Flat-spectrum Radio Quasars and BL Lacs Dominate the Anisotropy of the Unresolved Gamma-Ray Background

We summarize here the parameterizations of the gamma-ray luminosity function (GLF) and the SED of the BLLs and FSRQs adopted in our analysis. For more details, see Ajello et al. (2012) and Ajello et al. (2014). The GLF Φ(L_γ , z, Γ) = d³ N/dL_γ dVdΓ—defined as the number of sources per unit of luminosity L_γ , the comoving volume V at redshift z, and the photon spectral index Γ—is typically decomposed in terms of its expression at z = 0 and a redshift evolution function:

$$\begin{eqnarray}&&{\rm{\Phi }}({L}_{\gamma },z,{\rm{\Gamma }})={\rm{\Phi }}({L}_{\gamma },0,{\rm{\Gamma }})\times e({L}_{\gamma },z),\end{eqnarray} \tag{ 1 }$$

where L_γ is the rest-frame luminosity in the energy range (0.1–100) GeV, i.e., ${L}_{\gamma }={\int }_{0.1\,\mathrm{GeV}}^{100\,\mathrm{GeV}}{\rm{d}}{E}_{r}\,{ \mathcal L }({E}_{r})$, with

$$\begin{eqnarray}&&{ \mathcal L }({E}_{r})=\displaystyle \frac{4\pi {d}_{L}^{2}(z)}{(1+z)}E\,\displaystyle \frac{{\rm{d}}N}{{\rm{d}}E},\end{eqnarray} \tag{ 2 }$$

with E being the observed energy, related to the rest-frame energy E_r as E_r = (1 + z)E. The comoving volume element in a flat homogeneous universe is given by d² V/dΩdz = c χ²(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 a ΛCDM cosmology, with parameters from the final full-mission Planck measurements of the cosmic microwave background anisotropies (Aghanim et al. 2020).

At redshift z = 0, the parameterization of the GLF is

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Phi }}({L}_{\gamma },0,{\rm{\Gamma }}) & = & \displaystyle \frac{A}{\mathrm{ln}(10){L}_{\gamma }}{\left[{\left(\displaystyle \frac{{L}_{\gamma }}{{L}_{0}}\right)}^{{\gamma }_{1}}+{\left(\displaystyle \frac{{L}_{\gamma }}{{L}_{0}}\right)}^{{\gamma }_{2}}\right]}^{-1}\\ & & \times \exp \left[-\displaystyle \frac{{\left({\rm{\Gamma }}-\mu ({L}_{\gamma })\right)}^{2}}{2{\sigma }^{2}}\right],\end{array}\end{eqnarray} \tag{ 3 }$$

where A is a normalization factor, the indices γ₁ and γ₂ 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 the 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 μ*:

$$\begin{eqnarray}&&\mu ({L}_{\gamma })={\mu }^{* }+\beta \left[\mathrm{log}\left(\displaystyle \frac{{L}_{\gamma }}{\mathrm{erg}\,{{\rm{s}}}^{-1}}\right)-46\right].\end{eqnarray} \tag{ 4 }$$

On the other hand, the FSRQs have a narrower distribution, meaning that the inclusion of this effect would not affect the fit, so we can fix μ(L_γ ) = μ*.

We adopt a luminosity-dependent density evolution (LDDE):

$$\begin{eqnarray}\begin{array}{rcl}e({L}_{\gamma },z) & = & \left[{\left(\displaystyle \frac{1+z}{1+{z}_{c}({L}_{\gamma })}\right)}^{-{p}_{1}}\right.\\ & & +{\left.{\left(\displaystyle \frac{1+z}{1+{z}_{c}({L}_{\gamma })}\right)}^{-{p}_{2}}\right]}^{-1},\end{array}\end{eqnarray} \tag{ 5 }$$

with ${z}_{c}{({L}_{\gamma })={z}_{c}^{* }\times ({L}_{\gamma }/{10}^{48}\mathrm{erg}\,{{\rm{s}}}^{-1})}^{\alpha }$, ${p}_{1}({L}_{\gamma })={p}_{1}^{* }+\tau \times (\mathrm{log}({L}_{\gamma })-46)$, and ${p}_{2}({L}_{\gamma })={p}_{2}^{* }+\delta \times (\mathrm{log}({L}_{\gamma })-46)$. We set δ to 0.64 for both populations (Ajello et al. 2015), 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:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}N}{{\rm{d}}E}\sim {\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{-{\rm{\Gamma }}}.\end{eqnarray} \tag{ 6 }$$

For definiteness, the spectral 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

$$\begin{eqnarray}&&S({E}_{\min },{E}_{\max })={\int }_{{E}_{\min }}^{{E}_{\max }}\displaystyle \frac{{\rm{d}}N}{{\rm{d}}E}{e}^{-\tau (E,z)}\,{\rm{d}}E,\end{eqnarray} \tag{ 7 }$$

where τ(E, z) describes the attenuation by the extragalactic background light (Finke et al. 2010). Unless explicitly stated otherwise, the flux in the following computations of the dN/dS and the associated figures always refers to the energy bin from 1 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 to hand, we can compute the differential number of blazars per integrated flux and solid angle as

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}N}{{\rm{d}}S}={\int }_{0.01}^{5.0}{\rm{d}}z{\int }_{1}^{3.5}{\rm{d}}{\rm{\Gamma }}\,{\rm{\Phi }}[{L}_{\gamma }(S,z,{\rm{\Gamma }}),z,{\rm{\Gamma }}]\,\displaystyle \frac{{\rm{d}}V}{{\rm{d}}z}\,\displaystyle \frac{{\rm{d}}{L}_{\gamma }}{{\rm{d}}S},\end{eqnarray} \tag{ 8 }$$

and the size of the gamma-ray intensity fluctuations between the energy bins i and j can be cast in the following form (assuming that the Poisson noise term is the dominant contribution):

$$\begin{eqnarray}\begin{array}{rcl}{C}_{{\rm{P}}}^{{ij}} & = & {\displaystyle \int }_{0.01}^{5.0}{\rm{d}}z\displaystyle \frac{{\rm{d}}V}{{\rm{d}}z}{\displaystyle \int }_{1}^{3.5}{\rm{d}}{\rm{\Gamma }}{\displaystyle \int }_{{L}_{\min }}^{{L}_{\max }}{\rm{d}}{L}_{\gamma }\,{\rm{\Phi }}({L}_{\gamma },z,{\rm{\Gamma }})\\ & & \times {S}_{i}({L}_{\gamma },z,{\rm{\Gamma }})\,{S}_{j}({L}_{\gamma },z,{\rm{\Gamma }})\left[1-{\rm{\Omega }}(S({L}_{\gamma },z,{\rm{\Gamma }}),{\rm{\Gamma }})\right].\end{array}\end{eqnarray} \tag{ 9 }$$

The term Ω(S, Γ) accounts for the Fermi-LAT sensitivity to detect a source, and it is modeled through a step function 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}_{{\rm{P}}}}$, which accounts for the uncertainty in its description. We checked that a smoother, more realistic sensitivity function only has a negligible effect on C_P. The bounds in the L_γ integration are L_min = 7 × 10⁴³ erg s⁻¹ and L_max = 1 × 10⁵² erg s⁻¹, for BLLs, taken from Ajello et al. (2014), and L_min = 1 × 10⁴⁴ erg s⁻¹ and L_max = 1 × 10⁵² erg s⁻¹, for FSRQs, taken from Ajello et al. (2012).

The C_Ps of the BLLs and FSRQs are additive, i.e., ${C}_{{\rm{P}}}={C}_{{\rm{P}}}^{\mathrm{BLL}}+{C}_{{\rm{P}}}^{\mathrm{FSRQ}}$. We neglect the (multipole-dependent) clustering term (discussed below in the case of DM), since we checked that, in the multipole range of interest, it is a few orders of magnitude smaller than the C_P term.

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 can be taken as 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 the flux, photon spectral index, and redshift. The latter has been estimated for some of the Fermi-LAT resolved sources—those for which association with a source from another catalog was possible. Therefore, we 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:

$$\begin{eqnarray}&&\begin{array}{l}{\left\langle \displaystyle \frac{{\rm{d}}N}{{\rm{d}}S}\right\rangle }_{{ijk}}=\displaystyle \frac{1}{{S}_{i,\max }-{S}_{i,\min }}\\ \,\,\times \,{\displaystyle \int }_{{S}_{i,\min }}^{{S}_{i,\max }}{\rm{d}}S{\displaystyle \int }_{{z}_{j,\min }}^{{z}_{j,\max }}{\rm{d}}z{\displaystyle \int }_{{{\rm{\Gamma }}}_{k,\min }}^{{{\rm{\Gamma }}}_{k,\max }}{\rm{d}}{\rm{\Gamma }}\,{\rm{\Phi }}({L}_{\gamma },z,{\rm{\Gamma }})\,\displaystyle \frac{{\rm{d}}V}{{\rm{d}}z}\,\displaystyle \frac{{\rm{d}}{L}_{\gamma }}{{\rm{d}}S}.\end{array}\end{eqnarray} \tag{ 10 }$$

This is, however, the true theoretical prediction of the dN/dS. In practice, we have to consider the uncertainties of the source parameters, as detected by 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 a function of S, a smoothly broken power law as a function of z, and a Gaussian as a function of Γ. While the power-law behavior is very smooth, the impact of the 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 a function of S, we determine the mean and rms of σ_Γ in 10 flux bins logarithmically spaced between 5 × 10⁻¹¹ cm⁻² s⁻¹ and 1 × 10⁻⁷ cm⁻² s⁻¹ (the blue data points in Figure 1), and then fit a power law (the red line) to those data points:

$$\begin{eqnarray}&&{\sigma }_{{\rm{\Gamma }}}(S)=A\,{S}^{-\gamma }.\end{eqnarray} \tag{ 11 }$$

The best-fit values are A = 2.82 × 10⁻⁶ 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). This inclusion of the experimental resolution is a key ingredient to obtaining 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:

$$\begin{eqnarray}&&\begin{array}{l}{\left\langle \displaystyle \frac{{\rm{d}}N}{{\rm{d}}S}\right\rangle }_{{ijk}}={\displaystyle \int }_{{z}_{j,\min }}^{{z}_{j,\max }}{\rm{d}}z{\displaystyle \int }_{{{\rm{\Gamma }}}_{k,\min }^{(\mathrm{obs})}}^{{{\rm{\Gamma }}}_{k,\max }^{(\mathrm{obs})}}{\rm{d}}{{\rm{\Gamma }}}^{(\mathrm{obs})}\,\displaystyle \int {\rm{d}}{\rm{\Gamma }}\,\displaystyle \frac{{\rm{d}}V}{{\rm{d}}z}\,\displaystyle \frac{{\rm{d}}{L}_{\gamma }}{{\rm{d}}S}\\ \,\,\times \,\displaystyle \frac{\exp \left(-\tfrac{1}{2}{\left(\tfrac{{\rm{\Gamma }}-{{\rm{\Gamma }}}^{(\mathrm{obs})}}{{\sigma }_{{\rm{\Gamma }}}({S}_{i})}\right)}^{2}\right)}{\sqrt{2\pi }{\sigma }_{{\rm{\Gamma }}}({S}_{i})}\,{\rm{\Phi }}\left({L}_{\gamma }({S}_{i},z),z,{\rm{\Gamma }}\right).\end{array}\end{eqnarray} \tag{ 12 }$$

Zoom In Zoom Out Reset image size

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, as encoded in the ${C}_{{\rm{P}}}^{{ij}}$ and measured in Ackermann et al. (2018). We perform two sets of fits: one on the 4FGL catalog alone and one that 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 functions 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 χ²s:

$$\begin{eqnarray}&&\begin{array}{l}-2\mathrm{log}{ \mathcal L }({{\boldsymbol{\theta }}}_{\mathrm{BLL}},{{\boldsymbol{\theta }}}_{\mathrm{FSRQ}},{{\boldsymbol{\theta }}}_{{\rm{n}}})\\ ={\chi }_{4\mathrm{FGL}}^{2}({{\boldsymbol{\theta }}}_{\mathrm{BLL}},{{\boldsymbol{\theta }}}_{\mathrm{FSRQ}})+{\chi }_{{C}_{{\rm{P}}}}^{2}({{\boldsymbol{\theta }}}_{\mathrm{BLL}},{{\boldsymbol{\theta }}}_{\mathrm{FSRQ}},{{\boldsymbol{\theta }}}_{{\rm{n}}}).\end{array}\end{eqnarray} \tag{ 13 }$$

Here, θ _BLL and θ _FSRQ denote the parameters of the BLL and FSRQ models, while θ _n marks the 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}_{{C}_{{\rm{P}}}}$.

2.3.1. Fit to the 4FGL Catalog Data

We extract the source count distributions of four different source classes from the 4FGL catalog, then test our model against these distributions. In more detail, we use the following four source classes:

• 1.  
  BLL: the sum of the identified and associated BLLs (the sources labeled BLL or bll in the 4FGL catalog);
• 2.  
  FSRQ: the sum of the identified and associated FSRQs (the sources labeled FSRQ or fsrq in the 4FGL catalog);
• 3.  
  BLZ: the sum of the identified and associated blazars (the sources labeled BLL, BCU, FSRQ, bll, bcu, or fsrq in the 4FGL catalog); and
• 4.  
  ALL: all sources in the catalog.

We expect our model to respect the following constraints:

• 1.  
  The sum of the source count distributions 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 χ² is labeled ${\chi }_{\mathrm{ALL}}^{2}$.
• 2.  
  The sum of the source count distributions of our models for BLLs and FSRQs should be at least as large as the sum of the 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, the BLZ provides a lower bound for our models (labeled ${\chi }_{\mathrm{BLZ}}^{2}$).
• 3.  
  Our BLL model has to explain at least the observed source count distribution of BLLs, and thus also provides a lower bound (labeled ${\chi }_{\mathrm{BLL}}^{2}$).
• 4.  
  In analogy to BLL, the FSRQ model also receives a lower bound from the observed source count distribution (labeled ${\chi }_{\mathrm{FSRQ}}^{2}$).

Each of these constraints would give rise to a contribution to the 4FGL χ². However, a naive sum of the χ² from the three lower bounds would lead to double counting, since the BLL and FSRQ sources also appear in the BLZ class. To avoid this, we consider only the most constraining lower bound between the BLZ case and the combination of BLLs and FSRQs, i.e.:

$$\begin{eqnarray}&&{\chi }_{4\mathrm{FGL}}^{2}={\chi }_{\mathrm{ALL}}^{2}+\max \left({\chi }_{\mathrm{BLZ}}^{2},{\chi }_{\mathrm{BLL}}^{2}+{\chi }_{\mathrm{FSRQ}}^{2}\right).\end{eqnarray} \tag{ 14 }$$

The upper bound from ALL is implemented as:

$$\begin{eqnarray}&&{\chi }_{\mathrm{ALL}}^{2}=\sum _{i}{\left[\max \left(\displaystyle \frac{{\left\langle \tfrac{{\rm{d}}N}{{\rm{d}}S}\right\rangle }_{\mathrm{BLZ},i}-{\left(\tfrac{{\rm{d}}N}{{\rm{d}}S}\right)}_{\mathrm{ALL},i}}{{\sigma }_{\mathrm{ALL},i}},0\right)\right]}^{2}.\end{eqnarray} \tag{ 15 }$$

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 the BLLs and FSRQs, 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.

Table 1. The Binning of the 4FGL Fit

Source Class	Variable	Min	Max	Number of Bins	Scaling
ALL	S [cm−2 s−1]	10−10	10−7	10	log
	Γ	1.0	3.5	1	linear
BLZ	S [cm−2 s−1]	10−10	10−7	10	log
	Γ	1.0	3.5	1	linear
	z	0.0	4.0	6	log in (1+z)
BLL	S [cm−2 s−1]	10−10	10−7	10	log
	Γ	1.6	2.4	5	linear
	z	0.0	4.0	6	log in (1+z)
FSRQ	S [cm−2 s−1]	10−10	10−7	10	log
	Γ	2.1	2.9	5	linear
	z	0.0	4.0	6	log in (1+z)

Download table as:  ASCIITypeset image

For the lower bounds, which relate to identified or associated blazars, we would also like to consider the redshift information, which is not directly provided 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 in the 4LAC catalog is incomplete, i.e., not all sources have a redshift measurement. Since we only consider the identified or associated blazars as a lower bound, this does not represent a problem, but it might not be the most constraining option. Hence, we again consider 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

$$\begin{eqnarray}&&{\chi }_{\mathrm{BLZ},{Sz}}^{2}=\sum _{i,j}{\left[\min \left(\displaystyle \frac{{\left\langle \tfrac{{\rm{d}}N}{{\rm{d}}S}\right\rangle }_{\mathrm{BLZ},{ij}}-{\left(\tfrac{{\rm{d}}N}{{\rm{d}}S}\right)}_{\mathrm{BLZ},{ij}}}{{\sigma }_{\mathrm{BLZ},{ij}}},0\right)\right]}^{2}.\end{eqnarray} \tag{ 16 }$$

In the second case, we disregard the redshift information by integrating over the redshift. We define

$$\begin{eqnarray}&&{\chi }_{\mathrm{BLZ},S}^{2}=\sum _{i}{\left[\min \left(\displaystyle \frac{{\left\langle \tfrac{{\rm{d}}N}{{\rm{d}}S}\right\rangle }_{\mathrm{BLZ},i}-{\left(\tfrac{{\rm{d}}N}{{\rm{d}}S}\right)}_{\mathrm{BLZ},i}}{{\sigma }_{\mathrm{BLZ},i}},0\right)\right]}^{2}.\end{eqnarray} \tag{ 17 }$$

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 χ²s, due to double counting, so again we choose the most constraining one:

$$\begin{eqnarray}&&{\chi }_{\mathrm{BLZ}}^{2}=\max \left({\chi }_{\mathrm{BLZ},{Sz}}^{2},{\chi }_{\mathrm{BLZ},S}^{2}\right).\end{eqnarray} \tag{ 18 }$$

As described in Equation (14), we select the lower bounds by comparing ${\chi }_{\mathrm{BLZ}}^{2}$ with the ones from the analysis of the individual BLL and FSRQ source classes. In this latter case, we can use the full information from the catalogs and compare the models with the data from a three-dimensional grid of flux, redshift, and photon spectral index. Again, since the redshift information is incomplete, we define the χ² as the maximum of the two cases:

$$\begin{eqnarray}&&{\chi }_{{\rm{M}}}^{2}=\max \left({\chi }_{{\rm{M}},{Sz}{\rm{\Gamma }}}^{2},{\chi }_{{\rm{M}},S{\rm{\Gamma }}}^{2}\right),\end{eqnarray} \tag{ 19 }$$

where M stands for either BLL or FSRQ. The individual χ²s in the two cases are defined by

$$\begin{eqnarray}&&{\chi }_{{\rm{M}},{Sz}{\rm{\Gamma }}}^{2}=\sum _{i,j,k}{\left[\min \left(\displaystyle \frac{{\left\langle \tfrac{{\rm{d}}N}{{\rm{d}}S}\right\rangle }_{{\rm{M}},{ijk}}-{\left(\tfrac{{\rm{d}}N}{{\rm{d}}S}\right)}_{{\rm{M}},{ijk}}}{{\sigma }_{{\rm{M}},{ijk}}},0\right)\right]}^{2},\end{eqnarray} \tag{ 20 }$$

when redshift information is included, and by

$$\begin{eqnarray}&&{\chi }_{{\rm{M}},S{\rm{\Gamma }}}^{2}=\sum _{i,k}{\left[\min \left(\displaystyle \frac{{\left\langle \tfrac{{\rm{d}}N}{{\rm{d}}S}\right\rangle }_{{\rm{M}},{ik}}-{\left(\tfrac{{\rm{d}}N}{{\rm{d}}S}\right)}_{{\rm{M}},{ik}}}{{\sigma }_{{\rm{M}},{ik}}},0\right)\right]}^{2}\end{eqnarray} \tag{ 21 }$$

in the case with the flux and spectral index bins only.

To avoid a bias from the Galactic plane, we exclude small latitudes with ∣b∣ < 30° from our analysis.

2.3.2. Fit to the C_P Data

The fit of the APS is performed on the autocorrelation (i = j) and cross-correlation (i ≠ j) measurements, where i and j denote the energy bins. The ${\chi }_{{C}_{{\rm{P}}}}^{2}$ is defined as

$$\begin{eqnarray}&&{\chi }_{{C}_{{\rm{P}}}}^{2}=\sum _{i\leqslant j}\displaystyle \frac{{\left[{\left({C}_{{\rm{P}}}^{{ij}}\right)}_{\mathrm{meas}}-{\left({C}_{{\rm{P}}}^{{ij}}\right)}_{\mathrm{th}}\right]}^{2}}{{\sigma }_{{C}_{{\rm{P}}}^{{ij}}}^{2}}.\end{eqnarray} \tag{ 22 }$$

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, ${\sigma }_{{C}_{{\rm{P}}}^{{ij}}}^{2}$ are the uncertainties of the measured C_P, again taken from Ackermann et al. (2018).

As anticipated above, we perform two fits in this work. 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 Section 2.5, this extrapolation agrees well with the actual C_P measurement. As such, we can go one step further and also perform the second fit, which combines the resolved point sources (4FGL) with the C_P. This combined fit provides GLF and SED models that are consistent with the gamma-ray observations below the flux threshold of the 4FGL catalog. When we derive the 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.

As a result of our fits, we obtain the GLFs and SEDs of the 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 the 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 for correctly assessing the uncertainty of the full blazar model, as, for example, in Section 3, where blazars constitute 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 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 models are well compatible with the parameter values of Ajello et al. (2012) and Ajello et al. (2014), respectively. This means that the GLF of the sources belonging to the fainter regime probed in this work closely follow the one of 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}^{* }$.

Table 2. The Fit Results of the GLFs

	4FGL fit		4FGL+CP fit
	FSRQ	BLL	FSRQ	BLL
${\mathrm{log}}_{10}\left(A\,[{\mathrm{Mpc}}^{-3}]\right)$	$-{9.35}_{-0.37}^{+0.62}$	$-{9.80}_{-1.11}^{+0.40}$	$-{9.65}_{-0.47}^{+0.61}$	$-{9.01}_{-1.11}^{+1.21}$
${\mathrm{log}}_{10}\left({L}^{* }\,[\mathrm{erg}/{\rm{s}}]\right)$	${48.36}_{-0.66}^{+0.31}$	${47.85}_{-0.52}^{+0.79}$	${48.53}_{-0.60}^{+0.42}$	${47.26}_{-1.09}^{+0.75}$
γ1	${0.57}_{-0.09}^{+0.15}$	${1.03}_{-0.07}^{+0.12}$	${0.72}_{-0.09}^{+0.13}$	${0.92}_{-0.07}^{+0.18}$
γ2	${1.93}_{-0.43}^{+0.14}$	${1.95}_{-0.45}^{+0.16}$	${1.97}_{-0.42}^{+0.21}$	${1.88}_{-0.38}^{+0.12}$
${z}_{c}^{* }$	${0.93}_{-0.27}^{+0.20}$	${1.05}_{-0.53}^{+0.16}$	${0.87}_{-0.24}^{+0.14}$	${1.06}_{-0.56}^{+0.16}$
${p}_{1}^{* }$	${5.86}_{-5.02}^{+2.15}$	${7.48}_{-4.51}^{+3.97}$	${8.37}_{-3.35}^{+3.78}$	${4.01}_{-3.54}^{+0.77}$
${p}_{2}^{* }$	$-{0.88}_{-0.14}^{+0.77}$	$-{1.97}_{-0.43}^{+1.83}$	$-{0.77}_{-0.11}^{+0.67}$	$-{0.93}_{-0.19}^{+0.83}$
α	${0.20}_{-0.16}^{+0.07}$	${0.28}_{-0.13}^{+0.14}$	${0.11}_{-0.11}^{+0.02}$	${0.31}_{-0.07}^{+0.17}$
μ*	${2.50}_{-0.04}^{+0.03}$	${2.03}_{-0.04}^{+0.04}$	${2.49}_{-0.04}^{+0.03}$	${2.05}_{-0.04}^{+0.03}$
σ	${0.19}_{-0.04}^{+0.02}$	${0.22}_{-0.05}^{+0.03}$	${0.20}_{-0.04}^{+0.02}$	${0.19}_{-0.03}^{+0.02}$
β		${0.06}_{-0.05}^{+0.02}$		${0.03}_{-0.03}^{+0.01}$
${k}_{{C}_{{\rm{P}}}}$			${1.09}_{-0.12}^{+0.16}$
${\chi }_{\mathrm{ALL}}^{2}$	1.7		2.9
${\chi }_{\mathrm{BLZ}}^{2}$	3.0		2.7
${\chi }_{\mathrm{BLL}}^{2}$	11.4		12.7
${\chi }_{\mathrm{FSRQ}}^{2}$	7.4		8.0
${\chi }_{4\mathrm{FGL}}^{2}$	20.5		26.3
${\chi }_{{C}_{{\rm{P}}}}^{2}$			81.0
χ2	20.5		107.3

Download table as:  ASCIITypeset image

Figure 2 shows the best fits and uncertainties of the dN/dS from the combined fit to 4FGL+C_P in comparison to the dN/dS extracted from the 4FGL catalog. The four different panels correspond to the different contributions to ${\chi }_{4\mathrm{FGL}}^{2}$. 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°). The sum of our models for the BLLs and FSRQs is in agreement with those data. Because of the statistical technique that 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 the 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 the 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 the same 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 4FGL-only fit.

Zoom In Zoom Out Reset image size

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 autocorrelation, 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 the BLL and FSRQ models provides a good fit to the data. This was also the case in the analysis of unresolved sources after 2 yr of Fermi-LAT data (Di Mauro et al. 2014). In the meantime, however, 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): for energies below 200 GeV, bright point sources are masked from the 4FGL catalog, but in the last two bins at high energies, bright point sources are masked using the 3FHL catalog. This leads to a change in 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 (so when an energy bin below 200 GeV is present, the mask is mostly provided by the point sources of the 4FGL catalog). Note also from Table 2 that the nuisance parameter ${k}_{{C}_{{\rm{P}}}}$, introduced to allow for a possible rescaling of the flux threshold sensitivity from the reference model, is within uncertainties compatible with the default value of 1.

Zoom In Zoom Out Reset image size

It is also interesting to look at the 4FGL-only setup, and to 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 of 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 FSRQs.

Zoom In Zoom Out Reset image size

Figure 4 already hints that both BLLs and FSRQs are required to describe the C_P data. We better quantify this statement in the following. Using the results of the 4FGL-only fit, we calculate the posterior distribution of the ${\chi }_{{C}_{{\rm{P}}}}^{2}$, assuming (i) the sum of the FSRQs and BLLs, labeled all, and (ii) only the BLLs. An FSRQ-only 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 Figure 5 show that the sum of the FSRQs and BLLs is preferred. Using the full posterior of the 4FGL-only fit, and defining ${{ \mathcal L }}_{{C}_{{\rm{p}}}}=\exp (-{\chi }_{{C}_{{\rm{p}}}}^{2}/2)$, we calculate the Bayes factors of hypotheses (i) and (ii), obtaining 4.0 × 10³. More details about the calculation of the Bayes factors are given in Appendix B. In this sense, the two physical populations (BLLs and FSRQs) are clearly preferred over a single population of only BLLs or only FSRQs. We note, however, that the preference for these two populations is based on the catalog prior. As discussed in Appendix A, the C_P data are 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.

Zoom In Zoom Out Reset image size

Finally, Figure 6 shows a triangle plot with posterior distributions for the 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 the BLLs is μ* ∼ 2.0, with a width of σ ∼ 0.2. As expected, the FSRQs follow a softer energy spectrum, with an average index of μ* ∼ 2.5 and a similar width. Also, the shape of the GLF at small L is reasonably constrained. The index γ₁ lies between 0.4 and 1.0 for the FSRQs and between 0.5 and 1.2 for the BLLs, while the behavior at large L (see γ₂) is less constrained. We note the degeneracy between A and L*. This is because, in both cases, at first order, their main impact on the fit is to change the normalization of the GLF. The redshift dependence is only weakly constrained due to degeneracies with other parameters.

Zoom In Zoom Out Reset image size

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 the BLLs and FSRQs show some level of mutual correlation in the fit. It is only for the sake of clarity that we do not show the entire triangle plot in Figure 6.

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 mAGNs and SFGs) cannot account for the C_P data. Thus, FSRQs and BLLs provide an inescapable contribution to the UGRB intensity. In Figure 7, we compare the UGRB measurement of Fermi-LAT from Ackermann et al. (2015) with the prediction from our models. The measurement accounted for the 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%.

Zoom In Zoom Out Reset image size

The APS of the cross-correlation between a source field X in the energy bin i and a source field Y in the energy bin j reads (Fornengo & Regis 2014)

$$\begin{eqnarray}&&{C}_{{\ell }}^{{X}_{i}{Y}_{j}}=\int \displaystyle \frac{{\rm{d}}\chi }{{\chi }^{2}}{W}_{i}^{X}(\chi ){W}_{j}^{Y}(\chi ){P}^{{X}_{i}{Y}_{j}}\left(k=\displaystyle \frac{{\ell }}{\chi },\chi \right),\end{eqnarray} \tag{ 23 }$$

where W^X (χ) is the window function of the field X, P^XY (k, χ) is the three-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)

$$\begin{eqnarray}\begin{array}{rcl}{W}_{\mathrm{DM}}(z,E) & = & \displaystyle \frac{{\left({{\rm{\Omega }}}_{\mathrm{DM}}{\rho }_{c}\right)}^{2}}{4\pi }\displaystyle \frac{\langle {\sigma }_{\mathrm{ann}}v\rangle }{2{m}_{\mathrm{DM}}^{2}}{\left(1+z\right)}^{3}{{\rm{\Delta }}}^{2}(z)\\ & & \times \displaystyle \frac{{\rm{d}}{N}_{\mathrm{ann}}}{{\rm{d}}E}\left[E(1+z)\right]{e}^{-\tau \left[z,E(1+z)\right]},\end{array}\end{eqnarray} \tag{ 24 }$$

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, Δ²(z) is the clumping factor, and 〈σ_ann v〉 denotes the velocity-averaged annihilation cross section of DM particles, assumed here to be the same in all DM 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 the gamma-ray photons, which we model as in Finke et al. (2010).

To determine the DM autocorrelation, we compute the three-dimensional power spectrum with the so-called halo model approach (see, e.g., the review in Cooray & Sheth 2002). For X = Y = DM, we have

$$\begin{eqnarray}\begin{array}{rcl}{P}_{\mathrm{DM},\mathrm{DM}}^{1{\rm{h}}}(k,z) & = & {\displaystyle \int }_{{M}_{\min }}^{{M}_{\max }}{\rm{d}}M\ \displaystyle \frac{{\rm{d}}{n}_{{\rm{h}}}}{{\rm{d}}M}{\left(\displaystyle \frac{{\hat{u}}_{\kappa }(k| M,z)}{{{\rm{\Delta }}}^{2}(z)}\right)}^{2}\\ {P}_{\mathrm{DM},\mathrm{DM}}^{2{\rm{h}}}(k,z) & = & {\left[{\displaystyle \int }_{{M}_{\min }}^{{M}_{\max }}{\rm{d}}M\,\displaystyle \frac{{\rm{d}}{n}_{{\rm{h}}}}{{\rm{d}}M}\,{b}_{{\rm{h}}}\,\displaystyle \frac{{\hat{u}}_{\kappa }(k| M,z)}{{{\rm{\Delta }}}^{2}(z)}\right]}^{2}\\ & & \times {P}^{\mathrm{lin}}(k,z),\end{array}\end{eqnarray} \tag{ 25 }$$

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 ${\hat{u}}_{\mathrm{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 Navarro–Frenk–White (NFW) DM density profile (Navarro et al. 1996). All the ingredients in Equations (24) and (25) are modeled as in Ammazzalorso et al. (2020). The minimal and maximal halo masses are set at ${M}_{\min }={10}^{-6}\,{M}_{\odot }$ and ${M}_{\min }={10}^{18}\,{M}_{\odot }$, respectively.

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 as to what concerns the extrapolation of the concentration at low masses, leading 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 the 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 the blazars. The blazar window function can be phrased as W_BLA(z, E) = χ(z)² 〈f_S〉, with the mean flux defined as

$$\begin{eqnarray}\begin{array}{rcl}\langle {f}_{{\rm{S}}}\rangle & = & {\displaystyle \int }_{1}^{3.5}{\rm{d}}{\rm{\Gamma }}{\displaystyle \int }_{{L}_{\min }}^{{L}_{\max }}{\rm{d}}{L}_{\gamma }\,{{\rm{\Phi }}}_{{\rm{S}}}({L}_{\gamma },z,{\rm{\Gamma }})\\ & & \times S\left({L}_{\gamma },z,{\rm{\Gamma }}\right)(1-{\rm{\Omega }}).\end{array}\end{eqnarray} \tag{ 26 }$$

The three-dimensional power spectrum of the cross-correlation between annihilating DM and the blazars is given by

$$\begin{eqnarray}&&\begin{array}{l}{P}_{\mathrm{BLA},\mathrm{DM}}^{1{\rm{h}}}(k,z)={\displaystyle \int }_{1}^{3.5}{\rm{d}}{\rm{\Gamma }}{\displaystyle \int }_{{L}_{\min }}^{{L}_{\max }}{\rm{d}}{L}_{\gamma }\,\displaystyle \frac{{\rm{\Phi }}({L}_{\gamma },z,{\rm{\Gamma }})}{\langle {f}_{{\rm{S}}}\rangle }\\ \,\,\times \,S\left({L}_{\gamma },z,{\rm{\Gamma }}\right)(1-{\rm{\Omega }})\displaystyle \frac{{\hat{u}}_{\kappa }(k| M({L}_{\gamma },z)}{{{\rm{\Delta }}}^{2}(z)},\end{array}\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\begin{array}{l}{P}_{\mathrm{BLA},\mathrm{DM}}^{2{\rm{h}}}(k,z)={\displaystyle \int }_{{M}_{\min }}^{{M}_{\max }}{\rm{d}}M\,\displaystyle \frac{{\rm{d}}n}{{\rm{d}}M}{b}_{{\rm{h}}}\displaystyle \frac{{\hat{u}}_{\kappa }(k| M,z)}{{{\rm{\Delta }}}^{2}(z)}\\ \,\,\times \,{\displaystyle \int }_{1}^{3.5}{\rm{d}}{\rm{\Gamma }}{\displaystyle \int }_{{L}_{\min }}^{{L}_{\max }}{\rm{d}}{L}_{\gamma }\,{b}_{{\rm{S}}}\displaystyle \frac{{\rm{\Phi }}({L}_{\gamma },z,{\rm{\Gamma }})}{\langle {f}_{{\rm{S}}}\rangle }S\left({L}_{\gamma },z,{\rm{\Gamma }}\right)(1-{\rm{\Omega }})\\ \,\,\times \,{P}^{\mathrm{lin}}(k,z),\end{array}\end{eqnarray} \tag{ 28 }$$

where b_S is the bias of the blazars with respect to the matter density, for which we adopt b_S(L_γ , z) = b_h[M(L_γ , z)]. 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).

For the modeling of the signal expected from the Galactic subhalos, we follow 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 originating from 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 antibiased subhalo distribution, corresponding to fiducial model A1 of Ando (2009), with a boost factor for subhalos set to unity.

The subhalo number density as a function of the distance f from the center of the galaxy reads

$$\begin{eqnarray}\begin{array}{rcl}{n}_{\mathrm{sh}}(r) & = & f\displaystyle \frac{{M}_{\mathrm{vir},\mathrm{MW}}}{2\pi {r}_{-2}^{3}{M}_{\min }}\gamma {\left(\displaystyle \frac{3}{{\alpha }_{E}},\displaystyle \frac{2\,{c}_{-2}}{{\alpha }_{E}}\right)}^{-1}{\left(\displaystyle \frac{2}{{\alpha }_{E}}\right)}^{3/{\alpha }_{E}-1}\\ & & \times \exp \left[-\displaystyle \frac{2}{{\alpha }_{E}}{\left(\displaystyle \frac{r}{{r}_{-2}}\right)}^{{\alpha }_{E}}\right],\end{array}\end{eqnarray} \tag{ 29 }$$

where M_vir,MW is the Milky Way virial mass, α_E = 0.68, γ is the lower incomplete gamma function, r₋₂ = 0.81 r_200,MW, c₋₂ = r_vir/r₋₂, 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}_{\odot }$. 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

$$\begin{eqnarray}&&{\overline{n}}_{\mathrm{sh}}=\displaystyle \frac{1}{2}{\int }_{-1}^{1}{\rm{d}}\,\cos (\psi )\,{n}_{\mathrm{sh}}\left(r\left(\cos \psi \right)\right),\end{eqnarray} \tag{ 30 }$$

where $r\left(\cos (\psi )\right)=\sqrt{{r}_{0}^{2}+{s}^{2}-2\,{r}_{\oplus }\,s\,\cos (\psi )}$, 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 direction $\hat{n}$ 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

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\overline{n}}_{\mathrm{sh}}}{{\rm{d}}M}={\overline{n}}_{\mathrm{sh}}\left(r\right)\displaystyle \frac{{\alpha }_{0}-1}{{M}_{\min }}{\left(\displaystyle \frac{M}{{M}_{\min }}\right)}^{-{\alpha }_{0}},\end{eqnarray} \tag{ 31 }$$

where α₀ = 1.9.

The subhalo luminosity L for a subhalo with 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 with the channel under study. It reads

$$\begin{eqnarray}&&L=\displaystyle \frac{\langle {\sigma }_{\mathrm{ann}}v\rangle {M}^{2}}{24\pi {m}_{\chi }^{2}{r}_{s}^{3}}\int {\rm{d}}E\,\displaystyle \frac{{\rm{d}}{N}_{\mathrm{ann}}}{{\rm{d}}E}.\end{eqnarray} \tag{ 32 }$$

With the above ingredients, the APS for the Galactic subhalo contribution can then be expressed as

$$\begin{eqnarray}\begin{array}{rcl}{C}_{{\ell }}^{{ij}} & = & \displaystyle \frac{1}{16{\pi }^{2}{f}_{\mathrm{sky}}}{\displaystyle \int }_{{M}_{\min }}^{{M}_{\max }}{\rm{d}}M{\displaystyle \int }_{{s}_{\min }}^{{s}_{\max }}{\rm{d}}s\displaystyle \frac{{L}^{i}(M){L}^{j}(M)}{{s}^{2}}\\ & & \times \displaystyle \frac{{\rm{d}}{\overline{n}}_{\mathrm{sh}}}{{\rm{d}}M}\,{\tilde{u}}_{\mathrm{sh}}^{2}\left(\displaystyle \frac{{\ell }}{s},M\right),\end{array}\end{eqnarray} \tag{ 33 }$$

where i and j refer to energy bins. The integral along the line of sight is performed up to ${s}_{\max }={r}_{\mathrm{vir},\mathrm{MW}}=258\,\mathrm{kpc}$, and starts at ${s}_{\min }=\sqrt{L/(4\pi {S}_{\mathrm{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⁻¹⁰ cm⁻² s⁻¹ for the first 10 energy bins (the 4GFL threshold) and 2 × 10⁻¹⁰ cm⁻² s⁻¹ for the two highest energy bins (the 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 ${\tilde{u}}_{\mathrm{sh}}$ denotes its Fourier transform. The maximal halo mass for the subhalos in our galaxy M_max = 10¹⁰ M_⊙.

In summary, the APS involving DM is given by the sum of four terms—three extragalactic terms (the autocorrelation from the extragalactic DM halos and the cross-correlations with the BLLs and FSRQs) and the Galactic term, 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).

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 blazar-only case, on top of which the DM contribution is added.

The naive approach to considering the full uncertainty would be to perform an extended parameter scan, which would include the blazar parameters and the DM parameters at the same time. 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 by-product of the MultiNest scan is a set of parameter vectors that follows the multidimensional posterior distribution. This set is provided in the so-called equal-weights sample. We will apply importance sampling to obtain the posterior distribution of the full parameter space (blazars and DM) by using the equal-weights set of the fit to blazars only.

First, we note that we can approximate the integral of the product of the 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:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \int {\rm{d}}{{\boldsymbol{\theta }}}_{{\rm{b}}{lz}}\displaystyle \frac{{{ \mathcal L }}_{0}({{\boldsymbol{\theta }}}_{{\rm{b}}{lz}})\,{p}_{0}({{\boldsymbol{\theta }}}_{{\rm{b}}{lz}})}{{Z}_{0}}\,f({{\boldsymbol{\theta }}}_{{\rm{b}}{lz}})\\ \,\,\approx \,\displaystyle \frac{1}{N}\sum _{i=1}^{N}f({{\boldsymbol{\theta }}}_{\mathrm{blz},i}),\end{array}\end{eqnarray} \tag{ 34 }$$

where ${{ \mathcal L }}_{0}$ is the likelihood, Z₀ is the evidence, and p₀ is the prior. The subscript 0 indicates quantities that refer to the fit without DM. We note that the factor ${{ \mathcal L }}_{0}({{\boldsymbol{\theta }}}_{{\rm{b}}{lz}}){p}_{0}({{\boldsymbol{\theta }}}_{{\rm{b}}{lz}})/{Z}_{0}$ is by definition the posterior distribution. Furthermore, we know that the integral over the prior p₀( θ _blz ) is normalized to 1. If we set $f({{\boldsymbol{\theta }}}_{{\rm{b}}{lz}})={Z}_{0}/{{ \mathcal L }}_{0}({{\boldsymbol{\theta }}}_{{\rm{b}}{lz}})$, we see that the evidence is given by

$$\begin{eqnarray}&&{Z}_{0}=\displaystyle \frac{N}{\displaystyle \sum _{i=1}^{N}1/{{ \mathcal L }}_{0}({{\boldsymbol{\theta }}}_{\mathrm{blz},i})}.\end{eqnarray} \tag{ 35 }$$

We can obtain the marginalized likelihood by integrating over the background parameters:

$$\begin{eqnarray}&&\bar{{ \mathcal L }}({{\boldsymbol{\theta }}}_{\mathrm{DM}})=\int d{{\boldsymbol{\theta }}}_{{\rm{b}}{lz}}\,{ \mathcal L }({{\boldsymbol{\theta }}}_{{\rm{b}}{lz}},{{\boldsymbol{\theta }}}_{\mathrm{DM}})\,p({{\boldsymbol{\theta }}}_{{\rm{b}}{lz}}).\end{eqnarray} \tag{ 36 }$$

Here, θ _DM = {m_DM, 〈σ_ann v〉} denotes the DM parameters, ${ \mathcal L }({{\boldsymbol{\theta }}}_{{\rm{b}}{lz}},{{\boldsymbol{\theta }}}_{\mathrm{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₀( θ _blz )), we see that the integral of Equation (36) is approximately given by the sums

$$\begin{eqnarray}&&\bar{{ \mathcal L }}({{\boldsymbol{\theta }}}_{\mathrm{DM}})\approx \displaystyle \frac{{\sum }_{i=1}^{N}\tfrac{{ \mathcal L }({{\boldsymbol{\theta }}}_{\mathrm{blz},i},{{\boldsymbol{\theta }}}_{\mathrm{DM}})}{{{ \mathcal L }}_{0}({{\boldsymbol{\theta }}}_{\mathrm{blz},i})}}{{\sum }_{i=1}^{N}\tfrac{1}{{{ \mathcal L }}_{0}({{\boldsymbol{\theta }}}_{\mathrm{blz},i})}}.\end{eqnarray} \tag{ 37 }$$

In the next step, we can turn this equation into an expression for a marginalized χ² by using the definition of our likelihood (${ \mathcal L }=\exp (-{\chi }^{2}/2)$):

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Delta }}{\bar{\chi }}^{2}({{\boldsymbol{\theta }}}_{\mathrm{DM}})={\bar{\chi }}^{2}({{\boldsymbol{\theta }}}_{\mathrm{DM}})-{\bar{\chi }}_{0}^{2}\\ \,\,=\,-2\mathrm{log}\displaystyle \frac{{\sum }_{i=1}^{N}\exp \left(-\tfrac{{\chi }^{2}({{\boldsymbol{\theta }}}_{\mathrm{blz},i},{{\boldsymbol{\theta }}}_{\mathrm{DM}})-{\chi }_{0}^{2}({{\boldsymbol{\theta }}}_{\mathrm{blz},i})}{2}\right)}{N}.\end{array}\end{eqnarray} \tag{ 38 }$$

Finally, we obtain the DM limit at the 95% confidence level from the requirement ${\rm{\Delta }}{\bar{\chi }}^{2}({{\boldsymbol{\theta }}}_{\mathrm{DM}})\leqslant 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 contribution scales with 〈σ_ann v〉², while the DM × BLZ contribution scales linearly with 〈σ_ann v〉. So we can tabulate the C_P contributions of the background and DM at a reference value of 〈σ_ann v〉 = 3 × 10⁻²⁶ cm³/s. The equal-weights set contains ${ \mathcal 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 10, since the typical number of evaluations required for a full parameter scan of the blazar and DM parameters requires $\,{ \mathcal 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.

We derive the constraints on the annihilation of DM into a pair of $\bar{b}b$ quarks, which serve as an illustrative example. The limits for the other hadronic channels are expected to be at a similar level. The left panel of Figure 8 shows the marginalized Δχ² in the plane of the DM mass and the 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 ${\rm{\Delta }}{\bar{\chi }}^{2}$ is ∼3.5, which corresponds to a local significance of less than 2σ and an even smaller global significance. Consequently, we can derive the DM limits as a function of the DM mass. In the fiducial setup, namely, using the "LOW" model for the concentration–mass relation of the DM halos (see Section 3), we can place an upper limit of 〈σ_ann v〉 = 10⁻²⁵ cm³ s⁻¹ on the annihilation cross section at the DM mass of 10 GeV. The limit gradually weakens to 〈σ_ann v〉 = 3 × 10⁻²³ cm³ s⁻¹ at 4 TeV. In a more aggressive setting for the concentration parameter, i.e., the "HIGH" model, we obtain a DM limit that 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.

Zoom In Zoom Out Reset image size

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 being linear instead of quadratic (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), but 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 the DM only constitutes a subdominant part. The largest DM contribution stems from the extragalactic DM halos, closely followed by the contribution of Galactic DM subhalos, which is smaller roughly 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 being strongly dominated by the extragalactic term.

Zoom In Zoom Out Reset image size

Finally, we stress again that the absence of a DM signal in the C_P provides further confirmation that the sum of the FSRQs and BLLs fully explains the entire UGRB anisotropy, and that no additional weak component is required to match the data.

The source count distribution as a function of the photon flux S is modeled as a power law:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}N}{{\rm{d}}S}=A{\left(\displaystyle \frac{S}{{S}_{0}}\right)}^{-\gamma },\end{eqnarray} \tag{ A1 }$$

where A is an overall normalization, γ is the power-law index, and S₀ is a reference flux. In the following, S₀ is fixed to 1 × 10⁻¹⁰ cm² s⁻¹, and the fluxes S and S₀ always refer to the photon flux in the energy bin from 1 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:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}N}{{\rm{d}}E}=K{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{-{\rm{\Gamma }}}\,\exp \left(-\displaystyle \frac{E}{{E}_{c}}\right).\end{eqnarray} \tag{ A2 }$$

Here, K is the normalization, E₀ is a reference energy, Γ is the photon spectral index, and E_c is the energy of the exponential cutoff. The exponential cutoff allows the mimicking of the attenuation of gamma rays at high energies. By definition, the flux S is given by $S={\int }_{1\,\mathrm{GeV}}^{100\,\mathrm{GeV}}{\rm{d}}E\,{\rm{d}}N/{\rm{d}}E$. So the ratio of the fluxes S_i and S is given by:

$$\begin{eqnarray}&&{s}_{i}=\displaystyle \frac{{S}_{i}}{S}=\displaystyle \frac{{\int }_{{E}_{\min ,i}}^{{E}_{\max ,i}}{\rm{d}}E\,{E}^{-{\rm{\Gamma }}}\,\exp (-E/{E}_{c})}{{\int }_{1\,\mathrm{GeV}}^{100\,\mathrm{GeV}}{\rm{d}}E\,{E}^{-{\rm{\Gamma }}}\,\exp (-E/{E}_{c})}.\end{eqnarray} \tag{ A3 }$$

Finally, we allow for an intrinsic distribution of the photon spectral indices, following a Gaussian with mean μ and width σ:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}N}{{\rm{d}}{\rm{\Gamma }}}=\displaystyle \frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\displaystyle \frac{{\left(\mu -{\rm{\Gamma }}\right)}^{2}}{2{\sigma }^{2}}\right).\end{eqnarray} \tag{ A4 }$$

Then the C_P of the unresolved point sources is given by

$$\begin{eqnarray}&&{C}_{{\rm{P}}}^{{ij}}=\int {\rm{d}}S\,{\rm{d}}{\rm{\Gamma }}\,{S}_{i}{S}_{j}\displaystyle \frac{{{\rm{d}}}^{2}N}{{\rm{d}}S{\rm{d}}{\rm{\Gamma }}}\left[1-{\rm{\Omega }}(S,{\rm{\Gamma }})\right],\end{eqnarray} \tag{ A5 }$$

where Ω(S, Γ) is the detector efficiency for resolving point sources, which we model as a θ-function at a Γ-dependent flux threshold, S_thr (see the main text). By considering photon spectral indices between 1 and 3, the C_P is then calculated as

$$\begin{eqnarray}\begin{array}{rcl}{C}_{{\rm{P}}}^{{ij}} & = & {\displaystyle \int }_{1}^{3}{\rm{d}}{\rm{\Gamma }}{\displaystyle \int }_{0}^{{S}_{\mathrm{thr}}({\rm{\Gamma }})}{\rm{d}}S\,\displaystyle \frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\displaystyle \frac{{\left(\mu -{\rm{\Gamma }}\right)}^{2}}{2{\sigma }^{2}}\right)\\ & & \times {S}^{2}\,{s}_{i}{s}_{j}A{\left(\displaystyle \frac{S}{{S}_{0}}\right)}^{-\gamma }\\ & = & {\displaystyle \int }_{1}^{3}{\rm{d}}{\rm{\Gamma }}\displaystyle \frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\displaystyle \frac{{\left(\mu -{\rm{\Gamma }}\right)}^{2}}{2{\sigma }^{2}}\right)\\ & & \times {s}_{i}{s}_{j}\displaystyle \frac{A\,{S}_{\mathrm{thr}}^{3}({\rm{\Gamma }})}{3-\gamma }{\left(\displaystyle \frac{{S}_{\mathrm{thr}}({\rm{\Gamma }})}{{S}_{0}}\right)}^{-\gamma }.\end{array}\end{eqnarray} \tag{ A6 }$$

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.

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. The first analysis then employs a single source population, with the SED and GLF as specified in the previous section. 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}_{{C}_{{\rm{P}}}}$) that varies the value of the flux threshold (see the main text for more details). This fit has four 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 parameter, this amounts to six 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. The uncertainties are stated in the Bayesian statistical framework.

The results of our four fits are presented in Table 3 and Figure 10. The χ²/degrees of freedom 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 χ²s of the first and second fits, i.e., the fits without the Γ distribution, we see that including a second population reduces the total χ² by 10.3, which formally corresponds to a statistical shift by 2.8σ. So this could be interpreted as a small hint of two populations. However, if we look at the results with a Γ distribution, the conclusion changes. Here, the χ² 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 either the green or blue contours (population a or b of the two-population fit). Furthermore, both populations, a and b, are required to have a relatively high normalization of A > 10⁻¹¹ cm⁻² s⁻¹ sr⁻¹. On the other hand, if we allow for a Γ distribution, the red contours are fully included in both the green and blue contours. And, in the case of two populations, one of the normalizations A can be pushed to negligible values of 10⁻¹² cm⁻² s⁻¹ sr⁻¹. We compare the C_P energy autocorrelations of the best fit to the data in Figure 11. It is clearly visible that, in the case without the Γ distribution, two populations provide a better fit, while in the case with the Γ distribution there is almost no difference.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3. Fit Results of the Phenomenological Model

	w/o Γ Distribution			w/ Γ Distribution
	1 pop.	2 pops		1 pop.	2 pops
		pop. a	pop. b		pop. a	pop. b
${\mathrm{log}}_{10}\left(A\,[{\mathrm{cm}}^{-2}{{\rm{s}}}^{-1}\mathrm{sr}]\right)$	${12.18}_{-0.16}^{+0.17}$	${11.87}_{-0.13}^{+0.36}$	${11.88}_{-0.17}^{+0.34}$	${12.26}_{-0.18}^{+0.15}$	${11.71}_{-0.20}^{+0.81}$	${11.52}_{-0.42}^{+0.96}$
μ	${2.11}_{-0.03}^{+0.04}$	${1.86}_{-0.11}^{+0.19}$	${2.50}_{-0.29}^{+0.18}$	${2.18}_{-0.05}^{+0.06}$	${2.02}_{-0.07}^{+0.22}$	${2.44}_{-0.32}^{+0.17}$
σ	⋯	⋯	⋯	${0.34}_{-0.08}^{+0.07}$	${0.30}_{-0.09}^{+0.13}$	${0.30}_{-0.08}^{+0.18}$
${\mathrm{log}}_{10}\left({E}_{c}\,[\mathrm{GeV}]\right)$	${2.19}_{-0.13}^{+0.12}$	${1.98}_{-0.15}^{+0.12}$		${1.88}_{-0.12}^{+0.10}$	${1.88}_{-0.11}^{+0.10}$
${k}_{{C}_{{\rm{P}}}}$	${1.00}_{-0.50}^{+0.29}$	${1.00}_{-0.50}^{+0.29}$		${1.00}_{-0.24}^{+0.96}$	${0.97}_{-0.42}^{+0.38}$
χ2	85.2	74.9		74.7	74.5
Δχ2	10.3			0.2

Download table as:  ASCIITypeset image

Finally, we have a look at the C_P energy correlation matrix (the coefficients are defined as ${C}_{{\rm{P}}}^{{ij}}/\sqrt{{C}_{{\rm{P}}}^{{ii}}{C}_{{\rm{P}}}^{{jj}}}$). By definition, the diagonal coefficients are equal to one. If there was only one source population with a single and fixed photon spectral index, the offdiagonal coefficients would also be equal to 1. Instead, the C_P measurement (Ackermann et al. 2018) showed that the offdiagonal coefficients are ∼0.6, which was interpreted as an indication of two populations. In Figure 12, we also show that a single source population with a distribution of photon spectral indices provides the correct offdiagonal pattern.

Zoom In Zoom Out Reset image size

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 $\sigma ={0.34}_{-0.08}^{+0.07}$. 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, BLLs and FSRQs.