Constraining Extended Gamma-ray Emission from Galaxy Clusters
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 gamma-ray backgrounds and point sources from the Fermi 2-year 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 Fermi-LAT data we detect several new point sources (not present in the Fermi 2-year 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 clusters
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 self-interaction cross-section 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 N-body simulations show that cold dark matter halos contain a population of self-bound 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 gamma-ray 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 gamma-ray 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 line-like 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 large-scale distribution of dark matter while avoiding some of the uncertainties arising from the astrophysical modelling of galactic gamma-ray 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 gamma-ray 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 DM-dominated 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 (Geringer-Sameth & Koushiappas, 2011; Ackermann et al., 2011) resulted in no significant detection but have began to rule out the canonical annihilation cross-section 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 NFW-squared profile could miss most of the signal. In fact, just such a search using the 11-month Fermi-LAT data has yielded no significnat detection of emission from six clusters (Ackermann et al., 2010).
Using the 45-month data, we consider possible contributions from cosmic ray (CR) induced gamma-ray 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 semi-analytic 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 signal-to-noise 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 three-year 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 2-year 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 cross-section 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 non-standard 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
|Minimum Energy||100 MeV|
|Maximum Energy||100 GeV|
|Maximum zenith angle
|ABS (ROCK ANGLE)
|ROI-based zenith angle cut||yes|
We list the basic properties of the three clusters in Table 1.
|Enhancement due to subhalos within
1.2 Maximum-likelihood fitting
We use the pyLikelihood tool shipped with the Fermi Science
Tools software package (version v9r27p1-fssc-20120410) 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 degree-wide 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 best-fit 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,
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:
The factor comes from the constraint that the normalization parameter be non-negative. The significance of a detection can thus be quoted as (one-sided 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.
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 45-month 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 power-law 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 J-map 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 cross-section 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 number-flux relation which extrapolates smoothly from that measured for brighter sources. In this case the photon counts within a given pixel are no longer Poisson-distributed 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 Monte-Carlo simulations to re-evaluate the distribution of for the more realistic background model and provide corrections to the results of the standard analysis where needed.
2 Modeling gamma-ray emission in clusters
We model the observed gamma-ray 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 CR-induced 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 Fermi-LAT data server, while the point sources include those from the LAT 2-year point source catalogue, 2FGL (Nolan et al., 2012), as well as several, newly detected by us, in the 45-month data. In addition, an improved EG model which includes a population of un-detected sources is also analyzed. We now describe in detail our models for DM annihilation and CR emission.
2.1 Dark matter annihilation emission
The gamma-ray intensity along the line-of-sight due to DM annihilation is given by:
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 velocity-averaged cross-section (or annihilation rate) for that channel, which is predicted to be constant in the low velocity limit appropriate to present-day cold DM particles (see e.g., Jungman et al. (1996)). The line-of-sight integration of the density squared is often expressed in terms of a dimensionless factor,
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 Fermi-LAT 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
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 mass-concentration relation:
(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
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:
Below we fit the subhalo emission surface brightness beyond the virial radius and extrapolate to several times the virial radius using an exponential decay,
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 Cosmic-ray induced gamma-ray emission within clusters
The cosmic ray induced gamma-ray emission is calculated following a semi-analytic prescription, derived from high resolution numerical simulations of galaxy clusters, that models cosmic ray physics self consistently (Pinzke & Pfrommer, 2010). The gamma-ray photon production rate (or source function) from pion decay is found to be separable into a spatial and a spectral part:
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 power-law 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 X-ray observations in the Appendix.
The differential gamma-ray flux from this source function, , is simply the integral of along the line-of-sight. 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,
We take as our fiducial CR model and also consider the case when is fitted from the actual gamma-ray 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 cross-section, . 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.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 “CR-only” 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 CR-only 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)||
3.2 Constraints on DM annihilation
Given the low significance of the CR detection in the CR-only 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:
- Fiducial-CR 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.
- Optimal-CR model.
The CR level is taken as the best-fit value listed in Table 2.
- Free-CR model.
The normalization of the CR level is left as a free parameter in the fit.
- No-CR 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 Free-CR model and to less than in the Fiducial-CR model.
The value of for the no-CR 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 2-year 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 power-law, rather than a DM annihilation spectrum, result in a similarly high significance detection, , and a best-fit 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 cross-section.
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 No-CR and Fiducial-CR
The flux upper limits are translated into cross-section 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 cross-section 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 cross-section must be smaller for low mass particles. With an enhancement of order due to subhalos, a much lower cross-section 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 cross-section limits drop below the fiducial thermal cross-section
of for GeV. Of the three clusters, Virgo has the highest flux upper
limits but it still places the tightest constraints on the
annihilation cross-section. Our limits are much lower than those in
the 11-month Fermi-LAT 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
Geringer-Sameth & Koushiappas (2011).
As have been seen in section 3.2.1, the EXT model places tighter constraints on the cross-section 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 cross-section 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 cross-section 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.
3.3 Allowing for an undetected point source population
Although we have detected five new point sources in the 45-month 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 Monte-Carlo simulations to re-calibrate 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
where is the Hessian maxtrix. The 95% upper limit is calculated from , so that
Note that . Similar equations hold for the improved background model, with , and replaced by , and respectively, assuming that there is no bias in the best-fit parameters. The 95% upper limit is then corrected for the undetected point source fluctuations according to
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 cross-section 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 3-year Fermi-LAT 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 cross-section 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 cross-section for extended models are at least 100 times lower than those for point source models. Our cross-section constraints are much tighter than those from an analysis of clusters using the 11-month 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 (Geringer-Sameth & Koushiappas, 2011; Ackermann et al., 2011).
Our new limits exclude the thermal cross-section for GeV for and final states, and for GeV for final states. We note that the annihilation cross-section in dark matter halos need not be the standard thermal cross-section of supersymetric models. In cases where the cross-section is velocity dependent, for example, through p-wave contributions at freeze-out (see e.g., Jungman et al. (1996)), one can easily have a different average cross-section. 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 cross-section, 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 cut-off mass of , rather than our assumed , would be sufficient to increase the cross-section 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 X-ray observations. For Fornax, the zero-significance 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 45-month 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 number-flux relation extrapolates smoothly that of the detected sources. Using Monte-Carlo 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 Monte-Carlo 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 cross-section of Huang et al. (2012) for the test case of the Fornax cluster with 3-year 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 signal-to-noise ratio can potentially be enhanced by stacking many clusters. Such an analysis was recently carried out by Huang et al. (2012), but the signal-to-noise was degraded because of their assumption of an NFW annihilation profile. These authors considered an extended subhalo-dominated 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” subhalo-dominated profile. It is also tempting to extend the search for DM annihilation using multi-wavelength data, from the radio to very high energy gamma-rays and even in the neutrino channel(Dasgupta & Laha, 2012), where different systematics are expected for different bands.
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 Cosmo-Comp (PITNGA-2009-238356), and partially supported by NSFC (10878001, 11033006, 11121062) and by the CAS/SAFEA International Partnership Program for Creative Research Teams (KJCX2-YW-T23). 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 CM-203-2012 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 power-law 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
Appendix B Monte-Carlo 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 power-law 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
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 CR-only and the DM-only 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 Gamma-ray images for Coma and Fornax
In this Appendix we show gamma-ray images for the Coma and Fornax cluster regions. The corresponding image for Virgo is shown in Fig. 1.
Appendix D Semi-Analytic formula for the Cosmic Ray induced gamma-ray emission
Here we summarize the relevant equations for calculating the CR induced gamma-ray 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:
The spectrum is given as:
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
The proton cut-off energy is
The energy is related to the photon cut-off energy, , through the momentum relation . The effective cross-section for proton-proton interactions is given by:
The gas density is fitted with multiple beta-profiles:
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 X-ray 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 X-ray 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 cross-section upper limits are expected to be roughly inversely propotional to mass. This is indeed the case in Fig. E.1, where .
- See also a preliminary result by the Fermi-LAT collaboration (Vitale et al., 2009).
- We notice that several new point sources in Virgo are also identified in a concurrent paper(Macías-Ramírez et al., 2012) and are found to reduce the significance of DM-like emission in the cluster, consistent with what we find here.
- 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.
- DATA-QUAL: flag indicating the quality of the LAT data, where 1 = OK, 2 = waiting review, 3 = good with bad parts, 0 = bad
- LAT-CONFIG: 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://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 Fiducial-CR model is ruled out, the Optimal-CR model yields the lowest upper limit.
- If systematic uncertainties in the halo mass parameters assumed by Geringer-Sameth & Koushiappas (2011) are considered, the lower bounds of their derived limits become comparable to our limits.
- The “Fermi-1yr” 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.
- 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 e-prints, 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 e-prints, arXiv:1103.3481
- Collar J. I., 2011b, ArXiv e-prints, arXiv:1106.0653
- Dasgupta B., Laha R., 2012, ArXiv e-prints, arXiv:1206.1322
- Diemand J., Kuhlen M., Madau P., 2007, ApJ, 657, 262, arXiv:astro-ph/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 e-prints, arXiv:astro-ph/0409027
- Fouqué P., Solanes J. M., Sanchis T., Balkowski C., 2001, A&A, 375, 770, arXiv:astro-ph/0106261
- Gao L., Frenk C. S., Jenkins A., Springel V., White S. D. M., 2012, MNRAS, 419, 1721, arXiv:1107.1916
- Geringer-Sameth A., Koushiappas S. M., 2011, Physical Review Letters, 107, 241303, arXiv:1108.2914
- Gondolo P., Silk J., 1999, Physical Review Letters, 83, 1719, arXiv:astro-ph/9906391
- Han J., Frenk C. S., Eke V. R., Gao L., White S. D. M., 2012, ArXiv e-prints, arXiv:1201.1003
- Hoffman G. L., Olson D. W., Salpeter E. E., 1980, ApJ, 242, 861
- Hooper D., 2012, ArXiv e-prints, 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:hep-ph/9506380
- Karachentsev I. D., Nasonova O. G., 2010, MNRAS, 405, 1075, arXiv:1002.2085
- Linden T., Hooper D., Yusef-Zadeh F., 2011, ApJ, 741, 95, arXiv:1106.5493
- Macías-Ramí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:astro-ph/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 e-prints, 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:astro-ph/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 e-prints, 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 e-prints, arXiv:1012.0588
- Zimmer S., Conrad J., for the Fermi-LAT Collaboration, Pinzke A., 2011, ArXiv e-prints, arXiv:1110.6863