Self-consistent Modeling of Reionization in Cosmological Hydrodynamical Simulations

In our current standard model of the universe, hydrogen and helium account for 99% of the baryonic mass density (Planck Collaboration et al. 2016a). After the recombination epoch, these elements remain neutral until ultraviolet radiation from star-forming galaxies and active galactic nuclei reionizes them. Therefore, this ultraviolet background (UVB) governs the ionization state of intergalactic gas and plays a key role in its thermal evolution through photoheating. During the reionization of ${\rm{H}}\,{\rm{I}}$ and later $\mathrm{He}\,{\rm{II}}$, ionization fronts propagate supersonically through the intergalactic medium (IGM), impulsively heating gas to ∼10⁴ K (see, e.g., Abel & Haehnelt 1999; McQuinn 2012; Davies et al. 2016). As the universe evolves, it is well known that the balance between cooling due to Hubble expansion and inverse-Compton scattering of cosmic microwave background (CMB) photons and heating due to the gravitational collapse and photoionization heating give rise to a well-defined temperature–density relationship in the IGM (Hui & Gnedin 1997; McQuinn 2012):

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

where ${\rm{\Delta }}=\rho /\bar{\rho }$ is the overdensity with respect to the mean and T₀ is the temperature at the mean density. Immediately after the reionization of ${\rm{H}}\,{\rm{I}}$ ($z\lesssim 6$) or $\mathrm{He}\,{\rm{II}}$ ($z\lesssim 3$), T₀ is likely to be around $\sim 2\times {10}^{4}$ K and $\gamma \sim 1$ (Bolton et al. 2009; McQuinn et al. 2009); at lower redshifts, T₀ decreases as the universe expands, while γ is expected to increase and asymptotically approach a value of 1.62 (Hui & Gnedin 1997).

Another important physical ingredient to describe the thermal state of the IGM is the gas pressure support. At small scales and high densities, baryons experience pressure forces that prevent them from tracing the collisionless dark matter. This pressure results in an effective 3D smoothing of the baryon distribution relative to the dark matter, at a characteristic scale. known as the Jeans pressure smoothing scale, ${\lambda }_{{\rm{P}}}$. In an expanding universe with an evolving thermal state, this scale at a given epoch is expected to depend on the entire thermal history, because fluctuations at earlier times expand or fail to collapse depending on the IGM temperature at that epoch (Gnedin & Hui 1998; Kulkarni et al. 2015). Recently, Rorai et al. (2013) and A. Rorai et al. (2015, in preparation) have shown that an independent measurement of the pressure smoothing scale can be obtained using the coherence of Lyα forest absorption in close quasar pairs (Hennawi et al. 2006, 2010).

Lyα forest observations between $2\lt z\lt 6$ probe the moderate overdensities characteristic of the IGM and therefore are a crucial tool to understand the properties of the UVB. In the past decade, the precision of these measurements has continued to grow both in terms of their numbers (BOSS³ survey) and in quality (high signal-to-noise ratio spectrum from, e.g., O'Meara et al. 2015). However, while it seems that we keep learning more and more about the ionization history of the universe, for both ${\rm{H}}\,{\rm{I}}$ and $\mathrm{He}\,{\rm{II}}$ reionizations (e.g., Becker & Bolton 2013; Syphers & Shull 2014; Worseck et al. 2014; Becker et al. 2015) the thermal history of the universe is still far from certain. The statistical properties of the Lyα forest are sensitive to the thermal state of the gas, through both thermal broadening of lines and pressure support. When constraints on the thermal history are reviewed, they yield very puzzling results. Measurements of T₀ from different groups utilizing different methodology are in poor agreement (Schaye et al. 2000; Bolton et al. 2008, 2014; Lidz et al. 2010; Becker et al. 2011; Garzilli et al. 2012; Rudie et al. 2012; Boera et al. 2014). A similar problem appears when measurements of the slope of the temperature–density relation, γ, are compared. At $z\simeq 3$ some authors have even found that γ is either close to isothermal ($\gamma =1$) or even inverted ($\gamma \lt 1$; Bolton et al. 2008; Viel et al. 2009; but see Lee et al. 2015). Most studies of the thermal state of the IGM ignore uncertainties resulting from the unknown pressure smoothing scale (but see Becker et al. 2011; Puchwein et al. 2015), which produces a 3D smoothing that is difficult to disentangle from the similar but 1D smoothing resulting from thermal broadening (Peeples et al. 2010a, 2010b; Rorai et al. 2013). Therefore, ignoring this effect has probably contributed to the confusing and sometimes contradictory published constraints on T₀ and γ (Puchwein et al. 2015).

With the help of accurate models of the IGM, the statistics of the Lyα forest can be used to constrain its thermal parameters and ultimately cosmic reionization. Ideally one will run coupled radiative transfer hydrodynamical simulations that include extra physics governing the sources of ionizing photons (stars, quasars, etc.). Despite significant progress on this front (Gnedin 2014; So et al. 2014; Wise et al. 2014; Norman et al. 2015; Ocvirk et al. 2016; Pawlik et al. 2015), these simulations are still too costly for sensible exploration of the parameter space. For this reason, the dominant approach, implemented in the vast majority of hydrodynamical codes, is to assume that all gas elements are optically thin to ionizing photons, such that their ionization state can be fully described by a uniform and isotropic UV+X-ray background radiation field. Thus, the radiation field is encapsulated by a set of photoionization and photoheating rates that evolve with redshift for each relevant ion. The minimal set of ions is ${\rm{H}}\,{\rm{I}}$, $\mathrm{He}\,{\rm{I}}$, and $\mathrm{He}\,{\rm{II}}$ in order to track the most relevant ionization events, as well as the thermal heating associated with them. Of course, although this optically thin approximation is a valid assumption once the mean free path of ionization photons, ${\lambda }_{\mathrm{mfp},\nu }$, is large enough, it is certainly not true during cosmic reionization events. As such, this optically thin approach is not meant to provide an accurate description of reionization itself, but it should at least provide a reasonable description of the heat injection associated with reionization.

This is important since galaxies forming during the reionization epoch are sensitive to the thermal state of the gas, and even well after reionization gas elements can retain thermal memory of reionization heating (Gnedin & Hui 1998; Kulkarni et al. 2015). It is important to remark here that these UVB models have relevant consequences for galaxy formation and evolution models and hydrodynamical simulations. Several groups have already shown how important the UVB model is to determine the star formation of the first galaxies and their evolution by not only setting the minimum halo mass able to form stars (i.e., halos massive enough to overcome gas pressure forces; Rees 1986; Sobacchi & Mesinger 2013) but also regulating the gas accretion from the IGM into the more massive halos (Quinn et al. 1996; Simpson et al. 2013; Benítez-Llambay et al. 2015; Wheeler et al. 2015).

The standard approach is to adopt photoionization and photoheating rates from semianalytical synthesis models of the UVB (Haardt & Madau 1996, 2001; Faucher-Giguère et al. 2009; Haardt & Madau 2012). However, these UVB synthesis models surely break down during reionization events, and the validity of using them in optically thin simulations (during reionization) is questionable. Moreover, as we will show, these models are fundamentally inconsistent during reionization, leading to different reionization histories in the simulations than the ones given by the authors. Specifically, they reionize the universe too early, and as a result they produce spurious heating of the IGM at early times (see Section 2 and Figure 1). In this paper, we improve on the limitations of current UVB models to provide reliable ionization and thermal histories during reionization by developing a new method to model ionization and heating during reionization in hydrodynamical simulations. In the context of this method, we demonstrate how to run simulations with self-consistent ionization and thermal histories that agree with constraints from the CMB and IGM measurements. Moreover, we make these new tables publicly available in the default format used by most cosmological codes.

Zoom In Zoom Out Reset image size

The outline of the paper is as follows. In Section 2 we discuss in detail current standard methods that include the effect of the UVB in optically thin hydrodynamical simulations. We show that these models have problems reproducing the desired ionization and thermal histories. In Section 3 we present a new method to improve the current models of the UVB during reionization events. The different reionization models considered in this work, based on current observational constraints, are motivated in Section 4. We describe the basic details of the hydrodynamical cosmological code that we have used in this work, the analysis pipeline, and the properties of the simulations in Section 5. The ionization and thermal histories of the simulations using the new UVB models are shown and examined in Section 6. We explore the possibility of reproducing observational constraints on the ${\rm{H}}\,{\rm{I}}$ and $\mathrm{He}\,{\rm{II}}$ transmission in Section 7. In Section 8 we discuss the limitations of our new approach, provide a comparison to previous work, and discuss previous work using incorrect UVB models in galaxy formation simulations that likely overestimate the impact of photoionization feedback. We conclude in Section 9. In Appendix A we provide details on how the new photoionization and photoheating rates are derived in our method. In Appendix B we present the ionization and thermal histories of several widely used UVB models. The effects of cosmology on the new models are discussed in Appendix C. Finally, in Appendix D we present the photoionization and photoheating rates of the new models.

Current UVB models used in optically thin hydrodynamical simulations give the evolution of both photoionization, ${{\rm{\Gamma }}}_{\gamma }$ (${\rm{H}}\,{\rm{I}},\mathrm{He}\,{\rm{I}},\mathrm{He}\,{\rm{II}}$) and photoheating $\dot{q}$ (${\rm{H}}\,{\rm{I}},\mathrm{He}\,{\rm{I}},\mathrm{He}\,{\rm{II}}$) rates with redshift. Therefore, during reionizations the simulated IGM is photoheated everywhere by the same spectrum using a homogeneous value of $\dot{q}$. It is important to remark that the photoheating rates are given as energy per second per ion. Then, for example, in the case of ${\rm{H}}\,{\rm{I}}$ reionization, the photoheating rate per volume for each resolution element is ${n}_{{\rm{H}}{\rm{I}}}{\dot{q}}_{{\rm{H}}{\rm{I}}}$.

Haardt & Madau (1996) were the first to try to develop self-consistent UVB models in a cosmological context using radiative transfer methods and taking into account observations of the ionizing sources (namely, quasars and galaxy luminosity functions) and the absorption of the ionizing photons (column density distribution of neutral hydrogen, ${N}_{{\rm{H}}{\rm{I}}}$, absorbers, and ${\rm{H}}\,{\rm{I}}$ mean flux). An ionization front heats up the gas behind it (see, e.g., Abel & Haehnelt 1999; McQuinn 2012; Davies et al. 2016), and therefore it is crucial to include radiative transfer effects in the UVB modeling. Self-consistent methods use photoionization modeling codes that implement 1D radiative transfer (e.g., cloudy) to try to take this effect into account. Subsequent efforts have developed these models further (Faucher-Giguère et al. 2009; Haardt & Madau 2012). These UVB models are adopted in essentially all nonadiabatic cosmological hydrodynamic simulations to compute the ionization state and photoheating rates of intergalactic gas (Somerville & Davé 2015). These include simulations focusing on the properties of the IGM (e.g., Katz et al. 1996; Miralda-Escudé et al. 1996; Lukić et al. 2015), but also simulations modeling galaxy formation and evolution (e.g., Vogelsberger et al. 2013; Hopkins et al. 2014; Shen et al. 2014; Governato et al. 2015; Davé et al. 2016).

We first want to present the ionization and thermal histories obtained when one of the most widely used UVB models (e.g., tabulated photoionization and photoheating rates) is used (Haardt & Madau 2012, hereafter HM12). The upper panel of Figure 1 shows the ${\rm{H}}\,{\rm{II}}$ (solid black line) and $\mathrm{He}\,{\rm{III}}$ (dot-dashed black line) ionization history calculated by the HM12 model (black lines), which indicates that ${\rm{H}}\,{\rm{I}}$ reionization should finish at ${z}_{\mathrm{reion},{\rm{H}}{\rm{I}}}=6.7$ and $\mathrm{He}\,{\rm{II}}$ at ${z}_{\mathrm{reion},\mathrm{He}{\rm{II}}}=2.75$. These are given by the volume filling fraction evolution, ${Q}_{{\rm{H}}{\rm{II}}}(z)$ which can be thought of as the probability that the hydrogen in a given region is ionized (Madau et al. 1999), and ${Q}_{\mathrm{He}{\rm{III}}}(z)$, the analogous quantity for a doubly ionized helium region. In this panel we also show the ionization history of a hydrodynamical cosmological simulation using the HM12 UVB model (green lines). This simulation uses the standard methodology employed in other optically thin simulations that we describe in detail in Section 5. To obtain the ionization history from the simulations, we computed the ionization fraction of each volume element in the simulation, ${x}_{{\rm{H}}{\rm{II}}}={n}_{{\rm{H}}{\rm{II}}}/{n}_{{\rm{H}}}$ and ${x}_{\mathrm{He}{\rm{III}}}={n}_{\mathrm{He}{\rm{III}}}/{n}_{\mathrm{He}}$, and then calculate the average weighting in volume (i.e., averaging all of the cells).⁴ The first striking thing that we learn from this comparison is that using the HM12 photoionization rates effectively reionizes the universe much earlier than the reionization redshift reported by HM12. It appears that some aspect of the HM12 calculation is not internally consistent.

The middle panel of Figure 1 shows the integrated electron scattering optical depth, defined as

$$\begin{eqnarray}&&\tau e(z)=c{\sigma }_{{\rm{T}}}{\int }_{0}^{z}{n}_{{\rm{e}}}\displaystyle \frac{{(1+z^{\prime} )}^{2}{dz}^{\prime} }{H(z^{\prime} )},\end{eqnarray} \tag{ 2 }$$

where c is the velocity of light, ${\sigma }_{{\rm{T}}}$ is the Thomson cross section, H(z) is the Hubble parameter, and ${n}_{{\rm{e}}}$ is the proper electron density. We have computed the electron density in our simulation as ${n}_{{\rm{e}}}=\langle {n}_{{\rm{H}}}\rangle (1+\chi )\langle {x}_{{\rm{H}}{\rm{II}}}\rangle +\chi \langle {x}_{\mathrm{He}{\rm{III}}}\rangle$, where $\chi ={Y}_{{\rm{p}}}/(4{X}_{{\rm{p}}})$ and ${X}_{{\rm{p}}}$ and ${Y}_{{\rm{p}}}$ are the hydrogen and helium mass abundances, respectively (see Appendix A.1 for a detailed derivation of this equation). We also make the standard assumption that the reionization of $\mathrm{He}\,{\rm{I}}$ is perfectly coupled with that of ${\rm{H}}\,{\rm{I}}$. The observational constraints on ${\tau }_{{\rm{e}}}$ coming from the CMB (${\tau }_{{\rm{e}}}=0.078\pm 0.019$; Planck Collaboration et al. 2016a) are indicated by a gray band.⁵ The results of the simulation (green) not only differ from the expected results of the model (black) but also are in strong disagreement with the observational constraints.

Moreover, the lower panel of Figure 1 shows the thermal history of this simulation (green line) via the evolution of the temperature at mean density, T₀ defined in Equation (1). See Section 5.1 for a discussion of the procedure used to fit the ρ–T relation. We can see not only that reionization events occur too early, but also that the heating associated with them starts at much earlier times, $z\sim 15$. We have confirmed that this result is not due to resolution, assumed cosmological model, atomic rates considered, or some particular code characteristic (see Figures 3 and A1 in Puchwein et al. 2015, for the same effect but using a SPH Lagrangian code, GADGET-3, and using both an equilibrium and nonequilibrium ionization solver). Thus, why do the reionization histories in the simulation and the one calculated by these authors differ so much?

UVB semianalytical synthesis models—such as the one used by HM12—rely on two main assumptions: (1) the photoionizing background is everywhere uniform, and (2) radiative transfer is optically thin. Clearly, both assumptions break down during reionization. While a full solution requires a radiative transfer simulation, a large number of IGM studies are insensitive to the reionization details. Thus, the goal is to nevertheless have an approximate UVB background representing the mean UVB to adopt in an optically thin simulation. Therefore, it is important to state upfront that the rates obtained under these sets of assumptions and the ones obtained in a patchy reionization model (e.g., a radiation transfer simulation) can differ significantly during reionization. Lidz et al. (2007) illustrate this point in a simple way by considering two toy models. In the first case, representing patchy reionization, we imagine equal-sized ionized bubbles each with an interior neutral fraction ${x}_{{\rm{H}}{\rm{I}},\mathrm{IN}}\ll 1$, filling a fraction ${Q}_{{\rm{H}}{\rm{II}}}$ of the volume of the IGM, which is otherwise completely neutral. For simplicity, we neglect helium and consider an IGM with a uniform density and temperature. In the second model, representing uniform ionization, the neutral fraction is identical at each location within the IGM with $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle =1-{Q}_{{\rm{H}}{\rm{II}}}$. In each case, photoionization equilibrium tells us that $\langle {{\rm{\Gamma }}}_{{\rm{H}}{\rm{I}}}\rangle \sim \langle \alpha {n}_{{\rm{H}}}\rangle \langle {(1-{x}_{{\rm{H}}{\rm{I}}})}^{2}/{x}_{{\rm{H}}{\rm{I}}}\rangle$, where α here is the recombination factor. Now, in a uniformly ionized IGM we get $\langle {{\rm{\Gamma }}}_{{\rm{H}}{\rm{I}}}{\rangle }_{\mathrm{uniform}}\sim \langle \alpha {n}_{{\rm{H}}}\rangle {Q}_{{\rm{H}}\,{\rm{II}}}^{2}/(1-{Q}_{{\rm{H}}{\rm{II}}})$. In the toy patchy model, we have ${(1-{x}_{{\rm{H}}{\rm{I}}})}^{2}/{x}_{{\rm{H}}{\rm{I}}}\sim 1/{x}_{{\rm{H}}{\rm{I}},\mathrm{IN}}$ inside each ionized bubble and ∼0 outside of ionized regions. Hence, $\langle {{\rm{\Gamma }}}_{{\rm{H}}{\rm{I}}}{\rangle }_{\mathrm{patchy}}\,\sim {Q}_{{\rm{H}}{\rm{II}}}\langle \alpha {n}_{{\rm{H}}}\rangle /{x}_{{\rm{H}}{\rm{I}},\mathrm{IN}}$. The ratio of the volume-averaged photoionization rates is just $\langle {{\rm{\Gamma }}}_{{\rm{H}}{\rm{I}}}{\rangle }_{\mathrm{patchy}}/\langle {{\rm{\Gamma }}}_{{\rm{H}}{\rm{I}}}{\rangle }_{\mathrm{uniform}}\sim (1-{Q}_{{\rm{H}}{\rm{II}}})\,/{Q}_{{\rm{H}}{\rm{II}}}{x}_{{\rm{H}}{\rm{I}},\mathrm{IN}}$. This will typically be a very large number: for example, if 50% of the volume is filled by ionized bubbles ${Q}_{{\rm{H}}{\rm{II}}}=0.5$, each with an interior neutral fraction of ${x}_{{\rm{H}}{\rm{I}},\mathrm{IN}}\,={10}^{-4}$, the volume-averaged photoionization rate is a factor of 10⁴ times larger in the patchy reionization model than in the uniform model.

This is also relevant to understanding how the ionization history given by the volume filling factor can be compared with the one given by the volume-averaged ionization fractions and the intrinsic differences between the two. The Q formalism only knows about sources and sinks of ionizing photons and does not tell us anything about the value of the neutral fraction in highly ionized regions. In what follows we focus on hydrogen reionization, but analogous considerations also apply to helium. The volume filling factor evolution in UVB synthesis models is computed as

$$\begin{eqnarray}&&\displaystyle \frac{{{dQ}}_{{\rm{H}}{\rm{II}}}}{{dt}}=\displaystyle \frac{{\dot{n}}_{\mathrm{photon},{\rm{H}}{\rm{II}}}}{\langle {n}_{{\rm{H}}}\rangle }-\displaystyle \frac{{Q}_{{\rm{H}}{\rm{II}}}}{\langle {t}_{\mathrm{rec},{\rm{H}}}\rangle },\end{eqnarray} \tag{ 3 }$$

where ${t}_{\mathrm{rec},{\rm{H}}}$ is the hydrogen recombination time and ${\dot{n}}_{\mathrm{photon},{\rm{H}}{\rm{I}}}$ is the mean number of ionizing photons emitted by all radiation sources available per second.⁶ In the context of these models, ${\dot{n}}_{\mathrm{photon},{\rm{i}}}={\int }_{{\nu }_{L,i}}^{\infty }\tfrac{{\epsilon }_{\nu }}{h\nu }d\nu$ is considered to be the number of ionizing photons emitted into the IGM by all radiation sources (HM12; So et al. 2014), where ${\epsilon }_{\nu }$ is the total emissivity as a function of frequency obtained by the assumption of the sources (i.e., galaxies' and quasars' luminosity functions).⁷ This model uses observational constraints of the distribution of absorbers along the line of sight, and it is able to reproduce available measurements of the mean free path at 1 Ryd and the Lyα effective opacity.

To clarify what is happening at these high redshifts, we can use the equation of cosmological radiative transfer in its "source function" approximation, which allows us to write the following relation between emissivity, radiation intensity, and mean free path, ${\lambda }_{\nu }$, by ignoring photon redshifting effects (minimal at high redshifts): $4\pi {J}_{\nu }={\lambda }_{\nu }{\epsilon }_{\nu }$ if only local radiation sources contribute to the ionizing background intensity (HM12). From this relation it is now much more easy to see what went wrong in the UVB model. The ${\rm{H}}\,{\rm{I}}$ photoionization rates, ${{\rm{\Gamma }}}_{{\rm{H}}{\rm{I}}}$ given by HM12, are calculated as

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{H}}{\rm{I}}}=\int 4\pi \displaystyle \frac{{J}_{\nu }}{h\nu }{\sigma }_{\nu ,{\rm{H}}{\rm{I}}}d\nu ,\end{eqnarray} \tag{ 4 }$$

where ${J}_{\nu }$ is the radiation intensity. The ionization history, Q, was computed using the emissivity alone, using an analytical approximation that does not know anything about the mean free path assumed in the model. On the other hand, the photoionization rates depend on both the emissivity and the mean free path, indicating that the mean free path extrapolation done at high redshift was wrong and yields values that are systematically too high. To illustrate this, we show in Figure 2 the mean free path assumed in the model (green line) and observational constraints on the mean free path at 1 Ryd (black symbols). Notice that at high redshifts $z\gt 5$, there are no available constraints and the model thus corresponds to a blind extrapolation.

Zoom In Zoom Out Reset image size

The ${\rm{H}}\,{\rm{I}}$ photoheating rate (energy per time per ion), ${\dot{q}}_{{\rm{H}}{\rm{I}}}$, given by these models is calculated as

$$\begin{eqnarray}&&{\dot{q}}_{{\rm{H}}{\rm{I}}}=\int 4\pi \displaystyle \frac{{J}_{\nu }}{\nu }(\nu -{\nu }_{{\rm{H}}{\rm{I}}}){\sigma }_{\nu ,{\rm{H}}{\rm{I}}}d\nu ,\end{eqnarray} \tag{ 5 }$$

so the photoheating rates are overestimated for the same reason that the photoionization rates are. As noted above, the total heating rate does not depend on the amplitude of ${J}_{\nu }$, but only on its shape (Theuns et al. 2002b; McQuinn et al. 2009). However, as we argued above, in the standard UVB semianalytical synthesis modeling approach, Q and ${J}_{\nu }$ are not required to be internally consistent, so in practice the actual heating rate per volume approaches a constant at early times ($z\gt {z}_{\mathrm{reion}}$), whereas it should go to zero. Therefore, the total heat produced during ${\rm{H}}\,{\rm{I}}$ and $\mathrm{He}\,{\rm{II}}$ reionization in the simulation was also applied earlier than when it should be.

We have also tested other widely used UVB models in the literature (Haardt & Madau 1996, 2001; Faucher-Giguère et al. 2009, hereafter HM96, HM01, FG09, respectively). We present the full results for these other models in Appendix B, showing that they all share a similar problem, i.e., for ${\rm{H}}\,{\rm{I}}$ or $\mathrm{He}\,{\rm{II}}$ (or both) reionization the gas heating associated with the reionization event occurs much earlier than one will naively expect from these models. Although it was known that these UVB models do not properly model the reionization process, it is important that this problem has not been directly confronted in the literature (see, however, Puchwein et al. 2015). This oversight is likely due to the fact that, for the ionization and thermal histories of the IGM, most comparisons between simulations and observations have been performed in the context of the Lyα forest and focused on the redshift range for which such observations are available, i.e., between $2\lt z\lt 6$, or even lower redshifts. The thermal parameters studied so far in the literature, T₀ and γ, depend on the instantaneous values of these rates once we are far enough in time from reionization. So they quickly forget about early heating. Therefore, what happens at higher redshifts is not that relevant for these values. However, it is expected that this problem will have more important consequences for correctly modeling the pressure smoothing scale, as it depends on the full thermal history. This is because at IGM densities, the dynamical time that it takes the gas to respond to temperature changes at the Jeans scale (i.e., the sound-crossing time) is the Hubble time (Gnedin & Hui 1998).

To summarize, current UVB models result in different ionization and thermal histories than those quoted by the authors. This is because the low-z mean free path values have been blindly extrapolated to high z (see Figure 2), where they result in a photoionization rate far too high. Motivated by these results, we have developed a new method to build self-consistent effective ionization and thermal histories during reionization epochs in optically thin simulations, which we describe in the next section.

As stated in the previous section, current photoionization rates widely used in optically thin codes do not properly track the desired ionization histories and lead to incorrect thermal histories for the gas in these simulations, where the gas is heated much earlier than it should be. We present here a new way of creating self-consistent UVB models during reionization events (${\rm{H}}\,{\rm{I}}$, $\mathrm{He}\,{\rm{I}}$, and $\mathrm{He}\,{\rm{II}}$) to be used in optically thin hydrodynamical simulations. Different groups in the field have modified these tables, both ionization and photoheating rates (especially photoheating, and more in the context of $\mathrm{He}\,{\rm{II}}$ reionization where more observations are available), with different justifications: accounting for different physical effects as nonequilibrium ionization or radiative transfer effects, matching specific observables, or just exploring the parameter space (e.g., Haehnelt & Steinmetz 1998; Theuns et al. 2002a; Bolton et al. 2005; Jena et al. 2005; Pawlik et al. 2009; Wiersma et al. 2009; Lukić et al. 2015; Puchwein et al. 2015). However, this has generally been done by just applying a multiplying factor to the standard models assumed or by applying different simple cutoffs.⁸

Our approach will be based on building effective values of the photoionization and photoheating rates that can be substituted into the standard optically thin equations to yield the desired results (see FG09 for an early motivation of this approach). The main goal is to make sure that the heating due to reionization in the simulation is consistent with the reionization model itself. To enforce this, we have calculated the volume-averaged values of both the photoionization and photoheating rates that give us the desired ionization and total heat injection for an input reionization model. We give here a global overview of the method and the different assumptions made in our models. We explain in full detail how we derive the photoionization rates in Appendix A.1 and the photoheating rates in Appendix A.2. We will discuss the different caveats and limitations of our model in Section 8.

Each of our reionization models is defined by one free parameter, the total heat input ${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}$ during reionization, and one free function, the reionization history, which in this context we define as the volume-averaged ionization fraction evolution, $\langle {x}_{{\rm{H}}{\rm{II}}}\rangle (z)$. The reionization of $\mathrm{He}\,{\rm{II}}$ is analogously treated. Using these parameters, we will derive effective photoionization and photoheating rates that can be used in hydrodynamical simulations using the following assumptions:

• 1.  
  That all species are in ionization equilibrium at all times. This is done to be fully consistent with (most of) the codes that will be using these UVB models, but it could be changed in the future.
• 2.  
  That the gas composition can be approximated as primordial. Therefore, the evolution of the number density of electrons, ${n}_{{\rm{e}}}$, is given as ${n}_{{\rm{e}}}={n}_{{\rm{H}}{\rm{II}}}+{n}_{\mathrm{He}{\rm{II}}}+2{n}_{\mathrm{He}{\rm{III}}}$.
• 3.  
  That $\mathrm{He}\,{\rm{I}}$ reionization is perfectly coupled with ${\rm{H}}\,{\rm{I}}$ reionization.
• 4.  
  That $\mathrm{He}\,{\rm{II}}$ reionization is not relevant during ${\rm{H}}\,{\rm{I}}$ reionization and vice versa.
• 5.  
  That the heating due to reionization is perfectly coupled to the reionization process. Therefore, the heating can be written as a function of the total heat injection and the ionization history, i.e.,

  $$\begin{eqnarray}&&\displaystyle \frac{{{dT}}_{{\rm{H}}{\rm{I}}}}{{dt}}\propto {\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}\displaystyle \frac{d\langle {x}_{{\rm{H}}{\rm{II}}}\rangle }{{dt}}.\end{eqnarray} \tag{ 6 }$$

Finally, the new effective rates are only used during reionization, i.e., while $\langle {x}_{{\rm{H}}{\rm{II}}}\rangle \lt 1.0$. Once the reionization redshift, defined as when the input ionization history $\langle {x}_{{\rm{H}}{\rm{II}}}\rangle =1$, is reached, we can simply use the photoionization and photoheating rates of common UVB models (in our case HM12).

Let us first focus on how we obtain the new photoionization rates to be applied during reionization. We obtain the new effective photoionization rates by volume-averaging the ionization equilibrium equations and using the assumptions enumerated above. In particular, for the ${\rm{H}}\,{\rm{I}}$ photoionization rates we get

$$\begin{eqnarray}&&\langle {{\rm{\Gamma }}}_{\gamma ,{\rm{H}}{\rm{I}}}\rangle (z)={C}_{{\rm{H}}{\rm{II}}}\langle {n}_{{\rm{H}}}\rangle (z){\alpha }_{{\rm{r}},{\rm{H}}{\rm{II}}}(\langle T\rangle )(1+\chi )\displaystyle \frac{\langle {x}_{{\rm{H}}{\rm{II}}}{\rangle }^{2}(z)}{\langle {x}_{{\rm{H}}{\rm{I}}}\rangle (z)},\end{eqnarray} \tag{ 7 }$$

where ${C}_{{\rm{H}}{\rm{II}}}$ is a volume-averaged correction factor linked with the well-known clumping factor and ${\alpha }_{{\rm{r}},{\rm{H}}{\rm{II}}}(\langle T\rangle )$ is the recombination coefficient for which a specific volume-averaged temperature of the IGM, $\langle T\rangle$, has to be assumed. We refer the reader to Appendix A.1 for all the details.⁹

In addition to affecting abundances, photoionization injects energy into the gas when a high-energy ($h\nu \gt h{\nu }_{{\rm{T}}}$) photon transfers more energy to an electron than what is necessary to unbind it from the atom. In order to include this heat transfer to the IGM by the UVB, models also provide photoheating rates that are included in simulations as an effective energy source. The temperature of the IGM is much more affected by this photoheating and different atomic cooling processes than by the gravitational collapse of cosmic structure. These processes are accounted for as a global heating and cooling term that is added to the equation of energy for the gas:

$$\begin{eqnarray}\displaystyle \frac{\partial E}{\partial t} & = & -\displaystyle \frac{1}{a}{\boldsymbol{v}}\cdot {\rm{\nabla }}E-\displaystyle \frac{\dot{a}}{a}\left(3\displaystyle \frac{p}{\rho }+{{\boldsymbol{v}}}^{2}\right)\\ & & -\displaystyle \frac{1}{\rho a}{\rm{\nabla }}\cdot (p{\boldsymbol{v}})+\displaystyle \frac{1}{a}{\boldsymbol{v}}\cdot {\rm{\nabla }}{\rm{\Phi }}+\displaystyle \frac{{{\rm{\Lambda }}}_{\mathrm{HC}}}{a\rho },\end{eqnarray} \tag{ 8 }$$

where ${{\rm{\Lambda }}}_{\mathrm{HC}}$ represents the combined heating and cooling terms,¹⁰ which can be expanded as

$$\begin{eqnarray}{{\rm{\Lambda }}}_{\mathrm{HC}} & = & {a}^{4}{n}_{{\rm{H}}{\rm{I}}}{\dot{q}}_{{\rm{H}}{\rm{I}}}+{a}^{4}{n}_{\mathrm{He}{\rm{I}}}{\dot{q}}_{\mathrm{He}{\rm{I}}}+{a}^{4}{n}_{\mathrm{He}{\rm{II}}}{\dot{q}}_{\mathrm{He}{\rm{II}}}\\ & & -{a}^{4}{n}_{{\rm{e}}}{f}_{{\rm{c}}}({n}_{{\rm{H}}{\rm{I}}},{n}_{{\rm{H}}{\rm{II}}},{n}_{\mathrm{He}{\rm{I}}},\ldots )\\ & & -{a}^{4}{n}_{{\rm{e}}}{f}_{\mathrm{Compton}}(T,{T}_{\mathrm{CMB}})+...,\end{eqnarray} \tag{ 9 }$$

where the first three terms represent the photoheating rates of ${\rm{H}}\,{\rm{I}}$, $\mathrm{He}\,{\rm{I}}$, and $\mathrm{He}\,{\rm{II}}$, respectively.¹¹ We have also included an atomic cooling term (${f}_{{\rm{c}}}$) and an inverse-Compton scattering term (${f}_{\mathrm{Compton}}$).

We can obtain effective photoheating rates for the reionization heating by volume-averaging the relation between the heat per unit of time produced by a certain reionization model, $d{\rm{\Delta }}T/{dt}$, and the one produced by a photoheating rate. This will allow the use of our new models in standard hydrodynamical codes. For the ${\rm{H}}\,{\rm{I}}$ photoheating rate, and using the set of assumptions enumerated above, we obtain

$$\begin{eqnarray}&&{\dot{q}}_{{\rm{H}}{\rm{I}}}={C}_{\dot{q},{\rm{H}}{\rm{I}}}\displaystyle \frac{3{k}_{{\rm{B}}}}{2\langle \mu \rangle {X}_{{\rm{p}}}\langle {x}_{{\rm{H}}{\rm{I}}}\rangle (z)}\displaystyle \frac{d{\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}}{{dt}},\end{eqnarray} \tag{ 10 }$$

where ${X}_{{\rm{p}}}$ is the hydrogen mass abundance, $\langle \mu \rangle$ is the volume-averaged molecular weight, $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle (z)$ is the assumed reionization history of the model, ${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}$ is its total heat input, and ${C}_{\dot{q},{\rm{H}}{\rm{I}}}$ is a volume-averaged correction factor that we set to 1 at all redshifts. We refer the reader to Appendix A.2 for more details on how this equation was obtained.

Figure 3 shows an example of how a UVB model is built. In the left panel we show a model that should match a specific late ${\rm{H}}\,{\rm{I}}$ reionization history that assumes a reionization redshift of ${z}_{\mathrm{rei}}=6.55$ (late reionization model; see more details in Section 4) with a total heat input due to ${\rm{H}}\,{\rm{I}}$ reionization of ${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}=2\times {10}^{4}$ K. The upper left panel shows the ${\rm{H}}\,{\rm{I}}$ ionization history assumed to build the model. The middle panel shows the assumed temperature evolution of the ${\rm{H}}\,{\rm{I}}$ reionization event. As explained above, we have assumed that the temperature evolution follows the ionization history. The lower left panel shows the actual instantaneous heat input, $d{\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}/{dz}$, derived using Equation (33) for this model. The integral of this line gives the total heat input, ${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}$.

Zoom In Zoom Out Reset image size

The upper right panel gives the ${\rm{H}}\,{\rm{I}}$ photoionization rate (blue line) obtained using Equation (7) and the method described above. This rate is compared with the equivalent of the HM12 model (green line). We also plot observational constraints on the photoionization rates from different studies (Bolton & Haehnelt 2007; Faucher-Giguère et al. 2008; Calverley et al. 2011; Wyithe & Bolton 2011; Becker & Bolton 2013). Notice that all these constraints are below z = 6. The middle right panel shows the photoheating rates derived using Equation (10) for the ${\rm{H}}\,{\rm{I}}$. Notice that for the new model, below its reionization redshift (${z}_{\mathrm{reion},{\rm{H}}{\rm{I}}}=6.55$) we apply the same photoionization and photoheating rates used in HM12, as we think that they are reasonable after reionization. The rates of our model raise abruptly at around the reionization redshift because in the input ${\rm{H}}\,{\rm{I}}$ reionization model the transition from $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle \sim 0.1$ to $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle \sim {10}^{-4}$ is very fast. The lower right panel shows the reionization function input for the model (solid blue line; defined to be 1.0 at $z={z}_{\mathrm{reion},{\rm{H}}{\rm{I}}}$). It also shows the volume-averaged ${\rm{H}}\,{\rm{I}}$ ionization fraction obtained running a cosmological hydrodynamical simulation that uses this new photoionization rate (dashed blue line). The solid black line shows the reionization evolution of the HM12 model as calculated by their authors, while the dashed green line stands for the outcome of a hydrodynamical simulation using the HM12 photoionization rate. Current best observational constraints on the ${\rm{H}}\,{\rm{I}}$ fraction at different redshifts are also plotted using different symbols and error bars (Bolton et al. 2005; Fan et al. 2006; Bolton & Haehnelt 2007; McGreer et al. 2015). Notice again that all available constraints are below z = 6 so that both models are in agreement with these observations.

We can derive the evolution of the mean free path at 1 Ryd for the new model, making some simple assumptions, and compare it with the ones obtained by HM12. Our modeling has no shape information on the intensity, ${J}_{\nu }$; however, one can assume the HM12 shape and that all intensities differ by the same constant: ${J}_{\nu ,\mathrm{new}}=C\times {J}_{\nu ,\mathrm{HM}12}$. With this we can approximate the new radiation intensity at 1 Ryd as the ratio between the two photoionization rates ${J}_{912,\mathrm{new}}=({{\rm{\Gamma }}}_{{\rm{H}}{\rm{I}},\mathrm{new}}/{{\rm{\Gamma }}}_{{\rm{H}}{\rm{I}},\mathrm{HM}12})\times {J}_{912,\mathrm{HM}12}$. Then we obtain the mean free path using the source function approximation that has the same emissivity assumed by HM12: ${\lambda }_{\mathrm{mfp},912,\mathrm{new}}=4\pi {J}_{912,\mathrm{new}}/{\epsilon }_{912,\mathrm{HM}12}$. Figure 2 shows the evolution of the mean free path derived for this new model (blue line). The green solid line stands for the mean free path used in the HM12 model. This plot illustrates how the volume-averaged mean free path drops significantly above the reionization redshift.

We will now discuss the different reionization models that we want to simulate by using our new approach to compute photoionization and photoheating during reionization. In order to define a reionization model, we need to set the ${\rm{H}}\,{\rm{I}}$ and $\mathrm{He}\,{\rm{II}}$ reionization history via the volume-averaged ionization fractions.¹²

First, we need to define the shape of our reionization histories. We use the lower incomplete gamma function, g:

$$\begin{eqnarray}\langle {x}_{{\rm{H}}{\rm{II}}}\rangle =\left\{\begin{array}{ll}0.5+0.5\times g(1/{n}_{1},| z-{z}_{0}{| }^{{n}_{1}}), & z\lt =\,{z}_{0}\\ 0.5-0.5\times g(1/{n}_{2},| z-{z}_{0}{| }^{{n}_{2}}), & z\gt {z}_{0}\end{array}\right.,\end{eqnarray} \tag{ 11 }$$

where ${n}_{1}=50$, ${n}_{2}=1$, and z₀ is a free parameter that sets the redshift where ${x}_{{\rm{H}}{\rm{II}}}({z}_{0})=0.5$. This function defines a slower start for the reionization function but a fast finish. This specific shape was motivated by radiative transfer simulation results (Ahn et al. 2012; Park et al. 2013; Pawlik et al. 2015), which seem to favor rapid and sudden ${\rm{H}}\,{\rm{I}}$ reionization histories. In particular, n₁ and n₂ were fitted to mimic Pawlik et al. (2015) results. In this work we use only this shape and experiment with the redshift of reionization, but our function could be easily modified in the future to explore a wider range of models. We explore the relevant range of reionization parameters taking into account the CMB constraints on the integrated electron scattering optical depth, ${\tau }_{{\rm{e}}}=0.078\pm 0.019$ (Planck Collaboration et al. 2016a).¹³ We consider three models: an early, middle, and late ${\rm{H}}\,{\rm{I}}$ reionization history, which have reionization redshifts¹⁴ of ${z}_{\mathrm{reion},{\rm{H}}{\rm{I}}}=9.70$, 8.30, and 6.55, respectively, all of which give ${\tau }_{{\rm{e}}}$ within $1\sigma$ of the CMB measurements. These ${\rm{H}}\,{\rm{I}}$ reionization models are plotted as solid lines in the upper panel of Figure 4, and their associated ${\tau }_{{\rm{e}}}$ values are shown in the lower panel.

Zoom In Zoom Out Reset image size

For $\mathrm{He}\,{\rm{II}}$ reionization we consider a history based on FG09 (last column of their Table 2) results, which sets the reionization redshift to finish at ${z}_{\mathrm{reion},\mathrm{He}{\rm{II}}}=3.0$ (hereafter He_A). The specific evolution of the full ionized fraction can be well described by $\langle {x}_{\mathrm{He}{\rm{III}}}\rangle =1.0-\arctan (z-{z}_{\mathrm{reion}})$. This is in agreement with standard models of $\mathrm{He}\,{\rm{II}}$ reionization, which, as a result of the high-energy requirement to double ionize helium (${E}_{\nu }\gt 54.4$ eV), assume that quasars must be the main drivers of this process (e.g., Miralda-Escudé et al. 2000; Compostella et al. 2014; Worseck et al. 2016). This $\mathrm{He}\,{\rm{II}}$ reionization model is plotted as a dot-dashed line in the upper panel of Figure 4.

4.1. Total Heat Input due to H i and He ii Reionization

The simulations used in this work were performed with the Nyx code (Almgren et al. 2013). Nyx follows the evolution of dark matter simulated as self-gravitating Lagrangian particles, and baryons modeled as an ideal gas on a uniform Cartesian grid. The Eulerian gas dynamics equations are solved using a second-order-accurate piecewise parabolic method to accurately capture shock waves. We do not make use of adaptive mesh refinement (AMR) capabilities of Nyx in the current work, as the Lyα forest signal spans nearly the entire simulation domain rather than isolated concentrations of matter, where AMR is more effective. For more details of these numerical methods and scaling behavior tests, see Almgren et al. (2013).

Besides solving for gravity and the Euler equations, we also include the main physical processes fundamental to modeling the Lyα forest. First, we consider the chemistry of the gas as having a primordial composition with hydrogen and helium mass abundances of ${X}_{{\rm{p}}}$ and ${Y}_{{\rm{p}}}$, respectively. In addition, we include inverse-Compton cooling off the microwave background and keep track of the net loss of thermal energy resulting from atomic collisional processes. We used the updated recombination, collision ionization, dielectric recombination rates, and cooling rates given in Lukić et al. (2015). All cells are assumed to be optically thin, and radiative feedback is accounted for via a spatially uniform but time-varying UVB radiation field given to the code as a list of photoionization and photoheating rates that vary with redshift following the method described in Section 3.

In order to generate the initial conditions, we have used the music code (Hahn & Abel 2011) and a camb (Lewis et al. 2000; Howlett et al. 2012) transfer function. All simulations started at ${z}_{\mathrm{ini}}=159$ to be sure that nonlinear evolution is not compromised (see, e.g., Oñorbe et al. 2014, for a detailed discussion on this issue). Unless otherwise stated, all the simulations discussed in this paper assumed a ΛCDM cosmology with the following fundamental parameters: ${{\rm{\Omega }}}_{{\rm{m}}}=0.3192$, ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.6808$, ${{\rm{\Omega }}}_{{\rm{b}}}=0.04964$, h = 0.67038, ${\sigma }_{8}=0.826$, and ${n}_{{\rm{s}}}=0.9655$. These values are within 1σ agreement with last cosmological constraints from the CMB (Planck Collaboration et al. 2016a). The choice of hydrogen and helium mass abundances (${X}_{{\rm{p}}}=0.76$ and ${Y}_{{\rm{p}}}=0.24$, therefore $\chi =0.0789$) is in agreement with the recent CMB observations and big bang nucleosynthesis (Coc et al. 2013). Simulations were run down to z = 0.2, saving 32 snapshots¹⁵ from z = 20. Unless otherwise stated, all simulations presented here have a box size of length ${L}_{\mathrm{box}}=20$ $\mathrm{Mpc}/h$ and 1024³ resolution elements. This dynamical range guarantees that all the different physical parameters analyzed in this paper are converged with enough accuracy ($\lt 5 \%$ errors). We will discuss this issue in more detail in Section 8.1.

We first run one simulation using photoionization and photoheating values from the most widely used models (HM96, HM01, FG09, and HM12). We have already presented the reionization and thermal histories of the HM12 model in Section 2, but below we will further explore other properties of this simulation.¹⁶ We also run the three ${\rm{H}}\,{\rm{I}}$ reionization histories presented above, an early, middle, and late reionization model (EarlyR, MiddleR, and LateR; see above and Figure 4). All of them share the same heat input during ${\rm{H}}\,{\rm{I}}$ reionization, ${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}=2\times {10}^{4}$ K, and $\mathrm{He}\,{\rm{III}}$ reionization model: $\langle {x}_{\mathrm{He}{\rm{III}}}\rangle (z)$ and ${\rm{\Delta }}{T}_{\mathrm{He}{\rm{III}}}=1.5\times {10}^{4}$ K.

In order to study the effect of different total heat input during ${\rm{H}}\,{\rm{I}}$ reionization, ${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}$, we run three more simulations that share all parameters with the ${\rm{H}}\,{\rm{I}}$ middle reionization run, MiddleR, but varying this parameter: MiddleR-Hcold (${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}=1.5\times {10}^{4}\,{\rm{K}}$), MiddleR-Hwarm (${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}=3\times {10}^{4}\,{\rm{K}}$), and MiddleR-Hhot (${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}=4\times {10}^{4}\,{\rm{K}}$). We also explored the effects of different global net heating during $\mathrm{He}\,{\rm{II}}$ reionization by running a set of four more ${\rm{H}}\,{\rm{I}}$ middle reionization simulations (MiddleR, ${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}=2\times {10}^{4}$), in which we just changed the heat input during $\mathrm{He}\,{\rm{II}}$ reionization, ${\rm{\Delta }}{T}_{\mathrm{He}{\rm{II}}}$: MiddleR-noHe (no $\mathrm{He}\,{\rm{II}}$ reionization), MiddleR-Hecold (${\rm{\Delta }}{T}_{\mathrm{He}{\rm{II}}}=1\times {10}^{4}\,{\rm{K}}$), MiddleR-Hewarm (${\rm{\Delta }}{T}_{\mathrm{He}{\rm{II}}}=2\times {10}^{4}\,{\rm{K}}$), and MiddleR-Hehot (${\rm{\Delta }}{T}_{\mathrm{He}{\rm{II}}}=3\times {10}^{4}\,{\rm{K}}$). 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.

Table 1.  Summary of Simulations

Sim	$\langle {x}_{{\rm{H}}{\rm{II}}}\rangle (z)$	${z}_{\mathrm{reion},{\rm{H}}{\rm{I}}}$	${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}$	$\langle {x}_{\mathrm{He}{\rm{III}}}\rangle (z)$	${\rm{\Delta }}{T}_{\mathrm{He}{\rm{II}}}$
			(K)		(K)
HM12	HM12 a	⋯	⋯	⋯	⋯
FG09	FG09 b	⋯	⋯	⋯	⋯
HM01	HM01 c	⋯	⋯	⋯	⋯
HM96	HM96 d	⋯	⋯	⋯	⋯
LateR	Late reionization	6.55	2 × 104	He_A	1.5 × 104
MiddleR	Middle reionization	8.30	2 × 104	He_A	1.5 × 104
EarlyR	Early reionization	9.70	2 × 104	He_A	1.5 × 104
MiddleR-Hcold	Middle reionization	8.30	1.5 × 104	He_A	1.5 × 104
MiddleR-Hwarm	Middle reionization	8.30	3 × 104	He_A	1.5 × 104
MiddleR-Hhot	Middle reionization	8.30	4 × 104	He_A	1.5 × 104
MiddleR-noHe	Middle reionization	8.30	2 × 104	None	⋯
MiddleR-Hecold	Middle reionization	8.30	2 × 104	He_A	1 × 104
MiddleR-Hewarm	Middle reionization	8.30	2 × 104	He_A	2 × 104
MiddleR-Hehot	Middle reionization	8.30	2 × 104	He_A	3 × 104

Notes. Column (1): simulation code. Column (2): ${\rm{H}}\,{\rm{I}}$ ionization history assumed for the model. See Section 4 for details. Column (3): ${\rm{H}}\,{\rm{I}}$ reionization redshift. Column (4): total heating assumed for ${\rm{H}}\,{\rm{I}}$ reionization. Column (5): $\mathrm{He}\,{\rm{III}}$ ionization history assumed for the model. All models assume a ${z}_{\mathrm{reion},\mathrm{He}{\rm{II}}}=3.0$. See Section 4 for details. Column (6): total heating assumed for $\mathrm{He}\,{\rm{II}}$ reionization. Unless otherwise stated, all simulations have a box size of length ${L}_{\mathrm{box}}=20$ $\mathrm{Mpc}/h$ and 1024³ resolution elements.

^aUsing Haardt & Madau (2012) tabulated values. ^bUsing Faucher-Giguère et al. (2009) tabulated values (2011 December update). ^cUsing Haardt & Madau (2001) tabulated values. ^dUsing Haardt & Madau (1996) tabulated values.

Download table as:  ASCIITypeset image

Table 2.  Cosmological Parameters Used in Simulations

Model	${{\rm{\Omega }}}_{{\rm{m}}}$	${{\rm{\Omega }}}_{{\rm{b}}}$	${{\rm{\Omega }}}_{{\rm{\Lambda }}}$	h	${\sigma }_{8}$	${n}_{{\rm{s}}}$
A	0.320	0.0496	0.681	0.670	0.826	0.966
B	0.298	0.0477	0.702	0.686	0.873	0.974
C	0.333	0.0517	0.667	0.658	0.757	0.971
D	0.275	0.0460	0.725	0.702	0.816	0.960

Note. Column (1): model name. Column (2): total matter content. Column (3): total baryonic matter content. Column (4): cosmological constant. Column (5): normalized Hubble constant (0–1). Column (6): normalization of the power spectrum. Column (7): spectral index.

Download table as:  ASCIITypeset image

5.1. Analysis of the Simulations

Whenever Lyα forest spectra are created from the simulation, we compute the ${\rm{H}}\,{\rm{I}}$ optical depth at a fixed redshift, which can then be easily converted into a transmitted flux fraction, $F={e}^{-\tau }$. That is, we do not account for the speed of light when we cast rays in the simulation; we use the gas state at a single cosmic time. The simulated spectra are not meant to look like full Lyα forest spectra, but just recover the statistics of the flux in a small redshift window. Our calculation of the spectra accounts for Doppler shifts due to bulk flows of the gas, as well as for thermal broadening of the Lyα line. We refer to Lukić et al. (2015) for specific details of these calculations. This procedure results in the Lyα flux as a function of wavelength or equivalently time or distance. Following the standard approach, we then rescale the UV background intensity so that the mean flux of all the extracted spectra from the simulation matches the observed mean flux at the respective redshift (see Section 7 for more details on the specific value that we have chosen). We therefore have neglected noise and metal contamination in our skewers so far, but this will not be relevant in this paper.

Using these skewers, we have also calculated the curvature flux statistics, $\langle | \kappa | \rangle$, where

$$\begin{eqnarray}&&\kappa =\displaystyle \frac{F^{\prime\prime} }{{[1+{(F^{\prime} )}^{2}]}^{3/2}},\end{eqnarray} \tag{ 12 }$$

$F^{\prime}$ is the first derivative of the flux with respect to the velocity separation between pixels, and $F^{\prime\prime}$ is the second derivative. We have done this for each simulation following the method described in Becker et al. (2011).¹⁷

We measured the thermal parameters of the simulation at each snapshot by fitting the ${\rho }_{{\rm{b}}}\mbox{--}T$ relation with linear least squares in ${\mathrm{log}}_{10}\,{\rm{\Delta }}$ and ${\mathrm{log}}_{10}\,T$, fitting the range $-0.7\lt {\mathrm{log}}_{10}\,{\rm{\Delta }}\lt 0.0$ and ${\mathrm{log}}_{10}\,T/K\lt 4.5$.¹⁸

To characterize the gas pressure support in all our simulations, we have followed the recent work by Kulkarni et al. (2015) and use the real-space Lyα flux, ${F}_{{\rm{H}}{\rm{I}},\mathrm{real}}$. This quantity is defined as ${F}_{{\rm{H}}{\rm{I}},\mathrm{real}}=\exp (-{\tau }_{{\rm{H}}{\rm{I}},\mathrm{real}})$, where ${\tau }_{{\rm{H}}{\rm{I}},\mathrm{real}}$ is the real-space Lyα optical depth, which is identical to the observed Lyα optical depth except that the convolution integral that accounts for the redshift-space effects of the peculiar velocity field and thermal line broadening has not been included. This field naturally suppresses dense gas and is thus robust against the poorly understood physics of galaxy formation, revealing pressure smoothing in the diffuse IGM. The ${F}_{{\rm{H}}{\rm{I}},\mathrm{real}}$ 3D power spectrum is accurately described by a simple fitting function with a Gaussian cutoff at ${\lambda }_{{\rm{P}}}$, which is then defined as the pressure smoothing scale. This statistic has the added advantage that it directly relates to observations of correlated Lyα forest absorption in close quasar pairs, proposed as a method to measure this scale, and enables one to quantify it in simulations (Rorai et al. 2013; A. Rorai et al. 2015, in preparation).

We now present thermal properties of the IGM in LateR, MiddleR, and EarlyR simulations, which only differ in their redshift of ${\rm{H}}\,{\rm{I}}$ reionization. We first focus on the evolution with redshift of the temperature at mean density, T₀, and the slope of the temperature–density relation, γ. The left panel of Figure 5 shows the evolution of these parameters for the HM12 run (green line). It also shows the thermal history of the LateR (blue), MiddleR (magenta), and EarlyR (orange) simulations. In the upper left panel we plot the evolution of γ, which exhibits the expected convergence to a value close to ∼1.6 after all reionization events for all models, resulting from the balance of photoheating with adiabatic cooling (Hui & Gnedin 1997; McQuinn & Upton Sanderbeck 2016). The larger decrease of γ during $\mathrm{He}\,{\rm{II}}$ reionization in the EarlyR, MiddleR, and LateR runs seems to indicate a temperature increase more independent of density than in the HM12 run. Puchwein et al. (2015) have also shown that, for a fixed UVB model, using a nonequilibrium approach will tend to create a more pronounced feature (we will further discuss this in Section 8). In our modeling we use equilibrium photoionization; hence, the flattening occurs for different reasons. In fact, this is just because our ionization model injects more heat into the IGM than the model in HM12. One expects this type of effect during reionization when applying a uniform UVB model through the whole volume, as in that case we are applying a constant temperature increase at each resolution element. At higher initial temperature this corresponds to a lower increase in the logarithm of the temperature, so that the temperature–density relation flattens in log–log space. This effect gets magnified as we increase the amount of heat applied to the whole box. Therefore, our results show that the exact evolution of γ in simulations depends on the assumed shape for the ${\rm{H}}\,{\rm{I}}$ and $\mathrm{He}\,{\rm{II}}$ reionization histories and their total heat input. In this panel we also plot the value of γ of Bolton et al. (2014) at z = 2.4, derived from absorption-line profiles in the Lyα forest.¹⁹

Zoom In Zoom Out Reset image size

The lower left panel of Figure 5 shows the evolution of T₀ for the same set of simulations. As expected, in the new models ${\rm{H}}\,{\rm{I}}$ reionization produces a much later heating than in the HM12 run. In fact, it can be clearly seen that the temperature at mean density of the new runs rises following the ${\rm{H}}\,{\rm{I}}$ reionization of each model. The heating during $\mathrm{He}\,{\rm{II}}$ reionization also shows significant differences with the HM12 run. Although the heating happens at basically the same time as in the HM12 run, the new models also exhibit a larger and sharper temperature increase at lower redshifts ($3\lesssim z\lesssim 4$), along with the expected decrease of the γ value. This is due to the different ionization history assumed for $\mathrm{He}\,{\rm{II}}$, which rises steeply at these redshifts. We compare these models to the measurements of the temperature at mean density, T₀, from Lidz et al. (2010), based on the wavelet technique (Meiksin 2000; Theuns & Zaroubi 2000). We also plot constraints obtained by combining the γ measurement by Bolton et al. (2014) plotted in the upper panel with measurements of the temperature at the optimal density, $T({{\rm{\Delta }}}_{{\bigstar}})$, by Becker et al. (2011) and Boera et al. (2014) derived from the observed curvature of the Lyα forest transmitted flux. In this case we propagated errors from both measurements. Since observations of the Lyα forest are only available at $z\lesssim 6$, Figure 5 illustrates that it will be challenging to constrain ${\rm{H}}\,{\rm{I}}$ reionization from measurements of T₀ and γ at these redshifts because the IGM quickly loses thermal memory of reionization. Our new UVB models provide, by construction, a sharper temperature rise due to $\mathrm{He}\,{\rm{II}}$ reionization than the HM12 run. This again illustrates that the exact evolution of the IGM thermal parameters depends on the assumed shape for the ${\rm{H}}\,{\rm{I}}$ and $\mathrm{He}\,{\rm{II}}$ reionization histories, as well as their total heat input.

In the upper right panel of Figure 5 we plot another interesting property defining the thermal state of the IGM, which is the temperature at the optimal overdensity, $T({{\rm{\Delta }}}_{{\bigstar}})$. This "optimal" overdensity at each redshift is defined as the one for which curvature measurements of the Lyα forest are more sensitive (see Becker et al. 2011; Boera et al. 2014). The curvature measurements allow one to determine this parameter because they are not able to break the degeneracy between T₀ and γ. To calculate the $T({{\rm{\Delta }}}_{{\bigstar}})$ of our simulations, we have used the function fit to these optimal densities by Becker et al. (2011) from a suite of hydrodynamical simulations (${\mathrm{log}}_{10}{{\rm{\Delta }}}_{{\bigstar}}=A\times z+B$, where $A=-0.24596$ and B = 1.22218).²⁰ Inspecting Figure 5, we see that all simulations give the same temperature at the optimal density at $z\lt 6$, which is not surprising as we have already shown that at these redshifts all of them have roughly the same temperature–density relation (see the left panel of Figure 5). We also see the more pronounced rise in temperature due to $\mathrm{He}\,{\rm{II}}$ reionization in the new models compared with HM12 between $3\lesssim z\lesssim 4$.

The middle right panel of Figure 5 shows the curvature flux statistics, $\langle | \kappa | \rangle$, for all these simulations, and it is clear that they do not match at these redshifts. This is because, as explained above, the flux statistics depends not only on the temperature–density relation of the IGM (γ, T₀) but also on the pressure smoothing scale of the IGM, ${\lambda }_{{\rm{P}}}$. In fact, the lower right panel of Figure 5 shows the evolution of the pressure smoothing scale, ${\lambda }_{{\rm{P}}}$, with redshift for all the simulations. This panel clearly illustrates the dependence of the pressure smoothing scale on the full thermal history of the universe and not just on the instantaneous temperatures. That is, whereas the temperatures of all models agree at $z\lt 7$, differences in the pressure smoothing scale persist to much lower redshifts. This explains the differences between the curvature statistics for all the simulations and indicates that for fixed values of T₀ and γ the value of the optimal density, ${{\rm{\Delta }}}_{{\bigstar}}$, will be degenerate with the pressure smoothing scale, ${\lambda }_{{\rm{P}}}$.

We have confirmed this by recreating the same study done by Becker et al. (2011) to obtain the optimal densities using hydrodynamical simulations. With a similar set of simulations to that of these authors, we found almost identical results for the values of the optimal densities.²¹ However, including simulations that have identical values of T₀ and γ but different values of ${\lambda }_{{\rm{P}}}$ complicates the simple unique definition of ${{\rm{\Delta }}}_{{\bigstar}}$ by Becker et al. (2011) and instead adds scatter to the relationship between curvature and $T({{\rm{\Delta }}}_{{\bigstar}})$. The upper panel of Figure 5 also shows the determinations of the temperature at the optimal density using the curvature of the Lyα forest transmitted flux (Becker et al. 2011; Boera et al. 2014) compared with our simulations. Given the strong dependence of the curvature on the pressure smoothing scale ${\lambda }_{{\rm{P}}}$ resulting from the different reionization histories, it is clear that the error bars on $T({{\rm{\Delta }}}_{{\bigstar}})$ are likely underestimated. These issues pertaining to the thermal history are discussed in Becker et al. (2011; see also Puchwein et al. 2015) but were not included in the error budget. The difference in $\mathrm{log}\langle | \kappa | \rangle$ between our late reionization model (LateR) and early reionization model (EarlyR) at $z\sim 4$ is ∼0.05. From the results presented by Becker et al. (2011, Figures 1 and 10) this difference implies already 20%–25% error in temperature.

The aforementioned issues related to the pressure smoothing scale and thermal history can ease the $2\sigma$ level of disagreement between Lidz et al. (2010) measurements of T₀ and the Becker et al. (2011) measurements of $T({{\rm{\Delta }}}_{{\bigstar}})$ at $z\gt 4$ (Figure 5), although this does not seem to be enough to explain it fully. At lower redshifts a comparison of the two measurements is more challenging because it depends on what one assumes for the temperature–density relation slope, γ. Based on our results, it is clear that these conflicting measurements lead to different interpretations of the ${\rm{H}}\,{\rm{I}}$ and $\mathrm{He}\,{\rm{II}}$ reionization events. On the one hand, the Becker et al. (2011) results point toward a temperature of the IGM at mean density of ${T}_{0}\sim 1\times {10}^{4}$ K by $z\sim 4.7$ and a clear heating event later at $z\sim 3$, which they associated with $\mathrm{He}\,{\rm{II}}$ reionization. On the other hand, Lidz et al. (2010) higher measured temperatures require a higher-energy injection (${{\rm{\Delta }}}_{T}$) than we assumed in our models (${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}=2\times {10}^{4}$ K) and an earlier injection of heat that could be associated with $\mathrm{He}\,{\rm{II}}$ reionization at higher redshift than inferred by Becker et al. (2011) and Boera et al. (2014). Based on their measurements, Lidz et al. (2010) claim that the $\mathrm{He}\,{\rm{II}}$ reionization event should be completed by z = 3.4 and that the temperature at lower redshifts is consistent with the fall-off expected from adiabatic cooling. Although one can argue about the statistical significance of these discrepancies, especially given that the error bars are underestimated because neither study marginalized out the pressure smoothing scale ${\lambda }_{{\rm{P}}}$, we see no reason to prefer one set of measurements over the other. For this reason we will not attempt any further interpretation of these measurements with our numerical simulations, and we defer detailed data-to-model comparisons to future work.

We now want to discuss the evolution of the pressure smoothing scale in the different simulations, which, as can clearly be seen from Figure 5, retains memory of the reionization events. As we mentioned in the introduction, this is because, at IGM densities, the dynamical time that it takes the gas to respond to temperature changes at the Jeans scale (i.e., the sound-crossing time) is the Hubble time. The first thing to notice is that the HM12 model results in a much larger pressure smoothing scale than that of any our models, even the early reionization one. This is a direct result of the premature reionization of and the spurious associated heating (starting at $z\sim 15$) produced by this model. Our new models correct this issue, properly tying reionization heating to reionization history, resulting in later heating and a smaller overall pressure scale. In addition, significant differences are also found below z = 6 among the simulations with different ${\rm{H}}\,{\rm{I}}$ reionization histories (LateR, MiddleR, and EarlyR), even though these simulations share exactly the same photoionization and photoheating values at these redshifts and therefore have very similar thermal parameters, i.e., T₀, γ, and $T({\rm{\Delta }})$. These differences in ${\lambda }_{{\rm{P}}}$ arise because the IGM has had more time to respond to its hotter temperature when reionization occurs earlier, resulting in a larger pressure scale. The sensitivity of the pressure smoothing scale to the reionization history highlights the importance of constructing self-consistent models of reionization and applying them to optically thin simulations to better understand the thermal evolution of the IGM.

It is interesting to discuss the redshift evolution that we find for the pressure smoothing scale, using ${\lambda }_{{\rm{P}}}$. This parameter can be seen as the gas pressure scale at the density most sensitive for Lyα observations. Therefore, there are two physical effects that contribute to the value of this parameter. First, the IGM is heated as the universe evolves, so we expect the pressure smoothing scale to increase with time. On the other hand, there is not just one pressure smoothing scale in the IGM, but one for each density, and, just from linear theory, we expect it to be higher at lower densities, ${\lambda }_{{\rm{P}}}\propto {n}_{{\rm{H}}}^{-1/2}$ (Schaye 2001). As we go to lower redshifts, the neutral hydrogen density is being further diluted by the expansion of the universe, and observations start to be more sensitive to higher densities that have a smaller pressure smoothing scale. The combination of both processes can produce the somewhat surprising behavior of flattening of ${\lambda }_{{\rm{P}}}$ at $z\lesssim 4$ (decrease in the case of the HM12 model due to an $\mathrm{He}\,{\rm{II}}$ reionization with a lower heat injection). We also found that since all these models share the same UVB after reionization, they tend to converge to the same ${\lambda }_{{\rm{P}}}$ values at lower redshifts. Although the pressure smoothing scale of the IGM depends on the full thermal history, the thermal memory of past reionization events eventually fades as the gas evolves toward lower redshifts.

6.1. Heating during Hydrogen Reionization

Figure 6 shows the results from simulations (MiddleR-Hcold, MiddleR, MiddleR-Hwarm, and MiddleR-Hhot) using different total input heat during ${\rm{H}}\,{\rm{I}}$ reionization, i.e., different ${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}$ ($1\times {10}^{4}$, $2\times {10}^{4}$, $3\times {10}^{4}$, and $4\times {10}^{4}$ K, respectively). All these simulations share the same ${\rm{H}}\,{\rm{I}}$ ionization history (middle reionization, ${z}_{\mathrm{reion},{\rm{H}}{\rm{I}}}=8.3$) and exactly the same $\mathrm{He}\,{\rm{II}}$ ionization history and heating. The left panel of Figure 6 shows the evolution of γ and T₀ for these simulations. As expected, they differ significantly during ${\rm{H}}\,{\rm{I}}$ reionization, due to the different heat input applied in them. At lower redshifts, $z\lt {z}_{\mathrm{reion},{\rm{H}}{\rm{I}}}$, all these simulations share the same HM12 photoionization and photoheating rates; hence, eventually T₀ and γ thermal parameters tend to converge to the same values. This shows again that these thermal parameters depend more strongly on the instantaneous value of these rates. However, it is very interesting to remark that this convergence is not immediate, but that it takes some time for each simulation to converge after reaching the redshift at which they all have exactly the same rates (McQuinn & Upton Sanderbeck 2016). In any case, our simulations show that γ and T₀ have little sensitivity to the details of ${\rm{H}}\,{\rm{I}}$ reionization after $z\lt 5$.

Zoom In Zoom Out Reset image size

In the right panels of Figure 6 we show the evolution of the temperature at optimal density ($T({{\rm{\Delta }}}_{{\bigstar}})$; upper panel), the curvature ($\langle | \kappa | \rangle$; middle panel), and the pressure smoothing scale (${\lambda }_{{\rm{P}}}$; lower panel). As expected, $T({{\rm{\Delta }}}_{{\bigstar}})$ follows the same trend as T₀ and γ (see discussion above). However, the curvature statistics and the pressure smoothing scale for these simulations clearly show a different behavior at redshifts above $z\gtrsim 3$, while T₀, γ, and $T({{\rm{\Delta }}}_{{\bigstar}})$ have already forgotten reionization at $z\sim 5$. Simulations with a higher heat input at high redshift show higher-pressure smoothing scale values even at lower redshift, due to the dependence of this parameter on the full thermal history. Comparing this result with Figure 5, we can see that the ionization history and the total heat input are degenerate in terms of the pressure smoothing scale. That is, the earlier that H i reionization injects heat into the IGM, the larger the gas pressure scale, ${\lambda }_{P}$. However, a later but hotter (larger heat input) H i reionization also results in a larger pressure smoothing scale. This cautions one about the interpretation of curvature-based measurements of the IGM temperature (Becker et al. 2011) at $z\gt 3.5$, since the middle right panel of Figure 6 clearly illustrates that the curvature has a strong dependence on thermal history, even when the instantaneous temperature is the same in all models.

6.2. Heating during Helium Reionization

Finally, we discuss simulations for which we modified only the total heat input during $\mathrm{He}\,{\rm{II}}$ reionization. These are MiddleR-noHe, MiddleR-Hecold, MiddleR, MiddleR-Hewarm, and MiddleR-Hehot, and the specific values of ${\rm{\Delta }}{T}_{\mathrm{He}{\rm{II}}}$ used in them were 0 (no $\mathrm{He}\,{\rm{II}}$ reionization), $1\times {10}^{4}$, $1.5\times {10}^{4}$, $2\times {10}^{4}$, and $3\times {10}^{4}$ K, respectively. Apart from this, these simulations share exactly the same $\mathrm{He}\,{\rm{I}}$ ionization history and also the same ${\rm{H}}\,{\rm{I}}$ ionization and photoheating rates (middle reionization, ${z}_{\mathrm{reion},{\rm{H}}{\rm{I}}}=8.30$). Therefore, it is not surprising that the evolution of γ, T₀, and ${T}_{{{\rm{\Delta }}}_{{\bigstar}}}$ thermal parameters at redshift above $z\gtrsim 5$ is exactly the same for all simulations (see upper left, lower left, and upper right panels of Figure 7). It is only when $\mathrm{He}\,{\rm{II}}$ reionization starts heating the IGM that these simulations differ. In particular, runs with a higher total heat input result in a much steeper rise when $\mathrm{He}\,{\rm{II}}$ reionization commences in T₀, accompanied by a commensurate fall when $\mathrm{He}\,{\rm{II}}$ reionization is completed. With the slope of the temperature–density relation, γ, the effect is the opposite. The curvature statistic and the pressure smoothing scale also illustrate the effect of the different heating during $\mathrm{He}\,{\rm{II}}$ reionization. Simulations with a larger late heat input due to $\mathrm{He}\,{\rm{II}}$ reionization give rise to a larger pressure smoothing scale. Notice that the pressure smoothing scales of these models begin to diverge at $z\lt 3$, once $\mathrm{He}\,{\rm{II}}$ reionization has already completed, as there is a delay before the effect propagates. This delay is due to the dynamical time that it takes the gas to respond to temperature changes at the Jeans scale (i.e., the sound-crossing time), which, as discussed above for IGM densities, is close to the Hubble time.

Zoom In Zoom Out Reset image size

The simplest possible Lyα flux statistic is the mean transmitted flux $\langle {F}_{{\rm{H}}{\rm{I}}}\rangle$, or equivalently, the effective optical depth ${\tau }_{{\rm{H}}{\rm{I}}}=\mathrm{log}\langle {F}_{{\rm{H}}{\rm{I}}}\rangle$. It is commonly the case that simulations do not recover the observed mean flux, but that simulated fluxes are rescaled to match the observed mean. This rescaling is often understood as equivalent to adjusting the specific intensity of the ${\rm{H}}\,{\rm{I}}$ photoionization rate used in the simulation and is justified based on how poorly constrained the ionizing background is. Notice, however, that this rescaling is generally done directly in redshift space. Lukić et al. (2015) conducted a detailed study of the effect that this rescaling can have on different Lyα statistics. They found that for the large rescalings—those where optical depth has to be rescaled by a factor of 2 or more—the error on flux power spectrum is a few percent. Also, the larger the rescaling is, the larger the error that is introduced. The rescaling error is therefore small, but not negligible, and most importantly, it is puzzling why one should continue to run simulations that systematically produce mean flux values excluded by observations at the few sigma level and continue to compensate by rescaling the optical depth by a factor of few, as is currently required with the HM12 or FG09 UVB tables. For this reason we wish to correct for this error in our new UVB models by renormalizing the input ${\rm{H}}\,{\rm{I}}$ photoionization rate so that the post-processing correction will be minimal at all relevant redshifts. We want to emphasize here that the goal of this step is not to remove the need for future rescaling of the optical depth in simulations, but only to ensure that mean fluxes obtained by simulations are roughly consistent with current observations, therefore removing the need for large rescalings. Trying to do better than that would be a pointless exercise, as the change in cosmological parameters, as well as having different resolution or box size, will anyway change the mean flux at a few percent level. The resolution of the simulations discussed here is sufficient for obtaining mean flux converged at a few percent level (Lukić et al. 2015). In fact, we have run simulations of our LateR, MiddleR, and EarlyR UVB models using a larger box size, ${L}_{\mathrm{box}}=40$ Mpc h⁻¹, and 2048³ resolution elements to confirm that this is indeed the case.

Observational constraints coming from quasar absorption lines (Fan et al. 2006; Becker et al. 2007; Kirkman et al. 2007; Faucher-Giguère et al. 2008; Becker & Bolton 2013) show that the mean flux smoothly evolves from about 0.2 at z = 5 to about 0.9 at z = 2, as expansion gradually lowers the density and the UVB intensity slowly increases. Figure 8 shows a compilation of these observations using different symbols with $1\sigma$ error bars. We also plot suggested fits to the mean flux evolution by various authors (Fan et al. 2006; Kim et al. 2007; Viel et al. 2013a). However, none of these fits do a particularly good job of describing the full evolution of the mean flux, and we therefore opt for our own fit using all observational data points between $0.2\lt z\lt 5.85$. We found that the functional form

$$\begin{eqnarray}&&{\tau }_{{\rm{H}}{\rm{I}}}=A\times {e}^{(B\times \sqrt{z})}\end{eqnarray} \tag{ 13 }$$

provides an optimal fit, with A = 0.00126 and B = 3.294 as the best-fit parameters. This fit is also shown in Figure 8 as a solid black line.

Zoom In Zoom Out Reset image size

The left panel of Figure 9 shows the ${\rm{H}}\,{\rm{I}}$ mean flux evolution in simulations using different well-known UVB models. It is clear that several of them significantly underpredict the observed mean flux at all redshifts. It is also worth pointing out that, somewhat coincidentally, the HM01 UVB model is doing a good job in recovering the observed mean flux. Notice, however, that since our thermal history is different from theirs, we cannot just assume their photoionization rates. Hence, we have modified the ${\rm{H}}\,{\rm{I}}$ photoionization rate in our models after reionization so that the ${\rm{H}}\,{\rm{I}}$ mean flux in the simulations matches the fit to current observational constraints. Before ${\rm{H}}\,{\rm{I}}$ reionization, when we apply our new methodology, this does not apply. We also modified the ${\rm{H}}\,{\rm{I}}$ and $\mathrm{He}\,{\rm{I}}$ photoheating rate by the same factor so that the heat input at these redshifts is conserved, i.e., ${n}_{{\rm{H}}{\rm{I}},\mathrm{old}}{\dot{q}}_{{\rm{H}}{\rm{I}},\mathrm{old}}={n}_{{\rm{H}}{\rm{I}},\mathrm{new}}{\dot{q}}_{{\rm{H}}{\rm{I}},\mathrm{new}}$, and therefore we get exactly the same thermal histories. We have confirmed that this is the case by comparing the evolution of the thermal parameters in the new models versus the old ones, and we show the mean flux evolution of the MiddleR simulation as a dashed line in Figure 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Regarding the effect of changing the thermal history, we have compared the differences between the simulations in which we change the heat input due to $\mathrm{He}\,{\rm{II}}$ reionization. We found differences in the ${\rm{H}}\,{\rm{I}}$ mean flux up to $\sim 10 \%$ and $\sim 15 \%$ between our our fiducial MiddleR model and the two most extreme simulations MiddleR-Hehot and MiddleR-noHe, respectively. This is because through the recombination factor ($\alpha \propto {T}_{0}^{-0.7}$) the neutral fraction has a sensitivity to the gas temperature, not just the photoionization: ${n}_{{\rm{H}}{\rm{I}}}\propto {{\rm{\Gamma }}}_{{\rm{H}}\,{\rm{I}}}^{-1}{T}_{0}^{-0.7}{{\rm{\Delta }}}^{0.7* (\gamma -1)}$. This maximum difference corresponds to $z\sim 3$, the redshift at which the thermal parameters between these simulations are most different. Of course, simulations using thermal histories that deviate even further from these would increase these differences. We have confirmed that mean flux differences due to current uncertainties in the cosmological parameters have a much weaker effect on the mean flux than any of the above systematics in the UVB models described above (see Appendix C for a full discussion on cosmological parameters).

Finally, in the right panel of Figure 9 we show the most recent observations of the $\mathrm{He}\,{\rm{II}}$ transmission (Worseck et al. 2016) and compare them again with the mean flux derived from our hydrodynamical simulations using standard UVB models (HM96, HM01, FG09, HM12).²² We show the mean fluxes, $\langle {F}_{\mathrm{He}{\rm{II}}}\rangle ={e}^{-{\tau }_{\mathrm{He}{\rm{II}}}}$, and not the optical depths in order to focus on the low-redshift results, $z\lesssim 2.7$, where observations seem to indicate that $\mathrm{He}\,{\rm{II}}$ reionization has already finished. Results at higher redshifts are not that conclusive and may indicate that this reionization happened much more slowly than has been assumed (see Worseck et al. 2016, for a detailed discussion). Therefore, we want to focus here just on the low-redshift values, where reionization is completed and our method should be valid. For these it seems that FG09 $\mathrm{He}\,{\rm{II}}$ photoionization rates are doing the best job in reproducing the observations. For this reason we decided to use the $\mathrm{He}\,{\rm{II}}$ photoionization and photoheating rates of this model in our new UVB models after $\mathrm{He}\,{\rm{II}}$ reionization.

In this section we elaborate on the convergence of the results presented in this work, compare them with recent efforts done in the field, and finally discuss their implications for galaxy formation simulations.

8.1. Resolution and Convergence

8.2. Comparison to Previous Work

8.3. Implications for Galaxy Formation

As mentioned in the introduction, an ionizing UVB inhibits gas accretion and photoevaporates gas from the shallow potential wells of low-mass dwarf galaxies. This effect due to gas being heated up by photoionization can result in negative feedback: suppressing star formation inside reionized regions, thus impeding their continued growth (see, e.g., Rees 1986; Efstathiou 1992). In our current picture of galaxy formation this feedback is considered to be a good candidate for resolving the "missing satellite" problem (Klypin et al. 1999; Moore et al. 1999), which arises because in the cold dark matter (ΛCDM) framework the number of simulated dark matter subhalos is much larger than the number of observed dwarf galaxies (Babul & Rees 1992; Bullock et al. 2001). This picture seems consistent with recent observations that have seen uniformly old stellar populations in ultrafaint galaxies (Brown et al. 2014).

This topic has been a very active area of research in the past years, and using cosmological simulations is key to trying to understand these effects in their full context (Quinn et al. 1996; Thoul & Weinberg 1996; Gnedin 2000; Hoeft et al. 2006; McQuinn et al. 2007; Okamoto et al. 2008; Noh & McQuinn 2014; Benítez-Llambay et al. 2015). More recent simulations with much better resolution and a more complete feedback model also seem to point in this direction (Oñorbe et al. 2015; Wheeler et al. 2015). In this context, the work by Simpson et al. (2013) is particularly relevant. Using an enzo high-resolution cosmological hydrodynamical simulation, they found that turning on the UVB at z = 7.0 versus z = 8.9 resulted in an order-of-magnitude change in the final stellar mass of a 10⁹ ${M}_{\odot }$ dark matter halo. In addition, some very interesting constraints are starting to come from radiative transfer hydrodynamical simulations. Wise et al. (2014) find that very faint galaxies (${M}_{\mathrm{UV}}\sim -6$, ${M}_{* }\sim {10}^{3.59}$ in halos of ${M}_{h}=1\times {10}^{7}$) will still form at high redshift and contribute a significant amount to the ionizing photon budget during cosmic reionization. However, in order to avoid overproducing the observed abundance of classical satellites of the Milky Way, studies based on dark-matter-only simulations argue for a critical mass closer to $\sim {10}^{9}$ ${M}_{\odot }$ (Madau et al. 2008; Boylan-Kolchin et al. 2014).

For all this it is clear that the use of current standard UVB models that reionize and, more importantly here, heat the IGM at a much higher redshift than was desired will have a strong impact on the results of galaxy formation hydrodynamical simulations. For example, simulations that use the HM12 UVB background start spuriously heating up the IGM at $z\sim 15$ (see Figure 1). We can now perform a simple first analytical calculation to get an approximate idea of this effect. In linear theory, the instantaneous cosmological mass cooling threshold of the neutral IGM (sometimes also referred to as the Jeans mass) is given by (see, e.g., Iliev et al. 2007)

$$\begin{eqnarray}{M}_{\min } & = & 3.9\times {10}^{9}{M}_{\odot }{\left(\displaystyle \frac{{T}_{\mathrm{IGM}}}{{10}^{4}{\rm{K}}}\right)}^{3/2}{\left(\displaystyle \frac{1+z}{10}\right)}^{3/2}\\ & & \times {\left(\displaystyle \frac{{{\rm{\Omega }}}_{{\rm{b}}}{h}^{2}}{0.0223}\right)}^{-3/5}{\left(\displaystyle \frac{{{\rm{\Omega }}}_{{\rm{m}}}{h}^{2}}{0.15}\right)}^{-1/2},\end{eqnarray} \tag{ 14 }$$

where ${T}_{\mathrm{IGM}}$ is the temperature of the IGM. As has been found in the simulations mentioned above, the actual cooling mass differs somewhat from this instantaneous Jeans mass since the mass scale on which baryons succeed in collapsing out of the IGM along with the dark matter must be determined, even in linear theory, by integrating the differential equation for perturbation growth over time for the evolving IGM (Gnedin & Hui 1998). In fact, there is no single mass above which a collapsing halo retains all its gas and below which the gas does not collapse with the dark matter. Instead, simulations show that the cooled gas fraction in halos decreases gradually with decreasing halo mass (e.g., Okamoto et al. 2008). Using Equation (14) can give us a first estimate of the possible effects of using different UVB models. Figure 10 shows the results of this equation when we use the different thermal histories of the HM12, FG09, LateR, MiddleR, and EarlyR models. The differences between the models are quite significant.

Different studies of the reionization redshift of collapse structures have shown that the median reionization redshift of halos moves to lower redshift and the scatter increases (Weinmann et al. 2007; Alvarez et al. 2009; Li et al. 2014). Although both parameters depend substantially on the details of reionization (see, e.g., Ocvirk et al. 2013), massive halos, $\gt {10}^{15}$ ${M}_{\odot }$, can be reionized significantly earlier than the average region in the universe (${\rm{\Delta }}{z}_{\mathrm{reion},50 \% }\sim 2$; see Figure 2 of Li et al. 2014). For Milky Way halos (∼10¹² ${M}_{\odot }$, ) or dwarf galaxies (∼10¹⁰ ${M}_{\odot }$) the typical reionization redshift is expected to be much lower (${\rm{\Delta }}{z}_{\mathrm{reion},50 \% }\sim 0.5$ and ${\rm{\Delta }}{z}_{\mathrm{reion},50 \% }\lesssim 0$, respectively), making FG09 and HM12 models not optimal for these studies. This is even more true if one considers the last constraints on ${\tau }_{{\rm{e}}}$ from Planck Collaboration et al. (2016b). By overestimating the heat at high z, we are not only overestimating the pressure smoothing scale at lower redshifts but also overestimating the effect of the UVB in galaxy formation and evolution.²³ For all these reasons, a detailed review of the results of galaxy formation simulations using new consistent UVB models that fulfill all observational constraints, like the ones developed here, is certainly needed. This is another reason why we make the models presented in this work publicly available to the community (see Appendix D) and encourage colleagues to adopt them in future galaxy formation work and revisit previous calculations.

In this paper we have presented results from optically thin cosmological hydrodynamical simulations using the Nyx code (Almgren et al. 2013; Lukić et al. 2015). As commonly done in multiple IGM and galaxy formation studies, the UV background is modeled as a uniform and isotropic field that evolves with redshift. Operationally, the UVB determines the photoioinization and photoheating rates of ${\rm{H}}\,{\rm{I}}$, $\mathrm{He}\,{\rm{I}}$, and $\mathrm{He}\,{\rm{II}}$, which are inputs to the code.

We have demonstrated that when canonical models of the UVB, like that of HM12, are used in hydrodynamical simulations, the ionization of the IGM and, more importantly, the concomitant heating occur far too early, inconsistent with the reionization histories calculated by the respective authors, and in violation of current observational constraints on reionization. We argue that this results from the fact that these models dramatically overestimated the mean free path of ionizing photons at root at $z\gt 5$, resulting from the blind extrapolation of a model fit to lower-redshift ($z\lt 5$) measurements. As a result, the amplitudes of the photoionization and photoionization heating rates are far too high at $z\gt 6$. This premature heating spuriously heats the IGM to $\sim {10}^{4}$ K by $z\sim 13$, and because the IGM gas pressure smoothing scale depends on the full thermal history, it produces an erroneously large pressure smoothing scale at nearly all redshifts. We argue that a correct and consistent model of the reionization and thermal history is crucial for obtaining the correct pressure scale in simulations, which is necessary for interpreting Lyα forest statistics at $z\lt 6$—not doing so can bias estimates of the thermal state of the IGM. We also discussed the implications of this spurious early heating on galaxy formation simulations.

Motivated by these issues, we have developed a new method to generate UVB models for hydrodynamical simulations that allow one to self-consistently simulate different reionization models. We implement this by volume-averaging the photoionization and energy equations. In the model, each reionization event is defined by the ionization history with redshift and the total heat input of the reionization event. In this sense our new models provide a very promising tool to explore the parameter space of possible ionization and thermal histories. In this work, we investigated models in which we changed the redshift at which ${\rm{H}}\,{\rm{I}}$ reionization ends and the amount of heat input associated with both ${\rm{H}}\,{\rm{I}}$ and $\mathrm{He}\,{\rm{II}}$ reionization. We studied the effect of these changes on the thermal history of the IGM, in particular the temperature at mean density, T₀, the slope of the temperature–density relation, γ, the temperature at the optimal density probed by curvature measurements, $T({{\rm{\Delta }}}_{{\bigstar}})$, and the gas pressure smoothing scale, ${\lambda }_{P}$. We have shown how important the degeneracies between these parameters can be in order to derive the thermal parameters of the IGM using the curvature Lyα statistic. These rates have also been corrected to improve the agreement with measurements of the average ${\rm{H}}\,{\rm{I}}$ and $\mathrm{He}\,{\rm{II}}$ transmission after reionization.

The UVB plays a fundamental role in determining the star formation of the first galaxies and their evolution by not only setting the minimum halo mass able to form stars but also regulating the gas accretion from the IGM into more massive halos. Previous studies utilizing UVB models that suffer from the spurious early heating described in this paper have thus overestimated the effect of this photoheating feedback and the resulting suppression of star formation at high redshift. We therefore argue that galaxy formation simulations should be revisited using our new UVB models.

We make our new UVB models publicly available so that the community can better explore the consequences and effects of different ionization and thermal histories in all types of hydrodynamical cosmological simulations. Tables with the photoionization and photoheating rates of the new models can be found in Appendix D, which are in the standard "TREECOOL" file format, ready and easy to use with most cosmological codes, including gadget, arepo, and gizmo. We encourage anyone interested in implementing some other specific UVB model to contact the authors. We will also be happy to provide help incorporating these models into other codes by request.

Like all optically thin simulations, our approach misses UVB fluctuations that could produce scatter in reionization times and temperatures between different regions of the universe (e.g., Abel & Haehnelt 1999; Meiksin & White 2004; Davies & Furlanetto 2014, 2016; Gontcho A Gontcho et al. 2014; Pontzen 2014; D'Aloisio et al. 2015; Malloy & Lidz 2015). Of course, an immediate solution to this problem in optically thin simulations will be to add some extra dependence on a specific property of each resolution element in the simulation (e.g., density, distance from a halo that could host a galaxy/quasar, etc.). However, the computational challenge is to try to do this without making the simulation prohibitively expensive. In this sense, an interesting solution that is worth exploring will be to assign a specific reionization redshift to each resolution element of the simulation using, for example, an excursion set formalism (e.g., Furlanetto et al. 2004). In this picture, a resolution element will not see the UVB background until its reionization redshift. We plan to pursue this idea in the near future. Another very valuable piece of information that would improve current optically thin hydrodynamical simulations would be if future radiative transfer simulations (e.g., Gnedin 2014; So et al. 2014; Norman et al. 2015; Pawlik et al. 2015; Ocvirk et al. 2016) were to make the probability distribution function of their photoionization rates publicly available (and perhaps also its dependence on density), and not just the evolution of the mean and/or median values.

Finally, our new parameterization for the heating and ionization produced by the UVB allows us to explore a broader range of reionization models, as well as any other physical scenarios that could alter the thermal history of the IGM. This will allow us to better test the effect of such models in simulations of galaxy formation and the IGM. These models could include Population III stars (Manrique et al. 2015), X-ray pre-heating coming from starburst galaxies, supernova remnants or miniquasars (Oh 2001; Glover & Brand 2003; Madau et al. 2004; Furlanetto 2006), dark matter annihilation or decay (Liu et al. 2016) or cosmic rays (Samui et al. 2005), from the intergalactic absorption of blazar TeV photons (Chang et al. 2012; Puchwein et al. 2012) or from broadband intergalactic dust absorption (Inoue & Kamaya 2008). We expect more detailed studies on these physically motivated models in the future.

We thank the members of the ENIGMA group at the Max Planck Institute for Astronomy (MPIA) for helpful discussions. J.F.H. acknowledges generous support from the Alexander von Humboldt foundation in the context of the Sofja Kovalevskaja Award. The Humboldt foundation is funded by the German Federal Ministry for Education and Research. Z.L. was 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. The ASCR Leadership Computing Challenge (ALCC) program has provided NERSC allocation under the project "Cosmic Frontier Computational End-Station." This work made extensive use of the NASA Astrophysics Data System and of the astro-ph preprint archive at arXiv.org.

In this paper we have presented a new way of obtaining effective photoionization and photoheating rates for different reionization models. These can be used in optically thin hydrodynamical simulations to account for the emission of galaxies and quasars. Here we provide detailed derivation of these rates for future reference.

We remind the reader that the new effective rates are only computed during reionization, i.e., while $\langle {x}_{{\rm{H}}{\rm{II}}}\rangle \lt 1$. After reionization, our models use the rates from common UVB models (e.g., FG09, HM12; see Section 7). New values of the photoionization rates during reionization are forced to never exceed the values of the model plugged in after reionization. This is done to guarantee no numerical artifacts in the limit where reionization is almost complete, $\langle {x}_{{\rm{H}}{\rm{II}}}\rangle \sim 1$, where our new rates are not well defined.

A.1. Volume-averaged Optically Thin Ionization Equations

In the context of optically thin hydrodynamical codes, it is often assumed that the gas is of the primordial chemical composition, where the resulting reaction network includes six atomic species: ${\rm{H}}\,{\rm{I}}$, ${\rm{H}}\,{\rm{II}}$, $\mathrm{He}\,{\rm{I}}$, $\mathrm{He}\,{\rm{II}}$, $\mathrm{He}\,{\rm{III}}$, and e⁻. Codes generally evolve those species under the assumption of ionization equilibrium (see, however, Gnat & Sternberg 2007; Vasiliev 2011; Oppenheimer & Schaye 2013; Richings et al. 2014a, 2014b; Puchwein et al. 2015, for nonequilibrium treatments). The resulting system of algebraic equations is

$$\begin{eqnarray}&&({{\rm{\Gamma }}}_{{\rm{e}},{\rm{H}}{\rm{I}}}{n}_{{\rm{e}}}+{{\rm{\Gamma }}}_{\gamma ,{\rm{H}}{\rm{I}}}){n}_{{\rm{H}}{\rm{I}}}={\alpha }_{{\rm{r}},{\rm{H}}{\rm{II}}}{n}_{{\rm{e}}}{n}_{{\rm{H}}{\rm{II}}}\\ &&({{\rm{\Gamma }}}_{{\rm{e}},\mathrm{He}{\rm{I}}}{n}_{{\rm{e}}}+{{\rm{\Gamma }}}_{\gamma ,\mathrm{He}{\rm{I}}}){n}_{\mathrm{He}{\rm{I}}}=({\alpha }_{{\rm{r}},\mathrm{He}{\rm{II}}}+{\alpha }_{{\rm{d}},\mathrm{He}{\rm{II}}}){n}_{{\rm{e}}}{n}_{\mathrm{He}{\rm{II}}}\\ &&[{{\rm{\Gamma }}}_{\gamma ,\mathrm{He}{\rm{II}}}+({{\rm{\Gamma }}}_{{\rm{e}},\mathrm{He}{\rm{II}}}+{\alpha }_{{\rm{r}},\mathrm{He}{\rm{II}}}+{\alpha }_{{\rm{d}},\mathrm{He}{\rm{II}}}){n}_{{\rm{e}}}]{n}_{\mathrm{He}{\rm{II}}}\\ &&\quad ={\alpha }_{{\rm{r}},\mathrm{He}{\rm{III}}}{n}_{{\rm{e}}}{n}_{\mathrm{He}{\rm{III}}}+({{\rm{\Gamma }}}_{{\rm{e}},\mathrm{He}{\rm{I}}}{n}_{{\rm{e}}}+{{\rm{\Gamma }}}_{\gamma ,\mathrm{He}{\rm{I}}}){n}_{\mathrm{He}{\rm{I}}}.\end{eqnarray} \tag{ 15 }$$

In addition, there are three closure equations for the conservation of charge and hydrogen and helium abundances. Radiative recombination (${\alpha }_{{\rm{r}},{\rm{X}}}$), dielectronic recombination (${\alpha }_{{\rm{d}},{\rm{X}}}$), and collisional ionization (${{\rm{\Gamma }}}_{{\rm{e}},{\rm{X}}}$) rates are strongly dependent on the temperature, which itself depends on the ionization state through the mean mass per particle μ,

$$\begin{eqnarray}&&T=\displaystyle \frac{2}{3}\displaystyle \frac{{m}_{p}}{{k}_{{\rm{B}}}}\mu \ {e}_{\mathrm{int}},\end{eqnarray} \tag{ 16 }$$

where m_p is the mass of a proton, ${k}_{{\rm{B}}}$ is the Boltzmann constant, and ${e}_{\mathrm{int}}$ is the internal thermal energy per mass of the gas. For a gas composed of only hydrogen and helium, μ is related to the number density of free electrons relative to hydrogen by $\mu =(1+4\chi )/[1+\chi +({n}_{{\rm{e}}}/{n}_{{\rm{H}}})]$. The reaction network equations are iteratively solved together with the ideal gas equation of state, $p=2\rho {e}_{\mathrm{int}}/3$, to determine the temperature and equilibrium distribution of species. Above and throughout this paper we have assumed an adiabatic index of 5/3.

In order to produce consistent UVB models that reliably reproduce different reionization histories, we have derived the volume-averaged version of the ionization equilibrium equations presented in Equation (15). We start with the ${\rm{H}}\,{\rm{I}}$ reionization and derive here in detail the equation for ${\rm{H}}\,{\rm{I}}$ photoionization. We will address the $\mathrm{He}$ single and double reionization afterward, as the method is very similar. We start by doing the volume averages of Equation (15) for ${\rm{H}}\,{\rm{I}}$,

$$\begin{eqnarray}&&\langle {{\rm{\Gamma }}}_{{\rm{e}},{\rm{H}}{\rm{I}}}{n}_{{\rm{e}}}{n}_{{\rm{H}}{\rm{I}}}\rangle +\langle {{\rm{\Gamma }}}_{\gamma ,{\rm{H}}{\rm{I}}}{n}_{{\rm{H}}{\rm{I}}}\rangle =\langle {\alpha }_{{\rm{r}},{\rm{H}}{\rm{II}}}{n}_{{\rm{e}}}{n}_{{\rm{H}}{\rm{II}}}\rangle .\end{eqnarray} \tag{ 17 }$$

The first thing is that collisional ionization terms, ${{\rm{\Gamma }}}_{{\rm{e}}}$, are mainly relevant in shocks at high temperatures and densities. They will have a negligible effect on this volume-averaged calculation, so we can discard them. The important point here is that each term in this equation is nonlinear; hence, in principle, it is not possible to compute their volume averages unless the cross-correlation of the abundances of ${n}_{{\rm{H}}{\rm{I}}}$ and ${n}_{{\rm{H}}{\rm{II}}}$ with each other, as well as with the radiation and temperature fields, is known (since ${{\rm{\Gamma }}}_{{\rm{e}},{\rm{H}}{\rm{I}}}$ and ${\alpha }_{{\rm{r}},{\rm{H}}{\rm{II}}}$ are temperature dependent). A convenient way to encapsulate this unknown information is using correction factors. With all this we can write the volume-averaged equivalent of Equation (15) as

$$\begin{eqnarray}&&{C}_{\gamma ,{\rm{H}}{\rm{I}}}\langle {{\rm{\Gamma }}}_{\gamma ,{\rm{H}}{\rm{I}}}\rangle \langle {n}_{{\rm{H}}{\rm{I}}}\rangle ={C}_{{\rm{r}},{\rm{H}}{\rm{II}}}\langle {\alpha }_{{\rm{r}},{\rm{H}}{\rm{II}}}\rangle \langle {n}_{{\rm{e}}}\rangle \langle {n}_{{\rm{H}}{\rm{II}}}\rangle ,\end{eqnarray} \tag{ 18 }$$

where the correction factors are defined as

$$\begin{eqnarray}{C}_{\gamma ,{\rm{H}}{\rm{I}}} & = & \displaystyle \frac{\langle {{\rm{\Gamma }}}_{\gamma ,{\rm{H}}{\rm{I}}}{n}_{{\rm{H}}{\rm{I}}}\rangle }{\langle {{\rm{\Gamma }}}_{\gamma ,{\rm{H}}{\rm{I}}}\rangle \langle {n}_{{\rm{H}}{\rm{I}}}\rangle }\\ {C}_{{\rm{r}},{\rm{H}}{\rm{II}}} & = & \displaystyle \frac{\langle {\alpha }_{{\rm{r}},{\rm{H}}{\rm{II}}}{n}_{{\rm{e}}}{n}_{{\rm{H}}{\rm{II}}}\rangle }{\langle {\alpha }_{{\rm{r}},{\rm{H}}{\rm{II}}}\rangle \langle {n}_{{\rm{e}}}\rangle \langle {n}_{{\rm{H}}{\rm{II}}}\rangle }.\end{eqnarray} \tag{ 19 }$$

Then, we can rewrite Equation (18) as

$$\begin{eqnarray}&&\langle {{\rm{\Gamma }}}_{\gamma ,{\rm{H}}{\rm{I}}}\rangle =\displaystyle \frac{{C}_{{\rm{r}},{\rm{H}}{\rm{II}}}}{{C}_{\gamma ,{\rm{H}}{\rm{I}}}}\displaystyle \frac{\langle {\alpha }_{{\rm{r}},{\rm{H}}{\rm{II}}}\rangle \langle {n}_{{\rm{e}}}\rangle \langle {n}_{{\rm{H}}{\rm{II}}}\rangle }{\langle {n}_{{\rm{H}}{\rm{I}}}\rangle }.\end{eqnarray} \tag{ 20 }$$

Now, in order to obtain an ${\rm{H}}\,{\rm{I}}$ photoionization rate, we will make use of the following assumptions: (1) We assume that the $\mathrm{He}\,{\rm{I}}$ reionization occurs perfectly coupled with the ${\rm{H}}\,{\rm{I}}$ reionization process. This is a very common assumption in reionization models because ionization of hydrogen and the first ionization of helium require photons with similar energies (13.6 eV and 24.6 eV, respectively). Therefore, the same physical process responsible for ${\rm{H}}\,{\rm{I}}$ reionization must be also liable for the $\mathrm{He}\,{\rm{I}}$ reionization. We also consider that the helium second ionized state number density is negligible during ${\rm{H}}\,{\rm{I}}$ and $\mathrm{He}\,{\rm{I}}$ reionization, ${n}_{\mathrm{He}{\rm{III}}}\sim 0$, so $\langle {x}_{\mathrm{He}{\rm{III}}}\rangle =0$. This is a correct assumption for all the models discussed in this work but could be easily modified in the future if needed. (2) The evolution of number density of free electrons, ${n}_{{\rm{e}}}$, can be approximated by ${n}_{{\rm{e}}}={n}_{{\rm{H}}{\rm{II}}}+{n}_{\mathrm{He}{\rm{II}}}+2{n}_{\mathrm{He}{\rm{III}}}$. We can rewrite this equation as a function of the hydrogen density and the ionized fractions: ${n}_{{\rm{e}}}={n}_{{\rm{H}}}[(1+\chi ){x}_{{\rm{H}}{\rm{II}}}+\chi {x}_{\mathrm{He}{\rm{III}}}]$, where $\chi ={Y}_{{\rm{p}}}/4{X}_{{\rm{p}}}$ and we have assumed again that $\mathrm{He}\,{\rm{II}}$ reionization follows that of ${\rm{H}}\,{\rm{II}}$. The volume-averaged value can be written as

$$\begin{eqnarray}&&\langle {n}_{{\rm{e}}}\rangle =\langle {n}_{{\rm{H}}}\rangle [{C}_{{x}_{{\rm{H}}{\rm{II}}}}(1+\chi )\langle {x}_{{\rm{H}}{\rm{II}}}\rangle +{C}_{{x}_{\mathrm{He}{\rm{III}}}}\chi \langle {x}_{\mathrm{He}{\rm{III}}}\rangle ],\end{eqnarray} \tag{ 21 }$$

where C_x factors stand for the correction factors defined as ${C}_{{x}_{{\rm{H}}{\rm{I}}}}\langle {n}_{{\rm{H}}}\rangle \langle {x}_{{\rm{H}}{\rm{I}}}\rangle =\langle {n}_{{\rm{H}}{\rm{I}}}\rangle$, ${C}_{{x}_{{\rm{H}}{\rm{II}}}}\langle {n}_{{\rm{H}}}\rangle \langle {x}_{{\rm{H}}{\rm{II}}}\rangle =\langle {n}_{{\rm{H}}{\rm{II}}}\rangle$, etc. As we are assuming that $\mathrm{He}\,{\rm{II}}$ is not relevant during ${\rm{H}}\,{\rm{I}}$ reionization, we can discard the second term inside parentheses for our current derivation. Finally, putting together Equation (20) with Equation (21), we arrive at

$$\begin{eqnarray}&&\langle {{\rm{\Gamma }}}_{\gamma ,{\rm{H}}{\rm{I}}}\rangle (z)={C}_{{\rm{H}}{\rm{II}}}\langle {n}_{{\rm{H}}}\rangle (z){\alpha }_{{\rm{r}},{\rm{H}}{\rm{II}}}(\langle T\rangle )(1+\chi )\displaystyle \frac{\langle {x}_{{\rm{H}}{\rm{II}}}{\rangle }^{2}(z)}{\langle {x}_{{\rm{H}}{\rm{I}}}\rangle (z)},\end{eqnarray} \tag{ 22 }$$

where ${C}_{{\rm{H}}{\rm{II}}}$ encapsulates all correction parameters described above and can be written as

$$\begin{eqnarray}&&{C}_{{\rm{H}}{\rm{II}}}=\displaystyle \frac{{C}_{{\alpha }_{{\rm{r}},{\rm{H}}{\rm{II}}}}\times {C}_{{\rm{r}},{\rm{H}}{\rm{II}}}\times {C}_{{x}_{{\rm{H}}{\rm{II}}}}^{2}}{{C}_{\gamma ,{\rm{H}}{\rm{I}}}\times {C}_{{x}_{{\rm{H}}{\rm{I}}}}},\end{eqnarray} \tag{ 23 }$$

where ${C}_{{\alpha }_{{\rm{r}},{\rm{H}}{\rm{II}}}}$ is defined as ${C}_{{\alpha }_{{\rm{r}},{\rm{H}}{\rm{II}}}}=\langle {\alpha }_{{\rm{r}},{\rm{H}}{\rm{II}}}\rangle /{\alpha }_{{\rm{r}},{\rm{H}}{\rm{II}}}(\langle T\rangle )$ to clarify that in general the volume-averaged recombination factor is redefined as the recombination factor at a certain mean temperature value. It is customary to use the temperature at mean density, T₀, as this value. In the case of optically thin hydrodynamical simulations a constant photoionization rate is used throughout the whole volume, so ${C}_{\gamma }=1$.

In general, the unknown ratio of the IGM's true recombination rate to its hypothetical rate under the assumption of uniform density and temperature is often referred to as the IGM clumping factor (${C}_{\mathrm{IGM}}={C}_{{\rm{r}},{\rm{H}}{\rm{II}}}$). Notice also that, if for Equation (19) we assume that ${n}_{{\rm{e}}}={n}_{{\rm{H}}{\rm{II}}}$ and that ${\alpha }_{{\rm{r}}}$ is a constant, it is easy to redefine the IGM clumping in a much more familiar form: ${C}_{\mathrm{IGM}}=\langle {n}_{{\rm{H}}\,{\rm{II}}}^{2}\rangle /\langle {n}_{{\rm{H}}{\rm{II}}}{\rangle }^{2}$ (Kohler et al. 2007; Pawlik et al. 2009; Finlator et al. 2012; Kaurov & Gnedin 2014).

We can also obtain analogous equations for the volume-averaged helium photoionization rates from Equation (15) to the one we obtained above, Equation (22), for those of hydrogen. We have just made one extra assumption, which is that ${n}_{\mathrm{He}{\rm{III}}}$ is negligible during the single reionization of helium and that during the double reionization of helium ${n}_{\mathrm{He}{\rm{I}}}$ is negligible. In the context of all the $\mathrm{He}\,{\rm{II}}$ reionization models discussed in this paper this is a fair assumption. For the $\mathrm{He}$ first ionization we obtain

$$\begin{eqnarray}&&\langle {{\rm{\Gamma }}}_{\gamma ,\mathrm{He}{\rm{I}}}\rangle (z)={C}_{\mathrm{He}{\rm{II}}}\langle {n}_{{\rm{H}}}\rangle (z){\alpha }_{{\rm{r}},\mathrm{He}{\rm{II}}}(\langle T\rangle )(1+\chi )\displaystyle \frac{\langle {x}_{\mathrm{He}{\rm{II}}}{\rangle }^{2}(z)}{\langle {x}_{\mathrm{He}{\rm{I}}}\rangle (z)},\end{eqnarray} \tag{ 24 }$$

and ${C}_{\mathrm{He}{\rm{II}}}$ is analogous to ${C}_{{\rm{H}}{\rm{II}}}$,

$$\begin{eqnarray}&&{C}_{\mathrm{He}{\rm{II}}}=\displaystyle \frac{{C}_{{\alpha }_{{\rm{r}},\mathrm{He}{\rm{II}}}}\times {C}_{{\rm{r}},\mathrm{He}{\rm{II}}}\times {C}_{{x}_{\mathrm{He}{\rm{II}}}}^{2}}{{C}_{\gamma ,\mathrm{He}{\rm{I}}}\times {C}_{{x}_{\mathrm{He}{\rm{I}}}}}.\end{eqnarray} \tag{ 25 }$$

For the $\mathrm{He}$ double ionization,

$$\begin{eqnarray}&&\langle {{\rm{\Gamma }}}_{\gamma ,\mathrm{He}{\rm{II}}}\rangle (z)=\displaystyle \frac{1}{{C}_{\gamma ,\mathrm{He}{\rm{I}}}}\langle {n}_{{\rm{H}}}\rangle (z)[{C}_{{x}_{{\rm{H}}{\rm{II}}}}(1+\chi )\\ &&\qquad \quad \quad \quad +\chi {C}_{{x}_{\mathrm{He}{\rm{III}}}}\langle {x}_{\mathrm{He}{\rm{III}}}\rangle (z)]\\ &&\left[{C}_{{\alpha }_{{\rm{r}},\mathrm{He}{\rm{III}}}}\times {C}_{{\rm{r}},\mathrm{He}{\rm{III}}}\displaystyle \frac{{\alpha }_{{\rm{r}},\mathrm{He}{\rm{III}}}(\langle T\rangle )\langle {x}_{\mathrm{He}{\rm{III}}}\rangle (z)}{\langle {x}_{\mathrm{He}{\rm{II}}}\rangle (z)}\right.\\ &&\qquad -{C}_{{\alpha }_{{\rm{r}},\mathrm{He}{\rm{II}}}}\times {C}_{{\rm{r}},\mathrm{He}{\rm{II}}}\times {\alpha }_{{\rm{r}},\mathrm{He}{\rm{II}}}(\langle T\rangle )].\end{eqnarray} \tag{ 26 }$$

In order to describe the reionization history of one model, we will need to specify two functions, the evolution of ${\rm{H}}\,{\rm{II}}$ and $\mathrm{He}\,{\rm{III}}$ volume-averaged ionization fractions, $\langle {x}_{{\rm{H}}{\rm{II}}}\rangle (z)$ and $\langle {x}_{\mathrm{He}{\rm{III}}}\rangle (z)$. As $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle +\langle {x}_{{\rm{H}}{\rm{II}}}\rangle =1$ and $\langle {x}_{\mathrm{He}{\rm{I}}}\rangle +\langle {x}_{\mathrm{He}{\rm{II}}}\rangle +\langle {x}_{\mathrm{He}{\rm{III}}}\rangle =1$, we assume that $\mathrm{He}\,{\rm{I}}$ reionization is totally coupled with ${\rm{H}}\,{\rm{I}}$ reionization.

To calculate the average values of the radiative recombination rates (${\alpha }_{{\rm{r}},{\rm{i}}}$), which depend on temperature, we need to specify the evolution of the volume-averaged temperature, $\langle T\rangle (z)$, in our model. To do this, we first need to define the total heat input produced by each reionization event, ${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}$ and ${\rm{\Delta }}{T}_{\mathrm{He}{\rm{II}}}$, which will be two free parameters in our reionization models. We describe in the next paragraph the recent efforts to calculate these parameters from theoretical models. Thus, we will make another assumption in our model to describe the evolution of the volume-averaged temperature: (3) that the evolution of the volume-averaged temperature can be well approximated with the evolution of the reionization history times the total heat input, i.e., $\langle T{\rangle }_{{\rm{H}}{\rm{I}}}(z)={\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}\times \langle {x}_{{\rm{H}}{\rm{II}}}\rangle (z)$ and $\langle T{\rangle }_{\mathrm{He}{\rm{II}}}(z)={10}^{4}+{\rm{\Delta }}{T}_{\mathrm{He}{\rm{II}}}\times {x}_{\mathrm{He}{\rm{III}}}(z)$. This is a very rough estimation of the thermal history, as we are neglecting all cooling, but it serves as a first-order approximation of the gas temperature during reionization in order to compute the different rates that we need to compute Equation (22). We will show below that this assumption is accurate enough for our purposes and for the range of models covered in this work. More elaborate assumptions, including analytic estimations of the Compton cooling, could be implemented in the future.

Finally, to generate the models of this work, we used the following correction factors. For the ${\rm{H}}\,{\rm{I}}$ and $\mathrm{He}\,{\rm{I}}$ photoionization rates we used ${C}_{{\rm{H}}{\rm{II}}}={C}_{\mathrm{He}{\rm{II}}}=1.5$ for $z\geqslant 10$ and ${C}_{{\rm{H}}{\rm{II}}}={C}_{\mathrm{He}{\rm{II}}}=2.0$ for $6\lt z\lt 10$. Notice that we never go below redshift z = 6 when we compute the ${\rm{H}}\,{\rm{I}}$ rates, and that these values are in good agreement with results from radiative transfer simulations (clumping factor, C₁₀₀, of Pawlik et al. 2009) in this range of redshifts. For the $\mathrm{He}\,{\rm{II}}$ photoionization rates we used ${C}_{\mathrm{He}{\rm{III}}}=1.5$ at all redshifts, which seems to allow us to accurately recover the input ionization models.

A.1.1. First-order Estimation of the Ionization History in Hydrodynamical Simulations from the Photoionization Rates

Using the relation between volume-averaged quantities derived above, one can obtain an analytical expectation of the evolution of the volume-averaged ${\rm{H}}\,{\rm{I}}$ ionization fraction, $\langle {x}_{{\rm{H}}{\rm{I}}}\rangle$, in optically thin simulations once a certain photoionization rate is assumed. This can be useful to predict the approximate ionization history that a certain photoionization rate model will produce in an optically thin hydrodynamical simulation. In particular, for the ${\rm{H}}\,{\rm{I}}$ photoionization rate we can rewrite Equation (22) above as

$$\begin{eqnarray}&&{C}_{{\rm{H}}{\rm{II}}}\langle {n}_{{\rm{H}}}\rangle (z){\alpha }_{{\rm{r}},{\rm{H}}{\rm{II}}}(\langle T\rangle )(1+\chi ){(1-{C}_{x,{\rm{H}}{\rm{I}}}\langle {x}_{{\rm{H}}{\rm{I}}}\rangle (z))}^{2}\\ &&\quad -{C}_{x,{\rm{H}}\,{\rm{II}}}^{2}\langle {{\rm{\Gamma }}}_{\gamma ,{\rm{H}}{\rm{I}}}\rangle (z)\langle {x}_{{\rm{H}}{\rm{I}}}\rangle (z)=0,\end{eqnarray} \tag{ 27 }$$

where now we are using that ${C}_{x,{\rm{H}}{\rm{II}}}\langle {x}_{{\rm{H}}{\rm{II}}}\rangle =1-{C}_{x,{\rm{H}}{\rm{I}}}\langle {x}_{{\rm{H}}{\rm{I}}}\rangle$. If we define $A={C}_{{\rm{H}}{\rm{II}}}\langle {n}_{{\rm{H}}}\rangle (z){\alpha }_{{\rm{r}},{\rm{H}}{\rm{II}}}(\langle T\rangle )(1+\chi )$, we can rewrite this equation as

$$\begin{eqnarray}&&A\langle {x}_{{\rm{H}}{\rm{I}}}{\rangle }^{2}(z)-[2A+\langle {{\rm{\Gamma }}}_{\gamma ,{\rm{H}}{\rm{I}}}\rangle (z)]\langle {x}_{{\rm{H}}{\rm{I}}}\rangle +A=0.\end{eqnarray} \tag{ 28 }$$

Solving for the volume-averaged ${\rm{H}}\,{\rm{I}}$ ionization fraction, we get

$$\begin{eqnarray}&&\langle {x}_{{\rm{H}}{\rm{I}}}\rangle =1-\langle {x}_{{\rm{H}}{\rm{II}}}\rangle =\displaystyle \frac{2A+{{\rm{\Gamma }}}_{\gamma ,{\rm{H}}{\rm{I}}}\pm \sqrt{{{\rm{\Gamma }}}_{\gamma ,{\rm{H}}\,{\rm{I}}}^{2}+4A{{\rm{\Gamma }}}_{\gamma ,{\rm{H}}{\rm{I}}}}}{2A},\end{eqnarray} \tag{ 29 }$$

which gives us an estimation of the expected hydrogen neutral fraction once we discard the nonphysical solution. It is expected that ${C}_{{\rm{H}}{\rm{I}}}$ changes between simulations that implement very different physical processes. Therefore, this formula could also be used to obtain the value of ${C}_{{\rm{H}}{\rm{I}}}$ in optically thin hydrodynamical simulations using one run and then compute analytically the expected results with different UVB models.

In Figure 11 we plot the result obtained from Equation (29) using the photoionization rate from the HM12 model (green line) and the late reionization model (blue line). We assumed $T=2\times {10}^{4}$ K at all redshifts, and we used the clumping given by Pawlik et al. (2009, their fit to C₁₀₀).

Zoom In Zoom Out Reset image size

A.2. Volume-averaged Heating Equations

The next step is to include in a consistent way the heating occurring during different reionization epochs. To do this, we have implemented a variation of an idea suggested by FG09. In this work the authors also raised the problems pointed out in the previous section of incorporating the effects of a prescribed UVB in cosmological hydrodynamical simulations under the assumption of an optically thin plasma. Their work was motivated by $\mathrm{He}\,{\rm{II}}$ reionization, but the idea can be easily applied for any reionization event. They proposed, as a more physically motivated approach, to increase the temperature of each gas element by an amount $d{\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}(z)/{dz}$ (which is subsequently allowed to cool) at each time step ${\rm{\Delta }}z$ in the simulation.²⁴ Here ${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}(z)$ is the total heat input or cumulative temperature increase evolution owing to reionization and is precomputed given the desired ${\rm{H}}\,{\rm{I}}$ reionization history. Both the total heat input from ${\rm{H}}\,{\rm{I}}$ reionization and that from $\mathrm{He}\,{\rm{II}}$ reionization will be free parameters in our model. We discuss the specific values used in this work in Section 4.

Now, we want to go one step further from the $d{\rm{\Delta }}T/{dt}$ idea and relate this change in temperature to a specific photoheating rate, $\dot{q}$. This will allow the use of our new models in standard hydrodynamical codes. To do this, we can, again, volume-average the relation between the heat per unit of time in one cell produced by the dT/dz model and the one produced by a photoheating rate. The change of internal energy density due to a certain change in temperature of the gas can be related to a new photoheating rate in the following way:

$$\begin{eqnarray}&&\displaystyle \frac{{{de}}_{\mathrm{int}}}{{dt}}=\displaystyle \frac{3{k}_{{\rm{B}}}}{2{m}_{{\rm{p}}}\mu }\displaystyle \frac{d{\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}}{{dt}}=\displaystyle \frac{{n}_{{\rm{H}}{\rm{I}}}{\dot{q}}_{{\rm{H}}{\rm{I}}}}{{\rho }_{\mathrm{proper}}}.\end{eqnarray} \tag{ 30 }$$

Using the same assumptions considered to obtain the new photoionization rates, we can get a volume-averaged value for the heating rate. We need also to assume that the volume-averaged molecular weight is $\langle \mu \rangle =(1+4\chi )\,/(1+\chi +\langle {n}_{{\rm{e}}}/{n}_{{\rm{H}}}\rangle )$. Then we obtain

$$\begin{eqnarray}&&{\dot{q}}_{{\rm{H}}{\rm{I}}}={C}_{\dot{q},{\rm{H}}{\rm{I}}}\displaystyle \frac{3{k}_{{\rm{B}}}}{2\langle \mu \rangle {X}_{{\rm{p}}}\langle {x}_{{\rm{H}}{\rm{I}}}\rangle (z)}\displaystyle \frac{d{\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}}{{dt}},\end{eqnarray} \tag{ 31 }$$

where ${C}_{\dot{q},{\rm{H}}{\rm{I}}}$ is a correction factor defined as

$$\begin{eqnarray}&&\langle \displaystyle \frac{\mu }{{x}_{{\rm{H}}{\rm{I}}}}\displaystyle \frac{d{\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}}{{dt}}\rangle ={C}_{\dot{q},{\rm{H}}{\rm{I}}}\displaystyle \frac{\langle \mu \rangle }{\langle {x}_{{\rm{H}}{\rm{I}}}\rangle }\langle \displaystyle \frac{d{\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}}{{dt}}\rangle .\end{eqnarray} \tag{ 32 }$$

Note that this approach sets ${\dot{q}}_{\mathrm{He}{\rm{I}}}=0.0$, as we include the heating produced by $\mathrm{He}\,{\rm{I}}$ reionization in the ${\rm{H}}\,{\rm{I}}$ heating rate. An identical approach is used to obtain the $\mathrm{He}\,{\rm{II}}$ heating rate, ${\dot{q}}_{\mathrm{He}{\rm{II}}}$. To compute these values in our models, we have used a correction factor, ${C}_{\dot{q}}=1$.

Therefore, in our model, once a total heat input due to reionization (${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}$) is chosen, the exact photoheating rates will depend on the assumption made on $d{\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}/{dz}$. In order to simplify this, we assume that the evolution of the total heat input can be well approximated by the volume-averaged ionization fraction evolution: $d{\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}/{dt}\sim d\langle {x}_{{\rm{H}}{\rm{II}}}\rangle /{dt}$. From here we can derive the derivative of the total heat input evolution with redshift, which is used in Equation (31) to obtain the tabulated photoheating rates that will go into the code. For this reason in our reionization models, the thermal history is defined based on one free parameter, the total heat input ${\rm{\Delta }}T$, and one free function, the ionization history, which in this context we define as the volume-averaged ionization fraction evolution, $\langle {x}_{{\rm{H}}{\rm{II}}}\rangle (z)$. In the case of ${\rm{H}}\,{\rm{I}}$ reionization the heating term is defined as

$$\begin{eqnarray}&&\left|\displaystyle \frac{d{\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}(z)}{{dz}}\right|=\displaystyle \frac{{\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}\left|\tfrac{d\langle {x}_{{\rm{H}}{\rm{II}}}\rangle (z)}{{dz}}\right|}{{\int }_{{z}_{\mathrm{reion}}}^{{z}_{\mathrm{recomb}}}\left|\tfrac{d\langle {x}_{{\rm{H}}{\rm{II}}}\rangle (z)}{{dz}}\right|{dz}}.\end{eqnarray} \tag{ 33 }$$

The denominator factor is a normalization to guarantee that the integrated amount of heating injected at each time step corresponds to the total heat input.

Here we present the ionization and thermal histories produced by other UVB models, widely used in the literature. We have run simulations using the following UVB models (HM96, HM01, FG09), in addition to the HM12 model largely discussed in this work. The left panel of Figure 12 shows the ionization history for all these models. It is particularly interesting to first focus on the results of FG09, which use the same approximation as HM12 (Equation (3)) to calculate the reionization redshift of their model that obtained complete $\mathrm{He}\,{\rm{II}}$ reionization by $z\sim 3$ and ${\rm{H}}\,{\rm{I}}$ reionization by z = 6. However, it is clear that this model is producing a much higher ${\rm{H}}\,{\rm{I}}$ reionization, more close to $z\sim 10$, indicating that they have the same problems as the HM12prescription. In fact, the method to compute their expected ionization history is the same as the one used by HM12. In this case the effect is not as extreme, only because tabulated values start at a lower redshift and, although they produce what seems to be a very high redshift ${\rm{H}}\,{\rm{I}}$ reionization, the final outcome of the model is still within CMB observational constraints (Planck Collaboration et al. 2016a).²⁵ We do not have any information on what were the expected ionization histories for the other two older Haardt & Madau models, but it is still very interesting to show their reionization histories, as they are still used in the literature. First of all, notice that in these models, as well as in the FG09 run, the ${\rm{H}}\,{\rm{I}}$ reionization is described as a simple step function. The HM01model assumes a much more shallower $\mathrm{He}\,{\rm{II}}$ reionization history than any other model considered in this work, starting as early as $z\sim 10$. This will allow us to see the effect of these types of models on the thermal history. HM96 ${\rm{H}}\,{\rm{I}}$ reionization happens very late, and therefore this model does not fulfill the Planck Collaboration et al. (2016a) CMB constraints.²⁶

Zoom In Zoom Out Reset image size

The middle and right panels of Figure 12 illustrate the thermal histories of these models. The first thing that we want to point out is the effect of a shallower $\mathrm{He}\,{\rm{II}}$ reionization model on the HM01 model. This basically produces a smoother evolution of the temperature at mean density, T₀, eliminating any sharp behavior at lower redshifts. At lower redshifts, the instantaneous thermal parameters γ and $T({\rm{\Delta }})$ are converging to the same values, but pressure smoothing scale and curvature are not.

Figure 13 shows results for hydrodynamical simulations in which a simple redshift cutoff was applied to the HM12 model. In this approach the ${\rm{H}}\,{\rm{I}}$ and $\mathrm{He}\,{\rm{I}}$ UVB rates are set to zero above a certain redshift: z = 11 (red), z = 9 (blue), and z = 7 (brown). We also show the original HM12 model (green lines) for comparison. All the runs share the same $\mathrm{He}\,{\rm{II}}$ rates. The left panel of Figure 13 shows the ionization histories of all these models. The middle and right panels present their thermal histories. As expected, by applying a cutoff to the UVB rates, the reionization redshift and its thermal signatures move down to the cutoff redshift. The gas pressure scale shows very clearly the effect of producing the heating due to ${\rm{H}}\,{\rm{I}}$ reionization at lower redshift. As discussed in Section 6, this is because the pressure smoothing scale depends on the full thermal history of the universe and not just on the instantaneous temperatures.

Zoom In Zoom Out Reset image size

In this appendix we present the results of simulations that only differ in the cosmological parameters and have the same UVB model, HM12. The random seeds in the initial condition are also the same. Table 2 summarizes the cosmological parameters used in each run. Cosmological model A is the default cosmological model used in this work. Models B and C were selected from the posterior distribution of the Planck results (Planck Collaboration et al. 2016) in order to differ as much as possible in the expected matter power spectrum. These should maximize the difference between the models while keeping them within the limits allowed by CMB observations. Model D was used in Lukić et al. (2015) and is also in agreement with the last CMB results.

Figure 14 shows the ionization and thermal histories for these four simulations. Both ionization histories, as well as the evolution of thermal parameters, are almost identical for all the runs. This means that within current observational constraints, the cosmological structure formation does not change significantly the thermal evolution of the IGM, which is determined by the UVB model used. In fact, this result is expected, as from our analytical calculation of the photoionization and photoheating rates we can easily calculate how these rates will change with some cosmological parameters. The difference in photoionization rate values for different cosmologies will be the difference between ${{\rm{\Omega }}}_{{\rm{b}}}{h}^{2}$ in the models. The difference in effective photoheating values for different cosmologies will be driven by the difference between H(z).²⁷

Zoom In Zoom Out Reset image size

We want to emphasize here that even when two simulations/models share the same ionization and thermal evolution, that does not mean that the Lyα observables (probability density function, flux power spectrum, etc.) from these runs will not show differences between them, as the observables could have their own dependence on cosmological parameter or other parameters.

We are making our new UVB models freely available for public use, so that they can be used by the whole community in future hydrodynamical simulations. We provide here the photoionization and photoheating rates for the late ${\rm{H}}\,{\rm{I}}$ reionization (Table 3), middle ${\rm{H}}\,{\rm{I}}$ reionization (Table 4), and early ${\rm{H}}\,{\rm{I}}$ reionization (Table 5) models, assuming a total heat input during reionization of ${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}=2\times {10}^{4}$ K for ${\rm{H}}\,{\rm{I}}$ and ${\rm{\Delta }}{T}_{\mathrm{He}{\rm{II}}}=1.5\times {10}^{4}$ K for $\mathrm{He}\,{\rm{II}}$. The rates for these models are also shown in Figure 15. We refer to Section 3 for a careful explanation on how these rates were obtained. In our new models we have also applied a small correction to the ${\rm{H}}\,{\rm{I}}$ and $\mathrm{He}\,{\rm{I}}$ photoionization rates of HM12 once reionization is completed to ensure that the simulations match current best observations of the ${\rm{H}}\,{\rm{I}}$ mean flux at different redshifts. The goal is to reduce the effect of the current standard post-process rescaling approach done with simulations that aim to reproduce Lyα statistics. The photoheating rates have been corrected by the same factor in order to keep the heat input per volume element constant in the models and therefore keep exactly the same thermal histories. Section 8 elaborates on the limitations and applicability of these models. Tables 3–5 are published in their entirety in a machine-readable format. Just a portion is shown here, as a guidance to their format and content. We encourage anyone interested in running some other specific model to contact the authors.

Zoom In Zoom Out Reset image size

Table 3.  Tabulated UV Background for Late ${\rm{H}}\,{\rm{I}}$ Reionization (${z}_{\mathrm{reion},{\rm{H}}{\rm{I}}}=6.55$), ${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}=2\times {10}^{4}$ K; He_A Reionization, ${\rm{\Delta }}{T}_{\mathrm{He}{\rm{II}}}=1.5\times {10}^{4}$ K

${\mathrm{log}}_{10}(z+1)$	${{\rm{\Gamma }}}_{{\rm{H}}{\rm{I}}}$	${{\rm{\Gamma }}}_{\mathrm{He}{\rm{I}}}$	${{\rm{\Gamma }}}_{\mathrm{He}{\rm{II}}}$	${\dot{q}}_{{\rm{H}}{\rm{I}}}$	${\dot{q}}_{\mathrm{He}{\rm{I}}}$	${\dot{q}}_{\mathrm{He}{\rm{II}}}$
	(s−1)	(s−1)	(s−1)	(erg s−1)	(erg s−1)	(erg s−1)
0.0000	5.700e−14	3.100e−14	1.122e−16	3.561e−25	4.486e−25	5.008e−27
0.0212	7.131e−14	3.942e−14	1.291e−16	4.466e−25	5.632e−25	5.729e−27
0.0414	8.817e−14	4.882e−14	1.564e−16	5.546e−25	6.944e−25	6.874e−27
0.0607	1.081e−13	6.037e−14	1.892e−16	6.806e−25	8.499e−25	8.215e−27

Note. Column (1): redshift (logarithm). Column (2): ${\rm{H}}\,{\rm{I}}$ photoionization rate. Column (3): $\mathrm{He}\,{\rm{I}}$ photoionization rate. Column (4): $\mathrm{He}\,{\rm{II}}$ photoionization rate. Column (5): ${\rm{H}}\,{\rm{I}}$ photoheating rate. Column (6): $\mathrm{He}\,{\rm{I}}$ photoheating rate. Column (7): $\mathrm{He}\,{\rm{II}}$ photoheating rate.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Table 4.  Tabulated UV Background for Middle ${\rm{H}}\,{\rm{I}}$ Reionization (${z}_{\mathrm{reion},{\rm{H}}{\rm{I}}}=8.30$), ${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}=2\times {10}^{4}$ K; He_A Reionization, ${\rm{\Delta }}{T}_{\mathrm{He}{\rm{II}}}=1.5\times {10}^{4}$ K

${\mathrm{log}}_{10}(z+1)$	${{\rm{\Gamma }}}_{{\rm{H}}{\rm{I}}}$	${{\rm{\Gamma }}}_{\mathrm{He}{\rm{I}}}$	${{\rm{\Gamma }}}_{\mathrm{He}{\rm{II}}}$	${\dot{q}}_{{\rm{H}}{\rm{I}}}$	${\dot{q}}_{\mathrm{He}{\rm{I}}}$	${\dot{q}}_{\mathrm{He}{\rm{II}}}$
	(s−1)	(s−1)	(s−1)	(erg s−1)	(erg s−1)	(erg s−1)
0.0000	5.700e−14	3.100e−14	1.122e−16	3.561e−25	4.486e−25	5.008e−27
0.0212	7.131e−14	3.942e−14	1.291e−16	4.466e−25	5.632e−25	5.729e−27
0.0414	8.817e−14	4.882e−14	1.564e−16	5.546e−25	6.944e−25	6.874e−27
0.0607	1.081e−13	6.037e−14	1.892e−16	6.806e−25	8.499e−25	8.215e−27

Note. Column (1): redshift (logarithm). Column (2): ${\rm{H}}\,{\rm{I}}$ photoionization rate. Column (3): $\mathrm{He}\,{\rm{I}}$ photoionization rate. Column (4): $\mathrm{He}\,{\rm{II}}$ photoionization rate. Column (5): ${\rm{H}}\,{\rm{I}}$ photoheating rate. Column (6): $\mathrm{He}\,{\rm{I}}$ photoheating rate. Column (7): $\mathrm{He}\,{\rm{II}}$ photoheating rate.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Table 5.  Tabulated UV Background for Early ${\rm{H}}\,{\rm{I}}$ reionization (${z}_{\mathrm{reion},{\rm{H}}{\rm{I}}}=9.70$), ${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}=2\times {10}^{4}$ K; He_A reionization, ${\rm{\Delta }}{T}_{\mathrm{He}{\rm{II}}}=1.5\times {10}^{4}$ K

${\mathrm{log}}_{10}(z+1)$	${{\rm{\Gamma }}}_{{\rm{H}}{\rm{I}}}$	${{\rm{\Gamma }}}_{\mathrm{He}{\rm{I}}}$	${{\rm{\Gamma }}}_{\mathrm{He}{\rm{II}}}$	${\dot{q}}_{{\rm{H}}{\rm{I}}}$	${\dot{q}}_{\mathrm{He}{\rm{I}}}$	${\dot{q}}_{\mathrm{He}{\rm{II}}}$
	(s−1)	(s−1)	(s−1)	(erg s−1)	(erg s−1)	(erg s−1)
0.0000	5.700e−14	3.100e−14	1.122e−16	3.561e−25	4.486e−25	5.008e−27
0.0212	7.131e−14	3.942e−14	1.291e−16	4.466e−25	5.632e−25	5.729e−27
0.0414	8.817e−14	4.882e−14	1.564e−16	5.546e−25	6.944e−25	6.874e−27
0.0607	1.081e−13	6.037e−14	1.892e−16	6.806e−25	8.499e−25	8.215e−27

Note. Column (1): redshift (logarithm). Column (2): ${\rm{H}}\,{\rm{I}}$ photoionization rate. Column (3): $\mathrm{He}\,{\rm{I}}$ photoionization rate. Column (4): $\mathrm{He}\,{\rm{II}}$ photoionization rate. Column (5): ${\rm{H}}\,{\rm{I}}$ photoheating rate. Column (6): $\mathrm{He}\,{\rm{I}}$ photoheating rate. Column (7): $\mathrm{He}\,{\rm{II}}$ photoheating rate.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image