Large-scale Variation in Reionization History Caused by Baryon–Dark Matter Streaming Velocity

The Epoch of Reionization (EOR) is an exciting test site for theoretical models of structure formation. The evolution of the global ionization fraction of the intergalactic medium (IGM) is driven by competition between the production of ionizing photons from galaxies and the consumption of the photons by the intergalactic gas. Observational data from this era has accumulated over the last two decades, which have enable us to constrain the luminosity function of early galaxies (e.g., McLure et al. 2013; Finkelstein et al. 2015; Bouwens et al. 2015).

The first direct constraint on the reionization history came from the Gunn–Peterson (GP) trough feature in quasar spectra, where the Lyα opacity of the IGM rises steeply at z ∼ 6 marking the end of reionization (Songaila 2004; Fan et al. 2006). Measurements of the anisotropies of the cosmic microwave background (CMB) provide constraints on the integrated Thompson optical depth of the IGM. The latest results from the Planck Collaboration imply that reionization was in the midway at z = 7.7 ± 0.7 (Planck Collaboration et al. 2020). Also, the observed decline in the number density of Lyα emitters (Fontana et al. 2010; Pentericci et al. 2011; Schenker et al. 2012; Schmidt et al. 2016; Konno et al. 2018; Jung et al. 2018; Mason et al. 2018; Jung et al. 2020) and detailed analysis of the spectra of high-redshift quasars (Mortlock et al. 2011; Bañados et al. 2018; Davies et al. 2018) constrain the IGM neutral fraction, which provides indirect information on the nature of the sources of reionization (e.g., Robertson et al. 2013, 2015; Finkelstein et al. 2019).

The main interest of the present study is to explain the apparent variation in the Lyα opacity repeatedly reported by observations of distant quasars (Fan et al. 2006; Gallerani et al. 2008; Bolton et al. 2011; Becker et al. 2015). While it has been often regarded that the reionization ended before z = 6, some sightlines show a GP trough stretching well below that redshift, suggesting a substantial variation in the Lyα opacity at ∼50 Mpc scales. The most extreme case found is a ∼150 Mpc long trough stretching from z ≈ 6–5.5 in the sightline of quasar ULAS J0148 + 0600 reported by Becker et al. (2015). It is not understood how this quasar could coexist with other quasars that do not show GP troughs up to z = 6 because it is difficult to explain the large-scale fluctuations in the ultraviolet (UV) background and IGM temperature inferred from the observations.

The attempts to explain the Lyα opacity variation generally fall into two categories. One is to assume that reionization completed before z = 6 and the opacity variation comes from the variations in the IGM temperature (D'Aloisio et al. 2015) or in the UV background (Davies & Furlanetto 2016; D'Aloisio et al. 2017; Chardin et al. 2017; D'Aloisio et al. 2018). The other is to relax the assumption that reionization finished by z = 6 and consider the possibility that reionization ended later in void regions with low galaxy density (Kulkarni et al. 2019; Keating et al. 2020; Nasir & D'Aloisio 2020). Interestingly, Nasir & D'Aloisio (2020) have shown that both approaches could reproduce large-scale opacity variation under certain conditions. Montero-Camacho et al. (2019) and Montero-Camacho & Mao (2020) suggested that the variation in the IGM temperature leaves a measurable imprint on the Lyα forest at z ∼ 2, which will be measured by the Dark Energy Spectroscopic Instrument survey.

In this study, we explore the possibility that the relative streaming motion between baryons and dark matter causes the observed Lyα opacity variation by generating large-scale inhomogeneity in reionization. The streaming motion is a remnant of the baryonic acoustic oscillation (BAO), which froze out at the cosmic recombination with its amplitude around ∼30 km s⁻¹, ranging from 0–70 km s⁻¹ (Tseliakhovich & Hirata 2010). It is known to hinder the formation of Population III stars in minihalos (MHs) at the Cosmic Dawn (z ∼ 20–30) (Greif et al. 2011; O'Leary & McQuinn 2012; Naoz et al. 2013; Richardson et al. 2013; Asaba et al. 2016; Schauer et al. 2019) and generate large-scale imprints on the 21 cm fluctuations from the epoch (McQuinn & O'Leary 2012; Visbal et al. 2012; Muñoz 2019). However, its impact on the EOR (z ∼ 6) has been considered unimportant because the streaming motion does not affect the galaxies in atomically cooling halos with M > 10⁸ M_⊙, which produce most of the ionizing photons during the reionization era (Stacy et al. 2011; Fialkov et al. 2014).

Here, we pay attention to the photon consumption in structures below the Jeans mass of photoionized gas (∼10⁸ M_⊙). The recombination rate of photoionized gas during reionization—as sometimes expressed by a clumping factor by which the average recombination rate is enhanced to account for unresolved density inhomogeneity—has long been suggested to depend upon the inclusion of inhomogeneity on the smallest scales between the Jeans-filtered scale of the photoionized gas following its reionization and the Jeans-filtered scale of the pre-reionization gas (Shapiro et al. 2003, 2004). This has been demonstrated in a series of papers, starting from the first rad-hydro simulations of the photoevaporation of MHs by the intergalactic I-fronts that propagate during the EOR (Shapiro et al. 2004; Iliev et al. 2005b) and the application of those results to model the impact of this process on the progress of reionization (Iliev et al. 2005a).

Meanwhile, the clumping factor and the enhancement of the recombination rate associated with this small-scale structure in a three-dimensional cosmological density field, in the filamentary structure outside of MHs, as well, was subsequently described by Emberson et al. (2013) in an approximate way, by taking a snapshot of a tiny, highly resolved cosmological N-body + hydro simulation without radiative transfer or ionizing radiation, and then post-processing the static gas density field by introducing ionizing radiation, and relaxing the outcome to solve for the "Strömgren surface" separating the ionized exteriors of halos from the self-shielded interiors of halos. This work established the demands for spatial and mass resolution required to make the simulated value of the total recombination rate in the ionized regions converge.

Following that static calculation, Park et al. (2016, hereafter P16), developed the Gadget-RT code to perform this calculation with fully coupled radiation hydrodynamics. The results showed that the impact of a small-scale gas structure (SSGS) on the global progress of reionization could be significantly affected by the additional recombinations that result from these small-scale density fluctuations, including both the MHs and the fluctuating IGM around them, and parameterized the effect for application to future models. Recently, this problem has been revisited by D'Aloisio et al. (2020, hereafter D20), who confirmed the importance of the SSGS contribution in the ionization photon budget of reionization and showed a quantitative impact on the global reionization history. This small-scale recombination can be implemented in large-scale reionization simulations as a sub-grid correction, and the predictions for reionization observable are affected significantly by the correction (Mao et al. 2020).

In the meantime, Hirata (2018) showed that the streaming motion disrupts a large portion of SSGSs at z ∼ 6 despite its characteristic decay over time. Thus, our purpose here is to explore the impact of this effect on the clumping factor and recombination rate boost reported by P16 and D20. In particular, we focus our interest here on the possibility that the effect is modulated on large scales, correlated with the BAOs.

Another potentially significant effect is the pre-reionization heating by X-ray background. During the reionization, X-rays from accreting black holes are expected to have heated the IGM to some level (Ricotti & Ostriker 2004; Eide et al. 2018). The formation of SSGSs in the pre-reionization IGM is largely owed to its low gas temperature. Thus, intense X-ray preheating would likely suppress SSGSs. The exact value of the gas temperature in the pre-reionization phase is highly uncertain. Thus, we shall consider it as a free parameter between 10 and 1000 K in this work.

Therefore, our goal is to assess the impact of the streaming motion and the X-ray preheating on the recombination by the SSGSs and the resulting reionization history using a series of hydrodynamic simulations developed by P16. In a similar work, Cain et al. (2020, hereafter C20) investigated the signature of the streaming motion in the 21 cm signal. This study focuses on the impacts on the reionization history.

The rest of the paper is organized as follows. In Section 2, we explain how we initialize our simulations with and without the streaming motion. In Section 3, we explain how we evolve the initial conditions (ICs) of the IGM until it is exposed to the ionizing radiation. In Section 4, we report the simulated results on the clumpiness of the IGM after being exposed to ionizing background radiation. In Section 5, we calculate the reionization history from the clumping factor measured from our simulations to assess the impact of the streaming velocity and X-ray preheating on reionization. In Section 6, we summarize and discuss our results. For the rest of this paper, we assume a Lambda cold dark matter cosmology consistent with the Wilkinson Microwave Anisotropy Probe 9 yr results (Hinshaw et al. 2013): Ω_m,0 = 0.276, Ω_b,0 = 0.045, h = 0.703, σ₈ = 0.8, and n_s = 0.961.

We create cosmological ICs at z = 200 with 256³ particles for gas and dark matter each in a 200 h⁻¹ kpc box. In this setup, the gas and dark matter particle masses are 44 and 8.4 h⁻¹ M_⊙, respectively. P16 showed that this mass resolution well captures the internal gas structure of MHs during the evaporation. In a similar work, D20 also arrived at a similar conclusion in their convergence test (see appendix C of their work).

The simulation box size (200 h⁻¹ kpc) in this study is not large enough to capture the large-scale variation in density. The local overdensity of such a small volume can deviate substantially from zero, and P16 and D20 showed that the recombination rate depends sensitively on the local overdensity. The box size and resolution test in Appendix C of C20 showed that our box size would miss large-scale contributions to the clumping factor, but they are not affected much by the streaming motion. We limit the scope of this study to reporting the qualitative impacts of the streaming motion on the recombination in a mean-density volume.

With the setups described above, we create three sets of ICs with three different streaming velocity levels, V_bc,i = 0, 28, and 70 km s⁻¹. Since the streaming motion is highly coherent over a few-megaparsec scale, we implement the motion as a constant drift between gas and dark matter within our simulation box. The probability distribution of the streaming velocity at z = 1000 mostly ranges between 0 and 70 km s⁻¹, with its average near ∼30 km s⁻¹ (see Figure 2 of Ahn 2016). Thus, our choices should cover most of the possible values of V_bc,i .

To generate the ICs, we use the Baryon–CDM Cosmological Initial Condition Generator for Small Scales (BCCOMICS) package of Ahn & Smith (2018). ⁷ BCCOMICS calculates the perturbation equations for the density/velocity fluctuation amplitude with the streaming velocity terms to account for the streaming effect in a self-consistent manner. The details of the calculation are described in Ahn (2016) and Ahn & Smith (2018). We note that the streaming motion is often implemented in an approximate manner by adding a constant drift velocity to the gas density/velocity generated without accounting for the streaming motion. In this case, the ICs do not include the smoothing effect of the streaming motion on the gas density field taking place before the initialization of the simulation, substantially overestimating the gas density fluctuation amplitude at small scales (Park et al. 2020).

We show the initial density/velocity fields of the ICs in Figure 1. The same set of random phases is used for both sets of ICs to prevent the cosmic variance from affecting our analysis. In the case with a streaming velocity of 70 km s⁻¹, the gas is moving toward the lower right corner of the figure with respect to the underlying dark matter. The dark matter density is nearly unaffected by the streaming, while the gas density shows some smoothing effect toward the streaming direction.

Zoom In Zoom Out Reset image size

We evolve the initial gas/dark matter particle distribution at z = 200 using the publicly available smoothed particle hydrodynamics (SPH) code GADGET-2 (GAlaxies with Dark matter and Gas intEracT; Springel et al. 2001; Springel 2005) ⁸ down to z = 5.5. In this step, the gas cools adiabatically as expected for the IGM before reionization.

These treatments are based upon the assumption, justified by Haiman et al. (2000), that the typical MH was not yet forming stars before the external I-front overtook it. The rise of the UV ionizing radiation background responsible for the I-front that overtook each MH was led by an earlier rise of the Lyman–Werner background to the threshold level required to dissociate the H₂ molecules inside MHs to suppress their star formation.

We first simulate structure formation with a minimum gas temperature of ${T}_{\min }=10$ K down to z = 20. From z = 20, we run four cases with different T_min's to simulate the impact of different levels of preheating by X-ray sources. ⁹ In one case, we leave T_min at 10 K after z = 20, assuming a minimal level of preheating. In the other three cases, we raise T_min to 100, 300, and 1000 K. In combination with three levels of streaming velocity, V_cb,i = 0, 28, and 70 km s⁻¹, we create 12 cases of pre-reionization IGM in total. Figure 2 describes the time evolution of the density field for V_cb,i = 0 and 70 km s⁻¹.

Zoom In Zoom Out Reset image size

The density maps in Figure 2 show the lasting impacts of the streaming velocity and X-ray preheating on gas density. The case with zero streaming velocity and minimal preheating (${T}_{\min }=10$ K) depicted in the upper side of the upper panel is the most abundant in structures, and the cases with a higher V_cb,i or higher T_min's yields fewer structures. Notably, the streaming motion makes a substantial difference even at z ≲ 6, when it has decayed below 1 km s⁻¹. This result may appear to contradict previous studies that found that the streaming motion has a minor impact on halos that host star-forming galaxies during reionization. The mass scale for the SSGSs is below that for those galaxies (∼10⁸ M_⊙), and they are still susceptible to the streaming effect. Hirata (2018) reported similar findings from similar simulations.

We show gas density power spectra of four cases in Figure 3 to describe the impacts of streaming and preheating more quantitatively. Both suppress the small-scale end of the power spectrum, but the detailed response is somewhat different between the two factors. The suppression by preheating is more focused on small scales above a particular wavenumber as the increased T_min raises the filtering mass scale, below which the structures are removed (compare solid versus dashed lines). On the contrary, the suppression by the streaming motion appears in a wider range of scales (compare blue versus cyan lines).

Zoom In Zoom Out Reset image size

The EOR is thought to have happened in an extended period, presumably starting from density peaks and finishing void regions in the final stage. A small volume like our simulation box will be exposed to the ionizing background at different timings, depending on where it is located in the universe. The recombination rate is expected to depend on that timing because the SSGSs would gravitationally evolve over time. Thus, we treat the reionization redshift of our simulation box (z_re) as an extra parameter for the SSGS recombination rate and save the snapshots at ${z}_{\mathrm{re}}=5.5,\,6,\,7,\,9$, and 12 for each choice of T_min and V_cb,i . Therefore, we prepare 60 snapshots of pre-reionization IGM in total to simulate the impact of reionization on SSGSs.

We track the fate of SSGSs after the onset of the external ionizing background radiation (EIBR) in the snapshots produced from the previous steps. We use the GADGET-RT (GADGET with Radiative-Transfer) code developed from P16, where the EIBR is simulated with a self-shielding scheme and a nonequilibrium chemistry of the hydrogen and helium-related species described in Yoshida et al. (2006, 2007). Among the chemical reactions included in the code, the ionization (H → H⁺ + e⁻) and recombination (H⁺ + e⁻ → H) of hydrogen are relevant to our calculations. We turn on the blackbody radiation of 10⁵ K, with a photoionization rate of 9.2 × 10⁻¹² cm⁻² s⁻¹ Hz⁻¹ sr or Γ₋₁₂ = 9.2 in the pre-reionization IGM generated from the previous steps. The spectral energy distribution and radiation intensity of EIBR have been shown to have a minor impact on the evolution of SSGSs in P16, D20, and C20.

For each chemistry time step, we first list all the gas particles of which hydrogen is more than 50% neutral. For each of those target particles, we locate the neighboring particles within 200 physical parsecs utilizing the tree structures prepared for the gravity solver. We, then, calculate the average neutral hydrogen column density for six directions, ±x, ±y, and ±z, for each target particle. This way, we ionize neutral regions from the outskirt while shielding the inner part of the gas clump properly. P16 has shown that this radiative-transfer algorithm reproduces physical quantities of evaporating MH gas with high accuracy. The MH gas is practically the only neutral gas that is shielded against the EIBR for a significant amount of time in the simulation (>10⁶ yr).

Figure 4 depicts the evolution of gas density after the EIBR is turned on at ${z}_{\mathrm{re}}=6$ for three cases. Upon reionization, the Jeans mass increases to ∼10⁸ M_⊙, and the SSGSs begin to be destroyed by the increased pressure. Most of the structure is destroyed after a few times 10⁶ yr, while gas clumps in the MHs last up to ∼10⁸ yr due to self-shielding.

Zoom In Zoom Out Reset image size

P16 and D20 have shown that SSGSs substantially increase the ionization photon budget for reionization by raising the recombination rate. In the following section, we describe how we calculated the recombination rate from our simulation and how it depends on z_re, T_min, and V_cb,i .

4.1. Recombination Rate

We briefly summarize how we calculate the recombination rate from our simulations and report the results. For the fully detailed derivations, we refer the readers to Section 3 of P16.

The average recombination rate in a volume V is given by

$$\begin{eqnarray}&&\bar{{ \mathcal R }}={{ \mathcal C }}_{\mathrm{sim}}{\alpha }_{B}(\bar{T}){\bar{n}}_{e}{\bar{n}}_{\mathrm{HII}},\end{eqnarray} \tag{ 1 }$$

where ${\alpha }_{{B}}{(T)=2.6\times {10}^{-13}(T/{10}^{4}{\rm{K}})}^{-0.7}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-3}$ is the case B recombination coefficient at temperature T, ${\bar{n}}_{X}\equiv {\left\langle {n}_{X}\right\rangle }_{V}$ is the volume weighted averaged number density of a species X, and

$$\begin{eqnarray}&&{{ \mathcal C }}_{\mathrm{sim}}=\frac{{\left\langle {\alpha }_{B}(T){n}_{e}{n}_{\mathrm{HII}}\right\rangle }_{V}}{{\alpha }_{B}(\bar{T}){\bar{n}}_{e}{\bar{n}}_{\mathrm{HII}}}\end{eqnarray} \tag{ 2 }$$

is the clumping factor, which accounts for the boost in the recombination rate due to the spatial fluctuations in T, n_e , and n_HII. Here, the subscript S denotes that the results from our simulations only include clumping in SSGSs. The mean temperature $\bar{T}$ is given by

$$\begin{eqnarray}&&\bar{T}=({m}_{p}/{k}_{B}){\left\langle (\gamma -1)u\right\rangle }_{M}{\left\langle {\mu }^{-1}\right\rangle }_{M}^{-1},\end{eqnarray} \tag{ 3 }$$

where u is the specific internal energy, μ is the mean molecular weight, ${\left\langle \right\rangle }_{M}$ denotes the mass-weighted average. We note that, as was reported in P16, our definition of the clumping factor gives about a 10% higher recombination rate than the conventional one, $\left\langle {n}^{2}\right\rangle /{\left\langle n\right\rangle }^{2}$, does because a higher cooling rate in a higher density makes α_B (T) mildly correlate with density.

Nearly 100% of the hydrogen atoms are ionized shortly after the onset of EIBR except for the dense MH gas, which takes only a small fraction of total gas giving n_HII ≈ n_H and helium atoms are singly ionized giving 8% of extra free electrons: n_e ≈ (1.08)n_H. Thus, the evolution of the clumping factor is the major factor that dominates the recombination rate. In the case of the SPH particle data, the clumping factor can be computed from the following equation:

$$\begin{eqnarray}&&{{ \mathcal C }}_{\mathrm{sim}}={\bar{n}}^{-1}{N}_{\mathrm{ptl}}\left[\frac{{{\rm{\Sigma }}}_{i}{f}_{e,i}{f}_{\mathrm{HII},i}{\alpha }_{B}({T}_{i}){n}_{i}}{\left({{\rm{\Sigma }}}_{i}{f}_{e,i}\right)\left({{\rm{\Sigma }}}_{i}{f}_{\mathrm{HII},i}\right){\alpha }_{B}(\bar{T})}\right].\end{eqnarray} \tag{ 4 }$$

Here, the subscript i denotes the ith SPH particle in the simulation, N_ptl is the number of SPH particles, n ≡ ρ/m_p is density in the unit of the proton mass, f_X ≡ n_X /n is the number density of a species X divided by n, and $\bar{T}$ is given by averaging over the particle values: ${N}_{\mathrm{ptl}}^{-1}{{\rm{\Sigma }}}_{i}{T}_{i}$.

Due to cosmic expansion, volumes with the same clumping factor yield different recombination rates at different redshifts. To compare the recombination rate from different z_re's, it is convenient to weight the clumping factor by the cosmic mean density $\bar{\rho }$. Specifically, multiplying the clumping factor by $\bar{\rho }$ normalized at its z = 7.5 value (${\hat{{ \mathcal C }}}_{\mathrm{sim}}\equiv [\bar{\rho }/{\bar{\rho }}_{z=7.5}]{{ \mathcal C }}_{\mathrm{sim}}$) makes ${\hat{{ \mathcal C }}}_{\mathrm{sim}}=1$ correspond to the recombination rate of 1 per H atom per gigayear. We plot ${\hat{{ \mathcal C }}}_{\mathrm{sim}}$ for nine cases in Figure 5.

Zoom In Zoom Out Reset image size

On a log–log axis (left panel of Figure 5), the clumping factor as the function of time shows a rise and decay as reported in P16 and D20. Initially, the ionization fronts are in the supersonic R-type phase, where the ionization proceeds toward dense gas clumps before the gas responds to the photoheating. The clumping factor rises as denser gas is ionized. After several megayears, ionized gas begins to expand and lower the clumping factor. On the linear–linear axis (right panel of Figure 5), the early R-type phase is negligibly short and unimportant for the integrated recombination, which is proportional to the area under the curve.

Comparison of ${\hat{{ \mathcal C }}}_{\mathrm{sim}}$ in the nine cases shown in Figure 5 confirm that the clumping factor depends sensitively on three parameters, z_re, V_cb,i , and T_min: ${\hat{{ \mathcal C }}}_{\mathrm{sim}}$ is higher for lower z_re's, lower V_cb,i 's, and lower T_min's. These parameter dependences have been partially shown in previous studies. D20 ran one case with ${T}_{\min }=1000$ K to find a much lower clumping factor than in the ${T}_{\min }=10$ K case, C20 showed that the clumping factor is suppressed when the streaming velocity is higher, and these studies repeatedly showed the z_re dependence. This study focuses on computing the variation in the reionization history due to the spatial fluctuations in V_cb,i and T_min.

Here, we calculate the reionization history from a one-zone model similar to the one used by D20 to assess the impacts of the streaming motion and the X-ray preheating. Due to the long fluctuation scale of the streaming motion (∼150 Mpc), we expect volumes that are ≲50 Mpc across to form a zone where V_cb,i is more or less uniform. We intend to describe the difference in reionization histories in such zones with different local values of V_cb,i .

We emphasize that our results from the one-zone model should not be taken to be quantitatively accurate as it leaves out most of the complexities of reionization. For example, the spatial inhomogeneity of reionization due to clustering of ionizing sources cannot be modeled by design. We defer more accurate modelings of reionization to future studies with more sophisticated numerical simulations.

5.1. One-zone Model

The change in the global mean ionization fraction of the IGM (Q_HII) is given by the volume-averaged difference between the ionization rate ${ \mathcal I }$ and the recombination rate ${ \mathcal R }$:

$$\begin{eqnarray}&&\frac{{{dQ}}_{\mathrm{HII}}}{{dt}}={\bar{n}}_{{\rm{H}}}^{-1}(\bar{{ \mathcal I }}-\bar{{ \mathcal R }})\end{eqnarray} \tag{ 5 }$$

(Shapiro & Giroux 1987; Madau et al. 1999). The volume average here involves the entire universe, unlike in the previous section where we only averaged over the simulation box. In this case, Equation (1) is slightly changed to

$$\begin{eqnarray}\begin{array}{rcl}\bar{{ \mathcal R }} & = & {{ \mathcal C }}_{\mathrm{tot}}{\alpha }_{B}(\bar{T}){\bar{n}}_{e}{\bar{n}}_{\mathrm{HII}}\\ & = & {{ \mathcal C }}_{\mathrm{tot}}{\alpha }_{B}(\bar{T})(1.08){\bar{n}}_{{\rm{H}}}^{2}{Q}_{\mathrm{HII}},\end{array}\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal C }}_{\mathrm{tot}} & \equiv & [{{ \mathcal C }}_{L}-1]+[{{ \mathcal C }}_{S}-1]+1\\ & = & {{ \mathcal C }}_{L}+{{ \mathcal C }}_{S}-1\end{array}\end{eqnarray} \tag{ 7 }$$

combines clumping from large scales beyond our simulation box (${{ \mathcal C }}_{L}$) and that from SSGSs in the simulation (${{ \mathcal C }}_{S}$). The large-scale contribution comes from structures above the Jeans mass of ionized gas and can be obtained from large-scale simulations with the box size ≳10 Mpc. In this study, we adopt the fitting function from Pawlik et al. (2009): ${{ \mathcal C }}_{L}=1+43{z}^{-1.71}$. For the global SSGS contribution ${{ \mathcal C }}_{S}$, one must consider that different regions in the universe can be exposed to the EIBR at different z_re's. Thus, we average the simulated results of ${{ \mathcal C }}_{\mathrm{sim}}$ over the past reionization history:

$$\begin{eqnarray}&&{{ \mathcal C }}_{S}(z)=\frac{{\displaystyle \int }_{{z}_{i}}^{z}{\dot{Q}}_{\mathrm{HII}}({z}^{^{\prime} }){{ \mathcal C }}_{\mathrm{sim}}({\rm{\Delta }}{t}^{^{\prime} },{z}^{^{\prime} })(d{\rm{\Delta }}{t}^{^{\prime} }/{{dz}}^{^{\prime} }){{dz}}^{^{\prime} }}{{Q}_{\mathrm{HII}}(z)},\end{eqnarray} \tag{ 8 }$$

where ${\rm{\Delta }}{t}^{{\prime} }=t(z)-t({z}^{{\prime} })$ is the time difference between z and ${z}^{{\prime} }$. We assume that the reionization begins at z = 12 and interpolate the simulated results of ${{ \mathcal C }}_{\mathrm{sim}}({\rm{\Delta }}t,{z}_{\mathrm{re}})$ obtained for ${z}_{\mathrm{re}}=5.5,\,6,\,7,\,9,$ and 12 to evaluate the integral.

For the source term of Equation (5), ${\bar{n}}_{{\rm{H}}}^{-1}\bar{{ \mathcal I }}$, we construct our one-zone models motivated by Iliev et al. (2007) using two classes of ionizing sources, the low-mass sources (≲10⁹ M_⊙) whose formation is suppressed in H ii regions due to Jeans mass filtering and the high-mass sources (≳10⁹ M_⊙) that are not suppressed. Then, the ionization term is given by

$$\begin{eqnarray}&&{\bar{n}}_{{\rm{H}}}^{-1}\bar{{ \mathcal I }}={f}_{H}+(1-{Q}_{\mathrm{HII}}^{n}){f}_{L},\end{eqnarray} \tag{ 9 }$$

where f_L and f_H are the efficiencies of the low- and high-mass sources, respectively. Here, the factor, $1-{Q}_{\mathrm{HII}}^{n}$, accounts for the suppression of low-mass sources in ionized regions. The exponent n would equal 1 if the low-mass sources were homogeneously distributed in space, but, in practice, it is less than 1 because the low-mass sources are highly clustered around density peaks where the reionization begins, and galaxies are more efficiently suppressed at the early time than at the late time of the reionization. Iliev et al. (2007) empirically found that n = 0.1 fits their simulated results. They provide the fitting functions for the f_L and f_H as the following ¹⁰ :

$$\begin{eqnarray}&&{f}_{L}=0.335\left(\displaystyle \frac{{f}_{\gamma ,L}}{250}\right)\exp (0.227z-0.02546{z}^{2})\,{\mathrm{Myr}}^{-1},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{f}_{H}=4.6312\left(\displaystyle \frac{{f}_{\gamma ,H}}{250}\right)\exp (-0.107z-0.02463{z}^{2})\,{\mathrm{Myr}}^{-1}.\end{eqnarray} \tag{ 11 }$$

Here, the z dependence is from the evolving source formation rate over the cosmic history, and f_γ,L and f_γ,H are the efficiency parameters for the two classes of ionizing sources.

We consider three models for the source term ${\bar{n}}_{{\rm{H}}}^{-1}\bar{{ \mathcal I }}$ by modulating f_γ,L and f_γ,H as in Table 1. Descriptions of the models are as follows.

• 1.  
  Accelerating model: Only consider the contribution from high-mass sources (f_γ,E ) , which are not suppressed by ionizing radiation. The source term monotonically increases over time, thereby accelerating the reionization and is qualitatively similar to the Robertson et al. (2015) model.
• 2.  
  Decelerating model: Multiply a factor of (z/7.21)⁴ to the accelerating model to make the sources less efficient at a later time, thereby decelerating the reionization. This is motivated by the "early reionization" scenario of Finkelstein et al. (2019), where the escape fraction of ionizing photons decreases over time.
• 3.  
  Self-regulated model: Includes the contribution from low-mass sources (f_γ,L ), which are suppressed as the global ionization fraction Q_HII rises. This scenario assumes that the low-mass sources have a much higher ionizing radiation escape fraction than the high-mass ones do, allowing them to contribute significantly to the ionizing background (see, e.g., Hutter et al. 2020a).

For each model, the values of f_γ,L and f_γ,H are tuned so that the reionization finishes at z = 6 when V_cb,i = 70 km s⁻¹ and ${T}_{\min }=10$ K.

Table 1. Reionization Source Models and Variation in the End-of-reionization Redshift

Model	fS	fL	ze (Tmin = 10 K)	ze (Tmin = 100 K)	ze (Tmin = 300 K)	ze (Tmin = 1000 K)
Accelerating	0	1560	5.76–6 a	5.84–6.04	5.94–6.08	6.09–6.15
Decelerating	0	1560[z/7.21]4	5.27–6	5.53–6.10	5.83–6.20	6.22–6.33
Self-regulated	15,600	975	5.78–6	5.86–6.04	5.94–6.08	6.09–6.15

Note.

^a The range is given by bracketing the results from the no-streaming case and maximal-streaming case.

Download table as:  ASCIITypeset image

For the sink term ${\bar{n}}_{{\rm{H}}}^{-1}\bar{{ \mathcal R }}$, we consider three "zones" with different recombination rates:

• 1.  
  No-streaming zone: The small-scale clumping factor ${{ \mathcal C }}_{S}$ is obtained from simulations with V_cb,i = 0 km s⁻¹ to represent volumes with a minimal streaming motion.
• 2.  
  Typical-streaming zone: To represent the volumes with typical magnitude of streaming velocity, ${{ \mathcal C }}_{S}$ is obtained from simulations with V_cb,i = 28 km s⁻¹.
• 3.  
  Maximal-streaming zone: ${{ \mathcal C }}_{S}$ is obtained from simulations with V_cb,i = 70 km s⁻¹, which give a maximal impact by the streaming motion.
• 4.  
  No-SSGS zone: A hypothetical case where the SSGS contribution is ignored by setting ${{ \mathcal C }}_{S}=0$.

5.2. Variation in Reionization History Induced by the Streaming Motion

We first consider the minimal X-ray preheating cases with ${T}_{\min }=10\,{\rm{K}}$. We present the clumping factor ${{ \mathcal C }}_{\mathrm{tot}}$ (upper panels), the reionization history Q_HII(z) (middle panels), and the net ionization rate dQ_HII/dt (lower panels) in Figure 6 for the three cases of the sink term with three reionization sources models. We assume that the results from the no-streaming zone and the maximal-streaming zone bracket the possible scatter in reionization history caused by the streaming motion and show that range as the gray shading. The variation in the end-of-reionization redshift z_e for ${T}_{\min }=10\,{\rm{K}}$ is listed in Table 1.

Zoom In Zoom Out Reset image size

The no-SSGS case generally finishes reionization much earlier than the no-streaming zone does. In the decelerating model (Figure 6(b)), for example, reionization ends at z_e ≈ 5.3 in the no-streaming zone, while it ends at z_e ≈ 6.5 in the no-SSGS zone. This highlights the significant contribution from SSGSs to the ionization photon budget, which has been reported by P16 and D20.

The reionization history of the maximal-streaming zone falls near the midpoint between those of the other two cases (z_e = 6), showing that the streaming motion can undo a significant portion of the delay in reionization caused by SSGSs. The main interest of this work is the difference between the no-streaming zone and the maximal-streaming zone, which sets the maximum impact of the streaming motion on the reionization history. The reionization history of the typical-streaming zone falls near the midpoint between those of the no-streaming and the maximal-streaming zones, suggesting that the impact of the streaming velocity is more or less linearly proportional to its initial magnitude V_cb,i .

Since the maximal-streaming zone is tuned to give z_e = 6, z_e in the no-streaming zone sets the scatter in reionization history, which are 5.76, 5.26, and 5.78 in the accelerating, decelerating, and self-regulated reionization models. These are comparable to the z_e = 5.5 inferred from ULAS J0148 + 0600, suggesting that the streaming motion can help explaining the apparent variation in Lyα opacity.

The clumping factor ${{ \mathcal C }}_{\mathrm{tot}}$ in the upper panel of Figure 6 shows what fraction the recombination is suppressed by the streaming motion. The clumping factor contribution from SSGSs in the no-streaming (maximal-streaming) zone is given by the gap between cyan (blue) and black lines. The SSGS contribution in the maximal-streaming zone is roughly half of that in the no-streaming zone in all three source models, showing that the streaming motion cancels about half of the recombination by SSGSs.

5.2.1. Reionization Source Model Dependence

The characteristic of each source model is well described by the ionization source term $\left\langle { \mathcal I }\right\rangle$ (lower panels of Figure 6). The source term in the accelerating model starts at ∼0.1 Gyr⁻¹ and increases steeply with time, while it decreases after peaking at 2 Gyr⁻¹ at z = 8 in the decelerating model. In the self-regulated model, the ionizing efficiency rises slowly, falling between the other two models.

A notable feature of the self-regulated model is that the source term is larger when the sink term is larger, partially canceling out the difference in net ionization rate caused by different amounts of recombination. This is because the "self-regulation" mechanism suppresses more ionizing sources when the global ionization fraction is higher.

The characteristic of each source model explains the resulting reionization history (see middle panels of Figure 6). The global ionization fraction in the accelerating model rises later than in other models, reaching 50% at z ≈ 7, while it does at z ≈ 8 in the other two models. Due to the steep rise of the ionizing efficiency, the global ionization fraction in the accelerating model catches up with that in the other two models near the end of reionization. The reionization history of the accelerating models agrees well with the constraint from Davies et al. (2018) and Whitler et al. (2020), while the other two models agree better with the constraints from Inoue et al. (2018), Planck Collaboration et al. (2020), Whitler et al. (2020), and Jung et al. (2020). The constraint on Q_HII from Hoag et al. (2019) is lower than in all the three models considered in this work. The lower bounds on Q_HII at z < 6 from McGreer et al. (2015) are consistent with all the cases considered in this work.

The variation in reionization history caused by the streaming motion is notably larger in the decelerating model than in other models: z_e ranges between 5.27 and 6 in the decelerating model, while it ranges between ∼5.75 and 6 in the other two models (see Table 1). This shows how sensitively the variation depends on the characteristic of reionization sources. Below, we explain how the source model affects the variation.

We find that a hydrogen ion in the no-streaming zone generally recombines, on average, about 0.15 times more than in the maximal-streaming zone by z = 6. As a result, the global ionization fraction in the no-streaming zone is near 85% in the accelerating model and the decelerating model, while 92% in the self-regulated model, where the self-regulation partially cancels the difference in ionization rate. The steeply rising source term in the accelerating model finishes the remaining neutral gas in the no-streaming zone relatively quickly. The reionization in the self-regulated model is finished as early as in the accelerating model despite the lower ionizing rate because of the lower remaining neutral fraction. On the contrary, the end of reionization in the no-streaming zone is much delayed by the declining ionizing rate after z = 6 in the decelerating model. After z = 6, a few percent of IGM neutral fraction persists to z ≲ 5.3 in the decelerating model, consistent with the constraints in McGreer et al. (2015).

Thus, we conclude that the ionizing photon production rate after z = 6 is a critical factor in the variation in the reionization history caused by the streaming motion as well as the difference in local streaming velocity magnitude. The variation would increase further if we lowered the ionizing efficiency at z < 6 in our source model.

5.3. Impact of X-Ray Preheating

To assess the impact of X-ray preheating, we calculate reionization history for four different minimum gas temperatures of ${T}_{\min }=10,\,100,\,300,$ and 1000 K. The results are shown in Figure 7 for the decelerating source model. The z_e variations for these T_min's are listed in Table 1 for all the three reionization source models considered in this work.

Zoom In Zoom Out Reset image size

As T_min increases, the impact of the streaming motion diminishes: the reionization histories in both the maximal-streaming and no-streaming zones converge toward the no-SSGS case, giving Δz_e = 0.48, 0.38, 0.26, and 0.07 for ${T}_{\min }=10,\,100,\,300,$ and 1000 K, respectively. In the most extreme X-ray preheating case of T_min = 1000 K, the recombination from SSGSs is nearly wiped out.

The effect of the X-ray preheating is similar to that of the streaming motion: suppressing the SSGSs (see Figure 2). Therefore, raising the preheating level globally would only the reduce the impact of the streaming motion.

Since Tseliakhovich & Hirata (2010) pointed out the potentially large impact of the relative streaming motion between baryons and dark matter on early structure formation, a number of studies followed to assess its impact on early galaxy formation and reionization. In this study, we have investigated the possibility of the streaming motion affecting the reionization sink, the small-scale gas clumping in the pre-reionization IGM. We have shown that the streaming motion can erase a significant portion of the clumpiness in the pre-reionization IGM and thus make reionization happen significantly earlier under certain conditions. Because the streaming velocity varies over large length scales (≳100 Mpc), our results suggest that the streaming motion can contribute to the variation in Lyα opacity observed from different quasar sightlines, and help to explain the large scatter from recent observation when combined with fluctuations in the UV background and IGM temperature. A schematic description of the large-scale inhomogeneity caused by the streaming motion is given in Figure 8.

Zoom In Zoom Out Reset image size

A hydrogen ion in a no-streaming (V_cb,i = 0 km s⁻¹) region can recombine ∼0.15 times more than in a highly streaming (V_cb,i = 70 km s⁻¹) region throughout the reionization era. As a result, 15% of the IGM remains neutral in the no-streaming zone when the reionization is already finished in the streaming zone. In our one-zone model for reionization history, we find that the end-of-reionization redshift can vary by up to Δz_e ∼ 0.5, depending on the characteristic of ionizing sources. We find that the variation increases as the ionizing efficiency of sources declines toward the end of reionization. On the other hand, the variation is much smaller when the source efficiency does not evolve, and the ionizing background steeply rises at late time. The variation can also be suppressed if low-mass galaxies are self-regulated due to the photoionization feedback.

X-ray preheating is another critical factor that affects the reionization history. The current lack of knowledge about the pre-reionization IGM temperature adds large uncertainty to the recombination by SSGSs. Raising the minimum gas temperature T_min results in the suppression of SSGSs, and therefore, less recombination by SSGSs. In the most extreme case of T_min = 1000 K, the preheating almost completely suppresses the recombination by SSGSs. However, we note that this conclusion is based on the assumption that the X-ray preheating is spatially homogeneous. The X-ray background may have had a significant level of spatial fluctuations generated by the large-scale distribution of X-ray sources. In that case, the fluctuating X-ray background can add to the variation in z_e . An interesting possibility is that the streaming velocity induces intermediate-mass direct collapse black holes as suggested by Hirano et al. (2017), and the radiation from these black holes constitutes a significant portion of the X-ray background. In such cases, the X-ray background would correlate with the streaming velocity magnitude, possibly enhancing the variation in reionization history. Hence, X-ray preheating is another critical, yet uncertain factor that determines the reionization history.

We emphasize that this study aims to highlight the potential importance of the streaming velocity in a qualitative sense. The numerical values from our one-zone models with our parameter choices are subject to numerous uncertainties. The validity of our assumption that the reionization begins at z = 12 is largely subject to the MH contribution in reionization, which may ionize up to 15% of the IGM by that redshift (see, e.g., Ahn et al. 2012). Also, the two levels of streaming velocity that bracket the entire range, V_cb,i = 0 and 70 km s⁻¹, are unlikely extremes for the mean value of volumes that are ∼50 Mpc across. A more realistic comparison would involve values like V_cb,i ∼5 and 55 km s⁻¹, which would reduce Δz_e by ∼30%.

The results of our small-box simulations pose extra uncertainties due to cosmic variance and fluctuations in local density. C20 found that the impact of the streaming on recombination rate can be several times smaller in larger simulation boxes or with a lower ionizing intensity (e.g., Γ₋₁₂ = 0.3). Our choice of the ionization background intensity, Γ₋₁₂ = 9.2, is more relevant to the early phase of the reionization, where the HII regions exist only in the vicinity of ionization sources. As the HII regions expand to fill the entire space, the mean value should decrease, especially in the decelerating reionization scenario. The box-size test of C20 suggests that the velocity streaming effect tends to be smaller in larger boxes due to massive halos that are less sensitive to the streaming motion. Indeed, our simulation box lacks halos more massive than 10⁷ M_⊙. However, halos with 10⁷ ≲ M ≲ 10⁹ M_⊙ are known to form a certain amount of stars during the pre-reionization phase (e.g., Jeon et al. 2019), violating one of the key assumptions of this study: no star formation in the simulated volume. Thus, neglecting star formation inside these halos would likely overestimate the ionization photon budget. We shall explore this possibility in our future study.

Considering the patchy nature of reionization, which our one-zone model does not account for, the fluctuations in the UV background and IGM temperature caused by the clustering of ionizing sources should also contribute significantly to the Δz_e from quasar sightlines on top of the contribution from the streaming velocity fluctuations. Given that various factors mentioned above (Γ₋₁₂, T_gas, box size, the actual range of V_cb,i , and self-regulation of ionizing sources) can make Δz_e of several times smaller than in our most extreme scenario (Δz_e = 0.7), it is more realistic to consider a value like Δz_e ∼ 0.1−0.2, which would help to explain observed Δz_e when combined with the other contributions.

The long wavelength of the streaming velocity fluctuations would also leave a unique signature on the reionization observable like the 21 cm background, the near-infrared background, and the CMB via the kinetic Sunyaev–Zel'dovich effect. It is worth exploring these effects in the ionization field using seminumerical (e.g., Mesinger et al. 2011; Sobacchi & Mesinger 2014; Hutter et al. 2020b) or fully numerical simulations of reionization (e.g., CoDa simulation; Ocvirk et al. 2016, 2020).

The authors thank A. Mesinger, C. Cain, A. D'Aloisio, and the anonymous referee for helpful comments on this work. This work was supported by the World Premier International Research Center Initiative (WPI), MEXT, Japan. H.P. was supported by JSPS KAKENHI grant No. 19K23455. K.A. was supported by NRF-2016R1D1A1B04935414 and NRF-2016R1A5A1013277, and appreciates APCTP and IPMU for their hospitality during completion of this work. S.H. was supported by JSPS KAKENHI grant No. 18J01296. P.R.S. was supported in part by U.S. NSF grant No. AST-1009799, NASA grant No. NNX11AE09G, NASA/JPL grant No. RSA Nos. 1492788 and 1515294, and supercomputer resources from NSF XSEDE grant No. TG-AST090005 and the Texas Advanced Computing Center (TACC) at the University of Texas at Austin. Numerical computations were carried out on Cray XC50 and PC cluster at Center for Computational Astrophysics, National Astronomical Observatory of Japan.