# Mind the gap: The discrepancy between simulation and reality drives
interpretations of the Galactic Center Excess

Sascha Caron  , Christopher Eckner  , Luc Hendriks  , Guðlaugur Jóhannesson  ,
Roberto Ruiz de Austri  , Gabrijela Zaharijas

###### Abstract

The Galactic Center Excess (GCE) in GeV gamma rays has been debated for over a
decade, with the possibility that it might be due to dark matter annihilation
or undetected point sources such as millisecond pulsars (MSPs). This study
investigates how the gamma-ray emission model ($\gamma$EM) used in Galactic
center analyses affects the interpretation of the GCE’s nature. To address
this issue, we construct an ultra-fast and powerful inference pipeline based
on convolutional Deep Ensemble Networks. We explore the two main competing
hypotheses for the GCE using a set of $\gamma$EMs with increasing parametric
freedom. We calculate the fractional contribution ($f_{\mathrm{src}}$) of a
dim population of MSPs to the total luminosity of the GCE and analyze its
dependence on the complexity of the $\gamma$EM. For the simplest $\gamma$EM,
we obtain $f_{\mathrm{src}}=0.10\pm 0.07$, while the most complex model yields
$f_{\mathrm{src}}=0.79\pm 0.24.$ In conclusion, we find that the statement
about the nature of the GCE (dark matter or not) strongly depends on the
assumed $\gamma$EM. The quoted results for $f_{\mathrm{src}}$ do not account
for the additional uncertainty arising from the fact that the observed gamma-
ray sky is out-of-distribution concerning the investigated $\gamma$EM
iterations. We quantify the reality gap between our $\gamma$EMs using deep-
learning-based One-Class Deep Support Vector Data Description networks,
revealing that all employed $\gamma$EMs have gaps to reality. Our study casts
doubt on the validity of previous conclusions regarding the GCE and dark
matter, and underscores the urgent need to account for the reality gap and
consider previously overlooked “out of domain” uncertainties in future
interpretations.

##  1 Introduction

The all-sky gamma-ray emission in the GeV-TeV range contains a treasure trove
of information on the Galactic interstellar medium, the Galactic cosmic ray
(CR) population and CR accelerators. It primarily originates in the
interstellar emission (IE), produced by interactions between CRs and the
interstellar medium, in addition to the contribution of a large number of
individual sources. This emission was mapped in detail in the GeV energy range
by the Large Area Telescope (LAT) onboard the Fermi Gamma-Ray Space Telescope
[1], operating over the past decade, and at TeV energies, by the High-Altitude
Water Cherenkov Gamma-Ray Observatory (HAWC) [2] and High Energy Stereoscopic
System (H.E.S.S.) [3, 4, 5] telescopes. Some of the major issues, limiting the
amount of information that was so far extracted from the data, are the
relatively poor angular resolution of the detectors, limited statistics of the
observations, as well as degeneracy between various emission components that
is challenging to break with traditional analysis methods.

A prime example of a sky region that is particularly prone to these challenges
is the $\mathcal{O}$($10^{\circ}$) region surrounding the Galactic center
(GC). The Galaxy is mostly transparent to photons in this energy range except
at the highest energies [6, 7], and the emission is therefore integrated over
the entire diameter of the Galaxy, likely making it the most complex region of
the sky. This includes the GC which houses the supermassive black hole SgrA∗
[8] and has the highest star formation rate (SFR) in the Galaxy [9], both of
which are expected to impact the gamma-ray emission. In particular, the
population of CRs could differ from that in the rest of the Galaxy (due to,
e.g., freshly injected CRs [10]); their diffusion could be position-dependent
[11] and strong Galactic winds could be present in this region [12]. Due to
the high SFR a large number of individually unresolved sources is also
expected to be present in the region, with their cumulative emission
overlapping and possibly mimicking other emission components. The base of the
Fermi bubbles (FB), a large structure centered at the GC [13, 14], is also
expected to contribute to the measured emission, though the bubbles are
challenging to detect close to the Galactic plane. In addition, the GC is
expected to be the brightest place of a signal of dark matter (DM) particle
pair-annihilation, at the verge of the sensitivity of current telescopes for
WIMP DM models.

###  1.1 Understanding the GeV gamma-ray emission from the GC – a brief review

Excess emission over standard astrophysical backgrounds was detected in this
region and dubbed the Galactic center excess (GCE). First discovered in 2009
[15], this excess peaks at $1-3$ GeV with a spatial extension of $\sim
10^{\circ}$ (see also [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28]),
is robust against variations in the IE models (IEMs) [29], and is nearly
symmetric around the GC [30, 14, 31, 32, 33]. Possible explanations for the
GCE include a dim population of point sources, in particular, millisecond
pulsars (MSPs) [18, 20, 23, 34, 35, 36, 37, 38, 39, 40, 41], enhanced gamma-
ray emission due to cosmic-ray interactions in dense molecular clouds [42, 43]
or Galactic winds/past transient event(s) [10, 44]. Intriguingly, the GCE is
as well as consistent with the expectations for thermal dark matter pair-
annihilating into, e.g., quark final states and a mass around $40-60$ GeV
while the emission’s spatial profile is consistent with an adiabatically
contracted Navarro-Frenk-White (NFW) profile [45, 20, 23, 24, 29]. While there
is consensus in the community about the presence of the GCE, its properties
are less firmly established. Early works on the subject (e.g., [22, 30])
already pointed out that large uncertainties in the diffuse background
modeling impact the reconstruction of the GCE’s spatial and spectral
properties. As illustrated by the summary figure 1 in [46] that compares GCE
spectra obtained in a large number of works, it is not clear whether the GCE’s
spectrum extends to energies beyond $\sim 10$ GeV. There is a recent debate
about the preferred spatial morphology of the GCE based on the results of
template fitting. The authors of [31, 32, 33] reconstruct the GCE as entirely
consistent with spherical symmetry while – at least in the latter two analyses
– accounting for the uncertainty of the IE by varying the assumptions on
cosmic-ray propagation parameters. In contrast, [47] demonstrates that re-
simulating the gas distribution in the GC region leads to a strong preference
for a GCE following the asymmetric stellar density of the Galactic bulge,
without the need for DM as claimed in [31]. The preference for the stellar
bulge morphology has been pointed out many times in the literature by
independent groups and employing diverse inference methods [48, 49, 50, 51,
52].

The preference for a stellar bulge component can be interpreted as
corroborative evidence for the point-source nature of the GCE. The stellar
mass of the Galactic bulge of about $10^{10}\;M_{\odot}$ [53, 54] is dominated
by old stellar populations (> 5 Gyr). The innermost $\sim 200$ pc of the MW
feature an additional very-high-density region of stars, the nuclear stellar
cluster [55, 54]. Due to their age or stellar density, these components of the
MW are ideal places for the formation and presence of MSPs [56, 57, 58, 37].
However, the question of the point-source nature of the GCE has commonly been
treated fully decoupled from the question of the preferred spatial profile of
the GCE. Early works that leveraged the discriminating power of small-scale
photon clustering in the residual gamma-ray emission of the excess employed
wavelet transforms [34, 59] or, so-called, non-Poissonian template fitting
(NPTF) [35]. Both of these seminal papers found evidence for the point-source
nature of the GCE. Yet, the results of [34] were later questioned by an
independent study utilizing a dataset with a larger exposure and an updated
gamma-ray source catalog [60]. The initial results of the NPTF method were
followed by further analyses using the same approach, which cast doubt on the
robustness of the technique. In [61] an independent group published its
results that the NPTF approach always reconstructs an injected DM-like signal
– with an overall luminosity below a certain threshold – as point-source-like
emission. These claims were rebutted in [62, 63] showing explicitly that
diffuse background mis-modeling can induce the observed behavior of the NPTF
method. Employing instead an improved diffuse background model leads again to
the conclusion that the GCE is, at least partially, made out of point sources.
The response to these rebuttals was given in form of two publications [64, 65]
showing that an unmodeled asymmetry of the GCE emission can induce spurious
evidence for the point-source nature of the excess.

These recent developments in exploring the limitations of frequently employed
inference techniques, in particular template fitting and NPTF, revealed that
background mis-modeling plays a crucial role in examining the nature of the
GCE. To mitigate these uncertainties, several schemes have been proposed and
applied in the literature: Introducing nuisance parameters on small- and large
angular scales via spherical harmonic marginalization [63], varying diffuse
templates via Gaussian processes [66] or adaptive template-fitting with
SkyFACT [67, 49] that combines template fitting with image reconstructions
techniques. The latter method is utilizing a large number of nuisance
parameters $\mathcal{O}(>10^{5})$ and associated regularization terms to
prevent over-fitting and to optimize the compiled background model in a data-
driven way. The authors of [52] employed a SkyFACT-optimized IEM together with
the 1-point probability distribution function (1pPDF) technique [68] to
discern the contribution of sub-threshold sources to the gamma-ray emission of
the GC. The analysis showed that the 1pPDF’s sensitivity to dim sources
largely depends on the chosen IEM. However, the SkyFACT-optimized model is
consistent with 4FGL catalog sources and it exhibits the highest sensitivity
to point sources below the Fermi-LAT detection threshold.

Another means to extract the small-scale photon clustering in images of the GC
is via the help of deep learning and Convolutional Neural Networks (CNNs). A
range of new analysis methods was proposed to address this question [69, 70,
71, 72], with mixed results. Our previous work [69] pioneered the application
of CNN to this region of the sky and provided a proof of principle that such
methods could indeed successfully discern between the genuinely diffuse and
point source origin of the emission (albeit in a rather limited setting).
Building on that effort, we take a significantly more ambitious approach in
this work, improving both the network architecture and the modeling of
astrophysical emission, as detailed below.

###  1.2 The issue of discrepancies between model and reality

As discussed earlier, one of the biggest general issues when dealing with the
GCE is that its properties strongly depend on our knowledge (or lack thereof)
of the bright fore-/background emission originating from the rest of the
Galaxy. While we have a good broad description of the IE, detailed knowledge
of the distribution of CRs and their targets (gas and radiation) near the GC
is imperfect. This could be dealt with by either improving our models via
observational data or advanced inference techniques (see, e.g., [73, 74, 75,
76] or by adding flexibility to the IEM to parameterize away our imperfect
knowledge. Yet, there is always a balance between allowing for too much
freedom, which may explain any data, and not allowing for enough, which
results in spurious deficits and excesses that require additional, and
possibly erroneous, model components. This problem is inherent to both the
more traditional maximum likelihood and other methods (see e.g. [64, 63]) as
well as in our application of CNNs which rely on a gamma-ray emission model
($\gamma$EM) and simulations of reality from it.

The difficulty of harmonizing our incomplete knowledge of the full complexity
of observational data with the real world is often called the “reality gap”
and, as a result, a DM or other components might be needed to describe the
data only due to a wrong or missing component in the $\gamma$EM.

Figure 1: Visualisation of domain adaptation. The background model to generate
the GC images is expanded in dimensionality to make more complex (and
hopefully more accurate models). Model 1 is embedded into Model 2, which is
embedded in Model 3. Additionally, we define two scenarios with separate
parameter spaces but following the same domain adaptation procedures. This
results in six different background models. The output space of all background
models is always the same: a GC image with 5 energy bins. In each respective
case, this output is the input of the neural networks that predict the input
parameters of the models. For every background model, we trained a neural
network. Hence, there is a “Model 1, scenario A network” which is trained on
images generated from background Model 1, Scenario A, and also for Model 2 and
3 and Scenario B. Additionally, it is possible to apply a network trained on
simulations from Model 1 on simulations from Model 2 for example. The arrows
represent all experiments we have done: the six straight arrows represent
networks that are evaluated and trained on the same background model. The
curved dashed arrows represent experiments where the neural network of a
trained model is used to predict parameters from simulations of a more complex
model (i.e. a network trained on Model 1 which is applied on Model 2, and the
same for Model 2 and 3). The elbow-connected dashed lines represent
experiments where a neural network trained on Scenario A is applied to
Scenario B and the other way around.

In this work, we aim at answering the question of the nature of the GCE by
leveraging the small-scale photon clustering in gamma-ray data of the GC. At
the same time, we investigate how much the background model complexity
influences the results of our study. To this end, we compile three different
$\gamma$EM that are supposed to simulate the GC region. We dub these instances
Models 1, 2 and 3, which have 6, 17 and 24 free parameters respectively. Note
that Model 1A – the nomenclature becomes clear in the next paragraph – is
chosen to be similar to the $\gamma$EM adopted in Ref. [69], allowing us to
make a direct comparison with our previous work.

The main idea that we follow is to start with $\gamma$EMs exhibiting a minimal
number of parameters and then gradually increase the complexity of the model
through human intervention motivated by inspecting the residuals of the
resulting best fit, thus also testing the robustness of our determination of
the GCE properties when varying the background model. In machine learning (ML)
language such approaches are known as a change of the “domain” (or “domain
shift”). Methods known as “domain adaptation” aim to minimize the difference
between two domains, so that models trained in a source domain can generalize
(and make correct predictions) to a target domain. A typical example is that
an ML inference model for road signs trained from US images (e.g. a particular
GCE simulation) may not be applicable to European domains (the real data) and
may yield incorrect results that uncertainties may not entirely reflect. Here
both data sets (simulated data from a compiled $\gamma$EM and real data) share
the same feature space: diffuse emission following the distribution of gas and
the structure of radiation fields in the Galaxy, bright and dim gamma-ray
point sources, the Fermi Bubbles and the excess emission in the GC. In
particular, we try to test in what way such a gap between simulations and
reality is causing different interpretations of the GCE. Testing domain
changes uses here (training) data in Simulation Model 1 to solve the task of
interpreting the GCE in a different domain (e.g. Simulation Model 2). In our
case Model 1 is nested in Model 2 etc., see Fig. 1. Therefore, to adapt to the
different domains, we try to learn from the more general $\gamma$EM a
generalized distribution that is also applicable to the sub-model (and
hopefully to the real data).

In addition, we consider two independent (non-nested) $\gamma$EM-sets, dubbed
Scenario A and B (see Fig. 1 where scenario B is presented as ‘parallel’ to
Scenario A ). The two scenarios are motivated by uncertainties in the modeling
of the dim point source components and are therefore not nested. Such an
approach, of slowly adding parameters until a $\gamma$EM is found that leads
to featureless residuals, is complementary to techniques like the approach of
SkyFACT.

###  1.3 Contribution to the field

Our work presented here improves, extends and complements our previous proof-
of-principle study [69]. In addition, we advance the current state-of-the-art
of GCE analyses in the following aspects:

  * •

We create a series of increasingly complex template-based $\gamma$EMs for the
GC region, adjusting the flexibility of the IEM parameters while modeling the
GCE with a combination of DM and MSP templates. We generate Monte Carlo MSP
populations using observation-based priors [77].

  * •

We develop a high-speed inference pipeline using Deep Ensemble Networks (DENs)
that reconstructs best-fitting model parameters and associated uncertainties
for Monte Carlo data or real Fermi-LAT gamma-ray sky observations. We
demonstrate that the network can effectively separate gamma-ray data into
emission components, accurately determining the spatial shape and total flux
of the GCE.

  * •

We discover that the fractional contribution ($f_{\mathrm{src}}$) of a dim MSP
population to the GCE’s total luminosity depends on the complexity of the
$\gamma$EM used, indicating a reality gap between our models and real data. We
employ One-Class Deep Support Vector Data Description (Deep SVDD) networks to
quantify this reality gap, providing a powerful way to determine if the Fermi-
LAT gamma-ray sky is out-of-distribution for a specific $\gamma$EM. We
conclude that a reality gap persists across all $\gamma$EMs explored, leading
to an additional out-of-distribution uncertainty on $f_{\mathrm{src}}$
predictions not captured by our networks. This uncertainty is likely to be
present in most previous GCE studies.

###  1.4 Outline

This paper is organized as follows: In Sec. 2 we elaborate on the LAT data
selection and we describe the building blocks of our $\gamma$EMs, detailing
how the simulations are performed and how the first model iteration is
constructed. In Sec. 3 we describe the inference techniques utilized in this
work, i.e. DENs and a likelihood-based approach. Sec. 4 states our results
from the first step of the simulation with the simplest $\gamma$EM that we
introduced in Sec. 2. The following Sec. 5 is dedicated to a discussion of the
approach leading to an extension of our initial $\gamma$EM. We describe the
construction of the more complex $\gamma$EM iterations Model 2 and 3 and we
present the results obtained with respect to each of them. In Sec. 6 we assess
the robustness of the results presented in the previous sections. We initiate
Sec. 6.1 by describing how the likelihood-based analysis compares to the ML
approach. Subsequently, in Sec. 6.2, we engage in a more detailed exploration
of the reality gap. The paper is concluded in Sec. 7. In addition, we give
details about some of the analysis aspects in appendices: In Appendix A we
summarize the predictions of all six neural networks regarding the
astrophysical emission components. Appendix B offers a more detailed
discussion of the physical consequences of MSP scenario B in terms of sources
below and above the Fermi-LAT detection threshold. We provide in Appendix C
the full list of spatial and spectral residuals according to the DNN results.
Lastly, Appendix D supplements Sec. 6.2 by stating the results obtained for
those $\gamma$EMs and scenarios not discussed in the main text.

##  2 Data and gamma-ray emission models

###  2.1 Fermi data selection

We use Pass 8 (R3) data [78, 79] taken from the time interval from the 4th of
August 2008 to the 2nd of April 2018 ($\sim$10 years). To reduce the residual
CR background and increase the resolution, we select events in the
ULTRACLEANVETO event class of the FRONT type. We also apply the standard
zenith angle cut ($<100^{\circ}$) as well as the time cut filters
DATA$\\_$QUAL==1 $\&\&$ LAT$\\_$CONFIG==1. All _Fermi_ -LAT data selection,
preparation, simulation and manipulation is conducted via the _Fermi Science
Tools_ 111https://github.com/fermi-lat/Fermitools-conda (version 2.0.8).

Our Region Of Interest (ROI) is a $30^{\circ}\times 30^{\circ}$ square
centered on the Galactic center and we bin it using square spatial bins of
size $0.25^{\circ}\times 0.25^{\circ}$ using the plate carrée projection222The
size of the ROI we chose is constrained by the requirement of having a good
angular resolution while keeping the size of training images feasible for our
analysis (see Section 3). From a physics perspective, a larger region could
also potentially over-constrain some of the components, if their spatial
morphology is incorrectly modeled closer to the GC.. We consider the energy
range 500 MeV - 500 GeV, which we split into 5 energy bins (0.5-1 GeV, 1-2
GeV, 2-7 GeV, 7-20 GeV, 20-500 GeV). We do not use data below 500 MeV because
at lower energies the angular resolution, as determined by the 68% containment
of the point spread function (PSF), becomes worse than $1^{\circ}$. As the GCE
emission peaks around 2 GeV we chose finer binning around that energy and we
combine energy bins above 20 GeV in a single bin, due to the lower statistics
at higher energies.

The energy bins outlined above have been obtained by resuming a finer native
energy binning of 30 logarithmically spaced bins. This resummation
prescription is also applied to all simulated data in order to properly sample
the LAT’s instrument response functions. Moreover, we use energy dispersion
information for the simulation of our templates by using the flag
edisp_bins=-1. This flag means that the energy dispersion correction only
operates on the input spectra by adding an extra energy bin (of a width
consistent with the initial specification) to the lower and upper end of the
originally defined energy range.

###  2.2 Emission components and Model 1 set up

####  2.2.1 Conventional astrophysical gamma-ray emission components

Modeling the diffuse emission in the region of the Galactic center is a
complex task [14, 30, 63]. Here we adopt one of the IEMs from Fermi LAT’s 1st
SNR catalog [80] that was used to estimate the systematic uncertainty of the
IEMs in that work. The models are all based on GALPROP calculations that are
detailed in [1] but have been adjusted to match the LAT data. The chosen IEM
is based on a distribution of the CR sources which is assumed to be traced by
the distribution of pulsars (Lorimer [81]) with the height of the CR
propagation halo $z=10$ kpc and spin temperature of atomic hydrogen gas
$T_{S}=150$ K. In this model, the gas-related gamma-ray emission (CO and HI
components) is split into four Galactocentric rings (0–4 kpc, 4–8 kpc, 8–10
kpc and 10–30 kpc), allowing for additional freedom in the IEM by scaling them
independently. In addition, these IEMs include an inverse Compton (IC) gamma-
ray component induced by leptonic cosmic rays colliding with the cosmic
microwave background and interstellar radiation fields.

The IEM is supplemented by an isotropic gamma-ray background model (chosen
w.r.t. to the Fermi-LAT IRF file (P8R3$\\_$ULTRACLEANVETO$\\_$V2) and
conversion type
(FRONT)333iso$\\_$P8R3$\\_$ULTRACLEANVETO$\\_$V2$\\_$FRONT$\\_$v1.txt.

In addition to the IE and the isotropic background, we include emission from
the base of the Fermi bubbles [82]. The template of the non-uniformly
spatially structured, data-driven characterization of the bubbles is shown
among other components in Fig. 2. This emission is very uncertain close to the
Galactic plane and at the same time can have a significant impact on the
results, as already argued in [14, 70]. We adopt the same LogParabola spectrum
[82]

$\frac{\mathrm{d}N}{\mathrm{d}E}=F_{0}\left(\frac{E}{E_{0}}\right)^{-\alpha-\beta\ln\\!{\left(E/E_{0}\right)}}$
(2.1)

for the base of the FBs, with parameters $F_{0}=5\times
10^{-10}\;\mathrm{ph}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}\,\mathrm{MeV}^{-1}$,
$\alpha$ = 1.6, $\beta$ = 0.09, $E_{0}$ = 1 GeV.

We also explicitly include all sources listed in the fourth Fermi-LAT source
catalog DR2 (4FGL-DR2) [83], which is built using the same data period we
adopt here. We model sources within a radius of 20∘ around the GC via two
templates; one covering sources within a radius of 5 degrees from the GC and
another incorporating the rest of the sources. We use the best-fit 4FGL-DR2
source spectra as given in the catalog as our baseline. Note that we do not
supplement the 4FGL-DR2 sources with contributions from the unresolved part of
the respective population, which would add to the diffuse background modeling.
Such a contribution was studied in [84], where it was shown that it could be
sizeable, though its properties are highly uncertain. To answer our initial
science question, we however include the contribution from unresolved MSPs in
the GC region (as described below).

These templates are part of the definition of our Model 1-3 setups. We show
the data’s and templates’ two-dimensional spatial morphology and spectra in
Fig. 2 while Fig. 3 displays the templates’ longitude and latitude profiles in
comparison to the selected gamma-ray data set. We comment on the exact
composition of our $\gamma$EMs in Sec. 2.3 (Model 1), Sec. 5.1 (Model 2) and
Sec. 5.2 (Model 3) guided by the principle of increasing the $\gamma$EM’s
complexity after inspecting the performance of the neural network in terms of
its parameter prediction.

Figure 2: Compilation of the spatial (_upper panel_) and spectral (_lower
panel_) morphology of all gamma-ray templates used in our baseline setup to
model the gamma-ray emission in the GC region. The upper panel’s first image
displays the _Fermi_ -LAT data in our ROI between 1 GeV and 2 GeV, which is
the same energy bin chosen for the remaining templates. The templates are the
output of the Fermi Science Tools routine gtmodel and hence display the
expected events from the respective flux model for the given Fermi-LAT
observation time in the infinite statistics limit. The color indicates the
base-10 logarithm of the number of expected gamma-ray events per spatial
pixel. The spectral properties of the DM and MSP templates follow the best-fit
results for Model 2A as stated in the text (c.f. Fig. 10, which fixes
$\gamma=1.25$ in Eq. 2.2. The adjacent MSP template is based on the same
spatial profile whereas the spectral parameters read $\sigma_{L}=0.76$ and
$F_{\mathrm{MSP}}=4.1\times 10^{-7}\;$ph/cm2/s.

Figure 3: Longitude (_left panel_) and latitude profile (_right panel_) of the
selected _Fermi_ -LAT data set and the gamma-ray templates used in our
baseline setup to model the gamma-ray emission in the GC region. The templates
are the output of the Fermi Science Tools routine gtmodel and hence display
the expected events from the respective flux model for the given Fermi-LAT
observation time in the infinite statistics limit. The total contribution of
all four gas rings (CO, orange) and (HI, blue) is shown is a thick line, while
individual gas rings are displayed as dashed lines in shades of the respective
color. For definiteness, we have selected the second energy bin (1 GeV to 2
GeV) for all shown components. The relative normalisation parameters of all
$\gamma$EMtemplates as well as the spectral properties of the DM and MSP
templates in particular follow the best-fit results for Model 2A as stated in
the text (c.f. Fig. 7, which fixes $\gamma=1.25$ in Eq. 2.2. The adjacent MSP
template is based on the same spatial profile whereas the spectral parameters
read $\sigma_{L}=0.76$ and $F_{\mathrm{MSP}}=4.1\times 10^{-7}\;$ph/cm2/s. For
comparison, we show as a black line the longitude and latitude profiles of the
sum of all templates according to their inferred values in Model 2A.

####  2.2.2 Modelling of the GCE

On top of the conventional astrophysical emission, we also add templates for
the GCE. As we aim to distinguish between the DM and MSP (point source) origin
of the GCE we prepare two sets of GCE templates that either reflect the
expectations for DM pair annihilation or a population of faint millisecond
pulsars in the Galactic center. We assume the spatial profile of the GCE
emission is identical for both the DM and the MSP templates. We do this in an
attempt to keep the number of parameters minimal and focus on the algorithm’s
performance to distinguish between truly diffuse and templates comprised of
many point sources, all other ingredients being equal. We note that this
assumption might not be realistic. Indeed, there are indications that the
morphology of the GCE follows the box-shaped bulge emission closely, which
would be a strong argument for the MSP origin of the excess [49, 67], however,
we leave the exploration of this effect for follow-up work. Despite being a
significant signal in Fermi-LAT data of the GC, mismodeling of the GCE has
only a minor impact on the overall severity of the reality gap of our
$\gamma$EMs as it is mainly created by our insufficient knowledge of the
bright astrophysical components. Hence, we impose spherical symmetry of the
GCE and test whether the resulting signal should preferably be smooth or
granular given the Fermi-LAT data set as done in other works applying ML to
the GCE puzzle [71, 70, 72].

Spatial profile. Motivated by the results of previous studies of this region
we adopt a spherically symmetric profile and utilize the so-called generalized
NFW DM density profile [85, 86, 87] for the density distribution,

$\rho_{\mathrm{gNFW}}(r)=\frac{\rho_{s}}{(\frac{r}{r_{s}})^{\gamma}(1+\frac{r}{r_{s}})^{3-\gamma}},$
(2.2)

where $r_{s}=20$ kpc, and $\rho_{s}$ is such that
$\rho_{\mathrm{gNFW}}(r_{\odot})=0.4$ GeV cm-1 at $r_{\odot}=8.5$ kpc and we
treat the inner slope $\gamma$, as a free parameter. Previous fits to the GCE
have resulted in best-fit values $\gamma\in\left[1.1,1.3\right]$ [30, 31, 32],
however, we use a significantly wider range of this parameter
$\gamma\in\left[0.8,~{}1.3\right]$ to avoid potential biases. We do that
guided by the preliminary work [69], which demonstrated the sensitivity of the
outcomes with respect to the $\gamma$ parameter.

Dark matter template. For the DM gamma-ray spectrum, we assume the
$b{\bar{b}}$ annihilation channel (as implemented in the Fermi Science Tools
table). We assume a DM particle candidate of mass $m_{\mathrm{DM}}=42$ GeV and
thermally averaged, velocity-weighted self-annihilation cross section
$\langle\sigma v\rangle=2.2\times 10^{-26}$ cm3s-1 according to the best-fit
in [24]. However, we treat $\langle\sigma v\rangle$ as a free parameter by
introducing an additional normalization parameter that may vary between
$10^{-3}$ and 4 (using a log-uniform prior) to sufficiently sample the
expected total GCE luminosity given the assumed range of slope parameters
$\gamma$.

Millisecond pulsar template. To generate the MSP template together with the
spatial distribution (adopted as above) we need to specify the spectral shape
and the luminosity function of MSPs.
MSP gamma-ray spectrum: We chose the parameter shape as obtained in the
analysis of the stacked MSP spectrum derived in [88],

$\frac{dF}{dE}=F_{0}\left(\frac{E_{\gamma}}{1{\rm
GeV}}\right)^{-1.46}\exp{\left(-\frac{E_{\gamma}}{3.6\,{\rm GeV}}\right)}.$
(2.3)

$F_{0}$ is determined by the luminosity of the specific MSP drawn from the MSP
luminosity function, as specified in the following.

MSP luminosity function: We assume a log-normal luminosity distribution
(following [77], consistent also with earlier works, e.g. [89, 36, 90, 91,
92]), as

$\frac{dN}{dL}\sim\frac{1}{L}\exp\left[\frac{-\left(\log_{10}L-\log_{10}L_{0}\right)^{2}}{2\sigma_{L}^{2}}\right],$
(2.4)

where $L$ is the integrated luminosity in the 0.1 - 100 GeV range, and
$\log_{10}L_{0}=32.61$ (as found as the best fit in [77]). We simulate MSPs in
the luminosity range $L\in\left[10^{30},10^{37}\right]$ erg/s (which
encompasses the expected fluxes of objects placed at the Galactic center
distance that are well below the detection threshold of the Fermi-LAT). We
leave $\sigma_{L}$ as a free parameter with values $\sigma_{L}\in[0.3,1.2]$.
The selected range for $\sigma_{L}$ is larger than the $5\sigma$ confidence
interval around the best-fitting value reported in [77]. With this broad prior
we exhaustively sample the MSP luminosity function as to generate populations
that are almost entirely below the Fermi-LAT’s detection threshold or that
produce a large number of bright pulsars. Hence, our training data contains
instances of the fringe case where a population of many dim sources is
indistinguishable from smooth emission whose importance has been pointed out
in [72].
MSP map flux threshold: The treatment of close-to-threshold MSPs in the
training data set is obviously critical for our approach which aims to detect
faint (sub-threshold) sources. We, therefore, explore two different ways
(scenarios) in which we treat the high flux threshold in our MSP templates.

The main approach adopted in the text (‘Scenario A’), is to generate MSPs
below the flux value defined as the pulsars detection threshold in the 3PC
Fermi-LAT catalog [93], rescaled by the relative exposure. In particular, the
3PC catalogue is derived using the 4FGL procedure and data set (8 years of
data), while we use the 10-year data in this work. We therefore use the 3PC
threshold rescaled with a factor of $\sqrt{10/8}$, resulting in the flux
threshold values in the range $[2.0,4.7]\times 10^{-9}$ GeV s-1 cm-2, in our
region of interest. Since we do not match the MSP luminosity function to the
number of unassociated sources in 4FGL, the issue with this approach is that
it effectively breaks the luminosity function at the detection threshold,
making it impossible for the network to recover the corresponding parameters
(e.g. $\sigma_{L}$) and potentially biasing our results. In addition, the 3PC
sensitivity was derived using different data selections, IRFs and diffuse
models, and might not be fully applicable to our setup. Given how impactful
the treatment of subthreshold sources is expected to be on our results, we
also explore an alternative scenario (‘Scenario B’, further discussed in
Appendix B) in which we do not apply any flux threshold cut in the MSP
template maps. The issue with this approach is that we risk modeling (a
portion) of bright sources twice, as some of the bright MSPs might already be
part of the unidentified source population in 4FGL. As we will see, the
results in both cases are quite similar, strengthening our general
conclusions.

MSP population generation: We apply the following scheme to generate a
realization of individual MSPs:

  * •

Randomly draw a total GCE luminosity $F_{\mathrm{MSP}}$ compatible with the
GCE flux uncertainty range as visualised by the blue band in Fig. 15 of [14]
as the cumulative flux of the MSP population. Consequently, the total range of
the MSP template’s cumulative luminosity is
$F_{\mathrm{MSP}}\in\left[0.02,1.58\right]\times
10^{-6}\;\mathrm{ph}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}$.

  * •

Draw individual pulsars: (i) draw ($\ell$, $b$) from a gNFW profile, (ii) draw
an MSP luminosity from the MSP gamma-ray luminosity function, (iii) derive
$F_{0}$ corresponding to the drawn luminosity, (iv) check if the resulting
flux is less than the Fermi-LAT detection threshold (this step is omitted in
Scenario B).

  * •

Stop the MSP population generation when the randomly drawn total GCE
luminosity is saturated by the cumulative emission from all accepted
_subthreshold_ MSPs.

Fig. 4 illustrates the spread of the expected characteristics of the generated
MSP populations. For definiteness, we fix $F_{\mathrm{MSP}}=0.591\times
10^{-6}\;\mathrm{ph}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}$ and $\gamma=1.28$ to
display the number of sources per photon flux $S$ (0.1 GeV - 100 GeV) for
multiple values of $\sigma_{L}$. The shown source-count distributions
highlight the fact that our choice of priors guarantees to generate enough
variation in the MSP populations. The figure illustrates that for a
sufficiently large part of the $\sigma_{L}$ parameter space, there is a
significant number of MSPs just below the MSP detection threshold (denoted by
the grey band taken from the 3PC catalog in our ROI). At the same time, our
approach covers the extreme case of an MSP population consisting entirely of
dim sources below the detection threshold, which is almost degenerate with a
smooth, DM-like GCE emission. Further details about the generated MSP
populations are given in App. B.

Figure 4: Numbers of generated MSPs as a function of their photon flux $S$
(0.1 - 100 GeV) with respect to the GCE parameters
$F_{\mathrm{MSP}}=0.591\times
10^{-6}\;\mathrm{ph}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}$ and $\gamma=1.28$.
The grey band marks the MSP detection threshold (taken from 3PC catalog and
rescaled to our data set).

###  2.3 Definition and generation of Model 1 training data

Our first $\gamma$EMiteration is comprised of a minimal set of templates
derived from those described in the previous Sec. 2.2.1, i.e. it exhibits a
rather limited number of free parameters, to closely follow the approach of
our proof-of-principle work [69]. In what follows, we refer to our $\gamma$EMs
and scenarios as ”Model XY“ where
$\mathrm{X}\in\\{\mathrm{1},\mathrm{2},\mathrm{3}\\}$ (denoting the
$\gamma$EMiteration) and $\mathrm{X}\in\\{\mathrm{A},\mathrm{B}\\}$ (denoting
the MSP scenario).

We derive the best-fit diffuse model via a maximum likelihood fit utilizing
the astrophysical templates from Sec. 2.2.1 and confront them with the gamma-
ray data in the chosen ROI centered on the GC (c.f. Sec. 3.2 for the technical
details of the likelihood method). To prevent the astrophysical gamma-ray
components degenerate with the GCE from partially accounting for it – thus
reducing its overall luminosity and/or spatial shape – we add and fit the GCE
component during the likelihood maximization. The GCE is represented solely by
a DM-like template whose normalization and slope parameter $\gamma$ are free
to vary.

We combine the templates of the diffuse gamma-ray emission (CO, HI, IC,
isotropic background and FBs) according to their best-fit normalization
parameters in a single template, $N_{\mathrm{DIFF}}$, and repeat the same
procedure for the two templates containing the 4FGL-DR2 sources, resulting in
the template $N_{\mathrm{PS}}$. Consequently, Model 1 features six distinct
parameters denoted as
$\left\\{A_{\mathrm{DIFF}},A_{\mathrm{PS}},A_{\mathrm{DM}},\gamma,F_{\mathrm{MSP}},\sigma_{L}\right\\}$,
whose values are sampled within the ranges reported in Tab. 1 from the prior
distributions found in the last column of the same table and combined
according to Eq. 2.5

$\displaystyle N_{\mathrm{mock}}$ $\displaystyle=$ $\displaystyle
A_{\mathrm{DIFF}}N_{\mathrm{DIFF}}+A_{\mathrm{PS}}N_{\mathrm{PS}}+A_{\mathrm{DM}}N_{\mathrm{DM}}(\gamma)+N_{\mathrm{MSP},\gamma}(F_{\mathrm{MSP}},\sigma_{L},\gamma),$
(2.5)

where $N_{X}$ represents the template associated to the respective
$\gamma$EMparameters. Since we derived $N_{\mathrm{DIFF}}$ and
$N_{\mathrm{PS}}$ from a fit to the data, their priors are chosen closely
around 1 while the choice of the GCE parameters’ prior distributions is
explained above in Sec. 2.2.2. This way, we can generate mock sky images
$N_{\mathrm{mock}}$ that serve as training data for our neural network.

As a further step in the training data generation we impose a constraint that
the total photon counts in the first four energy bins of our mock sky images
_without_ the GCE components must be within $\pm 10\%$ of the total counts of
the _Fermi_ -LAT data in the respective energy bin. In total, we generated 1
million images for both scenarios A and B within this setup of Model 1.

Table 1: Summary of the parameters of Model 1 and their assumed ranges. These parameters are sampled from the prior distributions stated in the last column of this table and combined according to Eq. 2.5 to form a set of training data for our neural network approach. Note that an additional constraint on the total number of counts is applied a posteriori to the templates (see text for more details).  Parameter | Minimum | Maximum | Prior
---|---|---|---
$A_{\mathrm{DIFF}}$ | 0.5 | 2.0 | uniform
$A_{\mathrm{PS}}$ | 0.9 | 1.1 | uniform
$A_{\mathrm{DM}}$ | $10^{-3}$ | 4 | log-uniform
$\gamma$ | $0.8$ | $1.3$ | uniform
$F_{\mathrm{MSP}}$ | $0.02\times 10^{-6}$ | $1.58\times 10^{-6}$ | log-uniform
$\sigma_{L}$ | $0.3$ | $1.2$ | uniform

In this work we aim to distinguish between a DM-like (or “smooth”) and an MSP-
like (or “granular”) GCE. We, therefore, introduce an additional quantity,
$f_{\mathrm{src}}\in\left[0,1\right]$, on which we also train the neural
network. It helps characterize the composition of the emission associated with
the GCE. It is defined as the fraction of the GCE that is due to MSPs, i.e.

$f_{\mathrm{src}}=\frac{F_{\mathrm{MSP}}}{F_{\mathrm{MSP}}+F_{\mathrm{DM}}}\rm{,}$
(2.6)

where $F_{\mathrm{DM}}$ is the integrated gamma-ray flux (0.1 GeV - 100 GeV)
associated with the DM template for a given value of the normalization
parameter $A_{\mathrm{DM}}$. The value of $F_{\mathrm{DM}}$ hence follows from
integrating the differential gamma-ray flux formula for pair-annihilating
(Majorana) DM (see e.g. [94])

$\frac{d\Phi_{\gamma}}{d\Omega\,dE}(E,\psi,\gamma)=\left(\frac{1}{4\pi}\int_{\mathrm{l.o.s}}d\ell(\psi)\rho_{\mathrm{gNFW}}^{2}(\mathbf{r},\gamma)\vphantom{\frac{\langle\sigma
v\rangle_{\mathrm{ann}}}{2m_{\chi}^{2}}\sum_{f}B_{f}\frac{dN_{\gamma}^{f}}{dE}}\right)\left({\frac{\langle\sigma
v\rangle_{\mathrm{ann}}}{2m_{\mathrm{DM}}^{2}}\sum_{f}B_{f}\frac{dN_{\gamma}^{f}}{dE}}\right)$
(2.7)

with respect to the chosen energy range and ROI size. The ROI size is relevant
for the first term in parenthesis, the so-called $J$-factor, which is the
line-of-sight integral of the squared DM density distribution from a selected
portion of the sky. The second term in parenthesis is mostly fixed by our
choice of DM mass and annihilation channel while $\langle\sigma
v\rangle_{\mathrm{ann}}$ follows from the best-fit value of the normalization
parameter $A_{\mathrm{DM}}$. Since $f_{\mathrm{src}}$ is entirely computed
from a subset of the six $\gamma$EM parameters, it does not contain additional
information that the network may use to improve its predictions. However,
including $f_{\mathrm{src}}$ in the training parameter set has the advantage
that the network is able to make an assessment of the uncertainty of this
quantity on its own without further human intervention in contrast to, for
instance, error propagation using the estimated uncertainties of the
ingredients of Eq. 2.6.

The power of $f_{\mathrm{src}}$ to distinguish between a DM-like and MSP-like
GCE depends on the number of MSPs just below the detection threshold, which is
determined by $\sigma_{L}$. If the bulk of the flux comes from very faint
sources, the two templates become indistinguishable as, for example, explored
in the context of neural networks in [72]. As visualized in Fig. 4 our choices
for the MSP-related prior ranges are such that the generated populations cover
the latter extreme case as well as realization with enough sources near the
Fermi-LAT detection threshold. This is also visually demonstrated in Figure 2
which shows an example of an MSP template using the network’s predicted
parameters.

##  3 Inference techniques

While in this work we focus on the analysis technique based on computer
vision, we also compare our findings with a traditional statistical inference
method, namely the maximum likelihood template fitting technique.

In this section, we present the general framework of both analysis setups.

###  3.1 Neural network setup

For our setup, we use the deep ensembles technique [95], which enables us to
estimate the uncertainty in the network predictions. This is done through
random initializations of neural networks that allow for exploring different
modes in the network model space. This method has been shown to perform very
well in terms of accuracy and uncertainty [95].

The network architecture is as follows:

  * •

The input is normalised using the batch normalisation algorithm (see [96]),
which performs a normalisation step (scaling to zero mean and unit variance)
as part of the neural network architecture. This leads to more stable and
faster network training and less dependence on preprocessing;

  * •

Five convolutional “blocks", where one block comprises the following layers: a
convolutional layer with $X$ channels and a kernel size of 3, a convolutional
layer with $2X$ channels, a max-pooling layer with a pooling value of 2 and
the number of channels $X$ of the first convolutional layer is 8;

  * •

The last layer (an image that is 3x3 pixels and has 256 channels) is flattened
to a vector with 2304 elements;

  * •

For every output prediction, a subnetwork contains three dense layers of 128,
64 and 32 neurons and then two output neurons: one that predicts the output
mean, and one that predicts the logarithm of the output variance.

####  3.1.1 Training setup

The neural network was trained on Tesla V100 GPUs splitting the data in 80%
for training, 10% for validation and 10% for testing. The average training
time for this setup is around a few days, while a single prediction of the
trained network is in the order of hundreds of milliseconds. The network is
trained using the Adam optimizer with a learning rate of 0.01 [97], where the
learning rate is halved when there is no improvement in loss on the validation
set for three epochs. Training is stopped after ten of these halvings. The
source code to train the network is available on
https://github.com/rruizbazan/GalacticCenterNet.

####  3.1.2 Deep Ensembles

As mentioned above, to estimate the uncertainties of the neural network
predictions, we adopt the deep ensembles method [95]. To do this, we train a
set of five neural networks, each with a random initialisation of the network
parameters, along with a random shuffling of the data following the
prescription in [95].

To capture uncertainty, the networks predict two outputs per output parameter:
the predicted mean and its variance. Assuming that the data come from a
Gaussian distribution, we define the loss function as follows

$\mathcal{L}=\frac{1}{2}\cdot\sum_{i}^{T}e^{-\hat{z_{i}}}\cdot\left(\hat{y_{i}}-y_{i}\right)^{2}+\frac{1}{2}\cdot\hat{z_{i}},$
(3.1)

where $T$ is the number of output predictions, $\hat{z}$ is the predicted
logarithm of the variance, $\hat{y}$ is the predicted output value, and $y$ is
the true output value. Note that this loss function is equivalent to the
logarithm of a Gaussian (by substituting $z=\textrm{log}(\sigma^{2})$). When
$\sigma=1$, this equation reduces to a standard mean squared error loss
function.

For the final ensemble predictions, the mean and variance are derived from a
mixture of Gaussian distributions. Therefore, the combined mean and variance
of the predictions are $\hat{y}_{*}=N^{-1}\sum_{j}\hat{y}_{j}$ and
$\sigma^{2}_{*}=N^{-1}\sum_{j}(\sigma^{2}_{j}+\hat{y}^{2}_{j})-\hat{y}_{*}$
respectively, where $N$ is the number of ensembles which in our case is five.

###  3.2 Maximum likelihood approach

We employ the standard Poisson likelihood function

$\mathcal{L\\!}\left(\left.\bm{\mu}\right|\bm{n}\right)=\prod_{i,j}\frac{\mu_{ij}^{n_{ij}}}{\left(n_{ij}\right)!}e^{-\mu_{ij}}$
(3.2)

for binned $\gamma$EM data (our input templates, except the MSP one)
$\bm{\mu}$ and experimental data $\bm{n}$ (where the index $i$ runs over the
energy bins while the index $j$ enumerates the spatial pixels of our
templates). We take the best-fit parameters of our $\gamma$EM to be the
maximum likelihood estimators obtained via minimization of
$-\log{\mathcal{L}}$ based on the iminuit python package [98] and the Minuit2
minimization algorithm it provides.

In what follows, we will present the results of the likelihood approach in
combination with the corresponding predictions by the network. To stress it
again, the preparation of the diffuse and point-like source templates in the
definition of Model 1 relies on a prior maximum likelihood fit to obtain the
best-fit parameters for the templates (except for the MSP template) shown in
Fig. 2.

##  4 Results with the Model 1 set-up

###  4.1 Application to simulated data

We assess the trained network’s performance by running it on the validation
data set as a subset of the training data. The results are shown in Fig. 5.
The uncertainty bands mark the (1, 2) $\sigma$ containment intervals, while
the dark grey line marks the median of the prediction $y_{p}$. The uncertainty
bands are solely calculated from the scatter of the network’s (mean)
prediction for image samples that share the same true value $y_{t}$ of a
particular $\gamma$EM parameter (see Fig. 17 for the corresponding results in
scenario B). The uncertainty of predictions (as given by the DE networks) is
visualized for a set of randomly chosen data points (black) with error bars.
The size of the error bar reflects the epistemic and aleatoric uncertainty
summed in quadrature.

The results demonstrate that the network is able to predict the bright
astrophysical gamma-ray components with high precision. This precision is at
the per cent level for the diffuse template (comprising all contributors to
the diffuse emission in the GC). Such a result is expected on Monte Carlo data
since the bulk of the observed photons is either due to detected sources or
the Galactic diffuse emission. Moreover, the fact that we use Monte Carlo data
for the validation ensures perfect agreement between the $\gamma$EM and mock
data.

However, this accuracy is not matched by the predictive power regarding the
GCE components, which are reconstructed with substantial uncertainty while
retaining some predictive power in parts of the parameter space. We find that
the network is capable of consistently retrieving the injected normalization
of the DM and MSP components on average. This statement also applies to most
of the tested range of the gNFW’s inner slope parameter $\gamma$ except for
values at the lower boundary of the sampled prior distribution. As anticipated
however, the width of the MSP’s luminosity function $\sigma_{L}$ cannot be
reliably reconstructed by the network, as scenario A only features MSPs below
the nominal LAT detection threshold, which leaves very little information to
the network to establish a relation between the drawn MSPs in a template and
their underlying luminosity function. The corresponding results on Model 1
Scenario B corroborate this interpretation. In this scenario, $\sigma_{L}$ can
be predicted (with large uncertainties) for at least a fraction of the probed
prior range444This may be put in perspective using Fig. 16 in Appendix B where
we show the number of expected MSPs below and above the threshold for a
selection of GCE model parameters. Values of $F_{\mathrm{MSP}}$ where the
reconstruction of the width of the luminosity function becomes feasible seem
to be already at odds with the experimentally observed number of point-like
sources in the GC region so that the approach we employ here is most likely
not suited to shed light on the physical parameters of a (potential) MSP
population in the MW’s bulge..

The prediction of $f_{\mathrm{src}}$ is also characterized by sizeable
uncertainties, especially for values in the mid-range of this parameter. In
comparison to the preliminary work in [69], we obtain reasonable and even
comparable results in terms of the predicted uncertainty. Overall, the
containment bands derived in the present work are somewhat wider but they have
been obtained with the increased complexity of our training data, i.e.
introducing an explicit dependence on $\gamma$ and $\sigma_{L}$. Hence, the
results of this work corroborate our previous findings. Besides, we can
provide an estimate of the intrinsic uncertainty for each predicted value, a
piece of information that was inaccessible in [69]. We find large errors of
$y_{p}$ except at both endpoints of the probed parameter range. There are
further arguments why the proof-of-principle study was able to achieve a
smaller uncertainty on $f_{\mathrm{src}}$: _(i)_ The previous work used the
3FGL catalog based on 4 years of data, but applied it on data collected over 7
years. This inconsistency in the data could affect the flux of the MSP
template because many of the sources that were unresolved in 4 years of data
become resolvable with the increased exposure. _(ii)_ The prior range for
$\sigma_{L}$ was also different and included more samples with many bright
MSPs. These two effects might have affected the predictive power of the
$f_{\mathrm{src}}$ parameter, which was also directly predicted by the network
as the only trained parameter in the previous work.

Figure 5: Per panel: Prediction of the neural network $y_{p}$ (y-axis)
compared to the true parameter value $y_{t}$ (x-axis) based on the validation
data set as part of the training data generated for Model 1A. The observed
scatter of the predicted values $y_{p}$ is portrayed as grey-shaded bands
whose opacity (dark to light) denotes the $1\sigma/2\sigma$ containment while
the solid black line traces the respective median value. The black data points
are randomly chosen from the validation data set with the length of the error
bar reflecting the uncertainty of the prediction including the epistemic
uncertainty of the network derived via the deep ensembles method and the
aleatoric uncertainty summed in quadrature. In addition, we display the
network’s prediction on the real Fermi-LAT data as red points while the best-
fit results of the likelihood analysis on the real data are shown as orange
points. The green dotted line marks the case of perfect reconstruction; a
deviation from this line indicates a bias in the network’s prediction.

###  4.2 Application to the Fermi-LAT data

We applied the network trained on Monte Carlo data generated from Model 1A on
the Fermi-LAT data and show the findings in Fig. 5 as red points555Since the
true value, in this case, is not known, the point is placed on the diagonal.
(see also Fig. 17). The uncertainty of the predicted values is determined by
the network as described in Section 3.1. Because the uncertainty in general
depends on the explicit position in the high-dimensional parameter space, the
error bar is not identical in all cases to the size of the 68% containment
band of the validation set. In agreement with the validation results, the
algorithm generally finds accurate results for the astrophysical templates.
This is in line with the traditional, likelihood fits (orange points, see
Section 6.1 for details).

To further gauge the quality of the $\gamma$EM selected by the algorithm, we
calculated the residuals, i.e. the difference between the selected $\gamma$EM
and the LAT data, and present them in Fig. 6. To this end, we compute the
Poisson statistical significance $\sigma$, which is given by

$\sigma=\mathrm{sgn}\\!\left(\bm{n}-\bm{\mu}\right)\sqrt{\left|2\left(\bm{n}\ln{\frac{\bm{n}+\varepsilon}{\bm{\mu}}}-\bm{n}+\bm{\mu}\right)\right|}\mathrm{,}$
(4.1)

where $\bm{n}$ denotes the LAT data set whereas $\bm{\mu}$ refers to the
$\gamma$EM data following from the network’s best-fit parameters for Model 1A.
The value of $\varepsilon$ has been set to $10^{-100}$ to avoid numerical
inconsistencies when the value of a pixel in the LAT data is zero. Besides, we
resort to a different formula

$\sigma=\frac{\bm{n}-\bm{\mu}}{\sqrt{\bm{\mu}}}$ (4.2)

when the difference between $\gamma$EM and real data is small
($\left|\bm{n}-\bm{\mu}\right|\ll 1$) but $\bm{n}>0$ to reduce the impact of
possibly occurring numerical instabilities.

Figure 6: Left: Poisson statistical significance (c.f. Eq. 4.1) per energy bin
regarding the network’s prediction (c.f. Fig. 5) based on training data
comprised of the templates of Model 1 and scenario A. All images have been
convolved with a Gaussian smoothing kernel of $0.75^{\circ}$ width. Right:
Best-fit $\gamma$EM spectra and fractional spectral residuals regarding the
network’s prediction. Note that the contributions of all four rings per gas
component have been combined in a single item. The error bars of the total
spectral residuals have been computed using standard error propagation rules
using the statistical errors and network prediction errors for the relevant
parameters.

Significant negative residuals are present in the region of the Central
Molecular Zone (CMZ), left from the GC centered approximately on
$(\ell,b)=(2.0^{\circ},0.0^{\circ})$. This is a known issue and is a
consequence of the rapid variations in the $X_{CO}$ factor in the inner Galaxy
that is caused by tidal forces acting on the molecular gas near the GC. These
show up as elongated features in velocity in the CO emission, resulting in an
artificial increase in the integrated emission that is used to estimate the
molecular hydrogen column density.

There are also positive residuals away from the disk and towards the border of
our ROI. These indicate that the templates contributing at high latitudes
might need more freedom. Of those, the FBs template as a part of the diffuse
template $N_{\mathrm{DIFFUSE}}$ is the most uncertain and may warrant the
inclusion of alternatives to the currently assumed shape of the FBs in Model
1.

While the spectral residuals in the first three energy bins are almost
negligible – hence, pointing to a good agreement between LAT data and network
prediction – the consistency is less striking in the last two energy bins.
Here, the best-fitting $\gamma$EM overestimates the gamma-ray emission,
although the discrepancy is only at the per cent level. At these energies, the
majority of detected photons originated from the Galactic plane so these
regions of the sky are providing the bulk of the constraining power. As a
consequence, the normalization of the 4FGL source template is dominated by the
bright sources in the Galactic disc which may deteriorate the goodness-of-fit
regarding point-like sources at higher latitudes. In combination with the
relative scarcity of events at latitudes below and above the plane, such a
finding is within the expectations.

Since the reality gap, i.e. the fact that the original set of templates might
not capture all features of the real data presents the real limitation that
might bias our results, section 5 is devoted to extending the background model
in an attempt to shrink the reality gap and explore the effect it has on our
results for the GCE nature.

##  5 Extending the gamma-ray emission model

The spatial and spectral residuals shown in Fig. 6 indicate that our model of
the gamma-ray emission within our chosen ROI has well-defined limitations. We
therefore iteratively increase the complexity of our $\gamma$EMs to improve
the agreement with the LAT data and to test the robustness of our results.
Here we summarize our modification of the $\gamma$EM which are discussed and
reasoned in more detail in dedicated subsections. We end up iterating twice,
defining the following two new $\gamma$EMs:

  * •

Model 2 is created by freeing all normalization parameters related to diffuse
emission in the ROI previously combined in $N_{\mathrm{DIFFUSE}}$ as well as
those two of the single 4FGL-DR2 template $N_{\mathrm{PS}}$. The treatment of
the GCE components remains unaltered. Consequently, the overall normalization
of all templates shown in Fig. 2 are now individually varied so that the
number of parameters grows to 17666Note that, while we do have five spectral
bins, the spectral shape of each of the templates is fixed.. We create a
training data set analogously to Model 1.

  * •

Model 3 is created in the spirit to limit the extent by which the spectrum of
the IEM templates biases our results. To this end, we add a power-law
normalization for the gas components per ring and likewise for the IC
template. We use the same power-law adjustment for both gas components and
still allow for independent normalization of the HI and CO templates. This
adds five new parameters with respect to Model 2. To reduce the pronounced
overfitting of the interstellar emission to the left of the GC (see Fig. 6) –
which is known to be related to the preparation of the gas maps – we introduce
two new templates as well as a second template describing the FBs. More
details are provided in Sec. 5.2. In total, Model 3 features 24 parameters.

The $\gamma$EMs are summarized in table 2. In addition to applying the
$\gamma$EMs to the LAT data, we also test cross-domain applications of the
trained networks, applying the network trained on Model 1 data to Model 2 data
and also the network trained on Model 2 data to Model 3 data. This allows us
to determine any possible bias due to the additional freedom in the more
complex $\gamma$EMs in a controlled way. Note that, as discussed in Section 1,
we made sure that Model 1 is embedded in Model 2, which again is embedded in
Model 3, to guarantee an appropriate realization of the domain adaptation
approach (see Fig. 1). We stress, again, that our goal is to use the more
complex $\gamma$EMs (rather than many different low-dimensional $\gamma$EMs)
to learn those distributions and representations shared by possible
sub-$\gamma$EMs by means of the network itself.

Scenario A |  Sources above the detection threshold of 4FGL are excluded from the MSP template.
---|---
Scenario B |  All sources included in the MSP template, even those above the detection threshold.
Model 1 |  Reference $\gamma$EM, see Sec. 2.3 for details (6 parameters, 4 linear and 2 non-linear).
Model 2 |  Based on the templates used in the creation of Model 1. All templates are varied independently without prior constraints obtained via a maximum likelihood fit (17 parameters).
Model 3 |  The templates of Model 2 are kept, but now the IEM templates are normalized with a power law and a scale factor. In addition, the $CO_{1}$ template is split into two and there are now two templates for the FBs (24 parameters, see Sec. 5.2 for all details).
Table 2: Summary of the $\gamma$EMs that we explore in this work, see Fig. 1.

###  5.1 Model 2

As discussed above, Model 2 introduces a level of complexity by freeing all
normalization parameters of the components that constitute the astrophysical
background model. This training data generated within the framework of Model 2
no longer include prior information obtained via established statistical
inference methods like a maximum likelihood fit and enables an independent
comparison of the performance of a neural network to traditional methods. Note
that, to obtain the best-fit $\gamma$EM parameters, the network does not
minimize a Poisson log-likelihood function but rather an internal loss
function (see Eq. 3.1) whose relation to the likelihood function is not a
priori known. It is therefore not expected to recover the same best-fit
values.

To train the neural network with Model 2 data we generated 1.8 million images
for both scenarios A and B. The GCE parameter space is sampled from the same
prior distributions specified in Tab. 1 whereas all parameters related to
conventional astrophysics are drawn from uniform priors between 0 and 3. We
show the results of a neural network trained on Model 2 data and applied to
the real data in Fig. 7 for scenario A. The corresponding results for scenario
B are displayed in Fig. 18. The color coding of the panels follows the
respective plots for Model 1. For completeness, the observed spatial and
spectral residuals with respect to both scenarios of Model 2 are shown in
Appendix C since they do not change significantly and thus provide only a
marginal additional value to the discussion.
Despite the potentially different reasoning of a maximum likelihood fit and
supervised deep learning algorithms to derive best-fit parameters, the
correspondence between both methods is quite striking for the majority of the
astrophysical parameters not related to the GCE. The exceptions from this
general trend are those components whose main contribution is not necessarily
falling into the chosen ROI of this study, namely gamma-ray components that
are present at high Galactic latitudes, i.e. the third and fourth CO ring, the
fourth HI ring and the isotropic gamma-ray background. These components are
likewise predicted with a larger uncertainty by the network.
Regarding the GCE parameters, the predictions for Model 2 exhibit larger
uncertainties both with respect to the validation data set and the real LAT
data. Since the complexity of the $\gamma$EM has increased, a growth of the
network’s uncertainty on its own predictions can be understood in this regard.
On top of this, this deterioration of the precision of the network is
reinforced by the aforementioned astrophysical components that are difficult
to learn.
The predictions for the GCE parameters of the likelihood and ML approach are
not the same. As already witnessed in the case of Model 1, the ML approach is
capable of detecting a non-Poissonian component like the MSP population in the
MW’s bulge whereas it may seem impossible for all practical matters for the
likelihood fit777To illustrate this statement: The images have a size of
$120\times 120=14400$ pixels. A typical number of MSPs needed to saturate the
GCE given the reconstructed parameters of Model 2 is around 3000. Hence, there
are $\left(\begin{array}[]{c}17399\\\ 3000\end{array}\right)$ ways to place
these sources in the image (with repetition). In numbers, there are $\sim
10^{3470}$ ways to distribute these sources in the GC images while further
variations like the intrinsic brightness of each of the MSPs are not yet taken
into account.. In comparison to Model 1, the network finds a larger MSP-like
contribution to the GCE while at the same time, the spatial morphology of the
signal is predicted with a steeper slope $\gamma>1.2$ than obtained via the
maximum likelihood method. This in return increases the parameter
$f_{\mathrm{src}}$, which is about $50\%$ in this $\gamma$EM and scenario.
The differences between the likelihood and ML approach are the same in
scenario B but there, we find a smaller MSP contribution to the GCE of around
$20\%$ (although being higher than in Model 1B). Hence, dropping constraints
on the template normalizations imposed in Model 1 leads to an increase in the
reconstructed fractional contribution of unresolved point-like sources to the
GCE.

Figure 7: The same as Fig. 5, but in the context of Model 2 (Scenario A).

###  5.2 Model 3

In our final $\gamma$EM iteration, we enlarge our template selection (see Fig.
8) to address the aforementioned issue with the central CO ring. To further
alleviate the positive residuals at high latitudes, we add a template
quantifying the emission of the FBs because this component is the least
constrained and potentially one of the most degenerate with the DM template.
In particular,

  * •

We split the innermost CO ring – whose emission shows significant left/right
asymmetry – in two independent templates. For this rather ad-hoc solution, we
define a gamma-ray flux threshold at 25$\%$ of the maximal value found in the
input gamma-ray flux model of the first CO ring. Each pixel above this
threshold is part of a $N_{\mathrm{CO},1,\mathrm{high}}$ template whereas the
remainder of pixels is cast into an independent template to which we assign
the already used label $N_{\mathrm{CO},1}$.

  * •

We add a second FBs template $N_{\mathrm{FBcatenary}}$, taken from [99]. This
template is built by extracting residual intensity maps from the LAT data,
after the application of the Fermi diffuse model. The edges of the bubbles’
residuals are found to be well reproduced by two catenary curves (see [99] for
parameter values), which also reproduce correctly the structures observed
close to the GC in the ROSAT X-ray observations. While this template is quite
similar in morphology to the original FB template from Model 1, it is centered
closer to the GC. An additional difference is that we assume that it exhibits
a uniform intensity while the spectrum is the one that we use for the FB
template in Model 1.

  * •

We vary the spectral indices of four gas components, in such a way that we
keep the same spectral index per ring for CO and HI, but allow for their
normalization to change independently. Similarly, we also allow for variations
of the spectral index of the IC component. In all of the cases we allow for
the change of the spectral index in the $\Delta\gamma\in\left[-0.3,0.3\right]$
range, where we define the new spectral index as $\gamma-\Delta\gamma$, i.e.
negative values of $\Delta\gamma$ imply spectral softening.

The normalizations of these two templates yield two additional parameters in
our fits, while the spectral indices add another five parameters, see Table 2.

Figure 8: Compilation of the additional gamma-ray templates used in Models 2 and 3. The templates are the output of the Fermi Science Tools routine gtmodel and hence display the number of expected events from the respective flux model for the given Fermi-LAT observation time in the infinite statistics limit. The color indicates the decadic logarithm of the number of expected gamma-ray events per spatial pixel. For definiteness, we have selected the second energy bin (1 GeV to 2 GeV) for all shown templates.  Parameter | Minimum | Maximum | Prior
---|---|---|---
$A_{\mathrm{CO,1,high}}$ | 0 | 3 | uniform
$A_{\mathrm{FBcatenary}}$ | 0 | 3 | uniform
$PL_{1}$ | -0.3 | 0.3 | uniform
$PL_{2}$ | -0.3 | 0.3 | uniform
$PL_{3}$ | -0.3 | 0.3 | uniform
$PL_{4}$ | -0.3 | 0.3 | uniform
$PL_{\mathrm{IC}}$ | -0.3 | 0.3 | uniform
Table 3: Random sampling prior distributions for the extra parameters of Model
3.

To train the neural network with Model 3 data we generated 2.4 million images
for both scenarios A and B. The priors for the additional parameters are
stated in Tab. 3 while the remaining parameters are sampled as in Model 2. The
predictions of the neural network trained on data created within the Model 3A
framework are shown in Fig. 9 together with the corresponding maximum
likelihood results, while the residuals can be found in Fig. 21 in Appendix C
(see also Figs. 19 and 24, for Scenario B). To conduct the maximum likelihood
fit, we have used astrophysical templates corrected by the network’s best-fit
spectral indices.

The left-right asymmetry of the central residuals as well as those in the
Galactic plane, are reduced in this iteration (compare Figs. 6 and 21), with
the new CO template normalization ($A_{\mathrm{CO},1,\mathrm{high}}$) found to
be smaller than one, thus reducing the $\gamma$EM intensity in that region.
The algorithm finds the two FBs templates to be present in the data, with
comparable intensities. Regarding the reconstructed spectral indices, we note
that despite the added freedom, most $\Delta\gamma$ values are predicted close
to zero, meaning little spectral correction is required. A noteworthy
exception from this observation is the spectral correction for the first gas
ring, which is close to $\Delta\gamma=0.1$ for both scenarios indicating a
hardening of the default gamma-ray spectrum. Such a spectral hardening might
explain the increase of the residuals in Model 3 with respect to the last two
energy bins (see Figs. 21 and 24).

The parameters for the GCE are still compatible with those from the network
trained on Model 1 data, although there is a small trend towards higher values
of $F_{\mathrm{MSP}}$ and therefore a higher $f_{\mathrm{src}}$ value (Model
2A: $f_{\mathrm{src}}=(46\pm 36)\%$; Model 3A: $f_{\mathrm{src}}=(93\pm
20)\%$).

Figure 9: The same as Fig. 5, but in the context of Model 3 (Scenario A).

##  6 Discussion

###  6.1 Comparison with traditional likelihood template fitting

The main caveat is that when using the likelihood fitting, it is challenging
to include the MSP template(s) in the fitting stage due to their probabilistic
character. In addition, we limit our fitting to linear parameters only (except
for the $\gamma$ parameter), while the rest of the nonlinear parameters are
marginalized over. The nature of uncertainties in the two analysis methods is
different as well: while the likelihood analysis is based on the statistical
(Poisson) uncertainty, the mapping of the network’s aleatoric uncertainty of
the input to the predicted uncertainty is complex and depends on several
parameters (for example on the loss function of the network, see discussion in
Section 3.1). Despite a lack of a one-to-one comparison, we find that the
results obtained from these two methods are consistent, except in a few cases
discussed below.

To gauge the impact of an unmodeled MSP contribution in the likelihood
analysis, we use simulated data based on Model 3 with $\Delta\gamma=0.0$ and
draw in total 500 Poisson realizations where all other parameters are set to
1, except that $F_{\mathrm{MSP}}=5\times 10^{-7}$ ph/cm2/s. In Table 4 we
compare the likelihood analysis results in which the MSPs are present in the
mock data but not in the simulated gamma-ray data (the first column, labeled
’DM+MSP’), with the case where that template is not present in the mock data
(the second column, ’DM only’). In addition, we show the results of the
likelihood fit when the mock data contain only an MSP component which is then
fitted via the smooth DM template. We observe that the majority of parameters
corresponding to bright templates are unaffected by the omission of the MSP
template in the modeling step. However, the flux of the MSP template,
calculated from the mean best-fit results under the assumption that the
original DM setup is known (i.e. $A_{DM}=1$ and $\gamma=1$), and shown in the
last row, is recovered only partially.

| DM + MSP | DM only | MSP only
---|---|---|---
$A_{\mathrm{CO},1,\mathrm{high}}$ | $0.99\pm 0.01$ | $1.0\pm 0.01$ | $1.0\pm 0.01$
$A_{\mathrm{CO},1}$ | $1.0\pm 0.02$ | $1.0\pm 0.02$ | $1.0\pm 0.02$
$A_{\mathrm{CO},2}$ | $1.0\pm 0.01$ | $1.0\pm 0.01$ | $1.0\pm 0.01$
$A_{\mathrm{CO},3}$ | $0.99\pm 0.03$ | $1.0\pm 0.03$ | $0.98\pm 0.03$
$A_{\mathrm{CO},4}$ | $0.44\pm 0.38$ | $0.92\pm 0.44$ | $0.5\pm 0.39$
$A_{\mathrm{HI},1}$ | $1.0\pm 0.02$ | $1.0\pm 0.02$ | $1.0\pm 0.02$
$A_{\mathrm{HI},2}$ | $1.0\pm 0.01$ | $1.0\pm 0.01$ | $1.01\pm 0.01$
$A_{\mathrm{HI},3}$ | $1.01\pm 0.01$ | $1.0\pm 0.01$ | $1.01\pm 0.01$
$A_{\mathrm{HI},4}$ | $1.02\pm 0.08$ | $0.99\pm 0.08$ | $1.04\pm 0.07$
$A_{\mathrm{PS},1}$ | $1.0\pm 0.0$ | $1.0\pm 0.0$ | $1.0\pm 0.0$
$A_{\mathrm{PS},2}$ | $1.0\pm 0.01$ | $1.0\pm 0.01$ | $1.0\pm 0.01$
$A_{\mathrm{ISO}}$ | $0.74\pm 0.17$ | $0.97\pm 0.14$ | $0.92\pm 0.15$
$A_{\mathrm{IC}}$ | $1.03\pm 0.03$ | $1.01\pm 0.02$ | $0.99\pm 0.02$
$A_{\mathrm{FB}}$ | $1.0\pm 0.03$ | $1.0\pm 0.02$ | $1.0\pm 0.03$
$A_{\mathrm{FBcatenary}}$ | $1.06\pm 0.06$ | $1.01\pm 0.06$ | $1.0\pm 0.06$
$A_{\mathrm{DM}}$ | $1.26\pm 0.13$ | $0.98\pm 0.09$ | $0.51\pm 0.1$
$\gamma$ | $1.03\pm 0.02$ | $1.0\pm 0.02$ | $1.0\pm 0.03$
$F_{\mathrm{MSP}}$ |  $2.29\times 10^{-7}$ ph/cm2/s | / |  $2.71\times 10^{-7}$ ph/cm2/s
Table 4: Likelihood fit results for a mock data set based on Model 3 where all
normalization constants in the mock data are equal to unity, except for the
MSP template which features a total flux of $5\times 10^{-7}$ ph/cm2/s as well
as all spectral corrections $\Delta\gamma$ being set to zero. In the header
row of the table, we list the GCE templates that are used to create the mock
data utilized to perform the fit with the results listed in the respective
column. In all cases, the spatial GCE profile follows a gNFW density
distribution with $\gamma=1.0$. Note that the reported statistical error is
the standard deviation of the sample of the respective best-fit parameter
resulting from 500 Poisson realisations of the prepared gamma-ray mock data.
The ’recovered’ MSP total flux value shown in the last row is calculated from
the mean best-fit results for $\gamma$ and the DM normalization under the
assumption that the original DM setup is known (i.e. $A_{\mathrm{DM}}=1$ and
$\gamma=1.0$).

The optimal parameters derived from our likelihood fit to the real LAT data
are incorporated into our findings in Figs. 5, 7, and 9 as orange points. We
note that the likelihood outcomes align with those generated by the network,
particularly for the luminous gas emission components. However, for certain
components, the network predictions fall outside the one-sigma statistical
uncertainty, which is typically narrower in the likelihood outcomes.

The concurrence is less satisfactory for the fainter and less organized IC,
FBs, and ISO components, where the likelihood outcomes reveal a systematically
elevated normalization for the ISO component and a consistently reduced
normalization for the IC and FBs components. While not flawless, the degree of
agreement illustrates that the networks can effectively disentangle the data
into the astrophysical components they were trained on.

In terms of the GCE decomposition (i.e. DM normalization), there is a
significant discrepancy, which is expected, since the MSP component is not
included in the likelihood approach. In particular, the value of
$A_{\mathrm{DM}}$ is systematically larger to accommodate for the missing MSP
component, but also because the value of $\gamma$ is slightly smaller,
affecting the total flux from the smooth DM component.

###  6.2 Robustness of the results

In Fig. 10 we summarize our networks’ results for the parameters controlling
the GCE excess, plotting them for all $\gamma$EMs and scenarios. The parameter
values are mostly consistent considering the large uncertainty, but there
seems to be a trend towards lower values for $A_{DM}$ (or $F_{DM}$,
respectively), and higher values for $F_{MSP}$, $\gamma$, and $f_{src}$ when
going from the lower complexity of Model 1 to the higher complexity present in
Model 3. The value of $f_{src}$ is lower for Scenario B by construction. This
observation is in line with our expectations since in Scenario B a higher MSP-
like luminosity is linked to a larger number of MSPs above the detection
threshold in the simulated gamma-ray images (c.f. Fig. 16 in App. B). The room
in the real data for such artificially injected, bright sources in the GC is
rather limited, thus lowering the value of $f_{src}$. The network thus
consistently finds an excess over our astrophysical background models while
the characterization of the results into a smooth and point-source-like
component depends on the selected background model. It is tempting to conclude
that the trend of an increasing $f_{src}$ will continue as the $\gamma$EM’s
complexity increases, but we caution against that because there is no
guarantee that the increased complexity of the background model did not
adversely affect the characterization of the GCE.

To make a direct comparison with our previous work [69], we retrain the Model
3 networks to predict merely $f_{src}$ while treating the remaining set of
parameters as nuisance parameters, and label such networks ’1D’. We observe
that the predicted network uncertainty significantly shrinks in this case (as
shown in Fig. 10 middle panel), consistently with what is observed in [69].
This finding explains why the old results were more precise but not
necessarily more accurate.

The trend of increasing $f_{src}$ from Model 1 to Model 3 rather suggests that
the astrophysical background modeling is the dominating source of systematic
uncertainty, an issue we focus on in the rest of this section.

Figure 10: Top: Results of the network prediction on the real data in our
three Models and two Scenarios (circle: A; triangle: B) with respect to the
parameters of the GCE. Bottom: Comparison of the network’s prediction for
$f_{\rm src}$ per $\gamma$EM and scenario. 3-1D shows the result of networks
trained only on one ($f_{\rm src}$) parameter (treating the rest as nuisance
parameters) when applied to Models 3.

Outline of the robustness checks. To quantify its impact, we apply the network
trained on simulated data created by one $\gamma$EM to data simulated by
another $\gamma$EM. We do that in two steps. Firstly, we focus on cross-domain
application, for various $\gamma$EMs and scenario combinations of our
$\gamma$EMs and scenarios (see Fig. 1). After that, we extend our test by
applying the networks trained on our $\gamma$EMs, to an ’independent’
verification data set that is created under a range of substantially different
assumptions. For this independent $\gamma$EM, we chose the Galactic
Interstellar Emission Model for the 4FGL Catalog Analysis (so-called Pass 8
Fermi diffuse model, gll_iem_v07.fits)888See
https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/aux/4fgl/Galactic_Diffuse_Emission_Model_for_the_4FGL_Catalog_Analysis.pdf,
which is not based on GALPROP calculation of CR propagation but uses a
different, data-driven approach to calculate gas emissivity (see [100] for
details). It also uses ’patches’ of selected sky regions when deriving the
model, which is derived from the residual $\gamma$-ray intensity from fitting
the overall $\gamma$EM in an iterative procedure. As a consequence, the large-
scale residual emission (including the (smoothed) GCE and FBs) is part of the
diffuse model itself, in stark difference from the IEM our network was trained
on. Despite these differences, the overall emission in the verification data
set (that is built by summing up our usual PS, ISO, FB, and GCE templates
together with this IEM), contains the emission components the network is
expected to be broadly familiar with and should be able to quantify. We note
that the diffuse model gll_iem_v07.fits is not meant to be used for the
analysis of large-scale gamma-ray sources (like the GCE) in real Fermi-LAT
data999See
https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/aux/4fgl/Galactic_Diffuse_Emission_Model_for_the_4FGL_Catalog_Analysis.pdf
for a list of caveats compiled by the Fermi-LAT collaboration.. However, this
caveat is not relevant for the purpose of our robustness check since it is
solely based on Monte Carlo data.

We detail the outcomes of both of these approaches below.

Figure 11: Results of the cross-domain application with respect to validation
data prepared within the framework of Model 2A applied to a fully trained
network whose parameter reconstruction is tuned to Model 1A. (_Top_ :) Same as
the results for Model 2A shown in Fig. 7 reduced to the five parameters
characterizing the GCE. (_Bottom_ :) Performance of the network trained on
Model 1A data ($y_{p}$) when given the validation data set of Model 2A
($y_{t}$). The color coding is adopted from Fig. 7.

Cross domain applications. Since our $\gamma$EMs are nested, it is only
necessary to consider a small subset of all the possible combinations of
$\gamma$EMs and scenarios. For us, it suffices to explore the cross-domain
results obtained from using a network trained on a low number of parameters to
data that exhibits a higher complexity.

The outcome of this cross-check – focusing here on the GCE parameters – is
visualized in Fig. 11. The top row of the plot re-states the results of
network A2 on validation data of the same $\gamma$EM and scenario shown in
Fig. 7. The cross-domain application of the same validation data set to a
network trained on Model 1A data is displayed in the lower row. We observe
that in the lower panel (cross-domain results) the network results are clearly
biased as all prediction means depart from the dotted, green diagonal line
representing perfect reconstruction. Moreover, the predicted uncertainties are
underestimated in a significant number of cases. This shows very clearly that
the networks are very sensitive to the background model in the training data,
demonstrating nicely the relevance of the reality gap. We stress that this
observation is not an inherent limitation of our inference method based on NNs
but rather the impact of the reality gap. A template-based approach using a
maximum likelihood fit would be likewise affected.

Besides the explicit cross-model robustness check in Fig. 11, we have also
quantitatively examined the impact of cross-domain analyses on data from
different MSP scenarios (Scenarios A and B). Here, the induced bias of the
network predictions is much less pronounced than before. The cross-domain
results imply that a $\gamma$EM trained with data from scenario B is fairly
robust against changing the MSP scenario to A. Almost all parameters are
reconstructed without any additional biases. The reverse exercise, i.e.
letting a network trained on scenario A data reconstruct parameters of
scenario B data sets, leads to stronger biases within the $1\sigma$
containment band of the scatter of the predicted values. Hence, the systematic
model uncertainty introduced via the two MSP scenarios has substantially less
impact on the accuracy and precision of the network’s predictions.

New domain applications. We confirm this trend also by applying a network to a
significantly different $\gamma$EM, namely the Pass 8 Fermi diffuse model,
where we observe that the performance of the networks degrades when it comes
to estimating the GCE components. Fig. 12 shows the results for Model 3B,
while the rest of the cases are shown in Appendix D. The bright emission
components remain well predicted also in this case and predictions for some of
the faint templates, namely that of the FBs behave also as expected. In
particular, our networks tend to overpredict the normalization of the FB
template, which is reasonable given the fact that the FBs are now both, part
of the diffuse model itself and also added as a separate template. At the same
time, the networks have a significantly harder time predicting the DM and MSP
components, indicating possible internal biases. In Fig. 12 we observe, for
example, that the MSP and ISO components are significantly over-predicted
while the DM component is under-predicted. This was a priori not expected
because the residual component of the Pass 8 Fermi diffuse model was
deliberately smoothed so it would not contain emission from individual point
sources. We expected that the MSP component would be accurately predicted
while $A_{DM}$ would be over-predicted. A possible reason is that the smoothed
residual in the Pass 8 Fermi model is not perfectly symmetric around the GC
and therefore looks more similar to the MSP component. This corroborates our
conclusion of a lack of robustness of the features that our DENs can
reconstruct due to the presence of the reality gap. This lack of robustness
was also present in all other attempts to determine the nature of GCE as
evident by the large and sometimes contradicting body of scientific articles
on the nature of the GCE. In particular, the history of NPTF results described
in Sec. 1.1 illustrates the impact of background mis-modeling and the presence
of a reality gap in general.

The work we present in this paper goes beyond the exploration of previous
analyses that are subject to the reality gap/background mis-modeling. We lend
more freedom to the astrophysical contribution – a nested multi-parametric
modeling with progressively increasing complexity – to the GCE compared to
[61, 62] where the astrophysical background composition and flexibility are
essentially the same for the entire analysis (except for a few sanity checks).
As seen in the results of the ‘new domain application’, an artificially
injected signal (here: smooth, DM-like) may be mis-attributed to the MSP
contribution in presence of a substantial reality gap depending on the
modeling of the diffuse astrophysical backgrounds. However, we did not observe
a prominent influence of the 4FGL-DR2 templates on this matter. In the
following section, we propose a method to quantify and explore the severity of
the reality gap of a gamma-ray model of the GC region to minimize its impact
in future works on the GCE.

Figure 12: Robustness check of the network trained on Model 3B data by
applying a validation data set that contains the Fermi diffuse background
model instead of the gas and IC templates of Model 3’s set of diffuse
components. The line and color style is adopted from Fig. 5. We only show
those parameters of Model 3B that are in common with the validation data set
containing the Fermi diffuse background model.

###  6.3 Quantifying the reality gap

To quantify the gap between our $\gamma$EMs and a more accurate $\gamma$EM of
the GCE, we used an ensemble of network models [101, 102]. The Deep SVDD model
is an encoding network that is trained to map the input, in our case a [120 x
120 x 5] tensor, to a constant vector with a certain dimensionality. This
vector essentially defines a manifold onto which we compress the data.

The loss is defined as

$s(x)=\left[O_{n}^{d}-{\rm Model}~{}(x)\right]^{2},$ (6.1)

where the model maps the input x to the same tensor shape as the manifold $O$.
In our case, it is a vector of identical scalar values, with the subscript $n$
defining the scalar value and superscript $d$ the number of elements in the
vector.

When optimizing the SVDD model, the loss function between the predictions for
the training data set and the target is minimised. Therefore, data coming from
the same joint distribution can be expected to have a lower reconstruction
loss than data coming from a different distribution. In other words, a loss
out of distribution with respect to the test data would quantify how different
the real data is with respect to our $\gamma$EM.

To verify this, we follow [102] and train an ensemble of SVDD models with
different combinations of target output values $n=\left[10,25,50,100\right]$
and dimensionalities $d=\left[5,89,144,233\right]$. The SVDD architecture used
is the network described in Sec. 3, where the last layer was replaced by a
layer with as many neurons as the dimensionality of the output vector. The ML
model was trained with the Adam optimiser with a learning rate of 0.01,
halving the learning rate if the loss in the validation set did not improve
over five epochs.

Finally, we used the trained SVDD models to predict the distance in a
d-dimensional latent space between the encoded vector for the test data and
the Fermi diffuse model data with respect to the target vector in the case of
Model 3 101010This is equivalent to compute the $L_{2}$ norm or Euclidean
distance between two vectors in an n-dimensional space.. We combined the
results of the trained SVDD models resulting from the combination of n and d
values, by adding them together, as shown in Fig. 13 in the form of
histograms. There one can see how the real data (red) are clearly out of
distribution not only with respect to Model 3, but also to the Fermi diffuse
model. We also checked the consistency of our results by using noise as input,
which leads to the much larger prediction of distances (green line). Note that
the plots are normalized in such a way that the distance to the maximally
different $\gamma$EM is set to the value of 1.

Figure 13: SVDD prediction of the distance in latent space between the encoded
vector for the test data (blue, Model 3A (_left_) and Model 3B (_right_)), the
Fermi diffuse model (orange), the real data (red line) and an image of the
same dimensions filled only with Gaussian noise (green line). Note that the
normalized distance in both inset figures is adjusted to better reflect the
proximity of the real data to the training data compared to the pure noise.

##  7 Summary and Future Prospects

We use DeepEnsemble Networks to perform a detailed gamma-ray analysis of the
complex Galactic center region. By applying the trained networks on the
verification data set we demonstrate that they are capable of recovering
components of background emission and the presence of the GCE. In particular:

  * •

Bright components associated with conventional astrophysical processes are
detected robustly and consistently between our $\gamma$EMs, see Fig. 15. They
are also largely consistent with the prediction from the traditional
likelihood method and exhibit comparable uncertainties for most templates.

  * •

The networks robustly detect the presence of the GCE in all our $\gamma$EMs,
with the properties (flux and spatial distribution) consistent with other
works (see Fig. 14).

  * •

The nature of the GCE, however, while well predicted within each $\gamma$EM,
does not appear to be robust when the networks are applied outside of their
training domain. We conduct a dedicated study of the limitation of our
networks to successfully generalize the GCE nature and conclude that the
reality gap remains the final obstacle in addressing the nature of GCE in our
framework.

  * •

We further quantify the reality gap by adopting Deep SVDD architecture, which
allows us to compare the distance between the encoded vector for our
$\gamma$EMs and the real data rendering it possible to conclude that the real
data is out-of-distribution even for the most complex $\gamma$EM.

  * •

While this lack of robustness proved to be present in many other attempts to
determine the nature of GCE – as evidenced by the large and sometimes
contradicting body of scientific articles on the nature of the GCE – our work
provides the first detailed study of such limitations for nested $\gamma$EMs
aimed at extracting the small-scale clustering of photons in gamma-ray data.
More generally, whenever the NN application is strongly dependent on training
data based on theoretical or data-driven models, it is imperative to study and
quantify the reality gap before the results can be claimed robust.

In summary, this work represents a dedicated study of the often encountered
issue of $\gamma$EMs being an, at best, incomplete description of the full
physics of the GC region. We analyzed the consequences of this reality gap in
the setting of nested $\gamma$EMs with increasing complexity to derive a
quantitative statement about the contribution of point sources to the GCE; an
approach never attempted in the literature. We obtain
$f_{\mathrm{src}}=0.10\pm 0.07$ for Model 1B while Model 3B yields
$f_{\mathrm{src}}=0.79\pm 0.24$. We consider Model 3 as the most reliable
description of reality in our study. However, in light of the persistent
reality gap, we are unable to faithfully quantify the error on
$f_{\mathrm{src}}$ because our network’s predictions do not include the
additional uncertainty arising from the fact that the real data is outside of
the investigated $\gamma$EM space. While our analysis employs NN as the means
to infer the properties of the GCE – as it allows us to efficiently examine a
range of $\gamma$EMs and to account for various levels of uncertainties
associated with background and signal (GCE) components – the implications of
our work hold in more generality irrespective of the chosen inference
technique. We, therefore, argue that bridging the reality gap or understanding
the limitations and biases it introduces, is an imperative step before being
able to derive a verdict on the nature of the GCE. We invite the authors of
all previous studies aiming at assessing the nature of the GCE to verify that
their $\gamma$EM can bridge the reality gap to the LAT gamma-ray sky.

Figure 14: Reconstructed spectrum of the GCE (purple) according to the
network’s prediction for the DM (blue) and MSP (light red) components based on
training data of Model 3 + Scenario A (Left) and Model 3 + Scenario B (Right)
compared with the spectrum as derived in previous analyses [14, 30, 31]. The
colour-shaded bands around each prediction denote the network’s error
estimate. For definiteness, we have fixed $\gamma=1.28$, which is the best-fit
value of our network regarding Model 3 in both scenarios.

## Acknowledgments

The author(s) gratefully acknowledges the computer resources at Artemisa,
funded by the European Union ERDF and Comunitat Valenciana as well as the
technical support provided by the Instituto de Fisica Corpuscular, IFIC (CSIC-
UV). R. RdA is supported by PID2020-113644GB-I00 from the Spanish Ministerio
de Ciencia e Innovación. G. Z. and C. E. acknowledge the financial support
from the Slovenian Research Agency (grants P1-0031, I0-0033, J1-1700 and the
Young Researcher program). C. E. further acknowledges support by the “Agence
Nationale de la Recherche”, grant n. ANR-19-CE31-0005-01 (PI: F. Calore). The
work of C. E. has been supported by the EOSC Future project which is co-funded
by the European Union Horizon Programme call INFRAEOSC-03-2020, Grant
Agreement 101017536.

The Fermi LAT Collaboration acknowledges generous ongoing support from a
number of agencies and institutes that have supported both the development and
the operation of the LAT as well as scientific data analysis. These include
the National Aeronautics and Space Administration and the Department of Energy
in the United States, the Commissariat à l’Energie Atomique and the Centre
National de la Recherche Scientifique / Institut National de Physique
Nucléaire et de Physique des Particules in France, the Agenzia Spaziale
Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry
of Education, Culture, Sports, Science and Technology (MEXT), High Energy
Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency
(JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research
Council and the Swedish National Space Board in Sweden.

Additional support for science analysis during the operations phase is
gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy
and the Centre National d’Études Spatiales in France. This work performed in
part under DOE Contract DE-AC02-76SF00515.

## References

  * [1] Fermi-LAT collaboration, M. Ackermann et al., _Fermi-LAT Observations of the Diffuse Gamma-Ray Emission: Implications for Cosmic Rays and the Interstellar Medium_ , _Astrophys. J._ 750 (2012) 3, [1202.4039].
  * [2] HAWC collaboration, H. Zhou, C. D. Rho and G. Vianello, _Probing Galactic Diffuse TeV Gamma-Ray Emission with the HAWC Observatory_ , _PoS_ ICRC2017 (2018) 689, [1709.03619].
  * [3] H.E.S.S. collaboration, F. Aharonian et al., _Discovery of very-high-energy gamma-rays from the galactic centre ridge_ , _Nature_ 439 (2006) 695–698, [astro-ph/0603021].
  * [4] H.E.S.S. collaboration, A. Abramowski et al., _Diffuse Galactic gamma-ray emission with H.E.S.S_ , _Phys. Rev._ D90 (2014) 122007, [1411.7568].
  * [5] HESS collaboration, H. Abdalla et al., _Characterising the VHE diffuse emission in the central 200 parsecs of our Galaxy with H.E.S.S_ , _Astron. Astrophys._ 612 (2018) A9, [1706.04535].
  * [6] I. V. Moskalenko, T. A. Porter and A. W. Strong, _Attenuation of Very High Energy Gamma Rays by the Milky Way Interstellar Radiation Field_ , _ApJ_ 640 (Apr., 2006) L155–L158, [astro-ph/0511149].
  * [7] T. A. Porter, G. P. Rowell, G. Jóhannesson and I. V. Moskalenko, _Galactic PeVatrons and helping to find them: Effects of galactic absorption on the observed spectra of very high energy $\gamma$ -ray sources_, _Phys. Rev. D_ 98 (Aug., 2018) 041302, [1808.07596].
  * [8] R. Genzel, F. Eisenhauer and S. Gillessen, _The Galactic Center Massive Black Hole and Nuclear Star Cluster_ , _Rev. Mod. Phys._ 82 (2010) 3121–3195, [1006.0064].
  * [9] A. T. Barnes, S. N. Longmore, C. Battersby, J. Bally, J. M. D. Kruijssen, J. D. Henshaw et al., _Star formation rates and efficiencies in the galactic centre_ , _Monthly Notices of the Royal Astronomical Society_ 469 (Apr, 2017) 2263–2285.
  * [10] J. Petrović, P. D. Serpico and G. Zaharijaš, _Galactic Center gamma-ray ”excess” from an active past of the Galactic Centre?_ , _JCAP_ 10 (2014) 052, [1405.7928].
  * [11] D. Gaggero, D. Grasso, A. Marinelli, M. Taoso and A. Urbano, _Diffuse cosmic rays shining in the Galactic center: A novel interpretation of H.E.S.S. and Fermi-LAT $\gamma$-ray data_, _Phys. Rev. Lett._ 119 (2017) 031101, [1702.01124].
  * [12] E. Carlson, T. Linden and S. Profumo, _Improved Cosmic-Ray Injection Models and the Galactic Center Gamma-Ray Excess_ , _Phys. Rev. D_ 94 (2016) 063504, [1603.06584].
  * [13] G. Dobler, D. P. Finkbeiner, I. Cholis, T. Slatyer and N. Weiner, _The Fermi Haze: A Gamma-ray Counterpart to the Microwave Haze_ , _ApJ_ 717 (July, 2010) 825–842, [0910.4583].
  * [14] Fermi-LAT collaboration, M. Ackermann et al., _The Fermi Galactic Center GeV Excess and Implications for Dark Matter_ , _Astrophys. J._ 840 (2017) 43, [1704.03910].
  * [15] L. Goodenough and D. Hooper, _Possible Evidence For Dark Matter Annihilation In The Inner Milky Way From The Fermi Gamma Ray Space Telescope_ , _arXiv e-prints_ (Oct., 2009) arXiv:0910.2998, [0910.2998].
  * [16] Fermi-LAT collaboration, V. Vitale and A. Morselli, _Indirect Search for Dark Matter from the center of the Milky Way with the Fermi-Large Area Telescope_ , in _2009 Fermi Symposium_ , 12, 2009.  0912.3828.
  * [17] D. Hooper and L. Goodenough, _Dark Matter Annihilation in The Galactic Center As Seen by the Fermi Gamma Ray Space Telescope_ , _Phys.Lett._ B697 (2011) 412–428, [1010.2752].
  * [18] K. N. Abazajian, _The Consistency of Fermi-LAT Observations of the Galactic Center with a Millisecond Pulsar Population in the Central Stellar Cluster_ , _JCAP_ 03 (2011) 010, [1011.4275].
  * [19] D. Hooper and T. Linden, _On The Origin Of The Gamma Rays From The Galactic Center_ , _Phys. Rev. D_ 84 (2011) 123005, [1110.0006].
  * [20] K. N. Abazajian and M. 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 (2012) 083511, [1207.6047].
  * [21] D. Hooper, I. Cholis, T. Linden, J. Siegal-Gaskins and T. Slatyer, _Pulsars Cannot Account for the Inner Galaxy’s GeV Excess_ , _Phys. Rev. D_ 88 (2013) 083009, [1305.0830].
  * [22] C. Gordon and O. Macias, _Dark Matter and Pulsar Model Constraints from Galactic Center Fermi-LAT Gamma Ray Observations_ , _Phys. Rev. D_ 88 (2013) 083521, [1306.5725].
  * [23] O. Macias and C. Gordon, _The Contribution of Cosmic Rays Interacting With Molecular Clouds to the Galactic Center Gamma-Ray Excess_ , _Phys.Rev._ D89 (2014) 063515, [1312.6671].
  * [24] T. Daylan, D. P. Finkbeiner, D. Hooper, T. Linden, S. K. N. Portillo et al., _The Characterization of the Gamma-Ray Signal from the Central Milky Way: A Compelling Case for Annihilating Dark Matter_ , 1402.6703.
  * [25] K. N. Abazajian, N. Canac, S. Horiuchi and M. Kaplinghat, _Astrophysical and Dark Matter Interpretations of Extended Gamma Ray Emission from the Galactic Center_ , 1402.4090.
  * [26] B. Zhou, Y.-F. Liang, X. Huang, X. Li, Y.-Z. Fan, L. Feng et al., _GeV excess in the Milky Way: The role of diffuse galactic gamma-ray emission templates_ , _Phys. Rev. D_ 91 (2015) 123010, [1406.6948].
  * [27] X. Huang, T. Enßlin and M. Selig, _Galactic dark matter search via phenomenological astrophysics modeling_ , _JCAP_ 04 (2016) 030, [1511.02621].
  * [28] Fermi-LAT collaboration, M. Ajello et al., _Fermi-LAT Observations of High-Energy $\gamma$-Ray Emission Toward the Galactic Center_, _Astrophys. J._ 819 (2016) 44, [1511.02938].
  * [29] F. Calore, I. Cholis, C. McCabe and C. Weniger, _A Tale of Tails: Dark Matter Interpretations of the Fermi GeV Excess in Light of Background Model Systematics_ , _Phys. Rev. D_ 91 (2015) 063003, [1411.4647].
  * [30] F. Calore, I. Cholis and C. Weniger, _Background Model Systematics for the Fermi GeV Excess_ , _JCAP_ 03 (2015) 038, [1409.0042].
  * [31] M. Di Mauro, _The characteristics of the Galactic center excess measured with 11 years of Fermi-LAT data_ , 2101.04694.
  * [32] I. Cholis, Y.-M. Zhong, S. D. McDermott and J. P. Surdutovich, _Return of the templates: Revisiting the Galactic Center excess with multimessenger observations_ , _Phys. Rev. D_ 105 (2022) 103023, [2112.09706].
  * [33] S. D. McDermott, Y.-M. Zhong and I. Cholis, _A Phantom Menace: On the Morphology of the Galactic Center Excess_ , 2209.00006.
  * [34] R. Bartels, S. Krishnamurthy and C. Weniger, _Strong support for the millisecond pulsar origin of the Galactic center GeV excess_ , _Phys. Rev. Lett._ 116 (2016) 051102, [1506.05104].
  * [35] S. K. Lee, M. Lisanti, B. R. Safdi, T. R. Slatyer and W. Xue, _Evidence for Unresolved $\gamma$-Ray Point Sources in the Inner Galaxy_, _Phys. Rev. Lett._ 116 (2016) 051103, [1506.05124].
  * [36] H. Ploeg, C. Gordon, R. Crocker and O. Macias, _Consistency Between the Luminosity Function of Resolved Millisecond Pulsars and the Galactic Center Excess_ , _JCAP_ 08 (2017) 015, [1705.00806].
  * [37] C. Eckner et al., _Millisecond pulsar origin of the Galactic center excess and extended gamma-ray emission from Andromeda - a closer look_ , _Astrophys. J._ 862 (2018) 79, [1711.05127].
  * [38] G. Fragione, F. Antonini and O. Y. Gnedin, _Disrupted Globular Clusters and the Gamma-Ray Excess in the Galactic Centre_ , _Mon. Not. Roy. Astron. Soc._ 475 (2018) 5313–5321, [1709.03534].
  * [39] G. Fragione, V. Pavlík and S. Banerjee, _Neutron stars and millisecond pulsars in star clusters: implications for the diffuse $\gamma$-radiation from the Galactic Centre_, _Mon. Not. Roy. Astron. Soc._ 480 (2018) 4955–4962, [1804.04856].
  * [40] P. L. Gonthier, A. K. Harding, E. C. Ferrara, S. E. Frederick, V. E. Mohr and Y.-M. Koh, _Population syntheses of millisecond pulsars from the Galactic Disk and Bulge_ , _Astrophys. J._ 863 (2018) 199, [1806.11215].
  * [41] H. Ploeg, C. Gordon, R. Crocker and O. Macias, _Comparing the Galactic Bulge and Galactic Disk Millisecond Pulsars_ , _JCAP_ 12 (2020) 035, [2008.10821].
  * [42] W. de Boer, L. Bosse, I. Gebauer, A. Neumann and P. L. Biermann, _Molecular clouds as origin of the Fermi gamma-ray GeV excess_ , _Phys. Rev. D_ 96 (2017) 043012, [1707.08653].
  * [43] E. Carlson and S. Profumo, _Cosmic Ray Protons in the Inner Galaxy and the Galactic Center Gamma-Ray Excess_ , _Phys. Rev. D_ 90 (2014) 023015, [1405.7685].
  * [44] D. Gaggero, M. Taoso, A. Urbano, M. Valli and P. Ullio, _Towards a realistic astrophysical interpretation of the gamma-ray Galactic center excess_ , _JCAP_ 12 (2015) 056, [1507.06129].
  * [45] L. Goodenough and D. Hooper, _Possible Evidence For Dark Matter Annihilation In The Inner Milky Way From The Fermi Gamma Ray Space Telescope_ , 0910.2998.
  * [46] J. T. Dinsmore and T. R. Slatyer, _Luminosity functions consistent with a pulsar-dominated Galactic Center excess_ , _JCAP_ 06 (2022) 025, [2112.09699].
  * [47] M. Pohl, O. Macias, P. Coleman and C. Gordon, _Assessing the Impact of Hydrogen Absorption on the Characteristics of the Galactic Center Excess_ , _Astrophys. J._ 929 (2022) 136, [2203.11626].
  * [48] O. Macias, C. Gordon, R. M. Crocker, B. Coleman, D. Paterson, S. Horiuchi et al., _Galactic bulge preferred over dark matter for the Galactic centre gamma-ray excess_ , _Nature Astron._ 2 (2018) 387–392, [1611.06644].
  * [49] R. Bartels, E. Storm, C. Weniger and F. Calore, _The Fermi-LAT GeV excess as a tracer of stellar mass in the Galactic bulge_ , _Nature Astron._ 2 (2018) 819–828, [1711.04778].
  * [50] O. Macias, S. Horiuchi, M. Kaplinghat, C. Gordon, R. M. Crocker and D. M. Nataf, _Strong Evidence that the Galactic Bulge is Shining in Gamma Rays_ , _JCAP_ 09 (2019) 042, [1901.03822].
  * [51] K. N. Abazajian, S. Horiuchi, M. Kaplinghat, R. E. Keeley and O. Macias, _Strong constraints on thermal relic dark matter from Fermi-LAT observations of the Galactic Center_ , _Phys. Rev. D_ 102 (2020) 043012, [2003.10416].
  * [52] F. Calore, F. Donato and S. Manconi, _Dissecting the Inner Galaxy with $\gamma$-Ray Pixel Count Statistics_, _Phys. Rev. Lett._ 127 (2021) 161102, [2102.12497].
  * [53] L. Cao, S. Mao, D. Nataf, N. J. Rattenbury and A. Gould, _A New Photometric Model of the Galactic Bar using Red Clump Giants_ , _Mon. Not. Roy. Astron. Soc._ 434 (2013) 595–605, [1303.6430].
  * [54] M. Portail, O. Gerhard, C. Wegg and M. Ness, _Dynamical modelling of the galactic bulge and bar: the Milky Way’s pattern speed, stellar and dark matter mass distribution_ , _Mon. Not. Roy. Astron. Soc._ 465 (2017) 1621–1644, [1608.07954].
  * [55] R. Launhardt, R. Zylka and P. G. Mezger, _The nuclear bulge of the galaxy. 3. Large scale physical characteristics of stars and interstellar matter_ , _Astron. Astrophys._ 384 (2002) 112–139, [astro-ph/0201294].
  * [56] M. B. Davies and B. M. S. Hansen, _Neutron star retention and millisecond pulsar production in globular clusters_ , _MNRAS_ 301 (Nov., 1998) 15–24.
  * [57] C. Y. Hui, K. S. Cheng and R. E. Taam, _Dynamical Formation of Millisecond Pulsars in Globular Clusters_ , _Astrophys. J._ 714 (2010) 1149–1154, [1003.4332].
  * [58] N. Mirabal, _Dark matter vs. Pulsars: Catching the impostor_ , _Mon. Not. Roy. Astron. Soc._ 436 (2013) 2461, [1309.3428].
  * [59] R. Bartels and C. Weniger, _Millisecond Pulsars in the Galactic Bulge? An Extended Discussion on the Wavelet Analysis of the Fermi-LAT data_ , _IAU Symp._ 322 (2016) 193–196.
  * [60] Y.-M. Zhong, S. D. McDermott, I. Cholis and P. J. Fox, _Testing the Sensitivity of the Galactic Center Excess to the Point Source Mask_ , _Phys. Rev. Lett._ 124 (2020) 231103, [1911.12369].
  * [61] R. K. Leane and T. R. Slatyer, _Revival of the Dark Matter Hypothesis for the Galactic Center Gamma-Ray Excess_ , _Phys. Rev. Lett._ 123 (2019) 241101, [1904.08430].
  * [62] L. J. Chang, S. Mishra-Sharma, M. Lisanti, M. Buschmann, N. L. Rodd and B. R. Safdi, _Characterizing the nature of the unresolved point sources in the Galactic Center: An assessment of systematic uncertainties_ , _Phys. Rev. D_ 101 (2020) 023014, [1908.10874].
  * [63] M. Buschmann, N. L. Rodd, B. R. Safdi, L. J. Chang, S. Mishra-Sharma, M. Lisanti et al., _Foreground Mismodeling and the Point Source Explanation of the Fermi Galactic Center Excess_ , _Phys. Rev. D_ 102 (2020) 023023, [2002.12373].
  * [64] R. K. Leane and T. R. Slatyer, _Spurious Point Source Signals in the Galactic Center Excess_ , _Phys. Rev. Lett._ 125 (2020) 121105, [2002.12370].
  * [65] R. K. Leane and T. R. Slatyer, _The enigmatic Galactic Center excess: Spurious point sources and signal mismodeling_ , _Phys. Rev. D_ 102 (2020) 063019, [2002.12371].
  * [66] S. Mishra-Sharma and K. Cranmer, _Neural simulation-based inference approach for characterizing the Galactic Center $\gamma$-ray excess_, _Phys. Rev. D_ 105 (2022) 063017, [2110.06931].
  * [67] E. Storm, C. Weniger and F. Calore, _SkyFACT: High-dimensional modeling of gamma-ray emission with adaptive templates and penalized likelihoods_ , _JCAP_ 1708 (2017) 022, [1705.04065].
  * [68] H.-S. Zechlin, A. Cuoco, F. Donato, N. Fornengo and A. Vittino, _Unveiling the Gamma-ray Source Count Distribution Below the Fermi Detection Limit with Photon Statistics_ , _Astrophys. J. Suppl._ 225 (2016) 18, [1512.07190].
  * [69] S. Caron, G. A. Gómez-Vargas, L. Hendriks and R. Ruiz de Austri, _Analyzing $\gamma$-rays of the Galactic Center with Deep Learning_, _JCAP_ 05 (2018) 058, [1708.06706].
  * [70] F. List, N. L. Rodd, G. F. Lewis and I. Bhat, _The GCE in a New Light: Disentangling the $\gamma$-ray Sky with Bayesian Graph Convolutional Neural Networks_, 2006.12504.
  * [71] S. Mishra-Sharma and K. Cranmer, _Semi-parametric $\gamma$-ray modeling with Gaussian processes and variational inference_, in _34th Conference on Neural Information Processing Systems_ , 10, 2020.  2010.10450.
  * [72] F. List, N. L. Rodd and G. F. Lewis, _Extracting the Galactic Center excess’ source-count distribution with neural nets_ , _Phys. Rev. D_ 104 (2021) 123022, [2107.09070].
  * [73] P. Mertsch and A. Vittino, _Bayesian inference of three-dimensional gas maps. I. Galactic CO_ , _A &A_ 655 (Nov., 2021) A64, [2012.15770].
  * [74] P. Mertsch and V. H. M. Phan, _Bayesian inference of three-dimensional gas maps: II. Galactic HI_ , 2202.02341.
  * [75] C. M. Karwin, A. Broughton, S. Murgia, A. Shmakov, M. Tavakoli and P. Baldi, _Improved Modeling of the Discrete Component of the Galactic Interstellar $\gamma$-ray Emission and Implications for the $Fermi$ –LAT Galactic Center Excess_, 2206.02809.
  * [76] A. Shmakov, M. Tavakoli, P. Baldi, C. M. Karwin, A. Broughton and S. Murgia, _Deep learning models of the discrete component of the Galactic interstellar $\gamma$-ray emission_, _Phys. Rev. D_ 107 (2023) 063018, [2206.02819].
  * [77] R. Bartels, T. Edwards and C. Weniger, _Bayesian model comparison and analysis of the Galactic disc population of gamma-ray millisecond pulsars_ , _Mon. Not. Roy. Astron. Soc._ 481 (2018) 3966–3987, [1805.11097].
  * [78] W. B. Atwood, L. Baldini, J. Bregeon, P. Bruel, A. Chekhtman, J. Cohen-Tanugi et al., _New Fermi-LAT Event Reconstruction Reveals More High-energy Gamma Rays from Gamma-Ray Bursts_ , _ApJ_ 774 (Sept., 2013) 76, [1307.3037].
  * [79] P. Bruel, T. H. Burnett, S. W. Digel, G. Johannesson, N. Omodei and M. Wood, _Fermi-LAT improved Pass ~8 event selection_, _arXiv e-prints_ (Oct., 2018) arXiv:1810.11394, [1810.11394].
  * [80] Fermi-LAT collaboration, F. Acero et al., _The First Fermi LAT Supernova Remnant Catalog_ , _Astrophys. J. Suppl._ 224 (2016) 8, [1511.06778].
  * [81] D. R. Lorimer, A. J. Faulkner, A. G. Lyne, R. N. Manchester, M. Kramer, M. A. McLaughlin et al., _The Parkes Multibeam Pulsar Survey - VI. Discovery and timing of 142 pulsars and a Galactic population analysis_ , _MNRAS_ 372 (Oct., 2006) 777–800, [astro-ph/0607640].
  * [82] L. Herold and D. Malyshev, _Hard and bright gamma-ray emission at the base of the Fermi bubbles_ , _Astron. Astrophys._ 625 (2019) A110, [1904.01454].
  * [83] Fermi-LAT collaboration, J. Ballet, T. H. Burnett, S. W. Digel and B. Lott, _Fermi Large Area Telescope Fourth Source Catalog Data Release 2_ , 2005.11208.
  * [84] Fermi-LAT collaboration, F. Acero et al., _Fermi Large Area Telescope Third Source Catalog_ , _Astrophys. J. Suppl._ 218 (2015) 23, [1501.02003].
  * [85] J. F. Navarro, C. S. Frenk and S. D. M. White, _The Structure of cold dark matter halos_ , _Astrophys. J._ 462 (1996) 563–575, [astro-ph/9508025].
  * [86] J. F. Navarro, C. S. Frenk and S. D. M. White, _A Universal density profile from hierarchical clustering_ , _Astrophys. J._ 490 (1997) 493–508, [astro-ph/9611107].
  * [87] J. Zavala and C. S. Frenk, _Dark matter haloes and subhaloes_ , _Galaxies_ 7 (2019) 81, [1907.11775].
  * [88] A. McCann, _A stacked analysis of 115 pulsars observed by the Fermi LAT_ , _Astrophys. J._ 804 (2015) 86, [1412.2422].
  * [89] Fermi-LAT collaboration, M. Ajello et al., _Characterizing the population of pulsars in the inner Galaxy with the Fermi Large Area Telescope_ , 1705.00009.
  * [90] M. Winter, G. Zaharijas, K. Bechtol and J. Vandenbroucke, _Estimating the GeV Emission of Millisecond Pulsars in Dwarf Spheroidal Galaxies_ , _Astrophys. J. Lett._ 832 (2016) L6, [1607.06390].
  * [91] D. Hooper and G. Mohlabeng, _The Gamma-Ray Luminosity Function of Millisecond Pulsars and Implications for the GeV Excess_ , _JCAP_ 03 (2016) 049, [1512.04966].
  * [92] Fermi-LAT collaboration, F. Acero et al., _Fermi Large Area Telescope Third Source Catalog_ , _Astrophys. J. Suppl._ 218 (2015) 23, [1501.02003].
  * [93] B. Limyansky, _The Third Fermi Pulsar Catalog_ , in _AAS/High Energy Astrophysics Division_ , AAS/High Energy Astrophysics Division, p. 109.32, Mar., 2019.
  * [94] T. Bringmann and C. Weniger, _Gamma Ray Signals from Dark Matter: Concepts, Status and Prospects_ , _Phys. Dark Univ._ 1 (2012) 194–217, [1208.5481].
  * [95] B. Lakshminarayanan, A. Pritzel and C. Blundell, _Simple and scalable predictive uncertainty estimation using deep ensembles_ , 2017.
  * [96] S. Ioffe and C. Szegedy, _Batch normalization: Accelerating deep network training by reducing internal covariate shift_ , ICML’15, p. 448–456, JMLR.org, 2015.
  * [97] D. P. Kingma and J. Ba, _Adam: A method for stochastic optimization_ , 2014\.
  * [98] H. Dembinski, P. Ongmongkolkul, C. Deil, D. M. Hurtado, M. Feickert, H. Schreiner et al., _scikit-hep/iminuit: v1.5.2_ , .
  * [99] Fermi-LAT collaboration, 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 (2016) 26, [1602.07246].
  * [100] Fermi-LAT collaboration, 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 (2016) 26, [1602.07246].
  * [101] L. Ruff, R. Vandermeulen, N. Goernitz, L. Deecke, S. A. Siddiqui, A. Binder et al., _Deep one-class classification_ , in _Proceedings of the 35th International Conference on Machine Learning_ (J. Dy and A. Krause, eds.), vol. 80 of _Proceedings of Machine Learning Research_ , pp. 4393–4402, PMLR, 10–15 Jul, 2018.
  * [102] S. Caron, L. Hendriks and R. Verheyen, _Rare and different: Anomaly scores from a combination of likelihood and out-of-distribution models to detect new physics at the LHC_ , _SciPost Physics_ 12 (feb, 2022) .

##  Appendix A Astrophysical parameters

In this section, we present a complementary view of the results of the network
findings in terms of astrophysical components. In Fig. 15 we show the
network’s prediction for parameters not directly related to the GCE excess,
that were not part of Fig. 10 in the main text while retaining the style of
this plot. Note that the two astrophysical parameters of Model 1 are not shown
since they do not fit reasonably into this scheme.

We notice a slight anti-correlation between the evolution per model of the
normalizations for the isotropic and inverse-Compton templates. Moreover, both
scenarios indicate a trend of the 4FGL-DR2 template normalizations to decrease
with increasing model complexity. Lastly, the agreement between scenarios A
and B regarding the astrophysical parameters of our models are best for Model
3 lending further credence to the fact that it represents the best
configuration among the tested cases.

Figure 15: Same as Fig. 10 (upper panel) displaying a comparison of the
network’s prediction for relevant astrophysical parameters as a function of
the model number. Note that the template normalisation parameters associated
with the FBs templates have been converted into a photon flux over the full
energy range of our analysis using the assumed FBs spectral shape (see Sec.
2.2.1). We do not provide corresponding numbers for Model 1 since it does not
exhibit the same structure as the other two model iterations.

##  Appendix B MSP template without flux detection threshold - Scenario B

In this section, we show the results of the network in the case where the MSP
template generation proceeds without a high-flux cut-off, i.e. we generate
MSPs both below and above the LAT detection threshold. This is an approach
dubbed Scenario B in Section 2 and illustrated schematically in Fig. 1.

Figure 16: The color bar marks the number of MSPs in our mock catalogue that
have a flux above the 4FGL-DR2 detection threshold, for a given choice of
$\sigma_{L}$ and total luminosity of the MSP template $F_{MSP}$ (see Section 2
for details). The spatial distribution is fixed to $\gamma=1.28$. With the
dashed vertical line, we plot the best-fit value of $\sigma_{L}$ obtained in
[77], while the gray zone marks $1-3$ $\sigma$ uncertainty range. The
horizontal, dashed yellow line denotes the best-fit MSP cumulative flux with
respect to Model 3B while the height of the yellow box displays its 1$\sigma$
error margin. Analogously, the vertical shaded yellow line refers to the best-
fit value of $\sigma_{L}$ in this model and scenario; the width of the box
shows the 1$\sigma$ uncertainty of this parameter.

To gauge how many sources we add above the 4FGL detection threshold in this
procedure (risking the ’double counting’ issue) we show in Figure 16 the
number of sources above the 4FGL threshold that stay in our MSP template in
this scenario. We see that the number of bright sources is subdominant with
respect to the total number of 4FGL sources in our ROI which is 243 4FGL point
sources, within the angular distance to the GC of $|l,b|\leq 15^{\circ}$ (this
includes 22 pulsars and 153 unidentified objects).

All other procedure steps are identical as in the main text. Results are shown
in Figures 17, 18, 19 and the corresponding residuals of the models selected
by the network in Figs. 22, 23, 24.

Figure 17: The same as Fig. 5 (Model 1), but for Scenario B.

We observe that

  * •

All main results are unchanged in Scenario B. Indeed the network seems to
predict quite reliably the total flux of the MSP template, independently of
the scenario (see Fig. 10).

  * •

As expected, the network has higher sensitivity on the $\sigma_{L}$ parameter,
which is found to be in the $\sigma_{L}~{}\epsilon~{}[0.4-0.8]$ range,
consistent with the findings in [77]. This parameter was instead fully
unconstrained in Scenario A, as expected since in that case, the luminosity
function was discontinuous at the 4FGL detection threshold.

Figure 18: Model 2, Scenario B.  Figure 19: Model 3, Scenario B.

##  Appendix C Spatial and spectral residuals for the examined models and
scenarios

In this section, we list all remaining plots of spatial and spectral residuals
that have not been shown in the main text.

While the structure of the spatial and spectral residuals for Model 1A was a
good indicator and motivation for several model extensions, the corresponding
figures for other model iterations do not change significantly. Moreover,
these residuals are not the result of standard inference methods like a
maximum likelihood fit and should therefore not be interpreted as such. The
inherent purpose of a maximum likelihood analysis is to find those model
parameters that minimize the residual differences between the final best-fit
model and the input data. The machine learning approach employed in this work
rather aims at deriving a relation between the model parameters and the images
it is trained on, i.e. its task is to establish a parameter reconstruction
scheme. The resulting residuals on the image level as shown in, for example,
Fig. 6 are not part of the network’s loss function in Eq. 3.1. Therefore, it
is by no means guaranteed that even a perfectly unbiased and accurate
parameter reconstruction leads to the same result as a maximum likelihood
approach. The presence of a pronounced reality gap between the model and true
gamma-ray data may enhance this effect. In the explicit example of models and
scenarios probed in our work, we find that the spectral residuals in the last
two energy bins tend to become worse with increased model complexity whereas
the first three energy bins – those bins that contain the bulk of the GCE
emission – are fit fairly well. Such a behavior seems odd in the context of a
maximum likelihood fit but it might be reasonable for the optimization measure
created internally by training the neural network. Moreover, the error bars of
the spectral residuals are potentially overestimated because we lack
information on correlations among the parameters predicted by the network.
Retrieving such correlations is a field of active research in the machine-
learning community.

Figure 20: The same as Fig. 6 but in the context of Model 2.

Figure 21: The same as Fig. 6 but in the context of Model 3.

Figure 22: The same as Fig. 6 (Model 1), but for Scenario B.

Figure 23: Spatial and spectral residuals for Model 2B.

Figure 24: Spatial and spectral residuals for Model 3B.

##  Appendix D Consistency check with an alternative diffuse background model

In this appendix, we list the results of our consistency checks invoking the
Pass 8 Fermi diffuse background model used in the context of the 4FGL source
catalog supplementing the discussion in the main text (see Sec. 6.2).

Figs. 25 to 27 illustrate the bias of the network’s predictions induced by
applying a validation data set with a diffuse model the network has not been
initially trained on. For the comparison, we only select those model
parameters that can be reasonably confronted with each other, i.e. parameters
that are shared between the Fermi diffuse validation data set and the original
validation data set of the respective model and scenario.

The qualitative interpretation of these results is identical to the one
discussed in the main text using the example of Model 3B. Since the Fermi
diffuse model incorporates an additional (smoothed) contribution from the
Galactic center region that comprises the GCE and parts of the FBs, we obtain
biased network predictions. Each individual network tries to assign the excess
gamma-ray emission to several components that may absorb it according to the
training data it has seen. What components are affected strongly depends on
the respective model and scenario. If the model contains templates for the
FBs, the predicted normalization factor is always biased towards higher values
in full agreement with the expectations. Another common feature is the bias of
the gNFW inner slope parameter $\gamma$ towards larger values since an
additional GCE component is present in the validation data. However, the
excess emission of the Fermi diffuse model is not consistently attributed to
specific templates among models and scenarios. This again illustrates that
overcoming the reality gap is essential to _(i)_ robustly establish the
properties of the GCE and _(ii)_ identify the leading process that generates
it in the first place.

Figure 25: Same as Fig. 12 for Model 1 (Top: A); (Bottom: B).

Figure 26: Same as Fig. 12 for Model 2 (Top: A); (Bottom: B).  Figure 27: Same
as Fig. 12 for Model 3A.