Reionization with Simba: How much does astrophysics matter in modeling cosmic reionization?

Traditional large-scale models of reionization usually employ simple deterministic relations between halo mass and luminosity to predict how reionization proceeds. We here examine the impact on modelling reionization of using more detailed models for the ionizing sources as identified within the $100~{\rm Mpc/h}$ cosmological hydrodynamic simulation {\sc Simba}, coupled with post-processed radiative transfer. Comparing with simple (one-to-one) models, the main difference of using {\sc Simba} sources is the scatter in the relation between dark matter halos and star formation, and hence ionizing emissivity. We find that, at the power spectrum level, the ionization morphology remains mostly unchanged, regardless of the variability in the number of sources or escape fraction. In particular, the power spectrum shape remains unaffected and its amplitude changes slightly by less than 5-10\%, throughout reionization, depending on the scale and neutral fraction. Our results show that simplified models of ionizing sources remain viable to efficiently model the structure of reionization on cosmological scales, although the precise progress of reionization requires accounting for the scatter induced by astrophysical effects.


INTRODUCTION
Many upcoming surveys, such as the Square Kilometer Array (SKA, Mellema et al. 2013), the Hydrogen Epoch of Reionization Array (HERA, DeBoer et al. 2017), the Low Frequency Array (LOFAR, van Haarlem et al. 2013), are devoted to detecting reionization in the near future via its impact on HI 21cm emission topology. On the other hand, the James Webb Space Telescope (JWST, Gardner et al. 2006) the Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST, Ivezić et al. 2019), and Nancy Grace Roman Space Telescope (Roman, Spergel et al. 2015), are devoted to de-tecting the primary sources driving reionization at highredshift (galaxies and quasars). With these growing observational efforts, we are moving towards an era of detailed observations of high-redshift galaxy formation. It is therefore important to test and quantify the uncertainties in our current models of reionization, by exploring the interplay between galaxy formation and cosmology, in order to prepare for extracting the maximum amount of information from the future surveys.
The main sources of uncertainties in reionization models arise from the ionizing source population (Stark 2016), the clumpiness of the ionized gas (which sets the number of recombinations; D' Aloisio et al. 2020), and approximations in performing the radiation transport (RT) (Wu et al. 2021). The impact of RT has been investigated in detail in several studies (e.g. Ma-jumdar et al. 2014;Molaro et al. 2019;Zahn et al. 2011;Mesinger et al. 2011), who found that, using the same ionizing sources, the Excursion-Set Formalism (Press & Schechter 1974;Bond et al. 1991), which is an algorithm to identify maximally ionized bubbles within spherical regions around sources, produces similar topology to the post-processing radiative transfer of cosmological hydrodynamic simulations. Hence, it appears justified to use excursion-set based models to study the reionization signal on large scales. Meanwhile, the impact of including inhomogeneous recombination and clumping effects -the sink model -has also been studied in several works (e.g. Sobacchi & Mesinger 2014;Hassan et al. 2016), which concluded that inhomogeneous recombination reduces the ionized region size, delays reionization, and suppresses the large scale 21 cm power. In contrast, the large scale impact of including more physically motivated recipes to model the ionizing source population has not been studied. This is the goal of this work.
In Hassan et al. (2016), we took a first step to study the impact of modelling ionizing sources more accurately within semi-numerical simulations, and showed that including power-law mass-dependent ionizing emissivity motivated from hydrodynamical simulations as opposed to the usual linear dependence on mass (through an efficiency parameter ζ ≡ R ion /M h of ionizing photon emissivity R ion vs. halo mass M h ) increases the duration of reionization and boosts the small scale 21cm power spectrum. This indicates that a modest change in the ionizing emissivity-halo mass (R ion -M h ) relation has a major impact on reionization globally. This impact owes to the more massive sources in the super-linear scaling being more clustered and hence producing larger HII regions (Furlanetto et al. 2006;McQuinn et al. 2007) Currently, the most accurate way to model ionizing sources in a cosmological volume is through state-ofthe-art cosmological hydrodynamical galaxy formation simulations, in which star formation is modelled following, for instance, the sub-grid multiphase model of Springel & Hernquist (2003) or the H 2 -regulated model of Krumholz et al. (2009a). Such simulations include Illustris (Genel et al. 2014), IllustrisTNG (Pillepich et al. 2018), Eagle (Schaye et al. 2015), Horizon-AGN (Kaviraj et al. 2017 Ocvirk et al. 2020), SPHINX (Katz et al. 2021), THE-SAN (Kannan et al. 2022) using the AREPO-RT (Kannan et al. 2019), to name a few. It has recently become possible to run galaxy formation simulations on cosmological scales ( 100 Mpc) while resolving halos down to M h < 10 10 M , thanks to advances in computing capa-bilities. Hence, it is now possible to ask the question: Does detailed astrophysics matter for modeling cosmic reionization?
The main difference in using sources identified within these galaxy formation simulations consists in including the scatter around the R ion -M h relation. This scatter, which is due to stellar feedback, merger history and/or photo escape fraction, might have a direct correlation with the density field and/or other properties such as metallicity, and hence its impact on large scale reionization remains non-trivial. As a simple example, by assuming a lognormal scatter around sources identified within semi-numerical excursion-set based models, Whitler et al. (2020) showed that the scatter reduces the Lyα visibility for brightest galaxies.
In this work, we use ionizing sources identified within the state-of-the-art Simba simulation, that has been shown to reproduce a wide range of observations, at both low-z (see Davé et al. 2019) and during the reionization epoch (Wu et al. 2020;Leung et al. 2020). We consider a broad range of source models that uses the SIMBA star formation rates as inputs, but considering different minimum halo masses for sources owing to resolution limits in the simulations and varied escape fraction pictures. We then run large scale radiative transfer reionization simulations, with the radiative transfer done in post-processing using ATON (Aubert & Teyssier 2008, 2010) via a moment-based method with an M1 closure relation. We compare these results with a simple powerlaw model of the ionizing emissivity as a function of halo mass, as commonly used in large-scale semi-numerical models. Compared to previous works (e.g. Whitler et al. 2020), our approach has several advantages and improvements, namely the scatter in ionization rate naturally emerges from the detailed sub-grid physics implemented within Simba, and the impact on large scale topology is determined through accurate modelling of radiative transfer with ATON.
This paper is organized as follows: we briefly describe the Simba simulation in § 2.1, introduce the source models in § 2.2, explain how the reionization realization is generated in § 2.3, present our results in § 3, and conclude in § 4. Throughout this work, we adopt a ΛCDM cosmology with parameters Ω m = 0.3, Ω Λ = 0.7, Ω b = 0.048, h = 0.68, σ 8 = 0.82, and n s = 0.97, consistent with Planck Collaboration et al. (2016) constraints. All quantities are in comoving units unless otherwise noted.

Simba Simulations
We here briefly describe Simba and refer the reader to Davé et al. (2019) for more detailed information about the physics implemented in the simulation.
The Simba model was introduced in Davé et al. (2019). Simba is a follow-up to the Mufasa ) cosmological galaxy formation simulation using Gizmo's meshless finite mass hydrodynamics solver (Hopkins 2015(Hopkins , 2017. Radiative cooling and photo-ionisation heating are implemented using the updated Grackle−3.1 library (Smith et al. 2017). A spatially-uniform ionizing background from Haardt & Madau (2012) is assumed, in which self-shielding is accounted for following Rahmati et al. (2013). The chemical enrichment model tracks nine elements (C,N,O,Ne,Mg,Si,S,Ca,Fe) arising from Type II supernovae (SNe), Type Ia SNe, and Asymptotic Giant Branch (AGB) stars. Type Ia SNe and AGB wind heating are also included. The star formation-driven galactic winds are kinetically launched and decoupled into hot and cold phase winds, and the winds are metalloaded owing to local enrichment from supernovae, with overall metal mass being conserved. The mass rate entering galactic outflows is modelled with a broken power law following Anglés-Alcázar et al. (2017b). The quasi-linear scaling of wind velocity with escape velocity from Muratov et al. (2015) is adopted. Simba further implements black hole growth via a torque-limited accretion model (Anglés-Alcázar et al. 2017a), and the feedback from Active Galactic Nuclei (AGN) is implemented in both radiative and jet modes. Simba also accounts for the X-ray AGN feedback in the surrounding gas following Choi et al. (2012). Dust production and destruction are modeled on-the-fly, leading to predictions of dust abundance that broadly agree with a variety of constraints over cosmic time (Li et al. 2019).
Simba employs a molecular gas-based prescription following Krumholz et al. (2009b, hereafter KMT) to form stars. KMT is a physically motivated recipe to model star formation as seen in local disk galaxies, where a strong correlation is seen between SFR and molecular gas content (e.g. Leroy et al. 2008). Star formation is only allowed in the dense gas phase (ISM gas) above a hydrogen number density n H ≥ 0.13 cm −3 , though in practice the H 2 fraction forces star formation to occur at much higher densities (n H 1 cm −3 ). Simba assumes a Chabrier (2003) initial mass function (IMF) throughout.
Simba is an attractive platform to use in this study due its ability to match a wide range of observations as show in Davé et al. (2019). These observations include the stellar mass function evolution, main sequence, evolution in gas fractions (both neutral and molecular), the gas-phase and stellar mass-metallicity relations, galaxy photometric projected sizes, hot halo gas fractions, black hole mass-stellar mass relation, and the dust mass function. On the other hand, Simba fails to produce a sharp knee in the steller mass function at z=0, and produces larger low-mass quenched galaxy sizes, to name few examples. Despite these small disagreements, Simba remains a state-of-the-art platform for studying the scatter impact on reionization morphology.

Source Models
The ionization rate (R ion ) from each galaxy is computed using a stellar metallicity-dependent parameterization of the ionizing photon flux (Q(Z)), which reproduces the equilibrium values compiled by Schaerer (2003). This Q(Z) parameterization is provided by Finlator et al. (2011a): where Q is in units of s −1 (M yr −1 ) −1 and the last term converts to a Chabrier (2003) IMF. It is straightforward to convert the Q to R ion using the SFR as follows: We directly compute the ionizing emissivity from the SIMBA halo catalogs. We use the stellar mass-weighted metallicity Z of all star particles and the total SFR of all gas particles within each halo to evaluate Equation 2 in order to provide the emissivity required for the radiative transfer calculations.
As mentioned earlier, our aim is to use these realistic sources, as identified within Simba, to study the impact on reionization on cosmological scales, as compared to simplified models. We use the largest Simba run available which has a volume of 100h −1 Mpc with 1024 3 dark matter and gas particles each. While this simulation has a sufficiently large volume suitable to study large scale reionization (Iliev et al. 2014), the dark matter particle mass resolution is about 9.6×10 7 M , and hence the minimum halo mass resolved would be of order 10 9 M assuming 16 particles as minimum to identify halos in the friends-of-friends finder. However, if the resolution is high enough to resolve low mass halos down to 10 8 M , the number of sources increases, and in this case the scatter impact would be even smaller since the ionizing emissivity average in pixels is now taken over a larger number of sources. Later we include models with different number of sources to address this point.
To construct a simple model of the average R ion as a function of halo mass -a comparison model that averages out the fluctuations of individual halos -we follow our earlier approach presented in Hassan et al. (2016) where we obtained a parameterization for the ionizing emissivity R ion as a function of halo mass M h and redshift z from hydrodynamic simulations. We here repeat the same exercise to obtain the best fitting parameters for the following simple model of R ion : where A, B, C and D are free parameters. This fitting function, which is similar to the Schechter function, is motivated by the fact the R ion behaves as a power-law at high halo masses and exhibits a turnover at the low mass end. In Hassan et al. (2016), we found the turnover occurs around 10 8 M (see Figure 1 therein). Since the largest Simba run used in this study resolves halos down to ∼ 10 9−10 M , we refer to models using this simple relation (Equation 3) as power-law.
In order to estimate the slope of R ion more accurately from Simba, we include a higher resolution run which has a smaller volume of 25h −1 Mpc with the same number of particles (2 × 1024 3 ). Hence, this run resolves halos down to 10 8.5 M ; halos below this mass are not expected to contribute significantly (Finlator et al. 2011b). We combine these two simulations at different redshifts (z = 6 − 10) to obtain the R ion best fitting parameters as shown in Figure 1. Conservatively, one may obtain the best-fit parameters using catalogs from high (e.g. z=20) to low (e.g. z=6) redshifts. However, restricting the fitting to redshifts lower than z=10 would not impact the results since reionization, in models considered in this study, begins only afterward. This figure shows R ion as a function of M h at different z from Simba m100n1024 (light green) and m25n1024 (dark green) runs using Equation 2. The fitting function, Equation 3, is shown in red, and the best fitting parameters are found using a non-linear least square fit as follows: A = 3.51 × 10 39 ± 5.6 × 10 37 [sec −1 ], B = 3.3 × 10 8 ± 2.6 × 10 6 [M ], C = 0.234 ± 0.00078 and D = 3.2 ± 0.0082. These parameters are similar to those found earlier by Hassan et al. (2016), based on full radiative hydrodynamic simulations (with much smaller volume). It is evident that these parameters produce a good fit to the ionizing emissivity at different redshifts, and hence we will fit this as our power-law model of R ion to compare with using a more realistic R ion directly from Simba.
To perform a consistent comparison and to carefully investigate the scatter impact, we tune the photon escape fraction so that all models would produce roughly the same number of ionizing photons. In addition, we ensure that all models do produce the expected ionizing emissivity evolution (e.g. Becker & Bolton 2013) as well as the neutral fraction measurements by the end of reionization (e.g. Fan et al. 2006;Bouwens et al. 2015). It is worth noting that, for the purpose of this study, it is not important to match or reproduce any measurements. The most important part of this analysis is to ensure that the comparison is performed at the same redshift and neutral fraction in order to exclude any other cosmological and astrophysical effects.
Using the same density field and recombination rate, we consider five different source models, that are tuned to emit similar amount of ionizing photons and obtain same reionization history as follows: • Simba: uses all sources with emissivities computed using Equation 2, and f esc = 0.25 (blue in Figure 2). • Simba-random: uses only 25% of all sources, randomly selected within different halo mass bin of 0.25 dex. Similar to Simba, the emissivity is computed using Equation 2, but with f esc = 1.0 (cyan in Figure 2).
• Simba-brightest: uses only the brightest halos (the most luminous halos) in each mass bin of 0.25 dex (halos with the highest R ion ). These halos are about 7% out of all halos. Similar to Simba-random in terms of emissivity and unity f esc (dark blue in Figure 2).
• power-law: uses all sources with emissivities computed using Equation 3, and f esc = 0.25 (red line in Figure 2).
• power-law with scatter: Same as power-law, except with the addition of scatter assuming that R ion may be represented by a lognormal distribution with a mean given by Equation 3 and standard deviation set to 0.3 dex 1 . This model uses f esc = 0.2 (green in Figure 2).
A visual summary of these different models is provided in Figure 2, where solid lines represent the 50th percentile of the different R ion distributions, and the dark and light shaded areas are the corresponding 1-2 σ levels, respectively. As a reference, we show in red the power-law model in all panels. All these models include different sources of scatter. For instance, the scatter in Simba-random and Simba-brightest models is based on the variability of the photon escape fraction, and consider an extreme case where some halos (whether randomly or the brightest) have f esc = 1 and the rest are completely obscured (f esc = 0). The Simba-brightest is motivated by the fact that galaxies with high SFR would have stronger feedback (more gas ejection) and hence higher f esc . The other models, Simba or power-law with scatter, include different forms of scatter either from Simba's detailed subgrid physics or by assuming a simple lognormal distribution, respectively, while assuming the same f esc . It is worth noting that the Simba and power-law models are found to obtain same reionization history using the same photon escape fraction of 25%. This shows that the power-law model is a good representative of the average R ion in Simba. In fact, the power-law approximately represents the average of all considered models, expect the Simba-brightest where the models differ by a factor of 2-3. This shows that comparisons with these different models would reveal the impact of variability within either the number of sources or f esc . While the scatter in Simba increases the total R ion by ∼ 20% at the massive end, this increase has a negligible effect on the reionization history since the number density of massive sources (> 10 11 M ) is small as compared to the faint ones. In addition, it is also seen in the top panel of Figure 2 that the power-law has a higher R ion slope at low mass end than the average in Simba (blue solid line), which boosts the total R ion in the powerlaw model to a similar level as compared to Simba. By contrast, the power-law with scatter requires a 20 Figure 3. Comparison between ionization maps for all models, at different redshifts and stages during reionization. All models produce similar maps, indicating that the scatter, either in the ionizing emissivity (Simba vs power-law models) or in the photon escape fractions (comparing all Simba variants), has a minimal impact on the large scale morphology. % lower f esc as compared to the power-law. This is a consequence of assuming that the R ion can be modelled using a log-normal distribution, where the median is somewhat less than the mean. This shows that f esc is highly sensitive to the assumptions used to build the source model. We compare all these models in terms of their large scale ionization morphology in §3.

Radiative Transfer
All reionization realizations from different source models are obtained in post-processing using the radiative transfer module ATON (Aubert & Teyssier 2008, 2010, which is a moment based method with M1 closure relation. We run reionization from z ∼ 20 down to z ∼ 5, using all snapshots available in this redshift range, which are about 35 snapshots in the time interval of about 20 Myr. Single frequency has been used, and temperature has been evolved. We assume a black body spectrum with temperature T = 5 × 10 4 K (E=20.27 eV), and full speed of light (c=1). A similar version of ATON has been used in Keating et al. (e.g. 2020) and Kulkarni et al. (2019), for instance, to study the large scale opacity fluctuations and Lyman-α forest by the end of reionization. The main difference in our analysis lies in considering different source models as explained in the previous section. We first smooth densities of gas particles on a grid with number of cells N=256 3 , resulting in a resolution of 390.625 ckpc/h. While this is somewhat coarse density grid, in Hassan et al. (2016), we have performed a convergence test using the same source model R ion for different resolutions and configurations, and found that the change in the power spectrum remains minimal. Hence, we do not expect the number of cells would alter the conclusion. Using the gas density grids with the ionizing emissivities as inputs, we use ATON to self-consistently compute the ionization field, gas temperature and the photoionization rate. We run ATON with all source models, and compare the results in the following section.

RESULTS: REIONIZATION MORPHOLOGY
We here compare the effect of using the different source models described in the previous sections on reionization morphology. The main difference between these models lies in the scatter in the R ion -M h relation, and hence this study shows the impact of including different forms of scatter on reionization.
In Figure 3 and 4, we show a comparison between all models in terms of maps of the ionization field and their corresponding power spectra at different redshifts, respectively. We can see that the morphology remains mostly unaffected as seen in the maps as well as the power spectra. This is also seen in the ratio between all models to the power-law model (bottom panels in Figure 4), which is close to unity in all cases. However, it is expected for the Simba-random and Simbabrightest models to have a lower small scale power, due to the smaller number of sources as compared to other models. This effect is clearly seen at the early stages of reionization where bubbles are small. As reionization proceeds and bubbles grow, all models seem to have the same power ratio, except for the Simba-brightest model, due to the difference in the global neutral fraction. In addition, the source clustering and bias only change in the Simba-brightest model, which is the reason for its higher large scale power as compared to other models, particularly at the late stages of reionization. All these effects are too small to cause noticeable differences in the maps. It is worth mentioning that we have also checked the opposite scenario to the Simbabrightest model, where only the faintest halos are used with unity f esc . In that case (i.e. Simba-faintest), nearly 50% of halos are required to yield similar emissivity and reionization history to other models. All results remain unchanged whether random, brightest or faintest halos are used. This minimal impact might be due to the fact that hundreds if not thousands of galaxies exist within large ionized bubbles, and hence the scatter would average out. In this case, the scatter impact would be reduced by a factor of ∼ 1/ √ N , where N is the number of sources. Since these comparisons are performed at the same redshift and the 21cm brightness temperature field traces the ionization topology, the same conclusion is valid for the 21cm fluctuations on large scales. This shows that including different forms of scatter in the R ion -M h relation has a minimal impact on reionization morphology. However, we have checked the bi-spectra between these models and found bigger differences (scatter induces a more negative bi-spectrum, hence larger skewness). We leave investigating the impact of scatter on higher order summary statistics to future works.

CONCLUSIONS
We have studied the impact of scatter on reionization morphology by considering a broad range of models with different variations in the photon escape fractions and emissivities, using sources identified within the state-of-the-art galaxy formation simulation Simba . The main improvement over previous works is the inclusion of a realistic scatter in the R ion -M h relation. The radiative transfer is performed using the ATON routine, which implements a moment based method with M1 closure relation. We post-process the radiative trans- Figure 4. Comparison between ionization power spectra for all models, at different redshifts and stages during reionization (top panels). The bottom panels show the power ratio of all models to the power-law model. As seen in Figure 3, all models produce similar power spectra on small and large scales, irrespective of the variability in the emissivity or in the photon escape fraction.
fer on Simba's largest available simulation run, which has a comoving volume of 100 Mpc/h, sufficient to yield convergent reionization histories and capture the bubble overlap epoch by the end of reionization (Iliev et al. 2014).
Besides models with scatter, we constructed a simple source model (power-law) using a fitting function (Equation 3) that behaves as a power-law for halos larger than Simba's resolution limit (> 10 9 M ). This model adopts a one-to-one R ion -M h relation. We show in Appendix A that those missing halos below 10 9 M do not change the reionization morphology but are important to accurately estimate the ionization history and astrophysical parameters such as the photon escape fraction.
To investigate the impact of the scatter, we tune the number of sources and photon escape fractions in all models to produce similar evolution in the comoving ionizing emissivity density and reionization history. Our main key finding can be summarized as follows: The scatter in the R ion -M h relation does not significantly alter the reionization morphology as seen in the ionization maps (see Figure 3) or the power spectra (see Figure 4) if computed at the same global ionization fraction.
It is worthwhile to mention several limitations associated with this work. The radiative transfer is modelled using a moment based method with M1 closure that is known to over-ionize Lyman Limit Systems in post-reionization and does not accurately capture shadowing (e.g. Wu et al. 2021). The radiative transfer also has been run in post-processing and hence these results do not include the coupling between the hydrodynamics and radiative transfer. We also do not account for contribution from binary stars that would induce earlier reionization (e.g. Rosdahl et al. 2018;Doughty & Finlator 2021). Additionally, due to the low resolution in the current Simba run, we have relaxed the requirement to identify dark matter halos, assuming only 16 particles as a minimum to capture smaller halos. While we have checked explicitly that the results remain the same when we require a larger number of dark matter particles (e.g. 64) to identify halos, the comparison between Simbarandom and Simba-brightest with Simba, where a smaller number of sources were used (e.g. 25% and 7%), provides an evidence that, irrespective of the number of sources, models (with/without scatter) at approximately the same neutral fraction, would produce similar reionization morphology on cosmological scales. Recent works (e.g. Ma et al. 2020;Ocvirk et al. 2021) have shown that the photon escape fraction depends on halo mass, and varies by ∼ 2 orders of magnitude between low and high mass halos. Although most of our models adopt a constant f esc , Simba-brightest adopts an extreme case in f esc variability, where 7% of all halos (the brightest) are given unity f esc , and the rest (93%) have f esc = 0. This implies that, regardless of f esc variability with halo mass, the reionization morphology remains relatively similar (comparing Simba-brightest to other models). In Appendix B, we attempt to gain deeper insights on the minimal impact of scatter on the large scale morphology by comparing the ionizing emissivity history for randomly selected halos in terms of the dark matter halo evolution versus the star formation history. Iyer et al. (2020) have performed a detailed comparison between different galaxy formation models (Illus-tris, IllustrisTNG, Mufasa, Simba, EAGLE), and found that the variability of star formation histories within these simulations broadly show a similar power spectral density that is well described by a broken power-law. This indicates that similar results would be obtained using any of these galaxy formation models, since the R ion history follows closely the SF history. In addition, most of radiative transfer hydrodynamic simulations use either the OTVET or M1 approximation to solve the moment-based radiative transfrer eqaution, which have been shown to perform similar during the early and intermediate stages of reionization, but over-ionize selfshielding absorbers by the same amount by end of reionization (see e.g. Wu et al. 2021). This suggests that regardless of the radiative transfer method, similar conclusion might be drawn. As seen in Figure 1, the scatter decreases towards higher masses, which can be resolved within larger box sizes. Since the massive halos are very rare, we do not expect their presence to have a significant impact on the results but a slight modification to f esc might be required to obtain the same reionization history. It is worth mentioning that the scatter would increase the number of sources above a specific luminosity threshold, and hence future surveys with JWST surveys will be able to place tighter constrains on the scatter. We leave performing detailed JWST forecasting using Simba to future works.
Our results suggest that simplified models of ionizing sources, with some calibration of the parameters such as the photon escape fraction, can be used as a viable computationally efficient tool to model reionization on cosmological scales. for guidance on how to use it. We thank the anonymous referee for their constructive comments which have improved the manuscript significantly. SH and RSS acknowledge support from the Simons Foundation. RD  Since our results are based on the Simba simulation, which does not resolve halos below < 10 9 M , we here aim to quantify whether these halos are expected to change the reionization signal on cosmological scales. For this test, we make use of a semi-numerical approach, SimFast21 (Santos et al. 2010). We use the same cosmology, volume and number of cells, to generate the density field grids and halos using SimFast21. We set the minimum halo mass in the excursion set formalism to M h = 10 8 M . We then consider the power-law model with/without scatter to run ATON with the whole population of halos above > 10 8 M , and compare results in Figure 5. We compare these two runs in terms of the reionization topology at relatively the same ionized fraction and redshift. We find that the topology remains unaffected as seen in the ionization maps and in the corresponding ionization power. This indicates that the results presented in the main text remain valid, irrespective of these missing halos (< 10 9 M ).

B. VARIABILITY IN R ion HISTORY
To gain deeper insights on the scatter impact due to the evolution in dark matter or SF, we here focus on the variability in R ion history for individual halos as shown in Figure 6. We first compute the R ion history (solid lines) from Simba using the star formation history, which is a measure of the total star formation occurring over a time interval (∆t), and the metallicity history, which quantifies the mass-weighted metallicity over ∆t. We use the ages, masses and metallicities of all stellar particles associated with halos to compute these histories with ∆t = 10 Myr, and the R ion history is then computed using Equation 2. We next compute the R ion history differently in terms of the dark matter evolution over time (dashed lines), in which we trace halos back in time to identify their progenitors (a) Ionizing emissivity density,Ṅ ion , evolution as a function of redshift.
(b) Reionization history as compared with various measurements.
(c) Ionization maps and power spectra at different redshifts. Figure 5. Test with SimFast21 to examine the contribution of halos with M h < 10 9 M that are missing in Simba due to resolution. This figure indicates that those halos have no impact on the large scale morphology. Figure 6. Examples of the log 10 Rion history for randomly selected halos as a function of lookback time using the star formation history (Simba, solid) and the Rion fitting formula (power-law, dashed) in 10 Myr intervals. Halo masses and redshifts are quoted in the figure legends above the individual panels. The redshift quoted above represents the zero point in the x-axis. The value of Rion in Simba varies by 1 order of magnitude over time. The smooth evolution in the power-law model, which traces the evolution in dark matter halo mass, is a good approximation to the average Rion evolution in Simba. In all cases, the difference in instantaneous Rion between these two models is minimal, and hence similar amounts of photons are emitted. This figure provides insights on why the impact of star formation variability on the global reionization signal is found to be minimal. and use their masses and redshifts in the power-law formula (Equation 3) to compute the R ion . We compare these different forms of R ion histories in Figure 6 for randomly selected halos, and quote their masses and redshift in the figure legend. The R ion from Simba (solid), which traces the star formation history, fluctuates by 1 order of magnitude over intervals of 10 Myr. On the other hand, the smooth evolution in the power-law model, which traces the evolution in dark matter (dashed lines), is a good approximation of the average R ion evolution in Simba. In all cases, the difference in the instantaneous R ion between these different R ion histories is small, and hence similar amount of photons are emitted in each R ion model. If we instead use a larger bin size in time (e.g. ∆t = 30 or 40 Myr), then the star formation history becomes more smooth and similar to the dark matter evolution as seen in the power-law model. This indicates that the scatter might have a minimal impact on the large scale reionization signal, since the average of R ion over time is approximately similar to the smooth evolution of R ion as computed from the evolution in dark matter. This exercise provides insights on why the impact of star formation variability on the global reionization signal is found to be minimal.