SEARCH FOR EXTENDED GAMMA-RAY EMISSION FROM THE VIRGO GALAXY CLUSTER WITH FERMI-LAT

Galaxy clusters are the largest gravitationally bound objects in the universe. Observations of galaxy clusters, in particular estimations of the gravitation potential based on measurements of the velocity dispersions of member galaxies, suggested that Galaxy cluster masses were much larger than the values inferred by summing contributions from luminous matter, and lead to the proposal of non-luminous (i.e., dark) matter (DM) (Zwicky 1937). Since then, many observations, as well as N-body cosmological simulations, suggest that galaxy clusters contain a high amount of DM, making them prime targets to search for indirect DM signals (see, e.g., Clowe et al. 2006; Colafrancesco et al. 2006; Jeltema et al. 2009; Pinzke et al. 2009; Klypin et al. 2011; Pinzke et al. 2011; Sánchez-Conde et al. 2011; Gao et al. 2012; Hellwing et al. 2015).

Weakly interacting massive particles (WIMPs) constitute promising particle DM candidates. Among prominent WIMP candidates is the neutralino, which in many supersymmetric models is the lightest stable supersymmetric particle, allowing it to account for the observed relic DM density in the universe. In many of these models, the neutralino can self-annihilate into particle–antiparticle pairs which subsequently may produce other particles, including γ-rays. The γ-rays will propagate undeflected by the interstellar magnetic fields and thus reveal the location of the DM annihilation (see, e.g., Bertone et al. 2005; Feng 2010; Bringmann & Weniger 2012; Conrad et al. 2015, for a review of searches for indirect DM searches using γ-rays). Several predictions of the expected DM annihilation rate in cosmological environments, such as galaxy clusters, and the associated signals in γ-ray data show that current space-borne γ-ray detectors like the Large Area Telescope (LAT) on board the Fermi satellite (Ackermann et al. 2009) will be unlikely to detect this signal in case the DM is smoothly distributed (see, e.g., Pinzke et al. 2011).

However, a smooth DM distribution is indeed not expected and recent cosmological N-body simulations predict instead that DM virialized regions, known as halos, contain a large number of smaller virialized and highly concentrated substructures called sub-halos (Springel et al. 2005, 2008a; Diemand et al. 2008). Since the DM annihilation signal is proportional to the DM density squared, these highly concentrated sub-halos are expected to significantly boost the annihilation signal relative to the purely smooth DM scenarios (e.g., Kuhlen et al. 2008; Lavalle et al. 2008; Pieri et al. 2008; Springel et al. 2008a; Martinez et al. 2009). The exact signal enhancement depends on the abundance, distribution, and internal structural properties of the sub-halos. However, sub-halo properties are uncertain below the halo mass resolution of state-of-the-art N-body cosmological simulations, ∼O(10⁵ ${M}_{\odot }$) for Milky Way size halos (Diemand et al. 2008; Springel et al. 2008b; Hellwing et al. 2015) and 10⁸ ${M}_{\odot }$ for simulations of galaxy clusters (e.g., Gao et al. 2012). Thus, extrapolations of the relevant properties are required over several orders of magnitude in halo mass below the mass resolution limit in order to account for the whole halo mass range that is predicted to exist in the universe, and more specifically in clusters. Here, the fractional enhancement of γ-ray flux due to sub-halos is called the sub-halo boost, or boost-factor for short. As recently discussed, e.g., in Sánchez-Conde & Prada (2014), sub-halo boosts are very sensitive to the way these extrapolations are performed and boost estimates can vary drastically depending on the assumptions (e.g., Kamionkowski et al. 2010; Pinzke et al. 2011; Sánchez-Conde et al. 2011; Gao et al. 2012; Kuhlen et al. 2012; Nezri et al. 2012; Zavala & Afshordi 2014). A debate is ongoing as to whether the extrapolation with a power law to lower-mass halos (Pinzke et al. 2011; Gao et al. 2012) is justified or too optimistic (see, e.g., Diemand et al. 2008; Sánchez-Conde et al. 2011). ⁶²

A further challenge arises from the fact that cosmic-ray (CR) interactions in the intracluster medium (ICM) may also give rise to MeV–GeV γ-rays, which if observed, would be difficult to distinguish from a DM-induced signal (see, e.g., Pinzke et al. 2011; Ando & Nagai 2012). Despite intensive efforts, to date no γ-rays from clusters have been detected aside from those attributed to individual active galaxies (see, e.g., Ackermann et al. 2010a, 2014a; Huang et al. 2012; Huber et al. 2013).

The closest galaxy cluster is Virgo at a distance of about 15.4 ± 0.5 Mpc, subtending several degrees on the sky (Fouqué et al. 2001). The cluster consists of several sub-clusters which are located around giant elliptical galaxies, most prominently M87 and M49 as well as around the two smaller clusters associated with M100 and M60 (Schindler 2002). These sub-clusters are in the process of merging with one another, while the system is dominated by the most massive sub-cluster centered on M87. For the remainder of this work we refer to the sub-cluster centered on M87 as Virgo-I and the sub-cluster centered on M49 as Virgo-II. Also relevant for an analysis of γ-rays is the fact that M87 harbors a known active galactic nucleus (AGN, Abdo et al. 2009; Hada et al. 2014) which dominates the emission both in X-rays and γ-rays. Due to its proximity, Virgo is also an interesting target for the search of a DM-induced γ-ray signal with the Fermi-LAT. Earlier studies, concentrating on Virgo-I, tested for point-like or mildly extended (≤12) emission toward the center of the cluster and yielded upper limits on the integrated γ-ray flux of 14.1 × 10⁻⁹ ph cm⁻² s⁻¹ and 17.1 × 10⁻⁹ ph cm⁻² s⁻¹ respectively (Ackermann et al. 2010a), assuming a power-law spectrum of the γ-ray emission with photon index Γ = 2 in an energy band from 200 MeV to 100 GeV.

Claims of γ-ray emission induced by DM annihilation were put forward by Han et al. (2012b), using a very extended DM-induced emission profile considering only Virgo-I. ⁶³ Later studies attribute this putative signal to an incomplete point source model of the region (Han et al. 2012a; Macías-Ramírez et al. 2012).

Here we present a comprehensive analysis of the Virgo region searching for very extended emission and discuss the various systematic effects relevant for this analysis. In Section 2 we discuss our data selection and analysis. Sections 3 and 4 elaborate on the details of finding and characterizing extended excess emission, while Section 5 is devoted to the discussion of the uncertainties associated with the interstellar emission model (IEM). The aforementioned possibility of additional point sources constituting a putative signal is discussed in Section 6. As a result of these studies we devise an improved background model for the Virgo region. We use this model to derive new limits on the WIMP DM annihilation cross section (Section 7) and, in the case of CR-induced γ-ray production, on its flux. We compare our results to theoretical predictions in order to constrain relevant CR quantities (Section 8). We conclude and summarize our work in Section 9.

The main instrument on board the Fermi satellite, the LAT, is a pair-conversion telescope sensitive to γ-rays in the energy range from 20 MeV to >300 GeV. For a more detailed description the reader is referred to Ackermann et al. (2009) and for the characterization of the on-board performance to Ackermann et al. (2012).

We analyzed archival Fermi-LAT data between MJD 54682.7 (2008 August 04) and MJD 55789.5 (2011 August 16) corresponding to roughly three years of Pass 7 data. We chose Pass 7 data because the predicted γ-ray signal from possible DM annihilation stretches over several degrees in the case of the Virgo cluster. With the release of the reprocessed P7 data (P7REP), the LAT collaboration also released a new template to describe the Galactic foreground emission which is tuned to P7REP data. This model contains a component which is derived from re-injecting residual emission above a scale radius of about 2° (Casandjian, J.-M. for the Fermi-LAT Collaboration 2015). Residuals larger than this scale radius are absorbed into the model and thus the model cannot be used to search for emission with larger extension. The release of P7REP data also implied a switch in the data processing pipeline limiting the available Pass 7 data to the three years we used in this analysis. The reprocessing primarily impacts the energy reconstruction as it accounts for a time-dependent change of the calorimeter calibration constants (Bregeon & others for the Fermi-LAT collaboration 2013), which in case of a signal may result in a moderate (2%–3%) shift toward higher photon energies with respect to Pass 7 data as well as an improvement in the high energy LAT point-spread function (PSF), which however does not constitute any significant impact on our analysis. We also would like to point out that the recently released Pass 8 event selection provides ∼25% increased effective area above 1 GeV (Ackermann et al. 2015a) and increases the data live time by a factor of two compared to the data set used in our analysis of the Virgo region. We will show in subsequent sections that our findings and interpretation are limited by systematic uncertainties in the IEM. These uncertainties can not be overcome by statistics but require better understanding of the CR distribution and its interaction in our Galaxy which is used to derive the IEM. Moreover, the IEM released with Pass 8 contains a data-driven residual component at the scale of 2°. Hence, using this model is not suitable in the case of the spatially extended Virgo cluster.

The data were processed using the Fermi ScienceTools version v9r28p0. ⁶⁴ We selected events with high probability of being γ-rays by choosing the Source event class. In order to evade γ-ray contamination generated by CRs interacting with the atmosphere of the Earth, we removed events with an LAT zenith angle >100°. We excised time periods around bright solar flares and gamma-ray bursts and applied a rocking angle cut of 52°. Furthermore, we restricted our analysis to the 100 MeV to 100 GeV energy range and used the P7SOURCE_V6 instrument response functions.

To model the Galactic foreground emission caused by CRs interacting with the gas and radiation fields in our Galaxy, we use the gll_iem_v02.fit model. This IEM is the standard IEM provided by the Fermi-LAT collaboration for point source analysis in the un-reprocessed flavor. The isotropic γ-ray emission is accounted for by the isotropic_iem_v02.txt model. For simplicity we refer to this set of IEM and isotropic model as our standard model. We chose a 20° × 20° region of interest (ROI) centered on the center of Virgo-I (α₂₀₀₀ = 18771 and δ₂₀₀₀ = 1239) and performed a binned likelihood analysis with 01 spatial bins and 30 logarithmic bins in energy. ⁶⁵ A counts map of the ROI together with the position of all sources from the LAT 2-year catalog (Nolan et al. 2012) sources is shown in Figure 1. In addition to the two aforementioned diffuse model components, our background model contains all sources within a 30° radius around the Virgo-I center that are listed in the 2FGL catalog (Nolan et al. 2012). The Virgo ROI contains mostly extragalactic sources which may be variable and thus the two-year source parameters in the 2FGL catalog might be bad approximations for the three-year data. Another challenge arises by performing an analysis down to 100 MeV. At this energy the Fermi-LAT PSF with 68% containment radius is about 7° (Ackermann et al. 2012) and thus even far away but strong sources which are not modeled correctly might easily increase the significance of a very extended profile located at the cluster center. To account for this we free the normalization and spectral index of all sources within 5° from the center coordinates of either Virgo-I or II in addition to the bright sources in the ROI (see below). The sources left free in our fit are marked by crosses in Figure 1.

Zoom In Zoom Out Reset image size

Using a likelihood analysis we construct a test statistic (TS) following Mattox et al. (1996) to evaluate the improvement of the likelihood fit to the ROI when adding a new source by defining $\mathrm{TS}=-2(\mathrm{log}\;{{\mathcal{L}}}_{0}-\mathrm{log}\;{\mathcal{L}}).$ ${{\mathcal{L}}}_{0}$ refers to the maximum likelihood value for the null hypothesis and ${\mathcal{L}}$ for the value for the alternative hypothesis (including an additional source such as the cluster itself). In the case of one additional degree of freedom the significance can be written as $\sigma =\sqrt{\mathrm{TS}}$ (see Section 4 for details regarding the gauging of this quantity). We look for new point-like excesses and derive TS maps by placing a test point source at the location of each pixel of the map and maximizing the likelihood using the gttsmap tool.

To estimate the systematic uncertainties on source significance and flux properties caused by our limited knowledge of the IEM, we compare the standard IEM results to the results obtained with eight alternative models seeded with γ-ray emission maps generated by GALPROP ⁶⁶ (Vladimirov et al. 2011) and additional templates (as described in detail in Section 5) (de Palma & others for the Fermi-LAT collaboration 2013). We would like to stress that these models do not cover the entire systematic uncertainty associated with the IEM and the results are not expected to bracket the standard IEM results. Despite the potentially small coverage, using these models demonstrates the influence the IEM has on our result.

Instead of starting our search for extended emission from the cluster center with a specific physical model, we consider a more generic model for any potential excess, namely a simple uniform disk with a given radius, r_disk, centered on Virgo I. The uniform disk profile is very successful in finding weak extended emission and is less prone to degeneracy with strong point sources such as M87. Moreover, a disk profile is usually sufficient to find sources of various shapes because, even for very strong sources, a discrimination between different emission profiles is usually not possible (see, e.g., Lande et al. 2012). For these reasons, we use the disk profile with a power-law spectral model for all our systematic studies of the extended excess (Sections 3 and 4). To further characterize the extended excess, we perform a TS versus radius scan. The scan is performed in 05 steps of r_disk and shows a peak at r_disk = 3° with an associated TS-value of 14.2 as shown in Figure 2. This finding is in agreement with Han et al. (2012b) who attribute most of their found emission to the innermost 3° of their profile.

Zoom In Zoom Out Reset image size

Next, we considered the case of a DM-induced signal as proposed by Han et al. (2012b) which would lead to a γ-ray contribution peaked at the cluster center, again considering Virgo-I as the center of the cluster. To this end we substitute the disk profile with the DM profile used in Han et al. (2012b) in our source model of the Virgo region. We reproduce their results, finding TS ≃ 23, only if we fix the bright sources outside of 5° of Virgo-I to their catalog values as reported in 2FGL. However, when performing the analysis according to the description in Section 2, the TS-value drops considerably (TS ≃ 17).

We scan the inner 10° × 10° of the Virgo region with a disk with a radius of r_disk = 3° on a 05 grid to find the best position for the origin of the putative extended excess. To obtain this TS map we add a disk-like test source at each grid position and maximize the likelihood using gtlike. We leave all sources within a radius of 5° from the position of the test source and all bright sources (TS > 1000) free to vary in normalization and spectral index. This approach is needed since in the presence of a real extended source, the degeneracy between the extended source and overlapping or very close-by point sources will decrease the significance of the extended sources. The resulting TS map has peaks significantly offset from the centers of Virgo-I and Virgo-II. All of them are located in the lower right quadrant of our TS map shown in Figure 3. To obtain a more detailed description we made a finer grid of 02 spacing in the 4° × 4° region encompassing the highest TS values. This finely binned TS map is shown together with the coarse map in Figure 3. From the fine map it is evident that a large region of the Virgo ROI yields TS values above 25. In particular, there are two broad maxima seen which are spatially distinct from one another. Note that the typical 1σ localization contour for a source near threshold is about 01. Consequently, the peaks shown in Figure 3 appear to be larger than what would be expected from a point source. Repeating the study of r_disk versus TS at the positions of the two maxima we again find a clear peak at r_disk = 3° and thus continue all subsequent analysis using a disk with r_disk = 3°.

Zoom In Zoom Out Reset image size

In order to study the relationship between TS-values and significance, we performed an extended source search at 288 randomly selected sky positions (blank fields) to estimate the significance of finding a disk-shaped excess with r_disk = 3°. All test positions were selected so that their inner 5° do not overlap and hence each position can be treated as statistically independent. Furthermore, we only used positions with $| b| \gt 20^\circ$ and therefore exclude regions with bright Galactic γ-ray emission. The resulting TS distribution is shown in Figure 4 and is reasonably described by a χ² distribution with one degree of freedom. We note that the number of tested blank fields is too small to sample any probability density function at TS > 8. However, the absence of TS > 16 along with the fact that the majority of blank fields yield TS < 8 indicate that the significance of the excess we find here (up to TS ∼ 32) is larger than 3σ. We stress that this particular statement is only valid at high Galactic latitude from where we extracted our blank fields.

Zoom In Zoom Out Reset image size

The interstellar γ-ray emission is the dominant γ-ray background in Fermi-LAT analyses in the vicinity of the Galactic plane, while far from the Galactic plane the largest γ-ray contribution usually stems from discrete sources and isotropic γ-ray emission. However, in certain high-latitude regions the interstellar emission can contribute significantly to the γ-ray background and one of these regions is Virgo. Several features present in the IEM are contained within the Virgo ROI. The most striking features are the spatially uniform patch that is used to model the γ-ray emission from Loop I (Large et al. 1962) and some filaments of H i around the center of the Virgo cluster as shown in Figure 5. ⁶⁷ Loop I is only approximately modeled in the standard IEM and could influence our analysis of the region. For most point sources in the Virgo region the IEM does not play a significant role since it contributes on average only about 110 counts per square degree above 100 MeV. Yet, it must be considered carefully when analyzing weak and very extended emissions like our disk model in the previous section (accounting for about 800 predicted photons and about 28 counts per square degree). The isotropic γ-ray emission has the largest contribution to the diffuse γ-ray emission in the ROI and accounts for about 200 counts per square degree. We note that the spatial variation in the IEM γ-ray counts are up to a factor of two from the minimum value. These spatial variations in the IEM are about three times larger than the potential γ-ray counts from the disk profile. The weak disk emission could easily be caused by a missing feature in the IEM. Such a missing component would necessarily not be traceable by tracers of the interstellar gas (21 cm H i line and 2.6 mm CO line), and makes the systematic evaluation of IEM influences on our analysis mandatory.

Zoom In Zoom Out Reset image size

To assess the influence of the uncertainty in our knowledge of the IEM, besides the standard IEM, we use eight additional models. These models are seeded with simulated γ-ray intensities obtained from GALPROP, assuming CR halo heights of 4 and 10 kpc, two different CR source populations, referred to as "pulsars" and "supernova remnants" (SNRs), and two different H i spin temperatures (optically thin and 150 K). Afterwards each model's spectral component is multiplied with a log-parabola function to provide better model-data agreement and thus is not a direct output of GALPROP. For further details on the alternative models see de Palma & others for the Fermi-LAT collaboration (2013). Most of the parameters varied in the alternative models are only expected to have very slight influence on the results for the Virgo region since they should not influence the local CR density very much and that is what is mostly sampled when looking in the direction of Virgo. For each of these alternative models an individual modified isotropic diffuse contribution is used.

Each of these models is comprised of maps that are inferred from H i and CO tracers along with components modeling the large-scale diffuse residual structures, such as Loop I (Casandjian & others for the Fermi-LAT collaboration 2009) and the Fermi bubbles (Su et al. 2010; Ackermann et al. 2014c). The H ii emission that is inferred from GALPROP and based on the NE2001 model is added to the H i map. In addition each model also includes a model of emission from inverse Compton (IC) scattering of CR electrons on interstellar radiation also calculated by GALPROP. For these models the spectral line shifts of the H i and CO lines were used to derive maps for separate ranges of Galactocentric distance. Inside the Virgo ROI only H i is present and is located within the second and third ring H i templates (these rings correspond to Galactocentric distance ranges 4–8 kpc and 8–10 kpc, respectively). Besides these two H i components only the IC model and the Loop I template are included in our fit. Each normalization of the alternative IEMs components is left free in the likelihood fit. The second H i ring only contributes to the Virgo ROI because of the H ii emission that is added to the H i template. We note that there is considerable uncertainty in the estimation of the H ii emission. The Loop I template is a geometrical template adapted from Wolleben (2007) and based on a polarization survey at 1.4 GHz and it is modeled as two expanding shells centered on two local OB associations. Our standard IEM on the other hand uses a uniform-patch Loop I template whose shape was derived by visual inspection of the gamma-ray residuals when building the standard IEM. ⁶⁸ There is currently no template of Loop I available from observations at other wavelengths that adequately traces the gamma-ray emission observed by Fermi-LAT in the direction of Loop I. By following two different approaches to define a Loop I template into our models we can gain some insights into the influence of the Loop I modeling on our results.

We randomly select one position from the bottom panel of Figure 3 where TS ≥ 25 and repeat our likelihood calculations using these alternative diffuse models.

The results do not have consistent TS values between the individual IEMs as can be seen in Table 1. While the standard, I, V, and VII IEMs show high TS values for an extended source, the other models yield considerably lower values. Leaving the individual components of the IEM free to vary in the fit allows for a higher sensitivity to features in the IEM which might affect only one component map. Such a mis-modeling of a single component in the IEM could cause an extended excess roughly corresponding in shape with the mis-modeled or missing emission.

Table 1. Parameters for the Nine Interstellar Emission Models Used

Model	Sources	Halo height	TS	$\mathrm{log}({L}_{0})$	$\mathrm{log}({L}_{\mathrm{disk}})$	TSdisk	F(E > 1 GeV)	Γ
		(kpc)	(K)				× 10−9 cm−2 s−1
I	Pulsars	10	105	−342483.4	−342472.1	21.7	1.5 ± 0.1	1.85 ± 0.02
II	Pulsars	4	105	−342484.0	−342476.3	15.4	1.3 ± 0.3	1.8 ± 0.1
III	Pulsars	10	150	−342480.3	−342474.0	12.6	1.3 ± 0.3	1.8 ± 0.1
IV	Pulsars	4	150	−342481.3	−342474.0	14.6	1.3 ± 0.3	1.8 ± 0.1
V	SNR	10	105	−342481.9	−342471.3	21.2	1.5 ± 0.4	1.9 ± 0.1
VI	SNR	4	105	−342484.0	−342476.3	15.4	1.3 ± 0.3	1.8 ± 0.1
VII	SNR	10	150	−342479.2	−342469.4	19.7	1.5 ± 0.3	1.8 ± 0.1
VIII	SNR	4	150	−342481.1	−342474.1	14.1	1.2 ± 0.4	1.8 ± 0.1
Standard	...	...	...	−342494.3	−342480.4	27.4	1.7 ± 0.3	1.9 ± 0.1

Note. The models I–VIII have the normalization of each individual component left free in the fit while the standard model has fixed ratios between individual components. The columns contain, from left to right, the model name, the source population causing the γ-ray emission, the CR halo height, the spin temperature, the reference log-likelihood, the log-likelihood with the added disk, the γ-ray flux, and spectral index of the disk profile.

Download table as:  ASCIITypeset image

The large observed variation in TS-values of the disk emission for the alternative models provides an indication that the observed γ-ray emission may be due at least partially to an incomplete modeling of the IEM. In the direction of the Virgo cluster the interstellar gas is mostly local (within ∼1 kpc of the Sun) and thus only the local CRs should contribute to the IEM. The γ-ray emission caused by the local CR density is not expected to have large dependencies on the CR source distribution, CR halo height, or spin temperature and thus we would expect rather similar results for all models. However the alternative IEMs that have a relatively large TS for the disk are associated with large predicted photon counts, overemphasising the contribution from H i ring 2 by increasing its amplitude by a factor O(10–30). Such an implausibly large increase in the contribution is explained by the shape of the model components of each IEM. H i ring 2 covers a region in the projected sky that is similar to the Loop I template in the standard model and by upscaling the normalization beyond the physically viable bounds of the H i component the overall fit is improved. Note that the extension of H i ring 2 in the ROI depends on the CR halo height and thus introduces some dependency of the TS of the disk on this IEM parameter that is also not expected from the local CR density. The high normalization of the H i ring 2 demonstrates that some γ-ray emission is coming from its large region within the ROI that is not traced by the usual H i and CO tracers but can be partially compensated by overestimating their contribution in the ROI. The disk emission might just compensate a part of this large-area diffuse gamma-ray emission that is not overlapped by any other IEM component that could compensate it.

In general it is very difficult to obtain a conclusive picture as to the possible emission in the ROI. The IEM study suggests that while there may be some additional emission, the ROI contains several IEM components whose predicted emission overlap with the excess, making it difficult to disentangle their contribution in a log-likelihood fit to the γ-ray-data. A precise modeling of Loop I on the other hand poses a considerable challenge since its γ-ray emission is not well traced by the radio emission. Hence, we used a geometric model for Loop I in the alternative IEMs compared to the γ-ray residual-inferred template in the standard IEM. The low significance of the emission found in this work makes it extremely difficult to identify its origin due to the aforementioned issues but is likely caused by inaccuracies in the IEM. To describe the location of the extended emission we devised a double disk patch that accounts for the observed residual γ-ray emission. This double disk consists of two disks with r_disk = 3°, one at each of the maxima in our TS map shown in Figure 3 and a uniform single power law emission of the whole profile. In this way the model covers most of the extended emission and we can easily trace its position in sky maps and include it in all relevant plots discussing possible counter parts. In all further analysis when we are deriving upper limits on γ-ray emission from DM annihilation or CR interaction we leave the newly found extended emission unmodeled to be conservative as we can not make a definitive statement about the origin of the extended emission. This results in weaker upper limits compared to when including the aforementioned double-disk model for the emission when deriving upper limits.

Previously undetected (variable) point sources that are not contained in the 2FGL catalog provide a possibility to account for the extended excess we find. Such sources can be significant background sources in the three-year data set and should be included in the background model. This possibility was first addressed by Macías-Ramírez et al. (2012) and later by Han et al. (2012a). However, neither group searched for possible extended emission (not from DM annihilation) offset from the cluster center. Therefore we performed our own point source search in the Virgo ROI to be able to compare these results with the same analysis set up as our extended source search. We note that we search for new point source candidates without the double disk profile described in the previous section. This allows us to construct a new alternative background model including the newly found point source candidates but not considering any additional extended source profile.

We divide the Virgo ROI in 01 × 01 grid positions and fit all 2FGL sources in the region to obtain the log-likelihood reference value. We then add one point source at one grid point and calculate the TS of this source candidate and repeat this procedure at each grid point. We identify source candidates corresponding to grid points with TS > 15 and use them as seeds for the gtfindsrc localization algorithm. Finally we add all source candidates into our Virgo ROI model and fit the region again to obtain a better reference model for the extended source investigation.

Following this procedure, we find eight new source candidates, three of which appear clustered in close vicinity to one another. We note that when including any of the three, the remaining candidates are not statistically significant. Thus we only keep the brightest candidate in our background model and discard the other two. In summary, our new background model contains six new source candidates listed in Table 2. Four of the six sources lie in the vicinity of the center of Virgo-I and three of them are in the region where we find the extended emission (substantially offset from the center of Virgo-I). While adding these six source candidates to our model yields a much better description of the ROI, the closest candidate sources to the center of Virgo-I or Virgo-II are below the conventional source detection threshold (TS > 25). Note that five out of these six sources are contained in the latest four-year source catalog (3FGL, see Acero et al. 2015). To identify possible multi-wavelength counterparts to the γ-ray sources we searched in the 95% error circle around each source in the NASA/IPAC Extragalactic Database (NED) for potential γ-ray emitters. We do not find convincing counterparts for the new point source candidates; however, that does not perclude the possibility that they are real sources since a considerable fraction of Fermi-LAT sources are not associated with multi-wavelength counterparts (Nolan et al. 2012; Acero et al. 2015). The locations of the new source candidates are shown in the model map ⁶⁹ in Figure 6. While we do get slightly different results compared to the work of Macías-Ramírez et al. (2012) and Han et al. (2012a) the discrepancies can be attributed to the larger data set used in Macías-Ramírez et al. (2012) and Han et al. (2012a) and the slightly different analysis procedures, such as the number of source parameters freed in the search or the different definitions of the ROIs. Considering the very weak emission of all source candidates the consistency of the findings among these works is reasonable. We also note that when including the disk profile in addition to the new point source candidates in the ROI we obtain a TS value of only ∼7.5 for the disk emission.

Zoom In Zoom Out Reset image size

With this improved background model in hand, which contains new point source candidates but not the disk profile, we devote the remainder of this paper to study the Virgo cluster as a γ-ray emitter, either via DM annihilation or via CR interactions (Sections 7 and 8).

The γ-ray flux from annihilating DM particles of mass m_χ can be written as

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\gamma }(E,\psi )=\displaystyle \frac{1}{4\pi }\displaystyle \frac{\langle \sigma v\rangle }{2{m}_{\chi }}\displaystyle \sum _{j}{B}_{j}\displaystyle \frac{{{dN}}_{j}(E)}{{{dE}}_{j}}\times J(\psi ).\end{eqnarray} \tag{ 1 }$$

In the above equation, we have $\langle \sigma v\rangle$ as the thermally averaged product of DM self-annihilation cross section times velocity, and the sum runs over the final states the DM particles annihilate into with their specific γ-ray annihilation yields ${{dN}}_{j}(E)/{{dE}}_{j}$ and branching fraction B_j per final state j. We define the astrophysical J-factor as the line-of-sight integral of the squared DM density toward the observational direction, ψ, integrated over a solid angle ΔΩ:

$$\begin{eqnarray}&&J(\psi )={\displaystyle \int }_{{\rm{\Delta }}{\rm{\Omega }}}\left\{{\displaystyle \int }_{l.o.s.}{\rho }^{2}[l(\psi )]{dl}\right\}d{{\rm{\Omega }}}^{\prime }.\end{eqnarray} \tag{ 2 }$$

Density profiles ρ(r) for plausible DM distributions can be expressed in terms of a generalized Hernquist profile (Hernquist 1990):

$$\begin{eqnarray}&&\rho (r)={\rho }_{0}{\left(\displaystyle \frac{r}{{r}_{s}}\right)}^{-\gamma }{\left[1+{\left(\displaystyle \frac{r}{{r}_{s}}\right)}^{\alpha }\right]}^{\frac{\gamma -\beta }{\alpha }},\end{eqnarray} \tag{ 3 }$$

where α, β, and γ are shape parameters. High-resolution cosmological N-body simulations of cold DM halos indicate that their density profiles are well described by a Navarro–Frenk–White (NFW) profile where α = 1, β = 3, γ = 1 (Navarro et al. 1996, 1997). The quantity r_s is the characteristic scale radius of the profile with c as the concentration parameter such that r_s = r₂₀₀/c. r₂₀₀ is defined to be the virial radius for which within the cluster mass M₂₀₀ is: ${M}_{200}=4\pi {\rho }_{c}/3\times 200\times {r}_{200}^{3},$ where ρ_c is the critical density of the universe and ρ₀ is the characteristic density of the profile, ${\rho }_{0}=\displaystyle \frac{200}{3}\frac{{c}^{3}{\rho }_{c}}{\mathrm{ln}(1+c)-c/(1+c)}.$

As for the concentration–mass relation, we adopt here the one proposed by Prada et al. (2012) for the WMAP5 cosmology. This relation was derived from the ∼10 billion particle Bolshoi and MultiDark large scale structure N-body cosmological simulations (Klypin et al. 2011). ⁷⁰ By evaluating the integral ${M}_{500}={\displaystyle \int }_{0}^{{r}_{500}}{d}^{3}r\;\rho (r)$ with ρ(r) as defined in Equation (3), we can numerically determine the value for M₂₀₀ for which this equation is satisfied, using the reported values for M₅₀₀ and r₅₀₀ from Chen et al. (2007). We find M₂₀₀ = 5.6 × 10¹⁴ ${M}_{\odot }$ as the mass of Virgo-I with c = 4.21.

The integrated J-factor, for r < r₂₀₀ and assuming an NFW DM density profile, at a large angular diameter distance D_a from Earth to the center of the cluster can be approximated by:

$$\begin{eqnarray}&&{J}^{{\rm{NFW}}}=\displaystyle \frac{4\pi }{3}{\rho }_{0}^{2}{r}_{s}^{3}\displaystyle \frac{1}{{{D}_{a}}^{2}}.\end{eqnarray} \tag{ 4 }$$

7.1. Contribution from Substructure

Some fraction of the DM should reside in sub-halos within the NFW-like primary halo. The presence of sub-halos implies both a flattening of the surface brightness profile (see, e.g., Sánchez-Conde et al. 2011; Gao et al. 2012) and an enhancement (boost, denoted b) of the J-factor which may increase the total annihilation signal by orders of magnitudes. Here b = 0 corresponds to the case of the smooth NFW halo without the inclusion of additional substructure.

For the normalization of the DM substructure signal strength, we adopt a fiducial substructure model that follows the works by Sánchez-Conde & Prada (2014), assuming a moderate total boost factor of b = 33 as given by their proposed parametrization of the boost for the Virgo mass (DM-I). We contrast this conservative model with an optimistic model that implicitly adopts a power-law extrapolation of the mass–concentration relation to the smallest (unresolved in simulations) halo masses, yielding b = 1200 (DM-II, Gao et al. 2012).

For the spatial morphology of the DM-induced gamma-ray emission, including the predicted DM substructure signal, we adopt the form from a recent study of high resolution cosmological DM simulations of cluster-size halos (Gao et al. 2012). The projected luminosity profile from the substructure is approximated by

$$\begin{eqnarray}&&{I}_{\mathrm{sub}}(\theta )=\displaystyle \frac{16b\times {L}_{\gamma }^{\mathrm{NFW}}}{\pi \mathrm{ln}(17)}\displaystyle \frac{1}{{\theta }_{200}^{2}+16{\theta }^{2}},\quad \mathrm{for}\;\theta \leqslant {\theta }_{200},\end{eqnarray} \tag{ 5 }$$

where θ is the distance from the cluster center in degrees and θ₂₀₀ is the angle subtending the virial radius, given by ${\theta }_{200}=\mathrm{arctan}({r}_{200}/{D}_{a})\times 180^\circ /\pi ,$ and ${L}_{\gamma }^{\mathrm{NFW}}$ is the total γ-ray luminosity of the halo within the virial radius (Gao et al. 2012), defined as ${L}_{\gamma }=4\pi \times {{\rm{\Phi }}}_{\gamma }\times {D}_{a}^{2}.$

Note that the work by Sánchez-Conde & Prada (2014) does not address the change in the spatial morphology of the annihilation signal due to the presence of sub-halos. However, in previous works, Sánchez-Conde et al. (2011) have shown that moderate values for b also lead to a significant flattening of the annihilation profile in clusters. We found this flattening to agree reasonably well with the one implied by Equation (5) for moderate values of b. Thus, from here on, we assume this approximation to be a good representation of the spatial morphology of the DM substructure signal in both the DM-I and DM-II setups, with only the value of b differing from one to another substructure scenario. The resulting expected surface brightness profiles for DM annihilation are shown in Figure 7.

Zoom In Zoom Out Reset image size

Our choice of models for the substructure is motivated by assuming a common mass scale for both setups, M_cut = 10⁻⁶ ${M}_{\odot }$, at which the matter power-spectrum is truncated. Below this scale, no DM halos are formed and thus no halos will contribute to the expected DM signal. What the minimal DM halo size is and what properties these substructures have is, to a large extent uncertain and may depend on the specific DM particle model. The mass range for M_cut may vary from 10⁷ to 10⁻¹² ${M}_{\odot }$, where the upper limit comes the fact that we observe dwarf galaxies with that mass (McConnachie 2012). The lower limit is more uncertain and depends on the DM particle properties (Profumo et al. 2006; Bringmann 2009). ⁷¹

Recalling the discussion on the morphology of the cluster from Section 4, we remark that while the mass of Virgo-II is only about ∼13% that of Virgo-I, its concentration parameter, c = 5.58, is larger than that of Virgo-I (see Table 3). Since the DM annihilation flux is proportional to the third power of the concentration, the predicted DM-induced γ-ray flux of Virgo-II corresponds to about one-third of the predicted DM-induced γ-ray flux from Virgo-I. In order to account for this in our DM modeling, we consider the cluster as a merging system where each sub-cluster is modeled individually according to the description given in this section. We then co-add the two resulting templates to form a composite spatial map which we use in the likelihood analysis. The projected annihilation flux maps are shown for the different models we choose in Figure 8. In this figure we added a contour that indicates the spatial position of the excess we report in this paper and we stress that a DM origin is unlikely because of the large offset between the predicted DM annihilation profile and the contours of the γ-ray excess (see Section 4 for details).

Zoom In Zoom Out Reset image size

Table 3. Virgo Subclusters and Derived DM Density Profiles

Sub-cluster	M200	r200	θ200	c a	JNFW	JDM-I a	JDM-II b
	(×1014 M⊙)	(Mpc)	(°)		(×1017)	(×1018)	(×1020)
M87 (Virgo-I)	5.60	1.70	6.3	4.21	2.56	6.50	3.33
M49 (Virgo-II)	0.72 c	0.88	3.8	5.58	1.85	5.36	0.75

Notes. Shown are the characteristic quantities used to derive the resulting J-factors for the Virgo cluster modeled as a merging system between the sub-clusters associated with M87 and M49. Columns from left to right are name, mass, virial radius, angular radius θ₂₀₀, concentration parameter c, as well as J-factors for NFW and the DM models used in this analysis for each of the sub-clusters. All J-factors are given in units of GeV² cm⁻⁵ and have been computed over a solid angle subtending the virial radius of each sub-cluster.

^a Sánchez-Conde & Prada (2014). ^b Gao et al. (2012). ^c Chen et al. (2007).

Download table as:  ASCIITypeset image

We summarize the main characteristics of our chosen models in Table 3. We take the sum of the J-factors for Virgo-I and Virgo-II to be the total J-factor of the Virgo system.

7.2. DM Flux Enhancement due to IC Scattering

7.3. Limits on $\langle \sigma v\rangle$

We derive upper limits on $\langle \sigma v\rangle$ using the profile likelihood method (Rolke et al. 2005) as implemented in the MINOS-subroutine of the MINUIT package (James & Roos 1975) which is available through the Fermi Science Tools. We define the 95% upper limit on $\langle \sigma v\rangle$ as the value of $\langle \sigma v\rangle$ for a given mass m_χ where twice the difference in the log-likelihood, $2\times {\rm{\Delta }}{\mathcal{L}}=2.71$ with respect to the value of the log-likelihood for the best fit value. ⁷³

Figure 9 shows the dependency of the upper limits on the chosen DM annihilation channel for our fiducial models. The most constrained channels are $\chi \chi \to b\bar{b}$ and $\chi \chi \to {\tau }^{+}{\tau }^{-}.$ Accounting for IC emission in the leptonic channels e^± and μ^± improves the constraints we obtain from the prompt emission by two to three orders of magnitude, above DM masses of 50 GeV and 100 GeV, respectively. The limits for e^± are the most constraining for DM masses above ∼110 GeV, due to the enhanced flux predictions from IC.

Zoom In Zoom Out Reset image size

In Figure 10 we show our derived upper limits on $\langle \sigma v\rangle$ and their associated TS values for the $\chi \chi \to b\bar{b}$ channel and contrast our standard IEM with results obtained from using the alternative diffuse models as discussed in Section 5. Our optimistic limits exclude thermal WIMP cross sections below 40 GeV. The limits derived from the more conservative assumptions are a factor 20 weaker across the entire probed mass range. Even with the inclusion of additional point sources as done in this work, there is a residual TS ∼ 4 if we consider the more extended and elongated profile as predicted by our optimistic model (DM-II). For the DM-I case, this value is reduced even further. Considering the alternative diffuse models, the resulting limits are generally weaker, associated with residual TS < 5 except for WIMP masses ≲20 GeV. For lower masses, the alternative models give rise to residual TS peaking at TS ∼ 9.5 or ∼3.1σ. Low-mass DM models associated with relatively high TS-values for one diffuse model show a large spread (δTS ≃ 5) in TS for the alternative models.

Zoom In Zoom Out Reset image size

Recalling the description of the alternative diffuse models in Section 5, these differ from the standard IEM by having various large-scale components fit freely to the data (e.g., Loop I, IC, etc.). The extent of these large-scale components is comparable to the spatial extension of our cluster template which causes a degeneracy between the fit parameters for the diffuse components and Virgo. As a consequence we find that soft photons (E ≲ 10 GeV), which would otherwise be attributed to the background IEM, are now included in the number of predicted photon counts from Virgo for a light WIMP model. ⁷⁴ Note that this effect appears to be even more pronounced as the spatial template for the Virgo cluster is even more extended than the disk used in our previous study (refer to Section 4 for a detailed discussion). Finally, we also remark that this issue is by construction less apparent for the standard IEM, since here all components are fixed to their relative best-fit contributions obtained from a likelihood fit to the entire γ-ray-sky.

An alternative production mechanism of γ-rays originating from the Virgo region may be due to CR interactions. γ-rays are mainly produced in IC interactions of relativistic electrons or via hadronic pp-collisions producing pions and γ-rays through ${\pi }^{0}\to 2\gamma$ (Brunetti et al. 2012). The dominant production mechanism of γ-rays from CRs in the ICM is still debated: either cosmic-rays are accelerated directly in structure formation shocks (including the effect of AGNs and supernovae) through diffusive shock acceleration or an aged population of cosmic-ray are reaccelerated in the turbulent plasma generated by, e.g., merging clusters (see, e.g., Brunetti & Jones 2014, for a review).

Since there is no giant radio halo associated with the Virgo cluster and the central part of the cluster has properties similar to a cool-core cluster (Urban et al. 2011), we expect the γ-rays from a population of reaccelerated cosmic-rays (see, e.g., ZuHone et al. 2013) to be too faint to be detectable by the Fermi-LAT throughout its lifetime. However, there is a strong dependence on the uncertain turbulent profile. Indeed, Pinzke et al. (2015) showed that for a flatter turbulent profile than what was previously assumed, the γ-ray emission could be in reach with Fermi-LAT in the coming years. To keep the CR analysis simple, we neglect these aforementioned models as well as other leptonic models (Kushnir & Waxman 2009). Instead, we focus on constraining the γ-rays produced in a pure hadronic scenario in that region. Specifically, we adopt a simple but realistic model for the predicted universality of the CR spectra built up from diffusive shock acceleration in large-scale structure formation shocks (Pfrommer 2008; Pinzke & Pfrommer 2010). Based on these considerations, in this section we derive constraints on the CR-induced γ-ray flux and related CR quantities from Virgo.

8.1. Modeling and Results

Following earlier works in Ackermann et al. (2014a), we consider two different (hadronic) models for the CR distribution, the simulation-based approach by Pinzke et al. (2011), which predicts a γ-ray surface brightness which closely follows the X-ray emitting gas in the ICM, and a model in which the CRs are confined within the cluster virial radius but evenly distributed with no dependence on the ICM gas (flat model). The latter can thus be seen as a simplified proxy for CR-streaming models which can lead to more extended γ-ray brightness profiles (Enßlin et al. 2011; Wiener et al. 2013; Zandanel & Ando 2014). While the expected γ-ray morphology varies, we assume the spectrum to be approximated by the universal model as detailed in Pinzke et al. (2011) (the interested reader is referred to Figure 1 of that paper). Analogously with the results presented in Section 7.1, we construct a model that takes into account the merging state of the cluster by overlaying the spatial template inferred from X-ray profiles from Virgo-I with that of Virgo-II. We show the predicted flux maps in Figure 11. Outside r₂₀₀ we take the predicted flux to be negligible.

Zoom In Zoom Out Reset image size

In analogy with the results in the previous section, we use the profile likelihood method to derive 95% upper limits on the CR-induced γ-ray flux. Our results are shown in Table 4.

Table 4. CR-models and Derived Limits

CR Model	Fγ,pred(E > 100 MeV)	Fγ,95(E > 100 MeV)	$\langle {X}_{\mathrm{CR}}\rangle$	${\zeta }_{p,\mathrm{max}}$
	(×10−9 ph cm−2 s−1)	(×10−8 ph cm−2 s−1)
Simulation-based	15.0	1.2	6%	40%
Flat CR a	0.4	1.8	...	...

Notes. Shown are both predicted and observed integrated fluxes above 100 MeV for our models for CR-induced γ-rays as discussed in the text. For the simulation-based model, the remaining columns denote the volume-averaged CR-to-thermal pressure ratio and the maximum acceleration efficiency for CR protons, respectively. ^b

^a In order to provide a consistent description, we normalize each profile to the total CR number within r₂₀₀ for the simulation-based model. ^b As the observed flux for the flat CR-model is a factor ∼45 above the predictions, these limits cannot be used to constrain $\langle {X}_{\mathrm{CR}}\rangle$.

Download table as:  ASCIITypeset image

We exclude γ-ray integral fluxes above 1.2 × 10⁻⁸ ph cm⁻² s⁻¹ for the simulation-based CR model over the energy range 100 MeV–100 GeV, which is about a factor ∼1.4 stronger than previously published (Ackermann et al. 2010a). Using the flat model yields an integral flux limit of 1.8 × 10⁻⁸ ph cm⁻² s⁻¹ which is above the value that was published previously. This can however be explained by the fact that flat CR models are generically less constrained by current γ-ray data (Ackermann et al. 2014a; Zandanel & Ando 2014).

8.2. Constraints on ${\zeta }_{p,\mathrm{max}}$ and $\langle {X}_{\mathrm{CR}}\rangle$

8.3. Systematic Uncertainties due to IEM Modeling

8.4. Degeneracy of Results with M87

We find no strong evidence for extended emission associated with the Virgo cluster center. Yet, using the standard IEM we find a statistically significant extended excess from a disk profile with radius 3° clearly offset from the cluster center. Our TS map reveals two well-separated maxima, both clearly offset from the two main sub-clusters associated with the giant ellipticals M87 and M49. This signature makes a DM origin unlikely. Also, as there is no indication of accelerated CRs, evidenced by either radio or X-ray emission, an astrophysical origin due to, e.g., accelerated CRs in the ICM is questionable. We thus report upper limits on CR scenarios and DM-induced γ-rays.

Similar to previous studies, we carry out a search for new point sources in order to account for the increased data volume with respect to the employed source catalog. We find six new candidates in accordance with similar studies by Macías-Ramírez et al. (2012) and Han et al. (2012a). These new candidates, however, have no reported counterparts in other wavebands. Five of them are contained in the 3FGL-catalog. We carry out an alternative IEM study which is essential for estimating systematic uncertainties associated with the search for γ-ray emission from very extended sources. In our case the inconsistency between the IEMs demonstrates that the Virgo region is an especially difficult section of the sky. The proximity to poorly understood Galactic foregrounds emitting γ-rays, like Loop I, makes the search for extended emission from this region very challenging. Our study also reveals the challenges of searches for such low photon density sources even at high Galactic latitudes.

Accounting for the complex dynamics of the cluster, we model its emission by co-adding the contributions from the major sub-clusters centered on M87 and M49, respectively. In particular for very extended models, as predicted if considering large amounts of DM halo substructure, the spatial morphology departs from spherical symmetry. Resulting limits for either DM- or CR-induced γ-rays are generally weaker than that of other targets, e.g., dwarf spheroidal for the case of DM annihilation (Ackermann et al. 2011) and from collective cluster studies (Ackermann et al. 2014a). The DM limits from the Virgo analysis here, for instance, are about an order of magnitude above the thermal WIMP cross-section when assuming a realistic model for the sub-halo boost.

Finally we would like to stress that the main findings in this paper are expected to remain unchanged even if more data were to be included, as the uncertainties in the results are dominated by systematics associated with the IEM modeling. We emphasize that the improved source model used with the analysis roughly corresponds to the model presented in the current deepest γ-ray catalog, 3FGL. Also, while the predicted constraints on DM annihilation and CR processes improve by up to a factor of ≃1.4 if all available data are considered (six years instead of three years), targets other than the Virgo cluster may be better suited for analysis, e.g., the Coma cluster for CR processes and Fornax for DM prospects (see, e.g., Pinzke et al. 2011; Ando & Nagai 2012, for a discussion). ⁷⁵ Farther away, the predicted γ-ray emission from both clusters is expected to be within the detection reach of the LAT and their apparent extensions on the sky is significantly less than Virgo, which helps reduce the uncertainties associated with the foreground IEM modeling, thus allowing for a more robust analysis.

A.P. is grateful to the Swedish Research Council for financial support.

The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden.

Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France.

Facility: Fermi - Fermi Gamma-Ray Space Telescope (formerly GLAST)