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

X-ray photons can penetrate deep into the intergalactic medium (IGM), leading to preheating of the IGM prior to cosmic reionization. X-ray preheating wipes out some of the small-scale structures that would otherwise be present prior to the passage of an ionization front. Accurate modeling of the small-scale structure is vital to the post-reionization IGM since the small-scale structure is ultimately the dominant source of long-lasting relics from hydrogen reionization. However, the precise impact of X-ray preheating in the fossils from hydrogen reionization is highly uncertain. In this work, we explore and establish for the first time, the long-lasting impact of X-ray preheating in the post-reionization IGM via hydrodynamic simulations with high-mass resolution. We find that the addition of X-ray preheating astrophysics leads to an overall lesser impact of the effect of inhomogeneous reionization in the Lyman-$\alpha$ forest -- depending on specific X-ray prescription -- at low redshifts ($z \sim 2$) with respect to a model with no X-ray preheating. However, at high redshifts ($z \sim 4$), our results indicate a strengthening of the relics of reionization in the Lyman-$\alpha$ forest because the IGM becomes more transparent compared to the scenario with no preheating. Thus, the absence of X-ray preheating in Lyman-$\alpha$ modeling can lead to a biased inference of cosmological parameters. Nevertheless, optimistically, the inclusion of X-ray preheating emerges as a promising novel avenue to probe the astrophysics of cosmic dawn.


INTRODUCTION
The Universe is currently believed to transition from primarily neutral to mainly ionized around  ∼ 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 ∼ 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 ( < 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 forest at low redshifts ( ≲ 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 ∼ 3 × 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  ∼ 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 (∼ 10 6 M ⊙ ) (Shapiro et al. 2004) but also suppress the growth of structure for cosmological time scales, eventually reaching a filtering scale of ∼ 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 signa-tures in the Ly forest at low redshifts ( ≤ 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 forest, i.e. the absorption features due to neutral hydrogen absorption of Lyman- 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)  1 .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 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 crosssections 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 ( ∼ 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- quasars and active galactic nuclei (Venkatesan & Benson 2011;Madau & Haardt 2015;Ross et al. 2019) 2 , 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 1 However, 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. 2 Nevertheless, 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  > 10.
should occur faster (Hirata 2018), which would lead to an overall suppression of the large-scale enhancement in both the Ly 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 ∼ 31 K as early as  = 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 ( ∼ 12).Regarding the latter effect, MM20 demonstrated that the acceleration or deceleration of the initial stages of reionization, i.e. increase or decrease of  HII ( ∼ 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 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.

PREHEATING OF THE IGM
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 Xray preheating since it is degenerate with reionization astrophysics (Ewall-Wice et al. 2016).However, the reionization relics in the postreionization 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 forest power spectrum via the impact of Xrays on the imprints from H i reionization that survive in the IGM to low redshifts ( < 4).Before covering our implementation of Xray 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 forest power spectra (M19; Oñorbe et al. 2019;Wu et al. 2019;Molaro et al. 2022)  3 .However, accurately modeling the imprint of reionization relics on the low-redshift post-reionization IGM ( ≤ 4) requires great mass resolution because the gas in underdense regions proves crucial.Gas in minivoids gets heated by the UV ionization front to ∼ 2 × 104 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 ( ∼ 2) Ly forest.
To compute the Ly 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 forest using a slightly unrealistic X-ray preheating model, where the transition from absorption to emission of the 21 cm brightness temperature occurs at  < 16, while the EDGES collaboration could have discovered the absorption peak of the signal centered around  = 17 (Bowman et al. 2018) 4 .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 Δ X gas () = 300 10 1 +  6.5 K . (1) 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  ′′ and  ′′ +  ′′ is in units of s −1 , where   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  ′′ ( ′′ +  ′′ ). * is the fraction of baryons that becomes stars,  coll is the collapsed fraction, Ω  is the baryon density,  crit,0 is the critical density (today), and R ′′ is the comoving distance between  ′ , the redshift of interest, and  ′′ . R ′′ nl is the evolved -with respect to redshift  ′′ -density smoothed on scales of R ′′ .The X-ray efficiency,   , 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  seen at (,  ′ ), from sources within  ′′ and  ′′ + ′′ , 21cmFAST assumes that the X-ray luminosity of sources follows a power law of the form,  ∝ (/ 0 ) −  , with  being a constant and  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  0 = ℎ 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  ′′ (see the detailed derivation in §3.1.2 of Mesinger et al. 2011).
The energy threshold,  0 , is the other free parameter in this Xray prescription that is allowed to vary in this work.A larger value corresponds to inefficient X-ray heating of the IGM due to more Xray 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 Figure 1.Gas temperature as a function of redshift for all X-ray preheating models considered in this work.We also include the CMB temperature evolution (dashed blue line).Moreover, the Hirata-300 model is also shown (magenta dotted line).The heating transition occurs between 10.7 ⪅  ⪅ 13.0 for most models.
IGM temperatures that we consider throughout this work.In summary, the "Fast-Fid" model has  0 = 500 eV and   = 2×10 56 M −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  0 = {100, 1000, 1500} eV and three additional variations for the X-ray ionization efficiency   = {1, 4, 8} × 10 56 M −1 ⊙ .For the nomenclature, we refer to the "E2" model as the one with  0 = 1000 eV, the "  3" model has   = 8 × 10 56 M −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   3 model has the earliest transition from 21 cm absorption to emission ( ≈ 15).The bulk of models experiences the heating transition from 10.7 ⪅  ⪅ 13.0, while the latest model to surpass the CMB temperature is the E3 model where it happens at  = 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  = 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   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  ≈ 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 +  > 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, Δ X gas is one (two) order of magnitude lower than  CMB at  = 50 ( = 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.

SIMULATIONS
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(Springel et al. , 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 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 × 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 × 10 3 M ⊙ , and a gas-particle mass of 2.27 × 10 2 M ⊙ .The Type II simulations have a box size of 1275 kpc on each side, 2 × 192 3 total number of particles, streaming velocity of 33 km/s , a dark matter particle mass of 9.72 × 10 3 M ⊙ , and a gas-particle mass of 1.81 × 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 forest survey.In contrast, the Type II boxes are large enough to form multiple Ly 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  B is the Boltzmann's constant,  is the total number density of particles,  is the gas temperature,   is the net volumetric heating rate in erg cm −3 s −1 from process ,  =  −5/3 is the particle entropy in the code, and   / is the fractional rate of increase of thermal energy from process .
To include X-ray preheating into the Gadget-2 simulations, it is sufficient to use Eq. ( 5) to account for the extra heating Δ due to X-rays, i.e.   / ∝ Δ  /.Thus, one needs to replace Δ  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   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   3 model.In contrast, the underdense areas have been severely diminished.Regarding the differences between the two Xray 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 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 Δ 2   () where   is the Case A recombination rate and Δ = 1 +  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  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  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  1 value as a function of the redshift of observation and local redshift of reionization, and what we aim to identify is how  1 changes with the different heating models described in Table 1.

RESULTS WITH TYPE I SIMULATIONS: IMPACT OF X-RAY PREHEATING ON IGM
MM20 demonstrated that while X-ray preheating can affect the initial stages of the reionization process, e.g. by raising the global  HII to ∼ 0.1 before UV heating takes over, this effect has a minimal impact on the imprints from inhomogeneous reionization in the postreionization IGM, and therefore the resulting memory of reionization in the Ly 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 forest because the temperature-density relation becomes bimodal due to the underdense regions suffering shockheating and compression to mean density in addition to UV heating (Hirata 2018).The impact of X-ray preheating on the temperaturedensity 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 flux transmission, including the effect of reionization, are given by (M19) where  F is the flux bias, which is negative since the presence of matter implies absorption,  F is the redshift-space distortion parameter,  is the cosine of the angle of the  vector with respect to the line z = 8 Fiducial z = 6 Fiducial Hirata-300 Hirata-300 Hirata-300 Often in Ly simulations, one varies the normalization -in this case,  1 -which is equivalent to varying the ionizing background, until one obtains the correct mean transmitted flux in our simulations.Specifically,  1 is the optical depth that must be assigned to a patch of gas with mean density and temperature  = 10 4 K in order for the mean transmitted flux to match observations. parametrizes the variations of transparency for an opacity cube that suddenly reionizes at  re , and is observed at  obs , relative to an opacity cube in the reference simulation with no X-ray preheating and that reionizes at  re = 8.Note that  < 0 implies that the reference simulation is more transparent or equivalently has a higher F. We plot the relative transparency of the IGM for our Type I highresolution 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  X−ray <  Fid 5 for early reionization scenarios ( 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  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 clouds inside them leading to not-so-accurate sampling.Reassuringly, the different Xray models, including the simple Hirata-300 model and the Fiducial model, converge for very early reionization ( = 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  re = 12 has the shortest response time.
As seen in Figure 3, the same general trend from the early reionization scenarios,  X−ray <  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  re ≈ 7.99, while the last model is the   3 model that crosses the threshold at  re ≈ 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 B1, including redshifts that were omitted in the figure.From Table B1, the largest discrepancy between the X-ray prescriptions and the Fiducial model occurs at  re = 6 and  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  re = 8), compared to 3.489% more transparent in the   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 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.
5 This is easier to appreciate in Table B1.

3D Flux power spectrum
The Ly flux power spectrum can be computed from Eq. ( 6) by where the first term of the r.h.s. is the conventional Ly forest power spectrum  3D F , the second of the r.h.s. is the contribution due to the memory of reionization  3D re , and we neglect higher order terms of .Here,  ,  is the cross-power spectrum of matter and transparency of the IGM (see below for the explicit form),  L m is the linear matter spectrum, which we compute using CLASS (Lesgourgues 2011;Blas et al. 2011), and the non-linear correction  NL is given by (Arinyoi-Prats et al. 2015) ln which takes care of the isotropic non-linear enhancement (first term on the r.h.s. of Eq. ( 9)), Jeans smoothing6 , i.e. gas pressure isotropically suppresses the power to below the Jeans scale (third term therein), and the suppression due to peculiar velocities along the line of sight (second term therein).Δ 2  is the dimensionless linear matter power spectrum and the different free parameters are taken from Tables 8 and 9 of Arinyo-i-Prats et al. ( 2015), corresponding to simulations with box size  = 60 ℎ −1 Mpc.They also run various simulations to test the convergence of parameters when varying box size and resolution.According to their Tables 2 and 3, when doubling the box size, the changes are within 5 percent for  F and  F and within 10 percent for non-linear correction parameters except Jeans scale   (about 20 percent).The behaviour is similar when increasing the resolution by a factor of (5/3) 3 , but due to the limited variation of simulations, their results do not show clear convergence with resolution.However, our analysis of the impact of X-ray preheating on the Ly power spectrum is unlikely to be significantly affected by the parameter uncertainties.A rudimentary estimate indicates that the fractional uncertainty is approximately the fractional uncertainty of the parameters, thus about several percent since the X-ray preheating impact mainly arises on large scales.The main source of uncertainty should still be the transparency results from our small-box simulations7 and the uncertainty in the timeline of reionization.
The range of redshift in the Arinyo-i-Prats et al. ( 2015) tables is from 2.2 to 3.0, so for  > 3 calculations, we introduce a redshift evolution factor [(1+)/(1+ * )] 3.55 (Palanque-Delabrouille et al. 2013) and use a pivot redshift  * = 3.For the conventional Ly forest power spectrum, we have  3D F ( > 3) = [(1 + )/(1 +  * )] 3.55  3D F ( * ).Analogously, to calculate the memory of reionization term  3D re , the redshift evolution of the flux bias and redshift-space distortion parameter, for  > 3, is given by Note that similar bias evolution strategies are often used when dealing with Ly forest observations (see e.The second term of Eq. ( 8) corresponds to the memory of reionization, i.e. the imprints due to H i inhomogeneous reionization in the Ly forest, and the key physical quantity is the cross-power spectrum of matter and transparency of the IGM, which parametrizes the role of inhomogeneous reionization -primarily, the void response to the reionization process.It can be written as (see M19 for the derivation) One final caveat before computing the long-lasting effect of X-rays in the IGM is that we include a biasing procedure when computing  ,  at the largest scales.The need for this procedure is clear because our 21cmFAST simulations only have a box size of 400 Mpc, and so there are just so few modes that can fit in the first few bins when the separation between the bins is Δ = 2/400 Mpc −1 .Therefore we implement a cutoff at  cut = 0.044 Mpc −1 .At  <  cut , we model  ,  using a linear biasing model, which should be accurate at scales larger than the ionization bubble scale.
We chose to quantify the impact of X-ray preheating in the 3D Ly power spectrum in terms of the fractional difference of model "x" with respect to the Fiducial model, The change of pressure smoothing with respect to X-ray preheating is negligible for non-small scales.Consequently, only the reionization term appears in the numerator of the fractional difference.In Figure 4, we plot the fractional difference using Eq. ( 13) for the different X-ray models at different angles with respect to the line of sight and at two different redshifts.(Note that we drop the Hirata-300 model from the analysis at this point since there is no simple way of implementing a 21cmFAST model that has an analogous evolution of the gas temperature.)At low redshifts, Figure 4 indicates the maximum deviation from the Fiducial model occurs for the E1 model with a roughly nine percent decrease in the signal, although there is intense competition with the   3 model.Meanwhile, the minimum variation is encountered in the E3 model at large  where the impact is quantified at less than a percent.The effect of X-ray preheating is also more substantial in directions perpendicular to the line of sight, which remarkably aligns with the () sweet spot to recover the memory of reionization in the Ly forest because the reionization term -second term in Eq. ( 8) -has only one Kaiser term compared to the conventional Ly term that has two Kaiser factors.Overall, the hierarchical structure between the X-ray preheating models can also be appreciated in the impact of the 3D flux power spectrum, with the E3 model being the one with the lowest variation and the   3 and E1 models competing for the largest one.Note that this trend is consistent with that observed in the transparency of the IGM at  obs = 2 (see Figure 3) and consistent with the gas temperature (Figure 1) once one accounts for the fact that the extra heating in the photon-abundant model (E1) has a lower X-ray efficiency compared to that of the   3 0.0 0.2 0.4 0.6 0.8 1.0 0.10 0.08 0.06 0.04 0.02 0.00 model.Hence the constant early rise of gas temperature is slightly less important than the sudden late rise since there is not enough time to completely influence the transparency results, especially at higher redshifts.Interestingly, the situation "normalizes" when observation is at a low enough redshift.Naturally, the scenario is likely to differ if one observes close to the tail end of reionization because the impact of the large amplitude of the gas temperature in the late-rise scenario might become dominant, which can already be seen in the fractional difference at high redshifts of Figure 4.
Analogously, at high redshift, we see the same trend with respect to .Now the most significant difference clearly happens for the E1 model and corresponds to a decrease of roughly 10 percent.However, there is fierce competition and the   3 model has the largest deviations at some small scales.This sudden takeover at large scales is likely due to the large impact of the last snapshot (for the local reionization scenarios, i.e.  re = 6) in the integrand of Eq.( 11) -see also Figure 1 of MM20.The smallest change still corresponds to the E3 model in directions close to the line of sight, even though the E2 model can take the last spot at the smallest wavenumber, given the error bars this last feature is likely not statistically significant.
Although the amplitude of the fractional difference seen in Figure 4 for our Type I boxes is somewhat expected because the Fiducial scenario disregards any X-ray preheating physics, the spreads between the X-ray preheating prescriptions are surprising.For instance, at  = 0.1 Mpc −1 and  = 0.1 the models span the range of ∼ 0.05  3D,Fid at  obs = 2 and ∼ 0.08  3D,Fid at  obs = 4.However, note that the significance of the spread is slightly diminished if one were to include X-ray preheating physics in the Fiducial model.In this case, the Fast-Fid model would be the reference point and therefore the difference between the X-ray models will be slightly lesser but still quite significant.For now, the choice of this new Fiducial scenario is a good compromise between modeling accuracy and computational resources for future studies aiming at quantifying the response of the IGM to the reionization process.The sensitivity to X-ray preheating opens a novel window into the astrophysics that governs the cosmic dawn billions of years after it happened.

1D Flux power spectrum
As the complement to the 3D flux power, the Ly forest 1D power spectrum requires the cross-correlation between pixels along the same line of sight.Mathematically, the expression is given by integrating the 3D counterpart over the perpendicular direction to the line of sight (Palanque-Delabrouille et al. 2013; see also Eq. 10 of M19) and can be written as the sum of the conventional Ly forest power spectrum  1D F and the contribution due to the impact of inhomogeneous reionization  1D re .Here, we implement a cutoff at  max = 7 Mpc −1 .The arbitrary choice of this cutoff is not expected to change our results since the H i reionization relics are coupled to the ionization bubble scales and quickly decline as the wavenumber increases, and so it is close to zero at  = 1 Mpc −1 .
Analogous to Eq. ( 13), the 1D fractional difference is We assume that the impact of X-ray preheating on the pressure smoothing is small at the redshift range and scales of interest8 .Hence, the fractional difference for the 1D power spectrum has the same form as the 3D counterpart.In Figure 5, we plot the fractional difference in the 1D flux power spectrum for the different X-ray preheating models considered by our Type I simulations.At low redshifts, the 1D fractional difference recovers the hierarchical structure observed in Figure 4 with the largest deviation (< 4% in absolute value) caused by the   3 model while the smallest effect (< 1%) is present in the E3 model.At high redshifts, we observe a similar order of the models -albeit with fierce competition between E1 and   3 models at large scales -and a larger spread in the values of the fractional differences that is, crucially, comparable to the strength observed in the 3D findings.In particular, we obtain a maximum deviation larger than 6 percent for both the E1 and   3 models.From an optimistic point of view, this is a promising avenue for learning about the astrophysics that governs the cosmic dawn since medium and high-resolution P1D measurements and spectra already exist (Chabanier et al. 2019;Murphy et al. 2019;O'Meara et al. 2021).On the other hand, this implies that the Ly forest can produce biased results in the absence of proper countermeasures against this broadband systematic.However, as shown in Montero-Camacho et al. (2023), this effect will not impact the accuracy of the location of the BAO peak at any redshift but can affect measurements dependent on the shape quite drastically, which could, for instance, jeopardize recent efforts to weaponize the AP effect (Alcock & Paczynski 1979) to extract additional cosmological information from Ly two-point statistics (Cuceu et al. 2022(Cuceu et al. , 2021)).
Furthermore, the 1D fractional differences have a distinct morphology with respect to the 3D analogs.This trend is easily understood via Eq.( 14) --the integration over perpendicular directions to the line of sight leads to mode-mixing; hence small wavenumbers, which suffer more greatly from the impact of reionization and therefore of X-ray preheating, can be mapped into large wavenumbers.This shape is also encouraging, particularly for future efforts aimed at constraining X-ray preheating astrophysics (e.g.properties of the X-ray sources in the high- IGM) through the impact of inhomogeneous reionization in the 1D flux power spectrum.This shape is also the reason why the   3 model has an easier time at both redshifts and most scales compared to its rival the E1 model which tends to rise at large scales in  3D .
Given that  1D data already exists, future work will assess the potential implications of considering X-ray preheating as part of the fiducial model for inference of cosmological parameters.For now, we highlight that the values seen in Figure 5 are larger than the statistical errors of recent Ly forest surveys (Chabanier et al. 2019).In particular, the statistical error at  = 0.1 Mpc −1 is 0.0033 Mpc at  = 2.2 and 0.0731 Mpc at  = 4.In comparison, the effect of X-ray preheating in the IGM (i.e.no X-ray preheating vs the   3 X-ray prescription) at the same wavenumber is 0.0055 Mpc at  = 2.2 and 0.0654 Mpc at  = 4.
Finally, we recognize that the results obtained in §4.1 and in §4.2, which use our Type I simulations, suffer from poor sampling given the small box size.Nonetheless, they are sufficient to inform us of the general trends between the different X-ray prescriptions, particularly to highlight the ability to discern between different cosmic dawn astrophysics.To accurately capture the long-lasting impact of X-ray preheating in the flux power spectra, we include the analysis of the larger, Type II boxes in §5, where we only consider the Fast-Fid X-ray prescription and compare it with the Fiducial model due to limited availability in computational resources.

RESULTS WITH TYPE II SIMULATIONS: ACCURATE CALCULATION FOR THE FAST-FID MODEL
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 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 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  results for the Type II transparencies.Thus, there is no modification for the "heating models" from MM20 that form our large box simulation suite.

Impact on IGM evolution
We tabulate the relative transparency as a function of the redshift of observation and of the redshift of local reionization, which are the data used in Figure 6, in Table B2.Note that the reference scenario is now a Fiducial model Type II box where the local reionization occurs at  re = 8.Furthermore, we showcase the transparency results for the  Type II boxes in Figure 6. Figure 6 shows that the Fast-Fid model, which uses the default prescription for X-ray preheating in 21cmFAST (Mesinger et al. 2011), is more transparent than the Fiducial model at high redshift of observation and for both late and early reionization scenarios; nevertheless, for early reionization scenarios, the Fiducial and the Fast-Fid model are both less transparent than the reference, which is the same trend observed for the Type I boxes in Figure 3.
In contrast to the Type I results, Figure 6 showcases a transition from more transparent to less transparent for the Fast-Fid model as time evolves, at least within the error bars.Hence, we see a recovery of the trend found in Figure 3 at low redshift of observation and for late reionization scenarios, albeit with slightly lower amplitude for the Fiducial model.For early local reionization, the tendency is obscured due to the Monte Carlo scatter between the different realizations.
Interestingly, the Type II boxes reach considerably larger values than the Type I counterparts at the high redshift of observation.Besides, / re achieves considerably larger values (in absolute magnitude) than for the Type I boxes at late  re .This feature is of significant importance because of the redshift evolution nature of the reionization relics in the forest, where the highest redshifts lead to a stronger signal since the HEMD gas has not had enough time to relax the additional injected energy during the reionization process (Hirata 2018).Furthermore, the Type II boxes have different behavior at the high end of the local reionization redshift tail.We attribute these features to the increased statistics and more reliable sampling of the linear regime that the Type II boxes achieve compared to that of the Type I analogs.In particular, better sampling allows the Type II boxes to more accurately capture the impact of the initial larger (high  obs ) gas temperature and, consequently, the wiping-out of small-scale structures.In contrast, at low  obs , the initial higher temperature -and changed densities -has an easier time relaxing to the temperature-density relation because there should be a lesser HEMD phase.
The maximum deviation, quantified as the absolute difference between the two models, present in our relative transparency results is 3.465% and occurs at  re = 6 and  obs = 4, which is smaller than the one predicted for the Type I models (4.518%).However, this was achieved by the   3 model at  re = 6 and  obs = 3, and thus provides an unfair comparison.On the other hand, the minimum deviations tend to occur at early reionization scenarios (and at low redshifts of observation), which is also the trend seen in Table B1.

Impact on flux power spectra
Because of the enhancement at high  obs in Figure 6, we now expect the impact on the Ly forest power spectrum lead to a rise on the X-ray heating signatures at the high redshift of observation, and eventually result in a decrease on the strength of the effect at low redshifts, which is a similar see-saw pattern to the one that is expected in the 21 cm power spectrum (Pritchard & Furlanetto 2007).The cause of these see-saw mechanisms is essentially the same: X-ray photons will first increase the power and raise the temperature; however, this would eventually lead to a decrease in power since it comes at the cost of erasing fluctuation in small scales.
We plot the fractional difference of the 3D flux power spectrum for the fiducial 21cmFAST model ("Fast-fid" model) with respect to the model with no X-ray preheating in Figure 7. Figure 7 includes the standard deviation from the mean due to the use of four different realizations of the small-scale high-mass resolution simulations.
At low redshift of observation, Figure 7 shows the trend of increasing the impact (in absolute value) at small values of , which is the same trend observed in Type I boxes (see Figure 4).Nonetheless, Type II boxes have a slightly diminished impact, with a maximum deviation -at large scales and along the perpendicular direction -  3 but for the Type II boxes and including the Monte Carlo scatter due to different realizations.We only consider the overall fiducial model with no X-ray preheating (blue solid line) and the standard 21cmFAST X-ray prescription, "Fast-Fid" model (green dotted line).Interestingly, in contrast to Figure 3, there is now a turnover with the X-ray model leading to larger transparency initially at high redshifts.
around three percent.Do note that a potential Type II   3 model (or even an E1 model) would give rise to an enhancement in the relative strength of the signal.However, the purpose of this section is to establish what can be expected for a fiducial model of the long-lasting impact of X-ray preheating in the Ly forest.Thus, the priority is not on the intriguing ability of this effect to constrain the astrophysics that governs cosmic dawn.
On the other hand, at high redshift, Figure 7 shows a strengthening of the memory of reionization -the imprints from H i reionization in the Ly forest -when accounting for the effect of X-ray preheating in the IGM.The effect is reinforced at large scales and for small values of , which corresponds exactly to the best window of opportunity to directly capture the epoch of reionization astrophysics from the Ly forest (Montero-Camacho & Mao 2021).In particular, we obtain an enhancement of more than six percent in this window.This feature is unique to Type II boxes since it originates on the larger transparency of the X-ray preheating model with respect to the reference scenario at high  obs and for late local reionization.
Likewise, Figure 8 illustrates the effect of X-ray preheating in the 1D Ly flux power spectrum as a function of wavenumber and both for the low and high redshift of observation.Similar to Figure 7 the Monte Carlo scatter from different realizations of the highresolution Gadget2 simulations has been included.The mean result at low redshift is consistent with the values obtained for Type I boxes (see Figure 5), although it is also diminished, analogously to the 3D effect.We see a decrease in the signal, with respect to a scenario with no X-ray preheating, peaking at less than two percent.
In addition, at high redshift, we observe a moderate increase in the signal at large scales of less than four percent while a more modest enhancement is also perceptible at smaller scales (slightly larger than two percent).Again, this morphology occurs because of the integration with respect to the directions perpendicular to the line of sight, which can map small  ⊥ modes into large  and vice versa.In comparison with Figure 5, the relative impacts seen in Figure 8 may be lesser but still significant, particularly due to the enhancement of the memory of reionization in the regime where it is the strongest and due to its relative strength when compared with current statistical error (Chabanier et al. 2019) and compared to that of near-term measurements (e.g., Ravoux et al. 2023;Karaçaylı et al. 2023).

SUMMARY
The high-redshift range of the Ly 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 forest, and highlight its sensitivity to different X-ray preheating models using numerical simulations.
X-ray preheating affects the Ly 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 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 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 ≤  ≤ 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 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 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 forest, X-ray preheating leads to a reinforcement of the signal.This enhancement may imply yet more optimistic outcomes when the  Ly 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   is homogeneous at  = 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 forest but also the long-lasting effects of X-ray preheating in the fiducial model of Ly 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).

APPENDIX A: TEMPERATURE-DENSITY RELATION
Here we focus on the effect of X-ray preheating as seen in the temperature-density relation and its evolution.
It is expected that the heated gas from cosmic reionization would eventually recover a tight temperature-density relation  ∝ Δ −1 (Hui & Gnedin 1997;McQuinn & Upton Sanderbeck 2016), where  is the gas temperature and Δ = / ρ is the ratio of the density and mean density.However, the bimodal nature of the temperaturedensity relation, caused by the way the small-scale structure reionized (Hirata 2018), requires cosmological timescales to relax the additional injected energy during the reionization process.Thus, the recovery of the temperature-density relation, even for models with X-ray preheating, is delayed until  ∼ 2 9 .Given the strength of Xrays on the transparency of the IGM, we expect that differentiating between different X-ray prescriptions would be difficult.Thus we opt to showcase only the Fiducial model here and compare it with the   3 model since the latter showcases the largest deviations for late reionization scenarios.
Figure A1 and Figure A2 illustrate the impact of X-ray preheating in the bi-modal nature of the temperature-density relation for patches of gas that have late local reionization ( re = 6) and local reionization consistent with Planck Collaboration et al. 2018 ( re = 8), respectively.In both figures, the Fiducial model (no X-ray heating) leads to a more spread out distribution of the gas, i.e. a less tight blue contour, hence a stronger High-Entropy Mean-Density (HEMD; Hirata 2018; MM20) phase of the temperature-density relation, than the   3 model (shown in the bottom panel).As time evolves, the HEMD gas, which originates in the underdense regions that were heated by 9 Naturally He ii reionization will also disturb the temperature-density relation.However, no high-mass resolution analog study has been done to investigate the way the small-scale structure reionizes.Hence, the study of the memory of He ii reionization and its impact on the H i relics will be left to explore in future work.shocks and UV radiation (and compressed to mean density), traces the evolution of gas in denser regions.Thus, we see a partial recovery of the temperature-density relation at  obs = 2 for the late local reionization scenario and a more complete recovery for the Planck local reionization scenario.This recovery is more pronounced in the model with X-rays because some small-scale structures that would otherwise be present were smoothed out.

APPENDIX B: IGM TABLES
We tabulate the relative transparency, as a percentage, of our Type I and Type II boxes with respect to a reference scenario with local reionization occurring at redshift 8 in Tables B1 and B2, respectively.We highlight that there can be significant deviations between the results for Δ ln  1 for the Type I and II boxes.For instance, at  re = 9 and  obs = 2.5 the Fast-Fid and Fiducial models report an approximately 80 percent difference (ignoring the error bars).In contrast, at  re = 7 and  obs = 2.5 the Fast-Fid and Fiducial models have a deviation of 6 and 15 percent, respectively.The discrepancy between boxes arises from the difficulty the Type I box has to capture the forest due to poor sampling.However, we acknowledge that convergence for the transparency result has not been achieved even with the larger Type II boxes.
Our Type II boxes correspond to Hirata (2018) II-C boxes and have a missing variance of 1.32 at  = 2.5 (see their §5.3),which is an improvement from the 2.67 missing variance of the Type I boxes that have a factor of 9 smaller volume.Nevertheless, convergence in Δ ln  1 for even the II-F box ( = 2551 kpc, 216 times the Type I volume, and a missing variance of 0.77) used in M19 has not been achieved.The necessary computational resources to achieve satisfactory convergence in Δ ln  1 are too restrictive, particularly given the significant number of simulations needed to compute our variance, which is why we have chosen to use our Type II boxes for this analysis.
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 2 .
Figure 2. Log-normal gas density distribution of representative snapshots from the high-resolution Type I simulations, which have a box size of 425 kpc.The left, middle, and right columns correspond to redshifts  = 6, 8, and 10, respectively.The top, middle, and bottom rows correspond to the Fiducial model (no X-ray preheating of the IGM), the Hirata-300 model (with X-ray preheating prescription in Eq. 1), and the   3 model (with the X-ray preheating obtained using the 21cmFASTv1.3prescription), respectively.

Figure 3 .
Figure 3. Transparency of the IGM  = Δ ln  1 for our Type I simulations, as a function of local reionization redshift, relative to the scenario where local reionization redshift is 8 in the Fiducial model.Left, middle, and right columns correspond to different redshifts of observation  obs = 2.0, 3.0, and 4.0, respectively.For visual reference, we highlight the level of the reference scenario with a dashed magenta line.

Figure 4 .
Figure 4.The impact of X-ray preheating in the 3D Ly power spectrum -the fractional difference of the Ly power spectrum in different X-ray preheating models with respect to that in the Fiducial model.Shown are the results at three different  (the cosine of the angle with respect to the line of sight), at low redshift  obs = 2 (top) and at high redshift  obs = 4 (bottom), respectively.The largest deviation occurs at high redshifts and in the direction nearly perpendicular to the line of sight (i.e.small ).

Figure 5 .
Figure5.The impact of X-ray preheating in the 1D Ly power spectrum -the fractional difference of the Ly power spectrum in different X-ray preheating models with respect to that in the Fiducial model.Shown are the results at low redshift  obs = 2 (left) and at high redshift  obs = 4 (right), respectively.

Figure 6 .
Figure 6.Same as Figure3but for the Type II boxes and including the Monte Carlo scatter due to different realizations.We only consider the overall fiducial model with no X-ray preheating (blue solid line) and the standard 21cmFAST X-ray prescription, "Fast-Fid" model (green dotted line).Interestingly, in contrast to Figure3, there is now a turnover with the X-ray model leading to larger transparency initially at high redshifts.

Figure 7 .
Figure 7. Same as Figure 4 but for the Type II boxes with only the standard 21cmFAST X-ray prescription, "Fast-Fid" model and including the Monte Carlo scatter due to different realizations.The green line represents the mean result from the different realizations.

Figure 8 .
Figure 8. Same as Figure 5 but for the Type II boxes with only the standard 21cmFAST X-ray prescription, "Fast-Fid" model and including the Monte Carlo scatter due to different realizations.The green curve shows the mean result.

Figure A1 .
Figure A1.The evolution of temperature-density relation in the post-reionization era for the latest local reionization redshift ( re = 6), which is consistent with recent constraints for the tail end of reionization (see, e.g.Keating et al. 2020).The gas particles have been smoothed using a Gaussian density estimator with an FWHM of 0.05 dex on each axis.The color scale shows the mass-weighted distribution of gas in units of log 10 probability per log 10 Δ per log 10 .Contours are shown for every decade in probability density.Shown are the results at the redshift of observation  obs = 2.0 (left), 4.0 (middle), and 5.5 (right), for different X-ray preheating models -Fiducial model (i.e.no X-ray preheating; top row), and the   3 model (i.e.X-ray preheating is taken from the 21cmFAST prescription; bottom row), respectively.

Figure A2 .
Figure A2.Same as Figure A1 but with the redshift of local reionization  re = 8 (in line with the constraints from Planck Collaboration et al. 2018).
,  (,  obs ) = −   is the growth rate,  ,  HI is the cross-power spectrum between matter and neutral fraction field, which reflects the ionized bubble spatial structure and is needed to account for the patchy nature of reionization.Here, / re quantifies the response of the IGM to the reionization process, i.e. how the transparency of the IGM changes as the local redshift of reionization changes.Note that each preheating model gives us a unique  ,  (,  obs ) that accounts for the evolution of both  ,  HI (,  re ) and / re ( re ,  obs ) during the epoch of reionization.