Searching for gamma-ray emission from galaxy clusters at low redshift

TL;DR

This paper detects a significant cross-correlation between unresolved gamma-ray emission and galaxy clusters at low redshift, suggesting potential astrophysical sources like AGNs or intracluster medium emission.

Contribution

It introduces a robust method for estimating the cross-correlation covariance matrix using mock realizations, applicable to future cross-correlation analyses.

Findings

Detected a 3.5 sigma cross-correlation signal.

Identified potential gamma-ray sources within galaxy clusters.

Developed generalized mock data techniques for covariance estimation.

Abstract

We report the identification of a positive cross-correlation signal between the unresolved $\gamma$-ray emission, measured by the \emph{Fermi} Large Area Telescope, and four different galaxy cluster catalogues. The selected catalogues peak at low-redshift and span different frequency bands, including infrared, optical and X-rays. The signal-to-noise ratio of the detected cross-correlation amounts to 3.5 in the most significant case. We investigate and comment about its possible origin, in terms of compact $\gamma$-ray emission from AGNs inside clusters or diffuse emission from the intracluster medium. The analysis has been performed by introducing an accurate estimation of the cross-correlation power-spectrum covariance matrix, built with mock realisations of the gamma and galaxy cluster maps. Different methods to produce the mock realizations starting from the data maps have been investigated and compared, identifying suitable techniques which can be generalized to other cross-correlation studies.

Tables

Table 1. Table 1: Energy bins in GeV used in our analysis. E min subscript 𝐸 min E_{\rm min} and E max subscript 𝐸 max E_{\rm max} denote the lower and upper bound of the bins, while ℓ min subscript ℓ min \ell_{\mathrm{min}} and ℓ max subscript ℓ max \ell_{\mathrm{max}} show the interval in multipole ℓ ℓ \ell over which the fit of the angular power spectrum is performed: the lower bound is chosen in order to exclude possible galactic-foreground residual contamination, the upper limit is driven by the Fermi -LAT PSF, whose 68% containment angle θ cont subscript 𝜃 cont \theta_{\rm cont} is reported.

Bin	$E_{\min }$ [GeV]	$E_{\max }$ [GeV]	${\ell }_{\min }$	${\ell }_{\max }$	${\theta }_{\mathrm{cont}}$(deg)
1	0.631	1.202	40	251	0.50
2	1.202	2.290	40	316	0.58
3	2.290	4.786	40	501	0.36
4	4.786	9.120	40	794	0.22
5	9.120	17.38	40	1000	0.15
6	17.38	36.31	40	1000	0.12
7	36.31	69.18	40	1000	0.11
8	69.18	131.8	40	1000	0.10
9	131.8	1000	40	1000	0.10

Table 2. Table 2: Cluster catalogues used in our analysis. For each catalogue, we show the total number of clusters after the selection discussed in the text, the signal-to-noise ratio for cross-correlation with the Fermi -LAT maps and the difference in the best-fit χ 2 superscript 𝜒 2 \chi^{2} ’s of a flat angular power spectrum in multipole and the physical AGN model.

Catalogue	Type	Reference	Total	SNR	$\Delta {\chi }_{\mathrm{FLAT}-\mathrm{AGN}}^2$
WHY18	infrared	Wen et al. (2018)	8999	2.1	0.6
SDSSDR9	optical	Banerjee et al. (2018)	2582	2.3	0.2
MCXCsub	X-ray	Reiss & Keshet (2018)	109	3.5	4.5
HIFLUGCS	X-ray	Reiprich & Böhringer (2002a)	105	3.2	0.6

Equations

$$F_{\Delta E}^{\gamma }\exp \left(-\frac{R^2}{2{\theta }_{\Delta E}^2}\right)>\frac{F_{\Delta E,\mathrm{faintest}}^{\gamma }}{5},$$

$$C_{\ell ,\mathrm{FLAT}}^{c_j{\gamma }_i}=\left(\frac{1}{N_{c_j}}\sum _{k=0}^{N_{c_j}}{F}_{\Delta E_i,k}^{\gamma }\right)W_{\ell }^{\Delta E_i},$$

$$C_{\ell ,\mathrm{FLAT}}^{c_j{\gamma }_i}=C_P^j{\bar{E}}_i^{-{\alpha }_{1,j}}\Delta E_i{W}_{\ell }^{\Delta E_i},$$

$$C_{\ell }^{c_j{\gamma }_i}=\int \frac{d\chi }{{\chi }^2}{W}_{c_j}(\chi )W_{{\gamma }_i}(\chi )P_{c_j{\gamma }_i}(k=\ell /\chi ,\chi ),$$

$$W_{c_j}(z)=\frac{4\pi \chi (z{)}^2}{N_{c_j}}\int dM\frac{d^2{n}_{c_j}}{dMdV},$$

$$N_{c_j}=\int dMdV\frac{d^2{n}_{c_j}}{dMdV}.$$

$$W_{{\gamma }_i}(E,z)=\chi (z{)}^2{\int }_{{\mathcal{L}}_{\min }}^{{\mathcal{L}}_{\max }(F_{\mathrm{sens}},z)}d\mathcal{L}{\Phi }_{\mathrm{S}}(\mathcal{L},z,E)\frac{d{N}_{\mathrm{S}}}{dE}\left(\mathcal{L},z\right)\times e^{-\tau \left[E(1+z),z\right]},$$

$$P_{c_j,{\gamma }_j}^{1h}(k,z)={\int }_{{\mathcal{L}}_{\min ,i}(z)}^{{\mathcal{L}}_{\max ,i}(z)}d\mathcal{L}{\varphi }_i(\mathcal{L},z)\times \frac{\mathcal{L}}{⟨f_{{\gamma }_i}⟩}\frac{⟨N_{c_j}(M(\mathcal{L}))⟩}{{\bar{n}}_{c_j}}$$

$$P^{2h}_{c_{j},\gamma_{i}}(k,z)=\biggl{[}\int^{\mathcal{L}_{\mathrm{max},i}(z)}_{\mathcal{L}_{\mathrm{min},i}(z)}\ d\mathcal{L}\ \Phi_{i}(\mathcal{L},z)b_{h}(\mathcal{M(L)})\frac{\mathcal{L}}{\langle f_{\gamma_{i}}\rangle}\biggr{]}\times\biggl{[}\int_{M_{\mathrm{min}}}^{M_{\mathrm{max}}}\ dM\ \frac{dn}{dM}b_{h}(M)\frac{\langle N_{c_{j}}\rangle}{\bar{n}_{c_{j}}}\biggr{]}\ P^{\mathrm{lin}}(k)\,,$$

$$C_{\ell,\mathrm{AGN}}^{c_{j}\gamma_{i}}\ =\ C_{\ell,{\mathrm{FLAT}}}^{c_{j}\gamma_{i}}(C_{p}^{j},\alpha_{1,j})+\ A_{j}\ \Bigl{(}\frac{\bar{E}_{i}}{\bar{E}_{0}}\Bigr{)}^{-\alpha_{2,j}}\ \Delta E_{i}\ C_{\ell}^{\Delta E_{0}}\ W_{\ell}^{\Delta E_{i}}\,;$$

$$\mathrm{SNR}=\frac{C_{\ell ,\mathrm{DATA}}{\Gamma }_{\ell {\ell }^{\prime }}^{-1}{C}_{\ell ,\mathrm{MODEL}}}{\sqrt{C_{\ell ,\mathrm{MODEL}}{\Gamma }_{\ell {\ell }^{\prime }}^{-1}{C}_{\ell ,\mathrm{MODEL}}}},$$

$$\chi^{2}=\sum_{i=1}^{N_{\rm bin}}\sum_{\ell=\ell_{\rm min}}^{\ell_{\rm max}}\sum_{\ell^{\prime}=\ell_{\rm min}}^{\ell_{\rm max}}\ \Bigl{[}(C_{\ell,\mathrm{DATA}})^{i}-(C_{\ell,\mathrm{MODEL}})^{i}\Bigr{]}\ (\textbf{$\Gamma$}_{\ell\ell^{\prime}}^{-1})^{i}\ \Bigl{[}(C_{\ell^{\prime},\mathrm{DATA}})^{i}-(C_{\ell^{\prime},\mathrm{MODEL}})^{i}\Bigr{]}\,,$$

$$-2\ln \mathcal{L}={\chi }^2.$$

$$\mathrm{AIC}=2k-2\ln \mathcal{L},$$

$${\Gamma }_{\ell {\ell }^{\prime }}^{c{\gamma }_i}\equiv \text{cov}[C_{\ell }^{c{\gamma }_i},C_{{\ell }^{\prime }}^{c{\gamma }_i}]=\sum _{k=1}^{N_{\gamma }}\sum _{j=1}^{N_c}(C_{\ell }^{k,c{\gamma }_i}-{\overline{C_{\ell }}}^{c{\gamma }_i})(C_{{\ell }^{\prime }}^{j,c{\gamma }_i}-{\overline{C_{{\ell }^{\prime }}}}^{c{\gamma }_i}),$$

$${\overline{C_{\ell }}}^{c{\gamma }_i}=\frac{1}{N}\sum _{a=1}^{\mathrm{N}}{C}_{\ell }^{a,c{\gamma }_i}$$

$$\Gamma_{\ell\ell^{\prime}}^{c\gamma_{i}}=\frac{1}{2}\Bigl{(}\Gamma_{\ell\ell^{\prime}}^{\hat{c}\gamma_{i}}+\Gamma_{\ell\ell^{\prime}}^{c\hat{\gamma_{i}}}\Bigr{)}\,.$$

$${\Gamma }_{\ell {\ell }^{\prime }}^{c{\gamma }_i}=\frac{C_{\ell }^{cc}{C}_{\ell }^{{\gamma }_i{\gamma }_i}+(C_{\ell }^{c{\gamma }_i}{)}^2}{(2\ell +1)\Delta \ell f_{\mathrm{sky},\mathrm{i}}}{\delta }_{l{l}^{\prime }}^{\mathrm{K}},$$

$$r=\frac{{\text{Γ}}_{\ell {\ell }^{\prime }}^{c{\gamma }_i}}{\sqrt{{\text{Γ}}_{\ell \ell }^{c{\gamma }_i}{\text{Γ}}_{{\ell }^{\prime }{\ell }^{\prime }}^{c{\gamma }_i}}}.$$

$${\Gamma }_{\ell \ell }^{c{\gamma }_i{\gamma }_j}\equiv \mathrm{cov}[C_{\ell }^{c{\gamma }_i},C_{\ell }^{c{\gamma }_j}]=\frac{C_{\ell }^{cc}{C}_{\ell }^{{\gamma }_j{\gamma }_i}+C_{\ell }^{c{\gamma }_i}{C}_{\ell }^{c{\gamma }_j}}{(2\ell +1)f_{\mathrm{sky}}},$$

$$r_{\mathrm{E}}=\frac{{\Gamma }_{\ell \ell }^{c{\gamma }_i{\gamma }_j}}{\sqrt{{\Gamma }_{\ell \ell }^{c{\gamma }_i{\gamma }_i}{\Gamma }_{\ell \ell }^{c{\gamma }_j{\gamma }_j}}}$$

$$P_{E_i}=\frac{1}{\Delta \ell }\sum _{{\ell }_{\min }}^{{\ell }_{\max }}{C}_{\ell }^{c{\gamma }_i},$$

$$\sigma^{2}_{P_{E_{i}}}=\frac{1}{\Delta\ell}\Bigl{(}\frac{1}{\Delta\ell}\sum_{\ell_{\rm min}}^{\ell_{\rm max}}\sigma^{2}_{C_{\ell}^{c\gamma_{i}}}\Bigr{)}\,,$$

$$\widehat{{\Gamma }_{\mathrm{B}}}=\frac{1}{N_r-1}\sum _{k=1}^{N_k}(x_i^k-{\bar{x}}_i)(x_j^k-{\bar{x}}_j)$$

$${\bar{x}}_i=\frac{1}{N_r}\sum _{k=1}^{N_r}{x}_i^k.$$

$$\widehat{{\Gamma }_{\mathrm{J}}}=\frac{N_{\mathrm{r}}-1}{N_{\mathrm{r}}}\sum _{k=1}^{N_{\mathrm{r}}}(x_i^k-{\bar{x}}_i)(x_j^k-{\bar{x}}_j)$$

$$f(\theta ,\varphi )=\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }{a}_{\ell m}{Y}_{\ell m}(\theta ,\varphi ).$$

$$C_{\ell }=\frac{1}{2\ell +1}\sum _{m=-1}^{\ell }|a_{\ell m}{|}^2.$$

$$a_{\ell m}⟶a_{\ell m}{e}^{i{\phi }_{\ell m}}{\phi }_{\ell m}\in \mathbb{R};$$

$$\text{cov}[\text{X}]=\frac{1}{N_{\mathrm{r}}}\sum _{k=1}^{N_{\mathrm{r}}}(x_i^k-{\bar{x}}_i)(x_j^k-{\bar{x}}_j)$$

Full text

1 Introduction

Galaxy clusters are one of the most important tracers of the Large Scale Structure (LSS) of the Universe, being the largest virialized objects formed by the gravitational instability. Because of their large dimension, mass and formation history, they represent a unique cosmological probe. They are also fundamental from an astrophysical point of view: they host galaxies, but also ionized hot gas thermalized via collisionless virial shocks, dark matter (DM) and relativistic cosmic rays (CRs) accelerated by the shocks present at the edge of the clusters. If we focus on $\gamma$-ray emission, the CRs can produce photons via inverse Compton, non-thermal bremsstrahlung and decay of $\pi^{0}$. DM particles can also induce $\gamma$-rays through the same mechanisms arising from the products of DM annihilation or decay.

Within the hierarchical structure formation process, clusters form in the node of the cosmic web, where also different populations of astrophysical objects, as the Active Galactic Nuclei (AGN), are present. These astrophysical sources, that can be found within the clusters themselves, contribute to the total $\gamma$-ray flux that we observe. On top of this emission, clusters could host a spatially extended $\gamma$-ray contribution. The detection of a signal of $\gamma$-rays from CRs in clusters could help in understanding the origin of the radio halos (e.g. van Weeren et al. (2019) for a review). The observation of $\gamma$-rays originated from annihilation or decay of DM particles would confirm their existence, whilst a non-detection can help to reduce the size of the bucket of DM particle models.

The study of the $\gamma$-rays from galaxy clusters requires a double effort from an observational point of view: we need accurate galaxy cluster catalogues (in terms of mass, dimension and position), in order to be able to distinguish between the DM halo, the Intra-Cluster Medium (ICM) and the sources (such as AGN) inside the cluster, and precise measurements of the $\gamma$-ray photon flux. From the cluster catalogues side there are several surveys in the literature, obtained with different telescopes and at different wavelengths, that can be used for this purpose, while from the $\gamma$-ray side the most detailed sky recognition comes from the Fermi Large Area Telescope (LAT) Ackermann et al. (2010). Most recent analyses of the Fermi-LAT  data looking for a $\gamma$-ray signal from galaxy clusters using different techniques include (Branchini et al., 2017; Reiss et al., 2018; Brunetti et al., 2017; Lisanti et al., 2018; Hashimoto et al., 2019).

In this work we focus on the information we can derive from the joint analysis of the two observables, the $\gamma$-ray photon flux and the galaxy clusters distribution; in other words we proceed with a cross-correlation analysis. A similar approach was undertaken by Branchini et al. (2017) and Hashimoto et al. (2019), but using different cluster samples. Several analyses in the literature have been studying the cross-correlation of $\gamma$-ray with other LSS tracers (Ando et al., 2014; Shirasaki et al., 2014; Fornengo et al., 2015; Regis et al., 2015; Cuoco et al., 2015; Shirasaki et al., 2015; Ando & Ishiwata, 2016; Shirasaki et al., 2016; Feng et al., 2017).

Here we study the cross-correlation angular power spectrum (APS) between four selected galaxy cluster catalogues obtained in different bands (optical, infrared and X-ray) and $\gamma$-rays from the Fermi-LAT in different energy bins. We focus on low redshift catalogues and high-mass clusters, since our main goal is to disentangle a possible (and long-sought after) extended $\gamma$-ray emission from clusters. Such a signal would be originated from ICM or DM, since AGNs and galaxies have much more compact emissions. Improving from the previous works listed above, we introduce an accurate estimation of the power spectrum covariance matrix. This is built with mock realisations of the gamma and galaxy cluster maps and allows a precise statistical evaluation of the significance of the measured APS.

2 Data

The cross-correlation analyses we have carried out in this paper are based on: (i) The full-sky $\gamma$-ray intensity emission measured by the Fermi-LAT , for which we consider 9-years of data in the energy range between 630 MeV and 1 TeV; (ii) A series of low-redshift galaxy cluster catalogues built in different electromagnetic bands.

2.1 Fermi-LAT $\gamma$-rays maps

The Fermi-LAT pair-conversion telescope achieved remarkable results for $\gamma$-ray astronomy during its 10 years of operations. Its large energy coverage (20 MeV - 1 TeV) and the capability of rejecting the contamination from charged cosmic rays make the instrument particularly suitable to investigate the nature of the unresolved extra-galactic $\gamma$-ray background (UGRB). The angular resolution of the instrument is energy dependent and reaches $\sim$ 0.1∘ above 10 GeV. The photon and exposure maps adopted in the present analysis are produced with the Fermi Tools111We used the LAT Science Tools version https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/. The current version of the Fermi Tools subdivide the photon events into quartiles of angular resolution, from PSF0 to PSF3, corresponding to a progress from the worst to the best Point Spread Function (PSF). We find a trade-off between photon statistics and angular resolution, by selecting the best quartile PSF3 for energies below 1.2 GeV (where we have the highest photon counts, but worst PSF) and PSF1+2+3 for higher energies. In order to have the lowest contamination from cosmic-rays, we selected the Pass8 ULTRACLEANVETO event class222See http://www.slac.stanford.edu/exp/glast/groups/canda/lat_ Performance.htm for further details on photon event classes., that is recommended for diffuse emission analysis.

In this work we used 108 months of data (from mission week 9 to week 476). The analyses are performed on photon intensity maps, obtained by dividing the count maps by the exposure maps and the pixel area $\Omega_{\rm pix}=4\pi/\rm N_{\rm pix}$. We adopt a HEALPix (Gorski et al., 2005) pixelation format with resolution parameter $\rm N_{\rm side}=$ 1,024, which corresponds to a total number of pixel $\rm N_{\rm pix}$ = 12,582,912 and a mean spacing of $\sim 0.06\degree$, similar to the best angular resolution of the $\gamma$-ray data. We derived the intensity maps in 100 energy bins, evenly spaced in logarithmic scale between 100 MeV and 1 TeV. These micro-bins were subsequently re-binned into 9 larger energy bins, from 631 MeV to 1 TeV, which are listed in table 1. The data selection and the pre-processing steps follow the same procedure described in Ammazzalorso et al. (2018).

Since there is no clear sign of resolved extended emission from clusters, we adopt a cross-correlation technique which focusses on the UGRB component. To this aim, we need to mask resolved point sources. Moreover, we mask the Galactic emission, which acts as a foreground and, while not correlating with the extra-galactic cluster distribution, nevertheless contributes a sizeable source of noise to the error budget. We therefore build a set of masks for the $\gamma$-ray maps by adopting the following criteria:

• •

We mask resolved sources by adopting to the 4FGL catalogue (The Fermi-LAT collaboration, 2019a) that contains 5523 sources. Above 10 GeV, we include also additional sources present in the 3FHL catalogue (Ajello et al., 2017), that is specific to high energies. Each source is masked taking into account both the source brightness and the PSF resolution in the specific bin, as done in Ammazzalorso et al. (2018).

• •

The galactic disk emission is masked by means of a latitude cut that excludes the portion of the sky with $|b|<30^{\circ}$.

For point sources, the masking radius around each catalogues sources is defined as:

[TABLE]

where $F_{\Delta E}^{\gamma}$ is integral flux of the source in the specific energy bin $\Delta E$ under consideration, $F_{\Delta E,\rm faintest}^{\gamma}$ is the flux of the faintest source in the same energy bin, and $\theta_{\Delta E}$ is the 68% containment angle in that energy bin, as provided by the Fermi-LAT PSF. The resulting energy-dependent masks aim at properly covering resolved sources and avoiding artefacts due to source flux leakage outside the mask, but at the same time maintaining a good sky coverage. For further details and impact of the mask, see also Ackermann et al. (2018).

Even though we adopt a mask to cover the brightest part of the galactic plane emission, we nevertheless additionally adopt a procedure of foreground removal at higher latitudes, in order to reduce the contribution of this component to the error budget. This foreground removal is obtained by subtracting a model of the galactic foreground contribution, for which we use the template maps provided by the Fermi-LAT  Collaboration with the Galactic emission model gll_iem_v06.fits333https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html. The foreground template is projected in HEALPix maps keeping the same $N_{\rm side}$ of the intensity maps and the same micro binning. A free normalisation is assigned to the template map and a free constant is added, the latter representing the average UGRB emission and a possible cosmic-ray contamination of the $\gamma$-ray maps. The resulting foreground model is fitted with a Poissonian likelihood to the masked intensity maps. The best-fit normalisations obtained with this procedure are all well compatible with unity, showing a successful description of the foreground emission. We normalise the foreground templates with the obtained parameters and then subtract the foreground maps from the corresponding intensity maps: this procedure is performed after having re-binned them in the macro energy-bins then used in the cross-correlation analysis. For additional information about the foreground removal, see Ackermann et al. (2018) and Ammazzalorso et al. (2018), where the same procedure is adopted. In figure 1, we show as an example the $\gamma$-ray masked map in the ($2.3-4.8$) GeV energy bin before (left panel) and after (right panel) the galactic foreground removal.

2.2 Galaxy cluster catalogues

We base our analysis on four galaxy cluster catalogues obtained from observations in different frequency bands: one in the infrared (WHY18), one in the optical (SDSSDR9) and two in the X-ray band (MCXCsub and HIFLUGCS). We use these catalogues because they contain clusters which are massive and located at relatively low-redshift, thus offering an expected enhanced cross-correlation signal. Being massive and relatively close, some of them might also extend beyond the angular resolution of the Fermi-LAT detector, allowing us to investigate if a cross-correlation signal due to diffuse emission in clusters is present on large scales.

WHY18 is obtained by combining photometric galaxies from 2MASS, the Wide-field Infrared Survey Explorer (WISE, Wright et al., 2010) and the SuperCOSMOS Sky Survey (Hambly et al., 2001). The main selection applied to the sources is such to include clusters with ${\rm M_{500}}>3\times 10^{14}\rm M_{\odot}$ 444In our analysis we consider ${\rm M_{500}}$ corresponding to an overdensity $\Delta_{c}=500$. The mass in this case is defined as ${\rm M_{500}}=4\pi\,r_{500}^{3}\,500\,\rho_{c}/3=250\,r_{500}^{3}\,H(z)^{2}/G$, where $\rho_{c}$ is the critical density, $r_{500}$ is the virial radius, $H(z)$ the Hubble parameter and $G$ the gravitational constant. and redshift between 0.025 and 0.3. The final number of clusters in this catalogue is 47,500 (Wen et al., 2018).

SDSSDR9 is part of the full Sloan Digital Sky Survey (SDSS, York et al., 2000) catalogue, obtained using an Adaptive Matched Filter (AMF, Kepner et al., 1999) technique, based on a model of galaxy distribution. This filter is built using the cluster density radial profile, the galaxy luminosity function and the redshift. After applying the filter, the catalogues contains 49,479 clusters with redshift between 0.045 and 0.691 (Banerjee et al., 2018);

MCXC is an X-ray catalogue (Piffaretti et al., 2011) obtained by collecting 1,743 clusters from two main types of X-ray observations: (i) contiguous area survey ROSAT (Voges et al., 1999); (ii) deeper pointer X-ray observations. MCXCsub is built from the full MCXC catalogue by selecting a sub-set of 112 clusters with ${\rm M_{500}}>10^{13}\rm M_{\odot}$, angular diameter larger than 0.2°, latitude larger than 20°and in positions of the sky to avoid contamination from bright $\gamma$-ray point-sources (Reiss et al., 2018).

HIFLUGCS contains 63 clusters selected from ROSAT with latitude larger than 20°and flux between 0.1 and 2.4 keV larger than $2\times 10^{-11}\rm ergs\ s^{-1}\ cm^{-2}$ (Reiprich & Böhringer, 2002b).

Starting from the catalogues described above, we perform a selection by considering only clusters at redshift smaller than 0.2 and with mass ${\rm M_{500}}>10^{13}\rm M_{\odot}$, since we aim at investigating clusters having an angular size in the sky at least comparable with the Fermi-LAT  PSF. Moreover, in order to consider only clusters with robust identification, we applied a further cut on WHY18 by retaining only clusters with richness larger than 5. Table 2 reports the number of clusters selected for our analysis in each catalogue. It is worth to mention that not all the original cluster catalogues report the cluster mass in terms of ${\rm M_{500}}$. When a different mass is provided, we have converted it into ${\rm M_{500}}$, following the formalism described in the Appendix of Hu & Kravtsov (2003). The final redshift distribution of the cluster catalogues adopted in our analysis is shown in figure 2.

For each of the cluster catalogue, the cross-correlation analysis requires a proper mask that takes into account the portion of sky covered by the catalogue and possible systematic effects due to Galactic contamination or misidentification that derives from high stellar number densities. For the infrared catalogues WHY18 we adopt a mask derived from Alonso et al. (2015); for the optical catalogues SDSSDR9 the mask just refers to removing the portions of sky not covered by the SDSS survey; for the X-ray catalogues MCXCSub and HIFLUGCS the mask is a sharp cut for $|b|>20\degree$ with the additional removal of the Virgo region Schellenberger & Reiprich (2017). For all the catalogues, we assume uniform coverage, at least on the angular scales probed in our correlation analysis (i.e., below a few degrees).

3 Models

In this section we describe how we model the expected cross-correlation APS between the $\gamma$-ray sky and galaxy clusters at low redshift.

The first expected contribution for the cross APS is provided by $\gamma$-ray emissions from the center of clusters. This emission can arise from compact sources located inside the cluster like AGNs or from the innermost concentration of the ICM. Given the size of the Fermi-LAT PSF, such correlation between the $\gamma$-ray emission and clusters is well represented by correlation at zero angular separation (i.e., the emission comes from a region much smaller than the Fermi-LAT  PSF around the cluster center): this implies that the APS has a flat behaviour in multipole (also called shot-noise term), except for the correction due to the Fermi-LAT  beam window function:

[TABLE]

with $N_{c_{j}}$ being the number of objects in the cluster catalogue $j$ and $F^{\gamma}_{\Delta E_{i},k}$ the $\gamma$-ray flux coming from the center of the cluster $k$ in the energy bin $\Delta E_{i}$. The beam window function $W_{\ell}^{\Delta E}$ accounts for the finite angular resolution of the Fermi-LAT  instrument that suppresses high multipoles. It is computed from $W_{\ell}(E)=2\pi\int_{-1}^{1}d\cos\theta\,P_{\ell}(\cos\theta)\,{\rm PSF}(\theta,E)$ and then averaging over the energy bin by using an energy spectrum of index $-2.3$, which represents the mean spectral index of the UGRB (Ackermann et al., 2015).

Since gamma-ray sources typically have energy spectra that can be well approximated by a power-law, the energy dependence of the cross APS is modelled as a power law as well, and eq. 2 can be approximated as:

[TABLE]

where $\bar{E}_{i}=\sqrt{E_{i}\,E_{i+1}}$ is the geometric mean of the lower ($E_{i}$) and upper ($E_{i+1}$) bounds of the energy bin, $\Delta E_{i}=E_{i+1}-E_{i}$ is the size of the bin, $C_{P}^{j}$ is the normalization of the power-law relation and $\alpha_{1,j}$ is the spectral index. The index $j$ on the parameters refers to the fact that we fit the cross APS separately for each cluster catalogue.

In addition to the shot-noise term a correlation at larger angular separation angles might also be expected. We model this APS following the halo-model approach, where the correlation can be decomposed into the so-called 1-halo and a 2-halo terms (“1h” and “2h”, in the notation below). The former refers to the correlation between two points residing in the same physical halo, the latter to the case of two points belonging to two different halos. Correlations on angular scales larger than the Fermi-LAT  PSF can be due either to a 1-halo term associated to an extended ICM (or DM halo) or to a 2-halo contribution, involving e.g. two AGNs residing in different structures. For definiteness, since the modeling of $\gamma$-ray emission from ICM suffers of large uncertainties (Zandanel & Ando, 2014), we limit our modelling of the 2-halo term to the case of the more robust AGN emission.

The general expression defining the theoretical cross APS is:

[TABLE]

where $W_{c_{j}}$ and $W_{\gamma_{i}}$ are the window functions of the cluster catalogue $j$ and $\gamma$-ray source population $i$, $P_{c_{j}\gamma_{i}}$ is the cross-correlation 3D power spectrum, and $\chi$ is the comoving radial distance. To model these quantities we follow (Branchini et al., 2017).

The window function of cluster catalogues can be written as:

[TABLE]

with the number of objects in the cluster catalogue given by

[TABLE]

We empirically derive the cluster mass function from the catalogues themselves, considering the estimated redshift and $M_{500}$ masses provided by the cluster catalogues.

The window function of the $\gamma$-ray emission from a given source population is

[TABLE]

where $\mathcal{L}$ is the $\gamma$-ray rest-frame luminosity in the energy interval $0.1$ to $100\,\mathrm{GeV}$, $\Phi_{\rm S}$ is the $\gamma$-ray luminosity function of the source class $i$ of astrophysical emitters included in our analysis, and $dN_{\rm S}/dE$ is its observed (unabsorbed) energy spectrum. The upper bound $\mathcal{L}_{\rm max}(F_{\rm sens},z)={\rm min}[\mathcal{L}(F_{\rm sens},z),\hat{\mathcal{L}}_{\rm max}]$ with $F_{sens}$ being the flux above which an object is resolved in the FL8Y and 3FHL catalogues and consequently masked in our analysis. The minimum (maximum) intrinsic luminosity $\mathcal{L}_{\rm min}$ ($\hat{\mathcal{L}}_{\rm max}$) depends on the properties of the source class under consideration. For definiteness, we focus our analysis on misaligned AGNs, which are modeled as in Branchini et al. (2017).

Finally, the last ingredient we need is the three dimensional power spectrum, here decomposed in its 1-halo and 2-halo terms:

[TABLE]

[TABLE]

where $\langle f_{\gamma_{i}}\rangle$ is the mean luminosity density defined as $\langle f_{\gamma_{i}}\rangle=\int\ d\mathcal{L}\ \mathcal{L}\Phi_{i}(\mathcal{L},z)$ and $\bar{n}_{c_{j}}$ the mean cluster number density defined as $\bar{n}_{c_{j}}=\int\ dM\ \langle N_{c_{j}}dn/dM\rangle$. The halo mass function $dn/dM$ and bias $b_{h}$ are derived from Sheth & Tormen (1999), while the halo mass-luminosity relation is taken from Camera et al. (2015).

All these ingredients allow us to model a possible large-scale correlation, resulting in a non-flat dependence on the multipole angular scale. Although we use a model tuned on AGNs, we allow a certain level of flexibility in order to intercept possible large-scale deviations from the specific AGN emission. This is obtained by allowing a free power-law index for the energy dependence and a free overall normalisation parameter. This second model is therefore defined as:

[TABLE]

it contains the expected small-scale shot-noise term already introduced in eq. 3 to which we add the AGN-like model discussed above, with the fudging free parameters $\alpha_{2}$ and $A$. $C_{\ell}^{\Delta E_{0}}$ is the theoretical APS of eq. 4 calculated in a specific energy bin $\Delta E_{0}$, that we choose to be the $(1,2)$ GeV energy interval, with the spectral behaviour of the signal then carried by the free spectral index $\alpha_{2}$ and the size by the free normalisation parameter $A$. $\bar{E}_{i}$ is the geometric mean of the boundaries of the $i$ energy bin (reported in Table 1). $\bar{E}_{0}$ is the same mean energy for the reference $(1,2)$ GeV energy interval. While here we explicitly indicate the index $j$ (which labels the galaxy catalogue), for the sake of brevity in the rest of the paper we will omit the index $j$.

In conclusion, the FLAT model has two free parameters, the spectral index $\alpha_{1}$ and the normalisation factor $C_{P}$. It is embedded in the more complete physical AGN-like model, which is endowed with 2 additional free parameters, the spectral index $\alpha_{2}$ and the normalisation factor $A$ of the non-flat large-scale behaviour.

4 Cross-correlation signal

The goal of the analysis is to investigate the presence of a cross-correlation signal between $\gamma$-rays and low-redshift clusters, and then to study on which scales this signal originates, especially if a large-scale effect can be identified that could be related to the presence of a diffuse emission from the intra-cluster medium. The quantity we measure is the APS $C_{\ell}$ of the cross-correlation between the Fermi-LAT maps and the galaxy cluster catalogues described in Section 2. The APS is estimated by using Polspice (Chon et al., 2004), a tool to statistically analyse data pixeled on the sphere: it measures the two point auto- or cross-correlation APS, it is based on the fast spherical harmonic transforms allowed by iso-latitude pixelation such as HEALPix and corrects for the effects of the masks. The statistical method then adopted to quantify the presence of a signal against a null hypothesis is the Signal to Noise Ratio (SNR) defined as (see, e.g., Becker et al. (2016)):

[TABLE]

where $C_{\ell,\mathrm{DATA}}$ is the measured cross APS, Cℓ,MODEL is a model that grabs the physical features expected for the cross-correlation signal and $\Gamma^{-1}_{\ell,\ell^{\prime}}$ is the inverse of the cross APS covariance matrix. Theoretical models are described in the previous section, and the estimation of the covariance matrix is described below in a dedicated section. The model employed to assess the significance of the presence of a signal is $C_{\ell,\mathrm{MODEL}}=C_{\ell,\mathrm{FLAT}}$, where $C_{\ell,\mathrm{FLAT}}$ is defined in eq. 3, i.e., a model independent of the multipole $\ell$ (as expected from the shot-noise of a population of $\gamma$-ray sources) with an energy dependence similar to the one measured for the UGRB spectrum. The model parameters entering the computation of the SNR are determined as the best-fit parameters obtained as discussed below. The covariance matrix is carefully estimated using mocks, as described in detail in Section 5 and Appendix A.

As a second analysis we investigate whether a more refined and complete model for the cross APS is preferred over the simpler FLAT case. The model we use refers to the case where the $\gamma$-ray emission originates dominantly from AGNs. The model is outlined in Section 3 and its statistical preference over the FLAT case is determined by means of a $\chi^{2}$ analysis, where the $\chi^{2}$ function is defined as:

[TABLE]

where index $i$ denotes the energy bins and the sum over the multipole $\ell$ is limited to a range $(\ell_{\mathrm{min}},\ell_{\mathrm{max}})$ whose lower bound is chosen in order to exclude possible galactic-foreground residual contamination, while the upper bound is driven by the energy-dependent Fermi-LAT PSF, as discussed in Ackermann et al. (2018). The values for $\ell_{\mathrm{min}}$ and $\ell_{\mathrm{max}}$ for the different energy bins are reported in Table 1. The best fit FLAT and AGN models are determined by maximixing the Gaussian likelihood $\mathcal{L}$ defined as:

[TABLE]

The preference for the AGN model over the FLAT case is then determined by means the $\Delta\chi^{2}=\chi_{\mathrm{FLAT}}^{2}-\chi^{2}_{\mathrm{AGN}}$. This new quantity behaves approximately like a $\chi^{2}$ distribution with a number of degree of freedom (DOF) defined by the difference of the DOF between the two models that enter in the comparison. For the models used in our analysis and described in Section 3, we have DOF = 2.

We will also look at the Akaike information criterion (AIC), Akaike (1974), defined by the following expression:

[TABLE]

where $k$ is the number of parameters of the model and $\mathcal{L}$ is the likelihood. With this criterion one can estimate the relative quality of models in terms of information lost by the given model: the smaller is the loss of information, the better the model performs in reproducing the data.

5 Signal covariance through mocks

In order to determine the presence of a signal and its significance, we need a faithful determination of the covariance matrix, which is then used also to build a robust likelihood function adopted to constrain the models. To this aim, we investigated and compared different options and determined the optimal choice for our analysis. We performed an extensive and detailed investigation of the sources of covariance in the measurement of the cross APS between galaxy cluster catalogues and unresolved $\gamma$ ray emission. The main result we found is that the full covariance, along the two directions of multipole $\ell$ and energy, can be well approximated by a Gaussian estimate, diagonal in both dimensions, especially when the study is performed on binned (in $\ell$ and energy) data. While we retain in our analysis the full covariance for $\ell$, we can conclude that a Gaussian estimate of the error budget for this type of cross APS is a good approximation. We will substantiate these findings below. At the same time, we derived a general method under the approximation, valid in this case, that the cross-correlation contribution is smaller than the product of the auto-correlations of galaxies and gamma-rays, to provide and compare different estimates of the cross APS covariances, that can be implemented and adopted in any analysis of this kind. Details are provided in Appendices A and B.

The method adopted in this paper to build a reliable covariance matrix makes use of mock maps, obtained by randomising the true cluster and $\gamma$-ray data. The idea is to generate a large number of independent mock maps endowed with the same statistical properties of the original true maps. From these originated maps we then derive the covariance matrix in a direct way. The complete set of techniques investigated to produce mock maps, for both clusters and $\gamma$ rays, is described in Appendix A. For the final analysis, we decided to use the phase randomisation method to produce the mock $\gamma$ ray maps and the FLASK log-normal simulation for producing synthetic galaxy cluster catalogues. Even though the full description of these two methods is given in the Appendix, it is worth to provide some details here.

The concept behind the phase randomisation method is that every field defined on a sphere (like the $\gamma$ ray intensity maps measured by the Fermi-LAT) can be linearly decomposed in terms of spherical harmonics. The spherical harmonics are weighted by a set of coefficients from which the APS is defined. The APS of the map is invariant under rotations (i.e., under the phase shift on the coefficients described in eq. 29). This means that starting from the measured power spectrum of the true map, we can build mock maps by arbitrarily changing the phases of the spherical harmonic coefficients. All these synthetic maps conserve the APS, but they will give a non zero covariance.

For what concerns FLASK , all the information regarding the code and its implementation can be found in Xavier et al. (2016). FLASK  produces log-normal realisations of maps, starting from an input power spectrum. We feed FLASK with the measured APS of the cluster catalogues. All the synthetic maps generated by FLASK  possess the same power spectrum, but are otherwise randomised with respect to the original map that provides the input power spectrum. We stress that the code allows the user to adopt either a Gaussian or a log-normal probability distribution function: in our case we adopted a log-normal distribution, in order to include effects due to non-gaussianity in the covariance matrix.

To estimate the cross correlation covariance matrix, we produce 2000 realisations of the $\gamma$ maps in each of the 9 energy bin listed in Table 1 and 2000 mocks maps for each of the cluster catalogues. The large number of realisations is required to have numerical control on the off-diagonal terms of the covariance matrix: from our tests, we found that 2000 is a good compromise between statistics and computing time. We performed a large number of tests, which are summarised here in their main features, and discussed in the Appendix.

From these mock realisations, we can then construct the full covariance between different multipoles (for each energy bin, labelled by index $i$), obtained as:

[TABLE]

where $N_{\gamma}(N_{c})$ is the total number of $\gamma$-ray (cluster) mocks, $C_{\ell}^{a,c\gamma_{i}}$ is the APS measurement performed on the $a$-th mock realisations and:

[TABLE]

is the mean of the cross APS over all the realisations. In Appendix B we demonstrate that the covariance matrix can be obtained by averaging two estimates: $\Gamma_{\ell\ell^{\prime}}^{\hat{c}\gamma_{i}}$ obtained by correlating the $N_{\gamma}$ realisations of the $\gamma$ ray mocks with the true cluster map, and $\Gamma_{\ell\ell^{\prime}}^{c\hat{\gamma_{i}}}$ obtained by correlating the $N_{c}$ realisations of cluster mocks with the true $\gamma$ ray map:

[TABLE]

Each of the “half” covariance is obtained by using eq. 15. This technique allows us to reduce significantly the computing time, since we need $2N$ combinations instead of $N^{2}$.

A selection of the performed tests is shown in Fig. 3 and Fig. 4, where results are reported in terms of the binned covariance (the variance and covariance are binned over intervals of size $\Delta l=60$ with $\ell_{\textrm{min}}=40$ and $\ell_{\textrm{max}}=1000$ 555It is worth to mention the fact that binning the covariance matrix means to take a block sub-matrix in which the value in each block is given by the average over the covariance values in that block: this means that the diagonal of the binned covariance matrix includes some effect from the off-diagonal contribution of the full covariance matrix. ). Fig. 3 compares the results obtained through the mocks with the covariance estimate provided by Polspice and with the theoretical Gaussian prediction for a given energy bin $i$:

[TABLE]

where $C_{\ell}^{cc}$ and $C_{\ell}^{\gamma_{i}\gamma_{i}}$ denote the cluster and $\gamma$-ray autocorrelations, $f_{\rm sky,i}$ is the fraction of the sky probed by the survey in the energy bin $i$, $\Delta\ell$ is the multipole bin width and $\delta^{\rm K}_{ll^{\prime}}$ is the Kronecker symbol (the gaussian covariance is in fact diagonal, for which reason Fig. 3 shows the comparison for the variance at each multipole $\ell$). Polspice is a non minimum-variance estimator, while instead the theoretical estimate is not valid in presence of non-gaussianities, in which case it represents an underestimate of the true variance. For the estimation of the Polspice covariance matrix we refer to Efstathiou (2004) where the procedure to compute the so called pseudo-$C_{l}$ covariance matrix is described in detail. Here we just list the main steps. The code: (i) computes the $C_{l}$ from the auto-correlation function; (ii) corrects the auto-correlation function for the effect of the mask; (iii) computes the $C_{l}$ back; (iv) computes the $V$-matrix (see eq. 15a of Efstathiou (2004)); (v) estimates the final covariance matrix, corrected for the mask, beam and pixel effects, as in eq. 17 of Efstathiou (2004).

From Figure 3 we notice, as expected, that the Polspice variance over-estimates the gaussian prediction, as well as the variance from mocks. At the same time, we find that the variance obtained using the mocks is quite close to the Gaussian prediction. The differences between the two are of the order of $\sim 10\%$ at very large scales ($\ell<100$) for WHY18 and SDSSDR9, and always smaller than 10% for MCXCsub and HIFLUGCS; at smaller scales ($\ell>400$) they are slightly larger than $10\%$ for SDSSDR9, while smaller for all the other cluster catalogues; on intermediate scales ($100<\ell<300$) we observe differences between 5% and 10% for all the catalogues.

The “gaussianity” of the covariance matrix was not an expected result as it was not expected a large over-estimation of the variance from Polspice. We can confirm that the estimated covariance matrix is nearly (although not exactly) Gaussian by looking at the off-diagonal terms of the covariance matrix, that are small compared with the diagonal ones. In figure 4 we show the cross-correlation coefficient defined as:

[TABLE]

One can note that the off-diagonal terms of the mocks covariance matrix are always smaller than 5% with respect to the diagonal terms. The Polspice correlation coefficient, shown for comparison, is even more diagonal, as expected. Even though we obtain that the covariance matrix is significantly diagonal, nevertheless in our analyses we adopt the full (binned) covariance matrix.

Concerning the covariance between different energy bins, we again find that the cross APS between galaxy catalogues and $\gamma$ rays is quite diagonal, especially for energy bins of the size of those adopted in our analysis (reported in Table 1). This can be seen by evaluating the Gaussian estimator for the covariance in energy (at fixed multipole, for convenience):

[TABLE]

where $i$ and $j$ identify the different energy bins and by determining the corresponding correlation coefficient:

[TABLE]

which is analogous to the coefficient defined in Eq. 19 to investigate the off-diagonal terms of the covariance matrix for what concerns the multipole. Fig. 5 shows $r_{E}$ for some selected angular scales and for MCXCsub. For most of the angular scales, the off-diagonal elements of the covariance matrix are well below a 5% deviation from the diagonal elements. Only for smaller angular scales the effect reaches deviations of the order of 10%. Results are similar at different angular scales and for the other cluster catalogues.

6 Results

The measured cross-correlation APS between the galaxy clusters and the unresolved $\gamma$-ray intensity are shown in Fig. 6 and Fig. 7. The cross APS have been obtained by means of the Polspice estimator and the (co)variances have been derived as discussed in Section 5.

Fig. 6 shows a representative case of the the binned angular power spectrum $C_{\ell}^{c\gamma_{i}}$ as a function of the multipole $\ell$ (third energy bin of Table 1), for each of the four cluster catalogues. The error bars are the diagonal entries of the covariance matrix obtained from the mock analysis. Fig. 7 instead shows the energy dependence of the mean cross APS, defined as the average with respect to $\ell$ of the $C_{\ell}^{c\gamma_{i}}$ in each energy bin:

[TABLE]

where $\ell_{\rm min}$ and $\ell_{\rm max}$ are shown in Table 1 and $\Delta\ell=\ell_{\rm max}-\ell_{\rm min}$. The errors on $P_{E}$ are defined as:

[TABLE]

where $\sigma^{2}_{C_{\ell}^{c\gamma_{i}}}=\Gamma_{\ell\ell}^{c\gamma_{i}}$. To easy the visualization in the plot, the data of eq. (22) are multiplied by $E^{2.2}$ (expected behaviour of the UGRB (Ackermann et al., 2015)) and divided by $\Delta E$ (the width of the energy bin).

In order to determine the presence of a positive cross-correlation signal, we adopt the SNR defined in Eq. 11. Cℓ,MODEL is set at the best fit obtained with the featureless FLAT model defined in Section 3. We perform all our fits by employing a MCMC technique to determine the likelihood of Eq. 13. We specifically adopt a pure-Python implementation of Goodman and Weareas Affine Invariant Markov chain Monte Carlo Ensemble sampler (EMCEE) Foreman-Mackey et al. (2013). Once the best fit Cℓ,MODEL model is obtained, we determine the SNR, whose results are reported in Table 2. The SNR analysis shows that the clusters in the WHY18 and SDSSDR9 catalogues exhibit a mild preference for a positive cross-correlation signal, while those in MCXCSub and HIFLUGCS provide a larger SNR, in excess of 3. Therefore, although not large, an evidence of $\gamma$-ray emission from those clusters appears to be present.

In order to look for a possible large scale contribution, we then fit the measured cross APS by adopting this time a physical model which follows the features of an AGN-like $\gamma$-ray emission, as described in Section 3. This model, in fact, possesses a large-scale 2-halo term. We test whether the AGN-like model is preferred over the FLAT model by mean of a $\Delta\chi^{2}=\chi_{\mathrm{FLAT}}^{2}-\chi_{\mathrm{AGN}}^{2}$ test. Table 2 shows that for MCXCsub a large-scale contribution is preferred at the $2.1\sigma$ level, whilst the other catalogues do not show a preference for the $\mathrm{AGN}$ model over the $\mathrm{FLAT}$ one.

Thus MCXCsub, with a SNR of 3.5 and a $\Delta\chi^{2}=4.5$ pointing to a large-scale contribution, turns out to be the most interesting catalogue. We would expect HIFLUGCS and MCXCsub samples to provide similar results, since the clusters in the two catalogues share similar mass function and flux distributions. On the other hand, a statistical difference of less than $2\sigma$ between two different samples of (possibly) the same population is nevertheless plausible.

The AIC test confirms the preference for a large-scale contribution in MCXCsub: in particular $\mathrm{AIC}_{\mathrm{FLAT}}=152.05$ against $\mathrm{AIC}_{\mathrm{AGN}}=151.35$ is the indication that the AGN model allows for a smaller loss of information, even though with a somewhat smaller evidence than in the $\Delta\chi^{2}$ analysis. All the other three catalogues show instead a preference for the FLAT model, namely $\mathrm{AIC}_{\mathrm{FLAT}}<\mathrm{AIC}_{\mathrm{AGN}}$ .

In Fig. 8, we show the allowed regions obtained for the model parameters in the MCXCsub case. The contours and shaded areas refer to the 2-dimensional allowed regions at 68% (dark blue) and 95% (light blue) confidence levels. From Fig. 8 we see that the constraints are somewhat loose, but consistent with an AGN-like model, with the spectral indices close the mean blazar $\gamma$ ray emission (The Fermi-LAT collaboration (2019b)), but notably somewhat smaller. A spectral index lower than $\sim 2$ is indicative of a hardening of the spectrum for the unresolved population of blazars, a result compatible with the findings of Ref. (Ackermann et al., 2018). The AGN normalisation, instead, turns out to be much larger than expected, even though with a sizeable uncertainty. The model adopted in our analysis is normalised such that the integral of the window function over the redshift provides approximately the measured UGRB intensity (see (Ammazzalorso et al., 2018) for details). The value of the normalisation we obtain here from the MCMC is $A=71.5^{+19.9}_{-29.9}$. We verified that such conclusion is independent on the specific model of AGN or blazars adopted in eq. 10.

A value this large could therefore exceed significantly the UGRB intensity, if due to a 2-halo emission: in the halo model, it is the 2-halo term which is directly related to the total $\gamma$ ray emission, while instead the 1-halo term can be large without necessarily inducing an exceedingly large total emission. Indeed, both the gamma-ray intensity and the two halo term are set by the window function. This can be seen from their definitions: $I_{\gamma}=\int d\chi\,W_{\gamma}(\chi)$, and $C_{\ell}^{c\gamma,2h}=\int d\chi/\chi^{2}\,W_{c}(\chi)\,W_{\gamma}(\chi)\langle b_{c}(\chi)\rangle\,\langle b_{\gamma}(\chi)\rangle\,P^{lin}(k=\ell/\chi,\chi)$, where $\langle b_{i}\rangle$ is the bias of source $i$ with respect to matter. For a generic population emitting gamma-rays, $\langle b_{\gamma}\rangle\sim 1$ at low-$z$ (i.e., in the range we are considering). Therefore a normalization different from one for the cross-correlation APS could only be re-absorbed in the window function, thus affecting the intensity in the same way. For what concerns the one-halo term, there is instead an additional ingredient, that is poorly constrained, and can significantly change the strength of the correlation without affecting $W_{\gamma}$ and in turn $I_{\gamma}$, which is the average gamma-ray luminosity from a cluster of a given mass and redshift. Making this function steeply increasing with the cluster mass can boost the one-halo term.

The result we found might thus indicate that the model we implemented is able to effectively capture a large-scale contribution, but such correlation is not due to a 2-halo term involving $\gamma$-rays from AGNs (or other galactic sources) in 2 different halos at large physical distances. On the contrary it might be seen as a potential indication in favour of a diffuse emission from the ICM (which would instead be a 1-halo term). Indeed, a relevant 1-halo term providing correlation on scales around $0.5-1$ degree and provided by $\gamma$-rays from the ICM can be obtained (Branchini et al., 2017; Reiss et al., 2018), with no obvious violation of other existing bounds.

Clearly, if such a signal is present, it must be provided by the clusters with a size of their diffuse emission significantly larger than the Fermi-LAT  PSF. The 1-halo signal from the clusters with an angular dimension below/around the Fermi-LAT  PSF would instead be described by a featureless APS, like in the FLAT case. In order to investigate more deeply this issue, we subdivided the MCXCsub cluster catalogues in two sub-catalogues by looking at the angular size of the clusters. The selection is done according to the angular dimension of the clusters $\theta_{500}=r_{500}/d_{A}(z)$, where $r_{500}$ is the virial radius relative to an overdensity of 500 (as defined in footnote 4) and $d_{A}(z)$ is the angular diameter distance, which depends on redshift $z$. We compute the average $\bar{\theta}_{500}$ over all the MCXCsub clusters and we focus this new analysis on all those clusters with $\theta_{500}>\bar{\theta}_{500}$. They are expected to be the main contributors to the extended 1-halo correlation if the hypothesis of ICM emission is correct. We have that $\bar{\theta}_{500}=0.267$ and the total number of clusters with angular size smaller or larger than $\bar{\theta}_{500}$ are 72 and 37, respectively. The angular size distribution of the MCXCsub is shown in Figure 9. We can see that most of the MCXCsub clusters have a size larger than the Fermi-LAT PSF in most energy bins. The latter is reported in Table 1. The vertical solid line in the figure refers to an angle of 0.2 deg, which is an approximate illustration of the Fermi-LAT beam.

We then perform the fit of only the MCXCsub which are larger than the average $\bar{\theta}_{500}$. This analysis requires to build a new set of cluster mocks (2000) produced in the same way as discussed above, from which we determine the APS covariance matrix. The best fit results are shown in 10 and the results for the $\Delta\chi^{2}_{AGN-FLAT}$ turns out to be only 2.1, which corresponds to $0.96\sigma$. Contrary to the expectations for an ICM signal, the statistical significance does not increase when only the largest clusters are considered: instead, it rather decreases as compared to the full MCXCsub sample. This significantly weakens a possible ICM interpretation, and leaves open the alternative between an unresolved blazar population (with a slightly harder energy spectrum as compared to the resolved ones) and a diffuse emission from the cluster itself, like in the case of the intra-cluster medium emission. The value of the best fit parameters are similar to what is found for the full MCXCSub case, with a rather large normalisation parameter: $A=65.0^{+24.8}_{-35.7}$. This time, the parameter is consistent at $1.8\sigma$ with an interpretation in terms of LSS from AGN emission, although the large uncertainty does not allow to make firm conclusions.

For illustrative purposes, in order to visualize the angular scales at which this excess occurs, we show the cross correlation function (CCF) also in configuration space. The CCF for the subset of MCXCsub including most extended clusters (those with $\theta_{500}>\bar{\theta}_{500}$) is reported in Fig. 11. The plot refers to the energy range 1-10 GeV, where the photon count statistics is large and the angular resolution not too poor. The grey area indicates the size of Fermi-LAT PSF.

The CCF exhibits a significant ”noise” term at small angular scales, compatible with unresolved AGN point-like emission. A peak is present on angular scales of the order of $0.7\degree$, although not so statistically significant to determine a preference for a large-scale term. This excess occurs at similar scales as those obtained in Ref. Reiss & Keshet (2018) by means of a stacking analysis of the $\gamma$-ray emission around galaxy clusters and is there interpreted as due to the presence of virial shocks in the clusters.

7 Conclusions

We analysed the cross-correlation angular power spectrum between the unresolved extra-galactic $\gamma$-ray background measured by the Fermi-LAT and the large scale structure of the Universe at low redshift traced by four galaxy clusters identified in three different bands: WHY18 (infrared band), SDSSDR9 (optical band), MCXCsub and HIFLUGCS (X-ray band). The main motivation was to investigate whether the cross-correlation technique could identify the presence of an extended $\gamma$-ray emission possibly compatible with an intra-cluster medium emission.

For all the four catalogues, the analysis confirmed that the unresolved $\gamma$-ray emission observed by Fermi-LAT correlates with the large-scale clustering in the Universe as observed by Fornengo et al. (2015) with the LSS tracer given by the CMB lensing, by Xia et al. (2015) with galaxies and by Branchini et al. (2017) specifically with clusters. We found that the largest significance occurs for the galaxy clusters identified in the X-ray band, i.e. MCXCsub and HIFLUGCS, for which the SNR is 3.5 and 3.2, respectively. When compared with a theoretical model which contains an explicit term referring to a large-scale $\gamma$-ray emission, MCXCsub exhibits a clear preference for this type of emission as compared to the a model containing only “shot-noise” emission from unresolved point sources, like sub-threshold AGNs. The energy spectrum of this latter component is found to be slightly harder that the mean spectral behaviour of resolved blazars, possibly indicating differences between the resolved and unresolved components of the AGN population, as observed also in (Ackermann et al., 2018). Further investigation of the extended emission could not disentangle between the two options offered by a large-scale 1-halo term, possibly linked to an intra-cluster medium emission, and a large-scale 2-halo contribution, like it would occur if the correlation is due to the large-scale distribution of point-like AGNs.

However, the analysis in angular space shows a peak in the correlation function on angular scales of the order of $0.7\degree$, which appears compatible with the results of Reiss & Keshet (2018) obtained by means of a stacking analysis and where the peak is associated to the $\gamma$-ray emission in virial shocks. In our analysis, we confirm the presence of a fluctuation on similar angular scales, although we do not have the sensitivity to determine whether the peak in the correlation function has a physical origin or instead just reflects a statistical fluctuation.

In developing the analysis, we also derived and tested technical tools specifically designed to determine reliable multidimensional covariance matrices, which are a key ingredient for the study of cross correlation signals. These methods are summarised in the Appendices and refer to development, test and comparison of general numerical techniques for the massive and efficient production of mock realisations of the sky for cross-correlation studies, like the correlation of catalogues of galaxies or clusters with $\gamma$-ray maps. We then developed a semi-analytic framework that allowed us to properly join the information coming from the two pieces of the covariance matrix (galaxies/clusters catalogues from one side, and $\gamma$-rays from the other side) without overestimating the error matrix: this is clearly important for the estimation of the significance of the presence of a signal and for the inference of model parameters. These techniques are general enough such that they can be used for any distribution of objects and can be adapted easily to different statistical and astrophysical analyses.

Acknowledgements

We warmly thank Enzo Branchini for very useful and interesting discussions. This work is supported by the following grants: Departments of Excellence (L. 232/2016), awarded by the Italian Ministry of Education, University and Research (MIUR); The Anisotropic Dark Universe, Number CSTO161409, funded by Compagnia di Sanpaolo and University of Torino; TAsP (Theoretical Astroparticle Physics) project, funded by the Istituto Nazionale di Fisica Nucleare (INFN); PRIN 2017 project (Progetti di ricerca di Rilevante Interesse Nazionale) The Dark Universe: A Synergic Multimessenger Approach, Number 2017X7X85K, funded by MIUR. MR acknowledges support by “Deciphering the high-energy sky via cross correlation” funded by Accordo Attuativo ASI-INAF n. 2017-14-H.0. J.-Q.X. is supported by the National Science Foundation of China under grant No.11633001 and No.11690023.

Appendix

As pointed out in the main section of the paper, an accurate estimate of the covariance matrix of the cross-correlation angular power spectrum along its different dimensions (multipole, energy, redshift) is necessary to infer the statistical significance of our results. The method we are using to derive the cross-correlation APS is based on Polspice, which is a non-minimum variance algorithm. We therefore investigated methods based on the production of mock realisations of the two maps which enter the cross-correlation analyses, from which we derive estimates of the covariances of the cross-APS, in order to determine which combination for the generation of the two maps is more suitable for cross-correlation analyses, and to verify if the results are stable across different techniques for the generation of mocks. The method we devise starts from the true maps under investigation, from which the mocks are produced, and is quite general: this makes it directly adaptable to derive the covariance also along those directions (like energy and redshift) for which Polspice cannot be used. Since cross-correlations deal with two observables, the generation of a suitably large number $N$ of maps for each set of observables implies the crossing of $N\times N$ maps: this can make the computation of the covariance quite demanding from the computational point of view, since $N$ needs to be large enough to make all the procedure reliable. We therefore devised, tested and validated a method which allows the computing resources to grow linearly with $N$ rather than quadratically: this is obtained by deriving two estimates of the covariance by correlating the $N$ mocks of the first observable with the true map of the second observable, and then performing the opposite: we demonstrate below that the total covariance can then just be obtained as the average of these two “half” covariances. The results thus obtained are a faithful estimate of the global $N\times N$ crossing.

In the reminder of this Appendix we discuss the various techniques adopted to generate mock maps for both galaxy/cluster catalogs and for $\gamma$ ray maps. In Appendix B we derive a method which allows to drastically reduce the computing power necessary for the mock analysis and in Appendix C we show results of the methods, applied to some specific cases based, for definiteness, on the 2MPZ galaxy catalog and the Fermi-LAT $\gamma$ ray maps.

Appendix A Generation of mocks

The methods we use to generate mock maps starting from a tue map are the following:

• •

Bootstrap

• •

Jackknife

• •

Phase Randomization

• •

Gaussian Realizations (synfast)

• •

Lognormal Realizations (Flask)

We can group these five methods in two categories: resampling procedures (bootstrap and jackknife), that allow us to build mocks just with a reorganisation of the original data sample (galaxes/clusters or $\gamma$-ray map); generated fields procedures (phase randomization, synfast and Flask), that use the statistical distribution of the original data sample to build mocks.

We pre-process our data sets (either galaxies/clusters or $\gamma$ rays) in HEALPix format: this allows us to adopt the same procedure for both type of observables. Galaxies/cluster maps are produced in terms of number counts per pixels, $\gamma$ rays maps in terms of photon intensity per pixel. We adopt a HEALPix pixelation format with resolution parameter $N_{\rm side}=1024$, which correspond to a total number of pixels $N_{\rm pix}\,=12,582,912$ and mean spacing of $\sim 0.06^{\circ}$. Being interested in the unresolved component of the $\gamma$ ray emission, we apply masks for resolved point sources and galactic emission, as described in Section 2.1. Masks may apply also to the galaxies/clusters catalogs. The presence of (typically quite different) masks for the two fields make the determination of the covariance matrix quite complex and involved, for which the mock techniques becomes especially useful.

For each method used to produce the mock realisations of both the galaxies/clusters and $\gamma$ ray maps, we test that the ensuing auto-correlation APS is recovered from the mocks: we measure the auto APS for each mock map and verify that the average of these APS recover the APS of the corresponding data maps.

A.1 Bootstrap

A map in HEALPix format is an array of pixels where each element represents the intensity of the specific pixel. To make one bootstrap realisation we follow these steps:

• •

We divided the full array in $N_{\rm sub}\,$ sub-arrays, so that each of them has $N_{\rm pix}\,/N_{\rm sub}\,$ pixels;

• •

We label each of the sub-array;

• •

We randomly pick $N_{\rm sub}\,$ sub-arrays with replacements and form a new resampled HEALPix map

Reiterating these three steps would produce $N_{\rm r}\,$ bootstrap realisations. The new HEALPix map originated in this way is characterised by the same number of pixels of the original data set, but each of the sub-array can be selected more than one times or not selected at all; for this reason we have to weight each sub-array by the number of times it is selected. This procedures is called bootstrap with replacements (Norberg et al., 2009).

In this case the estimator of the APS covariance matrix is given by:

[TABLE]

where $x_{i}$ is the i-th bootstrap realisation and $N_{\rm r}\,$ is the number of realisations and,

[TABLE]

A.2 Jackknife

As for the bootstrap technique, also for the jackknife method the map is divided in $N_{\rm sub}\,$ sub-arrays with the same number of pixels. Each of the sub-arrays is labelled, but in this case a realisation is obtained by systematically omitting one of the sub-array in each realisation. The resampling of the data-set consists of $N_{\rm sub}\,-1$ remaning sub-arrays with volume $(N_{\rm sub}\,-1)/N_{\rm sub}\,$ times the volume of the original data-set (Norberg et al., 2009). By definition there are only $N_{\rm r}\,=N_{\rm sub}\,$ different copies of the data set that are created in this way. In this case the APS covariance matrix estimator reads:

[TABLE]

where $\bar{x}_{i}$ is given by eq. (25). The factor $(N_{\rm r}\,-1)$ accounts for the lack of independence between the $N_{\rm r}\,$ copies of the data set.

A.3 Phase randomization

The implementation of this procedure is described in De Domenico & Lyberis (2012). It based on the fact that it is always possible to write an intensity map as a linear combination of spherical harmonics:

[TABLE]

from which the angular power spectrum is obtained:

[TABLE]

It is clear that Eq. (28) is invariant under a phase rotation on the harmonic amplitudes $a_{lm}$:

[TABLE]

Taking advantage of this symmetry, we can build independent realisations of the initial intensity map, each sharing the same APS. Since we determine the true-map APS from a masked sky, we need to correct for it, in order to produce a mock map that contains the correct statistical properties of the original map. The procedure we adopt is:

• •

Measure the auto APS ($C_{\ell}^{\gamma\gamma}$ or $C_{\ell}^{gg}$ for $\gamma$ rays or galaxies/clusters) from the masked data maps

• •

Transform: $a_{\ell m}\longrightarrow\tilde{a}_{\ell m}=a_{\ell m}e^{i\varphi_{\ell m}}$

• •

Construct a full-sky mock map: $\tilde{f}(\theta,\phi)=\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell}\ \tilde{a}_{\ell m}\ Y_{\ell m}(\theta,\phi)$

• •

Correct the mock map for incomplete sky: $\tilde{f}(\theta,\phi)\longrightarrow\tilde{f}(\theta,\phi)\ \times W(\theta,\phi)\ \times f^{-1/2}_{\rm sky}$

The $(f_{\rm sky})^{-1/2}$ accounts for the fact that the original map was masked and therefore the obtained harmonic amplitudes has reduced power as compared to the true one. $W(\theta,\phi)$ restores the mask on the mock map. Let us notice that this method looses information on the shot-noise, and therefore it can produce underestimate of the covariance in situations where the shot-noise is large.

The evaluation of the APS covariance matrix is finally done with eq. (26)

A.4 Gaussian realisations (Synfast)

Synfast666https://healpix.jpl.nasa.gov/html/facilitiesnode14.htm is a HEALPix routine that allows to generate realisations of a Gaussian random fields on a sphere, starting from an input APS. The procedure We therefore start from the APS describing the statistical distribution of the data sample we want to replicate:

• •

Measure the auto APS ($C_{l}^{\gamma\gamma}$ or $C_{l}^{gg}$ for $\gamma$ rays or galaxies/clusters) from the masked data maps

• •

The obtained APS is fed to synfast, which outputs a full-sky mock map: $\tilde{f}(\theta,\phi)$

• •

Mask the mock map: $\tilde{f}(\theta,\phi)\longrightarrow\tilde{f}(\theta,\phi)\ \times W(\theta,\phi)$

The evaluation of the APS covariance matrix is finally done with eq. (26)

A.5 Lognormal realisations (FLASK)

Flask777http://www.astro.iag.usp.br/~flask is a C++ code, parallelised with OpenMP, based on the work of Xavier et al. (2016) and created to generate mock realisations of galaxy distributions starting from their 3D power spectrum. Like synfast, it generates multiple correlated fields on spherical shells, after providing the power spectrum describing the distribution to be replicated. Differently from synfast, the generated maps are obtained from a lognormal distribution. The tomographic approach used by Flask slices the three dimensional space into spherical shells (redshift slices), each one discretized in Healpix maps. After generating the fields, Flask can apply selection functions and noise to them. The output can be in the form of a source catalogue and/or Healpix maps, among others.

We use Flask to generate $N_{\rm r}\,$ independent realizations. The evaluation of the APS covariance matrix is finally done with eq. (26). Although the code is thought to work for galaxy distributions in different redshift bins, we tried to use it also for $\gamma$ ray maps, by using the APS instead of the 3D power spectrum.

A.6 Covariance estimators’ relations

As we have shown in the previous sections, we can set a different APS covariance matrix estimator for each of the methods we use to produce mocks. It can be useful to have a look to the relation between the different estimators.

Given a data vector $\textbf{X}=\{x_{1},x_{2},...,x_{N_{\rm r}\,}\}$ we can write the generic expression for the covariance matrix of X as:

[TABLE]

with $\bar{x}_{i}$ is given by eq. (25). Eq. (30) can also be rewritten as:

[TABLE]

The unbiased definition of the sample covariance when the mean is derived from the sample itself is:

[TABLE]

We can relate the estimator of eq. (LABEL:eqn:covU) to the ones obtained with the jackknife (26) and bootstrap (24) techniques as:

[TABLE]

Appendix B Semi-analytic prediction of the cross-correlation covariance

The derivation of the covariance matrix for the cross-correlation APS requires to combine the information arising from the generation of many different set of maps. In order to obtain stable results, the required number of realisations for each of the two observables can be large (in our analysis on the cross-correlations between clusters and $\gamma$ rays, we produced 2000 mocks for each set), and can require to produce maps in several energy bins (for $\gamma$ rays) and redshift bins (for galaxies/clusters). In this section we show that we can obtain a reliable estimate of the full covariance matrix by performing a simpler combination, namely we can construct two partial estimates of the covariance matrix by combining separately: (i) the galaxies/clusters mock with the measured $\gamma$ ray map; (ii) the $\gamma$ ray mocks with the measured galaxies/clusters maps. The final estimate of the covariance is obtained as the average of these two “half” covariances. This reduces the number of combinations from $(N_{\rm r}\,)^{2}\times(n_{E}\times n_{z})$ to $2(N_{\rm r}\,\times n_{E}\times n_{z})$ (where $N_{\rm r}\,$ denotes the number of mock maps produced for each of the $n_{E}$ energy bins and $n_{z}$ redshift bins), which is $N_{\rm r}\,/2$ times faster. We show that this approach is correct in the limit of a large number $(N_{\rm r}\,)$ of realisations, by following a theoretical derivation based on the Gaussian prediction for the APS covariance matrix. We have numerically verified that the result shown here hold in a more general situation where gaussianity is not necessarily present. The method we are going to describe is valid in the case where the cross-correlation term is small when compared to the product of the auto-correlations for galaxies and gamma-rays.

Let us start with the Gaussian prediction for the APS covariance matrix: (Hu & Jain, 2004):

[TABLE]

where $\Delta l$ is the bin width, $f_{\rm sky}$ is the fraction of the sky probed by the surveys and $\delta^{\rm K}$ is the Kronecker symbol.

Let us denote with a hat symbol quantities which are measured on the real maps, while quantities obtained from mocks do not have the hat symbol. For instance, $\hat{C}_{l}^{gg}$ and $\hat{C}_{l}^{\gamma\gamma}$ are the galaxies/clusters the $\gamma$ ray auto-correlation APS measured on the true data maps. Let us define a covariance term obtained by computing the cross-correlation between the real galaxy distribution and the mock $\gamma$ rays realisations. From Eq. (35) and considering that we construct the covariance from the mock by averaging over the $N_{\rm r}\,$ realisations:

[TABLE]

Let us then define the corresponding counterpart term obtained by using the mocks for galaxies/clusters and data for $\gamma$ rays:

[TABLE]

When instead we use only mocks, Eq. (35) gives:

[TABLE]

If we take the average of the expressions 36 and 37, we obtain:

[TABLE]

Since the measurements of the APS ($\hat{C}_{l}^{gg}$ and $\hat{C}_{l}^{\gamma\gamma}$) obtained using the real map are well reproduced by the APS measurements from mock maps:

[TABLE]

we obtain:

[TABLE]

We then observe that in the limit of large $N_{\rm r}\,$:

[TABLE]

In this case, Eq. (38) and Eq. (41) give the same result. Therefore we can obtain a reliable estimate the covariance by simply averaging over the sum of the two “half” contributions:

[TABLE]

Appendix C Comparison of the methods to produce map mocks

In this Appendix we report some results for the determination of the cross-correlation covariance matrix obtained using the methods described in Appendix A and using the relation shown in Eq. 43 to join the separate information coming respectively from the $\gamma$ and galaxy mock maps. For definiteness, in this analysis we use the the 2MASS Photometric Redshift catalogue (2MPZ, Bilicki et al. (2014)) that is a galaxy catalogue built by cross-matching 2MASS XSC, WISE and SuperCOSMOS all-sky samples with galaxy photometric redshift reconstructed via an artificial neural network. The employed algorithm is the one described in Collister & Lahav (2004) and trained on several redshift surveys (2MRS, SDSS, 6dFGS, 2dFGRS and ZCAT). The all-sky accuracy of the redshift reconstructed by the network is close to $\sigma_{z}=0.015$ for nearly all the dataset with few outliers. The resulting 2MPZ sample contains almost 1 million galaxies with a median redshift of $z=0.07$. In the left panel of fig. 12 we show the 2MPZ catalogue in HEALPix projection and in the right panel we show its redshift distribution.

We use this catalogue to test our methods to exploit the large statistics in terms of number of galaxies (934,175), that allow us to reduce the galaxy shot-noise contribution.

We produce 256 galaxy mock maps starting from the 2MPZ galaxy catalogue and 256 $\gamma$-ray maps in each energy bins used also in the main text analysis (see Table 1) using the five methods described in Appendix A. With this relative small number of mocks, per energy bin, we can test the accuracy in the estimation of the diagonal of the cross-correlation covariance matrix and compare it with its Gaussian prediction of Eq. 35). For the non diagonal terms, we need to produce a much larger number of mocks, and in the analysis of the main text we use 2000 realizations for each observable.

The results are shown in figures 13-17. In each figure we show the variance estimated using Eq. 41 (coloured lines) for a low (left), intermediate (center) and large (right) energy bin; in all plots, the black line represents the Gaussian prediction. Each plot is accompanied by a lower panels, where we show the ratio between the result for each combination and the Gaussian variance.

In Fig. 13 the method to produce $\gamma$ ray mocks is fixed to phase randomisation, and the galaxy mocks methods are rotated among the five options discussed in Appendix A. We notice that phaseGAM+bootstrapGAL tend to systematically underestimate the covariance, producing results even smaller by $\sim 30\%$ than the Gaussian prediction. This happens for all energy bins. The combination phaseGAM+phaseGAL slightly underestimante by a few percent the Gaussian variance in the first energy bin, but it is almost Gaussian in the other cases; the other combinations produce estimates in excess of the the Gaussian prediction of about $\sim 30\%$ for phaseGAM+phaseGAL/flaskGAL and about $60\%$ for phaseGAM+jackknifeGAL.

A similar behaviour is observed in Figure 14, where synfast is used to produce the $\gamma$ ray mocks. All the combinations show almost the same trend observed in Fig. 13. We expected similar results between the phase randomisation and synfast technique, since the two technique are quite similarly implemented in the generation of mock maps.

Figure 15 shows the result for the bootstrap method applied to $\gamma$ rays, in combination with all methods for the galaxy mocks. In this case we observe that all the combinations except bootstrapGAM+jackknifeGAL underestimate the Gaussian prediction of about $\sim 30\%$ in the case of bootstrapGAM+synfastGAL/flaskGAL, and by more than $50\%$ for bootstrapGAM+phaseGAL/bootstrapGAL.

The jackknife applied to $\gamma$ rays in combination with all methods for galaxies is shown in Fig. 16. In this case, all the combinations largely overestimate the gaussian prediction with just jackknifeGAM+bootstrapGAL/phaseGAL remaining below a $50\%$ difference.

Finally, we show in figure 17, the adoption of Flask to produce $\gamma$ ray mocks: it is clear that all the combinations show an inconsistent behaviour, with peculiar fluctuations for multipole scales smaller than 800. This behaviour is likely due to the fact that Flask is built to reproduce a galaxy distribution and is not general enough to be used for $\gamma$ ray maps. We decided to test its use also to generate $\gamma$ ray mocks: although the APS is well reproduced, the behaviour of the covariance does not give results which look trustable, especially when compared with all the other methods shown above.

In conclusion, from the extensive analysis of the different combinations, we found that FlaskGAL + PhaseGAM represents a good options for estimating the covariance matrix for the cross-correlation APS of galaxies/clusters with $\gamma$ rays. This combination produces covariance in slight excess of the gaussian prediction for almost all situation tested.