Constraining Reionization with the $z \sim 5-6$ Lyman-$\alpha$ Forest Power Spectrum: the Outlook after Planck

TL;DR

This paper investigates how the timing and heat injection of reionization affect the Lyman-alpha forest power spectrum at redshifts 5-6, using hydrodynamical simulations to constrain reionization models in light of Planck data.

Contribution

The study demonstrates that the Lyman-alpha forest power spectrum can reveal reionization history details, helping distinguish between galaxy-driven and quasar-driven reionization scenarios.

Findings

Models consistent with Planck's optical depth favor galaxy-driven reionization.

The small-scale cutoff in the power spectrum is sensitive to the IGM thermal history.

High-resolution spectra can differentiate between reionization scenarios at z~6.

Abstract

The latest measurements of CMB electron scattering optical depth reported by Planck significantly reduces the allowed space of HI reionization models, pointing toward a later ending and/or less extended phase transition than previously believed. Reionization impulsively heats the intergalactic medium (IGM) to $\sim10^4$ K, and owing to long cooling and dynamical times in the diffuse gas, comparable to the Hubble time, memory of reionization heating is retained. Therefore, a late ending reionization has significant implications for the structure of the $z\sim5-6$ Lyman-$\alpha$ (ly$\alpha$) forest. Using state-of-the-art hydrodynamical simulations that allow us to vary the timing of reionization and its associated heat injection, we argue that extant thermal signatures from reionization can be detected via the Ly$\alpha$ forest power spectrum at $5< z<6$. This arises because the small-scale cutoff in the power depends not only the the IGM temperature at these epochs, but is also particularly sensitive to the pressure smoothing scale set by the IGM full thermal history. Comparing our different reionization models with existing measurements of the Ly$\alpha$ forest flux power spectrum at $z=5.0-5.4$, we find that models satisfying Planck's $\tau_e$, constraint favor a moderate amount of heat injection consistent with galaxies driving reionization, but disfavoring quasar-driven scenarios. We explore the impact of different reionization histories and heating models on the shape of the power spectrum, and find that they can produce similar effects, but argue that this degeneracy can be broken with high enough quality data. We study the feasibility of measuring the flux power spectrum at $z\simeq 6$ using mock quasar spectra and conclude that a sample of $\sim10$ high-resolution spectra with an attainable signal-to-noise ratio will allow distinguishing between different reionization scenarios.

Figures

Tables

Table 1. Table 1 Summary of Simulations.

Sim	H i reionization	$z_{\mathrm{reion},\text{H}\text{i}}$	$z_{\mathrm{reion},\text{H}\text{i}}^{0.5}$	$\Delta z$	${\tau }_e$	$\Delta T_{\text{H}\text{i}}$	$u_0\left(z=4.9\right)$
						(K)	(eV $m_p^{-1}$)
VeryEarlyR	Very Early	$9.70$	$10.75$	$2.59$	$0.0917$	$2\times {10}^4$	$7.86$
EarlyR	Early	$7.75$	$8.80$	$2.59$	$0.0698$	$2\times {10}^4$	$5.52$
MiddleR	Middle	$6.55$	$7.60$	$2.59$	$0.0574$	$2\times {10}^4$	$4.24$
LateR	Late	$6.00$	$7.05$	$2.59$	$0.0520$	$2\times {10}^4$	$3.85$
MiddleR-fast	Fast middle	$6.60$	$7.10$	$0.89$	$0.0498$	$2\times {10}^4$	$4.18$
MiddleR-cold	Middle	$6.55$	$7.60$	$2.59$	$0.0574$	$1\times {10}^4$	$3.16$
MiddleR-warm	Middle	$6.55$	$7.60$	$2.59$	$0.0574$	$3\times {10}^4$	$5.60$
MiddleR-hot	Middle	$6.55$	$7.60$	$2.59$	$0.0574$	$4\times {10}^4$	$7.08$

Equations

$$⟨x_{\text{H}\text{textsc}ii}⟩=\{\begin{array}{l}0.5+0.5\times g(1/n_1,|z-z_{\mathrm{reion},\text{H}\text{textsc}\mathrm{i}}^{0.5}{|}^{n_1}),\\ z<=z_{\mathrm{reion},\text{H}\text{textsc}\mathrm{i}}^{0.5}\\ 0.5-0.5\times g(1/n_2,|z-z_{\mathrm{reion},\text{H}\text{textsc}\mathrm{i}}^{0.5}{|}^{n_2}),\\ z>z_{\mathrm{reion},\text{H}\text{textsc}\mathrm{i}}^{0.5}\end{array}$$

Full text

1. Introduction

How and when the first luminous sources reionized diffuse baryons in the intergalactic medium (IGM) is one of the most fundamental open questions in cosmology. Recently, the Planck collaboration have released new tighter constraints on reionization from cosmic microwave background (CMB) observations (Planck Collaboration et al. 2016b, c). Using for the first time the low-multipole $EE$ data from Planck-HFI, the Planck team has significantly improved our constraints on the cosmic reionization optical depth, $\tau_{\rm e}$, finding $\tau_{\rm e}=0.058\pm 0.012$ (Planck Collaboration et al. 2016c).

The reionization of H i by the UV background from galaxies and/or quasars results in the highly ionized IGM probed at $z\lesssim 6$ by observations of the Ly$\alpha$ forest (McQuinn 2016). Indeed, observations of complete Gunn-Peterson absorption in the spectra of many of the highest $z\sim 6$ quasars, along with the steep rise of both the Ly$\alpha$ optical depth and its sightline-to-sightline scatter with redshift, has led to the consensus that we are witnessing the end of reionization only at $z\sim 6$ (Fan et al. 2006; McGreer et al. 2015; Becker et al. 2015). However, Ly$\alpha$ opacity can only set lower limits on the redshift of reionization $z\gtrsim 6$, because the overly sensitive Ly$\alpha$ transition saturates for volume-averaged neutral fractions $\langle x_{\textrm{H}\,\textsc{i}}\rangle\gtrsim 10^{-4}$, which is far too small to pinpoint the redshift of reionization. While new constraints have emerged from the possible presence of a Ly$\alpha$ IGM damping wing in the highest redshift known quasar at $z=7.1$ (Mortlock et al. 2011; Simcoe et al. 2012; Greig et al. 2016), and the decreasing strength of Ly$\alpha$ emission lines in $z\sim 6-7$ galaxies (Caruana et al. 2014; Schmidt et al. 2016; Sadoun et al. 2017), the resulting constraints on $\langle x_{\rm HI}\rangle$ are degenerate with the intrinsic properties of the high-$z$ quasars and galaxies that they have been deduced from. We are in need of another technique to probe when reionization occurred.

During reionization, ionization fronts propagate supersonically through the IGM, impulsively heating gas to $\sim 10^{4}$ K (Abel & Haehnelt 1999; Davies et al. 2016). The integrated energy balance of heating and inverse-Compton and adiabatic cooling then gives rise to a power law temperature-density relation, $T=T_{0}(\rho/\bar{\rho})^{\gamma-1}$ (Miralda-Escudé & Rees 1994; Hui & Gnedin 1997; Hui & Haiman 2003; Meiksin 2009; McQuinn et al. 2009; McQuinn & Upton Sanderbeck 2016). Another important physical ingredient to describe the thermal state of the IGM is the gas pressure support that produces an effective three-dimensional smoothing of the baryon distribution relative to the dark matter at a characteristic scale, $\lambda_{\rm P}$. In an expanding universe with an evolving thermal state at a given epoch, this scale depends on the entire thermal history of IGM because fluctuations at earlier times expand or fail to collapse depending on the IGM temperature at that epoch (Gnedin & Hui 1998; Rorai et al. 2013; Kulkarni et al. 2015; Oñorbe et al. 2017; Rorai et al. 2017a). At redshift $z$, the level of pressure smoothing depends not on the prevailing pressure/temperature at that epoch, but rather on the temperature of the IGM in the past. The IGM pressure-smoothing scale, $\lambda_{\rm P}$, thus provides an integrated record of the thermal history of the IGM, and is sensitive to the timing of and heat injection by reionization events.

Measurements of the statistical properties of the Ly$\alpha$ forest are sensitive to the thermal state of the IGM through the thermal Doppler broadening of absorption lines, as well as the pressure-smoothing. The standard approach has been to compare measurements of different statistics to cosmological hydrodynamical simulations (Zaldarriaga et al. 2001; Theuns et al. 2002b; Viel et al. 2009; Lidz et al. 2010; Becker et al. 2011; Garzilli et al. 2012; Rorai et al. 2013; Iršič & Viel 2014; Rorai et al. 2017b, a) to deduce the thermal parameters (e.g., $T_{0}$, $\gamma$ or $\lambda_{\rm P}$) that best describe the IGM thermal state. As a larger number of high-resolution spectra of quasars have become available at higher redshifts $z\gtrsim 4$, the same approach has been applied to study the thermal state of the IGM at these redshifts, where the Ly$\alpha$ forest is more sensitive to the timing and nature of hydrogen reionization (Theuns et al. 2002b; Hui & Haiman 2003; Furlanetto & Oh 2009; Cen et al. 2009; Becker et al. 2011; Viel et al. 2013a; Lidz & Malloy 2014; Garzilli et al. 2015; Nasir et al. 2016).

In the light of the new Planck constraints on reionization, as well as the increasing number of quasars discovered at high-$z$ (e.g. Bañados et al. 2014; Matsuoka et al. 2016; Bañados et al. 2016), it is pertinent to revisit what Ly$\alpha$ observations with current and upcoming facilities can tell us about H i reionization. This work aims to explore in detail the possibilities of using the Ly$\alpha$ forest 1D flux power spectrum at high-$z$ to constrain HI reionization. For this we have used a new method that has recently been introduced by Oñorbe et al. (2017). It builds on the Haardt & Madau (2012) model, but enables one to vary the redshifts of H i and He ii reionization, as well as their associated heat injection, allowing one to consistently simulate a more diverse range of reionization histories. This allows for a much comprehensive and consistent exploration of the space of thermal parameters than previously done. We present here the results of the 1D flux power spectrum at high-$z$ of a new set of hydrodynamical simulations that also improve the resolution used in previous studies at these redshifts. Additionally, we compare these power spectra with some recent measurements at these redshifts (Viel et al. 2013a).

The structure of this paper is as follows. In Section 2 we describe the characteristics of our hydrodynamical simulations and the different H i reionization models studied in this work. Section 3 presents the thermal evolution of the different reionization models obtained from the cosmological hydrodynamical simulations. In Section 4 we compare the 1D flux power spectrum of each model at $5\leq z\leq 6$ as well as with the best available observations. We discuss in Section 5 the relevance of our findings in the context of current observational and theoretical limitations. We conclude by presenting a summary of our results and an outlook in Section 6. In Appendix A we perform a set of convergence test for the optical depth and the 1D flux power spectrum at $5\leqslant z\leqslant 6$.

Throughout this paper we assumed a flat $\Lambda$CDM cosmology with the following fundamental parameters: $\Omega_{\rm m}=0.3192$, $\Omega_{\rm\Lambda}=0.6808$, $\Omega_{\rm b}=0.04964$, $h=0.67038$, $\sigma_{\rm 8}=0.826$ and $n_{\rm s}=0.9655$. These values agree within one sigma with the latest cosmological parameter constrains from the CMB (Planck Collaboration et al. 2016a, c). The mass abundances of hydrogen and helium ($X_{\rm p}=0.76$ and $Y_{\rm p}=0.24$) were chosen to be in agreement with the recent CMB observations and Big Bang nucleosynthesis (Coc et al. 2013).

2. Simulating Reionization Histories in Light of New Planck Constraints

The cosmological hydrodynamical simulations used in this work were performed using the Nyx code (Almgren et al. 2013). Application of Nyx to studies of the Ly$\alpha$ forest, and its convergence and resolution requirements are discussed in Lukić et al. (2015). We refer to these two works for more details of the numerical implementation, accuracy, and code performance. To generate the initial conditions for the simulations, we have used the music code (Hahn & Abel 2011), with the transfer function for our cosmological model obtained from camb (Lewis et al. 2000; Howlett et al. 2012). All simulations discussed in this work used the same initial conditions and have a box size of length $L_{\rm box}=20$ Mpc h*-1* and $1024^{3}$ resolution elements.

As is standard in hydrodynamical simulations that model the Ly$\alpha$ forest, all cells are assumed to be optically thin to radiation. Thus, radiative feedback is accounted for via a spatially uniform, but time-varying ultraviolet background (UVB) radiation field, input to the code as a list of photoionization and photoheating rates that vary with redshift (e.g. Katz et al. 1996). In order to simulate each reionization model discussed here, we have used the method presented in Oñorbe et al. (2017), which allows us to vary the timing and duration of reionization, and its associated heat injection, enabling us to simulate a diverse range of reionization histories. This method allows us to create the H i, He i and He ii photoionization and photoheating rates, which are inputs to the Nyx code, by volume averaging the photoionization and energy equations. This method requires that each reionization event is defined by the ionization history with redshift, e.g. $x_{\textrm{H}\,\textsc{i}}(z)$, and its associated total heat injection, $\Delta T$, which depends on the spectral shape and abundance of the ionizing sources, and the opacity of the IGM (Abel & Haehnelt 1999; McQuinn 2012; Davies et al. 2016; Park et al. 2016). We direct the reader to Oñorbe et al. (2017) for the details of this method.

In order to determine the reionization histories for our simulations, we explore the relevant range of reionization models considering the last $\tau_{\rm e}$ measurements by Planck and the Ly$\alpha$ optical depth at high redshift, which set a lower limit for H i reionization at $z=6$ (Fan et al. 2006; McGreer et al. 2015; Becker et al. 2015). In all simulations we also assumed the same He ii reionization model ending at $z=3$ and which does not produce a significant increase in the IGM temperature until $z<5$. This model is the same as assumed in the Faucher-Giguère et al. (2009) UVB model, and since we will be comparing with data at $z\geq 5$ this assumption will not impact our results (see Oñorbe et al. 2017, for more details). In this work we consider five model H i reionization histories constructed using the analytical formula presented in Oñorbe et al. (2017) chosen to match the results of radiative transfer simulations (Pawlik et al. 2009),

[TABLE]

where $g$ is the incomplete gamma function, $n_{1}=50$, $n_{2}=1$ and $z_{\rm reion,\textrm{H}\,\textsc{i}}^{0.5}$ is a free parameter that sets the redshift where $x_{\textrm{H}\,\textsc{ii}}(z_{0})=0.5$.

We run an early, middle, and late H i reionization history (EarlyR, MiddleR, LateR), which have a specific reionization redshifts (defined as the redshift where $\langle x_{\textrm{H}\,\textsc{ii}}\rangle=0.999$) of $z_{\rm reion,\textrm{H}\,\textsc{i}}=7.75$, $6.55$, and $6.0$ respectively, and are within $1\sigma$ of the Planck CMB measurements. We also run two more models, one with a very early reionization (VeryEarlyR, $z_{\rm reion,\textrm{H}\,\textsc{i}}=9.70$ which is $3\sigma$ discrepant with the Planck measurement) and a faster reionization (MiddleR-fast, $z_{\rm reion,\textrm{H}\,\textsc{i}}=6.55$). A summary of all the relevant parameters used in the runs presented in this work is shown in Table 1 along with the naming conventions we have adopted.

The full reionization history of each of our models is plotted in the upper left panel of Figure 1, shown as the redshift evolution of the electron fraction given by $x_{\rm e}=n_{\rm e}/n_{\textrm{H}}=(1+\chi)\langle x_{\textrm{H}\,\textsc{ii}}\rangle+\chi\langle x_{\textrm{He}\,\textsc{iii}}\rangle$ where $\chi=Y_{p}/(4X_{p})$ and $X_{p}$ and $Y_{p}$ are the hydrogen and helium mass abundances, $x_{\textrm{H}\,\textsc{ii}}(z)$ is the hydrogen ionized fraction, and $x_{\textrm{He}\,\textsc{iii}}(z)$ is the fraction of helium that is doubly ionized111Throughout this paper we made the standard assumption that He i reionization is perfectly coupled with that of H i.. The lower left panel of Figure 1 shows the evolution of the cosmic reionization optical depth, $\tau_{\rm e}$, for each of these models. In the right panel of Figure 1 we compare the newest Planck measurements of the CMB polarization EE power spectrum222Throughout this paper we adopt the convention $D_{\ell}=\ell(\ell+1)C_{\ell}/2\pi$. low multipoles (Planck Collaboration et al. 2016b) to these reionization models, where we have computed the EE power spectrum using the class Boltzmann code (Blas et al. 2011). For the observed low multipoles we show the unbiased QML2 results from Planck Collaboration et al. (2016b, black circles), illustrating the impressively high precision achieved by these CMB polarization measurements, which significantly reduces the allowed range of models.

To build the reionization models, we also need to assume the associated total heat injection, $\Delta T$, during H i reionization, which depends on the spectral shape and abundance of the ionizing sources, and the opacity of the IGM (Abel & Haehnelt 1999; Tittley & Meiksin 2007; McQuinn 2012; Davies et al. 2016; Park et al. 2016). To run all the simulations with the reionization model described above, we assumed $\Delta T_{\textrm{H}\,\textsc{i}}=2\times 10^{4}$ K which is the standard value obtained in galaxy-driven H i reionization models using 1D radiative transfer simulations (e.g. McQuinn 2012). Quasar-driven scenarios give higher heat injection values, $\Delta T_{\textrm{H}\,\textsc{i}}\sim 4\times 10^{4}$. Thus, in order to study the effect of different total heat input during H i reionization, we run three more simulations with the same H i reionization model as MiddleR, but varying the $\Delta T_{\textrm{H}\,\textsc{i}}$ parameter: MiddleR-cold ($\Delta T_{\textrm{H}\,\textsc{i}}=1\times 10^{4}$ K), MiddleR-warm ($\Delta T_{\textrm{H}\,\textsc{i}}=3\times 10^{4}$ K), and MiddleR-hot ($\Delta T_{\textrm{H}\,\textsc{i}}=4\times 10^{4}$ K).

3. Reionization-dependent Thermal Histories

Changing the timing and duration of reionization and its associated heat injection will manifest as changes in the evolution of the parameters governing the thermal state of the IGM. In Figure 2 we present the resulting thermal histories for all of these simulations. The upper panel shows the evolution of temperature $T(\Delta_{\bigstar})$ at the ‘optimal’ overdensity $\Delta_{\bigstar}$ probed by curvature measurements of the Ly$\alpha$ forest (see Becker et al. 2011; Boera et al. 2014), where we calculate the optimal density at each redshift using the functional form of $\Delta_{\bigstar}(z)$ given by Becker et al. (2011). The evolution of thermal parameters, $\gamma$ and $T_{0}$, governing the density-temperature relation, are shown in the second and third panel from the top, determined by fitting the distribution of densities and temperatures in the simulation following the linear least-squares method described in Lukić et al. (2015)333 Changing the thresholds used to do the fit within reasonable IGM densities produces differences just at a few percent level (see Lukić et al. 2015, for similar conclusions) and in any case it does not affect the conclusions presented in this work. We also found no relevant effects in the main results of this paper when we employed the fitting approach used in Puchwein et al. (2015)..

The evolution of the pressure-smoothing scale, $\lambda_{\rm P}$, with redshift is shown in the bottom panel. To characterize the pressure-smoothing scale in all our simulations we have followed the approach described by Kulkarni et al. (2015). These authors define a pseudo real-space Ly$\alpha$ flux field, which is the same as the true Ly$\alpha$ forest flux, but without redshift space effects such as peculiar velocities and thermal Doppler broadening. This field naturally suppresses the dense gas that would otherwise dominate the baryon power spectrum, making it robust against the poorly understood physics of galaxy formation and revealing the pressure-smoothing in the diffuse IGM.

Inspection of the left panel of Figure 2 reveals that simulations with different reionization histories but the same heat injection during reionization, $\Delta T$, all share a very similar $T_{0}$, $T(\Delta_{\bigstar})$ and $\gamma$ evolution at $z=5-6$. This is because once reionization is completed, the IGM thermal state asymptotes to a tight power-law temperature-density relation driven mainly by the photoheating rate and accelerated by Compton and adiabatic cooling. The time to converge to these asymptote values is around $\Delta z\sim 1-2$ (a few hundred Myr) once reionization is completed and mainly depends on the amount of heat injected during reionization (McQuinn & Upton Sanderbeck 2016). Since all these models share the same photoionization and photoheating values, once reionization is completed they all converge to the same thermal state at lower redshift $z\lesssim 6$. However, their pressure-smoothing scale, $\lambda_{\rm P}$, remains very different at these and lower redshifts. Models in which reionization happened at earlier times have a larger pressure-smoothing scale. As discussed above, this results from the dependence of the IGM pressure-smoothing scale on the full thermal history (Hui & Haiman 2003; Kulkarni et al. 2015; Oñorbe et al. 2017) and not just on the instantaneous temperature, and it will have important consequences for the statistics of the Ly$\alpha$ forest (Nasir et al. 2016; Oñorbe et al. 2017).

The right panels of Figure 2 show the thermal histories for simulations where we fixed the reionization history, but varied the total heat injection during H i reionization. Models with more heat injection give rise to higher temperatures and slightly lower $\gamma$ values not only during reionization, but also at later times while the IGM is still reaching their asymptote values. Note that similarly as for the different reionization history models, these models also produce a larger pressure-smoothing scale, and these differences persist even at lower redshifts long after the other thermal parameters, $\gamma$ and $T_{0}$, have relaxed to their asymptotic values. Once reionization is completed all these models share the same photoionization and photoheating rates and therefore $\gamma$, $T_{0}$, and $T(\Delta_{\bigstar})$ thermal parameters asymptote to the same values much faster than the pressure-smoothing scale, which retains memory of the thermal history.

The symbols with error bars in Figure 2 indicate recent observational constraints on the parameters governing the thermal state of the IGM at high redshift. In particular, the purple square is the Lidz et al. (2010) measurement of $T_{0}$ using wavelets at $z=4.20$ and black squares are the Becker et al. (2011) measurements of $T(\Delta_{\bigstar})$ based on the curvature statistic. Gray circles and black diamonds represent the joint fits to $\gamma$ and $T_{0}$ given by Viel et al. (2013a) and Garzilli et al. (2015) respectively, using the 1D flux power spectrum.444These measurements are marginalized over the mass of a warm dark matter particle. Viel et al. (2013a) and Garzilli et al. (2015) used different fitting approaches, but both used the same grid of hydrodynamical simulations in which the standard reionization redshift for the runs was $z_{\rm reion}=12$ and the lowest reionization redshift considered in the grid, by including one simulation, was $z_{\rm reion}=8$. While the thermal parameters measured by Becker et al. (2011), Viel et al. (2013a) and Garzilli et al. (2015) appear consistent with the models discussed here, the Lidz et al. (2010) $T_{0}$ measurement suggests a significantly hotter IGM at $z\sim 4$. The origin of this disagreement is unclear, but may result from differences in the methods used by these authors and/or the different hydrodynamical simulations compared to the data. Our work aims to shed more light on this issue by comparing the 1D power spectrum measurements at high redshift with an improved set of simulations.

4. The 1D Flux Power Spectrum at High Redshift $z\sim 5-6$

In order to explore the possibility of discriminating between the reionization models presented here with Ly$\alpha$ forest measurements, we calculate the 1D flux power spectrum, $P(k)$, for each simulation at $z=5.0$, $5.4$, and $6.0$. The 1D power spectrum of the Ly$\alpha$ forest is sensitive to the parameters governing the thermal state of the IGM. Pressure smoothing, $\lambda_{\rm P}$, damps out small-scale fluctuations in the gas, while random thermal motions (sensitive to temperature, or $T_{0}$ and $\gamma$) Doppler broadens Ly$\alpha$ forest lines, further reducing the amount of small-scale structure. Both of these effects combine to produce a prominent small-scale (high-k) cutoff in the flux power spectrum $P(k)$ (Zaldarriaga et al. 2001; Peeples et al. 2010; Rorai et al. 2013; Nasir et al. 2016). Therefore, by carefully studying this cutoff, we expect to be able not only to constrain the thermal state of the gas, but also its full thermal history.

We have created Ly$\alpha$ forest spectra from the simulation computing the H i optical depth at a fixed redshift, which can then be easily converted into a transmitted flux fraction, $F_{\textrm{H}\,\textsc{i}}=e^{−\tau_{\textrm{H}\,\textsc{i}}}$. We refer to Lukić et al. (2015) for specific details of these calculations. We computed the power spectrum, $P(k)$, of the fractional contrast, $\delta F$, at each redshift defined as $\delta F$=$F/\langle F\rangle-1$. A total of $1024^{2}$ skewers were used at each redshift. We computed the power spectrum of each skewer and then calculated the average value at each mode, $k$.

The overall level and precise shape of the fractional contrast power is still sensitive to the mean flux because it changes the density-flux mapping. This is to say, a lower mean flux will shift Ly$\alpha$ observations sensitivity toward lower densities. For this reason, when computing the 1D power spectrum of the fractional contrast from a simulation, it is still very important which mean flux was considered. Following the standard approach, we rescaled the mean flux of each simulation to match a fixed mean flux value. Of course, this rescaling does not affect measurements since one directly measures a flux contrast field. While the mean flux value has been precisely measured at lower redshift, the measurements at $z\gtrsim 5$ are more uncertain. At $z=5.0$ and $5.4$, the current best measurements for the mean flux are the binned values computed by D’Aloisio et al. (2016) from the Becker et al. (2015) high-redshift quasar opacity measurements. These are $\langle F\rangle=0.14\pm 0.01$ for $z=5$, and $\langle F\rangle=0.08\pm 0.006$ for $z=5.4$. These values are consistent with the analytic formula presented by Viel et al. (2013a) derived from their own quasar sample, which are $\langle F\rangle=0.14603,0.071$ respectively. The Fan et al. (2006) measurements of the mean flux using a sample of high-$z$ quasars discovered in the Sloan Digital Sky Survey were $\langle F\rangle=0.1224\pm 0.03$ at $z=5.025$ and $\langle F\rangle=0.074^{+0.03}_{-0.06}$ at $z=5.450$. These values are also in good agreement with the current best measurements considering their larger errors. In fact, the global fit suggested by Fan et al. (2006) based on their own measurements gives $\langle F\rangle=0.1659,0.071$, respectively at these redshifts. However, since we will compare our models with the Viel et al. (2013a) observations of the 1D flux power spectrum in order to give a better qualitative idea of the results found in this work, we have considered the mean mean flux values that give a better overall normalization to these observations and assumed $\sim 7.5$% relative measurement error, which reflects the quoted errors in the results above: $\langle F\rangle=0.16\pm 0.01$ for $z=5$, and $\langle F\rangle=0.055\pm 0.004$ for $z=5.4$. While these values seem to be slightly far from the current best observations by D’Aloisio et al. (2016) at $z=5.4$, they are within the $1\sigma$ C.L. found by Viel et al. (2013a, see Table II) when they performed a marginalized fit of the 1D flux power spectrum for several parameters that included the mean flux ($\langle F\rangle=0.148^{+0.024}_{-0.007}$ at $z=5$, and $\langle F\rangle=0.045^{+0.02}_{-0.001}$ at $z=5.4$)555This is also in agreement with the marginalized fit made by Garzilli et al. (2015, see Table I) to the same dataset. These authors found $\langle F\rangle=0.142^{+0.023}_{-0.017}$ at $z=5$ and $\langle F\rangle=0.054^{+0.014}_{-0.01}$ at $z=5.4$.. Although the lower mean flux measurements from the 1D flux power spectrum could just indicate some fluctuation due to the small number of quasars used to compute the power spectrum at high-$z$, it definitely highlights the relevance of taking into account the mean flux degeneracy when any astrophysical or cosmological information from the 1D flux power spectrum is to be extracted. We return to this issue in our discussion of different degeneracies of the 1D flux power spectrum in Section 5.

The two panels of Figure 3 show the simulated dimensionless 1D flux power spectrum, $kP(k)/\pi$, computed at $z=5.0$ and $5.4$ for the models where we changed the H i reionization history but kept the total heat input constant: EarlyR, MiddleR, LateR, MiddleR-fast and VeryEarlyR (see Table 1). We first note in the two panels that the overall power level increases with redshift, which reflects the fact that as the average mean flux decreases toward higher $z$, density fluctuations are exponentially amplified (e.g. Viel et al. 2004; Palanque-Delabrouille et al. 2013; Viel et al. 2013a). Note also that at all redshifts the difference between models for the low-$k$ modes is very small and therefore the use of high-resolution spectra probing to $k\sim 0.1$ s km*-1* is key.666Detailed studies of the impact of metal absorption features in current state-of-the-art high-resolution spectra have shown that it can increase the 1D flux power spectrum at k$>0.1$ s km*-1* (McDonald et al. 2000, 2005; Lidz et al. 2010; Viel et al. 2013a). For this reason, this is typically the maximum k considered in high-resolution power spectra studies. The models separate at high-$k$ because their disparate reionization histories result in different levels of pressure-smoothing and thermal broadening (see Fig. 2), changing the shape of the small-scale (high-$k$) cutoffs in the power spectra. The color bands for each model show the variation in the 1D flux power spectrum due to the one-sigma uncertainty in the mean flux value at the corresponding redshift. Note that currently these errors seems to translate into a $\sim 10\%$ scatter in the 1D flux power spectrum, which is smaller than the current error bars on the $z=5.4$ measurements (Viel et al. 2013b).

The differences between these models at $z\gtrsim 5$ are particularly interesting because they result primarily from differences in the pressure-smoothing scale, $\lambda_{\rm P}$, as the other parameters governing the thermal state of the IGM, $\gamma$ and $T_{0}$, are very similar (see left panel of Figure 2). These models share exactly the same photoionization and photoheating rates at the redshifts considered, and differ solely in the timing of H i reionization heat injection. These results highlight that the 1D Ly$\alpha$ forest power spectrum is sensitive to the details of H i reionization history even at lower redshifts due to the different pressure-smoothing scale. Any attempts to derive astrophysical or cosmological parameters using high-$z$ Ly$\alpha$ forest observations that do not take this issue into account could obtain biased results (see also Puchwein et al. 2015; Oñorbe et al. 2017).

Figure 4 show the simulated dimensionless 1D flux power spectrum for simulations in which the timing of H i reionization is identical, but which have different amounts of heat injection, $\Delta T$: Middle-cold ($\Delta T=1\times 10^{4}$ K), MiddleR ($\Delta T=2\times 10^{4}$ K), MiddleR-warm ($\Delta T=3\times 10^{4}$ K), MiddleR-hot ($\Delta T=4\times 10^{4}$ K). As expected, the power spectrum shows a larger small-scale cutoff (i.e. toward lower $k$) for simulations with a higher heat input during reionization. It is clear that the effect in the 1D flux power spectrum of a high heat injection during H i reionization (MiddleR-hot) is degenerate with a reionization model with a lower heat input but that completes at higher redshift. Both physical processes produce a higher pressure-smoothing scale at lower redshift.

At higher redshifts, $z\gtrsim 5$ the differences in the power spectrum between simulations shown in Figure 4 are due not only to the effect of the pressure-smoothing scale, $\lambda_{\rm P}$, but also due to the differences in the other thermal parameters, $\gamma$ and $T_{0}$. This is because at these redshifts the IGM in these models are still reaching the asymptotic temperature-density relation after H i reionization (see Figure 2). This highlights the other physical process affecting the Ly$\alpha$ forest lines which is the thermal broadening along the line of sight that also affects the cutoff, in the 1D flux power spectrum. In fact, it is relevant to point out that the differences between models with different $T(\Delta_{\bigstar})$ produce larger differences at $k<0.04$ s km*-1* in the 1D flux power spectrum than thermal models that just differ in the pressure-smoothing scale, $\lambda_{\rm P}$ (see Figure 4). This could open a possibility to distinguish between both physical effects using different $k$-mode ranges of the 1D flux power spectrum, provided that H i reionization happens at low enough redshift to still see these effects.

4.1. Comparison with Observations

Viel et al. (2013a) made measurements of the 1D Ly$\alpha$ forest flux power spectrum at $z=5.0$ and $z=5.4$, using a sample of 25 high-resolution quasar spectra. The redshift bins had width $dz=0.4$ and contained data from $\sim 10$ quasars per bin. In Figure 3 and Figure 4 we also compare the measurements of Viel et al. (2013a, black circles) to our simulation results at the same redshift bins. From these figures it is clear that these measurements already have sufficient precision to begin distinguishing between different reionization models, once the degeneracy due to the mean flux has been taken into account. In what follows we report a first qualitative comparison of this data set with our simulations. A detailed quantitative analysis of these observations using a larger hydrodynamical grid of high-resolution large-volume simulations that expands the full parameter space of the thermal parameters and takes into account relevant degeneracies will be presented in another paper.

We first note from this comparison is that in the context of our current models and with the caveat that the mean fluxes have been chosen to best fit the power at each redshift, the $z=5$ and $z=5.4$ measurements appear to be in agreement with our fiducial model (MiddleR, green line) that uses the Planck $\tau_{e}$ value and $\Delta T=2\times 10^{4}$ K. This picture is consistent with the conventional wisdom that galaxies reionized hydrogen (Robertson et al. 2015). Reionization models driven by active galactic nuclei (AGN; see e.g. Chardin et al. 2015; Madau & Haardt 2015; Khaire et al. 2016) have recently gained traction in light of the discovery of an abundant population of faint AGN at high redshift $z\sim 4-6$ (Giallongo et al. 2015). Such models have higher photoelectric heating of H i and would also doubly ionize helium at these high redshifts (McQuinn 2012) increasing the amount of heat injection in addition to the one associated with H i reionization and moving the cutoff of the 1D flux power spectrum to lower $k$ modes.777In the simulations discussed in this work, we do not consider any high-$z$ He ii reionization model, but this could be easily added using the same formalism applied to H i reionization (see Oñorbe et al. 2017). In any case, current observations of the 1D flux power spectrum at $z=5$ and $z=5.4$ disfavor high-redshift $z\gtrsim 9$ reionization models, far away from Planck constraints on $\tau_{e}$, even for standard galaxy-driven heat injection. We have shown that this is due to the dependence of the 1D flux power spectrum cutoff on the timing of reionization because the pressure-smoothing scale retains memory of the thermal history.

4.2. Prospects for Measuring the Power Spectrum at $z\simeq 6$

Motivated by the results in the previous subsection for the 1D flux power spectrum at $z=5.0,5.4$ and the increasing number of high-$z$ quasars that are uncovered by recent surveys (e.g. Bañados et al. 2014; Matsuoka et al. 2016; Bañados et al. 2016), we also wish to explore the 1D flux power spectrum at $z=6$, and study the feasibility of making a power spectrum measurement at this redshift with current facilities. Whereas, as discussed previously, the power spectrum signal increases toward increasing redshift, the mean flux also begins to drop precipitously, lowering the signal-to-noise (${S/N}$) ratio level of the quasar spectra, thus increasing the importance of noise in the power spectrum measurement. At $z=6.0$ we have assumed a mean flux value of $\langle F\rangle=0.011$ ($\langle\tau_{\rm eff}\rangle=4.5$), which is consistent with the latest measurements of the effective optical depth at this redshift (Becker et al. 2015; D’Aloisio et al. 2016). We have also assumed that the $1\sigma$ error on the mean flux at this redshift will be $10\%$ of its value, which is a reasonable assumption given the large numbers of $z\sim 6$ quasars recently discovered, $\sim 150$.

To this end, we computed mock observations for the MiddleR simulation assuming high-resolution spectra (${S/N}$ of 35 per resolution element of $8$ km s*-1*) of 10 quasars at $z=6.3$ employing a path length equivalent of $\Delta z=0.5$ per quasar. For this we calculate the corresponding path length per quasar in cMpc at this redshift. We created random samples for the 10 quasars from the simulations and added noise realizations to each skewer. We then computed the mean power spectrum and subtracted off the average noise level. Results of one of these mock observations are shown as black squares in Figure 5 along side with the error bars computed from a set of 50 mock observations. At the $S/N$ considered, the measurement is dominated by cosmic variance. The noise in the quasar spectra is very significant at this redshift, but the increase in overall power due to the decrease of the mean flux still allows us to measure of the power using a sample of 10 quasars.

The two panels of Figure 5 show the 1D flux power spectrum for the same simulations as discussed in Figure 3 (left panel) and Figure 4 (right panel), but now at $z=6$. The color bands for each model indicate the variation in the 1D flux power spectrum due to $1\sigma$ changes in the mean flux. We can see that the overall scale of the flux power has increased compared to the values at lower redshift due to the decrease in mean flux. As discussed above, the power increases with redshift as the mean flux goes down because this amplifies the fluctuations. Note that the differences between the models have increased relative to $z\sim 5$, especially for high-$k$ modes and for simulations where the heat injection is varied. This is because at $z=6$, we are closer to reionization and therefore not only $\lambda_{\rm P}$, but also the other thermal parameters, $T_{0}$ and $\gamma$ (see Figure 2), are still affected by the details on how H i reionization happened (McQuinn & Upton Sanderbeck 2016).

Data with a size and $S/N$ comparable to our assumed mock are clearly within reach. For example, about $5$ such quasar spectra already exist in public telescope archives (Becker et al. 2015), so a sample of $10$ would be accessible with modest allocations of an $8$ m class telescope time. We have shown that the differences between models are greater at $z=6$ than at lower redshift, compensating for the possibly lower $S/N$. In order to understand how reionization heated the IGM and constrain both the reionization history and the heat injected, we therefore need to push as far back into their reionization epoch as possible, where not only $\lambda_{\rm P}$, but also $\gamma$ and $T_{0}$ still could have memory of reionization.

5. Discussion

5.1. Convergence of the Results

While the resolution and method employed in our Ly$\alpha$ simulations are currently state-of-the-art for this type of analysis (Lukić et al. 2015; Oñorbe et al. 2017), several aspects of our simulations call for caution. First, one should consider the level of convergence of the 1D flux power spectrum at these high redshifts for our grid of simulations that have $L_{\rm box}=20$ Mpc h*-1* and $N=1024^{3}$. Lukić et al. (2015) reported a careful convergence analysis of Nyx simulations of the Ly$\alpha$ forest similar to those employed in this work, but they only explored $2\leq z\leq 4$. At $z=4$, the 1D flux power spectrum of their simulations were converged to $<5\%$ in terms of spatial resolution (missing power in large k modes) but $<8\%$ due to box size effects (overestimated power at large $k$ modes). We show in the Appendix A resolution and box size convergence results at $z=5$ and $z=6$ that reach similar conclusions, although approaching $\lesssim 10\%$ convergence level at the $k\sim 0.1$ s km*-1* modes, as typically used for power spectrum measurements, and much better as we move to lower $k$ modes (larger scales). In any case, although convergence issues have to be taken into account, they do not seem to change the general conclusions of this work. Similar results at these redshifts but for simulations using the Gadget code can be found in Bolton et al. (2009) and Bolton et al. (2017).

Another relevant issue is that our simulations do not model galaxy formation. Therefore we neglect any impact that the effect of stellar or black hole feedback could have on the IGM (Theuns et al. 2002a; Kollmeier et al. 2006; Desjacques et al. 2006; Tepper-García et al. 2012). Viel et al. (2013b) showed that at $z\sim 3$ this effect could lead to differences of $\sim 10\%$ in the 1D flux power spectrum. This is a very relevant issue that should be explored in more detail with current state-of-the-art feedback models. However, it is expected that the effects of feedback on the IGM should only decrease at higher redshifts for two reasons: first, the Ly$\alpha$ forest is tracing lower gas densities at average locations of the universe (see, e.g., Fig. 7 in Lukić et al. 2015), and the smaller number of galaxies or black holes at these redshifts make it difficult for any feedback to alter the thermal state of the IGM.

Recent measurements of the Ly$\alpha$ optical depth at high redshift have found enhanced scatter at $z>5.5$ that exceeds what can be attributed to density fluctuations alone (Fan et al. 2006; Becker et al. 2015). It has been argued that they are driven by fluctuations in the radiation field (Davies et al. 2016; D’Aloisio et al. 2016), or the temperature field (D’Aloisio et al. 2015), both of which may be inevitable byproducts of a patchy, extended, and late-ending reionization process. Still others have interpreted these fluctuations as evidence that reionization was actually driven by rare AGN (Chardin et al. 2015; Madau & Haardt 2015; Khaire et al. 2016). In any case, the possible effect of UVB or temperature fluctuations is currently neglected in standard optically thin simulations like those we used in this work. Previous studies have shown that these effects are manifest on much larger scales, $\gtrsim 10$ Mpc h*-1*, that is,so much lower $k$ modes than those most sensitive to the thermal state of the IGM (e.g. Abel & Haehnelt 1999; Meiksin & White 2004; Cen et al. 2009; Pontzen 2014; Gontcho A Gontcho et al. 2014; Malloy & Lidz 2015; D’Aloisio et al. 2016). For example, D’Aloisio et al. (2016) computed semi-numerical models with and without temperature fluctuations and showed that they did not generate small-scale power. Considering both temperature and UVB fluctuations is clearly the direction in which the modeling needs to move forward. We plan to study this in detail in the near future with self-consistent hydrodynamical simulations.

5.2. Degeneracy with Cosmological Parameters and Warm Dark Matter

The 1D flux power spectrum depends not only on the thermal parameters of the IGM, but also on cosmological parameters. Here we have focused our analysis on one specific cosmological model, but we have checked that when one considers the range of parameters allowed by Planck Collaboration et al. (2016a) the differences at the $k$ modes studied in this work are never larger than $8\%$ and the cutoff is unaffected (see Oñorbe et al. 2017, for more details on these models).

While changes of cosmological parameters within $1\sigma$ of the Planck constraints result in minimal changes of the $z\sim 5-6$ power spectrum relative to the parameters governing reionization and its heating effect, this is not true when one considers dark matter particle properties, such as warm dark matter. Their free-streaming horizon leads to a suppression of the small-scale power and therefore to a degenerate effect on the power spectrum with the thermal parameters (Viel et al. 2013a; Garzilli et al. 2015). This is equivalent to a 3D smoothing at high redshift that continues to decrease as nonlinearities increase. Therefore the small-scale cutoff in warm dark matter models moves from very low k modes (large scales) to higher k modes (smaller scale) at progressively lower redshifts. This happens until the IGM becomes hotter and the IGM temperature (i.e. the $T_{0}$ and $\gamma$) determines the position of the cutoff. In order to provide reliable WDM constraints, it is therefore essential to marginalize out reionization nuisance parameters. Garzilli et al. (2015) highlighted the degeneracy between different warm dark matter masses and the temperature of the IGM. However, these authors did not discuss the extra degeneracy resulting from the unknown redshift of reionization and its associated heat injection, as we demonstrate in this work. This could be relevant given that the fiducial value in their models is $z_{\rm reion}=12$ and the lowest reionization redshift considered in their grid, by including one simulation, was $z_{\rm reion}=8$. From their Bayesian analysis, Viel et al. (2013b) and Garzilli et al. (2015) obtained very low temperature values at $5<z<5.4$ (see Figure 2). Although this is not in complete disagreement with our qualitative comparison of the measurements to our models, we caution that their simulations use a total of $2\times 512^{3}$ dark matter and gas particles within a periodic box of $L_{\rm box}=20$ cMpc h*-1888Note that Viel et al. (2013b) use 4 N=$2\times 768^{3}$ and $L_{\rm box}=20$ cMpc h-1* simulations to correct for resolution convergence.. Based on the resolution convergence tests presented in the Bolton et al. (2017, see their Figure A4) spectrum using the same code, these simulations could underestimate the power with $\sim 20\%$ error at the k modes most relevant to study the cut off of the power spectrum ($0.05\lesssim$ k $\lesssim 0.1$ s/km). This could produce an artificial increase of the cutoff just due to resolution.

5.3. Comparison to Previous Work

Nasir et al. (2016) used hydrodynamical simulations with a total of $2\times 512^{3}$ dark matter and gas particles within a periodic box of $L_{\rm box}=10$ cMpc h*-1* to discuss the possibility of constraining the thermal history of the IGM during H i reionization by studying the 1D flux power spectrum at $z=5$ from cosmological hydrodynamical simulations. To simulate different reionization histories, they adopted an approach different from ours. Namely they applied a multiplying factor to the Haardt & Madau (2001) photoheating rates ($A\times\dot{q}\times\Delta^{B}$, where $\Delta$ is the specific overdensity of that cell). They also used simple cutoffs of the Haardt & Madau (2001) UVB rates at various redshifts to model different reionization timing. Using these simulations, they showed the effect on the 1D flux power spectrum of different thermal histories in which they only changed the timing of H i reionization and tried to study its degeneracy with the temperature of the IGM, $T_{0}$. Interestingly, these authors found that this degeneracy between the timing of reionization and the temperature can in fact be broken using different scales of the 1D flux power spectrum. They find that the $z=5$ 1D flux power spectrum is more sensitive to the timing at $0.03<k<0.13$ s km*-1* scales. Our analysis also shows that these scales are the most sensitive to the timing of reionization (see Figures 3), but it also indicates that information about the lower $k$ modes will be crucial in order to break the degeneracy between the timing of reionization and other parameters, the temperature at mean density, $T_{0}$ or the mean flux at that specific redshift.

Nasir et al. (2016) attempt to quantify the effect of a different reionization timing using a new parameter, $u_{0}$, which is an integral in time of the heating rate per proton mass of the simulation at mean density.999It is defined as $u_{0}(z)=\int_{z}^{z_{\rm reion}}\frac{\sum_{i}n_{i}\dot{q}_{i}}{\bar{\rho}}\frac{dz}{H(z)(1+z)}$, where $\bar{\rho}$ is the mean background baryon density and $n_{i}$ and $\dot{q}$ stand for the density and photoheating rates of the following species: $i=[\textrm{H}\,\textsc{i},\textrm{He}\,\textsc{i},\textrm{He}\,\textsc{ii}]$. With this parameter, the authors try to quantify the pressure-smoothing scale in their models. In order to facilitate comparison with their work we have computed the value of $u_{0}$ at $z=4.9$ for our simulations and included them in Table 1. However, we caution about using this parametrization as $u_{0}$ measures the heating just at mean density. Therefore this parametrization will be valid as long as the the reionization models do not have density-dependent heating. Models with heating rates that depend on the density but normalized at mean density (e.g. Becker et al. 2011; Becker & Bolton 2013) share the same $u_{0}$ value, but have a different power spectrum at the $k$ modes that are more sensitive to the thermal history, $k>0.03$ s km*-1*, as they in fact have different pressure-smoothing scales.

6. Conclusions

In this paper we have used state-of-the-art hydrodynamical simulations that allowed us to self-consistently model different reionization models. We present an ensemble of simulations consistent with the latest measurements of the Thompson-scattering optical depth, $\tau_{\rm e}$ recently reported by Planck (Planck Collaboration et al. 2016b, c). These models are defined by when reionization happened, $z_{\rm reion}$, and how much heat was injected into the IGM during reionization, $\Delta T$. Our simulations show that although by $z\sim 6$ the temperature of IGM gas has mostly forgotten about reionization heat injection, the pressure-smoothing scale at these redshifts depends sensitively on how and when reionization occurred. This is because both the cooling and dynamical times in the rarefied IGM are long, comparable to the Hubble time, and therefore memory of H i reionization is retained (Rorai et al. 2013; Kulkarni et al. 2015; Oñorbe et al. 2017). We have found a degeneracy in the pressure-smoothing scale at $z<6$ between when reionization occurred and the amount of heat injected during reionization. For a fixed reionization history, the pressure-smoothing scale increases as we increase the heat injection. Similarly, the pressure-smoothing scale increases with the redshift of reionization in models with a fixed amount of heat injection.

In order to investigate the effects of these different thermal histories on the properties of the Ly$\alpha$ forest, we compute the Ly$\alpha$ 1D flux power spectrum at $z\sim 5-6$ for our simulation ensemble. Pressure smoothing damps out small-scale fluctuations in the IGM, while thermal vibrations of IGM gas Doppler broadens Ly$\alpha$ forest lines, further reducing the amount of small-scale structure. Both of these effects combine to produce a prominent small-scale (high-k) cutoff in the Ly$\alpha$ 1D flux power spectrum (Zaldarriaga et al. 2001; Peeples et al. 2010). We have found that at these high redshifts, the 1D flux power spectrum is especially sensitive to the pressure-smoothing scale of the IGM and not only its temperature. Therefore extant thermal signatures from reionization can be detected by analyzing the Ly$\alpha$ forest power spectrum at these redshifts.

We have also conducted a first qualitative comparison of the 1D flux power spectrum measurements at $z=5-5.4$ made by (Viel et al. 2013a) with our simulation ensemble. Taking Planck constraints on reionization at face value, we have shown that models with a fiducial heat input during H i reionization consistent with standard galaxy-driven reionization models are sufficient to explain the observations. We work on a more complete analysis of this in the near future, with a larger simulation grid, marginalizing out all the different relevant parameters, including the mean flux, and improving upon the reionization modeling.

We have also presented a feasibility study of performing a similar measurement at $z=6$ creating mock observations that assumed a realistic sample of quasars at this redshift both in terms of sample size and $S/N$. We found that combining 10 quasars should be enough to distinguish between different thermal histories of the IGM. Our results indicate that quasar spectra at high redshift can not only be useful to constrain when reionization happened via the study of Ly$\alpha$ opacity measurements (e.g. Fan et al. 2006; Becker et al. 2015), but also to understand the thermal history of the universe. This is the small-scale structure measured from high-resolution spectra can be used to understand the thermal history of the universe, further constraining the timing and heat injection by reionization. Taking into account that there are only a few direct observational probes of reionization currently available, we think that it is important to push in this direction in the near future.

In this regard, pushing these measurements at this and higher redshifts ($z\sim 6$) will be crucial to improve the power of the Ly$\alpha$ forest to constrain H i reionization. The new recent fivefold increase in the number of bright quasars at $z>5$ is obtained from deep wide-field optical/IR surveys like CFHQS (Willott et al. 2010), dark energy survey (DES, Reed et al. 2015), ESO public surveys (VST/KiDS and VISTA/VIKING Venemans et al. 2015a), and Pan-STARRS1 (Bañados et al. 2014; Venemans et al. 2015b). Currently, the total number $z>5.5$ quasars available for study is $\sim 173$, so if all could be used to study the cutoff of the Ly$\alpha$ forest power spectrum would reduce the errors by a factor of $\sim 4$ compared with current measurements. Moreover, large datasets of medium- and high-resolution Ly$\alpha$ spectra are already available at lower redshifts ($z\lesssim 4$, e.g. O’Meara et al. 2015; Iršič et al. 2017) that combined with the large redshift data can provide a more comprehensive constraint of reionization. Therefore this requires starting to focus now on improving the theory to exploit this increased precision. While modeling the 1D flux power spectrum with hydrodynamical simulations at sufficiently high accuracy is an incredible computational challenge, the advent of high-performance computing power and the high scalability of Nyx has allowed us to significantly improve the accuracy of our predictions in recent years.

Finally, given that there are few observables that are sensitive to the thermal state of baryons at the earliest redshifts, the 1D flux power spectrum at $z\gtrsim 5$ offers a unique opportunity to explore not only H i reionization, but also constrain other physical scenarios that alter the thermal history of the IGM, for example, models that alter the thermal state of the IGM via X-ray pre-heating coming from starburst galaxies, supernova remnants or miniquasars (Oh 2001; Glover & Brand 2003; Madau et al. 2004; Furlanetto 2006; Madau & Fragos 2017), dark matter annihilation or decay (Liu et al. 2016), cosmic rays (Samui et al. 2005), blazar heating (Chang et al. 2012; Puchwein et al. 2012), broadband intergalactic dust absorption (Inoue & Kamaya 2008) or high-$z$ exotic reionization scenarios driven by Population III stars (e.g. Manrique et al. 2015).

We thank M. White, J. Miralda-Escudé and the members of the ENIGMA group at the Max Planck Institute for Astronomy (MPIA) for helpful discussions. ZL was partially supported by the Scientific Discovery through Advanced Computing (SciDAC) program funded by U.S. Department of Energy Office of Advanced Scientific Computing Research and the Office of High Energy Physics. Calculations presented in this paper used the hydra cluster of the Max Planck Computing and Data Facility (MPCDF, formerly known as RZG) MPCDF is a competence center of the Max Planck Society located in Garching (Germany). We also used resources of the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. This work made extensive use of the NASA Astrophysics Data System and of the astro-ph preprint archive at arXiv.org. This work is based on observations obtained with Planck (http://www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada.

Appendix A Numerical Convergence

In this section we discuss the convergence tests with spatial resolution and box size for the two fundamental quantities studied in this paper, the mean optical depth and the 1D flux power spectrum. To study the effects of spatial resolution, we have run 4 simulations with the same thermal history (e.g., same UVB; MiddleR) and box size of $L_{\rm box}=10$ Mpc h*-1* but increasing numbers of resolution elements: $128^{3}$ (dot-dashed green), $256^{3}$ (dotted blue), $512^{3}$ (dashed red), and $1024^{3}$ (black). To simplify the comparison, simulations performed in the same box size share the same large-scale modes, the only difference being that higher resolution runs have more modes sampled on small scales. These simulations have a cell size of $78$, $39$, $20$, and $10$ kpc h*-1*, respectively, and therefore the $512^{3}$ run has the same spatial resolution as the simulations discussed in this work ($L_{\rm box}=20$ Mpc h*-1* and $1024^{3}$ cells). We also ran one more simulation with the same thermal history (MiddleR), a box size of $L_{\rm box}=40$ Mpc h*-1* and $2048^{3}$ cells in order to study box size effects. We compare this simulation with two other runs with the same spatial resolution, but decreasing box size: the $L_{\rm box}=20$ Mpc h*-1* - $N_{\rm cell}=1024^{3}$ run, which correspond to the simulations used in this work, and the $L_{\rm box}=10$ Mpc h*-1* - $N_{\rm cell}=512^{3}$ run also used in the spatial resolution study.

The evolution of the mean optical depth $\langle\tau_{\rm eff}\rangle$ for all these simulations is shown in Figure 6. We computed its evolution directly from the simulation mean flux, $\langle\tau_{\rm eff}\rangle=-\ln\langle F\rangle$, without any rescaling of the photoionization rate. Thus, all simulations use exactly the same photoionization rates at all redshifts. The left panel shows the convergence of the mean optical depth as we increase the spatial resolution, while the right panel shows the convergence for different box sizes. The simulations discussed in this work (red dashed lines in both panels) show a convergence level below $<5\%$ between $4\leqslant z\leqslant 6$ both in terms of spatial resolution and box size.

Figure 7 shows the convergence tests of the 1D flux power spectrum at redshift $z=5$ (left column) and $z=6$ (right column) for the same simulations. For this test, we rescaled the mean flux of all the simulations to the same value. We used the fit between mean flux and redshift suggested in Oñorbe et al. (2017), obtained using a wide range of data sets between $0<z<6$, but the exact values employed do not change our conclusions. For the resolution tests (upper row) and the box size tests (lower row), we find a $\lesssim 10\%$ level of convergence for $k$ modes lower than $\sim 0.04$ s km*-1* in the simulations with the same resolution and box size used in this work (red dashed lines). The error in these modes is mainly driven by box size effects as the resolution tests show a better convergence. However, for modes more relevant to study the thermal cutoff ($0.04<k<0.1$ s km*-1*) we find a $10\%$ convergence level at $z=5$ and $15\%$ at $z=6$ that is mainly driven by spatial resolution effects. Note that this quoted convergence level is the worst level at the highest $k$ mode, $0.1$ s km*-1*, but it decreases as we move to lower $k$ values. It is also very relevant to indicate that spatial resolution effects can only move the cutoff of the 1D flux power spectrum to higher $k$ values as we increase the resolution. Similar results at these redshifts but for simulations using the Gadget code can be found in Bolton & Becker (2009) and Bolton et al. (2017). Convergence results at lower redshifts for the same code used in this paper along with a more detailed discussion can be found in Lukić et al. (2015).