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

# 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

## Abstract

We have previously calculated the intergalactic background light (IBL) as a function of redshift from the Lyman limit in the far ultraviolet to a wavelength of 5 m near infrared range, based purely on data from deep galaxy surveys. Here we utilize similar methods to determine the mid- and far-infrared IBL from 5 m to 850 m. Our approach enables us to constrain the range of photon densities, by determining the uncertainties in observationally determined luminosity densities and spectral gradients. By also including the effect of the 2.7 K cosmic background photons, we determine upper and lower limits on the opacity of the universe to -rays up to PeV energies within a 68% confidence band.

Our direct results on the IBL are consistent with those from complimentary -ray analyses using observations from the Fermi -ray space telescope and the H.E.S.S. air Čerenkov telescope. Thus, we find no evidence of previously suggested processes for the modification of -ray spectra other than that of absorption by pair production alone.

diffuse radiation – galaxies: observations – gamma-rays: theory

## 1 Introduction

We have previously employed a fully empirical approach to calculating the intergalactic background light (IBL) to wavelengths up to 5 m by using observational data from deep galaxy surveys (Stecker, Malkan & Scully 2012 (hereafter SMS12); and Scully, Malkan & Stecker 2014 (hereafter SMS14)); for a similar approach, see also Helgason & Kashlinsky (2012) and Khaire & Srianand (2015). Here we extend our previous results from SMS12 and SMS14, determining both the intergalactic background light (IBL) at longer wavelengths out to 850 m and the subsequent -ray opacity of the Universe out to multi-TeV energies. The spectra of astrophysical sources of such high energy -rays are being studied by ground-based air Čerenkov telescopes and linked air Čerenkov telescope arrays.

We accomplish the goals of this paper by using very recent deep galaxy survey data at far infrared wavelengths where galaxy emission is produced by dust re-radiation rather than starlight. We also include the -ray opacity from photons of the 2.7 K cosmic background radiation (CBR). This enables us to extend our calculations of the -ray opacity of the Universe to energies much greater than a TeV.

Observations at wavelengths greater than 24 m, have been covered by the Multiband Imaging Photometer on the Spitzer space telescope (MIPS), now being dramatically advanced by the availability of data from the Photoconductor Array Camera and Spectrometer (PACS) and Spectral Photometric Imaging Receiver (SPIRE) instruments on the Herschel space telescope, the Planck space telescope, and by ground-based observations from the Atacama Large Millimeter Array (ALMA) and observations from the Balloon-borne Large Aperture Submillimeter Telescope (BLAST). By using only empirical data rather that modeling, our approach is, by definition, model independent. We use published luminosity functions and interpolations of luminosity densities between observed wavebands, observations from high-redshift galaxy surveys being now sufficient for this purpose.

Most past approaches to determining the IBL and its present-day spectral energy distribution, referred to as the ”extralgalactic background light (EBL)”, require assumptions about how the galaxy luminosity functions evolve (e.g. Malkan & Stecker 1998, 2001; Kneiske, Mannheim & Hartmann, 2002: Stecker, Malkan & Scully 2006 (SMS06); Franceschini et al. 2008; Finke, Razzaque, & Dermer, 2010; Kneiske & Dole 2010). Other approaches make use of semi-analytic models that require various assumptions regarding galaxy evolution, stellar population synthesis modeling, or star formation rates the properties of and dust attenuation, particularly for redshifts greater than 1 (e.g., Salamon & Stecker 1998; Gilmore et al. 2009; Somerville et al. 2012; Dominguez et al. 2011; Inoue et al. 2013; Driver et al. 2016). In contrast, our observationally-based approach is superior to model-based methods, since it enables a determination of both the IBL and its observational uncertainties. This is because we use observationally determined errors. Approaches that rely on modeling cannot determine such uncertainties.

Observational studies of blazar -ray spectra have also been used to probe the IBL (Ackermann et al. 2012; Abramowski et al. 2013; Biteau & Williams 2015) through its opacity effect caused by electron-positron pair production. This approach was originally suggested by Stecker, De Jager & Salamon (1992) when the infrared EBL was unknown. Our present method of determining the expected -ray opacity of the Universe by first using observational data to determine the IBL and then calculating its effect on -ray spectra of extragalactic sources complements the technique of using direct -ray observations. This is because the later approach requires a thorough knowledge of the intrinsic (unabsorbed) emission spectra of the -ray sources. This requirement introduces unavoidable uncertainties.

Our final results give the -ray opacity as a function of energy and redshift to within a 68% confidence band that is based on observational data. This opacity is solely caused by pair-production interactions of IBL photons with extragalactic -rays . Thus, a direct comparison of this effect with -ray spectra of extragalactic sources enables an assessment of possible additional spectral modifications. One such possible spectral modification has been suggested to be caused by secondary -ray production from cosmic-ray interactions along the line of sight to the source (Essey & Kusenko 2014). Another such modification might be caused by photon-axion oscillations during propagation from the source to the Earth (e.g., De Angelis et al. 2007; Mayer & Horns 2013).

In Section 2 we give the details of our calculations of upper and lower limits on IR galaxy luminosity densities for wavelengths between 5 m and 850 m as a function of redshift, paying particular attention to the effects of the emission spectra of polycyclic aromatic hydrocarbons (PAHs). In Section 3 we calculate the EBL based on the results from Section 2 and compare it with both other work and observational limits. In Section 4 we compute the resulting opacity of the universe to -rays out to a redshift of z = 5, including the total opacity from interactions with both IBL and thermal cosmic background radiation (CBR) photons. In our conclusion section, Section 5, we compare our results with those obtained by other methods in some of the papers mentioned above. We also discuss the implications of our results.

## 2 Calculating the Infrared IBL

We have previously calculated the intergalactic background light (IBL) as a function of redshift in the far ultraviolet to near infrared range, based purely on data from deep galaxy surveys (SMS12; SMS14). Here we utilize similar methods to extend the calculation into the mid- and far infrared IBL out to a wavelength of 850 m.

### 2.1 Determining the Luminosity Densities from Empirical Luminosity Functions

In our previous work the observationally tolerated ranges of photon densities were determined from the luminosity densities (LDs), , 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:

 ρLν=∫LmaxLmindLνLνΦ(Lν;z) (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 luminosity functions (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 luminosity densities (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 holding at high luminosities is usually referred to as the knee of the LF. The luminosity at the knee is usually designated as . 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:

 Φ(L)=cL∗((LL∗)a+(LL∗)b) (2)

Here, the parameters and 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 , which has the dimension of luminosity per Mpcdex. To compute the LDs, we integrate equation (1) over the galaxy luminosities between 4 10 L and 10 L.

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 SSM14, was to compute the observationally tolerated ranges of luminosity densities. 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 order of 10%. In determining the LDs, for consistency, all observational LFs are scaled to a Hubble parameter of

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 . Using published results derived from observations by the Spitzer and Herschel space telescopes, sufficient redshift coverage was found for wavelength bands of and 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 particularly for redshifts greater than . At 12 m we combined the results of Perez-Gonzalez et al. (2005) and Rodighiero et al. (2010) to compute LDs in redshift bins centered on redshifts of and . In order to achieve sufficient redshift coverage at m we combined some of the higher redshift m data from Lapi et al. (2011) with that of the m data from Gruppioni et al. (2013). This gains us additional coverage at redshifts of and .

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 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.

Our calculations of LDs at mid-IR to far-IR wavelengths presented here is an extension of the work done in SMS12 and SSM14. 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 to . 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 SSM14, 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 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 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 . In the lower limit case, we adopt a steeper decline proportional to . 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.

### 2.2 Taking account of PAH emission

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 and 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, Spoon & Tielens 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 has been shown by Lagache et al. (2005).

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.

## 3 The Extragalactic Background Light

At wavelengths above 100 m the FIRAS and DIRBE instruments aboard the Cosmic Background Explorer (COBE) have measured the EBL (Fixsen et al. 1998; Lagache et al. 1999). This diffuse background has now been in large part resolved by recent galaxy count studies using the Herschel telescope (Berta et al. 2010; Béthermin et al. 2012; Viero et al. 2015) along with ground based studies using ALMA (Fujimoto et al. 2016) and BLAST (Marsden et al. 2011). These more recent studies strongly support the COBE results.

However, there are no direct measurements of the EBL in the infrared range below 100 m owing to the predominance of the foreground radiation of Zodiacal light from interplanetary dust re-radiation. The flux of Zodiacal light is approximately two orders-of-magnitude larger than the EBL flux (Spiesman et al. 1995). In this region only lower limits obtained from galaxy counts exist.

Using our calculations of the LDs in the mid-IR and far-IR as a function of wavelength, we have constructed a 68% confidence band for the spectral energy distribution of the cosmic diffuse infrared background light (the IBL at ). Figure 4 shows our new results, combined with our previous results at shorter wavelengths, taken from SMS12 and SMS14. The light shaded band shows the maximum effect of PAH emission. It can be seen that taking account of the details of the PAH spectrum does not significantly affect our EBL results. Figure 4 also shows the observational lower limits on the EBL obtained from galaxy counts (in blue), extrapolations of mid-IR galaxy counts from Spitzer, and direct measurements (in black).

## 4 The Optical Depth from γ+γ→e++e− Interactions with IBL and 2.7 K CBR Photons

The co-moving radiation energy density for wavelength at redshift , , where , is the time integral of the co-moving luminosity density ,

 uν(z)=∫zmaxzdz′ρν′(z′)dtdz(z′), (3)

where and is the redshift corresponding to initial galaxy formation (Salamon & Stecker 1998), and

 dtdz(z)=[H0(1+z)√ΩΛ+Ωm(1+z)3]−1, (4)

with and .

The upper and lower limits on our co-moving energy densities, derived using equation (3), are shown in Figures 5 and 6.

In calculating the -ray opacities we use the relations for the photon energy, , and the photon density, .

The cross section for photon-photon annihilation to electron-positron pairs was first calculated by Breit & Wheeler (1934) as a solid result of quantum electrodynamics. The threshold for this interaction follows from the Lorentz invariance of the square of the four-momentum vector. This reduces to the square of the threshold energy in the c.m.s., , that is necessary to produce twice the electron rest mass:

 s=2ϵEγ(1−cosθ)=4m2e (5)

The Lorentz invariant quantity given by has been determined to hold to better than one part in (Stecker & Glashow 2001; Jacobson, Liberati, Mattingly & Stecker 2004).

The optical depth for -rays caused by electron-positron pair production interactions with photons of the stellar radiation background is given by

 τ(E0,ze)=c∫ze0dzdtdz∫20dxx2∫∞0dν(1+z)3[uν(z)hν]σγγ[s=2E0hc/λx(1+z)], (6)

where is the co-moving energy density of the photon field, (Stecker, De Jager, & Salamon  1992).

In equations (5) and (6), is the observed -ray energy at redshift zero, is the wavelength at redshift , is the redshift of the -ray source at emission, ,
being the angle between the -ray and the soft background photon, and the pair production cross section is zero for center-of-mass energy , being the electron mass. Above this threshold, the pair production cross section is given by

 σγγ(s)=316σT(1−β2)[2β(β2−2)+(3−β4)ln(1+β1−β)], (7)

where is the Thompson scattering cross section and (Jauch & Rohrlich 1955).

It follows from equation (5) that the pair-production cross section has a threshold at .

The optical depth of the universe to the CBR is given by

 τCBR=5.00×105√1.11PeVEγ∫z0dz′ (1+z′) e−(1.11PeVEγ(1+z′)2)√ΩΛ+Ωm(1+z′)3 (8)

(Stecker 1969; SMS06).

Figure 7 shows the 68% confidence opacity bands for interactions with IBL photons given for sources at and , calculated using the methods described above, along with the opacity produced by interactions of -rays with photons of the 2.7 K cosmic background radiation. Note that at the higher energies and redshifts where the opacity is dominated by interactions with CRB photons, the uncertainty band becomes a very thin line, since the CBR-dominated opacity is exactly determined by equation (8).

## 5 Discussion and Conclusions

### 5.1 Comparison with our previous backward evolution model

Ten years ago we made estimates the diffuse infrared background when hardly any mid-IR and far-IR luminosity functions 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 60m, 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 and W/Hz at 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 account purely by the evolution of , 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 and a high-luminosity slope of , provided that we compensate by slightly decreasing the normalization of the LF at and to Mpc dex. We then match the observed LFs at higher redshifts, by assuming that increases as for , with 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 8m observations of Huang et al. (2007) at . 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 Pozzi et al. (2004), Le Floch et al (2005), 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 12m 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 8m 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.

### 5.2 Our present results and comparison with other work

Figure 8 shows our 68% confidence band computed for the 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 and , in comparison with the opacity curves of Franceschini et al. (2008) and 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 -band given by Cirasuolo et al. (2010).

### 5.3 Conclusion

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 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 luminosity functions 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 cosmic background radiation (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 .

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 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).

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; Mayer & 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 aboard the James Webb Space Telescope.

### References

1. Abdo, A., Ackermann, M., Ajello, M. et al. 2010, ApJ 723, 1082
2. Abramowski, A., Acero, F., Aharonian, F. et al. 2013, A&A 550, A4
3. Ackermann, M., Ajello, M., Allafort, A. et al. 2012, Science 338, 1190
4. Ajello, M., Albert, A., Anderson, B. et al. 2016, Phys. Rev. Letters 116, 161101
5. Babbedge, T.S.R., Rowan-Robinson, M., Vaccari, M. et al. 2006, MNRAS 370, 1159
6. Béthermin, M., Le Floc’h, E., Ilbert, O. et al. 2012, A&A 542, A58
7. Berta, S., Magnelli, B., Lutz, D. et al. 2010, A&A 518, L30
8. Bouwens, R.J., Illingworth, G. D., Oesch, P. A. 2015, ApJ 803, 34
9. Biteau, J. & Williams, D. A. 2015, ApJ 812:60
10. Breit, G and Wheeler, J. A. 1934, Phys. Rev. 46, 1087
11. Caputi, K.I., Lagache, G., Yan, L. et al. 2007, ApJ 660, 97
12. Cirasuolo, M. McLure, R. J., Dunlop, J. S. et al. 2010, MNRAS401, 1166
13. De Angelis, A., Roncadelli, M. & Mansutti, O. 2007, Phys. Rev. D 76, 121301
14. Domínguez, A., Primack, J. R., Rosario, D. J. et al. 2011, MNRAS 410, 2556
15. Driver, S.P., Andrews, S.K., Davies, L.J. et al. 2016, arXiv:1605.01523.
16. Dye, S., Dunne, L., Eales, S. et al. 2010, A&A 518, L10
17. Dwek, E. & Krennrich, F. 2013, Astropart. Phys. 43, 112
18. Eales, S., Raymond, G.; Roseboom, I. G. et al. 2010, A&A 518, L23
19. Essey, W. & Kusenko, A. 2014, Astropart. Phys. 57, 30.
20. Fazio, G. G. & Stecker, F. W. 1970, Nature 226, 135
21. Finke, J. D., Razzaque, S. & Dermer, C. D. 2010, ApJ 712, 238
22. Finkelstein, S. L. Ryan, R. E., Jr., Papovich, C. et al. 2015, ApJ 810:71
23. Fixsen, D.J. Dwek, E., Mather, J. C., Bennett, C. L. & Shafer, R. A., ApJ 508, 123 (1998)
24. Franceschini, A., Rodighiero, G. & Vaccari, M. 2008, A&A,  487, 837
25. Fujimoto, S., Ouchi, M., Ono, Y., et al. 2016, ApJ Supp. 222, 1
26. Gilmore, R. C. and Madau, P., Primack, J. R. et al. 2009, MNRAS 399, 1694
27. Goto, T., Oi, N., Ohyama, Y., et al. 2015, MNRAS 452, 1684
28. Gruppioni, C., Pozzi, F., Rodighiero, G. et al. 2013, MNRAS 432, 23
29. Huang, J.-S., Ashby, M. L. N., Barmby, P. et al. 2007, ApJ 664, 840
30. Helgason, K. & Kashlinsky, A. 2012, ApJ 758, L13
31. Inoue, Y., Inoue, S., Kobayashi, M. A. R. et al. 2013, ApJ 768:197
32. Jacobson, T., Liberati, S. Mattingly, D. & Stecker, F.W. 2004, Phys. Rev. Letters 93, 021101
33. Jauch, J. M. & Rohrlich, F. 1955, The Theory of Photons and Electrons (Cambridge, MA: Addison-Wesley)
34. Khaire, V. & Srianand, R. 2015, ApJ 805:33
35. Kneiske, T. M., Mannheim, K. & Hartmann, D. H. 2002, A&A 386, 1
36. Kneiske, T. M. & Dole, H. 2010, A&A 515, A19
37. Lagache, G., Abergel, A., Boulanger, F. et al. 1999, A&A 344, 322
38. Lagache, G. Dole, H., Puget, J.-L. et al. 2004, ApJ Suppl. 154:112
39. Lapi, A. et al. González-Nuevo, J., Fan, L. et al. 2011, ApJ 742:24
40. Le Floc’h, E., Papovich, C., Dole, H. et al. 2005, ApJ 632, 169
41. Malkan, M. A. & Stecker, F. W. 1998, ApJ 496, 13
42. Malkan, M. A. & Stecker, F. W. 2001, ApJ 555, 641
43. Marchetti, L., Vaccari, M., Franceschini, A.et al. 2016, MNRAS 456, 1999
44. Marsden, G., Chapin, E. L., Halpern, M. et al. 2011, MNRAS 417, 1192
45. Meyer, M. & Horns, D. 2013, Phys. Rev. D 87, 035027
46. Negrello, M., Clemens, M., Gonzalez-Nuevo, J. et al. 2013, MNRAS 429, 1309
47. Patel, H., Clemens, D.L., Vaccari, M., et al. 2013, MNRAS 428, 291
48. Peeters, E., Spoon, H.W.W. & Tielens, A.G.G.M. 2004, ApJ 613, 986
49. Pérez-González, P.G., Pablo G., Rieke, G. H. et al. 2005, ApJ 630, 82
50. Pozzi, F., Gruppioni, C., Oliver, S. et al. 2004, ApJ 609, 122
51. Rodighiero, G., Vaccari, M., Franceschini, A. et al. 2010, A&A 518, A8
52. Salamon, M. H. & Stecker, F. W. 1998, ApJ 493, 547
53. Sano, K., Kawara, K., Matsuura, S. et al. 2016, ApJ 818, 72
54. Scully, S.T., Malkan, M.A. and Stecker, F.W. 2014 (SMS14), ApJ 784:138
55. Shipley, H.V., Heath V.; Papovich, C. et al. 2016, ApJ 818:60
56. Smith, J.D.T., Draine, B. T., Dale, D. A. et al. 2007, ApJ 656, 770
57. Smith, D.J.B., Dunne, L., da Cunha, E. et al. 2012, MNRAS 427, 703
58. Somerville, R. S., Gilmore, R. C., Primack, J. R. & Domínguez, A. 2012, MNRAS 423, 1992
59. Spoon, H.W.W., Marshall, J. A., Houck, J. R. et al. 2007, ApJL 654, L49
60. Spiesman, W.J., Hauser, M. G.; Kelsall, T. et al. 1995, ApJ 442, 662.
61. Stecker, F.W. 1969, ApJ 157, 507
62. Stecker, F. W., De Jager, O.C. & Salamon, M. H. 1992, ApJ 390, L49
63. Stecker, F.W. & Glashow, S.L. 2001, Astropart. Phys. 16, 97
64. Stecker, F.W., Malkan, M. A. & Scully, S. T. 2006 (SMS06), ApJ 648, 774
65. Stecker, F.W., Malkan, M. A. & Scully, S. T. 2012 (SMS12), ApJ 761, 128
66. Vaccari, M., Marchetti, L., Franceschini, A. et al. 2010, A&A 518, L20
67. Viero, M.P., Moncelsi, L., Quadri, R. F. et al. 2015, ApJ 809, L22
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters