THE COSMOLOGICAL IMPACT OF LUMINOUS TeV BLAZARS. II. REWRITING THE THERMAL HISTORY OF THE INTERGALACTIC MEDIUM

The local heating rate can be estimated in at least two ways, both of which give similar results. First, since the high-energy gamma rays deposit their energy locally, we can identify the local heating rate with the high-energy gamma-ray luminosity density of TeV blazars. In Paper I we showed that this is roughly 2.1 × 10⁻³ that of quasars, and thus approximately ∼(0.5–1.4) × 10³⁸ erg Mpc⁻³ s⁻¹ = (3.4–9.7) × 10⁻⁸ eV cm⁻³ Gyr⁻¹, depending on the minimum luminosity at which the heating mechanism operates.

Alternatively, given a sufficiently complete sample, we can estimate the local heating rate using that implied by the fluxes of the observed TeV blazars. In the present epoch, f(F_E, E, z) ≃ 1 at the relevant E, and thus the heating rate takes the particularly simple form

$$\begin{equation} \left.\dot{Q}\right|_{z=0}\simeq \sum _{\rm AGN} \int dE \theta \frac{F_{E,i}}{D_{\rm pp}}\,. \end{equation} \tag{ 11 }$$

To evaluate this sum, we have collated the 46 extragalactic VHEGR sources currently known.¹² Of these 46, only 28 have published measurements of their VHEGR flux from a combination of VERITAS, HESS, and MAGIC observations. For these 28 sources, we have extracted the parameterized spectra assuming the form

$$\begin{equation} \frac{dN}{dE} = f_0 \left(\frac{E}{E_0}\right)^{-\alpha }, \end{equation} \tag{ 12 }$$

where f₀ is the normalization in units of cm⁻² s⁻¹ TeV⁻¹. The gamma-ray energy flux is trivially related to dN/dE by F_E = EdN/dE∝E^{1 − α}, from which we obtain a VHEGR flux,

$$\begin{equation} F = E_0 f_0 \int _{100\,{\rm G}{\rm eV}}^{10\,{\rm T}{\rm eV}} dE\,\left(\frac{E}{E_0}\right)^{1-\alpha }, \end{equation} \tag{ 13 }$$

and for sources with a measured redshift a corresponding isotropic-equivalent luminosity, L = 4πD²_LF, where D_L is the luminosity distance. The resulting f₀, E₀, α, F, and L are collected in Table 1. In addition, we list the redshift, inferred distance, and absorption-corrected intrinsic spectral index at E₀, obtained via

$$\begin{equation} \hat{\alpha } = {-}\!\left.\frac{d\ln E^{-\alpha } e^{\tau \left(E,z\right)}}{d\ln E}\right|_{E_0} \simeq \alpha - \tau \left(E_0,z\right), \end{equation} \tag{ 14 }$$

where τ(E, z) is the optical depth accrued by a VHEGR emitted at redshift z and with observed energy E.¹³ For high-redshift sources $\hat{\alpha }$ can be less than 2, implying that an intrinsic spectral upper cutoff must exist. Here we conservatively take this to be at E ≃ 10 TeV, which is well justified given the distances to the two sources that dominate the observed local TeV flux, Mkn 421 and 1ES 1959+650, though as we shall see below, the heating rate is relatively insensitive to the particular values of the lower and upper spectral cutoffs.

Table 1. List of TeV Sources with Measured Spectral Properties in Decreasing 100 GeV–10 TeV Flux Order

Name	z	DCa	f0b	E0c	αd	Fe	log10Lf	$\hat{\alpha }$g	$\dot{q}$h	Classi	Reference
Mkn 421	0.030	129	68	1	3.32	1.7 × 103	45.6	3.15	44	H	Chandra et al. (2010)
1ES 1959+650	0.047	201	78	1	3.18	1.6 × 103	45.9	2.90	47	H	Aharonian et al. (2003)
1ES 2344+514	0.044	190	120	0.5	2.95	2.3 × 102	45.0	2.82	9.2	H	Albert et al. (2007c)
Mkn 501	0.034	150	8.7	1	2.58	85	44.4	2.39	6.0	H	Huang et al. (2009)
3C 279	0.536	2000	520	0.2	4.11	68	46.9	2.53	1.0	Q	MAGIC Collaboration et al. (2008)
PKS 2155-304	0.116	490	1.81	1	3.53	64	45.4	2.75	1.4	H	HESS Collaboration et al. (2010)
PG 1553+113	>0.09	>380	46.8	0.3	4.46	41	>44.9	<4.29	<3.88	H	Aharonian et al. (2008b)
W Comae	0.102	430	20	0.4	3.68	31	44.9	3.41	0.6	I	Acciari et al. (2009c)
3C 66A	0.444	1700	40	0.3	4.1	28	46.3	2.43	0.4	I	Acciari et al. (2009b)
1ES 1011+496	0.212	870	200	0.2	4	26	45.5	3.66	0.4	H	Albert et al. (2007b)
1ES 1218+304	0.182	750	11.5	0.5	3.07	24	45.4	2.37	0.8	H	Acciari et al. (2010c)
Mkn 180	0.045	190	45	0.3	3.25	20	44.0	3.17	0.6	H	Albert et al. (2006)
1H 1426+428	0.129	540	2	1	2.6	20	45.0	1.71	1.4	H	Aharonian et al. (2002)
RGB J0710+591	0.125	520	1.36	1	2.69	15	44.8	1.83	0.9	H	Acciari et al. (2010b)
1ES 0806+524	0.138	580	6.8	0.4	3.6	10	44.7	3.21	0.2	H	Acciari et al. (2009a)
RGB J0152+017	0.080	340	0.57	1	2.95	8.5	44.1	2.45	0.3	H	Aharonian et al. (2008a)
1ES 1101-232	0.186	770	0.56	1	2.94	8.2	44.9	1.50	0.3	H	Aharonian et al. (2007c)
1ES 0347-121	0.185	770	0.45	1	3.1	8.2	44.9	1.67	0.3	H	Aharonian et al. (2007a)
IC 310	0.019	83	1.1	1	2.0	8.1	42.8	1.90	0.1	H	Aleksić et al. (2010)
PKS 2005-489	0.071	300	0.1	1	4.0	8.0	44.0	3.56	0.1	H	Aharonian et al. (2005)
MAGIC J0223+430	–	–	17.4	0.3	3.1	7.6	–	<3.1	0.2	R	Aliu et al. (2009)
1ES 0229+200	0.140	590	0.7	1	2.5	6.4	44.5	1.51	0.5	H	Aharonian et al. (2007b)
PKS 1424+240	<0.66	<2400	51	0.2	3.8	6.3	<46.1	>1.42	0.1	I	Acciari et al. (2010a)
M87	0.0044	19	0.74	1	2.31	5.9	41.4	2.29	0.6	R	Acciari et al. (2008)
BL Lacertae	0.069	290	0.3	1	3.09	5.4	43.8	2.67	0.2	L	Albert et al. (2007a)
H2356-309	0.165	690	0.3	1	3.09	5.4	44.6	1.86	0.2	H	Aharonian et al. (2006b)
PKS 0548-322	0.069	290	0.3	1	2.86	4.0	43.7	2.44	0.2	H	Aharonian et al. (2010)
Centaurus A	0.0028	12	0.245	1	2.73	2.8	40.7	2.72	0.2	R	Raue et al. (2010)

Notes. ^aComoving distance in units of Mpc. ^bNormalization of the observed photon spectrum that we assume to be of the form dN/dE = f₀(E/E₀)^−α, in units of 10⁻¹² cm⁻²s⁻¹TeV⁻¹. ^cEnergy at which we normalize the spectrum in units of TeV. ^dObserved spectral index at E₀. ^eIntegrated energy flux between 100 GeV and 10 TeV in units of 10⁻¹² erg cm⁻²s⁻¹. ^fInferred isotropic integrated luminosity between 100 GeV and 10 TeV in units of erg s⁻¹. ^gInferred intrinsic spectral index at E₀. ^hLocal plasma instability heating rate in units of 10⁻¹⁰eV cm⁻³ Gyr⁻¹. ⁱH, I, L, Q, and R correspond to high-energy, intermediate-energy, low-energy-peaked BL Lac objects, flat-spectrum radio quasars, and radio galaxies of Faranoff–Riley Type I (FR I), respectively.

Download table as:  ASCIITypeset image

Using the observed VHEGR source in Table 1, the heating rate associated with the 23 IBL/HBL blazars (labeled I and H in the table; see below) is

$$\begin{equation} \dot{Q}_{\rm obs} \simeq \sum _{\rm obs AGN} \frac{E_0^2 f_0}{D_{\rm pp}(E_0,0)} \int _{100\,{\rm G}{\rm eV}}^{10\,{\rm T}{\rm eV}} \frac{dE}{E_0} \left(\frac{E}{E_0}\right)^{2-\alpha }\,. \end{equation} \tag{ 15 }$$

Note that for α ≃ 3, characteristic of the two dominant TeV sources, this depends logarithmically on the lower and upper spectral cutoffs. The resulting $\dot{Q}_{\rm obs}$ is shown in Figure 1 as a function of the upper spectral cutoff, E_max, for a handful of different lower cutoffs, ranging from 50 GeV to 200 GeV. For our fiducial values, $\dot{Q}_{\rm obs}=1.2\times 10^{-8}\,{\rm eV}\,{\rm c}{\rm m}^{-3}\,{\rm G}{\rm yr}^{-1}$, with this only weakly depending on the various cutoffs, changing at most by a factor of two over the reasonable ranges. Hence, for the remainder of this work, we choose a lower spectral cutoff of 100 GeV.

Zoom In Zoom Out Reset image size

To complete this estimate, we now correct for the selection effects of VHEGR observations. To correct for the pointed nature of VHEGR observations, we rely on the all-sky GeV gamma-ray observations from the Fermi satellite (Abdo et al. 2010b) of TeV blazars. Those belong to two subclasses of blazars, namely, high-energy-peaked BL Lac (HBL) objects and the somewhat less efficient accelerators, intermediate-energy-peaked BL Lac (IBL) objects, which are, in some cases, also able to reach energies beyond 100 GeV. Outside of the Galactic plane, Fermi observes 118 high-synchrotron-peaked (HSP) blazars and a total of 46 intermediate-synchrotron-peaked (ISP) blazars.¹⁴ Roughly half of the latter are likely to emit VHEGRs as indicated by their flat spectral index between 0.1 and 100 GeV, Γ ≲ 2 (see the spectral index distribution of Figure 14 in Abdo et al. 2009a). Of these potential 141 TeV blazars, only 22 have been coincidentally identified as TeV blazars (out of a total of 28 coincident TeV sources), while there are a total of 33 known TeV blazars (29 HBL, 4 IBL). If these 141 sources are all TeV emitters but have not been detected due to incomplete sky coverage of current TeV instruments, then the selection factor is η_sel = 141/33 = 4.3. In addition, the duty cycle of coincident GeV and TeV emission is η_duty = 33/22 = 1.5. Here we assume that the luminosity distribution of observed VHEGR sources reflects the true distribution after correcting for the effects of flux limitations in the observations. That this may be done is justified empirically in Section 5.1.2 of Paper I and demonstrated explicitly in Figure 5 of Paper I. Thus, we may use constant correction factors (independent of luminosity) to estimate the true distribution. Finally, by excluding the Galactic plane for galactic latitudes b < 10°, this is an underestimate by roughly η_sky = 1.17. Taken together, the true heating rate in our "standard" model is then

$$\begin{eqnarray} \left.\dot{Q}\right|_{z=0} &\simeq& 7\times 10^{-8} \left(\frac{\eta _{\rm sys}}{0.8}\right)\nonumber\\ &&\times \left(\frac{\eta _{\rm sel}}{4.3}\right) \left(\frac{\eta _{\rm duty}}{1.5}\right) \left(\frac{\eta _{\rm sky}}{1.2}\right) \,{\rm eV}\,{\rm c}{\rm m}^{-3}\,{\rm G}{\rm yr}^{-1},\quad \end{eqnarray} \tag{ 16 }$$

where η_sys is a remaining coefficient of order unity correcting for systematic uncertainties. These include corrections to the spectral model we currently use and the corrections to the completeness of the Fermi sample of HBLs and IBLs in accounting for the complete population of VHEGR sources. The fact that there are already four radio galaxies of Faranoff–Riley Type I, two flat-spectrum radio quasars, and two starburst galaxies emitting VHEGRs implies that we are probably too conservative in accounting for all the VHEGR sources. To explore the effects of the variation in the heating rate, we also adopt an "optimistic" heating model that is normalized to fit the observed IGM inverted temperature–density relation (Viel et al. 2009). This "optimistic" model has η_sys = 1.6, yielding a heating rate of $\dot{Q}|_{z=0} = 1.4\times 10^{-7}\,{\rm eV}\,{\rm c}{\rm m}^{-3}\,{\rm G}{\rm yr}^{-1}$.

In estimating $\dot{Q}$ from Equation (11), we have made a number of implicit assumptions. First, we have assumed that the observed distribution of TeV blazars locally is representative of the average distribution at any given point in the universe. Since the blazars that have been observed locally go out to z ≈ 0.5, we believe that this is a relatively safe assumption. Second, we have assumed that the current sample of blazars is sufficiently flux-complete to dominate the heating rate at Earth. That is, we have assumed that the effective flux limit of the current generation of Cerenkov telescopes is low enough that it captures the bulk of the sources responsible for the local heating of the IGM. Whether this is the case is less clear and will be the subject of Section 3.2 and the Appendix in some detail. Unfortunately, since the constant of proportionality relating the local TeV blazar and quasar luminosity densities is obtained using the observed set of blazars, our two estimates are strongly correlated.

Nevertheless, this provides a convenient lower limit on $\dot{Q}$, and we use this heating rate, which we denote the "standard" heating rate, as the mean heating rate of any fluid element in the present-day universe. Over a Hubble time, the total heat deposited per unit volume is then roughly

$$\begin{equation} \left. u_{\rm tot}\right|_{z=0} \simeq \frac{\left.\dot{Q}\right|_{z=0}}{H_0} \approx 10^{-6}\,{\rm eV}\,{\rm c}{\rm m}^{-3}, \end{equation} \tag{ 17 }$$

which is sufficient to raise the temperature of 1 + δ = 0.1 regions to approximately 1.7 × 10⁵ K and thus dramatically alter the thermal history of voids.

Extrapolating the above estimate to z > 0 requires understanding how the TeV blazar properties and population have evolved, as well as the evolution in f(F_E, E, z). Generally, the average heating rate is given by

$$\begin{equation} \dot{Q}= \int dV\,d\log _{10}L\,d\alpha \,d\Omega \, \tilde{\phi }_B(z;L,\alpha ,\Omega) \frac{\Omega }{2\pi } \dot{q}, \end{equation} \tag{ 18 }$$

where $\tilde{\phi }_B(z;L,\alpha ,\Omega ,z)$ is the physical number density of blazars at a given redshift per unit logarithmic isotropic-equivalent luminosity, spectral index, and blazar jet opening angle (we have assumed that all jets are symmetric). We make the simplifying assumption that $\tilde{\phi }_B$ is separable into components describing the evolving luminosity density distribution, $\tilde{\phi }_B(z,L)$, and a static, unit-normalized spectral distribution, φ_B(α, Ω), where $\tilde{\phi }_B(z;L,\alpha ,\Omega)=\tilde{\phi }_B(z,L)\varphi _B(\alpha ,\Omega)$. This produces

$$\begin{eqnarray} \dot{Q}&=& \int d\log _{10} L\, L \tilde{\phi }_B(z,L)\nonumber\\ &&\times \int D^2 dD\,d\alpha \,d\Omega \, \int dE\, \frac{\theta }{D_{\rm pp}} f \frac{L_E e^{-D/D_{\rm pp}}}{L D^2}, \end{eqnarray} \tag{ 19 }$$

where the bulk of the redshift and luminosity dependence is now contained in the outermost integral over log₁₀L, and we have used τ_E ≃ D/D_pp (since the heating is dominated by nearby objects, for simplicity we do not distinguish between different cosmological distance definitions, though see the Appendix for a more careful treatment). If there were no redshift or flux dependence in the remaining terms, it would be possible to simply normalize the heating rate by $\dot{Q}|_{z=0}$ and the estimated TeV blazar luminosity density, which we define to be

$$\begin{equation} \tilde{\Lambda }_B(z)\equiv \int _{\log _{10}L_{\rm min}}^\infty d\log _{10}L\, L \tilde{\phi }_B(z,L), \end{equation} \tag{ 20 }$$

where L_min ≃ 3 × 10⁴² erg s⁻¹ is chosen so that the plasma instabilities operate efficiently at all redshifts of interest. However, this is not entirely the case since f(F_E, E, z) retains some dependence on z, and thus we set

$$\begin{equation} \dot{Q}= \frac{\tilde{\Lambda }_B(z)}{\tilde{\Lambda }_B(0)} \left.\dot{Q}\right|_{z=0} \dot{\mathcal {Q}}_{\rm corr}(z), \end{equation} \tag{ 21 }$$

where $\dot{\mathcal {Q}}_{\rm corr}$ provides a correction due to the changing strength of the pair beam dissipation mechanism in comparison to inverse Compton cooling.

The VHEGR flux at Earth is dominated by a handful of very bright sources for which the absorption-corrected α ≃ 3 at 1 TeV and isotropic-equivalent luminosity is 5 × 10⁴⁵ erg s⁻¹. This simplifies the estimation of $\dot{\mathcal {Q}}_{\rm corr}$ significantly, giving

$$\begin{eqnarray} \dot{\mathcal {Q}}_{\rm corr}(z) &\simeq& \int\!\! \int dE dD \frac{\theta E^{-2}}{D_{\rm pp}} e^{-D/D_{\rm pp}} \nonumber\\ &&\times f\left(\frac{L_E e^{-D/D_{\rm pp}}}{4\pi D^2}, E, z \right)\bigg / \int dE \theta E^{-2},\quad \end{eqnarray} \tag{ 22 }$$

in which the form and magnitude of L_E are fixed. The resulting correction factors are shown in Figure 2 for isotropic-equivalent luminosities above 100 GeV ranging from 10⁴⁴ erg s⁻¹ to 10⁴⁶ erg s⁻¹. Generally, $\dot{\mathcal {Q}}_{\rm corr}$ makes a small (≲ 15%) correction following the peak of the quasar luminosity function (z = 2), though it can grow to as much as 80% by z = 6 for dim objects. However, since we will be primarily interested in bright blazars at z ≲ 4, $\dot{\mathcal {Q}}_{\rm corr}$ modifies the heating rate at relevant redshifts by ≲ 50% (see also Section 3.3). Thus, given the comparatively larger uncertainties associated with the TeV blazar luminosity function, plasma heating mechanism, and additional potential VHEGR sources, we neglect $\dot{\mathcal {Q}}_{\rm corr}$ in what follows.

Zoom In Zoom Out Reset image size

Therefore, obtaining the heating rate as a function of z is reduced to estimating $\tilde{\Lambda }_B(z)/\tilde{\Lambda }_B(0)$, the normalized evolution of the blazar luminosity density. However, no VHEGR emitting blazars are known beyond z ≃ 0.7, presumably due to the pair-creation absorption associated with propagation through the EBL. As a consequence, nothing is known about $\tilde{\Lambda }_B(z)$ at high redshifts directly. Nevertheless, in Paper I we showed that a close relationship between $\tilde{\phi }_B$ and the quasar luminosity function of Hopkins et al. (2007), $\tilde{\phi }_Q$, exists at low z:

$$\begin{equation} \tilde{\phi }_B(0.1,L) \simeq 3.8\times 10^{-3}\,\tilde{\phi }_Q(0.1,1.8L)\,. \end{equation} \tag{ 23 }$$

The insert in Figure 3 shows a comparison between the observed TeV blazar luminosity function, calculated from the distribution of the sources in Table 1, and applying an empirically determined flux limit of 4.19 × 10⁻¹² erg cm⁻² s⁻¹, with ϕ_B(0.1, L) (for a full description of how the luminosity function was formed and how it compares to ϕ_Q, see Section 5 of Paper I). Furthermore, we showed that once the ICCs are suppressed by the plasma beam instabilities (or some analogous mechanism), extending this relationship to high z was in excellent agreement with the best current constraints on the high-z TeV blazar population: the Fermi TeV blazar statistics and the Fermi measurement of the EGRB (between 100 MeV and 100 GeV). Thus, we estimate Λ_B(z) ≃ 2.1 × 10⁻³Λ_Q(z), shown in Figure 3, where Λ_Q(z) is the luminosity density of quasars.

Zoom In Zoom Out Reset image size

Figure 3 shows that employing the quasar luminosity density to estimate the heating rates has profound consequences. In particular, the inferred luminosity density of TeV blazars (solid line) rises rapidly with increasing z, with the comoving density increasing by roughly a factor of ∼10 by z = 2 (see also Figure 8 of Hopkins et al. 2007). In physical units this corresponds to a increase by a factor of nearly 300. Thus, we expect an increase in the local heating rate by a similar factor over the currently estimated value near z = 2.

We begin by estimating the instantaneous comoving mean separation of visible blazars above various luminosity thresholds,

$$\begin{equation} D_B(L_m,z) = (1+z)\left[\int _{\log _{10}L_m}^{\log _{10}L_M} \tilde{\phi }_B(L,z) d\!\log _{10}L\right]^{-1/3}, \end{equation} \tag{ 24 }$$

where the precise value of L_M is unimportant as long as it is above the peak of the luminosity function; here we choose L_M = 2 × 10⁴⁶ erg s⁻¹, consistent with the theoretically expected upper limit of TeV blazar luminosities. Since we have determined $\tilde{\phi }_B$ from the observed blazar population and are making use of the isotropic-equivalent luminosities, this is independent of the blazar jet opening angle; smaller opening angles will result in a correspondingly larger number of objects such that the total number seen is unchanged, and therefore D_B is fixed. A comparison between D_B and the mean free path of VHEGRs is shown in Figure 4. Because D_pp varies dramatically from 100 GeV to 10 TeV, for this purpose, we define a spectrally averaged mean free path, $\bar{D}_{\rm pp}$, determined implicitly by

$$\begin{equation} \int _{100\,{\rm G}{\rm eV}}^{10\,{\rm T}{\rm eV}} dE E^{1-\alpha } e^{-\bar{D}_{\rm pp}(z)/D_{\rm pp}(E,z)} = e^{-1} \int _{100\,{\rm G}{\rm eV}}^{10\,{\rm T}{\rm eV}} dE E^{1-\alpha }, \end{equation} \tag{ 25 }$$

in which we choose α = 3 based on the local TeV blazar sample. This is roughly the e-folding distance of the entire spectral band, shown in Figure 4 by the gray solid line, and in practice is quite close to D_pp(225 GeV, z). If the intrinsic TeV blazar spectra peak above 100 GeV, our estimate of the typical VHEGR mean free path could be significantly too large; thus, we also compare D_B to D_pp at 1 TeV (the gray dashed line in Figure 4), providing an extreme lower bound on the spectrally averaged mean free path in practice. The rapid increase of the density of EBL photons with z, associated with the larger star formation rate in the recent past, results in a substantially reduced D_pp by z = 1. Prior to z = 1, the physical number density of EBL photons remains nearly constant, and thus D_pp∝(1 + z) in comoving units.

Zoom In Zoom Out Reset image size

In the present epoch, $\bar{D}_{\rm pp}$ is quite large, and thus despite their sparsity, the mean separation of even bright blazars (10⁴⁶ erg s⁻¹) is less than $\bar{D}_{\rm pp}$. This remains true for objects with luminosities ⩽10⁴⁵ erg s⁻¹, which include the local luminosity-weighted median TeV blazar luminosity, L_0.5(z), defined such that

$$\begin{equation} \int _{\log _{10}L_{0.5}(z)}^{\log _{10}L_M} L \tilde{\phi }_B(L,z) d\!\log _{10}L= 0.5\tilde{\Lambda }_B(z). \end{equation} \tag{ 26 }$$

Hence, locally, we expect blazar heating to be quite uniform.

At high redshift matters change somewhat. Until z = 1, $\bar{D}_{\rm pp}(z)$ and D_pp(1 TeV, z) both decrease more rapidly than D_B. For z > 1, in comoving units $\bar{D}_{\rm pp}(z)$ and D_pp(1 TeV, z) increase slowly, though at a marginally larger rate than D_B. As a result, near z ∼ 1 the mean separation of TeV blazars is largest in comparison to the VHEGR mean free path. Thus, the mean separation between objects more luminous than 10⁴⁶ erg s⁻¹ is larger than $\bar{D}_{\rm pp}$ near z ∼ 1. Nonetheless, when the lower luminosity limit is dropped to 10⁴⁵ erg s⁻¹ or below, generally $D_B(z)<\bar{D}_{\rm pp}(z)$ for all z < 6. Similarly, the mean separation between blazars more luminous than L_0.5(z) is also smaller than our estimate of $\bar{D}_{\rm pp}$ for the relevant redshift range. Hence, we may expect that nearly all patches of the universe will be illuminated by at least one luminous TeV blazar following z ∼ 6.

Closely related to the mean separation is the number of blazars within the $\bar{\tau }\simeq D/\bar{D}_{\rm pp}=1$ surface (defined explicitly in the Appendix) above some luminosity limit. Figure 5 shows these for objects with luminosities above 10⁴⁵ erg s⁻¹ (blue short-dash) and L_0.5 (red dot). While these may be roughly inferred from the associated mean separations in Figure 4, they differ slightly from $(4\pi /3)(\bar{D}_{\rm pp}/D_B)^3$ due to the rapidly evolving $\bar{D}_{\rm pp}$ and blazar population combined with the finite look-back time in the integral. At z = 0 both are well above unity, with the luminosity-weighted median value substantially exceeding our estimate of roughly 170 visible extragalactic TeV blazars.¹⁵ However, since generally $\bar{D}_{\rm pp}(z)\simeq 4.4D_{\rm pp}(1\,{\rm T}{\rm eV},z)$, not all of the blazars that contribute to the local heating are expected to appear as strong TeV sources.

Zoom In Zoom Out Reset image size

The low-z behavior of the number of blazars is dictated by the rapid evolution in D_pp, associated with the recent rapid variation in the star formation rate (and thus the number density of EBL photons). Prior to z = 1, the EBL density, and hence D_pp, is roughly constant in physical units, and the evolution number of visible blazars becomes indicative of the underlying evolution of blazar population. At all z < 6 these estimates of the number of blazars responsible for the bulk of the heating, ${\mathcal N}_B$, exceed unity, implying only a small fractional spatial variation in the blazar heating rate. However, while the luminosity-weighted median estimate of ${\mathcal N}_B$ does give some idea of which population is responsible for most of the TeV blazar luminosity density, neither encapsulates the population responsible for the majority of the local heating. Thus, if the local heating were dominated by a handful of luminous sources, it is possible for the stochasticity to be much larger.

Therefore, we show a third estimate of ${\mathcal N}_B$, corresponding to the number of sources with $\dot{q}$ values larger than the heating-rate-weighted median value (see the Appendix for a precise definition). That is, at a given redshift, the ${\mathcal N}_B$ sources with the highest $\dot{q}$ values generate half of the total heating rate. Thus, if a single source were to dominate the local heating rate, ${\mathcal N}_B$ would be considerably less than unity, indicating a large degree of variability. However, below z ∼ 3.5 this is not the case; the fractional heating-rate-defined ${\mathcal N}_B$ is considerably larger than unity. At high redshifts the heating rate becomes increasingly dominated by more distant, luminous, and rarer objects, departing from the previous estimate of ${\mathcal N}_B$.

Nevertheless, while for z ≳ 3.5 few sources contribute nearly half of the heating rate, the total heating rate must be comparable to the total VHEGR luminosity density of the blazars. Therefore, the total number of contributing objects must be similar to the number of visible blazars above the luminosity-weighted median L (i.e., the red dotted line in Figure 5). As a consequence, the number of relevant sources must rapidly rise with decreasing fractional heating rate. This is explicitly indicated by the dashed line in Figure 5, which shows the number of sources responsible for 75% of the local heating rate.

In any case, it is clear that for the redshifts at which blazar heating is likely to be important, z ≲ 4, the heating rate will be relatively uniform. At z ∼ 4 it may experience order 50% fluctuations, and by z ∼ 6 it will exhibit order unity deviation. However, we note that all of our estimates of the heating inhomogeneity are predicated on the assumption that the $\tilde{\phi }_B$ accurately represents the distribution of blazars and has the associated considerable uncertainties. For a more thorough discussion of different estimates of the number of contributing blazars and a discussion of which blazars dominate the heating rate at a given redshift, see the Appendix.

In addition to the statistics of blazars, the homogeneity also depends on the intrinsic variability of TeV blazars. During the time that the TeV blazars have been observed, their VHEGR emission has remained remarkably stable. However, this provides only a weak lower limit (∼4 yr) on the variability timescale in these objects (Dermer et al. 2011). In addition to depending on the properties of the source itself, the variability of TeV blazars depends on the primary emission mechanism responsible for the VHEGR component of TeV blazars. There are two classes of inverse Compton models to explain the VHEGR emission.

• 1.  
  Synchrotron self-Compton (SSC) model. In this model, the synchrotron radiation field is Compton upscattered to TeV energies by a relativistic electron population (Jones et al. 1974; Ghisellini & Maraschi 1989). Recent work has led to the conclusion that a simple homogeneous, one-zone, SSC model cannot explain the spectral energy distribution (SED) of the majority of blazars (see Figure 36 of Abdo et al. 2010a). However, models with multiple SSC components are consistent with the blazar SEDs. Typically these invoke a steady, primary SSC component that peaks at the IR/optical (S) and γ-ray band (IC) and a second, more energetic, and usually more variable component, which peaks in the UV or X-ray band (S) and at GeV/TeV energies (IC). The variability of this energetic component increases the probability that a given patch of the IGM will see a TeV blazar during its history, resulting in larger homogeneity in the blazar heating.
• 2.  
  External radiation Compton (ERC) scenario. This model proposes that the relativistic jet electrons Compton scatter an external radiation field (Sikora et al. 1994; Dermer & Schlickeiser 2002) from the accretion disk or the dusty torus surrounding it. In the first case, the disk generates UV seed photons, which are then reflected toward the jet by the broad-line region within a typical distance from the accretion disk of the order of 1 pc. In the second case, the dusty torus could provide IR seed photons that are emitted at larger distances from the jet. In any of these ERC models, the VHEGR emission is expected to be very steady, implying long-lived TeV blazars, resulting in more patchy blazar heating.

We leave detailed studies of the inhomogeneity in the blazar heating rate resulting from intrinsic inhomogeneity and variability in the spatial distributions of the blazars themselves for future work.

Starburst galaxies are characterized by the presence of rapid star formation and hence are sites of numerous supernovae. Supernovae are known to be efficient particle accelerators and are presumed to be the primary source of the galactic cosmic rays. As cosmic rays can produce VHEGRs, starburst galaxies are potentially bright VHEGR sources.

While many starburst galaxies contain active AGNs that may also emit VHEGRs, unambiguous VHEGR emission from the starburst itself has been seen in at least two cases, M82 and NGC 253, with E > 100 GeV luminosities of 5 × 10³⁹ erg s⁻¹ and 2 × 10³⁹ erg s⁻¹, respectively (VERITAS Collaboration et al. 2009; Acero et al. 2009). The luminosities of these two sources are dwarfed by the local TeV blazars, and thus the VHEGR flux from M82 and NGC 253 does not contribute to the heating of the local IGM. In addition, their VHEGR luminosity produces a pair beam that is too dilute for plasma processes such as the "oblique" instability to beat cooling via inverse Compton. Hence, they do not contribute to the local VHEGR flux.

However, very high star-forming systems such as the ultraluminous infrared galaxies (ULIRGs), where L_IR > 10¹² L_☉ and L_IR is the bolometric infrared luminosity (5 μm < λ < 1000 μm), may evade this constraint. If the VHEGR emission from starburst galaxies is due to cosmic rays accelerated by supernovae, the VHEGR luminosity above 100 GeV, L, is then proportional to the star formation rate. For starbursts, this is linearly related to the continuum infrared luminosity (Kennicutt 1998). Thus, normalizing by M82 and NGC 253, we have

$$\begin{equation} L\simeq 6\times 10^{40} \left(\frac{L_{\rm IR}}{10^{12}\,L_\odot }\right)\,{\rm erg}\,{\rm s}^{-1}\,. \end{equation} \tag{ 27 }$$

While this relationship is extremely uncertain (the normalization varies by a factor of two between M82 and NGC 253), it suggests that ULIRGs, or perhaps hyperluminous infrared galaxies (L_IR > 10¹³ L_☉), may be sufficiently bright to contribute to the heating of the IGM.

While the fluxes from the brightest ULIRGs remain much smaller than those associated with the typical TeV blazars, there are many more starburst galaxies than AGNs. At all redshifts ULIRGs constitute the high-luminosity tail of the star-forming galaxy luminosity functions (Le Floc'h et al. 2005; Caputi et al. 2007; Magnelli et al. 2009; Goto et al. 2010). In the present epoch, the density of ULIRGs is roughly 4 × 10⁻⁷ Mpc⁻³, and thus the corresponding VHEGR luminosity density of these objects, ∼2 × 10³⁴ erg s⁻¹ Mpc⁻³ (see Table 8 of Caputi et al. 2007), is negligible in comparison to the TeV blazars, roughly 5 × 10³⁷ erg s⁻¹ Mpc⁻³.

However, due to the steep decline of the luminosity function in the ULIRG range, small changes in the location of the break luminosity result in large changes in the comoving luminosity density of ULIRGs. As a consequence, the comoving luminosity density of ULIRGs grows much faster than that of quasars (and thus presumably the TeV blazars). Nevertheless, even at z = 2, roughly the redshift of peak star formation, the comoving number density of ULIRGs remains below ∼2 × 10⁻⁴ Mpc⁻³ (see Table 8 of Caputi et al. 2007), corresponding to a comoving VHEGR luminosity density of ≲ 10³⁷ erg s⁻¹ Mpc⁻³, more than two orders of magnitude smaller than the contemporaneous TeV blazar population. For this reason, we neglect the starburst contribution to heating the IGM here, though they may represent a secondary source class.

Stellar-mass objects, such as magnetars and X-ray binaries, generally have difficulty reaching the flux limits required for plasma cooling to dominate the pair beam evolution. Even at Eddington-limited VHEGR luminosities, reaching isotropic-equivalent luminosities of 10⁴² erg s⁻¹ requires beaming factors of roughly 10⁴, corresponding to jet opening angles of roughly 2°. In practice, the formation of radio jets is believed to be associated with substantially sub-Eddington accretion flows in these objects, thus exacerbating the beaming requirement. More importantly, these sources may suffer compactness problems: the VHEGR emission regions in stellar-mass jets must necessarily be very far from the central object to avoid in situ pair production. Finally, were X-ray binaries and magnetars generally strong, persistent VHEGR emitters with sufficiently large fluxes, we would expect that many would have been already detected.

Gamma-ray bursts (GRBs) are natural candidates due to their large luminosities and strong inferred beaming, and we consider them as an example of the general class of energetic transient sources. Unfortunately, little is known about the VHEGR emission of GRBs. Currently there is a single report of a TeV signal associated with a GRB (GRB 970417A, Atkins et al. 2000, 2003), though due to the large distances at which they can be observed and comparative rarity, this is not unexpected. However, such high-energy emission is possible in principle, presumably due to inverse Comptonization of the prompt emission and/or X-ray afterglow (see Section VIII of Piran 2004, and references therein). Fermi observations of GRBs have shown that for many events the Band spectrum can be extended to ∼100 GeV (Abdo et al. 2009b, 2009c, 2009d; Ackermann et al. 2010), though in at least one case a spectral break below 10 GeV has been observed (GRB 090926A, Ackermann et al. 2011). Thus, it remains unclear if in practice the high-energy emission is attenuated within the emission region. Moreover, since GRBs are inherently short-lived events, the luminosity limits described in Paper I, which require the VHEGR emitting phase to last for a plasma cooling timescale (roughly 10²–10³ yr), are not directly applicable. A similar analysis, obtained by limiting the beam growth time to the GRB duration, gives VHEGR isotropic-equivalent energies of 10⁵⁴ erg. This is comparable to the total prompt and afterglow emission for only the brightest bursts, composing roughly 5% of GRBs observed by Swift (Gehrels et al. 2009). Nevertheless, even assuming that all GRBs produce the requisite high-energy emission, at an optimistic present local rate of roughly 0.5 Gpc⁻³ yr⁻¹, produces a comoving luminosity density of ≲ 10³⁶ erg s⁻¹ Mpc⁻³, roughly three orders of magnitude less than that due to TeV blazars at z = 1.

Our first definition is the simplest one can imagine; choose an intrinsic luminosity range and use the local density to determine the number within a volume defined by the spectrally average mean free path. That is,

$$\begin{equation} {\mathcal N}_{B,I}(z_o;L_m) = \frac{4\pi }{3} \bar{D}_{\rm pp}^3(z_o) \tilde{\Phi }_B(z_o;L_m,L_M)\,. \end{equation} \tag{ A9 }$$

This is an approximation of the number of sources with intrinsic luminosities in the specified range within an approximation of a single optical depth. This is shown for L_m = 10⁴⁵ erg s⁻¹ by the cyan dash-dotted line in Figure 11. Typically, where the blazar population is rapidly evolving it tends to be a poor estimate of the number of sources, overestimating this number by nearly an order of magnitude at z ≳ 3.

As long as $\tilde{\phi }_B(z,L)$ evolves slowly and $\bar{D}_{\rm pp}\ll c/H_0$, it is unnecessary to perform the relevant redshift integral to get the volume element. However, this is not always the case, and thus a more accurate estimate is obtained by integrating the blazar number density within the volume specified by unit optical depth. That is,

$$\begin{equation} {\mathcal N}_{B,II}(z_o;L_m) = \int _{z_o}^{z_1} 4\pi D_A^2(z;z_o) \frac{dD_P}{dz} \tilde{\Phi }_B(z;L_m,L_M) \,dz, \end{equation} \tag{ A10 }$$

where z₁ is defined implicitly by $\bar{\tau }(z_1;z_o)=1$.

For this definition, we may choose L_m in a variety of ways. In Figure 11 we show two in particular: that from setting L_m = 10⁴⁵ erg s⁻¹ (blue short-dashed line), which may be directly compared with the case shown for ${\mathcal N}_{B,I}$, and setting L_m = L_0.5 as defined by Equation (26), i.e., the luminosity-weighted median luminosity, above which sources produce half of the total local luminosity density (red dotted line). For z ≲ 3 this approximation for ${\mathcal N}_B$ gives similar results to ${\mathcal N}_{B,I}$ for a fixed L_m. At high redshifts, where $\tilde{\phi }_B$ is rapidly decreasing, the two diverge substantially. The difference between the ${\mathcal N}_{B,II}(z;L_{0.5})$ and the previous two is more striking and a consequence for the assumed evolving luminosity distribution of blazars, which is clearly important.

While a fixed intrinsic luminosity limit may be useful conceptually, even idealized surveys do not directly probe such a population. Instead, most surveys are flux limited, and thus we also define a flux-limited definition of ${\mathcal N}_B$. In this case we set the minimum luminosity via

$$\begin{equation} L_m(z;z_o,F_m) = 4\pi \left(\frac{1+z}{1+z_o}\right)^{\alpha -2} D_L^2(z;z_o) e^{\bar{\tau }(z;z_o)} F_m , \end{equation} \tag{ A11 }$$

where F_m is a fixed flux limit, from which we obtain

$$\begin{eqnarray} {\mathcal N}_{B,III}(z_o;F_m) &=& \int _{z_o}^{\infty } 4\pi D_A^2(z;z_o) \frac{dD_P}{dz}\nonumber\\ &&\times\, \tilde{\Phi }_B\left[z;L_m(z;z_o,F_m),L_M\right] \,dz.\qquad \end{eqnarray} \tag{ A12 }$$

This corresponds to the number of objects a flux-limited survey (in the 100 GeV–10 TeV band) performed by an observer at redshift z_o would detect. As such, it is the most directly comparable to the number of TeV blazars that have been observed. This is shown for a flux limit comparable to that inferred for the TeV sample in Table 1, 4.19 × 10⁻¹² erg cm⁻² s⁻¹ in Figure 11 by the green long-dash-dotted line. Note that the number of objects found at z = 0 corresponds nicely to the number observed after correcting for incompleteness of present TeV surveys. Of course, this is by construction since $\tilde{\phi }_B$ was obtained from the observed population. Nevertheless, it is striking that many more objects would have been observed above this flux limit during earlier epochs, vastly exceeding any of the preceding approximations of ${\mathcal N}_B$. However, it does not necessarily follow that all of these sources will have contributed substantially to the heating rate.

The most relevant approximation for ${\mathcal N}_B$, and the one we adopt as our primary definition, is set by the heating rates directly. The general idea is to do what we do naturally at Earth: arrange all of the sources visible by an observer at z_o by the local heating rate they induce, from largest to smallest, and count until a fixed fraction of the total heating rate is reached. That is, set ${\mathcal N}_B$ to be the minimum number of sources (on average) required to produce a given fraction of the total heating rate. Given its direct connection with the heating rate, this provides the most natural definition of the number of sources "responsible for the bulk of the heating." To do this, however, we first must explicitly define the heating rate in terms of the appropriate functions.

Given a fixed spectrum, there is a linear relationship between the local flux and the heating rate, defined by Equation (9):

$$\begin{eqnarray} \dot{q} &= \int _{E_m}^{E_M} \frac{F_E e^{-\tau (E,z;z_o)}}{D_{\rm pp}\left[E(1+z)/(1+z_o),z\right]} dE = \chi (z;z_o) F(E_m,E_M)\nonumber\\ &= \chi (z;z_o) \left(\frac{1+z}{1+z_o}\right)^{2-\alpha } \frac{L}{4\pi D_L^2(z;z_o)}, \end{eqnarray} \tag{ A13 }$$

where

$$\begin{equation} \chi (z;z_o) \equiv\! \int _{E_m}^{E_M}\! \frac{ E^{1-\alpha } e^{-\tau (E,z;z_o)} }{ D_{\rm pp}\left[E(1+z)/(1+z_o),z\right] } dE\! \bigg /\!\! \int _{E_m}^{E_M} E^{1-\alpha } dE \end{equation} \tag{ A14 }$$

is a function of the shape of the spectrum and the redshifts, similar to $\bar{\tau }(z;z_o)$. Thus, a given heating rate defines an intrinsic luminosity limit for a source at a given redshift:

$$\begin{equation} L_m(z;z_o,\dot{q}_m) = 4\pi \left(\frac{1+z}{1+z_o}\right)^{\alpha -2} \frac{D_L^2(z;z_o)}{\chi (z;z_o)} \dot{q}_m \,. \end{equation} \tag{ A15 }$$

With this, we may compute the heating rate as a function of $\dot{q}_m$, i.e., the heating rate associated with sources that produce a local heating larger than some limit:

$$\begin{eqnarray} \dot{Q}(z_o;\dot{q}_m) &=& \int _{z_o}^{\infty } 4\pi D_A^2(z;z_o) \frac{dD_P}{dz}\nonumber\\ &&\times \int _{\log _{10}L_m(z;z_o,\dot{q}_m)}^{\log _{10}L_M}\hspace{-19.91684pt} \dot{q} \,\tilde{\phi }_B(z,L)\,d\!\log _{10}L\,dz\nonumber\\ &=& \int _{z_o}^{\infty } 4\pi D_A^2(z;z_o) \frac{dD_P}{dz} \left(\frac{1+z}{1+z_o}\right)^{2-\alpha }\nonumber\\ && \times \frac{\chi (z;z_o)}{4\pi D_L^2(z;z_o)} \tilde{\Lambda }_B\left[z;L_m(z;z_o,\dot{q}_m),L_M\right] \,dz.\nonumber\\ \end{eqnarray} \tag{ A16 }$$

The heating from all gamma-ray blazars is obtained simply by setting $\dot{q}_m=0$. On the other hand, we may set $\dot{q}_m$ implicitly via

$$\begin{equation} \dot{Q}(z_o;\dot{q}_m) = \mathcal {Q} \, \dot{Q}(z_o;0), \end{equation} \tag{ A17 }$$

where $\mathcal {Q}$ ranges from 0 to 1, yielding a heating rate limit at each redshift that we shall call $\dot{q}_{\mathcal {Q}}$. From this, we may then obtain a number of contributing blazars:

$$\begin{eqnarray} {\mathcal N}_{B,IV}(z_o;\mathcal {Q}) &=& \int _{z_o}^{\infty } 4\pi D_A^2(z;z_o) \frac{dD_P}{dz}\nonumber\\ &&\times \tilde{\Phi }_B\left[z;L_m(z;z_o,\dot{q}_{\mathcal {Q}}),L_M\right] \,dz.\nonumber\\ \end{eqnarray} \tag{ A18 }$$

This is shown for $\mathcal {Q}=0.5$ and $\mathcal {Q}=0.75$ in Figure 11 by the black solid and long-dashed lines, respectively. While both are similar to the other measures of ${\mathcal N}_B$ below z ∼ 1, above this redshift they fall much more rapidly. This is due to the shift of $\tilde{\phi }_B(z,L)$ toward higher luminosities, and thus the luminosity density (and therefore heating rate) becomes dominated by fewer, more luminous sources. Nevertheless, ${\mathcal N}_B$ is a rapidly increasing function of $\mathcal {Q}$, as evidenced by the fact that ${\mathcal N}_{B,IV}$ increases by approximately an order of magnitude when $\mathcal {Q}$ increases from 0.5 to 0.75.

We begin with a toy problem in which we consider a fixed power-law luminosity function in a static, Euclidean universe. Specifically, we choose

$$\begin{equation} \tilde{\phi }(L) = \tilde{\phi }_* \left(\frac{L}{L_*}\right)^\xi \Theta \left(\frac{L_*}{L}\right), \end{equation} \tag{ A19 }$$

where Θ(x) is the Heaviside function, $\tilde{\phi }_*$ is the overall normalization of the luminosity function, L_* is a maximum luminosity (approximating a break), and ξ is an arbitrary constant. We will also assume that D_pp is independent of energy, i.e., there is some characteristic value for which τ = D/D_pp and we may bring it out of the energy integral in the definition of $\dot{q}$. As a consequence, a fixed limit in $\dot{q}_m$ corresponds directly to a fixed flux limit, F_m.

In a Euclidean universe the number of objects takes the particularly simple form

$$\begin{eqnarray} {\mathcal N}_E &=& \int _0^\infty dD 4\pi D^2 e^{-D/D_{\rm pp}} \int _{\log _{10}L_m}^\infty \tilde{\phi }(L) \,d\!\log _{10}L\nonumber\\ &=& \int _0^\infty dD 4\pi D^2 e^{-D/D_{\rm pp}}\, \Theta \left(\frac{L_*}{L_m}\right)\nonumber\\ &&\times \int _{\log _{10}L_m}^{\log _{10}L_*} \tilde{\phi }_* \left(\frac{L}{L_*}\right)^\xi \,d\!\log _{10}L. \end{eqnarray} \tag{ A20 }$$

The flux limit, F_m, gives a luminosity limit of L_m = 4πD²F_m, where we have ignored the optical depth. Since we will be most concerned with how peaked the various integrands are at nearby distances, this is not a significant oversight (including it would serve to make them only more so). The definition of L_m implies a maximum distance as well, with $D_*=\sqrt{L_*/4\pi F_m}$, and therefore L_m/L_* = (D/D_*)² ≡ x². Thus, we have

$$\begin{eqnarray} {\mathcal N}_E &=& \frac{4\pi \tilde{\phi }_* }{\ln 10} \int _0^{D_*} dD D^2 e^{-D/D_{\rm pp}} \xi ^{-1}\left[ 1 - \left(\frac{L_m}{L_*}\right)^\xi \right]\nonumber\\ &=& \frac{4\pi D_*^3 \tilde{\phi }_* }{\ln 10} \int _0^1 dx \,\xi ^{-1} x^2 (1 - x^{2\xi } ) e^{-x/x_{\rm pp}}. \end{eqnarray} \tag{ A21 }$$

From this we trivially obtain

$$\begin{equation} \frac{d{\mathcal N}_E}{dD} = \frac{1}{D_*} \frac{d{\mathcal N}_E}{dx} = \frac{4\pi D_*^2 \tilde{\phi }_* }{\ln 10} \xi ^{-1} x^2 (1-x^{2\xi }) e^{-x/x_{\rm pp}}, \end{equation} \tag{ A22 }$$

providing some notion of the location of the most numerous sources. Typical values of ξ range from −1 to 0, in practice, and $d\log {\mathcal N}_E/dx$ is shown in the top panel of Figure 12 for a variety of choices of ξ within this range.

Zoom In Zoom Out Reset image size

The integral for $\dot{Q}$ is similarly simple,

$$\begin{eqnarray} \dot{Q}_E &= \int _0^\infty dD 4\pi D^2 e^{-D/D_{\rm pp}} \int _{\log _{10}L_m}^\infty \frac{L}{4\pi D^2 D_{\rm pp}} \tilde{\phi }(L) \,d\!\log _{10}L\nonumber\\ &= \frac{\tilde{\phi }_* L_* D_*}{\ln 10 D_{\rm pp}} \int _0^{1} dx\, \frac{1-x^{2\xi +2}}{1+\xi }\, e^{-x/x_{\rm pp}}, \end{eqnarray} \tag{ A23 }$$

and thus

$$\begin{equation} \frac{d\dot{Q}_E}{dD} = \frac{\tilde{\phi }_* L_*}{\ln 10 D_{\rm pp}}\, \frac{1-x^{2\xi +2}}{1+\xi }\, e^{-x/x_{\rm pp}}, \end{equation} \tag{ A24 }$$

giving an idea of the location of the sources responsible for the bulk of the heating. For a variety of ξ, $d\log \dot{Q}_E/dx$ is shown in the middle panel of Figure 12.

Generally, we find that for all but the largest ξ, $d\log {\mathcal N}_E/dx$ and $d\log \dot{Q}_E/dx$ are peaked at small distances. In particular, both are typically dominated by x < x_pp. The bottom panel of Figure 12 shows the flux-limited L_m with (dashed) and without (solid) absorption included. Including absorption suppresses the contributions at large x/x_pp, forcing $d\log {\mathcal N}_E/dx$ and $d\log \dot{Q}_E/dx$ to be even more strongly peaked at small distances.

Since L_m is a strong function of D, contributions from different distances have different luminosity distributions. It is possible to roughly characterize this by defining a typical luminosity, L_eff, associated with contributions at a given D:

$$\begin{equation} L_{\rm eff} \equiv 4\pi D^2 D_{\rm pp}\frac{d\dot{Q}_E/dD}{d{\mathcal N}_E/dD} = L_* \frac{\xi }{1+\xi } \frac{1-x^{2\xi +2}}{1-x^{2\xi }} \,. \end{equation} \tag{ A25 }$$

With L_m, this is shown in the bottom panel of Figure 12. Generally L_eff is larger than L_m, and for small x substantially so. Thus, even for nearly flat $\tilde{\phi }(L)$ (ξ ∼ −1) the objects that contribute to the heating are not dominated by the numerous, intrinsically dim objects with luminosities L ∼ L_m.

Alternatively, we may perform the integral over D first, in which the flux limit implies a maximum distance to which a given object can be seen, $D_M=\sqrt{L/4\pi F_m}$, providing some insight into luminosity of the sources responsible for the heating. Doing so yields

$$\begin{eqnarray} \frac{d\dot{Q}_E}{d\log _{10}L} &= \frac{D_*}{D_{\rm pp}} L\tilde{\phi }(L) \int _0^{D_M/D_*} e^{-x/x_{\rm pp}} dx \nonumber\\ &= \frac{D_*}{D_{\rm pp}} L\tilde{\phi }(L) \big(1-e^{-\sqrt{L/L_*}/x_{\rm pp}}\big). \end{eqnarray} \tag{ A26 }$$

This is shown in Figure 13 for a variety of ξ. In all cases the luminosity at the break in $\tilde{\phi }(L)$, i.e., L_*, contributes most significantly to the heating rate. However, the relative importance of lower-luminosity objects does depend on the luminosity function; flat luminosity functions (i.e., ξ = −1) have many low-luminosity sources and thus induce more broad $d\log \dot{Q}_E/d\log _{10}L$. Nevertheless, it appears that the heating rate in our simple toy model is generally dominated by nearby objects near the peak in the luminosity function.

Zoom In Zoom Out Reset image size

We now return to the physically relevant case: heating due to TeV blazars with an evolving luminosity function in an evolving universe. Here we specify the distributions of the sources responsible for producing a fraction $\mathcal {Q}$ of the total heating rate. In this case, we may immediately read off $d{\mathcal N}_B/dz$ and $d\dot{Q}/dz$, the analogs of $d{\mathcal N}_E/dD$ and $d\dot{Q}_E/dD$, from Equations (A18) and (A16), respectively, yielding

$$\begin{eqnarray} \frac{d{\mathcal N}_B}{dz}(z;z_o) &=& 4\pi D_A^2(z;z_o) \frac{dD_P}{dz} \tilde{\Phi }_B[z;L_m(z;z_o,\dot{q}_{\mathcal {Q}}),L_M]\nonumber\\ \frac{d\dot{Q}}{dz}(z;z_o) &=& 4\pi D_A^2(z;z_o) \frac{dD_P}{dz} \left(\frac{1+z}{1+z_o}\right)^{2-\alpha }\nonumber\\ &\times&\frac{\chi (z;z_o)}{4\pi D_L^2(z;z_o)} \tilde{\Lambda }_B[z;L_m(z;z_o,\dot{q}_{\mathcal {Q}}),L_M].\qquad \end{eqnarray} \tag{ A27 }$$

These are shown, normalized by their integrated values, for z_o ranging from 0 to 4 in Figure 14 for $\mathcal {Q}=0.5$. In addition, we show an analogously defined characteristic luminosity,

$$\begin{eqnarray} L_{\rm eff}(z;z_o) &\equiv \frac{4\pi D_L^2(z;z_o)}{\chi (z;z_o)} \frac{d\dot{Q}/dz}{d{\mathcal N}_B/dz}\nonumber\\ &= \left(\frac{1+z}{1+z_o}\right)^{2-\alpha } \frac{ \tilde{\Lambda }_B[z;L_m(z;z_o,\dot{q}_{\mathcal {Q}}),L_M] }{ \tilde{\Phi }_B[z;L_m(z;z_o,\dot{q}_{\mathcal {Q}}),L_M] }.\nonumber\\ \end{eqnarray} \tag{ A28 }$$

Zoom In Zoom Out Reset image size

The generic features of our static toy model are also apparent here. At all observer redshifts the heating rate is dominated by the nearest sources. At low z_o, where $L\tilde{\phi }_B$ is nearly flat, the number of objects is also heavily weighted toward nearby objects, well within the redshift at which $\bar{\tau }=1$. However, $d{\mathcal N}_B/dz$ and $d\dot{Q}/dz$ evolve with observer redshift due to both the intrinsic evolution of $\tilde{\phi }_B$ and the background universe. As a consequence, by z_o ∼ 0.5 the peak of $d{\mathcal N}_B/dz$ has moved to the $\bar{\tau }=1$ redshift, implying that z₁ (where $\bar{\tau }(z_1,z_o)=1$) is not a particularly accurate estimate of the redshifts that contribute significantly to the heating. This is further supported by the high-z_o behavior of $d{\mathcal N}_B/dz$, which once again is heavily weighted at redshifts inside of z₁ due to onset of the decline in the blazar population.

The typical luminosities of objects responsible for the heating also evolve. At z_o ∼ 0 these are roughly 10⁴⁴ erg s⁻¹, rising to 3 × 10⁴⁵ erg s⁻¹ by z_o ∼ 2. We also compute the heating rate per logarithmic decade in luminosity:

$$\begin{eqnarray} \frac{d\dot{Q}}{d\log _{10}L} &=& \int _{z_o}^{z_m} 4\pi D_A^2(z;z_o) \frac{dD_P}{dz} \left(\frac{1+z}{1+z_o}\right)^{2-\alpha }\nonumber\\ &&\times \frac{\chi (z;z_o)}{4\pi D_L^2(z;z_o)} L\tilde{\phi }_B(z,L), \end{eqnarray} \tag{ A29 }$$

(where z_m is determined implicitly by $L_m(z_m;z_o,\dot{q}_{\mathcal {Q}})=L$) shown in Figure 6 for a number of z_o. At z_o ∼ 0 the distribution of heating rates is a relatively broad function of L centered near 10⁴⁵ erg s. Until z_o ∼ 2, as z_o grows $d\dot{Q}/d\log _{10}L$ becomes increasingly peaked and centered on increasingly larger luminosities. Above z_o ∼ 2 this trend reverses, though the distribution of luminosities that contribute appreciably to the heating rate never becomes comparable to that in the present epoch. Thus, generally it appears that the heating is due predominantly to nearby objects with luminosities comparable to 10⁴⁵ erg s⁻¹.