TY - JOUR
AU1 - Colavincenzo,, Manuel
AU2 - Tan,, Xiuhui
AU3 - Ammazzalorso,, Simone
AU4 - Camera,, Stefano
AU5 - Regis,, Marco
AU6 - Xia,, Jun-Qing
AU7 - Fornengo,, Nicolao
AB - ABSTRACT We report the identification of a positive cross-correlation signal between the unresolved gamma-ray emission, measured by the 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 realizations 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. cosmology: observations, cosmology: theory, large-scale structure of Universe, gamma-rays: diffuse background 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 π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 (AGNs), are present. These astrophysical sources, which 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 haloes (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; Brunetti, Zimmer & Zandanel 2017; Lisanti et al. 2018; Reiss, Mushkin & Keshet 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, Benoit-Lévy & Komatsu 2014; Shirasaki, Horiuchi & Yoshida 2014; Cuoco et al. 2015; Fornengo et al. 2015; Regis et al. 2015; Shirasaki, Horiuchi & Yoshida 2015; Ando & Ishiwata 2016; Shirasaki et al. 2016; Feng, Cooray & Keating 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 realizations 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 yr 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 yr 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 extragalactic gamma-ray background (UGRB). The angular resolution of the instrument is energy dependent and reaches ∼ 0.1° above 10 GeV. The photon and exposure maps adopted in the present analysis are produced with the Fermi Tools.1 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 class,2 which 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 /\ \mathit{ N}_{\rm pix}$|⁠. We adopt a HEALPix (Gorski et al. 2005) pixellation format with resolution parameter |$N_{\rm side} =$| 1024, which corresponds to a total number of pixel |${N}_{\rm pix}$| = 12582 912 and a mean spacing of ∼0.06°, 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 microbins were subsequently re-binned into nine 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). Table 1. Energy bins in GeV used in our analysis. Emin and Emax denote the lower and upper bound of the bins, while ℓmin and ℓmax show the interval in multipole ℓ 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 percent containment angle θcont is reported. Bin . Emin (GeV) . Emax (GeV) . ℓmin . ℓmax . θ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 Bin . Emin (GeV) . Emax (GeV) . ℓmin . ℓmax . θ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 Open in new tab Table 1. Energy bins in GeV used in our analysis. Emin and Emax denote the lower and upper bound of the bins, while ℓmin and ℓmax show the interval in multipole ℓ 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 percent containment angle θcont is reported. Bin . Emin (GeV) . Emax (GeV) . ℓmin . ℓmax . θ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 Bin . Emin (GeV) . Emax (GeV) . ℓmin . ℓmax . θ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 Open in new tab Since there is no clear sign of resolved extended emission from clusters, we adopt a cross-correlation technique which focuses 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 extragalactic 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 disc emission is masked by means of a latitude cut that excludes the portion of the sky with |b| < 30°. For point sources, the masking radius around each catalogues sources is defined as: $$\begin{eqnarray} F_{\Delta E}^\gamma \, \exp {\left(-\frac{R^2}{2 \theta _{\Delta E}^2}\right)} \gt \frac{F_{\Delta E,\rm faintest}^\gamma }{5}\, , \end{eqnarray}$$ (1) where |$F_{\Delta E}^\gamma$| is integral flux of the source in the specific energy bin ΔE under consideration, |$F_{\Delta E,\rm faintest}^\gamma$| is the flux of the faintest source in the same energy bin, and θΔE is the 68 per cent 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.fits.3 The foreground template is projected in HEALPix maps keeping the same Nside of the intensity maps and the same microbinning. A free normalization 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-fitting normalizations obtained with this procedure are all well compatible with unity, showing a successful description of the foreground emission. We normalize the foreground templates with the obtained parameters and then subtract the foreground maps from the corresponding intensity maps: this procedure is performed after having rebinned them in the macroenergy 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 Fig. 1, we show as an example the gamma-ray masked map in the (2.3−4.8) GeV energy bin before (left-hand panel) and after (right-hand panel) the galactic foreground removal. Figure 1. Open in new tabDownload slide Gamma-ray maps in the (2.3−4.8) GeV energy bin, masked with the procedure described in the text. The panels show the flux map before (left-hand panel) and after (right-hand panel) galactic-foreground subtraction. For illustration purposes, the maps have been rescaled to Nside = 128 and smoothed with a Gaussian beam of size σ = 0.4°. Figure 1. Open in new tabDownload slide Gamma-ray maps in the (2.3−4.8) GeV energy bin, masked with the procedure described in the text. The panels show the flux map before (left-hand panel) and after (right-hand panel) galactic-foreground subtraction. For illustration purposes, the maps have been rescaled to Nside = 128 and smoothed with a Gaussian beam of size σ = 0.4°. 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, Irwin & MacGillivray 2001). The main selection applied to the sources is such as to include clusters with |${\rm \mathit{ M}_{500}}\gt 3\times 10^{14}\rm M_\odot$|4 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 1743 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 \mathit{ M}_{500}} \gt 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 \mathit{ M}_{500}} \gt 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 M500. When a different mass is provided, we have converted it into M500, 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 Fig. 2. Figure 2. Open in new tabDownload slide Redshift distribution of the galaxy cluster catalogues used in our analysis: WHY18 (limited to richness larger than 9), SDSSDR9, MCXCsub, and HIFLUGCS. Figure 2. Open in new tabDownload slide Redshift distribution of the galaxy cluster catalogues used in our analysis: WHY18 (limited to richness larger than 9), SDSSDR9, MCXCsub, and HIFLUGCS. 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-fitting χ2’s of a flat angular power spectrum in multipole and the physical AGN model. Catalogue . Type . Reference . Total . SNR . |$\Delta \chi ^2_{\mathrm{FLAT-AGN}}$| . WHY18 Infrared Wen, Han & Yang (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 Catalogue . Type . Reference . Total . SNR . |$\Delta \chi ^2_{\mathrm{FLAT-AGN}}$| . WHY18 Infrared Wen, Han & Yang (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 Open in new tab 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-fitting χ2’s of a flat angular power spectrum in multipole and the physical AGN model. Catalogue . Type . Reference . Total . SNR . |$\Delta \chi ^2_{\mathrm{FLAT-AGN}}$| . WHY18 Infrared Wen, Han & Yang (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 Catalogue . Type . Reference . Total . SNR . |$\Delta \chi ^2_{\mathrm{FLAT-AGN}}$| . WHY18 Infrared Wen, Han & Yang (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 Open in new tab 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° 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 centre 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 centre): 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: $$\begin{eqnarray} C_{\ell ,{\mathrm{FLAT}}}^{c_j \gamma _i}=\left(\frac{1}{N_{c_j}}\sum _{k=0}^{N_{c_j}}\, F^\gamma _{\Delta E_i,k}\right)W_\ell ^{\Delta E_i} \, , \end{eqnarray}$$ (2) 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 centre of the cluster k in the energy bin ΔEi. 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 equation (2) can be approximated as: $$\begin{eqnarray} 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}\, , \end{eqnarray}$$ (3) where |$\bar{E}_i=\sqrt{E_i\, E_{i+1}}$| is the geometric mean of the lower (Ei) and upper (Ei + 1) bounds of the energy bin, ΔEi = Ei + 1 − Ei is the size of the bin, |$C_P^j$| is the normalization of the power-law relation, and α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 haloes. 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 modelling 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: $$\begin{eqnarray} C_\ell ^{c_j \gamma _i} = \int \frac{\mathrm{ d}\chi }{\chi ^2} \ W_{c_j}(\chi) W_{\gamma _i}(\chi) \ P_{c_j \gamma _i}(k=\ell /\chi ,\chi)\, , \end{eqnarray}$$ (4) 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 χ 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: $$\begin{eqnarray} W_{c_j}(z) = \frac{4\pi \chi (z)^2}{N_{c_j}} \ \int \ \mathrm{ d}M \ \frac{\mathrm{ d}^2n_{c_j}}{\mathrm{ d}M \ \mathrm{ d}V}\, , \end{eqnarray}$$ (5) with the number of objects in the cluster catalogue given by $$\begin{eqnarray} N_{c_j} = \int \ \mathrm{ d}M \ \mathrm{ d}V \ \frac{\mathrm{ d}^2n_{c_j}}{\mathrm{ d}M \ \mathrm{ d}V} \ . \end{eqnarray}$$ (6) We empirically derive the cluster mass function from the catalogues themselves, considering the estimated redshift and M500 masses provided by the cluster catalogues. The window function of the gamma-ray emission from a given source population is $$\begin{eqnarray} W_{\gamma _i}(E,z) &=& \chi (z)^2 \int _{\mathcal {L}_{\rm min}}^{\mathcal {L}_{\rm max}(F_{\rm sens},z)} \mathrm{ d}\mathcal {L} \, \Phi _{\rm S}(\mathcal {L},z,E) \, \frac{\mathrm{ d}N_{\rm S}}{\mathrm{ d}E}\left(\mathcal {L},z\right)\nonumber\\ &&\times \, e^{-\tau \left[E(1+z),z\right]} \ , \end{eqnarray}$$ (7) where |$\mathcal {L}$| is the gamma-ray rest-frame luminosity in the energy interval 0.1 to |$100\, \mathrm{GeV}$|⁠, ΦS is the gamma-ray luminosity function of the source class i of astrophysical emitters included in our analysis, and dNS/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 Fsens 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 modelled as in Branchini et al. (2017). Finally, the last ingredient we need is the 3D power spectrum, here decomposed in its 1-halo and 2-halo terms: $$\begin{eqnarray} P^{1h}_{c_j,\gamma _j}(k,z) = \int ^{\mathcal {L}_{\mathrm{max},i}(z)}_{\mathcal {L}_{\mathrm{min},i}(z)} \ \mathrm{ d}\mathcal {L} \ \phi _i(\mathcal {L},z) \times \frac{\mathcal {L}}{\langle f_{\gamma _i}\rangle }\frac{\langle N_{c_j}(M(\mathcal {L}))\rangle }{\bar{n}_{c_j}}\, \end{eqnarray}$$ (8) $$\begin{eqnarray} P^{2h}_{c_j,\gamma _i}(k,z) &=& \biggl [\int ^{\mathcal {L}_{\mathrm{max},i}(z)}_{\mathcal {L}_{\mathrm{min},i}(z)}\ \mathrm{ d}\mathcal {L} \ \Phi _i(\mathcal {L},z)b_{h}(\mathcal {M(L)})\frac{\mathcal {L}}{\langle f_{\gamma _i}\rangle }\biggr ] \nonumber\\ &&\times \, \biggl [\int _{M_\mathrm{min}}^{M_\mathrm{max}} \ \mathrm{ d}M \ \frac{\mathrm{ d}n}{\mathrm{ d}M}b_h(M)\frac{\langle N_{c_j}\rangle }{\bar{n}_{c_j}}\biggr ] \ P^{\mathrm{lin}}(k)\, , \end{eqnarray}$$ (9) where |$\langle f_{\gamma _i}\rangle$| is the mean luminosity density defined as |$\langle f_{\gamma _i}\rangle = \int \ \mathrm{ 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 \ \mathrm{ d}M \ \langle N_{c_j} \mathrm{ d}n / \mathrm{ d}M \rangle$|⁠. The halo mass function dn/dM and bias bh 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 normalization parameter. This second model is therefore defined as: $$\begin{eqnarray} C_{\ell ,\mathrm{AGN}}^{c_j \gamma _i} \ &=& \ C_{\ell ,{\mathrm{FLAT}}}^{c_j \gamma _i}(C_p^j,\alpha _{1,j}) \nonumber\\ &&+ \ 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}; \end{eqnarray}$$ (10) it contains the expected small-scale shot-noise term already introduced in equation (3) to which we add the AGN-like model discussed above, with the fudging free parameters α2 and A. |$C_\ell ^{\Delta E_0}$| is the theoretical APS of equation (4) calculated in a specific energy bin ΔE0, 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 α2 and the size by the free normalization 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 α1 and the normalization factor CP. It is embedded in the more complete physical AGN-like model, which is endowed with two additional free parameters, the spectral index α2 and the normalization 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 intracluster medium. The quantity we measure is the APS Cℓ 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 pixelled 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 pixellation 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): $$\begin{eqnarray} \mathrm{SNR} = \frac{C_{\ell ,\mathrm{DATA}} \ \Gamma ^{-1}_{\ell \ell ^{\prime }} \ C_{\ell ,\mathrm{MODEL}}}{\sqrt{C_{\ell ,\mathrm{MODEL}}\ \Gamma ^{-1}_{\ell \ell ^{\prime }} \ C_{\ell ,\mathrm{MODEL}}}}\, , \end{eqnarray}$$ (11) where Cℓ, 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ℓ, MODEL = Cℓ, FLAT, where Cℓ, FLAT is defined in equation (3), i.e. a model independent of the multipole ℓ (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-fitting 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 χ2 analysis, where the χ2 function is defined as: $$\begin{eqnarray} \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 ] \ ({\bf \Gamma }_{\ell \ell ^{\prime }}^{-1})^i \nonumber\\ &&\times \, \Bigl [(C_{\ell ^{\prime },\mathrm{DATA}})^i - (C_{\ell ^{\prime },\mathrm{MODEL}})^i\Bigr ] \, , \end{eqnarray}$$ (12) where index i denotes the energy bins and the sum over the multipole ℓ is limited to a range (ℓmin, ℓ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 ℓmin and ℓmax for the different energy bins are reported in Table 1. The best-fitting FLAT and AGN models are determined by maximixing the Gaussian likelihood |$\mathcal {L}$| defined as: $$\begin{eqnarray} -2 \ln {\mathcal {L}} = \chi ^2 . \end{eqnarray}$$ (13) The preference for the AGN model over the FLAT case is then determined by means of the |$\Delta \chi ^2=\chi _{\mathrm{FLAT}}^2-\chi ^2_{\mathrm{AGN}}$|⁠. This new quantity behaves approximately like a χ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: $$\begin{eqnarray} \mathrm{AIC} = 2k - 2 \ \mathrm{ln} \ \mathcal {L} \, , \end{eqnarray}$$ (14) 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 ℓ and energy, can be well approximated by a Gaussian estimate, diagonal in both dimensions, especially when the study is performed on binned (in ℓ and energy) data. While we retain in our analysis the full covariance for ℓ, 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 autocorrelations 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 randomizing 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 randomization method to produce the mock gamma ray maps and the FLASK lognormal 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 randomization 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 equation A6). 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, Abdalla & Joachimi (2016). FLASK produces lognormal realizations 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 randomized 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 lognormal probability distribution function: in our case we adopted a lognormal distribution, in order to include effects due to non-Gaussianity in the covariance matrix. To estimate the cross-correlation covariance matrix, we produce 2000 realizations of the gamma maps in each of the nine energy bin listed in Table 1 and 2000 mocks maps for each of the cluster catalogues. The large number of realizations 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 summarized here in their main features, and discussed in the Appendix. From these mock realizations, we can then construct the full covariance between different multipoles (for each energy bin, labelled by index i), obtained as: $$\begin{eqnarray} \Gamma _{\ell \ell ^{\prime }}^{c\gamma _i} &\equiv& \textrm{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} - \bar{C_\ell }^{c\gamma _i})\nonumber\\ &&\times \, (C_{\ell ^{\prime }}^{j,c\gamma _i} - \bar{C_{\ell ^{\prime }}}^{c\gamma _i})\, , \end{eqnarray}$$ (15) where Nγ(Nc) 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 realizations and: $$\begin{eqnarray} \bar{C_\ell }^{c\gamma _i} = \frac{1}{N} \sum _{a=1}^{\rm N} C_\ell ^{a,c\gamma _i} \end{eqnarray}$$ (16) is the mean of the cross APS over all the realizations. 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γ realizations of the gamma-ray mocks with the true cluster map, and |$\Gamma _{\ell \ell ^{\prime }}^{c\hat{\gamma _i}}$| obtained by correlating the Nc realizations of cluster mocks with the true gamma-ray map: $$\begin{eqnarray} \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)\, . \end{eqnarray}$$ (17) Each of the ‘half’ covariance is obtained by using equation (15). This technique allows us to reduce significantly the computing time, since we need 2N combinations instead of N2. A selection of the performed tests is shown in Figs 3 and 4, where results are reported in terms of the binned covariance (the variance and covariance are binned over intervals of size Δl = 60 with ℓmin = 40 and ℓmax = 1000 5). 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: $$\begin{eqnarray} \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,i}}} \delta _{ll^{\prime }}^{\rm K} \, , \end{eqnarray}$$ (18) where |$C_\ell ^{cc}$| and |$C_\ell ^{\gamma _i\gamma _i}$| denote the cluster and gamma-ray autocorrelations, fsky, i is the fraction of the sky probed by the survey in the energy bin i, Δℓ 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 ℓ). 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-Cl covariance matrix is described in detail. Here we just list the main steps. The code: (i) computes the Cl from the autocorrelation function; (ii) corrects the autocorrelation function for the effect of the mask; (iii) computes the Cl back; (iv) computes the V-matrix (see equation 15a of Efstathiou (2004)); (v) estimates the final covariance matrix, corrected for the mask, beam and pixel effects, as in equation (17) of Efstathiou (2004). From Fig. 3 we notice, as expected, that the Polspice variance overestimates 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{{\ \rm per\ cent}}$| at very large scales (ℓ < 100) for WHY18 and SDSSDR9, and always smaller than 10 per cent for MCXCsub and HIFLUGCS; at smaller scales (ℓ > 400) they are slightly larger than |$10{{\ \rm per\ cent}}$| for SDSSDR9, while smaller for all the other cluster catalogues; on intermediate scales (100 < ℓ < 300) we observe differences between 5 per cent and 10 per cent for all the catalogues. Figure 3. Open in new tabDownload slide Comparison of the binned cross-APS variance estimated from mocks (blue) with the variance computed with Polspice (magenta) or estimated through the Gaussian prediction (black), in the third energy bin of Table 1. The upper panel shows the results for the WHY18 and SDSSDR9 catalogues while the lower panel stands for the MCXCsub and HIGLUCS catalogues. For each panel, the ratio between the result obtained with the mock catalogues and the Gaussian prediction is shown. The shaded areas highlight the levels of 5 per cent (dark grey) and 10 per cent (light grey) deviations. Figure 3. Open in new tabDownload slide Comparison of the binned cross-APS variance estimated from mocks (blue) with the variance computed with Polspice (magenta) or estimated through the Gaussian prediction (black), in the third energy bin of Table 1. The upper panel shows the results for the WHY18 and SDSSDR9 catalogues while the lower panel stands for the MCXCsub and HIGLUCS catalogues. For each panel, the ratio between the result obtained with the mock catalogues and the Gaussian prediction is shown. The shaded areas highlight the levels of 5 per cent (dark grey) and 10 per cent (light grey) deviations. Figure 4. Open in new tabDownload slide Cross-correlation coefficient r, as defined in equation (19), that shows the size of the off-diagonal terms of the covariance matrix for the multipole dimension. Top left (right) panels show the WHY18 (SDSSDR9) case, bottom left (right) the MCXCsub (HIFLUGCS) case. Each panel refers to a specific multipole value ℓ2 and shows the coefficient r as a function of a different multipole ℓj, reported on the horizontal axis. The blue lines stand for the analysis done with the mock catalogues, while the magenta lines refer to the covariance obtained with Polspice. The relative size of the off-diagonal terms of the covariance matrix as compared to the diagonal ones is therefore always below the few per cent level. The peaks occur when the coefficient r sits on the diagonal (ℓj = ℓ2), where r = 1, by definition. The various panels show results for some representative multipoles ℓ2 = 56, 89, 141, 224, 447, 981, for each of the four cluster catalogues. Figure 4. Open in new tabDownload slide Cross-correlation coefficient r, as defined in equation (19), that shows the size of the off-diagonal terms of the covariance matrix for the multipole dimension. Top left (right) panels show the WHY18 (SDSSDR9) case, bottom left (right) the MCXCsub (HIFLUGCS) case. Each panel refers to a specific multipole value ℓ2 and shows the coefficient r as a function of a different multipole ℓj, reported on the horizontal axis. The blue lines stand for the analysis done with the mock catalogues, while the magenta lines refer to the covariance obtained with Polspice. The relative size of the off-diagonal terms of the covariance matrix as compared to the diagonal ones is therefore always below the few per cent level. The peaks occur when the coefficient r sits on the diagonal (ℓj = ℓ2), where r = 1, by definition. The various panels show results for some representative multipoles ℓ2 = 56, 89, 141, 224, 447, 981, for each of the four cluster catalogues. The ‘Gaussianity’ of the covariance matrix was not an expected result as it was not expected a large overestimation 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 Fig. 4 we show the cross-correlation coefficient defined as: $$\begin{eqnarray} r = \frac{{\bf \Gamma }_{\ell \ell ^{\prime }}^{c\gamma _i}}{\sqrt{{\bf \Gamma }_{\ell \ell }^{c\gamma _i} \ {\bf \Gamma }_{\ell ^{\prime }\ell ^{\prime }}^{c\gamma _i}}}\, . \end{eqnarray}$$ (19) One can note that the off-diagonal terms of the mocks covariance matrix are always smaller than 5 per cent 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): $$\begin{eqnarray} \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}}}\, , \end{eqnarray}$$ (20) where i and j identify the different energy bins and by determining the corresponding correlation coefficient: $$\begin{eqnarray} 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}}} , \end{eqnarray}$$ (21) which is analogous to the coefficient defined in equation (19) to investigate the off-diagonal terms of the covariance matrix for what concerns the multipole. Fig. 5 shows rE 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 per cent deviation from the diagonal elements. Only for smaller angular scales the effect reaches deviations of the order of 10 per cent. Results are similar at different angular scales and for the other cluster catalogues. Figure 5. Open in new tabDownload slide Cross-correlation coefficient rE, as defined in equation (21), that shows the size of the off-diagonal terms of the covariance matrix for the energy dimension. The four panels refer to four arbitrary energy slices of the covariance matrix, representative of the general behaviour. Each line, identified by a different colour, stands for a different multipole l, while the horizontal scale refers to the integer index labelling the energy bins. The labels Ei = 1, 3, 5, 7 on top of each panel refer again to the energy bin. The relative size of the off-diagonal terms of the covariance matrix as compared to the diagonal ones is therefore almost always below the 10 per cent level. The peaks occur when the coefficient rE sits on the diagonal (i = 1 in the first panel, and so on), where rE = 1, by definition. Figure 5. Open in new tabDownload slide Cross-correlation coefficient rE, as defined in equation (21), that shows the size of the off-diagonal terms of the covariance matrix for the energy dimension. The four panels refer to four arbitrary energy slices of the covariance matrix, representative of the general behaviour. Each line, identified by a different colour, stands for a different multipole l, while the horizontal scale refers to the integer index labelling the energy bins. The labels Ei = 1, 3, 5, 7 on top of each panel refer again to the energy bin. The relative size of the off-diagonal terms of the covariance matrix as compared to the diagonal ones is therefore almost always below the 10 per cent level. The peaks occur when the coefficient rE sits on the diagonal (i = 1 in the first panel, and so on), where rE = 1, by definition. 6 RESULTS The measured cross-correlation APS between the galaxy clusters and the unresolved gamma-ray intensity are shown in Figs 6 and 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 binned angular power spectrum |$C_\ell ^{c \gamma _i}$| as a function of the multipole ℓ (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 ℓ of the |$C_\ell ^{c \gamma _i}$| in each energy bin: $$\begin{eqnarray} P_{E_i} = \frac{1}{\Delta \ell }\sum _{\ell _{\rm min}}^{\ell _{\rm max}} C_\ell ^{c\gamma _i}\, , \end{eqnarray}$$ (22) where ℓmin and ℓmax are shown in Table 1 and Δℓ = ℓmax − ℓmin. The errors on PE are defined as: $$\begin{eqnarray} \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)\, , \end{eqnarray}$$ (23) where |$\sigma ^2_{C_\ell ^{c\gamma _i}} = \Gamma _{\ell \ell }^{c\gamma _i}$|⁠. To easy the visualization in the plot, the data of equation (22) are multiplied by E2.2 (expected behaviour of the UGRB; Ackermann et al. 2015) and divided by ΔE (the width of the energy bin). Figure 6. Open in new tabDownload slide Binned angular power spectrum of the gamma-ray–cluster cross-correlation for the four cluster catalogues adopted in our analysis and in the third gamma-ray energy bin (chosen as a representative example). The error bars are obtained from the binned angular-power-spectrum covariance matrix estimated from 2000 mocks. Figure 6. Open in new tabDownload slide Binned angular power spectrum of the gamma-ray–cluster cross-correlation for the four cluster catalogues adopted in our analysis and in the third gamma-ray energy bin (chosen as a representative example). The error bars are obtained from the binned angular-power-spectrum covariance matrix estimated from 2000 mocks. Figure 7. Open in new tabDownload slide Energy spectrum PE of the cross-correlation angular power spectrum, for each of the four cluster catalogues. The plot shows PE rescaled by E2.2/ΔE, where ΔE is the width of the corresponding energy bin. The error bars are obtained from the angular-power-spectrum covariance matrix estimated from 2000 mocks. Figure 7. Open in new tabDownload slide Energy spectrum PE of the cross-correlation angular power spectrum, for each of the four cluster catalogues. The plot shows PE rescaled by E2.2/ΔE, where ΔE is the width of the corresponding energy bin. The error bars are obtained from the angular-power-spectrum covariance matrix estimated from 2000 mocks. In order to determine the presence of a positive cross-correlation signal, we adopt the SNR defined in equation (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 an MCMC technique to determine the likelihood of equation (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 means 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σ level, whilst the other catalogues do not show a preference for the AGN model over the FLAT one. Thus MCXCsub, with an SNR of 3.5 and a Δχ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σ 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 AICFLAT = 152.05 against AICAGN = 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 Δχ2 analysis. All the other three catalogues show instead a preference for the FLAT model, namely AICFLAT < AICAGN . 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 2D allowed regions at 68 per cent (dark blue) and 95 per cent (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 ∼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 normalization, instead, turns out to be much larger than expected, even though with a sizeable uncertainty. The model adopted in our analysis is normalized 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 normalization 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 equation (10). Figure 8. Open in new tabDownload slide Triangular plot for the bounds on the AGN model parameters obtained from the fit to the cross-correlation with the MCXCsub catalogue. The contours refer to the 68 per cent (dark blue) and 95 per cent (light blue) confidence levels. Figure 8. Open in new tabDownload slide Triangular plot for the bounds on the AGN model parameters obtained from the fit to the cross-correlation with the MCXCsub catalogue. The contours refer to the 68 per cent (dark blue) and 95 per cent (light blue) confidence levels. 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 \mathrm{ d}\chi \, W_\gamma (\chi)$|⁠, and |$C_\ell ^{c\gamma ,2h}=\int \mathrm{ 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 〈bi〉 is the bias of source i with respect to matter. For a generic population emitting gamma-rays, 〈bγ〉 ∼ 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γ and in turn Iγ, 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 two different haloes 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 θ500 = r500/dA(z), where r500 is the virial radius relative to an overdensity of 500 (as defined in footnote 4) and dA(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} \gt \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 Fig. 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. Figure 9. Open in new tabDownload slide Distribution of θ500 for MCXCsub. The vertical solid line is set at θ500 = 0.2°, as an indication of the average size of the Fermi-LAT PSF (the values of the 68 per cent containment angles in the different energy bins are reported in Table 1), while the dashed vertical line shows the average angular size |$\bar{\theta } = 0.267$| of the clusters in the MCXCsub catalogue. Figure 9. Open in new tabDownload slide Distribution of θ500 for MCXCsub. The vertical solid line is set at θ500 = 0.2°, as an indication of the average size of the Fermi-LAT PSF (the values of the 68 per cent containment angles in the different energy bins are reported in Table 1), while the dashed vertical line shows the average angular size |$\bar{\theta } = 0.267$| of the clusters in the MCXCsub catalogue. 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-fitting results are shown in 10 and the results for the |$\Delta \chi ^2_{\mathrm{ AGN}-\mathrm{ FLAT}}$| turns out to be only 2.1, which corresponds to 0.96σ. 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 intracluster medium emission. The value of the best-fitting parameters are similar to what is found for the full MCXCSub case, with a rather large normalization parameter: |$A = 65.0^{+24.8}_{-35.7}$|⁠. This time, the parameter is consistent at 1.8σ with an interpretation in terms of LSS from AGN emission, although the large uncertainty does not allow to make firm conclusions. Figure 10. Open in new tabDownload slide The same as in Fig. 8 but for the MCXCsub clusters selected according the criterion |$\theta _{500}\gt \bar{\theta }_{500}$|⁠. Figure 10. Open in new tabDownload slide The same as in Fig. 8 but for the MCXCsub clusters selected according the criterion |$\theta _{500}\gt \bar{\theta }_{500}$|⁠. 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}\gt \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. Figure 11. Open in new tabDownload slide Angular two-point cross-correlation function for the most extended MCXCsub clusters (those with |$\theta _{500}\gt \bar{\theta }_{500}$|⁠) and gamma-ray energies in the range (1,10) GeV. The shaded blue area is an estimate of the error obtained from the diagonal of the Polspice covariance matrix; the grey vertical indicates the region where the Fermi PSF effects are not negligible. Figure 11. Open in new tabDownload slide Angular two-point cross-correlation function for the most extended MCXCsub clusters (those with |$\theta _{500}\gt \bar{\theta }_{500}$|⁠) and gamma-ray energies in the range (1,10) GeV. The shaded blue area is an estimate of the error obtained from the diagonal of the Polspice covariance matrix; the grey vertical indicates the region where the Fermi PSF effects are not negligible. 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°, 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 extragalactic 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 intracluster 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 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 than 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 intracluster 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°, 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 summarized in the Appendices and refer to development, test, and comparison of general numerical techniques for the massive and efficient production of mock realizations 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-QX is supported by the National Science Foundation of China under grant No.11633001 and No.11690023. Footnotes 1 We used the LAT Science Tools version https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/ 2 See http://www.slac.stanford.edu/exp/glast/groups/canda/lat_ Performance.htm for further details on photon event classes. 3 https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html 4 In our analysis we consider M500 corresponding to an overdensity Δc = 500. The mass in this case is defined as |${\rm \mathit{ M}_{500}} = 4\pi \, r_{500}^3 \, 500\, \rho _c/3 = 250\, r_{500}^3\, H(z)^2/G$|⁠, where ρc is the critical density, r500 is the virial radius, H(z) the Hubble parameter, and G the gravitational constant. 5 It 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. 6 https://healpix.jpl.nasa.gov/html/facilitiesnode14.htm 7 http://www.astro.iag.usp.br/ flask REFERENCES Ackermann M. et al. . , 2010 , ApJ , 717 , L71 10.1088/2041-8205/717/1/L71 Crossref Search ADS Crossref Ackermann M. et al. . , 2015 , ApJ , 799 , 86 10.1088/0004-637X/799/1/86 Crossref Search ADS Crossref Ackermann M. et al. . , 2018 , Phys. Rev. Lett. , 121 , 241101 10.1103/PhysRevLett.121.241101 Crossref Search ADS PubMed Crossref Ajello M. et al. . , 2017 , ApJS , 232 , 18 10.3847/1538-4365/aa8221 Crossref Search ADS Crossref Akaike H. , 1974 , in Parzen E. , Tanabe K. , Kitagawa G. , eds, Selected Papers of Hirotugu Akaike . Springer , p. 215 Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Alonso D. , Salvador A. I. , Sánchez F. J. , Bilicki M. , García-Bellido J. , Sánchez E. , 2015 , MNRAS , 449 , 670 10.1093/mnras/stv309 Crossref Search ADS Crossref Ammazzalorso S. , Fornengo N. , Horiuchi S. , Regis M. , 2018 , Phys. Rev. , D98 , 103007 10.1103/PhysRevD.98.103007 Crossref Ando S. , Ishiwata K. , 2016 , J. Cosmol. Astropart. Phys. , 2016 , 045 10.1088/1475-7516/2016/06/045 Crossref Search ADS Crossref Ando S. , Benoit-Lévy A. , Komatsu E. , 2014 , Phys. Rev. D , 90 , 023514 Banerjee P. , Szabo T. , Pierpaoli E. , Franco G. , Ortiz M. , Oramas A. , Tornello B. , 2018 , New Astron. , 58 , 61 10.1016/j.newast.2017.07.008 Crossref Search ADS Crossref Becker M. R. et al. . , 2016 , Phys. Rev. , D94 , 022002 10.1103/PhysRevD.94.022002 Crossref Bilicki M. , Jarrett T. H. , Peacock J. A. , Cluver M. E. , Steward L. , 2014 , ApJS , 210 , 9 10.1088/0067-0049/210/1/9 Crossref Search ADS Crossref Branchini E. , Camera S. , Cuoco A. , Fornengo N. , Regis M. , Viel M. , Xia J.-Q. , 2017 , ApJS , 228 , 8 10.3847/1538-4365/228/1/8 Crossref Search ADS Crossref Brunetti G. , Zimmer S. , Zandanel F. , 2017 , MNRAS , 472 , 1506 10.1093/mnras/stx2092 Crossref Search ADS Crossref Camera S. , Fornasa M. , Fornengo N. , Regis M. , 2015 , J. Cosmol. Astropart. Phys. , 2015 , 029 10.1088/1475-7516/2015/06/029 Crossref Search ADS Crossref Chon G. , Challinor A. , Prunet S. , Hivon E. , Szapudi I. , 2004 , MNRAS , 350 , 914 10.1111/j.1365-2966.2004.07737.x Crossref Search ADS Crossref Collister A. A. , Lahav O. , 2004 , PASP , 116 , 345 10.1086/383254 Crossref Search ADS Crossref Cuoco A. , Xia J.-Q. , Regis M. , Branchini E. , Fornengo N. , Viel M. , 2015 , ApJS , 221 , 29 10.1088/0067-0049/221/2/29 Crossref Search ADS Crossref De Domenico M. , Lyberis H. , 2012 , preprint (arXiv:1204.1365) Efstathiou G. , 2004 , MNRAS , 349 , 603 10.1111/j.1365-2966.2004.07530.x Crossref Search ADS Crossref Feng C. , Cooray A. , Keating B. , 2017 , ApJ , 836 , 127 10.3847/1538-4357/836/1/127 Crossref Search ADS Crossref Foreman-Mackey D. , Hogg D. W. , Lang D. , Goodman J. , 2013 , PASP , 125 , 306 10.1086/670067 Crossref Search ADS Crossref Fornengo N. , Perotto L. , Regis M. , Camera S. , 2015 , ApJ , 802 , L1 10.1088/2041-8205/802/1/L1 Crossref Search ADS Crossref Górski K. M. , Hivon E. , Banday A. J. , Wandelt B. D. , Hansen F. K. , Reinecke M. , Bartelmann M. , 2005 , ApJ , 622 , 759 10.1086/427976 Crossref Search ADS Crossref Hambly N. C. , Irwin M. J. , MacGillivray H. T. , 2001 , MNRAS , 326 , 1295 10.1111/j.1365-2966.2001.04661.x Crossref Search ADS Crossref Hashimoto D. , Nishizawa A. J. , Shirasaki M. , Macias O. , Horiuchi S. , Tashiro H. , Oguri M. , 2019 , MNRAS , 484 , 5256 10.1093/mnras/stz321 Crossref Search ADS Crossref Hu W. , Jain B. , 2004 , Phys. Rev. D , 70 , 043009 Hu W. , Kravtsov A. V. , 2003 , ApJ , 584 , 702 10.1086/345846 Crossref Search ADS Crossref Kepner J. , Fan X. , Bahcall N. , Gunn J. , Lupton R. , Xu G. , 1999 , ApJ , 517 , 78 10.1086/307160 Crossref Search ADS Crossref Lisanti M. , Mishra-Sharma S. , Rodd N. L. , Safdi B. R. , Wechsler R. H. , 2018 , Phys. Rev. D , 97 , 063005 10.1103/PhysRevD.97.063005 Crossref Search ADS Crossref Norberg P. , Baugh C. M. , Gaztañaga E. , Croton D. J. , 2009 , MNRAS , 396 , 19 10.1111/j.1365-2966.2009.14389.x Crossref Search ADS Crossref Piffaretti R. , Arnaud M. , Pratt G. W. , Pointecouteau E. , Melin J.-B. , 2011 , A&A , 534 , A109 10.1051/0004-6361/201015377 Crossref Search ADS Crossref Regis M. , Xia J.-Q. , Cuoco A. , Branchini E. , Fornengo N. , Viel M. , 2015 , Phys. Rev. Lett. , 114 , 241301 10.1103/PhysRevLett.114.241301 Crossref Search ADS PubMed Crossref Reiprich T. H. , Böhringer H. , 2002a , ApJ , 567 , 716 10.1086/338753 Crossref Search ADS Crossref Reiprich T. H. , Böhringer H. , 2002b , ApJ , 567 , 716 10.1086/338753 Crossref Search ADS Crossref Reiss I. , Keshet U. , 2018 , J. Cosmol. Astropart. Phys. , 2018 , 010 10.1088/1475-7516/2018/10/010 Crossref Search ADS Crossref Reiss I. , Mushkin J. , Keshet U. , 2018 , J. Cosmol. Astropart. Phys. , 1810 , 010 10.1088/1475-7516/2018/10/010 Crossref Search ADS Crossref Schellenberger G. , Reiprich T. H. , 2017 , MNRAS , 469 , 3738 10.1093/mnras/stx1022 Crossref Search ADS Crossref Sheth R. K. , Tormen G. , 1999 , MNRAS , 308 , 119 10.1046/j.1365-8711.1999.02692.x Crossref Search ADS Crossref Shirasaki M. , Horiuchi S. , Yoshida N. , 2014 , Phys. Rev. D , 90 , 063502 Shirasaki M. , Horiuchi S. , Yoshida N. , 2015 , Phys. Rev. D , 92 , 123540 Shirasaki M. , Macias O. , Horiuchi S. , Shirai S. , Yoshida N. , 2016 , Phys. Rev. D , 94 , 063522 The Fermi-LAT collaboration , 2019a , preprint (arXiv:1902.10045) The Fermi-LAT collaboration , 2019b , preprint (arXiv:1905.10771) van Weeren R. J. , de Gasperin F. , Akamatsu H. , Brüggen M. , Feretti L. , Kang H. , Stroe A. , Zandanel F. , 2019 , Space Sci. Rev. , 215 , 16 Voges W. et al. . , 1999 , A&A , 349 , 389 Wen Z. L. , Han J. L. , Yang F. , 2018 , MNRAS , 475 , 343 10.1093/mnras/stx3189 Crossref Search ADS Crossref Wright E. L. et al. . , 2010 , AJ , 140 , 1868 10.1088/0004-6256/140/6/1868 Crossref Search ADS Crossref Xavier H. S. , Abdalla F. B. , Joachimi B. , 2016 , MNRAS , 459 , 3693 10.1093/mnras/stw874 Crossref Search ADS Crossref Xia J.-Q. , Cuoco A. , Branchini E. , Viel M. , 2015 , ApJS , 217 , 15 10.1088/0067-0049/217/1/15 Crossref Search ADS Crossref York D. G. et al. . , 2000 , AJ , 120 , 1579 10.1086/301513 Crossref Search ADS Crossref Zandanel F. , Ando S. , 2014 , MNRAS , 440 , 663 10.1093/mnras/stu324 Crossref Search ADS Crossref 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 realizations 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 × 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 × N crossing. In the reminder of this Appendix we discuss the various techniques adopted to generate mock maps for both galaxy/cluster catalogues and for gamma-ray maps. In Appendix B we derive a method which allows us 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 catalogue 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 reorganization 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 pixellation format with resolution parameter Nside = 1024, which correspond to a total number of pixels |$N_{\rm pix}\, = 12582\,912$| and mean spacing of ∼0.06°. 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 catalogues. 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 realizations of both the galaxies/clusters and gamma-ray maps, we test that the ensuing autocorrelation 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. A1 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 realization 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 realizations. The new HEALPix map originated in this way is characterized 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 procedure is called bootstrap with replacements (Norberg et al. 2009). In this case the estimator of the APS covariance matrix is given by: $$\begin{eqnarray} \widehat{\Gamma _{\rm B}} = \frac{1}{N_r-1} \ \sum _{k=1}^{N_k} \ (x^k_i - \bar{x}_i)(x^k_j - \bar{x}_j)\, , \end{eqnarray}$$ (A1) where xi is the i-th bootstrap realization and |$N_{\rm r}\,$| is the number of realizations and, $$\begin{eqnarray} \bar{x}_i = \frac{1}{N_r} \ \sum _{k=1}^{N_r} x_i^k \ . \end{eqnarray}$$ (A2) A2 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 realization is obtained by systematically omitting one of the sub-array in each realization. The resampling of the data set consists of |$N_{\rm sub}\, -1$| remaining 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: $$\begin{eqnarray} \widehat{\Gamma _{\rm J}} = \frac{N_{\rm r}\, -1}{N_{\rm r}\, } \ \sum _{k=1}^{N_{\rm r}\, } \ (x^k_i - \bar{x}_i)(x^k_j - \bar{x}_j)\, , \end{eqnarray}$$ (A3) where |$\bar{x}_i$| is given by equation (A2). The factor |$(N_{\rm r}\, -1)$| accounts for the lack of independence between the |$N_{\rm r}\,$| copies of the data set. A3 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: $$\begin{eqnarray} f(\theta ,\phi) = \sum _{\ell =0}^\infty \sum _{m=-\ell }^\ell \ a_{\ell m}\ Y_{\ell m}(\theta ,\phi)\, . \end{eqnarray}$$ (A4) from which the angular power spectrum is obtained: $$\begin{eqnarray} C_\ell = \frac{1}{2\ell +1} \ \sum _{m=-1}^\ell \ |a_{\ell m}|^2 \ . \end{eqnarray}$$ (A5) It is clear that equation (A5) is invariant under a phase rotation on the harmonic amplitudes alm: $$\begin{eqnarray} a_{\ell m}\quad \longrightarrow \quad a_{\ell m}e^{i\varphi _{\ell m}} \qquad \varphi _{\ell m} \in \mathbb {R} \, . \end{eqnarray}$$ (A6) Taking advantage of this symmetry, we can build independent realizations 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 (fsky)−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(θ, ϕ) 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 equation (A3) A4 Gaussian realizations (Synfast) Synfast6 is a HEALPix routine that allows us to generate realizations 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 equation (A3) A5 Lognormal realizations (FLASK) Flask7 is a C+ + code, parallelized with OpenMP, based on the work of Xavier et al. (2016) and created to generate mock realizations 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 3D 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 equation (A3). 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. A6 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 |${\bf X} = \lbrace x_1,x_2,...,x_{N_{\rm r}\, }\rbrace$| we can write the generic expression for the covariance matrix of X as: $$\begin{eqnarray} \textrm{cov}[{\bf X}] = \frac{1}{N_{\rm r}\, } \ \sum _{k=1}^{N_{\rm r}\, } \ (x_i^k-\bar{x}_i)(x_j^k-\bar{x}_j)\, \end{eqnarray}$$ (A7) with |$\bar{x}_i$| is given by equation (A2). Equation (A7) can also be rewritten as: $$\begin{eqnarray} \textrm{cov}[{\bf X}] = \frac{1}{N_{\rm r}\, } \ \sum _{k=1}^{N_{\rm r}\, } \ x_i^kx_j^k - \frac{1}{N_{\rm r}\, ^2} \ \sum _{k=1}^{N_{\rm r}\, }x_i^k\sum _{k=1}^{N_{\rm r}\, }x_j^k \ . \end{eqnarray}$$ (A8) The unbiased definition of the sample covariance when the mean is derived from the sample itself is: $$\begin{eqnarray} \widehat{\Gamma }_{\rm U} &\equiv & \frac{N_{\rm r}\, }{N_{\rm r}\, -1} \ \times \ \textrm{cov}[{\bf X}] \nonumber \\ &=&\frac{1}{N_{\rm r}\, -1} \ \sum _{k=1}^{N_{\rm r}\, } \ x_i^kx_j^k - \frac{1}{N_{\rm r}\, }\frac{1}{N_{\rm r}\, -1} \ \sum _{k=1}^{N_{\rm r}\, }x_i^k\sum _{k=1}^{N_{\rm r}\, }x_j^k\, . \end{eqnarray}$$ (A9) We can relate the estimator of equation (A9) to the ones obtained with the jackknife (A3) and bootstrap (A1) techniques as: $$\begin{eqnarray} \textrm{Jackknife:} \quad \widehat{\Gamma }_{\rm J} &=& \frac{(N_r-1)^2}{N_r} \ \times \ \widehat{\Gamma }_{\rm U} \end{eqnarray}$$ (A10) $$\begin{eqnarray} \textrm{Bootstrap:} \quad \widehat{\Gamma }_{\rm J} &=& \widehat{\Gamma _{\rm U}} \end{eqnarray}$$ (A11) APPENDIX B: SEMI-ANALYTIC PREDICTION OF THE CROSS-CORRELATION COVARIANCE The derivation of the covariance matrix for the cross-correlation APS requires us to combine the information arising from the generation of many different set of maps. In order to obtain stable results, the required number of realizations 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 us 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 nE energy bins and nz 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 realizations, 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 autocorrelations for galaxies and gamma-rays. Let us start with the Gaussian prediction for the APS covariance matrix: (Hu & Jain 2004): $$\begin{eqnarray} \Gamma _{\ell \ell }^{g\gamma _i} \equiv \textrm{cov}[C_\ell ^{g\gamma _i},C_\ell ^{g\gamma _i}]= \frac{C_\ell ^{gg}C_\ell ^{\gamma _i\gamma _i} + (C_\ell ^{g\gamma _i})^2}{(2\ell +1)\Delta \ell \ f_{\rm sky}} \delta _{\ell \ell ^{\prime }}^{\rm K} \, , \end{eqnarray}$$ (B1) where Δl is the bin width, fsky is the fraction of the sky probed by the surveys, and δ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 autocorrelation 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 realizations. From equation (B1) and considering that we construct the covariance from the mock by averaging over the |$N_{\rm r}\,$| realizations: $$\begin{eqnarray} \Gamma _{\ell \ell ^{\prime }}^{\hat{g}\gamma _i} &\equiv& \textrm{cov}[C_\ell ^{\hat{g}\gamma _i},C_{\ell ^{\prime }}^{\hat{g}\gamma _i}] \ \propto \ \hat{C}_\ell ^{gg}\frac{1}{N_r}\sum _{n=1}^{N_r} \ C_\ell ^{\gamma _i\gamma _i,n} \nonumber\\ &&+\frac{1}{N_r}\sum _{n=1}^{N_r} \ (C_\ell ^{\hat{g}\gamma _i,n})^2 . \end{eqnarray}$$ (B2) Let us then define the corresponding counterpart term obtained by using the mocks for galaxies/clusters and data for gamma-rays: $$\begin{eqnarray} \Gamma _{\ell \ell ^{\prime }}^{g\hat{\gamma _i}} \propto \ \hat{C}_l^{\gamma _i\gamma _i}\frac{1}{N_r}\sum _{n=1}^{N_r} \ C_l^{gg,n} + \frac{1}{N_r}\sum _{n=1}^{N_r} \ (C_l^{g\hat{\gamma _i},n})^2 \, . \end{eqnarray}$$ (B3) When instead we use only mocks, equation (B1) gives: $$\begin{eqnarray} \Gamma _{\ell \ell ^{\prime }}^{g\gamma _i} \ &\propto& \ \left(\frac{1}{N_r}\sum _{n=1}^{N_r} \ C_l^{gg,n} \right) \left(\frac{1}{N_r}\sum _{n=1}^{N_r} \ C_l^{\gamma _i\gamma _i,n} \right) \nonumber\\ &&+ \frac{1}{N_r^2}\sum _{n=1}^{N^2_r} \ (C_l^{g\gamma _i,n})^2 . \end{eqnarray}$$ (B4) If we take the average of the expressions (B2) and (B3), we obtain: $$\begin{eqnarray} (\Gamma _{\ell \ell ^{\prime }}^{g\gamma _i})_{\rm ave} &=& \frac{1}{2} \Bigl (\Gamma _{\ell \ell ^{\prime }}^{\hat{g}\gamma _i}+\Gamma _{\ell \ell ^{\prime }}^{g\hat{\gamma _i}}\Bigr) \ \propto \ \hat{C}_l^{\gamma _i\gamma _i}\frac{1}{N_r}\sum _{n=1}^{N_r} \ C_l^{gg,n}\nonumber\\ && +\hat{C}_l^{gg}\frac{1}{N_r}\sum _{n=1}^{N_r} \ C_l^{\gamma _i\gamma _i,n} + \frac{1}{N_r}\sum _{n=1}^{N_r} \ (C_l^{g\gamma _i,n})^2 . \end{eqnarray}$$ (B5) 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: $$\begin{eqnarray} \hat{C}_l^{gg} & \simeq & \frac{1}{N_r}\sum _{n=1}^{N_r} \ C_l^{gg,n} \nonumber \\ \hat{C}_l^{\gamma _i\gamma _i} & \simeq & \frac{1}{N_r}\sum _{n=1}^{N_r} \ C_l^{\gamma _i\gamma _i,n} \, . \end{eqnarray}$$ (B6) we obtain: $$\begin{eqnarray} (\Gamma _{\ell \ell ^{\prime }}^{g\gamma _i})_{\rm ave} \ &\propto& \ \left(\frac{1}{N_r}\sum _{n=1}^{N_r} \ C_l^{gg,n} \right) \left(\frac{1}{N_r}\sum _{n=1}^{N_r} \ C_l^{\gamma _i\gamma _i,n} \right) \nonumber\\ &&+ \frac{1}{N^2_r}\sum _{n=1}^{N^2_r} \ (C_l^{g\gamma _i,n})^2 \, . \end{eqnarray}$$ (B7) We then observe that in the limit of large |$N_{\rm r}\,$|⁠: $$\begin{eqnarray} \frac{1}{N^2_r}\sum _{n=1}^{N^2_r} \ (C_l^{g\gamma _i,n})^2 \simeq \frac{1}{N_r}\sum _{n=1}^{N_r} \ (C_l^{g\gamma _i,n})^2\, . \end{eqnarray}$$ (B8) In this case, equation (B4) and equation (B7) give the same result. Therefore we can obtain a reliable estimate of the covariance by simply averaging over the sum of the two ‘half’ contributions: $$\begin{eqnarray} \Gamma _{\ell \ell ^{\prime }}^{g\gamma _i} = \frac{1}{2} \Bigl (\Gamma _{\ell \ell ^{\prime }}^{\hat{g}\gamma _i}+\Gamma _{\ell \ell ^{\prime }}^{g\hat{\gamma _i}}\Bigr)\, . \end{eqnarray}$$ (B9) 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 equation (B9) to join the separate information coming respectively from the gamma and galaxy mock maps. For definiteness, in this analysis we use 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 σz = 0.015 for nearly all the data set with few outliers. The resulting 2MPZ sample contains almost 1 million galaxies with a median redshift of z = 0.07. In the left-hand panel of Fig. C1 we show the 2MPZ catalogue in HEALPix projection and in the right-hand 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. Figure C1. Open in new tabDownload slide Left-hand panel: All-sky map of the 2MPZ catalogue in HEALPix projection with Nside = 128. Right-hand panel: Redshift distribution of the 2MPZ catalogue. Figure C1. Open in new tabDownload slide Left-hand panel: All-sky map of the 2MPZ catalogue in HEALPix projection with Nside = 128. Right-hand panel: Redshift distribution of the 2MPZ catalogue. 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 equation (B1). 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 Figs C2–C6. In each figure we show the variance estimated using equation (B7) (coloured lines) for a low (left-hand panel), intermediate (centre), and large (right-hand panel) energy bin; in all plots, the black line represents the Gaussian prediction. Each plot is accompanied by a lower panel, where we show the ratio between the result for each combination and the Gaussian variance. Figure C2. Open in new tabDownload slide Diagonal of the cross-correlation covariance matrix estimated using 256 mocks; in black the Gaussian predictions; the coloured lines are the variances estimated using the phase randomization method for gamma maps + all the other methods for galaxy maps. In the lower panels it is shown that the ratio between each of the coloured lines and the Gaussian variance; the grey shaded area represents a 50 per cent interval. Figure C2. Open in new tabDownload slide Diagonal of the cross-correlation covariance matrix estimated using 256 mocks; in black the Gaussian predictions; the coloured lines are the variances estimated using the phase randomization method for gamma maps + all the other methods for galaxy maps. In the lower panels it is shown that the ratio between each of the coloured lines and the Gaussian variance; the grey shaded area represents a 50 per cent interval. In Fig. C2 the method to produce gamma-ray mocks is fixed to phase randomization, 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{{\ \rm per\ cent}}$| than the Gaussian prediction. This happens for all energy bins. The combination phaseGAM + phaseGAL slightly underestimate by a few per cent 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 Gaussian prediction of about |$\sim 30{{\ \rm per\ cent}}$| for phaseGAM + phaseGAL/flaskGAL and about |$60{{\ \rm per\ cent}}$| for phaseGAM + jackknifeGAL. A similar behaviour is observed in Fig. C3, where synfast is used to produce the gamma-ray mocks. All the combinations show almost the same trend observed in Fig. C2. We expected similar results between the phase randomization and synfast technique, since the two technique are quite similarly implemented in the generation of mock maps. Figure C3. Open in new tabDownload slide Same as Fig. C2 but using the Synfast methods for the gamma maps. Figure C3. Open in new tabDownload slide Same as Fig. C2 but using the Synfast methods for the gamma maps. Fig. C4 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{{\ \rm per\ cent}}$| in the case of bootstrapGAM + synfastGAL/flaskGAL, and by more than |$50{{\ \rm per\ cent}}$| for bootstrapGAM + phaseGAL/bootstrapGAL. Figure C4. Open in new tabDownload slide Same as Fig. C2 but using the Bootstrap methods for the gamma maps. Figure C4. Open in new tabDownload slide Same as Fig. C2 but using the Bootstrap methods for the gamma maps. The jackknife applied to gamma-rays in combination with all methods for galaxies is shown in Fig. C5. In this case, all the combinations largely overestimate the Gaussian prediction with just jackknifeGAM + bootstrapGAL/phaseGAL remaining below a |$50{{\ \rm per\ cent}}$| difference. Figure C5. Open in new tabDownload slide Same as Fig. C2 but using the Jackknife methods for the gamma maps. Figure C5. Open in new tabDownload slide Same as Fig. C2 but using the Jackknife methods for the gamma maps. Finally, we show in Fig. C6, 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. Figure C6. Open in new tabDownload slide Same as Fig. C2 but using the Flask methods for the gamma maps. Figure C6. Open in new tabDownload slide Same as Fig. C2 but using the Flask methods for the gamma maps. 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. © 2019 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Searching for gamma-ray emission from galaxy clusters at low redshift
JF - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/stz3263
DA - 2020-01-21
UR - https://www.deepdyve.com/lp/oxford-university-press/searching-for-gamma-ray-emission-from-galaxy-clusters-at-low-redshift-QXLFiyDx7b
SP - 3225
VL - 491
IS - 3
DP - DeepDyve
ER -