<figure><img src="image_1.png" /></figure>

# SCORCH. III. Analytical Models of Reionization with Varying Clumping Factors

Nianyi Chen¹, Aristide Doussot², Hy Trac¹, and Renyue Cen³

¹McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA; nianyic@andrew.cmu.edu

²Sorbonne Université, Observatoire de Paris, PSL Research University, CNRS, LERMA, F-75014 Paris, France

³Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

Received 2020 May 7; revised 2020 October 27; accepted 2020 November 7; published 2020 December 22

## Abstract

In the Simulations and Constructions of the Reionization of Cosmic Hydrogen (SCORCH) project, we compare analytical models of the hydrogen ionization fraction with radiation-hydrodynamic simulations. We derive analytical models of the mass-weighted hydrogen ionization fraction from the local ionization balance equations as a more accurate alternative to the widely adopted model based on the volume filling factor. In particular, our model has a recombination term quadratic in the ionization fraction, which is consistent with the two-body interaction nature of recombination. Then, we use the radiation-hydrodynamic simulations to study the clumping factors needed to solve the analytical equations and provide accurate fitting functions. We find that the ionized hydrogen clumping factors from our radiative transfer simulations are significantly different than those from other simulations that use a uniform photoionization background. In addition to redshift dependence, we also see the dependence of the ionized hydrogen clumping factor on ionization fraction, and we incorporate this into our fits. We calculate the reionization histories using our analytical models and clumping factors and compare with widely adopted models, and all of our models achieve <7% difference from simulation results, while the other models have >20% deviations. The Thomson optical depths from reionization calculated from our analytical models result in <5% deviation from simulations, while the previous analytical models have a >20% difference and could result in biased conclusions of the intergalactic medium reionization.

*Unified Astronomy Thesaurus concepts: Reionization (1383); Radiative transfer simulations (1967); Large-scale structure of the universe (902)*

## 1. Introduction

The epoch of reionization (EoR) is a period when the first stars, galaxies, and quasars emit UV photons and ionize the neutral hydrogen in the universe. These photons have a large impact on the state and temperature of the baryonic gas through photoionization and photoheating and hence also affect the structure formation in the late universe. Therefore, we can gain information about the first generation of luminous sources and the status of the intergalactic medium (IGM), as well as constrain astrophysics and cosmology, by tracing the detailed history of reionization.

Due to the limited observational data at this relatively high redshift, knowledge about the sources of ionizing photons and the evolution of the IGM is still incomplete, but some progress has already been made. For example, high-redshift galaxy observations show that galaxies most likely provided the bulk of the ionizing photons (e.g., Bouwens et al. 2015a; Finkelstein 2016; Livermore et al. 2017), but quasars could make some contribution toward the end of reionization (e.g., Madau & Haardt 2015). Planck Collaboration et al. (2020) recently inferred a Thomson optical depth $\tau = 0.054 \pm 0.007$ from measurements of the cosmic microwave background temperature and polarization angular power spectra, implying a late reionization midpoint at redshift $z = 7.7 \pm 0.6$ (e.g., Glazer et al. 2018). Becker et al. (2015) found evidence of a dark Ly$\alpha$ trough extending down to $z \sim 5.5$ in the spectrum of a high-redshift quasar, suggesting that reionization could have ended at $z < 6$ (e.g., Keating et al. 2020; Nasir & D’Aloisio 2020), later than previously assumed.

On the theoretical side, there are three main approaches to studying the EoR. The most accurate and expensive are the cosmological simulations combining N-body, hydro, and

radiative transfer (RT) algorithms to solve the coupled evolution of the dark matter, baryons, and radiation (e.g., Trac et al. 2008; Gnedin 2014; Norman et al. 2015; Semelin et al. 2017; Finlator et al. 2018; Doussot et al. 2019). On the next level of accuracy, there are semianalytical/numerical methods providing an approximate and efficient approach to solving both the spatial and temporal evolution of the reionization process. They are especially useful for making mock observations on large scales (e.g., Furlanetto et al. 2004; Zahn et al. 2007; Alvarez et al. 2009; Santos et al. 2010; Mesinger et al. 2011; Battaglia et al. 2013). The most convenient but least accurate are fast analytical calculations and models, which are preferred for exploring the large parameter space in forecasting or inference studies (e.g., Madau et al. 1999; Miralda-Escudé et al. 2000; Barkana & Loeb 2004; Kaurov & Gnedin 2014).

One of the most commonly used analytical models is a differential equation for the time evolution of the volume filling factor of ionized hydrogen (H II) by Madau et al. (1999, hereafter M99). This model has recently been used to constrain the reionization history and infer the properties of the radiation sources, such as the ionizing emissivity and radiation escape fractions (e.g., McQuinn et al. 2011; Haardt & Madau 2012; Kuhlen & Faucher-Giguère 2012; Bouwens et al. 2015b; Robertson et al. 2015; Price et al. 2016; Madau 2017; Ishigaki et al. 2018). Despite its widespread use, there are several concerns about this model. First, it was deriving Strömgren sphere analysis under the simplifying assumption of isolated H II regions; furthermore, it has generally been applied assuming a constant clumping factor in the recombination term, which is not true in reality.

To study and use the analytical models of the reionization fraction, one would often need to make use of the clumping
factors in order to simplify the calculations. Among different clumping factor definitions, the ionized hydrogen clumping factor is of most interest in the literature. The ionized hydrogen clumping factor accounts for the distribution of ionized hydrogen in the IGM. Most numerical and semianalytical simulations of the EoR choose to incorporate constraints on the clumping factors into their simulations (e.g., Bouwens et al. 2015b; Ishigaki et al. 2015; Robertson et al. 2015; Greig & Mesinger 2017). In addition, many recent studies tend to directly consider this quantity as independent of the redshift and roughly constant, which is yet unproven. On the contrary, some studies have discussed its likely variability (Gorce et al. 2018), while others tried to look for systematic errors in the computation or emphasize a possible scale dependency (Kaurov & Gnedin 2015).

The Simulations and Constructions of the Reionization of Cosmic Hydrogen (SCORCH) project is designed to provide radiation–hydrodynamic simulations, theoretical predictions, and mock observations to facilitate more accurate comparisons with current and future observations. In SCORCH I (Trac et al. 2015), we probe the connection between observed high-redshift galaxies and simulated dark matter halos to better understand the primary source of ionizing photons. By abundance matching galaxy UV luminosities to halo mass accretion rates, we construct a fiducial model for the galaxy luminosity functions that can be extrapolated to fainter magnitudes and higher redshifts. Building on this work, Price et al. (2016) used both parametric and nonparametric statistical methods to constrain the radiation escape fraction $f_{esc}(z)$ from high-redshift galaxies using Hubble Space Telescope and Planck observations. Their inferred results favor increasing $f_{esc}$ toward higher redshift in an approximately power-law relation. With a better understanding of the evolving abundance of high-redshift galaxies and the production of ionizing radiation, in SCORCH II (Doussot et al. 2019), we run and analyze three radiation–hydrodynamic simulations with the same fiducial galaxy luminosity functions but different radiation escape fraction models. The simulations are designed to have the same $\tau_e = 0.06$ and similar midpoints of reionization $7.5 < z < 8$ but different ionization histories. Recently, D’Aloisio et al. (2019) also used these simulations to study the heating of the IGM by hydrogen reionization.

In this paper, we derive an analytical model for the global hydrogen ionization fraction from the local ionization balance equation and study different clumping factors related to reionization with the RadHydro simulations in the context of these analytical models. We highlight the dependency of the clumping factors not only on the redshift but also on the ionization fraction. We also compare the current analytical models of reionization and our own models derived from simulation data in terms of the solved reionization history and the resulting Thomson optical depth. This paper is organized as follows. In Section 2, we first summarize the RadHydro simulations used to conduct our analysis. Then we review the widely adopted M99 model for reionization history before proposing and deriving our own analytical model, and we end this section by describing details about our measurements of the clumping factors from the simulation. In Section 3, we first present simulation data and fits for clumping factors for ionized hydrogen, recombination, ionization, and total hydrogen. Then we compare the reionization histories solved from different analytical models and clumping factors. Finally, we show the

comparison between the Thomson optical depths calculated from different analytical models and the simulation results. We adopt the following cosmological parameters in the simulations: $\Omega_m = 0.30$, $\Omega_\Lambda = 0.70$, $\Omega_b = 0.045$, $h = 0.70$, $n_s = 0.96$, and $\sigma_8 = 0.80$.

## 2. Method

### 2.1. Radiation–Hydrodynamic Simulations

In SCORCH II (Doussot et al. 2019), we present three radiation–hydrodynamic simulations with the same cosmic initial conditions and galaxy luminosity functions but different radiation escape fraction models. The simulations are run with the RadHydro code, which combines N-body and hydrodynamic algorithms (Trac & Pen 2004) with an adaptive ray-tracing algorithm (Trac & Cen 2007). It directly and simultaneously solves collisionless dark matter dynamics, collisional gas dynamics, and RT of ionizing photons. The ray-tracing algorithm uses adaptive splitting and merging to improve resolution and scaling. This code was previously used to simulate both hydrogen and helium reionization (e.g., Trac et al. 2008; Battaglia et al. 2013; La Plante & Trac 2016).

The RadHydro simulations have $2048^3$ dark matter particles, $2048^3$ gas cells, and up to 12 billion adaptive rays in a $50 h^{-1}$ Mpc comoving box. We track five frequencies (15.7, 21.0, 29.6, 42.9, and 74.1 eV) above the 13.6 eV hydrogen ionizing energy. We then compute the incident radiation flux and use it in the computation of the photoheating and photoionization rates needed by the nonequilibrium solvers to solve the ionization and energy equations. The same initial conditions, generated at a starting redshift of 300, are used in all of the simulations. All three simulations are run down to redshift $z = 5.5$.

Using an updated subgrid approach to model the radiation sources, we are able to both populate dark matter halos with galaxies by matching the galaxy luminosity functions and accurately compute the spatial distribution of ionizing sources. Following SCORCH I (Trac et al. 2015), the luminosity–accretion rate $L_{UV}(\dot{M})$ relation is inferred from the halo mass accretion rate and the abundance matching performed by equating the number density of galaxies to the number density of halos.

To generate halo and galaxy catalogs, a particle-particle-mesh (P$^3$M) N-body simulation with $3072^3$ dark matter particles is run using a high-resolution version of the same initial conditions as the RadHydro simulations. A hybrid halo finder is run on the fly every 20 million cosmic years to locate dark matter halos and build merger trees. With a particle mass resolution of $3.59 \times 10^5 h^{-1} M_\odot$, we can reliably measure halo quantities such as mass and accretion rate down to the atomic cooling limit.

These simulations are consistent with the latest Planck observations (Planck Collaboration et al. 2016, 2020), as they have been designed to have a fixed Thomson optical depth $\tau_e = 0.06$. They start with the same initial conditions and have the same galaxy populations but use different radiation escape fraction models $f_{esc}$. Following Price et al. (2016), we choose a two-parameter single power law,

$$f_{esc}(z) = f_8 \left( \frac{1+z}{9} \right)^{a_8}, \qquad (1)$$
<table><caption>**Table 1**
RadHydro Simulations</caption><thead><tr><th>Model</th><th>L [h<sup>-1</sup> Mpc]</th><th>N<sub>dm</sub></th><th>N<sub>gas</sub></th><th>N<sub>RT</sub></th><th>f<sub>8</sub></th><th>a<sub>8</sub></th><th>τ</th><th>z<sub>mid</sub></th><th>Δ<sub>z</sub></th></tr></thead><tbody><tr><td>Sim 0</td><td>50</td><td>2048<sup>3</sup></td><td>2048<sup>3</sup></td><td>512<sup>3</sup></td><td>0.15</td><td>0</td><td>0.060</td><td>7.95</td><td>4.68</td></tr><tr><td>Sim 1</td><td></td><td></td><td></td><td></td><td>0.13</td><td>1</td><td>0.060</td><td>7.91</td><td>5.45</td></tr><tr><td>Sim 2</td><td></td><td></td><td></td><td></td><td>0.11</td><td>2</td><td>0.060</td><td>7.83</td><td>6.54</td></tr></tbody></table>

where $f_8$ is the value of the escape fraction at $z = 8$ and $a_8$ is the exponent that changes between our three simulations. With our three runs, we test $a_8 = 0, 1$, and $2$.

Table 1 summarizes the parameters for the three RadHydro simulations. Here $z_{\text{mid}}$ is the redshift at which 50% of the hydrogen is ionized, and $\Delta_z$ is the redshift interval between 5% and 95% ionization. From the table, we see that the main difference in the three simulations is the reionization histories resulting from different treatment of the escape fraction. Sim 0 has a constant $f_{\text{esc}}$, and reionization starts latest but ends earliest out of the three models. Sim 1 has $f_{\text{esc}}(z)$ varying linearly with $1+z$ and is an intermediate model. Sim 2 has $f_{\text{esc}}(z)$ varying quadratically, and reionization starts earliest but ends latest. For more details on how different models affect other aspects of reionization, please see SCORCH I and II (Trac et al. 2015; Doussot et al. 2019).

## 2.2. Volume Filling Factor Model

We begin the discussion of analytical models of reionization history by reviewing the reference model proposed by M99, which has been widely used in analytical calculations (e.g., Bolton & Haehnelt 2007; Kuhlen & Faucher-Giguère 2012; Bouwens et al. 2015b; Price et al. 2016),

$$ \frac{dQ_{\text{HII}}}{dt} = \frac{\langle \dot{n}_{\gamma} \rangle_{\text{V}}}{\langle n_{\text{H}} \rangle_{\text{V}}} - \frac{Q_{\text{HII}}}{\bar{t}_{\text{rec}}}, \quad (2) $$

where $Q_{\text{HII}}$ is the volume filling factor of ionized hydrogen, $\langle n_{\text{H}} \rangle_{\text{V}}$ is the volume-averaged total hydrogen number density, $\langle \dot{n}_{\gamma} \rangle_{\text{V}}$ is the volume-averaged photon production rate, and the effective recombination time is given by

$$ \bar{t}_{\text{rec}} = \frac{1}{(1 + Y/4X)\alpha(T_0)\langle n_{\text{H}} \rangle_{\text{V}} C}, \quad (3) $$

where $X = 0.76$ and $Y = 0.24$ are the hydrogen and helium mass fractions, respectively; $C$ is the clumping factor; $\alpha$ is the recombination coefficient; and $T_0$ is the temperature of the IGM at mean density, fixed to be $T_0 = 2 \times 10^4$ K throughout this paper in order to match the expectations from star-forming galaxy spectra (e.g., Hui & Haiman 2003; Trac et al. 2008).

It is worth noting that Equation (2) is derived from a constant density hypothesis. It thus implies that the mass- and volume-weighted ionization fractions are equal.

## 2.3. Mass-weighted Ionization Fraction Model

As an alternative to the volume filling factor, in this paper, we propose an analytical model for solving the mass-weighted ionization fraction because it comes from the widely used reionization balance equations and is frequently used when calculating observables. The mass-weighted global ionization fraction can be calculated as the ratio between volume-

weighted number densities,

$$ \begin{aligned} \langle x_{\text{HII}} \rangle_{\text{M}} &= \frac{\sum \rho_{\text{H},\text{r}} x_{\text{HII},i}}{\sum \rho_{\text{H},i}} = \frac{\sum \rho_{\text{HII},i}}{\sum \rho_{\text{H},i}} \\ &= \frac{\langle \rho_{\text{HII}} \rangle_{\text{V}}}{\langle \rho_{\text{H}} \rangle_{\text{V}}} = \frac{\langle n_{\text{HII}} \rangle_{\text{V}}}{\langle n_{\text{H}} \rangle_{\text{V}}}, \end{aligned} \quad (4) $$

where the summation is over all cells in the simulation, and we have used the relation $\rho_{\text{H}} = m_{\text{H}}n_{\text{H}}$ in each of the Eulerian cells.

In order to obtain a self-consistent, rigorous derivation of the global reionization equation for the mass-weighted ionization fraction, we start with the local ionization balance equation of hydrogen (e.g., Gnedin & Ostriker 1997),

$$ \frac{dn_{\text{HII}}}{dt} = \Gamma n_{\text{HI}} + \gamma_{\text{coll}} n_{\text{e}} n_{\text{HI}} - \alpha(T) n_{\text{e}} n_{\text{HII}} - 3Hn_{\text{HII}}, \quad (5) $$

where $n_{\text{HII}}$ is the physical number density of ionized hydrogen, $n_{\text{e}}$ is the physical number density of free electrons, $\gamma_{\text{coll}}$ is the collisional ionization rate, $\Gamma$ is the photoionization rate, $\alpha$ is the recombination coefficient, and $H$ is the Hubble parameter. The right-hand side contains four effects that govern the local ionization of hydrogen: the first two terms are the photoionization and collisional ionization, respectively, which increase the ionized hydrogen number density; the third term is the recombination, which decreases the ionized hydrogen number density; and the last term accounts for the decrease in the physical number density due to the universal expansion. Because the collisional ionization is insignificant compared with photoionization in the low-density IGM regions we are interested in, in subsequent derivations, we will drop the second term and focus on photoionization only.

Taking a volume-weighted average on both sides of Equation (5) and dividing both sides by $\langle n_{\text{H}} \rangle_{\text{V}}$, we obtain the global equation for the mass-weighted ionization fraction,

$$ \frac{d \langle x_{\text{HII}} \rangle_{\text{M}}}{dt} = \frac{\langle \Gamma n_{\text{HI}} \rangle_{\text{V}}}{\langle n_{\text{H}} \rangle_{\text{V}}} - \frac{\langle \alpha(T) n_{\text{e}} n_{\text{HII}} \rangle}{\langle n_{\text{H}} \rangle_{\text{V}}}, \quad (6) $$

where the left-hand side follows from Equation (4).

Instead of the volume filling factor $Q$, we focus on the evolution of $\langle x_{\text{HII}} \rangle_{\text{M}}$, as this quantity is more relevant for calculating observables such as the Thomson optical depth. Note that the fourth, universal expansion term on the right-hand side of Equation (5) has been canceled by the same term on the left-hand side that results from the time derivative of the denominator of Equation (4).

One issue with solving for global reionization directly with Equation (6) is that the last term is difficult to compute. Thus, the recombination term is often parameterized by the clumping factor for quick calculations. If we relate the free electron number density $n_{\text{e}}$ to the density of ionized hydrogen $n_{\text{HII}}$,
assuming helium is singly ionized, by

$$n_e = \left(1 + \frac{Y}{4X}\right)n_{\text{HII}}, \qquad (7)$$

we can then rewrite Equation (6) in terms of the ionization and recombination clumping factors as

$$\frac{d\langle x_{\text{HII}} \rangle_{\text{M}}}{dt} = C_{\text{I}}\langle\Gamma\rangle_{\text{V}}(1 - \langle x_{\text{HII}} \rangle_{\text{M}}) - C_{\text{R}}\langle\alpha(T)\rangle_{\text{V}}\langle n_{\text{H}} \rangle_{\text{V}}\left(1 + \frac{Y}{4X}\right)\langle x_{\text{HII}} \rangle_{\text{M}}^2, \quad (8)$$

where $C_{\text{I}}$ is the photoionization clumping factor (e.g., Kohler et al. 2007) that occurs due to spatial fluctuations in the radiation field,

$$C_{\text{I}} = \frac{\langle \Gamma n_{\text{HI}} \rangle_{\text{V}}}{\langle \Gamma \rangle_{\text{V}} \langle n_{\text{HI}} \rangle_{\text{V}}}, \qquad (9)$$

and $C_{\text{R}}$ is the recombination clumping factor,

$$C_{\text{R}} = \frac{\langle \alpha(T) n_{\text{HII}} n_e \rangle_{\text{V}}}{\langle \alpha(T) \rangle_{\text{V}} \langle n_{\text{HII}} \rangle_{\text{V}} \langle n_e \rangle_{\text{V}}}. \qquad (10)$$

Here we use the most general way of defining the clumping factor; namely, we take the spatial temperature variation into account when taking the global average of the recombination rate instead of assuming a constant temperature.

So far, we have defined the clumping factors and shown how the global ionization fraction can be iteratively solved from the differential equation once the clumping factors are known. Before discussing in more detail the measurements of the clumping factors, we would like to point out that the ionization rate $\Gamma$ is often difficult to compute without directly using the simulations. A more convenient quantity to use in the photoionization term is the ionizing photon production rate, $n_e$. This quantity can be computed from galaxy luminosity functions and the escape fraction without running RT simulations. If we use the photon production rate $n_e$ in the place of the photoionization rate $\Gamma$, Equation (8) becomes

$$\frac{d\langle x_{\text{HII}} \rangle_{\text{M}}}{dt} = \frac{\langle n_e \rangle_{\text{V}}}{\langle n_{\text{H}} \rangle_{\text{V}}} - C_{\text{R}} \langle \alpha(T) \rangle_{\text{V}} \langle n_{\text{H}} \rangle_{\text{V}} \left(1 + \frac{Y}{4X}\right) \langle x_{\text{HII}} \rangle_{\text{M}}^2. \quad (11)$$

Notice that only the first term on the right-hand side changes, and there is no need of an ionization clumping factor $C_{\text{I}}$ in this case. In the following sections, we will use both Equations (8) and (11) to solve for the reionization history and will evaluate the performance of both.

## 2.4. Clumping Factors

We have defined the clumping factors and shown their usefulness in solving for reionization history in the previous section. Now we want to provide details about the physical meaning of clumping factors and how we calculate them from our simulations.

In simulation subgrid modeling and global analytical models, clumping factors are used to account for the excess of recombination or photoionization due to fluctuations in gas density and radiation field, respectively, when solving for ionization fractions in the IGM. When calculating clumping

factors from the simulations, the intrahalo gas is often excluded, because the ionization within the halo is already accounted for in the escape fraction. Since we model $f_{\text{esc}}$ explicitly in our RadHydro simulations, we also exclude the halo gas in our clumping factor computation, and we do so by applying empirical upper limits on the gas overdensity. If the overdensity of a region is greater than this limit, the region is considered to be a part of the intrahalo gas and is not taken into account during the computation.

In order to find a robust functional form and make our models broadly applicable, we choose to study clumping factors under three different density upper limits. We use $\Delta = 50, 100$, and $200$ in units of the global mean gas density as our three density cuts, and the resulting clumping factors are named $C_{\Delta<50}$, $C_{\Delta<100}$, and $C_{\Delta<200}$, respectively. The clumping factors are computed while the RadHydro simulations are running instead of in postprocessing. Doing the calculations at many time steps instead of for a few saved snapshots allows us to better study the redshift evolution.

We have not done extensive convergence tests for the clumping factors due to limited computational resources. In addition to the fiducial high-resolution simulations, we have run low- and medium-resolution versions that have four and two times lower spatial resolution, respectively. However, we have not attempted ultrahigh-resolution simulations with smaller volumes, as our current box is already about the minimum size required to accurately capture larger ionized regions. Furthermore, the calculations and models for physical processes such as cooling, star formation, and feedback are quite sensitive to resolution, especially in high-density regions. While we do see smaller variations in the clumping factors between the two highest-resolution runs compared to the two lowest-resolution runs, it is difficult to do a fair comparison for the reasons stated earlier. In the following sections, we provide estimated uncertainties due to resolution effects from the differences between our mid- and high-resolution runs. Readers can treat our clumping factor measurements as lower limits and use the upper limits provided in the figures as guidance for higher-resolution results.

# 3. Results

Having formulated how to use the clumping factors to solve for reionization history and shown how we define and calculate the clumping factors in the IGM from the RadHydro simulations, in this section, we will provide measurements, fitting formulae, and parameters for the clumping factors. We will also show that with our fits and models for the global reionization balance equation, we can recover the evolution of ionized hydrogen accurately in comparison to direct radiation hydrodynamics simulations.

## 3.1. Total Hydrogen Clumping Factor

The total hydrogen clumping factor $C_{\text{H}}$ is calculated from all hydrogen gas, both neutral and ionized, from the simulation. It is defined as

$$C_{\text{H}} = \frac{\langle n_{\text{H}}^2 \rangle_{\text{V}}}{\langle n_{\text{H}} \rangle_{\text{V}}^2}. \qquad (12)$$

The total hydrogen clumping factor is not used in the ionization equations mentioned in the previous section, but we provide our measurements and fits for it here for completeness.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Total hydrogen clumping factor $C_H$ as a function of redshift for the three definitions given by the simulation (dashed) and the fitting (solid). The shaded bands show the uncertainty estimation from the resolution effect. The relative errors are shown in the lower panel. The fitting function is given in Equation (13), and parameters are given in Table 2.</figcaption></figure>

Figure 1 shows the evolution of the total hydrogen clumping factor $C_H$ for our three definitions of $C_H$ given by the simulation. The total hydrogen clumping factor increases with time, in agreement with the commonly acknowledged fact that the collapse gas fraction increases with time during the EoR. Furthermore, when we apply a higher upper limit in density, we are including more gas around the halos in our calculation. As the gas density around the halos grows at lower redshift, we would expect a higher $C_H$ for the high-density-cut models such as $\Delta < 200$, as can be seen in the plot. To show the uncertainties in $C_H$ from numerical resolution, we also plot in Figure 1 the estimated error bands of the clumping factors based on the differences between the mid- and high-resolution simulations. There is an ~20% difference between different resolutions.

As both ionized and neutral hydrogen atoms are considered indifferently, $C_H$ does not depend on the ionization fraction $x_{HII}$, so we consider its evolution only as a function of the redshift. Moreover, because of its independence of $x_{HII}$, the value of the total hydrogen clumping factor is almost the same for Sims 0, 1, and 2. It does, however, depend on the maximum overdensity $\Delta$ where measurements are taken. Hence, we fit three sets of parameters for different maximum overdensities but not different simulations.

For each density cut, we fit the hydrogen clumping factor by a single power law with a running exponent. The fitting formula is

$$C_H(z) = 1 + C_0 \left( \frac{1+z}{9} \right)^{\alpha(z)}, \quad (13)$$

where $C_0$ is a constant controlling the overall amplitude, and $\alpha$ is a linear function in redshift with $\alpha(z) = \alpha_0 + \alpha_1(z - 8)$.

To fit this function to our simulation data, we use the basin-hopping algorithm (Wales & Doye 1997) and minimize the

<table><caption>**Table 2**
$C_H$: Total Hydrogen Clumping Factor</caption><thead><tr><th rowspan="2">Model</th><th colspan="5">C<sub>H</sub></th></tr><tr><th>C<sub>0</sub></th><th>α<sub>0</sub></th><th>α<sub>1</sub></th><th>δ<sub>rms</sub></th><th>δ<sub>max</sub></th></tr></thead><tbody><tr><td>C<sub>H,Δ&lt;50</sub></td><td>1.62</td><td>-2.27</td><td>-0.107</td><td>0.4%</td><td>2.1%</td></tr><tr><td>C<sub>H,Δ&lt;100</sub></td><td>2.02</td><td>-2.89</td><td>-0.079</td><td>0.5%</td><td>2.6%</td></tr><tr><td>C<sub>H,Δ&lt;200</sub></td><td>2.53</td><td>-3.55</td><td>-0.056</td><td>1.1%</td><td>3.9%</td></tr></tbody></table>

<table><caption>**Table 3**
Uncertainty Range Fitting for $C_H$</caption><thead><tr><th>Model</th><th>C<sub>0</sub></th><th>α<sub>0</sub></th><th>α<sub>1</sub></th></tr></thead><tbody><tr><td>C<sub>H,Δ&lt;50</sub></td><td>[1.1,2.1]</td><td>[-2.6,-2.1]</td><td>[-0.064,-0.13]</td></tr><tr><td>C<sub>H,Δ&lt;100</sub></td><td>[1.3,2.8]</td><td>[-3.2,-2.8]</td><td>[-0.017,-0.11]</td></tr><tr><td>C<sub>H,Δ&lt;200</sub></td><td>[1.5,3.6]</td><td>[-3.7,-3.5]</td><td>[0.030,-0.091]</td></tr></tbody></table>

relative error:

$$\epsilon = \sum_{z=z_0}^{z_n} \frac{|C_{\text{sim}} - C_{\text{fit}}|}{C_{\text{sim}}}. \quad (14)$$

Here we choose $z_0 = 14$ and $z_n = 5.5$ to be the redshift range of our fitting for $C_H$. We optimize three sets of parameters for the three density threshold measurements.

The results of our fitted parameters are shown in Table 2. For low redshifts ($z < 10$), the $\alpha_0$ term dominates in the exponent, as $\alpha_1$ in the linear term is relatively small, and high-$\Delta$ models have more negative $\alpha_0$, which leads to a steeper slope that is physically due to the coincidence of the rapid rise of mass fraction in halos that are both sources and sinks of the reionization process at the last stage of the reionization. As we go to higher redshifts, the linear term in the exponent begins to dominate, and the values for three density cuts cross with each other. We also fit the upper and lower limits of the clumping factors estimated from different resolutions and provide fitting parameters in Table 3, in order to show how the range of parameters corresponds to the range in clumping factors. The fitted curves are plotted alongside the simulation data in Figure 1, and their differences are shown in the lower panel.

## 3.2. Ionized Hydrogen Clumping Factor

One key term in the global reionization equation is the clumping factor in the recombination term. In Section 2.3, we have defined the recombination clumping factor $C_R$ with the spatial variation of the recombination rate. In practice, the recombination rate $\alpha$ is sometimes treated as a constant, and its value is fixed at the average temperature of the medium and independent of the recombination case. In this case, the recombination clumping factor is reduced to the ionized hydrogen clumping factor:

$$C_{\text{HII}} = \frac{\langle n_{\text{HII}}^2 \rangle_{\text{V}}}{\langle n_{\text{HII}} \rangle_{\text{V}}^2}. \quad (15)$$

Previous work has also provided fits for the ionized hydrogen clumping factor and used them to solve for the evolution of the ionization fraction. For example, Shull et al. (2012) used hydrodynamic simulations with a uniform ionizing background to fit $C_{\text{Shull}}$, corresponding to our $C_{\text{HII}}$, and their fitting is used in Madau (2017) for computing the volume filling factor $Q$, which we will discuss in more detail in Section 3.5. Kaurov & Gnedin (2015) also made measurements
<table><caption>**Table 4**
C<sub>HII</sub>: Ionized Hydrogen Clumping Factor</caption><thead><tr><th rowspan="2">Model</th><th colspan="6">C<sub>HII</sub></th></tr><tr><th>C<sub>0</sub></th><th>α<sub>0</sub></th><th>α<sub>1</sub></th><th>β<sub>0</sub></th><th>δ<sub>rms</sub></th><th>δ<sub>max</sub></th></tr></thead><tbody><tr><td>C<sub>HII,Δ&lt;50</sub></td><td>2.07</td><td>-1.43</td><td>-0.213</td><td>-1.32</td><td>5.0%</td><td>9.7%</td></tr><tr><td>C<sub>HII,Δ&lt;100</sub></td><td>2.69</td><td>-1.92</td><td>-0.196</td><td>-1.33</td><td>2.2%</td><td>8.3%</td></tr><tr><td>C<sub>HII,Δ&lt;200</sub></td><td>3.48</td><td>-2.50</td><td>-0.185</td><td>-1.34</td><td>3.1%</td><td>21%</td></tr></tbody></table>

of the ionized hydrogen clumping factor using RadHydro simulations and applied a lower limit of ~0.9 on the ionization fraction in their calculations to only account for ionized gas. Another common practice in the literature is to use a constant C = 3 model for the recombination term (e.g., Kuhlen & Faucher-Giguère 2012; Robertson et al. 2015; Bouwens et al. 2015b), mostly based on the two simulation measurements mentioned above. In this section, we show our measurements from the RadHydro simulations as an alternative to the existing models, which use uniform photoionizing backgrounds.

We first measure the ionized hydrogen clumping factor $C_{HII}$ defined in Equation (15). As said in Section 2.4, we make our measurements with three overdensity thresholds $\Delta = (50, 100, 200)$ to exclude the dense regions around the halos. However, unlike many previous works (e.g., Jeeson-Daniel et al. 2014; Kaurov & Gnedin 2015), we do not impose a lower limit on the overdensity, nor do we have a lower limit on the ionization fraction of the cell. This is because we are using the clumping factors for the purpose of global analytical modeling of reionization history, rather than subgrid models within simulations.

Inspired by the fitting function of $C_H$, we use the same functional form, a single power law with a running exponent, for the redshift dependence. While $C_H$ only depends on redshift, we include an additional dependence on $x_{HII}$ for $C_{HII}$,

$$C_{HII}(z, x_{HII}) = \left[ 1 + C_0 \left( \frac{1+z}{9} \right)^{\alpha(z)} \right] \langle x_{HII} \rangle_M^{\beta_0}, \quad (16)$$

where $\alpha(z) = \alpha_0 + \alpha_1(z - 8)$, and $C_0, \alpha_0, \alpha_1$, and $\beta_0$ are parameters that we will find via optimization. Here the ionized hydrogen clumping factor depends on both redshift and mass-

weighted ionization fraction because we want to allow for different reionization modeling when using our fits.

Similar to the fitting method introduced in the previous section, we use the basin-hopping algorithm and minimize the relative error. Here, instead of fitting the data from one simulation, we now optimize against Sims 0, 1, and 2 simultaneously to account for the ionization fraction dependence, and we fit different sets of parameters for the three density cuts. Here we choose $z_0$ in Equation (14) to be the redshift when $\langle x_{HII} \rangle_M = 1\%$ and $z_n$ when $\langle x_{HII} \rangle_M = 99.9\%$, which results in different redshift ranges for the three simulations. The best-fit parameters are listed in Table 4. From the table, we see that the general trend in the three parameters for the redshift dependency part is similar to that in $C_H$. As for the ionization fraction dependency, we see that the models for the three density cuts all have $C_{HII} \propto x_{HII}^{-1.3}$, so the dependence on ionization fraction is not affected much by the density upper limit.

The resulting evolution of $C_{HII}$ matches the simulation results within 20% for any models and definitions and has mean square errors within 5%. With all parameters monotonic in the density cut $\Delta$, one could interpolate these parameters for intermediary definitions of the clumping factor $C_{HII}$ for further simulations.

Figure 2 shows the ionized hydrogen clumping factors from the three simulations and the fitted curves for the three simulations with the $\Delta>100$ density cut. Here $C_{HII}$ starts off at a relatively high value of ~$10^3$ due to the early ionization of the high-density regions around the sources and then decrease rapidly with redshift as larger regions get ionized. Toward the end of reionization, almost all of the neutral hydrogen is ionized, so the value of $C_{HII}$ becomes close to the total hydrogen clumping factor at order unity.

To demonstrate uncertainties in $C_{HII}$, we also plot in Figure 2 the estimated error bands to the clumping factors based on the differences between the mid- and high-resolution simulations. From the plots, we see that there is an ~20% difference between different resolutions. We also fit the upper and lower limits of the clumping factors and provide fitting parameters in Table 5, in order to show how the range of parameters corresponds to the range in clumping factors.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Ionized hydrogen clumping factors from the three RadHydro simulations and the fitted curves for the $\Delta < 100$ density cut. We show the simulation data (dashed) and estimated resolution uncertainties (shaded regions) along with the fitted data from Equation (16) (solid) for Sims 0 (blue), 1 (purple), and 2 (orange). The lower panels show the deviations of the fits from the simulation data.</figcaption></figure>

<figure><img src="image_4.png" /></figure>
<table><caption>**Table 5**
Uncertainty Range Fitting for $C_{\text{HII}}$</caption><thead><tr><th>Model</th><th>C<sub>0</sub></th><th>α<sub>0</sub></th><th>α<sub>1</sub></th><th>β<sub>0</sub></th></tr></thead><tbody><tr><td>C<sub>HII,Δ&lt;50</sub></td><td>[1.7,2.4]</td><td>[-1.7,-1.3]</td><td>[-0.28,-0.17]</td><td>[-1.3,-1.3]</td></tr><tr><td>C<sub>HII,Δ&lt;100</sub></td><td>[2.1,3.2]</td><td>[-2.2,-1.9]</td><td>[-0.27,-0.17]</td><td>[-1.3,-1.4]</td></tr><tr><td>C<sub>HII,Δ&lt;200</sub></td><td>[2.4,4.4]</td><td>[-2.7,-2.5]</td><td>[-0.22,-0.18]</td><td>[-1.3,-1.4]</td></tr></tbody></table>

Figure 3 shows the dependence of $C_{\text{HII}}$ on the density cuts. Similar to $C_{\text{H}}$, lower density cuts (e.g., $\Delta < 50$) exclude more regions around the halo and result in lower clumping, which also leads to a lower total recombination rate. The lower panels show the deviations of the fits from the simulation data, and we can see that with the fitting formula given in Equation (16), we can fit the three simulations within 5% error during most of reionization, although the fitting can get larger than the data values at high redshift due to the limitation of the functional form.

Comparing to previous works on $C_{\text{HII}}$ (e.g., Pawlik et al. 2009; Shull et al. 2012; Jeeson-Daniel et al. 2014; Kaurov & Gnedin 2015), our measured values are higher at $z > 8$, and the reason is twofold. First, our simulations use RT to track the evolution of photons and gases, so the process is patchy throughout reionization. This leads to higher patchiness and a higher $C_{\text{HII}}$ compared with the simulations that turn on a uniform ionizing background around $z = 8$. Second, as mentioned before, we did not apply any lower limit on either $\Delta$ or $x_{\text{HII}}$. This means that we have included all cells outside halos, leading to a larger $C_{\text{HII}}$ at high redshift when the ionization fraction is very inhomogeneous.

### 3.3. Recombination Clumping Factor

As mentioned in Section 2.3, instead of assuming a constant recombination coefficient, it is more accurate to take the spatial variation of $\alpha$ into account when calculating the recombination clumping factor. Hence, we also calculated $C_{\text{R}}$, defined in Equation (10) from our three simulations, and we provide our results and fits here. In our calculation, one slight difference from the definition in Equation (10) is that instead of using $\langle \alpha \rangle$ in the denominator, we choose a fiducial value of $\alpha = \alpha_{\text{B}}(2 \times 10^4 \text{ K})$, such that our renormalized recombination

clumping factor $C_{\text{rec}}$ is defined as

$$C_{\text{rec}} = \frac{\langle \alpha(T)n_e n_{\text{HII}} \rangle}{\alpha_{\text{B}}(2 \times 10^4 \text{ K}) \langle n_e \rangle \langle n_{\text{HII}} \rangle} \quad (17)$$

In practice, this is only a matter of multiplication by a constant, and we find it more convenient to choose a widely used value. Readers can easily rescale $C_{\text{rec}}$ differently if they want.

Figures 4 and 5 show our measurements of $C_{\text{rec}}$ from the RadHydro simulations for three simulations and three density cuts, respectively. Same as $C_{\text{HII}}$, we also plot the uncertainty estimation from resolution effects with the shaded bands. The general trend and magnitude resembles the ionized hydrogen clumping factor. Therefore, we use the same functional form (Equation (16)) as in $C_{\text{HII}}$ to fit the recombination clumping factor,

$$C_{\text{rec}}(z, x_{\text{HII}}) = \left[ 1 + C_0 \left( \frac{1+z}{9} \right)^{\alpha(z)} \right] \langle x_{\text{HII}} \rangle_{\text{M}}^{\beta_0}, \quad (18)$$

where $\alpha(z) = \alpha_0 + \alpha_1(z - 8)$, and $C_0$, $\alpha_0$, $\alpha_1$, and $\beta_0$ are parameters we want to optimize.

Like in $C_{\text{HII}}$, we also fit three sets of parameters for three different density cuts, and each set of parameters is fitted to three simulations simultaneously. The parameters from the fits are presented in Table 6, with uncertainty ranges shown in Table 7. From the tables, we see that the dependence of $x_{\text{HII}}$ on redshift is similar for $C_{\text{HII}}$ and $C_{\text{rec}}$ and almost the same for the three density cuts. The exponent in the redshift dependence has larger linear terms, leading to larger slope at high redshifts.

For comparison, we add the data of $C_{\text{R}}$ from Finlator et al. (2012) to Figure 4. The data were originally in $x_{\text{HII},\text{V}}C_{\text{R}}$ and used $\alpha(10^4 \text{ K})$ for the averaged recombination coefficient (note that their $C_{\text{R}}$ is defined similarly to our $C_{\text{rec}}$, although the normalization is different), but to show a direct comparison with our results, we divide it by $x_{\text{HII},\text{V}}$ and rescale to $\alpha(2 \times 10^4 \text{ K})$. After rescaling, we see from the plot that our data and fits agree with the data from Finlator et al. (2012), although there are some discrepancies due to our different reionization histories and box sizes. Reionization is later in our simulations, and the difference of ionization fractions enters into the difference in the clumping factor, where later reionization can lead to higher clumping at the same redshift. Furthermore, the

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Ionized hydrogen clumping factors with estimated uncertainties from Sim 1 for different density cuts, $\Delta < 50$ (blue), 100 (purple), and 200 (orange). With a higher density cut, we are including more high-density gas around the halos, and the overall clumpiness is higher.</figcaption></figure>
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Measurements of the recombination clumping factors from the three RadHydro simulations as a function of redshift with the $\Delta < 100$ density cut. We show the simulation data (dashed lines) together with the estimated uncertainties (shaded regions) and our fitting results (solid lines) for Sims 0 (blue), 1 (purple), and 2 (orange). We also show the results from Finlator et al. (2012; black dotted) for comparison. The lower panels show the percentage difference between the fitted and simulation results. In most periods throughout reionization, our fits are within 10% of the simulation results.</figcaption></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Effect of different density cuts on the recombination clumping factor within Sim 1. The effect of density cuts on $C_{\text{rec}}$ is not very prominent compared to that on $C_{\text{HII}}$.</figcaption></figure>

<table><caption>**Table 6**
$C_{\text{rec}}$: Normalized Recombination Clumping Factor</caption><thead><tr><th>Model</th><th>C<sub>0</sub></th><th>α<sub>0</sub></th><th>C<sub>rec</sub><br>α<sub>1</sub></th><th>β<sub>0</sub></th><th>δ<sub>rms</sub></th><th>δ<sub>max</sub></th></tr></thead><tbody><tr><td>C<sub>rec</sub>,Δ&lt;50</td><td>2.96</td><td>-1.31</td><td>-0.310</td><td>-1.39</td><td>4.9%</td><td>14%</td></tr><tr><td>C<sub>rec</sub>,Δ&lt;100</td><td>3.38</td><td>-1.50</td><td>-0.341</td><td>-1.41</td><td>2.0%</td><td>8.6%</td></tr><tr><td>C<sub>rec</sub>,Δ&lt;200</td><td>3.57</td><td>-1.58</td><td>-0.354</td><td>-1.42</td><td>3.7%</td><td>19%</td></tr></tbody></table>

<table><caption>**Table 7**
Uncertainty Range Fitting for $C_{\text{rec}}$</caption><thead><tr><th>Model</th><th>C<sub>0</sub></th><th>α<sub>0</sub></th><th>α<sub>1</sub></th><th>β<sub>0</sub></th></tr></thead><tbody><tr><td>C<sub>rec</sub>,Δ&lt;50</td><td>[2.4, 3.5]</td><td>[-1.5, -1.2]</td><td>[-0.31, -0.32]</td><td>[-1.4, -1.4]</td></tr><tr><td>C<sub>rec</sub>,Δ&lt;100</td><td>[2.6, 4.1]</td><td>[-1.7, -1.4]</td><td>[-0.32, -0.36]</td><td>[-1.4, -1.4]</td></tr><tr><td>C<sub>rec</sub>,Δ&lt;200</td><td>[2.7, 4.5]</td><td>[-1.7, -1.6]</td><td>[-0.32, -0.38]</td><td>[-1.4, -1.5]</td></tr></tbody></table>

box size from Finlator et al. (2012) is eight times smaller than ours, which leads to differences in the source distribution and can also result in our higher clumping factors compared with

theirs. Also from Finlator et al. (2012), we can see that if one intends to use our models for a simulation subgrid, it is possible to use $x_{\text{HII},V}C_R$ instead of applying an ionization fraction lower limit at $x_{\text{HII}} = 0.95$, as the evolution of the two values is close to each other.

### 3.4. Photoionization Clumping Factor

To use Equation (8) with the ionization rate $\Gamma$, we also need to compute the ionization clumping factor $C_1$ defined in Equation (9). The measurement of $C_1$ was previously computed by Kohler et al. (2007) using a small ($4 h^{-1}$ Mpc) box. In their work, the focus was toward the very end of reionization due to the interest in $\text{Ly}\alpha$ lines, and the evolution of the ionization clumping factor in the long redshift interval before the end was ignored in the fit. In our measurement, we want to track the evolution of $C_1$ throughout reionization.

In Figure 6, we show the evolution of $C_1$ as a function of redshift and mass-weighted neutral fraction, together with our uncertainty estimations. From the left panel, we can see a turnover at $z = 6-7$, when reionization is mostly complete.
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Evolution of the photoionization clumping factor $C_I$ as a function of redshift (left) and mass-weighted neutral fraction (right) for three RadHydro simulations. Same as in Figure 1, the estimated uncertainties are shown by the shaded regions.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** The M99 model reionization histories in comparison with the simulation reionization histories. From left to right are the reionization histories from Sims 0, 1, and 2. We calculate the volume filling factor $Q$ with the fixed recombination factor $C = 3$ (red dashed), $C_{\text{Shull}}$ from Shull et al. (2012) following Madau (2017, yellow dashed), and our fitted recombination clumping factor $C_{\text{rec}}$ (green dashed). We plot the calculated results of $Q$ against the mass-weighted (blue solid) and volume-weighted (purple solid) reionization fractions from the RadHydro simulations. In the lower panels, we show the fractional difference of the M99 histories from the simulation results.</figcaption></figure>

Before this turnover redshift, the value of $C_I$ was decreasing with time; this is because the ionization front propagates to the larger regions in the IGM where the process of photoionization became less concentrated in space. Then, toward the end of reionization, the photonization rate was once again dominated by the remaining neutral regions in the IGM close to halos where the density is high, so the ionization clumping factor has a steep increase. The fluctuations seen in the plots are due to the episodic star formation in a finite-sized box.

### 3.5. Ionization Fraction

Having fitted the clumping factors needed in solving the differential equation for the evolution of the hydrogen ionization fraction, we now compare our models (Equations (8) and (11)) calculated with the recombination clumping factor to the M99 formalism calculated with constant $C = 3$, the Shull et al. (2012) fit, and our fit for the recombination clumping factor $C_{\text{rec}}$.

First, we calculate the evolution of the volume filling factor $Q$ by solving Equation (2) following Madau (2017), where we assume an average temperature of $T = 2 \times 10^4$ K and use $C = 3$ and $C_{\text{Shull}} = 2.9((1+z)/6)^{-1.1}$ to model the

recombination rate. We choose these two clumping factors because the $C = 3$ model is a widely adopted simplification and $C_{\text{Shull}}$ reproduces the Madau (2017) calculation. For comparison purposes, we also use our fitted clumping factor $C_{\text{HII}}$ together with the M99 model to solve for the volume filling factor $Q$, although we caution the reader that this is not a self-consistent solution. The photoionization term is approximated by $\dot{n}_{\gamma}$ as in M99, where $\dot{n}_{\gamma}$ is directly measured from the simulations (see Trac et al. 2015 for details).

In Figure 7, we show the M99 model reionization histories with $C = 3$, $C_{\text{Shull}}$, and $C_{\text{rec}}$ in comparison with the simulation reionization histories. From left to right are the reionization histories from Sims 0, 1, and 2. We plot them against the mass- and volume-weighted reionization fractions from the RadHydro simulations. In the lower panels, we show the percentage difference of the M99 histories from the simulation results.

To better isolate the effect of the analytical modeling from the use of different clumping factors, let us first look at the volume filling factors solved with our fitted $C_{\text{rec}}$. We can see from the plots that solving the M99 equations with $C_{\text{rec}}$ results in later reionization compared to the simulations, and this is due to the fact that the M99 recombination term, which is linear
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Evolution of $\langle x_{HII} \rangle_M$ from the analytical model proposed in this work, solved with ionizing photon production rate $\dot{n}_\gamma$ (red dashed) and photoionization rate $\Gamma$ (yellow dashed). We plot them against the mass-weighted (blue solid) and volume-weighted (purple solid) reionization fractions from the RadHydro simulations. The resulting ionization fractions follow those from the simulations much closer than the M99 models.</figcaption></figure>

in $Q$, overestimates the recombination rate and thus postpones the overall ionization of hydrogen. Now if we instead use the $C = 3$ and $C_{\text{Shull}}$ clumping factors, as was done in Madau (2017), we will see the opposite effect of an earlier reionization compared to both our simulation results and the $z_{\text{rei}} = 7$ value in Shull et al. (2012). The main reason for this difference is that the recombination clumping factors measured from our RT simulations are much larger than those measured from a uniform ionizing background, especially during the early stage of reionization. This underestimation of the recombination clumping factor overcorrects for the overestimation of recombination rates due to the linear term in $Q$. In either case, the solved $Q$ can deviate from the simulated reionization histories by up to 20%.

To use the volume filling factor model to solve for the reionization history, it is necessary to redefine the clumping factor used in the recombination term. The clumping factor should be $C_{R,Q} = C_R Q$, where $C_R$ follows the definition in Equation (10), in order for the recombination term to match that in Equation (11). Based on the comparison with Finlator et al. (2012) in Section 3.3, we can see that $C_R Q$ is close to a recombination clumping factor with a >95% ionization fraction cut, so one may also use the $C_R(x_{HII} > 95\%)$ calculated within highly ionized regions together with the M99 model.

Next, we use our model in Equations (8) and (11) together with our fits for $C_{\text{rec}}$ in Equation (18) to solve for the mass-weighted ionization fraction. Here we choose the $C_{\text{rec}}$ with $\Delta < 100$ to solve for reionization history because it is consistent with our definition of the escape fractions in the RadHydro simulations. This density cut also falls between the mean overdensity of spherical top-hat halos and the overdensity at the virial radius of an isothermal DM halo during reionization (Pawlik et al. 2009). Using $C_{\text{rec},\Delta<50}$ instead will lead to slightly earlier reionization as the recombination rate becomes lower, and using $C_{\text{rec},\Delta<200}$ will lead to slightly later reionization. Note that $C_{\text{rec},\Delta<50}$ and $C_{\text{rec},\Delta<200}$ are not consistent with our reionization rates in Equations (8) and (11), although in either case, the effect of density cuts on solving the reionization history is small, with a maximum difference of <5%.

Figure 8 shows the evolution of $\langle x_{HII} \rangle_M$ from the analytical model proposed in this work, solved with an ionizing photon production rate $\dot{n}_\gamma$ and a photoionization rate $\Gamma$. The resulting ionization fractions follow those from the simulations much closer than the M99 models. In particular, the ionization fraction solved with the photoionization rate $\Gamma$ shows almost no deviation from the simulation results. The model with $\dot{n}_\gamma$ also does not deviate much from the simulation data, with an at most 5% difference in Sim 0 toward the end of reionization.

In the above solutions, we have used $\dot{n}_\gamma$ and $\Gamma$ directly from our simulations. However, for future applications, people might have different $\dot{n}_\gamma$ and $\Gamma$ from ours. Our analytical model is derived from a first-principle local ionization balance equation, so it has the generality to be applicable to different reionization scenarios. In our three simulations, for example, we also have different ionization rates, and the model works fine in all three cases. One has to be cautious and make sure that the photon production rate is either defined consistently with our density cuts or interpolated from our values shown in the tables. Also, since there is a redshift and ionization fraction dependence in the clumping factors, when using different ionization rates, one has to recalculate the clumping factors at each redshift iteratively when solving for the reionization history. As long as there is consistency in the definition of all variables in the equation, it is fine to use our global analytical equations in other reionization models.

As was mentioned in Section 3.3, our clumping factors are subject to uncertainties. Therefore, we test the sensitivity of our model to different variations in the clumping factors within the uncertainties. We use the upper- and lower-limit fittings of the recombination clumping factors shown in Table 7 to solve for the mass-weighted ionization fraction and find that using the upper limits will result in a maximum of 4% lower ionization fraction, while using the lower limits will result in a maximum of 4% higher ionization fraction when we keep the photon production rate the same. Thus, we conclude that our analytical model is not very sensitive to the variation in clumping factors from the resolution limits.

### 3.6. Thomson Optical Depth

One goal of computing the ionization fraction of hydrogen through the analytical models is to infer the Thomson optical
<table><caption>**Table 8**
Thomson Optical Depth</caption><thead><tr><th rowspan="2">Model</th><th>M99,</th><th>M99 + Shull12</th><th colspan="2">τ<sub>e</sub></th><th rowspan="2">RadHydro</th></tr><tr><th>C = 3</th><th></th><th>Our Model</th><th>Our Model</th></tr><tr><th></th><th></th><th></th><th>w/ Γ</th><th>w/ ˙n<sub>γ</sub></th><th></th></tr></thead><tbody><tr><td>Sim 0</td><td>0.065</td><td>0.067</td><td>0.059</td><td>0.060</td><td>0.060</td></tr><tr><td>Sim 1</td><td>0.066</td><td>0.069</td><td>0.059</td><td>0.060</td><td>0.060</td></tr><tr><td>Sim 2</td><td>0.067</td><td>0.070</td><td>0.061</td><td>0.059</td><td>0.060</td></tr></tbody></table>

depth $\tau_e$, which is an important observable for reionization and parameter for cosmological surveys. It is defined as

$$ \tau_e(z) = \sigma_T \int_0^z dz' \frac{dt}{dz'} \langle n_e(z') \rangle_v, \quad (19) $$

where the volume-averaged free electron number density,

$$ \begin{aligned} \langle n_e \rangle_v &= \langle n_{HII} \rangle_v + \langle n_{HeII} \rangle_v + 2 \langle n_{HeIII} \rangle_v, \\ &= \langle x_{HII} \rangle_M \langle n_H \rangle_v + \langle x_{HeII} \rangle_M \langle n_{He} \rangle_v \\ &\quad + 2 \langle x_{HeIII} \rangle_M \langle n_{He} \rangle_v, \end{aligned} \quad (20) $$

is related to the volume-averaged number densities and mass-weighted ionization fractions for hydrogen and helium. The redshift integration is performed from the present up until the beginning of reionization.

We calculate $\tau_e$ using different models for the ionization fractions and compare how the choice of reionization models affects the value of $\tau_e$. In our calculation, we assume single ionization of helium until redshift $z = 3$, after which helium gets fully ionized. For each reionization history model, we calculate $n_e$ from $\langle x_{HII} \rangle_M$ or $Q$ (note, though, that the correct calculation of $n_e$ should use the mass-weighted ionization fraction) and integrate from the beginning of reionization down to redshift $z = 0$.

Table 8 shows the optical depths from four different models in comparison with the $\tau_e$ values from the RadHydro simulations. Not surprisingly, just like the reionization history, the calculation based on Equation (8) follows the simulation results most closely, but the approximated method with $\dot{n}_\gamma$ in Equation (11) also resembles the simulations in terms of $\tau_e$. With the two models following M99, the values of $\tau_e$ are systematically higher by 5%–15% due to the earlier reionization of the gas, as already shown in Figure 7.

In addition, note that Table 8 shows the total optical depth, which includes both the reionization and post-reionization contributions. Out of the two, the post-reionization term contributes ~0.03 to all of the models equally, so if we only consider the difference between the reionization optical depth when constraining reionization models, the deviation from the simulation results could be as large as 20%–30% with the M99 models.

Given the latest $\tau_e$ measurement from Planck Collaboration et al. (2020) of $0.054 \pm 0.007$, which has a 13% uncertainty, conclusions about the IGM using the optical depth with the M99 model could be systematically biased, while the error in $\tau_e$ calculated with our model is well within the tolerance of current observations.

## 4. Conclusion

In this paper, we address issues in the analytical models proposed in Madau et al. (1999) based on Strömgren sphere

analysis and propose an alternative model to solve for the mass-weighted ionization fraction based on the local ionization balance equations. There are two main differences between our model and the M99 model. First, our derivation leads to a recombination term quadratic in the ionization fraction, while the M99 model uses a linear term. Second, we do not assume constant clumping factors of ionized hydrogen. We also do not use a globally averaged temperature to calculate the recombination coefficient $\alpha(T)$, although this has a minor effect in solving for the reionization history compared with the previous differences.

To study the evolution of clumping factors and include spatial fluctuations in the recombination term, we use the results of the RadHydro simulations from the SCORCH project. The measurements are based on three RadHydro simulations described in detail in Doussot et al. (2019). These simulations assume identical initial conditions and parameters, except for the evolution of the ionizing photon escape fractions. The escape fractions $f_{esc}$ are constant, linear, and quadratic in the three simulations, leading to different reionization histories. To match the idea of temporal and spatial variation of clumping factors and recombination coefficients in the analytical model we propose, those three simulations do not compel the clumping factor to be constant but rather record its free evolution as another output result.

Our first key result is the time evolution of the ionized hydrogen clumping factor $C_{HII}$. We show that $C_{HII}$ depends on both redshift and the mass-weighted ionization fraction. We measured $C_{HII}$ from the three RadHydro simulations using three different definitions, setting the overdensity cut at $\Delta < 50, 100$, and $200$ to exclude the regions within the halos. We clearly see the strong dependence of $C_{HII}$ on redshift, where $C_{HII}$ starts off at a very high value due to the patchiness of reionization and drops down quickly to order unity toward the end of reionization. Based on the data from all three simulations, we provide empirical fitting functions of $C_{HII}(z, \langle x_{HII} \rangle)$ and show that it fits the measured $C_{HII}$ to within 20% throughout reionization.

Next, we take into account the spatial variation of temperature and measure the normalized recombination clumping factors $C_{rec}$ in the same way as $C_{HII}$. We also provide fits for $C_{rec}$, and the fits are within a 20% relative error from the simulation results. With caution, our fit can be interpolated for other empirically set density limits and, to some extent, other reionization histories.

In addition to the clumping factors in the recombination term, we also show the redshift dependence of the ionization clumping factor $C_I$ and total hydrogen clumping factors. In particular, we pay attention to the entire redshift range throughout reionization, instead of focusing on the end. Our measurements of both result in higher values compared to previous works, and we observe a turning point in the evolution of $C_I$ with the neutral fraction.

Finally, we use both the Madau (2017) methods and our models and fits to solve for the evolution of the mass-weighted ionization fraction and compare both to our simulation results. While the M99 model results in a >20% difference from all three simulations, our model fits the simulation results to within 5%. The low clumping factor used in Madau (2017) results in an earlier end of reionization, and it also has a >10% difference in the Thomson optical depth $\tau_e$ and a >20% difference in $\tau_{e,\text{reion}}$. Our model, on the other hand, matches both the
ionization history and the value of $\tau_e$ from the simulations much better, with a <5% difference from the simulations.

We thank the anonymous reviewer for asking thoughtful questions and providing suggestions to make our work more comprehensive. We thank Alexander Kaurov, Francois Lanusse, and Michelle Ntampaka for helpful discussions. We are grateful for the data provided by Kristian Finlator in his paper Finlator et al. (2012). N.C. is supported by the John Peoples Jr. Research Fellowship at Carnegie Mellon University. A.D. acknowledges the McWilliams Center for Cosmology for hosting his internship. H.T. acknowledges support from STScI grant HST-AR-15013.002-A. R.C. is supported in part by NASA grant 80NSSC18K1101. Simulations were run at the NASA Advanced Supercomputing (NAS) Center.

## ORCID iDs

Nianyi Chen https://orcid.org/0000-0001-6627-2533

Aristide Doussot https://orcid.org/0000-0002-4586-7586

Hy Trac https://orcid.org/0000-0001-6778-3861

Renyue Cen https://orcid.org/0000-0001-8531-9536

## References

Alvarez, M. A., Busha, M., Abel, T., & Wechsler, R. H. 2009, ApJL, 703, L167

Barkana, R., & Loeb, A. 2004, ApJ, 609, 474

Battaglia, N., Trac, H., Cen, R., & Loeb, A. 2013, ApJ, 776, 81

Becker, G. D., Bolton, J. S., Madau, P., et al. 2015, MNRAS, 447, 3402

Bolton, J. S., & Haehnelt, M. G. 2007, MNRAS, 382, 325

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015a, ApJ, 811, 140

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015b, ApJ, 803, 34

D'Aloisio, A., McQuinn, M., Maupin, O., et al. 2019, ApJ, 874, 154

Doussot, A., Trac, H., & Cen, R. 2019, ApJ, 870, 18

Finkelstein, S. L. 2016, PASA, 33, e037

Finlator, K., Keating, L., Oppenheimer, B. D., Davé, R., & Zackrisson, E. 2018, MNRAS, 480, 2628

Finlator, K., Oh, S. P., Özel, F., & Davé, R. 2012, MNRAS, 427, 2464

Furlanetto, S. R., Zaldarriaga, M., & Hernquist, L. 2004, ApJ, 613, 1

Glazer, D., Rau, M. M., & Trac, H. 2018, RNAAS, 2, 135

Gnedin, N. Y. 2014, ApJ, 793, 29

Gnedin, N. Y., & Ostriker, J. P. 1997, ApJ, 486, 581

Gorce, A., Douspis, M., Aghanim, N., & Langer, M. 2018, A&A, 616, A113

Greig, B., & Mesinger, A. 2017, MNRAS, 472, 2651

Haardt, F., & Madau, P. 2012, ApJ, 746, 125

Hui, L., & Haiman, Z. 2003, ApJ, 596, 9

Ishigaki, M., Kawamata, R., Ouchi, M., et al. 2015, ApJ, 799, 12

Ishigaki, M., Kawamata, R., Ouchi, M., et al. 2018, ApJ, 854, 73

Jeeson-Daniel, A., Ciardi, B., & Graziani, L. 2014, MNRAS, 443, 2722

Kaurov, A. A., & Gnedin, N. Y. 2014, ApJ, 787, 146

Kaurov, A. A., & Gnedin, N. Y. 2015, ApJ, 810, 154

Keating, L. C., Weinberger, L. H., Kulkarni, G., et al. 2020, MNRAS, 491, 1736

Kohler, K., Gnedin, N. Y., & Hamilton, A. J. S. 2007, ApJ, 657, 15

Kuhlen, M., & Faucher-Giguère, C.-A. 2012, MNRAS, 423, 862

La Plante, P., & Trac, H. 2016, ApJ, 828, 90

Livermore, R. C., Finkelstein, S. L., & Lotz, J. M. 2017, ApJ, 835, 113

Madau, P. 2017, ApJ, 851, 50

Madau, P., & Haardt, F. 2015, ApJL, 813, L8

Madau, P., Haardt, F., & Rees, M. J. 1999, ApJ, 514, 648

McQuinn, M., Oh, S. P., & Faucher-Giguère, C.-A. 2011, ApJ, 743, 82

Mesinger, A., Furlanetto, S., & Cen, R. 2011, MNRAS, 411, 955

Miralda-Escudé, J., Haehnelt, M., & Rees, M. J. 2000, ApJ, 530, 1

Nasir, F., & D'Aloisio, A. 2020, MNRAS, 494, 3080

Norman, M. L., Reynolds, D. R., So, G. C., Harkness, R. P., & Wise, J. H. 2015, ApJS, 216, 16

Pawlik, A. H., Schaye, J., & van Scherpenzeel, E. 2009, MNRAS, 394, 1812

Planck Collaboration, Adam, R., Aghanim, N., et al. 2016, A&A, 596, A108

Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6

Price, L. C., Trac, H., & Cen, R. 2016, arXiv:1605.03970

Robertson, B. E., Ellis, R. S., Furlanetto, S. R., & Dunlop, J. S. 2015, ApJL, 802, L19

Santos, M. G., Ferramacho, L., Silva, M. B., Amblard, A., & Cooray, A. 2010, MNRAS, 406, 2421

Semelin, B., Eames, E., Bolgar, F., & Caillat, M. 2017, MNRAS, 472, 4508

Shull, J. M., Harness, A., Trenti, M., & Smith, B. D. 2012, ApJ, 747, 100

Trac, H., & Cen, R. 2007, ApJ, 671, 1

Trac, H., Cen, R., & Loeb, A. 2008, ApJL, 689, L81

Trac, H., Cen, R., & Mansfield, P. 2015, ApJ, 813, 54

Trac, H., & Pen, U.-L. 2004, NewA, 9, 443

Wales, D. J., & Doye, J. P. K. 1997, JPCA, 101, 5111

Zahn, O., Lidz, A., McQuinn, M., et al. 2007, ApJ, 654, 12