Measurement of the Extragalactic Background Light Spectral Energy Distribution with VERITAS

The universe is permeated by diffuse photon fields, tracking light emission from the epoch of recombination to the present day. At ultraviolet (UV), optical and infrared wavelengths, the extragalactic background light (EBL) is primarily due to stellar radiation in the UV/optical band, and reradiation of starlight into the infrared band after absorption by matter (primarily dust). Secondary contributions to the UV/optical intensity come from energy release associated with accretion of matter in active galactic nuclei. More exotic contributions have also been hypothesized, including dark-matter decay in the early universe and radiation from primordial stars (Dwek et al. 2005). A number of comprehensive reviews discuss the EBL (Hauser & Dwek 2001; Dwek & Krennrich 2013); only key points are discussed here.

The EBL intensity as a function of wavelength and its evolution with redshift have been predicted by a number of authors (Finke et al. 2010; Domínguez et al. 2011; Gilmore et al. 2012; Franceschini & Rodighiero 2017, 2018). While a detailed description of these models is beyond the scope of this paper, these approaches require information about stellar and galactic formation and evolution, and models of the matter (dust) distribution in galaxies and of radiative transfer. Comparing the measured spectral energy distribution (SED) of the EBL against theoretical predictions thus tests the current understanding of the density of and radiation by visible matter since the epoch of reionization.

While direct measurements of the EBL SED have been attempted, most recently by Matsuura et al. (2017) and Zemcov et al. (2017), correct subtraction of contaminating foregrounds is challenging, and the measurements tend to be limited by systematic uncertainties. Robust lower bounds on the EBL SED, on the other hand, have been generated by summing the light expected based on resolved galaxies. Deep surveys in recent years have improved the precision of these measurements (Madau & Pozzetti 2000; Fazio et al. 2004; Keenan et al. 2010; Driver et al. 2016). It is worth noting that by design, most EBL models are in good agreement with the limits from galaxy counts measurements. However, this approach is by construction insensitive to contributions from unresolved sources and diffuse emission.

Spectral measurements of distant gamma-ray emitters provide an alternative method to measure the EBL SED, with an independent set of uncertainties and sensitivity to diffuse contributions. The interaction of very-high-energy (VHE; >100 GeV) gamma-rays with EBL photons has been discussed since the early days of VHE gamma-ray astronomy (Nikishov 1962; Gould and Schréder 1967). At center-of-mass energies above the energy threshold for electron-positron pair creation, the VHE gamma-ray flux will be attenuated by interaction with the EBL. The probability of an interaction increases with VHE gamma-ray energy and with the distance the gamma-ray traverses. Consequently, EBL attenuation results in differences between the EBL-absorbed and intrinsic spectra of VHE gamma-ray emitters. For a given energy and source distance, these differences are expected to evolve in the same way in blazar spectra, independent of the source.

High-frequency-peaked blazars (HBLs) are canonically used in measurements of EBL-induced spectral modification, owing to their abundance (∼50 HBLs were detected in the VHE range at the time of writing²⁸ ), their detection at VHE to large redshifts (for example, z = 0.604 for PKS 1424+240), and their emission of gamma-rays up to high energies (detected to tens of TeV for nearby sources). Intermediate-frequency-peaked blazars (IBLs), though less abundant in the VHE catalog, can also be used in EBL measurements. As discussed in Section 4, it is necessary to make some assumptions about the intrinsic source spectra to extract the EBL imprint from the observed spectra. However, biases resulting from such choices are mitigated by making conservative assumptions and utilizing observations of a large sample of HBLs and IBLs located at a range of distances.

VERITAS utilizes the imaging atmospheric Cerenkov technique, measuring Cerenkov emission induced by particle showers in Earth's atmosphere. The imaging of such showers allows the reconstruction of the energy and arrival direction of their gamma-ray progenitors. VERITAS consists of an array of four telescopes (Holder et al. 2006) located at the Fred Lawrence Whipple Observatory in southern Arizona, USA (+31°40'3021, −110°57'777, 1268 m above sea level). Each telescope includes a 12 m reflector of Davies–Cotton design (Davies & Cotton 1957) and a photomultiplier tube camera with 499 pixels. The instrument has been upgraded several times: in 2009, one of the telescopes was relocated, improving sensitivity with a more symmetric array, and in 2012, the cameras and trigger system were upgraded, improving sensitivity at the lowest energies (Kieda et al. 2013). The instrument is sensitive to gamma-rays with energies from ∼85 GeV to greater than 30 TeV, has a field of view of 35, an energy resolution of 15%–25%, and an angular resolution (given as the 68% containment radius) of <01 at 1 TeV (Park et al. 2015).

The data used in this measurement were collected with the camera center offset by 05 from the source position (wobble pointing mode). This observation mode facilitates background estimation, as source and background regions can be defined within the same field of view (Fomin et al. 1994). The data included in this analysis span from the beginning of VERITAS operations in 2007 September to 2016 June. Data were collected under good weather conditions, and include observations taken under both dark sky conditions, and for a fraction of the data (∼20%), under moonlight (Archambault et al. 2017). The increased energy threshold of moonlight observations was considered acceptable, as the spectral behavior at the highest energies is of interest for EBL measurements. The data were processed with one of the standard VERITAS calibration and event reconstruction pipelines (Krause et al. 2016; Maier & Holder 2017), and the results were cross-checked with an independent chain (Krawczynski et al. 2006; Cogan 2007; Daniel et al. 2007). A gamma/hadron selection corresponding to an energy threshold (defined by the average energy bias falling below 10%) of ∼170 GeV (205 GeV for a few cases) was used to preserve sensitivity to low-energy gamma-rays.

The blazars considered in this analysis are listed in Table 1. The source sample comprises 13 HBLs and 1 IBL, 3C 66A. The source sample was selected with the goal of smoothly covering a broad redshift range, with a focus on objects with high detection significance to minimize statistical uncertainties. As discussed in Section 4, the photon spectrum is the primary input to the EBL analysis, and the uncertainties on the spectral points determine the uncertainty of the EBL measurement.

A distinguishing characteristic of blazars is their flux variability, which in some cases has been shown to correlate with spectral variability (Albert et al. 2007; Fossati et al. 2008; Abramowski et al. 2010). Averaging together hard and soft spectral shapes can distort the shape of the time-averaged spectrum, leading to the appearance of spectral hardening at high energy if a spectrally hard bright state and a spectrally soft dim state are averaged. To avoid this and the consequent risk of biases, the flux and spectral variability of the sources were studied. For sources that showed evidence for flux variability (as shown in Table 1, all sources other than 1ES 0414+009), the data sets were divided by flux level and spectral fits were performed for each flux bin. The spectra for sources that showed possible spectral variability (>2σ variation in a spectral parameter other than the flux normalization) are divided into high- and low-state spectra and treated independently in the EBL analysis. As shown in Table 1, only 1ES 1959+650 and 3C 66A show evidence of spectral variability. For 1ES 0502+675 and MS 1221.8+2452, the sources experienced bright states outside of which the sources were not detected: for the former over a single observing period, and for the latter during a two night flare. For these two sources, spectra were constructed only from data collected during these bright periods.

Three spectral models are considered: a simple power law (${dN}/{dE}={N}_{0}{(E/{E}_{0})}^{-{\rm{\Gamma }}}$), a power law with an exponential cutoff (${dN}/{dE}={N}_{0}{(E/{E}_{0})}^{-{\rm{\Gamma }}}\exp (-E/{E}_{{\rm{c}}})$), and a log parabola (${dN}/{dE}={N}_{0}{(E/{E}_{0})}^{-a-b\mathrm{log}(E/{E}_{0})}$). A binned χ² fit is utilized. The majority of the spectra are well-described by a simple power law, in which case a more complex fit function is not considered. Table 2 summarizes the fit results for all observed spectra that are well described by a power law, while Tables 3 and 4 show the fit results for observed spectra best described by a power law with an exponential cutoff and a log parabola, respectively. For several cases, both a power law with an exponential cutoff and a log parabola provided an acceptable fit, in which case the better fit is shown.

The approach taken here for extracting the EBL significance has been employed successfully in the literature (Mazin & Raue 2007; Orr et al. 2011; Biteau & Williams 2015). Rather than evaluating the consistency of the observed blazar spectrum with a model (such as Franceschini et al. 2008; Finke et al. 2010; Domínguez et al. 2011; Gilmore et al. 2012), the gamma-ray extinction predicted by generic EBL shapes and its impact on the blazar spectra described in Section 3 are studied.

An example of several EBL shapes is shown in Figure 1. The shapes are generated by defining a set of 12 grid points in the wavelength of EBL photons, λ_EBL, within the range 0.1–100 μm. The distance between grid points is set to be of the order of the FWHM of the EBL kernel (the pair-production cross section integrated along the line of sight). The EBL kernel width imposes a physical limit on the narrowest structure that can be probed. At each grid point, a value for the EBL intensity is set via random number generation between 1.0 and 50.0 nW m⁻² sr⁻¹, producing a flat initial distribution in intensity. The intensity range was chosen to fully contain the upper and lower limits from complementary measurements. The 12 intensity values (one random value for each grid point in λ_EBL) define an EBL shape. The EBL intensity was not allowed to change by more than a factor of 2.5 between consecutive grid points, eliminating EBL shapes with unresolvably sharp structures. This restriction spoils the initially flat EBL intensity probability distribution. The flat starting distribution is regained by attaching a weight at each grid point for each EBL shape. This ensures that the results are not biased toward any particular EBL intensity due to a concentration of EBL shapes at that intensity. It should also be emphasized that external constraints from direct measurements or galaxy counts estimates are not considered in the generation of the EBL shapes. An ensemble of 480,000 EBL shapes is considered in the analysis described below, smoothly covering the intensity space from 1.0 to 50.0 nW m⁻² sr⁻¹.

Zoom In Zoom Out Reset image size

The generated EBL shapes are taken to represent possible EBL SEDs at redshift z = 0. They are represented as second-order splines. As has been discussed in detail in the literature (e.g., Abdalla et al. 2017; Dwek & Krennrich 2013; Biteau & Williams 2015), an EBL number density at z = 0 can be translated into an optical depth τ(E, z), quantifying the attenuation of a gamma-ray flux of energy E traveling a distance z. For each EBL shape i, τ_i(E, z) is numerically computed for a range of energies spanning 0.1–20.0 TeV and redshifts z = 0.03–1.0. The energy and redshift points at which the calculations are made are spaced evenly in logarithmic space in energy, and semi-logarithmically in redshift. Flat ΛCDM cosmology is assumed in the calculation, with dark-energy density Ω_Λ = 0.73, matter density Ω_M = 0.27, and Hubble constant H₀ = 70 km s⁻¹ Mpc⁻¹. Fractional numerical integration errors are stored as σ_i(τ_i(E, z)). Integration errors increase with energy and redshift, but on average are a few percent.

The evolution of the EBL with redshift must be considered; both the universe's expansion and the evolution of radiation sources affect the intensity of the EBL. The evolution is handled empirically, as proposed by Madau & Phinney (1998): the EBL number density η_EBL is scaled with redshift as ${(1+z)}^{3-{f}_{\mathrm{evo}}}$ rather than by the purely volumetric (1 + z)³. Multiple values of f_evo were tested, and a value of f_evo approximating the redshift evolution following the models of Domínguez et al. (2011), Franceschini et al. (2008), and Gilmore et al. (2012) was selected. The selected value is 1.7. This value does not accurately reproduce the evolution to all redshifts or account for the differences in evolution between the models. The impact on the measurement is discussed in Section 5. It should be noted that selecting f_evo to reproduce the model-based redshift evolution of the EBL introduces a model dependence into what is otherwise a model-independent measurement.

Following the computation of τ_i(E, z) for the 480,000 EBL shapes, the observed photon spectra are corrected for the EBL absorption predicted by each shape, by multiplying each differential flux point by ${e}^{\tau (E^{\prime} ,z^{\prime} )}$, where E' and z' reflect the energy of the differential flux point and the redshift of the source. A Delaunay interpolation is performed to extract the optical depth at the correct energy and redshift, given that τ_i(E, z) is calculated on a grid in energy and redshift. The numerical uncertainty σ_i(τ_i(E, z)) is propagated into the total uncertainty on the differential flux. The resulting EBL-corrected spectra, dN_i/dE_i, are subjected to a binned χ² fit. As with the observed spectra, three fit models are considered: a power law, a power law with an exponential cutoff, and a log parabola. In the case in which the observed spectrum is well-described by a power law, only a power law is tested for the EBL-corrected spectra, in order to avoid overfitting the data. For EBL-corrected shapes fit with three models, the power-law fit is selected if it results in an acceptable p-value (p ≥ 0.05). If a curved spectral shape is required to obtain an acceptable fit, the model with the better p-value is selected.

Two assumptions are made about intrinsic blazar spectra, and hence on the possible shapes of the EBL-corrected spectra. The first is that the spectra must be power law or concave: the allowed range on the fit parameters describing curvature is restricted, disallowing convex shapes (Dwek & Krennrich 2013). Second, the spectral indices of the EBL-corrected spectra are restricted to be softer (larger) than Γ = 1.0 for the power law and power law with exponential cutoff fits. This is a more conservative restriction than the common assumption that the intrinsic spectral index should correspond to that measured by Fermi-LAT in the GeV energy range, where EBL attenuation is minimal for the redshifts considered here. The choice of Γ = 1.0 is somewhat arbitrary. It is significantly harder than the average spectral indices measured by Fermi-LAT (Ackermann et al. 2015). However, the second Fermi-LAT flare catalog (2FAV) records a small number of flares with Γ ∼ 1.5, with large uncertainties (Abdollahi et al. 2017), motivating the selection of Γ = 1.0.

The EBL shapes are generated with λ_EBL between 0.1 and 100 μm, however, the range of λ_EBL depends on the source redshift and the energy range covered by the photon spectrum. The sensitive wavelength range is defined by the width of the pair-production cross section after integration over the line of sight, and can be approximated as

$$\begin{eqnarray}&&{\lambda }_{\mathrm{EBL}}\simeq [0.5\,\mu {\rm{m}}\mbox{--}5\,\mu {\rm{m}}]\times \left(\displaystyle \frac{{E}_{\gamma }}{1\,\mathrm{TeV}}\right)\times {(1+z)}^{2},\end{eqnarray} \tag{ 1 }$$

taking E_γ as the energy of the lowest-energy and highest-energy differential flux points in the source's photon spectrum. The highest-energy flux point is taken to be either one flux point beyond the last point with ≥2σ significance, or the last point with the number of source region counts >0 by 2σ, whichever point is of higher energy.

A projection of the EBL intensity probability distribution at each λ_EBL grid point is made. Initially a flat distribution, each EBL shape is weighted by exp(−χ²/2), with χ² taken from the fit to the EBL-corrected spectrum. This choice of weighting disfavors EBL shapes whose EBL-corrected spectra are poorly described by the physically motivated spectral models considered.

The constraints from multiple blazars are combined by weighting the EBL model by $\exp {(-({\rm{\Sigma }}{\chi }^{2})}_{n}$/2), where n is an index over the sources. It is worth noting that the spectra are produced with comparable energy binning, ensuring that sources do not receive a greater or lesser weight in the combination based on the choice of binning. At each grid point in λ_EBL, only sources contributing according to Equation (1) are included in the sum. Two example distributions, after χ² weighting, are shown in Figure 2.

Zoom In Zoom Out Reset image size

A measurement of the EBL intensity is made by integrating the desired (68% or 95%) containment. As the probability distributions are in general asymmetric, the quantiles are extracted by integrating inward from where the distribution tails to zero, rather than outward from a central value. Examples of 68% containment bands are shown by the shaded regions in Figure 2. In the right panel, only an upper limit on the EBL intensity can be set, as opposed to the left panel, where both the upper and lower tails of the probability distribution fall to zero.

The accuracy of the method is tested by taking spectra from a bright, nearby source (ten spectra derived from short observations of the Crab Nebula, where event counts are comparable to those of the EBL sources) and calculating the expected spectra at z = 0.09, 0.2, 0.3, 0.4, 0.5, and 0.6 using the model of Gilmore et al. (2012). The extracted 68% containment band contains the input value of the EBL intensity for all redshifts.

Several systematic uncertainties must be considered in the measurement of the EBL SED with the technique described above:

• 1.  
  the number of EBL shapes considered;
• 2.  
  the energy-scale uncertainty and finite energy resolution of VERITAS;
• 3.  
  uncertainties in the redshift evolution of the EBL;
• 4.  
  the inclusion of sources with uncertain redshifts in the analysis (namely PG 1553+113 and 3C 66A).

The choice of 480,000 as the number of generic EBL shapes was made to balance smooth coverage with affordable computing time. Randomly removing 10% of shapes affected the extracted EBL bounds by less than 1%, indicating that the precision of measurement is not limited by the number of shapes.

The uncertainty on the absolute energy scale of VERITAS is ∼20% (Madhavan 2013), while the energy resolution is finite and varies with energy. The combination of these two uncertainties propagates into an uncertainty on the photon index Γ, which becomes relevant when the EBL-corrected spectra are fit with power laws (or power laws with exponential cutoffs) with a limit of Γ = 1.0. A Monte Carlo study was performed per source using the corresponding event counts, effective areas, and exposure times, and assuming a symmetric ∼20% energy-scale uncertainty and an energy-dependent energy resolution. The study verified that the energy-scale/resolution uncertainties contribute a ∼10% uncertainty on Γ. By shifting the limiting Γ by ±0.10 × Γ_obs and recomputing the containment bands on the EBL intensity, the uncertainty is conservatively accounted for. The maximum impact on the measured EBL intensity is 10%.

As discussed in Section 4, for the redshift range considered in this analysis, the EBL evolution of Gilmore et al. (2012) is reproduced within uncertainties by using a scaling factor f_evo = 1.7 to treat the redshift evolution of the EBL number density. However, the evolution of the EBL with redshift is not well known, and different models predict different evolution. Figure 3 shows the ratio of τ predicted by two EBL models, Franceschini et al. (2008) and Gilmore et al. (2012), as a function of energy for different redshifts. The difference in the predicted values of τ between the models grows with redshift and energy. The theoretical uncertainty from the spread in the model predictions can be propagated as an additional uncertainty to the EBL SED measurement. The curves shown in Figure 3 are used to obtain an uncertainty as a function of energy by taking the difference in the curve for z = 0.01 and the curve for the redshift of the source under consideration. The τ_i(E, z) associated with each spectral point of the source's photon spectrum is then varied by the uncertainty, and the most conservative containment band on νI_ν is taken. While this does not perfectly account for the spread in the theoretical models—only two models are considered—it attempts to quantify the impact of an uncertainty that has typically been ignored in the literature. The maximum impact on the extracted EBL intensity is 12%, comparable in size to the impact of the energy-scale uncertainty.

Zoom In Zoom Out Reset image size

Two of the sources included in the measurement, PG 1553+113 and 3C 66A, do not have well-measured redshifts. The redshift ranges for these sources are shown in Table 1. The redshift uncertainty is propagated to the EBL SED measurement as follows: the analysis is performed with the redshifts set to the upper and low bounds of the ranges, and the most conservative containment band on νI_ν is taken.

To summarize, the energy-scale uncertainty and the assumed redshift evolution of the EBL are the largest contributors to the total systematic uncertainty. Section 6 presents the results with and without systematic uncertainties due to the energy scale and redshift evolution. The uncertain redshifts of PG 1553+113 and 3C 66A, however, are accounted for in all presented results.

The measured EBL SED is shown in Figure 4. The 68% and 95% containment bands are shown, with the dark and light shaded regions denoting the bands with and without systematic uncertainties, respectively. The filled circles indicate the median values of the EBL intensity probability distributions. Upward-facing arrows indicate the lower limits from galaxy counts, downward-facing arrows show the upper limits from direct measurements (summarized in Dwek & Krennrich 2013). The results are further summarized numerically in Table 5, including systematic uncertainties. The measured band is in good agreement with estimates from galaxy counts, indicating no hint of a diffuse component, although the uncertainties are relatively large. The VERITAS results are consistently lower than the limits from direct measurements.

Zoom In Zoom Out Reset image size

Table 5.  Measured 68% and 95% Containment Bands on the EBL Intensity as a Function of EBL Wavelength, including Systematic Uncertainties

λ	νIνmin(95% CI)	$\nu {I}_{\nu }^{\min }$(68% CI)	$\nu {I}_{\nu }^{\max }$(68% CI)	$\nu {I}_{\nu }^{\max }$(95% CI)
(μm)	(nW m−2 sr−1)	(nW m−2 sr−1)	(nW m−2 sr−1)	(nW m−2 sr−1)
0.56	⋯	⋯	15.2	25.2
1.00	1.5	2.5	9.4	15.4
1.80	2.5	4.4	12.5	15.9
3.20	2.6	3.5	7.5	9.8
5.60	1.4	1.9	4.7	6.4
10.00	⋯	⋯	2.2	3.5
18.00	⋯	⋯	2.2	3.6
32.00	⋯	⋯	4.7	6.8
56.00	⋯	⋯	5.6	10.4

Note. A dash indicates that no lower limit could be set. The treatment of the redshift uncertainty for two of the sources is discussed in Section 5.

Download table as:  ASCIITypeset image

A test of the importance of individual sources was made, shown in Figure 5. Each of the sources listed in Table 1 was removed from the combined result, and the containment bands recalculated. The most conservative values of $\nu {I}_{\nu }^{\min }$ and $\nu {I}_{\nu }^{\max }$ at each value of λ are shown. The resulting 68% containment band is shown is in Figure 5, as is the 68% containment band considering all spectra (without systematic uncertainties). While the precision of the measurement is reduced when the full source sample is not used, the results are not dominated by any individual source. Indeed, at different wavelengths, the source whose removal most impacts $\nu {I}_{\nu }^{\min }$ and $\nu {I}_{\nu }^{\max }$ varies. The exclusion of the 1ES 1959+650 high state, 1ES 1215+303, 1ES 0229+200, 1ES 1218+304, 1ES 1011+496, and 1ES 0502+675 spectra contribute to the broadening of the containment bands.

Zoom In Zoom Out Reset image size

The measurement shown here is in good agreement with other recent measurements. Figure 6 shows a comparison between the VERITAS 68% containment band and the results of Biteau & Williams (2015) (gamma-ray only, 1σ error bars), Abdalla et al. (2017; 1σ error bars, including systematic uncertainties), and Acciari et al. (2019; 1σ error bars, including systematic uncertainties). The model of Gilmore et al. (2012) is included for reference. All gamma-ray measurements are consistent with each other and with Gilmore et al. (2012), which for this wavelength region is consistent with other state-of-the-art models. Taken together, these measurements point to an EBL SED, that is well described by the lower limits from galaxy counts, although the uncertainties remain too large to rule out a diffuse component. It is clearly seen that the approach taken here has more power to constrain high EBL intensities than low EBL intensities, which has been remarked upon in other EBL measurements using imaging atmospheric Cerenkov telescopes (Abdalla et al. 2017; Acciari et al. 2019).

Zoom In Zoom Out Reset image size

Nine years of VERITAS observations of 14 hard-spectrum blazars, located from redshift z = 0.044 to z = 0.604, were used to extract a measurement of the EBL intensity as a function of wavelength. The approach taken made minimal assumptions about the EBL shape and evolution, and made conservative assumptions about the intrinsic spectra of blazars. No individual source is dominant in the resulting joint measurement, and the systematic uncertainties are small in comparison to the statistical uncertainties. This indicates that the VERITAS results can be improved by deeper exposures and better-measured spectra. This is particularly true at large values of λ_EBL, where the state-of-the-art EBL models begin to disagree in their predictions. The agreement between the extracted EBL SED and lower limits from galaxy counts indicates that most, if not all, of the EBL photons can be attributed to resolved galaxies in the wavelength range covered by this measurement. Looking beyond currently operating instruments, a more complete picture of the EBL SED, particularly at the longest wavelengths, will be provided by the forthcoming Cerenkov Telescope Array (Acharya et al. 2017).

This research is supported by grants from the U.S. Department of Energy Office of Science, the U.S. National Science Foundation, and the Smithsonian Institution, by NSERC in Canada. This research used resources provided by the Open Science Grid, which is supported by the National Science Foundation and the U.S. Department of Energy's Office of Science, and resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231. E.P. acknowledges a Marie Curie Intra-European Fellowship within the 7th European Community Framework Programme and the Young Investigators Program of the Helmholtz Association for support during the period of the project, the DJEI/DES/SFI/HEA Irish Centre for High-End Computing (ICHEC) for the provision of computational facilities and support, and thanks Jonathan Biteau for useful discussion. We acknowledge the excellent work of the technical support staff at the Fred Lawrence Whipple Observatory and at the collaborating institutions in the construction and operation of the instrument.

Facility: VERITAS. -