The long-lasting effect of X-ray preheating in the post-reionization intergalactic medium

The Universe is currently believed to transition from primarily neutral to mainly ionized around $z\sim 7.7$ (Planck Collaboration et al., 2018). This phase transition, known as cosmic reionization, occurs via the passage of ionization fronts through the intergalactic medium (IGM) resulting in a violently perturbed IGM that gets heated up to $\sim 10^{4}$ K (Miralda-Escudé & Rees, 1994; Furlanetto & Oh, 2009; Zeng & Hirata, 2021). The heating of the IGM is primarily driven by ultraviolet (UV) photons that originate from the first stars and galaxies, and thus reionization happens in a patchy fashion where denser regions reionize first due to a larger abundance of UV sources (see e.g. McQuinn, 2016, for a review). Although the physics of reionization are – to a certain extent – understood, finer details like the timeline of the reionization process (e.g. Morales et al., 2021; Bosman et al., 2022; Lu et al., 2022; Jin et al., 2022) and the nature and evolution of the ionizing sources (Lidz et al., 2021; Kostyuk et al., 2022; Garaldi et al., 2022) are currently the focus of ongoing efforts.

Furthermore, the violent injection of energy into the IGM during reionization coupled with the inhomogeneous nature of reionization leads to relics in the post-reionization IGM that can survive to low redshifts ($z<4$) (Hirata, 2018). These imprints of reionization, which affect the evolution of the post-reionization IGM (Trac et al., 2008; Wu et al., 2019; Oñorbe et al., 2019), can be categorized into two types – the “void response” and the “halo response”. The void response is responsible for the imprints of reionization that affect the Ly$\alpha$ forest at low redshifts ($z\lesssim 4$) (Montero-Camacho et al. 2019, hereafter referred to as M19). The reason the (mini) voids play a crucial role is that they not only get ultraviolet heated by the ionization sources but also experience shock heating, and accompanying compression to mean density, from shocks that originate in the surrounding denser regions. Overall, this results in a long-lasting impact since they can get pushed to a higher adiabat than the denser gas and can reach temperatures of $\sim 3\times 10^{4}$ K (Hirata, 2018). On the other hand, we define the halo response to encompass the effects that affect the denser regions. Gas in these regions can likely recover from the reionization “energetic kick” by $z\sim 3.5-4$ (Long et al., 2022a). The halo response, i.e. the imprints left by the passage of ionization fronts on shallow potential wells, not only destroy small-scale baryonic structure ($\sim 10^{6}$ M_⊙) (Shapiro et al., 2004) but also suppress the growth of structure for cosmological time scales, eventually reaching a filtering scale of $\sim 100$ kpc after hundreds of millions of years (D’Aloisio et al., 2020). After cosmic reionization, the thermal evolution of the IGM was believed to follow a stringent power law, the so-called temperature-density relation (Hui & Gnedin, 1997), which arises as a consequence of the scaling of the recombination rate with temperature and it is facilitated by both adiabatic expansion and Compton cooling (McQuinn & Upton Sanderbeck, 2016). Although both types of reionization relics affect the relaxation/recovery of the temperature-density relation, only the void response leads to signatures in the Ly$\alpha$ forest at low redshifts ($z\leq 4)$. On the other hand, the evolution of the filtering of masses, i.e. how many baryons can sit on a given small halo host, is greatly impacted by the halo response to the passage of ionization fronts (Long et al., 2022b) and it is vital to compute the impact of inhomogeneous reionization in the post-reionization H i 21 cm intensity mapping. Note that both void and halo responses are somewhat artificial separations made here for an easier understanding. Besides, both effects need to be considered while simultaneously accounting for the inhomogeneous nature of reionization for an accurate model of the post-reionization IGM.

These imprints from hydrogen reionization leave strong ripples that cosmological probes of the post-reionization IGM can recover. For instance, for the Ly$\alpha$ forest, i.e. the absorption features due to neutral hydrogen absorption of Lyman-$\alpha$ radiation emitted from distant quasars, its power spectra have been shown to suffer an enhancement at large scales due to the fossils from H i reionization (M19; Molaro et al. 2022; Puchwein et al. 2022)¹¹1However, Mishra & Gnedin (2022) found no enhancement at large scales; nevertheless, relics from cosmic reionization in the post-reionization era were still present. The reason for this disagreement is still under investigation.. In addition, the line intensity mapping (LIM) of neutral hydrogen 21 cm hyperfine transition also suffers a similar effect, which is driven by the second type of relic, i.e. the hydrodynamic response (Long et al., 2022a).

Naturally, this novel broadband signal, which originates during the epoch of reionization, depends on the astrophysics that governs cosmic reionization (Montero-Camacho & Mao 2020, hereafter referred to as MM20). The ability to learn about the timeline of reionization (Montero-Camacho & Mao, 2021) or to separate cosmology from reionization astrophysics (Montero-Camacho et al., 2023) will depend on a clear holistic understanding of the reionization process and on accurate handling of other astrophysical and instrumental systematics that will be present in both 21 cm LIM and Ly$\alpha$ forest. Here in this paper, we will focus on one piece of the puzzle, namely the impact of X-ray preheating in the post-reionization IGM.

During a period preceding cosmic reionization, X-rays heat the IGM up and raise the gas temperature from below to above the cosmic microwave background (CMB) temperature, leading to a transition of observing this epoch of heating from absorption to emission through the 21 cm line relative to the CMB. X-rays have small cross-sections for interactions with atoms. Hence, they penetrate deep into the neutral IGM and subsequently preheat the IGM through photoelectric absorption. However, the impact of X-ray preheating in the high-redshift IGM ($z\sim 15$) is highly uncertain due to a lack of observational constraints. Thus, theoretical models of the IGM often consider medium IGM temperatures ranging from tens to hundreds of Kelvins (see Figure 4 of Ross et al., 2017). Furthermore, other sources of uncertainty are the extrapolation of low-redshift X-ray luminosity functions to high redshifts (Lehmer et al., 2022), uncertainty regarding the abundance, clustering, and nature of the heating sources, including X-ray binaries (see, e.g. Xu et al., 2014; Madau & Fragos, 2017; Kaur et al., 2022; Sartorio et al., 2023), high-$z$ quasars and active galactic nuclei (Venkatesan & Benson, 2011; Madau & Haardt, 2015; Ross et al., 2019)²²2Nevertheless, Eide et al. (2018) claimed that, based on their seeding prescription, nuclear black holes are too scarce and thus do not contribute significantly to heating or ionization at $z>10$., cosmic rays from supernovae in Population III stars (Gessey-Jones et al., 2023), and hot interstellar medium emission (Pacucci et al., 2014; Ma et al., 2021).

As a result of X-ray preheating, some structures at small scales “puff up” before the passage of an ionization front. Therefore, the relaxation into the usual temperature-density relation of the IGM should occur faster (Hirata, 2018), which would lead to an overall suppression of the large-scale enhancement in both the Ly$\alpha$ forest and on the 21 cm LIM imprint in the post-reionization era. In addition, the presence of an anisotropic radiation field due to the sources of X-ray photons impacts the high-redshift IGM in several ways. For instance, X-rays should first increase the 21 cm power spectrum on large scales (Pritchard & Furlanetto, 2007; Ewall-Wice et al., 2016), followed by an overall reduction of the strength of 21 cm fluctuations due to X-rays wiping off temperature fluctuations in small scales. Moreover, the presence of such energetic sources is expected to leave clear imprints on the 21 cm bispectrum (Watkinson et al., 2019; Ma et al., 2021) and to significantly increase non-Gaussianity (Ross et al., 2017; Ross et al., 2019). Furthermore, X-ray preheating yields a sizable contribution to the early stages of reionization, which can lead to an earlier, albeit slightly extended, epoch of reionization and a more uniform H ii morphology (Mesinger et al., 2013).

Recently, HERA Collaboration et al. (2022) has placed constraints on the X-ray luminosity and on the effect of X-ray preheating in the cosmic dawn, disfavoring the extrapolation of the local relationship between soft X-ray luminosity and star formation (Mineo et al., 2012) to high redshifts. On the other hand, their results are consistent with X-rays being produced by a population of low-metallicity high-mass X-ray binaries (Fragos et al., 2013). Moreover, HERA Collaboration et al. (2022) finds that the intergalactic medium must have been heated above the CMB temperature of $\sim 31\,{\rm K}$ as early as $z=10.4$.

In principle, preheating of the IGM could leave observable signatures in the post-reionization era through two main effects. First, X-rays will wipe out some of the small-scale structures that would otherwise be present. Second, X-ray preheating sets up the initial conditions for cosmic reionization since X-ray photons heat up the neutral gas before the UV ionizing sources, thus leading to a contribution to the global ionized hydrogen fraction prior to the start of proper cosmic reionization ($z\sim 12$). Regarding the latter effect, MM20 demonstrated that the acceleration or deceleration of the initial stages of reionization, i.e. increase or decrease of $\overline{x}_{\rm HII}(z\sim 12)$, has a minimal impact on the long-lasting relics from reionization in the IGM (as long as non-extreme scenarios of X-ray preheating are considered). In this work, we focus on the former effect, i.e. we investigate the impact of the “puffing-up” of structures during the early stages of the reionization process. We highlight that assessing the dependence of the impact of inhmogeneous reionization on cosmic dawn astrophysics is a vital step toward future Bayesian frameworks for the inference of unbiased cosmological parameters from cosmological probes of the post-reionization era (Montero-Camacho et al., 2023).

The rest of this paper is organized as follows. We introduce the different models for X-ray preheating of the IGM considered throughout this work in §2 and describe our simulations in §3. In §4, we establish the impact of preheating on the IGM. In particular, we quantify the expected impact on the 3D and 1D Ly$\alpha$ power spectrum in §4.1 and §4.2, respectively. Note that in order to facilitate the exploration of several different X-ray preheating models, the calculations in §4 sacrifice accuracy due to the use of small box lengths. Thus, we make a separate analysis with larger boxes but with only a single X-ray preheating model in §5. We summarize our findings in §6. We show the impact of preheating on the IGM temperature-density relation in Appendix A. Some intermediate, tabulated results are left to Appendix B.

Throughout this work, we assume the following fiducial cosmology: $h=0.6774$, $\Omega_{\rm c}h^{2}=0.1417$, $\Omega_{\rm b}h^{2}=0.02230$, $A_{\rm s}=2.142\times 10^{-9}$, $n_{\rm s}=0.9667$. Our chosen cosmology is close to Planck 2015 “TT + TE + EE + lowP + lensing + ext” result (Planck Collaboration et al., 2016). Our fiducial model has no X-ray preheating.

During the epoch of X-ray heating, X-ray photons are expected to raise the temperature of the IGM to higher levels than the CMB temperature, thereby leading to an observational signature in the 21 cm signal. Although the precise degree of heating and its redshift evolution is uncertain, it is nonetheless vital to understand X-ray preheating since it is degenerate with reionization astrophysics (Ewall-Wice et al., 2016). However, the reionization relics in the post-reionization IGM were found to be insensitive to X-ray preheating (MM20), thus allowing for a possible window for constraining the astrophysics of reionization without worrying about these degeneracies (Montero-Camacho & Mao, 2021). Here we revisit that claim, using high-resolution simulations capable of tracking the way the small-scale structure reionizes, and using realistic X-ray preheating models.

Throughout this paper we are interested in the observational signatures that cosmic dawn astrophysics (in this case X-ray preheating) leaves on the 3D Ly$\alpha$ forest power spectrum via the impact of X-rays on the imprints from H i reionization that survive in the IGM to low redshifts ($z<4$). Before covering our implementation of X-ray preheating, we briefly describe the fundamental ideas behind the reionization relics that make it possible to have any X-ray preheating signatures.

Inhomogeneous reionization seeds temperature fluctuations, fluctuations in the pressure smoothing, and affects the ionizing radiation field (see, e.g. Puchwein et al., 2022). The resulting smoking gun is the large-scale enhancement present in the Ly$\alpha$ forest power spectra (M19; Oñorbe et al. 2019; Wu et al. 2019; Molaro et al. 2022)³³3Molaro et al. (2023) presented evidence for a possible measurement at $z\geq 4$.. However, accurately modeling the imprint of reionization relics on the low-redshift post-reionization IGM ($z\leq 4$) requires great mass resolution because the gas in underdense regions proves crucial. Gas in minivoids gets heated by the UV ionization front to $\sim 2\times 10^{4}$ K. In addition, this underdense gas can be shock-heated by shocks that originate in denser regions. As a result, the underdense gas is compressed to the median density and gets heated to even higher temperatures than gas in denser regions (Hirata, 2018). This additional heating, which originates on the way the small-scale structure reionizes and takes cosmological timescales to relax, is the reason for the memory of reionization in the low-redshift ($z\sim 2$) Ly$\alpha$ forest.

To compute the Ly$\alpha$ forest 3D flux power spectrum, including the imprints from H i inhomogeneous reionization, we will use the model described in §4.1, which consists of using the linear matter power spectrum, the biasing model and non-linear correction from Arinyo-i-Prats et al. (2015), and a simulation-based hybrid method to obtain the reionization relics. The reason for the hybrid methodology is that the dynamical range is too large for a single simulation because we need not only the statistical power to capture the patchy nature of reionization but also great mass resolution to accurately track the way the small-scale structure reionizes.

Hirata (2018) established the impact of X-ray preheating on the BAO signal of the Ly$\alpha$ forest using a slightly unrealistic X-ray preheating model, where the transition from absorption to emission of the 21 cm brightness temperature occurs at $z<16$, while the EDGES collaboration could have discovered the absorption peak of the signal centered around $z=17$ (Bowman et al., 2018)⁴⁴4The cosmological nature of this result remains controversial (see, e.g., Singh et al. 2022).. For reference in our analysis, we include the X-ray model introduced in Hirata (2018), hereafter referred to as “Hirata-300 model”, in which the additional injected energy into the IGM due to X-rays is given by

In order to investigate the impact of X-ray preheating in the imprints from patchy reionization, we utilize the X-ray prescription used in 21cmFASTv1.3 (Mesinger & Furlanetto, 2007; Mesinger et al., 2011), where the heating rate per particle is computed by summing the contributions from sources located in concentric spherical shells. We refer the interested reader to Mesinger et al. (2011) for the explicit prescription. Here we only recap the steps that highlight the free parameters in 21cmFAST.

The total X-ray emission rate is a necessary ingredient to compute the X-ray heating per baryon in the pre-reionization era. This rate per redshift interval from luminous sources between $z^{\prime\prime}$ and $z^{\prime\prime}+dz^{\prime\prime}$ is

in units of ${\rm s}^{-1}$, where $\zeta_{X}$ is the number of X-ray photons per solar mass, which manage to escape their parent galaxies, and the remaining terms in the right-hand-side of Eq. (2) altogether correspond to the total star formation rate inside a spherical shell of inner (outer) radius at the comoving redshift $z^{\prime\prime}$ ($z^{\prime\prime}+dz^{\prime\prime}$). $f_{*}$ is the fraction of baryons that becomes stars, $f_{\rm coll}$ is the collapsed fraction, $\Omega_{b}$ is the baryon density, $\rho_{\rm crit,0}$ is the critical density (today), and ${\rm R^{\prime\prime}}$ is the comoving distance between $z^{\prime}$, the redshift of interest, and $z^{\prime\prime}$. $\delta_{\rm nl}^{\rm R^{\prime\prime}}$ is the evolved – with respect to redshift $z^{\prime\prime}$ – density smoothed on scales of ${\rm R^{\prime\prime}}$. The X-ray efficiency, $\zeta_{X}$, is one of the free parameters of the X-ray prescription of 21cmFASTv1.3, and its role is to control the degree of X-ray heating prior to reionization. Higher values would eventually cause reionization to occur earlier.

To compute the arrival rate, i.e. the number of X-ray photons per unit time per unit frequency bin with frequency $\nu$ seen at $(\boldsymbol{x},z^{\prime})$, from sources within $z^{\prime\prime}$ and $z^{\prime\prime}+dz^{\prime\prime}$, 21cmFAST assumes that the X-ray luminosity of sources follows a power law of the form, $L\propto(E/E_{0})^{-\alpha}$, with $\alpha$ being a constant and $E_{0}$ being the energy threshold for the lowest energy X-ray photons not absorbed by galaxies, and hence escaping into the IGM. The arrival rate can then be written as a function of the total X-ray emission rate (defined in Equation 2) as follows

where $E_{0}=h\nu_{0}$ and the exponential term accounts for the IGM attenuation for X-ray photons. Finally, the code calculates the X-ray heating per baryon by integrating Eq. (3) over frequency and $z^{\prime\prime}$ (see the detailed derivation in §3.1.2 of Mesinger et al. 2011).

The energy threshold, $E_{0}$, is the other free parameter in this X-ray prescription that is allowed to vary in this work. A larger value corresponds to inefficient X-ray heating of the IGM due to more X-ray photons being absorbed by the host galaxies, i.e. fewer photons preheat the IGM, and thus there is a slight delay in the reionization process.

We use the “heating models” of MM20 to establish the range of IGM temperatures that we consider throughout this work. In summary, the “Fast-Fid” model has $E_{0}=500$ eV and $\zeta_{X}=2\times 10^{56}$ M${}_{\odot}^{-1}$. This should not be confused with the overall Fiducial (“Fid”) model that has no X-ray preheating whatsoever. Furthermore, three variations on the energy threshold are considered $E_{0}=\{100,1000,1500\}$ eV and three additional variations for the X-ray ionization efficiency $\zeta_{X}=\{1,4,8\}\times 10^{56}$ M${}_{\odot}^{-1}$. For the nomenclature, we refer to the “E2” model as the one with $E_{0}=1000$ eV, the “$\zeta_{X}3$” model has $\zeta_{X}=8\times 10^{56}$ M${}_{\odot}^{-1}$, and similarly for other models. We summarize the nomenclature of models in Table 1.

In Figure 1, we plot the gas temperature evolution with respect to the redshift with different prescriptions of X-ray preheating considered herein, based on the simulations presented in §3. Note that the $\zeta_{X}3$ model has the earliest transition from 21 cm absorption to emission ($z\approx 15$). The bulk of models experiences the heating transition from $10.7\lessapprox z\lessapprox 13.0$, while the latest model to surpass the CMB temperature is the E3 model where it happens at $z=9$. This feature is expected for the E3 model since most X-ray photons are absorbed by their host galaxies in this photon-starved scenario, therefore leading to a delayed epoch of heating.

From Figure 1, we observe that, at $z=10$, the gas temperature ranges from $19$ to $564$ K. We highlight that the high-temperature end is consistent with the findings of Ross et al. (2017) and the models with the early heating transition of Fialkov et al. (2014), while the low end is roughly consistent, albeit slightly earlier for several models, with the current constraints from HERA Collaboration et al. (2022). Interestingly, the E1 model showcases some slightly deviant behavior with respect to the other models during the epoch of reionization. We attribute this feature to the photon-abundant nature of the X-ray model; however, due to the lower X-ray efficiency compared to that of model $\zeta_{X}3$ it takes a longer time for those X-ray photons to make their impact. Besides, they get a boost on their escape efforts thanks to the UV ionization going on during reionization.

The 21cmFAST simulation suite of MM20 computes the gas temperature up to $z\approx 35$. Hence, to ensure that X-rays do not permeate the whole dark age (Pritchard & Loeb, 2012), the high-redshift evolution of the gas temperature contribution due to X-rays is assumed to follow a power law (Hirata, 2018) as follows. For $1+z>40$,

Note that Eq. (4) does not have the same redshift evolution as the Hirata-300 model and, because of the low values of the added gas temperature, the impact on the hydrodynamics outside of the epoch of reionization and epoch of heating will be negligible. For instance, $\Delta T^{\rm X}_{\rm gas}$ is one (two) order of magnitude lower than $T_{\rm CMB}$ at $z=50$ ($z=100$). The reason for choosing a different redshift evolution than the Hirata-300 model is because Figure 1 shows that the 21cmFAST models do not decay strongly during the cosmic dawn, hence we chose a slightly less sudden drop. Given the large difference between the injected temperature of the gas (due to X-rays) and the CMB temperature at those high redshifts, the choice of knee is unlikely to play any role in our analysis — as long as there is one to limit the role of X-rays in the dark age.

We follow M19; MM20; Long et al. 2022a on simulating the impact of inhomogeneous reionization in the post-reionization IGM via a hybrid strategy that uses small box high-resolution simulations, which can accurately track down how gas reionizes below the Jeans mass prior to the passage of an ionization front, and large box semi-numerical simulations, which account for the patchy nature of reionization. We use a modified version of Gadget-2 (Springel et al., 2001, 2005; Hirata, 2018) for the small boxes and 21cmFAST (Mesinger & Furlanetto, 2007; Mesinger et al., 2011) for the large boxes.

The large box simulation suite was taken from MM20 and consists of 21cmFASTv1.3 simulations with a box size of 400 Mpc on each side, with $768^{3}$ cells for the matter field, and $256^{3}$ cells for the H i field. For each heating model considered in this work, we run 21cmFAST to generate the IGM temperature, neutral hydrogen field, and density contrast evolution during the epoch of reionization and cosmic dawn. The neutral hydrogen field and density contrast are needed to incorporate the inhomogeneous nature of reionization when calculating the imprint of reionization on the Ly$\alpha$ forest (see Equation 11 below) while the gas temperature is needed as an input for the Gadget-2 simulations.

Our Gadget-2 simulations include adiabatic expansion (including Hubble expansion), shock heating, Compton heating and cooling for neutral gas (accounting for residual ionization). For ionized gas, we also include Compton cooling, He ii cooling, recombination cooling, photoionization heating, and free-free cooling. For details, see §4.3 of Hirata (2018).

We use two different configurations for the small box simulations — “Type I”, a smaller configuration used to swiftly explore the qualitative impact of different X-ray models (see §4 below), and “Type II”, a larger configuration for comparison purposes using a couple of selected models (see §5 below). Both configurations use high-mass resolution simulations which were extensively tested and described in Hirata (2018). The Type I simulations have a box size of 425 kpc on each side, $2\times 128^{3}$ total number of particles, streaming velocities between baryons and dark matter of 33 km/s (r.m.s. at decoupling), a dark matter particle mass of $1.21\times 10^{3}$ M_⊙, and a gas-particle mass of $2.27\times 10^{2}$ M_⊙. The Type II simulations have a box size of 1275 kpc on each side, $2\times 192^{3}$ total number of particles, streaming velocity of 33 km/s , a dark matter particle mass of $9.72\times 10^{3}$ M_⊙, and a gas-particle mass of $1.81\times 10^{3}$ M_⊙. Note that the Type II simulations have been included in this analysis since the Type I boxes are not large enough to reliably sample the linear regime at the redshift of interest for a Ly$\alpha$ forest survey. In contrast, the Type II boxes are large enough to form multiple Ly$\alpha$ clouds and thus provide enhanced statistical power. Naturally, the ideal scenario would have been to have a single simulation cover both the small-scale structure response to the reionization process and the inhomogeneous nature of reionization. Nevertheless, the dynamical range would be too large to be computationally affordable.

We ran four different realizations of each Type I and Type II model considered throughout this work. We neglect the variance on the 21cmFAST simulations since previous work found that the error budget is dominated by the Gadget simulations (M19). Although we have computed the sample variance on the mean, for most of our Type I results we do not include this in our plots since the plots are too clustered due to the many models being considered. In contrast, we illustrate the Monte Carlo scatter due to the different realizations of the Fiducial and Fast-Fid model in §5.

For a description of the heating and cooling processes that are implemented in the small box simulations, we refer readers to §4 of Hirata (2018). Here we only recap the inclusion of X-ray preheating.

In Gadget-2, all heating and cooling rates are written as

where $k_{\rm B}$ is the Boltzmann’s constant, $n$ is the total number density of particles, $T$ is the gas temperature, $\dot{E}_{i}$ is the net volumetric heating rate in erg cm^-3 s^-1 from process $i$, $A=p\rho^{-5/3}$ is the particle entropy in the code, and $\dot{A}_{i}/A$ is the fractional rate of increase of thermal energy from process $i$.

To include X-ray preheating into the Gadget-2 simulations, it is sufficient to use Eq. (5) to account for the extra heating $\Delta T$ due to X-rays, i.e. $\dot{A}_{X}/A\propto\dot{\Delta T_{X}}/T$. Thus, one needs to replace $\Delta T_{X}$ using Eq. (1) in the Hirata-300 model or use results generated with 21cmFAST in other models. Note that Eq. (4) is used here for every 21cmFAST model to guarantee that X-rays do not permeate the dark age, i.e. predominant X-ray sources at very high redshifts that could compete with the CMB should not exist.

In Figure 2, we plot the gas distribution for our overall Fiducial model, the Hirata-300 model, and the $\zeta_{X}3$ model at different redshifts, prior to the start of the local reionization process. It is straightforward even by eye to grasp the impact of X-ray preheating — without X-ray preheating, the Fiducial model exhibits a much richer structure at all the shown redshifts when compared to the latter two X-ray models. Note that the skeleton of the denser regions remains although there is significant smearing in both the Hirata-300 and the $\zeta_{X}3$ model. In contrast, the underdense areas have been severely diminished. Regarding the differences between the two X-ray prescriptions, while they can be found by eye, particularly at the lower redshifts, we can easily quantify the differences in terms of the transparency of the intergalactic medium (see §4 below).

As described in §4.4 of Hirata (2018), to compute the Ly$\alpha$ transmission from our small-box snapshots we first select one of the axes as the line-of-sight direction (we later average over the results for all three axes). We assign the H i abundance proportional to $\Delta^{2}\alpha_{A}(T)$ where $\alpha_{A}$ is the Case A recombination rate and $\Delta=1+\delta_{\rm b}$ is the gas density in units of the mean baryon density in the Universe. The neutral hydrogen particles are then interpolated into a grid based on their redshift-space positions. A Gaussian with thermal width corresponding to the temperature $T$ and the mass of the hydrogen atom is used for the line-of-sight interpolation and we allow the Gaussian to wrap up by including periodic-box images due to the small box size. This procedure results in an opacity cube with arbitrary normalization. We obtain the normalization $\tau_{1}$ by matching our results to the observed mean flux from Kim et al. (2007).

Thus, the relevant observable from the Gadget-2 simulations is the $\tau_{1}$ value as a function of the redshift of observation and local redshift of reionization, and what we aim to identify is how $\tau_{1}$ changes with the different heating models described in Table 1.

MM20 demonstrated that while X-ray preheating can affect the initial stages of the reionization process, e.g. by raising the global $\overline{x}_{\rm HII}$ to $\sim 0.1$ before UV heating takes over, this effect has a minimal impact on the imprints from inhomogeneous reionization in the post-reionization IGM, and therefore the resulting memory of reionization in the Ly$\alpha$ forest, i.e. the large-scale enhancement in the transmitted flux power spectrum, would be capable of probing the astrophysics of reionization without dealing with the degeneracy introduced by cosmic dawn astrophysics.

However, X-rays can penetrate deeply into the IGM, which can also effectively wipe out some of the small-scale structures prior to the beginning of the reionization process. As seen in Figure 2, this impact appears to be even larger in the underdense regions. Underdense regions are crucial to establish the strength of the memory of reionization in the Ly$\alpha$ forest because the temperature-density relation becomes bimodal due to the underdense regions suffering shock-heating and compression to mean density in addition to UV heating (Hirata, 2018). The impact of X-ray preheating on the temperature-density relation of the IGM was left unexplored in MM20. Armed with the simulations herein, we can now investigate in this paper the long-lasting impact of X-ray preheating in the post-reionization IGM via the impact of inhomogeneous reionization (M19), through its effect on the temperature-density relation of the IGM. See Appendix A for the effect of X-ray preheating in the temperature-density relation.

The fluctuations in Ly$\alpha$ flux transmission, including the effect of reionization, are given by (M19)

where $b_{\rm F}$ is the flux bias, which is negative since the presence of matter implies absorption, $\beta_{\rm F}$ is the redshift-space distortion parameter, $\mu$ is the cosine of the angle of the $\boldsymbol{k}$ vector with respect to the line of sight, $\delta_{m}$ is the matter density contrast, $b_{\Gamma}$ is the radiation bias (Arinyo-i-Prats et al., 2015; Hirata, 2018) needed here to convert from optical depth fluctuation to flux. The transparency of the IGM, $\psi$, is a function of both redshift of observation $z_{\rm obs}$ and the redshift of local reionization $z_{\rm re}$, and is defined by the relative strength to a fiducial scenario where reionization happens at $z_{\rm re}=8$ in the Fiducial model (no X-ray preheating), i.e.

Often in Ly$\alpha$ simulations, one varies the normalization — in this case, $\tau_{1}$ — which is equivalent to varying the ionizing background, until one obtains the correct mean transmitted flux in our simulations. Specifically, $\tau_{1}$ is the optical depth that must be assigned to a patch of gas with mean density and temperature $T=10^{4}$ K in order for the mean transmitted flux to match observations. $\psi$ parametrizes the variations of transparency for an opacity cube that suddenly reionizes at $z_{\rm re}$, and is observed at $z_{\rm obs}$, relative to an opacity cube in the reference simulation with no X-ray preheating and that reionizes at $z_{\rm re}=8$. Note that $\psi<0$ implies that the reference simulation is more transparent or equivalently has a higher $\bar{F}$.

We plot the relative transparency of the IGM for our Type I high-resolution small-box simulations in Figure 3. Irrespective of the redshift of observation, some models like the Fast-Fid and Hirata-300 exhibit a small preference for $\psi^{\rm X-ray}<\psi^{\rm Fid}$ ⁵⁵5This is easier to appreciate in Table 2. for early reionization scenarios ($z_{\rm re}>8$). Therefore, the consequence of including any X-ray prescription is a less transparent scenario relative to that of the Fiducial model at any $z_{\rm obs}$ for these scenarios. However, as detailed in §5, this trend could be artificial due to the Type I boxes not having a large enough size for enough Ly$\alpha$ clouds inside them leading to not-so-accurate sampling. Reassuringly, the different X-ray models, including the simple Hirata-300 model and the Fiducial model, converge for very early reionization ($z=12$). This can be attributed to two factors. First, the very early start of the reionization process results in less time for X-rays to permeate the IGM. Second, the strong UV radiation field dominates over the X-ray radiation in this case, further limiting the impact of different X-ray prescriptions in this regime. Besides, the gas does not react instantly to X-ray preheating, i.e. the puffing up of structures takes time, and $z_{\rm re}=12$ has the shortest response time.

As seen in Figure 3, the same general trend from the early reionization scenarios, $\psi^{\rm X-ray}<\psi^{\rm Fid}$, is present for late reionization scenarios. In this case, the trend is more pronounced and is still generally insensitive to the redshift of observation. The consequence of adding X-rays to our Type I simulations is that the IGM becomes less transparent compared to the Fiducial scenario at a given local redshift of reionization. However, these trends are partially hindered by the fact that the Type I simulations do not accurately sample the forest (see §5). In Figure 3, the magenta-dashed line indicates where the equality with the reference scenario occurs. The first X-ray model to approach equality is the E3 model at $z_{\rm re}\approx 7.99$, while the last model is the $\zeta_{X}3$ model that crosses the threshold at $z_{\rm re}\approx 7.5$. These values vary slightly depending on the redshift of observation; however, the order of the models does not change.

For clarity, we tabulate the data used in Figure 3 in Table 2, including redshifts that were omitted in the figure. From Table 2, the largest discrepancy between the X-ray prescriptions and the Fiducial model occurs at $z_{\rm re}=6$ and $z_{\rm obs}=2.5$, where the Fiducial model predicts that the IGM is 7.840% more transparent than the reference scenario (i.e. Fiducial model with $z_{\rm re}=8$), compared to 3.489% more transparent in the $\zeta_{X}3$ model than the reference scenario, i.e. a reduction by more than a factor of two. On the other hand, the smallest deviations tend to occur in early reionization scenarios due to the less prominent role of X-ray preheating in these settings.

An interesting feature of Figure 3 is the larger transparency of the fiducial model compared to that of the X-ray models. In principle, one would have expected that the increase in temperature should smooth things out initially since the Ly$\alpha$ optical depth decreases with increasing temperature (McDonald & Miralda-Escudé, 2001; Bolton et al., 2005). We note that the Type-II boxes do recover this expected behavior as seen in §5.

We now move on to quantify the effect of X-ray preheating on the observables that we are considering in this work — the 3D and 1D flux power spectrum.

In this section, we showcase the impact of X-ray preheating in our Type II boxes, which have a box size of 1275 kpc (a factor of 3 larger than the Type I boxes). These Type II boxes are large enough to form multiple Ly$\alpha$ clouds and thus enhance our statistical power. Thus, the Type II simulations more accurately sample the linear regime for the redshifts of interest for Ly$\alpha$ forest surveys. Motivated by our findings in §4, i.e. a significant deviation from the fiducial model but a harder-to-resolve discrepancy between X-ray prescriptions, here we focus only on the Fast-Fid X-ray preheating model, which we recommend future work should do when accounting for X-ray preheating signatures present in the post-reionization IGM, namely to consider the fiducial X-ray preheating prescription of 21cmFAST. This focus herein is also due to the limitation in computational resources. For reference, we compare with a Type II version of our fiducial model, i.e. no X-ray preheating whatsoever.

For clarity, the only piece of our hybrid strategy being changed in this section is the substitution of the Type I $\psi$ results for the Type II transparencies. Thus, there is no modification for the “heating models” from MM20 that form our large box simulation suite.

The high-redshift range of the Ly$\alpha$ forest is often seen as an attractive probe of the nature of dark matter, particularly to constrain warm dark matter models (see e.g. Palanque-Delabrouille et al., 2020; Puchwein et al., 2022). However, our findings suggest yet another hurdle on the way for the quest for the nature of dark matter. Namely, the thermal state of the high redshift IGM not only depends on the relics from cosmic reionization, but through them, it also depends on the mechanism that set the initial conditions for cosmic reionization, i.e. X-ray preheating.

The precise mechanism through which the early intergalactic medium is preheated prior to cosmic reionization is highly uncertain, which arises, due to a lack of direct observational constraints, from the extrapolation of our understanding of the low-redshift Universe, particularly the X-ray luminosity functions (Lehmer et al., 2022), and the uncertainty regarding the abundance, clustering, and nature of the X-ray sources (Ma et al., 2021). Consequently, new avenues of probing the cosmic dawn are extremely valuable opportunities for discerning the astrophysics that governs this process. In this work, we establish, for the first time, the long-lasting impact of X-ray preheating in the Ly$\alpha$ forest, and highlight its sensitivity to different X-ray preheating models using numerical simulations.

X-ray preheating affects the Ly$\alpha$ forest in two main ways: (1) X-rays wipe out some of the small-scale structures that would be otherwise present, leading to a less transparent IGM at low redshifts. (2) Besides the “puffing up” of structures during the early reionization process, X-rays will penetrate deeply into the IGM raising its temperature and therefore setting up an earlier start to cosmic reionization. This early start also eventually leads to a decrease in the cross-correlation of matter and ionized bubble spatial structure.

Overall, the inclusion of X-ray preheating in the Ly$\alpha$ forest leads to a decrease in the memory of reionization, i.e. a decrease on the large-scale enhancement due to the impact of inhomogeneous H i reionization in the forest, at low redshifts with the total change depending on the specific X-ray preheating prescription. Furthermore, we find that the deviation from a model with no X-ray preheating can almost reach the $10\%$ mark for the 3D flux power spectrum of our Type I boxes at high redshift and for directions perpendicular to the line of sight, which also coincides with the sweet spot for the detection of the memory of reionization in the Ly$\alpha$ forest (Montero-Camacho & Mao, 2021). Likewise, the deviations in the already observed 1D flux power spectrum can lead to a change larger than $6\%$ at high redshifts. However, these quantitative findings are subject to the inability of Type I boxes to reliably sample the linear regime in the redshift range $2\leq z\leq 4$. Thus, we make a separate analysis with the aim to quantify the changes we should expect due to X-ray preheating in what we suggest should be implemented into future fiducial models of the Ly$\alpha$ forest interested in accurate IGM evolution in the post-reionization era.

Regarding the larger Type II boxes, our findings for the flux power spectra are encouraging for the current and next generation of Ly$\alpha$ forest surveys. At low redshifts, where the reionization relics are lesser and the forest is subject to other astrophysical systematics like UV clustering (Pontzen, 2014; Gontcho A Gontcho et al., 2014; Tie et al., 2019; Long & Hirata, 2023) and He ii reionization (Compostella et al., 2013; La Plante et al., 2017; Upton Sanderbeck & Bird, 2020), the long-lasting impact of X-ray preheating results in a decrease in the strength of the memory of reionization. In contrast, at large redshifts, where the impact of reionization is a dominant effect in the Ly$\alpha$ forest, X-ray preheating leads to a reinforcement of the signal. This enhancement may imply yet more optimistic outcomes when the Ly$\alpha$ forest’s ability to discern different X-ray preheating scenarios is forecasted or eventually measured.

During the completion of this work, Muñoz (2023) identified a potential issue in 21cmFAST simulations, namely that it underpredicts adiabatic fluctuations by approximately an order of magnitude due to the assumption that $T_{k}$ is homogeneous at $z=35$ (see their Figure 13). Although these findings can result in large effects on the 21 cm signal, we do not expect significant changes in the memory of reionization since X-ray preheating eventually takes over.

Our results underscore the necessity of including not only the memory of reionization in the Ly$\alpha$ forest but also the long-lasting effects of X-ray preheating in the fiducial model of Ly$\alpha$ analyses. Nevertheless, our findings also raise the same warning for other sensitive probes of the IGM in the post-reionization era. Thus, future work will investigate the impact of X-ray preheating in the imprints of H i reionization present in 21 cm intensity mapping (Long et al., 2022a).

We thank the anonymous referee for their insightful comments. This work is supported by the National SKA Program of China (Grant No. 2020SKA0110401), the Major Key Project of PCL, NSFC (Grant No. 11821303, 12050410236), and the National Key R&D Program of China (Grant No. 2018YFA0404502). PMC was partially supported by the Tsinghua Shui Mu Scholarship during the completion of this work. We are grateful to Catalina Morales-Gutierrez and Chris Hirata for their useful suggestions and comments. The authors acknowledge the Tsinghua Astrophysics High-Performance Computing platform at Tsinghua University for providing computational and data storage resources that have contributed to the research results reported within this paper. This work made extensive use of the NASA Astrophysics DataSystem and the following open-source python libraries: matplotlib (Hunter, 2007), numpy (Harris et al., 2020), and scipy (Virtanen et al., 2020).

The data underlying this article will be shared on reasonable request to the corresponding authors.