# The Galactic center excess at the highest energies:
morphology and photon-count statistics

Silvia Manconi<EMAIL_ADDRESS>Laboratoire d’Annecy-le-Vieux de
Physique Théorique (LAPTh), CNRS, USMB, F-74940 Annecy, France  Francesca
Calore  Laboratoire d’Annecy-le-Vieux de Physique Théorique (LAPTh), CNRS,
USMB, F-74940 Annecy, France  Fiorenza Donato  Department of Physics,
University of Torino, via P. Giuria, 1, 10125 Torino, Italy  Istituto
Nazionale di Fisica Nucleare, via P. Giuria, 1, 10125 Torino, Italy

###### Abstract

The nature of the GeV gamma-ray Galactic center excess (GCE) in the data of
Fermi-Large Area Telescope (LAT) is still to be unveiled. We present a new
analysis of the inner Galaxy Fermi–LAT data at energies above 10 GeV, based on
an innovative method which combines the skyFACT adaptive template fitting with
and the 1pPDF pixel-count statistics. We find a strong evidence for the GCE
also at high energies, $\sigma>5$ regardless of the GCE spatial template.
Remarkably, our fits prefer the bulge morphological model over the dark matter
one at high significance, and show no evidence for an additional dark matter
template on top of the bulge component. Through the 1pPDF analysis, we find
that the model best describing the gamma-ray data requires a smooth, diffuse
GCE following a bulge morphology, together with sub-threshold point sources.
The 1pPDF fit reconstructs a consistent population of faint point sources down
at least to $10^{-12}$ ph cm-2 s-1 . Between $10^{-12}$ ph cm-2 s-1 and
$10^{-11}$ ph cm-2 s-1 the 1pPDF measures a number of point sources
significantly higher than the ones in the Fermi 4FGL catalog. The robustness
of our results brings further support to the attempt of explaining, at least
partially, the high-energy tail of the GCE in terms of a population of point
sources, likely corresponding to millisecond pulsars.

††preprint: LAPTH-008/24

##  I Introduction

The nature of the GeV gamma-ray Galactic center excess (GCE) in the data the
he Fermi-Large Area Telescope (LAT) [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] is still
to be unveiled. This unexpected component, detected significantly on top of
the background coming from known sources, could be hinting at GeV thermal
relic dark matter particles annihilating in the Galactic dark matter halo [5,
9, 11, 12, 13, 14]. Albeit intriguing, the dark matter interpretation is
currently challenged by a growing corpus of evidence that the morphology of
the excess is better described by what expected from a population of
millisecond pulsar (MSP)-like sources in the Galactic bulge [15, 16, 17, 18].
While the GCE photon flux is peaked at about few GeV, a high energy component
extending up to tens of GeV has been reported by various studies. In
particular, the high energy tail is found to be statistically significant no
matter the nature of the spatial template used to describe the GCE emission,
i.e. either when using a dark matter-like or stellar bulge-like excess
template [8, 17, 16, 14]. The GCE high energy tail has been further
characterized in e.g. [11, 10, 19]. The spectrum was found to be consistent
among different the inner Galaxy (IG) [11], but the reconstructed properties
resulted hampered by large systematic uncertainties, due to the intrinsic
lower statistics and diffuse emission modeling. Recently, a few template
fitting analysis on updated Fermi–LAT datasets found high-energy spectra
consistent with earlier results [14, 20, 21, 22].

If confirmed, such high-energy photons could be naturally explained by the
inverse Compton scattering of electrons and positrons emitted by a population
of MSPs in the Galactic bulge [23, 24]. For this putative population of
unresolved MSPs, it would be crucial to know their phenomenological
characteristics, such as the luminosity function, so to understand how many
sources are expected to shine at different photon fluxes [25]. Population
studies often rely on observational-driven quantities, such as the source-
count distribution $dN/dS$, i.e. the distribution of the number of sources $N$
as a function of their observed gamma-ray flux $S$, where $S$ is the detected
flux integrated in a given energy range [$E_{\rm min},E_{\rm max}$]. As for
the Fermi–LAT, the number of detected sources typically decreases for
$S\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise
0.86108pt\hbox{$<$}\hss}}\lower
3.87495pt\hbox{\kern-1.90002pt$\sim$}}10^{-11}$ ph cm-2 s-1 , where the
detection efficiency drops below one [26, 27].

The pioneering work of Ref. [28] indicated how to measure the $dN/dS$ below
catalog flux threshold, and to decompose the different contributions to the
gamma-ray sky by analyzing the statistics of gamma-ray photon counts in a
pixelized sky map. To date, two implementations of the mathematical formalism
introduced by Ref. [28] have been developed and applied to the Fermi–LAT data:
the so-called Non Poissonian Template Fitting (NPTF) [29, 30] and the 1pPDF
[31, 32, 33]. The main applications of the method have led to measure the
$dN/dS$ of isotropic, extragalactic sources [31, 27, 34, 35], and to
contextually constrain dark matter annihilations at high latitudes [33].

When these methods are applied to low Galactic latitudes, the modeling for the
bright diffuse emissions could bias the results, as demonstrated in Refs. [36,
37, 38, 39], casting doubts on early evidence for unresolved point sources to
explain the GCE [29]. In [40] (PaperI in what follows) we pioneered a new
approach combining 1pPDF and the adaptive-constrained template fitting skyFACT
algorithm [41, 16]. More specifically, we used photon-count statistical
methods within the 1pPDF framework to measure faint point sources, well below
the Fermi-LAT detection threshold in the IG, between 2 and 5 GeV. We found
that all point sources and the bulge-correlated diffuse emission each
contributes O(10%) of the total IG emission, thus disclosing a potential sub-
threshold point-source contribution to the GCE at low energies. In addition,
we showed that mis-modeling of the Galactic diffuse emission causing residual
at low angular scales can lead to spurious evidence for new components, such
as point sources, and we overcame this limitation by optimizing the background
models through adaptive template fitting with skyFACT. Complementary
approaches have been recently explored, using spherical harmonics [39] and
Gaussian processes [42], as well as aiming at improving physics inputs
relevant for the computation, e.g. the reconstructed gas maps [20, 43].

At the highest energies, the nature of the GCE is even more challenging to
probe. A recent dedicated analysis of the high-energy tail of the GCE [44]
analyzed 7 years of Pass8 Fermi–LAT data with template fitting and NPTF. Using
standard template fitting techniques, the authors find a statistically
significant high-energy tail at $>10$ GeV with properties similar to the peak
GCE emission at few GeV. An emission correlated with a GCE template was
detected up to $\sim 50$ GeV, even if above about $30$ GeV the morphology of
the signal was found to be unconstrained. The NPTF analysis captured a mild
evidence for a point-source explanation of the excess at energies above $\sim
5$ GeV, although hampered by the same systematic uncertainties of the NPTF
analysis in the 2–20 GeV energy bin [29], as subsequently demonstrated in [36,
37, 38].

The goal of the current work is to assess the GCE significance and measure the
characteristics of IG gamma rays at energies larger than 10 GeV. We explore
the source-count distribution of point-like sources and its spatial
morphology. We will do so by combining state-of-the-art methods as pioneered
in PaperI. Our investigation of the GCE at high energies is novel in multiple
aspects: (i) we adopt adaptive-constrained template fitting to probe the
evidence of the GCE and reconstruct its morphology at the highest energies;
(ii) we combine this technique with pixel count statistical methods in order
to assess the role of sub-threshold point sources to the GCE at $E>10$ GeV,
while minimizing the mis-modelling of Galactic diffuse emission backgrounds;
(iii) we characterize the properties of the gamma-ray emission in the IG at
$E>10$ GeV by reconstructing the flux distribution of point sources well below
the Fermi-LAT detection threshold, and by estimating their spatial density as
compared to other regions in the sky.

The paper is organized as follows. The methodology is introduced in section
II, followed by the definition of the Fermi–LAT dataset and the region of
interest used in section III. The results of the skyFACT and 1pPDF analysis
are illustrated in sections IV and V, before concluding in section VI.
Additional systematic tests are summarized in the appendix.

##  II Methodology

The gamma-ray emission observed by Fermi–LAT is usually separated in truly
diffuse emissions and point-like (or extended) sources. Among diffuse
emissions, the one of our own Galaxy, i.e. the Galactic diffuse emission,
represents the dominant contribution to the observed photons, in particular
towards the Galactic plane. It originates from hadronic and leptonic cosmic
rays interacting with the interstellar gas and radiation fields. Other diffuse
structures extended few tens of degrees in the sky are the Fermi Bubbles [45,
46] and the Loop I [47].

Individual sources of gamma rays, connected to e.g. compact Galactic objects
such as pulsars, or distant active galaxies are collected in the source
catalogs compiled by the LAT collaboration. The last catalog, the 4FGL-DR4
[48, 49] contains 7194 sources significantly detected between energies of 50
MeV to 1 TeV with a significance larger than $5\sigma$ over a data-taking
period of about 14 years. At high latitudes, sources that are too faint to be
detected with a statistical significance of $5\sigma$ during a given survey,
together with a residual cosmic-ray contamination, make up the diffuse gamma
ray background [50] (DGRB), which exhibits anisotropies that could be used to
characterize the different populations of dim sources contributing to it [51].
These so-called unresolved sources can also considerably contribute to the
gamma-ray sky in the inner part of the Galaxy.

The methods, based on the skyFACT and the 1pPDF tools, both model diffuse and
individual sources of gamma rays using complementary approaches. They are
based on the use of templates, i.e. pixelized maps describing the energetic
characteristics of different sky components as a function of the sky
coordinates. To this, the photon-count method adds a statistical description
of bright and faint point sources. While the methods and procedures generally
follow what introduced in PaperI for the energy bin 2-5 GeV, we present here
the developments that turned out to be necessary to extend our analysis to
higher energies.

Figure 1:  The Fermi–LAT counts per pixel in the energy bin [10, 300] GeV for
CLEAN events as measured within the region of interest for the 1pPDF analysis
in the Inner Galaxy (upper left panel), Outer Galaxy (upper right) and
extragalactic sky (lower panel). All maps are in Mollweide projection and
using HEALPix resolution parameter $\kappa=7$.

###  II.1 Adaptive-constrained template fitting with SkyFACT

The skyFACT method was introduced in Ref. [41] and applied to the GCE in Ref.
[16]. By introducing a large number of parameters, and using regularization
condition for the likelihood, it accounts for intrinsic uncertainties in
spectral and spatial predictions of the different templates entering the fit
of Fermi–LAT data.

The modeling of the gamma-ray sky and statistical analysis closely follows
Ref. [16] and PaperI. The gamma-ray sky is interpreted as the combination of
several model components, each of which is characterized by a spectral and a
spatial (i.e. morphology) input template. The Galactic diffuse emission
follows the implementation in Ref. [16]: It is composed by an inverse Compton
contribution whose spectrum and morphology are computed for a typical scenario
of cosmic-ray sources and propagation parameters with the DRAGON code [52],
and a contribution from neutral pion decay split into three different
templates for more flexibility. Each gas template corresponds, spatially, to a
different gas ring (0 – 3.5 kpc, 3.5 – 6.5 kpc, 6.5 – 19 kpc), and it is the
sum of atomic and molecular hydrogen distribution within that ring, as
available in the GALPROP public release111https://galprop.stanford.edu/ [53].
The gamma-ray input spectrum of the gas component corresponds to the pion
decay spectrum, it is taken from Ref. [54], and it is the same for the three
rings. On top of the Galactic diffuse emission, the model takes contribution
from: an isotropic spatial component with the best-fit DGRB spectrum from Ref.
[50] and a uniform spatial template; all point-like and extended 4FGL sources;
the Fermi bubbles with input spectrum from Ref. [45] and a uniform geometrical
template input morphology. Finally, we consider, on top of this _background_
model, a template for the GCE.

In our setup, the baseline choice for the GCE input spectrum is a power law
with index of -2.5. We will show that the choice of the spectral index does
not affect the final results since we leave full freedom for spectral
modulation, see below. In PaperI baseline run’s configuration, the GCE input
spectral model was instead an MSP-motivated spectrum, i.e. a power-law with
exponential cutoff [55]. However, this spectrum, being zero at high energies,
cannot be modulated such to account for additional flux at $E>10$ GeV, so that
is not adapted for studying the presence of a high-energy tail of the GCE. We
discuss systematic uncertainties related to the GCE input spectrum in appendix
A. As in PaperI, we will test two different morphologies of the GCE: (i) a
dark-matter inspired template following a generalised Navarro-Frenk-White dark
matter density squared with inner slope 1.26 (NFW126); (ii) a bulge template,
made up by a boxy-bulge component from Ref. [56] and a nuclear bulge as in
Ref. [16]. While a more recent determination of the bulge from Ref. [57] has
been proposed to better fit the data [58, 59], we adopt the same bulge model
as in PaperI, for the sake of comparing fairly with our previous results. We
will nonetheless test also this new bulge model in appendix A. We refer to
Refs. [41, 16, 40] for more details about the gamma-ray model.

skyFACT implements adaptive-constrained template fitting by enabling each
model component to vary both spectrally and spatially. Starting from the input
spectra and morphologies of each model component, as mentioned above, skyFACT
minimizes a penalized likelihood composed by a standard Poissonian term and a
regularization term which is governed by hyperparameters preventing
overfitting and effectively constraining the modulation of the different
components. Modulation parameters for the spectral and spatial parts of each
component are therefore adjusted in order to maximize the likelihood. The
result of the fit for each component is therefore a modulated, output spectrum
and spatial template. How much this output differs from the input templates
depends on the freedom allowed during the fit. The minimization algorithm is
the L-BFGS-B (Limited memory BFGS with Bound constraints) one. For more
details about the technical implementation, we refer the reader to Ref. [41].
For this work, the baseline setup of hyperparameters follows run5 of Ref. [41]
for all components but the additional GCE. Similarly to what done in PaperI,
the GCE spectrum is allowed to fully vary in each energy bin independently,
while its morphology is fixed to the input template, so that there is no
additional freedom enabled for spatial modulation parameters.

The statistical framework implemented by skyFACT is purely frequentist. The
output of the fit is the maximum likelihood value of the specific model
considered. As such, we cannot properly do model comparison – which would
require a Bayesian approach – within skyFACT. On the other hand, we can test
the preference for additional model components by considering a nested model,
i.e. a model which includes a new sky contribution (e.g. GCE) on top of
another, less complex, model. The augmented model reduces to the other (i.e.
to the null hypothesis) for some choice of the additional component model
parameters. For nested models the calculation of the evidence for any
additional component is well defined, and follows the modified
$\delta-\chi^{2}$ statistics as described in Ref. [16]. More specifically, we
start by fitting the gamma-ray sky without a GCE diffuse template (sF-noGCE).
We proceed adding a GCE diffuse template following the bulge morphology (sF-B)
or the NFW126 (sF-NFW126) morphology, with spectral properties defined above.
Finally, we test the evidence for a NFW126 on top of a bulge-GCE diffuse
template, and vice versa through nested model comparison.

Errors are estimated from computing an approximated covariance matrix as the
inverse of the Fisher information matrix from mock data of the best-fit model.
See details in Ref. [41].

We notice that the scope of this work is the robust determination of the IG
source-count above 10 GeV. To this end, as demonstrated in PaperI, we do not
need to achieve a perfect description of the gamma-ray sky, rather to reduce,
effectively, the fit residuals. skyFACT can satisfactorily reach this goal, by
optimizing in a data-driven way the gamma-ray diffuse model. While a better
implementation of the Galactic emission is indeed possible (see e.g. [20]),
deriving the best Galactic diffuse gamma-ray model is beyond our scope. We
therefore refrain from making model comparison statements, which are in any
case prevented by the limitations of the maximum likelihood approach.

###  II.2 Photon count statistics technique

Our goal is to measure the flux distribution of point sources in the IG at
high energies, and compare their density in different sky regions. The 1pPDF
method [31, 32, 33] measures this quantity by a statistical analysis of the
probability distribution $p_{k}^{(p)}$ of the photon counts $k^{(p)}$ in each
pixel $p$ of a pixelized map. The core of the method is based on the fact that
different classes of photon sources are expected to contribute to the
$p_{k}^{(p)}$ with different statistics. The $p_{k}^{(p)}$ of truly diffuse,
isotropic emissions follows a Poissonian distribution. Non-Poissonian
contributions from point sources and other complex diffuse structures can
alter the probability distribution $p_{k}$, permitting to investigate these
components through the 1pPDF of the observed gamma rays.

By exploiting the statistics of photon counts in the data, the 1pPDF extracts
the $dN/dS$ in relation to the number of $k$-photon sources $x^{(p)}_{k}$ in
each pixel, see [31] for the details of the mathematical formulation. The
population of point sources is modeled through a $dN/dS$ shaped as a multiple
broken power law (MBPL):

$\frac{\mathrm{d}N}{\mathrm{d}S}=A_{\rm
S}\cdot\begin{cases}\left(\frac{S}{S_{0}}\right)^{-n_{1}}\;\;\;\;S>S_{\mathrm{b}1}\,;\\\
\left(\frac{S_{\mathrm{b}1}}{S_{0}}\right)^{-n_{1}+n_{2}}\left(\frac{S}{S_{0}}\right)^{-n_{2}}\;\;\;\;S_{\mathrm{b}2}<S\leq
S_{\mathrm{b}1}\,;\\\
\mathmakebox[\widthof{{}\left(\frac{S_{\mathrm{b}1}}{S_{0}}\right){}}][c]{\vdots}&\\\
\left(\frac{S_{\mathrm{b}1}}{S_{0}}\right)^{-n_{1}+n_{2}}\left(\frac{S_{\mathrm{b}2}}{S_{0}}\right)^{-n_{2}+n_{3}}\cdots\
\left(\frac{S}{S_{0}}\right)^{-n_{N_{\mathrm{b}}+1}}\\\ \hskip
108.12054ptS\leq S_{\mathrm{b}N_{\mathrm{b}}},\\\ \end{cases}$ (1)

where the free parameters are the normalization $A_{\rm S}$, the position of
the flux break, the indices $n_{i}$ of the broken power law, and we fix
$S_{0}=5\cdot 10^{-9}$ ph cm-2 s-1 . We utilize a MBPL with two free breaks
$N_{b}=2$, by sampling directly the parameters in Eq. (1). Since we do not
mask the bright sources, the measured $dN/dS$ is expected to follow the one of
Fermi–LAT detected sources in the bright regime. Typically, the 1pPDF measures
point sources down to (about one order of magnitude) lower fluxes [31, 32,
40]. We remind that the 1pPDF method in its current implementation does not
allow to include in the statistical fit the spatial distribution of the point
source component. Thus, the measured point sources are meant to be isotropic
in the considered region of interest. In PaperI, we collected insights on the
spatial distribution of point sources by subdividing our ROI in sub-regions,
and computing the source density as a function of the position in the sky.
Given the limited statistics of the gamma-ray sky at energies larger than 10
GeV, we anticipate that there would not be enough photons to constrain
multiple source-count distributions in different sub-regions.

Similarly to PaperI, we consider within the 1pPDF-fit the following model
components for the gamma-ray sky: (i) A population of isotropic point sources,
characterized by their source-count distribution, (ii) an isotropic diffuse
emission (with free normalization), (iii) the Galactic diffuse emission
template (as derived by skyFACT or from other existing models), and (iv) a GCE
smooth template as derived by skyFACT and following the bulge or DM
morphology, both with free normalization. We give more details on each
component in what follows. The number of counts in each map pixel
$x^{(p)}_{\mathrm{diff}}$ coming from diffuse templates (ii)-(iv) is defined
as:

$x_{\mathrm{diff}}^{(p)}=A_{\mathrm{gal}}x_{\mathrm{gal}}^{(p)}+A_{\mathrm{GCE}}x_{\mathrm{GCE}}^{(p)}+\frac{x_{\mathrm{iso}}^{(p)}}{F_{\mathrm{iso}}}F_{\mathrm{iso}}\,.$
(2)

The quantity $F_{\mathrm{iso}}$ represents the integral flux of the isotropic
diffuse emission, which we use directly as a sampling parameter in the fit.
The other diffuse templates describe the Galactic diffuse emission and a GCE
smooth template. As for the GCE smooth template, they correspond to the best
fits obtained within each skyFACT run (see Sec. II.1) for the bulge and DM
components. We allow for a rescaling factor $A_{\mathrm{GCE}}$ relative to the
best-fit normalization found by skyFACT. Similarly, when we employ the
Galactic diffuse emission templates as optimized by skyFACT, we allow for a
$A_{\mathrm{gal}}$ normalization factor with respect to the best fit. To
demonstrate the effect of diffuse emission model systematics, we employ other
widely used templates, also entering with a $A_{\mathrm{gal}}$ normalization
factor. Specifically, the official spatial and spectral template released by
the Fermi–LAT Collaboration for Pass 8 data (Official P8) (gll_iem_v06.fits
[60]), and the model labeled A (modA) optimized for the study of the DGRB in
Ref. [50] are tested. We note that in modA the Fermi bubbles are not modeled.

We fit the gamma-ray sky using the pixel-dependent likelihood function
$\mathcal{L}({\bf\Theta})$ as defined within the L2 method in Ref. [31], thus
taking into account the spatial morphology of the diffuse templates. The free
parameters ${\bf\Theta}$, together with their prior intervals are summarized
in the appendix B. The posterior distribution is defined as
$P({\bf\Theta})=\mathcal{L}({\bf\Theta})\pi({\bf\Theta})/\mathcal{Z}$, where
$\pi({\bf\Theta})$ is the prior and $\mathcal{Z}$ is the Bayesian evidence.
The $P({\bf\Theta})$ is sampled using MultiNest [61] in its standard
configuration, setting 800 live points with a tolerance criterion of 0.2. To
get prior-independent frequentist maximum likelihood parameter estimates [62],
we build one-dimensional profile likelihood functions for each parameter using
the final posterior sample. To perform Bayesian model comparison, the Bayes
factors are computed by using the nested sampling global log-evidence
$\ln(\mathcal{Z})$ provided by MultiNest.

###  II.3 Combining SkyFACT and photon-count statistics insights

Our hybrid approach, already presented in PaperI, relies on a two-step
procedure. First, we perform an adaptive-constrained template fit of the
gamma-ray IG with skyFACT. Second, for each skyFACT run, the corresponding
skyFACT-optimized output models are used as input for the 1pPDF fit. By using
the Galactic diffuse emission models as optimized through a skyFACT fit, in
PaperI, indeed, we demonstrated that the residuals are reduced and the results
of the 1pPDF are more stable against diffuse mismodeling. We adopt the
skyFACT-optimized components for the Galactic emission, and also for the GCE
template, modeled either as the bulge or NFW126 model. We then proceed with
the 1pPDF fit to Fermi–LAT data using the different skyFACT-optimized Galactic
diffuse emission models, the GCE diffuse templates, consistently for each
case. We extract the Bayesian evidence $\mathcal{Z}$ and compute the Bayes
factors between each model $i$ and $j$ as
$B_{ij}=\exp(\ln\mathcal{Z}_{i}-\ln\mathcal{Z}_{j})$, since the compared
models are not nested anymore.

The interpretation of the results, and specifically the statistical statements
about the GCE morphology are thus based on (i) the nested model comparison
done within the skyFACT optimization and (ii) the Bayesian model comparison
done fitting again the skyFACT-optimized templates with the 1pPDF method,
which adds the modeling of faint point sources. Furthermore, the measure of
the $dN/dS$ is extracted within the main ROI, and its characteristics compared
to the source count distribution found at high latitudes and in the outer
Galactic disk.

##  III Fermi–LAT data

Consistently with the dataset analyzed in PaperI, we consider 639 weeks of
Fermi–LAT data [63] until 2020-08-27 and apply standard quality cuts. For the
skyFACT fit, we employ all FRONT+BACK P8R3 ULTRACLEANVETO events (evtype=3) in
order to maximize the statistics against the huge number of free parameters in
the fit. Our main region of interest (ROI) is the so-called IG defined as
40${}^{\circ}\times 40^{\circ}$ around the Galactic center. We binned the map
into 0.5∘ pixels. No mask is applied neither on point sources, nor on the
Galactic plane. The full energy range (0.2 – 500 GeV) is divided into 30
logarithmically-spaced bins. The fit is performed either in the energy range
$0.3-300$ GeV (27 bins), similarly to PaperI, or in the high-energy range 10 –
300 GeV (13 bins).

The benchmark choice for the 1pPDF fit is to consider all P8R3 CLEAN events
and focusing in a single bin covering the 10–300 GeV energy range. Given the
limited number of photons in the considered high-energy range, this choice
increases the statistics of photon counts with which the 1pPDF can measure the
$dN/dS$ down to low fluxes. Since the point spread function (PSF) of Fermi–LAT
at these energies is expected to be of the order of 0.1 – 0.5∘, i.e. smaller
or comparable to the pixel size, we consider all PSF quartile events. The
1pPDF performance has been tested using CLEAN Fermi–LAT high-energy events
elsewhere [27]. Even less stringent quality cuts (SOURCE events) have been
used in past analysis of the IG to increase the sensitivity of photon-count
statistical methods [44]. Nevertheless, to test the stability of our
conclusions against the event selection, we discuss results using P8R3
ULTRACLEANVETO events with and without further selecting the quartile of best
angular resolution (PSF3), see appendix B. All the components entering the
1pPDF fits are corrected for the PSF effect as detailed in Ref. [31]. The
photon events are binned spatially using the HEALPix equal-area pixelation
scheme [64] with resolution parameter $\kappa=7$ [64].

Different ROIs in the sky are considered in this work. The IG ROI is shared by
both the skyFACT and 1pPDF analysis. Cuts at latitudes $|b|>0.5^{\circ}$ or
$1^{\circ}$ are introduced to check the stability of 1pPDF results. We note
that, in an effort to maximize the statistics for the high-energy fit, we
increased the ROI for the 1pPDF analysis with respect to the $20^{\circ}\times
20^{\circ}$ IG region used in PaperI. We tested that the main results are
unchanged when reducing the ROI to the one of PaperI. As in PaperI, two
control regions are introduced to compare the source density with respect to
the IG: the outer Galaxy (OG, $|b|<20^{\circ}$, $60^{\circ}<|l|<90^{\circ}$,
cutting the inner $|b|<1^{\circ}$) and the extragalactic region (EG,
$|b|>40^{\circ}$, $|l|>90^{\circ}$). Figure 1 shows the Fermi–LAT data counts
map for the IG, OG, EG ROI (left to right) in the energy bin 10–300 GeV. The
total number of counts within each region for CLEAN events is 92529 (IG),
50619 (OG) and 52324 (EG).

##  IV Results: excess significance and morphology with skyFACT

In this section, we discuss the results of the fits to the gamma-ray sky with
skyFACT with the aim of determining the excess significance and its morphology
at high energies.

First, we run the fit in the OG ROI in the 10 – 300 GeV energy range to serve
as a control region. We check the evidence for an additional GCE component on
top on the SF-noGCE model. For the GCE, we consider either an NFW126 or a
Bulge template. In none of the two cases, we find significant evidence for the
GCE ($\sigma\sim 0$). This result is consistent with the fit in the full
energy range, 0.3 – 300 GeV, performed in PaperI. Therefore, the baseline
skyFACT model in the OG results from the optimization of the sky components
without GCE and it will be used as input for the subsequent 1pPDF fit.

Second, we analyze the IG ROI. We first consider the full energy range, 0.3 –
300 GeV, and test the evidence for an additional GCE component, as explained
in section II.1. The GCE spatial template is fixed and assumed to follow
either NFW126 or Bulge. The input spectrum is a power law with index -2.5, and
it is allowed to freely vary energy bin by energy bin. No matter what GCE
template is considered, we found strong evidence for this additional
component: $\sigma\sim 13$ (9.9) for the additional Bulge (NFW126) template.
We then check the preference for Bulge on top of NFW126 and, viceversa, of
NFW126 on top of Bulge. Analogously to PaperI, we find that it is necessary to
add Bulge on top of NFW126 ($\sigma=$ 9.1), while there is no need for an
additional NFW126 on top of the Bulge ($\sigma<4$). We notice that in the full
energy range and adopting a power-law input spectrum for the GCE, the evidence
for the GCE and the preference of the Bulge on top of NFW126 is consistent
with the results found in PaperI for a MSP-like input spectrum, i.e. a power
law (index -1.46) with exponential cutoff at 3.6 GeV, see also discussion in
appendix A. We can understand these results by considering that, in the full
energy range, the fit is driven by the low-energy data (1 – 5 GeV). Therefore
the possible high-energy tail of the excess, if there, does not strongly
weight in determining the evidence for the GCE and Bulge. In figure 2, we
report the best-fit spectra of the SF-B run (i.e. including the Bulge GCE) in
the full energy range 0.3 – 300 GeV. The best-fit spectrum for the Bulge (run
with power-law input spectrum) is consistent with the one in PaperI (run using
an input MSP-like spectrum for the Bulge) for energies below 10 GeV, while at
higher energies the Bulge shows significant emission, see also discussion in
appendix A.

Since the best-fit results in the full energy range are driven by the low-
energy data, in order to properly test the presence of the GCE at higher
energies we need to run the fit in a narrower energy range and evaluate the
significance of additional components therein. We therefore consider the 10 –
300 GeV energy range with same setup as above for the different sky
components. We summarize the results in table 1: The first result is that we
find strong evidence for the GCE (either modeled with a Bulge or an NFW126
template) also at high energies, $\sigma>5$ regardless of the GCE spatial
template. Secondly, our fit prefers the Bulge model over the NFW126 at high
significance ($\sigma=5.5$). While the evidence is lower than the one found in
the full energy range, it still strongly indicates the need of the GCE when
fitting the high-energy Fermi sky. Similarly, at high energies, there is no
evidence for an additional NFW126 template on top of the Bulge component
($\sigma=0.8$). The integrated flux of the Bulge component in this case is
about a factor of two higher than the best-fit integrated flux of the Bulge
above 10 GeV from the fit in the full energy range. We will demonstrate that
this difference does not impact the 1pPDF results which are stable against
variation of the overall Bulge input intensity. In appendix A, we discuss
further systematics that may affect the evidence of the GCE in the high-energy
interval, especially the modeling of the Fermi bubbles, and show that this
evidence is robust against them.

We remind the reader that our scope is not to find the overall best-fit gamma-
ray sky model, rather to provide an input model for the 1pPDF analysis which
maximally reduces residuals in an effective way and allows us to derive robust
conclusions on the source-count distribution. The skyFACT fit results for the
SF-B model in the 10 – 300 GeV energy interval will be used as input for the
1pPDF analysis in section V. The baseline templates are shown in figure 3.

Table 1: Excess evidence in the IG. Results of the skyFACT gamma-ray fit in the IG in the high-energy interval, showing the evidence (in units of $\sigma$) of the GCE itself, and of the Bulge and NFW126 templates on top of the baseline background model. In all runs, the GCE morphology is fixed, while the spectrum is allowed to freely vary energy bin by energy bin. The input GCE spectrum is a power law with index -2.5. | GCE Bulge (NFW126) | Bulge over NFW126 | NFW126 over Bulge
---|---|---|---
Evidence ($\sigma$) | 8.1 (5.6) | 5.5 | 0.8
Figure 2:  skyFACT best-fit spectral components. Results of the IG analysis in
the 0.3 – 300 GeV energy range, for the preferred SF-B model composed by the
background model, and an additional GCE Bulge component (i.e. BB+NB). See
section II.1 for more details on the model components.

Figure 3:  The skyFACT-optimized templates for the Galactic diffuse emission
(left), the sum of nuclear and boxy stellar bulge (the Bulge, center), and
NFW126 dark matter as integrated in the 10 – 300 GeV energy bin. Note the
different scales in each panel.

##  V Results: source–count distribution

In this section we discuss the results for the source–count distribution
$dN/dS$ as obtained analyzing Fermi–LAT data with the 1pPDF.

Figure 4:  Extragalactic and Outer Galaxy $dN/dS$. Source-count distribution
of the EG ROI (left panel) and of the OG (right panel) as obtained from the
1pPDF analysis. The $dN/dS$ result using the Official P8 (blue lines) model,
modA (magenta lines) and the skyFACT Galactic diffuse emission model without
any component modeling the GCE (yellow line, only right panel) are reported.
The $1\sigma$ uncertainty bands are indicated with colored areas. The count
distribution of 4FGL sources is illustrated with black points for comparison
in the bright flux regime.

###  V.1 Extragalactic and Outer Galaxy sky regions

The result for the $dN/dS$ of the EG (OG) sky is illustrated in the left
(right) panel of figure 4. We report the result obtained with the 1pPDF
analysis of CLEAN all PSF events in the 10–300 GeV energy range when using the
Official Pass8 diffuse emission template (blue) and the modA (magenta). As for
the OG ROI, we also show the results obtained using an optimized skyFACT
diffuse template (yellow) without including a GCE template. The best fit
results are reported together with the $1\sigma$ uncertainty bands. The
source-count distribution from the sources in the 4FGL is shown for reference
with black points, where uncertainties include statistical uncertainties on
the source count only. As expected, the 1pPDF recovers the $dN/dS$ of bright
sources in agreement with the 4FGL number counts down to fluxes of about
$S\sim 10^{-11}$ ph cm-2 s-1 . For lower fluxes, the $dN/dS$ of 4FGL sources
decreases quickly, while the 1pPDF measures the collective contribution of
faint point sources down to at least one order of magnitude at $S\sim
10^{-12}$ ph cm-2 s-1 . At even lower fluxes, the $1\sigma$ uncertainty band
increases significantly, and the measured $dN/dS$ becomes fully degenerate
with the diffuse isotropic emission (in other words, the position of the
second break of the $dN/dS$ and the normalization of the diffuse isotropic
emission are highly correlated). This measurement of the $dN/dS$ is robust
against possible systematics coming from the $dN/dS$ parametrisation as MBPL
or using the Hybrid approach, as well as the event selection, see the appendix
B.

In the EG ROI, consistent $dN/dS$ are measured when using the Official Pass8
and the modA diffuse models, confirming the robustness of the results against
diffuse emission mismodeling at high latitudes, similarly to what observed in
the 2–5 GeV energy range in PaperI. The measured $dN/dS$ is in overall
agreement with previous 1pPDF results obtained within similar event selection,
energy bin and region of interest in the high latitude sky [32]. As expected,
it becomes unconstrained below about $6\cdot 10^{-13}$ ph cm-2 s-1 .

Using the Official P8 diffuse model and the skyFACT optimized one, the $dN/dS$
of the OG is found to be well described by a power law of index
$n_{2}=2.05^{+0.17}_{-0.08}$, from $5\cdot 10^{-10}$ ph cm-2 s-1 down to
$5\cdot 10^{-12}$ ph cm-2 s-1 , where the uncertainty band opens by few orders
of magnitude. These $dN/dS$ results are compatible with the source counts of
bright sources within the 4FGL catalog, and reveal unresolved sources down to
at least $5\cdot 10^{-12}$ ph cm-2 s-1 . Instead, we observe effects of
mismodeling the diffuse emission when using modA, while the results using the
skyFACT and the Official P8 model are consistent within statistical
uncertainties. Similarly to what observed in PaperI, in this scenario the
1pPDF looses its sensitivity and fails to measure faint sources at fluxes
lower than about $10^{-11}$ ph cm-2 s-1 . The $dN/dS$ measured by the 1pPDF
with modA suggests that residuals in the fitted gamma-ray sky are interpreted
as bright point sources around the sensitivity threshold of the 4FGL catalog.
The measure of the $dN/dS$ of gamma-ray sources in the OG ROI in the 10–300
GeV energy range is a novel result of this work.

Figure 5:  Inner Galaxy $dN/dS$. Source-count distribution of the IG ROI
measured by the 1pPDF using different skyFACT-optimized diffuse models (left
panel, for $|b|>1^{\circ}$) and compared to non-optimized models from the
literature (right panel, for $|b|>0.5^{\circ}$). The $1\sigma$ uncertainty
bands are indicated with colored areas. The count distribution of 4FGL sources
is illustrated with black points for comparison in the bright flux regime.
Gray points represent the count distribution of 4FGL sources without any
analysis flag, see [65].

###  V.2 Inner Galaxy

We proceed to measure the source count distribution of bright and faint point
sources within the IG ROI in the 10–300 GeV energy bin, and considering
different Galactic diffuse emission models and GCE diffuse models as obtained
in the dedicated skyFACT fit described in Sec. II. When using the skyFACT-
optimized diffuse models, we consistently include in the 1pPDF fit the GCE
diffuse model following a NFW126 or a bulge morphology, with a free overall
normalization, apart for the noGCE setup. This permits us to evaluate again,
this time with Bayesian statistics as detailed in Sec. II.3, which model is
preferred when fitting the IG including also faint point sources. The results
for the 1pPDF fits using the diffuse models obtained within the SF-noGCE, SF-
NFW126 and SF-B fits are summarised in Table 2 for the 1pPDF fits of CLEAN
allPSF data. The columns illustrate the latitude cut, the Bayesian log-
evidence, the point source fluxes for $S<10^{-10}$ ph cm-2 s-1 and the
normalization of the diffuse GCE for the different configuration explored.
Regardless of the latitude cut, we find that the model best describing the
gamma ray data is the one in which the Galactic diffuse emission corresponds
to the skyFACT fit and that also include a smooth, diffuse GCE following a
bulge morphology with Bayes factor $\ln B_{ij}\gtrsim 30$. Adding a diffuse
GCE template is also preferred with respect to the SF-noGCE case with similar
Bayes factor. Both the Galactic diffuse emission and the GCE smooth templates
are found in the 1pPDF fit with normalizations consistent with one. We tested
that similar results are obtained using UCV allPSF and PSF3 data selections.
This corroborates the result obtained within skyFACT: The IG gamma-ray sky at
energies larger than 10 GeV is best described when a smooth GCE following a
bulge morphology is included in the model, also when modeling the faint point
sources.

The PS flux is found to be consistent among the different optimized Galactic
diffuse models employed, and slightly higher when cutting only the inner
$0.5^{\circ}$, see next section for further discussion on the source density.
By integrating all the PS fluxes below $S<10^{-10}$ ph cm-2 s-1 (and thus in
an intermediate regime between bright and faint), we obtain a total PS flux
which is comparable to the total flux measured by skyFACT in the full IG ROI
within the fit for $E>10$ GeV, equal to $2.37\cdot 10^{-8}$ ph cm-2 s-1 .

Table 2: Results for the 1pPDF analysis of the IG LAT data, CLEAN, all PSF. First three columns: setup of the analysis and latitude mask of the IG. The $\ln(\mathcal{Z})$ is the nested sampling global log-evidence extracted from Multinest. Last two columns: flux from IG point sources (in units of $10^{-8}$ ph cm-2 s-1 , and for $S<10^{-10}$ ph cm-2 s-1 ), and normalization of smooth GCE template in the 1pPDF fit when relevant. Description |  $|b|$ cut [∘] | $\ln(\mathcal{Z})$ | PS flux | $A_{\rm B/NFW126}$
---|---|---|---|---
No GCE | 0.5 | $-19056.81$ | $3.42^{+0.78}_{-1.58}$ | -
NFW126 | 0.5 | $-19060.14$ | $3.76^{+0.41}_{-1.75}$ | $0.84^{+0.10}_{-0.16}$
Bulge | 0.5 | $-19027.22$ | $3.94^{+3.89}_{-1.98}$ | $0.97^{+0.06}_{-0.08}$
No GCE | 1 | $-18085.2$ | $2.75^{+0.76}_{-1.39}$ | -
NFW126 | 1 | $-18072.6$ | $2.84^{+0.7}_{-1.3}$ | $1.45^{+0.05}_{-0.15}$
Bulge | 1 | $-18054.34$ | $2.68^{+1.07}_{-1.08}$ | $1.02^{+0.08}_{-0.04}$

In addition, we find that a consistent population of faint point sources is
reconstructed no matter the skyFACT-optimized Galactic diffuse emission model.
This is illustrated in the left panel of figure 5. The $dN/dS$ as measured by
the 1pPDF and cutting the inner $|b|<1^{\circ}$ is reported for the SF-noGCE,
SF-NFW126, and SF-B setups with different colors, together with the $1\sigma$
uncertainty band. No matter the setup of diffuse templates used, the 1pPDF
consistently measures faint point sources down at least to $10^{-12}$ ph cm-2
s-1 . At brighter fluxes and up to $10^{-11}$ ph cm-2 s-1 , the reconstructed
$dN/dS$ is found compatible with the source count in the 4FGL, in particular
when flagged sources are excluded. Between $10^{-12}$ ph cm-2 s-1 and
$10^{-11}$ ph cm-2 s-1 the 1pPDF measures a number of point sources
significantly higher than the ones in the 4FGL. Between $10^{-9}$ ph cm-2 s-1
and $10^{-12}$ ph cm-2 s-1 , the $dN/dS$ follows a power law with index
$n_{2}\sim 2$. At fluxes lower than $10^{-12}$ ph cm-2 s-1 the uncertainty
band opens significantly and the second break in the $dN/dS$ becomes fully
degenerate with the isotropic diffuse background normalization.

In the right panel of figure 5 we illustrate the effect of Galactic diffuse
emission mis-modeling for $|b|<0.5^{\circ}$. When using non-optimized Galactic
diffuse emission models such as the Official P8 (blue line) or the modA
(magenta), the 1pPDF method is not able to reconstruct bright point sources
even at high fluxes of $10^{-10}$ ph cm-2 s-1 , and spurious peaks in the
$dN/dS$ are found with high significance. We attribute this effect to the
residuals left when analyzing gamma rays at low latitudes, which are
misattributed by the 1pPDF to a bright population of point sources.

###  V.3 Spatial distribution: Source density

In order to characterize the spatial distribution of the point source
population measured by the 1pPDF we compute the source density in the
different ROIs analyzed. We remind that, in the current implementation, the
1pPDF method measures an isotropic population of point sources in the region
of interest, and thus is not sensitive to the pixel-by-pixel spatial
distribution of sources. Thus the measured $dN/dS$, and in turn the source
density are to be considered as average quantities within the ROI. We compute
the source density by integrating the measured $dN/dS$ (and the corresponding
$1\sigma$ uncertainty) in the flux interval [$3\cdot 10^{-12},1\cdot
10^{-11}$] ph cm-2 s-1. This flux interval corresponds to the regime in which
the 4FGL is incomplete for the source counts, and thus characterizes the
unresolved point sources collectively measured by the 1pPDF. The source
density obtained integrating the $dN/dS$ as measured in the EG ROI is of $\sim
0.12^{+0.02}_{-0.01}$ sources/deg2. The OG source density is found to be
slightly higher, at the level of $0.20^{+0.09}_{-0.03}$ sources/deg2,
suggesting that Galactic sources contribute to the $dN/dS$ in addition to the
extragalactic sources measured in the EG.

The source density in the full IG ROI is $0.31^{+0.11}_{-0.09}$ sources/deg2
when cutting the inner $|b|<1$ degree, and $0.28^{+0.10}_{-0.09}$ sources/deg2
when cutting the inner $|b|<0.5$ degree. Both results are compatible within
errors, and reveal a source density higher with respect to the EG and OG.
However, the reduced photon statistics in the 10–300 GeV energy bin prevents
to robustly measure the radial and longitude profiles of point sources using
multiple regions as done in PaperI. We thus provide the source density within
the North and South hemispheres, corresponding to $0.35^{+0.20}_{-0.24}$
sources/deg2 and $0.33^{+0.22}_{-0.29}$ sources/deg2 respectively, as well as
for positive and negative longitudes, equal to $0.47^{+0.30}_{-0.32}$
sources/deg2 and $0.18^{+0.15}_{-0.15}$ sources/deg2. The source density is
found to be slightly higher than the OG and the average IG level within
positive longitudes. We note that, if the GCE consists of MSPs within the
stellar bulge of the Milky Way, and if they contribute to the unresolved point
sources in the considered flux interval, an asymmetry among negative and
positive longitudes is indeed expected, being the stellar bulge more luminous
at positive longitudes, see e.g. figure 3. However, we leave the
interpretation of these measurements to forthcoming work.

##  VI Conclusion and Outlook

Figure 6:  The GCE energy spectrum is illustrated in linear scale to highlight
the $E>10$ GeV tail. Our results obtained with skyFACT for energies larger
than $10$ GeV are shown in magenta. The measured spectra from the template-
based analysis of Ref. [11] is reported with blue points, where the shaded
band encloses the systematic uncertainties on the GCE spectra when varying the
Galactic diffuse emission modeling. Model interpretations assuming dark matter
annihilations (Ref. [21], yellow dot-dashed line, rescaled by a factor 1.5) or
MSP prompt (purple dotted) plus IC emission (dashed) from Ref. [24] are
overlaid for comparison.

In this work, we presented a new study of the IG Fermi–LAT data at energies
above 10 GeV with the goal of investigating the robustness of the high-energy
tail of the GCE, and assessing the role of sub-threshold point sources. Our
analysis is based on an innovative method which combines adaptive template
fitting and pixel count statistical methods, while minimizing the mis-
modelling of Galactic diffuse emission backgrounds. The present paper extends
our previous results of PaperI focused on the 2-5 GeV energy bin, where the
GCE is prominent.

We performed adaptive-constrained template fits with skyFACT, firstly in the
full energy 0.3 – 300 GeV, and then in the high-energy range 10 – 300 GeV. The
gamma-ray sky was interpreted as the combination of several components, each
one modeled by a specific spectral and morphological input template. In
particular, we considered a NFW126 or a Bulge template for the GCE. We found a
strong evidence for the GCE also at high energies, $\sigma>5$ regardless of
the GCE spatial template. Remarkably, our fit preferred the Bulge model over
the NFW126 at high significance ($\sigma=5.5$) and showed no evidence for an
additional NFW126 template on top of the Bulge component ($\sigma=0.8$).

The skyFACT fit results, which maximally reduce residuals in an effective way,
were then employed as input model for the 1pPDF analysis. After measuring the
$dN/dS$ in the OG and EG control regions, we proceeded to measure the source
count distribution of bright and faint point sources within the IG ROI in the
10–300 GeV energy bin, and considered different Galactic diffuse emission
models and GCE diffuse models as obtained in the dedicated skyFACT fit. We
found that the model best describing the gamma-ray data is the one that
contains the skyFACT fit for the Galactic diffuse emission with  the inclusion
of a smooth, diffuse GCE following a bulge morphology. Adding a diffuse GCE
template is also preferred with respect to the SF-noGCE case. Both the
Galactic diffuse emission and the GCE smooth templates are found in the 1pPDF
fit with normalizations consistent with one. Notwithstanding the low
statistics characterizing the emission at energies $>10$ GeV, the IG gamma-ray
sky at energies larger than 10 GeV is best described when a smooth GCE
following a bulge morphology is included in the model, also when modeling the
faint point sources. A further result reached by the 1pPDF is that a
consistent population of faint point sources is reconstructed. No matter the
skyFACT-optimized Galactic diffuse emission model, the 1pPDF consistently
measures faint point sources down at least to $10^{-12}$ ph cm-2 s-1 . Between
$10^{-12}$ ph cm-2 s-1 and $10^{-11}$ ph cm-2 s-1 the 1pPDF measures a number
of point sources significantly higher than the ones in the 4FGL.

In the current implementation, the 1pPDF method is not sensitive to the pixel-
by-pixel spatial distribution of sources. In an attempt to characterize the
spatial distribution of the point-source population measured by the 1pPDF, we
computed the source density in the EG, OG and IG ROIs by integrating the
measured $dN/dS$ in the flux interval [$3\cdot 10^{-12},1\cdot 10^{-11}$] ph
cm-2 s-1. We found an IG source density higher with respect to the EG and OG.
Similar source densities are found within the North and South hemispheres,
while a higher density was found for positive longitudes with respect to
negative ones. If GCE were explained by point sources such as MSPs, one would
expect an asymmetry among negative and positive longitudes, being the stellar
bulge indeed more luminous at positive longitudes.

The spectrum of the GCE as found in our skyFACT fit to the Fermi–LAT data at
energies above 10 GeV is reported in figure 6, along with earlier results. The
linear $x$-axis emphasizes the high-energy tail of the GCE, where we confirm
and support the existence of a significant high-energy tail in the GCE
spectrum at energies of $E>10$ GeV, consistently with the results from
template fitting by e.g. Refs. [11]. For illustrative purpose, we overlay two
model interpretations picked from the literature. Dark matter annihilation of
weakly interacting massive particles with $m_{\rm DM}=40$ GeV in the
$b\bar{b}$ channel [21] (yellow dot-dashed line, rescaled by a factor 1.5 to
match the normalization of GCE spectrum found in Ref.[11]) can explain the
peak of the GCE at around few GeV, but fails to explain the high-energy tail.
Annihilation into a combination of different final states as suggested by
specific dark matter models could explain the spectrum until few tens of GeV,
but they are subject by other constraints, see e.g. [66, 67, 68, 69].
Conversely, the cumulative flux from MSPs in the Galactic bulge, including the
prompt (purple dotted) and inverse Compton (IC, dashed) emissions as modeled
by Ref. [24] have a spectrum compatible with both the GCE peak and the high-
energy tail, depending on the spectral properties of the $e^{\pm}$ population
emitted [23].

The robustness of our results confirming the existence of a GCE at high
energies, demonstrating the preference of the Fermi–LAT data for a Bulge
morphology and measuring a population of point sources below the detector
threshold, motivates further studies towards the understanding of the GCE in
terms of point sources and specifically of MSPs. This exploration, along with
other viable and complementary interpretations, are left for future study.

###### Acknowledgements.

We warmly acknowledge C. Eckner for discussion on skyFACT performance and
careful reading of the manuscript, and D. Malyshev and N. Rodd for careful
reading of the manuscript and useful comments. SM acknowledges the European
Union’s Horizon Europe research and innovation program for support under the
Marie Sklodowska-Curie Action HE MSCA PF–2021, grant agreement No.10106280,
project VerSi. FC acknowledges support by the “Agence Nationale de la
Recherche”, grant n. ANR-19-CE31-0005-01. F.D. acknowledges the support of the
Research grant TAsP (Theoretical Astroparticle Physics) funded by Istituto
Nazionale di Fisica Nucleare.

## References

  * Goodenough and Hooper [2009] Lisa Goodenough and Dan Hooper, “Possible Evidence For Dark Matter Annihilation In The Inner Milky Way From The Fermi Gamma Ray Space Telescope,” arXiv preprint (2009), arXiv:0910.2998 [hep-ph] .
  * Vitale and Morselli [2009] Vincenzo Vitale and Aldo Morselli (Fermi-LAT), “Indirect Search for Dark Matter from the center of the Milky Way with the Fermi-Large Area Telescope,” in  _2009 Fermi Symposium_ (2009) arXiv:0912.3828 [astro-ph.HE] .
  * Hooper and Goodenough [2011] Dan Hooper and Lisa Goodenough, “Dark Matter Annihilation in The Galactic Center As Seen by the Fermi Gamma Ray Space Telescope,” Phys. Lett. B 697, 412–428 (2011), arXiv:1010.2752 [hep-ph] .
  * Hooper and Linden [2011] Dan Hooper and Tim Linden, “On The Origin Of The Gamma Rays From The Galactic Center,” Phys. Rev. D 84, 123005 (2011), arXiv:1110.0006 [astro-ph.HE] .
  * Abazajian and Kaplinghat [2012] Kevork N. Abazajian and Manoj Kaplinghat, “Detection of a Gamma-Ray Source in the Galactic Center Consistent with Extended Emission from Dark Matter Annihilation and Concentrated Astrophysical Emission,” Phys. Rev. D86, 083511 (2012), arXiv:1207.6047 [astro-ph.HE] .
  * Gordon and Macias [2013] Chris Gordon and Oscar Macias, “Dark Matter and Pulsar Model Constraints from Galactic Center Fermi-LAT Gamma Ray Observations,” Phys. Rev. D88, 083521 (2013), [Erratum: Phys. Rev.D89,no.4,049901(2014)], arXiv:1306.5725 [astro-ph.HE] .
  * Zhou  _et al._ [2015] Bei Zhou, Yun-Feng Liang, Xiaoyuan Huang, Xiang Li, Yi-Zhong Fan, Lei Feng, and Jin Chang, “GeV excess in the Milky Way: The role of diffuse galactic gamma-ray emission templates,” Phys. Rev. D 91, 123010 (2015), arXiv:1406.6948 [astro-ph.HE] .
  * Calore  _et al._ [2015a] Francesca Calore, Ilias Cholis, and Christoph Weniger, “Background Model Systematics for the Fermi GeV Excess,” JCAP 03, 038 (2015a), arXiv:1409.0042 [astro-ph.CO] .
  * Daylan  _et al._ [2016] Tansu Daylan, Douglas P. Finkbeiner, Dan Hooper, Tim Linden, Stephen K. N. Portillo, Nicholas L. Rodd, and Tracy R. Slatyer, “The characterization of the gamma-ray signal from the central Milky Way: A case for annihilating dark matter,” Phys. Dark Univ. 12, 1–23 (2016), arXiv:1402.6703 [astro-ph.HE] .
  * Ajello  _et al._ [2016] M. Ajello _et al._ (Fermi-LAT), “Fermi-LAT Observations of High-Energy $\gamma$-Ray Emission Toward the Galactic Center,” Astrophys. J. 819, 44 (2016), arXiv:1511.02938 [astro-ph.HE] .
  * Calore  _et al._ [2015b] Francesca Calore, Ilias Cholis, Christopher McCabe, and Christoph Weniger, “A Tale of Tails: Dark Matter Interpretations of the Fermi GeV Excess in Light of Background Model Systematics,” Phys. Rev. D91, 063003 (2015b), arXiv:1411.4647 [hep-ph] .
  * Agrawal  _et al._ [2015] Prateek Agrawal, Brian Batell, Patrick J. Fox, and Roni Harnik, “WIMPs at the Galactic Center,” JCAP 05, 011 (2015), arXiv:1411.2592 [hep-ph] .
  * Murgia [2020] Simona Murgia, “The Fermi–LAT Galactic Center Excess: Evidence of Annihilating Dark Matter?” Ann. Rev. Nucl. Part. Sci. 70, 455–483 (2020).
  * Di Mauro [2021] Mattia Di Mauro, “Characteristics of the Galactic Center excess measured with 11 years of $Fermi$-LAT data,” Phys. Rev. D 103, 063029 (2021), arXiv:2101.04694 [astro-ph.HE] .
  * Abazajian [2011] Kevork N. Abazajian, “The Consistency of Fermi-LAT Observations of the Galactic Center with a Millisecond Pulsar Population in the Central Stellar Cluster,” Journal of Cosmology and Astroparticle Physics 1103, 010 (2011), arXiv:1011.4275 [astro-ph.HE] .
  * Bartels  _et al._ [2018] Richard Bartels, Emma Storm, Christoph Weniger, and Francesca Calore, “The Fermi-LAT GeV excess as a tracer of stellar mass in the Galactic bulge,” Nature Astronomy 2, 819–828 (2018), arXiv:1711.04778 [astro-ph.HE] .
  * Macias  _et al._ [2018] Oscar Macias, Chris Gordon, Roland M. Crocker, Brendan Coleman, Dylan Paterson, Shunsaku Horiuchi, and Martin Pohl, “Galactic bulge preferred over dark matter for the Galactic centre gamma-ray excess,” Nature Astronomy 2, 387–392 (2018), arXiv:1611.06644 [astro-ph.HE] .
  * Macias  _et al._ [2019] Oscar Macias, Shunsaku Horiuchi, Manoj Kaplinghat, Chris Gordon, Roland M. Crocker, and David M. Nataf, “Strong Evidence that the Galactic Bulge is Shining in Gamma Rays,” JCAP 09, 042 (2019), arXiv:1901.03822 [astro-ph.HE] .
  * Ackermann  _et al._ [2017] M. Ackermann _et al._ (Fermi-LAT), “The Fermi Galactic Center GeV Excess and Implications for Dark Matter,” Astrophys. J. 840, 43 (2017), arXiv:1704.03910 [astro-ph.HE] .
  * Pohl  _et al._ [2022] Martin Pohl, Oscar Macias, Phaedra Coleman, and Chris Gordon, “Assessing the Impact of Hydrogen Absorption on the Characteristics of the Galactic Center Excess,” Astrophys. J. 929, 136 (2022), arXiv:2203.11626 [astro-ph.HE] .
  * Cholis  _et al._ [2022] Ilias Cholis, Yi-Ming Zhong, Samuel D. McDermott, and Joseph P. Surdutovich, “Return of the templates: Revisiting the Galactic Center excess with multimessenger observations,” Phys. Rev. D 105, 103023 (2022), arXiv:2112.09706 [astro-ph.HE] .
  * McDermott  _et al._ [2023] Samuel D. McDermott, Yi-Ming Zhong, and Ilias Cholis, “On the morphology of the gamma-ray galactic centre excess,” Mon. Not. R. Astron. Soc. 522, L21–L25 (2023), arXiv:2209.00006 [astro-ph.HE] .
  * Petrović  _et al._ [2015] Jovana Petrović, Pasquale D. Serpico, and Gabrijela Zaharijas, “Millisecond pulsars and the Galactic Center gamma-ray excess: the importance of luminosity function and secondary emission,” JCAP 02, 023 (2015), arXiv:1411.2980 [astro-ph.HE] .
  * Gautam  _et al._ [2022] Anuj Gautam, Roland M. Crocker, Lilia Ferrario, Ashley J. Ruiter, Harrison Ploeg, Chris Gordon, and Oscar Macias, “Millisecond pulsars from accretion-induced collapse as the origin of the Galactic Centre gamma-ray excess signal,” Nature Astron. 6, 703–707 (2022), arXiv:2106.00222 [astro-ph.HE] .
  * Dinsmore and Slatyer [2022] Jack T. Dinsmore and Tracy R. Slatyer, “Luminosity functions consistent with a pulsar-dominated Galactic Center excess,” JCAP 06, 025 (2022), arXiv:2112.09699 [astro-ph.HE] .
  * Ackermann  _et al._ [2016] M. Ackermann _et al._ (Fermi-LAT), “Resolving the Extragalactic $\gamma$-Ray Background above 50 GeV with the Fermi Large Area Telescope,” Phys. Rev. Lett. 116, 151105 (2016), arXiv:1511.00693 [astro-ph.CO] .
  * Di Mauro  _et al._ [2018] Mattia Di Mauro, Silvia Manconi, Hannes-S. Zechlin, Marco Ajello, Eric Charles, and Fiorenza Donato, “Deriving the contribution of blazars to the Fermi-LAT Extragalactic $\gamma$-ray background at $E>10$ GeV with efficiency corrections and photon statistics,” Astrophys. J. 856, 106 (2018), arXiv:1711.03111 [astro-ph.HE] .
  * Malyshev and Hogg [2011] Dmitry Malyshev and David W. Hogg, “Statistics of Gamma-Ray Point Sources below the Fermi Detection Limit,” Astrophys. J. 738, 181 (2011), arXiv:1104.0010 [astro-ph.CO] .
  * Lee  _et al._ [2016] Samuel K. Lee, Mariangela Lisanti, Benjamin R. Safdi, Tracy R. Slatyer, and Wei Xue, “Evidence for Unresolved Gamma-Ray Point Sources in the Inner Galaxy,” Phys. Rev. Lett. 116, 051103 (2016), arXiv:1506.05124 [astro-ph.HE] .
  * Mishra-Sharma  _et al._ [2017] Siddharth Mishra-Sharma, Nicholas L. Rodd, and Benjamin R. Safdi, “NPTFit: A code package for Non-Poissonian Template Fitting,” Astron. J. 153, 253 (2017), arXiv:1612.03173 [astro-ph.HE] .
  * Zechlin  _et al._ [2016a] Hannes-S. Zechlin, Alessandro Cuoco, Fiorenza Donato, Nicolao Fornengo, and Andrea Vittino, “Unveiling the Gamma-ray Source Count Distribution Below the Fermi Detection Limit with Photon Statistics,” Astrophys. J. Suppl. 225, 18 (2016a), arXiv:1512.07190 [astro-ph.HE] .
  * Zechlin  _et al._ [2016b] Hannes-S. Zechlin, Alessandro Cuoco, Fiorenza Donato, Nicolao Fornengo, and Marco Regis, “Statistical Measurement of the Gamma-ray Source-count Distribution as a Function of Energy,” Astrophys. J. Lett. 826, L31 (2016b), arXiv:1605.04256 [astro-ph.HE] .
  * Zechlin  _et al._ [2018] Hannes-S. Zechlin, Silvia Manconi, and Fiorenza Donato, “Constraining Galactic dark matter with gamma-ray pixel counts statistics,” Phys. Rev. D 98, 083022 (2018), arXiv:1710.01506 [astro-ph.HE] .
  * Lisanti  _et al._ [2016] Mariangela Lisanti, Siddharth Mishra-Sharma, Lina Necib, and Benjamin R. Safdi, “Deciphering Contributions to the Extragalactic Gamma-Ray Background from 2 GeV to 2 TeV,” Astrophys. J. 832, 117 (2016), arXiv:1606.04101 [astro-ph.HE] .
  * Manconi  _et al._ [2020] Silvia Manconi, Michael Korsmeier, Fiorenza Donato, Nicolao Fornengo, Marco Regis, and Hannes Zechlin, “Testing gamma-ray models of blazars in the extragalactic sky,” Phys. Rev. D 101, 103026 (2020), arXiv:1912.01622 [astro-ph.HE] .
  * Leane and Slatyer [2019] Rebecca K. Leane and Tracy R. Slatyer, “Revival of the Dark Matter Hypothesis for the Galactic Center Gamma-Ray Excess,” Phys. Rev. Lett. 123, 241101 (2019), arXiv:1904.08430 [astro-ph.HE] .
  * Leane and Slatyer [2020a] Rebecca K. Leane and Tracy R. Slatyer, “Spurious Point Source Signals in the Galactic Center Excess,” Phys. Rev. Lett. 125, 121105 (2020a), arXiv:2002.12370 [astro-ph.HE] .
  * Leane and Slatyer [2020b] Rebecca K. Leane and Tracy R. Slatyer, “The enigmatic Galactic Center excess: Spurious point sources and signal mismodeling,” Phys. Rev. D 102, 063019 (2020b), arXiv:2002.12371 [astro-ph.HE] .
  * Buschmann  _et al._ [2020] Malte Buschmann, Nicholas L. Rodd, Benjamin R. Safdi, Laura J. Chang, Siddharth Mishra-Sharma, Mariangela Lisanti, and Oscar Macias, “Foreground Mismodeling and the Point Source Explanation of the Fermi Galactic Center Excess,” Phys. Rev. D 102, 023023 (2020), arXiv:2002.12373 [astro-ph.HE] .
  * Calore  _et al._ [2021] F. Calore, F. Donato, and S. Manconi, “Dissecting the Inner Galaxy with $\gamma$-Ray Pixel Count Statistics,” Phys. Rev. Lett. 127, 161102 (2021), arXiv:2102.12497 [astro-ph.HE] .
  * Storm  _et al._ [2017] Emma Storm, Christoph Weniger, and Francesca Calore, “SkyFACT: High-dimensional modeling of gamma-ray emission with adaptive templates and penalized likelihoods,” JCAP 08, 022 (2017), arXiv:1705.04065 [astro-ph.HE] .
  * Mishra-Sharma and Cranmer [2022] Siddharth Mishra-Sharma and Kyle Cranmer, “Neural simulation-based inference approach for characterizing the Galactic Center $\gamma$-ray excess,” Phys. Rev. D 105, 063017 (2022), arXiv:2110.06931 [astro-ph.HE] .
  * Shmakov  _et al._ [2023] Alexander Shmakov, Mohammadamin Tavakoli, Pierre Baldi, Christopher M. Karwin, Alex Broughton, and Simona Murgia, “Deep learning models of the discrete component of the Galactic interstellar $\gamma$-ray emission,” Phys. Rev. D 107, 063018 (2023), arXiv:2206.02819 [astro-ph.HE] .
  * Linden  _et al._ [2016] Tim Linden, Nicholas L. Rodd, Benjamin R. Safdi, and Tracy R. Slatyer, “High-energy tail of the Galactic Center gamma-ray excess,” Phys. Rev. D 94, 103013 (2016), arXiv:1604.01026 [astro-ph.HE] .
  * Ackermann  _et al._ [2014] M. Ackermann _et al._ (Fermi-LAT), “The Spectrum and Morphology of the $Fermi$ Bubbles,” Astrophys. J. 793, 64 (2014), arXiv:1407.7905 [astro-ph.HE] .
  * Herold and Malyshev [2019] Laura Herold and Dmitry Malyshev, “Hard and bright gamma-ray emission at the base of the Fermi bubbles,” Astron. Astrophys. 625, A110 (2019), arXiv:1904.01454 [astro-ph.HE] .
  * Wolleben [2007] M. Wolleben, “A New Model For The Loop-I (The North Polar Spur) Region,” Astrophys. J. 664, 349–356 (2007), arXiv:0704.0276 [astro-ph] .
  * Abdollahi  _et al._ [2022] Soheila Abdollahi _et al._ (Fermi-LAT), “Incremental Fermi Large Area Telescope Fourth Source Catalog,” Astrophys. J. Supp. 260, 53 (2022), arXiv:2201.11184 [astro-ph.HE] .
  * Ballet  _et al._ [2023] J. Ballet, P. Bruel, T. H. Burnett, and B. Lott (Fermi-LAT), “Fermi Large Area Telescope Fourth Source Catalog Data Release 4 (4FGL-DR4),” arXiv preprint (2023), arXiv:2307.12546 [astro-ph.HE] .
  * Ackermann  _et al._ [2015] M. Ackermann _et al._ , “The Spectrum of Isotropic Diffuse Gamma-Ray Emission between 100 MeV and 820 GeV,” Astrophys. J. 799, 86 (2015), arXiv:1410.3696 [astro-ph.HE] .
  * Ackermann  _et al._ [2018] M. Ackermann _et al._ (Fermi-LAT), “Unresolved Gamma-Ray Sky through its Angular Power Spectrum,” Phys. Rev. Lett. 121, 241101 (2018), arXiv:1812.02079 [astro-ph.HE] .
  * Maccione  _et al._ [2011] Luca Maccione, Carmelo Evoli, Daniele Gaggero, and Dario Grasso, “DRAGON: Galactic Cosmic Ray Diffusion Code,” (2011), ascl:1106.011 .
  * Porter  _et al._ [2022] Troy A. Porter, Gudlaugur Johannesson, and Igor V. Moskalenko, “The GALPROP Cosmic-ray Propagation and Nonthermal Emissions Framework: Release v57,” Astrophys. J. Supp. 262, 30 (2022), arXiv:2112.12745 [astro-ph.HE] .
  * Ackermann and others [2012] M. Ackermann and others, “Fermi-LAT Observations of the Diffuse $\gamma$-Ray Emission: Implications for Cosmic Rays and the Interstellar Medium,” Astrophys. J. 750, 3 (2012), arXiv:1202.4039 [astro-ph.HE] .
  * McCann [2015] A. McCann, “A stacked analysis of 115 pulsars observed by the Fermi LAT,” Astrophys. J. 804, 86 (2015), arXiv:1412.2422 [astro-ph.HE] .
  * Cao  _et al._ [2013] Liang Cao, Shude Mao, David Nataf, Nicholas J. Rattenbury, and Andrew Gould, “A new photometric model of the Galactic bar using red clump giants,” Mon. Not. Roy. Astron. Soc. 434, 595–605 (2013), arXiv:1303.6430 [astro-ph.GA] .
  * Coleman  _et al._ [2020] B. Coleman, D. Paterson, C. Gordon, O. Macias, and H. Ploeg, “Maximum entropy estimation of the Galactic bulge morphology via the VVV Red Clump,” Mon. Not. Roy. Astron. Soc. 495, 3350–3372 (2020), arXiv:1911.04714 [astro-ph.GA] .
  * Zhong and Cholis [2024] Yi-Ming Zhong and Ilias Cholis, “Robustness of the Galactic Center Excess Morphology Against Masking,” arXiv e-prints , arXiv:2401.02481 (2024), arXiv:2401.02481 [astro-ph.HE] .
  * Song  _et al._ [2024] Deheng Song _et al._ , “Robust inference of the Galactic center gamma-ray excess spatial properties,” In preparation (2024).
  * Acero  _et al._ [2016] F. Acero _et al._ , “Development of the Model of Galactic Interstellar Emission for Standard Point-source Analysis of Fermi Large Area Telescope Data,” Astrophys. J. Suppl. 223, 26 (2016), arXiv:1602.07246 [astro-ph.HE] .
  * Feroz  _et al._ [2009] F. Feroz, M. P. Hobson, and M. Bridges, “MULTINEST: an efficient and robust Bayesian inference tool for cosmology and particle physics,” Mon. Not. Roy. Astron. Soc. 398, 1601–1614 (2009), arXiv:0809.3437 .
  * Rolke  _et al._ [2005] W. A. Rolke, A. M. López, and J. Conrad, “Limits and confidence intervals in the presence of nuisance parameters,” Nuclear Instruments and Methods in Physics Research A 551, 493–503 (2005), physics/0403059 .
  * [63] Publicly available at https://heasarc.gsfc.nasa.gov/
FTP/fermi/data/lat/weekly/photon/. More details are found at
https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Data/LAT_DP.html.

  * Gorski  _et al._ [2005] K. M. Gorski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann, “HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere,” Astrophys. J. 622, 759–771 (2005), astro-ph/0409513 .
  * Abdollahi  _et al._ [2020] S. Abdollahi _et al._ (Fermi-LAT), “$Fermi$ Large Area Telescope Fourth Source Catalog,” Astrophys. J. Suppl. 247, 33 (2020), arXiv:1902.10045 [astro-ph.HE] .
  * Di Mauro and Winkler [2021] Mattia Di Mauro and Martin Wolfgang Winkler, “Multimessenger constraints on the dark matter interpretation of the Fermi-LAT Galactic center excess,” Phys. Rev. D 103, 123005 (2021), arXiv:2101.11027 [astro-ph.HE] .
  * Calore  _et al._ [2022] Francesca Calore, Marco Cirelli, Laurent Derome, Yoann Genolini, David Maurin, Pierre Salati, and Pasquale Dario Serpico, “AMS-02 antiprotons and dark matter: Trimmed hints and robust bounds,” SciPost Phys. 12, 163 (2022), arXiv:2202.03076 [hep-ph] .
  * Balan  _et al._ [2023] Sowmiya Balan, Felix Kahlhoefer, Michael Korsmeier, Silvia Manconi, and Kathrin Nippel, “Fast and accurate AMS-02 antiproton likelihoods for global dark matter fits,” JCAP 08, 052 (2023), arXiv:2303.07362 [hep-ph] .
  * Foster  _et al._ [2023] Joshua W. Foster, Yujin Park, Benjamin R. Safdi, Yotam Soreq, and Weishuang Linda Xu, “Search for dark matter lines at the Galactic Center with 14 years of Fermi data,” Phys. Rev. D 107, 103047 (2023), arXiv:2212.07435 [hep-ph] .

##  Appendix A skyFACT fits: systematics and fit ranges

In this section, we test our skyFACT fit results against a number of
systematics.

Figure A7:  Systematic tests on GCE spectrum. Left: GCE best-fit spectra
(BB+NB components) reconstructed with the skyFACT fit in the range 0.3 – 300
GeV, adopting different input spectra: a power law with index -2.5 (deep pink
solid, baseline choice), a power law with exponential cutoff with parameters
$\alpha=-1.46$ and $E_{0}=3.6$ GeV (purple, typical MSP-like “prompt”
spectrum), or $\alpha=-2$ and $E_{0}=50$ GeV (violet). Dotted lines correspond
to the input spectra. Right: GCE best-fit spectra (BB+NB components)
reconstructed with the skyFACT fit in the range 10 – 300 GeV (orange solid,
baseline), testing uncertainties related to the input spectrum of the
innermost gas ring (green, gas I); spectral freedom on the Fermi bubbles
(royal blue, 25% and navy, 50%); spatial morphology of the Fermi bubbles
(teal, structured morphology). We also show the GCE Bulge spectrum when the
new BB determination in Ref. [57] is adopted (red line, BB Coleman).

First, we consider the fit in the 0.3 – 300 GeV. In figure A7 (left), we show
the best-fit spectra of the GCE Bulge component resulting from the skyFACT
when assuming different GCE input spectra: (i) a power law with index
$\alpha=-2.5$ (baseline of this work); a power law with exponential cutoff
with parameters (ii) $\alpha=-1.46$ and $E_{0}=3.6$ GeV (as assumed in
previous analyses as a proxy for the cumulative MSP spectrum, and also adopted
in PaperI) or (iii) $\alpha=-2$ and $E_{0}=50$ GeV. The latter input spectrum
is meant to provide significant input flux at high energies so to mimic the
behavior of the pure power law. From figure A7 (left), we can conclude that
the GCE spectral behavior is consistent regardless of the assumed input
spectrum up to $E\sim 10$ GeV. At higher energies, the (ii) input spectrum is
basically zero and cannot be properly modulated despite the full spectral
freedom. Instead, input spectra with non-zero high-energy values, i.e. like
(i) and (iii), can be modulated also at higher energies. The best-fit GCE
spectra we find for (i) and (iii) are compatible above 10 GeV, showing that,
as long as the input spectrum is non-zero in the energy range of interest
skyFACT successfully modulates it providing robust results against the choice
of the input spectral parametrization. In all cases, there is strong evidence
for the GCE, and the Bulge is needed on top of the NFW126 component.

Secondly, we consider the high-energy interval 10 – 300 GeV, and perform
skyFACT fits therein. As mentioned in section IV, the evidence of the GCE at
higher energies is still strong, and a Bulge component is still required on
top of NFW126. We here perform a series of runs to test different systematics.
First, we noticed that the innermost gas ring (gas I) resulted to be strongly
suppressed with respect to the best-fit in the 0.3 – 300 GeV skyFACT run. We
can explain that because of the fact that the higher energy data possess less
constraining power on this component given the reduced statistics and the fact
that the gas emission is the stronger at low energies, where the gamma-ray
production from neutral pion decay is peaked. We setup a run in which the gas
I component is therefore initialized to same best-fit value of the 0.3 – 300
GeV skyFACT run. The spectral behavior of the GCE component is stable against
this variation, and the evidence for GCE and preference for Bulge at high
energies are unaffected. We display the high-energy GCE Bulge best-fit
spectrum in figure A7 (right).

The modeling of the Fermi bubbles [46] was found to represent an important
source of systematic uncertainty in the evidence of a high-energy tail of the
GCE. In particular, Ref. [19] found that whenever a low-latitude Fermi bubbles
template was added to the fit, the GCE high-energy tail ($E>10$ GeV)
disappeared, see their figure 9, right panel. We would like to stress here
that the Fermi bubbles are an effective gamma-ray model component, defined in
skyFACT through their spectrum. We constrain the spectrum to be compatible
with the best-fit spectral behavior at high-latitude found in Ref. [45], which
is approximately $\propto E^{-1.9}$, and allow it to vary only mildly. On the
other hand, the spatial morphology is completely reconstructed by skyFACT
through modulation of the spatial coefficients which are fully free to vary,
see details in Ref. [41, 40]. Therefore, by construction, for all skyFACT runs
the Fermi bubbles do possess a low-latitude component.

Nevertheless, given the possible strong degenerancies between Fermi bubbles
and GCE, we test some possible systematics related to modeling of the Fermi
bubbles. In our framework, given the degenerancies at play, we cannot obtain
convergent fit if both Fermi bubbles and GCE are left free spatially and
spectrally, so that we need to constrain partly their model in order to get
meaningful results. We therefore test the following setups: (a) We
progressively increase the spectral freedom of the Fermi bubbles from the 1%
variation in the baseline setup to 25% and 50%; (b) Instead of leaving the
morphology of the Fermi bubbles free to readjust during the fit, we fix their
spatial part to the best-fit structured template from Ref. [18] and allow full
spectral freedom. In figure A7 (right), we show that the high-energy GCE
spectrum is stable against variations of Fermi bubbles spectrum and
morphology. The evidence for GCE is always present for these model variations,
albeit slightly reduced (to 5.3 $\sigma$) when using a fixed, structured Fermi
bubbles morphology, which however provides overall a worse fit than our
baseline setup. We therefore find strong evidence for the GCE at high energy
in the presence of the Fermi bubbles at low latitudes.

Finally, we also test a more recent model of the stellar bulge [57], which has
been found to better fit the gamma-ray data in the full energy range [58, 59].
The best-fit Bulge spectrum at high energy is fully consistent when different
bulge templates are used.

Figure A8:  Systematic tests for the IG $dN/dS$ reconstruction using different
skyFACT diffuse emission models.  The source count distribution of the IG
(masking $b<|1|^{\circ}$) as obtained from the 1pPDF analysis is measured
varying various ingredients of the skyFACT Galactic diffuse emission
derivation. The ’Full E range’ (magenta dashed) corresponds to the case in
which all the diffuse templates are taken from the skyFACT fit in the 0.3–300
GeV interval. See the right panel of figure A7 and text for more details.

##  Appendix B 1pPDF fits: systematics and fit ranges

To demonstrate that the most crucial systematics on the skyFACT fit explored
in the previous section have no impact on the 1pPDF results in the IG, we show
in figure A8 the $dN/dS$ reconstructed for various cases. The results of the
Baseline fit obtained in the 10–300 GeV are found to be consistent with the
ones obtained using the skyFACT diffuse models fitted to the full energy range
(magenta dashed), as well as the ones using the structured FB template (blue
dotted), no matter the latitude cut. Additionally, the $\ln(\mathcal{Z})$ for
the fit using the skyFACT diffuse models obtained in the full 0.3–300 GeV
range are respectively -19118.29 and -18119.29 for the latitude cut at
$0.5^{\circ}$ and $1^{\circ}$, indicating that the skyFACT diffuse models
optimized for the high energy fit better describe the gamma-ray sky fitted by
the 1pPDF.

Table AIII: 1pPDF analysis of the IG: parameter priors. See text for details.  Method | Parameter | Prior | Range
---|---|---|---
| $A_{\rm gal}$ | flat | [0.5,1.5]
| $F_{\rm iso}$ | log-flat | [$4\cdot 10^{-9}$, $5\cdot 10^{-8}$]
| $A_{\rm B/NFW126}$ | flat | [$0.5$, $1.5$]
MBPL | $A_{S}$ | log-flat | [$10^{9}$, $5\cdot 10^{13}$]
| $S_{\rm b1}$ | log-flat | [$3\cdot 10^{-11}$, $10^{-9}$]
| $S_{\rm b2}$ | log-flat | [$1\cdot 10^{-14}$, $3\cdot 10^{-11}$]
| $n_{\rm 1}$ | flat | [-1, 3]
| $n_{\rm 2}$ | flat | [1.5, 2.5]
| $n_{\rm 3}$ | flat | [-2 , 2 ]
Hybrid | $A_{S}$ | log-flat | [$10^{7}$, $10^{10}$]
| $S_{\rm b1}$ | log-flat | [$3\cdot 10^{-11}$, $10^{-9}$]
| $S_{\rm b2}$ | log-flat | [$10^{-12}$, $3\cdot 10^{-11}$]
| $n_{\rm 1}$ | flat | [2.05 ,7]
| $n_{\rm 2}$ | flat | [1.7,2.3]
| $n_{\rm 3}$ | flat | [-2,2]
| $A_{\rm nd1}$ | log-flat | [$10^{8}$, $10^{15}$]
| $S_{\rm nd1}$ | fixed | $10^{-13}$
| $n_{\rm f}$ | fixed | -10

Figure A9:  Systematic tests on $dN/dS$ reconstruction. Source count
distribution of the extragalactic ROI (upper panels), of the outer Galaxy
(center panels) and of the IG (masking $b<|1|^{\circ}$, and using the SF-BB
setup) as obtained from the 1pPDF analysis varying the 1pPDF fit method (left)
and the Fermi–LAT event selection (right). Consistent results are obtained in
all ROIs for all explored configurations.

The panels in figure A9 illustrate further systematics test performed on the
reconstructed $dN/dS$ within the 1pPDF fits. The upper, middle and lower rows
depict the results for the EG, OG and IG ROI, respectively. In each row, the
left panel shows that the 1pPDF results presented in the main text are robust
with respect of the modeling of faint point sources in the unresolved regime.
Indeed, modeling the $dN/dS$ within the Hybrid method [31], thus adding nodes
at fluxes $<10^{-12}$ ph cm-2 s-1 leaves the $dN/dS$ results compatible with
the MBPL fit down to the methods’s sensitivity. The right panels of figure A9
show instead how the $dN/dS$ results are robust against the choice of the
Fermi–LAT event selection. The measured $dN/dS$ is fully consistent within
$1\sigma$ uncertainties among CLEAN and UCV events, and also restricting to
the best quartile of PSF events (PSF3). The only observable effect is an
increased $1\sigma$ uncertainty band for event types including less photons,
specifically in small ROIs such as the OG and the IG.

We summarize in Table AIII the parameter priors for the 1pPDF analysis of the
IG. The table contains three blocks: in the first, parameters in common with
all the fit setups are given, such as the normalization of the diffuse
emission templates and of the isotropic emission. The second block refers to
the parameters for the $dN/dS$ fit when using the MBPL approach, specifically
the normalization of the $dN/dS$ the break positions and indexes, see Eq. (1).
The last block refers to the Hybrid approach, where the MBPL is extended with
a node. The normalizations $A_{S},\,A_{\rm nd1}$ are in units of s cm2sr-1,
and the break positions $S_{\rm bn,snd1}$ are in units of ph cm-2 s-1 .
Finally, $F_{\rm iso}$ is given in units of cm2s-1sr-1. All other parameters
are dimensionless.