New Constraints on IGM Thermal Evolution from the Lyα Forest Power Spectrum

To compare to existing measurements, which did not apply masking of spectral regions, but instead treated metal contamination statistically by comparing to lower redshift data where most metals are outside the Lyα forest, we compute the power spectrum based on ∼50,000 noiseless, high-resolution skewers from our simulation. We will refer to this as the "perfect model."

However, to account for the window function introduced on the power spectrum by masking parts of the data fully, when comparing to our measurement from Paper I, we compute the power spectrum based on the skewers for each combination of parameters applying the full forward modeling technique described in Paper I to our hydrodynamical simulations. Henceforth we will call this the "forward model." This technique consists of several steps of post-processing the hydrodynamical simulation outputs followed by a power spectrum computation in the same way as for the data. To forward model an individual quasar spectrum we first merge randomly selected skewers (without repetition) to cover the same path length as the data, then convolve the spectra with a Gaussian smoothing kernel reducing the resolution of the models to match that of the data, rebin the models onto the pixels of the observed spectra, and add noise drawn from a Gaussian distribution for each individual pixel with a standard deviation equal to the 1σ uncertainty of the corresponding quasar spectral pixel reported by the data reduction pipelines. Finally and most importantly, we mask the forward modeled spectrum in exactly the same way as the data to account for the windowing effects resulting from gaps in the data and our metal masking procedure. We then compute the power spectrum by utilizing ∼50,000 skewers from our simulation (see Paper I for a more detailed description of the individual steps). Note that while the full forward modeling of noise and resolution might not be completely necessary as they have been corrected in the measurement (and are corrected in the same way inside the forward modeling procedure as well), there might be subtle effects on the masking correction. We therefore want to make the model spectra as similar to the data as possible. Note that this does not change our model precision, which is dominated by data set size rather than noise or resolution.

To perform a fit to the data and infer the thermal state at a particular redshift we need to be able to compute power spectra on a continuous range of parameters. Therefore, we need to interpolate between the discrete and sparse outputs of the THERMAL grid. To perform this task we follow the emulation approach of Heitmann et al. (2006) and Habib et al. (2007). For details, we refer the reader to their papers (and references therein) as well as Paper I; in the following we summarize the main steps of the approach. First, we decompose the simulated logarithmic power spectra onto a principal component analysis (PCA) basis. We save the PCA vectors as well as the coefficients ${A}_{i}({{\boldsymbol{\Theta }}}_{j})$ at each thermal model location ${{\boldsymbol{\Theta }}}_{j}$. We then use a Gaussian process (GP) to interpolate the coefficients ${A}_{i}({{\boldsymbol{\Theta }}}_{j})$ onto any arbitrary location in parameter space ${\boldsymbol{\Theta }}$. Taking the dot product of the PCA vectors with these interpolated coefficients then gives the power spectrum evaluated at any parameter location.

We thus calculate a GP for each principal component coefficient (using George, see Ambikasaran et al. 2016) using a squared exponential kernel plus an additional white noise contribution

$$\begin{eqnarray}&&K({{\boldsymbol{\Theta }}}_{1},{{\boldsymbol{\Theta }}}_{2})=\exp (-0.5({{\boldsymbol{\Theta }}}_{1}-{{\boldsymbol{\Theta }}}_{2}){C}_{l}^{-1}({{\boldsymbol{\Theta }}}_{1}-{{\boldsymbol{\Theta }}}_{2}))+{\sigma }_{n}{\delta }_{\mathrm{ij}}\end{eqnarray} \tag{ 3 }$$

for parameter values ${{\boldsymbol{\Theta }}}_{i}$, a chosen distance metric C_l (which is defined by a smoothing length l for each parameter, i.e., its diagonal) and a noise contribution σ_n (for an in-depth introduction to GP techniques, see Rasmussen & Williams 2005).

As the hydrodynamical grid consists of far less models (∼50)¹¹ than the previous DM-based grid (∼500) used in Paper I, we must be more careful about the interpolation errors resulting from our emulation procedure. Instead of just using a kernel with a fixed hand-tuned smoothing length, which was our approach in Paper I, we additionally optimized our kernel parameters by maximizing GP-likelihood using the scipy.optimize (Jones et al. 2001) package and the so-called L-BFGS-B (Zhu et al. 1997) method.¹² We then performed the analysis using the optimal smoothing lengths l and noise σ_n for the kernel for each GP emulator.

We estimate the emulation uncertainties using a cross-validation scheme to propagate interpolation errors. To do this we generate the emulator, but leave one simulation out of the training set.¹³ We denote emulators with a model (defined by parameters ${\boldsymbol{\Theta }}$) left out as $\mathrm{emu}\setminus {\boldsymbol{\Theta }}$. We then compare the actual models (with power P_model) for this simulation to the emulator (with power ${P}_{\mathrm{emu}\setminus {\boldsymbol{\Theta }}}$) at the parameters ${\boldsymbol{\Theta }}$ of this model:

$$\begin{eqnarray}&&{\rm{\Delta }}{P}_{\mathrm{emu}}({\boldsymbol{k}},{\boldsymbol{\Theta }})={P}_{\mathrm{model},{\boldsymbol{\Theta }}}({\boldsymbol{k}})-{P}_{\mathrm{emu}\setminus {\boldsymbol{\Theta }}}({\boldsymbol{k}},{\boldsymbol{\Theta }}).\end{eqnarray} \tag{ 4 }$$

We show the accuracy of the emulation in Figure 3. This shows quantiles of the deviations ${\rm{\Delta }}{P}_{\mathrm{emu}}$ from the true underlying model inside our cross-validation sample. We see that for most models in our parameter space the emulator works to better than 1%. However, emulation uncertainty can increase to the 5% level (with a preference for underestimation at $k\gt 0.06\,{\rm{s}}\ {\mathrm{km}}^{-1}$) for some models. As the uncertainty in our power spectrum measurements is ∼2% (for the 68% quantile) on large scales (k ≲ 0.01 $\,{\rm{s}}\ {\mathrm{km}}^{-1}$) and ≳5% on smaller scales, measurement errors are much larger than these interpolation errors. Nevertheless, we opted to add the covariance matrix for the interpolation process to our likelihood. This covariance matrix can be obtained by performing:

$$\begin{eqnarray}&&{C}_{\mathrm{emu},\mathrm{ij}}=\langle {{\rm{\Delta }}}_{\mathrm{emu}}({{\boldsymbol{k}}}_{{\rm{i}}},{\boldsymbol{\Theta }}){{\rm{\Delta }}}_{\mathrm{emu}}({{\boldsymbol{k}}}_{{\rm{j}}},{\boldsymbol{\Theta }})\rangle \end{eqnarray} \tag{ 5 }$$

with the average performed over all possible combinations of model parameters inside our grid for each redshift bin.

Zoom In Zoom Out Reset image size

Due to the variety of thermal histories in the THERMAL suite some simulations can have extremely close values of their thermal parameters at some specific redshifts. In order to avoid possible problems in the emulator due to this issue we removed models from the THERMAL grid that did not satisfy a distance threshold¹⁴ and are left with 45–65 models per redshift.

We perform a Bayesian MCMC analysis on the power spectrum data at each individual redshift using the emcee package (Foreman-Mackey et al. 2013) based on the affine invariant sampling technique (Goodman & Weare 2010) and assuming the multivariate Gaussian likelihood:

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal L } & \equiv & P(\mathrm{data}| \mathrm{model})\\ & \propto & \displaystyle \prod _{\mathrm{data}\,\mathrm{sets}}\displaystyle \frac{1}{\sqrt{\det (C)}}\exp \left(-\displaystyle \frac{{{\boldsymbol{\Delta }}}^{{\rm{T}}}{C}^{-1}{\boldsymbol{\Delta }}}{2}\right)\\ {\boldsymbol{\Delta }} & = & {{\boldsymbol{P}}}_{\mathrm{data}}-{{\boldsymbol{P}}}_{\mathrm{emu}}\\ C & = & {C}_{\mathrm{data}}+{C}_{\mathrm{emu}}.\end{array}\end{eqnarray} \tag{ 6 }$$

with C_emu being the covariance of the interpolation procedure and C_data being the covariance of an individual measurement. For these covariances we use published values if available. For our own data set from Paper I as well as the Viel et al. (2013) data set, we used the published uncertainties (i.e., the diagonal covariance elements) and combined them with the correlation matrix of the model closest in parameter space to obtain an estimate of the covariance, i.e., we perform nearest neighbor interpolation between covariance matrices obtained at every point (see Paper I for details on this approach).

Our modeling so far depends on four parameters, T₀ and γ describing the thermal state, λ_P for the pressure smoothing depending on the full thermal history, and $\bar{F}$ for the mean transmission that corresponds to a given UVB amplitude. There is, however, one additional parameter that we input in our models for each data set¹⁵ to generate the observed correlation between $\mathrm{Si}\,{\rm{III}}$ and Lyα (see McDonald et al. 2006; Palanque-Delabrouille et al. 2013). Finally, because of significant uncertainties in the resolution of the XQ-100 data (see the detailed discussion in Appendix B of Paper I), we also marginalize over the resolution of the XQ-100 measurement whenever we use this data, giving us another parameter. Therefore, we have a total of five (in the case of high-resolution data only) to eight (in the case of fitting 3 data sets of which one comes from XQ-100) parameters. We assume flat priors on $\mathrm{log}{T}_{0},\mathrm{log}{\lambda }_{{\rm{P}}},\gamma$. We now go into further detail about the modeling and assumptions for the other parameters.

We add $\mathrm{Si}\,{\rm{III}}$ correlations to the model analytically by multiplying the model power spectrum with an oscillating signal as correlations inside a spectrum correspond to oscillations of the corresponding power spectrum:

$$\begin{eqnarray}&&{P}_{\mathrm{tot}}=(1+{a}_{\mathrm{Si}\,{\rm{III}}}^{2}+2a\cos (k\,{\rm{\Delta }}v)){P}_{{\rm{H}}{\rm{I}}}\end{eqnarray} \tag{ 7 }$$

with ${a}_{\mathrm{Si}{\rm{III}}}$ being a free nuisance parameter for the strength of the correlation. In previous works this was typically expressed as ${a}_{\mathrm{Si}{\rm{III}}}={f}_{\mathrm{Si}{\rm{III}}}/(1-\bar{F})$ with ${f}_{\mathrm{Si}{\rm{III}}}$ being a redshift independent quantity that was fit using the entire data set. We adopt this same parametrization but opt to fit for a unique value of ${f}_{\mathrm{Si}{\rm{III}}}$ at each redshift and for each data set because of the different metal treatment in the data sets and as we do not perform a joint fit of different redshifts here. We assume a flat prior on each ${f}_{\mathrm{Si}{\rm{III}}}$ and demand correlations to be positive.

We modeled the resolution of the XSHOOTER spectrograph R_new by multiplying the measured XQ-100 power spectrum with the resolution dependent part of the window function:

$$\begin{eqnarray}&&{W}_{R}(k,R)=\exp \left(-\displaystyle \frac{1}{2}{\left({kR}\right)}^{2}\right),\end{eqnarray} \tag{ 8 }$$

using the resolutions quoted in Iršič et al. (2017a) and dividing by ${W}_{R}(k,{R}_{\mathrm{new}})$. Note that the resolution of the instrument depends on two different factors: the resolution for a fully illuminated slit (or "slit resolution") and the seeing, which gives rise to higher spectral resolution if smaller than the slit size. We assume two limits for the resolving power of the XQ-100 data set. The lower limit assumes the slit resolutions quoted in the XQ-100 data release paper (López et al. 2016) as well as a fully illuminated slit (leading to ${R}_{\mathrm{UVB}}=4350$ for the UVB arm of the instrument, ${R}_{\mathrm{VIS}}=7410$ for the VIS arm¹⁶ ). The upper limit assumes a seeing of 065 (smaller than the slit) and higher values for the slit resolution¹⁷ (leading to R_UVB = 8230 and R_VIS = 12184). We assume a flat prior between these two limits. As z = 3.6 is using both spectral arms we use the lowest and highest of the four resolution values above as the limits here. Note that this choice of priors on spectroscopic resolution is an extremely conservative choice that will significantly weaken the constraints that can be obtained from this XQ-100 data set. This is most acute in the UVB arm because of its intrinsically lower resolution. A more careful analysis of the XQ-100 resolution would allow us to adopt a far stronger prior on these values, which would increase the precision of constraints deduced from power spectra measured from such moderate resolution spectra.

Note that most previous measurements (exceptions to this are, e.g., Lidz et al. 2010; Iršič et al. 2017b) of the IGM's thermal state did not attempt to marginalize over the uncertainty in the mean flux estimate. Instead, typically simulations that match the mean flux of the data assuming perfect knowledge of this quantity are used (e.g., in Voigt profile fitting or curvature analyses). For $\bar{F}$ we used both a flat prior (corresponding to performing a joint measurement of $\bar{F}$ and the thermal state) and a Gaussian-shaped prior. For the Gaussian prior we assumed a mean based on the fit by Oñorbe et al. (2017a) to a compilation of recent measurements (Fan et al. 2006; Kirkman et al. 2007; Faucher-Giguère et al. 2008b; Becker et al. 2013) and a standard deviation based on the uncertainties for the most recent measurements at z ≤ 4.0: Becker et al. (2013) for $2.2\leqslant z\leqslant 4.0$, Faucher-Giguère et al. (2008a) for z = 2.0, and Kirkman et al. (2005) for z = 1.8. For z ≥ 4.2 we use ${\sigma }_{\bar{F}}=0.03$, which is loosely based on the discrepancy between Fan et al. (2006) for z ≥ 4.6 and the measurements by Becker et al. (2011) in the range 4.1 ≤ z ≤ 4.7 (see also Bosman et al. 2018; Eilers et al. 2018, for more recent mean flux measurements that are discrepant by a similar amount for 5.0 ≤ z ≤ 5.4).

To avoid extrapolating from our model grid we additionally require that all thermal parameters lie inside the convex hull of our model grid (see Figure 2), i.e., the smallest convex shape, including all THERMAL grid points. The convex hull is evaluated numerically by triangulating the model grid (using scipy.spatial.Delaunay) and for each MCMC sample we test whether it is inside the triangulation when evaluating the prior. Otherwise the prior is set to zero. To see the effective prior resulting from only using this non-rectangular region where we have models, we performed an MCMC run assuming a completely uninformative data set, i.e., using only the priors in our fit and a constant likelihood. The results of this procedure are shown in Figure 4 for z = 2.8. In some contours, e.g., T₀ and λ_P, we can see that the parameters are highly correlated already since our grid is non-rectangular. We argue, however, that these correlations are physically motivated as models perpendicular to these correlations are hard to produce (see Section 3) and that this behavior actually constitutes prior information for our analysis.

Zoom In Zoom Out Reset image size

We performed fits of the parameters governing the thermal state using combinations of all data sets discussed in Section 2 in 16 individual redshift bins with 1.8 < z < 5.4, where we used a bin size Δz = 0.2 for z ≤ 4.2 and Δz = 0.4 for z ≥ 4.6.

The power spectra of each data set are summarized and compared to models based on our posterior MCMC chains in Figure 5. Note that for visualization purposes we only compare window function, $\mathrm{Si}\,{\rm{III}}$ correlation, and resolution corrected data to the perfect model. The window function due to masking was taken out of the UVES/HIRES data by multiplying measurement points with the median ${P}_{\mathrm{emu},\mathrm{perfect}}/{P}_{\mathrm{emu},\mathrm{forward}}$ for our MCMC chain and propagating its uncertainties using Gaussian error propagation for each individual mode (see Paper I for a more detailed description of this process).¹⁸ Analogously, we rescaled the XQ-100 power to use the "best-fit" resolution correction, i.e., we renormalize with ${W}_{R}(k,{R}_{\mathrm{new}})/{W}_{R}(k,R)$ (see Equation (8)) from the posterior and removed $\mathrm{Si}\,{\rm{III}}$ correlations from the data applying Equation (7). We can see that satisfactory fits have been achieved at all redshifts.

Zoom In Zoom Out Reset image size

In Figure 6, we further illustrate the posterior distribution is inferred via our MCMC at z = 2.8 with a so-called "corner plot." We can see that the data strongly constrains all parameters (e.g., compare to Figure 4 or the blue curves in the 1d histograms, for which the likelihood is assumed to be completely uninformative). The most important feature we see is that there are strong degeneracies between some parameters, e.g., the diagonal contours between permutations of T₀, γ, and $\bar{F}$. Note that the strong correlation between T₀ and γ is well understood and results from the IGM not probing the mean density, but instead mild overdensities at these redshifts (see, e.g., Lidz et al. 2010; Becker et al. 2011). We also infer a low mean transmitted flux $\bar{F}=0.69\pm 0.01$ compared to the Becker et al. (2013) measurement of $\bar{F}=0.727\pm 0.009$ (green band). It is interesting to note that this low value however agrees well with the joint constraint on mean transmission evolution by Palanque-Delabrouille et al. (2015) obtained from the BOSS power spectrum yielding $A=0.0028\pm 0.0002,\eta =3.67\pm 0.02$ for $\bar{F}{(z)=\exp (-A(1+z)}^{\eta })$ resulting in $\bar{F}(z=2.8)\approx 0.687\,\pm 0.020$. Note that the data set used in this analysis overlaps with the one we used here, but the simulations and inference procedure are independent and our analysis has additional higher resolution data available. Independent of the BOSS data, we also obtain similarly low $\bar{F}$ values when performing fits on the high-resolution data from Paper I alone.

Zoom In Zoom Out Reset image size

Additionally, the posterior distribution for γ shows a clear preference for values γ ≈ 2.1, far above the expected value of ∼1.6 for IGM gas in photoionization equilibrium long after reionization events (Hui & Gnedin 1997; McQuinn & Upton Sanderbeck 2016) Note again that there is a strong anticorrelation between γ and $\bar{F}$, so while our analysis prefers a high value of γ and a low value of $\bar{F}$, this is a movement along the degeneracy direction. We will further discuss this issue in Section 5.2.

The redshift evolution of individual parameters, determined from the 1d-marginalized posteriors, is illustrated in Figure 7. For 3.0 ≤ z ≤ 3.4 we also performed fits including the XQ-100 data, and fully marginalized over our lack of knowledge of the exact spectroscopic resolution (see discussion in Section 4.4). As including this data set did not significantly change our results, we decided to leave those points off the plot for clarity. Numerical values for the marginalized parameters are tabulated in Table 3 (see the Appendix for further details).

Zoom In Zoom Out Reset image size

There are several noteworthy features in Figure 7. First, the disagreement that we saw at z = 2.8 between our inferred value of $\bar{F}$ and recent measurements is also present at all other redshifts z < 3 (green and blue data points compared to the pink-shaded region in the lower panel). At the same time γ reaches very high values in the same redshift range. Also, T₀ drops strongly from z = 3.0 to lower redshifts, but due to the degeneracies between T₀, γ, and $\bar{F}$, these measurements are all strongly correlated and this effect is therefore expected. Note that these trends—high γ, low $\bar{F}$, and low T₀—persists if we fit the high-resolution data alone, as the BOSS data alone do not individually constrain all of these parameters due to the lack of high k-modes (resulting from limited spectral resolution).

Second, for z ≥ 3 we can see that γ shows little evolution and the mean transmitted flux $\bar{F}$ is consistent with the Oñorbe et al. (2017a) fit to recent measurements. We can also see that T₀ increases from ≈5100 K at z = 5.0 to $\approx {\rm{15,000}}\,\,{\rm{K}}$ at z = 3.4. This rise could be explained by the onset of $\mathrm{He}\,{\rm{II}}$ reionization, which we discuss in more detail in Section 5.5 where we compare our inferred parameter values to models of IGM thermal history that treat reionization heating.

In summary, we can see that the power spectrum analyzed here can in principle achieve high-precision constraints on IGM thermal parameters and the mean transmission, but the high values of γ ≃ 2 inferred at z < 3 and concomitant discrepancies between our inferred mean flux and the Becker et al. (2013) measurements might indicate systematics in our procedure. We consider this issue in detail in the next section.

In the previous section we found low values of $\bar{F}$ compared to Becker et al. (2013) and possibly unphysically high values of γ. While both parameters are degenerate and the degeneracy direction matches with our discrepancy this might point toward some problem within the analysis. To investigate this scenario we want to isolate the change in the power spectrum when moving along the degeneracy direction of our posterior distributions. Due to the dimensionality of the parameter space and correlations between different parameters this cannot be achieved by simple cuts along a parameter direction. Therefore, we designed the following procedure to generate model curves tracking the degeneracy direction for different values of γ (also see the illustration in Figure 8):

• 1.  
  We take the posterior of our MCMC analysis (i.e., the Markov chain) and define bins such that the median of γ inside a bin is equal to a desired quantile of the marginalized γ distribution (which are chosen to be equivalent to ±1σ, ±2σ). These bins are shown as colored bars in the left panel of Figure 8.
• 2.  
  For γ values in our chain within a given bin, we then compute the median of all other parameters. Because of the way we chose our γ bins, this yields the quantile of interest for γ, whereas the other parameters will track their corresponding degeneracy direction with respect to γ. This can be seen in the colored squares in the right panel of Figure 8.
• 3.  
  For the set of parameters at each of the quantiles (e.g., the 84% quantile in γ and the median in all other parameters for the corresponding bin) we can then generate a model using our GP emulator.

The result of this procedure is shown in Figure 9 for the power spectrum at redshift z = 2.8, which is the highest redshift showing a high γ value. We compare models generated in this way to the measured power spectra shown as the blue and green points in the figure. Bands show the 68% confidence interval at each k for models generated using our emulator with random draws from the posterior distribution. Note that the forward model (due to both masking and forward modeling of noise and resolution) can generate slightly more converged model power spectra than the perfect model using the same parameters. The latter band is therefore actually a prediction for k ≳ 0.02 $\,{\rm{s}}\ {\mathrm{km}}^{-1}$ and its slightly larger extent is not surprising. Also note that due to the way we chose to produce curves with different thermal parameters and the dimensionality of the space the range spanned by the dashed curves is typically smaller than the colored bands. This is expected as the band shows the actual spread in the five/six (depending on the number of data sets used) dimensional parameter space, whereas the lines are based on a quantile for one of the parameters and values at the center of the distribution close to that quantile for all others, which will lead to a point inside the respective hypersurface, e.g., parameters of the purple/blue curve fall inside the five/six dimensional 68% surface, where the band corresponds to the actual surface).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We can see that all five models shown basically lead to the same power except for the highest k-values measured k ≥ 0.07 (smallest scales). At those scales a higher γ and lower $\bar{F}$ indeed seem to provide a better fit to the data, whereas at larger scales (smaller k) the model does not seem to be strongly affected by the parameters when moving along the degeneracy.

However, for other redshift bins (see Figure 10) the sensitivity of the power spectrum toward changes in γ for a region around the median value shifts to different scales. For example, at z ≤ 2.0 the most dominant effect seems to be on large scales, but note that we do not have the high-precision BOSS measurement and therefore both the range in allowed power spectra and the range of parameters in the 2σ region of γ are larger. All other redshifts seem to suggest a highest sensitivity to γ at scales k ∼ 0.05 $\,{\rm{s}}\ {\mathrm{km}}^{-1}$, different from both the lowest redshifts and z = 2.8. While we note that differences between models of different γ along the degeneracy direction are typically small compared to our measurement errors for an individual k-bin, it is clear that the data of all bins combined has the precision to distinguish between these models, and that our inference is producing sensible fits. One might argue that the fact that the k-modes that are driving the fits to high γ and low $\bar{F}$ change for different redshift bins is a source of concern, but we caution that the degeneracies in this multi-dimensional parameter space are complex and not always easy to visualize. We are confident that these results are not spurious, since this high γ, low $\bar{F}$ combination persists consistently across all redshift bins with z ≤ 2.8, and both measurements and our inference of different redshift bins are completely independent. We will return to this issue of discrepant γ and $\bar{F}$ values in Section 6 when we discuss possible systematic errors in our hydrodynamical simulations.

Zoom In Zoom Out Reset image size

Given that independent precise constraints on the mean transmission exist we now consider the effect of applying a Gaussian prior on the mean transmission based on these measurements (see discussion in Section 4.4 for details). Henceforth we will refer to these fits as the "strong prior" results, and we will designate them as our fiducial measurements (as opposed to the joint fits for thermal parameters and $\bar{F}$ described in previous sections). Note that most previous analyses of the IGM thermal properties have simply assumed perfect knowledge of the mean transmission (see Lidz et al. 2010; Iršič et al. 2017b for exceptions), such that this "strong prior" approach is more consistent with previous efforts.

We present the redshift evolution of posterior parameter degeneracies assuming the strong prior in Figure 11. Each panel in these figures shows the 2d-marginalized 68% and 95% confidence regions of T₀ versus γ. While γ and T₀ are strongly anticorrelated at low redshifts z ≤ 3.4, i.e., the contours are close to diagonal, this correlation gets weaker at higher redshifts (especially at z ≥ 4.2), i.e., contours become aligned with the axes due to lower overdensities probed by the power spectrum. Likewise, the γ versus $\bar{F}$ confidence regions are shown in Figure 12. Note that these properties are correlated independent of redshift, in stark contrast to the thermal parameter degeneracy, while still changing shape and direction due to the different precision of the measurements. Therefore, a change of prior for the mean transmission measurements propagates into γ at high redshifts (z ≥ 4.2), but does not affect T₀ significantly. At lower redshifts (especially for z ≤ 3.4), however, γ is strongly correlated with both T₀ and $\bar{F}$, so a change in priors for any of the three quantities always affects the results on the other two quantities as well. Consequently, the change in our mean flux prior affects lower redshifts (especially z < 3) more strongly than higher ones.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We show the fully marginalized posterior constraints on thermal parameters as a function of redshift in Figure 13 (for tabulated values, see Table 4 with further details in the Appendix). We can see that now the values of γ cover the theoretically expected value of γ ≈ 1.6 at low redshifts z < 3, while the values of T₀ obtained are higher than in the fit using a flat prior on mean transmission because of the degeneracies between T₀ and γ. This is more clearly illustrated in Figure 14 where we compare T₀ and γ evolution for the different prior assumptions. We can see that indeed the changes between both the two fits are strongly anticorrelated between T₀ and γ and that the change in marginalized parameters between the two cases can be large, particularly for γ where the differences at 2.2 ≤ z ≤ 2.6 are ≳2σ and as high as 3σ at z = 2.8.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

However, this is the first thermal evolution measurement performed over the whole epoch of $\mathrm{He}\,{\rm{II}}$ reionization and beyond based on the power spectrum. Using the strong mean flux prior we also obtained reasonable results, including physically possible measurements of γ, a rise in temperature for z ≳ 3 as time progresses (or redshift decreases), and the first measurement of the IGM cooling down thereafter. In the next sections we compare our strong prior results to recent thermal parameter measurements from different methods as well as models of IGM thermal evolution.

A comparison of our results to recent measurements of thermal parameters is shown in Figure 15. We discuss the various data sets involved and elaborate on the comparison to our new measurement below.

Zoom In Zoom Out Reset image size

The phase angle PDF of quasar pairs (Rorai et al. 2013) measures the smoothness of the 3d distribution of IGM gas and therefore directly constrains the pressure smoothing scale λ_P independent of the instantaneous thermal state of the IGM (i.e., T₀ and γ). Rorai et al. (2017b) measured λ_P from a sample of quasar pairs in four redshift bins between $2.0\leqslant z\leqslant 3.6$. Figure 15 shows that our inferred values of λ_P are fully consistent with the Rorai et al. (2017b) measurement. We also see a smaller uncertainty in our power-spectrum-based measurement. Part of the explanation for these small error bars lies in how our model grid, and therefore our prior probability, is set up. As discussed in Section 4.4, the degeneracy between T₀ and λ_P within our simulated models combined with our approach of not extrapolating to regions outside our model grid result in a positive correlation in the prior probability between these parameters. However, the power spectrum cutoff is sensitive to a degenerate combination of both λ_P and thermal broadening (Peeples et al. 2010; Rorai et al. 2013), leading to an anticorrelation in the likelihood. So the correlations inside our prior (Figure 4) and the degeneracy direction of the likelihood due to the aforementioned effect are nearly perpendicular, and as the posterior distribution is the product of these two, the resulting constraints appear very tight. However, we argue that it is hard to generate physical models without imprinting the correlation between thermal state and λ_P that depends on the integrated thermal history of the IGM. While the uncertainties in λ_P might still be somewhat underestimated, we note that our prior grid degeneracy has a strong physical motivation (see also Section 4.4).

The orange points in Figure 15 show the T₀ and γ measurements from Lidz et al. (2010), who decomposed the Lyα forest of the Dall'Aglio et al. (2008) data set into wavelets and analyze the PDF of their squared amplitudes to derive constraints on the thermal state of the IGM. We note that this data is a subset of that used to compute the power spectrum in Paper I and analyzed here. Note that their γ constraint is often limited to the boundaries of their fits.¹⁹ While the wavelet analysis results at z = 2.6 are consistent with our measurement, we disfavor the z = 2.2, z = 3.0, z = 4 and especially z = 3.4 wavelet results, which seem to indicate a far hotter IGM than our measurement. The origin of this discrepancy is unclear, but it has also been noted before by Becker et al. (2011).

Another method for obtaining constraints on the thermal state of the IGM is by decomposing the Lyα forest into individual absorption lines, assuming that a cutoff in the distribution of column densities ${N}_{{\rm{H}}{\rm{I}}}$ versus Doppler parameter b exists and can be attributed to lines that are only thermally broadened (see, e.g., Schaye et al. 2000; Rudie et al. 2012; Bolton et al. 2014; Hiss et al. 2018; Rorai et al. 2018). Especially the new Hiss et al. (2018; which is based on the same data set as Paper I and is using a subset of the same simulation grid) thermal evolution result seems to hint toward a period of heating until z ∼ 2.8 that could be attributed to $\mathrm{He}\,{\rm{II}}$ reionization.

For both T₀ and γ we see broad agreement between our measurements and the line-fitting results at most redshifts. Of particular interest are z = 2.4 and z = 2.8, where several line-fitting measurements exist. At z = 2.4 we do reproduce the result from Hiss et al. (2018; blue points) as well as Bolton et al. (2014; green) in both T₀ and γ. At z = 2.8 we agree with the Rorai et al. (2018; brown point), but obtain higher precision. However, agreement with Hiss et al. (2018) at this redshift seems to be poor as they measure both higher T₀ and lower γ (which is along the degeneracy direction for line-fitting analyses as well as the power spectrum). Part of this discrepancy might come from systematics in the Voigt profile analysis, depending on the cutoff fitting algorithm chosen, as Hiss et al. (2018; see Appendix B) find either a multimodal posterior probability distribution for T₀, γ with a similar 68% confidence interval as Rorai et al. (2018), or a unimodal distribution with the values shown here depending on the cutoff fitting algorithm used. Whether this multimodal behavior results from systematics in the measurement procedure or is a real physical effect from, e.g., a real multimodal IGM TDR is not yet clear, but we do not see such behavior in our power spectrum analysis. For the other overlapping redshifts (except z = 2.2 and z = 2.4, which match very well) we generally measure a lower T₀ and higher γ compared to Hiss et al. (2018).

The most precise measurements of temperature in the IGM so far are based on the mean curvature in the Lyα forest $\langle \kappa \rangle$ (Becker et al. 2011; Boera et al. 2014). These measurements constrain T(Δ_⋆) at an optimal overdensity Δ_⋆ at which a one-to-one relation between the mean curvature $\langle \kappa \rangle$ of the Lyα forest and T(Δ_⋆) exists independent of the slope γ of the TDR (note again Figure 11, which shows the corresponding degeneracy for the power spectrum). However this method is not able to measure γ or T₀ independently. To compare to curvature-based measurements, we compute ${T}_{{{\rm{\Delta }}}_{\star }}={T}_{0}{{\rm{\Delta }}}_{\star }^{\gamma -1}$ (using the values for ${{\rm{\Delta }}}_{\star }$ given by Becker et al. 2011) for each sample in our MCMC chain and evaluate the 68% confidence interval. This approach allows us to directly compare to what the curvature results measure.

The agreement with the curvature analysis seems to be generally good for the largest part of the overlapping redshift range, but we seem to measure overall slightly higher temperatures. There are some redshifts $z=2.6,2.8,3.2,3.4$ where our analysis gives significantly higher temperatures than implied by the curvature measurements. Note in particular that at z = 2.8, where we see the strongest discrepancy between our results and the curvature measurements, multiple measurements of the thermal state have been performed via several different methods, and these results do not fully agree with each other. We argue that the overall agreement is still good given the significantly different data sets, statistical approaches, and models used for both types of analysis. For example, the difference in measured thermal state might potentially arise due to the different sensitivity of both statistics to metal contamination. While the power spectrum is only weakly affected by residual metal lines on the very smallest scales we cover here (see the comparison in Walther et al. 2018), the averaged squared curvature is basically measuring ${\int }_{{k}_{m}{in}}^{\infty }{k}^{5}P(k)d\,\mathrm{ln}\,k$ (see Appendix D in Puchwein et al. 2015) and thus enhances the weight of residual small-scale contamination in the Lyα forest. At the same time small-scale contaminants, like, e.g., leftover metal lines, would decrease the obtained IGM temperatures as there is now too much small-scale power, thus leading to a colder IGM in curvature than in power spectrum analyses. Additionally, these measurements did not marginalize over the mean flux in the simulations, thereby, e.g., potentially underestimate their errors.

Finally, we also show the Garzilli et al. (2017) measurement of T₀ at $4.2\leqslant z\leqslant 5.4$ based on the same as the Viel et al. (2013) data set we use here (gray points, the limit is at the 1σ level), but using a different analysis pipeline and including a WDM particle mass as an additional free parameter. We can see that for $z\leqslant 5$ the agreement is good, but for z = 5.4 we seem to get slightly higher values of T₀ than their 1σ upper limit. Part of that difference can be attributed to the additional freedom in their model.

Overall, we conclude that the agreement between our data and previous results is reasonably good. Our measurement is a strong advancement with regard to previous analyses, especially due to the large range of uniformly covered redshifts and due to jointly constraining T₀ and γ over this full range.

In the previous sections we performed a self-consistent measurement of thermal evolution in the IGM from z = 5.4 to z = 1.8, corresponding to 3 Gyr of cosmic history. In this section, we compare to simulations to thermal evolution due to $\mathrm{He}\,{\rm{II}}$ reionization, as this is expected to be the dominant process setting the thermal state of the IGM at this epoch.

In Figure 16, we show comparisons between our thermal evolution measurement and models based on different approaches. The solid curves show the "explicit reionization" simulations from our model grid for which hydrogen reionizes (to a level ${x}_{{\rm{H}}{\rm{II}}}=99.9 \%$, note that for our models this point is typically reached with a delay of ${\rm{\Delta }}z\approx 1$ compared to the corresponding ${z}_{\mathrm{reion},50}$ ²⁰ ) at ${z}_{\mathrm{reion},{\rm{H}}{\rm{I}}}=6.5$ in agreement with the Planck Collaboration et al. (2016a; and also Planck Collaboration et al. 2018) constraints, but for which the parameters governing $\mathrm{He}\,{\rm{II}}$ reionization are varied (see Table 2). The red dashed-dotted curve is showing an extreme version of these models for which $\mathrm{He}\,{\rm{II}}$ was never reionized.²¹ The measured temperature at mean density is significantly higher (by ≳3σ for each of the seven individual redshift bins with 2.2 ≤ z ≤ 3.4) than this "no $\mathrm{He}\,{\rm{II}}$" scenario for z ≤ 4.6 suggestive of a period of $\mathrm{He}\,{\rm{II}}$ reionization taking place.

Zoom In Zoom Out Reset image size

Table 2 .  Thermal Evolution Models Used in Comparisons to Existing Measurements: Parameters are the Reionization Redshifts and the Total Heat Input during Reionization for ${\rm{H}}\,{\rm{I}}$ and $\mathrm{He}\,{\rm{II}}$, see Oñorbe et al. (2017a) for Details

Model Name	${z}_{\mathrm{reion},{\rm{H}}{\rm{I}}}$	${z}_{\mathrm{reion},\mathrm{He}{\rm{II}}}$	${\rm{\Delta }}{T}_{{\rm{H}}{\rm{I}}}[{\rm{K}}]$	${\rm{\Delta }}{T}_{\mathrm{He}{\rm{II}}}[{\rm{K}}]$
No $\mathrm{He}\,{\rm{II}}$	7.3	⋯	20000	⋯
Cold $\mathrm{He}\,{\rm{II}}$	6.55	3.0	20000	10000
Standard $\mathrm{He}\,{\rm{II}}$	6.55	3.0	20000	15000
Warm $\mathrm{He}\,{\rm{II}}$	6.55	3.0	20000	20000
Hot $\mathrm{He}\,{\rm{II}}$	6.55	3.0	20000	30000
Late $\mathrm{He}\,{\rm{II}}$	6.55	2.8	20000	15000

Download table as:  ASCIITypeset image

Table 3.  Fiducial Evolution of Thermal Parameters Assuming a Flat Prior on $\bar{F}$

z	λP	T0	γ	$\bar{F}$	${T}_{{{\rm{\Delta }}}_{\mathrm{power}}}$	${T}_{{{\rm{\Delta }}}_{\star }}$	${{\rm{\Delta }}}_{\mathrm{power}}$
	(kpc)	(${10}^{4}\,\,{\rm{K}}$)			(${10}^{4}\,\,{\rm{K}}$)	(${10}^{4}\,\,{\rm{K}}$)
1.8	${79.0}_{-11.9}^{+16.0}$	${0.684}_{-0.122}^{+0.180}$	${1.97}_{-0.26}^{+0.16}$	${0.872}_{-0.018}^{+0.020}$	${1.160}_{-0.239}^{+0.239}$	${4.288}_{-1.525}^{+2.162}$	1.705
2.0	${93.0}_{-17.4}^{+8.3}$	${0.734}_{-0.071}^{+0.093}$	${2.15}_{-0.26}^{+0.09}$	${0.831}_{-0.011}^{+0.033}$	${1.096}_{-0.121}^{+0.122}$	${5.749}_{-2.290}^{+1.011}$	1.437
2.2	${91.0}_{-6.4}^{+6.3}$	${0.789}_{-0.068}^{+0.085}$	${2.13}_{-0.13}^{+0.09}$	${0.796}_{-0.009}^{+0.010}$	${1.369}_{-0.106}^{+0.120}$	${4.942}_{-0.771}^{+0.773}$	1.638
2.4	${87.2}_{-5.1}^{+5.4}$	${0.831}_{-0.078}^{+0.112}$	${2.07}_{-0.18}^{+0.13}$	${0.772}_{-0.012}^{+0.013}$	${1.593}_{-0.123}^{+0.143}$	${3.995}_{-0.631}^{+0.717}$	1.841
2.6	${88.3}_{-4.5}^{+3.7}$	${1.000}_{-0.090}^{+0.146}$	${1.93}_{-0.17}^{+0.15}$	${0.745}_{-0.013}^{+0.012}$	${1.936}_{-0.084}^{+0.095}$	${3.449}_{-0.345}^{+0.445}$	2.012
2.8	${93.8}_{-4.2}^{+4.2}$	${1.000}_{-0.087}^{+0.112}$	${2.16}_{-0.13}^{+0.09}$	${0.688}_{-0.010}^{+0.013}$	${1.982}_{-0.149}^{+0.163}$	${3.911}_{-0.408}^{+0.426}$	1.818
3.0	${80.6}_{-5.6}^{+6.0}$	${1.429}_{-0.271}^{+0.313}$	${1.47}_{-0.24}^{+0.26}$	${0.694}_{-0.015}^{+0.009}$	${2.027}_{-0.143}^{+0.155}$	${2.347}_{-0.226}^{+0.269}$	2.110
3.2	${84.9}_{-6.3}^{+4.7}$	${1.115}_{-0.149}^{+0.230}$	${1.85}_{-0.25}^{+0.21}$	${0.623}_{-0.019}^{+0.018}$	${1.910}_{-0.150}^{+0.167}$	${2.465}_{-0.253}^{+0.278}$	1.882
3.4	${90.1}_{-5.8}^{+4.7}$	${1.330}_{-0.215}^{+0.295}$	${1.82}_{-0.27}^{+0.24}$	${0.569}_{-0.023}^{+0.021}$	${2.202}_{-0.214}^{+0.206}$	${2.592}_{-0.277}^{+0.283}$	1.791
3.6	${79.4}_{-9.6}^{+10.0}$	${1.010}_{-0.296}^{+0.360}$	${1.74}_{-0.36}^{+0.28}$	${0.512}_{-0.021}^{+0.022}$	${1.160}_{-0.338}^{+0.394}$	${1.704}_{-0.561}^{+0.602}$	1.194
3.8	${79.4}_{-6.7}^{+8.4}$	${1.029}_{-0.246}^{+0.287}$	${1.74}_{-0.39}^{+0.29}$	${0.433}_{-0.026}^{+0.025}$	${1.320}_{-0.273}^{+0.355}$	${1.548}_{-0.338}^{+0.441}$	1.442
4.0	${72.3}_{-5.8}^{+7.6}$	${0.863}_{-0.187}^{+0.271}$	${1.42}_{-0.34}^{+0.37}$	${0.387}_{-0.022}^{+0.017}$	${0.942}_{-0.201}^{+0.288}$	${1.090}_{-0.262}^{+0.340}$	1.218
4.2	${77.0}_{-6.0}^{+3.6}$	${0.905}_{-0.082}^{+0.122}$	${1.73}_{-0.40}^{+0.33}$	${0.355}_{-0.031}^{+0.025}$	${1.051}_{-0.082}^{+0.087}$	${1.246}_{-0.165}^{+0.118}$	1.215
4.6	${73.7}_{-5.8}^{+4.9}$	${0.910}_{-0.117}^{+0.119}$	${1.54}_{-0.39}^{+0.37}$	${0.278}_{-0.028}^{+0.023}$	${0.966}_{-0.109}^{+0.126}$	${1.037}_{-0.124}^{+0.153}$	1.134
5.0	${57.3}_{-4.3}^{+4.0}$	${0.535}_{-0.092}^{+0.117}$	${1.54}_{-0.33}^{+0.31}$	${0.159}_{-0.020}^{+0.018}$	${0.555}_{-0.095}^{+0.119}$	${0.580}_{-0.102}^{+0.122}$	1.067
5.4	${54.4}_{-4.5}^{+4.3}$	${0.597}_{-0.132}^{+0.152}$	${1.55}_{-0.29}^{+0.29}$	${0.060}_{-0.008}^{+0.009}$	${0.551}_{-0.118}^{+0.138}$	${0.613}_{-0.139}^{+0.158}$	0.868

Download table as:  ASCIITypeset image

To allow a comparison between different $\mathrm{He}\,{\rm{II}}$ reionization scenarios, the gray, blue, green, and purple curves show models with ${z}_{\mathrm{reion},\mathrm{He}{\rm{II}}}=3.0$, assuming different amounts of heat being injected ${\rm{\Delta }}{T}_{\mathrm{He}{\rm{II}}}$ into the IGM, varying from 10,000 K (cold) to 30,000 K (hot); whereas the orange curve shows a model with ${\rm{\Delta }}{T}_{\mathrm{He}{\rm{II}}}={\rm{15,000}}\,\,{\rm{K}}$ but ${z}_{\mathrm{reion},\mathrm{He}{\rm{II}}}=2.8$ (late). We can see that the models predict an extended period of heating (i.e., increasing T₀) until ${z}_{\mathrm{reion},\mathrm{He}{\rm{II}}}$ followed by the IGM cooling down due to the expansion of the universe whose effects on the thermal state cannot be fully counteracted by ionizations anymore. Overall, for our measurement this rise and fall in T₀ lies between the standard and warm $\mathrm{He}\,{\rm{II}}$ evolution models for 2.2 ≤ z ≤ 4.6, disfavoring particularly hot or late phases of $\mathrm{He}\,{\rm{II}}$ reionization.

We also compare to the analytical thermal evolution model by Upton Sanderbeck et al. (2016) and the fiducial non-equilibrium reionization model by Puchwein et al. (2018; see their Figure 6). We note that while the general shape of thermal evolution looks similar to both models for z ≤ 4.6, we seem to obtain a slightly less pronounced peak in T₀. Overall, the temperature evolution we see in this redshift range is indeed well modeled by an $\mathrm{He}\,{\rm{II}}$ reionization event followed by photoionization equilibrium in an adiabatically expanding IGM. While it has been argued that this effect has been seen before (Becker et al. 2011), previous work did not break the degeneracy between γ and T₀. Note that also the cooldown of the IGM after reionization has never been conclusively observed due to this degeneracy.

Models of $\mathrm{He}\,{\rm{II}}$ reionization typically also show a dip in γ resulting from the IGM being more isothermal during reionization events (see also, e.g., McQuinn et al. 2009). We can also see this effect by comparing γ for our $\mathrm{He}\,{\rm{II}}$ models with the no-$\mathrm{He}\,{\rm{II}}$ model. Note that while the Upton Sanderbeck et al. (2016) model also shows a dip (albeit at later times and with a more strongly isothermal γ), the fully non-equilibrium simulation by Puchwein et al. (2018) does show an intrinsically smaller γ and no strong dip. The reason for this is that the non-equilibrium model reached γ = 1 at z = 7 due to ${\rm{H}}\,{\rm{I}}$ reionization and is still recovering from this feature, i.e., it did not yet forget about the timing of ${\rm{H}}\,{\rm{I}}$ reionization. The "dip" for this model therefore manifests in the near constant evolution from z ∼ 5 to z ∼ 3 compared to an otherwise expected rise in γ.

We can see this dip in γ for the measurement at z ∼ 3.9 aligned in redshift with the expected decrease due to $\mathrm{He}\,{\rm{II}}$ reionization in our explicit $\mathrm{He}\,{\rm{II}}$ reionization models. Note that the dip is only ∼2σ significant compared to the no-$\mathrm{He}\,{\rm{II}}$ region model, but overall a slightly higher value for γ than this model is preferred. Also note that on the data side this feature is currently dominated by XQ-100 data (which is the highest resolution data available at 3.6 ≤ z ≤ 4.0), which we strongly degraded by marginalizing over resolution. Additional high-resolution data or an accurate determination of the XQ-100 data set resolution at these redshifts and adopting a prior based on those results could therefore lead to additional constraints on $\mathrm{He}\,{\rm{II}}$ reionization due to its signature in the slope of the TDR.

Note that this feature also strongly relies on precise knowledge of $\bar{F}$, as the expected decrease is very shallow. Additionally, γ values for z > 4.2 might have a significant uncertainty as measurements of the mean transmission get less accurate for this range due to the smaller amounts of data available and stronger fluctuations in the ionization state of the IGM. Thus, there are currently several discrepant measurements of $\bar{F}$ (as discussed in Section 4.4), which consequently lead to a high uncertainty in γ.

At early times (z ≥ 5, we call those points the highest redshift measurements) we can see that the measured T₀ is lower than in any of the models. Note again that similarly low temperatures were also obtained by Garzilli et al. (2017) based on the same data set in a fully independent analysis. While one could in principle think that an earlier redshift of ${\rm{H}}\,{\rm{I}}$ reionization gives the IGM more time to cool thereafter leading to lower temperatures at these times, models suggest that is not the case, and T₀ has essentially forgotten about the timing of reionization by z = 5.4 (see, e.g., Oñorbe et al. 2017b who present models for a range of different $6.0\lt {z}_{\mathrm{reion},{\rm{H}}{\rm{I}}}\lt 9.7$). Instead, the post-reionization thermal state mostly depends on the spectral shape of the UVB (McQuinn & Upton Sanderbeck 2016) and a low temperature at z ∼ 5.4 requires lower photoheating rates, i.e., a particularly soft spectral energy distribution for ionizing sources is needed, which is not favored by current models of the UVB (Faucher-Giguère et al. 2009; Haardt & Madau 2012; Stanway et al. 2016; D'Aloisio et al. 2018b; Puchwein et al. 2018).

While it may be that the thermal state at z ≃ 5.4 would still be sensitive to the reionization redshift for the particularly late z ≲ 6 reionization scenario (which now seems to be allowed regarding the newest cosmic microwave background (CMB) results from Planck Collaboration et al. 2018), this would nevertheless need to be in conjunction with very low reionization heat injection. The recent results by D'Aloisio et al. (2018b) who used radiative transfer to simulate photoheating by ionization fronts during ${\rm{H}}\,{\rm{I}}$ reionization suggest that such low levels of IGM heating are unlikely. Finally, note that the onset of $\mathrm{He}\,{\rm{II}}$ reionization can only increase model temperatures and therefore worsen the disagreement, as none of the models shown exhibits any $\mathrm{He}\,{\rm{II}}$ reionization before z = 4.8.

Therefore, the small temperatures we (and also other groups using the same data set) obtain at z ≥ 5 are challenging to fit with current models of reionization. Consequently, models fitting the low-T measurements would also lead to a colder IGM at later times without additionally increasing heating due to, e.g., $\mathrm{He}\,{\rm{II}}$ reionization. However, as current constraints at the highest redshifts rest upon the single data set by Viel et al. (2013) based on a handful of objects, future measurements based on larger samples of quasar spectra obtained might change those low-T₀ results.