# Deciphering Contributions to the Extragalactic Gamma-Ray Background from 2 GeV to 2 TeV

## Abstract

Astrophysical sources outside the Milky Way, such as active galactic nuclei and star-forming galaxies, leave their imprint on the gamma-ray sky as nearly isotropic emission referred to as the Extragalactic Gamma-Ray Background (EGB). While the brightest of these sources may be individually resolved, their fainter counterparts contribute diffusely. In this work, we use a recently-developed analysis method, called the Non-Poissonian Template Fit, on up to 93 months of publicly-available data from the *Fermi* Large Area Telescope to determine the properties of the point sources that comprise the EGB. This analysis takes advantage of photon-count statistics to probe the aggregate properties of these source populations below the sensitivity threshold of published catalogs. We measure the source-count distributions and point-source intensities, as a function of energy, from 2 GeV to 2 TeV.
We find that the EGB is dominated by point sources, likely blazars, in all seven energy sub-bins considered. These results have implications for the interpretation of IceCube’s PeV neutrinos, which may originate from sources that contribute to the non-blazar component of the EGB.
Additionally, we comment on implications for future TeV observatories such as the Cherenkov Telescope Array. We provide sky maps showing locations most likely to contain these new sources at both low ( GeV) and high ( GeV) energies for use in future observations and cross-correlation studies.

MIT/CTP-4808

## 1 Introduction

The Extragalactic Gamma-Ray Background (EGB) is the nearly isotropic all-sky emission that arises from sources outside of the Milky Way. The OSO-3 [37, 65] and SAS-2 satellites [50, 49] were the first to see hints of the EGB and have since been followed by EGRET [87, 94] and, most recently, the *Fermi* Large Area Telescope^{1}

Most recently, *Fermi* presented a measurement of the EGB intensity from 100 MeV to 820 GeV [15]. The total EGB intensity is the sum of all resolved point sources (PSs) and smooth isotropic emission. The smooth emission, referred to as the Isotropic Gamma-Ray Background (IGRB), arises from PSs that are too faint to be resolved individually as well as other truly diffuse processes. It is also important to note that both the EGB and IGRB may be contaminated by cosmic rays that are mis-identified as gamma rays; this emission is expected to be smoothly distributed across the sky.
Of the known gamma-ray emitting PSs at high latitudes, which are captured by Fermi’s 3FGL [9] catalog from 0.1–300 GeV and the more recent 2FHL [16] catalog from 50–2000 GeV, the dominant source class is blazars.

In this work, we use a novel analysis method, called Non-Poissonian Template Fitting (NPTF), to study the source populations that contribute to the EGB in a data-driven manner. The method relies on photon-count statistics to illuminate the aggregate properties of a source population, even when its constituents are not individually resolvable [72, 67, 68]. This allows us to constrain the contribution of PSs to the EGB whose flux is too dim to be detected individually. While at very low fluxes the NPTF also loses the ability to distinguish PSs from smooth emission, the threshold for PS detection is lower for the NPTF than it is for other techniques that rely on finding individually-significant sources. This is because the NPTF only measures the aggregate properties of a PS population.

Using the NPTF, we are able to recover, for the first time, the source-count distribution (*e.g.*, flux distribution) for isotropically distributed PSs at high Galactic latitudes, as a function of energy from 1.89 GeV to 2 TeV. This builds on previous studies that use related methods to obtain the source-count distributions in single energy bins from 2–12 GeV [104, 68] and from 50–2000 GeV [17].

The source-count distribution for a given astrophysical population convolves information about its cosmological evolution. For a flat, non-expanding universe, a uniformly distributed population of galaxies has a differential source-count distribution , where is the source flux at Earth and is the differential number of sources [84]. This is the well-known Euclidean limit. However, the power-law index changes when one takes the standard CDM cosmology and more realistic assumptions for the redshift evolution of source-dependent observables such as luminosity. Therefore, the features of the source-count distribution—especially, its power-law indices and/or flux breaks—encode information about the number of source classes contributing to the EGB as well as their cosmological evolution.

These source-count distributions provide the keys for interpreting the GeV–TeV sky. For example, they enable us to obtain the intensity spectrum for PSs, down to a certain flux threshold, as a function of energy. We find that while the EGB is dominated by PSs, likely blazars, in the entire energy range from 1.89–2000 GeV, there is also room for other source classes, which contribute flux more diffusely, to produce a sizable fraction of the EGB.
Our findings may therefore leave open the possibility that IceCube’s PeV neutrinos [3, 2, 4, 5] can be explained by hadronic interactions in *e.g.*, SFGs [75, 96, 24] or mAGN [61], which—as we show in Sec. 3—show up as smooth isotropic emission under the NPTF. Additionally, the high-energy source-count distributions allow us to make predictions for the number of blazars, which dominate the high-energy data, that will be resolved by upcoming TeV observatories such as the Cherenkov Telescope Array (CTA) [38, 45].
While our analysis does not let us conclusively identify the locations of these sources, we provide maps showing the locations on the sky where, statistically, there are most likely to be PSs.

This paper is organized as follows. We begin in Sec. 2 by reviewing the analysis methods. Sec. 3 then applies these methods to simulated sky maps. We cannot stress the importance of these simulated data studies enough; they are crucial for proving the stability of the analysis methods and laying the foundation for the data results that follow. Our data study is divided into two separate analyses for low (1.89–94.9 GeV) and high (50–2000 TeV) energies, described in Sec. 4 and 5, respectively. The global fits to the full energy range, as well as their implications, are discussed in Sec. 6. Further details on the creation of the simulated data maps and supplementary analysis plots are provided in the Appendix. The main results of this work are summarized in a few key figures. In particular, the source-count distributions for the low and high-energy analyses are shown in Figs. 4, 6, and 8, respectively, while Fig. 9 presents a spectral fit to the PS intensity from 2 GeV to 2 TeV.

## 2 Methodology

In this paper, we make use of both Poissonian and non-Poissonian template-fitting techniques.
Poissonian template fitting is a standard tool in astrophysics for decomposing a sky map into component “templates” with different spatial morphologies. The NPTF builds upon this technique by allowing for the addition of templates whose spatial morphology traces the distribution of a PS population, even if the exact position of the sources that make up that population are not known. More precisely, in both template-fitting procedures one starts with a data set that consists of counts in each pixel .^{2}

(1) |

where denotes the probability of observing photons in pixel with model parameters .

In Poissonian template fits, the probabilities are Poisson distributions, with the model parameters only determining the means of the distributions. That is, the mean expected number of photon counts at each pixel may be written as

(2) |

where the sum is over template components and denotes the mean of the component for model parameters . The may parameterize, for example, the overall normalization of the templates or the shapes of the templates. Then, the probability is simply given by the Poisson distribution with mean .

In the NPTF, the situation is more complicated because we do not know where the PSs are. As a result, if we want to calculate the probability of observing photons in a given pixel , we must first calculate the probability that a PS (or a collections of PSs) exists in the vicinity of the pixel , with a given flux (or set of fluxes). Then, for that PS population, we calculate the probability of photons being produced in pixel . Convolving these two calculations together leads to distinctly non-Poissonian probabilities. In particular, the probability distributions in the presence of unresolved PSs tend to be broader than Poisson distributions, if both distributions have the same mean expected number of photon counts. The intuition behind this fact is that relative to a diffuse source, a collection of PSs leads to more “hot” pixels with many photons (where there are PSs) and more “cold” pixels with very few photons (where there are no PSs).

### 2.1 The Templates

We include three Poissonian templates for (1) diffuse gamma-ray emission in the Milky Way, assuming the *Fermi* p8r2 (*gll_iem_v06.fits*) foreground model, (2) uniform emission from the *Fermi* bubbles [95], and (3) smooth isotropic emission. Each of these templates is associated with a single model parameter describing its overall normalization. Variations to the choice of foreground model and bubbles template will be discussed in Sec. 4.2.

The model parameters specific to the isotropic-PS population enter into the source-count distribution , which we characterize as a triply-broken power law: {widetext}

(3) |

In particular, there are three breaks, , along with four indices, , and the overall normalization, .^{3}

It is important to stress that the photon-count probabilities are non-Poissonian in the presence of unresolved PSs because their locations are unknown. Once we know where a PS is, we can fix its location and describe it through a Poissonian template with a free parameter for the overall normalization of the source. However, even resolved sources with known locations may be characterized by the non-Poissonian template if we do not also put down Poissonian templates at their locations. This is the approach that we take throughout this work; that is, we model both the resolved (in the 3FGL and 2FHL catalogs) and unresolved PS populations through a single distribution, without individually specifying the locations of any sources.

The point-spread function (PSF) must be properly accounted for in the template-fitting procedure. The diffuse models are smoothed according to the PSF using the Fermi Science Tools routine gtsrcmaps. The bubbles template is smoothed with a Gaussian approximation to the PSF, with width set to give the correct 68% containment radius in each energy bin. We follow the prescription developed in [72] to account for the PSF in the calculation of the non-Poissonian photon-count probabilities; for this, we use the King function parameterization of the PSF provided with the instrument response function for the given data set. In Sec. 4.2, however, we show that consistent results are obtained when using a Gaussian approximation to the PSF instead.

### 2.2 Bayesian Fitting Procedure

The formalism developed in [72, 67, 68] (see also [104] and [69]) is used to calculate the photon-count probability distributions in each pixel as a function of the Poissonian and non-Poissonian model parameters . Then, Bayesian techniques are used to construct a posterior distribution for the parameters and the likelihood function in (1). We construct the posterior distribution numerically using the MultiNest package [48, 34] with 700 live points, importance nested sampling and constant efficiency mode disabled, and sampling efficiency set for model-evidence evaluation.

All prior distributions are taken to be flat except for , which is taken to be log-flat. The prior ranges for the model parameters are shown in Tab. 1.

Parameter | Prior Range | Parameter | Prior Range | Parameter | Prior Range |
---|---|---|---|---|---|

[-10, 20] | |||||

ph | |||||

ph | |||||

ph |

These prior ranges successfully reconstruct the source-count distributions of simulated data sets, as discussed in Sec. 3. Variations to the prior ranges in Tab. 1 are considered in Sec. 4.2.

In Tab. 1, the parameter denotes the normalization of the template, which is defined in terms of a baseline value. The baseline value is obtained by first performing a Poissonian template fit over 17 (10) log-spaced energy sub-bins between and GeV (50 and 2000 GeV) for the low (high)-energy analysis. When this procedure is applied to the low-energy analysis where the known PSs are very bright, we mask the 300 brightest and most variable 3FGL sources, at 95% containment. At both high and low energies, we include a PS model constructed from the 3FGL catalog.^{4}

The actual energy bins used for the NPTF studies presented in this work are larger than the sub-bins described above. Therefore, the baseline normalizations used to define the NPTF priors in the energy range are found by applying the best-fit Poissonian normalizations from the individual sub-bins to the corresponding templates, which are then combined.^{5}

### 2.3 Exposure Correction

While the source-count distribution is defined in terms of flux, , with units of , the priors for the breaks in Tab. 1 are written in terms of counts, . To convert from flux to counts, we multiply by the exposure of the instrument, with units of cm s. However, the relation between flux and counts is complicated by the fact that the exposure of the instrument varies both with energy and position in the sky. Below, we describe how we deal with both complications, starting first with the energy dependence.

The exposure map in the energy sub-bin is given by . To construct the exposure map in the larger energy range from , which contains multiple energy sub-bins, we average over the of the individual sub-bins, weighted by a power-law spectrum , as this is generally consistent with the isotropic spectrum over most of our energy range.
This procedure introduces a source of systematic uncertainty in going from counts to flux, as not all source components have an energy spectrum consistent with this spectrum. However, we have checked that variations to this procedure—such as weighting the exposures in the sub-bins by power laws of the form , with varying between and —do not significantly change the results.^{6}

The breaks in Tab. 1, with units of counts, are defined relative to the mean exposure , averaged over all pixels in the region of interest (ROI). Because the NPTF is performed at the level of counts and not flux, we must also convert the source-count distribution to a distribution , which is unique to each pixel :

(4) |

Then, the photon-count probability distribution must be computed uniquely at each pixel. In practice, however, it is numerically expensive to perform this procedure for every pixel in the ROI. Instead, we follow [104] and break the ROI up into regions by exposure. Within each region, we assume that all pixels have the same exposure, which is taken to be the mean over all pixels in the sub-region. The likelihood function is then computed uniquely in each exposure region, and the total likelihood function for the ROI is the product of the likelihoods across exposure regions. In practice, we find that our results are convergent for . We will take throughout this work, though we have checked that our main results are consistent with those found using .

### 2.4 Data Samples

We run the NPTF analysis, as described above, on *Fermi* data, considering low (1.89–94.9 GeV) and high (50–2000 GeV) energies separately. The former is discussed in Sec. 4, while the latter is the focus of Sec. 5. The primary difference between the data sets used in these studies is the data-quality cuts; moving to higher energies requires loosening these criteria to avoid being limited by statistics. The overlap in energy between the two studies allows us to compare the consistency of the results when transitioning between analyses.

The low-energy study uses the Pass 8 *Fermi* data from August 4, 2008 to June 3, 2015. The primary studies use the top quartile of the *ultracleanveto* event class (PSF3) as ranked by angular resolution, although the top-three quartiles (PSF1–3) are also studied separately.^{7}*ultracleanveto* event class is the cleanest event class released with the Pass 8 data and is recommended for studies of the EGB. However, the source event class has an enhanced exposure and thus may be advantageous at high energies where statistics become limited. On the other hand, we expect the source data to have additional cosmic-ray contamination relative to the *ultracleanveto* data.

The recommended^{8}*nside*=128, which corresponds to pixels roughly to a side. We consider four separate energy bins: , , , and GeV.

In the low-energy analysis with ultracleanveto PSF3 data, the means of the weighted exposure maps in the four increasing energy bins are cm s over the region of interest with . The 68% containment radii for the PSF, averaged over the isotropic spectra in the energy sub-bins, are degrees. Going to PSF1–3 data, the exposures increase to cm s, while the 68% containment radii of the PSF degrade to degrees. Going to source data with PSF1–3, the exposures ( cm s) increase further, while the 68% containment radii ( degrees) are essentially the same as in the ultracleanveto case.

The high-energy analysis uses the Pass 8 *Fermi* data from August 4, 2008 to May 2, 2016 and all PSF quartiles of either the *ultracleanveto* or *source* event class. The ROI is also extended to . We include more data in the high-energy analysis as there are far fewer photons than at lower energies. We employ the recommended event-quality cuts as in the low-energy analysis and also choose *nside*=128 HEALPix pixelation. Results are presented for the three energy bins , and GeV. With ultracleanveto data, the weighted exposures in the energy bins are cm s, while with source data the exposures become cm s. For both data sets, the 68% containment radii are approximately degrees. We will also discuss results of analyses performed over a single wide-energy bin from GeV.

## 3 Simulated Data Studies

To study the behavior of the NPTF, we apply it to simulated data sets of the gamma-ray sky. These results are crucial both for understanding systematics associated with the NPTF as well as for interpreting the results of the NPTF in terms of evidence for or against the existence of these source populations.

A simulated data map can be created starting from a particular source population that contributes to the EGB. Using a theory model for the energy spectrum and luminosity function, the source-count distribution for that population can be derived in a specified energy range—see Appendix A.1 for further details on this procedure. The appropriate number of sources is then drawn from this function and randomly distributed across the sky, with counts chosen to follow the intensity spectrum. Sources are then smeared with the appropriate Gaussian PSF to mimic the desired Fermi data set
bin-wise in energy, and Poisson counts are drawn to obtain the simulated map for the population. This is then combined with the simulated contribution of the p8r2 foreground model and the *Fermi* bubbles, whose normalizations are determined from the Poissonian template fits to the real data, as described in Sec. 2.

For most of this section, we simulate data corresponding to the PSF3 event type (best PSF quartile) of the *ultracleanveto* event class and focus on the following four energy bins: [1.89, 4.75], [4.75, 11.9], [11.9, 30], and [30, 94.9] GeV.^{9}

### 3.1 Blazars

Active galactic nuclei (AGN) are the highly luminous central regions of galaxies where emission is dominated by accretion onto a supermassive black hole [100]. If the black hole is spinning, then relativistic jets may also form. Blazars are a subclass of AGN in which the jet is oriented within of the line-of-sight [25]. The spectral energy distribution of these objects is bimodal with a peak in the ultraviolet due to synchrotron radiation of electrons in the jet, and another peak in the gamma band from inverse Compton scattering of the same electrons [54, 56, 57, 31]. There is also the possibility of a hadronic contribution to blazar gamma-ray spectra, although this is likely to be sub-dominant [97, 35, 102]. Blazars may be further classified as either BL Lacertae (BL Lacs) or Flat Spectrum Radio Quasars (FSRQs), which are characterized by the absence or presence of broad optical/ultraviolet emission lines, respectively.

Before *Fermi*, few blazars had been identified in gamma rays, and to estimate the size of this population, one had to extrapolate based on those observed at lower frequencies. However, *Fermi* brought the discovery of many more blazars in the gamma-ray band, making it possible to study their properties directly [1, 20, 21, 43, 58, 81]. Most recently, 403 blazars (with ) from the First LAT AGN Catalog [7] were studied [22]. FSRQs and BL Lacs were considered together in the same sample to improve statistics. We use the best-fit luminosity and spectral energy distributions given in [22] (specifically, the luminosity-dependent density evolution, or LDDE, scenario) to model the blazar component in our simulated data and refer to it as the “Blazar–1” model. Alternatively, we also consider BL Lacs and FSRQs separately, adding up their respective contributions using the LDDE1 model from [21] and the LDDE model from [20], which we refer to as the “Blazar–2” model. This model predicts a much flatter source-count distribution below the *Fermi* detection threshold, with more low-flux sources. The two source-count models approximately bracket the current theoretical uncertainty in the faint-end slope of blazars, and we use them to study the response of different blazar models to the NPTF, although this is meant to be purely illustrative and by no means exhaustive.

Figure 1 shows the best-fit source-count distributions recovered when the NPTF analysis is run on the Blazar–1 simulated data map, assuming the PSF3 instrument response function. In each panel, the dark (light) red band is the 68% (95%) credible interval for the isotropic-PS source-count distribution as recovered from the posterior and constructed pointwise in flux. The red line shows the median source-count distribution, constructed in the same way. The dashed red curve, on the other hand, indicates the source-count distribution of the blazar model used to generate the simulated data. A flux histogram of the simulated PSs for the particular realization shown here is given by the red points, with vertical error bars indicating the 68% credible interval associated with Poisson counting statistics on the number of sources in that bin. Notice that these error bars become large at high fluxes because there are very few sources per flux bin.

In general, the reconstructed source-count distribution is in good agreement with the input source-count distribution at intermediate fluxes, with uncertainties becoming large at low and high fluxes. At high flux, this is due to the fact that it is unlikely to draw a bright source from the underlying source-count distribution. At low fluxes, it is difficult to distinguish PS emission from genuinely isotropic emission. To illustrate this point, we also mark the flux that corresponds to a single photon on average (in the particular energy range, region-of-interest, and event class) with the vertical dot-dashed black line. At fluxes corresponding to counts near or below 1 photon, it is difficult to distinguish PS emission from smooth emission with the NPTF, as evidenced by the growing uncertainties. In this low-flux regime, we do not expect that the NPTF will be able to fully recover the properties of the input source-count distribution.

The vertical dotted green lines in Fig. 1 correspond to the fluxes above which 90%, 50%, and 10% (from left to right) of the photon counts are accounted for, on average, by sources with larger flux. Note that in the lowest energy bin, 90% of the flux arises from sources that contribute more than one photon. Moving towards higher energies, a larger fraction of the flux arises from sources that contribute less than one photon. In all energy bins, more than 50% of the flux is accounted for by sources that contribute more than a single photon each.

The corresponding energy spectra for the various templates are shown on the top left panel of Fig. 2. As is evident, these blazars show up as PSs under the NPTF; indeed, the smooth isotropic flux (blue) is sub-dominant in each energy bin. Overlaid in dashed red is the spectrum for the simulated Blazar–1 sources. The sum of the smooth and PS isotropic components—which is simply the EGB intensity—is consistent with the simulated spectrum for the blazar model. The green curve shows the median of the posterior for the galactic diffuse model spectrum. The energy spectrum of the diffuse model is softer than that for blazars, so that the diffuse model dominates more at low energies than at high. The sum of the components (yellow band) is consistent with the total flux in the simulated data (black lines) at 68% confidence.

As a contrasting example, we also simulate the Blazar–2 model, which predicts more low-flux sources than the previous example we considered. The best-fit source-count distributions for the Blazar–2 simulated maps are shown in Fig. LABEL:fig:bl2dnds. Once again, we see good agreement between the input data and the recovered source-count distribution above the single-photon sensitivity threshold. In this case, however, the reference model predicts a larger fraction of flux coming from sources below this threshold. For example, about 50% of the flux comes from sub-single photon sources in the second energy bin, and this fraction only increases further at higher energies. The corresponding energy spectrum is shown in the top right panel of Fig. 2. As expected, an increasing amount of flux is absorbed by the Poissonian isotropic template. However, the EGB spectrum, shown by the purple band, is still consistent with the input spectrum for the Blazar–2 model.

To further quantify the ability of the NPTF to reconstruct the blazar flux as PS emission, it is convenient to consider the ratios in each energy bin, where is the PS intensity found by the NPTF and is the blazar intensity in the simulation. For the Blazar–1 model, we find^{10}

in each of the four respective energy bins, while for the Blazar–2 model, we find

for the particular Monte Carlo realizations shown.^{11}

Next, we explore whether including more quartiles of the ultracleanveto data, as ranked by PSF, increases our ability to reconstruct the blazar flux as PSs under the NPTF. When including more quartiles of data, there are two competing effects that determine our ability to constrain the PS flux: on the one hand, we increase the effective area, but on the other hand, we worsen the angular resolution of the data set. We investigate these effects by repeating the Monte Carlo tests described above using the PSF1–3 instrument response function. The results of the PSF1–3 tests are described in Appendix A.3, and here we simply quote the fractions

for a generic realization of the Monte Carlo simulations for the Blazar–2 model. The PSF1–3 event type increases our ability to distinguish between the blazar emission and smooth emission compared to the PSF3 event type.

### 3.2 Star-Forming Galaxies

Star-forming galaxies (SFGs) like our own Milky Way are individually fainter, though much more numerous, than blazars. The modeling of SFGs in the gamma-ray band is highly uncertain, as *Fermi* has only detected eight SFGs thus far [53]. However, SFGs could still contribute a sizable fraction of the total flux observed by *Fermi*. Even though SFGs are PSs, their flux is expected to be dominated by a large population of dim sources degenerate with smooth isotropic emission. Under the NPTF, therefore, we expect that the majority of their emission will be absorbed by the smooth isotropic template.
To illustrate this point, we simulate SFGs using the luminosity function and energy spectrum from [96]. In that work, input from infrared wavelengths was used to construct a model for the infrared flux from SFGs.
Then, a scaling relation was used to convert from infrared to gamma-ray luminosities. The contributions from quiescent and starburst SFGs were considered separately, along with SFGs that host an AGN. Note, however, that other models predict less emission from SFGs than this particular case—see *e.g.*, [71, 63, 12].

The results of the SFG-only simulations are described in Appendix A.4. We find that while the NPTF does detect a small PS component in the first few energy bins, as the result of a few SFGs above the sensitivity threshold of the NPTF in those energy bins, by far most of the SFG emission is detected as smooth isotropic emission, with the ratio in all energy bins, where is the intensity of smooth isotropic emission. Moreover, the intensity is consistent with the simulated EGB (SFG flux) in all energy bins, at 68% confidence.

### 3.3 Blazar and SFG combination

A perhaps more realistic scenario for testing the NPTF is to consider a scenario where both SFGs and blazars contribute to the EGB. Therefore, we create simulated maps that include both components and test them on the NPTF. The recovered energy spectra for the SFG + Blazar–1 (Blazar–2) example is shown in the bottom left (right) panel of Fig. 2. In both cases, the PS spectrum is consistent with that found in the blazar-only simulations, which are shown in the top panels in that figure. The reconstructed source-count distributions for these examples are not shown, as they are consistent with those found in the blazar-only cases.

In the case of of the Blazar–1 model, the spectra of the smooth isotropic emission and the PS emission trace the spectra of the input SFG population and blazar population, respectively. In the case of the Blazar–2 model, the PS flux is further below the input blazar spectrum, as was found in the blazar-only simulations. However, the smooth isotropic emission is further above the simulated SFG spectrum. In both cases, the sum of the smooth isotropic emission and PS emission (EGB) is consistent with the simulated blazar plus SFG flux.

There is, in fact, a subtle difference between the PS distribution recovered with and without the addition of a SFG population. The difference becomes noticeable when comparing the fractions

for SFG + Blazar–1 and

for SFG + Blazar–2 to the corresponding values for the blazar-only simulations. In the simulations with SFGs, the fractions are generally higher and have larger uncertainties. The reason for this is that the SFG emission is degenerate with an enhanced sub-threshold component to the PS source-count distribution.

Simulating data with the PSF1–3 instrument response function, we find that the ratios are somewhat closer to unity than in the PSF3 case. In particular, for the SFG + Blazar–2 model simulations,

The improved exposure allows the NPTF to probe lower fluxes and to therefore recover a larger fraction of the isotropic-PS emission.

## 4 Low-Energy Analysis: 1.89–94.9 GeV

The findings from the previous section illustrate that the NPTF procedure is able to set strong constraints on the PS (*e.g.*, blazar) and smooth Poissonian (*e.g.*, SFGs, mAGN) contributions to the EGB.
In this section, we focus on the energy range from 1.89–94.9 GeV, and begin by presenting the results of our benchmark analysis on the real *Fermi* data. This is followed by a detailed discussion of potential systematic uncertainties and their effects on the conclusions.

### 4.1 Pass 8 ultracleanveto Data

#### Top PSF Quartile

We begin by analyzing the ultracleanveto PSF3 data for , using the p8r2 foreground model. This is referred to as the “benchmark analysis” throughout the text. Table 2 provides the best-fit intensities for each spectral component, as a function of energy, and the best-fit spectra are plotted in the left panel of Fig. 3. The p8r2 diffuse model is shown in green (median only), while the smooth isotropic and isotropic-PS posteriors are shown by the blue and red bands, respectively. The best-fit spectrum for PSs with in the 3FGL catalog [9] is shown by the dashed black line in Fig. 3; the spectrum as plotted should be treated with care as systematic uncertainties are not properly accounted for. In particular, the 3FGL catalog includes sources between 0.1–300 GeV. At the high end of this range, the spectrum is driven to a large extent by extrapolations from lower energies, where the statistics are better. The potential errors associated with such extrapolations are difficult to quantify and are not shown in Fig. 3.
As a result, a direct comparison between the 3FGL spectrum and our results is difficult to make, especially in the highest energy bins. For this reason, we have a dedicated NPTF study for energies greater than 50 GeV in Sec. 5. Those results are compared to the *Fermi* 2FHL catalog [17], which is explicitly constructed at higher energies and is likely a more faithful representation of above-threshold PSs in this regime.

Energy | |||||
---|---|---|---|---|---|

1.89–4.75 | |||||

4.75–11.9 | |||||

11.9–30.0 | |||||

30.0-94.9 |

*Fermi*bubbles template intensity is defined relative to the interior of the bubbles, while the intensities of the other templates are computed with respect to the region . The best-fit EGB intensity, which is the sum of the smooth and PS isotropic contributions, is also shown. The posterior distributions for the template intensities are provided in Fig. S5–LABEL:fig:p8triangle4.

Energy | |||||||
---|---|---|---|---|---|---|---|

1.89–4.75 | |||||||

4.75–11.9 | |||||||

11.9–30.0 | |||||||

30.0-94.9 |

The source-count distributions reconstructed from the NPTF are shown in Fig. 4, with best-fit parameters provided in Tab. 3. For comparison, the binned 3FGL source-count distributions are also plotted; the vertical error bars represent 68% statistical uncertainties and do not account for systematic uncertainties. A few trends are clearly visible. First, each flux break tends to have large uncertainties. This may be a reflection of the fact that the real source-count distribution is not a simple triply-broken power law, but rather a more complicated function, as in the blazar simulations of Sec. 3. Therefore, the best-fit values for each of these parameters, when viewed independently, may be somewhat deceptive. As is evident in Fig. 4, the posteriors for the breaks and indices are distributed in such a way as to describe a smooth concave function for .

At very high and very low flux, the uncertainties on the indices ( and , respectively) become large. At high flux, this is simply due to the fact that there are very few sources, so the source-count distribution falls off rapidly. At low flux, the large uncertainties on arise from the difficulty in distinguishing the isotropic-PS contribution from its smooth counterpart. Indeed, below the single-photon boundary (dot-dashed black line), the NPTF analysis starts to lose sensitivity. The posterior distributions for the slopes above (below) the highest (lowest) break are highly dependent on the priors and so the quoted values in Tab. 3 should be treated with care.

The presence of any distinctive breaks encodes information about the number of source populations as well as their evolutionary properties. In all energy bins, we see that the NPTF places the lowest break, , close to the one-photon sensitivity threshold and the highest break, , in the vicinity of the highest-flux 3FGL source (see Tab. 3 for the exact values). The evidence for an additional break, , at intermediate fluxes varies depending on the energy bin. From 1.89–4.75 GeV, there is strong indication for a break at fluxes ph cm s, with the index above the break hardening to below the break. In the two subsequent energy bins, up to GeV, we also find evidence that the source-count distribution hardens as we move from high fluxes to below the second break, with the index below the second break 1.9-2.0 in both cases. In the last bin, the uncertainties are too large to determine if the source-count distribution changes slope at any flux above the lowest break .

#### Top Three PSF Quartiles

The benchmark analysis described in the previous section used only the top quartile (PSF3) of the Pass 8 ultracleanveto data set. This restriction selects events with the best angular resolution, but at the price of reducing the total photon count. In Sec. 3, we showed that including the top three quartiles of the Pass 8 ultracleanveto data may help constrain the source-count distribution at low fluxes. With that in mind, we now investigate how the results of the benchmark analysis change when using the PSF1–3 ultracleanveto data set.

In general, the best-fit intensities for the individual spectral components are consistent within uncertainties with those obtained using only the top quartile of data. (See Tab. S1 for specific values.) The PS flux does increase slightly in going from PSF3 to PSF1–3 in the upper energy bins due to the increased exposure. More specifically, the ratios of the median PS intensities measured with ultracleanveto PSF1–3 data to those measured with PSF3 data are in the four increasing energy bins. This can also be seen in the associated spectral intensity plot (right panel of Fig. 3), where the red bands are further above the 3FGL line in the last energy bins than in the corresponding plot for the PSF3 analysis (left panel). The intensity of the EGB is seen to increase slightly, in all energy bins, when going from PSF3 to PSF1–3 data, potentially suggesting additional cosmic-ray contamination with the looser photon-quality cuts, though the increases in EGB intensities are within statistical uncertainties.

The best-fit source-count distributions recovered by the NPTF with PSF1–3 data are shown in Fig. 5. For reference, the blue curve shows the best-fit for the PSF3–only analysis. The most important difference between the PSF3 and PSF1–3 results is that the source-count distributions extend to lower flux with PSF1–3 data. This is due to the fact that the exposure in each energy bin, averaged over the region of interest, is larger for the top three quartiles compared to the top quartile alone. As a result, the flux corresponding to single-photon detection is lower (compare the vertical dot-dashed line in Fig. 5 with that in Fig. 4), which improves the NPTF reach. Thus, the PSF1–3 analysis is sensitive to more sub-threshold sources. Note that the same trend was observed in the simulation tests in Sec. 3 in going from PSF3 to PSF1–3 data sets.

Other than the location of the lowest break, which is lower due to the increased exposure, all other source-count distribution parameters are consistent, within uncertainties, between analyses. At the lowest energy, the break at photons cm s is even more pronounced, with an index above the break and below the break. In the highest energy bin, the structure observed in the source-count distribution for the benchmark analysis has smoothed out.

### 4.2 Systematic Tests

The previous subsection illustrated how the results of the NPTF change when additional ultracleanveto PSF quartiles are included in the analysis. We also tested the stability of our analysis to variations in the region of interest, *Fermi* event class, foreground modeling, *Fermi* bubbles, PSF modeling, and choice of priors. These systematic tests are described in detail in Appendix B.2.

Figure LABEL:fig:systematicsplot briefly summarizes the results. The EGB intensity as measured by *Fermi* is shown by the gray band. To obtain this band, we use the best-fit power-law spectrum with exponential cut-off provided in [15]; the width of the gray band is found by varying between best-fit values for the three foreground models considered in that paper (Models A/B/C) and does not include statistical uncertainties, which become increasingly important at high energies. The smooth isotropic intensity, and thus the intensity of the EGB, is subject to large systematic uncertainties. As expected, the variation in smooth isotropic intensity is most pronounced when using the source event class, which contains more cosmic-ray contamination. However, the spectrum of emission from PSs as captured by the NPTF appears robust to all the systematic effects considered here. This is the primary conclusion of this subsection.

## 5 High-Energy Analysis: 50–2000 GeV

We now consider the NPTF results at high energies from 50–2000 GeV. The number of photons available decreases when moving to higher energies, so we loosen the restrictions on the PSF quartiles to maximize the sensitivity potential of the NPTF. In this section, the majority of the analyses are done using all quartiles of the ultracleanveto data, though we also show results using all quartiles of source data. For the same reason, we widen the ROI to rather than , although the results are not sensitive to this cut, as we will show.

The best-fit energy spectra recovered by the NPTF analysis for the high-energy study of ultracleanveto data is shown in the bottom right panel of Fig. 6.^{12}*Fermi*’s 2FHL catalog [16] (dashed black line). This recently-published catalog is based on 80 months of data and focuses on hard sources in the range from 50–2000 GeV. Statistical and systematic uncertainties are not accounted for in the determination of the 2FHL spectrum in Fig. 6; these are likely non-negligible, especially at the highest energies.

The best-fit source-count distributions for the three energy bins are also shown in Fig. 6, in the top row and bottom left panel. The black points in those panels denote the 2FHL source-count distributions, with vertical error bars indicating 68% Poisson errors. The statistical errors on the 2FHL sources are large due to the fact that there are not many sources. In all energy bins, the NPTF places the lowest break close to the single-photon sensitivity threshold (vertical dot-dashed line) and the highest break in the vicinity of the brightest 2FHL source, just as in the low-energy analysis. Most notably in the 50–151 GeV bin, the NPTF probes unresolved sources with fluxes nearly an order-of-magnitude below the apparent 2FHL threshold. We find no evidence for an additional intermediate-flux break in any of the energy bins, although it is difficult to make conclusive statements due to the large uncertainties in the individual source-count distributions.

We have completed a number of systematic tests of the high-energy analyses that include looking at all quartiles of the source data, requiring for both event classes, and using the third alternate prior choice, with . The results are summarized in Fig. 7, and the source-count distributions for each of the systematic tests are given in Appendix C. Importantly, the isotropic-PS intensity is consistent across all the tests. However, the EGB intensities recovered by the NPTF are, in general, higher than those measured by *Fermi*. This discrepancy is likely due to increased cosmic-ray contamination above 100 GeV, as suggested by the high IGRB intensities recovered by the NPTF at these energies. Indeed, the *Fermi* EGB study on Pass 7 data [15] used dedicated event classes with specific data cuts to minimize such contributions. Such an analysis is beyond the scope of this work, as our primary focus is on the PS populations. We simply caution the reader that the derived intensity for the smooth isotropic component in the high-energy analyses is subject to potentially large contamination.

It is possible to make stronger statements about the best-fit source-count distribution at high energies if we consider the wide-energy bin from 50–2000 GeV. The results are shown in the left panel of Fig. 8. Due to the improved statistics, the uncertainties on the source-count distribution are smaller than those for the three sub-bins. Other than the low-flux sensitivity break, the NPTF finds no preference for an additional break. The intermediate-flux break, , is essentially unconstrained as a result, and the power-law slope above (below) it are consistent within uncertainties: and , respectively. We compare this result to the best-fit source-count distribution (blue line) published by *Fermi* for sources in this same energy range [17]. There are important differences between the two analyses. In the *Fermi* study, simulated maps were created using several different source-count distributions, parametrized as singly broken power laws. The histogram of the photon-count distribution for each of these maps, averaged over the full region of interest, was compared to the actual data, and a fit was done to select the simulated maps that most closely resembled the data.
This method is related to but in many ways distinct from the NPTF. The NPTF considers the difference between Poissonian and non-Poissonian photon probability distributions at the pixel-by-pixel level, instead of averaging the distributions over the full region. Moreover, in our analysis we rely on semi-analytic techniques to calculate the photon-count probability distributions as we scan over the space of model parameters, instead of relying on Monte Carlo samples to numerically construct these distributions. As a result, we are able to consider source-count distributions with additional degrees of freedom and also scan over the normalizations of all of the background templates, which tend to be well determined given the pixel-by-pixel nature of the fit. In contrast, the intensity of all Poissonian models in [17], including the smooth isotropic emission, was kept fixed while scanning over the source-count distribution degrees of freedom.

The cumulative source-count plot is provided in the right panel of Fig. 8. Our result is in good agreement with the 2FHL sources above the catalog sensitivity threshold ph cm s. In the first few flux bins above this threshold, there appear to be more 2FHL sources than what is predicted by the NPTF, although the results are still consistent within uncertainties. This may be due to the Eddington bias [46] where extra sources are observed above threshold due to upward statistical fluctuations from sources immediately below.

Based on the results in Fig. 8, we can project the expected number of these sources that may be observed by the Cherenkov Telescope Array (CTA) [38, 45]. For energies above 50 GeV, the CTA flux sensitivity is cm s for 50 hours of observation per field-of-view (5 detection).^{13}*Fermi* study using the blue source-count distribution illustrated in Fig. 8.

## 6 Discussion and Conclusions

The primary focus of this paper is to characterize the properties of the PSs contributing to the EGB in a data-driven manner. To achieve this, we use a novel analysis method, referred to as Non-Poissonian Template Fitting (NPTF), which takes advantage of photon-count statistics to distinguish diffuse and PS contributions to gamma-ray maps with non-trivial spatial variations. We presented the NPTF results on *Fermi* Pass 8 data at low (1.89–94.9 GeV, ) and high (50–2000 GeV, ) energies. For the first time, the intensity and source-count distributions for the isotropic PSs have been obtained as a function of energy, up to 2 TeV. The best-fit source-count distributions probe fluxes below the current detection threshold for the *Fermi* 3FGL and 2FHL catalogs, providing information on the unresolved populations.

Through extensive studies of how the NPTF responds to simulated populations, we have shown that the analysis procedure reproduces the properties of input source classes. Therefore, the features of the best-fit source-count distributions obtained from the data provide a potential wealth of information about the source populations of the EGB. While a detailed interpretation of the source-count distributions in terms of particular theoretical models is beyond the scope of this paper, several important trends were observed.

In this work, the source-count distributions are parametrized as triply-broken power laws in the NPTF. At all energies, a break is fit at low (high) fluxes, below (above) which the analysis method loses sensitivity. Of particular interest is whether an additional break, , is preferred at intermediate flux. We find a break in the lowest energy bin (1.89–4.75 GeV) at cm s with slope above and below. In the subsequent two energy bins, 4.75–11.9 GeV and 11.9–30.0 GeV, there is a mild indication that the source-count distribution hardens below the intermediate flux break, though the change in slope is not as robust and significant as in the lowest energy bin. At higher energies, above 30 GeV, there is no indication that the source-count distribution changes slope at the intermediate break. This trend is in line with the expectations from the blazar simulations in Sec. 3. For example, in both Figs. 1 and LABEL:fig:bl2dnds (see also Fig. S2), which show the results of the NPTF run on simulated data with the Blazar–1 and Blazar–2 models, we find evidence for curvature in the source-count distribution at intermediate fluxes in the lowest energy bins, while at higher energies the recovered source-count distribution appears as a single power law at fluxes above the sensitivity threshold of the NPTF. In the energy bin from 50–2000 GeV the best-fit value for is essentially unconstrained and the slopes above and below it are consistent within uncertainties: and .

The NPTF also provides the best-fit intensities for the isotropic-PS populations as a function of energy. Figure 9 illustrates this spectrum for analyses done using the ultracleanveto event class. The filled red circles (open red boxes) show the results for the dedicated low (high)-energy analysis, with PSF1–3 data used at low energies and PSF0–3 data at high energies. For comparison, the *Fermi* EGB spectrum is shown by the black line [15]. This corresponds to the best-fit intensity using the Model A diffuse background from that study. To illustrate the systematic uncertainty on this curve, we also plot the spectra for diffuse models B and C (dashed and dotted, respectively).

Low-Energy Analysis | High-Energy Analysis | |||||||
---|---|---|---|---|---|---|---|---|

1.89–4.75 | 4.75–11.9 | 11.9–30 | 30–94.9 | 50–151 | 151–457 | 457–2000 | 50–2000 | |

Scenario A | ||||||||

Scenario B |

*Fermi*EGB intensity from [15], with foreground Model A (‘Scenario B’). Although the

*Fermi*analysis uses a different foreground model, it takes advantage of a dedicated event selection above 100 GeV that mitigates effects of additional contamination.

The PS fraction, defined as , is provided in Tab. 4 for each energy bin. While using the EGB intensity derived in this work (‘Scenario A’) is the most self-consistent comparison, this may underestimate the PS contribution above 100 GeV, where the NPTF appears to recover too much smooth isotropic emission due to increased cosmic-ray contamination in the data sets used, as already discussed. Therefore, we also show the PS fractions calculated relative to the *Fermi* EGB intensity from [15] for diffuse model A (‘Scenario B’). The comparison to the EGB as measured in [15] is not fully self consistent, since, for example, the foreground modeling and data sets in [15] differ from those used in this work to measure . However, the advantage of this comparison is that the *Fermi* analysis uses special event-quality cuts to mitigate contamination, and thus their measure of is likely more faithful than that presented in this work. These results are shown in the second row of Tab. 4. For the low-energy analysis, the PS fractions are consistent, within uncertainties, when is taken from our work or *Fermi*’s.^{14}*Fermi*’s EGB intensity.

In general, we find that approximately 50–70% of the EGB consists of PSs in the energy ranges considered. To interpret these results, we use the ratios obtained in the simulation studies of Sec. 3. In that section, we showed that the efficiency for the NPTF to recover the flux for the Blazar–2 model (with PSF1–3) is 100% in the first energy bin and drops to 60% in the fourth energy bin. For the Blazar–1 model, the efficiencies are consistently higher than the Blazar–2 scenario. These two blazar models are meant to illustrate extreme scenarios, with the Blazar–1 model having a significant fraction of the total flux arising from high-flux sources, while low-flux sources dominate instead in the Blazar–2 case. The high efficiency of the NPTF to recover the blazar component at low energies, combined with the PS fractions observed in the data (Tab. 4), clearly suggests that there is a substantial non-blazar component of the EGB up to energies 30 GeV. The interpretation of the results in the energy bin from 30.0–94.9 GeV is less clear. A proper interpretation of the results at higher energies in terms of evidence for or against a non-blazar component of the EGB requires dedicated blazar simulations, which we leave to future work.

Our results tend to predict fewer PSs (and photons from PSs) where we do overlap with previous studies. For example, a similar photon-count analysis was used by [104] to study 1–10 GeV energies in the Pass 7 Reprocessed data. They found an 80% PS fraction at these energies. At the lowest energies that we probe—which admittedly do not extend down as low as 1 GeV—we only find a 54% PS fraction (relative to Model A). Systematic uncertainties, as shown in Fig. LABEL:fig:systematicsplot, can affect the recovered PS intensities at the level, which can partially alleviate the tension between our results.

Above 50 GeV, the NPTF procedure predicts that of the EGB consists of PSs, with systematic uncertainties estimated at approximately . This fraction is smaller, and in slight tension, with the predicted value obtained in previous work [17]. The fact that our results suggest that there is more diffuse isotropic emission at high energies may help alleviate the tension between [17] and the hadronuclear () interpretation of IceCube’s PeV neutrinos [75]. Some models suggest, for example, that these very-high-energy neutrinos are produced in hadronuclear interactions, along with high-energy gamma-rays that would contribute to the IGRB [75, 96, 24, 62]. If the smooth isotropic gamma-ray spectrum (*i.e.*, the non-blazar spectrum) is suppressed above 50 GeV in the *Fermi* data, it could put such scenarios in tension with the data [26, 77]; however, that does not necessarily appear to be the case given the results of our analysis [79]. With that said, and as already mentioned, dedicated blazar simulations at high energies are needed to properly interpret our results at these energies.

The PS spectrum in Fig. 9 is well-modeled () as a power law with an exponential cut-off:

(5) |

where GeVcmssr, , and GeV are the best-fit parameters.^{15}

The NPTF allows us to make statistical statements about the properties of source populations contributing to the EGB, but at the expense of identifying the precise locations of these sources. However, it is still possible to make probabilistic statements about these locations. To do so, we compare the observed photon count in a given pixel, , to the mean expected value, , without accounting for PSs. To determine we include the diffuse background, smooth isotropic emission, and the Fermi bubbles templates, with normalizations as determined from the NPTF. The pixel-dependent survival function is defined as

(6) |

where CDF is the Poisson cumulative distribution function. The smaller the value of (or, conversely, the larger the value of ), the more probable it is that the pixel contains a PS. Figure 10 shows full-sky maps of for both low (1.89–94.9 GeV) and high (50–2000 GeV) energies.^{16}

Deciphering the constituents of the EGB remains an important goal in the study of high-energy gamma-ray astrophysics, with broad implications extending from the production of PeV neutrinos to signals of dark-matter annihilation or decay. The *Fermi* LAT has already played an important role in the discovery of many new sources in the GeV sky. By taking advantage of the statistical properties of unresolved populations, our results provide a glimpse at the aggregate properties of the sources that lie below the detection threshold of these published catalogs and suggest a wealth of detections for future observatories.

## Note Added

As this work was being completed, we became aware of a complementary analysis that also takes advantage of photon-count statistics to obtain the source-count distributions in the energy range from 1–171 GeV [103]. That study used the clean event class of the Pass 7 Reprocessed data to study the source-count distributions in five energy bins from 1–171 GeV. A direct comparison between the two results is challenging given the different data sets, foreground models, priors, energy sub-bins, pixelation, and regions of interest that were used between the two studies. In the energy range where the two studies overlap, we tend to find smaller PS fractions. (For a discussion of why we do not consider energies below 1.89 GeV, see Appendix A.2.) Specifically, [103] finds PS fractions of , , , and in the energy ranges [1.99, 5.0], [5.0, 10.4], [10.4, 50.0], and [50.0, 171] GeV. (PS fractions are quoted relative to the *Fermi* EGB intensity for foreground Model A [15].) These numbers can be compared to our results, summarized in Tab. 4. The best-fit source-count distributions recovered by [103] closely resemble the Blazar–2 scenario that we studied in simulations (see Sec. 3). Our studies on simulated data show that the NPTF can successfully recover the energy spectrum and source-count distributions for this blazar population, both for the case where it singularly contributes to the EGB, as well as the case where it contributes along with SFGs, modeled according to [96]. However, we observe different features in the best-fit source-count distributions when running the NPTF on the actual data.

## Acknowledgements

We thank the *Fermi* Collaboration for the use of the *Fermi* public data and the Fermi Science Tools. We also thank J. Balkind, K. Bechtol, K. Blum, M. Di Mauro, N. Rodd, and T. Slatyer for helpful conversations. LN is supported by the U.S. Department of Energy under grant Contract
Number DE-SC00012567. ML is supported by the DoE under grant Contract Number DE-SC0007968. BRS is supported by a Pappalardo Fellowship
in Physics at MIT.

## Appendix A Simulations of Extragalactic Gamma-ray Point Sources

In this Appendix, we provide further details about the simulations and analyses of extragalactic point sources.

### a.1 Energy-Binned Source-Count Distributions

We generate simulated maps directly from the source-count distribution . To obtain this, we need two inputs: the gamma-ray luminosity function, , and the source energy spectrum, [41]. Typically, the luminosity function (LF) is given by

(S1) |

where is the comoving volume, is the photon spectral index, is the redshift, is the number of sources, and is the rest-frame luminosity for energies from 0.1–100 GeV in units of GeV s.

The photon flux in this energy range, , is defined in terms of the source energy spectrum,

(S2) |

where the units are cm s, and GeV.

The source-count distribution is then given by

(S3) |

which can be accurately estimated as

(S4) |

where is the full-sky solid angle, is the comoving volume slice for a given redshift and is sufficiently small. To calculate , we need the following expression, which relates the luminosity to the energy flux:

(S5) |

where is the luminosity distance. For a given and , one can use (S2) to solve for the normalization of , which can be substituted into (S5), along with and , to obtain the associated value of the luminosity. The photon flux, , is related to the photon count, , via the mean exposure , which is averaged over 0.1–100 GeV and the ROI. This allows us to finally obtain from (S4).

The procedure outlined above allows one to obtain the source-count distributions based on models of luminosity functions and spectral energy distributions provided in the literature. For the AGN and SFG examples we consider in detail in this work, the luminosity functions correspond to photon energies from 0.1–100 GeV. However, we also need the source-count distributions in subset energy ranges corresponding to our energy bins of interest, with GeV. We rescale the fluxes for these individual energy bins of interest to those in the provided 0.1–100 GeV range using a procedure similar to [41]. Denoting quantities associated with this energy bin with a prime, we can write the new source-count distribution as

(S6) |

where is again sufficiently small—we set , and verify that the answer is robust to this choice. Note that the integral must still be done over (unprimed) because the luminosity function is explicitly defined in terms of it. So, we must solve for the photon flux over the full energy, , in terms of the value in the sub-bin, . The two are related via a proportionality relation

(S7) |

where the exponential factor accounts for the attenuation due to extragalactic background light (EBL) [60, 47, 89, 55, 13, 8, 44]. It arises from pair annihilation of high-energy gamma-ray photons with other background photons in infrared, optical, and/or ultraviolet, and is described by the optical depth, . We use the EBL attenuation model from [52].

Additionally, the expected gamma-ray spectrum can be calculated from the luminosity function as

(S8) |

We use this equation to appropriately weight the number of photons per energy sub-bin for the individual sources when creating simulated maps. This ensures that the variations in PSF and exposure within the larger energy bins used in the NPTF analyses are properly accounted for in the simulation procedure.

### a.2 Simulations at Energies 1.89 GeV

The main analyses presented in the text do not consider energies below 1.89 GeV because we find that systematic effects may become more important at these low energies. In particular, simulations done with the set of priors presented in Sec. 2 show that the NPTF can over-estimate the PS intensity at very low energies. That is, in simulations with both blazars and SFGs, the isotropic-PS template tends to absorb more emission than is simulated in blazars, while the smooth isotropic template absorbs less emission than is simulated in SFGs.

As an illustration, we show the results when the NPTF is run on a simulated map of Blazar–1 and SFG sources. Figure S1 shows the best-fit source-count distributions for the energy ranges 0.475–0.753 and 0.753–1.89 GeV. The PS fractions in these bins are , illustrating that the PS template is absorbing more PS emission than is simulated in blazars. From Fig. S1, we see that this is likely the result of the fit typically predicting more sources than it should at intermediate and low fluxes. This could be due to a variety of factors. For example, at very low energies the angular resolution of the detector quickly degrades, and the PS flux also becomes more and more sub-dominant compared to foreground emission. For these reasons—and out of an abundance of caution—we do not present results of the NPTF on data below 1.89 GeV. It is certainly possible that analyses at low energies could provide a wealth of interesting observations. However, having confidence in the NPTF results at such low-energies requires a more careful consideration of systematics, which goes beyond the scope of the present work.

### a.3 Ultracleanveto PSF1-3 Simulation Analysis

The simulated-data studies in the main text used the PSF3 instrument response function for the Pass 8 ultracleanveto data set. Here, we show what happens when using the PSF1–3 instrument response function instead. Including the top three quartiles of data increases the total exposure, though at the cost of decreased angular resolution. Figure S2 illustrates the result for the Blazar–2 simulated data set. The best-fit source-count distribution extends to lower fluxes than the corresponding distribution for top-quartile data in Fig. LABEL:fig:bl2dnds.

### a.4 Simulations of SFGs

Here, we show the results of running the NPTF analysis on simulated data where the EGB arises solely from SFGs. The resulting best-fit source-count distributions are shown in Fig. S3. In the first energy bin, the brightest SFGs, which contribute little more than 1 photon, are detected as PSs by the NPTF. At higher energies, the best-fit source-count distributions are consistent with zero and no significant evidence for a PS population is found. The energy spectrum plot in Fig. S4 shows that the SFG flux is absorbed by the smooth isotropic template. In comparison, the intensity absorbed by the isotropic-PS template is completely subdominant and is several orders of magnitude lower than its Poissonian counterpart at all energies.

## Appendix B Supplementary Results: Low-Energy Analysis

### b.1 Best-Fit Intensities and Posterior Distributions

This section includes supplementary information pertaining to the low-energy analysis presented in Sec. 4.1. In particular, Tabs. S1 and S2 present the best-fit intensities and source-count parameters for the NPTF analysis of the top-three quartiles of ultracleanveto data. Figures S5–LABEL:fig:p8triangle4 show the posterior distributions, in each energy bin, for the benchmark analysis.

Energy | |||||
---|---|---|---|---|---|

1.89–4.75 | |||||

4.75–11.9 | |||||

11.9–30.0 | |||||

30.0-94.9 |

Energy | |||||||
---|---|---|---|---|---|---|---|

1.89–4.75 | |||||||

4.75–11.9 | |||||||

11.9–30.0 | |||||||

30.0-94.9 |