Long troughs in the Lyman- α forest below redshift 6 due to islands of neutral hydrogen
Laura C. Keating, Lewis H. Weinberger, Girish Kulkarni, Martin G. Haehnelt, Jonathan Chardin, Dominique Aubert
MMNRAS , 1–11 (2019) Preprint 31 October 2019 Compiled using MNRAS L A TEX style file v3.0
Long troughs in the Lyman- α forest below redshift 6 dueto islands of neutral hydrogen Laura C. Keating (cid:63) , Lewis H. Weinberger , , Girish Kulkarni ,Martin G. Haehnelt , , Jonathan Chardin and Dominique Aubert Canadian Institute for Theoretical Astrophysics, 60 St. George Street, University of Toronto, ON M5S 3H8, Canada Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK Department of Theoretical Physics, Tata Institute of Fundamental Research, Homi Bhabha Road, Mumbai 400005, India Observatoire Astronomique de Strasbourg, 11 rue de l’Universite, 67000 Strasbourg, France
31 October 2019
ABSTRACT
A long (110 cMpc/ h ) and deep absorption trough in the Ly α forest has been observedextending down to redshift 5.5 in the spectrum of ULAS J0148+0600. Although noLy α transmission is detected, Ly β spikes are present which has led to claims that thegas along this trough must be ionized. Using high resolution cosmological radiativetransfer simulations in large volumes, we show that in a scenario where reionizationends late ( z ∼ . β . We havealso modelled the Ly α emitter population around the simulated troughs, and showthat there is a deficit of Ly α emitters close to the trough as is observed. Key words: galaxies: high-redshift – quasars: absorption lines – intergalactic medium– methods: numerical – dark ages, reionization, first stars
The Lyman- α (Ly α ) forest has long been recognised as a use-ful probe of the tail end of the epoch of reionization. Gunn-Peterson troughs observed towards redshift 6 are interpretedas a sign the intergalactic medium (IGM) is becoming in-creasingly neutral towards these redshifts. However, the Ly α forest saturates for relatively low H i fractions meaning thatother methods must be used to place constraints on reion-ization. One such observation is the large scale spatial fluc-tuations in the opacity of the Ly α forest. These were firstnoted by Fan et al. (2006), however Lidz et al. (2006) claimedthat the fluctuations could be explained by fluctuations inthe density field alone. More recently, larger sets of higherquality observations that push to higher redshifts have beenobtained (Becker et al. 2015; Bosman et al. 2018; Eilers et al.2018). These observations show even larger fluctuations inthe opacity of the Ly α forest at a given redshift. Further-more, Becker et al. (2015) found a 110 cMpc/ h trough inthe Ly α forest that showed complete absorption in Ly α , al-though Ly β transmission was observed. These observationscould not be explained by variations in the density field and (cid:63) E-mail: [email protected] were claimed to be a signature of patchy hydrogen reioniza-tion.Several competing models were proposed to explain howthe H i fraction could vary over such large scales. Chardinet al. (2015, 2017) proposed a contribution to the UV back-ground (UVB) from rare, bright sources such as AGN.Davies & Furlanetto (2016) suggested that the mean freepath of ionizing photons may be very short in underdenseregions, allowing for large scale fluctuations in the mean freepath which in turn drive fluctuations in the amplitude of theUVB. D’Aloisio et al. (2015) suggested that the temperaturedependence of the recombination rate may allow for a differ-ence in H i fractions between regions that ionized early andhave had time to cool compared with regions that ionizedlate and are still hot.Although each of these models can potentially explainthe distribution of Ly α forest opacities, they each requiresome key assumption which has been challenged. The modelrequiring rare, bright sources could lead to an early He ii reionization inconsistent with observations (e.g., D’Aloisioet al. 2017; Garaldi et al. 2019a). The fluctuating UVBmodel requires a mean free path at z > c (cid:13) a r X i v : . [ a s t r o - ph . C O ] O c t L.C. Keating et al. reach very high temperatures as it is ionized which appearto be inconsistent with measurements of the IGM tempera-ture above redshift 5 (e.g., Bolton et al. 2012; Walther et al.2019).There is however a way to discriminate between thesemodels, as pointed out by Davies et al. (2018). In the casewhere the fluctuations are primarily in the gas temperature,the most overdense regions of the IGM will have the highestopacity as they were ionized first. In the case where the fluc-tuations in the amplitude of the UV background dominate,the high opacity regions will instead be the low density re-gions that are far from sources. This was tested by Beckeret al. (2018), who used a narrowband filter on Subaru Hyper-Suprime Cam (HSC) to search for Ly α emitters (LAEs) inthe vicinity of the long trough of ULAS J0148+0600. Theyfound an underdensity of LAEs around the trough, suggest-ing that the sightline lay in a low density environment. Thisis consistent with the Davies & Furlanetto (2016) model,where a strong density dependence of the mean free pathresults in a rather low amplitude of the UVB in low densityregions.Kulkarni et al. (2019a) recently showed that anotherway to explain the fluctuations in the opacity of the Ly α forest is if reionization ended later than previously thought(see also Lidz et al. 2007; Mesinger 2010). Using radiativetransfer simulations of the reionization of the IGM, theyfound a reionization history that was consistent with Ly α forest constraints on the ionization state of the IGM be-low redshift 6 (McGreer et al. 2015; Bosman et al. 2018;Eilers et al. 2018) as well as constraints from the low opti-cal depth to Thomson scattering measured with the CMB(Planck Collaboration 2018). In this model, there were con-tributions from fluctuations in the UVB before reionizationends, as well as from temperature fluctuations in the IGM.However to reproduce the sightlines with the highest effec-tive optical depths, it was crucial that residual islands ofneutral hydrogen persisted down below redshift z (cid:54) . α forest below redshift 6.In section 3, we show that, as well as matching the distri-bution of mean fluxes measured along individual sightlines,these models can reproduce troughs as long as observed. Insection 4, we model the distribution of LAEs around oneof our simulated troughs and in section 5 we present a testto distinguish between a late reionization model and a sce-nario with a fluctuating UVB. Finally in section 6 we stateour conclusions. To model the properties of the IGM at the end of reion-ization, we present here a carefully calibrated simulation ofhydrogen reionization. We ran mono-frequency cosmologi-cal radiative transfer simulations in post-processing using aton (Aubert & Teyssier 2008, 2010), using a setup similarto Chardin et al. (2015). aton is a GPU-based code thatsolves the radiative transfer equation using a moment-basedmethod with M1 closure for the Eddington tensor (e.g., Lev- F Becker+15Bosman+18This workKulkarni+19
Figure 1.
The evolution of the mean flux (blue points) in oursimulation down to z ∼ ermore 1984; Gnedin & Abel 2001; Aubert & Teyssier 2008).We perform this post-processing on a set of density fields ob-tained from a hydrodynamic simulation run with Tree-PMSPH code p-gadget3 (last described by Springel 2005). Oursimulation has a box size 160 cMpc/ h and 2 × gas anddark matter particles. The initial conditions were taken fromthe Sherwood simulation suite (Bolton et al. 2017) and as-sume cosmological parameters from Planck CollaborationXVI. (2014) with Ω m = 0 . Λ = 0 . h = 0 . b = 0 . σ = 0 .
829 and n = 0 . m gas = 6 . × M (cid:12) / h and the dark matterparticles have m DM = 3 . × M (cid:12) / h . The gravitationalsoftening length we assume is 3.1 ckpc/ h . We run these sim-ulations using a simplified star formation scheme, wherebyany gas particles with a density greater than 1,000 times themean cosmic baryon density and temperature less than 10 K are immediately converted into stars (Viel et al. 2004).This allows us to lower the time required to run the sim-ulation without significantly affecting the properties of theintergalactic medium.As we are performing the radiative transfer in post-processing, we neglect the hydrodynamical back-reaction ofthe gas to the photoheating. To compensate for this byadding 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 downto redshift 4. To allow us to update our input density field in
MNRAS , 1–11 (2019) roughs in the Ly α forest and H i islands the radiative transfer simulations frequently, we take snap-shots of the simulation every 40 Myr starting just belowredshift 20. We search for haloes in the simulation on thefly using a friends-of-friends algorithm. As aton requiresthe input density fields to be described on a cartesian grid,we map the SPH particles onto a 2048 grid (with cell size78.125 kpc/ h ) using the SPH kernel. To assign sources inthe simulation, we assume that each halo has an emissivityproportional to its mass (Iliev et al. 2006). The emissivityof each source is determined by summing the total mass ofsources within our volume and dividing the assumed totalemissivity within our volume among them, proportional totheir masses. We take a minimum halo mass 10 M (cid:12) / h . Wefind changing this minimum mass has hardly any effect onour results (once we have recalibrated the simulations tomatch the mean flux of the Ly α forest).We calibrate the ionizing emissivity assumed in the sim-ulation by running many simulations until the properties ofour simulated Ly α forest (described later) agree well withthe observations. As described above, the emissivity of eachsource is completely determined by the mass of the hosthalo and the volume emissivity we assume. When we cali-brate the simulations, we only change the volume emissiv-ity. This is equivalent to assuming some evolution in theproperties of the sources, such as their escape fraction, butwe do not explicitly treat these parameters. We make somechanges with respect to the simulation presented in Kulka-rni et al. (2019a), namely using a lower photon energy and asomewhat different emissivity evolution. The starting pointfor the emissivity evolution, which we modify as needed, isa galaxy-like component from Puchwein et al. (2019) andan increasing AGN-like contribution at lower redshifts. Al-though we do not model ionization by AGN directly (weonly model one class of source, with a spectrum assumedto be galaxy-like), we find that we need an increasing vol-ume emissivity at lower redshifts to match the mean fluxof the Ly α forest (described further below). We do not tryto separate our volume emissivity into contributions fromgalaxy-like and AGN-like components, but only account forthe total emissivity that we require to reproduce the Ly α forest statistics. We therefore leave a more detailed compar-ison of the relative contributions of galaxies and AGN to thehigh-redshift UVB to future work. During reionization, weassume that all of our sources emit as 30,000 K blackbodies.This spectrum was chosen to provide reasonable agreementwith measurements of the IGM temperature below redshift6. This is equivalent to an average photon energy of 17.1 eVand a photoionization cross-section 3 . × − cm − . Wearrive at these numbers assuming the optically thick limit,assuming all photons are absorbed locally (Pawlik & Schaye2011; Keating et al. 2018). This is likely a good approxi-mation during reionization but will cause us to overestimateour photoheating after reionization has ended. For this spec-trum, gas that is instantaneously ionized would reach a max-imum temperature of approximately 13,000 K.At lower redshifts, in a crude attempt to mimic the ef-fects of the heating of the IGM by He ii reionization (e.g.,Becker et al. 2011), we increase the photon energy in oursimulation. We do not explicitly treat two classes of sources(galaxies and AGN) separately, but instead just try to cap-ture the average effect of the evolution of the shape and am-plitude 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, andby redshift 4 we use an average photon energy of 23.8 eV.For the intermediate redshifts we interpolate between thevalues at z = 5 . t < ∆ x/c where c is the speed of light and ∆ x is the cell size of thegrid. The chemical network only follows the photoionization,photoheating and cooling of hydrogen. The cooling processeswe include are Hubble cooling, collisional ionization and ex-citation cooling (Cen 1992), recombination cooling (Hui &Gnedin 1997), Compton cooling (Haiman et al. 1996) andfree-free cooling (Osterbrock & Ferland 2006).As we run the simulation, we extract the ionization stateand temperature of the gas on the fly along a lightcone takenat a 30 degree angle through our grid. This allows us to ac-count for the rapid evolution of the IGM along a line of sight,rather than having to interpolate between individual snap-shots of the simulation as in Kulkarni et al. (2019a). We thenconstruct mock Ly α forest spectra along sightlines along ouroutput lightcone. The optical depth was computed using theapproximate fit to Voigt profiles presented in Tepper-Garc´ıa(2006). We show the evolution of the mean flux in our simu-lations in Figure 1 compared with observations from Beckeret al. (2015) and Bosman et al. (2018). In this work, we onlyfocus on analysis of the Ly α forest down to z ∼ . β transmission along our sim-ulated troughs, it is important to get the mean flux of theLy α forest close to the observed values at lower redshiftsfor adding in the contamination from the foreground Ly α forest.We show the cumulative distribution function (CDF)of effective optical depths in Figure 2 compared with thedistribution from Bosman et al. (2018). Here we use theconventional definition τ eff = − log( ¯ F ), where τ eff is the ef-fective optical depth and ¯ F is the mean flux measured in50 cMpc/ h segments of the Ly α forest. We emphasise thatthe simulations were calibrated to match the observed meanflux, and once this was reproduced the agreement betweenthe observed and simulated τ eff CDFs came out naturally.Another compilation of effective optical depths is presentedin Eilers et al. (2018), however due to systematic differencesbetween the two datasets it is not possible to find a singlereionization history that fits both measurements at once.We therefore limit our study here to a comparison with the
MNRAS000
MNRAS000 , 1–11 (2019)
L.C. Keating et al. P ( < e ff ) z=5.0 z=5.200.5 P ( < e ff ) z=5.4 z=5.60 2 4 6 8 eff P ( < e ff ) z=5.8 0 2 4 6 8 eff z=6.0Bosman et al. 2018 Figure 2.
The cumulative distribution of effective optical depths in our model (blue region), compared to the observations from Bosmanet al. (2018) (grey region). For the model, the shaded region denotes the 15 th and 85 th percentiles of 500 realisations of the CDFs. Forthe observations, the shaded area represents the “optimistic” and “pessimistic” bounds for the data. n [ s c M p c ] Kulkarni+19 Q H II McGreer+15Mason+18,19Greig+17,19Davies+18 H I [ s ] Haardt & Madau 12Becker & Bolton 13Calverley+11D'Aloisio+18 m f p [ p M p c / h ] Worseck+14 T [ K ] Walther+19Boera+19Becker+11, = 1.5Bolton+12
Planck 18
Figure 3.
Properties of our radiative transfer simulation (blue line) compared with the simulation presented in Kulkarni et al. (2019a)(grey dashed line). Top left: The evolution of the ionizing emissivity with redshift. Top middle: The filling fraction of ionized gas. Alsoshown are estimates from dark gap statistics of the Ly α forest (McGreer et al. 2015), quasar damping wing studies (Greig et al. 2017,2019; Davies et al. 2018) and estimates from the fraction of Lyman break galaxies that show Ly α emission (Mason et al. 2018, 2019).Top right: The evolution of the H i photoionization rate measured in ionized regions compared with measurements from Becker & Bolton(2013), Calverley et al. (2011) and D’Aloisio et al. (2018). The black dotted line is the Haardt & Madau (2012) model. Bottom left: Theevolution of the mean free path for photons at 912 ˚A as well as measurements from Worseck et al. (2014). Bottom middle: The meantemperature at mean density compared with estimates from Walther et al. (2019), Boera et al. (2019), Becker et al. (2011) and Boltonet al. (2012). Bottom right: The Thomson optical depth compared with the latest Planck results (Planck Collaboration 2018).
Bosman et al. (2018) compilation. To measure the τ eff dis-tributions, we create 500 realisations of the Bosman et al.(2018) sample and draw different sightlines from our sam-ple in the same redshift distribution. We show the 15 th and85 th percentiles of the different CDFs here. As in Kulkarniet al. (2019a), in this late reionization model we find goodagreement with the observed sample, although we still fail to reproduce the high opacity end of the distribution in the z = 5 . MNRAS , 1–11 (2019) roughs in the Ly α forest and H i islands h ]10 F r a c t i o n o f s i g h t li n e s Becker+15z max = 5.879z max = 5.979 0 50 100 150Longest trough along sightline [cMpc/ h ]456789101112 e ff a l o n g t r o u g h
50 0 50cMpc/ h c M p c / h T r o u g h l e n g t h [ c M p c / h ] N u m b e r o f s i g h t li n e s Figure 4.
Left panel: Distribution of the longest trough we measure along a sample of sightlines through our simulation. The filled bluehistogram only shows troughs with a maximum redshift z (cid:54) . z = 0.1. The vertical blackdotted line marks the length of the long trough presented in Becker et al. (2015). Middle panel: The effective optical depth integratedover the trough length. The colour of the pixels reflects the number of points in each bin in log scale. Also shown is the 2 σ limit onthe effective optical depth ( τ eff > . require a slightly different evolution to the assumed emis-sivity in that work. First, we assume a lower photon en-ergy than in that work to try and recover temperatures thatare closer to the Ly α forest measurements. Using this lowerphoton energy requires us to have a slightly faster reioniza-tion to recover the observed mean flux – since the gas doesnot reach as high a temperature as it is ionized, we requireit to be ionized somewhat later so it has less time to cooldown. Second, since we construct the spectra along a light-cone through the simulation, there are slight changes to themean flux we recover from the simulations and therefore werequire a different emissivity evolution. Third, as mentionedabove, we also would like this simulation to match the meanflux of the Ly α forest below redshift 5. We find that thisrequires a rapid upturn of the ionizing emissivity just aboveredshift 5 (which could perhaps be attributed to the onsetof an increasing contribution of AGN to the photoionizingbackground 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 madealong a lightcone taken in one direction through our simula-tion volume. If this exercise was repeated for a larger sampleof sightlines, the required emissivity evolution would likelybe somewhat different. Furthermore, as already mentioned,the different measurements of the mean flux in this redshiftrange by Bosman et al. (2018) and Eilers et al. (2018) arealso not fully consistent within the quoted respective errors.We therefore caution against over-interpreting the emissiv-ity evolution presented here.Although a decline in the galaxy contribution to theUVB is predicted in models such as Haardt & Madau (2012)and Puchwein et al. (2019), it is difficult to find a physical in-terpretation for why the galaxy-like component should dropas rapidly as required here. It is perhaps also concerningthat this rapid drop is required just as reionization ends inour simulations. However, the only other way to decreasethe mean flux would be if the IGM was much colder, andit would be equally difficult to explain why it should coolso quickly since it needs to be reasonably warm to achievethe observed mean flux at redshift 6. The AGN-like com- ponent that we invoke to match the mean flux at lowerredshifts becomes dominant over the galaxies much ear-lier than in other, more physically motivated models (e.g.,Puchwein et al. 2019) and is in tension with the shape ofthe UV background inferred from metal lines (Boksenberget al. 2003) and measurements of the ionizing emissivity ofgalaxies (Steidel et al. 2018) at somewhat lower redshifts.We leave a more careful investigation of the relation of themean flux of the Ly α forest and the ionizing emissivity tofuture work, preferably using higher resolution simulationswith more careful modelling of the gas temperature, and animproved measurement of the mean flux and the flux PDF.The other panels of Figure 3 compare our simulationto other properties of the IGM. In this model, the volumeis 50 per cent ionized at z = 6 . z = 5 .
2. In general, the agreement between the modeland the observed IGM properties is good. However, as dis-cussed above, since we not self-consistently modelling theQSO-like component of our emissivity, our temperatures be-low redshift 5 in the voids maybe be underestimated andour photoionization rate will be overestimated. We empha-sise though that information from this redshift range is onlyincluded as the foregrounds to our Ly β forest modellingand does not affect our main conclusions. The main differ-ences between the model presented here and Kulkarni et al.(2019a) are driven by the temperature of the gas. Our choiceof a lower photon energy here brings the temperature atmean density in the simulations closer to the observations,however as discussed above this is at the expense of intro-ducing a rather sharp drop in the evolution of the emissivityat z ∼
6. To achieve the same mean flux, the lower tempera-tures require a higher photoionization rate. This also seemsto be closer to the observed values. We also show for com-parison the photoionization rate from the Haardt & Madau(2012) model, as this was used when running the underly-ing hydrodynamic simulation. We caution though that thepublished temperature measurements rely heavily on cali-bration with numerical simulations that assume a homoge-neous UV background, and do not account for the effects ofinhomogeneous reionization and spatial fluctuations in the
MNRAS000
MNRAS000 , 1–11 (2019)
L.C. Keating et al. 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.
We next measure the longest trough along 250,000 sightlinesalong our extracted lightcone. We post-process the spectrato match the observations, by rebinning onto coarser pix-els, convolving with an appropriate instrument profile andadding noise. We measure the length of the troughs by mea-suring out to a maximum redshift z max = 5 .
879 (the highredshift point of the Becker et al. 2015 trough). We thensearch for spikes in transmission below this redshift which weassume are regions of the spectrum that show transmissiongreater than the 1 σ noise level for at least four consecutivepixels and that have a combined significance of more than5 σ . We then define the regions in between these transmissionspikes as the troughs, and record the longest trough mea-sured for each sightline. As shown in the left panel of Figure4, our simulation does produce troughs in the Ly α forest aslong as the 110 cMpc/ h trough of ULAS J0148+0600, how-ever these are very rare in our simulation (approximately 1in 10,000). In a recent paper, Giri et al. (2019) also presentedan investigation of neutral islands at the end of reionization,and noted that they find islands long enough to reproducethe long trough of ULAS J0148+0600. Their simulation vol-ume is much larger (714 cMpc vs.
160 cMpc/ h in this work)and is therefore more suited for studying the statistics ofthese long troughs at the end of reionization. Detailed com-parisons with observations will however also require that thesimulations are calibrated to match the properties of the Ly α forest, such as in this work.The long troughs we do find are as dark as measured byBecker et al. (2015), who estimated that the effective opticaldepth was greater than 7.2 along the absorption trough ofULAS J0148+0600 (middle panel of Figure 4). We find thattroughs with length up to ∼
80 cMpc/ h have a relativelyhigh incidence rate in the simulation, but the incidence ratefalls rapidly above this point. The most natural explana-tion for this is that reionization is still ending too early inour simulations. As discussed above, there are 50 cMpc/ h segments of the Ly α forest that are observed to be totallyopaque at z = 5 . z = 0 . h . However, we note that we were notsuccessful in finding a model where reionization ended laterwhile still matching the observed mean flux.Another point worth making is that the size of our sim-ulation volume is still small in the context of reionizationstudies (e.g., Iliev et al. 2014). We therefore do not have alarge sample of neutral islands at low redshift in our simula-tion. This is shown in the right panel of Figure 4 – all of thelong troughs are clustered in one section of the box. Thismeans that the deficit of long troughs in our simulationscould be caused by an unfavourable geometry in these re- maining neutral islands, i.e. they may just happen to not bevery extended in the direction in which we extract our light-cone. A stronger statement about the incidence rate of longtroughs for a given ionization history will therefore requirelarger volumes while still maintaining resolution comparableto the simulation presented here. We therefore leave a moredetailed study of the incidence rate of these long troughs forfuture work.An example of a long (98.2 cMpc/ h ) trough is shown inFigure 5, as well as the associated physical quantities of thegas along the trough. We find that the region spanned bythe trough is coincident with an island of neutral hydrogen.It is clear from this figure that the shape of the neutralislands is directly linked to the trough length, and one couldimagine that a sightline taken at a small angle to the lineof sight shown here would produce a longer trough. As seenin Figure 5, the UVB is also lower in the vicinity of theneutral island, as it is only being illuminated by sources inone direction. This means that the completely neutral gasis surrounded by gas with an H i fraction of about 10 − ,still neutral enough to saturate the Ly α transmission. Wealso find Ly β transmission along the trough, as is observed.This can be explained by small pockets of lightly ionized gaswithin this neutral island. The photoionization rate insidethese pockets is still low, however the gas is ionized enoughto allow transmission of Ly β , since Ly β has an oscillatorstrength ∼ α . α EMITTERS
A convincing model of the long trough in ULAS J0148+0600must also explain the observed deficit of LAEs around thatsightline. Since we construct our spectra along a planar light-cone taken at an angle wrapped through our simulation, wecannot make use of periodic boundary conditions to centrea simulated trough in the middle of the lightcone. We there-fore pick the sightline we are interested in (in this case, thesightline shown in Figure 5) and rerun the radiative trans-fer simulation to extract a new lightcone centred on thissightline. This means that our analysis of the LAEs aroundthe trough is restricted to a single sightline. However, sincethe size of our volume is close to the field of view of HSCand also our longest troughs are co-spatial (Figure 4) ananalysis of more than one trough in the simulation wouldnot really remove the effects of cosmic variance. We chosethis trough because, although it is slightly shorter than thetrough of ULAS J0148+0600, it showed Ly β emission withsimilar properties to the observed long trough, which showsonly one Ly β transmission spike in the redshift range cov-ered by the narrowband filter used in Becker et al. (2018).To model the LAEs around our simulated trough, wefollow the method described in detail in Weinberger et al.(2019). We note that this method is similar to the LAEmodelling in Davies et al. (2018) and Becker et al. (2018),although as we will describe later we also take the effectof the local IGM opacity into account. We first assign UVluminosities to the dark matter haloes by abundance match-ing to the UV luminosity function at redshift 5.9 (Bouwenset al. 2015). We assume a 50 Myr duty cycle such that onlya fraction of our galaxies are UV bright at a given time,reflecting the bursty nature of star formation (Trenti et al. MNRAS , 1–11 (2019) roughs in the Ly α forest and H i islands l o g f H I -40040 c M p c / h l o g ( / ) -40040 c M p c / h l o g ( H I / s ) -40040 c M p c / h l o g ( T / K ) -40040 c M p c / h L y F l u x -40040 c M p c / h L y F l u x c M p c / h l o g f H I l o g ( / ) l o g ( H I / s ) l o g ( T / K ) L y e ff L y e ff Figure 5.
Properties of a sightline containing a long (98.2 cMpc/ h ) trough. The left panel shows a 1D sightline and the right panelshows 2D maps. The white dashed line on the right column corresponds to the redshift range where we measure the Ly α absorptiontrough. From top to bottom: the neutral hydrogen fraction, gas density, H i photoionization rate, gas temperature, Ly α flux (left) andLy α effective optical depth (right) and Ly β flux (left) and Ly β effective optical depth (right) measured over 50 cMpc/ h segments. Thegrey shaded region on the left panels shows the redshift range corresponding to the FWHM of the narrowband filter used to search forLAEs. α emission to each halobased on its UV luminosity using the relation P (REW | M UV ) ∝ exp (cid:18) − REWREW c ( M UV ) (cid:19) , (1)where REW is the rest-frame equivalent width, M UV is theUV magnitude and REW c is defined asREW c ( M UV ) = 23 + 7( M UV + 21 .
9) + 6( z − . (2) The normalization is defined such that the population hasREW min (cid:54) REW (cid:54)
REW max , where we choose REW max =300 ˚A. The minimum REW is defined in terms of M UV asREW min = −
20 ˚A M UV < − . , . M UV > − . , −
20 + 6( M UV + 21 . ˚A otherwise . (3)We note that the assignment of UV magnitude and REW MNRAS000
20 + 6( M UV + 21 . ˚A otherwise . (3)We note that the assignment of UV magnitude and REW MNRAS000 , 1–11 (2019)
L.C. Keating et al. × P ( R E W ) Shibuya+18 (D)Shibuya+18 (UD) 42.5 43.0 43.5 44.0 L Ly [erg/s]10 ( L L y ) [( l o g L ) M p c ] IntrinsicAfter attenuationKonno+18 50 0 50cMpc/ h c M p c / h N B Figure 6.
Left: Rest-frame equivalent width distribution of our simulated LAEs (blue line) compared to observations from Shibuya et al.(2018) in Deep (D) and UltraDeep (UD) fields. Middle: The LAE luminosity function of our simulated LAEs (blue line) compared toobservations from Konno et al. (2018). In both panels, the grey line shows the intrinsic distribution and the blue line shows the distributionafter the IGM attenuation is taken into account. The solid lines are the median values of 50 realisations of the LAE population andthe shaded regions represent the 15 th /85 th per cent confidence intervals. Right: The spatial distribution of one realisation of our LAEscentred around the trough shown in Figure 5. The colour of the points represents their NB816 magnitude. The dashed black lines showcircles centred around the trough with radii increasing by 10 cMpc/ h for each circle. h ]0.000.010.020.030.040.050.060.07 L A E [ h / c M p c ] Becker+18 0 20 40 60 80Radius [cMpc/ h ]0.600.650.700.750.800.850.90 T I G M h ]0.000.010.020.030.040.050.060.07 U V [ h / c M p c ] M UV < -21M UV < -20.5M UV < -20 Figure 7.
Left: The surface density of LAEs measured in annuli centred around the trough. The black points are the Becker et al.(2018) blue line is the median surface density measured in our LAE realisation. The shaded regions represent the 15 th /85 th per centconfidence intervals. Middle: The averaged IGM transmission measured in annuli centred around the trough, which we define as the ratioof the observed and intrinsic LAE luminosities. Right: The surface density of LBGs with M UV < −
20 (red), M UV < − . M UV < −
21 (blue) around the same sightline. is a random process, since only a fraction of the haloes areUV bright at a given time and the REW of each halo isdrawn from a distribution. We therefore repeated this pro-cess 50 times, to generate 50 different realisations of LAEpopulations.We next accounted for the attenuation of Ly α emissiondue to neutral gas around the haloes. The default resolutionof our radiative transfer simulations is too coarse to resolvethe circumgalactic medium (CGM) of halos, so we re-extracta higher resolution skewer of gas density and peculiar veloc-ity in a 10 cMpc/ h radius around each halo. The resolutionof these skewers is 9.7 ckpc/ h . We insert this short highresolution sightline into the longer lower resolution sight-line to model the attenuation of Ly α emission due to boththe CGM and IGM. Since we do not know the ionizationstate of the gas in the high resolution sightline, we assumethat the gas is in ionization equilibrium calculated using thephotoionization rates from our radiative transfer simulationplus 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 mag- nitude to get the flux at rest-frame 1600 ˚A, and assumingthat F λ ∝ λ − . To model the Ly α emission, we add a gaus-sian emission line with a REW given by the prescriptionoutlined above. As in Weinberger et al. (2019), we also as-sume some velocity shift between the redshift of the hosthalo and the redshift of the Ly α emission, which accountsfor the complex radiative transfer of the Ly α line out ofthe halo. We take this to be a constant factor times thecircular velocity of the host halo, with the constant factorchosen to ensure good agreement with the rest-frame equiv-alent widths and luminosities of observed LAEs. We thenattenuate the spectrum by multiplying by the normalisedLy α transmission we calculate along each sightline. In thiswork we found that a velocity shift of 1 . v circ was a goodchoice. Using these spectra, we then measure the magnitudeof the LAEs in the i NB
816 bands and select LAEsthat have NB (cid:54) . i − NB (cid:62) . r MNRAS , 1–11 (2019) roughs in the Ly α forest and H i islands minosity and rest-frame equivalent width for each LAE byintegrating over our constructed LAE spectra. As in Wein-berger et al. (2019), the agreement between the observedand simulated LAEs is very good, i.e. our modelled LAEpopulation agrees with the observed LAE luminosity func-tion and rest-frame equivalent width distribution from theSILVERRUSH survey (Figure 6), however we note that ourbox is too small to reproduce the brightest LAEs (seen bythe rapid cutoff of our modelled luminosity function).We next compare the spatial distribution of LAEsaround our simulated trough (right panel of Figure 6). Weshow here the projected position of LAEs coloured by their NB
816 magnitude. As in Becker et al. (2018), we find a largedeficit of LAEs in the immediate vicinity of the trough. Thisis to be expected, since the low density regions will be thelast to ionize. We therefore expect our remaining neutral is-lands to trace underdense regions of the IGM. Furthermore,since the Ly α optical depth along the trough is high, wealso expect the attenuation of Ly α emission to be strongerin the region around the trough. In the left panel of Fig-ure 7, we show the LAE surface density computed in annulicentred around the trough. The results from our simulatedLAE are shown in blue (again, the shaded regions mark the15 th /85 th percentiles of our 50 different LAE realisationsand the solid line is the median value). The results fromBecker et al. (2018) are shown in black. Similar to the ob-servations, we find that the surface density of LAEs beginsto decline rapidly in the inner 30 cMpc/ h around the trough.This is predominantly due to the trough passing through anunderdense region, although we note that that this effect iscomplemented by a decrease in the IGM transmission frac-tion of about 10 per cent close to the trough (middle panelof Figure 7). Here we define the IGM transmission fractionas T IGM = L obs / L int where L is the observed/intrinsic LAEluminosity (Mesinger et al. 2015). This lower transmissionfraction is in agreement with Sadoun et al. (2017), who sug-gest that the observed number of LAEs will be lower inregions with a lower photoionizing background, due to theincreased filling factor of self-shielded systems in the CGMof the host galaxies. However, it does not appear to be thedominant effect here.In the right panel of Figure 7, we show the surfacedensity of UV-selected galaxies that have M UV < − − . −
21. For the cut with M UV < −
20, we findthat the surface density far from the trough is well matchedto the LAEs. However, the decline in surface density iseven stronger for this selection, already beginning to de-cline at a radius 50 cMpc/ h from the trough and highlight-ing that this sightline passes through an underdense region.The median mass of the LAEs that pass the colour cutsused here is 4 . × M (cid:12) / h . For the UV selected galax-ies, the median masses are (1 . , . , . × M (cid:12) / h for M UV < ( − , − . , − As shown in Becker et al. (2018), the fluctuating UVB modelpresented in Davies & Furlanetto (2016) can also reproducethe deficit of LAEs in high τ eff regions. In this case, a shortmean free path leads to large fluctuations in the UVB. Thisproduces regions that, while still ionized, have a H i fraction high enough to saturate transmission of Ly α . We do not findthat this is the case in our simulations, as once the ionizedregions have percolated the mean free path becomes longand 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 thatproduce the high τ eff sightlines. In our case, however, this isbecause reionization has not yet ended and the voids hostthe residual neutral islands. In both models, the underdenseregions host the high τ eff sightlines because they are far fromthe regions where the space density of ionizing sources ishighest. Linking the Davies & Furlanetto (2016) model toa late reionization scenario removes the main issues withthat explanation. First, that the mean free path needs tobecome much shorter at z > α forest. We alsotune the emissivity in our model to match the observed meanflux, however once this is fixed then the reionization historyof our preferred model can be used to make predictions forhigher redshifts. In this sense, our late reionization modelcan be thought of as a dynamic version of the Davies &Furlanetto (2016) model.Ideally one would still like to discriminate between thetwo models though, to answer if indeed the high τ eff regionsare mildly ionized or completely neutral. One difference be-tween the two scenarios is that the fluctuating UVB modelpredicts that the underdense regions will always coincidewith high τ eff sightlines, and the low τ eff sightlines will liein overdense regions (Davies et al. 2018). In a late reioniza-tion model though, the underdense regions will align withhigh τ eff sightlines only in the time before they have beenionized. Once the neutral islands have vanished, the voidsshould instead correspond to the lowest τ eff regions, as mostof the Ly α forest transmission is dominated by the voids.These recently ionized regions should also be very hot, sup-pressing recombinations and further enhancing transmission(D’Aloisio et al. 2015).We present an example of two such sightlines in Fig-ure 8. The left panel shows a map of τ eff measured in 50cMpc/ h segments centred on the redshift corresponding topeak transmission of the NB
816 filter. The blue and redcircles mark two regions of high and low τ eff . In the middlepanel, we show the mean surface density of LAEs in thatregion. The surface density is taken to be inversely propor-tional to the square of the distance to the tenth nearestneighbour and the map has been smoothed with a gaus-sian kernel. The red and blue circles correspond to the sameregions as in the left panel. It is clear that both of these re-gions show a deficit of LAEs, suggesting that both of theseare underdense regions. This is shown explicitly in the rightpanel, where we compare the surface density of LAEs as afunction of radius for the two lines of sight. In the Davies& Furlanetto (2016) model, you would instead expect thatthe low τ eff region should show an overdensity of LAEs. Thissuggests that future observations searching for LAEs aroundsightlines that show significant Ly α forest transmission maybe interesting for differentiating between these two models. MNRAS000
816 filter. The blue and redcircles mark two regions of high and low τ eff . In the middlepanel, we show the mean surface density of LAEs in thatregion. The surface density is taken to be inversely propor-tional to the square of the distance to the tenth nearestneighbour and the map has been smoothed with a gaus-sian kernel. The red and blue circles correspond to the sameregions as in the left panel. It is clear that both of these re-gions show a deficit of LAEs, suggesting that both of theseare underdense regions. This is shown explicitly in the rightpanel, where we compare the surface density of LAEs as afunction of radius for the two lines of sight. In the Davies& Furlanetto (2016) model, you would instead expect thatthe low τ eff region should show an overdensity of LAEs. Thissuggests that future observations searching for LAEs aroundsightlines that show significant Ly α forest transmission maybe interesting for differentiating between these two models. MNRAS000 , 1–11 (2019) L.C. Keating et al. h c M p c / h h c M p c / h h ]0.000.010.020.030.040.050.060.07 L A E [ h / c M p c ] High eff
Low eff e ff L A E / L A E Figure 8.
Left: Map of the effective optical depth measured in a 50 cMpc/ h segment centred on the peak transmission wavelength of thenarrowband filter. The red and blue circles mark out regions that have low and high effective optical depths. Middle: Smoothed averageLAE surface density map for the same field as in the right panel. Right: Surface density of LAEs for the low and high effective opticaldepth regions. We have shown here that long Gunn-Peterson troughs ob-served in the Ly α forest above redshift 5 can be explained byislands of neutral hydrogen in the late reionization scenariopresented in Kulkarni et al. (2019a). We have demonstratedthat this model also explains the significant deficit of LAEsdetected around these high opacity sightlines by Becker et al.(2018). While these observations can also be explained bya model in which the amplitude of the UVB is low in un-derdense regions, we argue that the model presented here ismore consistent with the observed evolution of the mean freepath for ionizing photons. We further present predictions forthe distribution of LAEs around the low opacity sightlineswhich should allow future observations to discriminate be-tween these two scenarios. This however will still not bea definitive test of neutral gas in the IGM, so other testsshould also be pursued (Malloy & Lidz 2015).We note that the reionization model presented here stilldoes not reproduce the most opaque sightlines in the range5 . < z < . th /85 th interval. Searching formore saturated sightlines at low redshifts would be useful inconstraining whether these sightlines are rare. If they turnout not to be outliers, this may call for a reionization historythat ends even later than shown here. This would also likelyincrease the incidence rate of long ( >
100 cMpc/ h ) troughs.We note also that the temperature of the IGM is still a sig-nificant uncertainty in our models, as it can have a largeeffect on the mean flux (as the recombination rate and thusthe effective optical depth for a given UVB amplitude de-pend rather strongly on temperature). More accurate mod-els of the ionization history of the IGM would benefit frommore estimates of the temperature at high redshift, perhapsalso using methods that account for the pressure effects forspatial fluctuations of the temperature-density relation dueto reionization (e.g., O˜norbe et al. 2019, Puchwein et al. inprep).As more high signal-to-noise Ly α forest sightlines areaccumulated and measurements of the mean flux and fluxPDF further improve, there will be benefit in moving fromstudies of integrated quantities such as the effective opticaldepth to more detailed studies of individual features. Fu-ture work modelling the occurrence rate and distributionof these troughs, as well as the transmission peaks between them, should provide more detailed constraints on the lat-ter half of reionization (Chardin et al. 2018a; Garaldi et al.2019b). Tying these Ly α forest measurements to the galaxypopulation may also yield additional information such asthe sources that contribute to reionization (Kakiichi et al.2018), especially once JWST is operational. Studies utilis-ing also the Ly β forest should allow to push constraints onthe reionization history outside of the immediate vicinity ofhigh redshift quasars to z ∼ . ACKNOWLEDGEMENTS
REFERENCES
Aubert D., Teyssier R., 2008, MNRAS, 387, 295Aubert D., Teyssier R., 2010, ApJ, 724, 244Becker G. D., Bolton J. S., 2013, MNRAS, 436, 1023Becker G. D., Bolton J. S., Haehnelt M. G., Sargent W. L. W.,2011, MNRAS, 410, 1096Becker G. D., Bolton J. S., Madau P., Pettini M., Ryan-WeberE. V., Venemans B. P., 2015, MNRAS, 447, 3402Becker G. D., Davies F. B., Furlanetto S. R., Malkan M. A., BoeraE., Douglass C., 2018, ApJ, 863, 92MNRAS , 1–11 (2019) roughs in the Ly α forest and H i islands Boera E., Becker G. D., Bolton J. S., Nasir F., 2019, ApJ, 872,101Boksenberg A., Sargent W. L. W., Rauch M., 2003,arXiv:0307557,Bolton J. S., Becker G. D., Raskutti S., Wyithe J. S. B., HaehneltM. G., Sargent W. L. W., 2012, MNRAS, 419, 2880Bolton J. S., Puchwein E., Sijacki D., Haehnelt M. G., Kim T.-S.,Meiksin A., Regan J. A., Viel M., 2017, MNRAS, 464, 897Bosman S. E. I., Fan X., Jiang L., Reed S., Matsuoka Y., BeckerG., Haehnelt M., 2018, MNRAS, 479, 1055Bouwens R. J., Illingworth G. D., Oesch P. A., Caruana J., Hol-werda B., Smit R., Wilkins S., 2015, ApJ, 811, 140Calverley A. P., Becker G. D., Haehnelt M. G., Bolton J. S., 2011,MNRAS, 412, 2543Cen R., 1992, ApJS, 78, 341Chardin J., Haehnelt M. G., Aubert D., Puchwein E., 2015, MN-RAS, 453, 2943Chardin J., Puchwein E., Haehnelt M. G., 2017, MNRAS, 465,3429Chardin J., Haehnelt M. G., Bosman S. E. I., Puchwein E., 2018a,MNRAS, 473, 765Chardin J., Kulkarni G., Haehnelt M. G., 2018b, MNRAS, 478,1065D’Aloisio A., McQuinn M., Trac H., 2015, ApJ, 813, L38D’Aloisio A., Upton Sanderbeck P. R., McQuinn M., Trac H.,Shapiro P. R., 2017, MNRAS, 468, 4691D’Aloisio A., McQuinn M., Davies F. B., Furlanetto S. R., 2018,MNRAS, 473, 560Davies F. B., Furlanetto S. R., 2016, MNRAS, 460, 1328Davies F. B., Becker G. D., Furlanetto S. R., 2018, ApJ, 860, 155Dijkstra M., Wyithe J. S. B., 2012, MNRAS, 419, 3181Eilers A.-C., Davies F. B., Hennawi J. F., 2018, ApJ, 864, 53Fan X., et al., 2006, AJ, 132, 117Furlanetto S. R., Oh S. P., 2009, ApJ, 701, 94Garaldi E., Compostella M., Porciani C., 2019a, MNRAS, 483,5301Garaldi E., Gnedin N., Madau P., 2019b, ApJ, 876, 31Giri S. K., Mellema G., Aldheimer T., Dixon K. L., Iliev I. T.,2019, Monthly Notices of the Royal Astronomical Society,p. 2192Gnedin N. Y., Abel T., 2001, New Astron., 6, 437Greig B., Mesinger A., Haiman Z., Simcoe R. A., 2017, MNRAS,466, 4239Greig B., Mesinger A., Ba˜nados E., 2019, MNRAS, 484, 5094Haardt F., Madau P., 2012, ApJ, 746, 125Haiman Z., Thoul A. A., Loeb A., 1996, ApJ, 464, 523Hui L., Gnedin N. Y., 1997, MNRAS, 292, 27Iliev I. T., Mellema G., Pen U.-L., Merz H., Shapiro P. R., AlvarezM. A., 2006, MNRAS, 369, 1625Iliev I. T., Mellema G., Ahn K., Shapiro P. R., Mao Y., Pen U.-L.,2014, MNRAS, 439, 725Kakiichi K., et al., 2018, MNRAS, 479, 43Keating L. C., Puchwein E., Haehnelt M. G., 2018, MNRAS, 477,5501Konno A., et al., 2018, PASJ, 70, S16Kulkarni G., Keating L. C., Haehnelt M. G., Bosman S. E. I.,Puchwein E., Chardin J., Aubert D., 2019a, MNRAS, 485,L24Kulkarni G., Worseck G., Hennawi J. F., 2019b, Monthly Noticesof the Royal Astronomical Society, 488, 1035Levermore C. D., 1984, J. Quant. Spectrosc. Radiative Transfer,31, 149Lidz A., Malloy M., 2014, ApJ, 788, 175Lidz A., Oh S. P., Furlanetto S. R., 2006, ApJ, 639, L47Lidz A., McQuinn M., Zaldarriaga M., Hernquist L., Dutta S.,2007, ApJ, 670, 39Malloy M., Lidz A., 2015, ApJ, 799, 179 Mason C. A., Treu T., Dijkstra M., Mesinger A., Trenti M., Pen-tericci L., de Barros S., Vanzella E., 2018, ApJ, 856, 2Mason C. A., et al., 2019, MNRAS,McGreer I. D., Mesinger A., D’Odorico V., 2015, MNRAS, 447,499Mesinger A., 2010, MNRAS, 407, 1328Mesinger A., Aykutalp A., Vanzella E., Pentericci L., Ferrara A.,Dijkstra M., 2015, MNRAS, 446, 566O˜norbe J., Davies F. B., Luki´c Z., Hennawi J. F., Sorini D., 2019,MNRAS, 486, 4075Osterbrock D. E., Ferland G. J., 2006, Astrophysics of gaseousnebulae and active galactic nuclei. University Science BooksPawlik A. H., Schaye J., 2011, MNRAS, 412, 1943Planck Collaboration 2018, arXiv:1807.06209,Planck Collaboration XVI. 2014, A&A, 571, A16Puchwein E., Haardt F., Haehnelt M. G., Madau P., 2019, MN-RAS, 485, 47Sadoun R., Zheng Z., Miralda-Escud´e J., 2017, ApJ, 839, 44Shibuya T., et al., 2018, PASJ, 70, S14Springel V., 2005, MNRAS, 364, 1105Steidel C. C., Bogosavljevi´c M., Shapley A. E., Reddy N. A.,Rudie G. C., Pettini M., Trainor R. F., Strom A. L., 2018,ApJ, 869, 123Tepper-Garc´ıa T., 2006, MNRAS, 369, 2025Trenti M., Stiavelli M., Bouwens R. J., Oesch P., Shull J. M.,Illingworth G. D., Bradley L. D., Carollo C. M., 2010, ApJ,714, L202Viel M., Haehnelt M. G., Springel V., 2004, MNRAS, 354, 684Walther M., O˜norbe J., Hennawi J. F., Luki´c Z., 2019, ApJ, 872,13Weinberger L. H., Haehnelt M. G., Kulkarni G., 2019, MNRAS,485, 1350Worseck G., et al., 2014, MNRAS, 445, 1745This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS000