AN EMPIRICAL DETERMINATION OF THE INTERGALACTIC BACKGROUND LIGHT FROM UV TO FIR WAVELENGTHS USING FIR DEEP GALAXY SURVEYS AND THE GAMMA-RAY OPACITY OF THE UNIVERSE

In our previous work the observationally tolerated ranges of photon densities were determined from the LDs, ${\rho }_{{L}_{\nu }}$, themselves, with errors provided by the various authors at wavelengths ranging from the far-UV to the 5 μm wavelength in the near-IR. The LDs are computed by integrating fits to the observationally determined luminosity functions:

$$\begin{eqnarray}&&{\rho }_{{L}_{\nu }}={\int }_{{L}_{\min }}^{{L}_{\max }}{{dL}}_{\nu }\,{L}_{\nu }{\rm{\Phi }}({L}_{\nu };z).\end{eqnarray} \tag{ 1 }$$

Cases where authors did not directly compute the LDs were excluded because properly estimating their error requires knowledge of the covariance of the errors in the fit parameters of the LFs and also knowledge of any observational biases. To provide comprehensive redshift coverage of the LDs, SMS12 and SMS14 made use of continuum colors between the wavelength bands to fill in any gaps.

At wavelengths greater than 5 μm very few studies provide determinations of LDs directly, as most authors are more concerned with calculating the total IR density, integrated over all wavelengths. This is because the total IR density is an observable that is correlated with the star formation rate. Therefore, in this paper, we used observer-given analytic fits to the LFs at various wavelengths. When those fits were not provided, we ourselves fit the observed infrared galaxy LFs at various wavelengths and redshifts, and then use Equation (1) to obtain the LDs.

Galaxy LFs typically have characteristic shapes that are flatter at lower luminosities but fall off more steeply at the highest luminosities. However, at wavelengths greater than 5 μm the Schechter function, commonly used at optical wavelengths, does not provide a good fit to the observed LFs, as its exponential decrease falls off too quickly at the bright end compared with the measured LFs. Most observers instead fit their LFs to a broken power-law function that describes the data better (e.g., Marchetti et al. 2016). The transition region between the power law that holds at low luminosities and that holds at high luminosities is usually referred to as the knee of the LF. The luminosity at the knee is usually designated as L_*. Most of the luminosity density from an LF is produced by galaxies with luminosities in the vicinity of the knee; much fainter galaxies do not contribute much to the LD and the much brighter ones are quite rare. Thus, it is critical to determine the location of the knee in order to estimate the LDs with any accuracy.

For the cases where we were required to determine our own fits, we chose a double power-law fitting function of the form

$$\begin{eqnarray}&&{\rm{\Phi }}(L)=\displaystyle \frac{c}{{L}_{* }\left({\left(\tfrac{L}{{L}_{* }}\right)}^{a}+{\left(\tfrac{L}{{L}_{* }}\right)}^{b}\right)}.\end{eqnarray} \tag{ 2 }$$

Here, the parameters a and b are the indices of the power-law fits in the low-luminosity and high-luminosity ranges, respectively. The overall normalization of the LF is given by the parameter c, which has the dimension of luminosity per Mpc³ dex. To compute the LDs, we integrate Equation (1) over the galaxy luminosities between 4 × 10⁷ ${L}_{{}_{\odot }}$ and 10¹⁴ ${L}_{{}_{\odot }}$.

In order to accurately compute the errors in the LDs as derived from the fit parameters, which is essential for our calculation, we must do so in a way that includes terms involving the off-diagonal elements of the error matrix of the fits so as to account for covariance. In cases where the authors provided there own fits, we have therefore re-derived the fits of the LFs to generate the error matrix, retaining the same choice of fitting function and fixed parameters, if given.

Our goal, as in SMS12 and SMS14, was to compute the observationally tolerated ranges of LDs. This requires that we represent the error on these quantities as best we can. The statistical error is determined by properly propagating the fit parameter errors accounting for covariance, lest we overestimate the error. To compute the total error on each LD, we further added in quadrature an additional systematic error that accounts for cosmic variance. We compute cosmic variance based on the field sizes of the individual studies. For the AKARI Wide Field IR Survey Explorer (WISE) and GOODS fields, this value is typically of the order of 10%. In determining the LDs, for consistency, all observational LFs are scaled to a Hubble parameter of h = 0.7.

Our new additions to our observationally determined LDs extend our coverage of rest frame galaxy photon production from the near-IR to 850 μm in the far-IR, with enough determinations at each wavelength band to span the redshift range $0\leqslant z\leqslant 2\mbox{--}3$. Using published results derived from observations by the Spitzer and Herschel space telescopes, sufficient redshift coverage was found for wavelength bands of $8,12,15,24,35,60,90,$ and 250 μm.

There were two cases where we were required to combine LFs from different observational studies in order to provide enough coverage to discriminate the location of L_* particularly for redshifts greater than 1. At 12 μm we combined the results of Pérez-González et al. (2005) and Rodighiero et al. (2010) to compute LDs in redshift bins centered on redshifts of $1.2,1.6,$ and 2.0. In order to achieve sufficient redshift coverage at $90\,\mu {\rm{m}}$ we combined some of the higher redshift 100 μm data from Lapi et al. (2011) with that of the 90 μm data from Gruppioni et al. (2013). This gains us additional coverage at redshifts of $1.4,2.2,$ and 3.0.

At 160, 350, 500, and 850 μm, LF data only exists for the very nearby redshifts. We therefore assumed that their redshift evolution closely follows that of the 250 μm band (Lapi et al. 2011; Marchetti et al. 2016). We use local LDs calculated in those bands as a normalization to this evolution. This assumption is justified because the emission in this wavelength region is dominated by warm dust. At 160 μm we used the local LF given be Patel et al. (2013). At 350 and 500 μm, we used the combined local LF data of Marchetti et al. (2016) from Herschel/SPIRE, the Herschel/SPIRE estimate of Vaccari et al. (2010), and the Planck satellite data from Negrello et al. (2013). At 850 μm, we computed the local LD from the LF provided by Negrello et al. (2013). Figures 1 and 2 show the resulting derived values for ${\rho }_{{L}_{\lambda }}(z)$ and their errors together with the observationally determined ±1σ confidence bands for the 8, 12, 15, 24 μm bands and those of 35, 60, 90, and 250 μm, respectively.

Zoom In Zoom Out Reset image size

Our calculations of LDs at mid-IR to far-IR wavelengths presented here is an extension of the work done in SMS12 and SMS14. Thus, the results given in those papers for wavelengths less than 5 μm is almost unchanged from the results presented here. However, we have updated the far-UV calculations to include the more recent work from Bouwens et al. (2015) for LDs in the redshift range of 4 to 7. Even though the shape of the far-UV band can affect other bands when filling in for redshift gaps using colors, the overall calculation yields results that are qualitatively the same as those presented in SMS12 and SMS14, as the newer data do not significantly change the general trend.

In order to place 68% upper and lower limits by using observational data on ${\rho }_{{L}_{\nu }}$ we make as few assumptions about the luminosity density evolution as possible. In SMS12 and SMS14 we utilized a robust rational fitting function in the form of a broken power law dependent on $(1+z)$ in order to generate the confidence bands. We take the 68% confidence ranges of these fits in each waveband as the ±1σ confidence bands for the LDs. At redshifts beyond the redshift at the peak of star formation, the LDs decline with redshift. In accord with recent studies of the evolution of rest frame LFs in the UV that trace the star formation rate (Finkelstein et al. 2015), we conservatively assume upper and lower limit power-law functions in redshift to represent the rate of this decline as the highest redshifts. In the upper limit case, we assume a decline proportional to ${(1+z)}^{-2}$. In the lower limit case, we adopt a steeper decline proportional to ${(1+z)}^{-4}$. These assumptions have almost no impact on the derivation of the opacity confidence bands that we determined. Figures 1 and 2 show our results for LDs as a function of redshift at various wavelengths.

Zoom In Zoom Out Reset image size

At wavelengths greater than ∼20 μm, the LDs between the bands can be determined by smoothly interpolating between our observationally based LDs at specific wavelengths, since the spectral energy distributions (SEDs) from galaxies in this wavelength range are smooth modified blackbody spectra produced by dust re-radiation. However, in the 5–20 μm range, the situation is more complex.

In star-forming galaxies, the spectra between 7 and 13 μm are dominated primarily by PAH emission. These PAH molecules are found in very small dust grains in intergalactic media. They absorb the UV photons emitted by hot young O and B stars and reemit them in molecular emission bands in the mid-IR. They are thus a strong signature of active star formation in galaxies (Peeters et al. 2004). In this regard, we note that luminous star forming galaxies at higher redshift have more prominent PAH emission features. The importance of PAH features in the mid-IR at redshifts $0.5\leqslant z\leqslant 2.5$ has been shown by Lagache et al. (2004).

The average SED of nearby star-forming galaxies (from Spoon et al. 2007, see also Smith et al. 2007) is shown in Figure 3 normalized to our best-fit low-redshift LD confidence band. One can see that there is a relative "valley" between 9 and 11 μm. A simple direct interpolation between our 8 and 12 μm bands would therefore obtain an incorrectly high value for the LD in this wavelength range. Since we do not have wavelength coverage in this regime, we take this feature into account by lowering our interpolated LDs by a factor of 3 for the upper limit and a factor of 5 for the lower limit at 10 μm. The factor of 5 is chosen as a lower limit based on the difference between the peak at 8 μm and the depth of the valley near 10 μm from the SED. For the upper limit, we relax the depth of this feature because, while star forming galaxies make up the bulk of the IBL, contributions from other galaxy types have a less pronounced PAH feature. Since, at high redshifts, PAH emission correlates with the star formation rate (Shipley et al. 2016), we assume that our relative PAH shape factors of 3 and 5 coevolve with redshift.

Zoom In Zoom Out Reset image size

Ten years ago we made estimates of the diffuse infrared background when hardly any mid-IR and far-IR LFs had been observationally determined at high redshifts (SMS06). That work was therefore based on the assumptions of a "backward evolution" model. Starting from the well-determined local (z = 0) LF at 60 μm, we assumed that the average locally determined transformations between different mid-IR and far-IR wavelengths applied, unchanged, at all redshifts. The luminosity function used was a double power law, similar to that of Equation (2) with the parameters $a=-1.35,b=-3.6$, and ${L}_{* }=8.5\times {10}^{23}$ W Hz⁻¹ at z = 0 as determined at 60 μm, which was the wavelength for which the most complete galaxy LF existed. We then assumed that the effect of redshift evolution could be taken into account purely by the evolution of L_*, as described in SMS06.

We have compared the predictions of the SMS06 backwards evolution model with the recently observed IR LFs as used in this paper. If we make relatively small improvements in the assumed LF parameters, the backwards evolution model agrees with the new observations surprisingly well. The data favor a slightly flatter LF, with a low-luminosity slope of a = 1.25 and a high-luminosity slope of b = 3.25, provided that we compensate by slightly decreasing the normalization of the LF at L_* and z = 0 to $3.5\times {10}^{-3}$ Mpc⁻³ dex⁻¹. We then match the observed LFs at higher redshifts by assuming that L_* increases as ${(1+z)}^{3.0}$ for $0\lt z\lt 2.0$, with L_* being constant at higher redshifts. In this way, we can obtain a good fit between the backward evolution model and the observational data.

Although this slightly modified backwards evolution reproduces the observed LFs at all redshifts well, there are a few discrepancies, in either direction. The most substantial disagreement is with the 8 μm observations of Huang et al. (2007) at z = 0.15. Below the knee of the LF the observed LF the data are up to 0.3 dex higher than the LF of our backwards evolution model. On the other hand, our best-fitting model overpredicts the 15 μm LFs by up to 0.25 dex at luminosities above the knee as compared with the data of Le Floc'h et al. (2005), Pozzi et al. (2004), and Rodighiero et al. (2010). The likely explanation of both of these discrepancies is that our previously proposed simple model SEDs are based on average 12 μm luminosities as derived from the 60 μm observations. Therefore, they do not include the strong contributions from the PAH bands (See Section 2.2). Thus our SEDs will overpredict the dust continuum on either side of the PAH bump feature as shown in Figure 3. Correspondingly, if we lower the normalization to better fit the shorter and longer wavelengths, we under-estimate the broadband fluxes at 8 μm rest wavelength.

Of course, our new observationally based calculation presented in this paper avoids the problems of backward evolution models, since here we use directly observed LFs at 8, 15, and 24 μm. These data include the effect of the PAH emission features and show their significance.

Figure 8 shows our 68% confidence band computed for the z = 0 IBL (EBL) along with EBL SEDs obtained from the models of Franceschini et al. (2008) and Domínguez et al. (2011). Figure 9 shows our 68% opacity bands for $z=0.1,0.5,1,3$, and 5, in comparison with the opacity curves of Franceschini et al. (2008), Domínguez et al. (2011). We note that Franceschini et al. (2008) used data only up to 8 μm that were available at the time and did not include a PAH component in their model. The model of Domínguez et al. (2011) assumes a redshift evolution at redshifts greater than ∼1 that follows the evolution in the K band given by Cirasuolo et al. (2010).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In our previous papers (SMS12 and SMS14), we presented observationally based results for the IBL as a function of wavelength and redshift for wavelengths below of 5 μm. Based on those results, we computed the γ-ray opacity of the universe caused by electron–positron pair production up to a γ-ray energy of $1.6/(1+z)$ TeV. In this paper, we have extended our determinations of the IBL within 68% confidence bands. This determination defines upper and lower limits on the IBL to longer wavelengths that extend to 850 μm. Our model-independent results are based on observationally derived LFs from recent galaxy survey data from both local and high-redshift surveys. These data include results from Spitzer, Herschel, and Planck. We then use these results to calculate the opacity of the universe to γ-rays out to PeV energies. In doing so, we also take account of the redshift dependence of interactions of γ-rays with photons of the 2.7 K CBR (Stecker 1969), since the opacity from interactions with CBR photons dominates over that from interactions with IBL photons at the higher γ-ray energies and redshifts.

In Figure 10, we give an energy-redshift plot showing the highest-energy photons from extragalactic sources as a function of redshifts as determined from Fermi data (Abdo et al. 2010) plotted with our 68% confidence band for $\tau =1$.

Zoom In Zoom Out Reset image size

Our direct results on the IBL are consistent with those from complimentary γ-ray analyses using observations from the Fermi-LAT γ-ray space telescope and the H.E.S.S. air Čerenkov telescope. Figure 11 indicates how well our opacity results for z = 1 overlap with those obtained by the Fermi collaboration (Ackermann et al. 2012). Our results are also compatible with those obtained from higher-energy γ-ray observations using H.E.S.S. (Abramowski et al. 2013). This overlap of results from two completely different methods strengthens confidence that both techniques are indeed complimentary and supports the concept that the spectra of cosmic γ-ray sources can be used to probe the IBL (Stecker et al. 1992).

Zoom In Zoom Out Reset image size

Thus, we find no evidence for modifications of γ-ray spectra by processes other than absorption by pair production, either by cosmic-ray interactions along the line of sight to the source (Essey & Kusenko 2014) or line-of-sight photon–axion oscillations during propagation (e.g., De Angelis et al. 2007; Meyer & Horns 2013). In this regard, we note that the Fermi Collaboration has very recently searched for irregularities in the γ-ray spectrum of NGC 1275 that would be caused by photon–axion oscillations and reported negative results (Ajello et al. 2016).

We conclude that modification of the high-energy γ-ray spectra of extragalactic sources occurs dominantly by pair production interactions of these γ-rays with photons of the IBL. They therefore support the concept of using the future Čerenkov Telescope Array instruments to probe the cosmic background radiation fields at infrared wavelengths. This method can be used in conjunction with future deep galaxy survey observations using the near-infrared and mid-infrared instruments onboard the James Webb Space Telescope.