Characterizing γ-ray maps of the Galactic Centerwith neural density estimation

Siddharth Mishra-Sharma
The NSF AI Institute for Artificial Intelligence and Fundamental Interactions;
Massachusetts Institute of Technology;
Harvard University;
New York University
smsharma@mit.edu
Kyle Cranmer
New York University
kyle.cranmer@nyu.edu

Abstract
Machine learning methods have enabled new ways of performing inference on
high-dimensional datasets modeled using complex simulations. We leverage recent
advancements in simulation-based inference in order to characterize the contribution of various modeled components to γ-ray data of the Galactic Center recordedby the Fermi satellite. A specific goal here is to differentiate “smooth” emission, as
expected for a dark matter origin, from more “clumpy” emission expected for a
population of relatively bright, unresolved astrophysical point sources. Compared
to traditional techniques based on the statistical distribution of photon counts, ourmethod based on density estimation using normalizing flows is able to utilize moreof the information contained in a given model of the Galactic Center emission,and in particular can perform posterior parameter estimation while accounting for
pixel-to-pixel spatial correlations in the γ-ray map.

1

Introduction

Dark matter (DM) represents one of the major unsolved problems in particle physics and cosmology
today. The traditional Weakly-Interacting Massive Particle (WIMP) paradigm envisions productionof dark matter in the early Universe through freeze-out of dark sector particles weakly coupled to theStandard Model (SM) sector. In this scenario, one of the most promising avenues of detecting a dark
matter signal is through an observation of excess γ-ray photons at ∼ GeV energies from DM-richregions of the sky. The Fermi γ-ray Galactic Center Excess (GCE), first identified over a decadeago [1–4] using data from the Fermi Large Area Telescope (LAT) [5], is an excess of photons in
the Galactic Center with properties—such as energy spectrum and spatial morphology—broadly
compatible with expectation due to annihilating DM [6, 7].
The high dimensionality of γ-ray data has traditionally necessitated a description of the photon map
in terms of hand-crafted summary quantities e.g., the probability distribution of photon counts [8, 9]
or a wavelet decomposition of the photon map [10–13], in order to enable computationally tractable
analyses. While effective, this reduced description necessarily involves loss of information comparedto that contained in the original γ-ray map. On the other hand, recent developments in machinelearning have enabled analysis techniques that can extract more information from high-dimensional
datasets. Machine learning methods have recently shown promise for analyzing γ-ray data [14] andspecifically for understanding the nature of the Fermi GCE [15–17].
Fourth Workshop on Machine Learning and the Physical Sciences (NeurIPS 2021).

Your Paper
You
October 8, 2021
Abstract
Your abstract.
!

1

Normalizing flow

Introduction

Base distribution

Posterior distributionu ⇠ ⇡(u) = N (u; 0, 1)✓ ⇠ p(✓ | x)
Your introduction goes here! Simply start writing your document and use the Recompile button✓
to view the updated PDF preview. Examples of commonly used commands andPS
features are listed
12 params.
below, to help you get started.
Once you’re familiar with the editor, you can find various project setting in the Overleaf menu,
accessed via the button in the
very
top
left
of
the
editor.
To
view
tutorials,
user guides, and further
Flow transformations
documentation, please visit our help library, or head to our plans page to choose your plan.

·

···

··

✓Poiss6 params.

Conditioning context

How to create Sections and Subsections

···

2.1

·

···

Some examples to get started··

f8 (u7 ; s) · · · f2 (u1 ; s) f1 (u; s)2

Simply use the section
and subsection commands, as in this example document! With
Overleaf, all······
the formatting and numbering is handled automatically according to the template you’ve chosen. If
you’re using Rich Text mode, you can also create new section and subsections via the buttons in theeditor toolbar.

e=

···

2.3

32 ch.

···

Inferred quantities

How to add Tables

Component intensities
GCE PSGCE DM
Disk PS
Isotropic
Bubbles

Bubbles

···

id

64
e=
id
ns

How to include Figures

First you have to upload the image file from your computer using the upload link in the file-tree menu.
Then use the includegraphics
command to include it in your document.
Use the distributions
figure environment
1
Source-count=2 caption commande=
PS 1 in
anddetheto add a number and a caption to your figure. See the code for GCE
Figureid
i
DM
GCE PS
thisns section for anns example.
GCE
DM theDisk
PS
Note that your figure will automatically be placed in the most appropriate place for it,
given
Disk
PS
Isotropic
surrounding text and taking into account other figures or tables that may be close by. You
canfind
Isotropic
Bubblesout more about adding images to your documents in this help article on including images onOverleaf.

Number density of PSs

ns

8
12

2.2

Feature extractor

Summaries s(x)Input map x

Use the table and tabular environments for basic tables — see Table 1, for example. For more information, please see this help article on tables. 2 R128

256 ch.

256 ch.

Convolutional (7 layers)

Forward model x ⇠ p(x | ✓)Item
Fully-connected
Widgets
Gadgets

Quantity
42
13

PS flux

0

1 Intensity2

3

Table 1: An example table.

Figure 1: A schematic overview of the inference framework used in this work. A normalizing flow
is used to model posterior distribution of the parametersof interest characterizing the contribution
1of point source populations as well as diffuse (“smooth”) components to the γ-ray data. The flowtransformation from the base distribution to the posterior is conditioned on learned summaries of the
γ-ray map extracted using a convolutional neural network. The normalizing flow and feature-extractor
neural networks are trained simultaneously using maps simulated from the forward model. Once
trained, samples from the flow can be generated conditioned on a new dataset of interest in order toobtain an estimate of the corresponding parameter posteriors, which can be used to infer physicalquantities of interest such as source-count distributions of modeled PS populations as well as fluxes
associated with the diffuse components.
Here, we showcase a complementary approach that leverages recent developments in simulationbased inference (SBI, also referred to as likelihood-free inference; see, e.g., Ref. [18] for a recent
review) in order to weigh in on the nature of the GCE. In particular, we use conditional density
estimation techniques based on normalizing flows [19, 20] to characterize the contributions of various
modeled components, including “clumpy” PS-like and “smooth” DM-like emission spatially tracing
the GCE, to the γ-ray photon sky at ∼ GeV energies in the Galactic Center region. Rather than using
hand-crafted summary statistics, we employ a graph-based convolutional neural network architecture(previously utilized in Refs. [15, 16]) in order to extract summary statistics from γ-ray maps optimized
for the downstream task of estimating the distribution of parameters characterizing the contributionof modeled components to the GCE. Unlike traditional approaches based on the statistics of photon
counts, this approach lets us capture more of the information contained in a model of the Galactic
Center emission, and in particular implicitly uses the distribution of correlations between pixels as
an additional discriminating handle. As we show in our extended paper [? ] alongside more details
on the analysis pipeline and validation tests on simulated data, this fact makes our method more
resilient to certain systematic uncertainties associated with model misspecification in real Fermi data.
A schematic illustration of our method is presented in Fig. 1.

2

Model and inference

The forward model We use the datasets and spatial templates from Refs. [21, 22] to create
simulated maps of Fermi data in the Galactic Center region. The maps are spatially binned using
the HEALPix [23] pixelization scheme with resolution parameter nside=128, roughly corresponding
to pixel area ∼ 0.5 deg2 . The inner region of the Galactic plane, where the observed emission isespecially difficult to model, is masked at latitudes |b| < 2◦ , and a radial cut r < 25◦ defines our
region of interest (ROI).
2

The simulated maps are a combination of diffuse (alternatively referred to as smooth or Poissonian) and PS contributions. The smooth contributions include (i) the Galactic diffuse foregroundemission [24], (ii) spatially isotropic emission accounting for, e.g., uniform emission from unresolved sources of extragalactic origin, (iii) emission from resolved PSs included in the Fermi
3FGL catalog [25], and (iv) lobe-like emission associated with the Fermi bubbles [26]. Finally,(v) Smooth DM-like emission is modeled using a line-of-sight integral of the (squared) generalized−γ−3+γNavarro-Frenk-White (NFW) [27, 28] profile, ρgNFW (r) ∝ (r/rs ) (1 + r/rs )with innerslope γ = 1.2 motivated by previous GCE analyses [29, 6, 30]. The total smooth component isobtained as a Poisson realization of a linear combination of these spatial templates.
Assuming the locations of individual PSs are not known a-priori, the statistics of multiple PS
populations can be completely specified through (i) their spatial distribution, described by templates
T p discretized over pixels p, (ii) the distribution of expected photon counts S contributed by each PS,
p(S), and (iii) the distribution of the number of PSs for each population. Additionally, the modeled
instrumental point-spread function quantifies the spatial distribution of photon counts sourced byindividual PS around its location due to the finite angular resolution of the instrument. Here, weparameterize the distributions of photon counts S contributed by each PS through a doubly-broken
power law specified by the break locations {Sb,1 , Sb,2 }, spectral indices (slopes) {n1 , n2 , n3 }, andappropriately normalized to unity. Together, we denote these parameters by θPS .
The PS components of the simulated Fermi map are created as follows, practically implemented
using the code package NPTFit-Sim [31]. The total number of PSs to be simulated is drawn asn ∼ Pois(n | npix λ), where npix is the number of pixels in the ROI. The sample of PS angular
positions is drawn from a PDF constructed by linearly interpolating the relevant pixel-wise spatial
template T p ; {rn } ∼ p(r) ∝ T (r). The expected number of photons emitted by each PS, indexed
by i, is drawn by sampling from the mean source-count distribution, S ∼ p (S | θPS ), and scaling
to correct for non-uniform exposure of the satellite. The actual sample of photon counts emitted
by the simulated PSs, {xn }, is taken to be a Poisson realization of this expectation. The procedure
is repeated for each PS population, and the final simulated PS map is constructed by binning thesampled photon positions within the ROI according to the pixelization scheme used. The total map isobtained by combining the simulated diffuse and PS components. The inclusion of PSs in the forward
model introduced a large number of latent variable—the positions and fluxes associated with each
PS—and renders the full likelihood of the model intractable.
Modeled PS populations are often compactly described through the so-called source-count distribution(SCD) dN/dS, which quantifies the differential number density of sources per unit angular area
(more formally d2 N/dSdΩ, although we leave the area dependence implicit) emitting S photons
in expectation. The source-count distribution jointly describes the distribution of photon counts
from individual PSs p(S | θPS ) and their mean per-pixel abundance λ, and is related to these asdN /xdS = λ p(S | θPS )/Ωpix where the the pixel area Ωpix is used to convert the per-pixel source
count to per-area, agnostic to pixel size. We will present our results in terms of the source fluxes
dN/dF , with the conversion S = hiF where hi is the mean exposure in the region considered.
Two PS populations are modeled—(i) those correlated with the GCE, following an NFW profile, andii those tracing the Galactic disk, spatially modeled using a doubly exponential profile.
The forward model is thus specified by a total of 18 parameters—6 for the overall normalizations of
the Poissonian templates, and 6 × 2 parameters modeling the source-count distributions associatedwith GCE-correlated and disk-correlated PS populations {hS PS i, n1 , n2 , n3 , Sb,1 , Sb,2 }. hS PS idenotes the mean per-pixel counts contributed by a given PS population, and parameterizes theiroverall abundance.
Inference with likelihoods based on simplified data representations The 1-point PDF (probability distribution function) framework, first introduced in the context of γ-ray analyses in Ref. [32]and extended in Refs. [8, 9] under the name of non-Poissonian template fitting (NPTF), considers
a simplification of the problem by computing the pixel-wise likelihood assuming each pixel tobe statistically independent (1-point then referring to values over individual, independent spatial
positions in the sky). This significantly reduces the latent space dimensionality by eliminating thepositions of individual PSs as latent variables, localizing them within a pixel and modulating their
expected abundance by the modeled spatial template (e.g., GCE-correlated or disk-correlated in ourcase). Here, we use this method as a comparison point, and sample the posterior associated withparameters of interest with dynesty [33] using the likelihood from NPTFit [22].
3Component intensities

‘1-σ’

10−11

10−11

10−10

10−9 0

1

2

4

5

10

20

5

0

30

0

Component intensities

5
10DM-like [%]15

GCE flux fractions

15

‘1-σ’

Source-count distributions

Fermi data, Baseline
NPTF
PS-like [%]1-ph

F 2 dN/dF [ph cm−2 s−1 deg−2]310Intensity [10−7 ph cm−2 s−1 sr−1]

F [ph cm−2 s−1]

10−10Fermi data, Baseline
Simulation-based inference

π 0 /bremICS

GCE PS
Disk PS

10−12

GCE flux fractions

15

GCE PSGCE DM
Disk PSIso
Bub

PS-like [%]1-ph

F 2 dN/dF [ph cm−2 s−1 deg−2]

Source-count distributions
10−10

10−11

10

5

10−12
10−11

10−10

F [ph cm−2 s−1]

10−9 0

1

2

3

4

5

Intensity [10−7 ph cm−2 s−1 sr−1]

10

20

30

0

0

5
10
DM-like [%]15

Figure 2: Results of the baseline analysis on real Fermi data. (Top row) Analysis using neural
simulation-based inference with normalizing flows, and (bottom row) using the 1-point PDF likelihood
implemented in the non-Poissonian template fitting (NPTF) framework. While moderate preference
for a PS-like origin of the GCE is seen in the case of the NPTF analysis (bottom), the simulation-based
inference analysis attributes a smaller fraction of the GCE to PS-like emission (top).

Extracting representative features from γ-ray maps Rather than relying on hand-crafted datasummaries, a neural network is used to extract representative features sϕ (x) of the data x optimized
for the downstream density estimation task. We use the DeepSphere architecture [34–36] with aconfiguration similar to and inspired by that employed in Ref. [15]. DeepSphere is a graph-based
spherical convolutional neural network architecture tailored to data sampled on a sphere, and inparticular is able to leverage the hierarchical structure of data in the HEALPix representation. This
makes it well-suited for our purposes. The architecture consists of graph convolutional layers which,following a ReLU nonlinearity, coarsens the pixel representation by a factor of 4 with max poolingwhile doubling the number of feature dimensions until a maximum of 256. The output of the final
convolution layer is passed through a fully connected neural network with 1024 hidden units before
outputting 128 summaries.

Simulation-based inference with normalizing flows Simulation-based inference (SBI) refers to aclass of methods for performing inference when the data-generating process does not have a tractable
likelihood. This is the case for the forward model used here, where the presence of a large numberof PSs leads to a large latent space. We approximate the joint posterior over the parameters ofinterest θ given a γ-ray map x through a distribution p̂φ (θ | sϕ (x)) conditioned on summaries sϕ (x)from simulated samples {x}. The conditioned posterior distribution is parameterized by φ anddefined via normalizing flows [19, 20], which are a class of models that provide an efficient way ofmodeling high-dimensional probability distributions. Specifically, we use Masked Autoregressive
Flows (MAFs) [37] to define the flow transformation. We use 8 MAF transformations, each made up
of a 2-layer masked autoregressive neural network [38] with 128 hidden units and tanh activations.
Each transformation is conditioned on summaries sϕ (x) by including these as additional inputs into
the transformation blocks.
Normalizing flows allow for tractable density evaluation, and log p̂φ (θ | sϕ (x)) is used as the training
objective to simultaneously optimize parameters {φ, ϕ} associated with the convolution and flow
neural networks, respectively. 106 samples from the forward model are produced, with 15% of
these held out for validation. The model is trained for up to 30 epochs with early stopping, using abatch size of 256. The AdamW optimizer [39, 40] is used with initial learning rate 10−3 and weightdecay 10−5 , with learning rate decayed through cosine annealing. Experiments were performed on
RTX8000 GPUs at the NYU Greene computing cluster.
4

3

Application to Fermi data

We apply our neural simulation-based inference pipeline to the real Fermi dataset. As a point ofcomparison, we also run the NPTF method on the data using the same spatial templates and priorassumptions as those used in the corresponding SBI analyses. The results of the NPTF analysis
are shown in the bottom panel of Fig. 2. The left column shows the median (solid lines) as well
as middle-68/95% containment (dark/light shaded regions) of the posteriors on the source-countdistributions F 2 dN/dF of GCE-correlated (red) and disk-correlated (blue) PSs, evaluated point-wisein flux F . The dashed grey vertical lines correspond to the flux associated with a single expected
photon count per source (below which Poissonian and PS-like emission is expected to be perfectly
degenerate) and the approximate 1-σ threshold for detecting individual sources (below which thedegeneracy is often observed in practice [41, 24]). The middle column shows the posteriors on
various modeled emission components. The right column shows the joint posterior on the fraction of
DM- and PS-like emission in proportion to the total inferred flux in the ROI.
Consistent with previous 1-point PDF studies using a similar configuration, a significant fractionof the GCE—55.0+8.8−22.9 %–is attributed to PS-like emission. The top panel of Fig. 2 shows results
using the neural simulation-based analysis pipeline introduced in this paper. Although posteriors
for the astrophysical background templates are seen to be broadly consistent with those inferred in
the NPTF analysis, the preference for PSs is somewhat reduced in this case, with 37.9+8.9−19.2 % of the
GCE emission being PS-like.

4

Discussion

We have leveraged recent advances in neural simulation-based inference in order to jointly characterize
a putative DM-like signal and PS population associated with the observed Fermi Galactic Center
Excess. While broadly consistent with results of the traditional method, our method shows a reduced
preference for PS-like emission correlated with the GCE. In our extended work Ref. [? ], we
present additional details of our analysis, including a validation of the pipeline on simulated data
as well as a discussion of the impact of model misspecification within our framework. We show
there that, owing to the fact that it can extract more information from the forward model, ourmethod can be less sensitive to certain forms of model misspecification compared to traditional
approaches. Although a direct comparison is difficult, our results are broadly consistent with andcomplementary to those obtained in Ref. [16], which used a DeepSphere-based architecture which
was, in contrast to our parametric approach, combined with a novel neural network-based nonparameteric approach to infer the counts distributions associated to PS populations using histogramswith modeled uncertainties [42].
As in any Galactic Center γ-ray analysis, we caution of the potential of unknown systematics, suchas mismodeling on the scale of the size of the LAT point-spread function, to bias the results andconclusions of our analysis. Although machine learning-based analyses can utilize more of theinformation encoded in the forward model, and in particular in the present case can take advantage of
pixel-to-pixel correlations, this can also make them more susceptible to specific modeled featurescompared to traditional techniques based on data reduction to hand-crafted data summaries. We leave
a more detailed investigation of the impact of these effects to future work.
Code used for reproducing the results presented in this paper is available at https://github.com/
smsharma/fermi-gce-flows.

Acknowledgments and Disclosure of Funding
We thank Johann Brehmer and Tracy Slatyer for helpful conversations. We are grateful to Florian
List and Nick Rodd for their many helpful comments. SM would like to thank the Center for
Computational Astrophysics at the Flatiron Institute for their hospitality while this work was being
performed. This work was performed in part at the Aspen Center for Physics, which is supported by
National Science Foundation grant PHY-1607611. The participation of SM at the Aspen Center for
Physics was supported by the Simons Foundation. SM is supported by the NSF CAREER grant PHY1554858, NSF grants PHY-1620727 and PHY-1915409, and the Simons Foundation. KC is partially
supported by NSF awards ACI-1450310, OAC-1836650, and OAC-1841471, the NSF grant PHY5

1505463, and the Moore-Sloan Data Science Environment at NYU. This work is supported by the
National Science Foundation under Cooperative Agreement PHY-2019786 (The NSF AI Institute for
Artificial Intelligence and Fundamental Interactions, http://iaifi.org/). This material is based
upon work supported by the U.S. Department of Energy, Office of Science, Office of High Energy
Physics of U.S. Department of Energy under grant Contract Number DE-SC0012567. We thankthe Fermi-LAT Collaboration for making publicly available the γ-ray data used in this work. This
work made use of the NYU IT High Performance Computing resources, services, and staff expertise.
This research has made use of NASA’s Astrophysics Data System. This research made use of theastropy [43, 44], dynesty [33], getdist [45], IPython [46], Jupyter [47], matplotlib [48],MLflow [49], nflows [50], NPTFit [22], NPTFit-Sim [31], NumPy [51], pandas [52], PyGSP [53],Pyro [54], PyTorch [55], PyTorch Geometric [56], PyTorch Lightning [57], seaborn [58],sbi [59], scikit-learn [60], SciPy [61], and tqdm [62] software packages. We acknowledge
the use of data products and templates from the code repository associated with Ref. [15].1 We
acknowledge use of the DeepSphere graph-convolutional layer from the code repository associatedwith Ref. [63].2

References[1] L. Goodenough and D. 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].
[2] D. Hooper and L. Goodenough, “Dark Matter Annihilation in The Galactic Center As Seen by the Fermi
Gamma Ray Space Telescope,” Phys. Lett. B 697, 412 (2011), arXiv:1010.2752 [hep-ph].
[3] A. Boyarsky, D. Malyshev, and O. Ruchayskiy, “A comment on the emission from the Galactic Center as
seen by the Fermi telescope,” Phys. Lett. B 705, 165 (2011), arXiv:1012.5839 [hep-ph].
[4] D. Hooper and T. Linden, “On The Origin Of The Gamma Rays From The Galactic Center,” Phys. Rev. D
84, 123005 (2011), arXiv:1110.0006 [astro-ph.HE].
[5] W. B. Atwood et al. (Fermi-LAT), “The Large Area Telescope on the Fermi Gamma-ray Space Telescope
Mission,” Astrophys. J. 697, 1071 (2009), arXiv:0902.1089 [astro-ph.IM].
[6] T. Daylan, D. P. Finkbeiner, D. Hooper, T. Linden, S. K. N. Portillo, N. L. Rodd, and T. R. Slatyer, “Thecharacterization of the gamma-ray signal from the central Milky Way: A case for annihilating dark matter,”Phys. Dark Univ. 12, 1 (2016), arXiv:1402.6703 [astro-ph.HE].
[7] F. Calore, I. Cholis, and C. Weniger, “Background Model Systematics for the Fermi GeV Excess,” JCAP
03, 038 (2015), arXiv:1409.0042 [astro-ph.CO].
[8] S. K. Lee, M. Lisanti, and B. R. Safdi, “Distinguishing Dark Matter from Unresolved Point Sources in the
Inner Galaxy with Photon Statistics,” JCAP 05, 056 (2015), arXiv:1412.6099 [astro-ph.CO].
[9] S. K. Lee, M. Lisanti, B. R. Safdi, T. R. Slatyer, and W. Xue, “Evidence for Unresolved γ-Ray Point
Sources in the Inner Galaxy,” Phys. Rev. Lett. 116, 051103 (2016), arXiv:1506.05124 [astro-ph.HE].
[10] R. Bartels, S. Krishnamurthy, and C. 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].
[11] B. Balaji, I. Cholis, P. J. Fox, and S. D. McDermott, “Analyzing the Gamma-Ray Sky with Wavelets,” Phys.
Rev. D 98, 043009 (2018), arXiv:1803.01952 [astro-ph.HE].
[12] S. D. McDermott, P. J. Fox, I. Cholis, and S. K. Lee, “Wavelet-Based Techniques for the Gamma-Ray Sky,”
JCAP 07, 045 (2016), arXiv:1512.00012 [astro-ph.HE].
[13] 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, 231103 (2020), arXiv:1911.12369 [astro-ph.HE].
[14] S. Caron et al., “Identification of point sources in gamma rays using U-shaped convolutional neuralnetworks and a data challenge,” (2021), arXiv:2103.11068 [astro-ph.HE].
[15] F. List, N. L. Rodd, G. F. Lewis, and I. Bhat, “The GCE in a New Light: Disentangling the γ-ray Sky with
Bayesian Graph Convolutional Neural Networks,” Phys. Rev. Lett. 125, 241102 (2020), arXiv:2006.12504[astro-ph.HE].
12

https://github.com/FloList/GCE_NN
https://github.com/deepsphere/deepsphere-pytorch

6[16] F. List, N. L. Rodd, and G. F. Lewis, “Dim but not entirely dark: Extracting the Galactic Center Excess’
source-count distribution with neural nets,” (2021), arXiv:2107.09070 [astro-ph.HE].
[17] S. Caron, G. A. Gómez-Vargas, L. Hendriks, and R. Ruiz de Austri, “Analyzing γ-rays of the Galactic
Center with Deep Learning,” JCAP 05, 058 (2018), arXiv:1708.06706 [astro-ph.HE].
[18] K. Cranmer, J. Brehmer, and G. Louppe, “The frontier of simulation-based inference,” Proc. Nat. Acad.
Sci. 117, 30055 (2020), arXiv:1911.01429 [stat.ML].
[19] G. Papamakarios, E. Nalisnick, D. J. Rezende, S. Mohamed, and B. Lakshminarayanan, “Normalizing
flows for probabilistic modeling and inference,” (2019), arXiv:1912.02762 [cs.LG].
[20] D. J. Rezende and S. Mohamed, “Variational inference with normalizing flows,” in Proceedings of the
32nd International Conference on Machine Learning, ICML 2015, Lille, France, 6-11 July 2015, JMLR
Workshop and Conference Proceedings, Vol. 37, edited by F. R. Bach and D. M. Blei (2015) pp. 1530–1538.
[21] S. Mishra-Sharma, N. L. Rodd, and B. R. Safdi, “Supplementary material for NPTFit,” (2016).
[22] S. Mishra-Sharma, N. L. Rodd, and B. R. Safdi, “NPTFit: A code package for Non-Poissonian Template
Fitting,” Astron. J. 153, 253 (2017), arXiv:1612.03173 [astro-ph.HE].
[23] K. M. Gorski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelman,“HEALPix - A Framework for high resolution discretization, and fast analysis of data distributed on thesphere,” Astrophys. J. 622, 759 (2005), arXiv:astro-ph/0409513.
[24] M. Buschmann, N. L. Rodd, B. R. Safdi, L. J. Chang, S. Mishra-Sharma, M. Lisanti, and O. 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].
[25] F. Acero et al. (Fermi-LAT), “Fermi Large Area Telescope Third Source Catalog,” Astrophys. J. Suppl.
218, 23 (2015), arXiv:1501.02003 [astro-ph.HE].
[26] M. Su, T. R. Slatyer, and D. P. Finkbeiner, “Giant Gamma-ray Bubbles from Fermi-LAT: AGN Activity or
Bipolar Galactic Wind?” Astrophys. J. 724, 1044 (2010), arXiv:1005.5480 [astro-ph.HE].
[27] J. F. Navarro, C. S. Frenk, and S. D. M. White, “The Structure of cold dark matter halos,” Astrophys. J.
462, 563 (1996), arXiv:astro-ph/9508025 [astro-ph].
[28] J. F. Navarro, C. S. Frenk, and S. D. White, “A Universal density profile from hierarchical clustering,”
Astrophys.J. 490, 493 (1997), arXiv:astro-ph/9611107 [astro-ph].
[29] C. Gordon and O. 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].
[30] B. Zhou, Y.-F. Liang, X. Huang, X. Li, Y.-Z. Fan, L. Feng, and J. 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].
[31] N. L. Rodd and M. W. Toomey, NPTFit-Sim (2017).
[32] D. Malyshev and D. W. Hogg, “Statistics of gamma-ray point sources below the Fermi detection limit,”
Astrophys. J. 738, 181 (2011), arXiv:1104.0010 [astro-ph.CO].
[33] J. S. Speagle, “dynesty: a dynamic nested sampling package for estimating bayesian posteriors andevidences,” Monthly Notices of the Royal Astronomical Society 493, 3132 (2020).
[34] M. Defferrard, M. Milani, F. Gusset, and N. Perraudin, “Deepsphere: a graph-based spherical cnn,” arXiv
preprint arXiv:2012.15000 (2020).
[35] N. Perraudin, M. Defferrard, T. Kacprzak, and R. Sgier, “DeepSphere: Efficient spherical Convolutional
Neural Network with HEALPix sampling for cosmological applications,” Astron. Comput. 27, 130 (2019),
arXiv:1810.12186 [astro-ph.CO].
[36] M. Defferrard, N. Perraudin, T. Kacprzak, and R. Sgier, “DeepSphere: towards an equivariant graphbased spherical CNN,” in ICLR Workshop on Representation Learning on Graphs and Manifolds (2019)arXiv:1904.05146 [cs.LG].

7

[37] G. Papamakarios, T. Pavlakou, and I. Murray, “Masked autoregressive flow for density estimation,” inProceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17(Curran Associates Inc., Red Hook, NY, USA, 2017) pp. 2335–2344.
[38] M. Germain, K. Gregor, I. Murray, and H. Larochelle, “Made: Masked autoencoder for distribution
estimation,” in International Conference on Machine Learning (PMLR, 2015) pp. 881–889.
[39] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in 3rd International Conference on
Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings,edited by Y. Bengio and Y. LeCun (2015).
[40] I. Loshchilov and F. Hutter, “Decoupled weight decay regularization,” in 7th International Conference on
Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019 (2019).
[41] 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, 023014 (2020), arXiv:1908.10874 [astro-ph.CO].
[42] F. List, “The Earth Mover’s Pinball Loss: Quantiles for Histogram-Valued Regression,” in Proceedings ofthe 38th International Conference on Machine Learning (2021).
[43] A. M. Price-Whelan et al., “The Astropy Project: Building an Open-science Project and Status of the v2.0
Core Package,” Astron. J. 156, 123 (2018), arXiv:1801.02634.
[44] T. P. Robitaille et al. (Astropy), “Astropy: A Community Python Package for Astronomy,” Astron. Astrophys.
558, A33 (2013), arXiv:1307.6212 [astro-ph.IM].
[45] A. Lewis, “GetDist: a Python package for analysing Monte Carlo samples,” (2019), arXiv:1910.13970[astro-ph.IM].
[46] F. Perez and B. E. Granger, “IPython: A System for Interactive Scientific Computing,” Computing in
Science and Engineering 9, 21 (2007).
[47] T. Kluyver et al., “Jupyter notebooks - a publishing format for reproducible computational workflows,” in
ELPUB (2016).
[48] J. D. Hunter, “Matplotlib: A 2d graphics environment,” Computing In Science & Engineering 9, 90 (2007).
[49] A. Chen et al., “Developments in MLflow: A System to Accelerate the Machine Learning Lifecycle,” inProceedings of the Fourth International Workshop on Data Management for End-to-End Machine Learning,
DEEM’20 (Association for Computing Machinery, New York, NY, USA, 2020).
[50] C. Durkan, A. Bekasov, I. Murray, and G. Papamakarios, “nflows: normalizing flows in PyTorch,” (2020).
[51] C. R. Harris et al., “Array programming with NumPy,” Nature 585, 357 (2020).
[52] W. McKinney, “Data structures for statistical computing in python,” in Proceedings of the 9th Python in
Science Conference, edited by S. van der Walt and J. Millman (2010) pp. 51 – 56.
[53] M. Defferrard, L. Martin, R. Pena, and N. Perraudin, “PyGSP: Graph Signal Processing in Python,”(2017).
[54] E. Bingham et al., “Pyro: Deep universal probabilistic programming,” J. Mach. Learn. Res. 20, 28:1
(2019).
[55] A. Paszke et al., in Advances in Neural Information Processing Systems 32, edited by H. Wallach,
H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019) pp.
8024–8035.
[56] M. Fey and J. E. Lenssen, “Fast graph representation learning with PyTorch Geometric,” in ICLR Workshop
on Representation Learning on Graphs and Manifolds (2019) arXiv:1903.02428 [cs.LG].
[57] W. Falcon et al., “Pytorchlightning/pytorch-lightning: 0.7.6 release,” (2020).
[58] M. Waskom et al., “mwaskom/seaborn: v0.8.1 (september 2017),” (2017).
[59] A. Tejero-Cantero et al., “sbi: A toolkit for simulation-based inference,” Journal of Open Source Software
5, 2505 (2020).

8

[60] F. Pedregosa et al., “Scikit-learn: Machine learning in python,” Journal of Machine Learning Research 12,2825 (2011).
[61] P. Virtanen et al., “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” Nature
Methods (2020), https://doi.org/10.1038/s41592-019-0686-2.
[62] C. da Costa-Luis et al., “tqdm: A fast, Extensible Progress Bar for Python and CLI,” (2021).
[63] M. Defferrard, M. Milani, F. Gusset, and N. Perraudin, “Deepsphere: a graph-based spherical CNN,”in 8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April
26-30, 2020 (2020).

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper’s
contributions and scope? [Yes](b) Did you describe the limitations of your work? [Yes] See Sec. 4(c) Did you discuss any potential negative societal impacts of your work? [N/A] Potential
negative societal impacts were considered, and we believe this work does not present
any issues in this regard.
(d) Have you read the ethics review guidelines and ensured that your paper conforms tothem? [Yes]2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [N/A] No theoreticalresults were obtained in this work.
(b) Did you include complete proofs of all theoretical results? [N/A]
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)? [Yes] The code
repository associated with this paper and needed to reproduce all the results is linked
in Sec. 4.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they
were chosen)? [Yes] These are described in Sec. 3.
(c) Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)? [Yes] Uncertainties are reported, as we obtain the parameter
posterior distributions.
(d) Did you include the total amount of compute and the type of resources used (e.g., typeof GPUs, internal cluster, or cloud provider)? [Yes] Relevant details are provided inSec. 2.
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes] All code used for thisproject is cited in the Acknowledgments section.
(b) Did you mention the license of the assets? [N/A] Licenses are mentioned in the links
associated with individual code packages.
(c) Did you include any new assets either in the supplemental material or as a URL? [N/A]
No new assets (excluding the code used to reproduced the experiments) were producedin this work.
(d) Did you discuss whether and how consent was obtained from people whose data you’reusing/curating? [N/A](e) Did you discuss whether the data you are using/curating contains personally identifiable
information or offensive content? [N/A] No personal information is included in theassets utilized in this paper.
5. If you used crowdsourcing or conducted research with human subjects...
9

(a) Did you include the full text of instructions given to participants and screenshots, if
applicable? [N/A](b) Did you describe any potential participant risks, with links to Institutional Review
Board (IRB) approvals, if applicable? [N/A](c) Did you include the estimated hourly wage paid to participants and the total amountspent on participant compensation? [N/A]10