  The nature of the Fermi gamma-ray Galactic Center Excess (GCE) has remained a
persistent mystery for over a decade. Although the excess is broadly compatible
with emission expected due to dark matter annihilation, an explanation in terms
of a population of unresolved astrophysical point sources e.g., millisecond
pulsars, remains viable. The effort to uncover the origin of the GCE is
hampered in particular by an incomplete understanding of diffuse emission of
Galactic origin. This can lead to spurious features that make it difficult to
robustly differentiate smooth emission, as expected for a dark matter origin,
from more "clumpy" emission expected for a population of relatively bright,
unresolved point sources. We use recent advancements in the field of
simulation-based inference, in particular density estimation techniques using
normalizing flows, in order to characterize the contribution of modeled
components, including unresolved point source populations, to the GCE. Compared
to traditional techniques based on the statistical distribution of photon
counts, our machine learning-based method is able to utilize more of 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 gamma-ray map. This makes the method
demonstrably more resilient to certain forms of model misspecification. On
application to Fermi data, the method generically attributes a smaller fraction
of the GCE flux to unresolved point sources when compared to traditional
approaches. We nevertheless infer such a contribution to make up a
non-negligible fraction of the GCE across all analysis variations considered,
with at least $38^{+9}_{-19}\%$ of the excess attributed to unresolved point
sources in our baseline analysis.
