Simulating the gamma-ray emission from galaxy clusters

Simulating the gamma-ray emission from galaxy clusters: a universal cosmic ray spectrum and spatial distribution


Entering a new era of high-energy -ray experiments, there is an exciting quest for the first detection of -ray emission from clusters of galaxies. To complement these observational efforts, we use high-resolution simulations of a broad sample of galaxy clusters, and follow self-consistent cosmic ray (CR) physics using an improved spectral description. We study CR proton spectra as well as the different contributions of the pion decay and inverse Compton emission to the total flux and present spectral index maps. We find a universal spectrum of the CR component in clusters with surprisingly little scatter across our cluster sample. When CR diffusion is neglected, the spatial CR distribution also shows approximate universality; it depends however on the cluster mass. This enables us to derive a semi-analytic model for both, the distribution of CRs as well as the pion-decay -ray emission and the secondary radio emission that results from hadronic CR interactions with ambient gas protons. In addition, we provide an analytic framework for the inverse Compton emission that is produced by shock-accelerated CR electrons and valid in the full -ray energy range. Combining the complete sample of the brightest X-ray clusters observed by ROSAT with our -ray scaling relations, we identify the brightest clusters for the -ray space telescope Fermi and current imaging air Čerenkov telescopes (MAGIC, HESS, VERITAS). We reproduce the result in Pfrommer (2008), but provide somewhat more conservative predictions for the fluxes in the energy regimes of Fermi and imaging air Čerenkov telescopes when accounting for the bias of ‘artificial galaxies’ in cosmological simulations. We find that it will be challenging to detect cluster -ray emission with Fermi after the second year but this mission has the potential of constraining interesting values of the shock acceleration efficiency after several years of surveying. Comparing the predicted emission from our semi-analytic model to that obtained by means of our scaling relations, we find that the -ray scaling relations underpredict, by up to an order of magnitude, the flux from cool core clusters.

magnetic fields, cosmic rays, radiation mechanisms: non-thermal, elementary particles, galaxies: cluster: general, Galaxy: fundamental parameters

1 Introduction

1.1 General background

In the cold dark matter (CDM) universe, large scale structure grows hierarchically through merging and accretion of smaller systems into larger ones, and clusters are the latest and most massive objects that had time to virialise. This process leads to collisionless shocks propagating through the intra-cluster medium (ICM), accelerating both protons and electrons to highly relativistic energies (Drury, 1983; Blandford & Eichler, 1987; Sarazin, 1999). High resolution X-ray observations by the Chandra and XMM-Newton satellites confirmed this picture, with most clusters displaying evidence for significant substructures, shocks, and contact discontinuities (e.g., Rosati et al., 2002; Voit, 2005; Markevitch & Vikhlinin, 2007). In addition, observations of radio halos and radio relics demonstrate the presence of synchrotron emitting electrons with energies reaching 10 GeV in more than 50 clusters (Feretti, 2003; Ferrari et al., 2008), although their precise origin in radio halos is still unclear. Similar populations of electrons may radiate -rays efficiently via inverse Compton (IC) upscattering of the cosmic microwave background photons giving rise to a fraction of the diffuse -ray background observed by EGRET (Loeb & Waxman, 2000; Totani & Kitayama, 2000; Miniati, 2002, 2003; Inoue & Nagashima, 2005). Although there is no clear observational evidence yet for a relativistic proton population in clusters of galaxies, these objects are expected to contain significant populations of relativistic protons originating from different sources, such as structure formation shocks, radio galaxies, and supernovae driven galactic winds. The ICM gas should provide ample target matter for inelastic collisions of relativistic protons leading to -rays (Völk et al., 1996; Enßlin et al., 1997; Miniati, 2003; Pfrommer & Enßlin, 2003, 2004; Pfrommer et al., 2008; Pfrommer, 2008; Kushnir & Waxman, 2009) as well as secondary electron injection (Dennison, 1980; Vestrand, 1982; Blasi & Colafrancesco, 1999; Dolag & Enßlin, 2000; Miniati et al., 2001b; Pfrommer & Enßlin, 2004; Pfrommer et al., 2008; Kushnir & Waxman, 2009; Kushnir et al., 2009). These hadronic collision processes should illuminate the presence of these elusive particles through pion production and successive decay into the following channels:

This reaction can only unveil those cosmic ray protons (CRs) which have a total energy that exceeds the kinematic threshold of the reaction of  GeV. The magnetic fields play another crucial role by confining non-thermal protons within the cluster volume for longer than a Hubble time, i.e. any protons injected into the ICM accumulates throughout the cluster’s history (Völk et al., 1996; Enßlin et al., 1997; Berezinsky et al., 1997). Hence, CRs can diffuse away from the production site, establishing a smooth distribution throughout the entire ICM which serves as efficient energy reservoir for these non-gravitational processes (Blasi & Colafrancesco, 1999; Dolag & Enßlin, 2000; Miniati et al., 2001a).

There is only little known theoretically about the spectral shape of the CR population in the ICM. It is an interesting question whether it correlates with injection processes or is significantly modified by transport and re-acceleration processes of CRs through interactions with magneto-hydrodynamic (MHD) waves. The most important processes shaping the CR spectrum as a function of cluster radius are (1) acceleration by structure formation shock waves (Quilis et al., 1998; Miniati et al., 2000; Pfrommer et al., 2006), MHD turbulence, supernova driven galactic winds (Acciari et al., 2009), or active galactic nuclei (AGN), (2) adiabatic and non-adiabatic transport processes, in particular anisotropic diffusion, and (3) loss processes such as CR thermalization by Coulomb interactions with ambient electrons and catastrophic losses by hadronic interactions. The spectral distribution of CRs that are accelerated at structure formation shocks should be largely described by a power-law with a spectral index of the one-dimensional distribution given by


where is the shock compression factor. Strong (high Mach number) shocks that inject a hard CR population occur either at high redshift during the formation of the proto-clusters or today at the boundary where matter collapses from voids onto filaments or super-cluster regions. In contrast, merger shocks show weak to intermediate strength with typical Mach numbers in the range of (Ryu et al., 2003; Pfrommer et al., 2006; Skillman et al., 2008). AGNs or supernova remnants are expected to inject CRs with rather flat spectra, (Völk et al., 1996; Schlickeiser, 2002; Enßlin, 2003), but it is not clear whether they are able to build up a homogeneous population of significant strength.

The CRs offer a unique window to probe the process of structure formation due to its long cooling times. While the thermal plasma quickly dissipates and erases the information about its past history, the CR distribution keeps the fossil record of violent structure formation which manifests itself through the spectrum that is shaped by acceleration and transport processes. The cluster -ray emission is crucial in this respect as it potentially provides the unique and unambiguous evidence of a CR population in clusters through observing the pion bump in the -ray spectrum. This knowledge enables determining the CR pressure and whether secondary electrons could contribute to the radio halo emission. In the -ray regime, there are two main observables, the morphological appearance of the emission and the spectrum as a function of position relative to the cluster center. The morphology of the pion induced -ray emission should follow that seen in thermal X-rays albeit with a slightly larger extent (Pfrommer et al., 2008). The primary electrons that are accelerated directly at the structure formation shocks should be visible as an irregular shaped IC morphology, most pronounced in the cluster periphery (Miniati, 2003; Pfrommer et al., 2008).

1.2 The -ray spectrum of a galaxy cluster

Figure 1: Upper Panel: CR spectral distribution within the virial radius of a Coma-like cluster. It shows three distinct features: a cutoff around GeV energies, a concave shape, and a steepening at high energies due to diffusive losses. Lower Panel: we show the intrinsic -ray number flux weighted by the photon energy that does not take into account photon propagation effects. The arrows indicate the spectral mapping from the CR spectrum to the photon spectrum. The pion decay flux is denoted in blue color, the secondary inverse Compton (sIC) in red, and the primary inverse Compton (pIC) emission in green. Due to the large uncertainty in the diffusion coefficient , we demonstrate how varying by a factor of two changes the corresponding -ray spectrum at high energies (dotted lines). The green band shows how the pIC emission changes if we vary the maximum electron injection efficiency from 0.05 (top) to 0.01 (bottom).

How do the spectral electron and proton distributions map onto the -ray spectrum? We show the CR spectrum within the virial radius of a simulated Coma-like galaxy cluster in the upper part of Fig. 1. It is shaped by diffusive shock acceleration at structure formation shocks, adiabatic transport and the relevant CR loss processes4. Three distinct features are visible in the spectrum: a cutoff close to the proton rest mass at  GeV, a concave shape for proton energies above and a steepening due to diffusive losses at energies , where is the value of the diffusion coefficient at 1 GeV. The dotted lines represents different values of the diffusion coefficient which is varied by a factor two from its fiducial value. The low energy cutoff is due to a balance of Coulomb and hadronic losses at energies around a GeV (Jubelgas et al., 2008). As shown in the present paper and in an upcoming work by Pinzke & Pfrommer (in prep.), the concave curvature is a unique shape5 that is caused by the cosmic Mach number distribution in combination with adiabatic transport processes. These features are mapped onto the pion decay -ray emission spectra as a consequence of hadronic CR interactions.

This can be seen in the lower part of Fig. 1, where the arrows indicate the spectral mapping from the CR spectrum to the photon spectrum. In a hadronic interaction, CRs produce pions that decay into photons with an energy that is on average smaller by a factor eight compare to the original CR energy (see Section 2.2). At CR energies that are larger than the hadronic reaction threshold, the CR power-law behavior is linearly mapped onto the pion decay induced -ray spectrum (solid blue). This emission component clearly dominates the total photon spectrum and therefore shapes the total emission characteristics in the central parts of the cluster, where the densities are high. Note that this spectrum is an intrinsic spectrum emitted at the cluster position and converted to a flux while assuming a distance of 100 Mpc without taking into account photon propagation effects. Depending on the cluster redshift, the finite mean free path of high-energy -rays to -pair production on infra red (IR) and optical photons limits the observable part of the spectrum to energies  TeV for clusters with redshifts , and smaller energies for higher redshift objects (Franceschini et al., 2008).

Secondary CR electrons and positrons up-scatter cosmic microwave background (CMB) photons through the IC process into the -ray regime, the so-called secondary inverse Compton emission (sIC). This emission component originates from the flat high-energy part of the CR spectrum and produces a rather flat sIC spectrum up to the Klein-Nishina (KN) regime. At large electron energies, we enter the KN regime of IC scattering where the electron recoil effect has to be taken into account. It implies less efficient energy transfer in such an elastic scattering event compared to the Thomson regime and leads to a dramatic steepening of the sIC spectrum at -ray energies around 100 TeV (solid red line). The dash-dotted red line shows the hypothetical sIC spectrum in the absence of the KN effect (which is never realized in Nature). However, it clearly shows that the diffusive CR break is not observable in the sIC component for large clusters (while it can move to energies below the KN break for small enough clusters, causing a faster steepening there). The spectrum shown in green color represents the energy weighted photon spectrum resulting from the IC process due to electrons accelerated at structure formation and merger shocks, the primary inverse Compton emission (pIC). The exponential cutoff is due to synchrotron and IC losses which lead to a maximum energy of the shock-accelerated electrons. The green pIC band shows the effect of the maximum electron injection efficiency, where we use an optimistic value of (see e.g. Keshet et al. 2003) in the top and a value of at the bottom. This more realistic value is suggested to be the theoretically allowed upper limit for the injection efficiency that is consistent with the non-thermal radiation of young supernova remnants (Zirakashvili & Aharonian, 2010).

This work studies the spectral and morphological emission characteristics of the different CR populations in the -ray regime. We concentrate on observationally motivated high-energy -ray bands. (1) The energy regime accessible to the Fermi -ray space telescope with a particular focus on and (2) the energy regime accessible to imaging air Čerenkov telescopes (IACTs) assuming a lower energy limit of . In Section 2 we describe the setup of our simulations, explain our methodology and relevant radiative processes considered in this work. In Section 3, we study emission profiles and maps, as well as spectral index maps. We then present the CR spectrum and spatial distribution and show its universality across our simulated cluster sample in Section 4. This allows us to derive a semi-analytic framework for the cluster -ray emission in Section 5 which we demonstrate on the Perseus and Coma galaxy clusters. Furthermore, we study the mass-to-luminosity scaling relations (Section 6) and predict the -ray flux from a large sample of galaxy clusters for the GeV and TeV energy regimes in Section 7. We compare our work to previous papers in this field and point out limitations of our approach in Section 8. We conclude our findings in Section 9. Throughout this work we use a Hubble constant of , which is a compromise between the value found by the Hubble key project (, Freedman et al., 2001) and from that one inferred from baryonic acoustic oscillation measurements (, Percival et al., 2009).

2 Setup and formalism

We follow the CR proton pressure dynamically in our simulations while taking into account all relevant CR injection and loss terms in the ICM, except for a possible proton production from AGN and supernova remnants. In contrast, we model the CR electron population in a post-processing step because it does not modify the hydrodynamics owing to its negligible pressure contribution. We use a novel CR formalism that allows us to study the spectral properties of the CR population more accurately.

2.1 Adopted cosmology and cluster sample

The simulations were performed in a CDM universe using the cosmological parameters: , and . Here, denotes the total matter density in units of the critical density for geometrical closure today, . , and denote the densities of baryons, dark matter, and the cosmological constant at the present day. The Hubble constant at the present day is parametrized as , while denotes the spectral index of the primordial power-spectrum, and is the rms linear mass fluctuation within a sphere of radius Mpc extrapolated to .

Our simulations were carried out with an updated and extended version of the distributed-memory parallel TreeSPH code GADGET-2 (Springel, 2005; Springel et al., 2001). Gravitational forces were computed using a combination of particle-mesh and tree algorithms. Hydrodynamic forces are 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). Our simulations follow the radiative cooling of the gas, star formation, supernova feedback, and a photo-ionizing background (details can be found in Pfrommer et al., 2007).

The clusters have originally been selected from a low-resolution dark-matter-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 by adding short-wavelength modes within the Lagrangian regions in the initial conditions that will evolve later-on into the structures of interest. We analyzed the clusters with a halo-finder based on spherical overdensity followed by a merger tree analysis in order to get the mass accretion history of the main progenitor. The spherical overdensity definition of the virial mass of the cluster is given by the material lying within a sphere centered on a local density maximum, whose radial extend is defined by the enclosed threshold density condition . We chose the threshold density to be a constant multiple of the critical density of the universe . In the remaining of the paper, we use the terminology instead of . Our sample of simulated galaxy clusters consists of 14 clusters that span a mass range from to where the dynamical stages range from relaxed cool core clusters to violent merging clusters (cf. Table 1). Each individual cluster is resolved by to particles, depending on its final mass. The SPH densities were computed from the closest 48 neighbors, with a minimum smoothing length () set to half the softening length. The Plummer equivalent softening length is in physical units after , implying a minimum gas resolution of approximately (see also Pfrommer et al., 2007).

Cluster sim.’s dyn. state
[] [Mpc] [keV]
1 g8a CC 2.9 13.1
2 g1a CC 2.5 10.6
3 g72a PostM 2.4 9.4
4 g51 CC 2.4 9.4
5 g1b M 1.7 4.7
6 g72b M 1.2 2.4
7 g1c M 1.2 2.3
8 g8b M 1.1 1.9
9 g1d M 1.0 1.7
10 g676 CC 1.0 1.7
11 g914 CC 1.0 1.6
12 g1e M 0.93 1.3
13 g8c M 0.91 1.3
14 g8d PreM 0.88 1.2

(1) The dynamical state has been classified through a combined criterion invoking a merger tree study and the visual inspection of the X-ray brightness maps. The labels for the clusters are M–merger, PostM–post merger (slightly elongated X-ray contours, weak cool core (CC) region developing), PreM–pre-merger (sub-cluster already within the virial radius), CC–cool core cluster with extended cooling region (smooth X-ray profile). (2) The virial mass and radius are related by , where denotes a multiple of the critical overdensity . (3) The virial temperature is defined by , where denotes the mean molecular weight.

Table 1: Cluster sample.

2.2 Modeling of CR protons and induced radiative processes

Our simulations follow cosmic ray physics in a self-consistent way (Pfrommer et al., 2006; Enßlin et al., 2007; Jubelgas et al., 2008). We model the adiabatic CR transport process such as compression and rarefaction, and a number of physical source and sink terms which modify the cosmic ray pressure of each CR population separately. The most important source considered6 for acceleration is diffusive shock acceleration at cosmological structure formation shocks, while the primary sinks are thermalization by Coulomb interactions, and catastrophic losses by hadronization. Collisionless structure formation shocks are able to accelerate ions and electrons in the high-energy tail of their Maxwellian distribution functions through diffusive shock acceleration (for reviews see Drury, 1983; Blandford & Eichler, 1987; Malkov & O’C Drury, 2001). In the test particle picture, this process injects a CR distribution with a power-law in momentum and a slope that depends on the instantaneous sonic Mach number of the shock. The overall normalization of the injected CR distribution depends on the adopted sub-resolution model of diffusive shock acceleration (e.g., Enßlin et al., 2007); in particular it depends on the maximum acceleration efficiency which is the maximum ratio of CR energy density that can be injected relative to the total dissipated energy density at the shock. We assume that in the saturated regime of shock acceleration, 50 percent of the dissipated energy at strong shocks is injected into cosmic ray protons. 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. This high efficiency rapidly decreases for weaker shocks (decreasing Mach number) and eventually smoothly approaches zero for sonic waves (Enßlin et al., 2007). Our paper aims at providing a quantitative prediction of the -ray flux and hence the associated CR flux that we expect in a cluster depending on our adopted acceleration model. Non-detection of our predicted emission will limit the CR acceleration efficiency and help in answering these profound plasma astrophysics questions about particle acceleration efficiencies.

We significantly revised the CR methodology and allow for multiple non-thermal cosmic ray populations of every fluid element (Pinzke & Pfrommer, in prep.). Each CR population is a power-law in particle momentum7,


characterized by a fixed slope , a low-momentum cutoff , and an amplitude that is a function of the position of each SPH particle through the variable . For this paper we have chosen five CR populations with the spectral index distribution for each fluid element (a convergence study on the number of CR populations is presented in the appendix A.2). This approach allows a more accurate spectral description8 as the superposition of power-law spectra enables a concave curvature of the composite spectrum in logarithmic representation. Physically, more complicated spectral features such as bumps can arise from the finite lifetime and length scale of the process of diffusive shock acceleration or incomplete confinement of CRs to the acceleration region. These effects imprint an upper cutoff to the CR population locally that might vary spatially and which translates into a convex curvature in projection. Additionally, interactions of pre-existing CRs with MHD waves can yield to more complex spectral features. Future work will be dedicated to study these topics.

In addition to the spectral features mentioned above, we model in the post-processing the effect of high-energy CR protons that are no longer confined to a galaxy cluster as these are able to diffuse into the ambient warm-hot intergalactic medium (WHIM). In this paper we define WHIM to be the region within , which is a subset of the entire WHIM (Davé et al., 2001). Assuming particle scattering off magnetic irregularities with the Kolmogorov spectrum, we obtain the characteristic scaling of the diffusion coefficient , where we normalize at . One can estimate the characteristic proton energy at which the spectrum steepens (Völk et al., 1996; Berezinsky et al., 1997),


For the reminder, we adopt a value of the diffusivity that is scaled to for each cluster, as this volume is expected to fall within the virialised part of the cluster past the accretion shock region (Pfrommer et al., 2007; Molnar et al., 2009) and traps CRs in a cluster for time scales longer than a Hubble time. This choice also has the property that the diffusion break is at energies  GeV; hence it does not interfere with the pion decay as well as secondary IC emission in the energy regime accessible to IACTs as we will show in the following. The momentum of a photon that results from pion decay is given by


This approximate relation is derived using the inelasticity and multiplicity for the channel together with the two photons in the final state. Secondary electrons that are injected in hadronic CR interactions Compton up-scatter CMB photons. A break in the parent CR spectrum would imprint itself in the sIC spectrum if there are no other effects that modify the spectrum at lower energies. Compared to the pion decay emission, this break manifests at slightly higher energies (for parameters adopted in Fig. 1). The momentum of the electrons depends on the proton momentum through the relation given by hadronic physics


Here we used the channel together with the four particles in the final state of the charged pion decay (, , , ). Combining the classical inverse Compton formulae from CR electrons with energies


with the energy relation in equation (5) we obtain a break in the sIC spectrum. This steepening of the CR spectrum take place at high photon energies where we choose CMB photons with the energy  meV as source for the inverse Compton emission using Wien’s displacement law. It turns out that these energies are deeply in the Klein-Nishina regime. This means that in the rest frame of the energetic electron, the Lorentz boosted photon energy is comparable to or larger than the electron rest mass, , so that the scattering event becomes elastic. This implies a less efficient energy transfer to the photon and manifests itself in a break in the resulting IC spectrum. While the number flux scales as in the Thomson-limit for  TeV, it steepens significantly to in the extreme KN-limit for  TeV, where is the spectral index of the (cooled) CR electron distribution (Blumenthal & Gould, 1970).

2.3 Magnetic fields

High energy CR electrons with loose their energy by means of IC scattering off CMB photons as well as through interactions with cluster magnetic fields which results in synchrotron emission. The relative importance of these two emission mechanisms depends on the rms magnetic field strength, , relative to the equivalent field strength of the CMB, , where denotes the redshift. In the peripheral cluster regions, where , the CR electrons loose virtually all their energy by means of IC emission. In the central cluster regions, in particular in the dense centers of cool cores, the magnetic energy density is probably comparable or even larger than the energy density of the CMB (Vogt & Enßlin, 2005), . Hence in these regions, the radio synchrotron emission carries away a fraction of the CR electrons’ energy losses; an effect that reduces the level of IC emission. We model the strength and morphology of the magnetic fields in the post-processing (Pfrommer et al., 2008) and scale the magnetic energy density field by the thermal energy density through the relation


where and denote the core values. If not mentioned otherwise, we use the magnetic decline and the central magnetic field throughout this paper. The central thermal energy density , is calculated by fitting the modified -model


to the radial pressure . The parametrization in equation (7) is motivated by both cosmological MHD SPH simulations (Dolag et al., 1999) and radiative adaptive mesh refinement MHD simulations (Dubois & Teyssier, 2008). Rather than applying a densities scaling as those simulations suggest, we use a scaling with thermal gas energy density which is not affected by the over-cooled centers in radiative simulations that do not take into account AGN feedback.

2.4 CR electron acceleration and inverse Compton emission

Modeling diffusive shock acceleration

Collisionless cluster shocks are able to accelerate ions and electrons through diffusive shock acceleration (for reviews see Drury, 1983; Blandford & Eichler, 1987; Malkov & O’C Drury, 2001). Neglecting non-linear shock acceleration and cosmic ray modified shock structure, the process of diffusive shock acceleration uniquely determines the spectrum of the freshly injected relativistic electron population in the post-shock region that cools and finally diminishes as a result of loss processes. The -ray inverse Compton emitting electron population cools on such a short time scale  yrs (compared to the long dynamical time scale  Gyr) that we can describe this by instantaneous cooling.9 In this approximation, there is no steady-state electron population and we would have to convert the energy from the electrons to inverse Compton and synchrotron radiation. Instead, we introduce a virtual electron population that lives in the SPH-broadened shock volume only; this is defined to be the volume where energy dissipation takes place. Within this volume, which is co-moving with the shock, we can use the steady-state solution for the distribution function of relativistic electrons and we assume no relativistic electrons in the post-shock volume, where no energy dissipation occurs. Thus, the cooled CR electron equilibrium spectrum can be derived from balancing the shock injection with the IC/synchrotron cooling: above a GeV it is given by


Here, we introduced the dimensionless electron momentum , where is the electron momentum, is the spectral index of the equilibrium electron spectrum, is the gas density, is the magnetic energy density, and denotes the photon energy density, taken to be that of CMB photons. The primary CRe distribution in equation (9) is calculated in the post-processing with a spectrum reflecting the current Mach number of the shock (without the assumption of spectral bins). Superposing the individual spectra of a large number of SPH particles, each with a spectrum reflecting the accelerating shock, produces a well defined total spectrum with a running index in general. A more detailed discussion of this simplified approach can be found in Pfrommer et al. (2008).

Once the radiative cooling time due IC and synchrotron emission becomes comparable to the diffusive acceleration time scale, the injection spectrum experiences a high-energy cutoff (Webb G.M., 1984). The electrons start to pile-up at this critical energy; the super-exponential term describing the maximum energy of electrons reached in this process, however, effectively cancels this pile-up feature which results in a prolonged power-law up to the electron cutoff momentum (Zirakashvili & Aharonian, 2007). We account for this effect by using the following parametrization of the shock injected electron spectrum,


where is the distance from the shock surface, and describe the characteristic momentum and shape of the pile-up region. The continuous losses cause the cutoff to move to lower energies as the electrons are transported advectively with the flow downstream. Integration over the post-shock volume causes the cutoffs to add up to a new power-law that is steeper by unity compared to the injection power-law (equation 9). Hence, the shock-integrated distribution function – as defined in equation (60) and displayed in Fig. 22 – shows three regimes: (1) at low energies (but large enough in order not to be affected by Coulomb losses) the original injected power-law spectrum, (2) followed by the cooled power-law that is steeper by unity, , and (3) an ultimate cutoff that is determined from the magnetic field strength and the properties of the diffusion of the electron in the shock (we refer the reader to Section B.2 for a detailed discussion). Note that only the last two regimes are important for -ray IC emission.

The fact that we observe X-ray synchrotron emission at shocks of young supernova remnants (Vink et al., 2006; Slane et al., 1999) necessarily requires the existence of high-energy CR electrons with  TeV. The non-thermal synchrotron emission generated by CR electrons with energies is given by


where (TeV) and magnetic fields of order G are required to reach X-ray energies of order 10 keV. To keep the highly relativistic electrons from being advected downstream requires efficient diffusion so that they can diffuse back upstream crossing the shock front again. We therefore use the most effective diffusion, refereed to as Bohm diffusion limit, as the electron propagation model at the shock. Balancing Bohm diffusion with synchrotron/IC cooling of electrons enables us to derive a maximum energy of the accelerated CR electrons at the position of the shock surface (derived in appendix, see also Webb G.M., 1984; Enßlin et al., 1998)


where diffusion parallel to the magnetic field has been assumed. Here denotes the flow velocity in the inertial frame of the shock (), and the electron loss time scale due to synchrotron and inverse-Compton losses reads


where is the Thompson cross-section. denotes the shock compression and is given by


where denotes the sonic Mach number. Inserting typical numbers for different cluster regions show that the electron cutoff energy, , only varies within a factor two inside (Table 2). The equivalent cutoff energy in the IC spectrum can easily be derived from equation (6), yielding .

Cluster variables WHIM
60 120 240
4 0.4 0.04 0.004
1.3 2.3 3.0 3.7
1.2 2 3 6
15 9 5 0.2
2400 3100 3500 1400
10 2.5 0.6 0.04
50 100 65 6.5

(1) Other constants used: mean molecular weight .

Table 2: Electron energy cutoff in different regions of a typical cluster.

IC emission

Following Rybicki & Lightman (1979), we calculate the inverse Compton emission from electrons that up-scatter CMB photons. It should be noted that we neglect the inverse Compton emission induced by starlight and dust, which might contribute significantly in the inner 10 kpc of the cluster. The inverse Compton source density in units of produced photons per unit time interval and volume for a simple power-law spectrum of CRes (equation 9) scales as


where . When we account for the competition between radiative cooling and diffusive acceleration of electrons in the shock region, the shape of in the high-energy regime changes. Using the cooled electron distribution of equation (60), we construct an effective integrated source function for primary inverse Compton emission (see Section B.3 for a self-consistent and extensive description),


where Bohm diffusion has been assumed. The normalization constants and are derived in equations (67) and (66), respectively. The KN suppression of the IC spectrum is captured by (equation LABEL:eq:KN_fit), and the shape of the transition region from the Thompson to the KN regime is given by (equation 71). Following recent work that carefully models the non-thermal radiation of young supernova remnants (Zirakashvili & Aharonian, 2010), we typically adopt a maximum electron injection efficiency of . We note that this value seems to be at the upper envelope of theoretically allowed values that match the supernova data.

2.5 Multiphase structure of the ICM

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 that has not yet dissociated into the ICM (Dolag et al., 2009). All of these phases contribute to the -ray emission from a cluster. Physically, the stripped ISM should dissociate after a time scale that depends on many unknowns such as details of magnetic draping of ICM fields (Dursi & Pfrommer, 2008; Pfrommer & Dursi, 2009) on galaxies or the viscosity of the ICM. In SPH simulations, this dissociation process is suppressed or happens only incompletely in our simulations leaving compact galactic-sized point sources that potentially biases the total -ray luminosity high. On the other hand, once these stripped compact point sources dissociate, the CRs diffuse out in the bulk of the ICM, and produce -rays by interacting with protons of the hot dilute phase. This flux, however, is negligible since


Here () denotes the gas number density in the ICM (ISM), the CR number density in the ISM, and the CR number density of the CRs that diffused out of their dense ISM environment into the ambient ICM that is in pressure equilibrium with the ISM. In the second step of equation (17) we accounted for the adiabatic expansion that these CRs would experience as they diffused out. In the last step we assumed , i.e. that the CR luminosity of all compact galactic sources in a cluster is of the same order as the CR luminosity in the diffuse ICM; a property that is at least approximately true in our simulations as we will show later on. Leaving all gaseous point sources would definitively be too optimistic, removing all of them would be too conservative since cluster spiral galaxies should contribute to the total -ray emission (which defines our so-called optimistic and conservative models). Hence, we perform our analysis with both limiting cases, bracketing the realistic case. The effect from the gaseous point sources is largest in low mass clusters, where they constitute a few percent of the total ICM mass. For high-mass clusters this fraction is lower, and constitutes only about one percent. In practice, we cut multiphase particles with an electron fraction and a gas density above the star forming threshold . If nothing else is stated, we use our conservative model without galaxies throughout the paper.

3 Characteristics of -ray emission

From surface brightness maps that are obtained by line-of-sight integration of the source functions, we study the -ray emission to characterize the morphology of clusters. Additionally, we use emission profiles to compare pIC, sIC, and pion decay induced emission for different clusters.

Figure 2: The -ray emission above of our Coma-like cluster g72a is shown. We show the pion decay -ray emission that originates from hadronic CR interactions with ambient gas protons on the left. On the right, we plot the inverse Compton emission from both, primary and secondary CR electrons. Comparing the different -ray emission components, we note that the pion decay has a very regular morphology and clearly dominates the cluster region. In contrast, the emission from primary electrons shows an irregular morphology that traces the structure formation shock waves and dominates in the virial periphery and the warm-hot intergalactic medium.
Figure 3: Projected energy dependent photon index for the galaxy cluster g72a. On the left, we show the photon index between and on the right, the photon index between . Note that the central part (where pion decay dominates) shows little variations which reflects the spectral regularity of the CR distribution. At the periphery and beyond, there are larger fluctuations due to the inhomogeneity of the pIC emission.

3.1 Morphology of -ray emission

The left side of Fig. 2 shows the morphology of the -ray emission above that results from hadronic CR interactions with ambient gas protons. The right side shows the primary and secondary IC emission for the post-merger cluster g72a. The comparison of the two panels shows that the central parts are dominated by the pion induced -ray emission. It has a very regular morphology that traces the gas distribution. There is a transition to the pIC emission as the dominant emission mechanism outside the cluster in the WHIM at a level depending on the dynamical state of the cluster. The pIC emission is very inhomogeneous which can be easily understood since it derives from primary CR electrons that are directly accelerated at structure formation shocks. Structure formation is not a steady process, it rather occurs intermittently. The morphology of the -ray emission above for g72a was investigated in Pfrommer et al. (2008). It shows a very similar morphology, indicating a similar power-law CR spectrum.

Figure 4: Azimuthally averaged profiles of the -ray emission components above for two different clusters: large cool core (CC) g8a (left panel) and large merger g72a (right panel). Shown are the dominating pion decay induced emission (dash-dotted lines), the IC emission of primary CR electrons accelerated directly at structure formation shocks (dotted lines), the secondary IC electrons resulting from CRp-p interactions (dashed lines) and the sum of all emission components (thick line). We cut the galaxy emission from our individual emission components in both panels and additionally show the total emission including galaxies (grey solid line).

The projected energy dependent photon index is a well defined continuous quantity as it is defined through


where is the surface brightness of -rays with energies at the projected radius . Using maps, such as the ones in Fig. 2, we can extract the photon index between the energies to . Close to the pion bump (see Fig. 1) at (energy of a photon originating from a decaying at the threshold of the hadronic p-p reaction) the photon spectrum has a convex shape. This is characterized by a flatter photon index compared to the asymptotic limit, . Calculating at two discrete energies results in a slightly steeper value for than its continuous counterpart at .

Since we are interested in comparing the morphology of clusters to spectra, we calculate the projected photon index for g72a (see Fig. 3). We compare for an energy regime accessible to the Fermi space telescope of (left panel) to accessible for IACTs (right panel). In the Fermi regime, we find a median value of in the central regions of the cluster and a value in the periphery. The reason is that the total emission in the central regions of the cluster is dominated by the pion decay emission at 100 MeV with a lower spectral index than pIC due to the pion bump. In the periphery and the WHIM, where the pIC contributes substantially to the total emission, intermediate shocks with are typical (Pfrommer et al., 2006). Using the spectral index of the electron equilibrium spectrum, , where for and , results in the observed steepening of in the periphery.

We now turn to the energy region important for IACTs with the photon index . In the central regions the photon index traces the proton spectral index since this spatial region is dominated by the asymptotic regime of the pion emission. Moving towards the periphery, the photon index steepens to , despite the presence of strong external shocks as well as accretion shocks that efficiently accelerate electrons on these large scales. The reason for this steepening is the exponential cutoff in the pIC emission. Increasing the energy or the distance from the cluster results in an even steeper photon index.

3.2 Emission profiles

The profiles of different non-thermal -ray emission processes without galaxies are shown in Fig. 4 for, a large CC cluster (g8a, left) and large post-merging cluster (g72a, right). The secondary IC emission traces the dominating pion decay emission because to zeroth order, both components depend on , where is the gas number density and the CR number density. This would be exactly true if the magnetic field was smaller than G. In this case, the CRe population would exclusively cool by means of IC emission. Since we assume the central magnetic field to be larger than , a fraction of the CRe energy is radiated through synchrotron emission into the radio, causing a larger discrepancy of the sIC emission compared to the pion emission in the center.

In contrast to the centrally peaked secondary emission components, the average primary IC emission shows a rather flat surface brightness profile which can be nicely seen in our cool core cluster g8a. This is because we see the strong accretion shocks that efficiently accelerate CRes (in terms of the injected energy density) in projection. There are noticeable exceptions in the centers of both clusters: accreting small sub-clumps dissipate their gravitational energy through weak shocks in the larger core regions of clusters. However, once these weak shock waves encounter the (over-cooled) centers of our simulated clusters they transform into strong (high Mach number) shock waves. These inject a hard population of primary CR electrons which causes the centrally peaked and bright pIC emission. We also observe an excess pIC emission at a radius of  Mpc in g72a. This can be traced back to a prominent merger shock wave with a Mach number up to that accelerates primary CRes (see also Fig. 2).

The total emission flattens out in the cluster exterior close to because of two reasons: (1) there the pIC emission contributes significantly to the total emission because of strong merger shocks as well as accretion shocks, and (2) subhalos in the periphery that have not yet merged with the larger halo contribute to the pion decay induced emission. In this regime, the halo-halo correlation term starts to dominate the average density profile of a cluster with its characteristic flattening (Hayashi & White, 2008). This naturally translates to the pion emission profile that tightly correlates with the gas density distribution.

Figure 5: The -ray number flux weighted by photon energy of our Coma-like cluster g72a: core region (left panel) and WHIM region (right panel). The pion decay flux is shown in blue color, the secondary IC in red, and the primary IC emission in green. The solid IC lines assume a central magnetic field strength of while the dashed lines assume . The lower panels show the photon spectral index for each emission component, where the colors are the same as in the upper panel. Note that in the core region () the pion decay induced flux dominates the emission. In the WHIM region () the pion decay and sIC emission components are much smaller compared to the central part. Contrary, the pIC component is boosted in the WHIM due to accretion shocks onto our simulated cluster. Nevertheless, its emission falls slightly short of the pion decay.

The ratio between pion decay and sIC emission can be estimated analytically by calorimetric considerations. To this end we compare the CR energy spectrum, , at the respective CR energies which give rise to -ray emission at some specific photon energy, say  GeV. Additionally we have to include a factor, , that accounts for the possibility that the CRes do not only emit IC -rays but also radio synchrotron radiation. Using the hadronic branching ratios for the production of pions, and , as well as the multiplicities for the decay products in the respective decay channels, and , we obtain


Here is the CR proton spectral index between the CR energy that give rise to pion decay flux at 1 GeV and the CR energy that gives rise to sIC flux at 1 GeV. The -ray source function for pion decay and sIC ((Pfrommer & Enßlin, 2004)) is denoted by and , respectively. Finally, the factor , for magnetic fields much smaller than the CMB equivalent magnetic field, otherwise . In the region close to , the magnetic field is about 2.4  in our model (Table 2), which implies an emission ratio of about 50. For smaller photon energies than 1 GeV, the pion decay -ray emission falls below the asymptotic power-law due to the characteristic pion bump. This effect implies a lower ratio of about 40 (instead of the expected ratio of about 100) for the emission above 100 MeV in Fig. 4.

We also show the total surface brightness profile with galaxies in Fig. 4 (grey line). The resulting profile shows a boosted emission by about a factor two compared to the one where we exclude galaxies from the surface brightness. The entire population of these galaxies takes up only a negligible volume so that the volume weighted CR pressure is almost the same in either case, when taking these galaxies into account or not. We note that this bias needs to be addressed when deriving average CR pressure contributions from the cluster’s -ray emission (Aleksić et al., 2010). Especially in the inner parts, the profile is very inhomogeneous. Since galaxies follow an approximate Poisson distribution and since the inner radial bins of the profile sample only few galaxies, we naturally obtain a larger Poissonian scatter across the inner radial bins.

3.3 Emission spectra from the cluster core and WHIM

The central parts of clusters are characterized by high gas and CR densities, and magnetic fields – at least compared to average values of the ICM. Even though the cluster core region only makes up a fraction of the total volume of the ICM, the high densities result in a significant -ray flux contribution to the total flux from the cluster. In contrast to the cluster center, the WHIM is characterized by on average low gas and CR densities, and magnetic fields. The low densities cause a smaller total -ray flux from this region compared to the cluster core regions. However, the WHIM of the super-cluster region contains a large number of galaxies and groups that are accreted onto the cluster. This generates more shocks compared to the cluster core region. The cluster characteristics in the two regions give rise to different normalizations of the individual -ray emission components, but with a similar shape. The shape of the emission components from the different regions agrees with that of the entire cluster as shown in Fig. 1. The emission can be summarized as follows. The -decay is characterized by the so-called pion bump followed by a concave curvature and a diffusive break. The pIC emission component shows a power-law followed by an exponential cutoff while the sIC component has a power-law with similar index that is however followed by the Klein-Nishina break.

In Fig. 5 we show the -ray number flux weighted by photon energy from different regions of our Coma-like cluster g72a that we place at the distance of 100 Mpc. The left figure shows the -ray number flux within the core region (), where the -decay dominates over the sIC component that itself is sub-dominant to the pIC component. The surface brightness profile of the -decay is sufficiently flat in the core region so that the -ray flux is dominated by the outer scale around where most of the volume is. Hence the -decay flux is largely insensitive to numerical inaccuracies of our modeling of the physics at the very center of the cluster. Both the pIC and sIC emission components have a larger -ray flux in our models with weak central magnetic fields (, ) compared to our models with strong central magnetic fields (, ). The sIC emission with is characterized by on scales around the core radius which is the region that dominate the flux. Relative to the pion decay emission, the sIC is suppressed by a factor of which can be understood by considering hadronic decay physics and the fact that the CR energy spectrum, , is decreasing as a function of proton energy (see equation 19). Even though the pIC emission component is sub-dominant, it shows a rather flat spectral index. This implies only a few strong shocks that are responsible for the electron acceleration. These merging shock waves are traversing the cooling core region in the cluster center. We caution the reader that the over-cooling of the cluster centers in our simulations possibly overestimates the true shock strengths and numbers which also results in an artificially enhanced pIC emission. At high energies, the electron cooling time is smaller than the time scale for diffusive shock acceleration which causes an exponential cutoff in the electron spectrum which is passed on to the pIC spectrum. The energy scale of the cutoff (combining equations 6 and 12) scales with the magnetic field which causes the low magnetic field model to turn down faster than the large magnetic field model. Note that the second cutoff in the figure for the small central magnetic fields is caused by a small fraction of the particles that have unusually high electron energy cutoff. The lower panel shows the photon spectral index defined in equation (18), where for the -decay emission after the pion bump that slowly flattens out with energy. For both the pIC and sIC emission the photon spectral index above the MeV regime up to at about 100 GeV. The reason for the flat spectra is that the pIC emission is generated by a few strong shocks in the center, while the sIC emission is caused by protons in the flat high energy part of CR spectrum.

Now we turn to the -ray spectra in the WHIM which are shown in the right panel of Fig. 5 for our g72a cluster. Here we define the WHIM by the emission in the region as seen in 2D projection of the cluster. We see that the pion decay and sIC are suppressed by a factor that is larger than 10 compared to the flux within the core region since the emission correlates with the CR- and gas densities. The suppression of -ray flux emitted by the intergalactic medium is expected to be much greater due to the large density decrease. The presence of small groups in the super-cluster region with densities that are much larger than the average density in the WHIM partially counteracts the flux suppression. Contrary to the pion decay and sIC component, the pIC emission is boosted by a factor of a few compared to the center because of the larger spatial region in the WHIM that contains a greater number of shocks. This leads to comparable flux from the pIC and pion decay emission above the energy of the pion bump in the super-cluster region. Note however, that this flux is still sub-dominant compared to the pion decay flux emitted by the cluster core region. The different central magnetic fields do not play any significant role in the power-law regime in the WHIM since , which implies that the CRes mainly cool through IC emission. However, note that the pIC cutoff is shifted towards lower energies for these smaller magnetic fields since (as derived by combining equations 6 and 62). In the lower panel we show that the photon spectral index is steeper in the WHIM for all three emission components compared to the core region. The photon index is about 1.3 for both pIC and sIC below 10 GeV and about 1.4 for the pion decay above the energy of the pion bump. The reason for the steeper pIC is because most of the energy that is injected into primary electrons comes from multiple intermediate-strength shocks (accretion shocks), while in the cluster center the pIC emission is build up from a few strong shocks in the over-cooled center (merger shocks). The steeper sIC and pion decay photon indices are caused by the slightly steeper CR spectrum present in the WHIM of the g72a cluster (cf. Fig. 7).

4 The CR proton distribution

In this section we investigate the CR proton spectrum that we obtain from our simulations and discuss the relevant physics that gives rise to it. We explore the variance of the spectrum across our cluster sample and within individual clusters and show that it obeys a universal spectral shape. In addition, we study the spatial profile of the CRs within a cluster as well as across our cluster sample and find it to be approximately universal. This universal behavior enables us to construct a semi-analytic CR spectrum and to compute the -ray spectrum as well as other secondary decay spectra (electrons, neutrinos) semi-analytically.

Figure 6: Spectral universality across our cluster sample: we show the median CR proton spectrum (solid line) as a function of dimensionless CR momentum for our sample of 14 galaxy clusters in Table 1 together with the 68 percentiles (dotted lines). The CR spectrum of every cluster, , has been obtained by volume weighting the individual spectra of our SPH particles and normalizing them at . Note the small variance between different clusters below the diffusive break, which indicates a universal CR spectrum for galaxy clusters. The large scatter in the momentum regime is a consequence of the strong radial dependence of the diffusive break. The lower panel shows the CR spectral index of the cluster sample.

4.1 A universal CR spectrum

In Fig. 6 we show the median CR spectrum of all 14 galaxy clusters as well as the associated spectral index alpha. The CR spectrum of every cluster, , has been obtained by volume weighting the individual spectra of our SPH particles which have been normalized at the dimensionless proton momentum . Before discussing the spectral shape we note that it shows a remarkably small variance across our cluster sample which indicates a universal CR spectrum for galaxy clusters. There are three important features in the spectra.

  1. The spectrum shows a low-momentum cutoff due to efficient Coulomb cooling at these low momenta with : the CR energy is transferred into the thermal energy reservoir through individual electron scatterings in the Coulomb field of the CR particle (Gould, 1972). The Coulomb time scale of a mono-energetic CR population is very short, , where we show a momentum scaling that is valid only for the relevant non-relativistic regime. The Coulomb time scale for a power-law population of CRs can be significantly longer, , where we assumed a low-momentum cutoff and (Enßlin et al., 2007). This, however, is still short compared to the dynamical time scale . Hence, we expect the formation of an equilibrium cutoff of our CR spectrum around for the typical number densities that we encounter at the cluster cores. Note that in the presence of cooling processes only, the equilibrium cutoff is determined from the competition between Coulomb cooling and hadronic losses and converges around q = 1.685 if the cooling time is sufficiently short (see Jubelgas et al. (2008) for a detailed discussion). The reason for that is that Coulomb cooling shifts the cutoff to higher momenta, as the CRs with low momenta are transferred to the thermal reservoir. At high momenta, the cooling time due to hadronic interactions is shorter than the Coulomb cooling time. Hadronic cooling effectively removes the CRs with high energy and moves the cutoff towards lower momenta.

  2. In the momentum range between , the spectrum has a concave shape in double-logarithmic representation, i.e. it is a decreasing function with energy in the GeV/TeV energy regime. This is quantified by the momentum dependent spectral index (shown in the lower panel of Fig. 6) which ranges from at energies above a GeV to around 100 TeV. This spectral shape is a consequence of the cosmological Mach number distribution that is mapped onto the CR spectrum (Pfrommer et al., 2006). This mapping depends on the shock acceleration efficiency as a function of shock strength as well as on the property of the transport of CRs into galaxy clusters. Nevertheless, we can easily understand the qualitative features: the typical shocks responsible for CR acceleration are stronger at higher redshift since they encounter the cold unshocked inter-galactic medium.10 This implies the build-up of a hard CR population. Since the forming objects have been smaller in a hierarchically growing Universe, their gravity sources smaller accretion velocities which results in smaller shock velocities, . Hence the energy flux through the shock surfaces, , that will be dissipated is much smaller than for shocks today. This causes a lower normalization of this hard CR population. With increasing cosmic time, more energy is dissipated in weaker shocks which results in a softer injection spectrum. Despite the lower acceleration efficiency, the normalization of the injected CR population is larger that that of the older flat CR population which yields an overall concave spectral curvature. We will study the details of the CR acceleration and transport that leads to that particular spectrum in a forthcoming paper (Pinzke & Pfrommer, in prep.).

  3. There is a diffusive break in the spectrum at high momenta where the spectral index steepens by 0.3. The CRs above these energies can diffusively escape from the cluster within a Hubble time. The particular value of the steepening assumes that the CRs scatter off Kolmogorov turbulence. Using twice the virial radius of each cluster, we find that the diffusive break varies between the momenta , dependent on the characteristic size of a cluster (equation 3).

Figure 7: Spectral universality within a cluster – variance and radial dependence: we show the CR proton spectrum of a large merging cluster, g72a, as a function of dimensionless CR momentum . Top: median of the volume weighted and normalized CR spectrum (solid) and particle-by-particle variance as indicated by the 68 percentiles (dotted line). Bottom: volume weighted CR proton spectrum (median and 68 percentiles) for different radial bins. Those are represented by different colors: core region of the intra-cluster medium in red, the intermediate cluster scales in orange, the periphery of the ICM in green, and the WHIM in blue. The spectrum for the different radial bins has been offset vertically by a constant factor of four, with decreasing amplitude for increasing radius. The shift of the cutoff to a higher momentum for smaller radius is due to the smaller CR cooling time scales in denser regions. Note that the cosmic ray spectrum inside the virial radius is almost independent of radius. This behavior is not only observed for this cluster, but persistent across our cluster sample.

We now turn to the question on how universal is our CR spectrum within a cluster. In Fig. 7 we find that the intrinsic scatter of our CR spectra within a cluster is larger than the scatter among clusters. The reason being that formation shocks are intermittent as mass is accreted in clumps an not continuously. Hence there is a high intrinsic variance of the CR spectrum for similar fluid elements that end up in the same galaxy cluster. However, averaging over all fluid elements that accrete onto a galaxy cluster results in a very similar spectrum since different galaxy cluster experience on average the same formation history.11 We note that the spectral variance of g72a is representative for all the CR spectra in our sample.

We study the radial dependence of the CR spectrum for g72a at the bottom of Fig. 7. Inside the cluster, the spectral shape does not strongly depend on the radius. This is a crucial finding as it enables us to separate the spectral and the spatial part of the CR distribution. The level of particle-by-particle variance is similar to that of the total cluster spectrum. Also noticeable is the increasing low-momentum cutoff as we approach the denser cluster center. This is due to enhanced Coulomb losses and to a lesser extent increased adiabatic compression that the CR distribution experienced when it was transported there. Outside the cluster g72a, for radii , we observe a considerably steeper CR spectrum at CR energies of  TeV compared to the cluster region. We note that this behavior is not universal and we observe a large scatter of the CR spectral indices among our cluster sample at these large radii. This behavior might be caused by the recent dynamical activity of the cluster under consideration but a detailed characterization goes beyond the scope of this work and will be postponed (Pinzke & Pfrommer, in prep.).

4.2 The spatial distribution of CRs

We have shown that the CR spectrum is separable in a spectral and spatial part. To this end, we quantify the spatial part of the CR spectrum by taking the volume weighted CR spectrum in each radial bin (see Section A.1 in the appendix for derivation),


Here we assume , the subscripts M and V denote mass- and volume weighted quantities, respectively, and we introduced the dimensionless normalization of the CR spectrum,


where is the CR number density. We now have to take into account physical effects that shape the profile of . Those include acceleration of CRs at structure formation shocks, the subsequent adiabatic transport of CRs during the formation of the halos as well as non-adiabatic CR cooling processes; primarily hadronic interactions. At the same time we have to consider the assembly of the thermal plasma and CRs within the framework of structure formation that is dominated by CDM. Hierarchical structure formation predicts a difference in the halo formation time depending on the halo mass, i.e. smaller halos form earlier when the average mass density of the Universe was higher. This leads to more concentrated density profiles for smaller halo masses; an effect that breaks the scale invariance of the dark matter halo profile (Zhao et al., 2009). The central core part is assembled in a regime of fast accretion (Lu et al., 2006). This violent formation epoch should have caused the CRs to be adiabatically compressed. In the further evolution, some cluster have been able to develop a cool core which additionally could have caused the CRs to be adiabatically contracted. On larger scales, the gas distribution follows that of dark matter (at least in the absence of violent merging events that could separate both components for time scales that are of order of the dynamical time). On these scales, the CR number density roughly traces the gas distribution (Pfrommer et al., 2008). Overall, we expect the spatial CR density profile relative to that of the gas density to scale with the cluster mass. If non-adiabatic CR transport processes have a sufficiently weak impact, these considerations would predict flatter -profiles in larger halos as these halos are less concentrated.

Our goal is to characterize the trend of -profiles with cluster mass and to test whether we can dissect a universal spatial CR profile. To this end, we adopt a phenomenological profile that allows for enough freedom to capture these features as accurately as possible. Hence, we parametrize with shape parameters that include a flat central region given by , a transition region where the location is denoted by and the steepness by , and finally a flat outer region denoted by , through the equation


The core regions in our radiative simulations show too much cooling, and we possibly lack of important CR physics such as anisotropic CR diffusion that could smooth out any strong inhomogeneity at the center. Hence we adopt a conservative limit of the central profile value of .

Figure 8: The top panel shows the profile of the dimensionless normalization of the CR spectrum, . We show the mean and the standard deviation across our cluster sample which has been subdivided into two different mass intervals: large- (red), and low-mass clusters (blue) representing the mass range , and . The solid lines show the best fit to equation (22). The lower three panels show the mass dependence of the quantities which parametrize for low mass clusters (blue circles) and large mass clusters (red circles). The top small panel shows the asymptotic for large radii (), where each cross shows at for each cluster. The middle panel shows the transition radius , and the bottom panel shows the inverse transition width denoted by .

In the top panel of Fig. 8 we show the mean of the profiles from the cluster simulation sample. We subdivide or sample in two different mass intervals: large mass clusters (top red), and small mass clusters (bottom blue), in the mass range , and , respectively.12 The error-bars show the standard deviation from sample mean and the solid lines the best fit that will be discussed in more detail in Section 4.3. The profile of our large mass clusters is almost flat and shows only a very weak radial dependence. In contrast, the profile of our small clusters has a rather steep and long transition region and is increasing towards the center. The difference in normalization, transition width, and transition radius between low mass and large mass clusters indeed suggests that scales with the cluster mass in a way that we anticipated. We quantify the mass scaling of these shape parameters by a power-law fit to the small and large clusters in the three lower panels of Fig. 8: the value of in the asymptotically flat regime in the cluster periphery, (top); the transition radius, (middle); and the steepness of the transition, (bottom). As expected, we find that all three quantities increase slowly with mass. We additionally show the values of at for each cluster (black crosses). Within the scatter, these values are consistent with the mass scaling found by the profile fits in our two cluster mass bins. We have shown that there exists an almost universal spatial CR profile after taking into account the weak trends of the profile with cluster mass. Note that the particle-by-particle variance of the spatial CR profile within a cluster (that we address in Section A.3 in the appendix) additionally supports the conclusion of a universal spatial CR profile.

4.3 A semi-analytic model for the spatial and spectral CR distribution

In our simulations we use a CR spectral description which follows five different CR proton populations; each being represented by a single power-law. The CR populations are chosen such that we accurately can capture the total CR spectrum in the entire momentum space (a convergence study on the number of CR populations is presented in the appendix A.2). We want to model this spectrum analytically with as few CR populations as possible, but at the same time preserve the functional form of a power-law of each population. In that way we can easily obtain the total CR spectrum by superposition and apply a simple analytic formula to estimate the induced radiative processes.

In detail, we use a total CR proton spectrum that is obtained by summing over the individual CR populations (equation 2); each being a power-law in particle momentum with the total CR amplitude ,


where we assume universal spectral shape throughout the cluster. The normalization of depends on maximum shock acceleration efficiency , where and (functional will be studied in Pfrommer in prep.). Denoting the amplitude of each CR population by , where the number of spectral bins might be different from the five chosen in the simulations, we construct the relative normalization for each CR population


Here we have assumed that is only a weak function of radius. This is explicitly shown in Fig. 7, where the functional form is almost independent of the radius. The two energy breaks in the CR spectrum are represented by


The first term in equation (25) ensures the low-momentum slope as appropriate from the phase space volume that is populated by CRs (Enßlin et al., 2007) and the last term accounts for diffusive losses of CRs that steepen the spectrum by assuming Kolmogorov turbulence. The shape parameter determines the size of the transition in momentum space. Our choice of ensures a fast transition from one regime to the other. The low-momentum break is determined to from Fig. 9, and the high-momentum break is derived from equation (3).

Figure 9: Cosmic ray proton spectrum as a function of dimensionless CR momentum for our sample of 14 galaxy clusters (Table 1). The main panel shows median of the cosmic ray spectrum across our cluster sample. The spectrum of each cluster has been normalized at (black crosses). The blue solid line shows the best fit triple power-law to the simulation data. The light blue area shows the 68 percentile of the data points across our cluster sample. The bottom panel shows the difference between the relative fit and the simulation data (blue solid line) which amounts to less than 2 percent.

To capture the spectral universality of our cluster sample, we fit in Fig. 9 the median of our cluster sample of the CR spectrum with a triple power-law function. Because we normalize the spectrum at , this ensures that the normalized spectrum becomes independent of radius (cf. equation 23). The triple power-law fit represented by the blue line in the upper panel of Fig. 9 resulted in


The data from the simulation is denoted by black crosses, together with the 68 percentiles spread shown by the light blue area. In the lower panel, we show the relative difference between the simulation and the fit which indicates a difference of less than two percent from the GeV to PeV energy regime.

The spatial part of the CR spectrum is derived from the fit to the mean in the top panel of Fig. 8. The solid lines show the best fit to using equation (22). We find that the central value is the most conservative value that still provides a good fit for both mass intervals. Note that there is a large uncertainty in this value due to insufficient data in the center13, implying that it should be treated as a lower limit. However, the gamma-ray flux depends only weakly on the exact value of since most of the flux reside from the region outside the transition region. The mass dependence of is quantified in the three lower panels of Fig. 8, where we fit a power-law in mass for the normalization (top), the transition radius (middle), and the steepness of the transition region (bottom). We obtain the following relations,


5 Semi-analytic model for the -ray emission

In this section we derive a semi-analytic formula for the integrated -ray source function that is based on our semi-analytic CR distribution. Using the gas density profile of the cluster along with its virial radius and mass, this formula enables us to predict the dominating -decay induced -ray emission from that cluster. The density profiles can either be inferred from simulations or from X-ray data. To test the semi-analytic source function, we use density profiles from our cluster simulations. We also make -ray flux predictions for the Coma and Perseus cluster using their true density profiles as inferred from X-ray observations.

5.1 Schematic overview and semi-analytic formulae

Figure 10: Schematic overview of how our semi-analytic framework relates to the simulated CR and -ray quantities. From the cluster simulations we derive the CR spectrum for each SPH particle, . Going clockwise, the integrated -ray production rate for each SPH particle is derived within the framework of Pfrommer et al. (2008). The radial profile of the integrated -ray production rate, , is obtained by volume weighting this quantity of individual SPH particles. Returning to and going counter-clockwise instead, we can directly perform a volume-weighting of the CR spectrum in each radial bin, yielding . Fitting both the spatial and spectral part provides us with a semi-analytic model for the spectrum, . The semi-analytic model for the integrated -ray production rate, , is derived within the framework of Pfrommer et al. (2008). The explicit form of can be found in equations (30)-(34), while and are implicitly given by equation (32). Comparing the flux from the simulations and our semi-analytic model shows good agreement for our cluster sample.

In Fig. 10 we show a schematic overview of the simulated CR spectrum and the integrated source function together with the mapping to our semi-analytic model. From the cluster simulations we derive the CR spectrum describing the phase space density of CRs. When taking the spherical (volume weighted) average of , we obtain . We fit the spectral and spatial part of this function separately. The semi-analytic model for the CR distribution in clusters that we provide in equation (23) can now be used to predict the secondary radio synchrotron emission and the hadronically induced neutrino and -ray emission from galaxy clusters. Following the formalism in Pfrommer et al. (2008), we obtain the volume weighted and energy integrated and omnidirectional (i.e integrated over the solid angle) -ray source function due to pion decay14,