Dark matter searches in the \gamma rays with galaxy catalogues

Dark matter searches in the γ-ray extragalactic background via cross-correlations with galaxy catalogues

Alessandro Cuoco Jun-Qing Xia Marco Regis Enzo Branchini Nicolao Fornengo Matteo Viel Dipartimento di Fisica, Università di Torino, via P. Giuria 1, I–10125 Torino, Italy Istituto Nazionale di Fisica Nucleare, Sezione di Torino, via P. Giuria 1, I–10125 Torino, Italy Department of Astronomy, Beijing Normal University, Beijing 100875, China Key Laboratory of Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of Science, P. O. Box 918-3, Beijing 100049, P. R. China. Dipartimento di Matematica e Fisica, Università degli Studi “Roma Tre”, via della Vasca Navale 84, I-00146 Roma, Italy INFN, Sezione di Roma Tre, via della Vasca Navale 84, I-00146 Roma, Italy INAF Osservatorio Astronomico di Roma, Osservatorio Astronomico di Roma, Monte Porzio Catone, Italy INAF Osservatorio Astronomico di Trieste, Via G. B. Tiepolo 11, I-34141, Trieste, Italy INFN, Sezione di Trieste, via Valerio 2, I-34127, Trieste, Italy
Abstract

We compare the measured angular cross-correlation between the Fermi-LAT -ray sky and catalogues of extra-galactic objects with the expected signal induced by weakly interacting massive particle (WIMP) dark matter (DM). We include a detailed description of the contribution of astrophysical -ray emitters such as blazars, misaligned AGN and star forming galaxies, and perform a global fit to the measured cross-correlation. Five catalogues are considered: SDSS-DR6 quasars, 2MASS galaxies, NVSS radio galaxies, SDSS-DR8 Luminous Red Galaxies and SDSS-DR8 main galaxy sample. To model the cross-correlation signal we use the halo occupation distribution formalism to estimate the number of galaxies of a given catalogue in DM halos and their spatial correlation properties. We discuss uncertainties in the predicted cross-correlation signal arising from the DM clustering and WIMP microscopic properties, which set the DM -ray emission. The use of different catalogues probing objects at different redshifts reduces significantly, though not completely, the degeneracy among the different -ray components. We find that the presence of a significant WIMP DM signal is allowed by the data but not significantly preferred by the fit, although this is mainly due to a degeneracy with the misaligned AGN component. With modest substructure boost, the sensitivity of this method excludes thermal annihilation cross sections at 95% level. for WIMP masses up to few tens of GeV. Constraining the low-redshift properties of astrophysical populations with future data will further improve the sensitivity to DM.

Subject headings:
cosmology: theory – cosmology: observations – cosmology: dark matter – cosmology: large-scale structure of the universe – gamma rays: diffuse backgrounds

1. Introduction

The last few years have seen a tremendous improvement in our understanding of the -ray sky, mostly thanks to the observations performed by the Large Area Telescope (LAT) on board of the Fermi satellite [15]. Among the main issues that have been investigated, an important one is the understanding of the origin of the Isotropic Gamma-Ray Background (IGRB) [3, 32], i.e. the fraction of the extra-galactic -ray background (EGB) that has not been resolved into individual sources. The nature of the extragalactic emission is a recurrent issue which arises each time a new observational window of the electromagnetic spectrum is opened on the Universe. A good example is the quest for the origin of the soft X-ray background, with the important difference that the latter has now been largely resolved (see, e.g., [40]) whereas a significant fraction of the -ray flux is still diffuse, leaving large room for potential new discoveries.

The interest in the IGRB also stems from the consideration that the -ray band is a potential “golden channel” for the indirect detection of particle Dark Matter (DM). In fact, among the conventional astrophysical sources that contribute to the IGRB there is the possibility that a characteristic signal from DM annihilation or decay may also be present. After its first detection and early attempts to shed light on the origin of the IGRB (see e.g. [47, 31, 51, 56, 69, 68, 46, 70]), a significant step forward has recently been possible thanks to Fermi-LAT that is resolving an ever growing number of sources [55, 1, 5], most of which have been identified as blazars, almost equally split into Flat Spectrum Radio Quasars (FSRQs) and BL Lacs sub-classes.

The properties of the resolved sources can be used to extrapolate their contribution to the IGRB [9, 8, 29, 28, 7]. These population studies suggest that unresolved blazars account for only about 20% of the unresolved IGRB integrated above 100 MeV, while they can be the dominant component above few GeV. The remaining IGRB fraction is thought to be contributed by star-forming galaxies (SFGs) and misaligned AGN (mAGN), two types of sources that can contribute 10–50% each to the extragalactic -ray emission [42, 4, 27]. The contribution from additional potential sources like the millisecond pulsars located in our Galaxy at high Galactic latitudes turned out to be small [67, 20]. The contribution from known astrophysical sources to the IGRB has thus significant uncertainties and leaves room for an additional contribution by more exotic sources like DM.

Additional constraints on the origin of the IGRB can be obtained by analyzing the angular correlation properties of its fluctuations. Refs. [4, 26, 2] have confirmed the conclusions derived from the IGRB mean-intensity and source populations studies: a blazar population that contributes, at energy below 10 GeV, about 20% of the unresolved IGRB can account for the whole measured angular power, thus providing an independent confirmation that the IGRB is not dominated by emission from blazars in the low energy part of the spectrum. Constraints on the DM contribution have been derived in [13, 38, 33].

The accuracy in the analysis of the IGRB and its fluctuations is limited by the presence of Galactic foregrounds and bright sources. If incorrectly subtracted they can induce spurious contributions to both the mean IGRB intensity and its anisotropies. An effective way of dealing with this problem and filter out contaminations is to cross-correlate the IGRB with maps of sources (observed in other wavelengths or by other means) that trace the same structures where the actual IGRB sources reside but do not correlate with Galactic foregrounds. Basically, all catalogs of extragalactic objects at any redshift satisfy these conditions. In the framework of the IGRB investigation, the cross-correlation strategy has been proposed in [25, 14] and recently revisited in [11, 12]. The measurement was pioneered in [75] using the first 21 months of Fermi data. In that case, no statistically significant signal was observed. The analysis has then been recently updated using 60-months Fermi maps [76]. This time a significant (more than 3.5 C.L.) cross-correlation signal has been detected. The signal is present on angular scales smaller than 1 in the cross-correlation between the diffuse -ray emission cleaned by the Galactic foregrounds and four types of Large Scale Structure (LSS) tracers: radio galaxies in the NRAO VLA Sky Survey (NVSS) [19], near infra-red selected galaxies in the Two Micron All Sky Survey (2MASS) [43] , optically selected galaxies in the Sloan Digital Sky Survey (SDSS)-DR8 catalog [6] and quasi stellar objects (QSO) in the SDSS-DR6 catalogs (QSO-DR6) [61]. No significant correlation was observed with Luminous Red Galaxies (LRG) also from SDSS. The analysis further confirms that blazars provide a minor contribution to the IGRB as found in the IGRB mean-intensity and source populations studies, while a mixture of SFGs and mAGNs can in principle contribute to the majority of the IGRB.

A promising, possibly more effective in the context of DM, way to apply the cross-correlation technique is to use weak gravitational lensing maps (cosmic shear) instead of catalogs of LSS tracers [21, 35, 34, 22, 66]. This alternative approach, originally proposed in [21], has the advantage of probing directly the matter distribution, therefore avoiding the so-called ‘biasing’ issue, i.e., the fact that the mapping between the spatial distribution of extragalactic sources and that of the underlying mass density field is ill-known, and needs to be modelled. Cross-correlation of -rays with cosmic shear will become available in the next years with the release of cosmic-shear maps from wide area surveys, like, e.g., the Dark Energy Survey (DES) [73] and, in the next decade, by the satellite-based Euclid survey [48]. Finally, a similar technique, based on the cross-correlation of -rays with Cosmic Microwave Background (CMB) lensing maps, has been recently adopted in [34], where an evidence of has been reported providing a further direct evidence of the extragalactic origin of the IGRB, and of a subdominant role of Galactic sources.

In this paper we investigate the implications of the recent measurement of a cross-correlation between -rays and LSS tracers by [76] for both the DM and the main astrophysical contributors to the IGRB. This work builds upon the results obtained by [59] in which we concentrated on the low-redshift 2MASS catalog as a tracer of the LSS in the local Universe, and we have assumed that DM-induced -rays provide the dominant source of cross-correlation for such a low redshift observations. That approach has been motivated by the fact that the DM contribution to the cross-correlation is dominated by -rays emission at low-redshift (see e.g. [35] or Appendix A), which is where 2MASS galaxies mostly reside. In that analysis we found that the observed cross-correlation signal can indeed be explained by a DM emission, while its contribution to the total mean intensity is significantly below the IGRB intensity measured by Fermi. This implies that the cross-correlation technique can be a powerful probe of the particle nature of DM, even when the DM contribution to the IGRB is subdominant, which is what we expect in a realistic scenario. In [59] we found that the cross-correlation signal can be explained by a DM particle with mass in the tens to hundreds GeV range (depending on the -rays production channel) and, once the uncertainties in the DM distribution modeling is properly accounted for, a “thermal” value for the annihilation cross section ( cm s), which is the most appealing case for a weakly interacting massive particle (WIMP) DM. From the same analysis we have obtained upper bounds on the DM annihilation cross section and decay rate that turn out to be quite competitive with those obtained with different techniques, based either on local (Galactic halo, dwarf galaxies) or extra-galactic -ray emission. We point out that those constraints are conservative precisely because the DM is assumed to be the only source of the -ray signal.

In this follow-up paper we extend the study of [59] to the inclusion of astrophysical -ray emitters, and to the whole set of LSS-tracers catalogs. As it will be discussed in the next Sections, the redshift distributions of the -ray signal is a fingerprint that characterises the contribution of different astrophysical sources and of the DM. For this reason, the possibility to use catalogues of objects whose distributions peak at different redshifts is an effective way to extract the information encoded in the -rays maps and remove degeneracies. To this aim, in addition to DM, we account here for contributions from blazars (of both BL Lac and FSRQs types), SFGs and mAGNs. Consequently, we do not limit our cross correlation study to the 2MASS catalog but consider NVSS, SDSS-DR6 QSO, LRGs and SDSS-DR8 Main Galaxy samples as well. The approach will be similar to [76] where, however, only astrophysical sources have been fitted to the observed correlations. Here, beside including DM in the fit, we will use an improved description of the cross-correlation modeling between astrophysical sources and LSS tracers based on the Halo Occupation Distribution (HOD) formalism. As for the DM, we shall use the halo model to trace its spatial distribution and predict its cross-correlation with LSS (see e.g., [24, 13, 35]).

The paper is organized as follows. In Section 2 we present the theoretical estimate of the angular cross-correlation function and angular power spectra. Section 3 describes the statistical techniques employed in the determination of the parameters of the -rays emitters (DM and astrophysical sources) from the measured cross-correlation reported in [76]. Section 4 then shows our results, and finally Section 5 summarizes our conclusions. The technical aspects of our theoretical modeling are presented in a set of three Appendices. Appendix A introduces the modeling of the window functions of DM and astrophysical -rays sources and of catalogs of LSS tracers. Appendix B discusses the HOD of galaxies for the various catalogs. Appendix C describes the derivation of three-dimensional (3D) power spectra (PS). These are the ingredients used in Section 2.

In this work we assume a fiducial flat CDM model with the cosmological parameters derived by the Planck Collaboration in [57]: matter density parameter , baryon density parameter , reduced Hubble constant , rms matter fluctuations in a comoving sphere of 8 Mpc and spectral index of primordial scalar perturbations .

2. Formalism

To quantify the cross-correlation between -ray sources and the LSS tracers in the various catalogues, we consider both the 2-point angular cross-correlation function (CCF) and its Legendre transform, i.e. the cross angular power spectrum (CAPS). In the Limber approximation [50, 44, 45], the CAPS can be obtained by integrating the 3D-PS of cross-correlation :

 C(γg)ℓ=∫dχχ2Wγ(χ)Wg(χ)Pγg(k=ℓ/χ,χ), (1)

where denotes the radial comoving distance, is the so-called window function that characterizes the distribution of objects and -ray emitters along the line of sight, is the modulus of the wavenumber and is the multipole. We relate the cosmological redshift to the radial comoving distances through the differential relation, valid in a flat cosmology, , where is the expansion rate of the Universe.

The indices  and denote -ray emitters and extragalactic objects in different catalogs, respectively. We consider five types of -ray sources: three different flavours of AGNs (BL Lacs, FSQRs, mAGN), SFGs and DM. We will consider both the case of annihilating and decaying DM particles. For the LSS tracers, we consider five different catalogues: quasars in SDSS-DR6, 2MASS galaxies, NVSS radio sources, SDSS-DR8 Luminous Red Galaxies and SDSS-DR8 “main” galaxies.

Denoting the density fields of an LSS tracer with , and that of the gamma ray emitter with , where indicates the position in comoving coordinates and labels time (given the one-to-one correspondence between time and distance), the cross-power spectrum is defined as:

 ⟨^fγ(z,k)^f∗g(z′,k′)⟩=(2π)3δ3(k+k′)Pγg(k,z,z′), (2)

where is the Fourier transform of , indicates the average over the survey volume and the explicit dependence on and highlights the possibility that the two populations under study (-ray emitters and extragalactic LSS tracers) are located at two different redshifts. From the Limber approximation one gets , so, in practice, only is used. The modeling of the various power spectra used in our analysis is derived in Appendix C. Objects in the catalogs are described in terms of their halo occupation distribution (HOD), which is discussed in Appendix B.

The window function appearing in Eq. (1) weights the contribution of objects at different redshifts to the cross-correlation signal. In the case of LSS tracers it coincides with the redshift distribution of the objects, . More precisely, such that for a redshift distribution normalized to unity. The expressions of for the different types of LSS tracers that we consider here are the same as in [76] (see also Appendix A.3).

For a -ray emitter the window function can be defined in term of the -ray intensity integrated along the line-of-sight, as function of the direction in the sky , which can be written as:

 Ifγ(n)=∫dχfγ(χ,r)⟨fγ(χ)⟩Wγ(χ) (3)

so that . We will use a coordinate system centered on the observer so that . The expression of the density fields and window functions for the different classes of -ray sources are provided in Appendix A.

In the Appendices we also define our reference models for the astrophysical and DM -rays emitters, and for their cross-correlations with the LSS tracers. The mean intensity as function of energy of the different -ray emitters for our reference models is shown in the left panel of Fig. 1. The various curves in color indicate the contribution of each component, as indicated by the labels, while the black line indicates the sum of all astrophysical contributions. The predicted total energy spectrum matches the recent Fermi-LAT measurements [3] (solid dots with error bars). Similarly, as shown in the right panel of Fig. 1, we also verified that our reference model matches the observed angular power spectrum of the diffuse extragalactic gamma ray background measured in the 1–2 GeV energy band by the Fermi-LAT (grey strip) [2]. The different curves in color show the predicted angular power spectra of the various emitters that contribute to the total angular spectrum (solid black line). The model angular power spectra for the various gamma ray emitters have been derived by using in Eq. (1) the power spectrum of the source instead of the cross-spectrum , and instead of the product . With respect to the more accurate procedure used in [28], here we use the simplifying assumption that all the sources of a given population have the same photon spectral index (see Appendix A).

In the next Section we will fit the theoretical predictions to the measured cross-correlations and will estimate the free parameters of the models. For the astrophysical components, the results will be in the form of deviations from reference models, which we adopt from the literature as updated benchmarks. We will therefore allow variations only in their normalization, plus a correction term as specified in the next Section.

A technical remark to take into account when comparing the model with the data is that the experimental CAPS determined from the data are not deconvolved from the effect of the point spread function (PSF) of the instrument and the effect of map pixelization. To account for these effects we thus convolve our model prediction in Eq. (1) with the PSF and pixelization using the same procedure described in [76]. Formally, this is implemented defining the new quantity directly comparable with the data as where the effective beam window function parameterizes the PSF and pixelization effects (see [76] for more details).

Finally, in the following, we perform our analyses in terms of the cross correlation function rather than the cross-angular power spectra . To obtain the CCF we perform a Legendre transformation on our CAPS as follows:

 CCF(γg)(θ)=∑ℓ2ℓ+14π~CγgℓPℓ[cos(θ)], (4)

where is the angular separation in the sky and are the Legendre polynomials.

3. Statistical Analysis

In order to assess the possible presence of a DM signal in the IGRB, and its robustness to the presence of astrophysical emitters, we perform a statistical analysis fitting the observed cross-correlation data of [76] with a combination of both DM and astrophysical source models. Specifically, we define a statistic from the data , i.e., the observed CCF between the Fermi maps and the number of sources in catalogues [76], and , i.e. the model CCF calculated for the different types of -ray emitters as introduced in the previous Section and detailed in the Appendices. The is defined as:

 χ2=5∑p=13∑n=1∑θiθj(D(p,n)θi−M(p,n)θi(A))[C(p,n)]−1θiθj(D(p,n)θj−M(p,n)θj(A)), (5)

where the index runs over the five different catalogues of extragalactic sources (2MASS, NVSS, SDSS-DR6 QSO, SDSS-DR8 Main Sample Galaxies and SDSS-DR8 Luminous Red Galaxies), the index runs over three -rays energy ranges ( GeV, GeV and GeV), whereas the indices and run over 10 angular bins logarithmically spaced between and . is the covariance matrix that quantifies the errors on the CCFs in each angular bin and the covariances among different bins, and denotes the vector of free parameters which the CCF model depends upon (specified below). Both the covariance matrix and the measured CCFs are taken from [76]. In Eq. (5) the total is obtained by adding up the individual computed in three overlapping energy bands. There is, thus, in principle, a statistical dependence among the different energy bands that should be accounted for. Nonetheless, such dependence is expected to be small since photon counts are heavily dominated by events near the lower end of each energy interval because of the steep IGRB energy spectrum [3]. For this reason we will treat the CCFs estimated in the three energy intervals as statistically independent in the analysis.

For any given catalog of LSS tracers, energy band and angular bin (i.e. for a given choice of , , and ) the theoretical CCF can be expressed as a sum of different contributions:

 M(p,n)θi=5∑α=1Aαc(p,n)α(θi)+A(p)1hc(n)1h(θi). (6)

The sum runs over the five different -ray emitters: BL Lacs, FSRQs, SFGs, mAGNs and DM. The terms denote the benchmark theoretical model CCFs described in the Appendix and is a free normalization parameter that quantifies the individual contribution to the observed cross-correlation. Values thus denote models equal to our benchmarks, while values would correspond to deviations from the benchmarks.

Besides the normalization of each component, we have introduced in the fit also a further free parameter, , which we dub 1-halo correction-term. This term is introduced as a correction for possible inaccuracies in the modeling of the 1-halo contribution of the power spectrum (and thus mainly to the small scale cross-correlation signal) of the -ray sources (hence its name), as discussed after Eq. (C6) in Appendix C. For simplicity we model it as a constant term added in the CAPS, which is a good approximation of the 1-halo term itself for astrophysical components, except at very high multipoles , which we do not considered in our analysis. In real space, and taking into account the modulation introduced by the PSF of the instrument, the 1-halo correction-term explicitly reads:

 M(p,n)1h(θi)=A(p,n)1h∑ℓ2ℓ+14πWBnℓPℓ[cos(θi)], (7)

where is the (energy dependent) window function of the PSF, also introduced in the previous section. For definiteness we have assumed that has the same energy dependence as the IGRB spectrum . With this assumption we can separate the energy dependence from that on the source catalog where for GeV (we take GeV as the normalization energy). In this way, combining Eq.s (6) and (7) we have: . Notice that, in principle, each candidate -ray emitter has its own 1-halo correction-term for each catalog and energy band. However they are degenerate and only their sum can be constrained. We have thus grouped them together so that in Eq. (6) the reported 1-halo correction-term actually represents the sum of the 1-halo correction-terms from all the components, for a given catalog and energy band.

All five values are treated as free parameters in the analysis. Notice that they can be either positive or negative since we intend them as possible correction to our benchmarks, and the natural expectation would thus be = 0 if the benchmarks are correct. The whole set of and coefficients (plus the additional parameter represented by the DM mass, upon which the DM signal depends) defines the parameter vector of Eq. (5), which represents the full set of parameters over which our analysis is performed.

For what concerns the particle DM contribution, we consider both the case where -rays are produced through DM particle annihilation and the case of DM decay. The DM mass is varied from 10 GeV to 5 TeV and we will show the results for a DM which dominantly annihilates/decays into one of the following -rays production channels: , , and . For annihilating DM, the signal strongly depends on the clustering at small scales and, in particular, on the amount of substructures. As discussed in Appendix A, results will be shown for three DM substructure models: high, low and ns. This will bracket the uncertainty on the reconstructed DM parameters arising from DM structure modeling. The high scenario provides a more optimistic case with the largest boosting factor for the -ray annihilation flux. The ns scheme, where DM substructures are absent, provides a lower limit to to the annihilation signal and therefore represent the most conservative scenario. Finally, the low scheme represents an intermediate case which can be (currently) considered as the most realistic one and that we regard as our reference model. Each one of these three scenarios predicts a different CCF, with in Eq. (6). Since the intensity of the DM signal is proportional to the DM annihilation cross section (or DM decay rate, ), we normalize our calculations to the reference values , i.e. the so-called “thermal value” which correspond to a DM particle thermally produced in the early Universe which, alone, would account for the observed DM relic abundance. For decaying DM we normalize the models to a decay rate of , which is the decay rate which would produce a DM signal equal to the one of an annihilating DM with thermal cross section in the low substructure scheme (for DM masses around 100 GeV). The parameter for can be thus seen as the annihilation or decay rate in units of or , respectively.

In summary, the global fit will be performed in a 11-dimensional parameter space, with the parameter vector given by . All the parameters in the fit are linear, except which enters non-linearly in the fit through . Beside the above fit, we will also consider different configurations, namely different parameter vectors , to cross-check the robustness of the results. In particular, we will consider the case where the are set equal to zero. These additional analyses will be described in more details in the next Section.

In order to efficiently scan the multi-parameter space we adopt the Markov Chain Monte Carlo (MCMC) strategy publicly available in the cosmomc package [49]. We will use linear priors limited to positive values for the normalization of the astrophysical components , although we will also check log priors. For the parameters we allow for linear priors with negative values since the 1-halo correction-term can either correct for over-estimation or under-estimation of the small-scale cross-correlation. Finally we will use a logarithmic prior for and since, theoretically, the possible values of the DM mass and signal normalization can span several orders of magnitude.

Notice that in our analysis we consider only the cross-correlation signal and ignore the intensity and the auto-correlation of the IGRB. These additional observational inputs will be used to perform an independent a posteriori check on the validity of our results. While the total intensity and the -rays autocorrelation calculated from the derived best-fit configurations must not exceed the measured values (this will be a sanity check), if they fall short of accounting for the data this might indicate either that the measured IGRB contains an unaccounted contribution which does not correlate with extragalactic tracers (possibly of Galactic origin) or that the modeling of the known components is imperfect/incomplete. We will discuss more in detail these aspects in section 4.1 and in the conclusions.

4. Results and Discussion

The triangle plot shown in Fig. 2 summarizes the results of our analysis for a benchmark annihilating DM case with final state and low substructure scheme. The plot shows the posterior marginal distributions of the normalization parameters . The two-dimensional plots refer to the and credible regions for each pair of parameters, while the diagonal shows the one-dimensional posterior distribution for each parameter. The parameter space is eleven-dimensional, but for clarity we show only a part of the full triangle plot (without including here the parameters and ). The two-dimensional posterior of (, ) is shown separately in the left panel of Fig. 3. The posterior probability for is instead displayed in the right panel of Fig. 3, while the one-dimensional posteriors for the parameters are shown in Fig. 4. Finally, Fig. 5 shows the best-fit results compared with the measured cross correlation functions.

A noticeable result from Fig. 2 is the fact that all the posteriors seem to peak at or close to it, except (but with a low significance) for the DM and mAGN constributions. This does not necessarily imply that the best fit is found when the contributions from all components is zero. Instead, it is an indication that strong degeneracies are present. A two-dimensional analogy is given by a case with only two parameters related by a simple relation const. While the degeneracy would be clearly seen in the two-dimensional posterior, both the one-dimensional and posteriors peak at zero, although and are never both zero at the same time. This is precisely the results we find here, although the high dimensionality of our parameter space prevents us from clearly trace the parameter degeneracy even in the two-dimensional posteriors plots.

The degeneracy of the different astrophysical components can be traced to the behavior of their respective window functions , which possess a relatively similar evolution as a function of redshift, and to a similar behavior of their cross-correlation 3D power spectra. As can be seen in Fig. 13 in Appendix A, apart from the DM case for which the -rays emission is concentrated at low redshift with a fast decrease for increasing distances, astrophysical sources possess a relatively broad kernel. Fig. 14 and Fig. 16 instead show some examples of 3D cross-spectra between LSS tracers and the various astrophysical sources considered here or DM: we notice that, for a given LSS tracer (e.g. 2MASS in Fig.14), the behaviors are quite similar for all astrophysical sources (while, instead, differences can be appreciated for the DM case). These facts, together with the relatively large error bars makes astrophysical-component separation currently difficult. On the other hand, given the somewhat different 3D cross-spectra, perspective to separate the DM component are, perhaps, brighter.

An exception in this line of reasoning are mAGNs and SFGs which exhibit a significant degree of degeneracy with DM in the 2D posterior contours. The main features of the DM signal is that it peaks at low redshift and that is mostly contributed by massive halos. To mimic such a signal an astrophysical source must then preferentially be hosted in large halos at low . Both SFGs and mAGNs meet the redshift requirement while the blazars do not, since their window peaks at higher . However, only mAGNs are believed to be hosted in large halos, while SFGs typically populate galaxy-size halos. Objects in large halos at low redshifts are expected to have a large bias and, more importantly, their correlation properties at the Mpc scale is dominated by a large one-halo term. This introduces a characteristic feature in the cross-PS that differentiates mAGNs from SFGs, making their contribution more similar to the DM one at Mpc scales (see left panel of Fig. 15). At the lowest redshift considered (namely, in the cross correlation with 2MASS), the Mpc scale corresponds to a sub-degree scale in the CCF. Nonetheless, given the present still large error bars, the above feature is only weakly constrained and thus a further degeneracy of both components with SFGs still remains on top of the mAGN-DM main degeneracy. Further investigation of this issue is reported later below. Instead, further differences between the mAGNs and and the DM cases are expected at smaller angles which, unfortunately, cannot be investigated given the size of the Fermi-LAT PSF.

Difficulties in modeling the 1-halo term in the HOD framework described in Appendix B propagates into uncertainties in predicting the cross-power at small-angles. To account for this potential source of systematic errors we introduced in Eq. (7) the 1-halo correction-terms . The 1D marginalized posteriors of the associated extra five parameters are shown in Fig. 4 as black solid curves. The various datasets are consistent with the case with different confidence levels, except NVSS. In this case we find a strong and statistically significant deviation from zero. This can also be appreciated in the fit to the observed CCF in Fig. 5 where the presence of a prominent 1-halo correction-term is required to fit the data at small angles. There is a likely explanation for this additional contribution: the presence in the NVSS catalog of -ray point sources (i.e., AGN) that are just below Fermi detection threshold. These sources would add their auto-correlation signal at zero-lad that, because of the PSF, spreads out to deg scale. This effect requires some fine tuning of the parameters defining the 1-halo term in Eq. (C6) which the benchmark model fails to catch, thus requiring a large correction term. The effect is also discussed in [76] to which we refer the reader for further discussion. The relevance of this term in the fit to NVSS data is expected to affect our constraints of the DM properties. To investigate this issue we use three further fitting procedures in addition to the one adopted so far. The four fitting procedures are as follows:

• NVSS-10, . All the 10 NVSS data points are fitted and the 1-halo-correction terms are free parameter of the fit. This is the standard fitting procedure used to obtain the results shown in Fig. 2.

• NVSS-10, . All the 10 NVSS data points are fitted and all the 1-halo-correction terms are set equal to zero.

• NVSS-6, . The first 4 NVSS data points at small angles are excluded from the fit. The 1-halo-correction terms are used as free parameters in the fit.

• NVSS-6, . The first 4 NVSS data points are excluded from the fit. All 1-halo-correction terms are set equal to zero.

Fig. 4 shows the posteriors for the NVSS-10, (black solid curves) and the NVSS-6, (red dashed curves) cases. It can be seen that there are no significant differences between the two fitting schemes, except for the NVSS case in which becomes obviously unconstrained when the first four data points, where the fit is guaranteed by the 1-halo-correction term, are ignored. Fig. 6 quantifies the impact of the four fitting schemes on the posterior probabilities of all parameters. The plots show that the fitting procedure does have an impact on some parameter. In particular the results obtained with the NVSS-10, fit deviates from the others in most of the cases. However, this is also the scheme that provides the worst fit to the various datasets, as illustrated by the comparatively larger values listed in Table 1, so that the results from this case are likely somewhat biased. Instead, all the three remaining schemes provide reasonably good fits to the different datasets. The NVSS-6, case provides a slightly worse fit to the data, particularly to the LRG sample, than the other schemes. It is interesting to notice that this fitting scheme favours a non-zero mAGN component (see the panel in Fig. 4) which is absorbed by the 1-halo-correction term when a fitting scheme with is adopted. This indicates that a degeneracy between the and the parameter is present. We notice that in all three cases the is lower than total number of degrees of freedom. This is partly due to some unaccounted correlation between the three energy bins and between the different catalogues, and partly to the fact that the error bars are probably slightly over-estimated. Indeed, it is known that the algorithm implemented in the PolSpice software which is used in [76] is not a minimum variance estimator of the error bars, i.e., it does not provide the smallest error possible [30].

Let us now discuss in more details the implications for the DM component. From the CCF plot in Fig. 5 we see that DM provides a significant contribution to the fit. Yet, the posterior probability of in the bottom right panel of Fig. 2 does not provide a clear indication for a DM component. As discussed above, this is an indication that a DM signal may indeed be there but is degenerate with some other component, in particular the mAGN one. In practice, the present datasets and our cross-correlation analysis cannot distinguish between the case of a large DM contribution with sub-dominant mAGN signal and that of a mAGN signal that dominates over the DM contribution. In fact, the situation is further complicated by the aforementioned degeneracy between mAGN and the 1-halo-correction terms. Before discussing the degeneracy issue more in detail, it is worth pointing out that i) our results are robust to the choice of the fitting strategy, as shown in the right panel of Fig. 3 and in the bottom panel of Fig. 6, and that ii) despite the DM vs. mAGN degeneracy we are able to set constraints on the annihilation cross-section, as shown in the left panel of Fig. 3, able to exclude the thermal value at for DM masses up few tens of GeV (in the low substructure scheme) that, again, are robust against the adopted fitting scheme. We also verified the robustness of the DM constraints with respect to the choice of the priors for the astrophysical components. Specifically, we considered the case of log-flat priors instead of a linear-flat ones, and we found that the posteriors of the DM parameters are unaffected. Some small variations are present in the constraints of the astrophysical parameters, which is expected since at the moment the significance of the measurement is still not very high, and in this regime some prior dependence is typically still present.

One thing to notice about the DM vs. mAGN degeneracy is that few mAGNs have been detected in -rays so far. As a consequence their model contribution to the IGRB and its anisotropies is rather uncertain. One key quantity is the relation between the -ray luminosity of these objects, and the mass of their host halo . Varying this relation within its uncertainty range, which is rather large (see e.g. [22]), modifies the predicted cross-correlation signal. Fig. 15 illustrates this point. In the right panel we show the cross-power spectrum mAGN-2MASS galaxies (solid line) and how it changes when the relation is varied within its uncertainty band (dashed curves). Considering halo masses in the lower bound of the uncertainty strip significantly decreases the amplitude of the 1-halo term and reduces the amplitude of the PS on Mpc scales. As a result the PS contributed by mAGNs will be very similar to that contributed by SFG, as shown in the left panel of Fig. 15, and since their window functions are also very similar (see Fig. 13), their contributions to the cross-power become fully degenerate.

We are therefore entitled to consider a scenario in which, due to this degeneracy, we set the mAGN contribution equal to zero and assume that it is absorbed by the SFG one. To explore this situation we consider four additional fitting schemes. In all of them we ignore the first data-points of the NVSS dataset (i.e. we use the NVSS-6 scheme) and set . The four schemes are obtained from all possible combination of and that are either set equal to zero or let free to vary. The four combinations are explicitly shown in the first column of Table 2 in which we summarize the results of the analysis.

The inclusion of the 1-halo-correction terms improves the fit appreciably, although the improvement is mainly driven by the LRG and QSO datasets. The inclusion of DM with two extra-parameters also improves the fit decreasing the best fit by 4.3 and 5.8 for the fits with and without 1-halo-correction terms, respectively, with improvement mainly coming from a better fit to the 2MASS data. No scheme provides a good fit to the CCF with the QSO for 500 MeV. This is possibly an indication of an imperfect modeling of the energy spectrum in the QSO correlation.

Fig. 7 shows the triangle plot for the case 0 and 0. When the mAGN contribution is suppressed (0) a non-vanishing DM component provides quite a good fit, with about a deviation from zero. Furthermore, the best fit values in Table 2 are very similar to those in Table 1 in which the mAGN component was included in the model (71.6 vs 72.5 and 89.5 vs 85.8 for the case with and without 1-halo-correction terms, respectively), a fact that corroborates the evidence for a degeneracy between DM and mAGNs. Fig. 8 illustrates the robustness of these results to the inclusion of the 1-halo-correction term, 0, for two different final annihilation states: (left set of plots) and and (right plots). In all cases the best fits preference is for the presence of a DM component, even when the extra degree of freedom is included. Furthermore, these plots reveal a degeneracy between and which is to be expected since, as can be seen in Eq. (A1), the WIMP signal approximately scales with . Indeed, contains a factor plus an integral over the energy which introduces a contribution roughly proportional to (being the spectrum integrated up its endpoint, which is ). Inspection of Fig. 7 also shows more clearly the degeneracy between DM and the SFG components, left from the mAGN-DM-SFG degeneracy after removing the mAGN component. One consequence of this is that performing a fit excluding either the SFG or the mAGN component would enhance the strength of the DM signal, without affecting appreciably the value of the best fit .

In conclusion, our analysis indicates that a significant DM contribution to cross-correlation is entirely plausible. However, the degeneracy with other astrophysical sources, namely mAGNs and SFGs, largely originating from the current observational uncertainties, prevents us from drawing a definitive conclusion. Future analyses with increased -ray statistics and improved angular resolution in which the cross correlation is extended to other catalogues of extragalactic objects will help to break the present degeneracies and to pinpoint the correct scenario.

As a final remark, we also mention that, as a cross-check, we have performed the analysis employing the astrophysical models adopted in [76]. The constraints on the -ray astrophysical contributions are different, which is expected given the different modeling. Regarding DM, using the NVSS-6, , fit configuration, which is the closest to the one used in [76], we compare in Fig. 9 the DM posteriors from 3 different fits using 3 different SFG models: the one adopted in this work (black solid curve), the SFG1 model from [76] (red dashed curve) and a modified version of SFG1 (blue dot-dashed curve) with redshift-dependent bias equal to the the bias of the present model (while the original SFG1 model has bias equal to 1 for all ). The SFG2 model of [76] is very similar to present SFG model and is not considered. The plot indeed shows that the DM results are not significantly dependent from the SFG model adopted.

With no unambiguous indication for a DM components we can nevertheless set constraint on the properties of the DM candidates. To this purpose we perform, for any given mass of the DM particle candidate, an individual 10 parameter fit and set the 95% bound on from the posterior distribution. The results for the annihilating DM are summarised in Fig. 10. In the left panel we focus on the annihilation channel. The solid line with different colours refer to constraints obtained from each of the three energy band separately ( GeV) as well as the ones obtained by their combination (red line). All theee results refer to the low substructures model and are obtained with the reference NVSS-10 0 fitting scheme. The upper and lower dot-dashed curves show how the bounds change in the the ns and high substructure model, respectively.

In the right panel we compare the results of different final annihilation channels (, and ) to the original case (red curve) shown in the left plot, for the low substructures scenario and combining all energy bands. All the results refer to the benchmark NVSS-10 0 case, but the other fitting schemes provide nearly indistinguishable constraints. The black curve is taken from [59] and refers to the case in which we assumed that all the 2MASS -ray correlation is produced by DM, with no astrophysical contribution. As expected, including the astrophysical sources makes the constraints stronger, of about a factor of 4. The gain is significant and will further improve once the DM-mAGN-SFG degeneracies discussed above will be removed.

As expected, uncertainties on the bounds driven by the substructure model are significant. The left panel of Fig. 10 shows that assuming the high model would strengthen the constraints on the cross section by about one order of magnitude, whereas in the ns scenario, the bounds would weaken by about a factor of 5. This implies that the thermal annihilation rate is excluded at the 95 % level up to masses of 6, 25, 250 GeV in the ns, low and high scenarios, respectively.

In Fig. 11 we instead show the 95% lower bounds on the lifetime of a decaying DM particle, for various decay final states. Bounds on DM decay, being proportional to the DM density (and not DM density squared, as instead the annihilation signal) depend on the total DM mass in structures and are not affected by the different substructure modeling. As for the annihilation case, including the astrophysical sources in the analysis improves the constraints, again by about a factor of 4, with respect to those obtained by ignoring the astrophysical components [59].

Finally, to test the robustness of our DM constraints we have repeated the analysis using the same astrophysical models used in [76] and we found that they are very similar to the ones obtained in the present analysis.

4.1. Self consistency tests: mean intensity and auto-correlation of the IGRB

As anticipated in Section 3 instead of including the mean IRGB intensity and its auto-correlation in the fit, we use these additional observational inputs a posteriori as a self-consistent test for our best fitting model.

We define, , the fractional mean IGRB intensity predicted by the cross-correlation fit, as follows

 InTOTAnIGRB = AFSRQInFSRQ+ABLLacInBLLac+ (8) AmAGNInmAGN+ASFGInSFG+ADMInDM,

where are the integrated -ray intensities of our reference models for the five -ray emitters considered here and shown in Fig. 1 and identifies the energy band The total intensity is defined as , where the sum runs over the five types of emitters. In our model = for the energy ranges GeV, respectively, which are consistent with the measured IGRB [3]. We thus expect that the have values close to unity to match observations. Note that the parameters need not to be the same in each energy band since the total signal and the individual contributions have different scaling in energy. However, the difference is not very large.

Similarly, we define the IGRB auto-correlation predicted from the cross-correlation fit, as a fraction of the measured one, in terms of the parameter as:

 CP,TOTfCP,IGRB = A2FSRQCP,FSRQ+A2BLLacCP,BLLac+ (9) A2mAGNCP,mAGN+A2SFGCP,SFG,

where , , and , all of them in units of (cmssrsr. are the predicted average auto-correlation signals in the multipole range and in the energy band - GeV. We have neglected the DM contribution since it is largely subdominant with respect to FSRQs, BLLacs and mAGNs (see Fig. 1). The SFG contribution, which is also subdominant, is considered for the sake of completeness. In the above equation we made the assumption that the amplitude of the auto-correlation signal scales with the square of the normalization parameters of the individual components. Unlike the cross-correlation case we did not include any 1-halo-correction term since we model astrophysical emitters as point sources for which no additional small-scale power is expected to contribute to the auto-correlation signal. Like the mean intensity, the value characterizes a model which saturates the measured IGRB auto-correlation.

In Fig. 12 we show the posterior probabilities for in the three energy bands considered in our analysis, and . We find that the typical value is between 20% and 50% in the two lower energy bands (upper panels) whereas for GeV is in the broader range 10% – 80%. These results are robust to the details of the fitting procedure, as demonstrated by the similarity of the various curves. They imply that the extragalactic sources considered in our model (BL Lac, mAGN, SFG, FSRQ and DM) which, as we have seen, provide a good match to the observed cross correlation with LSS tracers, also contribute to a significant fraction of the IGRB, although possibly not to the whole signal. This result is interesting but should also be taken with a grain of salt given the complexity of our cross-correlation model. For example [76], using a different model for the SFG emission and bias, was able to account for a larger fraction of the IGRB, although again not 100%. It should be also noted that the measurement of the IGRB in [3] is affected by systematic errors induced by the imperfect model of the foreground Galactic emission, even if the size of this systematic uncertainty does not seem to be large enough to saturate our models to 100% of the total emission. If indeed it turns out that additional -ray sources are required to explain the total intensity of the IGRB, then the results of our analysis set a rather sever constraint: their correlation with LSS tracers must be weak. This would imply that they should be local, possibly of Galactic origin, like the millisecond pulsar or, perhaps, diffuse inverse Compton photons from cosmic-ray electrons scattering on the optical/infrared Galactic inter-stellar radiation field. Future analyses with newer and additional datasets will help to clarify this interesting issue.

The posterior for is instead consistent with unity, although its probability distribution actually spans several orders of magnitude from to 10, meaning that the measured auto-correlation does not provide a very stringent cross-check.

5. Summary and Conclusions

In this paper we have used the cross-correlations recently measured in [76] between Fermi-LAT diffuse -ray maps and different catalogs of LSS-tracers to investigate the origin of the IGRB and the nature of the various sources that may contribute to it, including DM annihilation or decay. This work extends that of [59] which used only the -ray-2MASS correlation and considered DM as the only source of the extragalactic -ray signal. Our main results are as follows:

• Our theoretical models provide a good fit to the cross-correlation measured in all employed catalogs of extragalactic tracers, namely SDSS-DR6 quasars, 2MASS, NVSS, SDSS-DR8 LRGs and SDSS-DR8 MG. The quality of the fit is quantified by means of a analysis in which we account for covariance among the errors in different angular bins whereas we ignore the covariance among energy bins and among the different catalogs. The first approximation is justified by the photon statistics, which is dominated by low energy event, making each of the energy bins considered in our analysis effectively independent. The second approximation is justified by the spatial distributions of the objects in the different catalogs that, with the partial exception of the NVSS one, do not significantly overlap with each other.

• In our cross-correlation function (CCF) models we consider four different types of astrophysical sources (two flavours of blazars, FSRQs and BL Lacs, SFGs and mAGNs) and, in addition, annihilating/decaying DM. The rationale behind the choice of these astrophysical sources is that previous analyses have shown that they are the main contributors to the IGRB and its angular auto-correlation. These two observational constraints are not considered in our fit. Instead, we use them a posteriori to check the consistency of our best fitting models which are based solely on the measured CCF. We find that models that provide a good match to the cross-correlation fall short of accounting for the mean -ray intensity. The discrepancy is not large, less than a factor of two, especially in the high energy band, and could be accounted for by a combination of model uncertainty and imperfect subtraction of the Galactic foreground. However, it may also indicate that additional types of sources that do not cross-correlate with the LSS, like -ray sources within our Galaxy, are required to account for the whole IGRB intensity.

• Including DM among the possible IGRB sources does not significantly improve the quality of the fit, and does not indicate a preference for a particular DM mass or annihilation cross-section/decay rate. We find that the reason for the low statistical significance on the presence of a DM component does not lie in the fact that the fit rejects this component, while it is rather due to the presence of a model degeneracy with other types of astrophysical sources, mainly mAGNs and SFGs. In other words, a significant DM contribution gives an equally good fit as a case with a negligible DM contribution and a larger mAGNs and SFGs emission. Neglecting the mAGN component in the fit partially breaks this degeneracy and provides a small () preference for DM. The best fit is found for a rather canonical WIMP DM candidate with GeV that annihilates into at a rate which is of the order of the thermal value for the benchmark low DM clustering scenario considered. A candidate with a slight smaller mass of about 30 GeV that annihilates into provides an equally good fit.

• Breaking this degeneracy is the main goal of future cross-correlation analyses similar to the present one. Fortunately, this is a realistic goal. One of the main reason for this degeneracy is the uncertainty on the mAGN and SFG luminosity in -ray which, to date, has been directly measured for a handful of very nearby objects. However, their number is bound to increase thanks to the fact that Fermi-LAT will keep taking data in the next few years. In addition, the quality of the Fermi maps is also expected to increase both in terms of photon statistics, which will allow to better sample the energy behaviour and to improve the sensitivity to characteristic DM spectral features, and angular resolution, which would allow us to push the correlation analysis to smaller angular scales where the 1-halo term dominates.

• We turn the non-detection of DM into limits on the annihilation cross-section/decay rate as function of the DM mass. Our derived constraints are comparable in strenght to most of the current indirect detection method that exploits the -ray sky [59]. These constraints are rather robust to the astrophysical details of the models but, as expected, do depend on the detail of the DM substructure and small-scale clustering. For this reason and with the aim of bracketing current theoretical uncertainties, in addition to the low scenario which represents the current, somewhat conservative, benchmark substructure model, we have explored two additional, rather extreme cases: the ns case in which we completely ignore substructures and that provides extremely conservative constraints of the DM properties, and the high scenario in which substructures are more numerous and have an higher density concentration. In the most conservative ns scenario our method excludes, at a credible level larger than 95%, that DM particles with masses smaller than 10 GeV annihilating entirely into could have a thermal cross section. In the optimistic scenario, the same statement applies to particles lighter than GeV. The bounds are a factor of stronger than the most conservative case considered in [59] in which only DM is used in order to saturate the 2MASS cross-correlation. Constraints on DM decay time for DM decaying into are s, roughly independently from the DM mass.

All in all we are confident that the results obtained in our analysis, which are already quite remarkable considering that this is the first time that a genuine cross-correlation signal is detected in the Fermi-LAT -ray maps, will soon improve significantly. In this respect this work also represents a proof of concept that illustrates the potential of the cross-correlation analysis. We base our optimism on the fact that the new Pass8 data, with improved effective area and angular resolution, will soon be released by the Fermi-LAT Collaboration and that additional catalogs of objects al relatively low redshifts with wide, almost all-sky, angular coverage and well determined redshift distribution are already available [18] and some new ones are being compiled [17].

Acknowledgements

This work is supported by the research grant Theoretical Astroparticle Physics number 2012CPPYP7 under the program PRIN 2012 funded by the Ministero dell’Istruzione, Università e della Ricerca (MIUR), by the research grants TAsP (Theoretical Astroparticle Physics) and Fermi funded by the Istituto Nazionale di Fisica Nucleare (INFN), and by the Strategic Research Grant: Origin and Detection of Galactic and Extragalactic Cosmic Rays funded by Torino University and Compagnia di San Paolo. MV and EB are supported by PRIN MIUR and IS PD51 INDARK grants. MV is also supported by ERC-StG cosmoIGM, PRIN INAF JX is supported by the National Youth Thousand Talents Program, the National Science Foundation of China under Grant No. 11422323, and the Strategic Priority Research Program, The Emergence of Cosmological Structures of the Chinese Academy of Sciences, Grant No. XDB09000000.

Appendix A Window Functions

In this Appendix we discuss the modeling of the window functions adopted for the calculation of the cross-correlation angular power spectrum of Eq. (1), which are in turn the ingredient for the determination of the cross correlation function defined in Eq. (4).

a.1. Dark matter

a.1.1 Annihilating dark matter

DM annihilations in haloes and in their substructures produce -ray photons. This emission traces the DM density squared : therefore the density field responsible for the correlation signal is . The window function reads:

 Wδ2(χ)=(ΩDMρc)24π⟨σav⟩2mDM2[1+z(χ)]3Δ2(χ)∫Eγ>EmindEγdNadEγ[Eγ(χ)]e−τ[χ,Eγ(χ)], (A1)

where is the cosmological abundance of DM, is the critical density of the Universe, is the mass of the DM particle, and denotes the velocity-averaged annihilation rate, assumed here to be the same in all haloes. indicates the number of photons produced per annihilation event, and sets the -ray energy spectrum. We will consider annihilation into quarks as representative of a typical soft annihilation spectrum (with -rays mostly arising from production and decay of neutral pions), and into leptons as representative of a hard-spectrum channel (where -rays mostly arising from final state radiation), with and final states as intermediate possibilities. is the energy threshold of the Fermi-LAT maps considered in the analysis, namely: GeV. The factor accounts for absorption due to the extra-galactic background light, and we model the optical depth as in [36].

A crucial quantity in Eq. (A1) is the so-called clumping factor :

 (A2)

The clumping factor involves the integral of the halo number density above the so-called minimal halo mass , multiplied by the total number of annihilations produced in the generic haloes of mass at redshift with density profile and with subhalos providing a “boost” to the emission given by . We assume a reference value of for , which corresponds to a typical free-streaming mass in the WIMP DM scenario. We adopt the halo mass function from [65] and we assume that the halos are characterized by the so-called Navarro-Frenk-White (NFW) universal density profile [53]. The profile is completely determined by the total mass of the halo and by its size. We express the latter in terms of the concentration parameter , taken from [58] (see also [64] for an analytic fit of of [58]).

Concerning the boost provided by subhaloes hosted in the main haloes, we consider three scenarios (high, low, and ns) as extreme cases bracketing the effect. In Eq. (A2) this is indeed the most uncertain quantity [64, 37, 33, 54]. In the low scenario, the function is computed following [64] (see, in particular, their Eq. (2) assuming a subhalo mass function ). The high scenario stems instead from the found in [37], and assuming no redshift dependence. In the ns case, we simply set : this can be considered as the most conservative approach.

The blue curves in Fig. 13a show a quantity related to the window function of Eq. (A1), defined as the kernel of the -ray emission entering in the computation of the angular power spectrum discussed in Section 2. Specifically, we define the kernels as:

 Ka=√cz/H(z)Wa(z)/χ(z), (A3)

such that:

 C(γigj)ℓ=∫dln(z)Kγi(z)Kgj(z)Pγigj(k,z). (A4)

In Fig. 13a we choose a reference particle-physics model with GeV, and annihilation channel. The three clustering scenarios (ns, low and high for dotted, solid and dashed curves, respectively) share approximately the same redshift dependence, but they correspond to different sizes of the clumping factor and consequently of the intensity of the DM-induced -ray flux. Note that a comparison with previous works in the literature can be non-trivial, as different groups employ different prescriptions for the ingredients of the DM clustering and, in particular, for the boost factor. Fig. 13 can be useful also as a normalization test when confronting the results presented in the rest of the paper with other works.

a.1.2 Decaying dark matter

If instead of being stable, the DM particles decay, while having a negligible self-annihilation rate, the produced -rays traces the DM density linearly, i.e., , The window function in this case reads:

 Wδ(χ)=ΩDMρc4πΓdmDM∫Eγ>EmindEγdNddEγ[Eγ(χ)]e−τ[χ,Eγ(χ)], (A5)

where is the decay rate. The photon yield, , is assumed to be the same as for annihilating DM, but with the energy of the process given by instead of . In other words, with the kinematic end-point being at . The kernel in the case of decaying DM is shown as a cyan curve in Fig. 13a. In the plot we report reference particle-physics model with GeV, and decays into quarks. Note that for decaying DM, the window function does not depend on the details of the DM clustering. We notice also that DM kernels peak at low redshifts, both for annihilating and decaying DM, and have a relative fast decrease with distance.

a.2. Astrophysical sources

For astrophysical sources, we adopt as the characterizing parameter the source -ray luminosity in the energy interval (0.1 – 100) GeV. For a power-law energy spectrum with spectral index , the window function takes the form:

 WSi(χ)=(αi−2)⟨fSi(χ)⟩4πE20[1+z(χ)]2∫Eγ>EmindEγ(EγE0)−αie−τ[χ,Eγ(χ)], (A6)

where MeV is just the normalization energy, and stands for each of the -rays sources adopted in our analysis: BL Lac, FSRQ, mAGN and SFG. The mean luminosity produced by an unresolved class of objects located at a distance from us is denoted by and is given by:

 (A7)

where is the -ray luminosity function for the source class . The upper bound, , is the luminosity above which an object can be resolved, given the detector sensitivity for which we assume the value above 1 GeV [55, 1]. The precise value depends slightly on and on the catalogue of resolved point sources, although varying within these different values has only a weak impact of the window function. Conversely, the minimum luminosity depends on the properties of the source class under investigation. The four populations of astrophysical -ray emitters (i.e., BL Lac, FSRQ, mAGNs and SFGs) are discussed in the following. For each of them we describe the choice of and of the -ray luminosity function.

a.2.1 Blazars

We consider BL Lacertae (BL Lacs) and flat-spectrum radio quasars (FSRQ) separately. The -ray luminosity function of BL Lacs and FSRQ is taken from [8] and [9], respectively, where it is derived from a parametric fit of the redshift and luminosity distributions of resolved blazars in the Fermi-LAT catalogue. The lower limit of the integral in Eq. (A7) is set to (BL Lac) and (FSRQ). For the energy spectrum, we consider a simple power-law with a spectral index taken from the average spectral index in [8, 9], namely, we assume and .

The kernels of unresolved blazars are shown by the solid red (BL Lac) and magenta (FSRQ) lines in Fig. 13a. Note that they strongly decrease at low since Fermi-LAT has already detected a large number of the closest (brightest) emitters of these classes.

a.2.2 Misaligned AGNs

In the case of mAGN, we follow [27], which studied the correlation between the -ray luminosity and the core radio luminosity at 5 GHz, and derived the GLF from the radio luminosity function. We consider their best-fit vs. relation and assume an average spectral index of 2.37. The solid green line in Fig. 13a indicate the contribution of unresolved mAGNs.

a.2.3 Star-forming galaxies

As done in [4], we assume that the -ray and infrared (IR) luminosities are correlated in the case of SFG. We adopt the best-fit vs. relation from [4] while for the IR luminosity function we adopt the one from [39] , (adding up spiral, starburst, and SF-AGN populations of their Table 8), as considered in [72]. The spectral index is taken to be for all the 3 components although starbursts galaxies would require in principle a somewhat harder spectrum. Nonetheless, this component is subdominant in the total SFG contribution except for high energies and at high redshift (i.e., in the ranges which are less relevant for the analyses in our work). The above choice has thus no practical effects on our results. The kernel associated to unresolved SFGs is the solid orange line in Fig. 13a. All the different single peaked sub-populations provide sizable contributions and this gives raise to different peaks.

The -ray emission produced by the four extragalactic astrophysical populations described above accounts for approximately the whole IGRB and autocorrelation angular power spectrum (see Fig. 1). As described in the main text we however introduced a normalizing constant for each population to be determined by the fit. Apart from the extragalactic DM-induced emission described in Secs A.1.1 and A.1.2, there may be a contribution associated with annihilations/decays in the DM halo of the Milky Way. This is not included since it does not correlate with the LSS tracers.

a.3. Galaxy catalogues

For galaxies, we take the redshift distributions reported in [76]. The associated kernels are shown in Fig. 13b. The 2MASS kernel peaks at low-redshift, and a comparison with the -rays kernels shown in the left panel of the same Fig. 13 indicates that the 2MASS catalogue is the most suitable for investigating a DM signal in the cross-correlation analysis, followed by the SDSS Main Galaxy Sample catalogue.

Appendix B Halo occupation distribution of galaxies

In this work, we compute the angular cross-correlation between the unresolved -ray sky and the number of galaxies in specific catalogues. In order to estimate the latter from a theoretical point of view (and since we adopt the halo model description for the structure clustering), we need to describe how galaxies populate halos. Namely, we need to model how many galaxies of a certain catalogue are present in a halo of mass and how they are spatially distributed. To this aim, we employ the halo occupation distribution (HOD) formalism.

We follow the approach described in [78] (for review on HOD, see also [16, 24]), where the HOD is parameterized by distinguishing the contributions of central and satellite galaxies, (since different formation histories typically imply different properties for galaxies residing at the centers of halos with respect to satellite galaxies). These can be modeled with the following functional forms:

 ⟨Ncen(M)⟩ = 12[1+erf(logM−logMthσlogM)] (B1) ⟨Nsat(M)⟩ = (MM1)αexp(−McutM) (B2)

With this formalism, we need five parameters for each galaxy population: denotes the approximate halo mass required to populate the halo with the considered type of galaxies, with the transition from 0 to 1 central galaxy modeled by means of Eq. (B1), and set by the width . The satellite occupation is described by a power law (with index and normalization set by the mass ), with an exponential cutoff at low masses. The value of the five HOD parameters for each of the considered galaxy population is discussed in the following. For some catalogues, we will also consider similar but slightly different functional forms.

We selected those galaxy samples with available HOD which more closely resemble the catalogues considered in the cross-correlation analysis of this work. We caution, however, that since the matching of the two samples is not perfect some differences in the associated HODs might be expected. Nevertheless, this should not affect our results in a dramatic way.

Eqs. (B1) and (B2) provide the number of galaxies in a halo of mass . Concerning the spatial distribution, we treat central and satellite galaxies separately. The former is taken as a point-source located at the center of the halo (the point-source approximation is expected to break down only for ). Satellite galaxies are instead described in an effective way with a spatial distribution following the host-halo profile. In other words, we express the density field of galaxies with:

 gg(x−x′|M)=⟨Ncen(M)⟩δ3(x−x′)+⟨Nsat(M)⟩ρh(x−x′|M)/M. (B3)

Note that:

 ∫d3xgg(x)=⟨Ncen(M)⟩+⟨Nsat(M)⟩=⟨N(M)⟩. (B4)

b.0.1 2mass Hod

A determination of the HOD for the 2MASS galaxies is not present in the literature (to the knowledge of the authors). In [77], a sample of about 200,000 SDSS galaxies mostly residing in the redshift range and with r-band magnitude was analyzed. Such ranges of redshift and magnitude are analogous to the ones of the adopted 2MASS catalogue [43], with the cross-identification of the latter with SDSS found to be successful for about 90% of the sources [52]. We can thus exploit the HOD results of [77]. They considered a step function ( for and for ) instead of Eq. (B1), and set in Eq. (B2). The analysis was performed by splitting the sample in luminosity bins, but for our purposes we can consider the averaged best-fit parameters weighted over the number of galaxies in each bin. We found , , and .

b.0.2 Nvss Hod

The NVSS sub-sample considered in this work ([19]) contains sources brighter than 10 mJy at 1.4 GHz. The vast majority of them is associated to bright AGNs. To model the AGN HOD, we follow [23]. For bright objects, they found , , , , and .

b.0.3 SDSS-DR8 Main Galaxy Sample HOD

In [63], the clustering of more than three million photometrically selected SDSS galaxies was analyzed. In particular, the sample was defined requiring de-reddened r-band magnitudes and absolute magnitudes , in the redshift range and masking objects with Galactic extinction . This galaxy sample is very similar to the one adopted in this work for the cross-correlation analysis ([6]), except for a more limited redshift (and magnitude) range. Since the peak of the redshift distribution is at , such difference is not expected to play a major role. After averaging the best-fit values of HOD parameters over the different redshift bins considered in [63] (weighted for the number of galaxies in each bin), we obtain , , , and . The functional form of the satellite HOD considered in [63] is: