{'fraction_non_alphanumeric': 0.03714447257798864, 'fraction_numerical': 0.016494997418815605, 'mean_word_length': 4.378106821787414, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Testing gamma-ray models of blazars in the extragalactic sky

The global contribution of unresolved gamma-ray point sources to the extragalactic gamma-ray background has been recently measured down to gamma-ray fluxes lower than those reached with standard source detection techniques, and by employing the statistical properties of the observed gamma-ray counts. We investigate and exploit the complementarity of the information brought by the one-point statistics of photon counts (using more than 10 years of Fermi-LAT data) and by the recent measurement of the angular power spectrum of the unresolved gamma-ray background (based on 8 years of Fermi-LAT data). We determine, under the assumption that the source-count distribution of the brightest unresolved objects is dominated by blazars, their gamma-ray luminosity function and spectral energy distribution down to fluxes almost two orders of magnitude smaller than the threshold for detecting resolved sources. The different approaches provide consistent predictions for the gamma-ray luminosity function of blazars, and they show a significant complementarity.


I. INTRODUCTION
Since the start of operations of the Large Area Telescope (LAT) on board the Fermi satellite, the gamma-ray sky has become an extremely powerful tool to test the nature of high-energy emissions at all latitudes. The all-sky gamma-ray emission measured by Fermi -LAT is typically described in terms of: (i) Galactic and extragalactic resolved point (and few extended) single sources [1]; (ii) a Galactic diffuse emission, caused by the interaction of cosmic rays with the interstellar gas and radiation fields [2]; (iii) the unresolved gamma-ray background (UGRB), which is what remains of the total measured gamma-ray emission after the subtraction of (i) and (ii) (sources that are too faint to be detected individually are defined as unresolved ). The Extragalactic Gamma-ray Background (EGB) instead includes all the sources of gamma rays outside the Galaxy which have been resolved, plus the UGRB.
The UGRB is statistically isotropic, with tiny anisotropies that have been detected at small angular scales [3][4][5]. The anisotropies in the UGRB can be ascribed to the global contribution of one (or more) unresolved populations of point sources. The UGRB could also hide signals of annihilation or decay of dark matter particles in our Galactic halo, or in outer galaxies [6][7][8]: however, these searches are hampered by significant uncertainties [9], among which the ones connected to the contribution from astrophysical unresolved point sources playing a major role.
Several source populations contribute to the EGB. At high latitudes, Fermi -LAT has detected blazars, radio galaxies, star forming galaxies (SFGs) and milli-second pulsars [1,10]. Blazars, a class of active galactic nuclei (AGN), are the most numerous population of individual extragalactic gamma-ray sources [1,[11][12][13][14][15]. Depending on the orientation of the relativistic jet of the active galaxy with respect to the observer, AGNs are divided in blazars and misaligned AGNs (mAGNs) [16]. Blazars are in turn sub-divided into two categories, depending on the radio luminosity and the morphology of the emission: BL Lacs have low radio luminosity, which comes predominantly from the center and the jets, while flat spectrum radio quasars (FSRQs) have high radio luminosities which are concentrated in the edge-elongated radio lobes [17]. They also have different photon indices, softer for the FSRQs (∼ 2.4) and harder (∼ 2.1) for the BL Lacs [18]. In addition to AGN emission, the mechanism causing the diffuse emission in the Milky Way, such as the interaction of cosmic rays in the interstellar gas and with interstellar radiation fields, is expected to produce gamma rays in SFGs. Only few galaxies of this type have been detected so far in gamma rays, e.g. the M31 and M33 [19]. They are intrinsically faint but numerous [20], and the extrapolation of models suggests that they can possibly contribute to the observed UGRB in a significant way [8,15,18,[20][21][22][23][24]. However, the extrapolation of the gamma-ray source count to the unresolved flux regime is guided by correlations to the source count observed at different wavelengths, and consequently suffers from significant uncertainties [25].
The contribution to the EGB from these gamma-ray source populations can be characterized by their differential source count distribution dN/dS. This is the areal source density per solid angle element 1 , where N is the number of sources in a given flux interval (S, S + dS), and S is the integral gamma-ray flux of a source in an energy bin. The dN/dS for each source class, in the resolved regime, can be determined through the cataloged point sources. The number of resolved sources is however limited by the detection threshold of the instrument, which needs to be estimated for each catalog [26]. The dissection of the EGB composition is currently complemented by statistical methods, able to dig deeper into the unresolved regime. In fact, analyses employing the statistical properties of the observed gamma-ray counts map have recently measured the contribution from individual sources, and the diffuse EGB components, down to lower gamma-ray fluxes with respect to standard sourcedetection methods [27][28][29][30][31][32]. In particular, in Refs. [30] and [33] it was shown that the 1-point probability distribution function (1pPDF) of counts maps serves as a unique tool for precise measurements of the contribution from unresolved point sources to the gamma-ray sky and the EGB's composition. Within the 1pPDF analysis, the contribution from unresolved point sources to the EGB has been characterized by fitting the non-Poissonian contribution of sources to the photon counts, and with the prediction computed from a description of the dN/dS with a generic multiply broken power law (MBPL). As shown in [33], the 1pPDF analysis, performed with the generic MBPL approach, has the sensitivity to probe the extrapolation of the dN/dS blazar models in the unresolved flux regime.
The 1pPDF method can be generalized to include a more physical parametrization of the dN/dS. We perform here, for the first time, a fit of Fermi -LAT data at latitudes |b| > 30 deg with the 1pPDF method using a specific phenomenological model for describing the gamma-ray emission from the blazar population as the dominant contributor. In combination with this analysis, we consider the two-point angular power spectrum of the UGRB, recently measured on 8 years of Fermi-LAT data [5]. Also in this case, blazars are expected to dominate the anisotropy signal [34], and it has been shown that the gamma-ray angular power spectrum (APS) has the power to constrain the modeling of the unresolved blazar component [35].
In summary, in this paper we combine the investigation of the blazar component in the gamma-ray extragalactic emission, showing that the 1pPDF and two-point APS offer complementary information in the determination of the parameter of blazar models. We then confront these results with the characterization of blazar features we obtain in the resolved regime by using the most recent catalogs of Refs. [1,36].
The paper is organized as follows. Section II describes the model we adopt for the gamma-ray luminosity function and spectral energy distribution of blazars. The computation of the relevant observables is outlined in Section III. Results on the combined analysis of the APS and 1pPDF are presented in Section IV. In Section V we summarize our results and main conclusions.

II. MODEL FOR THE BLAZAR POPULATIONS
The main aim of the paper is to constrain the model of the gamma-ray emission of blazars at all redshift by applying the 1pPDF and APS analyses. These two observables can be computed from the gamma-ray luminosity function (GLF) and spectral energy distribution (SED) of blazars. We consider here the model for the GLF and SED derived in Ref. [24]. The authors of Ref. [24] do not differentiate the two blazar classes (BL Lacs and FS-RQs) since the adoption of a larger sample allows for a better determination of the integrated emission from the whole population in the regime of overlapping luminosities. Specifically, we adopt the following decomposition of the GLF Φ(L γ , z, Γ) = dN/dL γ dV dΓ (defined as the number of sources per unit of luminosity L γ , co-moving volume V at resdhift z and photon spectral index Γ) in terms of its expression at z = 0 and a redshift-evolution function: where L γ is the rest-frame luminosity in the energy range 0.1 − 100 GeV, given by L γ = 100 GeV 0.1 GeV dE r L(E r ), with: E being the observed energy, related to the rest-frame energy E r as E r = (1 + z) E. The co-moving volume element in a flat homogeneous Universe is given by d 2 V /dΩdz = c χ 2 (z)/H(z), where χ is the co-moving distance (related to the luminosity distance d L by χ = d L /(1 + z)), and H is the Hubble rate. At redshift z = 0, the parametrization of the GLF model is [24]: where A is a normalization factor, the indexes γ 1 and γ 2 govern the evolution of the GLF with the luminosity L γ and the Gaussian term takes into account the distribution of the photon indices Γ around their mean µ(L γ ), with a dispersion σ. The mean spectral index is allowed to slightly evolve with the luminosity from a value µ * as [24]: Following the results obtained in Ref. [24], we adopt the luminosity-dependent density evolution (LDDE): Concerning the SED, we model it through a double power law: where we assume that E b correlates with Γ according to log(E b /(GeV) = 9.25 − 4.11 · Γ. Given a SED, the flux S(E min , E max ) in a given energy interval is obtained as: where τ (E, z) describes 2 the attenuation by the extragalactic background light (EBL) [37]. The energy flux S E (E 1 , E 2 ) in a given energy interval is instead: The GLF and SED models have a large number of free parameters, which in Ref. [24] have been determined by fitting the Fermi -LAT catalog data, and follow-up observations of blazars. In our analysis we will adopt as free parameters those which grab the dominant dependencies, i.e. the GLF normalization parameter A, the central value µ * for the photon spectral index Γ, the powerlaw index γ 1 that governs the dependence of the GLF at high luminosity and the central values of the powerlaw indices p * 1 and p * 2 which set the redshift dependence of the LDDE. All other parameters have been fixed at the values obtained in Ref. [24], for definiteness. We have checked both larger (including e.g. also z * c ) and smaller sets of free parameters, obtaining that our method is sensitive dominantly to the stated parameters and we will therefore report the results on this set.

III. THE TECHNIQUES FOR DISSECTING THE BLAZAR MODELS
As mentioned above, in this paper we analyze the 1pPDF, APS and the most recent gamma-ray catalogs and their combined constraining power. In this section, we describe each of these techniques.
A. The 1pPDF photon-count statistics technique The 1pPDF method relies on defining a probability generating function -generically derived from a superposition of Poisson processes -for the photon count maps. The mathematical formulation of the 1pPDF method, its implementation, and its application to Fermi -LAT data are discussed in [28,30,33], to which we refer for any detail (see also [38]). In this method, the expected number of point sources in map pixel p contributing exactly m photons to the total pixel photon content is given by the dN/dS, being S the integral photon flux of a source in the energy band [E min , E max ] (observed energies) as defined in Eq. (10): where Ω pix is the solid angle of the pixel, and C (p) (S) denotes the average number of photons by a source with flux S which contribute to the pixel p. The isotropic distribution of gamma-ray point sources dN/dS was generically parameterized with a MBPL in Refs. [30,33], with the overall normalization, a number of N b break positions and therefore N b + 1 power-law components connecting the breaks as free parameters. In the current analysis we instead progress beyond this generic description, and assess if the dN/dS of high latitude, extra-galactic sources required to fit the data can be described by a blazar population, described by the physical model of the previous section. We concentrate our analysis to photon energies in the interval from 1 GeV to 10 GeV, which is where we have at the same time large statistics and a good angular resolution of the Fermi-LAT detector. The differential number of blazars per integrated flux and solid angle dN/dS can be computed from the model described in Section II as: where L γ (S E , z, Γ) is the luminosity of a source endowed with an energy flux S E , located at redshift z and with spectral index Γ, being S E the flux in a specific energy bin. The integration bounds for Γ in Eq. (13) are such to properly cover the distribution of observed blazars, while the integration bounds for the redshift z cover the interval in which we expect the vast majority of their emission [24].
Within the 1pPDF method applied here, the total gamma-ray emission is described by summing an isotropic distribution of point-like blazars and two diffuse background components, the Galactic foreground emission and an additional isotropic component, which are described by 1-photon source terms. The total diffuse contribution x (p) diff is then given by: For the isotropic component x (p) iso , we use the integral flux F iso as a sampling parameter, in order to have physical units of flux 3 . The first term accounts for the Galactic foreground emission, described with an interstellar emission model (IEM). Further details on the considered IEMs are given below. The global normalization of the IEM template A gal is taken as a free fit nuisance parameter. The second term describes all contributions indistinguishable from purely diffuse isotropic emission. The diffuse isotropic background emission is assumed to follow a power law spectrum (photon index Γ = 2.3), with its integral flux F iso serving as the free normalization parameter.
Concerning the data-set, we analyzed all-sky Fermi -LAT gamma-ray data from 2008 August 4 (239,557,417 s MET) through 2018 December 10 (566,097,546 s MET). We used Pass 8 data 4 , along with the corresponding instrument response functions. The Fermi Science Tools (v10r0p5) 5 were employed for event selection and data processing. The data selection referred to standard quality selection criteria (DATA_QUAL==1 and LAT_CONFIG==1), to values of the rocking angle of the satellite smaller than 52 • , and maximum zenith angle of 90 • . We selected events passing the ULTRACLEANVETO event class, and we use the corresponding instrument response functions. A correction for the source-smearing effects coming from the finite detector point-spread function (PSF) has been also applied in Eq. (12), as detailed in Ref. [30]. To avoid significant PSF smoothing, the event sample is restricted to the PSF3 quartile (see [30,33]). Data are analyzed in the energy range from 1 GeV to 10 GeV, and binned using the HEALPix equalarea pixelization scheme [39] with a resolution parameter κ = 7, being N pix = 12N 2 side the number of pixels, with N side = 2 κ . To avoid significant contamination from the diffuse emission of our Galaxy, we analyzed the data for |b| > 30 deg. The 1pPDF likelihood function is defined as the L2 method in [30]. The nested sampling algorithm included in the MultiNest framework [40] is used to sample the parameter space, with 1500 live points together with a tolerance criterion of 0.2. The IEM has been fixed according to the official spatial and spectral template as provided by the Fermi -LAT Collaboration for the Pass 8 analysis framework (gll_iem_v06.fits, see Ref. [41]).

B. The angular power spectrum technique
The APS of the gamma-ray intensity fluctuations is defined as: where the indexes i 3 We note that that the ratio x (p) iso /F iso does not depend on F iso . 4 Publicly available at https://heasarc.gsfc.nasa.gov/ \FTP/fermi/data/lat/weekly/photon/.
More details are found at https://fermi.gsfc.nasa.gov/ssc/data/analysis/ documentation/Cicerone/Cicerone_Data/LAT_DP.html 5 https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/ and j label here the energy bins. The coefficients a m are the amplitudes of the expansion in spherical harmonics of the intensity fluctuations, and n identifies the direction in the sky. The sum defines an average over the modes m for each multipole . For i = j the APS describes the energy auto-correlation, while for i = j the APS describes the cross-correlation of the fluctuations in two different energy bins.
If the population that dominates the APS is composed of point-like, relatively bright and non-numerous sources, its anisotropy signal is dominated by the so-called Poisson noise term and the APS does not depend on the angular multipole , i.e. C ij C ij P . One can check that, at the level of fluxes probed by the Fermi -LAT telescope, this is the case for blazars [34]. In our physical model, the blazar APS can be computed as: The upper and lower bounds in the L γ integration are set to L min = 10 43 erg/s and L max = 10 52 erg/s [24]. The term Ω(S, Γ) accounts for the Fermi -LAT sensitivity to resolve a source, which depends on its photon flux S and spectral index Γ. Our reference model for Ω is a step function: sources with a given spectral index Γ are considered to be unresolved if their flux in the energy range 1-100 GeV (10-1000 GeV) is below the detection threshold S thr (Γ) of the fourth Fermi -LAT catalog (4FGL) of gamma-ray sources [1] (the third catalog of hard Fermi -LAT sources (3FHL) [42]) catalog. The 4FGL catalog is based on eight years of data in the energy range from 50 MeV to 1 TeV and contains 5065 sources which are detected with a confidence level (C.L) above 4σ. On the other hand the 3FHL catalog is focussed on energies above 10 GeV. It is based on 7 years of data and contains 1556 objects. The threshold S thr is defined such that > 98% of sources in the catalog with spectral index Γ have a flux larger than S thr (Γ). We checked (as discussed in the following) that the exact procedure adopted for implementing the Fermi -LAT sensitivity threshold does not alter our final results 6 . The computation of the C ij P requires the same ingredients as in the dN/dS case: the GLF and SED. One can interpret the C P as the second moment of the dN/dS, as can be seen by comparing Eqs. (13) and (15). This allows us to combine the constraining power of the 1pPDF method and the anisotropy analysis in the determination of the free parameters characterizing the blazar model.
The measured C P 's adopted in our analysis are taken from Ref. [5], where the measurement is performed on Pass 8 data 4 of the P8R3_SOURCEVETO_V2 event class with PSF1+2+3 type events. The data selection comprises 8 years, binned in 12 energy bins between 524 MeV and 1 TeV. The contribution from the resolved sources in the energy range (0.5 − 14.5) GeV, (14.5 − 120) GeV, (120 − 1000) GeV is masked using the source list of the FL8Y, FL8Y+3FHL, 3FHL catalogs, respectively. The low latitude Galactic interstellar emission is masked, and a Galactic diffuse template based on the model gll_iem_v6.fits [41] has been subtracted in order to reduce the contamination from high-latitude Galactic contribution. For a full description of methods and results, we refer the reader to Ref. [5].
The fit of the APS is performed on the auto-and crosscorrelation energy bins. The χ 2 APS is defined as: Here the subscripts meas denotes the measured C P from Ref. [5] and the subscript th denotes the C P calculated from Eq. (15). Furthermore, σ 2 C ij P is the uncertainty of the measured C P . The likelihood L = exp(−χ 2 APS /2) is sampled using the MultiNest package in a configuration with 2000 live points, an enlargement factor of efr=0.7, and a stopping parameter of tol=0.1. The results in the next section will be discussed within the frequentist framework.

C. The 4FGL and 4LAC catalogs for resolved blazars
As a further technique, relevant for resolved sources, we analyze the most recent source catalogs to constrain the blazar model [24]. The 4FGL catalog [1] is now available, as well as an early release of the fourth catalog of AGNs (4LAC) [36], both obtained with eight years of data. In addition to the blazar type classification, the 4LAC collects also the spectral features, variability and redshift estimates, the last being crucial to constrain the GLF. The constraints to the GLF obtained from the catalogs of resolved blazars will also be used in Sec IV C as a prior for the APS fit.
We use the source count distributions extracted from the 4FGL catalog in a χ 2 -fit in which we vary the same five GLF parameters as for the 1pPDF and APS fits: A, µ * , γ 1 , p * 1 , and p * 2 . The total χ 2 4FGL can be divided into three contributions arising from the total number of observed point sources, the number of associated blazars 7 , 7 In this paper associated blazars refers to the sum of identified and associated sources classified as BL Lacs, FSRQs, or blazars of uncertain type (BCU), namely, the 4FGL source classes are BLL, BCU, FSRQ, bll, bcu, fsrq. and blazars with redshift measurements: In the following we will define each contribution. For the first term, we extract the source-count distribution of all sources, (dN/dS) all,i , in 12 flux bins i ranging equally spaced in log(S) from 10 −12 cm −2 s −1 to 10 −7 cm −2 s −1 , where S = S(1 GeV, 100 GeV). To avoid a strong contamination of Galactic sources, we restrict the analysis to sources at latitudes with |b| > 30 deg. We compare the extracted source count distribution to the average source count distribution from the blazar GLF dN/dS th,i , which is the integral of dN/dS (see Eq. (13)) in the energy bin i divided by ∆S i . Among the unassociated sources in the 4FGL catalog, we expect that some of them are not blazars. Therefore, we use (dN/dS) all,i only as an upper limit in the fit. In terms of the χ 2 definition this means: The contribution of the associated sources is defined with a very similar procedure. There are only two small differences: (i) instead of extracting the total source count distribution, we extract the source count distribution of associated blazars, (dN/dS) as,i , and (ii) we use (dN/dS) as,i as a lower limit in the fit since the association in the catalog might be incomplete. The χ 2 as is defined by: Note that by taking χ 2 all as an upper limit and χ 2 as as a lower limit, there is no double counting in Eq. (17).
In the last term of Eq. (17) The source count distribution, (dN/dS) z,ij , is extracted equivalently to the procedure described above. The only difference is that the number count is restricted to the redshift in each bin. The corresponding source count distribution of the GLF model, dN/dS th,ij , is obtained by restricting the integration range of z in Eq. (13) to the redshift bin. Since the redshift measurements in the catalog are incomplete, the source count distributions extracted from the 4LAC catalog are taken as lower limits: Let us observe that the combination of χ 2 all + χ 2 z allows the determination of p * 1 and p * 2 , while the combination of χ 2 all + χ 2 as constrains A, γ 1 and, mildly, µ * such that effectively we avoid any double counting.
Finally, we have to discuss one subtlety which is connected to the detection efficiency for point sources. As discussed above, we model the sensitivity to detect point sources in the 4FGL catalog by a step function with a threshold flux, S thr (Γ). Below the flux threshold the observed source count distribution in the catalog drops and can no longer be used in the fit. Consequently, we restrict the sum over i in Eq. (18) to those flux bins which are above the maximal threshold flux at Γ ∼ 2.3 of S thr = 1.1 × 10 −10 cm −2 s −1 . Note that we do not require to restrict the sums in Eqs. (19) and (20) because they serve as lower limit and a decrease of the observed source count distribution only weakens the limit.
In order to sample the 5-dimensional parameter space we use the MultiNest package. We adopt 2000 live points, an enlargement factor of efr=0.7, and a stopping parameter of tol=0.1. The results presented in the next section are interpreted in the frequentist approach.

IV. RESULTS
In this section, we first present the results obtained by applying the 1pPDF analysis to the specific blazar dN/dS model introduced in Section II. Then, we probe the blazar model through the APS analysis, and combine the two methods. Finally, we check the compatibility of our results with the 4FGL catalog.

A. Results from the photon-count statistics analysis
The results on the determination of the dN/dS for high latitude blazars, obtained with the 1pPDF analysis, is shown in Fig. 1: the red solid line is the result obtained by using the blazar model of Section II, while the dashed gray line refers to the results obtained by employing a MBPL, as done in [30]. The shaded areas of corresponding color denote the 1σ frequentist uncertainty. For the physical blazar model of Section II we vary the parameters A, µ * , γ 1 , p * 1 and p * 2 and marginalize over two nuisance parameters, the normalization of the Galactic foreground emission A gal , and the flux of the isotropic gamma-ray emission, F iso . In the case of the MBPL, we adopt a mode with three nodes (see [30] for details) and we obtain the following results: for the normalization parameter A S = 2.31 +7.67  sponding index n 4 at very low fluxes, is not statistically significant. Finally, Fig. 1 also shows the counts for the resolved sources listed in the 4FGL catalog. The MBPL result shows that the 1pPDF is able to determine the behavior of the dN/dS more than one order of magnitude in flux lower than the catalog threshold (S ∼ 2 − 3 × 10 −10 cm −2 s −1 ), namely at S ∼ 8 × 10 −12 cm −2 s −1 , below which the uncertainty band increases significantly. When this is translated to the physical blazar model, it allows to determine and trust the behavior of the dN/dS down to the same flux level, therefore extending the understanding of the blazar model in the unresolved regime. Let us also notice that the fact that the results obtained with the physical blazar model are very well consistent with those obtained with the generic MBPL analysis and with the 4FGL catalog sources, reinforces our assumption that point sources emitting photons at high latitudes in the energy range from 1 GeV to 10 GeV are consistent with a blazar origin even in the unresolved regime.
The best-fit values of the relevant parameters of the GLF blazar model, together with their uncertainties, are reported in Tab. I. We obtain values which are largely compatible (except for p * 2 , where compatibility is present only at about the 2σ level) with the reference model of Ref. [24], which was adapted to the resolved component and to a source catalog predating the 4FGL. In Tab. I we also show the results for the same parameters, obtained by fitting the 4FGL catalog (see Sec. IV C and Fig. 4), in which case the agreement between our results and Ref. [24] is well inside 1σ for all parameters. These results indicate that the unresolved blazar component (down to fluxes of the order of about 8 × 10 −12 cm −2 s −1 ) has similar properties as those which are currently resolved, I: Best-fit parameters for each of the techniques investigated in this paper. The first column lists the free parameters, while the following four columns contain the corresponding best fits. The last column reports the reference values from Ref. [24].

B. Results from the angular correlation analysis
As for the 1pPDF analysis, we compute the autocorrelation APS in the [1,10] GeV energy bin and the cross-correlation APS of this bin with the other twelve bins from 0.5 GeV to 1 TeV adopted in Ref. [5]. In this case, in addition to the A, µ * , γ 1 , p * 1 and p * 2 parameter, we have nuisance parameters which allow us to change the flux threshold of the point-source detection by a factor of k = 0.5 to 2.0 relative to S thr (more comments are provided at the end of this subsection).
The results are reported in Fig. 2. The left-panel refers to the auto-correlation APS amplitude C P as a function of the energy, while the right panel stands for one case of cross-correlation, specifically the cross-correlation of the [8.3, 14.5] GeV energy bin with all the other bins. We note that the best-fit model well reproduces the measurement obtained in Ref. [5], demonstrating that the blazar model is compatible also with the APS of the photon field fluctuations, and that the study of the unresolved components by means of two different methods (the 1pPDF and the APS) provide consistent results, as quantified below. The best-fit values for the parameters and their errors are reported in Tab. I: the results are well compatible with those obtained in the 1pPDF analysis, including the value obtained for the p * 2 parameter. While the 1pPDF 8 There is a subtlety connected to the counting of the degrees of freedom (dof) in the fit of the blazar model to the 4FGL+4LAC catalog data. We use the total number of point sources as upper limit and the number of BLL+BCU+FSRQ as lower limit in 7 flux bins. The number of fit parameters is 5. Using this information gives a dof of 2. The subtlety is that, on top of the mentioned constraints, we use redshift information as lower limit in the fit, in effectively 28 bins. However, the χ 2 at the best-fit point is only marginally affected by these lower limits. So, we decided not to count this information in the dof stated in the table.
and APS results are well compatible with the catalog results, the fact that p * 2 turns out somehow lower for both analyses (sensitive to the unresolved blazar component) might be indicative that the fainter blazar emission starts to point toward a slightly different regime.
Previous analyses of gamma-ray APS found evidence for two populations instead of a single population [4,5,35]. We also test here this hypothesis, following a strategy already used in Ref. [35]. On top of the blazar physical model, we add an additional soft and faint component for which we assume dN/dS = A PWL (S/S 0 ) −βPWL (where S refers to the flux in the energy bin 1-100 GeV) and an energy spectrum given by dN/dE ∼ E −ΓPWL . We then perform a fit with the sum of the blazar physical model plus such additional generic power-law component. In total, this fit involves 8 free parameters: the 5 parameters already used in our reference analysis, plus A PWL , β PWL , Γ PWL . We find a slight improvement in the χ 2 , but not statistically significant, being smaller than at the 2σ C.L. This then justifies the adopted procedure to fit the C P with a single blazar population: namely, the underlying assumption of our analysis that blazars are the dominant contributor to the unresolved gamma-ray sky, in the regime just below the Fermi-LAT detection threshold. Notice that we are adopting a different approach as compared to Ref. [5], where a preference for 2 populations was instead present: we describe the gamma ray emission in terms of a physical blazar model and we allow for a distribution of their spectral indices Γ with a dispersion of σ = 0.28 [24] (see Eq. (3)), instead of adopting a given spectral index as done in Ref. [5]. In this case, the single-blazar model is able to describe the measured APS. We leave for a future work the investigation of the possible presence of subdominant additional unresolved populations.
Let us conclude this section by elaborating more on the issue of the flux threshold sensitivity. The measurement of angular correlations is performed masking sources from the FL8Y and 3FHL catalogs. Therefore, the measured C P depends on the efficiency of the Fermi-LAT source detection, see Eq. (15). An exact estimate of such efficiency is challenging, and a typical assumption when    calculating the C P in the blazar model is that this efficiency Ω Θ abruptly changes from 0 to 1 at a certain flux denoted as S thr . We adopt a Θ-like cut as benchmark. Note that in contrast to many previous analysis we take also the Γ dependence of S thr into account. In order to test the impact on the blazar fit of our modeling of the efficiency, we replace S thr by k S thr , and marginalize over k. Furthermore, as a test we replace the Θ-like cut by a smooth function: The parameter η is varied from 2.5 to 4. We have verified that the results of the physical parameters (A, γ 1 , p * 1 , p * 2 , µ * ) are stable against these changes of the functional form of the efficiency function, with a value for the nuisance parameter k close to 1.

C. Complementarity of 1pPDF, C P and 4FGL catalog
The two methods adopted to investigate the unresolved side of the gamma-ray emission (1pPDF and APS) produce compatible results, but also provide complementary information. This can be seen by analyzing the full parameter space, reproduced in Fig. 3, which shows the 1-dimensional and 2-dimensional χ 2 distributions.
The preferred regions obtained with the two techniques always exhibit overlap within a 2σ C.L, demonstrating compatibility. However, the APS analysis significantly constrains the central value of the blazar spectral index µ * , while being much less effective on the other parameters. This occurs because the APS analysis involves several energy bins (through the cross-correlation in energy) and this allows us to characterize the blazar SED. On the other hand, the 1pPDF method is more constraining on the other GLF parameters, especially the normalization A and the parameter γ 1 which governs the luminosity evolution. Clearly, since in the 1pPDF we are adopting a single energy bin, we have small sensitivity on the SED.
The results of the blazar model fit to the 4FGL catalog are shown in Fig. 4. The lower and upper black triangles mark the source count distribution of all point sources (dN/dS) all,i and point sources associated as blazars (dN/dS) as,i , respectively. The best fit of the blazar model lies between the two source count distributions, which serve as upper and lower limit in the fit. The colored triangles show the source count distribution in four redshift bins (dN/dS) z,ij . Those data points are a lower limit to the blazar model, since the redshift catalog is incomplete. We observe that the best-fit model fulfills this requirement, and lies above the colored data points for all the redshift ranges. The corresponding best fit parameters for this fit are reported in Tab. I.
The results obtained by fitting the source count distribution of the 4FGL catalog are also provided in Fig. 3 (green contours). The results are well compatible with those obtained with the 1pPDF and APS analyses. As expected, there is very good agreement to the 1pPDF analysis, since the catalogs and the 1pPDF analysis directly probe the number of point sources, although in two different regimes (resolved for catalog, resolved and unresolved for 1pPDF). We note that the catalog fit pro- vides the strongest constraints on the parameter p * 2 , by excluding values smaller than about −7. To interpret this constraint, we remind that p * 2 changes the shape of the LDDE at z > ∼ z * c = 1.25. The other two methods cannot exclude small values of p * 2 since, in contrast to the catalog fit, they do not contain explicit redshift information.
As a further result, we show in Figs. 5 and 6 how the different observables would be reconstructed if only the best-fit from one of the techniques is used. In Fig. 5 we show that the source count distribution provided by the best-fit parameters of the APS analysis is in good agreement with the 1pPDF and 4FGL catalog analyses for what concerns the unresolved regime. On the other hand, the APS study would over-predict the measured dN/dS in the resolved part. The lack of precision of the C P analysis in this regime is somewhat expected, since it is based only on data below Fermi-LAT detection threshold. If one attempts to describe a complete model of both resolved and unresolved blazars, this has to be complemented by other techniques, as we show at the end of this section.
The prediction that would be obtained for the APS as a function of energy by using only the information coming from the 1pPDF or from the 4FGL catalog analyses is  shown in Fig. 6. Since they are obtained in a single energy bin, they cannot be very predictive for what concerns the blazar SED. This becomes manifest if one compares the precision obtained from the APS analysis (blue regions) in the reconstruction of the energy spectrum to what is predicted by the 1pPDF (red) or the 4FGL catalog (green) analyses. Therefore, Figs. 5 and 6 reinstate the complementarity of the different probes in cornering the blazar model. We note that the prediction of the dN/dS from the C P and vice versa show deviations above the 1σ level. A similar deviation is visible also in the parameter contours shown in Fig. 3. We checked explicitly that at the 3σ level all the bands are compatible with the data points. We also checked explicitly the compatibility between the C P and 1pPDF predictions and the dN/dS of the catalog in all our 4 redshift bins.
The contours from 1pPDF cannot be simply combined with APS or 4FGL analyses without computing the appropriate co-variance. Indeed, the 1pPDF uses data both in the resolved and unresolved regimes. The combination can be instead performed between APS and 4FGL analyses, since they rely on separate data-sets. To demonstrate again the complementarity between the C P measurement and the information in the 4FGL catalog, we perform a further joint fit to both observables, in which the joint χ 2 C P +4FGL is defined as sum of the two individual χ 2 s defined in Eqs. (16) and (17), respectively. We obtain a good fit with a minimal joint χ 2 /dof of 94.5/79 which  can be separated into a contribution from the C P fit of 86.2 and the 4FGL fit and 8.3. The combination of both observables guarantees that both, the measured dN/dS (Fig. 5) in the resolved part and the measured C P (Fig. 6) in the unresolved regime, are properly reproduced at 1σ. Furthermore, we observe that the nuisance parameter, k, is very well constrained by the combination of the two methods, since the 4FGL information fixes the dN/dS above S thr .
Results are shown in Fig. 7 and the best-fit values are reported into Tab. I. One can explicitly note the striking complementarity already mentioned above, namely, the best-fit regions shrink to the overlap of the two individual fits. We recommend to use the values of the C P +4FGL fit to obtain a good agreement or the blazar model in the resolved and unresolved regime.

V. CONCLUSIONS
In this paper we adopted and compared different statistical methods to constrain the gamma-ray emission from blazars. Based on the most recent Fermi -LAT data at high Galactic latitudes, we derived the description of the blazar luminosity function and spectral energy distribution, with best-fit parameters provided in Tab  statistical properties of the observed gamma-ray counts. We analyzed the 1pPDF and two-point APS, and compared the results to the characterization provided by the analysis of resolved sources in the 4FGL catalog. We found that the 1pPDF and APS can indeed extend our knowledge of the blazar GLF and SED to the unresolved regime, and are able to determine the dN/dS of blazars down to fluxes almost two orders of magnitude smaller than the Fermi -LAT detection threshold for resolved sources.
The different approaches provide predictions that are generically in good agreement with each other. Moreover, they show a significant complementarity. The APS analysis better characterizes the blazar SED, since it involves several energy bins (and their cross-correlation). The 1pPDF is more constraining for what concerns the normalization and luminosity evolution of the GLF. The analysis of the redshift distribution of the resolved sources in the catalogs allows a more refined determination of the GLF redshift evolution. The complementarity of the different techniques in constraining the parameters of the GLF and SED models of blazars can be appreciated in Figs. 3 and 7.
We plan for future works to further constrain the GLF of blazars, and potentially other source populations, by investigating the 1pPDF in different energy bins, and by performing a two-point correlation analysis with catalogs of blazars at different wavelengths.