Cosmic Reionization May Still Have Started Early and Ended Late: Confronting Early Onset with Cosmic Microwave Background Anisotropy and 21 cm Global Signals

It is a usual approximation that PPR is proportional to the star formation rate density (SFRD: SFR per co-moving volume, in units of M_⊙ s⁻¹ cMpc⁻³). Rigorously, this becomes true only if t_*, the lifetime of stars that we take as a constant for a given stellar species (e.g., Population II stars), is infinitesimal compared to the ionization time, i.e., x/(dx/dt). Even more rigorously, for PPR at time t to be the ionization rate at t, the photon travel time inside H ii regions should also be very small. In this paper, we take this usual assumption, PPR ∝ SFRD, which has been used extensively in semi-analytical modeling. The ionization rate ξ (in units of s⁻¹) will then be given by

$$\begin{eqnarray}&&\xi ={f}_{\mathrm{esc}}\mathrm{PPR}={f}_{\mathrm{esc}}{N}_{\mathrm{ion}}\displaystyle \frac{\mathrm{SFRD}}{{m}_{b}{n}_{b,0}},\end{eqnarray} \tag{ 1 }$$

where f_esc, N_ion, m_b , and n_b,0 are the ionizing-photon escape fraction out of halos, the total number of ionizing photons emitted per stellar baryon during the lifetime of a star, the baryon mass, and the average co-moving baryon number density, respectively.

We also modify the typical assumption of semi-analytical calculations to incorporate a nonzero duty cycle of star formation. Let us first briefly review the typical assumption in semi-analytical calculations: PPR and SFRD are assumed proportional to df_coll/dt, the growth rate of the halo-collapsed fraction f_coll, such that

$$\begin{eqnarray}\begin{array}{rcl}\xi (t) & = & {\xi }_{{\rm{a}}}(t)={f}_{* }{f}_{\mathrm{esc}}\displaystyle \frac{{N}_{\mathrm{ion}}}{{t}_{* }}{\displaystyle \int }_{t-{t}_{* }}^{t}{dt}^{\prime} \,\displaystyle \frac{{{df}}_{\mathrm{coll}}}{{dt}^{\prime} }\\ & \simeq & {f}_{* }{f}_{\mathrm{esc}}{N}_{\mathrm{ion}}\displaystyle \frac{{{df}}_{\mathrm{coll}}}{{dt}}(t)={f}_{\gamma }\displaystyle \frac{{{df}}_{\mathrm{coll}}}{{dt}}(t),\end{array}\end{eqnarray} \tag{ 2 }$$

where ξ_a denotes the ionization rate due to newly "accreted" matter, f_* is the star formation efficiency in such an episode (e.g., Loeb & Barkana 2001; Furlanetto 2006), and f_γ ≡ f_* f_esc N_ion is the proportionality coefficient between ξ and df_coll/dt. The second equality in Equation (2) is a rigorous form that considers only the newly accreted matter, whose maximal lookback time in the integration should be t_*. Therefore, we can denote this assumption as the "mass-accretion dominated star formation" scenario (MADSF). In the limit of infinitesimal t_* or slowly varying f_coll, the last approximation becomes valid. Because df_coll in Equation (2) is the amount of matter that has been newly accreted into halos per time interval dt, considering ξ_a(t) only is identical to having a null duty cycle such that any gas that has once formed stars can never form stars again. Note also that df_coll does not include halos whose masses have increased by a merger event at t, even though in nature "wet merger" events occur quite frequently. Because this null duty cycle assumption is somewhat extreme, we need to consider a more general case of nonzero duty cycle.

We quantify the nonzero duty cycle as the time fraction that stars emit radiation, f_DC ≡ t_*/(t_* + t_dorm), where t_dorm is the duration that a stellar baryon remains dormant after the death of the host star until a new star formation episode occurs. The same account on t_dorm has been addressed by Mirocha et al. (2018). Both t_* and t_dorm are average quantities, and consequently f_DC will be a global parameter. Ergodicity is assumed such that this temporal duty cycle is the same as the spatially averaged fraction of gas that has once resided inside a star and now resides in a post-generation star. In this work, we further assume that f_DC is constant over time for a given stellar species. With nonzero f_DC, we obtain the ionization rate ξ_d due to "duty cycle":

$$\begin{eqnarray}\begin{array}{rcl}{\xi }_{{\rm{d}}}(t) & = & {f}_{* }{f}_{\mathrm{esc}}{f}_{\mathrm{DC}}{f}_{\mathrm{coll}}(t-{t}_{* })\displaystyle \frac{{N}_{\mathrm{ion}}}{{t}_{* }}\\ & \simeq & {f}_{* }{f}_{\mathrm{esc}}{f}_{\mathrm{DC}}{f}_{\mathrm{coll}}(t)\displaystyle \frac{{N}_{\mathrm{ion}}}{{t}_{* }}={f}_{* ,{\rm{d}}}{f}_{\mathrm{esc}}{f}_{\mathrm{coll}}(t)\displaystyle \frac{{N}_{\mathrm{ion}}}{{t}_{* }},\end{array}\end{eqnarray} \tag{ 3 }$$

where only "old" gas at lookback time t_* is allowed to regenerate stars, and the approximation holds when t_* is treated as infinitesimal. The last equality absorbs f_DC into the effective f_* of halo gas, f_*,d. Equation (3) obviously does not account for any newly accreted gas. Note also that Equation (3) holds due to ergodicity: f_*,d is a spatially averaged quantity over many galaxies, which will be identical to the time average of star formation episodes on a single galaxy. We denote a scenario based on ξ = ξ_d as the "all-halo replenished star formation" scenario (AHRSF).

The most generic form for ξ should implement both contributions from old gas and newly accreted gas, or combining the episodes of MADSF and AHRSF, such that

$$\begin{eqnarray}&&\xi ={\xi }_{{\rm{d}}}+{\xi }_{{\rm{a}}}={f}_{* ,{\rm{d}}}{f}_{\mathrm{esc}}{f}_{\mathrm{coll}}\displaystyle \frac{{N}_{\mathrm{ion}}}{{t}_{* }}+{f}_{* ,{\rm{a}}}{f}_{\mathrm{esc}}{N}_{\mathrm{ion}}\displaystyle \frac{{{df}}_{\mathrm{coll}}}{{dt}},\end{eqnarray} \tag{ 4 }$$

where f_* of newly accreted gas is denoted by f_*,a, to be distinguished from f_*,d. The condition f_*,a = 0, corresponds to the case where newly accreted gas waits longer than t_* to start star formation activity, and practically there is no theoretical constraint on the relative strengths of f_*,a over f_*,d. With Equation (4), both cases of ξ ∝ f_coll/t_* and ξ ∝ df_coll/dt can be accommodated in terms of special cases of this generic form with {f_*,d ≠ 0, f_*,a = 0} and {f_*,d = 0, f_*,a ≠ 0}, respectively. Usually, the relation ξ ∝ f_coll/t_* has been used in numerical radiation transfer simulations (e.g., Iliev et al. 2007) and the relation ξ ∝ df_coll/dt in semi-analytical calculations (e.g., Loeb & Barkana 2001; Furlanetto 2006). In this paper, we only consider these two special cases, denoting the former by "F" and the latter by "dF" as the nomenclature for star formation scenarios. Because we use f_*,a and f_*,d mutually exclusively for these special cases, we drop the subscripts "a" (accretion) and "d" (duty cycle) from f_* for simplicity.

As long as cosmology is fixed, the evolution of the volume ionized fraction will be uniquely determined for a given f_γ in MADSF and for a given g_γ ≡ f_* f_esc N_ion/(t_*/10 Myr) in AHRSF, because df_coll/dt and f_coll are determined solely by cosmological parameters. Therefore, reionization models with MADSF and AHRSF will be parameterized by f_γ (e.g., Furlanetto 2006) and g_γ (e.g., Ahn et al. 2012), respectively. Of course, f_γ and g_γ need not be constant in time. f_esc, f_*, t_*, and N_ion can be time-varying individually and collectively (in terms of f_γ and g_γ ). Allowing f_γ and g_γ to change in time can yield even more variants of reionization scenarios, especially by reshaping dx/dz. For simplicity, we do not explore this possibility. Nevertheless, future CMB observations will become more accurate, and models inferred from the data could not be matched well by the constant values of f_γ and g_γ . In such a case, time-varying f_γ and g_γ may have to be considered in our models.

We consider three types of models: (1) the vanilla model, (2) the self-regulated model type I (SRI), and (3) the self-regulated model type II (SRII). From type (1) to type (3), these models become progressively sophisticated in the physical processes considered: the vanilla model does not consider any feedback, SRI considers the photo-heating feedback, and SRII considers both photo-heating and LW feedback effects. In each model, we allow both cases of ξ ∝ df_coll/dt and ξ ∝ f_coll/t_* that were described in Section 2.1. In all of these models, the source term is the ionization rate, and the sink term is the recombination rate, such that the change rate of the global (volume) ionization fraction x is

$$\begin{eqnarray}&&\displaystyle \frac{{dx}}{{dt}}=\xi -\alpha {{Cn}}_{{\rm{e}}}x,\end{eqnarray} \tag{ 5 }$$

where α is the hydrogen recombination coefficient (in units of cm³ s⁻¹), C is the average clumping factor, and n_e = n_H + n_He is the (proper) number density of electrons inside H ii regions (e.g., Furlanetto 2006; Iliev et al. 2007). In this work, we use the case B recombination coefficient for α, or α = α_B(T = 10⁴ K) = 2 × 10⁻¹³ cm³ s⁻¹, and adopt a fitting formula for C given by

$$\begin{eqnarray}&&C=\max (3,17.6\exp [-0.1z+0.0011{z}^{2}]),\end{eqnarray} \tag{ 6 }$$

which is a conservative combination of work by Iliev et al. (2005), Pawlik et al. (2009), and So et al. (2014). While it is possible that the numerical resolution limit of previous numerical simulations may have led to underestimation of C (Mao et al. 2020), we do not consider this possibility in this paper. The model variance is mainly caused by the variance in ξ in Equation (5), which will be described in the following subsections.

Physical parameters governing ξ, such as f_esc, f_*, N_ion, etc., would depend on halo properties. One of the crucial halo properties is the virial temperature T_vir. The natural borderline between MHs and ACHs is the temperature endpoint of the atomic line cooling (if dominated by hydrogen Lyα line cooling), ∼ 10⁴ K. ACHs have T_vir > 10⁴ K, and MHs have T_vir < 10⁴ K. This classification can be cast into the mass criterion

$$\begin{eqnarray}&&M\gtrless {10}^{8}\,{h}^{-1}{M}_{\odot }{\left[1.98\left\{{{\rm{\Omega }}}_{m}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}{\left(1+z\right)}^{-3}\right\}\left(\displaystyle \frac{1+z}{10}\right)\right]}^{-\tfrac{3}{2}},\end{eqnarray} \tag{ 7 }$$

where we assume the mean molecular weight μ = 0.6 for fully ionized gas (Barkana & Loeb 2001). One important subtlety is that this redshift-dependent mass criterion is not commonly respected in numerical radiative transfer (RT) simulations, because the minimum mass of halos resolvable in accompanying N-body or uniform-grid simulations tends to be constant in time. ³ Therefore, many RT simulations take a constant-mass classification scheme, such that ACHs and MHs correspond to halos with M > M_th and M < M_th, respectively, with a constant-mass threshold M_th. When a simulation domain is increased to a few ∼10 cMpc, the numerical resolution limit quickly reaches ∼10⁸ M_⊙, a common value of M_th for such classification. For example, in many large-volume simulations, MHs are ignored, and all halos that are numerically resolved, or, e.g., those with M > 10⁸ M_⊙, are taken as ACHs. However, this scheme even loses track of those halos with T_vir > 10⁴ K and M < 10⁸ M_⊙ at z ≳ 18. Therefore, one should be wary of the negligence when the high-redshift (z ≳ 20) astrophysics, for example, the early stage of cosmic reionization, is investigated with such a scheme. It may ignore not only MHs but also a substantial amount of ACHs, as long as we believe in the constant-temperature criterion as a natural distinction.

In order to have a fair comparison of the semi-analytical models to the previous numerical RT simulation results based on this constant-mass criterion for ACH/MH classification, we adopt a constant M_th (=10⁸ M_⊙) in this work. Comparison will be made to numerical simulations that (1) further split ACHs into low-mass and high-mass species with a constant-mass classification (Section 2.2.2), and (2) cover the full dynamic range of halos, namely MHs, LMACHs, and HMACHs but again with a constant-mass classification (Section 2.2.3).

2.2.1. Vanilla Model

This model assumes a very simplified form for PPR: PPR is typically assumed to be proportional to df_coll/dt, and no feedback effect on star formation is considered. The essential ingredient is not the relation ξ ∝ df_coll/dt but the lack of feedback, and thus we also allow the relation ξ ∝ f_coll/t_*. We therefore have

$$\begin{eqnarray}&&\xi ={\xi }_{{\rm{a}}}=\sum _{i}{f}_{* }^{\left(i\right)}{f}_{\mathrm{esc}}^{\left(i\right)}{N}_{\mathrm{ion}}^{\left(i\right)}\displaystyle \frac{{{df}}_{\mathrm{coll}}^{\left(i\right)}}{{dt}}=\sum _{i}{f}_{\gamma }^{\left(i\right)}\displaystyle \frac{{{df}}_{\mathrm{coll}}^{\left(i\right)}}{{dt}}\end{eqnarray} \tag{ 8 }$$

for ξ ∝ df_coll/dt ("dF") and

$$\begin{eqnarray}&&\xi ={\xi }_{{\rm{d}}}=\sum _{i}{f}_{* }^{\left(i\right)}{f}_{\mathrm{esc}}^{\left(i\right)}{f}_{\mathrm{coll}}^{\left(i\right)}\displaystyle \frac{{N}_{\mathrm{ion}}^{\left(i\right)}}{{t}_{* }^{\left(i\right)}}\end{eqnarray} \tag{ 9 }$$

for ξ ∝ f_coll/t_* ("F"). These equations are generalized from Equations (2) and (3) with summation to accommodate different halo species i, with the superscript "(i)" denoting physical quantities of the halo species i. The simplicity of Equations (8) and (9) is the essence of vanilla models that ignore any feedback effects.

2.2.2. Self-regulated Model Type I

The work by Iliev et al. (2007) is among the first 3D RT simulations of self-regulated reionization based on multispecies halo stars, but with radiation sources restricted to ACHs. This simulation starts with classifying ACHs into two mass categories, namely (1) the HMACHs and (2) the LMACHs, with M ≳ 10⁹ M_⊙ and 10⁸ M_⊙ ≲ M ≲ 10⁹ M_⊙, respectively. Again, even though it is more natural to classify halos in terms of T_vir, to make a direct comparison to Iliev et al. (2007), we adopt this convention in this work. While the accurate boundary does not exist, LMACHs defined this way roughly correspond to those halos that are subject to the "Jeans-mass filtering": if these halos are formed inside regions that have been already ionized, accretion of baryonic gas will not be efficient enough to form stars inside due to the high temperature (T ≳ 10⁴ K) of the regions (Efstathiou 1992; Shapiro et al. 1994; Thoul & Weinberg 1996; Navarro & Steinmetz 1997; Gnedin & Hui 1998; Gnedin 2000; Dijkstra et al. 2004). The suppression of star formation is likely to be not as abrupt as assumed in Iliev et al. (2007) but rather gradual in halo mass (Efstathiou 1992; Navarro & Steinmetz 1997; Dijkstra et al. 2004). Nevertheless, in this work we simply adopt the self-regulation scheme of Iliev et al. (2007).

We denote such a type as SRI. Because HMACHs are unaffected by the feedback and LMACHs are suppressed inside H ii regions, the ionization rate will be given by

$$\begin{eqnarray}&&\xi ={f}_{* }^{{\rm{H}}}{f}_{\mathrm{esc}}^{{\rm{H}}}{N}_{\mathrm{ion}}^{{\rm{H}}}\displaystyle \frac{{{df}}_{\mathrm{coll}}^{{\rm{H}}}}{{dt}}+{f}_{* }^{{\rm{L}}}{f}_{\mathrm{esc}}^{{\rm{L}}}{N}_{\mathrm{ion}}^{{\rm{L}}}\displaystyle \frac{{{df}}_{\mathrm{coll}}^{{\rm{L}}}}{{dt}}(1-{x}^{\eta })\end{eqnarray} \tag{ 10 }$$

for ξ ∝ df_coll/dt ("dF") and

$$\begin{eqnarray}&&\xi ={f}_{* }^{{\rm{H}}}{f}_{\mathrm{esc}}^{{\rm{H}}}{f}_{\mathrm{coll}}^{{\rm{H}}}\displaystyle \frac{{N}_{\mathrm{ion}}^{{\rm{H}}}}{{t}_{* }^{{\rm{H}}}}+{f}_{* }^{{\rm{L}}}{f}_{\mathrm{esc}}^{{\rm{L}}}{f}_{\mathrm{coll}}^{{\rm{L}}}\displaystyle \frac{{N}_{\mathrm{ion}}^{{\rm{L}}}}{{t}_{* }^{{\rm{L}}}}(1-{x}^{\eta })\end{eqnarray} \tag{ 11 }$$

for ξ ∝ f_coll/t_* ("F"), where superscripts "H" and "L" denote HMACH and LMACH, respectively, and 0 < η < 1. The reason why such an amplified suppression term (1 − x^η ) is used instead of (1 − x) is that LMACHs are clustered more strongly inside H ii regions than in neutral regions. The value η = 0.1 is the empirical one found in 3D numerical simulations by Iliev et al. (2007), such that x^η roughly estimates the mass fraction of LMACHs occupied by H ii regions when the global ionization fraction is x. The effect of such a small value of η is to make star formation in LMACHs, and consequently the reionization history, more strongly regulated than the unrealistic case with η = 1 where LMACHs are uniformly distributed in space. The assumption for the Jeans-mass filtering in Equation (11) is that even those halos that have collapsed earlier cannot host new star formation episodes. Even though one can generalize ξ into a combined form of Equations (10) and (11), as we mentioned in Section 2.1, we restrict our models to these two categories of "dF" and "F."

2.2.3. Self-regulated Model Type II

The work by Ahn et al. (2012) is a unique 3D RT simulation of reionization in that (1) a full dynamic range of halos, from MHs to HMACHs, is treated as radiation sources, (2) both the Jeans-mass filtering and the LW feedback are considered, and (3) the simulation box is large (∼150 Mpc) enough to provide a reliable statistical significance. This simulation does not neglect MHs as most other large-box simulations do. MHs are subject both to the Jeans-mass filtering and the LW feedback. Jeans-mass filtering of MHs is obvious due to the smallness of the the virial temperature (T ≲ 10⁴ K). The LW feedback occurs due to the fact that H₂ is the main cooling agent in the primordial environment, and MHs are usually formed first in the primordial environment. Stars born in this environment will be Population III stars. Even inside MHs, after a few episodes of star formation, the chemical environment can gain metallicity beyond the critical value Z ≃ 10⁻³. However, dynamical feedback from supernova explosion inside MHs is believed to be very destructive (Greif et al. 2007; Yoshida et al. 2007), such that it may take longer than, e.g., the halo merger time for post-generation star formation episodes to occur in the same MH. It is possible that the supernova feedback in the massive MHs, which remained neutral even after photoionization from stars inside, could have been confined inside the halo and led to the next episode of star formation (Whalen et al. 2008). However, the dominant contribution to the number of MHs is from the least massive ones, and so it is appropriate to assume the destructive feedback. If one assumed the most destructive feedback effect of the first episode of star formation inside MHs, then it would be equivalent to assuming that only the newly forming MHs form stars. This assumption was taken in Ahn et al. (2012), which we implement in our modeling here as well.

We denote such a type as SRII. The ionization rate will be given by

$$\begin{eqnarray}\begin{array}{rcl}\xi & = & {f}_{* }^{{\rm{H}}}{f}_{\mathrm{esc}}^{{\rm{H}}}{N}_{\mathrm{ion}}^{{\rm{H}}}\displaystyle \frac{{{df}}_{\mathrm{coll}}^{{\rm{H}}}}{{dt}}+{f}_{* }^{{\rm{L}}}{f}_{\mathrm{esc}}^{{\rm{L}}}{N}_{\mathrm{ion}}^{{\rm{L}}}\displaystyle \frac{{{df}}_{\mathrm{coll}}^{{\rm{L}}}}{{dt}}(1-{x}^{\eta })\\ & & +\,\displaystyle \frac{{f}_{\mathrm{esc}}^{{\rm{M}}}{M}_{\mathrm{III}}{N}_{\mathrm{ion}}^{{\rm{M}}}}{\mu {m}_{{\rm{H}}}{n}_{b,0}}\displaystyle \frac{{{dn}}^{{\rm{M}}}}{{dt}}\left[1-\min \left\{\max \left({x}_{\mathrm{LW}},x\right),1\right\}\right]\end{array}\end{eqnarray} \tag{ 12 }$$

for ξ ∝ df_coll/dt ("dF") and

$$\begin{eqnarray}\begin{array}{rcl}\xi & = & {f}_{* }^{{\rm{H}}}{f}_{\mathrm{esc}}^{{\rm{H}}}{f}_{\mathrm{coll}}^{{\rm{H}}}\displaystyle \frac{{N}_{\mathrm{ion}}^{{\rm{H}}}}{{t}_{* }^{{\rm{H}}}}+{f}_{* }^{{\rm{L}}}{f}_{\mathrm{esc}}^{{\rm{L}}}{f}_{\mathrm{coll}}^{{\rm{L}}}\displaystyle \frac{{N}_{\mathrm{ion}}^{{\rm{L}}}}{{t}_{* }^{{\rm{L}}}}(1-{x}^{\eta })\\ & & +\,\displaystyle \frac{{f}_{\mathrm{esc}}^{{\rm{M}}}{M}_{\mathrm{III}}{N}_{\mathrm{ion}}^{{\rm{M}}}}{\mu {m}_{{\rm{H}}}{n}_{b,0}}\displaystyle \frac{{{dn}}^{{\rm{M}}}}{{dt}}\left[1-\min \left\{\max \left({x}_{\mathrm{LW}},x\right),1\right\}\right]\end{array}\end{eqnarray} \tag{ 13 }$$

for ξ ∝ f_coll/t_* ("F"), where the superscript M denotes MH, M_III is the mass of Population III stars per MH, n^M is the co-moving number density of MHs, n_b,0 ≃ 2 × 10⁻⁷ cm⁻³ is the co-moving number density of baryons, μ = 1.22 is the mean molecular weight of MH gas, and x_LW is the LW intensity J_LW normalized by the threshold intensity J_LW,th given by

$$\begin{eqnarray}&&{x}_{\mathrm{LW}}={J}_{\mathrm{LW}}/{J}_{\mathrm{LW},\mathrm{th}},\end{eqnarray} \tag{ 14 }$$

to implement the suppression of star formation inside MHs (see Equations (12) and (13)) in a similar fashion with Ahn et al. (2012).

Because we only take newly forming MHs as sources, the last terms in Equations (12) and (13) are identical. Having $\xi \propto {{dn}}^{{\rm{M}}}/{dt}$ instead of $\xi \propto {{df}}_{\mathrm{coll}}^{{\rm{M}}}/{dt}$ for MHs is to accommodate the tendency found in numerical simulations of Population III star formation: Population III stars under the primordial environment will form mostly in isolation (Abel et al. 2000; Bromm et al. 2002) or as a few binary systems at most (Turk et al. 2009; Stacy et al. 2010), and the total mass of Population III stars in MHs is determined by the atomic physics (Hirano et al. 2014, 2015) and is not strongly dependent on the mass of MHs. The actual mass of Population III stars can vary substantially according to physical properties, such as the angular momentum, the local LW intensity, and the mass accretion rate of their host MHs (Hirano et al. 2014, 2015), and thus M_III should be taken as the average mass of Population III stars per MH.

We limit the MH mass to 10⁵ M_⊙ ≤ M ≤ 10⁸ and use the Sheth–Tormen mass function (Sheth & Tormen 1999) based on the typical linear matter density perturbation obtained from the linear Boltzmann solver CAMB (Lewis et al. 2000) to calculate n ^M. Our choice of the minimum mass of MHs, ${M}_{\min }={10}^{5}\,{M}_{\odot }$, is also used in Ahn et al. (2012) and is indeed reasonable due to the following reasons. This value is about half of the Jeans mass of the neutral IGM when the baryon–dark matter streaming velocity is considered (Tseliakhovich et al. 2011). The actual minimum mass to host Population III stars could be somewhat larger than this Jeans mass and also redshift-dependent (see, e.g., Figure 1 of Glover 2013; Hirano et al. 2015). Even though some claim that ${M}_{\min }\gtrsim$ a few 10⁶ M_⊙ and thus the impact of MH stars on cosmic reionization is negligible (e.g., Kimm et al. 2017), other high-resolution numerical simulations (e.g., Hirano et al. 2015) find that massive Population III stars are hosted mostly by halos in the mass range ${M}_{\min }\simeq [{10}^{5},{10}^{6}]\,{M}_{\odot }$.

We also note that the ionization history x(t) would depend on ${M}_{\min }$ (or similarly on M_III) more weakly than J_LW,th and f_esc, and therefore determining an accurate value of ${M}_{\min }$ would not be too crucial. This is because star formation in MHs, during the time when MH stars are the dominant radiation sources, is regulated in a way to maintain J_LW ≃ J_LW,th (Ahn et al. 2012; see also Section 3.1). If ${M}_{\min }$ had been larger than our fiducial value and thus n^M had been smaller, then MH stars would have produced less ionizing and LW radiation in the beginning and drive resulting suppression weaker (or (1 − x_LW) larger). In this case, because of reduced suppression, MH star formation will soon be expedited until J_LW reaches J_LW,th and produce an x(t) evolution similar to the fiducial case thereafter. Similarly, any additional change in n^M due to the baryon–dark matter streaming effect is likely to be unimportant in determining x(t).

The way we implement the LW feedback as a multiplicative factor (1 − x_LW) in Equations (12) and (13) roughly follows the work of Yoshida et al. (2003) and O'Shea & Norman (2008), where they find a gradual increase in ${M}_{\min }$, the minimum mass of halos that can form stars, as J_LW increases if J_LW < J_LW,th. O'Shea & Norman (2008) also find that when J_LW ≥ J_LW,th ∼ 0.1 × 10⁻²¹ erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹, ${M}_{\min }$ and ${T}_{\mathrm{vir},\min }$ (the minimum virial temperature of halos that can form stars) jump to those of ACHs such that star formation in MHs are fully suppressed. Yoshida et al. (2003) find practically the same result, namely the full suppression of star formation inside MHs when J_LW ≥ J_LW,th ∼ 0.1 × 10⁻²¹ erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹ based on a combination of semi-analytical analysis and numerical simulation. Therefore, the exact functional form of the suppression is not important once J_LW reaches J_LW,th and the increasing number of MHs afterward try to produce more photons but fail to do so due to the self-regulation. As long as the condition that MHs become fully devoid of star formation when the condition J_LW ≥ J_LW,th is met, the formalism will correctly predict the self-regulation by LW feedback. This is exactly how the LW feedback is implemented in Equations (12)–(14). Other work (e.g., Mirocha et al. 2018; Qin et al. 2021) implementing the LW feedback on MH stars and forecasting the 21 cm background does not usually take this approach, which will be discussed in detail in Section 3.3.

The self-regulation of MH stars is predominantly governed by the LW feedback. In all of the SRII models we tested (see the detailed model parameters in Section 3.1), x_LW > x, and thus the regulation factor $(1-\min \{\max ({x}_{\mathrm{LW}},x),1\})$ in Equations (12) and (13) is practically identical to $(1-\min \{{x}_{\mathrm{LW}},1\})$. We also note that we do not use biased suppression of star formation for MH stars, as quantified by η for LMACHs in Equations (10)–(13). This is because star formation inside MHs is likely to occur much more diffusively than that inside ACHs. This tendency is indeed observed in the numerical simulation by Ahn et al. (2012): ionized regions generated by MH stars are almost uniformly spread in space, in contrast to those by ACH stars (Figure 2 of Ahn et al. 2012). Therefore, we simply take an unbiased regulation factor $(1-\min \{{x}_{\mathrm{LW}},1\})$ for MH stars.

We take ${N}_{\mathrm{ion}}^{{\rm{M}}}=$ 50,000, which is a value suitable for very massive stars (M_* ≳ 100 M_⊙). As was noted in Fialkov et al. (2013), quantifying suppression of SFR by LW intensity is not very straightforward when J_LW < J_LW,th, if, e.g., one considers the temporal evolution of LW intensity during halo formation. Our suppression scheme described by Equations (12)–(14) could instead perfectly mimic the full suppression by any threshold intensity J_LW,th, and our ignorance of any other details is parameterized by the value of J_LW,th. How J_LW is evaluated is described in Section 2.3.2.

There are a few radiation backgrounds that determine the reionization history and the 21 cm background. Any unprocessed background intensity (erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹) at observing frequency ν and redshift z is given by

$$\begin{eqnarray}&&{J}_{\nu }=\displaystyle \frac{c}{4\pi }{\left(1+z\right)}^{3}{\int }_{z}^{\infty }\displaystyle \frac{h\nu ^{\prime} {{ \mathcal N }}_{\nu ^{\prime} }(z^{\prime} )}{(1+z^{\prime} )H(z^{\prime} )}{{\rm{e}}}^{-{\tau }_{\nu }}{dz}^{\prime} ,\end{eqnarray} \tag{ 15 }$$

where h is the Planck constant, ${{ \mathcal N }}_{\nu ^{\prime} }(z^{\prime} )$ is the photon-number luminosity density (s⁻¹ Hz⁻¹ cMpc⁻³) at source frequency $\nu ^{\prime}$ and redshift $z^{\prime}$, and τ_ν is the optical depth from $z^{\prime}$ to z at ν. The soft-UV, LW, and X-ray backgrounds are all unprocessed types and thus given by Equation (15). The Lyα background is a processed type, and thus is not given by Equation (15); an appropriate description will be given instead in Section 2.3.4.

2.3.1. H-ionization by UV Background

2.3.2. H₂-dissociation by Lyman–Werner Background

The frequency-averaged, global LW intensity is given by

$$\begin{eqnarray}\begin{array}{rcl}{J}_{\mathrm{LW}}(z) & = & {\left\langle \displaystyle \frac{c}{4\pi }{\left(1+z\right)}^{3}{\displaystyle \int }_{z}^{\infty }\displaystyle \frac{h\nu ^{\prime} {{ \mathcal N }}_{\nu ^{\prime} }(z^{\prime} )}{(1+z^{\prime} )H(z^{\prime} )}{e}^{-{\tau }_{\nu }}{dz}^{\prime} \right\rangle }_{\nu }\\ & = & \displaystyle \frac{c}{4\pi }{\left(1+z\right)}^{3}{\displaystyle \int }_{z}^{\infty }\displaystyle \frac{{\left\langle h\nu ^{\prime} {{ \mathcal N }}_{\nu ^{\prime} }(z^{\prime} )\right\rangle }_{\nu ^{\prime} }}{(1+z^{\prime} )H(z^{\prime} )}{f}_{\mathrm{mod}}(z,z^{\prime} ){dz}^{\prime} ,\end{array}\end{eqnarray} \tag{ 16 }$$

where ${\left\langle \right\rangle }_{\nu }$ and ${\left\langle \right\rangle }_{\nu ^{\prime} }$ are frequency-averages in the observed band at z and in the emitted band at $z^{\prime}$, respectively, $\nu ^{\prime}$ is limited below the Lyman limit (LL, $h\nu ^{\prime} \lt 13.6\,\mathrm{eV}$), and ${f}_{\mathrm{mod}}$ is the "picket-fence modulation factor" that accounts for the trimming of bands of radiation from a source at $z^{\prime}$ due to the redshifting of continuum into Lyman resonance lines (Ahn et al. 2009). As seen in Equation (16), ${f}_{\mathrm{mod}}$ replaces the attenuation factor ${e}^{-{\tau }_{\nu }}$ and is given approximately by

$$\begin{eqnarray}&&\begin{array}{l}{f}_{\mathrm{mod}}(z,z^{\prime} )\\ =\,\left\{\begin{array}{ll}1.7\exp \left[-{\left(\tfrac{{r}_{\mathrm{cMpc}}}{116.29\alpha }\right)}^{0.68}\right]-0.7 & \mathrm{if}\,\tfrac{{r}_{\mathrm{cMpc}}}{\alpha }\leqslant 97.39\\ 0 & \mathrm{if}\,\mathrm{otherwise},\end{array}\right.\end{array}\end{eqnarray} \tag{ 17 }$$

where $\alpha \,\equiv \,{(h/0.7)}^{-1}{({{\rm{\Omega }}}_{m}/0.27)}^{-1/2}{[(1+z^{\prime} )/21]}^{-1/2}$ and r_cMpc is the co-moving distance that light has traveled from $z^{\prime}$ to z, in units of cMpc. ${f}_{\mathrm{mod}}$ is useful when calculating the inhomogeneity of the LW background if inhomogeneous source distribution is given (Ahn et al. 2009), while in this work, it only serves as the relative weight that sources at $z^{\prime}$ contributes to J_LW(z).

Because J_LW is simply the frequency-averaged intensity, actual H₂-dissociation rates by individual LW lines should be further implemented. This could be achieved by some multiplication factor weighted by line-wise dissociation rates, which would change the effective weight of source-contribution from a smooth form (${f}_{\mathrm{mod}}$) to a discrete form (f_LW in Fialkov et al. 2013), or by interpreting J_LW,th as the threshold intensity weighted by the same line-wise dissociation rates. We take the latter option in this work. We also note that in this one-zone model, Equation (16) is equivalent to the LW intensity that is averaged over the sawtooth-modulated spectrum (Haiman et al. 1997). Then, suppression of SFR by dissociation of H₂ is implemented in the form of Equations (12)–(14).

2.3.3. Ionization and Heating by X-Ray Background

Global X-ray intensity determines the heating rate and the ionization rate of the IGM outside H ii regions (or "bulk IGM" as in Mirocha 2014). X-ray photon-number intensity is given by

$$\begin{eqnarray}&&{N}_{\nu }(z)=\displaystyle \frac{{J}_{\nu }(z)}{h\nu }=\displaystyle \frac{c}{4\pi }{\left(1+z\right)}^{2}{\int }_{z}^{{z}_{f}}\displaystyle \frac{{{ \mathcal N }}_{\nu ^{\prime} }(z^{\prime} )}{H(z^{\prime} )}{{\rm{e}}}^{-{\tau }_{\nu }}{dz}^{\prime} ,\end{eqnarray} \tag{ 18 }$$

where we used Equation (15) and the fact that $\nu =\nu ^{\prime} (1+z)/(1+z^{\prime} )$. Once N_ν (z) is known, we can calculate the photoionization rate

$$\begin{eqnarray}&&{\xi }_{{\rm{X}}}(z)=4\pi {n}_{{\rm{H}}{\rm\small{I}}}(z){\int }_{{\nu }_{\min }}^{{\nu }_{\max }}{N}_{\nu }(z){\sigma }_{\nu ,{\rm{H}}{\rm\small{I}}}d\nu ,\end{eqnarray} \tag{ 19 }$$

the secondary ionization rate

$$\begin{eqnarray}&&{\tilde{\xi }}_{{\rm{X}}}(z)=4\pi {n}_{{\rm{H}}{\rm\small{I}}}(z){\int }_{{\nu }_{\min }}^{{\nu }_{\max }}{N}_{\nu }(z){\sigma }_{\nu ,{\rm{H}}{\rm\small{I}}}\displaystyle \frac{h\nu -h{\nu }_{\mathrm{th}}}{h\nu }d\nu ,\end{eqnarray} \tag{ 20 }$$

and the heating rate (erg s⁻¹)

$$\begin{eqnarray}&&{\epsilon }_{{\rm{X}}}(z)=4\pi {n}_{{\rm{H}}{\rm\small{I}}}(z){\int }_{{\nu }_{\min }}^{{\nu }_{\max }}{N}_{\nu }(z){\sigma }_{\nu ,{\rm{H}}{\rm\small{I}}}\left(h\nu -h{\nu }_{\mathrm{th}}\right)d\nu ,\end{eqnarray} \tag{ 21 }$$

where h ν_th = 13.6 eV is the hydrogen Lyman-limit energy, and σ_{ν,H I} is the photoionization cross section of the hydrogen atom at frequency ν. The ionization rate equation for the bulk IGM is then given by

$$\begin{eqnarray}&&\displaystyle \frac{{{dx}}_{{\rm{b}}}}{{dt}}=\left({\xi }_{{\rm{X}}}+{\tilde{\xi }}_{{\rm{X}}}\right)(1-{x}_{{\rm{b}}})-\alpha {{Cx}}_{{\rm{e}}}^{2}{n}_{{\rm{H}}},\end{eqnarray} \tag{ 22 }$$

where x_b is used to denote the ionized fraction of the bulk IGM and to be distinguished from the volume ionized fraction x. Even though this will affect the volume ionization rate (Equation (5)) as well, we take the approximation that 1 − x_b ≃ 1 until the reionization ends. This holds true for reionization scenarios we consider in this work. The energy rate equation is given by

$$\begin{eqnarray}&&\displaystyle \frac{3}{2}\displaystyle \frac{d}{{dt}}\left(\displaystyle \frac{{k}_{{\rm{B}}}{T}_{k}{n}_{{\rm{b}}}}{\mu }\right)={\epsilon }_{{\rm{X}}}(z)+{\epsilon }_{\mathrm{comp}}-{ \mathcal C },\end{eqnarray} \tag{ 23 }$$

where k_B is the Boltzmann constant, T_k is the kinetic temperature of gas, n_b is the proper number density of baryons, _comp is the Compton heating (when T_K < T_CMB, and cooling when T_K > T_CMB) rate, and ${ \mathcal C }$ is the cooling rate. We only include the adiabatic cooling by cosmic expansion, which is dominant over the recombination cooling and the collisional excitation+ionization cooling.

2.3.4. Lyα Background

The hydrogen Lyα background is crucial in determining the 21 cm background by decoupling the spin temperature from the CMB temperature through the Lyα pumping process, or the Wouthuysen–Field mechanism (Wouthuysen 1952; Field 1958). The photon-number intensity of the Lyα background is given by

$$\begin{eqnarray}&&{N}_{\alpha }=\displaystyle \frac{c}{4\pi }{\left(1+z\right)}^{2}\sum _{n=2}^{{n}_{\max }}{f}_{\mathrm{rec}}(n){\int }_{z}^{z{{\prime} }_{n}}{dz}^{\prime} \displaystyle \frac{{{ \mathcal N }}_{\nu ^{\prime} }(z^{\prime} )}{H(z^{\prime} )},\end{eqnarray} \tag{ 24 }$$

where ${n}_{\max }$ (=23) is the effective maximum principal quantum number of Lyman resonances, f_rec(n) is the probability for an Lyn photon (Ly1 ≡ Lyα, Ly2 ≡ Lyβ, Ly3 ≡ Lyγ, etc.) to be converted to an Lyα photon, and $z{{\prime} }_{n}$ is the redshift satisfying (Pritchard & Furlanetto 2006)

$$\begin{eqnarray}&&\displaystyle \frac{1+z{{\prime} }_{n}}{1+z}=\displaystyle \frac{1-{\left(n+1\right)}^{-2}}{1-{n}^{-2}}.\end{eqnarray} \tag{ 25 }$$

The Lyα background can also be generated by the collisional excitation of H atoms induced by energetic electrons generated by the X-ray background (e.g., Ahn et al. 2015b). However, we do not include this mechanism here because it is usually negligible when the X-ray efficiency is not extremely high. As seen in Section 2.4, we only impose a minimal level of X-ray background in this work.

The spin temperature T_s is a parameter representing the ratio of up to down states of the hyperfine structure (Field 1958):

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{1}}{{n}_{0}}=3\exp \left[-\displaystyle \frac{{T}_{* }}{{T}_{s}}\right],\end{eqnarray} \tag{ 26 }$$

where n₁ and n₀ are the number of hydrogen atoms in the up (triplet) state and the down (singlet) state, respectively, and T_* = 0.0628 K is the energy difference of the two states in terms of temperature. In the absence of Lyman resonance photons, T_s is driven to the CMB temperature T_CMB radiatively. Otherwise, the absorption and re-emission of Lyman resonance photons can drive T_s to the color temperature of Lyman lines. The dominant radiative coupling is by the Lyα photons, and the repeated scattering of Lyα photons against thermalized gas brings the Lyα color temperature T_α into T_k . The mechanical pumping by collision, separately, drives T_s into T_k . Then the spin temperature becomes

$$\begin{eqnarray}&&{T}_{s}^{-1}=\displaystyle \frac{{T}_{\mathrm{CMB}}^{-1}+({x}_{c}+{x}_{\alpha }){T}_{k}^{-1}}{1+{x}_{c}+{x}_{\alpha }}\end{eqnarray} \tag{ 27 }$$

where x_c is the collisional coupling coefficient and x_α is the Lyα pumping coefficient. x_c is given by

$$\begin{eqnarray}&&{x}_{c}=\displaystyle \frac{4\kappa (1-0){n}_{{\rm{H}}}{T}_{* }}{3{A}_{10}{T}_{\mathrm{CMB}}},\end{eqnarray} \tag{ 28 }$$

where A₁₀ = 2.85 × 10⁻¹⁵ s⁻¹ is the spontaneous emission coefficient, and κ(1 − 0) is the collisional deexcitation coefficient defined and tabulated as a function of T_k in Zygelman (2005). x_α is given by Pritchard & Furlanetto (2006)

$$\begin{eqnarray}\begin{array}{rcl}{x}_{\alpha } & = & \displaystyle \frac{16{\pi }^{2}{T}_{* }{e}^{2}{f}_{\alpha }}{27{A}_{10}{T}_{\mathrm{CMB}}{m}_{e}c}{S}_{\alpha }{N}_{\alpha }\\ & = & \displaystyle \frac{{S}_{\alpha }{N}_{\alpha }}{1.165\times {10}^{-10}\left[(1+z)/20\right]\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}\,{\mathrm{sr}}^{-1}},\end{array}\end{eqnarray} \tag{ 29 }$$

where m_e is the electron mass and S_α is a correction factor of the order of unity that accounts for the distortion of the line profile by thermalized atoms and peculiar motion (Chen & Miralda-Escudé 2004; Chuzhoy & Shapiro 2006; Hirata 2006). In this paper, we adopt the functional form of S_α suitable for co-moving gas without peculiar motion (Chuzhoy & Shapiro 2006):

$$\begin{eqnarray}&&{S}_{\alpha }=\displaystyle \frac{\exp \left[-0.37{\left(1+z\right)}^{1/2}{T}_{k}^{-2/3}\right]}{1+0.4{T}_{k}^{-1}}.\end{eqnarray} \tag{ 30 }$$

We cover a limited but representative set of parameters in calculating reionization histories. For the vanilla model, we use parameters that are sampled similarly to Furlanetto (2006) and Bernardi et al. (2015). For SRI, we use parameters including those of Iliev et al. (2007), which first suggested the self-regulation scheme of the model. For SRII, we use parameters including those of Ahn et al. (2012), which provide the physical basis of the model. Parameters and some characteristics of reionization are listed in Tables 1–3. For the vanilla model and SRI, we accommodate both dF and F star formation scenarios (Section 2.1). The resulting x(z)'s are plotted in Figures 2–5, and overlaid on the 68% and 95% constraints from the Planck Legacy Data ("PLD"; Planck Collaboration et al. 2020). Note that in these figures we show the ionized volume fraction in terms of the electron fraction ${x}_{{\rm{e}}}\equiv \left\langle {n}_{{\rm{e}}}\right\rangle /\left\langle {n}_{{\rm{H}}}\right\rangle$, with $\left\langle \,\,\right\rangle$ being the volume average, such that ${x}_{{\rm{e}}}=(1+\left\langle {n}_{\mathrm{He}}\right\rangle /\left\langle {n}_{{\rm{H}}}\right\rangle )x\,=\,1.079\,x$ if helium atoms are assumed singly ionized in H ii regions. The PLD constraints are in fact on x_e, and they are shown as shaded regions in Figures 2–5.

Table 1. Model Parameters and Resulting Characteristics of the Vanilla and the SRI

Model	fγ	gγ	${f}_{\gamma }^{{\rm{H}}}$	${f}_{\gamma }^{{\rm{L}}}$	${g}_{\gamma }^{{\rm{H}}}$	${g}_{\gamma }^{{\rm{L}}}$	zbegin	zend	Δz3−97	τes	τes(15–30)	χ2/ν
V-L_dF (P)	19.62	⋯	⋯	⋯	⋯	⋯	15.87	5.50	8.18	0.05998	3.97E-4	0.81
V-M1_dF	23.55	⋯	⋯	⋯	⋯	⋯	16.32	6.21	7.98	0.06638	5.11E-4	0.86
V-M2_dF	34.88	⋯	⋯	⋯	⋯	⋯	17.01	7.27	7.66	0.07647	7.57E-4	1.19
V-H_dF	56.69	⋯	⋯	⋯	⋯	⋯	17.85	8.52	7.33	0.08899	1.23E-3	2.31
V-L_F	⋯	0.4077	⋯	⋯	⋯	⋯	12.21	5.49	5.06	0.04794	2.73E-5	0.86
V-M1_F (P)	⋯	0.7042	⋯	⋯	⋯	⋯	13.00	6.25	5.10	0.05555	4.71E-5	0.81
V-M2_F	⋯	1.472	⋯	⋯	⋯	⋯	14.06	7.32	5.11	0.06651	9.85E-5	0.87
V-H_F	⋯	3.140	⋯	⋯	⋯	⋯	15.13	8.43	5.10	0.07852	2.10E-4	1.32
SRI-L0_dF (P)	⋯	⋯	26.02	0	⋯	⋯	12.82	5.81	5.56	0.05436	2.01E-5	0.82
SRI-LL_dF	⋯	⋯	26.02	91.93	⋯	⋯	16.04	6.05	7.96	0.06380	4.14E-4	0.83
SRI-LH_dF	⋯	⋯	26.02	919.3	⋯	⋯	19.11	6.82	10.41	0.09311	2.77E-3	2.93
SRI-HL_dF	⋯	⋯	91.93	91.93	⋯	⋯	16.22	8.25	6.15	0.08057	4.60E-4	1.46
SRI-HH_dF	⋯	⋯	91.93	919.3	⋯	⋯	19.12	8.78	8.54	0.10083	2.81E-3	4.54
SRI-0H_dF	⋯	⋯	0	919.3	⋯	⋯	19.11	⋯	⋯	0.07659	2.75E-3	1.80
SRI-L0_F	⋯	⋯	⋯	⋯	0.8673	0	10.81	5.78	3.84	0.04756	1.53E-6	0.87
SRI-LL_F (P)	⋯	⋯	⋯	⋯	0.8673	8.673	14.41	6.24	6.49	0.06140	1.16E-4	0.81
SRI-LH_F	⋯	⋯	⋯	⋯	0.8673	86.73	17.09	6.93	8.34	0.08816	8.63E-4	2.22
SRI-HL_F	⋯	⋯	⋯	⋯	6.938	8.673	14.59	8.27	4.82	0.07562	1.25E-4	1.16
SRI-HH_F	⋯	⋯	⋯	⋯	6.938	86.73	17.10	8.79	6.74	0.09310	8.71E-4	2.95
SRI-0H_F	⋯	⋯	⋯	⋯	0	86.73	17.09	4.10	10.12	0.08597	8.62E-4	1.98

Note. SRI-0H_dF does not reach the end of reionization, not even x = 0.97, until z = 0. Denoted by "(P)" are models that best fit the observed ${C}_{l}^{\mathrm{EE}}$ of the PLD.

Download table as:  ASCIITypeset image

All models are described by simple ordinary differential equations (ODEs), and thus can be easily integrated with ODE solvers. We start numerical integration from z = 40, when the contribution of any type of halo to reionization and heating is believed to be negligible. The initial value of x is set to an arbitrarily small value, because the volume occupied by H ii regions at z = 40 must be negligible. Ionization rate equations to solve are not stiff, and we use a fourth-order, adaptive Runge–Kutta integrator with both the relative tolerance and the absolute tolerance of x set to 10⁻⁵. For SRII, we need an extra effort to calculate J_LW(t) at any time t, because J_LW(t) regulates SFR inside MHs at t via x_LW and impacts PPR (Equations (12) and (13)). Therefore, we calculate x(t) and J(t) at each incrementally increasing time step, by integrating Equation (5) with Equation (13) (or Equation (12) if dF assumed) and using Equations (14) and (16).

The vanilla model shows a smooth and monotonic evolution of x (Figure 1). The monotonic behavior of x(t) is easily explained by the fact that ξ is proportional to the monotonically increasing f_coll (F scenario) or df_coll/dt (dF scenario). This characteristic makes τ_es and z_end tightly correlated, and thus lacks the "leverage" to accommodate reionization scenarios that are different in x(z) but degenerate in τ_es and z_end, as long as one type of star formation scenario is chosen from F or dF. One can of course make a somewhat more sophisticated variant of this model by, e.g., allowing multiple species of halos with different f_γ 's (a mixture of Population II and Population III stars; Furlanetto 2006). Nevertheless, due to the lack of any self-regulation, the resulting reionization histories of such variants would still remain similar to the original vanilla model. The duration of reionization is in general more extended in the dF scenario than the F scenario. This is due to the fact that df_coll/dt grows more slowly than f_coll. Therefore, for given z_end, the dF scenario produces larger τ_es than the F scenario. This tendency is clearly presented in Figure 1 and Table 1.

Zoom In Zoom Out Reset image size

The SRI model adds a little more complexity to the characteristics compared to the vanilla model (Figure 2). Due to the existence of self-regulation, the SRI is expected to have more extended reionization histories than the vanilla model. This has indeed been shown to be the case for a consistent halo-selection criterion for HMACHs and LMACHs (Iliev et al. 2007). However, one subtlety in our modeling scheme complicates such an expectation. Because we use a constant-mass criterion in SRI, while a constant-temperature criterion in the vanilla model, we find that in some cases, the duration of SRI models can be shorter than that of vanilla models. Had we used the same halo-selection criterion, SRI models would have larger Δz than vanilla models, which we actually tested and confirmed. Aside from this complication, which is not essential, the general trend is that (1) the larger the value of ${g}_{\gamma }^{{\rm{L}}}/{g}_{\gamma }^{{\rm{H}}}$ is, the larger the duration of reionization becomes, and (2) the addition of LMACHs to HMACH-only scenarios extends the duration of reionization. Also, cases with very large ${g}_{\gamma }^{{\rm{L}}}/{g}_{\gamma }^{{\rm{H}}}$ (e.g., the SRI-0H_F case in Figure 2) slow down reionization significantly at the end of reionization, producing histories as symmetric as the tangent-hyperbolic model that has been used extensively in the analysis of the CMB data. It is easy to understand this behavior: ${g}_{\gamma }^{{\rm{L}}}/{g}_{\gamma }^{{\rm{H}}}$ is a rough measure of the relative contribution of LMACHs to reionization to that of HMACHs, and the self-regulation becomes stronger as x becomes larger. One very extreme case is SRI-0H_dF, which never finishes reionization due to a strong self-regulation. The trend that Δz is larger in the dF scenario than in the F scenario is the same as in the vanilla model. The general trend of the vanilla and SRI models can also be seen in Figure 3 in terms of ξ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The SRII model has features richer than the vanilla and the SRI models (Figures 4 and 5). The most notable feature is the existence of the early, extended, and slowly increasing phase in x. This is due to the self-regulation of star formation, even stronger than that in the SRI model, which takes place inside MHs. The star formation inside MHs is mainly regulated by the LW background J_LW, which quickly builds up to reach J_LW,th. Continuum photons below the Lyman limit and emitted at redshift z travel a cosmological distance (∼100{(1 + z)/21}^−0.5 cMpc), in contrast to the hydrogen-ionizing UV photons that travel up to the ionization front and are then absorbed. Therefore, any newly forming MHs will be under the influence of the LW background long before being exposed to the ionizing photons, and any pre-ionized region would have been under the overcritical LW intensity (J_LW > J_LW,th). We find that this is indeed the case: when tested with $\max ({x}_{\mathrm{LW}},x)$ replaced by x_LW in Equation (13), the resulting x(z) was not affected.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Another notable feature of SRII is that, in some cases where the contribution of ionizing photons by MHs is as significant as to drive x beyond ≳10%, there exists a phase where x decreases in time. ⁴ This is mainly due to the fact that (1) the LW feedback renders J_LW ≃ J_{LW th} when MHs dominate as the main radiation sources (Figure 11), and (2) the large difference in the number of soft-UV (h ν = ∼ 11 − 13.6 eV) photons per ionizing photon of Population II and Population III stars (see the detail in Section 2.4), and (3) the drop of f_esc from MH values (f_esc ∼ 1) to ACHs (f_esc ∼ 0.2). Then, LMACHs (assumed to host Population III stars) and HMACHs (assumed to host Population II stars) start to generate soft-UV photons to make up J_LW near J_{LW th}. This is achieved at the expense of ACHs putting out a much smaller amount of ionizing photons than MHs (assumed to host Population III stars). Consequently, the IGM gains a chance to recombine faster than ionization (Figure 6: see dips in ξ). In practice, however, this recombination is slight and not as dramatic as the "double reionization" that has been suggested by Cen (2003). Given the same set of $\left\{{g}_{\gamma }^{{\rm{H}}},{g}_{\gamma }^{{\rm{L}}}\right\}$ as in SRI, the SRII model has z_end that is practically identical to that of SRI while it has τ_es and Δz that are both boosted from those of SRI (Tables 1–3). This is simply due to the existence of additional photon sources, or MH stars, that ionize the IGM only to a limited extent (x ≲ 15% at most under our parameter range, but may be increased if M_III and J_LW,th are pushed to higher values) such that z_end is not much affected but can increase τ_es and extend the duration of reionization substantially by strongly regulated ionization history.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The 2018 PLD provides, among improvements from the 2015 data, so far the most precise measurement of the large-scale CMB polarization anisotropy. The large-scale E-mode autocorrelation angular power spectrum, ${C}_{l}^{\mathrm{EE}}$, is strongly affected by the history of reionization. The quadrupole moment of the CMB anisotropy generates linear polarization after Thomson scattering from the viewpoint of an electron, and the polarization signal is observed after being modulated by the relevant wave-modes, resulting in affecting ${C}_{l}^{\mathrm{EE}}$ mostly in the low-l (≲20) regime (Hu & White 1997; Haiman & Knox 1999; Dodelson 2003).

The extended high-redshift (z ≳ 15) ionization tail predicted by Ahn et al. (2012) and SRII models here has been advocated by Miranda et al. (2017) and Heinrich & Hu (2018), based on their principal component analysis (PCA) of the Planck 2015 ${C}_{l}^{\mathrm{EE}}$ data observed through the LFI. They claimed that a specific SRII model with a substantial high-redshift tail, corresponding roughly to the L2M1J2 case of Ahn et al. (2012), was favored over the vanilla model at ∼ 1σ level. Later, Millea & Bouchet (2018) used both the LFI and the HFI (proprietary at the time) data, with a well-handled physicality (x > 0) prior, to claim against too much contribution from the z ≳ 15 epoch. They constrained the optical depth from 15 ≤ z ≤ 30, or τ(15, 30), to τ(15, 30) < 0.015 at the 2σ level. The Planck 2018 analysis based only on the low-l E-mode polarization further reduces this value to τ(15, 30) < 0.007 at the 2σ level and ${\tau }_{\mathrm{es}}={0.0504}_{-0.0079}^{+0.0050}$ at the 1σ level (Planck Collaboration et al. 2020).

In light of the constraints described above, we compare our model x(z)'s from Section 3.1 to PLD. The main purpose of this task is to (1) understand whether any class of our models are preferred by observation and (2) whether the degeneracy of models in τ_es and z_ov can be broken. For example, Ahn et al. (2012) showed that SRII models with an extended tail in x(z) could be distinguished from the vanilla- or SRI-type models, even when the models have the same τ_es (=0.085) and z_ov (=6.8). The pictorial comparison of x(z) to the Planck constraint (1σ and 2σ constraints shown in the shaded regions) is shown in Figures 1–5. We test the relative goodness of several selected models by calculating the reduced chi square,

$$\begin{eqnarray}&&{\chi }^{2}/\nu =\displaystyle \frac{1}{29}\sum _{l=2}^{30}\displaystyle \frac{{\left({D}_{l}^{\mathrm{EE}}-{\tilde{D}}_{l}^{\mathrm{EE}}\right)}^{2}}{{\sigma }_{l}^{2}},\end{eqnarray} \tag{ 31 }$$

where ${D}_{l}^{\mathrm{EE}}\equiv l(l+1){C}_{l}^{\mathrm{EE}}/(2\pi )$ and ${\tilde{D}}_{l}^{\mathrm{EE}}\equiv l(l+1){\tilde{C}}_{l}^{\mathrm{EE}}/(2\pi )$ are the E-mode power spectra corresponding to a model x(z) and PLD, respectively, and σ_l is the standard deviation of ${\tilde{D}}_{l}^{\mathrm{EE}}$ due to the cosmic variance and the noise. ⁵ ${C}_{l}^{\mathrm{EE}}\equiv \left\langle {\left|{a}_{{lm}}^{{\rm{E}}}\right|}^{2}\right\rangle$ for given l averaged over m = [−l, l], with the spherical-harmonics decomposition of the E-mode anisotropy $E(\theta ,\phi )={\sum }_{{lm}}{a}_{{lm}}^{{\rm{E}}}{Y}_{{lm}}(\theta ,\phi )$. In calculating ${D}_{l}^{\mathrm{EE}}$, we use a version of the Boltzmann solver CAMB that was modified to allow for a generic shape of x(z) (Mortonson & Hu 2008; downloadable from http://background.uchicago.edu/camb_rpc/). For the base cosmology, we use the best-fit parameter set of PLD. While this is not a full likelihood analysis including other data products such as the temperature anisotropy, the value of χ²/ν from Equation (31) can indicate the relative goodness of models because the impact of reionization histories is the strongest in the E-mode (see, e.g., an identical approach by Qin et al. 2020b). E-mode power spectra of selected models against the PLD are plotted in Figure 8.

Zoom In Zoom Out Reset image size

It is interesting to note that the constraint on reionization by the Planck observation provides a good match to the observed Γ (Section 2.3.1). If the observed Γ (e.g., Calverley et al. 2011) is translated into ${\dot{N}}_{\mathrm{ion}}$ with a reasonable set of physical parameters (λ_mfp = 10 cMpc: the IGM mean free path to H-ionizing photons at z ≃ 6; and α_s = α_b = 2: the spectral hardness of H-ionizing photons in Equation (21) of Bolton & Haehnelt 2007), PLD-favored models have a good agreement with the observed Γ at z ≃ 6 (Figures 3 and 6). This consistency between the two independent observations, even though uncertainties are large, is encouraging.

We now claim that some SRII models with substantial high-redshift tails are still among those highly favored by PLD. Because χ²/ν is a measure of the goodness of fit, we can use the value of χ²/ν to find the PLD-favored models. We indeed find many models can explain the PLD low-l ${C}_{l}^{\mathrm{EE}}$ fairly well even though the variance in τ_es of such models is substantial. Models that fit the PLD ${C}_{l\leqslant 30}^{\mathrm{EE}}$ best are marked by arrows in Figure 7 and denoted by "(P)" in model names in Tables 1–3: all of these models have almost the same likelihood with χ²/ν = 0.80–0.82, but with a substantial spread on τ_es with τ_es ≃ [0.0544, 0.0643]. If allowance is extended to models with χ²/ν ≤ 1, then the allowed optical depth becomes τ_es ≃ [0.044, 0.072]. This indicates that some of our SRII models are still well within the PLD constraint and are favored as much as those models without high-redshift ionization tails. As seen in Tables 1–3, the most favored models with the smallest values of χ²/ν, χ²/ν = 0.80, are indeed SRII-L0-100-e1.0-J0.1 and SRII-L0-300-e1.0-J0.1, which have a substantial high-redshift tail that reaches maximum x = 0.12 at z = 13.6 and x = 0.15 at z = 13.6, respectively. Such tails contribute to τ_es(15–30) substantially: τ_es(15–30) = 0.0063 and 0.0078 for SRII-L0-100-e1.0-J0.1 and SRII-L0-300-e1.0-J0.1, respectively. We also note that these maximum-likelihood models are clustered around τ_es = 0.06, substantially different from the inferred value τ_es = 0.0504 by PLD-${C}_{l\leqslant 30}^{\mathrm{EE}}$. The reason why such a difference occurs is unclear; this is nevertheless a very important issue and a further investigation is warranted.

The model with the strongest ionization tail of all, SRII-L0-300-e1.0-J0.1, is worthy of close attention. Compared to the 2σ constraint of PLD, τ_es(15–30) < 0.007, SRII-L0-300-e1.0-J0.1 actually violates this constraint with τ_es(15–30) = 0.0078 but is still the best fit (to PLD ${C}_{l\leqslant 30}^{\mathrm{EE}}$) of all of the models we tested. This model also has τ_es = 0.0596, which is about 2σ away from the E-mode-only best-fit estimate by PLD, τ_es = 0.0504. If we do not consider other CMB observables and assume a flat prior, we can conclude that this model is as good as or just slightly better than other tailless models with χ²/ν = 0.81–0.82. It is interesting to see that there exists a weak tension between τ_es's estimated by the low-l E-mode polarization and the CMB lensing: the low-l E-mode data of PLD prefers such a low ${\tau }_{\mathrm{es}}(={0.0504}_{-0.0079}^{+0.0050})$, while the CMB lensing of PLD prefers higher τ_es at around τ_es ≃ 0.08 (Planck Collaboration et al. 2020). This tendency of the CMB lensing favoring large values of τ_es, even though the uncertainty is large, is on par with favoring two-stage reionization models with substantial ionization tails. Our findings are in slight disagreement with the PLD constraint that was constructed using nonparametric Bayesian inference. Based on our forward modeling and goodness-of-fit approach, we argue that a family of models with a substantial high-redshift ionization tail reaching ${x}_{\max }\sim 0.12$ is still just as strong a contender at the moment as the tailless models.

Will there be a chance to probe a high-redshift ionization tail in the future? The high-redshift tail tends to boost ${C}_{l}^{\mathrm{EE}}$ at 14 ≲ l ≲ 24 (e.g., Ahn et al. 2012; Miranda et al. 2017). SRII models in Figure 8 produce ${C}_{14\leqslant l\leqslant 24}^{\mathrm{EE}}$ larger than that of the rest of models, and in particular, SRII-L0-300-e1.0-J0.1 has the strongest ${C}_{l}^{\mathrm{EE}}$ at 8 ≤ l ≤ 24. In principle, models degenerate in τ_es can have different ${C}_{l\lt 30}^{\mathrm{EE}}$'s due to the variance in x(z). Ahn et al. (2012), using a PCA, had indeed predicted that high-precision CMB observations could break the degeneracy in τ_es and probe (or disprove) the existence of the high-redshift ionization tail. We stress this point again through Figure 9, showing two models from Ahn et al. (2012) that are degenerate both in τ_es and z_end but are clearly different in x(z), especially in the existence of the high-z ionization tail, and in the resulting ${C}_{l}^{\mathrm{EE}}$. From Figures 8 and 9, we observe that the boost of ${C}_{14\leqslant l\leqslant 24}^{\mathrm{EE}}$ in the two-stage reionization models with the high-z ionization tail against the tailless models is a universal effect. As seen in Figure 7, the relation between τ_es and χ²/ν is not exactly monotonic but instead there exists some scatter in χ²/ν for the same τ_es and vice versa. Such a scatter increases as τ_es increases, which is due to the increased freedom in constructing x(z) for given τ_es. However, because the PLD E-mode power spectrum prefers such a low τ_es, as of now the leverage of having a pronounced tail has somewhat diminished from that prediction. Nevertheless, it is possible that observation by a more accurate apparatus might find a preference for a higher τ_es than Planck that is still hampered by the large noise in measuring the polarization anisotropy. Therefore, we need a better apparatus than Planck to (1) see whether or not τ_es could get larger than the estimate by PLD to allow for more pronounced two-stage reionization models and (2) break the model degeneracy in τ_es better than Planck to probe the ionization tail even when the tail is weak.

Zoom In Zoom Out Reset image size

We also briefly describe another type of constraint from CMB observations. The kinetic Sunyaev–Zel'dovich effect can arise from the peculiar motion of H ii bubbles during EoR and can affect the small-scale (l ∼ a few thousand) temperature anisotropy power spectrum ${C}_{l}^{\mathrm{TT}}$. Measurement of ${C}_{l}^{\mathrm{TT}}$ by the South Pole Telescope, especially ${C}_{l=3000}^{\mathrm{TT}}$ (Reichardt et al. 2012), was used by Zahn et al. (2012) to constrain the duration of reionization Δz ≡ z(x = 0.25) − z(x = 0.99) to Δz < 4 − 7 at the 2σ level (depending on the assumed correlation between the thermal Sunyaev–Zel'dovich effect and the cosmic infrared background; see also the similar assessment by Mesinger et al. 2012 and Battaglia et al. 2013). Without MH stars, reionization occurs always in a patchy way and thus any addition of electrons, or equivalently extension of Δz, increases ${C}_{l=3000}^{\mathrm{TT}}$ monotonically, as was assumed in Zahn et al. (2012), Mesinger et al. (2012), and Battaglia et al. (2013). However, Park et al. (2013) readdressed this issue with a variety of reionization scenarios including the SRII-type, and found that the added duration of reionization beyond this limit could still be accommodated by the measured ${C}_{l=3000}^{\mathrm{TT}}$. As claimed in Park et al. (2013), H ii regions by MHs are distributed almost uniformly (Ahn et al. 2012), and thus the increase in Δz in SRII models does not guarantee an increase in ${C}_{l=3000}^{\mathrm{TT}}$. Therefore, the largeness of Δz₃₋₉₇ ≡ z(0.03) − z(x = 0.97) of many SRII models (see, e.g., those denoted by "(P)" in Tables 2 and 3) should not be considered as a violation of such a constraint. Instead, constraining Δz using the small-scale ${C}_{l}^{\mathrm{TT}}$ should be restricted to only a limited set of models without MHs.

Table 2. Model Parameters and Resulting Characteristics of SRII, with M_III = 100 M_⊙

Model	${g}_{\gamma }^{{\rm{H}}}$	${g}_{\gamma }^{{\rm{L}}}$	MIII	${f}_{\mathrm{esc}}^{{\rm{M}}}$	JLW,th	zbegin	zend	Δz3−97	τes	τes(15–30)	χ2/ν
SRII-L0-100-e0.1-J0.05	0.8673	0	100	0.1	0.05	10.86	5.78	3.86	0.04801	3.45E-4	0.86
SRII-L0-100-e0.1-J0.10	0.8673	0	100	0.1	0.10	11.03	5.78	3.88	0.04850	6.18E-4	0.86
SRII-L0-100-e0.5-J0.05	0.8673	0	100	0.5	0.05	11.18	5.78	3.90	0.04983	1.72E-3	0.85
SRII-L0-100-e0.5-J0.10	0.8673	0	100	0.5	0.10	22.15	5.79	4.07	0.05226	3.08E-3	0.83
SRII-L0-100-e1.0-J0.05	0.8673	0	100	1.0	0.05	22.91	5.79	3.96	0.05210	3.44E-3	0.83
SRII-L0-100-e1.0-J0.10 (P)	0.8673	0	100	1.0	0.10	25.38	5.79	4.51	0.05696	6.17E-3	0.80
SRII-LL-100-e0.1-J0.05 (P)	0.8673	8.673	100	0.1	0.05	14.48	6.24	6.50	0.06168	3.75E-4	0.81
SRII-LL-100-e0.1-J0.10 (P)	0.8673	8.673	100	0.1	0.10	14.90	6.24	6.53	0.06199	6.37E-4	0.81
SRII-LL-100-e0.5-J0.05 (P)	0.8673	8.673	100	0.5	0.05	15.26	6.24	6.54	0.06284	1.44E-3	0.82
SRII-LL-100-e0.5-J0.10	0.8673	8.673	100	0.5	0.10	22.14	6.24	6.76	0.06446	2.78E-3	0.83
SRII-LL-100-e1.0-J0.05 (P)	0.8673	8.673	100	1.0	0.05	22.90	6.24	6.59	0.06432	2.79E-3	0.82
SRII-LL-100-e1.0-J0.10	0.8673	8.673	100	1.0	0.10	25.38	6.24	13.93	0.06762	5.50E-3	0.86
SRII-LH-100-e0.1-J0.05	0.8673	86.73	100	0.1	0.05	17.11	6.93	8.35	0.08832	1.01E-3	2.23
SRII-LH-100-e0.1-J0.10	0.8673	86.73	100	0.1	0.10	17.19	6.93	8.35	0.08847	1.15E-3	2.24
SRII-LH-100-e0.5-J0.05	0.8673	86.73	100	0.5	0.05	17.22	6.93	8.36	0.08895	1.60E-3	2.29
SRII-LH-100-e0.5-J0.10	0.8673	86.73	100	0.5	0.10	22.05	6.93	8.41	0.08973	2.37E-3	2.37
SRII-LH-100-e1.0-J0.05	0.8673	86.73	100	1.0	0.05	17.46	6.93	8.37	0.08976	2.35E-3	2.36
SRII-LH-100-e1.0-J0.10	0.8673	86.73	100	1.0	0.10	25.38	6.93	8.50	0.09135	3.91E-3	2.55
SRII-HL-100-e0.1-J0.05	6.938	8.673	100	0.1	0.05	14.62	8.27	4.82	0.07586	3.48E-4	1.17
SRII-HL-100-e0.1-J0.10	6.938	8.673	100	0.1	0.10	14.71	8.27	4.83	0.07610	5.70E-4	1.18
SRII-HL-100-e0.5-J0.05	6.938	8.673	100	0.5	0.05	14.79	8.27	4.84	0.07684	1.26E-3	1.20
SRII-HL-100-e0.5-J0.10	6.938	8.673	100	0.5	0.10	22.13	8.27	4.91	0.07807	2.40E-3	1.26
SRII-HL-100-e1.0-J0.05	6.938	8.673	100	1.0	0.05	15.35	8.27	4.86	0.07807	2.41E-3	1.25
SRII-HL-100-e1.0-J0.10	6.938	8.673	100	1.0	0.10	25.38	8.27	5.04	0.08057	4.70E-3	1.38
SRII-HH-100-e0.1-J0.05	6.938	86.73	100	0.1	0.05	17.12	8.79	6.74	0.09325	1.01E-3	2.96
SRII-HH-100-e0.1-J0.10	6.938	86.73	100	0.1	0.10	17.19	8.79	6.75	0.09339	1.16E-3	2.98
SRII-HH-100-e0.5-J0.05	6.938	86.73	100	0.5	0.05	17.22	8.79	6.75	0.09387	1.60E-3	3.04
SRII-HH-100-e0.5-J0.10	6.938	86.73	100	0.5	0.10	22.05	8.79	6.79	0.09464	2.35E-3	3.15
SRII-HH-100-e1.0-J0.05	6.938	86.73	100	1.0	0.05	17.43	8.79	6.76	0.09467	2.34E-3	3.14
SRII-HH-100-e1.0-J0.10	6.938	86.73	100	1.0	0.10	25.38	8.79	6.88	0.09622	3.87E-3	3.38
SRII-0H-100-e0.1-J0.05	0	86.73	100	0.1	0.05	17.11	4.10	10.12	0.08612	1.01E-3	1.99
SRII-0H-100-e0.1-J0.10	0	86.73	100	0.1	0.10	17.19	4.10	10.13	0.08627	1.15E-3	2.00
SRII-0H-100-e0.5-J0.05	0	86.73	100	0.5	0.05	17.22	4.10	10.13	0.08676	1.60E-3	2.04
SRII-0H-100-e0.5-J0.10	0	86.73	100	0.5	0.10	22.05	4.10	10.18	0.08754	2.37E-3	2.12
SRII-0H-100-e1.0-J0.05	0	86.73	100	1.0	0.05	17.46	4.10	10.15	0.08756	2.35E-3	2.11
SRII-0H-100-e1.0-J0.10	0	86.73	100	1.0	0.10	25.38	4.10	10.27	0.08916	3.92E-3	2.27

Note. In the table, M_III and J_LW,th are in units of M_⊙ and 10⁻²¹ erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹, respectively. Denoted by "(P)" are models that fit the observed ${C}_{l}^{\mathrm{EE}}$ of the PLD best.

Download table as:  ASCIITypeset image

Table 3. Model Parameters and Resulting Characteristics of SRII, with M_III = 300 M_⊙

Model	${g}_{\gamma }^{{\rm{H}}}$	${g}_{\gamma }^{{\rm{L}}}$	MIII	${f}_{\mathrm{esc}}^{{\rm{M}}}$	JLW,th	zbegin	zend	Δz3−97	τes	τes(15–30)	χ2/ν
SRII-L0-300-e0.1-J0.05	0.8673	0	300	0.1	0.05	10.89	5.78	3.86	0.04817	4.50E-4	0.86
SRII-L0-300-e0.1-J0.10	0.8673	0	300	0.1	0.10	11.10	5.78	3.89	0.04875	7.75E-4	0.86
SRII-L0-300-e0.5-J0.05	0.8673	0	300	0.5	0.05	11.45	5.78	3.92	0.05063	2.29E-3	0.84
SRII-L0-300-e0.5-J0.10 (P)	0.8673	0	300	0.5	0.10	23.26	5.79	4.14	0.05350	3.88E-3	0.82
SRII-L0-300-e1.0-J0.05 (P)	0.8673	0	300	1.0	0.05	23.50	5.79	4.02	0.05365	4.50E-3	0.82
SRII-L0-300-e1.0-J0.10 (P)	0.8673	0	300	1.0	0.10	27.15	5.79	15.28	0.05960	7.79E-3	0.80
SRII-LL-300-e0.1-J0.05 (P)	0.8673	8.673	300	0.1	0.05	14.51	6.24	6.50	0.06177	4.49E-4	0.81
SRII-LL-300-e0.1-J0.10 (P)	0.8673	8.673	300	0.1	0.10	15.26	6.24	6.54	0.06215	7.73E-4	0.82
SRII-LL-300-e0.5-J0.05 (P)	0.8673	8.673	300	0.5	0.05	20.00	6.24	6.55	0.06330	1.83E-3	0.82
SRII-LL-300-e0.5-J0.10	0.8673	8.673	300	0.5	0.10	23.26	6.24	6.86	0.06528	3.45E-3	0.83
SRII-LL-300-e1.0-J0.05	0.8673	8.673	300	1.0	0.05	23.50	6.24	6.65	0.06527	3.57E-3	0.83
SRII-LL-300-e1.0-J0.10	0.8673	8.673	300	1.0	0.10	27.15	6.24	14.77	0.06944	6.96E-3	0.88
SRII-LH-300-e0.1-J0.05	0.8673	86.73	300	0.1	0.05	17.12	6.93	8.35	0.08836	1.04E-3	2.23
SRII-LH-300-e0.1-J0.10	0.8673	86.73	300	0.1	0.10	17.21	6.93	8.36	0.08856	1.23E-3	2.25
SRII-LH-300-e0.5-J0.05	0.8673	86.73	300	0.5	0.05	17.26	6.93	8.36	0.08917	1.77E-3	2.30
SRII-LH-300-e0.5-J0.10	0.8673	86.73	300	0.5	0.10	23.22	6.93	8.42	0.09022	2.78E-3	2.41
SRII-LH-300-e1.0-J0.05	0.8673	86.73	300	1.0	0.05	17.69	6.93	8.38	0.09024	2.74E-3	2.40
SRII-LH-300-e1.0-J0.10	0.8673	86.73	300	1.0	0.10	27.15	6.93	8.54	0.09235	4.76E-3	2.64
SRII-HL-300-e0.1-J0.05	6.938	8.673	300	0.1	0.05	14.63	8.27	4.82	0.07594	4.10E-4	1.17
SRII-HL-300-e0.1-J0.10	6.938	8.673	300	0.1	0.10	14.75	8.27	4.83	0.07623	6.85E-4	1.18
SRII-HL-300-e0.5-J0.05	6.938	8.673	300	0.5	0.05	14.91	8.27	4.85	0.07722	1.58E-3	1.22
SRII-HL-300-e0.5-J0.10	6.938	8.673	300	0.5	0.10	23.26	8.27	4.93	0.07875	2.97E-3	1.28
SRII-HL-300-e1.0-J0.05	6.938	8.673	300	1.0	0.05	23.50	8.27	4.88	0.07880	3.01E-3	1.28
SRII-HL-300-e1.0-J0.10	6.938	8.673	300	1.0	0.10	27.15	8.27	5.12	0.08198	5.90E-3	1.45
SRII-HH-300-e0.1-J0.05	6.938	86.73	300	0.1	0.05	17.12	8.79	6.74	0.09329	1.05E-3	2.97
SRII-HH-300-e0.1-J0.10	6.938	86.73	300	0.1	0.10	17.20	8.79	6.75	0.09349	1.24E-3	2.99
SRII-HH-300-e0.5-J0.05	6.938	86.73	300	0.5	0.05	17.27	8.79	6.75	0.09411	1.78E-3	3.06
SRII-HH-300-e0.5-J0.10	6.938	86.73	300	0.5	0.10	23.22	8.79	6.81	0.09513	2.77E-3	3.21
SRII-HH-300-e1.0-J0.05	6.938	86.73	300	1.0	0.05	17.58	8.79	6.77	0.09512	2.69E-3	3.19
SRII-HH-300-e1.0-J0.10	6.938	86.73	300	1.0	0.10	27.15	8.79	6.91	0.09721	4.70E-3	3.50
SRII-0H-300-e0.1-J0.05	0	86.73	300	0.1	0.05	17.12	4.10	10.12	0.08617	1.04E-3	1.99
SRII-0H-300-e0.1-J0.10	0	86.73	300	0.1	0.10	17.21	4.10	10.13	0.08637	1.24E-3	2.01
SRII-0H-300-e0.5-J0.05	0	86.73	300	0.5	0.05	17.27	4.10	10.13	0.08698	1.78E-3	2.06
SRII-0H-300-e0.5-J0.10	0	86.73	300	0.5	0.10	23.22	4.10	10.19	0.08803	2.79E-3	2.16
SRII-0H-300-e1.0-J0.05	0	86.73	300	1.0	0.05	17.69	4.10	10.15	0.08802	2.71E-3	2.14
SRII-0H-300-e1.0-J0.10	0	86.73	300	1.0	0.10	27.15	4.10	10.31	0.09016	4.77E-3	2.36

Note. The unit convention is the same as that in Table 2. Denoted by "(P)" are models that fit the observed ${C}_{l}^{\mathrm{EE}}$ of the PLD best.

Download table as:  ASCIITypeset image

The main variants determining δ T_b are the X-ray heating efficiency and the Lyα intensity, which determine T_K and x_α , respectively. The X-ray efficiency is not a direct product of the stellar radiation and is thus the main cause of the uncertainty in δ T_b . The Lyα intensity, on the other hand, is almost solely determined by the stellar radiation and is closely related to the ionizing PPR and the SED. We do not consider the creation of Lyα photons due to the excitation of H atoms by the X-ray-induced electrons, which is a good approximation unless the X-ray efficiency is extremely high (f_x ≫ 1 with f_x in Equation (32)). For the X-ray efficiency, we use the common parameter f_X (Furlanetto 2006; Mirocha 2014), defined as the fudge parameter connecting the co-moving X-ray luminosity density ${{ \mathcal L }}_{\nu }$ ($=h\nu {{ \mathcal N }}_{\nu };$ in erg s⁻¹ Hz⁻¹ cMpc⁻³) to SFRD (Section 2.1):

$$\begin{eqnarray}&&{f}_{X}=\displaystyle \frac{\int d\nu {{ \mathcal L }}_{\nu }}{{c}_{X}{\rm{SFRD}}},\end{eqnarray} \tag{ 32 }$$

where the additional proportionality coefficient c_X is fixed to c_X = 3.4 × 10⁴⁰ erg s⁻¹(M_⊙ yr⁻¹)⁻¹, an extrapolation of the 2–10 keV relation between ${{ \mathcal L }}_{\nu }$ and SFRD (or equivalently between L_ν and SFR on average galaxies) by Grimm et al. (2003) to h ν ≥ 0.2 keV. Here, we limit ${{ \mathcal L }}_{\nu }$ to the energy range h ν = [0.2, 30] keV and a power-law SED ${{ \mathcal L }}_{\nu }\propto {\nu }^{-1.5}$. For the Lyα intensity and the LW intensity, we take a simple distinction between Population II and Population III stars. Population III stars are assumed to have N_ion = 50,000, N_{α L} = 4800, and N_{β L} = 2130, where N_{α L} and N_{β L} are the number of photons emitted by a stellar baryon during the stellar lifetime in the energy range from Lyα to LL and from Lyβ to LL, respectively. Population II stars are assumed to have N_ion = 6000, N_{α L} = 9690, and N_{β L} = 3170. This makes N_{α L} /N_ion and N_{β L} /N_ion of Population III stars about an order of magnitude smaller than those of Population II stars, respectively. N_{α L} /N_ion and N_{β L} /N_ion strongly affect J_LW (Equation (16)) and N_α (Equation (24)) for a given PPR. We use the following SED conventions for each category of models:

• 1.  
  Vanilla model: Population II SED.
• 2.  
  SRI model: Population III SED for LMACH; Population II SED for HMACH.
• 3.  
  SRII model: Population III SED for LMACH and MH; Population II SED for HMACH.

The claimed detection of a ∼500 mK absorption dip around ν ≃ 78 MHz by EDGES has been a matter of debate, mainly due to the fact that it is impossible to explain such a large amplitude in the standard ΛCDM framework, if the background after a successful foreground removal is composed only of the CMB and the 21 cm background. In the ΛCDM universe, the kinetic temperature of the IGM is limited to the adiabatically cooled value (T_K ∼ 10.2 K[(1 + z)/21]²), and even at the maximum Lyα coupling, $| \delta {T}_{b}|$ is limited to $\sim 200\,{\rm{mK}}$. Another difficulty faced by the EDGES result is the existence of a peculiar spectral shape in δ T_b , a flat trough of δ T_b from ν = 72 to 85 MHz and lines connecting to the ends of the trough from ν = 65 and 92 MHz, which is in contrast with a smooth dip predicted by models in the ΛCDM.

We show our model predictions on δ T_b for a selected set of models, and compare these with the EDGES result. The selection criterion is the goodness of model fits to the PLD, and we use those minimum-χ²/ν models marked by arrows in Figure 7. We also tune f_X to produce the largest absorption dip for each model but under the condition δ T_b > 0 at z ≲ 9, to (1) comply with the EDGES result with the deepest absorption possible and (2) compensate for our ignorance of the Lyα heating, which might naturally turn the 21 cm background into emission before the end of reionization (Chuzhoy & Shapiro 2007; Ciardi & Salvaterra 2007; Mittal & Kulkarni 2021 ⁶ ) and some hints of the IGM heating at z ∼ 9 (Monsalve et al. 2017; Singh et al. 2018; Ghara et al. 2020; Mertens et al. 2020). Chuzhoy & Shapiro (2007) first showed that the Lyα-recoil heating can solely increase T_K beyond T_CMB before reionization is completed, correcting the estimate by Chen & Miralda-Escudé (2004). In this way, we show how far off each reionization model is from the EDGES result even when the maximum absorption is achieved in each model. One can be more inclusive in model selection because future CMB observations will probe the CMB polarization with better accuracy; nevertheless we stick to this choice here.

Different model categories show distinctive features in δ T_b (Figures 10 and 11, with the shaded regions indicating the redshift bin of the EDGES absorption trough of ∼500 mK), as follows.

• 1.  
  The Planck-favored vanilla models show the familiar ∼180 mK absorption dip. The moment of the absorption dip and the start of the absorption due to the Lyα pumping (to be distinguished from the absorption due to collisional pumping at z ≳ 25) are delayed in the dF case compared to the F case. The dip resides at z ≃ 16 and 14.5 for the dF and F cases, respectively. The start of the absorption is at z ≃ 26 and 21 for the dF and F cases, respectively.
• 2.  
  The Planck-favored SRI models show weaker absorption dips, with $\delta {T}_{b,\min }\simeq [-120,-140]\,\mathrm{mK}$, than the vanilla models. The dF case shows delays in the moments of the absorption dip and the start of the absorption compared to the F case, just as in the vanilla models. Compared to the vanilla models, both the dip and the start of the absorption are delayed: the dip is at z ≃ 12 (F)–13 (dF), and the absorption starts at z ≃ 18 (F)–20 (dF).
• 3.  
  The Planck-favored SRII models show a very slowly deepening absorption "slope" during 28 ≳ z ≳ 14, which is the epoch that is about the same as the full EDGES-low observational window, before the absorption dip at z ≃ 11 − 12 occurs. The absorption dip is with $\delta {T}_{b,\min }\simeq [-90({\rm{L}}0),-120(\mathrm{LL})]\,\mathrm{mK}$, located away from the EDGES trough window. Where the EDGES trough exists, the models have the limited differential brightness temperature, δ T_b ≃ [ − 30, − 60] mK.

Several details of these results are noteworthy. (1) The ∼ 180 mK $\left|\delta {T}_{b,\min }\right|$ of the vanilla model is roughly the maximum amplitude allowed in the ΛCDM cosmology (Furlanetto 2006; Mirocha 2014; Bernardi et al. 2015). (2) The main cause for the overall delays in the dip and the start of the absorption for the SRI model relative to the vanilla model is the fact that LMACHs, which dominate the early phase of reionization, are assumed to host Population III stars. Because ${({N}_{\alpha L}/{N}_{\mathrm{ion}})}_{\mathrm{Pop}\mathrm{III}}$ is much smaller than ${({N}_{\alpha L}/{N}_{\mathrm{ion}})}_{\mathrm{Pop}\mathrm{II}}$ at a given level of x, N_α , and x_α of the SRI models, while LMACHs dominate, they also become much smaller than those of the vanilla models that use the Population II SED. Such a smallness of x_α then renders T_S to couple only weakly to T_K, bringing down the amplitude of δ T_b (Figure 10). (3) For the reason same as in detail (2), had we used the Population III SED for the vanilla model, we would have gotten delays in the dip and the start of the absorption just similar to the SRI case. Accordingly, the amplitude of the absorption dip would have been reduced. Furlanetto (2006) tested the vanilla model using both the Population II and Population III SEDs to observe this tendency. (4) The Lyα pumping efficiency of SRII in the EDGES absorption trough window is almost constant at the restricted value x_α ≃ 0.1 − 0.5, and x_α is bound to < 1 at z ≳ 14. This is mainly caused by the combined effect of the strongly self-regulated SFRD by the LW feedback and the smallness of ${({N}_{\alpha L}/{N}_{\mathrm{ion}})}_{\mathrm{Pop}\mathrm{III}}$, both being relevant to MHs that dominate this epoch.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Compared to the signal interpreted by the EDGES team, none of the global δ T_b 's of our models can match that of EDGES in amplitude and spectral shape. Among the tested models, V-L_dF provides δ T_b closest to the EDGES data. However, this is as close as one can get to the amplitude of the EDGES absorption trough in the ΛCDM, because the case is tuned to produce the coldest T_K and the largest x_α , with x_α ≃ 5–20 in the EDGES window, with a reionization history well under the PLD constraint. The SRI models are terrible in matching the EDGES data, mainly due to the shift of the absorption dip caused by the nature of the Population III SED of LMACH stars dominating this era. The SRII models are as bad as the SRI models in matching the EDGES data within the trough window, mainly due to the dominance of MH stars with Population III SEDs. Note that all of the cases are tuned to produce maximum possible $\left|\delta {T}_{b}\right|$ in absorption with ${f}_{X}^{{\rm{L}}}=0$ (SRI) and ${f}_{X}^{{\rm{L}}}={f}_{X}^{{\rm{M}}}=0$ (SRII) with nonzero ${f}_{X}^{{\rm{H}}}$. Any addition of nonzero X-ray heating will increase T_K from these null-heating cases to reduce $\left|\delta {T}_{b}\right|$ and worsen the mismatch between the EDGES' interpretation and the theory.

The SRII models produce δ T_b that is the most peculiar in spectral shape, because across the full EDGES window (14 ≲ z ≲ 24), δ T_b (z) is almost featureless without much variation. It would even be possible that the process of foreground removal, utilizing the spectral smoothness of the foreground, from the observed signal could completely remove the true EoR (or dark ages) signal at z ≳ 14 or ν ≲ 95 MHz to yield only a null result. The limited amplitude of the absorption depth, $\left|\delta {T}_{b}(z\gtrsim 14)\right|\lesssim 60\,\mathrm{mK}$, and the featureless spectral shape designate our SRII model category as the one that disagrees with the EDGES data most. Considering the excellent agreement of many two-phase reionization models in the SRII category with the PLD polarization data, rather extreme alternative explanations to the standard model are required to explain the EDGES data in order to accept SRII. If SRII were the right model, the excess radio background or other alternatives should (1) almost solely contribute to the absorption trough of ∼500 mK because even the small (≲60 mK) absorption signal is likely to be removed in the foreground-removal process and (2) offset the possible absorption depth of δ T_b ∼ − 100 mK at z ∼ 12, in case X-ray heating is inefficient. Obviously, independent observations such as the 21 cm intensity mapping by radio interferometers will help to settle this issue.

We note that the limited amplitude and the featureless spectral shape of the global δ T_b (z) in the SRII models are a novel result, and such features are in large disagreement with other studies that also implement the LW feedback on Population III stars inside MHs to investigate its impact on δ T_b (Mirocha et al. 2018; Mirocha & Furlanetto 2019; Mebane et al. 2020; Qin et al. 2020a, 2021). Let us explain the major reason for such a discrepancy. A big difference lies in this work and others in the manner that LW feedback is implemented. As described in Section 2.2.3, we assume that there exists a threshold value of J_LW such that star formation inside MHs is fully suppressed as long as J_LW > J_LW,th, based on the observed abrupt change in ${M}_{\min }$ (the minimum mass of star-forming halos) and ${T}_{\mathrm{vir},\min }$ (the minimum virial temperature of star-forming halos) of star-forming MHs as J_LW varies across the value J_LW ∼ 0.1 × 10⁻²¹ erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹ in Yoshida et al. (2003) and O'Shea & Norman (2008). In contrast, other studies usually adopt a prescription where ${M}_{\min }$ and consequently ${T}_{\mathrm{vir},\min }$ are smooth functions of J_LW. While numerical coefficients vary somewhat in the literature (e.g., see difference between Fialkov et al. 2013 and Schauer et al. 2020), the commonly used functional form for ${M}_{\min }$ is either a redshift-independent one (advocated by Machacek et al. 2001 and Wise & Abel 2007),

$$\begin{eqnarray}&&{M}_{\min }/{M}_{\odot }=2.5\times {10}^{5}\left[1+6.8{\left(4\pi {J}_{21}\right)}^{0.47}\right]\end{eqnarray} \tag{ 33 }$$

where J₂₁ ≡ J_LW/(10⁻²¹ erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹), or a redshift-dependent one (Fialkov et al. 2013; Visbal et al. 2020; Qin et al. 2021),

$$\begin{eqnarray}&&\begin{array}{l}{M}_{\min }/{M}_{\odot }\\ =\,2.5\times {10}^{5}{\left(\displaystyle \frac{1+z}{26}\right)}^{-1.5}f({v}_{\mathrm{bc}};z)\left[1+6.96{J}_{21}^{0.47}\right]\end{array}\end{eqnarray} \tag{ 34 }$$

where f(v_bc; z) ≥ 1 is an additional factor accounting for the suppression of star formation due to the baryon–dark matter streaming velocity (Tseliakhovich & Hirata 2010; Ahn 2016). If one uses one of Equations (33) and (34), the effective J_LW,th is much larger than J₂₁ = 0.1 and thus star formation inside MHs will become much stronger than our prescription at a given J_LW. Therefore, studies adopting Equation (34), e.g., Fialkov et al. (2013), Qin et al. (2021), and Visbal et al. (2020), find Lyα intensity much stronger than our prediction, x_α (z ≳ 14) ≲ 0.5, in SRII models. Consequently, these studies find that MH stars cause the absorption dip to occur much earlier at z ∼ 16–18 than in cases without MH stars (Mebane et al. 2020; Qin et al. 2020a) and with amplitude easily reaching $\left|\delta {T}_{b}\right|\gtrsim 100\,\mathrm{mK}$. However, we stress that if one focuses on ${T}_{\mathrm{vir},\min }$, O'Shea & Norman (2008) clearly shows ${T}_{\mathrm{vir},\min }\simeq 8000\,{\rm{K}}$ when J₂₁ ≃ 0.1. If we take this result and extrapolate to any other redshift, one can instead conclude that J_LW,th should lie around J₂₁ = 0.1. In this paper, we parameterize J_LW,th but to a limited value of J_LW,th,21 ≤ 0.1, respecting the results of Yoshida et al. (2003) and O'Shea & Norman (2008). Of course, one cannot exclude one LW feedback scheme against another at the moment, because there is practically no observational constraint on ultra-high-z (z ≳ 14) radiation sources.

The success of our SRII models in producing two-stage reionization models, which are somewhat favored by PLD against the vanilla and SRI models, can be taken as a hint that SRII models may represent the reality, including how the LW feedback operates in nature. If this were true, the gross disagreement of the global δ T_b (z) of SRII models with the EDGES result would be hardly conceivable in the standard model, or the EDGES result could be an incorrect claim. Further study is warranted.