{'fraction_non_alphanumeric': 0.03860460229063781, 'fraction_numerical': 0.02356835137123043, 'mean_word_length': 4.3233583174851775, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Unveiling the Gamma-ray Source Count Distribution Below the Fermi Detection Limit with Photon Statistics

The source-count distribution as a function of their flux, dN/dS, is one of the main quantities characterizing gamma-ray source populations. We employ statistical properties of the Fermi Large Area Telescope (LAT) photon counts map to measure the composition of the extragalactic gamma-ray sky at high latitudes (|b|>30 deg) between 1 GeV and 10 GeV. We present a new method, generalizing the use of standard pixel-count statistics, to decompose the total observed gamma-ray emission into (a) point-source contributions, (b) the Galactic foreground contribution, and (c) a truly diffuse isotropic background contribution. Using the 6-year Fermi-LAT data set (P7REP), we show that the dN/dS distribution in the regime of so far undetected point sources can be consistently described with a power law of index between 1.9 and 2.0. We measure dN/dS down to an integral flux of ~2x10^{-11} cm^{-2}s^{-1}, improving beyond the 3FGL catalog detection limit by about one order of magnitude. The overall dN/dS distribution is consistent with a broken power law, with a break at 2.1^{+1.0}_{-1.3}x10^{-8} cm^{-2}s^{-1}. The power-law index n_1=3.1^{+0.7}_{-0.5} for bright sources above the break hardens to n_2=1.97+-0.03 for fainter sources below the break. A possible second break of the dN/dS distribution is constrained to be at fluxes below 6.4x10^{-11} cm^{-2}s^{-1} at 95% confidence level. The high-latitude gamma-ray sky between 1 GeV and 10 GeV is shown to be composed of ~25% point sources, ~69.3% diffuse Galactic foreground emission, and ~6% isotropic diffuse background.


INTRODUCTION
The decomposition of the extragalactic gamma-ray background (EGB; see Fornasa & Sánchez-Conde (2015) for a recent review) is pivotal for unveiling the origin of the nonthermal cosmic radiation field. The EGB comprises the emission from all individual and diffuse gamma-ray sources of extragalactic origin, and thus it originates from different mechanisms of gamma-ray production in the Universe. The EGB can be dissected by resolving the various point-source contributions, characterized by their differential source-count distribution dN/dS as a function of the integral source flux S (see, e.g., Abdo et al. 2010b;Singal 2015). Conventionally, the EGB emission that is left after subtracting the resolved gamma-ray sources is referred to as the isotropic diffuse gamma-ray background (IGRB; Ackermann et al. 2015b). The Large Area Telescope (LAT) on board the Fermi satellite (Ackermann et al. 2012) has allowed the discovery of more than 3,000 gamma-ray point sources, collected in the 3FGL catalog (Acero et al. 2015). Resolved sources amount to about 30% of the EGB (Ackermann et al. 2015b) below ∼100 GeV (while above ∼100 GeV this percentage can rise to about 50%).
For resolved point sources listed in catalogs the dN/dS distributions of different source classes can be charac-terized. Among these, blazars represent the brightest and most numerous population, and, consequently, their dN/dS is the best-determined one. Blazars exhibit two different subclasses: flat-spectrum radio quasars (FSRQs), with a typically soft gamma-ray spectrum characterized by an average power-law photon index of ∼2.4, and BL Lacertae (BL Lac) objects, with a harder photon index of ∼2.1. The dN/dS distribution of blazars has been studied in detail in several works (Ajello et al. , 2014Broderick et al. 2014a;Di Mauro et al. 2014c;Harding & Abazajian 2012;Inoue & Totani 2009;Stecker & Venters 2011;Stecker & Salamon 1996). Besides blazars, the EGB includes fainter sources like misaligned active galactic nuclei (mAGN; Di Mauro et al. 2014a; Inoue 2011) and star-forming galaxies (SFGs; Ackermann et al. 2012b;Fields et al. 2010;Lacki et al. 2014;Tamborra et al. 2014; Thompson et al. 2006). A contribution from Galactic sources such as millisecond pulsars (MSPs) located at high Galactic latitude is possible, although it has been constrained to be subdominant (Calore et al. 2014;Gregoire & Knodlseder 2013). Finally, pure diffuse (not point-like) components can contribute, for instance caused by pair halo emission from AGN, clusters of galaxies, or cascades of ultra-high-energy cosmic rays (UHECRs) on the CMB (see Fornasa & Sánchez-Conde (2015) and references therein).
In the usual approach, the dN/dS distributions of different populations (inferred from resolved sources) are extrapolated to the unresolved regime and used to investigate the composition of the IGRB (i.e., the unresolved EGB). This approach has revealed that the abovementioned three main components well explain the observed IGRB spectrum, constraining further contribu-tions to be subdominant, including a possible exotic contribution from dark matter (DM) annihilation or decay Cholis et al. 2014;Di Mauro & Donato 2015). While the above-mentioned approach is very useful, a clear drawback is caused by the fact that it relies on the extrapolation of dN/dS distributions. In this work, we will focus on a method to overcome this problem by conducting a direct measurement of the dN/dS in the unresolved regime.
Detection capabilities for individual point sources are intrinsically limited by detector angular resolution and backgrounds. This makes in particular the IGRB a quantity that depends on the actual observation (Ackermann et al. 2015b). The common approach of detecting individual sources (Acero et al. 2015;Ackermann et al. 2016) can be complemented by decomposing gamma-ray skymaps by statistical means, using photon-count or intensity maps. One of the simplest ways of defining such a statistic is to consider the probablity distribution function (PDF) of photon counts or fluxes in pixels, commonly known as P (D) distribution in the radio (e.g., Condon 1974;Scheuer 1957;Vernstrom et al. 2014Vernstrom et al. , 2015, and references therein) and X-ray (e.g., Hasinger et al. 1993; So ltan 2011, and references therein) bands. Recently, this technique has been adapted to photon-count measurements in the gamma-ray band; see Malyshev & Hogg (2011), henceforth MH11, for details. Various theoretical studies have also been performed (Baxter et al. 2010;Dodelson et al. 2009;Feyereisen et al. 2015;Lee et al. 2009). In addition, this method has been used to probe unresolved gamma-ray sources in the region of the Galactic Center (Bartels et al. 2016;Lee et al. 2015Lee et al. , 2016, as well as to constrain the source-count distribution above 50 GeV (Ackermann et al. 2015a).
As argued above, this method has the advantage of directly measuring the dN/dS in the unresolved regime, thus not relying on any extrapolation. A difference with respect to the use of resolved sources is that in the PDF approach only the global dN/dS, i.e., the sum of all components, can be directly measured: since no individual source can be identified with this method, counterpart association and the separation of dN/dS into different source components become impossible. The PDF approach nonetheless offers another important advantage with respect to the standard method: the use of the dN/dS built from cataloged sources close to the detection threshold of the catalog is hampered by the fact that the threshold is not sharp but rather characterized by a detection efficiency as a function of flux (Abdo et al. 2010b;Ackermann et al. 2015a). The dN/dS thus needs to be corrected for the catalog detection efficiency, which, in turn, is a nontrivial quantity to determine (Abdo et al. 2010b). On the contrary, the PDF approach treats all the sources in the same way, resolved and unresolved, and can thus determine the dN/dS in a significantly larger flux range, without requiring the use of any efficiency function.
In the following, we will measure the high-latitude dN/dS with the PDF methodology using 6 years of gamma-ray data collected with the Fermi-LAT. We will show that for the 1 GeV to 10 GeV energy band we can measure the dN/dS down to an integral flux of ∼ 10 −11 cm −2 s −1 , which is a factor of ∼ 20 lower than the nominal threshold of the 3FGL catalog.
This article is structured as follows: In Section 2 we introduce the mathematical framework of the analysis method, supplemented by a detailed description of our extensions to previous approaches, the modeling of source and background components, and the fitting procedure. The gamma-ray data analysis is addressed in Section 3. Section 4 is dedicated to details of the statistical analysis approach and the fitting technique. The resulting global source-count distribution and the composition of the gamma-ray sky are considered in Section 5. Section 6 addresses the angular power of unresolved sources detected with this analysis. Possible systematic and modeling uncertainties are discussed in Section 7. Eventually, final results are summarized in Section 8.

THE STATISTICS OF GAMMA-RAY PHOTON COUNTS
In the present analysis, we assume the gamma-ray sky at high Galactic latitudes to be composed of three different contributions: • A population of gamma-ray point sources. Given that the analysis is restricted to high Galactic latitudes, this source population is considered to be dominantly of extragalactic origin. Sources can thus be assumed to be distributed homogeneously across the sky.
• Diffuse gamma-ray emission from our Galaxy, mostly bright along the Galactic plane but extending also to the highest Galactic latitudes. We will refer to this component as Galactic foreground emission. The photon flux in map pixel p from this component will be denoted as F gal .
• Gamma-ray emission from all contributions that are indistinguishable from diffuse isotropic emission, such as extremely faint sources. We will include in this component possible truly diffuse emission of extragalactic or Galactic origin, such as, for example, gamma rays from cosmological cascades from UHECRs, or possible isotropic subcomponents of the Galactic foreground emission. In addition, the component comprises the residual cosmic-ray background. All together this emission will be denoted as F iso .
A more detailed account of the individual components is given in Section 1 and later in this section. Following the method of MH11, we considered the celestial region of interest (ROI) to be partitioned into N pix pixels of equal area Ω pix = 4πf ROI /N pix sr, where f ROI is the fraction of sky covered by the ROI. The probability p k of finding k photons in a given pixel is by definition the 1-point PDF (1pPDF). In the simplest scenario of purely isotropic emission, p k follows a Poisson distribution with an expectation value equal to the mean photon rate. The imprints of more complex diffuse components and a distribution of point sources alter the shape of the 1pPDF, in turn allowing us to investigate these components by measuring the 1pPDF of the data.
The usual way in which the 1pPDF is used requires us to bin the photon counts of each pixel into a histogram of the number of pixels, n k , containing k photon counts, and to compare the p k predicted by the model with the estimator n k /N pix . This method is the one adopted by MH11. By definition, this technique does not preserve any spatial information of the measurement or its components (for example, the uneven morphology of the Galactic foreground emission), resulting in an undesired loss of information. We will instead use the 1pPDF in a more general form, including pixel-dependent variations in order to fully exploit all the available information.

Generating Functions
An elegant way of deriving the 1pPDF including all the desired components exploits the framework of probability generating functions (see MH11 and references therein for details). The generating function P (p) (t) of a discrete probability distribution p (p) k , which may depend on the pixel p and where k = 0, 1, 2, . . . is a discrete random variable, is defined as a power series in an auxiliary variable t by (1) The series coefficients p (p) k can be derived from a given P (p) (t) by differentiating with respect to t and evaluating them at t = 0, The method of combining individual components into a single P (p) (t) makes use of the summation property of generating functions, i.e., the fact that the generating function for the sum of two independent random variables is given by the product of the generating functions for each random variable itself. In our case, the general representation of P (p) (t) for photon-count maps can be derived from considering a superposition of Poisson processes; see Appendix A and MH11 for a more detailed explanation. The generating function is therefore given by where the coefficients x (p) m are the expected number of point sources per pixel p that contribute exactly m photons to the total photon count of the pixel, and m is a positive integer. In the derivation of Equation (3), it has been assumed that the x m are related to the differential source-count distribution dN/dS, where S denotes the integral photon flux of a source in a given energy range [E min , E max ], by The number of counts C (p) (S) expected in pixel p is given as a function of S by for sources with a power-law-type energy spectrum ∝ E −Γ , where Γ denotes the photon index and the pixeldependent exposure 5 as a function of energy is denoted by E (p) (E). In Equation (4), we have assumed that the PDF for a source to contribute m photons to a pixel p follows a Poisson distribution with mean C (p) (S). Gammaray sources have been assumed to be isotropically distributed across the sky, i.e., dN/dS is pixel independent, while, in principle, Equation (4) allows for an extension of the method to spatially dependent dN/dS distributions. The generating functions for diffuse background components correspond to 1-photon source terms, with x (p) m = 0 for all m except m = 1: where x (p) diff denotes the number of diffuse photon counts expected in pixel p for a given observation. 6 This quantity is given by with f (p) diff (E) being the differential flux of the diffuse component as a function of energy.
The relation in Equation (4) allows measuring the source-count distribution dN/dS from pixel-count statistics. Furthermore, we can observe that the 1pPDF approach may allow the detection of point-source populations below catalog detection thresholds: if the sourcecount distribution implies a large number of faint emitters, pixels containing photon counts originating from these sources will be stacked in an n k -histogram, increasing the statistical significance of corresponding k-bins. The average number of photons required from individual sources for the statistical detection of the entire population will therefore be significantly smaller than the photon contribution required for individual source detection.
The simple 1pPDF approach refers to a measurement of p k which is averaged over the considered ROI. The generating function for the 1pPDF measurement therefore reduces to a pixel average, where we made use of the fact that the total generating function factorizes in the point-source component and the diffuse component, P (p) (t) = P  (3) and (6)). 5 The experiment exposure, which depends on energy and position, is discussed in Section 3.
6 Equation (6) can be derived from Equation (1)  The numerical implementation of Equation (8) in its most general form is computationally complex (see MH11; Lee et al. 2015). In the ideal situation of an isotropic point-source distribution and homogeneous exposure, P S (t) ≡ P (p) S (t) factorizes out of the sum, reducing the pixel-dependent part of Equation (8) to the diffuse component, which is easy to handle. The exposure of Fermi -LAT data is, however, not uniformly distributed in the ROI (see Section 3) and requires appropriate consideration.
To correct the point-source component for exposure inhomogeneities, we divided the exposure map into N exp regions, separated by contours of constant exposure such that the entire exposure range is subdivided into N exp equally spaced bins. In each region, the exposure values were replaced with the region averages, yielding N exp regions of constant exposure. The approximation accuracy is thus related to the choice of N exp . In this case, Equation (8) reads where P i = {p|p ∈ R i } denotes the subset of pixels belonging to region R i . In this way, P S (t) becomes independent of the inner sum and factorizes, significantly reducing the required amount of computation time.
The probability distributions p k or p (p) k can eventually be calculated from P(t) or P (p) (t), respectively, by using Equation (2).

Source-count Distribution
The source-count distribution dN/dS characterizes the number of point sources N in the flux interval (S, S+dS), where S is the integral flux of a source in a given energy range. The quantity N actually denotes the areal source density per solid angle element dΩ, which is omitted in our notation for simplicity. In this analysis, we parameterized the source-count distribution with a power law with multiple breaks, referred to as multiply broken power law (MBPL) in the remainder. An MBPL with N b breaks located at S bj , j = 1, 2, . . . , N b , is defined as where S 0 is a normalization constant. The n j denote the indices of the power-law components. The dN/dS distribution is normalized with an overall factor A S , which is given by A S = dN/dS (S 0 ) if S 0 > S b1 . We required a finite total flux, i.e., we imposed n 1 > 2 and n N b +1 < 2.

Source Spectra
The whole population of gamma-ray sources is disseminated by a variety of different source classes (see Section 1 for details). In particular, FSRQs and BL Lac objects contribute to the overall dN/dS at high Galactic latitudes. The spectral index distribution of all resolved sources in the energy band between 100 MeV and 100 GeV (assuming power-law spectra) is compatible with a Gaussian centered on Γ = 2.40 ± 0.02, with a half-width of σ Γ = 0.24 ± 0.02 (Abdo et al. 2010b). We thus used an index of Γ = 2.4 in Equation (5).

Galactic Foreground and Isotropic Background
The Galactic foreground and the diffuse isotropic background were implemented as described in Equation (6). The total diffuse contribution was modeled by with A gal being a normalization parameter of the Galactic foreground component x gal . For the isotropic component x (p) iso the integral flux F iso was directly used as a sampling parameter, in order to have physical units of flux.
Galactic Foreground -The Galactic foreground was modeled using a template (gll iem v05 rev1.fit) developed by the Fermi-LAT collaboration to compile the 3FGL catalog (Acero et al. 2015) 7 . The Galactic foreground model is based on a fit of multiple templates to the gamma-ray data. The templates used are radioderived gas maps splitted into various galactocentric annuli, a further dust-derived gas map, an inverse Compton emission template derived with the GALPROP code, 8 and some patches designed to describe observed residual emission not well represented by the pervious templates, such as the Fermi bubbles and Galactic Loop I.
The Galactic foreground template comprises predictions of the differential intensity at 30 logarithmically spaced energies in the interval between 50 MeV and 600 GeV. The spatial map resolution is 0.125 • , which was resampled to match the pixelization scheme and spatial resolutions used in our analysis. The predicted number of counts per pixel x (p) gal was obtained from integration in the energy range [E min , E max ] as described in Section 2.1.
In order to include the effects caused by the point spread function (PSF) of the detector, we smoothed the final template map with a Gaussian kernel of 0.5 • . We checked that systematics of this coarse PSF approximation (see Section 3) were negligible, by comparing kernels with half-widths between 0 • and 1 • . Figure 1 shows the model prediction for the diffuse Galactic foreground flux between 1 GeV and 10 GeV and Galactic latitudes |b| ≥ 30 • . The complex spatial morphology of the Galactic foreground emission is evident. The intensity of Galactic foreground emission significantly decreases with increasing latitude. The integral flux predicted by the model in the energy range ∆E between 1 GeV and 10 GeV is F gal (∆E) 4.69 × 10 −5 cm −2 s −1 for the full sky and F gal (∆E; |b| ≥ 30 • ) 6.42×10 −6 cm −2 s −1 for high Galactic latitudes |b| ≥ 30 • . gal is plotted in Galactic coordinates (l, b) using a Mollweide projection of the sphere. The Galactic Center is in the middle of the map. The Galactic plane has been masked for latitudes |b| < 30 • (in gray). The color mapping is log-linear.
Since the model reported in gll iem v05 rev1.fit was originally normalized to best reproduce the whole gamma-ray sky, we allowed for an overall different normalization parameter A gal in our analysis, given that we explored different ROIs. Nonetheless, A gal is expected to be of order unity when considered a free fit parameter.
Isotropic Background -The expected counts for the diffuse isotropic background component F iso were derived assuming a power-law spectrum with spectral index Γ iso = 2.3 (Ackermann et al. 2015b). We verified that using the specific energy spectrum template provided by the Fermi -LAT Collaboration (iso clean front v05.txt) had no impact on our results.

PSF Smearing
The detected photon flux from point sources is distributed over a certain area of the sky as caused by the finite PSF of the instrument. Photon contributions from individual point sources are therefore spread over several adjacent pixels, each containing a fraction f of the total photon flux from the source. Apart from being a function of the pixel position, the fractions f depend on the location of a source within its central pixel. A smaller pixel size, i.e., a higher-resolution map, decreases the values of f , corresponding to a relatively larger PSF smoothing. Equation (4) must therefore be corrected for PSF effects. Following MH11, the PSF correction was incorporated by statistical means, considering the average distribution of fractions ρ(f ) among pixels for a given pixel size. To determine ρ(f ), we used Monte Carlo simulations distributing a number of N fiducial point sources at random positions on the sky. The sources were convolved with the detector PSF, and the fractions f i , i = 1, . . . , N pix , were evaluated for each source. The sums of the fractions f i were normalized to 1. We used the effective detector PSF derived from the data set analyzed below, corresponding to the specific event selection cuts used in our analysis. The effective detector PSF was obtained by averaging the detector PSF over energy and spectral index distribution. This is further explained in Section 3. The average distribution function ρ(f ) is then given by where ∆N (f ) denotes the number of fractions in the interval (f, f + ∆f ). The distribution obeys the normalization condition df f ρ(f ) = 1 .
The expected number of m-photon sources in a given pixel corrected for PSF effects is given by (14) Figure 2 depicts the distribution function ρ(f ) derived for the effective PSF of the data set for two different pixel sizes. The function ρ(f ) is also shown assuming a Gaussian PSF with a 68% containment radius resembling the one of the actual PSF. Compared to the Gaussian case, the more pronounced peak of the detector PSF reflects in a strongly peaked ρ(f ) at large flux fractions. Reducing the pixel size, i.e., effectively increasing PSF smoothing (in the sense of this analysis), shifts the peak of ρ(f ) to smaller f . The impact of the large tails of the detector PSF becomes evident at small fractions.

Data Fitting
To fit the model (H) to a given data set (D), we used the method of maximum likelihood (see, e.g., Olive & Particle Data Group (2014) for a review). We defined the likelihood L(Θ) ≡ P (D|Θ, H) in two different ways, which we refer to as L1 and L2 in the following. The likelihood function describes the probability distribution function P of obtaining the data set D, under the assumption of the model (hypothesis) H with a given parameter set Θ.
For a source-count distribution following an MBPL with N b breaks and the previously defined background contributions, the parameter vector is given by (15) containing N Θ = 2N b + 4 free parameters.

Likelihood L1
The L1 approach resembles the method of the simple 1pPDF(see MH11). Given the probability distribution p k for a given Θ, the expected number of pixels containing k photons is ν k (Θ) = N pix p k (Θ). The probability of finding n k pixels with k photons follows a Poissonian (if pixels are considered statistically independent), resulting in the total likelihood function where k max denotes the maximum value of k considered in the analysis.

Likelihood L2
The simple 1pPDF approach can be improved by including morphological information provided by templates. The L2 approach defines a likelihood function that depends on the location of the pixel. The probability of finding k photons in a pixel p is given by p (p) k for a given parameter vector Θ. We emphasize that now the data set comprises the measured number of photons k p in each pixel p, instead of the n k -histogram considered in L1. For clarity, the function p (p) k is therefore denoted by P (k p ) ≡ p (p) k in the following. The likelihood function for the entire ROI is then given by It should be noted that the L2 approach is a direct generalization of the L1 approach. The 1pPDF approach already provides the PDF for each pixel, and it is thus natural to use the appropriate PDF for each pixel instead of using the average one and comparing it with the n khistogram. The L2 approach can then be seen as building a different n k -histogram for each pixel, comparing it with the appropriate p k distribution and then joining the likelihoods of all the pixels together in the global L2 one. The fact that for each pixel the n k -histogram actually reduces to a single count does not pose a matter-ofprinciple problem.

Bayesian Parameter Estimation
The sampling of the likelihood functions L 1 (Θ) and L 2 (Θ) is numerically demanding and requires advanced Markov Chain Monte Carlo (MCMC) methods to account for multimodal behavior and multiparameter degeneracies. We used the multimodal nested sampling algorithm MultiNest 9 (Feroz & Hobson 2008;Feroz 9 Version v3. 8, 20148, October et al. 20098, , 2013 to sample the posterior distribution P (Θ|D, H). The posterior is defined by Bayes's theorem as P (Θ|D, H) = L(Θ)π(Θ)/Z, where Z ≡ P (D|H) is the Bayesian evidence given by and π(Θ) is the prior. MultiNest was used in its recommended configuration regarding sampling efficiency. For our analysis setups, we checked that sufficient sampling accuracy was reached using 1,500 live points with a tolerance setting of 0.2. Final acceptance rates typically resulted in values between 5% and 10%, while the final samples of approximately equal-weight parameter space points consisted of about 10 4 points. From the marginalized one-dimensional posterior distributions, for each parameter we quote the median, and the lower and upper statistical uncertainties were derived from the 15.85% and 84.15% quantiles, respectively. In the case of log-flat priors (see below), we assumed the marginalized posterior distribution to be Gaussian for deriving single-parameter uncertainty estimates in linear space. The derivation of uncertainty bands of the dN/dS fit exploited the same method but using the full posterior.
Priors were chosen to be flat or log flat, depending on the numerical range required for a parameter. Details are discussed in Section 4.3.3.

Frequentist Parameter Estimation
Bayesian parameter estimates from the posterior distributions are compared to parameter estimates employing the frequentist approach. The MCMC method intrinsically provides samples of a posterior distribution that depends on the prior. Nonetheless, if the number of samples is sufficiently high such that also the tails of the posterior are well explored, it can be assumed that the final sample reasonably explored the likelihood function. Profile likelihood functions (see, e.g., Rolke et al. 2005) can be built from the posterior sample. In particular, we built the profile likelihood of the dN/dS fit and one-dimensional profile likelihoods for each parameter. We quote the maximum likelihood parameter values and 68% confidence level (CL) intervals derived under the assumption that the profiled −2 ln L follows a chi-squared distribution with one degree of freedom, i.e., we quote the values of the parameters for which −2∆ ln L = 1. 10 The advantage of profile likelihood parameter estimates is that they are prior independent.

FERMI -LAT DATA
The analysis is based on all-sky gamma-ray data that were recorded with the Fermi -LAT 11 within the first 6 years of the mission. 12 Event selection and processing were performed with the public version of the Fermi Science Tools (  We used Pass 7 Reprocessed (P7REP) data along with P7REP V15 instrument response functions. The application of the analysis method presented here is restricted to the energy bin between E min = 1 GeV and E max = 10 GeV. The lower bound in energy was motivated by the size of the PSF, which increases significantly to values larger than 1 • for energies below 1 GeV (Ackermann et al. 2012). The significant smoothing of point sources caused by a larger PSF may lead to large uncertainties in this analysis (see Section 2.3). Effects of a possible energy dependence of dN/dS are mitigated by selecting an upper bound of 10 GeV.
Data selection was restricted to events passing CLEAN event classification, as recommended for diffuse gammaray analyses. We furthermore required FRONT-converting events, in order to select events with a better PSF and to avoid a significant broadening of the effective PSF. Contamination from the Earth's limb was suppressed by allowing a maximum zenith angle of 90 • . We used standard quality selection criteria, i.e., DATA QUAL==1 and LAT CONFIG==1, and the rocking angle of the satellite was constrained to values smaller than 52 • . The data selection tasks were carried out with the tools gtselect and gtmktime.
The resulting counts map was pixelized with gtbin using the equal-area HEALPix pixelization scheme (Górski et al. 2005). The resolution of the discretized map is given by the pixel size, θ pix = Ω pix . For the statistical analysis employed here, the optimum resolution is expected to be of the order of the PSF: while undersampling the PSF leads to information loss on small-scale structures such as faint point sources, oversampling increases the statistical uncertainty on the number of counts per pixel. We thus compared two choices for the map resolution, where κ denotes the HEALPix resolution parameter: κ = 6 (N side = 64), 14 corresponding to a resolution of ∼ 0.92 • , and κ = 7 (N side = 128), corresponding to a resolution of ∼ 0.46 • . These choices slightly undersample or oversample the actual PSF, respectively.
We used gtltcube and gtexpcube2 to derive the exposure map as a function of energy. The lifetime cube was calculated on a spatial grid with a 1 • spacing. The 14 The number of pixels of the all-sky map is given by N pix = 12N 2 side ; N side can be obtained from the resolution parameter by N side = 2 κ . 8.35e+10 1.26e+11 cm 2 s Figure 4. Fermi-LAT exposure map of the 6-year data set, averaged over the energy band between 1 GeV and 10 GeV. The exposure map is divided into 20 equally spaced regions. For each region, the mean exposure is plotted. The coordinate system matches the one used in Figure 3. The color mapping is linear.
exposure map imposed a spatial grating of 0.125 • (in Cartesian projection) and the same energy binning as used in the Galactic foreground template. The map was projected into HEALPix afterwards. The statistical analysis requires a careful correction for effects imposed by the PSF; see Section 2.3 for details. The PSF of the data set was calculated with gtpsf for a fiducial Galactic position (l, b) = (45 • , 45 • ) as a function of the displacement angle θ and the energy E. We checked that changes of the PSF at other celestial positions were negligible. Given that the PSF strongly depends on energy, analyzing data in a single energy bin requires appropriate averaging. The effective PSF of the data set was calculated by weighting with the energydependent exposure and power-law type energy spectra where E(E) = E (p) (E) ROI denotes the exposure averaged over the ROI. An average spectral index Γ = 2.4 was assumed.
The analysis presented in this article was carried out for high Galactic latitudes |b| ≥ 30 • , aiming at measuring the source-count distribution and compositon of the extragalactic gamma-ray sky. For |b| ≥ 30 • (corresponding to f ROI = 0.5), the photon counts map comprises 862,459 events distributed in 24,576 pixels (κ = 6). The counts map, with a minimum of 5 events per pixel and a maximum of 4,101 events, is shown in Figure 3.
The energy-averaged exposure map of the data set is shown in Figure 4 for the full sky, divided into 20 equal-exposure regions (see Equation (9)). The fullsky (unbinned) exposure varies from 8.22 × 10 10 cm 2 s to 1.27 × 10 11 cm 2 s. The mean of the energy-averaged exposure is 9.18 × 10 10 cm 2 s for |b| ≥ 30 • .

ANALYSIS ROUTINE
The following section is dedicated to details of the analysis method and to the analysis strategy developed in this article. The analysis aims at measuring (i) the contribution from resolved and unresolved gamma-ray point sources to the EGB, (ii) the shape of their source-count distribution dN/dS, and (iii) the resulting total composition of the gamma-ray sky, in the energy band between 1 GeV and 10 GeV. The restriction to Galactic latitudes |b| ≥ 30 • provides a reasonable choice for ensuring that the dominant source contributions are of extragalactic origin.

Expected Sensitivity
The source-population sensitivity of the method can be estimated from the theoretical framework discussed in Section 2. By definition, the total PDF incorporates background components as populations of 1-photon sources (see Equation 6). Sources contributing on average two photons per pixel should be clearly distinguishable from background contributions. The limiting sensitivity on the point-source flux is thus given by the inverse of the average exposure, yielding a value of S sens 2.31×10 −11 cm −2 s −1 for a pixel size corresponding to resolution κ = 6. This value gives a back-of-theenvelope estimate of the sensitivity to the point-source population, while the actual sensitivity additionally depends on quantities such as the unknown shape of the source-count distribution, the relative contribution from foreground and background components, and the number of evaluated pixels N pix (i.e., the Galactic latitude cut). The actual sensitivity will be determined from a data-driven approach in Section 5, as well as from simulations in Appendix D.
In comparison, the sensitivity of the 3FGL catalog drops at a flux of ∼ 2.2 × 10 −10 cm −2 s −1 for the energy band between 1 GeV and 10 GeV. 15 Additional sensitivity can be achieved to lower fluxes by correcting for point-source detection efficiency. However, determining the point-source detection efficiency is nontrivial. The catalog detection procedure needs to be accurately reproduced with Monte Carlo simulations and the method is not completely free from assumptions regarding the properties of the unresolved sources. A clear advantage of the method employed here is, instead, that no detection efficiency is involved. As indicated by the value of S sens , we will see that this analysis increases the sensitivity to faint point-source populations by about one order of magnitude with respect to the 3FGL catalog.

Analysis Setup
The L2 approach emerged to provide significantly higher sensitivity than the L1 approach, as a consequence of the inclusion of spatial information. We will thus use the second method L 2 (Θ) as our reference analysis in the remainder. We will nonetheless present in the main text a comparison of the two approaches, showing that they lead to consistent results.
All pixels in the ROI were considered in the calculation of the likelihood. The upper bound on the number of photon counts per pixel, k max , as used in Equation (16) was always chosen to be slightly larger than the maximum number of counts per map pixel.

Source-count Distribution Fit
15 See Section 4.2 in Acero et al. (2015). The catalog threshold has been rescaled to the 1 GeV to 10 GeV energy band assuming an average photon index of 2.4.
The source-count distribution dN/dS was parameterized with the MBPL defined in Equation (10). For readability, the following terminology will be used in the remainder: the source-count distribution is subdivided into three different regimes, defined by splitting the covered flux range S into three disjoint intervals, 0, 10 −10 cm −2 s −1 : faint-source region, 10 −10 , 10 −8 cm −2 s −1 : intermediate region, The quantity S cut corresponds to a high cutoff flux of the source-count distribution. The observational determination of S cut is limited by cosmic variance, and a precise value is therefore lacking. Unless stated otherwise, we chose a cutoff value S cut = 10 −6 cm −2 s −1 , which is almost one order of magnitude higher than the flux of the brightest source listed in the 3FGL catalog within the ROIs considered in this work (see Section 5.1). The stability of this choice was checked by comparing with S cut = 10 −5 cm −2 s −1 .
In the following, we describe our strategy to fit the dN/dS distribution to the data. A validation of the analysis method with Monte Carlo simulations is described in Appendix D.

Parameters of dN/dS
Normalization -The reference normalization flux S 0 was kept fixed during the fit. A natural choice for S 0 would be the flux where the uncertainty band of the dN/dS reaches its minimum (pivot flux). In this way, undesired correlations among the fit parameters are minimized. We refrained from a systematic determination of the pivot point, but we instead fixed S 0 to a value of S 0 = 3 × 10 −8 cm −2 s −1 after optimization checks. We checked for robustness by varying S 0 within the range [0.1 S 0 , S 0 ], obtaining stable results. 16 Remaining parameter degeneracies were handled well by the sampling.
Number of Breaks -Previous works investigating the gamma-ray dN/dS distribution with cataloged sources concluded that the dN/dS distribution above |b| > 10 • is well described by a broken power law down to a flux of ∼ 5 × 10 −10 cm −2 s −1 , with a break at (2.3 ± 0.6) × 10 −9 cm −2 s −1 (Abdo et al. 2010b). The following analysis increases the sensitivity to resolving point sources with a flux above ∼ 2 × 10 −11 cm −2 s −1 and provides a significantly smaller statistical uncertainty. We therefore parameterized dN/dS with up to three free breaks (N b ≤ 3), in order to find the minimum number of breaks required to properly fit the data. In the case of N b = 3, one break was placed in the bright-source region, a second in the intermediate region, and the last one in the faint-source region; see Section 4.3.3 for details. We compared these results with setups reducing N b to one or two free breaks, to investigate stability and potential shortcomings in the different approaches.

Fitting Techniques
16 Given that the choice of S 0 turns out to be larger than the position of the first break, we note that increasing the interval to larger fluxes is not required.
We employed three different techniques of fitting the dN/dS distribution to the data, in order to investigate the stability of the analysis and to study the sensitivity limit. The third technique, which we refer to as the hybrid approach, is a combination of the two other techniques. This hybrid approach proved to provide the most robust results.
MBPL Approach -The MBPL approach comprises fitting a pure MBPL with a number of N b free break positions. The total number of free parameters is given by N Θ = 2N b + 4 (including free parameters of the background components). The parameters of the MBPL are sampled directly.
Node-based Approach -The complexity of the parameter space, including degeneracies between breaks and power-law indices, can be reduced by imposing a grid of N nd fixed flux positions, which we refer to as nodes S ndj , where j = 0, 1, . . . , N nd − 1. Nodes are counted starting from the one with the highest flux, in order to maintain compatibility with the numbering of breaks in the MBPL described in Equation (10). The free parameters of the source-count distribution correspond to the values of dN/dS at the positions of the nodes, i.e., A ndj = dN/dS (S ndj ). The index of the power-law component below the last node, n f , is kept fixed in this approach.
The parameter set {A ndj , S ndj , n f } can then be mapped to the MBPL parameters using Equation (10), i.e., the dN/dS distribution between adjacent nodes is assumed to follow power laws. Technically, it should be noted that S cut ≡ S nd0 in this case. A choice of N nd nodes therefore corresponds to choosing an MBPL with N nd − 1 fixed breaks. The quantity A S is to be calculated at a value close to the decorrelation flux to ensure a stable fit. The total number of free parameters is given by N Θ = N nd + 2.
While this technique comes with the advantage of reducing the complexity of the parameter space, the choice of the node positions is arbitrary. This can introduce biases between nodes and can thus bias the overall dN/dS fit. The node-based approach is further considered in Appendix C.
We note that a similar approach has been recently used by Vernstrom et al. (2014) for measuring the sourcecount distribution of radio sources.
Hybrid Approach -The hybrid approach combines the MBPL approach and the node-based approach. Free break positions as used in the MBPL approach are required to robustly fit the dN/dS distribution and to determine the sensitivity; see Section 5 for details. Fitting a pure MBPL, however, was found to underestimate the uncertainty band of the fit at the lower end of the faint-source region. In addition, the fit obtained from the Bayesian posterior can suffer a bias for very faint sources, as demonstrated by Monte-Carlo simulations in Appendix D. We therefore chose to incorporate a number of nodes around the sensitivity threshold of the analysis, resolving the issues of the MBPL approach.
The hybrid approach is characterized by choosing a number N h b of free breaks, a number N h nd of nodes, and the index of the power-law component below the last node, n f . We note that the lower limit of the prior of the last free break S bN h b technically imposes a fixed node S nd0 , given that the first free node S nd1 is continuously connected with a power law to the MBPL component at higher fluxes. The setup corresponds to choosing an MBPL with N h b + N h nd + 1 breaks, with the last ones at fixed positions. The total number of free parameters in the hybrid approach is

Priors
We used log-flat priors for the normalization A S , the nodes A ndj , the breaks S bj , and the isotropic diffuse background flux F iso , while the indices n j and the normalization of the Galactic foreground map A gal were sampled with flat priors. Prior types and prior ranges are listed in Table 1 for the MBPL and hybrid approaches. In general, priors were limited to physically reasonable ranges. Prior ranges were chosen to cover the posterior distributions well.
In particular, data from the 3FGL catalog motivate that S 2 dN/dS 10 −11 cm −2 s −1 deg −2 in the intermediate region; see Section 5. The range of the prior for A S was therefore adjusted to cover the corresponding interval between 3 × 10 −12 cm −2 s −1 deg −2 and 8 × 10 −11 cm −2 s −1 deg −2 at least (assuming an index of 2). The ranges of the priors for the node normalizations were chosen similarly, but reducing the lower bound to a value of ∼ 10 −12 cm −2 s −1 deg −2 .
The ranges of the priors for the breaks were chosen to connect continuously and not to overlap, preserving a well-defined order of the break points. For both the MBPL and the hybrid approach, the upper bound of the first break S b1 approximately matched the bright end of the 3FGL data points (excluding the brightest source). It is advantageous to keep the prior range for the first break sufficiently small, in order to reduce a possible bias of the intermediate region by bright sources (mediated through the index n 2 ). For the MBPL approach, the lower bound of the last break was chosen almost two orders of magnitude below the sensitivity estimate of S sens 2 × 10 −11 cm −2 s −1 , to fully explore the sensitivity range. In the case of three breaks, the lower bound of the intermediate break was selected to match the sensitivity estimate. For the hybrid approach, the lower bound of the last free break was set to ∼ S sens /2 . We comment on the choice of the nodes in Section 5.2.
Index ranges were selected according to expectations, allowing enough freedom to explore the parameter space. The stability of these choices was checked iteratively. For the MBPL approach, the lower bound of the last index allowed for a sharp cutoff of the dN/dS distribution. For the hybrid approach, the index n f was fixed to a value of −10, introducing a sharp cutoff manually. This choice will be motivated in Section 5.1.
As discussed in Section 2.2.3, A gal is expected to be of order unity. The selection of the prior boundaries for F iso was based on previous measurements (see Ackermann et al. 2015b) and was further motivated iteratively.
Prior ranges reported in Table 1 are further discussed in Section 5.

Exposure Correction: Nexp
The results were checked for robustness with respect to variations of N exp (see Section 2.1). We found that Table 1 Prior Ranges Note.
-Prior types and ranges used for the different dN/dS and background parameterizations investigated in this work. Either the dN/dS was parameterized with a pure MBPL, or the MBPL was extended with a node (hybrid approach). For the node-based approach see Appendix C. The table lists the prior ranges used for the MBPL and hybrid approaches, given a number of N b or N h b free breaks, respectively. Priors listed in the first panel "Generic" were used in both setups identically. Ellipses indicate parameters not present in the specific model. a The normalization A S is given in units of 10 7 s cm 2 sr −1 , while the normalization of the node A nd1 is given in 10 14 s cm 2 sr −1 . The normalizations refer to S 0 = 3 × 10 −8 cm −2 s −1 . The breaks S b1 , S b2 , and S b3 , as well as the node position S nd1 , are given in units of cm −2 s −1 . The diffuse flux component F iso is given in units of 10 −7 cm −2 s −1 sr −1 . All other quantities are dimensionless. b The lower prior limit poses an exception to our requirement of n 1 > 2 (see Section 2.2.1). Because only a single break is allowed in this model, a suitable prior coverage must include the cases of a faint break with a hard index n 1 < 2, and a break in the bright-source region with a consequently softer index n 1 > 2.
the choice of this parameter is critical for a correct recovery of the final result, and it is closely related to the sensitivity of the analysis. In particular, small values 5 were found insufficient. Results were stabilized by using at least N exp = 15 contours, and we tested that increasing up to N exp = 40 did not have further impact. Insufficient sampling of the exposure (i.e., small values of N exp ) was seen to affect the faint end of the dN/dS by introducing an early cutoff and attributing a larger flux to the isotropic component. At the same time, the bestfit likelihood using small N exp values was significantly smaller than the one obtained choosing larger values, indicating that indeed the sampling was insufficient. As a final reference value we chose N exp = 20.

3FGL Catalog Data
The results are compared to the differential (dN/dS) and integral (N (> S)) source-count distributions derived from the 3FGL catalog for the same energy band and ROI. The method of deriving the source-count distribution from catalog data is described in Appendix B.

APPLICATION TO THE DATA
In this section, a detailed description and discussion of the data analysis and all setups chosen in this article are given. Final results are summarized in Section 8.
The data were fit by employing the MBPL approach and the hybrid approach consecutively. The use of the hybrid approach was mostly chosen to inspect the uncertainties in the faint-source region. It should be emphasized that the prior of the last free break and the position of the node depend on the results obtained with the MBPL approach.
All analyses were carried out using two different pixel sizes, i.e., HEALPix grids of order κ = 6 (N side = 64) and κ = 7 (N side = 128). Details are discussed in Section 3. We chose κ = 6 as a reference, due to the expected sensitivity gain. All parameters were stable within their uncertainty bands against changes to κ = 7. Results using κ = 7 are shown in Section 7.1.

MBPL Approach
The MBPL fit was employed using the priors as discussed in Section 4.3.3. The results are shown in Table  2 and Figure 5.
The source-count distribution was parameterized with one, two, and three free breaks. Table 2 lists all bestfit values and statistical uncertainties obtained for individual fit parameters, in addition to the corresponding likelihoods of the best-fit solutions. Single-parameter uncertainties can be large in general, given that correlations were integrated over. Comparing Bayesian (posterior) and frequentist (profile likelihood 17 ) parameter   S b3 · · · · · · · · · · · · 4.6 +4. Note. -Best-fit values and statistical uncertainties (68.3% CL) obtained for a pure MBPL fit to the data (MBPL approach). The table compares dN/dS fits with one, two, and three free breaks. Both the parameter values obtained from the Bayesian posterior and the values derived from the profile likelihood (PL) are given. The last three rows list the values of the L 1 and L 2 likelihoods for the best-fit results. The value of the Bayesian evidence Z is given in addition. Ellipses indicate parameters not present in the specific model. a A S is given in units of 10 7 s cm 2 sr −1 . We remind that the values correspond to a flux normalization constant of S 0 = 3 × 10 −8 cm −2 s −1 . The breaks S b1 , S b2 , and S b3 are given in units of 10 −8 , 10 −10 , and 10 −12 cm −2 s −1 , respectively. The units of the diffuse flux component F iso are 10 −7 cm −2 s −1 sr −1 . All other quantities are dimensionless.
estimates, best-fit values match within their uncertainties. Figure 5 shows the best-fit results and corresponding statistical uncertainty bands for the dN/dS distributions parameterized with two and three free breaks. We can see that there is good agreement between the dN/dS distributions derived from the Bayesian posterior (solid tainties are given in the caption of Figure 6(b).
black line and green band) and the dN/dS fits derived from the profile likelihood (dashed black line and blue band): they match well within their uncertainty bands. The uncertainty given by the profile likelihood is larger than the band from the posterior in all cases. The frequentist uncertainty estimates can therefore be considered more conservative. In common, the statistical uncertainty bands of the dN/dS fits obtained here are small compared to fits employing catalog points only (see Abdo ). This directly reflects the fact that the method is independent of source-detection or binning effects. The smallest statistical uncertainty appears to be around a flux of ∼ 10 −9 cm −2 s −1 . As shown in Table 2, the fit of the simplest dN/dS model with only a single break prefers a break at low fluxes, i.e., at ∼ 10 −11 cm −2 s −1 . Below that break, the dN/dS cuts off steeply. The source-count distribution in the entire flux range above that break was fit with the single power-law component, with an index of n 1 = 2.03 ± 0.02. We found that adding a break at higher fluxes, i.e., parameterizing dN/dS with two free breaks, instead improved the fit with a significance of ∼3σ. Here the bright-source region is resolved with a break at ∼ 2 × 10 −8 cm −2 s −1 . The region between the two breaks (faint-source region and intermediate region) is compatible with an index of n 2 = 1.97±0.03, while the index in the bright-source region n 1 = 3.1 +0.7 −0.6 is softer (see Figure 5).
The intermediate region is populated with numerous sources contributing a comparably large number of photons. Given the high statistical impact of these sources, it was found that a fit of the faint-source and intermediate regions with only a single power-law component can be significantly driven by brighter sources of the intermediate region. We therefore extended the dN/dS model to three free breaks, properly investigating possible features in the faint-source region below ∼ 3 × 10 −9 cm −2 s −1 . We found that the model comprising three free breaks is not statistically preferred against the two-break model (see Table 2). Furthermore, the three-break dN/dS distribution is consistent with the previous scenario within uncertainties (see Figure 5). Differences between the best fit from Bayesian inference and the best fit given by the maximum likelihood are not statistically significant.
It can be seen in Figure 5 that the source-count distribution as resolved by the 3FGL catalog (red data points; see Section 4.3.5) in the intermediate and the brightsource regions is well reproduced with both the two-break and the three-break fits. Again, we emphasize that this analysis is independent of catalog data, which are shown in the plot for comparison only.
From the MBPL approach, we therefore conclude that parameterizing dN/dS with two free breaks is sufficient to fit the data. The index n 2 = 1.97±0.03, characterizing the intermediate region of dN/dS, is determined with exceptionally high precision (∼2%), originating from the high statistics of sources populating that region. The accuracy of the Galactic foreground normalization A gal fit is at the per mil level.
We found that the fit prefers a source-count distribution that continues with an almost flat slope (in S 2 dN/dS representation) in the regime of unresolved sources, i.e., faint sources not detected in the 3FGL catalog. A strong cutoff was found at fluxes between ∼ 5 × 10 −12 cm −2 s −1 and ∼ 10 −11 cm −2 s −1 . This cutoff, however, falls well within the flux region where this method is expected to lose sensitivity and where the uncertainty bands widen. It should thus be considered with special care. Indeed, Monte Carlo simulations were used to demonstrate that such a cutoff can originate either from the sensitivity limit of the analysis or from an intrinsic end of the source-count distribution (see Appendix D for details). In the former case, possible point-source contributions below the cutoff are consistent with diffuse isotropic emission, and the fit therefore attributes them to F iso .
It was found that the uncertainty band below the cutoff can be underestimated due to lacking degrees of freedom in the faint-source end. Moreover, simulations revealed that the fit obtained from the Bayesian posterior can be biased in the regime of very faint sources. We therefore chose to improve the fit procedure by using the hybrid approach in Section 5.2.
Sampling -The triangle plot of the Bayesian posterior and the corresponding profile likelihood functions are shown in Figure 6 for parameterizing dN/dS with three free breaks.
It can be seen that the marginalized posterior distributions are well defined. We attenuated strong parameter degeneracies by adapting the normalization constant S 0 to the value quoted in the previous section. It becomes evident from the posteriors that the breaks S b2 and S b3 tended to merge to a single break; this is supplemented by the flatness of their profile likelihoods. It therefore explains the previous observation that adding a third break is not required to improve the fit of the data.

Hybrid Approach
We improved the analysis by applying the hybrid approach consecutively. Priors are discussed in Table 1: In particular, the region around the sharp cutoff revealed by the MBPL approach was parameterized with a node placed at A nd1 = 5×10 −12 cm −2 s −1 . 18 The lower bound of the prior of the last free break was set to ∼ S sens /2 . The cutoff was introduced manually by fixing the index of the power-law component describing fluxes smaller than A nd1 to n f = −10.
The fit was carried out with dN/dS parameterizations comprising one, two, and three free breaks. Figure 7 and Table 3 summarize the results. The differential dN/dS distributions fitting the data best are shown in the left column of the figure. In the right column, the corresponding integral source-count distributions N (> S) are compared to 3FGL catalog data, providing another reference for investigating the precision of the fit.
In the bright-source and intermediate regions, the results obtained with the MBPL approach and with the hybrid approach are consistent among each other within their uncertainties. As expected, the determination of the uncertainty bands in the faint-source region improved in the hybrid fit, given the further degree of freedom allowed. In all three scenarios (N h b = 1, 2, 3), the fits reproduce well the differential and the integral source-count distributions from the 3FGL catalog within uncertainties.
Comparing the three dN/dS models, we find that none are statistically preferred by the data; see Table 3. The fit of the model with only a single free break consistently placed the break in the bright-source region, given that the cutoff in the faint-source region is effectively accounted for by the node. As argued in the previous section, in this case the fit of the intermediate and 18 The value approximates the faint cutoff positions obtained from the posterior of the MBPL fit. Note. -Best-fit values and statistical uncertainties (68.3% CL) obtained with the hybrid approach. The table compares dN/dS fits with one, two, and three free breaks. Both the parameter values obtained from the Bayesian posterior and the values derived from the profile likelihood (PL) are given. The last three rows list the values of the L 1 and L 2 likelihoods for the best-fit results. The value of the Bayesian evidence Z is given in addition. Ellipses indicate parameters not present in the specific model. a A S is given in units of 10 7 s cm 2 sr −1 . We remind that the values correspond to a flux normalization constant of S 0 = 3 × 10 −8 cm −2 s −1 . The breaks S b1 , S b2 , and S b3 are given in units of 10 −8 , 10 −10 , and 10 −12 cm −2 s −1 , respectively. The normalization A nd1 of the node is given in units of 10 14 s cm 2 sr −1 . The diffuse flux component F iso is listed in units of 10 −7 cm −2 s −1 sr −1 . All other quantities are dimensionless.
faint-source regions of dN/dS was driven by the high statistical impact of the relevant brighter sources, yielding a small uncertainty band also for faint sources (see Figure 7(a)). To address this issue, we extended the model with two additional free breaks (N h b = 2, 3), leading to consistent uncertainty bands that were stabilized by the additonal degrees of freedom added in the intermediate and faint-source regions (see Figures 7(c) and 7(e)). Because the three-break fit is not statistically preferred against the two-break fit, we conclude that two free breaks and a faint node are sufficient to fit the data properly. A comparison with the maximum likelihood values for the MBPL fits in Table 2 reveals also no statistical preference for the hybrid result over the MBPL result, confirming that the data are not sensitive enough to distinguish point sources below the last node from a purely diffuse isotropic emission. Figure 8 compares the best-fit model 1pPDF distributions to the actual pixel-count distribution of the data set. We plot the results for both the Bayesian posterior and the maximum likelihood fits. The residuals (data − model)/ √ data are shown in addition. It can be seen that the pixel-count distribution is reproduced well. A comparison with a simple chi-squared statistic, evaluating the best-fit results using the binned histogram only, leads to reduced chi-squared values (χ 2 /dof) between 0.89 and 0.92 .
The triangle plot of the Bayesian posterior and the single-parameter profile likelihood functions are shown in Figure 9 for the dN/dS fit with two free breaks and a node.
The stability of the MBPL and hybrid approaches can be further demonstrated by comparing the respective triangle plots (see Figures 6(a) and 9(a)): the posteriors of parameters corresponding to each other in both approaches are substantially equal, with the exception of n 3 . It can be seen that the choice of the node in the hybrid approach stabilized the posterior of n 3 . We have therefore shown that the MBPL and hybrid approaches lead to comparable results except in the faint-source flux region, where the latter improves the determination of the uncertainty bands.

How many breaks?
Both the MBPL approach and the hybrid approach single out a best-fit dN/dS distribution that is consistent with a single broken power law for integral fluxes in the resolved range above S sens 2 × 10 −11 cm −2 s −1 . Although two breaks are preferred to properly fit the entire flux range, the second break found with the MBPL approach in the faint-source region is consistent with a sensitivity cutoff. Instead, in the hybrid approach, the second break is needed for a viable determination of the uncertainty band.
To further describe the physical dN/dS distribution at low fluxes, we therefore derived an upper limit on the position of a possible intrinsic second break S b2 . The uncertainty band obtained with the hybrid approach for N h b = 2 was used. In general, an intrinsic second break would have been present if the power-law indices n 2 and n 3 changed significantly by a given difference |n 2 − n 3 | > ∆n 23 . We exploited the full posterior to derive upper limits on S b2 by assuming given ∆n 23 values between 0.1 and 0.7, in steps of 0.1. In detail, the upper limits S UL b2 at 95% CL were obtained from the marginalized posterior P (S b2 |D, H), after removing all samples not satisfying  the given |n 2 − n 3 | constraint: (20) where π L (S b2 ) = 10 −11 cm −2 s −1 is the lower bound of the prior for S b2 . Frequentist upper limits were calculated from the profile likelihood, constructed from the same posterior as used in Equation (20), by imposing −2∆ ln L = 2.71 for 95% CL upper limits. The upper limits are shown in Figure 10. In consistency with the uncertainty bands derived in the previous section, S UL b2 decreases monotonically as a function of ∆n 23 , until the sensitivity limit of the analysis is reached. Assuming a fiducial index change of ∆n 23 = 0.3, we find that a possible second break of dN/dS is constrained to be below 6.4 × 10 −11 cm −2 s −1 at 95% CL. The corresponding frequentist upper limit is 1.3 × 10 −10 cm −2 s −1 .

Composition of the Gamma-ray Sky
The method allows decomposing the high-latitude gamma-ray sky (|b| ≥ 30 • ) into its individual constituents. The integral flux F ps contributed by point sources was derived by integrating the posterior sam-ples of S dN/dS in the range [0, S cut ], which effectively corresponds to the interval [S nd1 , S cut ] due to the steep cutoff below the node S nd1 . Results are presented in Table 4, comparing both Bayesian and frequentist estimates. The profile likelihood for F ps is shown in Figure 11. The integral flux from point sources is determined as F ps = 3.9 +0.3 −0.2 × 10 −7 cm −2 s −1 sr −1 , thus with an uncertainty less than 10%. 19 The contribution from Galactic foreground emission F gal was obtained accordingly by integrating the template (see Section 2.2.3), including the fit results for the normalization A gal (see Figure 9). The isotropic background emission F iso was sampled directly.
For convenience, individual components can be expressed as fractions q of the total map flux F tot . The fractions are listed in Table 4. We found that the highlatitude gamma-ray emission between 1 GeV and 10 GeV is composed of (25 ± 2)% point-source contributions, (69.3 ± 0.7)% Galactic foreground contributions, and (6 ± 2)% isotropic diffuse background emission.
Even if not indicated by Figures 6(a) and 9(a), remain-  Figure 11. Profile likelihood of the integral point-source flux Fps, obtained using the hybrid approach with two free breaks and κ = 6. Line styles and colors are as in Figure 6(b). Note. -The table lists Bayesian (posterior) and frequentist (PL) estimates for the three flux contributions discussed in the text. Fluxes are given in units of 10 −7 cm −2 s −1 sr −1 . The ratios qps, q gal , and q iso refer to the total map flux Ftot = i F i , which has been consistently derived from the posterior and profile likelihood, respectively. For comparison, the integral fluxes of all sources listed in the 2FGL and 3FGL catalogs for |b| ≥ 30 • are also given. The last row lists an estimate of the flux contributed by residual cosmic rays.
ing degeneracies between an isotropic Galactic component accounted for in the template and the F iso parameter considered in this analysis might be present.
The flux contribution from point sources can be compared to the flux of all sources resolved in the 3FGL catalog (for |b| ≥ 30 • ; see Table 4). From the difference F ps − F 3FGL cat we conclude that a flux of 1.4 +0.3 −0.2 × 10 −7 cm −2 s −1 sr −1 between 1 GeV and 10 GeV can be attributed to originate from so far unresolved point sources. With regard to the IGRB flux measured by Ackermann et al. (2015b), we could therefore clarify between 42% and 56% of its origin between 1 GeV and 10 GeV. 20 Residual Cosmic Rays -The sum of the values F iso = (0.9 ± 0.2) × 10 −7 cm −2 s −1 sr −1 and F ps = 3.9 +0.3 −0.2 × 10 −7 cm −2 s −1 sr −1 listed in Table 4 can be compared with the EGB derived in Ackermann et al. (2015b). In the energy range between 1 GeV and 10 GeV this amounts to values between 4.7 × 10 −7 cm −2 s −1 sr −1 and 6.4 × 10 −7 cm −2 s −1 sr −1 , including systematics in the Galactic diffuse modeling; these values compare well with the total F iso + F ps found here.
However, the truly diffuse isotropic background emission F iso incorporates residual cosmic rays (CRs) not rejected by analysis cuts (see Ackermann et al. 2015b), while for the EGB derived in Ackermann et al. (2015b) the CR contamination has been accounted for and subtracted. The level of residual CR contamination in the P7REP CLEAN selection used in this work has been estimated to be between 15% and 20% of the measured IGRB flux above 1 GeV (see Figure 28 in Ackermann et al. 2012), thus amounting to about 5-7 × 10 −8 cm −2 s −1 sr −1 . 20 The IGRB obtained by Ackermann et al. (2015b) in the 1 GeV to 10 GeV energy band is between ∼ 3.2 × 10 −7 cm −2 s −1 sr −1 and ∼ 4.3 × 10 −7 cm −2 s −1 sr −1 , including systematic uncertainties of the Galactic foreground modeling. Note that this measurement refers to the 2FGL catalog, which has been used for subtracting resolved sources from the EGB. We therefore attribute a flux of 1.8 +0.3 −0.2 × 10 −7 cm −2 s −1 sr −1 to unresolved point sources in this IGRB measurement (using F 2FGL cat as quoted in Table 4).

ANISOTROPY
Complementary to the 1pPDF, the anisotropy (or autocorrelation) probes unresolved point sources (Ackermann et al. 2012a;Cuoco et al. 2012;Di Mauro et al. 2014b;Ripken et al. 2014). The two observables can thus be compared. The anisotropy in a given energy band can be calculated from the dN/dS distribution by where S th is the flux threshold of detected point sources, assumed to be 'sharp' and independent of the photon spectral index of the sources. Indeed, the previous assumption is a good approximation for the 1 GeV to 10 GeV energy band (Cuoco et al. 2012). We thus calculated the predicted anisotropy from the dN/dS distribution measured in this work (hybrid approach, N h b = 2) as a function of the threshold flux S th . Results are shown in Figure 12. To derive the uncertainty band of C P , we sampled the dN/dS from the posterior and calculated C P from each sampling point of the dN/dS parameter space. The uncertainty on C P was then derived using both the Bayesian and the frequentist approaches; see Sections 2.4.3 and 2.4.4. The predicted C P can be compared to the value (1.1 ± 0.1) × 10 −17 (cm −2 s −1 sr −1 ) 2 sr measured in Cuoco et al. (2012) and Ackermann et al. (2012a), using a threshold of about 4-6 × 10 −10 cm −2 s −1 suitable for sources detected in the 1FGL catalog (Abdo et al. 2010a). It can be seen in Figure 12 that the predicted anisotropy is slightly higher than the measured value. This can in part be explained by the approximation of the threshold as a sharp cutoff, as well as a possible systematic underestimate of the measured anisotropy itself (Broderick et al. 2014b). In addition, a possible clustering of point sources at angular scales smaller than the pixel size could in principle be degenerate with the inferred dN/dS distribution, leading to systematically higher anisotropies. The anisotropy of clustering effects is, however, expected to be rather small as compared to the C P values found here, i.e., C cluster >200 10 −20 (cm −2 s −1 sr −1 ) 2 sr for multipoles corresponding to angular scales smaller than the pixel size (e.g., Ando et al. 2007;Cuoco et al. 2015). Clustering can thus be neglected in this analysis. For the moment, we deem the agreement reasonable, and we wait for an updated anisotropy measurement for a more detailed comparison.

SYSTEMATICS
The following section is dedicated to systematic and modeling uncertainties of the analysis framework. In particular, we extensively investigated possible uncertainties due to the chosen pixel size (Section 7.1), statistical effects imposed by bright point sources (Section 7.2), and the Galactic foreground modeling (Section 7.3).

Pixel Size
The results discussed in Section 5 were cross-checked using smaller pixels, i.e., HEALPix order κ = 7, slightly oversampling the effective PSF (see Section 3). All results were stable against the resolution change, given the corresponding uncertainty bands. However, it was found that the enhanced PSF smoothing increased the uncertainty in determining the first break. An example is given in Figure 13, showing the dN/dS distribution obtained with the hybrid approach considering three free breaks and a node. It is demonstrated in Section 7.2 that the increased uncertainty in the bright-source region in turn led to a small bias in determining the indices n 2 and n 3 . Table 5 summarizes fit results that do not become evident in Figure 13. The integral point-source flux F ps slightly decreased with respect to the value obtained for κ = 6, with a corresponding increase of the isotropic background emission F iso , while the sum F ps + F iso remained constant within (single-parameter) statistical uncertainties. This is consistent with resolving fewer point sources due to reduced sensitivity, given that the value of A gal stayed almost the same as found for κ = 6.

Point-source Masking
posterior profile likelihood 3FGL Figure 13. Differential source-count distribution dN/dS obtained with the hybrid approach for a HEALPix grid of order κ = 7. The hybrid approach was carried out allowing three free breaks and a node. Line styles and colors are as in Figure 5.
The presence of bright point sources and the corresponding shape of their source-count distribution may influence the overall fit of the intermediate region and the faint-source region. The strength of a possible bias may also depend on the pixel size.
The level of systematics caused by bright point sources was investigated with point-source masks. To eliminate the influence of bright sources, we removed all pixels including sources with an integral flux larger than or equal to S mask 10 −8 cm −2 s −1 . The value of S mask was chosen to be slightly below the first break determined from the overall fit (see Section 5). Source positions and fluxes were retrieved from the 3FGL catalog. For each source, all pixels included in a circle with a radius of 2.5 • (corresponding to ∼ 6σ psf ) 21 around the cataloged source position were masked in the counts map. We checked that the mask area was sufficiently large by comparing radii between 3σ psf and 7σ psf . Remnant effects became negligible for radii larger than ∼ 5σ psf .
The masked data were fit using the hybrid approach with three free breaks, in order to retain full sensitivity to a possible break in the faint-source region. Priors were chosen as listed in Table 1, with the exception of changing the upper bound of the prior of the first break to S mask . The prior of n 1 was changed accordingly to sample the interval [1.7, 2.3], substantially covering the intermediate region. In addition, the flux normalization constant was fixed to S 0 = 3 × 10 −9 cm −2 s −1 and the upper flux cutoff of dN/dS was set to S cut ≡ S mask .
The results are shown in Figure 14 for a pixelization with resolution parameters κ = 6 and κ = 7. It can be seen that the results are consistent with what was found in Section 5. For κ = 7 we find that the uncertainty band is slightly down-shifted as compared to κ = 6, but bestfit results match well within uncertainties. The value of A gal was determined to be 1.071 +0.004 −0.004 (1.072 +0.005 −0.005 ) for κ = 6 and 1.075 +0.004 −0.004 (1.073 +0.006 −0.004 ) for κ = 7, using the posterior (profile likelihood). The integral flux of the isotropic diffuse background emission F iso was obtained to be 0.9 +0.2 −0.2 (0.8 +0.5 −0.3 ) × 10 −7 cm −2 s −1 for κ = 6 and 1.2 +0.1 −0.2 (1.4 +0.2 −0.4 ) × 10 −7 cm −2 s −1 for κ = 7. The larger value of F iso in the latter case is consistent with the fact of resolving fewer point sources for κ = 7.
We conclude that systematic effects due to bright point sources are dominated by statistical uncertainties. Bright point sources do not affect the determination of the dN/dS broken power-law indices in the intermediate and faint-source regions. For κ = 7, comparing the analyses of full data ( Figure 13) and masked data (Figure 14(b)) indicates that systematic effects slightly increased with enhanced PSF smoothing, but effects on the indices n 2 and n 3 remain rather small.

Galactic Foreground
We checked our results for systematic uncertainties of the Galactic foreground model, considering three different approaches: • Dependence on the Galactic latitude cut. We selected different ROIs, covering regions |b| ≥ b cut . The parameter b cut was varied between 10 • and 70 • , in steps of 10 • .
• Extended Galactic plane mask (GPLL mask). The GPLL mask was generated from the Galactic foreground emission model discussed in Section 2.2.3, by merging mask arrays for |b| < 30 • , a Galactic plane mask removing all pixels above a flux threshold 22 of 10 −6 cm −2 s −1 sr −1 , and mask arrays for the Fermi bubbles and Galactic Loop I (Ackermann et al. 2014;Casandjian et al. 2009;Su et al. 2010).
The GPLL mask is shown in Figure 15(a).
• Dependence on the Galactic foreground model. Given systematic uncertainties of the Galactic foreground model in its entirety, we incorporated a different foreground model as derived for the preceeding Fermi -LAT data release Pass 7, named gal 2yearp7v6 v0.fits 23 . Although mixing different versions of data releases and diffuse models is not generally recommended, the purpose here is to gauge the effect of a model differing in intensity as well as in morphology. The deviations between the two models are shown in Figure 15 Figure 15. (a) Same as Figure 1, but overlaying the GPLL mask (in gray). (b) Difference map P7 − P7REP of the Galactic foreground emission models derived for P7 and P7REP data between 1 GeV and 10 GeV. The difference is given in units of counts pixel −1 , derived with the exposure map as discussed in Section 2.2.3. Projection and mask are the same as in Figure 1. The color bar has been clipped at −3 and 1, in order to visualize small differences. The entire range covered by the data is (−6.5, 12). The color mapping is linear.
The hybrid approach was employed for all setups, choosing three free breaks and a node. The prior setup resembled the one used in Section 5, but prior ranges were extended in particular cases to cover the posterior sufficiently well. The results of the analyses are summarized in Figure 16 and Table 6. We found that all results were stable against the systematic checks. In addition, it should be noted that the catalog (3FGL) data points derived for comparison were well reproduced in all cases.
In the bright-source region, the error band increases almost monotonically with increasing Galactic latitude cut, due to the decreasing number of bright sources present in the ROI. We note that for the 10 • cut the index n 1 = 2.58 +0.23 −0.14 matches well within uncertainties the index deduced by the Fermi Collaboration from 1FGL catalog data (n 1 = 2.38 +0.15 −0.14 , Abdo et al. 2010b) for the same latitude cut and energy band. The first break position, however, was found to be a factor of 2 to 3 larger than in the 1FGL analysis. The index below the first break is n 2 2.
The fits of the faint-source region were stable against changing the Galactic latitude cut. The slopes of the corresponding dN/dS fits match well within uncertaintites for increasing latitude. Uncertainties grow for higher Galactic latitude cuts given less statistics. For lower latitude cuts, Figure 16 indicates an upturn for very faint sources, which is, however, not significant. The stability against the Galactic latitude cut is further supplemented by the integral point-source flux F ps (see Table 6), which remains stable within uncertainties. Table 6 shows that the normalization of the Galactic foreground model, A gal , increases with the latitude cut by ∼10% from 10 • to 50 • , while the integral flux of the isotropic background emission remains constant (∼ 9 × 10 −8 cm −2 s −1 sr −1 ). 24 The increase of A gal thus indicates a gradual mismatch between foreground model and data. Likewise, it can also indicate the presence of  Figure 16. Galactic foreground systematics. The left column shows the dN/dS distributions obtained for different Galactic latitude cuts (10 • , 30 • , and 50 • ), using the hybrid approach with N h b = 3 breaks and a node. Besides the 70 • Galactic latitude cut in the top panel, the right column depicts the fit results using the GPLL mask and the P7 Galactic foreground template model. Line styles and colors are the same as in Figure 5.  Tables 3 and 4. a new component not covered by our analysis setup. We note that a similar behavior has been found in other analyses, including the 3FGL catalog (see, e.g., Figure 25 in Acero et al. 2015). The stability of the results obtained in this article is supplemented by comparing with the GPLL mask and the Pass 7 foreground model (P7 model). The GPLL mask in particular removes the Galactic lobes and Galactic Loop I, known as regions potentially affected by large systematic model uncertainties. Employing the P7 model introduces a different Galactic foreground model in its entirety. As demonstrated in Figure 15(b), the differences between the models exhibit a nontrivial morphology. The pixel distribution of photon-count differences extends to ∼3 for the dominating part of the ROI, i.e., systematics can be expected at the flux level of the sensitivity estimate S sens . The resulting dN/dS distributions and the integral point-source fluxes F ps are consistent within uncertainties. It is to be noted, however, that the integral isotropic background flux F iso increased by a factor of ∼2 for the P7 model. At the same time, F gal decreased, maintaining a stable sum F gal + F iso . We therefore remark that modeling uncertainties can cause F iso to depend on the Galactic foreground model.

CONCLUSIONS
In this article, we have employed the pixel-count distribution (1-point PDF) of the 6-year photon counts map measured with Fermi -LAT between 1 GeV and 10 GeV to decompose the high-latitude gamma-ray sky. This 10 -11 10 -10 10 -9 10 -8 10 -7 10 -6 10 -13 10 -12 10 -11 10 -10 posterior profile likelihood 3FGL Figure 17. Differential source-count distribution dN/dS obtained from 6-year Fermi-LAT data for high Galactic latitudes greater than 30 • . The fit was carried out by employing the hybrid approach with two free breaks and a node at the faint end of the distribution. Line styles and colors are as in Figure 5. The shape of the dN/dS distribution for very faint sources can be further constrained by the fact that the sum of the integral point-source flux and the Galactic foreground contribution must not exceed the total map flux Ftot. Corresponding constraints have been derived assuming the dN/dS distribution obtained from the Bayesian posterior down to the best-fit position of the last free break, i.e., 3 × 10 −11 cm −2 s −1 . Below that value, dN/dS has been extrapolated with a power-law component of varying index. At the boundary of the gray-shaded region the total point-source flux equals Ftot − F gal , requiring a break.
statistical analysis method has allowed us to dissect the gamma-ray sky into three different components, i.e., point sources, diffuse Galactic foreground emission, and a contribution from isotropic diffuse background. The analysis of the simple pixel-count distribution has been improved by employing a pixel-dependent approach, in order to fully explore all the available information and to incorporate the morphological variation of components such as the Galactic foreground emission. A summary of the main results obtained with this analysis follows. The distribution of point sources dN/dS has been fit assuming a multiply broken power law (MBPL approach) with one, two, and three free breaks. A possible bias in obtaining the correct statistical uncertainty band for faint-source contributions has been mitigated by extending the setup with a node, what we called the hybrid approach. Figure 17 summarizes the resulting dN/dS distribution at high Galactic latitudes b greater than 30 • .
We have found that both the MBPL approach and the hybrid approach single out a best-fit source-count distribution for |b| ≥ 30 • that is consistent with a single broken power law for integral fluxes S in the resolved range. Although two-break models are preferred to properly fit the entire flux range covered by the data, the second break found in the MBPL approach in the faint-source region is consistent with a sensitivity cutoff. Instead, in the hybrid approach, the second break is needed for a viable determination of the uncertainty band. The MBPL and hybrid approaches have led to comparable results except in the faint-source flux region, where the latter improved the uncertainty band. For bright sources with an integral flux above the first break at 2.1 +1.0 −1.3 × 10 −8 cm −2 s −1 the dN/dS distribution follows a power law with index n 1 = 3.1 +0.7 −0.5 . Below the first break, the index characterizing the intermediate region and the faint-source region of dN/dS hardens to n 2 = 1.97 +0.03 −0.03 . It is determined with exceptionally high precision (∼2%) thanks to the high statistics of sources populating that region. The fit is consistent with the distribution of individually resolved sources listed in the 3FGL catalog. We have measured dN/dS down to an integral flux of ∼ 2 × 10 −11 cm −2 s −1 , improving beyond the 3FGL catalog detection limit by about one order of magnitude. To further constrain the physical dN/dS distribution at low fluxes, we have derived an upper limit on a possible intrinsic second break from the uncertainty band obtained with the hybrid approach. We have found that a possible second break of dN/dS is constrained to be below 6.4 × 10 −11 cm −2 s −1 at 95% CL, assuming a change of ∆n ≥ 0.3 for the power-law indices below and above that break.
We have checked our results against a number of possible systematic and modeling uncertainties of the analysis framework. Likewise, the behavior of dN/dS has been investigated as a function of the Galactic latitude cut. We have considered Galactic latitude cuts in the interval between 10 • and 70 • . We have found that the faint-source and the intermediate regions of dN/dS are not altered, while the uncertainty band in the bright end becomes larger due to the decreasing number of bright sources in the ROI. At the same time, fitting the overall normalization of the Galactic foreground template has revealed that it significantly increases with higher latitude cuts. This indicates a possible gradual mismatch between the Galactic foreground model and the data at high latitudes, or a missing component not accounted for in our analysis setup. Note, however, that this increase does not affect the obtained dN/dS distribution, which is instead stable.
We have found that the high-latitude gamma-ray sky above 30 • is composed of (25 ± 2)% point sources, (69.3±0.7)% Galactic foreground, and (6±2)% isotropic diffuse background emission. Both the integral pointsource component and the sum of the Galactic foreground and diffuse isotropic background components were stable against Galactic latitude cuts and changes of the Galactic foreground modeling. The choice of the Galactic foreground can, however, affect the integral value of the diffuse isotropic background component itself.
With respect to the recent IGRB measurement by Ackermann et al. (2015b), this analysis allowed us to clarify between 42% and 56% of its origin between 1 GeV and 10 GeV by attributing it to unresolved point sources.
Using the definition of µ in the limit ∆S → dS and rewriting in terms of x m as defined in Equation (4) eventually gives the representation of the generating function quoted in Equation (3), i.e., This section describes our approach of deriving the source-count distribution dN/dS (uncorrected for detection efficiency) from the 3FGL catalog. The dN/dS distribution was derived self-consistently for each ROI considered in the article. We first selected all 3FGL sources contained in a given ROI. For each source we adopted the best-fit spectral model (power law, log-parabola, power law with exponential or super-exponential cutoff) indicated in the catalog, using the reported best-fit parameters. The source photon flux in the energy range of interest was calculated by integrating this spectrum. The dN/dS was built as a histogram from the above-mentioned flux collection, using appropriate binning and normalizing it to the solid angle covered by the ROI.

C. NODE-BASED APPROACH
The node-based approach as introduced in Section 4.3.2 serves as an independent cross-check for the complementary approach of keeping the positions of breaks as free fit parameters. We applied the node-based approach to the |b| ≥ 30 • data between 1 GeV and 10 GeV. The choice of the node positions was driven by two criteria, i.e., (a) to reasonably approximate the bright-source and intermediate regions covered by catalog data, and (b) to approximate possible features in the faint-source region without overfitting the data. We therefore chose seven nodes: 5 × 10 −7 , 10 −8 , 10 −9 , 3 × 10 −10 , 3 × 10 −11 , 10 −11 , and 5 × 10 −12 cm −2 s −1 . Remaining parameters and priors were chosen in the same way as discussed in Section 4.3.3 for the hybrid approach.
The dN/dS fit employing the node-based approach is shown in Figure 18. The fit matches well the results found in Section 5 within statistical uncertainties.

D. MONTE CARLO SIMULATIONS
The analysis method and the techniques of fitting the pixel-count distribution were validated with Monte Carlo simulations. We used the gtobssim utility of the Fermi Science Tools package to simulate realistic mock maps including a point-source contribution, the Galactic foreground, and a diffuse isotropic background component. Mock maps were analyzed with the same analysis chain as used for the real data. D.1. Setup Mock data were simulated for a time period of 5 years, using P7REP instrumental response functions and the Fermi -LAT spacecraft file corresponding to the real data set. Data selection resembled the procedure applied for real data. Accordingly, an energy range between 1 GeV and 10 GeV was chosen, and the effective PSF was derived in compliance with the simulated data set.
To demonstrate the applicability of the analysis and to investigate the sensitivity, we simulated realizations of four different toy source-count distributions, tagged A1, A2, B, and C. In all four cases, dN/dS was modeled with a broken power law, where n 1 denotes the index above the break and n 2 the index below the break: (A1) no break, with n 1 ≡ n 2 = 2.0, (A2) break at 10 −10 cm −2 s −1 , with n 1 = 2.0, n 2 = 1.6, (B) break at 10 −10 cm −2 s −1 , with n 1 = 2.3, n 2 = 1.6, and (C) break at 10 −10 cm −2 s −1 , with n 1 = 1.6, n 2 = 2.5. In particular, model A1 approximates what was found in the real data (see Section 5). Model A2 was chosen to investigate the sensitivity of the analysis in the faint-source region, while models B and C impose two extreme scenarios.
Point-source fluxes were simulated according to the given dN/dS model, and positions were distributed isotropically across the sky. Realized sources were passed to gtobssim individually. The flux range covered by the dN/dS distributions was limited to the interval [10 −12 , 10 −8 ] cm −2 s −1 . The lower bound of this interval was chosen to be sufficiently small to investigate the sensitivity limit. At the same time, the upper bound ensures a setup that is reasonably simple to study, while resembling the real data in all flux regions except the bright-source region. Flux spectra of individual point sources were modeled with power laws with a fixed power-law index of Γ = 2.0. In addition, models A1 and A2 were simulated incorporating a distribution of point-source spectral indices. We assumed a Gaussian distribution centered on Γ = 2.4, with a half-width σ Γ = 0.2.
The Galactic foreground was modeled using the template discussed in Section 2.2.3. The isotropic background emission was modeled with respect to the analysis cuts. The model is given by the corresponding analysis template iso clean front v05.txt 26 . The simulated background emission between 1 GeV and 10 GeV was normalized to an integral flux of ∼ 3 × 10 −7 cm −2 s −1 sr −1 in the case of the fixed-index simulations and to ∼ 1.5 × 10 −7 cm −2 s −1 sr −1 otherwise. To investigate a possible bias caused by a distribution of spectral indices, model A2 was simulated without any backgrounds (source-only), increasing sensitivity.

D.2. Results
The mock data were analyzed applying the procedure established in Section 5. The MBPL approach was conducted allowing three free breaks. Priors were adjusted to cover the intermediate and faint-source regions appropriately. The hybrid approach was carried out choosing two free breaks and a node at 5 × 10 −12 cm −2 s −1 (10 −12 cm −2 s −1 ) in the case of simulations with a fixed (variable) point-source spectral index. The node was placed at the faint cutoff deduced from the MBPL fit.
The results of the analyses are depicted in Figure 19 for the fixed-index simulations and in Figure 20 for the simulations including the spectral-index distribution. Figures 19(a) and 20(a) demonstrate that the MBPL approach recovered well the simulated dN/dS distributions (red data points) in the intermediate and faint-source regions. It can also be seen that the dN/dS fit follows statistical fluctuations around the model within allowed degrees of freedom. The position of the break, corresponding to parameter S b2 of the model fit, is well constrained and in good agreement with the simulated input. However, uncertainty bands are biased for very faint sources; in particular, for model C a sensitivity cutoff before the faint end of the simulated source distribution was found. The mismatch increases for the results obtained from the Bayesian posterior, while the profile likelihood fit is comparably more accurate. This behavior becomes most pronounced for model C.
The bias of the fit in the faint-source region can be significantly reduced with the hybrid approach; see Figures 19(b) and 20(b). The hybrid approach resolved the sampling issues affecting the Bayesian posterior. The data points are well covered by the derived uncertainty bands.
Possible systematics caused by a distribution of point-source spectral indices are addressed by Figure 20. The data sets with dN/dS realizations of models A1 and A2, each simulated incorporating the Gaussian distribution of spectral indices, were analyzed with the same analysis chain as used for real data, i.e., assuming a constant spectral index of 2.4 . Figure 20 shows that no evidence for a systematic effect on the dN/dS fit was found for S S sens . Below the sensitivity limit S sens , the uncertainty bands shift slightly downward in comparison to model A1 in Figure 19. The high statistics of the source-only simulation of model A2 indeed increased the sensitivity (see bottom row of Figure 20), as expected. We found that the break was recovered well, again indicating no important systematic effect.
The Galactic foreground normalization parameter A gal was found to be ∼1.05 in all considered scenarios, with no evidence for a dependence on the Galactic latitude cut. For the realization of model A1 for fixed spectral indices, for instance, the value of A gal obtained from the posterior was 1.050 ± 0.002, 1.055 ± 0.005, and 1.066 ± 0.014 for Galactic latitude cuts of 10 • , 30 • , and 50 • , respectively. Profile likelihood parameter estimates were similar, with slightly larger uncertainties. The overall effect of obtaining A gal larger than 1 can be attributed to remaining degeneracies between the Galactic foreground model and the diffuse isotropic background component. However, a slight dependence on the Galactic latitude cut cannot be excluded within statistical uncertainties.
In conclusion, all toy distributions were well reproduced with the hybrid approach within statistical uncertainties. The mock data indicate that the actual sensitivity depends on the source-count distribution and the background components, matching our expectation (see Section 4). One can nevertheless conclude from the two extreme scenarios (models B and C) that the sensitivity estimate S sens constitutes a conservative benchmark for the energy band between 1 GeV and 10 GeV.  Figure 19, but for two realizations of simulations including a point-source spectral index distribution (see text for details). The analysis has been carried out assuming a fixed spectral index of 2.4 .