A Methods

Evidence for Unresolved Gamma-Ray Point Sources in the Inner Galaxy


We present a new method to characterize unresolved point sources (PSs), generalizing traditional template fits to account for non-Poissonian photon statistics. We apply this method to Fermi Large Area Telescope gamma-ray data to characterize PS populations at high latitudes and in the Inner Galaxy. We find that PSs (resolved and unresolved) account for 50% of the total extragalactic gamma-ray background in the energy range 1.9 to 11.9 GeV. Within 10 of the Galactic Center with , we find that 5–10% of the flux can be accounted for by a population of unresolved PSs, distributed consistently with the observed GeV gamma-ray excess in this region. The excess is fully absorbed by such a population, in preference to dark-matter annihilation. The inferred source population is dominated by near-threshold sources, which may be detectable in future searches.


Dark-matter (DM) annihilation in the Galactic halo can contribute to the flux of high-energy gamma rays detected by experiments such as the Fermi Large Area Telescope [1]. Currently, an excess of GeV gamma rays has been observed by Fermi near the Galactic Center (GC) [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]. The signal extends 10 off the plane, is approximately spherically symmetric, and has an intensity profile that falls as with  [12, 14]. The morphology and energy spectrum of the signal is consistent with DM annihilation. There is some possible tension between the DM interpretation and other searches, especially in dwarf galaxies [17]; alternate explanations include a new population of millisecond pulsars (MSPs) [18, 19, 20, 11, 21, 22, 23, 24, 25, 26] or cosmic-ray injection [27, 28].

This Letter addresses the potential contribution of unresolved point sources (PSs) to the excess through the use of a new statistical method, called a non-Poissonian template fit (NPTF). Our approach is model-independent, in that we remain agnostic about the nature of the PSs. To verify the method, we use it to characterize unresolved gamma-ray PSs at high Galactic latitudes. These findings represent one of the most precise measurements of the contribution of PSs to the extragalactic gamma-ray background (EGB) and have important implications for characterizing its source components.

The main focus of this Letter is to use the NPTF to search for a population of unresolved gamma-ray PSs in the Inner Galaxy (IG) with a morphology consistent with that of the excess. We find that the NPTF strongly prefers a PS origin for the excess over a DM-like (smooth diffuse) origin. The Supplementary Material provides further details on the method, as well as additional cross-checks that support these conclusions.

This study analyzes the Extended Pass 7 Reprocessed Fermi data from August 4, 2008 to December 5, 2013 made available by [29]. A HEALPix [30] pixelization of the data with is used, corresponding to pixels 0.5 to a side. We emphasize that our study focuses on data in a single energy bin from 1.893–11.943 GeV and does not rely on or extract spectral information for the excess. The choice of this energy range keeps the signal-to-background ratio in the region of interest (ROI) high, maintains a sufficiently large number of photons over the full sky, and keeps the point-spread function relatively small and energy-independent.

The analysis utilizes the photon-count probability distribution in each pixel. In general, a given model for the gamma-ray flux, with parameters , predicts a probability of observing  photons in a pixel . Several source components, each modeled by a spatial template, can contribute photons in a pixel. To date, analyses using templates have assumed Poisson statistics for the photon-count distribution—specifically, that is the Poisson probability to draw counts with mean given by the sum of the template components in pixel .

To account for unresolved PSs, the standard template-fitting procedure must be generalized to include non-Poissonian photon counts. In the NPTF procedure, depends on a potentially pixel-dependent PS source-count function . The source-count function determines the average number of PSs within pixel that contribute photon flux between and . In this work, the source-count function is assumed to follow a broken power law, , with pixel-dependent normalization and indices () above (below) the break that are constant between pixels. For isotropically distributed PSs, is constant between pixels. To model a population of PSs that mimics a DM annihilation signal, must instead follow the DM annihilation template. Semi-analytic methods for calculating the with a broken power law source-count function were developed in  [31, 32].

Figure 1: The source-count function for high-latitude point sources, derived from applying non-Poissonian template fits to data with 3FGL sources [33] unmasked (green band) and masked (orange band). The colored bands indicate 68% confidence intervals, which are computed pointwise in from the posteriors for the source-count–function parameters, while the solid and dashed black lines show the median source-count functions. The black points show the source-count function of the detected point sources in the 3FGL catalog. The vertical error bars indicate 68% confidence intervals; the horizontal bars denote bins in .

We include templates for up to seven different components in the NPTF analysis: (1) diffuse background, assuming the Fermi p6v11 diffuse model; (2) Fermi bubbles, assuming uniform emission within the bubbles [34]; (3) isotropic background; (4) annihilation of Navarro-Frenk-White (NFW) [35, 36]-distributed DM, assuming no substructure; (5) isotropic PSs; (6) disk-correlated PSs, and (7) NFW-distributed PSs.2 Templates (1) through (4) are specified by a single normalization parameter each. Templates (4) and (7) assume a generalized NFW distribution with inner slope . Template (6) corresponds to a doubly exponential thin-disk source distribution with scale height 0.3 kpc and radius 5 kpc. The PS templates each have four parameters describing their respective source-count functions.

Bayesian methods (implemented with MultiNest [37, 38]) are used to extract posterior distributions for the model parameters. The prior distributions of all parameters are flat, except for those of the DM and PS normalization factors, which are log flat. Unless otherwise stated, the prior ranges of all parameters are sufficiently large so that the posterior distribution is well-converged.

We begin by applying the NPTF to data at high Galactic latitudes (). Fermi’s third source catalog (3FGL) [33] identifies 1307 gamma-ray PSs in this region of the sky, with 55% associated with Active Galactic Nuclei and 24% associated with pulsars, supernova remnants and other known gamma-ray sources. The remaining 21% are unassociated. Figure 1 shows the source-count function in terms of the flux of the 3FGL PSs in our energy bin (black points). The observed source-count function is suppressed below  photons, where it is hard to detect PSs over the diffuse background with the current exposure.

The NPTF is performed in this high-latitude region, including templates for the diffuse background, Fermi bubbles, isotropic emission, and isotropic PSs. The best-fit source-count function values are given in Tab. 1.3 The pointwise % confidence interval for the source-count function is shown in Fig. 1, shaded green. The source-count function matches the 3FGL data well above  photons.

The best-fit intensities obtained from the NPTF can be compared to those obtained from a standard template fit that neglects PSs. The diffuse-background and Fermi-bubbles intensities (averaged over the ROI) are consistent between both procedures. When the PS template is included, the isotropic-background intensity is photons/cm/s/sr and the isotropic PS intensity is photons/cm/s/sr. With no PS template, the isotropic-background intensity is over twice as high, photons/cm/s/sr. Thus, the PS intensity is absorbed by the isotropic-background template in the standard procedure.

ROI Template 3FGL Bayes Factor Bayes Factor NFW DM
[photons/cm/s] (Data) (Simulation) (95% confidence)
HL Iso. PS unmasked
IG \makecellNFW PS
Disk PS
unmasked 10 10 %
IG \makecellNFW PS
Disk PS
masked 10 10 %
Table 1: Best-fit values (16, 50, and 84 percentiles) for the source-count functions associated with the PS templates in the High Latitude (HL) and Inner Galaxy (IG) ROIs. The source-count function is fit by a broken power-law, where is the slope above (below) the break in , given by . The source-count function for the isotropic PS component in the IG is not included, as its flux fraction is subdominant. Depending on the analysis, the Fermi 3FGL sources may either be masked or unmasked. Where appropriate, we provide the Bayes factor in preference for including the NFW PS component, in both the real data and in simulations, as well as the constraint on the flux fraction (calculated as in Fig. 2) attributed to NFW DM.

The averaged intensity of the observed 3FGL PSs is 9.32 photons/cm/s/sr at high latitudes. Using the result of the NPTF described above and neglecting systematic uncertainties in modeling the 3FGL PSs, we predict that the intensity of unresolved PSs is photons/cm/s/sr. This can be checked explicitly by repeating the NPTF with all 3FGL PSs masked (at a 1 radius). The results of this fit are given in Tab. 1 and illustrated by the orange band in Fig. 1. The source-count function for the unresolved PSs agrees with that computed from the unmasked sky at low flux. This suggests that the NPTF is sensitive to contributions from unresolved PSs below Fermi’s detection threshold. The intensity of the isotropic background is photons/cm/s/sr, which agrees with that from the 3FGL-unmasked NPTF, within uncertainties. The intensity of the isotropic PSs is photons/cm/s/sr, which is slightly lower than the value inferred from the 3FGL-unmasked NPTF.

corresponds to the intensity of the isotropic gamma-ray background (IGRB), while gives the intensity of the EGB. While the PS template does absorb some contribution from Galactic PSs, extragalactic PSs are expected to dominate. From the 3FGL-unmasked NPTF, we infer that % of the EGB in this energy range is associated with both resolved and unresolved PS emission; from the 3FGL-masked NPTF and using the intensities of the 3FGL PSs, we find that % of the EGB is due to PS emission. These estimates appear to be consistent with those in [39, 31], though a direct comparison is made difficult by the fact that these analyses cover a different energy range and only use the first 11 months of Fermi data. Our estimates for agree with the most recently published results from Fermi [40].

Figure 2: (Left) Best-fit source-count functions within of the GC and , with the 3FGL sources unmasked. The median and 68% confidence intervals are shown for each of the following PS components: NFW (dashed, orange), thin-disk (solid, blue), and isotropic (dotted, green). The number of observed 3FGL sources in each bin is indicated. The normalization for the diffuse emission in the fit is consistent with that at high latitudes, as desired. (Right) Posteriors for the flux fraction within of the GC with arising from the separate PS components, with 3FGL sources unmasked. The inset shows the result of removing the NFW PS template from the fit. Dashed vertical lines indicate the , , and percentiles.

Figure 3: Same as Fig. 2, except with 3FGL sources masked.

Next, we use the NPTF procedure to determine the fraction of flux from unresolved PSs in the IG. These analyses include templates for the diffuse background, the Fermi bubbles, isotropic background, and NFW-distributed DM, in addition to isotropic, disk-correlated, and NFW-distributed PSs. While the prior ranges for the isotropic, isotropic PS, Fermi bubbles, and diffuse background template parameters are not constrained by the high-latitude fit, restricting these parameters to their high-latitude values does not significantly affect the results.4

The ROI consists of all pixels within of the GC with , masking out the plane. As above, we perform two analyses, one on the full ROI and another with all 3FGL PSs masked. For both cases, the source-count functions and flux fractions are quoted with respect to the region within of the GC and , with no PSs masked. The source-count function of the Galactic and unassociated 3FGL PSs in the IG is given by the black points in the left panel of Fig. 2, with the number of PSs in each bin indicated. The majority (90%) of these PSs are unassociated.

Consider, first, the case where the 3FGL sources in the IG are unmasked. The left panel of Fig. 2 shows the best-fit source-count function for the NFW PS (dashed, orange), isotropic PS (dotted, green), and disk PS (solid, blue) populations. The disk-correlated PS template accounts for the high-flux 3FGL sources. Below  photons/cm/s, the NFW PS template accounts for nearly all the PS emission; its source-count function has a steep cutoff just below the source sensitivity threshold. It is worth noting that there is no externally imposed threshold for the PS population in this case, as the 3FGL sources are not masked.

Figure 4: Results obtained by applying the NPTF to simulated data. (Left column) The source-count functions for the PS templates in the fit when NFW PSs are included in the simulated data (top) or not (bottom). Note that when NFW PSs are not simulated, an NFW DM component is instead. (Right column) The associated posteriors for the fraction of flux absorbed by the different templates in the fit. The inset plots show the results of analyzing the simulated data without an NFW PS template in the fit. All plots are relative to the region within 10 of the GC with and 3FGL sources unmasked. For the flux-fraction plots, the fractions are computed relative to the total number of counts observed in the real data.

The most pressing question to address is whether the excess flux in the IG is better absorbed by the NFW PS or NFW DM template. The right panel of Fig. 2 shows the respective flux fractions, computed relative to the total photon count in the inner region with and the 3FGL sources unmasked. The disk and isotropic PSs contribute and of the flux, respectively. In contrast, the best-fit flux fraction for the NFW PS component is , while the best-fit DM flux fraction is consistent with zero. The normalization of the diffuse model remains consistent (to within 1%) with expectations from the fit at high latitudes, suggesting that the NFW PS template is absorbing the excess, and only the excess, and corresponds to a source population distinct from the more disk-like population of resolved sources. When the NFW PS template is omitted (inset), the fraction of flux absorbed by the disk PS population is essentially unchanged at , and the DM template absorbs of the flux. The DM flux obtained in absence of an NFW PS template is consistent with other estimates in the literature [12, 14]. The model including the NFW PS contribution is preferred over that without by a Bayes factor .5

When the 3FGL sources are masked, the NPTF procedure yields a best-fit source-count function given by the orange band in the left panel of Fig. 3. Below the break, the source-count function agrees well with that found by the unmasked fit. In this case, the contributions from the isotropic and disk-correlated PS templates are negligible. The flux fraction attributed to the NFW PS component is , while the NFW DM template absorbs no significant flux.

In the masked analysis, the Bayes factor for a model that contains an NFW PS component, relative to one that does not, is 10, substantially reduced relative to the result for the unmasked case. Masking the 3FGL sources removes most of the ROI within of the GC, reducing photon statistics markedly, especially for any signal peaked at the GC. Furthermore, in the masked ROI, non-NFW PS templates can absorb a substantial fraction of the excess. For example, if only disk and isotropic PS templates are added, the flux fraction attributed to the disk template is %, while that attributed to NFW DM is % (the flux attributed to isotropic PSs is negligible). When no PS templates are included in the fit, the NFW DM template absorbs % of the total flux. As we will discuss later, this behavior agrees with expectations from simulated data. In this statistics-limited regime, the fit does not distinguish different models for the PS distribution at high significance,6 but there is still a strong preference for unresolved PSs. The Bayes factor for a model with disk and isotropic PSs, relative to one with no PSs, is 10, while the Bayes factor for a model with NFW, disk and isotropic PSs, relative to one with no PSs, is 10. The Bayes factors, best-fit source-count function parameters, and DM flux fractions for the 3FGL masked and unmasked analyses are summarized in Tab. 1.

To validate the analysis procedure, we generate simulated data using the best-fit parameters from the unmasked IG analysis; we include isotropic, Fermi bubbles, and Fermi p6v11 diffuse emission, as well as isotropic, disk and NFW-distributed PSs. The simulated data is then passed through the 3FGL-unmasked IG analysis pipeline described above. Details for how we perform the simulations may be found in the Supplementary Material.

The top row of Fig. 4 shows the source-count functions that are recovered from the NPTF (left), as well as the posterior distributions for the flux fractions of the separate components of the fit (right). The fitting procedure attributes the correct fraction of flux to NFW-distributed PSs, within uncertainties, and finds no evidence for NFW DM. When no NFW PS template is included in the fit (inset, top right), the NFW DM template absorbs the excess. Both the source-count functions and the flux fractions are consistent with the results obtained using real data. Additionally, we recover a Bayes factor of 10 in preference for NFW PSs when using the simulated data, which is similar to what we found for the actual analysis.

For comparison, the bottom row of Fig. 4 shows the result of running the NPTF on a simulated data set that does not include NFW-distributed PSs but does include NFW DM. The model parameters used to generate the simulated data are taken from the best-fit values of the analysis without NFW PSs on the real data. In this case, the fitting procedure finds no evidence for NFW PSs, as it should, and the Bayes factor in preference for NFW PSs is much less than 1. The source-count functions recovered for disk-correlated and isotropic PSs are consistent with those used to generate the simulated data.

The source-count function that we recover for NFW PSs in the IG differs at low flux from those previously considered in the literature, which were motivated by population models and/or data for disk MSPs [41, 23, 24, 19]. In particular, our source-count function seems to prefer an increasing below the break, implying most sources lie close to the cutoff luminosity, while previously-considered source-count functions tend to be flatter or falling in . If confirmed, this may suggest novel features of the source population; however, our results are also consistent with a flat or falling within uncertainties.

The results of the NPTF analyses presented here predict a new population of PSs directly below the PS-detection threshold in the IG. We estimate from the 3FGL unmasked (masked) analysis that half of the excess within of the GC with may be explained by a population of () unresolved PSs, with flux above () photons. The entire excess within this region could be explained by () PSs, although this estimate relies on extrapolating the source-count function to very low flux, where systematic uncertainties are large. Detecting members of this PS population, which appears to lie just below the current Fermi PS-detection threshold, would be convincing evidence in favor of the PS explanation of the GeV excess.

Acknowledgements.— We thank S. Ando, K. Blum, D. Caprioli, I. Cholis, C. Dvorkin, D. Finkbeiner, D. Hooper, M. Kaplinghat, T. Linden, S. Murgia, N. Rodd, J. Thaler, and N. Weiner for useful discussions. We thank the Fermi Collaboration for the use of Fermi public data and the Fermi Science Tools. B.R.S was supported by a Pappalardo Fellowship in Physics at MIT. This work is supported by the U.S. Department of Energy under grant Contract Numbers DE-SC00012567, DE-SC0013999 and DE-SC0007968. B.R.S. and T.R.S. thank the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1066293, for support during the completion of this work.

Note added in proof.—Recently, we became aware of Ref. [42], which studied the distribution of sub-threshold PSs in the inner Galaxy using a wavelet technique.

Evidence for Unresolved Gamma-Ray Point Sources in the Inner Galaxy

Supplementary Material

Samuel K. Lee, Mariangela Lisanti, Benjamin R. Safdi, Tracy R. Slatyer, and Wei Xue

The supplementary material is organized as follows. The analysis methods are described in Sec. A. Section B estimates the Fermi PS-detection threshold. Extended results for the analyses described in the main text are then given in Sec. C. Section D investigates a variety of systematic issues that may affect the Inner Galaxy analysis. In Sec. E, we make a detailed comparison to previous work on the luminosity function for PSs in the Inner Galaxy. Section F describes an additional test of the unresolved PS models utilizing the survival function. Lastly, Sec. G describes how the NPTF is validated using simulated data in the Inner Galaxy.

Appendix A Methods

This section describes in detail the analysis framework for the NPTF and the dataset used in this work.

a.1 Analysis framework

In the standard template-fitting procedure, the photon-count probability distribution of observing photons in a given energy bin and pixel is assumed to be Poissonian,


where is the expected number of photons in the energy bin, summed over all the templates that contribute in that pixel. In particular,


where is the best-fit normalization of the template in the given energy bin, and specifies the spatial dependence of the template. That is, gives the expected number of events in the energy bin and pixel (accounting for an energy-dependent exposure, which may vary across the sky) for that template. Past studies of the GeV excess have included templates that account for the spatial dependence of the diffuse and isotropic backgrounds, the Fermi bubbles, and an NFW-distributed DM component. We have repeated the standard template-fitting procedure and verified that our analysis pipeline reproduces results from those studies (see Sec. C.2 for details).

Unlike the standard template-fitting procedure, the NPTF does not assume that the photon counts are Poissonian. We compute the non-Poissonian photon-count probabilities using the method of generating functions. In particular, the total generating function for the photon-count probability distribution in each pixel is


where is an auxiliary variable, and and are the pixel-dependent generating functions for the diffuse and non-Poissonian PS contributions, respectively. The probability of observing photons in pixel is then


If , then (S4) reduces to the simple Poissonian form found in (S1). Explicit analytic formulae for these generating functions are derived in [31, 32].

The benefit of the NPTF procedure is that it allows one to include additional spatial templates in the fitting procedure that are not properly characterized by Poissonian statistics—in particular, it allows for PS templates. In this work, we include up to seven different templates, with up to 16 free parameters. These parameters are listed in Tab. S1, along with the prior ranges over which they are scanned. The Bayesian analysis results in posterior probability distributions for each of these parameters. We now explain each of the templates and the corresponding parameters in detail.

  • Isotropic Background
    The isotropic-background template, which is uniform over the sky before correcting for exposure, is specified by the normalization parameter in the NPTF. We define such that the number of counts arising from the isotropic background in pixel is given by ; that is, when the NPTF and the standard template fit agree exactly.7 The prior for is taken to be linear flat; the posterior is well-converged and is typically peaked around unity.

    The best-fit pixel-averaged intensity for the isotropic background in a given ROI is


    where is the pixel solid angle, is the exposure in pixel , and the angle brackets denote an average over the pixels in the ROI. Note that for the isotropic-background template.

    Parameter Prior Ranges
    High Latitude Inner Galaxy
    Table S1: Parameters and associated prior ranges for the high-latitude and IG analyses. Note that there can be up to three copies of the PS parameters when isotropic, disk-correlated, and NFW-distributed PS templates are all included in the fit. The break is scanned up to , the maximum number of photons observed in a given pixel within the ROI; in general, photons.
  • Diffuse Background
    The spatial dependence of the diffuse-background template is based on the Fermi p6v11 model for the majority of our analyses; however, we do perform a number of cross-checks with other diffuse models. As with the isotropic-background template, the diffuse-background template is specified by a single normalization parameter, , also defined such that when the NPTF agrees exactly with the standard template fit.8 We assume a linear-flat prior, and the resulting best-fit value is typically close to unity. The best-fit pixel-averaged intensity is computed as in (S5).

  • Fermi Bubbles
    The Fermi-bubbles template assumes uniform emission within the bubbles [34] and is also specified by a single normalization parameter, , defined analogously as above. We similarly assume the prior for is linear flat; the best-fit value is again typically close to unity. The best-fit pixel-averaged intensity is computed as in (S5).

  • NFW Dark Matter
    The spatial dependence of the DM template is determined by , which gives the integrated photon-count intensity at the center of pixel from DM annihilation in the Galactic halo. The intensity is computed by the line-of-sight integral


    where is the DM density profile and gives the distance from the GC. This work assumes a generalized Navarro-Frenk-White (NFW) density profile [35, 36]


    where  kpc is the scale radius. We take for the majority of the analyses.

    With the spatial dependence of fixed by the NFW intensity profile and exposure map, the DM template can be specified by a single normalization parameter, , defined analogously as above. We assume a prior that is log flat and covers a broad range (see Tab. S1). The normalization of is such that it is equal to unity for 35 GeV DM annihilating into with cross section  cm/s, and a local DM density at the solar circle of 0.3 GeV/cm.

  • NFW Point Sources
    Previous standard template fits have included the above Poissonian templates. However, in the NPTF, PS templates have non-Poissonian statistics, which can be derived from the source-count function. In the majority of this work, the PS source-count function in a given pixel is modeled as a broken power law


    where is the break, are the slopes above and below the break, and is a pixel-dependent normalization. We require and so that the total number of photons contributed from PSs is finite. The priors for the parameters , , and are all linear flat.

    For the NFW PS template, the pixel-dependent normalization of the source-count function, , is assumed to be proportional to , as would be needed for the PSs to mimic a DM signal.9 We can thus specify the NFW PS template by the three source-count function parameters and an overall normalization parameter, , defined such that


    As with , the prior for is log flat and covers a broad range.

    We also quote the pixel-averaged intensity , which is given by


    with the expected number of counts in pixel arising from the PS population:

  • Isotropic Point Sources
    The source-count function for the isotropic PS template is also modeled as a broken power law, as in (S8), except with . As above, the isotropic PS template can be specified by either its overall normalization parameter or the pixel-averaged intensity .

  • Disk-correlated Point Sources
    The thin-disk template is the projection along the line-of-sight of the source spatial distribution given by:


    where , are cylindrical polar coordinates. The source-count function associated with this template is modeled as a broken power law, as in (S8), except with .

Given these templates and their associated generating functions, the overall photon-count probability distribution can be written as a function of the 16 parameters


Then, for a data set consisting of the set of photon counts in each pixel , the likelihood function for observing a particular photon-count distribution over all pixels in the ROI is


With the priors specified above, this likelihood function can be used in the standard framework of Bayesian inference to compute both the posteriors and the evidence for models that include various subsets of the parameters . We use the MultiNest package for the Bayesian calculations [37, 38]. All MultiNest runs use 500 live points with importance nested sampling disabled, constant efficiency mode set to false, and sampling efficiency set for model-evidence evaluation; the typical number of posterior samples generated for each run is 10.

a.2 Data Selection Criteria

The NPTF analysis was performed using the Extended Pass 7 Reprocessed Fermi data from August 4, 2008 to December 5, 2013 made available by [29]. Ultraclean front-converting events with zenith angle less than 100 and “DATA_QUAL==1 && LAT.CONFIG==1 && ABS(ROCK.ANGLE) < 52” are selected, and a Q2 cut on the CTBCORE parameter is used to remove events with poor directional reconstruction. We have used the original CTBCORE-cut dataset throughout this work [29], but we have tested the effect of including Fermi data up to March 8, 2015 (mission week 353) and the results are consistent with those presented here. Additionally, we have rerun the analysis with Pass 8 Fermi data up to June 3, 2015, using the new Ultracleanveto class and the top quartile of events ranked by PSF (denoted PSF3). We adopt the recommended data quality cuts for this analysis,10 which are slightly different to those in our main analysis (e.g., the zenith angle cut has been reduced from to and the rocking angle cut has been removed). Again, the results are consistent with those presented in this work.

The main body of the Letter focused primarily on two regions of interest: a high-latitude analysis with and an IG analysis that included all pixels within of the GC, with . These regions are shown in Fig. S1. When masking identified PSs from the Fermi 3FGL catalog [33], all pixels within of the source are excluded. This mask is sufficiently large to completely contain the flux from the majority of the PSs; the results do not qualitatively change as the mask size is varied, for example, to .

Figure S1: The counts map for the high-latitude () analysis (clipped at 15 counts) is shown in the left panel. The IG analysis focuses on the region within of the GC, with . The associated counts map (clipped at 50 counts) is shown in the right panel. All pixels within 1 of known Fermi 3FGL sources are masked.

If the detector’s point spread function (PSF) is comparable to or larger than the pixel size, photons from a PS in a given pixel can leak into neighboring pixels. As a result, the PSF must be properly accounted for in the calculation of the photon-count probability distributions for the PS templates [31, 32]. To model the PSF for the CTBCORE-cut data, we use the instrument response functions made available by [29] to approximate the PSF as a two-dimensional Gaussian with energy-dependent width, (see Sec. D.4 for details).11 The NPTF analysis includes data from a single energy bin from 1.893–11.943 GeV. Small modifications to the choice of energy bin do not affect the conclusions. For these energies, the PSF parameter varies from at low energies to at high energies. In calculating the photon-count probability distributions for the PS components, we use the energy-averaged value for assuming the spectrum of the excess (high-latitude 3FGL PSs) for all PS templates in the IG (high-latitude) analysis. We further explore systematic uncertainties arising from the PSF in Sec. D.4.

Appendix B The Fermi-LAT Point-Source Detection Threshold

The analysis presented in this Letter does not rely on any prescription for the sensitivity of Fermi to PSs, except implicitly via masking out the known sources. However, it is still worth considering whether the PS population we infer should already have been resolved as individual sources, given the nominal sensitivity of Fermi.

To some degree, the detection threshold can be read off directly from figures such as Fig. 2; the threshold may be responsible for the turnover in the histogram of known sources between 2–4 photons/cm/s. However, our analysis is restricted to a single high-energy bin; the “faintest” resolved sources in our sample include sources that are not intrinsically faint but have very soft spectra, and so have very few photons at high energy. The analysis by which the 3FGL catalog was created used a wider range of energies, and so it is entirely possible that a source might be detected by its bright emission at low energies, yet be too faint to be independently detected in our 1.893–11.943 GeV energy range, thus giving the false appearance that very faint sources can be detected. The spectrum of the excess is rather hard relative to the known 3FGL sources in our ROI, so PSs comprising the excess will tend to become less detectable relative to the resolved 3FGL sources in lower energy bands.

Estimates for the PS sensitivity in the literature generally assume a specific source spectrum and are based on fewer years of Fermi data than our current analysis (although they may have similar statistics, because we are using front-converting Ultraclean events and have applied further cuts on the CTBCORE parameter). The Fermi Collaboration has presented longitude-averaged sensitivity estimates based on the integral flux from 0.1–100 GeV with three years of data, for sources with a pulsar-like spectrum [43]; for Galactic latitudes with , the mean sensitivity varies from erg/cm/s, with the percentile sensitivity ranging from erg/cm/s. The average over longitude probably leads to a somewhat lower threshold than is accurate in the neighborhood of the GC, but on the other hand, the 3FGL catalog is based on four years of data rather than three.

Let us assume the spectrum of the excess, and all sources comprising it, is approximately given by the best-fit broken power-law spectrum of [14], with indices 1.42 and 2.63 below and above the break, respectively, and a break energy of GeV. Then, a luminosity of erg/cm/s in the 0.1–100 GeV band corresponds to photons/cm/s in our energy band. Thus, a threshold of erg/cm/s in integrated energy flux translates to a photon-flux threshold photons/cm/s in the relevant energy range, for the percentile sensitivity 10 from the Galactic plane. This is directly above the break in our inferred source-count functions.

Appendix C Extended results

This section includes extended discussion of the results presented in the main text.

c.1 High latitudes

Fig. S2 shows the posterior probabilities for the photon intensity associated with each template in the high-latitude analysis (), with 3FGL sources unmasked and masked (left and right columns, respectively). The top row shows the results when only the isotropic, diffuse, and bubbles templates are included in the fit. The best-fit intensity values for the diffuse and bubbles components are essentially constant, regardless of whether or not the 3FGL sources are masked. However, the isotropic template is sensitive to the presence of the 3FGL sources, with its median best-fit intensity increasing from to  photons/cm/s/sr when going from the right to left panels.

When the isotropic PS template is included in the NPTF at high-latitudes (bottom row in Fig. S2), there is no effect on the flux absorbed by the diffuse and bubbles template. In this case, the isotropic template is insensitive (within uncertainties) to the masking of the 3FGL sources, with a photon intensity that remains unchanged between the left and right panels. Instead, the isotropic PS template accounts for the resolved PSs when the 3FGL sources are unmasked, and picks up flux (presumably from unresolved PS emission) when the 3FGL sources are masked.

Figure S2: The posterior probabilities for the photon intensity associated with a given template in the high-latitude analysis with 3FGL sources unmasked (left column) and masked (right column). The top row shows the result without a PS template (i.e., the standard template fit), while the bottom row shows the result of the NPTF with an isotropic PS template included. Whether or not the 3FGL sources are masked, the isotropic template in the standard template fit clearly absorbs intensity that is assigned to the isotropic PS template in the NPTF. Dashed vertical lines indicate the , , and percentiles.

The intensity of the IGRB that is predicted by the NPTF can be compared with published values from the Fermi collaboration [40]. Applying a standard template analysis to 50 months of Fermi LAT data, [40] determined the IGRB spectrum from 100 MeV to 820 GeV. The brightest PSs from the 2FGL catalog (the predecessor of the 3FGL catalog) were fitted individually, while a standard template was used to model the other identified sources. Three foreground models—labeled A, B, C—were studied to better understand the effects of variations in the diffuse emission. For each foreground model, the spectrum of the IGRB intensity was fit to a broken power law with an exponential cutoff. According to our estimates, the results in [40] predict IGRB intensities of , ,  photons/cm/s/sr in the energy range from 1.893–11.943 GeV for the three models, with systematic uncertainties on the order of in each case. The intensity of the IGRB computed using the NPTF (e.g.,  photons/cm/s/sr for the 3FGL-masked fit) is generally in agreement with the results of [40]. However, a direct comparison is difficult to make because, for example, the analyses treat the known PSs differently, use different PS catalogs, use different data sets, and use different foreground models. It would be interesting to repeat the analysis in [40] including a non-Poissonian PS template to more accurately extract the contribution of unresolved PSs to the total EGB.

The central analyses in this work assume the source-count function parameterization in (S8), however we also consider the effect of forcing the source-count function to have a sharp cutoff at high flux, while still allowing a break at lower flux. For example, we can impose a high-flux cutoff consistent with the flux of the brightest 3FGL PS in the ROI: photons. In this case, the NPTF gives the following best-fit values: , , and photons. Additionally, photons/cm/s/sr, which—subtracting the intensity of the identified 3FGL PSs—predicts that the unresolved PSs have an intensity photons/cm/s/sr. This is in good agreement with the estimate obtained from the 3FGL-masked NPTF described in the main body of the Letter (e.g., photons/cm/s/sr). Indeed, the agreement is better than that obtained from estimating the intensity of the unresolved PS population using the source-count function without the high-flux cutoff (e.g., photons/cm/s/sr).

c.2 Inner Galaxy

Figs. LABEL:Fig:_IG_triangle_nomask and LABEL:Fig:_IG_triangle_mask5 show the posterior probabilities for the IG analysis (within of the GC, with ), with 3FGL sources unmasked and masked, respectively. In both cases, the parameters are all well-constrained within the prior ranges, with one exception.12 Namely, the slope above the break of the PS source-count functions may have large error bars. This tends to happen in the masked analyses, where the source-count functions fall off steeply near the detection threshold.

As shown in Figs. LABEL:Fig:_IG_triangle_nomask and LABEL:Fig:_IG_triangle_mask5, the NPTF finds that the average intensity of the diffuse emission in this region is () photons/cm/s/sr with 3FGL sources unmasked (masked); the predicted intensities are similar in both scenarios, as desired. Additionally, these intensities change by less than from the respective values obtained by the standard template procedure (i.e., when only an NFW DM template is included). The fact that the intensity of the diffuse emission is essentially constant between the NPTF and standard template analysis highlights that the addition of the PS templates does not affect the flux absorbed by the diffuse background template in the IG region.

One cross-check of the results is to compare the predicted fractions of flux from DM in the template fits that do not include PSs—where the excess appears to be absorbed by the NFW DM template—to results found in previous template studies. Ref. [12] performed two analyses, with different ROIs, that are relevant for this comparison.13 The ROI for the first analysis was a region around the GC, while the second analysis used the full sky. In both cases, the plane was masked () in addition to the 300 brightest and most variable PSs in the 2FGL PS catalog, using estimated 95% containment masks. A spectral and spatial model was also included for the remaining 2FGL sources, based on their positions and spectra in the catalog. Ref. [12] used the same Poissonian templates as we do, so we can compare their results to our fits that do not include NFW PSs.

Using the results of [12], we can compute the predicted fraction of flux from DM relative to the total number of observed counts within of the GC, with , no PSs masked, and within the energy range considered in this Letter. For the restricted (full-sky) ROI, we find that [12] predicts the fraction of flux from DM to be (%). The difference between these values gives a sense for the systematic uncertainty that comes from changing ROIs. Given the fact that our analyses use different ROIs from [12], we consider our results to be consistent with theirs within systematic uncertainties. But importantly, regardless of the ROI, we find that when we include an NFW PS template in addition to an NFW DM template, the excess flux is preferentially absorbed by the PS template.

For the results presented in this work, we have used a double-exponential thin-disk template with scale height and radius of 0.3 and 5 kpc, respectively. However, we find that our results are not sensitive to variations in the disk template. For example, we have repeated the masked and unmasked analyses using a thick disk with scale height and radius of 1 and 5 kpc, respectively. This thick-disk distribution has been shown to be a good model for the distribution of MSPs identified by Fermi [19, 23, 44]. We have also considered the case of a thin disk with scale height and radius of 0.3 and 10 kpc, respectively. In both these cases, the results of our analyses remain essentially the same.

One substantial difference between our work and previous studies of the excess is that we include no energy information, using only a single large energy bin. In future work, it would be very useful to incorporate energy dependence into the NPTF likelihood function in order to extract a spectrum for the PS population. It would also be useful to understand whether DM substructure could give rise to all or part of the excess. In this Letter, we have modeled DM emission as smooth emission, but it may be the case that the DM emission is more PS-like because of, for example, DM subhalos.

Appendix D Systematic Uncertainties in the NPTF

There are various systematic uncertainties that may influence the IG analysis, which we examine in some detail in this section.14

d.1 Broadening the ROI

As a first cross check of our results, we perform the NPTF on the full sky with . This analysis can be used to test the North-South symmetry of the excess by masking each hemisphere in turn. The results for the full-sky and hemisphere analyses are similar, so only the latter are presented here. Additionally, we only show results for the 3FGL-unmasked analysis, since, when when masking 3FGL sources, the fraction of the sky masked near the GC is different in the North versus the South.

Fig. S3 shows the best-fit source-count function and flux fractions for the Northern and Southern hemispheres in the top and bottom rows, respectively. The source-count functions for the two regions are consistent within statistical uncertainties. In both cases, the orange band cuts off steeply around photons/cm/s. The inferred flux fraction of NFW PSs in the Northern (Southern) hemisphere analysis is % (%) in the region within 10 of the GC, with , while the Bayes factor in favor of the model with NFW, disk, and isotropic PSs over the model without NFW PSs is () in the Northern (Southern) hemisphere analysis. We conclude that there are no significant asymmetries in the inferred NFW PS population between the Northern and Southern hemispheres.

Figure S3: Results obtained by applying the NPTF to a full-sky map with either the Northern (top row) or Southern (bottom row) hemispheres masked. The templates included are: isotropic, diffuse, bubbles, NFW DM, NFW PSs, disk PSs, and isotropic PSs. The best-fit source-count functions (with 68% confidence intervals shaded) are shown in the left column, and the posterior probabilities for the flux fractions are shown in the right column. The source-count functions are plotted with respect to the region within 10 of the GC with () for the Northern (Southern) analysis. The number of 3FGL sources in these regions is indicated. The flux-fraction plots are calculated for the region within 10 of the GC with in both cases. The inset shows the posterior probabilities for the flux fractions when the NFW PS template is not included.

d.2 Varying the Diffuse Model

Figure S4: Same as Fig. 2, except using the Fermi p7v6 diffuse background model.

A potentially serious source of systematic uncertainty is due to limitations in modeling the diffuse gamma-ray background arising from the propagation of high-energy cosmic rays in the Galaxy. The primary contributions come from bremsstrahlung of high-energy cosmic rays passing through the interstellar gas, inverse Compton scattering of Cosmic Microwave Background, interstellar, and infrared radiation off high-energy electrons, and the decay of boosted pions produced in inelastic proton collisions with the interstellar gas. Modeling the cosmic-ray emission depends on many factors, including the location and spectrum of source injection, the gas distribution, magnetic fields, and the interstellar radiation field, as well as the diffusion parameters. Repeating the NPTF with different diffuse-background models can help to quantify the effects of mis-modeling the diffuse background.

The primary results presented in this work used the Fermi p6v11 diffuse model. This model is not the most recent to be released by Fermi. The more recent Fermi p7v6 diffuse model includes contributions from large-scale diffuse substructures such as Loop 1 and the bubbles. As found in [12], repeating the template analysis with the p7v6 model does not qualitatively affect the results for the GeV excess, except for the fact that the overall flux absorbed by the DM template is reduced. This may be due to the largely data-driven p7v6 diffuse model having absorbed part of the GeV excess.

Figure S5: Same as Fig. 3, except using the Fermi p7v6 diffuse background model.

Figure S6: Same as Fig. 3, except showing the median source-count functions for the additional thirteen diffuse models described in the text in dashed black; the median function for the p6v11 model is shown in solid blue, for comparison. The orange band is the overlap of the 68% confidence bands for all thirteen models. For computational reasons, the only PS template included in these tests is that for NFW PSs.

The left panel of Fig. S4 shows the best-fit source-count functions obtained for the IG study using the p7v6 diffuse model, with all other templates the same. Comparing to Fig. 2, we see that the overall features are recovered. The right panel of Fig. S4 shows the corresponding flux fraction plot. In this case, the flux fraction absorbed by the NFW (disk) PS template is ; when the NFW PS template is removed from the fit (inset, right panel), the majority of its corresponding flux is absorbed by the NFW DM template instead. The main difference with the p6v11 results is that the flux associated to NFW PSs is lower with p7v6, which is consistent with previous findings of a smaller flux fraction for the GeV excess in p7v6 studies. As a result, the Bayes factor in preference for the NFW PS template is reduced to in this case.15

Figure S5 shows the corresponding results for p7v6 when the 3FGL sources are masked. Now, the NPTF finds no excess flux in the ROI; the flux fractions for the PS templates, as well as the NFW DM template, are all consistent with 0%. The fact that no flux is picked up by the NFW DM template is different from what was previously observed in [12]. We have verified that this is due to the larger PS mask implemented here, which removes a considerable region close to the Galactic Center; repeating the analysis with the PS masking of [12] recovers their result.

In addition to studying the p7v6 diffuse model, we also consider thirteen other diffuse models from [45]: L6201000005, L6201502, L6201505, L6201000002, L10201505, O10201505, Y10201505, S10201505, S8201505, S8201505, S4201505, O10301505, L10301505. These models are chosen to span variations in the source distributions, the diffusion scale height and radius, the gas spin temperature, and cuts on the magnitude of E(B-V). For these models, we include separate templates for the , bremsstrahlung, and ICS components. The prior ranges for each of these three templates is dealt with in the same way as the diffuse-model prior ranges for the p6v11 and p7v6 diffuse models; in all cases, the three component templates are well-converged within the prior ranges. For brevity, we only show results for the 3FGL-masked analyses. Additionally, for computational reasons, we only include a non-Poissonian template for NFW PSs. When only this PS template is included, the Bayes factor in preference for the model with NFW PSs relative to that without is () for the p6v11 (p7v6) diffuse model.

Fig. S6 summarizes the results, showing the spread in the median source-count functions (dashed black lines) obtained by running the NPTF in the IG with each of these diffuse models. The orange band indicates the combination of the thirteen 68% confidence intervals. While there is some spread in the predicted source-count function, the NPTF consistently finds a non-zero flux for the NFW PS contribution in all cases, ranging from 3% to 9%; also, the DM flux is always consistent with zero when the NFW PS template is included in the fit. Specifically, the model including an NFW PS template is preferred over the model without such a template by Bayes factors in the range for all diffuse emission scenarios considered. In this sense, these thirteen diffuse models appear more similar to p6v11 than to p7v6.

For the model L6201000005, which provided the formal best fit in a previous analysis of the IG [14], we also test the effect of adding an independent diffuse template with free normalization, corresponding to the predicted gas-correlated emission (pion decay and bremsstrahlung) within the innermost Galactocentric ring. The addition of this template does not significantly alter the source-count function or flux fraction attributed to the NFW PSs.

d.3 Scan Along the Galactic Plane

Figure S7: Same as Fig. 2, except repeating the NPTF in the region within of , with .

To address the concern that the PS templates may be erroneously absorbing contributions from the diffuse background, we repeat the NPTF in regions centered at different points along the Galactic plane. Previous studies have reported additional bright excesses along the plane, with the most significant at  [12, 14]. The residual emission in this region is roughly similar to that at the GC; however, its energy spectrum is softer, suggesting a different origin. Repeating the NPTF on another region of sky—with an excess that is unlikely to arise from PSs—allows us to test whether the fitting procedure can adequately distinguish between extra diffuse emission on top of a mis-modeled diffuse background and extra PS emission.

We apply the NPTF to regions of the sky within of the central points , where , requiring throughout. The same templates are included as in the IG analysis, except the NFW templates are centered at the middle of each ROI. The most significant excess is found for the region centered at , consistent with previous findings. Fig. S7 shows the best-fit source-count functions for the different PS templates in the 10 region centered at this point. The disk template (solid, blue) successfully recovers the known PSs in this region of sky (numbering 13). Most notably, the error band on the source-count function for the NFW PS population (dashed, orange) is consistent with no unresolved PSs in this ROI.

The right panel of Fig. S7 shows that the flux fractions of the NFW DM and disk PS components are comparable in this region, while the isotropic and NFW PS templates pick up negligible flux. Excluding the NFW PS template (inset) does not significantly affect the DM and disk PS flux fractions. The Bayes factor in favor of the model with NFW PSs relative to that without is 0.1. We note that repeating the analysis with the exact masking criteria and NFW inner-slope value used by [12], we reproduce the (slightly smaller) flux fraction of the excess found by those authors.

d.4 Point Spread Function

Figure S8: Same as the left panel of Fig. 2, except setting the PSF parameter (left) and (right). When masking identified 3FGL sources, all pixels within are excluded.

An accurate PSF is a crucial component of the NPTF procedure, as an incorrect parametrization can lead to over- or under-estimation of the PS component. The instrument response functions for the Fermi CTBCORE-cut data were made available by [29]. For a given energy range, the authors of [29] provide the 68% and 95% containment radii, and .

We model the PSF as a two-dimensional Gaussian


setting such that the 68% containment radius of (S15) is equal to . The value of varies from 0.198 at the lowest end to 0.0492 at the highest end of the energy bin (1.893-11.943 GeV). As described earlier, the isotropic, Fermi bubbles, and NFW DM templates are smeared with this PSF. (The diffuse-background template is smeared with the exact energy-dependent PSF using the Fermi Science Tools.) The PSF must also be properly accounted for in the generating-function formalism that is used to obtain the non-Poissonian likelihood function (see [32] for further details).

It is known that the real PSF has power-law tails that are not captured by (S15[29]. One might rightfully be concerned that ignoring these power-law tails can lead to mischaracterization of the PS population. To illustrate the effect of mis-modeling the PSF, Fig. S8 shows the best-fit source-count function for the IG analysis assuming extreme values for the Gaussian-PSF parameter : in the left (right) panel. When the assumed PSF is too wide, photons in neighboring pixels arising from diffuse emission or unresolved sources may instead be erroneously interpreted as photons from a single PS. Thus, using a PSF that is wider than the true PSF will underestimate the number of unresolved PSs and overestimate the number of observed sources, relative to the true source-count distribution. This effectively shifts the source-count function to large flux values, as observed in the left panel of Fig. S8.

In comparison, using a PSF with a width that is too small underestimates the number of bright PSs. For example, in the extreme limit, assuming a delta-function PSF would require all photons from a given PS to be contained in a single pixel, while the photons from bright PSs are, in actuality, smeared over several pixels by the true PSF. In this case, as may be seen in the right panel of Fig. S8, the fitting procedure predicts fewer high-flux PSs, while at the same time predicting a larger population of unresolved sources.

The fact that the energy-averaged PSF used in Fig. 2 results in a source-count function that matches the high-flux observed PSs reasonably well gives us confidence in our treatment of the PSF. Additionally, we have illustrated here that extreme variations to the PSF-parameter do not significantly affect the results. In Sec. G, we further verify our treatment of the PSF using simulated data.

d.5 Radial Distribution Profile

The analyses presented thus far have fixed the inner slope of the generalized NFW density profile to be for the DM and PS components. Due to computational limitations, we do not scan over this parameter. Previous studies of the GeV excess using standard template-fitting methods, such as [12, 14], have scanned over this parameter and found best-fit values in the range between and . Variations in have negligible effects on our final result. Fig. S9 shows the best-fit source-count function in the inner 10 with , assuming and (left and right panels, respectively). As is evident, the source-count function for NFW PSs is not sensitive to the value of the inner slope in this range. The posterior probabilities for the flux fractions (not shown) are also not significantly changed from Fig. 2.

Figure S9: Same as the left panel of Fig. 2, except using (left) and (right) for the generalized NFW distribution.

d.6 Diffuse-Correlated Point-Source Template

Figure S10: The Fermi p6v11 diffuse background in the IG region (smoothed using Fermi Science Tools), with masked. Counts are clipped at 30.

The conclusions described in this Letter provide evidence that the GeV excess can be accounted for by a population of unresolved PSs. One possible systematic issue, however, is that the preference for unresolved PSs is driven by localized structure in the diffuse gamma-ray background that is not captured by our background model. As a test of this hypothesis, we repeat the NPTF in the IG adding a PS template that traces the diffuse model. Fig. S10 shows the diffuse background in the ROI. Because the flux from the diffuse emission is larger closer to the plane, a diffuse-correlated (diff-corr) PS template is not only sensitive to localized structure in the diffuse model, but also to the presence of a disk-like population of PSs. For example, if the unresolved PSs are preferentially located close to the plane, then they may be absorbed by the diff-corr PS template. Breaking the degeneracy between these two interpretations requires more careful study, which we postpone to future work. However, the preliminary results are illuminating, so we share them here.

The left panel of Fig. S11 shows the best-fit source-count functions obtained when doing the NPTF in the IG region (without masking the 3FGL sources), including a diff-corr PS template in addition to disk-correlated and NFW PS templates, along with the standard Poissonian templates. Here, the diff-corr PS template only extends to from the GC. To simplify the analysis, we do not include the subdominant isotropic PS template. The best-fit flux normalization for the diffuse background is consistent with that obtained from the high-latitude analysis. In addition, the best-fit NFW DM normalization is consistent with zero, and the recovered source-count function parameters for disk and NFW PSs are consistent with those found when not including the diff-corr PS template, since the diff-corr PS template does not absorb a significant fraction of the flux.

Figure S11: Best-fit source-count functions for PSs within of the GC and , with the 3FGL sources unmasked, for models with NFW PS (dashed, orange), diffuse-correlated PS (dotted, pink), and thin-disk PS (solid, blue) components. For this analysis, the NPTF includes an additional template corresponding to diffuse-correlated PSs. This new template is taken to have support either in the inner (left) or over the full ROI ( from the GC with ) (right). For the latter case, parameters for the thin-disk PS component are fixed to the best-fit values found in the standard unmasked IG analysis.

As a point of comparison, we repeat the procedure letting the diff-corr PS template have support over the full IG region. Now, there is a potential degeneracy between the diff-corr PS template, the disk-correlated PS template, and the diffuse template. To break some of this degeneracy, we fix the disk-correlated PS parameters to their best-fit values from the scan including disk-correlated, isotropic, and NFW PSs. These results are shown in the right panel of Fig. S11. Again, the best-fit NFW DM normalization is consistent with zero, however the source-count function for the NFW PSs changes. In particular, the source-count function for NFW PSs is shifted to lower flux, potentially suggesting that some of the near-threshold sources could either be more disk-like in morphology or associated with mis-modeling the diffuse background. However, the preference for NFW PSs remains high, with the model including NFW PSs preferred over that without by a Bayes factor . Unlike the previous analysis that used a truncated diff-corr PS template, here the best-fit normalization of the diffuse-background template is lower than its best-fit value at high latitudes. When the NFW PS template is not included, the normalization can be shifted down by as much as 20%; the best-fit normalization of the diffuse template is still shifted down by 10% when the NFW PS template is included.

We caution the reader that these results are subject to considerable uncertainty due to the large number of parameters in the fitting procedure and the potential degeneracies between them. In addition, this analysis does not appear to be robust to changing the size of the ROI; as the ROI is increased, the best-fit normalization of the diffuse background approaches its value from high latitudes, and the flux absorbed by the diff-corr PS template decreases.

d.7 Binned Source-Count Functions

Figure S12: Best-fit source-count function for NFW PSs within of the GC and , with the 3FGL sources masked, where the number of sources in each flux bin is allowed to float freely. Blue (pink) regions indicate () confidence intervals. The orange band is the 68% confidence interval from the 3FGL-masked NPTF with a broken power-law source-count function. The median number of sources attributed to each bin is indicated, along with the 68% confidence range.

The imposition of a broken power law for the source-count function might be over-constraining in some cases. For example, it could lead to the apparent exclusion of a DM component because the extrapolation to low flux of the source-count function from high-flux sources is already too large. Furthermore, the broken-power law prescription may yield misleading results for the true uncertainty on the source-count function at low flux. To address these concerns, we present preliminary results from an alternate analysis where the number of sources in each flux bin is allowed to float independently; we use the seven logarithmically-spaced flux bins shown in Fig. S12. Within each bin, is constant as a function of . The source-count function model parameters are seven normalization parameters, one for each bin, which are taken to have log-flat priors over the range indicated for in Tab. S1.

Using the binned source-count function, we perform the NPTF on the 3FGL-masked IG ROI. For simplicity, we leave out isotropic and disk-correlated PS templates from this analysis. This is justified by the fact that leaving out these two templates from the broken power-law analyses does not qualitatively affect the results for the NFW PS and NFW DM components. We find that the DM flux fraction is consistent with zero, while the NFW PS template absorbs the excess flux.

As shown in Fig. S12, we recover a source-count function broadly consistent with the original analyses assuming broken power laws. The blue (pink) bands indicate the 68% (95%) confidence intervals for the source-count function in each bin, while the solid black (dashed red) line shows the median (mean) of the distribution. The orange band shows the 68% confidence interval from the masked NPTF using the broken power-law formalism.

At low flux, the mean and median of the binned result differ by multiple orders of magnitude, reflecting the fact that the posterior distributions for these parameters are skewed. This is related to the fact that the low-flux bins are not well-constrained by the fit, so the posterior distributions for these parameters are heavily influenced by the log-flat prior distributions. In the broken power-law fit, only the overall normalization parameter was taken to have a log-flat prior range, while in the binned fit all PS parameters have log-flat prior ranges. This point, combined with the fact that the broken power-law is more constrained than the binned source-count function, leads to larger uncertainties at low (and, to some extent, high) flux in the binned analysis than in the broken power-law analysis.

An additional challenge with this method is that the source counts in neighboring bins are generally highly correlated, leading to large errors in individual bins; however, the total flux is still well constrained. For example, we find that % of the flux (in the inner 10 region with ) is associated with the NFW PS template. These results are consistent with the broken power-law results, within uncertainties. Similarly, the binned fit prefers a large number of additional unresolved sources with fluxes only slightly below the PS-detection threshold, but the exact bin in which these sources appear is more uncertain. In future work, we plan to refine the binned analysis to make it more robust to changes in bin size and prior assumptions.

Appendix E Comparison with Luminosity Functions in the Literature

Figure S13: Same as Fig. 2, but including comparisons to example source-count functions motivated by MSPs. In particular, the red lines assume an NFW spatial distribution and a gamma-ray luminosity function consistent with that of observed MSPs in the Milky Way [41]. Three scenarios are considered: a maximum luminosity cutoff of  erg/s (dotted red) or  erg/s (dashed red), as well as no cutoff (solid red). These curves are normalized to the flux of the excess in the ROI. The purple line (dot-dashed) shows the expected source-count function assuming the same MSP-motivated luminosity function, but a thick-disk spatial distribution; it is normalized to match the number of bright photons/cm/s MSPs and MSP candidates at high latitudes ().

The best-fit source-count function recovered by the NPTF for the unresolved NFW-correlated sources is quite different to those previously considered in the literature—see e.g.[32, 41, 23, 24, 19]. With our source-count function, both the number of sources and the flux are dominated by sources within a factor 2 of the break. For sources 8.5 kpc from the Earth and with the spectrum of the excess, this break corresponds to a luminosity above 0.1 GeV of erg/s.

We construct models of the expected source-count function using a luminosity function derived from observed MSPs in the nearby field of the Milky Way [41], for both NFW-distributed sources and the thick-disk model described by (S12) with scale radius and height of 5 kpc and 1 kpc, respectively. For the former case, we also consider the possibility that the luminosity function possesses a cutoff at or erg/s, following [23]. The thick-disk model has been normalized to produce the correct number of high-latitude () bright MSPs and MSP candidates, with photons/cm/s assuming the average pulsar spectrum, as presented in [19]. Because this model has the wrong spatial morphology to explain the excess, the purpose of showing this curve is to illustrate the likely contribution from a disk population of MSPs within the ROI. The models with NFW-distributed sources have instead been normalized to match the flux attributed to the NFW PS template in our analysis.

Fig. S13 shows the expected source-count functions, averaged over our standard ROI, for these scenarios, for the 3FGL-unmasked fit (Fig. 2). In agreement with [23], we find that for NFW-distributed sources with the luminosity function of [41] and the correct normalization to explain the excess, too many sources are predicted above the Fermi detection threshold when there is no luminosity cutoff (solid red) or when erg/s (dashed red). The case with a cutoff at erg/s (dotted red) evades this constraint, as expected, but requires new sources to explain the excess. Using the luminosity functions of [41], we find good agreement with the number of sources required to fit the excess as stated in [23].

The purple dot-dashed line in Fig. S13 shows the predicted source-count function for the thick-disk distribution. It is remarkably similar to the best-fit source-count function for the thin-disk PS model extracted from the data. In particular, we estimate the slope of the purple dot-dashed line at low flux in Fig. S13 to be , while the NPTF predicts the slope of the thin-disk PS source-count function below the break to be . In agreement with [19], the unresolved sources associated with such a population should not produce enough photons to explain the excess (as may be seen by comparison to the red lines, which have the correct total flux to explain the excess).

The source-count function for the NFW PSs rises sharply above the red curves at fluxes photons/cm/s. As this source-count function is very shallow below the flux threshold, the total number of sources is dominated by this relatively high-flux region, as is the total flux. For this reason, only sources are needed to account for the excess in the ROI, in contrast to the quoted in [23].16 With that said, we emphasize that estimates of the total number of NFW-distributed PSs based on the NPTF are highly uncertain and subject to large systematic uncertainties since the low-flux PSs are hard to constrain with our method.

Appendix F Survival Function

As a further cross-check that the PS identification is working self-consistently, we attempt to identify pixels that are likely to contain unresolved sources in the 3FGL-masked sky. To do so, we perform a standard template fit in the ROI (excluding the PS templates) and determine the best-fit reference model by taking the median values of the posterior distributions for the appropriate fit parameters. Using this reference model, the expected mean number of counts, , in a pixel can be obtained. From the observed counts map one can then determine the Poissonian cumulative probability to observe the real count, , in that pixel given the expectation of the reference model. The survival function for pixel is defined as


where denotes the cumulative probability of observing counts for a Poisson function with mean . For example, if the best-fit reference model predicts a total of photons in a given pixel and the observed number of counts is , then . In general, the smaller the value of in a pixel, the more likely it is that it contains an unresolved PS.

Figure S14: The fraction of pixels with , with defined in (S16). (Left) The result for the high-latitude analysis, with the 3FGL sources masked. The red crosses show the survival-function distribution for the observed data set, given the reference background model (diffuse+isotropic+bubbles). The bars indica