Impact of Self-shielding Minihalos on the Lyα Forest at High Redshift

Dense gas in minihalos with masses of 106−108 M ⊙ can shield themselves from reionization for ∼100 Myr after being exposed to the UV background. These self-shielded systems, often unresolved in cosmological simulations, can introduce strong absorption in quasar spectra. This paper is the first systematic study on the impact of these systems on the Lyα forest. We first derive the H i column density profile of photoevaporating minihalos by conducting 1D radiation–hydrodynamics simulations. We utilize these results to estimate the Lyα opacity from minihalos in a large-scale simulation that cannot resolve self-shielding. When the ionization rate of the background radiation is 0.03 × 10−12 s−1, as expected near the end of reionization at z ∼ 5.5, we find that the incidence rate of damped Lyα absorbers increases by a factor of ∼2−4 compared to at z = 4.5. The Lyα flux is, on average, suppressed by ∼3% of its mean due to minihalos. The absorption features enhance the 1D power spectrum up to ∼5% at k ∼ 0.1 h Mpc−1 (or 10−3 km−1 s), which is comparable to the enhancement caused by inhomogeneous reionization. The flux is particularly suppressed in the vicinity of large halos along the line-of-sight direction at separations of up to 10 h −1 Mpc at r ⊥ ≲ 2 h −1 Mpc. However, these effects become much smaller for higher ionizing rates (≳0.3 × 10−12 s−1) expected in the post-reionization Universe. Our findings highlight the need to consider minihalo absorption when interpreting the Lyα forest at z ≳ 5.5. Moreover, the sensitivity of these quantities to the ionizing background intensity can be exploited to constrain the intensity itself.


Introduction
Roughly between 10 8 and 10 9 yr after the Big Bang, the intergalactic medium (IGM) was reionized by UV radiation from early galaxies (see reviews by Loeb & Barkana 2001;Wise 2019).During this period, the initially cold and neutral IGM transitioned to a hot plasma of ∼20,000 Kelvin, making it impossible for dark matter halos with masses less than 10 8 M ⊙ to gravitationally accrete gas and form baryonic structures.Since then, 10 8 M ⊙ has served as the minimum mass scale for the formation of large-scale structures within the IGM.
Before reionization, however, the cold neutral IGM was capable of forming structures below 10 8 M ⊙ .The first stars and galaxies are believed to have formed in minihalos (MHs) with masses ∼ 10 6 M ⊙ at z ∼ 20 − 30 (Tegmark et al. 1997;Abel et al. 2002;Bromm et al. 2002;Yoshida et al. 2006;Bromm & Yoshida 2011;Loeb & Furlanetto 2013).After these first galaxies photodissociated hydrogen molecules in the IGM (Haiman et al. 1997), subsequent MHs could still accrete gas but were unable to form stars unless they grew massive enough (≳ 10 8 M ⊙ ) to excite atomic hydrogen and enable cooling (e.g., Benitez-Llambay & Frenk 2020).Consequently, baryons within MHs that formed after the first stars likely existed as non-star-forming gas clouds.During reionization, numerous such clumps were present in the intergalactic space.
Once reionization occurs, the Jeans mass increases to ∼ 10 8 M ⊙ , leading to the gradual destruction of smallscale structures due to the increased pressure of photoionized gas (e.g., Park et al. 2016;D'Aloisio et al. 2020;Park et al. 2021a;Puchwein et al. 2023;Chan et al. 2023).However, dense cores of MHs can remain neutral and serve as the sinks of the ionizing background for more than 100 Myr (Shapiro et al. 2004;Iliev et al. 2005a,b;Nakatani et al. 2020).These self-shielded systems restrict the mean free path of the ionizing photons, impeding the growth of ionized bubbles (e.g., Nasir et al. 2021), and influencing the large-scale morphology of reionization (Gnedin & Fan 2006;Choudhury et al. 2009;Alvarez & Abel 2012;Bianco et al. 2021;Cain et al. 2023).Moreover, the Lyα opacity of these systems can attenuate the Lyα emission originating from Park et al. nearby star-forming galaxies (Park et al. 2021b;Smith et al. 2022).
The evaporation time of the MHs is highly dependent on their mass.Relatively massive MHs with a mass around 10 8 M ⊙ retain a significant portion of their baryons even after reionization (Iliev et al. 2005a;Nakatani et al. 2020), and they are associated with Lyman limit systems in the post-reionization Universe (Storrie-Lombardi et al. 1994;Miralda-Escudé 2003;Songaila & Cowie 2010;Prochaska et al. 2010Prochaska et al. , 2015)), while MHs with less than 10 6 M ⊙ lose most of their gas within ∼ 10 7 yr.The intermediate systems, ranging from 10 6 to 10 8 M ⊙ , undergo photoevaporate over a period of ∼ 10 8 yr.The number density of these selfshielded systems is believed to have evolved rapidly between z = 5.5 and 4.5 as less massive ones evaporate earlier than the more massive ones do.
The Lyα forest, spectral features in spectra of distant quasars due to intervening Lyα absorbers, is the most effective means of probing intergalactic structures (for a recent review, see McQuinn 2016).Therefore, in this study, we will explore the potential impact of these shielded MHs on the statistical properties of the Lyα forest near the end of reionization.Without star formation, the MHs would follow the truncated isothermal sphere profile until they are exposed to the ionizing background.The neutral hydrogen column density of these objects can even exceed 2 × 10 20 cm −2 at the cores, resulting in the damping-wing opacity of atomic hydrogen casting an extended shadow of ≳ 10 Mpc in quasar spectra.Given the high number density of the intermediate-mass (10 6 − 10 8 M ⊙ ) MHs and the challenges in subtracting damped Lyα absorbers (DLAs) from the high-z Lyα forests due to low average flux, the damping-wing absorption by the self-shielded systems can significantly impact the statistics of the high-z Lyα forest.While the absorption features of these MHs in the 21cm forest were explored by Furlanetto & Loeb (2002), this study is the first to focus on their impact on the Lyα forest.
Recently, there has been growing attention toward the high-z Lyα forest as a promising future probe of reionization (Fan et al. 2006a).In addition to UV background fluctuations, inhomogeneous reionization induces largescale thermal fluctuations, which leave more lasting observable signatures on the distribution and power spectrum of the Lyα flux even after the end of the reionization (z ≳ 5; e.g., Songaila & Cowie 2002;Nasir et al. 2016;Oñorbe et al. 2017aOñorbe et al. , 2019;;Wu et al. 2019;Molaro et al. 2022Molaro et al. , 2023;;Puchwein et al. 2023).The number of quasars discovered at such high redshifts has been steadily increasing in recent yr (Fan et al. 2006b;Becker et al. 2007;Willott et al. 2010;Mortlock et al. 2011;Venemans et al. 2013;Becker et al. 2015;Wu et al. 2015;Jiang et al. 2016;Bañados et al. 2018;Yang et al. 2020), and their spectra can be utilized to constrain the details of reionization if reionization is modeled accurately (e.g., Lidz et al. 2006;Bolton & Haehnelt 2007;Oñorbe et al. 2017b;Qin et al. 2021).Additionally, Hirata (2018) found these signals can theoretically be detected down to z ∼ 2. Thus, a substantial number of quasars to be observed at z ≲ 4.5 by large-scale experiments such as the Dark Energy Spectroscopic Instrument (DESI) survey could also prove useful in this context (Montero-Camacho et al. 2019;Montero-Camacho & Mao 2021).
If the absorption by MHs affects the statistics of the high-z Lyα flux, it is crucial to consider this factor when analyzing the Lyα forest to constrain reionization.The amount of remaining neutral gas in MHs depends on the timing and intensity of the ionizing background.Therefore, the presence of MH absorption in the Lyα forest may provide extra information about reionization.In this study, we shall explore these possibilities.
Due to the significant difference in length scale between MH photoevaporation (∼ ckpc) and the Lyα forest (∼ 100 cMpc), it is computationally prohibitive to simulate the Lyα forest while simultaneously resolving MH self-shielding.Instead, we will employ two simulation methods at different scales to address this issue.In Section 2, we describe our 1-dimensional radiationhydrodynamics calculation for modeling the HI column density of individual MHs after exposure to the ionizing background radiation.In Section 3, we outline our largescale simulation for modeling the intergalactic medium and explain how we account for the small-scale absorption by MHs based on the results from Section 2. In Section 4, we present our main results on the impact of MH absorption on the Lyα forest.Finally, in Section 5, we summarize our work.Throughout this work, we adopt the following cosmology parameters: h = 0.675, Ω M = 0.31, Ω M = 0.31, Ω b = 0.0487, σ 8 = 0.82, and n s = 0.965, which are consistent with measurements by the Planck satellite (Planck Collaboration et al. 2020).

Simulation of Minihalo Photo-evaporation
The process of photoevaporation in MHs involves radiative transfer and hydrodynamics at sub-kpc scales, which is generally not feasible to resolve in simulations of Lyα forest covering ∼ 100 Mpc.Therefore, we utilize the 1D radiation-hydrodynamics code, as fully described in Ahn & Shapiro (2007) and further developed by Park et al. (2016), to calculate the HI column density in individual MHs exposed to the ionizing background radiation.This 1D code has been extensively tested for various cases with analytic solutions (See Appendix C of Ahn & Shapiro (2007)).
In this work, we extrapolate the fitting formula for the truncated isothermal sphere (TIS) profile to ten times the truncation radius, r t , to have a reasonable description for the outskirt of the MH.We keep track of N sh = 10,000 radial shells linearly spaced from r = 10 −3 r t to 10 r t .We bound the outermost shell with the pressure of that shell at the initial time step, which becomes practically negligible as soon as the ionization of outer shells photoheats the gas above 10,000 K.This setup is identical to that in the appendix of Park et al. (2016).
For the initial conditions, we adopt the TIS profile, which represents the equilibrium state of MHs in the absence of star formation (Shapiro et al. 1999).The fitting formula for this TIS profile is provided in Appendix A of Shapiro et al. (1999).To describe the outskirts of MHs, we extend this TIS profile to 10 times the truncation radius r t using the power-law slope at r t .The mass inside r t is similar to the virial mass.In this calculation, we consider 2, 3, 5, 7, 10, 20, 30, and 50 million M ⊙ to cover a mass range that would survive the ionizing background for a significant amount of time (≳ 10 8 yr).The code tracks 10,000 radial shells between 10 −3 and 10 r t .While the profile shape depends on the redshift of collapse, z col , previous works have shown that the photoevaporation rate is relatively insensitive to this quantity (Shapiro et al. 2004;Nakatani et al. 2020).
For the ultraviolet background (UVB) radiation, we adopt the blackbody spectrum of T bb = 100,000 K.We truncate the UVB spectrum above 54.4 eV to account for absorption by HeII.The radiation is thought to con- tain starlight from O-and B-type stars with a temperature range from 30,000 to 50,000 K, as well as EUV and X-ray radiation from accreting neutron stars and quasars (e.g., Haardt & Madau 2012).However, the exact spectral shape of the radiation remains unknown.The photoevaporation of MHs depends on the hardness of the ionizing radiation background, as higher-energy photons result in a thicker ionizing layer on MHs due to their longer mean free path (Iliev et al. 2005b;Nakatani et al. 2020).The hardness of the spectrum we adopted falls between the hardness of these two components.
For simplicity, we assume a constant and isotropic ionizing radiation with J −21 = 0.01 and 0.1.Here, J −21 represents the angle-averaged intensity, I ν dΩ/(4π), at the Lyman limit in units of 10 −21 erg cm −2 s −1 Hz −1 sr −1 .Another commonly used quantity in the literature is where σ ν is the cross section of hydrogen atoms to ionizing photons.This quantity represents the ionization rate per H atom in units of 10 −12 s −1 .Our choices of ionizing intensity, J −21 = 0.01 and 0.1, correspond to Γ −12 = 0.03 and 0.3, respectively.We will refer to this value when referring to each case in this paper.
Measurements from the quasar proximity zone suggest a steep increase of Γ −12 near the end of reionization, with the value rising to Γ −12 ∼ 1 from a lower, yet unknown, value (see Gaikwad et al. (2023) for the latest measurement).As a typical value for reionization, we will consider Γ −12 = 0.03 based on the reionization model by Ocvirk et al. (2021, see their Fig. 4).
Figure 1 shows the time evolution of the radial HI column density profile for halo masses of 2 × 10 6 and 2 × 10 7 M ⊙ and ionizing intensity of Γ −12 = 0.03 and 0.3.In the initial conditions given by the TIS profile, the HI column density exceeds 2 × 10 20 cm −2 at the core for both masses.When exposed to the lower ionizing intensity Γ −12 = 0.03 of the reionization era, the HI column density remains nearly unchanged for both halo masses even after 300 Myr.However, with the higher intensity Γ −12 = 0.3 of the post-reionization era, the core of the 2 × 10 7 M ⊙ halo remains as a DLA even after 300 Myr, while the core of the 2 × 10 6 M ⊙ halo loses most of its HI gas after 100 Myr and no longer appear as a DLA.
The HI column density profiles show that MHs with several million solar masses can withstand the ionizing radiation during reionization and manifest as DLAs for ∼ 100 Myr after reionization.More massive MHs above 10 7 M ⊙ would survive the stronger ionizing background of the post-reionization Universe and remain as DLAs for much longer durations.These findings align with the results of Shapiro et al. (2004) and Nakatani et al. (2020) using similar simulations.In Section 3, we will utilize these HI density profiles to estimate the HI column density in FoF halos of the large-scale cosmological simulation.

Simulation of Lyα Forest
We employ Nyx as our simulation tool for modeling the evolution of the intergalactic medium.Nyx is a cosmological N -body/hydrodynamics code specifically designed for efficient execution on massively parallel computers (Almgren et al. 2013).It solves the hydrodynamic equations for gas using a grid-based approach, while dark matter is represented by collisionless particles.Nyx exclusively relies on the particle-mesh calculation for gravitational dynamics.This allows Nyx to focus on simulating mildly overdense structures that are relevant to the Lyα forest while saving computational time by excluding particle-particle gravity calculations for highly overdense structures where Lyα transmissivity is negligible.
To ensure convergence in the Lyα flux statistics, it is necessary to simulate a volume of ≳ 100 Mpc while resolving structures down to scales of ∼ 10 kpc (Lukić et al. 2015).We initialize the matter density field within an 80 h −1 Mpc box at z i = 199 using 4096 3 dark matter particles with 6.5 × 10 5 h −1 M ⊙ each.Once the simulation starts, the gas density and velocity are initialized on a 4096 3 mesh, assuming baryons trace dark matter.= 9, 7, 6, 5.5, 5, and 4  To model inhomogeneous reionization, we use a hybrid scheme that implements reionization physics into largescale cosmological hydrodynamics simulations based on a small set of phenomenological input parameters (see Oñorbe et al. 2017a).During the onset of reionization, each cell is photoheated by ∆T = 20, 000 K according to the input reionization redshift (z re ) field, following Oñorbe et al. (2019).After the cell has been exposed to reionization, the UVB model of Oñorbe et al. (2017a, "Middle HI Reionization") is assumed to calculate the temperature evolution.

Inhomogeneous Reionization
We calculate z re on a 128 3 grid according to the method presented by Battaglia et al. (2013), where z re is calculated by smoothing and rescaling the initial density field with a Gaussian filter of 1 h −1 Mpc.Following the approach of Trac et al. (2022), we rescale the smoothed density field to match a desired global reionization history, xe (z), while preserving the order among the cell values (i.e., higher z re for higher smoothed density).In this study, we adopt a simple analytic model with free parameters representing the beginning (z begin ), middle (z middle ), and end of reionization (z end ).The model is given by where X b = 0.5(tanh[A(z − z begin )] + 1), X e = 0.5(tanh[A(z − z end )] + 1), and A = 1.1.The reionization history from this toy model is shown in Figure 2. We choose z begin = 9.5, z middle = 7, and z end = 5.5 in this work.The choice of z end = 5.5 is consistent with the recent measurement of the Lyman limit mean free path at z ≥ 5.5, suggesting that reionization ended at z ≲ 6 (Becker et al. 2021).In future applications, these parameters will be varied to explore different reionization histories and study the dependence of the Lyα forest on these parameters.Also, the input reionization model can easily be replaced with a more sophisticated one such as 21CMFAST (Mesinger et al. 2011).To assess the impact of inhomogeneous reionization, we also conducted another simulation with the same initial conditions but with z re = 7 assigned everywhere in the simulation volume.
The gas temperature maps in Figure 3 provide an overview of the simulation and impact of the inhomogeneous reionization at large scales 1 .Between z = 9 and 5.5, reionization heats the IGM from 100 K or below to 20,000 K, starting from overdense regions and progress- ing toward underdense regions.By z = 5.5, reionization is complete, but there is still a mild large-scale variation in temperature due to regions that reionized later having less time to adiabatically cool.This reionization relic is visible at z = 5, but it becomes mostly diluted by z = 4.

Lyα Opacity of MHs
First, we calculate the Lyα opacity in redshift space from the output HI density/velocity/temperature field of the Nyx simulation, assuming that the z-direction is the line-of-sight direction.However, this calculation does not include the opacity from self-shielded MHs since Nyx assumes that the cells are optically thin to the ionizing background.Therefore, we assign the HI column density to the MHs based on the 1D simulation results described in Section 2.
Next, we identify MHs in the snapshots of the Nyx simulation at z = 4.5, 5, and 5.5 using the Friends-of-Friends (FoF) algorithm implemented in the Nbodykit package2 .We set the linking length to 0.2 times the mean particle separation.Among the FoF halos identified, we consider those with less than 5 × 10 7 M ⊙ as MHs.We identify FoF groups down to 6 particles or 4 × 10 6 M ⊙ in mass.It should be noted that the accuracy of halo findings is not guaranteed at such small numbers of particles, but the identified particle groups give a reasonable estimate for the locations of MHs.We compare the mass function of the FoF halos at z = 5.5 to the Sheth-Tormen mass function in Figure 4.Although the two cases show a reasonable agreement, our simulation underproduces halos at around 10 7 M ⊙ by up to a factor of 3.This mass is in the middle of MH mass range considered in this work (4 × 10 6 − 5 × 10 7 M ⊙ ).While we do not aim for a precise recovery of the Sheth-Tormen function given that the mass function (MF) is poorly constrained at such low mass and high redshift, we shall consider a possible underestimation of MH absorption due to this difference when analyzing our results.
For every identified MH, we calculate the exposure time to the ionizing radiation based on the reionization redshift field.We assume the MH had a constant mass since its exposure to the radiation.We also assume a constant UVB intensity during exposure.Using the halo mass and exposure time of each identified MH, we calculate the HI column density profile based on the 1D simulation results described in Section 2. This allows us to determine the HI column densities of intervening MHs along any given sight line and for two ionizing intensities: Γ −12 = 0.03 and 0.3.When calculating the Lyα opacity due to the MHs, we treat them as point-like absorbers with a temperature of 10,000 K.This choice is made considering that Lyα opacity is insensitive to these parameters in the damping-wing regime.By doing so, we properly account for the absorption by MHs in the Lyα flux data obtained from our simulation.The UVB in our Lyα forest simulation evolves over time, as should be the case for the real Universe (see Fig. 15 of Oñorbe et al. 2017a).Thus, the constant Γ −12 assumed in our MH opacity calculation needs to be close to the time average of the evolving UVB over the photoevaporation time scale to yield realistic results.Reionization occurs between z = 9 and 5.5 in our model, spanning nearly 500 Myr (Fig. 2), and photoevaporation can also take several hundred Myr when the UVB is weak (Γ −12 ∼ 0.03) or the minihalo is on the massive end (M h > 10 7 M ⊙ ), as shown in Figure 1.At z = 5.5, the time-averaged Γ −12 can be much lower than the observed post-reionization value (Γ −12 ≳ 0.3), as photoevaporation of MHs may have taken place under a much weaker UVB of the reionizing Universe for most of the time, or it can also be as high as the post-reionization value if the UVB intensity does not fall steeply toward high-z.Thus, we mainly consider the Γ −12 = 0.03 case at z = 5.5, but also consider the Γ −12 = 0.3 case to explore the dependence of the Lyα forest on UVB intensity.For z = 5 and 4.5, Γ −12 = 0.3 should be a realistic choice as they are well into the post-reionization era.
We also note that the Lyα forest is insensitive to the UVB model if the mean opacity is rescaled to a certain target value as demonstrated in Lukić et al. 2015 (see their Sec.7).Since the mean flux will be constrained  1. Volume-averaged normalized Lyα flux for no-MH cases and MH-included cases with Γ−12 = 0.03 and 0.3 calculated for z = 5.5, 5, and 4.5.For the MH-included cases, we also provide the fractional flux decrement with respect to the no-MH case of the corresponding redshift in parentheses next to the transmission value.We have boldfaced the results for the realistic combinations of Γ−12 and z.
through observations in practice, the actual value of the UVB intensity used in the Lyα forest simulation is unimportant.

Impact of Minihalo Absorption on Lyα Forest
In Table 1, we tabulate the volume average of the normalized Lyα flux from our simulation for z = 5.5, 5, and 4.5 and for three cases of MH absorption: absorption not included, and absorption included with Γ −12 = 0.03 and with Γ −12 = 0.3.We find that MH opacity decreases mean flux by 2−3% when Γ −12 = 0.03 and by 0.1−0.5% when Γ −12 = 0.3.The impact of MH absorption subsides toward low redshift as their HI gas photoevaporates over time.The absorption may be a factor of few stronger if we consider the underproduction of MHs in our simulation.In any case, the MH absorption effect on the mean opacity would be negligibly small for z = 5 and 4.5, where Γ −12 ≳ 0.3.
In Figure 6, we present the Lyα flux along an example sight line in the z = 5.5 snapshot, demonstrating a notable absorption caused by self-shielded MHs assuming Γ −12 = 0.03.Panels (c) and (d) focus on the two neighboring MHs with 8.6 × 10 6 and 1.6 × 10 7 M ⊙ at x ≈ 16.2 h −1 Mpc, where several MHs cluster along a filamentary structure.These two MHs exhibit N HI = 10.0 × 10 20 and 4.1 × 10 20 cm −2 for this sight line and appear as a single DLA with N HI = 14.1 × 10 20 cm −2 in the spectrum due to their proximity.This example qualitatively highlights the impact of a self-shielded system on the Lyα forest.MHs typically form in overdense regions with n HI > 10 −8 cm −3 , where the Lyα flux is negligible regardless of the absorption by MHs.However, the absorption by the MHs can extend to ∼ 10 h −1 Mpc, significantly changing the average flux and large-scale shape of the Lyα forest.Panel (d) demonstrates the clustering of MHs within an overdense structure, suggesting that sight lines are more likely to encounter self-shielded systems near massive objects.Each aspect of MH absorption seen in this example will be investigated further in this section.

DLA Incidence Rate
The latest measurement of the DLA incidence rate comes from the SDSS BOSS survey.The survey finds the incidence rate gradually increases from dN/dX = 0.033 at z = 2 to 0.069 − 0.106 (68% limit) at z ∼ 4.5, where X ≡ (1 + z) 2 H 0 /H(z)dz is the absorption distance (Ho et al. 2021).These DLAs are considered to be associated with galactic disks hosted by atomically cooling halos in the post-reionization Universe.Since MHs are more numerous than galaxies, self-shielded MHs can significantly contribute to the DLA number density at higher redshifts toward the end of reionization.
We calculate the DLA incidence rate contributed by self-shielded MHs by summing the projected area of their DLA portion, illustrated as blue-filled circles in Figure 5.This result is shown in Figure 7 as a function of the minimum halo mass considered, describing the cumulative contribution of MHs toward the low-mass limit.To address the impact of the difference between the simulated MF and the Sheth-Tormen function, as shown in Figure 4, we also calculate the incidence rate, correcting for this difference, and show the result in the right panel of Figure 7.
When considering all the simulated MHs down to 4 × 10 6 h −1 M ⊙ , the additional incidence rate due to MHs before correction for the difference in MF (left panel of Fig. 7) is 0.3(0.05) at z = 5.5 for Γ −12 = 0.03(0.3).The incidence rate decreases to a small value (≲ 0.02) at z = 5 and 4.5 for Γ −12 = 0.3.Considering that the measurement by BOSS indicates dN/dX = 0.1 at most at z = 4.5, our results suggest a steep increase in the DLA incident rate toward z ∼ 5.5 in the case of Γ −12 = 0.03 or a mild increase in the case of Γ −12 = 0.3.Thus, finding a rise in dN/dX toward high-z can give a clue for the evolution history of the UVB near the end of reionization.
The right panel of Figure 7 shows the same quantities after correcting for the difference in MF in Nyx versus the Sheth-Torman profile.The incidence rate increased by a factor of ∼ 2 with similar trends, making the impact of MHs even more pronounced.
The incidence rate as a function of the MH mass shows the differential contribution from different halo masses.In the cases of z = 5 and 4.5, the curves flatten toward the low mass, below M h ∼ 1 × 10 7 M ⊙ .The flattening indicates that the MHs below the flattening mass are fully evaporated and do not contribute to the DLA population, and therefore, dN/dX has fully converged at this mass resolution.The flattening does not occur for Γ −12 = 0.03 and z = 5.5 (black solid line), suggesting that the photoevaporation of low-mass (∼ 10 6 M ⊙ ) MHs is still ongoing.This finding aligns with the 1D simulation result presented in Section 2 (Fig. 1), which indicates that the ionizing background with Γ −12 = 0.03 cannot completely photoevaporate MHs even after several hundred million years of exposure.
The results here show MHs can substantially contribute to the DLA population due to their large number.Considering that dN/dX has not fully converged at 4 × 10 6 h −1 M ⊙ , MHs below this mass may also con-tribute significantly to the Lyα opacity at z ∼ 5.5.On the contrary, the contribution of MHs appears minimal at z ≤ 5 with Γ −12 = 0.3.In this regime, neutral gas from larger halos, which we do not account for in this calculation, is likely the main constituent of DLAs.

Probability Distribution of Lyα Opacity
The absorption by MHs weakens the overall Lyα flux by introducing extra opacity of sight lines that encounter self-shielded MHs.To quantify this effect, we calculate the effective optical depth for approximately 26,000 equally spaced skewers in our simulation.The effective optical depth is defined as τ eff ≡ − log(⟨F⟩), where ⟨F⟩ represents the mean flux along 50 h −1 Mpc line segments.This quantity is robust and unaffected by the spectral resolution, enabling a reliable comparison be- tween different surveys or simulations.In Figure 8, we present the cumulative probability distribution (CDF) of τ eff with and without accounting for MH absorption.
We compare the CDF of τ eff with and without accounting for MH absorption assuming Γ −12 = 0.03 at z = 5.5 in Figure 8.The CDF rises more slowly up to τ eff = 2.5 when MH absorption is included, indicating that a significant fraction of line segments experienced a increase in τ eff up to 2.5 from a lower opacity.The impact of MH absorption becomes negligibly small for Γ −12 = 0.3.Thus, we do not show this case and lower-redshift cases (z = 5 and 4.5), where we expect Γ −12 ≳ 0.3.
We also calculate and present the flux CDF for the instantaneous reionization case without MH absorption to assess the impact of inhomogeneous reionization in Figure 8.These two cases exhibit agreement at the low flux limit, but the inhomogeneous reionization case converges to unity more slowly, indicating a larger number of high flux segments.This is attributed to a wider distribution of the IGM temperature in the inhomogeneous reionization case, resulting in a broader flux distribution.These findings are consistent with the findings of previous works (D 'Aloisio et al. 2018;Ishimoto et al. 2022).
In practical calculations, the Lyα optical depth is rescaled to match the observed mean flux.To examine if the impact of MHs on CDF remains in such cases, we rescale the flux of one of the two cases (with and without MHs) to match the mean flux and assess the impact on the flux CDF.For instance, the mean flux of the case without MH absorption is ∼ 3.2% higher than that with the absorption for Γ −12 = 0.03 cases (see Table 1).We globally increase the optical depth field of the MH-included case to align the mean flux with the no-MH case.After rescaling the flux, the flux CDFs of the two cases match, suggesting the effect of MHs cannot be distinguished from the flux CDF in practice, unlike the inhomogeneity of reionization.

Flux Power Spectrum
At z ≳ 4, temperature fluctuations arising from inhomogeneous reionization (as illustrated in Fig. 3) significantly impact the flux power spectrum at k ∼ 0.1 h −1 Mpc, enabling to constrain the details of reionization with high-z Lyα forest.In this regard, we evaluate the influence of MH absorption on power spectrum statistics to assess the importance of considering them in the analysis of high-z quasar spectra.While at z ≲ 4 DLAs can be readily identified and subtracted due to their rarity and higher mean flux of the Lyα forest, the same is difficult at z ≳ 5 where the Lyα opacity is saturated in most of the cells.As a result, MH absorption features are likely to affect power spectrum analysis at those redshifts.
We begin by rescaling the optical depth of the MH absorption-included cases to align the mean flux with that of the no-MH absorption case.Subsequently, we compute the 1D flux power spectrum for 512 2 equally spaced skewers and present the median spectra at z = 5.5 in Figure 9. Additionally, we plot the case of instant reionization without MH absorption and fractional dif- ference to facilitate a comparison between the effects of inhomogeneous reionization and MH absorption.
In the case of MH absorption with Γ −12 = 0.03, the flux power is constantly lower than in the no-MHabsorption case by 2% at k ≳ 0.3 h −1 Mpc, but it rises toward low-k from 0.3 h −1 Mpc becoming 3% higher at 0.1 h −1 Mpc (blue solid line in the lower panel).The extended absorption by MHs, demonstrated in Figure 6, adds to the large-scale power, but instead reduces the small-scale power by removing the forest in the absorbed line segment.The impact of MH absorption is negligible for Γ −12 = 0.3 (blue dashed line), and it stays so at the other redshifts considered in this work.Thus, we do not show plots for those cases.
The wavenumber range affected by the MH absorption coincides with that affected by inhomogeneous reionization.In the comparison of the inhomogeneous reionization case (black solid) to the instant reionization case (red solid) in Figure 9, we observe that the inhomogeneity of reionization increases the flux power by 50% at k ∼ 0.1 h −1 Mpc at z = 5.5.This is fairly larger than the −2 to +3% modulation caused by MH absorption with Γ −12 = 0.03 at the same redshift.Recalling that the MH effect on the DLA incidence rate is ∼ 2 times stronger when corrected for the underproduction of MHs in our simulation, we can speculate that the impact on the flux power spectrum could be as large as ∼ 5%.We also note that the instant reionization scenario is an extreme case, which likely gives a maximal difference from our inhomogeneous reionization scenario.Therefore, the MH absorption can be a subdominant but nonnegligible effect, which can introduce a bias when constraining the reionization history parameters such as the duration and midpoint with the Lyα forest.Given the negligible impact of MH absorption with Γ −12 = 0.3, this complication can be avoided by utilizing the lower-redshift (z ≲ 5) Lyα forest, where the MH effect should have subsided.

Cross correlation with Galaxies
Park et al.
In this section, we explore whether absorption by MHs can be observed in the galaxy-Lyα correlation signals.
As the sample size of the high-z galaxies and quasar sight lines continues to increase, we anticipate a significant improvement in this measurement in the coming decades.
As depicted in panel (d) of Figure 6, the MHs are highly clustered in overdense structures, implying that the absorption due to MHs would also correlate with massive structures.In Figure 10, we compare the flux field around the most massive FoF halo with that at a random location along a plane parallel to the LOS direction in the z = 5.5 snapshot for the Γ −12 = 0.03 case.The flux map reveals a higher density of MH absorption near the halo compared to at the random location.The Lyα flux around the halo is already lower than the average without the MH absorption up to a few Mpc from the halo, but the extended Lyα absorption by DLAs extends beyond that local trend along the LOS direction up to ∼ 10 h −1 Mpc from the halo.
To quantify the spatial anticorrelation between flux and galaxies, we stack the Lyα flux around 2000 largest FoF halos on a grid of LOS separation (r ∥ ≡ ∆r z ) and perpendicular-to-LOS separation (r ⊥ ≡ ∆r 2 x + ∆r 2 y ; i.e., impact parameter).The majority of the FoF halos have masses around 10 11 M ⊙ , corresponding to M UV between −19 and −21.The results are shown in Figure 11 for Γ −12 = 0.03 and z = 5.5.We do not show the result for the Γ −12 = 0.3 case as the impact is negligible and, therefore, for lower redshifts where the ionizing background should be much stronger (Γ −12 ≳ 0.3).
The flux contours resemble the typical galaxy twopoint correlation function as both galaxy and flux trace the underlying density.The contours are globally compressed along the LOS direction due to the linear gravitational motion (a.k.a. the Kaiser effect), but they are locally stretched along the LOS direction for perpendicular separation smaller than 0.5 h −1 Mpc, extending up to r ∥ ∼ 3 h −1 Mpc or ∼ 300 km s −1 in LOS velocity.This stretching is caused by the strong nonlinear peculiar motion of matter near the halo, also referred to as the finger-of-god (FoG) effect.
Comparing the solid line to the dotted line illustrates the impact of MH absorption.MH absorption stretches the contours further, similar to the FoG effect, with the effects extending to larger LOS separations beyond r ∥ ∼ 3 h −1 Mpc, reaching 10 h −1 Mpc.
As in the case of the DLA incidence rate, the Lyαgalaxy correlation signal can also be ∼ 2 times stronger if we consider the underproduction of MHs in our simulation.In that scenario, the signal can be detected provided a sufficient amount of data is collected for small impact parameters (r ⊥ < 1 h −1 Mpc) at high enough redshift z ≳ 5.5.Additionally, the sensitivity of this cross-correlation signal to UVB intensity can be exploited to constrain Γ −12 from the correlation signal.
The most relevant observation comes from Meyer et al. (2020), where the authors measured the correlation for 21 LAEs and 13 LBGs at z ∼ 6 in the proximity of eight quasar sight lines from z ≳ 6.The data within 10 cMpc is not sufficient to detect the signature of MHs described in Figure 11.However, they do find a significant decrease in the flux toward the galaxy at r < 10 cMpc, which broadly aligns with our results.With an increasing sample size, the cross-correlation signal would also be detectable in the future.

Limitations
Our method, as described in Section 3.2, unavoidably relies on several simplifying assumptions to capture subkiloparsec-scale scale physics in a volume of ∼ 100 Mpc.While most of these assumptions are not expected to significantly alter our qualitative findings, further investigations are necessary to assess how these assumptions may affect the results.
First, we assume the identified MHs to have the same mass since exposure to ionizing radiation.This assumption could lead to an overestimation of the MH absorption if the MH absorbing the quasar light were exposed to reionization much earlier.The time difference between the redshifts that we calculate in the Lyα forest (z = 5.5, 5, and 4.5) and the reionization redshift in our model (z re = 5.5 − 9) can span several hundreds of megayears, which is longer than the typical growth time-scale for MHs (e.g., see the discussion in Sec.5.6 of Nakatani et al. 2020).In such cases, the remaining neutral gas in MHs would be overestimated.Our results in Section 4.4 are most likely to be affected by this assumption as massive halos are typically located in overdense regions where reionization happens earlier than in other regions.We find that the surroundings of the 2000 massive halos were ionized before z = 8 in most cases when less than 20% of the simulation volume was ionized.For the same reason, the calculated MHs absorption at z = 4.5 is likely exaggerated, which would corroborate our finding that the MHs absorption is negligible at that time.
Our calculations do not take into account any star formation inside MHs.This is considered reasonable during reionization, as the MHs are sterilized by the Lyman Werner (LW) background emitted by star-forming galaxies.However, it is possible that PopIII stars from z ≳ 20 affected the MHs during reionization via feedback.The metal enrichment resulting from the deaths We have not accounted for the influence of the X-ray background on MHs.Hard X-rays have the ability to penetrate into the core of MHs and partially ionize the gas.This could lead to two potential effects: (a) the MH core might expand in response to increased pressure, or (b) it could contract by promoting H2 formation via the free electrons followed by H2 cooling3 .A separate investigation is required to quantify this effect using methods similar to our MH evaporation simulation.
Our simulation did not take the inhomogeneity in the ionizing background intensity into account.The back-ground radiation is expected to be stronger at overdensities, where galaxies are clustered.The galaxy-Lyα correlation measurement by Meyer et al. (2020) shows that the Lyα flux is above the mean at separations between 10 and 20 cMpc because the IGM density is close to the cosmic mean, while the ionizing intensity is above the mean at those separations.We do not find this trend in our calculation because the Lyα flux is calculated assuming a uniform ionizing background intensity.On the one hand, a higher mean Lyα flux at overdensities could lead to a stronger impact from MH absorption clustered at these locations.On the other hand, the stronger UVB background of overdense regions may result in earlier photoevaporation of MHs, thereby reducing their impact.
Lastly, as mentioned in Section 3.2, the constant UVB intensity assumed for the MH photoevaporation simulation is another simplification made during our calculation.The global UVB intensity is expected to rise steeply toward the end of reionization.Additionally, the intensity can locally evolve in the vicinity of massive galaxies according to their star formation rates.There-Park et al. fore, the MH photoevaporation process can differ from our results, even if the fixed Γ −12 matches the time average of the global UVB intensity.

Summary and Conclusion
The Lyα forest at z ≳ 5 is emerging as an effective probe of cosmology and reionization.During these high-redshift epochs, it is expected that the intergalactic medium contains MHs that formed prior to reionization and have yet to complete their photo-evaporation.Despite their potential significance, the impact of selfshielded MHs on the Lyα forest has not been systematically explored in the literature.This is partially because studying them is challenging due to their small size and the complex interplay of physical processes involved in their formation and evolution.This work addresses the gap using a hybrid scheme incorporating the HI column density of MHs, obtained from small-scale 1D simulations, into the cosmological simulation from Nyx to calculate the opacity arising from the MHs in the Lyα forest.We furthermore include the effect of inhomogeneous reionization in the IGM using a simple parametric model based on the ansatz from Battaglia et al. (2013) and Trac et al. (2022).
Our results are based on several simplifying assumptions that need to be tested in future studies.Ideally, one would directly resolve the MH evaporation in cosmological simulations rather than relying on those assumptions.Zoom-in or Lagrangian techniques may be required to cover the large dynamic range required for such simulations.
The impacts of MHs in our study are summarized as follows.
• The incidence rate of DLAs, dN/dX, increases steeply toward high-z from z ∼ 4.5 to 5.5 by ≳ 0.1.
• The Lyα flux is decreased by up to ∼ 5%.
• Flux around massive halos at short perpendicularto-LOS separation (r ⊥ ≲ 2 h −1 cMpc) can be particularly more suppressed due to clustered MHs in dense environments.
• These impacts are pronounced when assuming Γ −12 = 0.03 for MH photoevaporation (possible near the end of reionization) but diminish to negligible levels for Γ −12 = 0.3, anticipated at lower redshifts.
In conclusion, MHs can significantly influence the statistics of the high-z Lyα forest in the period shortly after the end of reionization.However, this effect is expected to diminish within ∼ 10 8 yr as the photoevaporation of HI gas in MHs progresses, particularly from the low-mass end, with the rapid rise of the ionizing background intensity after the end of reionization.Consequently, it is advisable to account for these effects when analyzing quasar spectra at z ≳ 5.5, whereas they are less impactful at z ≲ 5. Conversely, measuring the signature of MH absorption at z ≳ 5.5 would provide insights into the ionization status of MHs and the UVB intensity during that era.

Acknowledgement
The authors thank K. Kakiichi, N. Yoshida and the anonymous referee for the helpful comments on this work.This

Figure 2 .
Figure 2. Volume-averaged global reionization history of the inhomogeneous reionization model used in this work.

Figure 3 .
Figure 3. Temperature map of a slice from the Nyx simulation with the inhomogeneous reionization model.The six panels depict the snapshots at z = 9, 7, 6, 5.5, 5, and 4 with the color representing the temperature.Reionization takes place between z = 9 and 5.5 in the model used in this work.
Figure 3. Temperature map of a slice from the Nyx simulation with the inhomogeneous reionization model.The six panels depict the snapshots at z = 9, 7, 6, 5.5, 5, and 4 with the color representing the temperature.Reionization takes place between z = 9 and 5.5 in the model used in this work.

Figure 4 .
Figure 4. Cumulative halo mass function from the z = 5.5 snapshot of our Nyx simulation run (black) compared to that from the Sheth-Tormen function (blue).
Figure 5 displays the projected distribution of MHs on the xy plane, covering a region of 0.35×0.35(h −1 Mpc) 2 plane along the line of sight throughout the simulation box (80 h −1 Mpc).The size of the circle indicates the dense core where damped Lyα absorption occurs (i.e., N HI > 2 × 10 20 cm −2).If a sight line intersects with one of these circles, the corresponding MH will appear as a DLA in the Lyα forest.The DLA cores have sizes of ∼ 1 h −1 kpc, which is smaller than the cell size of our simulation (20 h −1 kpc), described as the grid in the figure.When calculating the Lyα flux on the mesh data with MH absorption, we assume that all the sight lines are passing through the center of each cell.

Figure 6 .
Figure 6.Panel (a): Lyα flux of an example sight line in the z = 5.5 snapshot of the simulated IGM.We compare the flux without MH absorption (black dotted) vs. with MH absorption assuming Γ−12 = 0.03 (black solid).The red solid line shows the absorption by the MHs alone.The red vertical lines mark the location of MHs intersecting with the sight line.Panel (b): HI density map of a slice with the sight line shown as a black straight line.Panel (c): the zoom-in HI density map of the squared region in panel (b).Panel (d): the zoom-in HI density map of the squared region in panel (c).We show the MHs on the slice as blue circles.The lighter blue circle marks the virial radius of the MHs, and the dark ones mark the DLAs.The color scale given in panel (c) also applies to panels (b) and (d).

Figure 7 .
Figure 7. Incidence rate of MH DLAs per absorption distance as a function of the minimum MH mass considered.The black, blue, and red lines describe the z = 5.5, 5, and 4.5 cases, respectively, and the solid and dashed lines describe the cases in which MH absorption is included assuming Γ−12 = 0.03 and 0.3, respectively.In the left panel, the results are directly from the simulated MHs without any correction, while in the right panel, the results are corrected for the difference between the simulated MF and the Sheth-Tormen MF.

Figure 8 .
Figure 8. Cumulative probability distribution function of τ eff ≡ − log ⟨F ⟩ for z = 5.5.The black and red solid lines describe the case that we do not account for the MH absorption in the inhomogeneous and instantaneous reionization models, respectively.The blue solid line describes the case in which MH absorption is introduced to the inhomogeneous reionization model assuming Γ−12 = 0.03.The cyan dashed line from rescaling the MH-included case to match the mean opacity with the no-MH case.

Figure 9 .
Figure 9. Upper panel: dimensionless power spectra of 1D flux, kP (k)/π, for inhomogeneous reionization run without MH absorption (black) with MH absorption assuming Γ−12 = 0.03 (blue) and instant (zre = 7) reionization case without MH absorption (red) at z = 5.5.Lower panel: fractional difference to the inhomogeneous reionization with no-MH absorption case.The line colors correspond to the same cases as in the upper panels, and we additionally show the result for the Γ−12 = 0.3 case with a blue dashed line.

Figure 10 .
Figure 10.Lyα flux map around the most massive FoF halo with M h = 4 × 10 12 M⊙(upper panel) vs. a random location in the simulation (lower panel) at z = 5.5.White regions are where the flux is significantly transmitted (30% or above), and gray regions are where the flux is completely absorbed by the IGM without MHs.The blue stripes are where the flux is absorbed by intervening MHs in the Γ−12 = 0.03 case.The LOS direction is along the x-axis, and the y-axis points to a direction perpendicular to the LOS.The red circle at the center of the upper panel marks the virial radius of the halo.

Figure 11 .
Figure 11.Stacked Lyα flux around 2000 most massive FoF halos shown with line contours for the z = 5.5.The x and y axes are perpendicular and parallel to the LOS direction.The solid and dotted contours describe no-MH and MH-included cases with Γ−12 = 0.03, respectively.Each case has nine line contours marking 10, 20, ..., 90% of the mean Lyα flux from inside to outside.
work was partially supported by the DOE's Office of Advanced Scientific Computing Research and Office of High Energy Physics through the Scientific Discovery through Advanced Computing (SciDAC) program.The development of Nyx as an AMReX application was supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics program under contract number DE-AC02005CH11231, and by the Exas-cale Computing Project (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy Office of Science and the National Nuclear Security Administration.This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DEC02-05CH11231.This research also used resources from the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357.
Park et al.