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

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³ dark matter particles, 2048³ gas cells, and up to 12 billion adaptive rays in a 50 h⁻¹ 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}_{\mathrm{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³M) N-body simulation with 3072³ 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 × 10⁵ h⁻¹ M_⊙, 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 τ_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,

$$\begin{eqnarray}&&{f}_{\mathrm{esc}}(z)={f}_{8}{\left(\displaystyle \frac{1+z}{9}\right)}^{{a}_{8}},\end{eqnarray} \tag{ 1 }$$

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

Table 1 summarizes the parameters for the three RadHydro simulations. Here z_mid is the redshift at which 50% of the hydrogen is ionized, and Δ_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_esc, and reionization starts latest but ends earliest out of the three models. Sim 1 has f_esc(z) varying linearly with 1 + z and is an intermediate model. Sim 2 has f_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).

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),

$$\begin{eqnarray}&&\displaystyle \frac{{{dQ}}_{\mathrm{HII}}}{{dt}}=\displaystyle \frac{\langle {\dot{n}}_{\gamma }{\rangle }_{{\rm{V}}}}{\langle {n}_{{\rm{H}}}{\rangle }_{{\rm{V}}}}-\displaystyle \frac{{Q}_{\mathrm{HII}}}{{\bar{t}}_{\mathrm{rec}}},\end{eqnarray} \tag{ 2 }$$

where Q_HII is the volume filling factor of ionized hydrogen, 〈n_H〉_V is the volume-averaged total hydrogen number density, $\langle {\dot{n}}_{\gamma }{\rangle }_{{\rm{V}}}$ is the volume-averaged photon production rate, and the effective recombination time is given by

$$\begin{eqnarray}&&{\bar{t}}_{\mathrm{rec}}=\displaystyle \frac{1}{(1+Y/4X)\alpha ({T}_{0})\langle {n}_{{\rm{H}}}{\rangle }_{{\rm{V}}}C},\end{eqnarray} \tag{ 3 }$$

where X = 0.76 and Y = 0.24 are the hydrogen and helium mass fractions, respectively; C is the clumping factor; α is the recombination coefficient; and T₀ is the temperature of the IGM at mean density, fixed to be T₀ = 2 × 10⁴ 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.

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{eqnarray}\begin{array}{rcl}\langle {x}_{\mathrm{HII}}{\rangle }_{{\rm{M}}} & = & \displaystyle \frac{\displaystyle \sum {\rho }_{{\rm{H}},i}{x}_{\mathrm{HII},i}}{\displaystyle \sum {\rho }_{{\rm{H}},i}}=\displaystyle \frac{\displaystyle \sum {\rho }_{\mathrm{HII},i}}{\displaystyle \sum {\rho }_{{\rm{H}},i}}\\ & = & \displaystyle \frac{\langle {\rho }_{\mathrm{HII}}{\rangle }_{{\rm{V}}}}{\langle {\rho }_{{\rm{H}}}{\rangle }_{{\rm{V}}}}=\displaystyle \frac{\langle {n}_{\mathrm{HII}}{\rangle }_{{\rm{V}}}}{\langle {n}_{{\rm{H}}}{\rangle }_{{\rm{V}}}},\end{array}\end{eqnarray} \tag{ 4 }$$

where the summation is over all cells in the simulation, and we have used the relation ρ_H = m_Hn_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),

$$\begin{eqnarray}&&\displaystyle \frac{{{dn}}_{\mathrm{HII}}}{{dt}}={\rm{\Gamma }}{n}_{\mathrm{HI}}+{\gamma }_{\mathrm{coll}}{n}_{{\rm{e}}}{n}_{\mathrm{HI}}-\alpha (T){n}_{{\rm{e}}}{n}_{\mathrm{HII}}-3{{Hn}}_{\mathrm{HII}},\end{eqnarray} \tag{ 5 }$$

where n_HII is the physical number density of ionized hydrogen, n_e is the physical number density of free electrons, γ_coll is the collisional ionization rate, Γ is the photoionization rate, α 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 〈n_H〉_V, we obtain the global equation for the mass-weighted ionization fraction,

$$\begin{eqnarray}&&\displaystyle \frac{d\langle {x}_{\mathrm{HII}}{\rangle }_{{\rm{M}}}}{{dt}}=\displaystyle \frac{\langle {\rm{\Gamma }}{n}_{\mathrm{HI}}{\rangle }_{{\rm{V}}}}{\langle {n}_{{\rm{H}}}{\rangle }_{{\rm{V}}}}-\displaystyle \frac{\langle \alpha (T){n}_{{\rm{e}}}{n}_{\mathrm{HII}}\rangle }{\langle {n}_{{\rm{H}}}{\rangle }_{{\rm{V}}}},\end{eqnarray} \tag{ 6 }$$

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

Instead of the volume filling factor Q, we focus on the evolution of 〈x_HII〉_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_e to the density of ionized hydrogen n_HII, assuming helium is singly ionized, by

$$\begin{eqnarray}&&{n}_{{\rm{e}}}=\left(1+\displaystyle \frac{Y}{4X}\right){n}_{\mathrm{HII}},\end{eqnarray} \tag{ 7 }$$

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

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d\langle {x}_{\mathrm{HII}}{\rangle }_{{\rm{M}}}}{{dt}} & = & {C}_{{\rm{I}}}\langle {\rm{\Gamma }}{\rangle }_{{\rm{V}}}(1-\langle {x}_{\mathrm{HII}}{\rangle }_{{\rm{M}}})\\ & & -{C}_{{\rm{R}}}\langle \alpha (T){\rangle }_{{\rm{V}}}\langle {n}_{{\rm{H}}}{\rangle }_{{\rm{V}}}\left(1+\displaystyle \frac{Y}{4X}\right)\langle {x}_{\mathrm{HII}}{\rangle }_{{\rm{M}}}^{2},\end{array}\end{eqnarray} \tag{ 8 }$$

where C_I is the photoionization clumping factor (e.g., Kohler et al. 2007) that occurs due to spatial fluctuations in the radiation field,

$$\begin{eqnarray}&&{C}_{{\rm{I}}}=\displaystyle \frac{\langle {\rm{\Gamma }}{n}_{\mathrm{HI}}{\rangle }_{{\rm{V}}}}{\langle {\rm{\Gamma }}{\rangle }_{{\rm{V}}}\langle {n}_{\mathrm{HI}}{\rangle }_{{\rm{V}}}},\end{eqnarray} \tag{ 9 }$$

and C_R is the recombination clumping factor,

$$\begin{eqnarray}&&{C}_{{\rm{R}}}=\displaystyle \frac{\langle \alpha (T){n}_{\mathrm{HII}}{n}_{e}{\rangle }_{{\rm{V}}}}{\langle \alpha (T){\rangle }_{{\rm{V}}}\langle {n}_{\mathrm{HII}}{\rangle }_{{\rm{V}}}\langle {n}_{e}{\rangle }_{{\rm{V}}}}.\end{eqnarray} \tag{ 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 Γ 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, $\dot{{n}_{\gamma }}$. This quantity can be computed from galaxy luminosity functions and the escape fraction without running RT simulations. If we use the photon production rate $\dot{{n}_{\gamma }}$ in the place of the photoionization rate Γ, Equation (8) becomes

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d\langle {x}_{\mathrm{HII}}{\rangle }_{{\rm{M}}}}{{dt}} & = & \displaystyle \frac{\langle {\dot{n}}_{\gamma }{\rangle }_{{\rm{V}}}}{\langle {n}_{{\rm{H}}}{\rangle }_{{\rm{V}}}}\\ & & -{{\boldsymbol{C}}}_{{\rm{R}}}\langle \alpha ({\boldsymbol{T}}){\rangle }_{{\rm{V}}}\langle {n}_{{\rm{H}}}{\rangle }_{{\rm{V}}}\left(1+\displaystyle \frac{Y}{4X}\right)\langle {x}_{\mathrm{HII}}{\rangle }_{{\rm{M}}}^{2}.\end{array}\end{eqnarray} \tag{ 11 }$$

Notice that only the first term on the right-hand side changes, and there is no need of an ionization clumping factor C_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.

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_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 Δ = 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_Δ<50, C_Δ<100, and C_Δ<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.

The total hydrogen clumping factor C_H is calculated from all hydrogen gas, both neutral and ionized, from the simulation. It is defined as

$$\begin{eqnarray}&&{C}_{{\rm{H}}}=\displaystyle \frac{\langle {n}_{{\rm{H}}}^{2}{\rangle }_{{\rm{V}}}}{\langle {n}_{{\rm{H}}}{\rangle }_{{\rm{V}}}^{2}}.\end{eqnarray} \tag{ 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 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 Δ < 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.

Zoom In Zoom Out Reset image size

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 Δ 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

$$\begin{eqnarray}&&{C}_{{\rm{H}}}(z)=1+{C}_{0}{\left(\displaystyle \frac{1+z}{9}\right)}^{\alpha (z)},\end{eqnarray} \tag{ 13 }$$

where C₀ is a constant controlling the overall amplitude, and α is a linear function in redshift with α(z) = α₀ + α₁(z − 8).

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

$$\begin{eqnarray}&&\epsilon =\displaystyle \sum _{z={z}_{0}}^{{z}_{{\rm{n}}}}\displaystyle \frac{| {C}_{\mathrm{sim}}-{C}_{\mathrm{fit}}| }{{C}_{\mathrm{sim}}}.\end{eqnarray} \tag{ 14 }$$

Here we choose z₀ = 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 α₀ term dominates in the exponent, as α₁ in the linear term is relatively small, and high-Δ models have more negative α₀, 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.

Table 2.  C_H: Total Hydrogen Clumping Factor

		CH
Model	C0	α0	α1	δrms	${\delta }_{\max }$
CH,Δ<50	1.62	−2.27	−0.107	0.4%	2.1%
CH,Δ<100	2.02	−2.89	−0.079	0.5%	2.6%
CH,Δ<200	2.53	−3.55	−0.056	1.1%	3.9%

Download table as:  ASCIITypeset image

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 α 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:

$$\begin{eqnarray}&&{C}_{\mathrm{HII}}=\displaystyle \frac{\langle {n}_{\mathrm{HII}}^{2}{\rangle }_{{\rm{V}}}}{\langle {n}_{\mathrm{HII}}{\rangle }_{{\rm{V}}}^{2}}.\end{eqnarray} \tag{ 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_Shull, corresponding to our C_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 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 Δ = (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,

$$\begin{eqnarray}&&{C}_{\mathrm{HII}}(z,{x}_{\mathrm{HII}})=\left[1+{C}_{0}{\left(\displaystyle \frac{1+z}{9}\right)}^{\alpha (z)}\right]\langle {x}_{\mathrm{HII}}{\rangle }_{{\rm{M}}}^{{\beta }_{0}},\end{eqnarray} \tag{ 16 }$$

where α(z) = α₀ + α₁(z − 8), and C₀, α₀, α₁, and β₀ 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₀ in Equation (14) to be the redshift when 〈x_HII〉_M = 1% and z_n when 〈x_HII〉_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}_{\mathrm{HII}}\propto {x}_{\mathrm{HII}}^{-1.3}$, so the dependence on ionization fraction is not affected much by the density upper limit.

Table 4.  C_HII: Ionized Hydrogen Clumping Factor

			CHII
Model	C0	α0	α1	β0	δrms	${\delta }_{\max }$
CHII,Δ<50	2.07	−1.43	−0.213	−1.32	5.0%	9.7%
CHII,Δ<100	2.69	−1.92	−0.196	−1.33	2.2%	8.3%
CHII,Δ<200	3.48	−2.50	−0.185	−1.34	3.1%	21%

Download table as:  ASCIITypeset image

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 Δ, 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 Δ>100 density cut. Here C_HII starts off at a relatively high value of  ∼10³ 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.

Zoom In Zoom Out Reset image size

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 3 shows the dependence of C_HII on the density cuts. Similar to C_H, lower density cuts (e.g., Δ < 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.

Zoom In Zoom Out Reset image size

Comparing to previous works on C_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_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 Δ or x_HII. This means that we have included all cells outside halos, leading to a larger C_HII at high redshift when the ionization fraction is very inhomogeneous.

As mentioned in Section 2.3, instead of assuming a constant recombination coefficient, it is more accurate to take the spatial variation of α into account when calculating the recombination clumping factor. Hence, we also calculated C_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 〈α〉 in the denominator, we choose a fiducial value of α = α_B(2 × 10⁴ K), such that our renormalized recombination clumping factor C_rec is defined as

$$\begin{eqnarray}&&{C}_{\mathrm{rec}}\equiv \displaystyle \frac{\langle \alpha (T){n}_{{\rm{e}}}{n}_{\mathrm{HII}}\rangle }{{\alpha }_{{\rm{B}}}(2\times {10}^{4}{\rm{K}})\langle {n}_{{\rm{e}}}\rangle \langle {n}_{\mathrm{HII}}\rangle }.\end{eqnarray} \tag{ 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_rec differently if they want.

Figures 4 and 5 show our measurements of C_rec from the RadHydro simulations for three simulations and three density cuts, respectively. Same as C_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_HII to fit the recombination clumping factor,

$$\begin{eqnarray}&&{C}_{\mathrm{rec}}(z,{x}_{\mathrm{HII}})=\left[1+{C}_{0}{\left(\displaystyle \frac{1+z}{9}\right)}^{\alpha (z)}\right]\langle {x}_{\mathrm{HII}}{\rangle }_{{\rm{M}}}^{{\beta }_{0}},\end{eqnarray} \tag{ 18 }$$

where α(z) = α₀ + α₁(z − 8), and C₀, α₀, α₁, and β₀ are parameters we want to optimize.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Like in C_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_HII on redshift is similar for C_HII and C_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.

Table 6.  C_rec: Normalized Recombination Clumping Factor

			Crec
Model	C0	α0	α1	β0	δrms	${\delta }_{\max }$
Crec,Δ<50	2.96	−1.31	−0.310	−1.39	4.9%	14%
Crec,Δ<100	3.38	−1.50	−0.341	−1.41	2.0%	8.6%
Crec,Δ<200	3.57	−1.58	−0.354	−1.42	3.7%	19%

Download table as:  ASCIITypeset image

For comparison, we add the data of C_R from Finlator et al. (2012) to Figure 4. The data were originally in x_HII,VC_R and used α(10⁴ K) for the averaged recombination coefficient (note that their C_R is defined similarly to our C_rec, although the normalization is different), but to show a direct comparison with our results, we divide it by x_HII,V and rescale to α(2 × 10⁴ 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 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_HII,VC_R instead of applying an ionization fraction lower limit at x_HII = 0.95, as the evolution of the two values is close to each other.

To use Equation (8) with the ionization rate Γ, we also need to compute the ionization clumping factor C_I defined in Equation (9). The measurement of C_I was previously computed by Kohler et al. (2007) using a small (4 h⁻¹ Mpc) box. In their work, the focus was toward the very end of reionization due to the interest in Lyα 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_I throughout reionization.

In Figure 6, we show the evolution of C_I 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. 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 photonionization 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.

Zoom In Zoom Out Reset image size

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_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 × 10⁴ K and use C = 3 and ${C}_{\mathrm{Shull}}=2.9{\left((1+z)/6\right)}^{-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_Shull reproduces the Madau (2017) calculation. For comparison purposes, we also use our fitted clumping factor C_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_Shull, and C_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.

Zoom In Zoom Out Reset image size

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_rec. We can see from the plots that solving the M99 equations with C_rec results in later reionization compared to the simulations, and this is due to the fact that the M99 recombination term, which is linear in Q, overestimates the recombination rate and thus postpones the overall ionization of hydrogen. Now if we instead use the C = 3 and C_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_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_RQ, 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_RQ 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_rec in Equation (18) to solve for the mass-weighted ionization fraction. Here we choose the C_rec with Δ < 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_rec,Δ<50 instead will lead to slightly earlier reionization as the recombination rate becomes lower, and using C_rec,Δ<200 will lead to slightly later reionization. Note that C_rec,Δ<50 and C_rec,Δ<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 〈x_HII〉_M from the analytical model proposed in this work, solved with an ionizing photon production rate ${\dot{n}}_{\gamma }$ and a photoionization rate Γ. 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 Γ 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.

Zoom In Zoom Out Reset image size

In the above solutions, we have used $\dot{{n}_{\gamma }}$ and Γ directly from our simulations. However, for future applications, people might have different $\dot{{n}_{\gamma }}$ and Γ 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.

One goal of computing the ionization fraction of hydrogen through the analytical models is to infer the Thomson optical depth τ_e, which is an important observable for reionization and parameter for cosmological surveys. It is defined as

$$\begin{eqnarray}&&{\tau }_{{\rm{e}}}(z)={\sigma }_{T}{\int }_{0}^{z}{dz}^{\prime} \displaystyle \frac{{dt}}{{dz}^{\prime} }\langle {n}_{{\rm{e}}}(z^{\prime} ){\rangle }_{{\rm{V}}},\end{eqnarray} \tag{ 19 }$$

where the volume-averaged free electron number density,

$$\begin{eqnarray}\begin{array}{rcl}\langle {n}_{{\rm{e}}}{\rangle }_{{\rm{V}}} & = & \langle {n}_{\mathrm{HII}}{\rangle }_{{\rm{V}}}+\langle {n}_{\mathrm{HeII}}{\rangle }_{{\rm{V}}}+2\langle {n}_{\mathrm{HeIII}}{\rangle }_{{\rm{V}}},\\ & = & \langle {x}_{\mathrm{HII}}{\rangle }_{{\rm{M}}}\langle {n}_{{\rm{H}}}{\rangle }_{{\rm{V}}}+\langle {x}_{\mathrm{HeII}}{\rangle }_{{\rm{M}}}\langle {n}_{\mathrm{He}}{\rangle }_{{\rm{V}}}\\ & & +2\langle {x}_{\mathrm{HeIII}}{\rangle }_{{\rm{M}}}\langle {n}_{\mathrm{He}}{\rangle }_{{\rm{V}}},\end{array}\end{eqnarray} \tag{ 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 τ_e using different models for the ionization fractions and compare how the choice of reionization models affects the value of τ_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 〈x_HII〉_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 reionzation down to redshift z = 0.

Table 8 shows the optical depths from four different models in comparison with the τ_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 τ_e. With the two models following M99, the values of τ_e are systematically higher by 5%–15% due to the earlier reionization of the gas, as already shown in Figure 7.

Table 8.  Thomson Optical Depth

			τe
Model	M99, C = 3	M99 + Shull12	Our Model w/ Γ	Our Model w/ ${\dot{n}}_{\gamma }$	RadHydro
Sim 0	0.065	0.067	0.059	0.060	0.060
Sim 1	0.066	0.069	0.059	0.060	0.060
Sim 2	0.067	0.070	0.061	0.059	0.060

Download table as:  ASCIITypeset image

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 τ_e measurement from Planck Collaboration et al. (2020) of 0.054 ± 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 τ_e calculated with our model is well within the tolerance of current observations.