© ESO, 2017

1. Introduction

The extragalactic background radiation at various electromagnetic frequencies, from radio to gamma-rays, is a fundamental constituent of the Universe, and was demonstrated to permeate it quite uniformly (e.g., Longair 2000). Such radiations have a key role during most of the history of universal expansion and the formation of cosmological structures. Given their ubiquity, radiations are a fundamental source of opacity for the propagation of high-energy cosmic-ray particles and photons throughout space-time (Nikishov 1962; Gould & Schreder 1966).

The most intense and cosmologically relevant among the diffuse radiation components is the Cosmic Microwave Background. As discussed below (Sect. 2), these numerous photons make a wall to the propagation of ultra-energetic particles and photons in the PeV regime (the Greisen–Zatsepin–Kuzmin effect), that has been confirmed by many experiments (see e.g., the Pierre Auger Collaboration 2010).

Another particularly important component of this radiation is the extragalactic background light (EBL) in the wavelength interval between 0.1 and 1000 μm, produced by cosmic sources during the evolutionary history of the Universe. Interactions between these and very high-energy photons from astrophysical sources, and the consequent pair production, produce strong and observable opacity effects (e.g., Stecker et al. 1992, among others). The corresponding high-energy exponential cutoffs are customarily identified in the high- and very-high energy (HE [0.2 to 100 GeV] – and VHE [above 100 GeV]) spectra of several blazars observed with imaging atmospheric Cherenkov telescopes, now operating between a few tens of GeV to tens of TeV (e.g., HEGRA, HESS, MAGIC, see e.g., Dwek & Slavin 1994; Stanev & Franceschini 1998; Aharonian et al. 2006; Albert et al. 2006). Of course, not only photons from cosmic sources interact with EBL photons, but any particles of sufficiently high energy, from cosmic rays to neutrinos, in principle. Therefore, the issue of a high-precision determination of the time-evolving EBL is of central importance.

Several efforts to model the EBL and its time evolution have been published, including those based on physical prescriptions about the formation and evolution of cosmic structures, in particular the semi-analytical models of galaxy formation. This is the so-called “forward evolution” approach, advocated by Gilmore et al. (2012), among others. The alternatives are the more empirical approaches essentially based on local properties of cosmic sources and educated extrapolations back in cosmic time (e.g., Stanev & Franceschini 1998; Pei et al. 1999; Stecker et al. 2006; Franceschini et al. 2008, henceforth AF2008), the “backward evolution” modeling.

Once the background photon densities are defined, the calculation of the cosmic photon-photon opacity and its effects on very high-energy spectra of blazars are standard physical practice (e.g., Stecker et al. 1992). Many papers have compared the predictions for such attenuated spectra with multi-wavelength blazar observations, and the results appear often rather consistent, for both local and high-redshift objects. An interesting example, illustrative of what is currently feasible, is reported in Ackermann et al. (2012), who analyzed years of Fermi Large Area Telescope (LAT) observations of a large sample of blazars over a wide range of redshifts to investigate the effects of the pair-production opacity in various redshift bins, assuming smooth spectral extrapolations towards the highest energies and the EBL evolutionary models by AF2008. These LAT-normalized blazar spectra show indeed cutoff features increasing with redshift just as expected by the model, providing remarkable proof of the overall validity of our understanding of the local EBL photon density and its time evolution.

Similarly, the multi-wavelength high-energy spectra of both BL Lac and Flat-spectrum radio quasars (the Active Galactic Nuclei – AGN – populations emitting at the highest energies, blazars henceforth) are usually fairly well reproduced in terms of standard blazar emission models (like the Synchrotron Self Compton or the External Compton processes), once corrected for photon-photon absorption based on the most accurate EBL models.

There are however still some areas of concern related with observations of the highest energy blazar spectra. One is the reported relative independence (or only moderate dependence) of the observed spectral indices of blazars in the limited redshift interval currently probed at TeV energies, while a faster increase might be expected due to the strong dependence of the opacity on the source distance (De Angelis et al. 2009; Persic & de Angelis 2008). Another related consideration was concerning the detection of variable emission in the energy range 70 to 400 GeV from the flat spectrum radio quasar PKS 1222+216, where such high-energy photons would be expected to be absorbed in the broad line region (Tavecchio et al. 2012).

Another source of possible concern are some apparent upturns of the spectra of blazars at the highest energies, once the observed spectra are corrected for EBL absorption. A characteristic instance would have been shown by the HEGRA observations of Mkn 501, whose EBL-corrected VHE spectrum appeared to show a statistically significant hardening above 10 TeV, as reported for example in Costamante (2012, 2013).

All the above considerations may have important implications for fundamental physics, as they have raised an important question about photon propagation in space. It has been suggested that photons can oscillate into axion-like particles (ALPs), which are a generic prediction of extensions of the Standard Model of elementary particles, in particular the super-string theories. Photons can oscillate with axions in a similar way as massive neutrinos do, which would require an external field, for example a nano-Gauss magnetic field, to compensate for the photon and ALP spin difference. Then a fraction of VHE photons could escape absorption, because ALPs do not interact with ambient photons (De Angelis et al. 2007; Tavecchio et al. 2012). ALPs may even be an important constituent of the cosmos because, while not contributing to dark matter, they are a possible candidate for the quintessential dark energy.

It has also been proposed that photons and particle interactions at very high energies might reveal deviations from predictions of the Standard Model of particle physics. A violation of the Lorentz transformation of special relativity might happen, for example, in the framework of quantum gravity theories (Jacobson et al. 2003), but would only manifest itself at energies comparable to the Planck energy. As suggested by many, this violation of the Lorentz invariance (LIV) might be testable for photons of energies above ~10 TeV interacting with the EBL, by either increasing or decreasing the standard photon-photon opacity. The famous luminous outburst of the local blazar MKN 501 in 1997 has offered the opportunity to test LIV effects in the spectrum above 10 TeV, corresponding to interactions with EBL photons of wavelengths >10 μm, if we consider the relation between interacting photon energies at the maximum of the cross section: (1)Based on the above considerations, several experiments are now being planned or under development, targeting the highest-energy photons, such as the ASTRI mini array (Di Pierro et al. 2013; Vercellone et al. 2013), HAWC (Abeysekara et al. 2013); HiSCORE (Tluczykont et al. 2014), LHAASO (Cui et al. 2014), all well sensitive in the photon energy interval 1–100 TeV, and above. Ultimately, the Cherenkov Telescope Array (CTA) will offer vast, factor 10 sensitivity improvement compared to current instrumentation over a very wide energy range (Acharya et al. 2013).

All this calls for a revision of the EBL modeling, with a particular consideration of the longer wavelength infrared and submillimetric regime. Our previous effort in AF2008 has revealed remarkable success in reproducing absorption spectral features for sources over a wide redshift interval and up to photon energies of a few TeV. However that EBL model was only partly optimized for wavelengths above several μm, including only a preliminary account of the Spitzer MIPS mid-IR data. The present paper is dedicated to reconsidering the issue of the propagation in space-time of the highest-energy cosmic photons. In particular, we update AF2008 by including new very extensive data in the far-IR and submillimeter obtained with the Herschel Space Observatory, in a spectral region where the EBL and its evolution with time cannot be directly measured. The updated model also makes full use of deep survey data from the Spitzer telescope.

Section 2 summarizes the new data needed for the improved EBL model. Section 3 reports on our data-fitting scheme and multi-wavelength spectral corrections. Section 4 describes the updated EBL model and the cosmic photon density. This accounts only for the contributions by known sources, like galaxies, active galactic nuclei and quasars, and does not consider various published attempts to measure the total-light EBL in the near-IR. Section 5 summarizes our new results on the cosmic photon-photon opacities and applications to a few known blazars. Section 6 contains some discussions on our results and prospects for future observations, including some detailed simulations for observations with the Cherenkov Telescope Array (CTA) concerning in particular the effects of a truly diffuse background from primeval sources.

We use in the following a standard cosmological environment with H₀ = 70 km s^-1 Mpc^-1, Ω_m = 0.3, Ω_Λ = 0.7. We indicate with the symbol S₂₄ the flux density in Jy at 24 μm (and similarly for other wavelengths). L_λ/L_⊙ is the νL_ν luminosity in solar units.

2. The new data

With reference to the analysis of AF2008, the new data made available during the recent years for the estimate of the contributions of extragalactic sources to the EBL concern in particular the source’s infrared emission. Instead, the available data in the optical and near-IR portion of the galaxy spectrum have largely remained unchanged. Considering also that the cosmic opacity originating from UV-optical background photons in AF2008 have proved so far quite successful in reproducing the blazar spectra at about 1 TeV or less, we have chosen not to modify the model at wavelengths shorter than 5 μm (corresponding approximately to the spectral coverage of the Spitzer IRAC surveys).

Many of the data at longer wavelengths that we will exploit in the present analysis were already discussed in detail in Franceschini et al. (2010, AF2010). We concentrate here on the new data only, that have mostly come from two infrared observatories in space, the NASA Spitzer telescope and the ESA Herschel observatory.

Fig. 1 Euclidean-normalized differential number counts of extragalactic sources at 24 μm compared with our model fit. The red circles are from the analysis of SWIRE survey data by Shupe et al. (2008), black squares from Papovich et al. (2004). The contribution by type-I AGNs is shown as green dot-dashed line, moderate-luminosity starbursts (the LIRGs) make the cyan short-dash line (type-II AGNs and starbursts are included in the same population on the assumption that in both classes the IR spectrum is dominated by starburst emission). The red long-dashed line corresponds to the population of high-luminosity sources dominating the IR emissivity at high redshifts. The dotted line is the separate contribution of normal spirals, while the continuous line is the total model counts.

Fig. 2 Euclidean-normalized differential number counts of extragalactic sources at 70 μm compared with our model fit. The blue open circles are from Berta et al. (2011), black triangles from Bethermin et al. 2010). Other datapoints are as reported in AF2010. The cyan short-dashed line is the LIRG and the red long-dashed line the ULIRG populations. Black dotted and green dot-dash lines are the normal spiral and type-1 AGN sources (see Sect. 3 for our population decomposition).

2.1. The Spitzer MIPS number counts and luminosity functions

The Spitzer MIPS imaging camera has been used to map several deep sky regions with the 24 μm channel and, to much shallower depths, at 70 μm. The former observations benefited from the excellent instrumental sensitivity and good spatial resolution, while the latter observations were limited in both senses. From the local EBL viewpoint, what matters particularly are the source number counts at both wavelengths that are reported in differential normalized Euclidean units in Figs. 1 and 2. MIPS observations at 70 μm (and those at 160 μm) are limited by sensitivity and source confusion due to the poor angular resolution. Galaxy number counts based on MIPS data by Frayer et al. (2006) and Bethermin et al. (2010) are reported in differential Euclidean-normalized units in Figs. 3. Figure 2 also includes data from the Herschel observatory that will be presented below. All these are compared with the model predictions discussed later in the paper.

Important to note is that both these counts, as well as those discussed in following sections, show a flat, roughly Euclidean behavior at bright fluxes, a well-defined maximum at a given flux density value, and then a rather quick convergence towards faint fluxes. This has the implication that the source contribution to the total EBL intensity at the two wavelengths is mostly resolved at the limiting fluxes of the deepest surveys (see Madau & Pozzetti 2000, for a similar consideration about the optical counts), and the residual contribution by the fainter sources is negligible. We estimate the latter, in any case, from our model fit to the data in Sect. 3, easily extrapolated to the faintest fluxes.

Since we are interested not only in the local EBL intensity, but also in its time evolution, we need information on the source spatial density and emissivity as a function of redshift. The classically used descriptors for this are the redshift-dependent luminosity functions, giving the source number density per comoving volume as a function of luminosity and cosmic time. Various flux-limited Spitzer samples have been used for deriving redshift surveys with good degrees of completeness, all suitable for the calculation of redshift-dependent luminosity functions. Rodighiero et al. (2010) have worked on a complete 24 μm selected sample and used it to calculate galaxy luminosity functions at the rest-frame wavelengths of 15 μm, shorter than the observation wavelengths to minimize the K-corrections for the typical source redshift (z ~ 1). We report these redshift-dependent luminosity functions in Fig. 4, together with the predictions of our evolutionary model in Sect. 3. Functions at longer wavelengths are discussed in the following section.

Fig. 3 Top: euclidean-normalized differential number counts of extragalactic sources at 100 μm from Berta et al. (2010). The cyan short-dashed line are the LIRG and the red long-dashed line the ULIRG populations. Black dotted and green dot-dash lines are the normal spiral and type-1 AGN sources. Bottom: differential number counts of extragalactic sources at 160 μm compared with our model fit. The red open square and blue open circle datapoints are from Berta et al. (2010) and Berta et al. (2011), respectively, while the small green open squares are a re-analysis of the Spitzer/MIPS data by Bethermin et al. (2010). Black triangles are from Magnelli et al. (2013). The lines are our model fits and population decomposition (see Sect. 3).

Fig. 4 Rest-frame 15 μm luminosity functions of IR-selected galaxies in redshift bins, in comoving volume units. The datapoints are from a Spitzer-MIPS sample selected at 24 μm by Rodighiero et al. (2010), and computed at 15 μm to minimize the K-corrections. The units in abscissa are the luminosity νL(ν) at 15 μm in solar luminosities. In ordinate these are the number dN/dlogL of galaxies per cubic comoving Mpc per unit logarithmic interval of νL(ν). The lines represent our multi-population fit to the data as discussed in Sect. 3. In particular, green dot-dashed line: type-I AGNs; cyan short-dashed line: evolving moderate-luminosity star-forming galaxies (LIRG); red long-dashed line: high-luminosity starbursts (ULIRG); lower dotted black line: quiescent spiral population. The upper black continuous line is the total predicted function.

Fig. 5 From top to bottom panel: euclidean-normalized differential number counts of extragalactic sources at 250μm, 350μm, and 500μm, from the Herschel SPIRE Hermes survey by Oliver et al. (2010). Green datapoints come from the P(D) fluctuation analysis by Glenn et al. (2010). Lines are as in the previous figures.

2.2. The Herschel multi-wavelength number counts and luminosity functions

The two imaging instruments onboard Herschel, PACS (Poglitsch et al. 2010), and SPIRE (Griffin et al. 2010), offer capabilities of deep imaging in the far-IR and sub-millimeter at the wavelengths of 70, 100, 160 μm and 250, 350, 500 μm, respectively, with large improvement in sensitivity and angular resolution compared to previous space facilities. SPIRE has carried out cosmological surveys in the submillimeter over tens of square degrees below the confusion limits (the HerMES survey, Oliver et al. 2010, 2012), and over hundreds of square degrees in a contiguous fashion to sensitivities of S₂₅₀ ~ 30–50 mJy (HerMES and H-ATLAS surveys, Eales et al. 2010). PACS observations (the PEP programme, Lutz et al. 2011, and H-GOODS, Elbaz et al. 2011), in particular, have exploited the diffraction-limited imaging capability of the instrument to reduce the confusion noise and ease source de-blending and identification. PACS imaging, however, requires long integrations and has been performed to such faint fluxes only in small areas, such as GOODS and COSMOS.

We report in Fig. 3 number count data at the effective wavelengths of 100 and 160 μm from Berta et al. (2010) and Magnelli et al. (2013), which account for all various corrections of sampling incompleteness, effective area, and fraction of spurious identifications. The 160 μm number counts include also data collected on wider areas by the Spitzer telescope and reported in Bethermin et al. (2010), these latter quite consistent with the Herschel-PACS data.

The SPIRE number count data at longer wavelengths are reported in Fig. 5. Although they are not directly related to our photon-photon opacity estimate, they are relevant to further constrain our model of the IR emissivity of cosmic sources and its evolution with redshift. Because we are interested in estimating absorption effects at energies up to a few tens of TeV, from Eq. (1) and from consideration of the pair-production cross-section, this means that we require knowledge of the EBL from a wavelength of approximately 10 up to approximately 100 μm, that is, the interval where the EBL will never be measured because of the overwhelming brightness of the foreground interplanetary dust emission. For a precise determination of the EBL intensity, we then need to know both the local IR source emissivity and its evolution with redshift, which is constrained by the deep number counts in the far-IR and sub-mm.

As we see, in all cases, such differential normalized counts display a fast increase when going from flux densities of a few hundred mJy to ~10 mJy. The slope of the counts is monotonically increasing with the effective wavelength.

The SPIRE counts are limited at the faint fluxes by the source confusion. A significant extension to fainter limits is possible based on the analysis of the fluctuations of the background integrated intensities (the probability distribution of deflections, P(D), see, for example, Franceschini 1982), as recently discussed in AF2010 and Glenn et al. (2010). An application to the deep Herschel images is reported by Glenn et al. (2010), who analyzed three fields included in the HerMES programme at all three SPIRE bands (250, 350, and 500 μm). The differential number counts determined in this way of the number counts are reported in Fig. 5 (green open squares). These estimates are consistent with those based on individually detected SPIRE sources, and reveal clear evidence for a break in the slope of the differential counts, with a maximum and fast convergence at low flux densities. Obviously, these data from the P(D) analysis are less reliable than those from individual source detections (in red). Nevertheless, they provide us with useful information about the convergence of the SPIRE counts at faint fluxes.

We note that only the green data points in Fig. 5, and the black asterisks and red arrows in Fig. 6 below, are based on this kind of stacking analysis; all others in the present paper are from individually detected sources at the respective wavelengths.

With this extension, the total counts in the three SPIRE bands account for as much as 64, 60, and 43% of the far-infrared background intensity observed by COBE (red arrows in Fig. 6). Indeed, in the wavelength interval covered by the SPIRE observations, the COBE/FIRAS surveys have allowed direct determination of the EBL intensity, after subtraction of the Galactic and Inter-Planetary dust emissions (Puget et al. 1996; Hauser et al. 1998; Lagache et al. 1999; see also AF2008). All this information about the integrated number counts and the direct EBL estimates, when available, is summarized in Fig. 6, including our updated modeling reported as a thick continuous line. The data points and limits appearing in this figure are also reported in the review by Dwek & Krennrich (2013), to which the reader is referred if interested in the tabulated values.

Fig. 6 Comparison of data on the CIRB intensity at long wavelengths and our best-fit model predictions (thick line). Most of the data are as reported and discussed in AF2008. The cyan arrowheads mark the fraction of the CIRB intensity resolved into sources by Herschel PACS (Berta et al. 2010, 2011) at 100 and 160 μm, a fraction corresponding to roughly half of the COBE intensity at these peak wavelengths. Red arrowheads at longer wavelengths mark the CIRB resolved fractions with Herschel SPIRE (Oliver et al. 2010; Glenn et al. 2010). The three upper blue datapoints in the far-IR are from Hauser et al. (1998), the three lower (red) ones are from a re-analysis of the DIRBE data from the all-sky COBE maps by Lagache et al. (1999), the shaded areas from Fixsen et al. (1998, black shade) and Lagache et al. (2004, red shade). The three black asterisks at 70, 100, and 160 μm are the background intensities from PACS resolved sources and from a background fluctuation analysis by Berta et al. (2011). The green data point comes from a background estimate from a MIPS deep survey by Frayer et al. (2006, 2009).

The Herschel cosmological surveys have also allowed us to make fundamental progress towards the determination of the redshift-dependent luminosity functions at long wavelengths. Particularly relevant are the deep multi-wavelength PACS and SPIRE observations of the GOODS-N, GOODS-S and COSMOS fields, including also deep coverage with MIPS at 70 μm. The former two cover an area of approximately 300 square arcminutes each, and the latter an area of approximately 2 sq. deg., all including a vast amount of complementary data in all e.m. bands (Giavalisco et al. 2004; Scoville et al. 2007; Sanders et al. 2007), optical spectroscopic surveys (Lilly et al. 2007) and highly precise photometric redshifts (Ilbert et al. 2010, 2013). Berta et al. (2010, 2011) reported a detailed description of the Herschel source extractions and data catalogs. Further deep IR surveys have been performed in the Extended Groth Strip and Extended Chandra Deep Field South for a total of more than 1000 square arcmins (Magnelli et al. 2009).

Gruppioni et al. (2013) and Magnelli et al. (2009) published independent multi-wavelength luminosity functions at the rest-frame 35 μm wavelength, over a wide redshift interval from local to z ≃ 4, well consistent with one another. These results are reported in Fig. 7 as red datapoints from Gruppioni et al. (2013) and black datapoints from Magnelli et al. (2013), both compared with our model fit.

Luminosity functions at longer wavelengths were also calculated by many authors based on the deep fields (Gruppioni et al. 2010, 2013; Magnelli et al. 2013; Eales et al. 2010). In addition, Herschel SPIRE data over a total of 39 deg² within five high-latitude fields have been used by Marchetti et al. (2016) to compute luminosity functions at 250/350/500 μm, and provide us with a census of the luminosity density in the low-redshift Universe at these wavelengths.

Fig. 7 Rest-frame 35 μm luminosity functions of IR-selected galaxies in redshift bins, in comoving volume units. The black datapoints are from a Spitzer-MIPS sample selected at the effective wavelength of 70 μm by Magnelli et al. (2009), and computed at 35 μm to minimize the K-corrections. Red datapoints are from Gruppioni et al. (2013). The units in ordinate and abscissa are as in Fig. 4. The lines show our multi-population fit to the data (see caption to Fig. 4). The upper black continuous line is the total predicted function best-fitting the luminosity function data. We limit our analysis to z = 1.3 because higher-redshift sources do not contribute significantly to the EBL.

All these data are relevant to constrain our models of the IR source emissivity discussed in the following section.

3. The data-fitting scheme

To interpret the previously discussed statistical data, we have adopted the multi-wavelength modeling scheme of IR galaxy evolution reported in AF2010 and used in AF2008, with few changes and some parameter optimization. Our aim consists of obtaining a phenomenological fitting scheme with the best possible adherence to this large variety of data in order to use it to accurately quantify the source IR emissivity. For this reason, and considering its good ability to reproduce the data, we did not attempt to substantially update our evolutionary model, by, for example, removing our distinction between LIRG and ULIRG populations (see below) that might be seen as somewhat schematic.

3.1. Population components

Our model considers four population components characterized by different physical and evolutionary properties. The first class are normal spirals, dominating the local source populations and assumed to contribute only at z< 1. The presence of this non-evolving population has already been identified in AF2010 in order to explain the observed number counts at bright fluxes and the low-redshift luminosity functions.

A second class, already introduced by Franceschini et al. (2001) among others, is a population of star-forming galaxies of moderate luminosities (LIRG). In order to explain the fainter number counts at all IR wavelengths, this population needs to evolve fast in redshift up to z ~ 1. The LIRG’s evolution rates in both luminosity and number density are slightly modified compared to AF2010 to account for the new Herschel data. The local fraction of the moderately-luminous star-forming population is assumed to be ~10% of the total galaxy population, in line with the observational constraints.

A third evolutionary population introduced by AF2010 to explain the Spitzer 24 μm and SCUBA data at high-redshifts are the ultra-luminous infrared galaxies (ULIRGs), dominating the cosmic IR emissivity above z ≃ 1.5 and, as such, less relevant for our present analysis. This class of sources is essential to reproduce the statistics, particularly the number counts, in the submillimeter (Herschel) and millimeter (JCMT, IRAM, CSO, APEX). We note that, as explained in AF2010, our used terminology for LIRG and ULIRG indicates objects of approximately 10¹¹L_⊙ at z ~ 1 and ~10¹²L_⊙ at z> 1.5, respectively.

Finally, we consider the contributions of AGNs, in particular type-I AGNs that are easily identified with simple combinations of optical-IR colors (Polletta et al. 2006; Hatziminaoglou et al. 2005). The type-II category, instead, is much more difficult to disentangle among the star forming galaxy population, and we do not treat them separately here.

Our choice to distinguish two separate populations of star forming galaxies, LIRGs and ULIRGs, with different luminosity functions and evolutionary histories, was suggested in order to fit the large observational dataset in the IR, that appear to require at least two components, one of moderate luminosity dominating at low z and a second high-luminosity class at higher redshifts. This choice allowed us the flexibility to fit statistical data in the mid-IR (e.g., the 24 μm counts) and in the sub-millimeter simultaneously. Simpler realizations, for example, including only one component of star forming galaxies, would not allow us to achieve acceptable fits (see also Rowan-Robinson 2009).

Fig. 8 Our adopted IR spectra of various galaxy populations. The short-dashed cyan line corresponds to our adopted spectrum for the moderate-luminosity LIRG star-forming population, while the red long-dashed curve is the spectrum of high-luminosity ULIRG sources. In both cases the spectra are similar to that of the prototype star forming galaxy M 82 (in the range from 5 to 18 μm it is precisely the ISOCAM CVF spectrum of M 82). The lower dotted line corresponds to a low-luminosity inactive spiral (see Sect. 3.2 and AF2010), while the upper dotted line is closer to that of ULIRGs. The lower dot-dash green line is the average type-I AGN spectrum. The boundaries of the MIPS 24 μm filter are also shown in the source rest-frames at various source redshifts. The open squares are the average fluxes estimated by Fadda et al. (2010) for a faint sample of LIRGs at z ~ 1 and ULIRGs at z ~ 2. The two blue filled circles at 127 and 182 μm and the red filled datapoint at 100 μm are the average fluxes for LIRGs and ULIRGs from Herschel observations.

3.2. The spectral model

In consideration of the multi-wavelength nature of the data used, a key ingredient in our description of the IR-selected galaxy emission is the spectral model to be adopted for the various considered sources. This is needed to calculate both the K-corrections and the transformations and interpolations of the luminosity functions at various wavelengths, as well as their volume emissivities.

Fig. 9 Data on the comoving IR bolometric luminosity density from 8 to 1000 μm as a function of redshift for the IR-selected galaxy population, compared with our model prediction discussed in Sect. 3. The luminosity density here is expressed in solar luminosities per cubic Mpc. Line types are as in the previous figures. Blue filled squares are estimates of the bolometric IR emissivity based on the Herschel data by Gruppioni et al. (2013) and the cyan circles by Dunlop et al. (2017). Red filled square datapoints are from the analysis of the 24 μm luminosity functions by Rodighiero et al. (2010), green triangles are from Perez-Gonzalez et al. (2005), all expressed in terms of bolometric luminosity density (left ordinate axis). The filled green circles are from Vaccari et al. (2010) and Marchetti et al. (2016). The black dashed line corresponds to the Madau & Dickinson (2014) fit to the star-formation rate density versus redshift, scaled from M⊙/yr to IR luminosity in L⊙ with the factor 1.7 × 10-10, as in Kennicutt (1998) for their adopted Salpeter’s IMF.

Since AGNs of type-I mostly contribute in the mid-IR, and have little influence on the far-IR and sub-mm statistics, we have simply imposed a priori the spectral shape for them. We assume for type-I AGNs a spectral energy distribution (SED) corresponding to an emission model by a pure face-on dusty torus from Fritz et al. (2006), as reported with a green line in Fig. 8.

As for the normal spiral population, we have simplified the treatment by AF2010 of a luminosity dependent spectral shape, by assuming a single constant spectrum, shown as a dotted line in Fig. 8, and dominated by cold dust.

Concerning the LIRGs and ULIRG populations, modeling their IR spectra is critical to achieve good fits to the data. Our adopted SEDs for both classes, that we assumed to be independent of redshift and luminosity, are reported in Fig. 8 as the cyan and red dashed lines. They have been obtained by slightly modifying the spectrum of the starburst galaxy M 82 (see e.g., Polletta et al. 2007): while in the range from 5 to 18 μm they are fixed to the mid-IR spectrum of M 82 observed by ISOCAM (Franceschini et al. 2001), at longer wavelengths the two spectra have been modified so as to best-fit all the multi-wavelength data reported in Sect. 2. The two average SEDs for LIRGs and ULIRGs needed for reproducing the data are not very dissimilar to one another, but the ULIRG spectrum appears to be broader, with a significant enhancement at λ ~ 60 μm and in the submillimeter. This is likely due to significantly different physical conditions in galaxies at different epochs.

Now, it is important to confirm our adopted spectral shapes with further direct observational constraints. Fadda et al. (2010) published ultra-deep mid-IR spectra of a representative and complete sample of 48 infrared-luminous galaxies obtained with the Spitzer IRS spectrograph. Half of these are LIRGs at z< 1 and the other half are ULIRGs at z> 1. The average spectra for the two classes of sources, as obtained by Fadda et al., are reported in Fig. 8 as small open squares, covering the ranges of 4–10 and 5–12 μm. As already remarked by Fadda et al. and shown in the figure, Spitzer IRS observations are in good agreement with our adopted mid-IR spectral shapes for both LIRGs and ULIRGs.

To complement the Spitzer’s data with longer wavelength data, Lo Faro et al. (2013) analyzed deep Herschel SPIRE and PACS photometric data in the GOODS-South field for 10 LIRGs and 21 ULIRGs of the Fadda et al. sample. We show in Fig. 8 the measurement of the average flux for the two subsamples as the blue filled circles for the LIRGs and red filled circles at 100 μm for the ULIRGs. Overall, we find excellent agreement between these average fluxes and our synthetic spectra, which supports our adopted spectral model.

Our model luminosity functions are calibrated to primarily fit those estimated with Herschel in the far-IR and sub-mm, as discussed in Sect. 2. The spectral model is particularly critical to reproduce well LFs in the mid-IR, like those from the IRAS 12 μm survey, and in the mid- and far-IR by Spitzer and Herschel, that are reported above.

3.3. Model to data comparison

A partial comparison of our model fitting with the statistical dataset is reported as thick lines in Figs. 1 to 7. As we see, the fits are from “acceptable” to “excellent” in all cases. The same happens to data at longer wavelengths, that are less relevant to our analysis however.

For completeness, we also summarize in Fig. 9 our best guess evolution model for the time-dependence of the comoving IR galaxy emissivity (black continuous line), compared with existing literature data, where we see a fast increase of it from z = 0 to 1. Our estimated IR luminosity density shown as a thick black line uses all IR data reported in the review by Madau & Dickinson (2014), plus some other independent determinations. The black dashed line in the figure is their best-fit to the star-formation rate density translated back to IR bolometric luminosity with the Kennicutt (1998) recipes (see figure caption). Further to a fair overall match between our analysis and Madau & Dickinson’s, there are some small differences that are apparent in the figure, in particular our estimate is slightly higher than theirs at 0.5 < z < 1.5 (and slightly lower above). This results from the fact that Madau & Dickinson’s fit, trying to intercept both the IR and UV data points, suffers some tension between the two sets of data (the IR density shows a maximum at slightly lower redshifts than the UV, even corrected for dust extinction). Our analysis reproduces more closely the IR data of our interest here.

Fig. 10 Proper number density of EBL photons (energy-weighted) as a function of the energy ϵ. The various curves correspond to different redshifts, as indicated in the figure insert.

Fig. 11 Optical depth by photon-photon collision as a function of the photon energy for sources located at z = 0.003, 0.01, 0.03, 0.1, 0.3, 0.5, 1, 1.5, 2, 2.5, 3, 4, from bottom to top. Black and red lines are from the present work and from AF2008, respectively. For the reader’s convenience, the thick black line marks the z = 1 line.

Fig. 12 Same as in Fig. 11, but including the contributions of CMB photons, producing a fast increase at the highest photon energies. Redshifts increase from bottom to top figure.

Fig. 13 Energies corresponding to optical depth values of τ = 0.1,1 and 10 for photon-photon collisions, as a function of the redshift distance of the source.

The result reported in Fig. 9 is important to further constrain our modeled EBL intensity between 10 and 300 μm .

While we have adopted a standard χ² test to evaluate the goodness-of-fit to individual datasets (like e.g., the number counts in a given waveband or the luminosity function at a rest-frame wavelength and redshift) and to verify the occurrence of possible serious misfits, we have not attempted to run global χ² minimization and automatic χ² searches of the model parameter space to look for degeneracies in the solutions and to provide quantitatively defined uncertainty intervals. This was because of the huge number of datapoints with largely inhomogeneous error estimates, and hidden and unaccounted-for systematic uncertainties in specific data. The latter is particularly the case for the redshift-dependent LFs, for which a reliable error estimate free of systematic uncertainties is virtually impossible. In some other cases, like for the FIRAS-CIRB intensity (see Fig. 6), errors are not even defined except for an uncertainty range, not suited for a χ² test. The consequence is that we will not offer in the following a real discussion of uncertainties and parameter degeneracies. Note that our analysis benefits from an extreme adherence to a vast amount of data, at the cost of an incapacity to provide a well-defined quantitative confidence range. However, considering only the constraints set by the source multi-wavelength number counts, we can estimate a global uncertainty of about ±10% on them, which directly translates onto a ±10% uncertainty on the photon optical depth up to z ≃ 1 and increasing above. The linear relation between number counts, background intensity, and photon-photon optical depth can be seen in Eqs. (5) to (13) in AF2008.

Fig. 14 Top left: photon-photon absorption correction (exp [τγγ]) for the source MKN 501 at z = 0.034, based on our improved EBL model. Bottom left: observed (open black) and absorption-corrected (filled red) spectrum. Data are taken from Aharonian et al. (2001; see also Aharonian et al. 1999). The intrinsic spectrum is fitted by a broken power-law with relatively flat photon spectral indices, as indicated in the insert. Top right: absorption correction for the source MKN 421 at z = 0.03 . Bottom right: observed (open black) and absorption-corrected (filled red) spectrum. Data are taken from Aharonian et al. (2002). The intrinsic spectrum was fitted by a steep broken power-law. Spectral indices are indicated.

Fig. 15 Spectral analysis of VHE observations of the blazar PKS 1441+25. Full cyan squares are from Fermi and open cyan squares are from MAGIC observations by Ahnen et al. (2015). The data show a spectral steepening in the Fermi regime at ϵ ≃ 3 GeV (the best-fit power laws are shown in black and the spectral indices in the insert) and an exponential cutoff in the range covered by MAGIC data. We fit these data with an intrinsic spectrum as for an LSP-ISP blazar, red curve, while the best-fitting green line includes an EBL absorption as from the present paper. Blue data points are the corrected spectrum assuming an optical/near-IR EBL intensity a factor 1.5 larger than our best-guess estimate. The panel on top reports the absorption correction exp [τγγ].

Fig. 16 Similar to Fig. 15, but including simulations of the performances of future observations with the CTA array (50 h observations) and different spectral corrections for photon-photon absorption. Cyan and red datapoints as in the previous figure. Here we assume the existence of a diffuse excess background in the near-IR from very high redshift sources (Population III objects), on top of our estimated EBL in the present paper. The figure includes a simulation of future observations of a similar gamma-ray source observed at the redshift z = 1.8 with the CTA array, to show how sensitive the instrument will be in constraining truly diffuse background signals from un-resolved high-redshift source population.

4. EBL model improvements and the cosmic photon density

Our current modification of the EBL model by AF2008 is based on new data at long wavelengths, λ> 10μm obtained from the Herschel deep and wide-area surveys, improving previous determinations based on the Spitzer data. Indeed, the latter required some substantial extrapolations from the observed mid-IR to longer wavelengths, that are not confirmed by the new far-IR data. Elbaz et al. (2010) and Rodighiero et al. (2010) found, in particular, that the highest-luminosity galaxies over a wide redshift interval have a relatively lower IR luminosity with respect to previous estimates. Consequently, both the new estimated contribution of sources to the EBL intensity at IR wavelengths and the photon number densities at λ> 8 μm are slightly lower than reported by AF2008.

Figure 6 shows our current assessment of the local EBL intensity (background intensity at redshift z = 0) after the above mentioned modifications. We see that all the spectral range from the UV to approximately 8 μm has remained unchanged as a consequence of the unchanged modeling of the galaxy populations at these wavelengths. Instead, the mid-IR portion from ~8 to 40 μm is now lower because of the slightly reduced IR emissivity of galaxies. More significant is the decrease of the whole far-IR and submillimetric peak of dust emission from galaxies; the peak EBL emission at 160 μm, while estimated to exceed 20 nW/m²/ sr by AF2008 based on Spitzer data, is now lowered to νI(ν) ~ 14 nW/m²/ sr, in good agreement with several other independent determinations (e.g., Lagache et al. 1999; Berta et al. 2011).

As already stressed, our work accounts only for the known source contributions, like those from galaxies, active galactic nuclei and quasars. We do not attempt to consider here various published measurements of the total-light EBL in the near-IR based on COBE (Gorjian et al. 2000; Wright & Reese 2000; Cambresy et al. 2001; Hauser & Dwek 2001; Dwek & Krennrich 2005, 2013), and IRTS (Matsumoto et al. 2005; Kashlinsky 2006) observations, and rocket experiments (Zemkov et al. 2014). These measurements suffer substantial uncertainties due to the poorly understood contribution of the Zodiacal Sun-reflected light. The occurrence of a truly diffuse additional component will only be discussed in Sect. 6 as signals potentially reachable by future observations with Cherenkov telescopes.

Unfortunately, the whole interval from 24 to 70 μm was not covered by any one of the space observatory missions during the last 20 yr (ISO, Spitzer, Herschel), and this is reflected in the lack of the corresponding data in Fig. 6. However, thanks to our spectral analysis in Sect. 3.2, the excellent available information at 24, 70, and 100 μm is easily interpolated inside that wavelength range, as we have done for the luminosity functions in Fig. 7. Note that our local EBL model at λ = 70 μm fits well the two published independent estimates in Fig. 6. In conclusion, we are completely confident that we control in an accurate and robust way the background intensity produced by sources over the whole IR domain, both the local EBL and its time evolution.

Figure 10 illustrates our model prediction for the time-dependent number density of diffuse cosmic photons. Again the difference with AF2008 is confined to photon energies ϵ < 0.2 eV. At the energy of 0.01 eV, corresponding to the peak of the photon density produced by cosmic sources (see Fig. 4 in AF2008), this difference is approximately a factor 2.5, in the sense of lower density values from the present analysis. For completeness, in our Fig. 10 we have added the contribution of the Cosmic Microwave Background (CMB) photons, showing up at wavelengths λ > 100 μm and ϵ < 0.01 eV as a very steep rise in photon density with wavelength and redshift.

5. The cosmic photon-photon opacity

Our previously discussed modifications of the local EBL intensity and its evolution impact on the estimate of the cosmic opacity for photon-photon interaction and pair-production, and on the HE and VHE photon horizon. All detailed formalism and calculations of the optical depths as a function of photon energies and cosmic distances are reported in AF2008, to which we defer the reader.

The detailed dependences of the photon-photon optical depth τ_γγ on the photon energy for sources over a range of distances are reported in Fig. 11, and compared there to those calculated by AF2008. As we see, differences due to our modified EBL concern mostly very high-energy observations of cosmic sources, from approximately 1 TeV to few tens of TeV photon energies, in agreement with the usual rule-of-thumb of Eq. (1), where our currently reduced background photon densities at long wavelengths imply reduced opacities. At lower energies the differences with AF2008 become essentially null.

The effects of the inclusion of CMB photons in the optical depth calculations are instead illustrated in Fig. 12. These effects start to be prominent at E_γ ≥ 100 TeV for local sources, and at progressively lower energies at increasing source distances.

The cosmological horizons for different reference values of the photon-photon optical depth τ_γγ and of the photon energy are reported in Fig. 13. At an energy of E₀ ≃ 1 TeV there is a maximum rate of the decrease of E₀ with redshift for a fixed τ because of the maximum in the EBL at λ ~ 1 μm. The flattening at the highest E₀ values is due to the CMB, while that at the highest redshifts is due to the decreased density of photons produced by galaxies at z > 1.

Based on the constraints from the number counts only, we estimate at roughly ±10% the uncertainties on τ(E_γ,z) for z < 1, and larger for z ≥ 1.

Our improved EBL modeling at the long IR wavelengths is particularly relevant for the interpretation of the highest-energy observations of blazars with Cherenkov telescopes. Because of the strong dependence of τ_γγ on distance, such high-energy multi-TeV photons can be better observed from local or low-redshift sources. Data on the two best known such examples, Markarian 501 (Mkn 501, z = 0.034) and Markarian 421 (Mkn 421, z = 0.031), are reported in Fig. 14. In the top panels, our best-guess extinction factors corresponding to the two source distances are reported as a function of the photon energy. In the bottom panels, we compare the observed (black) and extinction-corrected (red) spectral data, together with the corrected best-fit slopes. The data on the former source, particularly relevant, come from a re-analysis of HEGRA atmospheric Cherenkov observations of the famous 1997 flares of Mkn 501, based on algorithms providing improved energy resolution to best constrain the highest energy portion of the spectrum (Aharonian et al. 2001). As for the Mkn 421 data, these were taken during the strong TeV gamma-ray outbursts registered during the observational season of 2000–2001 (Aharonian et al. 2002).

In neither case do we see evidence in the absorption-corrected spectra for excess flux or upturns at the highest energies that might be indicative of an improper spectral correction. Some high-energy upturns for the two sources that were indicated by the analysis in AF2008 (their Fig. 9), and more significantly in Costamante (2013), are not confirmed by the present analysis. While the differences in τ_γγ between this latter and those reported in AF2008 are small at such low redshifts (see Fig. 11), because of the exponential dependence of the absorption correction on τ, these moderately significant upturns are removed.

Aharonian et al. (2002) have performed a systematic comparison of the VHE spectra of the two sources and, based on the fact that they are essentially at the same distance, concluded that the exponential convergence of the two spectra is not consistent with being due to the EBL photon-photon absorption only, and an intrinsic cutoff is required. Our present analysis suggests that the pair-production effect predicted by our re-evaluated EBL intensity is sufficient to explain the VHE spectra for both sources, assuming intrinsic power-law shapes (that in the case of Mkn 421 are one unit, in the spectral index, steeper than for Mkn 501).

6. Discussion and conclusions

6.1. Possible indications for new physics

One of the important topics still open after several published analyses and long debate concerns the detailed spectral forms of the few local blazars whose Cherenkov spectra have been observed up to and above E₀ ~ 10 TeV. These analyses (see e.g., Horns & Meyer 2012; Costamante 2013) reported evidence for substantial upturns in the VHE spectra of Mkn 501 and Mkn 421, once corrected for pair-production absorption according to standard EBL models – like if such absorption corrections would be too large. Even our previous analysis in AF2008 did not exclude such an effect in the two examined local blazars. These inferred spectral upturns at several TeV to tens of TeV energies are not consistent with standard blazar emission models (except assuming rather ad-hoc effects of photon absorption and pile-up). Indications for anomalies in the photon propagation and pair-production have been claimed by comparing the energy distributions for VHE photons received from distant cosmic sources (Horns & Meyer 2012).

As discussed in several papers (e.g., De Angelis et al. 2007), this is an important issue for fundamental physics because, considering the very high energies of photons and particles involved, these effects might reveal deviations from the Standard Model and might require new physics.

The most obvious possibility for increasing the transparency of the universe would be to decrease the photon-photon collision cross-section as a consequence of a violation of the Lorentz invariance, LIV. This may be a natural possibility in the framework of Quantum-Gravity theories, which predict a breakdown of classical physics at energies on the scale of the Planck energy (Jacob & Piran 2008). For example, a modification of the photon dispersion relation might lead to a shift of the energy threshold for pair production at the highest energies, and a decrease of the optical depth at increasing photon energy.

An alternative possibility considered for increasing the photon mean free path is to assume that photons oscillate into and from hypothetical axion-like particles (ALPs) when traveling in the presence of a magnetic field (e.g., De Angelis et al. 2007). During the ALP phase in their path to the Earth, photons would not suffer pair-production effects and would travel undisturbed, then increasing the universal transparency.

There may, however, no longer be a convincing need for new physics emerging from VHE observations of distant sources. This request appears now to be alleviated by the present analysis because, compared with the results in AF2008, our re-analysis of the IR source emissivities at long IR wavelengths, with the corresponding slightly lower inferred EBL intensity values, imply somewhat reduced photon-photon opacities at the highest photon energies at >1 TeV, relaxing the claimed potential inconsistencies.

The pair-production extinction corrections to the observed spectra of Mkn 421 and Mkn 501 reported in Fig. 14 based on our best-guess EBL are consistent with simple spectral extrapolations from the lower photon energies. The high-energy upturns in the corrected spectrum indicated in Costamante (2012) are not apparent here. Our current results are in line with those reported as fiducial models in Primack et al. (2011) and in Stecker et al. (2016). In our case, the modeling of the extragalactic source contribution to the EBL is purely empirical and based on deep photometric imaging data and the source luminosity functions at all wavelengths from UV to the sub-millimeter. Overall, our improved EBL modeling reduces, or even removes, the tension between VHE observations of blazars and standard physical interpretations of their spectra.

6.2. Revising the optical and near-IR EBL (COB) intensity

Interactions between the highest-energy photons produced by cosmic sources and the background light will likely remain a hot topic in coming years. While the pair production effect offers interesting constraints on the EBL intensity in the UV to the near-IR range, the exact level of this background is still subject to significant uncertainties.

Zemcov et al. (2014), for example, report an observed excess signal of near-IR background fluctuations at 1.1 and 1.6 μm on scales of approximately 10 arcmin from a sounding rocket experiment, a signal that may correspond to diffuse light produced by stars in the intergalactic space outside galaxies, potentially implying an optical-NIR EBL around two times higher than our estimated galaxy contribution in Fig. 4.

Abramowski et al. (2013) present a joint analysis of the spectra of a sample of blazars totaling approximately 10⁵ VHE gamma-ray events. Further to the evidence of a clear, highly significant signature of the EBL absorption, the team reports an overall best-guess estimate for the EBL intensity over almost two decades of wavelengths from UV to near-IR with a peak value of νI(ν) ≃ 15 nW/m²/ sr at 1.4 μ, with essentially the same spectral shape of AF2008. This figure is slightly (~40%) higher than our predicted value contributed by resolved sources as in Fig. 6. The same analysis reveals that this excess flux corresponds to a pair-production optical depth larger than reported in AF2008 by a factor ~1.34.

One of the interesting sources where one can look for EBL absorption effects is the high-redshift blazar PKS1441+25 (z = 0.94) that was detected with high significance by MAGIC during an outburst event in 2015 (Ahnen et al. 2015). We show in Fig. 15 results of our analysis of the VHE spectrum of this source. The cyan data-points are the observed fluxes from Fermi and MAGIC, red datapoints are corrected for EBL absorption based on our best-guess optical depth. As can be seen in the figure, the latter are well fit by a theoretical blazar spectrum (red curve, in E²dN/dE units) with a peak energy at approximately 1 GeV, a spectrum that is typical for the Intermediate or Low Synchrotron Peaked (ISP, LSP) objects. High-luminosity blazars at large redshifts tend indeed to be of this category.

In the same figure, we also show (in blue) the MAGIC flux data corrected assuming an EBL intensity in the whole optical/near-IR a factor 1.5 larger than our best-guess estimate, following the suggestion by Abramowski et al., among others. These corrected data would be largely consistent with a power-law extrapolation (black line) of the lower-energy Fermi and MAGIC points, however implying a serious misfit of the highest-energy datum. How seriously this bad fit should be taken, and how physically consistent a power-law continuation of this kind, like the black line, could be (note that this is a high luminosity object, like ISP or LSP kind), cannot be decided with the present data. Furthermore, it is difficult to be conclusive about the exact level of the EBL within the above mentioned boundaries.

6.3. Future prospects

Although atmospheric Cherenkov observatories, joining efforts with the Fermi all-sky surveys, are accumulating VHE photon detections at good rates from sources over a substantial range of redshifts, up to z ~ 1 as we have seen, it seems unlikely that this subject of photon propagation across the universe and EBL characterization will experience significant progress in coming years, further to already published results. This is due to the limited sensitivity, spectral resolution, and energy coverage of current instrumentation, on one side, and the complex interplay between physics and astrophysics on the other, as discussed in the previous subsection.

New-generation instrumentation is clearly needed to make the next step. One such occasion will be offered by the Cherenkov Telescope Array (CTA) in a few years from now (Actis et al. 2011). The full array will offer over factor 10 improvement in sensitivity and spatial and spectral resolution over current telescopes. Among the countless contributions to physics and astrophysics that can be envisaged from such an instrument, we certainly expect substantially improved constraints on the total EBL intensity (sources + diffuse emissions) over a wide interval of wavelengths from UV to the far-IR, as well as relevant constraints on its evolution with cosmic time.

As a particularly significant example of such an opportunity, we have considered using the pair-production absorption effect with an instrument like CTA to constrain the integrated emission of very high-redshift sources responsible for the cosmological re-ionization at redshifts of z ~ 8–10 (that we name Population III, or Pop III sources henceforth). While partly dedicated to this purpose, the new generation space telescope JWST itself will not directly detect them if they are single stars or stellar aggregates.

Instead CTA might see their signatures in the spectra of high-redshift blazars. This also takes advantage of the fact that, while the EBL contribution by known sources, such as galaxies, starts fading away at z> 1 (see Figs. 9 and 10), the Pop III background proper photon number density increases very quickly and proportionally to (1 + z)³ up to z ~ 7–10.

We have then performed simulations of future observations of a few tens of hours with CTA using the MAGIC and Fermi data of the z = 0.94 blazar PKS1441+25 as a reference. We then assumed, on top of our modeled galaxy background, the existence of a truly diffuse signal, like that due to Pop III sources active at very high redshifts. For the latter, we adopted a spectrum as measured by Matsumoto et al. (2005), with a peak at 1.4 μm, but scaled down in flux, following a similar procedure as in AF2008 (see also Raue et al. 2009).

We have simulated CTA observations of PKS1441+25, first of all assuming its measured redshift z = 0.94, a luminosity as observed during the 2015 outburst, and a spectral shape as our best-fit red line in Fig. 15. The simulation results are reported in Fig. 16: the upper continuous red line is the intrinsic source spectrum, while the upper sequence of red datapoints are the simulated data calculated assuming our best-guess EBL model for the photon-photon absorption. The corresponding 1σ error-bars are based on the predicted sensitivity and spectral resolution of the full CTA array for a 50-h integration¹.

The blue data-points correspond to the simulated spectral data corrected for absorption assuming, in addition to the best-guess known-source contribution to the EBL, that of a Pop III component with an intensity equal to only 5% of the IRTS background measured by Matsumoto et al. At the peak intensity wavelength of 1.4 μm, this corresponds to an intensity of 3.5 nW/m²/ sr, added to the known source contribution of 11.4 nW/m²/ sr (see Fig. 6). This Pop III contribution would amount to 1.5 nW/m²/ sr in bolometric units when integrated between 1 and 4 μm, corresponding to a comoving star-formation rate density of approximately 1 M_⊙/ yr/Mpc³ at z ~ 10, according to the analysis of Raue et al. (2009). As we see in Fig. 16, this small truly diffuse component implies a moderate hardening of the (blue) corrected spectrum, however with an overall spectral shape not overly dissimilar to a simple power-law extrapolated from the lower-energy data shown as a black line.

The relative incidence of such a Pop III additional component would become much more evident for a 50-h observation with CTA of the same source, assumed to be observed at z = 1.8 with the same luminosity. To simulate this occurrence, we have scaled the spectrum of PKS1441 to that redshift and performed the same calculation as before.

The lower continuous red line in Fig. 16 is the predicted intrinsic spectrum for z = 1.8, while the lower red datapoints simulate the CTA observation including the EBL absorption for the known-source contribution only. As we see, the source becomes fainter with z and, after application of the EBL absorption, it is detectable by CTA only up to 150 GeV, while the lower-redshift counterpart could be observed up to almost 400 GeV.

We have then corrected back these simulated spectral data assuming the small additional contribution of Pop III sources as above (5% of the IRTS background signal). The resulting absorption-corrected spectrum is shown as green datapoints in Fig. 16. We see that, in this case, the inclusion of a small Pop III signal dramatically increases the τ_γγ for such a high-z source, as shown also by the comparison of the green and blue lines in the upper panel. At the photon energy of 200 GeV, the inclusion of the Pop III excess increases the τ_γγ only by a factor 2 for a source at z = 0.9 (see blue line in Fig. 16 versus the black one in Fig. 15). Instead, for the same source observed at z = 1.8, the inclusion of the small Pop III excess dramatically increases the optical depth, as shown by the green line compared to the blue one in the upper panel of Fig. 16. As mentioned, this is a consequence of the very fast (1 + z)³ increase of the Pop III photon density above z = 1. This exercise illustrates how the low-z known source contributions to the EBL and residual high-z components could be disentangled by comparing CTA observations of blazars at different redshifts.

Our simulation exercise shows that new-generation Cherenkov arrays, with the vast expansion of the VHE cosmological horizon that they will offer, will provide us with new unique opportunities for sampling the history of radiant energy production in the universe.

Acknowledgments

We are glad to acknowledge useful discussions and exchanges with A. De Angelis, Mose’ Mariotti, Michele Doro, Sara Buson, Abelardo Moralejo, and Daniel Mazin among others. This work was mostly supported by the University of Padova.

References

All Tables

All Figures

Fig. 1 Euclidean-normalized differential number counts of extragalactic sources at 24 μm compared with our model fit. The red circles are from the analysis of SWIRE survey data by Shupe et al. (2008), black squares from Papovich et al. (2004). The contribution by type-I AGNs is shown as green dot-dashed line, moderate-luminosity starbursts (the LIRGs) make the cyan short-dash line (type-II AGNs and starbursts are included in the same population on the assumption that in both classes the IR spectrum is dominated by starburst emission). The red long-dashed line corresponds to the population of high-luminosity sources dominating the IR emissivity at high redshifts. The dotted line is the separate contribution of normal spirals, while the continuous line is the total model counts.

In the text

Fig. 2 Euclidean-normalized differential number counts of extragalactic sources at 70 μm compared with our model fit. The blue open circles are from Berta et al. (2011), black triangles from Bethermin et al. 2010). Other datapoints are as reported in AF2010. The cyan short-dashed line is the LIRG and the red long-dashed line the ULIRG populations. Black dotted and green dot-dash lines are the normal spiral and type-1 AGN sources (see Sect. 3 for our population decomposition).

In the text

Fig. 3 Top: euclidean-normalized differential number counts of extragalactic sources at 100 μm from Berta et al. (2010). The cyan short-dashed line are the LIRG and the red long-dashed line the ULIRG populations. Black dotted and green dot-dash lines are the normal spiral and type-1 AGN sources. Bottom: differential number counts of extragalactic sources at 160 μm compared with our model fit. The red open square and blue open circle datapoints are from Berta et al. (2010) and Berta et al. (2011), respectively, while the small green open squares are a re-analysis of the Spitzer/MIPS data by Bethermin et al. (2010). Black triangles are from Magnelli et al. (2013). The lines are our model fits and population decomposition (see Sect. 3).

In the text

Fig. 4 Rest-frame 15 μm luminosity functions of IR-selected galaxies in redshift bins, in comoving volume units. The datapoints are from a Spitzer-MIPS sample selected at 24 μm by Rodighiero et al. (2010), and computed at 15 μm to minimize the K-corrections. The units in abscissa are the luminosity νL(ν) at 15 μm in solar luminosities. In ordinate these are the number dN/dlogL of galaxies per cubic comoving Mpc per unit logarithmic interval of νL(ν). The lines represent our multi-population fit to the data as discussed in Sect. 3. In particular, green dot-dashed line: type-I AGNs; cyan short-dashed line: evolving moderate-luminosity star-forming galaxies (LIRG); red long-dashed line: high-luminosity starbursts (ULIRG); lower dotted black line: quiescent spiral population. The upper black continuous line is the total predicted function.

In the text

	Fig. 5 From top to bottom panel: euclidean-normalized differential number counts of extragalactic sources at 250μm, 350μm, and 500μm, from the Herschel SPIRE Hermes survey by Oliver et al. (2010). Green datapoints come from the P(D) fluctuation analysis by Glenn et al. (2010). Lines are as in the previous figures.
In the text

Fig. 6 Comparison of data on the CIRB intensity at long wavelengths and our best-fit model predictions (thick line). Most of the data are as reported and discussed in AF2008. The cyan arrowheads mark the fraction of the CIRB intensity resolved into sources by Herschel PACS (Berta et al. 2010, 2011) at 100 and 160 μm, a fraction corresponding to roughly half of the COBE intensity at these peak wavelengths. Red arrowheads at longer wavelengths mark the CIRB resolved fractions with Herschel SPIRE (Oliver et al. 2010; Glenn et al. 2010). The three upper blue datapoints in the far-IR are from Hauser et al. (1998), the three lower (red) ones are from a re-analysis of the DIRBE data from the all-sky COBE maps by Lagache et al. (1999), the shaded areas from Fixsen et al. (1998, black shade) and Lagache et al. (2004, red shade). The three black asterisks at 70, 100, and 160 μm are the background intensities from PACS resolved sources and from a background fluctuation analysis by Berta et al. (2011). The green data point comes from a background estimate from a MIPS deep survey by Frayer et al. (2006, 2009).

In the text

Fig. 7 Rest-frame 35 μm luminosity functions of IR-selected galaxies in redshift bins, in comoving volume units. The black datapoints are from a Spitzer-MIPS sample selected at the effective wavelength of 70 μm by Magnelli et al. (2009), and computed at 35 μm to minimize the K-corrections. Red datapoints are from Gruppioni et al. (2013). The units in ordinate and abscissa are as in Fig. 4. The lines show our multi-population fit to the data (see caption to Fig. 4). The upper black continuous line is the total predicted function best-fitting the luminosity function data. We limit our analysis to z = 1.3 because higher-redshift sources do not contribute significantly to the EBL.

In the text

Fig. 8 Our adopted IR spectra of various galaxy populations. The short-dashed cyan line corresponds to our adopted spectrum for the moderate-luminosity LIRG star-forming population, while the red long-dashed curve is the spectrum of high-luminosity ULIRG sources. In both cases the spectra are similar to that of the prototype star forming galaxy M 82 (in the range from 5 to 18 μm it is precisely the ISOCAM CVF spectrum of M 82). The lower dotted line corresponds to a low-luminosity inactive spiral (see Sect. 3.2 and AF2010), while the upper dotted line is closer to that of ULIRGs. The lower dot-dash green line is the average type-I AGN spectrum. The boundaries of the MIPS 24 μm filter are also shown in the source rest-frames at various source redshifts. The open squares are the average fluxes estimated by Fadda et al. (2010) for a faint sample of LIRGs at z ~ 1 and ULIRGs at z ~ 2. The two blue filled circles at 127 and 182 μm and the red filled datapoint at 100 μm are the average fluxes for LIRGs and ULIRGs from Herschel observations.

In the text

Fig. 9 Data on the comoving IR bolometric luminosity density from 8 to 1000 μm as a function of redshift for the IR-selected galaxy population, compared with our model prediction discussed in Sect. 3. The luminosity density here is expressed in solar luminosities per cubic Mpc. Line types are as in the previous figures. Blue filled squares are estimates of the bolometric IR emissivity based on the Herschel data by Gruppioni et al. (2013) and the cyan circles by Dunlop et al. (2017). Red filled square datapoints are from the analysis of the 24 μm luminosity functions by Rodighiero et al. (2010), green triangles are from Perez-Gonzalez et al. (2005), all expressed in terms of bolometric luminosity density (left ordinate axis). The filled green circles are from Vaccari et al. (2010) and Marchetti et al. (2016). The black dashed line corresponds to the Madau & Dickinson (2014) fit to the star-formation rate density versus redshift, scaled from M⊙/yr to IR luminosity in L⊙ with the factor 1.7 × 10-10, as in Kennicutt (1998) for their adopted Salpeter’s IMF.

In the text

	Fig. 10 Proper number density of EBL photons (energy-weighted) as a function of the energy ϵ. The various curves correspond to different redshifts, as indicated in the figure insert.
In the text

	Fig. 11 Optical depth by photon-photon collision as a function of the photon energy for sources located at z = 0.003, 0.01, 0.03, 0.1, 0.3, 0.5, 1, 1.5, 2, 2.5, 3, 4, from bottom to top. Black and red lines are from the present work and from AF2008, respectively. For the reader’s convenience, the thick black line marks the z = 1 line.
In the text

	Fig. 12 Same as in Fig. 11, but including the contributions of CMB photons, producing a fast increase at the highest photon energies. Redshifts increase from bottom to top figure.
In the text

	Fig. 13 Energies corresponding to optical depth values of τ = 0.1,1 and 10 for photon-photon collisions, as a function of the redshift distance of the source.
In the text

Fig. 14 Top left: photon-photon absorption correction (exp [τγγ]) for the source MKN 501 at z = 0.034, based on our improved EBL model. Bottom left: observed (open black) and absorption-corrected (filled red) spectrum. Data are taken from Aharonian et al. (2001; see also Aharonian et al. 1999). The intrinsic spectrum is fitted by a broken power-law with relatively flat photon spectral indices, as indicated in the insert. Top right: absorption correction for the source MKN 421 at z = 0.03 . Bottom right: observed (open black) and absorption-corrected (filled red) spectrum. Data are taken from Aharonian et al. (2002). The intrinsic spectrum was fitted by a steep broken power-law. Spectral indices are indicated.

In the text

Fig. 15 Spectral analysis of VHE observations of the blazar PKS 1441+25. Full cyan squares are from Fermi and open cyan squares are from MAGIC observations by Ahnen et al. (2015). The data show a spectral steepening in the Fermi regime at ϵ ≃ 3 GeV (the best-fit power laws are shown in black and the spectral indices in the insert) and an exponential cutoff in the range covered by MAGIC data. We fit these data with an intrinsic spectrum as for an LSP-ISP blazar, red curve, while the best-fitting green line includes an EBL absorption as from the present paper. Blue data points are the corrected spectrum assuming an optical/near-IR EBL intensity a factor 1.5 larger than our best-guess estimate. The panel on top reports the absorption correction exp [τγγ].

In the text

Fig. 16 Similar to Fig. 15, but including simulations of the performances of future observations with the CTA array (50 h observations) and different spectral corrections for photon-photon absorption. Cyan and red datapoints as in the previous figure. Here we assume the existence of a diffuse excess background in the near-IR from very high redshift sources (Population III objects), on top of our estimated EBL in the present paper. The figure includes a simulation of future observations of a similar gamma-ray source observed at the redshift z = 1.8 with the CTA array, to show how sensitive the instrument will be in constraining truly diffuse background signals from un-resolved high-redshift source population.

In the text