Constraints on dark matter and astrophysics from tomographic γ -ray cross-correlations

,


I. INTRODUCTION
The unresolved γ-ray background (UGRB) is the collective emission observed after subtracting the diffuse Galactic contribution and detected extragalactic sources.It originates from unresolved astrophysical sources, including star-forming galaxies [1,2], misaligned active galactic nuclei (AGNs) [3], blazars [4,5], millisecond pulsars [6], and the interaction of cosmic rays with the extragalactic background light [7] (see [8] for a review).More interestingly for fundamental physics, the UGRB may also be used for the indirect detection of particle dark matter.
Weakly-interacting massive particles (WIMPs) are one of the most promising candidates to make up a substantial fraction of the dark matter (DM) needed to explain cosmological observations [9][10][11][12][13].With masses in the range of m χ ≃ 1 GeV − 1TeV 1 , if produced thermally in the early Universe and then decoupled from the plasma, their relic density would be naturally of the order of the measured global dark matter abundance [14,15].In most scenarios, decay or annihilation of WIMPs into Standard Model particles leads to the emissions of γ-ray photons.For this reason, data from γ-ray observatories, such as the Fermi-LAT telescope, has been widely used to place constraints on WIMPs as a dark matter candidate [16,17].
Both the extragalactic astrophysical sources of the UGRB listed above, as well as any potential contribution from DM processes, would trace the same large-scale structures, and hence cause anisotropies in this background.These anisotropies have been detected and studied via measurements of the UGRB power spectrum [18][19][20][21], revealing important information about its energy spectrum and potential composition.However, significant additional information can be obtained from the cross-correlations of the UGRB with other tracers of the same large-scale structure.Firstly, crosscorrelations with high signal-to-noise tracers of structure are able to tease out low-significance signatures in the data [e.g.22,23].Secondly, the cross-correlation with tracers at different redshifts allows us to reconstruct the joint dependence of the UGRB on redshift, energy, and physical scales, which is vital to separate the astrophysical sources of the signal from its potential DM contributions.
In this sense, various probes of structure offer different advantages.The spatial distribution of galaxies is generally the highest signal-to-noise tracer of the matter fluctuations and, if redshift information is available, allows for an accurate tomographic reconstruction of the signal [24][25][26][27][28][29][30].The main difficulty, however, is ensuring an accurate modelling of the relationship between galaxy and matter overdensities, particularly on small scales [31][32][33].Cosmic shear, caused by the weak gravitational lensing of background galaxies, is a direct tracer of the matter fluctuations, and is therefore immune to this challenge.However, shape measurement noise significantly reduces the sensitivity of this tracer, even for the densest galaxy samples and, as a cumulative effect along the line of sight, separating the γ-ray signal into its contributions at different redshifts is less straightforward [34][35][36][37][38][39].Another promising cross-correlation tracer is the positions of galaxy clusters.As the most massive objects in the Universe, clusters of galaxies contain significant amounts of dark matter, and hence are promising environments to detect the associated γray signal [40][41][42][43].Cross-correlation studies with other less standard tracers have been carried out in the literature, including Cosmic Microwave Background anisotropies (CMB) [44], lensing of the CMB [45], the Cosmic Infrared Background [46], the late-time 21cm signal [47], the thermal Sunyaev-Zel'dovich effect [48], and high-energy neutrino events [49].Finally, some of the tightest constraints on WIMPs have been obtained from the study of the UGRB around local structures [50][51][52][53].
In this paper, we study the cross-correlation of maps of the UGRB from 12 years of Fermi-LAT data with tomographic maps of the galaxy overdensity constructed from the 2MASS Photometric Redshift survey (2MPZ) and the WISE-SuperCOSMOS photometric survey (WI-SC), covering the redshift range z ≲ 0.4.We interpret the associated signal in terms of both astrophysical γ-ray sources and DM processes.The main improvements from our analysis with respect to previous similar works (e.g.[25,28]) is the use of newer, more sensitive Fermi-LAT data (12-year dataset as opposed to the 8-year release used in previous works, together with its updated point-source mask), with correspondingly enhanced constraints on DM decay and annihilation, and the adoption of an agnostic modelling approach.This allows us to compress our data into the measurement of a few model-independent quantities 2 .As we show, this has key advantages, enabling an easy interpretation of our measurements in the context of arbitrary DM particle interactions, as well as providing robust methods to identify the presence of astrophysical contributions in the data.
This paper is structured as follows.Section II presents the theoretical background used to model the UGRB and its crosscorrelation with galaxies.The datasets used in our analysis are described in Section III.The data analysis methodology we employ is presented in Section IV.Section V then presents the results of our analysis, and the corresponding constraints on DM interactions and γ-ray astrophysics.We summarise these and conclude in Section VI.

II. THEORETICAL MODEL A. Power spectra and the halo model
Our two observables, the angular galaxy overdensity δ g (see Section II B), and the γ-ray intensity I γ (see Section II C), are sky projections of 3-dimensional fields, and can be written in general as where U is the 3D counterpart of the projected field u, χ is the radial comoving distance, and q u (χ) is the radial kernel defining the projection.The angular power spectrum of two such fields, u and v, can be calculated as where P UV (k, z) is the power spectrum of the associated 3D fields, and k ℓ ≡ (ℓ + 1/2)/χ.The equation above assumes the applicability of Limber's approximation [54,55], which holds when the radial kernels of the fields involved are significantly broader than their typical correlation length (this is the case for the observables studied here).The 3D power spectra can be described using the halo model [56][57][58].In this formalism, the power spectrum receives two contributions: the "two-halo" term, corresponding to the contribution from pairs of matter elements in different haloes, and the "one-halo" term, corresponding to matter elements belonging to the same halo: where we have omitted the redshift dependence for brevity.These two contributions are given by where n(M) is the halo mass function, b h (M) is the halo bias, U(k|M) is the Fourier transform of the halo profile for field U around a halo of mass M, P L (k) is the linear matter power spectrum, and ⟨• • • ⟩ denotes averaging over all possible realisations of haloes of the same mass.The bias-weighted observable ⟨bU⟩ will be of particular interest in Section II E.

B. Galaxy Clustering
We cross-correlate the γ-ray intensity maps with projected maps of the galaxy distribution in a set of tomographic redshift bins.Concretely, we produce maps of the projected galaxy overdensity δ g ( n) ≡ n g ( n)/ ng − 1, where n g is the angular number density of galaxies, and ng is its sky average.This quantity is related to the underlying 3D galaxy overdensity ∆ g via where p g (z) is the sample's redshift distribution, and dz/dχ = H(z) is the expansion rate at redshift z.
For simplicity, we assume that galaxies are a biased tracer of the underlying dark-matter overdensity where b g is the galaxy bias.We also assume that dark matter is distributed in haloes according to a truncated Navarro-Frenk-White (NFW) profile [59] (see Appendix A 1 for details).Thus, within this simple model, the radial kernel and halo profile associated with galaxy clustering are: where ρM,0 is the mean matter density today.To test that our final results are not sensitive to the details of the galaxy-halo relation, we repeat parts of our analysis for a more sophisticated model, using the so-called halo occupation distribution approach (HOD).The model is described in more detail in Appendix A 2.

C. Gamma ray intensity
γ-ray observations can be used to construct maps of the UGRB intensity, defined as the number N o of observed photons detected per unit time t o , detector area A o , photon energy ε o , and solid angle dΩ o along direction n: where the subscript o denotes quantities measured in the observer's frame.The contribution to the total observed intensity from sources in a distance interval dl e along the line of sight is given by where dχ is an interval of comoving distance.The ratio of observed to emitted photons can for now be written as where the optical depth τ accounts for the absorption of γ-ray photons by the extragalactic background light [60].
For simplicity, in what follows, we will label energies measured in the observer's frame with an uppercase E, and restframe energies simply with a Greek lowercase ε.
The total intensity can then be calculated by integrating along the line of sight, obtaining: where ṅγ ε is the emissivity: the physical The form of ṅγ ε (e.g. its dependence on cosmological and astrophysical parameters) depends on the process giving rise to γ-ray emission.In the next two sections we will present models to describe γ-ray emission from DM decay and annihilation, and from regular astrophysical processes. 3As opposed to comoving.

Dark matter annihilation
The probability per unit time for one DM particle to annihilate against another is where n DM is the number density of DM particles, and ⟨σ v⟩ is the annihilation cross section times the relative velocity between particles averaged over the velocity distribution.The total number of annihilations taking place per unit volume and time is therefore where ρ DM is the DM density, and m DM is the DM particle mass.The factor of 1/2 must be included to avoid doublecounting unique pairs of annihilating particles.
To obtain the emissivity due to annihilation we simply multiply the number density rate of annihilations by the spectrum dN ann γ /dε (i.e. the number of photons emitted per unit energy) for each annihilation:

Dark matter decay
Calculating the emissivity for DM decay is even simpler.The probability per unit time for a DM particle to decay d p dec /dt is simply given by the decay rate Γ.Thus, the number of decays per unit volume and time is To obtain the emissivity we simply multiply this by the decay spectrum dN dec γ /dε:

γ-ray intensity maps from DM decay and annihilation
Combining the results for decay and annihilation, the γ-ray intensity can be written schematically as where ε ≡ E(1 + z) is the rest-frame energy of the emitted photons.Above, C(z) is a radial kernel involving only cosmological quantities, F(ε) is a function depending only on particle physics properties of DM, and ∆ γ (x) is a threedimensional field tracing the large-scale structure.The specific form of these quantities depends on the emission process (decay or annihilation): • Annihilation: • Decay: Here, δ (x) ≡ ρ DM (x)/ ρDM − 1 is the DM overdensity, ρ c,0 ≡ 3H 2 0 /8πG is the critical density today, with H 0 the Hubble constant, Ω DM ≡ ρDM,0 /ρ c,0 is the fractional dark matter density, and we have used the fact that the DM density evolves as ρDM (z) = ρDM,0 (1 + z) 3 .
Finally, we will make use of intensity maps integrated over a finite energy bin These are then related to the quantities defined above via where With this, the radial kernel associated with the γ-ray intensity maps is: At the redshifts and energies under study, the optical depth τ can be safely ignored [60], and we will do so in what follows.
The halo profile associated with DM decay is simply the NFW profile, described in Appendix A 1: The case of annihilation is more complicated.Since annihilation is proportional to the mean of the squared DM density, the signal is highly sensitive to the amount of substructure, i.e. fluctuations around the mean NFW profile for haloes of a given mass, normally manifested in the form of subhaloes.We will present constraints from annihilation assuming three different models to describe the boost factor to the annihilation halo profile associated with substructures.Ordered by the amplitude of the associated boost factor, these are the substructure models of [61], [62], and [63] (labelled SC14, M16, and G12 hereon, respectively).Details of these four models are described in Appendix A 3, and follow the description in [64].Unless otherwise stated, our fiducial constraints on annihilation will assume the M16 substructure model.
Another important aspect of the model for annihilation is the minimum mass over which haloes contribute effectively to the signal [64] (see Eqs. 4 and 6).As in previous works (e.g.[25]), we will set this to be M min = 10 −6 M ⊙ , corresponding to a typical WIMP free-streaming mass.We explore the uncertainty associated with this choice in Appendix B.
4. Modeling F F(ε) depends solely on the fundamental properties of the DM particles: mass, cross section/decay rate, and the associated spectra.The specific form of F thus depends on the particle physics model used to describe DM and its interactions.Here we will take an agnostic approach and instead attempt to reconstruct the form of F directly from the data.To do so, we will model F(ε) simply as a step-wise function: where ), and the amplitudes F n are free parameters of the model.Inserting this into Eq.29, the model for the intensity map in the i-th bin takes a particularly simple form: where and Above, the two last Heavyside functions ensure that the integral is zero when there is no overlap between both energy intervals.
For simplicity, and since we will only explore relatively low redshifts (z ≲ 0.4), and hence the effects of redshifting from the source to the observer frame are mild relative to the energy bin widths, we will use the same bin boundaries used to construct the intensity maps, {E i }, as the edges {ε n } used to model F in Eq. 33.
Finally, note that the step-wise parametrisation of F(ε) (Eq.33) allows us to write the cross-correlation between the g-th galaxy sample and the i-th γ-ray intensity map as where The template power spectra C g,i,n ℓ thus depend solely on cosmological quantities (distances, redshifts, and power spectra of various powers of the matter overdensity), while all the particle-physics properties are compressed into the linear amplitudes F n .Our aim will therefore be to reconstruct F n from the measured cross-correlations.

E. Astrophysical γ-ray emission
Since astrophysical γ-ray sources trace the same large-scale structure as the dark matter structures that could contribute to the UGRB through decay or annihilation, their contribution may dominate the cross-correlation with the galaxy overdensity, and is the main contaminant for this type of study.Here we develop a simple and generic scheme, based on the halo model, to interpret our measurements purely in terms of an astrophysical signal.Our treatment follows closely that used by [65] to constrain the star formation rate density from crosscorrelations of the Cosmic Infrared Background.
Consider Eq. 5 in the large-scale limit (i.e. on scales larger than the typical size of a halo).In this regime, we can write the 2-halo contribution to the galaxy-γ-ray cross power spectrum as where b g is the galaxy bias, and we have defined the biasweighted mean γ-ray emissivity with d Ṅ/dε the specific luminosity (total number of photons per unit time and energy interval emitted by the halo).Assuming we know the galaxy bias b g , and the matter power spectrum P mm (k), the amplitude of the cross-correlation is therefore sensitive to the mean γ-ray emissivity (weighted by halo bias).The 1-halo term is more difficult to interpret since it depends on the covariance between the galaxy density and the γ-ray emissivity within haloes.However, following [28,65], we can assume the scale dependence of this term to be effectively flat in harmonic space.This is an acceptable approximation since the Fermi point-spread function (PSF) is too broad to resolve dark-matter haloes in detail.
Using this approximation, and assuming the galaxy samples under study have relatively narrow redshift bins (compared to the typical variation of ⟨b ṅγ ε ⟩), we can write a remarkably simple model for the cross-power spectrum between the g-th sample of galaxies and the i-th γ-ray intensity map, at energy E i : where z g is the mean redshift of the sample, C g,i 1h is the 1-halo contribution to this power spectrum, and the template Q g ℓ is (see Eq. 14) Note that ⟨b ṅγ ε ⟩ in Eq. 41 is evaluated at the rest-frame energy E(1 + z g ), and that itself is an intrinsic function of redshift (which we do not state explicitly above for brevity).
We can thus reconstruct the energy and redshift dependence of the γ-ray emissivity by fitting the model in Eq. 41 to the measured galaxy-γ-ray correlations using ⟨b ṅγ ε ⟩ and C g,E 1h as free parameters.

III. DATA A. The Fermi-LAT 12-year data
In this work we use 12 years of γ-ray observations from the Fermi Large Area Telescope (Fermi-LAT), which we process using the Fermi Tools and FermiPY [66]4 , following a similar procedure to Ackermann et al. [21].We reject the lowest quartile of photons according to their PSF (PSF0) and divide the remaining photons into 100 logarithmically spaced energy bins in the range [0.5248, 1000] GeV.To project and bin these observations spatially onto a sky map, we use the healpy package [67] and the HEALPIX pixelation scheme.More information on HEALPIX can be found in [68].The angular resolution of the data (parametrised by the HEALPIX resolution parameter nside), is nside = 1024 which corresponds to a pixel size of ∼ 3.4 arcminutes.
We impose masks that remove all bright γ-ray sources detected by Fermi-LAT and recorded in the LAT 12-year Source Catalog (4FGL-DR3) 5 , comprising of more than 6000 sources.The masks are different in each energy bin due to the energy-dependence of the Fermi point-source response function, which is summarised as follows.The reconstruction of the direction of photons detected by Fermi is imperfect, and the uncertainty depends strongly on energy.This effectively leads to a smearing of the γ-ray maps, which can be quantified by the PSF.In real space, the PSF corresponds to the probability of measuring an incoming photon direction that differs from its true direction by an angle θ .This leads to a suppression of power in the γ-ray maps on scales smaller than the typical extent of the PSF.In harmonic space, this is quantified by the harmonic transform of the real-space PSF, which we plot in Fig. 1.We see that b ℓ is consistently 1 for a larger range of ℓ values at higher energy bins than lower bins.
In particular, there is significant suppression on medium-large ℓ scales, with relatively low energy bins (energies 0.52 GeV to 15 GeV) suffering the most from this effect.We include the PSF in our fiducial analysis by multiplying all theoretical power spectra templates by the corresponding PSF function in harmonic space.
As well as masking all point sources, we mask the region of the sky for which the fiducial galactic diffuse emission template (gll iem v07) exceeds three times the isotropic template.Repeating our analysis using a factor four instead, we verified that the results presented here are unaffected by this choice.The isotropic template is spatially constant, with a spectrum given by the Fermi Isotropic Spectral Template 6 .The galactic template is calibrated using a model of inverse Compton emission and spectral line surveys of HI and CO and infrared tracers of dust column density and is described in more detail in Acero et al. [69].
After applying these masks, we refit the amplitudes of the isotropic and galactic templates in each of the 100 energy bins, maximising a Poisson likelihood.The residuals to this fit are corrected for the Fermi exposure and then summed into 12 logarithmically spaced energy bins.The edges of these bins are listed in the second column of Table II.These summed residual maps are the maps used in the remainder of the analysis. 6https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html
When analysing the 2MPZ and WI-SC data, we follow closely the treatment of [77,78].We divide the sample into 6 redshift bins, with 2MPZ comprising the lowest bin, and WI-SC being split into 5 bins with equal photometric redshift width (∆z ph = 0.05) in the range 0.1 ≤ z ph ≤ 0.35.These are the same redshift bins used in [77,78], and further details about the resulting galaxy samples are provided in those references.Unlike in [77,78], we characterised the redshift distribution of each of the six redshift bins employed using the direct calibration (DIR) method of [79], using a cross-matched spectroscopic sample including sources from SDSS DR14 [80], 2dFGRS [81], WiggleZ [82], GAMA [83], SHELS [84], VIPERS [85], and AGES [86].Details of the method can be found in [87].Galaxy weights were found through a search of the 20 nearest neighbors to each spectroscopic source in the multi-dimensional space of observed magnitudes (8 dimensions for 2MPZ, 4 dimensions for WI-SC).The resulting redshift distributions are shown in Fig. 2.
Galaxy overdensity maps were constructed from the number counts of galaxies in pixels, and corrected for contamination from extinction and stars as described in [78].The residual contamination is limited to large scales (ℓ ≲ 10), which we remove from our analysis (see Section IV A).The associated sky mask was constructed as described in [78].

A. Power spectrum estimation
To compute all power spectra we use the MASTER algorithm [88] implemented in Namaster7 [89].The estimator can be summarised as follows (refer to [89] for further details).Observations of a given field are usually on an incomplete sky.We can relate the true map of the field u( n), to the observed map ũ( n) via ũ( n) = g u ( n)u( n), where, in the simplest case, g u ( n) is a binary mask such that g u ( n) = 1 if a pixel at n has been observed and g u ( n) = 0 otherwise.In general, g u may be understood as a local weight, which can be tuned to optimise the precision of the estimated power spectra.To compute the C ℓ in terms of the observed maps, we use the "pseudo-C ℓ " estimator: where ũℓm and ṽ * ℓm are the harmonic coefficients of the masked maps ũ( n) and ṽ( n) respectively.Multiplying a true map by a fixed sky mask results in statistical coupling between the different harmonic coefficients, which leads to coupling between different power spectrum multipoles, thus making the pseudo-C ℓ estimator biased with respect to the true C ℓ : where is the "modecoupling matrix" (MCM), which can be computed entirely in terms of the pseudo-C ℓ of the masks.In principle, we would invert the MCM to obtain an unbiased estimator of the true C ℓ .In general, the MCM is not guaranteed to be invertible but can be approximately inverted by binning the pseudo-C ℓ into bandpowers.The procedure is summarised as follows: 1. We bin the pseudo-C ℓ to obtain the mode-coupled bandpowers Cuv q = ∑ ℓ∈q B ℓ q Cuv ℓ and the binned MCM , where B q is a binning operator for the q-th bandpower.
2. We invert the binned MCM to obtain the decoupled, bandpowers Ĉuv In principle, in order to relate our estimated power spectra to their theoretical prediction, we would need to apply the same binning operation to the latter.This can be achieved by convolving it with the so-called "bandpower window functions", which encode the effects of mode-coupling, binning, and MCM-inversion: In practice, since the power spectra under study are noisy and rather flat, we can simply approximate where ℓ q is the central multipole in band q.We use a mixed binning scheme, consisting of linear bin widths of ∆ℓ = 30 between 0 ≤ ℓ ≤ 240 and then logarithmic bins until ℓ max = 3071 with ∆log 10 (ℓ) = 0.055.
To avoid any potential systematics in the galaxy overdensity maps on large angular scales (e.g.extinction or star contamination), we remove the first bandpower (ℓ < 30).To avoid numerical inaccuracies in the spherical harmonic transforms (see Appendix A in [90]) we limit the smallest scale to ℓ = 2N side = 2048.This leaves a total of 24 bandpowers for each cross-correlation.The total data vector, containing 6 × 12 cross-correlations, thus has 1728 elements.
To estimate the statistical uncertainties of these power spectra, we made use of the analytical approach outlined in [91].The method assumes that all fields involved are Gaussiandistributed, and accurately accounts for the impact of modecoupling caused by the presence of sky masks.The method requires an estimate of the power spectra of all fields involved.For this, we use the pseudo-C ℓ estimate corrected by the sky fraction of the masks involved: Here C ab,Cov ℓ is the power spectrum between fields a and b used to estimate the covariance matrix, Cab ℓ is the corresponding pseudo-C ℓ , w a is the mask of field a, and ⟨• • • ⟩ denotes averaging over all pixels in the sky.[91] found that using this recipe was able to accurately account for the impact of mode-coupling in the covariance matrix.To test the validity of the analytical covariance, we compared it with the statistical uncertainties found via jackknife resampling in some of the cross-correlations, obtaining a reasonable agreement between both approaches.

B. Galaxy bias modelling
Our fiducial analysis assumes a linear galaxy bias relation.We fix the linear bias parameter of each redshift bin to the values measured by [77,78], corrected for the fiducial cosmology used here.Given the large uncertainties of the crosscorrelations analysed in this work, propagating the small statistical error in the linear bias measurement, obtained from the measurements of the galaxy auto-correlation, has a negligible effect on our final constraints.
For the same reasons, we expect the simplifying assumption of a linear bias relation to be sufficiently accurate for the analysis presented here.To test this, we repeat our analysis assuming a more sophisticated model of the galaxy-halo connection, in the form of a Halo Occupation Distribution (HOD) parametrisation.In particular, we use the HOD model of [92], as implemented in [78], for the analysis of the same galaxy samples used here.In this model, the number of galaxies in haloes of a given mass is determined in terms of two free parameters: • M min : the halo mass for which the mean number of central galaxies is 0.5.
• M 1 : the typical mass of haloes hosting one satellite galaxy.
All other parameters of the model (described in detail in Appendix A 2) were fixed to the values used in [78].As we will show, the simpler linear bias parametrisation recovers results that are compatible with those found using this HOD model given the measurement uncertainties.The values of the linear galaxy bias, as well as the best-fit values of (M min , M 1 ), found in [78] for each redshift bin, are shown in Table I.

C. Likelihood analysis
The models described in Section II depend on a set of primary parameters ⃗ θ .These are either the dark matter parameters F n , describing the function F(ε) (see Sections II D 3 and II D 4), or the astrophysical parameters {⟨b ṅγ ε ⟩,C g,E 1h } (see Section II E).We constrain these parameters from our data vector d, consisting of all 6 × 12 galaxy-γ-ray crosscorrelations.To do so, we assume that d follows a Gaussian likelihood where t is the theory prediction for our data vector, C is the covariance matrix, and K is a normalisation constant.
When the theory prediction is a linear function of the model parameters, which is the case for both F n and {⟨b ṅγ ε ⟩,C g,i 1h }, and assuming wide, uniform priors (as we will do here), the posterior distribution is Gaussian in those parameters by construction, and their mean and covariance can be calculated analytically.Specifically, consider the case: where T is an model-independent N d × N p matrix (with N d and N p equal to the number of data points and free parameters, respectively).The posterior mean of ⃗ θ (which coincides with its maximum a-posteriori -MAP), and its covariance are simply given by: (50) When constraining F(ε), ⃗ θ ≡ {F n }, and the elements of T are (see Eqs. 37 and 38) When constraining the astrophysical parameters, we can write the parameter vector as θ α,g,i , with α ∈ {1, 2}, and The elements of T, in turn, are (see Eqs. 41 and 42): After having measured these linear model parameters from a given set of power spectra, we make use of these measurements to constrain other secondary parameters (e.g.DM annihilation cross section or decay rate, or the functional dependence of the γ-ray emissivity on energy and redshift).In those cases, we still use a Gaussian likelihood, as in Eq. 48, using the measured linear parameters as data (i.e. with d and C given by ⃗ θ MAP and Cov θ in Eq.50 for the linear parameters).This is mathematically consistent, since the linear dependence of the theory on these primary parameters ensures that they follow a Gaussian distribution.Since these secondary parameters are in general not linearly related to the primary ones, in this case we make use of Markov-Chain Monte-Carlo (MCMC) techniques in order to obtain parameter constraints.For this, we use the emcee package8 [93].
All theoretical predictions were computed with the Core Cosmology Library (CCL [94]), making use of the CAMB Boltzmann solver [95] to compute the linear matter power spectrum.The non-linear power spectrum, where necessary, was computed using the HALOFit implementation of [96].We fixed all cosmological parameters to the best-fit values found by Planck [11].In halo model calculations, we made use of the halo mass function parametrisation of [97], the halo bias model of [98], and the concentration-mass relation of [99].We use a spherical overdensity halo mass definition with overdensity parameter ∆ = 200 with respect to the critical density (see Eq. A2).We verified that, changing the mass function and concentration parametrisations (to that of [100] and [101], respectively) leads to only small changes in our results, at the level of 10-20%, in line with the usual accuracy of most halo model ingredients [102].V. RESULTS

A. Power spectrum measurements
We compute the angular power spectrum for the galaxyclustering maps in 6 redshift bins with the γ-ray maps in 12 energy bins, resulting in 72 cross-correlations in total.We shall refer to an individual galaxy-γ-ray cross-correlation with the notation: redshift bin × energy bin.
As a preliminary exploration of the data vector, we compute a rough estimate of the signal-to-noise ratio (SNR) in each power spectrum as where N d = 24 is the number of data points in any given power spectrum, and χ 2 0 is the χ 2 statistic (see Eq. 48) for a null model (t = 0).Since we are not fitting the data to a particular model at this stage (we do not consider any noise models), we shall take Eq.54 as a crude estimate of the detection significance for our cross-correlations.The rationale behind Eq. 54 is that the expectation value of the χ 2 for purely noise-dominated Gaussian noise is N d , and thus the equation is a measure of the departure with respect to this expectation.Note that we will use a more principled definition for SNR when adopting a given model in the rest of the paper.
In Fig. 3, we present a visualisation of SNR for the 6 × 12 power spectra analysed in this work.From this qualitative estimate, we observe that the signal measurement is dominated by the intermediate energy bins (3 GeV −120 GeV), and that the 2-3 highest-and lowest-energy bins are largely noise-dominated.This is as expected, as Fermi-LAT is most sensitive at intermediate energies [103].The highest SNR is achieved in the 4 × 4 power spectrum, with a tentative ∼ 5σ detection of the signal.
Since each of the 72 power spectra explored carries only a relatively small part to the total SNR of the data, determining whether a signal is being consistently detected is not straightforward by analysing any individual spectrum.As we will see in Section V C, the data supports a model in which the γ-ray emissivity scales with rest-frame energy approximately as ε α , with α ∼ −2.3.We can use this to produce measurements of the cross-correlation between galaxy redshift bin g coadded over all 12 energy bins as follows: Note that we only do this here for visualisation purposes.All steps of our analysis in the next sections are based on the raw set of 6 × 12 power spectra.Fig. 4 shows the six coadded cross-correlations, as well as the approximate SNR of each bin.We see that we obtain consistently positive cross-correlations, detected between 2.7σ and 6.7σ .The figure also shows the best-fit theoretical prediction for the astrophysical model described in Section II E (solid blue line), and for the DM decay and annihilation models of Section II D (dotted red and dashed orange lines).These results are presented in more detail in Sections V B and V C.Although all models are able to provide a reasonably good fit to the data, these preliminary results already show that the DM models are not flexible enough to fully reproduce the redshift evolution of the signal.We will quantify this more accurately in Sections V B 1 and V C.
Using the best-fit theoretical predictions for the three models explored in the next sections, we can produce more accurate estimates of the detection significance of the signal, given by where, as before, χ 2 0 is the χ 2 with respect to a null model, χ 2

M
is the χ 2 with respect to the best-fit prediction within model M, and N θ is the number of free parameters of the model.Using this definition, together with the best-fit models presented in the next sections (DM decay and annihilation from the model-independent reconstruction of F(ε), and astrophysical sources), we obtain: show the best-fit predictions for the astrophysical model presented in Section V C, while the dashed orange, and dotted red lines show the best predictions for DM annihilation and decay, respectively, described in Section V B 1 (obtained from the model-independent reconstruction of F(ε), and hence independent of the specific decay/annihilation channel and WIMP mass).
We thus find that the cross-correlation between the diffuse γray emission as measured by Fermi, and the 2MPZ + WI-SC galaxy samples is detected at the ∼ 8-10σ level.

Measuring F
We produce measurements of the step-wise amplitudes F n characterising the energy dependence of the modelindependent quantity F(ε) (see Section II D 4) from our power spectrum measurements, using the analytical solution for linear parameters described in Section IV C. We do so for each redshift bin individually (combining all cross-correlations with the 12 energy bins), and for the full data vector containing all 6 × 12 cross-correlations.
The results are shown in Fig. 5 for annihilation (top panel) and decay (bottom panel).We can see that, within each redshift bin, as well as for the coadded measurements of F n , the function F(ε) seems to be well described by a power-law behaviour.Focusing on the energy bins with highest detection significance (e.g.ε ≃ 10 GeV), we further observe a coherent evolution of the signal with redshift, with the amplitude of F growing with z.Since, by construction, F(ε) depends only on particle-physics quantities, and should therefore not be redshift-dependent, this already provides qualitative evidence that the observed signal is not consistent with a purely DM origin.
The coadded measurements of F n are listed in Table II for DM decay and annihilation, including all 4 models of substructure explored.Since the model used to parametrise the cross-correlations is fully described by the F n (having fixed the cosmological and galaxy bias parameters), these measurements compress all the particle physics information contained in our 1728-element data vector into only 12 numbers.Furthermore, since these measurements have been obtained without assuming any specific decay/annihilation channels, they can be used directly to place constraints on arbitrary DM particle physics models.It is worth noting that the uncertainties on the measured F n in different energy bins are not entirely uncorrelated.The full covariance matrix is made publicly available together with the measurements of F.
We find that our measurements of F are robust to the choice of a simple linear bias parametrisation used here, with results changing by less than 10% and 40% of the statistical uncertainties for annihilation and decay, respectively, when using the HOD model described in Section IV B.
As an example, we fit the measured F n in each redshift bin to a simple power-law model of the form We choose a pivot frequency ε 0 = 20 GeV, and fit for the amplitude F 0 and spectral index α as free parameters.We imposed flat priors on both parameters: α ∈ [−3, 0] for both decay and annihilation, and F dec 0 ∈ [0.001, 3] × 10 −29 GeV −2 s −1 and F ann 0 ∈ [0.01, 3] × 10 −29 cm 3 GeV −3 s −1 .We sample the resulting posterior distribution using emcee.
We tabulate the constraints on F 0 and α for each redshift bin, and for the coadded measurements of F n in Table III.For the coadded measurements, we obtain the following spectral indices: α ann = −2.31± 0.08 and α dec = −2.19± 0.08 for annihilation and decay respectively.[21] modelled the un-  TABLE II.Coadded measurements of F n for annihilation and decay, for each of the 12 γ-ray energy bins weighted by the energy spectrum with an index of α = −2.3,together with their 68% uncertainties.The annihilation constraints are presented for the 4 different models of substructure described at the end of Section II D 3 and in Appendix A 3.  As noted in [21], these spectral indices are compatible with blazar-like sources and mis-aligned active galactic nuclei at energies above and below a few GeV, respectively.Our measurement for the annihilation and decay spectral indices lie in between both of these measurements, with the annihila-tion spectral index being compatible with the first measurement.Our measurements lie somewhat in between both values.[30] obtained a spectral index of −2.75 +0.71 −0.46 , in agreement with our results.[104] obtained a slightly steeper index of −2.41 ± 0.05 for a single power-law by modelling detected Fermi-LAT sources and the diffuse Galactic γ-ray emission for energy ranges of 0.2 to 100 GeV (probing most of the energy range considered in this work).Our annihilation spectral index is compatible with [104].Our results are in excellent agreement with those of [28] (see e.g.their Table 5), TABLE III.Constraints on the free parameters of the power-law model for F(ε) (see Eq. 57) for DM annihilation and decay, in each of the 6 redshift bins (first 6 rows), and for the coadded measurements (last row).

⟨z⟩
F ann 0 × 10 30 0.06 0.149 ± 0.047 0.137 ± 0.038 2.32 ± 0.19 2.24 ± 0.16 0.13 0.116 ± 0.031 0.138 ± 0.036 2.37 ± 0.16 2.28 ± 0.15 0.19 0.163 ± 0.032 0.237 ± 0.046 2.36 ± 0.13 2.22 ± 0.12 0.24 0.144 ± 0.032 0.245 ± 0.055 2.37 ± 0.13 2.27 ± 0.12 0.29 0.126 ± 0.033 0.302 ± 0.070 2.20 ± 0.16 2.06 ± 0.12 0.34 0.153 ± 0.039 0.430 ± 0. who determined the spectral index of diffuse γ-ray emission through cross-correlations of earlier Fermi data with a large set of galaxy samples, including 2MPZ and WI-SC.Our measurements of F 0 for decay and annihilation are shown in Fig. 6 as a function of redshift.As already anticipated in Figs. 4 and 5, we observe a trend for the amplitude of F(ε) to grow with redshift, especially for DM decay, which should not be the case if the detected signal were caused by DM processes.Quantifying the evidence for this trend is not entirely straightforward a priori, given the non-negligible overlap between the six redshift bins explored here.We will do so in Section V C.This trend indicates that at least a sizeable fraction of the measured signal is likely caused by astrophysical, baryonic sources, rather than DM processes.Thus, without a reliable prediction for what this fraction is, in what follows we will treat our measurements as contributing to the upper bound on the contribution to the diffuse γ-ray background from DM decay and annihilation.Specifically, when quoting an upper bound on a given DM property µ (e.g.⟨σ v⟩ of Γ), we will quote μ + 2σ µ , where μ and σ µ are the mean and standard deviation of the inferred quantity.Note that, since we will fix the WIMP mass when constraining Γ and ⟨σ v⟩, these parameters are still linearly related to the F n measurements, and to our original data vector of cross-correlations, and thus their statistical uncertainties are Gaussianly distributed.

Constraints on DM parameters
To obtain the constraints on DM annihilation and decay properties, we follow the prescription outlined in Eq.IV C and employ Eq. 50.The data vector d is now our measurements F(ε) and we set ⃗ θ MAP to be the DM quantities we wish to constrain: the velocity-averaged annihilation rate ⟨σ v⟩ and the decay rate Γ for DM annihilation and decay respectively, as a function of the WIMP mass m DM .We set the T matrix of the linear regression (Eq.50) to be (recall Eq. 23 and Eq.26): To determine T, we use the Fermi Tools and FermiPy to obtain the theoretical spectrum dN/dε for both annihilation and decay.Specifically, we use DMFitFunction9 in FermiPy to calculate the spectrum dN/dε as a function of DM mass, and decay channel.Since the signal follows roughly a power law, with the spectral indices quoted in the previous section, instead of convolving the spectrum with the energy bandpass of each channel, we simply evaluated the function at effective energies weighted by this power-law spectrum: using the values of α for decay and annihilation found above.We compute ⟨σ v⟩ and Γ at fixed DM mass, considering logarithmically spaced values within the range [2, 400] GeV.In Fig. 7, we present the 95% confidence limit constraints on ⟨σ v⟩ and Γ for different Standard Model (SM) particleanti-particle channels.The constraints derived here are contingent on the following implicit assumptions/conditions: 1.The only progenitor of γ-rays are DM particles, i.e we have not taken into account any astrophysical models that could potentially source the γ-rays .
2. Annihilation occurs non-relativistically with orbital angular momentum L = 0 (so-called s-wave annihilation).This allows us to treat the decay spectrum to be equivalent to the annihilation spectrum but with half of the DM mass.
3. The DM particle decays/annihilates into single SM particle-anti-particle channels with an branching ratio of 1, i.e. we omit the possibility of decaying/annihilating into multiple channels per event.
4. We impose kinematic constraints, such that a decay can only occur if m DM ≥ 2m SM and similarly m DM ≥ m SM for annihilation, where the subscript SM denotes the corresponding SM particle.Hence some particle channels (e.g t t) do not probe the full range of DM masses.
With these assumptions in mind, we find that the constraints derived for both Γ and ⟨σ v⟩ are competitive with previous studies that have considered cross-correlations between γ-rays and other cosmological tracers: weak-lensing [35][36][37], galaxy surveys [25,26], as well as constraints derived from local structures [50-52, 105, 106].In the following section, we shall only consider the b b, µ + µ − and τ + τ − channels for comparative purposes and for sake of clarity.For Γ, we find that across all three channels, we obtain constraints that are of the same order of magnitude at m DM ∼ 100 GeV when compared to constraints found by [26].For the b b channel and at a DM mass of ∼ 10 GeV, we obtain constraints that are less stringent than those of [37] for their conservative (DM + Astro) and optimisitc (only DM) models.But, as we go towards higher masses, at around ∼ 100 GeV, we obtain constraints that are slightly tighter than both these models.For both the µ + µ − and τ + τ − channels, our constraints are nearly an order of magnitude tighter for all masses greater than ∼ 10 GeV than both DM and DM + Astro models.At ∼ 10 GeV, all three channels are roughly at the same order of magnitude as [30] with the b b constraint derived in this work being weaker by a factor of 6.The constraints on the b b channel are up to two orders of magnitude weaker than those derived from other cross-correlations with galaxy surveys [25,26], with µ + µ − and τ + τ − channels being rouhgly the same order of magnitude.
Likewise, we find that the constraints on ⟨σ v⟩ derived in this work are competitive with those found in the literature.In Fig. 8, we present an exclusion plot for b b, µ + µ − and τ + τ − annihilation channels where the shaded areas represent regions that can be ruled out from our 95 % limit constraints in three substructure regimes: G12, M16, and SC14.The constraints on ⟨σ v⟩ across all substructure models have been computed by extrapolating the halo mass down to 10 −6 M ⊙ for the sake of concurrence with constraints presented in current literature.We explore the effect of changing this lower bound of the mass in Appendix B. When comparing equivalent substructure models, we find that the constraints derived in this work are up to roughly two orders of magnitude more stringent than [37], obtained from cross-correlation with weak lensing, across all three channels for all three substructure cases.In particular, the G12 constraints obtained are only marginally stronger than both the conservative (DM + Astro) and optimistic (DM only) models for the b b channel.In comparison to constraints derived via cross-correlations with galaxy catalogues for equivalent substructure regimes,  [30] from a field-level analysis of low-redshift structures.The dashed red and blue lines show the constraints obtained from dwarf spheroidals (dSph), as computed by [52] and [51], respectively.The solid black line shows the thermal relic limit [107].
our constraints are slightly less stringent than [25] by an order of magnitude at a mass of ∼ 10 GeV.This may be because our measurements are limited by the amplitude of the detected astrophysical signal, and hence our upper bound is not determined solely by the statistical error achieved.Although weaker for the b b channel, we find that the constraints derived in this work for both the τ + τ − and µ + µ − channels for equivalent substructure regimes are up to an order of magnitude more stringent than those found in [25].More in detail, our constraints for the τ + τ − channel achieves roughly the same constraining power with only slightly tighter constraints in the mass range 10 GeV − 100 GeV for the SC14 substructure regime.For constraints derived through DM-only models [26], we find that we can provide roughly the same constraining power for the b b channel.We obtain roughly an order of magnitude tighter constraints for µ + µ − and τ + τ − .
It is worth noting, however, that the detailed modelling of the impact of substructure on γ-ray emission from DM annihilation is a large source of theoretical uncertainty in this analysis.In [26] the difference in the most optimistic case (referred to as LOW in their study, and corresponding to our M16 model) and the most conservative case (HIGH, corresponding to G12) is roughly an order of magnitude.By comparison, we find the difference between these two models to be more modest.Note, however, that, our M16 model [62] is a refinement of the M16 model used in their work [61].The difference between these two cases is a more detailed modelling of halo concentrations.This comparison highlights the sensitivity of the constraints to the substructure modeling parameters.The constraints found assuming the potentially optimistic model of G12 for the subhalo boost factor gives rise to a factor of ∼ 3 difference in the constraints from the SC14 substructure case.The difference in both extreme cases is modest in comparison to [26], who report a difference of an order of magnitude.
In the case of b b, we are able to exclude the thermal relic cross section around m DM ∼ 10 GeV, in a range of WIMP masses that depend heavily on the substructure model assumed.In the most optimistic case (G12), for both the b b and τ + τ − channels, the constraints obtained are comparable to those found through the study of dwarf spheroidal galaxies (dSph) [51,52].Finally, our constraints are more stringent than those found by [30] from a field-level analysis of lowredshift structures, by a factor that ranges between ∼ 2 and ∼ 10 depending on the substructure model used.

C. Constraints on the diffuse astrophysical γ-ray background
In Section V B 1 we presented evidence that the DM kernel F(ε) evolves with redshift, in a manner that would be incompatible with DM decay or annihilation as the sole origin of the cross-correlation between the diffuse γ-ray background and the positions of galaxies in 2MPZ and WI-SC.This prevents our interpretation of the measured signal in terms of fundamental DM properties, and limits our ability to obtain a tighter upper limit on them.
In this section, we instead take a more agnostic approach, and use our measurements to place constraints on the energy and redshift dependence of the γ-ray emissivity from our data, regardless of its physical origin.To do this, we follow the methodology outlined in Section II E, which allows us to determine the bias-weighted γ-ray emissivity ⟨b ṅγ ε ⟩(z) from the 10 0 10 1 10 2 10 3 [GeV] two-halo regime of the galaxy-γ-ray cross-correlation.Using the linear reconstruction method of Section IV C, we obtain the measurements shown in Fig. 9.It is worth emphasizing that we obtain a single measurement of ⟨b ṅγ ε ⟩ for each individual cross-correlation: the cross-correlation between galaxies at mean redshift z g and the γ-ray map with mean observed energy E i , provides a measurement of ⟨b ṅγ ε ⟩ at redshift z g and at the rest-frame energy ε gi ≡ E i (1 + z g ).The position of each measurement along the x axis of the figure corresponds to the associated rest-frame energy (hence the displacement between points at different redshifts).
As in the case of F(ε), we find that the γ-ray emissivity follows a power-law-like dependence on energy.To quantify this further, we fit the measured values of ⟨b ṅγ ε ⟩ in each redshift bin g to a power-law of the form: with ε 0 = 20 GeV as before.As in Section V B 1, since the measured values of ⟨b ṅγ E(1+z g ) ⟩ are linearly related to the cross-power spectra, they also follow a Gaussian likelihood.We assume flat priors for both ⟨b ṅγ ε 0 ⟩ g and α g , and sample the posterior distribution using emcee.The results are shown in Table IV.In all cases, the power-law model is a good fit to the data, with reasonable χ 2 values and associated probabilities.These results are shown as black points with error bars in Fig. 10.We see that all redshift bins recover compatible values for the spectral index which, as expected, is also in agreement with the spectral index measurement found for F(ε) in Section V B 1. Furthermore, there is a clear evolution of the γ-ray emissivity with redshift.
To quantify this redshift evolution, and compare it with the evolution expected for DM decay and annihilation, we fit a global model of the form to the values of ⟨b ṅγ ε ⟩ measured from all pairs of redshift and energy bins.The amplitude and spectral index parameters (⟨b ṅγ ε 0 ⟩ 0 , α) have the same interpretation as before, with the power-law index β parametrising the redshift evolution.As before, we use a Gaussian likelihood and assume flat priors for all parameters, obtaining the following constraints: The model has a best-fit χ 2 of 89.5 for 69 degrees of freedom and thus provides an adequate fit to the data.The shaded bands in Fig. 10 show our 1σ constraints on ⟨b ṅγ 20 GeV ⟩(z) for this model, obtained from the MCMC chains.The posterior distribution on the redshift evolution parameter β is shown in Fig. 11.
We can use our constraints on β to quantify the evidence against a purely DM-related origin for the signal we measured based on its redshift dependence.This is straightforward for DM decay since, in this case, the emissivity is directly proportional to the dark matter overdensity.The value of β for DM decay can then be found by simply comparing the radial kernels of both models.The astrophysical kernel is proportional to (1 + z) −3 (see Eq. 14), whereas the DM kernel is constant (see Eqs. 21 and 25).Hence, for decay β decay = 3.This value is shown as a vertical red dashed line in Fig. 11.
For annihilation, the comparison is less straightforward.The annihilation radial kernel is proportional to (1 + z) 3 (see Eqs. 21 and 22).However, annihilation is proportional to (1 + δ ) 2 instead of (1 + δ ), and the power spectrum of both quantities evolves differently with redshift.The effective radial kernel for annihilation is therefore proportional to Here P δ ,δ (k, z) is the power spectrum of the matter overdensity, whereas P δ ,δ 2 (k, z) is the cross-spectrum between 1 + δ and (1 + δ ) 2 .The value of ∆β ann depends on scale, as well as on the model used to describe substructure.Over scales relevant for this analysis (0.01 Mpc −1 < k < 1 Mpc), and for the various models of substructure explored here, ∆β ann varies in the range [−2.8, −1.7].The corresponding allowed range of β ann is shown as a vertical blue.band in Fig. 11.
Considering the highest value of β compatible with DM annihilation, we thus find that our measurements are incompatible with a purely DM-related origin to the diffuse γ-ray background at the 2.8σ level (considering the most extreme annihilation case).This justifies our interpretation of the detected signal as an upper bound of the emission from DM processes, rather than an indirect detection of DM decay or annihilation.

VI. CONCLUSIONS
Observations of the UGRB present a unique probe to improve our understanding of high-energy astrophysics, as well as a window to constrain, or potentially detect, WIMP dark matter.In this work, we have analysed cross-correlations between γ-ray intensity maps, covering the energy range E ∈ [0.5 GeV, 1 TeV], with the galaxy overdensity at 6 different redshift ranges covered by the 2MPZ and WI-SC surveys.
We detect a positive cross-correlation between both datasets at the level of 8 − 10σ , confirming the extragalactic nature of the UGRB.The sensitivity of this measurement, and the availability of redshift data, allows us to reconstruct the dependence of the signal on rest-frame energy and redshift, enabling us to interpret it in the context of both dark matter and astrophysics.
In the context of WIMP searches, we make use of linear regression methods to reconstruct, from our cross-correlation measurements, a function F(ε), defined in Eqs.23 and 26, that depends solely on the particle physics parameters governing DM decay and annihilation.We do so in a modelindependent way that compresses all the information in our cross-correlations into a set of 12 numbers, {F n }, characterising the dependence of this function of rest-frame energy.The reconstructed function has a power-law energy dependence (F(ε) ∝ ε α ) with a spectral index α ≃ −2.3.Repeating this for each redshift bin independently, we find evidence that the amplitude of F(ε) evolves monotonically with redshift.This would not be possible for a signal sourced only by dark matter processes, implying the presence of astrophysical contamination in our measurements.Nevertheless, we may use the measured value of F(ε) to place an upper-bound constraint on the WIMP decay rate and the annihilation crosssection.Although, due to this contamination, the constraints are systematics-limited, we find bounds that are competitive with those obtained by other groups targeting similar largescale structure cross-correlations, as well as constraints from local structures.As in all other γ-ray WIMP searches, our annihilation constraints are hampered by the theoretical uncertainty on the impact of substructure.Constraints may vary by up to two orders of magnitude between a model with no substructure, and the most extreme substructure model studied here [61].In the most optimistic case, we can rule out, at the 95% confidence level, the thermal relic bound for WIMP masses of a few tens of GeVs assuming complete decay into b b quarks or τ + τ − leptons.Better control over the impact of substructures on the expected annihilation signal must be achieved before such a claim can be undisputed.
Arguably the main source of uncertainty in our ability to constrain DM physics, is our inability to characterise and clean the contamination from astrophysical sources in the signal, which, as the analysis of Section V C shows, must be present.This forces us to interpret the cross-correlation signal, detected at the ∼ 8 − 10σ level, as contributing to the upper bound on the DM contribution.Together with the uncertainty in the theoretical model from substructures, this constitutes arguably the main impediment in obtaining reliable constraints on DM annihilation.This could be improved in the future by incorporating (and marginalising over) a physicsbased model of all potential astrophysical sources.Removing these contaminants at the map level (e.g. by detecting and masking out further extragalactic point sources), would likely lead to larger improvements.Finally, other sectors of the data, including higher-order statistics (e.g. through a fieldlevel analysis of all tracers over the same volume [30]), as well as including γ-ray auto-correlations and measurements of the isotropic signal, could significantly enhance the con-straints presented here.Future Fermi-LAT data (e.g.masking with the new 14-year data source catalog 4FGL-DR4), and its combination with denser and more reliable galaxy catalogs from next-generation surveys, such as DESI [108], the Rubin Observatory LSST [109], Euclid [110], or the Roman Space Telescope [111], have the potential to improve the constraints presented here using similar analysis methods, assuming that the modelling challenges outlined above are tackled.

FIG. 1 .
FIG. 1.The point-spread function (PSF) in harmonic space b ℓ for the 12 Fermi-LAT energy bins considered in this work.The vertical markers indicate the mean of the multipole bandpowers used in this work.

34 FIG. 2 .
FIG.2.The redshift distributions of the six redshift bins considered in this work.The mean redshift of each respective bin is tabulated with the plot.

FIG. 3 .
FIG. 3. The signal-to-noise ratio ( SNR) estimates given a null model of each respective 6 × 12 galaxy-clustering and γ-ray crosscorrelations.The negative values arise from the sign function in Eq. 54

SNR decay = 8 . 2 ,FIG. 4 .
FIG.4.Cross-correlations between Fermi-LAT and the 6 2MPZ and WI-SC galaxy samples.Each panel shows the power spectrum coadded over γ-ray energies assuming a power-law spectrum with spectral index α = −2.3 and inverse-variance weighting.The measured power spectra are shown as black circles with error bars.Empty circles show the absolute value of a given measurement when negative.The solid blue lines show the best-fit predictions for the astrophysical model presented in Section V C, while the dashed orange, and dotted red lines show the best predictions for DM annihilation and decay, respectively, described in Section V B 1 (obtained from the model-independent reconstruction of F(ε), and hence independent of the specific decay/annihilation channel and WIMP mass).

FIG. 5 .
FIG.5.The tomographic and coadded measurements of F n for annihilation (top panel) and decay (bottom panel).The horizontal shift in data points is for visualisation purposes.The points marked by a cross are the negative values of F n , and the points marked by a dot are the positive values.

FIG. 6 .
FIG.6.Constraints on the amplitude of the F(ε) function, in the power-law model of Eq. 57, for DM annihilation (top panel) and decay (bottom panel).The points with error bars show constraints in each redshift bin, with the corresponding mean redshift shown in the x axis.The shaded horizontal bands show the constraint obtained from the coadded measurement of F(ε).

FIG. 7 .
FIG. 7. 95% upper bounds on the DM annihilation cross section (left) and decay rate (right) as a function of WIMP mass.The different lines show constraints assuming a single decay/annihilation channel into Standard Model particle-antiparticle pairs (see legend).
FIG.8.95% upper bounds on the DM annihilation cross-section for the b b, µ + µ − , and τ + τ − channels.The constraints found from our data are shown as shaded band.Constraints found for the four different substructure models described in Appendix A 3 are shown as purple bands with different shadings.The yellow dashed line shows the constraints found by[30] from a field-level analysis of low-redshift structures.The dashed red and blue lines show the constraints obtained from dwarf spheroidals (dSph), as computed by[52] and[51], respectively.The solid black line shows the thermal relic limit[107].

34 FIG. 9 . 2 αFIG. 10 .
FIG.9.Constraints on the bias-weighted mean γ-ray emissivity as a function of rest-frame energy ε and redshift (see legend).Obtained from the analysis of the galaxy-Fermi-LAT cross-correlations using the model described in Section II E.

FIG. 11 .
FIG.11.Constraints on the power-law index characterising the redshift evolution of the γ-ray emissivity using the global model in Eq. 61.The black solid line shows the posterior distribution obtained from the set of cross-correlations studied here.The vertical red line and shaded bands show the values of β expected for DM decay and annihilation.The measured signal is incompatible with a DM origin at the ∼ 3σ level.