PATCHY BLAZAR HEATING: DIVERSIFYING THE THERMAL HISTORY OF THE INTERGALACTIC MEDIUM

The intergalactic medium (IGM) constitutes the bulk of the baryons forming the cosmic web (Bond et al. 1996) as well as the reservoir of baryons available for the formation of galaxies and clusters (Rauch et al. 1997). Observations of metal-enriched gas link its evolution to the formation of galaxies and stars (see, e.g., Simionescu et al. 2009; Werner et al. 2010). Therefore, the thermal history of the IGM plays a central role in determining the development of structure in the visible universe.

The temperature of the IGM is mostly set by the photoionization of hydrogen and helium, competing with adiabatic cooling. As a result, the universe is slowly cooling once reionization is completed. The IGM temperature is expected to increase with density because denser regions experience a lower amount of adiabatic cooling, as well as more recombination-induced photoheating. Gas in the IGM that has not been shock-heated provides a lower envelope to the temperature–density distribution. Observations of Lyα absorption lines with 4 ≤ z ≤ 2 show that lower column density gas exhibits the lowest line widths, which is commonly interpreted as an indication that lower-density gas is colder than high-density gas (Kirkman & Tytler 1997; Ricotti et al. 2000; Schaye et al. 2000; Rudie et al. 2012). A temporary flattening of the slope of the inferred temperature–density relation around z ∼ 3 indicates an additional source of heating at that redshift.

The latter is expected due to He ii reionization (e.g., Tittley & Meiksin 2007; McQuinn et al. 2009; Compostella et al. 2013; Puchwein et al. 2015). Observations suggest that the tail-end of He ii reionization occurred between z ≃ 3.2 and 2.7 (Kriss et al. 2001; Syphers & Shull 2014). However, the bulk of He ii has likely ionized earlier, potentially as early as $z\;\geqslant \;4$ (Worseck et al. 2014). Its patchiness is probably related to the rarity of ionizing sources and the inhomogeneity of the IGM. Including full radiative transfer, simulations show enhanced scatter in the temperature–density distribution around z = 3 (McQuinn et al. 2009). Whether reionization may result in an inverted temperature–density distribution (Meiksin & Tittley 2012) is not firmly established yet (Compostella et al. 2013).

On top of that, recent measurements found that underdense regions may be warmer than predicted (Bolton et al. 2008; Viel et al. 2009), as well as exhibiting unexpectedly high temperatures at z < 2 (Boera et al. 2014), which is harder to explain solely as stemming from the ionization of He ii. Although an inverted temperature–density relation in underdense regions has not been firmly established yet (Bolton et al. 2014), these observations suggest that the thermal history of the IGM may be more complex than initially assumed.

Broderick et al. (2012) recently suggested a complementary heating mechanism, through TeV blazars. TeV blazars are active galactic nuclei (AGNs) emitting very high-energy gamma rays (E ≥ 100 GeV). They belong to the radio-loud subgroup of AGNs that have relativistic jets that are pointed toward us. About 50 of these sources have been significantly discovered so far (http://tevcat.uchicago.edu/) by ground-based Cerenkov telescopes such as MAGIC, H. E. S. S., and VERITAS. Those pointed observations have only skimmed the surface of a much larger population that manifests itself as a hard-spectra gamma-ray blazar, as observed with the space-based Fermi-LAT telescope (Broderick et al. 2014a).

The universe is mainly opaque to very high-energy gamma rays; they interact with the extragalactic background light (EBL), producing electron/positron pairs (Gould & Schréder 1967; Stecker et al. 1992). It is commonly assumed that the electron/positron pairs inverse Compton scatter photons of the cosmic microwave background, resulting in a distribution of photons with energies between 0.1 and 100 GeV. Such an emission component toward TeV blazars has not been observed so far (Aleksić et al. 2010), despite extended searches (H. E. S. S. Collaboration et al. 2014). One solution would be pair deflection due to the intergalactic magnetic field, thus lowering the surface brightness of the formed pair halo at GeV energies (Dermer et al. 2011; Vovk et al. 2012; Durrer & Neronov 2013).

However, the cascaded GeV emission would still contribute to the extragalactic gamma-ray background (EGRB). Fermi-LAT was able to resolve more sources than its predecessor EGRET, thereby limiting the isotropic component of the unresolved EGRB and severely constraining the redshift evolution of hard blazars in the presence of inverse Compton cascades (e.g., Venters 2010; Inoue & Ioka 2012; Murase et al. 2012). As a result of these constraints, it is now well-established that in this picture, the comoving number density of gamma-ray blazars, and by extension the TeV blazars, must either be constant or decreasing with redshift (Kneiske & Mannheim 2008; Venters 2010; Abazajian et al. 2011; Inoue & Ioka 2012). This lies in stark contrast to other classes of AGNs specifically, and all other tracers of the cosmological history of accretion onto galactic halos generally (e.g., star formation). Furthermore, after careful modeling of gamma-ray and X-ray survey selection effects there appears to be no evidence for a significantly different evolution of blazars in comparison to their radio-loud analogs such as radio galaxies (Giommi et al. 2012, 2013). Together, this casts doubt on the inverse Compton cascade picture and provides circumstantial evidence that pair beams instead transfer their energy directly to the IGM through plasma instabilities (Broderick et al. 2012; Schlickeiser et al. 2012, 2013; Chang et al. 2014). This naturally explains the EGRB with a blazar population that exhibits a redshift evolution in agreement with that of quasars (Broderick et al. 2014a, 2014b).

The pairs constitute a dilute, ultrarelativistic beam, which is subject to several plasma instabilities, of which the "oblique" instability (Bret et al. 2004) is the most powerful. Assuming that its efficiency in the linear regime extends to the nonlinear regime, Chang et al. (2012) show that it is responsible for increasing the temperature of the IGM by almost a factor of 10 in low-density regions. While this assumption is still being debated (see Miniati & Elyiv 2013; Sironi & Giannios 2014, and also Schlickeiser et al. 2012, 2013; Chang et al. 2014), throughout this paper we assume  that plasma instabilities are the dominant mechanism for the cooling of the pair beams.

Including TeV blazar heating in the thermal history of the IGM, Chang et al. (2012) were able to reproduce the inverted temperature–density relation for low-density regions. Pfrommer et al. (2012) found that TeV blazar heating is capable of creating a redshift-dependent entropy floor in clusters and galaxies, thus suppressing the formation of dwarf galaxies after the peak of blazar activity at redshift z ≃ 2 and potentially providing an explanation for the "missing satellite problem" and the "missing void dwarf problem" (Kravtsov 2010). Implementing uniform blazar heating, i.e., with a redshift-dependent but spatially homogeneous energy deposition rate per unit volume, in a cosmological hydrodynamical simulation of galaxy formation, Puchwein et al. (2012) find excellent agreement with the one- and two-point statistics of the Lyα forest, which is the main observational tracer of low-density regions in the universe.

However, the low thermal broadening of certain Lyα lines indicates the presence of cold gas at z = 2.4 (Rudie et al. 2012), which suggests that TeV blazar heating does not uniformly heat the whole universe. It is natural to expect a larger TeV flux close to higher-density regions where visible structures form. Conversely, in large underdense regions, far from massive black holes, heating is probably much lower. The goal of this paper is to go beyond the hypothesis of uniform heating and to link TeV blazar heating to the underlying clustered density field and take into account the bias of sources. This will lead to a more heterogeneous heating pattern and account for unheated regions while keeping the overall impact of blazar heating.

Self-consistently studying the evolution of the IGM from first principles involves modeling both the formation and evolution of galaxies at the largest scales of the universe. As this is still far beyond reach of current computers, we have determined a filter function that relates the heating fluctuations to the dark matter (DM) structure, similar to Pritchard & Furlanetto (2007), Barkana & Loeb (2005), and Pontzen (2014). Based on the hierarchical structure formation in a ΛCDM universe, it naturally selects the relevant length scales for TeV blazar heating (Section 2). We implement our blazar heating model in large-scale cosmological simulations (Section 3) in order to focus on the equation of state and thermal evolution of the IGM (Section 4). We then discuss how inhomogeneous heating could reconcile different observations (Section 5) and draw our conclusions (Section 6).

2.1. Intuitive Understanding

One zone models (Chang et al. 2012; Pfrommer et al. 2012) and numerical simulations (Puchwein et al. 2012) of blazar heating on the IGM assume that the heating is uniform. Because the heating rate depends on the local density of EBL and TeV photons, the assumption of uniform heating implies that the distributions of EBL and TeV photons are uniform. The EBL photon density in voids is only 2% less than the average value (Furniss et al. 2015). However, for TeV photons the mean free path compares with the separation between TeV sources for z ≥ 1, so the spatial fluctuations in the heating rate are likely non-trivial. Moreover, the sources of TeV photons tend to be clustered, so the IGM near these clustered regions will get an increased flux of photons in comparison to low-density regions because of the increased number of sources and the 1/r² flux dilution with increasing distance r from the source.

Our goal is to include a more realistic model for heating due to TeV blazars in numerical simulations. To properly calculate the heating fluctuations due to TeV blazar heating, the formation and evolution of accreting supermassive black holes must be modeled in a full self-consistent cosmological simulation. In addition, the TeV radiation from these systems must be ray-traced through the simulation volume. Such a task is computationally intractable. As a result, we have elected to model this TeV blazar heating in a statistical manner.

We assume that TeV blazars are associated with galaxies or alternatively quasars and that they roughly emit over 4π steradian. The latter assumption remains valid as long as the duty cycle of blazars is small enough such that jets point in all directions over cosmological timescales and the number of TeV blazars is large enough such that every spot in the universe is illuminated by at least a few TeV blazars, which is the case for z ≲ 6 (Chang et al. 2012). The heating rate at a given point ${\boldsymbol{x}}$ is determined by the received TeV flux from all the sources

$$\begin{eqnarray}&&\dot{Q}({\boldsymbol{x}},z)=\displaystyle \frac{1}{4\pi }{\displaystyle \int }_{{\rm{\Omega }}}d{\rm{\Omega }}{\displaystyle \int }_{0}^{\infty }{dr}^{\prime} \displaystyle \frac{{\mathcal{E}}({{\boldsymbol{r}}}^{\prime }+{\boldsymbol{x}},z^{\prime} )}{{D}_{\mathrm{pp}}(z)}{e}^{-\tau },\end{eqnarray} \tag{ 1 }$$

where $\dot{Q}$ and ${\mathcal{E}}$ are the heating rate density and emissivity of the sources (both in units of energy per unit time and per unit volume), ${D}_{\mathrm{pp}}$ is the mean free path for pair production on the EBL, τ is the associated optical depth along the line of sight, and z' is the redshift corresponding to the distance $| {{\boldsymbol{r}}}^{\prime }|$.

One can then express the resulting heating rate as a mean value $\bar{\dot{Q}}$ and a first order correction δ_H.

$$\begin{eqnarray}&&\dot{Q}({\boldsymbol{x}},z)=\bar{\dot{Q}}({\boldsymbol{x}},z)[1+{\delta }_{{\rm{H}}}({\boldsymbol{x}},z)].\end{eqnarray} \tag{ 2 }$$

The method is based on a Taylor expansion of the quantities describing the TeV sources and keeping only the first order corrections. Transforming to Fourier space yields the fluctuation amplitude

$$\begin{eqnarray}&&{\tilde{\delta }}_{{\rm{H}}}(k,z)={\tilde{W}}_{{\rm{H}}}(k,z)\tilde{\delta }(k,z),\end{eqnarray} \tag{ 3 }$$

where ${\tilde{W}}_{{\rm{H}}}$ is defined as the window function and maps the Fourier transform of the DM overdensity, $\tilde{\delta }$, to the Fourier transform of the heating fluctuations, ${\tilde{\delta }}_{{\rm{H}}}$. This naturally yields the relevant length scale for heating rate fluctuations. The detailed and pedagogical exposition of this method (in the Newtonian limit, in an expanding universe, and additionally accounting for a clustered source distribution) is given in the Appendices A–C. In the following section, we present the general method and highlight the underlying hypotheses of our work.

2.2. Window Function for TeV Blazar Heating

To determine the heating rate fluctuations we express the TeV emissivity in Equation (1) as a mean value and a first order correction. The heating rate at a given point is set by the received TeV flux from all the sources within a certain radius. We assume that the pairs lose their energy to the IGM at the point where they are created. As stated in Broderick et al. (2012) and Chang et al. (2014), this is a reasonable assumption, as the plasma instability length scales are many orders of magnitude smaller than the mean free path of these TeV photons.

The TeV gamma rays are emitted by accreting supermassive black holes at centers of galaxies, which cluster in overdense regions. Matter is tightly coupled to the underlying DM, the evolution of which is straightforward to model analytically within the linear approximation. The linear approximation is valid as long as the overdensity is small, which is true in the early universe and then breaks down at small scales as very dense structures form. Our computation takes into account the bias between baryonic matter and DM (Mo & White 1996), as we detail below.

To model cosmic distances, Equation (1) is integrated in redshift space and we take into account the resulting energy loss for the TeV photons as well as first order corrections due to a clustered blazar distribution and proper motions of the sources within the Hubble flow (Kaiser 1987). We integrate over the energy distribution of the TeV emission. Following the detailed derivation in the appendices, we obtain the window function ${\tilde{W}}_{{\rm{H}}}(\hat{k},z)$ as a function of comoving wavelength $\hat{k}$,

$$\begin{eqnarray}{\tilde{W}}_{{\rm{H}}}(\hat{k},z) & = & \displaystyle \frac{1}{N}{\displaystyle \int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}\displaystyle \frac{{dE}}{{D}_{\mathrm{pp}}(E,z)}{\displaystyle \int }_{z}^{\infty }{dz}^{\prime} \displaystyle \frac{{\bar{\hat{{\mathcal{E}}}}}_{{{\rm{E}}}^{\prime }}(z^{\prime} ){e}^{-\tau }}{(1+z^{\prime} )H(z^{\prime} )}\\ & & \times \displaystyle \frac{D(z^{\prime} )}{D(z)}[(b(z^{\prime} )+\displaystyle \frac{f}{3}){j}_{0}(\hat{k}\hat{r})-\displaystyle \frac{2f}{3}{j}_{2}(\hat{k}\hat{r})],\end{eqnarray} \tag{ 4 }$$

with

$$\begin{eqnarray}&&N={\displaystyle \int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}\displaystyle \frac{{dE}}{{D}_{\mathrm{pp}}(E,z)}{\displaystyle \int }_{z}^{\infty }{dz}^{\prime} \displaystyle \frac{{\bar{\hat{{\mathcal{E}}}}}_{{{\rm{E}}}^{\prime }}(z^{\prime} ){e}^{-\tau }}{(1+z^{\prime} )H(z^{\prime} )}.\end{eqnarray} \tag{ 5 }$$

Here, $\hat{r}$ is the comoving distance between the source and heated region, D denotes the linear growth factor defined in Equation (33), $f\equiv d\mathrm{log}\delta /d\mathrm{log}a$, E is the energy of the received TeV photon, E' is its initial energy, and ${\hat{{\mathcal{E}}}}_{{\rm{E}}}$ is the comoving blazar spectral luminosity density, b is the Eulerian bias, and j₀ and j₂ are spherical Bessel functions.

The blazar spectral luminosity density (in units of energy, per unit time, per volume, per energy) is determined by both the intrinsic source spectra and the luminosity function of TeV blazars. Here we adopt the model described in Broderick et al. (2014a), in which the spectra are assumed to be well described by a set of broken power laws and a luminosity function that depends on redshift, luminosity from 100 GeV to 10 TeV (L_TeV), and the low-energy spectral index (Γ_l). The assumed source spectral luminosity density is given by ${L}_{{\rm{E}}}({L}_{\mathrm{TeV}},{{\rm{\Gamma }}}_{{\rm{l}}})={\tilde{L}}_{{\rm{E}}}({{\rm{\Gamma }}}_{{\rm{l}}}){L}_{\mathrm{TeV}}$, where

$$\begin{eqnarray}&&{\tilde{L}}_{{\rm{E}}}({{\rm{\Gamma }}}_{{\rm{l}}})=\displaystyle \frac{{f}_{0}E}{{(E/{E}_{b})}^{{{\rm{\Gamma }}}_{{\rm{l}}}}+{(E/{E}_{b})}^{{{\rm{\Gamma }}}_{{\rm{h}}}}},\end{eqnarray} \tag{ 6 }$$

where the normalization f₀ is chosen such that ${\displaystyle \int }_{0.1\;\mathrm{GeV}}^{10\;\mathrm{TeV}}{\tilde{L}}_{{\rm{E}}}{dE}=1$, E_b = 1 TeV, and Γ_h = 3. The luminosity function is based on the spectral properties of the Fermi population of hard gamma-ray blazars and the quasar luminosity function reported by Hopkins et al. (2007). In terms of these, the blazar spectral luminosity density is

$$\begin{eqnarray}&&{\hat{{\mathcal{E}}}}_{{\rm{E}}}(z)=\int d{\mathrm{log}}_{10}{L}_{\mathrm{TeV}}\;d{{\rm{\Gamma }}}_{{\rm{l}}}\;{L}_{{\rm{E}}}({L}_{\mathrm{TeV}},{{\rm{\Gamma }}}_{{\rm{l}}})\hat{\phi }(z,{L}_{\mathrm{TeV}},{{\rm{\Gamma }}}_{{\rm{l}}}),\end{eqnarray} \tag{ 7 }$$

where $\hat{\phi }(z,{L}_{\mathrm{TeV}},{{\rm{\Gamma }}}_{{\rm{l}}})$ is the comoving analog of Equation (6) in Broderick et al. (2014a). This may be simplified further using the separability of $\hat{\phi }$ in Γ_l, i.e.,

$$\begin{eqnarray}&&\hat{\phi }(z,{L}_{\mathrm{TeV}},{{\rm{\Gamma }}}_{{\rm{l}}})=\hat{\phi }(z,{L}_{\mathrm{TeV}})\chi ({{\rm{\Gamma }}}_{{\rm{l}}}),\end{eqnarray} \tag{ 8 }$$

where $\int d{{\rm{\Gamma }}}_{{\rm{l}}}\chi ({{\rm{\Gamma }}}_{{\rm{l}}})=1$. As a result,

$$\begin{eqnarray}&&{\hat{{\mathcal{E}}}}_{{\rm{E}}}(z)=\langle {\tilde{L}}_{{\rm{E}}}\rangle \hat{{\rm{\Lambda }}}(z),\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&\langle {\tilde{L}}_{{\rm{E}}}\rangle =\int d{{\rm{\Gamma }}}_{{\rm{l}}}\;{\tilde{L}}_{{\rm{E}}}({{\rm{\Gamma }}}_{{\rm{l}}})\chi ({{\rm{\Gamma }}}_{{\rm{l}}})\end{eqnarray} \tag{ 10 }$$

which is now a function of energy alone, and

$$\begin{eqnarray}&&\hat{{\rm{\Lambda }}}(z)=\int d{\mathrm{log}}_{10}{L}_{\mathrm{TeV}}\;\hat{\phi }(z,{L}_{\mathrm{TeV}})\end{eqnarray} \tag{ 11 }$$

is the TeV luminosity density of TeV blazars.

In practice, for the TeV blazar luminosity function in Broderick et al. (2014a) $\langle {\tilde{L}}_{{\rm{E}}}\rangle$ is well fit over the energy range of interest, 100 GeV–10 TeV, by Equation (6) with Γ_l = 1.80 and E_b = 1.057, with an L₁ norm of 0.0048. This spectral shape peaks near 1 TeV and thus sets a characteristic gamma-ray energy.

Since the TeV blazar luminosity function in Broderick et al. (2014a) is simply the quasar luminosity function in Hopkins et al. (2007) rescaled in energy and luminosity, the comoving luminosity densities are also proportional. That is,

$$\begin{eqnarray}&&\hat{{\rm{\Lambda }}}(z)=\zeta {\hat{{\rm{\Lambda }}}}_{{\rm{Q}}}(z),\end{eqnarray} \tag{ 12 }$$

with ζ = 2.1 × 10⁻³ and the comoving quasar luminosity density ${\hat{{\rm{\Lambda }}}}_{{\rm{Q}}}$. The corresponding comoving TeV luminosity density, found by Chang et al. (2012) based on a fit to the comoving quasar luminosity density in Hopkins et al. (2007), is given by

$$\begin{eqnarray}&&\hat{{\rm{\Lambda }}}(z)={\hat{{\rm{\Lambda }}}}_{0}{10}^{1.18z-0.0812{z}^{2}-0.182{z}^{3}+0.0515{z}^{4}-0.00418{z}^{5}}\end{eqnarray} \tag{ 13 }$$

where ${\hat{{\rm{\Lambda }}}}_{0}=(1.7-4.8)\times {10}^{-36}\;\mathrm{erg}\;{{\rm{s}}}^{-1}\;{\mathrm{cm}}^{-3}$ is the blazar luminosity density at the current epoch.

The optical depth is set by the physical mean free path ${D}_{\mathrm{pp}}(E^{\prime} ,z)$ of TeV photons before they interact with an EBL photon to produce an electron–positron pair. Its redshift evolution is set by the density of the EBL. Despite careful studies that aim at constraining it, there remain uncertainties in the star formation history of the universe, as well as its metallicity and dust contents (see, e.g., Stecker et al. 2006; Franceschini et al. 2008). Following Chang et al. (2012) we use a prescription

$$\begin{eqnarray}&&{D}_{\mathrm{pp}}(E^{\prime} ,z^{\prime} )=35{(\displaystyle \frac{E^{\prime} }{1\;\mathrm{TeV}})}^{-1}{(\displaystyle \frac{1+z}{2})}^{-\xi }\;\mathrm{Mpc},\end{eqnarray} \tag{ 14 }$$

where ξ = 4.5 for z < 1 and ξ = 0 for z > 1 (Kneiske et al. 2004; Neronov & Semikoz 2009) and the mean free path is expressed in physical units. The proper mean free path is constant for z ≥ 1, however, the comoving mean free path increases for increasing redshift.

The fluctuations in TeV flux are related to the distribution of blazars, which is biased with respect to the distribution of DM halos. Luminous structures such as galaxies preferentially populate the high peaks of the DM density distribution. The square of the bias b of a given structure is the ratio between its power spectrum to the power spectrum of the DM halos and it is stronger for more massive objects, such as the host galaxies of quasars (see, e.g., Cooray & Sheth 2002 for a review).

Due to the small number of sources and TeV photon absorption, there is no observation of TeV blazar bias. Therefore, we use a model for quasar bias as well as a model for galaxy bias to illustrate the impact of the uncertainty on the value of the bias. Since AGNs are generally much more biased than galaxies, the latter provide a robust lower bound. However, TeV blazars may have a stronger bias than quasars, as radio-loud AGN are generally found in a more clustered environment than quasars (Mandelbaum et al. 2009; Simpson et al. 2012). Our model for quasar bias for z > 1 is based on a fit of data by Croom et al. (2005), Myers et al. (2007), and Shen et al. (2007). We use

$$\begin{eqnarray}&&{b}_{\mathrm{quasar}}(z)={10}^{0.27z-0.04},\end{eqnarray} \tag{ 15 }$$

which is very similar to Papageorgiou et al. (2012). Our model for galaxy bias is based on a fit of data by Marinoni et al. (2005), Steidel et al. (1998), and Kashikawa et al. (2006). We use

$$\begin{eqnarray}&&{b}_{\mathrm{galaxy}}(z)={10}^{0.174z}.\end{eqnarray} \tag{ 16 }$$

Figure 1 shows the corresponding evolution of quasar and galaxy bias as a function of redshift. Using the same data, Basilakos et al. (2008) found that the quasar and galaxy bias can be fitted using DM halo models of 10¹³ h⁻¹ M_⊙ and 10¹² h⁻¹ M_⊙, respectively.

Zoom In Zoom Out Reset image size

Substituting Equations (9) and (14) into Equation (4) then gives the complete window function for TeV blazar heating. Figure 2 shows the filter for z = 1, 2, and 4 for a model with galaxy and quasar bias. We have computed the window function using an embedded Runge–Kutta method that is able to capture the fast variation of the Bessel functions at large wave numbers while decreasing computing time at smaller wave numbers. The window function is integrated up to z = 8, which ensures that all of the sources are taken into account. Changes in the maximum redshift of the integration have no discernible impact on the window function.

Zoom In Zoom Out Reset image size

The window function describes how density fluctuations translate into heating fluctuations. The ${\hat{k}}^{-1}$ slope for large wavenumbers results from the r⁻² decrease in the TeV flux. The position of the break is set by a combination between the mean free path and the redshift evolution of the blazar luminosity density. One might expect that the comoving mean free path, which increases for increasing redshift at z > 1, would move the break to smaller wavenumbers at higher z. However, by integrating over the energy of the photons and including lower-energy photons that have longer mean free paths, the effect of the change in the mean free path is mitigated. This allows the rapidly decreasing blazar luminosity density for z ≳ 2 to move the break to a larger wavenumber: the luminosity density drops too steeply so that more distant sources within the photon mean free path do not contribute as much as sources that are closer. This latter effect moves the break to a larger wavenumber for increasing redshift when lower-energy photons with longer mean free paths are accounted for. This implicitly assumes that these low-energy photons would contribute to the local heating rate via pair production and plasma instabilities. A more rigorous calculation that accounts for the efficiency of blazar heating would vary with photon energy and redshift (see, for instance, Chang et al. 2014 for the effect of nonlinear Landau damping). However, this efficiency is not yet fully understood, so in this paper we assume that it is unity as a function of photon energy. As the normalization is set by the bias, the quasar bias model has more power at all scales and density fluctuations will lead to more enhanced heating. In both models, at high redshift, most of the power resides on large scales, where blazar heating traces the density fluctuations. At smaller scales, density fluctuations have no impact and blazar heating is uniform. At the current epoch, TeV blazar heating is close to uniform because there is power only at the largest scales (above 100 Mpc) where the universe is essentially uniform (see Clowes et al. 2013 and the references therein).

In both models, the window function remains positive at all scales; thus, underdense regions are the only areas where a lower than average heating rate is expected. Heating at rates larger than the average is possible for large-scale overdensities but also on small scales for highly nonlinear objects with $\delta \gg 1$. We caution that our simplified approach with a linear bias description starts to break down there and would have to be augmented with a nonlinear description of bias. However, the gas in these regions experiences shock heating in addition to photoheating such that the influence of blazar heating becomes negligible there. Moreover, those densities have little influence on the statistics of the high-redshift Lyα forest and thus shall not be the subject of the present work. After these first analytic estimates, we include the window function to model TeV blazar heating in cosmological simulations.

3.1. Cosmological Simulations

3.2. Including the TeV Blazar Heating Fluctuations

Figure 3 shows the heating rate fluctuations in the midplane of the L_x = 100 h⁻¹ Mpc simulation for z = 3 and 1 for both the galaxy and quasar bias model. Both simulation models are started from identical initial conditions such that any visible difference is solely due to the different bias assumptions in the blazar heating model. The corresponding density field is shown in the upper column and shows increasing structure formation as the redshift decreases. The heating map has a linear scale while the density scale is logarithmic. The heating rate fluctuations are on average much smaller and more spatially homogeneous than the density fluctuations. This is because the window function filters out small scales, which correspond to collapsed regions, where density fluctuations are the highest. To the zeroth order, one can thus consider TeV blazar heating to be uniform, as was assumed in Chang et al. (2012).

Zoom In Zoom Out Reset image size

Additional heating (i.e., δ_H > 0) occurs around clustered regions. This is expected, as the window function translates large-scale density fluctuations into heating rate fluctuations. Conversely, underdense regions such as the one around {x = 60, y = 70} (for z = 3) display heating that is below average, as they are isolated from sources, and their flux decreases as ${e}^{-\tau }{r}^{-2}$. This is most obvious at high redshift, and for the quasar model, where bias is the strongest. As the redshift decreases, heating rate fluctuations decrease because of the decreasing bias.

Figure 4 shows the volume-weighted ratio of the internal energy when blazar heating is included to the internal energy when blazar heating is not included as a function of the density. For both simulation models, we divided the simulation volume into equally sized small sub-volumes and computed the internal energy therein. Because both simulations have been started from identical initial conditions, the interference of the primordial density waves gives rise to a comparable morphology of the cosmic web so that we can compare thermodynamic quantities of (almost) identical sub-volumes of the different simulations. These maps clearly highlight that blazar heating has more impact in underdense regions, as the heating rate per baryon is higher. Even if these regions receive less heat than regions with higher density (see Figure 3), blazar heating increases the internal energy by a factor of two at z = 3 and about an order of magnitude at z = 1. At low redshift, even in this inhomogeneous model, most of the gas is heated up and there is minimal difference between the two models we used for the bias.

Zoom In Zoom Out Reset image size

The impact of inhomogeneous heating translates into a more complex temperature–density relation, as is shown on Figure 5. The color map shows the mass-weighted $T-{\rho }_{\mathrm{gas}}$ relation from our simulations and the gray contours show the case for uniform blazar heating with the same resolution (Puchwein et al. 2012). Temperature measures the integrated impact of TeV blazar heating over time. The left column shows the $T-{\rho }_{\mathrm{gas}}$ relation in a model with no blazar heating. In the latter case, the underdense gas follows a very narrow distribution, where the lowest density gas is the coldest. When blazar clustering is taken into account, the temperature–density relation has a significant scatter for underdense regions for z ≃ 2–3. For the quasar bias model, this scatter results in the lower envelope of the temperature–density that differs little from the case with no blazar heating. However, the mode of the temperature is very close to the uniform blazar heating case. At z = 1 both models show a very similar behavior and blazar heating can be considered nearly homogeneous, though the lower envelope sits at a lower temperature.

Zoom In Zoom Out Reset image size

Figure 6 shows the volume-weighted probability distribution functions of the temperature for all of the simulations. This provides a more quantitative view of the scatter in temperature and highlights the impact of clustering on TeV blazar heating. The simulations with inhomogeneous blazar heating show significant deviation from the uniform case for z ≥ 2. The higher value of the quasar bias results in the presence of more unheated gas than in the galaxy bias model. Here the effect of the lower envelope is clear, as there is a significant tail toward lower temperatures: the coldest zones have T ≤ 5 × 10³ K, while the warmest zones have T ≃ 10⁵ K for z = 3. Conversely, the warmest gas is only slightly warmer than the mode of the temperature. At lower redshift, the mode of the temperature is very similar in all models but the inhomogeneous models exhibit a larger amount of colder gas.

Zoom In Zoom Out Reset image size

To have a better understanding of the heating fluctuations with respect to density fluctuations we show the mass-weighted ${\delta }_{{\rm{H}}}-{\delta }_{\mathrm{gas}}$ distribution on Figure 7. The heating rate represents an instantaneous view of the impact of TeV blazar heating as opposed to temperature or internal energy, which probe the integrated injection history of non-gravitational heat. As in Figure 3, the quasar bias model stands out at high redshift (lower left panel). In this model, most of the particles receive slightly more heat than in the uniform case. In contrast, the additional heat is only a few times more than the uniform case, while certain regions receive orders of magnitude less heat than the mean value. These areas suffer from the decrease of the TeV flux and isolation from massive structures. This is consistent with the temperature probability distribution function presenting a prominent low temperature tail and a low probability for high temperatures for z ≥ 2. In the galaxy model, most of the gas is heated similarly to the uniform case, translating the lower bias for galaxies.

Zoom In Zoom Out Reset image size

Previous work on TeV blazar heating showed a significant increase in the temperature of the low-density IGM (Chang et al. 2012; Puchwein et al. 2012). These results, while in good agreement with observational data for the mean transmitted flux statistics as well as fits to individual lines, appear to be in potential tension with recent work by Rudie et al. (2012). This observation suggests that a significant presence of gas in the IGM has not been exposed to additional heating beyond photoheating. In the context of blazar heating, this suggests that the heating is not entirely uniform. In this work, we have studied the impact of inhomogeneous heating and find that when accounting for the clustering of sources, the variability in the heating rate is significant. In particular, we find that while most of the gas receives close to average heating, and has a temperature close to the uniform heating case, there is an important scatter in the temperature of the gas for z ≳ 2.

The lower envelope of the $T-\rho$ distribution we find is within the error bars of the power-law fit to the $T-\rho$ distribution derived in Bolton et al. (2014). Detailed comparison of the Lyα forest properties will be the subject of a forthcoming paper. If confirmed with the analysis of the resulting Lyα forest, our model would be able to account for both the observed increase in temperature in the low-density IGM (Viel et al. 2009; Becker et al. 2011; Boera et al. 2014) and the lower values of the line widths as suggested by Rudie et al. (2012). Taken together, both observational diagnostics appear to point to a more complicated temperature–density relation with a possibly broad, multi-valued distribution as opposed to a simple power-law relation.

Around 4 ≲ z ≲ 2.7, the impact of TeV blazar heating is qualitatively similar to the impact of late reionization of He ii. Various physical models and assumptions for the quasar distribution can result in highly variable temperature–density distributions in models of He ii reionization (Bolton et al. 2004; McQuinn et al. 2009). However, comparing the temperature–density distributions obtained from such models with our Figure 5, we notice two major differences. While our model and reionization models present a similar lower envelope of ≃3 × 10³ K at $\rho /\bar{\rho }=0.1$ (at z = 3), our model includes a significant fraction of gas above a few 10⁴ K, which is not present in reionization models. On top of this increased scatter, our model presents a warmer mean temperature than all the He ii models. The same comparison holds at z = 4.

Our model also predicts a warmer IGM at z = 2, which would not be expected from He ii reionization, as the dependence on the exact He ii reionization history will have largely faded at that time (Compostella et al. 2013). An increased temperature with a strong scatter in underdense regions could then be attributed to inhomogeneous blazar heating. Direct measurements of the low-density IGM are hard to obtain, as the Lyα forest is sensitive to overdensities of at least a few. The Lyα forest of He ii traces lower-density regions and is a promising observable, with more data becoming available (Worseck et al. 2014).

Our model neglects important physical aspects of TeV blazars such as their duty cycle and beaming of the gamma-ray emission. The gamma-ray duty cycle of BL Lac objects is of the order of 10% (Stecker & Salamon 1996). As TeV blazars have low accretion rates, their lifetime is long ≃5–7 Gyr (Cavaliere & D'Elia 2002). The gamma-ray sources are highly beamed, with observed opening angles peaking around 20° (Pushkarev et al. 2009). Parsec-scale, long-baseline interferometry observations of the inner parts of radio jets indicate variations of the position angle of around 20°, up to 120° over more than a decade (Lister et al. 2013). Assuming the gamma-ray emission follows the orientation of the radio jet, the rapid variability of the sources can potentially increase the scatter in the heating rate. Indeed, in the early universe, when few sources are present, individual variability increases the inhomogeneity of blazar heating. In such cases, our model is a lower limit for the true scatter. However, later on, as the number of sources at a given point increases, the variability of the heating rate will be reduced. On top of the intrinsic variability of the sources, shot noise, which becomes important when sources are rare, is not taken into account and will result in additional heating rate fluctuations. For instance, at z ∼ 1–2, the number of sources that contributes half the heating is ≈10 whereas at z = 3 it dwindles to ≈1 (Chang et al. 2012). While this would suggest that the shot noise is potentially large, this is mitigated somewhat by the fact that the next 25% of the heating is provided by ≈100 sources between z = 1–2 and ≈10 source at z = 3. Hence, the fluctuations that we calculate in this work are likely a lower limit.

Our model is based on quasar bias, which is likely lower than TeV blazar bias (Allevato et al. 2014). Radio loud AGNs are associated with red giant elliptical galaxies (Hyvönen et al. 2007) and are often at the center of small clusters and groups. Quasars have a wider distribution of host galaxies and are rarely found at the center of clusters. Therefore, the scatter found in our simulations at low redshift is probably a lower limit to the exact scatter. In our simulations, any dense source emits TeV photons, regardless of the type of galaxy, and the history of mergers, accretion, or star formation. However, based on the comparison of the galaxy and quasar bias models presented here with the uniform heating case, we expect a more detailed physical model to produce comparable results.

In this paper we have implemented inhomogeneous TeV blazar heating into cosmological simulations to study the impact of clustering on the thermal state of the IGM. This extends the work by Chang et al. (2012) and Puchwein et al. (2012), based on uniform blazar heating.

We developed a filtering function relating heating rate fluctuations to the linear DM fluctuations. Using this window function, we are able to model the relevant length scales for the blazar heating fluctuations and include them in a statistical fashion. Our method is a cost-effective alternative method to fully self -consistent simulations of blazar heating where modeling black hole formation and growth and complete radiative transfer would have been prohibitive.

Our model for the blazar heating fluctuations yields a temperature–density relation that can potentially reconcile the observed increased mean temperature of the IGM (Boera et al. 2014) while maintaining a fraction of the cold gas responsible for the lower envelope of the line width distribution (Rudie et al. 2012). Therefore, detailed modeling of the Lyα forest will be the subject of a forthcoming paper. If confirmed, this will clearly indicate a more complex thermal history of the IGM, potentially with an important impact on late forming structures.

A.L. and P.C. are supported by the UWM Research Growth Initiative, the NASA ATP program through NASA grant NNX13AH43G, and NSF grant AST-1255469. A.E B. and M.S. receive financial support from the Perimeter Institute for Theoretical Physics and the Natural Sciences and Engineering Research Council of Canada through a Discovery Grant. Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Research and Innovation. C.P. gratefully acknowledges the financial support of the Klaus Tschira Foundation. E.P. acknowledges support by the ERC grant "The Emergence of Structure during the epoch of Reionization." The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin and the NASA Advanced Supercomputing Division for providing HPC resources that contributed to the research results reported within this paper. The authors thank S. Furlanetto, A. Loeb, and M. Heahnelt for fruitful discussions.

We detail the derivation of the window function of Equation (4). To present a transparent derivation of the window function (and to clarify some confusion in the literature), we start with a simple toy model to which we add progressively more physics. A familiarized reader may directly switch to Appendix C. In Appendix A, we start from a purely Newtonian universe and assume that heating rate fluctuations perfectly trace density fluctuations and a monochromatic source of very high gamma rays. Then in Appendix B, we include the impact of cosmological expansion and a spectral energy distribution of TeV blazars that evolve with redshift as quasars, albeit with a significantly smaller normalization in the luminosity density. This impacts on the energy-dependent mean free path of TeV photons. Finally, in Appendix C we also account for the clustering of sources (which we model through the linear bias factor) and redshift distortions as a result of peculiar infall velocities onto cluster potentials.

A.1. Fluctuations with Respect to the Mean Heating Rate

The TeV flux received (in erg s⁻¹ cm⁻²) at position ${\boldsymbol{x}}$, at a given energy, is given by the sum over all the sources

$$\begin{eqnarray}J({\boldsymbol{x}}) & = & {\displaystyle \int }_{0}^{2\pi }d\phi {\displaystyle \int }_{0}^{\pi }\mathrm{sin}\theta d\theta \times {\displaystyle \int }_{0}^{\infty }{| {{\boldsymbol{x}}}^{\prime }-{\boldsymbol{x}}| }^{2}d| {{\boldsymbol{x}}}^{\prime }-{\boldsymbol{x}}| \displaystyle \frac{{\mathcal{E}}({{\boldsymbol{x}}}^{\prime })}{4\pi | {{\boldsymbol{x}}}^{\prime }-{\boldsymbol{x}}{| }^{2}}{e}^{-\tau }= & {\displaystyle \int }_{{\rm{\Omega }}}d{\rm{\Omega }}{\displaystyle \int }_{0}^{\infty }{dr}\displaystyle \frac{{\mathcal{E}}({\boldsymbol{r}}+{\boldsymbol{x}})}{4\pi }{e}^{-r/{D}_{\mathrm{pp}}},\end{eqnarray} \tag{ 17 }$$

where the emissivity ${\mathcal{E}}$ is an energy per unit time, per unit volume, $\tau =| {{\boldsymbol{x}}}^{\prime }-{\boldsymbol{x}}| /{D}_{\mathrm{pp}}$ is the optical depth along the line of sight, and ${D}_{\mathrm{pp}}$ is the mean free path for pair production, which is constant over the volume. In the last step, we introduced ${\boldsymbol{r}}={{\boldsymbol{x}}}^{\prime }-{\boldsymbol{x}}$, $d{\rm{\Omega }}=\mathrm{sin}\theta d\theta d\phi$. We have further assumed that emissivity is monochromatic such that:

$$\begin{eqnarray}&&{{\mathcal{E}}}_{{\rm{E}}}({\boldsymbol{x}})={\mathcal{E}}({\boldsymbol{x}})\delta (E-{E}_{0}),\end{eqnarray} \tag{ 18 }$$

where δ is the Dirac delta function and E₀ is the energy of the monochromatic photons. Similarly, the spectral flux is

$$\begin{eqnarray}&&{J}_{{\rm{E}}}({\boldsymbol{x}})=J({\boldsymbol{x}})\delta (E-{E}_{0}).\end{eqnarray} \tag{ 19 }$$

For simplicity, we have integrated Equation (17) over the monochromatic source.

The differential heating rate is defined as $d\dot{Q}({\boldsymbol{x}})\equiv {dJ}({\boldsymbol{x}})/{D}_{\mathrm{pp}}$ such that the volume-integrated heating rate $\dot{Q}$ and its average (in units of erg cm⁻³ s⁻¹) are given by

$$\begin{eqnarray}&&\dot{Q}({\boldsymbol{x}})=\displaystyle \frac{1}{4\pi }{\displaystyle \int }_{{\rm{\Omega }}}d{\rm{\Omega }}{\displaystyle \int }_{0}^{\infty }{dr}{\mathcal{E}}({\boldsymbol{r}}+{\boldsymbol{x}})\displaystyle \frac{{e}^{-r/{D}_{\mathrm{pp}}}}{{D}_{\mathrm{pp}}},\quad {\rm{and}}\quad \bar{\dot{Q}}=\displaystyle \frac{1}{4\pi }{\displaystyle \int }_{{\rm{\Omega }}}d{\rm{\Omega }}{\displaystyle \int }_{0}^{\infty }{dr}\bar{{\mathcal{E}}}\;\displaystyle \frac{{e}^{-r/{D}_{\mathrm{pp}}}}{{D}_{\mathrm{pp}}}=\bar{{\mathcal{E}}}.\end{eqnarray} \tag{ 20 }$$

The heating rate fluctuations at a given point ${\boldsymbol{x}}$ are then given by

$$\begin{eqnarray}&&{\delta }_{{\rm{H}}}({\boldsymbol{x}})=\displaystyle \frac{\dot{Q}({\boldsymbol{x}})-\bar{\dot{Q}}}{\bar{\dot{Q}}}=\displaystyle \frac{1}{4\pi \bar{\dot{Q}}{D}_{\mathrm{pp}}}{\displaystyle \int }_{{\rm{\Omega }}}d{\rm{\Omega }}{\displaystyle \int }_{0}^{\infty }{dr}[{\mathcal{E}}({\boldsymbol{r}}+{\boldsymbol{x}})-\bar{{\mathcal{E}}}]{e}^{-r/{D}_{\mathrm{pp}}}=\displaystyle \frac{1}{4\pi {D}_{\mathrm{pp}}}{\displaystyle \int }_{{\rm{\Omega }}}d{\rm{\Omega }}{\displaystyle \int }_{0}^{\infty }{dr}{\delta }_{{\rm{E}}}({\boldsymbol{r}}+{\boldsymbol{x}}){e}^{-r/{D}_{\mathrm{pp}}},\end{eqnarray} \tag{ 21 }$$

where δ_E are the fluctuations of the TeV photon emissivity.

A.2. Window Function

As the universe is infinite and asymptotically flat, we can expand the fluctuations into planar waves in order to find the length-scale dependence of heating rate fluctuations,

$$\begin{eqnarray}&&{\delta }_{{\rm{H}}}({\boldsymbol{x}})=\displaystyle \frac{1}{{(2\pi )}^{3}}{\displaystyle \int }_{-\infty }^{\infty }{d}^{3}{\boldsymbol{k}}^{\prime} {\tilde{\delta }}_{{\rm{H}}}({\boldsymbol{k}}^{\prime} ){e}^{-i{\boldsymbol{k}}^{\prime} \cdot {\boldsymbol{x}}},{\delta }_{{\rm{E}}}({\boldsymbol{r}}+{\boldsymbol{x}})=\displaystyle \frac{1}{{(2\pi )}^{3}}{\displaystyle \int }_{-\infty }^{\infty }{d}^{3}{\boldsymbol{k}}^{\prime} {\tilde{\delta }}_{{\rm{E}}}({\boldsymbol{k}}^{\prime} ){e}^{-i{\boldsymbol{k}}^{\prime} \cdot ({\boldsymbol{r}}+{\boldsymbol{x}})}.\end{eqnarray} \tag{ 22 }$$

Performing a convolution of the density fluctuations with the kernel ${{Cr}}^{-2}{e}^{-r/{D}_{\mathrm{pp}}}$ (with C as an arbitrary constant), and using the Fourier convolution theorem, Equation (21) becomes

$$\begin{eqnarray}&&{\tilde{\delta }}_{{\rm{H}}}(k)={\tilde{\delta }}_{{\rm{E}}}(k)\tilde{W}(k),\end{eqnarray} \tag{ 23 }$$

where the tilde denotes quantities in Fourier space and $\tilde{W}(k)$ is the Fourier transform of ${{Cr}}^{-2}{e}^{-r/{D}_{\mathrm{pp}}}$. This yields

$$\begin{eqnarray}&&{\tilde{\delta }}_{{\rm{H}}}({\boldsymbol{k}})=\displaystyle \frac{1}{4\pi {D}_{\mathrm{pp}}}{\displaystyle \int }_{{\rm{\Omega }}}d{\rm{\Omega }}{\displaystyle \int }_{0}^{\infty }{{dre}}^{-r/{D}_{\mathrm{pp}}}{\tilde{\delta }}_{{\rm{E}}}({\boldsymbol{k}}){e}^{-i{\boldsymbol{k}}\cdot {\boldsymbol{r}}}.\end{eqnarray} \tag{ 24 }$$

Introducing $\mu =\mathrm{cos}\theta$, where θ is the angle between the wave vector and the line of sight, and assuming statistical isotropy and homogeneity, Equation (24) can be simplified as

$$\begin{eqnarray}&&{\tilde{\delta }}_{{\rm{H}}}(k)=\displaystyle \frac{1}{2{D}_{\mathrm{pp}}}{\displaystyle \int }_{-1}^{1}d\mu {\displaystyle \int }_{0}^{\infty }{{dre}}^{-r/{D}_{\mathrm{pp}}}{\tilde{\delta }}_{{\rm{E}}}(k){e}^{-{ikr}\mu }={\tilde{\delta }}_{{\rm{E}}}(k)\displaystyle \frac{1}{{D}_{\mathrm{pp}}}{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{\mathrm{sin}({kr})}{{kr}}{e}^{-r/{D}_{\mathrm{pp}}}{dr}=\displaystyle \frac{{\tilde{\delta }}_{{\rm{E}}}(k)}{{{kD}}_{\mathrm{pp}}}\mathrm{arctan}({{kD}}_{\mathrm{pp}})=\tilde{\delta }(k)\displaystyle \frac{\mathrm{arctan}({{kD}}_{\mathrm{pp}})}{{{kD}}_{\mathrm{pp}}}.\end{eqnarray} \tag{ 25 }$$

In the last step, we have assumed for simplicity that the fluctuations in the TeV emissivity follow the DM density fluctuations, i.e., δ_E = δ. While we may justify this by the tight relation between the mass of the central black hole and the stellar mass of the bulge of the host galaxy (Häring & Rix 2004), we generalize the mapping between emissivity and DM fluctuations in Appendix C. Figure 8 shows window functions with different mean free paths in a purely Newtonian universe. When absorption is present, the impact of large-scale structure (i.e., low wavenumbers k) remains constant as distant sources are absorbed. For scales equal to or larger than the cutoff, heating fluctuations follow density fluctuations and overdense regions get more heat. At smaller scales, the density structure has less impact on the heating, unless a strong overdensity is present.

Zoom In Zoom Out Reset image size

The heating rate in an expanding universe at a comoving point $\hat{{\boldsymbol{x}}}$ at redshift z is given by

$$\begin{eqnarray}&&\dot{Q}(\hat{{\boldsymbol{x}}})={\displaystyle \int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}{dE}\displaystyle \frac{{J}_{{\rm{E}}}(\hat{{\boldsymbol{x}}})}{{D}_{\mathrm{pp}}(E,z)},\end{eqnarray} \tag{ 26 }$$

where the interval E_min to E_max is the energy range over which pair production is possible. Quantities indicated by a hat are measured in the comoving frame.

J_E is the specific intensity per unit energy (in erg s⁻¹ cm⁻² erg⁻¹). Assuming a flat geometry on spatial hypersurfaces of the Friedmann–Robertson–Walker metric, J_E is given by

$$\begin{eqnarray}&&{J}_{{\rm{E}}}(E,\hat{{\boldsymbol{x}}})=\displaystyle \frac{{(1+z)}^{3}}{4\pi }{\displaystyle \int }_{{\rm{\Omega }}}d{\rm{\Omega }}{\displaystyle \int }_{z}^{\infty }\displaystyle \frac{{dz}^{\prime} }{1+z^{\prime} }\;\displaystyle \frac{d\hat{r}}{{dz}^{\prime} }{\hat{{\mathcal{E}}}}_{{{\rm{E}}}^{\prime }}(\hat{{\boldsymbol{r}}}+\hat{{\boldsymbol{x}}},z^{\prime} ){e}^{-\tau (E,z,z^{\prime} )},\end{eqnarray} \tag{ 27 }$$

where the vector $\hat{{\boldsymbol{r}}}(\hat{r}(z,z^{\prime} ),\theta ,\phi )$, z is the absorption redshift, ${\hat{{\mathcal{E}}}}_{{\rm{E}}}={\epsilon }_{{\rm{E}}}/{(1+z^{\prime} )}^{3}$ is the comoving specific emissivity defined in Equation (7), $E^{\prime} =E(1+z^{\prime} )/(1+z)$ and τ is the approximate optical depth given by

$$\begin{eqnarray}&&\tau (E,z,z^{\prime} )\equiv \tau ={\displaystyle \int }_{z}^{{z}^{\prime }}{dz}^{\prime\prime} \displaystyle \frac{1}{{\hat{D}}_{\mathrm{pp}}(E,z^{\prime\prime} )}\displaystyle \frac{d\hat{r}}{{dz}^{\prime\prime} }.\end{eqnarray} \tag{ 28 }$$

An alert reader will recognize that the true optical depth must include the effect of a changing $E^{\prime\prime} =E(1+z^{\prime\prime} )/(1+z)$ as we integrate along z'' from z to z'. However, here we assume that the typical mean free paths for high-energy photons are much smaller than the Hubble length so that E'' ≈ E over the redshift interval, which is small compared to z. We then find that the energy-integrated and volume-integrated heating rate is given by

$$\begin{eqnarray}&&\dot{Q}[\hat{{\boldsymbol{x}}}(z)]=\displaystyle \frac{{(1+z)}^{3}}{4\pi }{\displaystyle \int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}\displaystyle \frac{{dE}}{{D}_{\mathrm{pp}}(E,z)}{\displaystyle \int }_{{\rm{\Omega }}}d{\rm{\Omega }}{\displaystyle \int }_{z}^{\infty }\displaystyle \frac{{dz}^{\prime} }{1+z^{\prime} }\;\displaystyle \frac{d\hat{r}}{{dz}^{\prime} }{\hat{{\mathcal{E}}}}_{{{\rm{E}}}^{\prime }}(\hat{{\mathcal{R}}},z^{\prime} ){e}^{-\tau (E,z,z^{\prime} )},\end{eqnarray} \tag{ 29 }$$

where $\hat{{\mathcal{R}}}=\hat{{\boldsymbol{x}}}+\hat{{\boldsymbol{r}}}$. Similar to Equation (20), the mean heating rate is given by

$$\begin{eqnarray}&&\bar{\dot{Q}}(z)=\displaystyle \frac{{(1+z)}^{3}}{4\pi }{\displaystyle \int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}\displaystyle \frac{{dE}}{{D}_{\mathrm{pp}}(E,z)}{\displaystyle \int }_{{\rm{\Omega }}}d{\rm{\Omega }}{\displaystyle \int }_{z}^{\infty }\displaystyle \frac{{dz}^{\prime} }{1+z^{\prime} }\;\displaystyle \frac{d\hat{r}}{{dz}^{\prime} }{\bar{\hat{{\mathcal{E}}}}}_{{{\rm{E}}}^{\prime }}(z^{\prime} ){e}^{-\tau (E,z,z^{\prime} )}.\end{eqnarray} \tag{ 30 }$$

This enables us to define the heating rate fluctuations

$$\begin{eqnarray}&&{\delta }_{{\rm{H}}}[\hat{{\boldsymbol{x}}}(z)]=\displaystyle \frac{\dot{Q}[\hat{{\boldsymbol{x}}}(z)]-\bar{\dot{Q}}(z)}{\bar{\dot{Q}}(z)}=\displaystyle \frac{{(1+z)}^{3}c}{4\pi \bar{\dot{Q}}}{\displaystyle \int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}\displaystyle \frac{{dE}}{{D}_{\mathrm{pp}}(E,z)}{\displaystyle \int }_{{\rm{\Omega }}}d{\rm{\Omega }}{\displaystyle \int }_{z}^{\infty }{dz}^{\prime} \\ &&\times \displaystyle \frac{[{\hat{{\mathcal{E}}}}_{{{\rm{E}}}^{\prime }}(\hat{{\mathcal{R}}},z^{\prime} )-{\bar{\hat{{\mathcal{E}}}}}_{{{\rm{E}}}^{\prime }}(z^{\prime} )]{e}^{-\tau (E,z,z^{\prime} )}}{(1+z^{\prime} )H(z^{\prime} )}=\displaystyle \frac{{(1+z)}^{3}c}{4\pi \bar{\dot{Q}}}{\displaystyle \int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}\displaystyle \frac{{dE}}{{D}_{\mathrm{pp}}(E,z)}{\displaystyle \int }_{{\rm{\Omega }}}d{\rm{\Omega }}{\displaystyle \int }_{z}^{\infty }{dz}^{\prime} \displaystyle \frac{{\bar{\hat{{\mathcal{E}}}}}_{{{\rm{E}}}^{\prime }}(z^{\prime} ){\delta }_{{\rm{E}}}(\hat{{\mathcal{R}}},z^{\prime} ){e}^{-\tau (E,z,z^{\prime} )}}{(1+z^{\prime} )H(z^{\prime} )}.\end{eqnarray} \tag{ 31 }$$

The TeV emission is related to the presence of supermassive black holes at the center of galaxies, which are located in collapsed DM halos. We can thus connect the fluctuations of the TeV emission, within a certain radius $\hat{r}$, to the underlying DM fluctuations δ. At this stage we assume that TeV fluctuations exactly match the DM fluctuations δ_E = δ. We explain in the next section that this is not exactly true and will take into account various corrections. The initial density fluctuations represent a Gaussian random field, whose exact properties depend on the earliest stages of the universe prior to recombination (Peebles 1982; Bardeen et al. 1986). They grow linearly between z' and z following (Heath 1977)

$$\begin{eqnarray}&&\delta (\hat{{\mathcal{R}}},z)=\delta (\hat{{\mathcal{R}}},z^{\prime} )D(z)/D(z^{\prime} ),\end{eqnarray} \tag{ 32 }$$

where the linear growth factor is given by

$$\begin{eqnarray}&&D(z)={D}_{0}H(z){\displaystyle \int }_{z}^{\infty }{dz}^{\prime} \displaystyle \frac{1+z^{\prime} }{{H}^{3}(z^{\prime} )}.\end{eqnarray} \tag{ 33 }$$

The linear approximation breaks down when the amplitude of the root mean square of the perturbations approaches unity. The evolution of the density field is then determined by the spherical collapse (Gunn & Gott 1972) and the virialization of halos. In the linear regime, the growth of the modes is independent of the wavenumber and we have

$$\begin{eqnarray}&&{\delta }_{{\rm{E}}}(\hat{{\boldsymbol{r}}},z^{\prime} )=\delta (\hat{{\boldsymbol{r}}},z^{\prime} )=\delta (\hat{{\boldsymbol{r}}},z)\displaystyle \frac{D(z^{\prime} )}{D(z)}=\displaystyle \frac{D(z^{\prime} )}{D(z)}\displaystyle \frac{1}{{(2\pi )}^{3}}\displaystyle \int {d}^{3}\hat{{\boldsymbol{k}}}^{\prime} \tilde{\delta }(\hat{{\boldsymbol{k}}}^{\prime} ){e}^{-i\hat{{\boldsymbol{k}}}^{\prime} \cdot \hat{{\boldsymbol{r}}}}.\end{eqnarray} \tag{ 34 }$$

And Equation (31) rewrites to

$$\begin{eqnarray}&&{\delta }_{{\rm{H}}}[\hat{{\boldsymbol{x}}}(z)]=\displaystyle \frac{{(1+z)}^{3}c}{4\pi \bar{\dot{Q}}}{\displaystyle \int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}\displaystyle \frac{{dE}}{{D}_{\mathrm{pp}}(E,z)}{\displaystyle \int }_{{\rm{\Omega }}}d{\rm{\Omega }}{\displaystyle \int }_{z}^{\infty }{dz}^{\prime} \displaystyle \frac{D(z^{\prime} )}{D(z)}\displaystyle \frac{{\bar{\hat{{\mathcal{E}}}}}_{{{\rm{E}}}^{\prime }}(z^{\prime} )\delta (\hat{{\mathcal{R}}},z){e}^{-\tau (E,z,z^{\prime} )}}{(1+z^{\prime} )H(z^{\prime} )}.\end{eqnarray} \tag{ 35 }$$

Fourier transforming yields $\tilde{\delta }(\hat{{\boldsymbol{k}}})$ on the left side while the right side transforms in a similar fashion to Equation (25). As the power spectrum of density fluctuations is statistically isotropic and homogeneous, this simplifies to

$$\begin{eqnarray}&&{\tilde{\delta }}_{{\rm{H}}}(\hat{k},z)=\tilde{\delta }(\hat{k},z)\displaystyle \frac{{(1+z)}^{3}c}{\bar{\dot{Q}}}{\displaystyle \int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}\displaystyle \frac{{dE}}{{D}_{\mathrm{pp}}(E,z)}{\displaystyle \int }_{z}^{\infty }{dz}^{\prime} \displaystyle \frac{D(z^{\prime} )}{D(z)}\displaystyle \frac{\mathrm{sin}(\hat{k}\hat{r}(z,z^{\prime} ))}{\hat{k}\hat{r}(z,z^{\prime} )}\displaystyle \frac{{\bar{\hat{{\mathcal{E}}}}}_{{{\rm{E}}}^{\prime }}(z^{\prime} ){e}^{-\tau (E,z,z^{\prime} )}}{(1+z^{\prime} )H(z^{\prime} )}.\end{eqnarray} \tag{ 36 }$$

Figure 9 shows the window function in an expanding universe at different redshifts. The position of the breaks results from a combination between the value of the mean free path and the redshift evolution of the blazar luminosity density. As discussed in the main text, the rapidly decreasing blazar luminosity density with increasing redshift for z ≳ 2, which moves the break toward larger wave vectors, dominates over the increasing comoving mean free path, which moves the break toward smaller wave vectors. This is due to the fact that the integration over the energy of the photons, which includes lower-energy photons with larger mean free paths, allows the former effect to dominate and moves the break toward smaller wave vectors with increasing redshift. The TeV emission fluctuations are not exactly equal to the DM density fluctuations. In the next section we will account for the various corrections that have to be taken into account to determine a more accurate window function.

Zoom In Zoom Out Reset image size

The TeV fluctuations are biased with respect to the DM fluctuations, as we detail in Section 2.2, yielding

$$\begin{eqnarray}&&{\delta }_{{\rm{E}}}(\hat{{\boldsymbol{x}}},z)=b(z)\delta (\hat{{\boldsymbol{x}}},z),\end{eqnarray} \tag{ 37 }$$

with b as the Eulerian bias. If galaxies were moving exactly with the Hubble flow, their redshift would yield their exact distance to an observer. However, structures falling toward a central potential of a cluster have an infall velocity toward the central overdensity, yielding different overdensities in redshift space ${\delta }_{{\rm{E}}}^{(\hat{s})}$ and real space. In our case, the blazar luminosity density evolution, the mean free path, and the bias as a function of redshift are all subject to these redshift space distortions, as measured by the z = 0 observer. We follow the derivation of the redshift space distortions of the density fluctuations of Mo et al. (2011). The key difference that we encounter is that redshift space distortions impact on $\hat{{\mathcal{R}}}=\hat{{\boldsymbol{x}}}+\hat{{\boldsymbol{r}}}$, whereas our computation involves the displacements $\hat{{\boldsymbol{r}}}$ between the heated point and a distant source. To this end, we define a comoving redshift space distance

$$\begin{eqnarray}&&\hat{{\boldsymbol{s}}}\equiv \displaystyle \frac{\dot{\hat{{\mathcal{R}}}}}{{aH}(z)}=\hat{{\mathcal{R}}}+\displaystyle \frac{{\hat{v}}_{{\mathcal{R}}}}{H(z)}{{\boldsymbol{e}}}_{{\mathcal{R}}},\end{eqnarray} \tag{ 38 }$$

which is given in units of length, and $a=1/(1+z)$ is the cosmological scale factor, ${\hat{v}}_{{\mathcal{R}}}$ is the radial component of the comoving peculiar velocity ($\hat{{\boldsymbol{v}}}={\boldsymbol{v}}/a$ with ${\boldsymbol{v}}$ denoting the physical peculiar velocity), and ${{\boldsymbol{e}}}_{{\mathcal{R}}}$ is the unit vector along the line of sight. The conservation of TeV blazar number then implies

$$\begin{eqnarray}&&{\hat{n}}_{{\rm{E}}}^{(\hat{s})}(\hat{{\boldsymbol{s}}}){d}^{3}\hat{{\boldsymbol{s}}}={\hat{n}}_{{\rm{E}}}(\hat{{\mathcal{R}}}){d}^{3}\hat{{\mathcal{R}}},\qquad {\rm{or}}\qquad [1+{\delta }_{{\rm{E}}}^{(\hat{s})}(\hat{{\boldsymbol{s}}})]{d}^{3}\hat{{\boldsymbol{s}}}=[1+{\delta }_{{\rm{E}}}(\hat{{\mathcal{R}}})]{d}^{3}\hat{{\mathcal{R}}},\end{eqnarray} \tag{ 39 }$$

where we employed the fact that the mean number density is the same in both spaces. Assuming statistical isotropy, we find

$$\begin{eqnarray}&&1+{\delta }_{{\rm{E}}}^{(\hat{s})}(\hat{{\boldsymbol{s}}})=[1+{\delta }_{{\rm{E}}}(\hat{{\mathcal{R}}})]\displaystyle \frac{{\hat{{\mathcal{R}}}}^{2}}{{(\hat{{\mathcal{R}}}+{\hat{v}}_{{\mathcal{R}}}/H(z))}^{2}}{(1+\displaystyle \frac{1}{H(z)}\displaystyle \frac{\partial {\hat{v}}_{{\mathcal{R}}}}{\partial \hat{{\mathcal{R}}}})}^{-1}.\end{eqnarray} \tag{ 40 }$$

Assuming $| {v}_{{\mathcal{R}}}| \ll \hat{{\mathcal{R}}}H(z)$ and only keeping  terms up to the lowest non-trivial order, we obtain

$$\begin{eqnarray}&&{\delta }_{{\rm{E}}}^{(\hat{s})}[\hat{{\boldsymbol{s}}}(\hat{{\mathcal{R}}})]={\delta }_{{\rm{E}}}(\hat{{\mathcal{R}}})-\displaystyle \frac{1}{H(z)}(\displaystyle \frac{\partial }{\partial \hat{{\mathcal{R}}}}+\displaystyle \frac{2}{\hat{{\mathcal{R}}}}){\hat{v}}_{{\mathcal{R}}}\approx {\delta }_{{\rm{E}}}(\hat{{\mathcal{R}}})-\displaystyle \frac{1}{H(z)}\displaystyle \frac{\partial {\hat{v}}_{{\mathcal{R}}}}{\partial \hat{{\mathcal{R}}}},\end{eqnarray} \tag{ 41 }$$

where we assumed in the second step that the scale of the fluctuations is small compared to the distance to the sources (i.e., we used the plane-parallel approximation according to Kaiser 1987). In the linear regime (${\delta }_{{\rm{E}}}(\hat{{\mathcal{R}}})\ll 1$), the peculiar velocity field is related to the density field via the Zel'dovich approximation (Mo et al. 2011),

$$\begin{eqnarray}&&{\delta }_{{\rm{E}}}(\hat{{\mathcal{R}}})=-\displaystyle \frac{b}{{fH}(z)}\displaystyle \frac{\partial \hat{{\boldsymbol{v}}}}{\partial \hat{{\mathcal{R}}}},\end{eqnarray} \tag{ 42 }$$

where $f\equiv d\mathrm{log}\delta /d\mathrm{log}a$. It can be shown in the linear regime that non-trivial solutions of the velocity field must be irrotational so that ${\boldsymbol{v}}={{\boldsymbol{\nabla }}}_{\hat{{\mathcal{R}}}}{{\rm{\Phi }}}_{v}$. Using this property, we can rewrite Equation (42) to yield

$$\begin{eqnarray}&&{\hat{v}}_{{\mathcal{R}}}=-\displaystyle \frac{{fH}(z)}{b}\displaystyle \frac{\partial }{\partial \hat{{\mathcal{R}}}}{{\boldsymbol{\nabla }}}_{\hat{{\mathcal{R}}}}^{-2}{\delta }_{{\rm{E}}}(\hat{{\mathcal{R}}}),\end{eqnarray} \tag{ 43 }$$

where ${{\boldsymbol{\nabla }}}_{\hat{{\mathcal{R}}}}^{-2}$ represents the inverse Laplacian. It follows that

$$\begin{eqnarray}&&{\delta }_{{\rm{E}}}^{(\hat{s})}[\hat{{\mathcal{R}}}(\hat{{\boldsymbol{r}}})]=[1+\displaystyle \frac{f}{b}\displaystyle \frac{{\partial }^{2}}{\partial {\hat{{\mathcal{R}}}}^{2}}{{\boldsymbol{\nabla }}}_{\hat{{\mathcal{R}}}}^{-2}]{\delta }_{{\rm{E}}}[\hat{{\mathcal{R}}}(\hat{{\boldsymbol{r}}})]=[1+\displaystyle \frac{f}{b}\displaystyle \frac{{\partial }^{2}}{\partial {\hat{r}}^{2}}{{\boldsymbol{\nabla }}}_{\hat{{\boldsymbol{r}}}}^{-2}]{\delta }_{{\rm{E}}}[\hat{{\mathcal{R}}}(\hat{{\boldsymbol{r}}})].\end{eqnarray} \tag{ 44 }$$

This demonstrates that redshift space distortions do not depend on the line of sight distance from the observer to the source but only on the shape of the local gravitational potential of the source! This property enables us to treat redshift space distortions locally by performing a Fourier decomposition of the density field into plane waves, which yields

$$\begin{eqnarray}&&{\tilde{\delta }}_{{\rm{E}}}^{(\hat{s})}[{\hat{{\boldsymbol{k}}}}^{(\hat{s})}(\hat{{\boldsymbol{k}}})]=(1+\displaystyle \frac{f}{b}{\mu }^{2}){\tilde{\delta }}_{{\rm{E}}}(\hat{{\boldsymbol{k}}})=(b+f{\mu }^{2})\tilde{\delta }(\hat{{\boldsymbol{k}}}),\end{eqnarray} \tag{ 45 }$$

where $\mu \equiv {k}_{z}/k$. Keeping only the first order corrections for density fluctuations, Equation (29) yields

$$\begin{eqnarray}&&\dot{Q}(\hat{{\boldsymbol{x}}},z)=\displaystyle \frac{{(1+z)}^{3}c}{4\pi }{\displaystyle \int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}\displaystyle \frac{{dE}}{{D}_{\mathrm{pp}}(E,z)}{\displaystyle \int }_{{\rm{\Omega }}}d{\rm{\Omega }}{\displaystyle \int }_{z}^{\infty }{dz}^{\prime} \displaystyle \frac{{\bar{\hat{{\mathcal{E}}}}}_{{{\rm{E}}}^{\prime }}(z^{\prime} ){e}^{-\tau (E,z,z^{\prime} )}}{(1+z^{\prime} )H(z^{\prime} )}[1+{\delta }_{{\rm{E}}}^{(\hat{s})}(\hat{{\mathcal{R}}},z^{\prime} )].\end{eqnarray} \tag{ 46 }$$

Subtracting the mean heating rate (Equation (30)) then yields the fluctuations

$$\begin{eqnarray}&&{\delta }_{{\rm{H}}}(\hat{{\boldsymbol{x}}},z)=\displaystyle \frac{1}{4\pi \bar{\dot{Q}}}{\displaystyle \int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}{dE}{\displaystyle \int }_{{\rm{\Omega }}}d{\rm{\Omega }}{\displaystyle \int }_{z}^{\infty }{dz}^{\prime} X(E,z,z^{\prime} ){\delta }_{{\rm{E}}}^{(\hat{s})}(\hat{{\mathcal{R}}},z^{\prime} ),\end{eqnarray} \tag{ 47 }$$

where we introduced

$$\begin{eqnarray}&&X(E,z,z^{\prime} )=\displaystyle \frac{{(1+z)}^{3}c}{{D}_{\mathrm{pp}}(E,z)}\displaystyle \frac{{\bar{\hat{{\mathcal{E}}}}}_{{{\rm{E}}}^{\prime }}(z^{\prime} ){e}^{-\tau (E,z,z^{\prime} )}}{(1+z^{\prime} )H(z^{\prime} )}\end{eqnarray} \tag{ 48 }$$

for convenience. Transforming to Fourier space yields

$$\begin{eqnarray}&&{\tilde{\delta }}_{{\rm{H}}}(\hat{{\boldsymbol{k}}},z)=\displaystyle \frac{1}{4\pi \bar{\dot{Q}}}{\displaystyle \int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}{dE}{\displaystyle \int }_{{\rm{\Omega }}}d{\rm{\Omega }}{\displaystyle \int }_{z}^{\infty }{dz}^{\prime} \displaystyle \frac{X(E,z,z^{\prime} )}{{(2\pi )}^{3}}\displaystyle \int {d}^{3}\hat{{\boldsymbol{x}}}{e}^{i\hat{{\boldsymbol{k}}}\cdot \hat{{\boldsymbol{x}}}}\displaystyle \int {d}^{3}\hat{{\boldsymbol{k}}^{\prime} }{e}^{-i\hat{{\boldsymbol{k}}^{\prime} }\cdot (\hat{{\boldsymbol{r}}}+\hat{{\boldsymbol{x}}})}\tilde{\delta }(\hat{{\boldsymbol{k}}^{\prime} },z)\displaystyle \frac{D(z^{\prime} )}{D(z)}[b(z^{\prime} )+f{\mu }^{2}].\end{eqnarray} \tag{ 49 }$$

Assuming statistical isotropy and homogeneity allows us to simplify this expression, and we arrive at

$$\begin{eqnarray}&&{\tilde{\delta }}_{{\rm{H}}}(z,\hat{k})=\displaystyle \frac{\tilde{\delta }(z,\hat{k})}{2\bar{\dot{Q}}}{\displaystyle \int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}{dE}{\displaystyle \int }_{{\rm{\Omega }}}d{\rm{\Omega }}{\displaystyle \int }_{z}^{\infty }{dz}^{\prime} X(E,z,z^{\prime} )\displaystyle \frac{D(z^{\prime} )}{D(z)}{\displaystyle \int }_{-1}^{1}d\mu [{e}^{-i\hat{k}\hat{r}\mu }(b+f{\mu }^{2})]\\ &&=\displaystyle \frac{\tilde{\delta }(z,\hat{k})}{\bar{\dot{Q}}}{\displaystyle \int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}{dE}{\displaystyle \int }_{z}^{\infty }{dz}^{\prime} X(E,z,z^{\prime} )\displaystyle \frac{D(z^{\prime} )}{D(z)}[(b(z^{\prime} )+\displaystyle \frac{f}{3}){j}_{0}(\hat{k}\hat{r})-\displaystyle \frac{2}{3}{{fj}}_{2}(\hat{k}\hat{r})],\end{eqnarray} \tag{ 50 }$$

where we used the spherical Bessel functions

$$\begin{eqnarray}&&{j}_{0}(x)=\displaystyle \frac{\mathrm{sin}(x)}{x},\end{eqnarray} \tag{ 51 }$$

$$\begin{eqnarray}&&{j}_{2}(x)=(\displaystyle \frac{3}{{x}^{2}}-1)\displaystyle \frac{\mathrm{sin}(x)}{x}-\displaystyle \frac{3\mathrm{cos}(x)}{{x}^{2}}.\end{eqnarray} \tag{ 52 }$$

We have further used

$$\begin{eqnarray}&&{\displaystyle \int }_{-1}^{1}{\mu }^{2}{e}^{-{ikr}\mu }d\mu =\displaystyle \frac{2\mathrm{sin}({kr})}{{kr}}+4\displaystyle \frac{\mathrm{cos}({kr})}{{({kr})}^{2}}-4\displaystyle \frac{\mathrm{sin}({kr})}{{({kr})}^{3}}.\end{eqnarray} \tag{ 53 }$$

The window function for the heating rate fluctuations is then given by

$$\begin{eqnarray}&&{\tilde{W}}_{{\rm{H}}}(z,\hat{k})=\displaystyle \frac{{(1+z)}^{3}c}{\bar{\dot{Q}}}{\displaystyle \int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}\displaystyle \frac{{dE}}{{D}_{\mathrm{pp}}(E,z)}{\displaystyle \int }_{z}^{\infty }{dz}^{\prime} \displaystyle \frac{{\bar{\hat{{\mathcal{E}}}}}_{{{\rm{E}}}^{\prime }}(z^{\prime} ){e}^{-\tau (E,z,z^{\prime} )}}{(1+z^{\prime} )H(z^{\prime} )}\displaystyle \frac{D(z^{\prime} )}{D(z)}[(b(z^{\prime} )+\displaystyle \frac{f}{3}){j}_{0}(\hat{k}\hat{r})-\displaystyle \frac{2}{3}{{fj}}_{2}(\hat{k}\hat{r})]\end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray}&&=({\displaystyle \int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}\displaystyle \frac{{dE}}{{D}_{\mathrm{pp}}(E,z)}{\displaystyle \int }_{z}^{\infty }{dz}^{\prime} \displaystyle \frac{{\bar{\hat{{\mathcal{E}}}}}_{{{\rm{E}}}^{\prime }}(z^{\prime} ){e}^{-\tau (E,z,z^{\prime} )}}{(1+z^{\prime} )H(z^{\prime} )}\displaystyle \frac{D(z^{\prime} )}{D(z)}[(b(z^{\prime} )+\displaystyle \frac{f}{3}){j}_{0}(\hat{k}\hat{r})-\displaystyle \frac{2f}{3}{j}_{2}(\hat{k}\hat{r})])/\\ &&({\displaystyle \int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}\displaystyle \frac{{dE}}{{D}_{\mathrm{pp}}(E,z)}{\displaystyle \int }_{z}^{\infty }{dz}^{\prime} {\bar{\hat{{\mathcal{E}}}}}_{{{\rm{E}}}^{\prime }}(z^{\prime} )\displaystyle \frac{{e}^{-\tau (E,z,z^{\prime} )}}{(1+z^{\prime} )H(z^{\prime} )}).\end{eqnarray} \tag{ 55 }$$

This equation is similar to Pritchard & Furlanetto (2007) and Barkana & Loeb (2005) but does not contain the unnecessary area correction factor, which would be present only for a Lagrangian description of bias.