© ESO, 2013

1. Introduction

The extragalactic background light (EBL) is a diffuse field of UV, optical, and infrared photons, with wavelengths in the range λ = 0.1 − 1000  μm, on which the integrated history of star formation in the Universe is imprinted. The spectral energy distribution (SED) of the EBL consists of two distinct components: the first, peaking in νF_ν around ≃ 1  μm and commonly referred to as the cosmic optical background (COB), was produced by thermal emission from stars since the big bang. The second component peaks at longer wavelengths (≃100  μm), has peak energy density comparable to the COB. It is referred to as the cosmic infrared background (CIB), and originates from the absorption and reemission of starlight by dust (see Hauser & Dwek 2001, for review).

Direct measurements of the EBL density are difficult due to local foregrounds, such as the zodiacal light and Galactic radiation (Hauser & Dwek 2001). They and are often interpreted as upper limits, while galaxy number counts in optical or infrared provide lower limits (Madau & Pozzetti 2000).

Since γ rays of observed energy E_γ can interact with EBL photons of energy E_EBL at a redshift z through γγ → e⁺ + e⁻ when E_γ/1   TeV > 0.26   eV/E_EBL(1 + z), the spectra of distant extragalactic sources measured in the very high-energy (VHE, E > 200 GeV) regime should differ from their emitted (intrinsic) spectra if the EBL density is nonzero¹. Since a large fraction of the emitted power in BL Lac-type blazars is in γ-ray band, this must be accounted for when SEDs are used to model their underlying physical properties (Coppi & Aharonian 1999).

Finding clear evidence for this EBL-induced attenuation has proven remarkably difficult to date. The fall-off in the EBL spectral density between the COB and CIB peaks (around 0.1   eV) should be visible as a kink in the measured VHE spectra around 1   TeV. This was sought, for example, in the blazar H 1426+428 by Aharonian et al. (2003), but results have been inconclusive given the large statistical errors. The signature of the EBL should also be evident in studies of the global population of VHE sources. This is because the energy-dependent attenuation increases with distance, such that the observed spectra are expected to become softer, i.e., the photon index, Γ, in power-law spectral fits should increase with redshift, z. Such studies have not been successful, with no evidence for a redshift-dependent effect being found by Mori (2009), De Angelis et al. (2009, the authors attributing this to varying spectral states inducing a large scatter in the data), or Orr et al. (2011, see in particular Fig. 13).

It has, however, been possible to constrain the EBL density in the energy range where they interact with observed γ rays. Using VHE spectra from distant BL Lac objects and a theoretically motivated conjecture that the photon index of the intrinsic spectrum cannot be harder than Γ_I ≃ 1.5, several authors have derived upper limits for the density close to the lower limits from galaxy counts (Dwek & Krennrich 2005; Aharonian et al. 2006). Recently, Ackermann et al. (2012) and Abramowski et al. (2013) measured the EBL density using its imprint in the spectra of BL Lac objects and found a density of EBL compatible with the best upper limits to date (Meyer et al. 2012).

Operating as an all-sky monitor, the Fermi LAT (Atwood et al. 2009) observes γ rays in the high energy (HE, 100   MeV < E < 300 GeV) range, where the effects of the EBL are much smaller than in the VHE. Sources detected both in the HE by Fermi and in the VHE by imaging atmospheric-Cherenkov telescopes (IACTs) then provide an opportunity to probe the effects of the energy-dependent attenuation from the EBL absorption across a much wider energy range. Here we present an updated list of GeV-TeV sources, building on the work of Abdo et al. (2009) and Ackermann et al. (2011).

2. Selection of the sources

Since the first detection of an extragalactic γ-ray source in the VHE range by Whipple (Punch et al. 1992), 50 active galactic nuclei (AGN) have been discovered in this energy band. The rate of detections has increased dramatically with the increased sensitivity of the latest generation of IACTs (VERITAS, MAGIC and H.E.S.S.) and an improved observation strategy using data from the Fermi LAT. An up-to-date view of the VHE sky can be found by browsing the TeVCat catalog² (Wakely & Horan 2008).

We selected AGN from TeVCat for which an HE and a VHE spectrum have been published and a firm redshift determined. Of the 50 AGN, 30 are BL Lacs with published spectral information. Three of these, namely, 3C 66A, PKS 1424+240, and PG 1553+113, do not have a firm redshift determination. In addition, there are two VHE-detected flat-spectrum radio quasars (FSRQs), 3C 279 and 4C +21.35. Despite their large redshift (see Table 1), internal absorption close to the emission region may strongly affect their spectra (Aharonian et al. 2008; Sitarek & Bednarek 2008); hence we include them in this study only for illustrative purposes. We also include the radio galaxies Centaurus A and M 87 in the sample.

The second Fermi catalog of AGN (2LAC, Ackermann et al. 2011) includes 1057 sources associated with AGN of many kinds. In this list, 36 out of the 37 VHE BL Lacs have a Fermi counterpart; the six that are not detected are SHBL J001355.9-185406, 1ES 0229+200, 1ES 0347-121, PKS 0548-322, HESS J1943+213, and 1ES 1312-423.

By merging the two lists, our sample contains 23 HE-VHE BL Lacs and two radio galaxies. The characteristics of these sources are listed in Table 1, including the photon indexes from published power-law spectral fits in the HE (from 2LAC) and VHE (from the reference given in the table), which form the dataset for this study.

AGN are observed to be variable in all wavelengths. In the VHE regime, flaring episodes have been observed from a number of such sources, in particular those that are bright and/or close by (such as Mrk 421, Mrk 501, and PKS 2155−304). However, most have not shown clear evidence of variability, many having been detected close to the sensitivity threshold of the instrument. It is not improbable that variability is a common feature of HBLs in the VHE regime, future instruments will be able to probe this in a larger population of fainter, more distant sources. In the HE range, BL Lacs, and especially high-frequency peaked BL Lacs (HBLs), are found to be the less variable (Abdo et al. 2010a). To reduce any bias introduced by the use of non-simultaneous observations, we used the VHE spectrum that has the lowest flux reported in the literature, which results in a generally good agreement between the overlapping energy ranges (Abdo et al. 2009).

3. Interpretation

3.1. Spectral evolution with the redshift

The mean HE and VHE indexes of our sample are ⟨ Γ_HE ⟩  = 1.86 and ⟨ Γ_VHE ⟩  = 3.18, respectively. For each source the photon index measured in the VHE range is greater than or compatible with that found in the HE, i.e., ΔΓ = Γ_VHE − Γ_HE ≳ 0. In the HE band, the standard deviation of the measured indexes is σ_HE = 0.26, and the excess variance, which accounts for the measurement errors³, is ${\mathit{\sigma }}_{\mathrm{HE}}^{\mathrm{XS}}\mathrm{=}\mathrm{0.24}$. In the VHE regime, the standard deviation is σ_VHE = 0.49, while the excess variance is ${\mathit{\sigma }}_{\mathrm{VHE}}^{\mathrm{XS}}\mathrm{=}\mathrm{0.10}$, showing that most of the sample variance can be ascribed to the errors on the individual measurements rather than to the intrinic distribution.

The points in Fig. 1 show ΔΓ versus the redshift z for our sample of sources. The two close-by radio galaxies do not show significant spectral breaks. For all other sources, ΔΓ ≳ 0.5, and those more distant than z = 0.1 exhibit a break of ΔΓ ≳ 1.0. A dependence of ΔΓ with z is apparent in our sample.

If we were to assume that the intrinsic spectrum of each object was well represented by a single power law across the entire HE and VHE domain, as seems to be the case for the two nearby radio galaxies for which ΔΓ ~ 0, we should expect that any significant break in the measured spectrum is the result of absorption on the EBL, i.e., ${\mathrm{\Gamma }}_{\mathrm{VHE}}\mathrm{=}{\mathrm{\Gamma }}_{\mathrm{Int}}\mathrm{+}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{EBL}}\mathrm{\left(}\mathit{E,z}\mathrm{\right)}\mathrm{\approx }{\mathrm{\Gamma }}_{\mathrm{Int}}\mathrm{+}\frac{\mathrm{d}\mathit{\tau }}{\mathrm{d}\log \mathit{E}}\mathrm{\left(}\mathit{E,z}\mathrm{\right)}$, where τ(E_γ,z) is the optical depth due to the attenuation by pair production (Abdo et al. 2009, Eq. (2)). To leading order, the EBL break increases linearly with the redshift, and we would expect ΔΓ to be correlated with z. The Pearson correlation factor for our dataset is ρ = 0.56 ± 0.11, more than 5σ away from 0, showing clear evidence that this dependency exists.

In Appendix B.2 we outline a simple method to evaluate the size of the break that might be expected from redshifting a curved intrinsic synchrotron self-Compton (SSC) spectrum within the HE and VHE observation windows, which gives rise to a K correction for the measured spectral indexes. This study shows that such a K correction would account for ~15% of the observed break for the most distant source in the sample at z = 0.3.

3.2. Expected EBL-induced spectral break

To further evaluate the data, we estimated the size of the spectral break as a function of redshift from the EBL density model of Franceschini et al. (2008, hereafter Fra08). We performed a simple simulation in which a hypothetical source with a flat spectrum (Γ_I = 0) is placed at a distance z and its flux attenuated by the EBL. The spectrum is simulated in 20 logarithmically spaced bins per decade. To evaluate the measured index a power law model if fitted to the bins above 200 GeV, with equal weighting per bin. The limited photon flux at the highest energies is accounted for in an ad hoc manner; the upper energy bound of the fit is chosen to be the point at which the differential flux is 1% of the flux at 200 GeV (up to a maximum of 10 TeV).

The predicted ΔΓ(z) is the black line shown in Fig. 1. The shaded gray areas show uncertainties on this calculation, which arise from:

• the ~10% energy resolution typical of IACTs, which is taken into account by shifting the energy bins by ± 10% (dark gray area),

• the threshold energy of the observations, which can vary from 100 GeV to more than 500 GeV (gray area), and

• the systematic error on the measured photon spectral index, typically 0.2 (light gray).

It is clear from Fig. 1 that, for the majority of sources, the observed break, ΔΓ, is systematically larger than that predicted by the EBL model, ΔΓ_EBL. This is most notably the case for 1ES 2344+514 (z = 0.044), PKS 2005-489 (z = 0.071), W Comae (z = 0.102), PKS 2155-304 (z = 0.116), and H 1426+428 (z = 0.129). The difference, ΔΓ − ΔΓ_EBL, is almost certainly the result of (convex) intrinsic curvature in the spectra of these objects, which is not unexpected (Perlman et al. 2005) and can have several interpretations, for instance, as being due to a turnover in the distribution of the underlying emitting particles (acceleration effects; e.g., Massaro et al. 2006) or to Klein-Nishina suppression (emission effects).

Nevertheless, it is striking that there are many sources for which ΔΓ ≃ ΔΓ_EBL. For these sources, the intrinsic broad-band γ-ray spectra are compatible (within errors) with single power laws. For those that additionally have Γ_HE ≲ 2.0, the high-energy peaks are not constrained by the current observations, even though they have a well-defined observational νF_ν peak. The most striking examples are H 2356-309 (z = 0.129), 1RXS J101015.9-311909 (z = 0.142), 1ES 1101-232 (z = 0.186), 1ES 0414+009 (z = 0.287), and S5 0716+714 (z = 0.310).

Fig. 1 Value of ΔΓ as a function of the redshift z. The black line is the theoretical break obtained with the Franceschini et al. (2008) model. The uncertainties due to the energy resolution of IACTs (dark gray), the different threshold energies (medium gray), and the systematic errors of 0.2 (light gray) are shown. The FSRQs and the BL Lacs with uncertain redshift are shown for illustration.

3.3. Constraining the EBL density

Since a significant fraction of the observed break can be directly attributed to the EBL, we attempted to constrain its density by applying a Bayesian model, which takes into account the effects discussed in the previous section. The model is described fully in Appendix A. We used two prescriptions to account for the effects of the EBL on the spectra, ΔΓ_EBL(E,z). In the first, the break is modeled as a linear function of the redshift, with coefficient a (as in Stecker & Scully 2006). We further assumed that the VHE measurements cover approximately the same energy range, so that the effect of energy threshold can be neglected: $\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{EBL}\hspace{0.17em}\mathit{i}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{EBL}}\mathrm{\left(}{\mathit{E}}_{\mathit{i}}\mathit{,}{\mathit{z}}_{\mathit{i}}\mathrm{\right|}\mathit{a}\mathrm{)}\mathrm{\approx }\mathit{a}{\mathit{z}}_{\mathit{i}}\mathrm{+}\mathit{O}\mathrm{(}{\mathit{z}}_{\mathit{i}}^{\mathrm{2}}\mathrm{)}\mathit{.}$(1)In the second model we applied a scaling factor, α, to the EBL model of Fra08, which resulted in an expected break of $\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{EBL}\hspace{0.17em}\mathit{i}}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}\mathrm{=}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{EBL}}\mathrm{\left(}{\mathit{E}}_{\mathit{i}}\mathit{,}{\mathit{z}}_{\mathit{i}}\mathrm{\right|}\mathit{\alpha }\mathrm{)}\mathrm{=}\mathit{\alpha }\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{Fra}\mathrm{08}}\mathrm{(}{\mathit{E}}_{\mathit{i}}\mathit{,}{\mathit{z}}_{\mathit{i}}\mathrm{)}\mathit{,}$(2)where ΔΓ_Fra08(E_i,z_i), was calculated for each source as in Sect. 3.2.

For both models, the posterior probability is computed and the results are given in Table 2. Figure 2 depicts the resulting ΔΓ for each model, using the mean value (i.e., ⟨ a ⟩ and ⟨ α ⟩) and the 95% confidence level (CL) lower limit (i.e., a_P < 95% and α_P < 95%).

The Bayesian model gives a value of ⟨ a ⟩  = 5.37 ± 0.65 and a 95% CL upper limit of 6.44 for the linear EBL model, which is significantly less than the value of 8.4 ± 1.0 reported by Yang & Wang (2010) using a simple χ² fit, which did not account for the intrinsic breaks. Stecker & Scully (2010) found that their baseline model can be approximated by a linear coefficient of 7.99, which cannot be reconciled with the results presented here. The null hypothesis, i.e., that there is no dependence of ΔΓ with the redshift (a = 0), is rejected at more than 8σ. The spectral break predicted using the model of Fra08 is in good agreement with the data; the mean scaling factor is ⟨ α ⟩  = 0.85 ± 0.10 and the 95% CL limit is α < 1.02.

Fig. 2 Value of ΔΓ as a function of the redshift z. Overlaid are the predicted breaks obtained from the Bayesian fit (simple line), as well as the 95% CL lower limits (hashed area). Left is the linear approximation (Eq. (1)) and right is the scaled model of Fra08 (Eq. (2)). For both models, none of the sources are significantly lower than the 95% CL lower limit. The symbols are descibed in Fig. 1.

4. Conclusions and perspectives

We have shown that broad-band γ-ray spectra (from ~ 100 MeV to a few TeV) carry the imprint of the EBL and provide a unique dataset to probe its properties. The redshift dependence of the difference of the photon indexes in VHE and HE range, ΔΓ, is found to be compatible with expectations from EBL attenuation. The Pearson correlation coefficient shows that ΔΓ and z are significantly correlated. We developed a Bayesian model to fit the dataset, accounting for intrinsic spectral softening, and find that the EBL density is consistent with the value predicted by Franceschini et al. (2008). Similar results have been found by Ackermann et al. (2012) and Abramowski et al. (2013), who modeled the EBL-absorbed spectra of AGNs detected in the HE and VHE regimes and found scaling factors of α_Fermi = 1.02 ± 0.23 and ${\mathit{\alpha }}_{\mathrm{HESS}}\mathrm{=}\mathrm{1.2}{\mathrm{7}}_{-0.15}^{\mathrm{+}\mathrm{0.18}}$, respectively. Their approach has the potential to be more powerful than that used here, since it can probe the features of the EBL-absorption signature as a function of energy. However, their approach is more reliant on the detailed modeling of the intrinsic spectra of the objects and of the EBL density and does not take advantage of the wider energy band available when the HE and VHE observations are combined. The two approaches are complementary and yield roughly compatible results within the combined statistical and systematic errors. Taking only the statistical errors, a χ² fit to the HESS, Fermi, and our results gives a mean value of α_combined ≈ 0.98 and a value of  χ² = 5.45 for two degrees of freedom, compatible with the hypothesis that the values are consistent at the 1.85σ level. Our model also offers a simple prescription for constraining the redshift of a GeV-TeV sources based on their measured value of ΔΓ (see Appendix A). Applying our findings to PG 1553+113 and 3C 66A leads to z < 0.64 and z < 0.55, respectively, in good agreement with the spectroscopic constraints.

No sources have breaks significantly smaller than ΔΓ_EBL. Significant deviations from the expected EBL-induced spectral breaks could indicate either that there is concave curvature in the intrinsic spectrum of the source or that other processes are at play during the propagation of γ rays, such as cosmic-ray interactions along the line of sight, which create spectral softening in high-redshift source spectra (Essey & Kusenko 2012). Our findings indicate that experimental uncertainties need to improve and firmer redshift estimations be established before the significance of this effect can be assessed.

The recent commissioning of the 28 m H.E.S.S. 2 telescope, the commissioning of the upgraded MAGIC telescopes, and the upgrade of the VERITAS cameras and trigger should increase the distance at which new blazars can be detected, while the planned CTA project (Actis et al. 2011) should reduce the uncertainties mentioned above due to its superior sensitivity. Finally, a useful feature of studies of spectral breaks vs. redshift, such as this one, is their capacity to provide distance estimations for BL Lacs (see, e.g., Abdo et al. 2010b; Prandini et al. 2010, 2012), since an estimated 50% of this population has an unknown or uncertain redshift.

Online material

Appendix A: Bayesian model

Appendix A.1: Development of the model

To extract information about the EBL density from the data presented in Fig. 1, we develop a “hierarchical” Bayesian model (using the terminology of Gelman et al. 2003), which is described by source-by-source spectral parameters and global parameters specifying, among other things, the EBL density of primary interest here. We use a Bayesian methodology to write the posterior density for the model parameters, marginalize over the source-by-source parameters that are not of interest, and produce estimates and confidence intervals for the EBL density.

In the development below, we make repeated use of the conditional probability rule (CPR) that the joint probability of two (sets of) events A and B, P(A,B), can be expressed as P(A,B) = P(A | B)P(B). If the two events are independent this becomes P(A,B) = P(A)P(B). We also frequently use the rule for marginalizing (or integrating) over unwanted parameters, P(A) = ^∫P(A,B)dB. Finally, we use the standard identity for the product of two Gaussians (see, e.g., Ahrendt 2005). In particular, if N(x | μ,σ²) denotes a Gaussian distribution with mean μ and variance σ², then $\begin{array}{ccc}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}\mathit{N}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{\mu }}_{\mathrm{1}}\mathit{,}{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{\left)}\mathit{N}\mathrm{\right(}\mathit{x}\mathrm{\left|}{\mathit{\mu }}_{\mathrm{2}}\mathit{,}{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{\right)}\mathrm{d}\mathit{x}\mathrm{=}& & \\ \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}\mathrm{(}{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}{\mathrm{)}}^{\frac{\mathrm{1}}{\mathrm{2}}}}\exp \begin{array}{c}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\mathrm{-}\frac{\mathrm{(}{\mathit{\mu }}_{\mathrm{1}}\mathrm{-}{\mathit{\mu }}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}\mathrm{(}{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}\mathrm{)}}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\mathit{.}& & \end{array}$(A.1)The dataset to be modeled consists of N GeV and TeV spectral measurements, ${\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{M}}$ and ${\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{M}}$ with measurement variances of ${\mathit{\sigma }}_{\mathrm{G}\mathit{i}}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{T}\mathit{i}}^{\mathrm{2}}$ and redshifts z_i (which are themselves not considered as measurement data). In what follows we refer frequently to the measured spectral break, $\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{M}}\mathrm{=}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{M}}\mathrm{-}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{M}}$. We write the dataset as $\mathit{Y}\mathrm{=}\mathrm{\left\{}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{M}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{M}}\mathrm{\right\}}$.

Each source is parametrized by four values, the intrinsic spectral indexes in the GeV and TeV regimes, ${\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}$ and ${\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}$, and the spectral indexes after absorption by the EBL, ${\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}$ and ${\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}$. In what follows we refer frequently to the intrinsic spectral break, $\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}\mathrm{=}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathrm{-}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}$. The global parameters that describe the EBL absorption itself are denoted abstractly as G and we write the set of all parameters of the model as $\mathrm{\Theta }\mathrm{=}\mathrm{\left\{}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}\mathit{,G}\mathrm{\right\}}$.

Bayes’ theorem allows us to write the posterior probability distribution of the parameters after the measurements have been made, P(Θ | Y), in terms of the standard likelihood of the data, P(Y | Θ), and the prior probability distribution of the model parameters, P(Θ): $\begin{array}{ccc}\mathit{P}\mathrm{\left(}\mathrm{\Theta }\mathrm{\right|}\mathit{Y}\mathrm{)}\mathrm{\propto }\mathit{P}\mathrm{(}\mathrm{\Theta }\mathrm{\left)}\mathit{P}\mathrm{\right(}\mathit{Y}\mathrm{\left|}\mathrm{\Theta }\mathrm{\right)}\mathit{.}& & \end{array}$The relation is written as a proportionality, instead of as an equality, because a global normalization factor has been neglected.

The model has four primary components that we discuss below. These are: (1) the likelihood for the GeV and TeV measurements, given the true absorbed spectral indexes of the sources; (2) a relationship between the absorbed index in the GeV regime and the intrinsic (unabsorbed) index; (3) a relationship for the analogous indexes in the TeV regime; and (4) a specification of how the intrinsic index in the GeV regime is related to the intrinsic index in the TeV regime.

We assume that the measurements of the individual indexes are independent (no correlation between measurements from different sources or between the GeV and TeV bands) and that the distribution for each measurement (${\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{M}}$ or ${\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{M}}$) is Gaussian with mean given by the appropriate absorbed index and with variance given by the measurement errors squared. Therefore, the likelihood is $\begin{array}{ccc}\mathit{P}\mathrm{\left(}\mathit{Y}\mathrm{\right|}\mathrm{\Theta }\mathrm{)}\mathrm{=}\underset{\mathit{i}}{\mathrm{􏽙}}\mathit{N}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{M}}\mathrm{|}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}\mathit{,}{{\mathit{\sigma }}_{\mathrm{G}\mathit{i}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{N}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{M}}\mathrm{|}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}\mathit{,}{{\mathit{\sigma }}_{\mathrm{T}\mathit{i}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{.}& & \end{array}$For the prior, we assume that the only link between the source-by-source index parameters for different sources comes through the global parameters G. Therefore, using the CPR we can write $\begin{array}{ccc}\mathit{P}\mathrm{\left(}\mathrm{\Theta }\mathrm{\right)}\mathrm{=}\mathit{P}\mathrm{\left(}\mathit{G}\mathrm{\right)}\underset{\mathit{i}}{\mathrm{􏽙}}\mathit{P}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}\mathrm{|}{\mathit{G}}^{\mathrm{)}}\mathit{.}& & \end{array}$It now remains only to describe how the intrinsic and absorbed indexes for each source are related. We assume that the absorbed index in each band depends only on the intrinsic index in that band and on the EBL parameters. Repeatedly applying the CPR, this gives $\begin{array}{ccc}\mathit{P}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}\mathrm{|}{\mathit{G}}^{\mathrm{)}}\mathrm{=}& & \\ \mathit{P}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}\mathrm{|}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}\mathit{,}{\mathit{G}}^{\mathrm{)}}\mathit{P}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}\mathrm{|}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathit{,}{\mathit{G}}^{\mathrm{)}}\mathit{P}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}\mathrm{|}{\mathit{G}}^{\mathrm{)}}\mathit{.}& & \end{array}$(A.2)We further assume that there is no absorption in the GeV regime and that the absorption in the TeV regime changes the index in a deterministic way. Specifically, we assume that (i) ${\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}\mathrm{=}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}$, and (ii) ${\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}\mathrm{=}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathrm{+}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{EBL}}\mathrm{\left(}{\mathit{E}}_{\mathit{i}}\mathit{,}{\mathit{z}}_{\mathit{i}}\mathrm{\right|}\mathit{G}\mathrm{)}$, where ΔΓ_EBL(E_i,z_i | G) is a function giving the change in TeV index for a source at redshift z_i measured at a TeV “threshold” energy of E_i. This can be expressed in terms of a probability using the Dirac δ-function: $\begin{array}{ccc}\mathit{P}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}\mathrm{|}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}\mathit{,}{\mathit{G}}^{\mathrm{)}}& \mathrm{=}& \mathit{\delta }\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}\mathrm{-}{{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}}^{\mathrm{)}}\mathrm{(}{\mathit{i}}^{\mathrm{)}}\\ \mathit{P}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}\mathrm{|}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathit{,}{\mathit{G}}^{\mathrm{)}}& \mathrm{=}& \mathit{\delta }\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}\mathrm{-}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathrm{-}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{EBL}}{\left({\mathit{E}}_{\mathit{i}}\mathit{,}{\mathit{z}}_{\mathit{i}}\mathrm{|}\mathit{G}\right)}^{\mathrm{)}}\mathrm{(}\mathit{i}{\mathit{i}}^{\mathrm{)}}\mathit{.}\end{array}$Equation (A.2) then reads $\begin{array}{ccc}\mathit{P}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}\mathrm{|}{\mathit{G}}^{\mathrm{)}}\mathrm{=}& & \\ \mathit{\delta }\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}\mathrm{-}{{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}}^{\mathrm{)}}\mathit{\delta }\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}\mathrm{-}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathrm{-}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{EBL}}{\left({\mathit{E}}_{\mathit{i}}\mathit{,}{\mathit{z}}_{\mathit{i}}\mathrm{|}\mathit{G}\right)}^{\mathrm{)}}\mathit{P}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}\mathrm{|}{\mathit{G}}^{\mathrm{)}}\mathit{.}& & \end{array}$The final and most interesting part of the model is to describe how the intrinsic indexes in the two bands are related. We expand this using the CPR to give $\begin{array}{ccc}\mathit{P}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}\mathrm{|}{\mathit{G}}^{\mathrm{)}}\mathrm{=}\mathit{P}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathrm{|}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}\mathit{,}{\mathit{G}}^{\mathrm{)}}\mathit{P}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}\mathrm{|}{\mathit{G}}^{\mathrm{)}}\mathit{.}& & \end{array}$We adopt a uniform prior for ${\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}$, $\mathit{P}\mathrm{\left(}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}\mathrm{\right)}\mathrm{=}\mathrm{1}$ and restrict ourselves to forms for the conditional probability that can be expressed as a function of the intrinsic break, $\mathit{P}\mathrm{\left(}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathrm{\right|}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}\mathit{,G}\mathrm{)}\mathrm{=}\mathit{P}\mathrm{(}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}\mathrm{\left|}\mathit{G}\mathrm{\right)}$. We assume that the break for each source is drawn from a single universal distribution that has, at most, some dependence on the global parameter set, G. The final expression for the prior for the parameters of each source is $\begin{array}{ccc}\mathit{P}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}\mathrm{|}{\mathit{G}}^{\mathrm{)}}\mathrm{=}& & \\ \mathit{\delta }\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}\mathrm{-}{{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}}^{\mathrm{)}}\mathit{\delta }\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}\mathrm{-}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathrm{-}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{EBL}}{\left({\mathit{E}}_{\mathit{i}}\mathit{,}{\mathit{z}}_{\mathit{i}}\mathrm{|}\mathit{G}\right)}^{\mathrm{)}}\mathit{P}\mathrm{(}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}\mathrm{|}{\mathit{G}}^{\mathrm{)}}\mathit{.}& & \end{array}$Putting everything together, recognizing that we have not yet discussed the parameters G, and leaving the exact choice of prior for the intrinsic break open for the present time, the full posterior density is $\begin{array}{ccc}\mathit{P}\mathrm{(}\mathrm{\left\{}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}\mathrm{\right\}}\mathit{,G}\mathrm{\left|}\mathrm{\right\{}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{M}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{M}}{\mathrm{\}}}^{\mathrm{)}}\mathrm{\propto }& & \\ \mathit{P}\mathrm{(}{\mathit{G}}^{\mathrm{)}}\underset{\mathit{i}}{\mathrm{􏽙}}\mathit{P}\mathrm{(}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}\mathrm{|}{\mathit{G}}^{\mathrm{)}}\mathit{N}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{M}}\mathrm{|}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}\mathit{,}{{\mathit{\sigma }}_{\mathrm{G}\mathit{i}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{N}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{M}}\mathrm{|}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}\mathit{,}{{\mathit{\sigma }}_{\mathrm{T}\mathit{i}}^{\mathrm{2}}}^{\mathrm{)}}& & \\ \mathrm{\times }\mathit{\delta }\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}\mathrm{-}{{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}}^{\mathrm{)}}\mathit{\delta }\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}\mathrm{-}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathrm{-}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{EBL}}{\left({\mathit{E}}_{\mathit{i}}\mathit{,}{\mathit{z}}_{\mathit{i}}\mathrm{|}\mathit{G}\right)}^{\mathrm{)}}\mathit{.}& & \end{array}$We marginalize over the parameters that are not of interest, ${\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}$, ${\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}$, ${\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}$, and ${\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}$ to give the posterior probability for the global parameters G. The integrals over the parameters for each source can be done separately: $\begin{array}{ccc}{\mathrm{ℐ}}_{\mathit{i}}\mathrm{=}\mathrm{\int }\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}\mathrm{\int }\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathrm{\int }\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}\mathrm{\int }\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}& & \\ \mathit{P}\mathrm{(}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}\mathrm{|}{\mathit{G}}^{\mathrm{)}}\mathit{N}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{M}}\mathrm{|}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}\mathit{,}{{\mathit{\sigma }}_{\mathrm{G}\mathit{i}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{N}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{M}}\mathrm{|}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}\mathit{,}{{\mathit{\sigma }}_{\mathrm{T}\mathit{i}}^{\mathrm{2}}}^{\mathrm{)}}& & \\ \mathrm{\times }\mathit{\delta }\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}\mathrm{-}{{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}}^{\mathrm{)}}\mathit{\delta }\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}\mathrm{-}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathrm{-}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{EBL}}{\left({\mathit{E}}_{\mathit{i}}\mathit{,}{\mathit{z}}_{\mathit{i}}\mathrm{|}\mathit{G}\right)}^{\mathrm{)}}\mathit{.}& & \end{array}$The integrals over ${\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{A}}$ and ${\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{A}}$ can be done immediately against the delta functions to give $\begin{array}{ccc}{\mathrm{ℐ}}_{\mathit{i}}\mathrm{=}\mathrm{\int }\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}\mathrm{\int }\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathit{P}\mathrm{(}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}\mathrm{|}{\mathit{G}}^{\mathrm{)}}& & \\ \mathrm{\times }\mathit{N}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{M}}\mathrm{|}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}\mathit{,}{{\mathit{\sigma }}_{\mathrm{G}\mathit{i}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{N}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{M}}\mathrm{|}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}\mathrm{+}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{EBL}}\left({\mathit{E}}_{\mathit{i}}\mathit{,}{\mathit{z}}_{\mathit{i}}\mathrm{|}\mathit{G}\right)\mathit{,}{{\mathit{\sigma }}_{\mathrm{T}\mathit{i}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{.}& & \end{array}$Making a change of integration variable from ${\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{I}}$ to $\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}$, the Gaussians can be manipulated to give $\begin{array}{ccc}{\mathrm{ℐ}}_{\mathit{i}}\mathrm{=}\mathrm{\int }\mathrm{d}\mathrm{(}\mathrm{\Delta }{{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}}^{\mathrm{)}}\mathit{P}\mathrm{(}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}\mathrm{|}{\mathit{G}}^{\mathrm{)}}\mathrm{\int }\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}& & \\ \mathrm{\times }\mathit{N}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}\mathrm{|}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{M}}\mathit{,}{{\mathit{\sigma }}_{\mathrm{G}\mathit{i}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{N}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{I}}\mathrm{|}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{M}}\mathrm{-}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}\mathrm{-}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{EBL}}\left({\mathit{E}}_{\mathit{i}}\mathit{,}{\mathit{z}}_{\mathit{i}}\mathrm{|}\mathit{G}\right)\mathit{,}{{\mathit{\sigma }}_{\mathrm{T}\mathit{i}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{.}& & \end{array}$The second integral can be evaluated using Eq. (A.1). Putting all the source integrals together gives the general expression $\begin{array}{ccc}\mathit{P}\mathrm{(}\mathit{G}\mathrm{\left|}\mathrm{\right\{}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{M}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{M}}{\mathrm{\}}}^{\mathrm{)}}\mathrm{\propto }\mathit{P}\mathrm{(}{\mathit{G}}^{\mathrm{)}}& & \\ \underset{\mathit{i}}{\mathrm{􏽙}}\mathrm{\int }\mathrm{d}\mathrm{(}\mathrm{\Delta }{{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}}^{\mathrm{)}}\mathit{P}\mathrm{(}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}\mathrm{|}{\mathit{G}}^{\mathrm{)}}\mathit{N}\mathrm{(}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}\mathrm{|}\mathrm{\Delta }{{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{M}}}^{}& & \\ \mathrm{-}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{EBL}}\left({\mathit{E}}_{\mathit{i}}\mathit{,}{\mathit{z}}_{\mathit{i}}\mathrm{|}\mathit{G}\right)\mathit{,}{\mathit{\sigma }}_{\mathrm{G}\mathit{i}}^{\mathrm{2}}\mathrm{+}{{\mathit{\sigma }}_{\mathrm{T}\mathit{i}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{.}& & \end{array}$(A.3)

Appendix A.2: Various priors for $\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}$

We examine three concrete cases for $\mathit{P}\mathrm{\left(}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}\mathrm{\right|}\mathit{G}\mathrm{)}$. The first two are based on simple assumptions and result in analytic expressions for the full posterior probability that are easy to understand. The final prior distribution is derived from Monte Carlo realizations of an SSC model, as described in Appendix B.1. The results from this final case are presented in Sect. 3.3. Here we develop the first two cases. A relatively weak assumption is that the TeV index can be no harder than the GeV index. It would seem reasonable to express this using the Heaviside step function, Θ(x), as $\mathit{P}\mathrm{\left(}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}\mathrm{\right|}{\mathit{G}}^{\mathrm{)}}\mathrm{=}\mathrm{\Theta }\mathrm{(}\mathrm{\Delta }{{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}}^{\mathrm{)}}\mathit{.}$(A.4)However, this is unsatisfactory, as it asserts that the mean intrinsic break is infinite, $\mathrm{⟨}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}\mathrm{⟩}\mathrm{\to }\mathrm{\infty }$, which is clearly not realistic. As seen below, this results in an unphysical posterior distribution. To remedy this failing, we instead assume that the prior distribution of the intrinsic break is given by a Gaussian, truncated at negative values: $\mathit{P}\mathrm{\left(}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}\mathrm{\right|}{\mathit{G}}^{\mathrm{)}}\mathrm{=}\mathrm{\Theta }\mathrm{(}\mathrm{\Delta }{{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}}^{\mathrm{)}}\mathit{N}\mathrm{\left(}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}\mathrm{\right|}{\mathit{\mu }}_{\mathrm{I}}\mathit{,}{{\mathit{\sigma }}_{\mathrm{I}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{.}$(A.5)In this case the prior is parametrized by two values, μ_I and σ_I, which must be either estimated in the problem (i.e., added to the global parameter set G) or specified externally. We will simply assume μ_I = 0 and derive the results for various reasonable values of σ_I.

Equation (A.3) can be used to calculate the posterior distribution in the three cases for $\mathit{P}\mathrm{\left(}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}\mathrm{\right)}$ discussed above. In the second case, with the prior given by Eq. (A.5), the final integration can be done to give, after applying the formula for the product of Gaussians, $\begin{array}{ccc}\mathit{P}\mathrm{(}\mathit{G}\mathrm{\left|}\mathrm{\right\{}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{M}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{M}}\mathrm{\}}\mathit{,}{\mathit{\mu }}_{\mathrm{I}}\mathrm{=}\mathrm{0}\mathit{,}{\mathit{\sigma }}_{\mathrm{I}}\mathrm{)}\mathrm{\propto }& & \\ \mathit{P}\mathrm{(}{\mathit{G}}^{\mathrm{)}}\underset{\mathit{i}}{\mathrm{􏽙}}\exp \begin{array}{c}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\mathrm{-}\frac{{\mathrm{(}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{EBL}}\left({\mathit{E}}_{\mathit{i}}\mathit{,}{\mathit{z}}_{\mathit{i}}\mathrm{|}\mathit{G}\right)\mathrm{-}\mathrm{\Delta }{{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{M}}}^{\mathrm{)}}}^{\mathrm{2}}}{\mathrm{2}\mathrm{(}{\mathit{\sigma }}_{{\mathrm{G}}_{\mathrm{i}}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{{\mathrm{T}}_{\mathrm{i}}}^{\mathrm{2}}\mathrm{+}{{\mathit{\sigma }}_{\mathrm{I}}^{\mathrm{2}}}^{\mathrm{)}}}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}& & \\ \mathrm{\times }\mathrm{erfc}\begin{array}{c}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\frac{{\mathit{\sigma }}_{\mathrm{I}}\mathrm{(}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{EBL}}\left({\mathit{E}}_{\mathit{i}}\mathit{,}{\mathit{z}}_{\mathit{i}}\mathrm{|}\mathit{G}\right)\mathrm{-}\mathrm{\Delta }{{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{M}}}^{\mathrm{)}}}{\sqrt{\mathrm{2}}{\mathrm{(}{\mathit{\sigma }}_{{\mathrm{G}}_{\mathrm{i}}}^{\mathrm{2}}\mathrm{+}{{\mathit{\sigma }}_{{\mathrm{T}}_{\mathrm{i}}}^{\mathrm{2}}}^{\mathrm{)}}}^{\frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{(}{\mathit{\sigma }}_{{\mathrm{G}}_{\mathrm{i}}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{{\mathrm{T}}_{\mathrm{i}}}^{\mathrm{2}}\mathrm{+}{{\mathit{\sigma }}_{\mathrm{I}}^{\mathrm{2}}}^{\mathrm{)}}}^{\frac{\mathrm{1}}{\mathrm{2}}}}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\mathit{,}& & \end{array}$(A.6)where $\mathrm{erfc}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}}{\sqrt{\mathit{\pi }}}{{}^{\mathrm{\int }}}_{\mathit{x}}^{\mathrm{\infty }}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}{\mathit{x}}^{\mathrm{\prime }\mathrm{2}}}$ is the complementary error function.

It is instructive to examine the two limiting cases of σ_I = 0 and σ_I → ∞. In the first case, we arrive at $\begin{array}{ccc}\log \mathit{P}\mathrm{(}\mathit{G}\mathrm{\left|}\mathrm{\right\{}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{M}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{M}}\mathrm{\}}\mathit{,}{\mathit{\mu }}_{\mathrm{I}}\mathrm{=}\mathrm{0}\mathit{,}{\mathit{\sigma }}_{\mathrm{I}}\mathrm{=}{\mathrm{0}}^{\mathrm{)}}\mathrm{\propto }& & \\ \sum _{\mathit{i}}\frac{{\mathrm{(}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{EBL}}\left({\mathit{E}}_{\mathit{i}}\mathit{,}{\mathit{z}}_{\mathit{i}}\mathrm{|}\mathit{G}\right)\mathrm{-}\mathrm{\Delta }{{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{M}}}^{\mathrm{)}}}^{\mathrm{2}}}{{\mathit{\sigma }}_{{\mathrm{G}}_{\mathrm{i}}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{{\mathrm{T}}_{\mathrm{i}}}^{\mathrm{2}}}\mathit{,}& & \end{array}$which is exactly the expression that would result from a simple least-squares fit to the measured spectral breaks. In the second case (σ_I → ∞), we have $\begin{array}{ccc}\mathit{P}\mathrm{(}\mathit{G}\mathrm{\left|}\mathrm{\right\{}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{M}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{M}}\mathrm{\}}\mathit{,}{\mathit{\mu }}_{\mathrm{I}}\mathrm{=}\mathrm{0}\mathit{,}{\mathit{\sigma }}_{\mathrm{I}}\mathrm{\to }{\mathrm{\infty }}^{\mathrm{)}}\mathrm{\propto }& & \\ \underset{\mathit{i}}{\mathrm{􏽙}}\mathrm{erfc}\begin{array}{c}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\frac{\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{EBL}}\left({\mathit{E}}_{\mathit{i}}\mathit{,}{\mathit{z}}_{\mathit{i}}\mathrm{|}\mathit{G}\right)\mathrm{-}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{M}}}{\sqrt{\mathrm{2}}{\mathrm{(}{\mathit{\sigma }}_{\mathrm{G}\mathit{i}}^{\mathrm{2}}\mathrm{+}{{\mathit{\sigma }}_{\mathrm{T}\mathit{i}}^{\mathrm{2}}}^{\mathrm{)}}}^{\frac{\mathrm{1}}{\mathrm{2}}}}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\mathit{,}& & \end{array}$(A.7)which is the same expression as would be derived starting from the prior given by Eq. (A.4). We therefore have an expression that transforms continuously between the two clearly identifiable extremities as a function of σ_I.

As described in Sect. 3.3, we attempt to derive constraints on two simple models for the EBL: (1) a linear approximation, given by Eq. (1), and (2) a scaling of the EBL model of Franceschini et al. (2008), as described by Eq. (2). In both cases, we assume a uniform prior in the positive region of the parameter space for a, P(a) = Θ(a) and α, P(α) = Θ(α), respectively.

Using either of these, the deficiencies in the model given by Eq. (A.4) is finally evident. Combining Eqs. (2) and (A.7), we get $\begin{array}{ccc}\mathit{P}\mathrm{(}\mathit{\alpha }\mathrm{\left|}\mathrm{\right\{}{\mathrm{\Gamma }}_{\mathrm{G}\mathit{i}}^{\mathrm{M}}\mathit{,}{\mathrm{\Gamma }}_{\mathrm{T}\mathit{i}}^{\mathrm{M}}\mathrm{\}}\mathit{,}{\mathit{\mu }}_{\mathrm{I}}\mathrm{=}\mathrm{0}\mathit{,}{\mathit{\sigma }}_{\mathrm{I}}\mathrm{\to }{\mathrm{\infty }}^{\mathrm{)}}\mathrm{\propto }& & \\ \mathrm{\Theta }\mathrm{(}{\mathit{\alpha }}^{\mathrm{)}}\underset{\mathit{i}}{\mathrm{􏽙}}\mathrm{erfc}\begin{array}{c}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\frac{\mathit{\alpha }\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{Fra}\mathrm{08}}\mathrm{(}{\mathit{E}}_{\mathit{i}}\mathit{,}{\mathit{z}}_{\mathit{i}}\mathrm{)}\mathrm{-}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{M}}}{\sqrt{\mathrm{2}}{\mathrm{(}{\mathit{\sigma }}_{\mathrm{G}\mathit{i}}^{\mathrm{2}}\mathrm{+}{{\mathit{\sigma }}_{\mathrm{T}\mathit{i}}^{\mathrm{2}}}^{\mathrm{)}}}^{\frac{\mathrm{1}}{\mathrm{2}}}}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\mathit{.}& & \end{array}$(A.8)The most probable value occurs at α = 0, which is consistent with the assertion in this case that $\mathrm{⟨}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{I}}\mathrm{⟩}\mathrm{\to }\mathrm{\infty }$.

The results in Sect. 3.3 are derived from combining Eqs. (1) (or (2)) and (A.3) with the prior calculated in Appendix B.1 and integrating numerically.

Appendix A.3: Results with different priors

Table A.1 presents the results obtained with the half-Gaussian prior of Eq. (A.5) (with μ_I = 0) for four values for the variance of the distribution of the intrinsic break, σ_I = 0.25,0.5,1.0,  and  2.0, and using the prior derived from SSC modeling, illustrated in Fig. B.1.

With increasing value of σ_I, the most probable value of a or α decreases, since a lower EBL-induced break is needed to reproduce the data. Results derived with values of σ_I = 0.5  and  1.0 are in good agreement with those from the SSC model.

Appendix A.4: Constraints on the redshift

The Bayesian methodology can also be used to constrain the redshifts of GeV-TeV blazars from their measured spectral breaks. In the case of a single source, Eqs. (A.3) and (2) can be adopted to express the posterior probability for the parameters G =  { α,z }, given their priors P(α) and P(z)⁴. This can then be marginalized over α to give $\begin{array}{ccc}\mathit{P}\mathrm{(}\mathit{z}\mathrm{|}{\mathrm{\Gamma }}_{\mathrm{G}}^{\mathrm{M}}\mathit{,}{{\mathrm{\Gamma }}_{\mathrm{T}}^{\mathrm{M}}}^{\mathrm{)}}\mathrm{\propto }\mathit{P}\mathrm{(}{\mathit{z}}^{\mathrm{)}}\mathrm{\int }\mathrm{d}\mathit{\alpha P}\mathrm{(}{\mathit{\alpha }}^{\mathrm{)}}\mathrm{\int }\mathrm{d}\mathrm{(}\mathrm{\Delta }{\mathrm{\Gamma }}^{\mathrm{I}}\mathrm{)}\mathit{P}\mathrm{(}\mathrm{\Delta }{\mathrm{\Gamma }}^{\mathrm{I}}\mathrm{)}\mathrm{\times }& & \\ \mathit{N}\mathrm{(}\mathrm{\Delta }{\mathrm{\Gamma }}^{\mathrm{I}}\mathrm{|}\mathrm{\Delta }{\mathrm{\Gamma }}^{\mathrm{M}}\mathrm{-}\mathit{\alpha }\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{Fra}\mathrm{08}}\left(\mathit{E,z}\mathrm{|}\mathit{G}\right)\mathit{,}{\mathit{\sigma }}_{\mathrm{G}}^{\mathrm{2}}\mathrm{+}{{\mathit{\sigma }}_{\mathrm{T}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{.}& & \end{array}$(A.9)The prior for the redshift, P(z), could be estimated from the redshift distribution of detected GeV-TeV detected blazars, i.e., from the values presented in Table 1. However, the true distribution is probably not well represented by the small number of sources detected, so we seek an alternative approach. Another option is to use a flat distribution, which is conservative but also unrealistic. We instead compromise and use the distribution of 2FGL BL Lac objects and unknown AGN. Since Fermi AGN are detected to higher redshifts than those at TeV energies, this should still be a conservative approach. For the prior on the EBL scaling, P(α), we use the positive portion of a Gaussian with mean 1.0 and RMS 0.3, P(α) = Θ(α)N(α | 1.0,0.3²). Finally, we use the prior on ΔΓ^I derived from our SSC simulations, as described above.

The mean redshift and upper limits for some values of ΔΓ, calculated using a typical value of ${\mathit{\sigma }}_{\mathrm{G}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{T}}^{\mathrm{2}}\mathrm{=}{\mathrm{0.2}}^{\mathrm{2}}$, are given in Table A.2. The corresponding relation between ⟨ z ⟩ and ΔΓ can be approximated by $\begin{array}{ccc}\mathrm{⟨}\mathit{z}\mathrm{⟩}\mathrm{\approx }\mathrm{0.024}\mathrm{+}\mathrm{0.079}\mathrm{\Delta \Gamma }\mathrm{+}\mathrm{0.022}\mathrm{\Delta }{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{-}\mathrm{0.0010}\mathrm{\Delta }{\mathrm{\Gamma }}^{\mathrm{3}}& & \end{array}$and the relation between z_P < 95% and ΔΓ by $\begin{array}{ccc}{\mathit{z}}_{\mathit{P}\mathit{<}\mathrm{95}\mathrm{\%}}\mathrm{\approx }\mathrm{0.081}\mathrm{+}\mathrm{0.081}\mathrm{\Delta \Gamma }\mathrm{+}\mathrm{0.080}\mathrm{\Delta }{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{-}\mathrm{0.011}\mathrm{\Delta }{\mathrm{\Gamma }}^{\mathrm{3}}\mathit{.}& & \end{array}$

Appendix B: Synchrotron self-Compton simulations

Appendix B.1: Determination of the intrinsic break properties

Fig. B.1 ΓVHE − ΓHE distribution (red histogram) obtained with the SSC simulations and application of cuts described in the text.

Fig. B.2 Value of ΔΓ as a function of the redshift z. The gray line is the break due to redshift effect as computed with the SSC simulation.

In order to derive a plausible prior probability density for the intrinsic break between HE and VHE for use in the Bayesian model, we produce a set of Monte Carlo simulations of hypothetical BL Lacs using a one-zone SSC model (Band & Grindlay 1985), which is often used to successfully reproduce the time-averaged SEDs from radio to TeV energies.

We assume a spherical emission zone, with a size R, moving at a bulk Doppler factor δ. This region is filled by a uniform magnetic field B, and a population of electrons with a density N_e(γ) is responsible for the synchrotron emission. The synchrotron photons are upscattered by the same population of electrons to produce γ rays.

The distribution of electrons is described by a power law with an exponential cutoff (Lefa et al. 2012), N_e(γ) = N₀γ^p·exp( − γ/γ_cut). The model therefore has three parameters to describe the electron population (N₀, p, and γ_cut) and four to describe the jet properties (z, R, δ, and B). Among these parameters, R, δ, and N₀ have only an achromatic effect, and z produces a small K correction, which is evaluated separately in Appendix B.2. The spectral break is also insensitive to the value of the index p of the electron distribution. The parameters that determine the intrinsic break are therefore B and γ_cut.

We perform 10⁵ simulations in which the values of B and γ_cut are uniformly drawn in the range 0.01 < B < 0.5 G and 3 × 10⁴ < γ_cut < 1 × 10⁷. The other parameters are kept fixed at the values given in Table B.1, which are typical for BL Lacs.

Since only the BL Lac-type SEDs are of interest in this study, simulations have been removed for which the synchrotron emission peaks at energies lower than 10 eV or higher that 10 MeV and which have an HE index of Γ > 2. The distribution of the intrinsic break ${\mathrm{\Gamma }}_{\mathrm{VHE}}^{\mathrm{I}}\mathrm{-}{\mathrm{\Gamma }}_{\mathrm{HE}}^{\mathrm{I}}$, depicted in Figure B.1, has a sharp rise below 0.2 and a long tail at higher break values.

Appendix B.2: Effects of the energy shift due to the distance

The SSC simulations can be used to evaluate the size of the K correction required to account for the redshifting of the intrinsic spectrum into the fixed HE and VHE observation windows. This effect produces a trend of increasing observed ΔΓ with z, even in the absence of EBL absorption.

As before, the SSC parameters are fixed to the values in Table B.1, but with B = 0.1   G and γ_cut = 1.6 × 10⁵. This produces an intrinsic spectrum with ${\mathrm{\Gamma }}_{\mathrm{HE}}^{\mathrm{I}}\mathrm{=}\mathrm{1.85}$ and ${\mathrm{\Gamma }}_{\mathrm{VHE}}^{\mathrm{I}}\mathrm{=}\mathrm{3.04}$, i.e., ΔΓ_I = 1.19. Simulating the same source with redshifts in the range 10^-4 < z < 0.7 and with no EBL absorption leads to an additional component in the observed ΔΓ, which increases with redshift, as depicted in Fig. B.2. This corresponds to the K correction for this intrinsic spectral shape and is too small to explain the trend observed in the data.

References

All Tables

All Figures

	Fig. 1 Value of ΔΓ as a function of the redshift z. The black line is the theoretical break obtained with the Franceschini et al. (2008) model. The uncertainties due to the energy resolution of IACTs (dark gray), the different threshold energies (medium gray), and the systematic errors of 0.2 (light gray) are shown. The FSRQs and the BL Lacs with uncertain redshift are shown for illustration.
In the text

Fig. 2 Value of ΔΓ as a function of the redshift z. Overlaid are the predicted breaks obtained from the Bayesian fit (simple line), as well as the 95% CL lower limits (hashed area). Left is the linear approximation (Eq. (1)) and right is the scaled model of Fra08 (Eq. (2)). For both models, none of the sources are significantly lower than the 95% CL lower limit. The symbols are descibed in Fig. 1.

In the text

	Fig. B.1 ΓVHE − ΓHE distribution (red histogram) obtained with the SSC simulations and application of cuts described in the text.
In the text

	Fig. B.2 Value of ΔΓ as a function of the redshift z. The gray line is the break due to redshift effect as computed with the SSC simulation.
In the text