THE EFFECT OF BLAZAR SPECTRAL BREAKS ON THE BLAZAR CONTRIBUTION TO THE EXTRAGALACTIC GAMMA-RAY BACKGROUND

The gamma-ray sky as currently observed by the Large Area Telescope (LAT) on board the Fermi Gamma-Ray Space Telescope consists of resolved point sources (such as active galactic nuclei (AGNs), pulsars, and star-forming galaxies), transient gamma-ray sources (e.g., gamma-ray bursts), and the diffuse gamma-ray radiation comprised of emission from the Galaxy and the extragalactic gamma-ray background (EGB). The origins of the EGB are, as yet, unknown; however, it is expected that unresolved, extragalactic point sources provide a sizable contribution to the EGB.

In its first year of taking data, Fermi observed 1451 resolved point sources with significance greater than 5σ (First Fermi-LAT catalog (1FGL); Abdo et al. 2010b), of which 573 were associated with blazars (AGNs for which the jet is closely aligned with the observer's line of sight (Blandford & Konigl 1979)). Thus, blazars constitute the largest class of astrophysical objects associated with gamma-ray sources, and it has long been suspected that unresolved blazars should provide a substantial contribution to the EGB (Padovani et al. 1993; Stecker et al. 1993; Salamon & Stecker 1994; Chiang et al. 1995; Stecker & Salamon 1996; Kazanas & Perlman 1997; Chiang & Mukherjee 1998; Sreekumar et al. 1998; Mukherjee & Chiang 1999; Mücke & Pohl 2000; Giommi et al. 2006; Narumoto & Totani 2006; Dermer 2007; Kneiske & Mannheim 2008; Inoue et al. 2008, 2010; Inoue & Totani 2009; Abdo et al. 2010c; Stecker & Venters 2011).

It should then, perhaps, come as no surprise that the Fermi-LAT, with its improved sensitivity enabling it to observe many more blazars, measured the integrated intensity of the EGB to be substantially lower than that measured by the Energetic Gamma-Ray Experiment Telescope (EGRET) aboard the Compton Gamma-Ray Observatory in the 1990s (1.03 (± 0.17) × 10⁻⁵ cm⁻² s⁻¹ sr⁻¹ for Fermi versus 1.45 (± 0.05) × 10⁻⁵ cm⁻² s⁻¹ sr⁻¹ for EGRET; Sreekumar et al. 1998; Abdo et al. 2010g). However, exactly how much of the change in the EGB integrated intensity is due to the ability to resolve many more blazars remains unclear. On the one hand, since the Sreekumar et al. (1998) determination of the EGRET EGB, models of the galactic foreground emission have been updated to reflect recent observations of the interstellar medium (Abdo et al. 2010g), accounting for at least some of the change. On the other hand, the determination of the distribution of blazars with respect to luminosity and redshift (the blazar gamma-ray luminosity function, GLF) and by extension, their contribution to the EGB directly from their observed flux distribution (such as that performed in Abdo et al. 2010c) is non-trivial, for four reasons.

First, blazars are variable at gamma-ray energies and without knowledge of the degree to which the gamma-ray luminosity of a blazar changes during the flaring state and the blazar duty cycle, the reconciliation of the observed blazar source counts with model GLFs is impossible (Stecker & Salamon 1996). This is because the increase in flux of a blazar (sometimes up to an order of magnitude; see, e.g., McLaughlin et al. 1996; Mukherjee et al. 1997; Nolan et al. 2003; Vercellone et al. 2004) during the flaring period introduces a selection effect in that flaring blazars are easier to detect than quiescent blazars, and quiescent blazars would have to be relatively bright and nearby to be detectable. The largest impact on the observed blazar source counts will be at fluxes just above the Fermi-LAT sensitivity as some blazars below the threshold make it into the sample because they are flaring. While the magnitude of the impact on the observed blazar source counts depends on the blazar GLF and the blazar duty cycle, the effect will be to flatten the faint-end slope of the observed blazar source counts (as blazars move from lower flux bins to higher ones) leading to an underestimation of the blazar contribution to the EGB as extrapolated from the observed source counts alone. Abdo et al. (2010c) concluded that since the peak-to-mean flux ratio for most Fermi sources is ∼ a factor of two, that there are no large systematic uncertainties due to blazar variability. However, a small peak-to-mean ratio is expected for blazars that are observed mostly in a flaring state, and as most of the observed blazars are fainter blazars that will be more represented by flaring blazars, then while it will be true that most blazars will have a small peak-to-mean ratio, it is not necessarily the case that the effect of blazar variability will be small.

Second, the large angular resolution of the Fermi-LAT at lower energies (∼5° at 100 MeV) could introduce significant source confusion for GLF models that predict high blazar densities resulting in fewer resolved blazars and a higher contribution to the EGB at lower energies (Stecker & Salamon 1999; Stecker & Venters 2011). The impact of source confusion on observed source counts depends on the source density as given by the blazar GLF. As such, accounting for source confusion requires a priori knowledge of the source density, which is exactly what the observer wishes to determine. Thus, source confusion further complicates analyses of the blazar contribution to the EGB based solely on observed source counts (for a detailed discussion of source confusion and the difference between the definition employed in this discussion and that employed by Abdo et al. 2010c, see Stecker & Venters 2011).

Third, as demonstrated in Abdo et al. (2010c), the reconstructed fluxes of the individual sources are subject to a considerable amount of uncertainty, particularly at the faint end (see their Figure 6). Finally, the model of the Fermi-LAT detection efficiency employed by Abdo et al. (2010c) assumed power-law spectra for the sources.⁵ As such, these uncertainties complicate the determination of the distribution of blazars with respect to luminosity and redshift (the GLF) and the determination of the blazar contribution to the EGB from their observed flux distribution. Notably, Stecker & Venters (2011) found that contrary to the conclusion of Abdo et al. (2010c), the observed flux distribution of Fermi blazars does not, as yet, rule out the possibility that the EGB is dominated by emission from unresolved blazars (see also Abazajian et al. 2011). Thus, the contribution of still unresolved blazars to the EGB remains in dispute.

Additional information about the contributions to the EGB can be obtained through studying the shape of its energy spectrum and the shapes of the collective intensity spectra of suspected contributors (Pavlidou & Venters 2008, hereafter PV08; Venters et al. 2009; Venters 2010). Analysis of the first year of Fermi data yielded an EGB spectrum consistent with a featureless power law with spectral index Γ ∼ 2.4 (Abdo et al. 2010g). Upon first reflection, the Fermi EGB spectrum appears to be consistent with the hypothesis that the EGB is dominated by emission from unresolved blazars since the mean spectral index for blazars, Γ₀, is also ∼2.4. However, several effects complicate this simple picture.

First, the astrophysical population of blazars is actually composed of two separate sub-populations (flat-spectrum radio quasars, FSRQs, and BL Lacertae type objects, BL Lacs) with distinct spectral properties (Γ₀ ∼ 2.45 for FSRQs and Γ₀ ∼ 2.2 for BL Lacs; Abdo et al. 2010c). Second, even within a given sub-population of blazars, the spectral indices of individual blazars form a distribution with some spread (Stecker & Salamon 1996; Venters & Pavlidou 2007), which causes the collective spectrum of unresolved blazars to curve as harder blazar spectra become more important at higher energies (Stecker & Salamon 1996; Pavlidou & Venters 2008). On the other hand, this effect is somewhat mitigated by the spectral bias introduced by the fact that blazars with harder spectral indices are easier to observe in a flux-limited survey⁶ and, thus, are more likely to be resolvable and not play as big a role in producing the EGB. Third, just as unresolved blazars are expected to contribute to the EGB, unresolved members of other known astrophysical gamma-ray emitters (such as star-forming and starburst galaxies) should also contribute, but with spectra that are substantially different from those of blazars (Pavlidou & Fields 2002; Fields et al. 2010; Lacki et al. 2011; Makiya et al. 2011; Stecker & Venters 2011) and may not even resemble power laws at gamma-ray energies. Finally, recent observations conducted by Fermi indicate that even blazar spectra can break at ∼ GeV energies and would no longer be describable by simple power laws (Abdo et al. 2009a; Inoue et al. 2008, 2010; Inoue & Totani 2009).

The intuitive conclusion would be that the combination of these effects should cause the energy spectrum of the EGB intensity to exhibit features (Pavlidou & Fields 2002). The observation instead of an EGB with a single, featureless power-law spectrum begs the question: is the lack of observed spectral features necessarily indicative of a single dominant source population with individual unbroken power-law spectra in the energy range of ⩾300 MeV? If not, then what can such a featureless power-law EGB spectrum tell us about the relative contributions of the separate populations and how do they depend on their collective spectra? The answer to these questions provide novel constraints on blazar GLFs and thus could have profound implications for the cosmological properties of known astrophysical gamma-ray emitters, their corresponding contributions to the EGB, and the general properties of their spectra at gamma-ray energies.

In this paper, we calculate the collective spectrum of unresolved blazars with individual spectra exhibiting broken power laws. In so doing, we seek to determine whether such a population of unresolved blazars can result in a featureless power-law spectrum resembling that of the EGB. If so, we could then investigate the implications for the relative contributions of FSRQs and BL Lacs, which could, in turn, have implications for the cosmological properties of blazars. In Section 2, we present the formalism of the calculation of the spectral shape of the collective unresolved blazar emission. In Section 3, we discuss the inputs of the calculation and their uncertainties. In Section 4, we present the results of the calculation, and we discuss them in Section 5.

To calculate the collective spectrum of unresolved blazars, we follow the formalism as outlined in PV08 with one major difference. Instead of taking blazar spectra to be simple power laws over gamma-ray energies (F∝E^−Γ, where Γ is the photon spectral index at gamma-ray energies), blazar spectra are taken to be smoothly broken power laws:

$$\begin{equation} F_E(E) = F_0\left[\left(\frac{E}{E_b}\right)^{\Gamma _1n}+\left(\frac{E}{E_b}\right)^{\Gamma _2n}\right]^{-1/n}, \end{equation} \tag{ 1 }$$

where F_E(E) is the differential photon flux in units of photons per unit area per unit energy per unit time, E_b is the break energy, Γ₁ is the low-energy slope, Γ₂ is the high-energy slope, and n quantifies the sharpness of the transition from the low-energy power law to the high-energy power law. For the purposes of this paper, we take n = 1. The total flux, F, of photons with energies greater than some fiducial energy, E_f, is found by integrating F_E(E) over energy,

$$\begin{equation} F = F_0\int _{E_f}^{\infty }\left[\left(\frac{E}{E_b}\right)^{\Gamma _1n}+\left(\frac{E}{E_b}\right)^{\Gamma _2n}\right]^{-1/n}{\it dE}. \end{equation} \tag{ 2 }$$

For the purposes of this paper, we take E_f = 100 MeV. Then, the contribution of a single unresolved blazar to the EGB is

$$\begin{equation} I_1 = \frac{F\left[(E/E_b)^{\Gamma _1n}+(E/E_b)^{\Gamma _2n}\right]^{-1/n}}{4\pi \int _{E_f}^{\infty }\left[(E/E_b)^{\Gamma _1n}+(E/E_b)^{\Gamma _2n}\right]^{-1/n}{\it dE}}, \end{equation} \tag{ 3 }$$

where the flux of one source is uniformly distributed over the entire sky in anticipation of an isotropically distributed cosmological population and I has units of photons per unit area per unit energy per unit time per unit solid angle.

Following PV08, we characterize the flux distribution of unresolved blazars as a function g(F) and the distribution of blazar spectral indices (or spectral index distribution, SID) as a function p(Γ), where Γ₁ = Γ − ΔΓ₁, Γ₂ = Γ + ΔΓ₂, and p(Γ), ΔΓ₁, and ΔΓ₂ are determined from observations (see Section 3). Then, the total contribution of unresolved blazars to the EGB is given by

$$\begin{equation} I(E) = \int _{-\infty }^{\infty }d\Gamma \int _{0}^{F_{\rm min}}dFg(F)I_1p(\Gamma), \end{equation} \tag{ 4 }$$

where F_min is the sensitivity of the gamma-ray telescope under consideration. For the first year of Fermi data, F_min ∼ 2 × 10⁻⁹ photons cm⁻² s⁻¹.

Equation (4) can be characterized in terms of factors that determine the overall magnitude of the unresolved blazar contribution to the EGB (flux terms) and factors that determine the overall shape of the blazar contribution (spectral index terms). For a carefully chosen definition of p(Γ) (determined from the analysis of spectral indices of the flux-limited sample of Fermi blazars that accounts for the spectral bias inherent in a flux-limited catalog; see Section 3 and Venters et al. 2009), the magnitude and shape terms decouple, and Equation (4) can be rewritten as

$$\begin{equation} I(E)\! =\! I_0\int _{-\infty }^{\infty }\! d\Gamma p(\Gamma)\frac{\left[(E/E_b)^{(\Gamma - \Delta \Gamma _1)n}+(E/E_b)^{(\Gamma + \Delta \Gamma _2)n}\right]^{-1/n}}{\mathcal {S}(E_f,\Gamma)}, \end{equation} \tag{ 5 }$$

where I₀ is a normalization constant depending on the flux distribution of unresolved blazars and

$$\begin{equation} \mathcal {S}(E_f,\Gamma) = \int _{E_f}^{\infty } {\it dE} \left[\left(\frac{E}{E_b}\right)^{(\Gamma - \Delta \Gamma _1)n}+\left(\frac{E}{E_b}\right)^{(\Gamma + \Delta \Gamma _2)n}\right]^{-1/n}. \end{equation} \tag{ 6 }$$

As demonstrated in PV08, the unresolved blazar contribution to the EGB is not just a question of magnitude, but also of spectral shape, and the spectral shape is sensitive to the distribution of blazar spectral indices at GeV energies. Both PV08 and Abdo et al. (2010c) assumed that blazar spectra at gamma-ray energies are power laws and did not account for possible breaks in the spectra. In this paper, blazar spectra take the forms of smoothly broken power laws (see Section 2).

For each sub-population of blazars, we determine the SID, p(Γ), from the likelihood analysis of Venters & Pavlidou (2007) fitting blazars from the Fermi-LAT First Catalog of AGNs (Abdo et al. 2010d) to Gaussian SIDs accounting for errors in measurement of individual spectral indices. Due to the survey bias toward harder spectral indices present in Fermi data, we applied the likelihood analysis only to the subset of blazars with photon fluxes ≳ 7 × 10⁻⁸ photons cm⁻² s⁻¹ and galactic latitudes ≳ 10° (as per Abdo et al. 2010c). We determined that the maximum-likelihood Gaussian SID can be characterized by a mean (Γ₀) and a spread (σ₀) with maximum-likelihood parameters determined to be Γ₀ = 2.45 and σ₀ = 0.16 for FSRQs and Γ₀ = 2.17 and σ₀ = 0.23 for BL Lacs, which are similar to the findings of Abdo et al. (2010c). Based on 1FGL spectra of two prominent blazars from the Fermi-LAT Bright AGN Sample (Abdo et al. 2009a; see Figure 1), we model the spectral breaks by taking ΔΓ₁ = 0.1 and ΔΓ₂ = 0.9. As in the measured blazar spectra, we treat the break energies distinctly for FSRQs and BL Lacs. We took the break energy, E_b, to be 4 GeV for FSRQs and 15 GeV for BL Lacs. In so doing, we consider two cases:

Scenario 1. Blazars within a population evolve such that their break energies are observed to be roughly the same.

Scenario 2. The break energies of the intrinsic spectra of blazars within a sub-population are the same, but because of redshift effects, they are observed to be different.

In order to perform the calculation in the Scenario 2, we separately model the contribution arising from different redshift bins from an assumed GLF model for a given sub-population of blazars and then determine the composite spectrum in the final step. For the purposes of this analysis, we assume the best-fit PLE (pure luminosity evolution) model of Narumoto & Totani (2006) for FSRQs and their best-fit LDDE (luminosity-dependent density evolution) model for BL Lacs.⁷ In this manner, we seek to investigate the impact of changing the GLF model(s) of the blazars. We should also note that alternative statistical treatments on more extensive data sets as in Abdo et al. (2010f) could result in different estimates of the break energies and changes in spectral index.

Another alternative treatment performed by Abdo et al. (2010c) consists of stacking the spectra of observed blazars. However, the analysis was performed on the flux-limited sample of Fermi blazars, which consists almost entirely of FSRQs (hence, the similarity between the stacked spectrum for the flux-limited sample of blazars and that of FSRQs). Moreover, while such an analysis is effective in determining the average spectral index of a given population of blazars, it reveals nothing about the spectral properties of the population of blazars (e.g., the spread in the SID, spectral breaks). Furthermore, as the Fermi-LAT survey is biased against sources with softer spectra, it is certainly biased against sources exhibiting significant breaks in their spectra. We also note that despite finding that the SIDs of blazars have non-zero spreads, Abdo et al. (2010c) did not account for these spreads in their calculated spectra of the blazar contributions to the EGB as evidenced by the fact that these spectra, as well as their estimated uncertainties, are still power laws (see Figures 18, 19, and 20 of Abdo et al. 2010c). Similarly, Abdo et al. (2010c) did not account for spectral breaks in their calculated spectra of the blazar contributions to the EGB. Inclusion of these effects would introduce curvature in the spectra of the blazar contributions to the EGB and not result in the power laws indicated by Abdo et al. (2010c).

As in PV08, we do not include information about the magnitude of the unresolved blazar contribution to the EGB. Instead, we normalize the collective spectrum of each sub-population in order to best fit (as determined using a χ² analysis) the measured spectrum of the EGB from the Fermi-LAT first year data (Abdo et al. 2010g). For the composite spectrum of FSRQs and BL Lacs, we normalize the collective spectrum of each sub-population such that their total composite spectrum is the best-fit spectrum to the measured spectrum.

We should note that while applying the aforementioned cuts to the sample accounts for much (though likely not all; see Venters et al. 2009) of the survey spectral bias, doing so also significantly reduces the sample size and will ultimately introduce more uncertainty in the determination of the likelihood parameters. Furthermore, while in the flux-limited sample of blazars, spectral index and flux do appear to decouple (see Figure 1 of Abdo et al. 2010c), we expect it to be only approximately correct. A more detailed calculation, however, requires a statistical analysis to determine the intrinsic blazar SID, which would, in turn, require knowledge of the blazar GLFs since the survey spectral bias depends on the redshift and luminosity distributions of blazars⁸ (Venters et al. 2009). In the case of FSRQs, applying the high-flux cut changes the SID very little (Abdo et al. 2010c; Stecker & Venters 2011), so we do not expect the collective spectrum of unresolved FSRQs to be appreciably different from that determined from the resolved FSRQs. On the other hand, in the case of BL Lacs, the sparseness of the population at fluxes above the cut introduces a considerable margin of error. As such, the collective spectrum of unresolved BL Lacs could be harder or softer than that expected from the resolved BL Lacs. As BL Lacs are, on average, harder than FSRQs, they will likely contribute most significantly at the higher energies, and hence, the differences in spectra between resolved and unresolved BL Lacs will mostly impact the collective unresolved blazar spectrum at higher energies.

Finally, we should note that recent multi-wavelength observations conducted by Fermi and other telescopes suggest that blazars can be further subdivided into the categories of low synchrotron peak, intermediate synchrotron peak, and high synchrotron peak (Abdo et al. 2010e). In principle, one could gain further insight by applying this procedure to each subdivision as there is some indication that these categories are also spectrally distinct at gamma-ray energies (Abdo et al. 2010e). Notably, as demonstrated in Abdo et al. (2010f), most high-synchrotron-peaked BL Lacs (HBLs) do not exhibit spectral breaks in the Fermi-LAT energy range while those HBL breaks that have been observed are different in character (going from soft to hard rather than hard to soft). However, in light of the cuts on galactic latitude and flux, further subdividing the sample of blazars could lead to small sample sizes in some of the categories, and consequently to poor SID parameter determination. Thus, for the purposes of this analysis, we retain the original classification scheme of FSRQs and BL Lacs. In any case, given that BL Lacs likely contribute most significantly at higher energies, the effect of the HBL spectra is likely to improve the fit at high energies, though a more detailed study is in order as more data become available. One might also be concerned about systematic changes in blazar spectral indices while flaring; however, analyses of EGRET blazar spectral indices found no evidence of such systematic changes in spectral index with flaring (Nandikotkur et al. 2007; Venters & Pavlidou 2007), and Fermi observations of individual blazars have thus far revealed no systematic changes in spectral index with time or flux (Abdo et al. 2009b, 2010a; Ackermann et al. 2010).

The spectral shapes of the contributions of unresolved blazars to the EGB as determined from the Gaussian SIDs and broken power-law gamma-ray spectra discussed in Section 3 are plotted in Figures 2–4. In Figure 2, we have plotted the best-fit spectral shapes of the collective intensities of unresolved FSRQs (χ²_red = 4.1) and BL Lacs (χ²_red = 3.3). In Figure 3, we have plotted the best-fit spectral shape of the total collective spectrum of unresolved blazars in Scenario 1 (solid line; χ²_red = 0.7; for explanation, see Section 3) and the individual contributions from FSRQs (dashed line) and BL Lacs (dot-dashed line). In Figure 4, we have plotted the best-fit spectral shape of the total collective spectrum of unresolved blazars in Scenario 2 (χ²_red = 0.6). For comparison, we have also plotted the best-fit spectral shape of the total collective spectrum assuming that FSRQ spectra break at E_b ∼ 1 GeV and BL Lac spectra break at E_b ∼ 5 GeV (assuming Scenario 1; Figure 5). Note that such spectra do not adequately fit the spectra in Figure 1. As such, this alternative is included merely for the sake of comparison.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As noted in Section 3, BL Lacs tend to be harder than FSRQs, and thus, the BL Lac collective spectrum is noticeably harder than that of the FSRQs. Also apparent is the effect of the higher break energy for the BL Lacs as compared with the FSRQs. As can be seen in Figure 2, the collective spectrum of neither sub-population of blazars reproduces well the spectrum of the EGB: the BL Lacs are too hard while the FSRQs are too soft. However, if the spectra are added together as in Figures 3 and 4 (renormalized as described in Section 3), the resulting spectrum (solid line) reproduces well that of the EGB (being within the statistical error bars of nearly all of the data points and within ∼1.5σ of the last data point). The harder spectra and higher break energies for the BL Lacs compensate for the suppressed intensity of FSRQs at higher energies. Notably, the redshift dependence of Scenario 2 shifts the breaks to lower energies resulting in the transition from FSRQ dominance to BL Lac dominance occurring at lower energies (∼1 GeV rather than ∼3 GeV). Thus, the relative contribution from BL Lacs to the total blazar collective spectrum is greater in Scenario 2 than in Scenario 1. In the case of the even lower break energies as shown in Figure 5, the total collective spectrum fits the EGB fairly well (χ²_red = 1.1), but the relative contributions of the FSRQs and BL Lacs are such that the contribution from BL Lacs dominates at all energies. This is in contrast with the scenarios presented in Figures 3 and 4 in which BL Lacs dominate only energies greater than ∼ few GeV. Also, the sharper contrast between Figures 3 and 5 than that between Figures 3 and 4 indicate that the transition and the resulting relative contributions of FSRQs and BL Lacs are more sensitive to the break energies than the model of the blazar GLF(s). This is due to the fact that each contribution is dominated by its closest and brightest (though still unresolved) members for which the observed break energies are similar to their intrinsic values for the GLF(s) considered.

We have calculated the spectral shapes of the contributions to the EGB arising from BL Lacs and FSRQs assuming blazar spectra can be described as broken power laws. We found that in the case that blazar spectral breaks are indeed common neither sub-population alone can adequately reproduce the spectrum of the EGB. However, in a combined spectrum, the harder spectra of the BL Lacs could compensate for the softer spectra of the FSRQs, resulting in a collective blazar spectrum that is similar to that of the EGB within uncertainties. Furthermore, we have found that the relative contributions of FSRQs and BL Lacs required to fit the EGB spectrum are sensitive to the nature of the spectral breaks and, to a lesser extent, the blazar GLFs.

The question of whether the relative contributions necessary to reproduce the overall spectrum of the EGB as determined in this method are reasonable given future Fermi measurements of the blazar GLFs (which account for the notable observational uncertainties discussed below and in Section 1) would provide insight into the overall blazar contribution to the EGB as well as the nature of blazar variability. According to current Fermi observations, the luminosities of FSRQs are greater than those of BL Lacs (L_{γ, FSRQ} ∼ 10² × L_{γ, BLL}; Abdo et al. 2009a, 2010d). On the other hand, BL Lacs are likely situated at lower, on average, redshifts than FSRQs (Dermer 2007). Thus, it is possible that the closer proximity of BL Lacs (or their numbers) compensates for their deficit in luminosity with respect to FSRQs resulting in their roughly comparable or slightly enhanced relative contribution to the collective blazar spectrum seen in Figures 3, 4, and 5. However, it should be noted that the closer proximity of BL Lacs and their harder spectral indices would also make them more easily observable by Fermi, limiting the contribution of unresolved BL Lacs to the EGB—though if the luminosity function for BL Lacs is broad enough, there could be many low-luminosity BL Lacs that would escape detection.

It thus remains to be seen how the collective intensities of the sub-populations of blazars actually compare with one another. In order for a true comparison to be drawn, the GLF for each sub-population of blazars needs to be measured. If it should be the case that the relative contributions are not reasonable given gamma-ray observations, then either the breaks included in the best-fit spectrum are not typical of blazars or blazars are not sufficient to explain the EGB. Already, Fermi studies of observed blazar flux distributions have suggested that emission from unresolved blazars may not comprise the dominant contribution to the EGB (Abdo et al. 2010c). However, given that the physics behind the gamma-ray emission of other known and speculated contributors (i.e., normal galaxies, cascades of ultra-high energy cosmic rays and TeV gamma rays, and dark matter annihilation) result in spectra that are quite distinct from that of the EGB as measured by Fermi (see, e.g., Ando et al. 2007; Kalashev et al. 2009; Siegal-Gaskins & Pavlidou 2009; Ahlers et al. 2010; Fields et al. 2010; Venters 2010; Berezinsky et al. 2011; Stecker & Venters 2011), a close resemblance of the collective spectrum of unresolved blazars could be striking.

The necessary reconciliation of the clues to the blazar contribution to the EGB provided by studies of the collective blazar spectrum with the clues provided by studies of the blazar source counts could thus provide insight in the gamma-ray emission properties of blazars. Blazar variability plays a substantial role as flaring blazars would be more easily observed by Fermi than quiescent blazars. As such, studies of the blazar GLF, observed blazar flux distributions, and ultimately, the blazar contribution to the EGB are largely a question of the blazar duty cycle. If blazars spend the majority of their lifetimes in the more observationally challenging quiescent state, then studies of the blazar contribution to the EGB based on observed flux distributions such as that presented in Abdo et al. (2010c) could underestimate the number of quiescent blazars resulting in an underestimation of the blazar contribution to the EGB. Thus, the apparent discrepancy of the predictions of the two analyses could be the result of the variability of gamma-ray emission in blazars, the study of which could have implications for gamma-ray emission in blazars.

A comparison of the two types of analyses could also provide insight in the cosmological properties of blazars. At 100 MeV, the Fermi-LAT angular resolution is ∼5° (Atwood et al. 2009); hence, for blazar GLFs that predict large number densities, many blazars, particularly those on the faint-end of the source count distributions, that are, in principle, observable would not be distinguishable from other blazars (and thus are unresolved) and would not be included in source count distributions. In such a scenario, the faint-end of the blazar flux distribution would underestimate the number of blazars in a given flux bin, and the faint-end slope might appear flatter than it should be,⁹ resulting in an underestimation of the blazar contribution to the EGB. Thus, the apparent discrepancy of the results of the two analyses could be the result of a large number density of blazars.

In conclusion, we have demonstrated that even with the inclusion of spectral breaks, the collective spectrum of unresolved blazars reproduces well the spectrum of the Fermi EGB for several break models. As such, we find that the possibility that the collective intensity of unresolved blazars dominates the EGB is not excluded on spectral grounds, even if spectral breaks are shown to be common among Fermi blazars. Given the remaining controversy concerning the blazar GLFs, we conclude that it is, as yet, premature to rule out blazar dominance of the EGB. Furthermore, we have shown that relative contributions of the sub-populations of blazars required to fit the EGB spectrum are sensitive to the nature of their breaks; hence, the methodology we present in this paper can be used to constrain the GLFs of blazars. As models for the GLFs of blazars (accounting for the uncertainties outlined in this paper) and more data on spectral breaks become available, the study of the spectral shape of the blazar contribution to the EGB in light of the breaks can provide insight into the high-energy physics of blazars and their cosmological properties.

We gratefully acknowledge enlightening discussions with Floyd Stecker. T.M.V. was supported by an appointment to the NASA Postdoctoral Program at the Goddard Space Flight Center, administered by Oak Ridge Associated Universities through a contract with NASA. V.P. acknowledges support for this work provided by NASA through Einstein Postdoctoral Fellowship grant number PF8-90060 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. This work was partially supported by NASA through the Fermi GI Program grant number NNX09AT74G.