The impact and response of mini-halos and the inter-halo medium on cosmic reionization

An ionization front (I-front) that propagates through an inhomogeneous medium is slowed down by self-shielding and recombinations. We perform cosmological radiation hydrodynamics simulations of the I-front propagation during the epoch of cosmic reionization. The simulations resolve gas in mini-halos (halo mass 10 4 ≲ M h [M ⊙ ] ≲ 10 8 ) that could dominate recombinations, in a computational volume that is large enough to sample the abundance of such halos. The numerical resolution is sufficient (gas-particle mass ∼ 20M ⊙ , spatial resolution < 0 . 1 ckpc ) to allow accurate modelling of the hydrodynamic response of gas to photo-heating. We quantify the photo-evaporation time of mini-halos as a function of M h and its dependence on the photo-ionization rate, Γ − 12 , and the redshift of reionization, z i . The recombination rate can be enhanced over that of a uniform medium by a factor ∼ 10 − 20 early on. The peak value increases with Γ − 12 and decreases with z i , due to the enhanced contribution from mini-halos. The clumping factor, c r , decreases to a factor of a few at ∼ 100 Myr after the passage of the I-front when the mini-halos have been photo-evaporated; this asymptotic value depends only weakly on Γ − 12 . Recombinations increase the required number of photons per baryon to reionize the Universe by 20-100 per cent, with the higher value occurring when Γ − 12 is high and z i is low. We complement the numerical simulations with simple analytical models for the evaporation rate and the inverse Strömgren layer. The study also demonstrates the proficiency and potential of sph-m1rt to address astrophysical problems in high-resolution cosmological simulations.


INTRODUCTION
The epoch of reionization (eor) refers to the time during which luminous sources (e.g.hot stars) transformed the inter-halo medium (hereafter ihm) 1 from mostly neutral to mostly ionized (see reviews by, e.g.Loeb & Barkana 2001;Madau 2017or Wise 2019).Following cosmic recombination at redshift z ∼ 1100, the Universe remained neutral up to redshifts z ∼ 15 − 30 when the first stars formed and emitted ionizing photons (Tegmark et al. 1997;Abel et al. 2002;Bromm et al. 2002).These stars and the first galaxies ionized their immediate surroundings.With more and brighter galaxies forming, cumulatively more than one ionizing photon per hydrogen was emitted, causing these cosmological HII regions (Shapiro & Giroux 1987) to percolate (Gnedin 2000a;Barkana & Loeb 2004;Furlanetto et al. 2004;Iliev et al. 2006a), signalling the end of the eor.After reionization, the ionizing photons emitted by galaxies and quasars keep the ihm highly ionized, with ionizations approximately balanced by recombinations in higher density regions (Haardt & Madau 2012;Bolton & Haehnelt 2007).
Observations of the polarization of the cosmic microwave background (cmb) due to Thomson scattering place the midpoint of the eor (where about half the ihm is highly ionized) at around z ∼ 8 (Planck Collaboration et al. 2020).This is consistent with the evolution of the fraction of galaxies that are detected in Lyman-α emission (Mason et al. 2018), and constraints on the neutral fraction in the ihm inferred from the spectra of high-redshift quasars (Fan et al. 2006;Mortlock et al. 2011;Bouwens et al. 2015).Measurements of the patchy kinetic Sunyaev-Zeldovich effect in the cmb place limits on the duration of the eor (Zahn et al. 2012;George et al. 2015), with Planck Collaboration et al. (2016) setting an upper limit to the width of the reionization period of ∆z ≲ 2.8.
There are three main challenges to our understanding of the eor: (i) determining the rate at which sources emit ionizing photons, (ii) computing the escape fraction, fesc, of ionizing photons that contribute to ionizing the ihm (rather than being absorbed in the immediate surroundings of where they were produced), and (iii) understanding the nature and evolution of photon sinks, where ionized gas recombines again.
The precise nature of the dominant ionizing sources remains controversial.These could be stars in very faint galaxies below the detection limit of the Hubble Space Telescope deep fields (with 1500Å magnitude MUV ≳ −13; Robertson et al. 2013;Finkelstein et al. 2019).Another possibility is slightly brighter galaxies that have higher values of fesc because they drive strong winds (Sharma et al. 2016(Sharma et al. , 2017;;Naidu et al. 2020).jwst promises to measure the slope of the faint-end luminosity function, which could constrain the contribution of faint galaxies to the emissivity.However, it will remain challenging for observations to constrain fesc directly, although the UV-slope and the dominance of nebular lines in the spectra of galaxies may be good proxies for fesc (e.g.Chisholm et al. 2018Chisholm et al. , 2022)).
In this paper, we focus on Challenge (iii): understanding the nature and evolution of photon sinks.Clumpy gaseous structures impede reionization by boosting the opacity (often described in terms of the mean free path of ionizing photons) and the recombination rate (which consumes ionizing photons) of the ihm.
The recombination rate in an inhomogeneous medium relative to that of a uniform medium is characterised by the 'clumping factor', c l,all ∼ ⟨n 2 HII ⟩/⟨nH⟩ 2 , where ⟨⟩ denotes a volume average and nHII and nH are the density of ionized and total gas.Calculating the clumping factor accurately is challenging.Its numerical value is uncertain, ranging from 2 to 30 (e.g.Gnedin & Ostriker 1997;Trac & Cen 2007;Iliev et al. 2007;Pawlik et al. 2009;McQuinn et al. 2011;Raičević & Theuns 2011;Finlator et al. 2012;Shull et al. 2012), depending on its definition and the simulation outcome.
To further exacerbate the issue, a significant number of recombinations occurs in gas inside the smallest gravitationally bound structures: mini-halos.These are low-mass halos, 104 ≲ M h /M⊙ ≲ 10 8 , not massive enough to form a galaxy (without metals or molecules) but able to retain their cosmic share of baryons before reionization.Haiman et al. (2001) claimed that mini-halos boost the number of photons that are required to reionize the universe by a factor of ∼ 10 (although they assumed mini-halos are optically thin).The semi-analytical model of Iliev et al. (2005b) and the sub-grid model of Ciardi et al. (2006) also suggest that mini-halos can boost reionization photon budget, based on high-resolution 2D radiation hydrodynamical (rhd) simulations (Shapiro et al. 2004).Emberson et al. (2013) performed a suite of high-resolution cosmological hydrodynamics simulations that resolve all mini-halos with at least 100 particles.Post-processing these simulations with radiative transfer (hereafter rt), they obtained values of c l,all ∼ 10, confirming the major impact of small-scale structures on reionization.But Emberson et al. (2013) might have overestimated the clumping factor by not including photoheating and photo-evaporation.Park et al. (2016) performed cosmological rhd simulations with photo-evaporation.These calculations resolve mini-halos, and model the ionizing background as a roughly uniform radiation field with Gadget-RT.But the simulated volumes are relatively small and may not capture the full mass range of mini-halos.D' Aloisio et al. (2020) performed rhd simulations in a simulation volume with linear extent ∼ 1 cMpc, using a uniform mesh with cell size ∼ 1 ckpc (which might be larger than the virial radius of low-mass mini-halos).Both Park et al. (2016) andD'Aloisio et al. (2020) concluded that c l,all ∼ 10 during the early stages (≪ 100Myr) of reionization, but then drops rapidly to lower values ∼ 2, due to photo-evaporation.But these studies do not explicitly study mini-halo photo-evaporation and its impact on reionization, which we investigate here.
Improving upon previous studies, we perform highresolution (dark matter particle mass ∼ 100 M⊙ and adaptive spatial resolution < 0.1 ckpc) cosmological rhd simulations in volumes of linear extent ∼ 1 cMpc to investigate the response of gas in small-scale structures on the passage of an ionization front and its impact on reionization.
The simulations are performed with the cosmological simulation code swift (Schaller et al. 2018;Schaller et al. 2023) 3 with smoothed particle hydrodynamics (sph ;Lucy 1977;Gingold & Monaghan 1977).We model radiation hydrodynamics with sph-m1rt, a novel sph two-moment method with local Eddington tensor closure (Chan et al. 2021).The simulated volume is large enough to sample the full mass range of mini-halos (see §3.2) and the resolution is high enough to resolve the small mini-halos 4 .
With this simulation suite, we investigate the minihalo photo-evaporation process and its dependence on redshift and photo-ionization rate.We study how high-density gas can stall the ionization front due to self-shielding.Secondly, we perform a detailed analysis of how the inhomogeneous medium impedes reionization, including the clumping factor, the number of recombination photons, and the role of mini-halos.We extend upon Chan et al. (2023), which presented only the clumping factor evolution of one of our representative simulations.
This paper is organized as follows.In §2, we describe the stages of reionization.We introduce our simulation suite in §3.The analysis of the simulations is described in §4, focusing on the impact of an ionization front on gas in cosmic filaments and mini-halos.We use these results to gauge the impact of these structures on the progress of reionization (in terms of recombinations and the clumping factor).We put our results in context and discuss them in §5.Finally, we conclude with a summary and an outlook for future research.A set of appendices provide the details of our numerical implementation and convergence tests.This is a long paper, and we recommend the reader to start by skimming through the following sections.First, the definition and discussion around mini-halos in §2.1 (Fig. 1) and the stages of reionization in the beginning of §2.2 (Fig. 3).The main results are Figs.13 and 14 in §4.4,which quantify the rate at which mini-halos photo-evaporate.Finally, the impact of small-scale structures on reionization is presented in §4.5 (Figs. 15,16,and 17).

Halos and their role in reionization
The impact of halos on the ihm pre-and post-reionization depends on their mass.Below some minimum mass, halos have too shallow potential wells to accrete or retain cosmic gas.The value of this minimum mass depends on the mean temperature of the gas, T ihm .Prior to reionization, T ihm depends on redshift z as (e.g.Furlanetto et al. 2006) . (1) Here, z dec is the redshift where the cosmic gas decouples from the cmb temperature.We consider the following three estimates of the (preionization) minimum mass: (i) the Jeans mass, MJ, MJ = 5kBT ihm GµmH (e.g.Mo et al. 2010), where µ ≈ 1.22 is the mean molecular weight per hydrogen atom of fully neutral primordial gas.We note that MJ decreases with cosmic time, therefore Properties of halos with different masses as a function of redshift, the yellow vertical band indicates the eor, assumed to occur around z ∼ 6. M J is the Jeans mass (black solid line, Eq. 2), M F is the filtering mass (black dashed line, Eq. 3), M min is the mass above which halos contain 50% of the cosmic baryon fraction on average as computed from the simulations (red dot-dashed line), M ACH is the halo mass above which gas at the virial temperature can cool atomically (blue dashed line, Eq. 4).Halos with mass M min < M h < M ACH are mini-halos (left magenta region).
After reionization, M char is the mass above which halos contain 50% of the cosmic baryon fraction on average as computed by low-mass halos need to accrete baryons to reach the cosmic baryon fraction.
(ii) the 'filtering mass', MF, introduced by Gnedin 2000b to account for the redshift dependence of T ihm , (3) see Appendix A for more details (iii) Mmin, the halo mass above which halos contain on average more than 50% of the cosmic baryon fraction, as determined from the simulations described in §3.25 .
The ihm temperature increases by almost three orders of magnitude during reionization, and hence so does the minimum mass.We use the simulations of Okamoto et al. (2008) to compute the (post-reionization) minimum mass and refer to it below as M char (the characteristic value of the postreionization minimum mass).We further define MACH as Figure 1 illustrates the evolution of these masses, with the vertical yellow band indicating the assumed eor.Pre-reionization, MJ and MF decrease with cosmic time, whereas Mmin remains approximately constant.Numerically, Mmin ≈ 2 × 10 4 M⊙, which is a factor of a few larger than MJ and a factor of a few smaller than MF .
Average baryon fractions of halos (pre-reionization) in units of the cosmic mean are shown in Fig. 2 at various redshifts; the horizontal dotted line indicates 50 per cent.The figure also illustrates that neither MJ nor MF provides accurate estimates of Mmin.
In Figure 1, the left magenta region represents minihalos, the main topic of this paper.They are halos with Mmin ≲ M h ≲ MACH, which contain their cosmic share of baryons but this gas cannot cool atomically6 .After reionization, M char > MACH implies that there are no more halos that contain gas which cannot cool atomically: all mini-halos are evaporated during reionization.Such gas-free halos occupy the grey region.In the lower part of the diagram, the ihm is highly ionized (black region) except in halos with mass above a characteristic mass, M char .The outskirts of the gas in these halos are highly ionized but the halos' potential wells are deep enough that the gas cannot photo-evaporate: these are Lyman-limit systems.The inner parts of the halos can selfshield, and the gas is neutral; these are damped Lyman-α systems (yellow, solid discs).In the middle part of the diagram, the ionization front has passed only recently, and the gas in mini-halos is still photo-evaporating (hollow circles).The ionization front has not yet reached the top of the diagram, where the ihm is still neutral (yellow region).Yellow vertical stripes are regions where gas remains neutral due to the shadowing of the ionizing radiation by an intervening self-shielded region.The I-front is slowing down and becomes increasingly corrugated rather than planar as a result of the density inhomogeneities in the ihm.
Finally, the right magenta region represents the 'Reionization-Limited HI clouds' (RElHIcs) in the mass range Llambay et al. (2017).They contain a significant amount of gas yet do not host a galaxy.
The main point to take away from Fig. 1 is this: before reionization, to capture all photon sinks requires resolving halos down to masses ∼ 10 4 M⊙7 .After reionization, halos with masses ≳ 10 8 M⊙ are the dominant photon sinks.This highlights the extended range in masses (from ∼ 10 4 M⊙ to ≳ 10 8 M⊙) of photon sinks during and after eor.

The progression of the reionization process
We divide the eor into three characteristic evolutionary phases, (i) the ionization front (I-front) propagation phase, (ii) the mini-halo evaporation phase, and (iii) the relaxed phase.These phases are illustrated in Fig. 3 and described below: • (I-front) Propagation Phase : the R-type8 I-front propagates at a small fraction of the speed of light, ionizing and heating the low-density ihm.Halos more massive than a minimum mass, Mmin, contain a large fraction of their cosmic complement of baryons, and this gas is dense enough to self-shield and hence initially remains neutral.Dense gas in such halos significantly increases the recombination rate compared to that of a uniform ihm.At the end of this phase, the neutral, high-density, self-shielded gas in halos is embedded in a highly ionized ihm.
• (Mini-halo) Evaporation Phase: the I-front propagates into halos, ionizing the gas in their outskirts.In halos more massive than a minimum mass, M char , most of this ionized gas does not photo-evaporate, and their inner parts remain self-shielded.In contrast, the gas in lower mass halos photo-evaporates at a speed slower than the sound speed, taking 10-100 of Myrs to complete9 .Such mini-halos are very numerous, and their photo-evaporation consumes a large number of ionizing photons, comparable to the number of photons needed to reionize the ihm (e.g.Shapiro et al. 2004;Iliev et al. 2005a).At the end of this phase, all mini-halos are photo-evaporated.Only halos more massive than M char contain significant amounts of gas.
• (Post-reionization) Relaxed Phase: The I-front also propagates into halos more massive than mini-halos, but their deeper potential wells prevent them from photoevaporating.These photon sinks correspond to Lyman-limit (the ionized outskirts of the halos) and damped-Lyman-α systems (the self-shielded inner parts, hereafter referred to as lls's and dla's, see, e.g.Theuns 2021).After reionization, these structures are mainly responsible for determining the opacity, the recombination rate and the clumping factor.

(I-front) Propagation phase
Consider the case of an initially planar I-front propagating through a (nearly) uniform ihm which is fully neutral at time t = 0.The speed of the I-front along its propagation direction10 , which we assume to be the Z-axis, is where ⟨nH⟩ is the (mean) hydrogen density by number, x is the neutral fraction at coordinate Z of the gas once it is ionized, and αB is the temperature-dependent recombination coefficient.The solution to this differential equation is provided we assume that the gas downstream of the I-front is highly ionized, x ≪ 1.Here, where tr is the recombination time and ZS is the onedimensional analogue of the Strömgren radius 11 .The photoionization rate, Γ, is related to F by and the case-B recombination rate is The numerical values in Eq. ( 8) assume that the ionizing spectrum is that of a black body with effective temperature TBB = 10 5 K, e.g.massive stars, which yields a frequencyaveraged photo-ionization cross section 12 of σHI = 1.63 × 10 −18 cm 2 .The temperature of the ihm is T ihm ≈ 2.2×10 4 K after being flash ionized by such an ionizing spectrum (e.g.Chan et al. 2021).
The I-front's propagation speed is considerably slower than the speed of light.Note that the location of the Strömgren layer, ZS, itself evolves as the Universe expands (see, e.g.Shapiro & Giroux 1987 for a more accurate calculation of the evolution of such Hii regions).Eq. ( 7) also shows that the recombination time at the mean density, tr, is longer than the age of the Universe, tU = 2/(3H(z)) = 0.65 (9/(1 + z)) 3/2 Gyr, at a redshift, z ≲ 8. Gas in minihalos is at a higher density and hence has a shorter recombination time, thus slowing down the I-front until the gas photo-evaporates.

(Mini-halo) Evaporation phase
The dense gas in mini-halos traps the I-front.The response of the gas, when overrun by the I-front, depends on ratios of three time-scales: (i) the sound-crossing time tsc, (ii) the I-front crossing time tic, and (iii) the recombination time t r,h .
The sound-crossing time is the time for a sound wave 11 Using the value of dZ/dt suggests that it would take of order ≳10 Myr for an I-front to cross 1 Mpc.In fact, it takes considerably longer than this, because of density inhomogeneities and recombination.
12 Our value differs from that of Emberson et al. (2013) because we use a different spectral shape.
to travel a distance equal to the virial radius, R h , of the halo, where M h is the virial mass of the halo and cs is the sound speed, which we evaluated at the temperature of the ihm after flash ionization (Chan et al. 2021).
The I-front propagating time, tic, is the time it takes for an I-front to cross a mini-halo, neglecting recombinations, Myr .
(11) Finally, the mean recombination time of the gas in a halo, t r,h , is: We evaluated the case-B recombination coefficient at a temperature, T = 2.2 × 10 4 K, ∆ is the gas over-density which we set to 200 (e.g.Mo et al. 2010).Furthermore, c mh is the clumping factor of the gas inside a mini-halo; in the simulations discussed below, we find typical values for c mh in the range 2-4.The values of these three times are similar for our default choice of parameters, e.g. a halo mass of M h = 10 5 M⊙ and a photo-ionization rate with Γ−12 ∼ 0.1.The dependence on halo mass is the same for tic and tsc, but the redshift dependence differs; tr does not depend on M h .Depending on the values of z and M h , we identify the three following regimes: • Sound-speed limited regime: When tic ≪ tsc, the I-front races through a mini-halo so quickly that its gas cannot recombine nor react hydrodynamically.The gas is ionized and heated and will photo-evaporate13 on a timescale, t ∼ tsc.This is typically the case for low-mass haloes, at low redshift, and when Γ−12 is large.
• Ionization limited regime: When tsc ≪ tic, the I-front moves slowly and the gas can photo-evaporate as soon as it ionizes: this also implies that the total gas density tracks the ionized gas density.
• I-front trapping regime: If tr ≪ tic, the gas is photoionized but recombines very quickly14 .The I-front moves at sub-sonic speeds and eventually nearly stalls at the inverse Strömgren layer (Shapiro et al. 2004).This regime occurs in halos of mass, M h ≳ 10 6 M⊙, and in the central dense region of low-mass halos.
In principle, the duration of the photo-evaporation phase is not simply established by the time it takes to photoevaporate the most massive mini-halos.It is because lowmass mini-halos may contribute more to the recombination rate since they are much more numerous.We will use the simulations described below to estimate the duration of this phase.

(Post-reionization) Relaxed phase
Long (≫ 100 Myr) after being overrun by the I-front, halos with M h < M char are photo-evaporated and unable to accrete gas (Okamoto et al. 2008): such halos no longer contribute to recombination.Gas in the outskirts of more massive halos is highly ionized, with the inner parts self-shielded and neutral.These lls's and dla's determine the attenuation length of ionizing photons, and thereby the relationship between the emissivity of ionizing photons and the photoionization rate (e.g.Faucher-Giguère et al. 2009;McQuinn et al. 2011;Haardt & Madau 2012).
Cosmological rhd simulations are required to capture the propagation and evaporation phases and the transition to the post-reionization phase.To resolve photon sinks during these stages, such simulations must resolve mini-halos above the Jeans mass in a computational volume that is large enough to sample the rarer lls's and dla's that determine the mean free path of ionizing photons in the postreionization phase.This is a tall order, even if we neglect the even more challenging calculation of resolving the nature of the ionizing sources and the thorny issue of determining the fraction of those photons that can escape their natal cloud.In this paper, we focus on photon sinks during the eor: we simply inject ionizing photons into our computational volume at a specified rate and follow how these photo-evaporate gas out of small halos.

Code and numerical set-up
We simulate a periodic cubic cosmological volume using the publicly-available15 sph code SWIFT (Schaller et al. 2016(Schaller et al. , 2018)).Among the various implementations of the sph algorithms (e.g.Borrow et al. 2022) included in the code, we select the entropy-based version described by Springel & Hernquist (2002) and Springel (2005).
Radiation hydrodynamics is solved with the sph-m1rt two-moment method with a modified m1 closure16 , as described by Chan et al. (2021) a simulation suite in which we vary the redshift when we start injecting photons, zi, and the intensity of the radiation, Γ−12.The spectrum of the radiation is that of a black body of temperature, TBB = 10 5 K, and is treated in a single frequency bin using the grey approximation; this also means that we neglect any spectral hardening of the radiation.The optically thin direction of the Eddington tensor is taken to be along the initial direction of propagation of the I-front; this improves the ability of the method to cast shadows and handle self-shielding.To reduce the computational cost, we propagate radiation at a reduced speed of light, c (Gnedin & Abel 2001).In our implementation, c scales with the value of the smoothing lengths of the sph particles 17 .The interaction of radiation with matter is calculated with a nonequilibrium thermo-chemistry solver with hydrogen only (as in Chan et al. 2021).We include helium when calculating the heat capacity, but we do not consider its interaction with radiation.Note that we also neglect molecular hydrogen and other elements (see §5 for a discussion on the caveats of our approach).Our original rt implementation did not account for the cosmological redshifting of radiation.We describe and test in Appendix C our choice of co-moving variables.
Our implementation accounts for the decrease in the proper density of photons as the Universe expands, but it does not account for the increase in the wavelengths of these photons.
The mean free path of ionizing photons is short in the case 17 See Appendix B for more details on this 'variable' speed of light approximation.
we simulate here.Therefore, this is a reasonable approximation.
Our simulation suite does not include feedback from evolving stars.However, as very dense gas particles severely limit the simulation time-step, and since we do not include the correct physics for these high-density regions anyway, we simply convert gas particles into stars once their density exceeds a physical density of 10 hydrogen atoms per cm 3 and an over-density of ∆ = ρ/⟨ρ⟩ = 10 3 .These criteria are similar to the 'quick-Lyα' approximation used by, e.g., Viel et al. (2004), who pointed out that the impact of this approximation on the Lyα flux power spectrum is small (less than 0.2%).Moreover, the density of gas particles that turn into stars is higher than that of the regions that give rise to Lyman-limit systems.This indicates that our approximation is unlikely to affect the I-front speed or the photo-evaporation timescales of mini-halos (see further discussions on this in §5).
We generate the initial conditions at redshift z = 127, using the publicly-available music code (Hahn & Abel 2011).The adopted cosmological parameters are: h = 0.678, Ωm = 0.307, ΩΛ = 0.693, and Ω b = 0.0455, σ8 = 0.811, and ns = 0.961, where symbols have their usual meaning.The hydrogen and helium mass fractions are X = 0.752 and Y = 1−X, respectively.We use Eq.(1) to compute the temperature at the mean density, T ihm , which gives the normalization of the adiabat describing the temperature-density  Reed et al. (2007).Colours correspond to z = 10.2 (blue squares, and solid line), z = 7.9 (green circle and dotted line ), and z = 6.0 (red triangles and dashed line).Empty symbols show mass bins with fewer than 5 halos per decade in halo mass.The thin vertical lines on the right show the minimum mass of atomic cooling halos (Eq.4); halos in the grey region are not well resolved, containing fewer than 300 dark matter particles.Our fiducial simulation resolves mini-halos with more than 300 particles in a volume sufficiently large to contain a few halos in which gas can cool atomically.
relation of the particles in the initial conditions, Tini(ρ):

Simulation Suite
Our main objective is to study the photo-evaporation of mini-halos and how this impacts the progression of reionization.This requires sampling and resolving mini-halos from the Jeans mass, Mmin, to the atomic cooling limit, MACH, (see Fig. 1).We vary the redshift at which we start injecting radiation and the flux of the injected ionizing photons and run individual simulations until several 100 Myrs after the radiation injection.Simulation parameters are listed in Table 1.We motivate the choices as follows.
We set the dark matter particle mass of the fiducial simulations to mDM ∼ 100 M⊙, so that a halo of mass Mmin (see §2) is resolved with ∼ 300 particles.The baryonic content of such halos is resolved with roughly 30% accuracy (Naoz et al. 2009) and the overall clumping factor of the simulated volume is accurate to ∼ 10 − 20% (Emberson et al. 2013).
To include massive halos of mass ∼ MACH, we consider a fiducial linear extent of 800 ckpc.We demonstrate in Fig. D8 that such a volume contains more than one halo of mass MACH at z = 8, and 10 at z = 6.This linear extent also yields approximately converged values for the mean free path of ionizing photons at the end of the eor (for unrelaxed gas; see Emberson et al. 2013)  We perform simulations with three choices of the 'reionization redshift', zi ∼ 10, 8, and 6, (zi is the redshift where we start injecting ionizing photons from two opposing faces of the cubic volume at constant flux).This range covers approximately current observational estimates for the start and tail-end of the eor (e.g.Fan et al. 2006;Planck Collaboration et al. 2020).The values of F correspond to photoionization rates of Γ−12 = 0.03 − 0.3 (where F and Γ−12 are related by the frequency-averaged photo-ionization rate, as in Eq. 8).This range in Γ−12 is motivated by observational estimates (e.g.Calverley et al. 2011;Wyithe & Bolton 2011;D'Aloisio et al. 2018) as well as simulation results (e.g.Rosdahl et al. 2018).
We take the value of the reduced speed of light to be proportional to that of the I-front: this allows us to use a lower value of c at lower Γ−12, decreasing the wall-clock time of the simulations (see Appendix D for tests of numerical convergence).Once the majority of the ihm is ionized, i.e. at the start of the evaporation phase of the eor, most I-fronts are expected to be D-type and hence propagate locally at a speed comparable to the sound speed, which is much slower than the speed of light.Therefore, we set c = c/100 once the mass-weighted neutral hydrogen fraction drops below 5% (see the similar approach and convergence tests in D 'Aloisio et al. 2020).
Finally, note that the simulated volumes are all relatively small and are not necessarily representative cosmological volumes during the eor.Rather we think of them as selected patches of the Universe that are overrun by an I-front due to sources outside of the simulated volume.These patches are ionized at various redshifts, zi, and with a range of values of the photo-ionization rate, Γi.

Halo identification
We use nbodykit (Hand et al. 2018) to identify halos using the friend-of-friend algorithm (Davis et al. 1985), as simply connected regions with a mean density of ∼ 200 times the average density of the Universe.The halo mass function of our simulations at redshifts before zi is compared to that computed with the colossus python package19 (Diemer 2018) in Fig. 4 (we selected the Reed et al. 2007 fit).The agreement between the simulation results and the fit indicates that our small volume contains the expected number of mini-halos up to the mass of atomic-cooling halos, as desired.We have performed additional runs in which we vary the box size and numerical resolution, see Appendix D.
heating and thus photo-evaporation.However, if photo-heating is included, halos with virial mass below the characteristic mass found by Okamoto et al. (2008) will eventually photo-evaporate, with more massive halos determining the clumping factor.An accurate calculation of the clumping factor then requires an even larger volume.To study this situation, our simulation suite includes larger volumes simulated at lower resolution.

Overview
We illustrate the first two reionization phases -I-front propagation and mini-halo evaporation -using the M512z8G03 run.In this simulation, ionizing photons are injected after redshift zi = 7.87, with a constant flux equivalent to a photoionization rate of Γ−12 = 0.3.The value of zi falls around the midpoint of the eor as inferred from the Thompson optical depth (Planck Collaboration et al. 2016), and the value of Γ−12 is typical of the expected mean value during this time (Γ−12 = 0.3 • • • 0.6, D' Aloisio et al. e.g. 2018).Therefore, this setup simulates the history of a typical patch of the universe during the eor.The evolution of the other simulations is qualitatively similar to that of M512z8G03.Simulations with higher values of zi have fewer self-shielded clouds since structure formation is less advanced, and those at lower zi or higher values of Γ−12 have more self-shielded clouds.
Fig. 5 shows the simulated volume at z = 7.84, which is ∼ 3 Myr after ionizing photons were first injected (propagation phase, upper row), when the I-front has propagated over approximately 1/6 th of the simulation volume.Using Eq. ( 5) with αB = 0 for the speed of the I-front in a homogeneous medium, we would expect an I-front to cross the computational volume in a time ∼ ⟨nH⟩ L /F ∼ 6 Myr (where L is the proper box size and F is photon flux) if recombination could be neglected.Clearly, density inhomogeneities, and in particular mini-halos, slow down the I-front by a factor ∼ 2−3 on average.The top row shows how the initial planar I-front becomes corrugated.This is because high-density regions locally slow down the front, as opposed to low-density regions which increase the front speed.Sufficiently dense gas in halos may even stop the I-front casting a shadow of gas that remains neutral behind them.Downstream of the I-front, the highly ionized gas is also photo-heated, reaching a temperature of T ∼ 2 × 10 4 K in a timescale of order Γ −1 (photo-ionization rate; see Eq.8).The middle panel shows that the I-front is sharp, with the distance from where the gas is mostly ionized to where it is mostly neutral of around Therefore the ihm during the eor is either highly ionized or mostly neutral, to a good approximation.
The lower row of Fig. 5 corresponds to z = 7.58 (∼ 30 Myr from the start of photon injection) when more than 99% of the volume has been highly ionized.When the I-front has already crossed the computational volume, it leaves behind filaments and halos with mostly neutral and cold gas.The lower-left panel shows that these filaments are photo-evaporating, with gas expanding out of the shallow potential well once it is ionized and heated (see also Bryan et al. 1999).The flow velocity is comparable to the local sound speed.These expansion waves compress and heat the gas in the filaments' outskirts, as seen in the lower right panel, where hotter gas surrounds the expanding cooler filaments.The denser gas in mini-halos takes longer to fully photo-evaporate, and sufficiently massive halos retain most of their baryons.
This discussion suggests assigning gas to three categories: (1) the low-density ihm, e.g.voids; (2) the filaments of the cosmic web; (3) collapsed halos including mini-halos.These structures stand out in the top left panel of Fig. 5.To investigate how these structures are impacted by reionization and vice versa how they affect reionization, we proceed as follows.In §4.2, we investigate the impact of I-front on the ihm and filaments, and in §4.3, study how self-shielding keeps the central parts of mini-halos neutral.In §4.4,we will turn to the photo-evaporation of mini-halos.In §4.5, we will quantify how these small-scale structures impede reionization and the role of photo-evaporation/relaxation.

Response of the ihm to reionization
We examine the volume density distribution of our simulations in Fig. 6.Here we plot ∆ 3 PV (∆) as a function of log ∆.We show these curves for the total gas density (solid blue) and ionised gas density (dashed blue).First, we study the response of the ihm to the passing of I-front for the simulation M512z8G03 (left panel) quantitatively.Upon the passage of I-front, low-density (∆ < 5) gas got highly ionised and heated to T ∼ 2×10 4 K (see also Fig. 7).But this gas changes relatively little in volume density (since it cannot expand further).On the other hand, self-shielded gas at high density (log ∆ > 3) remains mostly neutral and cold, whose PV also is not strongly affected by reionization due to self-shielding.Photo-heating significantly reduces PV of intermediate densities, log ∆ ∼ 1 − 3, since this gas is not self-shielded and expanded after photo-heating.More detailed analysis (including particle tracking) is presented in Appendix F.
Other curves in Fig. 6 show the total gas density pdf's for all gas in halos (dashed orange), halos in the mass range ) and log M h [M⊙] < 5 (grey dashed).The top panels correspond to the beginning of reionization when the I-front is still traversing the volume; the lower panels are at redshift z = 5.Panels from left to right correspond to different values for zi and Γ−12, as indicated.The straight black dashed line shows PV ∝ ∆ −2.5 , which corresponds to the pdf slope of an isothermal profile (see Miralda-Escudé et al. 2000 and discussions below).
In the left panel, we can see how halos with log M h < 7 photo-evaporate, with their contribution to the high-density pdf log ∆ > 2, decreasing dramatically between z = 6 and z = 5.The intermediate density gas, 1 < log ∆ < 3, is mostly associated with the more massive halos, log M h > 8.
The middle panels correspond to a case with zi = 8 and a lower photo-ionisation rate of Γ−12 = 0.03 as compared to zi = 6 and Γ−12 = 0.3.Nevertheless, the pdf's at z = 5 are quite similar, with the most striking difference being the location of the upturn in the pdf, which occurs around log ∆ = 2.6 in the left panel and log ∆ = 2 in the middle panel.We note that this upturn is also the location where the gas turns from mostly ionised to mostly neutral.
The right panels correspond to a model with zi = 10 and Γ−12 = 0.3.Although this model has the same value of the photo-ionisation rate as the model in the left panel, the pdf's look strikingly different, with, in particular, very little gas at high densities.The reason for this becomes clear by comparing the top panels: the more massive halos have not formed yet by z ∼ 10, and these halos do not accrete gas once it is photo-heated.This shows that reionization has Radiation is injected at redshift z i ≈ 8. Panels from left to right are hydrogen gas density (n H ) (Mach number in the lower left), the neutral fraction (x HI ), and the gas temperature (T ).Panels labelled 'projected' correspond to projections of the computational volume; the others correspond to infinitesimally thin slices through the mid-plane.As the I-front propagates through the volume, gas becomes ionized and photo-heated, and mini-halos photo-evaporate.Photo-heating causes the expansion of the filaments.
a bigger impact on more massive halos by preventing them from accreting hot gas rather than by photo-evaporating their gas.Miralda-Escudé et al. (2000) found that the gas pdf is described well by the form PV ∝ ∆ γ exp(−∆ 4/3 /σ 2 ∆ ), where γ(= −2.5) and σ∆ are fitting parameters.Their model is based on simulations with uniform UV background and the optically thin approximation (Miralda-Escudé et al. 1996).
With full radiation hydrodynamics here, our result (Fig. 6) does not agree with Miralda-Escudé et al. (2000) (compared with the γ = −2.5 scaling in Fig. 6).Our pdf has a stronger dip at log ∆ ∼ 2, followed by a stiff profile at higher density.With RHD, we capture the self-shielding of gas inside halos, which is less affected by the radiation, so the dense gas can maintain a stiff profile.The ionised gas, on the other hand, is photo-heated and evaporated, which explains the dip at log ∆ ∼ 2.  2020).These results support that the selfshielding of gas is responsible for the deviation from Miralda-Escudé et al. (2000).
The processes described here can be visualised in the temperature-density diagrams of Fig. 7, before (upper panel) and after (lower panel) the passage of the I-front.The gas in the simulations has an initial entropy, T /ρ γ−1 , at the start of the calculations, set by the initial temperature of the gas (Eq.13).This initial entropy is shown by the diagonal dotted line.But the gas can increase its entropy and move upwards from the dotted line through shocks during structure formation.
Before reionization, the gas temperature remains below ∼ 10 4 K due to Compton and atomic line cooling.Hence, pre-reionization gas exists in the triangular-shaped region of Figure 6.Probability distribution by volume multiplied by ∆ 3 , ∆ 3 P V , as a function of the logarithm of the over-density, log ∆.In the top row, the I-front is still sweeping the computational volume.The lower row shows post-reionization (in practice, we show the last simulation snapshot at znow ∼ 5).The 'Total' lines refer to all gas in the volume, whereas the 'Halo' lines refer only to gas in halos.We also distinguish the contribution of gas in halos according to halo mass (differently coloured lines).Those labelled 'Ionized' correspond to ∆ 3 P V (∆)x 2 HII such that the area under the curve is proportional to the recombination rate per unit volume per decade in ∆.The left shaded regions correspond to densities too low to resolve accurately, and the right shaded regions indicate gas above our artificial 'star formation threshold' where the simulation is unreliable due to missing physics.The short dashed black lines indicate the power-law P V ∝ ∆ γ , for γ = −2.5 for comparison.This power law describes the scaling of Miralda-Escudé et al. (2000)'s model at high densities (see the main text).The low-mass mini-halos are photo-evaporated between the upper and lower panels.
the upper panel.But we note that by mass, most of the gas remains close to the entropy floor, as shown by the contours.
After the I-front has crossed the computational volume, gas with log ∆ ⩽ 2 is ionized and photo-heated to T ∼ 2.2 × 10 4 K. Gas at higher densities self-shields, remains mostly neutral, and is at log T [K] ∼ 2 − 4. The lower density gas with log ∆ ∼ 0−1 is mostly at T ∼ 2.2×10 4 K, but some gas is hotter, and some gas is cooler.This results in a diamondshaped region in Fig. 7.As we explained previously, the origin of the hotter gas is due to adiabatic compression and shocking of under-dense gas by filaments 20 .But clearly, some gas can also cool adiabatically, explaining the lower diamond region.
The connection with filaments is illustrated in more detail in Fig. 8, which shows a close-up of two filaments that are nearly parallel to the y-axis (plotted horizontally in this figure).The middle panel shows that both filaments are expanding in the x direction with a speed |vx| of up 20 We tracked particles with T > 9 × 10 4 K and log ∆ ∼ 0 − 1 back in time to investigate the evolution of their entropy.After being photo-ionized, their entropy changes slightly, whereas their density may change by one order of magnitude.This means that their temperature change is mostly due to adiabatic compression/expansion, with a small contribution from shock.This is consistent with the scatter in temperature being approximately symmetric around the T ∼ 2.2 × 10 4 K. to 20 km s −1 (which is higher than the sound speed).This adiabatic expansion cools the central region of the filament and heats the surroundings ihm through compression (upper panel).Similar adiabatic changes are seen in the expansion of mini-halos and the resulting compression in their surroundings, e.g.Test 7 in Iliev et al. (2009).

The inverse Strömgren layer
When an I-front overruns neutral gas in a halo, its denser interior may be able to stall the front, provided the recombination rate is high enough.In mini-halos, the front will continue to move inwards as the gas in the outer parts photoevaporates, but in more massive halos, the centre may continue to self-shield.In this section, we estimate the density, nS, at the location of this 'inverse Strömgren layer', above which gas self-shields and remains mostly neutral21 .In the simulations, we determine nS as the minimum gas density where x > 1/2.The value of nS is of interest because it impacts the overall recombination rate and the clumping factor in a patch of Universe.A good understanding of what sets nS An analytical estimate for nS can be derived from the model of Theuns (2021) for damped Lyman-α systems.Assume that the density profile in a halo of virial radius R h is the power law where n H,h (= 200/3⟨nH⟩) is the hydrogen density at the virial radius.Assuming that the gas is in photo-ionization  equilibrium, the rate of change of the neutral density is zero, where x is the neutral fraction and we have neglected collisional ionizations and any contribution from helium or other elements.The photo-ionization rate Γ at radius R in the cloud is related to its value Γ h at R h by where the optical depth Combining these yields Taking the logarithm on both sides yields and taking the derivative with respect to R Figure 10.The transition to neutral gas in the model of Theuns (2021).Filled symbols denote the density at the inverse Strömgren radius, computed from Eq. ( 28).Open symbols denote the minimum density where x > 1/2, obtained from integrating Eq. ( 22).In both cases we used Eq. ( 30) to relate the photoionization rate at the virial radius, Γ h , to its volume-averaged value, Γ 0 ≡ 10 −12 Γ −12 s −1 .The symbol type encodes virial mass, M h in units M ⊙ , and its colour encodes redshift, as per the legends.The black dotted line shows the scaling ∝ Γ 2/3 0 .We have evaluated the recombination coefficient at T = 10 4 K and set σ HI = 1.62 × 10 −18 cm 2 .This is a differential equation for x(r), which we write in dimensionless form as where τ h ≡ σHI n H,h R h is a characteristic optical depth for the halo, and r ≡ R/R h .The boundary condition is the where the ionization and recombination time at r = 1 are t h ≡ Γ h −1 and t h ≡ (αBn H,h ) −1 .This differential equation has no closed-form solution, the numerical solution is plotted in Fig. 9 for a range of halo masses and two redshifts.In the ionized outskirts of the halo, we can take x ≪ 1 so that 1 − x ≈ 1, and the differential equation simplifies to with solution This approximate solution is also plotted in Fig. 9.Not surprisingly, it captures the increase in x with increasing nH very well when x ≪ 1, but it also captures rather well the density and location where x = 1/2 and even where x = 1.At a given redshift, nS (where x = 1) is higher for lower halo masses.For a given halo mass, nS increases with increasing redshift.
We obtain the scaling of nS with M h and z as follows.We start by computing the value RS of the inverse Strömgren layer by expressing that the recombination rate along a ray, from R h to RS, equals the impinging flux -this simply means that all photons impinging on the halo at R h have been used up by recombination between R h and RS: so that for our 1/r 2 density profile We note that tr depends on redshift but not on halo mass, whereas τ h depends both on M h and redshift.The approximation of neglecting the '1' in the round brackets applies to most cases of interest.In this approximation, we derive the following scaling relation for the hydrogen density nS at the Strömgren radius, with the numerical value taking the gas temperature to be T = 10 4 K when evaluating the recombination coefficient, and taking σHI = 1.62 × 10 −18 cm −2 .We note that the dependencies on halo mass and redshift are relatively weak.
So far, we characterised the photo-ionization rate by its value Γ h at the virial radius.It is likely that Γ h < Γ0, where Γ0 is the volume-averaged photo-ionization rate because the gas in the surroundings of the halo also causes absorption and hence suppresses the ionizing flux.We can estimate the importance of this effect by simply extrapolating the 1/r 2 density profile to infinity: This yields the following relation, We plot nS computed from Eq. ( 28) and use the previous equation to relate Γ h to Γ0 in Fig. 10.The transition from neutral to ionised gas is relatively sharp in the simulations, as can be seen by comparing the 'Ionised' and 'Neutral' lines in Fig. 6.This justifies the Miralda-Escudé et al. ( 2000) assumption for a characteristic density dividing mostly neutral and mostly ionised gas, which was also shown in other numerical studies, e.g.Mc-Quinn et al. (2011), Park et al. (2016) andD'Aloisio et al. (2020).Therefore, the value of nS is relatively well defined, and operationally we determine it for each halo as the minimum density for which x > 1/2.We plot this value in Fig. 11 as a function of Γ−12.The analytical relation of Eq. ( 28) captures the nS ∝ Γ 2/3 −12 scaling, but the simulated value of nS is smaller.We suspect this is due to the density structure in the accreting gas, which the analytical model does not account for.
We generalise the analytical calculation to the case of a halo overrun by a plane-parallel I-front in Appendix E. We use this to test the ability of the rt scheme to capture an I-front at the typical numerical resolution with which we simulate mini-halos.

Photo-evaporation of mini-halos
Here we examine the time-dependent effects of the ionizing radiation on individual halos.To do so, we need to calculate how long mini-halos have been irradiated with ionizing radiation.
In our simulations, we inject ionizing photons from two opposing sides of the periodic volume.As a consequence, mini-halos close to the injection side are affected by the radiation earlier (and for longer) than those further away.To account for this time difference, we first calculate the I-front position, defined as the location where the volume-weighted neutral fraction is ∼ 0.5 (we do this in cubic cells of extent 10 ckpc instead of considering a plane-parallel I-front).For each halo, we then record the time since it encountered the I-front.Hence we can compute the duration t dur between the current time and the time that the halo was first irradiated.
We select ∼ 20 mini-halos with mass M h ∼ 10 6 M⊙ and compute the median of their spherically averaged radial density profiles.We assume that the centre of the halo corresponds to the location where the density of neutral gas is highest 22 .
The density profiles of these halos are plotted in Fig. 12 for various values of t dur , and two values of the photoionization rate, Γ−12 = 0.3 and 0.03.
If radiative cooling is inefficient (as in mini-halos), the maximum gas density is limited by the entropy of the IGM after decoupling from radaiation (Visbal et al. 2014).Therefore the central density profiles are flat.Outside of the central regions, the profiles of the total density (solid lines) Figure 12.Gas over-density profiles of halos with M h ∼ 10 6 M ⊙ at different times, t dur , after they were over-run by an I-front.The solid lines are total gas density, whereas the dashed lines are neutral gas density.The arrows indicate the 15 th and 85 th percentiles of halos with similar values of t dur .The black dashed lines show power-laws for isothermal profiles, ρ(∆) ∝ ∆ −2 , that describe the run of the total density with radius in the range 0.1 ⩽ r/R h ⩽ 1 relatively well (R h is the virial radius of the halos).Gas in the shaded region is within a gravitational softening length from the centre.
are reasonably well approximated by nH(r) ∝ r −2 (black dot-dashed line), at least before the I-front has propagated significantly into the halo's gas.
Therefore halos are in the sound-speed limited regime for Γ−12 = 0.3 (upper panel, using the nomenclature of §2.2.2).The I-front propagates rapidly into halos, reaching down to 5 per cent of the virial radius within 50 Myr.Gas starts photo-evaporating, but the time it takes to leave the halo is longer than tic.Therefore the halo contains a large amount of highly ionized photo-heated gas.Given that the gas has density profile is nH ∝ r −2 , the neutral gas has a density profile approximately nHI(r) = x nH(r) ≈ αB/Γ−12n 2 H ∝ r −4 .In contrast, halos are in the ionization limited regime when Γ−12 = 0.03 (lower panel).The I-front propagates so slowly into the cloud that the photo-heated gas has time to photo-evaporate and leave the halo.Consequently, the neutral and ionized gas profiles almost trace each other.For this low value of Γ−12, even halos of this low mass can hold on to a large fraction of their gas for several 100Myr, and this affects the duration of the photo-evaporation phase of the eor.
We demonstrate how fast a halo photo-evaporates in A more quantitative measure of the photo-evaporation timescale is plotted in Fig. 14.We have computed the value Thin dotted lines indicate that evaporation is not yet complete when the simulation is stopped.The lower panels are similar, except that halos must lose 99.9 per cent of their gas.Triangles are the simulation results of Iliev et al. (2005a).Our evaporation times are slightly longer, but the trends with M h and Γ −12 are very similar: the evaporation time increases with M h and z i , and decreases with increasing Γ −12 .
of t dur after which a halo contains only 10 per cent of the cosmic baryon fraction, which we refer to as t ev,10% (upper panel).For Γ−12 = 0.3, we find the approximate scaling t ev,10% ∝ M 1/2.3 h (solid lines in the top right panel), which is somewhat shallower than the tic ∝ M 1/3 h dependence of the I-front crossing time on halo mass.That panel also shows that the evaporation time depends surprisingly little on redshift.The dependence is even weaker than that of the sound crossing time, tsc ∝ (1 + z) −1 .t ev,10% can be a few times tsc, because recombinations significantly delay ionization and hence photo-evaporation.
In the left and central panels, we notice that the dependence of t ev,10% on Γ−12 is weaker than that of tic ∝ Γ −1 −12 .At the lower values of Γ−12 = 0.03, the more massive halos have not yet reached t ev,10% by the end of the simulation run: this is indicated by the thin dotted lines.
The dashed lines indicate when the halos contain less than 10 per cent of neutral gas.When Γ−12 is low (e.g.Γ−12 = 0.03), there is little difference between the total and neutral gas lines because once ionized, gas quickly leaves the halo: these halos are in the ionization-limited regime.However if Γ−12 is larger (Γ−12 = 0.3, say), some of the ionized gas is still inside the halo because it has not had time yet to photo-evaporate: these halos are in the sound-speed limited regime.The difference between the two regimes is more pronounced at lower z and lower halo mass, as seen in the central panel.
The lower panel of Fig. 14 shows t ev,0.1% , the time after which halos contain less than 0.1 per cent of the cosmic baryon fraction.The panels compare our results to those of Iliev et al. (2005a) (triangles), and we find an agreement within a factor of two.Several differences between our simulations and theirs may explain the difference.Firstly, ours are 3D cosmological rhd simulations, in which halos grow in mass through mergers and accretion during photoevaporation.In contrast, theirs are 2D non-cosmological simulations.Secondly, our simulations include the effect of shadowing since we simulate a cosmological volume.Third, we consider radiation in one frequency bin, whereas their radiative transfer is multi-frequency.The multi-frequency treatment will allow high-energy radiation to penetrate further into halos and include effects of spectral hardening (which boost the photo-evaporation rates).Finally, our sph simulation is particle-based, whereas they used a uniform grid.Our spatial resolution is comparable to theirs at high density but is lower in lower-density regions.
We do not plot the evaporation times of the total gas (solid) in the lower panels.It is because, in our simulations, the total gas fraction does not drop below 0.1 per cent: even mini-halos with M h ∼ 10 4 M⊙ can hold on to a small fraction of (highly ionized) gas.

The clumping factor, cr
The recombination rate in an inhomogeneous medium is enhanced compared with that of a uniform medium of the same mean density by the clumping factor, cr23 , which we define here as with the extra subscript 'all' added for reasons explained below.Here, angular brackets, ⟨⟩, denote volume-averaging and T is the mass-averaged temperature.This definition captures the dependence of the recombination rate on temperature (see also Park et al. 2016).Other definitions have appeared in the literature: these clumping factors can differ by tens of per cent, but the qualitative trends are similar (see Appendix D for details).
Structures change the I-front speed and location; inserting Eq. (31) into Eq.( 5) yields the solution in the approximation that nH, tr, and F are all time independent.As expected, the I-front moves slower and stalls earlier when clumping is more pronounced (greater cr).
We cannot apply Eq. ( 31) directly to our simulations because (i) we inject radiation from the sides into the computational volume, and (ii) we include collisional ionization in the calculation.To account for (i), we only include gas downstream from the I-front in the calculation of cr: We approximate the I-front by a plane where the volumeweighted neutral fraction first reaches 50 per cent (and perpendicular to the injection direction).We also resolve (ii) this way since collisional ionization contributes relatively little to the recombination rate in the downstream region, where photo-ionization dominates.Furthermore, our small simulation volume does not contain many more massive halos where collisional ionizations are frequent.

Evolution of the clumping factor
Fig. 15 shows the evolution of the clumping factor as computed using Eq. ( 33).In all runs, cr increases rapidly to a peak value24 , and after ∼ 100 Myr, starts to decrease to an asymptotic value of ≲ 5.The height of the peak, c r,peak , increases with increasing Γ−12 and with decreasing zi.We find that approximately where the exponent γ ∼ 0.5−0.6.The increase in c r,peak with decreasing zi is simply due to structure formation: there are far fewer halos at higher z to boost the recombination rate.
The dependence on Γ−12 is a result of the increase of nS (the density where the recombination rate is sufficiently high that the gas remains significantly neutral, Eq. 27) with Γ−12 and the dependence of the photo-evaporation time scale on Γ−12.
We caution the reader that our small simulation volume does not contain sufficient massive halos that dominate cr post-reionization.This may affect the simulation with zi ∼ 10 most because photo-heating of the ihm at high z prevents halos from accreting gas.
The contribution to the recombination rate of gas at different densities is plotted as the blue line, labelled as 'Total' in Fig. 16.At low reionization redshift zi, high Γ−12, and close to the reionization redshift (z ∼ zi), the recombination rate is dominated by the denser gas with ∆ > 100 that is mostly inside of halos (upper left panel).At slightly later times (lower left panel), lower-density gas with ∆ of the order of a few starts to dominate because the denser gas has photo-evaporated (in fact, a considerable fraction of the gas with initial over-density ∆ ∼ 100 decreases in density to values of a few following its photo-evaporation, see the discussion of Fig. F1).
Our small simulation volume underestimates the contribution of halos to recombinations to some extent.Still, we find that similar trends (a rapid rise in cr to a peak value, followed by a decrease on a longer time scale) in a simulation performed in a larger computational volume (see Appendix D).
The central panels in Fig. 16 also demonstrate how photo-evaporation decreases the contribution of halo gas to the recombination rate.Due to photo-evaporation, gas around the mean density, ∆ ∼ 1, increasingly dominates the recombination rate at later times.This is even more striking in the right panels of the figure: early reionization prevents the accretion of gas onto more massive halos that form after zi so that halos now contribute much less to recombinations.
Figure 16.Total recombination rate per hydrogen, ṅrec/nH, as a function of gas over-density ∆, per decade in ∆.The area under the curve measures the contribution of gas at that over-density to the recombination rate in the simulation volume.Line styles and colours are the same as in Fig. 6.In a given column of panels, the upper row corresponds to z ⪅ z i , whereas the lower row corresponds to a lower redshift.Gas in halos contributes significantly to recombinations when z ⪅ z i and Γ −12 is high.As mini-halos photo-evaporate, the contribution of the ihm to the overall recombination rate increases.At even lower redshifts (not shown), recombinations in more massive halos not contained in our small simulation volume dominate eventually.
The other coloured lines in Fig. 16 quantify the contribution to the recombination rate in halos of a given mass range, as indicated by the legend.For our choice of photoionization rate, Γ−12, and reionization redshift, zi, we find that mini-halos with mass below 10 6 M⊙ (grey dashed line) do not contribute significantly to the recombination rate, and by extension to the clumping factor.This is partly because they photo-evaporate quickly, as seen by comparing the upper panel with the lower panel in any given column.Rather, the more massive mini-halos with masses ≳ 10 6−7 M⊙ dominate recombinations during the eor.
Our findings agree with those of Iliev et al. (2005a), who also concluded that mini-halos with M h ≪ 10 6 M⊙ contribute little to cr.They find that halos with mass less than ∼ 10 6 M⊙ increase25 the photon consumption during eor by less than a 10-20 per cent.Since low-mass mini-halos do not affect cr significantly, we could relax the required resolution of our simulations to mDM ≲ 10 3 M⊙ (see Appendix D).
We now return to the late-time behaviour of the clumping factor, as plotted in Fig. 15.The value of cr drops from its peak to an asymptotic value of ∼ 2 − 3 after several 100 Myrs, in agreement with previous simulations of the post-eor ihm (e.g.Pawlik et al. 2009;McQuinn et al. 2011).The duration of the characteristic drop is set by the photo-evaporation time of mini-halos, which is of the order of 100 Myrs for the more massive mini-halos and ≲ 50 Myr for mini-halos with mass M h = 10 5 M⊙ (see Fig. 14).
Interestingly, we find that the value of the clumping factor at the end of the eor is relatively insensitive to the photo-ionization rate.For example, in Fig. 15, changing Γ−12 by over an order of magnitude does not affect the late-time clumping factors by more than ≳ 10 per cents.This was also seen in the simulations of D' Aloisio et al. (2020).This finding appears unexpected at first: larger values of Γ−12 cause the ionized-neutral transition to occur at higher densities (Eq.27) and this increases the recombination rate.Therefore if gas in halos were to dominate the recombination rate (and hence cr), then we would expect that cr ∝ Γ 1/3 (e.g.McQuinn et al. 2011).However, we find that recombinations mostly occur in gas with density around the mean at these early times, as seen in the lower panels of Fig. 16.Since this gas is highly ionized in any case, changing the value of Γ−12 has little effect on cr.
We now compare our results to those obtained by D'Aloisio et al. ( 2020) and Park et al. (2016).Note that the evolution of the clumpy factor depends strongly on the simulation setup and the reduced speed of light, so we only make qualitative comparisons here.Our values of cr are apparently ∼ 50 per cent higher than those plotted in  factor of 3 lower than that used by them, Γ−12 = 0.92.Applying the scaling of Eq. ( 34) with γ = 0.5, we would expect c r,peak ≈ 4.7 × (0.92/0.3) 1/2 ≈ 8.2 for the case of zi = 10 and Γ−12 = 0.92, which is reasonably close to their value of c r,peak = 7.5 for that value of Γ−12 and zi in their run M_I-1_z10.We also agree that the suppression of cr due to photo-evaporation is more significant at lower Γ−12 because the gas in mini-halos starts to photo-evaporate before radiation ionizes the higher-density gas.

Photon consumption in the clumpy medium
In Fig. 17, we plot the ratio Nr ≡ Nrec/NH, where Nrec is the cumulative number of recombination in the simulation volume up to a given time, and NH is the total number of hydrogen atoms in that volume.We computed and recorded this ratio as the simulation progressed.We find that Nr increases with Γ−12 (different line styles) but is rather strikingly only weakly dependent on zi (different colours).
Since c r,peak ∝ Γ γ −12 , a higher photo-ionization rate increases the clumping factor and hence also Nr, as expected.The weak dependence of Nr on zi was also noted by Park et al. 2016.The mean instantaneous recombination rate per hydrogen atom is where the integral is over the computational volume.Since we found that approximately cr ∝ (1+z) −3 (Eq.34), whereas ⟨nH⟩ ∝ (1 + z) 3 , the redshift dependence of the mean recombination rate per hydrogen atom is weak, which may explain the surprisingly weak dependence of Nr on zi Finally, we conclude that recombinations in the ihm and mini-halos increase the number of ionizing photons per hydrogen atom required to ionize the Universe by a factor 1 + Nr ≈ 1.2 − 1.7.The lower value corresponds to Γ−12 = 0.03, the upper value to Γ−12 = 0.3.Hence the claim that recombinations increase 1 + Nr by 20-100 per cent, which we made in the abstract.

DISCUSSION
Our simulation setup does not account for the formation of molecular hydrogen (H2) and hence neglects cooling due to H2.The simulations of Shapiro et al. (2004), Emberson et al. (2013), Park et al. (2016), andD'Aloisio et al. (2020) made the same assumption.Molecular hydrogen cooling can trigger star formation in halos of masses as low as 10 6 M⊙ at z ∼ 10 (e.g.Tegmark et al. 1997), orders of magnitude below the mass of halos in which gas can cool atomically (MACH discussed in § 2).Feedback from the first Pop iii stars may remove gas from some mini-halos, thus reducing their impact on the propagation of the I-front and the value of cr.However, molecular hydrogen can be photo-dissociated by Lyman-Werner (LW) photons (Stecher & Williams 1967), which are emitted by those first stars (Haiman et al. 1997).An LW radiation background may thus have already turned off the molecular cooling channel at the redshifts we considered (z ≲ 10, e.g.Trenti & Stiavelli 2009), suppressing Pop.iii star formation and hence reducing its impact on mini-halos26 .
We treat ionizing radiation using a frequency-averaged photo-ionization cross-section rather than following photons with a range of frequencies.This approximation does not capture the effects of spectral hardening or preheating.If included, these effects would most probably decrease the photo-evaporation time of mini-halos; hence, it may well be that our simulations overestimate the impact of such halos on the eor.For example, Iliev et al. (2005a) found that a harder spectrum can reduce the photo-evaporation time by ∼ 50%.
We neglected sources of ionizing photons inside the computational volume.Instead, we injected photons from two opposing faces of the cubic simulation volume.Previous studies also inject radiation from boundaries and neglect the sources inside the volume (Emberson et al. 2013;Park et al. 2016;D'Aloisio et al. 2020).This approximation is reasonable for our purposes, given the small size of the computational volume.
Our simulation volume contains halos with mass greater than MACH that would plausibly be able to form stars and contribute to the ionizing flux.Counter-intuitively, including stellar feedback from evolving stars forming in these halos could potentially increase the net recombination rate.Indeed, McQuinn et al. (2011) found that the recombination rate is higher in a model that included stellar winds compared to a pure rhd run.Feedback can also avoid the overcooling problem, where too much gas turns into stars (Katz 1992).Simulations with radiation sources within the simulation volume would be an important future improvement.These would be necessary for simulations of larger cosmological volumes where spatial correlations between photon sources and sinks are important.
Our simulations do not include pre-heating by X-rays emitted from an early generation of accreting black holes, as envisioned by Ricotti & Ostriker (e.g. 2004).Pre-heating would increase the minimum mass of halos that contain gas before reionization, reducing the impact of mini-halos on the eor.D'Aloisio et al. ( 2020) presented a simulation where the ihm is pre-heated by X-rays below redshift z = 20, and found a factor of two suppression of c r,,peak for a simulation with zi = 8.They also concluded that X-ray pre-heating did not affect the value of cr at the end of the eor.Park et al. ( 2021) also studied the impact of X-ray preheating (starting at z = 20) and found a much larger decrease in the value of c r,peak .However, their small simulation volume (200 h −1 ckpc) does not sample the more massive mini-halos that may be less affected by X-ray pre-heating, which may lead them to overestimate the impact of pre-heating.Fialkov et al. (2014) suggested that X-ray pre-heating most likely occurs at lower redshifts, not long before zi, which would reduce its impact on mini-halos.
We have neglected the relative velocity between dark matter and baryons post recombination (dark matterbaryon streaming), as discussed by Tseliakhovich & Hirata (2010).The streaming velocity can be several times larger than the pre-reionization sound speed, and this effect suppresses the early formation of low-mass halos and affects their baryon contents.Studying this effect, Park et al. (2021) found a cr lower by a factor of two.In contrast, Cain et al. (2020) concluded that baryon streaming decreases cr by only 5-10 per cent (in regions with the root-mean-square streaming velocity).They claimed that the small simulation volume used by Park et al. (2021) exaggerates the streaming effects 27 .
In summary, the previous discussion highlights some of the various ways in which our simulations could be improved in the future.Nonetheless, these are unlikely to change our main conclusions substantially.

CONCLUSIONS
We presented cosmological radiation hydrodynamics simulations of the propagation of ionization fronts in a cosmological density field.Radiative transfer is performed with a two-moment method with local Eddington tensor closure, implemented in the swift smoothed particle hydrodynamics code, as described by Chan et al. (2021).Our simulations follow the ionization and heating of the gas, leading to the photo-evaporation of (gas in) low-mass halos (mini-halos; Tvir < 10 4 K; M h ≲ 10 8 M⊙), and self-shielding of gas in more massive halos.The simulation can resolve the smallest mini-halos that contain gas before reionization (gas mass resolution mgas ∼ 20 M⊙; spatial resolution ∼ 0.1 ckpc).We inject ionizing photons from two opposing sides of the computational volume at a given photo-ionization rate,Γ−12, and for a given specified 'reionization redshift', zi.We assume that the ionizing sources have a black body spectrum with temperature TBB = 10 5 K.
These simulations improve upon the 2D rhd simulations performed by Shapiro et al. (2004) and Iliev et al. (2005a), who studied a more idealised set-up of a single minihalo overrun by an I-front assuming cylindrical symmetry (see also Nakatani et al. 2020).In contrast to previous cosmological rhd simulations of the intergalactic medium (Park et al. 2016;D'Aloisio et al. 2020), we have either a larger volume or higher (spatial) resolution.We also quantify the relative contribution to the clumping and recombination of mini-halos and the ihm.The results are presented in §4 and summarised in §6.
Our main results are as follows.In terms of the impact of the passage of the I-front on the cosmological density distribution, we find that • Upon the passage of I-front, the low density ihm is rapidly heated to T ∼ 2 × 10 4 K.This value is expected when a black body spectrum of photons with temperature TBB = 10 5 K flash ionizes hydrogen when non-equilibrium effects are accounted for (Chan et al. 2021), but spectral hardening is neglected.Higher density filaments expand supersonically, heating the surrounding ihm to T ∼ 10 5 K through shocks and adiabatic compression.Gas also cools through adiabatic expansion.As a result, the temperature of ionized gas scatters around T ∼ 10 4 K. Gas at higher densities in mini-halos slows the I-front as dense gas remains neutral due to self-shielding.
• At later times, mini-halos photo-evaporate, which moves gas from an over-density of ∆ ∼ 100 to ∆ ∼ 5 (see §4.4).The photo-evaporation timescale becomes shorter with increasing Γ−12, lowering zi and the mini-halo mass.The evaporation time is of order 50 − 100 Myrs, which can be several times longer than the sound crossing time.This is the case for halos of mass M h ≳ 10 6 M⊙, which trap the I-front for a long time (≫ 10 Myr) due to the short recombination time of their gas.
• Our estimates of minihalo photo-evaporation times agree qualitatively with those of Shapiro et al. (2004) and Iliev et al. (2005a), who performed high-resolution idealized 2D radiation hydrodynamics simulations.
In terms of the impact of mini-halos on the overall recombination rate as quantified by the clumping factor, cr, we found in § 4.5 that: • cr increases rapidly up to a peak value of c r,peak ∼ 20 as the I-front overruns mini-halos.The value of c r,peak is higher when Γ−12 is higher or zi is lower ( § 4.5.2).As minihalos photo-evaporate over a time-scale of ∼ 100 Myrs, cr decreases to a value of ≈ 2 − 4. At this stage, most recombinations occur in the ihm rather than in halos.Consequently, at the final stages of the eor, the value of cr depends only weakly on the value of Γ−12.
• The inhomogeneous ihm increases the reionization photon budget by 20-100 per cent, depending on the value of the photo-ionization rate during the eor (higher values of Γ−12 increase the budget) and to a lesser extent the value of zi (see Fig. 17).
• Low-mass mini-halos of mass ≪ 10 6 M⊙ do not contribute significantly to cr or the reionization photon budget because they photo-evaporate very quickly.The relative contribution to the recombination rate as a function of halo mass is plotted in Fig. 16.
We envision the following avenues for further research.Firstly, investigate the effect of the photo-evaporating minihalos on the value and evolution of the mean free path of ionising photons during the eor.What is the nature of absorbers that dominate the opacity (see, e.g.Nasir et al. 2021)?Does this explain the rapid evolution claimed by Becker et al. (2021)?Do these absorbers leave a trace on Lyman α forest (Park et al. 2023)?Secondly, study the origin of the power-law-like gas density profile (ρ ∝ r −2 ) in mini-halos pre-reionization.Thirdly, predict the shape and evolution of the gas density probability distribution (pdf), as well as the evolution of the pdf of neutral gas before, during, and after the eor.Finally, it would be worthwhile to incorporate the results presented here in a sub-grid physics model for recombinations (see e.g.Cain et al. 2021Cain et al. , 2023 for effects along these lines).Such a model could then be applied to a simulation at a much lower resolution but in a much larger volume.This would allow for the incorporation of mini-halos in reionization simulations that can connect with observations of the eor in upcoming 21-cm surveys.
While this paper focuses on the reionization of a clumpy universe, it also serves as a benchmark of the two-moment method SPHM1RT in cutting-edge cosmological simulations.We are coupling our radiation hydrodynamics method to interstellar medium and galaxy models.It will be promising in addressing a multitude of astrophysical problems, ranging from active galactic nuclei, HII regions, to selfconsistent reionization.drey Kravtsov, and Paul Shapiro.
We would like to pay our gratitude and respects to the late Prof.Richard Bower, who passed away in January of 2023.Richard was a Professor in the Physics department at the Durham University, with an expertise in understanding cosmology and galaxies with semi-analytic methods and simulations.He played a pioneering and leading role in many world-class astrophysics projects, including GALFORM, EAGLE, and SWIFT.He was also an excellent supervisor and mentor for numerous students and postdocs, who continue his legacy of inspiration in academia and industry.
We are greatly indebted to Richard for his innovative ideas and advice on this and related projects.
This work was supported by Science and Technology Facilities Council (STFC) astronomy consolidated grant ST/P000541/1 and ST/T000244/1.We acknowledge support from the European Research Council through ERC Advanced Investigator grant, DMIDAS [GA 786910] to CSF.
TKC is supported by the 'Improvement on Competitiveness in Hiring New Faculties' Funding Scheme from the Chinese University of Hong Kong (4937210,4937211,4937212). TKC was supported by the E. Margaret Burbidge Prize Postdoctoral Fellowship from the Brinson Foundation at the Departments of Astronomy and Astrophysics at the University of Chicago.
ABL acknowledges support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (GA 101026328).
This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk).The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1.DiRAC is part of the National e-Infrastructure.
The research in this paper made use of the swift opensource simulation code (http://www.swiftsim.com,Schaller et al. 2018) version 0.9.0.This work also made use of matplotlib (Hunter 2007), numpy (van der Walt et al. 2011), scipy (Jones et al. 2001), swiftsimio (Borrow & Borrisov 2020), nbodykit (Hand et al. 2018), colossus (Diemer 2018), and NASA's Astrophysics Data System.We thank an anonymous referee for their careful reading of the paper and their valuable comments on its content.
Because the gas temperature in the ihm is finite, the gas density distribution is smoother than that of the dark matter.However, the gas smoothing depends on the entire temperature evolution.
The filtering scale is where the baryonic perturbation are smoothed, given the temperature evolution.The filtering scale can be understood in terms of an evolving Jeans length, as described by Gnedin & Hui (1998).First define the (comoving) Jeans wave-number, kJ , at a given scale factor (a) as, where ρ m,0 = Ωmρc is the average matter density at z = 0, and cS is the sound speed.The filtering wave-number, kF , is obtained by integrating Jeans wave-number over time: where D+ is the linear growth factor.
In a flat Universe, D+ is equal to the scale factor, a, so D+ ∝ a ∝ t 2/3 at redshifts z ≳ 2. Therefore where z dec = 1/a dec − 1 ≈ 130 is the redshift where the gas temperature decouples from the cmb temperature.Before reionization and the onset of any other heating sources, the ihm temperature is: Therefore the (co-moving) Jeans scale is independent of redshift before decoupling.
We combine Eqs.A3 and A6 to calculate the filtering wave-number: In the simulations, we have neglected the difference between the baryon and dark matter power spectra at the starting redshift.Thus to enable a fair comparison with the simulations, we only integrate Eq. (A3) from a dec to find, 3 ln(a/a dec ) − 6 + 6 a dec a The filtering mass (eq3) can then be inferred from the resultant kF and MJ .When generating initial conditions for the simulations, we assume that gas traces the dark matter.For consistency, we therefore replace z dec by the starting redshift of the simulations, zic = 127, to estimate MF in the simulations, with the result plotted in Fig. 1.The resulting value of MF is lower than that quoted by Gnedin & Hui (1998), who start the integration in the limit of z → ∞.Our estimate of MF is close to the more accurate calculation performed by Naoz & Barkana (2007), who take account of the fact that the gas and cmb temperatures do not decouple completely at z dec .Presumably this better agreement is mostly coincidental.

APPENDIX B: VARIABLE SPEED OF LIGHT
The time-step in simulations is limited by the Courant-Friedrichs-Lewy condition (CFL, Courant et al. 1928), ∆t < Ch/vsign, where h is the spatial resolution, vsign the signal speed in the simulated fluid, and C is a constant of order unity.For radiation, the signal speed is the speed of light (vsign = c), so this time step is often many orders of magnitude smaller than that for the fluid itself, where the signal speed is of the order of the sound speed.Therefore, Figure B2.Tests of the rsl and vsl implementation, in which a packet of radiation streams in the optically thin limit from left to right.Shown in all cases is the product of the energy density times c, which is independent of the choice of c.The black dashed line is the initial shape of the packet, with thin coloured lines the correct solution at later times.The corresponding numerical solution in the rsl case with c = 1 is shown as red downwards triangles.The blue stars and green circles correspond to the numerical solution in the vsl approximation, which uses c = 1 for x < 0 and c = 0.5 and c = 2.0) for x > 0, respectively.The numerical calculation uses 10 3 particles in 1 dimension, and in all cases is shown at time t = 1.The simulation reproduces the correct solution, apart from the inevitable diffusion associated with the SPH scheme.
the inclusion of rt can dramatically reduce the time step and, in turn, increases the compute time.
However, even though radiation travels with speed c, ionization fronts typically travel much slower, see for example Eq. ( 7).Gnedin & Abel (2001) first suggested that reducing the speed of light, c → c ≪ c, would improve the computational speed of the code without sacrificing its accuracy: this is the 'reduced speed of light' approximation (hereafter rsl).The essential idea is that the state of the gas, when it is in equilibrium with the radiation (e.g.photoionization equilibrium), is independent of the propagation speed.In fact, the rsl approximation may give accurate results even in some non-equilibrium situations, provided the light crossing time is much shorter than other time-scales (see e.g.Rosdahl et al. 2013;Chan et al. 2021).
The required c varies across the simulation volume.Indeed, I-front moves fast in low-density regions, and hence c should not be much smaller than c.However, such a large c results in very short time-steps in high-density regions28 , where h is small.Fortunately, the I-front speed is actually very low in dense regions, so we could apply a much lower value of c without loss of accuracy.Therefore, it is very tempting to allow c to vary in space and time.This is the variable speed of light approximation (hereafter vsl, e.g. one vsl version in Katz et al. 2017).
The first requirement (1) of a vsl scheme is to keep the equilibrium photo-ionization rate independent of c, so does the equilibrium neutral gas fraction.The photo-ionization rate is proportional to the number density of photons, the speed of light, and the photo-ionization cross-section, i.e. nγ c σHI.This implies that the product of the radiation density times the reduced speed of light, nγ c, and the radiation flux should not change when c is used29 .
The second requirement (2) is that the scheme should reduce to the rsl approximation when c is uniform.
For reference, the terms entering the rate of change of radiation variables that involve spatial derivatives of these variables (and hence are computed by summing over sph neighbours that may have a different value of c), are given by (see Chan et al. 2021, their Eqs. 5-7), where the subscript 'prop' means that we only write terms that involve the propagation of radiation across particles.
We are going to modify these equations for vsl, so that they can be applied consistently even if two particles that exchange radiation have different values of c.To be definite, let these particles be i and j, with rsl values ci and cj.
To satisfy requirements (1) and ( 2), we devise a new three-step vsl scheme, as illustrated in Fig. B1.In short, (A) we transform the equations to a system with a common rsl, c0.Next, (B) we solve the propagation equation with c0.Finally, (C) we transform the equations back to the ci and cj.
In the first step (A), we construct a reference frame with c0 and then transform the quantities from the local (particle i) frame (with the rsl ci) to the reference frame through: where we label the c0 reference frame with the subscript "0".We apply the same transform for particles j, the neighbourhood particles of i.In this transform, ρf and ρcξ are invariant, so the radiation energy density in the reference frame is equal to that in the original frame, i.e ρc0 ξ0 = ρciξi, so this satisfies the requirement (1).Hence, the radiation equations in the c0 frame become: In the second step (B), we solve Eq.B4 and Eq.B5 (B7) Note that we have to pick the "difference" SPH differential form (e.g.see Tricco &Price 2012 andChan et al. 2021) if the speed of light changes discontinuously.Otherwise, the estimated gradient will be inaccurate, which would lead to numerical reflections at discontinuities.In the final step (C), we convert back to the local (ci) frames, and similarly for particle j.
The solution will depend on our choice for c0.We take c0 = max(ci, cj) for each pairwise interaction term to ensure symmetry and recover the rsl limit (i.e.satisfying the requirement (2)).In experiments, we found that the choice of c0 makes little difference30 as long as (I) c0 is between ci and cj and (II) c0 is symmetric with respect to i and j.
We also apply these changes to the artificial diffusion (see § 2.5 in Chan et al. 2021).

D ξ0,i Dt
and artificial viscosity: where D f is given by Eq. (B10) has to be solved with the "difference" SPH form too.
We demonstrate the robustness of our scheme in Fig. B2, which shows a case of radiation propagation in the optically thin limit.We consider a rapid change of c at x = 0 and investigate whether our scheme can work in the following aspects.
First, both rsl(s) and vsl(f) keep the product E rad c approximately constant within the plateau region, so they have the same ionizing power.Second, both the vsl(s) and vsl(f) profiles follow c specified in different regions.For example, the vsl(s) line reaches x = 0.5 at t = 1 (with c = 0.5).
To achieve this speed, the vsl(s) profile is compressed by a factor of two.Third, we have explicitly checked that photon conservation is better than 1% in all runs even with sudden changes in c.Finally, our scheme is also robust during the transition.Even with a rapid change in c at x = 0, the radiation profile at x = 0 is remarkably smooth.
The second test is the 1D photo-ionization test in Fig. B3.A constant photon flux is injected from the left and recombinations balance photo-ionizations.In the vsl runs, the reduced speed of light increases or decreases by a factor of ten at x = 0. We find that the equilibrium neutral fraction profiles computed with the vsl approximation match well that computed with the rsl approximation, demonstrating the accuracy of the vsl implementation.
The final test is a 3D isothermal Strömgren sphere (Iliev et al. 2006b Test 1) in Fig. B4.A source at the centre is emitting ionizing radiation at a constant rate of Ṅγ = 5 × 10 48 photons s −1 .The volume has a linear extent of 20 kpc, filled with 32 3 gas particles in a glasslike distribution.The gas is pure hydrogen with density nH = 10 −3 cm −3 .The collisional ionization coefficient is β = 3.1 × 10 −16 cm 3 s −1 , whereas the recombination coefficient is αB = 2.59 × 10 −13 cm 3 s −1 (with the on-thespot approximation).The photo-ionization cross section is σγHI = 8.13 × 10 −18 cm 2 .We apply the variable speed of light approximation such that c = 0.05c for r < 4 pkpc and c = 0.1c for r > 4 pkpc.Even with a sudden change of c by a factor of two, the equilibrium neutral fraction profile still agrees well with the analytic solution (from Chan et al. 2021).
Our vsl scheme improves upon the original scheme described by Katz et al. (2017) as follows.Our scheme allows for adaptive time stepping without the need for sub-cycling.Our scheme can also be applied to any moment-based radiative transfer scheme, e.g.adaptive mesh refinement, moving mesh, or any meshless method.

APPENDIX C: COSMOLOGICAL RADIATIVE TRANSFER
To account for cosmological expansion, we convert from proper to co-moving variables in the usual way, where primed variables are co-moving and a is the expansion factor.Furthermore, r is position, χ opacity, ρ gas density, W is the kernel and ∇ the gradient operator.Furthermore, we transform the radiation quantities as follows, Here, ξ(= E rad /ρ) is the ratio of the radiation energy density over the gas density f is the ratio of the radiative flux over the gas density and P is the radiation stress tensor (see table 1 in Chan et al. 2021).Our choice accounts only for the cosmological dilution of photons but not the redshifting of the wavelengths of photons.As a result, both gas and photon densities have the same dependency on scale factor, so ξ = E rad /ρ = E ′ rad /ρ ′ = ξ ′ .The primed radiation equations are then unchanged from the original proper equations, if we neglect small cosmological terms when the scale factor changes slowly.These are small corrections for the current application 31 and ignored in many previous works, e.g. in Shapiro & Giroux (1987); Rosdahl et al. (2013).
We verify the implementation of these equations by simulating the evolution of a cosmological Hii region in a homogeneous, constant-temperature, matter-dominated flat universe, with 10 per cent of its mass density in hydrogen (and neglecting helium).The hydrogen is initially fully neutral and is ionized by a source that switches on at zi = 9, emitting ionizing radiation at a constant rate of Ṅγ = 5 × 10 48 photons s −1 at hν = 13.6eV(and keeping the temperature constant).The collisional ionization coefficient is β = 3.1 × 10 −16 cm 3 s −1 , whereas the recombination coefficient is αB = 2.59 × 10 −13 cm 3 s −1 (at T = 10 4 K); making the on-the-spot approximation).The photo-ionization cross section is σγHI = 8.13 × 10 −18 cm 2 .
In the non-cosmological case, provided the density and temperature are uniform and constant, the I-front initially moves at the speed of light, before slowing down to the speed derived in Eq. ( 7) which eventually approaches zero due to 31 One of the neglected terms is the cosmological redshifting of the wavelength of photons.This is relevant for those photons that can travel over distances comparable to the Hubble distance, but not for the ionizing photons that have a mean-free path of only a few Mpc in our set-up.recombinations.Its position rI as a function of time is given by Here, ti is the time that the sources switches on, tr is the recombination time, t −1 r = αB nH, and rS = 3tr Ṅγ/(4πnH) is the Strömgren radius.In the cosmological case, the I-front initially also moves at the speed of light before slowing down as described by Eq. ( 7).However, the Strömgren radius increases as the density decreases due to the expansion of the Universe, Rs ∝ a 2 (provided the temperature is taken to be constant, rather than decreasing with the gas density).As a consequence, the ratio rI /rS decreases again at late times.The analytical solution is given by (Shapiro & Giroux 1987) where E2 is the exponential integral of order 2, and the constant λ is the ratio of the initial time over the initial recombination time, λ ≡ ti/tr(tI ).
We simulate the set-up using 64 3 gas particles in a computational volume with a linear extent of 40 ckpc, using c = 0.01 c, and compare the outcome to the analytical solution in Fig. C1.The Hii region is not well resolved initially, but at later times the agreement between the simulation and the analytical result is excellent.The good agreement at late times demonstrates that our choice of co-moving radiation variables works correctly.

APPENDIX D: NUMERICAL CONVERGENCE D1 Choice of clumping factor
In this section we contrast various definitions of the clumping factor, cr, that appear in the literature.For the first three of those, we only consider gas that is downstream of the I-front (as we did in the main text).
The first definition is c l from Emberson et al. ( 2013): (D1) asymptotic value of cr does not dependent on the choice of c, provided c ⩾ 0.05.

D3 Convergence with particle mass
We have performed simulations with different mass resolutions in the same cosmological volume to verify the level of numerical convergence, see Table 1.In Fig. D4, we plot the baryon fraction (in units of the cosmological mean value, f b ) as a function of halo mass, for simulations with different particle masses32 .In our simulations, the baryonic mass fraction of halos converges to within ∼ 30 per cent for halo mass M h ≳ 10 4.3 M⊙ provided that the dark matter particle mass mDM ∼ 100 M⊙ (200 dark matter particles per halo).This agrees well with the findings of Naoz et al. (2009), who found that ∼ 300 particles are needed to reach convergence at the 30 per cent level.Given our fiducial resolution of mDM ∼ 100 M⊙, our simulations resolve the baryon fraction of halos of mass Mmin to within 30 per cent.Here, Mmin is the minimum mass of halos that retain 50 per or more of their baryons, as discussed in §2.1.We examine the impact of finite mass resolution on the photo-evaporation of mini-halos in Fig. D5.At a given halo mass, mini-halos evaporate faster at lower resolution.This is partially because the central gas density is lower when the profile is not well-resolved, which enables the I-front to propagate deeper into the halo.This effect is less pronounced   The simulations differ significantly for halos with with mass M h < 10 6 M ⊙ , but these do not contribute much to the over all recombination rate.For halos more massive than M h = 10 6 M ⊙ , the lower resolution simulation yields more recombinations at higher density.This causes the over all recombination rate to be higher at lower resolution.
photo-evaporate, cr is determined by more massive minihalos halos which can be sufficiently resolved even when mDM = 800 M⊙.
We verify this scenario in Fig. D7 in which we plot the recombination rate as a function of over-density at the instance when cr is near its peak.Although the contribution of halos with mass M h < 10 6 M⊙ depends on resolution, they do not contribute significantly to the recombination rate at this stage of the evolution (This conclusion will depend on zi, since lower-mass halos contribute more at higher z.).The value of cr is dominated by halos with mass in the range 10 6−8 M⊙ halos, which can be resolved sufficiently well at a particle resolution of mDM ≲ 10 3 M⊙ (at least in the ionized outer regions)33 .At low Γ−12 and/or late times, the clumping factor is dominated by gas at over-densities ∆ < 100, and the recombination rate of this gas is resolved well when mDM < 10 3 M⊙.
However, mini-halos need to be resolved spatially as well.A mini-halo with M h = 10 6 M⊙ has a virial radius of ∼ 0.3 pkpc at z ∼ 8.For a quasi-Lagrangian code such as swift, the smoothing length of a gas particle is less than 0.05 kpc at the critical self-shielding density (nS ∼ 0.1 cm −3 ) and even smaller at higher density, when mDM ≲ 800 M⊙ (and when using equal numbers of dark matter and gas particles).Therefore spatial resolution is not a limiting factor in our calculations.
In contrast, spatial resolution may be challenging in mesh-based simulations, e.g uniform mesh methods might  2007), the vertical line shows M crit , the halo mass above which gas can be retained after reionization (Okamoto et al. 2008), as discussed in Section 2.
have difficulties resolving the central parts of mini-halos (see e.g. the convergence test in the appendix of D'Aloisio et al. 2020).Artificial diffusion can be problematic when the gas moves at high speed across a static mesh (Robertson et al. 2010).This can suppress the formation of structures that are not well resolved (Pontzen et al. 2021).Whether and how this affects the value of cr in mesh simulations may require further study.
In summary, a dark matter particle resolution of mDM ≲ 100 M⊙ is needed to model the photo-evaporation of minihalos (provided that the adaptive spatial resolution is high enough) 34 .This can be relaxed to mDM ≲ 10 3 M⊙ for studies of the clumping factor, because the lower-mass halos that are not well resolved at this coarser resolution contribute little to cr.

D4 Convergence with simulation volume
We quantify the impact of simulation volume on the halo mass function in Fig. D8.We consider three cubic volumes, S, M and L, with linear extents Lc ∼ 400, ∼ 800 and ∼ 1600 ckpc respectively.The numerical resolution is the same in these three runs.Not surprisingly, the number of halos with mass M h > 10 8 M⊙ is severely underestimated in simulation S. Such halos are significant photon sinks at the end of the eor and later, and clearly the simulation volume needs to be large enough to sample them correctly (see 34 Here we assume that the ihm gas follows the adiabatic cooling limit (Eq.1), which is a lower limit to the ihm temperature.In reality the ihm will be hotter than this lower limit due to, e.g., X-ray pre-heating as suggested by 21cm observations(HERA Collaboration et al. 2023).Thus, M min might be higher and the mass resolution requirement less stringent.also Appendix C of Cain et al. 2020).Our production runs are performed in simulations with Lc = 800 ckpc, (see Table 1) and therefore only begin to sample these more massive halos.
Fig. D9 compares the evolution of cr for these three simulations.The (first) value of c r,peak is similar for all three runs (see also Park et al. 2016).This first maximum depends on mini-halos with mass ≳ 10 6 M⊙, which are well sampled even in simulation S. The second peak is higher when Lc is smaller.The reason is that the ionization fronts that propagate from two opposing sides of the simulation cube overlap earlier for smaller values of Lc.At this earlier stage, cr is still higher, as mini-halos have had less time to photo-evaporate.The late time clumping factor is 50% higher in the simulation L compared to M , and a factor of two compared to S. This suggests that even our fiducial volume (M ) does not contain enough (massive) collapsed structures after photoevaporation to compute cr accurately.
We examine the reason for the relatively poor convergence with Lc in more detail in Fig. D10.The ratio of more massive to lower mass halos is larger in simulation L compared to M (Fig. D7).The gas in and surrounding these more massive halos contributes significantly to the recombination rate and hence cr in simulation L. This effect causes the dependence of cr post-reionization on Lc.
After reionization, the minimum halo mass that can contribute significantly to cr is close to the critical mass M char from Okamoto et al. (2008) (as discussed in Section 2.1, the minimum mass above which halos retain 50 per cent or more of their cosmological baryon fraction).This mass is around 10 8 M⊙ at z = 8 and Fig. D8 shows that only simulation L samples a reasonable number of such halos.This is not surprising: the parameters for our production runs were a trade-off between mass resolution and simulation volume for a given number of simulation particles and hence computation time.Our value of cr ≲ 3 post-reionization in simulation L agrees well with the value of cr ∼ 3 in the Lc = 2 h −1 cMpc simulation of D' Aloisio et al. (2020), and the range of cr = 2.4 − 2.9 found at z = 6 in the simulations with Lc = 40 h −1 cMpc by McQuinn et al. (2011).
Our convergence requirement on Lc is more stringent than that of Emberson et al. (2013), who obtained differences of less than 5 per cent in cr for simulations with Lc = 1 cMpc and L = 0.5 cMpc.We suggest that the difference stems from the fact that Emberson et al. (2013) did not account for photo-evaporation and hence are unduly biased to low-mass mini-halos, which dominate the recombination rate.Such low-mass halos are already sampled reasonably well in simulations with smaller Lc.
It would be interesting to extend the types of simulations performed here to larger volumes and investigate the impact of even more massive halos.Such simulations could then also be used to calculate the evolution of the attenuation length of ionizing photons.One aspect that such simulations should take into account is that these more massive halos may host a galaxy, and the feedback from the galaxy's stars may affect the distribution of neutral gas inside and outside the halo.

APPENDIX E: I-FRONT TRAPPING IN A MINIHALO
Here we test the ability of our rt implementation to trap an I-front in a halo by comparing to analytical solutions.These are based on the model presented in Section §4.3, except that here we assume that an initially planar I-front over runs a halo, rather than that the halo is illuminated by ionizing radiation from all directions.
Consider therefore an initially planar I-front, propagating along the Z axis from Z > 0 to Z < 0, and encountering a mini-halo located at the origin of the coordinate system (Z = 0).The mini-halo gas is modelled as a singular isothermal profile of Eq. ( 15).We assume that the gas den- sity outside the halo's virial radius R h is low enough such that the attenuation is negligible (but see §4.3).All of the gas is assumed to be static and has a uniform and constant temperature (hence, the recombination coefficient αB is also uniform and constant).
There are two ways to calculate the equilibrium neutral hydrogen profile.First, we consider the inverse Strömgren Layer (ISL) approximation under which gas downstream from the I-front is fully ionized and upstream is fully neutral.With this approximation, we can compute the total where the optical depth is given by We follow Altay & Theuns (2013) and take the logarithm on both side and differentiate with respect to Z, We then simplify and obtain the final differential equation with boundary condition x → x h when Z → Z h .We can integrate this equation numerically to find x(Z, b).This more accurate expression is actually quite similar to the approximate, analytical solution of Eq. (E3) (see the lower panel of Fig. E1).
We now perform a simulation with the following numerical set-up.The computational volume has linear extent of Lc = 6 pkpc, and is filled with particles such that the mean density is ⟨nH⟩out = 2 × 10 −4 cm −3 .The sphere is located at the centre of the simulation volume, has radius R h = 0.8 pkpc and its density is normalized by n H,h = (200/3)⟨nH⟩out so that its mean density is 200 times the surrounding density (and hence mimicking a cosmological halo).This density profile is realised with ≈ 2000 particles, this is equivalent to the resolution of a halo of mass 2.4 × 10 5 M⊙ in our main production runs (Table 1).We neglect hydrodynamics and keep the temperature of the gas constant and uniform, using a recombination coefficient of αB = 2.59 × 10 −13 cm 3 s −1 everywhere.We inject radiation from one side of the cubic simulation volume, with constant photon flux of F = 1.5 × 10 5 cm −2 s −1 , and use σHI = 8.13 × 10 −18 cm 2 .Finally, we impose that the optically thin direction of the Eddington tensor is parallel to the Z-axis (Chan et al. 2021) (which is perpendicular to the side through which we inject radiation), and use the reduced speed of light value of c = 0.01c.
The results of the simulation are compared to the analytical solution in Fig. E1, at 20 Myr after radiation was first injected.The upper panel is a slice through the centre showing the neutral fraction in the simulation, with the solid red line the approximate location of the inverse Strömgren layer from Eq. (E1).The lower panel compares the simulation results (black circle symbols) to the more accurate expression found by integrating Eq. (E10) numerically (black solid line).The agreement between the simulation and the analytical solution is within a few percentages, demonstrating the accuracy of the numerical scheme.This close agreement also suggests that the simulation can accurately account for I-front trapping in mini-halos with mass M h = 2.5×10 5 M⊙ and higher.To clarify better how regions of different densities react, we select particles in three ranges of over-density ∆ ≡ ρ/⟨ρ⟩ 35 Since the horizontal axis is log ∆ rather than ∆, the area under the curve is not proportional the fraction of volume at a given density.Rather, that would be dP/d log ∆ ∝ ∆ dP/d∆.at z = 7.84, and track them to z = 7.58.We plot their pdf's as a function of ∆, xHI and T in Fig. F1.

APPENDIX F: DETAILED RESPONSE TO THE PASSING OF IONISATION FRONTS
The initially lowest-density gas (∆ < 5, brown dashed line) changes little in volume density due to reionization but becomes mostly highly ionized, xHI ⩽ 10 −4 , and is photoheated to a temperature, T ∼ 2.2 × 10 4 K.A small fraction of this gas is only partially ionized and partially heated (T < 10 4 K).The initially highest density gas (∆ > 200, grey dot-dashed line) also changes relatively little, staying mostly dense, neutral and cold, as it self-shields from the I-front.A small fraction by mass becomes ionized and heated, expanding to lower densities.The passage of the I-front impacts predominantly gas initially at intermediate densities (5 < ∆ < 200, purple dotted line): this gas expands significantly to lower densities as it is ionized and photo-heated.
Some further features of Fig. F1 are worth noting.Although most of the gas is photo-heated to T ∼ 2.2 × 10 4 K, a fraction is hotter, reaching T ∼ 10 5 K.This corresponds to initially low-density gas (∆ < 5) that is first photo-heated by I-front and subsequently further adiabatically or shockheated by the expansion of a nearby filament.This can be seen in Fig. 5, where the hot gas appears as red lines delineating filamentary structures.This comparison also shows the origin of the slightly cooler gas of log T [K] ∼ 3.8: this is initially intermediate density gas in filaments that cools adiabatically as the filaments expand.Finally, there is a 'knee' of gas at temperature T ⩽ 10 3 K below which gas has not been ionized and photo-heated.This paper has been typeset from a T E X/ L A T E X file prepared by the author.5).We select gas at the higher redshift, z = 7.84, in three density ranges ∆ < 5, 5 < ∆ < 200 and ∆ > 200, as illustrated by coloured top patches in the left panel, and trace it to redshift z = 7.58.The probability distribution by volume of this gas is plotted as brown dashed, purple dotted and grey dash-dotted lines in all panels.

Figure 1 .
Figure1.Properties of halos with different masses as a function of redshift, the yellow vertical band indicates the eor, assumed to occur around z ∼ 6. M J is the Jeans mass (black solid line, Eq. 2), M F is the filtering mass (black dashed line, Eq. 3), M min is the mass above which halos contain 50% of the cosmic baryon fraction on average as computed from the simulations (red dot-dashed line), M ACH is the halo mass above which gas at the virial temperature can cool atomically (blue dashed line, Eq. 4).Halos with mass M min < M h < M ACH are mini-halos (left magenta region).After reionization, M char is the mass above which halos contain 50% of the cosmic baryon fraction on average as computed by Okamoto et al. (2008) (red solid line), and M crit is the mass above which halos can host a galaxy according to Benitez-Llambay & Frenk (2020) (green dotted line).We moved the M char line down since halos below M char can still hold a small fraction of gas: such halos can host RElHIcs (right magenta region, Benítez-Llambay et al. 2017).Halos in the grey region cannot retain a significant fraction of their baryons, whereas halos in cyan region could host a galaxy.Hence, the halos above the grey region are photon sinks, whereas those in cyan region host photon sources.The figure also shows that all mini-halos are photo-evaporated during and after the eor.
Figure1.Properties of halos with different masses as a function of redshift, the yellow vertical band indicates the eor, assumed to occur around z ∼ 6. M J is the Jeans mass (black solid line, Eq. 2), M F is the filtering mass (black dashed line, Eq. 3), M min is the mass above which halos contain 50% of the cosmic baryon fraction on average as computed from the simulations (red dot-dashed line), M ACH is the halo mass above which gas at the virial temperature can cool atomically (blue dashed line, Eq. 4).Halos with mass M min < M h < M ACH are mini-halos (left magenta region).After reionization, M char is the mass above which halos contain 50% of the cosmic baryon fraction on average as computed by Okamoto et al. (2008) (red solid line), and M crit is the mass above which halos can host a galaxy according to Benitez-Llambay & Frenk (2020) (green dotted line).We moved the M char line down since halos below M char can still hold a small fraction of gas: such halos can host RElHIcs (right magenta region, Benítez-Llambay et al. 2017).Halos in the grey region cannot retain a significant fraction of their baryons, whereas halos in cyan region could host a galaxy.Hence, the halos above the grey region are photon sinks, whereas those in cyan region host photon sources.The figure also shows that all mini-halos are photo-evaporated during and after the eor.

Figure 2 .−and
Figure 2. The gas fraction inside halos at different redshifts before reionization in our highest resolution simulation (S512G00; see §3 for details), in units of the cosmic baryon fraction (f b ≡ Ω b /Ωm).Vertical dotted lines indicate M min , the mass above which the baryon fraction of a halo is f b /2 or more on average (differently coloured lines from left to right, for z = 7.4, 9.6 and 19.4).Vertical dashed lines indicate the halo mass above which gas at the virial temperature can cool atomically, M ACH (from right to left, for z = 7.4, 9.6 and 19.4).Horizontal lines indicate a baryon fraction of f b /2 (dotted) and f b (dashed).Mini-halos are halos with masses between M ACH (vertical dashed) and M min (vertical dotted).

Figure 3 .
Figure3.A schematic picture of the three stages of reionization, in which an initially planar ionization front (red dashed line) moves upwards in the diagram.In the lower part of the diagram, the ihm is highly ionized (black region) except in halos with mass above a characteristic mass, M char .The outskirts of the gas in these halos are highly ionized but the halos' potential wells are deep enough that the gas cannot photo-evaporate: these are Lyman-limit systems.The inner parts of the halos can selfshield, and the gas is neutral; these are damped Lyman-α systems (yellow, solid discs).In the middle part of the diagram, the ionization front has passed only recently, and the gas in mini-halos is still photo-evaporating (hollow circles).The ionization front has not yet reached the top of the diagram, where the ihm is still neutral (yellow region).Yellow vertical stripes are regions where gas remains neutral due to the shadowing of the ionizing radiation by an intervening self-shielded region.The I-front is slowing down and becomes increasingly corrugated rather than planar as a result of the density inhomogeneities in the ihm.

Figure 4 .
Figure 4. Symbols: halo mass function of the fiducial simulation (M512) at three different redshifts, lines show the fitting functions fromReed et al. (2007).Colours correspond to z = 10.2 (blue squares, and solid line), z = 7.9 (green circle and dotted line ), and z = 6.0 (red triangles and dashed line).Empty symbols show mass bins with fewer than 5 halos per decade in halo mass.The thin vertical lines on the right show the minimum mass of atomic cooling halos (Eq.4); halos in the grey region are not well resolved, containing fewer than 300 dark matter particles.Our fiducial simulation resolves mini-halos with more than 300 particles in a volume sufficiently large to contain a few halos in which gas can cool atomically.

Figure 5 .
Figure5.Impact and response of gas structures to the passage of ionization fronts, entering the computational volume simultaneously from the top and the bottom.Individual panels visualize slices/projections of the M512z8G03 simulation box at z = 7.84 (upper panels) and z = 7.58 (lower panels).Radiation is injected at redshift z i ≈ 8. Panels from left to right are hydrogen gas density (n H ) (Mach number in the lower left), the neutral fraction (x HI ), and the gas temperature (T ).Panels labelled 'projected' correspond to projections of the computational volume; the others correspond to infinitesimally thin slices through the mid-plane.As the I-front propagates through the volume, gas becomes ionized and photo-heated, and mini-halos photo-evaporate.Photo-heating causes the expansion of the filaments.

Figure 7 .
Figure 7. Temperature vs over-density diagrams for the M512z8G03 run at the start of reionization, z ≲ z i , (upper panel) and when the clumping factor cr is maximal (lower panel).The colour scale encodes the mass-weighted neutral fraction.The red dashed line is the temperature where the cooling due to HI is maximum.The blue dotted line is the initial entropy (Eq.13).Contours indicate the surface density of gas particles in this plot, evenly spaced in log surface density.

Figure 8 .
Figure 8. Slices through the M512z8G03 run at z = 7.58, after the I-front has crossed the simulation volume (see Fig.5).Panels show temperature (upper panel), velocity in the x-direction (central panel), and hydrogen density (lower panel).These particular slices illustrate the photo-evaporation of two filaments, indicated by dashed black lines in each panel.At this moment in time, the density in the filaments is higher than in their surroundings (lower panel), yet the gas is cooler due to adiabatic expansion.The gas expands at a speed close to the sound speed (central panel) and compresses and heats the gas in the surrounding ihm as seen in the upper panel.This type of compressional heating is the origin of the hotter gas, T > 10 4.2 K, seen in Fig.7.

Figure 9 .
Figure 9. Neutral fraction as a function of radius (lower horizontal axis) or density (upper horizontal axis) in halos with a spherical density profile of the form n H= n H,h (R h /R) 2, where R h is the virial radius, for different halo masses (different line-styles) and redshifts (colours) as per the legend, obtained by integrating Eq. (22).The horizontal dot-dashed line and solid line indicates x = 1/2 and x = 1.The dotted line is the approximation for x(R) from Eq. (25).Small dots and small diamonds indicate the value where x = 1/2 and x = 1 in the approximate model.The calculations take Γ h = 10 −12 s for the photo-ionization rate at the virial radius, αr = α B (T = 10 4 K) for the recombination rate, and σ HI = 1.62 × 10 −18 cm −2 for the photo-ionization cross section.The approximate model is fairly accurate, even in predicting the location and density where x = 1/2 and x = 1.The value of n H,h = 8.6 × 10 −3 and 4.1 × 10 −3 cm −3 at z = 8 and z = 6.

Figure 11 .
Figure 11.Hydrogen density, n S , where gas transitions to neutral, x > 1/2, in the simulations plotted against the volumeaveraged photo-ionization rate, Γ −12 , after the I-front has crossed the simulation volume.Large markers refer to gas in halos that contain more than 10 per cent of the cosmic baryon fraction, f b ; small markers are for halos containing less than 10 per cent of f b .Colours and markers encode different values of z i , as per the legend.The dashed line indicates the theoretical scaling (n S ∝ Γ 2/3 −12 ; Eq. 28).

Figure 13 .
Figure 13.Ratio between gas mass within the virial radius and the halo mass, in units of the cosmic baryon fraction, f b .Solid lines refer to total gas mass, dashed lines to neutral gas.Different colours correspond to different times, t dur , since the halo was overrun by I-front, from blue (recently) to purple (long ago).The upper panel has Γ −12 = 0.3, whereas the lower panel has the lower photo-ionization rate of Γ −12 = 0.03.Halos in the left shaded regions are resolved with fewer than three hundred dark matter particles.Halos in the right shaded region are susceptible to atomic cooling.

Fig. 13 .
We plot the halo baryon fraction in units of the cosmic mean, Mgas/(f b M h ) (where f b ≡ Ω b /Ωm), at various values of t dur .More massive halos can hold onto their gas for longer times.A M h = 10 5 M⊙ halo loses half of its gas in t dur = 20 Myr for Γ−12 = 0.3.On the other hand, a halo with M h = 10 7 M⊙ has only lost ≈ 20 per cent of its gas at t dur = 50 Myr at Γ−12 = 0.3.A lower Γ−12 can also slow down the photo-evaporation, e.g. the M h = 10 5 M⊙ halo takes twice as long to photo-evaporate at Γ−12 = 0.03 (than Γ−12 = 0.3).

Figure 14 .
Figure14.Evaporation time, tev, as a function of halo mass, M h .The upper panels show the time to lose 90 per cent of the total gas mass (solid lines) or the neutral gas mass (dashed lines) (relative to the cosmic share of baryons; after the halo was overrun by I-front).Thin dotted lines indicate that evaporation is not yet complete when the simulation is stopped.The lower panels are similar, except that halos must lose 99.9 per cent of their gas.Triangles are the simulation results ofIliev et al. (2005a).Our evaporation times are slightly longer, but the trends with M h and Γ −12 are very similar: the evaporation time increases with M h and z i , and decreases with increasing Γ −12 .

Figure 15 .
Figure 15.Clumping factor, cr (Eq.33), as a function of redshift, z, for simulations with different photo-ionization rates Γ −12 (solid, dashed, dot-dashed lines correspond to Γ −12 = 0.3, 0.15, and 0.03).The line colour indicates different redshifts of ionization (yellow, red and blue correspond to z i ∼10, 8 and 6).When the I-front enters the simulation volume, cr rapidly increases to a peak value, c r,peak .Then, cr decreases on a time scale of order 100 Myrs before reaching asymptotically an approximately constant value.The value of c r,peak increases with increasing Γ −12 and decreasing z i , but the asymptotic late-time value of cr is approximately independent of Γ −12 .
Fig. 7 of D'Aloisio et al. (2020), when comparing the runs with zi = 8, 6 and Γ−12 = 0.3.However, D'Aloisio et al. (2020) used the clumping factor definition of Eq. (D3), which is tens of per cent lower than our fiducial definition Eq. (33) (because Eq.D3 used a lower reference temperature for recombination rate; see Appendix D).Thus, our results are roughly consistent.Our simulations are not directly comparable to those of Park et al. (2016), since we used Γ−12 = 0.3, roughly a

Figure 17 .
Figure 17.Lower panel: the cumulative number of recombinations per hydrogen atom, Nrec/N H , as a function of time since ionizing photons were injected in the simulation volume.The thick lines represent the simulation results.The thin lines are computed for a homogeneous medium with the same median temperature.The ratio of these is plotted in the upper panel.Line colours and styles are as in Fig. 15.Although the clumping factor can be quite high, cr = 5 − 20 (upper panel and Fig. 15), recombinations increase the required number of ionizing photons per hydrogen to complete reionization by a much smaller fraction, 1.2-1.7,because mini-halos photo-evaporate quickly.This fraction increases with Γ −12 (different line styles) but is only weakly dependent on z i (different colours).

Figure B1 .
Figure B1.Schematic diagram illustrating our variable-speedof-light implementation.The left and right lower panels represent particles i and j with different specific radiation energies ξ and reduced speeds of light c respectively.To evaluate the energy and flux exchanges, particles i and j are transformed to a frame with a reduced speed of light c0 (upper panels).After evaluating the particle-particle interaction, radiation variables are transformed back to the original frame.

Figure B3 .
Figure B3.Neutral gas fraction, x HI = n HI /n H , in a 1D photoionization test of the vsl approximation.We use dimensionless variables, in which the simulation volume contains 100 gas particles with hydrogen density n H = 1.0, photo-ionization crosssection σ HI = 10 2 and recombination coefficient α B = 4.0.Radiation is injected at the location x = −8 with a constant radiation flux Fγ = 60.Different curves show the run of x HI with x for simulations with different values of rsl with cr the value of rsl at x > 0 in units of the rsl at x < 0. Small horizontal lines indicate the local value of the smoothing length.The profile is shown at time t = 10 2 when x HI is in equilibrium.The test shows that the equilibrium value does not depend on the choice of the rsl.

Figure B4 .
Figure B4.Isothermal Strömgren sphere from Iliev et al. (2006b) Test 1: a source of radiation, located at radius r = 0, photoionizes hydrogen gas, kept at constant density and constant temperature.The system is shown a time t = 500 Myrs after the source is switched on.Red points are the neutral hydrogen fraction, x HI = n HI /n H , of individual gas particles, with the thick red points showing binned values, with horizontal error bars indicating the bin width, and vertical error bars the standard deviation; the black line is the analytical solution.Small blue points, and thick blue points with error bars show the corresponding ionized fraction, 1 − x HI , with the yellow line the analytical solution.The vertical light blue line is the approximate analytic location of the Strömgren radius.In this vsl test, c increases by a factor two at r = 4 kpc.This has little effect on the numerical solution, which is almost identical to the case where c is constant (Fig. B4 in Chan et al. 2021).

Figure
Figure C1.I-front location, r I , divided by the Strömgren radius, r S , as a function of redshift; r I is operationally defined as the location where 50 per cent of the gas is ionized.The source switches on at redshift z i = 9.Magenta squares show the simulation result, with error-bars depicting the gas smoothing lengths divided by r S .The black line is the analytic solution from Shapiro & Giroux (1987).

Figure D1 .
Figure D1.Evolution of the clumping factor in simulation M512z8G03, where t = 0 is the instance the I-front enters the simulation volume.Different curves correspond to different definition of the clumping factor, as described in section D1.

Figure D4 .
Figure D4.Lower panel: Baryon fraction of halos in units of the cosmic mean, f b , for simulations with different particles masses (colours).Upper panel: the baryon mass divided by the baryon mass of the highest resolution simulation.In both panels, thick (thin) lines refer to halos resolved with more than (less than) 300 DM particles.Results are shown at redshift z = 8, before radiation was injected in the simulations.The baryon fraction is converged to within 20-30 per cent, provided the halo is resolved with more than ∼ 300 particles.

Figure D5 .
Figure D5.Ratio between gas and halo mass, in units of the cosmic baryon fraction, f b .It is similar to Fig. 4, but for simulations with different particles masses (plotted using different line styles).Different colours correspond to different times, t dur , since the halo was over run by I-front: blue: 1 Myr; green: 20 Myr; brown: 100 Myr.The simulation volume is 800 ckpc on a side.Higher resolution simulations yield longer photo-evaporation time scales at a given halo mass.

Figure D6 .
Figure D6.As Fig. D3, but for simulations with different particles masses.The simulation volume is 800 ckpc on a side, and we use Γ −12 = 0.3, and z i = 8.The evolution of the clumping factor differs by less than 5 per cent between these runs.

Figure D7 .
Figure D7.As Fig. 16, showing the contribution to the recombination rate as a function of over-density with coloured lines depicting the net contribution of halos of different mass.Results are shown at the instance that cr = c r,peak .The thick lines are for simulation M512z8G03 (particle mass m DM = 120 M ⊙ ), the thin lines are simulation M128z8G03 (with m DM = 9600 M ⊙ ).The simulations differ significantly for halos with with mass M h < 10 6 M ⊙ , but these do not contribute much to the over all recombination rate.For halos more massive than M h = 10 6 M ⊙ , the lower resolution simulation yields more recombinations at higher density.This causes the over all recombination rate to be higher at lower resolution.

Figure D8 .
FigureD8.Impact of the extent of the simulation volume on the halo mass function at redshift 7.9.Colours and symbols correspond to simulations performed in volumes of linear extent ranging from 400 to 1600 ckpc, as per the legend.The empty symbols show bins with fewer than 5 halos per log 10 (M ).The diagonal line is the fit to the mass function fromReed et al. (2007), the vertical line shows M crit , the halo mass above which gas can be retained after reionization(Okamoto et al. 2008), as discussed in Section 2.

Figure D9 .
Figure D9.Convergence of the clumping factor for different linear extents of the simulation volume.The mass resolution is m DM = 960 M ⊙ , and we assumed Γ −12 = 0.3, and z i = 8.The value of cr at its peak only depends weakly on volume, but the asymptotic value of cr increases from ∼ 2 to ∼ 3 when increasing the linear extent from 800 to 1600 pkpc.The figure also suggests that the asymptotic value of cr is underestimated, even in the largest volume.

Figure D10 .
Figure D10.As Fig. D7, but for simulations with different linear extent, Lc: Lc = 800 ckpc (simulation S, thin lines) and Lc = 1600 ckpc (simulation L, thick lines).The simulation particle mass is m DM = 800 M ⊙ for both.Simulation L contains many more halos with mass 10 8 M ⊙ compared to simulation M .Gas in those halos increases the recombination rate at ∆ ∼ 10 2 significantly, as can be seen by comparing the thick and thin green lines.

Figure E1 .
Figure E1.Ionization front trapping by a static gas cloud with a 1/R 2 density distribution: the initially planar ionization front moves vertically down.Upper panel: slice through the centre of the sphere, with colours indicating neutral fraction.The dashed blue line is the radius of the sphere, the solid red line is the approximate location of the inverse Strömgren layer (isl) from Eq. (E1).The arrow indicates the origin and the positive Z direction.Lower panel: neutral fraction x HI as a function of Z, where the Z axis runs vertically up from the centre of the sphere (as indicated in the upper panel).The black circles and cyan triangles are obtained from the simulation with impact parameters b = 0 pkpc and b = 0.35 pkpc respectively.The black solid line and cyan dotted line are obtained by numerically integrating Eq. (E10) with impact parameters b = 0 pkpc and b = 0.35 pkpc respectively.The vertical dashed black line is the location of the isl from Eq. (E1) for an impact parameters b = 0 pkpc.Finally, the black dashdot line represents the approximate neutral fraction (Eq.E3).The simulation follows the analytical result very closely.The text describes the simulation set-up in detail.

F
recombination rate, R, along a line parallel to the Z axis at impact parameter b.Equating R to the flux incoming flux yields the location ZISL(b) of the inverse Strömgren layer, where gas transits from ionized to neutral.ZISL(b) can be solved through the following equation,Z h = (R 2 h − b 2 ) 1/2.The sign of ZISL can be positive or negative.Ignoring the optical depth of the neutral gas downstream, we can find an approximate expression for the neutral fraction x ≡ nHI/nH in photo-ionization equilibrium, x h is the neutral fraction at R h .Second, we can drop the ISL approximation and calculate the neutral gas distribution considering the optical depth.Assuming photo-ionization equilibrium, we calculate the photon number density nγ:∂nγ ∂t = − ∂fγ ∂Z − nHIcσHInγ = 0 , (E4)where fγ is the photon flux which is approximately fγ = nγc if we ignore scattering.Then we solve nγ by integration nγ

Fig
Fig. F1 shows the response of the ihm to the passing of I-front for the simulation M512z8G03 quantitatively.The probability density distribution by volume, PV , is plotted at z = 7.84 (blue line) and z = 7.58 (orange line), as a function of 35 log ∆, log xHI, and log T .Photo-ionization reduces PV at intermediate densities, log ∆ ∼ 1 − 3, due to photo-evaporation of the outskirts of mini-halos.To clarify better how regions of different densities react, we select particles in three ranges of over-density ∆ ≡ ρ/⟨ρ⟩

Figure F1 .
Figure F1.Probability distribution by volume, P V , of the gas over-density, ∆(= ρ/⟨ρ⟩), neutral fraction, x HI , and gas temperature, T (left to right), for two different redshifts as labelled, for M512z8G03 (the simulation shown in Fig.5).We select gas at the higher redshift, z = 7.84, in three density ranges ∆ < 5, 5 < ∆ < 200 and ∆ > 200, as illustrated by coloured top patches in the left panel, and trace it to redshift z = 7.58.The probability distribution by volume of this gas is plotted as brown dashed, purple dotted and grey dash-dotted lines in all panels.

Table 1 .
Details of the simulation suite.From left to right: simulation identifier, linear extent of the simulated volume L box , redshift of reionization z i , photo-ionization rate in units of 10 −12 s −1 , Γ −12 , number of gas and dark matter particles, N , gas particle masses mgas, dark matter particle masses m DM , gravitational softening lengths l soft , reduced speed of light at the mean density, c, and mean hydrogen density, ⟨n H ⟩, at z i .The simulation identifier encodes the extent of the simulated volume, the number of particles, the value of z i and the value of Γ −12 .