Dark matter or millisecond pulsars? A deep learning-based analysis of the Fermi Galactic Centre Excess

The γ -ray Galactic Centre Excess (GCE) remains one of the few observed high-energy signals for which a dark matter (DM) origin is a plausible explanation. We present a deep learning-based analysis of the γ -ray sky in the Galactic Centre region, carefully accounting for the mathematical degeneracy between faint point-sources (PSs) such as millisecond pulsars (MSPs) and DM-like Poisson emission. Using recent models for the Galactic foregrounds, we find that relatively few bright PSs just below Fermi ’s detection threshold seem unlikely to explain the GCE, although we continue to find evidence for PSs. Looking ahead, further improvements in the modelling of the γ -ray sky will be crucial for distinguishing between a DM-like and point-like morphology of the signal.


Introduction
The Galactic Centre (GC) constitutes a prime target for indirect dark matter (DM) detection given its relative closeness to Earth and its high DM density (e.g. [1,2]). Shortly after the launch of the Fermi satellite in 2008, an excess of photon counts from the vicinity of the GC was reported, for which DM annihilation was proposed as a possible explanation [3,4]. In the past decade, this excess -known as the Galactic Centre Excess (GCE) -has been thoroughly scrutinised by numerous researchers, and its energy spectrum, spatial morphology, as well as luminosity seem to be largely consistent with an annihilation signature of weakly interacting massive particles (WIMPs), see e.g. [5][6][7][8][9][10][11][12][13][14]. However, astrophysical point-sources (PSs), each of them too faint to be individually resolved, provide an alternative explanation for the GCE (e.g. [15][16][17][18][19][20][21]). The arguments in favour of a PS explanation, most prominently a putative population of millisecond pulsars (MSPs), can generally be subdivided into two categories. First, there are claims that template fits exhibit a preference for the GCE to follow the stellar bulge rather than a squared Navarro-Frenk-White halo profile [22] as would be expected for DM annihilation [23][24][25] and second, some studies found a granular point-like GCE to be preferred over a smooth DM-like GCE morphology [26,27]. Regarding the preference for the stellar bulge over a halo-like shape of the GCE, no final conclusions have been reached yet; for example, the recent studies in Refs. [13,14] found good agreement with a DM-like distribution. For investigating the point-like vs. smooth (Poissonian) nature of the GCEindicative of a PS or DM origin of the GCE, respectively -two methods have traditionally been used: the Non-Poissonian Template Fit (NPTF) [26,28] (see also [29]) and wavelet analyses [27,30,31]. While early analyses using these methods obtained a strong preference for a point-like GCE [26,27], doubt has been cast on the robustness of these results in Refs. [32][33][34], which demonstrated that a spurious preference for a point-like GCE may arise with the NPTF as an artefact of mismodelling the signal and backgrounds (see also Refs. [35] and [36] for discussions about the robustness of the NPTF). Further, a reanalysis of the Fermi data with the wavelet method by Ref. [37] identified a GCE consistent with Poisson emission when using an updated PS mask for the known bright sources [38].
In these proceedings, we will provide a brief synopsis of our work in Refs. [39,40], where we demonstrated that convolutional neural networks (CNNs) [41] are a powerful alternative analysis method for the γ-ray sky, which appear robust to the systematic issues identified for the NPTF. In §2, we will introduce our analysis framework, followed by a brief summary of our main results in §3 and concluding remarks in §4.

Method
We model the γ-ray flux from the inner Galaxy region as a superposition of different emission components, each of which is described by an individual spatial template: (1) diffuse emission from pion decays and bremsstrahlung produced by the interaction of cosmic rays with the interstellar gas (π 0 + BS), (2) diffuse emission due to inverse Compton (IC) up-scattering of interstellar medium photons by cosmic-ray electrons, (3) isotropic emission originating from extragalactic sources, (4) the Fermi bubbles [42], (5) a PS population in the Galactic Disk, and (6) the GCE, which we model with a generalised NFW-squared profile with exponent γ = 1.2. For the disk PSs and the potentially point-like GCE, the exact locations of the individual subthreshold sources are unknown, for which reason the resulting photon counts follow non-Poissonian statistics (see e.g. [26] for more details). In Refs. [39,40], we considered the inner θ ≤ 25 • around the GC and masked latitudes |b| ≤ 2 • as well as the resolved PSs in the 3FGL catalogue [43]. Our full data selection criteria are detailed in those works.

034.2
For our analysis of the Fermi data, we employ CNNs, which have become vastly popular across the sciences due to their ability to detect patterns in complex image data such as template edges and gradients for our task of disentangling γ-ray emission. While CNNs are arguably less interpretable than traditional analysis methods, they also possess notable advantages: first, CNNs assess entire regions of different sizes in γ-ray maps -in contrast to the product-over-pixel likelihood maximised by traditional template fitting methods, which may produce overly narrow posteriors because the counts in neighbouring pixels are in fact not independent of each other due to the point-spread function (PSF), which is what a product likelihood implicitly assumes (see Ref. [44]). This conceptual difference leads to a different behaviour of the methods in the presence of modelling deficiencies, see Ref. [40]. Second, as a representative of supervised machine learning methods, CNNs learn to infer the underlying relation between inputs (given by photon-count maps) and outputs from a labelled training dataset. Hence, the CNN can learn the effect of the PSF without being restricted by the need for explicitly prescribing an analytic approximation of the true likelihood. We take the label of each photon-count map to be the flux fraction of each template as well as the source-count distributions (SCDs) F d N d log 10 F of the PS templates (disk and GCE), which express the expected flux F of the d N PSs within each infinitesimal logarithmic flux bin d log 10 F .
We generate the training data by drawing from Poissonian distributions and using a modified version of NPTFit-Sim [45] for the Poissonian and PS templates, respectively. Importantly, we model all emission components that are potentially point-like as non-Poissonian. The rationale is as follows: Poisson emission can be viewed as the limit of N → ∞ PSs that each contribute a negligible flux F → 0 to the map, where the limit is taken such that N F = µ = const. Thus, by choosing the priors for the SCDs in such a way that narrow SCDs peaking at very low fluxes corresponding to ≪ 1 count per PS occur sufficiently often, many training maps will contain PS emission that is virtually indistinguishable from Poissonian emission. Using this approach, the (technically artificial) separation between Poissonian and pointlike emission can be circumvented, acknowledging the fact that sufficiently faint PS emission is statistically degenerate with Poisson emission. 1 Our analysis framework consists of two separate CNNs. The first CNN learns to estimate the flux fraction belonging to each template, disregarding whether the emission is Poissonian or point-like. Then, a second CNN is trained to predict the SCDs of the PS templates. For this purpose, we discretise the SCDs as histograms and utilise the Earth Mover's Pinball Loss function [46], which quantifies the uncertainties for the SCD histograms in terms of quantiles.
Our CNN architecture is built on the graph-based DeepSphere framework [47,48].
Finally, by comparing the predicted SCD for the Fermi GCE with our CNN predictions for simulated best-fit maps that contain an entirely Poissonian GCE, constraints on the Poissonian flux fraction of the Fermi GCE can be derived: GCE flux that the CNN attributes to bright PSs with high confidence is unlikely to have been produced by a Poissonian emission mechanism that follows the smooth NFW-squared template, whereas flux at the faint end of the estimated SCD cannot confidently be attributed to either dim PSs or Poissonian emission given their inherent degeneracy. For the full details of our methodology, we refer to Refs. [39,40].

Results
Here, we present the main results derived from the Fermi data. Extensive experiments validating the performance of our CNN-based framework on simulated maps, as well as various  [23,49]) and has been shown to give a significantly better (albeit still suboptimal) fit to the Fermi data than p6v11. In particular, the overly hard IC component of p6v11 at energies ≥ 10 GeVwhich is fused with the diffuse foregrounds from pion decay and bremsstrahlung into a single template -may lead to over-subtraction of other emission components [11,36]. The results presented herein rely on the same priors for the data generation as in Ref. [39], but we binned the photon counts into HEALPix [50] pixels at resolution N side = 128 (instead of N side = 256 in that work). Notably, the flux assigned to the GCE increases by almost a factor of two when replacing p6v11 with Model O, and much more flux is attributed to disk PSs. For completeness, we remark that the flux fraction estimates obtained with the NPTFit implementation [28] of the NPTF are very similar to those of our CNNs for both diffuse models (see Ref. [39]), demonstrating that the modelling uncertainties outweigh the differences between the results of the two methods. Figure 1 shows our SCD estimates for the GCE, again for both diffuse models. Here, the difference is particularly striking: consistent with the original NPTFit-based analysis in Ref. [26], we identify a GCE PS population located not far below the Fermi detection threshold when using the diffuse model p6v11, implying that some hundreds to a few thousand relatively bright PSs would suffice to explain the GCE. However, with the more recent Model O, the SCD shifts to much fainter fluxes, peaking at a value corresponding to 3 − 4 counts per PS with the data selection cuts in this present analysis. Importantly, for such a SCD, we cannot exclude that its low-flux end in fact represents Poissonian emission in view of the above-mentioned Poisson vs. dim PS degeneracy. Yet, also in this case the GCE was not found to be consistent with 100% Poissonian flux, and a PS contribution was required (see Ref. [40]).

Conclusions and outlook
In these conference proceedings, we have provided a bird's eye view of our machine learningbased framework for the analysis of the γ-ray sky presented in Refs. [39,40]. While modelling uncertainties still preclude a definite answer as to whether the GCE is caused by DM annihilation or by PSs such as MSPs, an important avenue is the development of novel analysis methods that -once the modelling of the γ-ray sky becomes sufficiently accurate -are able to accurately and reliably infer the GCE's flux, SCD, and energy spectrum. Importantly, since different analysis methods exhibit different systematics in the presence of modelling deficiencies (which can of course never be completely eradicated), taking the results of multiple methods in aggregate provides a more complete and in particular robust picture of the GCE's nature. A valuable extension of our method would be to consistently incorporate modelling uncertainties into our analysis by endowing templates with additional degrees of freedom (as done by e.g. Ref. [51] for a likelihood-based approach). In this spirit, Refs. [52,53] recently explored the use of Gaussian Processes for making the diffuse emission model more flexible in the context of Normalising Flows.
Another possible refinement is to harness the spectral information alongside the photoncount statistics. So far, methods that attempt to distinguish DM from PSs based on photoncount statistics such as the NPTF or deep learning methods have analysed the Fermi data integrated over a certain energy range (e.g. 2 − 20 GeV), thereby discarding the spectral information contained in the data (which, of course, has been the subject of other templatebased analyses of the GCE, however, see e.g. Ref. [11]). This direction seems particularly promising because it would not only provide more nuanced information, but would also enable (1) the use of state-of-the-art templates that are tailored to specific energy ranges (e.g. [54]) and (2) a more accurate treatment of the Fermi PSF, which improves at higher energies. 2 With further progress on the modelling of the γ-ray sky, improvements in the analysis methods, and possibly the (non-)detection of the putative Galactic Bulge MSPs with upcoming radio telescopes, the 2020's will hopefully be the decade in which the GCE mystery is resolved.