Revealing Reionization with the Thermal History of the Intergalactic Medium: New Constraints from the Lyα Flux Power Spectrum

The epoch of hydrogen reionization represents one of the most dramatic phases of evolution of the universe. During this period, the UV radiation from the first luminous sources reionized the neutral hydrogen (and He i) atoms in the diffuse intergalactic medium (IGM), driving the transition from a neutral to a highly ionized universe. Understanding sources and timing of this transformation can reveal crucial information on the properties of the first objects and the environment in which they were formed. When and how reionization happened therefore remains a primary subject of interest in extragalactic astrophysics (for a review, see Becker et al. 2015a).

The most direct probes of the highly ionized IGM have been obtained from observations of intergalactic Lyα absorption along the lines of sight to high-redshift quasars. Measurements of Lyα transmission along some lines of sight suggest that reionization was largely complete by z ∼ 6 (e.g., McGreer et al. 2015). On the other hand, large fluctuations in IGM opacity remain at z ≲ 6, suggesting that lingering evidence of reionization may remain in the IGM to somewhat lower redshifts (Fan et al. 2006; Becker et al. 2015b; Bosman et al. 2018; Eilers et al. 2018). While current constraints from cosmic microwave background (CMB) observations are consistent with a rapid reionization at redshift z_rei ≃ 7.7 ± 0.7 (Planck Collaboration et al. 2018), measurements of the fraction of neutral hydrogen at high redshift have also been obtained from the presence of Lyα damping wings (Mortlock et al. 2011; Simcoe et al. 2012; Greig et al. 2017; Davies et al. 2018) and from the weakening of Lyα emission lines in z ∼ 6–8 galaxies (e.g., Caruana et al. 2014; Schmidt et al. 2016; Sadoun et al. 2017; Mason et al. 2018). The available data seem to generally support a late reionization scenario (with the bulk of reionization happening at z ∼ 6–8) but are still consistent with a relatively broad range of reionization histories.

The sources responsible for reionization also remain uncertain. Star-forming galaxies have commonly been considered the most likely candidate (e.g., Finkelstein 2016; Bouwens et al. 2015, 2016). Scenarios in which active galactic nuclei (AGNs) make a substantial contribution, however, continue to be considered (e.g., Giallongo et al. 2015; Madau & Haardt 2015; D'Aloisio et al. 2017; Parsa et al. 2018).

Further insight may be gained from using the IGM thermal history to constrain the reionization process. The temperature of the IGM should increase significantly via photoionization heating during hydrogen reionization, and because its cooling time is long, the low-density gas retains some useful memory of when and how it was reionized (Miralda-Escudé & Rees 1994; Abel & Haehnelt 1999; Upton Sanderbeck et al. 2016). At the mean density of the IGM the characteristic signature of reionization is expected to be an increase in temperature of tens of thousands of kelvin as an ionization front sweeps through (e.g., D'Aloisio et al. 2018), followed by cooling (over Δz ∼ 1–2) toward a thermal asymptote set primarily by the balance between photoheating by the UV background (UVB) and adiabatic cooling due to the expansion of the universe (e.g., McQuinn et al. 2009). The interplay among these effects is expected to lead to a power-law temperature–density (T–ρ) relation for the low-density gas (${\rm{\Delta }}=\rho /\bar{\rho }\lesssim 10$) of the form

$$\begin{eqnarray}&&T({\rm{\Delta }})={T}_{0}{{\rm{\Delta }}}^{\gamma -1},\end{eqnarray} \tag{ 1 }$$

where T₀ is the temperature at the mean density and γ – 1 is the slope of the relation (Hui & Gnedin 1997; Puchwein et al. 2015; McQuinn & Upton Sanderbeck 2016).

Following reionization, the increase in gas pressure due to the boost in temperature smooths out of the gas on small scales (e.g., Gnedin & Hui 1998; Rorai et al. 2013; Kulkarni et al. 2015). The degree of "Jeans smoothing" in the IGM prior to a given redshift is sensitive to timing and the total heat injection during and after reionization. Measurements of both the gas temperature evolution and the Jeans smoothing at redshifts approaching reionization (z ≳ 4) can therefore constrain the timing of this process and potentially provide information on the nature of the ionizing sources.

In the past two decades the Lyα forest in quasar spectra has been the main laboratory for the study of the thermal state of the IGM. In combination with cosmological hydrodynamical simulations, previous efforts have used a variety of statistical approaches to measure the IGM temperature–density relation at 1.5 ≲ z ≲ 5.4 from the shapes of the Lyα absorption lines (e.g., Ricotti et al. 2000; Schaye et al. 2000; McDonald et al. 2001; Theuns et al. 2002; Bolton et al. 2008, 2014; Lidz et al. 2010; Becker et al. 2011; Rudie et al. 2012; Boera et al. 2014, 2016; Hiss et al. 2018; Rorai et al. 2017a, 2018). However, the widths of these features are sensitive to both the instantaneous temperature of the gas (thermal broadening) and Jeans smoothing (which increases the Hubble broadening) due to the heat injection at higher redshifts. In previous works the impact of pressure smoothing either has generally not been included or has been considered a source of systematic error. For example, Viel et al. (2013a) and Iršič et al. (2017a) account for this effect by adding the redshift of reionization as a nuisance parameter for their warm dark matter constraints.

On the other hand, the first direct measurement of the characteristic filtering scale over which the gas is pressure smoothed (λ_P) has recently been obtained from the analysis of the Lyα absorption correlations using close quasar pairs at z ∼ 2–3 (Rorai et al. 2017b). This method largely disentangles the impacts of thermal broadening and pressure smoothing; however, the lack of known quasar pairs at higher redshifts prevents it from being used at redshifts closer to hydrogen reionization.

An alternative means of simultaneously constraining temperature and Jeans smoothing is presented by the Lyα flux power spectrum (Puchwein et al. 2015; Nasir et al. 2016, hereafter N16; Walther et al. 2019). N16 demonstrated, using hydrodynamical simulations, that the Lyα flux power spectrum exhibits different scale dependences for the temperature and Jeans smoothing. In particular, probing small scales (wavenumber log(k/km⁻¹ s) ≳ −1) increases the sensitivity to different reionization scenarios (see also Oñorbe et al. 2017 for an independent analysis). Although the 1D flux power spectrum statistic has been already explored in several works (e.g., Kim et al. 2004; McDonald et al. 2006; Viel et al. 2013a; Palanque-Delabrouille et al. 2015; Iršič et al. 2017a; Yèche et al. 2017), the lack of high-resolution, high signal-to-noise ratio (S/N) Lyα forest spectra has so far prevented these small scales from being measured at redshifts approaching reionization (but see Walther et al. 2018 for an analysis at z < 4).

In this paper we present a new measurement of the Lyα flux power spectrum at z ∼ 4−5.2 obtained from a sample of high-resolution, high-S/N spectra. We extend the measurement to previously unexplored small scales (log(k/km⁻¹ s) = −0.7). By comparing the data to predictions from a suite of hydrodynamical simulations, we investigate the IGM temperature evolution and, simultaneously, its integrated thermal history. We then demonstrate how the combined constraints offer new insights on the timing and sources of the hydrogen reionization process.

For this work we have adopted the parameterization of the Jeans smoothing effect described in N16. We characterize the integrated thermal history of the IGM using the cumulative energy per unit mass, u₀, injected into the gas at the mean cosmic density during and after the reionization process. As we demonstrate, this quantity can be directly used to constrain reionization models.

This paper is organized as follows. In Section 2 we introduce the observational sample of high-resolution spectra. An overview of the simulations used to test and interpret the measurements is presented in Section 3. In Section 4 we introduce the power spectrum method, discussing the effect of the most relevant thermal parameters. In Section 5 we present the observational power spectrum results and discuss the strategies applied to take into account and reduce systematic uncertainties. The calibration and analysis of the synthetic power spectrum models are described in Section 6. The Markov chain Monte Carlo (MCMC) analysis, comparing models with the observational measurements, is described in Section 7, where we also present our main results for the IGM temperature at the mean density and the integrated thermal history. As an example of how our thermal constraints can be used to test reionization histories, we apply our results to instantaneous reionization models in Section 8. We summarize our findings and conclude in Section 9. Tests for various systematic effects are described in the appendices.

We obtained high-resolution spectra of a sample of 15 quasars spanning emission redshifts 4.8 ≲ z_em ≲ 5.4. The quasars and their basic properties are listed in Table 1. The spectra for 11 of the objects were obtained with the Keck High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994), while the remaining four were taken with the Ultraviolet and Visual Echelle Spectrograph (UVES; Dekker et al. 2000) on the Very Large Telescope (VLT).

The spectra were reduced using a custom set of IDL routines that include optimal sky subtraction (Kelson 2003) and extraction techniques (Horne 1986). For each object a single 1D spectrum was extracted simultaneously from all exposures after individually applying telluric absorption corrections and relative flux calibration to the 2D frames. Telluric corrections were modeled based on the ESO SKYCALC Cerro Paranal Advanced Sky Model (Noll et al. 2012; Jones et al. 2013). For the UVES data we found that flux calibration derived from standard stars yielded sufficiently accurate agreement between overlapping orders. For HIRES, however, this approach produced well-known moderate (∼10%) interorder flux discrepancies. For all except one of our HIRES quasars, therefore, we used lower-resolution spectra from Keck/ESI, VLT/X-Shooter, or Gemini/GMOS to derive a custom response function for each exposure. The remaining object, J2111−0156, was calibrated using a response function from a standard star. We verified that our final flux power spectra remained essentially unchanged if standard-star flux calibration was used for every object. We therefore do not expect this aspect of the reduction to significantly impact our results.

The HIRES objects were observed using a 086 slit, giving a nominal resolution FWHM of ∼6 km s⁻¹. The UVES spectra were taken with a 10 slit, giving a nominal resolution of ∼7 km s⁻¹. The telluric models for the UVES data, however, indicated somewhat higher resolution consistent with a typical seeing of 08. Consequently, we adopt a resolution of 6 km s⁻¹ for the full data set, which is sufficient to resolve small-scale features in the Lyα forest. We therefore expect that even the smallest scale of the flux power spectrum measured in this work (log(k/km⁻¹ s) = −0.7, or Δv ∼ 30 km s⁻¹) will not be strongly affected by the finite spectroscopic resolution (but see Section 5.6.3). For all the quasars, the echelle orders were redispersed onto a common wavelength scale with a dispersion of 2.5 km s⁻¹ per pixel.

According to the analysis presented in N16 using mock observations with a redshift path Δz = 4, a continuum-to-noise ratio (C/N) of ∼15 per 3 km s⁻¹ pixel is necessary to break the degeneracy between thermal broadening and pressure smoothing and measure the thermal parameters with a statistical uncertainty of ∼20%. Conservatively, we have chosen our sample imposing this minimum threshold inside the Lyα forest region.

We have fitted the continuum in our spectra using spline fits guided by power-law extrapolations of the continuum redward of the Lyα emission line. Given the high levels of absorption at z ≳ 4, the continuum measurements are necessarily characterized by large uncertainties (∼10%–20%). We therefore use these estimations only to derive a rough estimate of the C/N level. In measuring the power spectrum, as described in Section 5.1, we adopt an approach that does not require a priori knowledge of the continuum. The redshift coverage and the median C/N for the Lyα forest region of our sample are reported in Table 1. The majority of the spectra have larger C/N than our cut, with a typical value per pixel in the range of 20–30. This high C/N assures that the power spectrum measurement at small scales will not be strongly affected by uncertainties in the noise modeling.

To test systematics associated with the observed power spectrum and to interpret our observational results, we used synthetic spectra derived from hydrodynamical simulations and processed to closely match the characteristics of the real data. We ran a large set of hydrodynamical simulations that span a range of thermal histories at z > 4. The simulations run following the Sherwood simulations suite (Bolton et al. 2017), which uses a modified version of the parallel smoothed particle hydrodynamics code P-GADGET-3, an updated and extended version of GADGET-2 (Springel 2005). The models adopt the cosmological parameters Ω_m = 0.308, Ω_Λ = 0.692, h = 0.678, Ω_b = 0.0482, σ₈ = 0.829, and n_s = 0.961, consistent with the CMB constraints of Planck Collaboration et al. (2014). Initial conditions were obtained using transfer functions generated by CAMB (Lewis et al. 2000). Because the vast majority of the absorption systems probed by the Lyα forest at z > 4 correspond to overdensities ${\rm{\Delta }}=\rho /\bar{\rho }\lesssim 10$, our analysis will not be affected by the star formation prescription (Viel et al. 2013b). Therefore, to increase the computational speed, gas particles with temperature T < 10⁵ K and overdensity Δ > 10³ are converted to collisionless particles (Viel et al. 2004).

The bulk of our simulations use a box size of 10 h⁻¹ cMpc and 2 × 512³ gas and dark matter particles, corresponding to a gas particle mass of 9.97 × 10⁴ h⁻¹ M_⊙. In addition, we use runs with larger box size and different mass resolution to test numerical convergence (see Appendix F).

We note that our simulations are not intended to self-consistently model reionization. Instead, we employ models with a wide variety of thermal histories so that our ultimate constraints on the temperature and integrated heating of the IGM are as general as possible. The gas in our models becomes optically thin at a redshift z_OT, after which it is photoionized and heated by a uniform UVB, which is a scaled version of the background from Haardt & Madau (2012). The thermal history of a given simulation is therefore determined by the choice of z_OT and UVB scaling factor.

The photoheating rates from Haardt & Madau (2012) (${\epsilon }_{i}^{\mathrm{HM}12}$) for the different species (i = [H i, He i, He ii]) have been rescaled proportionally by a constant factor ζ using the relation ${\epsilon }_{i}=\zeta {\epsilon }_{i}^{\mathrm{HM}12}$ (see Table 2). The combination of z_OT and ζ will determine both the instantaneous temperature and the total integrated heating per unit mass at the epoch where the power spectrum is measured. Models with larger z_OT and/or ζ will tend to have higher values of u₀.

Table 2.  List of Hydrodynamical Simulations Used in This Work

Model	L(h−1 cMpc)	Mgas(h−1M⊙)	zOT	ζ	${T}_{0}^{z=5}\,({\rm{K}})$	γz=5	${u}_{0}^{z=5}$($\mathrm{eV}\,{m}_{{\rm{p}}}^{-1}$)
S10_0.3z7	10	9.97 × 104	7	0.3	3162	1.52	1.3
S10_0.3z9	10	9.97 × 104	9	0.3	3388	1.49	2.3
S10_0.3z15	10	9.97 × 104	15	0.3	3388	1.51	5.0
S10_0.55z7	10	9.97 × 104	7	0.55	4553	1.51	1.9
S10_0.55z9	10	9.97 × 104	9	0.55	5086	1.48	3.3
S10_0.55z12	10	9.97 × 104	12	0.55	5110	1.51	5.2
S10_0.55z15	10	9.97 × 104	15	0.55	5093	1.52	7.4
S10_0.55z19	10	9.97 × 104	19	0.55	5074	1.52	10.0
S10_1z7	10	9.97 × 104	7	1.0	6607	1.50	2.7
S10_1z9	10	9.97 × 104	9	1.0	7413	1.51	4.7
S10_1z12	10	9.97 × 104	12	1.0	7510	1.51	7.6
S10_1z15 (10–512)	10	9.97 × 104	15	1.0	7413	1.50	10.6
S10_1z19	10	9.97 × 104	19	1.0	7457	1.52	14.7
S10_1.8z7	10	9.97 × 104	7	1.8	9725	1.49	3.9
S10_1.8z9	10	9.97 × 104	9	1.8	10866	1.50	6.8
S10_1.8z12	10	9.97 × 104	12	1.8	10900	1.51	10.9
S10_1.8z15	10	9.97 × 104	15	1.8	10865	1.51	15.5
S10_1.8z19	10	9.97 × 104	19	1.8	10827	1.52	21.4
S10_3.3z7	10	9.97 × 104	7	3.3	13803	1.48	5.5
S10_3.3z9	10	9.97 × 104	9	3.3	15488	1.50	9.9
S10_3.3z12	10	9.97 × 104	12	3.3	15821	1.48	16.3
S10_3.3z15	10	9.97 × 104	15	3.3	15488	1.52	23.2
S10_1z9_g1	10	9.97 × 104	9	1.0	7413	1.00	4.7
S20_1z15 (20–1024)	20	9.97 × 104	15	1.0	7413	1.50	10.6
S40_1z15 (40–2048)	40	9.97 × 104	15	1.0	7413	1.50	10.6
S40_1z9 (40–2048–zr9)	40	9.97 × 104	9	1.0	7413	1.51	4.7
S10_1z15_256	10	7.97 × 105	15	1.0	7413	1.50	10.6
S10_1z15_768	10	2.95 × 104	15	1.0	7413	1.50	10.6

Note. Entries in bold correspond to the Sherwood simulations first introduced in Bolton et al. (2017); the model names used in that work are given in brackets. For each simulation we report the name (Column (1)), box size (Column (2)), mass resolution (Column (3)), redshift at which the gas becomes optically thin (Column (4)), and constant factor used to rescale the photoheating rates for different thermal histories (Column (5)). The thermal parameters that describe the T–ρ relation at z = 5 are also listed: the temperature of the gas at the mean density (Column (6)) and the power-law index γ (Column (7)). Finally, the cumulative energy per unit mass deposited into the IGM at the mean density by z = 5 is given in Column (8) (see text for details). Further details on the simulation methodology are presented in Bolton et al. (2017).

Download table as:  ASCIITypeset image

A summary of the simulations used in this work is listed in Table 2. For each model we selected the simulation outputs between 4.0 ≲ z ≲ 5.4 with a redshift step Δz = 0.1. At each redshift, synthetic spectra of the Lyα forest were produced by choosing 5000 "lines of sight" throughout the simulation box. In Section 6.2 we describe how these lines of sight were combined to create realistic mock spectra.

Following N16, the integrated thermal history in our simulations is parameterized using u₀, the cumulative energy deposited per unit mass into the gas at the mean density. At each redshift u₀ is defined as

$$\begin{eqnarray}&&{u}_{0}(z)={\int }_{z}^{{z}_{\mathrm{OT}}}\displaystyle \frac{{\displaystyle \sum }_{i}{n}_{i}{\epsilon }_{i}}{\bar{\rho }}\displaystyle \frac{{dz}}{H(z)(1+z)},\end{eqnarray} \tag{ 2 }$$

where $\bar{\rho }$ is the mean mass density and n_i and _i represent, respectively, the number density and the photoheating rates for the species i = [H i, He i, He ii]. As shown in N16 (see their Figure 4), this parameter correlates with the density power spectrum of the cosmic gas in the simulations, with larger u₀ corresponding to a smoother distribution of gas for overdensities Δ < 10. These are the overdensities at which the Lyα forest is sensitive at z > 4 (e.g., Becker et al. 2011), suggesting that, at these redshifts, u₀ serves as a useful parameterization for the prior IGM thermal history. In Section 6.3.3 we will further consider the redshift range of integration over which u₀ optimally correlates with the flux power spectrum.

Examples of the evolution of u₀ in our models are presented in Figure 1, along with the corresponding evolution of the temperature at the mean density, T₀ (for the complete set of models, see Appendix I). The left panels show how increasing the photoheating rate in the simulations produces larger values in both the temperature and u₀. The right panels show models with the same photoheating rate but different z_OT. These converge to the same value of T₀ provided that sufficient time has elapsed after the onset of heating (Δz ∼ 1–2; e.g., McQuinn & Upton Sanderbeck 2016); however, they remain distinct in terms of u₀ values, reflecting differences in the total integrated thermal history and therefore in the amount of pressure smoothing.

Zoom In Zoom Out Reset image size

Both thermal broadening and pressure smoothing tend to reduce the amount of small-scale structure in the forest. Figure 2 shows the effect of thermal broadening (top panel) and pressure smoothing (bottom panel) on simulated Lyα forest spectra at z = 5. While the impacts of T₀ and u₀ are visually similar, the scale dependences of these effects make it possible to break the degeneracy (e.g., N16; Oñorbe et al. 2017).

Zoom In Zoom Out Reset image size

The top row of Figure 3 demonstrates how the shape of the 1D Lyα flux power spectrum at z = 5 varies for models with different instantaneous temperature (left panel) and integrated thermal histories (right panel). Similar results were shown in N16 but are expanded here to include a broader range of thermal histories. As described in Section 6.3.1, we use post-processing to vary the T₀ for a fixed u₀ (top left) or to impose the same T₀ for models with different u₀ (top right). We also demonstrate the impact of varying γ (bottom left) and the effective optical depth, τ_eff (bottom right). As noted by N16, the scale dependences of T₀ and u₀ differ somewhat. While the impact of pure thermal broadening increases continuously toward smaller scales, the effect of changing u₀ peaks near log(k/km⁻¹ s) ∼ −0.9 to −0.8.

Zoom In Zoom Out Reset image size

Comparing the two panels of the first row of Figure 3, it is clear that in order to distinguish models characterized by an early reionization (large u₀) from those with high T₀ values, it is necessary to probe the power spectrum down to log(k/km⁻¹ s) ∼ −0.7. Our effort to measure the power spectrum down to these scales is described in the following section.

In this section we describe our procedure for measuring the flux power spectrum from the observed spectra. The following strategies have been tested using synthetic spectra for two reasons: first, to detect and quantify systematic effects in the calculation of the power spectrum, and second, to guarantee a fair comparison between simulated models and the observed data.

5.1. Rolling Mean

We performed the power spectrum measurement on the flux contrast estimator

$$\begin{eqnarray}&&{\delta }_{F}=\displaystyle \frac{F-\mathop{F}\limits^{\mbox{--}}}{\mathop{F}\limits^{\mbox{--}}},\end{eqnarray} \tag{ 3 }$$

where F is the transmission in the Lyα forest and $\mathop{F}\limits^{\mbox{--}}$ is the mean flux. When computing δ_F, we need to first divide out the intrinsic shape of the quasar spectrum, which can impact the power spectrum at large scales (log(k/km⁻¹ s) ≲ −2; e.g., Kim et al. 2004; Viel et al. 2013a; Iršič et al. 2017a). However, directly estimating the continuum is difficult at z ≳ 4 owing to the high levels of absorption in the Lyα forest. We therefore used a rolling mean approach, wherein $\bar{F}$ is estimated locally by smoothing the observed spectrum using a boxcar average. We used a boxcar window of 40 h⁻¹ cMpc, which was chosen to be short enough to roughly capture relevant features in the quasar continua over the forest (see Appendix A.1 for details). Examples of this approach are presented in Figure 4.

Zoom In Zoom Out Reset image size

5.2. Proximity Regions

5.3. DLAs

5.4. Bad Pixels

5.5. Lyα Sections and Redshift Subsamples

5.6. Measuring the Power Spectrum

For each of the 20 h⁻¹ cMpc forest regions we calculate the power spectrum from the flux contrast δ_F defined in Equation (3). Our spectra are unevenly sampled because they are masked, so we use a Lomb–Scargle algorithm (Lomb 1976; Scargle 1982) to compute the raw power of each region (P_masked(k)). In all of our calculations we use k-bins logarithmically spaced with Δlog k = 0.1. To obtain the final power spectrum values, P_F(k), for each section, we first correct the raw P_masked(k) for the effect of masking. Second, we subtract from the corrected P_data(k) an estimate of the contribution to the power from noise, P_N(k). All these steps are described in the following sections.

5.6.1. Masking Correction Function

The masking procedure described in Sections 5.3 and 5.4, and in particular the masking of sky line residuals, impacts the power spectrum owing to the application of a complex window function. In order to correct for this, we apply a masking correction function, C_m(k), to the raw power obtained from each of the 20 h⁻¹ cMpc Lyα forest sections,

$$\begin{eqnarray}&&{P}_{\mathrm{data}}(k)={P}_{\mathrm{masked}}(k)\times {C}_{m}(k),\end{eqnarray} \tag{ 4 }$$

where P_data(k) is the corrected quantity used to infer the final power and P_masked(k) is the raw power initially computed from masked spectra.

We determine the effect of masking for each of the Lyα forest sections contributing to the analysis using the following procedure. First, we create hundreds of synthetic spectra with the same characteristics (i.e., size, noise, redshift) of each of the real 20 h⁻¹ cMpc sections with and without the same masking applied. The final correction is then obtained from the average of the ratio between the power of the unmasked (P_sim) and masked (${P}_{\mathrm{sim}}^{\mathrm{mask}}$) simulated spectra,

$$\begin{eqnarray}&&{C}_{m}(k)=\left\langle \displaystyle \frac{{P}_{\mathrm{sim}}(k)}{{P}_{\mathrm{sim}}^{\mathrm{mask}}(k)}\right\rangle .\end{eqnarray} \tag{ 5 }$$

Because the impact of masking on the power spectrum in principle depends on its underlying shape, possible systematics may arise from choosing a particular simulation run to compute C_m(k). We use the 20 h⁻¹ cMpc run of Table 2 for the final correction; therefore, we quantify these possible uncertainties in Appendix A.5.

5.6.2. Noise Subtraction

In principle, the noise power can be directly computed from the flux error array output by the data reduction pipeline (e.g., Iršič et al. 2017a; Walther et al. 2018). This approach, however, relies on the precision of the pipeline; underestimating or overestimating these uncertainties could significantly impact the final power, especially at the small scales we are interested in. We therefore estimate the amount of noise for each forest section directly from the raw power spectrum of the data. At the smallest scales (log(k/km⁻¹ s) ≳ −0.2) the power is dominated by noise fluctuations and, assuming that the noise in adjacent wavelength bins is uncorrelated, can be fitted with a constant value. We then assume that P_N is constant over all scales and subtract it from the total power obtaining the "noiseless" power spectra. This method is illustrated in Figure 5 and has been tested on synthetic data after adding the observational noise arrays to the simulated lines of sight (see Appendix A.3).

Zoom In Zoom Out Reset image size

5.6.3. Resolution Correction

In this work we forward-model the synthetic spectra generated from simulations to match the instrumental resolution and pixel size of the data. For reference, however, we include a version of the observed flux power spectrum that has been corrected for resolution as

$$\begin{eqnarray}&&{P}_{{\rm{F}}}(k)=\displaystyle \frac{{P}_{\mathrm{data}}(k)-{P}_{{\rm{N}}}(k)}{{W}_{R}^{2}(k,R,{{d}v}_{p})},\end{eqnarray} \tag{ 6 }$$

using the window function adopted in Palanque-Delabrouille et al. (2015),

$$\begin{eqnarray}&&{W}_{R}(k,R,{{d}v}_{p})=\exp \left(-\displaystyle \frac{1}{2}{\left({kR}\right)}^{2}\right)\displaystyle \frac{\sin ({{kd}v}_{p}/2)}{({{kd}v}_{p}/2)}.\end{eqnarray} \tag{ 7 }$$

Assuming our nominal resolution R = 2.55 km s⁻¹ (FWHM = 6 km s⁻¹) and pixel size ${dv}_{p}$ = 2.5 km s⁻¹, the correction for the smallest scale considered in this work (log(k/km⁻¹ s) = −0.7) is ${W}_{R}^{2}\sim 0.76$.

We note that the actual spectral resolution of the data will depend on the seeing of the observation and may be different from the nominal one. A possible error in the power spectrum due to uncertainties in the spectral resolution, even when forward-modeling the simulations, must therefore be taken into account. We estimate an error of 10% in the spectral resolution, corresponding to an uncertainty in the power of ≲5% at log(k/km⁻¹ s) ≲ −0.7. This correction is smaller than our statistical error, so we do not expect that uncertainties in the resolution will significantly affect the measurements (see Appendix A.4).

5.7. Metals

The flux power spectrum measured directly from the observational spectra contains both the power coming from the Lyα forest and a small contribution from intervening metal lines. These lines tend to show individual components significantly narrower than Lyα (b ≲ 15 km s⁻¹), which will increase the power on small scales (e.g., Lidz et al. 2010). This section describes our approach to quantifying and removing the effect of metals on our final power spectrum measurements.

The high level of Lyα absorption at high redshift makes it very challenging to directly identify all metal lines in the forest. We therefore estimate the metal power spectrum directly from regions of quasar spectra redward of the Lyα emission line, where only metal absorption systems are present (e.g., McDonald et al. 2005; Palanque-Delabrouille et al. 2013). The metal power measured in this way will not take into account transitions with rest-frame wavelength shorter than the Lyα line. Correlation features like the one observed for Si iii (λ1206) in McDonald et al. (2006), however, will tend to affect the power spectrum on scales larger than the ones considered in this work (log(k/km⁻¹ s) ≲ −2.5).

We measured the metal power spectrum from two samples of high-resolution quasar spectra. First, we use a subset of the spectra listed in Table 1 with emission redshift 4.5 ≲ z_em ≲ 5.3. Second, we use a sample of spectra of quasars with emissions redshifts 3.4 ≲ z_em ≲ 4.1 from Boera et al. (2014, 2016) (Table 3). The latter sample allows us to measure the metal power spectrum over observed wavelengths similar to those spanned by the Lyα forest at 4.0 ≲ z ≲ 4.4. While metals redward of Lyα are not a perfect estimate of those that appear in the forest at higher redshifts, analyzing multiple samples allows us to check for redshift evolution in the metal power spectrum.

Table 3.  List of Quasars Used for the Analysis of the Metal Power Spectrum at z ∼ 4.2

Name	zem
J010604−254651	3.36500
J162116−004250	3.70270
J132029−052335	3.70000
J124957−015928	3.63680
J014049−083942	3.71290
J115538+053050	3.47520
J014214+002324	3.37140
J123055−113909	3.52800
J110855+120953	3.67160
J005758−264314	3.65500

Note. For each object we report the name (Column (1)) based on the J2000 coordinates of the quasar and the emission redshift (Column (2)). All the spectra have been taken with the UVES spectrograph (see Boera et al. 2014 and Murphy et al. 2019 for details).

Download table as:  ASCIITypeset image

Figure 6 shows the comparison between our measurements of the metal power spectrum at z ∼ 4.2 (blue dashed line) and at z ∼ 5.4 (red dashed line), obtained using the same data analysis procedure described in the previous sections. The most significant difference in the metal power between these two redshifts is at large scales (log(k/km⁻¹ s) ≲ −1.5), where the contribution of metals to the final flux power spectrum measurement is less relevant. Note that our measurement of the metal power spectrum at z ∼ 4.2 is also in good agreement with the one computed from the XQ-100 data sample at the same redshift (Iršič et al. 2017a) (blue dotted line) even if the latter is slightly noisier.

Zoom In Zoom Out Reset image size

Given the weak evolution in redshift of the metal power (already noted in previous works; e.g., Palanque-Delabrouille et al. 2013), we corrected our final power spectrum measurements assuming a metal contribution constant with redshift and equal to the average between our two measured metal power spectra (black solid line in Figure 6). The effect of the metal correction on the final flux power spectrum is shown in Appendix A.6.

5.8. The New Power Spectrum Measurements

The main observational results of this work are presented in Figure 7, where we plot the final Lyα flux power spectrum measured from our data. The values are tabulated in Appendix K and reported as a function of scale for the three redshift bins centered at z = 4.2, 4.6, and 5.0. Solid colored lines and data points represent the power spectrum results without the resolution correction described in Section 5.6.3, while the corresponding dashed lines are the measurements corrected for finite resolution and pixel size. The 1σ errors are estimated from the bootstrap covariance matrix of the data, corrected and regularized following the procedure described in Section 5.9.

Zoom In Zoom Out Reset image size

Our measurements at scales log(k/km⁻¹ s) ≳ −1 are the first ones made at these redshifts (see Walther et al. 2018 for an analysis at z < 4). At larger scales, however, we can compare with the results derived from high-resolution spectra by Viel et al. (2013a) and medium-resolution data by Iršič et al. (2017a). These comparisons are presented in Appendices C and D.

5.9. Covariance Matrix

As demonstrated in previous works (e.g., Viel et al. 2013a; Iršič et al. 2017a), the covariance matrix obtained via bootstrapping of a limited data set is necessarily noisy. We therefore regularized the observed covariance matrix using the correlation coefficients estimated from the simulated spectra following an approach similar to the one used by Lidz et al. (2006). We first used the simulations to verify the ability of the bootstrapped errors to reproduce the real statistical variance. For this test, using the 40 h⁻¹ cMpc box simulation S40_1z15, we created hundreds of samples of simulated lines of sight that closely reproduce the characteristics of the observational data (see Section 6.2 for details), and we compare the variance computed directly from these realizations with the uncertainty obtained from the bootstrapping of only one synthetic sample randomly chosen. We verify that, as already shown by previous studies (e.g., Kim et al. 2004; Iršič et al. 2013; Palanque-Delabrouille et al. 2013), the bootstrapping technique underestimates the cosmic variance by up to ∼25%, with the discrepancy level increasing toward smaller scales. We therefore increased the observational bootstrapped error at all scales by ∼15%–25%, where the correction has been computed separately for each redshift.

The elements of the final covariance matrix, C_ij, are then computed as

$$\begin{eqnarray}&&{C}_{{ij}}={R}_{{ij}}\sqrt{{C}_{{ii}}^{\mathrm{data}}{C}_{{jj}}^{\mathrm{data}}},\end{eqnarray} \tag{ 8 }$$

with

$$\begin{eqnarray}&&{R}_{{ij}}=\displaystyle \frac{{C}_{{ij}}^{\mathrm{sim}}}{\sqrt{{C}_{{ii}}^{\mathrm{sim}}{C}_{{jj}}^{\mathrm{sim}}}},\end{eqnarray} \tag{ 9 }$$

where ${C}_{{ii}}^{\mathrm{data}}$ are the diagonal elements of the bootstrapped observational covariance matrix, corrected as previously described, and C^sim are the elements of the simulated covariance matrix obtained from the multiple realizations of synthetic lines of sight.

For comparison, Figure 8 shows the simulated and observational covariance matrices for the redshift bin at z = 5. As expected, the bootstrapped matrix is noisier, but both the matrices show a similar structure, with correlations increasing toward the smallest scales, log(k/km⁻¹ s) ≳ −1.5. This similarity gives us confidence that the simulation model used for the covariance matrix regularization is reasonably capturing the data properties. We note that the off-diagonal correlation structure will depend somewhat on the precise shape of the power spectrum and therefore on the thermal parameters characterizing the model. In Appendix B we verify that the particular choice of the simulation S40_1z15 for this analysis is not significantly affecting our final constraints.

Zoom In Zoom Out Reset image size

In this section we describe how we calibrate and analyze synthetic spectra to create power spectrum models that will be used to fit the observational measurements in Section 7.

6.1. Constructing Mock Lines of Sight

To ensure the correct comparison between simulations and observational data, we need to produce mock lines of sight with the same resolution and redshift coverage as our observed sample. We first resample and smooth the synthetic Lyα spectra produced from the simulations in Table 2 to match the spectral resolution and the pixel size of the real spectra. We then progressively merge multiple synthetic sections randomly selected from the Δz = 0.1 simulation snapshot closest to the Lyα forest redshift that we want to cover. We choose an arbitrary starting point along each section, taking advantage of the periodicity of the simulation box. We take into account the mild redshift evolution of the mean flux along the line of sight by rescaling the optical depths such that the global effective Lyα optical depth (${\tau }_{\mathrm{eff}}=-\mathrm{ln}(\bar{F})$) in the simulation box follows the relation

$$\begin{eqnarray}&&{\tau }_{\mathrm{eff}}=1.56{\left(\displaystyle \frac{1+z}{5.75}\right)}^{4.0}.\end{eqnarray} \tag{ 10 }$$

See Appendix G for details. We note that while we use Equation (10) to calibrate the optical depth evolution within each simulated line of sight, the overall Lyα mean flux of each redshift bin will be treated as a free parameter in our models.

When required for testing purposes (see Appendix A), the spectral noise at the same level of the corresponding observational line of sight is added to the synthetic spectra, as well as the same pixel masking. Because the power spectrum computed from the real data is corrected for these systematics, the final power spectrum models used for the MCMC fit are computed without noise or pixel masking.

6.2. Modeling the Power Spectrum

6.3. Varying Model Parameters

To be able to fit the power spectrum measurement of Section 5.8, we need a grid of models that cover the parameter space that we want to explore. In each redshift bin we consider four parameters to describe the power spectrum: the thermal parameters T₀, u₀, and γ and the effective Lyα optical depth, τ_eff. While the large set of simulations listed in Table 2 spans a wide range of thermal histories, by themselves they are not sufficient to evaluate the power spectrum in all the possible combinations of thermal parameters. We therefore use the interpolation scheme described below.

6.3.1. Varying T₀ and $\gamma$

In order to separate the impact of thermal broadening and Jeans smoothing in the power spectrum models, we applied a simple post-processing procedure to the simulated spectra. This is achieved by translating and rotating the entire T–ρ plane of the simulations to match the new T₀ and γ values. We recompute the optical depths in each of our models over an extended range of power-law T–ρ relationships, with T₀ = 3000–15,000 K in steps of 1000 K and γ = 0.7–1.7 in steps of 0.1.

We note that at z ≳ 4 the Lyα forest is mainly sensitive to gas close to the mean density (e.g., Becker et al. 2011). For this reason we do not expect to place strong constraints on γ. We nevertheless treat γ as a free parameter in our fitting code. The impacts of T₀ and γ on the flux power spectrum are demonstrated in Figure 3. Note that scales log(k/km⁻¹ s) ≳ −0.8 seem to be insensitive to variations in γ, while considerable changes in this parameter create minor shifts in the power for scales log(k/km⁻¹ s) ≲ −0.8.

6.3.2. Varying τ_eff

We rescaled the optical depths in our models to span a wide range of τ_eff values. At each redshift the reference value has been obtained from Equation (10), while the entire range of optical depths covered by our models is τ_eff = 0.6–2.2 in steps of 0.1. Different mean fluxes within each redshift bin have been obtained by multiplying Equation (10) by single, fine-tuned scalars when calibrating the simulated lines of sight. The impact of varying τ_eff on the power spectrum is shown in the bottom right panel of Figure 3.

6.3.3. Varying u₀

By post-processing our simulations to a common set of thermal parameters, we can isolate how the power spectrum depends on the integrated heating. N16 demonstrated that the flux power spectrum at z = 5.0 (averaged over scales −1.5 ≲ log(k/km⁻¹ s) ≲ −0.8) correlates with u₀. They further argued that the correlation is strongest when u₀ is integrated between z = 12 and 5, reflecting the timescales over which the Jeans smoothing is sensitive to heat injection. Here we reevaluate this redshift dependence using our more extended suite of models. For each of the redshift bins at which we compute the power spectrum and each of the scales sensitive to u₀ (−1.4 ≲ log(k/km⁻¹ s) ≲ −0.8) we empirically determine the "characteristic" redshift range of integration ($\bar{{\rm{\Delta }}{z}^{{u}_{0}}}$) for which the power is closest to a one-to-one function of u₀. The method is demonstrated in Figure 9. We first post-process all of the 10 h⁻¹ cMpc simulations of Table 2 to the same values of T₀, γ, and τ_eff. We then fit a power law to kP_k versus u₀, where u₀ is integrated using Equation (2) over a redshift interval ${\rm{\Delta }}{z}^{{u}_{0}}$. The preferred interval, $\bar{{\rm{\Delta }}{z}^{{u}_{0}}}$, is the one that minimizes the χ² for this fit.

Zoom In Zoom Out Reset image size

The characteristic $\bar{{\rm{\Delta }}{z}^{{u}_{0}}}$ computed for the scales −1.4 ≲ log(k/km⁻¹ s) ≲ −0.8 are reported in Figure 10 for the different redshift bins. As expected, the sensitivity of the power spectrum to the previous thermal history varies slightly with the redshift at which the power spectrum is measured. The Lyα structures observed at progressively lower redshifts seem to slowly lose sensitivity to earlier epochs; while the power spectrum measured at z = 5.0 still maintains sensitivity up to z ≳ 13, at z = 4.2 the forest traces the thermal history of the gas mainly for z ≲ 12. Interestingly, we find that the power spectrum at z = 5.0 is less sensitive to heating happening at z ≲ 6, even though the power spectra at z = 4.6 and 4.2 retain sensitivity all the way down to their respective redshifts. We generally expect that the gas density distribution will exhibit some delay in responding to changes in gas pressure. Further investigation revealed that a delay did appear for all three redshifts when peculiar velocities were turned off. The delay increased with the redshift at which the power spectrum was measured, with a delay at z = 5.0 that was larger than the one found with peculiar velocities turned on. This suggests that peculiar velocities may play a role by decreasing the delay between heat injection and a change in the power spectrum. Presumably this occurs because, as the gas is heated, redshift distortions created by accelerating the gas precede changes in the density field. This effect may partly explain the lack of a gap at z = 4.2 and 4.6. For now we adopt these relations as empirical and leave more detailed physical insights to future work.

Zoom In Zoom Out Reset image size

Because the $\bar{{\rm{\Delta }}{z}^{{u}_{0}}}$ values are generally constant among different scales within the same redshift bin, we adopt average values (black dashed vertical lines in Figure 10) as integration bounds in Equation (2). The fiducial redshift range over which u₀ is integrated is given in Table 4. We note that our u₀–kP_k relationship, while remarkably tight over scales sensitive to u₀, does exhibit scatter. In the final MCMC analysis, therefore, the amount of scatter about the u₀–kP_k fit at each scale has been included as systematic uncertainty. We note that while we chose our fiducial redshift ranges to maximize the sensitivity of the power spectrum to u₀, in principle we could constrain this parameter integrated within any reasonable redshift range if properly accounting for the systematic uncertainty in the u₀ versus kP_k fit.

Table 4.  Fiducial Redshift Ranges Used to Compute the u₀ Parameters for Fitting the Flux Power Spectrum at Different Redshifts

z	${\rm{\Delta }}{z}_{\mathrm{fid}}^{{u}_{0}}$
4.2	4.2–12
4.6	4.6–13
5.0	6.0–13

Note. Redshift bins are indicated in Column (1), with the corresponding fiducial u₀ redshift intervals in Column (2).

Download table as:  ASCIITypeset image

To obtain constraints on the IGM temperature and integrated thermal history from the observational power spectrum measurements obtained in Section 5, we adopted a Bayesian MCMC approach to measure T₀, u₀, γ, and τ_eff for each of the three redshift bins independently. In this section we present the method and the main findings of this analysis.

7.1. The MCMC Analysis

We constructed a grid of power spectrum models following the post-processing approach given above, where for a given choice of T₀, γ, and τ_eff the dependence of the power spectrum on u₀ is derived from the fits described in Section 6.3.3. We then perform a multilinear interpolation among the grid points of the four-dimensional parameter space. We implemented the interpolation scheme using a Bayesian MCMC approach. At each redshift, applying flat priors for all variables, we obtain the set of parameters that maximize a Gaussian multivariate likelihood function:

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=-\displaystyle \frac{1}{2}{{\boldsymbol{\Delta }}}^{T}{C}^{-1}{\boldsymbol{\Delta }}-\displaystyle \frac{1}{2}\mathrm{ln}\,\det (C)-\displaystyle \frac{N}{2}\mathrm{ln}2\pi ,\end{eqnarray} \tag{ 11 }$$

where ${\boldsymbol{\Delta }}$ is the residual vector between the power spectrum values of the data and the model and C is the N × N data covariance matrix (where N is the number of data points).

We tested the interpolation scheme by removing one of the models used for the interpolation and using it to generate mock data (see Appendix H). We found that the parameters u₀ and T₀ are recovered accurately. Small biases (within the 68% uncertainties) appear in the recovered values of γ and τ_eff owing to their intrinsic degeneracy at large scales and the poor sensitivity of the high-redshift power spectrum to γ. Fortunately, however, the relatively weak constraints on these parameters do not bias our results for T₀ and u₀.

To test the reliability of the best-fitting values, for each redshift we ran three independent chains of 2 × 10⁵ iterations (half of which are discarded as burn-in) from different randomly chosen initial parameters. We verify that all the chains were converged by comparing the between-chain and within-chain variances for each parameter using the Gelman–Rubin test.

7.2. Results

Figures 11–13 display the posterior likelihood distributions for the parameters T₀, u₀, γ, and τ_eff at redshifts 4.2, 4.6, and 5.0, respectively. While the inclusion of small scales (log(k/km⁻¹ s) ≳ −1.0) in the power spectrum allows relatively tight constraints on both T₀ and u₀, some degeneracy between these two variables is still noticeable at all redshifts. This is expected since both of these parameters act on intermediate scales in a similar way. Degeneracies between γ and τ_eff increase with redshift, with slightly weaker constraints on τ_eff obtained toward higher redshifts. As expected, γ shows broad bounds at all redshift (with the 1σ contours covering almost the entire parameter space), reaffirming that the Lyα forest at high redshifts mainly probes gas around the mean density and is not highly sensitive to the slope of the T–ρ relation. Figures 11–13 demonstrate that our measurements of u₀ and T₀ are not highly affected by degeneracies with γ and τ_eff.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The final results of the MCMC analysis are summarized in Table 5. The temperatures are constrained with ∼15% uncertainties at all redshifts, while the error on u₀ varies from ∼18% for the z = 4.2 and z = 4.6 redshift bins up to ∼30% at the highest redshift, in good agreement with the forecast presented by N16. Our results for τ_eff are highly consistent with the measurements of Becker et al. (2013) at z = 4.2. We are somewhat higher at z = 4.6, but altogether our constraints appear to bridge the evolution of τ_eff at z ≲ 4 measured by Becker et al. (2013) and at z ≳ 5 from Bosman et al. (2018) (see Appendix G).

Table 5.  Best-fitting Values and Marginalized 68% Confidence Intervals for the Fits to Our Power Spectrum Measurements

z	T0/103 (K)	u0(eV ${m}_{p}^{-1})$	γ	τeff
4.2	${8.13}_{-0.97}^{+1.34}$	${7.29}_{-1.35}^{+0.98}$	${1.21}_{-0.28}^{+0.23}$	${1.02}_{-0.04}^{+0.04}$
4.6	${7.31}_{-0.88}^{+1.35}$	${7.10}_{-1.45}^{+0.83}$	${1.29}_{-0.26}^{+0.19}$	${1.41}_{-0.09}^{+0.08}$
5.0	${7.37}_{-1.39}^{+1.67}$	${4.57}_{-1.16}^{+1.37}$	${1.33}_{-0.27}^{+0.18}$	${1.69}_{-0.11}^{+0.10}$

Note. The power spectrum redshift (Column (1)) is reported along with the best-fitting values of T₀ (Column (2)), u₀ (Column (3)), γ (Column (4)), and τ_eff (Column (5)).

Download table as:  ASCIITypeset image

In Figure 14 we show the best power spectrum models compared with the measurement. Visually there is good agreement with the data at all redshifts.

Zoom In Zoom Out Reset image size

The final temperature constraints at the IGM mean density are presented in Figure 15. Our new results (green points) are compared with the previous measurements of Becker et al. (2011) (gray triangles) at z > 3.5 obtained with the curvature method. We have added the systematic uncertainty for Jeans smoothing estimated by Becker et al. to those data. See also Appendix E for a more comprehensive comparison of recent temperature measurements at the redshifts covered by our analysis. Our temperature measurements show good agreements with this previous work in the overlapping redshift bins. This accord is significant because we analyzed a largely independent set of quasar spectra with a different method and using a new suite of hydrodynamical simulations. Most significantly, we now explicitly fit for u₀, removing the systematic uncertainty in T₀ related to Jeans smoothing. Overall our T₀ values are consistent with little evolution over 4.2 ≲ z ≲ 5.0. Given the known trend of increasing temperatures at z ≲ 4 (e.g., Becker et al. 2011; Boera et al. 2014; Walther et al. 2019), our measurement at z = 4.2 may include a contribution from the initial phase of IGM reheating due to the He ii reionization (e.g., Worseck et al. 2011; Syphers & Shull 2014). We consider this possibility in the final part of our analysis, when we use the new thermal constraints to evaluate hydrogen reionization scenarios.

Zoom In Zoom Out Reset image size

Finally, in Figure 16 we present our first constraints on the integrated thermal history of the IGM. Our u₀ measurements are plotted at the minimum redshift of the fiducial ranges given in Table 5. As expected, u₀ increases from z = 6 to z = 4.2, reflecting ongoing heat injection after reionization.

Zoom In Zoom Out Reset image size

Our observational constraints on T₀ and u₀ can, in principle, be used to test any reionization model for which a thermal history can be calculated. While we leave for future work the analysis of extended and more realistic reionization scenarios, in this section we demonstrate the potential of this approach using semianalytical models of instantaneous reionization.

8.1. Modeling Instantaneous Reionization

We model the thermal history of instantaneous hydrogen reionization using a semianalytical approach similar to the one adopted in Upton Sanderbeck et al. (2016). To obtain T₀ as a function of redshift, we solve the equation describing the temperature evolution of a Lagrangian fluid element at the cosmic mean density, i.e., with Δ = 1 (e.g., Miralda-Escudé & Rees 1994; Hui & Gnedin 1997; McQuinn & Upton Sanderbeck 2016),

$$\begin{eqnarray}&&\displaystyle \frac{{dT}}{{dt}}=-2{HT}+\displaystyle \frac{2T}{3{\rm{\Delta }}}\displaystyle \frac{d{\rm{\Delta }}}{{dt}}+\displaystyle \frac{2}{3{k}_{B}{n}_{\mathrm{tot}}}\displaystyle \frac{{dQ}}{{dt}},\end{eqnarray} \tag{ 12 }$$

where H is the Hubble parameter and n_tot is the total number density of particles (electrons and ions). Equation (12) is valid in the approximation that the number of particles remains fixed, describing well the post-reionization gas. The first term on the right-hand side of Equation (12) takes into account the cooling due to adiabatic expansion, while the second term gives the adiabatic heating and cooling due to structure formation. The thermal history at the mean density has been shown to depend weakly on this second term (McQuinn & Upton Sanderbeck 2016); we will therefore ignore it in our calculation. The differences among the models will depend instead on the third term, which encodes photoheating (the only heating source considered in these calculations) and radiative cooling processes. We can expand this term as

$$\begin{eqnarray}&&\displaystyle \frac{{dQ}}{{dt}}=\displaystyle \sum _{X}\displaystyle \frac{{{dQ}}_{\mathrm{photo},X}}{{dt}}+\displaystyle \frac{{{dQ}}_{\mathrm{Compton}}}{{dt}}+\displaystyle \sum _{i}\displaystyle \sum _{X}{R}_{i,X}{n}_{e}{n}_{X},\end{eqnarray} \tag{ 13 }$$

where ${{dQ}}_{\mathrm{photo},X}/{dt}$ is the photoheating rate of ion X, ${{dQ}}_{\mathrm{Compton}}{dt}$ is the Compton cooling rate, and R_i,X is the cooling rate coefficient for the ion X and cooling mechanism i. Because we are modeling the temperature of the gas at the end of hydrogen reionization, in Equation (13) we will consider only the H i (dominant) and He i photoheating. As for the cooling term, we include Compton and H ii recombination in our calculations. As discussed in McQuinn & Upton Sanderbeck (2016), these represent the relevant cooling processes that shape the temperature evolution. We compute these cooling terms using the rate coefficients provided by Hui & Gnedin (1997). The optically thin photoheating after reionization is modeled as (e.g., Upton Sanderbeck et al. 2016)

$$\begin{eqnarray}&&\displaystyle \frac{{{dQ}}_{\mathrm{photo},X}}{{dt}}\approx \displaystyle \frac{h{\nu }_{X}}{{\gamma }_{X}-1+{\alpha }_{{bk}}}{\alpha }_{A,X}{n}_{\tilde{X}}{n}_{e},\end{eqnarray} \tag{ 14 }$$

where ν_X is the frequency associated with the ionization potential of species X, and γ_X is the corresponding approximate power-law index of the photoionization cross section, for which we assume γ_X = 2.8 for H i and γ_X = 1.7 for He i. Equation (14) is valid in the approximation of photoionization equilibrium with an ionizing background that has a power-law specific intensity of the form ${J}_{\nu }\propto {\nu }^{-{\alpha }_{\mathrm{bk}}}$. In detail, the photoheating rate will also depend on α_A,X, the case A recombination coefficient associated with the transition from X i → X for species X ∈ [H i, He i]; on the number density of the species $\tilde{X}\in$ [H, He]; and on the electron number density n_e.

To compute the total energy deposited into the gas by photoheating, u₀, we just need to consider the third term of Equation (12) (and the first term of Equation (13)). The equation to solve for u₀ will then be

$$\begin{eqnarray}&&\displaystyle \frac{{dT}}{{dt}}=+\displaystyle \frac{2}{3{k}_{B}{n}_{\mathrm{tot}}}\displaystyle \sum _{X}\displaystyle \frac{{{dQ}}_{\mathrm{photo},X}}{{dt}}.\end{eqnarray} \tag{ 15 }$$

Because the specific internal energy can be expressed as $u=\tfrac{3}{2}{k}_{B}T\tfrac{1}{m}$, Equation (15) can be solved as

$$\begin{eqnarray}&&\displaystyle \frac{{{du}}_{0}}{{dt}}=\displaystyle \frac{1}{\bar{\rho }}\displaystyle \sum _{X}\displaystyle \frac{{{dQ}}_{\mathrm{photo},X}}{{dt}},\end{eqnarray} \tag{ 16 }$$

where $\bar{\rho }$ is the mean mass density.

8.2. Instantaneous Reionization Parameters

We parameterize our models using three numbers: the redshift of instantaneous reionization, z_rei, the temperature reached by the IGM during hydrogen reionization, T_rei, and the spectral index of the post-reionization ionizing background, α_bk. Figure 17 presents the effects on the evolution of T₀ (top row) and u₀ (bottom row) of these three parameters. The first column shows models with the same z_rei and α_bk but different reionization temperatures. Radiative transfer calculations suggest that temperatures during reionization should reach 17,000 K ≲ T_rei ≲ 25,000 K (e.g., Miralda-Escudé & Rees 1994; D'Aloisio et al. 2018); however, we explored T_rei down to 10,000 K and up to 30,000 K, where the upper range is similar to the hottest scenario of short and late reionization considered in D'Aloisio et al. (2018). The second column in Figure 17 demonstrates how changing the timing of reionization influences the histories of T₀ and u₀. We tested models spanning a range of redshifts from z_rei = 5.5 up to z_rei = 12. Finally, the third column shows the effect of changing the spectral index of the post-reionization ionizing background. The value of α_bk can be connected to the intrinsic spectral index of the sources, α_s, via the expression α_bk ≈ α_s − 3(β − 1), where β is the logarithmic slope of the column density distribution of intergalactic hydrogen absorbers (Upton Sanderbeck et al. 2016), which is valid at z ≳ 3 when the physical mean free path of 1 ryd photons λ_MFP ≪ cH⁻¹. The value of β may vary, but for this analysis we adopt β = 1.3 from Songaila & Cowie (2010).

Zoom In Zoom Out Reset image size

Here we focus on two cases: reionization driven by star-forming galaxies with a soft α_bk = 1.5 (α_s ∼ 2.4, within the commonly adopted range between 1 and 3; e.g., Bolton & Haehnelt 2007; Kuhlen & Faucher-Giguère 2012), and models of quasar-driven reionization with α_bk = 0.5 (corresponding to α_s ∼ 1.4; e.g., Telfer et al. 2002; Shull et al. 2012). We note that radiative transfer calculations have shown that the temperature increase from the passage of an ionization front may not depend strongly on the spectrum of the ionizing sources (D'Aloisio et al. 2018). We therefore constrain T_rei independently from α_bk.

8.3. Method

To be conservative, the constraints on instantaneous reionization models presented in this paper will be obtained using the thermal parameters in the lowest redshift bin (z = 4.2) only as upper limits because they may be affected by the extra heating due to the He ii reionization. For each combination of parameters (α_bk, z_rei, T_rei) we obtain three likelihood values corresponding to the redshifts of the observational constraints: ${L}_{\mathrm{model}}[\lt {T}_{0}^{z=4.2},\lt {u}_{0}^{z=4.2}]$, ${L}_{\mathrm{model}}[{T}_{0}^{z=4.6},{u}_{0}^{z=4.6}]$, and ${L}_{\mathrm{model}}[{T}_{0}^{z=5.0},{u}_{0}^{z=6.0}]$. These probabilities describe how well a model can simultaneously fit the observed values of T₀ and u₀ at each redshift. A given L_model is obtained by associating the model's T₀ and u₀ values with the probability derived from the corresponding full posterior distribution in Figures 11, 12, or 13 (middle panels).

The final probability of each model is computed by multiplying the independent likelihood values obtained at each redshift:

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{model}} & = & {L}_{\mathrm{model}}[\lt {T}_{0}^{z=4.2},\lt {u}_{0}^{z=4.2}]\\ & & \times \,{L}_{\mathrm{model}}[{T}_{0}^{z=4.6},{u}_{0}^{z=4.6}]\times {L}_{\mathrm{model}}[{T}_{0}^{z=5.0},{u}_{0}^{z=6.0}].\end{array}\end{eqnarray} \tag{ 17 }$$

8.4. Results

Before presenting the final results, we stress that, unlike previous attempts to constrain reionization using the instantaneous temperature alone (e.g., Theuns et al. 2002; Raskutti et al. 2012), the power of our approach relies on the simultaneous use of measurements of both T₀ and u₀. Figure 18 demonstrates how these measurements separately constrain the likelihood contours for our galaxy-driven reionization models. In the top panel the 68% and 95% probability contours are shown for the temperature constraints only, while in the middle panel they are given for the u₀ constraints only. The models that better fit the T₀ and u₀ data cover two different but intersecting regions in the T_rei versus z_rei parameter space. Applying both constraints simultaneously therefore reduces the allowed parameter space considerably (bottom panel of Figure 18).

Zoom In Zoom Out Reset image size

Figure 19 shows the final 68% and 95% 2D probability contours for models of instantaneous reionization driven by softer (α_bk = 1.5; green contours) and harder sources (α_bk = 0.5; blue contours). For softer sources the favored models are the ones with z_rei ∼ 8 and reionization temperature of 20,000 K ≲ T_rei ≲ 25,000 K. These temperatures are consistent with the values predicted by radiative transfer models (e.g., Miralda-Escudé & Rees 1994; McQuinn 2012; D'Aloisio et al. 2018). For harder sources the thermal data prefer earlier reionizations and lower temperatures (T_rei ≲ 20,000 K). These results are driven by the fact that for harder post-reionization ionizing backgrounds the IGM temperature needs more time to cool in order to match the relatively low values observed at z ≲ 5. We note that even lower T_rei would be needed to fit the observations if we included the contribution of He ii photoheating, which has been conservatively excluded. Radiative transfer calculations may disfavor T_rei ≲ 17,000 K, as this would imply reionization front speeds unexpectedly low even for the early stages of reionization (D'Aloisio et al. 2018). Some of the parameter space in Figure 19 preferred by harder sources may therefore be disfavored on physical grounds.

Zoom In Zoom Out Reset image size

In Figure 17, we also compare our constraints to the empirically calibrated UVB model recently presented by Puchwein et al. (2019). Their predictions for the instantaneous temperature are larger than our measurements and are inconsistent with our new constraints at around the 2σ level. Furthermore, the cumulative energy input into the IGM at mean density also exceeds our constraint at z = 6 and is again inconsistent at around 2σ. This suggests that there is slightly too much IGM heating at z > 6 in the fiducial Puchwein et al. (2019) model, under the assumption of a ΛCDM cosmology. This difference will be exacerbated further for models where AGN provide a substantial contribution to the photon budget for reionization.

Finally, we compare our results to constraints on the redshift of instantaneous reionization derived from the most recent Planck measurements of the Thompson scattering optical depth (Planck Collaboration et al. 2018). We marginalized over T_rei and α_bk (where α_bk was allowed to vary from 0.5 to 1.5) to obtain the 1D probability distribution on z_rei from the IGM thermal history. Figure 20 compares our marginalized constraints on z_rei (green solid line) to those derived for an instantaneous reionization from the Planck baseline optical depth constraint τ_e = 0.0544 ± 0.0073 (based on Planck TT, TE, EE+lowE+lensing; Planck Collaboration et al. 2018) (red dotted–dashed line). The two distributions are broadly consistent and have comparable constraining power, with ${z}_{\mathrm{rei}}\simeq {8.5}_{-0.8}^{+1.1}$ from the thermal history and ${z}_{\mathrm{rei}}\simeq {7.7}_{-0.7}^{+0.7}$ from the Planck results. The combined probability distribution (blue dashed line) gives ${z}_{\mathrm{rei}}\simeq {8.1}_{-0.5}^{+0.5}$. While we emphasize that the instantaneous reionization adopted in this section is simplistic, the proof of concept presented here demonstrates the potential of our observational constraints to inform reionization scenarios.

Zoom In Zoom Out Reset image size

In this work we have presented the first simultaneous constraints on the instantaneous temperature and integrated thermal history of the IGM at z > 4 and demonstrated how these results can be used to test different scenarios of hydrogen reionization. We have utilized a sample of 15 Keck/HIRES and VLT/UVES high-resolution and high-C/N spectra to obtain new measurements of the Lyα forest flux power spectrum over redshifts 4.0 ≲ z ≲ 5.2, for the first time pushing the measurement down to the smallest scales currently accessible to high-resolution quasar spectra at these redshifts (log(k/km⁻¹ s) ≤ −0.7). We fit the new flux power spectra to obtain robust constraints on the instantaneous IGM temperature, T₀, and integrated energy input per unit mass, u₀, marginalizing over the slope of the T–ρ relation and the effective optical depth, and assuming a ΛCDM cosmology.

In agreement with previous results from the curvature method (Becker et al. 2011), we find temperatures of T₀ ∼ 7000–8000 K and no strong temperature evolution over 4.2 ≲ z ≲ 5.0. Our first constraints on u₀ show a significant increase from u₀ ∼ 4.5 eV ${m}_{{\rm{P}}}^{-1}$ at z > 6 to 7.1 eV ${m}_{{\rm{P}}}^{-1}$ at z > 4.2–4.6. These results are consistent with a heating from reionization at z > 6 and residual photoionization heating over z ∼ 6–4.

Our constraints on T₀ and u₀ can be used to test any reionization scenario for which the temperature and the energy injection into the IGM can be calculated. As a proof of concept we analyzed simplistic, semianalytical models of instantaneous reionization. These toy models depend on three parameters: the IGM temperature reached during reionization (T_rei), the redshift of reionization (z_rei), and the spectral index of the post-reionization UVB (α_bk), which is related to the sources driving the reionization process. We find that our measurements prefer instantaneous reionization redshifts near z_rei ∼ 8 with T_rei ∼ 20,000 K for a relatively soft UVB dominated by ionizing photons from star-forming galaxies. Our fully marginalized constraints on the reionization redshift, ${z}_{\mathrm{rei}}\simeq {8.5}_{-0.8}^{+1.1}$, are moreover comparable to those from recent Planck results.

While tests of more realistic scenarios of reionization are left for future work, the proof of concept presented here demonstrates the potential of the IGM thermal history at high redshift to impose tight constraints on the timing—and possibly the sources—of reionization.

We thank Simeon Bird and Anson D'Aloisio for helpful conversations. E.B. and G.D.B. were supported by the National Science Foundation through grant AST-1615814. J.S.B. acknowledges the support of a Royal Society University Research Fellowship. This work is based on observations made at the W. M. Keck Observatory, which is operated as a scientific partnership between the California Institute of Technology and the University of California; it was made possible by the generous support of the W. M. Keck Foundation. It also includes observations made with ESO Telescopes at the La Silla Paranal Observatory under program ID 092.A-0770. E.B. thanks Michael Murphy and the Centre of Astrophysics and Supercomputing at Swinburne for granting the access to the Swinburne supercomputer facility during the preparation of this work. The hydrodynamical simulations used in this work were performed with supercomputer time awarded by the Partnership for Advanced Computing in Europe (PRACE) 8th Call. We acknowledge PRACE for awarding us access to the Curie supercomputer, based in France at the Tres Grand Centre de Calcul (TGCC). This work also made use of the DiRAC High Performance Computing System (HPCS) at the University of Cambridge. These are operated on behalf of the STFC DiRAC HPC facility. This equipment is funded by BIS National E-infrastructure capital grant ST/J005673/1 and STFC grants ST/H008586/1 and ST/K00333X/1. We thank Volker Springel for making P-GADGET-3 available.

In this appendix we review some of the steps of our analysis to check and quantify possible systematic uncertainties arising from the specific strategies adopted.

A.1. Rolling Mean

As described in Section 5.1, we computed the flux contrast δ_F of Equation (3) using a rolling mean along the entire Lyα forest region. The application of this technique on both simulated and observational lines of sight guarantees a fair comparison between models and the real measurement when the continuum level is unknown, but it may introduce possible bias when comparing our power spectrum with previous works in which the power was computed from continuum-normalized spectra. Using the simulations, we tested different averaging functions and window sizes in order to minimize the impact of the rolling mean on the power spectrum at the relevant redshifts and to verify the ability of the rolling mean to capture continuum fluctuations. As demonstrated below, we found that a 40 h⁻¹ cMpc boxcar rolling mean is able to recover the power at all relevant scales even in the presence of continuum fluctuations.

Figure 21 shows, for our three redshift bins, the comparison between the power spectrum computed from simulated data sets using the rolling mean technique (green dashed line) and using a fixed mean flux (black solid line). Both the synthetic samples of lines of sight used in this test have been created following the procedure described in Section 6. For the rolling mean model, we first imposed on each of the lines of sight a random continuum selected from the real continua fitted for the XQ-100 survey (López et al. 2016). We then run the 40 h⁻¹ cMpc boxcar rolling mean directly on the total Lyα + continuum flux. Differences between the two power spectra are shown in the bottom panel of each plot and compared with the statistical error characterizing our observational sample (green shaded region; see Section 5.9). We note that, at all redshifts, the discrepancies between the two models always lie well within the statistical error, with systematic uncertainties σ_roll typically ≲0.20σ_stat. We therefore do not expect our results to be sensitive to this averaging choice.

Zoom In Zoom Out Reset image size

A.2. Windowing Effects

As explained in Section 5.5, we compute the observational power spectrum in 20 h⁻¹ cMpc sections of Lyα forest. Dividing the spectra into many small regions may introduce artificial excess power at intermediate and small scales owing to a windowing effect. This effect does not not appear in the simulations because of their periodicity.

In Figure 22 we present the effect on the flux power spectrum of dividing the spectra into smaller sections. At each redshift bin, the power computed from the largest 40 h⁻¹ cMpc simulation box (baseline; solid black line) is compared with the power computed from the same simulation but dividing each of the native synthetic spectra into sections of 10 h⁻¹ cMpc (blue dot-dashed line) and 20 h⁻¹ cMpc (green dashed line). Changes in the power are reported as fractions of the baseline power spectrum values in the bottom panel of each plot. While cutting the spectra into 10 h⁻¹ cMpc sections introduces an excess of power at the small scales (log(k/km⁻¹ s) ≳ −1.1) of ∼20%–25%, the windowing effect for the 20 h⁻¹ cMpc sections is less significant, with variations in the power ≲8% at all scales. We therefore opted for the latter section size in our analysis.

Zoom In Zoom Out Reset image size

A.3. Noise Subtraction

In this appendix we test the noise power subtraction method described in Section 5.6.2 using synthetic data sets generated from simulations. For each synthesized line of sight we use the error and the flux arrays of one of our observed spectra to fit a linear correlation between the signal and the noise level. Using these correlations, we then add noise to the simulated samples in a flux-dependent way. We construct synthetic samples of lines of sight following the procedure described in Section 6.2 using our 20 h⁻¹ cMpc simulation box and compute the power spectrum with and without adding noise. In Figure 23 we show the noisy (red dot-dashed line) and noiseless (baseline; black solid line) power spectra obtained for each redshift bin.

Zoom In Zoom Out Reset image size

We finally applied the noise power subtraction method to the noisy model and compare the corrected power spectrum (green dashed line in Figure 23) to the noiseless baseline. The noiseless power is recovered with errors of ≲2%. This suggests that this step of the data analysis is not introducing relevant systematic effects in the final results.

A.4. Instrumental Resolution

In this appendix we review the effect of uncertainties in the instrumental resolution when smoothing the synthetic spectra to match the resolution of the observed data. We calibrate the synthetic Lyα forest spectra using different instrumental resolutions, taking our nominal values of FWHM = 6 km s⁻¹ as a baseline.

Figure 24 shows the variations in the power spectrum expected for a change of +10% (red dot-dashed line) and −10% (green dashed line) in the observed spectral resolution. We note that the power changes by ≲5% at all scales (see lower panel of each plot).

Zoom In Zoom Out Reset image size

A.5. Masking Correction Function

In this appendix we analyze the systematic uncertainties arising from the choice of a particular simulation in computing the masking correction, C_m(k), described in Section 5.6.1. For the fiducial masking correction we adopted the simulation S20_1z15 of Table 2. To demonstrate that this particular choice of model (with T₀ ∼ 7500 K and γ ∼ 1.5 at the redshifts of interest) does not relevantly affect the final power spectrum measurements, we compare our fiducial results with the measurements obtained when computing the masking correction from post-processed runs with different values of T₀ and γ. In particular, we tested two extreme cases: a colder model, ${C}_{m}^{\mathrm{test}1}(k)$, with T₀ ∼ 4000 K and γ ∼ 1.5, and an isothermal model, ${C}_{m}^{\mathrm{test}2}(k)$, with T₀ ∼ 7500 K and γ ∼ 1.0.

Figure 25 shows the power spectrum obtained using different masking correction functions for the redshift bins considered in our analysis. In all the cases, computing the masking correction using models with different thermal parameters affects the small scales only mildly, with variations in the power of ≲4% for scales log(k/km⁻¹ s) ≳ −1.1. We therefore do not expect that uncertainties in the masking correction function will relevantly affect our final constraints.

Zoom In Zoom Out Reset image size

A.6. Metal Correction

The Lyα forest region is affected by narrow metal line contaminants that may increase the flux power at the small scales. To correct for the metal contribution, we subtract the metal power evaluated in Section 5.7 from the final power spectrum measurements. The effect of the metal subtraction on the power spectrum measurements is presented in Figure 26 for our three redshift bins. As expected, only the small scales (log(k/km⁻¹ s) ≳ −1.0) are significantly affected. The variations are relatively small (<10%) and well within the statistical errors of our final measurements.

Zoom In Zoom Out Reset image size

In this appendix we test how strongly the choice of simulation model for the covariance matrix regularization affects the final constraints on T₀ and u₀. In principle, off-diagonal coefficients of the covariance matrix will mildly depend on the shape of the power spectrum and therefore on the thermal parameters characterizing the models.

Figure 27 shows the results for T₀ and u₀ obtained when using a covariance matrix derived from the fiducial model S40_1z15 (Covariance matrix_1) and from the S40_1z9 model (Covariance matrix_2). Both the nominal results and error estimates shows only modest (≲3%) changes.

Zoom In Zoom Out Reset image size

In this appendix we compare our power spectrum measurements to those from Viel et al. (2013a). While Viel et al. (2013a) included somewhat lower resolution spectra from Magellan/MIKE, for consistency we limit our comparison to their measurements obtained with the HIRES spectrograph.

We consider two main differences between our new estimates and the older Viel et al. (2013a) power spectra: flux contrast estimators and cosmic variance. First, Viel et al. (2013a) normalized the spectra using a spline continuum estimate and then computed the mean flux in large sections of data. To test whether this affected the results, we recomputed the Viel et al. (2013a) power spectra by applying the procedure described in Section 5 to the Viel et al. data, using the same sections of spectra. Note that, given the somewhat larger scales probed by Viel et al., we do not expect any relevant effect due to noise or resolution. We verified that we were able to reproduce the previous results with good precision and no significant bias was introduced by the different flux contrast estimators.

We next consider whether sample size may be playing a role. In each of our redshift bins the number of independent lines of sight contributing to our measurement is always more than double the number in Viel et al. (2013a). In particular, for the z = 4.2 bin we used 12 lines of sight versus 4 in Viel et al. (2013a), for the z = 5.0 bin we used 12 versus 5, while at z = 4.6, where the largest differences between the power spectra are seen, we used 15 quasars versus 5 in the previous work.

To determine whether cosmic variance can explain the discrepancy between our results and those of Viel et al. (2013a), we computed the power spectrum from subsamples of our data. Figure 28 shows the comparison between the Viel et al. (2013a) power spectra (black dashed lines) and the 68% (darker shaded regions) and 95% (lighter shaded regions) contours of the distribution of power spectrum realizations obtained from a Monte Carlo sampling of our lines of sight. For each redshift bin we randomly select from our sample the same number of lines of sight contributing to the Viel et al. (2013a) measurement and use them to compute the power spectrum. Repeating the process for many (∼200) realizations, we verify that the Viel et al. (2013a) results fall within this distribution. At z = 4.2 and z = 5.0 most of the Viel et al. (2013a) values fall within the 68% contours. For the z = 4.6 bin the agreement is slightly worse but still largely within the 95% region. We note that the errors in different k bins are correlated, as discussed in Section 5.9. We further note that our sample size is still relatively modest, and our Monte Carlo technique is likely to underestimate the cosmic variance at the 95% level.

Zoom In Zoom Out Reset image size

In this appendix we compare our power spectrum measurement at z = 4.2 to the one from Iršič et al. (2017a). These authors use spectra from the XQ-100 Legacy Survey (López et al. 2016), collected using the lower-resolution VLT/X-Shooter spectrograph (R ∼ 11 km s⁻¹ corresponding to an FWHM = 26 km s⁻¹ for the VIS arm). We therefore limited our comparison to the power spectrum scales where the resolution correction for the Iršič et al. (2017a) measurements is ≲20%, corresponding to log(k/km⁻¹ s) < −1.4.

Figure 29 presents the comparison between the Iršič et al. (2017a) power spectrum (black points) and our measurements (blue points) at z = 4.2. For both data sets we also include the final 68% uncertainties (colored shaded regions). Even if the X-Shooter power spectrum shows tendentially lower values than ours, the two measurements are consistent within the 68% uncertainties at all scales but one point (at the scale log(k/km⁻¹ s) = −1.92). A possible explanation for the offset in the power at scales log(k/km⁻¹ s) ≳ −1.9 may be differences in the spectra samples and in the redshift coverage of the z = 4.2 bin: while we include spectra falling within a broad Δz = 0.4 redshift bin, Iršič et al. (2017a) adopt a narrower Δz = 0.2 bin. Nevertheless, as for the Viel et al. (2013a) power spectrum, we verified that, using the same X-Shooter spectra, we were able to reproduce the Iršič et al. (2017a) results at the considered scales without introducing any significant bias due to possible differences in the analysis.

Zoom In Zoom Out Reset image size

In Figure 30 we compare our IGM temperature measurements to previous constraints from the literature over the redshift range covered by our analysis. We note that the temperature values from Iršič et al. (2017b) (pink stars) and Garzilli et al. (2017) (brown circles) and the z ≳ 4.2 results from Walther et al. (2018) (orange squares) are all obtained from the flux power spectrum measurements of Viel et al. (2013a), although calibrated with different sets of simulations; therefore, they are not fully independent. Among these constraints, the larger error bars reported for Iršič et al. (2017b) reflect a 1.5 correction factor applied by these authors to the nominal errors associated with the Viel et al. (2013a) measurements.

Zoom In Zoom Out Reset image size

While generally within the 68% uncertainties of our new results, previous estimates from the flux power spectrum statistic suggest a significant decrease in temperature at z ≳ 4.6. Our thermal constraints, however, obtained from a larger and higher-quality sample of high-resolution spectra, do not show any strong evolution in the temperature at these redshifts (see also Appendix C for a comparison between our recent flux power spectrum measurements and the ones from Viel et al. 2013a).

In this appendix we examine the convergence of the Lyα flux power spectrum in the simulations used in this work. We used multiple simulations with the thermal history model 1_z15 in Table 2. All of the synthetic Lyα forest lines of sight were produced using the procedure described in Section 6.1. The flux power spectrum was then computed as in Section 6.2.

The tests are shown in Figure 31, where the convergence with box size for a fixed mass resolution (M_gas = 9.97 ×10⁴ h⁻¹M_⊙) is reported in the left column and the convergence with mass resolution for a fixed box size (L = 10 h⁻¹ cMpc) is displayed on the right. The results show that a small correction for both box size and mass resolution needs to be applied to the power spectra derived from our nominal 10 h⁻¹ cMpc simulations. When increasing the box size, the power decrease up to ∼15%, particularly at very small scales (log(k/km⁻¹ s) > −1). In contrast, when increasing the mass resolution, the power toward small scales (log(k/km⁻¹ s) ≳ −1.4) increases progressively, reaching a correction of ∼15% at z = 5.0 for our nominal mass resolution. Because the corrections are in opposite directions, the final scale factor for box and mass resolution convergence is ≲5% at all scales. However, we note that in principle the mass resolution convergence may depend on the underlying IGM density structure and, consequently, on the choice of a particular thermal history model. We therefore verified that the entity of this possible systematic effect was negligible when compared to the 68% statistical uncertainties characterizing the observational data and that the final thermal constraints were not affected by the numerical corrections.

Zoom In Zoom Out Reset image size

We further verified, using lower mass resolution simulations, that increasing the box size up to L = 160 h⁻¹ cMpc for a fixed mass resolution does not introduce additional power at the scales considered in this work.

As explained in Section 6.1, when constructing mock samples, we account for the mild redshift evolution of the mean flux along the line of sight by initially rescaling the effective Lyα optical depth using Equation (10). In this appendix we show how the choice of this relation for the τ_eff evolution, while somewhat arbitrary, represents a reasonable transition between the measurements of Becker et al. (2013) and the newer results of Bosman et al. (2018). Figure 32 presents the comparison among the different τ_eff evolutions. The results from the analysis of Becker et al. (2013) (gray triangles) at z < 4 are smoothly connected to the most recent measurements of Bosman et al. (2018) (red squares) at z > 5 by the fiducial fit adopted in this work (green dashed line). For comparison, we also show the constraints on the optical depth obtained from our MCMC chains and reported in Table 5 (green points). We recover values broadly consistent with Equation (10) for the three redshift bins.

Zoom In Zoom Out Reset image size

In this appendix we describe the test performed to verify the interpolation scheme implemented in the MCMC analysis. For this we remove one model from the interpolation grid (in the example below we excluded the model S10_0.55z9 of Table 2) and test how well the thermal parameters for this model are recovered when it is used to generate artificial data.

Figure 33 displays the correct values (red squares) overlaid on the parameter constraints obtained from the MCMC analysis at z = 5.0. The thermal parameters T₀ and u₀ are recovered accurately by the analysis, with discrepancies ≲5%. The correct values of γ and τ_eff fall within the 1σ probability distribution, although the peaks of their posterior distributions are somewhat biased toward lower values. As expected, the power spectrum at these redshifts is not sensitive enough to break the degeneracy between γ and τ_eff. The poor constraints on γ, however, do not affect our constraints on T₀ and u₀.

Zoom In Zoom Out Reset image size

We have also tested how well our interpolation scheme was able to recover the thermal parameters of a completely independent model. We fit the power spectrum extracted from a simulation using the UVB model of Puchwein et al. (2019) assuming nonequilibrium ionization, and we verified that the values of T₀ and u₀ (see Figure 17) were recovered within the 68% uncertainties given by our MCMC method. This gives some reassurance that, as intended, our results do not significantly depend on the specific thermal histories adopted for the modeling in this work.

For illustrative purposes we show in Figure 34 the evolution of the thermal parameters u₀ and T₀ for all the simulations listed in Table 2. While these models are not meant to represent realistic reionization scenarios, they provide a wide range of thermal histories and can be used to explore the thermal state of the IGM in a relatively model-independent way.

Zoom In Zoom Out Reset image size

In this appendix we show the relationship between u₀ and the characteristic real-space flux power cutoff scale, λ_P, as defined by Kulkarni et al. (2015). At each redshift we compute λ_P for all the models of Table 2 following the method described in Kulkarni et al. We then fit a relationship between the corresponding u₀ computed over the fiducial redshift range. Figure 35 shows the best-fitting relationship between u₀ and λ_P for the redshifts considered in this work. While a certain level of scatter about the fit is present at all redshifts, there is clearly a positive correlation between the two variables. Using the current constraints on u₀, we can then attempt to obtain a rough estimate of λ_P (green square with error bars).

Zoom In Zoom Out Reset image size

We note that these estimates for λ_P are smaller than the recent constraints at 2 < z < 4 from quasar pairs presented by Rorai et al. (2017a), suggesting that the pressure smoothing scale increases toward lower redshift as the IGM is photoheated further (see the instantaneous temperature measurements presented by Becker et al. 2011)

In Tables 6–8 we report the power spectrum measurements obtained in this work for the three redshift bins centered at z = 4.2, 4.6, and 5.0. In each table the values of the power spectrum obtained with (Column (3)) and without (Column (2)) instrumental resolution and pixel size correction (R.C.) are reported for each scale (Column (1)). The corresponding 68% uncertainties are shown in Column (4). The covariance matrices for the power spectrum measurements may be found in the online version of this article.