SEARCH FOR GAMMA-RAY EMISSION FROM THE COMA CLUSTER WITH SIX YEARS OF FERMI-LAT DATA

The radiative yields of energetic cosmic rays (CRs) traversing intracluster gas in galaxy clusters and interacting with background radiation fields span a wide spectral range from radio to high-energy γ rays. Extended regions of radio emission (referred to as halos and relics) have already been observed in many clusters (e.g., Ferrari et al. 2008). Compton scattering of relativistic electrons by the cosmic microwave background (CMB) radiation is the dominant process for emission above 50 keV (where thermal emission from hot intracluster gas is sufficiently weak) up to O(100) $\mathrm{MeV}$ (e.g., Brunetti & Jones 2014). Searches for this non-thermal (NT) X-ray emission have not yielded conclusive results (e.g., Rephaeli et al. 2008; Ajello et al. 2009; Wik et al. 2009).⁵⁸ At energies higher than O(100) $\mathrm{MeV}$ the radiative decay of neutral pions (produced in energetic proton interactions with ambient protons, henceforth referred to as p–p) is expected to dominate, if CR interactions in the intracluster gas take place at non-negligible rates.

While there is not yet observational evidence for energetic protons in clusters, the presence of relativistic electrons is well established from radio observations. Clearly, strong radio galaxies in clusters are sources of relativistic electrons: e.g., M87 in Virgo (Bolton et al. 1949), NGC 4869 in Coma (Willson 1970), NGC 1275 in Perseus (van den Bergh 1961); see also Dutson et al. (2013) for a comprehensive search for the brightest cluster galaxies. It is also possible that energetic particles are (re)accelerated in the intracluster medium (ICM), by, e.g., merger and accretion shocks (for a recent review, see Brunetti & Jones 2014). From the large extent of cluster radio halos and typical synchrotron energy loss times, it is commonly thought that the emitting electrons are secondary, namely decay products of charged pions (produced in p–p interactions; see, e.g., Brunetti & Jones 2014, for a recent review). Since contemporary viable models predict that the dominant γ-ray emission process is π⁰ decay (e.g., Brunetti & Jones 2014) and since the electron lifetimes to radiative losses are generally significantly shorter than the source crossing time (Petrosian 2001), in this work we do not consider purely leptonic models.

An analysis of 50 clusters using four years of data from the Large Area Telescope (LAT) on board the Fermi satellite (Ackermann et al. 2014a) resulted in upper limits on the CR-induced γ-ray emission (see also Ackermann et al. 2010; Huber et al. 2013; Griffin et al. 2014; Prokhorov & Churazov 2014). The γ-ray upper limits from these analyses put stringent constraints on the energy density of energetic protons in clusters, severely constraining models for acceleration by structure formation shocks (reviewed by Brunetti & Jones 2014).

Due to substantial interest in exploring cluster NT emission, and the availability of a longer LAT data set with an improved event-level analysis (see below), a dedicated analysis of the most favorable cluster candidates is warranted. The Coma cluster is a natural choice for our search, being the nearest massive cluster with a bright radio halo. Coma is also located near the North Galactic pole where the diffuse gamma-ray intensity is at a minimum. The vast array of broadband observations, including detailed measurements of the radio halo in Coma, provide a sound basis for testing theoretical models for the gamma-ray emission (e.g., Gabici & Blasi 2004; Reimer et al. 2004; Berrington & Dermer 2005; Pinzke & Pfrommer 2010; Brunetti et al. 2012; Pinzke et al. 2015). Coma has been studied using the LAT (e.g., Ackermann et al. 2010; Arlen et al. 2012; Han et al. 2012; Prokhorov 2014; Zandanel & Ando 2014), its predecessor, EGRET (e.g., Reimer et al. 2003), as well as Cherenkov telescopes (e.g., Aharonian et al. 2009; Arlen et al. 2012), but has yet to be detected in MeV-to-TeV γ rays. We report here on the deepest observation of the Coma cluster covering the $\mathrm{MeV}$–$\mathrm{GeV}$ band to date, obtained using six years of LAT data analyzed with Pass 8 (Atwood et al. 2013).

We describe the data analysis procedure in Section 2 and present our results in Section 3. A comparison of the results with predicted emission levels is presented in Section 4; we briefly summarize our conclusions in Section 5. In this paper we use the Planck measurement of the Hubble constant, H₀ = 70 km s⁻¹ Mpc⁻¹ (Planck Collaboration et al. 2014). At a distance of ≃100 Mpc (z = 0.023), the virial radius of the cluster of 2.0 Mpc corresponds to a subtended angle on the sky, θ₂₀₀ = 123 (Ackermann et al. 2014a).

The LAT is a pair conversion telescope with a large field of view of ∼2.4 sr sensitive to γ rays from ∼20 $\mathrm{MeV}$ up to >300 $\mathrm{GeV}$ (Atwood et al. 2009; Ackermann et al. 2012c).⁵⁹ This work uses 6 years (MET 239557414–428903014) of public Pass 8 LAT data. Pass 8 is an extensive rewrite of the core reconstruction algorithms (Atwood et al. 2013) applied to the data taken by the LAT that is informed by in-flight performance and characteristics of the LAT, which results in an increase in γ-ray acceptance by 20%–40% with respect to source events reconstructed with Pass 7REP (Bregeon et al. 2013), depending on energy (see Appendix A for details). We select photons with energies from 100 $\mathrm{MeV}$ to 10 $\mathrm{GeV}$ within a region of interest (ROI) centered on the Coma cluster at α₂₀₀₀ = 194.95, δ₂₀₀₀ = 27.98 and select a square of 15° × 15° as our ROI.⁶⁰

Limiting the data selection to zenith angles less than 90° allows us to effectively remove photons originating from the Earth limb.⁶¹ We use gtmktime to select time intervals during which the LAT was in nominal science operations mode and excluded intervals coincident with γ-ray bursts and solar flares. We bin our data in 16 logarithmically spaced bins in energy and use a spatial binning of 01 per pixel.

In our analysis we include all sources listed in the latest Fermi catalog of point-like and extended sources (3FGL, Acero et al. 2015) along with the standard diffuse Galactic foreground emission model recommended for point source analysis. We include a spatially isotropic model accounting for the extragalactic diffuse γ-ray background and misclassified CRs. The normalizations of both Galactic and extragalactic diffuse emission models are left free in the likelihood fit. The model components are convolved with the parametrized detector response represented by the P8R2_SOURCE_V6 instrument response functions (IRFs).

We make use of recent developments of the Science Tools package that incorporate the finite energy dispersion of the LAT for the likelihood analysis by convolving the point sources in our ROI with the LAT energy dispersion.⁶² Note that the energy dispersion is already accounted for when creating the above-mentioned diffuse models. During the likelihood fit, we allow all sources that are separated from the cluster center by less than 65 to have a free normalization to allow for the broad PSF at the lowest energies (the 68% containment radius of photons at normal incidence with an energy of 100 $\mathrm{MeV}$ is roughly 4°). This choice ensures that 99.9% of the predicted γ-ray counts (integrated over energy for an E^−2.3 spectrum, which would roughly correspond to the predicted cluster emission) is contained within the chosen radius (that is substantially larger than θ₂₀₀ of the cluster, see Section 2.1 for details).

In order to assess the completeness of our reference emission model for the Coma region (without a cluster source), we use the gttsmap tool to search for any additional gamma-ray sources (Figure 1). For each pixel in the map we evaluate the test statistic (TS), defined as twice the likelihood ratio between the best-fit value of the alternative hypothesis, including a source and the best-fit value of the null-hypothesis (background only). In the TS-map a point source modeled as a power law with Γ = 2.3 is tested at each pixel of the map. We find four new point source candidates that coincide with regions with individual TS values >25. The details of the improved background model along with the best-fit values of all free parameters are given in Appendix B.

Zoom In Zoom Out Reset image size

We find two residual structures whose peak TS values in individual pixels (with a pixel size of 005) are ∼9 and ∼13, respectively. These structures, located within the virial radius, are separated from each other by about 1 Mpc at the distance of the Coma cluster, if within the cluster. We further investigate this potential excess emission by repeating the TS-map calculation in a soft (E < 1 $\mathrm{GeV}$) and a hard (E > 1 $\mathrm{GeV}$) energy band. While the hard band reveals at least three distinct areas of excess emission that spatially overlap with the coordinates of X Comae, and NGC 4839, respectively, the soft map does not permit an immediate association with discrete sources.⁶³ Moreover, when comparing our results with radio observations of the Coma cluster we find that the weak diffuse excess in the soft band roughly overlaps with that of the radio halo. In Section 4 we compare the level of this residual emission with predicted emission from likely sources in the cluster. We further investigate this apparent excess by creating spectral residuals in each energy bin i; for this we take the number of observed photons in each bin c_i and then evaluate the number of predicted photon counts given our background-only hypothesis m_i and calculate the spectral residual ${r}_{i}=({c}_{i}-{m}_{i})/{m}_{i}$. We evaluate the residuals over the full ROI as well as within 0.5 and 1.0 × θ₂₀₀ and show these results in Figure 2. While the residuals evaluated over the entire ROI are practically negligible, we note an indication of a slight excess within the cluster virial radius.

Zoom In Zoom Out Reset image size

2.1. Cluster Spatial Modeling

In this section we motivate our benchmark spatial template, hereafter referred to as the cored profile, which is motivated by observations of the radio halo in Coma.

Evidence for NT particle populations in clusters is currently limited to measurements of extended regions of radio emission. The diffuse halo component in clusters (mostly intracluster emission) that remains after subtracting compact sources, has a power law spectrum with spectral flux ${f}_{\nu }\propto {\nu }^{-\alpha }$, with α ∼ 1.0–1.5 (Ferrari et al. 2008), considerably steeper than typical Galactic spectra (α ∼ 0.60–0.75). The radio halo in Coma, which is centered on two dominant radio galaxies (NGC 4869 and NGC 4874), has α ≃ 1.3 (Schlickeiser et al. 1987; Thierbach et al. 2003). From the recent 352 MHz WSRT measurements of Brown & Rudnick (2011) the halo is ∼1 Mpc in size. The contour map in the latter paper, or, more conveniently, its conversion to an azimuthally averaged profile (as plotted in Figure 10 of Planck Collaboration et al. 2013) can be parametrized by the following distribution in terms of the angular distance θ from the cluster center

$$\begin{eqnarray}&&{f}_{r}(\theta )={\left[1+\displaystyle \frac{{\theta }^{2}}{{\theta }_{c}^{2}}\right]}^{-p},\end{eqnarray} \tag{ 1 }$$

with θ_c ≃ 15' and p ≃ 1.4.

The synchrotron spectral flux is proportional to the projected (along the line of sight) value of the product N_eB^α+1, where N_e and B are the emitting electron (column) density and locally averaged value of the magnetic field, respectively; α is the spectral index of the radio-synchrotron emission. Due to heavy energy losses of primary electrons (with the required energy of more than a few GeV) propagating out of their source regions, it is usually assumed that the dominant contributions to the (observed) radio and (predicted) γ-ray emissions are from charged (π^±) and neutral (π⁰) pion decays, respectively. The spectral intensity of the γ-ray emission is predicted to have the characteristic π⁰ decay bump shape, unless the emission is dominated by a Compton component due to scattering of the secondary electrons off the CMB, in which case a power law with an index nearly equal to α would be predicted.

Whereas the spatial profile of N_e is essentially the same as that of both π^± and π⁰, respective sources of the radio and γ-ray emission, the radio profile reflects also the spatial dependence of B^α+1. In the hot, fully ionized intracluster gas, electrical conductivity is high, and magnetic fields are expected to be frozen into the plasma. If so, the field spatial dependence is expected to scale as n^2/3, with n being the intracluster gas density (Rephaeli 1988); the index in this scaling is only somewhat lower, 1/2, if magnetic energy (rather than flux) is conserved. With a commonly used gas density profile of the form ${(1+{r}^{2}/{r}_{c}^{2})}^{-3\beta /2}$, where r_c ($\propto {\theta }_{c}$) is the gas core radius, and typically β = 2/3, it follows that the (projected) profile of the γ-ray emission is expected to be considerably shallower than the measured radio profile.⁶⁴

Taking the radio profile from Equation (1) and β = 2/3, the profile of the γ-ray flux is predicted to be at most moderately steep, ${f}_{\gamma }{(\theta )\propto (1+{\theta }^{2}/{\theta }_{c}^{2})}^{-1/4}$.

In addition to the cored profile, we consider models that bracket the two extremes of possible spatial models: one that is based on the cluster being modeled as a point source (point-like) and another in which we assume a uniform distribution of predicted photons out to the virial radius (disk-like). Note that the the tabulated data for the bin-by-bin likelihoods for disk emission of varying size (0.1 × R₂₀₀ − 1.0 × R₂₀₀) is available in a supplementary tar package.

2.2. Likelihood Analysis

The resultant SED for the Coma cluster is shown in Figure 3. The SED was deduced by first deriving the observed 95% confidence level one-sided upper limits, selected as the interval within which twice the difference in the log-likelihood with respect to the best-fit value of the alternative hypothesis (including the cluster) equals 2.71. These limits were then compared with those obtained by randomly selecting high Galactic latitude regions of the sky and repeating the analysis in these blank fields.⁶⁵ The observed limits are typically within the 68% containment band, while a few bins fall in the 95% containment band.

Zoom In Zoom Out Reset image size

The information from each bin is combined into the global log-likelihood that contains the spectral information over the whole energy range of the analysis:

$$\begin{eqnarray}&&{ \mathcal L }(\mu ,\hat{\theta }| D)=\displaystyle \prod _{i}{{ \mathcal L }}_{i}({\mu }_{i},\hat{\theta }| {D}_{i}).\end{eqnarray} \tag{ 2 }$$

In this equation we have the binned Poisson likelihood term ${{ \mathcal L }}_{i}$, in each energy bin i, where D_i corresponds to the observed counts in bin i and μ_i is the normalization (the only free parameter) of the power law with Γ = 2.3 that we have assumed in each bin. $\hat{\theta }$ indicates that the nuisance parameters, θ (i.e., normalizations of free point sources and diffuse background emission components) are taken from the global fit over the entire energy range. Irrespective of the origin of the dominant emission process (be it hadronic or leptonic) in the relevant energy range explored here, the spectrum can effectively be assumed to have a power-law form. Thus, we select a set of spectral indices for which we provide integral flux upper limits between 100 $\mathrm{MeV}$ and 10 $\mathrm{GeV}$(bottom row of Figure 3). The maximum TS value found in our analysis is ∼13 for the cored profile (Γ = 2.3). When using the blank fields to assess the null-hypothesis distribution, and considering the trials factor associated with testing both for disk and point-like emission, this TS value corresponds to a global significance of ∼1.8σ. Given the low statistical significance of the residual emission, we determine integral flux upper limits.

The deduced integral flux limits in the energy range from 100 $\mathrm{MeV}$ to 100 $\mathrm{GeV}$ are 5.2 × 10⁻⁹ ph cm⁻² s⁻¹ for our benchmark model. Modeling the cluster as a point source yields a limit that is almost a factor of two lower, 3.1 × 10⁻⁹ ph cm⁻² s⁻¹, whereas for the uniform disk model we deduce a slightly higher limit 5.8 × 10⁻⁹ ph cm⁻² s⁻¹. Note, however, that the observed limit in each case is above the blank-field median expectation, falling instead into the 68% containment band for the disk and the 95% containment band for point-like and cored profiles, respectively. Our main results are summarized in Table 1.

Table 1.  Summary of the Likelihood Analysis Results

Spatial Model	Power-law Index	${F}_{\mathrm{obs}}^{95\%}(E\gt 100\;\mathrm{MeV})$	$\bar{{F}_{\mathrm{exp}}^{95\%}}(E\gt 100\;\mathrm{MeV})$
		(×10−9 ph cm−2 s−1)	(×10−9 ph cm−2 s−1)
Cored Profile	2.3	5.2	1.7
Point Source	2.3	3.1	1.1
Disk	2.3	5.8	2.1
Cored Profile	2.0	3.2	1.1
Point Source	2.0	1.7	0.5
Disk	2.0	3.8	1.4

Note. From left to right: spatial model assumed along with the adopted power-law index followed by the observed integral flux upper limit. The last column denotes the median expected integral flux upper limit that we determined from blank fields. All limits are calculated at 95% C.L. Note that while we refer to Γ = 2.3 as our benchmark model, we also list the values for Γ = 2.0 for easier comparison with previous works (see the text for details).

Download table as:  ASCIITypeset image

The limits we present here appear to be weaker than previously reported ones, most notably when comparing our results with those from the most recent analysis of the Coma cluster using Fermi-LAT data (Zandanel & Ando 2014). However, the observed limits (solid black line in Figure 3) cannot be compared naively this way because our analysis uses a larger (and different) data set and an improved background model in terms of sources included in the ROI (see Appendix A for details regarding a comparison of five years of Pass 7 and Pass 8 data). Since Zandanel & Ando (2014) report limits with associated TS values of zero, these results should instead be compared to our expected sensitivity based on the blank field study. Doing so we find that the expected sensitivity for emission from a point-like source is similar to that reported in earlier works by Zandanel & Ando (2014), but the sensitivity for an extended source has improved significantly; the median expected limit for the disk-like emission for Γ = 2.0 is a factor of two lower than what was reported in earlier works.

3.1. Systematic Uncertainties.

To assess the robustness of our results we perform a number of systematic checks and vary fiducial parameters such as the spectral or spatial binning as well as the number of free sources and repeat our analysis for these choices. In particular we investigate how these changes affect the bin-by-bin likelihood evaluations and how these changes influence our likelihood analysis for soft and hard spectra. We tested our results against a set of IRFs that represent minimal or maximal boundaries in the computation of the effective area and PSF within the systematic uncertainties of Pass 8 that are chosen to maximize and minimize effective area and PSF, respectively. To explore the uncertainties associated with the diffuse foreground emission, we follow the method in Ackermann et al. (2014a) by employing a set of alternative models (Ackermann et al. 2012b; de Palma et al. 2013). The total systematic uncertainties of the analysis we present here typically do not exceed 22% for E ≳ 300 $\mathrm{MeV}$. The fractional uncertainty increases up to 54% if the full energy range (100 $\mathrm{MeV}$ ≲ E ≲ 10 $\mathrm{GeV}$) is considered. For extended sources with fluxes near the current sensitivity limit of the LAT, the uncertainty from diffuse foreground modeling is at least as large as the uncertainties from the other contributors, as shown in Table 2.

Table 2.  Budget of Systematic Uncertainties

Category	Variation	SED Impact	Integral (Hard)	Integral (Soft)
IRFs	bracketing MIN/MAXa	<10% (E ≤ 300 $\mathrm{MeV}$)	<3%	<5%
	(σAeff ∼ 3%, σPSF ∼ 5%)	<1% (300 $\mathrm{MeV}$ ≤ E ≤ 2 $\mathrm{GeV}$)	⋯	⋯
	⋯	<5% (E ≥ 2 $\mathrm{GeV}$)	⋯	⋯
Diffuse modeling	alt. diffuse modelsb	<50% (E ≤ 300 $\mathrm{MeV}$)	<12%	<28%
	⋯	<10% (E > 300 $\mathrm{MeV}$)	⋯	⋯
Free radiusc	10° (nominal 6.5)	<10%	<4%	<6%
Spectral binning	±50%	⋯	<10%	<10%
Spatial binning	±50%	<10%	<2%	<2%
SED index	±10%	<2% (E > 400 $\mathrm{MeV}$)	<5%	<5%
		<10% (E < 400 $\mathrm{MeV}$)	⋯	⋯
Totald	⋯	<54% (E ≤ 300 $\mathrm{MeV}$)	<21%	<42%
		<22% (E > 300 $\mathrm{MeV}$)	⋯	⋯

Notes. For the evaluation of systematic uncertainties we consider two cases: the effect on the bin-by-bin likelihood (here the uncertainties should be applied to the upper limit reported in each bin) and the effect on the broadband spectra assuming a hard (Γ = 1.6) and soft (Γ = 2.6) power-law spectrum when all spectral bins are combined (here the uncertainty is given with respect to the integral flux upper limit).

^aSee Sections 5.7 and 6.5 in Ackermann et al. (2012a) for a detailed discussion of the bracketing IRF approach. ^bSee Ackermann et al. (2012b) for the detailed model description. ^cCompared to the nominal choice of 65, there are 13 more sources that are left to vary in the fit. ^dValues have been added in quadrature.

Download table as:  ASCIITypeset image

Our analysis of six years of LAT Pass 8 data above 100 $\mathrm{MeV}$ for the Coma cluster region does not reveal excess emission at a significance level that is large enough to claim detection of diffuse emission in the ICM. The TS map in Figure 1 is nonetheless of appreciable interest: there seems to be residual (background-subtracted) emission from an area that overlaps partly with the Coma virial radius. It is therefore reasonable to compare the sensitivity of our analysis (as gauged through the blank field analysis) with the emission predicted from likely sources in the cluster.

The deduced median integral flux upper limits for our benchmark and disk models with Γ = 2.3, which corresponds to the measured index of the radio (spectral energy) flux specified above (α = 1.3), translate to the luminosity bounds of L(>100 $\mathrm{MeV}$) ∼ 1.4 × 10⁴² erg s⁻¹ and L(>100 $\mathrm{MeV}$) ∼ 1.8 ×  10⁴² erg s⁻¹, respectively, on the combined emission from all compact and extended sources in the Coma region.

The combined 408 MHz flux density from the two central radio galaxies NGC 4869 (5C 4.81) and NGC 4874 (5C 4.85) is 2.1 × 10⁻²³ erg cm⁻² s⁻¹ Hz⁻¹ (Willson 1970; Kim et al. 1994). An estimate of the combined γ-ray luminosity of these two galaxies can be readily made if it is assumed that the radio continuum emission is mostly due to Compton scattering of the radio-producing electrons by the CMB, a reasonable expectation in light of the conclusion reached by Abdo et al. (2010) in their analysis of Fermi-LAT measurements of Cen A, the nearby radio galaxy. If so, the total Compton luminosity, L_C, can be estimated from the measured radio (bolometric synchrotron) luminosity, L_S, by the simple scaling L_C ≃ L_Sρ₀/ρ_B, where ρ₀ and ρ_B = B²/8π are the (present) CMB and magnetic field energy densities, respectively. With a mean field strength of ∼5 μG across the the extended radio regions of NGC 4874 and NGC 4869 (Feretti & Giovannini 1987; Feretti et al. 1990), we compute L_C ∼ 6 × 10⁴⁰ erg s⁻¹, and L(0.1–10 $\mathrm{GeV}$) ∼ 2 × 10⁴⁰ erg s⁻¹. Additional emission from other radio galaxies in the cluster renders this estimate a lower limit on the total γ-ray emission from all radio sources in Coma.

The other discrete ("compact") sources are mostly from star-forming galaxies, whose contribution to the cluster γ-ray emission was estimated by Storm et al. (2012) based on a scaling relation between galactic IR and γ-ray emission. They estimated that the superposed 0.1–100 $\mathrm{GeV}$ luminosity of star-forming galaxies in Coma is in the wide range of ∼3 × 10⁴⁰–3 × 10⁴² erg s⁻¹. Together, the estimates of radio and star-forming galaxies yield a lower limit on the combined emission from galaxies in Coma, L_gal(>100 $\mathrm{MeV}$)  ∼ 5 × 10⁴⁰ erg s⁻¹. (We note in passing that the predicted γ-ray emission from X Comae, a background active galactic nucleus at z = 0.091 located about a degree north to the Coma center, is too weak to be of relevance for our discussion here, since its radio flux is much weaker than the two radio galaxies within the cluster.)

Quantitative estimates of the intracluster γ-ray emission can only be made in the context of specific models for intracluster energetic protons and electrons producing this emission in interactions with intracluster gas (via π⁰ decay) and Compton scattering off the CMB, respectively. A viable model must include the particle source distribution, the propagation mode, and the spatial distributions of intracluster gas and magnetic field. While models for cluster NT particles abound (reviewed recently by Brunetti & Jones 2014), and their predicted levels of γ-ray emission span a wide range, we can roughly estimate a minimal level of γ-ray emission from the measured radio halo luminosity, avoiding the need for a detailed model (and inherent untested assumptions). The halo flux is about twice (at 408 MHz) that of the two central radio galaxies, and the mean (volume-averaged) intracluster magnetic field is at least a factor five lower than the value adopted for the central radio galaxies (Bonafede et al. 2010). We would then expect that the total number of radio-emitting electrons in the halo to be considerably higher than that in the central radio galaxies. In the context of a secondary electron origin of the measured radio halo emission, it can be shown that the extended π⁰-decay γ-ray emission is comparable to or higher than that from the central radio galaxies. Indeed, this has recently been quantified in a detailed study by Rephaeli & Sadeh (2016).

The above considerations lead to the lower limit L(>100 $\mathrm{MeV}$) ∼ 1 × 10⁴¹ erg s⁻¹ on the combined galactic and extended γ-ray emission from Coma. Based on this estimate the predicted minimal γ-ray flux is at a level that is an order of magnitude lower than our deduced upper limit. Aside from the instrument response functions, the measurement sensitivity is mainly governed by two factors: the detected (and thus modeled) discrete sources observed by Fermi, i.e., the 3FGL, and the model for the Galactic and isotropic γ-ray foreground emission. Moreover, the residuals found in this work may be the first indication of cluster emission which could become more significant with further observations.

The extrapolated sensitivity (given by the median expectation from blank fields, scaled to a longer exposure time) of our analysis to 10 years of LAT exposure is comparable to recent predictions of diffusion models of CR protons with turbulent reacceleration of leptonic secondaries (Brunetti & Lazarian 2011; Brunetti et al. 2012). However, even if the exposure of the object could be doubled by the end of the Fermi mission, at least in the context of the analysis techniques presented here, it appears unlikely that Coma will be significantly detected in the LAT data given the current results. This forecast, based upon the enhanced sensitivity of Pass 8 and a comparison to blank fields, is in agreement with the conclusions of Zandanel & Ando (2014), who found no indications of residual emission (TS ∼ 0) in an analysis of 5 years of Pass 7 reprocessed LAT data.

The main conclusions from this paper can be summarized as follows.

• 1.  
  In an analysis of γ rays between 100 $\mathrm{MeV}$ and 10 $\mathrm{GeV}$ collected over a period of six years with the LAT and reprocessed with Pass 8, we find excess emission within the cluster virial radius; however, the statistical significance of this emission is well below the threshold to claim detection of γ-ray emission from the cluster.When spectrally analyzed, this excess emission above the background expectation can be separated into a soft and a hard component; the soft (E < 1 $\mathrm{GeV}$) appears to roughly spatially overlap with parts of the well measured giant radio halo, while the hard component (E > 1 $\mathrm{GeV}$) can be associated with parts of the halo.
• 2.  
  Using a γ-ray template derived from the combination of Planck and WSRT observations of the radio halo, we find a maximum TS value of ∼13 for power-law emission with a spectral index of ∼2.3 (after correcting for trial factors, this corresponds to a global significance of ∼ 1.8 σ). We derive limits on the integral CR-induced γ-ray flux, and our observed limit excludes fluxes above 5.2 × 10⁻⁹ ph cm⁻² s⁻¹. While we focus in this paper on hadronic models where the radio emission is of secondary origin, we emphasize that the results presented here are of a more universal nature; the relatively small variation of ∼60% between our benchmark model and the most extreme disk profile indicates that our results are robust with respect to the assumed spatial distribution of γ rays. Our results and the bin-by-bin likelihood profiles provide the basis for comparisons with a variety of specific models that will help improve our understanding of CR physics in the Coma cluster.⁶⁶
• 3.  
  Based on scaling considerations and radio measurements, we derive a robust lower limit on the γ-ray luminosity of Coma that is a factor of ∼few below the median sensitivity given our analysis of blank fields.

We thank L. Rudnick for the radio continuum maps of Coma that we used in our study. We acknowledge useful discussions with F. Zandanel and G. Brunetti.

Y.R. gratefully acknowledges the hospitality extended to him during a visit to KIPAC, where this work was initiated.

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.

We would like to thank the DOE SLAC National Accelerator Laboratory Computing Division for their strong support in performing the large amount of simulations and repeated analysis used in the blank field study that were necessary for this work.

This research made use of APLpy.⁶⁷ We acknowledge the use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration and the use of HEALPix http://healpix.jpl.nasa.gov/ (Górski et al. 2005).

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

Pass 8 is the successor of the previous event level analysis (Pass 7REP; Ackermann et al. 2012a; Bregeon et al. 2013). The improvements include a tree-based pattern recognition algorithm to identify and reconstruct tracks, a new energy reconstruction that uses minimum spanning trees and calorimeter clustering, which also extend the energy range at both high and low energies and an improved background rejection algorithm from the anti-coincidence shield. Pass 8 also improved on the implementation and training of the classification trees used in the gamma-ray selection algorithms, resulting in increased gamma-ray selection efficiency with respect to Pass 7, while keeping the same background rejection power (Atwood et al. 2013).⁶⁸ For high-level science analysis this presents a significant improvement in all metrics. Among the key improvements are a ∼25% increase in the acceptance (and more than 50% below 100 $\mathrm{MeV}$ and above 300 $\mathrm{GeV}$) along with an improved angular resolution. As a result the point source sensitivity in the energy range between 1 and 10 $\mathrm{GeV}$ is enhanced by 30%–40% with respect to Pass 7REP. Similarly, the improved PSF provides increased sensitivity toward spatially extended sources.

Due to the nature of the changes introduced with Pass 8, there is only partial overlap between the Pass 7REP and Pass 8 event samples, with the Pass 8 sample generally being the larger of the two. The reason for this is two-fold: the changes in the event selection cause different events to pass the source selection criteria in Pass 8 than in Pass 7REP as well as selecting different residual non-photon events. The increase in γ-ray acceptance on the other hand increases the number of reconstructed events. Figure 4 shows the fraction of shared events between our here reported Pass 8 data set and the older Pass 7REP event sample after subjecting both data sets to the same ROI-based selection criteria. The total number of photons in our data set is about 612k. Below ∼3 $\mathrm{GeV}$, the fraction of shared events drops considerably, and toward the lowest energies there are about four times as many events in the Pass 8 data set that can partially account for the differences between previously published results based on Pass 7 and the analysis we present here.

Zoom In Zoom Out Reset image size

The basis for our background model is the 3FGL catalog of sources along with a set of diffuse templates rescaled to be used with Pass 8. We find several grid positions in the significance map (described in Section 2) that correspond to statistically significant (TS > 25) γ-ray excesses over the initial background model. The coordinates of four new point sources added to the background model are given in Table 3, along with their parameters as derived from a broadband spectral fit. For two sources we find that a LogParabola (LP) of the form $\displaystyle \frac{{dN}}{{dE}}={N}_{0}{\left(\frac{E}{{E}_{b}}\right)}^{-(\alpha +\beta \mathrm{log}(E/{E}_{b}))}$ gives a better fit to the data than a power law (PL); the remainder is modeled using a simple PL instead. The normalizations of each source is left free to vary in the likelihood fit.

Table 3.  Point Source Candidates in ROI

Name	R.A.	decl.	Model	Parameters	TS
	(°)	(°)
xFGL J1904.00+3465	190.40	34.65	PL	Γ = 1.88, E0 = 3496 $\mathrm{MeV}$	112
xFGL J1927.34+3129	192.73	31.29	LP	α = 1.83, β = 0.39, Eb = 2407 $\mathrm{MeV}$	86
xFGL J2027.04+2955	202.70	29.55	LP	α = 2.57, β = 0.41, Eb = 1237 $\mathrm{MeV}$	280
xFGL J1914.49+2080	191.45	20.80	PL	Γ = 3.45, E0 = 1000 $\mathrm{MeV}$	1220

Note. From left to right: name, R.A., decl., and assumed model along with fixed spectral parameters and TS values obtained in the likelihood fit to the Coma ROI. All coordinates are given in J2000 epoch.

Download table as:  ASCIITypeset image

Table 4 shows the best-fit values for the background-only model, including our new candidate sources. The associated uncertainty refers to the 68% parabolic error reported by MIGRAD. Note that the background model also comprises sources that are outside the selected 15° × 15° ROI, such as the bright source 3FGL J1224.9+2122 (4C +21.35), whose parameters are fixed to the reported values in 3FGL.

Table 4.  Best-fit Parameters of Background Model

Source Name	Normalization, N0
	(cm−2 s−1 $\mathrm{MeV}$−1)
3FGL J1230.3+2519	(2.11 ± 0.09) ×10−12
3FGL J1231.7+2847	(0.93 ± 0.04) ×10−12
3FGL J1254.5+2210	(1.10 ± 0.10) ×10−13
3FGL J1258.1+3233	(1.40 ± 0.10) ×10−12
3FGL J1258.4+2123	(0.30 ± 0.10) ×10−12
3FGL J1301.5+3333	(0.29 ± 0.05) ×10−12
3FGL J1303.0+2435	(2.80 ± 0.10) ×10−12
3FGL J1310.6+2446	(0.35 ± 0.09) ×10−13
3FGL J1310.6+3222	(2.60 ± 0.05) ×10−11
3FGL J1314.8+2349	(0.53 ± 0.03) ×10−12
3FGL J1321.0+2215	(4.20 ± 0.10) ×10−12
3FGL J1323.0+2942	(0.60 ± 0.03) ×10−12
3FGL J1326.1+2931	(0.40 ± 0.30) ×10−14
3FGL J1332.8+2723	(0.66 ± 0.09) ×10−12
xFGL J1914.49+2080	(1.50 ± 0.20) ×10−13
xFGL J1904.00+3465	(1.40 ± 0.40) ×10−14
xFGL J1927.34+3129	(0.60 ± 0.10) ×10−13
xFGL J2027.04+2955	(2.10 ± 0.40) ×10−13
Extragalactic Diffusea	(1.04 ± 0.01)
Galactic Diffusea	(1.04 ± 0.02)

Note.

^aThe fitted value corresponds to the overall (unit-less) normalization of an all-sky template. The nominal value is 1.0.

Download table as:  ASCIITypeset image