A NOVEL APPROACH IN THE WEAKLY INTERACTING MASSIVE PARTICLE QUEST: CROSS-CORRELATION OF GAMMA-RAY ANISOTROPIES AND COSMIC SHEAR

Weak gravitational lensing refers to the small distortions of images of distant galaxies, produced by the distribution of matter located between galaxies and the observer (Kaiser 1992; Bartelmann & Schneider 2001; Munshi et al. 2008; Bartelmann 2010). A distorted image can be described by the so-called distortion matrix, normally parameterized in terms of the convergence κ (controlling modifications in the size of the image) and the shear γ (accounting for shape distortions). While the former is a direct estimator of matter density fluctuations along the line of sight, the latter is easier to measure, through correlations in the observed source ellipticities. In the flat-sky approximation, the two generate identical angular power spectra and we thus focus on the shear.

The autocorrelation between the gravitational shear in two different directions can provide information on the clustering of the large-scale structures responsible for the lensing effect. The technique has already been used in Jullo et al. (2012), Tereno et al. (2011), de la Torre et al. (2010), and Schrabback et al. (2010) where data of the COSMOS galaxy survey allow for a measurement of the two-point correlation function of the shear at angular scales between 0.1 and 10 arcmin. Future surveys, like Pan-STARRS, Dark Energy Survey (DES; The Dark Energy Survey Collaboration 2005), and Euclid (Laureijs et al. 2011; Amendola et al. 2012), due to their larger coverage and improved sensitivities, will be able to reconstruct two-dimensional shear maps, from which one can extract the autocorrelation angular power spectrum (PS).

The same dark matter (DM) structures that act as lenses can themselves emit light at various wavelengths, including the γ-ray range. While γ-rays can be produced by astrophysical sources hosted by DM halos (i.e., star-forming galaxies (SFGs) or active galactic nuclei (AGNs)), DM itself may be a source of γ-rays, through its self-annihilation or decay, depending on the properties of the DM particle. Those γ-rays emitted by DM should therefore have the potential to exhibit strong correlation with the gravitational lensing signal.

In this Letter, we propose to study the cross-correlation of gravitational shear with the extragalactic gamma-ray background (EGB), i.e., the residual radiation contributed by the cumulative emission of unresolved γ-ray sources, as a novel and potentially relevant channel of DM investigation.

The most recent measurement of the EGB was performed by the Fermi Large Area Telescope (LAT) in Abdo et al. (2010b), covering a range between 200 MeV and 100 GeV: the emission is obtained by subtracting the contribution of resolved sources (both point-like and extended) and the Galactic foreground (due to cosmic-ray interaction with the interstellar medium) from the whole Fermi-LAT data. Unresolved astrophysical sources like blazars (Abazajian et al. 2011; Stecker & Venters 2011; Singal et al. 2012), SFGs (Ackermann et al. 2012b; Lacki et al. 2012), or radio galaxies (Inoue 2011; Massaro & Ajello 2011) contribute to the EGB but the exact amount of their contribution is still unknown. The γ-rays produced by DM annihilation or decay can also contribute to EGB (Ullio et al. 2002; Zavala et al. 2010, 2011; Cirelli et al. 2011; Fornasa et al. 2013). However, the fact that the EGB energy spectrum is compatible with a power law, without any evident spectral feature, suggests that DM cannot play a leading role in the whole energy range (Abdo et al. 2010a; Calore et al. 2012). In the angular anisotropies of the EGB emission, the DM also plays a subdominant role: indeed, a detection of a significant autocorrelation angular PS has been recently reported (Ackermann et al. 2012a; for multipoles ℓ > 100, which is the range of interest for our analysis, since there the contamination of the Galactic foreground can be neglected), but the features of such a signal (in particular its independence on multipole and energy) seem to indicate an interpretation in terms of blazars (Cuoco et al. 2012; Harding & Abazajian 2012).

The contribution of unresolved astrophysical sources to the EGB can also be analyzed by cross-correlating the gamma-ray emission with available catalogs of resolved galaxies (Xia et al. 2011). Both cosmic shear and γ-ray emission depend on the large-scale structure of the universe: because this is what generates the lensing effect and because those same structures can produce γ-rays, either from the hosted astrophysical sources or directly from DM annihilation/decay. A certain level of cross-correlation between cosmic shear and γ-ray emission is therefore expected. The key point of our analysis is to understand whether the shear/γ-rays cross-correlation is within reach of future galaxy surveys, as DES and Euclid, and under which circumstances such signal can be proficiently used to disentangle a true DM signal from the other astrophysical γ-rays sources. It is the first time, to our knowledge, that such analysis is performed. If successful, this approach could provide direct evidence that what is measured by weak-lensing surveys is indeed due to DM and is not, e.g., a manifestation of alternative theories of gravity.

The source intensity along a given direction n can be written as

$$\begin{equation} I_g (\boldsymbol n) = \int d\chi \, g(\chi,\boldsymbol n)\,\tilde{W}(\chi), \end{equation} \tag{ 1 }$$

where χ(z) is the radial comoving distance, g is the density field of the source, and $\tilde{W}$ is the window function (which does not depend on n). We then define a normalized version $W=\langle g \rangle \,\tilde{W}$, so 〈I_g〉 = ∫dχ W(χ). Expanding the intensity fluctuations of two source populations i and j in spherical harmonics, one can compute the cross-correlation angular PS (here in the dimensionless form):

$$\begin{equation} C_\ell ^{(ij)} = \frac{1}{\langle I_i \rangle \langle I_j \rangle } \int \frac{d\chi }{\chi ^2} W_i(\chi)\, W_j(\chi) P_{ij}(k=\ell /\chi,\chi). \end{equation} \tag{ 2 }$$

The definition of the three-dimensional (3D) PS P_ij is $\langle \hat{f}_{g_i} (\chi,\boldsymbol k) \hat{f}_{g_j}^\ast (\chi ^\prime,\boldsymbol k^\prime) \rangle = (2\pi)^3 \delta ^3 (\boldsymbol k - \boldsymbol k^\prime) P_{ij}(k,\chi,\chi ^\prime)$, where $f_g \equiv [g(\boldsymbol x|m,z)/\bar{g}(z)-1]$ ($\hat{f}_g$ is its Fourier transform) and the Limber approximation (k = ℓ/χ) is assumed to hold. We consider the sources to be characterized by a parameter m (typically the mass), and g(x|m) is the density field of an object associated with m, while $\bar{g}(z)=\langle g(\boldsymbol n,z)\rangle$. P_ij can be computed following the so-called halo-model approach. The two-point correlation is given by the sum of two components, the one-halo and two-halo terms, i.e., $P_{ij}=P_{ij}^{1h}+P_{ij}^{2h}$ (Scherrer & Bertschinger 1991; Ando & Komatsu 2006; Ando et al. 2007a):

$$\begin{eqnarray} P_{ij}^{1h}(k) &=& \int dm\ \frac{dn}{dm} \hat{f}_i^\ast (k|m)\,\hat{f}_j(k|m) \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} P_{ij}^{2h}(k) &=& \left[\int dm_1\,\frac{dn}{dm_1} b_i(m_1) \hat{f}_i^\ast (k|m_1) \right] \nonumber \\ &&\times \left[\int dm_2\,\frac{dn}{dm_2} b_j(m_2) \hat{f}_j(k|m_2) \right]\,P^{\rm lin}(k), \quad \end{eqnarray} \tag{ 4 }$$

where dn/dm is the number density distribution of sources, P^lin is the linear matter PS, and b_i(m) is the linear bias between the object i and matter. Note that the average of 〈g〉 is given by

$$\begin{equation} \bar{g}(z)=\langle g(\boldsymbol n,z) \rangle =\int \ dm\,\frac{dn}{dm}\int \ d^3\boldsymbol x\,g(\boldsymbol x|m,z), \end{equation} \tag{ 5 }$$

which implies that at small k (where $\hat{f}\sim \int d^3\boldsymbol x\,g(\boldsymbol x|m)/\bar{g}$) the terms in the square brackets in Equation (4) are ∼1 (except in the case of a significant bias). The two-halo term is thus normalized to the standard linear matter PS at small k, which motivates the normalization of the window function introduced above.

We aim at cross-correlating the shear signal (source i in Equations (2)–(4)) with γ-rays emitted by DM, SFGs, and blazars (source j in Equations (2)–(4)). For what concerns weak lensing, W takes the form (see, e.g., Bartelmann 2010)

$$\begin{equation*} W^\kappa (\chi)=\frac{3}{2}{H_0}^2{\Omega _m}[1+z(\chi)]\chi \int _\chi ^\infty \!\!d\chi ^{\prime }\,\frac{\chi ^{\prime }-\chi }{\chi ^{\prime }}\frac{dN}{d\chi ^{\prime }}(\chi ^{\prime }), \nonumber \end{equation*}$$

and dN/dχ represents the redshift distribution of the sources, normalized to unity area (such that we can take 〈I_i〉 = 1 in Equation (2)). For Euclid $dN/dz=A_E\,z^2\,e^{-(z/z_0)^{1.5}}$, where z₀ = z_m/1.4 with z_m = 0.9 being the median redshift of the survey and A_E is fixed by the normalization ∫dzdN/dz = 1. For DES, dN/dz = A_D (z^a + z^ab)/(z^b + c), with a, b, and c provided in Table 1 of Fu et al. (2008), and A_D fixed by the normalization. Since gravitational lensing is sourced by the potential wells of the large-scale structure, whose Poisson equation relates to the matter distribution ρ, we have g_κ(x) = ρ(x), and f_κ(x) is given by the density contrast δ(x). For the bias in Equation (4), we use the estimates in Cooray & Sheth (2002). We adopt the halo mass function dn/dm of Sheth & Tormen (1999), the halo concentration from Muñoz-Cuartas et al. (2011), and a Navarro–Frenk–White halo density profile (Navarro et al. 1997).

For the case of γ-rays from decaying DM, we again have f_d(x) = δ(x) (we assume ρ_m ≃ ρ_DM). The window function is now given by

$$\begin{equation} W^{\gamma _d}(E_\gamma,z) = \frac{1}{4\pi } \frac{\Omega _{{\rm DM}}\rho _c}{m_\chi \tau _d} J_d(E_\gamma,z), \end{equation} \tag{ 6 }$$

where m_χ and τ_d are the mass and decay lifetime of the DM particle, $J_d=\int _{E_\gamma }^\infty dE ({dN_d(E(1+z))}/{dE}) e^{-\tau (E(1+z),z)}$ with dN_d/dE(E) being the number of γ-ray photons emitted per decay event in (E, E + dE), and τ being the optical depth for absorption (Stecker et al. 2007). Note that the factor Ω_DMρ_c comes from the normalization of W, since in this case $\langle g_d \rangle =\bar{\rho }$.

The DM annihilation signal scales with ρ², thus we have $\hat{f}_{a}\equiv \tilde{u}(k|m)$ given by the Fourier transform of $\rho ^2(\boldsymbol x|m)/\langle \rho ^2 \rangle$. In the literature, equations are often written in terms of the so-called clumping factor:

$$\begin{equation} \Delta ^2(z)=\frac{\langle \rho ^2 \rangle }{\bar{\rho }^2}=\int _{m_{\rm min}}^{m_{\rm max}} dm\,\frac{dn}{dm}\int \ d^3\boldsymbol x\,\frac{\rho ^2(\boldsymbol x|m)}{\bar{\rho }^2}, \end{equation} \tag{ 7 }$$

and the window function has the form

$$\begin{equation*} W^{\gamma _a}(E_\gamma,z) = \frac{(\Omega _{{\rm DM}}\rho _c)^2}{4\pi } \frac{(\sigma _a v)}{2 m_\chi ^2} (1+z)^3\,\Delta ^2(z)\,J_a(E_\gamma,z), \nonumber \end{equation*}$$

where (σ_av) is the velocity-averaged annihilation rate (which we assume to be the same in all halos) and $J_a=\int _{E_\gamma }^\infty dE ({dN_a}/{dE})(E(1+z)) \,e^{-\tau (E(1+z),z)}$ with dN_a/dE(E) being the number of γ-ray photons emitted per annihilation in the energy range (E, E + dE). In the annihilating DM case, the predictions for both the window function and the PS heavily depend on the (unknown) clustering at small masses (i.e., on the minimum halo mass, concentration below approximately 10⁶ M_☉, and on the amount of substructures). In our benchmark case, we consider m_min = 10⁻⁶ M_☉ (typical free-streaming mass for weakly interacting massive particles, WIMPs) and include unresolved subhalos following the scheme of Kamionkowski et al. (2010) with parameters tuned as in the HIGH scenario of Section 3.3 in Fornasa et al. (2013; within our halo model, it induces only moderate boost factor ∼2). In Figures 1 (left) and 2 (left), we estimate the impact of clustering uncertainties by considering a different subhalo scheme (from Via Lactea; Kamionkowski et al. 2010), a different extrapolation of the halo concentration (from Virgo Collaboration; Fornasa et al. 2013), and a conservative case with m_min = 10⁷ M_☉ (minimum halo mass currently inferred from dynamical measurements) and no substructures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The formalism sketched in Equations (1)–(5) can also be used for a population of astrophysical sources, by replacing the mass with the source luminosity $\mathcal {L}$ as the characterizing parameter. This leads to the replacement of dm dn/dm with $d\mathcal {L}\, \Phi$, where Φ is the γ-ray luminosity function (GLF). For the range of multipoles of interest (ℓ < 10³), both blazars and SFGs can be approximated as point sources and we have $g_S(\mathcal {L},\boldsymbol x-\boldsymbol x^{\prime })=\mathcal {L}\,\delta ^3(\boldsymbol x-\boldsymbol x^{\prime })$, which leads to

$$\begin{eqnarray} P_{\kappa \gamma _S}^{1h}(k,z) &=& \int _{\mathcal {L}_{\rm min}(z)}^{\mathcal {L}_{\rm max}(z)} d\mathcal {L}\,\Phi (\mathcal {L},z)\,\frac{\mathcal {L}}{\langle g_S \rangle } \,\tilde{v}(k|m(\mathcal {L})) \nonumber \\ P_{\kappa \gamma _S}^{2h}(k,z) &=& \left[\int _{\mathcal {L}_{\rm min}(z)}^{\mathcal {L}_{\rm max}(z)} d\mathcal {L}\,\Phi (\mathcal {L},z)\, b_S(\mathcal {L},z)\,\frac{\mathcal {L}}{\langle g_S \rangle } \right]\nonumber \\ &&\times \left[\int _{m_{\rm min}}^{m_{\rm max}} dm\,\frac{dn}{dm} \tilde{v}(k|m) \right] \,P^{\rm lin}(k,z), \qquad \end{eqnarray} \tag{ 8 }$$

with $\tilde{v}(k|m)$ being the Fourier transform of $\rho (\boldsymbol x|m)/\bar{\rho }$ and $\langle g_S \rangle = \int d\mathcal {L}\,\Phi \,\mathcal {L}$. In Equations (8) a relation between the source luminosity $\mathcal {L}$ and the host-halo mass m is required. We compute the source bias b_S through the halo bias by means of $b_S(\mathcal {L},z)=b_h(m_h(\mathcal {L}),z)$, for which we need again a relation between host-halo mass and source luminosity. On the other hand, since at low redshift and in the mass range of interest b_h ∼ 1, the two-halo term is only very mildly dependent on the description of $m_h(\mathcal {L})$. For a power-law spectrum with index α, the window function is

$$\begin{equation*} W^{\gamma _S}(E_\gamma,z) = \frac{A_S(z)\,\langle g_S(z) \rangle }{4\pi \,E_0^2}\int _{E_\gamma }^\infty dE \left(\frac{E}{E_0}\right)^{-\alpha } e^{-\tau (E,z)}, \nonumber \end{equation*}$$

where E₀ = 100 MeV and A_S is a factor that depends on which specific luminosity is chosen as the characterizing parameter (as we will describe below).

The GLF of blazars is computed following the model described in Inoue & Totani (2009) with the AGN X-ray luminosity function from Ueda et al. (2003) and with the numerical value of parameters derived in Harding & Abazajian (2012) by fitting Fermi-LAT data on EGB diffuse emission and anisotropies. The spectrum is taken to be a power law with α = 2.2, and $\mathcal {L}$ is the γ-ray luminosity at 100 MeV (which leads to A_S = (1 + z)^−α). We assume that no blazars fainter than the luminosity cutoff $\mathcal {L}_{\rm min}=10^{42}$ erg s⁻¹ can exist at any redshift, while $\mathcal {L}_{\rm max}(z)$ is the maximum luminosity above which a blazar can be resolved (for 5 yr Fermi-LAT, it is computed taking F_max = 2 × 10⁻⁹ cm⁻² s⁻¹ for E > 100 MeV). The relation between halo-mass and blazar luminosity can be described through $m_h=10^{11.3}\,M_\odot (\mathcal {L}/10^{44.7}\,{\rm erg\,s^{-1}})^{1.7}$ following Ando et al. (2007b), where the blazar γ-ray luminosity is linked to the mass of the associated supermassive black hole, which is in turn related to the halo mass. The description of $m_h(\mathcal {L})$ suffers from sizable uncertainties which propagate to the prediction of the one-halo term. However, as can be seen from Figures 1 (middle) and 2 (middle), where we introduce an alternative model (model B) which dramatically increases $m_h(\mathcal {L})$ with respect to our benchmark case (model A), the blazar contribution remains largely subdominant.

For the GLF of SFGs, we follow results from the Fermi-LAT Collaboration (Ackermann et al. 2012b), which are based on the infrared (IR) luminosity function derived in Rodighiero et al. (2010), and the rescaling relation between γ-ray and IR luminosity obtained analyzing resolved SFGs (Ackermann et al. 2012b). The spectrum is assumed to be a power law with α = 2.7, similar to the Milky Way case, and $\mathcal {L}$ is the γ-ray luminosity between 0.1 and 100 GeV (which leads to A_S = (α − 2)/(1 + z)²). The dependence of the SFG–shear PS on the $m(\mathcal {L})$ relation is milder than for blazars. In this case, the relation could, in principle, be computed from the relation between γ-ray luminosity and star formation rate (SFR; Ackermann et al. 2012b), the Schmidt–Kennicutt law (connecting SFR and gas density), and the ratio of gas to total galactic mass. This leads to different relations for each different sub-population of SFGs (e.g., ellipticals are much brighter than spirals of the same mass); on the other hand, we do not have γ-ray data to compute the specific GLF of the sub-populations, thus we have to derive an effective averaged relation. Assuming a power-law scaling $m={\cal A}\times 10^{12}\,M_\odot (\mathcal {L}/10^{39}\,{\rm erg\, s^{-1}})^{\cal B}$ and a maximum galactic mass of m_max = 10¹⁴ M_☉, we can find $\cal A$ and $\cal B$ using, e.g., the Milky Way data (m ≃ 10¹² M_☉ and $\mathcal {L}\simeq 10^{39}\,{\rm erg\, s^{-1}}$) and requiring that the mass associated with the maximum luminosity ∼10⁴³ erg s⁻¹ (this can be computed from the maximum observed IR luminosity (Rodighiero et al. 2010) rescaled to γ-ray frequency (Ackermann et al. 2012b)) not to exceed m_max. We found ${\cal A} \simeq 1$ and ${\cal B} \simeq 0.5$. This is just a simple benchmark model, and we estimated the impact of the associated uncertainty (by varying $\cal A$ and $\cal B$ within reasonable ranges) in Figures 1 (right) and 2 (right).

For the sake of clarity, we focus on a benchmark annihilating (decaying) DM scenario, where the WIMP has a mass of 100 GeV (200 GeV), annihilation (decay) rate of (σ_av) = 8 × 10⁻²⁶ cm³ s⁻¹ (τ_d = 3 × 10²⁶ s) and dominant final state $\bar{b}b$. The characteristics of the DM particle are chosen to saturate (at least in one particular energy range) the EGB emission, without violating the experimental constraints.⁴ In particular, we note that, although we take DM to be a significant component of the EGB at E ≳ 1 GeV in Figure 3 (left), it is basically impossible to obtain an evidence for DM from the angular PS of γ-rays alone because the latter is dominated by the blazar contribution.

Zoom In Zoom Out Reset image size

In Figure 4, we show the ingredients of Equation (2) for the computation of the shear/γ-ray cross-correlation angular PS: the window function for the cosmic shear signal nicely overlaps with the DM window function, both for annihilating and decaying DM, while this happens only at intermediate redshifts for the SFG window function and only at high redshifts for the case of blazars. This suggests that a tomographic approach could be a powerful strategy to further disentangle different contributions in the angular PS (this will be pursued in a future work; S. Camera et al. 2013, in preparation). The shear signal is stronger for larger DM masses. The same is also true for the γ-ray signal from DM and this fact gives a large one-halo contribution which dominates starting from k ≲ 1 h Mpc⁻¹ in Figure 4 (right). Galaxies have masses ≲ 10¹⁴ M_☉, thus they correlate with the shear signal of lower-mass halos and the one-halo contribution becomes important at slightly smaller scale k ≳ 1 h/Mpc⁻¹. Since the bulk of unresolved blazars in 5 yr Fermi-LAT will be hosted in relatively small halos at large redshift, the one-halo term of the blazar/shear PS is suppressed. Thus, an important result is that, since both the shear and DM-induced γ-ray signals are stronger for larger halos, their cross-correlation is more effective with respect to the case of astrophysical sources. This, together with the sizable overlapping of the DM γ-ray and shear window functions at low redshift, leads to the expectation of a sizable DM signal in the angular PS, which is indeed what we find in Figure 5. For ℓ ≲ 100, the two-halo term dominates for all the sources, thus the relative size is roughly given by the relative contribution in the total EGB emission. At ℓ ≳ 100, the one-halo term starts to be important in the DM case which grows more rapidly than the astrophysical sources. At ℓ ≳ 10³, the one-halo term also takes over in the SFG spectrum which is brought again close to the DM curve. Blazars are largely subdominant in the whole range of multipoles.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The observational forecasts for the cross-correlation between DES or Euclid and Fermi-LAT are shown for the benchmark models considered in this work (for error estimates, we take observational performances from Atwood et al. (2009), The Dark Energy Survey Collaboration (2005), and Laureijs et al. (2011)). Figure 5 shows that a DM signal can be disentangled in the angular PS at ℓ ≲ 10³. The same conclusion can be derived for DM models with different mass and annihilation/decay channels, provided the DM is a significant component of the total γ-ray EGB (at least in one energy bin) as in our assumptions.

In this Letter, we discussed the cross-correlation angular PS of weak-lensing cosmic shear and γ-rays produced by WIMP annihilations/decays and astrophysical sources. We showed that this method can provide novel information on the composition of the EGB. Since the shear signal is stronger for structures of larger masses and most of the γ-ray emission from decaying and annihilating DM is also produced in large mass halos, their cross-correlation is typically stronger than the case of astrophysical sources (which are associated with galactic-mass halos). The combination of Fermi-LAT with forthcoming surveys like DES and Euclid can thus potentially provide evidence for WIMPs.

We thank L. Amendola, A. Cuoco, and A. Green for useful comments. N.F. and M.R. acknowledge INFN grant FA51. N.F. acknowledges support of Spanish grant MULTIDARK CSD2009-00064. M.F. is supported by a Leverhulme Trust grant. S.C. acknowledges support from FCT-Portugal under grant PTDC/FIS/100170/2008. S.C. work is funded by FCT-Portugal under Post-Doctoral Grant SFRH/BPD/80274/2011.