A MEASUREMENT OF SMALL-SCALE STRUCTURE IN THE 2.2 ⩽ z ⩽ 4.2 Lyα FOREST

After a low-density gas element is photoheated during reionization, it will subsequently cool and gas elements with similar photoheating histories generally land on a "temperature–density relation" (Hui & Gnedin 1997):

$$\begin{eqnarray} &&T({\rm r}) = T_0({\rm r})\, [1 + \delta _\rho ({\rm r})]^{\gamma ({\rm r})-1}. \end{eqnarray} \tag{ 1 }$$

Here δ_ρ(r) denotes the fractional gas over-density (implicitly smoothed on the Jeans scale) at spatial position r. T₀ is the temperature of a gas element at the cosmic mean density, and the power-law index γ approximates the density dependence of the temperature field. The temperature that a gas element reaches at say, z = 3, depends on the temperature that it reaches during reionization and on its subsequent cooling and heating. The temperature attained by each gas element during reionization depends mostly on the shape of the spectrum of the sources that ionize it. The relevant spectrum is generally modified from the intrinsic spectral shape of an ionizing source, owing to intervening material between a source and the gas element in question, which tends to harden the ionizing spectrum. After a gas element is photoheated during reionization, adiabatic cooling owing to the expansion of the universe is the dominant cooling mechanism (for the bulk of the low-density gas that makes up the Lyα forest).⁷ When a gas element is significantly ionized during reionization it reaches photoionization equilibrium and receives only a small amount of additional photoheating as low levels of residual neutral material are ionized. During reionization, gas elements gain heat as hydrogen is ionized, as helium is singly ionized, and when helium is doubly ionized. If helium is doubly ionized significantly after hydrogen is ionized, two separate "reionization events" may be important in determining the thermal history of the IGM. As both hydrogen and helium reionization are extended inhomogeneous processes, T₀ and γ may be strong functions of spatial position following reionization events. However, once a sufficiently long time passes after reionization, gas elements reach a "thermal asymptote" and lose memory of the initial photoheating during reionization (Hui & Gnedin 1997). At this point the inhomogeneities in T₀ and γ should again be small.

In the absence of He ii photoheating, one expects the temperature of the IGM at z ≈ 3 to be T₀ ≲ 10, 000 K, with the precise temperature depending on the timing of hydrogen reionization and the nature of the ionizing sources (Hui & Haiman 2003). Sufficiently long after a reionization event, the slope of the temperature–density relation, γ, tends to γ ≈ 1.6, owing to the competition between adiabatic cooling and residual photoionization heating (Hui & Gnedin 1997). He ii reionization likely raises the temperature of the IGM by roughly 10, 000 K, with the precise increase depending on the spectrum of the ionizing sources and other factors. He ii photoheating and the spread in timing of He ii reionization flatten the temperature–density relation to γ ≈ 1.3 (McQuinn et al. 2008).

The temperature of the IGM has three separate effects on Lyα forest spectra. First, increasing the temperature of the absorbing gas increases the amount of Doppler broadening: thermal motions spread the absorption of a gas element out over a length (in velocity units) of $b = \sqrt{2 k T/m_p} \approx 13$ km s⁻¹ for T = 10⁴ K gas. Second, the gas pressure and Jeans smoothing scale increase with increasing temperature. Since it takes some time for the gas to move around and the gas pressure to adjust to prior heating, this effect is sensitive not to the instantaneous temperature, but to prior heating (Gnedin & Hui 1998). This effect is more challenging for simulators to capture because properly accounting for it requires re-running entire simulations after adjusting the simulated ionization/reheating history. The Jeans smoothing effect is not completely degenerate, however, with the Doppler broadening one because Jeans-smoothing smooths the gas distribution in three dimensions, while Doppler broadening smooths the optical depth in one dimension (Zaldarriaga et al. 2001).⁸ Finally, the recombination coefficient is temperature dependent, scaling as T^−0.7: hotter gas recombines more slowly and reaches a lower neutral fraction than cooler gas.

The first two of these effects mostly impact the amplitude of small-scale fluctuations in the Lyα forest (e.g., Zaldarriaga et al. 2001). For the range of models we are interested in presently, the first effect (Doppler broadening) should be the dominant influence on the small-scale power. At a given redshift, the small-scale structure in the Lyα forest is most sensitive to the temperature of absorbing gas at some characteristic density, with less dense gas giving very little absorption and more dense gas giving rise to mostly saturated absorption. At z ≈ 3 the forest is sensitive mostly to the temperature of gas a little more dense than the cosmic mean (McDonald et al. 2001; Zaldarriaga et al. 2001). At higher redshifts, the absorption is sensitive to the temperature of somewhat less dense gas, while at lower redshifts the absorption depends on more dense gas (Davé et al. 1999).

Next we describe our method for constraining T(r) (Equation (1)) from absorption spectra. Following earlier work (Theuns & Zaroubi 2000; Zaldarriaga 2002; Theuns et al. 2002b), we convolve Lyα transmission spectra with a filter that pulls out high-k modes across each spectrum. As mentioned above, Doppler broadening convolves the optical depth field with a Gaussian filter with a—temperature-dependent—width of tens of km s⁻¹. We hence desire a filter that extracts Fourier modes with wavelengths of tens of km s⁻¹ across each spectrum.

We have found that a very simple choice of filter accomplishes this task. In configuration space, the filter we use may be written as

$$\begin{eqnarray} &&\Psi _n(x) = A {\rm exp}(i k_0 x) {\rm exp}\left[-\frac{x^2}{2 s_n^2}\right]. \end{eqnarray} \tag{ 2 }$$

We fix the normalization, A, by requiring the filter to have unit power—i.e., after filtering a white-noise field with noise power spectrum P_N(k) = Δuσ², the filtered field has variance σ². (Δu denotes the size of a spectral pixel in velocity units.)⁹ With this normalization, the filter's Fourier transform in k-space is

$$\begin{eqnarray} &&\Psi _n(k) = \pi ^{-1/4} \sqrt{\frac{2 \pi s_n}{\Delta u}} {\rm exp}\left[-\frac{(k - k_0)^2 s_n^2}{2}\right]. \end{eqnarray} \tag{ 3 }$$

In configuration space this filter is simply a plane wave, damped by a Gaussian. In Fourier space, the filter is a Gaussian centered around k = k₀. We would like the filter to have zero mean. Throughout this work we choose k₀s_n = 6, in which case Equation (3) shows that the zero mode of the filter Ψ_n(k = 0) is extremely close to zero, satisfying closely the zero mean requirement. This filter clearly has the properties of being localized in both configuration space and Fourier space. These are among the defining properties of a "wavelet filter," and the filter of Equations (2) and (3) is known as a "Morlet Wavelet" in the wavelet literature.¹⁰ We plot its form in Figure 1 for s_n = 34.9 km s⁻¹, which, as we discuss further below, turns out to be one convenient choice. Note that the filters Ψ_n (Equation (2)) do not form an orthogonal set, but this is unnecessary for our present purposes. We do not expand the entire spectrum in terms of a wavelet basis in this work—the Morlet wavelet, with locality in real and configuration space, is simply a convenient filter.

Zoom In Zoom Out Reset image size

We then convolve each observed (or simulated) spectrum with the above filter. In this paper, we consider throughout the fractional Lyα transmission field, δ_F = (F − 〈F〉)/〈F〉. Here F = e^−τ is the Lyα transmission and 〈F〉 is the global average Lyα transmission. We label the flux field, δ_F, convolved with the filter Ψ_n as a_n:

$$\begin{eqnarray} &&a_n(x) = \int dx^\prime \Psi _n(x - x^\prime) \delta _F(x^\prime), \end{eqnarray} \tag{ 4 }$$

and compute the convolution using fast Fourier transforms. Note that a_n(x) is a complex number for our choice of filter, Ψ_n(x). A measure of small-scale power is then

$$\begin{eqnarray} &&A(x) = |a_n(x)|^2, \end{eqnarray} \tag{ 5 }$$

which for brevity of notation we sometimes refer to as "the wavelet-filtered field" or as "the wavelet amplitudes" (even though it is proportional to the transmission field squared). It is also useful to note that the average wavelet amplitude is just

$$\begin{eqnarray} &&\langle |a_n(x)|^2 \rangle = \int _{-\infty }^{\infty } \frac{dk^\prime }{2 \pi } [\Psi _n(k^\prime)]^2 P_F(k^\prime), \end{eqnarray} \tag{ 6 }$$

with P_F(k) denoting the power spectrum of δ_F. Hence, the mean wavelet amplitude is nothing more than the usual flux power spectrum for some convenient "band" of wavenumbers (see Figure 5 for further illustration). Additional statistics of A(x), beyond the mean, characterize the spatial variations in the small-scale transmission power.

We frequently find it convenient to smooth A(x) using a top-hat filter of smoothing length L:

$$\begin{eqnarray} &&A_L(x) = \frac{1}{L} \int _{-\infty }^{\infty } dx^\prime \Theta (|x - x^\prime |; L/2) A(x^\prime). \end{eqnarray} \tag{ 7 }$$

Here Θ(|x − x'|; L/2) = 1 for |x − x'| ⩽ L/2 and is zero otherwise. Smoothing the wavelet filtered field is desirable since the small-scale power is not a perfect indicator of the local temperature, and smoothing reduces the noisy excursions that the wavelet amplitudes can take. Since the hot regions are expected to be rather large during He ii reionization (McQuinn et al. 2008), we can smooth considerably without diluting any temperature inhomogeneities. We generally adopt L = 1000 km s⁻¹, corresponding to roughly ≈10 co-moving Mpc/h at z = 3. We discuss this choice further below.

Since thermal broadening smooths the optical depth field on tens of km s⁻¹ scales, A_L(x) should be a good tracer of the temperature for suitable choices of s_n. In order to illustrate this concretely, we apply the filter to a simulated spectrum from a simple hypothetical inhomogeneous temperature model, following a similar example from Theuns & Zaroubi (2000). Specifically, we splice together simulated lines of sight (generated from a cosmological simulation with the temperature adjusted in a post-processing step, see Section 4.1) with alternating portions of spectrum drawn from each of a "hot" temperature model with T₀ = 2 × 10⁴ K, and γ = 1.3, and a "cold" temperature model with T₀ = 1 × 10⁴ K, and γ = 1.3. We refer the reader to Section 4.1 and Appendix C for details regarding the simulated spectra. If the wavelet filtered field provides a good indicator of the temperature, regions with hot temperatures should tend to produce low wavelet amplitudes, while the cold regions should produce high wavelet amplitudes. The results of this test are shown in Figure 2, for smoothing scales of s_n = 34.9 km s⁻¹ and L = 1000 km s⁻¹. Cold regions tend to contain several narrow lines and produce a large response after filtering: the regions near Δv = 6000 km s⁻¹ and 15, 000 km s⁻¹ have A_L ≳ 0.02. The hot regions typically have A_L ≲ 0.005 and never reach the large amplitudes found in the cold regions. There is some variance in the wavelet amplitude from region to region—for example, A_L is not as large in the cold region near Δv = 10, 000 km s⁻¹ as it is at Δv = 6000 km s⁻¹ and 15, 000 km s⁻¹. Nonetheless, the smoothed wavelet amplitude is a fairly good tracer of the underlying temperature field.

Zoom In Zoom Out Reset image size

In order to quantify this further, we calculate the probability distribution function (PDF) of the smoothed wavelet amplitudes. We do this for the two choices of small-scale smoothing adopted in this paper (see Section 2.3): s_n = 34.9 km s⁻¹, and twice this, s_n = 69.7 km s⁻¹. The PDF of smoothed wavelet amplitudes will be the main statistic we consider in this paper. For now, we examine models with homogeneous temperature–density relations. The models we select for the temperature–density relation loosely correspond respectively to what one expects right after He ii reionization (T₀ ≈ 20–25, 000 K and γ = 1.3; McQuinn et al. 2008), and to what one expects if H i, He i, and He ii are all ionized much before z ≈ 3 (T₀ ≈ 7500–10, 000 K and γ = 1.6; Hui & Haiman 2003). The latter, cooler model might be expected if, for example, the IGM is reionized by abundant faint quasars which have sufficiently hard spectra to doubly ionize helium at the same time they reionize hydrogen or if high redshift galaxies have a surprisingly hard spectrum and can doubly ionize helium themselves. Note that the precise z ≈ 3 temperature in the early reionization models is determined by residual photoheating and depends on the reprocessed spectra of the post-reionization ionizing sources (Hui & Haiman 2003).

The PDFs in these models are shown for two choices of small-scale smoothing in Figure 3 (s_n = 34.9 km s⁻¹) and Figure 4 (s_n = 69.7 km s⁻¹). A larger range of models will be examined in Section 4. Considering first the smaller smoothing scale (Figure 3), one sees that the peak of the PDF in the T₀ = 20, 000 K, γ = 1.3 model is reached at a smoothed wavelet amplitude that is roughly a factor of 2 smaller than the peak location in the T₀ = 10, 000 K, γ = 1.6 model. The PDFs in the hotter T₀ ≈ 25, 000 K model and the colder T₀ = 7500 K model differ by even more. In the midst of He ii reionization, one expects an inhomogeneous temperature field and the true temperature–density relation may be a mix of the models shown here. At any rate, the wavelet PDFs differ significantly in the models with 20, 000 K and those with cooler temperatures. This further demonstrates—beyond the visual inspection of Figure 2—that the wavelet PDF is a useful statistic for constraining the thermal history and He ii reionization. The typical wavelet amplitude in each model is significantly larger at s_n = 69.7 km s⁻¹ (Figure 4), a consequence of the roughly exponential fall-off in flux power toward high k (Zaldarriaga et al. 2001). The PDFs still vary significantly with temperature–density relation at this larger smoothing scale, although the sensitivity is a little bit reduced.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Before we move on to analyze observational data, let us further consider the two smoothing scales, s_n and L, in our calculations. We make measurements for two choices of small-scale smoothing: s_n = 34.9 km s⁻¹ and s_n = 69.7 km s⁻¹.¹¹ For the former choice of smoothing scale |Ψ_n(k)|² is proportional to a Gaussian centered on k₀ = 6/s_n = 0.17 s km⁻¹, with width $\sigma _k = \sqrt{2}/s_n = 0.04$ s km⁻¹. The latter choice of smoothing scale centers the Gaussian on k₀ = 6/s_n = 0.086 s km⁻¹, with a width of $\sigma _k = \sqrt{2}/s_n = 0.02$ s km⁻¹. The range of scales probed by these filters is shown in comparison to simulated flux power spectra in Figure 5. As illustrated in Figures 3, 4, and 5, the wavelet PDFs are slightly less sensitive to the IGM temperature for the larger smoothing scale filter. On the other hand, the results at the larger smoothing scale are less sensitive to metal line contamination and other systematics. Increasing the smoothing by still another factor of 2 would almost completely remove the sensitivity to temperature (see Figure 5). Decreasing s_n by an additional factor of 2 (to s_n = 17.4 km s⁻¹) increases the fractional difference between model curves, but brings one very far out on the exponential tail of the power spectrum (Figure 5) and makes the results very sensitive to metal line contamination, detector noise, and pixelization effects. The two choices of filtering scale used here represent a compromise between discriminating power and systematic effects. Considering both choices of filtering scale gives a consistency check on the results and helps to protect against systematic effects.

Zoom In Zoom Out Reset image size

Let us now consider the large-scale smoothing, L. Naively, one would want to tune this filtering to precisely the scale on which the temperature field is inhomogeneous. Since the power spectrum of temperature fluctuations during He ii reionization has a relatively well-defined peak (McQuinn et al. 2008), one might expect the variance of the wavelet amplitudes to also show a clear maximum at some characteristic smoothing scale. However, in practice we find that this is washed out in Lyα forest spectra, which as one-dimensional skewers suffer owing to aliasing from high-k modes transverse to the line of sight (Kaiser & Peacock 1991). To illustrate this, consider the two-point function of the wavelet amplitudes (squared),

$$\begin{eqnarray} &&\xi _A(|v_1 - v_2|) = \frac{\langle A(v_1) A(v_2) \rangle - \langle A \rangle ^2}{\langle A \rangle ^2}, \end{eqnarray} \tag{ 8 }$$

and its Fourier transform, the power spectrum of wavelet-amplitudes squared, P_A(k). Here v₁ and v₂ are two points along a quasar spectrum and 〈A〉 is the globally averaged wavelet amplitude squared, and we have normalized this two-point function by the (square of the) mean wavelet amplitude squared. The power spectrum of wavelet amplitude squared fluctuations encodes how much the small-scale power spectrum fluctuates across a quasar spectrum as a function of scale. It involves a product of four values of δ_F and is hence a four-point function.

We show two simulated examples of P_A(k) in Figure 6 for s_n = 34.9 km s⁻¹. One can see that, except for the small-scale cut-off, the power spectra are quite flat as a function of scale. The flatness is a direct consequence of aliasing from high-k modes transverse to the line of sight, whose contribution to P_A(k) on large scales may swamp that from large-scale temperature inhomogeneities. Indeed, we have experimented with various inhomogeneous temperature models, including simulated models from McQuinn et al. (2008) and find similarly flat power spectra. Unfortunately P_A(k) does not directly indicate the scale dependence of temperature fluctuations, as one might naively hope. One might be able to get around this by using quasar pairs to measure the power spectrum of wavelet amplitude squared transverse to the line of sight. We defer, however, investigating this to future work. For the moment, our main conclusion is that, owing to the flatness of P_A(k), the precise smoothing scale L is relatively unimportant. Hence we generally stick to L = 1, 000 km s⁻¹ as a convenient choice. We nevertheless investigate the dependence on large-scale smoothing from observational and simulated data in Section 4.4.

Zoom In Zoom Out Reset image size

To summarize, by applying a very simple filter to a quasar spectrum, we can measure the small-scale power spectrum of transmission fluctuations as a function of position across each spectrum, and thereby constrain the temperature of the IGM. Note that our procedure does not involve identifying absorption lines and fitting profiles to identified lines (although we find in Section 3 that it is important to identify metal absorbers in the forest which does involve line fitting). It is instead within the spirit of treating the forest as a one-dimensional random field and measuring the statistics of this continuous field (e.g., Croft et al. 1998). This is more appropriate given the modern understanding that the forest arises from fluctuations in the line-of-sight density field, rather than discrete absorbing clouds (e.g., Hernquist et al. 1996; Miralda-Escudé et al. 1996; Katz et al. 1996b). In this way our approach is very similar to Theuns & Zaroubi (2000) and Zaldarriaga (2002), and somewhat resembles Zaldarriaga et al. (2001), but is rather different than Schaye et al. (2000), Ricotti et al. (2000), and McDonald et al. (2001).

Additionally, recall that the widths of most of the absorption lines in the Lyα forest are dominated by the Hubble expansion across an absorber, and not by thermal broadening (Hernquist et al. 1996; Weinberg et al. 1998). In order to determine the temperature with a line-fitting method, one typically looks for a low-end cut-off in the distribution of line widths (e.g., Schaye et al. 1999, 2000). One might worry that this throws out information as thermal broadening smooths the spectrum everywhere. In practice, though, it appears that most of the signal and information in our method also arise from deep narrow lines which produce a large response after wavelet filtering. Another possible issue is that the precise interpretation of the line width cut-off in the line-fitting studies is unclear when the temperature field is inhomogeneous. It would certainly be interesting to compare more closely the different methods, but we defer this to future work. For now, note that our method is very simple to apply.

We aim to estimate the small-scale power in a way that minimizes sensitivity to uncertainties in the quasar continuum. Dall'Aglio et al. (2008) carefully continuum fit the data we use here and used Monte Carlo simulations to check the accuracy of their fits. We can further mitigate uncertainties by considering fluctuations in the transmission around the mean, relative to the mean. This is helpful because the overall normalization of the continuum divides out. Provided that the continuum varies slowly across each spectrum in comparison with the fluctuations in the forest, we can additionally remove any slowly varying trend produced by the quasar continuum—or any slowly varying residuals in the case of data that has previously been continuum fitted—and obtain an unbiased estimate of the small-scale structure in the forest (Hui et al. 2001). For each spectrum, we estimate a running mean flux by filtering the data on large scales as in Croft et al. (2002), Kim et al. (2004), and Lidz et al. (2006). Our estimate of the fractional transmission is then

$$\begin{eqnarray} &&\hat{\delta }_F(\Delta v) = \frac{F(\Delta v) - F_R(\Delta v)}{F_R(\Delta v)}. \end{eqnarray} \tag{ 9 }$$

Here F(Δv) is the flux at velocity separation Δv, and F_R(Δv) is the spectrum smoothed with a large radius filter. We use here a Gaussian filter with radius R = 2500 km s⁻¹. One may form $\hat{\delta }_F$ using either the raw flux or a continuum-normalized flux. In this work, we use the continuum fitted data from Dall'Aglio et al. (2008) throughout. The large-scale filter removes any slowly varying trend owing to structure in the underlying quasar continuum from, e.g., weak emission lines, or slowly varying residuals in the case of continuum fitted data. It also means that we sacrifice measuring large-scale modes in the Lyα forest, but we presently focus on small-scale structure, and sufficiently large-scale modes are regardless dominated by structure in the quasar continuum. We refer the reader to Croft et al. (2002) and Lidz et al. (2006) for some tests illustrating the robustness of $\hat{\delta }_F$ to continuum-fitting uncertainties. As a double-check that the present results are insensitive to the precise $\hat{\delta }_F$ estimator, we also generated $\hat{\delta }_F$ with a different choice of large-scale smoothing for one of our redshift bins, R = 10, 000 km s⁻¹—i.e., close to the flat mean case—and found a nearly identical wavelet PDF.

We begin by estimating $\hat{\delta }_F$ across each spectrum, first re-binning, using linear interpolation, all of the data onto uniform pixels in velocity space with Δu = 4.4 km s⁻¹. We consistently use the same binning in constructing simulated spectra. This avoids effects from variable pixelization, while still preserving the scales of interest.¹² After forming $\hat{\delta }_F$ across each spectrum, we break the data into several (contiguous and non-overlapping) redshift bins of full-width Δz = 0.4, centered around $\bar{z} = 2.2, 2.6, 3.0, 3.4, 3.8$, and 4.2. Owing to uneven redshift sampling in the data set, the redshift bin at $\bar{z} = 3.8$ (Dall'Aglio et al. 2008) would be almost entirely empty and so we do not consider it further here. This occurs because most of the spectra in the Dall'Aglio et al. (2008) sample have emission redshift z_em ≲ 3.7, but the sample has two high quality spectra at emission redshift above z_em ≳ 4.6, which contribute extended (≳150 co-moving Mpc) stretches to our highest redshift bin at $\bar{z}=4.2$. We select only spectral regions that lie between rest frame wavelengths of λ_r = 1050 Å and λ_r = 1190 Å. This conservative cut serves to remove spectral regions that may be contaminated by either the proximity effect, by the Lyβ forest (and other higher Lyman series lines), or by Lyβ and OVI emission features. We then form the wavelet amplitude squared field, smoothed at L = 1000 km s⁻¹, using Equations (2)–(7). The resulting spectra and wavelet amplitudes are visually inspected. Regions impacted by DLAs, or with obvious spurious stretches, are removed from the data sample by hand.

It is instructive to examine a few example spectra visually before measuring their detailed statistical properties. In Figures 7–10, we show several spectra, along with the corresponding (smoothed) wavelet amplitudes squared for a few redshift bins. The most conspicuous change across the different redshift bins is the increasing average absorption with increasing redshift. Since we are considering fractional fluctuations, this manifests itself as an increase in the fraction of pixels with $\hat{\delta }_F$ close to −1, with occasional excursions to very large $\hat{\delta }_F$. The next impression provided by the spectra appears at first tantalizing: most regions have low A_L, but there are occasional upward excursions over portions of the spectrum. This behavior is especially apparent for the smaller of the two filtering scales, and is less apparent in the highest redshift case (Figure 10).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Consider for example the spectrum Q2139-44, in the $\bar{z}=3.0$ bin, convolved with a s_n = 34.9 km s⁻¹ Morlet filter, as shown in Figure 9. In this spectrum the regions near Δv = 5000, 7500, and 12, 500 km s⁻¹ all have relatively high wavelet amplitudes, A_L ≳ 0.02, while the rest of the spectral regions have low amplitude. Inspecting the simulated PDF of Figure 3, the low amplitude floor with A_L ≈ 0.005 seems to indicate hot T₀ ≈ 20, 000 K gas, while the regions with A_L ≳ 0.02 seem to require cooler gas T₀ ≲ 7500 K gas.

At first glance, these upward wavelet amplitude excursions seem to be cold regions embedded in an otherwise hot IGM. This is what one naively expects in the midst of He ii reionization: cool regions where H i and He i reionized long ago, and hotter regions where helium is doubly ionized. Before we dispel this fantasy—these regions are contaminated by metal absorbers (see Section 3.2)—let us add some sightlines from simulated models to further illustrate this (Figure 11).¹³ The sightlines show that the low wavelet amplitude floor in the observed spectrum roughly matches the hot IGM sightline. This implies that there are indeed significant quantities of hot ≈20, 000 K gas in the IGM at z = 3. However, the hot model fails to produce the high wavelet amplitude excursions seen in the data. Matching these seems, at first glance, to require a cooler model—one with roughly T₀ ≈ 7500 K, γ = 1.6, for example. (To be clear, note that the simulation and observational data are drawn from different realizations, so one does not expect the simulated case to match the observations region-by-region or feature-by-feature. The meaningful comparison is the overall number of regions with high or low wavelet amplitude.) At first it is tempting to conclude that we are detecting temperature inhomogeneities from incomplete He ii reionization.

Zoom In Zoom Out Reset image size

We need, however, to consider a very important systematic. A hot region that lands at the same wavelengths as a "clump" of prominent narrow metal lines may look to us like a cold region. The wavelet filter just tells us the total level of small-scale power from place to place and does not distinguish whether absorption arises from H i or some other element. To make a robust estimate of the IGM temperature, we need to identify metal line absorbers within the Lyα forest.¹⁴ We expect metal line contamination to be most severe in the low redshift bins, where the fractional contribution of metals to the overall opacity in the forest is highest (e.g., Faucher-Giguère et al. 2008b), and on the smaller of our two filtering scales (see Appendix B).

Naturally, distinguishing metal absorption lines and Lyα lines within the Lyα forest is a challenging and imperfect process. We do, however, have a few separate handles on distinguishing metal lines from Lyα lines within the forest. First, we identify all of the metal absorbers redward of Lyα and look for "partner" transitions. The partner transitions are additional transitions that lie at the same redshift as an identified red-side line, yet which land within the Lyα forest. Next, we search for doublets within the Lyα forest, which can be identified by their distinctive optical depth ratios and by the characteristic separation between a doublet's two components. For instance, C iv is a doublet with a strong component at λ_r = 1548.2 Å, and a weaker component at λ_r = 1550.8 Å, and the ratio of the absorption cross sections of the two components is 2. So C iv should stand out as a doublet with the two components separated by ≈640 km s⁻¹, with the lower wavelength line a factor of 2 stronger than its partner component. Mg ii is another prominent doublet. After identifying a doublet, one can use the estimated redshift of an identified doublet to search for additional transitions at the same redshift: we look for C ii/iii/iv, N ii/iii/v, O i/vi, Mg i/ii, Al ii/iii, Si ii/iii/iv, S vi, and Fe ii, and consider further transitions for DLAs. This approach already identifies a host of metal lines within the forest, but there are inevitably some remaining metal lines left within the forest. For example, there are sometimes absorbers where the doublet features are undetectable owing to line blending. To further mitigate metal line contamination, our final step is to mark extremely narrow lines (with b-parameter b ≲ 7 km s⁻¹) as metals. This final cut amounts to only 25% of the identified lines. Clearly, one needs to be careful about making cuts based on line width: doing so could bias us against detecting cold regions. However, for an H i line to have a linewidth of b ≲ 7 km s⁻¹ it needs to have an implausibly low temperature of T ≲ 3000 K. Hence, we are confident that this final cut does not bias our results, yet it helps protect against remaining unidentified metal lines within the forest. We will subsequently present tests to check how much the results depend on the precise way in which we excise metal line contaminated regions.

One further source of potential systematic bias is that it is difficult to identify metal lines at wavelengths where Lyα produces saturated absorption. This means that regions with strong Lyα absorption are preferentially left in the data sample that is analyzed. Owing to this, the data sample analyzed may not be entirely representative, which could bias the results. This systematic may be strongest at the highest redshifts in our sample, where the Lyα absorption is most saturated. On the other hand, the abundance of metals likely drops off toward high redshift which may make this systematic less important at high redshift.¹⁵

In this paper, we have identified metal lines for 11 of the 40 spectra in our data sample. The identified metals come entirely from portions of spectrum absorbing at z ≲ 3—where we expect the metal line contamination to be strongest—and not in the higher redshift bins. That is, we do not presently have estimates of metal line contamination in the redshift bins centered around $\bar{z}=3.4$ and $\bar{z}=4.2$. In these redshift bins, we will focus entirely on the larger (s_n = 69.7 km s⁻¹) filtering scale where the metal line contamination is less of an issue (Appendix B).

In order to check the influence of metal line contamination, we calculate the wavelet amplitudes as before and excise regions impacted by metal line contamination. An important assumption here is that gas absorbing in a metal line transition at a given wavelength is spatially uncorrelated with gas absorbing in Lyα at nearby wavelengths. If this assumption were violated, we could bias ourselves by preferentially removing regions of above average hydrogen absorption when excising metal contaminated regions. Fortunately, most of the metal line transitions have rest-frame wavelengths that are very different than that of Lyα and so the gas absorbing in a metal transition at a given wavelength is very widely separated (in physical space) from that absorbing in Lyα. Hence the metal and Lyα absorption are uncorrelated. This justifies our approach.¹⁶ Since the wavelet filter is not completely local, pixels with metal line absorption will contaminate neighboring pixels after filtering. Furthermore, we generally smooth the wavelet squared field over L = 1000 km s⁻¹. As a simple and conservative cut, we examine the fraction of contaminated pixels within a smoothing length L around each pixel, and discard a pixel if less than f_m = 95% of its neighbors (within a smoothing length) are metal free.

We find that metal line contamination can have a significant impact, especially for s_n = 34.9 km s⁻¹ and at z ≲ 3. We show a few example sightlines in Figures 12–14. It is striking that the most prominent peaks in the wavelet filtered field at s_n = 34.9 km s⁻¹, shown in the figures, correspond very closely to metal lines. Essentially, our filter was designed to look for temperature inhomogeneities, but it appears most effective at identifying metal-line contaminated regions. In fact, wavelet filtering may be a good way to quickly identify prominent metals in the forest. The metal line contamination is less severe for spectra passed through the larger wavelet filter (s_n = 69.7 km s⁻¹). The amplitude of fluctuations in the forest is much greater on this smoothing scale. The metals also generally contribute more power on the larger smoothing scale, but the amplitude of fluctuations from H i increases more strongly with smoothing scale, and so the metals are fractionally less important on larger scales. This is perhaps seen most easily in the example of Figure 13. In the s_n = 34.9 km s⁻¹ panel of this figure, all of the prominent peaks are metal contaminated regions. In the larger smoothing scale panel, there are some peaks from H i and some from metals, and the heights of the various peaks are comparable. The more significant contamination of the metals on the smaller smoothing scale likely results because the metal lines tend to be narrower than the H i lines. In Appendix B, we find qualitatively similar results by adding metal line absorbers, with empirically derived properties, to mock Lyα forest spectra.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since we can attribute many of the peaks observed in the wavelet amplitudes to metal lines, this does imply, however, that the temperature inhomogeneities cannot be too large. If temperature inhomogeneities were present and large, we would expect to see more high wavelet amplitude regions left over after excising the metals. In particular, consider Figure 11. In this example, we found that the low wavelet amplitude regions of the spectrum are consistent with hot 20, 000 K gas. While we have not identified metal lines for this particular spectrum, our results from other lines of sight clearly suggest that the high wavelet amplitude regions are metal-contaminated rather than genuine cold regions with T₀ ≈ 7500–10, 000 K. The lack of high wavelet amplitude regions after metal excision implies there are few such cold regions left, and that most of the volume of the IGM at z ≈ 3 is hot with T₀ ≈ 20, 000 K (although see Section 4.1 for a discussion regarding the dependence of our results on γ).

It is clear, however, that metal line contamination is a very important systematic for these measurements, although the contamination is less of an issue on the larger smoothing scale and for the high redshift measurements. This issue is not unique to our method, although the detailed impact of metal lines will depend on the precise algorithm for constraining the IGM temperature. For instance, measurements based on fitting the minimum width of absorption lines in the Lyα forest need to carefully avoid including metal lines in the sample of lines used to estimate the temperature. Power spectrum based temperature estimates need to account for the small-scale power contributed by metal absorbers or mask the metal absorbers before estimating the power spectrum.

Let us now move past mere visual inspection and measure statistical properties from the observed spectra. We focus mostly on the PDF of A_L for our fiducial choices of s_n = 34.9 km s⁻¹, s_n = 69.7 km s⁻¹, L = 1000 km s⁻¹, and f_m = 0.95. In each redshift bin, we find the minimum and maximum A_L and then choose ten evenly spaced logarithmic bins in A_L for the PDF measurement. We tabulate the average A_L and the differential PDF in each A_L bin for each redshift bin. The average redshift of pixels in a redshift bin is typically close to the redshift at bin center, and the error bars are still fairly large, so we ignore any issues associated with redshift evolution across each bin and quote all results at the bin center.

We use a jackknife resampling technique to calculate error bars for the PDF measurements. We first estimate the PDF from the entire data sample within a given redshift bin, $\hat{P}(A_i)$. Here $\hat{P}(A_i)$ is the PDF estimate for the ith A_L bin, and A_i is the average wavelet amplitude squared and smoothed within the bin. Next we divide the data set into n_g = 10 subgroups and estimate the PDF of the data sample omitting each subgroup. Let $\tilde{P}_k(A_i)$ represent the PDF estimate omitting the pixels in the kth subgroup. Then our estimate of the jackknife covariance between bins i and j, C(i, j), is

$$\begin{eqnarray} &&C(i,j) = \sum _{k=1}^{n_g} [\hat{P}(A_i) - \tilde{P}_k(A_i)][\hat{P}(A_j) - \tilde{P}_k(A_j)]. \end{eqnarray} \tag{ 10 }$$

In practice our estimates of the off-diagonal elements of the covariance matrix are very noisy. Consequently, we will be forced to ignore the off-diagonal elements of C(i, j). We have tested the jackknife error estimator with lognormal mocks (see McDonald et al. 2006; Lidz et al. 2006) that approximately mimic the properties of the current data set. We generate 10, 000 mock realizations of a z = 3 data set and compare error bars estimated from the dispersion across the mock realizations with the jackknife error estimates. In the mock data, we find that neglecting the off-diagonal elements in the covariance matrix increases the average value of χ² by ≈1 for 14 degrees of freedom (the mock PDFs had 15 rather than 10 A_L bins), and so ignoring the off-diagonal elements is likely a good approximation. The jackknife estimates of the diagonal elements of the covariance matrix agree with direct estimates of the dispersion across the mock data to better than 20% on average, although the jackknife estimator sometimes underpredicts the errors in the tails of the PDF more severely. We provide tables of the wavelet PDF measurements in Tables 1–5.

We plot the measured wavelet PDF in the next section, but pause to consider first the impact of shot noise. The observed Lyα forest spectra are contaminated by random noise owing to Poisson fluctuations in the discrete photon count and around the mean night sky background count, as well as by random read-out noise in the CCD detector (see, e.g., Hui et al. 2001 for discussion). We need to consider how this noise impacts the wavelet PDF measurements.

In Appendix A, we derive estimates of the noise bias for measurements of the first two moments of the wavelet amplitude PDF. We apply these formulae here to estimate the impact of noise on the present measurements. On the larger smoothing scale, s_n = 69.7 km s⁻¹, we find that shot-noise bias is unimportant for our present data set. For example, at z = 3, applying the formulae of Appendix A, we find that the noise contamination to the mean wavelet amplitude is less than one-third of the 1σ statistical error on this quantity for our present data sample. Similarly, in this redshift bin and for this smoothing scale, we find that the wavelet amplitude variance is biased by random noise only at the ≈1% level. However, the shot-noise bias is not negligible on the smaller smoothing scale, s_n = 34.9 km s⁻¹. For instance, a quasar spectrum with S/N ≈ 50 at the continuum contributes a mean wavelet amplitude owing to noise of 〈|a^noise_n|²〉 ≈ (N/S)²/〈F〉 ≈ 6 × 10⁻⁴ at z ≈ 3 (Appendix A; Hui et al. 2001). This is comparable to the wavelet amplitude signal in the tail of the PDF in the favored hot IGM models (see Figure 3). The more significant noise contamination on the smaller smoothing scale owes to the rapid decline in signal power toward small scales. To guard against noise bias at the smaller smoothing scale, we cut spectra with S/N ⩽ 50 redward of Lyα from the sample used in the smaller smoothing scale measurement. We cut based on the red side noise, rather than using a noise estimate in the forest, to avoid introducing any possible selection bias. We estimate the red side S/N using spectral regions which lie between rest-frame wavelengths of 1250 Å and 1350 Å that are free of absorption lines. It is imperfect to define our cut based on the S/N redward of Lyα since the detector sensitivity varies with wavelength, but we use regions close to Lyα and so this simple cut should be adequate for present purposes. Further, we add noise to the mock spectra when we compare them with the measurement on the smaller smoothing scale (Section 4.5).

For the purpose of this project and related Lyα forest work, we have run a new suite of cosmological smoothed particle hydrodynamics (SPH) simulations using the simulation code Gadget-2 (Springel 2005). The simulations adopt a LCDM cosmology parameterized by n_s = 1, σ₈ = 0.82, Ω_m = 0.28, Ω_Λ = 0.72, Ω_b = 0.046, and h = 0.7 (all symbols have their usual meanings), consistent with the Wilkinson Microwave Anisotropy Probe (WMAP) constraints from Komatsu et al. (2009). Each simulation was started from z = 299, with the initial conditions generated using the Eisenstein & Hu (1999) transfer function. We ran several simulations to test the convergence of our results with boxsize, as well as mass and spatial resolution (see Appendix C). From these tests, we determined that the best choice simulation for the present project has a boxsize of L_box = 25 Mpc/h and N_p = 2 × 1024³ particles, and this run is the fiducial simulation in what follows. This simulation represents a fairly significant improvement in boxsize and resolution compared to most previous work (see Section 4.10 for details). It has approximately the gas mass recommended for resolution convergence in a recent study by Bolton & Becker (2009) and tracks over an order of magnitude more particles than the simulations of these authors. In each run, the softening length was taken to be 1/20th of the mean inter-particle spacing. In order to speed up the calculation, we chose an option in Gadget-2 that aggressively turns all gas at density greater than 1000 times the cosmic mean density into stars (e.g., Viel et al. 2004). Since the forest is insensitive to gas at such high densities, this approximation should not effect our results.

The simulations were run using the Faucher-Giguère et al. (2009) photoionizing background, which is an update of the Haardt & Madau (1996) model (see also Katz et al. 1996a; Springel & Hernquist 2003).¹⁷ The ionizing background was turned on at z = 7 in the simulations. This model for the ionizing background determines the photoheating and gas temperature in the simulation. We would like, however, to explore a wide range of thermal histories. In order to do this, we make an approximation. The approximation is to fix the fiducial ionization history to the Faucher-Giguère et al. (2009) model for the purpose of running the simulation and accounting for gas pressure smoothing, but to vary the temperature–density relation (Equation (1)) when constructing simulated spectra. This "post-processed spectra" approximation neglects the dependence of Jeans smoothing on the detailed thermal history of the IGM, but correctly incorporates thermal broadening for a given temperature–density relation model, parameterized by T₀ and γ. It also neglects the inhomogeneities in T₀ and γ expected during He ii reionization. Finally, by assuming a perfect temperature–density relation in constructing mock absorption spectra, we also neglect the impact of shock heating—which adds scatter to the temperature density relation (Hui & Gnedin 1997)—on the amount of thermal broadening. We caution against taking the results of these first pass, homogeneous temperature–density relation calculations too literally: if the IGM temperature is significantly inhomogeneous, these calculations provide only a crude approximation. The calculations are meant only to get a sense for whether the IGM is mostly at T₀ ≈ 20, 000 K, or instead at T₀ ≈ 10, 000 K, and to check whether large temperature inhomogeneities are present. We intend to make more detailed theoretical calculations in future work.

Although our measurements of the wavelet amplitude PDF are robust to uncertainties in fitting the quasar continuum, our interpretation of the measurements still relies on estimates of the mean transmitted flux. Specifically, we follow the normal procedure of adjusting the intensity of the simulated ionizing background in a post-processing step, so that the simulated mean transmitted flux (averaged over all sightlines) matches the observed mean flux. We assume here the best-fit measurements of Faucher-Giguère et al. (2008b), and subsequently explore the impact of uncertainties in the mean flux (Section 4.3). We adopt their estimates in Δz = 0.2 bins, and use their measurements that include a correction for metal line opacity based on Schaye et al. (2003), and a continuum-fitting correction (which accounts for the rarity of regions with nearly complete transmission, F = 1, at high redshift). The corresponding Faucher-Giguère et al. (2008b) measurements in our redshift bins are $(\bar{z}, \langle F \rangle) = (2.2, 0.849), (2.6, 0.778), (3.0, 0.680), (3.4, 0.566),$ and (4.2, 0.346). We output simulation data at every Δ_z = 0.5 between z = 4.5 and z = 2. In order to generate model wavelet PDFs at redshifts in between two stored snapshots, we measure the simulated wavelet PDF from each stored snapshot and linearly interpolate to find the PDF at the precise desired redshift.

Let us first compare the measured PDF in the different redshift bins for s_n = 69.7 km s⁻¹. The results of these calculations are shown in Figures 15–19. We start with a qualitative "chi-by-eye" assessment and provide more quantitative constraints in Section 4.6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The blue histogram with error bars in Figure 15 shows the measured PDF at $\bar{z}=4.2$, uncorrected for metal line contamination. We have not identified metal lines in the high redshift spectra contributing to this redshift bin, but we expect metal line contamination to have only a small effect on the wavelet PDF at this redshift and smoothing scale (see Appendix B). The model curves with T₀ ≈ 7, 500–10, 000 K and γ = 1.6 correspond roughly to models in which H i is reionized early, and He ii is not yet ionized. One expects a similarly low temperature in models in which each of H i, He i, and He ii are all ionized early. Interestingly, these models produce too many large wavelet amplitude regions and too few small wavelet amplitude regions compared to the data. The model curves with T₀ = 15, 000 K, and each of γ = 1.3 and γ = 1.6 are fairly close to the measurements, but overproduce slightly the high amplitude tail. These two curves are almost completely degenerate because the wavelet amplitude PDF is sensitive to the temperature over only a limited range in overdensity. At this redshift the measurements appear most sensitive to the temperature at densities near the cosmic mean, and so the models depend sensitively on T₀ but not on γ. The model with T₀ = 2 × 10⁴ K, γ = 1.3 is the best overall match to the data of the models shown, although it over-predicts the point near A_L ≈ 0.4 by more than 2.5σ. Finally, the model with T₀ = 2.5 × 10⁴ K seems to produce too many low wavelet amplitude regions, and too few high amplitude pixels. It is also interesting that the measured PDF is not much wider than the model PDFs. Taken at face value, this argues against the temperature field being very inhomogeneous at this redshift.

The results are tantalizing because they suggest the IGM is fairly hot with T₀ ≈ 15–20, 000 K at z = 4.2. This requires some amount of early He ii photoheating and/or H i reionization to end late and heat the IGM to a high temperature. If metal line contamination is in fact significant, this only strengthens the argument for a high temperature at $\bar{z}=4.2$: metal lines can only add power and increase the number of high wavelet amplitude regions. Similarly, the finite resolution of our numerical simulations causes us to underestimate the IGM temperature (see Appendix C). While we show that our results are mostly converged with respect to simulation resolution in Appendix C, convergence is most challenging at high redshift and this may lead to a small systematic underestimate in this redshift bin. On the other hand, we show in Section 4.3 that a cooler IGM model (T₀ ≈ 10, 000 K) can match the PDF measurement at this redshift if the true mean transmitted flux is 2σ higher than the best-fit value estimated by Faucher-Giguère et al. (2008b).

The measurements in our next redshift bin ($\bar{z}=3.4$) suggest the presence of even hotter gas in the IGM (Figure 16). At this redshift the best overall match is the model with T₀ = 2.5 × 10⁴ K, γ = 1.3. Even a fairly hot model with T₀ ≈ 2 × 10⁴ K, γ = 1.3 produces too few low wavelet amplitude regions and too many high amplitude ones. Models with lower temperatures are clearly quite discrepant. At this redshift, the measured PDF is a bit wider than the simulated ones. This might owe to temperature inhomogeneities, or it may indicate some metal line contamination since, as with the $\bar{z}=4.2$ data, we have not identified and excised metal lines in this redshift bin. In either of these cases, the measurements may allow for some even hotter gas at T₀ ≈ 3 × 10⁴ K.

The measurements at $\bar{z}=3.0$ indicate similarly hot gas (Figure 17). By this redshift, the average absorption in the forest is increased and the wavelet PDF is most sensitive to gas a little more dense than the cosmic mean, at roughly 1 + δ = Δ ≈ 2 for our method. This means that models that have a lower temperature at mean density (T₀), yet a steeper temperature–density relation (γ) give similar wavelet PDFs to models with higher T₀ and flatter γ at this redshift. This explains why the model curves with T₀ = 2.5 × 10⁴ K, γ = 1.3, and T₀ = 2 × 10⁴ K, γ = 1.6 are nearly identical to each other, as are the models with T₀ = 2 × 10⁴ K, γ = 1.3, and T₀ = 1.5 × 10⁴ K, γ = 1.6. At this redshift, the metal line correction appears fairly important: it shifts the peak of the PDF to lower amplitude and narrows the histogram somewhat (as seen by comparing the black dashed histogram and the blue solid histogram in the figure). The error bars are significantly larger for the metal-cleaned measurement than for the full measurement. This is mostly because we only have metal line identifications for some of the spectra in this bin and the metal-cleaned measurement hence comes from a smaller number of spectra, and also because we use a smaller portion of each spectrum after metal cleaning. The mean wavelet amplitude changes by less than the 1σ error bar as we vary f_m between f_m = 0.8 and f_m = 1, and so f_m = 0.95 is a conservative choice, and we hence stick to this choice throughout. After accounting for metal contamination, the model curves with T₀ = 2 × 10⁴ K, γ = 1.3, and T₀ = 1.5 × 10⁴ K, γ = 1.6 are somewhat disfavored. Again, the cooler models with T₀ = 7500–10, 000 K differ strongly with the measurement, regardless of the metal correction. The hottest model shown with T₀ = 3.0 × 10⁴ K, γ = 1.3 produces too many low-amplitude and two few high amplitude regions. The models with T₀ = 2.5 × 10⁴ K, γ = 1.3 and T₀ = 2.0 × 10⁴ K, γ = 1.6 are strongly degenerate and each roughly match the measured PDF.

Proceeding to lower redshift, the data at $\bar{z}=2.6$ disfavor some of the hotter IGM models (Figure 18). At this redshift, the models shown with T₀ = 3.0 × 10⁴ K, γ = 1.3; T₀ = 2.5 × 10⁴ K, γ = 1.3; T₀ = 2.0 × 10⁴ K, γ = 1.6 all produce too many low amplitude regions, and too few high amplitude ones. The other models shown with T₀ = 2.0 × 10⁴ K, γ = 1.3; T₀ = 1.5 × 10⁴ K, γ = 1.6, and T₀ = 1.0 × 10⁴ K, γ = 1.6 are closer to the measurements, although none of the models are a great fit. The cooler model with T₀ ≈ 7500 K is again a poor match to the measurement. The preference for somewhat more moderate temperatures at this redshift may result from cooling after He ii reionization completes at higher redshift.

Finally, the measurement in the $\bar{z}=2.2$ bin is shown in Figure 19. The general features are similar to the results at $\bar{z}=2.6$: the hotter models are clearly a poor match to the data, and there is some preference for cooler temperatures, although none of the models are a great match to the data. The models with (T₀, γ) = (2.0 × 10⁴ K, 1.3), (1.5 × 10⁴ K, 1.6), and (1.0 × 10⁴ K, γ = 1.6) are the closest matches of the models shown. At this redshift, the mean transmission is high (〈F〉 = 0.849), and the method is sensitive to somewhat overdense gas as a result. The similarity between the models with T₀ = 3.0 × 10⁴ K, γ = 1.3 and T₀ = 2.0 × 10⁴ K, γ = 1.6 suggests that the PDFs are most sensitive to densities around Δ = 3.9 at this redshift. We expect scatter in the temperature–density relation from shock-heating to be most important at this low redshift, especially since the wavelet PDF is becoming sensitive to the temperature of moderately overdense gas. This may be part of the reason for the poorer overall match between simulations, where the effects of shocks on T are ignored in post-processing and observations at this redshift. We will investigate this in more detail in the future.

In summary, our measurements appear to support a picture where the IGM is being heated in the middle of the redshift range probed by our data sample, with the temperature likely peaking between z = 3.0–3.8, before cooling down toward lower redshifts. The favored peak temperature appears to be around T₀ ≈ 25–30, 000 K, somewhat hotter than found by most previous authors (see Section 4.10), although consistent with theoretical expectations from photoheating during He ii reionization, especially if the quasar ionizing spectrum is on the hard side of the models considered by McQuinn et al. (2008, see their Figure 12).

In the previous section, we showed model wavelet PDFs for varying temperature–density relations, but left the underlying cosmology and mean transmitted flux fixed. Here we consider how much the wavelet PDF varies with changes in these quantities. As far as the underlying cosmology is concerned, we restrict our discussion to uncertainties in the amplitude of density fluctuations. Note that there is some (2σ level) tension between the amplitude of density fluctuations determined from the Lyα forest and WMAP constraints (Seljak et al. 2006).

In order to gauge the impact of uncertainties in the amplitude of density fluctuations, we generate mock spectra for a given model using simulation outputs of varying redshift. In particular, we consider a model at $\bar{z}=4.2$ with 〈F〉 = 0.346, T₀ = 2 × 10⁴ K, and γ = 1.3, which roughly matches the measured PDF. We generate mock spectra in this model from outputs at z_o = 3.5, 4.0, and 4.5. For the prediction in our fiducial cosmology, we linearly interpolate between wavelet PDFs generated from the z = 4.0 and z = 4.5 outputs. Using instead the model PDFs at z_o = 3.5 or 4.0 (with the mean transmitted flux fixed at 〈F〉 = 0.346)—in which structure formation is more advanced—should mimic a model with a higher amplitude of density fluctuations, while using the z = 4.5 snapshot should correspond to a model with smaller density fluctuations. Our fiducial model has σ₈(z = 0) = 0.82, roughly in between the preferred values inferred from the forest alone and that from WMAP-3 alone (Seljak et al. 2006). Using the outputs at z_o = 3.5, 4.0, and 4.5 for the $\bar{z} = 4.2$ mock spectra should roughly correspond to models with σ₈(z = 0) = 0.95, 0.85, and 0.78, respectively. The results of these calculations, shown in Figure 20, illustrate that the wavelet PDF is only weakly sensitive to the underlying amplitude of density fluctuations. The mean small-scale power is exponentially sensitive to the temperature, which is uncertain at the factor of ≈2 level, and so it is unsurprising that ≈10% level changes in the amplitude of density fluctuations have relatively little impact. The small effect visible in the plot is that the wavelet PDF shifts to smaller amplitudes for the outputs in which structure formation is more advanced. This likely owes to the enhanced peculiar velocities in models with larger density fluctuations, which suppress the small-scale fluctuations in the forest via a finger-of-god effect (e.g., McDonald et al. 2006). The impact of uncertainties in the amplitude of density fluctuations on the wavelet PDF is similarly small at other redshifts, and so we do not discuss this further here.

Zoom In Zoom Out Reset image size

The amplitude of fluctuations in the forest, and the wavelet PDF, are sensitive to the mean transmitted flux and uncertainties in this quantity impact constraints on the temperature from the PDF measurements. The mean transmitted flux partly determines the "bias" between fluctuations in the transmission and the underlying density fluctuations, with the bias increasing as the mean transmitted flux decreases. This impacts the small-scale transmission power spectrum, and the wavelet PDF, as well as fluctuations on larger scales. When the gas density is sufficiently high, and/or the ionizing background sufficiently low—i.e., when the mean transmitted flux is small—even slight density inhomogeneities produce absorption features, yielding large transmission fluctuations on small scales.

In the previous section, we adopted the best-fit values of the mean transmitted flux from Faucher-Giguère et al. (2008b), but now consider variations around these values. These authors provide estimates of the statistical and systematic errors on their mean transmitted flux measurements. Their 1σ errors at our bin centers are (z, 〈F〉 ± 1 − σ) = (2.2, 0.849 ± 0.017), (2.6, 0.778 ± 0.017), (3.0, 0.680 ± 0.02), (3.4, 0.566 ± 0.022), and(4.2, 0.346 ± 0.042). Their systematic error budget accounts for uncertainties in estimating metal line contamination, and for uncertainties in corrections related to the rarity of true unabsorbed regions at high redshift, among other issues. Nonetheless, there is some tension between the measurements of different groups. We refer the reader to Faucher-Giguère et al. (2008b) for a discussion.

Below z ≲ 4 uncertainties in the mean transmitted flux have a noticeable yet fairly small impact on our constraints. A typical example, in the $\bar{z} = 3.4$ redshift bin, is shown in Figure 21. The solid black line in the figure shows a model with T₀ = 2.5 × 10⁴ K, γ = 1.3 that adopts the best-fit value for the mean transmission, 〈F〉 = 0.566. The blue dashed line is the same model, but with the mean transmitted flux shifted up from the central value by 1σ. This reduces the amplitude of transmission fluctuations in the model and shifts the wavelet PDF towards slightly lower amplitudes. Reducing the transmission by 1σ has the opposite effect of boosting the typical wavelet amplitudes slightly, as illustrated by the red dotted line in the figure. While the uncertainty in the mean transmitted flux can shift the preferred temperature around slightly, the effect at this redshift is relatively small and has little impact on our main conclusions. For example, a cooler IGM model with T₀ = 1.0 × 10⁴ K and γ = 1.6 still differs greatly from the PDF measurement, even after assuming a mean transmitted flux that is 2σ higher than the central value. This is demonstrated by the magenta line in Figure 21.

Zoom In Zoom Out Reset image size

The impact of uncertainties in the mean transmitted flux is more important in our highest redshift bin, at $\bar{z}=4.2$. The impact is larger at this redshift both because data samples are smaller and the fractional error on the mean transmitted flux is larger at this redshift, and because the wavelet amplitudes are more sensitive to the mean transmission once the transmission is sufficiently small. We repeat the exercise of the previous figure at $\bar{z}=4.2$ and present the results in Figure 22. In this case, the model that roughly goes through the PDF measurement with our best-fit mean transmitted flux has T₀ = 2.0 × 10⁴ K and γ = 1.3. After shifting the mean transmitted flux up in this model by 1σ it produces too many low wavelet amplitude regions, and too few high amplitude ones, in comparison to the measurement. Indeed, at this redshift, even the cooler IGM model with T₀ = 1.0 × 10⁴ K, γ = 1.6 will pass through the measurement after a 2σ upward shift in the mean transmitted flux. In other words, accounting for uncertainties in the mean transmitted flux, the cool IGM model with T₀ = 1.0 × 10⁴ K, γ = 1.6 can only be excluded at roughly the 2σ level.

Zoom In Zoom Out Reset image size

Furthermore, systematic concerns with direct continuum-fitting are most severe at high redshift, and the agreement between different measurements, while generally good at lower redshifts, is marginal above z ≈ 4 or so (Faucher-Giguère et al. 2008b). Direct continuum measurements must correct for the fact that there are few genuinely unabsorbed regions at high redshift, which can cause one to systematically underestimate the mean transmitted flux. Part of the disagreement can be traced to the fact that some of the measurements in the literature do not make this important correction. Since Faucher-Giguère et al. (2008b) make a correction using cosmological simulations, we consider their measurement to be more reliable than many of the other previous measurements. However, McDonald et al. (2006) constrain the mean transmitted flux based on a multi-parameter fit to their Sloan Digital Sky Survey (SDSS) power spectrum measurements, which should be immune to this concern. Their best fit to the redshift evolution of the mean transmitted flux gives 〈F〉 = 0.41 at z = 4.2. This disagrees with the Faucher-Giguère et al. (2008b) measurement at this redshift by 1.6σ. The overall disagreement is in fact more severe than this because there is a similar level of disagreement in neighboring redshift bins. Dall'Aglio et al. (2008) also perform a direct continuum fit, correct for the rarity of unabsorbed regions at high redshift with a different methodology, and find a best fit to the redshift evolution of the opacity of 〈F〉 = 0.40 at z = 4.2. Again, this measurement is in tension with the measurement we adopt. Adopting either of these measurements for the best-fit mean transmitted flux would favor a cooler IGM temperature. See Faucher-Giguère et al. (2008b) for further comparison and discussion regarding different measurements of the mean transmitted flux in the literature.

The measured PDF in the $\bar{z}=3.4$ redshift bin requires hot (T₀ ≳ 20, 000 K) gas. Interestingly, the PDF in this redshift bin is somewhat broader than the theoretical model curves, which assume a homogeneous temperature–density relation. This may be the result of uncleaned metal line contamination, but a more interesting possibility is that the wide measured PDF indicates temperature inhomogeneities from ongoing He ii reionization. We argued in Section 2.3 that the precise choice of large-scale smoothing, L, should be relatively unimportant. Nevertheless, to further explore the exciting possibility that the data indicate temperature inhomogeneities in this redshift bin, we measure the PDF for a few additional choices of L and compare with theoretical models.

The results of these calculations are shown in Figure 23. In addition to our usual large-scale smoothing of L = 1000 km s⁻¹, we also compare simulated and observational wavelet PDFs for L = 200, 2000, and 5000 km s⁻¹. Here we use 15 logarithmically spaced A_L bins for the PDF measurement, rather than 10 as in the previous sections, to increase our sensitivity to any bi-modality in the PDF. The mean of the model curves with different smoothing scales is of course fixed, while the width of the PDF increases with decreasing smoothing scale (see Section 2.3, Figure 6). At all smoothing scales, the simulated model with T₀ = 25, 000 K and γ = 1.3 is the best overall match to the data. The fit is poorest at L = 2000 km s⁻¹, but it is not clear precisely how to interpret this since the model is a formally poor fit at each smoothing scale. There does appear to be a slight, yet tantalizing, hint that the PDF is bimodal on large smoothing scales: this trend is most apparent at L = 2000 km s⁻¹ and L = 5000 km s⁻¹. This may be a first indication of temperature inhomogeneities from ongoing He ii reionization, or it may be the result of uncleaned metal line contamination, as the abundance of metals can vary significantly on large smoothing scales. It will be interesting to revisit this measurement with larger data samples in the future.

Zoom In Zoom Out Reset image size

In the previous sections we found that our results at s_n = 34.9 km s⁻¹ are quite susceptible to metal-line contamination and somewhat to shot-noise bias. Because of this, we will not presently use the results at this smoothing scale in constraining the thermal history of the IGM. Nevertheless, as a consistency check we compare here the measured wavelet PDF at this smoothing scale with simulated models.

As mentioned previously, to guard against shot-noise bias, we cut spectra with a (red-side) S/N ⩽ 50 and add random noise to the mock spectra. Provided we cut out the low S/N data, the random noise mainly impacts only the low wavelet amplitude tail by decreasing the number of very low amplitude wavelet regions. We add Poisson distributed noise to the mock spectra, assuming that the noise is dominated by Poisson fluctuations in the photon counts from the quasar itself. We have experimented with incorporating Poisson distributed sky noise, and Gaussian random read-noise, and find qualitatively similar results at fixed noise level. We estimate the average wavelet amplitude in the forest contributed by noise (after our S/N cut) as described in Appendix A, and find that it corresponds to S/N ≈ 70, per 4.4 km s⁻¹ pixel at the continuum for the $\bar{z}=3.0$ bin. In Figure 24, we compare some example model PDFs with the measurements and find results gratifyingly close to those at larger smoothing scale. In particular, the model with T₀ = 2.0 × 10⁴ K, and γ = 1.6 at $\bar{z}=3.0$ that roughly matched the measurement on larger scales, matches the PDF on this smaller scale as well. For contrast, we show a hotter and a colder IGM model which are again a poor match. At $\bar{z}=2.2$ and $\bar{z}=2.6$ the results are similar to the previous ones, suggesting a cooler IGM at these redshifts. Comparing the blue and black dashed histograms, it is clear that the metal contamination correction is quite important at this scale and we do not use these results in what follows.

Zoom In Zoom Out Reset image size

We have also compared the s_n = 34.9 km s⁻¹ wavelet PDF in the two highest redshift bins—where we have not identified metal lines—with model PDFs. The measured PDF at z = 3.4 looks similar to the T₀ = 2.5 × 10⁴ K, γ = 1.3 model that we previously identified as the best general match of our example models at s_n = 69.7 km s⁻¹, except with a fairly prominent tail toward high wavelet amplitudes. We expect more significant metal contamination at this smoothing scale (Appendix B), and so this is in line with our expectations. Indeed, the tail toward high wavelet amplitude looks similar to the one in the top panel of Figure 33. Similar conclusions hold at z = 4.2, except the agreement without excising metals is better, likely owing to the smaller impact of metals at this redshift (Appendix B).

In this section, we perform a preliminary likelihood analysis, in order to provide a more quantitative constraint on the thermal history of the IGM from the wavelet measurements. We confine our analysis to a three-dimensional parameter space, spanning a range of values for T₀, γ, and 〈F〉. The results of the previous section suggest that CDM models close to a WMAP-5 cosmology should all give similar wavelet PDFs, and so it should be unnecessary to vary the cosmological parameters in this analysis. In order to facilitate this calculation, we adopt here an approximate approach to cover the relevant parameter space. We generate the wavelet PDF for a range of models by expanding around a fiducial model in a first-order Taylor series (see Viel & Haehnelt 2006 for a similar approach applied to SDSS flux power spectrum data). In particular let p denote a vector in the three-dimensional parameter space. Then we calculate the wavelet PDF at a point in parameter space assuming that:

$$\begin{eqnarray} P(A_L,\bm {p}) = && P\left(A_L,{\bm { p} ^0}\right) \nonumber \\ && +\sum _{i=1}^3 \frac{\partial P(A_L, p_i)}{\partial p_i} \bigg |_{\bm {p} =\bm {p } ^0} \left(p_i - p^0_i\right). \end{eqnarray} \tag{ 11 }$$

Although inexact, this approach suffices to determine degeneracy directions, approximate confidence intervals, and the main trends with redshift. We use the results of the previous section to choose the fiducial model to expand around: at each redshift we choose the best match of the example models in the previous section as the fiducial model. Using the Taylor expansion approximation of Equation (11), we then estimate the wavelet PDF for a large range of models, spanning T₀ = 5000–35, 000 K, γ = 1.0–1.6, and 〈F〉 = F_c ± 3σ_F (subject to a 〈F〉 prior). Here F_c denotes the central value from Faucher-Giguère et al. (2008b), and σ_F denotes their estimate of the 1σ uncertainty on the mean transmitted flux. For each model PDF in the parameter space, we first compute χ² between the model and the wavelet PDF data, ignoring off-diagonal terms in the co-variance matrix. We then add to this χ² an additional term to account for the difference between the model mean transmitted flux and the best-fit value of Faucher-Giguère et al. (2008b). Finally, we marginalize over 〈F〉 (subject to the above prior based on the results of Faucher-Giguère et al. 2008b) to compute two-dimensional likelihood surfaces in the T₀–γ plane at each redshift, and marginalize over γ to obtain reduced, one-dimensional likelihoods for T₀. We assume Gaussian statistics, so that 1σ (2σ) two-dimensional likelihood regions correspond to Δχ² = 2.30(6.17), while one-dimensional constraints correspond to Δχ² = 1(4).

The best-fit models at z = 4.2, 3.4, 3.0, 2.6, and 2.2 have χ² = 9.5, 19.8, 5.7, 8.0, and 23.1 respectively for 7 degrees of freedom (10 A_L bins minus 1 constraint since the PDF normalizes to unity, minus two free parameters). The fits at z = 4.2, 3.0, and 2.6 are acceptable, while the χ² values in the z = 3.4 and z = 2.2 bins are high (p-values of 6 × 10⁻³ and 2 × 10⁻³ respectively). The poor χ² in these redshift bins results because the measured PDFs are broader than the theoretical models in these bins, as discussed previously. We will nevertheless consider how χ² changes around the best-fit models in these redshift bins, although we caution against taking the results too literally.

The constraints from these calculations are shown in Figures 25 and 26. They are qualitatively consistent with the example models shown in the previous section. The degeneracy direction of the constraint ellipses results because the z = 4.2 measurements are sensitive only to the temperature close to the cosmic mean density, while the lower redshift measurements start to constrain only the temperature of more overdense gas. The best-fit model at z = 4.2 has T₀ ≈ 20, 000 K, but uncertainties in the mean transmitted flux allow cooler models with T₀ ≈ 10, 000 K at ≈2σ, as discussed previously. The z = 3.4 measurements indicate the largest temperatures, and require that T₀ ≳ 20, 000 K at 2σ confidence. The lower redshift measurements, particularly that at z = 2.6, generally favor cooler temperatures although at only moderate statistical significance.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 26 shows ±1 and 2σ error bars on the temperature at mean density after marginalizing over γ and 〈F〉. We conservatively allow γ to vary over γ = 1.0–1.6, even though γ ≳ 1.2 is expected theoretically (McQuinn et al. 2008). If we enforced a prior that γ be steeper than 1.2, then the results at z ≲ 3.4 would disfavor some of the higher T₀ models. The T₀ results are consistent with the IGM temperature falling off as T₀ ∝ (1 + z)² below z = 3.4, i.e., below this redshift the temperature evolution appears consistent with simple adiabatic cooling owing to the expansion of the universe. Theoretically, we expect the temperature fall-off to be similar, but slightly slower, than the adiabatic case just after reionization with the temperature evolution eventually slowing owing to residual photoionization heating (Hui & Gnedin 1997). The statistical errors are however still large, and a flat temperature evolution is also consistent with the T₀ constraints, although this case is disfavored theoretically (see below). Note also that enforcing a γ ⩾ 1.2 prior would disfavor the high T₀ models that are otherwise allowed at z = 2.2 and z = 2.6, strengthening the case for cooling below z ≈ 3.4.

Moreover, the high temperatures at z = 3.0 and z = 3.4 suggest recent He ii photoheating. To illustrate this point, we show several example thermal history models in Figure 26, considering both cases without any He ii photoheating, and ones in which H i/He i/He ii are all reionized together at high redshift (z ⩾ 6). The upper blue dashed line is a late H i reionization model (z_r = 6), with a high temperature at reionization (T_r = 3 × 10⁴ K), and a hard spectrum near the H i/He i ionization thresholds (with a specific intensity near threshold of J_ν ∝ ν^−α and α = 0). This case should roughly indicate the highest possible temperature without He ii photoheating over the redshift range probed. Note that this is a rather extreme situation, since even if reionization completes as late as z = 6, much of the volume will be reionized significantly earlier (e.g., Lidz et al. 2007). The lower blue dashed line is an early reionization model (z_r = 12 and α = 2) that approximately indicates the lowest plausible temperature without He ii photoheating.

Finally, perhaps the most interesting case is the black dot-dashed line which shows a model in which H i/He i/He ii are all reionized together at z = 6. Here we assume that the temperature at reionization is T_r = 3 × 10⁴ K, since atomic hydrogen line cooling should keep the temperature less than this when all three species are ionized simultaneously (Miralda-Escudé & Rees 1994; Abel & Haehnelt 1999; A. Lidz et al. 2010, in preparation). The temperature after reionization depends on the ionizing spectrum, which determines the amount of residual photoheating. The curve here adopts a quasar-like spectrum, reprocessed by intervening absorption, to give α_H i = 1.5 near the H i ionization threshold and α_He ii = 0 near the He ii ionization threshold (Hui & Haiman 2003). This case is hence similar to the other z = 6 reionization model, except with the addition of residual He ii photoheating. Each of the examples considered gives too low a temperature in the z = 2.2, z = 3.0, and z = 3.4 redshift bins, particularly at z = 3.4 and z = 3. One can further ask how hard the post-reionization ionizing spectrum would need to be to give a thermal asymptote as large as ≈20, 000 K. For a power-law spectrum we find, using the thermal asymptote formula of Hui & Haiman (2003), that an implausibly hard spectrum with α≲−0.73 is required to match the 2σ lower limit on the z = 3.4 temperature. In fact, there is evidence that galaxies rather than quasars produce most of the ionizing background at z ≳ 3 (e.g., Faucher-Giguère et al. 2008b), and so assuming even a quasar-like spectrum likely overestimates residual photoheating for plausible early He ii reionization models. In summary, although the errors allow the possibility of a slow temperature evolution and T₀ ≈ 20, 000 K, this temperature is higher than expected from residual photoheating long after reionization.

The simplest interpretation is that He ii reionization heats the IGM, with the process completing near z ≈ 3.4, at which point there is relatively little additional heating and the universe expands and cools. The redshift extent over which the heat input occurs is, however, not well constrained by our present measurement. Clearly the large error bars on the measurements still leave room for other possibilities. For example, models in which He ii reionization completes a bit later at z ≈ 3—or perhaps even as late as z ≈ 2.7 as favored by a recent analysis of He ii Lyα forest data by Dixon & Furlanetto (2009)—or earlier at z ≈ 4 are likely consistent with our present measurements given the large error bars. We will consider this further in future work. Finally, other heating mechanisms may be at work in addition to photoionization heating.

Recently, Bolton et al. (2008), Becker et al. (2007), and Viel et al. (2009) have suggested that measurements of the Lyα flux PDF favor an inverted temperature density relation (γ < 1), i.e., situations where low density gas elements are hotter than overdense ones. Bolton et al. (2008) and Viel et al. (2009) construct simulated models with inverted temperature–density relations by adding heat into the simulations in a way that depends on the local density, i.e., on the density smoothed on the Jeans scale. This particular case for an inverted temperature–density relation seems unphysical to us since heat input from, e.g., reionization, should be coherent on much larger scales. Nonetheless, we can consider this as a phenomenological example that the flux PDF data favor and examine the implications of these models for the small-scale wavelet amplitudes. Theoretically, Trac et al. (2008) and Furlanetto & Oh (2009) find that hydrogen reionization does produce a weakly inverted temperature–density relation. This effect is driven by the tendency for large-scale overdensities to reionize hydrogen first, coupled presumably with the small correlation between the overdensity on large scales and that on the Jeans scale. On the other hand, McQuinn et al. (2008) find that He ii reionization leads to a non-inverted equation of state with γ ≈ 1.3 in the midst and at the end of He ii reionization. We refer the reader to this paper for further discussion.

To explore this, we generate mock spectra and measure wavelet amplitudes for several inverted temperature–density relation models and compare with our z = 3 measurements. As before, we are considering the impact of the temperature in a post-processing step, and so we are not accounting for differences between the gas pressure smoothing in the inverted models and that in the simulation. Likewise, we incorporate thermal broadening assuming a perfect temperature–density relation, and so the impact of scatter in the temperature-density relation is ignored in this part of the calculation. We consider inverted temperature density relations with a power-law index of γ = 0.5, close to the value suggested by Bolton et al. (2008) and Viel et al. (2009) from their flux PDF measurements near z = 3. The results of these calculations are shown in Figure 27. These cases also roughly match the observed PDF, but require a very high temperature at mean density of T₀ ≈ 40–45, 000 K. The reason for this is that the wavelet PDF measurements are sensitive mostly to the temperature around a density of Δ ≈ 2 at this redshift. In the previous section we found that models with, for example, T₀ ≈ 25, 000 K and γ = 1.3 roughly match the data. A model with an inverted temperature density relation (γ = 0.5) produces the same temperature at a density of Δ ≈ 2 only for a much higher temperature (at mean density) of T₀ ≈ 45, 000 K. The figure suggests that the expected degeneracy between T₀ and γ indeed extends to even these inverted temperature–density relations. Hence one can fit the measurements with a very high T₀, small γ model, although the inverted cases produce slightly wider PDFs. While these can fit the data, the high required temperatures seem unlikely to us, and we disfavor inverted models for this reason.

Zoom In Zoom Out Reset image size

Bolton et al. (2008) and Viel et al. (2009) found that inverted models with substantially smaller temperature at mean density match their flux PDF measurements. On the other hand, Viel et al. (2009) did a joint fit to the flux PDF and the SDSS flux power spectrum from McDonald et al. (2006). Recall that the SDSS measurements are sensitive only to the large-scale flux power spectrum (k ≲ 0.02 s km⁻¹), and thus depend on IGM parameters differently than the small-scale wavelet measurements explored here. Their joint fit requires high T₀ for cases with inverted temperature–density relations, similar to our conclusions from a different type of measurement. There thus appears to be some tension with the flux PDF measurement, which may reflect systematic errors in one or more of the measurements and/or the modeling. We intend to consider this further in future work.

Let us further consider the implications of our measurements for the presence or absence of temperature inhomogeneities in the IGM. In most redshift bins, the measured PDF has comparable width to the simulated PDFs, which assume a perfect temperature–density relation.¹⁸ The possible exceptions are the $\bar{z}=3.4$ bin (where metal contamination is a possible culprit) and the $\bar{z}=2.2$ bin (where scatter from shocks may be most important). One might wonder if the widths of the wavelet PDFs are too small to be compatible with ongoing or recent He ii reionization, which is presumably a fairly inhomogeneous process. A related question regards the precise meaning of our temperature constraints in the presence of inhomogeneities: which temperature do we measure exactly—the mean temperature, the minimum temperature, etc.? We intend to address these issues in detail in future work, but we outline a few pertinent points here. In this discussion, we draw on the results of McQuinn et al. (2008).

The first point is that temperature inhomogeneities during He ii reionization, while likely important, are smaller than one might naively guess. McQuinn et al. (2008) emphasized the importance of hard photons, with long mean free paths, for He ii photoheating: much of the heating during He ii reionization by bright quasars occurs far from sources, rather than in well-defined "bubbles" around ionizing sources. This is quite different than during H i/He i reionization by softer stellar sources, where the ionizing photons have short mean free paths and heating does occur within well-defined bubbles. Since the hard photons have long mean free paths, and a "background" radiation field from multiple sources needs to be built up before these photons appreciably ionize and heat the IGM, the heating is much more homogeneous than might otherwise be expected. The softer photons, typically absorbed in bubbles around the quasar sources, only heat the IGM by δT ≲ 7000 K. Consider the temperature PDFs in Figure 11 of McQuinn et al. (2008). This figure illustrates that by the time any gas is heated significantly, there are very few completely cold regions left over in the IGM: the temperature field is more homogeneous than might be expected.

Simplified models with discrete ≈30, 000 K bubbles around quasar sources and a cooler IGM outside (e.g., Lai et al. 2006) are hence not realistic, and overestimate the temperature inhomogeneities. In the McQuinn et al. (2008) simulations, the temperature inhomogeneities peak at a level of σ_T/〈T〉 ≈ 0.2, which is reached in the early phases of He ii reionization. For contrast, a hypothetical two-phase hot/cold IGM with hot bubbles that are three times as hot as a cooler background IGM, gives a more substantial peak fluctuation level of σ_T/〈T〉 = 0.58, reached when the hot bubbles fill 25% of the IGM. In the midst of He ii reionization, the McQuinn et al. (2008) simulations predict roughly 10% level temperature fluctuations on large scales from inhomogeneous He ii heating. This level of scatter may be hard to discern with our existing measurement.

To illustrate this, we compare the $\bar{z}=4.2$, s_n = 69.7 km s⁻¹ measurement to a simplified and extreme two phase model. This redshift bin probes extended stretches of spectrum along just two lines of sight. Imagine a model where one line of sight passes entirely through cold regions of the IGM with T₀ = 10⁴ K, γ = 1.6, while the other line of sight passes entirely through hot regions with T₀ = 2.5 × 10⁴ K, γ = 1.3. This is a contrived example since each sightline probes hundreds of co-moving Mpc, and so each sightline should in reality probe a mix of temperatures, but this simple case nonetheless illustrates the challenge of detecting temperature inhomogeneities. For simplicity, in this hypothetical model we imagine that each line of sight probes an equal stretch through the IGM so that the wavelet PDF is a fifty–fifty mix of the hot and cold models. In this scenario, the mean IGM temperature is 〈T〉 = 17, 500 K and the fluctuation level is σ_T/〈T〉 ≈ 0.43, i.e., substantially larger than we expect. The wavelet PDF in this model is shown in Figure 28. This simple model clearly produces too broad a PDF, but it is also apparent that smaller, likely more realistic, levels of inhomogeneity will be hard to distinguish with the existing data. For example, an inhomogeneous model with fewer cold regions than in the hypothetical two-phase model would agree with the measurement. Indeed, the data may even favor slightly inhomogeneous models, but we leave exploring this to future work. In particular, note that the homogeneous theoretical models have poor χ² values at z = 2.2 and z = 3.4 (Section 4.6). The z = 4.2 and z = 3.4 data, which may be in the midst of He ii reionization, and which are sensitive to the temperature near the cosmic mean density, are the best redshift bins for further exploration. Provided the inhomogeneities are relatively small, as suggested by the measurements in most redshift bins, ambiguities in which temperature we constrain precisely are unimportant, and our temperature estimates should be accurate.

Zoom In Zoom Out Reset image size

Another possible issue, related to the discussion in Section 2.3, is that the one-dimensional nature of the Lyα forest may obscure detecting temperature inhomogeneities from He ii reionization. Consider the three-dimensional power spectrum of temperature fluctuations in Figure 10 of McQuinn et al. (2008). There is a large-scale peak in the three-dimensional power spectrum, owing to inhomogeneous heating, and a prominent small-scale ramp-up that results from the temperature–density relation and small-scale density inhomogeneities. The large-scale peak in the power spectrum is essentially the signal we are after, while the small-scale ramp-up is noise as far as extracting inhomogeneities is concerned. However, the one-dimensional temperature power spectrum may be more relevant than the three-dimensional one for absorption spectra. In the one-dimensional temperature power spectrum, high-k transverse modes, which are dominated by the small-scale ramp-up, will be aliased to large scales, swamping the temperature inhomogeneities. This argument is imperfect though, since the one-dimensional temperature power spectrum is not exactly the relevant quantity either: absorption spectra are insensitive to the temperature of large overdensities, which regardless produce saturated absorption. It will be interesting to consider this further in the future, and to consider the potential gains from cross-correlating the wavelet amplitudes of pairs of absorption spectra.

A final issue, particularly relevant in the highest redshift bin, is that the temperature inhomogeneities may depend on the timing and nature of hydrogen reionization. The temperature contrast between regions with doubly ionized helium and those in which only H i/He i are ionized depends on when hydrogen (and He i) reionized. Specifically, the temperature contrast between H ii/He ii and H ii/He iii regions will be reduced if hydrogen is reionized late to a high temperature, and increased if hydrogen reionizes early to a smaller temperature. Moreover, heating from hydrogen reionization will itself be inhomogeneous (e.g., Cen et al. 2009). Extending the measurements in this paper to higher redshift can help disentangle the impact of hydrogen and helium photoheating. Further modeling will also be helpful.

As mentioned previously, a shortcoming of our modeling thus far is that we have run only a single simulated thermal history in describing the gas density distribution in the IGM: we vary the thermal state of the gas only as we construct mock absorption spectra and incorporate thermal broadening. Similar approximations are common in the Lyα forest literature (although see Schaye et al. 2000). The gas density distribution is sensitive to the full thermal history of the IGM (Gnedin & Hui 1998) and so properly accounting for a range of thermal histories requires running many simulations. We expect thermal broadening to be more important than Jeans smoothing for our measurements, since thermal broadening directly smooths the optical depth field and results in a roughly exponential decrease in small-scale flux power (Zaldarriaga et al. 2001), while Jeans smoothing acts on the three-dimensional gas distribution and has a less direct impact.

In order to investigate the impact of the post-processing approximation on our results we have run several additional simulations with varying thermal histories. Our present aim is to check that this approximation has only a minor impact on our main results, and we defer a more detailed study to future work. Our additional simulations track 2 × 256³ particles in a 12.5 Mpc/h simulation box. We compare these results with those from our fiducial thermal history in an identical boxsize with the same particle number and initial conditions. Our convergence tests (Appendix C) suggest that these runs are inadequate for a detailed comparison to data, but they should suffice to gauge the relative importance of Jeans smoothing. The Faucher-Giguère et al. (2009) ionizing background turns on in our simulation at z = 7 and yields a slowly varying temperature–density relation at z = 2–4 with T₀ ≈ 11, 000 K, γ ≈ 1.6 at z = 3. This calculation likely underestimates the photoheating during He ii reionization (e.g., Miralda-Escudé & Rees 1994; Abel & Haehnelt 1999; McQuinn et al. 2008; Faucher-Giguère et al. 2009), and produces a smaller temperature than suggested by our measurements, and so it may underestimate the amount of Jeans smoothing.

In order to check this, we artificially boost the photoheating rates (of H i/He i/He ii) at fixed photoionization rate in our simulations by a factor of 3 for all redshifts below a redshift z_heat, varying z_heat across different simulation runs. While this procedure does not produce a very realistic thermal history, it provides a simple way of making the gas hotter and increasing the amount of gas pressure smoothing. We have run simulations with z_heat = 3.5, 4.0, 4.5, 5.5, and 7.0. Varying z_heat is important because it takes the gas some time to respond to prior heating: our measurement is sensitive to gas near the cosmic mean density in which case the sound crossing time is comparable to the Hubble time (Gnedin & Hui 1998; Schaye 2001). The gas in runs with z_heat close to z = 3 has had less time to adjust to the boost in photoheating and so we expect less Jeans smoothing in these runs. The simulated temperature at mean density for z_heat = 7.0 is about a factor of 2 larger than in our fiducial run at z = 3, and the temperature–density relation has a similar slope. These values are comparable to that suggested by our z = 3 temperature measurement.

With these simulation runs in hand, we produce (z = 3) mock spectra in our usual post-processing step assuming T₀ = 2 ×  10⁴ K, γ = 1.6 so that any differences result from differences in the amount of Jeans smoothing in the different runs. The resulting wavelet PDFs (for s_n = 69.7 km s⁻¹) are shown in Figure 29. The impact is quite small for the z_heat = 3.5, 4.0 models, moderate for z_heat = 4.5, and more significant for the models with earlier heating boosts. The results converge when z_heat is sufficiently high, with the z_heat = 5.5 and 7.0 simulations giving similar wavelet PDFs. Our interpretation is that the gas in the recent heating models has not had enough time to respond to the enhanced photoheating and is hence distributed similarly to the model without enhanced photoheating.¹⁹ Visual inspection of the gas density field along lines of sight through the various runs also supports this interpretation.

Zoom In Zoom Out Reset image size

We can further use these calculations to estimate any bias in our post-processing approximation at z = 3. If the level of Jeans smoothing in the real universe is similar to our z_heat = 4.5 model then the temperature estimated in our fiducial calculation is biased high by ≈2500 K, and by ≈5000 K if the level of Jeans smoothing resembles that in the z_heat ⩾ 5.5 models. If the level of Jeans smoothing is more like the z_heat ⩽ 4. models then the bias should be quite negligible at z = 3. In our favored interpretation that He ii reionization completes near z ≈ 3.4 we would expect our post-processing approximation to be quite good at z ⩾ 3, unless the heating is very extended in redshift. We hence expect our measurement of the peak temperature near z ≈ 3–4 to be robust to our imperfect modeling of Jeans smoothing. At lower redshift, where the gas has had more time to respond to prior heating, our approximation should be less good, although likely still within the present measurement errors (e.g., a 5000 K bias at z = 2.2 would be a 1.3σ bias at this redshift). These conclusions are at least broadly consistent with previous works which have argued that Jeans smoothing is sub-dominant, although not entirely negligible, for flux power spectrum (Zaldarriaga et al. 2001; Peeples et al. 2010a) and linewidth measurements (Schaye et al. 2000). Observational studies of the absorption spectra of close quasar pairs may help disentangle the effects of thermal broadening and Jeans smoothing (e.g., Peeples et al. 2010b).

A detailed comparison with previous measurements is difficult since our methodology differs from that of most previous work. Instead, we will simply compare the bottom line, and make a few remarks about the differences. Figure 30 shows our constraints on T₀(z), compared to the results of Schaye et al. (2000), Ricotti et al. (2000), McDonald et al. (2001), and Zaldarriaga et al. (2001). It is encouraging that some of the main trends are similar across all of the measurements: for example, all of the measurements favor a fairly hot IGM near z ≈ 3. In this sense, our work reinforces the previous results. There are differences in the details, however: the peak temperatures in Ricotti et al. (2000) and Schaye et al. (2000) are reached at lower redshift than in our analysis. The McDonald et al. (2001) and Zaldarriaga et al. (2001) results are, on the other hand, flat as a function of redshift, although they adopt wide redshift bins and may average over any temperature increase. Our measurements are also fairly consistent with a flat temperature evolution given the large error bars on our measurements. Our results mostly favor higher temperatures than the previous measurements, particularly the high redshift points of Schaye et al. (2000), although the differences still lie within the error bars.²⁰

Zoom In Zoom Out Reset image size

Although our results broadly support the main conclusions of previous authors, a possible reason for some of the differences is related to improvements in simulations of the forest over the past decade or so. In Appendix A, we found that our method—and we suspect related methods—require fairly large simulation volumes and high mass and spatial resolution, particularly at high redshift (see also, e.g., Bolton & Becker 2009). The requisite particle number, while achievable today, was of course prohibitive for past studies. Indeed, this was one of our motivations for revisiting the temperature measurements. While some of the previous studies varied simulation resolution and boxsize, they often considered only a single additional run, which may have been inadequate to fully assess convergence. Finite resolution, in particular, can bias temperature estimates low.

It is instructive to compare our fiducial simulation with a boxsize of L_b = 25 Mpc/h and N_p = 2 × 1024³ particles to the main runs of previous work. Schaye et al. (2000) used a (L_b, N_p) = (2.5 Mpc/h, 2 × 64³) SPH simulation, Ricotti et al. (2000)'s main runs were (2.56 Mpc/h, 2 × 256³) HPM calculations, McDonald et al. (2001) used a Eulerian hydrodynamic simulation with L_b = 10 Mpc/h and 288³ cells, and Zaldarriaga et al. (2001) used a dark matter only simulation with L_b = 16 Mpc/h and 128³ dark matter particles. Given the differences between methods, we will not try to estimate the impact of systematic errors from finite boxsize and resolution on previous results. Indeed, many of these previous studies did at least partly test their results for convergence with simulation boxsize and resolution and so the requirements of these previous methods may be less stringent than ours. However, it is clear that increases in computing power allow us to do a much better job with respect to boxsize and resolution than previous work. Finally, improved estimates of the mean transmitted flux (Faucher-Giguère et al. 2008b), and improved masking of metal lines, may also contribute to some of the differences between our results and previous work.

In comparison to this work, a strength of some of the previous studies is that they break the T₀–γ degeneracy by examining the column density dependence of the b-parameter distribution (Schaye et al. 1999, 2000; Ricotti et al. 2000; McDonald et al. 2001). In principle we may be able to extend our analysis by examining the correlation between A_L and δ_F which is directly analogous to the b-parameter column density correlation, yet does not require line fitting. We defer this study to future work.