Flat spectrum radio quasars and BL Lacs dominate the anisotropy of the unresolved gamma-ray background

We analyze the angular power spectrum (APS) of the unresolved gamma-ray background (UGRB) emission and combine it with the measured properties of the resolved gamma-ray sources of the Fermi-LAT 4FGL catalog. Our goals are to dissect the composition of the gamma-ray sky and to establish the relevance of different classes of source populations of active galactic nuclei in determining the observed size of the UGRB anisotropy, especially at low energies. We find that, under physical assumptions for the spectral energy dispersion, i.e. by using the 4FGL catalog data as a prior, two populations are required to fit APS data, namely flat spectrum radio quasars (FSRQs) at low energies and BL Lacs (BLLs) at higher energies. The inferred luminosity functions agree well with the extrapolation of the FSRQ and BLL ones obtained from the 4FLG catalog. We use these luminosity functions to calculate the UGRB intensity from blazars, finding a contribution of 20% at 1GeV and 30% above 10 GeV. Finally, bounds on an additional gamma-ray emission due to annihilating dark matter are derived.


INTRODUCTION
The extragalactic gamma-ray sky has been surveyed by the Fermi Large Area Telescope (LAT) since the sum-  (Atwood et al. 2009). The outstanding capability of this instrument has been groundbreaking in several aspects of high energy astrophysics. One important result is the detection and cataloging of extragalactic gamma-ray sources. The 8-year Fermi-LAT source catalog, called the 4FGL catalog (Abdollahi et al. 2020) 1 , counts more than 3300 extragalactic sources, which is more than 60% of the entire catalog. Almost the totality of the extragalactic sources are blazars, a sub-class of active galactic nuclei (AGNs) with a jet pointing towards us: 35% are BL Lacs (BLL), about 22% are flat spectrum radio quasar (FSRQs) and about 41% are blazars of unknown type (BCU).
On top of the numerous detected extragalactic sources, even more numerous sub-threshold sources populate the unresolved gamma-ray background (UGRB). 2 The UGRB emission represents about 20% of the total gamma-ray emission and offers a unique observable of the extragalactic gamma-ray sky below the Fermi-LAT source detection threshold.
The UGRB is by definition a mission-time dependent component: The more the Fermi-LAT surveys the sky, the more sensitive it becomes to less bright sources, leaving only the faintest objects unresolved. Guaranteed contributors to the UGRB emission are sub-detectionthreshold blazars (Cuoco et al. 2012;Di Mauro et al. 2018), misaligned AGNs (mAGNs) (Di Mauro et al. 2013), and star-forming galaxies (SFG) (Roth et al. 2021;Tamborra et al. 2014). Additionally, we cannot exclude contributions from more exotic components, such as dark matter (DM) (Ando 2009;Bringmann et al. 2014;Ajello et al. 2015;Fornasa et al. 2016;Zechlin et al. 2018) The UGRB emission has been studied through three main observables: its energy spectrum (Abdo et al. 2011;Ackermann et al. 2015), its 1-point probability distribution function (1pPDF) through the photon counts statistics (Zechlin et al. 2016b,a;Lisanti et al. 2016;Di Mauro et al. 2018), and its angular power spectrum (APS) Fornasa et al. 2016;Ackermann et al. 2018). The latter two observables investigate fluctuations over the UGRB isotropic emission to infer the properties of the underlying sources at the subthreshold level. In this unresolved regime, mAGNs and SFGs, fainter than blazars but extremely more numerous, are expected to dominate the UGRB energy spectrum (Di Mauro et al. 2013;Roth et al. 2021). At the same time, at the current level of sensitivity, the blazars produce a higher level of spatial anisotropy than mAGNs and SFGs, and hence the former are expected to dominate the APS of the UGRB (Di Cuoco et al. 2012). SFGs and mAGNs could eventually emerge once the majority of the blazars have been resolved. Moreover, an improvement of the sensitivity is necessary in order to reveal the large-scale structure (LSS) of the Universe traced by gamma-ray sources (see e.g. Ando (2009)) that is encoded in a multipole-dependent APS. It is, therefore, crucial to update the UGRB anisotropy measurement in parallel to the detection of more sources in the LAT catalogs.
The latest UGRB anisotropy measurement has been performed by the Fermi-LAT Collaboration in 2018 (Ackermann et al. 2018). In that work, 8 years of Pass-8 (R3) data were analyzed, and it was consistently adopted the source catalog based on the same amount of observation time (FL8Y, a preliminary version of the 4FGL). The APS of the UGRB was measured in 12 energy bins between 500 MeV and 1 TeV. Additionally, the cross-correlation signal between the different energy bins (generically denoted by i and j) was derived. In all cases, the APS (above = 50) was compatible with a constant value, C ij P , with no hint of LSS signature in the multipole range considered. This result confirms that the UGRB intensity fluctuation field, at the current level of sensitivity of the detector to point sources is still dominated by a population of relatively bright and not very numerous sources, so that the isotropically distributed fluctuations from Poisson noise dominate over the correlation due to clustering. Additionally, the anisotropy energy spectrum revealed a preference for a double power-law trend (with a high energy exponential cutoff) over a single power law (with a high energy exponential cutoff), placing a spectral break around 5 GeV.
Previous interpretation works, based on antecedent measurements of the UGRB anisotropy energy spectrum, were devoted to determining the components that contribute to the measured signal. In particular, (Ando et al. 2017) studied the results of (Fornasa et al. 2016) and inferred the presence of a second steeper component, in addition to the blazar-only model, emerging below 2 GeV. However, the very soft spectral index implied by this analysis challenges the interpretation in terms of a known source population. Recently, (Manconi et al. 2020) combined the 1pPDF, using methods as in (Zechlin et al. 2016a), and the latest measurement of the anisotropy energy spectrum of the UGRB by (Ackermann et al. 2018) to test blazars models (yet not distinguishing between BLLs and FSRQs) as well as the source count distribution of blazars extracted from the 4FGL catalog. They found that the assumption of the UGRB fluctuation field being entirely dominated by blazars is in agreement with both observables, which appear to show remarkable complementarity. Past works focused also on the DM interpretation of the UGRB anisotropy, as (Fornasa et al. 2016), where numerical simulations were used to model the DM distribution and its uncertainty in order to constrain the contribution from weakly inter-acting massive particles (WIMPs) in Galactic and extragalactic structures. The derived bounds are in the same ball-park as for other UGRB probes, but still significantly above the so-called thermal WIMP scenario. For a comprehensive overview about UGRB-related measurements and interpretation works prior (Fornasa et al. 2016), we address the reader to the review of (Fornasa & Sánchez-Conde 2015).
In this work, we will focus on the latest measurement of the UGRB anisotropy energy spectrum (Ackermann et al. 2018). We investigate the contributions of different blazar types, distinguishing between BLLs and FSRQs. We find that BLLs and FSRQs can account for the totality of the UGRB anisotropy and also well reproduce the spectral features observed by (Ackermann et al. 2018). The analysis allows us to constrain many of the most relevant parameters of the blazar models in the unresolved regime. In a second step, we include the contribution to the UGRB arising from an annihilating DM particle and perform a global analysis to derive constraints on the particle DM parameters. We account for both Galactic and extragalactic DM contributions, under different assumptions of the DM subhalo contribution and by including cross-terms in the anisotropy APS, due to the cross-correlation of the blazars contribution with the DM halos hosting them.
The paper is structured as follows: Section 2 is devoted to blazars and we describe the blazar model adopted in our study, we introduce the fit procedure and we show the results. In Section 3, we discuss the DM constraints for both Galactic and extragalactic DM components. Finally, we conclude in Section 4. Additionally, we present a phenomenological approach to the interpretation of the UGRB anisotropy energy spectrum in Appendix A, while in Appendix C we relate our results to the findings of previous measurements.

MODELING BLAZAR POPULATIONS
In (Manconi et al. 2020) it has been pointed out that a single blazar model is sufficient to describe both the anisotropy level C P and the 4FGL catalog data, at the price of allowing for a relatively broad distribution of the spectral index. Such an approach can be seen as an effective description, where different sub-populations (with narrower spectral index distributions) are combined in a single model . We reproduce the finding of (Manconi et al. 2020), although in a more general way and by using a phenomenological model, in Appendix A. However, we note that the blazar model in (Manconi et al. 2020) was only compared to the catalog data in bins of flux and redshift, but not in bins of spec-tral index which constrains the SED. In contrast, here we intend to use the full catalog information.
In this section, we will therefore consider a physical description of the two populations of blazars that are more numerous in the 4FGL catalog, namely BLLs and FSRQs. We aim to assess their ability to explain the APS measurement. In other words, we test the possibility that FSRQs, with properties compatible with their cataloged sample, are the population that accounts for the low-energy anisotropy found in (Ackermann et al. 2018), while BLLs are at the origin of the high-energy anisotropy. Below we will show that two populations are preferred when the information of 4FGL catalog are included as a prior.

The gamma-ray luminosity function
We summarize here the parametrizations of the gamma-ray luminosity function (GLF) and spectral energy distribution (SED) of BLLs and FSRQs adopted in our analysis. For more details, see  and (Ajello et al. 2014). The GLF Φ(L γ , z, Γ) = d 3 N/dL γ dV dΓ, defined as the number of sources per unit of luminosity L γ , co-moving volume V at redshift z and photon spectral index Γ, is typically decomposed 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, i.e. 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 comoving distance (related to the luminosity distance d L by χ = d L /(1 + z)), and H is the Hubble parameter. We use ΛCDM cosmology with parameters from the final full-mission Planck measurements of the CMB anisotropies (Aghanim et al. 2020). At redshift z = 0, the parametrization of the GLF is: where A is a normalization factor, the indices γ 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 σ. It turns out that the GLF of BLLs has a relatively broad distribution in terms of luminosity. For this reason, we allow the mean spectral index to slightly evolve with luminosity from a value µ * : On the other hand, FSRQs have a narrower distribution, so that the inclusion of this effect would not affect the fit and we can fix µ(L γ ) = µ * . We adopt a luminosity-dependent density evolution (LDDE): . We set δ to 0.64 for both populations , while τ is fixed to 3.16 for FSRQs and to 4.62 for BLLs (Ajello et al. 2014).
The SED is modeled as a power law: For definiteness, the spetral index Γ will be taken to be in the range (1, 3.5) (see also Manconi et al. (2020)). Given the SED, the photon flux S(E min , E max ) in a given energy interval is obtained by: where τ (E, z) describes the attenuation by the extragalactic background light (EBL) (Finke et al. 2010a). Unless explicitly stated otherwise, the flux in the following computations of the dN/dS and the associated figures always refer to the energy bin from 1 GeV to 100 GeV. The free parameters of the model are summarized in Table 2, together with their best-fit values obtained as outlined below.
With the physical models of the GLF and SED at hand, we can compute the differential number of blazars per integrated flux and solid angle as: and the size of the gamma-ray intensity fluctuations between energy bins i and j can be cast in the following form (assuming the Poisson-noise term is the dominant contribution): The term Ω(S, Γ) accounts for the Fermi-LAT sensitivity to detect a source, and it is modeled through a stepfunction becoming equal to one at the flux threshold sensitivity S thr as described in (Manconi et al. 2020) (Section III.B.1). It depends on Γ and it includes a nuisance parameter k C P which accounts for the uncertainty in its description. We checked that a smooth, more realistic sensitivity function only has a negligible effect on the C P . The bounds in the L γ integration are L min = 7×10 43 erg/s and L max = 1×10 52 erg/s, for BLL, taken from (Ajello et al. 2014), and L min = 1 × 10 44 erg/s and L max = 1 × 10 52 erg/s, for FSRQ, taken from ). The C P 's of BLL and FSRQ are additive, i.e., C P = C BLL P + C FSRQ P . We neglect the (multipole-dependent) clustering term (discussed below in the case of DM), since we checked that, in the multipole range of interest, is a few orders of magnitude smaller than the C P term.

The source count distribution
The source count distribution, dN/dS, is defined as the number of sources per flux and solid angle and is a function of the flux S. In principle, there could also be a directional dependence, but blazars are observed up to relatively large distances, such that their distribution 10 10 10 9 10 8 10 7 S [cm 2 s 1 ] can be taken isotropic. On the other hand, populations of blazars are known to evolve with time, i.e. they depend on redshift. Furthermore, blazars do not have a unique SED: for this reason, we adopt a distribution for the photon spectral index Γ, which in equation (3) is assumed to be Gaussian.
The source count distribution (see equation (8)) depends on flux, photon spectral index, and redshift. The latter has been estimated for some of the Fermi-LAT resolved sources, for which the association to a source from another catalog was possible. We, therefore, build a three-dimensional grid for the average source count distribution in the flux bin i, the redshift bin j, and the photon spectral index bin k: This is, however, the true theoretical prediction of the dN/dS. In practice, we have to consider uncertainties of the source parameters as detected by the Fermi-LAT. This especially matters for those parameters for which the dN/dS shows a strong dependence. In our case, the first-order behavior of the dN/dS is a power law as function of S, a smoothly broken power law as function of z, and a Gaussian as function of Γ. While the power-law behavior is very smooth, the impact of resolution can be important in the Gaussian tails of Γ. We use a data-driven method to obtain the uncertainty on the determination of the spectral index. In Figure 1 we show the uncertainty of the spectral index, σ Γ , as stated in the 4FGL catalog against the flux value S. As expected, the average uncertainty increases at smaller flux values. To obtain the average uncertainty as function of S we determine the mean and RMS of σ Γ in 10 flux bins logarithmically spaced between 5 × 10 −11 cm −2 s −1 and 1 × 10 −7 cm −2 s −1 (blue data points in Figure 1) and then fit a power law (red line) to those data points: The best-fit values are A = 2.82 × 10 −6 and γ = 0.48. Assuming that to a first approximation the resolution behaves as a Gaussian, we can convolute the dN/dS of equation (10) with a Gaussian with the width σ Γ (S). The inclusion of the experimental resolution is a key ingredient to obtain a reasonable fit for small fluxes and extreme Γ bins. However, this would lead to a quadruple integral for the dN/dS which is computationally not feasible. So, instead of integrating over the flux bin, we simply evaluate the dN/dS at the (geometric) mean, S i , of the flux bin which gives a good approximation because of the smooth behavior of the dN/dS. So finally for the comparison with the 4FGL catalog, the dN/dS of our model is given by:

Fit procedure
Our data sets are the resolved sources collected in the Fermi-LAT 4FGL catalog (Abdollahi et al. 2020), from which we construct the source count dN/dS, and the anisotropy of the gamma-ray sky due to unresolved sources, encoded in the C ij P and measured in (Ackermann et al. 2018). We perform two sets of fits: one on the 4FGL catalog alone and one which combines the 4FGL catalog with the C P measurement. The fit strategy follows the procedure introduced in (Manconi et al. 2020).
In this paper, we model the luminosity function of BLLs and FSRQs separately. We assume that both blazar populations can be described by the functional form of equation 3, but with different parameter values. The log-likelihood of our fit to the 4FGL catalog data and the C P data is given by the sum of the two individual χ 2 s: Here θ BLL and θ FSRQ denote the parameters of the BLL and FSRQ models, while θ n marks nuisance parameters for the C P . In the following we also use a short notation: θ blz = {θ BLL , θ FSRQ , θ n }. We note that in practice there is only a single nuisance parameter, k CP .

Fit on the 4FGL catalog data
We extract the source count distribution of four different source classes from the 4FGL catalog and then test our model against these distributions. In more detail, we use the following four source classes: • BLL: the sum of identified and associated BL Lacs (sources labeled BLL or bll in the 4FGL catalog), • FSRQ: the sum of identified and associated FS-RQs (sources labeled FSRQ or fsrq in the 4FGL catalog), • BLZ: the sum of identified and associated blazars (sources labeled BLL, BCU, FSRQ, bll, bcu, or fsrq in the 4FGL catalog), and • ALL: all sources in the catalog.
We expect our model to respect the following constraints: • The sum of the source count distribution of our models for BLLs and FSRQs should not overshoot the observed total source count distribution of all sources (ALL). In this sense, the source count distribution ALL provides an upper bound for our models. This contribution to the χ 2 is labeled χ 2 ALL . • The sum of the source count distribution of our models for BLLs and FSRQs should be at least as large as the sum of identified and associated blazars (BLZ). However, our models are allowed to lie above the observed source count distribution because of unassociated sources. In this sense, BLZ provides a lower bound for our models (labeled χ 2 BLZ ). • Our BLL model has to explain at least the observed source count distribution BLL and thus also provides a lower bound (labeled χ 2 BLL ). • In analogy to BLL, also the FSRQ model receives a lower bound from the observed source count distribution (labeled χ 2 FSRQ ). Each of these constraints would give rise to a contribution of the 4FGL χ 2 . However, a naive sum of the χ 2 from the three lower bounds would lead to a double counting since BLL and FSRQ sources appear also in the BLZ class. To avoid this, we consider only the most constraining lower bound between the BLZ case and the combination of BLL and FSRQ, i.e., The upper bound from ALL is implemented as: Here the dN/dS in angle brackets denotes the model prediction from equation (12) and the dN/dS in round brackets is the source count distribution extracted from the 4FGL catalog. The model includes the sum of BLL and FSRQ, namely dN/dS BLZ,i = dN/dS BLL,i + dN/dS FSRQ,i . For this contribution, we integrate over all redshifts and all values of the photon spectral index. The remaining index i denotes the flux bin, as summarized in Table 1. If the flux of the bin is below the detection threshold S thr , the bin is excluded from the sum.
For the lower bounds, related to identified or associated blazars, we would like to consider also the redshift information, which is not directly provide in the 4FGL catalog. So, we extract the information from the 4LAC catalog (Ajello et al. 2020), which contains spectroscopic redshift measurements. However, the redshift information of the 4LAC catalog is incomplete, i.e. not all sources have a redshift measurement. Since we consider the identified or associated blazars only as a lower bound this does not represent a problem but it might not be the most constraining option. Hence, we consider again two cases. In the first case, we include redshift information and compare our model in bins of a two-dimensional grid in redshift and flux through: In the second case, we disregard the redshift information by integrating over the redshift. We define: Depending on the model parameter point, either equation (16) or equation (17) provides the stronger constraint. Similar to the discussion above we cannot use the sum of both χ 2 s because of double counting, and again we choose the most constraining one: As described in equation (14), we select the lower bounds by comparing χ 2 BLZ with the ones from the analysis of the individual source classes BLL and FSRQ. In this latter case, we can use the full information of the catalogs and compare models with data on a threedimensional grid in flux, redshift, and photon spectral index. Again, since the redshift information is incomplete, we define the χ 2 as the maximum of the two cases:  where M stands for either BLL or FSRQ. The individual χ 2 in the two cases are defined by when redshift information is included, and in the case with flux and spectral index bins only.
To avoid a bias from the Galactic plane we exclude small latitudes with |b| < 30 deg from our analysis.

Fit on the CP data
The fit of the APS is performed on the auto (i = j) and cross (i = j) correlation measurements, where i and j denote the energy bins. The χ 2 CP is defined as: The subscript meas refers the measured C P obtained in (Ackermann et al. 2018), while the subscript th denotes the theoretical estimation of C P calculated as in equation (9). Finally, σ 2 C ij P are the uncertainties of the measured C P , again taken from (Ackermann et al. 2018).

Analysis strategy
As anticipated above, in this work, we perform two fits. The first one utilizes only the 4FGL catalog while the second one additionally uses the C P measurement. The basic idea is as follows: From the first fit, we obtain constraints on the GLF and SED models of the two blazar populations in the flux regime of resolved point sources. From there we can extrapolate to the unresolved flux regime and calculate the C P . As will become clear in the next Section 2.5, this extrapolation agrees well with the actual C P measurement. So, we can go one step further and also perform the second fit that combines the resolved point sources (4FGL) with the C P . This combined fit provides GLF and SED models which are consistent with gamma-ray observations also below the below the flux threshold of the 4FGL catalog. When we derive DM constraints in Section 3 the combined fit serves as the baseline.
The large parameter space investigated in this work is sampled using MultiNest (Feroz et al. 2009). We use a configuration with 800 live points, an enlargement factor of efr = 0.7, and a stopping parameter of tol = 0.1. In the following, we present the results in the Bayesian statistical framework.

Results on the GLF and SED of BLLs and FSRQs
As a result of our fits, we obtain the GLFs and SEDs of FSRQs and BLLs. The 4FGL categorization of the blazars into the two classes is incomplete which leaves some degeneracy. We allow the fit to attribute the uncharacterized blazars either to BLLs or FSRQs. This is only possible because we fit both source classes at the same time. We note that this treatment leads to correlations of BLL parameters with FSRQ parameters and vice versa. These correlations are important to assess the uncertainty of the full blazar model correctly, as for example in Section 3 where blazars pose the background for our DM search.
The results are summarized in Table 2, where we state the mean values and the 1σ uncertainty derived from the marginalized posterior for each parameter. Results are provided for two setups: In the first setup, we only fit to the resolved point sources of the 4FGL catalog, while in the second setup we fit both the resolved sources and the APS data. The obtained parameter values of the two setups are compatible within their uncertainties. The GLF parameters of the FSRQ and BLL model are very well compatible with the parameters values of Ajello et al. (2012) and Ajello et al. (2014), respectively. This means that the GLF of sources belonging to the fainter regime probed in this work closely follows from the one in the brighter end. We note that we use a slightly different definition of the LDDE than Ajello et al. (2012) which leads to slightly different parameter values for p * 1 , p * 2 , and z * c . Figure 2 shows the best-fit and uncertainties of dN/dS from the combined fit of 4FGL+C P in comparison to the dN/dS extracted from the 4FGL catalog. The four different panels correspond to the different contributions to χ 2 4FGL . In more detail, the data points in the upper left panel show the dN/dS of all sources from the 4FGL catalog (at |b| > 30 deg). The sum of our models for BLLs and FSRQs is in agreement with those data. Because of the statistical technique we adopted (see above), it is expected to stay at the level or below those data points. The open white data point is below the flux threshold, so it is excluded from the analysis. The upper right panel shows the data points of the dN/dS for all identified or associated blazars. The different colors show the dN/dS in different redshift bins, while the black points are summed over all redshifts. Our model for the dN/dS of BLLs plus FSRQs lies as expected at the level or above the data points. In both upper panels, the source count distribution is integrated over all photon spectral indices, from 1.0 to 3.5. Furthermore, the upper panels only show constraints from the 4FGL catalog on the sum of the BLL and FSRQ models. The two lower panels, instead, look at the individual models for FSRQs and BLLs. Furthermore, they focus on the dN/dS for specific bins of the photon spectral index, corresponding to the peak of the distribution for each class. The lower left panel compares the dN/dS BLLs in the Γ bin from 1.9 to 2.1. Again the different colors correspond to different redshift bins and the black points contain the sum over all redshifts. Finally, the lower right panel is as the left panel but for FSRQs and a Γ bin from 2.4 to 2.6. All in all, we see that our model matches the constraints from the 4FGL catalog very well. We show these plots only for the 4FGL+C P fit but we note that they look very similar for the 4FGLonly fit. Figure 3 compares the C P measurement with the best fit model and uncertainty from the 4FGL+C P setup. The left panel shows the C P auto-correlation, while the right panel shows an example of cross-correlation, between the (8.3, 14.5) GeV energy bin and all other energy bins. The sum of BLL and FSRQ model provides a good fit to the data. This was also the case in the analysis of unresolved sources after 2 years of Fermi-LAT data Di . However, in the meantime many more point sources were resolved and the C P measurement decreased by a factor of about 10 making the latest data much more sensitive to faint populations. Still we find that blazars explain the entire latest C P measurement. The feature at 200 GeV in our model (and also, less visible, in the data) is related to a change in the analysis adopted for the measurement from (Ackermann et al. 2018). While for energies below 200 GeV bright point sources are masked from the 4FGL catalog, in the last two bins at high energies, bright point sources were masked using the 3FHL catalog. This leads to a change of the actual S thr which explains the feature. The feature does not appear in the right panel because for the cross-correlation the masks of the two energy bins involved in the measurement were joined (and so when an energy bin below 200 GeV is present, the mask is mostly provided by point sources of the 4FGL catalog). Note also from Table 2   threshold sensitivity from the reference model, is within uncertainties compatible with the default value of 1.
It is also interesting to look at the setup of 4FGL-only, and use the extrapolation of the GLF and SED model to predict the C P . As shown in Figure 4 our prediction agrees very well with the measurement. We note that the C P is dominated by FSRQs at low energies, below ∼ 2 GeV, while BLLs dominate at higher energies. The domination of BLLs at high energies is expected since they have a harder SED than FSRQs. The fact that there is a transition from FSRQs to BLLs in the C P at low energies introduces a softening in the spectral index at low energies. This softening has previously been interpreted as a possible hint for a new source population (Ando et al. 2017). The new C P data (Ackermann et al.  2018) and the detailed treatment presented here, allow us to interpret it in terms of FSRQ. Figure 4 already hints that both BLLs and FSRQs are required to describe the C P data. We quantify this statement a bit better in the following. Using the results of the 4FGL-only fit, we calculate the posterior distribution of the χ 2 CP assuming (i) the sum of FSRQs and BLLs, labeled all, and (ii) only BLLs. An only FS-RQs hypothesis is excluded since it cannot explain the high-energy C P data. In order to consider the systematic uncertainty on the exact flux threshold, we profile over the normalization of the C p . The results in Fig. 5 show that the sum of FSRQ and BLL is preferred. Using the full posterior of the 4FGL-only fit and defining  L Cp = exp(−χ 2 Cp /2), we calculate the Bayes factors of the hypotheses (i) and (ii) obtaining 4.0 × 10 3 . More details about the calculation of the Bayes factors are given in App. B. In this sense, the two physical populations (BLLs and FSRQs) are clearly preferred over a single population of only BLLs or FSRQs. We note, however, that the preference for these two population is based on the catalog prior. As discussed in the Appendix A, the C P data by itself is not sufficient to distinguish between a scenario with one or two populations. Without the catalog prior, a hypothetical and more general GLF and SED can provide a good fit to the C P data.
Finally, Figure 6 shows a triangle plot with posterior distributions for parameters of the GLF and SED for the BLL and FSRQ models. The posteriors correspond to the 4FGL+C P setup. The diagonal contains the marginalized one-dimensional posteriors for each individual parameter, while the panels in the lower half show the 1 and 2σ contours for each combination of two parameters. Since the BLL and FSRQ models have the same functional form we can combine them into the same triangle, using different colors. We observe that the SED parameters µ * and σ are well constrained for both populations. The average photon spectral index of BLLs is µ * ∼ 2.0 and a width of σ ∼ 0.2. As expected, FSRQs follow a softer energy spectrum with an average index with µ * ∼ 2.5 but a similar width. Also, the shape of the GLF at small L is reasonably constrained. The index γ 1 lies between 0.4 and 1.0 for FSRQs and between 0.5 and 1.2 for BLLs, while the behavior at large L (see γ 2 ) is less constrained. We note the degeneracy between A and L * . This is because, at first order, in both cases, their main impact in the fit is to change the normalization of the GLF. The redshift dependence is only weakly constrained due to degeneracies with other parameters.
We provide the covariance matrix of our fits in the ancillary files (arXiv version). This covers, to a first approximation, the degeneracies and correlations of the fit parameters. As a final comment, let us note that, clearly, the uncertainties on the parameters of BLLs and FSRQs show some level of mutual correlation in the fit. It is just for the sake of clearness that we do not show the entire triangle plot in Figure 6.

Blazar contribution to the UGRB
We see that the measured gamma-ray angular correlations require the presence of FSRQs and BLLs. Populations with a GLF peaked at lower luminosities (like misaligned AGN and star forming galaxies) cannot account for the C P data. Thus, FSRQs and BLLs provide an unavoidable contribution to the UGRB intensity. In   2015) with the prediction from our models. The measurement accounted for contribution of point sources from the 2FGL catalog (Nolan et al. 2012). To be consistent, we apply a flux threshold corresponding to the 2FGL catalog, taken from Ackermann et al. (2015), for the predictions in Figure 7. We conclude that blazars provide a significant contribution to the UGRB, accounting for about 30% between 10 and 100 GeV. At energies below 1 GeV the contribution decreases to about 20%.

BOUNDS ON WIMP DARK MATTER
The UGRB observed by the Fermi-LAT could conceal a signal from DM particles. We focus our analysis on annihilating dark matter. Since annihilation occurs universally in all dark matter structures, we need to consider the gamma-ray emission from both the halo of our Galaxy and from extragalactic structures. We, therefore, have two dark matter contributions to the anisotropy APS, one for the Galactic halo and one from the extragalactic dark matter distribution (the two contributions are not expected to cross-correlate). Moreover, extragalactic structures host the same astrophysical sources (BLL and FSRQ) which we have discussed in the previous sections. This induces a cross-correlation term in the APS between extragalactic dark matter and sources emission. All these terms are properly modeled and taken into account in our analysis, as outlined below. In the following, we will discuss the contribution of DM only to the anisotropy signal C P , since the contribution of dark matter halos to the source count dis-tribution is significantly suppressed as compared to the one arising from astrophysical sources in the resolved flux regime. A DM contribution to the dN/dS could in principle emerge on top of astrophysical sources only at very low fluxes.

Extragalactic Dark Matter modeling
The APS of the cross-correlation between a sourcefield X in the energy bin i and a source-field Y in the energy bin j reads (Fornengo & Regis 2014) where W X (χ) is the window function of field X, P XY (k, χ) the 3-dimensional cross power-spectrum of the fluctuations of the two fields and χ denotes the comoving distance, related to redshift by dχ = (c dz)/H(z).

The window function for annihilating DM is given by (Ando & Komatsu 2006; Fornengo & Regis 2014)
where Ω DM and ρ c are the present-day cosmological abundance of DM and the critical density of the Universe, respectively, m DM is the mass of the DM particle, ∆ 2 (z) is the clumping factor, and σ ann v denotes the velocity-averaged annihilation cross-section of dark matter particles, assumed here to be the same in all dark matter halos. dN ann /dE indicates the number of photons produced per annihilation as a function of energy and sets the gamma-ray energy spectrum and τ (E, z) denotes the optical depth of gamma-ray photons, which we model as in (Finke et al. 2010b).
To determine the DM auto-correlation, we compute the 3D power spectrum with the so-called halo-model approach (see, e.g., the review in Cooray & Sheth (2002)). For X = Y = DM, we have where dn h /dM is the halo mass function, P lin (k, z) is the linear matter power spectrum, b h (M ) is the linear bias, andû ann (k|M, z) denotes the Fourier transform of the density profile of the DM halos (see, e.g., see the Appendix of Cuoco et al. (2015)). We assume the NFW DM density profile (Navarro et al. 1996). All the ingredients in equations (24) and (25)  To characterize the halo profile and the subhalo contribution, we need to specify their mass concentration. The description of the concentration parameter c(M, z) at small masses and for subhalos is still an open issue and provides our largest source of uncertainty. In the following, we consider two models that we name 'LOW', where we take the description of c(M, z) from (Correa et al. 2015), and 'HIGH', where we follow (Neto et al. 2007). They differ for what concerns the extrapolation of the concentration at low masses and lead to a difference of about one order of magnitude in the final bounds on the annihilation cross section, as shown in Section 3.4.
Clearly, the distribution of extragalactic DM halos (and in turn of the annihilation signal) has some level of correlation with the blazar distribution, therefore inducing a cross-correlation signal between DM and blazars. The blazar window function can be phrased as: W BLA (z, E) = χ(z) 2 f S , with the mean flux defined as: The 3D power spectrum of the cross correlation between annihilating DM and blazars is given by: where b S is the bias of blazars with respect to the matter density, for which we adopt The relation M (L γ , z) between the mass of the host halo and the luminosity of the hosted blazar is taken from (Camera et al. 2015).

Galactic Dark Matter modeling
For the modeling of the signal expected from the Galactic subhalos, we followed the treatment of Ando (2009). In general, the prediction of the APS for the Galactic component is the sum of two contributions, one arising from the main halo and the other originated by substructures. The main (smooth) halo contribution is subdominant for multipoles above a few: since we are dealing with multipoles larger than 50, it is here neglected. For the substructure contribution, we consider an anti-biased subhalo distribution, corresponding to the fiducial model A1 of Ando (2009), with a boost factor of subhalos set to unity.
The subhalo number density as a function of the distance f from the center of the Galaxy reads where M vir,MW is the Milky Way virial mass, α E = 0.68, γ is the lower incomplete gamma function, r −2 = 0.81 r 200,MW , c −2 = r vir /r −2 and the fraction f of DM enclosed in subhalos is fixed to 0.2. The minimal subhalo mass is set at M min = 10 −6 M . We do not include a truncation for the subhalo distribution at large radii. The angle-average number density referred to our position in the Galaxy is where r (cos(ψ)) = r 2 0 + s 2 − 2 r ⊕ s cos(ψ), r ⊕ = 8.5 kpc is our distance from the center of the Galaxy, s represents the distance along the line of sight, and ψ is the angle between the directionn of observation and the direction to the Galactic center.
Numerical simulations suggest that the mass distribution of subhalos follows a power-law behavior with the mass M of the subhalo and can be written as where α 0 = 1.9. The subhalo luminosity L for a subhalo with a mass M depends on the particle properties of DM, σ ann v and m DM , as well as on the energy spectrum dN ann /dE associated to the channel under study. It reads With the above ingredients, the APS for the Galactic subhalo contribution can then be expressed as where i and j refer to energy bins. The integral along the line of sight is performed up to s max = r vir,MW = 258 kpc, and starts at s min = L/(4πS thr ). In this specific case, we use a simplified approach: Instead of the full, i.e., Γ-dependent flux threshold, we use a fixed threshold averaged over Γ. More specifically, the thresholds (corresponding to the flux from 1 to 100 GeV) are set to S thr = 10 −10 cm −2 s −1 for the first ten energy bins (4GFL threshold) and 2 × 10 −10 cm −2 s −1 for two highest energy bins (3FHL threshold). The value of s min is thus chosen such that only unresolved Galactic subhalos are considered in the determination of the APS. We adopted an NFW profile for the internal density distribution of the subhalos andũ sh denotes its Fourier transform. The maximal halo mass for subhalos in our Galaxy, M max = 10 10 M .
Summarizing, the APS involving DM is given by the sum of four terms, three extragalactic (the autocorrelation from extragalactic DM halos and the crosscorrelation with BLLs and FSRQs) and the Galactic one, which is not expected to cross-correlate with extragalactic source populations. Since the DM contributions are not flat in multipole but -dependent, they are averaged in the multipole range considered for the determination of the C P measurement in (Ackermann et al. 2018).

Statistical framework
We derive bounds on the DM annihilation cross section as a function of the DM mass by marginalizing over the uncertainties in the astrophysical background model. This means that we are not using the approximation of a fixed background model obtained from the best-fit of the blazars-only case, on top of which the DM contribution is added.
The naive approach to consider the full uncertainty would be to perform an extended parameter scan which would include at the same time the blazar parameters and the DM parameters. However, this is computationally very expensive. We instead use a method called importance sampling in order to recycle the information from background-only fits which significantly speeds up the calculation. The same approach has recently been used in a different context (Kahlhoefer et al. 2021).
One side-product of the MultiNest scan is a set of parameter vectors that follows the multidimensional posterior distribution. This set is provided in the socalled equal-weights sample. We will apply importance sampling to obtain the posterior distribution of the full parameter space (blazars and DM) by using the equalweights set of the fit to only blazars.  (38). The left panel shows the constraints in the "LOW" scenario, while the right panel reports the comparison between "LOW" and "HIGH" cases.
First, we note that we can approximate the integral of the product of background, i.e. the blazar, posterior and an arbitrary function f over the blazar parameters, θ blz , by a sum over the parameter points in the set of the equal-weights sample: where L 0 is the likelihood, Z 0 is the evidence, and p 0 is the prior. The subscript 0 indicates quantities that refer to the fit without DM. We note that the factor L 0 (θ blz ) p 0 (θ blz )/Z 0 is by definition the posterior distribution. Furthermore, we know that the integral over the prior p 0 (θ blz ) is normalized to 1. If we set f (θ blz ) = Z 0 /L 0 (θ blz ) we see that the evidence is given by We can obtain the marginalized likelihood by integrating over the background parameters: Here θ DM = {m DM , σ ann v } denotes the DM parameters, L(θ blz , θ DM ) is the likelihood of the full parameter space, and p(θ blz ) is the prior of the blazar parameters. Using equation (34) and assuming the same prior (p(θ blz ) = p 0 (θ blz )), we see that the integral of equa-tion (36) is approximately given by the sums In the next step, we can turn this equation into an expression for a marginalized χ 2 by using the definition of our likelihood (L = exp(−χ 2 /2)) Finally, we obtain the DM limit at the 95% C.L. from the requirement ∆χ 2 (θ DM ) ≤ 3.84.
In practice, we evaluate equation (38) on a grid of m DM with 14 grid points logarithmically spaced between 10 GeV and 4 TeV. We note that the annihilation cross section σ ann v only changes the normalization of the C P contributions but not the shape. The DM×DM contributions scales with σ ann v 2 while the DM×BLZ contributions scales linearly with σ ann v . So, we can tabulate the C P contributions of background and DM at a reference value of σ ann v = 3 × 10 −26 cm 3 /s. The equalweights set contains O(10 5 ) parameter vectors. After the tabulation, the evaluation of equation (38) takes about one second. The importance sampling reduces the computing time by about a factor of ten since the typical number of evaluations required for a full parameter scan of blazar and DM parameters requires O(10 6 ) evaluations. A further advantage of the importance sampling is that the tabulation can be parallelized to an arbitrary degree which is not possible for a Monte Carlo based parameter sampling with MultiNest. FSRQxDM BLLxDM Figure 9. Comparison of the CP from blazars and DM. The DM contributions are shown at 100 GeV (left) and 1 TeV (right), and computed in the "LOW" scenario, for annihilation intobb. In both panels we choose the value of σannv of DM at the limit derived in Figure 8.

Constraints on DM annihilation
We derive constraints on the annihilation of DM into a pair ofbb quarks which serves as an illustrative example. The limits for other hadronic channels are expected at a similar level. The left panel of Figure 8 shows the marginalized ∆χ 2 in the plane of DM mass and annihilation cross section as derived from equation (38). For DM masses between 15 and 140 GeV a small DM contribution slightly improves the fit of the C P data. However, it is statistically not significant. The maximal improvement of the ∆χ 2 is ∼ 3.5 which corresponds to a local significance of less than 2σ and an even smaller global significance. Consequently, we can derive DM limits as a function of the DM mass. In the fiducial setup, namely, using the "LOW" model for the concentration-mass relation of DM halos, see Section 3, we can place an upper limit of σ ann v = 10 −25 cm 3 /s on the annihilation cross section at the DM mass of 10 GeV. The limit gradually weakens to σ ann v = 3 × 10 −23 cm 3 /s at 4 TeV. In a more, aggressive setting for the concentration parameter, i.e., the "HIGH" model, we obtain a DM limit which is almost one order of magnitude stronger. The plot in the left panel of Figure 8 shows the limits for the "LOW" model, while the comparison between the two cases is shown in the right panel. In the "HIGH" scenario, we can exclude a thermal WIMP for m DM < 20 GeV.
As explained above the measurement of the APS is dominated by the Poisson noise term. For this reason, the DM bounds in Figure 8 are weaker than the ones from other probes of the UGRB (Di Mauro & Donato 2015;Charles et al. 2016), such as the total intensity energy spectrum and the cross-correlation APS with gravitational tracers, which are less affected by the noise be-ing linear instead of quadratic probes of the UGRB (see, e.g., Figure 4 in Regis et al. (2015)). There are also other strategies of indirect DM searches using gamma-rays, e.g. from the dwarf spheroidal galaxies or the Galactic Center, or using cosmic-ray antiprotons. Those DM limits are typically stronger by 1-2 orders of magnitude (see e.g. Leane (2020); Slatyer (2021) and references therein), however, those analyses are affected by different systematic uncertainties.
In Figure 9, we show the contribution of all the different components to the C P for two exemplary DM masses of 100 GeV (left panel) and 1 TeV (right panel), in the "LOW" scenario. The DM components are evaluated at the σ ann v values corresponding to the limit shown in Figure 8 for that DM mass. As expected, the blazar components are dominant and DM only provides a sub-dominant part. The largest DM contribution stems from the extragalactic DM halos, closely followed by the contribution of Galactic DM subhalos which is roughly smaller by a factor of 2, while the contribution of BLZ×DM is smaller by 1-2 orders of magnitude. The uncertainty bands in Figure 9 represent the 1σ uncertainty in the blazar background models. In the "HIGH" scenario, the picture is similar but with the DM contribution strongly dominated by the extragalactic term.
Finally, we stress again that the absence of a DM signal in the C P is a further confirmation that the sum of FSRQs and BLLs fully explains the entire UGRB anisotropy and no additional, weak component is required to match the data.

CONCLUSIONS
In this work, we compared models of the GLF and SED of blazars to the latest measurement of the energy spectrum of the UGRB anisotropies (Ackermann et al. 2018) and the properties of the resolved gamma-ray sources of the Fermi-LAT 4FGL catalog. We considered two different blazar populations, distinguishing between BL Lacs (BLLs) and flat spectrum radio quasars (FS-RQs). We found that BLLs and FSRQs can account for the totality of the UGRB anisotropy, with BLLs dominating the APS at high energies, while FSRQs being important at GeV energies. The derived models well reproduce the size and spectral features observed by (Ackermann et al. 2018), and the properties of source number counts of the 4FGL catalog. Our analysis significantly constrains the redshift and luminosity dependence of the blazar GLF and the spectrum of the SED in the unresolved regime. We also calculate the contribution of the unresolved population of FSRQs and BLLs to the the UGRB intensity spectrum finding a non-negligible contribution of about 30% between 10 and 100 GeV and about 20% at 1 GeV.
In the second part of the paper, we included a contribution to the UGRB arising from annihilating DM and performed a global fit to derive constraints on the particle DM parameters. We computed both Galactic and extragalactic DM contributions, and included cross-terms in the APS, due to the cross-correlation of blazars with the DM halos hosting them. The dominant term arises from extragalactic DM halos, which strongly depends on the poorly-known description of DM subhalos. To bracket the uncertainty, we considered two different scenarios, "LOW" and "HIGH", which lead to an upper limit of σ ann v = 10 −25 cm 3 /s and σ ann v = 1.5 × 10 −26 cm 3 /s, respectively, on the annihilation cross section at the DM mass of 10 GeV, for annihilation into bottom quarks.
The present analysis of the UGRB anisotropies is based on the measurement of (Ackermann et al. 2018), where no evidence for an -dependent APS was found. Further data, and more resolved sources, would allow to reduce the level of the Poisson noise C P and to measure an APS unveiling the large-scale clustering of gammaray sources. This would allow us to deepen our understanding of blazar populations, as well as exploit the APS observable in a much more powerful way in the context of DM bounds. In this paper, we have considered two blazar populations, BLLs and FSRQs. Both populations are important and their sum provides a good fit to angular correlations of the UGRB. In more detail, FSRQs give the largest contributions to the C P below a few GeV, while BLLs dominate at higher energies. As a consequence, we found that the FSRQs are responsible for the softening of the C P at low energies, which is an interesting result because this softening has also been interpreted as a possible hint for a new source population (Ando et al. 2017). In a similar spirit, also (Ackermann et al. 2018) claimed that the C P data indicates a hint for two populations. On the other hand, in (Manconi et al. 2020) a good fit of the C P data is obtained with one population using a model which combines FSRQs and BLLs into a single GLF and SED model.
Here, we explore the issue of one versus two populations more thoroughly by employing a phenomenological model for the SED and GLF of two hypothetical source populations. The phenomenological model is slightly simplified as compared to the physical models for the GLFs and SEDs of the BLL and FSRQ source populations used in the main text. The most important difference is that the phenomenological model does not include a redshift dependence. However, the C P data is not sensitive to redshift information because it is integrated along the line of sight. We find that our model is sufficient to explore the question about the number of source populations. Furthermore, it has the clear advantage that the large parts of the computations can be done analytically. We find that the width of the spectral index distribution is a key parameter to distinguish the scenarios of one or two populations.

A.1. Definition of the GLF and SED
The source count distribution as a function of the photon flux S is modeled as a power law: where A is an overall normalization, γ is the power-law index, and S 0 is a reference flux. In the following, S 0 is fixed to 1 × 10 −10 cm 2 s −1 and the fluxes S and S 0 always refer to the photon flux in the energy bin from 1 GeV to 100 GeV. To relate the flux S to the flux S i in a different energy bin i we need the SED, which we model as a power law with an exponential cutoff: Here K is the normalization, E 0 a reference energy, Γ the photon spectral index, and E c the energy of the exponential cutoff. The exponential cutoff allows mimicking the attenuation of gamma rays at high energies. By definition, the flux S is given by S = 100 GeV 1 GeV dE dN/dE. So, the ratio of the fluxes S i and S is given by: Finally, we allow for an intrinsic distribution of the photon spectral indices following a Gaussian with mean µ and width σ: Then, the C P of unresolved point sources is given by: where Ω(S, Γ) is the detector efficiency to resolve point sources, which we model as a θ-function at a Γdependent flux threshold, S thr (see main text). By considering photon spectral indices between 1 and 3, the C P is then calculated as: We note that all the energy-information of the C P is encoded in the factor s i s j such that to a first approximation the slope, γ, of the dN/dS is degenerate with the overall normalization, A. So, we fix γ to a benchmark value of 2.2 in the following analysis.   Figure 10. Triangle plots showing the best fit region of the phenomenological model fitted to the CP data for different fit setups. On the diagonal, we display the marginalized likelihood for all parameter and the remaining panels below the diagonal contain the 1 and 2 σ contours derived from the marginalized likelihood in two dimensions for each parameter combination. The red contours and lines correspond to a fit with a single population, while the green and blue contours correspond to the two different populations of a fit containing two populations (labeled a and b). Left panel: Phenomenological model with a fixed photon spectral index. Right panel: Phenomenological model with a Gaussian distribution of the photon spectral index.

A.2. Fits and results for the phenomenological model
In this section, we perform a total of four fits. The fits differ by the number of source populations (one or two) and by the inclusion of a distribution in the photon spectral index.
In the first two fits, we neglect the distribution of spectral indices, namely we force σ = 0. Then, the first analysis employs a single source population with the SED and GLF specified in the previous paragraph. The free parameters are the normalization of the GLF (A), the photon spectral index (µ), and the energy cutoff of the SED (E c ). Furthermore, we use a nuisance parameter (k CP ) that varies the value of flux threshold (see main text for more details). This fit has 4 free parameters.
In the second analysis, we use two source populations, labeled a and b. They have the same functional form, but different parameters: in particular, for the normalization of the GLF and the photon spectral index. However, we force them to have the same value for the cutoff, E c . Together with the nuisance, this amounts to 6 free parameters. To avoid an extra degeneracy in the fit we restrict the photon spectral indices to µ a < µ b .
The third and fourth fits are very similar to the first two, but we allow for a distribution in Γ. Consequently, this leads to one (σ) and two (σ a and σ b ) more parameters in the fits, respectively. We use the MultiNest code to perform the four fits and to obtain the posterior distributions presented in the following. Uncertainties are stated in the Bayesian statistical framework.
The results of our four fits are presented in Table 3 and in Figure 10. The χ 2 /dof of all four fits are close to 1 and we conclude that all of them give a good fit to the C P data. However, when comparing the χ 2 s of the first and CP energy-correlation matrix for the phenomenological model in the setup with only one population but including a Γ-distribution. second fit, i.e. the fits without the Γ distribution, we see that including a second population reduces the total χ 2 by 10.3, which formally corresponds to a statistical shift by 2.8σ. So, this could be interpreted as a small hint for two populations. However, if we look at the results with a Γ distribution the conclusion changes. Here the χ 2 only decreases by 0.2 when a second population is introduced, which is not significant. So, we cannot distinguish between one or two populations based only on the C P data. This can also be seen from the triangle plots in Figure 10. Without allowing for a Γ distribution, the red posterior contours (one population) are not fully included in the contours of either the green or blue (population a or b of two population fit). Furthermore, both populations, a and b are required to have a relatively high normalization of A > 10 −11 cm −2 s −1 sr −1 . On the other hand, if we allow for a Γ distribution, the red contours are fully included in both the green and the blue contours. And, in the case of 2 populations, one of the normalizations A can be pushed to negligible values of 10 −12 cm −2 s −1 sr −1 . We compare the C P energy autocorrelations of the best-fit to data in Figure 11. It is clearly visible that, in the case without Γ distribution, two populations provide a better fit, while in the case with the Γ distribution there is almost no difference.
Finally, we have a look at the C P energy-correlation matrix (coefficients defined as C ij P / C ii P C jj P ). By definition, the diagonal coefficients are equal to one. If there was only one source population with a single and fixed photon spectral index, also the off-diagonal coefficients would be equal to 1. Instead, the C P measurement (Ackermann et al. 2018) showed that the off-diagonal coefficients are ∼ 0.6, which was interpreted as an indication for two populations. In Figure 12, we show that also a single source population with a distribution of photon spectral indices provides the correct off-diagonal pattern.

A.3. Conclusions
We conclude that, based on the C P measurement itself, it is not possible to distinguish between two populations with narrow spectral index distributions and a single population with a broader distribution. In the latter case, the required value for the width of the Γ distribution is σ = 0.34 +0.07 −0.08 . Additional information can help to solve the riddle. In the analysis of the main text, we included a physical model and the constraining power from the 4FGL catalog. This allowed us to break the degeneracy and to determine the presence of two populations, BLL and FSRQ.

B. CALCULATION OF THE BAYES FACTOR
In the main text, we use a physical model with two population, BLL and FSRQ. Because of the 4FGL catalog prior, this model is significantly preferred over a model with only one BLL populations. This can be quantified by the calculation of Bayes factors, which is biefly summarized here. The Bayes factor is defined as the ratio of evidences between the two hypotheses: two populations (labeled all) and one population (labeled BLL). For flat priors, it becomes B = Z all Z BLL = dθ blz L 4FGL (θ blz )L CP,all (θ blz ) dθ blz L 4FGL (θ blz )L CP,BLL (θ BLL ) .
Here, L CP,all and L CP,BLL are the likelihoods for the C P data using the sum of BLLs plus FSRQs and only BLLs, respectively. We note that the C P likelihoods also depend on the nuisance parameters θ n which we suppress in the notation. In analogy to Eq. (34), we can approximately replace the integrals by sums, where i runs over the all points in the equal-weights sample of the 4FGL-only fit. Finally, we note one subtlety: The normalization of the C P model strongly depends on the exact value of the catalog threshold which is only known approximately. To take the uncertainty due to the threshold into account, we profile over the C P normalization, individually for each parameter point i. The mean renormalization factors are 1.04 and 1.16 for the hypotheses of two populations and one population, respectively.

C. RESOLVING PREVIOUS UGRB ANISOTROPY MEASUREMENTS
As mentioned in the main text, the UGRB emission is exposure-dependent: the more the LAT observes, the  Fornasa et al. 2016, PRD, 94, 123005 Ackermann et al. 2018. Comparison between the measurements from (Fornasa et al. 2016) (red points) and (Ackermann et al. 2018) (black points). The CP's of the populations of FSRQs, BLLs, and BCUs of the 4FGL not yet resolved in the 3FGL, are shown in orange, blue and green, respectively, and the best-fit to the CP of (Ackermann et al. 2018) is given by the black curve. Their sum is shown with the red curve, which well-reproduces the measurement of (Fornasa et al. 2016). The comparison is restricted to energies below 10 GeV, see text.
more sensitive is the survey to fainter sources, and consequently less unresolved emission contributes to the UGRB. We show this effect by comparing the UGRB anisotropy energy spectrum measured by (Fornasa et al. 2016), which masked sources from the 3FGL source catalog (Acero et al. 2015) based on 4 years of LAT survey, with the one in (Ackermann et al. 2018). In Figure 13 we report both measurements in red and black, respectively. The difference between the two measurements, which consider similar total mission times (7.5 years in (Fornasa et al. 2016) and 8 years in (Ackermann et al. 2018)), is mainly due to the difference in the number of resolved sources masked away from the analyzed maps: the 4FGL counts approximately 2000 more sources than the 3FGL. We test this by evaluating the anisotropy energy spectrum of the populations of FSRQs, BLLs, and BCUs of the 4FGL not yet resolved (i.e., not present) in the 3FGL (let us call it ∆C P,4FGL−3FGL ). We verify that the sum of the best-fit model of the measurement in (Ackermann et al. 2018) plus the ∆C P,4FGL−3FGL reproduces the anisotropy energy spectrum as measured by (Fornasa et al. 2016). In agreement with our expectations, the results in Figure 13 show that there is a transition between FSRQs and BLLs in the ∆C P,4FGL−3FGL when going from lower to higher energies. So it is plausible to observe a similar transition in the C P measurement of (Ackermann et al. 2018).