© The Authors 2026

1. Introduction

Dark matter (DM) constitutes a significant portion of the Universe, amounting to approximately 27% of the total mass-energy content according to constraints from Planck (Planck Collaboration VI 2020). Unlike baryonic matter, DM appears invisible, with its presence inferred from its gravitational effects on visible matter and the large-scale structure of the Universe. The first evidence of DM came from Fritz Zwicky’s 1930s study of the Coma Cluster (Zwicky 1937), which was followed by the discovery of flat galaxy rotation curves (Rubin & Ford 1970). Although precise measurements of the DM cosmological abundance and statistical properties of its distribution in the Universe have been achieved, its microscopic properties, such as its non-gravitational interactions and composition, are still mostly unknown.

One of the most studied candidates for DM is the weakly interacting massive particle (WIMP), a hypothetical particle that interacts only through the weak nuclear force and gravity (Bertone et al. 2005). With masses ranging from approximately 1 GeV to 1 TeV, thermally produced WIMPs in the early Universe would naturally exhibit a relic density that matches the observed DM abundance. This is due to their decoupling from the primordial plasma in the non-relativistic regime and their weak interactions with lighter standard model particles (Lee & Weinberg 1977; Gunn et al. 1978).

The weak coupling of WIMPs allowed us to test the hypothesis of WIMP DM and offers a promising method to indirectly probe it. Assuming WIMPs are the building blocks of the large-scale structure of the Universe, it is likely that WIMPs within DM haloes may annihilate or decay into detectable standard model particles. These particles, known as final state particles, can then interact to produce gamma rays, neutrinos, and cosmic rays. Searching for the primary and secondary radiation of standard model particles from the DM annihilation process is a known approach to indirectly detecting DM (see e.g. Cirelli 2012 for a review).

Various instruments have been utilised in efforts to detect DM annihilation products. These include Cherenkov telescopes such as the High Energy Stereoscopic System (H.E.S.S.; e.g. Abramowski et al. 2011, 2012; Abdallah et al. 2018), Very Energetic Radiation Imaging Telescope Array System (VERITAS; e.g. Archambault et al. 2017), and the Florian Goebel Major Atmospheric Gamma-ray Imaging Cherenkov (MAGIC) telescopes (e.g. Aleksić et al. 2011; MAGIC Collaboration 2016; Acciari et al. 2022; Abe et al. 2023); neutrino telescopes such as IceCube (Abbasi et al. 2023); the High Altitude Water Cherenkov Experiment (HAWC; e.g. Aartsen et al. 2018; albert2018dark); and the Fermi Large Area Telescope (Fermi-LAT; e.g. Ackermann et al. 2010; MAGIC Collaboration 2016; McDaniel et al. 2024; Totani 2025). These instruments study the DM annihilation process by targeting regions predicted to have high DM densities in the local Universe, including the Galactic centre (GC), and the Galactic halo and subhaloes, where dwarf spheroidal galaxies (dSphs) of the Milky Way reside. Among the DM annihilation products, gamma-rays are the most investigated because they travel without being deflected by magnetic fields, allowing them to be traced back to their source and providing crucial insights into the spatial distribution of DM. At high energies, however, gamma-rays can be partially attenuated during propagation, so their observed flux carries integrated information about the intervening cosmic medium.

The unresolved gamma-ray background (UGRB) is a significant focus in DM studies on larger scales, offering a potential probe for the indirect detection of DM across cosmological distances (e.g. Fornasa & Sánchez-Conde 2015). The UGRB emission due to DM annihilation is expected to be anisotropic, arising from the underlying large–scale distribution of DM haloes and subhaloes (see e.g. Ando & Komatsu 2006). In addition to potential DM signals, the UGRB comprises cumulative emissions from all unresolved astrophysical sources, including star-forming galaxies (SFGs; Ajello et al. 2020; Roth et al. 2021), misaligned active galactic nuclei (mAGNs; Di Mauro et al. 2013), blazars (Korsmeier et al. 2022), millisecond pulsars (Calore et al. 2014), and the interaction of cosmic rays with extragalactic background light (EBL; Kalashev et al. 2009).

Weak gravitational lensing is a powerful tool for cosmological studies (see, e.g. Bartelmann & Schneider 2001). It occurs when the gravitational field of the intervening large-scale structure induces slight distortions in the light paths of distant background objects. Weak lensing can serve as an unbiased tracer of the matter distribution, making it essential for studying the distribution and properties of DM (Camera et al. 2013). Weak gravitational lensing can be measured statistically from the apparent correlations between the ellipticities of distant galaxies. By comparing these correlations to model predictions, one can constrain cosmological parameters such as Ω_m and σ₈ (see, e.g. Secco et al. 2022; Amon et al. 2022; Dalal et al. 2023; Wright et al. 2025b; Stölzner et al. 2025).

In recent years, cross-correlations between UGRB maps and weak lensing data have become a powerful tool for probing DM properties. The first measurement of the gamma-ray shear cross-correlation was presented by Shirasaki et al. (2014). Shirasaki et al. (2016) performed a cross-correlation analysis using 7 years of Fermi-LAT gamma-ray data with weak lensing surveys from the Canada-France-Hawaii Telescope Lensing Survey (CFHTLenS) and the Red Cluster Sequence Lensing Survey (RCSLenS). Building on this, Tröster et al. (2017) cross-correlated the UGRB maps from 8 years of Fermi-LAT data with tomographic weak lensing data from CFHTLenS, RCSLenS, and the Kilo Degree Survey (KiDS-450). The work from Shirasaki et al. (2018) applied a similar approach using the first-year data from the Hyper Suprime-Cam Subaru Strategic Programme (HSC-SSP). These works are all consistent with non-detections and demonstrate that weak lensing data can provide significant information for constraining DM properties. Ammazzalorso et al. (2020), for the first time, achieved a 5.3σ detection of a cross-correlation signal between gamma-rays from 9 years of Fermi-LAT UGRB and tomographic weak lensing data from Dark Energy Survey Y1 (DES Y1). More recently, Thakore et al. (2025) detected an 8.9σ cross-correlation using 12 years of Fermi-LAT UGRB and DES Y3 lensing, dominated by large angular scales and consistent with blazar emission. While they also explored a DM component, they noted its degeneracy with astrophysical modelling and the tension with external gamma-ray limits and therefore did not interpret it as evidence of DM. We further investigate the potential reasons for the discrepancy between the DES and KiDS weak-lensing cross-correlation results in Appendix D.

In addition to weak lensing, other large-scale structure-based indirect detection approaches also provide valuable insights into DM properties. Spatial distributions of galaxies offer a high signal-to-noise ratio (S/N) for tracing matter fluctuations (Ando 2014; Cuoco et al. 2015; Regis et al. 2015; Shirasaki et al. 2015; Cuoco et al. 2017; Ammazzalorso et al. 2018; Bartlett et al. 2022; Kostić et al. 2026; Paopiamsap et al. 2024), and their cross-correlations with gamma-rays yield significant signals (e.g. Paopiamsap et al. 2024). Additionally, cross-correlation studies with galaxy clusters have great potential, as these massive structures contain large amounts of DM, and their cross-correlation with gamma-rays has also been explored as a means of detecting DM signatures (Branchini et al. 2017; Hashimoto et al. 2019; Colavincenzo et al. 2020; Di Mauro et al. 2023). Other tracers including anisotropies of the cosmic microwave background (CMB; Tan et al. 2020), CMB lensing (Fornengo et al. 2015), the cosmic infrared background (CIB; Feng et al. 2017), the late-time 21 cm signal (Pinetti et al. 2020), the thermal Sunyaev-Zel’dovich effect (Shirasaki et al. 2020), and high-energy neutrinos (Negro et al. 2023) are also used in cross-correlation studies.

Based on previous work, our study investigates the nature of DM by cross-correlating tomographic weak lensing data from the fifth data release of the Kilo Degree Survey (KiDS-Legacy, Wright et al. 2024) with the energy-binned intensity map of the UGRB from 15 years of Fermi-LAT observations (Atwood et al. 2009). Properly modelling the contribution of astrophysical sources allows the UGRB measurements to better constrain the DM component (see, e.g. Fornasa & Sánchez-Conde 2015 for a review).

This paper is structured as follows: Theoretical models are introduced in Sect. 2; we describe the data used for our analysis in Sect. 3; the data analysis methods and estimators are presented in Sect. 4; the results are discussed in Sect. 5; we summarise our conclusions in Sect. 6. We employed a flat Lambda cold dark matter (ΛCDM) cosmological model, with parameters taken from Planck Collaboration VI (2020): (h,Ω_DMh²,Ω_bh²,σ₈,n_s) = (0.677,0.119,0.022,0.810,0.967).

2. Theoretical models

Both observations, the weak gravitational lensing convergence (denoted as κ) and the gamma-ray emission intensity (denoted as g), are formalised as projections of three-dimensional fields. The general form of the projection is given by

$$\begin{array}{c}\hfill u(\widehat{\theta })={\displaystyle \int \mathrm{d}}\chi W^u(\chi )\delta (\chi \widehat{\theta },\chi ),\end{array}$$(1)

where u represents the projected field in the angular direction $\widehat{\theta }$; χ is the radial comoving distance; $\delta (\chi \widehat{\theta },\chi )$ denotes the 3D fluctuation of the field u at the coordinate $(\chi \widehat{\theta },\chi )$; W^u(χ) is the radial kernel that reflects the distribution of the field along the line-of-sight.

The theoretical model for the cross-correlation power spectrum between two such fields, u and v, can be well estimated by the Limber approximation (Limber 1953; Kaiser 1992; Kilbinger et al. 2017) when the radial kernels are significantly wider than their correlation length and at small scales (ℓ ≳ 10). It can be expressed as

$$\begin{array}{c}\hfill C_{\ell }^{\mathit{uv}}={\displaystyle \int \frac{\mathrm{d}\chi }{{\chi }^2}}{W}^u(\chi )W^v(\chi )P_{\mathit{uv}}(k=\frac{\ell +1/2}{\chi (z)},z(\chi )).\end{array}$$(2)

We assume a spatially flat ΛCDM cosmology, where z is the redshift, and P_uv(k_ℓ, z) is the cross power spectrum of the associated 3D fields u and v. It can be described as follows, with k and ℓ being the modulus of the wave number and the angular multipole, respectively:

$$\begin{array}{c}\hfill ⟨{\delta }_u(\mathit{k}){\delta }_v(\mathit{k}\prime )⟩={(2\pi )}^3{\delta }^3(\mathit{k}+\mathit{k}\prime )P_{\mathit{uv}}(k).\end{array}$$(3)

Here ⟨…⟩ indicates the ensemble average of the quantity inside the brackets. Therefore, the cross-correlation angular power spectrum for lensing convergence κ and gamma-ray intensity g can be expressed as

$$\begin{array}{c}\hfill C_{\ell }^{\mathrm{g}\kappa }={\displaystyle {\int }_{E_{\min }}^{E_{\max }}}\mathrm{d}E{\displaystyle {\int }_0^{\infty }}\mathrm{d}z\frac{c}{H(z)}\frac{1}{\chi {(z)}^2}{W}_{\mathrm{g}}(E,z)W_{\kappa }(z)P_{\mathrm{g}\kappa }(k,z),\end{array}$$(4)

where E represents the gamma-ray energy; E_min and E_max define the upper and lower limits of the energy bin being integrated; H(z) is the Hubble rate at redshift z; W_g(E, z) and W_κ(z) are the window functions for gamma-ray intensity and gravitational lensing, respectively; and P_gκ(k, z) denotes the 3D cross power spectrum.

2.1. Gravitational lensing

The window function for gravitational lensing is given by (see, e.g. Bartelmann 2010)

$$\begin{array}{c}\hfill W_{\kappa }(z)=\frac{3H_0^2{\mathrm{\Omega }}_{\mathrm{m}}}{2}(1+z)\chi (z){\displaystyle {\int }_z^{\infty }}\mathrm{d}{z}^{\prime }\frac{\chi (z^{\prime })-\chi (z)}{\chi (z^{\prime })}n(z^{\prime }),\end{array}$$(5)

where H₀ is the Hubble constant; Ω_m represents the matter density parameter of the Universe; n(z) is the redshift distribution of the source galaxies in each tomographic bin. The redshift distributions of the KiDS-Legacy gold-selected sample are shown in the top panel of Fig. 1 for the six bins used in this analysis. The integral term quantifies the lensing efficiency as a function of the source galaxies relative to the lens. The window functions for the KiDS-Legacy survey with six redshift bins are illustrated in the lower panel of Fig. 1.

Fig. 1. Top: Redshift distributions of the KiDS-Legacy gold-selected sample for tomographic bins, where the photometric redshift bin with z in [0.1, 2.0] is the sum weighted by the effective number densities listed in Table 1. The shaded vertical bands indicate the boundaries of the tomographic photometric redshift bins. Bottom: Weak lensing window functions for six tomographic bins in KiDS-Legacy.

2.2. Gamma-rays from DM

Proposed DM particles, such as WIMPS, may annihilate upon collision, creating standard model particles and high-energy photons. We assumed that DM is its own antiparticle (self-conjugate), which we take as a standard benchmark in the literature, and the annihilation window function is given by (Ando & Komatsu 2006; Fornengo & Regis 2014; Tröster et al. 2017)

$$\begin{array}{cc}\hfill W_{{\mathrm{g}}_{\mathrm{ann}}}(z,E)& =\frac{{({\mathrm{\Omega }}_{\mathrm{DM}}{\rho }_{\mathrm{c}})}^2}{4\pi }\frac{⟨{\sigma }_{\text{ann}}v⟩}{2m_{\mathrm{DM}}^2}{(1+z)}^3{\mathrm{\Delta }}^2(z)\hfill \\ \hfill & \times \frac{\mathrm{d}{N}_{\mathrm{ann}}}{\mathrm{d}E}[E(1+z)]{\mathrm{e}}^{-\tau [z,E(1+z)]},\hfill \end{array}$$(6)

where Ω_DMρ_c is the DM density in the Universe; ρ_c is the critical density of the Universe; m_DM is the rest mass of DM; ⟨σ_annv⟩ represents the thermally averaged annihilation cross-section; Δ²(z) represents the clumping factor; dN_ann/dE represents the energy spectrum of gamma-rays, indicating the number of photons emitted per unit energy per annihilation for different final states. It can be described as a function of the rest frame energy of the photons and m_DM. The optical depth τ, accounting for gamma-ray attenuation by the EBL, was taken from Franceschini et al. (2008).

We adopted the energy spectrum of the annihilation process from ‘A Poor Particle Physicist Cookbook for DM Indirect Detection’ (PPPC4DMID, Cirelli et al. 2011) with electroweak corrections (Ciafaloni et al. 2011). Our results are presented for three different final state particles: $b\overline{b}$, μ⁺μ⁻ and τ⁺τ⁻. In each of these three scenarios, we assumed that DM annihilates and decays into each final state with a branching ratio of 1. The continuous energy spectrum was obtained through two-dimensional linear interpolation for each state according to the DM mass and gamma-ray energy.

The clumping factor Δ²(z) in Eq. (6) reflects the concentration of DM within haloes and subhaloes and indicates the density distribution and the aggregation tendencies of DM in the cosmological structure. For the parameters involved in the clumping factor, we used the assumption in Tröster et al. (2017), which is defined as (Fornengo & Regis 2014)

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}^2(z)=\frac{⟨{\rho }_{\mathrm{DM}}^2⟩}{{\overline{\rho }}_{\mathrm{DM}}^2}& ={\displaystyle {\int }_{M_{\min }}^{M_{\max }}}\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}(M,z)[1+b_{\mathrm{sub}}(M,z)]\hfill \\ \hfill & \times {\displaystyle \int {\mathrm{d}}^3}\mathbf{x}\frac{{\rho }^2(\mathbf{x}|M,z)}{{\overline{\rho }}_{\mathrm{DM}}^2},\hfill \end{array}$$(7)

where M_min and M_max represent the minimal and maximal halo masses, assumed as 10⁻⁶ M_⊙ and 10¹⁸ M_⊙. As in previous works (e.g. Cuoco et al. 2015; troster2017cross), the minimum halo mass M_min = 10⁻⁶ M_⊙ was chosen to correspond to a typical WIMP free-streaming mass, below which haloes do not contribute effectively to the signal. The uncertainty associated with the selection of M_min is explored in Appendix B of Paopiamsap et al. (2024). The current mean DM density is denoted by ρ_DM, and dn/dM is the halo mass function (Sheth & Tormen 1999). The DM density profile within a halo of mass M and at redshift z is denoted as ρ(x|M, z), which is taken from the Navarro-Frenk-White profile (Navarro et al. 1997). The boost factor b_sub quantifies the enhancement of the halo emission attributable to the presence of subhaloes. We considered three scenarios with different concentration parameters for subhaloes: LOW (Sánchez-Conde & Prada 2014), MID (Moliné et al. 2017) and HIGH (Gao et al. 2012).

The process of DM decay involves the spontaneous transition of unstable DM particles into standard model particles, accompanied by the emission of gamma-rays. The window function for this decay process is described as follows (Ando & Komatsu 2006; Ibarra et al. 2013; Fornengo & Regis 2014; Tröster et al. 2017):

$$\begin{array}{c}\hfill W_{{\mathrm{g}}_{\mathrm{dec}}}(z,E)=\frac{{\mathrm{\Omega }}_{\mathrm{DM}}{\rho }_c}{4\pi }\frac{{\mathrm{\Gamma }}_{\mathrm{dec}}}{m_{\mathrm{DM}}}\frac{\mathrm{d}{N}_{\mathrm{dec}}}{\mathrm{d}E}[E(1+z)]{\mathrm{e}}^{-\tau [z,E(1+z)]},\end{array}$$(8)

where Γ_dec is the decay rate; dN_dec/dE represents the energy spectrum of the DM decay process. It is analogous to DM annihilation at twice the energy (Cirelli et al. 2011), which can be expressed as dN_dec/dE(E) = dN_ann/dE(2E). The window functions for both the annihilation and the decay of DM are shown in Fig. 2.

Fig. 2. Window functions for gamma-rays produced by DM annihilation, decay, and astrophysical sources across an energy range of 0.5 − 1000 GeV. The annihilation process assumes ⟨σannv⟩ = 3 × 10−26 cm3 s−1 (Steigman et al. 2012), where the three scenarios are considered: HIGH (purple), MID (orange), LOW (green), representing different amounts of DM concentration for subhaloes. The DM decay (black) assumes a particle decay rate of Γdec = 5 × 10−28 s−1, taken as a representative benchmark consistent with current limits. Both processes are modelled with the assumption of the final state in $b\overline{b}$ pairs and mDM = 100 GeV for this figure. The contribution from unresolved blazars, which were taken as the only astrophysical source of the UGRB in this analysis, is shown in pink.

2.3. Gamma-rays from astrophysical sources

The UGRB is dominated by unresolved astrophysical sources, which constitute the main source of contamination when probing DM signals. The astrophysical gamma-ray background is generally thought to be dominated by contributions from blazars, mAGN, and SFGs. Compared to mAGN and SFGs, blazars are less numerous but produce a higher level of spatial anisotropy at the current level of sensitivity. Consequently, blazars are expected to dominate the angular power spectrum of the UGRB (Cuoco et al. 2012; Di Mauro et al. 2014). The window function from the astrophysical sources is shown as the pink line in Fig. 2.

In this analysis, we adopted the assumptions from Korsmeier et al. (2022), considering blazars as the only contributors to the unresolved astrophysical gamma-ray background. Specifically, we used an averaged luminosity-dependent density evolution (LDDE) model based on two populations of blazars: BL Lacertae objects and Flat Spectrum Radio Quasars, which are more numerous in the 4GFL catalogue (Ajello et al. 2015). The window function for astrophysical sources is modelled in the general form:

$$\begin{array}{c}\hfill W_{{\mathrm{g}}_{\mathrm{S}}}(z,E)=\chi {(z)}^2⟨f_{\mathrm{S}}⟩,\end{array}$$(9)

where ⟨f_S⟩ is the mean flux. For blazars, it can be expressed as

$$\begin{array}{cc}\hfill ⟨f_{\mathrm{S}}⟩& ={\displaystyle {\int }_1^{3.5}}\mathrm{d}\mathrm{\Gamma }{\displaystyle {\int }_{L_{\min }}^{L_{\max }}}\mathrm{d}L\mathrm{\Phi }(L,z,\mathrm{\Gamma })\times \frac{\mathrm{d}F}{\mathrm{d}E}(L,z,\mathrm{\Gamma })(1-\mathrm{\Omega }(S,\mathrm{\Gamma })),\hfill \end{array}$$(10)

where Γ is the photon spectral index, where the integral range follows the assumption in Ajello et al. (2015); L is the rest-frame gamma-ray luminosity in the energy range [0.1, 100] GeV; the luminosity bounds are L_min = 10⁴³ erg/s and L_max = 10⁵² erg/s; Φ is the gamma-ray luminosity function (GLF), which provides the number density of sources per unit luminosity and per unit co-moving volume at redshift z and photon spectral index Γ; dF/dE is the spectral energy distribution (SED); the term Ω refers to the Fermi-LAT sensitivity to detect sources.

The GLF term taken from Ajello et al. (2015) can be decomposed into its expression at z = 0 and a redshift-evolution term in the form Φ(L, z, Γ) = Φ(L, z = 0, Γ) e(z, L), where

$$\begin{array}{cc}\hfill \mathrm{\Phi }(L,z=0,\mathrm{\Gamma })& =\frac{A}{ln(10)L_{\gamma }}{[{\left(\frac{L_{\gamma }}{L_{\ast }}\right)}^{{\gamma }_1}+{\left(\frac{L_{\gamma }}{L_{\ast }}\right)}^{{\gamma }_2}]}^{-1}{\mathrm{e}}^{-0.5{[\mathrm{\Gamma }-\mu (L_{\gamma })]}^2/{\sigma }^2}.\hfill \end{array}$$(11)

The parameters are taken from the LDDE best-fit values presented in Table 1 of Ajello et al. (2015).

The SED term is modelled as a power law:

$$\begin{array}{c}\hfill \frac{\mathrm{d}F}{\mathrm{d}E}={\displaystyle \int d}E\frac{\mathrm{d}{N}_{\mathrm{S}}}{\mathrm{d}E}{\mathrm{e}}^{-\tau (E,z)},\end{array}$$(12)

where the integral range is [0.1, 100] GeV; the term e^{−τ(E, z)} refers to the attenuation of photons by the EBL; and dN_S/dE is the observed energy spectrum, which can be described as (Ajello et al. 2015)

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\mathrm{S}}}{\mathrm{d}E}=K{[{\left(\frac{E}{E_{\mathrm{b}}}\right)}^{1.7}+{\left(\frac{E}{E_{\mathrm{b}}}\right)}^{2.6}]}^{-1}.\end{array}$$(13)

We note that the GLF term was defined with luminosities in the rest frame energy range of [0.1, 100] GeV. When predicting the unresolved gamma-ray flux, we extrapolated the corresponding source spectra to our analysis range of [0.5, 1000] GeV, including attenuation by the EBL, in order to match the gamma-ray energy bins used in the cross-correlation. Similar methods were applied to the flux regime, where GLF models calibrated on resolved Fermi-LAT sources were extrapolated to the unresolved population (e.g. Manconi et al. 2020; Korsmeier et al. 2022).

From simulations, the relation between break energy E_b and the LAT-measured power-law photon index can be approximated as log E_b(GeV)≈9.25 − 4.11Γ. The normalisation factor K is given by (Zhang et al. 2022)

$$\begin{array}{c}\hfill K=L{\left[4\pi d_L^{2}k{\displaystyle \int \mathrm{d}}EE\frac{\mathrm{d}{N}_{\mathrm{S}}}{\mathrm{d}E}(K=1)\right]}^{-1},\end{array}$$(14)

where k is the K-correction term, obtained as (Wang et al. 2015)

$$\begin{array}{c}\hfill k=\left[{\displaystyle {\int }_{E_{\min }/(1+z)}^{E_{\max }/(1+z)}}\mathrm{d}EE\frac{\mathrm{d}{N}_{\mathrm{S}}}{\mathrm{d}E}\right]{\left[{\displaystyle {\int }_{E_{\min }}^{E_{\max }}}\mathrm{d}EE\frac{\mathrm{d}{N}_{\mathrm{S}}}{\mathrm{d}E}\right]}^{-1}.\end{array}$$(15)

The term Ω in Eq. (10) refers to the Fermi-LAT sensitivity and is modelled as a step function (Manconi et al. 2020)

$$\begin{array}{c}\hfill \mathrm{\Omega }(\mathrm{\Gamma })=\mathrm{\Theta }[S_{E,100}-S_{\mathrm{thr}}(\mathrm{\Gamma })],\end{array}$$(16)

where $S_{E,100}={\displaystyle {\int }_{0.1\mathrm{GeV}}^{100\mathrm{GeV}}}\mathrm{d}EE\frac{\mathrm{d}F}{\mathrm{d}E}$. The sources with a given spectral index Γ are regarded as undetected (Ω = 0) if their flux falls below the effective detection threshold S_thr(Γ). We derive S_thr(Γ) from the fourth Fermi-LAT source catalogue (4FGL, Abdollahi et al. 2020) following Section. B.1 of Manconi et al. (2020), where we bin 4FGL sources in Γ and define S_thr in each bin such that 98% of the catalogue sources have S_E, 100 > S_thr(Γ), and then interpolate in Γ to obtain a smooth threshold function.

2.4. Halo model

The 3D power spectrum can be effectively described using the DM halo model (see Asgari et al. 2023, and references therein). Within this formalism, the power spectrum, P_uv(k), is decomposed into two components: the two-halo term, representing the correlation between different haloes, and the one-halo term, representing the correlations within the same halo:

$$\begin{array}{c}\hfill P_{\mathit{uv}}(k)=P_{\mathit{uv}}^{2\mathrm{h}}(k)+P_{\mathit{uv}}^{1\mathrm{h}}(k).\end{array}$$(17)

The general 3D cross-power spectrum of one- and two-halo terms, related to the fields u and v is given by

$$\begin{array}{cc}\hfill P_{\mathit{uv}}^{1\mathrm{h}}(k)& ={\displaystyle {\int }_{M_{\min }}^{M_{\max }}}\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}(M)⟨\widehat{u}(k|M)\widehat{v}(k|M)⟩,\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill P_{\mathit{uv}}^{2\mathrm{h}}(k)& =⟨b\widehat{u}⟩(k)⟨b\widehat{v}⟩(k)P^{\mathrm{lin}}(k),\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill ⟨b\widehat{u}⟩(k)& ={\displaystyle {\int }_{M_{\min }}^{M_{\max }}}\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}{b}_{\mathrm{h}}(M)⟨\widehat{u}(k|M)⟩,\hfill \end{array}$$(20)

where P^lin is the linear power spectrum; b_h is the linear bias; $\widehat{u}(k|M)$ and $\widehat{v}(k|M)$ represent the Fourier transforms of the halo profile for u and v with mass M.

We calculated the cross angular power spectrum of convergence and gamma-rays emitted from DM with the Core Cosmology Library (CCL, Chisari et al. 2019). The matter power spectrum was estimated via the CAMB library (Lewis et al. 2000), where the non-linear part of the spectrum was fitted using halofit (Smith et al. 2003; Takahashi et al. 2012). The prediction of the angular power spectrum without the Fermi-LAT PSF correction (as discussed in Sect. 3.2) is demonstrated in Fig. 3. We note that the modelling uncertainty in the transition between the 1- and 2-halo transition regime may lead to small shifts in the predicted cross-correlation amplitude and therefore to modest changes in the inferred upper limits.

Fig. 3. Model of $C_{\ell }^{\mathrm{g}\kappa }$ for three clumping cases of DM annihilation, decaying DM, and astrophysical background. The models and formatting are the same as in Fig. 2. The models assumed the redshift n(z) for z in the range of [0.1, 2.0] and energy range from 0.5 to 1000 GeV. The effect of the Fermi-LAT PSF was not considered. The blazar component shown in pink represents the unresolved blazar contribution to the astrophysical background adopted in this analysis.

For the analysis of cross-correlating the cosmic shear and unresolved astrophysical sources, we adopted the one- and two-halo models with the window function detailed in Sect. 2.3. The 3D cross-power spectrum of shear and astrophysical sources can be expressed as

$$\begin{array}{cc}\hfill P_{\mathrm{g}\mathrm{S}\kappa }^{1\mathrm{h}}(k,z)& ={\displaystyle {\int }_{L_{\min }}^{L_{\max }}}\mathrm{d}L\frac{\mathrm{\Phi }(L,z)}{⟨f_{\mathrm{S}}⟩}\frac{\mathrm{d}F}{\mathrm{d}E}(L,z){\widehat{u}}_{\kappa }(k|M(L,z),z),\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill P_{\mathrm{g}\mathrm{S}\kappa }^{2\mathrm{h}}(k,z)& =\left[{\displaystyle {\int }_{L_{\min }}^{L_{\max }}}\mathrm{d}L{b}_{\mathrm{S}}(L,z)\frac{\mathrm{\Phi }(L,z)}{⟨f_{\mathrm{S}}⟩}\frac{\mathrm{d}F}{\mathrm{d}E}(L,z)\right]\hfill \\ \hfill & \times \left[{\displaystyle {\int }_{M_{\min }}^{M_{\max }}}\mathrm{d}M\frac{\mathrm{d}{n}_{\mathrm{h}}}{\mathrm{d}M}{b}_{\mathrm{h}}(M,z){\widehat{u}}_{\kappa }(k|M,z)\right]\hfill \\ \hfill & \times P^{\mathrm{lin}}(k,z),\hfill \end{array}$$(22)

where b_S is the bias of gamma-ray astrophysical sources relative to the matter distribution, which can be computed in terms of the halo bias b_S(L, z) = b_h(M(L, z)) as assumed in Camera et al. (2015), and the linear bias b_h is taken from the model of Tinker et al. (2010); Φ is the GLF of blazars; dF/dE is the differential photon flux; and ⟨f_S⟩ is the mean flux of blazars.

We modelled the angular cross-power spectrum of astrophysical sources from the UGRB and weak lensing signals also using CCL. For the transformation of the gamma-ray luminosity of a blazar to the mass of its host DM halo, we adopted the ‘Model B’, which corresponds to an intermediate scaling between halo mass and blazar luminosity, following the relation introduced by Camera et al. (2013):

$$\begin{array}{c}\hfill M={10}^{13}{M}_{\odot }{\left(\frac{L}{{10}^{44}{\mathrm{erg}\mathrm{s}}^{-1}}\right)}^{1.7}.\end{array}$$(23)

3. Data

3.1. KiDS

We used the fifth data release from the Kilo-Degree Survey (KiDS-Legacy; Wright et al. 2024) for the weak gravitational lensing data sample. KiDS is a wide-field optical imaging survey conducted using the OmegaCAM camera on the VLT Survey Telescope (VST; Capaccioli & Schipani 2011), specifically designed for weak lensing applications. The footprint of the KiDS can be divided into two patches: KiDS-N, covering regions along the equator, and KiDS-S, covering the southern sky. In conjunction with the infrared data from the VISTA Kilo-degree INfrared Galaxy survey (VIKING; Edge et al. 2013), the KiDS-Legacy footprint spans 1347 deg², with an effective area of 967.4 deg² after masking. The KiDS-Legacy dataset provides photometry in nine optical and near-infrared bands: ugriZYJHK_s, along with a second i-band pass and extra 23 deg² of KiDS-like calibration observations of deep spectroscopic surveys (Wright et al. 2025b), allowing for reliable photometric redshift estimation and redshift distribution calibration in KiDS-Legacy.

The galaxy shape measurements in KiDS-Legacy were performed with lensfit (Miller et al. 2007, 2013) and were calibrated with SKiLLS simulations as detailed in Li et al. (2023), which demonstrated that residual shear-related systematic errors are negligible after calibration. The robustness of the KiDS shear measurements were further validated by an independent METACALIBRATION-based analysis (Yoon et al. 2025). Our study selected galaxies from the gold-selected sample of KiDS-Legacy, which offers a high-accuracy redshift distribution. We conducted a tomographic cross-correlation analysis by dividing the sample into six photometric redshift bins: (0.10, 0.42], (0.42, 0.58], (0.58, 0.71], (0.71, 0.90], (0.90, 1.14] and (1.14, 2.00], based on the Bayesian photometric redshift z_B inferred from the BPZ algorithm (Benitez 2000). The redshift distributions n(z) were estimated and calibrated via a combination of self-organising map colour-based direct calibration, and cross-correlation methods using nine-band KiDS+VIKING photometry, and deep spectroscopic reference samples as detailed in Wright et al. (2025a). We summarise the shape noise and effective galaxy number density of KiDS-Legacy for each tomographic bin in Table 1.

3.2. Fermi-LAT

We analysed 15 years of observational data from the Fermi-LAT, spanning up to 3 July 2023. The Fermi-LAT has been surveying the gamma-ray sky since 2008 (Atwood et al. 2009). Its mission is to conduct long-term, high-sensitivity observations of celestial sources, covering an energy range from ∼20 MeV to > 300 GeV. The LAT is a wide field-of-view imaging gamma-ray telescope with a large effective area, combined with excellent energy and angular resolution.

In this study, we used the Pass 8 event data and processed the data preparation with FERMI SCIENCE TOOLS version 2.2.0¹. Our selection criteria included ultracleanveto photons and applied cuts on the data for higher quality with filters DATA_QUAL>0 && LAT_CONFIG==1. Photons with the lowest PSF quartile (PSF0) were excluded, and the surviving events were binned into 100 logarithmically spaced bins between 0.5 and 1000 GeV. The generated photon counts and exposure maps were projected in HEALP IX² (Górski et al. 2005) format with nside=1024, corresponding to a pixel size of approximately 3.4 arcmin. Intensity maps were obtained by dividing the counts by the exposure and were subsequently co-added into the analysis bins as described below. The process is facilitated by the healpy package (Zonca et al. 2019).

To isolate the UGRB, we applied masks for both diffuse Galactic emission and point sources. We masked sky regions where the fiducial Galactic diffuse emission model (gll_iem_v07³) exceeds three times the intensity of the isotropic template. All gamma-ray sources from the 4FGL-DR3 catalogue⁴ were then masked. The mask size was set according to the LAT point spread function (PSF), which describes the energy-dependent smearing of photon directions due to the finite angular resolution of the instrument. We analysed the PSF in each energy bin using the P8R2_SOURCE_V6 instrument response functions. In harmonic space, the real-space PSF is represented by the beam window b_ℓ, which attenuates the observed angular power spectrum, particularly at high multipoles (see Fig. 4). For each energy bin, we set the source mask radius to be twice the LAT 68% containment radius, derived from the corresponding PSF parameterisation of the response functions. To assess the effectiveness of the point source masking we adopted, we estimated the fraction of point source flux leakage. Even in the lowest-energy bin (0.5–1 GeV), which suffered most from the PSF; the leakage fraction is ≲3%, indicating that residual contamination from resolved sources was negligible for our analysis.

Fig. 4. Fermi PSF in harmonic space bℓ in ten different energy bins. It shows an inverse relationship between the energy and the suppression of the PSF effect. The dashed black line represents ℓ = 1500, referring to the upper limit of ℓ in the cross-power spectrum measurement.

After applying the masks, we refitted the normalisations of the Galactic diffuse and isotropic templates in each of the 100 fine energy bins by maximising a Poisson likelihood. The residual photon counts were then corrected for the LAT exposure and combined into ten energy bins covering [0.5, 1000] GeV, with edges [0.5, 1.0, 1.99, 3.97, 7.93, 15.83, 31.59, 63.04, 125.81, 251.08, 1000] GeV, which were used in the subsequent cross-correlation analysis. We summarise the photon counts in the unmasked region for each energy bin in Table 2.

4. Methods

In this section, we describe the methodologies we employed to calculate the cross-correlation between the galaxy lensing data and the gamma-ray intensity map, to estimate the covariance, and to constrain parameters related to DM.

4.1. Power spectrum estimation

In this study, we created maps and masks for both KiDS-Legacy and Fermi-LAT data to estimate the cross-correlation angular power spectrum of weak lensing and gamma-ray intensity. We used the MASTER algorithm (Hivon et al. 2002) implemented in NaMaster⁵ (Alonso et al. 2019) to compute the cross power spectra. NaMaster provides a general framework to estimate pseudo-C_ℓ on masked fields with arbitrary spin and contaminants, and to deconvolve mask-induced mode coupling via the coupling matrix.

For the weak lensing component of the KiDS-Legacy data, we created HEALP IX maps for different redshift bins with nside=1024 covering the footprint of the KiDS-Legacy gold-selected sample catalogue. We took the weighted average of the galaxy ellipticities for each pixel, where the weights were determined by the shape measurement uncertainties, forming the spin-2 field {γ₁, γ₂}. The corresponding masks, matching the coverage and resolution of the maps, were generated as a sum of shape measurement weights for each pixel.

For the Fermi-LAT gamma-ray data, we constructed the gamma-ray intensity maps for different energy bins as described in Sect. 3.2. The mask used for the gamma-ray data was binary (1 for the unmasked regime, and 0 for the masked regime) and energy-dependent, with point sources and Galactic diffuse emission removed.

Using the three fields, {g, γ₁, γ₂}, we applied NaMaster to estimate the cross-power spectra between the weak lensing E/B-modes and the gamma-ray intensity I_γ. The E-mode is the lensing-induced component and thus carries the cosmological shear signal, which can correlate with gamma-ray emission from astrophysical sources, as well as any potential DM contribution. The B-mode of the shear is obtained by rotating galaxy shapes by 45°. Since gravitational lensing is parity-conserving, pure gravitational lensing is expected to produce only E-mode patterns, with no B-modes in the absence of systematic effects or intrinsic alignments. The B-mode therefore provides a valuable null test and diagnostic of residual systematics. The E- and B-modes are analogous to the decomposition used in CMB polarisation studies (see, e.g. Kamionkowski et al. 1997), where E/B-modes represent the curl-free and divergence-free components, respectively.

The resultant cross-power spectrum was organised into five linearly spaced multipole bins spanning from 200 to 1500. The upper limit was set by the Fermi–LAT PSF, which suppresses power on small angular scales (high ℓ; see Fig. 4), whereas the lower limit was chosen considering the survey sky coverage and residuals from foreground subtraction on large angular scales.

4.2. Covariance

We estimated the covariance of the angular cross-power spectra between tomographic weak-lensing data and energy-binned gamma-ray intensity maps using two approaches: a weighted jackknife resampling, used as a data-driven consistency check, and a Gaussian covariance computed with NaMaster, which is adopted as the fiducial covariance in the analysis.

Jackknife resampling is a non-parametric technique that estimates covariance by systematically excluding one or more tiles of observations from the dataset and recalculating the estimate across all subsets. In our analysis, we employed the delete-one weighted jackknife resampling to estimate the covariance of the cross-correlation between weak lensing data and gamma-rays.

For each redshift–energy bin, we divided the corresponding footprint into 50 approximately equal-area sub-regions using the HEALP IX scheme. The k-th sub-region was removed from the sample, and we recalculated the pseudo-C_ℓ of the remaining parts for all realisations. Realisations were weighted by the sum of galaxy shape measurement weights in the retained area. The covariance of the measurement was computed as follows (Mohammad & Percival 2022):

$$\begin{array}{c}\hfill {\mathbf{C}}^{\mathrm{JK}}(C_{\ell }^{\mathrm{g}\kappa },C_{\ell \prime }^{\mathrm{g}\kappa })=(n_s-1)\frac{{\sum }_k^{n_s}{w}_k(C_{\ell ,k}^{\mathrm{g}\kappa }-⟨C_{\ell }^{\mathrm{g}\kappa }⟩)(C_{\ell \prime ,k}^{\mathrm{g}\kappa }-⟨C_{\ell \prime }^{\mathrm{g}\kappa }⟩)}{{\sum }_k{w}_k},\end{array}$$(24)

where n_s = 50 is the number of subsamples used in the jackknife resampling, w_k is the weight for the k-th patch, $C_{\ell ,k}^{\mathrm{g}\kappa }$ is the cross-power spectrum between convergence and gamma-rays for the k-th patch, and $⟨C_{\ell ,k}^{\mathrm{g}\kappa }⟩$ is the weighted average of the cross-angular power spectrum for all jackknife subsamples. Because the gamma-ray masks were energy dependent, the jackknife procedure was applied independently to each redshift–energy bin, and the resulting estimates were used only for the corresponding diagonal covariance blocks. However, we note that the jackknife subsamples do not carry strictly equal statistical weight, since the effective area is modulated by the energy-dependent Fermi-LAT mask, particularly at low energies where the PSF-driven point-source masking is most complex.

Gaussian covariance estimation with NaMaster was employed under the Gaussian-field assumption and accounted for the mode-coupling induced by the survey masks through coupling matrices. In practice, for each pair of fields {g, κ}, we first measured the coupled pseudo-C_ℓ ${\stackrel{\sim }{C}}_{\ell }^{\mathrm{g}\kappa }$ from the masked maps and normalised them by the effective overlap of the corresponding masks:

$$\begin{array}{c}\hfill C_{\ell }^{\mathrm{g}\kappa ,\mathrm{Cov}}\simeq \frac{{\widehat{C}}_{\ell }^{\mathrm{g}\kappa }}{⟨w_{\mathrm{g}}{w}_{\kappa }⟩},\end{array}$$(25)

where $C_{\ell }^{\mathrm{g}\kappa ,\mathrm{Cov}}$ is the power spectrum between the gamma-ray field g and the E-mode shear (convergence) field κ used to estimate the covariance matrix. The term ${\stackrel{\sim }{C}}_{\ell }^{\mathrm{g}\kappa }$ is the corresponding coupled pseudo-C_ℓ. The gamma-ray masks w_γ are binary maps, and the convergence masks w_κ correspond to the weighting maps of the galaxy shape measurements. The factor ⟨w_gw_κ⟩ represents the average product of the masks over all pixels.

We assessed the consistency of the NaMaster Gaussian covariance by comparing its diagonal terms with jackknife estimates across all tomographic and energy bins, and found a good agreement. Given that the NaMaster covariance naturally accounts for mask-induced mode coupling and provides a self-consistent description of both variances and cross-bin correlations, we adopt it as the fiducial covariance in the likelihood analysis. The jackknife estimates are used as an independent validation of the covariance amplitude on the diagonal. We show the correlation matrix corresponding to the fiducial covariance in Appendix A.

For the forecast covariance, we assumed Gaussian-distributed fields and accounted for mask-induced mode-coupling using NaMaster. In this case, we modelled power spectra ${\widehat{C}}_{\ell }^{\kappa \kappa \prime }$, ${\widehat{C}}_{\ell }^{\mathrm{gg}\prime }$ and ${\widehat{C}}_{\ell }^{\mathrm{g}\kappa }$ as inputs, where ${\widehat{C}}_{\ell }^{\kappa \kappa \prime }$ and ${\widehat{C}}_{\ell }^{\mathrm{gg}\prime }$ denote the estimated auto power spectra of convergence κ and gamma-ray g, respectively. The cross-spectrum ${\widehat{C}}_{\ell }^{\mathrm{g}\kappa }$, defined in Sect. 2, includes the DM-induced contribution (from annihilation or decay, depending on the scenario considered) but does not include the blazar component.

The lensing auto power spectrum ${\widehat{C}}_{\ell }^{\kappa \kappa \prime }$ was modelled as the combination of the cosmic shear signal and shape noise:

$$\begin{array}{c}\hfill {\widehat{C}}_{\ell }^{\kappa \kappa \prime }=C_{\ell }^{\kappa \kappa \prime }+{\delta }_{\kappa \kappa \prime }\frac{{\sigma }_e^2}{n_{\mathrm{eff}}},\end{array}$$(26)

where C_ℓ^κκ′ is the cosmic shear signal computed for each tomographic bin; σ_e² and n_eff are the dispersion of the ellipticities, and effective galaxy number density of the survey, which only contribute to the auto spectra terms.

The gamma-ray power spectrum ${\widehat{C}}_{\ell }^{\mathrm{gg}\prime }$ was estimated from the gamma-ray flux maps used in the cross-correlation. We measured the auto- and cross-angular spectra among the ten energy bins using NaMaster in 15 logarithmically spaced ℓ-bins between 30 and 2000. We fitted the measured spectra with the form

$$\begin{array}{c}\hfill {\widehat{C}}_{\ell }^{\mathrm{gg}}=C_{\mathrm{Poisson}}+c{\ell }^{\alpha },\end{array}$$(27)

where the Poisson term C_Poisson captures the photon shot noise, and c and α account for a power-law contribution for large scales. For the auto-spectra, we found that the Poisson term dominates and the best-fitting amplitude c was consistent with zero, while for cross-spectra the Poisson term vanished and only the power-law component was retained. A more detailed comparison between the different covariance prescriptions, including Jackknife, NaMaster covariance constructed from the measured power spectra, and the NaMaster covariance based on analytical power spectra, is presented in Appendix B.

4.3. Likelihood

We constrained the model parameters p = {⟨σ_annv⟩,Γ_dec} for each DM mass bin as described in Sect. 2, assuming the measured power spectra follow a Gaussian likelihood:

$$\begin{array}{c}\hfill {\chi }^2=-2lnL(\mathit{d}|\mathit{p})=[\mathit{d}-\mathit{\mu }(\mathit{p}){]}^{\mathrm{T}}{\mathbf{C}}^{-1}[\mathit{d}-\mathit{\mu }(\mathit{p})],\end{array}$$(28)

where d is the data vector, comprising 6 × 10 cross-correlations with six redshift bins in weak lensing data and ten energy bins in gamma-ray data; μ(p) is the theoretical cross-correlation predicted by the model parameters; and C⁻¹ is the inverse of the covariance matrix. Uniform priors are assigned for the model parameters p.

Given that our data lack the constraining power to fully determine the DM parameters, we derived the 95% confidence upper limit constraints on ⟨σ_annv⟩ and Γ_dec from the contours of the likelihood surface. We followed the assumption from Tröster et al. (2017), and for a specified confidence interval p_c, the contours were defined by the parameter sets d, where

$$\begin{array}{c}\hfill {\chi }^2(\mathit{d},\mathit{p})={\chi }^2(\mathit{d},{\mathit{p}}_{\mathrm{ML}})+\mathrm{\Delta }{\chi }^2(p_{\mathrm{c}}),\end{array}$$(29)

where p_ML is the maximum likelihood estimation of the parameters, χ² is calculated from Eq. (28), and the Δχ²(p_c) corresponds to the quantile function of the χ² distribution. For 2σ contours, Δ_χ²(0.95) = 6.18. We calculated the limits on ⟨σ_annv⟩ and Γ_dec for different final states of standard model particles with m_DM in logarithmically spaced values in the range [10,1000] GeV.

5. Results

5.1. Cross-correlation measurement

We present our measurements of the cross-correlation angular power spectrum between the Fermi-LAT gamma-ray intensity maps (across ten energy bins) and KiDS-Legacy weak lensing data (across six redshift bins), utilising a combined covariance, as detailed in Sect. 4.2. We show the results in Fig. 5, with black and red points representing the cross-correlation of gamma-rays with the components of the E- and B-modes of shear, respectively.

Fig. 5. Cross angular power spectrum Cℓgκ between the Fermi-LAT gamma-ray intensity map and KiDS-Legacy weak lensing data when using a NaMaster covariance. E- and B-mode cross-spectra are shown as black and red points, respectively. Panels are labelled by the tomographic redshift bin zi and gamma-ray energy bin Ej. The reduced χν2 (ν = 5) with respect to null signal for the E-mode in each bin is also provided in each subplot.

The χ_ν² values with respect to the null signal with degrees of freedom ν = 5 are illustrated in subplots for all bins. The global χ² values for the E- and B-mode are 248.07 and 256.28, respectively, evaluated over 300 data points. These correspond to p-values of 0.987 for the E-mode and 0.956 for the B-mode, indicating that both are consistent with the null hypothesis and that no significant cross-correlation signal is detected. We note that this χ² is below the degrees of freedom, suggesting that the covariance could be somewhat overestimated. As discussed in Appendix B, the jackknife, NaMaster–measured, and NaMaster–analytic covariance estimates yield consistent diagonal variances, indicating internal consistency among the different estimation methods. However, all three approaches assume isotropic and spatially uniform noise, which is not strictly valid for the Fermi–LAT intensity maps. Because the gamma-ray maps are dominated by photon shot noise that depends on the local exposure time, and since the exposure varies across the sky, the effective noise level is spatially inhomogeneous. This can lead to mildly overestimated variances and hence a lower global χ². We tested this by splitting the KiDS footprint into North and South patches. We removed the area of KiDS–N with the strongest exposure variations to make both regions have comparable effective sky coverage. As a side effect, this also makes KiDS–N more uniform in exposure than KiDS–S, which helps isolate the effect of inhomogeneous photon noise. The result shows that KiDS–S, which exhibits stronger Fermi exposure variations, yields a lower χ² (241.81 in KiDS–S and 301.37 in KiDS–N, among 300 data points). This is consistent with the interpretation that anisotropic photon noise leads to a conservative covariance estimate. Overall, these results are consistent with the non-detection of signals across all cross angular power spectrum measurements for the six redshift bins and ten energy bins.

5.2. Constraints on DM parameters

We applied the formalism detailed in Sect. 4.3 to derive constraints on the DM thermally averaged annihilation cross-section ⟨σ_annv⟩ and decay rate Γ_dec for three final states: $b\overline{b},{\mu }^{+}{\mu }^{-}$, and τ⁺τ⁻, each state with a branching ratio of 1. In the modelling, we accounted for the blazar contribution to the UGRB as described in Sect. 2.3. Although the background of astrophysical sources dominates the mixed UGRB, the improvements in constraints after excluding blazar effects remain limited. Due to the non-detection of cross-correlation signals, we can only provide 95% upper bounds, without establishing both upper and lower limits.

We present our results on the 95% upper bounds of ⟨σ_annv⟩ and Γ_dec with m_DM in the range of [10, 1000] GeV in Fig. 6 and Fig. 7, respectively. Compared with the analysis of Paopiamsap et al. (2024), our constraints are less stringent. Their study focuses on the cross-correlation between the UGRB and the galaxy overdensity at low redshift (z ≲ 0.4), using wide-area galaxy surveys covering a large fraction of the extragalactic sky, which yields a high S/N detection at the 8–10 σ level. In contrast, our study relies on shear measurements over a smaller effective coverage.

Fig. 6. Upper bounds at 95% confidence on the DM annihilation cross-section ⟨σannv⟩ for $b\overline{b},{\mu }^{-}{\mu }^{+}$, and τ−τ+ states making use of six redshift bins and ten energy bins. The unresolved blazar contribution to the UGRB is included as the astrophysical background component in the analysis. For reference, dashed lines indicate the limits obtained without subtracting the blazar contribution. The thermal relic cross-section for WIMPs (Steigman et al. 2012) is shown as dashed-dot lines. The upper bounds are presented for three DM clumping scenarios: LOW (green), MID (orange) and HIGH (purple).

Fig. 7. Upper bounds at 95% confidence on the DM decay rate Γdec for $b\overline{b},{\mu }^{-}{\mu }^{+}$, and τ−τ+ states include the unresolved blazar contribution to the UGRB as the astrophysical background (solid curves). The dashed curves show the limits obtained without the blazar component.

Our results are simplified as we do not account for gamma-ray emission from inverse Compton scattering of CMB photons. Compared to previous cross-correlation analyses based on CFHTLenS, RCSLenS, and KiDS data (Tröster et al. 2017), our constraints are generally tighter for LOW and MID clumping scenarios, while differences at HIGH clumping are likely driven by modelling choices such as substructure prescriptions and the treatment of additional emission components. For heavy DM candidates, inverse Compton scattering can significantly contribute to the gamma-ray signal at high energies (Ando & Ishiwata 2016). In addition, we adopt a simplified model of the astrophysical background. Incorporating more complete modelling of astrophysical sources and additional emission processes is therefore expected to yield more accurate and potentially tighter constraints in future analyses.

We also compared the 95% upper bounds of ⟨σ_annv⟩ and Γ_dec from large-scale cosmological cross-correlation (for the τ⁺τ⁻ state from this work and Paopiamsap et al. 2024) with constraints derived from local structures as shown in Fig. 8. Local probes such as IceCube (neutrinos), H.E.S.S., and MAGIC (gamma-rays) focus on regions of high DM density, such as the Galactic Centre and dwarf spheroidal galaxies, where annihilation or decay signals are expected to be largest. By contrast, cosmological analyses constrain DM by averaging over a much larger volume and are less tied to assumptions about individual targets and substructure. As shown in Fig. 8, cosmological methods provide competitive constraints on low m_DM for both ⟨σ_annv⟩ and Γ_dec.

Fig. 8. Left: Upper bounds at 95% confidence on the DM thermally averaged cross-section ⟨σannv⟩ of HIGH clumping case for the τ+τ− state, compared to constraints from local probes. This includes results from neutrino detectors such as IceCube (Abbasi et al. 2023) and combined results from ANTARES and IceCube (Albert et al. 2020); the gamma-ray telescope H.E.S.S. (Abdallah et al. 2016), combined results from Fermi-LAT and MAGIC (MAGIC Collaboration 2016); and the red dotted lines represent thermal relic cross-section. Right: Upper bounds on the DM decay rate Γdec for the τ+τ− state, compared to limits from IceCube (Aartsen et al. 2018) for neutrino surveys, HAWC (Albert et al. 2018) and searches of cosmic rays and the interstellar medium with Fermi-LAT (Ackermann et al. 2012) from gamma-ray telescopes. The cosmological constraints from this work (black solid lines) and work from Paopiamsap et al. (2024) (black dashed lines) are also shown. GC: Galactic centre; dSph: dwarf spheroidal galaxies; Cosmo: cosmological probe.

5.3. Forecast of DM parameters with Euclid

Future surveys aiming to measure cosmic shear, covering a wider area and more extensive depth, such as LSST and Euclid, are expected to offer improved sensitivity to DM signals. In this section, we present a forecast of constraints on the DM annihilation cross-section and decay rate using a Euclid-like survey. The Euclid survey observes about 15 000 deg² of the extragalactic sky. For this analysis, we adopted ten tomographic bins, together with the survey area, redshift distributions, intrinsic ellipticity dispersion, and expected source number densities specified in Table 4 of Euclid Collaboration: Blanchard et al. (2020).

We employed a Fisher matrix formalism to forecast the uncertainties of DM parameters ⟨σ_annv⟩ and Γ_dec as a function of the logarithmically spaced DM mass m_DM. The Fisher matrix with Gaussian likelihood is defined as

$$\begin{array}{c}\hfill {\mathbf{F}}_{\mathit{ij}}=\langle \frac{{\partial }^{2}L(\mathit{d}|\mathit{p})}{\partial p_i\partial p_j}\rangle ={\displaystyle \sum _{\ell ,\ell \prime }}\frac{\partial {\mathit{\mu }}_{\ell }(\mathit{p})}{\partial p_i}{\mathbf{C}}_{\ell \ell \prime }^{-1}\frac{\partial {\mathit{\mu }}_{\ell \prime }(\mathit{p})}{\partial p_j},\end{array}$$(30)

where L(d | p) is the likelihood of observing the data d given the model parameters p, C_ℓℓ′ is the covariance, and ${\mathit{\mu }}_{\ell }(\mathit{p})$ is the theoretical prediction for the cross power spectrum, which represents the stacked vector containing all tomographic–energy combinations at ℓ. The marginalised uncertainties on the parameters p_i are given by

$$\begin{array}{c}\hfill \sigma (p_i)=\sqrt{{({\mathbf{F}}^{-1})}_{\mathit{ii}}}.\end{array}$$(31)

We used the same Fermi–LAT gamma-ray intensity maps as in the measurements. The covariance was calculated using the NaMaster-based analytical estimation detailed in Sect. 4.2, and we compared the forecasts using different covariance prescriptions in Appendix C, given that we employed the combined covariance for the measurements. The redshift distributions for ten equally populated tomographic bins with edges {0.001, 0.418, 0.560, 0.678, 0.789, 0.900, 1.019, 1.155, 1.324, 1.576, 2.50} were modelled as in Sect. 3.3.1 in Euclid Collaboration: Blanchard et al. (2020). We roughly estimated the effective number density for each redshift bin to be the galaxy number density divided by the number of redshift bins, n_gal/N_z, which was 3 arcmin⁻² per bin. The shape noise was set to σ_ϵ = 0.3, the total intrinsic ellipticity dispersion for Euclid, where σ_ϵ² = 2σ_e² relates it to the per-component ellipticity dispersion. We adopted the same multipole range as in the KiDS-Legacy measurement (ℓ_max = 1500), which is considered as the pessimistic choice of ℓ_max (Euclid Collaboration: Blanchard et al. 2020). The Fisher forecast for the cross-correlation between gamma-rays from Fermi-LAT and weak lensing data from Euclid-like survey is shown in Fig. 9, presenting the 95% bounds of ⟨σ_annv⟩ and Γ_dec for three different final state particles, where we assume a fiducial DM signal corresponding to ⟨σ_annv⟩ = 3 × 10⁻²⁶ cm³ s⁻¹ for annihilation and Γ_dec = 5 × 10⁻²⁸ s⁻¹ in the decay case. The exclusion of the thermal relic cross-section in the $b\overline{b}$ and τ⁺τ⁻ HIGH scenario should not be interpreted as a detection, but rather as a result of the small forecasted covariance in this case. Since these limits are derived within a Fisher forecast framework, they only reflect the statistical sensitivity of the experiment under the assumed model. Relative to KiDS-Legacy, the Euclid-like forecasts tighten the limits by factor of 2–4 for both ⟨σ_annv⟩ and Γ_dec depending on the final states and DM masses. This should be regarded as a conservative estimate, because the sensitivity of the cross-correlation does not scale solely with the weak-lensing data, but is also limited by the properties of the gamma-ray maps. At low energies, the large Fermi–LAT PSF requires larger point-source masking, which reduce the effective sky overlap, while at high energies the limited photon counts lead to large statistical uncertainties. As a result, the improvement in Euclid coverage does not propagate into a proportional tightening of the final constraints. Within this framework, surveys such as LSST, which achieve comparable effective source number densities and sky coverage, are expected to yield constraints of a similar order.

Fig. 9. Upper: Forecast of 95% upper and lower bounds on the DM annihilation cross-section ⟨σannv⟩ for $b\overline{b},{\mu }^{-}{\mu }^{+}$, and τ−τ+ states from cross-correlation of ten energy bins from Fermi-LAT and ten redshift bins from Euclid-like survey. The expected value ⟨σannv⟩ = 3 × 10−26 cm3 s−1 is shown as dot-dashed lines for all cases. Lower: Forecast of 95% upper and lower bounds on the DM decay rate Γdec for $b\overline{b},{\mu }^{-}{\mu }^{+}$, and τ−τ+ states. The dot-dashed lines represent the expected value for the forecast with Γdec = 5 × 10−28 s−1.

6. Conclusions and discussions

In this study, we performed a cross-correlation analysis between the gamma-ray intensity map from 15 years of Fermi-LAT data and weak lensing data from KiDS-Legacy, covering 1347 deg² of the sky. Our aim was to constrain the DM annihilation cross-section ⟨σ_annv⟩ and the decay rate Γ_dec for WIMP masses between 10 GeV and 1 TeV for three different final states: $b\overline{b}$, μ⁺μ⁻, and τ⁺τ⁻. We measured the cross-correlation angular power spectrum using six tomographic redshift bins for the weak lensing data and ten energy bins for the gamma-ray data, with an energy range from 0.5–1000 GeV. We found no significant detection in the multipole range 200 ≤ ℓ < 1500.

We modelled the unresolved astrophysical background in the UGRB with a blazar-dominated component, in order to isolate any potential DM contribution. Our analysis demonstrates the 95% upper bounds for ⟨σ_annv⟩ and Γ_dec. Compared to previous findings, our constraints exhibit greater stringency for the LOW and MID clumping scenarios. In comparison to the upper bounds from local structure observations, such as neutrino searches by IceCube and gamma-ray observations by H.E.S.S., MAGIC and Fermi-LAT, the cosmological cross-correlation approach is competitive and provides complementary, large-scale constraints on DM properties, particularly in the low-mass region.

A notable aspect of our study is that we do not find a statistically significant UGRB–lensing cross-correlation in the KiDS-Legacy × Fermi-LAT analysis. By contrast, using galaxy clustering as the large-scale structure tracer, Paopiamsap et al. (2024) report an ∼8–10σ detection when cross-correlating with the UGRB. As a robustness check, our pipeline recovers a comparable signal (∼8σ), consistent with the higher S/N of density tracers relative to shear. For weak lensing, DES Y1 found a 5.3σ detection (Ammazzalorso et al. 2020), and DES Y3 found an 8.9σ detection that cross-correlates with Fermi-LAT (Thakore et al. 2025). We recovered a DES Y3 UGRB-shear signal at ∼6σ in real space. It corresponds to a null-hypothesis χ² test using a jackknife-based covariance within our pipeline. In addition, our analysis uses ULTRACLEANVETO Fermi-LAT events, while DES Y3 adopted CLEANVETO events, which can further affect the quoted significance. Applying the methodology to KiDS-Legacy weak lensing and 15-year Fermi-LAT UGRB yields no significant detection. This difference can be explained by the effective sky area and the sensitivity to low-ℓ modes after gamma-ray processing and masking. Our stricter LAT selections might reduce photon counts and sky fraction and tend to down-weight the lowest multipoles, where DES Y3 derived much of its significance. Applying the methodology to KiDS-Legacy weak lensing and 15-year Fermi-LAT UGRB yields no significant detection. We also tested alternative multipole binning, and our conclusions were not changed. A summary of the additional validation tests is provided in Appendix D. The more fragmented geometry of the KiDS–Fermi-LAT overlap may increase mask-induced mode coupling. This reduces the number of truly independent low-ℓ modes and, therefore, weakens our leverage on large-scale correlations. Differences in tomography and noise properties further modify the effective kernel overlap between lensing and the gamma-ray intensity field. Variations in the redshift distributions, the effective source number density, and shape noise may affect the expected cross-correlation amplitude and uncertainty. In addition, choices in estimators and modelling, such as the treatment of the Fermi-LAT beam function and the level of diffuse-foreground masking, set the effective weight assigned to different angular scales. Those settings that down weight the lowest multipoles would reduce sensitivity to large-scale correlations. Another contributing factor is that the global χ² values are slightly lower than expected. As discussed in Sect. 5.1, this likely reflects a conservative covariance estimation arising from the spatially varying photon noise in the Fermi intensity maps, which should be accounted for in future analyses with higher signal-to-noise data.

A cross-correlation detection between the UGRB and large-scale structure tracers does not necessarily imply a DM detection. Without sufficiently complete modelling of unresolved astrophysical emission and residual large-scale foregrounds, the measured correlation can be explained by conventional source populations. In this context, DES Y3 × Fermi (Thakore et al. 2025) interprets their detection mainly as information on the unresolved blazar contribution. Likewise, Paopiamsap et al. (2024) models the astrophysical component using UGRB-galaxy cross-correlations and reports DM constraints as upper limits, obtained by conservatively attributing all the measured signal to a possible DM contribution. Although galaxy clustering currently yields a generally higher S/N in UGRB cross-correlation measurements, weak lensing provides a complementary probe. Unlike galaxy tracers, lensing is sensitive to the total matter distribution and does not depend on assumptions about galaxy bias. This makes UGRB–lensing cross-correlations less sensitive to modelling uncertainties in the tracer population, though the statistical uncertainties are larger at present.

We also forecast the potential improvements that future surveys can bring to DM research. Future surveys, such as the LSST and Euclid, with larger coverage and better photometry, will provide higher S/N and more precise cross-correlation measurements. For a Euclid-like configuration, our forecasts indicate about ∼2–4 times tighter constraints on both ⟨σ_annv⟩ and Γ_dec relative to KiDS-Legacy. On the gamma-ray side, the ongoing observations of Fermi-LAT will increase the total exposure time and reduce photon shot noise, which in turn increases the S/N of the cross-correlation and tightens the DM constraints.

In our analysis, we adopted fixed spectra and modelled the UGRB with blazars only, estimating limits via a maximum-likelihood approach without marginalising over astrophysical-background parameters. Given the low S/N in our measurements and the relatively simple parameter space of the DM model from our assumption, this approach simplified the analysis while maintaining computational efficiency. Although this method does not fully account for potential uncertainties in the blazar power spectra, it is sufficient within the scope of our assumptions, as the coupling between the blazar contribution and the DM-related parameters is negligible. Furthermore, focusing on blazars also constitutes a conservative modelling choice for the astrophysical background, avoiding overfitting in the absence of a significant detection.

Future studies, particularly those with higher sensitivity or more significant signals, would benefit from a more detailed treatment of astrophysical contributions. This involves accounting for additional sources of astrophysical background, such as SFGs, mAGN, and millisecond pulsars, or adopting Bayesian-based techniques and marginalising over the parameters related to astrophysical contributions. These advancements will be critical for refining constraints and deepening our understanding of DM properties.

Acknowledgments

We gratefully acknowledge the anonymous referee for their valuable comments that helped us to improve the manuscript. We thank Bhashin Thakore, David Alonso, Stefano Camera, Marco Regis, Henk Hoekstra, Stefan Fröse, Paul-Simon Blomenkamp and Tristan Gradetzke for helpful discussions during the preparation of this manuscript. SZ acknowledges the support from the Deutsche Forschungsgemeinschaft (DFG) SFB1491. H. Hildebrandt is supported by a DFG Heisenberg grant (Hi 1495/5-1), the DFG Collaborative Research Center SFB1491, an ERC Consolidator Grant (No. 770935), and the DLR project 50QE2305. ZY, CH and BS acknowledges support from the Max Planck Society and the Alexander von Humboldt Foundation in the framework of the Max Planck-Humboldt Research Award endowed by the Federal Ministry of Education and Research. TT acknowledges funding from the Swiss National Science Foundation under the Ambizione project PZ00P2_193352. MA acknowledges the UK Science and Technology Facilities Council (STFC) under grant number ST/Y002652/1 and the Royal Society under grant numbers RGSR2222268 and ICAR1231094. DJB was supported by the Simons Collaboration on “Learning the Universe” and support was provided by Schmidt Sciences, LLC. MB is supported by the Polish National Science Center through grant no. 2020/38/E/ST9/00395. D.E. acknowledges the support from DFG Sonderforschungsbereich 1491 Project F5 and BMBF ErUM-Pro. CH acknowledges the UK Science and Technology Facilities Council (STFC) under grant ST/V000594/1. BJ acknowledges support by the ERC-selected UKRI Frontier Research Grant EP/Y03015X/1 and by STFC Consolidated Grant ST/V000780/1. LM acknowledges the financial contribution from the grant PRIN-MUR 2022 20227RNLY3 “The concordance cosmological model: stress-tests with galaxy clusters” supported by Next Generation EU and from the grant ASI n. 2024-10-HH.0 “Attività scientifiche per la missione Euclid – fase E”. DN acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant agreement No. 101053992). AP acknowledges the support from “la Caixa” Foundation (ID 100010434, code LCF/BQ/DI24/12070011) Funding for this work was partially provided by the “Center of Excellence Maria de Maeztu” award to the ICCUB CEX2024-001451-M funded by MICIU/AEI/10.13039/501100011033. RR is partially supported by an ERC Consolidator Grant (No. 770935). AHW is supported by the Deutsches Zentrum für Luft- und Raumfahrt (DLR), under project 50QE2305, made possible by the Bundesministerium für Wirtschaft und Klimaschutz, and acknowledges funding from the German Science Foundation DFG, via the Collaborative Research Center SFB1491 “Cosmic Interacting Matters – From Source to Signal”.

References

Appendix A: Visualisation of the covariance

In this appendix, we present a visualisation of the covariance structure used in the analysis. As described in Sect. 4.2, the fiducial covariance is given by the NaMaster Gaussian covariance, which accounts for mask-induced mode coupling and correlations between multipole, redshift, and energy bins. To illustrate the resulting correlations, we show the corresponding correlation matrix in Fig. A.1.

Fig. A.1. Correlation matrix for the cross-correlation between weak lensing data and gamma-ray intensity maps across different redshift and energy bins. The redshift bins (z1 – z6) are the six tomographic bins listed in Table 1. The energy bins (E1 – E10) span 0.5–1000 GeV from low to high energy.

Appendix B: Comparison of covariance prescriptions

The constraints presented in the main text rely on the NaMaster Gaussian covariance as the fiducial choice. Since different covariance prescriptions rely on different assumptions, it is useful to assess how the results depend on this choice. In this appendix, we compare three covariance estimates: the jackknife covariance, the NaMaster Gaussian covariance constructed from analytical power spectra, and the NaMaster covariance based on measured power spectra.

We compare the diagonal terms (variances) of the jackknife (green), the NaMaster-analytic (blue) and NaMaster-measured (black) estimates in Fig. B.1. Overall, the different approaches yield variances of comparable magnitude across most redshift and energy bins, indicating a broadly consistent description of the statistical uncertainties. However, at low energies, where the Fermi-LAT masks are most complex and the effective sky coverage is smallest because of the Fermi PSF, the jackknife estimates tend to be larger than those from NaMaster. This behaviour is expected, as the jackknife captures additional sources of variance related to spatial fluctuations and inhomogeneous masking that are not included in the Gaussian approximation.

Fig. B.1. Comparison of diagonal variances for the gamma–shear cross-power spectra. Each subplot corresponds to one tomographic redshift bin (columns, from z1 to z6) and one Fermi–LAT energy bin (rows, from e1 to e10). The diagonal elements Var(Cℓgκ) are shown in the five multipole bins used in the measurement. The jackknife estimates are shown in green, the NaMaster Gaussian covariance based on analytical power spectra in blue, and that based on measured power spectra in black.

The two NaMaster-based estimates generally agree well with each other, indicating that the covariance is relatively insensitive to the choice of input power spectra. We note that for the highest redshift bins, the analytical-based covariance tends to have slightly smaller variances, which may reflect differences between the assumed model and the measured spectra.

Appendix C: Forecast from different covariance

We adopted a Fisher forecast for a Euclid-like survey with wider sky coverage and better photometry to estimate the 95% upper bounds of the DM parameters, as detailed in Sect. 5.3. In the forecasts, we used the analytical Gaussian covariance described in Sect. 4.2 for the Euclid-like survey. To assess the validity of this approximation, we compare the fisher forecasts for a KiDS-like survey obtained using two covariance prescriptions: the NaMaster Gaussian covariance (used as fiducial in the measurements) and the analytical covariance used in the Fisher forecasts. This allowed us to quantify the impact of the covariance modelling on the inferred sensitivities and to evaluate the robustness of the Euclid-like forecasts.

The resulting constraints on the annihilation cross-section ⟨σ_annv⟩ and decay rate Γ_dec are shown in Fig. C.1 for the HIGH clumping scenario. For the KiDS-like case, the forecasts obtained with the NaMaster covariance constructed from measured and analytical power spectra are in good agreement over most of the parameter space, showing that the analytical approximation captures the main features of the covariance relevant for forecasting. The measurements derived from the KiDS-Legacy data, shown for comparison, are broadly consistent with the forecasted sensitivity. Small deviations are observed in some regions, particularly at low m_DM, but remain within the expected level of statistical variation. At low m_DM, the $b\overline{b}$ channel lies close to its kinematic thresholds, so forecasts near the threshold might be less robust, which accounts for the small data and forecast differences observed in that regime. The Euclid-like forecasts exhibit a clear improvement across all channels, reflecting the increased sky coverage and statistical power of the survey. This improvement is consistent between annihilation and decay scenarios.

Fig. C.1. Upper: Comparison of the forecasts from Euclid-like survey with analytical covariance, KiDS-Legacy with analytical and combined covariance and the KiDS-Legacy measurement of 95% bounds on the DM annihilation cross-section ⟨σannv⟩ for three final states of HIGH clumping scenario from combined covariance. The dot-dashed lines represent the expected cross-section $⟨{\sigma }_{\mathrm{ann}}v⟩=3\times {10}^{-26}{\mathrm{cm}}^3{\mathrm{s}}^{-1}$. Lower: Forecast for three cases and KiDS-Legacy measurement of 95% upper and lower bounds on the DM decay rate Γdec for three final states. The model value of decay rate is taken as Γdec = 5 × 10−28 s−1.

Overall, the comparison indicates that the NaMaster covariance constructed from analytical power spectra provides a reliable description of the uncertainties for Fisher forecasts, and the fiducial NaMaster covariance based on measured power spectra ensures a consistent treatment of mask-induced mode coupling in the measurements.

Appendix D: Tests of the discrepancy between DES Y3 and KiDS results

We performed a series of consistency tests to investigate the origin of the discrepancy between the highly significant cross-correlation signal reported for DES Y3 and the null detection found in this work using KiDS-Legacy shear data. Below, we briefly summarise these tests and the possible causes suggested by their outcomes.

Firstly, the difference between the DES and KiDS results cannot be explained by survey area or angular scale coverage alone. DES Y1, which covers an area comparable to KiDS-Legacy (∼1500 deg²), has already reported a significant detection, indicating that the signal is not driven solely by the larger DES Y3 footprint. Moreover, the DES analysis probes a wide range of angular scales using the tangential shear statistic, with a significant signal detected from arcminute to degree scales. These scales overlap substantially with those probed in our pseudo-C_ℓ analysis, suggesting that the discrepancy is not due to DES accessing scales unavailable in KiDS.

Secondly, we tested whether differences in the construction and processing of the Fermi–LAT maps could explain the discrepancy. Using the UGRB maps and pipeline choices from Paopiamsap et al. (2024), we reproduced their galaxy–gamma-ray cross-correlation signal when correlating the Paopiamsap et al. (2024) UGRB maps with their galaxy samples (2MPZ and WISE–SuperCOSMOS). Repeating the same cross-correlation using our fiducial Fermi–LAT maps yields a consistent signal level, with small differences attributable to the different energy binning adopted in this work. We further applied the same Fermi UGRB maps to weak-lensing shear tracers. When cross-correlated with DES shear, we recover a significant detection with ∼6σ using cross correlations. The corresponding χ² values for each cross bin of this validation test are reported in Table D.1, where we use the Paopiamsap et al. (2024) UGRB maps together with the DES Y3 shear catalogue. The lowest Fermi energy bin is excluded because of the large PSF masking. In contrast, when cross-correlating with KiDS shear, we obtain a null result, and this remains true whether we use the Paopiamsap et al. (2024) UGRB maps or our fiducial Fermi–LAT maps. Taken together, these tests indicate that our implementation is able to reproduce the literature galaxy–gamma-ray measurements and the DES shear–gamma-ray detection when using the corresponding external data products. Moreover, the KiDS shear null detection is not driven by the particular choice of the Fermi–LAT map but persists across independent gamma-ray map constructions.

Finally, differences between DES and KiDS shear measurements are unlikely to be the dominant cause. The consistency between DES Y3 and KiDS-1000 cosmic shear results has been extensively tested and discussed in the literature (see, e.g. Dark Energy Survey and Kilo-Degree Survey Collaboration 2023), with no evidence of significant systematic offsets that could account for the observed discrepancy in the gamma-ray cross-correlation.

Based on the differences observed between the KiDS-N and KiDS-S footprints, our current working hypothesis is that foregrounds and large-scale anisotropies in the Fermi–LAT sky, arising from Galactic emission and spatially inhomogeneous exposure, lead to an inhomogeneous cross-correlation signal across the sky. Under this inhomogeneity interpretation, the effective signal amplitude within the KiDS footprint could significantly differ from DES, as this signal includes more than just the homogeneous signal expected from WIMP annihilation or decay. This hypothesis will be tested more robustly with upcoming wide-area weak-lensing surveys such as Euclid and LSST, which will allow nearly full-sky cross-correlation analyses and a more direct assessment of spatial variations in the signal.

All Tables

All Figures

Fig. 1. Top: Redshift distributions of the KiDS-Legacy gold-selected sample for tomographic bins, where the photometric redshift bin with z in [0.1, 2.0] is the sum weighted by the effective number densities listed in Table 1. The shaded vertical bands indicate the boundaries of the tomographic photometric redshift bins. Bottom: Weak lensing window functions for six tomographic bins in KiDS-Legacy.

In the text

Fig. 2. Window functions for gamma-rays produced by DM annihilation, decay, and astrophysical sources across an energy range of 0.5 − 1000 GeV. The annihilation process assumes ⟨σannv⟩ = 3 × 10−26 cm3 s−1 (Steigman et al. 2012), where the three scenarios are considered: HIGH (purple), MID (orange), LOW (green), representing different amounts of DM concentration for subhaloes. The DM decay (black) assumes a particle decay rate of Γdec = 5 × 10−28 s−1, taken as a representative benchmark consistent with current limits. Both processes are modelled with the assumption of the final state in $b\overline{b}$ pairs and mDM = 100 GeV for this figure. The contribution from unresolved blazars, which were taken as the only astrophysical source of the UGRB in this analysis, is shown in pink.

In the text

Fig. 3. Model of $C_{\ell }^{\mathrm{g}\kappa }$ for three clumping cases of DM annihilation, decaying DM, and astrophysical background. The models and formatting are the same as in Fig. 2. The models assumed the redshift n(z) for z in the range of [0.1, 2.0] and energy range from 0.5 to 1000 GeV. The effect of the Fermi-LAT PSF was not considered. The blazar component shown in pink represents the unresolved blazar contribution to the astrophysical background adopted in this analysis.

In the text

	Fig. 4. Fermi PSF in harmonic space bℓ in ten different energy bins. It shows an inverse relationship between the energy and the suppression of the PSF effect. The dashed black line represents ℓ = 1500, referring to the upper limit of ℓ in the cross-power spectrum measurement.
In the text

Fig. 5. Cross angular power spectrum Cℓgκ between the Fermi-LAT gamma-ray intensity map and KiDS-Legacy weak lensing data when using a NaMaster covariance. E- and B-mode cross-spectra are shown as black and red points, respectively. Panels are labelled by the tomographic redshift bin zi and gamma-ray energy bin Ej. The reduced χν2 (ν = 5) with respect to null signal for the E-mode in each bin is also provided in each subplot.

In the text

Fig. 6. Upper bounds at 95% confidence on the DM annihilation cross-section ⟨σannv⟩ for $b\overline{b},{\mu }^{-}{\mu }^{+}$, and τ−τ+ states making use of six redshift bins and ten energy bins. The unresolved blazar contribution to the UGRB is included as the astrophysical background component in the analysis. For reference, dashed lines indicate the limits obtained without subtracting the blazar contribution. The thermal relic cross-section for WIMPs (Steigman et al. 2012) is shown as dashed-dot lines. The upper bounds are presented for three DM clumping scenarios: LOW (green), MID (orange) and HIGH (purple).

In the text

	Fig. 7. Upper bounds at 95% confidence on the DM decay rate Γdec for $b\overline{b},{\mu }^{-}{\mu }^{+}$, and τ−τ+ states include the unresolved blazar contribution to the UGRB as the astrophysical background (solid curves). The dashed curves show the limits obtained without the blazar component.
In the text

Fig. 8. Left: Upper bounds at 95% confidence on the DM thermally averaged cross-section ⟨σannv⟩ of HIGH clumping case for the τ+τ− state, compared to constraints from local probes. This includes results from neutrino detectors such as IceCube (Abbasi et al. 2023) and combined results from ANTARES and IceCube (Albert et al. 2020); the gamma-ray telescope H.E.S.S. (Abdallah et al. 2016), combined results from Fermi-LAT and MAGIC (MAGIC Collaboration 2016); and the red dotted lines represent thermal relic cross-section. Right: Upper bounds on the DM decay rate Γdec for the τ+τ− state, compared to limits from IceCube (Aartsen et al. 2018) for neutrino surveys, HAWC (Albert et al. 2018) and searches of cosmic rays and the interstellar medium with Fermi-LAT (Ackermann et al. 2012) from gamma-ray telescopes. The cosmological constraints from this work (black solid lines) and work from Paopiamsap et al. (2024) (black dashed lines) are also shown. GC: Galactic centre; dSph: dwarf spheroidal galaxies; Cosmo: cosmological probe.

In the text

Fig. 9. Upper: Forecast of 95% upper and lower bounds on the DM annihilation cross-section ⟨σannv⟩ for $b\overline{b},{\mu }^{-}{\mu }^{+}$, and τ−τ+ states from cross-correlation of ten energy bins from Fermi-LAT and ten redshift bins from Euclid-like survey. The expected value ⟨σannv⟩ = 3 × 10−26 cm3 s−1 is shown as dot-dashed lines for all cases. Lower: Forecast of 95% upper and lower bounds on the DM decay rate Γdec for $b\overline{b},{\mu }^{-}{\mu }^{+}$, and τ−τ+ states. The dot-dashed lines represent the expected value for the forecast with Γdec = 5 × 10−28 s−1.

In the text

	Fig. A.1. Correlation matrix for the cross-correlation between weak lensing data and gamma-ray intensity maps across different redshift and energy bins. The redshift bins (z1 – z6) are the six tomographic bins listed in Table 1. The energy bins (E1 – E10) span 0.5–1000 GeV from low to high energy.
In the text

Fig. B.1. Comparison of diagonal variances for the gamma–shear cross-power spectra. Each subplot corresponds to one tomographic redshift bin (columns, from z1 to z6) and one Fermi–LAT energy bin (rows, from e1 to e10). The diagonal elements Var(Cℓgκ) are shown in the five multipole bins used in the measurement. The jackknife estimates are shown in green, the NaMaster Gaussian covariance based on analytical power spectra in blue, and that based on measured power spectra in black.

In the text

Fig. C.1. Upper: Comparison of the forecasts from Euclid-like survey with analytical covariance, KiDS-Legacy with analytical and combined covariance and the KiDS-Legacy measurement of 95% bounds on the DM annihilation cross-section ⟨σannv⟩ for three final states of HIGH clumping scenario from combined covariance. The dot-dashed lines represent the expected cross-section $⟨{\sigma }_{\mathrm{ann}}v⟩=3\times {10}^{-26}{\mathrm{cm}}^3{\mathrm{s}}^{-1}$. Lower: Forecast for three cases and KiDS-Legacy measurement of 95% upper and lower bounds on the DM decay rate Γdec for three final states. The model value of decay rate is taken as Γdec = 5 × 10−28 s−1.

In the text