text:
# Gamma-Rays from Intergalactic Shocks

## 1 Introduction

In the past three decades, clear evidence has emerged indicating the existence of an extragalactic $`\gamma `$-ray background (EGRB). The origin of this radiation, however, has remained highly speculative. The first unambiguous detection of isotropic extragalactic $`\gamma `$-ray emission was obtained by the SAS-2 satellite at 1977. Subsequent experiments, especially the EGRET instrument aboard the *Compton Gamma Ray Observatory* (CGRO), have confirmed the existence of this radiation. These measurements indicate a generally isotropic (fluctuation amplitude $`20\%`$ on $`20^{}`$ scales) spectrum in the energy range 30 MeV-120 GeV, well fitted by a single power law, with photon number density per energy interval $`dJ/dϵϵ^{\left(2.10\pm 0.03\right)}`$ Sreekumar et al. (1998). Known extragalactic sources, such as blazars, account for less than 7% of the EGRB, and unidentified blazars probably contribute no more than 25% of the radiation (Mukherjee & Chiang 2000; but for a different model see Stecker & Salamon 2001). Thus, most of the EGRB flux originates from a source yet to be identified.

Recently, Loeb & Waxman (2000) have proposed a model which attributes most of the EGRB flux to emission from large-scale structure-forming regions in the universe. Strong collisionless, non-relativistic shock waves characteristic of these regions accelerate particles to highly relativistic energies by the Fermi-acceleration mechanism, first discussed as sources of ultra-high energy cosmic rays by Norman, Melrose & Achterberg (1995) and by Kang, Ryu & Jones (1996). Electrons, accelerated in such shocks to highly relativistic energies with Lorentz factors up to $`\gamma _{\text{max}}10^7`$, scatter a small fraction of the cosmic microwave background (CMB) photons up to $`\gamma `$-ray energies, as inverse-Compton radiation. The estimated energy flux from this process per decade in photon energy, $`ϵ^2(dJ/dϵ)1.5\text{ keV}\text{ cm}^2\text{ s}^1\text{ sr}^1`$, is in good agreement with the detected EGRB, for typical cosmological parameters, provided that the shocked matter accumulated a present-day mass-average temperature $`T10^7\text{ K}`$, and that a fraction $`\xi _e0.05`$ of the associated thermal energy was transferred to the relativistic electrons. In this model, nearby rich galaxy clusters are expected to be bright $`\gamma `$-ray sources, as they host strong, dense shocks Loeb & Waxman (2000); Totani & Kitayama (2000).

Kawasaki & Totani (2001) report a correlation between a sample of high-latitude steady unidentified EGRET sources and pairs or groups of galaxy clusters, where strong accretion and merger shocks are expected. However, this sample is too small and the correlation too weak for it to be statistically significant. Enßlin et al. (2001) suggest that various features of the asymmetric radio plasma ejected from the giant radio galaxy NGC 315 can provide the first indication of a structure-formation shock in a filament of galaxies. According to the above model, such a shock should be accompanied by $`\gamma `$-ray emission, yet to be detected.

The model presented above has several interesting implications, other than an identification of the origin of the EGRB: $`\gamma `$-ray emission could be the first tracer of the undetected warm-hot intergalactic medium (WHIM), gas at temperatures $`10^5\text{ K}T10^7\text{ K}`$ which is believed to contain a large fraction of the baryons in the low redshift universe (e.g. Davé et al. 2001a). The angular fluctuations in the EGRB should, according to the model, reflect the underlying forming structure, and thus serve as a diagnostic of intergalactic shocks. The relativistic electrons emit synchrotron radiation, as they gyrate in intergalactic magnetic fields. The resulting radio radiation, although weaker than the CMB, could have a fluctuating component dominating the CMB fluctuations at low ($`\nu 10\text{ GHz}`$) frequencies Waxman & Loeb (2000). Inverse-Compton emission from shock-accelerated electrons in the hard X-ray (HXR) and extreme ultra-violet (EUV) bands was recently studied by Miniati et al. (2001). They find that such emission could account for a small fraction of the measured EUV excess, and for the entire HXR excess in a few reported galaxy clusters, provided that the latter contain relatively strong merger or accretion shocks.

In order to test this model, calibrate its free parameters and extract useful predictions, one must follow the evolution of structure in the intergalactic medium (IGM) in detail. The complex, non-linear processes and many scales involved, make this a difficult problem, rendering numerical simulations a powerful and indispensable tool. In this paper, we apply cosmological simulations to the problem of non-thermal emission from shock–accelerated electrons, and demonstrate how various EGRB predictions can be efficiently extracted from them. For this purpose, we have developed the necessary tools to localize intergalactic shocks in the simulations, inject the corresponding shock-accelerated electrons into the gas, and calculate the resulting radiation. We consider a hydrodynamic simulation of a $`\mathrm{\Lambda }\text{CDM}`$ universe, currently considered the most successful theory of structure formation, which describes a flat universe containing dark matter and presently dominated by vacuum energy. The relativistic electrons are assumed to carry a fraction $`\xi _e`$ of the shock thermal energy. We adopt the value $`\xi _e=0.05`$, which we show is inferred from collisionless non-relativistic shocks in the interstellar medium.

Overall, we find an emitted inverse-Compton spectrum well-described by a mildly broken power law, $`dJ/dϵϵ^{2.04}`$ in the energy range $`10\text{ keV}1\text{ GeV}`$ and $`dJ/dϵϵ^{2.13}`$ in the energy range $`10\text{ GeV}5\text{ TeV}`$, in good agreement with the prediction of Loeb & Waxman (2000). The calculated energy flux (per decade in photon energy), $`ϵ^2(dJ/dϵ)=50160\text{ eV}\text{ cm}^2\text{ s}^1\text{ sr}^1`$, constitutes $`10\%`$ of the EGRB flux and is smaller by a factor of $`10`$ than estimated above, mainly due to the lower present-day temperature reached by the simulation, $`T_04\times 10^6\text{ K}`$, and its slower heating rate. We predict the existence of $`\gamma `$-ray halos, associated with large-scale structure. Although too faint to account for the EGRET unidentified $`\gamma `$-ray sources, several of these halos should be detectable by future experiments, such as the *Gamma-Ray Large-Area Space Telescope* (GLAST, planned 111See http://glast.gsfc.nasa.gov to be launched in 2006), provided that $`\xi _e0.03`$. The number of detectable sources is sensitive to the fraction of shock-energy transferred to the relativistic electrons, with $`15`$ well-resolved GLAST sources predicted for $`\xi _e=0.05`$. Detection of intergalactic shock-induced $`\gamma `$-ray sources may be used to calibrate the free parameter $`\xi _e`$. We examine images of such an object, as predicted to be produced by GLAST and by the *Major Atmospheric Gamma-ray Imaging Čherenkov* (MAGIC) telescope (expected 222See http://hegra1.mppmu.mpg.de/MAGICWeb to become operational during 2002). These images reveal an elliptic ring-like feature, of diameter corresponding to $`510\text{ Mpc}`$, surrounding a nearby ($`z0.012`$) rich galaxy cluster. This ring, tracing the cluster accretion shock, has luminous peaks along its circumference, probably indicating the locations of intersections with galaxy filaments, channeling gas into the cluster.

We describe the Loeb-Waxman model in §2. We first consider the gross characteristics of intergalactic shocks based on order of magnitude estimates. In this context, we parameterize the distribution of accelerated electrons, and discuss the inverse-Compton and synchrotron radiation they emit. The simulation we chose to study and its basic predictions are described in §3. We discuss the underlying theory and cosmological model, present some characteristics of the simulated universe and examine its structure. Section 4 describes how EGRB predictions can be extracted from cosmological simulations in general and then focuses on the method applied to the simulation discussed in §3. In §5, we present our results regarding the radiation predicted from the simulation. We describe the resulting inverse-Compton spectrum, analyze the corresponding $`\gamma `$-ray sky maps and study the features and distribution of the simulated $`\gamma `$-ray sources. In §6 we discuss the implications of our results for future experiments and models. Appendix A contains some formulae used to calculate and integrate the radiation emitted by the accelerated electrons. Appendix B presents an algorithm devised to identify point sources in our simulated maps of the $`\gamma `$-ray sky.

## 2 Theory

This section describes the model for non-thermal emission from structure formation shocks proposed by Loeb and Waxman (2000), using simple order-of-magnitude estimates. First, we estimate the parameters of the large-scale shock waves involved in structure formation. Next, we discuss the high-energy electrons accelerated by such shocks. We evaluate the power index and cutoff of their distribution, and evaluate its normalization by estimating the efficiency of electron acceleration in intergalactic shocks. We conclude by calculating the radiation emitted by these electrons, mainly through inverse-Compton scattering off CMB photons. The estimates presented in this section needed for our numerical simulation are discussed in more detail in §4.

### 2.1 Intergalactic Shock Waves

Large-scale structure in the universe is believed to have evolved by the gravitational collapse of initially over-dense regions. As such an initial seed accreted increasing amounts of mass, large converging flows of matter were accelerated towards it, inevitably brought to a violent halt. This process resulted in large-scale shock waves, forming as massive gas flows collided with opposite flows or with the newly formed object. In the following, we concentrate our attention on the ’recent’ stages of structure formation, taking place at moderate to low redshifts ($`z3`$), and deduce some order of magnitude estimates of the resulting shock wave properties, following Cen & Ostriker (1999).

Linear theory of structure formation predicts that the length and mass scales of an object, forming (becoming non-linear) at these low redshifts, are roughly given by:

$$\lambda _{NL}(z)=\lambda _0\text{ }\frac{f(z)}{1+z}\text{and}M_{NL}(z)=M_0\text{ }f(z)^3\text{ ,}$$
(1)
where $`\lambda _0`$ and $`M_0`$ are the length and mass scales of structures forming at the present time and $`f(z)`$ is a slowly varying function of the redshift $`z`$, satisfying $`f(z=0)=1`$, such that smaller scales become non-linear first. These parameters are sensitive to the cosmological model. We use a ’concordance’ $`\mathrm{\Lambda }\text{CDM}`$ model of Ostriker & Steinhardt (1995) - a flat universe with normalized vacuum energy density $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, matter energy density $`\mathrm{\Omega }_M=0.3`$, baryon energy density $`\mathrm{\Omega }_B=0.04`$, Hubble parameter $`h=0.67`$, and initial perturbation spectrum of slope $`n=1`$ and normalization $`\sigma _8=0.9`$ \- to find:

$$\lambda _010h_{67}^1\text{ Mpc}\text{and}M_01.5\times 10^{14}h_{67}^1\left(\frac{\mathrm{\Omega }_M}{0.3}\right)M_{}\text{ ,}$$
(2)
where $`h_{67}h/0.67`$. The order of magnitude estimates obtained from linear theory in this subsection should be regarded with caution, as non-linear effects become important in the forming objects, and the phases of the perturbations can no longer be ignored. The distribution of the forming objects according to linear theory may be better estimated using the approach of Press & Schechter (1974).

Geometric considerations dictate that the gas in an object of size $`\lambda `$, collapsing over a time period $`t_{\text{coll}}`$, will achieve velocities of order $`v\lambda /t_{\text{coll}}`$. Since a structure that became non-linear at redshift $`z`$ will collapse over a time period of order $`t_{\text{coll}}H(z)^1`$, where $`H(z)`$ is the Hubble parameter, the formation of such an object is expected to involve shocks with velocity:

$$v_{\text{sh}}(z)\text{ }k\text{ }\lambda _{NL}(z)\text{ }H(z)=v_0\text{ }\frac{f(z)h(z)}{1+z}\text{ ,}$$
(3)
where $`k`$ is a dimensionless number of order unity, $`v_0k\lambda _0H_0`$ is the typical shock velocity at $`z=0`$ and $`h(z)H(z)/H(z=0)`$. For the $`\mathrm{\Lambda }\text{CDM}`$ model described above we find:

$$v_0700\text{ km s}^1\text{and}h(z)=\sqrt{\mathrm{\Omega }_\mathrm{\Lambda }+\mathrm{\Omega }_M(1+z)^3}\text{ .}$$
(4)
Hence, the shock waves resulting from converging flows during structure formation in the IGM have a non-relativistic velocity similar to those of shocks surrounding supernova remnants (where accelerated particles are observed) and they are steadily growing in both scale and velocity.

The baryonic component may be treated as an ideal gas with an adiabatic index $`\mathrm{\Gamma }=5/3`$, and is sufficiently rarefied such that the shocks concerned are collisionless. Using the Rankine-Hugoniot shock adiabat Landau & Lifshitz (1959), we can relate the shocked gas temperature to the shock velocity,

$$T(z)T_0\left[\frac{f(z)h(z)}{1+z}\right]^2,$$
(5)
where $`T_0[m_p/(\mathrm{\Gamma }k_B)](\alpha v_0)^2`$ is the temperature of gas, shocked at $`z=0`$, $`k_B`$ is the Boltzmann constant, and $`\alpha c_{s,d}/v_{\text{sh}}`$ is the ratio between the post-shock (downstream) sound velocity and the shock velocity (relative to the unshocked gas), approaching $`0.56`$ for a strong shock. For the $`\mathrm{\Lambda }\text{CDM}`$ model above we find:

$$T_010^7\text{ K}\text{ .}$$
(6)
As we show in §3.2, this estimate of $`T_0`$ is indeed valid only to order of magnitude and appears to exceed the true value by a factor $`3`$.

### 2.2 Shock-Accelerated Electron Population

#### 2.2.1 Acceleration Mechanism

Collisionless, non-relativistic shock waves are generally known to accelerate a power-law energy distribution of high-energy particles. This phenomenon has been observed in astrophysical shock waves on various scales, such as in shocks forming when the supersonic solar wind collides with planetary magnetospheres, in shocks surrounding supernovae (SNe) remnants in the interstellar medium, and probably also in shocks in many of the most active extragalactic sources: quasars and radio galaxies Drury (1983); Blandford & Eichler (1987). These power law distributions extend up to energies $`100\text{ keV}`$ in the Earth’s bow shock, to $`100\text{ TeV}`$ in SNe remnants Tanimori et al. (1998) where shock velocities are similar to those of intergalactic shocks, and over $`10^6\text{ TeV}`$ in Galactic cosmic rays.

Magnetic fields of amplitude $`B0.11\mu \text{G}`$ are measured in halos of galaxy clusters Kronberg (1994) and in the Coma super-cluster Fusco-Femiano et al. (1999), suggesting that the IGM in which the shocks propagate is significantly magnetized. The measured magnetic field energy constitutes a fraction $`\xi _B0.01`$ of the typical post-shock thermal energy in cluster environments Waxman & Loeb (2000).

Electrons are expected to be accelerated in intergalactic shocks by the Fermi acceleration mechanism. The velocity difference between the two sides of a shock front acts as a first-order acceleration mechanism ($`dE/dtE`$), provided that slow-moving elastic scattering centers are embedded in the gas. Alfvén waves or magnetic inhomogeneities can function as such centers, scattering electrons through small angles and causing them to repeatedly diffuse through a shock front. Diffusion of the energetic electrons downstream, away from the shock front, provides a Poissonian escape route. These effects lead to a characteristic power-law distribution of the high-energy electrons steady-state number density:

$$\frac{dn_e}{dE_e}(E)AE^s\text{ ,}$$
(7)
where $`A`$ is a normalization constant. The vast size and perpetual duration of these shocks enable them to dissipate significant fractions of their energies in this process.

The resulting distribution of high energy electrons can be calculated using the test-particle approximation Drury (1983); Blandford & Eichler (1987), where the test particles, accelerated by the shock, are assumed to have no effect on its structure. This approximation, valid for the acceleration of a small population of high energy electrons, predicts that the power-law index of the resulting population depends only on the kinematic structure of the shock front, and is given by:

$$s=\frac{r+2}{r1}\text{ ,}$$
(8)
where $`r`$ is the compression factor of the shock. For a shock in an ideal gas, $`r`$ is given by:

$$\frac{1}{r}=\frac{\mathrm{\Gamma }1}{\mathrm{\Gamma }+1}+\frac{2}{(\mathrm{\Gamma }+1)M^2}\text{ ,}$$
(9)
where $`Mv_{\text{sh}}/c_{s,u}`$ is the Mach number, and $`v_{\text{sh}}`$ and $`c_{s,u}`$ denote the shock velocity and the pre-shock (upstream) sound velocity, respectively. The shock waves of interest are generally strong, i.e. involve high ($`M1`$) Mach numbers. Hence, we find $`r(\mathrm{\Gamma }+1)/(\mathrm{\Gamma }1)=4`$, and thus $`s2`$.

In order to completely define the power-law distribution discussed, we must specify its normalization and the energy range over which it extends. The Fermi acceleration mechanism per-se does not impose stringent limits on the electron energies achieved, as long as the total energy carried by the accelerated electrons is small compared to the shock energy and the shock structure remains intact. The e-folding time for electron acceleration is given Loeb & Waxman (2000) by:

$$\tau _{\text{acc}}\frac{r_\text{L}c}{v_{\text{sh}}^2}2\times 10^4\gamma _7(B_7T_{d,7})^1\text{ yr}\text{ ,}$$
(10)
where $`r_\text{L}`$ is the electron Larmor radius, $`\gamma _7`$ is the electron Lorentz factor - $`\gamma _e`$ \- in units of $`10^7`$, $`B_7`$ is the magnetic field in units of $`0.1\mu \text{G}`$, and $`T_{d,7}`$ is the downstream temperature of the gas, in $`10^7\text{ K}`$. This time-scale is much shorter than the shock lifetime. The limit on electron energy is therefore not determined by the acceleration process, but rather by cooling processes: electrons will stop accelerating once their cooling rate overcomes their acceleration rate. The dominant cooling process for the relevant parameters is inverse-Compton scattering of the electrons off CMB photons, with a characteristic cooling time:

$$\tau _{\text{IC}}=\frac{3m_ec}{4\sigma _T\gamma _eu_{\text{CMB}}}2\times 10^5\gamma _7^1(1+z)^4\text{ yr}\text{ ,}$$
(11)
where $`u_{\text{CMB}}`$ is the CMB energy density. Synchrotron cooling of the electrons is less efficient, with a characteristic cooling time $`\tau _{\text{syn}}=(u_{\text{CMB}}/u_\text{B})\tau _{\text{IC}}10^3(B_7)^2(1+z)^4\tau _{\text{IC}}`$, where $`u_\text{B}`$ is the magnetic field energy density. Thus, inverse-Compton cooling is faster, as long as the downstream magnetic field satisfies $`B<3(1+z)^2\mu \text{G}`$. Equating $`\tau _{\text{acc}}`$ and $`\tau _{\text{IC}}`$ leads to the maximal Lorentz factor:

$$\gamma _{\text{max}}=3.3\times 10^7\left(B_7T_{d,7}\right)^{1/2}(1+z)^2\text{ .}$$
(12)
Since the electron energy loss is dominated by scattering of CMB photons, the flux of high energy photons produced by shock accelerated electrons depends only weakly on the magnetic field strength. The magnetic field determines, nevertheless, the highest energy of accelerated electrons, and hence the energy to which the inverse-Compton spectrum extends. For $`B0.1\text{ }\mu \text{G}`$, as observed in cluster halos and which corresponds to a fraction $`\xi _B0.01`$ of the post-shock thermal energy, the highest energy electrons up-scatter CMB photons to energy $`ϵ1\text{ TeV}`$. It should be pointed out here, that the observed sub-$`\mu \text{G}`$ fields may be representative of only the downstream (post shock) magnetic field strength, and that the upstream field may be significantly lower. In this case, the electron acceleration time will be longer, and the highest energy of accelerated electrons and of inverse-Compton photons may be significantly lower than estimated from equation (12). The strength of the upstream magnetic field depends on the processes leading to IGM magnetization, which are not yet known. Magnetic fields may be produced by turbulence induced by the large scale structure shocks themselves Kulsrud et al. (1997), by “contamination” by the first generation of stars Rees (1987) or radio sources Daly & Loeb (1990); Furlanetto & Loeb (2001), or by galactic winds (Kronberg, Lesch, & Hopp 1999). Under all these scenarios, the upstream magnetic field strength is not expected to be much smaller than the downstream field strength. Moreover, near equipartition fields are expected to be produced in collisionless shocks through the growth of electromagnetic instabilities, as indicated, for example, by observed supernovae remnants Helfand & Becker (1987); Cargill & Papadopoulos (1988) and $`\gamma `$-ray bursts Gruzinov & Waxman (1999); Medvedev & Loeb (1999), and there is evidence from supernovae remnant observations for an enhanced level of magnetic waves ahead of the shock front (Achterberg, Blandford, & Reynolds 1994) We therefore adopt for our estimates of the maximum electron energy magnetic field strengths corresponding to fixed $`\xi _B0.01`$, which reproduce the observed large scale field strengths.

We can now parameterize the normalization of the relativistic electron distribution, by assuming that a fraction $`\xi _e`$ of the thermal energy density induced by the shock, $`\mathrm{\Delta }u_{\text{th}}`$, is transferred to these electrons:

$$u_e=\xi _e\mathrm{\Delta }u_{\text{th}}=\frac{3}{2}\xi _ek_B\left(n_dT_dn_uT_u\right)\text{ ,}$$
(13)
where $`u_em_ec^2_1^{\gamma _{\text{max}}}\gamma (dn_e/d\gamma )𝑑\gamma `$ is the total energy density of relativistic electrons and $`u`$ and $`d`$ subscripts indicate upstream (pre-shock) and downstream (post-shock) values, correspondingly. The gas number density, $`n`$, is related to the baryon number density, $`n_{\text{bar}}`$, by $`nm=n_{\text{bar}}m_p`$, where $`m0.59m_p`$ is the average particle mass in an ionized gas with a hydrogen mass fraction $`\chi =0.76`$. For strong shocks $`s2`$ and $`n_uT_un_dT_d`$, thus $`u_eA\mathrm{ln}\gamma _{\text{max}}`$, and we find:

$$A\frac{3\xi _en_dk_BT_d}{2\mathrm{ln}\gamma _{\text{max}}}\text{ .}$$
(14)
#### 2.2.2 Acceleration Efficiency

The efficiency of particle acceleration in structure formation shocks, parameterized above by $`\xi _e`$, is an important parameter bearing linearly on our results, and thus deserves a special discussion. Although no existing model credibly calculates the acceleration efficiency, we can evaluate $`\xi _e`$ in intergalactic shocks by estimating the acceleration efficiency of other astrophysical shock waves.

The best analogy to intergalactic shocks may be found in supernova remnant (hereafter SNR) shocks, drawing on the similarity between the velocities of intergalactic shocks and SNR shocks, both of the order of $`10^3\text{ km s}^1`$. The physics of shock waves is essentially determined by three parameters: the shock velocity $`v_{\text{sh}}`$, the pre-shock density $`n_u`$ and the pre-shock magnetic field $`B_u`$. Although the plasma density is very different in intergalactic shocks and in SNR shocks, the density may be extracted from the problem by measuring time in units of $`\nu _{p,i}^1`$, where $`\nu _{p,i}`$ is the plasma frequency of the ions. The pre-shock magnetic field, parameterized by the cyclotron frequency of the ions, $`\nu _{c,i}`$, may not be scaled out of the problem as well. However, by comparing $`\nu _{c,i}`$ to the growth rate of electromagnetic instabilities in the shocked plasma, $`\nu _{\text{e.m.}}\nu _{p,i}v_{\text{sh}}/c`$, we find that for strong shocks $`(\nu _{c,i}/\nu \text{e.m.})^2=\frac{B_u^2/8\pi }{(1/2)m_pn_uv_{\text{sh}}^2}1`$. We thus assume there is a well behaved limit when this ratio approaches zero, suggesting that the pre-shock magnetic field has little effect on the characteristics of strong shocks. With this assumption, we expect to find much similarity between strong shocks in different environments, as long as their shock velocities are comparable.

The energy of the relativistic electrons accelerated in SNR shocks was calculated by several authors, using the multi-frequency emission from SNR’s, preferably remnants with TeV detection such as SN1006 (G$`327.61.4`$; Tanimori et al. 1998). Dyer et al. (2001) have modeled the multi-frequency emission from SN1006, finding a shock efficiency of $`\xi _e=5.3\%`$. The total energy they find in relativistic electrons constitutes a fraction $`\eta _e1.4\%`$ out of the *total supernova explosion energy* ($`W5\times 10^{50}\text{ erg}`$, as often quoted in the literature). Other authors Mastichiadis & de Jager (1996); Aharonian & Atoyan (1999) have estimated $`\eta _e1\%2\%`$ in SN1006, which corresponds to $`\xi _e`$ values in the range of $`3.8\%7.6\%`$. Ellison, Slane and Gaensler (2001) have modeled the emission from supernova remnant $`\text{G}347.30.5`$, estimating the energy of relativistic electrons to constitute at least $`1.2\%`$ \- and more likely $`2.5\%`$ \- of the shock kinetic energy flux, corresponding to $`2.5\%\xi _e5\%`$.

An independent, less reliable method for estimating the electron acceleration efficiency relies on the observation, that relativistic ions accelerated in SNR shocks are dynamically significant Reynolds & Ellison (1992); Baring et al. (1999), and thus carry a large fraction (tens of percents, e.g. $`50\%`$ according to Ellison et al. 2001) of the post shock thermal energy. The efficiency of electron acceleration may then be found by estimating the ratio $`\eta `$ between the relativistic electron energy and the energy of relativistic ions. One may easily estimate the value of $`\eta `$ in interstellar medium cosmic rays: using the electron-to-proton ratio in a given energy, e.g. $`1\%`$ at GeV energies Barwick et al. (1998), correcting for electron cooling and integrating over energies, we find $`\eta 3\%4\%`$, i.e. $`\xi _e1\%2\%`$. However, since there is no clear relation between $`\eta `$ in the interstellar medium and its value immediately behind the $`10^3\text{ km s}^1`$ shocks found in SNR’s and in the intergalactic medium, this estimate of $`\xi _e`$ is controversial.

To summarize, a shock efficiency of $`\xi _e0.05`$ is found from studies of supernovae remnants, characterized by shock velocities similar to those of intergalactic shocks and thus expected to exhibit similar behavior. Since there is some uncertainty in this parameter, values of $`\xi _e`$ in the range $`0.020.10`$ are considered plausible.

### 2.3 Resulting Radiation

In this subsection we calculate the radiation emitted by the relativistic electron population described above. As mentioned earlier, most of the energy is produced when the electrons inverse-Compton scatter CMB photons, resulting in a spectrum that extends up to TeV energies. A secondary radiation process is synchrotron emission from the electrons, gyrating in the magnetic field, resulting in a similarly sloped spectrum extending up to the infrared. Other radiative processes experienced by the electrons, such as Bremsstrahlung and Čherenkov radiation, are less important.

We begin by calculating the diffuse radiation, resulting from the various intergalactic shocks at recent epochs (moderate to low redshifts). We then focus our attention on a single forming structure and the angular distribution of relevant sources.

#### 2.3.1 Diffuse Radiation

Inverse-Compton scattering describes the interaction of a free electron with radiation, when the electron kinetic energy exceeds the photon energy. In the ultra-relativistic limit, $`\gamma _e1`$, the average relative increase in photon energy is given by $`\mathrm{\Delta }ϵ_\gamma /ϵ_\gamma (4/3)\gamma _e^2`$. Hence, an electron with energy $`E_e=\gamma _em_ec^2m_ec^2`$ will scatter a CMB photon to an average energy:

$$ϵ_\gamma \frac{4}{3}\gamma _e^2ϵ_{\text{CMB}}=\zeta E_e^2\text{ ,}$$
(15)
where we have defined $`\zeta (4/3)(m_ec^2)^2ϵ_{\text{CMB}}`$, and $`ϵ_{\text{CMB}}3.83k_BT_{\text{CMB}}`$ is the average energy of a CMB photon. The extent of the spectrum is determined by the electron range of energies: The maximal electron Lorentz factor, $`\gamma _{\text{max}}3.3\times 10^7`$, is given by equation (12), and equation (11) imposes a minimal energy threshold, as electrons with Lorentz factors smaller than $`\gamma _{\text{min}}200`$ hardly radiate during a Hubble time $`10^{10}\text{ yr}`$. Plugging this range of Lorentz factors into equation (15) shows that the emitted photon spectrum extends from $`50\text{ eV}`$ up to $`\text{TeV}`$ energies.

The emitted spectrum may be found by assuming that each relativistic electron scatters a photon only once. This approximation, discussed in §4.3, gives qualitatively correct results and is exact for strong shocks. With this simplification, we find that the emitted energy density per unit photon energy is given by:

$$\frac{du_\gamma }{dϵ_\gamma }(ϵ)=\frac{1}{2\zeta }\frac{dn_e}{dE_e}\left(\sqrt{ϵ/\zeta }\right)=\frac{A}{2}\zeta ^{s/21}ϵ^{s/2}\text{ .}$$
(16)
For a strong shock, where the electron number density satisfies $`dn_e(E)/dE_eAE^2`$, this reduces to:

$$ϵ\frac{du_\gamma }{dϵ_\gamma }(ϵ)=A/2\text{ .}$$
(17)
For weak shocks, the approximation in equation (16) introduces a small numerical error ($`2`$), because the temporal evolution of the electron distribution affects the normalization of the resulting spectrum.

Determining the present-day photon spectrum requires integration over different shocks at various times. This, in turn, is determined by the redshift-dependence of the IGM parameters, and requires the use of non-linear structure formation theory. We note that the contribution of a shock that developed at redshift $`z`$ to the observed value of $`ϵ[du_\gamma (ϵ)/dϵ_\gamma ]`$, is roughly proportional to $`(1+z)^4`$. Thus, the main contribution to the spectrum observed today, especially at high photon energies, is expected to originate from recent shocks at $`z0.5`$. Assuming that a fraction $`f_{\text{sh}}`$ of the baryons in the universe were recently shocked to a temperature $`T10^7\text{ K}`$, when $`\gamma _{\text{max}}`$ already approached $`10^7`$, we find from equation (14) and equation (17):

$$ϵ\frac{du_\gamma ^{\text{IC}}}{dϵ_\gamma }(ϵ)6.4\times 10^7\text{ }T_7\left(\frac{\xi _ef_{\text{sh}}}{0.05}\right)\left(\frac{\mathrm{\Omega }_Bh_{67}^2}{0.04}\right)\text{ eV cm}^3\text{ .}$$
(18)
The inverse-Compton flux per decade in photon energy thus becomes:

$$ϵ^2\frac{dJ_\gamma ^{\text{IC}}}{dϵ_\gamma }(ϵ)1.5\text{ }T_7\left(\frac{\xi _ef_{\text{sh}}}{0.05}\right)\left(\frac{\mathrm{\Omega }_Bh_{67}^2}{0.04}\right)\text{ keV cm}^2\text{ s}^1\text{ sr}^1\text{ .}$$
(19)
This result is in excellent agreement with observational data, when supplemented by the contribution from unresolved point sources Loeb & Waxman (2000). One can employ the Press-Schechter mass function in order to carry out the time-integration over the radiation resulting from intergalactic shocks at various epochs. This approach gives, for slightly different cosmological parameters ($`\mathrm{\Omega }_\mathrm{\Lambda }=0.65`$, $`\mathrm{\Omega }_B=0.05`$, $`h=0.7`$), a similar result Loeb & Waxman (2000):

$$ϵ^2\frac{dJ_\gamma ^{\text{IC,PS}}}{dϵ_\gamma }(ϵ)1.5\left(\frac{\xi _e}{0.05}\right)\text{ keV cm}^2\text{ s}^1\text{ sr}^1\text{ .}$$
(20)
The synchrotron radiation emitted by the accelerated electrons may now be estimated. A relativistic electron with a Lorentz factor $`\gamma _e`$, gyrating in a magnetic field of amplitude $`B`$, emits photons with characteristic energy $`ϵ_\gamma [(\mathrm{}eB)/(m_ec)]\gamma _e^2`$. Synchrotron radiation is thus expected in photon energies ranging from $`10^{12}\text{ eV}`$ and up to $`1\text{ eV}`$. The emitted flux can be crudely estimated as:

$`ϵ^2{\displaystyle \frac{dJ_\gamma ^{\text{syn}}}{dϵ_\gamma }}(ϵ)`$ $``$ $`ϵ^2{\displaystyle \frac{dJ_\gamma ^{\text{IC}}}{dϵ_\gamma }}(ϵ){\displaystyle \frac{u_B}{u_{\text{CMB}}}}`$ (21)
$``$ $`1.5\text{ }T_7B_7^2\left({\displaystyle \frac{\xi _ef_{\text{sh}}}{0.05}}\right)\left({\displaystyle \frac{\mathrm{\Omega }_Bh_{67}^2}{0.04}}\right)\text{ eV cm}^2\text{ s}^1\text{ sr}^1\text{ ,}`$
whereas a more careful approach is to assume that the energy density of the post-shock magnetic field constitutes a fraction $`\xi _B`$ of the shock-induced thermal energy. This assumption, combined with the Press-Schechter mass function, leads to a similar (somewhat higher) flux:

$$ϵ^2\frac{dJ_\gamma ^{\text{syn,PS}}}{dϵ_\gamma }(ϵ)3\left(\frac{\xi _e}{0.05}\right)\left(\frac{\xi _B}{0.01}\right)\text{ eV cm}^2\text{ s}^1\text{ sr}^1\text{ ,}$$
(22)
where a typical value $`\xi _B=0.01`$ produces magnetic fields consistent with observations of the outer halos of X-ray clusters Waxman & Loeb (2000).

#### 2.3.2 Radiation From Forming Objects, Angular Statistics

We start with the general case of an astrophysical object formed recently, and calculate its inverse-Compton luminosity. We assume that the object accretes mass at a rate $`\dot{M}`$, heating it up to a temperature $`T`$ by strong shocks. The baryon number accretion rate is thus given by $`\dot{M}(\mathrm{\Omega }_B/\mathrm{\Omega }_M)m^1`$. We assume, as previously, that a fraction $`\xi _e`$ of the post-shock thermal energy is converted into relativistic electrons. The maximal electron Lorentz factor $`\gamma _{\text{max}}`$ is determined from equation (12), and depends on the assumed magnetic field. Electrons with cooling time longer than the virialization time of the object, $`t_{\text{vir}}`$, will not have sufficient time to radiate a significant fraction of their energy. This introduces a low-energy cutoff in the emitted spectrum, at photon energy $`ϵ_{\text{min}}=(4/3)\gamma _{\text{min}}^2ϵ_{\text{CMB}}`$, where $`\gamma _{\text{min}}2\times 10^{12}\text{ yr}/t_{\text{vir}}`$. This effect reduces the total luminosity by a factor $`\mathrm{ln}\left(\gamma _{\text{max}}/\gamma _{\text{min}}\right)/\mathrm{ln}\gamma _{\text{max}}`$. Combining these considerations, we find that the total and specific luminosities of the object are given by:

$$L^{\text{IC}}=\dot{M}\frac{\mathrm{\Omega }_B}{\mathrm{\Omega }_M}m^1\text{ }\xi _e\frac{3}{2}k_BT\text{ }\frac{\mathrm{ln}\left(\gamma _{\text{max}}/\gamma _{\text{min}}\right)}{\mathrm{ln}\gamma _{\text{max}}}\text{ ,}$$
(23)
and:

$$ϵ_\gamma \frac{dL^{\text{IC}}}{dϵ_\gamma }=\frac{L^{\text{IC}}}{2\mathrm{ln}\left(\gamma _{\text{max}}/\gamma _{\text{min}}\right)}\text{for energies}1\left(\frac{\gamma _{\text{min}}}{10^3}\right)^2\text{ keV}ϵ_\gamma 0.1\left(\frac{\gamma _{\text{max}}}{10^7}\right)^2\text{ TeV}\text{ .}$$
(24)
Structure formation at the current epoch produces dense, large-scale sheets and filaments, entangled in space. The densest regions appear at the intersections of these features, around the locations of galaxy clusters. As an example, we consider such a galaxy cluster. The typical gas mass involved in the formation of a cluster is $`M_{\text{gas}}10^{14}M_{}`$, accreting over a virialization time $`t_{\text{vir}}10^9\text{ yr}`$ and shocked to a temperature $`T_{\text{gas}}10^7\text{ K}`$. A magnetic field of amplitude $`B0.1\mu \text{G}`$, as suggested by observations, gives $`\gamma _{\text{min}}2000`$ and $`\gamma _{\text{max}}3.3\times 10^7`$. Whence,

$$L_{\text{cluster}}^{\text{IC}}2.2\times 10^{45}\left(\frac{\xi _e}{0.05}\right)\left(\frac{k_BT_{\text{gas}}}{5\text{ keV}}\right)\left(\frac{M_{\text{gas}}}{10^{14}M_{}}\right)\left(\frac{t_{\text{vir}}}{10^9\text{ yr}}\right)^1\text{ erg s}^1\text{ ,}$$
(25)
and:

$$ϵ_\gamma \frac{dL_{\text{cluster}}^{\text{IC}}}{dϵ_\gamma }0.05L_{\text{cluster}}^{\text{IC}}\text{for energies}5\text{ keV}ϵ_\gamma 1\text{ TeV}\text{ .}$$
(26)
A comprehensive discussion of particle acceleration in galaxy clusters and the resulting inverse-Compton radiation up to hard X-rays, may be found in Sarazin (1999).

The angular statistics of the radiation predicted by the model can be investigated by simplifying assumptions regarding the mass distribution. Waxman and Loeb (2000) calculated the angular correlation of the spectrum using the Press-Schechter mass function. On sub-degree angles they find fluctuations $`40\%`$ in the inverse-Compton spectrum and fluctuations $`>100\%`$ in the synchrotron spectrum.

## 3 Simulation

The previous section presented some order of magnitude estimates, leading from linear or non-linear structure formation theories to an EGRB spectrum, that roughly matches present-day observations. However, further progress along this route is severely compromised by the complexity and non-linearity of structure formation. This situation naturally calls for the use of cosmological simulations that incorporate hydrodynamical effects. For this project, we analyze a simulation of a $`\mathrm{\Lambda }\text{CDM}`$ cosmology. In the following, we describe the cosmological model and background theory of this simulation, and discuss some characteristics of the simulated universe and the predicted underlying structure.

### 3.1 Model and Theory

The cosmological model adopted in the simulation is the ’concordance’ $`\mathrm{\Lambda }\text{CDM}`$ model of Ostriker & Steinhardt (1995), presented in §2.1. An initial Harrison-Zel’dovich spectrum of fluctuations was assumed, and linear structure formation theory Eisenstein & Hu (1998) employed to predict the resulting spectrum at a certain redshift $`z_{\text{start}}`$, providing initial conditions for the simulation. The various parameters of the cosmological model are summarized in Table 1 (see e.g. Springel, White & Hernquist 2001b).

The simulation was performed using the TreeSPH code GADGET (Springel, Yoshida, & White 2001a). This code, designed for both serial and parallel computers, is intended for collisionless and gas-dynamical cosmological simulations.

Dark matter is modeled by the code as a self-gravitating collisionless fluid, represented by a large number $`N_{\text{DM}}`$ of particles. Newton’s gravitational potential is modified by introducing a spline softening at small separations, equivalent to treating each particle as a normalized spline mass distribution, instead of as a point mass (see, e.g. Hernquist & Katz 1989). This modification is necessary in order to suppress large angle scattering in two-body collisions. Poisson’s equation is solved using a hierarchical tree algorithm Barnes & Hut (1986).

The baryonic component of the universe is modeled by the simulation as an ideal gas, with an adiabatic index $`\mathrm{\Gamma }=5/3`$, and is represented by a large number of particles, similar to the dark matter. The gas is governed by the continuity equation:

$$\frac{d\rho }{dt}+\rho 𝐯=0\text{ ,}$$
(27)
the Euler (momentum) equation:

$$\rho \frac{d𝐯}{dt}=P\rho \mathrm{\Phi }+\overline{\mathrm{\Sigma }}\text{ ,}$$
(28)
and the thermal energy equation (the first law of thermodynamics):

$$\frac{du}{dt}=\frac{P}{\rho }𝐯\frac{1}{\rho }\overline{\mathrm{\Sigma }}_{ij}_i𝐯_j\frac{\mathrm{\Lambda }(u,\rho )}{\rho }\text{ ,}$$
(29)
where $`P`$, $`\rho `$, $`u`$ and $`𝐯`$ are the pressure, mass density, thermal energy per unit mass and velocity of a gas element; $`\mathrm{\Phi }`$ is the gravitational potential, found from Poisson’s equation:

$$^2\mathrm{\Phi }(r,t)=4\pi G\rho (r,t)\text{ ;}$$
(30)
and we have used Lagrangian derivatives:

$$\frac{d}{dt}=\frac{}{t}+𝐯\text{ .}$$
(31)
An artificial viscosity Steinmetz (1996) term $`\overline{\mathrm{\Sigma }}`$ was added to the above equations, in order to capture shocks. A cooling function $`\mathrm{\Lambda }(u,\rho )`$ may be included in the energy equation, to incorporate radiative cooling processes.

These hydrodynamical equations are solved using the SPH - Smoothed Particle Hydrodynamics - technique (Lucy 1977; Gingold & Monaghan 1977; Monaghan 1992), well suited for non-axisymmetric three-dimensional astrophysical problems. The version of the code used for the simulation employs adaptive SPH smoothing lengths, such that each particle has a constant number $`N_s`$ of neighbors lying within its kernel. The hydrodynamical equations and their solutions in the formalism of this version of SPH appear in Springel et al. (2001a). We note that the simulation analyzed in this paper employed a formulation of SPH which does not rigorously conserve entropy in adiabatic parts of the flow (for a discussion, see Hernquist 1993, Springel & Hernquist 2002a).

The code simulates the evolution of dark matter and gas (hereafter: SPH) particles within a finite simulation box of comoving size $`L_{\text{com}}^3`$, using periodic boundary conditions. Due to technical reasons, to be discussed in the next section, we have chosen to employ a simulation that did *not* include radiative cooling processes. Radiative cooling, dominated by Bremsstrahlung and line emission at the epoch of interest, is sub-dominant relative to inverse-Compton emission from relativistic electrons in all but the densest regions within the cores of clusters where the baryon number density exceeds $`10^3\text{ cm}^3`$, and in cold ($`T10^4\text{ K}`$), dilute regions. We discuss the significance of radiative cooling towards the end of this section. The simulation neglected feedback from star formation and supernovae. The code parameters used in the simulation are presented in Table 2.

### 3.2 Resulting Structure

The simulation chosen for the project, as most cosmological simulations, predicts a highly structured IGM in the present-day universe, as illustrated in Figure 1. Several hot and dense objects are spread throughout space, separated by large empty voids. The largest of these structures are sheet-like objects (pancakes) and filaments, entangled in space, some as long as $`100\text{ Mpc}`$. At the intersection of the filaments appear denser, spheroidal objects, with typical size $`1\text{ Mpc}`$, where galaxy clusters are formed. Simulations predict that the pancake-resembling objects are the first large structures to appear. At the intersections of these pancakes, filaments emerge, which in turn intersect and drain mass into spheroidal objects. At the large scales discussed, dark matter and baryons essentially trace the same structure.

The objects described above are significantly denser than their surroundings. Thus, the last stages of their evolution are dominated by non-linear processes. These processes are mainly gravitational, though other physical processes, such as non-linear gas dynamics and radiative cooling (when included), influence these stages of structure evolution as well.

It is instructive to examine the thermodynamic evolution of the universe, by applying some averaging schemes to the temperature and density at different epochs, as shown in Figure 2. The mass averaged temperature is notably higher than the volume averaged temperature, because hot regions are typically much denser, leaving most of the volume with a dilute, cooler gas (e.g. Davé et al. 1999). The present-day averaged temperature reached by the simulation, $`T(z=0)=3.9\times 10^6\text{ K}`$ (see Figure 2), is significantly lower than the order-of-magnitude estimate presented earlier. (Note that the approximation of eq.  is in rough agreement with the mass averaged temperature, although deviating towards higher temperatures at recent times.) Furthermore, the simulation result for the present-day mass averaged temperature is lower than that predicted by Cen & Ostriker (1999) by a factor $`2.5`$, but in good agreement with the results of e.g. Refregier et al. (2000), Croft et al. (2001), and Refregier & Teyssier (2000). R. Cen (private communication) informs us that revised estimates for the temperature evolution of the Cen & Ostriker (1999) simulation, correcting for an error in the integration of the thermal energy equation, now actually lie below the other simulation results.

The volume averaged density is proportional to $`(1+z)^3`$, due to the expansion of the universe, whereas the approximately constant mass averaged density indicates that much of the mass has virialized. The temperature averaged density is higher than the volume averaged density, due to the above mentioned correlation between hot and dense regions. This (temperature averaged) density decreases in time, growing closer to the volume averaged density, as the dilute gas surrounding virialized structures heats up.

We examine the effect of radiative cooling by super-imposing cooling-time contours on the distribution of mass in the temperature-density phase-space (Figure 3, left). The figure suggests that radiative cooling has important effects on the present-day gas, residing in the most dense, hot regions and in the most cold, dilute regions. The cooling time in these volumes due to Bremsstrahlung and line emission, $`t_{\text{cool}}`$, is comparable to or smaller than the dynamical time, $`t_\text{H}H^1`$. Such regions, where $`t_{\text{cool}}<t_\text{H}`$, presently contain approximately $`25\%`$ of the gas mass, and a progressively larger mass fraction at higher redshifts (Figure 3, right). The dense and hot regions, where radiative cooling is significant, are associated with the cores of clusters, well beyond the accretion shock fronts of interest and should thus have little effect on our results, although the weaker merger and flow shocks may be affected. Radiative effects in the cold, dilute regions cools the collapsing gas and may result in slightly stronger accretion shocks. However, as these regions are cold to begin with, this effect should be unimportant.

## 4 Method

Next we discuss how cosmological simulations can be exploited to yield various predictions regarding the EGRB, according to the Loeb-Waxman model. Parts of the discussion are specific to the cosmological simulation chosen for the project, introduced in the previous section, and to the EGRB extraction schemes used. Predictions of shock-induced radiation are obtained from the simulation in three basic steps: First, the simulation is used to keep track of the forming structures and locate the shocks involved in their virialization. Second, relativistic electron populations are injected into the shocked gas, according to the predictions of the Fermi acceleration mechanism. And third, the radiation emitted by the relativistic electrons is calculated and integrated to give measurable EGRB features. The organization of this section corresponds to these three steps.

### 4.1 Extracting Shocks

Since shock waves are discontinuities in the thermodynamic quantities, locating shocks involves identifying large gradients in the properties of the gas. Gradients of quantities such as pressure, temperature or velocity, can be calculated directly from a simulation. In our case, this requires SPH kernel interpolation, as described in the previous section. Next, one can search for gradients larger than some imposed thresholds. These would imply the presence of shocks, stronger than a corresponding shock strength. The shock front could then be traced out, yielding the shock topology and the jump conditions across the shock, which can in turn be checked to see if they satisfy the Rankine-Hugoniot adiabat.

Such an approach was adopted by Miniati et al. (2000), using simulations based on a version of the total variation diminishing scheme (TVD), designed to capture strong shocks Ryu et al. (1993). They identified very strong accretion shocks with Mach numbers $`M`$ of up to $`10^3`$ around the filaments and halos of their simulations. These halos exhibit, in addition, weaker shocks ($`M310`$), originating from the merger of structures or from accretion shocks into old structures, already embedded within newly formed halos. Combined, these shocks lead to a complicated pattern of inter-winding shocks. Miniati et al. find that the jump conditions, imposed by a shock, are altered by adiabatic gravitational compression, especially when the shock is close to the core of a cluster. Since the limited simulation resolution renders these two effects inseparable, shock parameters were extracted under the assumption that gravitational effects are isothermal.

#### 4.1.1 Shock-Extraction Scheme

Since we are interested more in the shocked gas than in the topology of shocks, we find SPH simulations sufficient for our needs, although they do provide limited shock front resolution. This enables us to adopt a simpler approach for locating shocks, appropriate for Lagrangian simulations such as SPH. This method involves following a gas element throughout the simulation, searching for rapid increases in its *entropy*. In cosmological simulations, entropy changes of the gas should result only from shock waves or from cooling processes. This is the reason we have chosen, for simplicity, to work with a simulation in which radiative cooling is neglected, shown in the previous section to yield reliable shocks. In such a simulation, the entropy of the gas should only increase, and do so in discrete steps, decisively indicating the presence of shocks. Ideally, one could then identify all entropy ’jumps’ as shock waves, determining their space-time coordinates. The thermodynamic jump conditions across a shock could then be found, by comparing the thermodynamic parameters of the gas element before ($`T_i`$, $`\rho _i`$) and after ($`T_f`$, $`\rho _f`$) passing through the shock, directly, or using the entropy-change, $`\mathrm{\Delta }S=C_V\mathrm{ln}\left[(T_f/T_i)(\rho _f/\rho _i)^{(1\mathrm{\Gamma })}\right]`$, induced by the shock. Using the Rankine-Hugoniot adiabat, we find that the compression factor $`r`$ of a shock - related to its Mach number by equation (9) - can be determined from $`\mathrm{\Delta }S`$ according to:

$$\mathrm{\Delta }S_{}C_V^1\mathrm{\Delta }S=\mathrm{ln}\left[\frac{r(\mathrm{\Gamma }+1)(\mathrm{\Gamma }1)}{(\mathrm{\Gamma }+1)r(\mathrm{\Gamma }1)}r^\mathrm{\Gamma }\right]\text{ .}$$
(32)
In SPH simulations, one can choose the SPH particles *themselves* as the gas elements to be followed, since they are of constant mass and the simulation is Lagrangian. This technically simplifies our shock-locating scheme, as kernel interpolation becomes unnecessary. It is instructive to study the distribution of mass - i.e. SPH particles - that accumulated various entropy-changes during recent epochs. Figure 4 (left) gives this distribution in terms of the mass fraction $`f(\mathrm{\Delta }S_{}^{})=m_{\text{tot}}^1dm(\mathrm{\Delta }S_{}^{})/d\mathrm{\Delta }S_{}`$ of SPH particles that accumulated a normalized entropy change $`\mathrm{\Delta }S_{}^{}`$ between $`z=2`$ and $`z=0`$, where $`m_{\text{tot}}`$ is the total mass of SPH particles. First, note that a small fraction of the gas experienced *negative* entropy changes, most probably as a result of lack of strict entropy conservation in the absence of shocks, owing to the particular formulation of SPH used for this simulation (see, e.g. Springel & Hernquist 2002a). The entropy-difference distribution among such particles resembles a normal distribution, consistent with accumulating noise. The shape of the remaining distribution has an interesting structure, suggesting a pronounced population of ’weak’ shocks, with $`M310`$, and a large population of stronger shocks, with Mach numbers up to a few $`10^3`$, in agreement with Miniati et al. (2000).

In order to determine the location and temporal duration of shock waves, one must examine the simulated gas in sufficiently short time intervals. This is also necessary in order to identify multiple shocks suffered by an SPH particle, as the effect of consecutive shocks is not additive, although the entropy change they induce is. We chose to examine snapshots of the simulated gas, separated by the simulation-box light-crossing time, as illustrated in Appendix A. For the mass resolution of the simulation, the entropy-increase experienced by most shocked particles is spread over several of these snapshots, indicating that shorter time intervals are not required. With the time intervals that were chosen, the entropy-difference distribution between any two snapshots (Figure 4, right) is fairly constant throughout the period of interest. Based on this distribution, we impose a minimal threshold $`\mathrm{\Delta }S_{}^{\text{min}}`$, such that only particles experiencing a larger entropy-increase between two consecutive snapshots are considered shocked. A threshold value $`\mathrm{\Delta }S_{}^{\text{min}}=0.4`$ was found to be reasonable, considering numerical noise, although it entails the danger of eliminating some of the mildly weak shocks ($`M10`$), if an SPH particle requires many snapshots to cross such a shock. The fact that only $`10\%`$ of the gas suffers entropy-changes larger than $`\mathrm{\Delta }S_{}^{\text{min}}`$ simplifies the computational work considerably. Demanding that the entropy increases at least by $`n_s\mathrm{\Delta }S_{}^{\text{min}}`$ during $`n_s=2`$ consecutive snapshots, and dismissing shocked particles that did not reach a temperature $`T_{\text{min}}10^5\text{ K}`$ by the present time, further simplifies the analysis. Our final results do not strongly depend upon the choice of these threshold values near those specified above. Note that with finite resolution, SPH particles may pass through a shock in a period longer than the snapshot temporal resolution. The simulation must thus be traced down to $`z<0`$ values in order to capture all recent shocks.

The temperature jump across a shock determines the thermal energy (per particle) it adds to the gas, $`\mathrm{\Delta }U_{\text{th}}/N=(\mathrm{\Gamma }1)^1k_B(T_dT_u)`$, and is thus vital for the calculation of shock-induced radiation. In principle, this jump could be found directly, by comparing the temperature of the shocked SPH particle before and after a shock takes place. However, with finite simulation resolution, this approach leads to errors, due to significant adiabatic heating or cooling of the particle in parallel to the shock. In some extreme cases, we even found the shocked particle cooler than it was prior to the shock. (A similar problem when radiative cooling is active has been identified by e.g. Hutchings & Thomas  and Martel & Shapiro .) A partial remedy for this problem is to assume that the entropy induced by the shock - and thus the compression factor - was correctly determined so that the state of the shocked particle can be deduced from its state prior to the shock, or vice versa. Indeed, neither the pre-shock state of the particle nor its post-shock state are accurately known, due to the inaccurate shock timing and ongoing adiabatic processes, leading to some inevitable error. We chose to regard the temperature of the gas, before the shock was detected, as the true pre-shock temperature: $`T_uT_i`$. Hence, we take:

$$T_d=T_ir^{(\mathrm{\Gamma }1)}e^{\mathrm{\Delta }S_{}}\text{ }$$
(33)
as the post-shock temperature, where $`r`$ and $`\mathrm{\Delta }S_{}`$ are related by equation (32). This choice yields more realistic results than the alternative, $`T_d=T_f`$, as discussed in §4.1.2.

#### 4.1.2 Shock-Extraction Results

The shock-extraction scheme, described above, identifies shocks in $`40\%`$ of the simulated gas mass during recent epochs ($`z<2`$), where about $`30\%`$ of the baryons were shocked at least once (Figure 5, left). Most identified shocks are concentrated around the halos in the simulation, although filaments and sheets are traced out as well. In some cases, but not always, shock fronts can be identified (Figure 6).

The distribution of entropy-differences, accumulated by *identified shocks* (Figure 5, right), suggests that most of the strong shocks with Mach numbers in the range $`M>10^2`$ were located. It confirms, however, our suspicion of loosing most of the ’weak’ shock population ($`M310`$). In such shocks, the entropy-change suffered by a single particle is spread over many snapshots, and is thus eliminated by the imposed entropy-change threshold. Figure 5 (right) suggests that some of the very large accumulating entropy-increases, suffered by an SPH particle, result from multiple shocks. Separating the shocks in such cases results in a population of weaker shocks, artificially concentrated at the threshold shock strength ($`M4`$).

As a ’sanity check’ of the proposed shock-location and energy-evaluation schemes, we compare the thermal energy gain of the simulation to the evaluated energy, induced by the identified shocks. We find that the energy, injected by our shocks, constitutes a fraction $`f_{\text{sh}}2/3`$ of the total thermal energy gained by the simulated gas (see Table 3 and Figure 5, left). For comparison, in the estimates of §2.1 we took $`f_{\text{sh}}=1`$. Miniati et al. (2000) find, for a slightly different $`\mathrm{\Lambda }\text{CDM}`$ model ($`\mathrm{\Lambda }=0.55`$, $`h=0.6`$), that the energy passing through shocks with Mach numbers $`M>10^2`$ between redshift $`z=1.5`$ and today, approximately equals the present-day thermal energy of their simulation.

The above results were obtained by estimating the energy of an SPH particle, gained when passing through a shock, based on its temperature *before* passing through the shock (i.e. taking $`T_u=T_i`$). Estimating this energy gain based on the temperature *after* passing through the shock (i.e. taking $`T_d=T_f`$) results in shock energies lower, on average, by a factor $`4.5`$. This result is a combined effect of the low resolution of the simulation and adiabatic cooling, taking place in parallel to the shock.

We find the shock-extraction results described above reassuring, as they indicate that the proposed shock location-and-evaluation recipe captures most of the strong shocks, and measures their energy reasonably well. Better performance of this algorithm may be achieved by using simulations with a higher mass resolution and with a more accurate handling of adiabatic flows (e.g. Springel & Hernquist 2002a). It should be stressed that the elimination of weak shocks by imposing an entropy threshold for shock identification has only a small effect on the predicted diffuse $`\gamma `$-ray radiation, especially at high energies. First, the consistency check described above verifies that the energy injected by weak shocks is small compared to the energy injected by strong shocks: accounting for all weak shocks may increase the integrated flux at most by a factor of $`1.5`$. Second, the contribution of weak shocks to the $`\gamma `$-ray emission, especially at high photon energies, is small due to the steep distribution of the electrons accelerated by such shocks, reducing their $`\gamma `$-ray efficiency. At $`10\text{ GeV}`$, for example, shocks with $`M=10`$ have $`66\%`$ efficiency (compared to strong shocks), dropping below $`50\%`$ for $`M8`$. Weak shocks could, nonetheless, increase the $`\gamma `$-ray brightness of individual nearby structures such as galaxy clusters experiencing merger events, mainly in low photon energies, thereby slightly altering the appearance of their inner regions, and as a result may somewhat enhance the predicted source number counts.

### 4.2 Accelerated Electron Population

Here we rely on the arguments of §2.2, deriving expressions more general and more precise than the order of magnitude estimates presented earlier. These results are used to calculate the relativistic electron distributions produced by shocks, and their temporal evolution. Appropriate electron populations can thus be ’injected’ into the shocked SPH particles to simulate radiative processes.

#### 4.2.1 Production

We postulate that a fraction $`\xi _e`$ of the post-shock density of the thermal energy induced by a shock, $`\mathrm{\Delta }u_{\text{th}}=[1/(\mathrm{\Gamma }1)]n_dk_B(T_dT_u)`$, is converted into a power law distribution of relativistic electrons, $`dn_e(E)/dE_e=AE^s`$. The power law index is related to the compression factor by $`s=(r+2)/(r1)`$, where $`\mathrm{\Gamma }=5/3`$ was used here and is assumed henceforth. Imposing a distribution cutoff at the maximal energy attainable by the electrons, $`E_{\text{max}}\gamma _{\text{max}}m_ec^2`$, fixes the normalization:

$$A=\xi _e\mathrm{\Delta }u_{\text{th}}\times \{\begin{array}{cc}1/\mathrm{ln}\gamma _{\text{max}}\hfill & \text{for }s=2\text{ ;}\hfill \\ (s2)(m_ec^2)^{s2}/\left[1\gamma _{\text{max}}^{(s2)}\right]\hfill & \text{for }s>2\text{ .}\hfill \end{array}$$
(34)
We estimate the maximal Lorentz factor of the electrons by equating their acceleration e-folding time to the cooling time due to inverse-Compton scattering, yielding:

$$\gamma _{\text{max}}^2=\frac{3(\mathrm{\Gamma }1)}{8}\frac{eB_d}{mc^2\sigma _Tu_{\text{CMB}}}\left(k_BT_u+r\frac{\mathrm{\Gamma }+1}{\mathrm{\Gamma }1}k_BT_d\right)\text{ .}$$
(35)
The average mass per particle is given by

$$m=\frac{4m_p}{3\chi +1+4\chi _e}0.59m_p\text{ ,}$$
(36)
where $`\chi `$ is the hydrogen mass fraction, $`\chi _e`$ is the ionization ratio (per nucleon) and the second equality holds for a fully ionized gas, where $`\chi _e=(1+\chi )/2`$, with the hydrogen mass fraction of the simulation, $`\chi =0.76`$, assumed henceforth.

We parameterize the post-shock magnetic field, $`B_d`$, by assuming that a fraction $`\xi _B`$ of the post-shock thermal energy is transferred into the energy of this field, giving:

$$B_d=4.3\times 10^8\text{ }\left[T_{d,7}\left(\frac{\xi _B}{0.01}\right)\left(\frac{n}{10n_{\text{av}}(z)}\right)\right]^{1/2}(1+z)^{3/2}\text{ Gauss ,}$$
(37)
where $`n_{\text{av}}(z)`$ is the average particle number density at redshift $`z`$. Hence,

$$\gamma _{\text{max}}=1.1\times 10^7\text{ }(rT_{d,7}+T_{u,7}/4)^{1/2}\left[T_{d,7}\left(\frac{\xi _B}{0.01}\right)\left(\frac{n}{10n_{\text{av}}(z)}\right)\right]^{1/4}(1+z)^{5/4}\text{ ,}$$
(38)
which for strong shocks, where $`r4`$ and $`T_uT_d`$, reduces to:

$$\gamma _{\text{max}}^{(\text{strong})}=2.2\times 10^7\text{ }T_{d,7}^{3/4}\left[\left(\frac{\xi _B}{0.01}\right)\left(\frac{n}{10n_{\text{av}}(z)}\right)\right]^{1/4}(1+z)^{5/4}\text{ .}$$
(39)
#### 4.2.2 Evolution

After the relativistic electrons diffuse away from the shock, they cool, mainly by inverse-Compton scattering off CMB photons. The electron distribution - initially a power law - may thus evolve into a different form, as energetic electrons cool faster:

$$P_{\text{IC}}\frac{dE}{dt}=CE^2\text{ ,}$$
(40)
where we have defined $`C4u_{\text{CMB}}\sigma _T/3m_e^2c^3`$. This leads to a partial differential equation for the electron distribution:

$$\frac{}{t}n_\text{E}(E,t)=CE^2\frac{}{E}n_\text{E}(E,t)+2CEn_\text{E}(E,t)\text{ ,}$$
(41)
where the notation $`n_\text{E}n_e/E_e`$ is introduced. We ignore adiabatic cooling due to the Hubble expansion since the cooling timescales relevant for our calculated $`\gamma `$–ray spectrum are all much shorter than the Hubble time. For the initial (power-law) condition we find the solution:

$$n_\text{E}(E,t)=\{\begin{array}{cc}AE^s(1CEt)^{(s2)}\hfill & \text{for }E<\text{ min }\{E_{\text{max}},E_{\text{cool}}(t)\}\text{ ;}\hfill \\ 0\hfill & \text{for }E>\text{ min}\{E_{\text{max}},E_{\text{cool}}(t)\}\text{ ,}\hfill \end{array}$$
(42)
where we have defined $`E_{\text{cool}}(t)1/Ct`$. This solution implies that the electron distribution steepens (if $`s>2`$) towards $`E_{\text{cool}}(t)`$, beyond which all electrons have already cooled off. For strong shocks, where $`s2`$, the electron distribution is stationary for $`E<E_{\text{cool}}(t)`$, as equal numbers of electrons enter and leave each energy bin at any given time.

The inverse-Compton cooling time of a relativistic electron is a function of its Lorentz factor $`\gamma `$, approximately given by $`t_{\text{cool}}(\gamma )2.3\times 10^{12}\gamma ^1(1+z)^4\text{ yr}`$. This implies that electrons with $`\gamma <\gamma _{\text{min}}(z)`$, where $`\gamma _{\text{min}}`$ is calculated in equation (A7), do not have sufficient time to cool between redshift $`z`$ and today, and will not significantly contribute to the emitted radiation.

### 4.3 Emitted Spectrum

#### 4.3.1 Inverse Compton Emission

In the previous subsection we have demonstrated how the electron distribution, carried by each simulated SPH particle, can be calculated at any given time. Thus, it is straightforward, if somewhat lengthy, to calculate the inverse-Compton spectrum, emitted by each SPH particle. The emitted energy density per unit photon energy is obtained by integrating the single particle emission function, over both electron and photon distributions Rybicki & Lightman (1979):

$$\frac{dP_\gamma ^{\text{IC}}}{dϵ_\gamma }(ϵ,t)=\frac{3}{4}\sigma _Tc\frac{ϵ}{ϵ_0}\nu (ϵ_0)𝑑ϵ_0\gamma ^2\frac{dn_e}{d\gamma _e}(\gamma ,t)f\left(\frac{ϵ}{4\gamma ^2ϵ_0}\right)𝑑\gamma \text{ ,}$$
(43)
where $`\nu (ϵ_0)`$ is the specific number density of the incident photon field, assumed constant (in the time period of interest) and isotropic, $`f(x)`$ is the emission function, given by:

$$f(x)=2x\mathrm{ln}x+x+12x^2\text{ ,}$$
(44)
and we have assumed Thomson scattering in the electron rest frame: $`\gamma ϵ_0m_ec^2`$. The incident CMB photon field is approximately blackbody radiation, whence:

$$\nu (ϵ_0)=\frac{8\pi ϵ_0^2/h^3c^3}{\mathrm{exp}\left(ϵ_0/k_BT_{\text{CMB}}\right)1}\text{and}T_{\text{CMB}}2.73(1+z)\text{ K .}$$
(45)
The assumptions made above are all satisfied for the ultra-relativistic electrons ($`\gamma 100`$) of interest: The CMB photon field is isotropic and varies slowly in time, thus it can be treated as constant for the fast emission of such electrons. The Thomson condition is satisfied for $`\gamma _757(1+z)^1`$, rendering Klein-Nishina corrections unnecessary.

A straightforward integration of equation (43) may be replaced by a simpler approach. The typical lifetime of a shock, comparable to the Hubble time $`t_{\text{sh}}t_\text{H}H^1=1.5\times 10^{10}h(z)^1\text{ yr}`$, is much longer than the cooling time of ultra-relativistic electrons. Hence, a population of such electrons, constantly produced by a well resolved shock, will resemble a *steady* source of radiation. This enables us to approximate the emission from these electrons as *instantaneous*, i.e. we assume that all ultra-relativistic electrons with $`\gamma >\gamma _{\text{min}}(z)`$ loose all their energy to radiation, *immediately* after passing through the shock front. The error thus introduced is small, for electrons satisfying $`t_{\text{cool}}(\gamma )t_{\text{sh}}\gamma 155h(z)/(1+z)^4`$. Further justification for this approximation may be found in the temporal resolution of the snapshots used. The time between consecutive snapshots examined is determined by the simulation-box light-crossing time, of order $`t_{\text{res}}6.5\times 10^8(1+z)^1\text{ yr}`$. Thus, it is obviously pointless to follow the temporal evolution of electrons, satisfying $`t_{\text{cool}}(\gamma )<t_{\text{res}}\gamma >3.5\times 10^3(1+z)^3`$.

In the ’instantaneous emission’ approximation described above, we are interested only in the *time-integrated* specific power, emitted by the accelerated electrons:

$$\frac{du_\gamma ^{\text{IC}}}{dϵ_\gamma }(ϵ)=_0^{t_{\text{sh}}}\frac{dP_\gamma ^{\text{IC}}}{dϵ_\gamma }(ϵ,t)𝑑t\text{ .}$$
(46)
Since the only time-dependence in the emitted power of equation (43) appears in the distribution of electrons, $`dn_e(\gamma ,t)/d\gamma _e`$, we may change the order of integration and first calculate the time-integrated electron distribution:

$$\frac{dn_e^{\text{int}}}{dE_e}(E)=_0^{t_{\text{sh}}}n_\text{E}(E,t)𝑑t\frac{A}{C(s1)}E^{(s+1)}\text{ ,}$$
(47)
where the second equality is obtained by taking the upper integration limit to infinity or to $`\text{max}\{t_{\text{sh}},t_{\text{cool}}(E)\}`$, which is well justified for ultra-relativistic electrons. Integrating the power emitted by the *evolving* electron distribution, as in equation (43), may thus be replaced by calculating the *instantaneous* emission from the *integrated* electron distribution of equation (47).

Since the integrated electron distribution is a simple power law, we may use the well known result for inverse-Compton emission of a power-law electron distribution, scattering a blackbody photon field Rybicki & Lightman (1979). This leads to:

$`{\displaystyle \frac{du_\gamma ^{\text{IC}}}{dϵ_\gamma }}(ϵ)`$ $`=`$ $`Q(s)A\left[k_BT_{\text{CMB}}/\left(m_ec^2\right)^2\right]^{s/21}ϵ^{s/2}`$ (48)
$`=`$ $`Q(s)Aϵ_{\text{IC}}^{(s/21)}(1+z)^{s/21}ϵ^{s/2}\text{ ,}`$
where $`Q(s)`$ is a numerical coefficient, defined by:

$$Q(s)\frac{135}{2\pi ^4}2^s\frac{s^2+6s+16}{(s+4)^2(s+6)(s+2)}\mathrm{\Gamma }\left(3+\frac{s}{2}\right)\zeta \left(3+\frac{s}{2}\right)\text{ ,}$$
(49)
$`\zeta (s)`$ is the Riemann zeta function, and we defined $`ϵ_{\text{IC}}\left(m_ec^2\right)^2/\left[k_BT_{\text{CMB}}(z=0)\right]=1.1\text{ PeV}`$. This approximate result holds for emitted photon energies $`ϵ`$ in the range $`4\gamma _{\text{min}}^2ϵ_0ϵ4\gamma _{\text{max}}^2ϵ_0`$. Writing this as:

$$3.6\left(\frac{\gamma _{\text{min}}}{10^3}\right)^2\text{ keV}ϵ3.2\left(\frac{\gamma _{\text{max}}}{3\times 10^7}\right)^2\text{ TeV}\text{ ,}$$
(50)
indicates that the approximation holds for the entire relevant $`\gamma `$-ray regime.

It is interesting to compare the above results to the approximation used in §2.3.1, where only the first photon emitted by each electron was treated, thus neglecting the evolution of the electron distribution. The approximate result of equation (16) agrees with the last equality in equation (48), up to a function of the power index $`s`$. The two would become identical, if one was to replace $`Q(s)`$ in the latter by $`1/2`$, and take $`ϵ_{\text{IC}}=0.22\text{ PeV}`$. For a strong shock, we indeed find that $`Q(s=2)=1/2`$, and precisely recover equation (17) - the early estimate for a strong shock. This is not surprising, recalling that the electron distribution for $`s=2`$ is mostly stationary, rendering evolutionary effects insignificant. The numerical mismatch between equation (16) and equation (48) varies as function of $`s`$, but remains of order $`2`$ for relevant values of $`s`$ ($`s8`$).

#### 4.3.2 Space-Time Integration

Two different space-time integration schemes were used. The first, *line of sight* integration, calculates the specific intensity $`I(ϵ)`$ at a given direction $`(\theta ,\varphi )`$ on the sky, where $`\theta `$ is its angle with the $`z`$-axis and $`\varphi `$ is the angle of its projection on the $`x`$-$`y`$ plane with the $`x`$-axis. Such a scheme may be used to study the statistics of the predicted radiation, by calculating moments of the specific intensity: the mean intensity $`J(ϵ)=I(ϵ)`$, the correlation of intensities at various angular separations, etc. The second scheme provides *direction bin* integration, calculating the flux predicted within a small region $`\sigma =(\theta _1,\theta _2;\varphi _1,\varphi _2)`$ on the sky, defined by $`\theta _1\theta <\theta _2`$ and $`\varphi _1\varphi <\varphi _2`$. The flux, related to the specific intensity by $`F(ϵ)=_\sigma I(ϵ)\mathrm{cos}\theta d\mathrm{\Omega }`$, is required for sky maps and for source identification and number counts.

Both integration schemes involve simulating an integration volume (IV, for short), backtracking the propagation of radiation towards an imaginary observer. The IV starts at an arbitrarily chosen observer location today, and is moved away from the observer and back in time, at the speed of light. The radiation from sources (i.e. shocked SPH particles) within the IV is calculated and summed up, to yield the total radiation detected by the observer today. Since the simulation-box light-crossing time is much shorter than the Hubble time, periodic boundary conditions must be used to simulate the motion of the IV.

We treat the crossing of a shock front by an SPH particle as a single source of radiation for our integration scheme and call it, for short, a shock. Let such a shock, with specific luminosity $`L_{\text{sh}}(ϵ)`$ and redshift $`z`$, lie at coordinates $`(r,\theta _0,\varphi _0)`$ of a comoving coordinate system centered upon the observer. The contribution of this shock to the specific flux $`F_{\text{sh}}(ϵ)`$, detected by the observer today, satisfies:

$$ϵF_{\text{sh}}(ϵ)=\frac{ϵ^{}L_{\text{sh}}(ϵ^{})}{4\pi d_L(z)^2}\text{ ,}$$
(51)
where $`d_L(z)=rR_0(1+z)`$ is the luminosity distance between source and observer, $`R_0`$ is the present-day value of the cosmological scale factor, and we have defined $`ϵ^{}(1+z)ϵ`$. With the spatial and temporal resolution used, many ($`10`$) shocks fall within the IV in each time step. This enables us to approximate the shocks as point sources, such that the contribution of a shock with space-time coordinates $`(𝐫,t)`$ to the detected flux is determined by:

$$\mathrm{\Delta }\left[F(ϵ,\sigma )\right]=\{\begin{array}{cc}F_{\text{sh}}(ϵ)\hfill & \text{if }(𝐫,t)\text{ IV}(\sigma )\text{ ;}\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$
(52)
The specific intensity contributed by the shock is calculated by approximating the latter as a sphere of uniform specific brightness $`B_{\text{sh}}(ϵ)`$. The proper radius of this sphere is determined by the mass density $`\rho (z)`$ of the SPH particle, according to $`d_{\text{pr}}(z)=\left[3M_{\text{SPH}}/4\pi \rho (z)\right]^{1/3}`$, where $`M_{\text{SPH}}`$ is the mass of an SPH particle. The brightness of the sphere is related to its detected flux by:

$$F_{\text{sh}}(ϵ)=\pi B_{\text{sh}}(ϵ)\left[\frac{d_{\text{pr}}(z)}{d_A(z)}\right]^2\text{ ,}$$
(53)
where $`d_A(z)=rR_0/(1+z)`$ is the angular diameter distance between the source and the observer. Equating this result and equation (51) gives:

$$ϵB_{\text{sh}}(ϵ)=\frac{ϵ^{}L_{\text{sh}}(ϵ^{})}{4\pi ^2d_{\text{pr}}(z)^2(1+z)^4}\text{ ,}$$
(54)
and the contribution of the source to the specific intensity at direction $`(\theta ,\varphi )`$ is given by:

$$\mathrm{\Delta }\left[I(ϵ,\theta ,\varphi )\right]=\{\begin{array}{cc}B_{\text{sh}}(ϵ)\hfill & \text{if the ray }(\theta ,\varphi )\text{ intersects the sphere;}\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$
(55)
We approximate the specific luminosity of a shock by:

$$L_{\text{sh}}(ϵ)\frac{1}{\mathrm{\Delta }t_{\text{sh}}}\left[\frac{M_{\text{SPH}}}{\rho (z)}\right]\frac{du_\gamma }{dϵ_\gamma }(ϵ)\text{ ,}$$
(56)
where $`\mathrm{\Delta }t_{\text{sh}}t_ft_i`$ is the period of time during which the shock was detected. Since our shock locating-and-timing scheme is limited by the temporal resolution of the snapshots, $`\mathrm{\Delta }t_{\text{sh}}`$ tends to overshoot the true duration of the shock. This error is partially compensated by the integration schemes, as they are limited by the same snapshot resolution: although the approximate luminosity will tend to undershoot its true value, the probability of including the shocked particle in the IV is enhanced proportionally.

## 5 Results

Next we present various features of the radiation, inverse-Compton emitted by the shock-accelerated electrons, according to the cosmological simulation studied. We begin by presenting the predicted spectrum and discussing its implications. Next, we present some maps of the simulated $`\gamma `$-ray sky and discuss the sources they are composed of. We conclude with source number-counts extracted from such maps.

### 5.1 Predicted Spectrum

The inverse-Compton spectrum, predicted by the simulation, extends from the optical band into the far $`\gamma `$-ray regime, with photon energies in the range $`1\text{ eV}<ϵ<10^4\text{ TeV}`$ (Figure 7). However, the high energy tail of this spectrum is modified by photon-photon pair-production interactions, rendering the universe non-transparent for $`ϵ1\text{ TeV}`$ photons. Photons of $`ϵ100\text{ TeV}`$ thus interact with CMB photons or, at much higher energies, with radio background photons. Lower energy ($`1\text{ TeV}<ϵ<100\text{ TeV}`$) photons interact mainly with the intergalactic infrared (IR) background. The $`e^+e^{}`$ pairs, produced in such interactions, may inverse-Compton scatter CMB or IR photons, thus initiating a cascade process Nikishov (1962); Gould & Schreder (1967); Protheroe & Stanev (1993). As a result, the $`1\text{ TeV}`$ spectrum is significantly attenuated, whereas at lower photon energies, piling-up of the cascading photons slightly enhances (flattens) the spectrum.

The resulting spectrum provides roughly equal energy flux per decade in photon energy, in the photon energy range $`10\text{ keV}ϵ1\text{ TeV}`$. This flux ranges from $`ϵ^2dJ_\gamma (ϵ)/dϵ_\gamma 160\text{ eV cm}^2\text{ s}^1\text{ sr}^1`$ at $`ϵ=10\text{ keV}`$, through $`100\text{ eV cm}^2\text{ s}^1\text{ sr}^1`$ at $`ϵ=1\text{ GeV}`$, and down to $`50\text{ eV cm}^2\text{ s}^1\text{ sr}^1`$ at $`ϵ=1\text{ TeV}`$. The number flux of photons in the energy range $`10\text{ keV}ϵ1\text{ GeV}`$ is well approximated by:

$$\frac{dJ_\gamma }{dϵ_\gamma }(ϵ)1.1\times 10^{14}\left(\frac{ϵ}{100\text{ MeV}}\right)^{2.04}\left(\text{eV cm}^2\text{ s sr}\right)^1\text{ ,}$$
(57)
whereas in the energy range $`10\text{ GeV}ϵ5\text{ TeV}`$ the spectrum is slightly steeper and may be approximated, before including the cascade process, by:

$$\frac{dJ_\gamma }{dϵ_\gamma }(ϵ)8.8\times 10^{19}\left(\frac{ϵ}{10\text{ GeV}}\right)^{2.13}\left(\text{eV cm}^2\text{ s sr}\right)^1\text{ .}$$
(58)
The slope of the resulting spectrum and the energy range over which it extends are in good agreement with the order of magnitude estimates, carried out in §2.3.1. The resulting flux, on the other hand, falls short of the early estimates of equation (19) and equation (20), by a factor $`10`$ (Figure 7).

The discrepancy between the above flux, calculated from the simulation, and the flux predicted in §2.3.1, is due to a combination of several factors. First, the temperature reached by the simulation is a factor of $`2.6`$ lower than the order-of-magnitude temperature $`10^7\text{ K}`$ assumed previously. The thermal energy gain in the epoch examined ($`0<z<2`$) is thus smaller than the thermal energy used in §2.3.1 by a factor of $`3`$. Second, the energy induced by the identified shocks constitutes a fraction $`f_{\text{sh}}2/3`$ of the thermal energy gain in the simulation, as opposed to the previous assumption $`f_{\text{sh}}=1`$. Third, the simulation identifies only approximately *half* of the shock-induced electron-energy in electrons, that had sufficient time to cool ($`\gamma >\gamma _{\text{min}}`$), in contrast to the factor $`1\mathrm{ln}(\gamma _{\text{min}})/\mathrm{ln}(\gamma _{\text{max}})2/3`$ of §2.3.1, indicating the presence of weak shocks and very recent shocks. Finally, in §2.3.1 we neglected the redshift of the emitted energy, caused by the expansion of the universe, assuming that most the detected radiation was emitted recently ($`z<0.5`$). However, the energy accumulated by the detected shocks increases almost steadily for $`z<2`$ (Figure 5, left), such that a significant fraction of the predicted flux originates at high redshifts: only $`40\%`$ of this flux was emitted at $`z<0.5`$. This effect reduces the predicted flux by a factor $`5/3`$. Summing up the effects listed above, found from Table 3, yields the discrepancy factor: $`3(3/2)(4/3)(5/3)=10`$.

### 5.2 Sky Maps, $`\gamma `$-ray Sources and Source Number-Counts

By integrating the predicted flux, arriving from different directions, one can construct maps of the simulated $`\gamma `$-ray sky. Maps of the entire sky (Figure 8) and of two selected regions (Figure 9) are presented in the following. They reveal a globally isotropic picture, with strong local fluctuations, e.g. flux variation over $`2`$ orders of magnitude at $`1^{}`$ resolution and over more than $`3`$ orders of magnitude at $`5^{}`$ resolution. Such maps of the predicted $`\gamma `$-ray sky can be directly compared to the maps produced by detectors such as EGRET Sreekumar et al. (1998) aboard the CGRO satellite, detectors on board GLAST, planned for launch in 2006, and Čherenkov telescopes such as MAGIC, VERITAS 333See http://veritas.sao.arisona.edu/ and HESS 444See http://mpi-hd.mpg.de/hfm/HESS/HESS.html. Since such devices have finite angular resolution, the simulated maps must be convolved with appropriate filter functions, as demonstrated in Figure 9. The relevant parameters of EGRET, GLAST and MAGIC are summarized in Table 4.

Various $`\gamma `$-ray sources, of different shapes and sizes, are observed in the simulated maps of the $`\gamma `$-ray sky. Some of these sources exhibit irregular shapes, revealing the underlying structure of their emitter. Two typical examples, depicted in Figure 9, are strong emission from a semi-spherical source, and weak emission from a filamentary object. Some interesting features of the semi-spherical $`\gamma `$-ray source are shown in Figure 10. This source is the signature of a nearby rich galaxy cluster, lying $`53\text{ Mpc}`$ from the simulated observer, at a redshift of $`z0.012`$. The $`\gamma `$-ray emission from the cluster appears approximately as an elliptic ring, corresponding to a diameter $`510\text{ Mpc}`$, surrounding the location of the cluster and indicating the position of its accretion shock. The luminosity has two evident peaks along the circumference of the ring, approximately at the locations of its intersections with two large galaxy filaments. This suggests that the strong emission from these two regions is due to shocked gas, channeled by the filaments into the cluster region. Other galaxy filaments may be responsible for less luminous peaks of emission along the circumference of the ring, on the right hand side of Figure 10 images.

The simulated sky maps produced may be decomposed into discrete $`\gamma `$-ray sources, yielding number counts of sources above given photon energies (Figure 11). This is performed using a greedy source-identification algorithm devised for the purpose, discussed in Appendix B. The resulting number counts of sources above $`100\text{ MeV}`$ and above $`1\text{ GeV}`$ fall short of the EGRET detection threshold. Our simulation thus predicts that none of the $`60`$ unidentified extragalactic EGRET sources Özel & Thompson (1996) can be attributed to intergalactic shock-induced emission. The number-counts are especially sensitive to the fraction of shock-energy transferred to the relativistic electrons. For $`\xi _e=0.05`$ our simulation predicts a few dozen sources detectable by GLAST at both photon energy ranges, among them $`16`$ well-resolved sources above $`1\text{ GeV}`$ associated with different objects in the sky (see Figure 8). For $`\xi _e=0.02`$, the simulation predicts only one source detectable by GLAST, whereas $`\xi _e=0.10`$ yields more than a hundred GLAST sources, well resolved above $`1\text{ GeV}`$, yet too faint for EGRET detection.

Since our source-identification algorithm was devised for spherical sources, very large, nearby objects tend to be broken down and identified as several separate, weaker ones. This implies that when comparing our source number counts with analytical predictions Waxman & Loeb (2000); Totani & Kitayama (2000), one should impose a threshold on the angular size of the sources. As evident from Figure 11, Our number counts fall short of those of Waxman & Loeb (2000) by a factor of $`6`$. This may be explained, at least partially, by the same factors listed above that reduce the simulated diffuse $`\gamma `$-ray flux compared to the analytic calculation: the lower average temperature of the simulation suggests less hot, massive clusters, and combined with the thermal efficiency of the shocks and the fraction of electrons that had sufficient time to cool, both lower than assumed in the analytical predictions, may account for the discrepancy factor.

## 6 Discussion

We have studied the $`\gamma `$–ray emission from intergalactic shock waves, using a hydrodynamic $`\mathrm{\Lambda }\text{CDM}`$ cosmological simulation, according to the model proposed by Loeb & Waxman (2000). A simple approach for extracting large-scale shocks from a Lagrangian simulation, based on following entropy changes of single simulated particles, was employed. This method was shown to identify most of the strong (Mach number $`M100`$) accretion shocks in the simulation containing a large fraction ($`f_{\text{sh}}2/3`$) of the gas thermal-energy gain, although some of the weaker ($`M10`$) flow or merger shocks might be missed. Relativistic electron populations, Fermi-accelerated by the shocks, were injected into the shocked gas, assuming that a fraction $`\xi _e=0.05`$ of the shock thermal energy was transferred into relativistic electrons. The inverse-Compton radiation, emitted as these electrons scatter CMB photons, was then calculated and integrated over space and time, to yield the radiation detected by an imaginary observer.

This radiation was found to extend from the optical band and into the far $`\gamma `$-ray regime, containing roughly equal energy flux per decade in photon energy, $`ϵ^2[dJ_\gamma (ϵ)/dϵ_\gamma ]50160\text{ eV cm}^2\text{ s}^1\text{ sr}^1`$, in the photon energy range $`10\text{ keV}ϵ1\text{ TeV}`$. The calculated spectral *slope*, $`dJ/dϵϵ^{2.04}`$ for $`10\text{ keV}ϵ1\text{ GeV}`$ and $`dJ/dϵϵ^{2.13}`$ for $`10\text{ GeV}ϵ5\text{ TeV}`$ (before accounting for photon-photon pair-production cascade), is consistent with the predictions of Loeb & Waxman (2000) and with the EGRET observations. The *energy flux* constitutes $`10\%`$ of the EGRB flux, or up to $`15\%`$ of the EGRB flux after subtracting the expected contribution from unresolved point sources based on empirical modeling of the luminosity function of blazars Mukherjee & Chiang (1999).

The simulated $`\gamma `$-ray sources, mostly spheroidal halos associated with large-scale structure, fall short of the EGRET detection threshold, thus producing a diffuse background. Hence, the simulation predicts that none of the $`60`$ unidentified EGRET extragalactic $`\gamma `$-ray sources can be attributed to intergalactic shock-induced radiation. Several sources do fall within the detection range of GLAST, planned for launch in 2006, provided that $`\xi _e0.03`$. The number of detectable sources is sensitive to the fraction of shock-energy transferred to the relativistic electrons: For $`\xi _e=0.05`$ we find several dozen sources detectable by GLAST, among them $`16`$ well-resolved sources (see Figures 8, 11), whereas $`\xi _e=0.10`$ is still insufficient for EGRET detection, but leads to more than a hundred well-resolved GLAST $`\gamma `$-ray sources. Such sources could be identified as emission from intergalactic shock-accelerated electrons, by association with large-scale structure, a characteristic $`dJ/dϵϵ^2`$ spectrum, extending up to the far $`\gamma `$-ray regime, and a low cutoff at photon energy $`ϵ_{\text{min}}5\left(t/10^9\text{ yr}\right)^2\text{ keV}`$, caused by low-energy electrons not having sufficient time to cool during the lifetime $`t`$ of the source. Detection of intergalactic shock-induced $`\gamma `$-ray sources may be used to calibrate $`\xi _e`$, once their temperature has been determined. The shock-induced $`\gamma `$-ray emission, a first direct tracer of the warm-hot IGM ($`10^5\text{ K}T10^7\text{ K}`$), may thus be used to determine the baryon density in the present-day universe.

We have examined the $`\gamma `$-ray signature of a nearby simulated rich galaxy cluster, lying at a redshift $`z0.012`$. An elliptic emission ring of a diameter corresponding to $`510\text{ Mpc}`$ surrounds the cluster, tracing the location of its accretion shock (see Figures 9, 10). The luminosity peaks along the circumference of this ring, probably at the locations of intersections with galaxy filaments, channeling gas into the cluster region. Similar features should be detected by GLAST or by the MAGIC telescope in nearby rich galaxy clusters, at photon energies above $`10\text{ GeV}`$, and by VERITAS and HESS above $`200\text{ GeV}`$, if the background signal is assumed flat. A ring-like feature resembling the simulated emission should thus be detected in such galaxy clusters lying at redshifts $`z0.01`$, whereas the brightest regions of emission along the circumference of the ring could be detected in clusters with redshifts up to $`z0.025`$. These bright regions of emission lie above the background signal and should be detected at redshifts $`z0.015`$ even for $`\xi _e=0.02`$, whereas $`\xi _e=0.10`$ will enable the detection of the emission ring itself at similar redshifts. Detection of such a ring will enable a direct study of the cluster accretion shock, as well as establish the validity of our model. Detected bright regions along the circumference of such a ring may be tested to comply with data regarding a nearby galaxy filament and be used to measure the density and velocity of the untraced gas in such a filament.

As mentioned above, our algorithm for extracting shocks succeeded in capturing strong shocks, but missed some of the weaker ($`M10`$) ones. Accounting for these weak shocks could only have a minute effect on the predicted diffuse $`\gamma `$-ray radiation, since most ($`2/3`$) of the thermal energy has been identified in strong shocks and because the steep spectrum of electrons accelerated by weak shocks reduces their $`\gamma `$-ray efficiency. Weak shocks could, however, increase the $`\gamma `$-ray brightness of nearby structures such as galaxy clusters encompassing ongoing merger processes, mainly in low photon energies, slightly altering the appearance of their inner regions and enhancing the predicted source number counts.

In our work we have neglected the radiation resulting from the protons accelerated in the intergalactic shocks. These protons lead to $`\gamma `$-ray emission as they scatter off thermal protons to produce pions, most of the $`\gamma `$-ray flux originating from $`\pi ^0`$ decay. Even with a large proton-to-electron energy ratio $`\eta `$, the radiation resulting from the relativistic protons is, on average, negligible compared to the emission from electrons. Only in the densest ($`n>10^3\text{ cm}^3`$) regions in the cores of clusters, the relativistic protons may alter the $`\gamma `$-ray spectrum. This effect may make the cores of clusters even brighter $`\gamma `$-ray sources, provided that is not eliminated by diffusion of the relativistic protons out of the over-dense regions. For the values of $`\xi _e`$ and $`\eta `$ discussed above, the relativistic protons are dynamically significant and play an important role in shock dynamics and in particle acceleration. These effects, possibly detected in the radio spectra of some supernovae remnants, are small Reynolds & Ellison (1992), suggesting that the linear approach adopted in this paper is sufficiently adequate.

One can confirm that at least part of the EGRB is associated with the large-scale structure of the universe, by cross-correlating $`\gamma `$-ray maps with other sources that trace the same structure, such as galaxies, X-ray gas and the Sunyaev-Zel’dovich effect. The intergalactic shock-origin of $`\gamma `$-ray radiation may be verified by cross-correlating the EGRB with radio maps, partly attributed to synchrotron emission from the same shock-accelerated electrons. The synchrotron radio radiation is sub-dominant to the CMB, although its strong fluctuations on sub-degree scales were shown to contaminate CMB fluctuations at low ($`\nu 10\text{ GHz}`$) frequencies Waxman & Loeb (2000). Detection of a radio – $`\gamma `$-ray cross-correlation signal from nearby clusters may be used to measure the intergalactic magnetic field. The auto-correlation of inverse-Compton or synchrotron radiation, as well as their cross-correlation, may be calculated directly from our simulation.

A study of the shock-induced component of the extragalactic radiation, combined with a calculation of shock-induced emission according to cosmological simulations as described in this paper, may provide an important test of structure-formation models. Angular statistics of the $`\gamma `$-ray radiation, both on its own and cross-correlated with other sources, may be measured and compared to the predictions of such simulations. This will enable a study of the large-scale structure of the universe, and provide insight to the structure-formation process itself. We plan to apply the tools, developed in this research, to other cosmological simulations, including high-resolution simulations incorporating radiative cooling processes and supernovae feedback (e.g. Springel & Hernquist 2002b). Such research will permit an extensive study of shock-induced radiation, in the framework of various cosmological structure-formation models. Other physical phenomena, such as the acceleration of cosmic-rays by intergalactic shocks, may thus be investigated as well.

Order-of-magnitude estimates, presented in §2.3.1, predict an inverse-Compton flux higher by a factor of $`10`$ than obtained from the simulation we have studied. This is primarily attributable to the lower present-day temperature reached by the simulation, $`T_04\times 10^6\text{ K}`$, and to its slower heating-rate, which seem to comply with other recent simulations. Our simulation neglected radiative cooling, suffered from relatively low mass resolution ($`10^{11}M_{}`$), and did not include feedback from supernovae explosions, photo-heating the surrounding matter. The effect of radiative cooling on the accretion shocks is expected to be small, as cooling is efficient in hot, dense regions, beyond these shocks, and in cold, dilute regions, where the temperatures are low to begin with. Cooling should have a non-negligible effect on merger and flow shocks, but since these are much weaker than the accretion shocks, this should not significantly alter our results. In order to examine the sensitivity of our results to the low resolution of the simulation, a preliminary check of a similar adiabatic $`\mathrm{\Lambda }`$CDM simulation, with identical cosmological parameters but with eight times the mass resolution, was carried out, exhibiting no substantial difference from the results presented above. Outflows from galaxies (due to starburst or quasar activities) could also result in shock acceleration of relativistic electrons as they impact on the surrounding IGM. At the same time, preheating of the collapsing matter by, for example, supernovae, may lead to weaker, cooler shocks, as recently pointed out by Totani & Inoue (2001), and thus result in softer, weaker photon emission. However, Davé et al. (2001b) find that the luminosity-temperature relation for galaxy clusters, which originally motivated the addition of preheating to structure formation models, may in fact be attributed to radiative cooling. Therefore, preheating may not have played an important role in structure formation after all, and thus have little effect on our results.

*Note added in proof*: Recently, Scharf & Mukherjee (2002) have pointed out a possible association of $`\gamma `$-ray emission with the locations of galaxy clusters, by cross-correlating high Galactic latitude EGRET data with Abell clusters. The correlation signal, identified at a confidence level of $`3\sigma `$, is broadly consistent with the extended emission and with the source number counts discussed above, suggesting good prospects for the detection of extended $`\gamma `$-ray emission from clusters by future observations.

*Acknowledgments:* This work was supported in part by grants from the Israel-US BSF (BSF-9800343) and NSF (AST-9900877). UK & EW thank the Harvard-Smithsonian Center for Astrophysics for its kind hospitality during a period where part of this work was done. AL thanks the Minerva Einstein Center for partial support during his visit to the Weizmann Institute. EW is the incumbent of the Beracha foundation career development chair.

## Appendix A Integration Parameters

In order to perform the space-time integration schemes described in section 4.3.2, one must find the redshift dependence of several parameters: the time that elapsed between redshift $`z`$ and the present, the coordinate distance, traversed by a photon in that time, and the minimal Lorentz factor of electrons that lost a significant fraction of their energy since redshift $`z`$. Formulae for these parameters are provided below.

Using the FRW equations, we find that the time that elapsed between redshift $`z`$ and the present is given by:

$$\mathrm{\Delta }t=t_0t(z)=_{R(z)}^{R(z=0)}\frac{1}{\dot{R}}𝑑R=H_0^1_0^z\frac{1}{(1+z)h(z)}𝑑z\text{ ,}$$
(A1)
where $`h(z)H(z)/H_0`$. The coordinate distance traversed by a photon in that period is given by:

$$r=g\left[_{R(z)}^{R(z=0)}\frac{c}{\dot{R}R}𝑑R\right]=g\left[\frac{c}{R_0H_0}_0^zh(z)^1𝑑z\right]\text{ ,}$$
(A2)
where $`g(x)`$ is a function of the curvature $`k`$, defined by:

$$g(x)\{\begin{array}{cc}\mathrm{sin}(x)\hfill & \text{for }k=1\text{ ;}\hfill \\ x\hfill & \text{for }k=0\text{ ;}\hfill \\ \mathrm{sinh}(x)\hfill & \text{for }k=1\text{ .}\hfill \end{array}$$
(A3)
For our $`\mathrm{\Lambda }\text{CDM}`$ model, where $`\mathrm{\Omega }_\mathrm{\Lambda }+\mathrm{\Omega }_M=1`$ and $`h(z)=\sqrt{\mathrm{\Omega }_\mathrm{\Lambda }+(1+z)^3\mathrm{\Omega }_M}`$, the above equations are solved to give:

$$\mathrm{\Delta }t=\frac{2}{3H_0\sqrt{\mathrm{\Omega }_\mathrm{\Lambda }}}\mathrm{ln}\left[\frac{(1+\sqrt{\mathrm{\Omega }_\mathrm{\Lambda }})(1+z)^{3/2}}{\sqrt{\mathrm{\Omega }_\mathrm{\Lambda }}+\sqrt{\mathrm{\Omega }_M(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }}}\right]\text{ ,}$$
(A4)
and:

$$r=\frac{c}{3R_0H_0}\frac{(1)^{2/3}}{\mathrm{\Omega }_M^{1/3}\mathrm{\Omega }_\mathrm{\Lambda }^{1/6}}\text{ }\beta [\frac{1}{\mathrm{\Omega }_\mathrm{\Lambda }},1+\frac{\mathrm{\Omega }_M}{\mathrm{\Omega }_\mathrm{\Lambda }}\left(1+z\right)^3,\frac{1}{2},\frac{1}{3}]\text{ ,}$$
(A5)
where $`\beta `$ is the incomplete beta function, defined by:

$$\beta [t_1,t_2,a,b]_{t_1}^{t_2}t^{a1}\left(1t\right)^{b1}𝑑t\text{ .}$$
(A6)
In a similar fashion, we find the redshift dependence of $`\gamma _{\text{min}}(z,p)`$, the minimal Lorentz factor of an electron, that looses a fraction $`p`$ of its energy between redshift $`z`$ and today by inverse-Compton scattering. This parameter, which determines the low cutoff of the resulting spectrum, is given by:

$$\gamma _{\text{min}}(z,p)=\frac{p}{(1p)}\frac{3m_ecH_0}{4\sigma _Tu_{\text{CMB}}(z=0)}\left[_0^z\frac{(1+z)^3}{h(z)}𝑑z\right]^1\text{ ,}$$
(A7)
which for our $`\mathrm{\Lambda }\text{CDM}`$ model can be solved to give:

$`\gamma _{\text{min}}(z,p)`$ $`=`$ $`400h_{67}{\displaystyle \frac{p}{(1p)}}\mathrm{\Omega }_M/\{(1+z)h(z)1+\sqrt{\mathrm{\Omega }_\mathrm{\Lambda }}F\left({\displaystyle \frac{\mathrm{\Omega }_M}{\mathrm{\Omega }_\mathrm{\Lambda }}}\right)`$ (A8)
$``$ $`(1+z)\mathrm{\Omega }_\mathrm{\Lambda }F\left[(1+z)^3{\displaystyle \frac{\mathrm{\Omega }_M}{\mathrm{\Omega }_\mathrm{\Lambda }}}\right]\}\text{ ,}`$
where $`F(x)`$ is the hyper-geometric function $`{}_{2}{}^{}F_{1}^{}(\frac{1}{3},\frac{1}{2},\frac{4}{3},x)`$, defined by $`{}_{2}{}^{}F_{1}^{}(\frac{1}{3},\frac{1}{2},\frac{4}{3},x)_{k=0}^{\mathrm{}}\left(\frac{1}{3}\right)_k\left(\frac{1}{2}\right)_k(x)^k/\left(\frac{4}{3}\right)_kk!`$.

The redshift dependence of the parameters calculated above, also demonstrating the criteria for snapshot selection (simulation-box light-crossing time), is depicted in Figure 12.

## Appendix B Source Identification Algorithm

In order to count discrete $`\gamma `$-ray sources, predicted by the simulation, we have devised a greedy algorithm for extracting point sources from a sky map. This algorithm is used to isolate such sources, determine their location and flux, and measure their angular size. The algorithm operates iteratively, identifying sources and removing them from the map, until the contrast of the remaining map is smaller than some imposed threshold (of order $`2`$). At each iteration, the algorithm searches for the brightest point in the sky, and attempts to identify the brightest source associated with it. Our candidate for the source is assumed to have its kernel positioned about this bright point. We determine the area around the kernel, associated with the source by examining concentric rings about the kernel and adding them to the ’growing’ source. This accumulation process lasts as long as the flux within such a ring is sufficiently higher than the background, throughout its circumference, and yet monotonically decreasing, at least at one point along the circumference. Once the entire source has been identified, the flux associated with it is carefully removed from the map. The flux subtracted from each ring is determined according to the faintest point along its circumference, in order to prevent upsetting adjacent sources.

The source number counts produced by the algorithm were analyzed in section 5.2, with respect to the experimental devices EGRET and GLAST. In order to do so, algorithm parameters were tuned to mimic the source identification process, used to analyze the experimental data. Hence, the angular resolution of the analyzed sky maps was set to match the point-source resolution of the experimental devices, and the angular size of each source was compared to the photon-position error of each experiment, as illustrated in Figure 11.

Figure 13 demonstrates the results of the algorithm, when applied to two regions in the simulated $`\gamma `$-ray sky. The sky maps displayed in this paper were produced by approximating the source, associated with each SPH particle, as a rectangle in $`\theta \varphi `$ space. The resulting image of a source, produced by a single SPH particle, thus assumes a trapezoidal-like shape when mapped to a surface, centered at the source. Similarly, The concentric ’rings’ discussed above are in fact rectangles in $`\theta \varphi `$ space. The identified sources are thus depicted in Figure 13 as ellipses instead of rectangles for illustrative purposes only.