The impact of inhomogeneous subgrid clumping on cosmic reionization II: Modelling stochasticity

: Small-scale density fluctuations can significantly affect reionization, but are typically modelled quite crudely. Unresolved fluctuations in numerical simulations and analytical calculations are included using a gas clumping factor, typically assumed to be independent of the local environment. In Paper I we presented an improved, local density-dependent model for the sub-grid gas clumping. Here we extend this using an empirical stochastic model based on the results from high-resolution numerical simulations which fully resolve all relevant fluctuations. Our model reproduces well both the mean density-clumping relation and its scatter. We applied our stochastic model, along with the mean clumping one and the Paper I deterministic model, to create a large-volume realisations of the clumping field, and used these in radiative transfer simulations of cosmic reionization. Our results show that the simplistic mean clumping model delays reionization compared to local density-dependent models, despite producing fewer recombinations overall. This is due to the very different spatial distribution of clumping, resulting in much higher photoionization rates in the latter cases. The mean clumping model produces smaller H II regions throughout most of reionization, but those percolate faster at late times. It also causes significant delay in the 21-cm fluctuations peak and yields lower non-Gaussianity and many fewer bright pixels in the PDF distribution. The stochastic density-dependent model shows relatively minor differences from the deterministic one, mostly concentrated around overlap, where it significantly suppresses the 21-cm fluctuations, and at the bright tail of the 21-cm PDFs, where it produces noticeably more bright pixels. ABSTRACT Small-scale density ﬂuctuations can signiﬁcantly aﬀect reionization, but are typically modelled quite crudely. Unresolved ﬂuctuations in numerical simulations and analytical calculations are included using a gas clumping factor , typically assumed to be independent of the local environment. In Paper I we presented an improved, local density-dependent model for the sub-grid gas clumping. Here we extend this using an empirical stochastic model based on the results from high-resolution numerical simulations which fully resolve all relevant ﬂuctuations. Our model reproduces well both the mean density-clumping relation and its scatter. We applied our stochastic model, along with the mean clumping one and the Paper I deterministic model, to create a large-volume realisations of the clumping ﬁeld, and used these in radiative transfer simulations of cosmic reionization. Our results show that the simplistic mean clumping model delays reionization compared to local density-dependent models, despite producing fewer recombinations overall. This is due to the very diﬀerent spatial distribution of clumping, resulting in much higher photoionization rates in the latter cases. The mean clumping model produces smaller H II regions throughout most of reionization, but those percolate faster at late times. It also causes signiﬁcant delay in the 21-cm ﬂuctuations peak and yields lower non-Gaussianity and many fewer bright pixels in the PDF distribution. The stochastic density-dependent model shows relatively minor diﬀerences from the deterministic one, mostly concentrated around overlap, where it signiﬁcantly suppresses the 21-cm ﬂuctuations, and at the bright tail of the 21-cm PDFs, where it produces noticeably more bright pixels.

The impact of inhomogeneous subgrid clumping on cosmic reionization II: modelling stochasticity

INTRODUCTION
The Epoch of Reionization (EoR) is an important period in the history of the Universe, which encompasses the creation of the first stars and galaxies that subsequently influenced the formation and evolution of latter-day structures. These luminous objects have produced enough UV-radiation to both alter their host galaxy composition and to propagate into the intergalactic medium (IGM), ultimately ionizing it for a second time (Furlanetto et al. 2006;Zaroubi 2012;Ferrara & Pandolfi 2014).
⋆ Contact e-mail: M.Bianco@sussex.ac.uk The key goal of reionization simulations is to provide numerical framework for constraining EoR observables, for example the detection of the 21cm hyperfine transition of neutral hydrogen fluctuations (Bowman & Rogers 2010;Paciga et al. 2013;Yatawatta et al. 2013;Parsons et al. 2014;Jelić et al. 2014;Jacobs et al. 2015;Dillon et al. 2015;Robertson et al. 2015;Ali et al. 2015;Pober et al. 2015;Patil et al. 2017;Mertens et al. 2020;Ghara et al. 2020) and Lyman-α damping wings (Davies et al. 2018;Greig et al. 2019). Such simulations require large volumes, of several hundreds cMpc size in order to correctly derive the cosmic reionization history, to account for abundance and clustering of expected sources and to sample vast regions of the universe for detection of the redshifted 21cm hyperfine transition of neutral hydrogen (Mellema et al. 2006b;Iliev et al. 2014), as relevant for current and upcoming experiments (e.g. LOFAR 1 , SKA 2 ). At the same time, EoR simulations need to include fluctuations in the density distribution down to the Jeans mass of the cold gas, which is in sub-kpc scale, so as to correctly model recombination effects and thus properly track the expansion of ionizing fronts throughout reionization (Park et al. 2016). Unfortunately, because of limited dynamic range, satisfying both of these requirements in a single fully-numerical simulation is currently unachievable, and will remain challenging in the future. Hence, in largescale simulations, the sources and sinks of ionizing radiation often act on scales much smaller than the resolution level and need to be treated using sub-grid prescriptions. Consequently, simulations may adopt incorrect values for various relevant quantities (e.g. density, temperature, gas pressure, etc.) smoothed on the (relatively coarse) grid scales and this could influence the predicted observational signatures. In this work, our focus is on how sub-grid density inhomogeneities are considered within the volume elements of large-scale simulations. Depending on how gas density fluctions vary in space and over time (local degree of "clumpiness"), the recombinations in the IGM can significantly affect the progress and nature of the reionization process. For every ionized atom that recombines with a free electron, an additional ionizing photon should be produced in order to ionize it again and keep the IGM highly ionized. In this way potentially a substantial portion of the sources photon budget could be depleted. In simulations, the recombination rate R is discrete, averaged on a mesh giving R = α B (T)x 2 i n 2 , where α B (T) is the (temperature-dependent) Case B recombination coefficient, x i is the ionized fraction, n is the number density and for simplicity we assumed pure hydrogen gas. This indicates the number of electron-proton recombination per second in a volume, for a given gas chemistry, within each grid cell. Early semi-analytical models have adopted a common methodology named Clumping Factor Approach, that defines the averaged recombination rate in terms of a clumping factor C = n 2 / n 2 , which corrects for the difference between the cell-averaged n 2 and the actual value, thereby accounting for unresolved small-scale (sub-grid) structure in simulations (Gnedin & Ostriker 1997;Tegmark et al. 1996;Ciardi & Ferrara 1997;Madau et al. 1999;Valageas & Silk 2004). If not correctly treated, this approach can underestimate the impact of sub-grid inhomogeneities on absorption of radiation. In some cases this term is just completely ignored, i.e.: C = 1 (Onken & Miralda-Escudé 2004;Kohler et al. 2007), but the more common and simplistic approaches consist in either a constant global term (Cen 2002;Zhang et al. 2007) or a time evolving global term (Iliev et al. 2005;Mellema et al. 2006b;Iliev et al. 2007;Pawlik et al. 2009), averaged on the entire box volume, also referred as the biased homogeneous or globally averaged clumping model. Recently we presented our first work (see Mao et al. 2019, for reference), hereafter Paper I, where we investigated the impact of a spatially varying, local density dependent subgrid clumping factor on reionization observables. In the 1 http://lofar.org 2 https://skatelescope.org present paper we extend the discussion and propose a more realistic and accurate treatment of the Clumping Factor Approach, that takes into account also the scatter around the mean clumping-density relation observed in high-resolution simulations.
We use a high-resolution N-body simulation of a small volume of side length 9 cMpc, with spatial and mass resolution of approximately 200 pc and 5000 M ⊙ , to statistically describe IGM density fluctuations down to the Jeans mass in the cold, pre-reonization gas and then to implement these sub-grid density fluctuations into two large volume (714 cMpc and 349 cMpc os side length) reionization simulations. By adapting the small-scale sub-grid to the resolution of larger boxes we then model the correlation between density and clumping factor, comparing three different models (details in §2.3), in order to infer the clumping factor from the coarse density grid of the large volume, see §2.4. Finally we perform a radiative transfer simulation to study the effect of this sub-grid inhomogeneity approach on observables of reionization.
This paper is organized as follows. In § 2 we present the N-body and radiative transfer (RT) simulation used, the numerical methods, §2.2 and our models in §2.3. In §2.4 we discuss the realisation of the clumping factor for large volumes from sub-grid inhomogeneity correlation. In §3 we analyse our RT simulation results and look into how our models influence the basic features of EoR: the reionization history in §3.1, the volume-averaged ionization fraction evolution, the integrated Thompson optical depth and then the Bubble size distribution in §3.3. To better understand the change in ionization morphology we describe a side-by-side comparison of box slice shot with zoom §3.2. In §3.4 we analyse the 21cm signal power spectra and the brightness temperature distribution. Our conclusions are summarized in §4.

Numerical Simulations
We use N-body simulations to follow the evolution of cosmic structures, performed with the CUBEP 3 M code (Harnois-Déraps et al. 2013). The code uses particle-particle on short-range and particle-mesh on long-range to calculate gravitational forces. We use set of three N-body simulations, whose parameters are summarized in Table 1. Our clumping factor modeling is based on small, high resolution volume box (6.3 h −1 Mpc= 9 Mpc, 1728 3 particles, labelled SB in Table 1). This has sufficient spatial and mass resolution to resolve the smallest halos that can hold cold, neutral gas. Our main larger-volume N-body simulation is referred to as LB-1 (500 h −1 Mpc= 714 Mpc, 6912 3 ≈ 330 billion particles). A smaller simulation, LB-2, (244 h −1 Mpc= 349 Mpc, 4000 3 = 64 billion particles) will be used as comparison to analyse possible influence of box size and resolution in the realisation of sub-grid clumping factor and prove the stability of our method. For both of the large-volume simulations the minimum halos mass resolved is 10 9 M ⊙ , while halos with 10 8 M ⊙ < M halo < 10 9 M ⊙ are implemented using a sub-grid model , thereby all atomicallycooling halos (ACHs) with minimum mass M halo 5×10 8 M ⊙ are included. We are using updated N-body simulations compared to Paper I, we illustrate this further in §A.
Subgrid clumping II 3 Table 1. N-body simulation parameters. Minimum halo mass is 10 5 M ⊙ , 10 9 M ⊙ and 10 9 M ⊙ , corresponding to 20, 40 and 25 particles, respectively in SB, LB-1 and LB-2. In all cases the force smoothing length is fixed at 1/20 of the mean inter-particle spacing. An on-the-fly spherical overdensity halo finder (Harnois-Déraps et al. 2013;Watson et al. 2013), with overdensity parameter ∆ = 130, creates an halo catalogues at given redshift, that is later used as inputs for the radiative transfer simulation. The remaining particles are categorized as part of the IGM. In this work we do not include any effects from minihaloes M halo < 10 8 M ⊙ . Even though these sources could have driven ionization in the early phase of EoR, their effect on later stage is expected to be minor because of molecular dissociation by UV background radiation from primordial luminous sources, up to a point that their contribution is negligible compared to heavier ACHs (Ahn et al. 2009). Initial conditions are generated using the Zel'dovich approximation and the power spectrum of the linear fluctuations is given by the CAMB code (Lewis et al. 2000). The SB N-body simulation starts at redshift z = 300, while LB-1 and LB-2 at z = 150, which gives enough time to significantly reduced non-linear decaying modes (Crocce et al. 2006), while at the same time fluctuations are small enough to ensure linearity of density field at the respective resolutions. The cosmological parameter are based on WMAP 5 years data observation and consistent with final Planck results, for a flat, ΛCDM cosmology with the following parameters, Ω Λ = 0.73, Ω m = 0.27, Ω b = 0.044, H 0 = 70 km s −1 cMpc −1 , σ 8 = 0.8, n s = 0.96 and the cosmic helium abundance η He = 0.074 (Komatsu et al. 2010). Our method is general and can be applied in any cosmological background, but the specific fitting parameters we provide are based on the above values.
We simulate the Epoch of Reionization using the C 2 -Ray code (Mellema et al. 2006b), a photon-conserving radiative transfer (RT) code based on short characteristic ray-tracing. The LB-1 and LB-2 N-body simulations provide the IGM density fields and halo catalogues with masses, velocities, position and other variables, for a total of 76 snapshots, equally spaced in time (∆t = 11.54 Myr) in the redshfit interval z ∈ [6; 50]. For computational feasibility, the density grid is coarsened for the radiative transfer simulation to 300 3 (LB-1), and 250 3 (LB-2). The high-resolution N-body simulation (SB) data input is initially interpolated onto a 1200 3 (SB) grid, which can then be coarsened to the required resolution as discussed in the next section. These grids correspond respectively to cell sizes of length 2.381 cMpc, 1.394 cMpc and 7.5 kpc. For brevity we will refer to these grids as the subgrid volumes for SB, and coarse volumes in LB-1 and LB-2, noted . cr s . Just as in Paper I, the interpolation of the particles onto a grid is performed with a Smoothed-Particle-Hydrodynamic-like method (SPH-like), which then yields coarse-grid density, velocity and clumping fields (see sect. 2.2 in Paper I for details).
Ionization sources for the radiative transfer simulations are characterised by the ionizing photon production rate per unit time N γ , given by where m p is the proton mass, M halo is the total halo mass within coarse-grid cell, ∆t s = 11.53 Myr, the lifetime of stars set equal to the time between N-body snapshots. f γ is the efficiency factor, defined as where f ⋆ is the star formation efficiency, f esc is the photons escape fraction and N i is the stars ionizing photon production efficiency per stellar atom, it depends on the initial mass function (IMF) of the stellar population, e.g. for Pop II (Salpeter IMF) N γ ∼ 4000, the value for f ⋆ and f esc are still uncertain, therefore these parameters can be tuned in order to match the observational constrain that we will discuss in §3. Here we adopt the partial suppression model of (Dixon et al. 2015), whereby for LMACHs located in a neutral cell the efficiency factor is set to f γ = 8.2, while in an ionized cell (above 10%) we set f γ = 5 to account for feedback. For HMACHs the efficiency factor has a constant value of f γ = 5, equivalent to e.g. N i = 5000, f ⋆ = 0.05 and f esc = 0.02.

Coarse-Grid Method
Our clumping factor calculations are based on N-body data and neglect any hydrodynamical effects on the clumping factor. Accounting for the gas pressure provides additional smoothing at small scales and therefore our clumping factor values should be considered as upper limits. Moreover, we are interested in the reionization of the IGM and therefore exclude the halos from our calculations. The contribution of recombinations inside haloes is already taken into account in Equation 1 through the photon escape fraction and should not be counted again.
In order to represent the N-body particles into a regular grid, we adopt the SPH-like smoothing technique described in §2.2 of Paper I, we refer the readers to e.g. Shapiro et al. (1996) for more general details on SPH smoothing methods. In LB-1 and LB-2 simulations we use regular meshes directly produced by the SPH code at the required resolution (the specific values used here are listed in Table 1). In the SB simulation we adopt a more flexible approach, whereby we first produce all quantities on a very fine mesh (here 1200 3 ), which is later coarsened as required in order to approximately match the cell sizes used in the LB simulations. A window mesh function smooths the SB mesh-grid on a coarser-grained mesh, with size defined by Equation 3. The method allows the windows function to overlap. The percentage of overlap N % is chosen in order to achieve the required resolution size of the LBs and at the same time obtain a large enough set of coarsened SB data, since Mesh 3 crs−gr gives the total number of data point then interpolated by the clumping models (see Figure 1).
where Res LB is the coarse grained resolution of the large box and BoxSize SB the box size of the small box. We employ the SB cell-wise quantities expressed with equations (4) and (7) to compute the parametrization of the correlation models.
Hereafter, we will refer to them as the sub-coarse-grid or SB data, whereas in the case of LBs we name them RT-mesh grid. In our case we have Mesh crs−gr = 8 with percentage overlap N % = 53% for 714 cMpc (LB-1) and Mesh crs−gr = 13 with N % = 50% for 349 cMpc (LB-2). We define the gas clumping factor based on the cell-wise averaged quantities (e.g. Iliev et al. 2007 The mean cell over-density is defined where n IGM is the global average of the IGM number density over the entire box volume (in this paper, we always refer to quantities in comoving units).

Modeling the Overdensity-Clumping Correlation
In this work we consider several models for the parametrisation of the correlation between the local coarse overdensity 1 + δ cell and the coarse clumping factor C IGM, cell .

I) Biased Homogeneous Subgrid Clumping (BHC)
The simplest approach is to set a constant (redshiftdependent) clumping factor C(z), for the entire simulation volume (e.g. Madau et al. 1999 ). In our case, we evaluate this globally averaged clumping for every SB simulation snapshot at the appropriate coarse resolution and then fit it with an exponential function of the form: where C 0 , c 1 and c 2 are the fitting parameters. We refer to this model as biased homogeneous clumping (Paper I) since that volume-averaged value is then multiplied by the local cell density to obtain the recombination rate, effectively biasing recombinations towards high-density regions.

II) Inhomogeneous Subgrid Clumping (IC) Model
This model, where the local gas clumping is set based on one-to-one, deterministic relation with the cell density, was first presented in Paper I. We include it here for comparison purposes. The relation of the clumping with the overdensity in Equation 4 is fit by a quadratic function: where x = log 10 (1 + δ cell ) and y = log 10 (C IGM, cell ), the cellwise quantity from SB simulation. For each snapshot z i we evaluate the fitting parameters a i , b i and c i using the coarsegrid field we derived in §2.2.

III) Stochastic Subgrid Clumping (SC) Model
This model, first presented here, aims to account for the natural stochasticity of the relation between local clumping and overdensity, as observed in full numerical simulations. This stochasticity is due to various environmental effects beyond the dependence of the clumping on the local density, and results in a significant scatter around the mean relation used in the IC model ( Fig. 1). We model this scatter from the simple one-to-one relation by binning the SB coarse-grained clumping in several (here five) wide bins of overdensity ∆δ j . In each bin we fit the scatter using a log-normal distribution.
where x = C IGM, cell . For each snapshot z i and bin ∆δ j we evaluate and record the parameters µ ij and σ ij . A stochastic process is then applied to generate lognormal random values from two dimensional uniformly distributed variable u 1 , u 2 ∈ [0, 1], by using a modified 3 Box-Müller transformation.
where µ ij and σ ij are the weighed log-normal parameters for LB-1 and LB-2. Finally, we note that the range of overdensities in the SB simulation is inevitably narrower due to the smaller volume compared to our target reionization volumes. For data beyond the SB limits, for high and low Subgrid clumping II 5 Figure 1. Sample correlation between local coarse IGM overdensity and coarse clumping factor at redshift z = 7.305 for LB-1 resolution (1.394 Mpc cells, left panels) and LB-2 resolution (2.391 Mpc, right). Shown are the coarsened SB N-body data at these resolutions (black crosses), the IC model (deterministic) fit (red line) and the globally-averaged clumping factor (horizontal dashed line). The (blue) points with error bars represent the expected value and standard deviation of the log-normal distribution (see text) in each overdensity bin. Vertical lines (solid grey) indicate the bin limits, whose sizes are adjusted so that each bin contains the same number of data points (approximately ∼ 400 (left) and 100 (right)). For each figure the right panel shows the log-normal distribution (solid line) of the clumping within each density bin vs. the actual data (shadow area), where we include in the legend below each panels a short description of the relevant parameters. [t] Figure 2. Comparison of mean clumping values for the three different models, the redshift evolution of the mean clumping factor for the three models, respectively, shows the range of the standard deviation. On the bottom plot we show the relative where the left image is for LB-1 and right for LB-2. On the upper plot we have the dashed black line is BHC, in red IC and in green SC model, the shadow error in percentage of the difference with BHC.
densities, we fix the mean value to the one given by the IC model, while standard deviation is fixed to the one obtained in the closest density bin.
These distributions are then sampled randomly to create realisations of the clumping in large-volume simulations. A similar approach, but in a different context, has been used previously by Tomassetti et al. (2014) and Lupi et al. (2018), motivated by observation of density distribution in giant molecular clouds.
In Figure 1 we show examples of the resulting parametrization obtained from the three models at redshift z = 7.305, applied at the LB-1 and LB-2 RT resolutions. We show the coarse-grained N-body data, along with the BHC and IC models, as well as the mean, E[X] = e µ+ 1 2 σ 2 , and the standard deviation, SD[X] = e µ+ 1 2 σ 2 e σ 2 − 1, of our proposed log-normal distribution of the stochasticity. On the side plot we show the coarse-data distribution (shadow histogram) and the resulting log-normal fit (solid line) with brief description of the density-bin limits and fitting parameters shown in the legend.

Clumping Implementation in Large-Scale Volumes
We used simulations LB-1 and LB-2 as examples of our method for creating large-volume simulation sub-grid clumping realisation. Results are shown in Figure 3. In the figures we show the N-body data upon which the model is based (black crosses), the volume-averaged clumping factor BHC (black horizontal line), the one-to-one quadratic fit IC (red solid line), the expectation value E[X] and the standard deviation SD[X] of the log-normal distribution in each density bin (blue error-bar points) with the relative bins limits also shown (dashed vertical line). Finally, our SC model . We plot the 38% (0.5 σ), 68% (1 σ) and 95% (2 σ) confidence interval to highlight the realization distribution. Cross point are the coarse SB data used to calibrate the model parameters. In the case of z = 7.305 they correspond to the one of Figure 1 (left column).
Subgrid clumping II 7 clumping realisation (green area) based on the density field of LB-2 is shown with contours corresponding to the 95% (outer), 68% (middle) and 38% (inner) confidence interval. Tables with parameters of the three models used in this paper can be found online 4 . The results illustrate the extend to which each subgrid clumping model reproduces the trends in the direct N-body data throughout the evolution. The BHC (mean-clumping) model roughly matches the peak of the contours and its evolution over time. The IC model (quadractic fit) captures well the general trend of the density-clumping relation and tracks well the highest density of data points. Finally, our new SC model realisation fully reproduces the data, including the scatter around the mean relation. The contours trace the majority of the simulation data quite closely, apart from a few outliers. However, a few things should be noted here.
First, as noted in § 2.3, the large volumes generally sample much wider range of environments than smaller ones used to produce the model, thus inevitably the large-volume realisation should extrapolate to over-densities outside the range sampled by the direct N-body data, for both larger and smaller over-densities. Second, again as discussed above, for statistical reasons we fixed the bin sizes so that they contain same number of data points, which inevitably results in quite uneven bin widths. These are very narrow near the peak density of points and are quite wide for extreme values of the over-density. The combination of these factors yields the 'flairing' of the realisaion (green) points at both large and small values of the over-density, and thus possible minor discrepancies with what we find in simulations. However, only small fraction of the points are in these regions, as demonstrated by the density of points, and therefore it is unlikely this will affect the results in any significant way. Obviously, the IC model is potentially affected in a similar way, since the quadratic fit is used beyond the range of the original data points.
As a consistency check, we compare the redshift evolution of the volume average of the clumping realizations based on the SC and IC models vs. the actual global mean C glob based on the simulated data ( Figure 2). The vertical line indicates the redshift at which the SB simulation was stopped, thus data beyond that is extrapolated. The relative errors of the mean values (bottom panels) are in agreement within the 6 − 7% for LB-2 (right) and within 10% for LB-1 (left), throughout the relevant redshift range 6 ≤ z ≤ 30. At the highest redshifts (z > 30) the errors appear larger, however over that redshift range the density fluctuations are small and thus all clumping factors converge to 1 and do not contribute to the recombination rate. Hence, the local density inhomogeneity does not significantly affect the global averages; however, we expect that the local clumping factor plays a greater role in the recombination and ionization at small scales (e.g. on the H II region size distribution, ionized bubble volume evolution, etc.). The proposed models are roughly consistent with results of previous papers Park et al. (2016), Iliev et al. (2007) and once our the RT-simulation are performed we expect to ob-4 table for model parameters: https://github.com/micbia/SubgridClumping tain similar confirmation from the work of Iliev et al. (2012).

CLUMPING MODEL EFFECT ON OBSERVATIONAL SIGNATURE
The sub-grid clumping model employed affects the local IGM recombination rates, which is then reflected in the derived observable signatures of reionization. In order to understand and try to quantify the importance of this choice, we perform three RT simulations where we fix the source production efficiencies of ionizing photons and vary solely the clumping model. At each time step the precomputed Nbody density fields are used to create a realisation of the corresponding gridded clumping factor, as described in §2.3. These clumping grids are then stored and provided as additional inputs to full radiative transfer simulations with the C 2 -RAY code (Mellema et al. 2006a). Specifically, the simulation used for this section is LB-1. The simulation redshifts span the range z = 40 to 6, for a total of 125 snapshots. The corresponding aperture on the sky vary from 3.6 to 4.7 deg per side, and covers the redshifted 21-cm frequency range from 26 to 45 MHz. The resolution evolves from 43.5 ′′ to 56 ′′ in the spatial direction, and from 0.08 to 0.15 MHz in frequency.

Reionization History
Our results on the reionization histories are presented in Fig. 4 and Table 2. Perhaps counter-intuitively, either of the more realistic, density-dependent clumping treatments (SC and IC) yield somewhat faster evolution and an earlier end of reionization compared to the BHC model. The former models diverge from BHC around z 12, and thereafter the mean reionization is accelerated with a maximum difference atx i = 70%, of ∆z ≃ 0.3 at z ≃ 7.5, corresponding to a time difference of approximately 36 Myr. is delayed by ∆z = 0.1, or 17 Myr. Here there is very little difference between the SC and IC models. Compared to the observational constraints, all three models reionize somewhat early, however these constraints are largely upper limits, and with significant uncertainties. Moreover, our main interest is the relative effect of different sub-grid clumping models, rather than a faithful reproduction of the constraints. During reionization, free electrons scatter CMB photons via inverse Compton scattering, suppressing CMB anisotropies on all scales and introducing polarization on large angular sizes. The contribution from free electron can be quantified by the integrated Thomson scattering optical depth along the line of sight, given by where σ T = 6.65 × 10 −25 cm 2 is the Thomson cross section, c the speed of light and n e is the electron density at a given redshift.
In Figure 5 we plot the volume mean of Equation 12, integrated back in redshift. In agreement with the global reionization histories, the inhomogeneity-dependent models are very similar to each other and are slightly optically thicker than the BHC case, due to the more advanced reionization in the latter. Regardless of this small difference, all three cases are in close agreement with the Planck-LFI 2015 results Ade et al. (2016), which found τ e = 0.066 ± 0.016 corresponding to an instantaneous reionization for redshift z r eion = 8.8 +1.7 −1.4 .
The importance of recombinations throughout reionization could be quantified by the (dimensionless) mean rate of recombinations per hydrogen atom per Hubble time: In Figure 6 (top left) we show the evolution of the mean of this quantity over the full simulation volume (solid lines), as well as averaged only over the over-dense (dashed lines) and under-dense (dot-dashed lines) regions. Colours indicate the model used, as per legend. We also show (bottom left) the relative percentage difference compared to the BHC model. As could be expected, the number of recombinations grows strongly over time, starting close to zero, departing from BHC model around z ∼ 12, and then all reaching ∼ 15 at late times, as more and more structures form over time. Although all models end up at similar values by z ∼ 6, the BHC model lags behind throughout the evolution. The IC and SC models yield very similar values at all times. The over-/under-dense volumes yield much higher/lower number of recombinations, respectively, demonstrating the wide variety of outcomes dependent on the local conditions. Interestingly, the over-dense average for the BHC model results in very similar recombinations to the full-volume averages of SC and IC models, showing that at least on average the clumping in these last models behaves the same way as the over-dense regions in BHC. Overall, the SC model shows a few percent higher recombination rate (∼ 1 − 5%) compared to the IC model. This is most likely due to the stochastic nature of the realization process, also related to the broader scatter in Figure 2 (shaded areas).
In Figure 6 (right panels) we compare the (nonequilibrium) photoionization rates Γ i computed during the run by the C 2 -Ray code. Just as above, all mean photoionization rates are essentially the same until z ∼ 12, after which the BHC model one rises more slowly, lagging behind the other two cases by about factor of 2.5 throughout most of the evolution, eventially catching up by z ∼ 6. The average rates in the overdense regions are higher than the mean (reflecting the inside-out nature of reionization) by a similar amount, while the mean photoionization rates in the under-dense regions lag behind by larger factors, up to several hundred, before again rising steeply and catching up with the mean by z ∼ 6. Interestingly, the mean rate in BHC overdense regions is again very close to the whole volume means of IC and SC models. The average values in the under-dense regions remain the same for all models until much later, z ∼ 9.5, indicating that the specific clumping model has little influence before that redshift.
At first glance, it seems somewhat counter-intuitive that reionization proceeds faster in the denity-dependent models IC and SC, despite their notably higher recombination rates. The reason for this is that in the former cases also the suppression of low-mass galaxies (LMACHs) due to radiative feedback is weaker than in the BHC case, as illustrated in Figure 7. In the BHC case essentially all such galaxies are suppressed by z ∼ 8.5, while in the density-dependent models the suppression is slowed down, allowing LMACHs to last longer in high density regions. This is further clarified in Figure 8   tion of ionized fraction of cells at five different reionization stages,x i = 0.1, 0.3, 0.5, 0.7, 0.9, approximately corresponding to redshift between z ≃ 12 − 6 (see Table 2). The vertical line indicates the partial suppression threshold for LMACHs. Early on (x i = 0.1) the gas clumping has yet had very litte effect, due to the still small ionized fraction and the short time available for recombinations, thus all models yield very similar results, with only BHC showing slightly fewer highly ionized cells. As reionization progresses (x i = 0.3), IC and SC models remain very similar, while BHC is gaining more ionized cells, and at the same time it is starting to show a lack of neutral regions. Starting from roughly mid-point of reionization (x i = 0.5), the dearth of neutral cells becomes ever more prominent whereas the peak of highly ionized cells stays roughly similar for all models. A faint difference between SC and IC is visible at late times, where slightly more cells remain neutral in SC.

Reionization morphology
The globally-averaged quantities discussed above (Figure 4, 5 and 6) give an overall idea of the reionization history. Next step is to understand how the sub-grid gas clumping model affects the propagation of radiation and the local features of reionization. In Figure 9 we show box slice of LB-2 and compare simulation snapshots with similar globally averaged ionized fraction and the three gas clumping models. From top to bottom row we havex i = 0.3, 0.5, 0.7, 0.9 (in Table 2 we list the corresponding redshift at which this occurs and its consequent time delay compared to the BHC model) and from left to right column we have the different models BHC, IC and SC. Red/crimson regions indicate highly ionized cells x i > 0.9, in dark blue neutral regions x i < 0.1, and in green/aquamarine the transition phase x i ≈ 0.5. Within each image we embed a zoom-in region, of 85 cMpc per side, to better appreciate the morphological changes of a randomly selected under-dense neutral clump, as ionized fronts expand (bluer blob, right column plots).
Our simulations reproduce the general reionization features found in previous simulations (e.g. Iliev et al. 2014;?). In high density regions LMACH are the first halos to form. In our simulations they make their first appearance at redshift z = 21, and by z ∼ 12 every volume element contains at least one ionizing source. At first, a modest number of isolated sources, highly clustered on small scale but homogeneously distributed on large scale, start to ionize their surrounding gas, forming small regions of a few Mpc size. The presence of sub-grid gas clumping slows down the propagation of the I-fronts and yields somewhat smaller, more fragmented H II regions. Throughout reionization, these HII bubbles grow and eventually overlap, at which point the ionization process accelerates and many of the smaller bubbles percolate to much larger connected volumes.
The side-by-side comparison shows some notable differences between BHC and the two density-dependent models, with the latter starting at a faster pace, with earlier local percolation, then slowing down compared to the former case. Modest differences appear between the three models in terms of large scale morphology, with a higher degree of ioniziation around early sources in the density-dependent models IC and SC (respectively central and right panel). From around the mid-point of reionization (50% ionization by volume, second row of images) we can see neighbouring growing regions connecting to each other and starting to highly ionize the linking filament. At this point, accordingly to Figure 8, all cells in BHC have surpassed the threshold limit x i = 0.1 for the partial suppression of low mass haloes. For IC (middle) and SC (right column) the degree of ionization around sources is visibly more intense compared to BHC, in fact we can distinguish highly ionized cells clustered around the high density peak, whereas under-dense regions are kept fairly neutral. This diversity is due to the higher recombination rate in inhomogeneity dependent model, shown in Figure 6 (left), that effectively reduces the number of photons able to escape the cells of origin and spread into the neighbour grid elements. This is not the case for BHC, to which clumping factor in high density regions is underestimated and ionizing photons are free to percolate and been absorbed elsewhere in the surrounding IGM, therefore interconnecting filament cells between sources clearly appears extended and in a more advanced neutral-ionized transition (blue/aquamarine).
Later on (x i = 0.7, third row of images), ionized regions have grown substantially and become strongly ionized. A first look suggests similar structure patchiness on large scale, although from the zoom-in we can observe that BHC has a wider and smoother transition between the ionized/neutral phases, whereas IC and SC show a narrower front, allowing more cells that host under-density to stay neutral. When the same transition region dwell across the three model, density dependent model show more irregularity with occasionally one or few cells appearing slightly more ionized then their surrounding.
The morphology differences are more evident at late times (x i = 90%, bottom row of images), whereby HII bubbles connect together to form one vast interconnected highly ionized region. At this stage the vast ionized IGM in IC and SC show variations that follow the higher recombinations due to density fluctuations, which is not the case in BHC model and therefore the same regions appear uniformly highly ionized, x ≈ 1. On the other hand there are no striking difference between IC and SC, except for small variations, of a few pixels of size, on the ionized/neutral boundaries. We suspect that this is numeric artefact due to the stochastic nature of SC. We are developing a more complete clumping model, that we will present in future work, to exclude this uncertainty.

Bubble Size Distribution
One of the key characteristics of reionization, which directly affects all observables is the normalized distribution of bubble sizes R dN/dR or volume sizes V 2 dN/dV of ionized regions (Furlanetto et al. 2004). A number of complementary approaches to calculate these distributions have been proposed (e.g. Friedrich et al. 2011;Lin et al. 2016;Giri et al. 2018a).
Here we employ the Mean-Free-Path (MFP) method to calculate R dN/dR, and the Friends-of-Friends (FOF) algorithm (Iliev et al. 2006) to obtain V 2 dN/dV bubble size distributions (BSD). For both methods we employ the TOOLS21CM 5 python package for EoR simulations analysis (Giri et al. 2020). In both cases, we apply a threshold value of x th = 0.9, since we want to highlight differences in distribution of highly ionized regions that develop around sources.
Results are shown in Figure 10 and Figure 11, respectively we see the typical traits of the percolation process, with volume ranges that roughly corresponding to what is expected from large simulated box . We present our results at four different reionization milestones, x i = 0.3, 0.5, 0.7, 0.9, see Table 2 for corresponding redshifts. In the case of MFP-BSD, we calculate the mean bubble size byR = ∫ (R dN/dR) dR, represented by the corresponding vertical lines for each simulation. The sharp cut-off at small scales 2.381 cMpc, for MFP-BSD, and 13.498 cMpc 3 for FOF-BSD correspond to the simulation cells size and volume respectively.
Early-on (x i = 0.3, top left panel, Figure 10), LB-2 hosts small H II bubbles with radius smaller then 10 cMpc. For inhomogeneity-dependent models IC and SC, distributions present many more highly-ionized regions, indication of a faster radiation propagation around sources. All three distributions peak at the size corresponding to one cell. The same Subgrid clumping II 11 Figure 9. Box slice comparison of LB-2 ionization fraction for different clumping models. In red/crimson highly ionized regions x i > 0.9, in green/aquamarine transition x i ≈ 0.5 and in dark blue neutral regions x i < 0.1. The zoom-in covers an area of 85 cMpc per side and each pixel represent a volume element of 2.381 cMpc per side. We compare slices at same global average ionization fraction, from top to bottom row we havex i = 0.3, 0.5, 0.7, 0.9 (see Table 2 for corresponding redshifts). From left to right column respectively we show the models BHC, IC and SC.   trend is confirmed by the topologically-connected FOF volumes (Figure 11), which are however typically larger than MFP, with volumes between 30 − 700 cMpc 3 for BHC and a wider distribution for IC and SC, from one cell up to a few thousand Mpc 3 . Even though the number of bubbles increase as reionization progress, atx i = 0.5 (top right), the MFP-BSD remain similar. However, the FOF-BSD shows a qualitative transition when the small H II regions start to percolate into much larger, connected one. Their sizes vary widely, with a broad flat distribution (plateau) at smaller scales (V < 10 5 − 10 6 cMpc 3 ). However, BHC and IC also show a bifurcated distribution, with a second peak at large scales, at 10 5 cMpc 3 for BHC and 10 6 cMpc 3 for IC, indicating that percolation process has started (Iliev et al. 2008. Compared to BHC, the IC distribution is shifted toward larger sizes, such that the limit for the plateau and the percolation cluster are up to one order of magnitude higher. A narrower separation between these two volume range indicates that the merging of ionized region in BHC has just started Giri et al. 2018a), whereas in the case of IC this process is already ongoing. On the other hand, IC and SC distribution show similarity at small scale but they differ for larger volumes. The former distribution shows a constant and continuous range of scales from large volumes V ∼ 10 6 cMpc 3 down to one cell, sign that ionized regions are in principle less interconnected and therefore the presence of one dominant super cluster has not yet occurred.
During the later stages of the reionization process (x i = 0.7, bottom left) this bifurcation of the FOF-BSD continues and strenghtens, with ever more small patches merging into the large one, while smaller patches become fewer and on average ever smaller. At this stage the three models present similar volume distributions, whereas their MFP-BSD varies. BHC distribution starts to show a clear characteristic size peak. Albeit of similar shape, the BHC size distribution is clearly shifted to smaller scales, with the average bubble size smaller by a few Mpc and the distribution peak at scale about a factor of 2 smaller (8 vs 15 Mpc) Towards late reionization (x i = 0.9, bottom right), the volume limit for isolated regions to grow before merging is further reduced to V ∼ 10 3 cMpc 3 , while the percolation cluster surpass volumes of 10 8 cMpc 3 (i.e. close to the full simulation volume) in all the three cases. In Figure 10, the sizes distribution in the BHC model has surpassed the other two, with average radius of 54.84 cMpc. IC and SC show again similar distribution but with an increasing, although still minor, difference in the mean radius. Volume distribution in Figure 11 present a similar situation, the only difference between IC and SC consists in the value of the volume merging limit, with a difference up to 1 cMpc 3 .

21-cm Signal Statistics and Power Spectra
The hyperfine transition of neutral hydrogen redshifed into meter wavelengths is a key observable of reionization. Its characteristic emission/absorption line has rest-frame wavelength λ 0 = 21.1 cm and corresponding frequency 1.42 GHz. Radio interferometry telescopes measure the intensity of this signal by quantifying the differential brightness temperature δT b ≡ T b − T CMB signal from patches of the sky, given as: here x HI is the fraction of neutral hydrogen and 1 + δ = n N,IGM /n N,IGM is the local IGM overdensity. The differential brightness is characterized by the relation between the CMB temperature T CMB and spin temperature T S (see e.g. Furlanetto et al. 2006 andZaroubi 2012 for extended discussion). Equation 14 saturates when the neutral hydrogen decouples from CMB photons and couples with the IGM gas heated by X-ray sources (e.g. Ross et al. 2019), so that T S ≫ T CMB , which is the approximation we adopt here. This is known as the heating-saturated approximation where the signal is for the majority observable in emission, δT b > 0, true only at low redshift z < 15. Thus in our simulation the approximated differential brightness is dependent on the density distribution of the neutral gas and redshift, such that δT b ∝ √ 1 + z (1 + δ) x HI . From the RT and N-body simulation outputs we calculate the differential brightness coeval cube at each time step. The cube is then smoothed in the angular direction by a Gaussian kernel with a FWHM of λ 0 (1 + z)/B, where B = 2 km corresponds to the maximum baseline of SKA1-Low core. Smoothing along the frequency axis is done by a top-hat kernel with the same width and the above Gaussian kernel. SKA1-Low will not observe the coeval cube. Instead it will observe a lightcone, in which the signal evolves along the line of sight direction. We construct lightcones from our simulation results using the method described in Giri et al. (2018a). This method is also implemented in TOOLS21CM.
In Figure 12 we show the smoothed lightcone for the three different clumping models, BHC, IC and SC, respectively from top to bottom. This type of data maps the 21-cm differential brightness evolution at the observed frequency ν obs = ν 0 /(1 + z), where ν 0 = 1.42 GHz is the rest frame frequency when the signal was emitted at redshift z. We then express the comoving box length in corresponding angular aperture of 4.65 • at z = 6.583.
Early on, the IGM remains mostly neutral, the average signal largely follows Eq. 14 (δT b > 30 mK) and the fluctuations are driven by the density distribution. The gas clumping also remains low and therefore at low frequencies, ν obs > 120 MHz, there is no visible difference between simulations. As radiation escapes the host halos, it starts to form small isolated transparent regions around sources and gradually suppresses the average signal. The H II regions are still small and thus are smoothed over by the observation beam. Figure 12 shows very similar evolution for the three simulations at frequency higher then 130 MHz (z < 10), but with different intensity of signal suppression. For example the appearances of the first transparent regions, due to lack of neutral hydrogen, at ν obs ≃ 147 MHz and angular position 3.2 • and 1.1 • shows that ionization around sources are more consistent for the simulation with inhomogeneity dependent clumping. This is the case even at higher frequency ν obs > 180 M Hz (z < 7), during the final phases of reionization the morphology and size of the percolation cluster  ∼ 3 mK that are extensively linked together. The IC model shows the same morphology but with considerably smaller and more isolated regions of signal. The SC model, in the other hand, shows a conspicuous lack of signal and regions of emission have only of a few Mpc size.
These differences between models are more clearly observed in the statistics of the 21-cm differential brightness temperature fluctuations. are significant variation in the statistics of the differential brightness temperature -rms, PDFs, skewness and power spectra -shown in Figure 13, Figure 14, and Figure 15. The low frequency cut-off is cho-sen for range where differences between models becomes noticeable. The high density peaks get ionized early, and the corresponding H II regions are smaller then the interferometer resolution, thus their effect on rms (Figure 13, top) is to diminish the averaged δT b without increasing fluctuations. At this stage the signal mostly follows the underlying density field, apart from the peaks and there is little difference between the models. The observed frequency of the RMS dip indicates the timing at which HII regions become larger then the interferometry smoothing scale and eventually start to overlap locally. This is the case at fre- Subgrid clumping II 15 Figure 13. Differential brightness statistic quantities derived from the lightcones data smoothed on the core baseline of SKA1-Low (B = 2 km). Plot on top shows the frequency evolution of the signal root mean squared (RMS). Bottom plot shows the skewness and an inset panel show the frequency evolution of the averaged differential brightness in logarithmic scale. Figure 14. Probability distribution functions of the differential brightness temperature at ionized fractions x i = 0.1, 0.3, 0.5, 0.7 and 0.9, for the three clumping models, as labelled.
quency larger then 120 MHz (z ≈ 11). For the IC and SC models, the turnover occurs earlier and with a steeper slope than the BHC model, indication that signal fluctuations increase faster and stronger. Moreover the peak value of the RMS fluctuations varies, in the case of IC and SC models the amplitude is 14% higher, despite having a lower averaged brightness temperature then the homogeneous case, indicating that the signal is sensitiive to a more physical treatment of the clumping factor. This is the consequence of a lower clumping factor values in under-dense regions, consistent with the conclusion in subsection 3.1. The faster propagation of I-fronts, in the vast low density regions, leads to a earlier second peak in the RMS of the two former approaches. In order of appearance at ν obs = 165 MHz (z = 7.56) for SC, 169 MHz (z = 7.34) for IC and slightly later at 176 MHz (z = 7.06) for BHC, respectively when the average neutral fraction isx n = 0.33, 0.28 and 0.25. The subsequent decline is the results of reionization reaching its final stage, with almost complete ionization. The averaged 21-cm fluctuations level at different scales is reflected in the power spectra (Figure 15), where we compare the results for models BHC, IC and SC at epochs at which the mean ionization fractions are x i = 0.1, 0.3, 0.5, 0.7, 0.9, as well as around reionization completion x i = 0.99. At first, the 21-cm signal follows the underlying density distribution of neutral hydrogen and the power spectra are very similar and approximatively a power law in all three cases. The flattening of the power spectra is an indication of the expanding ionized region, shifting the signal toward larger scales while suppressing small structures. Interestingly, this characteristic appears at the same scale regardless of the clumping model but modest difference in amplitude of signal. The BHC model yields systematically lower power at all scales and at all redshifts except close to overlap. The stochastic relation between local overdensity and clumping factor does not have a large effect throughout most of reionization, and is noticeable predominantly at small scales later on. The most significant differences between IC and SC models emerges at the end of reionization (x i = 0.99), where the SC model has less power on all scales, by factor of up to a few. In fact, at that time the SC model has less power than even the BHC, except at the small scales k > 0.3Mpc −1 .
The 21-cm signal fluctuations are strongly non-Gaussian (e.g. Mellema et al. 2006b;Giri et al. 2019) and therefore are not fully described by the power spectra. We therefore also present the 21-cm differential brightness temperature distribution moments of first (PDFs; Figure 14) and second order (skewness; Figure 13, lower panel). For all the models and all times, 21-cm PDFs are bimodal in nature, which is a clear signature of non-Gaussianity (e.g. Ichikawa et al. 2010;Giri et al. 2018b). Even though all the models show non-Gaussinity, there are significant variations between models. The SC and IC models are much more non-Gaussian, with many more pixels at both high low values. Particularly, they show a very strong tail at high values. This is somewhat stronger for the SC model at all redshifts, indicating that the clumping scatter yields more high brightness temperature peaks, by factor of a few. The signal skewness confirms these observations. It is going from negative to positive symmetry at ν obs ≃ 170 MHz, when the volume ionized fraction is close to x i = 0.6 − 0.7 and the RMS fluctuations reach maximum. Differences between models are noticeable  only later, once the simulation overpass the peak in fluctuations, at frequency larger then 180 MHz. At this point the skewness increases exponentially.

CONCLUSION
Studies of the large scale reionization morphology and its imprint on the observable signatures requires large simulated volumes of a several hundred cMpc per side. Due to computational limitations which limit the dynamic range, uniformly high resolution cannot be achieved in such a volume. Therefore no general model of the local recombinations on scale below the resolution of large numerical simulation exists. Typically a constant value of clumping factor is used, but recently we presented a more general model (Paper I), that depends on the local density, and we demonstrated how an over-simplistic treatment of the clumping factor can have a strong effect on the simulated reionization timescale, topology and size distributions of the ionized region.
In the current work we extend and improve this method by including an empirical stochastic subgrid gas clumping (SC) model (see §2.3) based on the results from highresolution N-body simulation, where the full range of relevant fluctuations is fully resolved. Our approach considers a novel parametrization of the correlation between local IGM overdensity and clumping factor, which take into accounts the scatter due to e.g. tidal forces. We employ a high resolution N-body simulation SB, of spatial resolution 260 pc per side, that resolve the Jeans length of the cold IGM and structure evolution on scale much smaller then the resolution of EoR simulations. The density-binned scatter is then modelled with a log-normal distribution. Those distributions are then randomly sampled to create a realization of the scatter. We then apply our method to the density fields of larger volumes LB-1 (714 cMpc per side) and LB-2 (349 Mpc) to infer its sub-grid clumping factor (see §2.4). Subsequently we post-process the large scale N-body snapshot with C 2 Ray radiative transfer cosmic reionization simulation code, in order to present the impact of various modeling approaches for gas clumping on reionization observables (see §3). We then compare our stochastic model SC with the inhomogeneous clumping model, IC, which is a simple deterministic densitydependent fit, and a globally redshift dependent averaged clumping factor BHC, whereby the subgrid gas clumping is independent of the local density.
We find that density-dependent clumping models, IC and SC, exhibit similar behaviour for globally averaged quantities, meanwhile there is a tangible difference when compared to the volume-averaged model BHC. For instance, the reionization history ( Figure 4) is delayed by as much as ∆z ∼ 0.3 at x i = 0.7 (z ∼ 7.5) and the average neutral fraction decrease swiftly for z < 10. The evolution of ionized regions in IC and SC models is a bit faster due to the on average lower gas clumping factor that decreases gas recombination in the under-dense regions. Meanwhile, as structure formation advance, the higher clumping factor C > 20 in high-density regions considerably increase the recombination rate, such that recombination is twice as effective as in the BHC model case for z < 12. We find that the increase of rate in these regions, due to the different density-dependent gas clumping approach, is responsible for the divergence in the simulated observables. Despite the fact that the over-dense medium constitute a minor fraction of the box volume, compared to the vast under-dense IGM, it is responsible for the majority of recombinations. Our model and the IC method Subgrid clumping II 17 behave similarly, with only 5% of relative error to each other. This difference is mainly due to the broad scatter at high density in the clumping-density correlation plot (Figure 3). The clumping factor for IGM in the proximity of sources, is extremely high C ∼ 100 and the introduced stochasticy can extend it to a factor of few hundreds more. Moreover, the simulated electron scattering optical depth is very similar in IC and SC models and the choice of the clumping model has little effect on the feedback of sources. The density-dependence of the subgrid gas clumping accelerates the propagation of ionizing fronts in the low density IGM (Figure 12), By z < 10 (ν obs > 130 MHz) the regions with low 21-cm signal around the sources are more pronounced than in the BHC case. The differences between the new stochastic approach and the IC model are minor, mostly appearing at late times (z < 7, ν obs > 175), where the SC scenario presents considerably less residual neutral gas then the other two models. These last region of neutral gas are mostly in large voids and distant from any ionizing sources, therefore our interpretation is that at lower redshift the empirical stochastic model becomes predominant in under-dense IGM, accelerating the propagation of ionizing radiation in these regions. Meanwhile, at early stages of reionization the gas recombination in high density region drives the reionization process, resulting in reduction of the ionizing photons propagating into the neutral surroundings.
We compared the simulation-derived observables at the same reionization milestones, x i = 0.1, 0.3, 0.5, 0.7, 0.9. Compared to our previous work, the bubble-size distributions (based on both mean free path and FOF methods) do not show large variation, as an indication that the SC model does not increase the recombination rate in a way that significantly alters the morphology and sizes of the ionized regions. The same conclusion can be deduced from statistics of the 21-cm differential brightness temperature. As we demonstrated in Paper I, the density-dependent model increase the amplitude and shift the fluctuations peak position to lower frequency with a difference of approximately 20 MHz compared to BHC model, and just a few MHz of difference when compared to the SC model. Hence, the peak occurs at stage of reionization that differ only of few percentagex n ≈ 0.3 for SC and IC models and 0.25 for BHC.
The PDFs of the redshifted 21-cm distributions show some notable differences between our models. While all distributions are non-Gaussian, the IC and SC yield significantly more non-Gaussianity, with long tail of bright pixels, which is very different from the BHC model. The bright tail is longer for the SC model compared to IC, predicting many more and brighter pixels at all redshifts.
The power spectra of the 21-cm signal (see Figure 15) show that in early phase of reionization, the BHC scenario yields a weaker signal, when compared to density dependent models on all scales. IC and SC differ somewhat at large scale k < 0.1 cMpc − 1 for x)i = 0.3 − 0.5. This largely disappears byx i = 0.7. Towards the final stages of reionization (x i = 99%) results for three models differ. The IC model predicts the highest signal at all scales, higher by a feactor of a few compared to SC. The BHC model signal is intermediate between them for most except the smallest scales.
The results presented here are not intended as a detailed prediction of the reionization observables, but rather a demonstration that an over-simplistic treatment of the clumping factor can have strong effect on the reionization morphology and thus on simulated observables. The widelyused BHC model, overestimates the rate at which the ionized IGM recombines, and therefore have a strong influence on the timescale of reionization, morphology of the ionized region and the intensity of the expected 21-cm signal. We demonstrated that density dependent model takes better account the cumulative effect of the clumping factor on the gas recombination rate. On the other hand, we have also shown that accounting for the scatter around the average, deterministic local density-clumping relation has only modest effects on the reionization morphology and observables, predominantly towards the end of the reionization process. This indicates that the deterministic IC model is usually sufficient except possibly around and after overlap.
The gas clumping factors presented here should be considered as an upper limit to the actual clumping since they are derived based on high-resolution N-body simulations and thus do not capture the photo-ionization feedback that would suppress small-scale density fluctuations. Consequently it overestimates the recombination rate throughout reionization. We leave a more realistic approach, that follows the feedback effects, and the complex physics of the cold gas (T < 10 4 K) in IGM, for future work. [h] Figure A1. Correlation between local coarse IGM over-density and coarse clumping factor at redshift z = 7.305 for the SB simulation. In red, the IC model interpolation ran with the version 1 of the N-body code, with the solid blue line the same quantity but with the updated code. Lower panel, the ratio between the old and new quantity. Downloaded from https://academic.oup.com/mnras/advance-article/doi/10.1093/mnras/stab787/6185052 by University Hospital Zurich user on 29 April 2021