A self-consistent 3D model of fluctuations in the helium-ionizing background

TL;DR

This paper develops a 3D semi-analytic model of helium-ionizing background fluctuations that accounts for spatially varying mean free paths, successfully reproducing observed optical depth variations in the He II Ly$ ext{alpha}$ forest.

Contribution

It introduces a novel model that incorporates local ionization-dependent mean free paths, improving the understanding of large-scale fluctuations in the helium-ionizing background.

Findings

Model reproduces observed optical depth distributions

Fluctuations are more severe when mean free path varies spatially

Highlights importance of opacity fluctuations in IGM modeling

Abstract

Large variations in the effective optical depth of the He II Ly$\alpha$ forest have been observed at $z\gtrsim2.7$, but the physical nature of these variations is uncertain: either the Universe is still undergoing the process of He II reionization, or the Universe is highly ionized but the He II-ionizing background fluctuates significantly on large scales. In an effort to build upon our understanding of the latter scenario, we present a novel model for the evolution of ionizing background fluctuations. Previous models have assumed the mean free path of ionizing photons to be spatially uniform, ignoring the dependence of that scale on the local ionization state of the intergalactic medium (IGM). This assumption is reasonable when the mean free path is large compared to the average distance between the primary sources of He II-ionizing photons, $\gtrsim L_\star$ quasars. However, when this is no longer the case, the background fluctuations become more severe, and an accurate description of the average propagation of ionizing photons through the IGM requires additionally accounting for the fluctuations in opacity. We demonstrate the importance of this effect by constructing 3D semi-analytic models of the helium ionizing background from $z=2.5$-$3.5$ that explicitly include a spatially varying mean free path of ionizing photons. The resulting distribution of effective optical depths at large scales in the He II Ly$\alpha$ forest is very similar to the latest observations with HST/COS at $2.5 \lesssim z \lesssim 3.5$.

Figures

Equations

$$\frac{d\bar{\tau }}{dz}(\nu ,z)={\int }_0^{\infty }d{N}_{\mathrm{H}\mathrm{I}}f(N_{\mathrm{H}\mathrm{I}},z)(1-e^{-{\tau }_{\nu }}),$$

$${\tau }_{\nu }=N_{\mathrm{H}\mathrm{I}}{\sigma }_{\mathrm{H}\mathrm{I}}(\nu )+N_{\mathrm{He}\mathrm{I}}{\sigma }_{\mathrm{He}\mathrm{I}}(\nu )+N_{\mathrm{He}\mathrm{II}}{\sigma }_{\mathrm{He}\mathrm{II}}(\nu )$$

$${\eta }_{\mathrm{thin}}\equiv \frac{N_{\mathrm{He}\mathrm{II}}}{N_{\mathrm{H}\mathrm{I}}}=\frac{{\Gamma }_{\mathrm{H}\mathrm{I}}}{{\Gamma }_{\mathrm{He}\mathrm{II}}}\frac{{\alpha }_{\mathrm{He}\mathrm{II}}^A}{{\alpha }_{\mathrm{H}\mathrm{I}}^A}\frac{Y}{4X},$$

$$J_{\nu }(\vec{r},t)=\sum _i\frac{L_i(\nu )e^{-\tau (\vec{r},{\vec{r}}_i,\nu )}}{(4\pi |\vec{r}-{\vec{r}}_i|{)}^2},$$

$$\tau (\vec{r},{\vec{r}}_i,\nu )={\int }_{{\vec{r}}_i}^{\vec{r}}\frac{d\tau }{dz}\frac{dz}{dl}d{r}^{\prime },$$

$${\Gamma }_{\mathrm{He}\mathrm{II}}(\vec{r},t)=4\pi {\int }_{{\nu }_{\mathrm{He}\mathrm{II}}}^{\infty }\frac{J_{\nu }(\vec{r},t)}{h\nu }{\sigma }_{\mathrm{He}\mathrm{II}}(\nu )d\nu ,$$

Full text

1 Introduction

The hydrogen and helium reionization epochs produce some of the largest scale features in the Universe. They reflect the cumulative ionizing photon output of galaxies and black holes, or in other words, the history of star formation and supermassive black hole accretion across cosmic time. Helium reionization is believed to finish at $z\sim 3$–$4$, when quasars start to dominate the ionizing photon production in the universe (e.g. Sokasian et al. 2002; Wyithe & Loeb 2003; Furlanetto & Oh 2008). Occurring near the peak era of star and black hole formation, this is the last global transition experienced by the intergalactic medium (IGM) and significantly heats that material (e.g. Theuns et al. 2002; Hui & Haiman 2003; Bolton et al. 2009). This heating indirectly affects the star formation history of the Universe (by influencing its fuel supply) as well as observables like the H I  Ly$\alpha$ forest. The He II-ionizing background also offers the prospect of shedding light on the sources that drive it, including the abundance, spectra, lifetimes, and emission geometries of quasars (e.g. Worseck & Wisotzki 2006; Worseck et al. 2007; Furlanetto & Lidz 2011).

Observations of substantial transmission through the He II  Ly$\alpha$ forest below $z\sim 2.7$ (Davidsen et al., 1996; Kriss et al., 2001; Zheng et al., 2004; Fechner et al., 2006) suggest that helium reionization is complete by that epoch. At $z\ga 2.7$, however, the He II  Ly$\alpha$ forest opacity rapidly increases and significant fluctuations in the He II  Ly$\alpha$ effective optical depth $\tau_{\mathrm{eff}}$ arise on scales ranging from $10$–$100$ comoving Mpc (e.g. Reimers et al. 1997; Heap et al. 2000; Shull et al. 2004). Recent observations with HST/COS have confirmed these features (Shull et al., 2010; Worseck et al., 2011; Syphers & Shull, 2013, 2014) and continue to increase the number of He II  Ly$\alpha$ forest sightlines (Syphers et al., 2012; Worseck et al., 2016). Straightforward models of the He II  ionizing background with a uniform mean free path of ionizing photons cannot explain the fluctuations observed at $z\ga 2.7$ (Furlanetto & Dixon, 2010), suggesting that the culprit is incomplete helium reionization (see also Worseck et al. 2011). Interestingly, the H I  Ly$\alpha$ forest above $z\sim 5.5$ shows a similar increase in $\tau_{\mathrm{eff}}$ fluctuations after the epoch of hydrogen reionization has finished (Fan et al., 2006; Becker et al., 2015), which may be due to a spatially varying mean free path (Davies & Furlanetto 2016, but see D’Aloisio et al. 2015). The incomplete reionization interpretation also appears to be in moderate tension with the existence of regions with significant He II  transmission at $z\sim 3.5$ (Worseck et al., 2016; Compostella et al., 2014). Unfortunately, the era immediately following He II  reionization has not been well-studied, as radiative transfer simulations (which have mostly focused on the reionization process) are expensive and hence limited to “best-guess” parameters.

Davies & Furlanetto (2014, henceforth DF14) showed that fluctuations in the ionizing background could play an important role in the evolution of the mean background by inducing spatial variations in the opacity to ionizing photons. Toward higher redshifts, as the opacity of the IGM to He II-ionizing photons increases, the fluctuations in the background similarly increase (Furlanetto, 2009), leading to a rapidly evolving ionizing background. However, once these fluctuations become an important source of additional opacity the assumption of uniform emission and absorption of ionizing photons in commonly used 1D cosmological radiative transfer models (Haardt & Madau, 1996, 2012; Faucher-Giguère et al., 2009) breaks down. A full treatment of these effects requires a 3D realization of ionizing background-dependent opacity and discrete ionizing sources. In this work we construct such a model and find that it can describe the observed evolution of $\tau_{\mathrm{eff}}$ and its variations between sightlines decently well even under the assumption of a post-reionization universe.

The structure of the paper is as follows. In Section 2, we outline our prescription for the fluctuating opacity of the IGM to He II-ionizing photons. In Section 3, we describe our numerical model for computing the ionizing background and qualitatively discuss the effect of a fluctuating mean free path. In Section 4, we employ a simple prescription for the He II  Ly$\alpha$ forest to relate our ionizing background model to observations. In Section 5, we discuss the uncertainties of our simplified approach. Finally, in Section 6 we conclude by discussing the implications of our model and future improvements that will enhance its predictive power.

In this work we assume a standard $\Lambda$CDM cosmology with $H_{0}=70$ km/s/Mpc, $\Omega_{m}=0.3$, and $\Omega_{\Lambda}=0.7$. Distance units should be assumed to be comoving unless otherwise specified.

2 The Opacity of the IGM to Ionizing Photons

The opacity of the IGM to ionizing photons under the assumption of Poisson-distributed absorbing clouds can be parameterized by the optical depth per unit redshift,

[TABLE]

where $f(N_{\mathrm{H\,I}},z)$ is the H I  column density distribution function (CDDF), and

[TABLE]

is the optical depth of an absorber with H I  column density $N_{\mathrm{H\,I}}$ at frequency $\nu$. We adopt the shape and normalization of the $z\approx 2.5$ CDDF from Prochaska et al. (2014) and assume that the normalization evolves as $(1+z)^{2.5}$ to roughly match the measured redshift evolution of the H I  effective optical depth (Becker et al., 2013). For He II-ionizing photons with $\nu>4$ $\nu_{\mathrm{H\,I}}=\nu_{\mathrm{He\,II}}$, the opacity of an individual absorber is dominated by its He II  content, which depends on the ratio of the H I  and He II  ionization rates ($\Gamma_{\mathrm{H\,I}}$, $\Gamma_{\mathrm{He\,II}}$). In the optically thin limit, $N_{\mathrm{He\,II}}$ can be related to $N_{\mathrm{H\,I}}$ by the parameter $\eta_{\mathrm{thin}}$,

[TABLE]

where $\alpha_{\mathrm{H\,I}}^{A}$ and $\alpha_{\mathrm{He\,II}}^{A}$ are the case-A recombination coefficients of H I  and He II  , and $X$ and $Y$ are the hydrogen and helium mass fractions, respectively. For optically thick absorbers, the relationship is more complicated due to self-shielding of He II  and H I  ionizing photons. Modern cosmological radiative transfer models (e.g. Faucher-Giguère et al. 2009; Haardt & Madau 2012) determine $\eta(N_{\mathrm{H\,I}})$ by computing simplified radiative transfer of the ionizing background assuming a slab geometry and Jeans ansatz for absorbers (Schaye, 2001). Because of the implied dependence of absorber properties on $\Gamma_{\mathrm{He\,II}}$, ionizing background calculations must iterate several times to achieve self-consistency between the radiation field and the He II  absorber distribution. The link between $\Gamma_{\mathrm{He\,II}}$ and absorbers leads to an enhanced sensitivity of $\Gamma_{\mathrm{He\,II}}$ to evolution in the emissivity of ionizing photons (DF14; see also McQuinn et al. 2011).

In DF14, we extended this idea to link the local opacity of the IGM to the local $\Gamma_{\mathrm{He\,II}}$, which fluctuates significantly due to the rarity of the dominant sources of He II-ionizing photons (Furlanetto, 2009). In effect, this means that the mean free path of ionizing photons $\lambda$ fluctuates along with the intensity of the ionizing background. We showed that if ionizing photons sample the distribution of ionization rates $f(\Gamma)$, the overall opacity of the IGM increases due to the skewed nature of the distribution. However, naive application of this effect to a one-dimensional cosmological radiative transfer model caused the ionizing background to vanish at all redshifts unless a somewhat ad-hoc correction due to proximity effects was applied. Such a one-dimensional model also cannot account for spatial coherence of ionizing background fluctuations on large scales. These limitations suggest that a one-dimensional model is insufficient to study the effect of opacity fluctuations on the He II  ionizing background, motivating the three-dimensional approach described in the following section.

In this work, we apply the Faucher-Giguère et al. (2009) model for $\eta(N_{\mathrm{H\,I}})$ to compute $d\bar{\tau}/dz$ as a function of the local intensity of the He II-ionizing background. This procedure assumes that the He II  fraction in the absorbers responds instantaneously to changes in $\Gamma_{\mathrm{He\,II}}$. In reality, the He II  fraction will change on a characteristic timescale $t_{\mathrm{eq}}=(\Gamma_{\mathrm{He\,II}}+\alpha_{\mathrm{He\,II}}n_{e})^{-1}\sim 10(\Gamma_{\mathrm{He\,II}}/10^{-14.5})^{-1}$ Myr. If the average quasar lifetime is shorter than this timescale then non-equilibrium ionization effects could be very important (McQuinn, 2009). All other analytic treatments of the He II-ionizing background have made the same assumption with respect to the absorbers of He II-ionizing photons (e.g. Fardal et al. 1998; Faucher-Giguère et al. 2009; Haardt & Madau 2012).

3 A Numerical Model of the He II Ionizing Background

In this work, we present a simple three-dimensional extension to the 1D model of DF14. The basic structure of the model is as follows. Quasars are randomly placed in a cosmological volume 500 Mpc on a side from $z=4$ to $z=2.5$ following the Hopkins et al. (2007) $B$-band quasar luminosity function (QLF). By placing quasars randomly we neglect their clustering; the effect of quasar clustering on fluctuations in the ionizing background is likely small (Dixon et al., 2014) but could play a role when the mean free path is very short (Desjacques et al. 2014; see Section 5.1). The ionizing spectrum of each quasar is determined by first converting the $B$-band luminosity to the luminosity at the H I  ionizing edge with the constant conversion factor from Hopkins et al. (2007) and assuming that the spectrum at $\nu>\nu_{\mathrm{H\,I}}$ is a power law $L_{\nu}\propto\nu^{-\alpha_{\mathrm{Q}}}$ with $\alpha_{\mathrm{Q}}=1.6$ in agreement with Telfer et al. (2002) and consistent with the recent estimate by Lusso et al. (2015). We assume isotropic emission of ionizing radiation and a “lightbulb” model for quasar light curves with a lifetime of 50 Myr, similar to Compostella et al. (2013, 2014). The ionizing background due to these quasars is then calculated at all points on a 503 grid, giving a spatial resolution of 10 Mpc. This resolution is sufficient because the fluctuations in the ionizing background that we hope to characterize will only manifest on scales larger than the typical mean free path inside the simulation volume ($\ga 20$ Mpc), with the dominant scale set by the distance between bright sources ($\sim 40$ Mpc, cf. DF14). Additionally, the simplified absorber physics from Section 2 is likely only a reasonable description on large-scales, so higher resolution would not add to our understanding.

The specific intensity of ionizing radiation $J_{\nu}$ at each point $\vec{r}$ on the grid is computed by adding together the contribution from every quasar $i$ with specific luminosity $L_{i}(\nu)$ at position $\vec{r}_{i}$ that is turned on at time $t-|\vec{r}-\vec{r}_{i}|/c$,

[TABLE]

where the optical depth of ionizing radiation from the quasar $\tau(\vec{r},\vec{r}_{i},\nu)$, which in previous work has been approximated as $\sim|\vec{r}-\vec{r}_{i}|/\lambda_{\mathrm{mfp}}$, is computed by integrating the fluctuating IGM opacity along the light cone,

[TABLE]

where $d\tau/dz$ is given by equation (1) evaluated at the redshift corresponding to $t^{\prime}=t-|\vec{r}^{\prime}-\vec{r}|/c$, and which depends on the local ionization rate $\Gamma_{\mathrm{He\,II}}(\vec{r}^{\prime},t^{\prime})$ through the absorber model discussed in Section 2. The ionization rate is then computed by integrating over frequency,

[TABLE]

where in practice this integral is computed discretely with 12 logarithmic bins in $\nu$ from $\nu_{\mathrm{He\,II}}$ to $10^{1.2}\nu_{\mathrm{He\,II}}$. We neglect the redshifting of photons as they travel from the quasar to the cell and use the proper distance to the source rather than the luminosity distance when computing $J_{\nu}$, but these should be relatively small effects. Because the opacity between the cell and the quasar depends on the history of the opacity between the source and cell, a fully self-consistent model requires a time-dependent fluctuating background. We adopt a fiducial time step of 5 Myr, corresponding to $\sim 6$–$7$ Mpc of light travel distance per step.

Because the IGM opacity depends on the local value of $\Gamma_{\mathrm{He\,II}}$, which in turn depends on the opacity, as discussed in Section 2, the calculation must be iterated in order to achieve self-consistency. We iterate until the average change of $\Gamma_{\mathrm{He\,II}}(\vec{r},t)$ across all grid cells between iterations is less than 1 per cent at $z=3$, which typically requires 7–8 iterations.

We have additionally run a “control” simulation that ignores spatial fluctuations in the mean free path. Instead, $d\tau/dz$ is spatially uniform and computed as a function of redshift from a standard 1D ionizing background model with the same input parameters (i.e. quasar emissivity, CDDF). The rest of the computation is done in a similar manner to the full model described above, including finite quasar lifetimes and light cone effects. We will refer the this model as the “uniform MFP” model in the rest of the paper.

3.1 Results

The solid curves in Figure 1 show the evolution of the distribution of ionization rates, $f(\Gamma)$, from the simulation volume at $z=3.1$–$2.5$ in steps of $dz=0.2$. The dotted curves show the distributions for the uniform MFP model, which are nearly identical to analytic models for $f(\Gamma)$ computed using the same mean free path (Meiksin & White, 2003; Furlanetto, 2009). The distributions in the full simulation are broader than the uniform MFP model due to the enhancement of large-scale features in the radiation field by the fluctuating mean free path, and the disagreement increases with increasing redshift as the mean free path fluctuations become more pronounced. In the fiducial model, at $z=$ (2.6, 2.8, 3.0, 3.2) the median and 16–84th percentile ionization rates are $\log{\Gamma_{\mathrm{He\,II}}}=$ (-14.33${}^{+0.22}_{-0.13}$, -14.52${}^{+0.29}_{-0.22}$, -14.75${}^{+0.38}_{-0.35}$, -14.91${}^{+0.47}_{-0.60}$), corresponding to $\lambda=$ (73${}^{+25}_{-13}$, 45${}^{+22}_{-12}$, 26${}^{+19}_{-10}$, 17${}^{+17}_{-10}$) Mpc.

The evolution of fluctuations in the helium-ionizing background is shown visually in Figure 2 as a light cone projection from $z=3.5$–$2.5$. The parabolic features are caused by the intersection of the light cone with the expanding 50 Mly-thick “light shell” from a particularly bright quasar near the sightline (similar to Figure 5 in Croft 2004). At high redshift, the mean background is dominated by the relatively transparent proximity zones around luminous quasars. By $z\sim 2.5$, the mean free path is $\ga 100$ Mpc, leading to modest background fluctuations of about a factor of two that are similar to observations (McQuinn & Worseck, 2014). Figure 3 shows a series of snapshots of a slice through the simulation volume, with $\Gamma$ from the uniform and fluctuating mean free path models and maps of the logarithmic difference between the two. The dominant features change from transparent proximity zones at $z>3$ (top row) to large-scale ($\sim 200$ Mpc) coherent structures at $z\sim 2.9$ (middle row) before finally the ionizing background becomes mostly uniform and unaffected by mean free path fluctuations at $z<2.7$ (bottom row). At early times, the fluctuating model has substantial differences from the uniform model principally because the mean free path is substantially larger near the bright sources, leading to larger proximity zones (or “fossils” once the quasar has turned off), and deeper opacity far from those sources where the mean free path becomes very small. At later times ($z=2.87$ in Figure 3), the amplitude of these effects decreases, but their large-scale coherence remains, thanks to the rarity of the sources.

From Figure 1 we see that the high-$\Gamma$ end of $f(\Gamma)$ ($\log{\Gamma}\ga-13.8$) changes very little with redshift. These high-$\Gamma$ regions typically lie in the transparent proximity zones of bright quasars (DF14), which do not evolve significantly because they are largely decoupled from the global ionizing background evolution (i.e. the local mean free path is long enough that $\Gamma(r)\propto r^{-2}$ is a good approximation). The small variation with redshift that remains reflects fluctuations in the (small) number of extremely luminous quasars present in the simulation volume.

4 Predictions for the He II Ly$\alpha$ Forest

The primary observable of He II  is the opacity of the He II  Ly$\alpha$ forest ($\lambda_{\mathrm{rest}}<304$ Å) in far-ultraviolet spectra of $z\sim 3$–$4$ quasars. In this section we describe a simple model of mock He II  Ly$\alpha$ transmission spectra through the fluctuating ionizing background computed in the previous section.

The foundation of our mock spectrum model is a Monte Carlo realization of Poisson-distributed absorbers drawn from the observed CDDF, computed in a manner similar to Worseck & Prochaska (2011). We draw H I  absorbers in the range $10^{11.3}<N_{\mathrm{H\,I}}<10^{19}$ cm*-2* from the same evolving CDDF used in Sections 2 and 3 to compute the ionizing continuum opacity and place them randomly along a sightline. The He II  content of each absorber is computed as in Section 2, assuming a constant uniform $\Gamma_{\mathrm{H\,I}}=10^{-12}$ s*-1* (Becker & Bolton, 2013) and $\Gamma_{\mathrm{He\,II}}$ extracted along the light cone from the fluctuating He II  ionizing background simulation. Doppler widths for each absorber are drawn in the range $10<b<100$ km/s from the distribution function $dN/db\propto b^{-5}\exp{[-b_{\sigma}^{4}/b^{4}]}$ (Hui & Rutledge, 1999) with $b_{\sigma}=24$ km/s (Kim et al., 2001), and we assume that $b_{\mathrm{H\,I}}=b_{\mathrm{He\,II}}$ (i.e. the Doppler widths are dominated by nonthermal motion). Finally, we compute transmission spectra assuming Voigt profiles for each absorber using the efficient approximation of Tepper-García (2006).

An example He II  (black) and H I  (grey) transmission spectrum is shown in the top panel of Figure 4. The transmission through the He II  forest is closely tied to the variations in $\Gamma_{\mathrm{He\,II}}$ along the sightline, shown in the bottom panel. Large-scale regions along the light cone with enhanced $\Gamma_{\mathrm{He\,II}}$ are apparent, which share a characteristic width corresponding to our chosen quasar lifetime. Our He II  Ly$\alpha$ forest model has effectively infinite resolution and shows narrow regions with transmission close to unity, but more closely resembles observations when smoothed to mimic the resolution of HST/COS G140L grating, shown by the thick red curve.

We then binned mock spectra from 3000 randomly directed sightlines into $\Delta z=0.04$ pieces to allow direct comparison to the observations of Worseck et al. (2016). We compute the effective optical depth of each bin as $\tau_{\mathrm{eff}}=-\ln(\sum T_{i}/N)$ where $T_{i}$ is the transmission in each of the $N$ pixels inside the bin. The solid curve in Figure 5 shows the evolution of the median effective optical depth in the full simulation. The median optical depth increases steadily from $\tau_{\mathrm{eff}}\sim 2$ at $z\sim 2.7$ to $\tau_{\mathrm{eff}}\sim 5$ by $z\sim 3.1$, similar to the observations by Worseck et al. (2016). The dark and light shaded regions in Figure 5 show the 16–84th and 2.5-97.5th percentiles of $\tau_{\mathrm{eff}}$ values. Across the entire redshift range the low-$\tau$ end of the distribution, corresponding to regions with high $\Gamma_{\mathrm{He\,II}}$, evolves slowly. In contrast, the high-$\tau$ end, corresponding to regions with low $\Gamma_{\mathrm{He\,II}}$, increases rapidly above $z\sim 2.8$.

The thin dashed and dotted curves in Figure 5 show the median and distributions of $\tau_{\mathrm{eff}}$ in the uniform MFP simulation, re-normalized slightly to match the median optical depth in the fiducial model at $z=2.55$. The median $\tau_{\mathrm{eff}}$ evolution of the uniform MFP simulation is very similar to the full model, but the distribution tends to be more narrow. Fluctuations in the mean free path are most important in the tails of the distribution, as seen by the growing 2.5–97.5th percentile width at $z\ga 2.8$. DF14 predicted that the evolution of $\Gamma$ should be accelerated by including fluctuations in the mean free path because the IGM opacity increases when averaged over $f(\Gamma)$. Instead, in our 3D model the evolution of $\Gamma$ from the uniform MFP model is strikingly similar to the fluctuating MFP model at $z\la 3$. This lack of excess opacity from fluctuations is likely due to the fact that ionizing photons are preferentially emitted from bright quasars with transparent proximity zones, so the distribution of $\Gamma$ seen by ionizing photons is biased towards higher values.

The mock distribution is similar to the observations in Worseck et al. (2016) (red points), with the exception of a handful of high-$\tau_{\mathrm{eff}}$ regions at $z\sim 2.75$ and a very low $\tau_{\mathrm{eff}}$ region at $z=3.44$. In particular, the well-measured $\tau_{\mathrm{eff}}\sim 5$ region seen towards HE2347-4342 is quite rare in our simulations, appearing only in $\sim 0.1$ per cent of sightlines. There is some evidence for disagreement at $z\sim 3.4$ – our model predicts that only $\sim 1$–$2/10$ $\Delta z=0.04$ segments should show detectable transmission (i.e. $\tau_{\mathrm{eff}}\la 5.5$) while the observed fraction is $5/10$. We leave a full statistical comparison of our modeled optical depth distributions to observations, including forward modeling of instrumental noise, to future work.

The upper and lower bounds to the $\tau_{\mathrm{eff}}$ distribution are driven by the volume of space far away from (i.e. $\tau(\vec{r},\vec{r}_{\mathrm{Q}},\nu_{\mathrm{He\,II}})>1$) and very close to luminous quasars, respectively. The former evolves very quickly with redshift – not only is the number density of luminous quasars decreasing above $z\sim 2.5$, but more importantly the mean free path also decreases very rapidly. Additionally, as one moves further away from luminous sources in real space, the weaker ionizing background causes the mean free path to shrink, causing a non-linear increase in the optical depth. The lower bound to the $\tau_{\mathrm{eff}}$ distribution evolves only due to the evolution of the bright end of the QLF, with some fluctuations due to Poisson variance in the number of very luminous quasars – the “bumpy” evolution in our model is due to this variance rather than poor sampling of the $\tau_{\mathrm{eff}}$ distribution.

Our model allows a considerable amount of “tuning” of parameters to reproduce the smoothly evolving median $\tau_{\mathrm{eff}}$ at $z\la 3$, namely through the emissivity of He II-ionizing photons (via the quasar spectral index $\alpha_{\mathrm{Q}}$), the ionization state of H I  absorbers (via $\Gamma_{\mathrm{H\,I}}$), and the exact form of the CDDF (via the shape, normalization, and minimum $N_{\mathrm{H\,I}}$). It is actually a remarkable coincidence that – without deliberate tuning – our fiducial set of model parameters provides good agreement to the He II  $\tau_{\mathrm{eff}}$. We discuss the effect of variations in model parameters further in Section 5.

In Figure 6 we compare the associated H I  and He II  $\tau_{\mathrm{eff}}$ on scales of $\Delta z=0.01$ ($\approx 10$ Mpc) for comparison to Figure 8 of Worseck et al. (2016). This scale is also coincidentally the spatial resolution of our simulations, and thus the smallest scale that we can make predictions for variations in $\tau_{\mathrm{eff}}$ due to the fluctuating ionizing background. The behavior of the median optical depth is very similar to that of the Worseck et al. (2016) data at all redshifts. However, at $z\sim 2.85$ our simulations fail to reproduce the substantial observed fraction of regions with $\tau_{\mathrm{H\,I}}\la 0.5$ and $\tau_{\mathrm{He\,II}}\ga 6$, indicating that the Universe has stronger (rare) downward fluctuations in $\Gamma_{\mathrm{He\,II}}$ at this epoch than our simulations produce. In addition, one might expect weaker fluctuations due to the correlation between absorber systems and quasars, which we explicitly ignore and likely plays a stronger role in the distribution of effective optical depths at these scales compared to $\Delta z=0.04$. Nevertheless, the agreement of our median relationship between H I  and He II  suggests that our modeling approach is not unreasonable. It is worth noting that for $z\ga 3$ a large fraction of sightlines have little-to-no detectable flux given current observational limitations (Worseck et al., 2016), so judging the agreement between our simulation and observations is difficult.

5 Model Uncertainties

We have shown that a relatively simple model of fluctuations in the helium ionizing background can reproduce most of the observed properties of the He II  Ly$\alpha$ forest across cosmic time. This implies that fluctuations in the radiation field should be considered a viable alternative to ongoing He II  reionization (e.g. Worseck et al. 2011) as an explanation of the observed $\tau_{\mathrm{eff}}$ variations, at least for $z\la 3.2$. We discuss some important caveats to our simple model assumptions below.

5.1 Caveats

(i) He II  Ly$\alpha$ forest prescription

Our model assumes that the He II  Ly$\alpha$ forest can be described as an ensemble of randomly distributed absorbers following the Prochaska et al. (2014) CDDF. While this kind of description has been successfully used to model the H I  effective optical depth in previous work (e.g. Haardt & Madau 2012), He II  effective optical depths greater than $\sim 2.0$ at $z\la 3$ require near-complete blanketing of the spectrum by numerous low $N_{\mathrm{H\,I}}$ ($\la 10^{12}$ cm*-2*) lines (Fardal et al., 1998). The abundance of such “absorbers” is completely unconstrained by observations, and physical interpretation via the Jeans ansatz (Schaye, 2001) becomes unrealistic with implied absorber sizes of $\ga 2$ Mpc. In effect, these lines act as a “smooth” component to the He II  Ly$\alpha$ absorption that is required to reproduce the observed $\tau_{\mathrm{eff}}\ga 1.5$ absorption. Our treatment of the He II  Ly$\alpha$ forest as a set of discrete absorption lines is a rough approximation, with the low $N_{\mathrm{H\,I}}$ systems representing the smooth low density environments of cosmic voids that provide the majority of the He II  Ly$\alpha$ opacity (Croft et al., 1997).

To judge whether our He II  Ly$\alpha$ forest model responds to a fluctuating ionizing background similar to a more realistic density field (at least on large scales), we performed a simple comparison to skewers through a NYX Eulerian hydrodynamical simulation (Almgren et al., 2013; Lukić et al., 2015), kindly provided to us by Z. Lukić. The simulation is 100 $h^{-1}$ Mpc on a side with 40963 gas elements and dark matter particles, and we use a single output of the density and velocity fields at $z=3$ for our comparison. To calculate He II  Ly$\alpha$ transmission spectra of the skewers, we computed the equilibrium He II  fraction of each cell assuming a range of uniform $\Gamma_{\mathrm{He\,II}}$ and included the effects of thermal broadening and redshift-space distortions due to peculiar velocities. We found that the normalization, scaling of $\tau_{\mathrm{eff}}$ with $\Gamma_{\mathrm{He\,II}}$, and variations in $\Gamma_{\mathrm{He\,II}}$ on scales of $\Delta z=0.04$ was very similar to our model over the range of relevant values. In detail, the median $\tau_{\mathrm{eff}}(\Gamma_{\mathrm{He\,II}})$ relationship in the hydrodynamical simulation was slightly steeper, $\tau_{\mathrm{eff,hydro}}\propto\Gamma_{\mathrm{He\,II}}^{-0.54}$ vs. $\tau_{\mathrm{eff,CDDF}}\propto\Gamma_{\mathrm{He\,II}}^{-0.47}$, with a $\sim 10$ per cent lower normalization. These differences are not large enough to significantly change our results.

(ii) Ionization equilibrium

Our simulations assume that the universe is in ionization equilibrium – that is, reionization is assumed to have completed some time in the past. Assuming the gas density probability distribution from Miralda-Escudé et al. (2000) and following their procedure to compute the neutral fraction of the IGM, the equilibrium He II  fraction in our simulation is below $1$ per cent by volume at $z\la 3$. For the lowest $\Gamma_{\mathrm{He\,II}}$ regions at $z\sim 3$ the equilibrium He II  fraction can be as high as $\sim 10$ per cent, so an equilibrium treatment is unlikely to be very accurate. This is the primary weakness of our model compared to full radiative transfer simulations (e.g. McQuinn et al. 2009; Compostella et al. 2013, 2014).

(iii) Clustering of sources and absorbers

In both the source and absorber models we neglect the effect of clustering. As mentioned previously, Dixon et al. (2014) showed that the impact of quasar clustering on ionizing background fluctuations was likely small. Recent work by Desjacques et al. (2014) suggests that clustering could have a significant effect on background fluctuations at $z\sim 3$ if the mean free path is comparable to the correlation length of quasars, $r_{\xi}\sim 15$ Mpc. The average mean free path in our simulations is substantially larger than this for the redshifts we are interested in ($\lambda_{\mathrm{mfp}}\ga 25$ Mpc), but it is reasonable to expect that the addition of mean free path fluctuations would enhance the effect of clustering to some degree.

We ran an additional ionizing background simulation with clustered sources to test this possibility. Using the semi-numerical simulation code DEXM (Mesinger & Furlanetto, 2007), we generated a realistic distribution of dark matter halos in a (500 Mpc)3 volume at $z=3$. We then re-ran one of our ionizing background simulations with the quasars randomly assigned to locations of the most massive $N_{\mathrm{QSO}}$ halos, corresponding to $M_{\mathrm{h}}\ga 1.3\times 10^{12}$ M⊙, and approximating the halo distribution as constant with redshift. The median and low-$\tau_{\mathrm{eff}}$ end of the resulting opacity distribution were nearly identical to the fiducial, uniformly-distributed source model. At $z\la 3$, the upper-$\tau_{\mathrm{eff}}$ end was boosted by $\sim 10$ per cent and $\sim 30$ per cent at the 84th and 97.5th percentiles, respectively, suggesting that the clustering of sources can increase the probability of high-$\tau_{\mathrm{eff}}$ outliers. While the model with clustering more closely reproduces the observed distribution of $\tau_{\mathrm{eff}}$ at $z<3$, detailed modeling of the highly uncertain connection between quasars and dark matter halos is outside of the scope of this work.

In constructing our mock Ly$\alpha$ forest spectra the only large-scale fluctuations in the density field we consider are the Poisson variations in the number of absorber systems, which are expected to be fairly small (Fardal et al., 1998). Cosmic variance in the large-scale density field at the $\Delta z=0.04\approx 40$ Mpc scale considered in this work should also be small compared to the fluctuations in the ionizing background, with $\sigma(R=40\,\mathrm{Mpc})\sim 0.1$. While the overall effect of Ly$\alpha$ forest clustering may be small, its largest effect would likely be to extend the tails of the $\tau_{\mathrm{eff}}$ distribution, which could also ease some tension with the highest $\tau_{\mathrm{eff}}$ regions at $z\sim 2.8$ in addition to including source clustering as described above. The correlations between the density field and the source field (i.e. the positions of quasars) could in principle lead to smaller fluctuations in $\tau_{\mathrm{eff}}$, but this is unlikely to be a strong effect – the typical excess overdensity associated with a $\sim 10^{12-13}$ M⊙ dark matter halo on scales comparable to the $\sim 40$ Mpc observed He II  Ly$\alpha$ forest bins should be $<10$ per cent (see, e.g., Faucher-Giguère et al. 2008).

5.2 Cosmic Variance

Despite the large volume probed by our simulations, the mean optical depth and fluctuations that we measure are still affected somewhat by cosmic variance. The non-linear response of the ionizing background due to mean free path fluctuations implies that a random upward boost in the number of bright quasars leads to large-scale regions with larger-than-average mean free path that persist until the radiation from those quasars leaves the simulation volume. This can be seen in our fiducial simulation which has a pair of $\sim 70L_{*}$ quasars turn on at $z\sim 3.4$. The “boost” in $\Gamma_{\mathrm{He\,II}}$ and scatter of optical depths to low values from $z\sim 3.2$–$3.4$ seen in Figure 5 is almost entirely attributable to these two sources. This demonstrates the sensitivity of our model to the number of the brightest sources, at least at early times when the average mean free path is small.

In Figure 7 we show the $\tau_{\mathrm{eff}}$ distributions for a set of nine realizations of the (500 Mpc)3 ionizing background model (including the fiducial one). The solid curves show the mean value of the median, 16–84th percentiles, and 2.5–97.5th percentiles, while the shaded regions correspond to the minimum and maximum values in the full set of realizations. Below $z\sim 3$ the $\tau_{\mathrm{eff}}$ distributions between realizations agree to within $\pm 10$ per cent, but at higher redshift the range increases to $\ga 20$ per cent. This suggests that simulations which hope to accurately predict the post-reionization $\tau_{\mathrm{eff}}$ distribution should be at least the size of our model (i.e. $\ga 500$ Mpc on a side), but even larger volumes may be required for convergence at the $10$ per cent level. Achieving accurate radiative transfer in such large volumes will require substantial computational resources, far more than our simplified approach.

5.3 Variations in Model Parameters

The primary input parameters to our model are: the quasar luminosity function, the quasar ionizing spectral index $\alpha_{\mathrm{Q}}$, the quasar lifetime $t_{\mathrm{Q}}$, the strength of the H I-ionizing background $\Gamma_{\mathrm{H\,I}}$, and the column density distribution of H I  absorbers $f(N_{\mathrm{H\,I}},z)$. We leave a full analysis of variations in these parameters to future work, but we qualitatively describe their effects and importance below.

We have applied the Hopkins et al. (2007) QLF to populate the quasars in our simulation. More recent measurements by Glikman et al. (2011) suggest that there are substantially more faint quasars than the Hopkins et al. (2007) QLF at $z\sim 4$, which would result in a more uniform ionizing background. These measurements appear to be in conflict with the $z\sim 4$ measurements by Masters et al. (2012), but they instead find evidence for a somewhat steeper faint-end slope at $z\sim 3$ (see also Ross et al. 2013). Accurate model predictions for the helium-ionizing background will require resolution of this discrepancy, which we do not attempt in this work.

The primary effect of adjusting $\alpha_{\mathrm{Q}}$ is on the normalization of the quasar ionizing output at $\nu_{\mathrm{He\,II}}$, because $L_{\mathrm{He\,II}}=4^{-\alpha_{\mathrm{Q}}}L_{\mathrm{H\,I}}$, where $L_{\mathrm{He\,II}}$ and $L_{\mathrm{H\,I}}$ are the specific quasar luminosities at $\nu_{\mathrm{H\,I}}$ and $\nu_{\mathrm{He\,II}}$, respectively. At fixed redshift, due to the self-consistency of absorbers with the ionizing background, the average ionization rate should roughly follow $\Gamma_{\mathrm{He\,II}}\propto\epsilon^{1/(2-\beta)}$ (McQuinn et al., 2011), where $\beta$ is the power law index of the CDDF in absorbers with He II  ionizing opacity $\tau_{\nu,\mathrm{He\,II}}\la 1$ ($N_{\mathrm{H\,I}}\la 10^{15.5}$ cm*-2*). In the Prochaska et al. (2014) CDDF model that we adopt, $\beta=1.6$–$1.75$ in the relevant range of column densities, so we expect $\Gamma_{\mathrm{He\,II}}\propto\epsilon^{2.5\mathrm{-}4.0}\propto 4^{-(2.5\mathrm{-}4.0)\alpha_{\mathrm{Q}}}$. In practice, we find that $\Gamma_{\mathrm{He\,II}}\propto\epsilon^{2.5}$, so an adjustment of $\alpha_{\mathrm{Q}}$ by $0.1$ results in a change of $\Gamma_{\mathrm{He\,II}}$ by a factor of $\sim 1.4$. This sensitivity is completely degenerate with any other adjustment in the ionizing emissivity, such as the conflict between the typical conversion of the Hopkins et al. (2007) rest-frame optical QLF and the UV luminosity density measured by Cowie et al. (2009). Changing the shape of the He II-ionizing spectrum alone leads to a minor adjustment of the relationship between $\Gamma_{\mathrm{He\,II}}$ and the mean free path of average-energy ionizing photons that regulates the fluctuations in the background, so a harder spectrum would result in weaker fluctuations and a softer spectrum would lead to stronger fluctuations at fixed $\langle\Gamma_{\mathrm{He\,II}}\rangle$.

In our model we assume that quasars emit ionizing photons isotropically at a fixed rate for a fixed lifetime of 50 Myr and fixed spectral index $\alpha_{\mathrm{Q}}$. This is an enormously simplified picture of quasars that ignores realistic light curves (e.g. Hopkins & Hernquist 2009), anisotropic/beamed emission, and variations in the ionizing spectral index (Vanden Berk et al., 2001; Telfer et al., 2002). We tested simulations with quasar lifetimes of 25 and 100 Myr and found that the main results of this work, the distribution of $\tau_{\mathrm{eff}}$ over $dz=0.04$, were qualitatively unchanged. The primary difference was in the characteristic width of features in $\Gamma_{\mathrm{He\,II}}$ along the light cone, as seen in Figures 2 and 4. If quasar lifetimes are significantly shorter than the light travel time corresponding to $dz=0.04$, the blending of features would likely reduce the amount of fluctuations observed. However, on those short time scales ($\Delta{t}\la 10$ Myr), the equilibration time of the gas (Section 2) would become comparable to the quasar lifetime.

In our model, $\Gamma_{\mathrm{H\,I}}$ is fixed at a constant value, consistent with measurements by Becker & Bolton (2013). Because $\lambda_{\mathrm{mfp}}\propto\eta^{-1}\propto\Gamma_{\mathrm{H\,I}}$, adjusting $\Gamma_{\mathrm{H\,I}}$ is roughly equivalent to changing the normalization of the mean free path at fixed $\Gamma_{\mathrm{He\,II}}$. In effect, we find that $\Gamma_{\mathrm{He\,II}}\propto\Gamma_{\mathrm{H\,I}}^{-1}$, so for a given absorber in our He II  Ly$\alpha$ forest model, $N_{\mathrm{He\,II}}\propto\Gamma_{\mathrm{H\,I}}^{-2}$, and so very roughly $\tau_{\mathrm{eff}}\propto\Gamma_{\mathrm{H\,I}}^{-2}$. If one assumed a smaller $\Gamma_{\mathrm{H\,I}}$ (e.g. Faucher-Giguère et al. 2009), a simple way to maintain the same $\tau_{\mathrm{eff}}$ would be to increase the He II-ionizing emissivity such that $\eta$ remained constant, or in other words, $\Delta\log{\Gamma_{\mathrm{He\,II}}}\sim 2\times\Delta\log{\Gamma_{\mathrm{H\,I}}}$.

The redshift evolution of the mean free path and the details of its dependence on the local background are sensitive to the evolution and shape of the CDDF. The CDDF of H I  absorbers has been well-measured at $z\sim 2.5$, but discrepancies between different works exist in the difficult (i.e. saturated) $N_{\mathrm{H\,I}}$ regime that is most important to cosmological He II  radiative transfer ($N_{\mathrm{H\,I}}\sim 10^{15}$–$10^{16}$ cm*-2*; see, e.g., Kim et al. 2013; Rudie et al. 2013; Prochaska et al. 2014). While recent observations have greatly increased our knowledge of the CDDF at $z\sim 2.5$, the evolution of both the shape and normalization of the CDDF are still very uncertain. The effects of different CDDF parameterizations on the evolution of the He II-ionizing background are discussed in detail in DF14.

Given the number of uncertain parameters in modeling both the sources and absorbers, there are important degeneracies in our model. For example, for a fixed CDDF and model of the physical nature of absorbers, the measured average He II  $\tau_{\mathrm{eff}}$ corresponds to a locus of reasonable combinations of the He II-ionizing emissivity and $\Gamma_{\mathrm{H\,I}}$. In that sense, the matching of our fiducial model to the $z\sim 2.5$ measurements of the He II  $\tau_{\mathrm{eff}}$ is something of a coincidence – other “solutions” exist in a reasonable range of parameter space. However, our results show that a self-consistent model of the He II-ionizing background, with parameters consistent with other measurements, is compatible with the vast majority of the forest observations. This conclusion does not result from a fortuitous choice of parameters but is true for a large swath of - though by no means the entirety of - parameter space.

6 Conclusion

Recent observations of excursions to high effective optical depths in the He II  Ly$\alpha$ forest at $z\la 3$ have been interpreted as evidence for ongoing He II  reionization (Shull et al., 2010; Worseck et al., 2011). We have shown that the majority of the scatter in $\tau_{\mathrm{eff}}$ measurements at $z\la 3.2$ can be explained by fluctuations in the ionizing background when mean free path fluctuations are included self-consistently. Our model consists of a 3D realization of randomly distributed quasars following the measured quasar luminosity function, with finite quasar lifetimes and a finite speed of light. We additionally let the mean free path vary depending on the strength of the local ionizing background in a manner analogous to standard 1D cosmological radiative transfer models (Faucher-Giguère et al., 2009; Haardt & Madau, 2012) that assumes ionization equilibrium throughout the IGM.

The resulting radiation field fluctuates strongly on large scales, leading to large variations between sightlines. These large-scale features are due to the additional coherence caused by mean free path fluctuations, which are present even though we have neglected the clustering of sources in our fiducial simulations. If bright sources randomly lie close to each other, the excess background is enhanced due to the locally transparent IGM. Similarly, regions that are far from bright sources suffer from a more opaque IGM, further decreasing the radiation they receive. The strength of these effects increase strongly with redshift, such that by $z\ga 3.2$ they dominate the structure of the ionizing background (see the top panels of Figure 3). At these relatively high redshifts it is likely that the progression of He II  reionization is important to the budget of He II-ionizing photons and the size of proximity zones, but it is currently unclear how to distinguish the two scenarios observationally. Despite the enhanced fluctuations of our model relative to a uniform ionizing background, a handful of observed regions with high optical depth at $z\sim 2.7$–2.9 are in modest tension with our model predictions, suggesting that we may still be missing an important piece of the puzzle.

Recently, significant transmission in the He II  Ly$\alpha$ forest has been observed to persist to $z\sim 3.5$ (Worseck et al., 2016), which appears to be in (model-dependent) tension with late ($z\la 3$) He II  reionization (McQuinn et al., 2009; Compostella et al., 2013, 2014). At face value, our model appears inconsistent with these observations because the likelihood of encountering large-scale $\tau_{\mathrm{eff}}\la 4$ regions is small. At these higher redshifts He II  reionization is likely still ongoing, so our assumption of ionization equilibrium is no longer valid – however, it seems likely that this would limit the size of ionized regions and exacerbate the tension. Resolution of this discrepancy could come from the addition of clustering to the quasar distribution (although our simple test of this scenario suggests otherwise) or changes in the evolution and shape of the QLF.

The large-scale variations in the He II-ionizing background we have described here show that He II  reionization is a rich and complex event. The interaction of the sources, IGM, and radiation field require careful modeling, which is made possible by our detailed understanding of the IGM at $z\sim 3$. Observers are discovering more and more lines of sight along which the He II  Ly$\alpha$ forest can be studied, and the large fluctuations in our model imply that a full understanding of the reionization event will require exploring these new lines of sight in detail.

Acknowledgements

We would like to thank G. Worseck for valuable discussions, and for sharing his He II  effective optical depth data with us prior to publication. We would also like to thank Z. Lukić for providing us with NYX hydrodynamical simulation output. SRF was partially supported by NASA grant NNX15AK80G, administered through the ATP program, and by a Simons Fellowship in Theoretical Physics. SRF also thanks the Observatories of the Carnegie Institute of Washington for hospitality while this work was completed. Some of these simulations were run on computers in the UK (Apollo at the University of Sussex). KLD was supported by the Science and Technology Facilities Council [grant number STI/L000652/1].