Semi-numeric simulations of helium reionization and the fluctuating radiation background

Abstract

INTRODUCTION

Throughout most of cosmic history, the metagalactic ionizing background evolved smoothly. The reionization of hydrogen (z ≳ 6) and of helium (fully ionized at z ∼ 3) were two major exceptions. During these phase transitions, the intergalactic medium (IGM) became mostly transparent to the relevant ionizing photons, allowing the high-energy radiation field to grow rapidly. Aside from its fundamental importance as a landmark event in cosmic history, the full reionization of helium has several important consequences for the IGM. Since helium makes up ∼24 per cent of the baryonic matter by mass, its ionization state significantly affects the mean free path of photons above 54.4 eV, the energy required to fully ionize helium. Additionally, reionizing helium should significantly alter the thermal history of the IGM (e.g. Hui & Gnedin 1997; Schaye et al. 1999, 2000; Zaldarriaga 2002; Furlanetto & Oh 2008b; McQuinn et al. 2009, hereafter M09; Garzilli et al. 2012). The resultant heating from reionization may also affect galaxy formation by increasing the pressure (and hence Jeans mass) of the gas, which may affect the star formation histories of galaxies (Wyithe & Loeb 2007). Determining the timing and duration of this epoch is, therefore, important for understanding the evolution of the IGM.

Because of its large ionization potential, He ii reionization is driven by high-energy sources: quasars and other active galactic nuclei. This process, therefore, offers an indirect probe of these populations and their distribution through space and even inside dark matter haloes. For example, the metagalactic background above 54.4 eV can even be used to constrain detailed properties of the sources, such as their lifetimes and emission geometries (Furlanetto & Lidz 2011).

The study of helium reionization has some advantages over its hydrogen counterpart (see Barkana & Loeb 2001 for a review) and may serve as an object lesson for our understanding of that earlier era. Unlike for the reionization of hydrogen, we have excellent data on the state of the IGM at z ∼ 3, including the temperature and density distribution (see Meiksin 2009 for a review). Furthermore, the main drivers of helium reionization are quasars (due to their hard spectra), which are much better understood than their counterparts that drive hydrogen reionization. In particular, we know the quasar luminosity function (QLF; e.g. Wolf et al. 2003; Richards et al. 2006; Masters et al. 2012) and quasar clustering properties (e.g. Ross et al. 2009; Shen et al. 2009; White et al. 2012).

Given the rarity of these sources and the inhomogeneity of the IGM, the process of helium reionization is expected to be patchy. This patchiness has important implications for the evolution of the IGM, and it is one of the crucial elements that must be modelled correctly. Many measurements make assumptions about the UV-background radiation and the temperature–density relation, ignoring the inhomogeneity and timing of He ii reionization. Therefore, quantifying the expected spread in He ii ionization states is crucial.

Recent and ongoing high-resolution measurements of the He ii Lyman α (Lyα) forest have started to test this picture. Dixon & Furlanetto (2009) showed that early data indicated a rapid decrease in the optical depth with cosmic time at z ≳ 2.7, which may be indicative of He ii reionization, while Furlanetto & Dixon (2010) showed that the large fluctuations observed in the optical depth were most likely due to ongoing reionization at z ∼ 2.8. Now, the Cosmic Origins Spectrograph on the Hubble Space Telescope is providing an even more powerful probe of the IGM at 2.3 ≲ z ≲ 3.7. The He ii optical depth along the explored sightlines varies significantly at z > 2.7 and generally increases rapidly at higher redshift (Shull et al. 2010; Worseck et al. 2011; Syphers et al. 2012; Syphers & Shull 2013). Four lines of sight also show significant variation in the He ii Lyβ optical depth (Syphers et al. 2011b). This large set of observations – and our detailed knowledge of the source populations and IGM at z ∼ 3 – make a detailed study of the He ii reionization process very timely.

Since quasars are rare, the ionizing background should exhibit significant fluctuations even after reionization is complete (Fardal, Giroux & Shull 1998; Bolton et al. 2006; Furlanetto 2009; Meiksin 2009). As mentioned above, the He ii Lyα forest measurements provide some direct evidence for this behaviour (Furlanetto & Dixon 2010). Additionally, the clumpy IGM (and radiative transfer through it) can induce additional fluctuations (Schaye 2001; Maselli & Ferrara 2005; Tittley & Meiksin 2007). During reionization, fluctuations in the UV radiation background are even greater, because some regions receive strong ionizing radiation while others remain singly ionized with no local illumination.

During the past decade, the reionization of helium and its impact on the IGM has received increased theoretical attention through a variety of methods. Both large-scale structures (spanning ∼1 Gpc comoving scales), given the rarity of quasars and their clustering properties, and the small scales of gas physics play a significant role in fully understanding this epoch, but accounting for both simultaneously is computationally challenging. Since this epoch is observationally constrained in a way that hydrogen reionization is not, a viable model must conform to these observations, such as the IGM density and temperature.

Several studies employ various combinations of analytic, semi-analytic, and Monte Carlo methods (e.g. Theuns, Schaye & Haehnelt 2000; Theuns et al. 2002; Gleser et al. 2005; Furlanetto & Oh 2008a; Furlanetto 2009) to approximate the morphology of ionized helium bubbles, heating of the IGM, He-ionizing background, etc. None of these methods can comprehensively include all relevant physics, especially spatial information like source clustering. Many numerical simulations of this epoch focus on scales ≤100³ comoving Mpc³ (Sokasian, Abel & Hernquist 2002; Paschos et al. 2007; Meiksin & Tittley 2012). These simulations miss the large ionized bubbles expected during helium reionization and fail to include many sources. M09 present large (hundreds of Mpc) N-body simulations with cosmological radiative transfer as a post-processing step. Note that M09 consider two classes of quasar models and vary the timing of reionization independently.

One major modelling hurdle is the implementation of quasars. The aggregate, empirical properties of quasars are well constrained with the QLF measured over a range of redshifts and to low optical luminosities (see Hopkins, Richards & Hernquist 2007). The spectral energy distribution is more uncertain and varies significantly from quasar to quasar (e.g. Telfer et al. 2002; Scott et al. 2004; Shull, Stevans & Danforth 2012). Furthermore, the lifetime of quasars (e.g. Kirkman & Tytler 2008; Kelly et al. 2010; Furlanetto & Lidz 2011) and their relation to dark matter haloes (e.g. Hopkins et al. 2006; Wyithe & Loeb 2007; Conroy & White 2013) are far from settled. The challenge is to satisfy the observed quantities while allowing for a range of possibilities for the more uncertain aspects.

To complement these numeric and analytic studies, here we present fast, semi-numeric methods that incorporate realistic source geometries to explore a large parameter space. We adapt the established hydrogen reionization code dexm (Mesinger & Furlanetto 2007) to the z = 3 era, as appropriate for helium. The code applies approximate but efficient methods to produce dark matter halo distributions. Using analytic arguments, ionization maps are generated from these distributions. The advantages of this approach are speed, as compared to cosmological simulations, and reasonably accurate spatial information (like halo clustering and a detailed, local density field), as compared to more analytic studies. We investigate a wide range of quasar properties.

These simulations are well suited to study many features of the helium reionization epoch, including interpreting the He ii Lyα forest and estimating the He-ionizing background. In this paper, we consider the spatial morphology of He ii reionization and the He ii photoionization rate during and after reionization. Given our ability to span a large parameter space, we aim to pinpoint the assumptions that most strongly impact the results and to evaluate the effect of various uncertainties.

We use semi-numeric methods, outlined in Section 2, to approximate the ionization morphology, density field, and sources of ionizing photons at z = 3. As a second layer, we place empirically determined active quasars in dark matter haloes, following several prescriptions in Section 3. With these inputs, we determine the He ii photoionization rate for post-reionization and several He ii fractions in Section 4. In Section 5, we compare our results and methods with those for hydrogen reionization. We conclude in Section 6.

In our calculations, we assume a cosmology with Ω_m = 0.26, Ω_Λ = 0.74, Ω_b = 0.044, H₀ = 100 h km s⁻¹ Mpc⁻¹) (with h = 0.74), n = 0.95, and σ₈ = 0.8. Unless otherwise noted all distances quoted are in comoving units.

SEMI-NUMERIC SIMULATIONS OF He ii REIONIZATION

To fully describe helium reionization and compute the detailed features of the Lyα forest, complex hydrodynamical simulations of the IGM, including radiative transfer effects and an inhomogeneous background, are required. Recent N-body simulations (e.g. M09) have advanced both in scale and in the inclusion of relevant physical processes. However, the major limitation is the enormous dynamic range, balancing large ionized regions with small-scale physics, required to model reionization. Generally, this problem is solved in two ways. First, most simulations follow only the dark matter, assuming that the baryonic component traces it according to some simplified prescription. Radiative transfer algorithms are applied to the dark matter field (perhaps with some modifications to reflect baryonic smoothing). Secondly, these codes add in the ionizing sources through post-processing, rather than through a self-consistent treatment of galaxy and black hole formation, which is not substantially better than simple analytic models. Also importantly, large simulations remain computationally intensive, limiting the range of parameter space that can be explored in a reasonable amount of time.

We adapt the semi-numeric code dexm¹ (Mesinger & Furlanetto 2007), originally designed to model hydrogen reionization at z ≳ 6, for lower redshifts and for helium reionization. This code provides a relatively fast semi-numeric approximation to more complicated treatments but remains fairly accurate on moderate and large scales. For full details, see Mesinger & Furlanetto (2007). In brief, we create a box with length L = 250 Mpc on each side, populated with dark matter haloes determined by the density field. Given these haloes, we generate a two-phase (singly and fully ionized) ionization field, where the He iii fraction is fixed by some efficiency parameter using a photon-counting method. For our high-level calculations, we additionally adjust the local density field to match recent hydrodynamical simulations.

There are drawbacks to our approach, of course. For one, this semi-numeric method does not follow the progression of He ii reionization through a single realization of that process. In order to ensure accurate redshift evolution, photon conservation requires a sharp k-space filter, whereas we use a real-space, top-hat filter, applied to the linear density field (see the appendix of Zahn et al. 2007, 2011). This filtering can result in ‘ringing’ artefacts from the wings of the filter response, in the limit of a rare source population (as is the case in He ii reionization). Moreover, during He ii reionization, the relevant ionizing sources are short-lived quasars. When a source shuts off, the surrounding regions can begin to recombine in a complex fashion (Furlanetto, Haiman & Oh 2008). Our method does not explicitly follow these recombinations, as we will discuss below. As a consequence, we perform all of our calculations at z = 3, near the expected mid-point of the He ii reionization according to current measurements. We expect our results to change only modestly throughout the redshift range z ∼ 2.5–3.5, because the source density does not change significantly over this range and the IGM density field evolution is similarly limited.

Our primary concern in adapting dexm to this period, as compared to the higher redshift regime for which it has been carefully tested, is the increasing non-linearity of the density field. There are two possible problems. First, the dramatically larger halo abundance may lead to overlap issues with the halo-finding algorithm. Secondly, the IGM density field will have modest non-linearities representing the cosmic web. We will address the latter problem in detail in Section 2.4.

Dark matter haloes

The backbone of our method is the dark matter halo distribution. In a similar way to N-body codes, we generate a 3D Monte Carlo realization of a linear density field on a box with N = 2000³ grid cells (L = 250 Mpc). This density field is evolved using the standard Zel'dovich approximation (Zel'Dovich 1970). Dark matter haloes are identified from the evolved density field using the standard excursion set procedure. Starting on the largest scale (the box size) and then decreasing, the density field is smoothed using a real-space top-hat filter at each point. Once the smoothed density exceeds some mass-dependent critical value, the point is associated with a halo of the corresponding mass. Each point is assigned to the most massive halo possible, and haloes cannot overlap. We choose our critical overdensity to match the halo mass function found in the Jenkins et al. (2001) simulations at z = 3, following the fitting procedure of Sheth & Tormen (1999). The final halo locations are shifted from their initial (Lagrangian) positions to their evolved (Eulerian) positions using first-order perturbation theory (Zel'Dovich 1970). To optimize memory use, we resolve the velocity field used in this approximation on to a lower resolution, 500³ grid.

The left-hand panel of Fig. 1 displays the resultant mass function as points, where the dashed line is the Sheth–Tormen (Sheth & Tormen 1999) result with Jenkins et al. (2001) parameters and the dotted line is the Press–Schechter analytic solution for the mass function (Press & Schechter 1974). We successfully reproduce the former mass function and, therefore, match large numerical simulations. The right-hand panel shows our halo power spectrum, defined as |$\Delta ^2_{\rm hh}(k,z) = k^3/(2\pi ^2V)\langle |\delta _{\rm hh} ({{\boldsymbol k}},z)|^2\rangle _k$|⁠. Here, |$\delta _{\rm hh} ({{\boldsymbol x}},z) \equiv M_{\rm coll}({{\boldsymbol x}},z)/\langle M_{\rm coll}(z)\rangle - 1$| is the collapsed mass field, and the volume V = L³. Also shown is the prediction from the halo model, which is a semi-analytic method for describing non-linear gravitational clustering. This model matches many statistical properties of large simulations (see Cooray & Sheth 2002). Specifically, all dark matter is assumed to be bound in haloes of varying sizes; then, given prescriptions for the distribution of dark matter within haloes and the number and spatial distribution of haloes, many large-scale statistical properties of dark matter can be calculated, including the halo power spectrum. At small scales, the measured power spectrum diverges from the halo-model equivalent, most likely because we do not include halo density profiles that would dominate at these scales.

Figure 1.

Left-hand panel: the calculated halo mass function shown as points. The dotted curve is the Press–Schechter analytic mass function. The dashed curve shows the Sheth–Tormen mass function, which agrees with our results by construction. Right-hand panel: halo power spectra. The solid line represents the halo power spectrum measured from our standard (L = 250 Mpc) box, and the dashed line is the halo model analytic equivalent.

Open in new tabDownload slide

Note that the haloes with which we are concerned are the most massive at z = 3, containing only a few per cent of the total mass. This fact means that overlap between mass filters is not very important, as it might be at lower masses: when a large fraction of the mass is contained inside collapsed objects, the halo filters could then overlap, in which case the ordering at which the cells were examined would affect the halo distribution and bias our results (Bond & Myers 1996). Fortunately, because we are concerned with quasars inside very massive dark matter haloes, we avoid this issue. This focus on large haloes also accounts for the large difference between the Sheth–Tormen and Press–Schechter mass functions in Fig. 1.

Ionizing sources

Because we focus on helium reionization, quasars are the relevant ionizing source.² Empirically, quasar properties are reasonably well constrained: the QLF has been measured to low optical luminosities over a range of redshifts (see the compilation in Hopkins et al. 2007). However, the spectral energy distribution – and hence the transformation of this observed luminosity function to the energies of interest to us – is more uncertain. Furthermore, the relation between quasars and their host dark matter haloes is not settled (e.g. Hopkins et al. 2006; Wyithe & Loeb 2007; Conroy & White 2013). These details affect the ionization maps, the photoionization rate, and the pertinent range of the halo mass function via the minimum mass. Given these factors, there is no standard approach to the inclusion of quasars in numeric and analytic models. To address this uncertain relationship, we considered a wide range of models that encompass the spectrum of popular theories.

The efficiency parameter ζ can be estimated from first principles given a model for quasars (see, e.g. Furlanetto & Oh 2008a). One could compute the average number of He ii ionizations per baryon inside of quasars (N_ion), estimate the efficiency with which baryons are incorporated into the source black holes (f_BH), and estimate the fraction of these photons that escape their host (f_esc); then ζ ∼ f_BHf_escN_ion, with additional corrections for the composition of the IGM and the number of photons ‘wasted’ through recombinations. However, we note that ζ cannot be inferred directly from the observed luminosity function, because it depends on the total number of ionizing photons emitted by a halo rather than the current luminosity. At minimum, we would require knowledge of the lifetime of the luminous phase of the quasar emission. Conversely, this implies that our choices for M_min and f_host enforce a quasar lifetime, which we will consider in Section 3.

This ζ also depends on the recombination history of the IGM, as (to the extent that they are uniform) such recombinations essentially cancel out some of the earlier ionizations (Furlanetto, Zaldarriaga & Hernquist 2004). In practice, recombinations are inhomogeneous and so much more complex to incorporate (Furlanetto & Oh 2005), especially because the small number of quasar sources can allow some regions to recombine substantially before being ionized again by a new quasar (Furlanetto et al. 2008). Our f_host parameter can approximately account for the latter effect (if one thinks of some of the ‘off’ haloes as having hosted quasars so long in the past that their bubbles have recombined). We also expect that inhomogeneous recombinations are relatively unimportant for the quasar case, at least until the late stages of reionization (Furlanetto & Oh 2008a; Davies & Furlanetto 2014).

Because we lack a strong motivation for a particular value of ζ, we will generally use it as a normalization constant. That is, we will specify |$x_{\rm {He\,\small {iii}}}$| and a model for the host properties (M_min and f_host), and then choose ζ to make equation (1) true. As described above, any such choice can be made consistent with the luminosity function if we allow the quasar lifetime to vary.

Throughout most of the paper, we consider two basic source models. Neither is meant to be a detailed representation of reality: instead, they take crude, but contrasting, physical pictures meant to bracket the real possibilities. In the first, our fiducial method, we imagine that some quasars existed long enough ago that any He iii ions they created have since recombined, so only a fraction of the haloes (which recently hosted quasars) contribute to the ionization field. We therefore assume that during reionization (1) only a fraction f_host of haloes above |$M_{\rm min} = 5 \times 10^{11} \,\mathrm{M}_{\odot }$| have ever hosted sources, and that (2) each halo with a source, regardless of its mass, produces the same total number of ionizing photons. The latter assumption effectively fixes our f_q parameter, which must vary inversely with mass.

We are left with two free parameters, f_host and ζ. For the first, we set |$f_{\rm host} \approx x_{\rm He\,\small {iii}}/Q_{\rm ion}$| as the fraction of haloes that contain sources, where Q_ion is the (assumed) total number of ionizing photons produced per helium atom (integrated over all past times). Physically, this means that at the moment the universe becomes fully ionized we would have f_host ≈ 1/Q_ion. In this case, the other haloes would have had ionizing sources so far in the past that their He iii regions have recombined since formation. Since we loosely assume z = 3 to be the mid-point of helium reionization, we then fix ζ so as to produce |$x_{\rm {He\,\small {iii}}}$| = 0.5 with f_host = 0.167, where Q_ion = 3 as expected at z = 3. This efficiency effectively determines the size of each ionized bubble: thus, their number density and size together fix the net ionized fraction. To set other ionized fractions, f_host is varied, where intuitively more sources produce more ionizations.

In detail, these parameters cannot generally be set by analytic arguments, as the total ionized fraction depends on the numerical method as well (because our algorithm does not conserve or track individual photons). A final adjustment (of order unity) is necessary in order to reproduce the proper He iii fraction. (Hence, |$f_{\rm host} \approx x_{\rm He\,\small {iii}}/Q_{\rm ion}$| above is not a strict equality.) Thus, the physical parameters that motivate each realization should be taken only as rough guides to the underlying source populations.

Our second source model (which we will refer to as the abundant-source method) is chosen primarily as a contrast to this fiducial one. We assume that (1) all haloes above |$M_{\rm min} = 5 \times 10^{11} \,\mathrm{M}_{\odot }$| have hosted sources, even when the global |$x_{\rm {He\,\small {iii}}}$| is small, and that (2) the ionizing efficiency of each halo is proportional to its mass (which demands that f_q = 1). These assumptions mean that (when averaged over large scales) the quantity ζf_hostf_q in equation (1) reduces to the simple form ζf_coll (the same form typically used for hydrogen reionization; Furlanetto et al. 2004). We then choose ζ to satisfy equation (1) for a specified |$x_{\rm {He\,\small {iii}}}$|⁠.

Neither of these approaches is entirely satisfactory, but they bracket the range of plausible models reasonably well and so provide a useful gauge of the robustness of our predictions. Our fiducial model fails to match naive expectations in two ways. First, it assumes that only a fraction of massive haloes host black holes that have ionized their surroundings. One possible interpretation is that some of these sources appeared long ago, so that their ionized bubbles have recombined. However, in most cases the required number of recombinations would be much larger than expected (Furlanetto et al. 2008). A simpler interpretation is that many of these haloes simply lack black holes. While this fails to match observations in the local Universe (see, e.g. Gültekin et al. 2009 and references therein), the situation at z ∼ 3 is less clear, and it is possible that many massive haloes have not yet formed their black holes.

On the other hand, this procedure does allow our calculations to crudely describe the evolution of the He iii field with redshift, even though the calculations are all performed at z = 3. The primary difference between our density field and one at, say z = 4, is the decreased halo abundance in the latter. Assuming that f_host < 1, therefore, roughly approximates the halo field at this higher redshift (though, of course, it does not properly reproduce their spatial correlations, etc.).

A second problem with our fiducial model is that, by assuming that each halo produces the same number of ionizing photons, the resulting black holes will likely not reproduce the well-known M–σ relation between halo and black hole properties (Gültekin et al. 2009). Our choice is motivated by two factors. First, any scatter in the high-redshift black hole mass–halo mass relation (as is likely if these haloes are still assembling rapidly) may dominate the steep halo mass distribution. Secondly, some models predict a characteristic halo mass for quasar activity (Hopkins et al. 2006). In practice, this inconsistency is probably not an enormous problem, because our high-mass threshold (and the resulting steep mass function shown in Fig. 1) means that only a fairly narrow range of haloes contribute substantially.

The abundant-source method at least qualitatively reproduces the relation between halo mass and black hole mass (and the expected ubiquity of quasars), but it leads to some other unexpected features. Most importantly, it demands that the ionized bubbles around many quasars be quite small when |$x_{\rm {He\,\small {iii}}}$| <1, as we will show later. Partly, this is because we work at z = 3, where the number of massive haloes is rather large. At higher redshift, an identical prescription would find fewer sources, which would allow for larger individual bubbles. Moreover, the characteristic luminosity of quasars does not vary dramatically throughout the helium reionization epoch, so we naturally expect that the average bubble size will also remain roughly constant (M09).

Instead, in our calculations within this model, the assumed efficiency parameter ζ varies with |$x_{\rm {He\,\small {iii}}}$| while the source distribution remains fixed. Thus, our sequence of ionized fractions should not be viewed as a sequence over time in the abundant-source model, as we do not allow the halo field to evolve. One could view the small sizes of the bubbles when |$x_{\rm {He\,\small {iii}}}$| < 1 as a result of strong ‘internal absorption’ within the bubbles, which would keep them compact, but that does not seem a likely scenario (Furlanetto & Oh 2008a).

Nevertheless, it is useful to consider a contrasting case from our fiducial model, in which reionization proceeds not by generating more He iii regions but by growing those that do exist. Our abundant-source model provides just such a contrast and so is useful as a qualitative guide to the importance of the source prescription.

Ionization field

First, the source field is smoothed to a lower resolution (⁠|$N_{\rm He\,\small {iii}} = 500^3$| grid cells) to speed up the ionized bubble search procedure, which takes considerable time at high ionized fractions relative to the other components of the code. We then filter this source field using a real-space, top-hat filter, starting at a (mostly irrelevant) large-scale R_max and decreasing down to the cell size. Every region satisfying equation (2) at each filtering scale is flagged as a He iii bubble, regardless of overlap (unlike the halo-finding case). The result is a two-phase map of the ionization field, dependent on the number density of haloes hosting ionizing sources.

The assumption of two phases with sharp edges is clearly a simplification. One-dimensional radiative transfer codes show that the ‘transition region’, or distance between fully ionized and singly ionized helium zones, is typically a few Mpc thick – more precisely, the radius over which the He iii fraction falls from ∼0.9 to ∼0.1 is ≲ 20 per cent of the size of the fully ionized zone (Davies, private communication; see also fig. 19 of M09). The relatively late timing of He ii reionization, the rarity of the ionizing sources, and the hard spectra of quasars are responsible for widening this region. In many applications, especially those involving the IGM temperature, this transition region is very important. However, for other applications – such as the He ii Lyα forest – the two-phase approximation is adequate, because even a small He ii fraction renders a gas parcel nearly opaque. Fortunately, most He iii regions are very large (many tens of Mpc across), where many are built through the overlap of bubbles from individual sources, so the relatively narrow transition regions do not dominate the ionization topology. Nevertheless, we must bear in mind this important simplification when applying our semi-numeric simulations to observables.

In the ionization field slices shown in Fig. 2 , we see the increasing ionized fraction through a static universe for the fiducial method. Displayed in the figure are |$x_{\rm {He\,\small {iii}}}$| = 0.20, 0.50, 0.80, and 0.90, corresponding to f_host = 0.049, 0.167, 0.400, and 0.560. The map is fairly bubbly in nature with noticeable clustering. The fact that a pixel residing in a bubble has an increased likelihood of being near another bubble has important consequences for our later calculations. Even at the highest |$x_{\rm {He\,\small {iii}}}$| displayed, there are large ‘neutral’ sections. Note the periodic boundary conditions.

Figure 2.

Slices from the ionization field at |$x_{\rm {He\,\small {iii}}}$| = 0.20, 0.50, 0.80, and 0.90 (from left to right) for the fiducial model. All slices are 250 Mpc on a side and 0.5 Mpc deep. White (black) indicates He iii (He ii) regions.

Open in new tabDownload slide

Figure 3.

Slices from the ionization field at |$x_{\rm {He\,\small {iii}}}$| = 0.21, 0.50, 0.80, and 0.90 (from left to right) for the abundant-source model. All slices are 250 Mpc on a side and 0.5 Mpc deep. White (black) indicates He iii (He ii) regions.

Open in new tabDownload slide

Beyond our two main methods, we consider some additional parameter choices, shown in Fig. 4. In the left three panels, we explore a range of minimum halo mass (effectively changing the quasar lifetime) for the abundant-source method. From left to right, |$M_{\rm {min}} = (0.5, 5, 10)\times 10^{11} \,\mathrm{M}_{\odot }$| at |$x_{\rm {He\,\small {iii}}}$| ≈0.80. As M_min is decreased, the number of sources increases, which implies that the quasar is ‘turned on’ for a shorter period of time, ionizing a smaller region given a fixed |$x_{\rm {He\,\small {iii}}}$|⁠. Similarly, a higher minimum mass means fewer sources and larger, rarer ionized bubbles, which is more pronounced at lower |$x_{\rm {He\,\small {iii}}}$|⁠. The rightmost panel of Fig. 4 shows the result for the fiducial method with f_host = 1 and |$x_{\rm {He\,\small {iii}}}$| ∼0.80. Importantly, the second panel (abundant-source model with |$M_{\rm {min}} = 5\times 10^{11} \,\mathrm{M}_{\odot }$|⁠) and this panel are very similar visually, indicating the exact ionization criterion is unimportant.

Figure 4.

Slices from the ionization field for variations on the abundant-source model at |$x_{\rm {He\,\small {iii}}}$| ≈ 0.80. From left to right, |$M_{\rm min} = (0.5, 5, 10)\times 10^{11} \,\mathrm{M}_{\odot }$|⁠. The rightmost panel shows the fiducial model with f_host = 1 at |$M_{\rm min} = 5\times 10^{11} \,\mathrm{M}_{\odot }$|⁠.

Open in new tabDownload slide

To quantify the statistical properties of our ionization field, we calculate the dimensionless ionized power spectrum |$\Delta ^2_{x_{\rm He\,\small {iii}}}(k) = k^3/(2\pi ^2V) \langle |\delta _{x_{\rm He\,\small {iii}}}(\boldsymbol {k})|^2\rangle _k$|⁠, where |$\delta _{x_{\rm He\,\small {iii}}}(\boldsymbol {x}) \equiv x_{\rm He\,\small {iii}}(\boldsymbol {x})/\bar{x}_{\rm He\,\small {iii}} - 1$|⁠. The results for the fiducial model at |$x_{\rm {He\,\small {iii}}}$| = 0.20, 0.50, 0.80, and 0.90 are shown by the dot-dashed, long-dashed, short-dashed, and solid lines, respectively, in the left-hand panel of Fig. 5. These curves are qualitatively similar to M09 results (except for the lowest |$x_{\rm {He\,\small {iii}}}$|⁠, which we expect to differ the most), with a large-scale peak and small-scale trough. Differences in the underlying quasar models affect the exact position of the peak, but our results are comparable: ∼30 Mpc scales compared to ∼60 Mpc in M09. Note that the positions of the peaks and troughs do not vary with the He iii fraction, which is consistent with M09. Since, in this fiducial model, we are only changing the number of haloes that contribute to the ionization at different |$x_{\rm {He\,\small {iii}}}$|⁠, the size of the ionized bubbles remain constant, leading to the lack of evolution in the peak scale with |$x_{\rm {He\,\small {iii}}}$|⁠. A similar phenomenon occurs in M09, because the typical luminosity of quasars remains roughly constant with redshift (as does the effective lifetime in their model), so the size of a typical ionized bubble remains the same throughout reionization.

Figure 5.

The power spectrum of |$x_{\rm {He\,\small {iii}}}$| at various ionized fractions. Left-hand panel: shown are |$x_{\rm {He\,\small {iii}}}$| = 0.20 (dot-dashed), 0.50 (long-dashed), 0.80 (short-dashed), and 0.90 (solid) for the fiducial model. Right-hand panel: the abundant-source method for |$x_{\rm {He\,\small {iii}}}$| = 0.21 (dot-dashed), 0.50 (long-dashed), 0.80 (short-dashed), and 0.90 (solid) is plotted. The dotted lines represent the inclusion of even more sources (⁠|$M_{\rm min} = 0.5\times 10^{11} \,\mathrm{M}_{\odot }$|⁠) for |$x_{\rm {He\,\small {iii}}}$| = 0.50 and 0.80.

Open in new tabDownload slide

The right-hand panel of Fig. 5 shows the results for the abundant-source model with |$x_{\rm {He\,\small {iii}}}$| = 0.21, 0.50, 0.80, and 0.90 (dot-dashed, long-dashed, short-dashed, and solid lines, respectively). The peaks are wider and evolve to smaller scales as |$x_{\rm {He\,\small {iii}}}$| decreases, since the lower ionization states are dominated by smaller and smaller bubbles. As expected from the ionization morphologies, the lower |$x_{\rm {He\,\small {iii}}}$| exhibit the largest differences from the fiducial model. The dotted lines assume the lowest minimum halo |$M_{\rm min} = 0.5 \times 10^{11} \,\mathrm{M}_{\odot }$| at |$x_{\rm {He\,\small {iii}}}$| = 0.50 and 0.80, as shown in the leftmost panel of Fig. 4. Since each ionized bubble is smaller than the |$M_{\rm min} = 5 \times 10^{11} \,\mathrm{M}_{\odot }$| case, the peak in the spectrum is shifted to the right and generally smoother. Omitted from this panel, the power spectra of the other variations outlined in Fig. 4 are nearly indistinguishable.

Density field

We use the density field not just to find dark matter haloes and ionized regions but also to track photon transfer through the IGM. Specifically, when calculating the photoionization rate for He ii later in this paper, we wish to follow the penetration of hard photons into neutral regions. This task requires a reasonably accurate small-scale density field appropriate to baryons. This density field is also crucial to investigating other observable quantities, such as the He ii Lyα forest. In this section, we describe our procedure for generating that distribution.

As mentioned above, the density field in dexm is generated through first-order perturbation theory. Specifically, matter particles on the high-resolution (2000³) grid are moved from their Lagrangian to their Eulerian coordinates using the displacement field (Zel'Dovich 1970), and then binned on to the lower resolution (500³) grid. The resulting field has an inherent particle mass quantization and can be smoothed with a continuous filter to obtain a continuous density field (as in all particle dynamics codes). Ideally, the filter scale should correspond to the Jeans mass. We therefore smooth the density field by applying a real-space, top-hat filter over 0.75 Mpc (or 1.5 pixel lengths). Note that the Jeans length in ionized gas at the mean density with an IGM temperature of ∼10⁴ K is ∼1 Mpc (near our smoothing size of 0.75 Mpc). The pre- and post-smoothing density probability distribution functions are shown in Fig. 6 with the dotted and dashed lines, respectively.

Figure 6.

The volume-weighted density distribution before, during, and after our mapping procedure with a reference curve. The dotted (short-dashed) line is the evolved (then smoothed) density distribution. The solid line shows the distribution used for subsequent calculations, which is mapped on to the BB09 fit, shown as the long-dashed line.

Open in new tabDownload slide

This approach tends to underpredict overdensities and overpredict voids, although Mesinger, Furlanetto & Cen (2011) show that during hydrogen reionization (z ≳ 6), the resulting density field agrees with one generated from a hydrodynamical simulation at the per cent level to about a dex around the mean density (spanning a large majority of the IGM). However, at the redshifts of interest for He ii reionization, the density field becomes more non-linear, and our first-order perturbation theory approach is less accurate. Rather than adding higher order corrections, we apply a simple mapping procedure to the smoothed particle hydrodynamics results of Bolton & Becker (2009, hereafter BB09), using the respective cumulative probably distributions of the overdensity p_DF(< Δ_ρ) from our calculation and p_BB(< Δ_ρ), the fit from BB09. Here, |$\Delta _{\rho } = \rho /\bar{\rho }$|⁠. The long-dashed curve in Fig. 6 is the BB09 fit. For each cell's overdensity |$\Delta _{\rho }^i$|⁠, we assign a new overdensity |$\Delta _{\rho }^j$| such that |$p_{\rm DF}(<\Delta _{\rho }^i) = p_{\rm BB}(<\Delta _{\rho }^j)$|⁠. In this way, we preserve the ordering of our underlying matter distribution, but we also achieve a result closer to the gas distribution, which is more appropriate for computing observables.

The final, mapped density distribution is the solid curve in the figure. The discreteness visible in the dotted line and a lack of very large overdensities render the matching procedure imperfect. The choice of filter scale has minimal effect on our later calculations. Compared to the smoothed field without mapping, the mapped distribution slightly narrows the photoionization curves in Section 4, since higher density regions effectively absorb more photons.

QUASAR MODELS

We are interested not only in how the source field affects the morphology of ionized gas but also in how it affects the metagalactic radiation background. This background directly affects the Lyα forest and also indirectly affects several interesting observables, like the temperature of the IGM and ionization states of metal line absorbers. Although obviously intertwined with the determination of the ionization field, we must take an additional step to determine this background, which depends on the instantaneous luminosities of the sources. Conversely, the ionization maps depend on the history of ionizing photons. We follow an empirical approach to determine the number and intrinsic properties of the ‘active’ quasars, using an observationally determined QLF. Then, we compare three basic models for placing these quasars in dark matter haloes.

Now that each quasar has a designated specific luminosity, L_ν (given by L_B and α), the three quasar models differ in the method for placing them inside haloes. Note that the two prescriptions for the ionization morphology only determine which host haloes have ‘turned on’ by z = 3 via f_host (and M_min), so active quasars are always placed in ionized bubbles.

QSO1. Haloes are randomly populated with no preference to host mass. This method is inspired by Hopkins et al. (2006). In this picture, the distribution of quasar luminosities is due primarily to evolution in the light curve of individual black holes rather than differences in the black holes themselves: as quasar feedback clears out the host halo's gas, less of the light is obscured, rapidly increasing the flux reaching the IGM (Hopkins et al. 2005a,b). This scenario is consistent with quasars residing in haloes of a characteristic size; by choosing random haloes, we preferentially populate the lowest halo mass, or 5 × 10¹¹ M_⊙, which is at the lower end of recent estimates (see, e.g. Ross et al. 2009; Trainor & Steidel 2012). As shown above, a larger M_min would produce rarer and larger bubbles, but our conclusions would be minimally affected as detailed below.

QSO2. In the second model, we assign the brightest quasars to the most massive haloes, which is more consistent with the well-known M–σ relation. Loosely, we are supposing that the luminosity of quasars increases with the halo mass, without fixing the exact form. This picture essentially assumes that the quasar's luminosity is proportional to its black hole's mass, and that the black hole mass increases monotonically with halo mass (see, e.g. Wyithe & Loeb 2002). Any other factors – such as the light curve of each quasar – are subdominant in setting the relationship.

QSO3. The last model assigns quasars to completely random positions within our box, irrespective of halo positions or ionized regions. This model does not include any clustering and serves mainly as an illustration.

PHOTOIONIZATION RATE

With the quasar configuration, ionization map, and density distribution, we can calculate the He ii photoionization rate Γ at any point in our box, even during reionization. Although the mean rate 〈Γ〉 determines the average behaviour of many observable quantities, we are particularly interested in the magnitude of fluctuations about this mean, which have important effects for scatter in observables and even for the evolution of the mean (Davies & Furlanetto 2014; Dixon, Davies, & Furlanetto, in preparation). In brief, we add up the contribution to Γ from all sources within our box, taking into account the local density and the ionization state.

Post-reionization

The mean free path is both uncertain and varies rapidly over the era of interest, with estimated values ranging from 15 to 150 Mpc at the ionization edge. We choose a fiducial value at this edge, λ₀ = 60 Mpc, to match more detailed calculations of λ_mfp that incorporate more sophisticated treatments of the Lyα forest (Davies & Furlanetto 2014). However, it is worth noting that the mean free path varies quite rapidly with redshift over the z ≈ 2.5–3.5 range and is subject to substantial uncertainties in the source population. We allow a range of λ₀ = 15, 35, 60, 80 Mpc in this work to evaluate the impact.

In detail, our assumption of a spatially constant mean free path is incorrect: because the ionization structure of the IGM absorbers is set by the (fluctuating) metagalactic radiation field, the amount of absorption – and hence the mean free path – fluctuates as well. We ignore this possibility here, but it has been approximately included in analytic models (Davies & Furlanetto 2014).

Figure 7.

Distribution of Γ (scaled to the mean) in the post-ionization regime. Left-hand panel: three models are shown: quasars placed in random haloes (black), most luminous quasars situated in most massive haloes (blue), and quasars placed randomly throughout the box (red). The dashed curves result from a fixed attenuation length, λ_mfp = 60 Mpc. The post-reionization analytic result is represented by the dotted line. Right-hand panel: from widest to narrowest curves, λ₀ = 15, 35, 60, and 80 Mpc.

Open in new tabDownload slide

Since the results for QSO1 and QSO3 are nearly identical, randomly placed quasars approximates quasars placed in random haloes, demonstrating that clustering of sources is unimportant. Predictably, the QSO2 result is wider, since clustering is more pronounced in this scenario. The largest dark matter haloes tend to be found near other large dark matter haloes, so the brightest quasars will tend to cluster. As noted in Furlanetto (2009), these brightest quasars dominate the distribution. The 〈Γ〉 are identical to within a few per cent for the three quasar models.

The left-hand panel of Fig. 7 also compares frequency-dependent (solid) and frequency-independent (dashed) mean free paths. The differences between quasar models are similar in both cases, but the curves are generally wider for a frequency-independent value. Photons at the highest energies can travel further in the frequency-varying model providing a higher minimum ionizing background, although the difference is small in the post-reionization regime. Along similar lines, the right-hand panel of Fig. 7 shows the impact of varying the mean-free-path normalization with λ₀ = 15, 35, 60, and 80 Mpc (widest to thinnest curves). As before, larger λ_mfp narrow the distribution. The larger λ₀ results are very similar, indicating that the quasars dominating the distribution are within ∼60 Mpc of each other; beyond that point, the intrinsic spread in luminosities dominates the distribution.

The analytic model of Furlanetto (2009), based on Poisson distributed sources (Zuo 1992; Meiksin & White 2003) and represented by the dotted line, with the a fixed λ_mfp = 60 Mpc is remarkably close to the semi-numeric result. Since it does not include frequency-dependence in λ_mfp, the analytic result unsurprisingly matches the fixed λ_mfp model best. The congruence of the analytic and semi-numeric distributions further indicates that the rarity of quasars dominates the distribution of the photoionization rate, not clustering.

Even in this post-reionization regime, Γ varies substantially about the mean, with very little dependence on the underlying details. The width of this distribution is essentially set by the QLF (constant for all models in this work) and the mean free path, where smaller values broaden the curve. In all cases, the high-Γ tail is especially robust. This segment of the distribution corresponds to regions near bright sources, or proximity zones, so model details, like the mean free path or exact quasar placement, are inconsequential.

During reionization

Next, we approximate the evolution of Γ during reionization (by varying |$x_{\rm {He\,\small {iii}}}$|⁠). High-energy photons can propagate through the He ii regions present in this era, because the optical depth experienced by photons decreases with frequency. Unlike previous analytic models, we can explicitly compute the ionizing background generated by these high-energy photons. In particular, with our density field and ionization map, we estimate the He ii column density |$N_{\rm {He\,\small {ii}}}$| for a given photon ray and, thereby, estimate the absorption.

Suppose we wish to calculate Γ at a particular point |${{\boldsymbol r}}$|⁠. Beginning there, we take a step of length |${\rm d}x_{\rm p} = \frac{1}{2} \Delta x_{\rm p}$|⁠, where Δx_p is the proper width of a pixel in our simulation. If the pixel contains He ii, the step's contribution to the intervening column density is |$\Delta N_{\rm {He\,\small {ii}}} = \Delta _{\rho } \times n_{\rm {He}} \times {\rm d}x_{\rm p}$|⁠, where Δ_ρ is the pixel's overdensity and n_He is the (proper) number density of helium atoms. Otherwise, there is no contribution to |$N_{\rm {He\,\small {ii}}}$|⁠. If the ionization state changes during a step, we use the |$x_{\rm {He\,\small {ii}}}$| present at the mid-point for the entire segment.

This component of the absorption typically dominates the opacity when |$x_{\rm {He\,\small {iii}}}$| is small: at |$x_{\rm {He\,\small {iii}}}$| ≈0.50, the ionized bubble size distribution, a reasonable proxy for the path length, peaks around 10 Mpc, which is much smaller than the mean free path we have assumed within the ionized gas. However, by |$x_{\rm {He\,\small {iii}}}$| ≈ 0.80, the peak is of the order of the mean free path. Thus, photons typically reach τ ∼ 1 while travelling through nominally ionized regions during these late stages, indicating that dense clumps within these ionized regions become important before the completion of reionization.

Fig. 8 displays the distribution of Γ for several |$x_{\rm {He\,\small {iii}}}$|⁠. In the left-hand panel, QSO1 (solid) and QSO2 (dashed) for the fiducial model are shown for |$x_{\rm {He\,\small {iii}}}$| = 0.20, 0.50, 0.80, 0.90, and post-reionization (from widest to thinnest distributions, respectively). A bimodal distribution develops at small ionized fractions due to the low-level background from hard photons escaping into the He ii regions. Because He ii reionization causes substantial IGM heating, whose amplitude is determined by the shape of the local ionizing background, this bimodality may have important consequences for this heating process (Furlanetto & Oh 2008a; Bolton, Oh & Furlanetto 2009; M09). Even with |$x_{\rm {He\,\small {iii}}}$| = 0.9, this hard photon background is still substantial. Interestingly, Mesinger & Furlanetto (2009) calculated the UV background using dexm during the H i reionization era. At that time, the background does not exhibit this behaviour. This difference is due to our inclusion of hard photons, which are more important for helium reionization (as the H i ionizing sources are expected to have had very soft spectra).

Figure 8.

Distribution of Γ (scaled to the mean fiducial value post-reionization) for several |$x_{\rm {He\,\small {iii}}}$| and varying λ₀. Left-hand panel: from lowest to highest at 〈Γ〉, |$x_{\rm {He\,\small {iii}}}$| = 0.20, 0.50, 0.80, 0.90, and post-reionization are shown for the fiducial model (solid) and QSO2 variation (dashed). Right-hand panel: the fiducial model for |$x_{\rm {He\,\small {iii}}}$| = 0.50, 0.80, and post-reionization is displayed with varied mean-free-path normalizations, λ₀ = 15, 35, 60, and 80 Mpc as long-dashed, short-dashed, solid, and dotted, respectively. Note that λ₀ = 60 Mpc (solid) is the fiducial value.

Open in new tabDownload slide

The choice of quasar model has only a minimal impact on the distributions, but it preferentially affects the high-end tail, which becomes increasingly important as |$x_{\rm {He\,\small {iii}}}$| decreases. Put another way, the high-Γ peak is shifted to the right for QSO2 for |$x_{\rm {He\,\small {iii}}}$| <1, while the low-Γ peak is unaffected. Since QSO2 clusters more bright sources together, this shift indicates that proximity to quasars trumps the radiation from more distant sources. This behaviour implies that detailed observations of He ii Lyα absorption inside of large ionized bubbles during reionization may help determine how quasars are distributed inside of dark matter haloes.

The right-hand panel of Fig. 8 illustrates the effect of changing λ₀ on the distribution of Γ, where λ₀ = 15, 35, 60, and 80 Mpc as long-dashed, short-dashed, solid, and dotted, respectively, for |$x_{\rm {He\,\small {iii}}}$| = 0.50, 0.80, and post-reionization. Note that the distribution is scaled to the post-reionization 〈Γ〉 for each λ₀. In the post-reionization regime, the general trend is that larger mean free paths translate into narrower distributions, as previously demonstrated. Similarly to post-reionization, the high-Γ tail is steeper for larger λ₀ during reionization. In contrast, the low-Γ peak does not shift to lower values for smaller λ₀, nor do the curves become significantly narrower. Since more photons can penetrate the He ii regions with larger mean free paths, the 〈Γ〉 of these regions (a proxy for the low-Γ peak) is higher. To roughly summarize, the low peak and steepness of the high-Γ tail are determined by how far photons can travel, and the high peak is shaped by the position of quasars (see Mesinger & Dijkstra 2008 for a similar conclusion for the hydrogen case).

The left-hand panel of Fig. 9 contrasts the abundant-source and fiducial methods. Here, |$x_{\rm {He\,\small {iii}}}$| = 0.21, 0.50, 0.69, 0.80, and post-reionization (widest to thinnest curves, respectively) with QSO1 (solid) and QSO2 (short-dashed). In the abundant-source case, the difference between QSO1 and QSO2 is even less pronounced than in the fiducial framework. With smaller ionized bubbles, fewer pixels reside in an ionized region with an active quasar, except when overlap dominates (higher |$x_{\rm {He\,\small {iii}}}$|⁠). Comparing the fiducial method with QSO1 (long-dashed curves), the abundant-source model are very similar for higher |$x_{\rm {He\,\small {iii}}}$|⁠. For the lowest ionization fraction, the high-Γ tail is substantially higher in the fiducial model, but the low-Γ peak remains quite similar. The decrease at high intensities occurs for two reasons. First, the abundant-source method enforces smaller ionized bubbles, so the proximity zones are smaller and less space is filled by the near regions (outside the He iii region, of course, the opacity rapidly increases). Secondly, the abundant-source method has many more empty He iii regions – that is, regions that have been ionized in the past but no longer host an active source. Because more of the IGM is inside of isolated bubbles, there is less chance for a nearby quasar to illuminate such a region. Overall, the difference between the two methods is surprisingly small, especially during the late stages of reionization.

Figure 9.

Distribution of Γ (scaled to the mean value post-reionization) for the abundant-source method. Left-hand panel: from lowest to highest at 〈Γ〉, |$x_{\rm {He\,\small {iii}}}$| = 0.21, 0.50, 0.80, 0.90, and post-reionization, where QSO1 (solid) and QSO2 (short-dashed). For comparison, the long-dashed curves are the fiducial model with QSO1. Right-hand panel: here, the minimum halo mass is varied with QSO1, and |$M_{\rm min} = (0.5, 5, 10)\times 10^{11} \,\mathrm{M}_{\odot }$| are the short-dashed, solid, and long-dashed curves, respectively. For clarity, only |$x_{\rm {He\,\small {iii}}}$| = 0.50, 0.80, and post-reionization (from widest to narrowest) are shown.

Open in new tabDownload slide

The right-hand panel of Fig. 9 shows the impact of the number of haloes hosting sources (in the past) via the minimum halo mass |$M_{\rm min} = (0.5, 5, 10)\times 10^{11} \,\mathrm{M}_{\odot }$| as the short-dashed, solid, and long-dashed curves, respectively. Note that these are variations on the abundant-source model, where the representative ionization map slices are shown in Fig. 4. The widest to thinnest curves are |$x_{\rm {He\,\small {iii}}}$| = 0.5, 0.8, and post-reionization. In the post-reionization regime, the results are nearly identical, meaning the exact placement of active quasars is relatively unimportant. With smaller ionized bubbles (i.e. lower M_min), the distributions are shifted towards lower Γ. Decreasing M_min reduces source clustering. Essentially, fewer regions are within an ionized bubble and close to a quasar.

Figs 8 and 9 span a wide range of quasar model possibilities. Importantly, the density distribution (as noted in Section 2.4), the exact ionization map, and even the basic quasar model do not greatly affect the photoionization rate distributions. Not only does this minimize the dependence on our fiducial choices, it also emphasizes the robustness of our general conclusions and the relative importance of uncertain quantities relating to quasars.

To tease out the cause of the bimodal nature of Γ, we contrast the distribution of Γ inside fully ionized bubbles with that outside them in Fig. 10 (short- and long-dashed curves, respectively). As expected, ionized areas around active quasars contain the highest Γ. The origin of the small-Γ distribution is less straightforward. The ionizing radiation that penetrates the singly ionized helium comprises the majority of this low-end segment. This radiation is mainly due to high-frequency photons that can travel larger distances before absorption. Ionized bubbles that do not host active quasars also make up a non-negligible fraction of this low-Γ tail (and especially in the transition region from one peak to the other). This contribution is more relevant early in reionization (when overlap of the bubbles is less significant, meaning neighbouring sources cannot illuminate such empty regions). We caution the reader that such empty regions are unlikely to be fully ionized, since the recombination time is relatively short (Furlanetto et al. 2008). Therefore, we likely overestimate the ionizing background here.

Figure 10.

The distribution of Γ within and without ionized regions. The solid, short-dashed, and long-dashed lines are the distributions of the entire box, He iii regions, and He ii regions, respectively. The widest to thinnest curves are |$x_{\rm {He\,\small {iii}}}$| = 0.50, 0.80, and post-reionization. The high-end Γ is determined almost exclusively by regions around quasars, which are necessarily inside ionized bubbles. The ionizing photons that escape the He iii dominate the low-end Γ, although ionized bubbles without an enclosed source also contribute.

Open in new tabDownload slide

Furlanetto (2009) explored the related f(J) through a combination of analytic calculations and a Monte Carlo model for the number of sources within a distribution of bubble sizes. The results during reionization show large fluctuations but do not exhibit the bimodal behaviour found in this work. The difference is due to our inclusion of local density effects and the improved treatment of hard photons via ν_min and a frequency-dependent mean free path.

COMPARISON WITH H i REIONIZATION

It is useful to consider the similarities and differences between our work here and analogous studies of H i reionization, especially given that we are adapting a hydrogen reionization code. In this section, we identify the major points of divergence between these two epochs.

Our semi-numeric approach was originally designed for calculations of H i reionization, to which it is undoubtedly better suited. Most importantly, the two-phase approximation (fully ionized or fully neutral) is very accurate for H i reionization, because the soft stellar sources most likely responsible for H i reionization have very short mean free paths ( ≲ a few kpc). The ‘photon-counting’ arguments inherent to the method are therefore almost exactly accurate for that case, whereas the transition region between He ii and He iii fills a non-negligible fraction of the volume. Our ionization maps must therefore be taken as approximate, although they are still useful for many purposes (such as calculations of the He ii Lyα optical depth, which is very large even in the transition region).

Another key difference between the ionization fields in H i and He ii reionization is the source abundance. We typically assume that any dark matter halo able to form stars will also produce ionizing photons during the earlier era, making those sources many times more common than luminous quasar hosts at z ∼ 3. During H i reionization, ionized bubbles typically have thousands of sources even early in the process (Furlanetto et al. 2004), but we have found that, for He ii reionization, even a single source can generate a large bubble. The ‘morphology’ of reionization (i.e. the distribution of ionized bubble sizes) is often taken as a key observable of the H i process, because the huge number of sources closely connects it to the underlying dark matter field. In the He ii case, however, the small numbers of sources make that connection mostly uninteresting: Poisson fluctuations instead play a key role (Furlanetto & Oh 2008a), and the distribution of bubbles relative to haloes is much more likely to tell us about the details of the black hole–host relation than it is the fundamentals of the reionization process. Quantitatively, there is little evolution in the shape of the power spectrum of the ionized fraction throughout helium reionization.

Additionally, because there are so many sources inside each H ii region, the ionizing background remains large at all times, preventing these bubbles from recombining substantially. When only one or two (short-lived) sources occupy each bubble, however, recombinations can substantially increase the He ii fraction inside even regions that are nominally fully ionized (Furlanetto et al. 2008). Fortunately, Fig. 10 shows that most such regions will be illuminated, at least late in the process. But this will, for example, increase the opacity of some He iii bubbles in the Lyα forest (and to ionizing photons). We have crudely modelled this effect through our source prescriptions (one can regard turning some haloes ‘off’ as instead allowing their ionized bubbles to recombine fully).

Our approach to describing the ionizing sources themselves is also rather different. The lack of data on high-redshift galaxies so far implies that naive prescriptions (with a constant ionizing efficiency per halo, or at most a systematic variation with halo mass) are more than adequate. On the other hand, the tight observational constraints on the QLF at z ∼ 3 require that it be directly incorporated into the calculation. This inclusion is easy to do for the instantaneous ionization rate but much harder for the cumulative number of ionizing photons, because the relationship between the visible sources and their dark matter halo hosts and the light curves of individual quasars remains uncertain. We have shown that the details of this association do not strongly affect many properties of reionization (such as the distribution of Γ), but they will show up in other observables. For example, one could imagine searching for massive haloes inside of large He iii regions in order to identify the descendants of bright quasar hosts.

A final key difference is the source spectrum, which is unimportant for stellar sources but crucial for the quasars responsible for He ii reionization. Although we have ignored the resulting partially ionized transition zones, we have included these spectra explicitly in evaluating the ionizing background. We are thus able to track the growth of a hard background in He ii regions, an effect likely to be much less important during H ii reionization (Mesinger & Furlanetto 2009).

DISCUSSION

We have introduced a new and efficient semi-numeric method for efficiently generating dark matter halo distributions and ionization maps relevant to the full reionization of helium. Specifically, we adapted an existing hydrogen reionization code (dexm) to identify haloes at lower redshift and changed the method for determining the ionization state to suit quasars as sources. The speed of our algorithm allows exploration of a wide range of possible quasar models not only for active quasars, but in the realization of the ionization map as well.

Our results broadly match previous simulations of related quantities, at least statistically speaking. We generated accurate halo mass functions as compared to N-body simulations. This congruence is important for determining the effect of quasar clustering. We also reproduced the He iii morphology, especially at high |$x_{\rm {He\,\small {iii}}}$|⁠, as demonstrated by the ionized power spectra, which compare favourably with M09 (modulo our different source prescriptions). We demonstrated that this method can easily produce simulation volumes hundreds of Mpc across to accurately capture the large-scale features of helium reionization. We also constructed a density field that matches (by design) the smoothed particle hydrodynamics simulations of BB09. A realistic density field is crucial to capturing at least some of the small-scale fluctuations in the observable quantities, which is lacking in semi-analytic models.

We have presented the first systematic study of fluctuations in the ionizing background during reionization. Our f(Γ) exhibit significant variation even post-reionization, in excellent agreement with analytic calculations (Furlanetto 2009). Furthermore, we calculate the more relevant f(Γ) as opposed to f(J) in the analytic case. We demonstrated that complications like a frequency-dependent mean free path and quasar clustering do not significantly affect the result in this regime.

We also find that a broader and bimodal distribution develops during reionization. This distribution differs from that during the hydrogen reionization era (Mesinger & Dijkstra 2008) and from analytic predictions for the He iii era (Furlanetto 2009). These fluctuations have important consequences, which we can quantify, for many observables, including the quasar spectra metal lines and IGM heating. In an upcoming paper, we will study the impact of these fluctuations on the He ii Lyα forest. The most important conclusion is that not every part of the IGM experienced the same reionization history or receives the same background radiation.

We included a frequency-varying mean free path, another crucial improvement on semi-analytic models. This frequency-dependent spectrum gives rise to the bimodal nature of f(Γ) we find during reionization. Specifically, hard photons can travel farther and escape into the He ii regions, creating a low-level background there. This frequency-dependence is important for studies of the He ii Lyα forest, metal lines, and IGM heating.

There are some important caveats to our methods as follows.

• We assign the quasars a simple – though frequency-dependent – mean free path with an uncertain normalization. This assumption ignores detailed radiative transfer effects and is important both during and after reionization. In particular, we prescribe the normalization of this mean free path, rather than allowing it to vary with the local radiation field (Davies & Furlanetto 2014).

• We assume a two-phase medium. In principle, our ionized bubble edges should not be abrupt, especially since helium recombines so quickly, but instead transition smoothly from fully ionized to singly ionized. This effect is less important during the late stages of reionization as the bubbles increasingly overlap (so that the ratio of their volume to their surface area decreases).

• We do not fully account for recombinations. As mentioned above, we do not explicitly include any He ii inside ionized bubbles, though likely some He iii would have recombined. Not only would this affect our effective mean free path (as determined by ν_min), we do not compensate for ionizing photons being ‘wasted’ on reionizing this gas when generating our ionization maps.

• Our simulations assume that quasars emit radiation isotropically. Many theoretical models pertaining to quasar behaviour predict otherwise. Although this effect is important during the early stages of reionization, once each ionized bubble contains (or has contained) multiple, randomly oriented quasars, beaming is no longer important.

Our conclusions are mostly insensitive to our underlying assumptions for active quasars and ionizing sources. For our wide range of explored models, the resultant f(Γ) varied minimally, except for λ₀. In particular, quasar clustering is a secondary effect, and a few rare, bright quasars dominate the spectrum. Although this makes our predictions robust, distinguishing between our models using observational data dependent on the ionizing background would be difficult if not impossible.

With the ionization maps, density field, and active quasar models, we are uniquely suited to study fluctuations in the He ii Lyα forest. In particular, we maintain spatial information that is important for considering correlations along a line of sight. We explore these fluctuations in an upcoming paper (Dixon, Davies, & Furlanetto, in preparation), where we show that the large variations in transmission indicate ongoing reionization at z ≳ 2.8. Similarly, these methods could apply to the He ii Lyman β forest, as measured by Syphers et al. (2011a).

Further applications that we intend to explore are the transverse proximity effect, as measured in the He ii Lyα forest (Jakobsen et al. 2003; Worseck & Wisotzki 2006; Worseck et al. 2007), whereby a quasar near another line of sight enhances the ionizing background. The size and nature of the apparent influence of the Lyα spectra may give important clues about the intrinsic properties of quasars (e.g. Furlanetto & Lidz 2011). Also, the spatial inhomogeneities of the UV radiation background should affect the ionization balance of the metals in the IGM. As an example, there may be a break in the ratio of C iv to Si iv at z ∼ 3 indicating helium reionization (Songaila 1998, 2005, though see Aguirre et al. 2004). The fluctuations we find in our He-ionizing background could have important consequences for these measurements. Combining past and upcoming observations with our theoretical tools will shed light not only on the timing and duration the helium reionization epoch, but also the nature of the IGM and the properties of quasars.

We thank F. Davies for helpful discussions and sharing preliminary results. This research was partially supported by the David and Lucile Packard Foundation, the Alfred P. Sloan Foundation, and NSF grant AST-0607470. KLD gratefully acknowledges support by the Science and Technology Facilities Council grant number ST/I000976/1.

REFERENCES