Mapping extragalactic dark matter structures through gamma-rays

Mapping extragalactic dark matter structures through gamma-rays


If dark matter is composed of neutralinos, the gamma-ray radiation produced in their annihilation offers an attractive possibility for dark matter detection. This process may contribute significantly to the extragalactic gamma-ray background (EGB) radiation, which is being measured by the FERMI satellite with unprecedented sensitivity. Using the high-resolution Millennium-II simulation of cosmic structure formation we have produced the first full-sky maps of the expected contribution of dark matter annihilation to the EGB radiation. Our maps include a proper normalization of the signal according to a specific supersymmetric model based on minimal supergravity. The new simulated maps allow a study of the angular power spectrum of the gamma-ray background from dark matter annihilation, which has distinctive features associated with the nature of the annihilation process. Our results are in broad agreement with analytic models for the gamma-ray background, but they also include higher-order correlations not readily accessible in analytic calculations and, in addition, provide detailed spectral information for each pixel. In particular, we find that color maps combining different energies can reveal the cosmic large-scale structure at low and intermediate redshifts.

I Introduction

Although dark matter accounts for most of the matter in the Universe, its nature remains unknown and its presence has only been inferred through its gravitational effects. However, if dark matter is made of neutralinos, a new particle predicted by Supersymmetry, it would also interact (although very weakly) with ordinary matter, and it might be detected soon in laboratories on Earth.

In addition, neutralinos can self-annihilate to produce ordinary particles like positrons Baltz and Edsjö (1999), neutrinos Berezinsky et al. (1996) and gamma-ray photons. If these byproducts of the annihilation are copious enough, they could be detected soon. In the present work, we will focus on the gamma-ray photons as a residual of the annihilation process.

This gamma-ray radiation is expected to be produced more abundantly in regions with high dark matter density. Thus, it seems best to look for it in very dense nearby regions, such as the center of our own Galaxy Bertone and Merritt (2005); Jacholkowska et al. (2006) and/or the centers of its satellite galaxies Wood et al. (2008). Actually, it turns out that the best prospects for the detection of gamma-rays from our Galaxy are obtained by looking slightly off-center to avoid confusion of the signal with other sources of gamma rays residing at the Galactic center Stoehr et al. (2003); Springel et al. (2008).

Outside our Galactic halo, gamma-rays are also produced in large quantities by the annihilation of dark matter in all the many halos and subhalos within our past light-cone, contributing to the so-called extragalactic gamma-ray background (EGB) radiation. The EGB has been measured by different satellites, in particular in the energy range between and by COMPTEL and EGRET Strong et al. (2004). Our understanding of the EGB is greatly improving by the Large Area Telescope aboard the recently launched FERMI satellite Atwood et al. (2009), which covers an energy range between and and features an improved sensitivity compared with its predecessor EGRET.

Although the EGB also receives contributions from other sources, such as blazars Ando et al. (2007a) and cosmic rays accelerated at structure formation shocks Jubelgas et al. (2008), the energy spectrum and angular power spectrum of the annihilation radiation have distinctive features that may open up effective ways for disentangling the signal Ando and Komatsu (2006); Siegal-Gaskins and Pavlidou (2009). Therefore, a detailed analysis of the EGB is a viable possibility for detecting dark matter.

Previous studies have analyzed this possibility using analytic approaches. However, it is not clear how accurately these methods capture the non-linear structures resolved in the newest generation of high-resolution cosmological simulations. Also, the previous analysis have so far been restricted to statistical statements about the power spectrum of the EGB, or its mean flux. For a full characterization of the signal, it would be useful to have accurate realizations of maps of the expected gamma-ray emission over the whole sky.

In our study Zavala et al. (2009), we therefore focus on predicting the EGB radiation directly from cosmological N-body simulations. To this end, we use the Millennium-II simulation (MS-II) Boylan-Kolchin et al. (2009). With particles in a homogeneously sampled volume of , it is one of the best resolved structure formation simulations to date. We use a ray-tracing technique to accumulate the signal from all halos and subhalos over the past light-cone of a fiducial observer, positioned at a plausible location of the Milky Way in the simulation box.

The values we calculate for each pixel of such a sky map are that of the specific intensity, the energy of photons received per unit area, time, solid angle and energy range:


where the integral is over the whole line of sight, is the comoving distance and is the energy measured by the observer at . The gamma ray emissivity (energy emitted per unit energy range, unit volume and unit time) associated with dark matter annihilation is given by:


where is the total differential photon spectrum summed over all channels of annihilation, and are the mass and density of neutralinos, and is the thermally averaged product of the annihilation cross section and the Møller velocity. Note that is evaluated at the blueshifted energy along the line-of-sight to compensate for the cosmological redshifting.

Ii Dark matter annihilation and the SUSY factor

Of all the quantities in Eq. (2), only depends on the spatial distribution of neutralinos, the rest is related to its intrinsic properties, which we can conveniently analyze separately by defining the so called SUSY factor:


In the present section we describe the particular SUSY model that we used to compute the gamma-ray emissivity produced by dark matter annihilation.

We restrict our analysis to the framework of the minimal supergravity model (mSUGRA), which is popular, among other reasons, due to its relative simplicity. The large number of free parameters in general SUSY is reduced to effectively four free parameters in mSUGRA: , , , which are the values of the gaugino and scalar masses and the trilinear coupling, all specified at the GUT scale, tan which is the ratio of the expectation values in the vacuum of the two neutral SUSY Higgs, and finally, the sign of , the Higgsino mass parameter.

The general 5-dimensional parameter space of the mSUGRA model is significantly constrained by various requirements: consistency with radiative electroweak symmetry breaking (EWSB) and experimental constraints on the low energy region. It is also normally assumed that if the major component of dark matter is the lightest neutralino, then its relic density should be equal to the observed abundance of dark matter . This condition reduces the allowed regions in the parameter space considerably.

The different regions in the mSUGRA parameter space that satisfy these constraints have been studied often in the past and have received generic names: (B) bulk region (low values of and ), (FP) focus point region (large values of ), (CA) co-annihilation region (low and , where is the lightest slepton) and (RAF) rapid annihilation funnel region (for large values of tan and a specific condition for 1). Instead of making a full scan of the allowed parameter space, we define several “benchmarks points” which fulfill all the experimental constraints discussed above and which can be taken as representative of the different regions that we just described. The benchmark points were selected following the analysis of Battaglia et al. (2004) and Gondolo et al. (2004). We use the numerical code DarkSUSY Gondolo et al. (2004, 2005) to analyze these benchmark points.

Figure 1: SUSY factor as a function of -ray emission energy for a selection of benchmark points representatives of the allowed regions in the mSUGRA parameter space: CA (black), FP (red), RAF (green ) and B (blue).

Fig. 1 shows the energy spectrum of computed for four of the benchmark points analyzed, which are representatives of the different regions: CA (black), FP (red), RAF (green ) and B (blue). The main features of Fig. 1 are determined by the photon spectrum which receives contributions from three mechanisms of photon production Bringmann et al. (2008): (i) gamma-ray continuum emission following the decay of neutral pions produced during the hadronization of the primary annihilation products; (ii) monoenergetic gamma-ray lines for neutralino annihilation in two-body final states containing photons; (iii) internal bremsstrahlung (IB), which leads to the emission of an additional photon in the final state; the IB contribution to the spectrum is typically dominant at high energies resulting in a characteristic bump. Processes (ii) and (iii) are subdominant to process (i), but they display distinctive spectral features intrinsic to the phenomenon of annihilation.

In Fig. 1, we have shown only 4 of all the benchmark points we have considered; the rest have very similar features and lie in between these 4 cases. Fig. 1 also indicates that the normalization of the spectrum of typically varies by three orders of magnitude among the different regions of the allowed parameter space.

For the rest of our analysis, we will use one particular benchmark point, which is the one we highlight with a thick blue line in Fig. 1. It is a model with  GeV that gives an upper limit on for the benchmark points in the bulk region and is also close to the maximum value of for all the benchmark points analyzed. The latter is the main reason to choose this particular model because it is desirable that the predicted gamma-ray flux has a large value among the different theoretical possibilities. There are other reasons as well. The mass of the neutralino in this model is “safely” larger than the lower mass bounds coming from experimental constraints ( according to Heister et al. (2004)), but low enough to be detectable in the experiments available in the near future. Also, a high value of the parameter tan seems to be favored by other theoretical expectations Núñez et al. (2008).

Iii The dark matter annihilation contribution to the EGB

The specific intensity (Eq. 1) describes the total emission from dark matter annihilations integrated over the full backwards light cone along a certain direction. We use the data of the MS-II to fill the whole volume contained in the past-light cone of an observer located at a fiducial position at . Since there are only 68 simulation outputs in total, we can approximate the temporal evolution of structure growth by using at each redshift along the past light-cone the output time closest to this epoch. To cover all space, we use periodic replication of the simulation box.

Our strategy is not to account for the emission at the level of individual dark matter particles, but instead to use entire halos and subhalos, allowing us to accurately correct for resolution effects. The list of dark matter particles at each output time is hence replaced by a catalog of dark matter substructures.

We further assume that these substructures are well represented by a spherically symmetric NFW density profile Navarro et al. (1997). Under this assumption the gamma-ray luminosity from a given substructure can be written as a scaling law Springel et al. (2008):


where is the maximum rotational velocity of the subhalo and is the radius where this maximum is reached. The most recent high-resolution N-body simulations support a slightly revised density profile that becomes gradually shallower towards the center, the so called Einasto profile. For this case, the coefficient in Eq. (4) changes to 1.87, increasing our results effectively by 50. This is an insignificant difference compared to other uncertainties in the analysis of the dark matter annihilation radiation. For the purposes of our work, the NFW profile is a good approximation and using it simplifies comparisons to a number of results in the literature.

Our map-making method then becomes a task to accurately accumulate the properly redshifted emission from these structures in discretized representations of the sky and to include a resolution correction for poorly resolved objects. All our maps use pixels, corresponding to an angular resolution of .

A given pixel in the simulated maps covers a solid angle . Our map-making code computes the average value of the specific intensity within the area subtended by this solid angle by conservatively distributing the emission of each substructure over the appropriate pixels. Combining Eqs. (24) we can write the relevant sum over all substructures in the light-cone as:


where the function is a weight function that distributes the luminosity of a given halo onto the pixels overlapping with the projected “size”of the halo; the latter depends on the distance of the observer to the halo and the transverse distance between the halo center and the center of the pixels it touches. Except for structures that are very nearby, the high central concentration of the emission of a subhalo and the limited angular resolution of our maps give most subhalos the character of unresolved point sources.

For the purpose of extending the predictions of our maps down to the damping scale limit of neutralinos (M), we divide the maps into separate components: (i) the contribution of resolved halos and subhalos with a minimum mass of (called “ReHS” in the following); (ii) the contribution of unresolved main halos with masses in the range (referred to as “UnH”); and finally, (iii) the contribution of unresolved subhalos in the same mass range as in the case (ii) (component “UnS”).

iii.1 Resolved halos and subhalos (ReHS)

In Fig. 2, we show full sky maps of the the -ray emission of all main halos and subhalos that are detected in the MS-II; this corresponds to an essentially perfectly complete sample of all halos above a mass limit of . In the top panel of Fig. 2, a partial map at low redshift (corresponding to the first shell at ) is shown in order to illustrate typical foreground structure, whereas the bottom panel gives an integrated map out to , which is approximately the full EGB from annihilation. The maps in Fig. 2 are for . In the color scale used for the maps, red corresponds to the highest and black to the lowest values of the specific intensity.

Figure 2: Upper panel: A partial map showing the EGB produced by dark matter annihilation in nearby structures. Only sources within 68 Mpc from the observer are considered for the map. The color scale gives a visual impression of the values of the specific gamma-ray intensity for each pixel in the map; the red color corresponds to the highest values of specific intensity. The observed energy of the simulated gamma-ray radiation is 10 GeV. Lower panel: The full gamma-ray sky map from dark matter annihilation containing sources up to . In both maps only the contribution of the smallest halos resolved by the simulation () has been taken into account.

iii.2 Unresolved halos (UnH)

The maps shown in Fig. 2 are complete down to the minimum mass in the MS-II that we can trust, . In order to make a prediction of the full EGB coming from dark matter annihilations, we extrapolate the -ray flux to account for the contribution of all missing dark matter halos down to the cutoff mass . For this purpose we use extrapolations of the power law behavior found for the total gamma-ray luminosity of host halos as a function of mass. This is an analysis that we present in detail in section 4.1 of Zavala et al. (2009).

The way we incorporate this extrapolated contribution in the -ray maps is the following. We assume that the EGB radiation from the missing halos in the mass range to is distributed on the sky in the same way as the one from the smallest masses we can resolve in the simulation, which we adopt as the mass range between and . Hence, we compute the value of a boost factor with which each resolved halo in the mass range needs to be multiplied such that the luminosity of the unresolved main halos is accounted for as well. We found that , and that its value is nearly independent of redshift up to the highest redshift we can reliable make the extrapolation, which is .

iii.3 Unresolved subhalos (UnS)

To add the contribution of unresolved substructures to the -ray maps, we first assume that the identified subhalo population is complete down to a mass limit of . Below this mass we compute the gamma-ray luminosity of subhalos using extrapolations of the power law behavior of the contribution of these subhalos to the luminosities of their host halos. The parameters of the power law are given in detail in section 4.2 of Zavala et al. (2009). We add this contribution in different ways for the following two cases.

In the case of resolved main halos with more than 100 particles, we distributed the extra luminosity coming from unresolved substructures among the subhalos of the main host, assuming in this way that unresolved subhalos are distributed in the same way as the resolved ones. If the main halos have no subhalos, the extra luminosity is given directly to them.

For the second case of all main halos with masses less than , we use the following strategy. From the results described in the last subsection we can get the total luminosity coming from all main halos with masses between the damping scale limit and . Thus, we can compute the boost factor for which this luminosity needs to be multiplied in order to include the full contribution of all subhalos of the main halos in this mass range. We found that is roughly independent of redshift. Using the range of values found for the power laws described before, we obtain that . We take the extreme values of this range in the following results as the minimum and maximum values that reflect the uncertainties in our extrapolation method and should bracket the true result.

iii.4 Isotropic and anisotropic components of the EGB from dark matter annihilation

The upper panel of Fig. 3 shows the mean value of the specific intensity divided by the comoving thickness of the shell at as a function of redshift. Recall that each partial map represents the total specific intensity in a shell of constant comoving thickness, so this provides information about where the signal is coming from. The black line is for the ReHS case, while the dashed line is for main halos only (resolved and unresolved, including the boost factor ). The remaining lines are for the resolved and unresolved main halos and their subhalos, boosted to include the contribution of structures all the way down to the damping scale limit; here the blue line is for the “minimum boost” with whereas the red line is for the “maximum boost” with , but both include the boost for unresolved main halos.

In all cases, the contribution from partial maps up until is almost constant, indicating that the increasing number of sources for shells at higher redshifts, due to the larger volume seen behind each pixel, is approximately compensated by the distance factor and the spectral effects from the cosmological redshifting and the intrinsic variation of the emission spectrum, as described by (see Fig. 1).

Figure 3: Upper panel: Mean annihilation intensity per comoving shell thickness, as a function of redshift for individual partial maps. The solid black line is for resolved halos and subhalos. The dashed black line is for main halos only, resolved and unresolved down to the cutoff mass using a boost factor . The blue and red solid lines show the contribution of all components, resolved and unresolved halos and subhalos, boosted with the extreme values in the interval . All lines are for , except for the black dotted line which is for the ReHS component with . Lower panel: The same as the upper panel but for the accumulated intensity, .

The lower panel of Fig. 3 shows the accumulated mean value of the specific intensity as a function of redshift. This clearly shows that the most relevant contributions come from maps up to . Overall, applying the maximum boost increases the mean value of of the ReHS case by three orders of magnitude. For most energies, the spectrum is monotonically decreasing. In these regions, the dependence on with redshift will have the generic form shown in Fig. 3. However, at the highest energies the shape is expected to change dramatically due to the importance of Internal Bremsstrahlung and/or of monochromatic lines (see Fig. 1) close to the rest mass energy of the dark matter particle. Such an effect is clearly visible in the black dotted line in Fig. 3, which is the ReHS case for an energy of and features a bump in the redshift distribution, a reflection of the importance of IB at the highest energies.

The angular power spectrum of the EGB is an important tool to study the statistical properties of its anisotropy on the sky. In fact, certain classes of -ray sources are expected to exhibit different power spectra, making this a potential means to identify the origin of the unresolved EGB.

Figure 4: Angular power spectrum of the EGB produced by dark matter annihilation as a function of the multipole for . The line styles are as in Fig. 3 except for the dotted lines that show the predictions from Ando et al. (2007b) for a subhalo-dominant contribution with and without considering tidal destruction (lower and upper dotted lines), respectively.

In Fig. 4, we show the angular power spectrum as a function of the multipole , at an observed energy of , for all the cases discussed in Fig. 3. At large scales, , the power spectrum is related to the clustering of dark matter halos. When only the main halo contribution to the EGB is considered (dashed black line), the normalization is lower because most of the -ray signal comes from low mass halos that are less clustered than more massive halos.

The blue line in Fig. 4, corresponding to the full extrapolation including subhalos but with the minimal boost , has exactly the same power spectrum than the case with main halos only, at all scales. This is because the signal from main halos is dominant in this case, and subhalos have a negligible effect. In contrast, in the case where subhalos have a significant contribution mediated by (red line), the normalization is larger at small scales because the signal is dominated by subhalos belonging to the most massive halos, and the latter are strongly clustered.

As the angular scale decreases, , the power spectrum depends more and more on the internal structures of halos. For the cases where substructures are ignored or are negligible (dashed-black and blue lines), the slope becomes steeper, with a slope close to 2. The power spectrum for the cases where substructures are relevant (red and black solid lines) behaves differently, however. In the range , it becomes slightly shallower, i.e. the signal is slightly more isotropic in this regime. Contrary to the strong central concentration of the matter in a halo, the number density profile of subhalos is considerably shallower than a NFW profile and produces a luminosity profile in projection which is essentially flat Springel et al. (2008). This effect continues until where the power spectrum becomes dominated by the low-mass main halos.

The dotted lines in Fig. 4 were taken from the analytical results of Ando et al. (2007b) and are shown here for comparison, see the las paragraphs of section 5.5 in Zavala et al. (2009) for a discussion on the differences between our results and those of Ando et al. (2007b).

Figure 5: Upper panel: Full-sky maps at energies 0.1 GeV and 32 GeV in the left and right, respectively. The maps were smoothed with a Gaussian beam with a FWHM of . At a single energy, a full-sky map is very smooth, nearby structures are only minimally visible. Lower panel: Ratio of the maps in the upper panel (left) and a partial map containing only nearby structures within 68 Mpc for an observed energy of 0.1 GeV (right). Creating difference maps (“color” maps) using different energy channels greatly enhances the signal of nearby structures.

iii.5 Energy dependence

An important result that we obtain from the analysis of our results is that there are differences in the shape of the power spectra for different energies. We can exploit these differences and enhance the power at small or large angular scales by taking the ratio of sky maps at different energies. The upper panel of Fig. 5 shows the total sky maps for and in the left and right, respectively, smoothed with a Gaussian beam with a FWHM of . Both maps show only small anisotropies, with the exception of a couple of prominent structures. The sky map in the lower left panel is the ratio of the two maps; it clearly enhances the signal of nearby structures. Convincing supporting evidence for this enhancement comes from the nearest partial map at for shown in the lower right corner of Fig. 5. Most of the prominent dark matter structures that can be seen in this map are also clearly present in the map with the energy ratio. We conclude that the -ray sky maps that the FERMI satellite will obtain at different energies could be used to construct difference maps which may then show enhanced correlations with nearby cosmic large-scale structure.


  1. The relation which should hold is: , where is the CP-odd Higgs boson; see for example Feng (2005) for details.


  1. E. A. Baltz and J. Edsjö, Phys. Rev. D 59, 023511 (1999).
  2. V. Berezinsky et al., Astroparticle Physics 5, 333 (1996).
  3. G. Bertone and D. Merritt, Modern Physics Letters A 20, 1021 (2005).
  4. A. Jacholkowska et al., Phys. Rev. D 74, 023518 (2006).
  5. M. Wood et al., ApJ 678, 594 (2008).
  6. F. Stoehr et al., MNRAS 345, 1313 (2003).
  7. V. Springel et al., Nature 456, 73 (2008).
  8. A. W. Strong, I. V. Moskalenko, and O. Reimer, ApJ 613, 956 (2004).
  9. W. B. Atwood et al., ApJ 697, 1071 (2009).
  10. S. Ando et al., MNRAS 376, 1635 (2007a).
  11. M. Jubelgas et al., A&A 481, 33 (2008).
  12. S. Ando and E. Komatsu, Phys. Rev. D 73, 023521 (2006).
  13. J. M. Siegal-Gaskins and V. Pavlidou, Physical Review Letters 102, 241301 (2009).
  14. J. Zavala, V. Springel, and M. Boylan-Kolchin, ArXiv e-prints (2009), eprint arXiv:0908.2428.
  15. M. Boylan-Kolchin et al., MNRAS 398, 1150 (2009).
  16. J. L. Feng, Annals of Physics 315, 2 (2005).
  17. M. Battaglia et al., European Physical Journal C 33, 273 (2004).
  18. P. Gondolo et al., JCAP 7, 8 (2004).
  19. P. Gondolo et al., New Astronomy Review 49, 149 (2005).
  20. T. Bringmann, L. Bergström, and J. Edsjö, Journal of High Energy Physics 1, 49 (2008).
  21. A. Heister et al., Physics Letters B 583, 247 (2004).
  22. D. Núñez et al., JCAP 5, 3 (2008).
  23. J. F. Navarro, C. S. Frenk, and S. D. M. White, ApJ 490, 493 (1997).
  24. S. Ando et al., Phys. Rev. D 75, 063519 (2007b).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description