SEARCH FOR COSMIC-RAY-INDUCED GAMMA-RAY EMISSION IN GALAXY CLUSTERS

The joint likelihood is a source stacking technique, which has been previously applied with LAT data in the search for DM (Ackermann et al. 2011) and to study the extragalactic background light (Ackermann et al. 2012b). In brief, if the goal is to constrain or estimate a single or a set of parameters common to a source class, then backgrounds and individual properties of each source can be modeled individually and treated as nuisance parameters in source-specific likelihoods. The source-specific likelihoods can then be multiplied to yield a joint likelihood function that is used for inference on the common parameter of interest. The joint likelihood function for our case can be written as

$$\begin{equation} \mathcal {L}(A_\gamma |\mathcal {D}) = \prod _{i}\mathcal {L}_{i} (A_\gamma,b_{i}|\mathcal {D}_{i}), \end{equation} \tag{ 1 }$$

where A_γ is a (dimensionless) universal scale factor which serves as the parameter of interest and $\mathcal {D}$ refers to the photon data for all regions-of-interest (ROIs). The physical interpretation of the universal scale factor in terms of CR-induced γ-ray emission is discussed further in Section 4.1. The b_i parameters correspond to the parameters describing the background components in the individual ROIs and are treated as nuisance parameters in the likelihood evaluation. The b_i s include the normalizations of the isotropic and Galactic diffuse components. We denote the photon data for each ROI as $\mathcal {D}_{i}$. The index i runs over the ROIs in the sample. Having constructed the likelihood function, we use the profile likelihood method (e.g., Rolke et al. 2005) to obtain the best-fit values and confidence intervals for the parameter of interest, A_γ. We constrain our parameter of interest and nuisance parameters to be positive. The joint likelihood function is implemented in the Fermi Science Tools using the Composite2 package, and profiling over the likelihood function is achieved by means of MINOS, which is part of the MINUIT package (James & Roos 1975). The 90% confidence interval, $[A^{{\rm LL}}_\gamma, A^{{\rm UL}}_\gamma ]$ is defined as the values of A_γ where the log-likelihood has changed by 2.71/2 with respect to its value at the maximum, $-2\Delta \log {\mathcal {L}}=2.71$ (Bartlett 1953; Rolke et al. 2005). These intervals can be reinterpreted as upper limits at the 95% confidence level (C.L.), if the parameter is unconstrained in the fit, which we do if the lower limit ⩽0.

For quantifying the significance of a potential excess, we employ the common likelihood ratio approach (Neyman & Pearson 1928):

$$\begin{equation} {\rm TS}=-2\log \left(\frac{\mathcal {L}(A_{\gamma }=0,\hat{\hat{b}})}{\mathcal {L}(\hat{A}_\gamma,\hat{b})}\right), \end{equation} \tag{ 2 }$$

where $\mathcal {L}(A_\gamma =0,\hat{\hat{b}})$ is the null hypothesis, i.e., it represents the likelihood evaluated at the best-fit $\hat{\hat{b}}$ under the background-only hypothesis (Section 5.1) and $\mathcal {L}(\hat{A}_\gamma,\hat{b})$ is the likelihood evaluated at the best-fit value of ${\hat{A}}_{\gamma }$, $\hat{b}$, when including our candidate γ-ray (cluster) source.

As we constrain the signal fit parameter A_γ ⩾ 0, the null distribution of the test statistic, TS, is given by (1/2)δ + (1/2)χ² for one degree of freedom (Chernoff 1952), which we verified by MC simulations. For details, the reader is referred to Appendix A. While these simulations agree well with the expectations from Chernoff's theorem, studies based on random ROIs that encapsulate systematic effects in the LAT data (e.g., imperfect diffuse modeling, unresolved background sources, and percent-level inconsistencies in the IRFs) indicate that the significance estimated by simulations is probably somewhat too high when compared to the asymptotic expectations (Ackermann et al. 2014).

Recalling the discussion in Section 2, sources and ROIs have to be selected such that the overlap is minimal since the joint likelihood function Equation (1) does not account for correlation terms between the different individual likelihood functions. For technical reasons, we cannot define a single ROI containing all 50 clusters, as this leads to an overflow in the number of free parameters that MINOS is not able to handle. We therefore construct non-overlapping ROIs, each containing one or more (non-overlapping) cluster sources.

In order to avoid overlaps between ROIs, we perform an iterative procedure in which we treat each cluster in our sample with its extension as listed in Table 1 as a seed source and construct a circular region around it. We require these regions to be at least 8° in radius in order to obtain a good fit to the background model. If we cannot accommodate clusters in separate ROIs, we mitigate any remaining overlap by enlarging the ROIs such that more than one cluster may be contained in one analysis region.

Using this method, we are able to construct 26 independent ROIs containing the clusters in our sample which are listed in Table 2. We show the locations and extensions of our selected clusters and the ROIs containing them in Figure 2.

In this paper, we focus on the CR-induced γ-ray signal and defer an analysis of γ rays originating from DM decay or annihilation to future studies. For a study of γ rays originating from star-forming galaxies in galaxy clusters, see, e.g., Storm et al. (2012). Hydrodynamic simulations of galaxy clusters in a cosmological framework that include CR physics show an approximately universal spectral and spatial CR distribution within clusters as a result of hierarchical structure growth (Pinzke & Pfrommer 2010). The γ-ray emission induced by decaying neutral pions dominates over the inverse Compton emission from primary shock-accelerated electrons or secondaries injected in hadronic CR p–p interactions within clusters (Miniati 2003; Pfrommer et al. 2008; Pinzke & Pfrommer 2010).

Using a very simplified analytic model that employs a CR power-law energy spectrum dn_CR/dε∝ε^−α with spectral index α = 2 (i.e., equal CR energy density per logarithmic energy interval), Kushnir & Waxman (2009) claim that the IC emission from primary accelerated electrons at accretion shocks dominates over pion decay γ rays by a factor of ≃ 150 (ζ_e/ζ_p) (kT/10 keV)^−1/2, where ζ_e and ζ_p denote the fraction of shock-dissipated energy that is deposited in CR electrons and protons, respectively. Instead of the centrally concentrated pion decay emission, this model would give rise to a spatially extended emission reaching out to the accretion shocks beyond the virial radius.

However, there are two simplifying model assumptions that conflict with results from numerical cluster simulations and observations of non-thermal emission of clusters and supernova remnant shocks, rendering the conclusions about this apparent dominance of the primary inverse Compton emission questionable.

First, the assumed spectral index is in conflict with numerical simulations of cosmological structure formation that have shown a continuous distribution of Mach numbers with weaker flow and merger shocks being more numerous in comparison to strong shocks with Mach numbers exceeding $\mathcal {M}\gtrsim 6$ (Ryu et al. 2003; Pfrommer et al. 2006; Skillman et al. 2008; Vazza et al. 2009, 2011). This necessarily implies a softer effective spectral injection index of primary shock-accelerated particles of α_inj ≃ 2.3, which is also consistent with observed spectral indices of radio relics. In fact, elongated relics show a mean radio spectral index of 〈α_ν〉 = 1.3 (Feretti et al. 2012), which translates to a cooling-corrected spectral index of CR electrons of α = 2α_ν = 2.6 (which is a lower limit since the uncorrected injection index could be as high as α = 2α_ν + 1 = 3.6). Hence phenomenologically, the relevant shock strengths for radio relic emission are on average characterized by small Mach numbers of $\mathcal {M}\lesssim 3$ and inconsistent with hard CR indices of α = 2. It can easily be seen that the different power-law indices are indeed the reason for the strongly differing conclusions on the importance of primary inverse Compton emission. Electrons with an energy of 500 GeV can Compton upscatter CMB photons to γ-ray energies of around 1 GeV. Adopting the effective spectral injection index of primary electrons, α_inj ≃ 2.3, and assuming a post-shock temperature at a typical accretion shock of kT ≃ 0.1 keV, we find a flux ratio of the Kushnir & Waxman (2009) model and that by Pinzke & Pfrommer (2010) of (500 GeV/0.1 keV)^0.3...0.6 ≃ 800...6 × 10⁵.

Second, ζ_e/ζ_p ∼ 1 is inconsistent with the observed ratio for Galactic CRs, which has a differential energy ratio of K_e/p ≃ 0.01 at 10 GeV (Schlickeiser 2002) and observations of supernova remnants constrain the ratio ζ_e/ζ_p ∼ 0.001 (Edmon et al. 2011; Morlino & Caprioli 2012). Assuming universality of the shock acceleration process, statistical inferences about the CR distributions at the solar circle and non-thermal modeling of individual supernova remnants indicate that electron acceleration efficiencies are very subdominant in comparison to that of protons.

Hence, as a baseline model, we apply the simulation-based analytical approach for the CR distribution (Pinzke & Pfrommer 2010), which only requires the gas density profile inferred from X-ray measurements as input. Note that these simulations account only for advective CR transport, where the CRs may be tied to the cluster plasma via small-scale tangled magnetic fields with cored CR profiles as a consequence. Additional CR transport such as CR diffusion and streaming have been neglected. Furthermore, we neglect the potentially important contribution from re-accelerated CRs (e.g., Brunetti & Lazarian 2007, 2011).

The pion decay γ-ray flux above energy E as a result of hadronic CR interactions, $F_{\gamma,\mathrm{exp}}^{\mathrm{CR}} ({>}E)$, can be parameterized as

$$\begin{equation} F_{\gamma,\mathrm{exp}}^{\mathrm{CR}} ({>}E) = A_\gamma \,\lambda _{\pi ^{0}\rightarrow \gamma }\left({>} E\right) \int _V dV \kappa _{\pi ^{0}\rightarrow \gamma }(R), \end{equation} \tag{ 3 }$$

where the integral extends over the cluster volume V, $\lambda _{\pi ^{0}\rightarrow \gamma }\left(>E\right)$ is the spectral γ-ray distribution, and $\kappa _{\pi ^{0}\rightarrow \gamma }(R)$ is the spatial distribution, both given in Pinzke & Pfrommer (2010). The parameter, A_γ, is a dimensionless universal scale factor, common to all the clusters in the (sub)sample. The predicted γ-ray flux above the minimum energy threshold 500 MeV, $F_{\gamma,\mathrm{exp}}^{\mathrm{CR}} (E>500\,\mathrm{MeV})$ is calculated using the formalism in Pinzke et al. (2011). We tabulate these values in Table 1. The quoted values of $F_{\gamma,\mathrm{exp}}^{\mathrm{CR}} (E>500\,\mathrm{MeV})$ correspond to a maximal efficiency ζ_{p, max} = 0.5 for diffusive shock acceleration of CR ions at structure formation shocks which translate into A_γ = 1 with correspondingly smaller values of A_γ for smaller efficiencies (obeying however a nonlinear relation). For completeness, we show the CR formalism in Appendix B.

The cluster brightness profile is used to fit the emission from each cluster and is derived from the line-of-sight (l.o.s.) integral of the γ-ray emissivity

$$\begin{eqnarray} \nonumber S_\gamma (\psi,>E) &=&A_\gamma \,\lambda _{\pi ^{0}\rightarrow \gamma }\left({>}E\right) \nonumber\\ &&\times \int _{\Delta \Omega } d\Omega \int _\mathrm{l.o.s} dl \,\frac{\kappa _{\pi ^{0}\rightarrow \gamma }\left(R(l)\right)}{4\pi }, \end{eqnarray} \tag{ 4 }$$

where l is the l.o.s. distance in the direction ψ that the detector is pointing and $R(l)=\sqrt{l^2+D_a^2-2 D_a l\cos \Psi }$ is the cluster radius. Here D_a is the angular diameter distance from the Earth to the center of the cluster halo and cos Ψ ≡ cos θcos ψ − cos φsin θsin ψ, with θ being the azimuthal angle and ϕ the polar angle. The angular integration dΩ = sin θdθ dφ is performed over a cone centered around ψ. The spatial features of our model are described in more detail in Appendix B.

Outside the very center (r > 0.03R₂₀₀), this model predicts a rather flat CR-to-thermal pressure profile, i.e., 〈X_CR〉 ∼ const. Most of the emission is contributed from the region around the core radius, which is well outside the central parts. In order to compare the chosen spatial CR profile, we contrast our analysis of the simulation-based model with two additional configurations in which the CR profile is derived from a constant X_CR profile (ICM model) and a constant P_CR profile (flat model). The normalizations of these CR profiles are fixed by assuming that the total CR number within R₂₀₀ in our baseline model is conserved. ⁶⁵ For illustration purposes, we show the surface brightness profiles using our three spatial emission models in Figure 3 for the case of the (massive) Coma cluster and the much less massive cluster A400.

Zoom In Zoom Out Reset image size

In our framework, the derivation of A_γ also allows us to constrain the CR-to-thermal pressure ratio, X_CR, in the ICM using the virial mass and virial radius of each cluster in our sample, since A_γ∝〈X_CR〉, where 〈X_CR〉 = 〈P_CR〉_V /〈P_th〉_V and the brackets indicate volume averages. To this end, we make use of the set of 14 galaxy cluster simulations presented in Pinzke & Pfrommer (2010), which span almost two orders of magnitude in mass and include different dynamical states ranging from relaxed to merging clusters. We show the CR-to-thermal pressure ratio X_CR as a function of radius and cluster mass in Figure 4. X_CR decreases for smaller radius approximately inversely with gas temperature since a composite of CRs and thermal gas favors the gas component over CRs upon adiabatic compression (e.g., Pfrommer & Enßlin 2004a). At a fixed radius, X_CR has a negative trend with mass, that is mainly driven by the virial temperature scaling of the thermal pressure distribution. We have $\langle X_\mathrm{CR}\rangle \propto \langle C\rangle / \langle P_\mathrm{th}\rangle \propto \tilde{C} / kT_{200}\sim M_{200}^{-0.23}$, where C is the normalization of the CR distribution function and $\tilde{C} = C m_p/\rho$ denotes the dimensionless normalization, which scales with cluster mass as $\langle \tilde{C}\rangle _V \propto M_{200}^{0.44}$ (Pinzke & Pfrommer 2010) and partially offsets the virial mass scaling of the cluster temperature $kT_{200}\propto M_{200}^{2/3}$.

Zoom In Zoom Out Reset image size

To formalize these considerations, we fit an empirical relation of X_CR(R/R₂₀₀, M₂₀₀) to the simulated data points and obtain

$$\begin{eqnarray} X_\mathrm{CR} &=& 0.023 A_\gamma \left(\frac{R}{R_{200}}\right)^{0.369} \!\left(\frac{M_{200}}{10^{15}\,{M}_\odot }\right)^{-0.239 \left(\frac{R}{R_{200}}\right)^{-0.258}}\!,\quad\ \ \ \end{eqnarray} \tag{ 5 }$$

where A_γ refers to our universal scale factor derived in the joint likelihood analysis. We tabulate the values for 〈X_CR〉 for two different integration radii, R₂₀₀ and the half-light radius R_HL assuming A_γ = 1 in Table 1. ⁶⁶

The background model for each ROI in this analysis includes templates for the diffuse Galactic and isotropic emission components as well as individual γ-ray sources reported in the 2nd Fermi-LAT catalog (Nolan et al. 2012; 2FGL). We model extended 2FGL sources according to the spatial templates provided by the Fermi Science Support Center. Unless stated otherwise, we use the standard diffuse and extragalactic γ-ray background templates that are recommended for performing data analysis of reprocessed LAT data. ⁶⁷ We note that the Galactic diffuse emission model we use includes a residual component of diffuse γ-ray emission that is not modeled by any template. This component is smoothly varying and does not contribute importantly to the intensity >500 MeV in the regions considered here. ⁶⁸

In the background model for each ROI, we include the union of 2FGL sources within the ROI radius enlarged by 5° and 2FGL sources located within 10° of any galaxy cluster in the ROI.

Since the average virial radius of clusters in our sample is less than 2°, we leave the normalizations of all sources free that are contained within a 4° radius around each cluster. In addition, we allow the normalization of the templates used to model the Galactic foreground and isotropic diffuse emission to vary freely.

One shortcoming of using the 2FGL catalog (based on 2 yr of LAT observations) to search within a data set covering 4 yr is that spectral parameters (in particular for variable sources) may have substantially changed. To account for this variability, we free the normalization of sources that coincide with bright spatial residuals, and we determine their values through performing a fit using a background-only model to obtain the best-fit for the null-hypothesis.

This procedure produces a large number of free parameters which are then fixed to their maximum likelihood values from the background-only model fit when maximizing the joint likelihood function introduced in Section 3.1 in order to avoid an overflow of free parameters and to ensure convergence of the maximization. The normalizations of the Galactic and isotropic diffuse components are left free in each ROI for the joint likelihood fit.

Our ROIs are well described individually by the null hypothesis, i.e., despite the increase in data volume, our results are consistent with emission from only previously detected individual Fermi-LAT sources and diffuse emission provided by the Galactic and extragalactic diffuse models, respectively, with the exception of three ROIs: the ROIs containing A400 and, less prominently A1367 and A3112 exhibit residual emission located within the virial radius of each respective cluster. We leave these excesses unmodeled in our baseline analysis and address the interpretation of these residuals in Section 5.3. It is also reassuring that the fitted normalizations of the two global diffuse backgrounds are narrowly scattered around the nominal value of 1 for both the extragalactic (isotropic) and the Galactic diffuse component (refer to Figure 5) across our entire sample. The isotropic component shows a slight bias toward normalizations >1 which however has negligible effect on our results. Figure 6 shows a stacked residual significance map for the full sample and CC and NCC sub-samples. These residual significance maps are created from theoretical model maps of predicted counts from the best-fit null hypothesis model summed over each cluster location. The combined model maps M are then subtracted from the stacked counts maps C from the same region and the residual significance R is computed as $R=(C-M)/\sqrt{M}$. For the spectral residuals, we refer the reader to Appendix E. No obvious excess is visible in either the spectral or spatial residuals.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We then repeat the fitting procedure including a model of our galaxy clusters with the predicted γ-ray flux F_exp(E > 500 MeV) and using the spatial template and spectral form proposed by Pinzke et al. (2011), leaving only A_γ to vary freely. We show the distribution of associated TS values for the respective samples in Figure 7.

Zoom In Zoom Out Reset image size

Assuming that backgrounds are properly modeled, and that CR physics governing the γ-ray emission of clusters is indeed universal, we calculate the best-fit value of the combined scale factor for the full sample along with the two morphological sub-samples. The global TS value of the scale factor for the full sample of 50 clusters is 7.3, corresponding to a 2.7σ evidence. We note that the largest contributors to this signal are however the aforementioned excesses which are spatially coincident with A400 but also less prominently from excesses toward A1367 and A3112. We discuss these special cases in the next section. Removing these clusters from the sample results in a drop of the significance below 2σ, yielding an upper limit to the common scale factor $A_\gamma ^{{\rm UL}}=0.29$ at 95% C.L. for the whole cluster sample containing the remaining 47 galaxy clusters. While individually the excess toward A400 yields a higher significance than the excesses toward either A3112 or A1367, in the combined limit, the contribution of this excess becomes negligible due to the lower flux prediction (which mainly determines the weight assigned to each cluster in the joint analysis), as compared to, e.g., the contribution from Coma that dominates the upper limit.

Since at present we can neither claim nor refute the origin of the observed excesses being due to γ rays from the ICM, we calculate upper limits on the universal scaling factor, leaving those respective excesses unmodeled. For our whole sample, we find a combined limit of $A^{{\rm UL}}_\gamma =0.41$ at 95% C.L. Considering the NCC/CC subsamples, we find weaker limits $A^{{\rm UL}}_\gamma =0.47$ on the scale factor for NCC systems as compared to $A^{{\rm UL}}_\gamma =0.49$ for the CC subsample.

We also calculated the limit on the combined scale factor for our two alternative spatial CR profiles. For the ICM model, we obtain $A^{{\rm UL}}_\gamma =0.48$, and for the flat model, the combined limit is roughly a factor of four larger with respect to the results from the baseline model, yielding $A^{{\rm UL}}_\gamma =1.78$, at 95% C.L. The associated global TS values are 7.2 and 9.7, respectively. ⁶⁹ In addition, a flat CR profile is preferred in the case of the Coma cluster which is also the cluster that contributes most to the constraints in the NCC cluster subsample. We provide the final values for our three setups in Table 3.

For the combined sample, we exclude A_γ = 1 at more than a 5 σ confidence level, while for CC/NCC, it is excluded at more than a 4 σ confidence level.

We obtained a TS value of 13.8, corresponding to a pre-trial factor significance of 3.7σ individually for each of the candidate γ-ray sources at the locations of A1367 (Leo cluster) and A3112. Conservatively, assuming a binomial distribution for the trial probability, we find that these excesses correspond to a post-trial significance of 2.6σ. For each of these regions, we calculated a TS map using an unbinned likelihood method in a 5° × 5° region centered around each cluster with a grid spacing of 01 between test positions. The TS value of a putative point-like source with a hadronic CR spectrum is evaluated at each test position on the grid to create a spatial map of the excess emission. We show these TS maps in Figure 8.

Zoom In Zoom Out Reset image size

In both cases, we find that the excess emission, albeit marginally offset from the center of the cluster (03 for A1367 and 01 for A3112) is still contained within the virial radius of the respective cluster. For A1367 the difference in TS obtained at the center of the cluster and the peak TS position as shown in Figure 8 is 15, while for A3112 the offset is smaller than the resolution used to create the TS maps. In order to test whether the emission is more appropriately modeled assuming an extended emission profile over a new point source, we follow the analysis presented in Lande et al. (2012), which gauges the spatial extent of the source by taking the difference, TS_ext between the TS for a point source signal hypothesis and the extended source signal hypothesis. Lande et al. (2012) found that sources with TS_ext > 16 are confidently ascertained to be spatially extended beyond the LAT point-spread function (PSF). ⁷⁰ For A1367 and A3112, we find TS_ext = 2.6 and TS_ext = 0.9, respectively. Assuming that the excesses toward all three clusters constitutes point-like emission, we first obtained a better estimate for the location of the excess using the gtfindsrc tool and then repeated the calculation of TS maps using a finer binning of 002. We report the best-fit values of these excesses along with their uncertainties in Table 4.

We also investigate the spectral behavior of the excesses toward A1367 and A3112, by replacing the hadronic CR spectrum with a featureless power-law, such that the flux becomes:

$$\begin{eqnarray} \nonumber F_{\gamma }^{\mathrm{PL}}\left(E\right)&=& A_\gamma \times F_{\gamma,\mathrm{exp}}^{\mathrm{CR}} ({>} 500\,\mathrm{MeV})\nonumber\\ &&\times \frac{(1-\Gamma)\times E^{-\Gamma }}{E_{\mathrm{max}}^{1-\Gamma }-E_{\mathrm{min}}^{1-\Gamma }}, \end{eqnarray} \tag{ 6 }$$

where Γ is the spectral index and E_min = 500 MeV and E_max = 200 GeV correspond to the energy range of the analysis. $F_{\gamma,\mathrm{exp}}^{\mathrm{CR}} (E>500\,\mathrm{MeV})$ denotes the expected integral flux above E_min. A_γ corresponds to the scale factor introduced in Section 5.2. For this test, we leave both A_γ and Γ free to vary. We report the best-fit values of these parameters in Table 5.

Table 5. Spectral Model Comparison

Cluster	Aγ	TS	Aγ	Γ	TSPL
	(Hadronic Model)	(Hadronic Model)	(Power Law)	(Power Law)	(Power Law)
A1367	$3^{+2}_{-1}$	13.8	1.7 ± 0.9	1.7 ± 0.3	17.3
A3112	3 ± 2	13.8	2.0 ± 1.5	1.7 ± 0.4	16.1
A400	$39^{+11}_{-10}$	52.7	43 ± 8	2.3 ± 0.2	52.8

Notes. Spectral model comparison of clusters that exhibit excess emission. Shown are the best-fit values for A_γ along with their associated TS values for the hadronic model and the corresponding values for A_γ when replacing the spectrum by a featureless power-law of index Γ, given by Equation (6). The last column indicates the obtained TS value with the power-law fit.

Download table as:  ASCIITypeset image

For both A1367 and A3112, we find that harder spectral indices are preferred over softer, with a best-fit value for Γ = 1.7 ± 0.3 and Γ = 1.7 ± 0.4, for A1367 and A3112, respectively.

While none of the aforementioned tests decisively exclude the attribution of the observed excesses to CR-induced γ-ray emission from the ICM, we note that both A3112 and A1367 are hosts to head-tail radio sources which may be the source of the observed γ-ray emission (see, e.g., Gavazzi & Jaffe 1982, 1987 for A1367 and Costa & Loyola (1998) for A3112). A discussion of supporting multifrequency arguments is given in Appendix D.

The excess found in A400 yields a TS value of 52.7 which nominally corresponds to a significance of 7.3σ pre-trial (6.7σ post-trial). We performed the same tests as in the case of A1367 and A3112 and find that the excess is contained within the cluster virial radius, although the TS map indicates that the emission is about 03 offset from the cluster center. The difference between the TS evaluated at the cluster center and at the position of the excess is 38. The test for extendedness (centered at the cluster center) yields TS_ext = 15.0, indicating a preference for an extended source, which is likely explained by the offset of the excess. This casts further doubt on the explanation of the excess being ICM emission that moreover would be concentrated toward the cluster center. Fitting the excess with a power-law spectral model as in the case of A3112 and A1367 yields a best-fit value for Γ = 2.3 ± 0.2 and A_γ = 43 ± 8 with an associated TS value of 52.8, which is similar to that obtained for the hadronic model.

In addition to these tests, we searched for source variability using aperture photometry and found no indications of variability on time scales of one month. We note, however, that the obtained scale factor for A400, $A_\gamma =39_{-10}^{+11}$, is in strong tension with baseline model expectations and the scale factor constraints derived from other clusters in our sample. If the excess toward A400 constitutes a signal, the calculated upper limit from A400 corresponds to the usual statement that scale factors larger than the upper limit are inconsistent with the data on the stated confidence level. If the excess instead stems from the unmodeled background, then the upper limit means that scale factors larger than the upper limit would make the model more inconsistent with the data than the background-only hypothesis allows at the stated confidence level. In both cases, the result is a conservative upper limit and in addition (as mentioned in Section 5.2) A400 does not affect the combined upper limit sizeably. Removing A400 from the sample yields a marginally stronger upper limit on the common scale factor for the combined sample of the remaining 49 clusters of $A_\gamma ^{{\rm UL}}=0.40$.

In the conservative CR model of Pinzke & Pfrommer (2010), A_γ = 1 corresponds to a pion decay γ-ray flux owing to CR protons accelerated at structure formation shocks with a maximal acceleration efficiency and neglecting active CR transport. ⁷¹ Allowing for CR streaming transport would cause a net CR flux to the dilute outer cluster regions and would reduce the γ-ray yield for that acceleration efficiency. Additionally, Pinzke & Pfrommer (2010) presented an optimistic CR model including effects which enhance the predicted γ-ray yield by a factor of 2–3 depending on the cluster mass. ⁷² However, that model does not account for CRs that are injected by AGNs over the cluster lifetime, which could also produce pions in inelastic collision with the ICM. The total energy dissipation by gravitational shocks exceeds that of AGNs for large clusters, making it unlikely that AGN-injected CRs dominate the diffuse γ-ray signal. However, this argument does not apply for smaller CCs (in particular in their central cluster regions) where the AGN appears to dominate the energy budget and could possibly also give rise to observable γ-ray emission (see, e.g., Pfrommer 2013, for the interpretation of the low-state of M87 in terms of diffuse pion decay emission while the high state is attributed to jet-induced emission). Hence, the most likely explanation for the possible signal with A_γ ≃ 39 in A400 is jet-related emission from point sources projected onto or within the cluster. A potential source for the emission in A400 may be the quadruple head–tail system 3C 75 which also shows X-ray core emission in the galaxies (Owen et al. 1985; Hudson et al. 2006). However, given that 3C 75 is located toward the center of A400 and the excess is 03 offset, which is about eight times larger than the 68% error radius for a pointsource, this possibility is unlikely.

Assuming that each cluster in our sample can be modeled according to our description in Section 4.1, one can also derive individual limits on A_{γ, i} where i refers to the individual cluster. These individual limits can be propagated into flux limits.

In addition, in order to assess the impact of modeling the clusters in our sample as extended sources, we have derived individual flux upper limits modeling the clusters as point sources.

In Figure 9, we show these two cases, contrasting the individually derived γ-ray flux limits for extended emission from those derived when assuming the cluster emission to be point-like. We tabulate the former for various energy bands along with their associated (pre-trial) TS values in Table 6. We note that the limit on A_γ, derived from the Coma cluster alone is comparable to the jointly derived limit from the full sample of all 50 clusters in this study, emphasizing its weight in the joint analysis. We investigated the potential dependence of the upper limits on the Galactic diffuse emission model. In the great majority of cases the limits vary by less than 30% for the range of models that we considered. Those clusters for which the dependence on the model was more sensitive are marked in Table 6. However, this sensitivity does not affect the combined limit.

Zoom In Zoom Out Reset image size

Table 6. Individual Flux Upper Limits

Cluster	$A_{\gamma,i}^{{\rm UL}}$	$F^{{\rm UL}}_{\gamma,500\,\mathrm{MeV}}$	$F^{{\rm UL}}_{\gamma,1\,\mathrm{GeV}}$	$F^{{\rm UL}}_{\gamma,10\,\mathrm{GeV}}$	TS	$A_{\gamma,i}^{{\rm UL}}$	$F^{{\rm UL}}_{\gamma,500\,\mathrm{MeV}}$	$F^{{\rm UL}}_{\gamma,1\,\mathrm{GeV}}$	$F^{{\rm UL}}_{\gamma,10\,\mathrm{GeV}}$	TS
		(× 10−10)	(× 10−11)	(× 10−12)			(× 10−10)	(× 10−11)	(× 10−12)
	Extended	Extended	Extended	Extended	Extended	Point-like	Point-like	Point-like	Point-like	Point-like
2A0335	0.52	1.5	6.7	3.6	0.0	0.45	1.3	5.9	3.1	0.0
A0085	1.35	3.9	18.0	9.5	1.0	1.26	3.7	16.9	9.0	1.3
A0119	3.54	5.0	23.0	12.2	2.7	2.65	3.8	17.2	9.1	0.8
A0133	2.32	1.7	7.6	4.0	0.0	2.23	1.6	7.3	3.9	0.0
A0262	1.67	2.0	9.3	4.9	0.0	1.07	1.3	5.9	3.1	0.0
A0400	50.11	22.2	101.7	53.6	52.7	41.20	18.3	83.6	44.1	37.7
A0478 a	1.14	2.8	12.7	6.7	0.0	0.88	2.1	9.8	5.2	0.0
A0496 a	1.94	5.5	25.2	13.3	1.9	1.46	4.1	18.9	10.0	0.5
A0548e	10.78	2.7	12.3	6.5	0.1	9.02	2.2	10.3	5.4	0.0
A0576	4.39	2.7	12.6	6.6	0.2	3.80	2.4	10.9	5.7	0.3
A0754	2.86	4.8	22.2	11.7	1.1	2.49	4.2	19.3	10.2	0.7
A1060	1.89	4.3	19.7	10.4	0.4	1.28	2.9	13.3	7.0	0.1
A1367	4.08	7.8	35.9	19.0	13.8	3.31	6.4	29.1	15.4	11.2
A1644	2.69	3.5	16.0	8.5	0.2	2.41	3.1	14.3	7.6	0.4
A1795 a	0.42	1.3	5.8	3.1	0.0	0.42	1.3	5.7	3.0	0.0
A2065	5.00	4.6	21.1	11.2	2.3	4.69	4.3	19.8	10.5	2.2
A2142	0.47	1.6	7.4	3.9	0.0	0.46	1.6	7.3	3.9	0.0
A2199	1.45	4.3	19.8	10.5	1.9	1.32	4.0	18.1	9.6	1.7
A2244	1.74	1.3	6.0	3.2	0.0	1.63	1.2	5.6	3.0	0.0
A2255	5.21	4.4	20.2	10.6	5.2	4.99	4.2	19.3	10.2	6.3
A2256	0.50	1.3	5.8	3.1	0.0	0.42	1.1	4.9	2.6	0.0
A2589 a	3.96	2.7	12.2	6.5	0.0	3.36	2.3	10.4	5.5	0.0
A2597	1.27	1.0	4.4	2.3	0.0	1.06	0.8	3.7	2.0	0.0
A2634	5.95	3.7	17.0	8.9	0.3	4.16	2.6	11.9	6.2	0.0
A2657 a	2.34	1.9	8.9	4.7	0.0	1.75	1.4	6.6	3.5	0.0
A2734	2.71	1.2	5.5	2.9	0.0	2.58	1.1	5.2	2.7	0.0
A2877	1.87	0.9	4.3	2.3	0.0	1.55	0.8	3.6	1.9	0.0
A3112	5.37	6.0	27.3	14.4	13.8	5.06	5.6	25.7	13.6	12.8
A3158	1.26	1.5	6.9	3.7	0.0	1.17	1.4	6.4	3.4	0.0
A3266	1.24	3.9	17.7	9.4	2.2	1.15	3.6	16.4	8.7	3.6
A3376	6.20	4.3	19.5	10.3	0.0	3.59	2.5	11.3	6.0	0.0
A3822	4.77	2.1	9.8	5.2	0.0	4.22	1.9	8.7	4.6	0.0
A3827	0.66	0.6	2.9	1.5	0.0	0.61	0.6	2.7	1.4	0.0
A3921 a	4.63	2.2	9.9	5.2	0.1	4.11	1.9	8.8	4.6	0.0
A4038	1.84	2.5	11.3	6.0	0.1	1.48	2.0	9.1	4.8	0.0
A4059	2.04	2.0	9.1	4.8	0.0	1.86	1.8	8.2	4.4	0.0
COMA a	0.35	4.0	18.5	9.8	0.7	0.22	2.5	11.3	6.0	0.1
EXO0422	0.70	0.5	2.3	1.2	0.0	0.65	0.5	2.2	1.1	0.0
FORNAX	3.73	3.1	14.3	7.6	0.0	1.29	1.1	4.9	2.6	0.0
HCG94	6.34	2.4	11.0	5.8	0.0	5.17	2.0	9.0	4.7	0.0
HYDRA-A	2.57	4.3	19.6	10.4	3.1	2.44	4.1	18.6	9.8	3.2
IIIZw54	10.97	5.0	22.7	12.0	0.5	9.74	4.4	20.2	10.6	0.3
IIZw108	3.31	1.3	6.1	3.2	0.0	3.25	1.3	6.0	3.2	0.0
NGC 1550	3.37	2.3	10.3	5.5	0.0	2.70	1.8	8.3	4.4	0.0
NGC 5044	4.96	3.0	13.8	7.3	0.0	4.29	2.6	11.9	6.3	0.0
RXJ2344	5.07	2.9	13.3	7.0	1.3	5.10	2.9	13.4	7.1	2.3
S405 a	3.11	1.2	5.7	3.0	0.0	2.36	0.9	4.3	2.3	0.0
S540	11.47	4.0	18.5	9.7	1.2	10.35	3.6	16.7	8.8	0.8
UGC03957	8.00	4.0	18.2	9.6	0.4	7.31	3.6	16.6	8.8	0.3
ZwCl1742	1.94	2.3	10.4	5.5	0.0	1.77	2.1	9.5	5.0	0.0

Notes. The columns contain from left to right the 95% upper limit on A_{γ, i} for each cluster, the derived flux upper limit above (500 MeV, 1 GeV, and 10 GeV), the associated TS value as well as the same quantities assuming the clusters to be modeled by a point source at the cluster position as given in Table 1. Fluxes are given in photon s⁻¹ cm⁻². ^aThe upper limits on A_{γ, i} derived using our alternative diffuse models (Section 6.3) for these clusters varied by more than 30% relative to those obtained using the standard diffuse emission model.

Download table as:  ASCIITypeset image

In Figure 10, we show individual γ-ray flux upper limits that are derived for the CR profile following: (1) a constant X_CR profile (ICM model) and (2) a constant P_CR profile (flat model). We tabulate these values along with their (pre-trial) TS values in Table 7. We find that the limits derived from the first model are very similar to those with our simulation-based model. On the other hand, the flat model (i.e., constant CR pressure profile) yields less stringent upper limits. We also note that the associated TS values for the flat CR profile are generally marginally higher than either those derived from the thermal ICM or our simulation-based approach. This is expected since the choice of a flat CR profile yields a substantially flatter γ-ray surface brightness profile, which in turn provides a better fit to the data compared to a cored profile, in particular, when considering that the excess emission we report in the previous section is offset from the cluster center. We also note that from the three excesses found, only those toward A3112 and A400 remain with TS > 9 when using the flat CR profile instead.

Zoom In Zoom Out Reset image size

The flux limits we derive substantially improve over previous limits (Ackermann et al. 2010b) due to the increase in data volume and improved modeling of the γ-ray sky as well as improved instrument understanding reflected by the use of reprocessed LAT data with on-orbit calibrations. Finally, we note that in particular for spatially extended clusters such as Coma or Fornax, the limits are substantially weakened with respect to when modeling them as point-like objects.

Motivated by the Fermi-LAT detection of few bright cluster galaxies (e.g., 2FGL J0627.1−3528 in A3392 and 2FGL J1958.4-3012 in RXC J1958−3011), a recent stacking study by Dutson et al. (2013) investigated a large sample of galaxy clusters that was selected according to the radio flux of bright cluster galaxies. Although based on a different scientific prior and methodology than our cluster analysis, the determined flux limits can be compared to the point-source upper limits reported here. About a dozen clusters are in common between both studies. However, Dutson et al. (2013) used the respective coordinates of the bright cluster galaxy, which are not necessarily consistent with the cluster center coordinates considered in our studies. Given the LAT PSF (Bregeon et al. 2013) and the considered ROIs this does not constitute a severe handicap for comparison. The use of different exposures (45 months in Dutson et al. 2013 and 48 months in this work), the use of different source models, and, perhaps most particularly noteworthy, the use of reprocessed LAT data with associated different Galactic and isotropic diffuse models (gal_2yearp7v6_v0 versus gll_iem_v05 and iso_p7v6source versus iso_clean_v05 respectively), as well as different analysis energy thresholds, render a strict comparison more problematic. For the majority of common clusters, the limits in Dutson et al. (2013) are marginally less sensitive, as expected regarding the slightly less exposure and the rather moderate changes in the diffuse background models.

However, there are two noticeable exceptions: A85 and A2634 appear to have more constraining upper limits in Dutson et al. (2013), besides less exposure and a lower analysis threshold. The discrepancies could be explained by differences in the construction of the ROI (treatment of variable sources, sources to be too faint to be in the 2FGL catalog, size of ROI). Taking the respective flux limits at face value, the differences do not amount to more than 35% between both studies. Our limits on extended cluster emission cannot be meaningfully compared to the point-source upper limits in Dutson et al. (2013) as they constitute alternative scientific priors for a different scientific problem. However, our individual limits on the γ-ray flux, while being specifically derived within the framework of the universal CR model by Pinzke & Pfrommer (2010), can, in principle, be used to constrain other classes of models.

We show the resulting upper limits on the CR-to-thermal pressure ratio, 〈X_CR〉 in Figure 11. These numbers were obtained by scaling the 〈X_CR〉 values in Table 1 with the limit on A_γ from Section 5.2. This procedure assumes universality of the CR distribution as suggested by hydrodynamical cosmological simulations of clusters (Pinzke & Pfrommer 2010) and implicitly asserts that the active CR transport does not appreciably modify the spatial distribution of the CRs. This is justified since the impact of CR streaming on the CR distribution of a cosmological cluster is not clear to date. Depending on the microscopic plasma physics that sets the CR streaming speed (i.e., competing damping mechanisms of the CR Alfvén waves) and the macroscopic distribution of cluster magnetic fields (Pfrommer & Dursi 2010; Ruszkowski et al. 2011, for evidence of radial bias of the magnetic geometry), CR streaming could either be a perturbation to the peaked advection-dominated CR distribution or cause a substantial flattening by a substantial net outward CR flux (Enßlin et al. 2011; Wiener et al. 2013).

Zoom In Zoom Out Reset image size

The median upper limit on the CR pressure ratio is 〈X_CR〉 < 0.006 for the combined sample and the NCC subsample within R_HL. Those constraints are relaxed to 〈X_CR〉 < 0.012 (combined sample) and 〈X_CR〉 < 0.013 (NCC) within R₂₀₀. For CC clusters, this limit is less stringent, yielding 〈X_CR〉 < 0.008 and 〈X_CR〉 < 0.014 within R_HL and R₂₀₀, respectively.

These limits are more constraining than the previous limits on the CR pressure that were obtained through flux upper limits on individual objects; in particular, they are constraining for individual limits on clusters using the initial 18 months of LAT data (Ackermann et al. 2010b); also, the limits improve those constraints that use four years of Fermi data on Coma, which yield 〈X_CR〉 < 0.017 (Arlen et al. 2012), provided the CR universality assumption holds. ⁷³ The most suitable cluster target for CR-induced γ-ray emission, the Perseus cluster, cannot be used to competitively constrain the CR pressure using Fermi-LAT data since the Fermi-LAT and MAGIC collaborations detected the central radio galaxy NGC 1275 in γ rays in the energy range from 300 MeV to >300 GeV (Fermi: 300 MeV–300 GeV, MAGIC: >200 GeV) (Abdo et al. 2009; Aleksić et al. 2012b). Non-observations of γ rays from the Perseus cluster above these energies by the MAGIC Collaboration (Aleksić et al. 2012a, 2010b) provide limits similar to those obtained from analyzing LAT observations of Coma, 〈X_CR〉 < 0.017 alone (Arlen et al. 2012; Zandanel & Ando 2013).

Our limits on X_CR probe the entire ICM and are much more constraining in comparison to the limits on the non-thermal pressure contribution of the central ICM in several nearby cD galaxies that have been derived by comparing the gravitational potentials inferred from stellar and globular cluster kinematics and from assuming hydrostatic equilibrium of the X-ray emitting gas (Churazov et al. 2008, 2010). Depending on the adopted estimator of the optical velocity dispersion, they find a non-thermal pressure bias of X_nt = P_nt/P_th ≈ 0.21–0.29, which probes the cumulative non-thermal pressure contributed by CRs, magnetic fields, and unvirialized motions. It is, however, conceivable that those central regions of CC clusters, which probe the enrichment of non-thermal components as a result of AGN feedback (e.g., Pfrommer 2013), are characterized by a larger CR pressure contribution in comparison to the bulk of the ICM that probes CRs accelerated by shocks associated with the growth of structure and magnetic fields that are reprocessed by these shocks.

Complementary limits on CRs are also derived from radio (synchrotron) observations (Pfrommer & Enßlin 2004a; Brunetti et al. 2007). Assuming a central cluster magnetic field of B ∼ 1 μG and CR spectral indices α = 2.1–2.4, this approach allows the CR energy to be constrained to a few percent while the limits are less stringent for steeper spectra and lower magnetic fields (Aharonian et al. 2009a).

The distributions of CRs within the virial regions of clusters are built up from shocks during the cluster assembly (Ryu et al. 2003; Pfrommer et al. 2006; Pinzke et al. 2013). While strong shocks are responsible for the high-energy population of CRs that could potentially be visible at TeV γ-ray energies, intermediate Mach-number shocks with $\mathcal {M}\simeq 3$–4 build up the CR population at GeV energies (Pinzke & Pfrommer 2010), which can be constrained through Fermi observations. The normalization of both the CR population and the γ-ray flux scale with the acceleration efficiency at shocks of corresponding strength. Our fiducial CR model (Pinzke & Pfrommer 2010) uses a simplified model for CR acceleration (Enßlin et al. 2007) in which the efficiency rises steeply with Mach number for weak shocks and saturates already at shock Mach numbers $\mathcal {M}\gtrsim 3$. Observations of supernova remnants (Helder et al. 2009) and theoretical studies (Jones & Kang 2005) suggest a value for the saturated acceleration efficiency of ζ_{p, max} ≃ 0.5, i.e., the fraction of shock-dissipated energy that is deposited in CR ions. Following these optimistic predictions, the model of Pinzke & Pfrommer (2010) assumes this value for the saturated efficiency, which serves as an input to our analysis. This model provides a plausible upper limit for the CR contribution from structure formation shocks in galaxy clusters since more elaborate models of CR acceleration predict even lower efficiencies for the Mach-number range $\mathcal {M}\simeq 3$–4 of relevance here (Kang & Ryu 2013, see also Appendix C for a more detailed discussion). If we scale the CR pressure contribution linearly with the maximum acceleration efficiency, using the previously derived limit on A_γ, we find $\zeta _{\mathrm{p,max}}^{{\rm UL}}=21\%$ for the combined sample while the CC- and NCC-subsamples yield 25% and 24% maximum acceleration efficiency, respectively. We note that this constrains the maximum acceleration efficiency only in the simplified model adopted here. For different acceleration models, these upper limits provide conservative constraints on the acceleration efficiency at intermediate strength shocks of Mach number $\mathcal {M}\simeq 3$–4, a regime complementary to that studied at supernova remnant shocks.

However, these conclusions rely on two major assumptions, namely, CR universality and the absence of efficient CR transport relative to the plasma rest frame. The latter assumption hypothesizes that CRs are tied to the gas via small-scale tangled magnetic fields, which implies that they are only advectively transported and that we can neglect the CR streaming and diffusive transport relative to the rest frame of the gas. Early work on this topic suggests that such CR transport processes are at work in clusters and cause a flattening of radial CR profiles that can significantly reduce the radio and γ-ray emission at high energies probed by Cherenkov telescopes but remains largely unaffected at lower energies probed by Fermi (Wiener et al. 2013). Moreover, different formation histories of clusters cause spatial variations of the CR distribution and hence a deviation from a universal distribution (Pinzke & Pfrommer 2010). To date, there is no consensus about the size of these effects, and more work is needed to fully quantify them.

6.1.1. Binning

6.1.2. Region Size

6.1.3. Free Sources

The IRFs consist of three separate parts (see Ackermann et al. 2012a, for details): the effective area which has an associated uncertainty of at most ∼10% in the energy range we consider, the PSF whose uncertainty can be conservatively estimated to be ∼15%, and the energy dispersion (whose uncertainty is negligible for this analysis). Using bracketing IRFs (see Sections 5.7.1 and 6.5.1 in Ackermann et al. 2012a) to quantify these uncertainties, we find that while individual ROIs may show variations of up to ∼21%, the effect on the combined scale factors and quantities derived from it is less than 7%.

Spatial residuals due to mismodeling of the large-scale Galactic diffuse foreground emission may be misinterpreted in terms of an extended γ-ray excess. We compare results derived using the standard diffuse emission model adopted for the baseline analysis (based on empirical fits of multiple spatial templates to γ-ray data) to results obtained when using a set of eight alternative diffuse emission models that were created using a different methodology with respect to the standard diffuse emission model. ⁷⁴

We chose these models to represent the most important parameters scanned in Ackermann et al. (2012c), in particular, CR source distribution, halo size, and spin temperature. We summarize the properties of the alternative models in Table 8, and we refer readers to de Palma et al. (2013) for details. The models we employed were tuned to the P7REP data. Although the models were created such that different components along the l.o.s. could be fit separately, we only adopted a free overall normalization, since at high Galactic latitudes the vast majority of the gas resides in the neighborhood of the solar system. Moreover, having different components as additional degrees of freedom in the fit makes a comparison of the TS values with the baseline analysis more difficult.

Table 7. Individual Flux Upper Limits (Alternative CR Profiles)

Cluster	$A_{\gamma,i}^{{\rm UL}}$	$F_{\gamma,\mathrm{exp}}^{\mathrm{CR}}$	$F^{{\rm UL}}_{\gamma,500\,\mathrm{MeV}}$	TS	$A_{\gamma,i}^{{\rm UL}}$	$F_{\gamma,\mathrm{exp}}^{\mathrm{CR}}$	$F^{{\rm UL}}_{\gamma,500\,\mathrm{MeV}}$	TS
Spatial Model	ICM	ICM	ICM	ICM	Flat	Flat	Flat	Flat
2A0335	0.85	1.8	1.5	0.0	5.14	0.4	2.2	0.0
A0085	1.64	2.4	3.9	0.6	6.15	0.8	4.6	1.1
A0119	3.91	1.3	5.2	2.7	9.55	0.6	6.1	3.0
A0133	3.83	0.5	2.0	0.0	15	0.1	1.9	0.0
A0262	2.39	1.0	2.3	0.0	5.63	0.7	3.8	0.0
A0400	59.23	0.4	22.5	52.9	95.1	0.3	24.5	61.4
A0478	1.57	1.9	3.0	0.0	6.9	0.4	2.6	0.0
A0496	2.74	2.1	5.7	2.3	14.66	0.6	9.0	5.4
A0548e	12.8	0.2	2.7	0.1	17.91	0.2	2.7	0.1
A0576	5.25	0.6	2.9	0.3	15.01	0.2	3.1	0.1
A0754	3.8	1.5	5.6	2.3	19.61	0.4	7.1	2.8
A1060	2.63	1.7	4.5	0.4	13.74	0.6	7.7	0.1
A1367	3.68	1.8	6.4	5.6	18.64	0.4	7.7	5.0
A1644	2.77	1.2	3.3	0.0	7.67	0.5	4.1	0.0
A1795	0.58	2.3	1.3	0.0	5.03	0.3	1.6	0.0
A2065	4.99	0.8	4.2	1.5	26.23	0.2	4.5	1.5
A2142	0.57	3.0	1.7	0.0	2.61	0.8	2.1	0.0
A2199	1.84	2.4	4.4	1.9	7.09	0.7	4.9	1.2
A2244	2.23	0.6	1.4	0.0	6.86	0.2	1.4	0.0
A2255	5.67	0.8	4.5	5.1	11.27	0.4	4.5	4.4
A2256	0.49	2.4	1.2	0.0	2.13	0.7	1.6	0.0
A2589	5.06	0.5	2.7	0.0	21.48	0.2	3.8	0.0
A2597	1.97	0.5	1.0	0.0	12.69	0.1	1.1	0.0
A2634	6.84	0.6	3.8	0.4	15.34	0.3	4.2	0.3
A2657	3.01	0.7	2.1	0.0	27.51	0.2	4.1	0.1
A2734	3.18	0.4	1.2	0.0	7.4	0.2	1.3	0.0
A2877	2.14	0.4	0.9	0.0	16.41	0.1	1.2	0.0
A3112	7.57	0.8	6.0	13.8	29.98	0.2	6.1	13.4
A3158	1.51	1.1	1.6	0.0	3.78	0.5	1.7	0.0
A3266	1.34	2.9	3.9	2.1	6.17	0.7	4.6	1.7
A3376	8.43	0.6	5.4	0.5	28.36	0.3	7.5	1.2
A3822	5.27	0.4	2.1	0.0	9.07	0.2	2.2	0.0
A3827	0.72	0.9	0.6	0.0	3.91	0.2	0.7	0.0
A3921	5.96	0.4	2.5	0.6	19.58	0.1	2.7	0.6
A4038	2.49	1.0	2.4	0.1	11.23	0.3	3.6	0.5
A4059	2.66	0.7	2.0	0.1	11.15	0.2	2.1	0.0
COMA	0.41	10.4	4.3	0.9	1.37	4.9	6.8	2.2
EXO0422	0.92	0.6	0.5	0.0	4.83	0.1	0.6	0.0
FORNAX	7.19	0.7	4.7	0.5	74.75	0.2	12.3	1.8
HCG94	8.44	0.3	2.6	0.0	23.26	0.2	3.7	0.0
HYDRA-A	3.68	1.2	4.2	3.0	18.86	0.2	3.9	1.2
IIIZw54	12.94	0.4	4.9	0.4	52.7	0.1	4.2	0.0
IIZw108	3.66	0.4	1.3	0.0	6.43	0.2	1.4	0.0
NGC 1550	5.38	0.5	2.5	0.0	25.85	0.2	4.9	0.0
NGC 5044	9.95	0.3	3.1	0.0	61.79	0.1	6.1	0.1
RXJ2344	5.79	0.5	2.9	1.3	30.83	0.1	2.8	0.4
S405	3.5	0.4	1.3	0.0	6.16	0.2	1.5	0.0
S540	14.76	0.3	4.1	1.2	48.07	0.1	4.5	1.0
UGC03957	10.74	0.4	4.1	0.5	75.34	0.1	5.3	0.6
ZwCl1742	2.34	1.0	2.3	0.0	12.66	0.2	3.0	0.0

Notes. Individual flux predictions and upper limits for alternative CR profiles for photon energies above 500 MeV. Columns 2–5 refer to the individual scale factor, the flux prediction above 500 MeV ($F_{\gamma,\mathrm{exp}}^{\mathrm{CR}}$), the derived upper limit on the γ-ray flux, and the individual TS value for the ICM model, while Columns 6–9 represent the values obtained for the flat CR model. All photon fluxes are given in 10⁻¹⁰ photon s⁻¹ cm⁻².

Download table as:  ASCIITypeset image

We emphasize that these eight models do not span the complete uncertainty of the systematics involved with interstellar emission modeling. They do not even encompass the full uncertainty in the input parameters that are varied. The resulting uncertainty should therefore only be considered as one indicator of the systematic uncertainty due to interstellar emission modeling. The tests we performed using these additional models considered a different energy range, from 500 MeV–100 GeV, because the alternative models were not derived for higher energies. However, since events at low energies dominate the fit, this difference is negligible. We have explicitly verified this by repeating the baseline analysis up to 100 GeV and we found that the differences between the computed combined scale factors are <1%.

In Figure 12, we show how the combined upper limit on the scale factor, $A_{\gamma }^{{\rm UL}}$ varies for the alternative models with respect to our standard model as well as how the significance of the excesses observed in A1367 and A3112, and A400 changes when using the alternative models.

Zoom In Zoom Out Reset image size

The spread in limits for the scale factor is rather small between the alternative diffuse models, but the choice of using the standard diffuse model versus any of the alternative diffuse models can affect the inferred limits for the different cluster samples by 20%–30%. Comparing the TS values for the three clusters indicates small variations across the alternative models for A400, A3112, and A1367. We repeated this procedure for the derivation of the individual flux limits (see Section 5.4) and found that the majority of our clusters follow the same trend. We have marked the clusters which show variations beyond 30% in Table 6.

We summarize the systematic uncertainties discussed in this section in Table 9, and we note that the main source of uncertainty is the accurate modeling of the foreground Galactic diffuse emission.

The joint likelihood approach discussed in Section 3.1 makes the assumption that the individual data samples are uncorrelated. Assuming two sources s₁ and s₂, in the uncorrelated case, the composite p-value is p_Σ = p₁ × p₂, where p₁ and p₂ are the p-values associated with s₁ and s₂, respectively. This case corresponds to two sources that are far away from each another. In this case, the difference between p_joint, which is the derived p-values from the joint likelihood, should be minimal, while for small distances, correlations impact the derivation of p_joint, and thus lead to an overestimation of the significance associated with this p-value. We assess this bias using an isotropic simulation with two identical power-law sources with Γ = −2.3 that were each modeled as an extended source with a disk of 2° in diameter. In Figure 13, we show the difference between p_joint and $p_\Sigma$. We find that at small distances (<3°), for all minimum energy thresholds, there is a substantial bias due to overlaps. Toward larger distances, this bias is reduced. In this toy model, 3° corresponds to R₂₀₀ + 1°. However, by requiring larger distances between sources, we reduce the number of viable cluster candidates for our search. Hence we decided to use E_min = 500 MeV as this threshold maximizes the number of clusters to be included (and thus the expected signal) while minimizing the bias on p_joint.

Zoom In Zoom Out Reset image size

In Figure 14, we show ROI-specific TS-distributions for a background only simulation (top). We find that the TS-distribution in the background-only case can be well described by $({1}/{2}) \delta + ({1}/{2}) \chi _{k}^{2}$ and obtain k = 1.1 ± 0.1 in an unbinned maximum likelihood fit, giving rise to the usual definition of significance as $\sqrt{{\rm TS}}$.

Zoom In Zoom Out Reset image size

In addition, we include the results from simulations including a (weak) putative CR-induced γ-ray signal corresponding to A_γ = 0.29 (nominal best-fit value for the combined sample) as well as A_γ = 1 (middle panel in Figure 14), i.e., assuming the predictions by Pinzke et al. (2011). Finally, we repeat the analysis assuming a strong signal, characterized by A_γ = 10 (bottom panel in Figure 14). For all simulations, the true value of A_γ is recovered in the combined result, validating our analysis approach, although the global significance varies with respect to the signal simulation. For the background-only case, only upper limits can be derived. Assuming the nominal best-fit value for A_γ from the combined sample, yields a combined significance of 2.1σ. For A_γ = 1 and A_γ = 10, these values are much higher, yielding a global TS value corresponding to a significance of 13.3σ and 79.9σ, respectively.

Given that a simulation with A_γ = 0.29 yields a combined significance comparable with what we have found in our analysis, we investigated how this signal could be studied further. To that end, we compare the expected γ-ray flux based on the ROI-specific scale factors A_γ, i with their associated ROI TS values and quantify the correlation through calculating Spearman-rank correlation coefficient (Spearman 1904), denoted by r_s . We find Spearman-rank coefficients >0.5 indicating a correlation. However, even for the background case, r_s = 0.6 is obtained, which is the same as what we find in the signal case for A_γ = 1. This illustrates that a correlation analysis with such a weak signal, or further, studying the excess we find, is difficult. Only for the strong signal of A_γ = 10 do we find a correlation coefficient r_s = 0.94 indicating a strong correlation between expected γ-ray flux and associated ROI TS value.