DETECTION OF THE COSMIC γ-RAY HORIZON FROM MULTIWAVELENGTH OBSERVATIONS OF BLAZARS

Very high energy (VHE; 30 GeV–30 TeV) photons do not travel unimpeded through cosmological distances in the universe. A flux attenuation is expected due to photon–photon pair production with the lower energy photons of the extragalactic background light (EBL) in the ultraviolet, optical, and infrared (IR; e.g., Nikishov 1962; Gould & Schréder 1967; Stecker et al. 1992; Salamon & Stecker 1998). The EBL is the radiation emitted by star formation processes (star light and star light absorbed/re-emitted by dust) plus a small contribution from active galactic nuclei (AGNs) integrated over redshift over all the cosmic star formation history (e.g., Hauser & Dwek 2001; Dwek & Krennrich 2013). Due to the properties of the interaction, a VHE photon of a given energy interacts mainly with an EBL photon within a well-defined and narrow wavelength range. Therefore, a signature of the EBL spectral distribution is expected in the observed VHE spectra of extragalactic sources (Ackermann et al. 2012b; Abramowski et al. 2013; Yuan et al. 2013).

An interesting feature in the observed VHE spectra of extragalactic sources as a consequence of EBL absorption is given by the cosmic γ-ray horizon (CGRH), which has not been clearly observed yet. The CGRH is by definition the energy, E₀, at which the optical depth of the photon–photon pair production becomes unity as a function of redshift. Therefore, it gives an estimate of how far VHE photons can travel through the universe. Due to the exponential behavior of the flux attenuation, an alternative definition is the energy at which the intrinsic spectrum is attenuated by the EBL by a factor of 1/e (see, e.g., Aharonian 2004). (The intrinsic source spectrum is the one that we would observe if there were no effect from the EBL, also known as the EBL-corrected spectrum.)

An extreme category of AGNs, known as blazars, has been shown to be the best target for extragalactic VHE detections. They are characterized by having their energetic γ-ray jets pointing toward us. In fact, most of the extragalactic sources already detected by imaging atmospheric Cherenkov telescopes (IACTs) such as H.E.S.S., MAGIC, and VERITAS (Hinton 2004; Lorenz 2004; Weekes et al. 2002, respectively) are blazars.¹⁰ The observations of broadband spectral energy distributions (SEDs) of blazars show that they are characterized by a double-peaked shape and that their emission covers all the electromagnetic spectrum from radio up to the most energetic γ-rays. The synchrotron/synchrotron self-Compton (SSC) model provides a successful explanation for this behavior for most cases (e.g., Abdo et al. 2011b, 2011d; Zhang et al. 2012). In this framework, a population of ultrarelativistic electrons causes the lower energy peak by synchrotron emission, while the second peak is then accounted for by inverse Compton production of γ-rays from the same population of high-energy electrons and photons in the low-energy peak.

The direct observation of the CGRH in the VHE spectra of blazars remains elusive. This is due to two main observational difficulties. First, the lack of knowledge of the intrinsic spectra at VHE. Extensive multiwavelength campaigns from radio up to γ-rays are needed in order to predict the unattenuated VHE emission from the synchrotron/SSC model with enough precision. These campaigns should preferably be simultaneous due to the short-time flux variability of blazars (e.g., Aleksić et al. 2011a, 2011b). In this situation, the typical procedure in the literature to estimate the intrinsic VHE spectrum is either to assume a limit for the hardness of the slope E^−Γ with Γ = 1.5 (Aharonian et al. 2006) or to extrapolate the Fermi-Large Area Telescope (LAT) spectrum up to higher energies (Georganopoulos et al. 2010; Orr et al. 2011).

Second, it has been possible only in recent years to detect a considerable number of blazars in the GeV energy range to allow us a statistical analysis. A large sample is necessary to reject intrinsic behaviors in the sources that mimic the effect of the CGRH. This improvement has been made thanks to the large data sets provided by the Fermi satellite (Ackermann et al. 2011) and the IACTs.

The only estimations of the CGRH so far are EBL-model-dependent lower limits from VHE observations of blazars (Albert et al. 2008), lower limits from Fermi-LAT observations of blazars (Abdo et al. 2010b), and the predictions from EBL models (e.g., Franceschini et al. 2008; Kneiske & Dole 2010; Finke et al. 2010; Domínguez et al. 2011a, hereafter D11; Gilmore et al. 2012; Stecker et al. 2012). (A table with a classification and description of the ingredients of these EBL models can be found in the proceeding by Domínguez 2012.) Indeed, an independent observation of the CGRH will also provide a completely independent and new test to the modeling of the EBL and consequently constraints on galaxy evolution. Furthermore, the CGRH measurement also can be useful to estimate the cosmological parameters with a novel and independent methodology (Blanch & Martinez 2005a, 2005b, 2005c; Domínguez & Prada 2013). The detection of the CGRH is a primary scientific goal of the Fermi Gamma-ray Telescope (Hartmann 2007).

This paper is organized as follows. Section 2 describes the blazar catalog used in our analysis. The methodology is explained in Section 3. Section 4 shows the results obtained from our analysis and in Section 5 the results are discussed. Finally, a brief summary of the main results is presented in Section 6.

Throughout this paper a standard ΛCDM cosmology is assumed, with Ω_m = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹ (Larson et al. 2011; Komatsu et al. 2011).

A catalog of quasi-simultaneous multiwavelength data from radio up to VHE for 15 blazars has been built. The data for energies lower than the Fermi-LAT regime (20 MeV–>300 GeV; Atwood et al. 2009) are taken from the data compilation presented by Zhang et al. (2012). We refer the reader to that paper for references and details of the different data sets and their simultaneity. That catalog is combined here with Fermi-LAT data from the second Fermi-LAT AGN catalog (2FGL; Ackermann et al. 2011) and the new Fermi-LAT hard-spectrum catalog (1FHL; M. Ackermann et al., in preparation), which include 2 and 3 yr of Fermi-LAT data, respectively. The 2FGL contains fluxes in the following six energy bins: 30–100 MeV, 100–300 MeV, 300 MeV–1 GeV, 1–3 GeV, 3–10 GeV, and 10–100 GeV. The 1FHL contains fluxes for the sources with the highest-energy emission in three bins that reach higher energy than the ones used in the second-year catalog (10–30 GeV, 30–100 GeV, and 100–500 GeV). Our final catalog includes data from IACTs as well. We use simultaneous VHE data from IACTs when available; otherwise the IACT observation closest in time to the lower-energy data is used following the suggestions by Zhang et al. (2012). Table 1 lists the 15 sources in the set of blazars that we study here. The catalog presented in Zhang et al. (2012) contains 24 blazars; however, we could not use all of them due to non-detections either by Fermi or the IACTs, which are essential for applying our methodology (see Section 3).

Our catalog covers a wide redshift range from z = 0.031 to z ∼ 0.5 and all 15 blazars are classified as BL Lac sources (which are typically characterized by rapid and large-amplitude flux variability and significant optical polarization). We note that in the cases of PG 1553+113 and 3C 66A the redshifts are uncertain, but still these sources are included in the analysis. A redshift in the range 0.395 < z ⩽ 0.58 is estimated for PG 1553+113 by Danforth et al. (2010; the upper limit is 1σ). The blazar 3C 66A typically is cited as having a redshift of z = 0.444, which is used here as well (cf. Bramel et al. 2005; Finke et al. 2008b).

Our methodology consists of finding the best-fitting homogeneous synchrotron/SSC models from multiwavelength data as simultaneous as possible from radio to the highest-energy γ-rays detected by the Fermi-LAT for the blazars in our catalog. These models predict the unattenuated VHE fluxes, which are compared with detections by IACTs. The ratios between the unattenuated and detected VHE fluxes give an estimate of the EBL optical depth. By means of a maximum likelihood technique that is independent of any EBL model and that is based only on a few physically motivated assumptions, we then derive the CGRH for each blazar.

3.1. Broadband Spectral-energy-distribution Fitting and Optical-depth Data Estimation

For every blazar in our sample, we built a quasi-simultaneous SED based on the data collected by Zhang et al. (2012) and on the LAT data from the 2FGL (Ackermann et al. 2011) and the 1FHL (M. Ackermann et al., in preparation). These data are shown in the insets of Figure 1. In many cases, these SEDs are not constructed from strictly simultaneous data. However, since we preferentially choose SEDs that are for a low state, we expect that the effects of variability are minimal in the γ-ray energy range. We then fit the unique multiwavelength data of each source with a synchrotron/SSC model using a χ² minimization technique. The model and fitting technique are fully described in Finke et al. (2008a; see also Mankuzhiyil et al. 2010). This model includes photoabsorption by photons internal to the blob. Since BL Lac objects and their presumed misaligned counterparts, FR I radio galaxies, usually lack optically thick dust tori (Donato et al. 2004; Plotkin et al. 2012), we do not include photoabsorption from this radiation source, although it could in principle be important. We note that for 1ES 1959+650, there is evidence for an optically thick dust component (Falomo et al. 2000). However, we did not obtain positive results with this source, as discussed in Section 5 below. The fitting technique is double nested, with the inner loop fitting the synchrotron component with a particular electron distribution. In this paper, we use a broken power law for the electron distribution with an exponential cutoff at high energies, that is, at γ_max. In some cases (Mrk 421, 1ES 2344+514, PKS 2005−489, H2356−309, 1ES 218+304, and 1ES 1101−232), we found that a single power law with exponential cutoff was sufficient to provide good fits. Our fits have as free parameters the electron indices p₁, p₂; the minimum, maximum, and break electron Lorentz factors, γ_min, γ_max, and γ_brk, respectively; and the overall normalization. Often γ_min and γ_max were kept constant during the fit. The outer loop fits the SSC model, the high-energy data, and has three free parameters: the Doppler factor, δ_D; the magnetic field strength, B; and the minimum variability timescale, t_{v, min}. We assume that the bulk Lorentz factor Γ_bulk = δ_D. In the fits, t_{v, min} was kept constant, with only δ_D and B as free parameters. We discuss in more detail our choices of t_{v, min} below. In the fits, we specifically leave off the IACT data; we only fit the IR through the LAT γ-rays. We use any radio points as upper limits, since this emission is likely from another region of the jet. We leave out the IACT data because for this fit, we are fitting data that are unaffected by EBL attenuation. In some of the more distant sources, the highest energy LAT point (at energy ∼224 GeV) can suffer significant attenuation. Therefore, we remove this data point for sources at z > 0.05. This choice is supported by a variety of observational evidences (see the proceedings by Domínguez 2012 and Primack et al. 2011). Once we have the resulting model curve from our fit that includes the extrapolation to VHE energies (the overall model is shown with a black line in the insets of Figure 1), we use this as the unattenuated/intrinsic spectrum for the source. We compare it with the flux observed from the IACT detection to calculate the absorption optical depth,

$$\begin{equation} \tau (E,z)=\ln (\frac{dF}{dE} {|}_{{\rm int}} / \frac{dF}{dE} {|}_{{\rm obs}}) \end{equation} \tag{ 1 }$$

for photons observed in an energy bin centered at energy E and a source at redshift z. Here F_obs is the observed differential flux and F_int is the intrinsic flux at the energies given by the IACT detection, i.e., the fluxes given by the synchrotron/SSC model evaluated at the energies sampled by the IACT. The uncertainties in τ come directly from the uncertainties in the IACT observations. The log₁₀(τ) data are shown with blue crosses in Figure 1. The method used here for measuring τ(E, z) is similar to the one described by Mankuzhiyil et al. (2010).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the synchrotron/SSC model, the radius of the spherical emitting region, R_blob, is determined from the minimum variability timescale, t_{v, min}, which is in turn constrained by the observed variability timescale through light travel time arguments so that t_{v, min} ⩽ t_v (e.g., Finke et al. 2008a). For consistency, we used the same t_{v, min} for all of our blazars, t_{v, min} = 10⁴ s and t_{v, min} = 10⁵ s. As we see in Section 4, the choice of variability time makes very little difference to the model curve, although it has a large effect on the model fit parameters that are not the focus of this paper. Thus, we are confident that the choice of t_{v, min} has little effect on our resulting measurement of τ(E, z). However, although the two SSC models are similar, they predict different VHE fluxes, which allows us to include the uncertainty in the variability timescale in our analysis.

It should be noted that several sources have been observed to have extremely rapid (∼10² s) variability timescales (e.g., Aharonian et al. 2007). These rapid flares are quite rare, and since we have chosen SEDs from blazars in a quiescent state, we think these short timescales (and the small emitting regions they imply) are unlikely. Furthermore, we have fit several of our objects with t_{v, min} = 10² s and found the resulting Doppler factors to be extremely high, δ ≳ a few hundred, and in two cases (1ES 1101−232 and 3C 66A) as high as 1000. We believe these high Doppler factors to be unreasonably high.

3.2. Maximum Likelihood Polynomial Fitting

We assume that log₁₀(τ) (as obtained from Equation (1)) may be described by a third-order polynomial in log₁₀(E),

$$\begin{equation} \log _{10}(\tau) = a_{0}+a_{1}\log _{10}(E)+a_{2}\log _{10}^{2}(E)+a_{3} \log _{10}^{3}(E), \end{equation} \tag{ 2 }$$

where E is in units of TeV. Lower order polynomials are not sufficient to describe the optical depth. Higher order polynomials introduce too many degrees of freedom and increase the computational time without increasing the precision of our analysis. This shape of the opacity (which is the integral of the EBL spectral intensity and the pair-production cross section; see, e.g., Domínguez et al. 2011a) is expected in the VHE range for two main reasons. First, the EBL SED must have a smooth shape as a consequence of the galaxy SEDs that produce the EBL and second because of the continuity of the cross section of the pair-production interactions (e.g., Dwek & Krennrich 2005). A maximum likelihood method that scans the parameter space is adopted to compute the likelihood of the estimated data given the different polynomials log₁₀(τ) (blue crosses in Figure 1) for each blazar. The four parameters of the third-order polynomial are explored by studying their probability density distributions. At this point, three physically motivated assumptions are made. First, that τ < 1 at E = 0.03 TeV, since EBL attenuation is expected to be significantly low at these energies. Second, that τ ⩽ UL(E, z), where UL is an upper limit calculated in the present work from the EBL upper limits presented in Mazin & Raue (2007); in particular, 1 ⩽ τ ⩽ UL(z) at E = 30 TeV. In Mazin & Raue (2007), two upper limits are presented coming from a realistic and a weaker assumption on the blazar emission that is called extreme and that provides the higher upper limits. The extreme case is used in our analysis since we want to keep our methodology conservative.

Third, we impose that τ should increase monotonically with energy. These assumptions will be discussed in more detail in Section 5. The E₀ value derived from each blazar, for both slow and fast variability timescales, is then estimated from the most likely polynomial in the four-dimensional parameter space (these E₀ values are given in each inset of Figure 1 and also in Table 2). The uncertainty is estimated by using a standard jackknife analysis (Wall & Jenkins 2003). The final (combined) E₀ for each blazar is then calculated as the geometric mean value for the two variability timescales.¹¹ We stress that the uncertainties in this value include the uncertainties derived from the two SSC modelings (physically the uncertainty in the minimum time variability t_{v, min}), which are bracketed by the two different predictions of the VHE fluxes. The lower and upper uncertainties of the combined E₀ are taken from the E₀ − ΔE₀ and E₀ + ΔE₀ of the state with the lowest and highest E₀, respectively. We stress that these uncertainties are more conservative than 1σ.

Table 2. The Measured and Predicted Cosmic γ-Ray Horizon

Source	Redshift	E0 ± ΔE0(fast/slow)	E0 ± (ΔE0)stat ± (ΔE0)sys	ED11 ± ΔED11	IACT Reference
		(TeV) (a)	(TeV) (b)	(TeV) (c)
Mkn 421	0.031	10.42 ± 7.72/11.91 ± 8.79	$11.14_{-8.44}^{+9.56}\pm 2.23$	$9.72_{-3.17}^{+1.85}$	Abdo et al. (2011d)
Mkn 501	0.034	11.85 ± 16.84/2.28 ± 1.02	$5.20_{-3.94}^{+23.49}\pm 1.04$	$8.75_{-3.31}^{+1.68}$	Acciari et al. (2011)
1ES 2344+514	0.044	None/6.65 ± 21.49	None	$6.01_{-3.23}^{+1.20}$	Albert et al. (2007a)
1ES 1959+650	0.048	None/none	None	$5.12_{-2.99}^{+1.02}$	Tagliaferri et al. (2008)
PKS 2005−489	0.071	1.99 ± 0.26/2.09 ± 0.25	$2.04_{-0.31}^{+0.30}\pm 0.41$	$1.83_{-1.06}^{+0.34}$	Kaufmann et al. (2009)
W Comae	0.102	None/none	None	$0.90_{-0.18}^{+0.09}$	Acciari et al. (2008)
PKS 2155−304	0.116	0.77 ± 0.17/0.88 ± 0.05	$0.82_{-0.22}^{+0.11}\pm 0.16$	$0.77_{-0.13}^{+0.07}$	Aharonian et al. (2009)
H1426+428	0.129	6.23 ± 7.64/13.24 ± 16.81	None	$0.68_{-0.11}^{+0.06}$	Aharonian et al. (2002)
1ES 0806+524	0.138	0.35 ± 0.04/0.85 ± 0.01	$0.55_{-0.24}^{+0.31}\pm 0.11$	$0.64_{-0.10}^{+0.05}$	Acciari et al. (2009a)
H2356-309	0.165	None/none	None	$0.54_{-0.07}^{+0.04}$	Abramowski et al. (2010)
1ES 1218+304	0.182	0.58 ± 0.02/0.46 ± 0.02	$0.52_{-0.08}^{+0.08}\pm 0.10$	$0.49_{-0.06}^{+0.04}$	Acciari et al. (2010)
1ES 1101−232	0.186	0.41 ± 0.02/0.39 ± 0.01	$0.40_{-0.02}^{+0.03}\pm 0.08$	$0.48_{-0.06}^{+0.04}$	Aharonian et al. (2006)
1ES 1011+496	0.212	None/none	None	$0.43_{-0.05}^{+0.03}$	Albert et al. (2007b)
3C 66A	0.444	0.29 ± 0.02/0.31 ± 0.02	$0.30_{-0.03}^{+0.03}\pm 0.06$	$0.23_{-0.02}^{+0.02}$	Acciari et al. (2009b); Aleksić et al. (2011a)
PG 1553+113	0.500$^{+0.080}_{-0.105}$	0.24 ± 0.01/0.23 ± 0.02	$0.23_{-0.03}^{+0.05}\pm 0.05$	$0.21_{-0.02}^{+0.02}$	Aleksić et al. (2010)

Notes. The 15 blazars in our catalog are listed with their estimated redshifts. The energy E₀ is shown for the two different variability timescales (t_{v, min} = 10⁴ s for the fast model and t_{v, min} = 10⁵ s for the slow model) and its combined value as described in the text in Columns (a) and (b), respectively. The combined E₀ is given with its statistical and systematic uncertainties (see the text for details). The E₀ estimated from the EBL model discussed in Domínguez et al. (2011a) (E_D11 ± ΔE_D11) is given as well in Column (c). None means that our methodology outputs no solution for the E₀.

Download table as:  ASCIITypeset image

The parameters that describe the synchrotron/SSC models from the fits to the quasi-simultaneous multiwavelength data are listed in Table 1. We provide two different fits to every blazar bracketing the expected intrinsic VHE fluxes. (These two fits are named slow and fast according to their variability timescale.) The methodology described in the previous section is applied to every blazar in our catalog. As we mentioned in Section 2, PG 1553+113 has an uncertain but well-constrained redshift (Danforth et al. 2010). Therefore, two different fits are provided for both redshift limits for this blazar.

Figure 1 shows all the fits for the 15 blazars used in our analysis (two fits per blazar for each minimum variability timescale, except four fits for PG 1553+113 to account for its redshift uncertainty). The fast minimum time variability fits (t_{v, min} = 10⁴ s) are shown on the left side of the figure whereas the slow minimum time variability fits (t_{v, min} = 10⁵ s) are shown on the right side. Each panel shows the log₁₀(τ) data derived from Equation (1) versus the log₁₀ of the energy in TeV. Figure 1 shows the upper limits of the optical depth calculated from the EBL upper limits provided by Mazin & Raue (2007) and the most likely polynomials. The E₀ value is calculated as the energy where log₁₀(τ) = 0 from the most likely polynomial in the maximum likelihood parameter space distribution. For comparison, the estimation of E₀ calculated from the EBL model based on observations presented by D11 is shown in every panel as a red dotted line. The uncertainties in the optical depth from the EBL modeling are shown as a red area. Every panel in Figure 1 also has an inset with the multiwavelength data and the best-fit synchrotron/SSC model (E²dN/dE versus log₁₀ of the frequency in Hz).

The final E₀ is then assessed combining the two E₀ results (slow/fast minimum time variability) from every blazar as the geometric mean value between these two estimates. The statistical uncertainties bracket the two E₀ values from every fit for every blazar with their uncertainties estimated from the maximum likelihood analysis using the likelihood distributions of the polynomial parameters. In our analysis, we also consider the systematic uncertainties in the LAT measurements. Their effect is estimated by artificially hardening and softening the overall SEDs of PKS 2005−489 (a blazar with low statistical uncertainties). First, the fluxes of the three lowest-energy LAT data for PKS 2005−489 are decreased by 10% (which is the typical Fermi-LAT systematic uncertainty of the effective area; Ackermann et al. 2012a) and the three highest-energy LAT bins are increased in flux by 10%. The overall SED fit is done with these new points, and the energy where τ = 1 is estimated, using the procedure described in Section 3. This procedure is repeated by increasing by 10% the three lowest-energy LAT data of PKS 2005−489 whereas decreasing by 10% the three highest-energy LAT data. This allows us to estimate an average systematic uncertainty of 20% in the energy where τ = 1 for PKS 2005−489. We thus assume that the systematic uncertainty from the uncertainty in the LAT is 20% for all sources. The observed CGRH is shown in Figure 2 with blue circles, the statistical uncertainties are shown with darker blue lines, and the statistical plus systematic uncertainties (added in quadrature) with lighter blue. A completely independent estimation of the CGRH from the EBL model described in D11 is also shown with its uncertainties, which are thoroughly discussed in D11. The uncertainties in the EBL modeling are larger in the far-IR region for the reasons discussed in D11. This leads to the larger uncertainties in the estimation of the CGRH from the EBL modeling at the lower redshifts. The reason is that this is the EBL region that mainly interacts with the higher-energy VHE photons that lead to determination of the CGRH in that redshift range.

Zoom In Zoom Out Reset image size

Our methodology offers more information on the optical depth than just the CGRH. Therefore, the same procedure followed to calculate the CGRH is applied to calculate the energies at which the optical depth is equal to 0.5, 2, and 3 (shown in Figure 3 with blue squares, green triangles, and magenta diamonds, respectively). The energies for those optical depths are plotted from the D11 model with their uncertainties as well (the same colors are used for each modeled optical depth as for the data). The uncertainties in these estimates are significantly larger since, in general, they are given by the extrapolation of the most likely polynomial outside the energy range of the Cherenkov detections.

Zoom In Zoom Out Reset image size

In this work, we present an estimation of the CGRH based on a multiwavelength compilation of blazars that includes the most recent Fermi-LAT data. We stress that our estimation of the CGRH is derived with only a few physically motivated constraints. These results represent a major improvement with respect to previous works. These previous works provide only lower limits for the CGRH such as the EBL-model-dependent limits estimated by Albert et al. (2008; which are based on a modified parameterization of the EBL models presented by Kneiske et al. 2002). Other CGRH limits are presented by Abdo et al. (2010b) using only Fermi-LAT observations.

The 1FHL (M. Ackermann et al., in preparation) is included in our analysis. The inclusion of this data set in our multiwavelength blazar catalog is essential for the right estimation of the CGRH since these measurements help to resolve the shape of the inverse Compton peak.

The optical depth is calculated using Equation (1), which describes the ratio between the intrinsic flux from the synchrotron/SSC models and the observed flux by IACTs. Then, these data are fitted to polynomials of third order imposing some constraints. We also require an increasing and monotonic behavior of the polynomials. Polynomials of order lower than three would not reproduce the expected optical-depth shape while larger order polynomials would introduce unnecessary parameters into the fits. The constraints are all physically motivated and EBL-model independent. As we said before, the first condition is that τ ⩽ 1 at E = 0.03 TeV, which means that the attenuation is rather weak at those low energies.

The second constraint is that 1 ⩽ τ ⩽ UL(z) at E = 30 TeV, where UL are the opacities calculated from the EBL upper limit in the local universe found in Mazin & Raue (2007). The upper limits of their so-called extreme case are used in our analysis. This extreme case represents the least constraining assumption on the blazar spectra since it allows us a wider range of spectral indices (i.e., this results in a rather conservative hypothesis for our analysis). For this same reason, we prefer to use as conservative upper limits the results by Mazin & Raue (2007) rather than the newer results by Meyer et al. (2012) that are based on a more constraining spectral condition. The EBL evolution is expected to affect the optical depth calculated at higher redshifts. To account for this effect, we evolve conservatively the EBL upper limits at all wavelengths as (1  +  z)⁵ (in the comoving frame) when calculating the optical depths from these EBL limits from Mazin & Raue (2007). We note that this is a robust limit given the fact that the maximum evolution (which is dependent on the wavelength) is (1 + z)^2.5 in a realistic model such as D11 for 0 ⩽ z ⩽ 0.6 (the redshift range of our blazar catalog).

The third constraint that we apply for our fits is to require only monotonically increasing functions for log₁₀(τ) as a function of log₁₀(E). This condition is also expected for any realistic EBL spectral intensity, which comes from galaxy emission, given the increasing behavior of the pair-production interaction with energy. Interestingly, we see in Figure 1 that in most cases the IACT observations are indeed detecting the flux decrement given by the CGRH feature (i.e., the Cherenkov observations span from negative to positive values of log₁₀(τ)).

We find that the CGRH derived from 9 out of 11 blazars where our maximum likelihood methodology can be applied is compatible with the expected value from the D11 model. The estimations from other EBL models such as Franceschini et al. (2008), Finke et al. (2010) (model C), and Somerville et al. (2012) are in agreement within uncertainties with the EBL model by D11. We note that the fit of 1ES 1101−232 has only one degree of freedom; see Table 1. The uncertainties of the two lowest redshift blazars (Mkn 501 and Mkn 421) are systematically higher because the optical depth for these cases becomes unity at energies larger than the energies observed by the Cherenkov telescopes. Therefore, in these cases τ = 1 is given by an extrapolation of the polynomials rather than an interpolation between observed energies (see Figure 1), leading to greater uncertainty. For the case of 1ES 2344+514 with fast flux variability timescale, a value of E₀ in agreement with the estimation by the D11 EBL model is derived. However, for this case the uncertainties are larger than E₀ and therefore no useful constraint can be derived. For the case of 1ES 2344+514 with slow flux variability timescale, the SSC-predicted flux is lower than the flux given by IACT data. For H1426+428, both flux variability timescales give uncertainties in the measurement of E₀ larger than E₀ and therefore no constraint can be derived. In both cases, the synchrotron/SSC model does not seem to correctly fit the multiwavelength data. Our maximum likelihood procedure cannot be applied to any flux state on four blazars (1ES 1959+650, W Comae, H2356−309, and 1ES 1011+496). There are different explanations for this fact. Some blazars have shown flux variability on the scale of minutes (e.g., Aharonian et al. 2007; Albert et al. 2008; Aleksić et al. 2011b; Arlen et al. 2013) and the IACTs tend to detect the sources in higher-flux states. In most cases, the LAT data are not simultaneous with the IACT and other multiwavelength data. We have tried to alleviate this problem by choosing SEDs that are based on a low, non-flaring state, where the variability seems to be small. In this way, the effects of variability from epoch to epoch have been minimized. We compare the long-term light curves in X-rays using the quick-look results from the All Sky Monitor (ASM) aboard the Rossi X-Ray Timing Explorer¹² (RXTE) with the time range of the IACT observation for those four blazars where our maximum likelihood procedure could not be applied (see the Appendix for more details). Clearly, 1ES 1011+496 was indeed detected by the IACTs in flaring states. The situation for 1ES 1959+650 is not clear, and the light curve of the H2356−309 observation was rather irregular. We could not find X-ray data for W Comae in the ASM database but this source was clearly detected in TeV in a strong flare (Acciari et al. 2008).

The synchrotron/SSC model is the standard model for fitting high-peaked TeV BL Lac objects and does seem to provide a good fit to their broadband SEDs (e.g., Zhang et al. 2012). However, there are some alternatives. High-peaked BL Lac objects are not thought to have a significant contribution to the γ-ray flux from scattering external photon sources, but there are some exceptions, such as the eponymous BL Lac (Abdo et al. 2011c). It has also been suggested that for some sources, a lepto-hadronic model provides a better fit, such as 1ES 0414+009 (Aliu et al. 2012). Non-variable TeV emission unrelated to the rest of the broadband SED could originate from Compton-scattering of the cosmic microwave background (CMB) by an extended jet (Böttcher et al. 2008), which would certainly complicate their modeling. Another way of creating non-variable TeV emission unrelated to the SED, which would also avoid much of the EBL attenuation, would be if the AGN produces a significant number of cosmic rays, which during propagation interact with the CMB and EBL to produce the observed γ-rays (Essey & Kusenko 2010; Essey et al. 2011). Finally, even if the synchrotron/SSC model is valid for the TeV blazars considered here, it is possible that the electron–positron pairs created by the VHE γ-ray interactions with EBL photons can Compton-scatter the CMB, producing γ-rays observable by the LAT (Neronov & Vovk 2010; Tavecchio et al. 2010; Taylor et al. 2011; Dermer et al. 2011; Vovk et al. 2012). This would complicate the modeling process since it would add other, poorly constrained parameters (Tavecchio et al. 2011). Nonetheless, the simple synchrotron/SSC is a very attractive model due to its success at fitting a large number of objects and its relatively small number of free parameters. The existence of axion-like particles could also allow the γ-rays to avoid the photoabsorption process, changing the expected VHE spectrum (Sánchez-Conde et al. 2009; Domínguez et al. 2011). A detailed broadband modeling including all these non-standard considerations is out of the scope of this work.

In summary, we built a catalog of 15 blazars from Zhang et al. (2012) by requiring a good multiwavelength spectral coverage and TeV detections. After fitting a synchrotron/SSC model to each source, there were four blazars where the TeV detection was at higher fluxes than the flux extrapolation from the models. As described above, the long-term X-ray data from RXTE were checked on the dates of the TeV observations for these blazars, hinting that H2356−309 and 1ES 1011+496 were flaring in X-ray and therefore probably in the TeV range as well. The situation is not clear for 1ES 1959+650. However, as described above potential problems with the synchrotron/SSC model cannot be ruled out. There are two other blazars (H1426+428 and 1ES 2344+514) for which the uncertainties that we derive for E₀ are too high to set any constraint. This is due to large uncertainties in the TeV measurements and/or the low number of TeV spectral points, which do not allow a reliable E₀ estimation.

The agreement between the CGRH observation presented in this work and the expected values from D11 indicates that these possibilities described above might not be relevant for many blazars. Furthermore, Figure 3 gives more information on the optical-depth shape derived from our methodology. This figure shows that the energies at which the optical depths are 0.5, 1 (the standard definition of CGRH), 2, and 3 are still compatible with the D11 model. The uncertainties are generally higher at τ different from 1 due to the fact that the Cherenkov detections do not span the energy range needed in order to derive a better estimation of those energies. As seen in Figure 1, the polynomials cut the horizontal lines of constant optical depths generally in wider energy ranges for τ values different from 1 (i.e., log₁₀(τ) = 0). The agreement between the observed CGRH and the expected CGRH from D11 is also consistent with 3C 66A being located at a redshift slightly lower than z ∼ 0.444¹³ and PG 1553+113 at 0.395 ⩽ z ⩽ 0.58. An independent confirmation of these redshifts will give support to both the current EBL knowledge and our methodology to derive the CGRH. By assuming an EBL model, it is possible to estimate redshifts using EBL attenuation in a realistic way considering the overall SED of the blazars (see the proceeding by Mankuzhiyil et al. 2011). However, we note that some previous estimation of the redshift is necessary in order to fit the synchrotron/SSC models (Abdo et al. 2010a, 2011a).

Orr et al. (2011) claimed that recent EBL models such as D11 are incompatible with IACTs' observations at more than 3σ. They based their conclusions on an analysis of ∼12 blazars using two different methods that they call Method 1 and Method 2. Their more constraining results are derived from Method 2 (see Section 3.2 in Orr et al. 2011). This approach is based on the expected difference between spectral indexes when the VHE spectrum is fitted by a broken power law. This spectral difference is attributed to EBL attenuation, setting limits on the intensity of the local EBL. We consider their results inconclusive. Their Method 2 relies on the assumption that the VHE spectra may be well fitted by broken power laws. We performed F-tests on all the fits of their blazar sample to test whether broken power laws (fitted by two different spectral indexes) could be actually considered better fits to the observed VHE spectra than simple power laws (fitted by a single spectral index).¹⁴ The F-tests performed for every one of their spectra show that for only 2 out of 12 cases (RGB J0152+017 and 1ES 1101-232) the spectra can be considered fitted better by broken power laws than simple power laws. Their results from Method 1 (see Section 3.1 in Orr et al. 2011) are inconclusive as well. This method relies on the assumption that the VHE intrinsic spectrum is described by a power-law extrapolation of the Fermi-LAT data points. As we see in the synchrotron/SSC fits of Figure 1, this is not a realistic assumption due to the shape of the inverse Compton peak. Furthermore, we show in the present work that a more sophisticated SSC-based analysis is compatible with the current EBL knowledge.

Some authors have treated the Fermi-LAT spectrum, extrapolated into the VHE regime, as an upper limit on the intrinsic spectrum and used this to compute upper limits on τ(E, z) (Georganopoulos et al. 2010; Meyer et al. 2012). This provides only upper limits on τ(E, z) rather than measurements, as we derive here. However, their techniques involve fewer assumptions about the blazar emission model and variability of the SED (see the discussion above). Thus, the two techniques are complementary.

Recently, Ackermann et al. (2012b) and Abramowski et al. (2013) claimed the detection of an imprint of the EBL in the blazar spectra. Ackermann et al. (2012b) base their analysis on Fermi-LAT data from z ∼ 0.2 to 1.6, whereas Abramowski et al. (2013) use H.E.S.S. data from blazars located at z ∼ 0.1. These works do not give any results in terms of the CGRH but we are able to estimate it from their results. We find that the results presented in our analysis are compatible with the results from these two independent works, which supports our conclusions.

From our results, we can conclude that the EBL data from direct detection by Cambrésy et al. (2001), Matsumoto et al. (2005), and Bernstein (2007) are likely contaminated by zodiacal light. This possibility has indirectly been proposed previously by several authors such as Aharonian et al. (2006), Mazin & Raue (2007), and Albert et al. (2008) using EBL upper limits, but we confirm these results using a more robust approach.

The CGRH horizon is detected in this work for the first time from a multiwavelength sample of blazars that includes the more recent Fermi-LAT data. Only a few general and physically motivated constraints on the EBL were necessary. As we see from our analysis, the observational estimation of the CGRH is compatible within uncertainties with the derivation by the observational EBL model described by Domínguez et al. (2011a), which is in agreement with the observational EBL model by Franceschini et al. (2008) and the theoretical methodology followed by Somerville et al. (2012) and Gilmore et al. (2012). All these EBL models are realistic representations of the current knowledge of the EBL (Domínguez et al. 2011b; Primack et al. 2011; Domínguez 2012). We have shown the ability of our methodology to study the opacity of the universe at different redshifts and to infer distances of blazars with unknown redshifts. Our methodology is sensitive to the total EBL, which includes light even from the faintest and most distant galaxies in the universe. This will allow us to set limits on the faint-end slope of the evolving galaxy luminosity function, which still remains controversial (see e.g., Reddy & Steidel 2009). The detection of the CGRH presented here will provide an independent test for cosmology and for the estimation of the cosmological parameters that will be presented in Domínguez & Prada (2013).

Our technique will benefit in the future with the improved statistics that Fermi will provide. The future Cherenkov Telescope Array is expected to provide VHE spectra with a better energy resolution, observed up to higher energies, and increase considerably the number of sources, which indeed will improve the CGRH determination. These prospects together with the increasing number of simultaneous multiwavelength observational campaigns (e.g., Abdo et al. 2011a, 2011b, 2011d) are promising for a better estimation of the optical depths due to EBL attenuation using our methodology and for the estimation of the CGRH to z > 0.5.

The authors thank Marco Ajello, M. A. Sánchez-Conde, Seth Digel, and David Williams for fruitful discussions. We thank Soebur Razzaque for useful comments on the draft and the anonymous referee for improving the manuscript. We acknowledge the support of the Spanish MICINN's Consolider-Ingenio 2010 Programme under grant MultiDark CSD2009-00064. J.R.P. acknowledges the support from a Fermi Guest Investigator grant and also from the NASA ATP grants NNX07AGG4G, NSF-AST-1010033, and NSF-AST-0607712.

The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States; the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France; the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy; the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan; and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden.

Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France.

In general, our multiwavelength data are taken from Zhang et al. (2012). Therefore, we refer to the reader to that paper for details. Here, we briefly discuss the SED variability for individual sources. Unless otherwise stated, the LAT data are from the 2FGL (Ackermann et al. 2011), which are integrated between 2008 August and 2010 August, and the 1FHL (M. Ackermann et al., in preparation), which are integrated between 2008 August and 2011 August.

• 1.  
  Mrk 421. The data for the SED of this source are entirely simultaneous, and little variability was detected, although a few small LAT flares are evident (Abdo et al. 2011d). There does not seem to be correlated variability between the X-rays and LAT γ-rays.
• 2.  
  Mrk 501. The SED is from 2006, taken from Anderhub et al. (2009), and is not simultaneous with the LAT data. Variability was of the order of a factor of two in VHE γ-rays, X-rays, and optical. Variability did not appear to be correlated, although the observations are quite sparse.
• 3.  
  1ES 2344+514. Most of the multiwavelength data are taken from 2005–2006 (Albert et al. 2007a) and are not variable, nor are they simultaneous with the LAT data. The LAT data showed minimal variability as well.
• 4.  
  1ES 1959+650. The data are from a multiwavelength campaign in 2006 where the optical and X-ray bands were in high state and showing significant variability. However, the VHE data were at the lowest fluxes ever detected for this source.
• 5.  
  PKS 2005−489. No significant variability on timescale less than a year was found by Aharonian et al. (2005). The X-ray data are from a high state of X-ray emission taken in 1998, whereas the VHE data are from 2009 when the source was in a similar X-ray state as in 1998 (Kaufmann et al. 2009).
• 6.  
  W Comae. The TeV data were taken during a strong VHE flare in 2008 (Acciari et al. 2008).
• 7.  
  PKS 2155−304. All the data including the LAT data are simultaneous, and variability was on the order of a factor of two in the X-rays and VHE γ-rays. The X-ray, VHE γ-ray, and optical variability appears to be approximately correlated.
• 8.  
  H1426+428. No simultaneous broadband data are found for this source.
• 9.  
  1ES 0806+524. No significant variability was found on a timescale of months (Acciari et al. 2009a).
• 10.  
  H2356−309. The TeV and X-ray data were taken simultaneously in 2005 (Abramowski et al. 2010), whereas the data for the optical and other X-ray bands were taken in 2004. In the VHE, significant variations of flux were detected on a timescale of months.
• 11.  
  1ES 1218+304. The optical as well as the X-ray data were taken quasi-simultaneously in 2005. The VHE data were taken between the end of 2008 and middle of 2009. The TeV emission for this source showed day-scale variability.
• 12.  
  1ES 1101−232. A multiwavelength campaign for this blazar was carried out from 2004 to 2006 covering different parts of the X-ray SED showing no variability. The VHE detections were taken at the beginning of 2006.
• 13.  
  1ES 1011+496. The VHE data were taken after an optical outburst in 2007 (Albert et al. 2007b), which may indicate a correlation between the optical and VHE flux. The X-ray data were taken in 2008 May.
• 14.  
  PG 1553+113. No variability for this source has been found, neither in the VHE regime nor in X-ray, by different instruments and campaigns (Reimer et al. 2008; Aleksić et al. 2010).
• 15.  
  3C 66A. The simultaneous observed broadband SED is from Abdo et al. (2011a). The VHE data were taken by MAGIC (Aleksić et al. 2011a) and VERITAS (Acciari et al. 2009b) in similar flux states.