MAGIC GAMMA-RAY TELESCOPE OBSERVATION OF THE PERSEUS CLUSTER OF GALAXIES: IMPLICATIONS FOR COSMIC RAYS, DARK MATTER, AND NGC 1275

In the course of this work, we used cosmological simulations of the formation of galaxy clusters to inform us about the expected spatial and spectral characteristics of the CR-induced γ-ray emission. A clear detection of the IC emission from shock-accelerated CR electrons will be challenging for Imaging Atmospheric Cherenkov Telescopes (IACTs) due to the large angular extent of these accretion shocks that subtend solid angles corresponding to up to six virial radii. For these instruments, the spatially concentrated pion-decay γ-ray emission resulting from hadronic CR interactions that dominates the total γ-ray luminosity (Pfrommer et al. 2008; Pfrommer 2008) should be more readily detectable than the emission from the outer region.

To address the question of universality and predictability of the expected γ-ray emission, we simulated a sample of 14 galaxy clusters that span one and a half decades in mass and show a variety of dynamical states ranging from relaxed cool core clusters to violent merging clusters (details are given in Section 4.1). In order to find the most promising target cluster in the local Universe for detecting the pion-decay emission, we computed the scaling relations between the γ-ray luminosity and cluster mass of our sample (Pfrommer 2008) and used these to normalize the CR-induced emission of all clusters in a complete sample of the X-ray brightest clusters (the extended HIFLUGCS catalog; Reiprich & Böhringer 2002). This favors a high-mass, nearby galaxy cluster with a scaling M^β₂₀₀/D²_lum, where M₂₀₀ is the virial mass,³⁴ D_lum is the luminosity distance, and β ≃ 1.32 is a weakly model-dependent scaling parameter that provides the rank ordering according to the brightness of each individual cluster (Pfrommer 2008). As a second criterion, we required low zenith angle observations, i.e., below 35°, that ensure the lowest possible energy thresholds and the maximum sensitivity for the detector. We carefully modeled the most promising targets, accounting for the measured gas density and temperatures from thermal X-ray measurements while assuming a constant CR-to-thermal gas ratio (Pfrommer & Enßlin 2004a). Cluster-wide extended radio synchrotron emission that informs about present high-energy processes was additionally taken into account before we selected the Perseus cluster as our most promising source. Although other clusters showed a somewhat higher γ-ray flux in our simulations (e.g., Ophiuchus), the facts that Perseus is observable at low zenith angles and that the expected emission is more spatially concentrated make it the best suited target for this observation.

Typically up to 80% of the total mass of a galaxy cluster is in the form of non-baryonic DM. Since the DM annihilation γ-ray signal is expected to be proportional to the integrated squared DM density along the line of sight (Evans et al. 2004; Bergström & Hooper 2006), it is obvious that galaxy clusters could be good candidates to look for DM as well. This is true despite the fact that they are located at much larger distances than other potential DM candidates, such as dwarf spheroidal galaxy satellites of the Milky Way or the Galactic Center. One obvious reason is the huge amount of DM hosted by clusters compared with the rest of candidates. Perseus, for example, is located ∼1000 times farther than Milky Way dwarfs, but it contains roughly six orders of magnitude more DM than the Willman 1 dwarf galaxy, one of the most promising DM candidates according to recent work (Strigari et al. 2007; Aliu et al. 2009a). Additionally, the presence of substructures could be of crucial importance. Substructures in clusters may significantly enhance the DM signal over the smooth halo, while we do not expect this to be of special relevance for dwarf galaxies since their outer regions are severely affected by tidal stripping (Pinzke et al. 2009; M. A. Sánchez-Conde et al. 2010, in preparation).

Essentially, the annihilation flux is proportional to the product of two parameters (see, e.g., Evans et al. 2004 for details): a first one that captures all the particle physics (DM particle mass, cross section, etc.), which we will label as f_SUSY, and a second one, J_astro, that accounts for all the astrophysical considerations (DM distribution, telescope point-spread function (PSF), etc.). The particle physics factor just acts as a normalization in the expected annihilation flux, so we can neglect it when performing a comparative study—as we are doing in this section. Concerning the astrophysical factor, the DM distribution is commonly modeled with radial density profiles of the form ρ(r) = ρ_s/[(r/r_s)^γ(1 + (r/r_s)^α)^(β−γ)/α], where ρ_s and r_s represent a characteristic density and a scale radius, respectively (Kravtsov et al. 1998). These density profiles are well motivated by high-resolution N-body cosmological simulations. Here we adopt the Navarro–Frenk–White (Navarro et al. 1997, hereafter NFW) DM density profile, with (α,β,γ) = (1, 3, 1). For an NFW profile, 90% of the DM annihilation flux comes from the region within r_s so that the corresponding integrated luminosity is proportional to r³_sρ²_s. We can derive r_s and ρ_s for Perseus, assuming M₂₀₀ = 7.7 × 10¹⁴M_☉ (as given in Table 1) and a concentration of ∼6 (as given by the Bullock et al. 2001 virial mass–concentration scaling relation). We obtain r_s = 0.384 Mpc and ρ_s = 1.06 × 10¹⁵ M_☉ Mpc⁻³, which translates into a total value of J_astro ∼ 1.4 × 10¹⁶ GeV² cm⁻⁵ for the scale radius region. In the case of Coma, although slightly (∼15%) more massive than Perseus, the fact that it is located significantly farther (101 Mpc) translates into a slightly lower annihilation flux. Virgo, only 17 Mpc away from us, gives a larger DM annihilation flux, but here the large extension of the region from which most of the annihilation flux is expected to come compared with Perseus (r_s ∼ 12 and r_s ∼ 03, respectively) could represent an obstacle from the observational point of view. Source extension is of special relevance for single-telescope IACTs, for which point-like sources (sources with an angular extension smaller than or similar to the telescope PSF) are more readily observable.

The central NGC 1275 radio galaxy is another strong motivation for γ-ray observations of the Perseus galaxy cluster. The detection at TeV energies of the radio galaxies M 87 (Aharonian et al. 2006) and Centaurus A (Aharonian et al. 2009c) has forced a substantial revision of the paradigm whereby very high energy (VHE) emission is a characteristic property of highly relativistic jets closely aligned with the line of sight, establishing radio galaxies as a new class of VHE γ-ray emitters. Note that NGC 1275 has various characteristics in common with Centaurus A which has also been interpreted as a possible source of ultra-high-energy CRs (Hardcastle et al. 2009).

The NGC 1275 radio galaxy is the brightest radio source in the northern sky. Its jet inclination angle seems to increase from 10°–20° at milliarcsecond scales up to 40°–60° at arcsecond scales (Dunn et al. 2006). Note that NGC 1275 was classified as a blazar by Angel & Stockman (1980) because of its optical polarization, and it has been seen to vary in the optical on time scales of a day (Geller et al. 1979). All these elements are promising from the point of view of the TeV detectability, since they suggest that the emission region is located at the base of the jet. In these conditions, in the scenario based on the structured jet model (Tavecchio & Ghisellini 2008), we expected VHE emission from the layer of the jet at a level detectable by MAGIC.

MAGIC observed the Perseus cluster for 33.4 hr during 2008 November and December, at zenith angles between 12° and 32°, which guarantees the lowest energy threshold. The observation was performed in the false-source tracking (wobble) mode (Fomin et al. 1994) pointing alternatively to two different sky directions, each at a 04 distance from the nominal target position.

The main background for Cherenkov telescopes is due to the hadronic CRs and the night sky background. Our standard analysis procedure is as follows (for a detailed description, see Albert et al. 2008c): data calibration and extraction of the number of photoelectrons per pixel are done (Albert et al. 2008a). This is followed by an image cleaning procedure using the amplitude and timing information of the calibrated signals. Particularly, the arrival times in pixels containing >6 photoelectrons (core pixels) are required to be within a time window of 4.5 ns and for pixels containing >3 photoelectrons (boundary pixels) within a time window of 1.5 ns from a neighboring core pixel. For the surviving pixels of each event, the shower parameters are reconstructed using the Hillas parameterization algorithm (Hillas 1985). Hadronic background suppression is achieved using a multivariate method called Random Forest (Breiman 2001; Albert et al. 2008b), which uses the Hillas parameters to define an estimator called hadronness (it runs from 0 for gammas to 1 for hadrons) by comparison with Monte Carlo (MC) γ-ray simulations. Moreover, the Random Forest method is used for the energy estimation of a reconstructed shower. The gamma/hadron (g/h) separation in the analysis was optimized on a sample of well-understood Crab Nebula data, which is commonly accepted as a standard reference source for VHE astronomy.

Part of the data has been rejected mainly due to the bad weather conditions during some observation days. The total data rejected amount to ∼27%, resulting in 24.4 hr effective observation time of very high data quality. Independent cross-checks were performed on the data giving compatible results.

Given the good data quality and the low zenith angles of observation, the analysis energy threshold turns out to be 80 GeV. Beyond this threshold, no significant excess of γ-rays above the background was detected in 24.4 hr of observation. In Figure 1, the α-plot for energies above 250 GeV, where the best integral sensitivity is obtained from a Crab Nebula data sample, is reported. The α parameter is defined as the angular distance between the shower image main axis and the line connecting the observed source position in the camera and the image barycenter. Background events are isotropic in nature and thus produce randomly oriented shower images. This results in a more or less smooth event distribution in the α-plot. The γ-ray events due to the source, on the other hand, are predominantly aligned to the observed position in the camera. For a detected source, this results in a significant excess of events at small α. A fiducial region α < 6° and a hadronness cut of 0.05 are chosen by optimizing the analysis on a Crab Nebula data sample.

Zoom In Zoom Out Reset image size

In Figure 2, the significance map for events above 150 GeV in the observed sky region is shown. The source-independent DISP method has been used. This implies the rise of the energy threshold from 80 GeV to around 150 GeV (see Domingo-Santamaria et al. 2005 for a detailed description). The significance distribution in the map is consistent with background fluctuations. In Figure 2, X-ray contours from the XMM-Newton observations (Churazov et al. 2003) are also shown.

Zoom In Zoom Out Reset image size

The significance was calculated according to Equation (17) of Li & Ma (1983), and upper limit estimation is performed using the Rolke method (Rolke et al. 2005). The upper limits in number of excess events are calculated with a confidence level of 95%. For the upper limit calculation, a systematic uncertainty of 30% in the energy estimation and effective area calculation is taken into account. Our systematic error budget is obtained by adding up the individual contributions in quadrature. The different sources of systematic uncertainties are mainly related to the differences between the real experimental conditions and the simulated ones (see Albert et al. 2008c for a detailed discussion on the systematic errors). The photon flux upper limit is finally reconstructed for a general γ-ray spectrum as described in Aliu et al. (2009a).

In Sections 4 and 5, we will discuss the implications of this observation for the CR and DM annihilation-induced γ-ray flux, respectively. Using the true density profile as obtained by X-ray measurements (Churazov et al. 2003), we will be able to model the spatial characteristics of the CR-induced γ-ray signal. Our simulations indicate that 60% of the total γ-ray flux are contained within a circle of radius r_0.6 = 015 (this angular scale corresponds to a physical radius of 200 kpc). We then compare the flux from within this region to the upper limits. As the characteristics of the considered emission region are close to a point source, we use point-like upper limits. The same conclusion is valid also for the DM annihilation signal. In this case, as explained in Section 2.2, 90% of the expected emission is coming from the scale radius region. For Perseus, we obtained r_s ∼ 03 which is somewhat extended compared to the telescope angular resolution. However, the fact that the NFW profile is very steep implies that the main DM emission comes from the core of the source that can be considered approximately point-like compared to our angular resolution.

To compute flux upper limits, we assume specific spectral indices that have been motivated by an astrophysical scenario (see the following sections). This "scenario-guided" approach allows us to provide the tightest limits on physically motivated parameters and underlying astrophysical models. In the next sections we will consider flux upper limits computed using a power-law γ-ray spectrum with spectral indices Γ of −1.5, −2.2, and −2.5. In Table 2, the corresponding integral flux upper limits for energies above 100 GeV are listed.

In Section 4, we will use an integral flux upper limit set above given energy thresholds in order to trace the energy range where we can better constrain the models. In Table 3, the obtained integral flux upper limits for Γ = −2.2 are shown. Note that we do not compute integral upper limits above 80 GeV (as we have not shown a cumulative α-plot for energies above this value). This is because the g/h separation for events below 100 GeV works in a substantially different way with respect to the higher energy events. Therefore, we analyze separately the events below 100 GeV and the events of higher energy, with different sets of analysis cuts.

Finally, for completeness, in Table 4 the differential flux upper limits for the assumed spectral indices are shown in different energy intervals. Spectral energy density (SED) upper limits can also be obtained from those differential flux upper limits, as shown in Section 6 discussing the observation implications for the radio galaxy NGC 1275.

There are few existing IACT observations of galaxy clusters (Perkins et al. 2006; Perkins 2008; Aharonian et al. 2009a, 2009b; Domainko et al. 2009; Galante et al. 2009; Kiuchi et al. 2009; Acciari et al. 2009). In Section 4.3, we will compare the limits on the CR-to-thermal pressure obtained by other IACTs with those derived in this work. However, there are two observations of the Perseus galaxy cluster made by WHIPPLE (Perkins et al. 2006) and VERITAS (Acciari et al. 2009) with which we can directly compare our upper limits.

The WHIPPLE Collaboration observed the Perseus galaxy cluster (Perkins et al. 2006) for ∼13 hr obtaining an integral upper limit above 400 GeV of 4.53 × 10⁻¹² cm⁻² s⁻¹ assuming a spectral index Γ = −2.1. We can compare this value with our integral upper limit above 400 GeV of 1.83 × 10⁻¹² cm⁻² s⁻¹ with Γ = −2.2 (see Table 3). Our upper limit is significantly lower than the WHIPPLE one; clearly, this is not a surprise as the MAGIC telescope belongs to a new generation of IACTs. More recently, the VERITAS Collaboration observed Perseus (Acciari et al. 2009) for ∼8 hr and obtained an integral upper limit above 126 GeV of 1.27 × 10⁻¹¹ cm⁻² s⁻¹ assuming Γ = −2.5. We can compare this value with our corresponding integral upper limit above 100 GeV of 7.52 × 10⁻¹² cm⁻² s⁻¹ (see Table 2). Despite the fact that the VERITAS sensitivity of about 1% of Crab Nebula (Otte et al. 2009) is better than the MAGIC one, our upper limit is slightly lower than that found by Acciari et al. (2009) as expected from the significant difference in observation time.

Simulations were performed using the "concordance" cosmological cold DM model with a cosmological constant (ΛCDM) motivated by First Year cosmological constraints of Wilkinson Microwave Anisotropy Probe (WMAP). The cosmological parameters of our model are Ω_m = Ω_dm + Ω_b = 0.3, Ω_b = 0.039, Ω_Λ = 0.7, h = 0.7, n = 1, and σ₈ = 0.9. Here, Ω_m denotes the total matter density in units of the critical density for geometrical closure today, ρ_crit(z = 0) = 3H²₀/(8πG). Ω_b and Ω_Λ denote the densities of baryons and the cosmological constant at the present day, respectively. The spectral index of the primordial power spectrum is denoted by n, and σ₈ is the rms linear mass fluctuation within a sphere of radius 8 h⁻¹Mpc extrapolated to z = 0.

Our simulations were carried out with an updated and extended version of the distributed-memory parallel TreeSPH code GADGET-2 (Springel et al. 2001; Springel 2005). Gravitational forces were computed using a combination of particle-mesh and tree algorithms. Hydrodynamic forces were computed with a variant of the smoothed particle hydrodynamics (SPH) algorithm that conserves energy and entropy where appropriate, i.e. outside of shocked regions (Springel & Hernquist 2002). We have performed high-resolution hydrodynamic simulations of a sample of galaxy clusters that span over one and a half decades in mass and show a variety of dynamical states ranging from relaxed cool core clusters to violent merging clusters. Our simulated clusters have originally been selected from a low-resolution DM-only simulation (Yoshida et al. 2001). Using the "zoomed initial conditions" technique (Katz & White 1993), the clusters have been re-simulated with higher mass and force resolution. In high-resolution regions, the DM particles had masses of m_dm = 1.61 × 10⁹ h⁻¹₇₀ M_☉ and SPH particles m_gas = 2.4 × 10⁸ h⁻¹₇₀ M_☉, so each individual cluster is resolved by 8 × 10⁴–4 × 10⁶ particles, depending on its final mass. The SPH densities were computed from 48 neighbors, allowing the SPH smoothing length to drop at most to half of the value of the gravitational softening length of the gas particles. This choice of the SPH smoothing length leads to our minimum gas resolution of approximately 1.1 × 10¹⁰ h⁻¹₇₀ M_☉. For the initial redshift, we chose 1 + z_init = 60. The gravitational force softening was of a spline form (e.g., Hernquist & Katz 1989) with a Plummer-equivalent softening length that is assumed to have a constant comoving scale down to z = 5 and a constant value of 7 h⁻¹₇₀ kpc in physical units at later epochs.

These simulations included radiative hydrodynamics, star formation, and supernova feedback and followed CR physics using a novel formulation that followed the most important injection and loss processes self-consistently while accounting for the CR pressure in the equations of motion (Pfrommer et al. 2006; Enßlin et al. 2007; Jubelgas et al. 2008). To obtain predictions of the GeV–TeV γ-ray emission from clusters, we used an updated version of the CR physics in our code. It is capable of following the spectral evolution of the CR distribution function by tracking multiple CR populations in each gaseous fluid element; each of these populations is described by an amplitude, a low-momentum cutoff, and a characteristic power-law distribution in particle momentum with a distinctive slope that is determined by the acceleration process at formation shocks or supernova remnants (A. Pinzke & C. Pfrommer 2010, in preparation). Adiabatic CR transport processes such as compression and rarefaction and a number of physical source and sink terms which modify the CR pressure of each particle are modeled. The most important sources considered are diffusive shock acceleration at cosmological structure formation shocks and optional injection by supernovae while the primary sinks are thermalization by Coulomb interactions and catastrophic losses by hadronic interactions. We note that the overall normalization of the CR distribution scales with the maximum acceleration efficiency at structure formation shock waves. Following recent observations at supernova remnants (Helder et al. 2009) as well as theoretical studies (Kang & Jones 2005), we adopt a realistic value of this parameter and assume that 50% of the dissipated energy at strong shocks is injected into CRs while this efficiency rapidly decreases for weaker shocks (Enßlin et al. 2007).

We computed the γ-ray emission signal and found that it obeys a universal spectrum and spatial distribution (A. Pinzke & C. Pfrommer 2010, in preparation). This is inherited from the universal concave spectrum of CRs in galaxy clusters that is caused by the functional form and redshift dependence of the Mach number distribution of structure formation shocks that are responsible for the acceleration of CRs (Pfrommer et al. 2006). The CR distribution has a spectral index of Γ ≃ −2.5 at GeV energies and experiences a flattening toward higher energies resulting in Γ ≃ −2.2 at energies above a few TeV. Hence, the resulting γ-ray spectrum from CR-induced pion decay shows a characteristic spectral index of Γ ≃ −2.2 in the energy regime ranging from 100 GeV to TeV. The spatial distribution of the CR number density is mainly governed by adiabatic transport processes (Pfrommer et al. 2007) and similarly attains an approximate universal shape relative to that of the gas density. These findings allow us to reliably model the CR signal from nearby galaxy clusters using their true density profiles as obtained by X-ray measurements that we map onto our simulated density profiles.

In addition to CR protons, we modeled relativistic electrons that have been accelerated at cosmological structure formation shocks (primary CR electrons) and those that have been produced in hadronic interactions of CRs with ambient gas protons (secondary CR electrons). Both populations of CR electrons contribute to the γ-ray emission through Compton up-scattering photons from the cosmic microwave background as well as the cumulative starlight from galaxies. It turns out that the pion-decay emission of the cluster dominates over the contribution from both IC components—in particular, for relaxed systems (Pfrommer 2008).

In our optimistic CR model (radiative physics with galaxies), we calculated the cluster total γ-ray flux within a given solid angle. In contrast, we cut the emission from individual galaxies and compact galactic-sized objects in our more conservative model (radiative physics without galaxies). In short, the ICM is a multiphase medium consisting of a hot phase which attained its entropy through structure formation shock waves dissipating gravitational energy associated with hierarchical clustering into thermal energy. The dense, cold phase consists of the true interstellar medium (ISM) within galaxies and at the cluster center as well as the ram-pressure-stripped ISM. These cold dense gas clumps dissociate incompletely in the ICM due to insufficient numerical resolution as well as so far incompletely understood physical properties of the cluster plasma. All of these phases contribute to the γ-ray emission from a cluster. To assess the bias associated with this issue, we performed our analysis with both limiting cases bracketing the realistic case.

In Figure 3, we compare the integral flux upper limits obtained in this work (see Table 3) with the simulated flux that is emitted within a circle of radius r_0.6 = 015 for our two models, with and without galaxies. The upper limits are a factor of 2 larger than our conservative model and a factor of 1.5 larger than our most optimistic model predictions implying consistency with our cosmological cluster simulations. We note however that our simulated flux represents a theoretical upper limit of the expected γ-ray flux from structure formation CRs; lowering the maximum acceleration efficiency would decrease the CR number density as well as the resulting γ-ray emission.

Zoom In Zoom Out Reset image size

In Figure 4, we show the simulated γ-ray surface brightness map of a cooling flow cluster of mass similar to Perseus. As the CR-induced γ-ray flux is a radially declining function, so is the CR pressure. A quantity that is of great theoretical interest is the CR pressure relative to the thermal pressure X_CR = P_CR/P_th, as it directly assesses the CR bias of hydrostatic cluster masses since the CR pressure enters in the equation of motion. On the right-hand side of Figure 4, we show the profile of the CR-to-thermal pressure (volume-weighted) of this simulated cluster. Moving from the periphery toward the center, this quantity is a steadily declining function until we approach the cooling flow region around the cD galaxy of this cluster (similar to NGC 1275) where the CR pressure rises dramatically relative to that of the thermal gas which cools on a short time scale (Pfrommer et al. 2006). The volume average is 〈X_CR〉 = 〈P_CR〉/〈P_th〉 = 0.02, dominated by the region around the virial radius, while the ratio of CR-to-thermal energy is given by E_CR/E_th = 0.032.³⁵ Perseus has a smaller mass and a corresponding temperature that is only half of that of our simulated cooling flow cluster. Noting that X_CR ∝ 1/P_th ∝ 1/kT,³⁶ we expect these values to be a factor of ≃2 larger in Perseus, yielding 〈X_CR〉 ≃ 0.04 for the entire cluster and 〈X_CR〉 ≃ 0.02 for the core region that we probe with the present observation.

Zoom In Zoom Out Reset image size

We have to scale our conservative model prediction by a factor of ∼2 to reach the upper limits (cf. Figure 3) which implies that this work constrains the relative pressure contained in CRs to <8% for the entire cluster and to <4% for the cluster core region. The presence of dense gas clumps potentially biases the simulated γ-ray flux high and hence the inferred limits on X_CR low. Another source of bias could be unresolved point sources inside the cluster such as an active galactic nucleus (AGN). In the presented simulation of the cool core cluster g51, the bias due to subclumps amounts to a factor of 1.5 but it could be as high as 2.4 which is the mean difference between our conservative and optimistic models across our scaling relations. We note however that the latter case is already excluded by our upper limits provided the maximum shock acceleration efficiency is indeed as high as 50%. While there are indications from supernova remnant observations of one rim region (Helder et al. 2009) as well as theoretical studies (Kang & Jones 2005) that support such high efficiencies, to date it is not clear whether these efficiencies apply in an average sense to strong collisionless shocks or whether they are realized for structure formation shocks at higher redshifts. Improving the sensitivity of the presented type of observations will help in answering these profound plasma astrophysics questions.

In Figure 4, we additionally compare a simulation where we only accelerate CRs at structure formation shocks with one where we additionally account for CRs that are injected through supernova feedback within the star-forming regions in our simulation. Outside the cD galaxy, there is no significant difference visible which suggests that the CRs injected into the ICM by supernova-driven winds are negligible compared with those accelerated by structure formation shocks. While this is partly an artifact of our simulations that neglect CR diffusion, we expect this behavior due to the adiabatic losses that CRs suffer as they expand from their compact galactic ISM into the dilute ICM. Assuming a conservative value for the density contrast of Δ = 10⁻³, the CR pressure is diluted by P_CR ∼ Δ^4/3 P_CR,ISM ∼ 10⁻⁴ P_CR,ISM.

As anticipated in Section 3.3, there are few existing IACT observations of galaxy clusters, some of which derived limits on the CR-to-thermal pressure contained in clusters, in particular the WHIPPLE observation of the Perseus cluster (Perkins et al. 2006) and the HESS observations of the Abell 85 (Aharonian et al. 2009a; Domainko et al. 2009) and Coma (Aharonian et al. 2009b) clusters. These works used simplifying assumptions about the spectral and spatial distribution of CRs. They typically assumed a single CR power-law distribution with a spectral index of Γ = −2.1 (that provides optimistic limits on the CR-to-thermal pressure) and assumed that the CR energy density is a constant fraction of the thermal energy density throughout the entire cluster. Based on these two assumptions, WHIPPLE and HESS found in Perseus and Abell 85 E_CR/E_th < 0.08, respectively, while HESS found E_CR/E_th < 0.2 in Coma.

To facilitate comparison with these earlier works, we repeated the data analysis with a spectral index Γ = −2.1 to obtain an integral upper limit $\mathcal {F}_{\mathrm{UL}} (>100\,\mbox{ GeV}) = 6.22\times 10^{-12}\,\mbox{ cm}^{-2}\,\mbox{ s}^{-1}$. Following the formula of Pfrommer & Enßlin (2004a), we compute the γ-ray flux of a CR population with Γ = −2.1 within a circular region of radius r_0.6 = 015 or equivalently 200 kpc. In our isobaric model of CRs, we assume that the CR pressure scales exactly as the thermal pressure and constrain E_CR/E_th < 0.053 which corresponds to an averaged relative pressure of 〈X_CR〉 = 〈P_CR〉/〈P_th〉 = 0.033. This would be the most stringent upper limit on the CR energy in a galaxy cluster.

In our adiabatic model of CRs, we account for the centrally enhanced CR number density due to adiabatic contraction during the formation of the cooling flow (Pfrommer & Enßlin 2004a). We assume that the CRp population scaled originally as the thermal population but was compressed adiabatically during the formation of the cooling flow without relaxing afterward (we adopted temperature and density profiles given by Churazov et al. 2003). In this model, we obtain an enhanced γ-ray flux level for virtually the same volume-averaged CR pressure or vice versa for a given flux limit; hence, we can put a tighter constraint on the averaged CR pressure. We constrain E_CR/E_th < 0.03 which corresponds to an averaged relative pressure of 〈X_CR〉 = 〈P_CR〉/〈P_th〉 = 0.019.

How can we reconcile these tighter limits with our simulation-based slightly weaker limit? We have to compare our simulated CR profile to a CR distribution that does not show any enhancement relative to the gas density. In the central region for r < 200 kpc, we derive an adiabatic compression factor of 1.7 that matches that in our simplified approach—suggesting that our simple adiabatic model captures the underlying physics quite realistically. Second, we have then to relate the pressure of a power-law spectrum with Γ = 2.1 to our simulated concave spectrum. Noting that the γ-rays at 100 GeV are produced by CR protons at ≃1 TeV, we normalize both spectra at 1 TeV and find that the simulated spectrum contains a larger pressure by a factor of 1.8. This factor brings the limit of our simplified adiabatic model into agreement with our simulation-based limit of the relative CR pressure 〈X_CR〉 < 4% for the cluster core region. Finally, since γ-ray observations are only sensitive to the cluster core regions (the emission is expected to peak in the center due to the high target gas densities), they cannot constrain the average CR-to-thermal pressure within the entire cluster. Hence, we have to use cosmological cluster simulations to address how much CR-to-thermal pressure could be additionally hidden in the peripheral cluster regions.

For clusters that host radio (mini-)halos, we are able to derive a minimum γ-ray flux in the hadronic model of CR interactions. The idea is based on the fact that a steady-state distribution of CR electrons loses all its energy to synchrotron radiation for strong magnetic fields (B ≫ B_CMB ≃ 3.2 μG) so that the ratio of γ-ray to synchrotron flux becomes independent of the spatial distribution of CRs and thermal gas (Pfrommer 2008). This can be easily seen by considering the pion-decay-induced γ-ray luminosity L_γ and the synchrotron luminosity L_ν of a steady-state distribution of CR electrons that has been generated by hadronic CR interactions:

$$\begin{equation} L_\gamma = A_\gamma \displaystyle \int {d}V\, n_{\rm CR}n_\mathrm{gas}, \end{equation} \tag{ 1 }$$

$$\begin{equation} L_\nu = A_\nu \displaystyle \int {d}V\, n_{\rm CR}n_\mathrm{gas} \frac{\displaystyle \varepsilon _B^{(\alpha _\nu +1)/2}}{\displaystyle \varepsilon _\mathrm{CMB}+\varepsilon _B} \end{equation} \tag{ 2 }$$

$$\begin{equation} \simeq A_\nu \displaystyle \int {d}V\, n_{\rm CR}n_\mathrm{gas}\qquad \mathrm{for}\quad \varepsilon _B \gg \varepsilon _\mathrm{CMB}. \end{equation} \tag{ 3 }$$

Here A_γ and A_ν are dimensional constants that depend on the hadronic physics of the interaction (Pfrommer et al. 2008; Pfrommer 2008) and α_ν ≃ 1 is the observed synchrotron spectral index. Hence, we can derive a minimum γ-ray flux in the hadronic model:

$$\begin{equation} \mathcal {F}_{\gamma,\mathrm{min}} = \frac{\displaystyle A_\gamma }{\displaystyle A_\nu }\frac{\displaystyle L_\nu }{\displaystyle 4\pi D_\mathrm{lum}^2}, \end{equation} \tag{ 4 }$$

where L_ν is the observed luminosity of the radio mini-halo and D_lum denotes the luminosity distance to the respective cluster. Lowering the magnetic field would require an increase in the energy density of CR electrons to reproduce the observed synchrotron luminosity and thus increase the associated γ-ray flux.

Using the values of Table 1, we obtain a minimum γ-ray flux in the hadronic model of the radio mini-halo of $\mathcal {F}_{\gamma, \mathrm{min}} ({>}100 \mbox{ GeV}) = 6\times 10^{-13}\, \mbox{cm}^{-2}\mbox{ s}^{-1}$, assuming a power-law CR distribution with Γ ≳ −2.3. This lower limit is independent of the spatial distribution of CRs and magnetic fields. We note that the spectral index is consistent with the radio data.³⁷ It turns out that the requirement of strong magnetic fields violates the energy conditions in clusters as it implies a magnetic energy density that is larger than the thermal energy density—in particular, at the peripheral cluster regions. The minimum γ-ray flux condition requires a constant (large) magnetic field strength throughout the cluster while the thermal energy density is decreasing by more than a factor of 100 from its central value. This would imply that the magnetic field eventually dominates the energy density at the virial regions—a behavior that is unstable as it is subject to Parker-like buoyancy instabilities. Additionally, such a configuration would be impossible to achieve in the first place as the magnetic energy density typically saturates at a fixed fraction of the turbulent energy density which itself is only a small fraction of the thermal energy density in clusters (Schuecker et al. 2004). Hence, these considerations call for lowering the assumed cluster magnetic fields which should strengthen the lower limits on the γ-ray flux considerably—however at the expense that these limits inherit a weak dependence on the spatial distribution of magnetic fields and CRs.

Estimates of magnetic fields from Faraday rotation measures (RMs) have undergone a revision in the last few years with more recent estimates typically in the order of a few μG with slightly higher values up to 10 μG in cooling flow clusters (Clarke 2004; Enßlin & Vogt 2006). For the Perseus radio mini-halo, Faraday RMs are available only on very small scales (Taylor et al. 2006), i.e. few tens of pc. RM estimates are of the order of ∼7000radm² leading to magnetic field values of ∼25 μG assuming that the Faraday screen is localized in the ICM. This, however, appears to be unlikely as variations of 10% in the RM are observed on pc scales (Taylor et al. 2002), while ICM magnetic fields are expected to be ordered on significantly larger scales of a few kpc (Taylor et al. 2006; Vogt & Enßlin 2005; Enßlin & Vogt 2006). Application of the classical minimum-energy argument to the Perseus radio mini-halo data leads to estimates for the central magnetic field strength of B₀ ≃ 7 μG or even B₀ ≃ 9 μG for the more appropriate hadronic minimum-energy argument (Pfrommer & Enßlin 2004b).

We select a cooling flow cluster of our sample that is morphologically similar to Perseus with a mass M₂₀₀ ≃ 10¹⁵ M_☉ (the simulated cluster g51 of Pfrommer et al. 2008). We adopt a conservative choice for the central magnetic field strength of ∼10 μG and parameterize the magnetic energy density in terms of the thermal energy density by ε_B ∝ ε^0.5_th which ensures ε_B < ε_th/3 in the entire cluster. This allows us to strengthen the physically motivated lower limit to $\mathcal {F}_{\gamma, \mathrm{phys.\,min}} ({>} 100 \mbox{ GeV}) = 8.5\times 10^{-13} \mbox{ cm}^{-2} \mbox{ s}^{-1}$ as shown by the dash-dotted line in Figure 3. In the hadronic model, this minimum γ-ray flux implies a minimum CR pressure relative to the thermal pressure. Figure 3 shows that the minimum flux $\mathcal {F}_{\gamma, \mathrm{phys.\,min}}$ is a factor of 3.6 lower than the simulated flux for Perseus in our conservative model. As seen in Section 4.2, this model corresponds to a relative CR pressure of 〈X_CR〉 = 〈P_CR〉/〈P_th〉 = 0.04 where the averages represent volume averages across the entire cluster. Hence we obtain a minimum relative CR pressure, 〈X_CR, min〉 = 〈P_CR, min〉/〈P_th〉/3.6 = 0.01. This minimum CR pressure corresponds to a minimum total CR energy of E_CR min = E_CR min/E_th × E_th = 1.6 〈X_CR, min〉 × E_th = 9 × 10⁶¹ erg where we integrated the temperature and density profiles from X-ray observations (Churazov et al. 2003) to obtain the total thermal energy of E_th = 5.7 × 10⁶³ erg. These considerations show the huge potential of combining future TeV γ-ray and radio observations in constraining physical models of the non-thermal cluster emission and of obtaining important insights into the average distribution of cluster magnetic fields.