Long troughs in the Lyman-$\alpha$ forest below redshift 6 due to islands of neutral hydrogen

TL;DR

This study uses high-resolution simulations to explain long Lyman-alpha forest troughs at redshift 5.5 as caused by neutral hydrogen islands, aligning with observed Ly$eta$ spikes and Ly$ extalpha$ emitter deficits.

Contribution

It demonstrates that late reionization scenarios can produce observed long troughs through neutral hydrogen islands, supported by detailed cosmological simulations.

Findings

Simulations reproduce observed long troughs at z~5.5.

Neutral hydrogen islands cause the troughs and allow Ly$eta$ transmission.

A deficit of Ly$ extalpha$ emitters is consistent with observations.

Abstract

A long (110 cMpc/$h$) and deep absorption trough in the Ly$\alpha$ forest has been observed extending down to redshift 5.5 in the spectrum of ULAS J0148+0600. Although no Ly$\alpha$ transmission is detected, Ly$\beta$ spikes are present which has led to claims that the gas along this trough must be ionized. Using high resolution cosmological radiative transfer simulations in large volumes, we show that in a scenario where reionization ends late ($z \sim 5.2$), our simulations can reproduce troughs as long as observed. In this model, we find that the troughs are caused by islands of neutral hydrogen. Small ionized holes within the neutral islands allow for the transmission of Ly$\beta$. We have also modelled the Ly$\alpha$ emitter population around the simulated troughs, and show that there is a deficit of Ly$\alpha$ emitters close to the trough as is observed.

Equations

$$P(\mathrm{REW}|M_{\mathrm{UV}})\propto \exp \left(\frac{-\mathrm{REW}}{{\mathrm{REW}}_{\mathrm{c}}(M_{\mathrm{UV}})}\right),$$

$${\mathrm{REW}}_{\mathrm{c}}(M_{\mathrm{UV}})=23+7(M_{\mathrm{UV}}+21.9)+6(z-4).$$

$${\mathrm{REW}}_{\min }=\{\begin{array}{ll}-20\mathring{\mathrm{A}}& M_{\mathrm{UV}}<-21.5,\\ 17.5\mathring{\mathrm{A}}& M_{\mathrm{UV}}>-19.0,\\ -20+6(M_{\mathrm{UV}}+21.5{)}^2\mathring{\mathrm{A}}& \mathrm{otherwise}.\end{array}$$

Full text

††pagerange: Long troughs in the Lyman-$\alpha$ forest below redshift 6 due to islands of neutral hydrogen–Long troughs in the Lyman-$\alpha$ forest below redshift 6 due to islands of neutral hydrogen††pubyear: 2019

1 Introduction

The Lyman-$\alpha$ (Ly$\alpha$) forest has long been recognised as a useful probe of the tail end of the epoch of reionization. Gunn-Peterson troughs observed towards redshift 6 are interpreted as a sign the intergalactic medium (IGM) is becoming increasingly neutral towards these redshifts. However, the Ly$\alpha$ forest saturates for relatively low H i fractions meaning that other methods must be used to place constraints on reionization. One such observation is the large scale spatial fluctuations in the opacity of the Ly$\alpha$ forest. These were first noted by Fan et al. (2006), however Lidz et al. (2006) claimed that the fluctuations could be explained by fluctuations in the density field alone. More recently, larger sets of higher quality observations that push to higher redshifts have been obtained (Becker et al., 2015; Bosman et al., 2018; Eilers et al., 2018). These observations show even larger fluctuations in the opacity of the Ly$\alpha$ forest at a given redshift. Furthermore, Becker et al. (2015) found a 110 cMpc/$h$ trough in the Ly$\alpha$ forest that showed complete absorption in Ly$\alpha$, although Ly$\beta$ transmission was observed. These observations could not be explained by variations in the density field and were claimed to be a signature of patchy hydrogen reionization.

Several competing models were proposed to explain how the H i fraction could vary over such large scales. Chardin et al. (2015, 2017) proposed a contribution to the UV background (UVB) from rare, bright sources such as AGN. Davies & Furlanetto (2016) suggested that the mean free path of ionizing photons may be very short in underdense regions, allowing for large scale fluctuations in the mean free path which in turn drive fluctuations in the amplitude of the UVB. D’Aloisio et al. (2015) suggested that the temperature dependence of the recombination rate may allow for a difference in H i fractions between regions that ionized early and have had time to cool compared with regions that ionized late and are still hot.

Although each of these models can potentially explain the distribution of Ly$\alpha$ forest opacities, they each require some key assumption which has been challenged. The model requiring rare, bright sources could lead to an early He ii reionization inconsistent with observations (e.g., D’Aloisio et al., 2017; Garaldi et al., 2019a). The fluctuating UVB model requires a mean free path at $z>5$ much shorter than observed at slightly lower redshift (Worseck et al., 2014). For the temperature fluctuations model, gas is required to reach very high temperatures as it is ionized which appear to be inconsistent with measurements of the IGM temperature above redshift 5 (e.g., Bolton et al., 2012; Walther et al., 2019).

There is however a way to discriminate between these models, as pointed out by Davies et al. (2018). In the case where the fluctuations are primarily in the gas temperature, the most overdense regions of the IGM will have the highest opacity as they were ionized first. In the case where the fluctuations in the amplitude of the UV background dominate, the high opacity regions will instead be the low density regions that are far from sources. This was tested by Becker et al. (2018), who used a narrowband filter on Subaru Hyper-Suprime Cam (HSC) to search for Ly$\alpha$ emitters (LAEs) in the vicinity of the long trough of ULAS J0148+0600. They found an underdensity of LAEs around the trough, suggesting that the sightline lay in a low density environment. This is consistent with the Davies & Furlanetto (2016) model, where a strong density dependence of the mean free path results in a rather low amplitude of the UVB in low density regions.

Kulkarni et al. (2019a) recently showed that another way to explain the fluctuations in the opacity of the Ly$\alpha$ forest is if reionization ended later than previously thought (see also Lidz et al., 2007; Mesinger, 2010). Using radiative transfer simulations of the reionization of the IGM, they found a reionization history that was consistent with Ly$\alpha$ forest constraints on the ionization state of the IGM below redshift 6 (McGreer et al., 2015; Bosman et al., 2018; Eilers et al., 2018) as well as constraints from the low optical depth to Thomson scattering measured with the CMB (Planck Collaboration, 2018). In this model, there were contributions from fluctuations in the UVB before reionization ends, as well as from temperature fluctuations in the IGM. However to reproduce the sightlines with the highest effective optical depths, it was crucial that residual islands of neutral hydrogen persisted down below redshift $z\leq 5.5$.

In this work we focus on explaining how the long trough of ULAS J0148+0600 can be explained in this late reionization model. In section 2, we present results from high resolution radiative transfer simulations that were tuned to match the mean flux of the Ly$\alpha$ forest below redshift 6. In section 3, we show that, as well as matching the distribution of mean fluxes measured along individual sightlines, these models can reproduce troughs as long as observed. In section 4, we model the distribution of LAEs around one of our simulated troughs and in section 5 we present a test to distinguish between a late reionization model and a scenario with a fluctuating UVB. Finally in section 6 we state our conclusions.

2 Modelling reionization

To model the properties of the IGM at the end of reionization, we present here a carefully calibrated simulation of hydrogen reionization. We ran mono-frequency cosmological radiative transfer simulations in post-processing using aton (Aubert & Teyssier, 2008, 2010), using a setup similar to Chardin et al. (2015). aton is a GPU-based code that solves the radiative transfer equation using a moment-based method with M1 closure for the Eddington tensor (e.g., Levermore, 1984; Gnedin & Abel, 2001; Aubert & Teyssier, 2008). We perform this post-processing on a set of density fields obtained from a hydrodynamic simulation run with Tree-PM SPH code p-gadget3 (last described by Springel, 2005). Our simulation has a box size 160 cMpc/$h$ and 2$\times$20483 gas and dark matter particles. The initial conditions were taken from the Sherwood simulation suite (Bolton et al., 2017) and assume cosmological parameters from Planck Collaboration XVI. (2014) with $\Omega_{\rm m}=0.308$, $\Omega_{\Lambda}=0.692$, $h=0.678$, $\Omega_{\rm b}=0.482$, $\sigma_{8}=0.829$ and $n=0.961$. The gas particles have mass $m_{\rm gas}=6.4\times 10^{6}$ $M_{\odot}$/$h$ and the dark matter particles have $m_{\rm DM}=3.4\times 10^{7}$ $M_{\odot}$/$h$. The gravitational softening length we assume is 3.1 ckpc/$h$. We run these simulations using a simplified star formation scheme, whereby any gas particles with a density greater than 1,000 times the mean cosmic baryon density and temperature less than 105 K are immediately converted into stars (Viel et al., 2004). This allows us to lower the time required to run the simulation without significantly affecting the properties of the intergalactic medium.

As we are performing the radiative transfer in post-processing, we neglect the hydrodynamical back-reaction of the gas to the photoheating. To compensate for this by adding some pressure smoothing to our input density fields, we also include photoheating by the Haardt & Madau (2012) uniform UV background in our hydrodynamic simulation. We start the simulations at redshift 99 and run them down to redshift 4. To allow us to update our input density field in the radiative transfer simulations frequently, we take snapshots of the simulation every 40 Myr starting just below redshift 20. We search for haloes in the simulation on the fly using a friends-of-friends algorithm. As aton requires the input density fields to be described on a cartesian grid, we map the SPH particles onto a 20483 grid (with cell size 78.125 kpc/$h$) using the SPH kernel. To assign sources in the simulation, we assume that each halo has an emissivity proportional to its mass (Iliev et al., 2006). The emissivity of each source is determined by summing the total mass of sources within our volume and dividing the assumed total emissivity within our volume among them, proportional to their masses. We take a minimum halo mass $10^{9}$ $M_{\odot}$/$h$. We find changing this minimum mass has hardly any effect on our results (once we have recalibrated the simulations to match the mean flux of the Ly$\alpha$ forest).

We calibrate the ionizing emissivity assumed in the simulation by running many simulations until the properties of our simulated Ly$\alpha$ forest (described later) agree well with the observations. As described above, the emissivity of each source is completely determined by the mass of the host halo and the volume emissivity we assume. When we calibrate the simulations, we only change the volume emissivity. This is equivalent to assuming some evolution in the properties of the sources, such as their escape fraction, but we do not explicitly treat these parameters. We make some changes with respect to the simulation presented in Kulkarni et al. (2019a), namely using a lower photon energy and a somewhat different emissivity evolution. The starting point for the emissivity evolution, which we modify as needed, is a galaxy-like component from Puchwein et al. (2019) and an increasing AGN-like contribution at lower redshifts. Although we do not model ionization by AGN directly (we only model one class of source, with a spectrum assumed to be galaxy-like), we find that we need an increasing volume emissivity at lower redshifts to match the mean flux of the Ly$\alpha$ forest (described further below). We do not try to separate our volume emissivity into contributions from galaxy-like and AGN-like components, but only account for the total emissivity that we require to reproduce the Ly$\alpha$ forest statistics. We therefore leave a more detailed comparison of the relative contributions of galaxies and AGN to the high-redshift UVB to future work. During reionization, we assume that all of our sources emit as 30,000 K blackbodies. This spectrum was chosen to provide reasonable agreement with measurements of the IGM temperature below redshift 6. This is equivalent to an average photon energy of 17.1 eV and a photoionization cross-section $3.9\times 10^{-18}$ cm*-2*. We arrive at these numbers assuming the optically thick limit, assuming all photons are absorbed locally (Pawlik & Schaye, 2011; Keating et al., 2018). This is likely a good approximation during reionization but will cause us to overestimate our photoheating after reionization has ended. For this spectrum, gas that is instantaneously ionized would reach a maximum temperature of approximately 13,000 K.

At lower redshifts, in a crude attempt to mimic the effects of the heating of the IGM by He ii reionization (e.g., Becker et al., 2011), we increase the photon energy in our simulation. We do not explicitly treat two classes of sources (galaxies and AGN) separately, but instead just try to capture the average effect of the evolution of the shape and amplitude of the UVB on the IGM. This change in the photon energy therefore applies to every photon in the simulation (both those emitted from sources and from recombinations). We begin increasing the energy at redshift 5.2 linearly, and by redshift 4 we use an average photon energy of 23.8 eV. For the intermediate redshifts we interpolate between the values at $z=5.2$ and 4 linearly with the cosmic expansion scale factor. Increasing the photon energy in this way is fully capturing the relevant physics, as any heating of the gas at these redshifts will come from recombinations, and we are not properly modelling e.g. the extra heating of the IGM by X-rays, so it is likely we will not recover an appropriate temperature-density relation as our voids will be under-heated. We also change the photon energy of photons that have already been emitted from sources, which will not conserve energy. We also neglect the corresponding change in the photoionization cross-section as we harden our spectra, as we wish to heat the gas without changing the photoionizing background. However, we found it was hard to recover the observed mean flux at redshift 4 if our gas was too cold and so decided to adopt this approximate treatment despite the shortcomings discussed above. We use the full speed of light in our simulations and take a timestep $\Delta t<\Delta x/c$ where $c$ is the speed of light and $\Delta x$ is the cell size of the grid. The chemical network only follows the photoionization, photoheating and cooling of hydrogen. The cooling processes we include are Hubble cooling, collisional ionization and excitation cooling (Cen, 1992), recombination cooling (Hui & Gnedin, 1997), Compton cooling (Haiman et al., 1996) and free-free cooling (Osterbrock & Ferland, 2006).

As we run the simulation, we extract the ionization state and temperature of the gas on the fly along a lightcone taken at a 30 degree angle through our grid. This allows us to account for the rapid evolution of the IGM along a line of sight, rather than having to interpolate between individual snapshots of the simulation as in Kulkarni et al. (2019a). We then construct mock Ly$\alpha$ forest spectra along sightlines along our output lightcone. The optical depth was computed using the approximate fit to Voigt profiles presented in Tepper-García (2006). We show the evolution of the mean flux in our simulations in Figure 1 compared with observations from Becker et al. (2015) and Bosman et al. (2018). In this work, we only focus on analysis of the Ly$\alpha$ forest down to $z\sim 5.5$ (the end-point of the Becker et al. (2015) long trough). However, as we also wish to look at the Ly$\beta$ transmission along our simulated troughs, it is important to get the mean flux of the Ly$\alpha$ forest close to the observed values at lower redshifts for adding in the contamination from the foreground Ly$\alpha$ forest.

We show the cumulative distribution function (CDF) of effective optical depths in Figure 2 compared with the distribution from Bosman et al. (2018). Here we use the conventional definition $\tau_{\rm eff}=-\log(\bar{F})$, where $\tau_{\rm eff}$ is the effective optical depth and $\bar{F}$ is the mean flux measured in 50 cMpc/$h$ segments of the Ly$\alpha$ forest. We emphasise that the simulations were calibrated to match the observed mean flux, and once this was reproduced the agreement between the observed and simulated $\tau_{\rm eff}$ CDFs came out naturally. Another compilation of effective optical depths is presented in Eilers et al. (2018), however due to systematic differences between the two datasets it is not possible to find a single reionization history that fits both measurements at once. We therefore limit our study here to a comparison with the Bosman et al. (2018) compilation. To measure the $\tau_{\rm eff}$ distributions, we create 500 realisations of the Bosman et al. (2018) sample and draw different sightlines from our sample in the same redshift distribution. We show the 15th and 85th percentiles of the different CDFs here. As in Kulkarni et al. (2019a), in this late reionization model we find good agreement with the observed sample, although we still fail to reproduce the high opacity end of the distribution in the $z=5.2$ redshift bin. This will be discussed later with reference to the incidence rate of long troughs in our simulations.

We present some further results from our simulation in Figure 3. The top left panel shows the assumed volume emissivity in our simulation compared with what was used in Kulkarni et al. (2019a). There are several reasons why we require a slightly different evolution to the assumed emissivity in that work. First, we assume a lower photon energy than in that work to try and recover temperatures that are closer to the Ly$\alpha$ forest measurements. Using this lower photon energy requires us to have a slightly faster reionization to recover the observed mean flux – since the gas does not reach as high a temperature as it is ionized, we require it to be ionized somewhat later so it has less time to cool down. Second, since we construct the spectra along a lightcone through the simulation, there are slight changes to the mean flux we recover from the simulations and therefore we require a different emissivity evolution. Third, as mentioned above, we also would like this simulation to match the mean flux of the Ly$\alpha$ forest below redshift 5. We find that this requires a rapid upturn of the ionizing emissivity just above redshift 5 (which could perhaps be attributed to the onset of an increasing contribution of AGN to the photoionizing background and the onset of He ii reionization, e.g. Haardt & Madau 2012; Puchwein et al. 2019; Kulkarni et al. 2019b). We further only measure the mean flux in spectra made along a lightcone taken in one direction through our simulation volume. If this exercise was repeated for a larger sample of sightlines, the required emissivity evolution would likely be somewhat different. Furthermore, as already mentioned, the different measurements of the mean flux in this redshift range by Bosman et al. (2018) and Eilers et al. (2018) are also not fully consistent within the quoted respective errors. We therefore caution against over-interpreting the emissivity evolution presented here.

Although a decline in the galaxy contribution to the UVB is predicted in models such as Haardt & Madau (2012) and Puchwein et al. (2019), it is difficult to find a physical interpretation for why the galaxy-like component should drop as rapidly as required here. It is perhaps also concerning that this rapid drop is required just as reionization ends in our simulations. However, the only other way to decrease the mean flux would be if the IGM was much colder, and it would be equally difficult to explain why it should cool so quickly since it needs to be reasonably warm to achieve the observed mean flux at redshift 6. The AGN-like component that we invoke to match the mean flux at lower redshifts becomes dominant over the galaxies much earlier than in other, more physically motivated models (e.g., Puchwein et al., 2019) and is in tension with the shape of the UV background inferred from metal lines (Boksenberg et al., 2003) and measurements of the ionizing emissivity of galaxies (Steidel et al., 2018) at somewhat lower redshifts. We leave a more careful investigation of the relation of the mean flux of the Ly$\alpha$ forest and the ionizing emissivity to future work, preferably using higher resolution simulations with more careful modelling of the gas temperature, and an improved measurement of the mean flux and the flux PDF.

The other panels of Figure 3 compare our simulation to other properties of the IGM. In this model, the volume is 50 per cent ionized at $z=6.7$ and 99.9 per cent ionized at $z=5.2$. In general, the agreement between the model and the observed IGM properties is good. However, as discussed above, since we not self-consistently modelling the QSO-like component of our emissivity, our temperatures below redshift 5 in the voids maybe be underestimated and our photoionization rate will be overestimated. We emphasise though that information from this redshift range is only included as the foregrounds to our Ly$\beta$ forest modelling and does not affect our main conclusions. The main differences between the model presented here and Kulkarni et al. (2019a) are driven by the temperature of the gas. Our choice of a lower photon energy here brings the temperature at mean density in the simulations closer to the observations, however as discussed above this is at the expense of introducing a rather sharp drop in the evolution of the emissivity at $z\sim 6$. To achieve the same mean flux, the lower temperatures require a higher photoionization rate. This also seems to be closer to the observed values. We also show for comparison the photoionization rate from the Haardt & Madau (2012) model, as this was used when running the underlying hydrodynamic simulation. We caution though that the published temperature measurements rely heavily on calibration with numerical simulations that assume a homogeneous UV background, and do not account for the effects of inhomogeneous reionization and spatial fluctuations in the temperature-density relation (Furlanetto & Oh, 2009; Lidz & Malloy, 2014; D’Aloisio et al., 2015; Keating et al., 2018) which is likely to bias the published measurements low.

3 Long troughs

We next measure the longest trough along 250,000 sightlines along our extracted lightcone. We post-process the spectra to match the observations, by rebinning onto coarser pixels, convolving with an appropriate instrument profile and adding noise. We measure the length of the troughs by measuring out to a maximum redshift $z_{\rm max}=5.879$ (the high redshift point of the Becker et al. 2015 trough). We then search for spikes in transmission below this redshift which we assume are regions of the spectrum that show transmission greater than the 1$\sigma$ noise level for at least four consecutive pixels and that have a combined significance of more than 5$\sigma$. We then define the regions in between these transmission spikes as the troughs, and record the longest trough measured for each sightline. As shown in the left panel of Figure 4, our simulation does produce troughs in the Ly$\alpha$ forest as long as the 110 cMpc/$h$ trough of ULAS J0148+0600, however these are very rare in our simulation (approximately 1 in 10,000). In a recent paper, Giri et al. (2019) also presented an investigation of neutral islands at the end of reionization, and noted that they find islands long enough to reproduce the long trough of ULAS J0148+0600. Their simulation volume is much larger (714 cMpc vs. 160 cMpc/$h$ in this work) and is therefore more suited for studying the statistics of these long troughs at the end of reionization. Detailed comparisons with observations will however also require that the simulations are calibrated to match the properties of the Ly$\alpha$ forest, such as in this work.

The long troughs we do find are as dark as measured by Becker et al. (2015), who estimated that the effective optical depth was greater than 7.2 along the absorption trough of ULAS J0148+0600 (middle panel of Figure 4). We find that troughs with length up to $\sim$ 80 cMpc/$h$ have a relatively high incidence rate in the simulation, but the incidence rate falls rapidly above this point. The most natural explanation for this is that reionization is still ending too early in our simulations. As discussed above, there are 50 cMpc/$h$ segments of the Ly$\alpha$ forest that are observed to be totally opaque at $z=5.2$ which our simulation does not reproduce. If reionization ended later, then the neutral islands would be longer and so would the associated troughs. We demonstrate this by moving the redshift at which we begin measuring the troughs at to a higher redshift ($\Delta z=0.1$). As shown in the left panel of Figure 4, this has the effect of increasing the length of our simulated troughs and we now find many more longer than 110 cMpc/$h$. However, we note that we were not successful in finding a model where reionization ended later while still matching the observed mean flux.

Another point worth making is that the size of our simulation volume is still small in the context of reionization studies (e.g., Iliev et al., 2014). We therefore do not have a large sample of neutral islands at low redshift in our simulation. This is shown in the right panel of Figure 4 – all of the long troughs are clustered in one section of the box. This means that the deficit of long troughs in our simulations could be caused by an unfavourable geometry in these remaining neutral islands, i.e. they may just happen to not be very extended in the direction in which we extract our lightcone. A stronger statement about the incidence rate of long troughs for a given ionization history will therefore require larger volumes while still maintaining resolution comparable to the simulation presented here. We therefore leave a more detailed study of the incidence rate of these long troughs for future work.

An example of a long (98.2 cMpc/$h$) trough is shown in Figure 5, as well as the associated physical quantities of the gas along the trough. We find that the region spanned by the trough is coincident with an island of neutral hydrogen. It is clear from this figure that the shape of the neutral islands is directly linked to the trough length, and one could imagine that a sightline taken at a small angle to the line of sight shown here would produce a longer trough. As seen in Figure 5, the UVB is also lower in the vicinity of the neutral island, as it is only being illuminated by sources in one direction. This means that the completely neutral gas is surrounded by gas with an H i fraction of about 10*-3*, still neutral enough to saturate the Ly$\alpha$ transmission. We also find Ly$\beta$ transmission along the trough, as is observed. This can be explained by small pockets of lightly ionized gas within this neutral island. The photoionization rate inside these pockets is still low, however the gas is ionized enough to allow transmission of Ly$\beta$, since Ly$\beta$ has an oscillator strength $\sim 5$ times lower that Ly$\alpha$.

4 Lyman-$\alpha$ emitters

A convincing model of the long trough in ULAS J0148+0600 must also explain the observed deficit of LAEs around that sightline. Since we construct our spectra along a planar lightcone taken at an angle wrapped through our simulation, we cannot make use of periodic boundary conditions to centre a simulated trough in the middle of the lightcone. We therefore pick the sightline we are interested in (in this case, the sightline shown in Figure 5) and rerun the radiative transfer simulation to extract a new lightcone centred on this sightline. This means that our analysis of the LAEs around the trough is restricted to a single sightline. However, since the size of our volume is close to the field of view of HSC and also our longest troughs are co-spatial (Figure 4) an analysis of more than one trough in the simulation would not really remove the effects of cosmic variance. We chose this trough because, although it is slightly shorter than the trough of ULAS J0148+0600, it showed Ly$\beta$ emission with similar properties to the observed long trough, which shows only one Ly$\beta$ transmission spike in the redshift range covered by the narrowband filter used in Becker et al. (2018).

To model the LAEs around our simulated trough, we follow the method described in detail in Weinberger et al. (2019). We note that this method is similar to the LAE modelling in Davies et al. (2018) and Becker et al. (2018), although as we will describe later we also take the effect of the local IGM opacity into account. We first assign UV luminosities to the dark matter haloes by abundance matching to the UV luminosity function at redshift 5.9 (Bouwens et al., 2015). We assume a 50 Myr duty cycle such that only a fraction of our galaxies are UV bright at a given time, reflecting the bursty nature of star formation (Trenti et al., 2010). Following Dijkstra & Wyithe (2012), we then assign a rest-frame equivalent width of Ly$\alpha$ emission to each halo based on its UV luminosity using the relation

[TABLE]

where REW is the rest-frame equivalent width, $M_{\rm UV}$ is the UV magnitude and REWc is defined as

[TABLE]

The normalization is defined such that the population has ${\rm REW}_{\rm min}\leq{\rm REW}\leq{\rm REW}_{\rm max}$, where we choose ${\rm REW}_{\rm max}=300\,{\rm\AA}$. The minimum REW is defined in terms of $M_{\rm UV}$ as

[TABLE]

We note that the assignment of UV magnitude and REW is a random process, since only a fraction of the haloes are UV bright at a given time and the REW of each halo is drawn from a distribution. We therefore repeated this process 50 times, to generate 50 different realisations of LAE populations.

We next accounted for the attenuation of Ly$\alpha$ emission due to neutral gas around the haloes. The default resolution of our radiative transfer simulations is too coarse to resolve the circumgalactic medium (CGM) of halos, so we re-extract a higher resolution skewer of gas density and peculiar velocity in a 10 cMpc/$h$ radius around each halo. The resolution of these skewers is 9.7 ckpc/$h$. We insert this short high resolution sightline into the longer lower resolution sightline to model the attenuation of Ly$\alpha$ emission due to both the CGM and IGM. Since we do not know the ionization state of the gas in the high resolution sightline, we assume that the gas is in ionization equilibrium calculated using the photoionization rates from our radiative transfer simulation plus the self-shielding prescription of Chardin et al. (2018b). We then calculate the optical depth along each sightline.

We construct mock LAE spectra by using the UV magnitude to get the flux at rest-frame 1600 Å, and assuming that $F_{\lambda}\propto\lambda^{-2}$. To model the Ly$\alpha$ emission, we add a gaussian emission line with a REW given by the prescription outlined above. As in Weinberger et al. (2019), we also assume some velocity shift between the redshift of the host halo and the redshift of the Ly$\alpha$ emission, which accounts for the complex radiative transfer of the Ly$\alpha$ line out of the halo. We take this to be a constant factor times the circular velocity of the host halo, with the constant factor chosen to ensure good agreement with the rest-frame equivalent widths and luminosities of observed LAEs. We then attenuate the spectrum by multiplying by the normalised Ly$\alpha$ transmission we calculate along each sightline. In this work we found that a velocity shift of $1.4v_{\rm circ}$ was a good choice. Using these spectra, we then measure the magnitude of the LAEs in the $i2$ and $NB816$ bands and select LAEs that have $NB816\leq 26.0$ and $i2-NB816\geq 1.2$ as in Becker et al. (2018). We note that Becker et al. (2018) also use $r2$ band magnitudes as part of their selection process, but we can not model that band here because our spectra do not extend down to sufficiently low redshift. We calculate the luminosity and rest-frame equivalent width for each LAE by integrating over our constructed LAE spectra. As in Weinberger et al. (2019), the agreement between the observed and simulated LAEs is very good, i.e. our modelled LAE population agrees with the observed LAE luminosity function and rest-frame equivalent width distribution from the SILVERRUSH survey (Figure 6), however we note that our box is too small to reproduce the brightest LAEs (seen by the rapid cutoff of our modelled luminosity function).

We next compare the spatial distribution of LAEs around our simulated trough (right panel of Figure 6). We show here the projected position of LAEs coloured by their $NB816$ magnitude. As in Becker et al. (2018), we find a large deficit of LAEs in the immediate vicinity of the trough. This is to be expected, since the low density regions will be the last to ionize. We therefore expect our remaining neutral islands to trace underdense regions of the IGM. Furthermore, since the Ly$\alpha$ optical depth along the trough is high, we also expect the attenuation of Ly$\alpha$ emission to be stronger in the region around the trough. In the left panel of Figure 7, we show the LAE surface density computed in annuli centred around the trough. The results from our simulated LAE are shown in blue (again, the shaded regions mark the 15th/85th percentiles of our 50 different LAE realisations and the solid line is the median value). The results from Becker et al. (2018) are shown in black. Similar to the observations, we find that the surface density of LAEs begins to decline rapidly in the inner 30 cMpc/$h$ around the trough. This is predominantly due to the trough passing through an underdense region, although we note that that this effect is complemented by a decrease in the IGM transmission fraction of about 10 per cent close to the trough (middle panel of Figure 7). Here we define the IGM transmission fraction as $T_{\rm IGM}=L_{\rm obs}$/$L_{\rm int}$ where $L$ is the observed/intrinsic LAE luminosity (Mesinger et al., 2015). This lower transmission fraction is in agreement with Sadoun et al. (2017), who suggest that the observed number of LAEs will be lower in regions with a lower photoionizing background, due to the increased filling factor of self-shielded systems in the CGM of the host galaxies. However, it does not appear to be the dominant effect here.

In the right panel of Figure 7, we show the surface density of UV-selected galaxies that have $M_{\rm UV}<-20$, $-20.5$ and $-21$. For the cut with $M_{\rm UV}<-20$, we find that the surface density far from the trough is well matched to the LAEs. However, the decline in surface density is even stronger for this selection, already beginning to decline at a radius 50 cMpc/$h$ from the trough and highlighting that this sightline passes through an underdense region. The median mass of the LAEs that pass the colour cuts used here is $4.1\times 10^{10}\,M_{\odot}$/$h$. For the UV selected galaxies, the median masses are $(1.1,1.7,3.1)\times 10^{11}\,M_{\odot}$/$h$ for $M_{\rm UV}<(-20,-20.5,-21)$.

5 Comparison with other work

As shown in Becker et al. (2018), the fluctuating UVB model presented in Davies & Furlanetto (2016) can also reproduce the deficit of LAEs in high $\tau_{\rm eff}$ regions. In this case, a short mean free path leads to large fluctuations in the UVB. This produces regions that, while still ionized, have a H i fraction high enough to saturate transmission of Ly$\alpha$. We do not find that this is the case in our simulations, as once the ionized regions have percolated the mean free path becomes long and any fluctuations in the UVB are quickly washed away. Our late reionization model is still similar to the Davies & Furlanetto (2016) in that it is the underdense regions that produce the high $\tau_{\rm eff}$ sightlines. In our case, however, this is because reionization has not yet ended and the voids host the residual neutral islands. In both models, the underdense regions host the high $\tau_{\rm eff}$ sightlines because they are far from the regions where the space density of ionizing sources is highest. Linking the Davies & Furlanetto (2016) model to a late reionization scenario removes the main issues with that explanation. First, that the mean free path needs to become much shorter at $z>5$ than observed by Worseck et al. (2014) at redshift 5. In our model, the mean free path does decrease sharply above redshift 5 due to the residual islands of neutral hydrogen. Second, there is no obvious way to predict the redshift evolution of the mean free path in the Davies & Furlanetto (2016) model – at each redshift it must be tuned to match constraints from the Ly$\alpha$ forest. We also tune the emissivity in our model to match the observed mean flux, however once this is fixed then the reionization history of our preferred model can be used to make predictions for higher redshifts. In this sense, our late reionization model can be thought of as a dynamic version of the Davies & Furlanetto (2016) model.

Ideally one would still like to discriminate between the two models though, to answer if indeed the high $\tau_{\rm eff}$ regions are mildly ionized or completely neutral. One difference between the two scenarios is that the fluctuating UVB model predicts that the underdense regions will always coincide with high $\tau_{\rm eff}$ sightlines, and the low $\tau_{\rm eff}$ sightlines will lie in overdense regions (Davies et al., 2018). In a late reionization model though, the underdense regions will align with high $\tau_{\rm eff}$ sightlines only in the time before they have been ionized. Once the neutral islands have vanished, the voids should instead correspond to the lowest $\tau_{\rm eff}$ regions, as most of the Ly$\alpha$ forest transmission is dominated by the voids. These recently ionized regions should also be very hot, suppressing recombinations and further enhancing transmission (D’Aloisio et al., 2015).

We present an example of two such sightlines in Figure 8. The left panel shows a map of $\tau_{\rm eff}$ measured in 50 cMpc/$h$ segments centred on the redshift corresponding to peak transmission of the $NB816$ filter. The blue and red circles mark two regions of high and low $\tau_{\rm eff}$. In the middle panel, we show the mean surface density of LAEs in that region. The surface density is taken to be inversely proportional to the square of the distance to the tenth nearest neighbour and the map has been smoothed with a gaussian kernel. The red and blue circles correspond to the same regions as in the left panel. It is clear that both of these regions show a deficit of LAEs, suggesting that both of these are underdense regions. This is shown explicitly in the right panel, where we compare the surface density of LAEs as a function of radius for the two lines of sight. In the Davies & Furlanetto (2016) model, you would instead expect that the low $\tau_{\rm eff}$ region should show an overdensity of LAEs. This suggests that future observations searching for LAEs around sightlines that show significant Ly$\alpha$ forest transmission may be interesting for differentiating between these two models.

6 Summary and Conclusions

We have shown here that long Gunn-Peterson troughs observed in the Ly$\alpha$ forest above redshift 5 can be explained by islands of neutral hydrogen in the late reionization scenario presented in Kulkarni et al. (2019a). We have demonstrated that this model also explains the significant deficit of LAEs detected around these high opacity sightlines by Becker et al. (2018). While these observations can also be explained by a model in which the amplitude of the UVB is low in underdense regions, we argue that the model presented here is more consistent with the observed evolution of the mean free path for ionizing photons. We further present predictions for the distribution of LAEs around the low opacity sightlines which should allow future observations to discriminate between these two scenarios. This however will still not be a definitive test of neutral gas in the IGM, so other tests should also be pursued (Malloy & Lidz, 2015).

We note that the reionization model presented here still does not reproduce the most opaque sightlines in the range $5.1<z<5.3$ within the 15th/85th interval. Searching for more saturated sightlines at low redshifts would be useful in constraining whether these sightlines are rare. If they turn out not to be outliers, this may call for a reionization history that ends even later than shown here. This would also likely increase the incidence rate of long ($>$ 100 cMpc/$h$) troughs. We note also that the temperature of the IGM is still a significant uncertainty in our models, as it can have a large effect on the mean flux (as the recombination rate and thus the effective optical depth for a given UVB amplitude depend rather strongly on temperature). More accurate models of the ionization history of the IGM would benefit from more estimates of the temperature at high redshift, perhaps also using methods that account for the pressure effects for spatial fluctuations of the temperature-density relation due to reionization (e.g., Oñorbe et al., 2019, Puchwein et al. in prep).

As more high signal-to-noise Ly$\alpha$ forest sightlines are accumulated and measurements of the mean flux and flux PDF further improve, there will be benefit in moving from studies of integrated quantities such as the effective optical depth to more detailed studies of individual features. Future work modelling the occurrence rate and distribution of these troughs, as well as the transmission peaks between them, should provide more detailed constraints on the latter half of reionization (Chardin et al., 2018a; Garaldi et al., 2019b). Tying these Ly$\alpha$ forest measurements to the galaxy population may also yield additional information such as the sources that contribute to reionization (Kakiichi et al., 2018), especially once JWST is operational. Studies utilising also the Ly$\beta$ forest should allow to push constraints on the reionization history outside of the immediate vicinity of high redshift quasars to $z\sim 6.5$ and beyond.

Acknowledgements

We thank George Becker for useful discussions and for the suggestion that motivated Section 5. This work was performed using resources provided by the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service (www.csd3.cam.ac.uk), provided by Dell EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences Research Council (capital grant EP/P020259/1), and DiRAC funding from the Science and Technology Facilities Council (www.dirac.ac.uk). This work used the COSMA Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility. This equipment was funded by a BIS National E-infrastructure capital grant ST/K00042X/1, DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E- Infrastructure. This work was supported by the ERC Advanced Grant 320596 “The Emergence of Structure During the Epoch of Reionization”.