COSMIC REIONIZATION ON COMPUTERS. II. REIONIZATION HISTORY AND ITS BACK-REACTION ON EARLY GALAXIES

Nickolay Y. Gnedin; Alexander A. Kaurov

doi:10.1088/0004-637X/793/1/30

1. INTRODUCTION

If cosmic reionization can be called the current frontier of extragalactic astronomy, then, in historic terms, we live in the middle of the 19th century, i.e., the frontier is being settled.

Ultra Deep Field campaigns with the Hubble Space Telescope pushed the search for the most likely reionization sources—young star-forming galaxies—to double-digit values of cosmic redshift (Bouwens et al. 2007, 2011; Oesch et al. 2012; Bradley et al. 2012; Schenker et al. 2013; Willott et al. 2013; Oesch et al. 2013, 2014; Bowler et al. 2014). Observations of Lyα emitters at z ∼ 7 (Hu et al. 2010; Ouchi et al. 2010; Pentericci et al. 2011; Kashikawa et al. 2011; Schenker et al. 2012; Ono et al. 2012; Caruana et al. 2012, 2014) indicate rapid change in their abundance as one rides deeper into frontier territory. Recent mind-blowing progress of the first generation experiments for detecting the redshifted 21cm signal from the epoch of reionization (Parsons et al. 2014; Dillon et al. 2014) promises a major observational breakthrough well before the end of this decade. Even along the well-trodden "Oregon Trail" of Lyα absorption spectroscopy of high-redshift quasars, new advances are expected in the very near future, as new discoveries of z > 6 quasars continue (Bañados et al. 2014; Venemans et al. 2013).

Theoretical studies did not stay behind the observational strides, rejuvenating a somewhat slowed-down progress of the second half of the last decade. The major push on the theory side was galvanized by the pioneering idea of Furlanetto et al. (2004), who realized that the standard lore of large-scale structure theory, Excursion Set formalism, can be applied to studying the reionization process. That idea generated a large following of semi-analytical and semi-numerical approaches for modeling reionization (Furlanetto & Oh 2005; Furlanetto et al. 2006; Mesinger & Furlanetto 2007; Alvarez & Abel 2007; Zahn et al. 2011; Mesinger et al. 2011; Alvarez & Abel 2012; Zhou et al. 2013; Battaglia et al. 2013; Kaurov & Gnedin 2013; Sobacchi & Mesinger 2014), and more traditional models were pursued as well (Choudhury & Ferrara 2005, 2006; Shull & Venkatesan 2008; Mitra et al. 2011; Venkatesan & Benson 2011; Mitra et al. 2012; Kuhlen & Faucher-Giguère 2012; Robertson et al. 2013). Unfortunately, on the numerical simulation front the progress was less dramatic, although important advances in simulation technology did take place (e.g., Iliev et al. 2006; Zahn et al. 2007; McQuinn et al. 2007; Trac et al. 2008; Shin et al. 2008; Croft & Altay 2008; Lee et al. 2008; Iliev et al. 2009; Aubert & Teyssier 2010; Friedrich et al. 2011; Ahn et al. 2012; Shapiro et al. 2012, for a complete review, see Trac & Gnedin 2011). However, the primary brake on the simulation progress—insufficient computing power—is finally being released, thanks to Moore's Law.

Modern High Performance Computing platforms have crossed an important threshold of "sustained peta-scale" performance. This level of performance, currently available on about a dozen or so (nonclassified) supercomputers across the globe, offers a unique opportunity for reionization theorists to make a substantial breakthrough in our ability to model cosmic reionization with high physical fidelity, and some of the most recent simulation work already took advantage of that opportunity (Iliev et al. 2014; So et al. 2014; Norman et al. 2013; Hutter et al. 2014).

The Cosmic Reionization On Computers (CROC) project is another effort in producing peta-scale simulations of reionization in sufficiently large volumes (above 100 Mpc in comoving units), with spatial resolution reaching down to 100 pc, and including most (if not all) of the relevant physical processes, from star formation and feedback to radiative transfer.

In the first paper in this series (Gnedin 2014, hereafter Paper I), we described in complete detail the simulation design and the calibration of numerical parameters. In this paper we explore the overall process of cosmic reionization as captured by CROC simulations and compare our theoretical predictions to several observational constraints.

We deliberately limit the scope of this paper to relatively easily computable quantities, which give only a global, broad-brush view of reionization, due to the limited human effort available for the analysis of the rich, but complex, simulation data. We intend to continue this paper series as more detailed, labor-intensive analyses get completed.

2. CROC SIMULATIONS

All CROC simulations are performed with the Adaptive Refinement Tree (ART) code (Kravtsov 1999; Kravtsov et al. 2002; Rudd et al. 2008). The ART code is capable of modeling a diverse set of physical processes, from the dynamics of dark matter and gas to star formation, stellar feedback, and radiative transfer. A detailed description of all physical processes followed in the CROC simulations is presented in Paper I.

CROC simulations performed in volumes with 20 h⁻¹ Mpc and 40 h⁻¹ Mpc on a side, and 80 h⁻¹ Mpc boxes will be added to the full data set as the project progresses. All simulations but one have the same mass resolution of 7 × 10⁶M_☉ (20 h⁻¹ Mpc boxes use 512³ dark matter particles, 40 h⁻¹ Mpc use 1024³ particles, etc.). One of the 20 h⁻¹ Mpc boxes (B20HR.uv2) has been run with 1024³ particles, achieving the eight times higher mass resolution of 9 × 10⁵M_☉; we use that simulation for testing numerical convergence and for some of the scientific results where high mass resolution is required. Using the full adaptive mesh refinement (AMR) functionality of ART, we reach the formal spatial resolution (smallest simulation cell size) of 125 pc at z = 6, and even higher resolution at earlier times, as our resolution remains constant in comoving units. The real resolution of the simulations is a factor of 2–3 worse.

Paper I describes the numerical parameters of the CROC simulations and how we calibrate their values. The only parameter that we vary in this paper is the "escape fraction of ionizing radiation up to the simulation resolution" epsilon _UV. In a numerical simulation with finite spatial resolution, not all absorptions of ionizing photons can be accounted for because some of the photons will be lost in structures that are not resolved in the simulation (like a parent molecular cloud). Hence, to account for those absorptions, we assign each stellar particle an ionizing luminosity

$L_{\rm ion} = \epsilon _{\rm UV}L^{\rm orig}_{\rm ion},$

where $L^{\rm orig}_{\rm ion}$ is the unattenuated luminosity of a single-age stellar population and the parameter epsilon _UV accounts for unresolved photon losses.

Since the ionizing output of our model galaxies is proportional to epsilon _UV, that parameter critically controls the whole process of reionization in the intergalactic medium (IGM).

For each value of the simulation parameters (box size, epsilon _UV, etc.) we perform a set of simulations that start from independent realizations of initial conditions and properly account for the fluctuations outside the box using the so-called DC mode (Gnedin et al. 2011). Hence, we can use a given simulation set to quantify the effect of cosmic variance on our results.

Table 1 lists the simulation sets that we use in this paper. Star formation and stellar feedback parameters in these simulations are calibrated so that the observed galaxy UV luminosity functions are matched to the observations at all redshifts from z = 6 to z = 10. Hence, for all simulation sets used in this paper stellar sources of cosmic reionization are followed accurately (at least for z ⩽ 10).

Table 1. Simulation Sets

Set Id	_UV	Stopping	Number of
Set Id	_UV	Redshift	Realizations
20 h⁻¹ Mpc boxes, 512³ particles
B20.uv1	0.1	5	6 (A–F)
B20.uv2	0.2	5	6 (A–F)
B20.uv4	0.4	5	3 (D–F)
20 h⁻¹ Mpc boxes, 1024³ particles
B20HR.uv2	0.2	5.7	1 (B)
40 h⁻¹ Mpc boxes, 1024³ particles
B40.uv1	0.1	5	3 (A–C)
B40.uv2	0.2	5.5	3 (A–C)

Download table as: ASCII Typeset image

3. RESULTS

3.1. Reionization History

In Paper I it was showed that the value of epsilon _UV between 0.1 and 0.2 provides the best match to the observed evolution of the Lyα forest at z < 6 (with 0.2 giving a better fit). This is again illustrated in Figure 1, where we show the evolution of the average mass- and volume-weighted hydrogen fractions for the best-fit 20 h⁻¹ Mpc set (B20.uv2) and both 40 h⁻¹ Mpc sets (in all cases we average over all independent realizations).

The shapes of the various curves in the figure are highly expected and are consistent between almost all previous simulations of reionization. Ionized fractions grow steadily with time, as ionized bubbles expand in the neutral IGM. Neutral fractions decrease in response until the moment of overlap, when the volume-weighted neutral hydrogen fraction decreases rapidly (Gnedin 2000). In the post-overlap stage the IGM is highly ionized, and the evolution of neutral fractions is governed by the mean free path of ionizing photons and the level of cosmic ionizing background (see also Trac & Gnedin 2011, for a general overview of reionization process).

The lack of numerical convergence, which we discussed in Paper I, is also visible in Figure 1—in the set B40.uv1 the overlap of ionized bubbles (indicated by the rapid drop in the average volume-weighted ${\scriptsize {HI}}$ fraction just before z = 6) occurs at about the same time as in the smaller box set B20.uv2, but the post-reionization Lyα forest is better matched by the set B40.uv2, which has a significantly earlier overlap.

As we discussed in Paper I, that lack of convergence is caused by cosmic variance⁴—having multiple independent realizations allows us to explore it well. Our (still unconverged) simulations sample a much larger number of independent sightlines than is available observationally, hence the lack of agreement between the simulations and observations at z > 6 cannot yet be taken seriously.

At z < 6 the situation is, however, completely different: simulations with the same value of the epsilon _UV parameter do converge, the cosmic variance is small (since the radiation field is dominated by the cosmic background—this is apparent from Figure 4 of Paper I), and the correct simulations should match the observational data. The mismatch between the simulations and the observations at z ≲ 5.3 is, therefore, real. In Paper I we also showed that, similarly, our simulations fail to match the galaxy UV luminosity function at z = 5 well enough. Both these discrepancies indicate that our simulations become inaccurate after z ≈ 5.3, most likely because our spatial resolution, being kept fixed in comoving units, degrades too much by z ≈ 5.

Another illustration of the role of cosmic variance is shown in Figure 2, where we plot the distribution of the average Gunn–Peterson optical depth 〈τ_GP〉_{Δz = 0.15} as a function of redshift for individual lines of sight. Lines of sight are generated randomly (i.e., starting at random locations and going in random directions) throughout each of the simulation boxes, and synthetic Lyα spectra are generated along each of them. The optical depth along each line of sight is computed as an average over a distance of 40 h⁻¹ Mpc (Δz ≈ 0.15), which is well-matched to the redshift interval used in observational studies (Fan et al. 2006). At z > 6 the distribution becomes exceptionally wide, and not just in the tails. At lower redshift, past overlap, the distribution narrows significantly, but still maintains a relatively long tail toward large values of 〈τ_GP〉_{Δz = 0.15}. Hence, we predict that, in a small fraction of all quasar sightlines, segments with 〈τ_GP〉_{Δz = 0.15} as high as 10 can be observed all the way down to z ≈ 5.5.

Figure 3 shows the same result in a more conventional way, where we plot a cumulative probability distribution for 〈τ_GP〉_{Δz = 0.15} in three redshift bins after reionization (such a distribution is, obviously, dependent on the redshift interval over which averaging is done). Even at z ∼ 5.5 there is a 0.4% probability to find a line of sight with no detectable Lyα flux. For such a large redshift interval to be devoid of any flux, it is not enough to just have a damped Lyα in that line of sight, rather it is a genuine feature of the cosmic variance due to large-scale density fluctuations.

3.2. Ionized Bubbles

Time evolution of the distribution of ionized bubbles is one of the most important characteristics of the reionization process. Unfortunately, in a realistic cosmological simulation the concept of an "ionized bubble" is not well-defined mathematically, especially at late times, as can be easily seen from Figure 4. Hence, in order to have a working definition that can also be compared with other studies, we closely follow the procedure from Zahn et al. (2007) to compute the probability that a given point in the simulation is located inside an ionized bubble of size R. The scale R is defined as the largest radius of a sphere in which the volume-weighted average ionized fraction is higher than the threshold value, which is chosen to be 90%. The only difference with Zahn et al. (2007) is normalization: we normalize the bubble size distribution to the total volume (an integral over the distribution is equal to the volume-weighted average ionized fraction), whereas Zahn et al. (2007) normalize their distributions to the total ionized volume (an integral over the distribution is equal to unity).

**Figure 4.** Slices through the computational domain for the first realization of the B40.uv1 set (run B40.uv1.A in the notation of Paper I) at z = 8 and z = 7. While at z = 8 individual ionized bubbles can still be identified reasonably well, by z = 7 the concept of a "bubble" becomes ill-defined.
Download figure:
Standard image High-resolution image

Distributions of sizes of ionized and neutral bubbles (shown as a differential volume function) are presented in Figure 5 for several values of redshift for both 20 h⁻¹ Mpc and 40 h⁻¹ Mpc simulation sets. Somewhat unexpectedly, we find good convergence in the bubble size distribution even for the 20 h⁻¹ Mpc box size all the way to z ∼ 7. At lower redshifts the convergence does break down, simply because some of the smaller boxes are going to be completely ionized before larger boxes (and the same applies to neutral "bubbles" at high redshifts). The volume fraction in such boxes is shown with horizontal arrows on both panels, and it also is reasonably consistent between the 20 h⁻¹ Mpc and 40 h⁻¹ Mpc simulation sets.

At face value, this result is inconsistent with the recent simulations of Iliev et al. (2014), who found incomplete convergence in simulation volumes as large as 114 h⁻¹ Mpc at all redshifts. Without detailed comparison, it is difficult to isolate the source of the discrepancy. We notice, however, that the spatial resolution of the radiative transfer solver in Iliev et al. (2014) simulations is only about 200 h⁻¹ kpc in comoving units, which is inadequate for resolving absorptions in galactic halos and Lyman limit systems, while our spatial resolution (0.6 h⁻¹ comoving kpc) is better suited for properly accounting for all absorptions of ionizing radiation.

3.3. Damping Wing of Lyα Absorption

Before the universe is completely reionized, patches of the still-neutral IGM can significantly absorb Lyα emission from high-redshift galaxies, as the damping wing of Lyα absorption extends far redward of the galaxy systemic velocity (Miralda-Escude 1998). In order to model that effect, we also generate synthetic Lyα spectra that originate at galaxy locations. The "sky" of each galaxy is sampled uniformly with 12 directions, corresponding to the 12 zero level cells of the HEALPix⁵ tessellation of a sphere (Górski et al. 2005). In order to exclude the local absorption from the galactic ISM, we start the line of sight 10 kpc away from the center of the galaxy.

An exact calculation of the effect of the damping wing on the Lyα emission line of a galaxy requires complex Lyα radiative transfer in the galactic ISM and surrounding IGM. Such a calculation is a separate research project in itself, and in any case our simulations do not have enough spatial resolution to perform such a calculation with sufficient accuracy. Instead, we approximate the effect of the damping wing by computing the absorption equivalent width of the red part of the synthetic spectrum,

$D_{\rm EW}= \int _{\lambda _0}^\infty e^{\displaystyle -\tau _\lambda } d\lambda,$

where λ₀ is the wavelength of Lyα at the systemic velocity of each model galaxy. Hence, we compute 12 values of D_EW for each galaxy, achieving dense sampling of the full distribution function for D_EW.

Comparing Figures 1 and 6, it is clear that the beginning of the overlap of ionized bubbles corresponds to a rapid decrease in the equivalent width of the damping wing—in our fiducial B20.uv2 run D_EW drops by an order of magnitude between z = 7 and z = 6.5. Such behavior is consistent with the observed sharp decline in the fraction of Lyα emitters between z = 6 and z = 7 (Pentericci et al. 2011; Schenker et al. 2012; Caruana et al. 2014), although, as we mentioned above, a more quantitative comparison would require a better model of Lyα emitters in the simulations.

**Figure 6.** Equivalent width of Lyα absorption D_EW as a function of redshift (top) or volume-weighted neutral hydrogen fraction (bottom). Colored lines show averages over all realizations of simulation sets B20.uv4 (red), B20.uv2 (green), and B20.uv1 (blue) for all galaxies with UV magnitudes between −22 and −18. Green semi-transparent bands give the distribution of D_EW around the mean (as 1–99 and 10–90 percentiles). Green dashed and dotted lines show D_EW for the subsets of galaxies with magnitudes in bins [−21.72, −20.25] and [−20.25, −18.75], respectively.
Download figure:
Standard image High-resolution image

We also notice that in our simulations there is little dependence of the damping wing equivalent width D_EW on galaxy luminosity—dotted and dashed lines in Figure 6 show the evolution of D_EW for two luminosity bins, but both lines trace similar behavior. The luminosity dependence of the fraction of Lyα emitters has been seen in some observational studies (cf. Schenker et al. 2012) but not in others (cf. Pentericci et al. 2011). If such dependence is confirmed by further observations, it would imply that brighter galaxies have an intrinsically higher probability of becoming a Lyα emitter than fainter ones.

The observed measurements of the fraction of galaxies that remain strong Lyα emitters have also been used to constrain the mean volume-weighted neutral fraction. The bottom panel of Figure 6 replaces the redshift axis with the (monotonically decreasing with time) volume-weighted ${\scriptsize {HI}}$ fraction. The sensitivity of $D_{\rm EW}(\langle X_{{\scriptsize {HI}}}\rangle _V)$ to the epsilon _UV parameter is much less than when D_EW is treated as a function of z, confirming the validity of the assumption that the decrease in the observed fraction of Lyα emitters at higher redshifts indicates a change of the volume-weighted average neutral fraction. In fact, it appears that the condition D_EW < 1000 km s⁻¹ corresponds to $\langle X_{{\scriptsize {HI}}}\rangle _V\lesssim 0.2$ , while the condition D_EW < 500 km s⁻¹ corresponds to $\langle X_{{\scriptsize {HI}}}\rangle _V\lesssim 0.1$ . This conclusion is in good agreement with other recent simulation studies (c.f. Taylor & Lidz 2014; Hutter et al. 2014).

3.4. Back-reaction of Reionization on Early Galaxies

The most immediate effect of reionization on the early galaxies—or, rather, galactic halos—is to expel photoionized gas from sufficiently low mass halos. This process, sometimes inaccurately called "photoevaporation," has been a focus of a large number of studies, reviewing which is beyond the scope of this paper. So far, the highest mass resolution (up to 3 × 10⁴M_☉—a critical numerical parameter for this question) has been achieved by Okamoto et al. (2008, they also review the previous works), who provided accurate fits for the average gas fraction as a function of halo mass and redshift. In Figure 7 we compare the gas fractions in our simulations with the fits from Okamoto et al. (2008). Our fiducial simulation sets barely have enough mass resolution to properly capture the characteristic mass below which halos start losing gas due to photoionization. The high resolution run B20HR.uv2 captures that effect well, and our results post-reionization (at z ≲ 6) agree with the Okamoto et al. (2008) fits extremely well. At earlier redshifts the difference is expected, as Okamoto et al. (2008) assumed an instantaneous reionization at z = 9, while we model the actual reionization history.

**Figure 7.** Average gas fractions (in units of the universal) for the high resolution run B20HR.uv2 at z = 6, 7, and 8 (solid lines) and for the fiducial B20.uv2 run at z = 5 (at earlier times the fiducial run lack enough mass resolution). A light orange band shows the rms scatter around the z = 6 line for the B20HR.uv2 run. Thin black lines are fits from Okamoto et al. (2008) that agree with our simulations after reionization (z ≲ 6) almost perfectly.
Download figure:
Standard image High-resolution image

A good agreement of our results with Okamoto et al. (2008) illustrates a simple but not widely appreciated fact that it takes only a few hundred million years for the effect of reionization on the gas fractions to be fully established (Iliev et al. 2005)—for example, our B20HR.uv2 run reionizes at z ≈ 7.3, and fully converges with the Okamoto et al. (2008) simulations by z = 6, only 210 Myr later.

A secondary effect of reionization is on the actual star formation rates in early galaxies. That effect can be large or nil, depending on the fraction of star formation in halos that are affected by photoionization. The simplest manifestation of such back-reaction is a change in the global cosmic star formation history. Studies of the response of the global star formation history to reionization were pioneered by Barkana & Loeb (2000), who found a large suppression in the global star formation history at reionization. This "Barkana & Loeb" effect has been revisited repeatedly in previous studies (Tassis et al. 2003; Wyithe & Loeb 2006; Barkana & Loeb 2006; Davé et al. 2006; Wyithe & Cen 2007; Pieri & Martel 2007; Yoshida et al. 2007; Muñoz & Loeb 2011; Duffy et al. 2014), with different groups disagreeing significantly regarding its strength. Global star formation histories of our three fiducial simulation sets as well as the high resolution run B20HR.uv2 are shown in Figure 10. No effect of reionization is noticeable in the figure, in contradiction with some of the previous studies.

In order to elucidate that disagreement, we show in Figure 9 the cumulative fraction of molecular gas in halos above a given mass at several redshifts. There is always at least a decade in mass difference between the halos that are affected by reionization and halos that contain a substantial amount of molecular gas. Since our model for star formation is crucially based on the (observationally motivated) paradigm that stars form primarily in the molecular gas, the negligible back-reaction on the star formation in early galaxies is naturally explained by the large difference in the two mass scales.

As a side note, we recall that our fiducial runs and the high resolution run B20HR.uv2 both reproduce observed galaxy UV luminosity functions at all redshifts z ≳ 6 (Figure 8 of Paper I), but they have global star formation histories that differ at z ≈ 10 (Figure 8) by more than the claimed observational errorbars from Bouwens et al. (2014). Hence, we conclude that the observational determination of the global star formation history is based on assumptions that do not always hold, and, hence, the systematic errors of such determinations are substantially larger than the formal statistical errors at z > 8.

**Figure 8.** Cosmic star formation histories for the three fiducial simulation sets with varied ionizing intensity and the high resolution run B20HR.uv2. Thin vertical lines mark reionization redshifts for each simulation set (defined as times when the volume-weighted neutral fraction falls below 10⁻³; see Gnedin 2000). Black points with error bars are the data from Bouwens et al. (2014). Reionization does not introduce any noticeable feature in the global star formation history.
Download figure:
Standard image High-resolution image

**Figure 9.** Average gas fractions (dotted lines—the same as shown in Figure 7) and cumulative molecular gas fractions in halos above a certain mass (solid lines) for the high resolution run B20HR.uv2 at several redshifts. Notice that halos affected by reionization contain little molecular gas and, hence, form no stars.
Download figure:
Standard image High-resolution image

While not affecting the global star formation history in any significant way, reionization may still leave more subtle signatures in the properties of early galaxies. For example, the sensitivity of the faint end slope of the galaxy luminosity function to the reionization history has been proposed as an important science goal for the James Webb Space Telescope (JWST; Gardner et al. 2006a). To explore the feasibility of such a test, we show in Figure 10 the galaxy luminosity functions for the three fiducial simulation sets with varying ionizing emissivity. In this case the effect of reionization is detectable, although it is not large. As the bottom panel of Figure 10 shows, a factor of two variation in the ionizing emissivity (which corresponds to about a Δz ≈ 0.5 change in the redshift of reionization—see Figure 6 and Figure 4 of Paper I) from our fiducial value epsilon _UV = 0.2 produces a change of about 0.05 in the faint end slope of the luminosity function for M₁₅₀₀ ≳ −17, although at brighter magnitudes the difference rapidly disappears. This variation is likely to be too small to be usable as a constraint on the reionization history.

**Figure 10.** Top: ultraviolet galaxy luminosity functions for the three simulation sets with varied ionizing intensity at six different redshifts (as shown in the legend). Circles with error bars are a compilation of recent observational measurements (Bouwens et al. 2007, 2011; Oesch et al. 2012; Bradley et al. 2012; Schenker et al. 2013; Willott et al. 2013; Oesch et al. 2013, 2014; Bowler et al. 2014). Different redshifts are shifted vertically by 1 dex for clarity. Bottom: differences between simulation sets with _UV = 0.1 (dotted lines) and _UV = 0.4 (solid lines) relative to the reference value of _UV = 0.2. Black dashed lines show the variation in the faint-end slope by ±0.05.
Download figure:
Standard image High-resolution image

epsilon — **Figure 10.** Top: ultraviolet galaxy luminosity functions for the three simulation sets with varied ionizing intensity at six different redshifts (as shown in the legend). Circles with error bars are a compilation of recent observational measurements (Bouwens et al. 2007, 2011; Oesch et al. 2012; Bradley et al. 2012; Schenker et al. 2013; Willott et al. 2013; Oesch et al. 2013, 2014; Bowler et al. 2014). Different redshifts are shifted vertically by 1 dex for clarity. Bottom: differences between simulation sets with _UV = 0.1 (dotted lines) and _UV = 0.4 (solid lines) relative to the reference value of _UV = 0.2. Black dashed lines show the variation in the faint-end slope by ±0.05.
Download figure:
Standard image High-resolution image

4. CONCLUSIONS

We present the reionization history and global characteristics of the reionization process from a suite of recent numerical simulations performed as part of the Cosmic Reionization On Computers (CROC) project. CROC simulations reproduce the observed evolution of the galaxy UV luminosity function between z = 10 and z = 6 well, and, hence, include realistic treatment of the dominant class of ionizing sources.

We find that, in order to match the observational constraints on the post-reionization Lyα forest at 5 < z < 6, we need to set the ionizing emissivity parameter epsilon _UV (which measures the escape fraction up to the resolution limit of our simulations) to just under epsilon _UV = 0.2. However, as we also emphasized in Paper I, cosmic variance increases sharply with redshift, and at z > 6 our simulations do not yet converge on the global properties of the IGM, such as the mean Gunn–Peterson optical depth or the volume-weighted ${\scriptsize {HI}}$ fraction. Since the statistical power of our simulations is much higher than the statistical reach of the existing absorption spectra of high-redshift quasars, we conclude that, unfortunately, the observations are unlikely to have reached convergence either.

In a further illustration of this, we show that the distribution of the Gunn–Peterson optical depth over the redshift intervals Δz ≈ 0.15 is extraordinarily wide at z > 6, but even at z < 6 the τ_GP distribution retains a relatively long tail toward high values.

The distributions of ionized and neutral bubbles during most of cosmic reionization is approximately flat, meaning that it is roughly equally likely for a random place of the universe to be in a large or a small bubble. We find good numerical convergence in bubble sizes down to z ∼ 7, at which point the finite sizes of our simulation boxes start biasing the distribution of ionized bubbles. That result illustrates the importance of achieving consistent numerical resolution between the gas dynamic solver and the radiative transfer solver—the mismatch between the two resolution likely results in erroneous overpropagation of ionizing radiation beyond the few mean free path lengths.

We show that the equivalent width of the damping wing of Lyα absorption increases rapidly from mellow values of D_EW ∼ 100 km s⁻¹ at z = 6 to a whopping D_EW ∼ 2000 km s⁻¹ by z = 7. While D_EW serves only as a rough proxy for the suppression of the galaxy Lyα emission line by the neutral IGM in front of it, this result is generally consistent with the observed sharp decline in the fraction of Lyα emitting galaxies at z = 7 as compared to z = 6. We also confirm conclusions from the previous simulation and analytical work that such suppression corresponds to a substantial, but not dominant, volume-weighted neutral fraction of about 0.2.

While our results on the reionization history are in good agreement with most previous studies, we find little back-reaction of reionization on the properties of early galaxies. Because galaxies that are affected by photoionization contain little molecular gas (and, hence, star formation), we find that the global star formation history is insensitive to the reionization history, i.e., the "Barkana & Loeb" effect does not exist. A more subtle effect of reionization is in modifying the faint end slope of the galaxy UV luminosity function, but such a modification is rather small (change in the slope of about 0.1 for a unit shift in the redshift of reionization). Since predicting the faint end slope to such precision theoretically would be extremely challenging, we conclude that, unfortunately, measuring the faint end slope by JWST will not be a useful constraint on reionization, contrary to expectations (Gardner et al. 2006b).

One observational constraint that we have ignored so far is the optical depth to Thompson scattering from the CMB observations by the WMAP mission. While the history of WMAP measurements of the Thompson optical depth is rocky, the latest value from the nine-year WMAP data is 0.089 ± 0.014 (or 0.081 ± 0.012 if other data are included in a joint fit, Bennett et al. 2013; Hinshaw et al. 2013). The value we get for fiducial sets B20.uv2 and B40.uv1 is 0.052 ± 0.003, and in the set B40.uv2, which reionizes earlier, the value for the Thompson optical depth only rises to 0.057 ± 0.004, which is only marginally (at a 2σ level) consistent with the WMAP values.

A large portion, if not all, of this discrepancy is due to incomplete numerical convergence of our simulations. In Paper I we compared our fiducial runs (an equivalent of 512³ particles in a 20 h⁻¹ Mpc box) with a single higher mass resolution run B20HR.uv2 that we were able to complete (an equivalent of 1024³ particles in a 20 h⁻¹ Mpc box). While numerical converge tests indicate that our fiducial runs account for 55% of all ionizing photons, the higher reslution B20HR.uv2 run accounts for 80% of them. As a result, the Thompson optical depth rises to 0.067 in that run. Simple linear extrapolation to the limit of 100% of ionizing radiation gives a value of 0.08 for the Thompson optical depth, fully consistent with the current observational measurements.

Determining whether incomplete numerical convergence is, indeed, the full story will have to wait for more powerful computers, however, as at present we are unable to run the whole ensemble of higher mass resolution simulations—for example, a higher mass resolution equivalent of our planned 80 h⁻¹ Mpc run would have 4096³ particles and will require of order of 200 million CPU hours, an amount not currently feasible to obtain for this kind of work.

We are grateful to George Becker for valuable comments on the early draft of this paper.

Simulations used in this work have been performed on the Joint Fermilab-KICP cluster "Fulla" at Fermilab, on the University of Chicago Research Computing Center cluster "Midway," and on National Energy Research Supercomputing Center (NERSC) supercomputers "Hopper" and "Edison."

COSMIC REIONIZATION ON COMPUTERS. II. REIONIZATION HISTORY AND ITS BACK-REACTION ON EARLY GALAXIES

Article metrics

Permissions

Author e-mails

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION

2. CROC SIMULATIONS