#  Robustness of the Galactic Center Excess Morphology Against Masking

Yi-Ming Zhong  Department of Physics, City University of Hong Kong, Kowloon,
Hong Kong SAR, China  Kavli Institute for Cosmological Physics, University of
Chicago, Chicago, IL 60637, USA<EMAIL_ADDRESS>Ilias Cholis
Department of Physics, Oakland University, Rochester, Michigan, 48309, USA
<EMAIL_ADDRESS>
###### Abstract

The Galactic Center Excess (GCE) remains an enduring mystery, with leading
explanations being annihilating dark matter or an unresolved population of
millisecond pulsars. Analyzing the morphology of the GCE provides critical
clues to identify its exact origin. We investigate the robustness of the
inferred GCE morphology against the effects of masking, an important step in
the analysis where the gamma-ray emission from point sources and the galactic
disk are excluded. Using different masks constructed from Fermi point source
catalogs and a wavelet method, we find that the GCE morphology, particularly
its ellipticity and cuspiness, is relatively independent of the choice of mask
for energies above $2-3$ GeV. The GCE morphology systematically favors an
approximately spherical shape, as expected for dark matter annihilation.
Compared to various stellar bulge profiles, a spherical dark matter
annihilation profile better fits the data across different masks and galactic
diffuse emission backgrounds, except for the stellar bulge profile from
Coleman et al. (2020), which provides a similar fit to the data. Modeling the
GCE with two components, one from dark matter annihilation and one tracing the
Coleman Bulge, we find this two-component model outperforms any single
component or combinations of dark matter annihilation and other stellar bulge
profiles. Uncertainty remains about the exact fraction contributed by each
component across different background models and masks. However, when the
Coleman Bulge dominates, its corresponding spectrum lacks characteristics
typically associated with millisecond pulsars, suggesting that it mostly
models the emission from other sources instead of the GCE that is still
present and spherically symmetric.

##  I Introduction

The Galactic Center Excess (GCE) has remained an unexplained component of
gamma-ray emission for a decade and a half since it was first identified in
the Fermi-LAT data by Ref. Goodenough and Hooper (2009). The GCE peaks in its
energy output around a GeV. It appears to be a diffuse and extended component
of gamma rays that can not be explained by steady-state galactic diffuse
emission, originating from the interaction with the interstellar medium of
cosmic rays, typically produced by astrophysical sources such as pulsars or
supernova remnants Hooper and Goodenough (2011); Vitale and Morselli (2009);
Abazajian (2011); Hooper and Linden (2011); Hooper and Slatyer (2013); Gordon
and Macias (2013); Abazajian  _et al._ (2014); Daylan  _et al._ (2016); Calore
_et al._ (2015a); Zhou  _et al._ (2015); Ajello  _et al._ (2016); Linden  _et
al._ (2016); Cholis  _et al._ (2022). The exact origin of the GCE is still a
topic of ongoing research debate, with the two most intriguing explanations
being dark matter annihilation Gordon and Macias (2013); Daylan  _et al._
(2016); Calore  _et al._ (2015b); Agrawal  _et al._ (2015); Berlin  _et al._
(2015); Karwin  _et al._ (2017); Ackermann  _et al._ (2017); Leane and Slatyer
(2019) and an unresolved population of millisecond pulsars (MSPs) Abazajian
and Kaplinghat (2012); Hooper  _et al._ (2013); Lee  _et al._ (2016); Petrović
_et al._ (2015); Cholis  _et al._ (2015a); Bartels  _et al._ (2018); Buschmann
_et al._ (2020); Gautam  _et al._ (2022); while also a sequence of cosmic-ray
burst events from the center of the galaxy provides a third alternative
Petrović  _et al._ (2014); Carlson and Profumo (2014); Cholis  _et al._
(2015b). If the dark matter explanation is true, the GCE will not only serve
as the first evidence that dark matter interacts with ordinary matter beyond
its gravitational pull, but it will also shed light on its production
mechanism. If the MSPs explanation is true, it could reveal a new population
of MSPs with luminosity function that is different from the known pulsars in
globular clusters or the galactic disk Cholis  _et al._ (2015a); Zhong  _et
al._ (2020). Finally, if the GCE comes from a sequence of burst events, it
would help us understand the conditions and recent history of the inner region
of the Milky Way.

Various efforts have been pursued to identify the origin of the GCE, which
include analyzing the details of its energy spectrum, its photon statistics,
and its morphology Daylan  _et al._ (2016); Calore  _et al._ (2015a); Bartels
_et al._ (2016); Ajello  _et al._ (2016); Macias  _et al._ (2018, 2019); Zhong
_et al._ (2020); Buschmann  _et al._ (2020); Di Mauro (2021); Cholis  _et al._
(2022); Pohl  _et al._ (2022); McDermott  _et al._ (2022). A common step
adopted in many GCE analyses is masking certain areas of the galactic center
(GC). Masks are applied to exclude pixels of the sky where conventional
sources of gamma-ray emission are expected to be so bright that even the
uncertainty on their expected flux dominates the possible signal. Those
include the galactic disk and the known point (or extended) sources that could
interfere with the GCE signal. Various methods exist for creating such masks.
One standard approach involves using the Fermi-LAT point source catalog to
exclude regions around the cataloged sources, combined with band-shaped masks
along $b=0^{\circ}$ to exclude the disk.

The choice of the mask and method to apply them to the data can potentially
affect the results of the GCE analysis. For example, earlier works have
demonstrated such masking effects on the Non-Poissonian Template Fitting Lee
_et al._ (2016), where switching from masks created using the Fermi-LAT 3FGL
catalog to the 4FGL catalog diminishes the statistical power of the method on
the GCE photon statistics Buschmann  _et al._ (2020). Refs. Zhong  _et al._
(2020); Cholis  _et al._ (2022) investigated how the GCE characteristics are
affected using alternative masks for the known Fermi-LAT point sources. They
found its energy spectrum and morphology robust to the different updates of
the Fermi-LAT catalog. Given the significance of the GCE morphology in its
interpretation, we perform an extensive analysis of the robustness of the GCE
morphology to alternative masking procedures in this work. We do that by
combining with the vast publicly available library of galactic diffuse
emission models produced in Ref. Cholis  _et al._ (2022) to account for the
astrophysical modeling uncertainties of the Milky Way conditions. Such an
analysis is especially timely given the increasing number of point sources
released in the latest Fermi-LAT catalogs Abdollahi  _et al._ (2020); Ballet
_et al._ (2020); Abdollahi  _et al._ (2022); Ballet  _et al._ (2023).

We begin in Sec. II by introducing the Fermi data, the templates, the masks
tested, and the statistical procedure used in our analysis. Sec. III presents
our injection tests for the template fitting method in a controlled
environment. Sec. IV shows the results of the GCE morphology under different
masks. We find that the GCE morphology remains largely spherical and resembles
that of a dark matter annihilation profile rather than that of a population of
dim gamma-ray point sources tracing known stellar populations. An additional
population of faint gamma-ray point sources tracing the stellar distribution
known as the “Coleman Bulge” Coleman  _et al._ (2020) may exist. However, the
spectrum associated with this component does not possess the characteristics
expected from MSPs. We conclude in Sec. V.

##  II Setup for the template fitting

Here, we describe the Fermi data and the templates used for the diffuse
emission and the GCE. We also describe the masks applied and the statistical
procedure used in our template fitting analysis. Except for the mask
construction, our analysis framework largely follows that of Refs. Cholis  _et
al._ (2022); McDermott  _et al._ (2022).

###  II.1 Data

We use Fermi Pass 8 data, version P8R3, recorded from August 4, 2008, to April
14, 2021, corresponding to weeks $9-670$ of Fermi-LAT
observations.111https://fermi.gsfc.nasa.gov/ssc/data/access/ We use Fermi
ScienceTools P8v27h5b5c8 for selection cuts and to calculate the relevant
exposure-cube files and exposure
maps,222https://fermi.gsfc.nasa.gov/ssc/data/analysis/ which allow us to pass
from fluxes to expected counts. The exposure is calculated for each pixel
using the Fermi ScienceTools P8v27h5b5c8.

We keep only FRONT-converted CLEAN data. In addition, we set the following
filters: zmax = $100^{\circ}$, DATA$\\_$QUAL==1, LAT$\\_$CONFIG==1, and
ABS(ROCK$\\_$ANGLE) < 52. Our data maps are centered at the galactic center
and cover a square window of $40^{\circ}$ per side in galactic coordinates in
Cartesian pixels of size $0.1^{\circ}\times 0.1^{\circ}$. Unlike a HEALPix
pixelization, our pixels do not have equal area, but we account for this in
our fits.

| $E_{\textrm{min}}-E_{\textrm{max}}{\rm\,[GeV]}$ | $\theta_{s}[^{\circ}]$ | $\theta_{l}[^{\circ}]$
---|---|---|---
0 | $0.275-0.357$ | 1.125 | 3.75
1 | $0.357-0.464$ | 0.975 | 3.25
2 | $0.464-0.603$ | 0.788 | 2.63
3 | $0.603-0.784$ | 0.600 | 2.00
4 | $0.784-1.02$ | 0.450 | 1.50
5 | $1.02-1.32$ | 0.375 | 1.25
6 | $1.32-1.72$ | 0.300 | 1.00
7 | $1.72-2.24$ | 0.225 | 0.750
8 | $2.24-2.91$ | 0.188 | 0.625
9 | $2.91-3.78$ | 0.162 | 0.540
10 | $3.78-4.91$ | 0.125 | 0.417
11 | $4.91-10.8$ | 0.100 | 0.333
12 | $10.8-23.7$ | 0.060 | 0.200
13 | $23.7-51.9$ | 0.053 | 0.175
Table 1: The energy bins and the energy-dependence of the fiducial radii
$\theta_{s}$ ($\theta_{l}$) used to mask known point sources with $\rm{TS}<49$
($\rm{TS}\geq 49$).

We bin the gamma-ray data in 14 energy bins spanning energies from 0.275 GeV
to 51.9 GeV, given in the first column of Tab. 1. The first eleven energy bins
have a constant log width; because the gamma-ray flux drops at higher energy,
the final three energy bins are wider, so each bin has a roughly similar
statistical impact in our fits.

###  II.2 Galactic Diffuse Emission and GCE Templates

We use the 80 astrophysical models that predict an equal number of sets of
high-resolution energy-dependent templates developed in Ref. Cholis  _et al._
(2022) for the galactic diffuse emission background.333We use the term
“background” to describe the combination of the galactic diffuse emission, the
isotropic and mostly extragalactic emission, and the emission from the Fermi
bubbles. If the GCE is the “signal”, the term “background” is used in this
context as the non-signal emission. For our region of interest (ROI),
depending on the latitude and energy, about 50–80$\%$ of the total emission
and 60–90$\%$ of the galactic diffuse emission is truly foreground emission,
i.e., emission generated between our location and where the inner galaxy
region begins (inner 3 kpc from the center of the Milky Way). See Ref. Calore
_et al._ (2015a) for a more detailed analysis of this. Each set includes
templates for $\pi^{0}$, Inverse Compton Scattering (ICS), and Bremsstrahlung,
calculated for each energy bin used in this analysis. In addition, we include
an energy-independent template for the isotropic background and a template for
the Fermi bubbles. In all cases, we projected the templates onto a
$0.1^{\circ}\times 0.1^{\circ}$ Cartesian pixel grid near the galactic center
that covers a square window of $40^{\circ}\times 40^{\circ}$. All templates
are subsequently smoothed by the energy-dependent Fermi-LAT point spread
function psf (2019).

Figure 1: The masks used in this work. Here, we show the masked regions (dark
purple) for the energy bin 5 ($1.02\,\text{GeV}<E<1.32\,\text{GeV}$) only.

For the GCE, we consider contracted, squared, and integrated along the line-
of-sight Navarro-Freak-White (NFW) profiles of the form,

$\rho(r)=\frac{\rho_{0}}{(r/r_{c})^{\gamma}(1+r/r_{c})^{3-\gamma}},$ (1)

where the cuspiness parameter $\gamma$ represents the inner slope of the
density profile, $r_{c}=20\,\text{kpc}$ is the scale radius and $\rho_{0}$ is
the scale density that is fixed by
$\rho(r=8.5\,\text{kpc})=0.4\,\text{GeV}/\text{cm}^{3}$. We further deform the
NFW profile by introducing an ellipticity (elongation) parameter $\epsilon$,
such that the opening angle from the GC, $\psi$, is related to the galactic
coordinates ($b$, $\ell$) by

$\cos(\psi)=\cos(b)\cos(\ell/\epsilon).$ (2)

A spherically symmetric NFW template has $\epsilon=1$ while an oblate
(prolate) NFW template has $\epsilon>1(<1)$. We construct the GCE template
with a combination of $\gamma$ and $\epsilon$ values. All the GCE templates
are projected to the same Cartesian grid as the galactic diffuse emission and
isotropic templates and smoothed by the same point spread functions.

###  II.3 Masks

We construct a variety of masks for the template fitting analysis.444The masks
are available at https://github.com/ymzhong/gce_mask. Most of the masks are a
combination of the point source mask and the disk mask.

_Point source masks_ : We take the location information of the point sources
from the 4FGL-DR1 Abdollahi  _et al._ (2020), DR2 Ballet  _et al._ (2020), and
DR3 Abdollahi  _et al._ (2022) catalogs. The point sources are classified into
those with a test statistic (TS) < 49 and those with TS$\geq 49$. We masked
the two classes of point sources with disk-shaped masks that have energy-
dependent radii $\theta_{s}$ and $\theta_{l}$ respectively (the same radii
used in Cholis  _et al._ (2022)), as listed in Tab. 1. We call those point
sources masks the “Standard 4FGLDRX” with X$=1,2,3$.

On top of the “Standard 4FGLDRX” masks, we introduce two variations on the
mask radii. For the “Large 4FGLDRX” masks, we increase the fiducial radii by a
factor of 1.5 for the point sources with TS$<49$ while keeping radii for the
sources with TS$\geq 49$ unchanged. For the “Small 4FGLDRX” masks, we decrease
the fiducial radii by a factor of 0.5 for the point sources with TS$<49$. The
radii for the sources with TS$\geq 49$ are again unchanged.

_Disk masks_ : We set up three band-shaped masks for the galactic disk. For
the standard disk mask (“L20”), we mask the region $|\ell|<20^{\circ}$ and
$|b|<2^{\circ}$, i.e., the entire disk region of ROI. For the “L8” (“L5”) disk
mask, we mask the region $|\ell|<8^{\circ}$ ($|\ell|<5^{\circ}$) and
$|b|<2^{\circ}$.

_Combined regular masks_ : Combing the point source masks and the disk masks
described above, we set up the following “regular” masks for our analysis:
“Standard 4FGLDRX + L20”, “Large 4FGLDRX + L20”, “Small 4FGLDRX + L20”,
“Standard 4FGLDRX + L8”, and “Standard 4FGLDRX + L5”. For “Standard 4FGLDRX +
L20”, we create masks for X=$1,2,3$. For other variations, the masks are only
built for X=$2,3$. We show some of those combined masks in the first seven
panels of Fig. 1.

_Wavelet-based masks_ : The wavelet method Bartels  _et al._ (2016); Zhong
_et al._ (2020) is a powerful way to identify potential point sources. Here,
we utilize the wavelet peaks found in Zhong  _et al._ (2020). We select peaks
with the statistic $\mathcal{S}>2$, $3$, and $4$ (as defined in Refs. Bartels
_et al._ (2016); Zhong  _et al._ (2020)), respectively, masking each peak with
a disk-shaped mask that has energy-dependent radii identical to those used for
the point sources with TS$<49$. We then combine these wavelet peak masks with
the “Standard 4FGLDRX” point source masks and leave out the disk masks. The
resulting masks are “Standard 4FGLDRX + Wavelet SY” with X=$2,3$ and
Y$=2,3,4$, where the Y values represent the thresholds of $\mathcal{S}$. We
show the masks in the last three panels of Fig. 1.

Figure 2: Masked pixel fraction for various masks used in the analysis. The
vertical bands indicate the setup of 14 energy bins.

For each mask, we compute the fraction of the masked pixels in each energy bin
to the total number of pixels in the inner $40^{\circ}\times 40^{\circ}$
galactic center region ($1.6\times 10^{5}$ pixels). The results are shown in
Fig. 2. Overall, the masked area decreases as the energy bin’s energy
increases. This reflects that the mask radii of the point sources decrease,
following the Fermi point spread function as the energy increases. Looking
into more detail, the top panel of Fig. 2 highlights the differences between
masks created by different point source catalogs as well as by different radii
criteria. As the number of point sources in the catalog increases, the masked
pixel fraction significantly increases. The difference is significant for
$E<2.24\,\text{GeV}$ and gets diminished for higher energy bins as the number
of masked pixels by the disk dominates over those by the point sources.
Increasing or decreasing the radii criteria introduces a significant variation
in the masked pixel fraction, which is the most significant for the lower
energy bins. The middle panel of Fig. 2 highlights the difference in how we
mask the galactic disk. Here, the difference between masking the entire
galactic disk or parts of it is most significant at the higher energy bins
($E>0.603\,\text{GeV}$). For lower energy bins, the two sets of masks yield a
similar masked pixel fraction, given that the point sources mask dominates the
lower latitudes. Combining the information from the two panels, we can
conclude that energy bins 0–3 are mostly affected by point sources masks,
energy bins 8–13 are mostly affected by the disk mask, while energy bins 4–7
are influenced by both types. The lower panel of Fig. 2 highlights the
difference between different wavelet-based masks and those with a band mask
for the disk. The wavelet-based masks, in general, mask a greater fraction of
the GC region than the “Standard 4FGLDRX + L20” masks at lower energy bins.
However, as energy increases and the mask radii shrink, the wavelet-based
masks are not blocked by the disk mask and have a sub 10% masked pixel
fractions for $E>3\,\text{GeV}$. Also, at energies around 1 GeV, typically we
only get 30% of the inner $40^{\circ}\times 40^{\circ}$ region masked.

###  II.4 Statistical Procedure

We construct the log-likelihood $\ln\mathcal{L}_{j}$ for a given energy bin
$j$ as,

$\displaystyle-2\ln\mathcal{L}_{j}={}$
$\displaystyle\Big{\\{}2\sum_{p}\left[\mathcal{M}_{j,p}\mathcal{C}_{j,p}+\ln\left[(\mathcal{M}_{j,p}\mathcal{D}_{j,p})!\right]\right.$
$\displaystyle\left.-(\mathcal{M}_{j,p}\mathcal{D}_{j,p})\ln(\mathcal{M}_{j,p}\mathcal{C}_{j,p})\right]\Big{\\}}$
$\displaystyle+\chi_{\textrm{Bubbles},j}^{2}+\chi_{\textrm{Iso},j}^{2}.$ (3)

The index $p$ runs over the pixels of the ROI. The data $\mathcal{D}_{j,p}$
and the mask $\mathcal{M}_{j,p}$ are described in detail in Sec. II.1 and Sec.
II.3, respectively. The expected counts
$\mathcal{C}_{j,p}=\mathcal{E}_{j,p}\sum_{i}c^{i}_{j}\Phi^{i}_{j,p}$ are
obtained from summing over all the diffuse emission component templates
$\Phi^{i}_{j,p}$ ($i$ runs over different components) with independent
normalizations $c_{j}^{i}$. We run fits without the GCE and with one or two
templates to describe it. In all cases, our fits assume independent
normalizations for each template. After linearly adding the templates, we
multiply their sum by the exposure $\mathcal{E}_{j,p}$. The “external
$\chi^{2}$” functions $\chi^{2}_{\textrm{Bubbles},j}$ and
$\chi^{2}_{\textrm{Iso},j}$ are constraints that act as penalties when the
Bubbles and isotropic normalizations deviate too much from their spectra
measured at high latitudes, see Ref. Ackermann  _et al._ (2014, 2015).

We use emcee Foreman-Mackey  _et al._ (2013), a Markov chain Monte Carlo
(MCMC) program, to sample the coefficient parameters $c^{i}_{j}$, which also
conveniently yields the parameters that minimize the negative log-likelihood,
Eq. 3. We use EnsembleSampler with 100 walkers and 1000 steps to assess
convergence. We then discard the first 300 steps of the resulting chain from
emcee and draw posterior probabilities for the parameters with ChainConsumer
Hinton (2016). The total negative log-likelihood is a sum over the value of
all the energy bins,

$-2\ln\mathcal{L}\equiv-2\sum_{j}\ln\mathcal{L}_{j}.$ (4)

Figure 3: Top two rows: The inferred energy spectra of the GCE from the mock
Fermi data for the four injection tests, utilizing all 80 background galactic
diffuse emission models of Ref Cholis  _et al._ (2022) and the “Standard
4FGLDR3 + L20” mask. Further details are provided in the text. Bottom panel:
the inferred GCE spectrum from the actual Fermi data for comparison. In the
fit, the GCE flux is permitted to be negative, although this occurs in only a
small subset of models and exclusively at energies below 0.7 GeV (indicated by
the termination of lines). For the five background models out of the 80 that
yield the best overall fit to the gamma-ray data, the best-fit GCE
normalization and its corresponding $2\sigma$ uncertainty range are
represented by purple lines and magenta bands, respectively. The green lines
give the GCE normalization from the five galactic diffuse emission models that
provide the worst fit to the gamma-ray data out of the 80 considered.

##  III Injection test

Injecting a GCE component into the Fermi data has been utilized to assess the
sensitivity and robustness of a method in determining the characteristics of
the GCE Leane and Slatyer (2019). In this study, we conduct an injection test
to evaluate how well the template fitting can identify the ellipticity
parameter $\epsilon$ of the GCE. We investigate two GCE signals: a spherical
NFW profile with parameters $\gamma=1.2,\epsilon=1.0$ and an oblate NFW
profile with $\gamma=1.2,\epsilon=1.4$. The GCE signals are injected into the
Fermi data with a flux twice the fiducial value. Note that the template fit of
the actual Fermi data indicates that the observed flux of the GCE is around
the fiducial values for both signals. Utilizing these mock Fermi data sets, we
conduct four separate template-fitting analyses as outlined in Tab. 2. In
these tests, scenarios (1) and (4) match the injected data with the correct
GCE template, whereas scenarios (2) and (3) involve a mismatch between the
injected data and the template.

| Mock data | Fitting GCE Template
---|---|---
(1) |  Fermi+2$\times$NFW($\epsilon=1.0$) | 80 GDEs, NFW($\epsilon=1.0$)
(2) |  Fermi+2$\times$NFW($\epsilon=1.0$) | 80 GDEs, NFW($\epsilon=1.4$)
(3) |  Fermi+2$\times$NFW($\epsilon=1.4$) | 80 GDEs, NFW($\epsilon=1.0$)
(4) |  Fermi+2$\times$NFW($\epsilon=1.4$) | 80 GDEs, NFW($\epsilon=1.4$)
Table 2: The mock Fermi data and templates used for the injection tests. For
the GCE templates, we use a spherical ($\epsilon=1.0$) or oblate
($\epsilon=1.4$) NFW profile with the cuspiness parameter of $\gamma=1.2$.
“GDE” stands for galactic diffuse emission model.

Figure 4: Comparison of the difference in the total negative log-likelihood,
$-2\Delta\ln\mathcal{L}$, for injection test scenario (1) vs. scenario (2)
(left), scenario (3) vs. scenario (4) (right). Colored lines represent some of
the galactic diffuse emission models that give the best fit to the Fermi data.

The resulting energy spectra of the GCE are shown in the top two rows of Fig.
3, where we apply the “Standard 4FGLDR3 + L20” mask. The magenta bands
indicate the $2\sigma$ fit ranges for the five best-fit background galactic
diffuse emission models. For comparison, the lower panel of Fig. 3 shows the
GCE spectrum obtained from fitting the actual Fermi data, applying the same
mask. Upon initial examination, all the injection tests, regardless of whether
the GCE template matches or mismatches, show fitted spectra that are
consistent with the injected GCE spectra (shown as black lines) up to 10 GeV
(except a small discrepancy in scenario (2)). Notably, the fluxes in these
tests are approximately three times higher than those obtained from the actual
Fermi data. That should be the case since the Fermi data already contains the
GCE component before injecting the additional GCE mock component into it.

We further look into the specifics by comparing the total negative log-
likelihoods in scenario (1) vs. scenario (2) and scenario (3) vs. scenario
(4). The results are shown in Fig. 4 in terms of the difference of total
negative log-likelihoods, $-2\Delta\ln\mathcal{L}$, relative to the best-fit
models that include both the galactic diffuse emission backgrounds and the
GCE. Our analysis shows that, in general, the cases where the GCE model
matches the injected data (scenarios (1) and (4)) result in a better fit
compared to their mismatched counterparts (scenarios (2) and (3)) across the
majority of the 80 galactic diffuse emission background models fitted. The
difference in $-2\Delta\ln\mathcal{L}$ ranges from approximately $50$ to $100$
for the five best-fit background models when comparing matched with mismatched
scenarios. Only a few background models that yield poorer fits show a
preference for the mismatched scenarios over the matched ones. Our results
suggest that the template fitting demonstrates robustness against
misinterpreting the ellipticity of the GCE and quantifies the relevant
statistical penalty.

##  IV The morphology of the GCE

In the subsequent sections, we provide a comprehensive list of tests focusing
on the morphology of the GCE. Our investigation covers the properties of the
GCE morphology, i.e., its cuspiness and ellipticity, across different energies
and under different masks as introduced in Sec. II.3 (see Fig. 1).
Additionally, we evaluate and compare various alternative models for the GCE
morphology proposed in the literature. These models encompass scenarios such
as a dark matter annihilation profile, a series of stellar bulge profiles, and
a combination of the dark matter annihilation profile and a stellar bulge.
Each model is examined under different masks to determine their respective
goodness of fit.

Figure 5: The evolution of the fit GCE morphology with energy using the
“Standard 4FGLDR3 + L20” mask. Left: how the inner slope (cuspiness) parameter
$\gamma$ changes with energy. We show the results for all 80 galactic diffuse
background models in gray lines, where we show what the best-fit $\gamma$
value was at a given energy range for each model. With the colored lines, we
show the same results for some of the best overall fit models for the entire
energy range. Right: Similarly, we show how the ellipticity parameter
$\epsilon$ changes with energy. At low energies, due to the large masks and
astrophysical background uncertainties, the GCE morphology, both in terms of
its cuspiness $\gamma$ and ellipticity $\epsilon$, can not be well
constrained. However, with increasing energy as the point sources mask
shrinks, the best-fit models converge to $1.1\leq\gamma\leq 1.3$ and
$\epsilon\simeq 1.0$. Some background models that provide a poor fit to the
data still allow for a wider range of $\gamma$ and $\epsilon$.

###  IV.1 The dependence of the GCE morphology with energy

One key question in determining the origin of the GCE is understanding the
dependence of its morphology on energy. We expect an energy-independent
morphology if the GCE originates purely from dark matter annihilation. This
scenario assumes that gamma rays from the underlying annihilation channel are
prompt emission photons rather than (ICS) photons up-scatted from high-energy
electrons and positrons that interact with the interstellar medium as they
propagate away from their point of origin. Dark matter particles give a
dominant prompt emission component when they annihilate mainly into quarks and
gluons. Instead, if the GCE is purely of stellar origin, it would consistently
show an elongation along the Galactic disk ($\epsilon>1$), regardless of
energy. Such a morphology is consistent with predictions from stellar bulge
profiles, where the gamma-ray emission is predominantly from the outer gap
regions of MSPs and thus directly tracks the location of the MSPs. Finally,
the Fermi bubbles might significantly affect the observed GCE emissions, or
even be related to it, being an earlier event of cosmic-ray burst activity
from the GC Carlson and Profumo (2014); Petrović  _et al._ (2014); Cholis  _et
al._ (2015b). Given the hard spectrum of the Fermi bubbles, their contribution
to the GCE is expected to increase with energy. In this scenario, the GCE will
exhibit a more pronounced elongation perpendicular to the galactic disk at
higher energies, i.e., give ellipticity parameter values of $\epsilon<1$.

In Fig. 5, we show how the GCE morphological properties change with energy.
The left panel, assuming the GCE is spherical, shows the evolution of the
cuspiness parameter $\gamma$ from Eq. 1 with energy. This parameter describes
the level of contraction in the NFW profile. We remind the reader that the GCE
morphology is proportional to the line-of-sight integral of $\rho(r)^{2}$. We
use mask “Standard 4FGLDR3 + L20” and include results for all 80 background
galactic diffuse emission models from Cholis  _et al._ (2022) and highlight in
color the results coming from background galactic diffuse emission models that
provide the best fit to the Fermi data. At low energies, the systematic
uncertainties of the background diffuse emission models prevent a definitive
conclusion about the value of $\gamma$. We remind the reader that various
alternative background models can give similar quality fits. However, for
gamma-ray energies above 2 GeV, there is a robust preference for the GCE to
exhibit a cuspiness of $\gamma\geq 1$, akin to a regular NFW profile or a
slightly adiabatically contracted version. Only a minority of background
models at high energies favor a cuspiness of $\gamma<1$, but these models
generally have poorer overall fit quality for the ROI. In App. A, we provide
the results from the same analysis for alternative ROIs (masks). For the
regular masks, our result on the cuspiness parameter being $\gamma\simeq 1.2$
remains valid. It should be noted that, in addition to the systematic
uncertainties of the galactic diffuse emission backgrounds, accurately
tracking the GCE morphology at lower energies is challenging due to the
significant fraction of the ROI being removed by point source masks, as
detailed in Sec. II.3 and Fig. 2.

In the right panel of Fig. 5, again using “Standard 4FGLDR3 + L20” mask, we
show how the ellipticity parameter $\epsilon$ of Eq. 2 changes with energy.
This parameter indicates the sphericity of the GCE. Again, we show results
from all 80 background diffuse emission models of Ref. Cholis  _et al._ (2022)
and highlight in color those from the best-fit background models. Similar to
the cuspiness analysis, at low energies, the combination of large
astrophysical uncertainties of backgrounds and the small fraction of the
unmasked ROI prevents us from conclusively determining a preference for
$\epsilon$. Nonetheless, at gamma-ray energies above 3 GeV, the background
models that provide a good fit to the data consistently suggest that the GCE
is approximately spherical. In App. A, we show results on the ellipticity
parameter for alternative masks. An ellipticity $\epsilon\simeq 1$ for
energies above 3 GeV is still the case for the regular masks.

###  IV.2 The cuspiness and ellipticity of the GCE for alternative masks

In this section, we test the general morphological properties of the GCE in
terms of its cuspiness and ellipticity under various mask choices.

Figure 6: The GCE cuspiness assuming the emission comes from annihilating dark
matter. Log-likelihoods are evaluated by summing all energy bins. The dark
matter profile is described by the $\gamma$ parameter of Eq. 1. We show
results for alternative mask choices. From the top left panel to the bottom
right, we use masks, “Small 4FGLDR3 + L20”, “Standard 4FGLDR3 + L20”, “Large
4FGLDR3 + L20”, “Standard 4FGLDR3 + L5”, “Standard 4FGLDR3 + Wavelet S4” and
“Standard 4FGLDR3 + Wavelet S2”. The alternative masks (ROIs) may affect the
ordering between the diffuse emission background models that provide the best
fit and the difference in $-2\Delta\ln\mathcal{L}$, but the common pattern is
a systematic preference for $1.0\leq\gamma\leq 1.4$.

Fig. 6 examines the cuspiness $\gamma$ of the GCE. We show the quality of fit
by comparing differences in total negative log-likelihood,
$-2\Delta\ln\mathcal{L}$. This statistic reflects how well the combination of
an astrophysical background and GCE model agrees with the Fermi data for each
of the six masks presented. To encompass the astrophysical uncertainties of
backgrounds, we analyze all 80 astrophysical galactic diffuse emission models
provided in Ref. Cholis  _et al._ (2022). In our plots of
$-2\Delta\ln\mathcal{L}$, we set baseline (zero point) as the combination of
background and GCE model that yields the best overall fit once summing over
all energy bins. Lines closer to zero indicate a higher quality fit. Since we
plot the $-2\Delta\ln\mathcal{L}$ values up to 8000, in some of the panels
(especially the lower ones of Fig. 6), only a subset of lines is visible. The
absent lines are models with poorer fits. For each panel, we color the
combination of background and GCE model that provides the best fit, along with
four more models that provide comparatively similar quality fits.

Direct comparison of log-likelihoods across different masks is not feasible,
as the number of pixels in the remaining ROIs where fitting occurs varies.
Also, the exact ranking of which combination of galactic diffuse emission
background models and GCE provides the best-fit changes slightly with
different ROIs (masks). This variation is expected, mainly because regions
dense in interstellar medium gas may be either revealed or masked by altering
the masks, and these regions are typically associated with high gamma-ray
fluxes. Nonetheless, our objective is to identify the properties of the GCE
morphology that are robust regardless of these varying assumptions. We achieve
this by searching for recurring patterns across the combination of different
choices in the diffuse emission backgrounds and masks.

In the case of Fig. 6, a clear pattern emerges: irrespective of the background
model and the mask used, the GCE consistently shows a preference for
$1.0\leq\gamma\leq 1.4$. This range corresponds to the morphology of a dark
matter annihilation signal, following a spatial distribution akin to a regular
NFW profile ($\gamma=1)$ or a more contracted profile that scales as $\sim
1/r^{\gamma}$ at the inner galaxy.

Figure 7: As in Fig. 6, we show the GCE ellipticity for alternative choices of
masks. Log-likelihoods are again evaluated by summing all energy bins. From
the masks (ROIs) shown in the top panels, the backgrounds giving the best fit
prefer $0.9\leq\epsilon\leq 1.5$. Using instead as a mask the “Standard
4FGLDR3 + L5”, we get $0.7\leq\epsilon\leq 1.2$. The two wavelet-based masks
of “Standard 4FGLDR3 + Wavelet S4” and “Standard 4FGLDR3 + Wavelet S2” give
for a small number of models, but also the ones providing the best fit, a
preference for $0.4\leq\epsilon\leq 0.8$.

These findings are not affected by using any of the remaining four masks
(ROIs) of Fig. 1: “Standard 4FGLDR1 + L20”, “Standard 4FGLDR2 + L20”,
“Standard 4FGLDR3 + L8”, “Standard 4FGLDR3 + Wavelet S3”. The “Standard
4FGLDR1 + L20” and “Standard 4FGLDR2 + L20” masks were tested in Ref. Cholis
_et al._ (2022), yielding the same conclusions as in this work. The “Standard
4FGLDR3 + L8” and “Standard 4FGLDR3 + Wavelet S3” masks represent intermediate
masking assumptions compared to those plotted here.

Fig. 7 focuses on examining the ellipticity $\epsilon$ of the GCE. Again, we
show the quality of fit in terms of $-2\Delta\ln\mathcal{L}$ for the combined
astrophysical background and GCE models against the Fermi data, applying the
same six masks as in Fig. 6. Lines closer to zero in $-2\Delta\ln\mathcal{L}$
represent better quality fits. In each case, we color the combinations of
background and GCE models that provide the best fit, along with four
additional background models yielding similar agreement to the observations.

The ellipticity of GCE appears to be more sensitive to the combination of
assumptions used, compared to its cuspiness. Under the masks shown in the top
panels (“Small 4FGLDR3 + L20”, “Standard 4FGLDR3 + L20”, and “Large 4FGLDR3 +
L20”), the diffuse emission background models that give the best fit prefer
$0.9\leq\epsilon\leq 1.5$. The “Standard 4FGLDR3 + L5” mask yields a similar
range of $0.7\leq\epsilon\leq 1.2$. However, wavelet-based masks, namely
“Standard 4FGLDR3 + Wavelet S4” and “Standard 4FGLDR3 + Wavelet S2”, give a
range of $0.4\leq\epsilon\leq 0.8$ for a few models that provide a better fit.
Such values of $\epsilon$ suggest that the GCE is elongated perpendicular to
the galactic disk, which would indicate the GCE comes from a burst of cosmic
rays similar or the same as the one given in the Fermi bubbles. Such an
ellipticity may also due to a triaxial dark matter profile Diemand  _et al._
(2008); Kuhlen  _et al._ (2008); Springel  _et al._ (2008); Dobler  _et al._
(2011). However, this ellipticity is unlikely to result from a population of
gamma-ray sources tracing the dense stellar regions of the Milky Way.
Additionally, it could also be due to a non-trivial selection effect in which
pixels are identified by the wavelet method as potential point sources.

###  IV.3 Comparing dark matter annihilation and stellar bulge profiles

Figure 8: The “Boxy Bulge” (1st col.), “BB+NB” (2nd col.), “Coleman Bulge”
(4th col.), “F98” (5th col.), and “X-shaped Bulge” (6th col.) profiles used in
this work, together with the “NFW $\gamma=1.2$” profile (3rd col.), for energy
bin 5 ($1.02\,\text{GeV}<E<1.32\,\text{GeV}$). The top row shows the profiles
in the full $40^{\circ}\times 40^{\circ}$ ROI. The bottom row shows the same
profiles with the “L5” disk mask, i.e., masking the region with
$|\ell|<5^{\circ}$ and $|b|<2^{\circ}$. All the profiles are smeared by the
Fermi PSF. The color represents the photon counts of the fiducial profiles
used in the template fitting.

In this section, we compare the hypothesis of dark matter annihilation being
responsible for the GCE with the alternative theory of unresolved gamma-ray
point sources, such as MSPs, being the origin. The latter hypothesis suggests
that the GCE, to some extent, follows one of the known stellar populations,
such as the galactic bulge or the nuclear stellar cluster (nuclear bulge) at
the center of the Milky Way. For the dark matter annihilation hypothesis, we
adopt the NFW profile with $\gamma=1.2$ and $\epsilon=1.0$ (referred to as
“NFW $\gamma=1.2$”), as this choice typically provides the best fit to the
gamma-ray data. We compare this with a sequence of alternative morphological
models proposed in the literature for the GCE.

We compare the dark matter annihilation template and the following stellar
bulge profiles:

  1. (i)

The “Boxy Bulge” profile, based on the stellar bulge measured at the infrared
Launhardt  _et al._ (2002).

  2. (ii)

A combination of the Boxy Bulge and the Nuclear Bulge “BB + NB”. The Nuclear
Bulge is measured at radio waves Launhardt  _et al._ (2002). This combination
has been proposed as the correct morphology for the GCE in Ref. Bartels  _et
al._ (2018).

  3. (iii)

The “Coleman Bulge” from Ref. Coleman  _et al._ (2020), which has been
recently suggested as an alternative for the GCE morphology Macias  _et al._
(2023).

  4. (iv)

The Freudenreich profile for the galactic bar, “F98”. It is measured at the
infrared Freudenreich (1998). This profile has been tested for the GCE
morphology in Ref. Macias  _et al._ (2018) and recently suggested as a
possible alternative for the GCE morphology Macias  _et al._ (2023).

  5. (v)

The “X-shaped Bulge”, observed in lower energy gamma rays and suggested as
having the correct morphology for the GCE in Refs. Macias  _et al._ (2018,
2019).

We showcase the morphology of the stellar bulge profiles, together with the
“NFW $\gamma=1.2$” profile, in Fig. 8. The top row shows the photon count
maps, smeared by the _Fermi_ PSF, in the $1.02-1.32\,\text{GeV}$ energy bin
(energy bin 5) of the fiducial templates. The bottom row shows these same
spatial maps after applying the “L5” disk mask, which covers the brightest
regions across all the profiles and is the largest common masking area used in
the regular masks. With the “L5” disk mask applied, the “Coleman Bulge”
profile bears the greatest similarity to the spherical “NFW $\gamma=1.2$”
among all the stellar bugle profiles.

Figure 9: Similarly to Figs. 6 and 7, for alternative choices of masks we show
how the assumption that the GCE comes from dark matter annihilation (“NFW
$\gamma=1.2$”) compares to various choices of bulges for the GCE (see text for
details). Log-likelihoods are again evaluated by summing all energy bins. We
also include the option that the GCE originates from a combination of dark
matter and a stellar population of sub-threshold point sources (“NFW
$\gamma=1.2$ $\&$ Coleman Bulge”). From the masks (ROIs) shown in the top
panels, the backgrounds giving the best fit prefer the GCE morphology to
follow the “NFW $\gamma=1.2$ $\&$ Coleman Bulge” choice. The two next best
options are either the purely “NFW $\gamma=1.2$” profile or the “Coleman
Bulge” profile (see text for further details). The same result comes for the
lower panel masks, with the difference that for the relevant background model
used to get the best fit–which changes between masks–the pure dark matter
annihilation profile for the GCE performs as well as the “NFW $\gamma=1.2$
$\&$ Coleman Bulge” profile.

The results from template fitting are shown in Fig. 9. We find that the dark
matter annihilation template, “NFW $\gamma=1.2$”, generally outperforms the
“Boxy Bulge”, “BB+NB” and the “X-shaped Bulge” morphologies, irrespective of
the mask or the astrophysical background model for the galactic diffuse
emissions.555For some background models the “Boxy Bulge” or “BB+NB” profiles
may seem to provide a better than “NFW $\gamma=1.2$”. However, these cases
involve background models that yield a significantly poorer quality of fit. It
highlights the importance of background model selection, as a suboptimal
choice can potentially give misleading conclusions about the morphological
properties of the GCE (see also Ref. McDermott  _et al._ (2022)). These
findings strengthen previous results from Ref. Cholis  _et al._ (2022), which
reached similar conclusions regarding the GCE morphology using the “Standard
4FGLDR1 + L20” and “Standard 4FGLDR2 + L20” masks.

Interestingly, the “F98” model’s performance compared to the “NFW
$\gamma=1.2$” varies depending on the specific combination of mask and
background models. The difference in the fit quality between the two models
for the GCE morphology is typically small. However, the “NFW $\gamma=1.2$”
morphology still slightly outperforms the “F98” profile, particularly with
regular masks such as “Standard 4FGLDR3 + L20”, “Small 4FGLDR3 + L20”, “Large
4FGLDR3 + L20”, and “Standard 4FGLDR3 + L5”, for the best-fit background
choices. The preference between “NFW $\gamma=1.2$” and “F98” varies more with
wavelet-based masks, but “NFW $\gamma=1.2$” still outperforms the best-fit
background choice.

When comparing the “NFW $\gamma=1.2$” to the “Coleman Bulge” profile, we find
that the exact choice of mask and galactic diffuse emission background affects
the preference between the two GCE templates. However, the difference in
$-2\Delta\ln\mathcal{L}$ is typically small. As shown in the fourth column of
Fig. 8, the “Coleman Bulge” profile has a morphology that is asymmetric peanut
shape, elongated along the galactic disk. However, when comparing to other
stellar bulge profiles, we note that the “Coleman Bulge” profile is
morphologically the closest to a spherical NFW profile among the bulge
profiles, particularly when masking the central “L5” region. Such similarity
makes the fit result less surprising. We also note that which of those (if
any) bulge profiles more faithfully represents the population of stellar
progenitors to the possible MSPs laying below the point source threshold
remains an open but important question.

Furthermore, we explore whether the GCE could be a combination of dark matter
annihilation and a stellar bulge profile, with free relative normalization for
each of the 14 energy bins we study. This can be justified as even if the GCE
is mostly a dark matter signal, we expect contribution from unresolved MSPs
and other gamma-ray sources at some level. We test combinations of the “NFW
$\gamma=1.2$” profile with (1) the “BB + NB” profile, (2) the “F98” profile,
and (3) the “Coleman Bulge” profile. For the first two cases, the NFW
$\gamma=1.2$ absorbs almost all the GCE-associated emission above 0.7 GeV
effectively in all background model combinations and masks. Furthermore, the
difference in log-likelihood between using the combined profile and the pure
“NFW $\gamma=1.2$” template is small. That is in agreement with earlier
results when we used the “Standard 4FGLDR2 + L20” mask in Refs. Cholis  _et
al._ (2022); McDermott  _et al._ (2022).

However, the combination of “NFW $\gamma=1.2$” with the “Coleman Bulge”
profile (“NFW $\gamma=1.2$ $\&$ Coleman Bulge”) presents a less definitive
picture. Practically in all scenarios, the combined profile “NFW $\gamma=1.2$
$\&$ Coleman Bulge” performs better or as well as either the “NFW
$\gamma=1.2$” or the “Coleman Bulge” profiles alone. The ranking again varies
depending on the mask and the background model. This leads us to explore
further the amount of the GCE emission that is absorbed by each of the two
components (“NFW $\gamma=1.2$” or “Coleman Bulge”) as a function of energy.

###  IV.4 The GCE spectrum from dark matter annihilation and a population of
point sources

In this section, we focus on the case that the GCE is attributed to two
overlapping components, each with its distinct morphology. One component is
modeled by the spherically symmetric dark matter annihilation template, “NFW
$\gamma=1.2$”. The other is described as the “Coleman Bulge” profile. We are
motivated to test the combined “NFW $\gamma=1.2$ $\&$ Coleman Bulge” based on
our findings in Sec. IV.3. This combination often provides the best fit to the
gamma-ray data across various choices of background models and masks (see also
Fig. 9).

Additionally, it is important to consider that gamma-ray point sources below
the detection threshold contribute to the data, and their influence is
generally captured by other templates in our fitting process. The isotropic
diffuse emission template accounts for the extragalactic sub-threshold
sources. The combination of galactic diffuse emission components tends to
absorb some flux from the galactic point sources, especially those correlated
with the gas-dense regions of the Milky Way. By incorporating the “Coleman
Bulge” template, we explicitly test the hypothesis that the GCE is associated
morphologically with the galactic bulge.

Figs. 10 and 11 show the flux decomposition of the GCE into these two
components for the six different masks previously used in Figs. 6-9. Each
component in our study has an independently adjustable normalization for each
of the 14 energy bins. The flux attributed to the dark matter component is
shown with a purple line for the best-fit normalization and a magenta band for
its 2$\sigma$ range. Similarly, the flux associated with the “Coleman Bulge”
is given by a light green line for the best-fit normalization and a dark green
band for its 2$\sigma$ range. For each mask, we show results from two out of
the 80 background models. We pick these two models among those providing a
good fit to the data and in a manner that is intended to envelope for the
reader the range of how much of the GCE each of the two components can absorb.

Figure 10: Assuming that the GCE comes from the combination of dark matter and
subthreshold gamma-ray point sources, we show the flux decomposition of the
GCE attributed to the spherically symmetric ($\epsilon=1.0$) dark matter
annihilation template with $\gamma=1.2$ (purple line and magenta band) and the
asymmetric and elongated along the galactic disk “Coleman Bulge” template
(green line and band). We show the evolution with the energy of the two
components contributing to the GCE. In each row, we show results from two
different gamma-ray background models keeping the mask fixed. The two
background models are selected to be among those providing the best agreement
to the data. From top to bottom, we use the following masks: “Small 4FGLDR3 +
L20”, “Standard 4FGLDR3 + L20”, “Large 4FGLDR3 + L20” and “Standard 4FGLDR3 +
L5”.

Figure 11: As in Fig. 10, we show the flux decomposition of the GCE attributed
to the dark matter annihilation template with $\gamma=1.2$ and $\epsilon=1.0$
(purple line and magenta band) and the “Coleman Bulge” template (green line
and band). We show results using masks “Standard 4FGLDR3 + Wavelet S4” and
“Standard 4FGLDR3 + Wavelet S2”.

In the analysis using regular masks, as shown in Fig. 10, the dark matter
component consistently comprises at least 50$\%$ of the GCE emission above 0.7
GeV, and in some instances, it accounts for over 90$\%$ of the GCE emission.
Moreover, even after accounting for systematic uncertainties, the spectrum of
the dark matter component is roughly the same as what has been proposed in
previous works as Hooper and Goodenough (2011); Vitale and Morselli (2009);
Abazajian (2011); Hooper and Linden (2011); Hooper and Slatyer (2013); Gordon
and Macias (2013); Abazajian  _et al._ (2014); Daylan  _et al._ (2016); Calore
_et al._ (2015a); Zhou  _et al._ (2015); Ajello  _et al._ (2016); Linden  _et
al._ (2016); Cholis  _et al._ (2022). This suggests that, within a factor of
$\simeq 2$ adjustments in the dark matter cross section and mass, the existing
explanations for the dark matter origin of the GCE remain valid, even
considering all tested background diffuse emission model systematics, mask
choices, and the potential independent contribution of a galactic bulge
component.

It is also notable that even when the “Coleman Bulge” component contributes
maximally (as seen in the right panels of Fig. 10), its spectrum is distinctly
different from that typically associated with MSPs Abazajian (2011); Abazajian
_et al._ (2014); Cholis  _et al._ (2014); Bartels  _et al._ (2016). Instead,
it more closely resembles the gamma-ray spectrum from the regular population
of galactic point sources Abdollahi  _et al._ (2022), which includes MSPs but
also contains the emission from younger gamma-ray pulsars, including pulsar
wind nebulae, sources in globular clusters, binaries, supernova remnants, and
star-forming regions. This suggests that the actual contribution of regular
MSPs to the GCE may be relatively minor (see e.g. Hooper  _et al._ (2013);
Cholis  _et al._ (2015a); Zhong  _et al._ (2020); Leane and Slatyer (2019)).
Results from the “Standard 4FGLDR1 + L20”, “Standard 4FGLDR2 + L20” and
“Standard 4FGLDR3 + L8” masks, lead to the same conclusions. In App. B, we
provide additional information on the statistical preference for either the
dark matter or the “Coleman Bulge” components of the GCE as a function of
energy.

Fig. 11, shows the results obtained using our wavelet-based masks “Standard
4FGLDR3 + Wavelet S2” and “Standard 4FGLDR3 + Wavelet S4”. The GCE spectra
derived from these masks are distinctively different from those derived when
portions of the galactic disk are masked. Depending on the specific case, the
dark matter component or the “Coleman Bulge” may dominate the GCE emission.
However, it is important to note that the gamma-ray spectrum associated with
the “Coleman Bulge”, again, does not resemble the known spectra of MSPs. This
discrepancy suggests a significant contribution from unresolved point sources
along the disk. Results from the “Standard 4FGLDR3 + Wavelet S3” mask agree
with those of the other two wavelet-based masks.

Figure 12: For the $40^{\circ}\times 40^{\circ}$ ROI with the “Standard
4FGLDR3 + L20” mask, we provide all galactic diffuse emission components.
Those are the Pi0+Bremsstrahlung in orange lines and band, the ICS in blue
lines and band, the isotropic in turquoise, Fermi bubbles in orange, the dark
matter component in purple line and magenta band, and the “Coleman Bulge”
component in green line and band. Left panel: We use galactic diffuse emission
model I of Ref. Cholis  _et al._ (2022). Right panel: We instead use the model
XLVIII from the same reference (see text for further details).

Finally, in Fig. 12, as a point of reference for the reader, we provide a
complete decomposition of gamma-ray emission from the $40^{\circ}\times
40^{\circ}$ ROI for one of the masks and two galactic diffuse emission
background models. The solid lines represent the original assumptions of the
Pi0+Bremsstrahlung and the ICS before fitting each energy bin. The dashed
lines show the fitted spectra for each component, along with their
corresponding 2$\sigma$ ranges as colored bands. On the left panel of Fig. 12,
we use the galactic diffuse emission model I of Ref. Cholis  _et al._ (2022),
where the dark matter template dominates the GCE. On the right panel, we use
model XLVIII from the same reference instead, in which case the dark matter
template contributes approximately half of the GCE flux above 1 GeV. Note that
in the case of the right panel, the flux component associated with the
“Coleman Bulge” profile does not have the hard spectrum followed by an
exponential cut-off expected by MSPs. Instead, the spectrum associated with
the “Coleman Bulge” is that expected from the population of all galactic
gamma-ray point sources. Furthermore, the inferred isotropic emission
component deviates from the true isotropic emission calculated at high
latitudes. On the left panel, the “Coleman Bulge” template component has a
similar spectrum to that on the right but an order of magnitude dimmer. We
find such spectral patterns indicative of the “Coleman Bulge” template are
likely some combination of extragalactic and galactic point sources and not an
unresolved millisecond pulsar population.

##  V Discussions and Conclusions

Masking is an essential step for the template-fitting analysis of the GCE,
crucial for minimizing the influence of the galactic disk and the known point
sources on the GCE emission, which are challenging to model accurately. In
this study, we constructed a variety of masks based on the information of the
point sources and the disk. These include regular masks, which combine a band-
shaped disk mask with point source masks based on the Fermi 4FGL catalogs, and
wavelet-based masks, which combine masks for high-significance peaks
identified by the wavelet method with masks for point sources also based on
4FGL. For the regular masks, we varied the band length for the disk mask and
the radii for the point source masks.

We then extensively investigated the impact of masking on the GCE morphology
analysis. Through injection tests, we demonstrated that our analysis method,
the template fitting, is sensitive and robust against the ellipticity of the
GCE. Through systematic analyses, we further showed that the choice of masks
does not significantly impact the inferred cuspiness and ellipticity of the
GCE, assuming the GCE follows contracted, squared, and integrated NFW
profiles, a common template to model dark matter annihilation. This robustness
in the characteristics of the GCE is evident at an energy bin-by-energy bin
level, except for the lowest energy bins, due to limited statistics after
masking and large astrophysical background uncertainties.

We also explored various hypotheses about the GCE’s morphology based on
different stellar bulge profiles, including “Boxy Bulge”, “BB+NB”, “Coleman
Bulge”, “F98”, and “X-shaped Bulge”. We compared these results to dark matter
annihilation (“NFW $\gamma=1.2$”). Across different choices of masks and
background diffuse emission models, we found that the “NFW $\gamma=1.2$” and
the “Coleman Bulge” profiles offer the best fits to the GCE, with the latter
profile giving slightly more often a better fit to the gamma-ray data, but
with a typically small difference in terms of $-2\Delta\ln\mathcal{L}$.

Additionally, we tested combinations of dark matter annihilation with the
“BB+NB”, the “F98”, and the “Coleman Bulge”. Among those combinations, we
found that “NFW $\gamma=1.2$ & Coleman Bulge” offers the best fit across
different choices of masks and background emission models. It also
consistently fits better than “NFW $\gamma=1.2$” or “Coleman Bulge” alone.
This prompted further exploration of the contributions from the two components
as a function of energy under different masks. We found that the dark matter
component consistently comprises a significant portion of the GCE for
$E>0.7\,\text{GeV}$. Although, in some cases, the “Coleman Bulge” can be an
important component, the corresponding spectra distinctly differ from that
typically associated with MSPs, resembling more the spectra of the regular
galactic point sources. Besides, when the “Coleman Bulge” is an important
component, the resulting spectrum associated with the isotropic emission
deviates from the measurement at high latitudes. The fact that the “Coleman
Bulge” does not yield a MSP-like spectrum is the case across both regular and
wavelet-based masks. All regular masks still give a prominent spherically
symmetric GCE component.

Our results using wavelet-based masks occasionally show a different GCE
morphology and spectrum from those obtained with regular masks, possibly due
to too much contribution from other galactic point sources. However, the
template fitting based on wavelet-based masks generally maintains a consistent
ranking of GCE hypotheses (dark matter annihilation, stellar bulges, and
combinations) as those based on regular masks. The wavelet-based masks do
accentuate the differences in $-2\Delta\ln\mathcal{L}$ between the best-fit
background models and poorer-fit ones. Further exploration with a broader pool
of background models may offer more suitable candidates for template fitting
under wavelet-based masks.

###### Acknowledgements.

We thank Samuel D. McDermott for participating in the early stage of the
project. We thank Matthew Buckley, Dan Hooper, Manoj Kaplinghat, Tracy
Slatyer, Deheng Song, and Yitian Sun for valuable discussions and Oscar Macias
and Deheng Song for providing us with the F98 and Coleman Bulge templates. We
also acknowledge the use of GALPROP Strong and Moskalenko (1998). This
research uses computational resources at the University of Chicago’s Research
Computing Center. We thank Edward W. Kolb for providing access to the
resources. YZ acknowledges the Aspen Center for Physics for its hospitality
during the completion of this study, which is supported by the National
Science Foundation under Grant PHY-1607611. YZ is supported by the Kavli
Institute for Cosmological Physics at the University of Chicago through an
endowment from the Kavli Foundation and its founder Fred Kavli. IC
acknowledges that this material is based upon work supported by the U.S.
Department of Energy, Office of Science, Office of High Energy Physics, under
Award No. DE-SC0022352.

##  Appendix A The GCE Morphology’s Energy Evolution for alternative masks

In Sec. IV.1 of the main text, we presented for the “Standard 4FGLDR3 + L20”
mask the evolution with energy of the GCE’s cuspiness and ellipticity. Here,
we show results for the other five masks, for which figures are presented in
the main text. Those are the, “Small 4FGLDR3 + L20”, “Large 4FGLDR3 + L20”,
“Standard 4FGLDR3 + L25”, “Standard 4FGLDR3 + Wavelet S4” and “Standard
4FGLDR3 + Wavelet S2” masks. For ease of comparison, we also include the
results of “Standard 4FGLDR3 + L20”.

In Fig. 13, we show the cuspiness evolution with energy for these six masks.
Using the regular masks (first four panels) and focusing on the background
diffuse emission models that give the best fits (colored lines), the cuspiness
converges to values of $1.1\leq\gamma\leq 1.4$ for energies 3 GeV or higher.
Only the wavelet-based masks provide inconclusive results on the GCE
cuspiness. As described in the main text, we think that to be the result of
the wavelet-based masks allowing too much contribution of other galactic (non-
MSP) point sources, whose spectrum has nothing to do with the GCE emission.

Figure 13: As in the left panel of Figure 5, the evolution of the cuspiness of
the GCE with energy for a sequence of masks (see text for more details).

In Fig. 14, we show the ellipticity evolution with energy for the same six
masks. Again, the regular masks, in combination with the best-fit background
models, converge to values of $0.8\leq\epsilon\leq 1.2$ for energies $>3$ GeV.
Instead, the wavelet-based masks provide inconclusive results on the GCE
ellipticity.

Figure 14: As in the right panel of Figure 5, the evolution of the ellipticity
of the GCE with energy for a sequence of masks (see text for more details).

##  Appendix B The energy evolution of the statistical preference for the GCE
components

In the main text, in Secs. IV.3 and IV.4 and Figs. 9, 10 and 11 we
investigated the possibility that the GCE is a combination of both a component
following the dark matter annihilation profile and a component that traces the
“Coleman Bulge” profile. As we showed in the main text, the component tracing
the dark matter annihilation profile is always present, with a normalization
that for energies $E>0.7\,\text{GeV}$ accounts for $\simeq 50-90\%$ of the
total GCE flux. Only when we use the “Standard 4FGLD3 + Wavelet S2” mask, we
find a lower flux associated with the dark matter annihilation profile.

In this appendix, we provide further results on the evolution with the energy
of the statistical preference for the GCE to trace either the dark matter
annihilation morphology or the “Coleman Bulge” one. We test the cases where
the GCE is purely tracing one or the other morphology. In Fig. 15, we show for
the same six masks as in section IV and for specific galactic diffuse
background models the energy bin-by-energy bin difference in negative log-
likelihood (in $-2\Delta\ln\mathcal{L}_{j}$) between the dark matter
annihilation “NFW with $\gamma=1.2$” profile and the “Coleman Bulge” profile.
The biggest differences show up at energies between 1 GeV and 5 GeV, with
little statistical difference between the two profiles for $E>10$ GeV where
the number of detected photons is small. Also, at $E<0.7$ GeV, there is
typically little statistical preference between the two options. At the lowest
energies, the background emission is too dominant.

Figure 15: For a sample of galactic diffuse models and masks, we show the
energy bin-by-energy bin level of statistical preference between the dark
matter annihilation “NFW with $\gamma=1.2$” profile vs. the “Coleman Bulge”
profile.

## References

  * Goodenough and Hooper (2009) Lisa Goodenough and Dan Hooper, “Possible Evidence For Dark Matter Annihilation In The Inner Milky Way From The Fermi Gamma Ray Space Telescope,” (2009), arXiv:0910.2998 [hep-ph] .
  * Hooper and Goodenough (2011) Dan Hooper and Lisa Goodenough, “Dark Matter Annihilation in The Galactic Center As Seen by the Fermi Gamma Ray Space Telescope,” Phys. Lett. B697, 412–428 (2011), arXiv:1010.2752 [hep-ph] .
  * Vitale and Morselli (2009) Vincenzo Vitale and Aldo Morselli, “Indirect Search for Dark Matter from the center of the Milky Way with the Fermi-Large Area Telescope,” arXiv e-prints , arXiv:0912.3828 (2009), arXiv:0912.3828 [astro-ph.HE] .
  * Abazajian (2011) Kevork N. Abazajian, “The Consistency of Fermi-LAT Observations of the Galactic Center with a Millisecond Pulsar Population in the Central Stellar Cluster,” JCAP 1103, 010 (2011), arXiv:1011.4275 [astro-ph.HE] .
  * Hooper and Linden (2011) Dan Hooper and Tim Linden, “On The Origin Of The Gamma Rays From The Galactic Center,” Phys. Rev. D 84, 123005 (2011), arXiv:1110.0006 [astro-ph.HE] .
  * Hooper and Slatyer (2013) Dan Hooper and Tracy R. Slatyer, “Two Emission Mechanisms in the Fermi Bubbles: A Possible Signal of Annihilating Dark Matter,” Phys. Dark Univ. 2, 118–138 (2013), arXiv:1302.6589 [astro-ph.HE] .
  * Gordon and Macias (2013) Chris Gordon and Oscar Macias, “Dark Matter and Pulsar Model Constraints from Galactic Center Fermi-LAT Gamma Ray Observations,” Phys. Rev. D 88, 083521 (2013), [Erratum: Phys.Rev.D 89, 049901 (2014)], arXiv:1306.5725 [astro-ph.HE] .
  * Abazajian  _et al._ (2014) Kevork N. Abazajian, Nicolas Canac, Shunsaku Horiuchi, and Manoj Kaplinghat, “Astrophysical and Dark Matter Interpretations of Extended Gamma-Ray Emission from the Galactic Center,” Phys. Rev. D90, 023526 (2014), arXiv:1402.4090 [astro-ph.HE] .
  * Daylan  _et al._ (2016) Tansu Daylan, Douglas P. Finkbeiner, Dan Hooper, Tim Linden, Stephen K. N. Portillo, Nicholas L. Rodd, and Tracy R. Slatyer, “The characterization of the gamma-ray signal from the central Milky Way: A case for annihilating dark matter,” Phys. Dark Univ. 12, 1–23 (2016), arXiv:1402.6703 [astro-ph.HE] .
  * Calore  _et al._ (2015a) Francesca Calore, Ilias Cholis, and Christoph Weniger, “Background Model Systematics for the Fermi GeV Excess,” JCAP 1503, 038 (2015a), arXiv:1409.0042 [astro-ph.CO] .
  * Zhou  _et al._ (2015) Bei Zhou, Yun-Feng Liang, Xiaoyuan Huang, Xiang Li, Yi-Zhong Fan, Lei Feng, and Jin Chang, “GeV excess in the Milky Way: The role of diffuse galactic gamma-ray emission templates,” Phys. Rev. D 91, 123010 (2015), arXiv:1406.6948 [astro-ph.HE] .
  * Ajello  _et al._ (2016) M. Ajello _et al._ (Fermi-LAT), “Fermi-LAT Observations of High-Energy $\gamma$-Ray Emission Toward the Galactic Center,” Astrophys. J. 819, 44 (2016), arXiv:1511.02938 [astro-ph.HE] .
  * Linden  _et al._ (2016) Tim Linden, Nicholas L. Rodd, Benjamin R. Safdi, and Tracy R. Slatyer, “High-energy tail of the Galactic Center gamma-ray excess,” Phys. Rev. D 94, 103013 (2016), arXiv:1604.01026 [astro-ph.HE] .
  * Cholis  _et al._ (2022) Ilias Cholis, Yi-Ming Zhong, Samuel D. McDermott, and Joseph P. Surdutovich, “Return of the templates: Revisiting the Galactic Center excess with multimessenger observations,” Phys. Rev. D 105, 103023 (2022), arXiv:2112.09706 [astro-ph.HE] .
  * Calore  _et al._ (2015b) Francesca Calore, Ilias Cholis, Christopher McCabe, and Christoph Weniger, “A Tale of Tails: Dark Matter Interpretations of the Fermi GeV Excess in Light of Background Model Systematics,” Phys. Rev. D 91, 063003 (2015b), arXiv:1411.4647 [hep-ph] .
  * Agrawal  _et al._ (2015) Prateek Agrawal, Brian Batell, Patrick J. Fox, and Roni Harnik, “WIMPs at the Galactic Center,” JCAP 05, 011 (2015), arXiv:1411.2592 [hep-ph] .
  * Berlin  _et al._ (2015) Asher Berlin, Stefania Gori, Tongyan Lin, and Lian-Tao Wang, “Pseudoscalar Portal Dark Matter,” Phys. Rev. D 92, 015005 (2015), arXiv:1502.06000 [hep-ph] .
  * Karwin  _et al._ (2017) Christopher Karwin, Simona Murgia, Tim M. P. Tait, Troy A. Porter, and Philip Tanedo, “Dark Matter Interpretation of the Fermi-LAT Observation Toward the Galactic Center,” Phys. Rev. D 95, 103005 (2017), arXiv:1612.05687 [hep-ph] .
  * Ackermann  _et al._ (2017) M. Ackermann _et al._ (Fermi-LAT), “The Fermi Galactic Center GeV Excess and Implications for Dark Matter,” Astrophys. J. 840, 43 (2017), arXiv:1704.03910 [astro-ph.HE] .
  * Leane and Slatyer (2019) Rebecca K. Leane and Tracy R. Slatyer, “Revival of the Dark Matter Hypothesis for the Galactic Center Gamma-Ray Excess,” Phys. Rev. Lett. 123, 241101 (2019), arXiv:1904.08430 [astro-ph.HE] .
  * Abazajian and Kaplinghat (2012) Kevork N. Abazajian and Manoj Kaplinghat, “Detection of a Gamma-Ray Source in the Galactic Center Consistent with Extended Emission from Dark Matter Annihilation and Concentrated Astrophysical Emission,” Phys. Rev. D 86, 083511 (2012), [Erratum: Phys.Rev.D 87, 129902 (2013)], arXiv:1207.6047 [astro-ph.HE] .
  * Hooper  _et al._ (2013) Dan Hooper, Ilias Cholis, Tim Linden, Jennifer Siegal-Gaskins, and Tracy Slatyer, “Pulsars Cannot Account for the Inner Galaxy’s GeV Excess,” Phys. Rev. D 88, 083009 (2013), arXiv:1305.0830 [astro-ph.HE] .
  * Lee  _et al._ (2016) Samuel K. Lee, Mariangela Lisanti, Benjamin R. Safdi, Tracy R. Slatyer, and Wei Xue, “Evidence for Unresolved $\gamma$-Ray Point Sources in the Inner Galaxy,” Phys. Rev. Lett. 116, 051103 (2016), arXiv:1506.05124 [astro-ph.HE] .
  * Petrović  _et al._ (2015) Jovana Petrović, Pasquale D. Serpico, and Gabrijela Zaharijas, “Millisecond pulsars and the Galactic Center gamma-ray excess: the importance of luminosity function and secondary emission,” JCAP 02, 023 (2015), arXiv:1411.2980 [astro-ph.HE] .
  * Cholis  _et al._ (2015a) Ilias Cholis, Dan Hooper, and Tim Linden, “Challenges in Explaining the Galactic Center Gamma-Ray Excess with Millisecond Pulsars,” JCAP 06, 043 (2015a), arXiv:1407.5625 [astro-ph.HE] .
  * Bartels  _et al._ (2018) Richard Bartels, Emma Storm, Christoph Weniger, and Francesca Calore, “The Fermi-LAT GeV excess as a tracer of stellar mass in the Galactic bulge,” Nature Astron. 2, 819–828 (2018), arXiv:1711.04778 [astro-ph.HE] .
  * Buschmann  _et al._ (2020) Malte Buschmann, Nicholas L. Rodd, Benjamin R. Safdi, Laura J. Chang, Siddharth Mishra-Sharma, Mariangela Lisanti, and Oscar Macias, “Foreground Mismodeling and the Point Source Explanation of the Fermi Galactic Center Excess,” Phys. Rev. D 102, 023023 (2020), arXiv:2002.12373 [astro-ph.HE] .
  * Gautam  _et al._ (2022) Anuj Gautam, Roland M. Crocker, Lilia Ferrario, Ashley J. Ruiter, Harrison Ploeg, Chris Gordon, and Oscar Macias, “Millisecond pulsars from accretion-induced collapse as the origin of the Galactic Centre gamma-ray excess signal,” Nature Astron. 6, 703–707 (2022), arXiv:2106.00222 [astro-ph.HE] .
  * Petrović  _et al._ (2014) Jovana Petrović, Pasquale Dario Serpico, and Gabrijela Zaharijaš, “Galactic Center gamma-ray ”excess” from an active past of the Galactic Centre?” JCAP 10, 052 (2014), arXiv:1405.7928 [astro-ph.HE] .
  * Carlson and Profumo (2014) Eric Carlson and Stefano Profumo, “Cosmic Ray Protons in the Inner Galaxy and the Galactic Center Gamma-Ray Excess,” Phys. Rev. D 90, 023015 (2014), arXiv:1405.7685 [astro-ph.HE] .
  * Cholis  _et al._ (2015b) Ilias Cholis, Carmelo Evoli, Francesca Calore, Tim Linden, Christoph Weniger, and Dan Hooper, “The Galactic Center GeV Excess from a Series of Leptonic Cosmic-Ray Outbursts,” JCAP 12, 005 (2015b), arXiv:1506.05119 [astro-ph.HE] .
  * Zhong  _et al._ (2020) Yi-Ming Zhong, Samuel D. McDermott, Ilias Cholis, and Patrick J. Fox, “Testing the Sensitivity of the Galactic Center Excess to the Point Source Mask,” Phys. Rev. Lett. 124, 231103 (2020), arXiv:1911.12369 [astro-ph.HE] .
  * Bartels  _et al._ (2016) Richard Bartels, Suraj Krishnamurthy, and Christoph Weniger, “Strong support for the millisecond pulsar origin of the Galactic center GeV excess,” Phys. Rev. Lett. 116, 051102 (2016), arXiv:1506.05104 [astro-ph.HE] .
  * Macias  _et al._ (2018) Oscar Macias, Chris Gordon, Roland M. Crocker, Brendan Coleman, Dylan Paterson, Shunsaku Horiuchi, and Martin Pohl, “Galactic bulge preferred over dark matter for the Galactic centre gamma-ray excess,” Nature Astron. 2, 387–392 (2018), arXiv:1611.06644 [astro-ph.HE] .
  * Macias  _et al._ (2019) Oscar Macias, Shunsaku Horiuchi, Manoj Kaplinghat, Chris Gordon, Roland M. Crocker, and David M. Nataf, “Strong Evidence that the Galactic Bulge is Shining in Gamma Rays,” JCAP 09, 042 (2019), arXiv:1901.03822 [astro-ph.HE] .
  * Di Mauro (2021) Mattia Di Mauro, “Characteristics of the Galactic Center excess measured with 11 years of $Fermi$-LAT data,” Phys. Rev. D 103, 063029 (2021), arXiv:2101.04694 [astro-ph.HE] .
  * Pohl  _et al._ (2022) Martin Pohl, Oscar Macias, Phaedra Coleman, and Chris Gordon, “Assessing the Impact of Hydrogen Absorption on the Characteristics of the Galactic Center Excess,” Astrophys. J. 929, 136 (2022), arXiv:2203.11626 [astro-ph.HE] .
  * McDermott  _et al._ (2022) Samuel D. McDermott, Yi-Ming Zhong, and Ilias Cholis, “A Phantom Menace: On the Morphology of the Galactic Center Excess,” (2022), arXiv:2209.00006 [astro-ph.HE] .
  * Abdollahi  _et al._ (2020) S. Abdollahi _et al._ (Fermi-LAT), “$Fermi$ Large Area Telescope Fourth Source Catalog,” Astrophys. J. Suppl. 247, 33 (2020), arXiv:1902.10045 [astro-ph.HE] .
  * Ballet  _et al._ (2020) J. Ballet, T.H. Burnett, S.W. Digel, and B. Lott (Fermi-LAT), “Fermi Large Area Telescope Fourth Source Catalog Data Release 2,” (2020), arXiv:2005.11208 [astro-ph.HE] .
  * Abdollahi  _et al._ (2022) Soheila Abdollahi _et al._ (Fermi-LAT), “Incremental Fermi Large Area Telescope Fourth Source Catalog,” Astrophys. J. Supp. 260, 53 (2022), arXiv:2201.11184 [astro-ph.HE] .
  * Ballet  _et al._ (2023) J. Ballet, P. Bruel, T. H. Burnett, and B. Lott (Fermi-LAT), “Fermi Large Area Telescope Fourth Source Catalog Data Release 4 (4FGL-DR4),” (2023), arXiv:2307.12546 [astro-ph.HE] .
  * Coleman  _et al._ (2020) Brenna Coleman, Dylan Paterson, Chris Gordon, Oscar Macias, and Harrison Ploeg, “Maximum Entropy Estimation of the Galactic Bulge Morphology via the VVV Red Clump,” Mon. Not. Roy. Astron. Soc. 495, 3350–3372 (2020), arXiv:1911.04714 [astro-ph.GA] .
  * psf (2019) “Fermi-LAT Performance,”  (2019).
  * Ackermann  _et al._ (2014) M. Ackermann _et al._ (Fermi-LAT), “The Spectrum and Morphology of the $Fermi$ Bubbles,” Astrophys. J. 793, 64 (2014), arXiv:1407.7905 [astro-ph.HE] .
  * Ackermann  _et al._ (2015) M. Ackermann _et al._ (Fermi-LAT), “The spectrum of isotropic diffuse gamma-ray emission between 100 MeV and 820 GeV,” Astrophys. J. 799, 86 (2015), arXiv:1410.3696 [astro-ph.HE] .
  * Foreman-Mackey  _et al._ (2013) Daniel Foreman-Mackey, David W. Hogg, Dustin Lang, and Jonathan Goodman, “emcee: The MCMC Hammer,” Publ. Astron. Soc. Pac. 125, 306–312 (2013), arXiv:1202.3665 [astro-ph.IM] .
  * Hinton (2016) S. R. Hinton, “ChainConsumer,” The Journal of Open Source Software 1, 00045 (2016).
  * Diemand  _et al._ (2008) J. Diemand, M. Kuhlen, P. Madau, M. Zemp, B. Moore, D. Potter, and J. Stadel, “Clumps and streams in the local dark matter distribution,” Nature 454, 735–738 (2008), arXiv:0805.1244 [astro-ph] .
  * Kuhlen  _et al._ (2008) M. Kuhlen, J. Diemand, P. Madau, and M. Zemp, “The Via Lactea INCITE simulation: galactic dark matter substructure at high resolution,” in  _Journal of Physics Conference Series_, Journal of Physics Conference Series, Vol. 125 (2008) p. 012008, arXiv:0810.3614 [astro-ph] .
  * Springel  _et al._ (2008) V. Springel, J. Wang, M. Vogelsberger, A. Ludlow, A. Jenkins, A. Helmi, J. F. Navarro, C. S. Frenk, and S. D. M. White, “The Aquarius Project: the subhaloes of galactic haloes,” Mon. Not. R. Astron. Soc. 391, 1685–1711 (2008), arXiv:0809.0898 [astro-ph] .
  * Dobler  _et al._ (2011) Gregory Dobler, Ilias Cholis, and Neal Weiner, “The Fermi Gamma-Ray Haze from Dark Matter Annihilations and Anisotropic Diffusion,” Astrophys. J. 741, 25 (2011), arXiv:1102.5095 [astro-ph.HE] .
  * Launhardt  _et al._ (2002) R. Launhardt, R. Zylka, and P. G. Mezger, “The nuclear bulge of the Galaxy. III. Large-scale physical characteristics of stars and interstellar matter,” Astron. Astrophys. 384, 112–139 (2002), arXiv:astro-ph/0201294 [astro-ph] .
  * Macias  _et al._ (2023) Oscar Macias, Deheng Song, and Manoj Kaplinghat, “Private communication,” (2023).
  * Freudenreich (1998) H. T. Freudenreich, “Cobe’s galactic bar and disk,” Astrophys. J. 492, 495–510 (1998), arXiv:astro-ph/9707340 .
  * Cholis  _et al._ (2014) Ilias Cholis, Dan Hooper, and Tim Linden, “A New Determination of the Spectra and Luminosity Function of Gamma-Ray Millisecond Pulsars,” (2014), arXiv:1407.5583 [astro-ph.HE] .
  * Strong and Moskalenko (1998) A. W. Strong and I. V. Moskalenko, “Propagation of cosmic-ray nucleons in the Galaxy,” Astrophys. J. 509, 212–228 (1998), astro-ph/9807150 .