A Detailed modelling of starlight and dust in clusters

Prospects of detecting gamma-ray emission from galaxy clusters: cosmic rays and dark matter annihilations

Abstract

We study the possibility for detecting gamma-ray emission from galaxy clusters. We consider 1) leptophilic models of dark matter (DM) annihilation that include a Sommerfeld enhancement (SFE), 2) different representative benchmark models of supersymmetric DM, and 3) cosmic ray (CR) induced pion decay. Among all clusters/groups of a flux-limited X-ray sample, we predict Virgo, Fornax and M49 to be the brightest DM sources and find a particularly low CR-induced background for Fornax. For a minimum substructure mass given by the DM free-streaming scale, cluster halos maximize the substructure boost for which we find a factor of . Since regions around the virial radius dominate the annihilation flux of substructures, the resulting surface brightness profiles are almost flat. This makes it very challenging to detect this flux with imaging atmospheric Cherenkov telescopes since their sensitivity drops approximately linearly with radius and they typically have linear resolution elements across a cluster. Assuming cold dark matter with a substructure mass distribution down to an Earth mass and using extended Fermi upper limits, we rule out the leptophilic models in their present form in 28 clusters, and limit the boost from SFE in M49 and Fornax to be . This corresponds to a limit on SFE in the Milky Way of which is too small to account for the increasing positron fraction with energy as seen by PAMELA and challenges the DM interpretation. Alternatively, if SFE is realized in Nature, this would imply a limiting substructure mass of —a problem for structure formation in most particle physics models. Using individual cluster observations, it will be challenging for Fermi to constrain our selection of DM benchmark models without SFE. The Fermi upper limits are, however, closing in on our predictions for the CR flux using an analytic model based on cosmological hydrodynamical cluster simulations. We limit the CR-to-thermal pressure in nearby bright galaxy clusters of the Fermi sample to and in Norma and Coma to . Thus we will soon start to constrain the underlying CR physics such as shock acceleration efficiencies or CR transport properties.

pacs:
95.35.+d, 95.85.Pw, 98.62.Gq, 98.65.-r, 98.70.Sa

I Introduction

Dark matter (DM) has been searched for in direct detection experiments Pato et al. (2011), at accelerators Ellis et al. (2002); Baer et al. (2007); Khachatryan et al. (2011) and also in indirect detection experiments looking for signals in the cosmic-ray (CR) spectra of antiprotons, positrons, neutrinos and all of the electromagnetic spectrum from radio waves to gamma-rays Bergstrom (2009). So far, the improvements in direct detection sensitivity have put this method into focus, but the situation may change considerably the coming few years as the CERN LHC experiments collect data, and new gamma-ray detectors are being planned, such as the CTA CTA (2010). In fact, it has recently been pointed out Bergstrom et al. (2011) that a dedicated ground-based gamma-ray detector would have potential that goes far beyond that of the other methods, depending on presently unknown parameters in the particle physics models for DM.

Among the astrophysical systems which will be very interesting to detect, and study, with gamma-ray detectors (Fermi, H.E.S.S., MAGIC, VERITAS, and eventually large detectors like CTA)1 belong galaxy clusters. The most promising directions in which to search for a gamma-ray annihilation signal (from the annihilation process itself, and also the accompanying bremsstrahlung and inverse Compton (IC) components coming from charged particles produced in the annihilations) are basically three:

1. The galactic center (g.c.). This is where all numerical simulations of cold dark matter (CDM) predict the highest density. However, the detailed DM density in the very central part is difficult to predict, due a to a possibly very complicated interplay between baryons, DM, and the central galactic black hole. Also, it is a very crowded region with many gamma-ray sources like pulsars, CR-illuminated molecular clouds, and other supernova remnants, which have to be subtracted from the data to extract the DM induced signal. In fact, there is a recent claim of an indication of a relatively light DM particle contribution to the gamma-ray flux from the g.c. Hooper and Goodenough (2011), but other hypotheses seem to work at least as well Boyarsky et al. (2010).

2. The dwarf spheroidal galaxies orbiting the Milky Way (MW), like Segue-1, Ursa Minor, Draco, Sagittarius, Sculptor, Carina or Willman-1 Bringmann et al. (2009); Acciari et al. (2010); Scott et al. (2010, 2010); Aleksić et al. (2011); Abramowski et al. (2011). The problem here is that the nature of many of these small, dark matter-dominated galaxies is not entirely clear, and the velocity dispersion estimates are based on rather small numbers of stars. Confusion with star clusters and tidal disruption are other complications. Once a satellite dwarf galaxy is accreted by the MW, the outer regions are severely affected by tidal stripping. The longer a satellite has been part of our Galaxy, and the closer it comes to the center during its pericentral passage, the more material is removed Gao et al. (2004). Thus the DM density profile is very uncertain for most of them, especially for radii larger than those probed by the stars. Nonetheless, by stacking the data together from many dwarf spheroidals these uncertainties can be made less severe, and preliminary results from Fermi-LAT shows this method to give quite promising results Garde ().

3. Galaxy clusters. This possibility has been less studied theoretically, however currently there is an ongoing observational campaign to detect gamma-ray emission from galaxy clusters Reimer et al. (2003); Perkins et al. (2006); Perkins (2008); Aharonian et al. (2009a); Galante and for the VERITAS Collaboration (2009); Domainko et al. (2009); Aharonian et al. (2009a, b); Kiuchi et al. (2009); Acciari et al. (2009); Aleksić et al. (2010a); Ackermann et al. (2010, 2010). In fact, we noted in a previous Letter Pinzke et al. (2009) that there are certain advantages that work in favour of this possible target for gamma-ray detection of DM annihilation. Galaxy clusters constitute the most massive objects in our Universe that are forming today. This causes their DM subhalo mass function to be less affected by tidal stripping compared to galaxy sized halos that formed long ago. The annihilation luminosities of the DM halo component for e.g. the Virgo cluster and the Draco dwarf scales in a way (see Pinzke et al. (2009)) that the ratio of gamma-ray luminosities from the smooth components is around 4, in favour of Virgo. In addition, there may be a further enhancement due to substructure, which to a large extent should be unaffected by tidal disruption, at least in the outer regions. According to a recent estimate Gao et al. (2011a), more massive haloes tend to have a larger mass fraction in subhalos. For example, cluster size haloes typically have 7.5 per cent of the mass within in substructures of fractional mass larger than , which is 25 per cent higher than for galactic haloes.2

In this paper, we will investigate in detail the potential of several of the most promising galaxy clusters to produce an annihilation gamma-ray yield which could be observable with present and planned gamma-ray detectors. Here, we sketch out the main arguments that need to be considered for maximizing the expected signal-to-noise ratio of a promising target cluster. First, we need to maximize the DM annihilation flux of an unresolved cluster, , where depends on the particle physics model of DM, is the smooth DM density profile, is the luminosity distance, and if we assume universality of the DM density profile. The dependence of on the halo formation epoch breaks the universality and slightly modifies the mass dependence, yielding Pinzke et al. (2009). Additionally, the presence of a hierarchy of substructures down to small scales Springel et al. (2008a) and the potential dependence of the particle physics cross section on the relative DM particle velocities Arkani-Hamed et al. (2009) may furthermore modify the scaling parameter and shall be one of the focus points of this work. Second, we need to minimize the expected noise which is a sum of instrumental noise, galactic and cluster-intrinsic foreground. While the galactic foreground varies across the sky, it is typically lower for increasing latitude (away from the galactic plane). It is thought that the gamma-ray foreground from clusters is dominated by CRs which are accelerated at cosmological formation shocks and transported over time although there could be a substantial contribution of CRs from AGNs and supernova-driven galactic winds. In this work we will also study (and constrain) the CR-induced emission from clusters to identify those objects that are expected to be especially dim which would imply a low CR-induced background. Since the CR-induced emission is expected to scale with the thermal X-ray emission of clusters Pfrommer and Enßlin (2004a), we estimate gamma-ray fluxes of a nearby X-ray flux complete sample of galaxy clusters from the Rosat all-sky survey (extended HIghest X-ray FLUx Galaxy Cluster Sample, HIFLUGCS, Reiprich and Böhringer (2002); Chen et al. (2007)).3 This complements and extends previous work related to DM in clusters Colafrancesco et al. (2006); Jeltema et al. (2009); Sanchez-Conde et al. (2011); Cuesta et al. (2011); Huang et al. (2011) and to CRs in clusters Pinzke and Pfrommer (2010); Brunetti and Lazarian (2011); Pfrommer et al. (2008); Kushnir and Waxman (2009); Donnert et al. (2010a).

In Section 2 we discuss the theory of DM and gamma-rays. In particular, we focus on the leptophilic (LP) and supersymmetric benchmark (BM) DM models in this work, and outline the framework for estimating the gamma-ray emission from various radiative processes. In Section 3 we calculate the gamma-ray fluxes and spectral distributions of four clusters that are identified as prime targets for DM observations (Fornax and Virgo) as well as to have a high CR-induced gamma-ray yield (Perseus and Coma). For the same clusters we derive the gamma-ray surface brightness profiles in Section 4. To extend the analysis to more clusters and increase the list of promising targets we estimate the fluxes from DM and CRs of all clusters in the extended HIFLUGCS sample in Section 5. In Section 6 we conclude and discuss our results.

Ii Theory

We start the section by discussing two different, but well motivated DM models; LP and supersymmetric DM. We then present the framework that is used to calculate the gamma-ray emission for these DM models using an Einasto DM density profile and the expected enhancement from DM substructures. This is followed by an outline of the framework of IC emission where we take into account photon fields of the cosmic microwave background (CMB), dust, and starlight. We end by summarizing our formalism of calculating the CR-induced pion decay gamma-ray emission which is thought to dominate the astrophysical gamma-ray signal.

ii.1 Detecting Particle Dark Matter

Besides the intrinsic interest in the gamma-ray flux from galaxy clusters generated by conventional hadronic and electromagnetic processes - which should be close to observability with the Fermi-LAT data Berezinsky et al. (1997); Enßlin et al. (2007); Pinzke and Pfrommer (2010) - the possible contribution from WIMP DM is of greatest interest. A WIMP (Weakly Interactive Massive Particle) which fulfils the WMAP bounds on the relic density of CDM of Komatsu et al. (2011)

will in many cases naturally give a gamma-yield which may be observable. In addition, there are possible enhancement effects known, such as the astrophysical boost from dense substructure of DM halos, or the particle physics boost from the Sommerfeld effect, which may increase the chances of detection further.

There are three methods for detecting DM candidates which are presently employed: First, at particle physics accelerators such as CERN’s LHC, one is now entering into a new energy regime which may allow the production of the heavy, electrically neutral and long-lived particles which may constitute DM. From these experiments one may get the first glimpse of the mass scale beyond the standard model where DM may reside. However, to really show that any of the hypothetical, new particles created at the LHC is actually the DM, one has to rely on the two other methods available for the DM search, namely direct and indirect detection. Second, direct detection methods, which are presently evolving rapidly, use the feeble interaction between DM particles of the galactic halo and nuclei such as Germanium, Sodium, Iodine, Argon or Xenon to infer the scattering cross section and a rough estimate of the mass of the DM particle (for the currently most sensitive search, see Aprile et al. (2010, 2011)). A characteristic of this method is that it basically only depends on the local DM density (which is rather well determined to be around GeV/cm) and the basic cross section for DM - nucleus scattering. A disadvantage is that it does not benefit from the two enhancement mechanisms mentioned above, boost from substructure and/or the Sommerfeld effect. Also, it cannot be excluded - although it seems at present improbable - that the local halo structure is such that the solar system happens to be in an underdense region.

Third, in indirect detection, in particular in the photon channel, one searches for products of DM annihilation in the galactic halo and beyond. Particularly interesting targets are the galactic center (for a recent possible indication of a signal, see Hooper and Goodenough (2011), however, see also Boyarsky et al. (2010)), the dwarf galaxies surrounding the MW Strigari et al. (2007); Essig et al. (2009); Scott et al. (2010), galaxy clusters (the focus of this work) Ghigna et al. (1998); Lewis et al. (2003); Boyarsky et al. (2008); Colafrancesco et al. (2006); Pinzke et al. (2009), and even the cosmological large scale structure Bergström et al. (2001); Ullio et al. (2002); Taylor and Silk (2003); Elsaesser and Mannheim (2005); Cuoco et al. (2011); Abazajian et al. (2010); Abdo et al. (2010); Zavala et al. (2011).

Leptophilic models

There is also a possibility to search for DM annihilation in the Milky Way indirectly through annihilation to antimatter channels, which however lacks the important directional signature of gamma-rays. Recently, it has however been much in focus due to the surprising findings of PAMELA Adriani et al. (2009) and Fermi Abdo et al. (2009). Viable particle physics models are rather constrained by other observations, however. For example, the non-observation of an enhanced antiproton flux by PAMELA Adriani et al. (2010) means that quark-antiquark final states have to be suppressed. This has led to the postulate of a whole class of leptophilic models with mainly or exclusively annihilation to leptons and antileptons. In addition, one has to have a sizeable fraction of final states containing muons or tau leptons as the shape of the spectrum disfavors direct annihilation to electrons and positrons only. A large enhancement, of the order of at least several hundred, of the annihilation rate is also needed, something that may be given by the Sommerfeld effect (for a pedagogical review of viable scenarios, see Arkani-Hamed et al. (2009)). As still another difficulty for these models, the lack of an IC signature from the Galactic center means that the DM density distribution must be cored rather than cuspy Bertone et al. (2009); Cirelli et al. (2009); Bergström et al. (2009a). (For a review of the DM modeling of these effects, see Bergstrom (2009). For a recent treatment, showing still viable models, see Finkbeiner et al. (2011).)

A more standard explanation of the PAMELA/Fermi excess would be pulsars or other supernova remnants (e.g., Atoyan et al. (1995); Hooper et al. (2009); Ahlers et al. (2009)). In that case, DM would more likely be explained by the conventional scenario of a WIMP, for example the thoroughly studied lightest supersymmetric neutralino (for reviews, see Jungman et al. (1996); Bergström (2000); Bertone et al. (2005)). We note that AMS-02 experiment on the International Space Station Ting et al. (2011) may give interesting clues to the origin of the antimatter excess.

In the scenario with leptophilic DM, one has to rely on some enhancement of the annihilation rate from the effects of DM halo substructure, a possibility which was realized long ago Silk and Stebbins (1993); Bergström et al. (1998); Moore et al. (1999). However, a boost factor as large as several hundred is difficult to achieve in the solar neighborhood, due to tidal stripping of subhalos in the inner part of the Galaxy. Sommerfeld enhancement is thus important for leptophilic models that claim to explain the PAMELA/Fermi result. This effect was computed for electromagnetism by Arnold Sommerfeld many years ago Sommerfeld (1931) and recently rediscovered in the quantum field theory of DM Hisano et al. (2005); Cirelli et al. (2007); Arkani-Hamed et al. (2009). In quantum mechanics describing electron scattering and electron positron annihilation, it is caused by the distortion of the plane wave describing the relative motion of the annihilating particle pair, due to the near formation of a bound state caused by photon exchange. In the ladder approximation for QED, one reproduces the Sommerfeld effect, and the square of the relative wave function at the origin (which enters into the probability for the short-distance process of annihilation) is increased by the factor Arkani-Hamed et al. (2009)

(1)

with the fine-structure constant, and the relative velocity. This amounts to for small velocities. For a Yukawa-like particle of mass , mediating a weak attractive force with coupling constant between DM particles with mass , the small-velocity limit of the enhancement becomes instead

(2)

In the general case, with mediators that may also excite virtual charged particles in the DM sector, one has to solve numerically a coupled system of differential equations using appropriate boundary conditions Hisano et al. (2005); Cirelli et al. (2007); Arkani-Hamed et al. (2009). In some cases, the enhancement factor can be as high as several hundred to a few thousand, depending however on the exact parameters of the theory. The effect is usually strongly velocity-dependent, depending on velocity as or, in the fine-tuning case of being very near resonance, as . This means that in a virialised system (such as a galaxy cluster) with large velocity dispersion the SFE will be smaller than the one expected in a single galaxy such as the MW, roughly by a factor . Note that the -scaling is valid only for ; at smaller velocities and outside resonances, the -enhancement saturates at Kamionkowski and Profumo (2008).

As LP models, apart from being slightly contrived from the particle physics point of view, are also rather limited by several sets of astrophysical data. We take advantage of the recent reanalysis Finkbeiner et al. (2011) to define a benchmark LP model that is still viable. It is found that the bounds from WMAP5 approximately imply a Sommerfeld boost factor which has to satisfy GeV, where is the fraction of energy from DM annihilation that ionizes the intergalactic medium, for annihilation to electrons and a factor of a few smaller for all other standard model final states.

Taking into account the whole cosmic history of LP models, also utilizing limit on spectral and polarization distortions of the CMB, the maximal value for the boost factor for a 1 to 2 TeV particle and sub-GeV force carriers is found to be Finkbeiner et al. (2011) between 400 and 800. Although lower than the first estimates (e.g., Bergström et al. (2009b); Meade et al. (2010)) this is still enough to explain the PAMELA and Fermi excess, given other astrophysical uncertainties.

We adopt our benchmark LP model from Finkbeiner et al. (2011) where we use a DM mass  TeV and a branching ratio of () into (). Furthermore, since most of the gamma-ray flux is expected to come from dense subhalos within clusters which have a velocity dispersion close to the velocity limit where the SFE saturates, we use for the cluster halos.4 However, for figures where we neglect the substructure boost, we adopt a velocity dependent SFE that is normalized to fit the electron and positron excess observed at Earth,

(3)

where the local boost factor required to explain the data is . The velocity dispersion for the MW, , and the mass dependent velocity dispersion for clusters is derived in Voit (2005).

Supersymmetric dark matter

The most studied models for DM are supersymmetric ones, especially models where the lightest supersymmetric particle, in most models the lightest neutralino, is stable. The stability is assured in viable supersymmetric models due to a discrete symmetry, -parity. This symmetry is needed from the particle physics point of view to avoid fast proton decay, and automatically makes the lightest supersymmetric particle stable, and therefore a good candidate for DM. In addition, as the coupling to ordinary matter is given by gauge couplings, the neutralino is automatically a WIMP candidate. Unfortunately, the breaking of supersymmetry (which has to be there since otherwise, e.g., there would exist a scalar electron with the same mass as the electron) means that a large number of essentially free parameters enter actual calculations. The neutralino is then a linear combination of the supersymmetric partners of the U(1) gauge boson, , one of the SU(2) gauge bosons, , and the two Higgs doublets in the Minimal Supersymmetric extension of the Standard Model (the MSSM), and , with mixing coefficients which depend on the supersymmetric breaking parameters Ellis et al. (1984). The parameter space is very large, too large in general to scan efficiently even with the most powerful computers. Therefore, often simplifying assumptions are made to limit the number of free parameters. Even then, a large number of viable models is found which give a relic density consistent with the WMAP data. This means that they are indeed very good templates for WIMP DM, something that explains their popularity besides the subjective statement that few more well-motivated models for DM have been put forward so far.

The wave function of the lightest neutralino can be written

(4)

with

(5)

where the gaugino fraction of a given neutralino is whereas is the higgsino fraction.

We use the DarkSUSYpackage Gondolo et al. (2004) to compute the mass, couplings and relic density for a given set of parameters. Actually, we take advantage of so-called benchmark (BM) models, which have been proposed in different contexts (in particular in Bringmann et al. (2009) giving predictions for imaging air Cherenkov telescopes, like MAGIC II and CTA). They give a good representation of models with high enough gamma-ray rates to become the first candidates for DM indirect detection in imaging air Cherenkov telescopes. We use the following supersymmetric BM models:

  • : This model was introduced in Battaglia et al. (2004), where its phenomenology at colliders was studied. Annihilation directly into lepton pairs is suppressed for neutralinos due to their Majorana nature and therefore helicity suppression for annihilation in the galactic halo Goldberg (1983). The higher order process (inner bremsstrahlung, IB) , which does not suffer from helicity suppression Bergström (1989); Bringmann et al. (2008), gives a considerable contribution due to light sleptons in this model.

  • : This model, also from Battaglia et al. (2004) is in the so-called coannihilation tail. The sleptons are nearly degenerate with the neutralino, which causes the large IB from the leptonic final states to give a high enhancement of the flux.

  • : A representative model for the funnel region, where the annihilation dominantly occurs in the -channel through exchange of the pseudoscalar Higgs boson Battaglia et al. (2004). Here IB contributions are not important.

  • : Annihilation in the so-called coannihilation region, introduced in Bringmann et al. (2008) as BM3, with a particularly large IB contribution.

The details for the DM BM models are summarized in Table 1.

BM
140 0.09
315 0.12
570 0.10
234 0.09
Table 1: Relevant parameters for the four benchmark models. The mass of the DM particles are given denoted by , the relic density by and the annihilation rate today by .

Final state radiation

There are two types of radiative processes which are important for our BM models, first the usual QED final state radiation for both the LP model and for supersymmetric models with charged final states. In addition, there may be direct emission from a virtual, charged, exchanged particle (such as a spin-0 slepton or squark), internal bremsstrahlung, IB. The latter is essentially only important when there is a helicity suppression for the lowest order annihilation process, which is the case for neutralinos, whose Majorana nature make the annihilation rate in the wave proportional to , for final state fermion . Interestingly, both the gamma-ray and fermion energy spectra are peaked at the highest energy for this process, which can in some models cause a rather spectacular bump in the spectra for these models.

The photon spectrum resulting from final state radiation is universal with only a weak dependence of the underlying particle physics model. The photon yield from this process is given by (see e.g. Bringmann et al. (2008))

(6)

Here, the normalized photon energy , is the fine-structure constant, and the charge and mass of the final state particle , respectively. The function depends on the spin of the final state and is given by

(7)

for fermions.

The differential energy spectrum for the IB process is more complicated (see Bergström (1989); Bringmann et al. (2008)), but can be computed using DarkSUSY. The differential photon energy spectra for all our BM models are shown in Fig. 1. The peaking at high photon energy and the resulting bump caused by inner bremsstrahlung is clearly seen for BM models , , and . Despite the different particle masses and cross sections, the spectrum of the continuum emission shows a remarkably similar shape which suggests that a simple rescaling of particle masses (responsible for the exponential cutoff in the spectrum) and cross sections (spectral amplitude) should yield a roughly scale-invariant continuum emission spectrum except for the presence of the final state emission feature which depends on the details of the specific decay channels. The right panel of Fig. 1 shows the cumulative number of leptons (electrons and positrons) above a given energy for our different DM models which is of interested for the high-energy IC emission of those leptons. Only LP models have energetic enough electrons such that the IC emission is powerful enough to either be constrained or detected at GeV energies and higher.

Figure 1: Source functions for different dark matter (DM) models. We show the simulated data from DarkSUSY Gondolo et al. (2004), and that generated by Ciafaloni et al. (2011); Cirelli et al. (2011); the solid lines show the fit to the data. Left panel: normalized differential continuum spectra for four different DM benchmark (BM) models; model (red circles), (orange squares), (green diamonds), and (blue triangles). We use Eq. (62) to fit the continuum spectra. Right panel: number of electron and positron per DM annihilation above the electron energy for different DM models; BM model (red circles), BM model (orange squares), BM model (green diamonds), and BM model (blue triangles), leptophilic (LP) DM annihilating indirectly into electrons and positrons (purple ) and into muons (black ) without considering Sommerfeld boosts. We use Eq. (63) and Eqs. (LABEL:eq:lp_mu-65) to fit the spectra of electrons and positrons from BM and LP models, respectively.

ii.2 Astrophysical modeling of DM induced emission

We now turn to the detailed modeling of the surface brightness profiles of DM annihilation emission and discuss the DM profiles for the smooth distribution and substructures.

General equations

The differential photon flux within a given solid angle along a line-of-sight (los) is given by

(8)

where

(9)

and denotes the surface brightness above the photon energy . The integration along the line-of-sight , in the direction that the detector is pointing, is parametrized such that the radius of the source , where is the distance from the Earth to the center of the cluster halo and . The angular integration is performed over a cone centered around and the opening angle is typically taken to be a few times the point spread function (PSF) . The limited angular resolution results in a probability that a photon coming from a direction ’ is instead reconstructed to a direction , where the underlying probability distribution follow a Gaussian:

(10)

We denote the total differential source function by , where we include contributions from five main processes; leptophilic DM annihilating through (neglecting the produced neutrinos) where the pairs IC upscatter background photons (LP-IC), leptophilic DM emitting final state radiation (LP-FSR), supersymmetric DM BM models where annihilating neutralinos generate pairs that upscatter background photons (BM-IC) and emit a continuum as well as final state radiation (BM-Cont), and CR proton induced that decays into gamma-rays (CR-). The source function is given by

(11)

where the differential CR to gamma-ray source function is denoted by (see Sec. II.4 and Pinzke and Pfrommer (2010) for further details). The subscript runs over the gamma-ray producing DM channels and the total differential boost factor for DM is given by:

(12)

It is the product of enhancement factors from SFE (see Eq. 3 and Sec. II.1.1) and from substructure enhancement over the smooth halo contribution 5 (see Eqs. 26-29 and Sec. II.2.3). The DM source function from the smooth halo for each process is written in the form:

(13)

where the annihilation rate density is given by

(14)

Here, the subscript runs over all kinematically allowed gamma-ray producing channels, each with the spectrum and annihilation cross section . We denote the DM mass with and the smooth DM density profile with .

The smooth DM density profile

Typically the universal Navarro-Frenk-White (NFW) density profile provides a good fit to both the observed and simulated clusters. It can be considered as a special case of the more general 5-parameter profile:

(15)

Here, denotes the inner slope, is the outer slope, and is the shape parameter that determines the profile shape at the scaling radius that characterizes the transition between the different power-law slopes. A cuspy NFW profile is given by . The characteristic overdensity for an NFW profile is given by

(16)

where the halo mass dependent concentration parameter is derived from a power-law fit to cosmological simulations with Macciò et al. (2008),

(17)

This mass scaling agrees well with Zhao et al. (2009) for cluster-mass halos after converting the concentration definitions according to Hu and Kravtsov (2003). In this work we choose to model the DM density by an Einasto density profile

(18)

that is slightly shallower in the center than the conventional Navarro-Frenk-White (NFW) profile, but provides a better fit to recent simulated high resolution DM halos Merritt et al. (2006); Navarro et al. (2010). It should also be noted that recent observations of the Abell 383 galaxy cluster find a density profile with a shallower inner slope of compared to an NFW profile, and can be ruled out with % confidence Newman et al. (2011). These observations are based on lensing and X-ray measurements as well as the stellar velocity dispersion of the central galaxy. In Eq. (18) we denote the density where the profile has a slope of by , and the radius by . We use that the density and determine through :

(19)

Here for and since is a slowly increasing function with concentration , we fix it for simplicity to .

In the recent dark matter simulation literature, it has become standard to characterize halos by the value, , where the circular velocity attains its maximum:

Especially for subhalos, this quantity seems to be more stable compared to the virial mass that is subject to tidal stripping. For comparison, we quote typical values for a large galaxy cluster and a galaxy group:

(21)
(22)

Substructures

High-resolution dissipationless DM simulations of MW type halos find substantial amount of substructures in the periphery of DM halos, while the substructures in the center suffer from dynamical friction and tidal effects depleting their central number densities. Since the DM annihilation rate depends on the density squared, the resulting flux from substructures is boosted compared to the smooth density distribution. While there is still a discrepancy in the literature of the exact value of the predicted boost factor of the DM annihilation luminosity from substructures in DM halos, this inconsistency starts to become resolved, apparently converging toward predictions at the high end Springel et al. (2008b, a); Anderson et al. (2010). Following recent high-resolution simulations of Galaxy-sized halos, we adopt a boost factor due to substructures of 230 for such a halo Springel et al. (2008b) (which needs to be scaled to cluster halos as we will discuss below).

The initial suggestions of a small total boost () of order unity by the Via Lactea simulations Kuhlen et al. (2008) made assumptions of computing the boost assuming that the substructure luminosity follows that of the smooth DM distribution. Physically, this would imply substructure to follow the smooth DM distribution and a radially independent concentration parameter of substructures. Both assumptions are in conflict with tidal mass loss of satellites which is at work in simulations and cause the radial number density to be anti-biased with respect to the host’s mass density profile yielding a significantly flatter subhalo distribution compared to the smooth DM distribution Anderson et al. (2010). This increases the substructure boost preferentially in the outer parts of DM halos with a factor that ranges between 20 and 1000 (for dwarfs and galaxy clusters), hence depending sensitively on the host halo mass (see Figure 2). Recent Phoenix simulations by Gao et al. Gao et al. (2011b) confirm the large enhancement due to substructures in clusters. In fact they find a boost of 1125 for a Coma like cluster which agrees within 10% with our estimated boost (see e.g. Table 2).

We use a double power-law function to fit the luminosity from the smooth component of substructures (i.e. substructures within substructures are not included) inside radius which is determined for the Aquarius simulations Springel et al. (2008b, a).6 Our best fit is given by

(23)

Here, , denote our fit variables, is the luminosity from the smooth halo without substructures within and . We derive the normalization function in Eq. (LABEL:eq:Csub) from the simulations in reference Springel et al. (2008a) using a value of for the mass of the MW halo in the simulation and for the mass of the smallest resolved subhalos in the MW simulation. The smallest mass of subhalos in reality, , is determined from the free streaming length of DM at decoupling—an effect that erases structure on scales smaller than the free streaming length. In applying Eq. (LABEL:eq:Csub), we implicitly assume that the mass power-law scaling relation is valid down to the free streaming mass of DM halos. In the CDM universe, this is conventionally taken to be Hofmann et al. (2001); Green et al. (2005) (see Bringmann (2009) for a discussion of the range expected in various DM models). Note that potentially the power-law could flatten towards smaller mass scales although current simulations show no hints of such a behavior and since Einstein’s gravity is a scale-free theory, we do not expect such a behavior on theoretical grounds either. For DM halos more (less) massive than the MW we expect a larger (smaller) boost from substructures, simply because of the larger (smaller) mass range down to the minimum mass .

To motivate the scaling of the substructure luminosity boost with limiting substructure mass in Eqs. (23-LABEL:eq:Csub), we show how it derives from the substructure mass function, Springel et al. (2008b, a) and substructure luminosity scaling with the limiting mass of satellites, . The total luminosity of substructures scales as

Tidal truncation is responsible in shaping the substructure luminosity scaling parameter . Here we only sketch out qualitative arguments and leave details and comparison to numerical simulations to future work. First, tidal effects truncate the subhalo profile primarily in the outer regions. As a result of this, the subhalo acquires a steeper logarithmic slope than the canonical 3 from the NFW-profile. Second, tidal stripping imposes a mass-dependency and should be stronger for smaller substructure. The stripping efficiency (for an effective cross section ) only depends on the ambient DM density of the host halo, , which implies a scaling of the mean free path of a subhalo of . Since there is less mass at any given density for smaller substructures, modest density inhomogeneities within the host halo seen by the orbiting substructures cause a stronger mass loss of these smaller satellites. Third, the tidal truncation radius should be a function of host halo radius due to the increasing smooth density profile for smaller radii. If we pick an effective host halo radius which dominates the luminosity contribution of the substructures and were to fit the tidally truncated subhalo density profile with an NFW-profile, the resulting would be biased toward lower values implying a higher concentration compared to an isolated halo. In order to reproduce our value of , we would need a concentration mass relation of (employing the formula derived by reference Pinzke et al. (2009)). As expected, this is steeper than the concentration-mass relation found for isolated halos on dwarf galaxy scales, Ishiyama et al. (2011), as well as galaxy and cluster scales, Macciò et al. (2008).

We now derive the squared density profile for the substructures using

(26)
(27)
(28)
(29)

Here, represents the particle physics factor, where is defined by Eq. (13). The different density profiles have some impact on the luminosity from annihilating DM, although the details of the density profile can be neglected compared to the dominating boost from substructures (assuming DM to be cold). In Fig. 2 we compare the radial dependence of the accumulative luminosity from different smooth cluster density profiles to the boosted luminosity due to substructure for different mass scales. We recalculate the overdensity, , for the shallower density profile with and rescale the concentration parameter in Eq. (17) with Newman et al. (2011) to account for the more centralized scale radius in cluster with a shallow inner slope. The emission of this profile with is about 30% larger within compared to a cuspy NFW (). This difference is built up within (i.e. close to ). Hence the slope of the central part of a cluster has little influence for the DM luminosity within as long as the degeneracy in the lensing mass measurements with (which decreases for decreasing inner slope) has been taken into account. The emission from an Einasto density profile is about 50% larger than the cuspy NFW profile in the periphery of the cluster, where the difference is mainly built up at a few percent of . The increase in luminosity due to substructures is negligible in the center of halos, but integrated out to the boost relative to the smooth emission profile amounts to approximately 20 for dwarf galaxies, for galaxies, and for galaxy clusters. We stress that these boost factors are only realized in the region around the virial radius of each respective halo which is mostly tidally stripped for dwarfs in the MW. Hence a more realistic boost from substructures is probably much smaller in satellite dwarf galaxies. In addition, these boost factors are only realized for direct annihilation emission (continuum emission or final state radiation) or IC scattering off homogeneous seed photon fields (CMB). For IC scattering of SD photons, the overlap of final state leptons and SD photons is smaller which causes the substructure boost over the smooth emission to be reduced by roughly two orders of magnitudes. It should also be noted that the substructure boosted fluxes from clusters are much more extended than for dwarf galaxies, hence more difficult to detect with Cherenkov telescopes. In fact, it was shown in Sanchez-Conde et al. Sanchez-Conde et al. (2011) that inside , where the sensitivity of Cherenkov telescopes is maximized, the expected DM flux including substructures from the brightest dwarf galaxy is about an order magnitude higher than the brightest cluster. However, that work assumed a substructure boost of about 50 in the most massive cluster halos which is more than an order magnitude smaller than what we use in this work, hence in projection we expect this difference to be a factor few smaller.

Figure 2: Radial dependence of DM annihilation luminosity of smooth halo and substructures. The solid lines show the accumulative smooth luminosity from a cluster with the mass for three different density profiles; an Einasto profile with (light blue), a cuspy NFW profile with (thick dark blue), and a core NFW profile with (thin purple). The dashed lines show the accumulative luminosity from substructures for three different mass scales; an galaxy cluster (thick red), an galaxy (orange), and an dwarf galaxy (thin green). All luminosities have been normalized with the luminosity within from a cuspy NFW profile. We have assumed the standard value for the limiting substructure mass of . Note the large expected boost from substructures in clusters (), and the relative small boost in dwarf galaxies ().

In Fig. 3 we show the radial regions that dominate the DM annihilation luminosity, in particular we show the differential contribution to the DM luminosity per logarithmic interval in radius for three different mass scales. Solving for the maximum of this curve, in combination with Eq. (27), we find that the luminosity from the smooth NFW profile peaks at . Despite the cuspy nature of the density profile, the luminosity is not dominated by the central region but by the transition region where the profile steepens because of the larger volume available there. For large clusters with typical concentrations of , the luminosity from the smooth profile is focused to the regime around 10% of . In contrast, emission from substructures is mainly contributed by the outer parts of DM halos. As shown in Fig. 3, the product of annihilation emissivity and emission volume increases towards and only starts to drop outside this radius. Note that, even though most of the substructure mass density has been erased in the central regions of DM halos, a cluster in projection has a significant enhancement due to substructures at a radius of just a few percent of .

Figure 3: Plot that shows where most of the DM annihilation luminosity originates. We show the differential contribution to the DM annihilation luminosity per logarithmic interval in radius which corresponds to where is the density at . The solid lines represent cuspy NFW density profiles for three different mass scales; an galaxy cluster (thick red), an galaxy (orange), and an dwarf galaxy (thin green). The dashed lines show the contribution from substructures for the same three mass scales. All luminosities have been normalized with the luminosity within from a cuspy NFW profile. We have assumed the standard value for the limiting substructure mass of . For the smooth profile, the majority of the flux is delivered by a region around as indicated by the maximum value of the curves. In contrast, for substructures the emission is dominated by regions around .

ii.3 Inverse Compton emission

In this section we outline the basics of inverse Compton (IC) emission. As target radiation fields we consider cosmic microwave background (CMB) photons, and the light from stars and dust (SD). We derive an analytic model from which we can estimate the spectral and spatial distributions of SD as a simple function of cluster mass.

The standard IC source function is given by Rybicki and Lightman (1979):

where is the energy of the upscattering electrons and is the energy of the background photon field. We represent the Thomson cross section with and the full differential Klein-Nishina (KN) cross section is captured by Blumenthal and Gould (1970):

(31)

where

(32)

The full KN cross section accounts for the less efficient energy transfer between the photon and electron once the energy of the Lorentz-boosted photon in the electron rest frame comes close to such that the scattering electron experiences a significant recoil. This results in a steepening of the IC gamma-ray spectrum. In the low energy Thomson regime the IC spectrum , however when the IC spectrum steepens due to the KN suppression to . Here we denote the steady state electron spectrum by .

We account for two major contributions to the number density of radiative background fields ; the CMB photons and the infra-red (IR) to ultra-violet (UV) light emitted by SD. The number density for the SD is given by where the specific SD energy density is given in the Appendix, Eq. (50). We model the CMB photon spectrum as a photon gas that is isotropically distributed and follows a black body spectrum with the temperature  K:

(33)

Note that the typical energy of a black body photon before scattering is given by , where and are the number- and energy-density derived by integrating and over the photon energy , respectively.

The electrons injected from annihilating DM also suffer from diffusive and radiative losses. Hence we have to calculate the equilibrium spectrum of the electrons plus positrons denoted by in Eq. (LABEL:eq:ICemiss). We derive this stationary solution using the cosmic ray transport equation (neglecting convection and re-acceleration effects):

(34)

where denotes the diffusion coefficient and the energy loss term. The source function yields the number of electrons and positrons produced per unit time, energy and volume element at the position :

(35)

where the annihilation rate density is defined in Eq. (14). The sum runs over the kinematically allowed annihilation final states , each with a branching ratio and a differential spectrum that represents the number of electrons plus positrons resulting from an annihilation event. We use the differential spectra derived from high-statistics simulations in Ciafaloni et al. (2011); Cirelli et al. (2011) to compute the cumulative number of electrons and positrons resulting from neutralinos annihilate indirectly into pairs as well as pairs in the LP model. We use DarkSUSYto compute the spectra from our four BM models where only a fraction of the annihilating neutralinos is converted into electrons and positrons (see Sec. II.1 and Fig. 1 for further details).

The electrons and positrons loose their energy on a timescale that is shorter than the diffusive timescale in the ICM of galaxy clusters for cosmic ray electrons which is larger than the Hubble time in our energy range Berezinsky et al. (1997); Enßlin et al. (2011). Hence, we neglect the first term of the r.h.s. in Eq. (34), and derive an expression for the equilibrium number density:

(36)
(37)
(38)

Here we include the three main radiative loss processes for the cosmic ray electrons and positrons: (1) IC losses on CMB photons with the equivalent field strength of the CMB of , where is the cosmological redshift. (2) Synchrotron losses on ambient magnetic fields where we parametrize the magnetic field in the galaxy cluster by . We adopt a magnetic decline of in this work which follows from flux frozen magnetic fields. (3) IC losses on starlight and dust with an energy density given by Eq. (A.2), where we outline the derivation in the following.

The emission of galaxy clusters at infra-red (IR) and ultra-violet (UV) wavelengths emerges from dust and starlight in both the galaxies and the ICM (e.g. Porter et al. (2006) and Gilmore et al. (2009)). Three distinctive components dominate these wavelengths: a central galaxy, the intra-cluster light (ICL), and individual cluster galaxies. We decided to use the accurately measured spectral shape of dust and starlight in the interstellar medium of our Galaxy to model the emission from clusters. We then normalize the two spectral components – far IR dust and starlight at wavelengths ranging from the near IR to UV – individually by using stacked cluster data and employ a measured mass-to-starlight luminosity scaling relation derived from observations of the brightest cluster galaxy (BCG) Haarsma et al. (2010).

In Fig. 4 we characterize the spectral shape through a fit to the galactic spectra presented in Porter et al. (2006). The figure is showing the spectra at , which is the radius there the SD energy density of a galaxy cluster equals the energy density of the CMB black-body distribution. Inside this central radius the SD component is dominating, which is shown in Fig. 5 where we compare the energy densities from different radiation fields in a galaxy cluster with the mass . For this figure we use two different profiles for the SD energy density, where the total profile includes the contribution from the ICL, the BCG and all the galaxies, while the galaxies are excluded in the smooth profile. To compute the IC emission from SD, we require a non-negligible overlap of the relativistic lepton distribution resulting from DM annihilation and SD. In fact, one can show with a simple order of magnitude calculation that the overlap, , of the photon field of individual galaxies (starlight and dust emission) and the smooth DM density is very small so that we can neglect the starlight contribution from galaxies to the IC emission for the reminder of this work,

(39)

Here we assume that the exponential scale height of the stellar light is less than the scale radius of the galaxy halo, and . We denote the number of galaxies within with and the concentration of a galaxy DM halo with .

We find that even for a cluster with a relative small mass the energy density of the SD components dominates over the CMB and the magnetic fields (with a central B field of 3 ) within about 10% of . Outside this radius the CMB is dominating the energy density of the cluster. Note that while the magnetic field is always subdominant in the cluster for our assumptions, we keep its contribution to the total energy density for consistency. Also note that we extract the spatial distribution of the SD light in clusters from a stacked emission analysis of Sloan Digital Sky Survey (SDSS) data at the redshift Zibetti et al. (2005) and do not attempt to correct for a potential evolution in this component (c.f. Fig. 22).

The spatial distribution of far-IR to UV light emitted by SD is quite different from what is expected for the IC upscattered SD photons. We illustrate this in Fig. 6, where the accumulative luminosity from SD is compared to the IC SD photons for a cluster where the boost from substructures is excluded. The gamma-ray luminosity from IC upscattered SD photons is dominated from the inner parts where there is the largest overlap of the cuspy DM profile and the peaked SD distribution. This is in marked contrast to the accumulative SD luminosity in the optical that is dominated by the outer parts of the cluster. Interestingly, it rises as a function of radius and follows the distribution of an NFW mass profile outside . In the center, the electrons and positrons mainly cool by IC upscattering SD photons, hence the spatial dependence of the cooling cancels the distribution of the SD source function, which results in a distribution that approaches a density square profile in the center. We also find that the luminosity within from the total SD model is a factor three larger compared to the more realistic smooth SD model.

Figure 4: Spectral dependence of radiation fields in a cluster of galaxies. The black dotted line of the left peak show the spectrum of CMB photons using a black body with a temperature of  K. The crosses of the middle and right peaks represent the measured spectra from stars extending from the near IR to UV and dust at far IR wavelengths (SD), respectively, and are derived in Porter et al. (2006) for a galaxy. We normalize the individual SD spectrum separately using the observed luminosity from SD in clusters. The SD luminosity is related to the cluster mass through Eqs. (52-56), where we use have used the mass in this figure. We renormalize the SD spectra to the radius , where the smooth energy density of the SD light (see Fig. 5) equals the energy density of the CMB. The red dashed lines show the fitted SD spectral model.
Figure 5: Spatial dependence of the energy density of radiation fields in a cluster of galaxies. The energy density of the CMB (blue dotted line) is isotropic with . The energy density of the light from stars and dust (SD) is denoted by the red dashed line and the solid orange line for the total SD light and the smooth SD light, respectively. For comparison we show the energy density of the stars and a low dust model (SLD) with the solid yellow line. The SD light has been renormalized to a cluster with the mass . Finally we show the energy density of two magnetic field models with a central magnetic field of 3 . The magnetic field scales with the gas density to the power ; dark green dash-dotted line () and light green dash-dotted line (). Note that the SD radiation is dominating the energy density inside .
Figure 6: Comparing the spatial dependence of emission from stars and dust (SD) to IC upscattered SD light. We show the luminosity without substructures inside radius . The solid lines show the optical/IR emission from SD and the IC upscattered SD light is shown by the dashed lines. The thick red lines show the total emission including the brightest cluster galaxy, the intra-cluster light and the additional intra-cluster galaxies, while these galaxies have been cut out in the smooth component shown by the thin orange lines. We normalize the SD light and the IC upscattered SD with the smooth luminosity within for each component. We find that the gamma-ray IC luminosities are dominated by the central regime, while the SD light is mainly build up in the outer parts of the cluster. For comparison, we show that the SD light traces the NFW mass profile (dotted blue) outside .

ii.4 Cosmic-ray induced gamma-ray emission

In supernovae remnants and on scales of galaxies, there are convincing evidences of non-thermal populations. Especially, in the MW, the cosmic rays are observed directly as well as indirectly through radio, X-ray, and gamma-ray emission. On larger scales of the order of few 100 kpc up to Mpc, there are currently a vast number of observations of radio emission coming from radio mini halos in the centers of cooling flow clusters, radio relics in the periphery of clusters Kempner et al. (2004), as well as giant radio halos amounting to a total number of more than 50 diffuse radio sources in clusters Feretti (2003); Ferrari et al. (2008). This type of emission is expected in clusters since the formation process of a galaxy cluster is a very energetic processes that induces both turbulence as well as frequently occurring merging and accretion shocks which both are thought to accelerate relativistic non-thermal protons and radio emitting electrons to high energies. The precise origin of these electrons, in especially relics and giant radio halos, is still not settled. One possible scenario for the production of these electrons are hadronic CR interactions with ambient gas protons which results in charged and neutral pions that decay into electrons, neutrinos, and gamma-rays (see Dennison (1980); Vestrand (1982); Blasi and Colafrancesco (1999); Dolag and Enßlin (2000); Miniati et al. (2001); Miniati (2003); Pfrommer and Enßlin (2004a, b); Pfrommer et al. (2008); Pfrommer (2008); Kushnir et al. (2009); Donnert et al. (2010b, a); Keshet and Loeb (2010); Enßlin et al. (2011) ).7 Supporting evidence comes from the smoothness of the extended radio emission that often resembles that observed in thermal X-rays. This can be easily explained in the hadronic model since the long cooling time of cosmic ray protons (CRs) of order the Hubble time allows for a cluster-filling population of CRs to build up over the formation history Völk et al. (1996); Ensslin et al. (1997); Berezinsky et al. (1997). The production of these secondaries depend both on the gas and CR densities in the cluster, where the CR density roughly traces the gas outside the core regime and is slightly enhanced in the center. This density scaling implies that clusters are great targets for Cherenkov telescopes with a high sensitivity for the central parts of nearby clusters. Detecting the cluster gamma-ray emission is crucial in this respect as it potentially provides the unique and unambiguous evidence of CR populations in clusters through observing the bump at about 100 MeV in the spectra.

We adopt the universal spectral and spatial gamma-ray model developed by Pinzke & Pfrommer Pinzke and Pfrommer (2010) to estimate the emission from decaying :s that dominates over the IC emission from primary and secondary electrons above 100 MeV in clusters. The gamma-ray formalism was derived from high-resolution simulations of clusters of galaxies that included radiative hydrodynamics, star formation, supernova feedback, and followed CR physics using a novel formulation that trace the most important injection and loss processes self-consistently while accounting for the CR pressure in the equation of motion Jubelgas et al. (2008); Enßlin et al. (2007); Pfrommer et al. (2006). We highlight two main uncertainties of the models, namely the acceleration efficiency at formation shocks and CR transport parameters. First, we note that the overall normalization of the CR and gamma-ray distribution scales with the maximum acceleration efficiency at structure formation shock waves. Following recent observations at supernova remnants Helder et al. (2009) as well as theoretical studies Kang and Jones (2005), we assume the maximally allowed acceleration efficiency for strong shocks that transfers 50% of the shock-dissipated energy (kinetic energy corrected for adiabatic compressional heating at the shock) to CRs. These efficiencies drop quickly for weaker shocks Enßlin et al. (2007) which dominate the gravitational energy dissipation inside galaxy clusters Pfrommer et al. (2006). These efficiencies will have to be corrected down if this value is not realized at strong structure formation shocks. Second, these simulations (and by extension the analytic model) neglect active CR transport such as streaming and diffusion relative to the gas, i.e., we assume that advective transport dominates and CRs are tightly coupled to the gas via magnetic fields tangled on sufficiently small scales which produces centrally enhanced profiles. However, CR diffusion and streaming tends to drive the CR radial profiles towards being flat, with equal CR number density everywhere. While the CR streaming velocity is usually larger than typical advection velocities and becomes comparable or lower than this only for periods with trans- and supersonic cluster turbulence during a cluster merger. As a consequence a bimodality of the CR spatial distribution may result with merging (relaxed) clusters showing a centrally concentrated (flat) CR energy density profile Enßlin et al. (2011). This translates into a bimodality of the expected diffuse radio and gamma-ray emission of clusters, since more centrally concentrated CR will find higher target densities for hadronic CR proton interactions Enßlin et al. (2011). As a result of this, relaxed clusters could have a reduced gamma-ray luminosity by up to a factor of five.

Iii Gamma-ray spectra

Spectrally resolved indirect DM searches have the advantage of probing different DM models through their characteristic spectral distributions. To make current and future DM searches more effective it is important to know in which energy band to focus the efforts in order to maximize potential DM signals over the expected background.

We focus in this section on the spectral distribution of gamma-rays from clusters. In the LP model, DM annihilation radiation includes final state radiation and gives rise to substantial amounts of electrons and positrons that IC upscatter background radiation fields to high energies. We also consider four supersymmetric DM models with a high gamma-ray yield in the form of continuum emission and IC induced emission. In addition to the annihilating DM, we estimate the gamma-ray flux induced by shock-accelerated CRs. Note that all fluxes in this section are derived within an angle corresponding to and are convolved with a PSF of .

We compare the calculated fluxes to gamma-ray upper limits (95% c.l.) set by Fermi-LAT after 18 months of observations Ackermann et al. (2010). In particular, to achieve a more reliable comparison we adopt, if nothing else is stated, the maximal spatially extended limits since the gamma-ray flux from our brightest clusters all have an angular extent on the sky when the boost from substructures is included. In fact, published Fermi upper limits are not derived for the kind of extended emission that we find with our treatment of substructures, hence these Fermi limits may not be adequate for some of these clusters. This and other recent work highlights that the improved substructure models imply that Fermi limits will have to be re-calculated to accommodate for the large source extensions. Note, however, that for most clusters the assumption of a point source is well justified. The flux upper limits, that we compare to, are a function of spectral index . However, in the relevant energy range,  GeV, the spectral index varies within . This changes the photon flux upper limits by , with Fermi-LAT being more sensitive to a hard spectrum Ackermann et al. (2010).

In Fig. 7 we show the differential flux from the Fornax cluster which is one of the best clusters for indirect DM searches due to its high DM annihilation fluxes and low CR induced fluxes. We show the emission of four different supersymmetric BM models and contrast it to the emission induced by CRs. Comparing this emission to the differential gamma-ray upper limits set by Fermi-LAT, we find that the upper limits are not violated. The predicted DM flux that is dominated by the continuum emission from and models (shown in the left panels) are about a factor few below the upper limits, making it hard for Fermi to probe these kind of DM models in the near future without a significant improvement in the analysis from e.g. stacking of clusters and improved methods for analyzing extended sources. Furthermore, the gamma-ray signal induced by CRs is expected to be about a factor 10 below the DM continuum flux from the and BM models at 10 GeV. However, the IC emission from upscattered CMB and SD photons is at least a factor of a few lower for and a factor 1000 lower for , than the expected flux from CRs above 100 MeV, making it very hard to distinguish such a signal from the foreground due to the similar spectral index. For energies below 100 MeV we expect IC emission by primary shock accelerated electrons to be dominating Pinzke and Pfrommer (2010) over the supersymmetric DM induced leptons. Hence in clusters, the IC emission from supersymmetric DM BM models can be neglected compared to the CR pion and DM continuum emission.

Figure 7: Comparing the differential flux from different models: we show the continuum emission from DM benchmark (BM) models (blue dashed), electrons and positrons from DM BM models that Compton-upscatter CMB photons (green dotted) as well as dust and star photons (orange dash-dotted), and CR induced gamma-ray emission (red solid). Each panel is associated with an individual DM BM model; upper left , lower left , upper right , and lower right . The emission is calculated for the Fornax cluster using a point spread function of . The substructures boost the gamma-ray flux from IC upscattered CMB and continuum emission by a factor of 890 while the IC upscattered SD photons are only boosted by 20. Extended upper limits () from Fermi-LAT after 18 months Ackermann et al. (2010) are also shown. In the near future we find it hard to detect even the brightest BM models, and , where the continuum emission is dominating the total emission in the GeV energy range.

For the LP DM models, however, the main contributions to the expected extended gamma-ray flux is coming from the IC upscattered photons on various photon background fields. Depending on the spatial and spectral distribution of electrons and positrons as well as the photon background fields, the resulting spectral distributions of gamma-rays can differ greatly. Hence it is interesting to understand which spectral features and energy regimes are dominating.

If the enhancement due to substructures in clusters is significant, then the distribution of electrons and positrons follows the radial profile of substructures outside the cluster center where the smooth DM density profile is dominating. In addition, the cooling of the steady state electron and positron distribution is dominated by SD photons in the center. Especially since the SD distribution is centrally peaked and the substructure distribution peaks in the outskirts around , the overlap between electrons/positrons and SD photons is small. This results in a suppression of the IC upscattered SD photons relative to the IC upscattered CMB photons, the final state radiation, and the continuum emission. However, if substructures are only marginally dominating over the smooth distribution in a cluster, the SD component becomes relatively more important.

In Fig. 8 we show the total IC emission as well as its individual contributions from different IC upscattered radiation fields. The left panel shows the gamma-ray emission from the LP model and the right panel from the BM model, and for comparison we overplot the emission expected from the CRs. Due to the flat electron and positron spectrum resulting from the LP DM model and the smaller mean energy of the CMB compared to the SD, we find that the upscattered CMB photons dominate the total DM IC emission in the energy regime below 100 GeV, while the SD dominate above this energy. For the BM models, this transition energy is shifted towards smaller energies since the electron and positron spectrum has a steeper spectral distribution (see Fig. 1). At the highest gamma-ray energies of about 100 GeV and above, the IC from starlight steepens because it probes the high energy tail of the electrons and more importantly, it suffers from the Klein-Nishina suppression. When substructures are present, most of their flux resides in the outer parts of clusters. However, if we remove the boost from substructures, the density profile of electrons and positrons is more centrally peaked, and the relative importance of the IC upscattered SD increases by a factor (see Fig. 8)8. A larger fraction of electrons (in the core) will now also cool by Compton-upscattering SD photons which suppresses the IC-upscattered CMB light. For comparison, we include the contribution of IC upscattered SD photons and to bracket the uncertainty in the SD model, with a smaller amount of dust. In this low-dust model, we have reduced the energy in dust by a factor 10. However, the resulting flux from the IC upscattered dust in this model only decreases with a factor that is slightly smaller than 10 since the IC cooling of the electrons and positrons also decreases.

Because of the large total boost factor for the LP model (), we overproduce the upper gamma-ray limit in the  GeV energy interval set by Fermi-LAT by about a factor 100. This strongly constrains the boost from both substructures and the SFE. Additionally, these constraints might improve with future more sensitive Cherenkov telescopes such as CTA. However, considering that an indirect detection of DM in clusters relies on the boost from substructures whose main contribution comes from the periphery, we conclude that the wide angular extent of clusters on the sky in gamma-rays suggest that these sources are not ideal for Cherenkov telescopes since their sensitivity drops linearly with source extension.

Figure 8: Comparing the flux from different inverse Compton upscattered radiation fields. We show the differential inverse Compton emission induced by leptophilic DM in the left panel and by the benchmark model in the right panel. The contribution from each individual radiation field from top line to bottom line; CMB (green dashed-triple-dotted), dust (purple dashed), low dust model (grey dashed), and stars (blue dashed-dotted). The sum of the three components are shown with the orange solid line. The red solid lines show the CR induced gamma-ray flux. The black arrow shows the spatially extended differential upper limit from Fermi Ackermann et al. (2010) indicating that the LP model assumptions such as Sommerfeld and/or substructure boost are in conflict with the upper limit. All fluxes are calculated for the Fornax cluster within using a point spread function of . For this cluster, the enhancement due to substructures from IC upscattered CMB and SD photons is 890. The saturated Sommerfeld boost is 530.

It is also interesting to compare the total contribution from the LP model, the brightest BM model (), and the CR induced emission. In Fig. 9 we show the integrated flux from Fornax for our different gamma-ray models and compare it to the unresolved integrated flux upper limit on Fornax set by Fermi-LAT where they averaged the flux over the energy range assuming a spectral index of 2. Again, the annihilation flux in the LP model is in conflict with the upper limits by Fermi-LAT, although only by a factor of 30 which is a factor few less constraining than the differential flux in the energy range  GeV. As shown in this figure, the LP model is dominating the entire gamma-ray energy range up to the DM rest mass energy in this model of about 1 TeV. Finally, we note that the flux the DM BM model is larger than the predicted emission from the CRs in the  GeV energy regime. Hence, for an experiment with a high sensitivity even for extended sources, the prospects for detecting the BM model over the expected gamma-ray background induced by CRs looks promising in clusters. Present day Cherenkov telescopes, however, have a trigger region that is smaller than the size of clusters, hence it has to be increased to several degrees to overcome problems with background estimation. In addition, even though the projected CTA point source sensitivity (, 50h) shows the potential of this experiment in constraining leptophilic models as well as BM models with a very large neutralino mass  TeV, we find that analysis techniques have to be developed that enable the detection of extended sources without too much degradation of sensitivity.

Figure 9: Comparing the energy integrated flux from different models. We show the emission from CR induced emission (red solid), a leptophilic (LP) model that includes both final state radiation and IC upscattered CMB, dust and starlight (green dash-dotted), and the benchmark model (BM) that includes continuum emission, and IC upscattered CMB, dust and starlight (blue dashed). The black arrow shows the spatially extended integrated flux upper limit set by Fermi-LAT again indicating challenges for the assumptions underlying the LP model. We also show projected CTA point source sensitivities (, 50h). The emission is calculated for the Fornax cluster using a point spread function of . The boost from Sommerfeld and substructures is about 530 and 890, respectively.

We continue by comparing the estimated differential flux from Fornax to three other clusters in Fig. 10; the close-by and well studied Virgo cluster, the X-ray bright Perseus clusters, and the massive merging Coma cluster. We summarize in Table 2 the expected gamma-ray flux in our DM and CR models and the corresponding boost factors. We find high DM induced gamma-ray fluxes from the Fornax and Virgo clusters, which confirms them as promising targets for indirect DM searches. Especially in the  GeV energy range it is quite striking how constraining the upper limits are for some of the clusters. At these energies, Fermi-LAT has its peak sensitivity due to a combination of increasing effective area and decreasing source spectra as a function of energy. The upper limits for Virgo and Perseus are unfortunately background contaminated by AGN activity from M87 and NGC 1275, respectively, and do not gain much from the increased sensitivity. In addition, the angular radius of Virgo is about , while the extended upper limit from Fermi-LAT is calculated assuming a radius. In fact, the ratio of the virial radius to the assumed extension for the upper limits for Fornax, Coma, Virgo are , respectively. Similarly, it is interesting to compare the radius that contains 68% of the flux to the assumed extension, , which is given by for Fornax, Coma, Virgo, respectively. This further motivates a re-calculation of the Fermi upper limits to account for the full extension of the galaxy clusters. We also find that the theoretical expectation of a high CR-induced flux in Coma and Perseus may significantly complicate indirect DM searches in those clusters if the boost factors are much lower than assumed in this work. The Fornax cluster, however, is a great target for indirect DM studies because the relative low gamma-ray flux from CRs, absents of an active AGN, and high DM gamma-ray flux.

What is the figure of merit for selecting the most promising cluster targets for indirect DM searches? To this end, we employ the luminosity-to-mass scaling relations. The gamma-ray luminosity from the smooth density profile is given by Pinzke et al. (2009)

(40)

and the total DM luminosity that includes boost factors for the LP and the BM models is given by

(41)
(42)

We note that the IC from upscattered SD photons scales slightly softer with mass, hence gives rise to a larger fraction of gamma-rays in low mass clusters compared to the total gamma-ray flux. Also note that we have not included the SFE for the LP model in Eq. (41) since it saturates to a constant value in the substructures. However, if only a small fraction of the DM resides in the subhalos, i.e. the scaling of the LP DM becomes shallower due to the mass dependence inferred from the velocity dispersion:

(43)
(44)
Cluster : :
DM-LP DM-BM- CR- DM-LP DM-BM- CR- [Mpc]
Coma
Perseus
Virgo
Fornax