Anomaly in the Opacity of the Post-reionization Intergalactic Medium in the Lyα and Lyβ Forest

Our original data sample comprises 34 quasar spectra at 5.77 ≤ z_em ≤ 6.54 and has been publicly released (Eilers et al. 2018) via the igmspec database (Prochaska 2017) and the zenodo platform.⁵ Previously, we have used this data set to analyze the sizes of quasar proximity zones (Eilers et al. 2017a), as well as the redshift evolution of the IGM opacity within the $\mathrm{Ly}\alpha$ forest (Eilers et al. 2018). For the study of the joint $\mathrm{Ly}\alpha$ and $\mathrm{Ly}\beta$ forest opacity that we present in this paper, we take a subset of 19 quasar spectra with S/N ≥ 10 per 10 km s⁻¹ pixel. All spectra have been observed at optical wavelengths (4000–10000 Å) with the Echellette Spectrograph and Imager (ESI; Sheinis et al. 2002) at the Keck II Telescope between the years of 2001 and 2016. The data were collected from the Keck Observatory Archive⁶ and complemented with our own observations. All observations used slit widths ranging from 075 to 10, resulting in a resolution of R ≈ 4000–5400.

The details of all individual observations as well as detailed information about the data reduction process can be found in Eilers et al. (2017a, 2018). We make further improvements on our data reduction in order to avoid biases in our analysis of the opacities within the $\mathrm{Ly}\beta$ forest relative to the opacities within the $\mathrm{Ly}\alpha$ forest, which could arise due to potential flux calibration issues or in the process of coadding individual exposures. Hence, we correct all final extracted and coadded spectra with a power law, if necessary, to match the observed photometry in the i- and z-band (see Appendix A for details on this procedure). The properties of all quasars in our data sample are listed in Table 1 in Eilers et al. (2018).

We normalize the quasar spectra by their continuum emission, which we estimate using a principal component analysis (PCA) reconstruction. The idea of the PCA is for each continuum spectrum to be represented by a mean spectrum plus a finite number of weighted principal component spectra (Suzuki et al. 2005; Suzuki 2006; Pâris et al. 2011). Because the quasar spectra in our data set experience substantial absorption blueward of the $\mathrm{Ly}\alpha$ emission line from intervening neutral hydrogen in IGM, we follow the procedure suggested by Pâris et al. (2011) and estimate the coefficients for the reconstructed continuum emission using only pixels redward of the $\mathrm{Ly}\alpha$ line, and then apply a projection matrix to obtain the coefficients for the continuum spectrum covering the whole spectral range between $1020\,\mathring{\rm A} \leqslant {\lambda }_{\mathrm{rest}}\leqslant 1600\,\mathring{\rm A}$. Because the PCA components do not cover bluer wavelengths, we extend the estimated continuum to bluer wavelengths by appending the composite quasar spectrum from Shull et al. (2012) continuously at λ_rest < 1020 Å; for details, see Eilers et al. (2018).

We estimate the opacity of the IGM by means of the effective optical depth, i.e.,

$$\begin{eqnarray}&&{\tau }_{\mathrm{eff}}=-\mathrm{ln}\langle F\rangle ,\end{eqnarray} \tag{ 4 }$$

where $\langle F\rangle$ denotes the continuum normalized transmitted flux averaged over discrete spectral bins along the line of sight to each quasar. In this study, we choose spectral bins of 40 cMpc, instead of a bin size of 50 cMpc$\,{h}^{-1}\,$ used in previous work (Becker et al. 2015; Bosman et al. 2018; Eilers et al. 2018). This smaller bin is chosen such that the $\mathrm{Ly}\beta$ forest region is better sampled. Note that whenever the average flux $\langle F\rangle$ is negative or detected with less than 2σ significance, we adopt a lower limit on the optical depth at the 2σ level, in accordance with previous work (Fan et al. 2006; Becker et al. 2015).

The Lyα forest lies between the $\mathrm{Ly}\alpha$ and $\mathrm{Ly}\beta$ emission lines at the rest-frame wavelengths λ_Lyα = 1215.67 Å and λ_Lyβ = 1025.72 Å, respectively, whereas the Lyβ forest covers the wavelength region between the $\mathrm{Ly}\beta$ and the $\mathrm{Ly}\gamma$ (at rest-frame wavelength λ_Lyγ = 972.54 Å) emission lines. However, we do not conduct the IGM opacity measurements within the whole wavelength region, but the proximity zones around each quasar are excluded as in Eilers et al. (2018), i.e., we mask the region around each quasar that is heavily influenced by the quasar's own radiation, resulting in enhanced transmitted flux. Additionally, we choose the rest-frame wavelengths 1030 and 975 Å as the respective minimum wavelengths for our measurements in the Lyα and Lyβ forest. Thus, for a typical quasar at z_em ∼ 6 with a proximity zone of R_p ≈ 5 pMpc, we measure the opacity of the Lyα and Lyβ forests between the observed wavelengths 7210–8389 Å and 6825–7078 Å, respectively.

Additionally, as described in detail in Eilers et al. (2018), we correct for small offsets in the zero level of each spectrum that might have been introduced due to improper sky subtraction, and mask the spectral regions around DLAs, in order to avoid any biases in the estimation of the IGM opacity.⁷

Note that we do not attempt to correct for any metal line contamination within the $\mathrm{Ly}\alpha$ and $\mathrm{Ly}\beta$ forests. Faucher-Giguère et al. (2008b) show in their Figure 8 that the relative metal correction to τ_eff in the $\mathrm{Ly}\alpha$ forest decreases with increasing redshift from 13% at z = 2%–5% at z = 4. Assuming a monotonic increase in the enrichment of the IGM with time, we expect to have a negligible metal contamination at z ≈ 6.

We first compare our measurements to predictions of the $\mathrm{Ly}\beta$ forest opacity from a published reionization model with a spatially inhomogeneous UVB by Davies et al. (2018b). The model consists of UVB fluctuations from an independent seminumerical simulation following Davies & Furlanetto (2016), which we summarize briefly here. First, using the 21cmFAST code, we generated a set of cosmological initial conditions in a (400 cMpc)³ volume. From these initial conditions, we derived the locations of dark matter halos following the excursion set method of Mesinger & Furlanetto (2007), and shifted their positions to z = 6 with the Zel'dovich approximation (Zel'dovich 1970). The UV luminosities of the halos were then computed via abundance matching (e.g., Vale & Ostriker 2004) to the Bouwens et al. (2015) UV luminosity function. Finally, the ionizing radiation field was computed assuming a halo mass-independent conversion from UV luminosity to ionizing luminosity, and a spatially varying mean free path of ionizing photons (Davies & Furlanetto 2016). We then impose the fluctuating UVB onto the Nyx skewers and recompute the ionization state of the gas assuming ionization equilibrium.

Because the UVB fluctuations were computed in an independent cosmological volume from the hydrodynamical simulation, this approach modestly overestimates the effect of UVB fluctuations on the IGM opacity because the anticorrelation between the UVB and the large-scale density field (Davies & Furlanetto 2016; Davies et al. 2018a) is lost. While small-scale correlations between the UVB and the density field may alter the relationship between $\mathrm{Ly}\alpha$ and $\mathrm{Ly}\beta$ forest opacity (Oh & Furlanetto 2005), the UVB fluctuations in the Davies et al. (2018a) model manifest on scales comparable to the size of our 40 cMpc bins, so we do not expect a substantial effect.

In order to model temperature fluctuations, we use a new hybrid method introduced in Oñorbe et al. (2019) that uses a small set of phenomenological input parameters, and combines a seminumerical reionization model (21cmFAST; Mesinger et al. 2011) to solve for the topology of reionization, as well as an approximate model of how reionization heats the IGM, with the cosmological hydrodynamical code Nyx. Instead of applying a uniform UVB to the whole simulation box, which was the procedure in the original Nyx simulations, the UVB is now applied to denser regions in the simulations at earlier reionization redshifts z_reion, whereas the underdense voids will be reionized at later times. At the same time, the gas temperature is heated up. The topology of the different reionization times z_reion are extracted from the seminumerical model.

This guarantees that the temperature evolution of the inhomogeneous reionization process and the small-scale structure of the diffuse gas of the IGM is resolved and captured self-consistently. In this work, we used the IR-A reionization model presented in Oñorbe et al. (2019), which is an extended reionization model consistent with the observations of the cosmic microwave background (CMB) by Planck Collaboration et al. (2018) with a median reionization redshift of z_reion = 7.75, a duration of the reionization process of Δz_reion = 4.82, and a heat injection due to reionization of ΔT = 20,000 K. We apply this model to a Nyx simulation with L = 40 Mpc h⁻¹ box side and 2048³ resolution elements. The full thermal evolution of this model at the redshifts of interest for this work can be seen in Figures 5 and 6 in Oñorbe et al. (2019).

We will see in Section 5 that neither of the two models described above provides a good fit to our opacity measurements in both the $\mathrm{Ly}\alpha$ and the $\mathrm{Ly}\beta$ forests. Because the ratio of these opacities is sensitive to the temperature–density relation of the IGM, we will also compare our measurements to models with a homogeneous IGM but different slopes γ − 1 of the temperature–density relation (see Equation (3)).

To this end, we simply ignore the temperatures of the extracted Nyx skewers and calculate new temperatures in post-processing for the gas densities along each skewer by means of Equation (3), assuming T₀ = 10,000 K. We then take the velocity skewers and map again the modified real space into redshift space to calculate mean fluxes and opacities.

We will compare our observations to the reionization models in three redshift bins with ${\rm{\Delta }}z=0.2$ at the central redshifts of z = 5.6, z = 5.8, and z = 6.0. At each redshift, we extract N_skew = 2000 skewers from the simulation box of the same size as our measurements, i.e., 40 cMpc. It is well known that the Lyα opacity depends on the unknown amplitude of the UVB Γ_{H i}, so we rescale the optical depth in the $\mathrm{Ly}\alpha$ forest to match the observations and determine this unknown parameter. For each skewer, we compute the mean flux, or equivalently, the effective optical depth according to

$$\begin{eqnarray}\begin{array}{rcl}\langle {F}^{\mathrm{Ly}\alpha }\rangle & = & \langle \exp \left[-{\tau }_{i}^{\mathrm{Ly}\alpha }\right]\rangle \\ & & =\langle \exp \left[-{A}_{0}\times {\tau }_{i}^{\mathrm{Ly}\alpha ,\mathrm{unscaled}}\right]\rangle ,\end{array}\end{eqnarray} \tag{ 6 }$$

where the angle brackets denote the average in the 40 cMpc bins. In the last equality, we introduce the scaling factor A₀, which provides the lever arm for tuning the photoionization rate of the UVB Γ_{H i} to match the observations.

Studies of the Lyα forest at lower redshift, i.e., 2 < z < 5, typically tune A₀ to match the mean flux. Because this study focuses on comparing the distribution of mean fluxes (effective optical depths) in 40 cMpc bins, we instead determine the value of A₀ by requiring the 25th percentile of the cumulative distribution of the mean flux to match the data. Specifically, for the data in any redshift bin, we can solve the equation $P(\gt \langle {F}^{\mathrm{Ly}\alpha }\rangle )=0.25$ for $\langle {F}^{\mathrm{Ly}\alpha }\rangle$. Our procedure for setting the unknown Γ_{H i} for any given model is then to vary A₀ until the model matches the $\langle {F}^{\mathrm{Ly}\alpha }\rangle$ determined from the data for that redshift bin. From the cumulative distribution of values in Table 3 in Appendix B, we measure $\langle {F}^{\mathrm{Ly}\alpha }{\rangle }_{25\mathrm{th}}(z=5.6)\approx 0.0713$, $\langle {F}^{\mathrm{Ly}\alpha }{\rangle }_{25\mathrm{th}}(z=5.8)\approx 0.0363$, and $\langle {F}^{\mathrm{Ly}\alpha }{\rangle }_{25\mathrm{th}}(z=6.0)\approx 0.0144$. Note that the choice of the 25th percentile is somewhat arbitrary, and we could have also chosen to match the observed mean or median flux, which would not change any of the conclusions of this work.

We then obtain the pure $\mathrm{Ly}\beta$ optical depth ${\tau }_{i}^{\mathrm{Ly}\beta }$ at each pixel from ${\tau }_{i}^{\mathrm{Ly}\alpha }$ by scaling according to the oscillator strengths of their transitions and the respective wavelengths, i.e.,

$$\begin{eqnarray}&&{\tau }_{i}^{\mathrm{Ly}\beta }=\displaystyle \frac{{f}_{\mathrm{Ly}\beta }}{{f}_{\mathrm{Ly}\alpha }}\,\displaystyle \frac{{\lambda }_{\mathrm{Ly}\beta }}{{\lambda }_{\mathrm{Ly}\alpha }}\,{\tau }_{i}^{\mathrm{Ly}\alpha }\end{eqnarray} \tag{ 7 }$$

with f_Lyα = 0.41641 and f_Lyβ = 0.079142 (Table 4 in Wiese & Fuhr (2009)). This scaling is correct, provided that there are no damping wing effects due to the presence of very dense or highly neutral gas. Given that our analysis is focused on the low-density IGM at z ≳ 5.5 probed by the $\mathrm{Ly}\alpha$ and $\mathrm{Ly}\beta$ forests at typical locations of the universe, this is a good approximation.

However, the observed $\mathrm{Ly}\beta$ effective optical depth ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\beta ,\mathrm{obs}}$ will be higher than the pure $\mathrm{Ly}\beta$ optical depth, due to the additional absorption from the $\mathrm{Ly}\alpha$ forest at the foreground redshift. Thus, we create additional skewers along different lines of sight with $\mathrm{Ly}\alpha$ forest absorption at the foreground redshift z_fg, i.e., ${\tau }_{i}^{\mathrm{Ly}\alpha ,\mathrm{fg}}$, using lower-redshift outputs from the Nyx simulation. These lower-redshift skewers are also rescaled, but to match the mean flux in order to be consistent with published measurements. We use the fitting formula presented in Oñorbe et al. (2017a) to the mean flux measurements by Faucher-Giguère et al. (2008b), i.e., $\langle {F}^{\mathrm{Ly}\alpha }\rangle ({z}_{\mathrm{fg}}=4.57)\approx 0.2371$, $\langle {F}^{\mathrm{Ly}\alpha }\rangle ({z}_{\mathrm{fg}}=4.74)\approx 0.1945$, and $\langle {F}^{\mathrm{Ly}\alpha }\rangle ({z}_{\mathrm{fg}}=4.91)\approx 0.1560$ for the redshift bins at z = 5.6, z = 5.8, and z = 6.0, respectively.

For each skewer of the high-redshift $\mathrm{Ly}\beta$ forest, we now draw a random $\mathrm{Ly}\alpha$ forest segment from a simulation output at the corresponding foreground redshift. In this way, we ensure that we account for the full distribution of foreground $\mathrm{Ly}\alpha$ opacities, and do not underestimate the scatter between different foreground sightlines. At each pixel, we then sum the pure $\mathrm{Ly}\beta$ forest optical depth with the foreground $\mathrm{Ly}\alpha$ optical depth to obtain the flux one would observe in the $\mathrm{Ly}\beta$ forest, i.e.,

$$\begin{eqnarray}&&{F}_{i}^{\mathrm{Ly}\beta ,\mathrm{obs}}=\exp \left[-({\tau }_{i}^{\mathrm{Ly}\beta }+{\tau }_{i}^{\mathrm{Ly}\alpha ,\mathrm{fg}})\right].\end{eqnarray} \tag{ 8 }$$

We calculate the observed effective $\mathrm{Ly}\beta$ optical depth by averaging the flux of all pixels in the spectral bin, i.e., ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\beta ,\mathrm{obs}}=-\mathrm{ln}\,\langle {F}_{i}^{\mathrm{Ly}\beta ,\mathrm{obs}}\rangle$.

Figure 3 shows the various distributions of effective optical depths in an example case for the uniform reionization model with γ = 1.5 at z = 6, but all other models look qualitatively similar. The dark blue histogram shows the distribution of ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\alpha }$ at z = 6, whereas the green histogram depicts the distribution of ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\alpha ,\mathrm{fg}}$ along different skewers at the foreground redshift z_fg ≈ 4.9. The yellow histogram shows the pure ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\beta }$ values obtained from rescaling according to Equation (6), while the light blue histogram represents the distribution of observed ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\beta \,\mathrm{obs}}$ values after accounting for the distribution of foreground $\mathrm{Ly}\alpha$ opacities.

Zoom In Zoom Out Reset image size

The spectral noise in the data influences the opacity measurements. This is because whenever the average flux is detected with less than 2σ significance, we adopt a lower limit on the effective optical depth measurement (see Section 2.3). To mimic the effects of the noisy data, for each mean flux value $\langle {F}_{\mathrm{sim}}\rangle$ from the simulated skewers, we randomly draw an uncertainty ${\sigma }_{\langle {F}_{\mathrm{data}}\rangle }$ from the measurements from our data set. Note that the noise level formally depends on the signal, but in the limit where the objects are fainter than the sky, which is almost always the case for the highly absorbed forests in high-redshift quasars, the noise is dominated by the sky background and read noise rather than the object photon noise, so this is a good approximation.

For this mock data with forward-modeled noise, we then apply the same criterion as for the real data and adopt a measurement only if $\langle {F}_{\mathrm{sim}}\rangle \geqslant 2{\sigma }_{\langle {F}_{\mathrm{data}}\rangle }$, or a lower limit at $2{\sigma }_{\langle {F}_{\mathrm{data}}\rangle }$ otherwise. This reduces the sensitivity to very high opacities in the simulated skewers, which we do not have in noisy data.

We compare our co-spatial measurements of ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\alpha }$ and ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\beta ,\mathrm{obs}}$ to predictions from the various models of the post-reionization IGM in three different redshift bins centered around z = 5.6, 5.8, and 6.0 in Figure 4. The contours in the middle and bottom panels show the predicted parameter space from models with spatial fluctuations in the UVB and the temperature field, respectively. Note that we did not yet forward model the spectral noise in the data to the simulations here, because we are overlaying the data with error bars and indicating limits. Thus, for this qualitative comparison, forward modeling the noise would amount to effectively double counting the noise.

Zoom In Zoom Out Reset image size

While our measurements fall within the predicted parameter space in the lowest-redshift bin, it is evident that the observed τ_eff measurements in the two higher-redshift bins are not well confined to the predicted parameter space. The measured ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\beta ,\mathrm{obs}}$ values show a larger scatter than expected by the models at z ∼ 5.8, as well as an offset in the mean at z ∼ 6.

Another way of illustrating this comparison is by means of the cumulative distribution function (CDF). In Figures 5 and 6, we present the CDFs of our opacity measurements as the red curves, as well as the predicted distributions from the post-reionization models including a fluctuating UVB and a fluctuating temperature field, respectively, as dotted black curves. Note that, whereas the measurements presented in Figure 4 are restricted to the spectral regions where we have co-spatial measurements of ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\alpha }$ and ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\beta ,\mathrm{obs}}$, the cumulative histograms for ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\alpha }$ in Figures 5 and 6 show all measurements from the full $\mathrm{Ly}\alpha$ forests along all quasar sightlines in our data set.⁸

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Both CDFs show clearly that, when fine-tuning the reionization models to match the observed transmission in the $\mathrm{Ly}\alpha$ forest, the effective optical depth in the $\mathrm{Ly}\beta$ forest is underestimated in the highest-redshift bin at z ∼ 6. Additionally, as noted before, we find a large increase in the scatter of the measurements, most notably at z ∼ 5.8, compared to the predictions from the models.

Note that, in the two highest-redshift bins in Figure 5, the difference between the observed mean Lyβ forest optical depth and the reionization model with a fluctuating UVB is Δτ^Lyβ_eff ≳ 1, which corresponds to a factor of ≳2.5 in the mean flux. Hence, systematic uncertainties in the data that might arise from poor quasar continuum estimates or issues in the data reduction process (see also Appendix A), would have to change the observed averaged flux by ≳250% to account for this offset. This seems highly unlikely, as there is no obvious reason why such errors, if they were present, would be so asymmetric as to cause the flux in the $\mathrm{Ly}\beta$ region to be systematically lower by this large factor. Furthermore, if such systematics were present, they would presumably also impact the $\mathrm{Ly}\alpha$ measurements. However, a comparison of independent measurements of the distribution of ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\alpha }$ between Bosman et al. (2018) and Eilers et al. (2018) (see Figure 7 in Eilers et al. (2018)) shows no evidence for systematic offsets of ${\rm{\Delta }}{\tau }_{\mathrm{eff}}^{\mathrm{Ly}\beta }\gtrsim 1$. Note that the difference in Figure 6 between the mean $\mathrm{Ly}\beta$ optical depth of the fluctuating temperature model predictions and the measurements is smaller, but still significant.

We know that the ratio of $\mathrm{Ly}\alpha$ and $\mathrm{Ly}\beta$ optical depths is sensitive to the ionization as well as the thermal state of the IGM (Oh & Furlanetto 2005; Trac et al. 2008; Furlanetto & Oh 2009). Thus, by varying the slope of the temperature–density relation of the intergalactic gas and adjusting it to be more isothermal (or inverted), we expect an increase in the predicted ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\beta ,\mathrm{obs}}$.

To this end, we also compare our observations to reionization models with a uniform IGM but different slopes of their temperature–density relations, i.e., the fiducial value of γ = 1.5, an isothermal temperature–density relation with γ = 1.0, and a highly inverted slope with γ = 0.0. This is shown in the top panels of Figure 4, as well as in Figure 7.

Zoom In Zoom Out Reset image size

The comparison shows that the model with a highly inverted temperature–density relation, i.e., γ = 0.0, predicts a higher $\mathrm{Ly}\beta$ opacity and seems to match the observations in the $\mathrm{Ly}\beta$ forest in the highest-redshift bin better than the previously discussed inhomogeneous reionization models, although it is still not fully capturing the whole parameter space of the observations. Additionally, none of the uniform reionization models show enough spatial fluctuations to predict the scatter in the $\mathrm{Ly}\alpha$ opacity.

How likely is it that the IGM follows an inverted temperature–density relation  at z ∼ 6? It has been argued that reionization events dramatically alter the thermal structure of the IGM, increasing the IGM temperature by at least an order of magnitude up to T = 25,000–30,000 K (Furlanetto & Oh 2009; D'Aloisio et al. 2019). Several hydrodynamical simulations that self-consistently model the reionization process have shown that a flat or inverted temperature–density relation arises naturally during reionization events, due to the different number densities of ionizing sources in over- and underdense regions (e.g., Trac et al. 2008; Finlator et al. 2018). The gas in dense regions ionizes first, due to the higher number of ionizing sources, and thus has more time to cool down, whereas low-density voids are ionized last and hence contain the hottest gas at the end of reionization, which leads to an at least partially inverted temperature–density relation of the IGM (Furlanetto & Oh 2009; Lidz & Malloy 2014; D'Aloisio et al. 2015, 2019). Oñorbe et al. (2019) modeled this process in detail, and their model (see Section 4.2) has an average slope of the temperature–density relation of γ ≈ 1.1, i.e., roughly isothermal, at z ∼ 6, but with a very large scatter (see their Figure 5), which arises due to a superposition of different heat injections and subsequent cooling at different times. Thus, their model predicts that at least some regions of the IGM do indeed follow an inverted temperature–density relation. This could explain why the model with a fluctuating temperature field (Figure 6) matches the ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\beta ,\mathrm{obs}}$ observations better than the model with fluctuations in the UVB (Figure 5).

Although Furlanetto & Oh (2009) claim that rapid adiabatic and Compton cooling will quickly erase any inversion of the thermal state of the IGM after the reionization event, we might still be observing the end stages of this process at z ∼ 5.5–6.0. This result complies with the conclusions from Trac et al. (2008). Contrary to Furlanetto & Oh (2009), they argue that the inversion of the IGM endures long past the reionization process (see also Finlator et al. 2018). Thus, our results provide support for a late and extended EoR, which is also in good agreement with several recent studies (Bosman et al. 2018; Eilers et al. 2018; Planck Collaboration et al. 2018; Kulkarni et al. 2019).

In previous studies, comparisons of τ_eff measurements to reionization models have been performed assuming an infinite S/N ratio of the simulated skewers. These noise-free predictions for the CDFs of τ_eff correspond to the dotted curves in Figures 5–7. However, in order to conduct a proper comparison between the noisy data of any realistic data set to these models, it necessary to model the spectral noise. Thus, as discussed in Section 4.5, we forward model our observations by adding spectral noise to the simulated skewers.

The models including the forward-modeled noise are shown as the dashed curves in Figures 5–7, which predict less scatter at high optical depths. As a result, the models with a spatially inhomogeneous reionization process that had been claimed previously to fully reproduce the scatter in the $\mathrm{Ly}\alpha$ opacity observations (Figures 5 and 6) do not contain enough spatial fluctuations once we account for the spectral noise in the data. This discrepancy is most obvious in the top right panel of each figure for ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\alpha }$ measurements at z = 6.0.

The reason for this discrepancy is quite simple: In any realistic data set, there will be a distribution in the S/N ratios of the spectra due to the various object magnitudes and exposure times. This noise manifests as noisy τ_eff measurements. In order to detect a very high τ_eff, one must not only encounter a rare fluctuation in the IGM along a quasar sightline, but this fluctuation must also be probed by a sufficiently bright quasar or deep enough spectrum to measure a high τ_eff. Thus, large τ_eff fluctuations in the IGM need to be more common than one naively expects from infinite S/N models.

We conclude that the forward modeling of spectral noise is essential to conduct a fair data model comparison, and our results suggest that stronger spatial inhomogeneities than previously assumed are needed to fully capture the observations. At lower opacities, the difference between perfect and noisy models is less significant, as most measurements are actually high-S/N detections rather than lower limits.