Detection of cross-correlation between gravitational lensing and gamma rays

TL;DR

This paper reports the first detection of a cross-correlation between gamma-ray emissions and gravitational lensing data, revealing insights into faint cosmic sources and dark matter.

Contribution

It introduces the first observational evidence of gamma-ray and gravitational lensing cross-correlation using DES Y1 and Fermi data, highlighting potential astrophysical and dark matter sources.

Findings

Detected cross-correlation signal with SNR of 5.3

Signal localized at small angular scales and high energies

Hints of correlation at larger separations

Abstract

In recent years, many gamma-ray sources have been identified, yet the unresolved component hosts valuable information on the faintest emission. In order to extract it, a cross-correlation with gravitational tracers of matter in the Universe has been shown to be a promising tool. We report here the first identification of a cross-correlation signal between gamma rays and the distribution of mass in the Universe probed by weak gravitational lensing. We use the Dark Energy Survey Y1 weak lensing catalogue and the Fermi Large Area Telescope 9-year gamma-ray data, obtaining a signal-to-noise ratio of 5.3. The signal is mostly localised at small angular scales and high gamma-ray energies, with a hint of correlation at extended separation. Blazar emission is likely the origin of the small-scale effect. We investigate implications of the large-scale component in terms of astrophysical sources and particle dark matter emission.

Tables

Table 1. Table 1: Energy bins over which the analysis is performed, 68% and 95% containment angles θ cont subscript 𝜃 cont \theta_{\rm cont} of the Fermi -LAT PSF, and photon counts in the unmasked Fermi area in each energy bin

	Bin number
	1	2	3	4	5	6	7	8	9
$E_{\min }$ [GeV]	$0.631$	$1.202$	$2.290$	$4.786$	$9.120$	$17.38$	$36.31$	$69.18$	$131.8$
$E_{\max }$ [GeV]	$1.202$	$2.290$	$4.786$	$9.120$	$17.38$	$36.31$	$69.18$	$131.8$	$1000.0$
${\theta }_{\mathrm{cont}}$ 68% [deg]	$0.50$	$0.58$	$0.36$	$0.22$	$0.15$	$0.12$	$0.11$	$0.10$	$0.10$
${\theta }_{\mathrm{cont}}$ 95% [deg]	$1.03$	$1.06$	$0.62$	$0.39$	$0.28$	$0.22$	$0.20$	$0.18$	$0.16$
Photon counts	$345230$	$444559$	$286209$	$102821$	$41148$	$16932$	$5250$	$1728$	$722$

Table 2. Table 2: Δ ​ χ mod 2 Δ subscript superscript 𝜒 2 mod \Delta\chi^{2}_{\rm mod} and SNR mod subscript SNR mod {\rm SNR}_{\rm mod} computed for the phenomenological  and physical  models, using either the full data set or the various subsamples discussed in the text. For dark matter in the physical  model, we consider the annihilation channel τ + ​ τ − superscript 𝜏 superscript 𝜏 \tau^{+}\tau^{-} . For the Low- z 𝑧 z case we selected the two first redshift bins, while for the High- z 𝑧 z case the last two bins, where the bins are defined as: 0.20 < ⟨ z ⟩ < 0.43 0.20 delimited-⟨⟩ 𝑧 0.43 0.20<\langle z\rangle<0.43 , 0.43 < ⟨ z ⟩ < 0.63 0.43 delimited-⟨⟩ 𝑧 0.63 0.43<\langle z\rangle<0.63 , 0.63 < ⟨ z ⟩ < 0.90 0.63 delimited-⟨⟩ 𝑧 0.90 0.63<\langle z\rangle<0.90 and 0.90 < ⟨ z ⟩ < 1.30 0.90 delimited-⟨⟩ 𝑧 1.30 0.90<\langle z\rangle<1.30 ; the Low- E 𝐸 E subset is defined for energies below 5 GeV, while the High- E 𝐸 E for energies above this value; finally, the Small- θ 𝜃 \theta /Large- θ 𝜃 \theta cases correspond to data points below/above 3 times the Fermi -LAT PSF.

	Data set
	Full	Low-$z$	High-$z$	Low-$E$	High-$E$	Small-$\theta$	Large-$\theta$
$\Delta {\chi }_{\mathrm{phe}}^2$	$27.5$	$17.5$	$10.4$	$1.7$	$21.7$	$17.0$	$6.0$
${\mathrm{SNR}}_{\mathrm{phe}}$	$5.3$	$4.2$	$3.2$	$1.3$	$5.1$	$4.2$	$2.5$
$\Delta {\chi }_{\mathrm{phys}}^2$	$27.0$	$4.8$	$12.5$	$1.5$	$16.2$	$12.5$	$4.8$
${\mathrm{SNR}}_{\mathrm{phys}}$	$5.2$	$2.2$	$3.5$	$1.2$	$4.2$	$3.5$	$2.1$

Table 3. Table 3: Δ ​ χ mod 2 Δ subscript superscript 𝜒 2 mod \Delta\chi^{2}_{\rm mod} computed for the physical  and phenomenological  models with respect to the null hypothesis for the various combinations of the blind analysis. Only the combination of true γ 𝛾 \gamma -ray and shear data (CX) and the low-noise mock γ 𝛾 \gamma -ray map with an injected artificial shear signal (AZ) show a large preference for the model vs the null hypothesis (no signal). The physical model adopted in the blind analysis is the reference model discussed in the text. For CX, CY and CZ, we show in parenthesis also the values considered in the blinding phase, before the post-unblinding correction mentioned in the text.

	Map combination
	AX	AY	AZ	BX	BY	BZ	CX	CY	CZ
$\Delta {\chi }_{\mathrm{phe}}^2$	$1.17$	$0.27$	$29.3$	$4.09$	$0.26$	$0.60$	$27.5\left(21.1\right)$	$1.56\left(1.27\right)$	$1.50\left(1.90\right)$
$\Delta {\chi }_{\mathrm{phys}}^2$	$-0.06$	$0.04$	$2.58$	$3.09$	$-0.02$	$3.33$	$18.2\left(9.91\right)$	$0.65\left(0.15\right)$	$2.06\left(1.72\right)$

Equations

$${\Xi }^{ar}(\theta )={\Xi }_{\Delta {\theta }_h,\Delta E_a,\Delta z_r}^{\mathrm{signal}}-{\Xi }_{\Delta {\theta }_h,\Delta E_a,\Delta z_r}^{\mathrm{random}}=\frac{{\sum }_{i,j}{e}_{ij,\mathrm{t}}^r{I}_j^a}{R{\sum }_{i,j}{I}_j^a}-\frac{{\sum }_{i,j}{e}_{ij,\mathrm{t}}^r{I}_{j,\mathrm{random}}^a}{R{\sum }_{i,j}{I}_{j,\mathrm{random}}^a},$$

$${\Xi }_{\mathrm{phe}}^{ar}(\theta )⟨I_a⟩=A_1\times E_a^{-{\alpha }_1}\times (1+z_r{)}^{{\beta }_1}\times {\hat{\Xi }}_{\text{PSF-like}}^a(\theta )+A_2\times E_a^{-{\alpha }_2+2.2}\times (1+z_r{)}^{{\beta }_2}\times {\hat{\Xi }}_{\text{2h-like}}^{ar}(\theta )$$

$${\chi }^2({\mathbf{P}}_{\mathrm{mod}})=[{\mathbf{\Xi }}_{\mathrm{data}}-{\mathbf{\Xi }}_{\mathrm{th}}({\mathbf{P}}_{\mathrm{mod}}){]}^{\mathsf{T}}{\mathbf{\Gamma }}^{-1}[{\mathbf{\Xi }}_{\mathrm{data}}-{\mathbf{\Xi }}_{\mathrm{th}}({\mathbf{P}}_{\mathrm{mod}})]$$

$$\mathrm{SNR}({\mathbf{P}}_{\mathrm{mod}})=\frac{{\mathbf{\Xi }}_{\mathrm{data}}^{\mathsf{T}}{\mathbf{\Gamma }}^{-1}{\mathbf{\Xi }}_{\mathrm{th}}({\mathbf{P}}_{\mathrm{mod}})}{\sqrt{{\mathbf{\Xi }}_{\mathrm{th}}^{\mathsf{T}}({\mathbf{P}}_{\mathrm{mod}}){\mathbf{\Gamma }}^{-1}{\mathbf{\Xi }}_{\mathrm{th}}({\mathbf{P}}_{\mathrm{mod}})}};$$

$${\Xi }_{\mathrm{phys}}^{ar}(\theta )⟨I_a⟩=A_{\mathrm{BLZ}}^{1h}\times {\hat{\Xi }}_{\mathrm{BLZ},1h}^{ar}(\theta )+A_{\mathrm{BLZ}}^{2h}\times {\hat{\Xi }}_{\mathrm{BLZ},2h}^{ar}(\theta )+A_{\mathrm{mAGN}}\times {\hat{\Xi }}_{\mathrm{mAGN}}^{ar}(\theta )+A_{\mathrm{SFG}}\times {\hat{\Xi }}_{\mathrm{SFG}}^{ar}(\theta )+A_{\mathrm{DM}}\times {\hat{\Xi }}_{\mathrm{DM}}^{ar}(\theta ;m_{\mathrm{DM}}).$$

$$\mathbf{e}=\mathbf{e}{|}_{\gamma =0}+\frac{\partial \mathbf{e}}{\partial \mathbf{\gamma }}{\mathbf{R}}_{\gamma =0}\mathbf{\gamma }+\dots =\mathbf{e}{|}_{\gamma =0}+{\mathbf{R}}_{\gamma }\mathbf{\gamma }+\dots .$$

$${\mathbf{R}}_{{\gamma }_{i,j}}=\frac{e_i^{+}-e_i^{-}}{\Delta {\gamma }_j},$$

$$⟨{\mathbf{R}}_{S_{i,j}}⟩=\frac{⟨e_i^{S+}⟩-⟨e_i^{S-}⟩}{\Delta {\gamma }_j},$$

$$⟨\mathbf{R}⟩=⟨{\mathbf{R}}_{\gamma }⟩+⟨{\mathbf{R}}_S⟩.$$

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

$$\widehat{\Gamma}_{ar\ell,bs\ell^{\prime}}=\frac{\delta^{\rm K}_{\ell\ell^{\prime}}}{(2\ell+1)\Delta\ell f_{\rm sky}}\left[C_{\ell}^{ar}C_{\ell^{\prime}}^{bs}+\big{(}C_{\ell^{\prime}}^{rs}+\mathcal{N}^{rs}\big{)}\big{(}C_{\ell}^{ab}+\mathcal{N}^{ab}\big{)}\right].$$

$$C_{\ell }^{ar}=\int \mathrm{d}E\mathrm{d}z\frac{1}{H(z)}\frac{W_{\text{gamma}}^a(E,z)W_{\mathrm{shear}}^r(z)}{\chi (z{)}^2}{P}_{\gamma \delta }\left[k=\frac{\ell }{\chi (z)},z\right],$$

$${\hat{\Xi }}^{ar}(\theta )=\sum _{\ell }\frac{2\ell +1}{4\pi \ell (\ell +1)}{C}_{\ell }^{ar}{W}_{\ell }^a{P}_{\ell }^{(2)}(\cos \theta ),$$

$$W_{\mathrm{shear}}^r(z)=\frac{3}{2}{H}_0^2{\Omega }_{\mathrm{m}}(1+z)\chi (z){\int }_z^{\infty }\mathrm{d}{z}^{\prime }\frac{\chi (z^{\prime })-\chi (z)}{\chi (z^{\prime })}{n}^r(z^{\prime }),$$

$$W_{\text{gamma,DM}}^a(E,z)=\frac{({\Omega }_{\mathrm{DM}}{\rho }_{\mathrm{c}}{)}^2}{4\pi }\frac{⟨{\sigma }_{\mathrm{ann}}v⟩}{2m_{\mathrm{DM}}^2}{\left(1+z\right)}^3{\Delta }^2(z)\frac{\mathrm{d}{N}_{\mathrm{ann}}}{\mathrm{d}E}\left[E(1+z)\right]e^{-\tau \left[E(1+z),z\right]},$$

$${\Delta }^2(z)\equiv \frac{⟨{\rho }_{\mathrm{DM}}^2⟩}{{\bar{\rho }}_{\mathrm{DM}}^2}={\int }_{M_{\min }}^{M_{\max }}\mathrm{d}M\frac{\mathrm{d}{n}_{\mathrm{h}}}{\mathrm{d}M}(M,z)\left[1+b_{\mathrm{sub}}(M,z)\right]\times \int {\mathrm{d}}^{3}x\frac{{\rho }_{\mathrm{h}}^2(\mathbf{x}|M,z)}{{\bar{\rho }}_{\mathrm{DM}}^2},$$

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

$$P_{{\gamma }_{\mathrm{DM}}\delta }^{1h}(k,z)$$

$$P_{{\gamma }_{\mathrm{DM}}\delta }^{2h}(k,z)$$

$$P_{{\gamma }_{\mathrm{S}}\delta }^{1h}(k,z)$$

$$P_{{\gamma }_{\mathrm{S}}\delta }^{2h}(k,z)$$

Full text

††preprint: FERMILAB-PUB-19-615-AE††preprint: DES-2019-0431

I Introduction

Astronomy at $\gamma$-ray frequencies represents a promising avenue for both astrophysics and particle physics. On one hand, the most violent phenomena in the Universe produce high-energy photons that travel all the way to Earth. Thus, they bring us information about the physics of rare events such as supernovae, and the behaviour of matter under extreme conditions, as in pulsars and active galactic nuclei (AGNs). On the other hand, the most elusive form of matter in the cosmos—dark matter, which represents about 25% of all the Universe’s energy content—is believed to consist of an exotic fundamental particle, which may annihilate or decay into standard-model particles and thus produce cosmic messengers including $\gamma$-ray photons. In the weakly interacting massive particle (WIMP) scenario, or for any hypothetical dark matter particle with mass in the GeV range or higher, dark matter particle annihilation/decay almost necessarily result in photons at $\gamma$-ray energies. Therefore, $\gamma$-ray astronomy represents a promising means to investigate the fundamental nature of dark matter. However, the faintness of the expected emission makes it very difficult to identify such a signal.

Since 2008, the Large Area Telescope (LAT) mounted on the Fermi satellite have been performing the most detailed observations of the extra-galactic $\gamma$-ray sky and resolved 5065 $\gamma$-ray sources in the energy range $50\,{\rm MeV}$ to $1\,{\rm TeV}$ 2019arXiv190210045T . Once the point sources and Galactic emission are removed, the remaining $\gamma$-ray photons form the so-called unresolved $\gamma$-ray background (UGRB).

A method to discriminate between non-thermal $\gamma$-ray emission due to astrophysical sources and possible dark matter annihilation/decay in the UGRB has been proposed in Ref. Camera:2012cj . This method relies on cross-correlations of UGRB maps with maps of other tracers of the underlying structure on cosmological scales, such as the weak gravitational lensing effect, or the clustering of galaxies and galaxy clusters (see also refs. Ando2014 ; Fornengo:2014 ) and CMB lensing Fornengo:2014cya ; Feng:2016fkl . These are direct gravitational probes of matter, most of which is thought to be dark matter. The energy, redshift and scale dependence of the aforementioned cross-correlations have the potential to disentangle signatures due to astrophysics from dark matter (see also Ref. Camera:2014rja ). More generally, the method can provide valuable information on the redshift distribution and on the clustering properties of the unresolved $\gamma$-ray source populations, including blazars, AGNs and star-forming galaxies.

Since cross-correlations of the UGRB with gravitational lensing have been proposed as a probe, several observational attempts have followed Shirasaki2014 ; Shirasaki:2016kol ; Troster:2016sgf ; Shirasaki:2018dkz , but none so far has detected the signal. Here, we report the first detection of such a cross-correlation. We used 108-month $\gamma$-ray data from Fermi-LAT and first year (Y1) shear measurements from the Dark Energy Survey (DES). In the following, we describe details of the analysis and discuss the results.

II Analysis and results

The observable we probe is the cross-correlation between the unresolved component of the $\gamma$-ray emission and gravitational shear. To this aim, the Fermi-LAT data have been pre-processed to produce the relevant energy-dependent response functions of the detector and full-sky maps of photon intensities in several energy bins. Resolved $\gamma$-ray sources and the bright Galactic plane emission have been masked with energy- and flux-dependent masks, in order to minimise the sky fraction removal. Furthermore, we have subtracted a model of the Galactic plane emission. Galactic foreground emission does not lead to false detection of a cross-correlation, since it does not correlate with the large-scale structure measured by gravitational shear, but it increases the variance of the measurements (see Supplemental Material and, e.g., Refs. Xia:2015wka ; Shirasaki2014 ; Shirasaki:2016kol ; Ammazzalorso:2018evf ; Shirasaki:2018dkz ). The weak lensing information is extracted by measuring the mean tangential ellipticity of source galaxies in the DES footprint around pixels weighted by their UGRB flux. The shear catalogue is divided in redshift bins in order to perform a tomographic analysis. As an illustration of the overlapping area bewteen DES and Fermi-LAT, Fig. 1 shows the DES footprint and the Fermi-LAT map for photon energies in the 1-10 GeV interval.

We measure the cross-correlation between the UGRB and gravitational shear through its 2-point angular correlation function. Specifically, we adopt the following estimator (see also Ref. Gruen:2017xjj ):

[TABLE]

where $\Xi^{\mathrm{signal}}_{\Delta\theta_{h},\Delta E_{a},\Delta z_{r}}$ is the correlation function in configuration space of the two observables measured in different angular ($\Delta\theta_{h}$), $\gamma$-ray energy ($\Delta E_{a}$) and lensing source-galaxy redshift ($\Delta z_{r}$) bins. The correlation is obtained by summing the products of tangential ellipticity of source galaxies $i$ relative to a pixel $j$, $e^{r}_{ij,\mathrm{t}}$, multiplied by the Fermi-LAT photon intensity flux in the $a$-th energy bin and in pixel $j$, $I^{a}_{j}$. The sum runs over all unmasked pixels $j$ and all sources $i$ in the redshift bin of the shear catalogue, and it is performed in each of the different photon energy bins (labelled by $a$) and source galaxies redshift bins (labelled by $r$). Lastly, $R$ is the mean response of ellipticity to shear for sources in the redshift bin, determined by the metacalibration algorithm Sheldon:2017szh ; Zuntz:2017pso to be between 0.54 and 0.73 for the source-galaxy redshift bins used here.

From the correlation function, we remove $\Xi^{\mathrm{random}}_{\Delta\theta_{h},\Delta E_{a},\Delta z_{r}}$, the measurement of tangential shear around random lines of sight. This is done by setting $I^{a}_{j,\rm random}=1$ anywhere within the sky region used for $\gamma$-ray measurements in that energy bin and 0 elsewhere. This reduces additive shear systematic effects, random very-large-scale structures, or chance shear alignments relative to the mask. The random subtraction, while not affecting the expected signal, lowers the variance at large angular separations (see also Refs. 2017MNRAS.471.3827S ; Gruen:2017xjj ).

We analyse the data in 12 logarithmically-spaced angular bins with radii between $5$ and $600\,\mathrm{arcmin}$, 9 photon energy bins between $0.631$ and $10^{3}\,\mathrm{GeV}$, and 4 redshift bins defined by $0.20<\langle z\rangle<0.43$, $0.43<\langle z\rangle<0.63$, $0.63<\langle z\rangle<0.90$ and $0.90<\langle z\rangle<1.30$, where $\langle z\rangle$ is the estimated expectation value of galaxy redshift from DES. The energy bins used in the analysis and the corresponding 68% and 95% containment angles of the Fermi-LAT PSF are shown in Tab. 1. These sum up to a total of 432 bins for the cross-correlation measurement. The analysis is performed blindly, i.e. on multiple variants of the measurements including artificial versions, in order to avoid experimental bias in measurement and interpretation of the signal. See the Supplemental Material for details.

The result of the measured cross-correlations, averaged over all energy and redshift bins, is shown in Fig. 2 in terms of the estimator $\Xi(\theta)$ defined in Eq. (1). Note that the data points reported on both panels are the same, although confronted with different models. A clear positive cross-correlation is observed, especially at small angular separations.

In order to determine the statistical significance of the signal, we test the deviation of the measurement from a null signal (null hypothesis of pure noise) by means of a phenomenological model, which aims at capturing the general expected features of the cross-correlation signal without resorting to any specific, detailed modelling of its physical origin (in the next Section, we will instead adopt a physical model to provide insights on the origin of the cross-correlation). In the halo-model approach, all mass in the large-scale structure of the Universe is associated with virialised dark matter halos, and the correlation function can thus be decomposed into the so-called 1-halo and a 2-halo terms (‘1h’ and ‘2h’, in formulae hereafter). The former refers to the correlation between two points in the same physical halo; the latter to the case in which the two points belong to two different halos. Point-like sources contribute at small angular scales with a 1-halo term, while at large scales they produce a 2-halo term resembling the large-scale structure matter distribution. In our case, we use the fact that the spatial extent of the 1-halo term is smaller than the beam window function of the Fermi-LAT. Then a phenomenological model can be constructed as:

[TABLE]

where $E_{a}$ and $z_{r}$ are the central values of the energy (measured in GeV) and redshift bins, and $\langle I_{a}\rangle$ is the measured photon flux. $\hat{\Xi}_{\textrm{PSF-like}}^{a}(\theta)$ is the Legendre transform of the beam window function (or point-spread function, PSF) integrated in the $a$th energy bin (in arbitrary units, being merely a template for the 1-halo term due to point-like $\gamma$-ray sources) and $\hat{\Xi}_{\textrm{2h-like}}^{ar}(\theta)$ is the Legendre transform of a generic 2-halo (i.e. large-scale) contribution, also convolved with the Fermi-LAT beam window function. Correlation functions with a hat have flux units, while those without a hat are normalised to the $\gamma$-ray flux as in Eq. (1), and therefore dimensionless. The two normalisations $A_{1}$ and $A_{2}$, spectral indices $\alpha_{1}$ and $\alpha_{2}$, and redshift evolution indices $\beta_{1}$ and $\beta_{2}$ are free parameters of the model 333For blazars, we expect an average spectral index of $2.2$, which is the reason for the term added to $\alpha_{2}$.. Gamma-ray sources typically have energy spectra that can be well approximated by a power-law, and so it is assumed in Eq. (2). For simplicity, we also assume a power-law scaling in redshift. Best fits and confidence intervals of the parameters are found in a Markov Chain Monte Carlo likelihood analysis.

The first statistical method adopted to quantify the presence of a signal, and its significance, against the null hypothesis relies on the $\Delta\chi^{2}$ test statistics, with the chi-squared defined in the usual way, i.e. :

[TABLE]

where ${\bm{\Xi}}_{\rm data}$ is the data vector, ${\bm{\Xi}}_{\rm th}$ is the theoretical cross-correlation for the models outlined above, described by the parameter set ${\bm{P}}_{\rm mod}$, and ${\bm{\Gamma}}$ is the data covariance matrix, detailed in the Appendix. (All angular, energy, and redshift bin indexes have been omitted for simplicity of notation.) The $\Delta\chi^{2}$ is defined as $\Delta\chi^{2}_{\rm mod}=\chi^{2}_{\rm null}-\chi^{2}({\bm{P}}_{\rm mod}^{\star})$, with $\chi^{2}({\bm{P}}_{\rm mod}^{\star})$ computed at the model parameter values ${\bm{P}}_{\rm mod}^{\star}$ that best fit the data, and $\chi^{2}_{\rm null}$ referring to no signal, i.e. ${\bm{\Xi}}_{\rm th}=0$. The second estimator is the matched filter signal-to-noise ratio (see e.g. Ref. Becker:2015ilr ),

[TABLE]

in analogy to $\Delta\chi^{2}_{\rm mod}$, we shall later refer to ${\rm SNR}_{\rm mod}\equiv{\rm SNR}({\bm{P}}_{\rm mod}^{\star})$.

In Table 2 we present the results on detection significance. The phenomenological model results for the full data set show a clear evidence for the presence of a cross-correlation signal, at the level of $\mathrm{SNR}_{\rm mod}=5.3$. Since the matched filter based on the phenomenological model captures the generic features of the cross-correlation signal, without committing to any specific physical description, this best assesses that indeed a cross-correlation between gravitational shear and unresolved $\gamma$-rays emission has been observed. In order to investigate the features of the signal in more detail, we repeat the tests by subdividing the data set according to redshift, energy, and angular separation. Specifically, Low/High-$z$ refers to the first two and second two redshift bins; Low/High-$E$ bins are defined by being below/above $5\,\mathrm{GeV}$; and Small/Large-$\theta$ separates angular scales below/above 3 times the Fermi-LAT PSF.

From Table 2 we infer that the signal is mostly concentrated at high energies and small angles. These results point towards an interpretation in terms of point-like sources with hard energy spectrum, broadly compatible with these sources being blazars. In fact, the best-fit for the spectral index of the PSF-like 1-halo component $\alpha_{1}=1.81^{+0.20}_{-0.24}$ is quite hard with respect to the spectral index of the average intensity of the UGRB Ackermann:2014usa , but compatible with BL Lac emission, which is the source population expected to be the most relevant in the range of fluxes probed by this analysis, just below the Fermi-LAT flux sensitivity threshold. Notice that this hard spectral index is in agreement with recent findings from $\gamma$-ray auto-correlation analysis Ackermann:2018wlo , possibly suggesting that BL Lac objects below the threshold have slightly harder spectra than those detected individually. The energy scaling of the 2-halo component is also compatible with a blazar origin, though this term shows lower statistical significance than the 1-halo component. Concerning the redshift dependence of the signal, the statistical significance is almost equally distributed among the lower and higher redshift bins. The allowed regions for the parameters of the *phenomenological * model are shown in Fig. 3, while the cross-correlation function for the best-fit of the phenomenological model are shown in the left panel of Fig. 2: the PSF-like term due to point-like sources well reproduces the behaviour of the measured cross-correlation up to about $1\,\deg$ scale. We note here that for the subset of High-$E$/Small-$\theta$, comprising 88 data-points, we do obtain a distinctive signal without application of the matched filter. The $\chi^{2}_{\rm null}=137$ for these points corresponds to a p-value of 0.0006, meaning that the null hypothesis is excluded at $3.5\sigma$ in this subset.

III Discussion

In the following we attempt a physical interpretation of the signal detected in the previous Section. Star-forming galaxies and misaligned AGNs are not expected to be able to produce a sufficiently hard energy spectrum, which thus points to a dominant blazar component. Particle dark matter in terms of WIMPs can also provide a hard spectrum, especially if the annihilation channel is predominantly leptonic or, in the case of a hadronic final state, if the dark matter mass is large (above a few hundred GeV).

Blazars are compact sources and, for our purposes, they can be considered as point-like—i.e. their size is, on average, much smaller than the Fermi-LAT PSF. Also the size of the halo hosting blazars rarely exceeds the Fermi-LAT PSF. This has a consequence that the angular correlation function for the 1-halo term essentially follows from the detector PSF. Conversely, the relevant dark matter halos have a larger angular extent, and the corresponding 1-halo correlation function thus drops more slowly with angular scale. On very large scales, the correlation functions of the two components have a similar angular behaviour, since the 2-halo power spectra differ only by the bias terms. The fact that our signal is detected with high significance only on small scales therefore points towards blazars as the dominant source. In order to investigate this interpretation, we perform the statistical tests discussed in the previous Section with a physical model, based on a detailed characterisation of the components expected to produce the cross-correlation signal: blazars (BLZ), misaligned active galactic nuclei (mAGN), star-forming galaxies (SFG) and possibly dark matter (DM). The physical cross-correlation function model reads:

[TABLE]

The model parameters are: free normalisations for the astrophysical sources, $A_{\rm BLZ}^{\rm 1h}$, $A_{\rm BLZ}^{\rm 2h}$, $A_{\rm mAGN}$, and $A_{\rm SFG}$; the mass of the dark matter particle, $m_{\rm DM}$; its velocity-averaged annihilation rate, $\langle\sigma_{\rm ann}v\rangle$, expressed in terms of the “thermal” cross-section $\langle\sigma_{\rm ann}v\rangle_{\rm th}=3\times 10^{-26}\,{\rm cm^{3}\,s^{-1}}$ through the normalisation $A_{\rm DM}\equiv\langle\sigma_{\rm ann}v\rangle/\langle\sigma_{\rm ann}v\rangle_{\rm th}$. Note that for blazars, which represent the astrophysical component expected to dominate the correlation signal at the current level of unresolved $\gamma$-ray emission, we allow the 1-halo and the 2-halo terms to be separately adjusted in the fit against the data. As for the phenomenological model, all terms depend on both energy and redshift, labelled by indices $a$ and $r$, respectively.

The results are shown in Table 2, where the overall significance of the presence of a signal, the preference for an origin at high energies and small angular scales are all confirmed. However, since in this case we have specific behaviours for the correlation functions as dictated by a physical model (different for each component, contrarily to the average generic case of the phenomenological model), we notice that a mild hint of large scale correlation is present—namely, in the Large-$\theta$ case. We note that both the physical and phenomenological models provide a good fit to the data according to their $\chi^{2}$ (see the Supplemental Material).

More details of the analysis are shown in Fig. 3, where the triangular plot of the profile likelihood distributions of the model parameters is reported. The likelihood exhibits a preference for a large 1-halo term of blazars with normalisation $A_{\rm BLZ}^{\rm 1h}=102^{+56}_{-57}$, while the normalisations of the blazar 2-halo term shows only a (weak) upper bound. The latter picture is shared also by the other astrophysical sources (SFG and mAGN) which are shown only in the Appendix for brevity.

The blazar-shear cross-correlation on small scales depends on the relation between the blazar $\gamma$-ray luminosity and the host-halo mass, a quantity which is rather uncertain. For our reference model this relation has been taken from Camera:2014rja , where it was derived by associating the $\gamma$-ray luminosity of blazars to the mass of the supermassive black hole powering the AGN and then relating the mass of the black hole to the mass of the dark matter halo. This procedure gives $M(\mathcal{L})=2\times 10^{13}M_{\odot}\left[\mathcal{L}/(10^{47}\,\mathrm{erg\,s^{-1}})\right]^{0.23}(1+z)^{-0.9}$, where $\mathcal{L}$ is the rest-frame luminosity of blazars in the energy range $0.1$ to $100\,\mathrm{GeV}$. We can therefore translate a value of $A_{\rm BLZ}^{\rm 1h}$ different from unity to a deviation from the reference $M(\mathcal{L})$ relation. The value we found implies that the average mass of a halo hosting an unresolved blazar is larger than the one adopted in Ref. Camera:2014rja , and most likely above $10^{14}M_{\odot}$. The cross-correlation signal with weak lensing seems therefore to be dominated by blazars residing in cluster-size halos.

The right panel of Fig. 2 shows that the cross-correlation at small angular scales requires a sizeable blazar 1-halo term. It also illustrates that the best-fit analysis exhibits a mild preference for some power at large scales. This can be accounted for either by the 2-halo term of blazars or by a DM contribution. The interplay of the angular, energy and redshift behaviours of the observed signal leads to a small preference for a DM component over a pure blazar contribution in the assumed model. Misaligned AGNs and star-forming galaxies are disfavoured since they do not meet the requirement of a hard energy spectrum.

To have a visual impression on the physical behaviours, we plot the energy and redshift dependence of the cross-correlation signal in Fig. 4. The average along the angular and redshift (energy) directions of each point of the energy (redshift) spectrum is performed by computing a matched filter amplitude ${\cal A}={\bm{\Xi}}^{\sf T}{\bm{\Gamma}}^{-1}\bar{\bm{\Xi}}_{\rm M}/(\bar{\bm{\Xi}}_{\rm M}^{\sf T}{\bm{\Gamma}}^{-1}\bar{\bm{\Xi}}_{\rm M})$, where $\bar{\bm{\Xi}}_{\rm M}$ is given by a sample model that we choose to be flat in energy and redshift, while scaling as $1/\theta$ in angle, to approximately reproduce the expected signal, and ${\bm{\Xi}}$ is either the measured data or the best-fit models introduced in the main text. The error on ${\cal A}$ is given by $\sigma^{2}_{\cal A}=({\bm{\Gamma}}^{-1}\bar{\bm{\Xi}}_{\rm M})^{\sf T}{\bm{\Gamma}}({\bm{\Gamma}}^{-1}\bar{\bm{\Xi}}_{\rm M})/(\bar{\bm{\Xi}}_{\rm M}^{\sf T}{\bm{\Gamma}}^{-1}\bar{\bm{\Xi}}_{\rm M})^{2}$. From Fig. 4 one can appreciate the small but noticeable differences in the energy and redshift scalings of the models of different physical components that have been just discussed.

Notice that the blazar model we are adopting, and which is outlined in the Supplemental Material, is based on the current understanding of blazars as derived from the Fermi-LAT resolved sources: the small preference in the fit for a contribution with features compatible with DM might be interpreted as an indication that unresolved blazars have different properties than the resolved ones.

We find that for a dark matter particle dominantly annihilating into the leptonic channel $\tau^{+}\tau^{-}$ the best-fit improves by $2.8\sigma$ as compared to a model where DM is not included. The lack of a degeneracy between dark matter and blazar amplitudes (see Fig. 3) indicates that the two components are supported by independent features of our cross-correlation data, in particular the small and large scale behaviour. The best-fit occurs for a dark matter particle of mass $m_{\rm DM}=(65\,\pm\,_{23}^{27})\,\mathrm{GeV}$ and annihilation rate $\langle\sigma_{\rm ann}v\rangle=(26\,\pm\,_{15}^{17})\times\langle\sigma_{\rm ann}v\rangle_{\rm th}$. In the case of a softer energy spectrum, as the one provided by a $\bar{b}b$ annihilation channel, the fit improvement is slightly lower, at $2.7\sigma$ level, with best-fit mass $m_{\rm DM}=(302\,\pm\,_{120}^{188})\,\mathrm{GeV}$ and annihilation rate $\langle\sigma_{\rm ann}v\rangle=(78\,\pm\,_{43}^{67})\times\langle\sigma_{\rm ann}v\rangle_{\rm th}$.

Let us remark that the main source of uncertainty concerning the dark matter signal described in this work is our ignorance on the impact of substructures. Specifically, the minimal halo mass and the amount and distribution of sub-haloes can significantly change the size of the signal. This, however, is common to all cosmological searches for a particle dark matter signature. Comparing our nominal model Moline2016 to a recent development Hiroshima:2018kfv based on both $N$-body simulations and analytical modelling, we found that the constraint on the annihilation cross-section gets shifted to higher annihilation rates by about one order of magnitude. This is due to a different amplitude of the expected signal, whereas the energy, redshift, and angular dependencies are only slightly modified, and the predicted cross-correlation function therefore just needs a larger normalisation, which in turn is directly provided by the annihilation rate.

We conclude by summarizing the major results of this analysis. We present the first detection of the cross-correlation between the $\gamma$-ray sky and the mass distribution in the Universe observed through gravitational lensing shear, with a significance of SNR$=5.3$. The bulk of this signal comes from the 1-halo term of point-like sources with a hard spectrum, most likely dominated by blazars. In addition, we find a hint for a cross-correlation on large scales with spectral and redshift behaviours that might imply that the population of blazars that are currently unresolved by Fermi-LAT has different characteristics than those obtained by extrapolation from observations of resolved blazars, or that an additional contributor to the $\gamma$-ray emission is present. The analysis of the cross-correlation of Fermi-LAT data with the forthcoming Year 3 and Year 5 DES data set, and improvements in modelling of the blazar population will likely clarify the source of the signal detected in this work, and characterise it more deeply.

IV Acknowledgments

We warmly acknowledge Mattia Fornasa, for his collaboration in the early phases of this project. This work is supported by: ‘Departments of Excellence 2018-2022’ grant awarded by the Italian Ministry of Education, University and Research (miur) L. 232/2016; Research grant ‘The Anisotropic Dark Universe’ No. CSTO161409, funded by Compagnia di Sanpaolo and University of Turin; Research grant TAsP (Theoretical Astroparticle Physics) funded infn; Research grant ‘The Dark Universe: A Synergic Multimessenger Approach’ No. 2017X7X85K funded by miur; Research grant ‘From Darklight to Dark Matter: understanding the galaxy/matter connection to measure the Universe’ No. 20179P3PKJ funded by miur; Research grant ‘Deciphering the high-energy sky via cross correlation’ funded by Accordo Attuativo ASI-INAF n.2017-14-H.0. This work was supported by the Department of Energy, Laboratory Directed Research and Development program at SLAC National Accelerator Laboratory, under contract DE-AC02-76SF00515 and as part of the Panofsky Fellowship awarded to DG. Support for DG was also provided by Chandra Award Number GO8-19101A, issued by the Chandra X-ray Observatory Center. SC is supported by miur through Rita Levi Montalcini project ‘prometheus – Probing and Relating Observables with Multi-wavelength Experiments To Help Enlightening the Universe’s Structure’. SA was supported by JSPS KAKENHI Grant Numbers JP17H04836, JP18H04578, and JP18H04340. This work was supported in part by the U.S. Department of Energy under contract number DE-AC02-76SF00515.

Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey.

The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l’Espai (IEEC/CSIC), the Institut de Física d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, Texas A&M University, and the OzDES Membership Consortium.

Based in part on observations at Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.

The DES data management system is supported by the National Science Foundation under Grant Numbers AST-1138766 and AST-1536171. The DES participants from Spanish institutions are partially supported by MINECO under grants AYA2015-71825, ESP2015-66861, FPA2015-68048, SEV-2016-0588, SEV-2016-0597, and MDM-2015-0509, some of which include ERDF funds from the European Union. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya. Research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Program (FP7/2007-2013) including ERC grant agreements 240672, 291329, and 306478. We acknowledge support from the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), through project number CE110001020, and the Brazilian Instituto Nacional de Ciência e Tecnologia (INCT) e-Universe (CNPq grant 465376/2014-2).

This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.

This paper has gone through internal review by the DES collaboration.

Appendix A Gravitational lensing data

The Dark Energy Survey (DES 2016MNRAS.460.1270D ) is a photometric survey performed with the Dark Energy Camera (DECam 2015AJ….150..150F ) on the Blanco $4\,\mathrm{m}$ telescope at Cerro Tololo Inter-American Observatory (CTIO) in Chile. DES has observed approximately $5000\,\deg^{2}$ of the Southern sky in five broadband filters, g, r, i, z, and Y, ranging from $400$ to $1060\,\mathrm{nm}$. The primary goal of DES is to study the origin of cosmic acceleration and the nature of dark energy through four key probes: weak lensing, clustering of the large-scale structure, galaxy clusters, and type Ia supernovae.

The first season of DES observations, from 31${}^{\textrm{st}}$ August 2013 to 9${}^{\textrm{th}}$ February 2014, resulted in the Y1 data stet analysed here Diehl:2014lea ; y1gold . During Y1, DES observed $\sim 1500\,\deg^{2}$ of the wide-area survey footprint with three to four dithered tilings per filter. The Y1 footprint consisted of two areas: one near the celestial equator including Stripe 82 (S822014ApJ…794..120A ), and a much larger area that was also observed by the South Pole Telescope (SPT2011PASP..123..568C ). In Fig. 1 of the main text we show with a solid white line the DES Y1 sky coverage in celestial coordinates; the complete sky coverage for the final survey is shown with a dashed line. Fig. 5 shows the DES Y1 footprint in Galactic coordinates.

In our analysis we make use of the metacalibration shear catalogue Zuntz:2017pso . The metacalibration catalogue yields a total of 35 million objects and the final number density of the selection is 5.5 galaxies per square arcminute. metacalibration Huff:2017qxu ; Sheldon:2017szh is a method to estimate weak lensing shear, calibrating associated biases directly from the imaging data, rather than from prior information about galaxy properties or simulations. metacalibration has been shown to be accurate at the part-per-thousand level in the absence of blending with other galaxies Sheldon:2017szh , and at the part-per-hundred level for the blending present in DES Y1 data Zuntz:2017pso .

The implementation of metacalibration consists of measuring the two-component ellipticity, ${\bm{e}}$, of a galaxy in the DES Y1 run by fitting a single Gaussian model to its single-epoch images in the riz bands. The galaxy images are then artificially distorted with a known shear, ${\bm{\gamma}}$, and the ellipticity is measured again to construct the response matrix, ${\bm{R}}_{\gamma}$, as the derivative of measured ellipticity w.r.t. shear. Thus, the ellipticity estimator can be written as the following Taylor expansion Huff:2017qxu ; Sheldon:2017szh ; Zuntz:2017pso :

[TABLE]

The elements of the response matrix are measured in metacalibration by the numerical derivative:

[TABLE]

where $e^{\pm}_{i}$ is the measurement of the $i$th ellipticity component made on an image sheared by $\pm\gamma_{j}$, and $\Delta\gamma_{j}=2\gamma_{j}$.

It is also necessary to correct for selection effects, i.e. shear biases that may occur when placing a cut on any quantities that change under shear, such as signal-to-noise ratio. This is accomplished by measuring the mean response of the estimator to the selection, repeating the selections on quantities measured on sheared images. The mean selection response matrix is given by:

[TABLE]

where $\langle e_{i}^{S\pm}\rangle$ represent the mean of the ellipticities measured on images without artificial shearing, but selected by properties measured on sheared images. The full response for the mean shear is given by the sum of the shear response and selection response:

[TABLE]

When measuring a shear statistics, such as mean shear or a shear two-point function, the measurement can be appropriately corrected for the mean response $R=\left(\langle\bm{R}\rangle_{11}+\langle\bm{R}\rangle_{22}\right)/2$ to produce a calibrated result.

DES galaxies were assigned to redshift bins using a re-implementation of the bpz algorithm Benitez2000 , which provides an estimate of the redshift probability density $p_{\rm BPZ}(z_{j})$ of each galaxy $j$. The assignment is done based on fluxes measured by metacalibration on the original and sheared images to correct for selection response of redshift binning. The fluxes in the bands griz of the metacalibration galaxies are used for estimating an expectation value of $\langle z\rangle=\int_{0}^{z}z_{j}{\rm d}z_{j}\,p_{\rm BPZ}(z_{j})$. Thus, four redshift bins are defined by $0.20<\langle z\rangle<0.43$, $0.43<\langle z\rangle<0.63$, $0.63<\langle z\rangle<0.90$, and $0.90<\langle z\rangle<1.30$. The corresponding four redshift distributions, $n^{i}(z)$, are taken as the stacked $p_{\rm BPZ}(z)$ of the galaxies, estimated from their improved joint-fit photometry y1gold ; Hoyle:2017mee . The mean values of each redshift bin stacked $p_{\rm BPZ}(z)$ are calibrated independently, and they are found to be consistent with the bpz estimate, with a joint uncertainty of $\sigma_{\langle z\rangle}\approx 0.015$ Hoyle:2017mee .

Since the purpose of this work is detection of the cross-correlation between the UGRB and shear, rather than accurate constraints on cosmological parameters, we do not explicitly account for the systematic uncertainties in shear and redshift calibration in our analysis. The moderate significance of our measurements and the accuracy of the DES Y1 calibrations justify this choice.

Appendix B Gamma-ray data

Fermi-LAT is a $\gamma$-ray pair-conversion telescope that has been operating for the last 10 years in space. Due to its wide energy range ($20\,\mathrm{MeV}$ to $1\,\mathrm{TeV}$) and its capability of rejecting the background of charged cosmic rays, it is an excellent instrument to investigate the UGRB. It scans the whole sky every three hours with a remarkable angular resolution for $\gamma$ rays ($\sim 0.1\,\deg$ above $10\,\mathrm{GeV}$).

In this work we used 108 months of data, from mission week 9 to week 476. The photon and exposure maps are produced with the LAT Science Tools version v10r0p5 444 https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/ . We selected the Pass8 ultracleanveto class 555See http://www.slac.stanford./glast/groups/canda/lat_Performance.htm , which has the lowest cosmic-ray contamination and is the most appropriate class for diffuse emission analysis. The Fermi Tools provide the possibility of choosing different angular resolutions, which are organised in four quartiles, from PSF0 to PSF3, corresponding to a transition from the worst to the best PSF. In order to have a balance between the photon count statistics and a good direction reconstruction, we selected the best quartile PSF3 for energies below $1.2\,\mathrm{GeV}$ (where we have the highest photon counts) and PSF1+2+3 for higher energies. The PSF is modelled according the Fermi-LAT specifications 666 https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_LAT_IRFs/index.html and for each energy bin in our analysis an effective PSF is determined by weighting the energy-dependent PSF by the intensity energy spectrum of the UGRB.

We produced 100 intensity maps in HEALPix projection with $N_{\rm side}=1024$, evenly spaced in logarithmic scale between $100\,\mathrm{MeV}$ and $1\,\mathrm{TeV}$, by dividing the count maps by the exposure and the pixel area. The size of such energy bins is small enough that the exposure can be approximated by its mean value within the energy bin when deriving the flux. The resulting flux maps are then re-binned into 9 larger energy bins between $631\,\mathrm{MeV}$ to $1\,\mathrm{TeV}$ by simply adding up fluxes from the smaller bins. We discard very low energies because the angular resolution is too poor for our purposes.

Since we are interested only in the UGRB, we benefit from excluding Galactic emission and resolved point sources. This is achieved by a process of masking and subtracting described below.

B.1 Masking $\gamma$-ray data

We build a set of masks according to the following two criteria:

Low latitudes, where the Galactic foreground is stronger, are removed by a flat cut excluding the region between $\pm 30\,\deg$ of latitude. 2. 2.

Sources identified in the list FL8Y are masked. FLY8 777 https://fermi.gsfc.nasa.gov/ssc/data/access/lat/fl8y/ has been recently released by the Fermi-LAT Collaboration as a preliminary version of the upcoming 4FGL catalogue. It contains 5523 sources. Above 10 GeV, we mask also the sources which are present in the 3FHL catalogue AjeloEtAll2017 , that is more accurate for high energy sources. Each source is masked taking into account both its source brightness and the detector PSF resolution in the specific energy bin. The masking radius $R$ has been defined by:

[TABLE]

where $F_{\Delta E}^{g}$ is the integral flux of the source in a given energy bin $\Delta E$, $F_{\Delta E,\rm faintest}^{g}$ is the flux of the faintest source in the same energy bin (and once divided by 5 provides an approximate estimate of the noise, i.e., the faintest source emission is approximately a measure of the $5\sigma$ level), and $\theta_{\Delta E}$ is the 68% containment angle in that energy bin, as provided by the Fermi-LAT PSF. We verified that the non-Gaussian tail of the PSF (not included in the Gaussian approximation in Eq. (10)) does not appreciably contaminate our maps.

This strategy aims at masking the Galactic plane and resolved sources over a sufficiently large area, in order to reduce the chance to have artefacts in the APS produced by source leakage and foreground emission. For further details and impact of the mask, see also ref. Ackermann:2018wlo .

B.2 Foreground removal

Despite the majority of the Galactic foreground has been removed by cutting out the Galactic plane, still high-latitudes might contain some contamination that causes noise in our analysis. While galactic emission does not contribute to the cross-correlation signal with gravitational shear, nevertheless a foreground contamination adds to the error budget. We therefore performed a galactic foreground subtraction by adopting the template maps provided by the Fermi-LAT Collaboration with the Galactic emission model gll_iem_v06.fits 888 https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html . The foreground template is projected in HEALPix maps with the same $N_{side}$ as the intensity maps and in the same 100 energy bins. Each template map is convolved with the Fermi-LAT PSF and assigned a free normalization. This component is added to a free constant (representing the UGRB and cosmic-ray contamination) and a Poissonian likelihood fit is performed globally on all the masked intensity maps. All obtained best-fit normalization parameters are of the order of unity, supporting a successful description of the foreground emission. The normalised foreground templates are then re-binned into the 9 larger energy bins and subtracted from the corresponding intensity maps.

In Fig. 5, we show an example of the Fermi-LAT $\gamma$-ray intensity map in the (1.2-2.3) GeV energy bin with the application of the mask described above and illustrating the effect of the Galactic foreground subtraction.

The impact of foreground removal in cross-correlation studies involving $\gamma$-rays has been discussed in Refs. Xia:2015wka ; Shirasaki2014 ; Shirasaki:2016kol ; Ammazzalorso:2018evf ; Shirasaki:2018dkz , where it was shown that the effect is marginal. As mentioned, the main effect of foreground residuals on our study would be a mild reduction of detection significance. This can be appreciated in Fig. 6, where we show the analogous of Fig. 2 of the main text, but without foreground removal. Since the amplitude $\Xi$ is dimensionless, in Fig. 6 it has been normalized as in Fig. 2 of the main text (i.e. relative to the mean intensity $\langle I_{\gamma}\rangle$ of the foreground-subracted maps), to allow a direct comparison of the two plots. Put in different words,Supplemental Material

Appendix C Estimation of the covariance matrix

The covariance matrix is estimated by combining the theoretical Gaussian large-scale structure covariance with realizations of the shape-noise term generated from the data. The latter term is expected to be dominating in the covariance, while the former should be small. To avoid notation conflict between the symbol $\gamma$ usually employed for both the weak lensing effect of gravitational shear, and for $\gamma$-rays, we shall here simply refer to shear through indexes $r,s$, labelling redshift bins, and to the UGRB anisotropies via indexes $a,b$, labelling energy bins.

In the Gaussian approximation, the generic element of the theoretical harmonic-space covariance matrix $\widehat{\bm{\Gamma}}$ reads:

[TABLE]

It represents the covariance between the measurements of two cross-correlation harmonic-space power spectra: one cross-correlating $\gamma$-ray map in energy bin $a$ with shear map in redshift bin $r$; and another considering $\gamma$-ray map in energy bin $b$ with shear map in redshift bin $s$. In the Gaussian approximation, these two cross-correlation power spectra, $C^{ar}_{\ell}$ and $C^{bs}_{\ell^{\prime}}$, have a covariance which is diagonal in $\ell\ell^{\prime}$, and which depends on both the cross-correlation signals, as well as the $\gamma$-ray–$\gamma$-ray and shear-shear auto-correlations, i.e. $C^{ab}_{\ell}$ and $C^{rs}_{\ell^{\prime}}$. Both auto and cross-correlation theoretical signals involving $\gamma$-rays have been corrected for the effect of the Fermi-LAT PSF beam function. In Eq. (11), $\mathcal{N}$ denotes the noise terms. Note that the noise does not depend on the angular scale. In Eq. (11) $f_{\rm sky}$ accounts for the incomplete sky coverage, with $f_{\rm sky}^{\rm DES}=0.042$ denoting the DES footprint (independent of the redshift bin) and $f_{\rm sky}^{a}=(0.315,0.289,0.404,0.449,0.468,0.475,0.476,0.477,0.477)$ being the fraction of unmasked pixels of the $\gamma$-ray map in the $a$-th energy bin. For the cross-correlation estimator with complex masks as we have here, an unambiguous definition of $f_{\rm sky}$ is not possible. We tested two options: a geometric mean $f_{\rm sky}=(f_{\rm sky}^{a}f_{\rm sky}^{\rm DES})^{1/2}$ and the $f_{\rm sky}$ relative to the overlap of the DES footprint with the Fermi-LAT unmasked sky. We comment more in the two alternatives below, but we anticipate that results on the SNR arising from the two choices exhibit only a minor change.

For details on the Gaussian covariance matrix, see Section 5.1 in ref. Camera:2014rja . Since our measurements are performed in real space, we compute the Legendre transform of Eq. (11) to return the real-space covariance matrix, $\bm{\Gamma}$, with entries $\Gamma_{ar\vartheta,bs\vartheta^{\prime}}$.

As mentioned above, the dominant contribution to the covariance is represented by the shape noise term. In the Gaussian approximation, it takes the form of $\mathcal{N}^{rs}(C_{\ell}^{ab}+\mathcal{N}^{ab})$. In order to estimate this contribution more accurately, without resorting to theoretical approximations, we produce 2000 realisations of the noise directly from the data in real space. To this end, we rotate each galaxy in our catalogue by an independent random angle. The shear signal as measured from these rotated source catalogues represents a random realisation of shape noise (see e.g. refs. Gruen:2017xjj ; Hikage:2018qbn ; Murata:2017zdo ).

We complete the construction of the total covariance by adding a theoretical estimate of the large-scale structure term, $C_{\ell}^{ar}C_{\ell}^{bs}+C_{\ell}^{rs}(C_{\ell}^{ab}+\mathcal{N}^{ab})$. We assume it subdominant in the covariance budget, and this is verified a posteriori. The $\gamma$-ray auto-correlation $C_{\ell}^{ab}$ entering in the theoretical estimate can be well-fitted by a simple model, given by a power-law plus a constant, i.e. $C_{\ell{\rm,mod}}^{ab}=A_{ab}\ell^{-\alpha_{ab}}+C_{\rm P}^{ab}$ (see e.g. ref. Ackermann:2018wlo ). We fit the three parameters $\{A_{ab},\,\alpha_{ab},\,C_{\rm P}^{ab}\}$ for each energy bin pair against the measurement of the auto-correlation power spectrum obtained by running the tool PolSpice (as in ref. Ackermann:2018wlo ). Finally, the $\gamma$-ray noise term, $\mathcal{N}^{ab}$, is computed with the analytic estimator of eq. 5 in ref. Fornasa:2016ohl .

The shear auto-correlation is derived using the galaxy redshift distributions described above, and adopting $\Lambda$CDM cosmology with parameter fiducial values from Planck Planck2015 . The model is compatible with previous measurements Troxel2018 .

The large-scale structure part of the covariance is added to the shape noise term by using the following procedure. We create a set of 2000 simulated datasets from a multivariate Gaussian distribution with zero mean and covariance equal to the large-scale structure part of the theoretical covariance. Then, we add this simulated data to the shape noise realisations discussed above, thus obtaining 2000 samples containing both terms. These 2000 samples are then used to obtain the covariance matrix. The inverse of the covariance estimated such is a biased estimator of the inverse covariance, with the bias depending on the number of realisations, number of bins, and parameters of the model chosen to fit the data anderson2009 ; Hartlap2007 . We apply the Anderson-Hartlap correction in order to de-bias the inverse covariance (as done e.g. in Section 3 of ref. Kilbinger2013 ).

In order to validate our procedure, we compare the shape noise term obtained via the simulations with the analogous term from the theoretical estimate, $\mathcal{N}^{rs}(C_{\ell}^{ab}+\mathcal{N}^{ab})$. In Fig. 7, we show a comparison between the theoretical shape noise standard deviation and the results from the simulations, in the case of a low-energy bin (left panel) and a high-energy bin (right panel) combined with the lowest DES redshift bin. We see that the variance obtained with the two techniques approximately agrees, yet can differ due to complex masking effects present in the data Troxel2018b . When deviations are present, the variance obtained with simulations is typically larger, as expected. By adopting the shape noise term directly derived from the data, we ensure that the quoted errors and values for the goodness of fit are correct.

We compare also the large-scale structure term with the covariance matrix obtained from simulations after rescaling the latter by $C_{\ell}^{rs}/\mathcal{N}^{rs}$. Results are in very good agreement with the choice of the effective parameter $f_{\rm sky}$ entering Eq. (11) as the geometric mean of DES and Fermi-LAT sky coverage. Another value for the $f_{\rm sky}$ parameter that is frequently used in the literature is the overlap between the two masks. We verified that, in this case, results would be just slightly modified, with the SNR becoming 4.8 instead of 5.3, as obtained in the main analysis.

Appendix D Theoretical models

The harmonic-space (often referred to as ‘angular’) power spectrum of the cross-correlation between a map of $\gamma$-rays in the $a$th energy bin and a map of gravitational shear in the $r$th redshift bin can be computed as:

[TABLE]

where $E$ is the $\gamma$-ray energy, and $\chi(z)$ is the comoving distance to redshift $z$, obeying $\mathrm{d}z/\mathrm{d}\chi=H(z)$ with $H(z)$ the Hubble rate. Formally, the integrals extend over the whole $\gamma$-ray energy spectrum and from $z=0$ to the horizon, but the window functions $W_{\textrm{\rm gamma}}^{a}(E,z)$ and $W_{\rm shear}^{r}(z)$ implement energy and redshift binning effectively reducing the integration range (for details, see next Subsections). Finally, $P_{\gamma\delta}$ is the three-dimensional cross-power spectrum between a given $\gamma$-ray population sourcing the UGRB emission and the matter density contrast $\delta$. It represents the three-dimensional correlation (in Fourier space) between what seeds the unresolved $\gamma$-ray emission and what sources the gravitational lensing effect, namely matter inhomogeneities. It is a function of both redshift and physical scale $k$, the modulus of the physical wavenumber. In the Limber approximation, $k$ and the angular multipole $\ell$ are linked by $k=\ell/\chi(z)$. This approximation is valid for $\ell\gg 1$, which is the case of the present work.

As mentioned before, the quantity measured from the data is the real-space cross-correlation of UGRB anisotropies with tangential shear, which is related to the harmonic-space cross-power spectrum of Eq. (12) by a Legendre transform:

[TABLE]

with $\theta$ being the angular separation on the sky, $P^{(2)}_{\ell}$ the Legendre polynomial of order two, and $W_{\ell}^{a}$ the Fermi-LAT PSF beam function in energy bin $a$, accounting for the finite resolution of the detector.

D.1 Gravitational lensing window function

The window function describes the mean distribution of the signal along the line of sight, in the given energy or redshift bin. For the gravitational lensing effect, the window function is given by (see e.g. ref. Bartelmann2010a ):

[TABLE]

where $H_{0}\equiv H(z=0)$ is the Hubble constant, $\Omega_{\rm m}$ is the matter abundance in the Universe (sum of the dark matter and the baryon abundances, $\Omega_{\mathrm{DM}}$ and $\Omega_{\mathrm{b}}$), and $n^{r}(z)$ is the redshift distribution of background galaxies in the lensing data set in bin $r$. The galaxy distribution depends on the data set and redshift cut, as described above.

D.2 WIMP-sourced $\gamma$-rays window function

The window function for UGRB anisotropies sourced by annihilations of dark matter particles reads Ando2006 ; Fornengo:2014 :

[TABLE]

where $\rho_{\rm c}$ is the critical density of the Universe, $m_{\rm DM}$ is the mass of the dark matter particle, and $\langle\sigma_{\rm ann}v\rangle$ denotes the velocity-averaged annihilation cross-section, assumed here to be the same in all dark matter haloes.

Among the other ingredients, we have $\mathrm{d}N_{\rm ann}/\mathrm{d}E$, indicating the number of photons produced per annihilation as a function of energy, and setting the $\gamma$-ray energy spectrum. We will consider it to be given by the sum of two contributions: prompt $\gamma$-ray production from dark matter annihilations (or decays); and inverse Compton scattering of dark matter produced electrons and positrons on CMB photons (we compute inverse Compton assuming negligible magnetic field and no diffusion for the produced electrons and positrons). Results of our analysis will be shown for two annihilation final states: $b\bar{b}$ pairs, which yield a relatively soft spectrum of photons and electrons, mostly associated to hadronisation into pions and their subsequent decay; $\tau^{+}\tau^{-}$, which provides a harder spectrum, mostly associated to final state radiation of photons and decay of the muons produced by $\tau$ decay into electrons (with subsequent $\gamma$-ray emission through inverse Compton on CMB), with an additional semi-hadronic decay into pions Fornengo2004 ; Cembranos2011 ; Cirelli2011 .

Then, the optical depth $\tau$ in Eq. (15) accounts for attenuation of $\gamma$-rays by the extra-galactic background light, and is taken from the literatureFinke2010 . Finally, the clumping factor $\Delta^{2}$ is related to the fact that the signal of annihilating dark matter depends on the average of the square of the dark matter density; it is defined as (see e.g. ref. Fornengo:2014 ):

[TABLE]

where: $\mathrm{d}n_{\rm h}/\mathrm{d}M$ is the halo mass function Sheth1999 ; $M_{\rm min}$ is the minimal halo mass, which we fix to $10^{-6}M_{\odot}$; $M_{\rm max}$ is the maximal mass of haloes, for which we use $10^{18}M_{\odot}$, although results are insensitive to the precise value assumed; $\rho_{\rm h}({\bm{x}}|M,z)$ is the dark matter density profile of a halo with mass $M$ at redshift $z$, taken to follow a Navarro-Frenk-White profile Navarro1997 ; and $b_{\rm sub}$ encodes the ‘boost’ to the halo emission provided by subhaloes. To characterise the halo profile and the subhalo contribution, we need to specify their mass concentration. For the main haloes we follow ref. Prada2012 . On the other hand, the description of the concentration parameter $c(M,z)$ at small masses and for subhaloes is still an open issue and provides a source of uncertainty. We considered the two most updated analyses, from Ref. Moline2016 , where we assume $\mathrm{d}n_{\rm sub}/\mathrm{d}M_{\rm sub}\propto M_{\rm sub}^{-2}$ (see also Ref. Sanchez-Conde2014 ), and from Ref. Hiroshima:2018kfv . We show only the former, since the two descriptions just differ by roughly a constant factor (around 10) in the derived final value of the annihilation cross section.

D.2.1 Astrophysical $\gamma$-ray sources window function

Besides the possible particle dark matter emission, $\gamma$-rays in the UGRB are certainly produced by a number of astrophysical sources. The most relevant $\gamma$-ray emitters include: blazars, misaligned active galactic nuclei, and star forming galaxies. Their window function is given by:

[TABLE]

where $\mathcal{L}$ is the $\gamma$-ray rest-frame luminosity in the energy interval $0.1$ to $100\,\mathrm{GeV}$, $\Phi_{\rm S}$ is the $\gamma$-ray luminosity function (GLF) of the source class S of astrophysical emitters included in our analysis, and $\mathrm{d}N_{\rm S}/\mathrm{d}E$ is its observed (unabsorbed) energy spectrum. The upper bound, $\mathcal{L}_{\rm max}(F_{\rm sens},z)$, is the luminosity above which an object is resolved in the FL8Y and 3FHL catalogues, and consequently masked in our analysis. As we are interested in the contribution from unresolved astrophysical sources, only sources with luminosity smaller than $\mathcal{L}_{\rm max}$ are included. Conversely, the minimum luminosity, $\mathcal{L}_{\rm min}$, depends on the properties of the source class under consideration.

We consider a unified blazar model combining BL Lacertae and flat-spectrum radio quasars as a single source class. The GLF and energy spectrum are taken from ref. Ajello2015 where they are derived from a fit to the properties of resolved blazars in the third Fermi-LAT catalogue Acero2015 . In the case of misaligned AGNs, we follow ref. Di-Mauro2013 , who built the GLF from the radio luminosity function of misaligned AGNs. We consider their best-fitting relation between the $\gamma$-ray and radio luminosities $\mathcal{L}$ -$L_{r,{\rm core}}$ and assume a power-law spectrum with index $\alpha_{\rm mAGN}=2.37$. To derive the GLF of star-forming galaxies, we start from the infrared luminosity function Gruppioni2013 (adding up spiral, starburst, and SF-AGN populations of their Table 8). Then we relate $\gamma$-ray and infrared luminosities using the best-fitting $\mathcal{L}$-$L_{\rm IR}$ relation from ref. Ackermann2012 . The energy spectrum is taken to be a power-law with spectral index $\alpha_{\rm SFG}=2.7$.

D.3 Three-dimensional seed power spectra

To compute the three-dimensional cross-power spectrum $P_{\gamma\delta}$ between the clustering of a given population of $\gamma$-ray emitters and the matter density field, we follow the halo model formalism (e.g. ref. Cooray2002 ), and write $P_{\gamma\delta}=P_{\gamma\delta}^{\rm 1h}+P_{\gamma\delta}^{\rm 2h}$. Below we derive the 1- and 2-halo terms for the various cases (see also ref. Fornengo:2014 ).

D.3.1 Dark matter $\gamma$-ray sources

The 3D cross power spectrum between $\gamma$-ray emission from particle dark matter and matter density is given by:

[TABLE]

where $P^{\rm lin}$ is the linear matter power spectrum, $b_{\rm h}$ is the linear bias (taken from the model of ref. Sheth1999 ), and $\hat{u}_{\delta}(k|M,z)$ is the Fourier transform of the matter halo density profile, i.e. $\rho_{\rm h}({\bm{x}}|M,z)/\bar{\rho}_{\rm DM}$. The Fourier transform of the $\gamma$-ray emission profile from annihilating dark matter is described by $\hat{v}_{\gamma_{\rm DM}}(k|M,z)$, and it is related to the square of the dark matter density profile. For its precise form, see the appendix of ref. Cuoco2015 .

D.3.2 Astrophysical $\gamma$-ray sources

The cross-correlation of the matter density with astrophysical $\gamma$-ray sources is given by the 3D power spectrum:

[TABLE]

where $b_{\rm S}$ is the bias of $\gamma$-ray astrophysical sources with respect to the matter density, for which we adopt $b_{\rm S}(\mathcal{L},z)=b_{\rm h}[M(\mathcal{L},z)]$. That is, a source with luminosity $\mathcal{L}$ has the same bias $b_{\rm h}$ as a halo with mass $M$, with the relation $M(\mathcal{L},z)$ between the mass of the host halo $M$ and the luminosity of the hosted object $\mathcal{L}$ taken from ref. Camera:2014rja . The mean flux $\langle f_{\rm S}\rangle$ is defined as $\langle f_{\rm S}\rangle=\int\mathrm{d}\mathcal{L}\mathrm{d}F/\mathrm{d}E\Phi_{\rm S}$.

Appendix E Blinding and unbliding procedure

The analysis has been performed by adopting a blinding technique, which relied on the execution of the cross-correlation analysis on nine combinations of data – one true and the other artificial. For both $\gamma$-rays and gravitational lensing we adopted 3 data variants, one of which was the true one. After the data were created, they were randomly assigned symbolic names. The analysis was then performed on all 9 combinations of data, without knowing which was the true version.

The aim of the blinding was to ensure that our analysis would not falsely detect a signal. Criteria to test the blind analysis were defined beforehand. After recording the results of the blind analysis and agreeing on its interpretation, the symbolic name assignments were revealed.

The construction of the versions of the data and the unblinding procedure is discussed below.

E.1 Blinded $\gamma$-ray maps

The $\gamma$-ray maps adopted in the blind analysis are:

• •

G0: The true Fermi-LAT $\gamma$-ray maps.

• •

G1: Poissonian random map with a constant expected photon count over all the unmasked pixels, i.e. in each pixel of the map we extract a random number from a Poissonian distribution with fixed mean. The mean was computed by taking the average counts of the real maps in the unmasked pixels in each energy bin, and then multiplying it by a factor of 10 in order to simulate improved statistics. The produced counts maps are then transformed into flux maps by the usual procedure of dividing them by the mean detector exposure in each energy bin and by the pixel area.

• •

G2: Random reshuffle of all unmasked pixels of the real map in each energy bin.

Once produced, the three sets of maps have been blindly and randomly assigned names A, B and C and the association stored.

E.2 Blinded shear maps

The shear signals adopted in the blind analysis are:

• •

K0: The true DES tangential component of galaxy shapes (this is the shear component for which we are looking for the cross-correlation signal when combined with the Fermi maps).

• •

K1: The cross-component (also known as $B$-modes or $\gamma_{\times}$) of galaxy shapes (which instead should yield pure shape noise and a null detection when correlated with the Fermi maps).

• •

K2: A linear combination of the null signal given by K1 and the cross-correlation signal of redMaGiC galaxies at redshift $z=0.2-0.45$ Prat2018 , with the linear combination coefficients chosen such that the signal should be neither plainly visible nor certainly undetectable. Since the correlation between redMaGiC galaxies and gravitational shear is significant, this mock data set is meant to inject in our analysis a situation potentially (but not necessarily) leading to what could be seen as a detection.

Once produced, the three versions of the data have been blindly and randomly assigned names X, Y and Z and the association stored. We note here that the case K1 provides the cross-shear null test investigated in previous attempts of measurement of the signal Shirasaki2014 ; Shirasaki:2016kol ; Troster:2016sgf ; Shirasaki:2018dkz .

E.3 Blind analysis and unblinding procedure

The nine versions of cross-correlations results were all processed and vetted, with no team member aware of which combination represented the true data vector. We proceeded to the unblinding only once a number of (previously defined) criteria where satisfied: essentially, we needed to have at least one set for each $\gamma$-ray map compatible with noise (because of K1) and one with some (possibly weak) signal (because of K2). This has been done by evaluating the $\chi^{2}$ differences defined in the main text as the statistical estimator. In order to evaluate the statistical significance of the obtained $\Delta\chi^{2}$, we derived the distribution of various cases from multivariate Gaussian realizations. In the following, we summarise the main and most relevant results.

We note one substantial change to the analysis after unblinding. Tests with versions of the Fermi-LAT data that differ in the subtraction of a Galactic foreground model revealed a bug in our treatment of weak lensing shear around pixels with negative foreground-subtracted flux. All results in the main text of this paper have been updated, leading to an increase in signal-to-noise ratio from $4.5$ to $5.3$. The error did not affect our unblinding choices, particularly because it did not affect shear measurements around the G1 and G2 maps, which have no pixels with negative flux.

The $\Delta\chi^{2}$ between null signal and models for all the combinations analyzed are reported in Table 3, for both the phenomenological and physical models. By looking at the phenomenological model analysis, the table shows that the majority of cases have a very low $\Delta\chi^{2}$, as expected from the combination of the data sets discussed in the previous Sections. Specifically, for each of the versions of the blinded $\gamma$-ray maps, at least one of the three shear data vectors is consistent with pure shape noise. Two cases exhibit a somewhat large $\Delta\chi^{2}$: AZ and CX. There are two possibilities in the subsets data combination that could provide such a situation: the combination of any version of the $\gamma$-ray map with the redMaGiC galaxy signal or the combination of the true Fermi and the true DES data, if a signal is in fact present (combination G0-K0). These results, as well as those discussed below, remain valid also after the post-unblinding correction mentioned above.

By looking at the physical case results, we notice that the largest $\Delta\chi^{2}$’s occur for CX, while AZ is consistent with no signal. The fact that AZ cannot be fitted by the physical model implies this case does not provide a real signal. Let’s comment that in this part of the analysis, for the physical case we used a reference model (further discussed in the next Section) which assumes for simplicity the normalisation of the 1-halo and 2-halo terms for blazars to be equal ($A_{\rm BLZ}^{1{\rm h}}=A_{\rm BLZ}^{2{\rm h}}$) and which refers to a DM particle annihilating into hadronic states, specifically into a $\bar{b}b$ pair. The null $\chi^{2}$ distribution and the distribution of expected $\Delta\chi^{2}$ for the reference physical model are reported in Fig. 8. They are obtained by drawing from a multivariate Gaussian with mean given by the reference model and covariance given by the same covariance used in the data analysis. The null $\chi^{2}$ distribution peaks around 433. The $\Delta\chi^{2}$ distribution is rather broad and implies that $\Delta\chi^{2}$ of all variants reported in Table 3 are potentially compatible with the true case, since their $\Delta\chi^{2}$ is well inside the distribution in Fig. 8. The peak of the distribution indicates that the expected $\Delta\chi^{2}$ for the reference model is around 28. Fig. 8 shows as shaded areas the $\chi^{2}$ and $\Delta\chi^{2}$ distributions obtained before the post-unblinding correction was applied. No difference is obtained for the $\chi^{2}$, while for the $\Delta\chi^{2}$ a distribution peaked at larger values of $\Delta\chi^{2}$ is obtained after applying the correction, as a consequence of the increased sensitivity to the presence of a signal. The results are consistent with their mutual data sets and show that the bug did not lead to erronoeus conclusions.

The unblinding revealed the following identifications: A $\xrightarrow{}$ G1, B $\xrightarrow{}$ G2, C $\xrightarrow{}$ G0, X $\xrightarrow{}$ K0, Y $\xrightarrow{}$ K1, Z $\xrightarrow{}$ K2. The coefficients of the linear combination in K2 were such that the injected signal was rather small, and therefore not easily identifiable. Moreover, being a signal due to the cross-correlation between the galaxy distribution and gravitational shear, it does not have to be compatible with our physical models (constructed specifically for the cross-correlation with $\gamma$-rays), while instead it could be well described by the phenomenological model (since it contains a generic 1-halo and 2-halo terms).

From all these considerations, we see first of all that the results shown in Table 3 are fully compatible with expectations: there is no “spurious” detection, while the presence of a signal occurs only for those cases for which this is potentially possible. In fact, the CX case (the one corresponding to the combination of the true Fermi and DES data) is the only one that presents a high $\Delta\chi^{2}$ for both the phenomenological and physical model. The $\Delta\chi^{2}$ is larger in the case of the phenomenological model since the latter has more freedom to adapt to data. CY and CZ are compatible with null signal. The smallness of the injected signal in Z makes this case essentially indistinguishable from pure noise. Finally, the null $\chi^{2}$ for CX is 468, compatible with the expectations of the distribution of Fig. 8.

None of the combinations involving A and B maps present a statistically significant $\Delta\chi^{2}$ either with respect to the physical or the phenomenological model, except for the phenomenological fit for the AZ case. We remind that map A is built from Poissonian $\gamma$-ray noise with enhanced photon count statistics. The latter means that the size of the errors is significantly reduced, and now the “small” injected signal is enough to provide a significant detection. On the other hand, the “artificial” nature of the signal is revealed by the fact that the physical model is not able to fit it. The freedom we left to the phenomenological model is instead large enough to make it able to include the signal of cross-correlation with galaxies.

In the next Section we outline the specifications of the reference model used for the blind analysis.

Appendix F Reference model

In the pipeline of the analysis before unblinding, we adopted a physical model with blazars, mAGN, SFG and DM annihilating into a $\bar{b}b$ pair, with a common normalization for the 1-halo and 2-halo terms for blazars, i.e., $A_{\rm BLZ}^{1{\rm h}}=A_{\rm BLZ}^{2{\rm h}}$. This was the first obvious options, since it adopts a nominal model for the blazar terms and minimises the number of free parameters. As seen from the triangle plot in Fig. 9, the parameter scan for the CX case provides upper bounds for the three astrophysical components and a hint for the presence of a DM signal is found at the $3.1\sigma$ C.L., with best-fit parameters $m_{\rm DM}=302$ GeV and normalization of the annihilation cross section relative to the natural scale $A_{\rm DM}=105$. Overall, the statistical significance of this reference physical model is ${\rm SNR}=4.2$, with \Delta\chi^{2}$$=18.2 (for brevity, in this Section, we only quote the significance obtained after the post-unbliding correction mentioned above).

In the right panel of Fig. 9, we show the cross-correlation signal for the best fit of this reference physical model, compared to the data points. While large scales are well fitted by the model, a clear lack of power is present at small angular scales. From Fig. 2 of the main text we instead know that in the phenomenological model, the PSF-like 1-halo term can account for the measured CCF at small scales properly. This implies that the reference model, in its simplest formulation, is excessively limiting the ability to adapt to the data: a larger 1-halo term is needed.

Considering that the data require a hard-spectrum, not compatible with astrophysical components other than blazars, we improve the modelling by allowing the 1-halo and 2-halo terms for blazars to be independently normalised. In this case, the small scale correlation is well fitted, and a peak in the posterior distribution for the 1-halo blazar component arises (see Fig. 3 of the main text). This improves the ${\rm SNR}$ from $4.2$ to $5.2$, while the dark matter component loses some of its significance with respect to the original reference case.

We have also investigated variations of the dark matter model in terms of considering different annihilation channels, as it is usually done in dark matter analyses. We found a preference for leptonic annihilation, with the preferred option being annihilation into $\tau$-lepton pairs: this increases the global ${\rm SNR}$ to $5.2$, with a preference for the presence of a dark matter signal at the $2.8\sigma$ confidence level. This is the model that we report in the main text. In the case the DM component is not included, we obtain $\Delta\chi^{2}=16.5$, as compared to $\Delta\chi^{2}=27$ for the case with DM.

Concerning the blazar 1-halo component, this depends on the relation between the blazar luminosity and the host-halo mass, as can be seen in Eq. 20. In the reference model used for the blind analysis this relation has been taken from ref. Camera:2014rja , where it was derived by linking the $\gamma$-ray luminosity of a blazar to the mass of the supermassive black hole powering the AGN and then relating the mass of the black hole to the mass of the DM halo. Combining these pieces together, this peocedure gives $M(\mathcal{L})=2\times 10^{13}M_{\odot}\,(\mathcal{L}/10^{47}\mathrm{erg\,s^{-1}})^{0.23}(1+z)^{-0.9}$, where $\mathcal{L}$ is the rest-frame luminosity of blazars in the energy range $0.1$ to $100\,\mathrm{GeV}$. Assuming this relation to be a power-law and fixing the redshift dependence as in the reference model (since, as we already mentioned, it is weakly constrained by data, and we therefore do not gain meaningful insight from the fit), we explore in the left panel of Fig. 10 which $M(\mathcal{L})$ relation would be needed in order to reproduce the 1-halo term of our refined model. The plot shows the best-fit relation and its 68% C.L. contours in the plane $(\alpha,M_{0})$ for a relation $M(\mathcal{L})=M_{0}\,(\mathcal{L}/10^{47}\mathrm{erg\,s^{-1}})^{\alpha}(1+z)^{-0.9}$. Our results suggest the average mass of a halo hosting a blazar is larger than the one adopted in Camera:2014rja , and most likely above $10^{14}M_{\odot}$. The cross-correlation signal with weak lensing seems therefore to be dominated by blazars residing in cluster-size halos. In the right panel of Fig. 10, we show with a few examples that modifying the $M(\mathcal{L})$ relation has a dramatic impact on the the 1-halo power spectrum, while it only mildly affects the 2-halo term.