UvA-DARE (Digital Academic Repository) Detection of Cross-Correlation between Gravitational Lensing and γ Rays

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

In recent years, many γ-ray sources have been identified, yet the unresolved component hosts valuable information on the faintest emission. In order to extract it, a cross-correlation with gravitational tracers of matter in the Universe has been shown to be a promising tool. We report here the first identification of a cross-correlation signal between γ rays and the distribution of mass in the Universe probed by weak gravitational lensing. We use data from the Dark Energy Survey Y1 weak lensing data and the Fermi Large Area Telescope 9-yr γ-ray data, obtaining a signal-to-noise ratio of 5.3. The signal is mostly localized at small angular scales and high γ-ray energies, with a hint of correlation at extended separation. Blazar emission is likely the origin of the small-scale effect. We investigate implications of the large-scale component in terms of astrophysical sources and particle dark matter emission. Introduction.-Astronomy at γ-ray frequencies represents a promising avenue for both astrophysics and particle physics. On one hand, the most violent phenomena in the Universe produce high-energy photons that travel all the way to Earth. Thus, they bring us information about the physics of rare events such as supernovae and the behavior of matter under extreme conditions, as in pulsars and active galactic nuclei (AGNs). On the other hand, the most elusive form of matter in the cosmos-dark matter, which represents about 25% of all the Universe's energy content-is believed to consist of an exotic fundamental particle, which may annihilate or decay into Standard Model particles and thus produce cosmic messengers including γ-ray photons. In the weakly interacting massive particle (WIMP) scenario, or for any hypothetical dark matter particle with mass in the GeV range or higher, dark matter particle annihilation or decay almost necessarily results in photons at γ-ray energies. Therefore, γ-ray astronomy represents a promising means to investigate the fundamental nature of dark matter. However, the faintness of the expected emission makes it very difficult to identify such a signal.
Since 2008, the large area telescope (LAT) mounted on the Fermi satellite has been performing the most detailed observations of the extragalactic γ-ray sky and resolved 5065 γ-ray sources in the energy range of 50 MeV-1 TeV [1]. Once the point sources and Galactic emission are removed, the remaining γ-ray photons form the so-called unresolved γ-ray background (UGRB).
A method to discriminate between nonthermal γ-ray emission due to astrophysical sources and possible dark matter annihilation or decay in the UGRB has been proposed in Ref. [2]. This method relies on cross-correlations of UGRB maps with maps of other tracers of the underlying structure on cosmological scales, such as the weak gravitational lensing effect, or the clustering of galaxies and galaxy clusters (see also Refs. [3,4]) and cosmic microwave background (CMB) lensing [5,6]. These are direct gravitational probes of matter, most of which are thought to be dark matter. The energy, redshift, and scale dependence of the aforementioned cross-correlations have the potential to disentangle signatures due to astrophysics from dark matter (see also Ref. [7]). More generally, the method can provide valuable information on the redshift distribution and on the clustering properties of the unresolved γ-ray source populations, including blazars, AGNs, and star-forming galaxies.
Since cross-correlations of the UGRB with gravitational lensing have been proposed as a probe, several observational attempts have followed [8][9][10][11], but none so far have detected the signal. Here, we report the first detection of such a cross-correlation. We used 108-month γ-ray data from Fermi-LAT and first year (Y1) shear measurements from the Dark Energy Survey (DES). In the following, we describe details of the analysis and discuss the results.
Analysis and results.-The observable we probe is the cross-correlation between the unresolved component of the γ-ray emission and gravitational shear. To this aim, the Fermi-LAT data have been preprocessed to produce the relevant energy-dependent response functions of the detector and full-sky maps of photon intensities in several energy bins. Resolved γ-ray sources and the bright Galactic plane emission have been masked with energy-and flux-dependent masks, in order to minimize the sky fraction removal. Furthermore, we have subtracted a model of the Galactic plane emission. Galactic foreground emission does not lead to false detection of a cross-correlation, since it does not correlate with the large-scale structure measured by gravitational shear, but it increases the variance of the measurements (see Supplemental Material [12] and, e.g., Refs. [8,9,11,50,51]). The weak lensing information is extracted by measuring the mean tangential ellipticity of source galaxies in the DES footprint around pixels weighted by their UGRB flux. The shear catalog is divided in redshift bins in order to perform a tomographic analysis. As an illustration of the overlapping area between DES and Fermi-LAT, Fig. 1 shows the DES footprint and the Fermi-LAT map for photon energies in the 1-10 GeV interval.
We measure the cross-correlation between the UGRB and gravitational shear through its two-point angular correlation function. Specifically, we adopt the following estimator (see also Ref. [52]): where Ξ signal Δθ h ;ΔE a ;Δz r is the correlation function in the configuration space of the two observables measured in different angular (Δθ h ), γ-ray energy (ΔE a ), and lensing sourcegalaxy redshift (Δz r ) bins. The correlation is obtained by summing the products of tangential ellipticity of source galaxies i relative to a pixel j, e r ij;t , multiplied by the Fermi-LAT photon intensity flux in the ath energy bin and in pixel j, I a j . The sum runs over all unmasked pixels j and all sources i in the redshift bin of the shear catalog, and it is performed in each of the different photon energy bins (labeled by a) and source-galaxy redshift bins (labeled by r). Lastly, R is the mean response of ellipticity to shear for sources in the redshift bin, determined by the METACALIBRATION algorithm [53,54] to be between 0.54 and 0.73 for the source-galaxy redshift bins used here.
From the correlation function, we remove Ξ random Δθ h ;ΔE a ;Δz r , the measurement of tangential shear around random lines of sight. This is done by setting I a j;random ¼ 1 anywhere within the sky region used for γ-ray measurements in that energy bin and zero elsewhere. This reduces additive shear systematic effects, random very-large-scale structures, or chance shear alignments relative to the mask. The random subtraction, while not affecting the expected signal, lowers the variance at large angular separations (see also Refs. [52,55]). We analyze the data in 12 logarithmically spaced angular bins with radii between 5 and 600 arc min, 9 photon energy bins between 0.631 and 10 3 GeV, and 4 redshift bins defined by 0.20 < hzi < 0.43, 0.43 < hzi < 0.63, 0.63 < hzi < 0.90, and 0.90 < hzi < 1.30, where hzi is the estimated expectation value of galaxy redshift from DES. The energy bins used in the analysis and the corresponding 68% and 95% containment angles of the Fermi-LAT pointspread function (PSF) are shown in Table I. These sum up to a total of 432 bins for the cross-correlation measurement. The analysis is performed blindly, i.e., on multiple variants of the measurements including artificial versions, in order to avoid experimental bias in measurement and interpretation of the signal. See the Supplemental Material for details [12].
The result of the measured cross-correlations, averaged over all energy and redshift bins, is shown in Fig. 2 in terms of the estimator ΞðθÞ defined in Eq. (1). Note that the data points reported on both panels are the same, although confronted with different models. A clear positive cross-correlation is observed, especially at small angular separations.
In order to determine the statistical significance of the signal, we test the deviation of the measurement from a null signal (null hypothesis of pure noise) by means of a phenomenological model, which aims at capturing the general expected features of the cross-correlation signal without resorting to any specific, detailed modeling of its physical origin (in the next section, we will instead adopt a physical model to provide insights on the origin of the cross-correlation). In the halo-model approach, all mass in the large-scale structure of the Universe is associated with virialized dark matter halos, and the correlation function can thus be decomposed into the so-called one-and twohalo terms (1h and 2h in formulas hereafter). The former refers to the correlation between two points in the same physical halo, the latter to the case in which the two points belong to two different halos. Pointlike sources contribute at small angular scales with a one-halo term, while at large scales they produce a two-halo term resembling the largescale structure matter distribution. In our case, we use the fact that the spatial extent of the one-halo term is smaller than the beam window function of the Fermi-LAT. Then a phenomenological model can be constructed as where E a and z r are the central values of the energy (measured in GeV) and redshift bins, and hI a i is the measured photon flux.Ξ a PSF−like ðθÞ is the Legendre transform of the beam window function (or PSF) integrated in the ath energy bin (in arbitrary units, being merely a template for the one-halo term due to pointlike γ-ray sources) and Ξ ar 2h-like ðθÞ is the Legendre transform of a generic two-halo (i.e., large-scale) contribution, also convolved with the Fermi-LAT beam window function. Correlation functions with a hat have flux units, while those without a hat are normalized to the γ-ray flux as in Eq. (1), and are therefore dimensionless. The two normalizations A 1 and A 2 , spectral indices α 1 and α 2 , and redshift evolution indices β 1 and β 2 are free parameters of the model [56]. γ-ray sources typically have energy spectra that can be well approximated by a power law, and so it is assumed in Eq. (2). For simplicity, we also assume a power-law scaling in redshift. Best fits and confidence intervals of the parameters are found in a Markov chain Monte Carlo likelihood analysis.
The first statistical method adopted to quantify the presence of a signal, and its significance, against the null hypothesis relies on the Δχ 2 test statistics, with the χ 2 defined in the usual way, i.e., where Ξ data is the data vector, Ξ th is the theoretical crosscorrelation for the models outlined above, described by the parameter set P mod , and Γ is the data covariance matrix, detailed in the Supplemental Material [12]. (All angular, energy, and redshift bin indexes have been omitted for simplicity of notation.) The Δχ 2 is defined as Þ computed at the model parameter values P ⋆ mod that best fit the data and χ 2 null referring to no signal, i.e., Ξ th ¼ 0. The second estimator is the matched filter signal-to-noise ratio (see, e.g., Ref. [57]), in analogy to Δχ 2 mod , we shall later refer to SNR mod ≡ SNRðP ⋆ mod Þ. In Table II we present the results on detection significance. The phenomenological model results for the full dataset show clear evidence for the presence of a crosscorrelation signal, at the level of SNR mod ¼ 5.3. Since the matched filter based on the phenomenological model captures the generic features of the cross-correlation signal, without committing to any specific physical description, this best assesses that indeed a cross-correlation between gravitational shear and unresolved γ-rays emission has been observed. In order to investigate the features of the signal in more detail, we repeat the tests by subdividing the dataset according to redshift, energy, and angular separation. Specifically, low (high) z refers to the first two (second two) redshift bins; low (high) E bins are defined by being below (above) 5 GeV. and small (large) θ separates angular scales below (above) 3 times the Fermi-LAT PSF.
From Table II we infer that the signal is mostly concentrated at high energies and small angles. These results point toward an interpretation in terms of pointlike sources with hard energy spectrum, broadly compatible with these sources being blazars. In fact, the best fit for the spectral index of the PSF-like one-halo component α 1 ¼ 1.81 þ0.20 −0.24 is quite hard with respect to the spectral index of the average intensity of the UGRB [58], but compatible with BL Lacertae (BL Lac) emission, which is the source population expected to be the most relevant in the range of fluxes probed by this analysis, just below the Fermi-LAT flux sensitivity threshold. Notice that this hard spectral index is in agreement with recent findings from γ-ray autocorrelation analysis [59], possibly suggesting that BL Lac objects below the threshold have slightly harder spectra than those detected individually. The energy scaling of the two-halo component is also compatible with a blazar TABLE II. Δχ 2 mod and SNR mod computed for the phenomenological and physical models, using either the full dataset or the various subsamples discussed in the text. For dark matter in the physical model, we consider the annihilation channel τ þ τ − . For the low-z case we selected the two first redshift bins, while for the high-z case the last two bins, where the bins are defined as 0.20 < hzi < 0.43, 0.43 < hzi < 0.63, 0.63 < hzi < 0.90, and 0.90 < hzi < 1.30; the low-E subset is defined for energies below 5 GeV, while the high-E is for energies above this value. Finally, the small (large) θ cases correspond to data points below (above) 3 times the Fermi-LAT PSF. origin, though this term shows lower statistical significance than the one-halo component. Concerning the redshift dependence of the signal, the statistical significance is almost equally distributed among the lower and higher redshift bins. The allowed regions for the parameters of the phenomenological model are shown in Fig. 3, while the cross-correlation function for the best fit of the phenomenological model are shown in the left panel of Fig. 2: the PSF-like term due to pointlike sources well reproduces the behavior of the measured cross-correlation up to about 1 deg scale. We note here that for the subset of high E and small θ, comprising 88 data points, we do obtain a distinctive signal without application of the matched filter. The χ 2 null ¼ 137 for these points corresponds to a p value of 0.0006, meaning that the null hypothesis is excluded at 3.5σ in this subset.
Discussion.-In the following, we attempt a physical interpretation of the signal detected in the previous section. Star-forming galaxies and misaligned AGNs are not expected to be able to produce a sufficiently hard energy spectrum, which thus points to a dominant blazar component. Particle dark matter in terms of WIMPs can also provide a hard spectrum, especially if the annihilation channel is predominantly leptonic or, in the case of a hadronic final state, if the dark matter mass is large (above a few hundred GeV).
Blazars are compact sources and, for our purposes, they can be considered as pointlike; i.e., their size is, on average, much smaller than the Fermi-LAT PSF. Also the size of the halo hosting blazars rarely exceeds the Fermi-LAT PSF.
This has a consequence that the angular correlation function for the one-halo term essentially follows from the detector PSF. Conversely, the relevant dark matter halos have a larger angular extent, and the corresponding onehalo correlation function thus drops more slowly with angular scale. On very large scales, the correlation functions of the two components have a similar angular behavior, since the two-halo power spectra differ only by the bias terms. The fact that our signal is detected with high significance only on small scales therefore points toward blazars as the dominant source. In order to investigate this interpretation, we perform the statistical tests discussed in the previous section with a physical model, based on a detailed characterization of the components expected to produce the cross-correlation signal: blazars (BLZs), misaligned active galactic nuclei (mAGN), starforming galaxies (SFGs) and possibly dark matter (DM). The physical cross-correlation function model reads The model parameters are free normalizations for the astrophysical sources, A 1h BLZ , A 2h BLZ , A mAGN , and A SFG , the mass of the dark matter particle m DM , and its velocityaveraged annihilation rate hσ ann vi, expressed in terms of the "thermal" cross section hσ ann vi th ¼ 3 × 10 −26 cm 3 s −1 through the normalization A DM ≡ hσ ann vi=hσ ann vi th . Note that, for blazars, which represent the astrophysical component expected to dominate the correlation signal at the current level of unresolved γ-ray emission, we allow the one-and two-halo terms to be separately adjusted in the fit against the data. As for the phenomenological model, all terms depend on both energy and redshift, labeled by indices a and r, respectively.
The results are shown in Table II, where the overall significance of the presence of a signal, the preference for an origin at high energies, and small angular scales are all confirmed. However, since in this case we have specific behaviors for the correlation functions as dictated by a physical model (different for each component, contrary to the average generic case of the phenomenological model), we notice that a mild hint of large-scale correlation is present, namely, in the large-θ case. We note that both the physical and phenomenological models provide a good fit to the data according to their χ 2 (see the Supplemental Material [12]).
More details of the analysis are shown in Fig. 3, where the triangular plot of the profile likelihood distributions of the model parameters is reported. The likelihood exhibits a preference for a large one-halo term of blazars with normalization A 1h BLZ ¼ 102 þ56 −57 , while the normalizations of the blazar two-halo term shows only a (weak) upper bound. The latter picture is shared also by the other astrophysical sources (SFG and mAGN), which are shown only in the Supplemental Material [12] for brevity.
The blazar-shear cross-correlation on small scales depends on the relation between the blazar γ-ray luminosity and the host-halo mass, a quantity which is rather uncertain. For our reference model this relation has been taken from [7], where it was derived by associating the γ-ray luminosity of blazars to the mass of the supermassive black hole powering the AGN and then relating the mass of the black hole to the mass of the dark matter halo. This procedure gives MðLÞ¼2×10 13 M ⊙ ½L=ð10 47 ergs −1 Þ 0.23 ð1þzÞ −0.9 , where L is the rest-frame luminosity of blazars in the energy range 0.1-100 GeV. We can therefore translate a value of A 1h BLZ different from unity to a deviation from the reference MðLÞ relation. The value we found implies that the average mass of a halo hosting an unresolved blazar is larger than the one adopted in Ref. [7] and most likely above 10 14 M ⊙ . The cross-correlation signal with weak lensing seems therefore to be dominated by blazars residing in cluster-size halos.
The right panel of Fig. 2 shows that the cross-correlation at small angular scales requires a sizeable blazar one-halo term. It also illustrates that the best-fit analysis exhibits a mild preference for some power at large scales. This can be accounted for either by the two-halo term of blazars or by a DM contribution. The interplay of the angular, energy, and redshift behaviors of the observed signal leads to a small preference for a DM component over a pure blazar contribution in the assumed model. Misaligned AGNs and star-forming galaxies are disfavored since they do not meet the requirement of a hard energy spectrum.
To have a visual impression on the physical behaviors, we plot the energy and redshift dependence of the crosscorrelation signal in Fig. 4. The average along the angular and redshift (energy) directions of each point of the energy (redshift) spectrum is performed by computing a matched filter amplitude given by a sample model that we choose to be flat in energy and redshift, while scaling as 1=θ in angle, to approximately reproduce the expected signal, and Ξ is either the measured data or the best-fit models introduced in the main text. The error on A is given by Fig. 4 one can appreciate the small but noticeable differences in the energy and redshift scalings of FIG. 4. (Left) Energy scaling of the measured signal and best-fit models, in terms of the matched filter amplitude A of the crosscorrelation between gravitational shear and γ rays (see text for its definition). The amplitude is divided by the size of the corresponding energy bin ΔE a and multiplied by the measured photon intensity hI a i in the same bin, to show the physical differential scaling in energy. (Right) Redshift scaling of the measured signal and models, again in terms of the matched filter amplitude A introduced in the text. the models of different physical components that have been just discussed.
Notice that the blazar model we are adopting, and which is outlined in the Supplemental Material [12], is based on the current understanding of blazars as derived from the Fermi-LAT resolved sources: the small preference in the fit for a contribution with features compatible with DM might be interpreted as an indication that unresolved blazars have different properties than the resolved ones.
We find that, for a dark matter particle dominantly annihilating into the leptonic channel τ þ τ − , the best fit improves by 2.8σ, as compared to a model where DM is not included. The lack of a degeneracy between dark matter and blazar amplitudes (see Fig. 3) indicates that the two components are supported by independent features of our cross-correlation data, in particular, the small-and large-scale behavior. The best fit occurs for a dark matter particle of mass m DM ¼ð65AE 27 23 Þ GeV and annihilation rate hσ ann vi ¼ ð26 AE 17 15 Þ × hσ ann vi th . In the case of a softer energy spectrum, as the one provided by abb annihilation channel, the fit improvement is slightly lower, at 2.7σ level, with best-fit mass m DM ¼ð302AE 188 120 Þ GeV and annihilation rate hσ ann vi ¼ ð78 AE 67 43 Þ × hσ ann vi th . Let us remark that the main source of uncertainty concerning the dark matter signal described in this Letter is our ignorance on the impact of substructures. Specifically, the minimal halo mass and the amount and distribution of subhalos can significantly change the size of the signal. This, however, is common to all cosmological searches for a particle dark matter signature. Comparing our nominal model [60] to a recent development [61] based on both N-body simulations and analytical modeling, we found that the constraint on the annihilation cross-section gets shifted to higher annihilation rates by about one order of magnitude. This is due to a different amplitude of the expected signal, whereas the energy, redshift, and angular dependencies are only slightly modified, and the predicted cross-correlation function therefore just needs a larger normalization, which, in turn, is directly provided by the annihilation rate.
We conclude by summarizing the major results of this analysis. We present the first detection of the crosscorrelation between the γ-ray sky and the mass distribution in the Universe observed through gravitational lensing shear, with a significance of SNR ¼ 5.3. The bulk of this signal comes from the one-halo term of pointlike sources with a hard spectrum, most likely dominated by blazars. In addition, we find a hint for a cross-correlation on large scales with spectral and redshift behaviors that might imply that the population of blazars that are currently unresolved by Fermi-LAT has different characteristics than those obtained by extrapolation from observations of resolved blazars, or that an additional contributor to the γ-ray emission is present. The analysis of the cross-correlation of Fermi-LAT data with the forthcoming year 3 and year 5 DES dataset and improvements in modeling of the blazar population will likely clarify the source of the signal detected in this Letter and characterize it more deeply.