Constraining Extended Gammaray Emission from Galaxy Clusters
Abstract
Cold dark matter models predict the existence of a large number of substructures within dark matter halos. If the cold dark matter consists of weakly interacting massive particles, their annihilation within these substructures could lead to diffuse GeV emission that would dominate over the annihilation signal of the host halo. In this work we search for GeV emission from three nearby galaxy clusters: Coma, Virgo and Fornax. We first remove known extragalactic and galactic diffuse gammaray backgrounds and point sources from the Fermi 2year catalog and find a significant residual diffuse emission in all three clusters. We then investigate whether this emission is due to (i) unresolved point sources; (ii) dark matter annihilation; or (iii) cosmic rays (CR). Using 45 months of FermiLAT data we detect several new point sources (not present in the Fermi 2year point source catalogue) which contaminate the signal previously analyzed by Han et al. Including these and accounting for the effects of undetected point sources, we find no significant detection of extended emission from the three clusters studied. Instead, we determine upper limits on emission due to dark matter annihilation and cosmic rays. For Fornax and Virgo the limits on CR emission are consistent with theoretical models, but for Coma the upper limit is a factor of 2 below the theoretical expectation. Allowing for systematic uncertainties associated with the treatment of CR, the upper limits on the cross section for dark matter annihilation from our clusters are more stringent than those from analyses of dwarf galaxies in the Milky Way. Adopting a boost factor of from subhalos on cluster luminosity as suggested by recent theoretical models, we rule out the thermal cross section for supersymmetric dark matter particles for masses as large as 100 GeV (depending on the annihilation channel).
keywords:
dark matter experiments, gamma ray experiments, galaxy clusters1 Introduction
The existence of dark matter (DM) in the universe has so far only been deduced from its gravitational effect, due to the lack of electromagnetic interactions of the DM with itself or with baryonic matter. There are several elementary particle candidates for DM in various extensions of the standard model of particle physics (Bertone et al., 2004). Weakly interacting massive particles or WIMPs (with a selfinteraction crosssection at roughly the weak scale) are one class of the popular dark matter candidates. These particles could be related to the electroweak symmetry breaking which is currently being explored by experiments at the LHC. For example, within the framework of the minimal supersymmetric standard model (MSSM), the lightest neutralino emerges as a candidate WIMP that is stable over cosmological timescales and can annihilate into standard model particles. WIMPs behave as cold dark matter since their primordial velocity dispersion is negligible. High resolution Nbody simulations show that cold dark matter halos contain a population of selfbound substructures (subhalos) whose number decreases with increasing subhalo mass as with (Diemand et al., 2007; Springel et al., 2008; Gao et al., 2012)
Much effort has been devoted to the search
for WIMPs either directly or indirectly. Direct detection involves
identifying the rare events of DM scattering off ordinary matter or
searching for new particles near the weak scale at the LHC. Indirect
detection involves looking for the annihilation or decay products of
dark matter in cosmic rays and gamma rays. In particular, pair
annihilation produces gammaray photons at a rate proportional to the
square of the dark matter density, which then propagate, almost
without absorption, to the observer. In this case, the Galactic centre
should be the brightest gammaray source on the sky (Springel et al., 2008, and
references therein). Extended emission (distinct from the
central point source) was reported from the central around
the Galactic centre by Hooper & Goodenough (2011); Hooper & Linden (2011b).
It has recently been reported that the ray emission from the region around the Galactic centre exhibits a linelike excess at energies GeV (Bringmann et al., 2012; Weniger, 2012; Tempel et al., 2012; Su & Finkbeiner, 2012). The intepretation of this signal as arising from dark matter particles, however, is controversial (see Boyarsky et al., 2012).
Targeting the entire sky rather than the Galactic centre in searching for annihilation radiation may seem a good strategy since this takes advantage of the largescale distribution of dark matter while avoiding some of the uncertainties arising from the astrophysical modelling of galactic gammaray sources. However, the fact that we are located near the centre of the Galactic halo and that most of the annihilation emission outside the Galactic centre is produced by dark matter substructures (Diemand et al., 2007; Springel et al., 2008) results in a gammaray map from annihilation that is almost uniform on large scales. This makes detection within the Milky Way halo a difficult task, exacerbated by the additional uncertainty of having to model the extragalactic background, which is more important on large scales (Zaharijas et al., 2010; Baxter & Dodelson, 2011).
Dwarf galaxies are the most DMdominated objects known, are relatively free from astrophysical contamination and appear compact on the sky. They are therefore promising targets to search for DM annihilation radiation. Recent joint analyses of eight to ten dwarf galaxies (GeringerSameth & Koushiappas, 2011; Ackermann et al., 2011) resulted in no significant detection but have began to rule out the canonical annihilation crosssection of for DM masses below GeV.
Galaxy clusters are the most massive virialized DM structures in the universe and are also good targets for indirect DM searches. The presence of a large population of DM substructures (or subhalos) predicted by numerical simulations further enhances the detectability of DM in clusters. Although the total mass within subhalos amounts to only 10 to 20 percent of the total halo mass, the density enhancement within subhalos can boost the total cluster annihilation luminosity by a factor as high as 1000 when extrapolated down to a subhalo mass limit of one Earth mass, the fiducial cutoff in the primoridal power spectrum of density fluctuations for a typical 100 GeV WIMP (Gao et al., 2012; Pinzke et al., 2011). As the distribution of subhalos is much less concentrated than that of the smooth main halo, the total annihilation emission from clusters is predicted to be extended. Thus, attempts to detect DM annihilation assuming a point source or NFWsquared profile could miss most of the signal. In fact, just such a search using the 11month FermiLAT data has yielded no significnat detection of emission from six clusters (Ackermann et al., 2010).
Using the 45month data, we consider possible contributions from cosmic ray (CR) induced gammaray emission and from DM annihilation. For the former (which can be as high as, or higher than the emission from cluster DM annihilation (Jeltema et al., 2009; Pinzke & Pfrommer, 2010; Pinzke et al., 2011), we adopt the semianalytic method developed by Pinzke & Pfrommer (2010). For the later, we adopt the model proposed by Gao et al. (2012) for the cluster DM annihilation profile. We provide constraints on both the CR and DM components for the three galaxy clusters analyzed by Han et al. (2011): Coma, which is predicted to have the highest signaltonoise according to Gao et al. (2012), and Fornax and Virgo which are predicted to have the lowest astrophysical contamination according to Pinzke et al. (2011).
The current paper replaces an earlier version by a subset of the
authors (Han et al., 2012, arXiv:1201.1003). After submission of that version, it was
pointed out to us that a number of point sources are present in the
full threeyear LAT data which were not detected significantly in the
data used for the ”official” Fermi point source catalogue available at
the time of our analysis, the LAT 2year point source catalogue (2FGL;
Nolan et al. (2012)). We now carry out our own point source detection in the
regions of interest and find several new point sources.
Huang et al. (2012) have recently reported a failure to detect significant DM annihilation emission from a combined analysis of eight galaxy clusters. Our work differs from theirs in several respects: firstly, we assume a DM annihilation profile based on high resolution cosmological simulations (Gao et al., 2012); secondly, we assess the impact of cosmic rays in the detection of dark matter; and finally, we include in our sample the Virgo cluster which turns out to be the best candidate. The constraints we set on the annihilation crosssection are consistent with those of Huang et al. (2012).
The paper is organized as follows. In section 1 we describe the data and provide an overview of the models of the Virgo, Fornax and Coma galaxy clusters regions used in the analysis (see Table 1). The specification of the nonstandard components of the models (dark matter and cosmic rays brightness profiles) is provided in Sec. 2. The constrains on CR emission and DM annihilation that we obtain are summarized in section 3 and discussed in Sec. 4.
The cosmological parameters used in this work are the same as those assumed by Gao et al. (2012): , , .
1.1 Data preparation
We analyze the first 45 months of data (04/08/2008 to 20/05/2012) from
the FermiLAT,
Minimum Energy  100 MeV 
Maximum Energy  100 GeV 
Maximum zenith angle 
100 degrees 
Event Class 
3 (P7CLEAN) 
DATAQUAL 
1 
LAT CONFIG 
1 
ABS (ROCK ANGLE) 
degrees 
ROIbased zenith angle cut  yes 
We list the basic properties of the three clusters in Table 1.
Coma  Fornax  Virgo (M87)  
RA (deg)  194.9468  54.6686  187.6958 
DEC (deg)  27.9388  35.3103  12.3369 
(Mpc) 
95.8  17.5  16.8 
() 
1.3e15  2.4e14  7.5e14 
(deg) 
1.3  4.1  6.2 

5.9e5  4.1e4  1.2e3 
Enhancement due to subhalos within 
1.3e3  6.5e2  1.0e3 
1.2 Maximumlikelihood fitting
We use the pyLikelihood tool shipped with the Fermi Science
Tools software package (version v9r27p1fssc20120410) to perform a
maximum likelihood (ML) analysis (Mattox et al., 1996). After applying
appropriate data cuts, as described in section 1.1, we bin
the data into 0.1 degreewide pixels
and 30 logarithmic energy bins within a radius of 10 degrees around
each cluster. This large radius is chosen to account for the large LAT
PSF size at low energies ( degrees at 100 MeV
In the standard Fermi likelihood analysis, the photon counts within each pixel are treated assuming Poisson statistics for each energy bin to calculate the likelihood. The bestfit parameters are obtained when the likelihood for the entire data set is maximized. The significance of a given component of interest (e.g. DM or CR) from the ML fitting is quantified by the likelihood ratio statistic,
(1) 
where is the maximum likelihood for the full model and is the maximum likelihood for the null hypothesis, i.e, the model without the component of interest. According to Wilk’s theorem, this test statistic, , approximately follows a distribution when the null hypothesis is true, with one degree of freedom for our case where the normalization is the only extra parameter in the alternative model. The probability that a given value of arises purely from fluctuations of the null hypothesis is:
(2) 
The factor comes from the constraint that the normalization parameter be nonnegative. The significance of a detection can thus be quoted as (onesided Gaussian confidence). Upper limits on the extra normalization parameter are obtained by searching for a null hypothesis where in the full model is constrained to be equal to the upper limit, , so that , corresponding to the 95% confidence interval.
1.3 Model
For the analysis we constructed a model to fit the data including all known foreground and background emission, as well as DM and CR components, as appropriate. We include all the point sources from 2FGL within a radius of 15 degrees from the cluster centre in the model, plus the most recent galactic (GAL) and extragalactic (EG) diffuse emission given by the template files gal_2yearp7v6_v0.fits and iso_p7v6clean.txt. Additionally, we have searched the 45month data for new point sources; we detect several of them within the cluster region (see Appendix A for more detail) and these are also included in our model. The normalization of the GAL and EG diffuse components are allowed to vary during the fitting. Within the cluster virial radius there are two 2FGL point sources and one newly detected point source in Fornax, six 2FGL, including the central AGN (M87; Abdo et al., 2009), plus four newly detected ones in Virgo. We allow the normalization and powerlaw spectral index of these thirteen point sources to vary freely. In addition, the parameters of all sources with variability index greater than 50 located within 10 degrees of the cluster centres are allowed to vary. Parameters for the other point sources are fixed as in the 2FGL catalog. From now on we refer to the model with GAL, EG and the known point sources as the “base model”.
A DM annihilation surface brightness template (given by the dimensionless factor J, see Eqn. 4 in Sec. 2.1) is generated for each cluster out to a 15 degree radius by summing up both the contribution from a smooth NFW profile and the contribution from subhalos. This Jmap is used to fit for extended cluster annihilation emission. For the point source model, the integrated factor (see Eqn. 5) is used to derive an annihilation crosssection from the fitted total flux. Similarly, a CR photon template is generated for each cluster out to three times the cluster virial radius, where the surface brightness has dropped to below of the central value and beyond which the model is not reliable. Images for various model components are shown in Fig. 1 taking Virgo as an example. We discuss these templates in more details in Sec. 2.
In the traditional Fermi analysis, the EG template is treated as a smooth component where all emission below the nominal point source detection limit is assumed to come from a smoothly distributed diffuse component. In this work, we also consider a more realistic one where a fraction is assumed to be contributed by fainter point sources with a numberflux relation which extrapolates smoothly from that measured for brighter sources. In this case the photon counts within a given pixel are no longer Poissondistributed since the photons arrive in packets. In principle, one can use the full distribution of photon counts from a population of randomly placed point sources to calculate the likelihoods and , but Eqns. 1 and 2, and the corresponding discussion, are not affected. However, since the full distribution of photon counts in this case (Han et.al., in prep.) is complicated and difficult to implement in the likelihood analysis, instead of recalculating and , in this work we use MonteCarlo simulations to reevaluate the distribution of for the more realistic background model and provide corrections to the results of the standard analysis where needed.
2 Modeling gammaray emission in clusters
We model the observed gammaray emission in clusters with several components as shown in Fig. 1: the galactic foreground (GAL), the extragalactic background (EG), emission from known point sources, DM annihilation and CRinduced emission. The GAL and EG diffuse emission are given by the most recent templates, gal_2yearp7v6_v0.fits and iso_p7v6clean.txt, which can be obtained from the FermiLAT data server, while the point sources include those from the LAT 2year point source catalogue, 2FGL (Nolan et al., 2012), as well as several, newly detected by us, in the 45month data. In addition, an improved EG model which includes a population of undetected sources is also analyzed. We now describe in detail our models for DM annihilation and CR emission.
2.1 Dark matter annihilation emission
The gammaray intensity along the lineofsight due to DM annihilation is given by:
(3) 
where is the DM particle mass; the density of DM; the particle model dependent term giving the differential number of photons produced from each annihilation event as a function of energy, , in a particular annihilation channel, ; and is the velocityaveraged crosssection (or annihilation rate) for that channel, which is predicted to be constant in the low velocity limit appropriate to presentday cold DM particles (see e.g., Jungman et al. (1996)). The lineofsight integration of the density squared is often expressed in terms of a dimensionless factor,
(4) 
If the source size is much smaller than the instrumental beam size, a point source approximation is applicable. In this case, the integration of over a large enough solid angle, , is used to determine the total flux for the point source, .
The cluster annihilation emission is modeled with the extended profile suggested by Gao et al. (2012). However, for part of the analysis and for comparison purposes, we will also use the point source approximation which, although inappropriate, has been employed in all previous analysis of FermiLAT data from clusters. We shall refer to models that assume these two profiles respectively as EXT and PT. If the cluster follows a smooth NFW profile, then its integrated factor which determines the total annihilation flux can be found as
(5) 
Here is the angular diameter distance to the cluster and and are the characteristic density and radius of the NFW profile. They are related to halo concentration, , and virial radius through the relations, and , with the critical density of the universe, the cluster virial radius within which the average density is and the concentration parameter, , is given by the following massconcentration relation:
(6) 
(Duffy et al., 2008). Here, is the mass enclosed within . Extrapolating to a cutoff mass of , the existence of subhalos will increase this flux by a factor
(7) 
Gao et al. (2012). Using the results of the simulations by these authors, the surface brightness profile of subhalo emission can be fitted within by the following formula:
(8) 
Below we fit the subhalo emission surface brightness beyond the virial radius and extrapolate to several times the virial radius using an exponential decay,
(9) 
The total annihilation profile is the sum of the contributions from a smooth NFW profile and the subhalo emission. This is completely dominated by subhalo emission except in the very centre of the cluster. We show the total annihilation profile and its decomposition into main halo and subhalo contributions in the left panel of Fig. 3, taking Virgo as an example. This profile is further inflated after convolution with the LAT point spread function.
We consider three representative annihilation channels, namely into ,
and final states. The annihilation spectrum
is calculated using the DarkSUSY package (Gondolo & Silk, 1999),
We note that the electroweak corrections recently proposed by Ciafaloni et al. (2011) (see also Cirelli et al. (2011)) can bring visible differences to the leptonic channel spectra at high WIMP masses before IC scattering. However, since IC photons dominate at the high mass end and the electroweak correction only significantly changes the positron yields at low energy, thus having little effect on the IC spectrum, the electroweak correction to the total spectrum is still negligible. The total photon yields are shown in Fig. 2. The almost flat spectrum with a cutoff around the energy corresponding to the WIMP mass comes from prompt annihilation emission including continuum secondary photons and final state radiation from charged final state particles. The low energy rise originates from IC scattered CMB photons.
2.2 Cosmicray induced gammaray emission within clusters
The cosmic ray induced gammaray emission is calculated following a semianalytic prescription, derived from high resolution numerical simulations of galaxy clusters, that models cosmic ray physics self consistently (Pinzke & Pfrommer, 2010). The gammaray photon production rate (or source function) from pion decay is found to be separable into a spatial and a spectral part:
(10) 
where the spatial part, , is proportional to the square of the gas density profile multiplied by a slowly varying radial function parametrized by cluster mass. The spectral part, , is almost independent of cluster mass and has a powerlaw form, , for the energy range GeV but flattens at low energies, as shown in Fig. 2. We summarize the detailed form of and plus the gas density profile for the three clusters derived from Xray observations in the Appendix.
The differential gammaray flux from this source function, , is simply the integral of along the lineofsight. This prescription is derived from the average emission profile for a sample of simulated clusters for a realistic choice of parameter values (e.g., for the maximum shock acceleration efficiency, ). In addition to the uncertainties in the model parameters there is also uncertainty in the observationally derived halo mass and gas density profile. In this work, we simply assume that the shape of is given by the model described above and account for the uncertainty in the model parameters, as well as sample variance with an additional normalization parameter, , so that,
(11) 
We take as our fiducial CR model and also consider the case when is fitted from the actual gammaray data as an optimal model. In the right panel of Fig. 3 we compare the CR profile for the fiducial model to the expected DM annihilation profile within our three clusters, assuming a fiducial DM particle model with particle mass, , annihilating through the channel with crosssection, . In general the CR emission is more centrally concentrated than the annihilation profile since the CR trace the gas profile. It can be readily seen that Fornax has a particularly low CR level while Coma is CR dominated. Coma has steeper profiles due to its larger distance and hence smaller angular size.
3 Results
3.1 Constraints on CR emission
With all the model components defined above, we first proceed with ML fitting for a model with no DM annihilation but with cosmic rays, the “CRonly” model hereafter. Note that the GAL and EG backgrounds, as well as the nearby point sources, are always included in the analysis, as described in section 1.2. The results for the CRonly model fits are listed in Table 2. The fitted CR levels all agree within a factor of three with the theoretical predictions. While Fornax is most consistent with no CR emission due to its intrinsically low CR level, the derived upper limit for Coma already rules out the fiducial value at 95% confidence.


(ph cm s) 




Coma  0.5  2.4e09  5.2  2.6  0.6  
Fornax  4.8  1.8e09  0.2  0.1  6.4  
Virgo  1.2  2.1e08  8.4  2.8  1.6 
3.2 Constraints on DM annihilation
Given the low significance of the CR detection in the CRonly model, it is not safe simply to adopt the best fit values for further extraction of the DM signal. Instead, we consider the following four families of cosmic ray models in the presence of a DM component:
 FiducialCR model.

The CR level is fixed to the theoretical expectation, . Since this value exceeds our derived upper limit for Coma, we exclude Coma from further discussion of this family.
 OptimalCR model.

The CR level is taken as the bestfit value listed in Table 2.
 FreeCR model.

The normalization of the CR level is left as a free parameter in the fit.
 NoCR model.

No CR emission is considered, only DM.
For each family, both point source (PT) and extended (EXT) profiles are considered for the DM component (the former merely for comparison with earlier work). Note that when calculating the for DM, the null hypothesis refers to the full model excluding only the DM component, or equivalently, to the base model plus a CR component modelled according to one of our four families of CR models. We show results for the , and DM annihilation channels.
For none of the combinations of DM and CR models considered here, do we obtain a detection of DM at high significance in any of the three clusters. The highest significance is obtained for Virgo for the channel in a DM model that has a particle mass of 30 GeV and the EXT profile, in the absence of CR. In this case, we find , corresponding to 3.4. This reduces to in the FreeCR model and to less than in the FiducialCR model.
The value of for the noCR model for Virgo can be compared with the value of reported in an earlier version of this paper (arXiv:1201.1003v.1) from a similar analysis of the 2year Virgo data (see Fig. 4). The decrease in significance is entirely due to the subtraction of the new point sources which we have detected in the Virgo region and which were not catalogued in the 2FGL. These previously undetected sources happen to lie within the virial radius of Virgo and can mimic the extended emission expected from DM annihilation. In fact, fits assuming an EXT profile but a powerlaw, rather than a DM annihilation spectrum, result in a similarly high significance detection, , and a bestfit spectral index . This is the typical spectral index of Fermi point sources (including the newly detected ones). The preference for a 30 GeV DM particle mass in the DM fits reflects a preference for a spectrum around 1 GeV, the energy scale from which most of the significance arises.
In fact, the significance of the Virgo detection is further reduced when we take account in the analysis of a possible undetected point source population, as we shall do in Appendix B. Thus, in what follows we use our analysis exclusively to set upper limits on the flux and annihilation crosssection.
The channel
In Fig. 5 we show the 95% confidence upper limits on
the DM annihilation flux and compare them to the CR levels. For each
cluster, the coloured stripes are defined by the minimum and maximum
upper limits corresponding to the four families of CR models. The optimal
CR levels in the three clusters are all comparable to the fitted DM
flux, and the DM flux upper limits for the four different CR models
vary only within a factor of two, with the NoCR and FiducialCR
The flux upper limits are translated into crosssection upper limits in Fig. 6, using Eqn. 3. These are also shown as coloured regions reflecting the variation in the different treatments of CR. Although the predicted flux upper limits decrease slowly with DM particle mass and remain within the same order of magnitude for the mass range considered, the resulting crosssection upper limits increase by a factor of 100 from low to high particle mass. This is because low mass particles correspond to higher DM number densities (the factor in Eqn. 3) for a given mass density, so to obtain the same flux level, the required crosssection must be smaller for low mass particles. With an enhancement of order due to subhalos, a much lower crosssection is needed (by a factor of at least 100) for extended annihilation models to achieve a slightly higher flux upper limit than point source models.
Our crosssection limits drop below the fiducial thermal crosssection
of for GeV. Of the three clusters, Virgo has the highest flux upper
limits but it still places the tightest constraints on the
annihilation crosssection. Our limits are much lower than those in
the 11month FermiLAT analysis by Ackermann et al. (2010), where the
tightest constraint came from Fornax for a much lower assumed subhalo
contribution of . Our limits are also tighter than that from
a joint analysis of the dwarf satellites of the Milky Way by
GeringerSameth & Koushiappas (2011).
The channel
As have been seen in section 3.2.1, the EXT model places tighter constraints on the crosssection than the PT model, and is the fiducial model expected from recent simulations. Therefore from now on we will only show results for the EXT model. The flux and crosssection upper limits for DM annihilating through the channel are plotted in Fig. s 7 and 8. The predicted flux upper limits for Coma and Virgo are still comparable to the CR level, with Fornax having much lower CR emission. The inferred crosssection falls below the canonical value for DM particle masses less than 10 GeV. Note the discontinuity in the upper limits around 100 GeV which reflect the transition from the prompt annihilation dominated regime to the IC emission dominated regime in the photon spectrum.
The channel
3.3 Allowing for an undetected point source population
Although we have detected five new point sources in the 45month data in the region of our three clusters, it is still necessary to account for population of still undetected point sources. When no unknown point sources are present, the probability of measuring a certain value of when the null hypothesis is true is given by the probability that Poisson fluctuations in the photon counts for the null model exceed some value. When a population of undetected point sources is present, the Poisson fluctuations become correlated and it is easier for the same amplitude of fluctuations to result in a given value of . In this case, the distribution of no longer follows the distribution, as predicted by Wilk’s theorem, because the data no longer follow a pure Poisson distribution from which the likelihood function is constructed.
Allowing for the presence of undetected point sources in the data will lead to weaker upper limits. We obtain these by performing MonteCarlo simulations to recalibrate the significance corresponding to a given value of . In the simulations we include the GAL and 2FGL sources, but we split the EG component into two parts: a population of undetected point sources and a residual smooth EG component, such that the sum of the two is consistent with the standard EG component. We consider a benchmark model for the undetected point source population which is close to the model derived by Abdo et al. (2010) and which contributes 14% of the EG background. Details of the simulations may be found in Appendix B. A standard likelihood analysis is then performed on the simulated data in order to derive appropriate values of for an assumed DM or CR component.
With the introduction of the undetected point source population, the distribution function for is found to be roughly described by . That is, the significance of a given value of is approximately reduced by a factor of compared to the significance of the same value of in the absence of the undetected point source population. For a DM component, we find for Coma and for Virgo and Fornax. The factor is not sensentive to the adopted DM spectrum. For the CR models, we find for Coma and Fornax, and for Virgo.
In order to obtain new limits from the corrected TS, let us first consider the likelihood function that has been maximized over all the nuisance parameters. Expanding around the maximum likelihood value of the parameter to leading order, we have
(12) 
where is the Hessian maxtrix. The 95% upper limit is calculated from , so that
(13) 
Note that . Similar equations hold for the improved background model, with , and replaced by , and respectively, assuming that there is no bias in the bestfit parameters. The 95% upper limit is then corrected for the undetected point source fluctuations according to
(14) 
where and are the upper limit and likelihood ratio from the standard analysis, while and are the corrected upper limit and likelihood ratio. For , the increase in the upper limit is at most 70%.
In Fig. 10, we show the corrected dark matter annihilation crosssection upper limits adopting for Coma and for Virgo and Fornax. The corrected TS and upper limits for CR models are listed in the last two columns of Table 2.
4 Discussion and conclusions
We have performed maximum likelihood fits to the 3year FermiLAT data for three galaxy clusters: Coma, Fornax and Virgo. We fit models which, in addition to point sources and galactic and extragalactic backgrounds, include emission due to dark matter (DM) annihilation and cosmic rays (CR). For the former, we assume both a point source and the theoretically predicted extended distribution of gamma rays in three generic annihilation channels, the , and channels. When searching for a dark matter signal, we experiment with different treatments of the CR component. In the traditional Fermi analysis, the extragalactic background (EG) is assumed to be a smooth component. In this work we have also investigated a more realistic EG model where a fraction of the EG emission comes from a population of undetected point sources.
Performing a standard likelihood analysis we obtain the following results:

In all three clusters and for the four different treatments of CR we have implemented, no significant detection of DM emission is obtained. We set upper limits on the flux and crosssection of DM annihilation in the three clusters we have investigated. Uncertainties in the CR component have only a mild effect on the upper limits: for the different CR models, the DM upper limit constraints agree to within a factor of two.
Models in which the DM annihilation emission has the extended profile predicted by cosmological simulations (Gao et al., 2012) have higher flux upper limits than models in which this emission is assumed to be a point source. Due to the large luminosity enhancement, of order of 1000, by emission from subhalos, the upper limits on the annihilation crosssection for extended models are at least 100 times lower than those for point source models. Our crosssection constraints are much tighter than those from an analysis of clusters using the 11month data (Ackermann et al., 2010), mostly because we take into account the effect of subhalos. Our constraints are also tighter than those from a joint analysis of Milky Way dwarf galaxies (GeringerSameth & Koushiappas, 2011; Ackermann et al., 2011).
Our new limits exclude the thermal crosssection for GeV for and final states, and for GeV for final states. We note that the annihilation crosssection in dark matter halos need not be the standard thermal crosssection of supersymetric models. In cases where the crosssection is velocity dependent, for example, through pwave contributions at freezeout (see e.g., Jungman et al. (1996)), one can easily have a different average crosssection. We emphasize that there is still a large uncertainty our adopted annihilation profile, which depends on a significant extrapolation of the resolved subhalo population by more than 10 orders of magnitude in mass. Taking this into account, the thermal crosssection, however, could still be reconciled with the data by assuming a larger cutoff mass in the WIMP power spectrum, thus reducing the contribution from subhalos and hence the factor. Since the total enhancement from subhalo emission scales as (Springel et al., 2008), a cutoff mass of , rather than our assumed , would be sufficient to increase the crosssection limits by a factor of 3.

Assuming no DM annihilation radiation, the gamma ray data for Coma and Virgo already set significant constraints on the CR level. For Virgo, the data are consistent with the predictions of the analytic CR model proposed by Pinzke & Pfrommer (2010) and Pinzke et al. (2011) while, for Coma, the data place an upper limit that is a factor of two below the analytical prediction, indicating either an uncertainty in model parameters such as halo mass, gas density and maximum shock injection efficiency, , or a peculiarity of the CR emission in Coma. If attributed to , the upper limit on the normalization parameter, , translates into an upper limit on of 0.3, assuming a linear form for . This is consistent with the estimates obtained independently by Zimmer et al. (2011) for Coma using Fermi data and by the Aleksić et al. (2012) for the Persus cluster using MAGIC observations. If interpreted as an error in the halo mass, a reduction in mass by a factor of is required to reconcile the model with the upper limits, assuming a simple CR luminosity scaling relation, (Pinzke & Pfrommer, 2010), or a factor of according to Eqn. 16 in the case when the gas density profile is fixed from Xray observations. For Fornax, the zerosignificance of a CR component is consistent with the low level predicted by the model.

Five new point sources with in Virgo and Fornax have been detected in the 45month data. Ignoring these new point sources results in a detection for a DM component in Virgo, in contrast to a detection when account is taken of these point sources.
In addition to the standard likelihood analysis, we have also investigated a model in which the EG component includes a population of undetected point sources whose numberflux relation extrapolates smoothly that of the detected sources. Using MonteCarlo simulations, we find that the standard Fermi likelihood analysis could overestimate the TS of extended emission by a factor of , and underestimate the upper limits by up to 70 percent. Adopting this more realistic EG model yields slightly looser upper limits, but does not quantitatively change any of the above conclusions. Still, it should be kept in mind that these corrections are derived from simulations assuming a particular distribution for the point source population. It is too computationally expensive to explore the parameter space of point source populations with MonteCarlo simulations. A more detailed and more general analytical study of the effect of undetected point sources will be presented elsewhere (Han et. al., in preparation).
In our analysis we have allowed the parameters of 2FGL point sources lying within the cluster virial radius to vary. This accounts for possible corrections to the 2FGL parameters in the presence of a DM or a CR component, while also avoiding the risk of refitting sources lying near the boundary of the data region with less accuracy. The parameters of highly variable sources are also kept free since the 2FGL parameters for these sources would be the average during a 2 year period whereas here we have 45 months of data. However, we also tried keeping all the point sources fixed or allowing the parameters of all the point sources within the data region to vary during the fitting. We find that this freedom in the treatment of the point sources has little impact on the DM model fits.
The cluster annihilation luminosity scales roughly linearly with halo mass, with the shape of the profile being almost independent of halo mass or concentration when expressed in terms of the normalized radius . We investigate the effect of mass uncertainties in Appendix E. We have also checked that the different energy cuts assumed in our analysis and in that of Huang et al. (2012) have no effect on the derived upper limits. We are able to reproduce the upper limits on the annihilation crosssection of Huang et al. (2012) for the test case of the Fornax cluster with 3year data, after adopting the same instrument response function and correcting for slightly different assumed subhalo contributions.
The CR model used in this analysis is still subject to improvement. This model is derived from simulations which, unavoidably, make simplifying assumptions. For example, the simulations only consider advective transport of CR by turbulent gas motions but there are other processes such as CR diffusion and streaming which may flatten the CR profiles (Enßlin et al., 2011). In particular, if the CR diffusion is momentum dependent this will entangle the spectral and spatial profile of CR and modify the morphology as well as the spectrum of the CR emission, thus invalidating our basic assumption that is the only free parameter. There could also be CR injected from AGN which are not accounted for in the current model.
Although we have not detected DM annihilation emission in our small cluster sample, the signaltonoise ratio can potentially be enhanced by stacking many clusters. Such an analysis was recently carried out by Huang et al. (2012), but the signaltonoise was degraded because of their assumption of an NFW annihilation profile. These authors considered an extended subhalodominated annihilation profile but only for individual clusters, not for the stack. Their stacked analysis placed looser constraints on DM annihilation emission than their analysis of individual clusters, presumably because the use of an inappropriate theoretical profile resulted in the different clusters yielding inconsistent results. Thus, it is clearly worth repeating the joint analysis with the “correct” subhalodominated profile. It is also tempting to extend the search for DM annihilation using multiwavelength data, from the radio to very high energy gammarays and even in the neutrino channel(Dasgupta & Laha, 2012), where different systematics are expected for different bands.
Acknowledgments
We thank Shaun Cole, Jie Liu, Yu Gao, John Lucey, Anders Pinzke, Christoff Pfrommer, Dan Hooper, Neal Weiner, Douglas Finkbeiner, Gregory Dobbler, Louie Strigari, Christoph Weniger, Savvas Koushiappas, and Fabio Zandanel for helpful discussions. JXH acknowledges the support on software issues from Tesla Jeltema and the Fermi science support team, especially Elizabeth C. Ferrara, Jeremy S. Perkins, Dave Davis and Robin Corbet. JXH is supported by the European Commissions Framework Programme 7, through the Marie Curie Initial Training Network CosmoComp (PITNGA2009238356), and partially supported by NSFC (10878001, 11033006, 11121062) and by the CAS/SAFEA International Partnership Program for Creative Research Teams (KJCX2YWT23). CSF acknowledges a Royal Society Wolfson research merit award and an ERC Advanced Investigator grant. The calculations for this work were performed on the ICC Cosmology Machine, which is part of the DiRAC Facility jointly funded by STFC, the Large Facilities Capital Fund of BIS, and Durham University. This work was supported in part by an STFC rolling grant to the ICC. The work of D.M. is supported in part from the SCOPES project IZ73Z0 128040 of Swiss National Science Foundation, grant No CM2032012 for young scientists of National Academy of Sciences of Ukraine, Cosmomicrophysics programme of the National Academy of Sciences of Ukraine and by the State Programme of Implementation of Grid Technology in Ukraine.
Appendix A Detection of New Point Sources
We model the new point sources assuming powerlaw spectra. For a given pixel, we calculate the value for an assumed new point source centered on that pixel. The calculation is performed using the binned method in the pyLikelihood tool, with a null model which includes the GAL and EG components and all the 2FGL sources within 15 deg of each cluster, but with the parameters of the 2FGL sources fixed. Around each cluster, we carry out a first scan of all the pixels within the cluster virial radius (and within 4 deg around Coma) using a pixel size of deg.
Regions with a peak are identified as potential locations of new point sources. We then scan each potential point source region using 10 times smaller pixels. The calculated map is then interpolated with cubic splines down to 0.002 degpixel. The value and location of the peak is taken as the and position for a new point source, if the peak . In case several peaks are clustered, we first extract the primary peak, then scan for lower peaks by including the newly detected sources into the null model. In our sample, no secondary peaks survive this iterative examination to be identified as new point sources.
The new point sources are listed in Table 3, and ploted in Figure 1 and C.1. Sources in Virgo and Fornax are prefixed by “V” and “F” in their names respectively. None of these new sources show significant variability when binned over monthly scale. The last column of Table 3 shows possible associations of astrophysical sources with these new detections, which are found to lie within the 2 confidence region of the detections.
Name  TS  RA (deg)  DEC (deg)  Flux ( ph cm s)  Spectral Index 
Seperation (deg) 
Possible Association 
V1  32.5  190.920  16.194  4.96  LBQS 1241+1624  
V2  31.8  185.698  11.116  2.31  [VV2006] J122307.2+110038  
V3  31.6  184.066  9.456  4.58  2MASX J12160619+0929096  
V4  30.5  185.894  8.286  4.42  SDSS J122321.38+081435.2  
F1  26.3  58.300  36.386  3.17  [VV98b] J035305.1362308 
Appendix B MonteCarlo simulation of undetected point source populations
To model the undetected population we adopt the following model based on the results of Abdo et al. (2010). Each point source is assumed to have a powerlaw spectrum defined by two parameters: flux and spectral index. The spectral index distribution is modeled as a Gaussian of mean and . The flux and spectral index are assumed to be independent. The differential number density of undetected point sources is assumed to be given by
(15) 
We adopt , and , as derived from Table 4 of Abdo et al. (2010). Since the total number of point sources diverges for , we cut off the flux distribution at . Due to the dependence of the detection efficiency on flux and spectral shape, there is no obvious cutoff in the maximum flux of undetected sources. We take as the detection threshold which corresponds to a detection completeness of , comparing 2FGL source counts and the model. This implies an undetected point source flux of 14% of the standard EG background, consistent with the results of Abdo et al. (2010). The synthetic spectrum of these undetected point sources is then subtracted from the standard EG template to yield a residual EG template for the simulation.
We perform 750 independent realizations of the 15 deg Virgo region in the presence of undected point sources. For each realization, we generate mock data in the following steps:

Generate a Poisson random number for the total number of undetected point sources within 15 deg.

For each point source, generate a random spectral index and a random flux according to the distributions specified above. Also, generate random coordinates for the point source according to a uniform distribution on the sky.

Feed these point sources and the 2FGL point sources within 15 deg, as well as the GAL and remaining EG components, to gtobssim.
The standard likelihood analysis is then applied to the simulated data without including any of the randomly generated point sources in the model. Here we only consider the CRonly and the DMonly models. In Fig. B.1 we show the cumulative probability distribution of values. Simple scaled versions of the standard distributions can roughly describe the distribution and provide the simplest way to convert the fitted to the standard distributed .
Appendix C Gammaray images for Coma and Fornax
In this Appendix we show gammaray images for the Coma and Fornax cluster regions. The corresponding image for Virgo is shown in Fig. 1.
Appendix D SemiAnalytic formula for the Cosmic Ray induced gammaray emission
Here we summarize the relevant equations for calculating the CR induced gammaray emission in galaxy clusters as derived by Pinzke & Pfrommer (2010) and Pinzke et al. (2011). The CR induced photon source function from pion decay can be decomposed as:
The spatial part is given by:
(16) 
with
(17)  
(18)  
(19)  
(20) 
The spectrum is given as:
(21)  
with , , . Here is the proton mass, the neutral pion mass and the speed of light. The maximum shock acceleration efficiency is chosen to be so that . The term describes the diffusive CR losses due to escaping protons as
(22) 
The proton cutoff energy is
(23) 
The energy is related to the photon cutoff energy, , through the momentum relation . The effective crosssection for protonproton interactions is given by:
(24) 
The gas density is fitted with multiple betaprofiles:
(25) 
where is the primordial hydrogen mass fraction and is the ratio of electron and hydrogen number densities in the fully ionized intracluster medium, with parameter values for , and listed in TABLE VI of Pinzke et al. (2011).
Appendix E Effect of mass uncertainties in Virgo
We adopt a virial mass for Virgo of , as estimated by Tully & Shaya (1984) from an analysis of the infall pattern of galaxies around Virgo. This value is consistent with other dynamical measurements (Smith, 1936; Hoffman et al., 1980; Tonry et al., 2000; Fouqué et al., 2001; Karachentsev & Nasonova, 2010). Mass estimates from Xray gas modelling tend to give somewhat lower values (Böhringer et al., 1994; Schindler et al., 1999; Urban et al., 2011). Thus, in addition to our adopted mass uncertainty from Tully & Shaya (1984), as an extreme case, we consider also a value of , obtained by scaling the Xray estimate to the virial radius (Urban et al., 2011). In Fig. E.1 we show the effect of adopting these different masses on the upper limits for DM annihilation in the channel. Since the flux upper limit is insensitive to slight changes in the profile shape and thus in the mass, while the luminosity (or integrated factor) scales linearly with mass, the crosssection upper limits are expected to be roughly inversely propotional to mass. This is indeed the case in Fig. E.1, where .
Footnotes
 See also a preliminary result by the FermiLAT collaboration (Vitale et al., 2009).
 We notice that several new point sources in Virgo are also identified in a concurrent paper(MacíasRamírez et al., 2012) and are found to reduce the significance of DMlike emission in the cluster, consistent with what we find here.
 http://fermi.gsfc.nasa.gov/cgibin/ssc/LAT/LATDataQuery.cgi
 We also tried using P7SOURCE_V6 IRF and Event Class 2 data. The results are consistent with those presented in this paper.
 ZENITH ANGLE (degrees): angle between the reconstructed event direction and the zenith line (originates at the centre of the Earth and passes through the centre of mass of the spacecraft, pointing outward). The Earth’s limb lies at a zenith angle of 113 degrees.
 EVENT CLASS: flag indicating the probability of the event being a photon and the quality of the event reconstruction.
 DATAQUAL: flag indicating the quality of the LAT data, where 1 = OK, 2 = waiting review, 3 = good with bad parts, 0 = bad
 LATCONFIG: flag for the configuration of the lat (1 = nominal science configuration, 0 = not recommended for analysis)
 ROCK ANGLE: angle of the spacecraft axis from the zenith (positive values indicate a rock toward the north).
 Angular diameter distance, from the NASA extragalactic database for Coma and Fornax, and from Tully & Shaya (1984) for Virgo.
 Cluster halo mass defined as the mass within the radius, , within which the average density equals 200 times the critical density of the universe. Values for Coma and Fornax are taken from Pinzke et al. (2011), while the value for Virgo is taken from Tully & Shaya (1984).
 Cluster halo mass defined as the mass within the radius, , within which the average density equals 200 times the critical density of the universe. Values for Coma and Fornax are taken from Pinzke et al. (2011), while the value for Virgo is taken from Tully & Shaya (1984).
 Integrated coefficient, , over the solid angle spanned by the cluster virial radius, assuming a smooth NFW density profile.
 Enhancement to the total annihilation luminosity within the virial radius due to substructures, extrapolated to a subhalo mass limit of . Note this factor scales with the minimum subhalo mass as (Springel et al., 2008).
 The LAT PSF size scales roughly as , so at 1 GeV it is deg
 http://www.darksusy.org.
 http://home.thep.lu.se/ torbjorn/Pythia.html
 Best fit normalization ( is the theoretical prediction)
 95% upper limit (UL) on the normalization
 TS after allowing for undetected point sources; see Section 3.3 for details
 Upper limit on the normalization factor after allowing for undetected point sources; see Section 3.3 for details
 In Coma, where the FiducialCR model is ruled out, the OptimalCR model yields the lowest upper limit.
 If systematic uncertainties in the halo mass parameters assumed by GeringerSameth & Koushiappas (2011) are considered, the lower bounds of their derived limits become comparable to our limits.
 The “Fermi1yr” constraint is only reproduced schematically, by reading out several data points from the original plot in the reference.
 Photon spectral index for .
 Distance to cluster centre.
References
 Aalseth C. E. et al., 2011a, Physical Review Letters, 106, 131301, arXiv:1002.4703
 Aalseth C. E. et al., 2011b, Physical Review Letters, 107, 141301, arXiv:1106.0650
 Abazajian K. N., Kaplinghat M., 2012, Phys. Rev. D, 86, 083511, arXiv:1207.6047
 Abdo A. A. et al., 2010, ApJ, 720, 435, arXiv:1003.0895
 Abdo A. A. et al., 2009, ApJ, 707, 55, arXiv:0910.3565
 Ackermann M. et al., 2011, Physical Review Letters, 107, 241302, arXiv:1108.3546
 Ackermann M. et al., 2010, J. Cosmology Astropart. Phys, 5, 25, arXiv:1002.2239
 Ahmed Z., et al., 2011, Physical Review Letters, 106, 131302, arXiv:1011.2482
 Aleksić J. et al., 2012, A&A, 541, A99, arXiv:1111.5544
 Angloher G. et al., 2012, European Physical Journal C, 72, 1971, arXiv:1109.0702
 Aprile E., et al., 2011, Physical Review Letters, 107, 131302, arXiv:1104.2549
 Baxter E. J., Dodelson S., 2011, Phys. Rev. D, 83, 123516, arXiv:1103.5779
 Bernabei R. et al., 2010, European Physical Journal C, 67, 39, arXiv:1002.1028
 Bertone G., Hooper D., Silk J., 2004, Phys. Rep., 405, 279
 Böhringer H., Briel U. G., Schwarz R. A., Voges W., Hartner G., Trümper J., 1994, Nature, 368, 828
 Boyarsky A., Malyshev D., Ruchayskiy O., 2011, Physics Letters B, 705, 165, arXiv:1012.5839
 Boyarsky A., Malyshev D., Ruchayskiy O., 2012, ArXiv eprints, arXiv:1205.4700
 Bringmann T., Huang X., Ibarra A., Vogl S., Weniger C., 2012, J. Cosmology Astropart. Phys, 7, 54, arXiv:1203.1312
 Ciafaloni P., Comelli D., Riotto A., Sala F., Strumia A., Urbano A., 2011, J. Cosmology Astropart. Phys, 3, 19, arXiv:1009.0224
 Cirelli M. et al., 2011, J. Cosmology Astropart. Phys, 3, 51, arXiv:1012.4515
 Collar J. I., 2011a, ArXiv eprints, arXiv:1103.3481
 Collar J. I., 2011b, ArXiv eprints, arXiv:1106.0653
 Dasgupta B., Laha R., 2012, ArXiv eprints, arXiv:1206.1322
 Diemand J., Kuhlen M., Madau P., 2007, ApJ, 657, 262, arXiv:astroph/0611370
 Duffy A. R., Schaye J., Kay S. T., Dalla Vecchia C., 2008, MNRAS, 390, L64, arXiv:0804.2486
 Enßlin T., Pfrommer C., Miniati F., Subramanian K., 2011, A&A, 527, A99, arXiv:1008.4717
 Finkbeiner D. P., 2004, ArXiv eprints, arXiv:astroph/0409027
 Fouqué P., Solanes J. M., Sanchis T., Balkowski C., 2001, A&A, 375, 770, arXiv:astroph/0106261
 Gao L., Frenk C. S., Jenkins A., Springel V., White S. D. M., 2012, MNRAS, 419, 1721, arXiv:1107.1916
 GeringerSameth A., Koushiappas S. M., 2011, Physical Review Letters, 107, 241303, arXiv:1108.2914
 Gondolo P., Silk J., 1999, Physical Review Letters, 83, 1719, arXiv:astroph/9906391
 Han J., Frenk C. S., Eke V. R., Gao L., White S. D. M., 2012, ArXiv eprints, arXiv:1201.1003
 Hoffman G. L., Olson D. W., Salpeter E. E., 1980, ApJ, 242, 861
 Hooper D., 2012, ArXiv eprints, arXiv:1201.1303
 Hooper D., Finkbeiner D. P., Dobler G., 2007, Phys. Rev. D, 76, 083012, arXiv:0705.3655
 Hooper D., Goodenough L., 2011, Physics Letters B, 697, 412, arXiv:1010.2752
 Hooper D., Linden T., 2011a, Phys. Rev. D, 83, 083517, arXiv:1011.4520
 Hooper D., Linden T., 2011b, Phys. Rev. D, 84, 123005, arXiv:1110.0006
 Huang X., Vertongen G., Weniger C., 2012, J. Cosmology Astropart. Phys, 1, 42, arXiv:1110.1529
 Jeltema T. E., Kehayias J., Profumo S., 2009, Phys. Rev. D, 80, 023005, arXiv:0812.0597
 Jungman G., Kamionkowski M., Griest K., 1996, Phys. Rep., 267, 195, arXiv:hepph/9506380
 Karachentsev I. D., Nasonova O. G., 2010, MNRAS, 405, 1075, arXiv:1002.2085
 Linden T., Hooper D., YusefZadeh F., 2011, ApJ, 741, 95, arXiv:1106.5493
 MacíasRamírez O., Gordon C., Brown A. M., Adams J., 2012, Phys. Rev. D, 86, 076004, arXiv:1207.6257
 Mattox J. R. et al., 1996, ApJ, 461, 396
 Nolan P. L., Abdo A. A., Ackermann M., Ajello M., Allafort A., et al., 2012, ApJS, 199, 31, arXiv:1108.1435
 Pinzke A., Pfrommer C., 2010, MNRAS, 409, 449, arXiv:1001.5023
 Pinzke A., Pfrommer C., Bergström L., 2011, Phys. Rev. D, 84, 123509, arXiv:1105.3240
 Schindler S., Binggeli B., Böhringer H., 1999, A&A, 343, 420, arXiv:astroph/9811464
 Smith S., 1936, ApJ, 83, 23
 Springel V. et al., 2008, Nature, 456, 73, arXiv:0809.0894
 Su M., Finkbeiner D. P., 2012, ArXiv eprints, arXiv:1206.1616
 Tempel E., Hektor A., Raidal M., 2012, J. Cosmology Astropart. Phys, 9, 32, arXiv:1205.1045
 Tonry J. L., Blakeslee J. P., Ajhar E. A., Dressler A., 2000, ApJ, 530, 625, arXiv:astroph/9907062
 Tully R. B., Shaya E. J., 1984, ApJ, 281, 31
 Urban O., Werner N., Simionescu A., Allen S. W., Böhringer H., 2011, MNRAS, 414, 2101, arXiv:1102.2430
 Vitale V., Morselli A., for the Fermi/LAT Collaboration, 2009, ArXiv eprints, arXiv:0912.3828
 Weniger C., 2012, J. Cosmology Astropart. Phys, 8, 7, arXiv:1204.2797
 Zaharijas G., Cuoco A., Yang Z., Conrad J., 2010, ArXiv eprints, arXiv:1012.0588
 Zimmer S., Conrad J., for the FermiLAT Collaboration, Pinzke A., 2011, ArXiv eprints, arXiv:1110.6863