The Fermi Galactic Center GeV Excess and Implications for Dark Matter

The analysis is based on 6.5 yr of Fermi LAT data recorded between 2008 August 4 and 2015 January 31 (Fermi Mission Elapsed Time 239,557,418–444,441,067 s). We select the standard good-time intervals, e.g., excluding calibration runs. In order not to be biased by residual backgrounds in all-sky or large-scale analysis, we select events belonging to the Pass 8 UltraCleanVeto class, which provides the highest purity against contamination from charged particles misclassified as gamma rays. Additionally, to minimize the contamination from emission from Earth's atmosphere, we select events with an angle θ < 90° with respect to the local zenith.

We use events with measured energies between 100 MeV and 1 TeV in 27 logarithmic energy bins, which is about seven bins per decade. For each energy bin events are binned spatially using HEALPix⁶⁹ (Górski et al. 2005) with a pixelization of order 6 (≈092 pixel size) or order 7 (≈046 pixel size). For all-sky fitting we use maps with adaptive resolution that have order 7 pixels in areas with large statistics (near the GP or close to bright sources) and order 6 in areas with fewer counts. The pixelization is determined based on the count map between 1.1 and 1.6 GeV at order 6, for which we further subdivide pixels with more than 100 photons to order 7. In the resulting maps we also mask 200 PSs with largest flux above 1 GeV from the Fermi LAT third source catalog (3FGL; Acero et al. 2015) within a radius of 1° (we mask pixels with centers within $1^\circ +a\sin (\pi /4)$ from the position of PSs, where a ≈ 046 is the size of pixels at order 7). Masking the bright PS effectively excludes pixels within about 2° from the GC (the fraction of masked pixels within 4° is about 50%; within 10° it is about 20%). We test the effect of the PS mask on the GC excess in Section 6.2, where we refit PS within 10° from the GC and only mask the 200 brightest PSs outside of 10°. The order 6 pixels at high latitudes are masked if any of the underlying order 7 pixels are masked.

We calculate the exposure and PSF using the standard Fermi LAT Science Tools package version 10-01-01 available from the Fermi Science Support Center⁷⁰ using the P8R2_ULTRACLEANVETO_V6 instrument response functions.

The emission measured by the LAT in any direction on the sky can be separated into individually detected sources, most of which are point-like sources, and diffuse emission. The majority of diffuse gamma-ray emission at GeV energies arises from inelastic hadronic collisions, mostly through the decay of neutral pions (π⁰). This component is produced in interactions of CR nuclei with interstellar gas; therefore, it is spatially correlated with the distribution of gas in the Milky Way. Another interstellar emission component, which becomes dominant at the highest and lowest gamma-ray energies, is due to inverse Compton (IC) scattering of leptonic CRs (electrons and positrons) interacting with the low-energy interstellar radiation field (ISRF). The ISRF can be considered to consist of three components: starlight, infrared light emitted by dust, and the cosmic microwave background (CMB) radiation. The IC contribution is expected to be less structured compared to the hadronic component. At energies ≲10 GeV, bremsstrahlung emission from electrons and positrons interacting with interstellar gas can become important. All of these three components of Galactic interstellar emission are brighter in the direction of the Galactic disk. Additionally, there is a diffuse emission component with approximately isotropic intensity over the sky. It is made of residual contamination from interactions of charged particles in the LAT misclassified as gamma rays, unresolved (i.e., not detected individually) extragalactic sources, and, possibly, truly diffuse extragalactic gamma-ray emission.

Throughout the paper, we will model Galactic interstellar emission from the large-scale CR populations in the Galaxy starting from predictions obtained through the GALPROP code⁷¹ (Moskalenko & Strong 1998; Strong et al. 2000, 2004; Ptuskin et al. 2006; Strong et al. 2007; Porter et al. 2008; Vladimirov et al. 2011). We use the GALPROP package v54.1 unless mentioned otherwise. GALPROP calculates the propagation and interactions of CRs in the Galaxy by numerically solving the transport equations given a model for the CR source distribution, a CR injection spectrum, and a model of targets for CR interactions. Parameters of the model are constrained to reproduce various CR observables, including CR secondary abundances and spectra obtained from direct measurements in the solar system, and diffuse gamma-ray and synchrotron emission. GALPROP is used to generate spatial templates for the gamma-ray emission produced in CR interactions with interstellar gas and radiation fields, which are then fitted to the data as described below.

The GALPROP models we employ in this paper assume CR diffusion with a Kolmogorov spectrum of interstellar turbulence plus reacceleration, and no convection. The diffusion coefficient is assumed to be constant and isotropic in the Galaxy. Additionally, unless otherwise mentioned, calculations assume azimuthal symmetry of the CR density with respect to the GC. For the Sample Model described in this section we chose one of the models from Ackermann et al. (2012). It assumes a CR source distribution traced by the measured distribution of pulsars (Lorimer et al. 2006, from now on referred to as Lorimer); the CR confinement volume has a height of 10 kpc and a radius of 20 kpc. It should be stressed that the parameters selected for the Sample Model represent only one of the possible choices. The goal of our analysis is not to find the best model, but rather to estimate the uncertainty in the GC excess due to the choice of parameters and analysis procedure. A study of the dependence of the results on propagation parameters and on the distribution of CR sources is presented in Section 4. Also, note that the Sample Model, as for most of the models considered in the paper, is derived from solving the transport equation in two dimensions (galactocentric radius, height over the GP). This speeds up the derivation of the model, but it is worth noting that cylindrical coordinates have a coordinate singularity at the GC. In our case, the modeling of interstellar emission is meant mainly for estimating the foreground/background emission (which dominates over emission from the region near the GC), and this is not a source of concern. In Section 4.4 we will employ some three-dimensional GALPROP models to address the case of CR sources near the GC.

The distribution of target gas is based on multiwavelength surveys. For the Sample Model described in this section, the calculation of gamma-ray fluxes from CR interactions employs the LAB survey (Kalberla et al. 2005) of the 21 cm line of H i and the survey of a 2.6 mm line of CO (a tracer of H₂) by Dame et al. (2001) to evaluate the distribution of atomic and molecular gas, respectively, in galactocentric annuli. The partitioning of the interstellar gas into galactocentric annuli based on the Doppler shifts of the lines is particularly uncertain at longitudes within about 10° of the GC, for latitudes within a few degrees of the Galactic equator. This is because the velocity from circular motion is almost perpendicular to the line of sight; therefore, Doppler shifts of the H i and CO lines are small relative to random and streaming motions of the interstellar medium in this range. The Sample Model taken from Ackermann et al. (2012) assumes H i column densities derived from the 21 cm line intensities for a spin temperature of 150 K. The dust reddening map of Schlegel et al. (1998) is used to correct the H i maps to account for the presence of dark neutral gas not traced by the combination of H i and CO surveys (Grenier et al. 2005; Ackermann et al. 2012). We neglect the contribution from ionized gas. The impact of the choice of input data for modeling the interstellar gas distribution is addressed in Section 4.3.

Since the distribution of CR densities in the Galaxy is not well constrained a priori, we fit the templates for the emission from interstellar gas split into galactocentric annuli to the gamma-ray data, independently for each energy bin, using the procedure described later. In the fit we use five independent annuli: three inner annuli spanning galactocentric radii 0–1.5 kpc, 1.5–3.5 kpc, and 3.5–8 kpc; a local ring spanning 8–10 kpc, and an outer ring spanning 10–50 kpc. Also, H i and CO maps are fitted independently to the data so that assuming an a priori CO-to-H₂ ratio is not required. In each ring, we add the bremsstrahlung and the hadronic components together in one template.

Furthermore, we separately fit to the data templates for the three IC components from CMB, dust infrared emission, and starlight, as models for the last two have significant uncertainties. Note that since we fit the model maps to the data in several independent energy bins and the morphology of gamma-ray emission from gas is determined mainly by the distribution of interstellar gas, the CR transport modeling in GALPROP in our case affects the analysis mainly through the morphology predicted for IC emission.

The LAT data revealed the presence of diffuse emission components extending over large fractions of the sky, which are not represented in the GALPROP templates that we use to model gamma-ray emission resulting from CR interactions in the Galaxy. Loop I is a giant radio loop spanning 100° on the sky (Large et al. 1962), which was also detected in LAT data (Casandjian & Grenier 2009). The origin of Loop I is an open question. It may be a local object, produced either by a nearby supernova explosion or by the wind activity of the Scorpio–Centaurus OB association at a distance of 170 pc (Wolleben 2007; Sun et al. 2015). Alternatively, it may be interpreted as the result of a large-scale outflow from the GC (Kataoka et al. 2013). In the Sample Model we account for Loop I using a geometric model (e.g., Figure 2 of Ackermann et al. 2014) based on a polarization survey at 1.4 GHz (Wolleben 2007). The geometric Loop I model assumes synchrotron emission from two shells. Each shell is described by five parameters: the center coordinates ℓ, b; the distance to the center d; and the inner (r_in) and outer (r_out) radius of the shell. The parameters are set to ℓ₁ = 341°, b₁ = 3°, d₁ = 78 pc, r_in,1 = 62 pc, r_out,1 = 81 pc, ℓ₂ = 332°, b₂ = 37°, d₂ = 95 pc, r_in,2 = 58 pc, and r_out,2 = 82 pc.

An additional large-scale extended emission component is represented by the so-called Fermi bubbles, two large gamma-ray lobes seen above and below the GC (Su et al. 2010; Ackermann et al. 2014). While the Fermi bubbles are well studied at high latitudes, a careful characterization of their properties close to the GP is complicated by large systematic uncertainties introduced by the modeling of other bright components of Galactic interstellar emission in this region. In the Sample Model we include a flat intensity spatial model of the Fermi bubbles at $| b| \gt 10^\circ$ (Figure 5 of Ackermann et al. 2014). In Section 5 we will derive an alternative template of the Fermi bubbles that includes emission at low latitudes.

We model the emission from the Sun (Moskalenko et al. 2006; Orlando & Strong 2007, 2008) and Moon, which is trailed along the ecliptic in our long data set, using templates derived with the Fermi Science Tools⁷² following the description in Johannesson et al. (2013).

For the Sample Model, we use the 3FGL catalog (Acero et al. 2015). We add all PSs in a single template in each energy bin. To construct the template, we use the parameterized spectra from the catalog and convolve with the PSF in each energy bin. Extended sources, except for the Large Magellanic cloud (LMC) and Cygnus region, are assembled in a separate template. Since the LMC and Cygnus are the brightest extended sources, we have independent templates. We describe extended sources based on the templates used in 3FGL. Unresolved Galactic sources, which may amount to up to ∼10% of the Galactic diffuse component (Acero et al. 2015), are not explicitly accounted for, but their spatial distribution is assumed to be similar to other Galactic components, and so the corresponding emission will be taken up by the other components in the fit (such as π⁰, bremsstrahlung, and IC).

We tentatively include in the Sample Model an additional component with the spatial distribution expected from annihilation of DM that follows in the inner Galaxy a generalized Navarro–Frenk–White (gNFW) profile with index γ = 1.25 and scaling radius r_s = 20 kpc. The NFW profile is an approximation to the equilibrium configuration of DM produced in simulations of collisionless particles (Navarro et al. 1997). Modified NFW profiles with γ > 1 are expected from numerical simulations of DM halos including baryons (e.g., Guedes et al. 2011). However, in our Sample Model the modified NFW profile is mainly motivated by earlier analyses of the GeV excess (Goodenough & Hooper 2009; Abazajian et al. 2014; Calore et al. 2015), which found it to be a good representation for the residual emission, with values of γ varying around ∼1.25. We will consider different γ values in Section 8.3.

The components of the Sample Model are summarized in Table 1.

We simultaneously fit the different components of diffuse emission (including an isotropic component) and a combined map of PSs to the Fermi LAT maps independently in each energy bin by maximizing the likelihood function based on Poisson statistics

$$\begin{eqnarray}&&\mathrm{log}{ \mathcal L }=\displaystyle \sum _{i}\,({d}_{i}\mathrm{log}{\mu }_{i}-{\mu }_{i}-\mathrm{log}({d}_{i}!)),\end{eqnarray} \tag{ 1 }$$

where d_i represents the photon counts in the spatial pixel with index i and μ_i represents the model counts in the same bin (since the fit is performed in each energy bin independently, we omit the energy bin index in this and the following equations). The model is constructed as a linear combination of templates

$$\begin{eqnarray}&&{\mu }_{i}=\displaystyle \sum _{m}\,{f}_{m}\,{P}_{i}^{(m)},\end{eqnarray} \tag{ 2 }$$

where m labels the components of emission and ${P}_{i}^{(m)}$ is the spatial template of component m in the appropriate energy bin corrected for exposure and convolved with the LAT PSF. The coefficients f_m are adjusted to maximize the likelihood.

Occasionally the best solution has negative normalization coefficients associated with some of the templates. Since this is unphysical, in such a case the corresponding template is removed and the fitting procedure is repeated. From comparing residual maps to the templates, it seems most likely that this behavior is due to incompleteness or imperfections of the models. The normalization of some of the templates is either over- or underestimated to compensate for such defects, and other templates' normalization in turn may react to this.

Our fitting strategy differs from previous examples in the literature (e.g., Calore et al. 2015; Ajello et al. 2016). We summarize below the main distinctive features:

• 1.  
  We fit the various emission templates independently in fine energy bins; this enables us to mitigate the impact on the results of several assumptions related to modeling background/foreground emission. We reiterate that the spectra of the various components, e.g., Figure 1, show that our procedure results in stable and physically plausible spectra.
• 2.  
  We perform an all-sky fit of all the component templates simultaneously; this provides us with a fast and simple procedure, which can also be consistently applied to other regions in the sky as a control sample to assess systematic uncertainties related to the GC results (Section 9). We characterize the impact of the analysis region choice on the excess spectrum in Section 3.2.

Zoom In Zoom Out Reset image size

The spectra of the components of the Sample Model fitted to the all-sky data are shown in Figure 1. The GC excess spectrum peaks around 3 GeV and extends up to about 100 GeV. The corresponding data maps, total model maps, and fractional residuals summed over several energy bins are shown in Figure 2. Although the Sample Model approximately reproduces the data, many excesses are evident. There is a clear residual associated with Loop I at energies below a few GeV in spite of including the geometrical template in the Sample Model. There are also residuals associated with substructures inside the Fermi bubbles. Furthermore, many excesses are seen along the GP. In Figure 3 we show the GC excess modeled by the gNFW annihilation template added back to the residual summed over energy bins between 1.1 and 6.5 GeV.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The analysis with the Sample Model also serves to confirm through inspection of the likelihood Hessian matrix that the large number of degrees of freedom does not create degeneracy between the model components, i.e., there is enough information in the gamma-ray data to separate them. However, some of the components that would be assigned negative fluxes are set to zero. As discussed before, this is most likely due to imperfections or incompleteness of the model. Although the procedure results in stable and physically plausible spectra for most of the various components, in a few instances this is not the case (e.g., in the Sample Model, for the gas rings between 1.5 and 3.5 kpc, the reason for which is discussed later in Section 4.3). This has limited impact on the determination of the GC excess properties, as the overall fore/background model is physically sound (Figure 1).

We start by testing the systematic uncertainty related to selection of the data sample. As an alternative to the sample analysis we use the Clean event class (P8R2_CLEAN_V6 instrument response functions) with a selection on zenith <100°. By considering the Clean class instead of the UltraCleanVeto, we estimate the magnitude of the residual CR contamination, which is larger for Clean class events compared to the UltraCleanVeto events. By using a larger zenith angle cut, we estimate a possible effect of emission from the Earth limb at ∼112°. In general, residual Earth limb emission, from gamma rays in the tails of the PSF, becomes more important at lower energies, where the tails are broadest. Using the Clean class events with the zenith angle cut <100° has relatively small influence on the GC excess spectrum (Figure 4, top left): the GC excess spectra are consistent within the statistical uncertainties.

We also subselect gamma-ray events with the best angular resolution: PSF classes 2 and 3 (Section 2.1). We then convolve the Sample Model components with the respective instrument response functions (notably, PSF) independently and perform a joint fit of the two data sets. The comparison of the GC excess flux with the Sample Model is again shown in the top left panel of Figure 4. There is a moderate effect on the spectrum at low energies only, where the LAT PSF gets worse.

One of the limitations of the template-fitting approach we use is that to model gamma-ray emission from gas we assume that the CR densities depend only on galactocentric radius and distance from the GP, and we rely on GALPROP to accurately predict the morphology of IC emission at each energy. Therefore, variations of the CR spectrum or mismodeling in one part of the Galaxy can lead to oversubtraction or unmodeled excesses in other regions.

One way to moderate this type of effect is to restrict the fitting procedure to a smaller ROI around the GC, so that there is more freedom to reproduce the features in the data for this specific part of the sky. To gauge the effect on the spectrum of the GC excess, we repeat the analysis in Section 2.2, restricting the ROI to some square regions: $| b| ,| {\ell }| \lt 10^\circ ,20^\circ ,30^\circ$. In this subsection we use maps with order 7 resolution (for all-sky fits we use adaptive resolution as discussed in Section 2.1), which gives more than 1000 pixels even for the $| b| ,| {\ell }| \lt 10^\circ$ case. This is generally sufficient to resolve the gas-correlated templates. However, the IC templates are rather smooth and may be degenerate in a small ROI. For this reason we combine the three IC templates in the Sample Model into a single template for fits in small ROIs. We also do not have the bubble template in the $| b| ,| {\ell }| \lt 10^\circ$ case, because it is defined only at $| b| \gt 10^\circ .$

The results are shown in the top right panel of Figure 4. We note that the gNFW cusp profile remains nondegenerate with the other components of emission even in the small ROI, because the degeneracies would result in large error bars, while the error bars on the GC excess flux remain reasonably small below 10 GeV. The intensity of the GC excess is generally reduced for the fits in smaller ROIs. For a 10° ROI the GC excess continues to be significant at energies below 400 MeV, while for a 30° ROI the excess cuts off below 1 GeV. The change in the GC excess flux for different ROI sizes is likely due to mismodeling of Galactic diffuse components.

Ackermann et al. (2012) explored the effects of varying several parameters of the GALPROP models that we use to create templates for interstellar gamma-ray emission. They concluded that the parameters with the largest impact on the predictions for gamma rays are (1) distribution of CR sources in the Galaxy, (2) height of the CR confinement halo, and (3) spin temperature used in deriving the atomic gas column densities from the 21 cm H i line intensities.

Our Sample Model in Section 2.2 uses the Lorimer pulsar distribution as a tracer of CR sources (supposedly SNRs, whose distribution is more difficult to determine from observations), a CR confinement height of 10 kpc, a radius of 20 kpc, and an H i spin temperature of 150 K. In order to quantify the impact of these choices on the spectrum of the GC excess, we use a subset of models in Ackermann et al. (2012). We have used different CR source distributions: an alternative pulsar distribution (Yusifov & Küçük 2004, hereafter referred to as Yusifov), the distribution of SNRs⁷³ (Case & Bhattacharya 1998), and the distribution of OB stars (Bronfman et al. 2000). Radial distributions of these CR source models are shown in Figure 5. We changed the CR confinement height from 10 to 4 kpc and its radius from 20 to 30 kpc. In addition, we derived the H i column densities from the 21 cm line intensities assuming an optically thin medium, which we formally modeled by setting the spin temperature to 10⁵ K.

Zoom In Zoom Out Reset image size

The resulting spectra for the GC excess are presented in the middle left panel of Figure 4. The largest effect is observed from the OB star source distribution model, which leads to an overall increase in the GC excess flux, while a decrease of the CR confinement height to 4 kpc leads to reduction of the flux at energies below a few GeV.

Zoom In Zoom Out Reset image size

IC emission is subdominant at GeV energies with respect to the gas-correlated components, especially near the GP. Its spatial distribution is expected to be smooth, but it depends on indirect knowledge of the ISRF and the calculated distribution of CR electrons. As a result, the IC emission is very difficult to model, especially near the GC, which can lead to significant uncertainties in the GC excess. In the Sample Model we use three IC components corresponding to the three seed ISRF components (CMB, starlight, and infrared), which are fitted to the gamma-ray data in independent energy bins. This procedure reduces the impact of our imperfect knowledge of the ISRF and CR electron spatial/spectral distribution. As an alternative approach, here we use a combined (i.e., summed over the three ISRF spectral bands) IC emission template, and we split the total IC emission into five galactocentric rings with the same boundaries as the gas templates.⁷⁴ The results are shown in the middle right panel of Figure 4. We find that there is a significant reduction of the GC excess flux between the models with combined ISRF components compared to models with ISRF components separated in different templates. We confirm that the IC emission can have a strong effect on the GC excess, which was previously discussed in, e.g., Ajello et al. (2016).

Uncertainties in the 3D models of the gas distribution in the Galaxy are important contributors to the uncertainties in the models of diffuse gamma-ray emission. In addition to the different values of the H i spin temperature considered in Section 4.1, in this section we explore (1) uncertainties due to the method used to partition the gas along the line of sight and (2) uncertainties related to the input interstellar tracer data and their angular resolution.

In the construction of the maps of interstellar gas used in Sample Model, Doppler shifts of the H i and CO lines are used to partition the gas column densities in annuli along the line of sight under the assumption of circular velocity around the GC. This method is not applicable toward the GC (and anticenter). For the Sample Model, the gas contents of the ring maps in this range are interpolated and renormalized as described in Appendix B of Ackermann et al. (2012). Also, in the case of CO, the line widths toward the GC are stretched by large noncircular motions (e.g., Dame et al. 2001). Therefore, all gas traced by CO at high velocities is assigned to the innermost ring based on the assumption that it is in the central molecular zone (CMZ). Furthermore, the H i maps in Ackermann et al. (2012) are augmented to incorporate dark neutral gas, i.e., neutral gas that is not traced by the combination of the H i and CO lines (Grenier et al. 2005). At low latitudes, where more than one ring map can have substantial column densities of interstellar gas, the inferred column densities of dark gas are distributed proportionally to the relative H i column densities in the rings. As noted in Ackermann et al. (2012), the correction for the dark gas component was limited to directions with $E(B-V)\lt 5$ mag, and special procedures were applied to handle large negative corrections. These considerations make the inferred column densities in the inner Galaxy relatively uncertain in the Sample Model.

Although Ackermann et al. (2012) assessed that this did not change the results of their large-scale analysis, we investigate here the impact on the properties of low-level residual emission seen toward the GC. Hence, we developed an alternative procedure to partition the gas column densities in the region at $| {\ell }| \lt 10^\circ$ (Appendix A) that employs complementary information from starlight (SL) extinction due to interstellar dust. Dust grains are thought to be well mixed with gas in the cold and warm phases of the interstellar medium, hence to be tracing the total (atomic and molecular) column densities (e.g., Bohlin et al. 1978). Marshall et al. (2006) proposed a method that combined infrared surveys with stellar population synthesis models to derive the distribution of dust in 3D. For our test we use the maps for the region of the Galactic bulge derived using this method by Schultheis et al. (2014) based on the VISTA Variables in the Via Lactea survey⁷⁵ (Minniti et al. 2010).

To further investigate the uncertainties related to the data sets used to build the gas maps, we considered alternative data sets that became recently available and, among other advantages, provide superior angular resolution. The H i maps in Ackermann et al. (2012) are based on the LAB survey (Kalberla et al. 2005), which has an effective angular resolution of ∼06, which is larger than the LAT angular resolution at energies ≳2 GeV. We produced alternative high-resolution maps by using the GASS survey (Kalberla et al. 2010). In the region of the sky where GASS data are available, including around the GC, they provide an angular resolution of ∼025.

Furthermore, the dark-gas correction in Ackermann et al. (2012) was based on the dust reddening map by Schlegel et al. (1998). The map by Schlegel et al. (1998) traces dust reddening based on dust thermal emission measured using IRAS and corrected for temperature variations using data from COBE/DIRBE. The latter has an angular resolution of ∼1°. For the alternative high-resolution maps we applied the same correction as in Ackermann et al. (2012), but based instead on the dust extinction map from Planck Collaboration et al. (2014) that is built using IRAS and Planck data (Planck public data release R1.20). The limit in reddening of 2 or 5 mag, above which no correction is applied in Ackermann et al. (2012), in our case was replaced by a comparable limit of 3 × 10⁻⁴ in dust optical depth at 353 GHz.

The effect of the alternative gas maps on the GC excess spectrum is shown in the bottom left panel of Figure 4. While the high-resolution gas maps only modestly change the excess spectrum, the alternative procedure to divide the gas into galactocentric rings yields an increase in the excess flux below a few GeV. The latter also results in more plausible spectra for the gas rings between 1.5 and 3.5 kpc that are set to zero in several energy bins in the Sample Model (2.4). The derivation of the gas distribution from dust extinction, however, has its own set of systematic and modeling uncertainties, which requires further investigation beyond the scope of the current analysis. Therefore, here we only use it to estimate the possible effect on the GC excess.

The large-scale distributions of CR sources that we consider peak at a few kiloparsecs from the GC, in correspondence with the so-called molecular ring and the main segment of the Scutum-Centaurus spiral arm. Some of them, notably the distribution used for the Sample Model, go to zero at the GC. However, there is evidence that a source of CRs up to PeV energies exists at the GC (HESS Collaboration et al. 2016; Gaggero et al. 2017). Gaggero et al. (2015) and Carlson et al. (2016) argued that taking into account additional sources near the GC may substantially change the significance and spectrum of the GC excess. We test two additional steady sources of CRs: one associated with the bulge/bar in the central kiloparsec of the Milky Way, and one associated with the CMZ in the innermost few hundred parsecs.

The stellar population of the Galactic bulge is older than 5 Gyr (e.g., Robin et al. 2012, and references therein). CRs that were accelerated by SNRs associated with the star formation activity in the bulge have either escaped the Galaxy or lost their energy (in the case of electrons). However, a possible source of CRs at the present time in the bulge is a population of MSPs that could potentially accelerate electrons and positrons to hundreds of GeV (e.g., Petrović et al. 2015). To model a possible population of MSPs in the bulge, we assume that their distribution is traced by the old stellar population in the bulge that we take from Robin et al. (2012). We parameterize the bulge as an ellipsoid with an orientation of 71 with respect to the Sun–GC direction (model B in Robin et al. 2012).

As an alternative, we consider a second possible population of CR sources distributed like the dense interstellar gas in the CMZ. The CMZ contains very dense molecular clouds that can host intensive star formation (e.g., Longmore et al. 2013) and, as a result, a significant rate of supernova explosions. The star formation rate (SFR) is rather uncertain in the CMZ and can vary from a few percent of the total SFR in the Galaxy, if traced by the free–free emission (Longmore et al. 2013), up to 10%–13%, if traced by young stellar objects (Yusef-Zadeh et al. 2009; Immer et al. 2012) or Wolf-Rayet stars (Rosslowe & Crowther 2015).

As a tracer of the CR production in the CMZ we use the distribution of molecular gas, which we model by a simplified axisymmetric version of Equation (18) in Ferrière et al. (2007). The radial distribution is described as

$$\begin{eqnarray}&&f{(R)\propto \exp (-(R/{L}_{c})}^{4})\times \exp (-| z| /{H}_{c}),\end{eqnarray} \tag{ 3 }$$

where R is the radial distance from the GC and z is the height above the GP. The two scaling factors were chosen to be L_c = 137 pc and H_c = 18 pc.

The additional source distributions, illustrated in Figure 5, are implemented in the GALPROP code to calculate the resulting gamma-ray emission. Owing to large uncertainty in their contributions, we treat these components independently from the rest of the templates. In the case of the bulge source, we add the IC emission from the additional electrons and positrons as an extra component, together with the components of the Sample Model in the fit to the data. In the case of the CMZ source, we add the IC emission template and the gas-correlated components associated with H i and H₂ in the first four rings in the Sample Model (Section 2.2), omitting the outer ring, i.e., we use four out of five rings in the Sample Model (nine additional parameters in each bin relative to the Sample Model).

Throughout our paper we are using GALPROP to model the particle propagation in two dimensions (for cylindrical symmetry in the Galaxy) with a 1 kpc resolution in radius and 100 pc resolution perpendicular to the disk. To have a more accurate description of the CR distribution in the CMZ case, we perform some 3D GALPROP runs for the CMZ source distribution with a resolution of 100 pc in all coordinates (the difference of the GC excess spectra for 2D and 3D GALPROP runs is less than about 2σ–3σ statistical uncertainties around a few GeV). For the CMZ source, we test different sizes of the propagation halo z = 2, 4, and 8 kpc and R = 10 kpc. The rest of the components are derived with a 2D GALPROP calculation with the same halo height and R = 20 kpc. The results are shown in the bottom right panel of Figure 4. The CMZ source of CRs with z = 2 and 4 kpc propagation halo height has little effect on the GC excess flux. The CMZ source of CRs with z = 8 kpc has a significant effect on the GC excess spectrum at energies below ∼4 GeV, while the bulge source of CR takes up a significant part of the GC excess around ∼10 GeV.

We derive the Fermi bubble template using the spectral component analysis (SCA) procedure used to extract the Fermi bubble component at high latitudes in Ackermann et al. (2014). In this derivation we will use the ROI $| b| \lt 60^\circ$, $| {\ell }| \lt 45^\circ .$

5.1.1. Subtraction of Gas-correlated Emission and Point Sources from the Data

The first step in modeling the Fermi bubbles is to subtract the gas-correlated emission and PSs from the data. We fit the data with a combination of gas-correlated emission components, the PS template obtained by adding 3FGL PSs (the overall normalization is free in each energy bin), and a combination of smooth templates. The smooth components are introduced as a proxy for the other components of emission, such as IC, Fermi bubbles, Loop I, extended sources, and GC excess. They are required to avoid biasing the determination of the contribution from the gas-correlated templates in the fit. As a basis of smooth functions we use spherical harmonics (calculated using the HEALPix package; Górski et al. 2005). The general basis of smooth functions makes it possible to model non-gas-related emission without a predefined template.

As a basis of smooth templates, we select the 30 spherical harmonics that provide the largest improvement in likelihood out of the first 100, i.e., ${Y}_{{lm}}(\theta ,\varphi )$ with degree l ≤ 9 (angular resolution ≈20°). In Section 5.1.3 we will test the consistency of the derivation by selecting the 60 most significant harmonics out of the first 225 (degree l ≤ 14, angular resolution ≈14°) and the 90 most significant harmonics out of 400 (degree l ≤ 19, angular resolution ≈10°). An example of data fitting by a combination of gas-correlated emission components, PSs, and spherical harmonics is shown in Figure 6.

Zoom In Zoom Out Reset image size

To speed up the calculations in this subsection, we use a quadratic approximation to the log likelihood in Equation (1),

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i}\,\displaystyle \frac{{({d}_{i}-{\mu }_{i})}^{2}}{{\sigma }_{i}^{2}},\end{eqnarray} \tag{ 4 }$$

where d_i is the photon counts in pixel i and μ_i is the model counts. The statistical uncertainty is calculated by smoothing the data count maps in each energy bin, ${\sigma }_{i}^{2}={\tilde{d}}_{i}$, to avoid bias in using either data or the model as an estimator of standard deviation (see, e.g., Appendix A in Ackermann et al. 2014). The smoothing radius R is chosen in each energy bin independently, such that there are at least 100 photons on average inside a circle of radius R. The minimum smoothing radius is 1°, while the maximum smoothing radius is 20°. For smoothing, photon counts inside the PS mask are approximated by an average of the neighboring pixels outside the mask.

5.1.2. Decomposition into Spectral Components

We use the results of the previous subsection to subtract the gas-correlated emission and PSs from the data. An example of the residual map summed over energy bins between 1.1 and 6.5 GeV is shown in the left panel of Figure 7. These residuals primarily consist of the Fermi bubbles, IC emission, isotropic background, and Loop I.

Zoom In Zoom Out Reset image size

At latitudes $| b| \gt 10^\circ$ the bubbles have an approximately uniform spectrum (e.g., Hooper & Slatyer 2013; Ackermann et al. 2014). We will assume that the spectrum of the bubbles at low latitudes is the same as at high latitudes and use the difference between the spectrum of the bubbles and that of other components to determine a template for the bubbles at low latitudes. The assumption about the bubble spectrum at low latitudes is a limitation of the current method, but, as we will see below, we will only need to use the spectrum in a relatively small energy range between 1 and 10 GeV. In this energy range the Fermi bubbles have a spectrum markedly different from the other gamma-ray emission components (see, e.g., Figure 1), and the LAT PSF is relatively good compared to energies below 1 GeV.

We further decompose the residuals, obtained by subtracting the gas-correlated emission and PSs from the data, as a combination of components correlated with the spectrum of the Fermi bubbles at high latitudes and the spectrum of the sum of IC, isotropic, and Loop I components. Between 1 and 10 GeV the high-latitude bubble spectrum is well fit by a power law ∝E^−1.9, while the sum of IC, isotropic, and loop I components obtained in the Sample Model (Section 2.2) is ∝E^−2.4. In this section we do not include a model for the GC excess component; we will take it into account as a separate spectral component in the next section. The model is determined as a combination of two spectral components

$$\begin{eqnarray}&&{M}_{\alpha i}={({E}_{\alpha }/{E}_{0})}^{-1.9}{H}_{i}+{({E}_{\alpha }/{E}_{0})}^{-2.4}{S}_{i},\end{eqnarray} \tag{ 5 }$$

where α is the energy bin index for energies between 1 and 10 GeV, i is the pixel index, and H and S are defined as a hard and a soft template. The reference energy is taken to be E₀ = 1 GeV; in this case the values of the H and S maps correspond to contributions at 1 GeV. The maps H_i and S_i are found by minimizing the χ² similar to the χ² in Equation (4):

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{1\,\mathrm{GeV}\lt {E}_{\alpha }\lt 10\,\mathrm{GeV}}\displaystyle \sum _{i}^{\mathrm{pixels}}\displaystyle \frac{{({R}_{\alpha i}-{M}_{\alpha i})}^{2}}{{\sigma }_{\alpha i}^{2}},\end{eqnarray} \tag{ 6 }$$

where ${R}_{\alpha i}$ are the residual maps obtained by subtracting gas-correlated emission and PSs from the data (Figure 7, left). The statistical uncertainty ${\sigma }_{\alpha i}^{2}$ is estimated from smoothed count maps (Section 2.3). We represent the maps H_i and S_i as linear combinations of residual maps: ${H}_{i}={\sum }_{\alpha }{f}_{{\rm{H}}}^{\alpha }{R}_{\alpha i}$ and ${S}_{i}={\sum }_{\alpha }{f}_{{\rm{S}}}^{\alpha }{R}_{\alpha i}$. The coefficients ${f}_{{\rm{H}}}^{\alpha }$ and ${f}_{{\rm{S}}}^{\alpha }$ are found by minimizing the χ² in Equation (6). The statistical uncertainties of H and S are calculated by propagating the uncertainties of the maps ${R}_{\alpha i}$, e.g., ${\sigma }_{{{\rm{H}}}_{i}}^{2}={\sum }_{\alpha }{{f}_{{\rm{H}}}^{\alpha }}^{2}{\sigma }_{\alpha i}^{2}$. The hard (H_i) and the soft (S_i) component maps are presented in Figure 7. The hard component primarily contains the Fermi bubbles and is used below to derive an all-sky model of the bubbles; the soft component contains isotropic background, IC emission, and Loop I.

5.1.3. Derivation of the Fermi Bubble Template

To derive the Fermi bubble template, we take the map of the hard spectral component ∝E^−1.9 in significance units smoothed with a 1° Gaussian kernel (Figure 8, top left) and cut in significance at the level of 2σ relative to the statistical uncertainty of the H_i map discussed after Equation (6). As one can see in Figure 8, the emission from the bubbles has a high significance and the cut at the 2σ level keeps most of the area of the bubbles. To eliminate the fluctuations and residuals outside of the Fermi bubbles, we select only the pixels that are above the threshold and are continuously connected to each other. The resulting Fermi bubble templates are shown in Figure 8.

Zoom In Zoom Out Reset image size

One of the main assumptions in this derivation of the Fermi bubbles is that their spectrum between 1 and 10 GeV below $| b| =10^\circ$ is the same as the spectrum above $| b| =10^\circ$. To reduce the dependence on this assumption, we split the derived Fermi bubble template into two templates: high latitude ($| b| \gt 10^\circ$) and low latitude ($| b| \lt 10^\circ$). The corresponding templates are shown in the bottom panels of Figure 8. The spectra of components derived with the new Fermi bubble templates in the Sample Model are shown in Figure 9. The spectrum of the low-latitude bubbles is similar to the spectrum of the bubbles at high latitudes between ∼100 MeV and ∼100 GeV, which supports the hypothesis of the homogeneous spectrum of the bubbles as a function of latitude. However, above 100 GeV the low-latitude bubble spectrum continues to be hard, while the high-latitude spectrum of the bubbles softens.

Zoom In Zoom Out Reset image size

The effect of the introduction of the low-latitude bubble template on the GC excess spectrum is shown in the right panel of Figure 9. Note that the Fermi bubble template in the Sample Model is determined only for $| b| \gt 10^\circ$. The GC excess above 10 GeV is taken up by the bubble template, while between 1 and 10 GeV the GC excess is reduced by a factor of 2 or more.

To test the robustness of the bubble template derivation and the effect on the GC excess flux, we also show the results for choosing different bases of smooth functions ℓ_max = 9, 14, and 19 (Section 5.1.1); different indices for the hard component n_hard = −1.8 and −2.0; different indices for the soft component n_soft = −2.3 and −2.5 (Section 5.1.2); and different significance thresholds in the derivation of the bubble template σ_cut = 1.8 and 2.2. The largest effect comes from the change in the soft component index n_soft = −2.3. The reason is that with the harder spectrum of the soft component a part of the bubble template is now attributed to the soft component. As a result, the bubble template has a smaller area, and it has a less significant influence on the GC excess flux.

In Figure 10 we show the residuals plus the GC excess modeled by the gNFW template with index γ = 1.25. We also show residuals in the model with all-sky bubbles without including a template for the GC excess. The excess remains in the presence of the all-sky bubble template, but it is reduced compared to the residuals in Figure 3. We note that Ajello et al. (2016) modeled the Fermi bubbles as an isotropic emission component within a 15° × 15° region around the GC, which led to a limited effect on the GC excess. This differs from our analysis, in which the Fermi bubbles have nonuniform intensity and become increasingly brighter near the GP, as derived from the SCA analysis. In conclusion, we find that the Fermi bubbles can significantly reduce the GC excess or even explain it completely above 10 GeV.

Zoom In Zoom Out Reset image size

In this section we apply the SCA technique to derive a template for the GC excess itself. Decomposition into spectral components was used previously by several groups to determine the morphology of the excess, in particular, de Boer et al. (2016) found that the excess emission resembles the distribution of molecular clouds near the GC, while Huang et al. (2016) argued that the excess morphology is spherical. Motivated by the possibility that the excess comes from a population of MSPs (Brandt & Kocsis 2015), we add the third spectral component with an average spectrum of observed MSPs $\propto {E}^{-1.6}{e}^{-E/4\mathrm{GeV}}$ (e.g., Cholis et al. 2014; McCann 2015). Consequently, we fit the residuals obtained after subtracting the gas-correlated emission and PSs in Section 5.1.1 between 1 and 10 GeV with three spectral components, hard $\propto {E}^{-1.6}{e}^{-E/4\mathrm{GeV}}$ (MSP like), medium ∝E^−1.9 (bubble like), and soft ∝E^−2.4. The derivation of the templates for the components is analogous to Equations (5) and (6), except that now there are three components instead of two. The maps of the templates are shown in Figure 11.

Zoom In Zoom Out Reset image size

The templates for the Fermi bubbles are derived by applying a cut in significance at 1.5σ to the medium spectral component. Owing to the presence of the third spectral component, the bubble component becomes relatively less significant; thus, we choose the 1.5σ cut rather than the 2σ cut used in the previous subsection. The template for the GC excess is derived by applying a 2σ cut in the hard component (Figure 12). The statistical uncertainties of the spectral component maps are derived by propagating the statistical uncertainties in the data maps (see discussion after Equation (6)). As before, we also split the Fermi bubble template into high- and low-latitude bubbles.

Zoom In Zoom Out Reset image size

The corresponding spectra for the GC excess and high- and low-latitude bubbles are shown in Figure 13. The spectra of the bubbles at high and low latitudes are consistent with each other between ∼1 and ∼100 GeV. At energies <10 GeV, the GC excess spectrum derived with the gNFW profile and the two-component SCA model of the bubbles is similar to the GC excess spectrum derived in the three-component SCA model (Figure 13, right).

Zoom In Zoom Out Reset image size

As an alternative derivation of the GC excess template, we use the spectrum $\propto {E}^{0.5}{e}^{-E/1.1\mathrm{GeV}}$ derived in Ajello et al. (2016) from the LAT data in the case of diffuse models with variable index and CR sources traced by distribution of pulsars. In this case the spectral shape is derived using a phenomenological spectral function to fit the LAT data and is not based on any specific scenario for the origin of the excess. The resulting spectrum for the alternative excess template is very similar to the spectrum derived with the template for the MSP-like spectrum.

In this section we present two PS search methods that were applied to the same data sets and diffuse models used in this work. The data selection is the same as for the Sample Model: 6.5 yr of UltraCleanVeto events with zenith angle cut θ < 90°. The goal is to test how much the selection of a PS detection algorithm can affect the inferred properties of the GC excess. Although both algorithms are based on a local likelihood method, there are differences in how the PSs are selected and localized. In both cases this is an iterative procedure from bright sources to faint ones, but the details are different, and we describe them in this subsection.

The first source detection algorithm is the same used in the production of the Fermi LAT source catalogs based on the pointlike package (Kerr 2010) and described in Acero et al. (2015). The data are binned in energy, in 14 bands from 100 MeV to 316 GeV (or four per decade), and separated into front and back event types. For each band and event type, the photons are binned using HEALPix, with pixel sizes selected to be small compared with the PSF. The log likelihood is then computed summing over the energy bands and event types. As described for 3FGL (Section 3.1.2 of Acero et al. 2015), the contributions to the likelihood function are "unweighted" for the lower energies, to account for systematics of the diffuse background spectrum. For the diffuse model, we use the Sample Model from Section 2.2 fitted to the data, but without adding the gNFW template.

The data are fitted using a diffuse model template, an isotropic template, and PSs in small ROIs covering the whole sky. The centers of the ROIs are determined by the centers of HEALPix pixels with nside = 12 (1728 tiles in total; average distance between the centers of ROIs is about 5°), and the radius of the ROIs is 5°. The pixel size depends on the energy. At energies below 10 GeV it is about 1/5 of the 68% containment radius, e.g., for back-entering events around 100 MeV the pixel size is about 1° (nside = 52). Note that nside = 12 and nside = 52 are nonstandard nside values, which are typically powers of 2. Above 10 GeV the analysis is unbinned. The likelihood for the data within the radius is optimized with respect to the spectral parameters of the sources located within the tile. Correlations with sources outside the tile, but contributing to the likelihood, are accounted for by iterating until changes of the likelihood for each tile are small.

In the search for new PSs, we calculate the likelihood ratio, expressed as a Test Statistic ($\mathrm{TS}=2{\rm{\Delta }}\mathrm{log}{ \mathcal L })$ for an additional PS, assuming a power-law spectrum with a fixed spectral index of 2.3 but variable flux, at each of the positions defined by nside = 512 [3.2M total]. This is done for all pixels within each ROI. A clustering analysis is applied to the resulting map of the pixels with TS > 25. All clusters with more than one pixel are used to define seeds for inclusion in the model. As for all sources, the spectral index is now optimized and the source is localized. If the power-law spectrum does not fit the data well, then the spectrum is described by a log parabola (e.g., Acero et al. 2015). A source candidate is accepted for inclusion if its optimized TS is greater than 25 and the localization process converged properly.

This source-finding procedure relies on the model being an accurate description of the data, given the set of sources and diffuse components. Thus, the set of sources needs to be fairly complete, so that new sources are weak and do not strongly affect the current model. We have found that it is necessary to rerun the procedure several times after adding new sources.

For the determination of the second new list of PS, we use the Fermipy package, a set of Python tools built around the Fermi LAT Science Tools that automate and enhance their functionalities.⁷⁶ In this case we use data between 300 MeV and 550 GeV binned with 5 bins per decade. The ROI is $| {\ell }| ,| b| \lt 22^\circ$, and we bin the data in 008 pixels (on a square grid). As a preliminary step, we start with the 3FGL PS and reoptimize their positions. We then perform a fit to the ROI with those sources and delete all 3FGL sources with TS < 49. Then we build a TS map (map of TS of a PS candidate at each position in the spatial grid) and select maxima with TS > 64 and separation from other sources greater than 05. The best positions and location uncertainties of source candidates are derived from the likelihood profile in the nine pixels around the TS map peak by fitting it with a paraboloid. We repeat the selection of TS maxima and PS localizations two more times for TS > 36 and separation greater than 04, and with TS > 25 and separation larger than 03. After this third iteration, we have a list of sources with TS > 25, and we perform a final fit to the ROI to determine the PS spectra.

Discussing the properties of the source candidates in the new lists is beyond the scope of this article. Instead, we focus in the following section on the effects that the choice of a PS list has on the determination of the GC excess spectrum.

We start the analysis of the effect of PS characterization near the GC by combining PSs within 10° from the GC into independent templates, using the spectra provided by the 3FGL catalog, pointlike or Fermipy. We use one of these templates at a time together with the GC excess template, the other diffuse emission components, and the PS template determined with sources outside of 10° from the GC. The effect on the GC excess spectrum is shown in the left panel of Figure 14. An independent template for sources inside 10° generally leads to a softer GC excess flux at energies below 1 GeV, while above 1 GeV the GC excess flux is not affected significantly by the introduction of the independent PS template for sources within 10°. To test the effect of the PS mask, we also include a model with 3FGL templates where we do not mask PSs within 10° from the GC. In this case we find an overall reduction of the flux attributed to the GC excess. The 1FIG list is not considered in this step because the spectra are provided only for energies larger than 1 GeV. The spectra for the Fermipy sources are extrapolated from the energies above 300 MeV used to determine the source list, which is likely the cause for the deviations of the GC excess spectrum at lower energies with respect to the other determinations.

Zoom In Zoom Out Reset image size

To allow more freedom in PS modeling, we also refit the spectra of the PSs within 10° from the GC with a free normalization for each source in each energy bin independently. The 3FGL Catalog has 76 sources within this region, and so our model has 76 additional free parameters in each energy bin compared to the Sample Model. The 1FIG catalog has 48 sources inside $| b| ,| {\ell }| \lt 7\buildrel{\circ}\over{.} 5.$ Pointlike provides 104 PS candidates, and Fermipy has 127 PS candidates inside R < 10°. Due to limited statistics at high energies, we restrict the energy range to 100 MeV–300 GeV (23 energy bins) for these fits. Similarly to the Sample Model, PSs more than 10° away from the GC (outside of $| b| ,| {\ell }| \lt 7\buildrel{\circ}\over{.} 5$ in the case of 1FIG) are combined in a single template based on the 3FGL catalog. Diffuse emission components, extended sources, the LMC, and Cygnus also have independent templates: all these templates have free overall normalization in each energy bin. We mask regions of 1° radius for the 200 brightest PSs outside of 10° from the GC.

The spectra derived for the GC excess are shown in the right panel of Figure 14. The spectrum below a few GeV is clearly dependent on the data set, the diffuse emission model, and the method used to derive the PS source list. In several cases the flux attributed to the GC excess is larger than in the Sample Model. In the case of the 3FGL Catalog, this could be attributed to an overestimation of the PS fluxes near the GC stemming from not having accounted for an excess at the GC explicitly. The largest effect is observed for refitting 1FIG sources that were derived using data above 1 GeV, without accounting in any way for a GC excess, and 1FIG has the fewest number of PSs near the GC, due to rather strict selection criteria (Ajello et al. 2016). On the other hand, refitting sources individually above 10 GeV reduces the flux attributed to the excess for all PS lists. This is consistent with the analysis of Linden et al. (2016), who find that the GC excess above ∼10 GeV can be explained by a contribution of PSs.

In the derivation of the new lists of PS, we did not account for large-scale residuals, as done, e.g., in Acero et al. (2016) for the derivation of the 3FGL catalog, or specifically for an excess near the GC. This is likely to result in several spurious sources that absorb positive residuals due to underprediction of diffuse emission. We use them as a conservative starting point in the derivation of the GC excess component, since some part of the emission from this component will be attributed to spurious PSs. By allowing free normalizations for both the gNFW template and the PSs, we have also tested whether the gNFW template is preferred relative to a combination of resolved PSs from a statistical point of view. We find that in some cases refitting PSs results in larger flux attributed to the GC excess at energies below a few GeV; however, for the pointlike algorithm, refitting the PSs leads to smaller GC excess flux for all energies. The conclusion is that although all PS detection algorithms are based on maximizing a local likelihood, the details on how the PSs are selected are important and may have a significant influence on the GC excess.

In this section, we derive the spectrum independently in sectors along and perpendicular to the GP by dividing the gNFW template into four templates centered at the GC with an opening angle of 90°: top, bottom, left, and right. The fit is all-sky, and the only difference from the Sample Model is that now we have four independent templates for the excess in the four quadrants instead of a single one. The spectra of the quadrant templates are shown in the left panel of Figure 16. The spectra of the top, bottom, and right quadrants are similar to each other, while the spectrum in the left quadrant is different from the other three.

Zoom In Zoom Out Reset image size

Let us recall that in the Sample Model, the Fermi bubble template is defined only for $| b| \gt 10^\circ .$ The difference in the spectra of the quadrant templates may be due to an asymmetry in emission from the Fermi bubbles near the GC (see Figure 8). To qualitatively investigate this hypothesis, we "correct" the spectrum of the left quadrant by adding a bubble-like spectrum (as derived in the Sample Model in Figure 1) and a background-like spectrum (∝E^−2.4, the same as the soft spectral component in the derivation of the all-sky bubble template) with free normalizations adjusted so that differences with respect to the bottom quadrant are minimized. We take the bottom quadrant as a reference since it is the least contaminated by emission from the GP and local gas to the north of the GC. The "corrected" left quadrant spectrum is compared to the others in the right panel of Figure 16. This shows that it is plausible to have a spherical excess on top of emission from the Fermi bubbles and other foregrounds. Another implication is that the bubbles are contributing mostly to the right, top, and bottom quadrants, i.e., the contribution is asymmetric with respect to the GC, which is consistent with the description of the bubbles by Acero et al. (2016).

Negative values of the GC excess spectrum below 1 GeV in, e.g., the right panel of Figure 16 are troubling. Throughout the paper, we restrict the values of the diffuse emission components to be non-negative, but since the presence of the GC excess is not known a priori, we treat it as a "residual" component and allow both positive and negative values. If the morphology of the foreground components is not modeled perfectly, the fit can try to compensate for imperfections by subtracting the GC excess template such that the correction can be larger than the otherwise positive flux from the excess component. Indeed, we notice that the contribution from the background-like component in the right panel of Figure 16 is negative, i.e., one needs to subtract the background-like spectrum from the left quadrant to minimize differences with the bottom quadrant, and the GC excess spectrum in the four quadrants is mostly positive with respect to this "corrected" zero level.

Figure 17 shows the spectra of the four quadrants once the all-sky bubble template derived in Section 5 is included in the fit. Similarly to Figure 13, the spectrum of the quadrant templates changes dramatically, and emission remains significant only below 10 GeV. The spectra in the four quadrants are closer to each other, but not yet consistent within the statistical uncertainties. This inconsistency, however, may be due to an imperfect derivation of the Fermi bubble template.

Zoom In Zoom Out Reset image size

In summary, we find that establishing whether the GC excess has spherical morphology is challenging, due to uncertainties in the contribution from low-latitude emission from the Fermi bubbles. However, at present we cannot exclude that a component with spherical morphology is present in addition to a continuation of the Fermi bubbles.

In Section 5.2, we derived templates for the emission near the GC correlated with the Fermi bubble spectrum at high latitudes and with an average MSP spectrum. Longitude and latitude profiles for the component with a bubble-like spectrum (Figure 11 middle) are presented in Figure 18. The profiles are shown at a reference energy of 2 GeV. The latitude profiles are relatively flat for $10^\circ \lesssim | b| \lesssim 50^\circ ,$ but the intensity increases by a factor of ∼5 near the GP. One can also see that the emission associated with the Fermi bubbles in this model is shifted to the right (negative longitudes) relative to the GC.

Zoom In Zoom Out Reset image size

Similarly, longitude and latitude profiles of the MSP-like component (Figure 11, right panel) are shown in Figure 19. The latitude and longitude profiles of this component are symmetric with respect to the GC, with a possible enhancement along the GP, which can be expected as a contribution from millisecond and regular pulsars in the Galactic disk (e.g., Faucher-Giguère & Loeb 2010; Grégoire & Knödlseder 2013).

Zoom In Zoom Out Reset image size

Finally, in Figure 20 we compare the profile as a function of radial distance from the GC at 2 GeV for the MSP-like spectral component with the total gamma-ray data and the gNFW profiles in the Sample Model, as well as for a standard NFW annihilation profile. The MSP-like profile is similar to the DM annihilation profiles (gNFW with γ = 1.25 in the Sample Model and the NFW profile) within ∼5° of the GC, but it flattens at a higher intensity than the gNFW profile, which is likely related to the positive values of the MSP-like component along the disk; cf., the longitude profile in the right panel of Figure 19.

Zoom In Zoom Out Reset image size

We also checked that using alternative PS templates within 10° from the GC derived for UltraCleanVeto data with pointlike and Fermipy tools (Section 6.2) does not significantly affect any of the profiles for the MSP-like component.

In summary, the profiles in latitude, longitude, and radial distance from the GC corroborate the hypothesis that the excess is not obviously consistent with expectations from DM annihilation with gNFW/NFW density profiles, but such a component may exist in addition to emission from the Fermi bubbles and from sources in the Galactic disk/bulge such as MSPs.

In this section, we assess the relative likelihoods of models in which we vary the centroid position of the gNFW annihilation template around the GC, as well as its radial index γ. Results from the scan in the position of the center of the gNFW template are shown in Figure 21. The spectra of the excess for cases with the component centered at b = 0° and with various longitudes are presented in the left panel of Figure 21, while the $-2{\rm{\Delta }}\mathrm{log}{ \mathcal L }$ values for different locations around the GC are shown in the right panel of Figure 21. The best-fit position is at l ≈ −1°. The spectrum of the excess depends on the location of the centroid. The spectra for the center at positive longitudes look similar to the GC excess in the left quadrant in the left panel of Figure 16, while the spectra for negative longitudes resemble more the spectrum of the Fermi bubbles with a less pronounced bump at a few GeV and the spectrum extending to lower energies. These findings are consistent with the possibility that the GC excess to the right (negative longitudes) from the GC is mixed with a contribution from low-latitude emission of the Fermi bubbles above 10 GeV. This is also consistent with the observation by Calore et al. (2015) that the best-fit longitude of the gNFW profile is at l ≈ −1° below 10 GeV and shifts to l ≲ −2° above 10 GeV.

Zoom In Zoom Out Reset image size

The $-2{\rm{\Delta }}\mathrm{log}{ \mathcal L }$ values for the variations of the gNFW center when the model includes the all-sky bubble template (i.e., including the component at low latitudes) are shown in Figure 22. In the left panel we show results for all-sky gNFW templates, while in the right panel we truncate the gNFW template at 10° from its center to test whether the difference in the best-fit location of the center is due to residuals away from the GC. Both the truncated and the all-sky gNFW templates have the best-fit position at b = 3°, ℓ = 1°. The range of $-2{\rm{\Delta }}\mathrm{log}{ \mathcal L }$ values is smaller when the all-sky bubble component is included in the model.

Zoom In Zoom Out Reset image size

The spectra of the GC excess and the $-2{\rm{\Delta }}\mathrm{log}{ \mathcal L }$ for different indices of the gNFW profile when all other components are accounted for as in the Sample Model are shown in Figure 23. The step in the index scan was 0.1. The best likelihood is obtained for the standard NFW profile with radial index γ = 1. Scanning the profile with other models from Section 4 and from Section 6 (including independent PS templates within 10° from the GC), we find that for most of the models the best-fit indices are between 0.9 and 1.2 (the scan step is 0.1 in each case), which overlaps with the range of indices found by Calore et al. (2015) and with the best-fit index of 1.2 found by Hooper & Slatyer (2013).

Zoom In Zoom Out Reset image size

The results of a scan of the all-sky gNFW profile in the presence of the all-sky bubble template derived in Section 5.1 are shown in the left panel of Figure 24. The best-fit index is equal to 0.4, which is significantly smaller than the best-fit index γ = 1 obtained in the scan with the Sample Model. The results of the scan for the truncated gNFW template are shown in the right panel of Figure 24. The best-fit index for the truncated template is 0.9.

Zoom In Zoom Out Reset image size

The rather large variations of the best-fit index and the gNFW centroid in the presence of the all-sky bubble template show that the inferred morphology of the excess critically depends on the model of the Fermi bubbles near the GC, which can bias the derivation of the morphology and, as we have shown in Section 5, the GC excess spectrum. Because of these large uncertainties, at present it is not possible to firmly associate the centroid of the excess with the GC itself or precisely determine its density profile.

Many groups have shown that the spectral energy distribution of the GC excess peaks at energies around a few GeV and can be fit with a model of either ∼40 GeV DM annihilating to $b\bar{b}$ or ∼10 GeV DM annihilating to τ⁺τ⁻ (Goodenough & Hooper 2009; Hooper & Goodenough 2011; Hooper & Linden 2011; Gordon & Macías 2013; Hooper & Slatyer 2013; Abazajian et al. 2014; Calore et al. 2015; Daylan et al. 2016).

We used DM annihilation spectra for a variety of DM masses and two representative DM annihilation channels, $b\bar{b}$ and τ⁺τ⁻, to model the GC excess spectrum that we found using our Sample Model (see Appendix C.1 for details). To estimate the uncertainty level of the DM-like signal, we repeat the analysis by placing the gNFW template at different locations along the GP instead of the GC. Since we compare fits from many regions across the GP with varying levels of γ-ray intensity, we quantify the best-fit DM component as a fraction of the effective background:

$$\begin{eqnarray}&&f=\displaystyle \frac{{N}_{\mathrm{sig}}}{{b}_{\mathrm{eff}}},\end{eqnarray} \tag{ 7 }$$

where N_sig is the number of signal counts integrated over the energy bins and ${b}_{\mathrm{eff}}$ is the "number of counts" in the effective background. If the signal were localized in a small region on the sky with expected intensity much smaller than the background intensity, then the statistical variance of the signal measurement would be proportional to the number of background counts in that region. In general, for a signal that covers a large portion of the sky (or possibly the entire sky) with a varying intensity, one can determine a weighted sum of the background counts, the effective background, so that the statistical variance of the signal is still proportional to this weighted sum (see Ackermann et al. 2015; Buckley et al. 2015; Caputo et al. 2016; Appendix C, for details about the evaluation of the effective background):

$$\begin{eqnarray}&&{b}_{\mathrm{eff}}=\displaystyle \frac{N}{{\sum }_{i,\alpha }\tfrac{{({P}_{i\alpha }^{(\mathrm{sig})})}^{2}}{{P}_{i\alpha }^{(\mathrm{bkg})}}-1},\end{eqnarray} \tag{ 8 }$$

where the sum is over energy indices α and pixel indices i, N is the total number of counts, and ${P}_{i\alpha }^{(\mathrm{sig})}$ (${P}_{i\alpha }^{(\mathrm{bkg})}$) are the signal (background) intensity distributions normalized to 1: ${\sum }_{i,\alpha }{P}_{i\alpha }^{(\mathrm{sig})}={\sum }_{i,\alpha }{P}_{i\alpha }^{(\mathrm{bkg})}=1$. In a particular case, if there is one energy bin and the background is uniformly distributed over the whole sky ${P}_{i}^{(\mathrm{bkg})}=1/{N}_{\mathrm{pix}}$, while the signal is uniformly distributed over a small number of pixels k_pix so that ${P}_{i}^{(\mathrm{sig})}=1/{k}_{\mathrm{pix}}$ and N_pix/k_pix ≫ 1, then Equation (8) gives ${b}_{\mathrm{eff}}\approx N\cdot {k}_{\mathrm{pix}}/{N}_{\mathrm{pix}}$, which is the number of background counts in the signal region.

Figure 25 (left) shows the GC excess spectrum in the Sample Model (Section 2.2) in count space. In this fit, the normalization of the GC excess template is fit to the data, together with the other templates in each energy bin. The first four energy bins had negative best-fit normalizations of (−3.4 ± 0.3) × 10⁴, (−2.3 ± 0.3) × 10⁴, (−1.5 ± 0.3) × 10⁴, and (−1.3 ± 1.8) × 10³ respectively. The errors given are the statistical errors only. Figure 25 also shows the best-fit DM annihilation spectrum for the $b\bar{b}$ channel for various values of the mass of the DM particle, m_DM. The DM annihilation count spectra were calculated from their corresponding flux spectra using the DMFitFunction within the standard Fermi Science Tools and the detector exposure for this data set. The signal counts N_sig for each DM model are simply the integral of the best-fit count spectrum. Therefore, we can evaluate the strength of the best-fitting DM model relative to the effective background ${b}_{\mathrm{eff}}$ (see Figure 25, right).

Zoom In Zoom Out Reset image size

We note that none of the DM fits to the GC excess spectrum are very good (the reduced χ² is >30 for all fits). In particular, the sample spectrum has a high-energy (>50 GeV) tail that cannot be explained by, e.g., <50 GeV DM annihilating to $b\bar{b}$. However, as we have shown, the high-energy tail may be consistent with a low-latitude extension of the Fermi bubbles (Section 5.1.3). Nevertheless, we consider the results of these fits as we are interested in quantifying the "DM-like" component of the spectrum for various DM models.

In addition to the gNFW excess template, we also fit the standard NFW (γ = 1.0) DM annihilation template since it is the best-fit template for the excess in our analysis (Section 8.3). Figure 26 shows the best-fit N_sig of various DM models to the GC excess spectra both in the Sample Model and in the model including the SCA bubbles (Section 5.1.3) as a fraction of ${b}_{\mathrm{eff}}$. The bubble template derived using the SCA method can account for a large amount of the GC excess, especially at high energies. Therefore, the reduction of the amplitudes of the best-fit DM models (especially at high masses) is expected.

Zoom In Zoom Out Reset image size

Because the GC is very bright in γ rays, many of the DM models we test have very small statistical errors in inferred N_sig (δN_sig < 0.01). However, we are not able to model the γ-ray sky to a similar level of precision (recall that the fractional residuals from our fits are typically in the range of −0.2 to 0.2; see Figure 3). Therefore, systematic uncertainties that may mimic or mask a DM signal need to be accounted for.

To assess these systematic uncertainties beyond what was already done with model variations, we estimate $\delta {N}_{\mathrm{sig},\mathrm{syst}}$ by fitting for DM-like signals in control regions along the GP, based on two assumptions: that the expected DM signal is approximately zero for 30° ≤ l ≤ 330°, and the systematic uncertainty scales with ${b}_{\mathrm{eff}}$ for effects that can induce or mask a DM-like signal. An excess may be a fraction of the background if it is caused by a single (or a few) errors in the modeling of the gamma-ray intensity, which are proportional to the "average" emission, or when the uncertainty is dominated by errors in a single component that also dominates the overall emission. Fluctuations due to several small effects, such as uncertainties in emission components where each component contributes a small fraction of the total emission, would be best estimated as a square root of the ${b}_{\mathrm{eff}};$ in this case the characteristic values would be ${N}_{\mathrm{sig}}/\sqrt{{b}_{\mathrm{eff}}}$. Fluctuations in emission that are caused by one or a few components that are not directly correlated with the overall gamma-ray emission, such as a local SNR or an AGN-like activity, would be best characterized by their absolute values. Since the gamma-ray emission toward the GC is the largest, taking the fractional excess as a figure of merit to estimate its significance is the most conservative assumption, which we will adopt for our analysis.

Control regions along the GP to estimate the modeling uncertainty were used before by Calore et al. (2015). They fit a DM-like spatial profile along the GP and represented the results as a covariance matrix in energy bins, which is used to determine the expected level of modeling uncertainty at the GC. Our approach is to fit DM-like excess along the GP including both the spatial profile and the energy spectrum of a DM annihilation channel. We then express the uncertainty as a ratio of the signal to the local effective background. Both of these differences are likely to increase the estimate of the modeling uncertainty, since we get the maximal possible DM-like signal in each location, and then we divide by the local background, which is smaller along the plane than at the GC.

We perform all-sky fits using the same diffuse emission components as in the Sample Model, but we shift the gNFW excess template in steps of 10° in longitude at b = 0 for 30° ≤ l ≤ 330°. Figure 27 shows the amplitudes of the best-fit DM model spectra (as a fraction of ${b}_{\mathrm{eff}}$) measured in the control regions.

Zoom In Zoom Out Reset image size

Though the fits were allowed to be negative, most longitudes preferred a positive DM normalization. Also, about half of the total fits had $| f| \lt 0.01$. We define δf_GP as the value of f for which we obtain larger values in only 16% of the fits along the GP, which corresponds to a one-sided 1σ exclusion in the case of a Gaussian distribution. We take δf_GP as representative of the amplitude of the components of the diffuse γ-ray residuals consistent with DM signals, and therefore a measure of the systematic uncertainty for our DM search in the GC:

$$\begin{eqnarray}&&\delta {f}_{\mathrm{syst}}=\delta {f}_{\mathrm{GP}},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\delta {N}_{\mathrm{sig},\mathrm{syst}}={b}_{\mathrm{eff}}\times \delta {f}_{\mathrm{syst}}.\end{eqnarray} \tag{ 10 }$$

See Appendix C for more on the motivation of this choice and the application of δf_GP in our fitting.

A comparison of the fits in the GC to the characteristic $\delta {N}_{\mathrm{sig},\mathrm{syst}}$ in Equation (10) is shown in Figure 28. The largest DM signal as a fraction of ${b}_{\mathrm{eff}}$ in the GC in the Sample Model is similar to the characteristic uncertainty level from the GP scans. Consequently, we cannot claim that the DM interpretation of the signal in the GC is robust when we compare it with the DM-like signals in the control regions along the GP. We note that the same conclusion holds when an NFW profile is adopted for the spatial distribution of the DM, and for the model that includes all-sky SCA bubbles. In the latter case, the GC signal is much smaller than those seen in the GP scan. Furthermore, the same conclusions are reached if we adopt higher energy thresholds of 300 MeV and 1 GeV for the analysis.

Zoom In Zoom Out Reset image size

As a corollary from these results, we conclude that the model variations discussed earlier in the paper may not capture the complete range of systematic uncertainties that can mimic a DM signal. Either the variations considered do not cover the full range of the associated model uncertainties, or there are other sources of gamma rays that play a significant role. A noteworthy candidate for the latter is a population of unresolved sources, such as MSPs, that other analyses (Bartels et al. 2016; Lee et al. 2016; Brandt & Kocsis 2015) are indicating as a likely cause of the GC excess.

In the previous subsection we have found that, although the GC excess is statistically significant and it is present in all models that we have considered, a similar level of fractional excesses is found at other locations along the GP, where no DM annihilation is expected. As a consequence, a DM interpretation of the GC excess is not robust. Since we cannot claim a detection of DM annihilation at the GC, we derive upper limits on the DM annihilation cross-section. In previous studies, e.g., Ackermann et al. (2015) and Buckley et al. (2015), the systematic uncertainty inferred from control regions was used also to set upper limits on the DM cross-section. In Appendix C.1 we use the modeling uncertainty derived from the scan along the GP in a derivation of upper limits on DM annihilation. This derivation assumes that the probability of having negative fluctuations (i.e., masking a real DM signal) is as large as the positive fluctuations seen in the GP scan. However, this is likely to be too conservative since in almost all fits in control regions along the GP we found positive DM-like excesses (Figure 27).

Therefore, here we proceed to set limits on $\langle \sigma v\rangle$ by requiring a tentative DM signal to not exceed the largest GC excess flux value in each energy bin for all the background models considered in this work (i.e., the upper edge of the blue band in Figure 15) at the 95% statistical confidence level. Results are shown in Figure 29 and compared with other relevant DM limits from the literature.

Zoom In Zoom Out Reset image size

Our limits are similar to the DM parameters found in works where the GC excess has been interpreted as DM. This is not surprising since all the independent GC excess analyses use similar data sets and interstellar emission models. We also see that the gNFW limits are stronger than the NFW limits since the gNFW profile has a steeper inner slope of the DM density toward the GC. Both gNFW profiles provide stronger limits than the dSph analysis at masses above a few TeV (few hundred GeV) for the $b\bar{b}$ (τ⁺τ⁻) channel. Although for the Sample Model the GC analysis also provides stronger than dSph DM limits below few tens of GeV for the $b\bar{b}$ channel, these limits are subject to larger modeling uncertainties. This is the regime where the analyses become statistics limited. Since the expected signal in our GC analysis is much larger than that of the dSph analysis, the statistical uncertainties are smaller, resulting in more constraining limits.

In Section 9, we fit a variety of DM models to the GC excess spectrum that we found using our sample background model. The differential γ-ray flux expected from DM annihilation in a given ROI (ΔΩ) is

$$\begin{eqnarray}&&\displaystyle \frac{d\phi }{{{dE}}_{\gamma }}=\left(\displaystyle \frac{\langle \sigma v\rangle }{8\pi }\displaystyle \frac{{{dN}}_{\gamma }}{{{dE}}_{\gamma }}\displaystyle \frac{1}{{m}_{\chi }^{2}}\right)\left({\int }_{{\rm{\Delta }}{\rm{\Omega }}}d{\rm{\Omega }}{\int }_{{\rm{l}}.{\rm{o}}.{\rm{s}}.}d{\ell }\,{\rho }_{\chi }^{2}({\ell })\right),\end{eqnarray} \tag{ 18 }$$

where dN_γ/dE_γ is the differential count spectrum of gamma rays from annihilation of a pair of DM particles, m_χ is the mass of each DM particle, $\langle \sigma v\rangle$ is the velocity-averaged annihilation cross-section, and ${\rho }_{\chi }^{2}$ is the density distribution of the DM. dN_γ/dE_γ depends on the relative strength of each annihilation channel for a specific DM model. In this work, we calculate limits on DM annihilation for two representative annihilation channels: $b\bar{b}$ and τ⁺τ⁻. The value in the first set of parentheses in Equation (18) depends on the particle nature of the DM, specifically on the particle mass, annihilation channel, and annihilation cross-section. The value in the second set of parentheses is the so-called "J factor," which depends on the distribution of DM. We consider both a standard NFW profile (γ = 1) and a slightly contracted generalized NFW profile (γ = 1.25). The J factors integrated over the whole sky (including the PS masking) assuming a local DM density of 0.4 GeV cm⁻³ are 2.15 × 10²² GeV² cm⁻⁵ and 1.53 × 10²² GeV² cm⁻⁵ for the gNFW and NFW profiles, respectively.

We fit using the χ² method to the GC excess spectrum obtained by fitting the excess spatial template (gNFW) in each energy bin independently (Section 2.3). The differential number of counts assuming a set of DM parameters is simply Equation (18) convolved with the PSF and multiplied by the instrument exposure. If we assume a specific mass, annihilation channel, and J factor, then the normalization of the DM component will be proportional to $\langle \sigma v\rangle$, which is left free in the fit. The χ² of each DM fit is

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{j}\,\displaystyle \frac{{({s}_{j}-{\lambda }_{j}(\langle \sigma v\rangle ,{\boldsymbol{p}}))}^{2}}{{\epsilon }_{j}^{2}}+{\left(\displaystyle \frac{\langle \sigma v\rangle }{\delta \langle \sigma v\rangle }\right)}^{2},\end{eqnarray} \tag{ 19 }$$

where s_j is the GC excess counts in energy bin j, λ_j is the expected photon counts for a given DM model, ${\boldsymbol{p}}$ are the assumed DM model parameters $[{m}_{\chi },{{dN}}_{\gamma }/{{dE}}_{\gamma },J \mbox{-} \mathrm{factor}]$, and _j is the statistical uncertainty on s_j from the spatial fit. The second term is the log likelihood of the prior probability on the annihilation cross-section. We conservatively assume that the prior $\langle \sigma v\rangle =0$. The uncertainty $\delta \langle \sigma v\rangle$ is derived by scanning the cusp profile along the GP; see Equations (17) and (18). If the value of the cross-section necessary to minimize the first term were large compared to the uncertainty in the denominator of the second term, then the χ² in Equation (19) would be large, which could have been interpreted as an exclusion of $\langle \sigma v\rangle =0$ prior. However, due to relatively large uncertainty on $\langle \sigma v\rangle$, the first term in Equation (19) can be minimized without a large increase in the second term, i.e., the value of $\langle \sigma v\rangle$ that minimizes the first term is consistent with the uncertainty. Based on this χ², accounting for systematic uncertainties from the GP in Section 9, we concluded that a DM interpretation of the GC excess at present is not robust.

We note that $\delta \langle \sigma v\rangle$ in Equation (19) equally accounts for systematic uncertainties that mask or induce DM-like signals. However, this was not what was seen in the control region fits, which were more like a one-sided Gaussian distribution with about half the normalizations being zero. Therefore, we chose $\delta \langle \sigma v\rangle$ by finding the 84th percentile of the fit profiles, i.e., the fractional ${b}_{\mathrm{eff}}$ for which only 16% of the fit results are larger. In the case of a Gaussian distribution this corresponds to a one-sided 1σ exclusion. The use of a nuisance parameter of this form is valid when calculating the TS since systematic uncertainties that would induce a DM-like signal, and so would reduce the TS of a DM annihilation component, are properly represented.

Then, in Section 9 we derive upper limits by requiring the DM signal to not exceed the largest value of the excess spectrum found in each energy bin under all the modeling scenarios considered. Alternatively, we can use Equation (19) to derive 95% confidence level limits based on the results of the GP scan by calculating when Δχ² = 3.84. We assume that δf_syst equally represents the amplitude of systematic uncertainties that mask or induce DM-like signals. Both versions of limits are shown in Figure 34. The limits based on the symmetric nuisance parameter from the GP scan are less constraining. However, the symmetric assumption results in limits that are likely too conservative, since a positive excess was found in the GC in every modeling scenario considered, and the vast majority of the fits in control regions along the GP found positive DM-like excesses.

Zoom In Zoom Out Reset image size