# High-Significance Detection of Correlation Between the Unresolved Gamma-Ray Background and the Large Scale Cosmic Structure

B Thakore, M Negro, M Regis, S Camera, D Gruen, N Fornengo, A Roodman, A Porredon, T Schutt, A Cuoco, et al.

► To cite this version:

B Thakore, M Negro, M Regis, S Camera, D Gruen, et al., High-Significance Detection of Correlation Between the Unresolved Gamma-Ray Background and the Large Scale Cosmic Structure. JCAP, 2025, 06, pp.037. 10.1088/1475-7516/2025/06/037. [hal-04910240]

HAL Id: hal-04910240

https://hal.science/hal-04910240v1

Submitted on 2 Nov 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# High-significance detection of correlation between the unresolved gamma-ray background and the large-scale cosmic structure

To cite this article: B. Thakore *et al* JCAP06(2025)037

View the [article online](#) for updates and enhancements.

## You may also like

- Performance of the CMS high-level trigger during LHC Run 2

A. Hayrapetyan, A. Tumasyan, W. Adam et al.

- Identification of low-momentum muons in the CMS detector using multivariate techniques in proton-proton collisions at sqrt(s) = 13.6 TeV

V. Chekhovsky, A. Hayrapetyan, V. Makarenko et al.

- Muon identification using multivariate techniques in the CMS experiment in proton-proton collisions at sqrt(s) = 13 TeV

A. Hayrapetyan, A. Tumasyan, W. Adam et al.
# High-significance detection of correlation between the unresolved gamma-ray background and the large-scale cosmic structure

**B. Thakore et al.**

Full author list at the end of the paper

E-mail: bhashinashish.thakore@unito.it

**ABSTRACT:** Our understanding of the γ-ray sky has improved dramatically in the past decade, however, the unresolved γ-ray background (UGRB) still has a potential wealth of information about the faintest γ-ray sources pervading the Universe. Statistical cross-correlations with tracers of cosmic structure can indirectly identify the populations that most characterize the γ-ray background. In this study, we analyze the angular correlation between the γ-ray background and the matter distribution in the Universe as traced by gravitational lensing, leveraging more than a decade of observations from the Fermi-Large Area Telescope (LAT) and 3 years of data from the Dark Energy Survey (DES). We detect a correlation at signal-to-noise ratio of 8.9. Most of the statistical significance comes from large scales, demonstrating, for the first time, that a substantial portion of the UGRB aligns with the mass clustering of the Universe as traced by weak lensing. Blazars provide a plausible explanation for this signal, especially if those contributing to the correlation reside in halos of large mass (~$10^{14} M_{\odot}$) and account for approximately 30–40% of the UGRB above 10 GeV. Additionally, we observe a preference for a curved γ-ray energy spectrum, with a log-parabolic shape being favored over a power-law. We also discuss the possibility of modifications to the blazar model and the inclusion of additional γ-ray sources, such as star-forming galaxies, misalting active galactic nuclei, or particle dark matter.

**KEYWORDS:** Bayesian reasoning, gamma ray experiments, weak gravitational lensing, galaxy surveys

ARXIV ePRINT: 2501.10506

<figure><img src="image_2.png" /></figure>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Theoretical framework</td><td>3</td></tr><tr><td>3</td><td>Data</td><td>6</td></tr><tr><td>3.1</td><td>Weak lensing data</td><td>6</td></tr><tr><td>3.2</td><td>Gamma-ray data</td><td>8</td></tr><tr><td>4</td><td>Analysis and results</td><td>9</td></tr><tr><td>5</td><td>Physical interpretation</td><td>14</td></tr><tr><td>6</td><td>Conclusion</td><td>20</td></tr><tr><td>A</td><td>Covariance matrix</td><td>22</td></tr><tr><td>B</td><td>Star forming galaxies</td><td>25</td></tr><tr><td>C</td><td>Misaligned AGNs</td><td>26</td></tr><tr><td>D</td><td>WIMP dark matter</td><td>27</td></tr><tr><td>E</td><td>Data verification and robustness checks</td><td>28</td></tr><tr><td></td><td>Author List</td><td>36</td></tr></table>

## 1 Introduction

High-energy astronomy is an endeavour that can potentially provide insights into the disciplines of astrophysics, cosmology and particle physics. Stemming from extremely violent events in the Universe, γ-rays act as messengers, providing information about the mechanisms of rare events as seen in supernovae, as well as those that occur due to matter under extreme stress, such as in the vicinity of pulsars or Active Galactic Nuclei (AGNs). In addition to their use in astrophysics, γ-ray frequencies can also shed light on an important cosmological component, dark matter (DM), representing approximately 25% of the Universe's energy budget [1]. DM has been theorized to consist of an exotic fundamental particle which may annihilate or decay into standard-model particles, thereby producing cosmic messengers. In the case of Weakly Interacting Massive Particles (WIMPs) or any other potential DM particle with a mass in the GeV range or higher, their annihilation or decay is likely to produce photons in the γ-ray spectrum [2]. However, due to the small cross-section of DM annihilation, the potential number of detectable events is limited. Furthermore, the γ-ray sky is filled with emissions from various astrophysical sources, such as pulsars, supernova remnants, AGN, and cosmic-ray interaction with Galactic interstellar medium and radiation fields. These emissions create background noise, masking potential signals from DM annihilation.

A method for distinguishing between the non-thermal γ-ray emissions originating from astrophysical sources and those that might be caused by DM annihilation or decay within
the Unresolved Gamma-Ray Background (UGRB) hinges on the concept of cross-correlating UGRB maps with various other maps that trace the underlying large-scale structure of the Universe. Such tracers include cosmic phenomena like the weak gravitational lensing effect [3–8], the clustering of galaxies [9–15] and galaxy clusters [16–20], and the lensing effect of the Cosmic Microwave Background (CMB) [21], which reflect the large-scale distribution of matter across cosmological distances (see also refs. [22–25]).

These techniques are especially promising because they serve as direct gravitational probes of the mass distribution in the Universe, the vast majority of which is expected to be DM [1]. By examining the angular, energy, and redshift behaviour of the cross-correlations, it becomes possible to disentangle signals that may arise from ordinary astrophysical sources from those associated with DM interactions. Astrophysical sources are hosted within DM halos, and thus exhibit correlations with other direct tracers of the underlying large-scale cosmic structure such as galaxy clustering or weak gravitational lensing. On the other hand, they are typically much smaller than DM halos (thus showing different angular spectra), and mainly follow the redshift dependence of star formation as opposed to the WIMP signals, which peak at very low redshifts. Astrophysical sources also have smooth, power-law like spectra in energy while WIMPs exhibit a cutoff and a more curved spectrum. These characteristics are handles to tell apart γ-ray emission from astrophysical and dark matter origin.

Beyond its utility in the search for DM, this approach also has broader implications for understanding the population of unresolved γ-ray sources, providing valuable insights into their redshift distribution and clustering properties. This can help to refine our understanding of a variety of γ-ray emitting objects, such as AGN, distant star-forming galaxies (SFGs), and other yet-to-be-resolved populations of astrophysical objects. A key contributor to the UGRB is the blazar population. Blazars are a type of AGN that is estimated to contribute significantly to the UGRB [26], dominating in particular at fluxes just below the Fermi-LAT source detection threshold [27]. Their influence must be therefore carefully accounted for when attempting to identify the composition of the signals in the aforementioned cross-correlations. Blazars exhibit distinct clustering patterns and spectra, which, when accurately modeled, can help distinguish their signals from other γ-ray sources, improving our understanding of their contribution to the UGRB.

In this work, we present a study involving cross-correlations using weak lensing. The weak lensing observable that we consider here is the tangential shear, which is a result of the distortion of background (source) galaxies as a result of the foreground lensing galaxies. This is known as galaxy-galaxy lensing.

Over the years, several observational attempts at disentangling cross-correlation signals for both their dark and visible nature have been carried out involving weak lensing [5–7]. In particular, a cross-correlation study in 2020 by Amazzalorso et al. [8] obtained the first identification of a cross-correlation signal between the UGRB, as seen by the Fermi-LAT 9 year data, and the distribution of mass in the Universe probed by weak gravitational lensing, measured by the DES Y1 datasets [28]. The cross-correlation signal was detected with a significance of 5.3σ, and found a preference for a DM-inclusive model over a purely astrophysical model at 2.8σ. Here, we build on the aforementioned detection using 12 years of γ-ray data from Fermi-LAT [29, 30], and three-year (Y3) weak lensing shear measurements from DES.
The paper is organized as follows. Section 2 describes the phenomenological and physical models used to interpret the cross-correlation signal. Section 3 describes the Fermi-LAT and DES data used for the cross-correlation analysis. Details of the analysis are reported in section 4, which also offers an interpretation of the results in terms of phenomenological models, outlining the main properties of the signal. Inferences in terms of astrophysical sources follow in section 5, followed by a conclusion in section 6. Appendix A provides details on the construction of the covariance used in the data analysis, while additional discussions on SFGs, misaligned AGNs (mAGNs), and DM annihilation can be found in appendices B, C, and D respectively. Finally in appendix E we elaborate upon the robustness checks that were performed to validate the results of the analysis.

## 2 Theoretical framework

In this work, we study the 2-point angular cross-correlation function between the UGRB and gravitational shear. The UGRB is obtained from photon counts measured by the Fermi-LAT in different energy bins, after removing the contribution from resolved sources and Galactic foreground. The gravitational shear is given by the tangential ellipticity of galaxies measured by DES, and trace the mass distribution in the Universe. Here we want to test whether the γ-ray background fluctuations are sourced by the matter distribution.

For a map of γ-rays in the a th energy bin and a shear catalog in the r th redshift bin, the real-space quantity depicting their cross-correlation $\hat{\Xi}^{ar}$ can be theoretically computed from the harmonic-space cross-power spectrum $C_{\ell}^{ar}$ by a Legendre transform [31]:

$$ \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), \qquad (2.1) $$

with $\theta$ being the angular separation on the sky, $P_{\ell}^{(2)}$ the Legendre polynomial of order two, and $W_{\ell}^{a}$ the beam window function. The latter is computed from a Legendre transform of the Fermi-LAT Point-Spread Function (PSF), see appendix II in [32], and accounts for the finite resolution of the detector (we neglect the DES and pixel smoothing since they act on much smaller scales than the Fermi-LAT PSF).

Correlations that occur at physical scales smaller than the Fermi-LAT PSF can be approximated with $C_{\ell}^{ar} = \text{const}$. On the other hand, correlations on very large scales can be well-described by the clustering in the linear regime, with the angular dependence dictated by the linear matter power spectrum [33], that we computed from the transfer function in [34]. In this study, we have adopted a halo model approach to describe both the phenomenological and physical characteristics of the cross-correlations [35]. In the halo-model framework, all mass in the Universe’s large-scale structure is assumed to reside within virialised DM halos. Consequently, the correlation function splits into two distinct contributions: the 1-halo and 2-halo terms (abbreviated as ‘1h’ and ‘2h’ in the succeeding equations). The 1-halo term captures correlations between two points within the same halo, dominating at small angular scales. On the other hand, the 2-halo term accounts for correlations between points in different halos, aligning with the broader matter distribution and prevailing at large scales.
<table><caption>**Table 1.** Gamma-ray energy bins over which the analysis is performed, 68% containment angles $\theta_{\text{cont}}$ of the Fermi-LAT PSF, and photon counts in the unmasked Fermi area in each energy bin along with the average measured intensity.</caption><thead><tr><th></th><th colspan="9">Bin number</th></tr><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th></tr></thead><tbody><tr><th>E<sub>min</sub> [GeV]</th><td>0.631</td><td>1.202</td><td>2.290</td><td>4.786</td><td>9.120</td><td>17.38</td><td>36.31</td><td>69.18</td><td>131.8</td></tr><tr><th>E<sub>max</sub> [GeV]</th><td>1.202</td><td>2.290</td><td>4.786</td><td>9.120</td><td>17.38</td><td>36.31</td><td>69.18</td><td>131.8</td><td>1000</td></tr><tr><th>θ<sub>cont</sub> 68% [deg]</th><td>1.00</td><td>0.58</td><td>0.36</td><td>0.22</td><td>0.15</td><td>0.12</td><td>0.11</td><td>0.10</td><td>0.10</td></tr><tr><th>Photon counts</th><td>360865</td><td>780820</td><td>551998</td><td>221181</td><td>89897</td><td>38277</td><td>11990</td><td>3757</td><td>1619</td></tr><tr><th>⟨I<sub>a</sub>⟩ [10<sup>-7</sup> cm<sup>-2</sup> s<sup>-1</sup> sr<sup>-1</sup>]</th><td>5.69</td><td>2.18</td><td>0.991</td><td>0.375</td><td>0.157</td><td>0.0670</td><td>0.0197</td><td>0.0062</td><td>0.0032</td></tr></tbody></table>

Our first theoretical model is a phenomenological model constructed as a Power-Law (PL) in energy as follows:

$$ \Xi_{\text{PL}}^{ar}(\theta) = \left[ A_1(E_a/E_0)^{-\alpha_1} \left(\frac{1+z_r}{1+z_0}\right)^{\beta_1} \hat{\Xi}_{\text{PSF-like}}^a(\theta) + A_2(E_a/E_0)^{-\alpha_2} \left(\frac{1+z_r}{1+z_0}\right)^{\beta_2} \hat{\Xi}_{\text{2h-like}}^{ar}(\theta) \right] \frac{\Delta E_a}{\langle I_a \rangle} \quad (2.2) $$

where $E_a$ and $z_r$ are the central values of the energy (measured in GeV) and redshift bins, $E_0$ is the pivot energy, chosen as the geometric mean of the energy bin centres at $E_0 = 13.7$ GeV, $\Delta E_a$ is the width of the energy bin and $\langle I_a \rangle$ is the measured photon flux, provided in table 1. Similarly to $E_0$, we have also defined a pivot redshift as a geometric mean of the centres of the four redshift bins at $z_0 = 0.64$.<fnref n="1" /> Note that the term in square bracket is differential in energy. $\hat{\Xi}_{\text{PSF-like}}^a(\theta)$ is the Legendre transform of the Fermi-LAT PSF (i.e., with $C_{\ell}^{ar} = 1$ in eq. (2.1)) and $\hat{\Xi}_{\text{2h-like}}^{ar}(\theta)$ is a generic large-scale contribution, obtained from linear theory. Correlation functions with a hat have flux units, while those without a hat are normalised to the $\gamma$-ray flux, and therefore dimensionless. We make this distinction because the comparison with the measurement is performed with dimensionless Cross-Correlation Function (CCF) as described below, while the comparison with physical models is easier for dimensional quantities. The two normalisations $A_1$ and $A_2$, the spectral indices $\alpha_1$ and $\alpha_2$, and the redshift evolution indices $\beta_1$ and $\beta_2$ are free parameters of the model. For the aforementioned parameters, we have considered a phenomenological model with a free relative amplitude between the PSF-like and large-scale terms that can capture the angular behaviour of the signal. As we will show in the next section, the redshift behaviour is not strongly constrained in our analysis and, for simplicity, we assume a power-law scaling. Gamma-ray sources typically have energy spectra that can be well approximated by a power-law, and so it is assumed in eq. (2.2).

<footnote n="1">¹The pivot energy and redshift values, while determined arbitrarily, are primarily chosen in order to decorrelate the uncertainties in the slope and the offset. By introducing a pivot energy and redshift, we can remove the degeneracy that arises from the interplay between the normalisation and the spectral and (in the case of a log-parabolic model) curvature indices, which makes the fit less stable. A pivot value removes this degeneracy and makes the fit more stable and easier to interpret by giving us cleaner uncertainties and their concomitant best fits.</footnote>
On the other hand, to explore the possibility of a “curved” spectrum, we devise a second phenomenological model, considering a Log-Parabola (LP) spectrum in energy (which also describes several γ-ray sources of the Fermi-LAT 4FGL catalog [29]):

$$ \Xi_{\text{LP}}^{\text{ar}}(\theta) = \left[ A_1(E_a/E_0)^{-\alpha_1+\gamma_1 \log_{10}(E/E_0)} \left(\frac{1+z_r}{1+z_0}\right)^{\beta_1} \hat{\Xi}_{\text{PSF-like}}^{a}(\theta) + A_2(E_a/E_0)^{-\alpha_2+\gamma_2 \log_{10}(E/E_0)} \left(\frac{1+z_r}{1+z_0}\right)^{\beta_2} \hat{\Xi}_{\text{2h-like}}^{ar}(\theta) \right] \frac{\Delta E_a}{\langle I_a \rangle}, \quad (2.3) $$

where $\gamma_1$ and $\gamma_2$ are the spectral indices providing the deviation from a power-law (henceforth referred to as curvature indices).

A physical model of the angular power spectrum of the cross-correlation between γ-ray sources and the gravitational shear can be derived as [23]

$$ C_{\ell}^{\text{ar}} = \int dE \, dz \frac{1}{H(z)} \frac{W_{\text{gamma}}^{a}(E, z) W_{\text{shear}}^{r}(z)}{\chi(z)^2} P_{\gamma\delta} \left[ k = \frac{\ell}{\chi(z)}, z \right], \quad (2.4) $$

where $\chi(z)$ is the comoving distance to redshift $z$, obeying $dz/d\chi = H(z)$ with $H(z)$ the Hubble rate, $W_{\text{gamma}}^{a}(E, z)$ and $W_{\text{shear}}^{r}(z)$ are the so-called window functions, providing the redshift distribution of the signals, and $P_{\gamma\delta}$ is the three-dimensional cross-power spectrum between a given γ-ray population sourcing the UGRB emission and the matter density contrast $\delta$. In the Limber approximation, the physical scale $k$ (the modulus of the physical wavenumber) and the angular multipole $\ell$ are related by $k = \ell/\chi(z)$. All the ingredients entering eq. (2.4) are detailed in ref. [8], with cosmological parameters from ref. [36].

In the main analysis, we will consider blazars (BLZ) to dominate the γ-ray source population. We discuss the addition of SFGs and mAGNs in appendices B and C. The window function is related to the γ-ray luminosity function (GLF) $\phi$, which depends on luminosity, energy and redshift. For BLZ we will follow ref. [27] and, unless otherwise specified, we set all the parameters describing $\phi$ to the best-fit of the angular auto-correlation found in ref. [27], see the BLL 4FGL+C_P fit in their table 2. Specifically, we adopt the following decomposition of the GLF $\Phi(L_\gamma, z, \Gamma) = dN/dL_\gamma dV d\Gamma$ (defined as the number of sources per unit of luminosity $L_\gamma$, co-moving volume $V$ at redshift $z$ and photon spectral index $\Gamma$) in terms of its expression at $z=0$ and a redshift-evolution function $e(L_\gamma, z)$:

$$ \Phi(L_\gamma, z, \Gamma) = \Phi(L_\gamma, 0, \Gamma) \times e(L_\gamma, z), \quad (2.5) $$

where $L_\gamma$ is the rest-frame luminosity in the energy range (0.1–100) GeV. At redshift $z=0$:

$$ \Phi(L_\gamma, 0, \Gamma) = \frac{A}{\ln(10)L_\gamma} \left[ \left(\frac{L_\gamma}{L_0}\right)^{\kappa_1} + \left(\frac{L_\gamma}{L_0}\right)^{\kappa_2} \right]^{-1} \times \exp \left[ -\frac{(\Gamma - \mu_{\text{BLZ}})^2}{2\sigma^2} \right], \quad (2.6) $$

where $A$ is a normalization factor, the indices $\kappa_1$ and $\kappa_2$ govern the evolution of the GLF with the luminosity $L_\gamma$ and the Gaussian term takes into account the distribution of the spectral indices $\Gamma$ around their mean $\mu_{\text{BLZ}}$, with a dispersion $\sigma$.

The redshift behaviour is given by

$$ e(L_\gamma, z) = \left[ \left( \frac{1+z}{1+z_c} \right)^{-p_1} + \left( \frac{1+z}{1+z_c} \right)^{-p_2} \right]^{-1}. \quad (2.7) $$
Then the γ-ray window function can be then written as [8]:

$$W_{\text{gamma}}^{a}(E, z) = \chi(z)^2 \int d\Gamma \int_{L_{\gamma}^{\min}}^{L_{\gamma}^{\max}} dL_{\gamma} \Phi_{\text{S}}(L_{\gamma}, z, \Gamma) \frac{dN}{dE} \times e^{-\tau[E(1+z),z]}, \quad (2.8)$$

where $dN/dE$ is the γ-ray spectrum, taken to be a power-law, $\Phi_S$ is the GLF of the unresolved source population, $\tau$ is the optical depth due to absorption, and $L_{\gamma}^{\max}$ ensures we are considering only unresolved sources. Let us note here that, in the following, when we mention the γ-ray spectrum, we will refer to the differential (in energy) photon spectrum, i.e., the quantity denoted by $dN/dE$ in eq. (2.8).

As already mentioned, we can separate the 1-halo and 2-halo contributions of the power spectrum, so we can write the physical model as:

$$\Xi_{\text{phys}}^{ar}(\theta) \langle I_a \rangle = A_{\text{BLZ}}^{1h} \hat{\Xi}_{\text{BLZ},1h}^{ar}(\theta, \mu_{\text{BLZ}}, p_1) + A_{\text{BLZ}}^{2h} \hat{\Xi}_{\text{BLZ},2h}^{ar}(\theta, \mu_{\text{BLZ}}, p_1). \quad (2.9)$$

The model parameters, that will be constrained through our Markov Chain Monte Carlo (MCMC) scan, are: two free normalizations for the 1-halo and 2-halo terms, $A_{\text{BLZ}}^{1h}$ and $A_{\text{BLZ}}^{2h}$, which effectively describe the normalization $A$ in eq. (2.6) and the normalization of the halo-matter bias (see e.g. ref. [35] for a complete description of halo biases), the spectral index $\mu_{\text{BLZ}}$ in eq. (2.6), governing the energy dependence, and the redshift parameter $p_1$ in eq. (2.7).<fnref n="2" /> All the other parameters of the GLF are set following the best-fit of BL Lacs<fnref n="3" /> from the analysis of γ-ray number counts and angular auto-correlation in [27].

As for the phenomenological model, all terms are computed in different energy and redshift bins, labeled by indices $a$ and $r$, respectively.

# 3 Data

## 3.1 Weak lensing data

The Dark Energy Survey is a six-year observing program that was carried out with a 570 megapixel camera, the Dark Energy Camera (DECam for short), mounted on the Blanco 4m telescope at the Cerro-Tololo Inter-American Observatory in Chile [38]. Spanning a total area of 4143 deg<fnref n="2" /> after masking, the Year 3 survey provides 2.76 times the weak lensing survey area compared to the 1500 deg<fnref n="2" /> for DES Y1 [39, 40]. The observations are carried out in five broadband filters, namely, the grizY, ranging from ~400 nm to ~1060 nm in wavelength. The DES Year 3 (Y3) analysis utilizes the data acquired over the first three years of observations, amounting to a total of 319 nights. In our analysis we make use of the METACALIBRATION shear catalogue, which is fully described in refs. [41, 42], and included in the DES Y3 shape catalogue of ref. [39]. The DES Y3 shape catalog that is used in our analysis is a subset of the objects in the GOLD catalog [43] that pass the METACALIBRATION cuts, and provides calibrated ellipticity measurements that describe the shapes of individual galaxies after correcting for

<footnote n="2">²We also investigated the possibility of using two different energy and redshift parameters for the 1- and 2-halo terms, but we found that this does not add more information.</footnote>

<footnote n="3">³BL Lacs, short for BL Lacertae, are a type of AGN specifically categorized as blazars. Calling a blazar a BL Lac separates it from the other blazar category of Flat Spectrum Radio Quasars (FSRQs). For more information about blazars see e.g. Mukherjee et al. [37].</footnote>
observational biases. It also includes redshift bin assignments for each galaxy, grouping them by estimated distance to enable studies of large-scale structure and galaxy-galaxy lensing. By introducing artificial shear to images and tracking the response of the estimator to the applied shear, one can address both model and noise biases through the introduction of a mean response factor *R*. The methodology indicates, therefore, that *METACALIBRATION* can be applied to calibrate any shear estimator, including the shapes that have been derived from model fitting or weighted moments. *METACALIBRATION* has been shown to be accurate at the part-per-thousand level in the absence of blending with other galaxies [41], and at the part-per-hundred level for the blending present in the DES Y3 data [39].

After accounting for the additive and multiplicative biases, the wide-field data for the weak lensing shape catalogue consists of approximately 100 million galaxies, having an effective source number density of $n_{\text{eff}} = 5.59 \text{ gal/arcmin}^2$, and a corresponding shape noise of $\sigma_e = 0.261$ [39]. The measurements were carried out in the *riz* bands, with the *g*-band excluded due to insufficient PSF modelling.

In order to establish reliable cross-correlation constraints, it is necessary to calibrate the redshift measurements for the source galaxies properly, allowing for a reliable estimate of the tangential shear that becomes one of the two key components of the two-point estimator that generates the cross-correlation values. For DES Y3, the method used to determine and calibrate the photometric redshift distribution of the wide-field galaxies is a combination of two methods: the Self-Organizing Maps $p(z)$ (abbreviated as SOMPZ) [44], and the clustering redshift technique WZ [45].

Some uncertainties arise from factors like the limited coverage of the Deep Fields and the finite number of simulated BALROG sources [46], which inject mock galaxies into real survey images to account for systematic biases (see, for example, ref. [47]). To mitigate these uncertainties, an ensemble of redshift distribution realisations, $n_i(z)$, is constructed for each redshift bin. By generating multiple realisations, the analysis captures a range of possible outcomes, ensuring more robust estimates of the underlying redshift distribution. Additionally, the clustering redshift technique (WZ) improves redshift constraints by cross-correlating the weak lensing sources with galaxies of known redshift.

The computation of the likelihood functions also relies on shear ratios (the full methodology is described in ref. [48]). Shear ratios provide additional constraining and validation power through the measurement of the galaxy-galaxy lensing signal of a lens galaxy redshift bin at small scales. They therefore reflect the ratio of mean lensing efficiencies of objects in those source bins with respect to the lens bin redshift. This, in turn, depends on the redshift distribution of the sources. Such a method is essentially independent from SOMPZ and clustering redshifts, because of its utilization of lensing signals.

The four DES Y3 redshift bins are depicted in figure 1. The mean redshift for each bin is chosen from ref. [49], with $\langle z_1 \rangle = 0.339$, $\langle z_2 \rangle = 0.528$, $\langle z_3 \rangle = 0.752$, and $\langle z_4 \rangle = 0.952$.<fnref n="4" />

The purpose of this work primarily concerns the detection of the cross-correlation between the UGRB and the tangential shear, as opposed to robust constraints on cosmological parameters. Due to the low Signal-to-Noise Ratio (SNR) of the work, there is an uncertainty

<footnote n="4">⁴The DES Y3 component of the analysis uses updated tomographic binning compared to older Y3 analyses. This is implemented due to the changes observed in the $\Delta\chi^2_{\text{min}}$ and cosmology results, as explained in footnote 5 of ref. [50].</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Normalized redshift distribution $n(z)$ of the DES Y3 source galaxies, with the data taken from ref. [36]. Reprinted figure with permission from [36], Copyright (2022) by the American Physical Society.</figcaption></figure>

in the shear bias due to blending (as discussed in ref. [40]) that is not propagated in this paper. The bias due to blending, however, is negligible here and will therefore not affect the interpretations arising from the analysis.

## 3.2 Gamma-ray data

Fermi-LAT is a γ-ray pair-conversion telescope that has been in operation since 2008. With its broad energy range (~30 MeV to > 1TeV) and effective rejection of charged cosmic-ray background, it is well-suited for studying the UGRB. The telescope scans the entire sky every three hours, achieving an angular resolution for gamma rays of about 0.1 degrees above 10 GeV. The γ-ray data used in the cross-correlation measurement rely on a statistic of 12 years of observation of the Fermi-LAT telescope, from August 4, 2008 to August 8, 2020 (Mean Elapsed Time (MET) = 239557417–618050000.0). The all-sky count and exposure maps are obtained with the Fermitools (v2.2.0).<fnref n="5" /> We use Pass8-R3 processed data and select SOURCEVETO_V2 event class<fnref n="6" /> and the combination of PSF1+2+3 event types, i.e., excluding the quartile of events with the worst angular reconstruction (labelled as PSF0) that would impact negatively on the measurement of the angular correlation at small scales. We have chosen SOURCEVETO_V2 as it provides the best compromise in terms of acceptance and background rejection. In fact SOURCEVETO_V2 has an acceptance comparable to P8R2_clean_v6, with a residual contamination almost equal to that of P8R2_ultracleanveto_v6 at all energies.<fnref n="7" /> This results in a clean photon dataset with sufficient event statistics, which is particularly critical when conducting cross-correlation analyses in order to ensure enough events when suppressing background is essential. The model of the PSF as a function of the energy

<footnote n="5">⁵https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/</footnote>

<footnote n="6">⁶See https://www.slac.stanford.edu/exp/glast/groups/canda/lat_Performance.htm.</footnote>

<footnote n="7">⁷See https://www.slac.stanford.edu/exp/glast/groups/canda/archive/pass8r3v2/lat_Performance.htm.</footnote>
is obtained using the gtpsf routine of the Fermitools, using the appropriate Instrument Response Function (IRFs) for the selected events.<fnref n="8" /> For each energy bin we determine an effective average PSF, by weighting the differential PSF(E) by the UGRB spectrum [52]. We generate γ-ray maps in HEALPix [53] format with $N_{side} = 1024$. Such pixeling is optimal for this analysis as it is comparable to the smallest angular scales probed by the DES Y3 weak lensing analyses, and is better than the Fermi-LAT PSF at any energy. As mentioned in section 2, we account for the Fermi-LAT PSF in our modelling.

To get the final flux maps, we follow the same approach as in refs. [8, 24, 32], which we summarize here for convenience. First, we produced counts maps in 100 micro, logarithmically spaced, energy bins between 100 MeV and 1 TeV; then we divide the count maps of each micro bin by the corresponding average-exposure maps computed in the same micro bin, and corrected by the pixel area, to obtain flux maps in units of photons/cm²/s/sr. The resulting flux maps are then summed up to produce maps for 9 macro energy bins between 631 MeV to 1 TeV, as reported in table 1. We discard $E < 0.6$ GeV because at low energy the angular resolution is too poor for our purposes (see following paragraphs).

To extract the UGRB component in the γ-ray maps, we exclude the majority of the Galactic foreground emission and the resolved point sources listed in the 4FGL-DR2 Catalog [29]. While the former is not expected to contribute to the cross-correlation signal with the (extragalactic) gravitational lensing, it nevertheless provides a noise term to our measurement. Therefore, in order to remove the majority of the bright emission from the Galactic plane, we apply a flat mask on latitudes $|b| < 30$ deg. Furthermore, to remove Galactic large-scale emission at higher latitudes, we perform a template fitting of the Galactic diffuse emission model and its subsequent subtraction at the micro-energy-bin level, following the procedure adopted in ref. [32] and utilizing the new Galactic emission model gll_iem_v07.fits.<fnref n="9" /> Finally, in order to extract the unresolved emission we are interested in, resolved point sources are removed from the maps following the masking approach described in ref. [32], which conservatively masks a region around each catalogue source determined by taking into account both the source's brightness and the PSF in the specific energy bin. Each energy bin has therefore a unique mask, and the fractions of sky available at different energies can be found in appendix A.

In figure 2, we show an example of the Fermi-LAT γ-ray intensity map in equatorial coordinates in the (9.120–17.378) GeV energy bin with the application of the mask described above and illustrating the effect of the Galactic foreground subtraction. In the maps we also illustrate the overlapping and the mass distribution of the DES footprint (as detailed in ref. [51]).

# 4 Analysis and results

To measure the cross-correlation between the UGRB and gravitational shear, we compute the following estimator of the 2-point correlation function (see also ref. [54]):

$$ \Xi^{ar}(\theta) = \Xi_{\Delta\theta_h, \Delta E_a, \Delta z_r}^{\text{signal}} - \Xi_{\Delta\theta_h, \Delta E_a, \Delta z_r}^{\text{random}} = \frac{\sum_{i,j} e_{t,ij}^r I_j^a}{R \sum_{i,j} I_j^a} - \frac{\sum_{i,j} e_{t,ij}^r I_{j,\text{random}}^a}{R \sum_{i,j} I_{j,\text{random}}^a}, \quad (4.1) $$

<footnote n="8">⁸https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_LAT_IRFs/index.html</footnote>

<footnote n="9">⁹https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Left: the DES Y3 sky coverage (in light gray) superimposed on the Fermi-LAT γ-ray intensity map (after Galactic foreground masking and subtraction) for photons in the (9.120-17.38) GeV energy range. The plot is in a Mollweide projection with equatorial coordinates, and has been downgraded to $N_{side} = 256$ and smoothed with a Gaussian beam of size $\sigma = 0.4^{\circ}$ for visualization purposes. The Galactic foreground has been masked and subtracted with the procedure described in section 3.2. Right: the map zoomed in on the surveys' overlapping regions, showing the DES Y3 mass map [51] along with the overlapping region in the Fermi-LAT 12 year γ-ray background. The map refers to the same energy bin as the image on the left.</figcaption></figure>

where $\Xi_{\Delta\theta_h, \Delta E_a, \Delta z_r}^{\text{signal}}$ is the correlation function in configuration space of the two observables measured in different angular $(\Delta\theta_h)$, γ-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$ in redshift bin $r$ relative to a pixel $j$, $e_{ij,t}^r$, multiplied by the Fermi-LAT photon intensity flux in the $a$-th energy bin and in pixel $j$, $I_j^a$. The sum runs over all unmasked pixels $j$ and all sources $i$ in the DES shear catalogue, and it is performed in each of the different photon energy bins and source galaxies redshift bins. The mean response $R$ is determined as described above.

From the correlation function, we remove $\Xi_{\Delta\theta_h, \Delta E_a, \Delta z_r}^{\text{random}}$, the measurement of tangential shear around random lines of sight. This is done by setting $I_{j,\text{random}}^a = 1$ anywhere within the sky region used for γ-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, e.g., refs. [54, 55]).

Using the estimator shown in eq. (4.1), we perform the cross-correlation measurement in 12 logarithmically-spaced angular bins, with radii between 5 and 600 arcmin, 9 photon energy bins, detailed in table 1, and the 4 redshift bins introduced in section 3.1, for a total of 432 bins. The cross-correlations have been computed using the `TreeCorr` [56, 57] package, with the MCMC scans performed using the affine-invariant ensemble sampler as configured in the Python package `emcee` [58], and the contour plots obtained using `ChainConsumer` [59].

For the MCMC analysis, the total number of walkers considered are twice the number of the free parameters of the model. We consider a uniform prior for each parameter and a Gaussian likelihood. For the spectral index $\mu_{\text{BLZ}}$, the prior range is [1.5, 3.0], based on findings in ref. [27]. We assume broad, uninformed priors for the remaining physical and phenomenological parameters. The convergence of the chains is assessed by computing the autocorrelation and the concomitant autocorrelation time as described in ref. [60].
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Measurement and models of the cross-correlation between γ-ray photons and gravitational tangential shear, for the log-parabola and power-law models, showing the angular behaviour (top left), along with the energy (top right) and redshift (bottom) dependences. The second angular bin and the second highest energy bins consisted of negative data points when averaged across energy and redshift/angle and redshift. They have thus been displayed with their respective 2-sigma upper limits, followed by a downward arrow. It is important to note here that these are one-dimensional projections of a multi-dimensional fit, shown to provide a visual representation of the models used. See text for details on the derivation of the reported quantities.</figcaption></figure>

For a qualitative assessment of the measured signal, we plot the angular, energy and redshift dependence of the cross-correlation measurements in figure 3, along with their 68% and 95% confidence intervals in figure 4. For the angular dependence, we simply average the cross-correlations and their concomitant phenomenological and physical model counterparts along energy and redshift. The energy and redshift dependencies, on the other hand, have been calculated using a matched filter amplitude $A = \Xi^\mathrm{T} \Gamma^{-1} \bar{\Xi}_\mathrm{M}/(\bar{\Xi}_\mathrm{M}^\mathrm{T} \Gamma^{-1} \bar{\Xi}_\mathrm{M})$,<fnref n="10" /> where $\Gamma$ is the covariance matrix (more information in appendix A), and $\bar{\Xi}_\mathrm{M}$ is given by a simple 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 $\Xi$ is the estimator defined in eq. (4.1), or the best fit phenomenological and physical models. Note that with boldface notation we indicate the full vector of the correlation function, while so far we have indicated the individual components of the vector. The simple model is then given as $\Xi(\theta) = A_s/\theta$,

<footnote n="10">¹⁰More information on the derivation of the matched filter SNR can be found on https://github.com/Bhash inT/Matched_F Filter_FermixDES.git.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Left: constraints on the parameters of the log-parabola phenomenological model. Right: constraints on the parameters of the power-law phenomenological model. In both panels, the 2D contours refer to the 68% and 95% credible regions, with the shaded areas in the 1-D subplots denoting the 68% credible interval for the associated posterior distributions.</figcaption></figure>

where $A_s$ is the normalization amplitude of the simple model. The error on $A$ is given by $\sigma_A^2 = (\Gamma^{-1}\bar{\Xi}_M)^T \Gamma (\Gamma^{-1}\bar{\Xi}_M)/(\bar{\Xi}_M^T \Gamma^{-1}\bar{\Xi}_M)^2$. In order to calculate the evolution of the energy (redshift) bins, we take all the data vectors (i.e., as cross-correlation signals or as their best fit phenomenological or physical models) and their covariances in each redshift (energy) bin, and calculate the amplitudes and standard deviations for the selected data vectors and models. We can clearly see a positive detection with a differential signal scaling in energy with a spectral index around two and some curvature, and a mildly increasing redshift behaviour.

To be more quantitative, and in order to determine the statistical significance of the signal, we test the deviation of the measurement from a null signal (pure noise) by means of the phenomenological model introduced above, and using two statistical methods.

First, we perform a $\Delta\chi^2$ test statistics, with the chi-squared defined as:

$$ \chi^2(P_{\text{mod}}) = [\Xi_{\text{data}} - \Xi_{\text{th}}(P_{\text{mod}})]^T \Gamma^{-1} [\Xi_{\text{data}} - \Xi_{\text{th}}(P_{\text{mod}})], \quad (4.2) $$

where $\Xi_{\text{data}}$ is the data vector, and $\Xi_{\text{th}}$ is the theoretical cross-correlation for the models outlined above, described by the parameter set $P_{\text{mod}}$. The $\Delta\chi^2$ is defined as $\Delta\chi^2_{\text{mod}} = \chi^2_{\text{null}} - \chi^2(P_{\text{mod}})$, with $\chi^2(P_{\text{mod}})$ computed at the model parameter values $P_{\text{mod}}$ that best fit the data, and $\chi^2_{\text{null}}$ referring to no signal, i.e. $\Xi_{\text{th}} = 0$. The best fits and confidence intervals of the parameters are found in an MCMC likelihood analysis. The second estimator of the significance of the signal is the matched filter signal-to-noise ratio (see e.g. ref. [61]),

$$ \text{SNR}(P_{\text{mod}}) = \frac{\Xi_{\text{data}}^T \Gamma^{-1} \Xi_{\text{th}}(P_{\text{mod}})}{\sqrt{\Xi_{\text{th}}^T(P_{\text{mod}}) \Gamma^{-1} \Xi_{\text{th}}(P_{\text{mod}})}}, \quad (4.3) $$

and we will evaluate $\text{SNR}_{\text{mod}} \equiv \text{SNR}(P_{\text{mod}}^*)$.
<table><caption>**Table 2.** $\Delta\chi^2_{\text{mod}}$ and $\text{SNR}_{\text{mod}}$ computed for the phenomenological and physical models, using either the full data set or the various subsamples discussed in the text. For the Low-$z$ case we selected the two first redshift bins ($z \lesssim 0.6$), while for the High-$z$ case the last two bins ($z \gtrsim 0.6$); the Low-$E$ subset refers to the first four energy bins, i.e., energies below 9 GeV, while the High-$E$ to the bins at higher energies; finally, the Small-$\theta$/Large-$\theta$ cases correspond to data points below/above 3 times 68% containment angle of the Fermi-LAT PSF.</caption><thead><tr><th rowspan="2"></th><th colspan="6">Data set</th><th rowspan="2">Model</th></tr><tr><th>Full</th><th>Low-z</th><th>High-z</th><th>Low-E</th><th>High-E</th><th>Small-θ</th><th>Large-θ</th></tr></thead><tbody><tr><td>Δχ<sup>2</sup><sub>lp</sub></td><td>78.9</td><td>3.40</td><td>75.3</td><td>23.3</td><td>55.4</td><td>8.09</td><td>71.26</td><td>Log-parabola</td></tr><tr><td>SNR<sub>lp</sub></td><td>8.89</td><td>2.45</td><td>8.70</td><td>4.89</td><td>7.49</td><td>2.86</td><td>8.45</td><td>Log-parabola</td></tr><tr><td>Δχ<sup>2</sup><sub>pl</sub></td><td>51.4</td><td>0.92</td><td>53.7</td><td>15.8</td><td>38.7</td><td>11.58</td><td>44.10</td><td>Power-law</td></tr><tr><td>SNR<sub>pl</sub></td><td>7.17</td><td>1.66</td><td>7.34</td><td>4.03</td><td>6.24</td><td>3.40</td><td>6.64</td><td>Power-law</td></tr><tr><td>Δχ<sup>2</sup><sub>phys-BLZ</sub></td><td>52.4</td><td>3.19</td><td>53.7</td><td>15.6</td><td>41.3</td><td>10.69</td><td>47.39</td><td>Physical</td></tr><tr><td>SNR<sub>phys-BLZ</sub></td><td>7.23</td><td>1.83</td><td>7.32</td><td>4.04</td><td>6.47</td><td>3.27</td><td>6.88</td><td>Physical</td></tr></tbody></table>

<table><caption>**Table 3.** The 68% credible interval (C.I.) and global best fit values of the parameters for the power-law and log-parabola phenomenological models, denoted as PL and LP respectively.</caption><thead><tr><th>Parameter</th><th>68% C.I. (PL)</th><th>Best fit (PL)</th><th>68% C.I. (LP)</th><th>Best fit (LP)</th></tr></thead><tbody><tr><td>A<sub>1</sub> (×10<sup>-12</sup>)</td><td>[0.2, 19.6]</td><td>17.3</td><td>[0.35, 22.84]</td><td>1.48</td></tr><tr><td>A<sub>2</sub></td><td>[0.058, 0.091]</td><td>0.077</td><td>[0.15, 0.21]</td><td>0.20</td></tr><tr><td>α<sub>1</sub></td><td>[1.82, 2.67]</td><td>2.13</td><td>[2.53, 4.98]</td><td>0.871</td></tr><tr><td>α<sub>2</sub></td><td>[1.93, 2.09]</td><td>2.01</td><td>[1.49, 2.13]</td><td>1.94</td></tr><tr><td>γ<sub>1</sub></td><td>N.A.</td><td>N.A.</td><td>[1.55, 4.80]</td><td>0.073</td></tr><tr><td>γ<sub>2</sub></td><td>N.A.</td><td>N.A.</td><td>[1.22, 2.50]</td><td>1.61</td></tr><tr><td>β<sub>1</sub></td><td>[0.76, 6.91]</td><td>4.63</td><td>[0.40, 7.91]</td><td>5.17</td></tr><tr><td>β<sub>2</sub></td><td>[3.45, 6.20]</td><td>4.83</td><td>[3.46, 6.20]</td><td>4.83</td></tr></tbody></table>

In table 2 we present the results on detection significance. The phenomenological model results for the full data show clear evidence for the presence of a cross-correlation signal at the level of $\text{SNR}_{\text{mod}} = 8.9$ for the log-parabola model, and $\text{SNR}_{\text{mod}} = 7.2$ for the power-law model. For both model types, we display in table 3 the best-fit values of all parameters, obtained from the maximum of the joint posterior distribution, and the 68% credible intervals, obtained from the 1D marginalized posterior distribution. 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 last two redshift bins; Low/High-$E$ bins are defined by being below/above 9 GeV, i.e., first four/last five energy bins; and Small/Large-$\theta$ separates angular scales below/above 3 times the 68% containment angle of the Fermi-LAT PSF, reported in table 1. From table 2 we infer that the signal is mostly concentrated at high energies, large angles, and high redshifts. Higher significance from higher redshift bins is somewhat expected (see also the
physical model below), because those bins have a higher lensing signal by integrating over longer line-of-sight distances. The evidence for correlation at large angles suggests that the measurement is not dominated by few very massive and very bright objects, but rather it comes from a clustered population of extragalactic sources. The low energy bins suffer from the poor Fermi-LAT angular resolution which prevents to obtain a signal at large significance. The evidence at large energies points towards an interpretation in terms of sources with a hard spectrum, namely with a relatively large amount of photons at high energy. While the obtained redshift and angular features are common to any γ-ray source population following large scale structures, the γ-ray spectrum might suggest a preference for blazars, rather than softer sources, such as SFGs and mAGNs. The best-fit for the spectral index of the 2-halo component for the log-parabolic model, $\alpha_2 = 1.94_{-0.45}^{+0.19}$ as well as the spectral index component for the power-law model $\alpha_2 = 2.01_{-0.08}^{+0.08}$ is quite hard with respect to the spectral index of the average intensity of the UGRB, which is $\alpha \approx 2.3$ [52], 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.

The spectral scaling has a noticeable curvature, with the log-parabola model having a significantly higher SNR compared to its power-law counterpart. The log-parabola model is strongly favored, at a $\Delta\chi^2 \sim 27$. This could point towards multiple physical phenomena, that will be discussed in the next section.

# 5 Physical interpretation

Having clearly assessed the presence of the cross-correlation signal, in this section, we attempt to disentangle the astrophysical γ-ray components of the cross-correlation signal. Different emitters can lead to different shapes of the cross-correlation signal as a function of the angular separation, energy, and/or redshift. In this work, we consider three astrophysical components, namely, BLZs, SFGs and mAGNs, and discuss also the possibility of a particle DM contribution.

We anticipate that when confronting data with expectations from astrophysical populations, we find that BLZ are preferred by the fit, with negligible contributions from SFG and mAGNs, as found also in [8]. This is essentially due to the harder γ-ray spectrum of the signal, compatible with BLZ but not with SFG and mAGNs, as already discussed above. We therefore start with the most economical choice by including only BLZ in our physical model, and come back to SFGs, mAGNs, and DM towards the end of this section and in the appendix.

For the purpose of this work, blazars can be considered as point-like sources — i.e. their size is much smaller than the Fermi-LAT PSF. Additionally, the size of the halo hosting blazars rarely exceeds the Fermi-LAT PSF. As a consequence, the angular correlation function for the 1-halo term essentially follows from the detector PSF. However, as seen from table 2 and figure 5, a cross-correlation on larger angular scales is required by the fit, and the statistical significance of the signal is actually driven by the 2-halo component. This means that we are able to probe the clustering of blazars, i.e., their large scale distribution. In order to investigate the BLZ properties needed to account for the measured cross-correlation signal, we perform the statistical tests discussed in the previous section, but now with a physical model, based on a characterisation of the 1- and 2-halo cross-correlations of weak lensing
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Similar to figure 3 but for the physical model with BLZ-only contribution. We display the angular (top left), energy (top right), and redshift dependence (bottom) of the cross-correlation.</figcaption></figure>

with γ-rays from blazars, as described in eq. (2.9). In ref. [27] it was found that two BLZ populations are needed to explain the angular power spectrum (APS) data, with FSRQs dominating at GeV energies and BL Lacs taking over from a few GeV. As discussed in the previous section, we have inconclusive results and little constraining power in the first two energy bins. Therefore we focus our modeling on BL Lacs.<fnref n="11" />

Note that we allow the 1-halo and the 2-halo terms to be separately adjusted in the fit against the data. The two normalization parameters are a simple effective way to account for a different luminosity dependence of the GLF (leading to a different amplitude of the signal) and for a different function for the host-halo mass versus γ-ray luminosity (which impacts differently the linear and non-linear bias, thus altering the ratio between 1- and 2-halo terms).

The results are shown in table 2 (along with the best fit parameter values shown in table 4 in a manner similar to that of table 3), where the overall significance of the presence of a signal, and its scalings as a function of energy, angular scales and redshift are all confirmed. Actually, the BLZ case provides results that are very similar to the phenomenological power-law model.

More details of the parameter constraints are shown in figure 6, where the triangular plot of the posterior distributions of the model parameters are reported. The posterior exhibits

<footnote n="11">¹¹The UGRB also contains unresolved Gamma-Ray Burst (GRB) components; however, it is highly unlikely that they are relevant in the energy bins below 100 GeV, i.e., the contributions from GRBs can only notably come from the final energy bin, and occupy a very small fraction of the UGRB (see e.g. refs. [62–64]). Keeping in mind that we do not expect significant UGRB contribution from GRB in the energy range we are looking at with integrated time of observation of 16 years, we can reliably conclude that the current GRB models cannot explain the spectral curvature that we observe in our work, and therefore have not been considered.</footnote>
<table><caption>**Table 4.** The 68% credible interval and global best fit values of the physical model parameters of eq. (2.9).</caption><thead><tr><td>Parameter</td><td>68% C.I.</td><td>Best-fit</td></tr></thead><tbody><tr><td>A<sub>1halo</sub></td><td>[7.51, 49.22]</td><td>34.17</td></tr><tr><td>A<sub>2halo</sub></td><td>[4.36, 6.70]</td><td>6.59</td></tr><tr><td>μ<sub>BLZ</sub></td><td>[1.91, 2.22]</td><td>2.07</td></tr><tr><td>p<sub>1</sub></td><td>[1.02, 5.46]</td><td>1.02</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Constraints on the parameters describing the physical model BLZ, based on the BL Lacs model in [27], and assuming a single population accounting for the measured cross-correlation. As before, the 2D contours refer to the 68% and 95% credible regions, with the shaded areas in the 1D subplots denoting the 68% credible interval for the 1D posterior distribution.</figcaption></figure>

a preference for a large 2-halo term of blazars with normalisation $A_{\text{BLZ}}^{2\text{h}} = 6.59_{-2.23}^{+0.11}$, while the normalisation of the blazar 1-halo term $A_{\text{BLZ}}^{1\text{h}} = 34.17_{-15.05}^{+26.66}$ shows somewhat weaker constraints. The evidence for a 2-halo term is much greater than what was observed in ref. [8]. This nicely meets expectations since the improvement of the DES Y3 dataset with respect to DES Y1 is the larger portion of the sky covered. Note that $A_{\text{BLZ}}^{2\text{h}}$ can be seen as the amplification of the signal due to a larger linear bias and/or GLF with respect to our reference choice. With both being $\mathcal{O}(1)$, one can account for our findings. Given the expected size of systematic uncertainties (see discussion below on $M(\mathcal{L})$ and the bias), we consider it as a good indication that our reference model is able to describe the physical picture.

The spectral index $\mu_{\text{BLZ}} = 2.07_{-0.16}^{+0.15}$ is consistent with BLZ emissions, as already mentioned when discussing the power-law phenomenological model from section 4. The redshift behaviour is also compatible with the BLZ model, see figure 5. The larger significance of high-$z$ bins with respect to ref. [8] might be again attributed to the enhanced sensitivity to large scale structures of the current dataset, which means that the signal is no longer dominated by a few bright (and so typically closer) sources. We tested this conclusion by measuring the correlation of Fermi-LAT maps with the 4FGL sources included, i.e., without masking the extragalactic sources. In this case, due to the presence of bright, low-redshift
sources, we observe a much more dramatic increase of the signal in the low-$z$ bin as compared to the higher-$z$ bins, which confirms the aforementioned explanation of the redshift behaviour. There are only two properties that do not exactly align with our expectations of blazar behaviour. The first has been already emphasized in the previous section, namely, that a curved $\gamma$-ray spectrum is preferred over a power-law (note that the absorbed blazar spectrum has a curvature, which however occurs at high energies, see discussion below). Second, the amplitude of the normalization parameters is surprisingly large, especially for the 2-halo case. To be compatible with the findings in ref. [27], concerning, e.g., the total contribution of BLZ to the UGRB intensity, we would need $A_{\text{BLZ}}^{2h} \lesssim 2$, which is about $3\sigma$ away from our results. Let us now try to address these two points.

A different $\gamma$-ray population might have a curved spectrum. The emission from SFGs has a significant contribution from $\pi_0$ decay, leading to a log-parabolic shape for the $\gamma$-ray spectrum, see, e.g., ref. [65]. On the other hand, the peak occurs below a few GeV, whilst, to fit our data, one needs a peak around 10 GeV, see figure 3. Therefore, the attempt to provide an explanation using SFGs as the log-parabolic component remains unsuccessful. More details are provided in appendix B.

A more appropriate $\gamma$-ray spectrum can be provided by particle DM in terms of WIMPs. Indeed, by choosing the appropriate DM mass and annihilation channel, a significant preference over the BLZ-only model can be obtained. We discuss this case in appendix D.

A less exotic explanation is to fine-tune the BLZ model. One could obtain a log-parabolic spectrum as a result of the combination of different BLZ sub-populations or by directly implementing a curved behavior. On the other hand, it is then not trivial to satisfy bounds from other probes, like the 1-point statistics or the auto-correlation, see, e.g., [27]. A comprehensive test of this hypothesis is beyond the scope of this work.

As a side note, we also emphasize that we detect the presence of a signal at energies above 100 GeV.<fnref n="12" /> This might suggest the presence of a very hard source population on top of the main BLZ contribution (although, currently, there is no statistically significant preference). For example, the magnetohydrodynamical simulations performed in [67] suggest that the majority of the UGRB at these energies originate from cosmic-rays in galaxy clusters, having a very hard spectral index.

Finally, instead of a curvature in the original $\gamma$-ray spectrum, a stronger absorption in the energy range between ~10 and 50 GeV can generate the required behaviour. Gamma-ray absorption is due to pair production from the interaction of $\gamma$-rays with optical and infrared light. To describe absorption, one has thus to estimate the extragalactic background light (EBL). This is not a trivial task, and requires to account for important systematic effects, e.g., the contamination from large zodiacal light associated with interplanetary dust in the Solar system (see e.g. ref. [68] for more information on EBL measurements). The curvature seen in the $\gamma$-ray spectrum of our correlation could be associated to a higher level of EBL at ultra-violet frequencies (that are the most relevant ones to generate pair production when interacting with 10–50 GeV photons) with respect to the reference case we are using [69, 70].

<footnote n="12">¹²To consider the energy dependence in the plots as the physical $\gamma$-ray spectrum is a bit misleading. Since at high energies Fermi-LAT is less sensitive to source detection, the masking becomes less effective (and in order to alleviate for it we added the 3FHL catalog [66], which is specifically devised to the detection of high-energy sources).</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Left: constraints on the parameters of the BLZ model, as in figure 6, but replacing the absorption optical depth of [69] with the one of [71]. The contour plot shown in figure 6 is also displayed below for comparison. Right: the energy behaviour of the 1- and 2-halo components, to show that the alternative absorption model can provide a slightly more curved spectrum in the 10–50 GeV range, compared to the case of figure 5.</figcaption></figure>

To this end, we implement the optical depth obtained from the EBL evaluation in ref. [71]. The latter differs by a factor of a few from the EBL in refs. [69, 70] in the far-UV.

Results are shown in figure 7. The spectrum is indeed slightly more curved than the case in figure 5. The SNR grows to 7.8 with a $\Delta\chi^2 \simeq 8$ with respect to the reference absorption model (compared to the previous SNR of 7.23). Therefore, a larger amplitude of the EBL in the UV can help in recovering the curvature we measured, as expected, even though additional modifications seem to be needed.

We now discuss the possible origin for the large value of $A_{\text{BLZ}}^{2\text{h}}$. The blazar-shear cross-correlation depends on the relation between the blazar $\gamma$-ray luminosity and the host-halo
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Constraints on the parameters describing the relation between the host halo mass and blazar luminosity, $M(\mathcal{L}) = M_0(L_{\text{BLZ}}/10^{47} \text{ erg/s})^{\alpha}(1+z)^{0.9}$. The amplitude $A_{\text{BLZ}}$ is in log-scale (in contrast to the linear scale in the previous plots) to illustrate better its relation with $M_0$.</figcaption></figure>

mass, a quantity which is rather uncertain. We followed ref. [4], 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 black hole mass to the mass of the DM halo: $M(\mathcal{L}) = 2 \times 10^{13} M_{\odot} [\mathcal{L}/(10^{47} \text{ 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 GeV. By considering larger values for the DM halo mass, one can obtain both larger 1-halo term and larger linear bias (i.e., larger 2-halo term), therefore reducing the amplitude needed for $A_{\text{BLZ}}^{\text{1h}}$ and $A_{\text{BLZ}}^{\text{2h}}$.

In order to test this scenario we allow for a more generic $M(\mathcal{L})$ relation: $M(\mathcal{L}) = M_0(L_{\text{BLZ}}/10^{47} \text{ erg/s})^{\alpha}(1+z)^{-0.9}$, with $M_0$ and $\alpha$ being free parameters. For this test, we set the energy and redshift parameters of the BLZ model in eq. (2.9) to the best-fit values found from the above analysis, and set a common normalization to both the 1- and 2-halo terms, $A_{\text{BLZ}}^{\text{1h}} = A_{\text{BLZ}}^{\text{2h}} = A_{\text{BLZ}}$. Along with the best fit value of $\alpha = 0.44^{+0.30}_{-0.34}$, we find the best-fit halo mass parameter value to be $M_0 = 10^{14.1^{+0.98}_{-0.61}} M_{\odot}$, suggesting that the average mass of the halo hosting the unresolved blazars we are probing with this analysis is larger than the average one derived in ref. [4], tailored to galaxy-size halos. Our result agrees with results in ref. [8] where $M_0 > 10^{14} M_{\odot}$ was found.

We can note from figure 8 that if the average halo mass is above $10^{14} M_{\odot}$, then a value of $A_{\text{BLZ}} \simeq 2$ is compatible with data. We have also found this to be true for the alternative BLZ model that we considered, with the model having a best-fit halo mass $M_0 = 10^{14.08^{+0.34}_{-0.81}} M_{\odot}$. Following ref. [27], one can see that $A_{\text{BLZ}} \lesssim 2$ is required in order to be consistent with number counts and auto-correlation analyses and makes the BLZ model providing about 30–40% of the UGRB above 10 GeV. Therefore, if the BLZ-only model is the correct interpretation of our measurement, the weak lensing signal has to be provided by cluster-size halos (i.e., with mass $\gtrsim 10^{14} M_{\odot}$) where the $\gamma$-ray blazars should reside.

The angular correlation between clusters of galaxies and the UGRB has been detected in refs. [16, 19] and the origin of the signal is compatible with $\gamma$-ray emission from blazars
hosted by the clusters. Relatively massive hosts are also found for the blazars above detection threshold, i.e., the ones in the Fermi-LAT catalogue [72].

## 6 Conclusion

In the present work, we measured and interpreted the angular correlation between the γ-ray sky and the matter distribution in the Universe. We employed 12 years of observation of the Fermi-LAT telescope and 3 years of data from DES. Building upon the results obtained with Fermi-LAT 9 yr and DES Y1 in ref. [8], where a signal of cross-correlation between γ-ray sources and gravitational lensing shear was found with an SNR of 5.3, we improved the statistical significance and found an SNR of 8.9.

Thanks to the larger portion of the sky available from DES Y3, we detect for the first time a large-scale correlation, i.e., a proof that a significant fraction of the UGRB follows the mass clustering of the Universe traced by weak lensing.

Blazars are a sufficiently good explanation of the signal, provided that the ones responsible for the measured correlation reside in halos of large mass (~$10^{14} M_{\odot}$), and that BLZ account for about 30–40% of the UGRB above 10 GeV.

We found a significant preference for a curved γ-ray spectrum, with the log-parabolic shape preferred over the power-law at $\Delta\chi^2 \sim 27$. This can have different physical interpretations, including a modification of the BLZ model with respect to the current state of the art, additional γ-ray sources, like SFG and DM, discussed in appendices B and D, or considering different UV EBL models to the default model implemented in this text.

In order to ascertain the properties of the γ-ray source population providing the reported signal, we plan to perform an analysis adding the cross-correlation of the γ-ray sky with the DES Y3 galaxy clustering. For this, we expect an even larger statistical significance with respect to the one reported here for the case of the lensing shear. Thus we should be able to work with sub-samples of galaxies that can emphasize the correlation with a given population, e.g., luminous red galaxies are more suited to explore the DM interpretation since they are not supposed to be SFGs or to host blazars, and to better characterize the redshift dependence of the signal to assess the impact of absorption.

## Acknowledgments

BT, MR, and NF acknowledge support from the Research grant TAsP (Theoretical Astroparticle Physics) funded by INFN. The work of MR and NF is supported by the Italian Ministry of University and Research (MUR) via the PRIN 2022 Project No. 20228WHTYC – CUP: D53C24003550006. SC acknowledges support from PRIN 2022 ‘EXSKALIBUR – Euclid-Cross-SKA: Likelihood InferenceBuilding for Universe’s Research’, Grant No. 20222BBYB9, CUP D53D2300252 0006, from the Italian Ministry of Foreign Affairs and International Cooperation (maeci), Grant No. ZA23GR03, and from the European Union – Next Generation EU. This research was also funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. 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, NSF NOIRLab, 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 NSF Cerro Tololo Inter-American Observatory at NSF NOIRLab (NOIRLab Prop. ID 2012B-0001; PI: J. Frieman), which is managed 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 MICINN under grants PID2021-123012, PID2021-128989 PID2022-141079, SEV-2016-0588, CEX2020-001058-M and CEX2020-001007-S, some of which include ERDF funds from the European Union. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya.

We acknowledge support from the Brazilian Instituto Nacional de Ciência e Tecnologia (INCT) do e-Universo (CNPq grant 465376/2014-2).

This document was prepared by the DES collaboration using the resources of the Fermi National Accelerator Laboratory (Fermilab), a U.S. Department of Energy, Office of Science, Office of High Energy Physics HEP User Facility. Fermilab is managed by Fermi Forward Discovery Group, LLC, acting under Contract No. 89243024CSC000002.

The Fermi LAT collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National
de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515.

**Author contributions:** B. Thakore performed the measurement with D. Gruen, and analyzed the comparison to models with M. Regis. M. Negro provided the Fermi-LAT gamma-ray data. They wrote the manuscript for this paper with constructive feedback throughout the process from the other core authors, namely, S. Camera, N. Fornengo and A. Roodman. The document has been through internal reviews within the DES collaboration, with A. Porredon and T. Schutt acting as internal reviewers, and within the Fermi-LAT collaboration with A. Cuoco acting as the internal reviewer. All of the reviewers have provided valuable expertise and notable feedback that contributed to an improvement in the quality of the paper. The remaining authors have made contributions to the DES Y3 Key Project analysis pipeline, including but not limited to the DES instruments, data collection, processing and calibration, and various analysis pipelines.

## A Covariance matrix

The covariance matrix used for the present analysis is obtained from a combination of the shape-noise term generated from random galaxy rotations, and the theoretically estimated Gaussian large-scale structure covariance. We expect the former term to dominate the covariance, with the latter expected to have a smaller contribution. As in the main text, the redshift bins for the shear are indicated with indices r, s while energy bins for the UGRB are denoted by indices a, b.

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

$$ \hat{\Gamma}_{ar\ell,bs\ell'} = \frac{\delta_{\ell\ell'}^K}{(2\ell + 1)\Delta\ell f_{\text{sky}}} [C_{\ell}^{ar} C_{\ell'}^{bs} + (C_{\ell}^{rs} + \mathcal{N}^{rs})(C_{\ell}^{ab} + \mathcal{N}^{ab})], \quad (A.1) $$

where the $C_{\ell}$'s denote (auto and cross) angular power spectra and $\mathcal{N}$ the noises. $\Delta\ell$ refers to the $\ell$-bin widths and $f_{\text{sky}}$ the sky coverage fraction. Eq. (A.1) can therefore be divided into two main parts: the large-scale structure component, taking the form $C_{\ell}^{ar}C_{\ell}^{bs} + C_{\ell}^{rs}(C_{\ell}^{ab} + \mathcal{N}^{ab})$, and the shape noise terms $\mathcal{N}^{rs}(C_{\ell}^{ab} + \mathcal{N}^{ab})$. All theoretical components involving $\gamma$-rays have been corrected for the Fermi-LAT PSF beam function prior to their inclusion in the covariance calculation. The noise terms $\mathcal{N}$ have no angular dependence. $f_{\text{sky}}$ accounts for the partial coverage of the sky, with $f_{\text{sky}}^{\text{DES}} = 0.11$ for DES Y3 (independent of the redshift bin), and $f_{\text{sky}}^{\text{Fermi}} = (0.073, 0.296, 0.433, 0.519, 0.535, 0.546, 0.554, 0.554, 0.554)$ denoting the $f_{\text{sky}}$ values for the Fermi sky. In order to be conservative, for the cross-correlations analysis we chose the overlap of the DES and Fermi-LAT sky, i.e., $f_{\text{sky}} = f_{\text{sky}}^{\text{DES}} \times f_{\text{sky}}^{\text{Fermi}}$ considering the smallest $f_{\text{sky}}$ value in cases where the energy bins (i.e., a and b) were different.
As mentioned previously, in a Gaussian approximation, the shape noise component can be estimated as $N^{rs}(C_{\ell}^{ab} + N^{ab})$. Expected to be the dominant noise contribution in the covariance, it is important to model this term reliably. To that end, i.e., to go beyond the Gaussian approximations, we generated 2000 realisations of the noise directly from the data in real space. This was done by rotating the ellipticities using an independent random angle between 0 and $\pi$ and calculating the cross-correlations using the ellipticities and the Fermi-LAT maps. The signal generated from this represents a random realization of the shape noise [39, 54, 73, 74].

With a major part of the covariance contribution addressed through the shape noise simulations, we complete the covariance calculations by adding the theoretical estimates of the large-scale structure covariance. The contribution of this term is subdominant in the covariance budget, which is verified both a posteriori, as well as in the previous iteration of this experiment of ref. [8]. This part of the covariance depends on both the cross-correlation signals the $\gamma$-ray-$\gamma$-ray and shear-shear autocorrelations $C_{\ell}^{ab}$ and $C_{\ell}^{rs}$. The details on how to measure the $\gamma$-ray autocorrelation have been described in ref. [32], with $C_{\ell}^{ab}$ being fitted with a simple model, given by a power-law added to a constant, written as $C_{\ell,\text{mod}}^{ab} = A_{ab}\ell^{-\alpha_{ab}} + C_{\text{P}}^{ab}$. The $\gamma$-ray noise term, $N^{ab}$, on the other hand, is computed using eq. (5) of ref. [24]. The shear autocorrelation is calculated using the galaxy redshift distributions described in section 3.1, assuming a $\Lambda$CDM cosmology, with the fiducial values of the model parameters taken from [36].

The large-scale structure part of the covariance is then added to the shape noise term through the creation of a set of 2000 simulated noise realisations from a multivariate Gaussian distribution. These simulations are in physical space and are meant to encode the fluctuations in the large-scale structure. They are added to the shape noise realisations, thereby obtaining 2000 samples, each one with different large-scale and shape noise realisations. These 2000 samples therefore provide a robust estimation of the full covariance matrix. The inverse of the covariance matrix thus obtained is a biased estimator of the inverse of the true covariance, with the bias depending on the number of realisations $N_s$ and bins $N_b$ of the data. The inverse covariance is corrected for bias by multiplying it with the Anderson-Hartlap factor, given as $(N_s - N_b - 2)/(N_s - 1)$ (see ref. [75] for more information, and section VI of ref. [76] for an example of its use in nondiagonal covariances). Using the simulated covariance, the cross-shear measurements returned a chi-squared value of 430, indicating a very close ($< 1\sigma$) agreement with the expected value of 432 for the noise component that is represented by the cross-shear.

In order to establish the reliability of a covariance obtained through simulations, we test it against another covariance obtained through the jackknife method, with each Fermi-LAT flux map divided into 100 jackknife patches. Since the Anderson-Hartlap factor in this case would become negative in the case of 432 bins, we selected an energy-redshift bin combination with a significant signal, and tested the noise modelling for this bin against that of the simulated covariance. For this purpose, we chose the highest redshift and energy bins, and generated 2000 mock datasets having a mean of zero from a multivariate Gaussian distribution using the jackknife covariance. The left side of figure 9 shows the results, depicting roughly the same chi-squared distribution, with the standard deviations of the simulated and jackknife covariances being approximately equal. This is also shown in a qualitative manner by the plot on the right, where we have compared the standard deviations for the jackknife and simulated covariances in the measured signal for the aforementioned bin combination. We have also
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Left: the chi-squared distribution for 2000 mock datasets in the highest energy and redshift bins as a comparison between the jackknife and simulated covariance matrices. Right: the measured cross-correlation signal for the highest energy and redshift bins, comparing the errors between the jackknife and simulated covariances.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Comparing the shape noise standard deviations generated by simulations against the shape noise standard deviations obtained using the Gaussian approximation, for the lowest (0.631–1.202 GeV, left) and the highest (131.825–1000 GeV, right) energy bins, taking the first redshift bin with a centre of 0.343 for the former and the final redshift bin with a centre of 0.964 for the latter.</figcaption></figure>

compared the shape noise obtained by the simulated covariance and the Gaussian covariance, as shown in figure 10 for the lowest energy and redshift as well as the highest energy and redshift bins. This is in contrast to the covariances that were compared in figure 9, where the components were not split into their contributions but were tested as a whole. Here, instead, we solely test the robustness of the shape noise obtained using the simulated covariance by using a Gaussian covariance as a reference. In order to depict the behaviour at larger angular scales in a clearer fashion, we have “flattened” the signal by multiplying the cross-correlation $\Xi(\theta)$ with its corresponding angular bin $\theta$. We find that the simulated covariance has a larger shape noise than its Gaussian counterpart, as expected, due to complex masking effects present in the data [77]. By quoting the shape noise term obtained directly from the data, we can ensure that the goodness of fit tests as well as the errors are reliable.
<figure><img src="image_13.png" /><figcaption>Figure 11. Constraints on the normalization parameters and spectral indices in a model containing blazar and SFGs on left, alongside a description of the angular evolution of the cross-correlations and the best-fit physical models on the right.</figcaption></figure>

To summarize, we have tested the simulated covariance in two different ways: first, by comparing it to an internally estimated covariance, we ensured that it was able to model the noise accurately using the null chi-squared test. Then we considered the dominant noise component of the simulated matrix, i.e., the shape noise, and tested it against its Gaussian counterpart, finding typically larger error bars for the former owing to masking effects present in the data from which the simulations were generated, which was expected. This allowed us to generate the shape noise errors from the data and ensure that the goodness of fit tests were reliable while giving us a full covariance matrix where the Anderson-Hartlap factor could be applied to de-bias the covariance. This validates the robustness of the simulated covariance, which has been used for the analysis in this work.

## B Star forming galaxies

The model of SFG is based on the GLF described in ref. [8]. As a first case, we add the contribution to the cross-correlation signal of such reference SFG model with a free normalization. The physical model of eq. (2.9) becomes:

$$ \Xi_{\text{phys}}^{ar}(\theta) \langle I_a \rangle = A_{\text{BLZ}}^{1h} \hat{\Xi}_{\text{BLZ},1h}^{ar}(\theta, \mu_{\text{BLZ}}, p_1) + A_{\text{BLZ}}^{2h} \hat{\Xi}_{\text{BLZ},2h}^{ar}(\theta, \mu_{\text{BLZ}}, p_1) + A_{\text{SFG}} \hat{\Xi}_{\text{SFG}}^{ar}(\theta), \quad (B.1) $$

where the parameter $A_{\text{SFG}}$ sets the amplitude of the SFG correlation. The constraints on the parameters are shown in the left panel of figure 11. The SNR does not vary with respect to the BLZ-only case reported in table 2, i.e., SNR = 7.2, and the contribution from SFGs is negligible, see right panel of figure 11. This was expected since the γ-ray spectrum of this reference SFG model is a power-law with spectral index of 2.7, something very different from the properties of the measurement described with the phenomenological model.

We then allowed more flexibility to the SFG model, including a log-parabolic spectrum and splitting the 1-halo and 2-halo components, so that eq. (B.1) can be re-written as:

$$ \Xi_{\text{phys}}^{ar}(\theta) \langle I_a \rangle = A_{\text{BLZ}}^{1h} \hat{\Xi}_{\text{BLZ},1h}^{ar}(\theta, \mu_{\text{BLZ}}, p_1) + A_{\text{BLZ}}^{2h} \hat{\Xi}_{\text{BLZ},2h}^{ar}(\theta, \mu_{\text{BLZ}}, p_1) \\ + A_{\text{SFG}}^{1h} \hat{\Xi}_{\text{SFG},1h}^{ar}(\theta, \mu_{\text{SFG}}, \mu_E) + A_{\text{SFG}}^{2h} \hat{\Xi}_{\text{SFG},2h}^{ar}(\theta, \mu_{\text{SFG}}, \mu_E), \quad (B.2) $$
<figure><img src="image_14.png" /><figcaption>Figure 12. Constraints on the normalization parameters and spectral indices in a model containing blazar and mAGNs on left, alongside a description of the angular evolution of the cross-correlations and the best-fit physical models on the right.</figcaption></figure>

where $A_{\text{SFG}}^{\text{1h}}$ and $A_{\text{SFG}}^{\text{2h}}$ are the normalizations on the 1-halo and 2-halo SFG components, while $\mu_{\text{SFG}}$ and $\mu_{\text{E}}$ refer to the spectral and curvature indices of the log-parabola in energy. Due to cosmic ray activity in dense interstellar matter, SFG can exhibit a log-parabolic behaviour in energy owing to $\pi_0$ decay. This peaks around the GeV (see e.g. ref. [65]) and so we set the pivot energy of the log-parabola to 1 GeV, which, together with the prior on $\mu_{\text{SFG}} \subset [1.5, 3]$, ensures the turnaround of the spectrum to be below a few GeV.

Performing the MCMC scan with the above model, we obtained similar results as with the simple power-law case introduced at the beginning of the section. The SNR remained unchanged at 7.2, the presence of SFGs comes out to be subdominant, and the SFG parameters are unconstrained. To obtain some preference for an SFG contribution, we should move the turnaround energy of the curved spectrum to above 10 GeV, something that lacks physical motivation.

## C Misaligned AGNs

For the mAGN model, we follow ref. [78], where the GLF is built from the radio luminosity function described in ref. [79]. In a process similar to the one taken for SFGs, we consider the contribution to the cross-correlation signal of this reference mAGN model including a free normalization $A_{\text{mAGN}}$, which alters eq. (2.9) as:

$$ \Xi_{\text{phys}}^{ar}(\theta) \langle I_a \rangle = A_{\text{BLZ}}^{1h} \hat{\Xi}_{\text{BLZ},1h}^{ar}(\theta, \mu_{\text{BLZ}}, p_1) + A_{\text{BLZ}}^{2h} \hat{\Xi}_{\text{BLZ},2h}^{ar}(\theta, \mu_{\text{BLZ}}, p_1) + A_{\text{mAGN}} \hat{\Xi}_{\text{mAGN}}^{ar}(\theta). \quad (C.1) $$

The constraints on the new set of parameters are reported in the left panel of figure 12. The inclusion of mAGN does not improve the SNR, which remains stable at 7.2 as in the BLZ-only case (see table 2). As discussed for SFGs, this result meets expectations considering that the reference mAGN model has a $\gamma$-ray energy spectrum described by a power-law with spectral index of 2.3. In fact, such an index, while providing a relatively harder spectrum with respect
<figure><img src="image_15.png" /><figcaption>Figure 13. Left: constraints on the normalization parameters and spectral indices in a model containing blazar and potential contributions from DM annihilation. Right: a description of the energy dependence of each component.</figcaption></figure>

to SFGs, describes a significantly softer spectrum compared to the best-fit value we obtain when describing the data with the phenomenological model. We conclude that both SFGs and mAGNs, as modelled under the aforementioned reference cases, can be considered to have negligible contributions to the cross-correlations signal compared to blazars. We note that also the UGRB small scale anisotropy in the LAT band is heavily dominated by blazars, as pointed out in ref. [27]. Future UGRB measurements extending at lower energies may open up new regimes in which SFG and/or mAGN can be expected to contribute more.

## D WIMP dark matter

Finally, we present a case including the astrophysical contribution from blazars as described in eq. (2.9) together with an additional term from particle DM annihilation. Thereby the model reads:

$$ \Xi_{\text{phys}}^{\text{ar}}(\theta) \langle I_a \rangle = A_{\text{BLZ}}^{\text{1h}} \hat{\Xi}_{\text{BLZ},1\text{h}}^{\text{ar}}(\theta, \mu_{\text{BLZ}}, p_1) + A_{\text{BLZ}}^{\text{2h}} \hat{\Xi}_{\text{BLZ},2\text{h}}^{\text{ar}}(\theta, \mu_{\text{BLZ}}, p_1) + A_{\text{DM}} \hat{\Xi}_{\text{DM}}^{\text{ar}}(\theta; m_{\text{DM}}), \quad (\text{D.1}) $$

where the DM contribution depends on the DM mass $m_{\text{DM}}$ along with the velocity-averaged annihilation rate $\langle\sigma_{\text{ann}}v\rangle$ expressed in terms of the “thermal” annihilation cross-section $\langle\sigma_{\text{ann}}v\rangle_{\text{th}} = 3 \times 10^{-26} \text{ cm}^{-2} \text{ s}^{-1}$, i.e., $A_{\text{DM}} = \langle\sigma_{\text{ann}}v\rangle/\langle\sigma_{\text{ann}}v\rangle_{\text{th}}$. As a benchmark final state of annihilation, we adopt the $b\bar{b}$ channel. All the ingredients entering the computation of the DM contribution are detailed in ref. [8]. The triangle plot is shown in figure 13. We find the best-fit annihilation cross-section and mass to be $\langle\sigma_{\text{ann}}v\rangle = 32_{-8}^{+10}\langle\sigma_{\text{ann}}v\rangle_{\text{th}}$, and $m_{\text{DM}} = 363_{-39.5}^{+138}$ GeV, which is also in agreement with results for the $b\bar{b}$ annihilation channel found in ref. [8]. We find that the SNR for this physical model is 8.9, which is analogous to the value obtained by the log-parabolic phenomenological model. Quantitatively, the $\Delta\chi^2 = 27$ between the DM-inclusive and a BLZ-only model would point towards a significant preference in favor of the former. However, there are two caveats related to such an inference.
<table><thead><tr><th></th><th>W</th><th>X</th><th>Y</th><th>Z</th></tr></thead><tbody><tr><td>Δχ<sup>2</sup></td><td>2.07</td><td>4.78</td><td>3.36</td><td>78.64</td></tr><tr><td>SNR</td><td>1.44</td><td>2.19</td><td>1.83</td><td>8.86</td></tr></tbody></table>

**Table 5.** The SNR and $\Delta\chi^2$ results for the log-parabola phenomenological model with respect to the null hypothesis for each data vector in the blind analysis. Only Z shows a large preference for the models, with X and Y compatible with null signals.

First, as discussed in the main text, the DM component might be effectively fitting a curved $\gamma$-ray spectrum which is actually produced by a different astrophysical mechanism. It is also quite unlikely that the astrophysical contribution becomes as subdominant as in figure 13, namely, this would require a peculiar redshift behaviour of the UGRB sources to suppress their correlation with lensing. Second, the fit requires a large $A_{DM}$, and such large values of the annihilation cross-section can be in tension with constraints from other probes, e.g., dwarf spheroidal galaxies [80] and the Large Magellanic Cloud [81]. On the other hand, the value of $A_{DM}$ has to be taken with a grain of salt, since there is an uncertainty in the overall amplitude of the cosmological DM signal, primarily associated to the poorly known contribution of subhalos (see for example ref. [82]).

A follow-up analysis, including DES Y3 galaxy clustering and more statistics will help in clarifying the viability of a particle DM contribution.

## E Data verification and robustness checks

In order to make sure that the cross-correlation described in eq. (4.1) is correct,<fnref n="13" /> we have tested the estimator using four datasets, three of which were dominated by noise while the remaining dataset contained the cross-correlation signal. The noise dominated datasets were chosen as: (i) the cross component of the galaxy shapes (known as B-modes or $\gamma_\times$), which should yield a null detection when cross-correlated with Fermi-LAT $\gamma$-ray maps; (ii) a realization of the shape noise that we obtained by randomly rotating the tangential components of the signal and cross-correlating it with the $\gamma$-ray maps; (iii) randomly reshuffled pixels of the masked $\gamma$-ray maps cross-correlated with the tangential ellipticity. Along with the true signal, the four datasets were blinded and randomly assigned the names W, X, Y, and Z. The results of a full MCMC analysis considering the log-parabola model as a fit for each of the blinded datasets (using the same covariance that was used in the data analysis) are reported in table 5. We find that W, X, and Y yield extremely low SNRs and $\Delta\chi^2$ values, while Z has values that correspond to those obtained in the data analysis in table 2. The results obtained here are compatible with our expectations: there are no spurious detections, and the presence of a signal occurs only where it is truly possible, with clear distinctions between the null and true signals. We note here that the estimator was also tested through an extensive blinding procedure in ref. [8], and has been consequently verified to be a robust method to determine the cross-correlation signal.

<footnote n="13">¹³The masked UGRB maps can be found on Github as: https://github.com/BhashinT/Fermi-LAT-Gamma-Ray-Maps-12-years-.</footnote>
<figure><img src="image_16.png" /><figcaption>Figure 14. Masks used for the “quadrant” test.</figcaption></figure>

We now turn our attention to the ninth energy bin in the dataset. As shown in table 1, the Fermi-LAT statistics are scarce in the final energy bin, leading to a large shot noise component, which can prevent a reliable signal detection. Therefore, the signal that we have obtained for this energy range could be the consequence of lurking pixels or sets of pixels. In order to test this, we conduct a “quadrant” test that proceeds as follows: we mask the γ-ray map using a combined mask of the DES Y3 footprint as well as the foreground subtracted Fermi-LAT map, and further divide the resultant mask into four parts, each of which covers half of the DES footprint, but oriented differently, as shown in figure 14. We then measured the amplitudes of the resultant cross-correlation signals as detailed in section 4. The covariance used for this test has been obtained using the same method described in appendix A, i.e., combining shape noise simulations with Gaussian realisations of the large-scale structure. The resulting amplitudes for each of these quadrants, as a function of the photon counts in each quadrant, are shown in figure 15. A mild scaling of the errors with the number of counts is present, with the strongest constraints found for the highest count number. The amplitudes, are very consistent in size, confirming that the signal we have obtained does not appear to be a spurious effect brought in by particular group of pixels, but appears to be a true cross-correlation signal.

Finally, in figure 16 we report the breakdown, in terms of energy and redshift bins, of the measured cross-correlations signal, together with the best-fit power-law phenomenological model. By looking at the different panels, we notice that the signal appears to become stronger as the energy and redshift increase, with the lowest energy-redshift bin seeming to possess only a visible 2-halo component, showing that the SNR for the total model (i.e. the sum of the 1-halo and 2-halo components) is higher at larger angular bins as depicted in table 2. We find here, visually, that as we move to larger energy and redshift bins, the 1-halo terms start to emerge, as the signal on the lower end of the angular scale becomes stronger relative to the lower energies and redshifts.
<figure><img src="image_17.png" /><figcaption>**Figure 15.** The amplitudes as a function of the γ-ray counts for each of the quadrants depicted in figure 14.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Phenomenological model fits for the cross-correlations between the tangential shear and photon flux intensities, based on the estimator described in eq. (4.1).</figcaption></figure>

References

[1] PLANCK collaboration, Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 (2020) A6 [Erratum ibid. 652 (2021) C4] [arXiv:1807.06209] [INSPIRE].

[2] P. Ullio, L. Bergström, J. Edsjö and C.G. Lacey, Cosmological dark matter annihilations into gamma-rays — a closer look, Phys. Rev. D 66 (2002) 123502 [astro-ph/0207125] [INSPIRE].
[3] S. Camera, M. Fornasa, N. Fornengo and M. Regis, *A Novel Approach in the Weakly Interacting Massive Particle Quest: Cross-correlation of Gamma-Ray Anisotropies and Cosmic Shear*, Astrophys. J. Lett. **771** (2013) L5 [arXiv:1212.5018] [INSPIRE].

[4] S. Camera, M. Fornasa, N. Fornengo and M. Regis, *Tomographic-spectral approach for dark matter detection in the cross-correlation between cosmic shear and diffuse γ-ray emission*, JCAP **06** (2015) 029 [arXiv:1411.4651] [INSPIRE].

[5] M. Shirasaki, S. Horiuchi and N. Yoshida, *Cross-Correlation of Cosmic Shear and Extragalactic Gamma-ray Background: Constraints on the Dark Matter Annihilation Cross-Section*, Phys. Rev. D **90** (2014) 063502 [arXiv:1404.5503] [INSPIRE].

[6] M. Shirasaki, O. Macias, S. Horiuchi, S. Shirai and N. Yoshida, *Cosmological constraints on dark matter annihilation and decay: Cross-correlation analysis of the extragalactic γ-ray background and cosmic shear*, Phys. Rev. D **94** (2016) 063522 [arXiv:1607.02187] [INSPIRE].

[7] T. Tröster et al., *Cross-correlation of weak lensing and gamma rays: implications for the nature of dark matter*, Mon. Not. Roy. Astron. Soc. **467** (2017) 2706 [arXiv:1611.03554] [INSPIRE].

[8] DES collaboration, *Detection of Cross-Correlation between Gravitational Lensing and γ Rays*, Phys. Rev. Lett. **124** (2020) 101102 [arXiv:1907.13484] [INSPIRE].

[9] J.-Q. Xia, A. Cuoco, E. Branchini and M. Viel, *Tomography of the Fermi-lat γ-ray Diffuse Extragalactic Signal via Cross Correlations With Galaxy Catalogs*, Astrophys. J. Suppl. **217** (2015) 15 [arXiv:1503.05918] [INSPIRE].

[10] M. Regis, J.-Q. Xia, A. Cuoco, E. Branchini, N. Fornengo and M. Viel, *Particle dark matter searches outside the Local Group*, Phys. Rev. Lett. **114** (2015) 241301 [arXiv:1503.05922] [INSPIRE].

[11] A. Cuoco, J.-Q. Xia, M. Regis, E. Branchini, N. Fornengo and M. Viel, *Dark Matter Searches in the Gamma-ray Extragalactic Background via Cross-correlations With Galaxy Catalogs*, Astrophys. J. Suppl. **221** (2015) 29 [arXiv:1506.01030] [INSPIRE].

[12] M. Shirasaki, S. Horiuchi and N. Yoshida, *Cross-Correlation of the Extragalactic Gamma-ray Background with Luminous Red Galaxies*, Phys. Rev. D **92** (2015) 123540 [arXiv:1511.07092] [INSPIRE].

[13] A. Cuoco, M. Bilicki, J.-Q. Xia and E. Branchini, *Tomographic imaging of the Fermi-LAT gamma-ray sky through cross-correlations: A wider and deeper look*, Astrophys. J. Suppl. **232** (2017) 10 [arXiv:1709.01940] [INSPIRE].

[14] S. Ammazzalorso, N. Fornengo, S. Horiuchi and M. Regis, *Characterizing the local gamma-ray Universe via angular cross-correlations*, Phys. Rev. D **98** (2018) 103007 [arXiv:1808.09225] [INSPIRE].

[15] A. Paopiamsap, D. Alonso, D.J. Bartlett and M. Bilicki, *Constraints on dark matter and astrophysics from tomographic γ-ray cross-correlations*, Phys. Rev. D **109** (2020) 103517 [arXiv:2307.14881] [INSPIRE].

[16] E. Branchini et al., *Cross-correlating the γ-ray sky with Catalogs of Galaxy Clusters*, Astrophys. J. Suppl. **228** (2017) 8 [arXiv:1612.05788] [INSPIRE].

[17] M. Shirasaki, O. Macias, S. Horiuchi, N. Yoshida, C.-H. Lee and A.J. Nishizawa, *Correlation of extragalactic γ rays with cosmic matter density distributions from weak gravitational lensing*, Phys. Rev. D **97** (2018) 123015 [arXiv:1802.10257] [INSPIRE].

[18] D. Hashimoto et al., *Measurement of redshift-dependent cross-correlation of HSC clusters and Fermi γ-rays*, Mon. Not. Roy. Astron. Soc. **484** (2019) 5256 [arXiv:1805.08139] [INSPIRE].
[19] M. Colavincenzo et al., *Searching for gamma-ray emission from galaxy clusters at low redshift*, Mon. Not. Roy. Astron. Soc. 491 (2020) 3225 [arXiv:1907.05264] [INSPIRE].

[20] X. Tan, M. Colavincenzo and S. Ammazzalorso, *Bounds on WIMP dark matter from galaxy clusters at low redshift*, Mon. Not. Roy. Astron. Soc. 495 (2020) 114 [arXiv:1907.06905] [INSPIRE].

[21] N. Fornengo, L. Perotto, M. Regis and S. Camera, *Evidence of Cross-correlation between the CMB Lensing and the Γ-ray sky*, Astrophys. J. Lett. 802 (2015) L1 [arXiv:1410.4997] [INSPIRE].

[22] S. Ando, A. Benoit-Lévy and E. Komatsu, *Mapping dark matter in the gamma-ray sky with galaxy catalogs*, Phys. Rev. D 90 (2014) 023514 [arXiv:1312.4403] [INSPIRE].

[23] N. Fornengo and M. Regis, *Particle dark matter searches in the anisotropic sky*, Front. Phys. 2 (2014) 6 [arXiv:1312.4835] [INSPIRE].

[24] M. Fornasa et al., *Angular power spectrum of the diffuse gamma-ray emission as measured by the Fermi Large Area Telescope and constraints on its dark matter interpretation*, Phys. Rev. D 94 (2016) 123005 [arXiv:1608.07289] [INSPIRE].

[25] C. Feng, A. Cooray and B. Keating, *Planck Lensing and Cosmic Infrared Background Cross-Correlation with Fermi-LAT: Tracing Dark Matter Signals in the Gamma-Ray Background*, Astrophys. J. 836 (2017) 127 [arXiv:1608.04351] [INSPIRE].

[26] M. Ajello et al., *The Origin of the Extragalactic Gamma-Ray Background and Implications for Dark-Matter Annihilation*, Astrophys. J. Lett. 800 (2015) L27 [arXiv:1501.05301] [INSPIRE].

[27] M. Korsmeier, E. Pinetti, M. Negro, M. Regis and N. Fornengo, *Flat-spectrum Radio Quasars and BL Lacs Dominate the Anisotropy of the Unresolved Gamma-Ray Background*, Astrophys. J. 933 (2022) 221 [arXiv:2201.02634] [INSPIRE].

[28] DES collaboration, *The Dark Energy Survey and Operations: Year 1*, in *Observatory Operations: Strategies, Processes, and Systems V*, proceedings of the SPIE Astronomical Telescopes + Instrumentation, Montréal, Quebec, Canada, 22–27 June 2014, pp. 332–346 [Proc. SPIE 9149 (2014) 91490V] [INSPIRE].

[29] *FERMI-LAT* collaboration, *Fermi Large Area Telescope Fourth Source Catalog*, Astrophys. J. Suppl. 247 (2020) 33 [arXiv:1902.10045] [INSPIRE].

[30] S. Abdollahi et al., *Incremental Fermi Large Area Telescope Fourth Source Catalog*, Astrophys. J. Suppl. 260 (2022) 53 [arXiv:2201.11184] [INSPIRE].

[31] DES collaboration, *Dark Energy Survey year 3 results: covariance modelling and its impact on parameter estimation and quality of fit*, Mon. Not. Roy. Astron. Soc. 508 (2021) 3125 [arXiv:2012.08568] [INSPIRE].

[32] *FERMI-LAT* collaboration, *Unresolved Gamma-Ray Sky through its Angular Power Spectrum*, Phys. Rev. Lett. 121 (2018) 241101 [arXiv:1812.02079] [INSPIRE].

[33] A. Cooray and R.K. Sheth, *Halo Models of Large Scale Structure*, Phys. Rep. 372 (2002) 1 [astro-ph/0206508] [INSPIRE].

[34] D.J. Eisenstein and W. Hu, *Power spectra for cold dark matter and its variants*, Astrophys. J. 511 (1997) 5 [astro-ph/9710252] [INSPIRE].

[35] M. Asgari, A.J. Mead and C. Heymans, *The halo model for cosmology: a pedagogical review*, Open J. Astrophys. 6 (2023) 39 [arXiv:2303.08752] [INSPIRE].

[36] DES collaboration, *Dark Energy Survey Year 3 results: Cosmological constraints from galaxy clustering and weak lensing*, Phys. Rev. D 105 (2022) 023520 [arXiv:2105.13549] [INSPIRE].
[37] R. Mukherjee et al., *EGRET observations of high-energy gamma-ray emission from blazars: An update*, Astrophys. J. 490 (1997) 116 [INSPIRE].

[38] B. Flaugher et al., *The Dark Energy Camera, Astron. J. 150 (2015) 150* [arXiv:1504.02900] [INSPIRE].

[39] M. Gatti et al., *Dark energy survey year 3 results: weak lensing shape catalogue*, Mon. Not. Roy. Astron. Soc. 504 (2021) 4312 [arXiv:2011.03408] [INSPIRE].

[40] DES collaboration, *Dark Energy Survey Y3 results: blending shear and redshift biases in image simulations*, Mon. Not. Roy. Astron. Soc. 509 (2022) 3371 [arXiv:2012.08567] [INSPIRE].

[41] E.S. Sheldon and E.M. Huff, *Practical Weak Lensing Shear Measurement with Metacalibration*, Astrophys. J. 841 (2017) 24 [arXiv:1702.02601] [INSPIRE].

[42] E. Huff and R. Mandelbaum, *Metacalibration: Direct Self-Calibration of Biases in Shear Measurement*, arXiv:1702.02600 [INSPIRE].

[43] I. Sevilla-Noarbe et al., *Dark Energy Survey Year 3 Results: Photometric Data Set for Cosmology*, Astrophys. J. Suppl. 254 (2021) 24 [arXiv:2011.03407] [INSPIRE].

[44] DES collaboration, *Phenotypic redshifts with self-organizing maps: A novel method to characterize redshift distributions of source galaxies for weak lensing*, Mon. Not. Roy. Astron. Soc. 489 (2019) 820 [arXiv:1901.05005] [INSPIRE].

[45] DES collaboration, *Dark Energy Survey Year 3 Results: clustering redshifts — calibration of the weak lensing source redshift distributions with redMaGiC and BOSS/eBOSS*, Mon. Not. Roy. Astron. Soc. 510 (2022) 1223 [arXiv:2012.08569] [INSPIRE].

[46] S. Everett et al., *Dark Energy Survey Year 3 Results: Measuring the Survey Transfer Function with Balrog*, Astrophys. J. Suppl. 258 (2022) 15 [arXiv:2012.12825] [INSPIRE].

[47] E. Suchyta et al., *No galaxy left behind: accurate measurements with the faintest objects in the Dark Energy Survey*, Mon. Not. Roy. Astron. Soc. 457 (2016) 786 [arXiv:1507.08336] [INSPIRE].

[48] DES collaboration, *Dark Energy Survey Year 3 results: Exploiting small-scale information with lensing shear ratios*, Phys. Rev. D 105 (2022) 083529 [arXiv:2105.13542] [INSPIRE].

[49] J. Myles et al., *Dark Energy Survey Year 3 results: redshift calibration of the weak lensing source galaxies*, Mon. Not. Roy. Astron. Soc. 505 (2021) 4249 [arXiv:2012.08566] [INSPIRE].

[50] DES collaboration, *Dark Energy Survey Year 3: Blue Shear*, arXiv:2410.22272 [INSPIRE].

[51] DES collaboration, *Dark Energy Survey Year 3 results: Curved-sky weak lensing mass map reconstruction*, Mon. Not. Roy. Astron. Soc. 505 (2021) 4626 [arXiv:2105.13539] [INSPIRE].

[52] M. Ackermann et al., *The spectrum of isotropic diffuse gamma-ray emission between 100 MeV and 820 GeV*, Astrophys. J. 799 (2015) 86 [arXiv:1410.3696] [INSPIRE].

[53] K.M. Górski et al., *HEALPix — A Framework for high resolution discretization, and fast analysis of data distributed on the sphere*, Astrophys. J. 622 (2005) 759 [astro-ph/0409513] [INSPIRE].

[54] DES collaboration, *Density Split Statistics: Cosmological Constraints from Counts and Lensing in Cells in DES Y1 and SDSS Data*, Phys. Rev. D 98 (2018) 023507 [arXiv:1710.05045] [INSPIRE].

[55] S. Singh, R. Mandelbaum, U. Seljak, A. Slosar and J. Vazquez Gonzalez, *Galaxy-galaxy lensing estimators and their covariance properties*, Mon. Not. Roy. Astron. Soc. 471 (2017) 3827 [arXiv:1611.00752] [INSPIRE].
[56] M. Jarvis, G. Bernstein and B. Jain, *The skewness of the aperture mass statistic, Mon. Not. Roy. Astron. Soc.* **352** (2004) 338 [astro-ph/0307393] [INSPIRE].

[57] M. Jarvis, *Treecorr: Two-point correlation functions*, in *Astrophysics Source Code Library*, ascl:1508.007 (2015).

[58] D. Foreman-Mackey, D.W. Hogg, D. Lang and J. Goodman, *emcee: The MCMC Hammer*, *Publ. Astron. Soc. Pac.* **125** (2013) 306 [arXiv:1202.3665] [INSPIRE].

[59] S.R. Hinton, *Chainconsumer*, *J. Open Source Softw.* **1** (2016) 00045.

[60] J. Goodman and J. Weare, *Ensemble samplers with affine invariance*, *Comm. Appl. Math. Comput. Sci.* **5** (2010) 65 [INSPIRE].

[61] DES collaboration, *Cosmic shear measurements with Dark Energy Survey Science Verification data*, *Phys. Rev. D* **94** (2016) 022002 [arXiv:1507.05598] [INSPIRE].

[62] D. Guetta et al., *TeV emission from Gamma Ray Bursts*, checking the hadronic model, *Eur. Phys. J. Web Conf.* **280** (2023) 01005 [INSPIRE].

[63] F.-S. Min, Y.-H. Yao, R.-Y. Liu, S. Chen, H. Lu and Y.-Q. Guo, *Contribution of γ-Ray Burst Afterglow Emissions to the Isotropic Diffuse γ-Ray Background*, *Astrophys. J.* **964** (2024) 195 [INSPIRE].

[64] Y.-H. Yao et al., *Contribution of High-energy GRB Emissions to the Spectrum of the Isotropic Diffuse γ-Ray Background*, *Astrophys. J.* **901** (2020) 106 [INSPIRE].

[65] M.A. Roth, M.R. Krumholz, R.M. Crocker and S. Celli, *The diffuse γ-ray background is dominated by star-forming galaxies*, *Nature* **597** (2021) 341 [arXiv:2109.07598] [INSPIRE].

[66] FERMI-LAT collaboration, *3FHL: The Third Catalog of Hard Fermi-LAT Sources*, *Astrophys. J. Suppl.* **232** (2017) 18 [arXiv:1702.00664] [INSPIRE].

[67] S. Hussain, R.A. Batista, E.M. de Gouveia Dal Pino and K. Dolag, *The diffuse gamma-ray flux from clusters of galaxies*, *Nat. Commun.* **14** (2023) 2486 [arXiv:2203.01260] [INSPIRE].

[68] A. Cooray, *Extragalactic Background Light: Measurements and Applications*, *Roy. Soc. Open Sci.* **3** (2016) 150555 [arXiv:1602.03512] [INSPIRE].

[69] J.D. Finke, S. Razzaque and C.D. Dermer, *Modeling the Extragalactic Background Light from Stars and Dust*, *Astrophys. J.* **712** (2010) 238 [arXiv:0905.1115] [INSPIRE].

[70] J.D. Finke et al., *Modeling the Extragalactic Background Light and the Cosmic Star Formation History*, *Astrophys. J.* **941** (2022) 33 [arXiv:2210.01157] [INSPIRE].

[71] F.W. Stecker, S.T. Scully and M.A. Malkan, *An Empirical Determination of the Intergalactic Background Light from UV to FIR Wavelengths Using FIR Deep Galaxy Surveys and the Gamma-ray Opacity of the Universe*, *Astrophys. J.* **827** (2016) 6 [Erratum ibid. **863** (2018) 112] [arXiv:1605.01382] [INSPIRE].

[72] V. Allevato, A. Finoguenov and N. Cappelluti, *Clustering of γ-ray selected 2LAC Fermi Blazars*, *Astrophys. J.* **797** (2014) 96 [arXiv:1410.0358] [INSPIRE].

[73] C. Hikage et al., *Cosmology from cosmic shear power spectra with Subaru Hyper Suprime-Cam first-year data*, *Publ. Astron. Soc. Jpn.* **71** (2019) 43 [arXiv:1809.09148] [INSPIRE].

[74] R. Murata et al., *Constraints on the mass-richness relation from the abundance and weak lensing of SDSS clusters*, *Astrophys. J.* **854** (2018) 120 [arXiv:1707.01907] [INSPIRE].
[75] J. Hartlap, P. Simon and P. Schneider, *Why your model parameter confidences might be too optimistic: Unbiased estimation of the inverse covariance matrix*, *Astron. Astrophys.* **464** (2007) 399 [astro-ph/0608064] [INSPIRE].

[76] DES collaboration, *Dark energy survey year 3 results: High-precision measurement and modeling of galaxy-galaxy lensing*, *Phys. Rev. D* **105** (2022) 083528 [arXiv:2105.13541] [INSPIRE].

[77] M.A. Troxel et al., *Survey geometry and the internal consistency of recent cosmic shear measurements*, *Mon. Not. Roy. Astron. Soc.* **479** (2018) 4998 [arXiv:1804.10663] [INSPIRE].

[78] F.W. Stecker, C.R. Shrader and M.A. Malkan, *The Extragalactic Gamma-Ray Background from Core Dominated Radio Galaxies*, *Astrophys. J.* **879** (2019) 68 [arXiv:1903.06544] [INSPIRE].

[79] Z. Yuan, J. Wang, D. Worrall, B.-B. Zhang and J. Mao, *Determining the core radio luminosity function of radio AGNs via copula*, *Astrophys. J. Suppl.* **239** (2018) 33 [arXiv:1810.12713].

[80] FERMI-LAT collaboration, *Searching for Dark Matter Annihilation from Milky Way Dwarf Spheroidal Galaxies with Six Years of Fermi Large Area Telescope Data*, *Phys. Rev. Lett.* **115** (2015) 231301 [arXiv:1503.02641] [INSPIRE].

[81] M. Regis et al., *The EMU view of the Large Magellanic Cloud: troubles for sub-TeV WIMPs*, *JCAP* **11** (2021) 046 [arXiv:2106.08025] [INSPIRE].

[82] FERMI-LAT collaboration, *Limits on Dark Matter Annihilation Signals from the Fermi LAT 4-year Measurement of the Isotropic Gamma-Ray Background*, *JCAP* **09** (2015) 008 [arXiv:1501.05464] [INSPIRE].
Author List

B. Thakore $^{a,b,c}$, M. Negro $^{d}$, M. Regis $^{a,b}$, S. Camera $^{a,b,e,f}$, D. Gruen $^{c,g}$, N. Fornengo $^{a,b}$,
A. Roodman $^{h,i}$, A. Porredon $^{j,k}$, T. Schutt $^{h,i,l}$, A. Cuoco $^{a,b}$, A. Alarcon $^{m}$, A. Amon $^{n}$,
K. Bechtol $^{o}$, M.R. Becker $^{p}$, G.M. Bernstein $^{q}$, A. Campos $^{r,s}$, A. Carnero Rosell $^{t,u,v}$,
M. Carrasco Kind $^{w,x}$, R. Cawthon $^{y}$, C. Chang $^{z,aa}$, R. Chen $^{ab}$, A. Choi $^{ac}$, J. Cordero $^{ad}$,
C. Davis$^{i}$, J. DeRose $^{ae}$, H.T. Diehl $^{af}$, S. Dodelson $^{z,aa,af}$, C. Doux $^{q,ag}$,
A. Drlica-Wagner $^{z,aa,af}$, K. Eckert $^{q}$, J. Elvin-Poole $^{ah}$, S. Everett $^{ai}$, A. Ferté $^{h}$,
M. Gatti $^{q}$, G. Giannini $^{aa,aj}$, R.A. Gruendl $^{w,x}$, I. Harrison $^{ak}$, W.G. Hartley$^{al}$,
E.M. Huff $^{am}$, M. Jarvis $^{q}$, N. Kuropatkin $^{af}$, P.-F. Leget $^{i}$, N. MacCrann $^{an}$,
J. McCullough $^{c,h,i,n}$, J. Myles $^{n}$, A. Navarro-Alsina $^{ao}$, S. Pandey $^{q}$, J. Prat $^{z,ap}$,
M. Raveri $^{aq}$, R.P. Rollins $^{ad}$, A.J. Ross $^{ar}$, E.S. Rykoff $^{h,i}$, C. Sánchez $^{q}$, L.F. Secco $^{aa}$,
I. Sevilla-Noarbe $^{k}$, E. Sheldon $^{as}$, T. Shin $^{at}$, M.A. Troxel $^{ab}$, I. Tutusaus $^{au}$, B. Yanny $^{af}$,
B. Yin $^{r}$, Y. Zhang $^{av}$, M. Aguena $^{u}$, D. Brooks $^{aw}$, J. Carretero $^{aj}$, L.N. da Costa $^{u}$,
T.M. Davis $^{ax}$, J. De Vicente $^{k}$, S. Desai $^{ay}$, P. Doel $^{aw}$, B. Flaugher $^{af}$, J. Frieman $^{aa,af}$,
J. García-Bellido $^{az}$, E. Gaztanaga $^{m,ba,bb}$, G. Gutierrez $^{af}$, S.R. Hinton $^{ax}$, D.L. Hollowood $^{bc}$,
K. Honscheid $^{ar,bd}$, D.J. James $^{be}$, K. Kuehn $^{bf,bg}$, O. Lahav $^{aw}$, S. Lee $^{am}$, M. Lima $^{u,bh}$,
J.L. Marshall $^{bi}$, J. Mena-Fernández $^{bj}$, R. Miquel $^{aj,bk}$, R.L.C. Ogando $^{bl}$, A. Palmese $^{r}$,
A. Pieres $^{u,bl}$, A.A. Plazas Malagón $^{h,i}$, S. Samuroff $^{aj,bm}$, E. Sanchez $^{k}$, D. Sanchez Cid $^{k}$,
M. Smith $^{bn}$, E. Suchyta $^{bo}$, G. Tarle $^{bp}$, V. Vikram $^{bq}$, A.R. Walker $^{av}$, N. Weaverdyck $^{ae,br}$

<footnote><sup>a</sup> Dipartimento di Fisica, Università degli Studi di Torino,
via P. Giuria 1, 10125 Torino, Italy</footnote>

<footnote><sup>b</sup> INFN — Istituto Nazionale di Fisica Nucleare, Sezione di Torino,
via P. Giuria 1, 10125 Torino, Italy</footnote>

<footnote><sup>c</sup> University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität,
Scheinerstr. 1, 81679 Munich, Germany</footnote>

<footnote><sup>d</sup> Department of Physics & Astronomy, Nicholson Hall, Tower Drive, Louisiana State University,
Baton Rouge, LA 70803, U.S.A.</footnote>

<footnote><sup>e</sup> Department of Physics & Astronomy, University of the Western Cape,
Robert Sobukwe Road, Bellville 7535, Cape Town, South Africa</footnote>

<footnote><sup>f</sup> INAF — Istituto Nazionale di Astrofisica, Osservatorio Astrofisico di Torino,
Strada Osservatorio 20, 10025 Pino Torinese, Italy</footnote>

<footnote><sup>g</sup> Excellence Cluster ORIGINS,
Boltzmannstr. 2, 85748 Garching, Germany</footnote>

<footnote><sup>h</sup> SLAC National Accelerator Laboratory,
2575 Sand Hill Road, Menlo Park, CA 94025 U.S.A.</footnote>

<footnote><sup>i</sup> Kavli Institute for Particle Astrophysics & Cosmology,
P.O. Box 2450, Stanford University, Stanford, CA 94305, U.S.A.</footnote>

<footnote><sup>j</sup> Astronomical Institute, Faculty of Physics and Astronomy, Ruhr University Bochum,
Universitätsstraße 150, 44780 Bochum, Germany</footnote>

<footnote><sup>k</sup> Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT),
Av. Complutense, 40, Moncloa, Aravaca, 28040, Madrid, Spain</footnote>

<footnote><sup>l</sup> Department of Physics, Stanford University,
382 Via Pueblo Mall, Stanford, CA 94305, U.S.A.</footnote>

<footnote><sup>m</sup> Institute of Space Sciences (ICE, CSIC), Campus UAB,
Carrer de Can Magrans, s/n, 08193 Barcelona, Spain</footnote>

<footnote><sup>n</sup> Department of Astrophysical Sciences, Princeton University,
Peyton Hall, Princeton, NJ 08544, U.S.A.</footnote>
o Physics Department, 2320 Chamberlin Hall, University of Wisconsin-Madison,
1150 University Avenue, Madison, WI 53706-1390, U.S.A.

p Argonne National Laboratory,
9700 South Cass Avenue, Lemont, IL 60439, U.S.A.

q Department of Physics and Astronomy, University of Pennsylvania,
209 South 33rd Street, Philadelphia, PA 19104, U.S.A.

r Department of Physics, Carnegie Mellon University,
5000 Forbes Avenue, Pittsburgh, PA 15213, U.S.A.

s NSF AI Planning Institute for Physics of the Future, Carnegie Mellon University,
5000 Forbes Avenue, Pittsburgh, PA 15213, U.S.A.

t Instituto de Astrofísica de Canarias,
E-38205 La Laguna, Tenerife, Spain

u Laboratório Interinstitucional de e-Astronomia — LINEA,
Rua Gal. José Cristino 77, Rio de Janeiro, RJ 20921-400, Brazil

v Universidad de La Laguna, Departamento de Astrofísica,
E-38206 La Laguna, Tenerife, Spain

w Center for Astrophysical Surveys, National Center for Supercomputing Applications,
1205 West Clark St., Urbana, IL 61801, U.S.A.

x Department of Astronomy, University of Illinois at Urbana-Champaign,
1002 W. Green Street, Urbana, IL 61801, U.S.A.

y Physics Department, William Jewell College,
500 College Hill, Liberty, MO 64068, U.S.A.

z Department of Astronomy and Astrophysics, University of Chicago,
5640 South Ellis Avenue, Chicago, IL 60637, U.S.A.

aa Kavli Institute for Cosmological Physics, University of Chicago,
933 East 56th Street, Chicago, IL 60637, U.S.A.

ab Department of Physics, Duke University,
120 Science Drive, Durham, NC 27708, U.S.A.

ac NASA Goddard Space Flight Center,
8800 Greenbelt Rd, Greenbelt, MD 20771, U.S.A.

ad Jodrell Bank Center for Astrophysics, School of Physics and Astronomy, University of Manchester,
Oxford Road, Manchester, M13 9PL, U.K.

ae Lawrence Berkeley National Laboratory,
1 Cyclotron Road, Berkeley, CA 94720, U.S.A.

af Fermi National Accelerator Laboratory,
Kirk Road and Pine Street, P.O. Box 500, Batavia, IL 60510, U.S.A.

ag Université Grenoble Alpes, CNRS, LPSC-IN2P3,
38000 Grenoble, France

ah Department of Physics and Astronomy, University of Waterloo,
200 University Ave W, Waterloo, ON N2L 3G1, Canada

ai California Institute of Technology,
1200 East California Blvd, MC 249-17, Pasadena, CA 91125, U.S.A.

aj Institut de Física d'Altes Energies (IFAE),
The Barcelona Institute of Science and Technology, Campus UAB,
08193 Bellaterra (Barcelona), Spain

ak School of Physics and Astronomy, Cardiff University,
The Parade, Cardiff CF24 3AA, U.K.

al Department of Astronomy, University of Geneva,
Chemin d'Écogia 16, CH-1290 Versoix, Switzerland

am Jet Propulsion Laboratory, California Institute of Technology,
4800 Oak Grove Dr., Pasadena, CA 91109, U.S.A.

an Department of Applied Mathematics and Theoretical Physics, University of Cambridge,
Wilberforce Road, Cambridge CB3 0WA, U.K.
ao Instituto de Física Gleb Wataghin, Universidade Estadual de Campinas,
13083-859, Campinas, SP, Brazil

ap Nordita, KTH Royal Institute of Technology and Stockholm University,
Hannes Alfvéns väg 12, SE-10691 Stockholm, Sweden

aq Department of Physics, University of Genova and INFN,
Via Dodecaneso 33, 16146, Genova, Italy

ar Department of Physics,
191 W Woodruff Ave, Columbus, OH 43210, U.S.A.

as Brookhaven National Laboratory,
Bldg 510, Upton, NY 11973, U.S.A.

at Department of Physics and Astronomy, Stony Brook University,
100 Nicolls Road, Stony Brook, NY 11794, U.S.A.

au Institut de Recherche en Astrophysique et Planétologie (IRAP),
Université de Toulouse, CNRS, UPS, CNES,
14 Av. Edouard Belin, 31400 Toulouse, France

av Cerro Tololo Inter-American Observatory,
NSF's National Optical-Infrared Astronomy Research Laboratory,
Casilla 603, La Serena, Chile

aw Department of Physics & Astronomy, University College London,
Gower Street, London, WC1E 6BT, U.K.

ax School of Mathematics and Physics, University of Queensland,
Brisbane, QLD 4072, Australia

ay Department of Physics, IIT Hyderabad,
Kandi, Telangana 502285, India

az Instituto de Física Teórica UAM/CSIC, Universidad Autónoma de Madrid,
Calle Nicolás Cabrera 13-15, 28049 Madrid, Spain.

ba Institut d'Estudis Espacials de Catalunya (IEEC),
08034 Barcelona, Spain

bb Institute of Cosmology and Gravitation, University of Portsmouth,
Dennis Sciama Building, Burnaby Road, Portsmouth PO1 3FX, U.K.

bc Santa Cruz Institute for Particle Physics,
Natural Sciences 2 Bldg, 552 Red Hill Rd, Santa Cruz, CA 95064, U.S.A.

bd Department of Physics, The Ohio State University,
Columbus, OH 43210, U.S.A.

be Center for Astrophysics, Harvard & Smithsonian,
60 Garden Street, Cambridge, MA 02138, U.S.A.

bf Australian Astronomical Optics, Macquarie University,
105 Delhi Rd, North Ryde NSW 2113, Australia

bg Lowell Observatory,
1400 Mars Hill Rd, Flagstaff, AZ 86001, U.S.A.

bh Departamento de Física Matemática, Instituto de Física, Universidade de São Paulo,
CP 66318, São Paulo, SP, 05314-970, Brazil

bi George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, and
Department of Physics and Astronomy, Texas A&M University,
College Station, TX 77843, U.S.A.

bj LPSC Grenoble,
53 Avenue des Martyrs, 38026 Grenoble, France

bk Institució Catalana de Recerca i Estudis Avançats,
Passeig de Lluís Companys, 23, L'Eixample, 08010 Barcelona, Spain

bl Observatório Nacional,
Rua Gal. José Cristino 77, Rio de Janeiro, RJ 20921-400, Brazil

bm Department of Physics, Northeastern University, Dana Research Center,
100 Forsyth St #111, Boston, MA 02115, U.S.A.
bn Physics Department, Lancaster University,
Bailrigg, Lancaster LA1 4YW, U.K.

bo Computer Science and Mathematics Division, Oak Ridge National Laboratory,
One Bethel Valley Rd., Oak Ridge, TN 37831, U.S.A.

bp Department of Physics, University of Michigan,
450 Church St, Ann Arbor, MI 48109, U.S.A.

bq Central University of Kerala, Tejaswini Hills,
Periye (PO), Kasaragod (DT), Kerala, 671316, India

br Department of Astronomy, University of California, Berkeley,
501 Campbell Hall, Berkeley, CA 94720, U.S.A.