Joint Bayesian Estimation of Quasar Continua and the Lyα Forest Flux Probability Distribution Function

The fluctuating Gunn–Peterson approximation (FGPA; see e.g., Croft et al. 1998, 1999; Weinberg et al. 1998) provides a relationship between the observed flux in the Lyα forest and the underlying matter density distribution. It represents a valuable tool for generating Lyα absorption spectra from cosmic matter density fields, which we will use to generate mock quasar spectra. Assuming that the intergalactic gas is in photoionization equilibrium and follows the temperature–density relation in Equation (1), the optical depth τ of the IGM scales with the matter density contrast as

$$\begin{eqnarray}&&\tau \propto {{\rm{\Delta }}}^{2-0.7(\gamma -1)}.\end{eqnarray} \tag{ 2 }$$

The normalized transmitted flux F in the Lyα forest is related to the intervening optical depth of the photoionized gas by

$$\begin{eqnarray}&&F=\exp (-\tau ).\end{eqnarray} \tag{ 3 }$$

Combining Equations (2) and (3) leads to a relationship between the transmitted flux F and the matter density contrast Δ.

McDonald et al. (2006) provide a semi-analytic model to generate absorption fields for the Lyα forest from matter density fields based on the log-normal model introduced by Bi et al. (1992) and Bi & Davidsen (1997). The goal of this semi-analytic model is to generate realistic flux fields F that approximately reproduce both the observed power spectrum as a function of wave number k and redshift z and the observed PDF of the Lyα forest. We introduce small modifications to this model to accommodate different values for γ in the temperature–density relation, and summarize its main features below.

The generation of each F field for the ith quasar spectrum begins by creating an initial Gaussian random field ${\delta }_{i,0}$ representing the underlying dark matter density contrast with a power spectrum given by:

$$\begin{eqnarray}&&{P}_{\delta }(k)=\displaystyle \frac{1+{\left(\tfrac{0.01{\rm{s}}{\mathrm{km}}^{-1}}{{k}_{0}}\right)}^{\nu }}{1+{\left(\tfrac{k}{{k}_{0}}\right)}^{\nu }}\exp [-{({{kR}}_{\delta })}^{2}],\end{eqnarray} \tag{ 4 }$$

with values of ${k}_{0}=0.001\,{\rm{s}}$ km⁻¹, $\nu =0.7$, ${R}_{\delta }=5.0$ km s⁻¹, which were chosen to approximately reproduce the dependence on wave number k in the observed flux power spectrum of the Lyα forest. An evolution of the amplitude A(z) of the density fluctuations with redshift is introduced by the transformation

$$\begin{eqnarray}&&{\delta }_{i}=A({z}_{i}){\delta }_{i,0}\end{eqnarray} \tag{ 5 }$$

with

$$\begin{eqnarray}&&{A}^{2}({z}_{i})=58.6{\left(\displaystyle \frac{1+{z}_{i}}{4}\right)}^{-2.82}.\end{eqnarray} \tag{ 6 }$$

The values for the amplitude and exponent here were chosen such that the final flux power spectrum would evolve like the observed one. In order to obtain the density field n_i of the baryonic diffuse matter of the IGM we do not use the squared log-normal transformation that has been used in previous work (see e.g., Bi et al. 1992; Bi & Davidsen 1997; McDonald et al. 2006), but instead use a log-normal transformation with an exponent of $2-\beta$ with $\beta =0.7(\gamma -1)$. This step introduces the dependence on the slope parameter γ of the temperature–density relation of the IGM that we are interested in. Hence the transformation of the density field n_i is computed as

$$\begin{eqnarray}&&{n}_{i}={n}_{0}{\left[\exp \left({\delta }_{i}-\left\langle \displaystyle \frac{{\sigma }_{i}^{2}}{2}\right\rangle \right)\right]}^{2-0.7(\gamma -1)},\end{eqnarray} \tag{ 7 }$$

where n₀ is the mean number density of the IGM. The factor ${\sigma }_{i}^{2}$ is computed from the input density power spectrum in Equation (4) and fixes the mean of the lognormal field to unity.

After the transformation the n_i field is smoothed with a Gaussian filter with standard deviation ${\sigma }_{\tau }=20$ km s⁻¹. Note that the quasar spectra that we will consider in this work have a resolution comparable to the achievable resolution of echelle spectrographs (see Section 2.3 for details), which is higher than the Gaussian filter applied here. However, since this smoothing is applied to the n_i field and not to the flux F itself, it accounts for the effects that give rise to the velocity width of the Lyα forest lines, such as thermal broadening, Jeans pressure smoothing (see Gnedin & Hui 1998; Rorai et al. 2013; Kulkarni et al. 2015; Oñorbe et al. 2017), or peculiar velocities and velocities due to the Hubble flow. We then multiply by a redshift evolution factor of $0.374{[(1+{z}_{i})/4]}^{5.1}$ (the form of this factor is again chosen to reproduce the observations) to produce a field τ. The transmitted flux can then be calculated using Equation (3).

Note that this simple toy model of the IGM only accommodates a dependence on the thermal parameter γ and ignores dependences on the temperature T₀ or the pressure smoothing scale ${\lambda }_{P}$. However, Bolton et al. (2008) used hydrodynamical simulations to show that the Lyα flux PDF is mostly sensitive to γ, whereas the other thermal parameters introduce only relatively modest changes to its shape. If we applied our method to real quasar spectra, the simple model adopted here would clearly be insufficient, but rather we would use hydrodynamical simulations that accommodate dependences on all thermal parameters. However, the main advantage of this model is that it is fast and a large grid of models can be generated quickly. We will comment more on this later and argue that future augmentation of the parameter space of the PDF should be easy to incorporate.

We re-scale the optical depth in order to obtain values for the mean flux of our mock Lyα forest spectra that are consistent with the measurements of Faucher-Giguère et al. (2008) at redshift z = 3, i.e., $\langle F\rangle \approx 0.680$. For an estimate of the mean flux at z = 5 we use the fitting formula presented in Oñorbe et al. (2017), because the measurements of the mean flux do not extend past $z\gt 4.85$, and obtain $\langle F\rangle \approx 0.189$ at redshift z = 5.

Figure 1 illustrates the effect of changing the slope parameter γ on the flux PDF at redshifts of z = 3 and z = 5 in the left and right panels, respectively. Decreasing γ increases the temperature in the underdense regions of the IGM, i.e., ${\rm{\Delta }}\lt 1.0$, according to the temperature–density relation (see Equation (1)). Higher temperatures cause the hydrogen recombination rate to decrease, therefore reducing the fraction of neutral H i and the Lyα optical depth, which in turn increases the transmitted flux in these regions. These differences are apparent at $F\gt 0.5$, where the peak of the flux PDF shifts toward higher values for lower values of γ. Accordingly, by measuring the Lyα flux PDF one can constrain the parameter γ.

Zoom In Zoom Out Reset image size

In order to generate mock quasar spectra we multiply Lyα absorption fields with quasar continua. For this purpose we use the estimated continua of 50 real quasar spectra taken by the Hubble Space Telescope (HST) Faint Object Spectrograph from Suzuki et al. (2005), which were collected and calibrated by Bechtold et al. (2002). All quasars in this sample have complete wavelength coverage from 1020 Å to 1600 Å in the rest frame. These moderate-resolution spectra ($R\sim 1300$) have been combined from several exposures, brought into the rest frame and re-binned into pixels of size 0.5 Å. The average signal-to-noise ratio of the chosen quasar sample is $\langle {\rm{S}}/{\rm{N}}\rangle =19.5$ per 0.5 Å pixel, and quasars with ${\rm{S}}/{\rm{N}}\lt 10$ per pixel were removed from the sample. All quasar spectra are renormalized to unity at 1280 Å in the rest-frame. The quasars in the sample are at very low redshifts in the range $0.14\lt z\lt 1.04$ with a mean redshift of $\langle z\rangle =0.58$. Quasars at these low redshifts show very few Lyα absorption lines, which makes it easier to correctly estimate their continuum level. Spectra with broad absorption lines or damped Lyα systems were excluded from the sample, because of the associated large uncertainties in placing their continuum level. In order to obtain precise and smooth continuum estimations for these spectra, Suzuki et al. (2005) fitted Chebyshev polynomials of different orders to each HST spectrum and applied further fine-tuning adjustments afterwards using a B-spline fit. For more details on this procedure we refer the reader to Suzuki et al. (2005). We use these smoothed continua as the true continua of our mock quasar spectra.

We multiply the empirical quasar continua with different realizations of the Lyα forest absorption field to produce mock quasar spectra. To simplify the analysis we assume no redshift evolution of the quasar continua. The absorption fields are also generated at a fixed redshift, and there is thus no redshift evolution in the Lyα forest along the quasar sightline. These assumptions should not alter the conclusions of this work. Thus we generate our own mock spectra by multiplying the Lyα forest F at different redshifts into the quasar continua in the wavelength region where Lyα forest is usually found, i.e., between Lyβ (1025.18 Å) and Lyα (1215.67 Å) emission.

We choose the pixel size and resolution of our mock spectra to be comparable to real high-resolution spectra taken with the Ultraviolet and Visual Echelle Spectrograph (UVES) at the Very Large Telescope and the High Resolution Echelle Spectrometer (HIRES) from the Keck telescope. Since the pixel size from these spectrographs is smaller than the pixel size of 0.5 Å used by Suzuki et al. (2005), we interpolate the HST continua onto a finer grid with the desired pixel scale, which we choose to be 2.5 km s⁻¹. We then add random white noise to our mock data, i.e., we ignore correlations between the noise and the transmission, with an S/N that is achievable with UVES or HIRES. The highest-quality quasar spectra at $z\sim 3$ from these instruments achieve an S/N in the Lyα forest of ${\rm{S}}/{\rm{N}}\approx 80$ per 6 km s⁻¹ resolution element (see e.g., Dall'Aglio et al. 2008; Lehner et al. 2014; O'Meara et al. 2015). For our purposes we assume a slightly larger resolution element of 10 km s⁻¹, which gives ${\rm{S}}/{\rm{N}}\approx 103$ per resolution element, assuming the same noise properties. A pixel size of 2.5 km s⁻¹ implies 4 pixels per resolution element, and thus we obtain ${\rm{S}}/{\rm{N}}\approx 51$ per pixel. For quasars at $z\sim 5$ we will assume a lower ${\rm{S}}/{\rm{N}}\approx 20$ per pixel, consistent with existing high-resolution quasar spectra at this redshift (Calverley et al. 2011). However, we will see in Section 3.2 that our results at z = 5 are fairly independent of the exact noise level, because the Lyα flux PDF at these redshifts is only mildly sensitive to noise. We add Gaussian-distributed white noise to our mock quasar spectra with a standard deviation ${\sigma }_{\mathrm{noise}}=\tfrac{\langle F\rangle }{{\rm{S}}/{\rm{N}}}$, where $\langle F\rangle$ is the mean transmitted flux and ${\rm{S}}/{\rm{N}}\approx 51$ or ${\rm{S}}/{\rm{N}}\approx 20$ for spectra at $z\sim 3$ and $z\sim 5$, respectively. Examples of our mock quasar spectra at redshift z = 3 and z = 5 are presented in Figure 2.

Zoom In Zoom Out Reset image size

The quasar continuum is modeled using a PCA, enabling us to describe the continuum shape via a simple linear model. Suzuki (2006) analyzed the shape of 50 HST quasar spectra and established a set of principal component spectra (PCS) that reproduce the continuum shape of each quasar with high accuracy. This accuracy in the continuum model is crucial in order to prevent biases when studying the statistical properties of the Lyα forest.

The idea of the PCA is to represent the continuum spectrum $| {q}_{i,\lambda }\rangle$ of the ith quasar by a reconstructed spectrum that consists of a mean spectrum $| {\mu }_{\lambda }\rangle$ and a sum of m weighted PCS $| {\xi }_{j,\lambda }\rangle$, where the index λ denotes the wavelength. Hence

$$\begin{eqnarray}&&| {q}_{i,\lambda }\rangle \approx | {\mu }_{\lambda }\rangle +\displaystyle \sum _{j=1}^{m}{\alpha }_{{ij}}| {\xi }_{j,\lambda }\rangle ,\end{eqnarray} \tag{ 8 }$$

where $| {\xi }_{j,\lambda }\rangle$ refers to the jth PCS and ${\alpha }_{{ij}}$ is its weight for quasar i.

Suzuki (2006) derived 10 PCS, which are shown in Figures 2 and 3 of their paper. We show their mean quasar spectrum and the first four PCS in Figure 3. The more components that are included in the reconstruction, the lower the variance of the individual spectra about the model fits. Table 1 in Suzuki (2006) summarizes the fraction of the variance accounted for depending on the number of principal components employed.

Zoom In Zoom Out Reset image size

In our analysis we will vary the number of PCS that we include when estimating the quasar continuum. Since we are only taking a finite number of components we cannot account for the total variance seen in the quasar spectra and thus we inevitably introduce an error in the continuum fit. In order to analyze the best precision possible for a chosen PCA basis, we estimate the continuum model coefficients by using an MCMC algorithm to sample the likelihood function ${ \mathcal L }=\exp (-{\chi }^{2}/2)$ with

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{\lambda }\displaystyle \frac{{({C}_{\mathrm{data},i\lambda }-{C}_{\mathrm{model},i\lambda }({\alpha }_{{ij}}))}^{2}}{{\sigma }_{i,\mathrm{noise}}^{2}}.\end{eqnarray} \tag{ 9 }$$

Here ${C}_{\mathrm{data},i\lambda }$ is the true quasar continuum determined by Suzuki et al. (2005) (see Section 2.2), ${C}_{\mathrm{model},i\lambda }$ is the model continuum from Equation (8), and ${\sigma }_{i,\mathrm{noise}}$ is the standard deviation of the white noise added to each spectrum. We use the means of the posterior probability distributions of each coefficient ${\alpha }_{{ij}}$ as our best continuum model for quasar i. Note that for this estimation we take all pixels in the whole available wavelength range, i.e., between 1020 Å and 1600 Å, into account and that we do not add artificial Lyα forest to the true quasar continua here. We then evaluate the precision of the estimated continuum model by calculating the continuum residuals, i.e., the relative continuum error,

$$\begin{eqnarray}&&{\rm{\Delta }}C/C=\displaystyle \frac{{C}_{\mathrm{data},\lambda }-{C}_{\mathrm{model},\lambda }}{{C}_{\mathrm{data},\lambda }},\end{eqnarray} \tag{ 10 }$$

only in the wavelength range that we are particularly interested in, i.e., the Lyα forest region between the Lyβ and Lyα emission peaks. However, in order to avoid biases due to the influence of proximity zones, and to guarantee that we exclude the Lyβ forest in the presence of quasar redshift errors, we only include the wavelengths between 1040 Å to 1190 Å. Note that all our mock spectra use the true quasar continua determined from the HST spectra.

The relative continuum error ${\rm{\Delta }}C/C$ is shown in the upper panel of Figure 4 for different numbers of PCS. The distribution for the ensemble of all 50 spectra follows a Gaussian distribution reasonably well. Note that this is not necessarily the case for a fit to a single quasar spectrum. The Gaussian distribution has a standard deviation of ${\sigma }_{\mathrm{PCA}}\approx 3.0 \%$ when all 10 components are used to reconstruct the continuum spectra (blue histogram). Fewer components, i.e., ${N}_{\mathrm{PCA}}=2$ (gray histogram) or ${N}_{\mathrm{PCA}}=4$ (yellow histogram), account for less individual variance in the spectra and thus decrease the precision of the estimated continuum and increase the width of the distribution of continuum residuals. The middle and lower panels of Figure 4 plot the mean ${\mu }_{\mathrm{PCA}}$ and standard deviation ${\sigma }_{\mathrm{PCA}}$ of the relative continuum error as a function of ${N}_{\mathrm{PCA}}$ used to reconstruct the continuum. Increasing the number of PCS implies a more precise, i.e., decreasing ${\sigma }_{\mathrm{PCA}}$, and less biased, i.e., decreasing $| {\mu }_{\mathrm{PCA}}|$, estimation of the quasar continua in the Lyα forest region.

Zoom In Zoom Out Reset image size

In order to calculate the PDF of the transmitted flux in the Lyα forest, the flux in the quasar spectra has to be normalized by its continuum level. However, the exact placement of the continuum level is very challenging because, particularly at higher redshifts, much of the continuum flux is absorbed by the Lyα forest. Thus one should consider the finite precision with which the continuum can be estimated. If one instead adopts a PCA model, as is done here, this continuum model has limited precision due to the finite number of PCS, and thus additionally has an intrinsic uncertainty which we quantified in the previous subsection (see Figure 4).

We will now analyze the effects of a misplaced or imprecise continuum level and the effects of noise in the spectra on the shape of the PDF. In the presence of noise and continuum errors the continuum normalized transmitted flux, which we denote as f, can take on values below 0 and above 1. We describe the level of uncertainty in the continuum normalized flux caused by the continuum error with ${\sigma }_{C}$. Note that we expect ${\sigma }_{C}$ to take on values comparable to or slightly larger than the intrinsic PCA continuum model error ${\sigma }_{\mathrm{PCA}}$ (see Figure 4), since ${\sigma }_{\mathrm{PCA}}$ represents the best continuum error one can achieve using the given PCA basis as the continuum model, but the continuum error ${\sigma }_{C}$ could generally be larger. We can approximate the transmitted flux f by modifying the perfectly normalized and noiseless flux F by adding white noise with a standard deviation ${\sigma }_{\mathrm{noise}}$ and introducing the effects of the continuum errors in the following way:

$$\begin{eqnarray}&&f=\displaystyle \frac{F+{ \mathcal N }(0,1){\sigma }_{\mathrm{noise}}}{1+{ \mathcal N }(0,1){\sigma }_{C}},\end{eqnarray} \tag{ 11 }$$

where ${ \mathcal N }(0,1)$ represents a random draw from a normal Gaussian distribution with mean $\mu =0$ and standard deviation $\sigma =1$. Note that the numerator in Equation (11) adds random white noise to the perfectly normalized transmitted flux, whereas the denominator introduces the effects of continuum errors.

In Figure 4 we have seen that the relative continuum errors ${\rm{\Delta }}C/C$ of the ensemble of all 50 quasars follows a Gaussian distribution. Hence, we can generate PDF models that incorporate the uncertainty in the continuum model by simply drawing Gaussian deviates with a standard deviation of the continuum error ${\sigma }_{C}$ and adding them to each flux pixel in the Lyα forest to mock up the continuum error. In this way we generate a new set of models for the transmitted Lyα forest flux PDF that incorporate the continuum errors as well as the effects of spectral noise. This new set of PDF models now has two free parameters: the slope parameter γ of the temperature–density relation of the IGM, where we consider values of $0.0\leqslant \gamma \leqslant 2.0$ in steps of ${\rm{\Delta }}\gamma =0.01$ for our PDF models and, additionally, the continuum error ${\sigma }_{C}$ that we allow to vary between $0.5 \%$ and $15.0 \%$ in steps of ${\rm{\Delta }}{\sigma }_{C}=0.1 \%$. The grid of PDF models is generated by averaging over many realizations of approximately normalized mock Lyα forests for every combination of the two free parameters γ and ${\sigma }_{C}$.

In Figure 5 we demonstrate the effects of continuum error and noise on the flux PDF at redshifts z = 3 (left panels) and z = 5 (right panels), for a fixed value of $\gamma =1.0$. The different colored curves represent PDF models constructed from mock Lyα forest spectra with different levels of continuum error added. The black curve shows the PDF taken from perfect mock Lyα forests without any continuum error or noise. There is an obvious difference between the PDFs at redshift z = 3 and those at z = 5. Whereas at redshift z = 3, one can clearly distinguish the different PDF models depending on the underlying continuum error, these differences are hardly noticeable at redshift z = 5. The continuum errors have only a weak impact on the shape of the PDF at the higher redshift, which can be understood as follows.

Zoom In Zoom Out Reset image size

The PDF is a function of the approximately normalized Lyα forest flux f, which itself is a function of the continuum error as shown in Equation (11). Thus:

$$\begin{eqnarray}&&{\rm{\Delta }}\mathrm{PDF}=\displaystyle \frac{d\mathrm{PDF}}{{df}}\displaystyle \frac{{df}}{d{\sigma }_{C}}{\rm{\Delta }}{\sigma }_{C}.\end{eqnarray} \tag{ 12 }$$

The first term $\tfrac{d\mathrm{PDF}}{{df}}\approx \tfrac{d\mathrm{PDF}}{{dF}}$ indicates that the change in the PDF is dependent on the slope of the PDF itself, thus the steeper the PDF, the larger the change with varying ${\sigma }_{C}$. This can be clearly seen in the PDF at z = 3 (left panel of Figure 5): at fluxes around $0.8\lesssim f\lesssim 1.0$ the PDF has a steep slope such that the term $\tfrac{d\mathrm{PDF}}{{df}}$ is large, as are the differences in the PDF models for the different values of ${\sigma }_{C}$, whereas at smaller flux values of $0.2\lesssim f\lesssim 0.8$ the slope of the PDF $\tfrac{d\mathrm{PDF}}{{df}}$ is smaller, as are the differences in the PDF models. But that explanation is not sufficient, because the PDF is also steep for very small flux values around $f\approx 0.1$ at z = 3 as well as at z = 5. At these small fluxes there is no change in the PDF models visible for different values of ${\sigma }_{C}$ because the second term of Equation (12) $\tfrac{{df}}{d{\sigma }_{C}}=\left|\tfrac{f}{1+{ \mathcal N }(0,1){\sigma }_{C}}\right|$ that indicates the change in the PDF is also dependent on the actual flux level f. Hence for small fluxes $f\approx 0.1$ there is a much smaller change in the PDF than at larger flux values $f\approx 1.0$. Thus at z = 5 where the flux values are lower overall (the mean flux is $\langle F\rangle =0.189$), there is very little change in the PDF when the continuum error ${\sigma }_{C}$ is increased. This implies that at lower redshifts the PDF is very sensitive to continuum errors, whereas at higher redshift the sensitivity is much lower, which is illustrated in Figure 5.

In order to estimate model parameters from measurements of the flux PDF, we introduce a likelihood function ${{ \mathcal L }}_{\mathrm{PDF}}$. The PDF models are explicitly dependent on two free parameters, the physical parameter γ and the continuum error ${\sigma }_{C}$. In addition there is an implicit dependence on the parameters of the quasar continuum model, i.e., the coefficients ${\alpha }_{{ij}}$ of the PCS for each quasar (see Equation (8)), since each spectrum needs to be divided by its estimated continuum level in order to normalize the transmitted Lyα forest flux before calculating the flux PDF. Thus the likelihood function has a multi-variate Gaussian form and is given by the following equation:

$$\begin{eqnarray}{{ \mathcal L }}_{\mathrm{PDF}} & = & { \mathcal L }(\gamma ,{\sigma }_{C},{\alpha }_{{ij}}| {\mathrm{PDF}}_{\mathrm{data}})\\ & = & \displaystyle \frac{1}{\sqrt{\det ({\boldsymbol{C}})}}\exp \left(-\displaystyle \frac{1}{2}{{\boldsymbol{d}}}^{{\rm{T}}}{{\boldsymbol{C}}}^{-1}{\boldsymbol{d}}\right).\end{eqnarray} \tag{ 13 }$$

Here ${\boldsymbol{C}}$ represents the covariance matrix governing the covariances between the different flux bins in the PDF. The vector ${\boldsymbol{d}}$ denotes the difference between the calculated flux PDF and the model PDF for each bin:

$$\begin{eqnarray}&&{\boldsymbol{d}}={\mathrm{PDF}}_{\mathrm{data}}-{\mathrm{PDF}}_{\mathrm{model}}(\gamma ,{\sigma }_{C},{\alpha }_{{ij}}).\end{eqnarray} \tag{ 14 }$$

Our new Bayesian formalism is based on this likelihood function ${{ \mathcal L }}_{\mathrm{PDF}}$ which we sample with an MCMC algorithm. The goal of an MCMC is to draw samples from a complex probability distribution for which direct sampling is challenging, thus providing an approximate representation of the true posterior probability distribution of the model parameters. For our MCMC algorithm we employ the python implementation emcee, which is described in detail by Foreman-Mackey et al. (2013). An important advantage of this algorithm compared to other codes is its affine-invariance property, first proposed by Goodman & Weare (2010). This makes the algorithm insensitive to any covariances among the parameters and requires hand-tuning of only one or two parameters compared to $\sim {N}^{2}$ for a traditional algorithm in an N-dimensional parameter space. Another useful property of the emcee implementation is that the parameter space is explored by a set of chains, so-called walkers, that evolve in parallel as an ensemble, rather than by a single chain. At every step each walker is randomly assigned to a partner walker, moving along lines which connect the single chains to each other. This feature makes the exploration of the parameter space extremely efficient and allows for parallel computing on multi-core processors.

There are various sources of covariance in the flux PDF. The flux in the Lyα forest contains spatial correlations that give rise to covariances in the flux PDF. The continuum residuals due to continuum estimation errors of the quasar continua will also result in non-diagonal entries in the covariance matrix of the flux PDF. In what follows we will investigate a set of models for the quasar continua that take different numbers of PCS into account, i.e., we vary ${N}_{\mathrm{PCA}}=0,1,2,3,4$. For each of the ${N}_{\mathrm{PCA}}$ that we explore we generate a separate covariance matrix.

To generate the covariance matrix we create a large set of mock quasar spectra as described in Section 2. For each mock spectrum, the underlying continuum, which is a draw from the 50 HST spectra, is fit with ${N}_{\mathrm{PCA}}$ PCS in the absence of Lyα forest absorption (see Section 3.1). The mock spectrum, which is a product of the true HST continuum times a realization of the Lyα forest, is divided by the approximate PCA continuum model, thus giving approximately normalized Lyα forest spectra with a continuum error that is intrinsic to the continuum model itself. We create many ensembles of ${N}_{\mathrm{QSO}}$ quasar spectra generated in this way, and calculate the PDF of each realization. We use this ensemble of PDFs to calculate the covariances between the different PDF flux bins and obtain a covariance matrix. This covariance matrix now contains the continuum error that is intrinsic to the continuum model with the chosen number of PCS. We choose a fixed value for $\gamma =1.0$, which we set to be our fiducial value for the mock data set, when generating the covariance matrices.

Note that we consider bins in the PDF in the flux range of $-0.5\leqslant f\leqslant 1.5$. In the absence of continuum error and noise only the bins between $0.0\leqslant f\leqslant 1.0$ would be populated. The inclusion of continuum error and noise will thus populate bins below 0 and above 1, but for any given γ and ${N}_{\mathrm{PCA}}$ not all flux bins between $-0.5\leqslant f\leqslant 1.5$ will be populated. Indeed the very low and very high end of this range rarely contain any flux pixels, and for a finite number of samples these bins may not be populated, such that the resulting covariance matrix is zero. However, models with values of γ significantly different from those used to generate the covariance matrix may actually populate these bins, such that we cannot just neglect the bins and reduce the allowed flux range. The correct approach would clearly be to adopt a more complicated procedure and compute a different covariance matrix for each model, rather than simply fixing γ and ${N}_{\mathrm{PCA}}$. However, we lack the ability to generate covariance matrices for models with different values of ${\sigma }_{C}$, other than those intrinsic to the continuum model, i.e., ${\sigma }_{\mathrm{PCA}}$, since by construction adding a random level of continuum error as described in Equation (11) has no covariance. In this way our covariance matrix is approximate, because we adopt a model with fixed ${N}_{\mathrm{PCA}}$ and an implied single value of ${\sigma }_{C}={\sigma }_{\mathrm{PCA}}$ to calculate it. This means that the range of values of ${\sigma }_{C}$ we can explore to create covariance matrices is limited and thus results in matrix entries with zero covariance at the edges of the considered flux range.

For simplicity, we opt to proceed with a single covariance matrix (for each employed value of ${N}_{\mathrm{PCA}}$). In order to prevent the covariance matrix from becoming singular, we apply a technique known as matrix shrinkage and add an identity matrix multiplied with a small value (e.g., Schaefer & Strimmer 2005) to our covariance, thus setting a floor on the diagonal elements. The size of the floor should be roughly the same size as the smallest calculated covariance and, after some experimentation, we set this value to ${10}^{-4}\times ({N}_{\mathrm{QSO}}/10)$, where ${N}_{\mathrm{QSO}}$ is the number of quasars in the data set that is analyzed. We will come back to this subtle point in Section 4 when discussing the effects of this floor on our results. Figure 6 shows the covariance and correlation matrices of the PDF for an ensemble of 10 quasar spectra at redshift z = 3 generated for a continuum model containing two PCS.

Zoom In Zoom Out Reset image size

We have shown previously that we will have to take continuum errors into account when analyzing the flux PDFs of the normalized quasar spectra. We have also shown how we modify the covariance matrix (Section 3.4) and PDF models (Section 3.2) to include the continuum error due to the intrinsic limitations of the continuum model. We will now demonstrate that these steps are necessary to obtain unbiased estimates of the thermal parameter γ that is our main interest. Note that—in this subsection only—we do not use our full PDF regulated continuum fitting algorithm for illustrating the importance of incorporating the continuum errors, but for the sake of simplicity we conduct an exercise with a simplified version of our method that only estimates the parameter γ and ${\sigma }_{C}$, but no continuum parameters. Rather, we analyze already normalized quasar spectra that have been divided by their best-fit continuum model with ${N}_{\mathrm{PCA}}=2$ that was estimated when no Lyα forest absorption was added to the spectra (see Section 3.1), thus resulting in approximately normalized Lyα forest spectra.

Our MCMC-based algorithm produces posterior probability distributions for the free parameters that we would like to estimate, in this particular case only γ and ${\sigma }_{C}$. We will now examine three different scenarios that we will explain in the following paragraph.

Figure 7 shows the resulting posterior probability distribution of the parameter γ when running our simplified algorithm for the three different situations for an ensemble of 10 quasar spectra. The right and left panels of the figure correspond to the analysis at redshift z = 3 and z = 5, respectively. The histogram shown in blue represents the first situation which is the ideal case, where there is neither spectral noise nor continuum error, i.e., ${\sigma }_{C}=0$, included in the Lyα forest spectra. Also the PDF models for this ideal case and the covariance matrix assume noise-free spectra without continuum error and thus only contain entries between $0\leqslant F\leqslant 1$. The covariance matrix used in this case was calculated from multiple PDFs assuming noise-free, perfectly continuum normalized spectra. In this ideal situation our algorithm produces an unbiased estimate for the parameter γ for both redshifts and recovers the fiducial input value of $\gamma =1.0$ shown as the black vertical line. The posterior probability distribution is not exactly centered around the fiducial value which is due to a random fluctuations in the realization of a finite number of quasar spectra.

Zoom In Zoom Out Reset image size

Of course this first scenario is very unrealistic, because real quasar spectra contain noise and we do not know their continuum level precisely and hence do not have perfectly normalized Lyα forests. Thus in the second scenario, shown as the yellow histogram, spectral noise and continuum error with ${\sigma }_{C}={\sigma }_{\mathrm{PCA}=2}\approx 8.4 \%$ are added to the Lyα forest spectra by dividing each quasar spectrum by its continuum model with ${N}_{\mathrm{PCA}}=2$. However, in the model fitting, we still use the PDF model and the covariance matrix corresponding to the ideal case of noise-free spectra with no continuum error. Thus the PDF of the approximately normalized Lyα forest spectra now contains flux pixels below $f\lt 0$ and above $f\gt 1$, but then all pixels outside the flux interval $f\in [0,1]$ are neglected in the likelihood function. Figure 7 shows that neglecting the noise and continuum error in the modeling clearly results in a very biased estimate for the parameter γ, with the true fiducial value of $\gamma =1.0$ lying far outside the resulting posterior probability distribution (yellow histogram).

This example demonstrates the need to incorporate the effects of spectral noise and continuum error into the PDF models and the covariance matrix. Hence we enlarge the range of fluxes that we consider and add these effects as described in Sections 3.2 and 3.4 to the models for the PDF and the covariance matrices (for ${N}_{\mathrm{PCA}}=2$, i.e., ${\sigma }_{C}={\sigma }_{\mathrm{PCA}=2}\approx 8.4 \%$). We again estimate the posterior distribution for γ, now marginalizing over ${\sigma }_{C}$ as a nuisance parameter. The Lyα forest spectra still contain noise and continuum errors but these effects are now included in the models of the flux PDF and the covariance matrix. The posterior probability distribution for γ in this third scenario is shown in Figure 7 as the gray histogram. The bias in the γ estimate that arose in the previous case is now removed, allowing us to recover the fiducial input value, albeit with slightly lower precision than in the idealized first case.

Incorporating continuum error into our PDF models introduces the free parameter ${\sigma }_{C}$ into our formalism. Thus the posterior probability distribution for γ in this last case is marginalized over this parameter. This marginalization and the reduced precision due to spectral noise causes the posterior probability distribution to broaden, which can be seen when comparing the blue and gray histograms. The resulting reduction in precision is greater at z = 3, but less obvious at z = 5, because continuum errors have a much smaller impact on the PDF models at higher redshift (see Figure 5). Indeed, the bias that arose in the second scenario when adopting the idealized (noiseless spectra and perfect continuum) models and covariance matrix at redshift z = 5 were largely due to the effects of ignoring the noise in modeling, rather than the continuum error.

Figure 8 shows the PDF models corresponding to the three different scenarios described in Figure 7. The black curves in each panel show the true PDF that ideally should have been recovered by our algorithm, whereas the colored curves show the actual recovered PDF. In the first and third case the recovered PDF models (blue dashed and gray curves) agree with the true PDF models (black curves) very well and lead to an unbiased estimate of the parameter γ (see Figure 7). The middle panels show the second scenario, where the quasar spectra contain noise and continuum error and thus the resulting PDFs (yellow curves) contain flux values outside of the considered flux range of $0\leqslant f\leqslant 1$, that are indicated by the black dashed lines. The resulting PDF models differ from the true PDF models (black curves) and prefer PDF models that correspond to biased values for γ (dotted black curves). Considering a wider range of flux values and modeling the effects of noise and continuum error in the PDFs and covariance matrices (bottom panels) resolves this discrepancy between the true and the recovered PDF models and the fiducial value of γ can be recovered again.

Zoom In Zoom Out Reset image size

We now have the tools to obtain unbiased estimates for thermal parameters of the IGM such as γ and can also estimate the uncertainty ${\sigma }_{C}$ inherent to imperfect continuum fits. In the following section we will apply our full method to ensembles of mock quasar spectra.

We investigate the efficacy of our algorithm for an ensemble of 10 high-resolution quasar spectra, i.e., ${N}_{\mathrm{QSO}}=10$, which is a realistic ensemble size corresponding to the path length of a typical redshift bin given current high-resolution spectroscopic quasar samples (M. Walther et al. 2017, in preparation). The continuum models that we analyze have different numbers of PCS, namely ${N}_{\mathrm{PCA}}=0,1,2,3,4$. No principal components ${N}_{\mathrm{PCA}}=0$, means that we divide each quasar by the mean quasar spectrum from Suzuki (2006) and only fit for the PDF parameters γ and ${\sigma }_{C}$. Thus the largest parameter space we are considering corresponds to ${N}_{\mathrm{QSO}}=10$ quasar spectra with each continuum modeled with ${N}_{\mathrm{PCA}}=4$, which results in a 42-dimensional parameter space. Due to the long computational time that it takes until convergence is achieved, we did not consider continuum models with more components. Also, Figure 4 shows that we do not expect a big improvement in the precision of the continuum estimation for continuum models with more than four PCS.

We show the results for a realization of an ensemble of ${N}_{\mathrm{QSO}}=10$ quasar spectra and a continuum model with ${N}_{\mathrm{PCA}}=2$ in Figure 10. We run our MCMC algorithm with 100 walkers for 10,000 steps and choose half the steps as the burn-in time until convergence of the MCMC chains is achieved. The second half of the chains is used as posterior probability distributions. The upper panels show the posterior probability distribution for the parameter γ at redshift z = 3 and z = 5 in the left and right panels, respectively. In both cases the fiducial input parameter of $\gamma =1.0$ with which the ensemble of mock quasar spectra was generated is recovered. The width (i.e., 68% region) of the posterior probability distribution, which gives the precision with which we can recover γ with this ensemble of 10 quasar spectra, is $\pm 8.6 \%$ at redshift z = 3 and $\pm 6.1 \%$ at z = 5. Note that we quote the average between the (50th–16th)-percentile and the (84th–50th)-percentile here, since the distributions are fairly symmetric.

Zoom In Zoom Out Reset image size

The middle panels show the posterior probability distribution of the continuum error ${\sigma }_{C}$ in the estimation of the continuum model of all quasar spectra in the ensemble. At redshift z = 3 the ensemble of 10 quasar spectra that we chose in this run constrains the continuum error to be ${\sigma }_{C}\approx 9.8\pm 2.1 \%$. At redshift z = 5 the continuum error ${\sigma }_{C}$ is only mildly constrained. This is not a surprising result given that the flux PDF at this high redshift is fairly insensitive to the level of uncertainty in the continuum estimation (see Figure 5). In this case the median value of ${\sigma }_{C}\approx 12.5\pm 2.2 \%$ should be considered with caution since it will also be dependent on prior assumptions for ${\sigma }_{C}$.

The lower panels of Figure 10 show the actual continuum residuals ${\rm{\Delta }}C/C$ for all quasars in the ensemble calculated with Equation (10), i.e., they show the relative continuum error between the true quasar continuum and the continuum model estimated via PDF regulation. The obtained median values of ${\rm{\Delta }}C/C\approx -2.1 \%$ at z = 3 and ${\rm{\Delta }}C/C\approx -2.8 \%$ at z = 5 indicate that there is a very small systematic overestimation of the continuum level. The width of this residual distribution, the actual PDF regulated continuum error ${\sigma }_{{\rm{\Delta }}C/C}$, i.e., the average of the (84th–50th) percentile and (50th–16th) percentile, shows that the actual continuum error for this ensemble of 10 quasar spectra is ${\sigma }_{{\rm{\Delta }}C/C}\approx 6.7 \%$ at z = 3 and ${\sigma }_{{\rm{\Delta }}C/C}\approx 10.5 \%$ at z = 5. These values are slightly lower than the estimated continuum error ${\sigma }_{C}$, i.e., ${\sigma }_{C}$ is slightly overestimated, but still within the $1\sigma$ region of the posterior probability distribution at z = 5 and slightly outside of it ($1.46\sigma$) at z = 3.

By means of the posterior probability distributions we can explore the degeneracy between the two parameters that govern the shape of the PDF. This is demonstrated in Figure 11 at redshift z = 3 and z = 5 in the left and right panels, respectively. The slope parameter γ is shown on the x-axis versus the continuum error ${\sigma }_{C}$ on the y-axis. The gray dashed lines indicate the fiducial parameter of $\gamma =1.0$ and the intrinsic error of the continuum model with two PCS, i.e., ${\sigma }_{C}={\sigma }_{\mathrm{PCA}=2}\approx 8.4 \%$ (see Figure 4).

Zoom In Zoom Out Reset image size

At redshift z = 3 we can obtain constraints for both parameters with our algorithm, since there is no degeneracy between the two parameters of interest. This result can be seen directly in Figures 1 and 5, since γ and ${\sigma }_{C}$ alter the shape of the PDF in different ways. At higher redshift at z = 5 the constraints on ${\sigma }_{C}$ are very weak, because the PDF of the transmitted flux in the Lyα forest is less sensitive to this parameter, as we have seen in Figure 5. However, at both redshifts we can constrain the slope parameter γ of the temperature–density relation of the IGM, which was our main goal.

In this subsection we investigate the robustness and reliability of our method for ensembles of 10 high-resolution quasar spectra. Hence we would like to determine whether our method results in possible biases in the model parameters, whether the estimated continuum error ${\sigma }_{C}$ is a good proxy for the actual PDF regulated continuum error ${\sigma }_{{\rm{\Delta }}C/C}$, and how well the continuum error compares to the intrinsic accuracy of the chosen PCA basis. Therefore we run our algorithm 100 times, each with a different realization of an ensemble of 10 quasar spectra. Each of the 100 runs entails an MCMC, which results in posterior probability distributions for γ and ${\sigma }_{C}$ as well as for the coefficients of the PCS for each quasar. Taking the mean of each posterior probability distribution of the PCS gives us a mean continuum model for each quasar spectrum, which enables us to calculate the continuum residual distribution ${\rm{\Delta }}C/C$ and its width, the actual PDF regulated continuum error ${\sigma }_{{\rm{\Delta }}C/C}$.

The results of this investigation are presented in Figure 12 for redshift z = 3 and z = 5 in the left and right panels, respectively. The upper panels show the estimated values for the slope parameter γ of the temperature–density relation as a function of ${N}_{\mathrm{PCA}}$ that were used to model the continuum of each quasar. From each of the 100 runs we obtain a posterior probability distribution for γ and adopt its median as the best estimate (see the upper panels of Figure 10). We then examine this distribution of estimates of γ from the 100 realizations of the quasar ensemble. The data points in the upper panels of Figure 12 indicate the median of this distribution, i.e., the 50th percentile, and the errorbars show the 16th and 84th percentile. In order to compare to the ideal but unrealistic situation when assuming ensembles of 10 noise-free spectra with no continuum error we also show the gray areas that show the median (gray dashed lines) and the 16th to 84th percentile region of estimates of γ. In these cases γ is the only free parameter to estimate with MCMC.

Zoom In Zoom Out Reset image size

The upper panels of Figure 12 show that we can recover the thermal parameter γ at both redshifts with an ensemble of 10 quasar spectra remarkably well. The medians of the posterior probability distributions for γ are mostly within the 16th to 84th percentile region of the estimates for ensembles of perfect noise-free spectra, and the fiducial input value of $\gamma =1.0$ lies within the 68% region of all estimates. For ${N}_{\mathrm{PCA}}=2$ which we will adopt as our fiducial model, the median estimate for γ and its $1\sigma$ fluctuations at z = 3 are $\gamma \approx 1.013$ with ${\sigma }_{\gamma }=8.3 \%$ and $\gamma \approx 1.029$ with ${\sigma }_{\gamma }=5.6 \%$ at z = 5. We expect $1\sigma$ fluctuations around the fiducial value of $\gamma =1.0$ of ${\sigma }_{\mathrm{med}\gamma }\simeq 8.3 \% /\sqrt{100}\approx 0.83 \%$ at z = 3 (${\sigma }_{\mathrm{med}\gamma }\simeq 5.6 \% /\sqrt{100}\approx 0.56 \%$ at z = 5) due to the fact that we only investigate 100 different realizations of our ensemble. Hence the estimated values of γ lie marginally outside the expected region. This may be a random fluctuation, however, since the other estimates of γ for different values of ${N}_{\mathrm{PCA}}$ all fluctuate in the same direction, which is suggestive that the values are slightly biased. Nevertheless the bias in the estimation of γ is still 6–7 times smaller at z = 3 (2–3 times smaller at z = 5) than its $1\sigma$ errorbars and thus we can tolerate this possible small bias. The marginalization over noise and continuum error leads to a modest reduction of precision in the measurement of γ, since the errorbars increase by a factor of ∼2 compared to the ideal noise-free case (gray areas).

The middle panels of Figure 12 illustrate how centered the continuum residuals obtained by PDF regulation are about zero. For each MCMC inference we compute the continuum residuals ${\rm{\Delta }}C/C$ in the region of the Lyα forest (i.e., rest-frame wavelengths 1040 Å $\leqslant \ \lambda \ \leqslant 1190$ Å) with Equation (10), and compute the median ${\mu }_{{\rm{\Delta }}C/C}$ of the distribution. This results in a distribution of 100 median values ${\mu }_{{\rm{\Delta }}C/C}$, and the median of this distribution of medians (points), as well as its 16th and 84th percentiles (errorbars), are plotted in the middle panels. The blue dashed lines and shaded regions indicate the corresponding ${\mu }_{{\rm{\Delta }}C/C}$ resulting from errors intrinsic to the PCA continuum model. Specifically, we MCMC fit only the continua of the the same 100 realizations of ${N}_{\mathrm{QSO}}=10$ quasars with the PCA model, and perform exactly the same procedure on the resulting ${\rm{\Delta }}C/C$ distributions, and 100 ${\mu }_{{\rm{\Delta }}C/C}$ median values.

A completely unbiased continuum estimation would give a continuum residual distribution with a median of ${\mu }_{{\rm{\Delta }}C/C}\approx 0$. We see that both at z = 3 and at z = 5, there is a small negative bias of a few percent in the estimates, i.e., we slightly overestimate the quasar continuum level. However, this bias is below 2% at z = 3 and below 3% at z = 5 for ${N}_{\mathrm{PCA}}=4$ and even smaller for ${N}_{\mathrm{PCA}}\lt 4$, and thus insignificant relative to the typical continuum residuals ${\rm{\Delta }}C/C$ of ∼7%–10% (lower panel, Figure 10). For ${N}_{\mathrm{PCA}}=2$, which represents our fiducial model, our algorithm overestimates the continuum by $\approx 1.5 \%$ at z = 3 and $\approx 2.4 \%$ at z = 5. Our measured ${\mu }_{{\rm{\Delta }}C/C}$ are marginally consistent, i.e., the 68% regions overlap, with the continuum residuals intrinsic to the continuum model, i.e., the shaded blue regions. For ${N}_{\mathrm{PCA}}\geqslant 2$ our PDF regulated continuum fits are just slightly more biased than that obtained by fitting the PCA model to the quasar spectra in the absence of the Lyα forest. However, the biases we obtain are much smaller than those estimated for hand-fitted quasar continua (see, e.g., Faucher-Giguère et al. 2008).

The lower panels of Figure 12 illustrate the behavior of the $1\sigma$ error of our continuum estimates. Recall that we characterize continuum errors in two different ways: in one we fit for the continuum error ${\sigma }_{C}$ and obtain an estimate via MCMC; in the other case we calculate the continuum residuals ${\rm{\Delta }}C/C$ using our knowledge of the the true continua of the mock spectra, and take the width of this distribution ${\sigma }_{{\rm{\Delta }}C/C}$ as the estimate for the continuum error. The black points with error bars represent the median and the 16th and 84th percentiles of the inferred continuum error ${\sigma }_{C}$ from our MCMC PDF continuum regulation. Specifically, for each of the 100 realizations of 10 quasar spectra, we performed MCMC PDF regulation, obtained a posterior probability distribution for ${\sigma }_{C}$ and adopted its median value as the best estimate. This results in 100 median ${\sigma }_{C}$ values, where the median and the 16th and 84th percentiles (errorbars) are shown in black. The blue data points correspond to the actual PDF regulated continuum error, i.e., the measured dispersion ${\sigma }_{{\rm{\Delta }}C/C}$ of the continuum residuals ${\rm{\Delta }}C/C$. For each run we compute the continuum residuals ${\rm{\Delta }}C/C$ for all pixels for the whole ensemble of quasars and take the average of the 16th and 84th percentile of this distribution as the estimate for ${\sigma }_{{\rm{\Delta }}C/C}$. We obtain 100 of these measurements for each realization of the quasar ensemble and plot the median value and the 16th and 84th percentiles. The blue shaded areas show the same distributions of median values for the ideal case with the noise-free spectra and no Lyα forest, i.e., the error intrinsic to the continuum model for a given ${N}_{\mathrm{PCA}}$.⁵

Comparing the actual error of the PDF regulated continuum residuals (blue points) with the intrinsic precision of the PCA basis (blue shaded region) in the lower left panel at z = 3, one sees that our algorithm recovers the continuum consistent with the best accuracy achievable with the PCA basis, at least for ${N}_{\mathrm{PCA}}\leqslant 2$. For higher values of ${N}_{\mathrm{PCA}}\gt 2$ the blue points lie slightly above the intrinsic accuracy of the PCA model shown by the shaded regions, and thus our PDF regulated continuum fits perform slightly worse than the underlying PCA model. One might expect the PDF regulated continuum error (blue points) to continue to decrease as ${N}_{\mathrm{PCA}}$ is increased, following the intrinsic error of the PCA basis (shaded region), but instead the PDF regulated continuum errors saturate around a value of ${\sigma }_{{\rm{\Delta }}C/C}=0.07$. We believe that this is due to the floor that we introduced into our covariance matrix (see Section 3.4) in order to prevent it from becoming singular. This floor sets a threshold below which we can no longer measure continuum errors. Once the chosen PCA basis becomes accurate enough, we can no longer recover the small continuum errors, since the floor reduces the sensitivity particularly in the high flux bins of the PDF, which are most sensitive to the continuum error (see Figure 5).

Our results at z = 5 in the lower right panel of Figure 12 tell a similar story. The PDF regulated continuum error (blue points) do not follow the intrinsic error of the PCA basis (shaded region) for ${N}_{\mathrm{PCA}}\geqslant 2$. This most likely results from the reduced sensitivity of the PDF at z = 5 to the continuum (see Figure 5), and possibly also due to the floor we imposed in the covariance matrix.

A comparison between the black data points showing the inferred continuum error ${\sigma }_{C}$ and the blue data points indicating the PDF regulated continuum error ${\sigma }_{{\rm{\Delta }}C/C}$ at z = 3 reveals that our algorithm gives a fairly good estimate for the actual continuum error. Nevertheless the inferred values are slightly biased. We underestimate the continuum error for ${N}_{\mathrm{PCA}}=0$ and ${N}_{\mathrm{PCA}}=1$ and slightly overestimate it for ${N}_{\mathrm{PCA}}\geqslant 2$. However, the inferred values of ${\sigma }_{C}$ usually lie within the 16th and 84th percentile region of the real measured value ${\sigma }_{{\rm{\Delta }}C/C}$ . The errorbars on the ${\sigma }_{C}$ estimates are very small, indicating that the MCMC always converges to the same range of values. We believe that this is also an artifact resulting from the floor imposed on our covariance matrix.

At redshift z = 5 the algorithm seems to recover the actual continuum error better, i.e., the inferred continuum error ${\sigma }_{C}$ (black data points) and real measured PDF regulated continuum error ${\sigma }_{{\rm{\Delta }}C/C}$ (blue data points) agree fairly well, but since the locations of the black inferred values are not well constrained because the posterior probability distribution for ${\sigma }_{C}$ is fairly flat (see middle right panel of Figure 10), this agreement is partly coincidental and dependent on our choice of the prior probability distribution on ${\sigma }_{C}$ in the MCMC.

In summary, the lower panels of Figure 12 illustrate that our PDF regulated continuum fitting method is producing continuum errors at z = 3 that roughly track the intrinsic precision of the PCA basis, since both the actual PDF regulated continuum error ${\sigma }_{{\rm{\Delta }}C/C}$ (blue points) as well as the inferred continuum error ${\sigma }_{C}$ (black points) lie within the blue shaded region for ${N}_{\mathrm{PCA}}\leqslant 2$. At redshift z = 5 the inferred continuum uncertainties are only weakly constrained but nevertheless reliable.

How many PCS should be used? At z = 3, the lower left panel shows a drop in the continuum error in both the inferred continuum error ${\sigma }_{C}$ and the real measured continuum error ${\sigma }_{{\rm{\Delta }}C/C}$ and thus an increase of the precision in the continuum model when two or more principal components are taken into account. At redshift z = 5 the drop in the continuum error is small, which is due to the fact that the PDF at this redshift barely constrains the quasar continuum level as we have seen before in Figure 5. At both redshifts the gains in precision for more ${N}_{\mathrm{PCA}}$, i.e., ${N}_{\mathrm{PCA}}=3$ or ${N}_{\mathrm{PCA}}=4$, are insignificant. Since we are only interested in the continuum residuals in the wavelength region of the Lyα forest, only the PCS that alter the shape of the quasar continuum in this region are influential for our method. The first principal component spectrum for example mainly addresses the shape of the Lyα emission line, and hence does not improve the continuum estimation in the Lyα forest region. This is why the continuum error is not reduced before taking at least two PCS, i.e., ${N}_{\mathrm{PCA}}\geqslant 2$, in the continuum model.

Thus we advocate using the continuum model with two PCS for this analysis, since we want to minimize the number of free parameters and thus the dimensionality of our parameter space which in turn limits the computational time. This implies that we have two nuisance parameters for the continuum estimation per quasar plus the global parameters that influence the shape of the PDF—${\sigma }_{C}$ and any parameters related to the physics of the IGM, which for our simplified IGM model is only γ. For an ensemble of 10 quasar spectra this results in a $2\times {N}_{\mathrm{QSO}}\,+2=22$-dimensional parameter space. However, in principle our method can accommodate larger ensembles of quasar spectra and a more complicated model of the IGM.

Although this study only considered a highly simplified single-parameter model of the thermal state of the IGM, in practice our method can easily accommodate a more complicated multi-parameter PDF model, the only limitations being computational time and effort. When applying this method to real quasar spectra one would use hydrodynamical simulations to generate the IGM model and also consider the temperature at mean density T₀, the Jeans pressure smoothing scale ${\lambda }_{{\rm{P}}}$ (Rorai et al. 2013; Kulkarni et al. 2015), and other possible nuisance parameters (Lee et al. 2015). Adding more thermal parameters or other nuisance parameters related to the IGM (e.g., Lyman limit systems) extends the computational effort only marginally. More parameters in the continuum model, on the other hand, add a more significant amount of time to the computation because they need to be estimated for each quasar in the ensemble and hence the number of continuum parameters multiplies by ${N}_{\mathrm{QSO}}$. In this work we compressed the continuum error of each quasar spectra into only two nuisance parameters. We tested our algorithm with an ensemble of 10 quasar spectra, but a larger ensemble with even double the number of quasar spectra is still easily feasible. A larger ensemble of quasar spectra would also further improve the precision on our estimates of the thermal parameters.

The mean flux of the Lyα forest has been held fixed in our analysis. Therefore, we computed 100 different realizations of the Lyα forest for each value of γ that we considered in our models and computed a scaling factor for the optical depth τ averaged over the 100 Lyα realizations, in order to match the mean flux in our mock spectra to measurements of the mean flux (see Section 2.1). While creating each mock quasar spectrum we re-scale τ with the scaling factor. In this way, the average mean flux of the Lyα forest is kept fixed, but individual mock spectra can have higher or lower mean flux values, which is expected to be the case also for real quasar spectra due to cosmic variance.

The mean flux could be considered as an additional nuisance parameter in the PDF models. This would not add significantly to the computation time (see Section 5.1), but might decrease the precision of the thermal parameter estimation slightly. However, given the precision of the most recent mean flux measurements (${\rm{\Delta }}\langle f{\rangle }_{z\sim 3}\approx 1 \%$ and ${\rm{\Delta }}\langle f{\rangle }_{z\sim 5}\approx 7 \%$ Becker et al. 2013), we opted to keep it fixed.

In our simplified model for Lyα forest absorption we do not take the redshift evolution within the Lyα forest along the line of sight into account, i.e., we assume a fixed value of γ. However, assuming a redshift evolution within the Lyα forest, and thus within the temperature–density relation of the IGM, would result in a redshift-dependent slope parameter, i.e., $\gamma (z)$. Our parameter space of the PDF models can be easily augmented with a redshift-dependent γ and we would only need to account for one more free parameter in the MCMC runs.

Rorai et al. (2016) discuss temperature–density relations with a break in the slope, i.e., a different slope for over- or underdensities. Even such nonlinear dependences can be easily accommodated in our algorithm by simply constructing a different set of PDF models that contain two different linear slope parameters γ. These models can be parameterized by two parameters ${\gamma }_{{\rm{o}}}$ and ${\gamma }_{{\rm{u}}}$ for over- and underdense regions, respectively, which simply adds another free parameter to our algorithm.

In principle the PDF models can be made arbitrarily complicated by adding more and more nuisance parameters that will be marginalized out in the end.

A more suitable set of PCS could further increase the precision of the continuum model and reduce the continuum error without increasing the number of free parameters in the continuum model and thus the computational time. The PCS that we use in this work from Suzuki (2006) cover the the whole spectral range from 1020 Å to 1600 Å of quasar spectra. However, for studies of the IGM, one is interested in a small range of wavelengths, i.e., the spectral range spanning the Lyα forest. Thus some of the PCS we use account for the different shapes in the Lyα and other emission lines, for example, but do not help improve the continuum model in the Lyα forest region. One could thus come up with a different PCA basis that would be more suited for our purposes. Other methods, such as an independent component analysis (e.g., Allen et al. 2011) or heteroscedastic matrix factorization (e.g., Tsalmantza & Hogg 2012), could also be useful to construct a more suitable continuum model while keeping the dimensionality of the PCA low.

When applying this algorithm to real quasar spectra, the chosen set of PCS might not be flexible enough to account for enough of the variance in the spectra. However, for high-z spectra, the PCA basis of Pâris et al. (2011) could be applied for instance, or a new set of PCAs could be created from a similar procedure by fitting the >10,000 spectra available with high S/N in the BOSS survey. If such an augmented and improved PCA were constructed and fit to real data, if the red side of each spectrum is modeled well by the PCA, there is no reason to believe that the blue side is not expected to match.

We used an approximation when constructing the covariance matrix, whereby we fixed the covariance matrix to have the intrinsic continuum error of a single continuum model. This fact required us to add a floor to the covariance matrix according to the matrix shrinkage approach, in order to prevent the matrix from becoming singular, since the very high and very low flux bins were rarely populated for the fiducial model used to construct the covariance. This approximation set a floor on the accuracy of our PDF regulated continuum, and limited our ability to obtain more precise continuum fits by increasing the number of PCS. Thus in the future it would be worthwhile to construct covariance matrices as a function of the model parameters γ, ${\sigma }_{C}$ and possibly others. This would remove the bias in the continuum residuals and increase the precision of the continuum estimation. However, constructing covariance matrices with varying ${\sigma }_{C}$ is tricky, since we introduced the continuum error into our PDF models by construction without any covariance (Equation (11)). One possible approach to solve this problem could be to interpolate between the covariance matrices with continuum errors intrinsic to the chosen continuum model. This will be the subject of future work.

Applying our new algorithm to real quasar spectra requires the spectra to be flux calibrated. Suzuki et al. (2003) showed how one can flux calibrate echelle spectra using lower-resolution spectra with accurate spectrophotometry. We do not expect any problems arising due to errors and uncertainties on this flux calibration, since we infer the continuum error ${\sigma }_{C}$ with MCMC, and any flux calibration errors would become a part of the continuum error that we can marginalize out.

In Section 3.1 we introduced a likelihood function ${ \mathcal L }=\exp (-{\chi }^{2}/2)$ (see Equation (9)) to find estimates for the continuum model to the data that we used to evaluate the intrinsic continuum error of the model. We will now use this likelihood function for the continuum model of each quasar on longer wavelengths than the Lyα emission line, i.e., on the red side of the spectrum. The likelihood function now reads:

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{cont}}={ \mathcal L }({\alpha }_{{ij}}| {C}_{\mathrm{data}})\\ &&=\displaystyle \prod _{i}^{{N}_{\mathrm{QSO}}}\displaystyle \prod _{\lambda }\displaystyle \frac{1}{\sqrt{2\pi {\sigma }_{i}^{2}}}\exp \left(-\displaystyle \frac{{({C}_{\mathrm{model},i\lambda }({\alpha }_{{ij}})-{C}_{\mathrm{data},i\lambda })}^{2}}{2{\sigma }_{i,\mathrm{noise}}^{2}}\right).\end{eqnarray} \tag{ 15 }$$

The free parameters of the continuum model ${\alpha }_{{ij}}$ are the coefficients of the PCS j of each quasar i as previously seen in Equation (8). Thus the likelihood function evaluates how well the continuum model given by Equation (8) describes the data on the red side of each quasar spectra.

It is important to note that only the wavelength range redward of the Lyα emission line will be taken into account to evaluate this likelihood function, i.e., 1216 Å–1600 Å, since there are no absorption lines present (we neglect possible contamination due to metals that are not incorporated into our mock data set—in real data these metals would be straightforward to mask), whereas the wavelength region blueward of Lyα will be considered by the previous likelihood function from Equation (13) as before.

Thus ${{ \mathcal L }}_{\mathrm{cont}}$ takes every pixel at wavelengths on the red side of the quasar spectrum into account, i.e., ∼40,000 data points. This likelihood function is clearly very different in character compared to the previous likelihood function ${{ \mathcal L }}_{\mathrm{PDF}}$ (Equation (13)), where we bin the data and evaluate the difference between the data and the model for every bin in the PDF and thus the likelihood function ${{ \mathcal L }}_{\mathrm{PDF}}$ only covers ∼25 data points. In the next subsection we will combine the two likelihood functions ${{ \mathcal L }}_{\mathrm{PDF}}$ and ${{ \mathcal L }}_{\mathrm{cont}}$ and will see that we have to slightly modify this new likelihood function ${{ \mathcal L }}_{\mathrm{cont}}$ due to its differences in character in order to make the two likelihood functions compatible.

We would like to join the two likelihood functions ${{ \mathcal L }}_{\mathrm{cont}}$ and ${{ \mathcal L }}_{\mathrm{PDF}}$ to make use of the information on the continuum level provided by the wavelength coverage on the red side of the Lyα emission peak, but nevertheless regulate the continuum level via the Lyα forest flux PDF on the blue side of each quasar spectrum. However, joining these two likelihood functions by multiplication is not as easy as one would naively expect. We have already mentioned that the character of the two functions is quite different: in one case we evaluate the likelihood for each pixel in the wavelength range between 1216 Å and 1600 Å, i.e., ∼40,000 data points, whereas in the other case we only evaluate the likelihood after binning the pixels, which results in more than three orders of magnitude fewer data points. We know that the mean of a ${\chi }^{2}$ distribution is roughly equal to the number of degrees of freedom, which corresponds to the number of our data points. Since both likelihood functions are of the shape ${ \mathcal L }=\exp (-{\chi }^{2}/2)$, this implies that the new likelihood function is much steeper than the other. Hence in a simple multiplication of the two likelihood functions, the steeper function will always dominate the joint likelihood function.

We will therefore modify the new likelihood function ${{ \mathcal L }}_{\mathrm{cont}}$ in order to make the two functions more compatible. We will augment the function with an additional uncertainty ${\sigma }_{\mathrm{red}}$, which reflects the inability of the red side of each spectrum to predict the continuum on the blue side. Note that this parameter is not one that we estimate with MCMC but rather one that we infer by conducting various tests.

The likelihood function ${{ \mathcal L }}_{\mathrm{cont}}$ then reads:

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{cont}}={ \mathcal L }({C}_{\mathrm{model}}({\alpha }_{{ij}})| {C}_{\mathrm{data}},{\sigma }_{\mathrm{red}})\\ &&\,=\displaystyle \prod _{i}^{{N}_{\mathrm{QSO}}}\displaystyle \prod _{\lambda }\displaystyle \frac{1}{\sqrt{2\pi {\sigma }_{i}^{2}}}\exp \left(-\displaystyle \frac{{({C}_{\mathrm{model},i\lambda }({\alpha }_{{ij}})-{C}_{\mathrm{data},i\lambda })}^{2}}{2({\sigma }_{i,\mathrm{noise}}^{2}+{\sigma }_{\mathrm{red}}^{2})}\right).\end{eqnarray} \tag{ 16 }$$

Now the two likelihood functions can be joined together by simply multiplying them:

$$\begin{eqnarray}&&{ \mathcal L }(\gamma ,{\sigma }_{C},{\alpha }_{{ij}})={{ \mathcal L }}_{\mathrm{cont}}({\alpha }_{{ij}}| {\sigma }_{\mathrm{red}})\times {{ \mathcal L }}_{\mathrm{PDF}}(\gamma ,{\sigma }_{C},{\alpha }_{{ij}}).\end{eqnarray} \tag{ 17 }$$

The larger the additional error ${\sigma }_{\mathrm{red}}$ is, the less influence has the likelihood function ${{ \mathcal L }}_{\mathrm{cont}}$ and the more influence has the likelihood function ${{ \mathcal L }}_{\mathrm{PDF}}$ in the joint likelihood function. In the extreme case of ${\sigma }_{\mathrm{red}}\to \infty$, the likelihood ${{ \mathcal L }}_{\mathrm{cont}}$ would be ignored and we are using only the blue wavelength coverage of each spectrum with ${{ \mathcal L }}_{\mathrm{PDF}}$ as we have done in the main text. In the other extreme of ${\sigma }_{\mathrm{red}}=0$, we basically fix the continuum level to the best fit of the continuum model to the red side of each quasar spectrum due to the steepness of ${{ \mathcal L }}_{\mathrm{cont}}$ and are unable to regulate the continuum on the blue side via the PDF any longer. The likelihood function ${{ \mathcal L }}_{\mathrm{PDF}}$ would then only be used to estimate γ and ${\sigma }_{C}$ with a fixed quasar continuum. A comparison of the two error terms ${\sigma }_{i,\mathrm{noise}}$ representing the Gaussian white noise that we add to each quasar spectrum with very high S/N with ${\sigma }_{\mathrm{red}}$, which reflects the limited ability to predict continuum on the blue side with an estimate of the continuum on the red wavelength side, reveals that the additional error ${\sigma }_{\mathrm{red}}$ dominates the error budget by roughly 2–3 orders of magnitude.

We illustrate the differences in the steepness of the two likelihood functions in Figure 13. We show the normalized probability distributions given by the the two functions ${{ \mathcal L }}_{\mathrm{cont}}$ and ${{ \mathcal L }}_{\mathrm{PDF}}$ dependent on one of the free parameters, the coefficient ${\alpha }_{i1}$ of the first principal component spectrum $| {\xi }_{1}\rangle$ for quasar i at redshift z = 3. All other principal component coefficients ${\alpha }_{{ij}}$ with $j\gt 1$ are thereby held fixed. The blue curve shows ${{ \mathcal L }}_{\mathrm{PDF}}$ from the Lyα forest PDF, whereas the gray curve represents ${{ \mathcal L }}_{\mathrm{cont}}$ from the continuum emission on the red side of the quasar spectrum with ${\sigma }_{\mathrm{red}}=0.0$. The red and yellow curves show also ${{ \mathcal L }}_{\mathrm{cont}}$, but now including an additional error on the red wavelength side ${\sigma }_{\mathrm{red}}$. Setting this additional error to ${\sigma }_{\mathrm{red}}=1.0$ (red curve) already causes ${{ \mathcal L }}_{\mathrm{cont}}$ to flatten and become considerably less steep. Its PDF widens compared to the gray curve with no additional error. Increasing ${\sigma }_{\mathrm{red}}$ will further flatten the likelihood function and thus increase the width of the probability distribution. This is shown as the yellow curve for ${\sigma }_{\mathrm{red}}=10.0$. This new likelihood function results now in a PDF that is of comparable width to that from ${{ \mathcal L }}_{\mathrm{PDF}}$. Multiplying these two functions ${{ \mathcal L }}_{\mathrm{PDF}}$ and ${{ \mathcal L }}_{\mathrm{cont}}$ with ${\sigma }_{\mathrm{red}}=10.0$ results in the probability distribution shown in black.

Zoom In Zoom Out Reset image size

Note that the true value for this parameter for this particular quasar is ${\alpha }_{i1}\approx -11.4$, which is strongly disfavored by fitting the pixels redward of the Lyα line if the ${\sigma }_{\mathrm{red}}$ is too small.

In Figure 14 we show the results of our analysis when including the information of each quasar spectrum redward of the Lyα emission line. We show the same quantities as in the previous Figure 12 for z = 3 and z = 5 in the left and right panels, respectively. The upper panels show the distribution of medians of 100 posterior probability distributions of γ estimated from different realizations of an ensemble of 10 quasar spectra. The middle panels show the distribution of medians of the continuum residual ${\rm{\Delta }}C/C$ distributions of each run, and the lower panels depict the continuum error—once measured as the widths of each continuum residual distribution, i.e., the actual PDF regulated continuum error ${\sigma }_{{\rm{\Delta }}C/C}$, shown as the blue data points and once inferred via MCMC as a free parameter influencing the shape of the PDF, i.e., the estimated continuum error ${\sigma }_{C}$. In this case the medians of the posterior probability distributions are used again as the best estimate. The medians and 68% regions of this distributions of best estimates are shown as the black data points.

Zoom In Zoom Out Reset image size

However, the x-axis of Figure 14 now shows the additional error ${\sigma }_{\mathrm{red}}$ that was added to the new likelihood function ${{ \mathcal L }}_{\mathrm{cont}}$. A value of ${\sigma }_{\mathrm{red}}\to \infty$ indicates that the likelihood ${{ \mathcal L }}_{\mathrm{cont}}$, and thus the red side of each quasar spectrum, has no influence in the joint likelihood function and only the spectral coverage blueward of the Lyα emission line is used to regulate the quasar continuum. This is the case we have discussed extensively in Section 4.2. Decreasing values of ${\sigma }_{\mathrm{red}}$ increase the influence of ${{ \mathcal L }}_{\mathrm{cont}}$ in the joint likelihood function. We run our algorithm as before 100 with different realizations of an ensemble of 10 quasar spectra in order to avoid statistical fluctuations due to the limited number of quasar spectra. Throughout this analysis we work with a continuum model with two PCS, i.e., ${N}_{\mathrm{PCA}}=2$. The gray and blue shaded areas in all panels depict the same ideal cases as in Figure 12, i.e., the gray shaded areas in the upper panels show the precision on γ achievable when assuming noise-free spectra and perfectly known quasar continua and the blue shaded areas in the middle and lower panels illustrate the intrinsic uncertainty of the chosen continuum model.

The behavior of our algorithm at z = 3 and z = 5 is similar. Surprisingly, the additional information of the wavelength coverage on the red side of each quasar spectrum does not seem to improve the estimation of the thermal parameter γ significantly. The opposite even seems to be the case at z = 3 for values of ${\sigma }_{\mathrm{red}}$ that are too small, i.e., ${\sigma }_{\mathrm{red}}=1.0$, where the 68% region is significantly larger than for higher values of ${\sigma }_{\mathrm{red}}$. However, the additional information on the red wavelength side does improve the quasar continuum estimation and reduces the small bias in the medians of the continuum flux residual distributions (middle panels) from ${\mu }_{{\rm{\Delta }}C/C}\approx -1.5 \%$ to ${\mu }_{{\rm{\Delta }}C/C}\approx 0.2 \%$ at z = 3 and from ${\mu }_{{\rm{\Delta }}C/C}\approx -2.3 \%$ to ${\mu }_{{\rm{\Delta }}C/C}\approx 0.3 \%$ at z = 5 for ${\sigma }_{\mathrm{red}}=10.0$ compared to the previous case from Section 4.2 when only wavelengths bluer than Lyα, i.e., ${\sigma }_{\mathrm{red}}\to \infty$, were taken into account.

The inferred estimates for the continuum error ${\sigma }_{C}$ (lower panels) also reflect this improvement in the quasar continuum model for ${\sigma }_{\mathrm{red}}=10.0$. At z = 3 the measured (blue data point) and inferred (black data point) estimates for the continuum error agree perfectly and show that our algorithm results in a reliable estimate of the actual continuum error. For lower values of ${\sigma }_{\mathrm{red}}$ the continuum level estimation is set entirely by the very steep red side continuum likelihood, which causes the continuum model on the blue wavelength side to no longer be a good fit. Hence the 68% region of the estimates for the continuum error increases in the actual PDF regulated continuum error ${\sigma }_{{\rm{\Delta }}C/C}$ as well as in the inferred values for ${\sigma }_{C}$. We mentioned already that the agreement of the measured and estimated continuum errors at z = 5 is probably only coincidental and reflects our choice of prior assumptions, since we do not have good constraints on ${\sigma }_{C}$ at these redshifts from the Lyα flux PDF.

We conclude that the information of the spectral coverage on the red wavelength side of each quasar spectrum does contain information for improving the quasar continuum model on the blue wavelength side of each spectrum. However, it is necessary to introduce an additional error ${\sigma }_{\mathrm{red}}$ in order to prevent the likelihood function ${{ \mathcal L }}_{\mathrm{cont}}$ to over-regulate the continuum model on the red side of the quasar spectrum. We experimentally determine the best value to be ${\sigma }_{\mathrm{red}}\approx 10.0$, i.e., we add an error that enlarges the noise on the red wavelength side of each quasar by multiple orders of magnitude compared to the Gaussian noise that was added to the spectra. Note that this choice for ${\sigma }_{\mathrm{red}}$ is purely based on the intent to match the widths of ${{ \mathcal L }}_{\mathrm{cont}}$ to ${{ \mathcal L }}_{\mathrm{PDF}}$, i.e., there is no physical meaning to its specific value.

However, the additional information of the red spectral coverage does not improve the estimation of the thermal parameter γ, regardless of the exact choice for the value of ${\sigma }_{\mathrm{red}}$. Since the computational time increases significantly when adding the likelihood ${{ \mathcal L }}_{\mathrm{cont}}$, we advocate using it only when interested in the continuum estimation of the quasars. If the thermal parameters are the main parameters of interest, ${{ \mathcal L }}_{\mathrm{cont}}$ provides no additional information and the computational time will be significantly shorter without including the red side likelihood function.

This study considered a highly simplified analytical model to simulate the Lyα forest absorption in high-redshift quasars. We will show now that the flux PDF of our analytical model behaves qualitatively similar to the PDF from hydrodynamical simulations. Therefore we choose different runs from the NyX simulations (Almgren et al. 2013) that differ in the slope parameter γ of the temperature–density relations, but other thermal parameters such as the temperature at mean density or the Jeans smoothing scale are kept fixed. The parameters of the chosen NyX simulations coming closest to our requirements had $\gamma =1.011$ (${T}_{0}=\mathrm{10,892}$ K, ${\lambda }_{P}=64.5$ ckpc) and $\gamma =1.575$ (${T}_{0}=\mathrm{10,753}$ K, ${\lambda }_{P}=66.5$ ckpc) at z = 3, and $\gamma =0.9436$ (${T}_{0}=7500$ K, ${\lambda }_{P}=69.66$ ckpc) and $\gamma =1.5219$ (${T}_{0}=7713$ K, ${\lambda }_{P}=67.4$ ckpc) at z = 5. Thus we compare these simulations to the PDF from the analytical models with $\gamma =1.01$ and $\gamma =1.58$ at z = 3, and $\gamma =0.94$ and $\gamma =1.52$ at z = 5. We calculate the difference in the PDF models for our analytical model and the hydrodynamical simulations and plot the relative change

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}\mathrm{PDF}}{\mathrm{PDF}}=\displaystyle \frac{\mathrm{PDF}\ (\gamma =1.58)-\mathrm{PDF}(\gamma =1.01)}{\mathrm{PDF}\ (\gamma =1.01)},\end{eqnarray} \tag{ 18 }$$

for z = 3 in the left panel of Figure 15, and correspondingly for z = 5 in the right panel. The PDF from hydrodynamical simulations differs from the analytical model at a 10%–20% level, but shows qualitatively the same behavior. Choosing a different smoothing kernel in our analytical model could reduce these differences (see Section B.2), i.e., the analytical model with slightly higher smoothing kernel will resemble more closely the PDF from the NyX simulations. Thus we do not expect any discrepancies arising in the parameter estimation with our new method due to the the simplified analytical model with which the PDF models are constructed.

Zoom In Zoom Out Reset image size

In order to test the influence of the smoothing kernel ${\sigma }_{\tau }$ that has been applied to the density field (see Section 2.1) we show in Figure 16 the relative differences between two PDF models with a smoothing kernel ${\sigma }_{\tau }=20$ km s⁻¹ and a smoothing kernel of ${\sigma }_{\tau }\approx 28$ km s⁻¹, which roughly corresponds to a change in temperature of a factor of two, i.e., ${\rm{\Delta }}T\approx 0.3$ dex. We plot

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}\mathrm{PDF}}{\mathrm{PDF}}=\displaystyle \frac{\mathrm{PDF}\ ({\sigma }_{\tau }\approx 20\,\mathrm{km}\,{{\rm{s}}}^{-1})-\mathrm{PDF}({\sigma }_{\tau }\approx 28\,\mathrm{km}\,{{\rm{s}}}^{-1})}{\mathrm{PDF}\ ({\sigma }_{\tau }\approx 20\,\mathrm{km}\,{{\rm{s}}}^{-1})}.\end{eqnarray} \tag{ 19 }$$

Changing the smoothing kernel influences the shape of the PDF at a 10% level; however, the slope parameter γ has a significantly stronger influence on its shape. In the future, the smoothing kernel ${\sigma }_{\tau }$ could be easily incorporated into our PDF models and be treated as a nuisance parameter.

Zoom In Zoom Out Reset image size