Strong New Constraints on the EBL

Strong New Constraints on the Extragalactic Background Light in the Near- to Mid-IR

M.R. Orr and F. Krennrich Department of Physics and Astronomy, Iowa State University, Ames, IA 50011 E. Dwek Observational Cosmology Lab, Code 665, NASA Goddard Space Flight Center, Greenbelt, MD 20771

Direct measurements of the extragalactic background light (EBL) in the near-IR to mid-IR waveband are extremely difficult due to an overwhelming foreground from the zodiacal light that outshines the faint cosmological diffuse radiation field by more than an order of magnitude. Indirect constraints on the EBL are provided by -ray observations of AGN. Using the combination of the Fermi Gamma-Ray Space Telescope together with the current generation of ground-based air Cherenkov telescopes (H.E.S.S., MAGIC, and VERITAS) provides unprecedented sensitivity and spectral coverage for constraining the EBL in the near- to mid-IR. In this paper we present new limits on the EBL based on the analysis of the broad-band spectra of a select set of -ray blazars covering MeV to several TeV. The EBL intensity at m is constrained to be . We find that the fast evolution and baseline EBL models of Stecker et al. (2006), as well as the model of Kneiske et al. (2004), predict significantly higher EBL intensities in the mid-IR (m) than is allowed by the constraints derived here. In addition, the model of Franceschini et al. (2008) and the fiducial model of Domínguez et al. (2011) predict near- to mid-IR ratios smaller than that predicted by our analysis. Namely, their intensities in the near-IR are too low while their intensities in the mid-IR are marginally too high. All of the aforementioned models are inconsistent with our analysis at the level.

diffuse radiation — galaxies: active — gamma rays: galaxies — infrared:diffuse background

1 Introduction

The extragalactic background light (EBL) contains all radiative energy releases, from nuclear and accretion processes, that have occurred since the epoch of recombination. While the EBL is the second most dominant cosmological radiation field in our universe, surpassed only by the cosmic microwave background (CMB), it has escaped detection through direct measurements at most wavelengths. Its spectrum is bimodal with one component peaking at , originating from radiation released through the formation of heavy elements and the accretion of matter onto black holes in active galactic nuclei (AGNs), and a second component peaking at , created through the absorption of UV/optical light that is re-radiated by dust at infrared wavelengths. The energy spectrum contains, therefore, important information about the cosmic evolution of these energy sources, the obscuring dust, and the relative contributions of starburst galaxies and AGNs to the energy released over cosmic time. In addition, the EBL contains information on the rate of core collapse supernovae which provides an important constraint on the normalization of the diffuse supernova neutrino background (Horiuchi et al., 2009). The trough in the EBL spectrum at is caused by the decrease in the stellar and AGN UV-optical emission towards the mid-IR, and the paucity of hot dust capable of emitting at these wavelengths. This regime is also the most difficult for direct measurements of the EBL due to an overwhelming foreground from zodiacal light. For detailed reviews of the origin, measurements, modeling, and cosmological implications of the EBL see Hauser & Dwek (2001) and Kashlinsky (2005).

Very high energy (VHE) gamma-rays (TeV) interact, via pair production, with photons in the near- to mid-IR regime of the EBL (i.e., ) (Gould & Schréder, 1967). Due to the energy dependence of the - optical depth, the observed VHE emission spectra of blazars are softer than the intrinsic spectra. Information regarding the EBL spectral energy distribution (SED) can be gleaned from this spectral absorption using a variety of approaches. A common technique employed assumes a theoretically based limit on the hardness of the intrinsic spectrum, and hence places a limit on the maximum level of EBL absorption given the softness of the observed spectrum (Aharonian et al., 2006, 2007d). This highest allowed level of absorption translates into an upper limit on the EBL intensity.

Another technique, being an expanded approach to that of Aharonian et al. (1999), utilized the TeV spectra of nearby blazars Markarian 421 and Markarian 501 () (Dwek & Krennrich, 2005). The authors identified EBL scenarios producing an exponential rise in luminosity with energy in absorption-corrected spectra. Such an exponential rise is highly unlikely, based on standard synchrotron self-Compton (Maraschi et al., 1992; Bloom & Marscher, 1996) and external inverse-Compton (Dermer & Schlickeiser, 1993; Sikora et al., 1994) models of blazars, and hence the authors conclude these scenarios are unrealistic.

In this paper we use two different approaches to constrain the EBL, one of which only recently became feasible as a result of Fermi measurements of blazar spectra in the 0.1 to 100 GeV energy regime. At these energies, blazar spectra are typically well characterized by a power-law and are a good representations of the intrinsic spectra, since the EBL causes a negligible amount of attenuation at these energies. The first approach, referred to hereafter as the spectral shape method, or Method 1, assumes that the same power-law characterizing the blazar spectrum at GeV energies extends up to TeV energies, while also allowing for curvature that softens the spectrum at very high energies. Consequently, only EBL SEDs yielding intrinsic blazar spectra consistent with the TeV-extended Fermi power-law are considered viable. The use of the power-law assumption is, of course, only valid over a limited energy regime since the inverse-Compton (IC) component shows a turnover close to its peak energy. However, spectra with IC peaks in the multi-TeV regime are generally well characterized by a power-law below the peak (Samuelson et al., 1998; Aharonian et al., 1999; Krennrich et al., 2001; Aharonian et al., 2002). Similar approaches have been used by Georganopoulos et al. (2010) and Mankuzhiyil et al. (2010) to constrain the EBL absorption optical depth for VHE photons and by Prandini et al. (2010) to constrain the distances of blazars with unknown redshifts.

The second approach, referred to as the TeV spectral break method, or Method 2, makes a less stringent assumption, namely that the intrinsic blazar spectrum at TeV energies is characterized by a single power-law, with no restrictions placed on the spectral index. We then search for evidence of a break at TeV in observed blazar spectra, caused by EBL absorption, and compare this with the spectral breaks predicted by different EBL scenarios. The EBL SEDs producing blazar spectral breaks consistent with observations are considered viable.

We characterize the family of EBL realizations by the m and m intensities, and the ratio between the two. As we show later on, the two approaches are complementary. The spectral shape method is most sensitive to the overall EBL intensity, whereas the TeV spectral break method is mostly sensitive to the relative difference between the m and m intensities, that is, the slope of the EBL between these two wavelengths. The two methods therefore probe different features of the EBL SED, ruling out different regions of parameter space, providing new constraints on the EBL. In Section 2, we discuss the range of EBL parameter space tested. Section 3 describes the two analysis methods used and how they can be combined to improve constraints on the EBL. The results of our analysis are presented in Section 4, with associated interesting and important caveats outlined in Section 5. Section 6 provides some suggestions for improving constraints on the EBL with future blazar measurements. Finally, a discussion of our results, along with concluding remarks, is provided in Section 7.

2 EBL Models

A variety of models exist describing the SED of the EBL (e.g., Kneiske et al. (2002); Primack et al. (2005); Stecker et al. (2006); Franceschini et al. (2008)). The techniques used to generate these models fall into three basic categories: forward evolution, backward evolution, and observed evolution over a particular redshift range. The forward evolution method (e.g., Primack et al. (2005)) makes use of early structure formation scenarios to predict the evolution of galactic luminosity functions forward in time. Backward evolution techniques (e.g., Stecker et al. (2006)) begin with the existing galaxy population and model their luminosity functions backward in time. The third approach to modeling uses deep galaxy surveys (e.g., Kneiske et al. (2002); Franceschini et al. (2008)) or tracers of chemical evolution to infer the cosmic star formation rate and compute the EBL.

The models tested here follow a more observational approach and were derived from what will be referred to throughout the text as the baseline model. This baseline model follows the shape outlined by lower limits derived from galaxy counts obtained with the Hubble Space Telescope (Gardner et al., 2000; Madau & Pozzetti, 2000), the Spitzer Space Telescope (Fazio et al., 2004; Papovich et al., 2004), and the Infrared Space Observatory (Elbaz et al., 2002). An approach similar to that of Mazin & Raue (2007), using third order splines, was implemented to produce this baseline EBL model. The advantage to this technique is that it allows one to easily adjust the shape of the SED through manipulation of the control points defining the spline curve.

Since the absorption of TeV gamma-rays by EBL photons is sensitive to the intensity at near- and mid-IR wavelengths, the regions between and were independently scaled to explore a variety of EBL shapes.111The near and mid-IR regions were scaled as a whole but future work will investigate scenarios using a finer resolution (e.g, including a near-IR excess due to Population III stars (Dwek et al., 2005)). The left panel of Figure 1 illustrates the range of scenarios investigated (shaded region). The baseline model, indicated by the thick solid curve, is shown along with two additional models (dashed and dotted curves) illustrating various levels of scaling in the two wavelength regimes. The right panel of Figure 1 shows the calculated optical depths for gamma-rays given each EBL scenario. All optical depths were calculated assuming the cosmological parameters , , and .

Figure 1: Left: EBL intensity versus photon wavelength. The shaded region indicates the range of scenarios tested. The thick solid line indicates the baseline shape used, from which all other scaled shapes are generated. For clarity, two additional models are shown (dotted and dashed) illustrating the independent scaling of the near- and mid-IR regions. Right: Optical depth (at ) versus gamma-ray energy in TeV for each EBL scenario tested. The optical depths for the baseline and two additional EBL models shown in the left panel are shown as well.

To cover the full range of intensities, 18 scaling factors were used in the near-IR regime and 29 in the mid-IR. This resulted in a total of 519 EBL scenarios.222The total number of scenarios would have been 522 if not for an imposed restriction that the mid-IR intensity be less than the near-IR intensity. This resulted in the exclusion of three models. The scaling factors were chosen to be linearly distributed in .

3 Determining the EBL

The analysis methods used here to constrain the EBL are described in Sections 3.1 and 3.2. In Section 3.3 we describe how these two methods were combined to derive even stronger constraints on the allowable EBL intensity parameter space.

3.1 Method 1 - Spectral Shape Method

With the advent of the Fermi Large Area Telescope (LAT) (Atwood et al., 2009), high energy observations of blazars are now possible in a regime where EBL attenuation is minimal. The operational energy range of the LAT spans GeV, overlapping observations with Imaging Atmospheric Cherenkov Telescopes (IACTs) ranging from TeV (Hinton, 2004; Albert et al., 2008; Holder et al., 2008). Since gamma-rays below GeV travel unimpeded by the EBL (Stecker et al., 2006; Franceschini et al., 2008; Kneiske et al., 2004; Gilmore et al., 2009), the first method described here used spectra measured by Fermi as a proxy for the intrinsic spectra in the TeV regime. Measurements by IACTs were used to calculate a de-absorbed spectrum for each EBL model discussed in Section 2.

The intrinsic TeV source spectra were determined using the relationship


where is the intrinsic spectrum, is the observed spectrum, and is the optical depth at energy and source redshift . The intrinsic spectrum was assumed to have a power-law form given by,333A power-law is generally a good approximation to blazar spectra over a limited energy range (e.g., or ) for measurements with limited statistics.


where is the normalization at energy , is the energy, and is the spectral index. The condition used to evaluate the consistency between the intrinsic TeV spectrum and the extrapolated Fermi spectrum was given by


where and are the calculated IACT intrinsic spectral index and variance, respectively, and are the Fermi spectral index and variance, respectively, and is the confidence level in units of standard deviations. An EBL scenario was considered viable if the de-absorbed spectrum satisfied the set criteria for consistency with the intrinsic spectrum predicted by the Fermi extrapolation. This criterion was that be within of (i.e., Equation 3), where is 1,2, or 3.

As indicated by Equation 3, it was assumed that hard spectrum blazars have no intrinsic spectral break between the Fermi and IACT energy regimes up to several TeV. This assumption was motivated by both observational and theoretical evidence. Firstly, observations of nearby hard spectrum blazars, such as Mrk 421 and Mrk 501, indicate that their observed IC peaks can be located at energies TeV (Dwek & Krennrich, 2005; Aleksić et al., 2010). Since EBL absorption softens the observed spectrum with respect the the intrinsic one, accounting for any absorption whatsoever will only move this peak to higher energies. In fact, Dwek & Krennrich (2005) show that most EBL realizations result in an intrinsic peak energy between and 5TeV for the blazar H 1426+428 (). It is also well known that the spectra of Mrk 421 (Krennrich et al., 2001; Aharonian et al., 2002) and Mrk 501 (Samuelson et al., 1998; Aharonian et al., 1999) have exponential cutoffs at energies TeV, further supporting the assumption of an IC peak located at TeV energies.

Population studies of blazars performed with Fermi show that the hardness of the Fermi spectrum is correlated with the synchrotron and IC peak frequencies (Abdo et al., 2010a, b). In other words, a harder Fermi spectrum generally means a higher IC peak energy. The sources used in Method 1 have some of the hardest spectra of all Fermi-detected blazars (even harder than Mrk 421 and Mrk 501) and hence are likely to have some of the highest energy IC peaks. This provides an additional motivation for the assumption that, for these particular sources, the intrinsic spectrum up to a few TeV is consistent with an extrapolation of the spectrum measured by Fermi.

Costamante et al. (2003), Katarzyński et al. (2006), Aharonian et al. (2007a), and Krennrich et al. (2008) have all shown that extremely hard intrinsic blazar spectra are unavoidable, for these already hard sources, using current EBL model estimates and observationally derived limits. Katarzyński et al. (2006) demonstrated that an IC peak at very high energies can be incorporated into the standard synchrotron self-Compton (SSC) emission model of blazars by introducing a low energy cutoff in the parent electron distribution responsible for the synchrotron and SSC emission. Furthermore, Tavecchio et al. (2009) showed that recent measurements from Swift provide supporting evidence for a very high energy IC peak in the hard spectrum blazar 1ES 0229+200.

The spectrum of 1ES 0229+200 is the only one used here extending to energies that could potentially span the blazar IC peak for the case of hard spectrum blazars. As such, the analysis of 1ES 0229+200 warranted a slightly different treatment.444It should be noted that this “different” treatment was also applied to the other sources used in Method 1. However, it had no impact on the final results. As shown in Figure 2(d), the spectrum of 1ES 0229+200 (Aharonian et al., 2007d) spans the energy range of TeV, nearly 1.5 orders of magnitude. This broad range, combined with the measurements above several TeV, necessitates the inclusion of a term to account for curvature in the intrinsic spectrum. This was done using a log-parabolic function of the form


where is the normalization at energy , is the energy, is the spectral index, and is the log-parabolic fit parameter. Additionally, an F-test was performed, along the lines of Dwek & Krennrich (2005), to test for the presence of an exponential rise in the intrinsic spectrum at the highest energies. This exponential fit had the form


where is the log-parabolic function from Equation 4, is the exponential fit parameter (with units of energy), and is the energy. If the resulting P-value from the F-test was , the exponential function was considered to yield a statistically significant improvement to the fit. The EBL model was then excluded based on the fact that it produced an unphysical intrinsic spectrum. A more detailed discussion of the F-test and its application is found in Dwek & Krennrich (2005).

Using the analysis described above, the hard spectrum blazars 1ES 0229+200 (Aharonian et al., 2007d), RGB J0710+591 (Acciari et al., 2010b), 1ES 1101-232 (Aharonian et al., 2006), and 1ES 1218+304 (Acciari et al., 2010a) have been studied. The combined Fermi and IACT spectra for these sources are shown in Figure 2. It is clear that the spectral breaks between the GeV and TeV regimes, as well as the energy ranges covered, vary between sources. This results in the derivation of a unique set of EBL constraints for each spectrum (Section 4.1). Hard spectrum blazars are ideal for these studies because they are more readily detected at TeV energies (Aharonian et al., 2007d; Krennrich et al., 2008). They also stand the greatest chance of being detected at energies upwards of 10TeV, broadening the EBL constraints into the mid-IR regime. Combining the analyses of multiple spectra further constrains the allowable parameter space. In the following, we refer to this approach as Method 1, or the spectral shape method.

(a) 1ES 1218+304 (VERITAS (Acciari et al., 2010a)).
(b) 1ES 1101-232 (H.E.S.S. (Aharonian et al., 2006)).
(c) RGB J0710+591 (VERITAS (Acciari et al., 2010b)).
(d) 1ES 0229+200 (H.E.S.S. (Aharonian et al., 2007d)).
Figure 2: Combined Fermi and IACT spectra in units of . The best fit Fermi spectrum is indicated by the shaded region. The Fermi flux points were calculated by fixing the spectral index to the value obtained from the fit over the full energy range and then fitting the integral flux over each energy bin. The IACT spectra are given by the red line and flux points in each plot. Note: Fermi has a weak detection of 1ES 0229+200 (the current result with months of data is ). For this analysis we have fixed its spectral index in the Fermi regime to a value of .

3.2 Method 2 - TeV Spectral Break Method

The absorption of gamma-rays by EBL photons may produce breaks in the observed blazar spectrum for some SEDs. This break is produced by changes in the slope of the - optical depth. As illustrated in the right panel of Figure 1, the optical depth calculated for certain EBL scenarios is nearly flat in the energy range of TeV resulting in an approximately energy independent absorption of gamma-rays. The observed spectral index in this energy range will be closer to the intrinsic value – producing a break in the spectrum at TeV. The magnitude of this break increases with the source redshift and depends on the near- to mid-IR ratio (Imran & Krennrich, 2008). With currently available data, no individual blazar spectrum is likely to show a statistically significant break in its spectrum. However, a large sample of blazars may reveal a redshift dependent trend. The presence/absence of such a trend in observations can be used to constrain the EBL. Table 1 lists the sample of blazars used to perform this study.

Source Name Redshift Spectral Index Method Reference # Spec. Points
Used l.t./g.t. TeV
1ES 2344+514 0.044 2 Albert et al. (2007a) 4/3
1ES 1959+650 0.048 2 Tagliaferri et al. (2008) 4/2
PKS 0548-322 0.069 - 2 Aharonian et al. (2010) 3/2
PKS 2005-489 0.071 2 Aharonian et al. (2005a) 6/3
RGB J0152+017 0.080 - 2 Aharonian et al. (2008) 4/2
PKS 2155-304 0.117 2 Aharonian et al. (2005b) 7/3
RGB J0710+591 0.125 1,2 Acciari et al. (2010b) 3/2
H 1426+428 0.129 2 Petry et al. (2002) 3/4
1ES 0229+200 0.140 1,2 Aharonian et al. (2007d) 3/5
1ES 1218+304 0.182 1,2 Acciari et al. (2010a) 7/2
1ES 1101-232 0.186 1,2 Aharonian et al. (2006) 9/4
1ES 0347-121 0.188 - 2 Aharonian et al. (2007b) 4/3
* Fermi analysis performed in this work.
probability of being a steady Fermi source.
Fermi has a weak detection of 1ES 0229+200 at . For this analysis we have assumed a Fermi spectral index of .
Table 1: Blazar sample used in the described analyses. The columns are as follows (left to right): source name, source redshift, GeV and TeV spectral indices, analysis technique performed with the source spectrum, reference to the TeV spectrum, number of TeV spectral points less than (l.t.) and greater than (g.t.) TeV, and the measured spectral break at TeV. All errors given are statistical only. Method 1 refers to that discussed in Section 3.1 while 2 refers to Section 3.2. Unless otherwise indicated, is taken from the Fermi 1 Year Catalog (

The shape of the EBL uniquely determines the spectral break versus redshift distribution and was calculated here for a multitude of scenarios using a test blazar spectrum. This test spectrum was chosen to represent the “average” spectrum of the blazars listed in Table 1. This was motivated, for example, by the fact that the spectrum of 1ES 0229+200 (Aharonian et al., 2007d) ranges from TeV whereas that of 1ES 1218+304 (Acciari et al., 2010a) spans TeV. Using a test blazar spectrum ranging from TeV would be unrealistic as there is no individual source in the sample covering such a wide energy range. To generate a spectrum representative of the overall sample of blazars, the average lowest energy bin, highest energy bin, and number of bins in the sample were calculated. This resulted in a test spectrum characterized by 8 energy bins ranging from 200GeV to 5TeV. The flux normalization and spectral index (intrinsic values) of the test blazar spectrum are irrelevant since only the magnitude of the spectral break is of concern. Nevertheless, values need to be chosen to perform the analysis. The Crab Nebula flux at 1 TeV was used for the normalization along with an intrinsic spectral index of 1.5. The underlying assumption in this method is that the intrinsic blazar spectrum is well described by a single power-law over the energy range considered.

Beginning with an exact power-law function, and calculating the spectrum due to EBL absorption, one is left with an observed spectrum that shows detailed features of the - optical depth and is decidedly not of a power-law (or broken power-law) form. This is, of course, not what is observed by IACTs. To produce an absorbed spectrum more consistent with observations, each flux point was assigned an error of 25% of the flux in that bin. This value was chosen to approximately coincide with the mean error of all flux points from the IACT measurements in the blazar sample. Gaussian fluctuations were added to the spectral points using a Normal distribution with a standard deviation equal to the 25% error bars. Optical depths for the test spectrum were calculated for all EBL scenarios over the redshift range 0.05–0.45 in steps of 0.05. No EBL evolution with redshift was implemented.

The absorbed spectra were calculated using the inverse relation of Equation 1. Each spectrum was fit with a broken power-law of the form


where is the normalization at the break energy , and are the spectral indices below and above , respectively, and is the energy. The spectral break was defined as


With this definition, a spectral hardening above the break energy would result in a positive . From here the distribution of as a function of redshift was determined for each EBL scenario. This process was repeated 100 times, mimicking multiple independent observations, each time yielding a different result due to the Gaussian fluctuations in the spectra. This facilitated the determination of the break versus redshift distribution to an arbitrary precision and accuracy.

Each distribution was fit with a line, i.e.,


where is the spectral break at redshift and and are free parameters. Using the best linear fit obtained from observational data, the , , and two dimensional contours in slope and intercept space (i.e., and , respectively, in Equation 8) were determined. The best fit results from the expected spectral break versus redshift distribution for each EBL scenario were then compared with the observationally derived contours. If the fit parameters for a particular scenario fell within the two dimensional contour, this scenario was included as part of the contour in EBL parameter space. The and contours in EBL parameter space were determined in the same fashion (Section 4.2). In the following, we refer to this approach as Method 2, or the TeV spectral break method.

3.3 Complementarity of Methods 1 & 2 in Constraining the EBL

The constraints placed on the EBL are presented in the form of a contour plot. This contour is drawn in the parameter space defined by the near-IR intensity at and the mid-IR intensity at . These two wavelengths are representative of the positions of the near-IR peak and mid-IR trough in the EBL SED. Given that the EBL models tested are relatively smooth, these two parameters serve as an adequate characterization of a given EBL scenario. An illustration of such a contour plot is shown in Figure 3. Each point represents one of the EBL scenarios tested. Shown in this fashion, the linear sampling of intensity in log-space is clearly seen. A constant ratio of intensities (i.e., ) would be given by a straight line with positive slope. It should be noted that parametrizing the results in this fashion does not mean the derived constraints only apply to the EBL intensity at and . This presentation is simply a way of characterizing the average slope of the EBL SED between these two wavelengths.

Figure 3: EBL intensity at plotted versus the intensity at . Each grey dot represents an EBL scenario tested. (a) As the intensities at and increase (moving up and to the right on the plot), the inferred intrinsic blazar spectrum becomes harder. The intrinsic spectrum softens as the two intensities decrease (down and to the left). This is indicated by the arrows in the figure. By placing restrictions on the maximum hardness and softness of the intrinsic spectrum, one obtains a contour in intensity space oriented approximately perpendicular to these two directions (shaded region). (b) Increasing the near-IR to mid-IR ratio increases the redshift dependence of the EBL-induced spectral break at TeV energies (i.e., becomes larger). Decreasing this ratio results in an increasingly negative dependence of spectral break versus redshift (i.e., becomes smaller to the point of being negative). This is indicated by the arrows in the figure. By comparing the redshift dependence of the TeV spectral break produced by various EBL models with observations, one obtains a contour in intensity space oriented as shown by the shaded region.

The shaded region in Figure 3 shows how a contour, derived from the analysis outlined in Section 3.1, might appear. As the near- and/or mid-IR intensities increase (decrease), the inferred intrinsic spectrum hardens (softens). A contour oriented as in Figure 3 is obtained when placing limitations on the maximum hardness/softness of the intrinsic spectrum.

An approximately orthogonal contour is obtained from the analysis outlined in Section 3.2. This is because, when increasing the near-IR to mid-IR intensity ratio, the dependence of the EBL induced spectral break with redshift becomes increasingly positive. That is to say, becomes larger. Conversely, decreasing the near-IR to mid-IR ratio makes smaller to the point where it can become negative. These directions are shown on the contour illustration in Figure 3. By comparing the redshift dependence of the TeV spectral break produced by different EBL models with that from observations, one can obtain a contour in EBL intensity space as represented in Figure 3 by the shaded region.

Comparing the two contours from Figure 3, it becomes clear that the EBL constraining methods outlined in Sections 3.1 and 3.2 are complementary to one another. Each method is sensitive to different characteristics of the EBL SED thereby improving the constraints at near- and mid-IR wavelengths.

4 Results

The techniques described in Sections 3.1 and 3.2 have been applied to a sample of blazars (Table 1). The four blazars selected for the spectral shape analysis were chosen for their redshifts (), hard spectra, and steady emission in the Fermi regime.

The blazar sample used in the TeV spectral break analysis was selected based on characteristics of the source spectra in the TeV regime. Most importantly, to test for a spectral break at TeV, each spectrum must have at least two data points both above and below the break energy. Spectra from large flaring events were excluded as their atypically high statistics skew the overall results. An example of this was the Crab flare event of PKS 2155-304 (Aharonian et al., 2007c). The two nearby blazars Mrk 421 and Mrk 501 were also excluded from this analysis due to the fact that the presence of a cutoff at TeV in their intrinsic spectra (Krennrich et al., 2001; Samuelson et al., 1998), combined with their low redshift (), greatly inhibits the search for a spectral break resulting from EBL absorption.

The analysis results for Methods 1 and 2 are presented in Sections 4.1 and 4.2, respectively, and are then combined in Section 4.3.

4.1 Spectral Shape Analysis

The source spectra used for the spectral shape analysis are shown in Figure 2. Each Fermi spectrum consists of approximately 19 months of data and was analyzed using the Science Tools v9r15p2 package distributed by the Fermi collaboration.555 Only events with zenith angles less than were included to avoid contamination from Earth’s gamma-ray albedo. All exposures and fluxes were calculated using the post-launch instrument response function P6_V3_DIFFUSE. The galactic and extragalactic diffuse gamma-ray backgrounds were modeled using data available on the Fermi Science Support Center website.666 The variability index of each source, as quoted in the Fermi 1 Year Catalog777, is less than 23.21 indicating a probability the source emission is steady. Spectral fitting was performed from MeV to the highest photon energy associated with the given source. TeV spectra were obtained from the references listed in Table 1.

(a) 1ES 1218+304
(b) 1ES 1101-232
(c) RGB J0710+591
(d) 1ES 0229+200
Figure 4: Constraints, provided by Method 1, on the EBL intensity, in units of , at m and m. The contours shown are for the 1 (red), 2 (blue), and 3 (green) sigma confidence intervals. Each grey dot represents one of the EBL scenarios tested. Note: Fermi has not obtained a significant detection of 1ES 0229+200. For this analysis we have assumed the spectral index in the Fermi regime to be .

The Fermi analysis of 1ES 0229+200 was handled differently than for the other three sources due to the fact that Fermi has only a weak detection of this source. Nevertheless, it was included since it has the hardest TeV spectrum of all the blazars studied here and provides a unique set of constraints. To obtain a spectrum from the Fermi data, an assumed spectral index of was used. This parameter was fixed in the source model so that the the best fit flux normalization could be determined. The overall test statistic (TS) for this source, between 1 and 28GeV, was 18.0 (). The TS values for the first (GeV) and second (GeV) energy bins were 2.6 and 15.4, respectively.

The , , and contours for each blazar, constraining the EBL intensity at m and m are shown in Figure 4. The reader is cautioned not to view these contours as strict exclusion regions for other authors’ models since the analysis is dependent on the shape of the EBL SED. However, given that most models have the same general shape, it is a safe assumption that scenarios falling outside the contours derived here are not consistent with this analysis.

Figure 4 shows that gamma-ray observations provide strong constraints on the m EBL intensity, but leave the m intensity largely undetermined. This is due to the lack of TeV data extending up to multi-TeV energies which provide the strongest constraints on the mid-IR region of the EBL. It is worthwhile to point out some particular properties of the different contours in Figure 4 as they provide insight into the constraints derived with each blazar.

One of the most noticeable features is the narrowness of the 1ES 1218+304 contours (Figure 4(a)). This is a result of the high precision measurement of the spectral index by Fermi and VERITAS in both the GeV and TeV regimes (see Table 1). The measurement errors are a factor of 2 or more smaller than those of the other blazars in the sample. Another distinguishing feature is the vertical orientation of the contours indicating that the TeV spectrum is minimally affected by EBL absorption in the mid-IR regime. This results from the fact that the spectrum only extends up to TeV.

The spectrum of 1ES 1101-232 extends up to TeV but is also insensitive to the mid-IR portion of the EBL (Figure 4(b)). This is most likely due to the somewhat large error on the spectral index measured by Fermi () as well as the large power-law fit probability () to the observed TeV spectrum (see Aharonian et al. (2006)).888A fit probability of 80% suggests the possibility of either overestimated errors or somewhat correlated data points. In either case, fluctuations in the highest energy spectral points (those sensitive to the mid-IR portion of the EBL) will have a limited impact on the overall power-law fit to the de-absorbed spectrum.

In contrast to 1ES 1218+304 and 1ES 1101-232, the contours for RGB J0710+591 (Figure 4(c)) are somewhat inclined, indicating sensitivity to the mid-IR EBL intensity. The spectrum of RGB J0710+591 extends up to TeV demonstrating that a modest increase in the maximum energy of the spectrum can have a significant impact on EBL studies.

The spectrum of 1ES 0229+200 is unique in that it starts at GeV and extends up to TeV. As a result, this source is the most sensitive, of all those analyzed here, to absorption by the mid-IR portion of the EBL. This is evident in Figure 4(d). While none of the aforementioned blazars strongly constrain the m EBL intensity, the spectral shape analysis of 1ES 0229+200 constrains the intensity to be .

4.2 TeV Spectral Break Analysis

The energy dependence of the - optical depth for some EBL scenarios may produce a break around TeV in blazar spectra (Section 3.2). Here, a spectral break energy of TeV is assumed. The choice of this value is easily arrived at using the the fact that the cross section peaks at , where is the EBL photon wavelength and is the gamma-ray energy (Dwek & Krennrich, 2005). Since the near-IR peak of the EBL occurs at m, a flattening of the - optical depth would be expected at TeV. The right panel of Figure 1 shows that the optical depth does, in fact, flatten out between TeV, further motivating this choice for the spectral break energy.

To determine the spectral break versus redshift distribution (), a sample of 12 blazars (Table 1) with spectra extending both above and below the break energy were fit with a broken power-law (Equation 6). Figure 5 shows the distribution of versus the source redshift for all blazars in the sample. Two linear fits were performed on the data – one leaving the intercept as a free parameter and one with the intercept fixed to 0. The results of these two fits were


In the absence of any spectral break one would expect . The data exclude this hypothesis at the level.

Figure 5: Distribution of the observed spectral break versus redshift assuming a break energy of 1.3TeV (the choice of which is motivated by the near-IR peak at m in the EBL scenarios used here). The sample includes the 12 blazars (see Table 1) available to date with spectra extending above and below the break energy. The spectral break is defined as the spectral index below the break energy minus the spectral index above the break energy (i.e., ) The black solid line represents the best linear fit given by and yields . The grey solid line shows the best linear fit assuming , giving with . The long-dashed line shows the distribution expected from the scenario of Aharonian et al. (2006), scaled by a factor of 0.55, yielding . The short-dashed line shows the distribution expected from the model of Franceschini et al. (2008), yielding . Fitting the data with a flat line at yields .

Figure 5 also shows the distributions produced by the models of Aharonian et al. (2006),999The EBL SEDs from Aharonian et al. (2006) are scaled versions of the model from Primack et al. (2001). scaled by a factor of 0.55, and Franceschini et al. (2008). A summary of how well these models describe the data, along with the other authors’ models considered here, is given in Table 2. Column 3 gives the best linear fit to the distribution produced by the given model. This is derived in the same way outlined in Section 3.2. However, in this case we use the EBL model from the reference listed rather than one of our own realizations discussed in Section 2. With the exception of the EBL scenario from (Aharonian et al., 2006), all models are inconsistent with TeV spectral break data at the level.

Model Reference    Best Fit
Aharonian Aharonian et al. (2006) 14.16 (22.95) ()
Aharonian Aharonian et al. (2006) 18.09 (21.97) ()
Aharonian Aharonian et al. (2006) 19.06 (23.24) ()
Stecker (fast) Stecker et al. (2006) 31.66 (48.54) ()
N/A 31.98
Dominguez Domínguez et al. (2011) 34.46 (31.49) ()
Franceschini Franceschini et al. (2008) 36.23 (32.67) ()
Kneiske Kneiske et al. (2004) 36.46 (34.54) ()
Stecker (baseline) Stecker et al. (2006) 61.87 (52.86) ()
Table 2: Summary of the comparison between TeV spectral break data and the EBL models from other authors considered here. The columns are as follows: Column 1: EBL model name. Column 2: Reference. Column3: Best linear fit, of the form , to the TeV spectral break distribution produced by the model (i.e., simulated, not observed data). Column 4: of the observed data assuming the best fit from Column 3. Column 5: Confidence level with which the model can be excluded given the data. Note: Parenthetical numbers in Columns 4 and 5 are obtained using the technique discussed in Section 5.2.

Figure 6 shows the , , and contours resulting from comparison between the distributions generated from the EBL models tested and observations (Figure 5). In addition, points taken from Aharonian et al. (2006), Kneiske et al. (2004), Stecker et al. (2006), Franceschini et al. (2008), and Domínguez et al. (2011) show where these other EBL models fall within the contours generated here.101010It should be noted that some of these EBL models have bumps at m. The data points should therefore only be used as a point of reference and not for a strict comparison. The compatibility of these models with our results was determined using the full EBL SED. The orientation of the contours are as anticipated from the discussion in Section 3.3. The limits on the near- to mid-IR ratio () become increasingly restrictive as the intensities at m and m increase. Using this technique, future direct measurements of the EBL at near-IR wavelengths (for instance by the James Webb Space Telescope111111 will help to further constrain the EBL in the mid-IR regime.

Figure 6: Constraints, provided by Method 2, on the EBL intensity, in units of , at m and m. The contours shown are for the 1 (red), 2 (blue), and 3 (green) sigma confidence intervals. The intensities for the models of Aharonian et al. (2006) (red stars), Stecker et al. (2006) (green triangles), Franceschini et al. (2008) (blue diamond), Kneiske et al. (2004) (orange square), and Domínguez et al. (2011) (purple hexagon) are also shown. The three red stars for the scenario of Aharonian et al. (2006) represent scaling factors of 1.0, 0.55, and 0.45 (right to left). For the model of Stecker et al. (2006), the two green triangles represent the so called fast and baseline evolution models (right to left).

4.3 Combining the Analyses

Improved constraints on the EBL can be obtained by combining the results of the analyses from Sections 4.1 and 4.2. Figure 7 shows the overlay of (a) and (b) contours from the analysis of RGB J0710+591, 1ES 1218+304, and . The contours from 1ES 1101-232 and 1ES 0229+200 are left off because they do not further restrict the constraints. Figure 8 shows the final allowed and regions of parameter space. The overlapping region indicates a very low EBL intensity at m. It should be reiterated that the assumption here was that the intrinsic TeV spectrum was approximately equal to the Fermi spectrum. If, in fact, there is an intrinsic break between the GeV and TeV spectra of blazars, the intrinsic TeV spectrum will be softer than is assumed here. This will move the Fermi-IACT derived EBL contours toward the left of Figure 8 and hence allow for lower EBL intensities at m. However, given the derived contour (Figure 6), the EBL intensity at m will then be constrained to even lower values.

(a) Overlayed contours.
(b) Overlayed contours.
Figure 7: Combined contours from the analysis of RGB J0710+591, 1ES 1218+304, and . The contours of 1ES 1101-232 and 1ES 0229+200 are left off because they do not further constrain the region of allowable parameter space beyond that shown in the figure. The separately colored shaded area in each panel indicates the region of overlap for the three contours.
Figure 8: Combined and contours from the analysis of RGB J0710+591, 1ES 1218+304, and . The contours of 1ES 1101-232 and 1ES 0229+200 are left off because they do not further constrain the region of allowable parameter space beyond that shown in the figure. Also shown for reference are the models of: Aharonian et al. (2006) (stars) scaled by 1.0, 0.55, and 0.45 (right to left); Stecker et al. (2006) (triangles) fast evolution (right) and baseline (left) models; Franceschini et al. (2008) (diamond); Kneiske et al. (2004) (square); Domínguez et al. (2011) (hexagon). The dotted line at indicates the total integrated EBL intensity of Hopwood et al. (2010). The dashed lines at and indicate the error on the value.

Figure 9 shows the confidence region for the EBL SED. The models of Aharonian et al. (2006), Stecker et al. (2006), Franceschini et al. (2008), and Kneiske et al. (2004) are also shown. This region is generated using the rectangular area encompassing the contour in Figure 8.

Figure 9: Approximate EBL SED confidence region (blue shaded area) derived from Figure 8. Also shown are the models of: Aharonian et al. (2006) (dotted lines) scaled by 1.0, 0.55, and 0.45; Stecker et al. (2006) (dot-dashed lines) fast evolution (upper) and baseline (lower) models; Franceschini et al. (2008) (dashed line); Kneiske et al. (2004) (upper solid line); Domínguez et al. (2011) (lower solid line and grey shaded area). The filled symbols indicate the galaxy count lower limits of (left to right): Gardner et al. (2000) (downward triangles); Madau & Pozzetti (2000) (upward triangles); Keenan et al. (2010) (squares with error bars); Levenson & Wright (2008) (sideways triangle); Fazio et al. (2004) (squares without error bars); Elbaz et al. (2002) (hexagon); Hopwood et al. (2010) (circle); Papovich et al. (2004) (diamond). The at m indicates the tentative detection of Dwek & Arendt (1998). Note: The partially transparent blue shaded area in the mid-IR shows the expansion of the confidence region when using the alternate test spectrum described in Section 5.2.

5 Caveats and Considerations Regarding Methods 1 & 2

There are a number of caveats pertinent to the analyses presented here. While each of these cannot be discussed exhaustively, we address the most relevant topics below. First, however, we make a few general remarks.

The EBL attenuation factors calculated throughout the paper were done so using the central energy (in logarithmic space) of each spectral energy bin. Another way of doing this would be to calculate the average attenuation across the entire energy bin. We have verified numerically that these two approaches result in absorption factors within of one another given the typical energy bin widths encountered here.

Not every possible shape of the SED was explored with the EBL parametrization used here. The aim of this work is to constrain the overall shape of the EBL and not to test for various bumps and wiggles. The reader is therefore cautioned not to view any individual data point, derived from galaxy counts, in Figure 9 as being inconsistent with our results if it lies outside the designated contour region. This analysis does not rule out the possibility that there are, for instance, bumps in the EBL spectrum (e.g., at m). What is demonstrated by this analysis is that the near- to mid-IR ratio of the EBL appears to be larger than is predicted by many models thereby requiring a relatively high intensity in the near-IR (see Section 7).

5.1 Hard Spectrum Blazars and TeV IC Peak Energies

The assumption made in Method 1 was that hard spectrum blazars have IC emission components peaking at TeV energies. We believe this assumption is well motivated (see Section 3.1) but, if it were to prove incorrect, and one instead allowed for some level of intrinsic curvature, the lower limits obtained here in the near-IR would be pushed to lower intensities. One could also interpret the upper boundary of the contour in Figure 9 as an upper limit. However, there is a very important caveat to this second interpretation. If the near-IR intensity decreases, the mid-IR intensity must also decrease so that the overall slope of the EBL SED between near- and mid-IR wavelengths remains within the allowed range determined with Method 2.

It should also be noted that the contours shown in Figure 4 do not all overlap in a common region of EBL parameter space. This could be interpreted as being due to intrinsic differences between some of the blazars (e.g., different intrinsic spectral breaks). However, given that this is a effect, and a common overlapping region for the contours can be identified, we do not consider this issue to be of significance.

5.2 Spectral Break Versus Redshift Calculation

There are multiple ways of to go about calculating the spectral break versus redshift distribution for a particular EBL scenario. The method implemented here (Section 3.2) used a test blazar spectrum representative of the average spectrum of the 12 blazar data sample. This determination governed both the energy range covered and the errors on each flux point. An absorbed spectrum was then calculated for a uniformly spaced set of redshifts.

An alternative to this approach is to use 12 different test blazar spectra, with each one being representative of a particular blazar in the dataset. In this case, one would calculate the absorbed spectrum, at each redshift of the observed blazars, using a test spectrum with energy bins and error bars equivalent to the associated observed blazar spectrum. This technique is most appropriate when testing EBL models that contain features in the SED (e.g., the bump at m in the model of Domínguez et al. (2011)). In cases such as this, the spectral break calculated, due to EBL absorption, can be sensitive to the energy range over which one performs the fit. Using the same energy range and binning of the spectral observations allows one to determine the expected spectral break given the specific dataset and an assumed EBL scenario. This technique is also preferable when the mid-IR EBL intensity is quite high (e.g., the Stecker et al. (2006) fast evolution model).

We have used the technique described above to compare the various EBL realizations of other authors with the spectral break versus redshift data. The exclusion probabilities obtained are shown in Table 2 along side those calculated using the simple best linear fit to the spectral break versus redshift distribution for each scenario. With the exception of the fast evolution model of Stecker et al. (2006), the values obtained using the two techniques are comparable. The reason for the difference in the results of the Stecker et al. (2006) fast evolution model is likely due to the fact that it has the highest EBL intensity, and therefore the largest amount of absorption, which makes it the most sensitive to the energy range over which the absorbed spectra are fit.

As an additional cross check of Method 2, we performed the exact same analysis described in Section 3.2 using a test blazar spectrum with 6 spectral points ranging in energy from GeV to TeV. This more limited energy range around TeV elucidates how the low and high energy portions of the spectrum affect the spectral break calculation. The contours obtained using this test spectrum are shown in Figure 10. The main effect of the reduced energy range is to shift the contours upward to slightly higher mid-IR intensities. This is due to the fact that, by excluding energies above TeV, the spectral fit is less impacted by the sharp rise in absorption, above a few TeV, that can occur for high mid-IR EBL intensities. This results in the calculation of a spectral break versus redshift distribution with a more positive, or less negative, slope which allows for higher mid-IR intensities. The corresponding contour region for the EBL SED is shown in Figure 9 as the partially transparent blue shaded region.

Figure 10: Constraints on the EBL intensity, provided by Method 2, using a test spectrum with 6 spectral points ranging from GeV to TeV. See Figure 6 for a full description of the contours and data points shown.

It is also important to understand the effects of each individual source on the observed spectral break versus redshift distribution shown in Figure 5. Namely, how do the results change if sources lying outside the fit (e.g., PKS 0548-322 and H 1426+428) are removed? By excluding H 1426+428 from the dataset, the fit result of Figure 5 changes by only . Removing PKS 0548-322 from the dataset results in a best fit of . In this case, both the slope and intercept of the fit change by less than while the exclusion significances of the various EBL models mentioned in the paper (e.g., Primack et al. (2005), Stecker et al. (2006), etc.) change only marginally.

5.3 Evolution of the EBL with Redshift

The evolution of the EBL as a function of redshift should be taken into account when calculating the gamma-ray absorption for sources with large redshifts. Since the quality of the data used, and the nature of the analysis performed, will impact exactly at what redshift EBL evolution becomes a necessary inclusion, we investigate its effects here.

The EBL is comprised of photons generated by sources that have been evolving since the early Universe. This evolution is, therefore, imprinted on the EBL but has been neglected in our analysis thus far. This means that, when calculating the absorption for a source at a redshift of, say, , we have been assuming that all of the target EBL photons are in place at this redshift and that there is no absorption and/or reemission of these photons (i.e., no evolution) between and . This assumption is reasonable for small redshifts, but one must determine at what redshift it is no longer warranted.

The evolution of the EBL has been examined in detail by many of the authors whose models are used here (e.g., Kneiske et al. (2004) and Franceschini et al. (2008)). There is also a simple first order approach to incorporating EBL evolution into absorption calculations. This method simply applies a correction factor to the cosmological photon numbers density scaling (e.g., Madau & Phinney (1996) and Raue & Mazin (2008)). This factor is chosen such that the resultant EBL evolution is in good agreement with more detailed models (e.g., Kneiske et al. (2002)). Raue & Mazin (2008) have shown that this approach can be used to match the EBL evolution well for redshifts of .

Using this approach, we have investigated the effects of EBL evolution on the results presented here. The inclusion of EBL evolution lowers the calculated optical depths by at a redshift of . For the blazars analyzed in Method 1, this translates to a softening of the intrinsic spectra by . This corresponds to a shift of , to higher EBL intensities, in the contours of Figure 4.

Method 2 is insensitive to the intrinsic spectral softening effects due to the incorporation of EBL evolution since this method relies on the measurement of the spectral break and not the spectral index itself. In other words, it depends on a relative measurement rather than an absolute one. It should be noted, however, that this is contingent upon the (reasonable) assumption that only the intensity of the EBL evolves over such small redshifts and not the shape of the SED. In other words, the relative contributions to the EBL from different emission components remain approximately the same over small redshifts.

5.4 Alternatives Scenarios Describing Blazar Spectral Observations

Alternative scenarios that could explain the observed spectra of blazars, while allowing for a wider range of EBL intensities, have been investigated by others. Some have proposed that the observed high energy gamma-ray emission from blazars could be dominated by secondary photon emission produced along the line of site by cosmic-ray interactions with EBL photons (Essey & Kusenko, 2010; Essey et al., 2010a, b). Essey et al. (2010a) have shown that such a scenario can be used to describe the observed spectra of many VHE blazars. Other authors have suggested that photons could oscillate into light axion-like particles (ALPs) via interactions with intergalactic magnetic fields (De Angelis et al., 2007, 2009). These ALPs, which travel unimpeded through the Universe, can then convert back into photons before reaching the Earth. Since ALPs are not affected by EBL absorption, this process would allow for an EBL intensity that is, in actuality, higher than would appear based on standard absorption calculations. These are just two of the many interesting avenues of research being pursued in relation to the EBL.

6 Improving Constraints on the EBL

The constraints on the EBL presented here can be improved through further observations with Fermi and IACTs. Fermi is continuously monitoring the entire sky and hence data on the sources used in Method 1 are being acquired all the time. Improved measurements in the TeV regime using current instruments will require either deep or flaring state observations of blazars and will predominantly impact the constraints obtained with Method 2. As can be seen from Figure 5, the majority of the blazar spectral breaks measured using current data are not, on an individual basis, statistically significant. The driving factor behind the large measurement errors is an insufficient knowledge of blazar emission at energies TeV for most of the distant sources. As deeper exposures are obtained, and source photon statistics continue to improve, measurements of many blazars will extend to ever higher energies until their spectra either intrinsically, or via EBL absorption, cut off.

As spectra span a larger energy range, and as the precision of each flux point improves, it is less likely that their intrinsic emission will be well described by a single power-law. When looking for evidence of a break at TeV, it would be necessary to modify the approach used here by either truncating the fit regime so as to exclude energies above or below where the IC peak is believed to be located, or use a more complicated fit function that includes (at least) a spectral break parameter along with a curvature parameter. It will also likely be necessary to use the more sophisticated spectral break versus redshift determination described in Section 5.2.

Another consideration, regarding the detection of blazars in flaring states, is that of spectral variability. Some blazars exhibit spectral hardening with increasing flux (Krennrich et al., 2002; Albert et al., 2007b, c) while others show no evidence for spectral variability when flaring (Acciari et al., 2010a). Any change in the overall slope of the spectrum would have a limited impact on Method 2, given that it relies on the relative change in spectral index (spectral break), but the limits obtained using the Fermi and IACT data as in Method 1 could be affected.121212The four sources used in Method 1 show no signs of spectral variability. If the TeV spectrum were to harden, while the Fermi spectrum remained unchanged, tighter constraints on the EBL could be obtained due to the smaller observed spectral break. However, if spectral variability takes place entirely above or below the IC peak, for hard spectrum blazars, where we have assumed that the IC peak is at multi-TeV energies, both Fermi and IACTs should see the same amount of spectral hardening thereby having no impact on the derived constraints. Simultaneous Fermi and IACT data are needed during both flaring and quiescent states to truly evaluate the effect (if it exists) of spectral variability on hard spectrum blazars.

7 Discussion & Conclusions

We have presented constraints on the EBL intensity using two distinct methods for the analysis of blazar spectra. The spectral shape method (Method 1) used Fermi observations of the three hard spectrum blazars RGB J0710+591, 1ES 1101-232, and 1ES 1218+304 as a proxy for the intrinsic source emission in the sub-TeV to TeV energy regime. In addition, 1ES 0229+200, which has a weak Fermi detection (), was used in this analysis, with an assumed Fermi spectral index of 1.5, as its elucidates the impact of sources detected above TeV on EBL constraints. Observed IACT spectra were corrected for attenuation due to different EBL realizations. Limits on the EBL were then placed by rejecting realizations that did not produce deabsorbed (i.e., intrinsic) IACT spectra consistent with the Fermi extrapolation.

The second method, the TeV spectral break method, involved two steps. The first consisted of characterizing the spectral break, occurring at TeV, using the IACT observations of 12 blazars located at different redshifts. The magnitude of the spectral break was then plotted as a function of redshift (Figure 5). We then assumed that the intrinsic blazar spectrum can be characterized as a power-law which, unlike Method 1, was not determined by the Fermi spectrum, and is therefore more general in nature. The second step consisted of comparing the data in Figure 5 with the redshift distribution of spectral breaks derived from the attenuation of an intrinsic test spectrum using different EBL realizations.

We showed that Method 1 primarily constrains the overall EBL intensity, whereas Method 2 primarily constrains the relative EBL intensities in the near- and mid-IR (m and m, respectively). The two approaches are therefore complementary and provide a closed confidence range for the EBL intensity. This is the first time these methods have been used to derive such a contour for the EBL SED.

The resulting constraints imposed by both methods on the the m and m intensity of the EBL are shown in Figure 8. Using ultra deep m mapping observations of the gravitational lensing cluster Abell 2218 with the AKARI infrared satellite, Hopwood et al. (2010) derived an intensity of for the integrated galaxy light (IGL) down to an intensity of mJy. A model fit to the differential contribution to the IGL, combining data from several fields, gives a total integrated EBL intensity of (Hopwood et al., 2010). The IGL seems to have totally resolved the EBL at m. Using our constraints, and the claimed m EBL detection, we derive a value of for the EBL intensity at m, the error representing a uncertainty. Comparison with the IGL (Madau & Pozzetti, 2000) suggests that , and perhaps as much as (Keenan et al., 2010), of the EBL has been resolved at this wavelength.

Figure 9 presents the confidence range of the EBL intensity derived in this paper as a function of wavelength (shaded blue). Our derived constraints are consistent with claimed EBL detections at m (Dwek & Arendt, 1998) and m (Levenson & Wright, 2008). Figure 9 also shows that our EBL constraints are consistent with some models, such as the scaled version of Aharonian et al. (2006) and the most recent model of Hopwood et al. (2010). However, the models of Stecker et al. (2006), Franceschini et al. (2008), Kneiske et al. (2004), and Domínguez et al. (2011) are disfavored by our analysis at confidence level.

By inspection of Figure 9, one may be tempted to conclude that the models of Franceschini et al. (2008) and Domínguez et al. (2011) are inconsistent with the analysis presented here predominantly because their intensities in the near-IR are too low. One might then suspect that this is purely due to constraints placed on the EBL by Method 1, which would be lower if the assumption of a single power-law between GeV and TeV energies was not made. This, however, would be an incorrect conclusion. The exclusion probabilities calculated in Table 2 are entirely independent of Method 1, using only the spectral break versus redshift distribution of Method 2. As such, what is truly driving the agreement or disagreement of a particular model with the data in this case is the near- to mid-IR ratio.

The minimum average slope between 2 and m, falling within our contour in Figure 9, is (). This value is in agreement with the limit placed by Aharonian et al. (2007d) of using 1ES 0229+200 data from H.E.S.S. Aharonian et al. (2007d) also calculate an upper limit on the EBL intensity at m of . This constraint is also compatible with our results, which places an upper bound at m of .

In summary, our analysis has shown the advantages of using these two complementary approaches in constraining the EBL. It has also provided evidence that the near- to mid-IR ratio of the EBL may be larger than previously thought. Future measurements with Fermi, as well as both current and next generation IACTs, such as the Cherenkov Telescope Array (CTA), will continue to improve these constraints.

FK acknowledges support from DOE. ED and FK acknowledge support from the Fermi Cycle-2 Guest Investigator program. The authors would also like to thank the anonymous referee for a thorough review of this manuscript which has improved both its clarity and completeness as well as revealed many subtleties regarding the analysis methods developed here.


  • Abdo et al. (2010a) Abdo, A. A., et al. 2010a, ApJ, 716, 30
  • Abdo et al. (2010b) —. 2010b, ApJ, 710, 1271
  • Acciari et al. (2010a) Acciari, V. A., et al. 2010a, ApJ, 709, L163
  • Acciari et al. (2010b) —. 2010b, ApJ, 715, L49
  • Aharonian et al. (2002) Aharonian, F., et al. 2002, A&A, 393, 89
  • Aharonian et al. (2005a) —. 2005a, A&A, 436, L17
  • Aharonian et al. (2005b) —. 2005b, A&A, 430, 865
  • Aharonian et al. (2006) —. 2006, Nature, 440, 1018
  • Aharonian et al. (2007a) —. 2007a, A&A, 470, 475
  • Aharonian et al. (2007b) —. 2007b, A&A, 473, L25
  • Aharonian et al. (2007c) —. 2007c, ApJ, 664, L71
  • Aharonian et al. (2007d) —. 2007d, A&A, 475, L9
  • Aharonian et al. (2008) —. 2008, A&A, 481, L103
  • Aharonian et al. (2010) —. 2010, A&A, 521, A69
  • Aharonian et al. (1999) Aharonian, F. A., et al. 1999, A&A, 349, 11
  • Albert et al. (2007a) Albert, J., et al. 2007a, ApJ, 662, 892
  • Albert et al. (2007b) —. 2007b, ApJ, 663, 125
  • Albert et al. (2007c) —. 2007c, ApJ, 669, 862
  • Albert et al. (2008) —. 2008, ApJ, 674, 1037
  • Aleksić et al. (2010) Aleksić, J., et al. 2010, A&A, 519, 32
  • Atwood et al. (2009) Atwood, W. B., et al. 2009, ApJ, 697, 1071
  • Bloom & Marscher (1996) Bloom, S. D., & Marscher, A. P. 1996, ApJ, 461, 657
  • Costamante et al. (2003) Costamante, L., Aharonian, F., Ghisellini, G., & Horns, D. 2003, NewA Rev., 47, 677
  • De Angelis et al. (2009) De Angelis, A., Mansutti, O., Persic, M., & Roncadelli, M. 2009, MNRAS, 394, L21
  • De Angelis et al. (2007) De Angelis, A., Roncadelli, M., & Mansutti, O. 2007, Phys. Rev. D, 76, 121301
  • Dermer & Schlickeiser (1993) Dermer, C. D., & Schlickeiser, R. 1993, ApJ, 416, 458
  • Domínguez et al. (2011) Domínguez, A., et al. 2011, MNRAS, 410, 2556
  • Dwek & Arendt (1998) Dwek, E., & Arendt, R. G. 1998, ApJ, 508, L9
  • Dwek & Krennrich (2005) Dwek, E., & Krennrich, F. 2005, ApJ, 618, 657
  • Dwek et al. (2005) Dwek, E., Krennrich, F., & Arendt, R. G. 2005, ApJ, 634, 155
  • Elbaz et al. (2002) Elbaz, D., Cesarsky, C. J., Chanial, P., Aussel, H., Franceschini, A., Fadda, D., & Chary, R. R. 2002, A&A, 384, 848
  • Essey et al. (2010a) Essey, W., Kalashev, O., Kusenko, A., & Beacom, J. F. 2010a, arXiv, 1011.6340v2
  • Essey et al. (2010b) Essey, W., Kalashev, O. E., Kusenko, A., & Beacom, J. F. 2010b, Phys. Rev. Lett., 104, 141102
  • Essey & Kusenko (2010) Essey, W., & Kusenko, A. 2010, Astropart. Phys., 33, 81
  • Fazio et al. (2004) Fazio, G. G., et al. 2004, ApJS, 154, 39
  • Franceschini et al. (2008) Franceschini, A., Rodighiero, G., & Vaccari, M. 2008, A&A, 487, 837
  • Gardner et al. (2000) Gardner, J. P., Brown, T. M., & Ferguson, H. C. 2000, ApJ, 542, L79
  • Georganopoulos et al. (2010) Georganopoulos, M., Finke, J. D., & Reyes, L. C. 2010, ApJ, 714, L157
  • Gilmore et al. (2009) Gilmore, R. C., Madau, P., Primack, J. R., Somerville, R. S., & Haardt, F. 2009, MNRAS, 399, 1694
  • Gould & Schréder (1967) Gould, R. J., & Schréder, G. P. 1967, Phys. Rev., 155, 1404
  • Hauser & Dwek (2001) Hauser, M. G., & Dwek, E. 2001, ARA&A, 39, 249
  • Hinton (2004) Hinton, J. A. 2004, NewA Rev., 48, 331
  • Holder et al. (2008) Holder, J., et al. 2008, in AIP Conf. Series, Vol. 1085, Proc. of the 4th International Meeting on High Energy Gamma-Ray Astronomy, 657
  • Hopwood et al. (2010) Hopwood, R., et al. 2010, ApJ, 716, L45
  • Horiuchi et al. (2009) Horiuchi, S., Beacom, J. F., & Dwek, E. 2009, Phys. Rev. D, 79, 83013
  • Imran & Krennrich (2008) Imran, A., & Krennrich, F. 2008, in Proc. of the 30th ICRC, Vol. 3, 981
  • Kashlinsky (2005) Kashlinsky, A. 2005, Phys. Rep., 409, 361
  • Katarzyński et al. (2006) Katarzyński, K., Ghisellini, G., Tavecchio, F., Gracia, J., & Maraschi, L. 2006, MNRAS, 368, L52
  • Keenan et al. (2010) Keenan, R. C., Barger, A. J., Cowie, L. L., & Wang, W.-H. 2010, ApJ, 723, 40
  • Kneiske et al. (2004) Kneiske, T. M., Bretz, T., Mannheim, K., & Hartmann, D. H. 2004, A&A, 413, 807
  • Kneiske et al. (2002) Kneiske, T. M., Mannheim, K., & Hartmann, D. H. 2002, A&A, 386, 1
  • Krennrich et al. (2008) Krennrich, F., Dwek, E., & Imran, A. 2008, ApJ, 689, L93
  • Krennrich et al. (2001) Krennrich, F., et al. 2001, ApJ, 560, L45
  • Krennrich et al. (2002) —. 2002, ApJ, 575, L9
  • Levenson & Wright (2008) Levenson, L. R., & Wright, E. L. 2008, ApJ, 683, 585
  • Madau & Phinney (1996) Madau, P., & Phinney, E. S. 1996, ApJ, 456, 124
  • Madau & Pozzetti (2000) Madau, P., & Pozzetti, L. 2000, MNRAS, 312, L9
  • Mankuzhiyil et al. (2010) Mankuzhiyil, N., Persic, M., & Tavecchio, F. 2010, ApJ, 715, L16
  • Maraschi et al. (1992) Maraschi, L., Ghisellini, G., & Celotti, A. 1992, ApJ, 397, L5
  • Mazin & Raue (2007) Mazin, D., & Raue, M. 2007, A&A, 471, 439
  • Papovich et al. (2004) Papovich, C., et al. 2004, ApJS, 154, 70
  • Petry et al. (2002) Petry, D., et al. 2002, ApJ, 580, 104
  • Prandini et al. (2010) Prandini, E., Bonnoli, G., Maraschi, L., Mariotti, M., & Tavecchio, F. 2010, MNRAS, 405, L76
  • Primack et al. (2005) Primack, J. R., Bullock, J. S., & Somerville, R. S. 2005, in AIP Conf. Proc., Vol. 745, 23
  • Primack et al. (2001) Primack, J. R., Somerville, R. S., Bullock, J. S., & Devriendt, J. E. G. 2001, in AIP Conf. Proc., Vol. 558, 463
  • Raue & Mazin (2008) Raue, M., & Mazin, D. 2008, IJMPD, 17, 1515
  • Samuelson et al. (1998) Samuelson, F. W., et al. 1998, ApJ, 501, L17
  • Sikora et al. (1994) Sikora, M., Begelman, M. C., & Rees, M. J. 1994, ApJ, 421, 153
  • Stecker et al. (2006) Stecker, F., Malkan, M., & Scully, S. 2006, ApJ, 648, 774
  • Tagliaferri et al. (2008) Tagliaferri, G., et al. 2008, ApJ, 679, 1029
  • Tavecchio et al. (2009) Tavecchio, F., Ghisellini, G., Ghirlanda, G., Costamante, L., & Franceschini, A. 2009, MNRAS, 399, L59

Appendix A Investigation of Systematics

To obtain a more thorough understanding of how the choice of certain parameter values (flux point error bars, number of iterations, redshift binning, etc.) affect this analysis, we have investigated the potential systematic effects introduced by each.

a.1 Test Blazar Spectrum Flux Point Errors

The application of 25% error bars to the test blazar spectrum (Section 3.2), while motivated by the mean percent error from VERITAS spectra, is a somewhat arbitrary choice. However, the systematic affect on the final result is limited, as illustrated in Figure 11. This figure shows the fit results for the TeV spectral break versus redshift distribution (Equation 8), after 100 iterations, assuming a particular EBL model. The left (right) panel shows the best fit slope (intercept) using six different size error bars, ranging from 15%–65%, for the test blazar spectrum. While the error in the best fit value increases with the error in the test blazar spectrum,131313This could be remedied by using a larger number of iterations (i.e., simulated observations). the standard deviation in the final result is of the error in the best fit determination from observations (Figure 5). From this one can be conclude that the systematic effect introduced by the choice of error bars in the test spectrum is minimal.

(a) Slope parameter.
(b) Intercept parameter.
Figure 11: Figures showing the slope (left) and intercept (right) fit results for the TeV spectral break versus redshift distribution (Equation 8) assuming a particular EBL model and using different size error bars for the test blazar spectrum (see Section 3.2). The six different choices of errors shown are 15%, 25%, 35%, 45%, 55%, and 65%. The standard deviations (unweighted) of the slope and intercept best fit parameters are 0.065 and 0.020, respectively. This spread in values is of the error in the best fit determination from observations (Figure 5).

a.2 Redshift Binning and the Number of Iterations for Determination

The choice of binning in redshift and the number of iterations used for the determination of are two potential sources of systematics. A redshift binning that is too course may result in a poor characterization of the spectral break versus redshift distribution. Many of the blazars used in the TeV spectral break analysis are separated in redshift from another source in the sample by (see Table 1). It may seem then that a redshift binning of 0.05 is insufficient for comparing the distributions for given EBL models with that from observation. Additionally, the choice of 100 iterations (i.e., simulated observations) for the determination of the distribution is not necessarily sufficient for the fit to converge to its true value.

Figure 12 addresses both these points. The four panels show the spectral break versus redshift distribution after 1, 25, 50, and 100 iterations using one particular EBL model. It can clearly be seen that the distribution has converged after 100 iterations. This is reflected not only in the data point error bars, but in the errors on the best fit parameters as well. It is also clear that, while the distribution is not perfectly linear, there is no evidence of features from small variations in redshift. In other words, a finer binning in redshift would have virtually no impact on the fit results.

(a) 1 iteration.
(b) 25 iterations.
(c) 50 iterations.
(d) 100 iterations.
Figure 12: distribution, using EBL absorbed test blazar spectra for a particular EBL scenario, calculated after 1, 25, 50, and 100 iterations.

a.3 Potential Pileup Effect at Very High Energies

The possibility exists that a redshift dependent spectral break at TeV could be the result of an effect intrinsic to the observing instrument. Perhaps the highest energy spectral bins contain a pileup of events that systematically hardens the spectrum. We have investigated this issue by examining two distributions. The first distribution of relevance is the spectral break (as shown in Figure 5) as a function of the spectral index obtained when fitting the blazars listed in Table 1 with a single power-law. This is plotted in Figure 13. A pileup effect would be most pronounced for soft sources due to fewer statistics at the highest energies. Consequently, if such an instrumental pileup were present, one would expect to see an increase in the measured spectral break as the spectral index increases (softens). The best linear fit to the data shown in 13 is given by , with a fit probability of 22%. This suggests that any pileup as a function of spectral index is approximately a effect.

The next important thing to establish is whether or not this effect could be responsible for the trend of increasing spectral break with redshift seen in Figure 5. This can be determined by looking at the single power-law spectral index distribution as a function of redshift (Figure 13). In order for a systematic pileup effect to be responsible for the spectral break versus redshift distribution observed, the spectral indices of blazars would have to be increasingly soft as the source redshift increases. As can be seen in Figure 13, this is not the case. The best linear fit to the spectral index versus redshift distribution is given by and is clearly consistent with being flat. This may seem surprising given that EBL absorption is expected to result in softer spectra as source redshifts increase. However, hard spectrum blazars are more readily detected by IACTs given that they generally have higher fluxes than their soft spectrum counterparts at very high energies. This results in the selection effect that blazars of higher redshift are more likely to be intrinsically hard.

One can conclude from these arguments that a pileup effect does not systematically produce a redshift dependent spectral break in IACT spectra. Furthermore, the spectra used in Method 2 come from three different IACTs all of which use different energy and spectral reconstruction techniques. This further reduces the likelihood that a redshift dependent spectral break could be purely an instrumental effect.

Figure 13: (a) Spectral break (as shown in Figure 5) plotted versus the spectral index obtained when fitting the blazar spectra with a single power-law. The best linear fit to the data is given by . (b) Single power-law spectral index plotted versus the blazar redshift. The best linear fit to the data is given by .
Comments 0
Request Comment
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
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description