Semi-analytic modelling of the extragalactic background light and consequences for extragalactic gamma-ray spectra

Mon. Not. R. Astron. Soc. 422, 3189–3207 (2012) doi:10.1111/j.1365-2966.2012.20841.x Semi-analytic modelling of the extragalactic background light and consequences for extragalactic gamma-ray spectra Rudy C. Gilmore,1,2Rachel S. Somerville,3,4Joel R. Primack2 and Alberto Dom´ ınguez5,6,7,8† 1SISSA, via Bonomea 265, 34136 Trieste, Italy 2University of California, Santa Cruz, CA 95064, USA 3Space Telescope Science Institute, Baltimore, MD 21218, USA 4Department of Physics and Astronomy, The Johns Hopkins University, Baltimore, MD 21218, USA 5Santa Cruz Institute for Particle Physics (SCIPP), University of California, Santa Cruz, CA 95064, USA 6Instituto de Astrof´ ısica de Andaluc´ ıa, CSIC, Apdo. Correos 3004, E-18080 Granada, Spain 7Departamento de F´ ısica At´ omica, Molecular y Nuclear, Universidad de Sevilla, Apdo. Correos 1065, E-41080 Sevilla, Spain 8Department of Physics and Astronomy, University of California, Riverside, CA 92521, USA Accepted 2012 February 28. Received 2012 January 31; in original form 2011 April 4 ABSTRACT Attenuation of high-energy gamma-rays by pair production with ultraviolet, optical and infrared (IR) extragalactic background light (EBL) photons provides a link between the history of galaxy formation and high-energy astrophysics. We present results from our latest semianalytic models (SAMs), which employ the main ingredients thought to be important to galaxy formation and evolution, as well as an improved model for reprocessing of starlight by dust to midand far-IR wavelengths. These SAMs are based upon a cold dark matter hierarchical structural formation scenario, and are successful in reproducing a large variety of observational constraints such as number counts, luminosity and mass functions and colour bimodality. Our fiducial model is based upon a Wilkinson Microwave Anisotropy Probe 5-year cosmology, and treats dust emission using empirical templates. This model predicts a background flux considerably lower than optical and near-IR measurements that rely on subtraction of zodiacal and galactic foregrounds, and near the lower bounds set by number counts of resolvable sources at a large number of wavelengths. We also show the results of varying cosmological parameters and dust attenuation model used in our SAM. For each EBL prediction, we show how the optical depth due to electron–positron pair production is affected by redshift and gamma-ray energy, and the effect of gamma-ray absorption on the spectra of a variety of extragalactic sources. We conclude with a discussion of the implications of our work, comparisons to other models and key measurements of the EBL and a discussion of how the burgeoning science of gamma-ray astronomy will continue to help constrain cosmology. The low EBL flux predicted by our fiducial model suggests an optimistic future for further studies of distant gamma-ray sources. Key words: cosmology: theory – diffuse radiation – gamma-rays: general – infrared: diffuse background. 1 INTRODUCTION The extragalactic background light (EBL) is the totality of light emitted by stars and active galactic nuclei (AGN) over the lifetime of the Universe. Today, this pervasive photon background consists of light emitted at all epochs, modified by redshifting and dilution E-mail: [email protected] †Visiting researcher at the Santa Cruz Institute for Particle Physics (SCIPP). due to the expansion of the universe. The bulk of the EBL occurs at wavelengths from the near-ultraviolet (UV) to the far-infrared (IR). In the UV, optical and near-IR most of the EBL is due to direct starlight,aswellasasubdominantcontribution from AGN (Schirber & Bullock 2003). From the mid-IR to submillimetre wavelengths, the EBL consists of reemitted light from dust particles, including both continuum thermal radiation and line emission from polycyclic aromatic hydrocarbons (PAH) molecules (Lagache, Puget & Dole 2005). The background at longer wavelengths is dominated by the cosmic microwave background, while shortward of the Lyman limit C 2012 The Authors Monthly Notices of the Royal Astronomical Society C 2012 RAS 3190 R. C. Gilmore et al. the background flux decreases rapidly due to attenuation by neutral hydrogen in stellar atmospheres and the interstellar and intergalactic media. Because the production of the EBL is directly linked to the star formation history of the universe, limits on the EBL can be used to provideconstraintson the historyofgalaxyformationand evolution. Observations of the extragalactic sky brightness can constrain the local background, however, they do not provide information about evolution of the background with redshift. Direct sky photometry has been attempted with a number of instruments, most notably DiffuseInfraredBackground Experiment (DIRBE)andFar-Infrared Absolute Spectrophotometer (FIRAS), and also ground-based optical telescopes (e.g. Mattila et al. 2011), but these type of measurements are subject to considerable uncertainties due to the large foreground sources that must be subtracted (Hauser & Dwek 2001). Integration over discrete sources seen in galaxy surveys (e.g. Madau & Pozzetti 2000; Keenan et al. 2010) is another way to estimate the EBL, but one that in principle can only provide a lower limit due to the possibility of unseen sources beyond the magnitude limits of the survey instrument or underestimation of true total luminosity of galaxies due to light in the faint outskirts (e.g. Bernstein 2007). Observations with highly sensitive satellite instruments have provided us with EBL lower limits from galaxy number counts across wide wavelength ranges. High-energy gamma-rays can interact with EBL photons in electron–positron pair production interactions (Nikishov 1962; Jelley 1966; Gould & Schreder 1967). By effectively removing these gamma-rays from view, this process has the potential to alter the observed spectra of extragalactic high-energy sources, and increasingly occlude those at higher redshifts. The rapid development of ground-based gamma-ray astronomy in the past 20yr has led to a number of attempts – e.g. Dwek & Krennrich (2005), Aharonian et al. (2006), Mazin & Raue (2007), Albert et al. (2008) – to constrain the EBL based on modification to gamma-ray spectra, a method that can provide a measurement of the EBL that is independent of direct observation. In principle, the cosmological history of the EBL could be reconstructed by comparing observations of high-energy sources at different redshifts to their intrinsic spectra. Unfortunately, the emission mechanisms and intrinsic spectra of GeV and TeV sources are still poorly understood. Understanding how the EBL is produced and how its spectral energy distribution (SED) evolves in redshift requires an understanding of the sources responsible for its production. This has been attempted by different authors using a variety of techniques. As enumerated in Dom´ ınguez et al. (2011, hereafter D11), calculations of the EBL fall into four general categories: (i) forward evolution beginning with initial cosmological conditions, such as the semi-analytic models (SAMs) used in this work; (ii) backwards evolution of the well-constrained present-day galaxy emissivity according to some prescription; (iii) evolution of galaxy properties that are inferred over some range in wavelength; (iv) direct observation of evolution in galaxy properties over the redshifts providing the major contribution to the background light, a category which describes the empirical method developed in D11. The last two of these have become much more powerful techniques in recent years due to large-scale surveys by groundand space-based instruments, especially at UV and IR wavelengths, where a great deal of progress has taken place in the last decade. Some of the first models to account for EBL production by the evolving galaxy population were Madau, Pozzetti & Dickinson (1998) and Franceschini (2001), using Hubble Space Telescope (HST)andInfrared Space Observatory (ISO) data, respectively, and Pei, Fall & Hauser (1999), who looked at chemical enrichment data in Lyαsystems. A two-part paper series by Kneiske and collaborators (Kneiske, Mannheim & Hartmann 2002; Kneiske et al. 2004) computed the EBL and subsequently predictions for attenuation of gamma-ray sources based on a parametrization of the star formation rate (SFR) density. These models separately include the contribution of the luminous infrared galaxy (LIRG)/ultraluminous infrared galaxy (ULIRG) population. An update to this work by Kneiske & Dole (2010) attempt to create a model using a similar method that produces a minimal background, and a similar method was employed in Razzaque, Dermer & Finke (2009) and Finke, Razzaque & Dermer (2010). The aforementioned work of D11 used observedevolutionuptoz=4intheK-band luminosity function combined with the evolving distribution of 25 galaxy SED types from a multiwavelength survey of galaxies to estimate the EBL and its evolution. Other authors have used backward evolution models to predict the EBL. These calculations begin with the present day galaxy luminosity function and attempt to trace this function backwards in time by assuming a functional form for the redshift evolution. In Malkan & Stecker (1998, 2001), IR luminosity functions from IRAS were extrapolated backwards in redshift using power-law functions. The model of Stecker, Malkan & Scully (2006) updated this work and computed the EBL below the Lyman limit (13.6eV) for two different cases of stellar evolution. The model of Rowan-Robinson (2001) also utilized a 60μm evolving luminosity function, and a four-component spectral model for IR and optical emission. One potential problem with this method is that it has difficulty accounting for the emissivity contribution of merger-triggered starbursts, believed to contribute an increasing fraction of the SFR density and IR emissivity with increasing redshift. Franceschini et al. (2001) made an attempt to account for this starburst phase in a backwards evolution model. Franceschini, Rodighiero & Vaccari (2008) published a sophisticated model using observed luminosity functions, and used it to calculate the EBL and gamma-ray attenuation. This model uses evolving luminosity functions in the near-IR up to z= 1.4 for two different galaxy populations (spiral and spheroidal) and local luminosity functions for the irregular/starbursting population, combined with synthetic SEDs to find the total emissivity. In forward evolution scenarios such as SAMs, predictions for the evolution of galaxy emissivities are made by beginning from the universe in its primordial state and simulating the process of galaxy formation. This is considerably more involved and challenging than the other methods of estimating the EBL, but can provide a degree of insight into the fundamental astrophysics processes that determine the emissivity that is lacking in other approaches. SAMs of structure formation based on cold dark matter (CDM) merger trees have been usedin several papersbyour grouptopredict theEBL.Primack etal. (1999) predicted the EBL using the SAM described in Somerville & Primack (1999) and Somerville, Primack & Faber (2001). Later work included improved treatment of absorption and reemission of starlight by dust and updated cosmological data (Primack, Bullock & Somerville 2005; Primack, Gilmore & Somerville 2008). In a companion paper to this work, Somerville et al. (2012, hereafter SGPD12), we present a new SAM based on galaxy formation in a Wilkinson Microwave Anisotropy Probe 5-year (WMAP5) cosmology. This model, which will be summarized in the next section, incorporates the physical processes thought to be most important in determining the evolution of these systems. From the luminosity density calculated in this model, we have predicted the evolving EBL out to high redshift. In this paper, we address the topic of gamma-ray attenuation and show how our estimated EBL affects C 2012 The Authors, MNRAS 422, 3189–3207 Monthly Notices of the Royal Astronomical Society C 2012 RAS Modelling of the EBL and gamma-ray spectra 3191 high-energy observations of extragalactic sources, and discuss the ability of gamma-ray telescopes to explore the distant universe. We also present a comprehensive comparison to the predictions for EBL and gamma-ray opacity that have been proposed by a number of recent authors. Our current work is an update to Primack et al. (2008),which presentsanearlier stageof ourresultsusing a‘concordance cosmology’ (CCDM), with parameters largely consistent with WMAP1. This work is also closely related and complementary to Gilmore et al. (2009), which used the CCDM model as the basis for a prediction of the UV background radiation out to high redshift, and therefore emphasized the calculation of optical depths for gamma-rays below 200 GeV. In that paper we included contributions to the UV emissivity from quasars, as well as an account of the attenuation of ionizing radiation escaping from galaxies and processing by neutral hydrogen in the intergalactic medium (IGM) using a radiative transfer code. In the following section, we briefly review the key elements of the SAM presented in SGPD12. Results are presented in Section 3, beginning with a review of key results related to the evolving backgroundradiationfromSGPD12 inSubsection 3.1.In Subsection3.2, we show how the population of photons in our calculated EBL impacts observations of extragalactic gamma-ray sources through pair production interactions. Subsection 3.3 deals with the comparison between our predicted gamma-ray optical depths and observations of very high energy (VHE) blazars, and constraints on the EBL that other authors have derived using high-energy data. In Section 4, we compare our EBL model with several others that have been proposed using a variety of techniques in recent years. We conclude in Section 5 with a summary of results and a discussion of how current and future high-energy observations will continue to constrain the EBL. 2MODEL This section summarizes the SAM that is used to predict the EBL in this work. This model is based upon the models that were first presented in Somerville & Primack (1999) and Somerville et al. (2001), with significant new ingredients as described in Somerville et al. (2008, hereafter S08). Readers should refer to SGPD12, as well as S08, for details. 2.1 Galaxy formation We assume a standard CDM universe and a Chabrier (Chabrier 2003) stellar initial mass function (IMF) that does not evolve in redshift. The model presented in this work uses cosmological parameters based on WMAP5, including a power spectrum normalization of σ8=0.82, a value that is intermediate between the previous findings of WMAP1andWMAP3. This value is within 1σof the recently published value from WMAP7 of 0.809 ±0.024 (Komatsu et al. 2011). The SAMs used here are based on Monte Carlo realizations of dark matter halo merger histories calculated using the modified Press–Schechter (Sheth & Tormen 1999) and extended Press–Schechter methods. Star formation occurs when gas is accreted by the galaxy and becomes available after cooling via atomic processes. Feedback from supernovae can heat and eject the cold gas within the galaxy. This gas will either be deposited in the hot reservoir connected with the dark matter halo, or returned to the IGM, depending on the wind velocity relative to the virial velocity of the halo. Star formation in our model occurs in two regimes, quiescent star formation in isolated galaxies and merger-driven starbursts. The former is treated using a recipe based on the empirical Schmidt–Kennicutt relation (Kennicutt 1989; Kennicutt et al. 1998). Mergers drive gas deep into galactic nuclei, fuelling black hole growth which power AGNdriven winds. Supermassive black holes can also produce radio jets that heat the hot halo gas and may eventually shut off cooling and eventually lead to a cessation of star formation. This ‘quenching’ of star formation tends to occur in massive galaxies, which are able to build massive black hole, and which are accreting gas from hot, tenuous haloes rather than via cold dense filaments. The chemical enrichment and star formation history of each galaxy are used to predict the total emission spectrum. We have adopted the Bruzual & Charlot (2003) stellar population models in this work. 2.2 Dust extinction and reemission Light emitted by stars can be absorbed and reemitted by dust. In the SAM, dust is modelled as a two-component distribution, using a modified version of the prescription of Charlot & Fall (2000), which treats separately the dense dust in giant molecular clouds that contain newborn stars and the much more diffuse cirrus in the interstellar medium (ISM). Extinction from the ISM component is proportional to the density and metallicity of cold dust, and a slab geometry is assumed. Stars younger than 107yr are enshrouded in a cloud of dust with optical depth τBC,V=μBC τV,0,whereμBC = 3andτV,0is the face-on, V-band extinction of the ISM component. To calculate extinction at other wavelengths, we have assumed a starburst attenuation curve (Calzetti et al. 2000) for the diffuse dustcomponent andapower-lawextinctioncurveAλ∝(λ/5500 Å)n, with n=0.7, for the birth clouds (Charlot & Fall 2000). We consider two possible normalizations for the extinction recipes. In our ‘WMAP5+fixed’ model, parameters are constant at all redshifts, and are adjusted to match observed relations between UV and IR luminosity for nearby galaxies. However, as discussed in SGPD12, we have found that this model has difficulty reproducing luminosity functions at higher redshift in the UV and optical bands. Motivated by this finding, we have created an ‘evolving’ model with redshift-dependent parameters tuned to match the observed UV and optical luminosity functions at all redshifts where they have been measured. In this model, total dust extinction is scaled by a factor (1 +z)−1at non-zero redshifts, and the opacity and lifetime of molecular clouds is scaled by a factor z−1above redshift 1. Because this evolving dust model is found to be more successful at matching high-redshift data, we favour this model and will refer to it as the ‘fiducial’ variant in the next section. The reemission of IR light by the dust due to thermal and PAH emission is estimated in our model using templates that describe the spectra of galaxies from the mid-IR to submillimetre as a function of the total IR luminosity, and are based on observations of galaxies in the local universe. Energy absorbed by dust from direct starlight is redistributed in the infrared according to a prescribed SED. These templates are embedded in our SAM, and account for emission at wavelengths from a few microns to the submillimetre, including the emission and absorption lines appearing in the PAH region. The dust emission templates we have used here are described in Rieke et al. (2009, hereafter R09), and are based on observations of 11 local LIRGs and ULIRGs combined with lower luminosity local systems. A comparison of these templates with those of Devriendt, Guiderdoni & Sadat (1999), Devriendt & Guiderdoni (2000), which were used in Primack et al. (2001, 2005, 2008), is available in SGPD12. Being observationally based, these templates suffer to some extent from starlight contamination at short wavelengths, and C 2012 The Authors, MNRAS 422, 3189–3207 Monthly Notices of the Royal Astronomical Society C 2012 RAS 3192 R. C. Gilmore et al. Figure 1. The SFR density in each of the two SAMs over cosmic time. The solid and dotted black curves are the predictions from our WMAP5and CCDM models, respectively. The dash–dotted red line is the SFR density inferred from the model of D11, which has been converted to a Chabrier IMF. Grey data points are from a compilation presented in S08; and the dashed line is the estimate of Hopkins & Beacom (2006). there is also a discontinuity in the galactic SEDs where the templates are joined with our stellar synthesis models, at about 4μm. In the wavelength range of 2–5 μm we have attempted to compensate for theseproblemsby fitting a power-lawextrapolationto the templates. The overall effect of this change to the integrated background light that we will present in the next section is minimal. 3 RESULTS 3.1 Overview of astrophysical results In this section, we present results for the evolving background light predicted by our SAMs. In addition to our fiducial model, which uses a WMAP5 cosmology combined with the evolving dust model and R09 templates described in the last section, we will also show results using employing the ‘fixed’ dust absorption model to facilitate comparison with our older work. Results from our CCDM model will also be shown when relevant. As in the previous section, much of the material here is a review of key results from SGPD12, and readers should refer to this work for further detail. We will also compare with the work of D11, a more observationally driven model that provides a useful contrast to our theoretical approach. Predictions for the evolving EBL and gamma-ray opacity in the fiducial and fixed models are available online in tabulated form.1 3.1.1 Star formation and luminosity density history The global SFR density arising in our WMAP5andCCDM models is shown in Fig. 1, along with observational estimates obtained using a variety of tracers. Our model makes predictions that are in agreement with the bulk of data for z<1, and tend to be slightly lower than observed at z∼2. All measured SFRs are subject to significant uncertainties, as seen in the scatter in results for the plotted data. Uncertainties in dust extinction impact all results relying on UV luminosity. Measurements of Hαand other spectral lines must take into account extinction as well as metallicity effects. Other 1http://physics.ucsc.edu/~joel/EBLdata-Gilmore2012 authors have attempted to measure SFRs based on 24μm and other mid-IR observations of warm dust. These results can be affected by AGN contamination, as well as PAH features that move in and out of the instrument bandpass with changing redshift. All of these uncertainties grow with increasing redshift, where our knowledge of dust distribution and galaxy SEDs becomes less reliable. The SFR density inferred from the UV and IR luminosity densities of D11 is also shown here. The larger value predicted in this work relative to ours is due to the considerably higher far-IR emission in D11; predictions in the optical and near-IR are similar. Above redshift 1, SFR density predictions from D11 are affected by the assumption made about the evolution of different galaxy spectral types. The luminosity density in our models is predicted by summing over the emission from individual galaxies. Results for total galactic emissivity as a function of wavelength and redshift are shown in the right-hand panel of Fig. 2. In the left-hand panel we show predicted luminosity density in the local universe, compared to constraints at a number of wavelengths. The local luminosity density has been extremely well measured in the optical and near-IR by large-scale surveys such as Sloan Digital Sky Survey (SDSS) and 2-Micron All Sky Survey (2MASS), and this provides a strong constraint on any model of the galaxy population. At longer wavelengths, we show how our models fit the local data at IR wavelengths, including data from IRAS and Submillimetre Common-User Bolometer Array (SCUBA). We have also compared our model with observational estimates of the evolving luminosity density at a number of different redshifts (Fig. 3). The peak emissivity redshift in our model changes depending on the wavelength considered. At UV bands, the emission closely follows the SFR, which peaks at z≈2.5 in our WMAP5 modelandz≈3intheCCDM model. Longer wavelengthsinclude significant contributions from progressively more evolved stellar populations, and therefore peak at later times. Recent evolutionary surveys such as DEEP2 and COMBO-17 allow us to compare the evolution of galaxy emissivity against accurate luminosity density data in several bands. Emissivity in the UV (1500 and 2800 Å) has been seen to increase out to nearly z=2 (Dahlen et al. 2007). Note the large discrepancy between the evolving and fixed dust attenuation models at high redshift in the Bband and UV; in the evolving model much less attenuation of starlight occurs in early star-forming galaxies. At high redshift our fiducial UV predictions are somewhat higher than the measurements of Bouwens et al. (2007); this is largely due to a substantial contribution to the total emissivity from faint galaxies which have very little extinction in our evolving dust model. In the Bband, Dahlen et al. (2005) find that emission increases out to at least z=1; this paper makes the claim that emissivity in the Band Rbands is consistent with being flat in the interval 1 < z<2. Results at the higher redshifts could be sensitive to the faint end slope assumed in calculating the luminosity density. In the K band we match well the local luminosity measurement of Kochanek et al. (2001). Available data at higher redshifts seem to suggest a falloff in emissivity beginning at about z1 which we do not find in our models. As discussed in SGPD12, our model does seem to overpredict the K-band luminosity of galaxies at and below L∗,beginningatredshifts 1–2. A corresponding overproduction of near-IR flux is not seen in the local luminosity density or K-band counts. 3.1.2 Local EBL flux and discrete sources As mentioned in the Introduction, measurements of the local (z= 0) EBL generally fall into two categories: direct sky photometry C 2012 The Authors, MNRAS 422, 3189–3207 Monthly Notices of the Royal Astronomical Society C 2012 RAS Modelling of the EBL and gamma-ray spectra 3193 Figure 2. Left: the luminosity density of the local universe. The solid black line is the WMAP5 model, and the dotted line is the CCDM model. Data at a number of wavelengths are shown from GALEX (blue circles), SDSS (red stars; Montero-Dorta & Prada 2009), 6dF (light blue squares; Jones et al. 2006), 2MASS (green stars; Cole et al. 2001; Bell et al. 2003). In the midand far-IR, the orange squares are from IRAS (Soifer & Neugebauer 1991), while blue stars are from an analysis of local emissivity using data from IRAS,ISO and SCUBA (Takeuchi et al. 2001). Right: three-dimensional representation of the evolution of the luminosity density in our WMAP5 model as a function of wavelength and redshift. and integrated counts of galaxies. Direct measurements provide an absolute measurement of the background light without regard to the sources responsible, but require subtraction of foreground sources present in the Milky Way and our Solar system in order to isolate the extragalactic signal. Integration of galaxy counts (galaxies per unit sky area at a given magnitude) is a way to set firm lower limits on the EBL, although the degree to which these measurements converge on the true value often remains controversial. The flux from faint sources will converge mathematically if the slope of the counts plotted on a log number versus flux diagram is flatter than unity, or in terms of magnitudes if α<0.4, for ln (N)∝αm. As expounded by Bernstein (2007), photometry of faint galaxies is fraught with difficulty in untangling the faint galactic fringes from the background, and it is possible to miss 50 per cent or more of the light associated with extended sources in simple aperture photometry. Large-scale surveys such as the SDSS, the 6-degree Field survey (6dF) and the 2MASS have provided us with an accurate accounting of the galaxies in the local universe, and surveys with the HST have complemented this data with extremely deep counts. Satellite instruments such ISOCAM, IRAC and MIPS provide data in the midand far-IR. A detailed presentation of galaxy number counts in our models compared with data can be found in SGPD12. Our prediction for the local EBL is generally in agreement with lower limits from integrated number counts. In the UV, limits from Gardner, Brown & Ferguson (2000) are considerably higher than the measurement from GALEX (Xu et al. 2005). This may be explained by the former’s use of data from the balloon-based FOCA experiment to find bright counts, which were in disagreement with those from GALEX at several magnitudes. Preliminary Herschel counts data from Berta et al. (2010) set only a weak lower limit on the far-IR background peak, and the author acknowledges that only about half the total IR background is likely being resolved. Absolute measurements of the EBL require the removal of foregroundsources, including stars,ISMemission and sunlight reflected from dust in the inner Solar system (often called ‘zodiacal’ light). The most robust direct measurements of the IR background to date come from DIRBE and FIRAS instruments on the Cosmic Background Explorer (COBE) satellite, though they are still fraught with uncertainty in sky subtraction (see fig. 2 in Hauser & Dwek 2001). The near-IR flux has been calculated from DIRBE observations by a variety of authors (Gorjian, Wright & Chary 2000; Wright & Reese 2000; Cambr´ esy et al. 2001; Wright 2001; Levenson, Wright & Johnson 2007) using foreground source subtraction techniques and modelling of the zodiacal light, and has generally yielded high estimates in this range compared to number counts. Another notable attempt to measure the near-IR background was Levenson & Wright (2008), which used IRAC data to calculate the best-fitting flux at 3.6 μm using a profile-fit to estimate the light from the unobservable faint fringes of galaxies. These results were 70per cent higher than those of the aperture method of Fazio et al. (2004), highlighting the large uncertainties that galaxy fringe issues can bring to EBL measurement. The present-day EBL obtained in each of our models is shown in Fig. 4. We also show results from D11 for comparison. The local EBL is calculated by integrating over the luminosity density at all wavelengths beginning at z=7.5, and accounting for the redshifting and dilution of photons as the universe expands. The EBL at a redshift z0and frequency ν0in proper coordinates can be written as (Peebles 1993) J(ν0,z 0)=1 4π∞ z0 dl dz (1 +z0)3 (1 +z)3(ν,z)dz, (1) where (ν,z) is the galaxy emissivity at redshift zand frequency ν=ν0(1 +z)/(1 +z0), and dl/dzis the cosmological line element, which is dl dz=c (1 +z)H0 1 m(1 +z)3+ (2) for a flat CDM universe. We assume here that the EBL photons evolve passively after leaving their source galaxies and are C 2012 The Authors, MNRAS 422, 3189–3207 Monthly Notices of the Royal Astronomical Society C 2012 RAS 3194 R. C. Gilmore et al. Figure 3. The luminosity density (integrated luminosity function of sources within a given redshift range) in our models versus redshift at 1500, 2800 Å, and in the B,zand Kbands (approximately 4500, 9130 Å and 2.2 µm, respectively). The final panel on the lower right shows the total amount of energy that is absorbed and reradiated by dust at IR energies; units on the y-axis for this panel are solar luminosities per Mpc3. The solid black line is our WMAP5 prediction with evolving dust, and the dotted line are prior results from our CCDM model. Dash–dotted violet shows the predictions from our WMAP5 model using fixed dust attenuation parameters. The long–short dashed red line is the prediction of D11. Observational data shown here are as follows. 1500Å: blue squares are from Dahlen et al. (2007), red stars are from Schiminovich et al. (2005), green stars are from Bouwens et al. (2007) and orange circles are from Reddy et al. (2008). The solid purple circle is a local measurement with GALEX by Wyder et al. (2005). 2800Å: blue squares and the purple circle are again from Dahlen et al. (2007) and Wyder et al. (2005), respectively. Red stars are from Gabasch et al. (2006). Bband: blue squares are from Dahlen et al. (2005). DEEP and COMBO-17 data from Faber et al. (2007) are shown as red stars and open red squares, respectively (these are very similar). This Combo-17 data is an update to that originally presented in Wolf et al. (2003), and we show the original points as green stars. The work of Marchesini et al. (2007) is shown as open purple hexes. zband: local measurements are provided by Montero-Dorta & Prada (2009) (red) and Blanton et al. (2003) (green). Blue squares are from Gabasch et al. (2006). Kband: the local determination is from Kochanek et al. (2001). High-redshift data are from Barro et al. (2009) (blue squares) and Arnouts et al. (2007) (open red circles). Total IR luminosity: observational estimates of the IR emissivity are from Caputi et al. (2007) (open blue pentagons), Reddy et al. (2008) (green circles), Rodighiero et al. (2010) (purple stars) and Le Floc’h et al. (2005) (red crosses). not affected by any further interactions except for cosmological redshifting. This is an acceptable approximation for photons at energies below the Rydberg energy of 13.61 eV. Photons above this energy are strongly attenuated by neutral hydrogen when leaving their galaxy of origin. At higher energies, photons are also capable of interacting with residual neutral hydrogen and, if sufficiently energetic, neutral and singly ionized helium in the IGM. The effect of these processes on the ionizing EBL is the topic of our previous work in Gilmore et al. (2009). The total flux of the integrated EBL as well as the contributions from the optical–near-IR and far-IR peaks and the mid-IR valley for each model are shown in Table 1. 3.1.3 Evolution of the background flux A correct determination of gamma-ray opacity at distances beyond the very nearby universe, z>0.05, requires accounting for the redshift-dependent evolution of the background at all wavelengths. The sharply increasing SFR density from z=0 back to z∼2, combinedwith the(z+1)4evolution ofproperflux inredshift means that the background was considerably more powerful in the recent past, a fact that can only be neglected in gamma-ray attenuation calculations for the closest extragalactic sources. With observations of VHE extragalactic sources now stretching out to redshifts of over C 2012 The Authors, MNRAS 422, 3189–3207 Monthly Notices of the Royal Astronomical Society C 2012 RAS Modelling of the EBL and gamma-ray spectra 3195 Figure 4. The predicted z=0 EBL spectrum from our fiducial WMAP5 model (solid black) and WMAP5+fixed (dash–dotted violet) dust parameters, and CCDM (dotted black) models, compared with experimental constraints at a number of wavelengths. D11 is shown for comparison in dashed–dotted red with the shaded area indicating the uncertainty region. Data: upward pointing arrows indicate lower bounds from number counts; other symbols are results from direct detection experiments. Note that some points have been shifted slightly in wavelength for clarity. Lower limits: the blue–violet triangles are results from HST and Space Telescope Imaging Spectrograph (STIS; Gardner et al. 2000), while the purple open triangles are from GALEX (Xu et al. 2005). The solid green and red triangles are from the Hubble Deep Field (Madau & Pozzetti 2000) and Ultra Deep Field (Dolch & Ferguson, in preparation), respectively, combined with ground-based data, and the solid purple triangle is from a measurement by the Large Binocular Camera (Grazian et al. 2009). In the near-IR J,Hand K bands, open violet points are the limits from Keenan et al. (2010). Open red triangles are from IRAC on Spitzer (Fazio et al. 2004), and the purple triangle at 15µm is from ISOCAM (Hopwood et al. 2010) on ISO. The lower limits from MIPS at 24, 70 and 160 µmonSpitzer are provided by B´ ethermin et al. (2010) (solid blue) and by Chary et al. (2004), Frayer et al. (2006) and Dole et al. (2006) (solid gold, open gold and open green, respectively). Lower limits from Herschel number counts (Berta et al. 2010) are shown as solid red triangles. In the submillimetre, limits are presented from the BLAST experiment (green points; Devlin et al. 2009). Direct detection: in the optical, orange hexagons are based on data from the Pioneer 10/11 Imaging Photopolarimeter (Matsuoka et al. 2011), which are consistent with the older determination of Toller (1983). The blue star is a determination from Mattila et al. (2011), and the triangle at 520nm is an upper limit from the same. The points at 1.25, 2.2 and 3.5 µm are based upon DIRBE data with foreground subtraction: Wright (2001, dark red squares), Cambr´ esy et al. (2001, orange crosses), Levenson & Wright (2008, red diamond), Gorjian et al. (2000, purple open hexes), Wright & Reese (2000, green square) and Levenson et al. (2007, red asterisks). In the far-IR, direct detection measurements are shown from DIRBE (Schlegel, Finkbeiner & Davis 1998; Wright 2004, solid red circles and blue stars) and FIRAS (Fixsen et al. 1998, purple bars). Blue–violet open squares are from IR background measurements with the AKARI satellite (Matsuura et al. 2011). Table 1. The integrated flux of the local EBL in our models (WMAP5 with evolving and fixed dust parameters, and the CCDM model) and the model of D11. Units are nW m−2sr−1. Wavelength range WMAP5 (fiducial) WMAP5+fixed CCDM D11 Optical–near-IR peak (0.1–8µm) 29.01 24.34 26.15 24.47 Mid-IR (8–50µm) 4.89 5.16 5.86 5.24 Far-IR peak (50–500 µm) 21.01 22.94 24.08 39.48 Total (0.1–500 µm) 54.91 52.44 56.09 69.19 C 2012 The Authors, MNRAS 422, 3189–3207 Monthly Notices of the Royal Astronomical Society C 2012 RAS 3196 R. C. Gilmore et al. Figure 5. The history of the EBL in each of our models. The top two plots show the background flux at past redshifts in the WMAP5 fiducial (left) and CCDM (right) models in standard units. Redshifts shown include z=0 (solid), z=0.2 (dotted), z=0.6 (short dashed), z=1 (long dashed), z=1.5 (dot–short dashed), z=2 (dot–long dashed) and z=2.5 (long and short dashed); also see the key in the upper left-hand panel. The bottom two plots show the same quantities, but now evolved to present-day (comoving), allowing easy comparison of the EBL in place at a particular time compared to the total in existence at z=0. 0.5, it is important in comparing different realizations of the EBL that we focus not only on the flux at z=0 but at higher redshifts as well, where behaviour may be quite different in a given model depending on how galaxies evolve. We show how the background develops in our models in two ways in Fig. 5. The top panels show the proper EBL SED from different redshifts in the rest frame, for each of our models. The bottom panels show the comoving EBL at those same redshifts; this is the background that would be seen today if all galaxy emissivity had been shut off below the indicated redshift. It can be seen in the top plots that the EBL photon density was considerably higher in the past at all wavelengths. The most striking increases from present day levels are in the midand far-IR, andintheUV. Complementary to Fig. 5, in Fig. 6 we show how the photons populating the IR EBL at various wavelengths today were produced as a function of redshift. As expected from our knowledge of obscured starbursting galaxies at high redshift, the midand far-IR parts of the EBL came into existence considerably sooner than the photons that are part of the optical–near-IR peak today. Our results are in reasonable agreement with a recent survey of submillimetre galaxies (Devlin et al. 2009), which has found that half of the background radiation at 250μm is produced at z>1.2, with this fraction increasing at longer wavelengths. The results for the WMAP5and CCDM models are qualitatively similar, however, due to earlier star formation in the latter, a greater percentage of photons were in place at a given redshift for all wavebands relative to the WMAP5 model. Because the measurements shown here are unavoidably incomplete at high redshift, the fact that our models overpredict the fraction of light in place at early times is not necessarily in conflict with these results. The rapid increase in flux at all wavelengths with increasing redshift means that the attenuation per unit distance increases a corresponding amount. Therefore, gamma-rays from more distant blazars suffer more attenuation than might be expected from the local EBL flux. In addition, the SED shape of the EBL changes, so a simple z-dependent scaling factor is not sufficient to allow accurate predictions of spectral modification for the more distant sources. C 2012 The Authors, MNRAS 422, 3189–3207 Monthly Notices of the Royal Astronomical Society C 2012 RAS Modelling of the EBL and gamma-ray spectra 3197 Figure 6. Build-up of the present-day background photon population in three different bands (observed frame), showing the fraction of the presentday comoving number density in place at a given redshift. Solid black is our fiducial WMAP5 model, dashed is the WMAP5+fixed and dash–dotted is the model of D11. We have compared against the measurements of Le Floc’h et al. (2009) (blue) at 24µm and Jauzac et al. (2011) at 70 (green) and 160µm (red). For the latter two, the extent of the shaded region is the 1σerror bound. 3.2 Gamma-ray attenuation The process of photon–photon scattering to electron–positron pairs is well understood from quantum electrodynamics. The basic kinematic requirement for this process is that there must be sufficient energy in the centre-of-mass frame of the two-photon system to create the pair. Including the effect of interaction angle as measured in the cosmological frame, this can be written as 2E1E2(1 −cos θ)≥2mec2,(3) where E1and E2are the photon energies and θis the angle of incidence. We are interested here in cases where the target background photon has energies from the far-IR (10−2eV) to the Lyman limit (13.6eV). The corresponding gamma-ray energies are therefore in the TeV or GeV range. We can rewrite equation (3) to define the minimum threshold energy Eth for a background photon to interact with a gamma-ray of energy Eγ, Eth =2m2 ec4 Eγ(1 −cosθ).(4) The cross-section for this process is (Gould & Schreder 1967; Madau & Phinney 1996) σ(E1,E 2,θ)=3σT 16 (1 −β2) ×2β(β2−2) +(3 −β4)ln1+β 1−β,(5) where β=1−2m2 ec4 E1E2(1 −cosθ),(6) and σTis the Thompson scattering cross-section. The cross-section is maximized for centre-of-mass energies of approximately twice the threshold energy 2mec2, and falls approximately as inverse energy for EEth. If we also account for θ,we find that the likelihood of absorption is maximized for photons at about four times the absolute threshold energy, with one factor of 2 from σand another in going from θ=π(‘head-on’ configuration) to the most probable angle of interaction θ≈π/2. If we assume θ=π/2, then we can define the characteristic energy or wavelength for the background photons that will most strongly affect a gamma-ray of energy Eγas Ebg =4m2 ec4 Eγ =1.044 TeV EγeV,(7) or equivalently, λbg =1.188 Eγ TeV μm.(8) Gamma-rays at a rest-frame energy above 1TeV are most attenuated by the nearand mid-IR range of the EBL, while those in the 200 GeV to 1 TeV regime are sensitive to light in the near-IR and optical peak in the EBL SED. Below 200 GeV it is mainly UV photons that have sufficient energy to cause the pair production interaction. Below 19 GeV only background photons with energies above the Lyman limit of 912 Å have sufficient energy to interact at any angle in the rest frame, and there is little attenuation due to the paucity of such photons (Oh 2001; Gilmore et al. 2009). To calculate the optical depth for a gamma-ray observedat energy Eγ, we perform the integral along the line of sight to the target at redshift z, τ(Eγ,z 0)=1 2z0 0 dzdl dz1 −1 du(1 −u) ×∞ Emin dEbg n(Ebg,z)σ(Eγ(1 +z),E bg,θ).(9) Wherewehave Emin =Eth (1 +z)−1=2m2 ec4 Eγ(1 +z)(1 −cosθ) to account for the redshifting of the gamma-ray energy. Here n(Ebg, z) is the proper density of target background photons as a function of energy Ebg and redshift z,anduis shorthand for cos θ.dl/dzis the cosmological line element, equation (2). For nearby sources, z0.05, it is sufficient to use the local EBL density n(Ebg,z=0). However, as we saw in Fig. 5, both the total power and SED of the EBL vary strongly with redshift, and in general it is therefore necessary to understand the evolution of the background to correctly compute gamma-ray opacities. The rapid increase in flux at all wavelengths with increasing redshift to z 2 means that the attenuation per unit proper distance increases a corresponding amount. This effect means that gammas from more distant blazars suffer more attenuation than might be expected from the local EBL flux. In addition, the functional form of the EBL changes, so a simple z-dependent scaling factor is not sufficient to allow accurate predictions of spectral modification for the more distant sources. Using the line-of-sight integral, equation (9), we show in Fig. 7 the optical depth versus gamma-ray energy for a variety of redshifts. A more general way to show EBL attenuation is to plot the ‘attenuation edge’ redshift where the optical depth reaches a certain value as a function of gamma-ray energy, and this is presented out to high redshift for three values of τin Fig. 8. This shows how telescopes with lower energy thresholds will allow us to peer deeper into the universe. See Gilmore et al. (2009) for a similar plot extending to lower gamma-ray energies and higher redshift. C 2012 The Authors, MNRAS 422, 3189–3207 Monthly Notices of the Royal Astronomical Society C 2012 RAS 3204 R. C. Gilmore et al. will not be attainable at elevated temperatures (van Dokkum et al. 2007). Our models lie below the level of direct detection of the absolute background by calculations using data from Hubble Wide Field and Planetary Camera 2 (WFPC2), DIRBE and IRTS (see the discussion for Fig. 4), however, we do find consistency with the newer optical results of Mattila et al. (2011), and also the Pioneer analysis of Matsuoka et al. (2011). The low significance and large error bars on the HST points of Bernstein (2007) mean that these results should not be considered inconsistent with an EBL at the level provided by resolved sources. All of our models are at least 1σbelow the flux from any of the near-IR direct detection calculations we have discussed. Limits from gamma-ray observations shown in Fig. 9 have strongly disfavoured the highest levels at this range. As discussed in Levenson et al. (2007), the present uncertainty in zodiacal light subtraction may be intractable without a new mission to directly study this foreground, such as the proposed ZEBRA experiment (see discussion in Zemcov et al. 2011). While our models are consistent within 1σwith number count measurements by MIPS in the midand far-IR, they are low compared to the DIRBE measurements of the far-IR peak. The zodiacal foreground is a sharply decreasing function of wavelength in this regime, and the DIRBE points are expected to suffer from less systematic error here than in the near-IR, especially at 240 μm, where our models lies beneath the data. However, our IR templates include a somewhat simplified treatment of blackbody emission at these long wavelengths, and we do not expect accurate reproduction in this regime. Fardal et al. (2007) compared the possible range of EBL flux measurements with observations of the fossil mass and SFR history of the universe. As our SAM reproduces these three observables, it is worth discussing our work in the context of their claim that a topheavy or ‘paunchy’ IMF can best fit simultaneously these data. This proposal is based on the argument that the low levels of estimated stellar mass are difficult to reconcile with the present-day EBL flux suggested as being intermediate between the direct detections and the lower limits from resolved galaxies, which is now known to be in conflict with blazar limits (see Fig. 9). Fardal et al. (2007) created three models of the EBL based on all available observational limits. Their minimal model, with total flux of 50nW m−2sr−1,issetby resolved number counts and is similar to our WMAP5andCCDM modelsin theopticaland near-IRoutto theKband. Theirbest-fitting model, based on a compromise between number counts compilations and the HST and DIRBE direct detection measurements, is substantially higher than our CCDM model. The K-band number counts are well measured by a number of surveys (see fig. 14 in SGPD12) which constrain the amount of stellar mass in the nearby universe. Two factors alleviate the discrepancy in our model. Our background fluxes are near the lowest levels considered in Fardal et al. (2007), with total fluxes of 56.09 and 54.91nW m−2sr−1for our CCDM and WMAP5 models, respectively (Table 1), and our global Chabrier IMF produces more high-mass stars than the dietSalpeter considered as the standard by these authors. For a near-IR flux much higher than our models to not overproduce the K-band counts, this flux would have to arise from a high-redshift population of sources unresolvable in our current surveys, which extend to >24 mag. As mentioned in the Introduction, there are reasons why the SFR measures we compare against at high redshift could be biased high. This interpretation favours our fiducial EBL model, which has slightly less flux than Fardal’s lowest model, and is in fairly good agreement with integrated star formation and observed Kcounts (see SGPD12). 5 DISCUSSION The EBL presents one of the primary barriers to extragalactic gamma-ray astronomy with ground-based instruments. Our determination of a fairly low extragalactic background across the optical and nearto mid-IR, supported by convergence with alternative methods such as Franceschini et al. (2008) and D11 as well as new direct measurement techniques and recent limits from gamma-ray experiments in the optical to mid-IR, is an optimistic prediction for the future of the field. The weight of this evidence also increasingly points to an EBL that is well known over this wavelength range, at a level near that of resolved number counts, though many questions remain about its redshift evolution. The ability of currentand next-generation experiments to detect blazars at larger distances is a function of several factors: the luminosity function and spectral evolution of these objects, the effective area (especially at the lowest energies) and duty cycles of these instruments, and the details of the increasingly uncertain non-local EBL at higher redshifts. The field of extragalactic VHE astronomy has grown considerably in the last 5yr, and ongoing progress on the instrumentation front suggests that many new detections may be coming. In Fig. 8, we have placed vertical lines at 50 and 100GeV to facilitate comparison with the gamma-ray attenuation edge curves. In the fiducial WMAP5 model, the universe does not become optically thick (τ>1) to 100 GeV gamma-rays until z∼0.9, and z∼1.6 for an observed energy of 50 GeV. Our fiducial model does predict slightly higher opacities at high redshift compared to the fiducial model of Gilmore et al. (2009); this is due to more UV light escaping high-redshift galaxies in our evolving dust model than in the fixed dust model used in the CCDM prediction. Nonetheless, the EBL does not become a significant barrier to VHE observations at these low energies in either model until redshifts considerably higher than those for which AGN have presently been detected by ground-based instruments. At the redshift of the current most distant confirmed source, 3C279, we find τ=0.38 for an observed energy of 100 GeV, and τ=1at 175 GeV. As mentioned, the approximation of a local EBL in optical depth calculations is only valid for nearby extragalactic observations. At redshifts above ∼0.3, differences in the evolution of star formation and galaxy emissivity begin to have a substantial effect on attenuation; two different EBL models with the same present-day normalization could have widely varying behaviour at these times. For instance, the Kneiske models and Stecker’s fast evolution models have star formation history peaks at a lower redshift than our models predict. In addition to predicting different results for the present-day EBL than our model, the evolution with redshift is also quite a bit different in these cases. As no direct observations of the EBL are possible at non-zero redshift, predicting attenuation from sources past these distances must be made on the basis of models of galaxy evolution, constrained by surveys of luminosity functions at high redshift. Recent surveys of the non-local universe such as DEEP2 and the multiwavelength GOODS and AEGIS surveys have become powerful tools in constraining the EBL at these distances. Observational methods such as Franceschini et al. (2008), Finke et al. (2010) and D11, which fall under types (ii), (iii) and (iv), respectively, in the classification scheme discussed in the Introduction, are complementary to our semi-analytic approach, which is of the first type. Beyond z∼1 to 2, uncertainties in available star formation and luminosity data become large, and theoretical models will continue to be required to understand the production and evolution of the background. C 2012 The Authors, MNRAS 422, 3189–3207 Monthly Notices of the Royal Astronomical Society C 2012 RAS Modelling of the EBL and gamma-ray spectra 3205 The other impact of the shift to higher redshift observations by lower energy-threshold instruments such as the HESS phase-II upgrade and the future Cherenkov Telescope Array (CTA) experiment is the change in the relevant absorbing photon population to UV wavelengths. Our models, and the others we reference in the previous section, predict a rapid falloff in transmission of gamma-rays above 500 GeV for blazars at the redshift of 3C 279. Detecting emission at or above 1 TeV from sources at this distance will require orders-of-magnitude gains in instrument effective area, or observations of flare events with similar increases in output. The energy range of primary interest for these types of sources is going to be 50 to 500 GeV for the next generation of IACTs, plus lower energy data from Fermi. Below 200GeV, it is the UV background that is primarily responsible for absorption. In Gilmore et al. (2009), we addressed the question of this background using four different models varying in high redshift star formation and quasar emissivity; the ‘fiducial’ model in that paper that we consider to be the most likely was based on the CCDM model in this work. The WMAP5 model with evolving dust parameters presented in this work is generally consistent with this model in its predictions for the evolution of UV emission, though it does predict a somewhat stronger UV background at high redshift. One significant weakness of our present approach is our use of templates to describe re-emission by dust at midand far-IR wavelengths, which is relevant to the interpretation of multi-TeV data from nearby blazars. This method makes the assumption that galaxies of a given bolometric luminosity emit light with a similar spectral distribution. As discussed in SGPD12, future progress in understandingthesewavelengthswill likelyrequiremovingpast this assumption. While there is much progress to be made modelling this part of the background SED, new models of dust will only have a substantial effect on our calculation of gamma-ray opacities for the nearest VHE sources. ACKNOWLEDGMENTS AD, RCG and JRP acknowledge support from a Fermi Guest Investigator Grant. RCG and JRP were also supported by NASA ATP grant NNX07AGG4G, NSF-AST-1010033 and NSF-AST0607712. RCG was supported during this work by a SISSA postdoctoral fellowship. RSS and RCG thank UCSC for hospitality during many visits. AD’s work has been supported by the Spanish Ministerio de Educaci´ on y Ciencia and the European regional development fund (FEDER) under projects FIS2008-04189 and CPAN-Ingenio (CSD2007-00042), and by the Junta de Andaluc´ ıa (P07-FQM-02894). We wish to thank Felix Aharonian, David A. Williams, Amy Furniss and A. Nepomuk Otte for helpful discussions concerning gamma–ray attenuation. Fabio Fontanot assisted in comparisons of different dust templates. The TeVCat (http://tevcat.uchicago.edu/) online TeV-source catalogue was used in the preparation of this paper. REFERENCES Abdo A. A. et al., 2010a, ApJ, 715, 429 Abdo A. A. et al., 2010b, ApJ, 723, 1082 Acciari V. et al., 2009a, ApJ, 690, L126 Acciari V. A. et al., 2009b, ApJ, 693, L104 Aharonian F. A., 2001, in Gladysheva O. G., Kocharov G. E., Kovaltsov G. A., Usoskin I. G., eds, Int. Cosmic Ray Conf. Vol. 27, 27th International Cosmic Ray Conference, Hamburg, Germany, 07-15th August 2001. Copernicus Gesellschaft, Berlin, p. 250 Aharonian F. A. et al., 1999, A&A, 350, 757 Aharonian F. et al., 2002, A&A, 384, L23 Aharonian F. et al., 2005a, A&A, 430, 865 Aharonian F. et al., 2005b, A&A, 436, L17 Aharonian F. et al., 2006, Nat, 440, 1018 Aharonian F. et al., 2007a, A&A, 473, L25 Aharonian F. et al., 2007b, A&A, 475, L9 Aharonian F. et al., 2008a, A&A, 477, 481 Aharonian F. et al., 2008b, A&A, 481, L103 Aharonian F. A., Khangulyan D., Costamante L., 2008c, MNRAS, 387, 1206 Albert J. et al., 2006a, ApJ, 639, 761 Albert J. et al., 2006b, ApJ, 642, L119 Albert J. et al., 2006c, ApJ, 648, L105 Albert J. et al., 2007a, ApJ, 654, L119 Albert J. et al., 2007b, ApJ, 662, 892 Albert J. et al., 2007c, ApJ, 666, L17 Albert J. et al., 2007d, ApJ, 667, L21 Albert J. et al., 2007e, ApJ, 669, 862 Albert J. et al., 2008, Sci, 320, 1752 Aleksi´ c J. et al., 2011a, ApJ, 730, L8 Aleksi´ c J. et al., 2011b, A&A, 530, 4 Arnouts S. et al., 2007, A&A, 476, 137 Atwood W. B. et al., 2009, ApJ, 697, 1071 Baldry I. K., Glazebrook K., 2003, ApJ, 593, 258 Barro G. et al., 2009, A&A, 494, 63 Bell E. F., McIntosh D. H., Katz N., Weinberg M. D., 2003, ApJS, 149, 289 Bernstein R. A., 2007, ApJ, 666, 663 Berta S. et al., 2010, A&A, 518, L30 B´ ethermin M., Dole H., Beelen A., Aussel H., 2010, A&A, 512, A78 Blanton M. R. et al., 2003, ApJ, 592, 819 Bouwens R. J., Illingworth G. D., Franx M., Ford H., 2007, ApJ, 670, 928 Bruzual G., Charlot S., 2003, MNRAS, 344, 1000 Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J., StorchiBergmann T., 2000, ApJ, 533, 682 Cambr´ esy L., Reach W. T., Beichman C. A., Jarrett T. H., 2001, ApJ, 555, 563 Caputi K. I. et al., 2007, ApJ, 660, 97 Chabrier G., 2003, PASP, 115, 763 Charlot S., Fall S. M., 2000, ApJ, 539, 718 Chary R. et al., 2004, ApJS, 154, 80 Cogan P., 2008, in Aharonian F. A., Hofmann W., Rieger F., eds, AIP Conf. Ser. Vol. 1085, 4th International Meeting on High Energy Gamma-Ray Astronomy Am. Inst. Phys., New York, p. 403 Cole S. et al., 2001, MNRAS, 326, 255 Cortina J., 2005, Ap&SS, 297, 245 da Cunha E., Charlot S., Elbaz D., 2008, MNRAS, 388, 1595 Dahlen T., Mobasher B., Somerville R. S., Moustakas L. A., Dickinson M., Ferguson H. C., Giavalisco M., 2005, ApJ, 631, 126 Dahlen T., Mobasher B., Dickinson M., Ferguson H. C., Giavalisco M., Kretchmer C., Ravindranath S., 2007, ApJ, 654, 172 Danforth C. W., Keeney B. A., Stocke J. T., Shull J. M., Yao Y., 2010, ApJ, 720, 976 Daniel M. K. et al., 2005, ApJ, 621, 181 Devlin M. J. et al., 2009, Nat, 458, 737 Devriendt J. E. G., Guiderdoni B., 2000, A&A, 363, 851 Devriendt J. E. G., Guiderdoni B., Sadat R., 1999, A&A, 350, 381 Dole H. et al., 2006, A&A, 451, 417 Dom´ ınguez A. et al., 2011, MNRAS, 410, 2556 (D11) Dwek E., Krennrich F., 2005, ApJ, 618, 657 Faber S. M. et al., 2007, ApJ, 665, 265 Fardal M. A., Katz N., Weinberg D. H., Dav´ e R., 2007, MNRAS, 379, 985 Fazio G. G. et al., 2004, ApJS, 154, 39 Finke J. D., Razzaque S., Dermer C. D., 2010, ApJ, 712, 238 Fixsen D. J., Dwek E., Mather J. C., Bennett C. L., Shafer R. A., 1998, ApJ, 508, 123 C 2012 The Authors, MNRAS 422, 3189–3207 Monthly Notices of the Royal Astronomical Society C 2012 RAS 3206 R. C. Gilmore et al. Fortin P., 2008, in Aharonian F. A., Hofmann W., Rieger F., eds, AIP Conf. Ser. Vol. 1085, 4th International Meeting on High Energy Gamma-Ray Astronomy Am. Inst. Phys., New York, p. 565 Franceschini A., 2001, in Harwit M., Hauser M. G., eds, Proc. IAU Symp. 204, The Extragalactic Infrared Background and Its Cosmological Implications. Astron. Soc. Pac., San Francisco, p. 283 Franceschini A., Aussel H., Cesarsky C. J., Elbaz D., Fadda D., 2001, A&A, 378, 1 Franceschini A., Rodighiero G., Vaccari M., 2008, A&A, 487, 837 Frayer D. T. et al., 2006, ApJ, 647, L9 Gabasch A. et al., 2006, A&A, 448, 101 Gardner J. P., Brown T. M., Ferguson H. C., 2000, ApJ, 542, L79 Gilmore R. C., Madau P., Primack J. R., Somerville R. S., Haardt F., 2009, MNRAS, 399, 1694 Gorjian V., Wright E. L., Chary R. R., 2000, ApJ, 536, 550 Gould R. J., Schreder G. P., 1967, Phys. Rev., 155, 1404 Grazian A. et al., 2009, A&A, 505, 1041 Hauser M. G., Dwek E., 2001, ARA&A, 39, 249 Hinton J. A., 2004, New Astron. Rev., 48, 331 Hopkins A. M., Beacom J. F., 2006, ApJ, 651, 142 Hopwood R. et al., 2010, ApJ, 716, L45 Jauzac M. et al., 2011, A&A, 525, A52 Jelley J. V., 1966, Phys. Rev. Lett., 16, 479 Jones D. H., Peterson B. A., Colless M., Saunders W., 2006, MNRAS, 369, 25 Katarzy´ nski K., Ghisellini G., Tavecchio F., Gracia J., Maraschi L., 2006, MNRAS, 368, L52 Keenan R. C., Barger A. J., Cowie L. L., Wang W., 2010, ApJ, 723, 40 Kennicutt R. C., Jr, 1989, ApJ, 344, 685 Kennicutt R. C., Jr, et al., 1998, ApJ, 498, 181 Kneiske T. M., Dole H., 2010, A&A, 515, A19 Kneiske T. M., Mannheim K., Hartmann D. H., 2002, A&A, 386, 1 Kneiske T. M., Bretz T., Mannheim K., Hartmann D. H., 2004, A&A, 413, 807 Kochanek C. S. et al., 2001, ApJ, 560, 566 Komatsu E. et al., 2011, ApJS, 192, 18 Konopelko A., Mastichiadis A., Kirk J., de Jager O. C., Stecker F. W., 2003, ApJ, 597, 851 Konopelko A., Cui W., Duke C., Finley J. P., 2008, ApJ, 679, L13 Lagache G., Puget J.-L., Dole H., 2005, ARA&A, 43, 727 Le Floc’h E. et al., 2005, ApJ, 632, 169 Le Floc’h E. et al., 2009, ApJ, 703, 222 Levenson L. R., Wright E. L., 2008, ApJ, 683, 585 Levenson L. R., Wright E. L., Johnson B. D., 2007, ApJ, 666, 34 Liu H. T., Bai J. M., Ma L., 2008, ApJ, 688, 148 Madau P., Phinney E. S., 1996, ApJ, 456, 124 Madau P., Pozzetti L., 2000, MNRAS, 312, L9 Madau P., Pozzetti L., Dickinson M., 1998, ApJ, 498, 106 Maier G. et al., 2008, in Caballero R., D’Olivo J. C., Medina-Tanco G., Nellen L., S´ anchez F. A., Vald´ es-Galicia J. F., eds, Int. Cosmic Ray Conf. Vol. 3, 30th International Cosmic Ray Conference, ICRC, p. 1457. Available online at http://www.icrc2007.unam.mx/proceedings Malkan M. A., Stecker F. W., 1998, ApJ, 496, 13 Malkan M. A., Stecker F. W., 2001, ApJ, 555, 641 Marchesini D. et al., 2007, ApJ, 656, 42 Matsumoto T. et al., 2005, ApJ, 626, 31 Matsuoka Y., Ienaka N., Kawara K., Oyabu S., 2011, ApJ, 736, 119 Matsuura S. et al., 2011, ApJ, 737, 2 Mattila K., Lehtinen K., Vaisanen P., von Appen-Schnur G., Leinert C., 2011, preprint (arXiv:1111.6747) Mazin D., Goebel F., 2007, ApJ, 655, L13 Mazin D., Raue M., 2007, A&A, 471, 439 Mazin D., Lindfors E., Berger K., Galante N., Prandini E., Saito T. (for the MAGIC collaboration), 2009, arXiv:0907.0366 Montero-Dorta A. D., Prada F., 2009, MNRAS, 399, 1106 Nikishov A. I., 1962, Soviet Phys. JETP, 14, 393 Nilsson K., Pursimo T., Sillanp¨ a¨ a A., Takalo L. O., Lindfors E., 2008, A&A, 487, L29 Oh S. P., 2001, ApJ, 553, 25 Peebles P. J. E., 1993, Principles of Physical Cosmology, Princeton Univ. Press, Princeton, NJ Pei Y. C., Fall S. M., Hauser M. G., 1999, ApJ, 522, 604 Prandini E., Dorner D., Mankuzhiyil N., Mariotti M., Mazin D. (for the MAGIC Collaboration), 2009, preprint (arXiv:0907.0157) Prandini E., Bonnoli G., Maraschi L., Mariotti M., Tavecchio F., 2010, MNRAS, 405, L76 Primack J. R., Bullock J. S., Somerville R. S., MacMinn D., 1999, Astropart. Phys., 11, 93 Primack J. R., Somerville R. S., Bullock J. S., Devriendt J. E. G., 2001, in Aharonian F. A., V¨ olk H. J., eds, AIP Conf. Ser. Vol. 558, High Energy Gamma-Ray Astronomy: International Symposium. Am. Inst. Phys., New York, p. 463 Primack J. R., Bullock J. S., Somerville R. S., 2005, in Aharonian F. A., V¨ olk H. J., Horns D., eds, AIP Conf. Ser. Vol. 745, High Energy Gamma-Ray Astronomy. Am. Inst. Phys., New York, p. 23 Primack J. R., Gilmore R. C., Somerville R. S., 2008, in Aharonian F. A., Hofmann W., Rieger F., eds, AIP Conf. Ser. Vol. 1085, 4th International Meeting on High Energy Gamma-Ray Astronomy. Am. Inst. Phys., New York,p.71 Razzaque S., Dermer C. D., Finke J. D., 2009, ApJ, 697, 483 Reddy N. A., Steidel C. C., Pettini M., Adelberger K. L., Shapley A. E., Erb D. K., Dickinson M., 2008, ApJS, 175, 48 Rieke G. H., Alonso-Herrero A., Weiner B. J., P´ erez-Gonz´ alez P. G., Blaylock M., Donley J. L., Marcillac D., 2009, ApJ, 692, 556 (R09) Rodighiero G. et al., 2010, A&A, 515, A8 Rowan-Robinson M., 2001, ApJ, 549, 745 Schiminovich D. et al., 2005, ApJ, 619, L47 Schirber M., Bullock J. S., 2003, ApJ, 584, 110 Schlegel D. J., Finkbeiner D. P., Davis M., 1998, ApJ, 500, 525 Sheth R. K., Tormen G., 1999, MNRAS, 308, 119 Soifer B. T., Neugebauer G., 1991, AJ, 101, 354 Somerville R. S., Primack J. R., 1999, MNRAS, 310, 1087 Somerville R. S., Primack J. R., Faber S. M., 2001, MNRAS, 320, 504 Somerville R. S., Gilmore R. C., Primack J. R., Dom´ ınguez A., 2012, MNRAS, doi:10.1111/j.1365-2966.2012.20490 (SGPD12) Somerville R. S., Hopkins P. F., Cox T. J., Robertson B. E., Hernquist L., 2008, MNRAS, 391, 481 (S08) Stecker F. W., Scully S. T., 2009, ApJ, 691, L91 Stecker F. W., Malkan M. A., Scully S. T., 2006, ApJ, 648, 774 Stecker F. W., Malkan M. A., Scully S. T., 2007, ApJ, 658, 1392 Superina G. et al., 2008, in Caballero R., D’Olivo J. C., Medina-Tanco G., Nellen L., S´ anchez F. A., Vald´ es-Galicia J. F., eds, Int. Cosmic Ray Conf. Vol. 3, 30th International Cosmic Ray Conference, ICRC, p. 913. Available online at http://www.icrc2007.unam.mx/proceedings Tagliaferri G. et al., 2008, ApJ, 679, 1029 TakeuchiT.T.,IshiiT. T.,HirashitaH.,YoshikawaK.,MatsuharaH.,Kawara K., Okuda H., 2001, PASJ, 53, 37 Tavani M. et al., 2008, Nuclear Instrum. Methods Phys. Res. A, 588, 52 Tavecchio F., Mazin D., 2009, MNRAS, 392, L40 Teshima M. et al., 2008, in Caballero R., D’Olivo J. C., Medina-Tanco G., Nellen L., S´ anchez F. A., Vald´ es-Galicia J. F., eds, Int. Cosmic Ray Conf. Vol. 3, 30th International Cosmic Ray Conference, ICRC, p. 1045. Available online at http://www.icrc2007.unam.mx/proceedings Toller G. N., 1983, ApJ, 266, 79 Treves A., Falomo R., Uslenghi M., 2007, A&A, 473, L17 van Dokkum P., Cooray A., Labb´ e I., Papovich C., Stern D., 2007, in StorrieLombardi L. J., Silbermann N. A., eds, AIP Conf. Ser. Vol. 943, The Science Opportunities of the Warm Spitzer Mission Workshop. Am. Inst.Phys.,NewYork,p.122 Wagner R. M., 2008, MNRAS, 385, 119 Wagner R. et al., 2008, in Caballero R., D’Olivo J. C., Medina-Tanco G., Nellen L., S´ anchez F. A., Vald´ es-Galicia J. F., eds, Int. Cosmic Ray Conf. Vol. 3, 30th International Cosmic Ray Conference, ICRC, p. 889. Available online at http://www.icrc2007.unam.mx/proceedings C 2012 The Authors, MNRAS 422, 3189–3207 Monthly Notices of the Royal Astronomical Society C 2012 RAS Modelling of the EBL and gamma-ray spectra 3207 Wolf C., Meisenheimer K., Rix H.-W., Borch A., Dye S., Kleinheinrich M., 2003, A&A, 401, 73 Wright E. L., 2001, ApJ, 553, 538 Wright E. L., 2004, New Astron. Rev., 48, 465 Wright E. L., Reese E. D., 2000, ApJ, 545, 43 Wyder T. K. et al., 2005, ApJ, 619, L15 Xu C. K. et al., 2005, ApJ, 619, L11 Yang J., Wang J., 2010, A&A, 522, A12 Zemcov M. et al., 2011, preprint (arXiv:1101.1560) This paper has been typeset from a T EX/L A TEX file prepared by the author. C 2012 The Authors, MNRAS 422, 3189–3207 Monthly Notices of the Royal Astronomical Society C 2012 RAS