Evidence for Unresolved GammaRay Point Sources in the Inner Galaxy
Abstract
We present a new method to characterize unresolved point sources (PSs), generalizing traditional template fits to account for nonPoissonian photon statistics. We apply this method to Fermi Large Area Telescope gammaray 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 gammaray 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 gammaray excess in this region. The excess is fully absorbed by such a population, in preference to darkmatter annihilation. The inferred source population is dominated by nearthreshold sources, which may be detectable in future searches.
Darkmatter (DM) annihilation in the Galactic halo can contribute to the flux of highenergy 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 cosmicray 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 nonPoissonian template fit (NPTF). Our approach is modelindependent, in that we remain agnostic about the nature of the PSs. To verify the method, we use it to characterize unresolved gammaray PSs at high Galactic latitudes. These findings represent one of the most precise measurements of the contribution of PSs to the extragalactic gammaray 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 gammaray 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 DMlike (smooth diffuse) origin. The Supplementary Material provides further details on the method, as well as additional crosschecks 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 signaltobackground ratio in the region of interest (ROI) high, maintains a sufficiently large number of photons over the full sky, and keeps the pointspread function relatively small and energyindependent.
The analysis utilizes the photoncount probability distribution in each pixel. In general, a given model for the gammaray 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 photoncount 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 templatefitting procedure must be generalized to include nonPoissonian photon counts. In the NPTF procedure, depends on a potentially pixeldependent PS sourcecount function . The sourcecount function determines the average number of PSs within pixel that contribute photon flux between and . In this work, the sourcecount function is assumed to follow a broken power law, , with pixeldependent 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. Semianalytic methods for calculating the with a broken power law sourcecount function were developed in [31, 32].
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 NavarroFrenkWhite (NFW) [35, 36]distributed DM, assuming no substructure; (5) isotropic PSs; (6) diskcorrelated PSs, and (7) NFWdistributed PSs.
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 wellconverged.
We begin by applying the NPTF to data at high Galactic latitudes (). Fermi’s third source catalog (3FGL) [33] identifies 1307 gammaray PSs in this region of the sky, with 55% associated with Active Galactic Nuclei and 24% associated with pulsars, supernova remnants and other known gammaray sources. The remaining 21% are unassociated. Figure 1 shows the sourcecount function in terms of the flux of the 3FGL PSs in our energy bin (black points). The observed sourcecount 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 highlatitude region, including templates for the diffuse background, Fermi bubbles, isotropic emission, and isotropic PSs. The bestfit sourcecount function values are given in Tab. 1.
The bestfit intensities obtained from the NPTF can be compared to those obtained from a standard template fit that neglects PSs. The diffusebackground and Fermibubbles intensities (averaged over the ROI) are consistent between both procedures. When the PS template is included, the isotropicbackground intensity is photons/cm/s/sr and the isotropic PS intensity is photons/cm/s/sr. With no PS template, the isotropicbackground intensity is over twice as high, photons/cm/s/sr. Thus, the PS intensity is absorbed by the isotropicbackground 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  —  —  —  
masked  
IG 
\makecellNFW PS
+ Disk PS 
unmasked  10  10  %  
IG 
\makecellNFW PS
+ Disk PS 
masked  10  10  %  
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 sourcecount 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 3FGLunmasked NPTF, within uncertainties. The intensity of the isotropic PSs is photons/cm/s/sr, which is slightly lower than the value inferred from the 3FGLunmasked NPTF.
corresponds to the intensity of the isotropic gammaray 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 3FGLunmasked NPTF, we infer that % of the EGB in this energy range is associated with both resolved and unresolved PS emission; from the 3FGLmasked 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].
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 NFWdistributed DM, in addition to isotropic, diskcorrelated, and NFWdistributed PSs.
While the prior ranges for the isotropic, isotropic PS, Fermi bubbles, and diffuse background template parameters are not constrained by the highlatitude fit, restricting these parameters to their highlatitude values does not significantly affect the results.
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 sourcecount functions and flux fractions are quoted with respect to the region within of the GC and , with no PSs masked. The sourcecount 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 bestfit sourcecount function for the NFW PS (dashed, orange), isotropic PS (dotted, green), and disk PS (solid, blue) populations. The diskcorrelated PS template accounts for the highflux 3FGL sources. Below photons/cm/s, the NFW PS template accounts for nearly all the PS emission; its sourcecount 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.
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 bestfit flux fraction for the NFW PS component is , while the bestfit 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 disklike 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 .
When the 3FGL sources are masked, the NPTF procedure yields a bestfit sourcecount function given by the orange band in the left panel of Fig. 3. Below the break, the sourcecount function agrees well with that found by the unmasked fit. In this case, the contributions from the isotropic and diskcorrelated 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, nonNFW 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 statisticslimited regime, the fit does not distinguish different models for the PS distribution at high significance,
To validate the analysis procedure, we generate simulated data using the bestfit parameters from the unmasked IG analysis; we include isotropic, Fermi bubbles, and Fermi p6v11 diffuse emission, as well as isotropic, disk and NFWdistributed PSs. The simulated data is then passed through the 3FGLunmasked 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 sourcecount 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 NFWdistributed 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 sourcecount 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 NFWdistributed PSs but does include NFW DM. The model parameters used to generate the simulated data are taken from the bestfit 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 sourcecount functions recovered for diskcorrelated and isotropic PSs are consistent with those used to generate the simulated data.
The sourcecount 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 sourcecount function seems to prefer an increasing below the break, implying most sources lie close to the cutoff luminosity, while previouslyconsidered sourcecount 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 PSdetection 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 sourcecount 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 PSdetection 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 DESC00012567, DESC0013999 and DESC0007968. B.R.S. and T.R.S. thank the Aspen Center for Physics, which is supported by National Science Foundation grant PHY1066293, for support during the completion of this work.
Note added in proof.—Recently, we became aware of Ref. [42], which studied the distribution of subthreshold PSs in the inner Galaxy using a wavelet technique.
Evidence for Unresolved GammaRay 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 PSdetection 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 templatefitting procedure, the photoncount probability distribution of observing photons in a given energy bin and pixel is assumed to be Poissonian,
(S1) 
where is the expected number of photons in the energy bin, summed over all the templates that contribute in that pixel. In particular,
(S2) 
where is the bestfit 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 energydependent 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 NFWdistributed DM component. We have repeated the standard templatefitting procedure and verified that our analysis pipeline reproduces results from those studies (see Sec. C.2 for details).
Unlike the standard templatefitting procedure, the NPTF does not assume that the photon counts are Poissonian. We compute the nonPoissonian photoncount probabilities using the method of generating functions. In particular, the total generating function for the photoncount probability distribution in each pixel is
(S3) 
where is an auxiliary variable, and and are the pixeldependent generating functions for the diffuse and nonPoissonian PS contributions, respectively. The probability of observing photons in pixel is then
(S4) 
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 isotropicbackground 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 wellconverged and is typically peaked around unity.The bestfit pixelaveraged intensity for the isotropic background in a given ROI is
(S5) 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 isotropicbackground template.
Parameter Prior Ranges High Latitude Inner Galaxy Table S1: Parameters and associated prior ranges for the highlatitude and IG analyses. Note that there can be up to three copies of the PS parameters when isotropic, diskcorrelated, and NFWdistributed 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 diffusebackground template is based on the Fermi p6v11 model for the majority of our analyses; however, we do perform a number of crosschecks with other diffuse models. As with the isotropicbackground template, the diffusebackground 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 linearflat prior, and the resulting bestfit value is typically close to unity. The bestfit pixelaveraged intensity is computed as in (S5). 
Fermi Bubbles
The Fermibubbles 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 bestfit value is again typically close to unity. The bestfit pixelaveraged intensity is computed as in (S5). 
NFW Dark Matter
The spatial dependence of the DM template is determined by , which gives the integrated photoncount intensity at the center of pixel from DM annihilation in the Galactic halo. The intensity is computed by the lineofsight integral(S6) where is the DM density profile and gives the distance from the GC. This work assumes a generalized NavarroFrenkWhite (NFW) density profile [35, 36]
(S7) 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 nonPoissonian statistics, which can be derived from the sourcecount function. In the majority of this work, the PS sourcecount function in a given pixel is modeled as a broken power law(S8) where is the break, are the slopes above and below the break, and is a pixeldependent 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 pixeldependent normalization of the sourcecount 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 sourcecount function parameters and an overall normalization parameter, , defined such that(S9) As with , the prior for is log flat and covers a broad range.
We also quote the pixelaveraged intensity , which is given by
(S10) with the expected number of counts in pixel arising from the PS population:
(S11) 
Isotropic Point Sources
The sourcecount 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 pixelaveraged intensity . 
Diskcorrelated Point Sources
The thindisk template is the projection along the lineofsight of the source spatial distribution given by:(S12) where , are cylindrical polar coordinates. The sourcecount 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 photoncount probability distribution can be written as a function of the 16 parameters
(S13) 
Then, for a data set consisting of the set of photon counts in each pixel , the likelihood function for observing a particular photoncount distribution over all pixels in the ROI is
(S14) 
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 modelevidence 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 frontconverting 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 CTBCOREcut 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,
The main body of the Letter focused primarily on two regions of interest: a highlatitude 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 .
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 photoncount probability distributions for the PS templates [31, 32]. To model the PSF for the CTBCOREcut data, we use the instrument response functions made available by [29] to approximate the PSF as a twodimensional Gaussian with energydependent width, (see Sec. D.4 for details).
Appendix B The FermiLAT PointSource 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 highenergy 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 frontconverting Ultraclean events and have applied further cuts on the CTBCORE parameter). The Fermi Collaboration has presented longitudeaveraged sensitivity estimates based on the integral flux from 0.1–100 GeV with three years of data, for sources with a pulsarlike 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 bestfit broken powerlaw 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 photonflux 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 sourcecount 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 highlatitude 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 bestfit 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 bestfit 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 highlatitudes (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.
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 3FGLmasked 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 nonPoissonian PS template to more accurately extract the contribution of unresolved PSs to the total EGB.
The central analyses in this work assume the sourcecount function parameterization in (S8), however we also consider the effect of forcing the sourcecount function to have a sharp cutoff at high flux, while still allowing a break at lower flux. For example, we can impose a highflux cutoff consistent with the flux of the brightest 3FGL PS in the ROI: photons. In this case, the NPTF gives the following bestfit 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 3FGLmasked 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 sourcecount function without the highflux 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 wellconstrained within the prior ranges, with one exception.
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 crosscheck 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.
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 (fullsky) 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 doubleexponential thindisk 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 thickdisk 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 PSlike 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.
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 NorthSouth symmetry of the excess by masking each hemisphere in turn. The results for the fullsky and hemisphere analyses are similar, so only the latter are presented here. Additionally, we only show results for the 3FGLunmasked 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 bestfit sourcecount function and flux fractions for the Northern and Southern hemispheres in the top and bottom rows, respectively. The sourcecount 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.
d.2 Varying the Diffuse Model
A potentially serious source of systematic uncertainty is due to limitations in modeling the diffuse gammaray background arising from the propagation of highenergy cosmic rays in the Galaxy. The primary contributions come from bremsstrahlung of highenergy cosmic rays passing through the interstellar gas, inverse Compton scattering of Cosmic Microwave Background, interstellar, and infrared radiation off highenergy electrons, and the decay of boosted pions produced in inelastic proton collisions with the interstellar gas. Modeling the cosmicray 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 diffusebackground models can help to quantify the effects of mismodeling 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 largescale 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 datadriven p7v6 diffuse model having absorbed part of the GeV excess.
The left panel of Fig. S4 shows the bestfit sourcecount 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.
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(BV). 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 diffusemodel prior ranges for the p6v11 and p7v6 diffuse models; in all cases, the three component templates are wellconverged within the prior ranges. For brevity, we only show results for the 3FGLmasked analyses. Additionally, for computational reasons, we only include a nonPoissonian 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 sourcecount 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 sourcecount function, the NPTF consistently finds a nonzero 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 gascorrelated emission (pion decay and bremsstrahlung) within the innermost Galactocentric ring. The addition of this template does not significantly alter the sourcecount function or flux fraction attributed to the NFW PSs.
d.3 Scan Along the Galactic Plane
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 mismodeled 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 bestfit sourcecount 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 sourcecount 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 innerslope value used by [12], we reproduce the (slightly smaller) flux fraction of the excess found by those authors.
d.4 Point Spread Function
An accurate PSF is a crucial component of the NPTF procedure, as an incorrect parametrization can lead to over or underestimation of the PS component. The instrument response functions for the Fermi CTBCOREcut 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 twodimensional Gaussian
(S15) 
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.89311.943 GeV). As described earlier, the isotropic, Fermi bubbles, and NFW DM templates are smeared with this PSF. (The diffusebackground template is smeared with the exact energydependent PSF using the Fermi Science Tools.) The PSF must also be properly accounted for in the generatingfunction formalism that is used to obtain the nonPoissonian likelihood function (see [32] for further details).
It is known that the real PSF has powerlaw tails that are not captured by (S15) [29]. One might rightfully be concerned that ignoring these powerlaw tails can lead to mischaracterization of the PS population. To illustrate the effect of mismodeling the PSF, Fig. S8 shows the bestfit sourcecount function for the IG analysis assuming extreme values for the GaussianPSF 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 sourcecount distribution. This effectively shifts the sourcecount 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 deltafunction 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 highflux PSs, while at the same time predicting a larger population of unresolved sources.
The fact that the energyaveraged PSF used in Fig. 2 results in a sourcecount function that matches the highflux observed PSs reasonably well gives us confidence in our treatment of the PSF. Additionally, we have illustrated here that extreme variations to the PSFparameter 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 templatefitting methods, such as [12, 14], have scanned over this parameter and found bestfit values in the range between and . Variations in have negligible effects on our final result. Fig. S9 shows the bestfit sourcecount function in the inner 10 with , assuming and (left and right panels, respectively). As is evident, the sourcecount 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.
d.6 DiffuseCorrelated PointSource Template
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 gammaray 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 diffusecorrelated (diffcorr) PS template is not only sensitive to localized structure in the diffuse model, but also to the presence of a disklike population of PSs. For example, if the unresolved PSs are preferentially located close to the plane, then they may be absorbed by the diffcorr 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 bestfit sourcecount functions obtained when doing the NPTF in the IG region (without masking the 3FGL sources), including a diffcorr PS template in addition to diskcorrelated and NFW PS templates, along with the standard Poissonian templates. Here, the diffcorr PS template only extends to from the GC. To simplify the analysis, we do not include the subdominant isotropic PS template. The bestfit flux normalization for the diffuse background is consistent with that obtained from the highlatitude analysis. In addition, the bestfit NFW DM normalization is consistent with zero, and the recovered sourcecount function parameters for disk and NFW PSs are consistent with those found when not including the diffcorr PS template, since the diffcorr PS template does not absorb a significant fraction of the flux.
As a point of comparison, we repeat the procedure letting the diffcorr PS template have support over the full IG region. Now, there is a potential degeneracy between the diffcorr PS template, the diskcorrelated PS template, and the diffuse template. To break some of this degeneracy, we fix the diskcorrelated PS parameters to their bestfit values from the scan including diskcorrelated, isotropic, and NFW PSs. These results are shown in the right panel of Fig. S11. Again, the bestfit NFW DM normalization is consistent with zero, however the sourcecount function for the NFW PSs changes. In particular, the sourcecount function for NFW PSs is shifted to lower flux, potentially suggesting that some of the nearthreshold sources could either be more disklike in morphology or associated with mismodeling 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 diffcorr PS template, here the bestfit normalization of the diffusebackground template is lower than its bestfit value at high latitudes. When the NFW PS template is not included, the normalization can be shifted down by as much as 20%; the bestfit 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 bestfit normalization of the diffuse background approaches its value from high latitudes, and the flux absorbed by the diffcorr PS template decreases.
d.7 Binned SourceCount Functions
The imposition of a broken power law for the sourcecount function might be overconstraining in some cases. For example, it could lead to the apparent exclusion of a DM component because the extrapolation to low flux of the sourcecount function from highflux sources is already too large. Furthermore, the brokenpower law prescription may yield misleading results for the true uncertainty on the sourcecount 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 logarithmicallyspaced flux bins shown in Fig. S12. Within each bin, is constant as a function of . The sourcecount function model parameters are seven normalization parameters, one for each bin, which are taken to have logflat priors over the range indicated for in Tab. S1.
Using the binned sourcecount function, we perform the NPTF on the 3FGLmasked IG ROI. For simplicity, we leave out isotropic and diskcorrelated PS templates from this analysis. This is justified by the fact that leaving out these two templates from the broken powerlaw 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 sourcecount function broadly consistent with the original analyses assuming broken power laws. The blue (pink) bands indicate the 68% (95%) confidence intervals for the sourcecount 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 powerlaw 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 lowflux bins are not wellconstrained by the fit, so the posterior distributions for these parameters are heavily influenced by the logflat prior distributions. In the broken powerlaw fit, only the overall normalization parameter was taken to have a logflat prior range, while in the binned fit all PS parameters have logflat prior ranges. This point, combined with the fact that the broken powerlaw is more constrained than the binned sourcecount function, leads to larger uncertainties at low (and, to some extent, high) flux in the binned analysis than in the broken powerlaw 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 powerlaw results, within uncertainties. Similarly, the binned fit prefers a large number of additional unresolved sources with fluxes only slightly below the PSdetection 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
The bestfit sourcecount function recovered by the NPTF for the unresolved NFWcorrelated sources is quite different to those previously considered in the literature—see e.g., [32, 41, 23, 24, 19]. With our sourcecount 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 sourcecount function using a luminosity function derived from observed MSPs in the nearby field of the Milky Way [41], for both NFWdistributed sources and the thickdisk 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 thickdisk model has been normalized to produce the correct number of highlatitude () 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 NFWdistributed sources have instead been normalized to match the flux attributed to the NFW PS template in our analysis.
Fig. S13 shows the expected sourcecount functions, averaged over our standard ROI, for these scenarios, for the 3FGLunmasked fit (Fig. 2). In agreement with [23], we find that for NFWdistributed 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 dotdashed line in Fig. S13 shows the predicted sourcecount function for the thickdisk distribution. It is remarkably similar to the bestfit sourcecount function for the thindisk PS model extracted from the data. In particular, we estimate the slope of the purple dotdashed line at low flux in Fig. S13 to be , while the NPTF predicts the slope of the thindisk PS sourcecount 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 sourcecount function for the NFW PSs rises sharply above the red curves at fluxes – photons/cm/s. As this sourcecount function is very shallow below the flux threshold, the total number of sources is dominated by this relatively highflux 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].
Appendix F Survival Function
As a further crosscheck that the PS identification is working selfconsistently, we attempt to identify pixels that are likely to contain unresolved sources in the 3FGLmasked sky. To do so, we perform a standard template fit in the ROI (excluding the PS templates) and determine the bestfit 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
(S16) 
where denotes the cumulative probability of observing counts for a Poisson function with mean . For example, if the bestfit 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.