# Gamma-Rays from Intergalactic Shocks

###### Abstract

Structure formation in the intergalactic medium (IGM) produces large-scale, collisionless shock waves, where electrons can be accelerated to highly relativistic energies. Such electrons can Compton scatter cosmic microwave background photons up to -ray energies. We study the radiation emitted in this process using a hydrodynamic cosmological simulation of a universe. The resulting radiation, extending beyond TeV energies, has roughly constant energy flux per decade in photon energy, in agreement with the predictions of Loeb & Waxman (2000). Assuming that a fraction of the shock thermal energy is transferred to the population of accelerated relativistic electrons, as inferred from collisionless non-relativistic shocks in the interstellar medium, we find that the energy flux of this radiation, , constitutes of the extragalactic -ray background flux. The associated -ray point-sources are too faint to account for the unidentified EGRET -ray sources, but GLAST should detect and resolve several -ray sources associated with large-scale IGM structures for , and many more sources for larger . The intergalactic origin of the shock-induced radiation can be verified through a cross-correlation with, e.g., the galaxy distribution that traces the same structure. Its shock-origin may be tested by cross-correlating its properties with radio synchrotron radiation, emitted as the same accelerated electrons gyrate in post-shock magnetic fields. We predict that GLAST and Čherenkov telescopes such as MAGIC, VERITAS and HESS should resolve -rays from nearby (redshifts ) rich galaxy clusters, perhaps in the form of a diameter ring-like emission tracing the cluster accretion shock, with luminous peaks at its intersections with galaxy filaments detectable even at .

^{†}

^{†}copyright: ©:

## 1 Introduction

In the past three decades, clear evidence has emerged indicating the
existence of an extragalactic -ray background (EGRB).
The origin of this radiation, however, has remained highly speculative.
The first unambiguous detection of isotropic extragalactic
-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 on
scales) spectrum in the energy range 30 MeV-120 GeV, well fitted by a
single power law, with photon number density per energy interval
Sreekumar98 . 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 , scatter a small fraction of the cosmic microwave background (CMB) photons up to -ray energies, as inverse-Compton radiation. The estimated energy flux from this process per decade in photon energy, , is in good agreement with the detected EGRB, for typical cosmological parameters, provided that the shocked matter accumulated a present-day mass-average temperature , and that a fraction of the associated thermal energy was transferred to the relativistic electrons. In this model, nearby rich galaxy clusters are expected to be bright -ray sources, as they host strong, dense shocks LoebWaxman2000 ; Totani2000 .

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 -ray emission, yet to be detected.

The model presented above has several interesting implications, other than an identification of the origin of the EGRB: -ray emission could be the first tracer of the undetected warm-hot intergalactic medium (WHIM), gas at temperatures 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 () frequencies WaxmanLoeb2000 . 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 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 of the shock thermal energy. We adopt the value , 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,
in the energy range and
in the energy range ,
in good agreement with the prediction of Loeb & Waxman (2000). The
calculated energy flux (per decade in photon energy), , constitutes of the EGRB flux and is smaller
by a factor of than estimated above, mainly due to the lower
present-day temperature reached by the simulation,
,
and its slower heating rate. We predict the existence of -ray
halos, associated with large-scale structure.
Although too faint to account for the EGRET unidentified -ray sources,
several of these halos should be detectable by future experiments, such as the
*Gamma-Ray Large-Area Space Telescope* (GLAST, planned
^{1}^{1}1See http://glast.gsfc.nasa.gov to be launched in 2006),
provided that .
The number of detectable sources is sensitive to the fraction of
shock-energy transferred to the relativistic electrons,
with well-resolved GLAST sources predicted for .
Detection of intergalactic shock-induced -ray sources may be used to
calibrate the free parameter .
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 ^{2}^{2}2See http://hegra1.mppmu.mpg.de/MAGICWeb to
become operational during 2002).
These images reveal an elliptic ring-like feature,
of diameter corresponding to , surrounding a nearby
() 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 -ray sky maps and study the features and distribution of the simulated -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 -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 (), 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:

(1) |

where and are the length and mass scales of structures forming at the present time and is a slowly varying function of the redshift , satisfying , such that smaller scales become non-linear first. These parameters are sensitive to the cosmological model. We use a ’concordance’ model of Ostriker & Steinhardt (1995) - a flat universe with normalized vacuum energy density , matter energy density , baryon energy density , Hubble parameter , and initial perturbation spectrum of slope and normalization - to find:

(2) |

where . 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 , collapsing over a time period , will achieve velocities of order . Since a structure that became non-linear at redshift will collapse over a time period of order , where is the Hubble parameter, the formation of such an object is expected to involve shocks with velocity:

(3) |

where is a dimensionless number of order unity, is the typical shock velocity at and . For the model described above we find:

(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 , and is sufficiently rarefied such that the shocks concerned are collisionless. Using the Rankine-Hugoniot shock adiabat Landau , we can relate the shocked gas temperature to the shock velocity,

(5) |

where is the temperature of gas, shocked at , is the Boltzmann constant, and is the ratio between the post-shock (downstream) sound velocity and the shock velocity (relative to the unshocked gas), approaching for a strong shock. For the model above we find:

(6) |

As we show in §3.2, this estimate of is indeed valid only to order of magnitude and appears to exceed the true value by a factor .

### 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 Drury83 ; Blandford87 . These power law distributions extend up to energies in the Earth’s bow shock, to in SNe remnants Tanimori98 where shock velocities are similar to those of intergalactic shocks, and over in Galactic cosmic rays.

Magnetic fields of amplitude are measured in halos of galaxy clusters Kronberg94 and in the Coma super-cluster Fusco99 , suggesting that the IGM in which the shocks propagate is significantly magnetized. The measured magnetic field energy constitutes a fraction of the typical post-shock thermal energy in cluster environments WaxmanLoeb2000 .

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 (), 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:

(7) |

where 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 Drury83 ; Blandford87 , 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:

(8) |

where is the compression factor of the shock. For a shock in an ideal gas, is given by:

(9) |

where is the Mach number, and and denote the shock velocity and the pre-shock (upstream) sound velocity, respectively. The shock waves of interest are generally strong, i.e. involve high () Mach numbers. Hence, we find , and thus .

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 LoebWaxman2000 by:

(10) |

where is the electron Larmor radius, is the electron Lorentz factor - - in units of , is the magnetic field in units of , and is the downstream temperature of the gas, in . 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:

(11) |

where is the CMB energy density. Synchrotron cooling of the electrons is less efficient, with a characteristic cooling time , where is the magnetic field energy density. Thus, inverse-Compton cooling is faster, as long as the downstream magnetic field satisfies . Equating and leads to the maximal Lorentz factor:

(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 , as observed in cluster halos and which corresponds to a fraction of the post-shock thermal energy, the highest energy electrons up-scatter CMB photons to energy . It should be pointed out here, that the observed sub- 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 Kulsrud97 , by “contamination” by the first generation of stars Rees87 or radio sources Daly90 ; Furlanetto01 , 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 Helfand87 ; Cargill88 and -ray bursts Gruzinov99 ; Medve99 , 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 , which reproduce the observed large scale field strengths.

We can now parameterize the normalization of the relativistic electron distribution, by assuming that a fraction of the thermal energy density induced by the shock, , is transferred to these electrons:

(13) |

where is the total energy density of relativistic electrons and and subscripts indicate upstream (pre-shock) and downstream (post-shock) values, correspondingly. The gas number density, , is related to the baryon number density, , by , where is the average particle mass in an ionized gas with a hydrogen mass fraction . For strong shocks and , thus , and we find:

(14) |

#### 2.2.2 Acceleration Efficiency

The efficiency of particle acceleration in structure formation shocks, parameterized above by , 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 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 . The physics of shock waves is essentially determined by three parameters: the shock velocity , the pre-shock density and the pre-shock magnetic field . 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 , where is the plasma frequency of the ions. The pre-shock magnetic field, parameterized by the cyclotron frequency of the ions, , may not be scaled out of the problem as well. However, by comparing to the growth rate of electromagnetic instabilities in the shocked plasma, , we find that for strong shocks . 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; Tanimori et al. 1998).
Dyer et al. (2001) have modeled the multi-frequency emission from SN1006,
finding a shock efficiency of .
The total energy they find in relativistic electrons constitutes a fraction
out of the *total supernova explosion energy*
(, as often quoted in the literature).
Other authors Mastichiadis96 ; Aharonian99 have estimated
in SN1006, which corresponds to values
in the range of .
Ellison, Slane and Gaensler (2001) have modeled the emission from supernova
remnant , estimating the energy of relativistic electrons
to constitute at least - and more likely - of the shock
kinetic energy flux, corresponding to .

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 Reynolds92 ; Baring99 , and thus carry a large fraction (tens of percents, e.g. 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 between the relativistic electron energy and the energy of relativistic ions. One may easily estimate the value of in interstellar medium cosmic rays: using the electron-to-proton ratio in a given energy, e.g. at GeV energies Barwick98 , correcting for electron cooling and integrating over energies, we find , i.e. . However, since there is no clear relation between in the interstellar medium and its value immediately behind the shocks found in SNR’s and in the intergalactic medium, this estimate of is controversial.

To summarize, a shock efficiency of 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 in the range 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, , the average relative increase in photon energy is given by . Hence, an electron with energy will scatter a CMB photon to an average energy:

(15) |

where we have defined , and 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, , is given by equation (12), and equation (11) imposes a minimal energy threshold, as electrons with Lorentz factors smaller than hardly radiate during a Hubble time . Plugging this range of Lorentz factors into equation (15) shows that the emitted photon spectrum extends from up to 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:

(16) |

For a strong shock, where the electron number density satisfies , this reduces to:

(17) |

For weak shocks, the approximation in equation (16) introduces a small numerical error (), 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 to the observed value of , is roughly proportional to . Thus, the main contribution to the spectrum observed today, especially at high photon energies, is expected to originate from recent shocks at . Assuming that a fraction of the baryons in the universe were recently shocked to a temperature , when already approached , we find from equation (14) and equation (17):

(18) |

The inverse-Compton flux per decade in photon energy thus becomes:

(19) |

This result is in excellent agreement with observational data, when supplemented by the contribution from unresolved point sources LoebWaxman2000 . 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 (, , ), a similar result LoebWaxman2000 :

(20) |

The synchrotron radiation emitted by the accelerated electrons may now be estimated. A relativistic electron with a Lorentz factor , gyrating in a magnetic field of amplitude , emits photons with characteristic energy . Synchrotron radiation is thus expected in photon energies ranging from and up to . The emitted flux can be crudely estimated as:

(21) | |||||

whereas a more careful approach is to assume that the energy density of the post-shock magnetic field constitutes a fraction of the shock-induced thermal energy. This assumption, combined with the Press-Schechter mass function, leads to a similar (somewhat higher) flux:

(22) |

where a typical value produces magnetic fields consistent with observations of the outer halos of X-ray clusters WaxmanLoeb2000 .

#### 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 , heating it up to a temperature by strong shocks. The baryon number accretion rate is thus given by . We assume, as previously, that a fraction of the post-shock thermal energy is converted into relativistic electrons. The maximal electron Lorentz factor is determined from equation (12), and depends on the assumed magnetic field. Electrons with cooling time longer than the virialization time of the object, , 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 , where . This effect reduces the total luminosity by a factor . Combining these considerations, we find that the total and specific luminosities of the object are given by:

(23) |

and:

(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 , accreting over a virialization time and shocked to a temperature . A magnetic field of amplitude , as suggested by observations, gives and . Whence,

(25) |

and:

(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 in the inverse-Compton spectrum and fluctuations 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 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’ model of Ostriker & Steinhardt (1995), presented in §2.1. An initial Harrison-Zel’dovich spectrum of fluctuations was assumed, and linear structure formation theory EisensteinHu employed to predict the resulting spectrum at a certain redshift , 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 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 BarnesHut .

The baryonic component of the universe is modeled by the simulation as an ideal gas, with an adiabatic index , and is represented by a large number of particles, similar to the dark matter. The gas is governed by the continuity equation:

(27) |

the Euler (momentum) equation:

(28) |

and the thermal energy equation (the first law of thermodynamics):

(29) |

where , , and are the pressure, mass density, thermal energy per unit mass and velocity of a gas element; is the gravitational potential, found from Poisson’s equation:

(30) |

and we have used Lagrangian derivatives:

(31) |

An artificial viscosity Steinmetz term was added to the above equations, in order to capture shocks. A cooling function 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 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
, 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
, and in cold (),
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 . At the intersection of the filaments appear denser, spheroidal objects, with typical size , 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, (see Figure 2), is significantly lower than the order-of-magnitude estimate presented earlier. (Note that the approximation of eq. [5] 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 , 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 , 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, , is comparable to or smaller than the dynamical time, . Such regions, where , presently contain approximately 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 . They identified very strong accretion shocks with Mach numbers of up to around the filaments and halos of their simulations. These halos exhibit, in addition, weaker shocks (), 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 (, ) and after (, ) passing
through the shock, directly, or using the entropy-change, , induced by
the shock. Using the Rankine-Hugoniot adiabat, we find that the
compression factor of a shock - related to its Mach number by
equation (9) - can be determined from according
to:

(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 of SPH particles that accumulated
a normalized entropy change between and ,
where 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 , and a large population of stronger shocks,
with Mach numbers up to a few , 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 , such that only particles experiencing a larger entropy-increase between two consecutive snapshots are considered shocked. A threshold value was found to be reasonable, considering numerical noise, although it entails the danger of eliminating some of the mildly weak shocks (), if an SPH particle requires many snapshots to cross such a shock. The fact that only of the gas suffers entropy-changes larger than simplifies the computational work considerably. Demanding that the entropy increases at least by during consecutive snapshots, and dismissing shocked particles that did not reach a temperature 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 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, , 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 [2000] and Martel & Shapiro [2001].) 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: . Hence, we take:

(33) |

as the post-shock temperature, where and are related by equation (32). This choice yields more realistic results than the alternative, , as discussed in §4.1.2.

#### 4.1.2 Shock-Extraction Results

The shock-extraction scheme, described above, identifies shocks in of the simulated gas mass during recent epochs (), where about 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 were located. It confirms, however, our suspicion of
loosing most of the ’weak’ shock population (). 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 ().

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 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 . Miniati et al. (2000) find, for a slightly different model (, ), that the energy passing through shocks with Mach numbers between redshift 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 ).
Estimating this energy gain based on the temperature *after* passing
through the shock (i.e. taking ) results in shock energies lower,
on average, by a factor .
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 -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 . Second, the contribution of weak shocks to the -ray emission, especially at high photon energies, is small due to the steep distribution of the electrons accelerated by such shocks, reducing their -ray efficiency. At , for example, shocks with have efficiency (compared to strong shocks), dropping below for . Weak shocks could, nonetheless, increase the -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 of the post-shock density of the thermal energy induced by a shock, , is converted into a power law distribution of relativistic electrons, . The power law index is related to the compression factor by , where was used here and is assumed henceforth. Imposing a distribution cutoff at the maximal energy attainable by the electrons, , fixes the normalization:

(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:

(35) |

The average mass per particle is given by

(36) |

where is the hydrogen mass fraction, is the ionization ratio (per nucleon) and the second equality holds for a fully ionized gas, where , with the hydrogen mass fraction of the simulation, , assumed henceforth.

We parameterize the post-shock magnetic field, , by assuming that a fraction of the post-shock thermal energy is transferred into the energy of this field, giving:

(37) |

where is the average particle number density at redshift . Hence,

(38) |

which for strong shocks, where and , reduces to:

(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:

(40) |

where we have defined . This leads to a partial differential equation for the electron distribution:

(41) |

where the notation is introduced. We ignore adiabatic cooling due to the Hubble expansion since the cooling timescales relevant for our calculated –ray spectrum are all much shorter than the Hubble time. For the initial (power-law) condition we find the solution:

(42) |

where we have defined . This solution implies that the electron distribution steepens (if ) towards , beyond which all electrons have already cooled off. For strong shocks, where , the electron distribution is stationary for , 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 , approximately given by . This implies that electrons with , where is calculated in equation (A7), do not have sufficient time to cool between redshift 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 :

(43) |

where is the specific number density of the incident photon field, assumed constant (in the time period of interest) and isotropic, is the emission function, given by:

(44) |

and we have assumed Thomson scattering in the electron rest frame: . The incident CMB photon field is approximately blackbody radiation, whence:

(45) |

The assumptions made above are all satisfied for the ultra-relativistic electrons () 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 , 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 , 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 loose all
their energy to radiation, *immediately* after passing through the
shock front. The error thus introduced is small, for electrons satisfying
. 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 . Thus, it is obviously pointless to follow the
temporal evolution of electrons, satisfying .

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

(46) |

Since the only time-dependence in the emitted power of equation (43) appears in the distribution of electrons, , we may change the order of integration and first calculate the time-integrated electron distribution:

(47) |

where the second equality is obtained by taking
the upper integration limit to infinity or to , 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 . This leads to:

(48) | |||||

where is a numerical coefficient, defined by:

(49) |

is the Riemann zeta function, and we defined . This approximate result holds for emitted photon energies in the range . Writing this as:

(50) |

indicates that the approximation holds for the entire relevant -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 . The two would become identical, if one was to replace in the latter by , and take . For a strong shock, we indeed find that , and precisely recover equation (17) - the early estimate for a strong shock. This is not surprising, recalling that the electron distribution for is mostly stationary, rendering evolutionary effects insignificant. The numerical mismatch between equation (16) and equation (48) varies as function of , but remains of order for relevant values of ().

#### 4.3.2 Space-Time Integration

Two different space-time integration schemes were used.
The first, *line of sight* integration, calculates the specific
intensity at a given direction on the sky,
where is its angle with the -axis and is the angle of
its projection on the - plane with the -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
,
the correlation of intensities at various angular separations,
etc. The second scheme provides *direction bin* integration,
calculating the flux predicted within a small region
on the sky, defined by and .
The flux, related to the specific intensity by
,
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 and redshift , lie at coordinates of a comoving coordinate system centered upon the observer. The contribution of this shock to the specific flux , detected by the observer today, satisfies:

(51) |

where is the luminosity distance between source and observer, is the present-day value of the cosmological scale factor, and we have defined . With the spatial and temporal resolution used, many () 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 to the detected flux is determined by:

(52) |

The specific intensity contributed by the shock is calculated by approximating the latter as a sphere of uniform specific brightness . The proper radius of this sphere is determined by the mass density of the SPH particle, according to , where is the mass of an SPH particle. The brightness of the sphere is related to its detected flux by: