IGM CONSTRAINTS FROM THE SDSS-III/BOSS DR9 Lyα FOREST TRANSMISSION PROBABILITY DISTRIBUTION FUNCTION

BOSS (Dawson et al. 2013) is part of SDSS-III (Eisenstein et al. 2011; the other surveys are SEGUE-2, MARVELS, and APOGEE). The primary goal of the survey is to carry out precision BAO measurements at z  ∼  0.5 and z  ∼  2.5, from the luminous red galaxy distribution and Lyα forest absorption field, respectively (see, e.g., Anderson et al. 2014; Busca et al. 2013; Slosar et al. 2013). Its eventual goal is to obtain spectra of ∼1.5 million luminous red galaxies and ∼170, 000 z  >  2.15 quasars over 4.5 yr of operation.

BOSS is conducted on upgraded versions of the twin SDSS spectrographs (Smee et al. 2013) mounted on the 2.5 m Sloan telescope (Gunn et al. 2006) at the Apache Point Observatory, NM. One thousand optical fibers mounted on a plug-plate at the focal plane (spanning a 3° field of view) feed the incoming flux to the two identical spectrographs, of which 160–200 fibers per plate are allocated to quasar targets (see Ross et al. 2012; Bovy et al. 2011 for a detailed description of the quasar target selection). Both spectrographs split the light into blue and red cameras that cover 3610–10140 Å, with the dichroic overlap region occurring at around 6000 Å. The resolving power R ≡ λ/Δλ ranges from 1300 at the blue end to 2600 at the red end.

Each plate is observed for sufficiently long to achieve the S/N requirements set by the survey goals; typically, 5 individual exposures of 15 minutes are taken. The data are processed, calibrated, and combined into co-added spectra by the "idlspec2d" pipeline, followed by a pipeline that operates on the 1D spectra to classify objects and assign redshifts (Bolton et al. 2012). However, as described later in this paper, we will generate our own co-added spectra from the individual exposures and other intermediate data products.

In this paper we use data from the publicly available SDSS Data Release 9 (DR9; Ahn et al. 2012). This includes 87,822 quasars at all redshifts that have been confirmed by visual inspection as described in Pâris et al. (2012). In Lee et al. (2013), we have defined a further subset of 54,468 quasars with z_qso ⩾ 2.15 that are suitable for Lyα forest analysis, and have provided in individual FITS files for each quasar various products such as sky masks, masks for damped Lyα absorbers (DLAs), noise corrections, and continua; these are designed to ameliorate systematics in the BOSS spectra and aid in Lyα forest analysis (see Table 1 in Lee et al. 2013 for a full listing). While we use this Lee et al. (2013) catalog as a starting point, in this paper we will generate our own custom co-added spectra and noise estimates.

The typical S/N of the BOSS Lyα forest quasars is low: 〈S/N〉 ≈ 2 pixel⁻¹ within the Lyα forest; this criterion is driven by a strategy to ensure a large number of sightlines over a large area in order to optimize the 3D Lyα forest BAO analysis (McDonald & Eisenstein 2007; McQuinn & White 2011), rather than increasing the S/N in individual spectra. However, for our analysis we wish to select a subset of BOSS Lyα forest sightlines with reasonably high S/N in order to reduce the sensitivity of our PDF measurement to inaccuracies in our modeling of the noise and continuum of the BOSS spectra. We therefore make a cut on the S/N, including only sightlines that have a median of 〈S/N〉 ⩾ 6 pixel⁻¹ within the Lyα forest,¹⁸ defined with respect to the pipeline noise estimate (see Lee et al. 2013)—this selects only ∼10% of the spectra with the highest S/N. The 1041–1185 Å Lyα forest region of each quasar must also include at least 30 pixels (Δv = 2071 km s⁻¹) within one of our absorption redshift bins of 〈z〉 = 2.3, 〈z〉 = 2.6, and 〈z〉 = 3.0, with bin widths of Δz = 0.3 (see Section 3).

We discard spectra with identified DLAs in the sightline, as listed in the "DLA Concordance Catalog" used in the Lee et al. (2013) sample. This DLA catalog (W. Carithers 2014, in preparation) includes objects with column densities $N_{{\rm H\,\scriptsize{I}}}\,{>}\,10^{20}\,\mathrm{cm^{-2}}$; however, the completeness of this catalog is uncertain below $N_{{\rm H\,\scriptsize{I}}} = 10^{20.3}\,\mathrm{cm^{-2}}$. We therefore discard only sightlines containing DLAs with $N_{{\rm H\,\scriptsize{I}}} \ge 10^{20.3}\,\mathrm{cm^{-2}}$, and take into account lower column-density absorbers in our subsequent modeling of mock spectra. At the relatively high S/N that we will work with (see below), the detection efficiency of DLAs is essentially 100% (see, e.g., Prochaska et al. 2005; Noterdaeme et al. 2012) and thus we expect our rejection of $N_{{\rm H\,\scriptsize{I}}}\ge 10^{20.3}\,\mathrm{cm^{-2}}$ DLAs to be quite thorough.

Measurements of the Lyα forest transmission PDF are known to be sensitive to the continuum estimate (Lee 2012), but in this paper we use an automated continuum-fitter, MF-PCA (Lee 2012), that is less susceptible to biases introduced by manual continuum estimation. Moreover, unlike the laborious process of manually fitting continua on high-resolution spectra, the automated continuum estimation can be used to explore various biases in continuum estimation. For this purpose, we will use the same MF-PCA continuum estimation used in Lee et al. (2013), albeit with minor modifications as described in Section 3.2. We select only quasars that appear to be well described by the continuum basis templates, based on the goodness-of-fit to the quasar spectrum redward of Lyα. This is flagged by the variable CONT_FLAG =1 as listed in the Lee et al. (2013) catalog (see Table 3 in that paper). Broad absorption line (BAL) quasars, the continua of which are difficult to estimate due to broad intrinsic absorption troughs, have already been discarded from the Lee et al. (2013) sample.

Another consideration is that the shape of the transmission PDF is affected by the resolution of the spectrum, especially since the BOSS spectrographs do not resolve the Lyα forest. The exact spectral resolution of a BOSS spectrum at a given wavelength varies as a function of both observing conditions and row position on the BOSS CCDs. The BOSS pipeline reports the wavelength dispersion at each pixel, σ_disp, in units of the co-added wavelength pixel size (binned such that ln (10) Δ(λ)/λ = 10⁻⁴). This is related to the resolving power by R ≈ (2.35 × 1 × 10⁻⁴ln 10 σ_disp)⁻¹. Palanque-Delabrouille et al. (2013) have recently found, using their own analysis of the width of the arc-lamp lines and bright sky emission lines, that the spectral dispersion reported by the pipeline had a bias that depended on the CCD row and increased with wavelength, up to 10% at λ ≈ 6000 Å. We will correct for this bias when creating mock spectra to compare with the data, as described in Section 4. Figure 1 shows the (uncorrected) pixel dispersions from 236 BOSS quasars from the 〈z〉 = 2.3, S/N = 6–8 bin, as a function of wavelength at the blue end (λ = 3700–4200 Å) of the spectrograph. At a fixed wavelength, there are outliers that contribute to the large spread in σ_disp, e.g., ranging from σ_disp ≈ 0.9–1.8 at 3700 Å. We therefore discard spectra with outlying values of σ_disp based on the following criterion: we first rank-order the spectra based on their σ_disp value evaluated at the central wavelength of each PDF bin (i.e., λ = [4012, 4377, 4863] Å at 〈z〉 = [2.3, 2.6, 3.0]), and then discard spectra below the 5th percentile and above the 90th percentile. This is illustrated by the dashed red lines in Figure 1.

Zoom In Zoom Out Reset image size

Finally, since our noise estimation procedure uses the individual BOSS exposures, we discard objects that have less than three individual exposures available.

Our final data set comprises 3373 unique quasars with redshifts ranging from z_qso = 2.255 to z_qso = 3.811, and a median S/N of S/N = 8.08 pixel⁻¹. This data set represents only a small subsample of the BOSS DR9 quasar spectra, but is over two orders of magnitude larger than high-resolution quasar samples previously used for transmission PDF analysis. Table 1 summarizes our data sample, and the statistics of the redshifts and S/N bins for which we measure the transmission PDF. Figure 2 shows histograms of the pixels used in our analysis, as a function of absorption redshift.

Zoom In Zoom Out Reset image size

Since we intend to model BOSS spectra with modest S/N, we need an accurate estimate of the pixel noise that also allows us to separate out the contributions from Poisson noise due to the background and sky as well as read noise from the detector. In this subsection, we will construct an accurate probabilistic model of the flux and noise of the BOSS spectrograph, based on the individual exposure data that BOSS delivers.

The basic BOSS spectral data consists of a spectrum of each raw exposure, f_λi (inclusive of noise), an estimate of the sky s_λi, and a calibration vector S_λi, where i indicates the exposure of the n_exp exposures taken.¹⁹ The quantity s_λi is the actual sky model that was subtracted from the fiber spectra in the extraction. The calibration vector is defined as $S_{\lambda i} \equiv f_{\lambda i}/f_{N i}$, with f_Ni being the flux of exposure i in units of photoelectrons. The idlspec2d pipeline then estimates the co-added spectrum of the true object flux, $\mathcal {F}_\lambda$, from the raw individual exposures, sky estimates, and calibration vectors.

The BOSS data reduction pipeline also delivers noise estimates in the form of variance vectors, which are, however, known to be inaccurate (McDonald et al. 2006; Desjacques et al. 2007; Lee et al. 2013; Palanque-Delabrouille et al. 2013).

To quantify the fidelity of the BOSS noise estimate, we used the so-called "side-band" method described in Lee et al. (2014b) and Palanque-Delabrouille et al. (2013), which uses the variance in flat, absorption-free regions of the quasar spectra to quantify the fidelity of the noise estimate. First, we randomly selected 10,000 BOSS quasars (omitting BAL quasars) from the Pâris et al. (2012) catalog in the redshift range 1.4  ⩽  z_qso  <  3.4, evenly distributed into 20 redshift bins of width Δz_qso = 0.1 (i.e., 500 objects per bin). We then consider the flat 1460 Å  <  λ_rest  <  1510 Å spectral region in the quasar rest-frame, which is dominated by the smooth power-law continuum and relatively unaffected by broad emission lines (e.g., Vanden Berk et al. 2001; Suzuki 2006) or absorption lines. The pixel variance in this flat portion of the spectrum should therefore be dominated by spectral noise, allowing us to examine whether the noise estimate provided by the pipeline is accurate. We then evaluate the ratio of σ_side, the pixel flux rms in the rest-frame 1460 Å  <  λ_rest  <  1510 Å region divided by the average pipeline noise estimate, σ_λ:

$$\begin{equation} \left\langle \frac{\sigma _{\rm side}}{\sigma _\lambda } \right\rangle = \frac{ \left[\sum f^2_\lambda - \bar{f}_\lambda ^2 \right]^{1/2}}{\sum \sigma _\lambda }, \end{equation} \tag{ 3 }$$

where the summations and average flux is evaluated in the quasar rest-frame 1460 Å  <  λ_rest  <  1510 Å.

In Figure 3, this quantity is averaged over the 500 individual quasars per redshift bin and plotted as a function of the observed wavelength corresponding to λ = (1 + 〈z_qso〉)1485 Å. With a perfect noise estimate, 〈σ_side/σ_λ〉 should be unity at all wavelengths, but we see that the BOSS pipeline underestimates the true noise in the spectra at λ  ≲  5000 Å, by up to ∼15% at the blue end of the spectra, with an overall tilt that changes over to an overestimate at λ ≳ 4500 Å. Lee et al. (2013) and Palanque-Delabrouille et al. (2013) provide a set of correction vectors that can be applied to the pipeline noise estimates to bring the latter to within several percent of the true noise level across the wavelength coverage of the blue spectrograph.

Zoom In Zoom Out Reset image size

Unfortunately, these noise corrections are inadequate for our purposes, since we want to generate realistic mock spectra that have different realizations of the Lyα forest transmission field from the actual spectra, i.e., a different $\mathcal {F}_{\lambda }$. We therefore require a method that not only accurately estimates the noise in a given BOSS spectrum, but also separates out the photon-counting and CCD terms in the variance that results from applying the Horne (1986) optimal spectral extraction algorithm:

$$\begin{equation} \sigma _{\lambda }^2 = S_{\lambda }\left(\mathcal {F}_{\lambda } + s_{\lambda }\right) + S_{\lambda }^2\sigma ^2_{\rm RN}, \end{equation} \tag{ 4 }$$

where σ_RN is the CCD read-noise.

To resolve this issue, we apply our own novel statistical method to the individual BOSS exposures to generate co-added spectra while simultaneously estimating the corresponding noise parameters for each individual spectrum. This procedure, which uses a Gibbs-sampled Markov-Chain Monte Carlo (MCMC) algorithm, is described in detail in the Appendix. Initially, we attempted to model the noise with just a single constant noise parameter that rescales the read-noise term of Equation (4), but this was found to be inadequate. This is likely because an optimal extraction algorithm weights by the product of the S/N and object profile, causing the corresponding variance to have a non-linear dependence on the flux and sky level. Furthermore, systematic errors in the reduction, sky-subtraction and calibration will result in additional noise contributions, which could depend on sky level, object flux, or wavelength, hence deviating from this simple model.

After considerable trial-and-error to find a model that best minimizes the bias illustrated in Figure 3, we settled on the form:

$$\begin{equation} \sigma _{\lambda i}^2 = A_1 \hat{S}_{\lambda i}\left(\mathcal {F}_{\lambda }+ s_{\lambda _i}\right) + A_2 \hat{S}_{\lambda i}^2\sigma ^2_{\rm RN, eff} \sigma _{\rm disp}(\lambda), \end{equation} \tag{ 5 }$$

where

$$\begin{equation} \hat{S}_{\lambda i}=S_{\lambda i} \left(1 - \exp (-A_3 \lambda + A_4) \right), \end{equation} \tag{ 6 }$$

where the A_j are free parameters in our noise model, while the σ_disp(λ) factor in the second term (the pixel dispersion) provides a rough approximation for the wavelength dependence of the spot size (i.e., the size of the raw CCD image in the spatial direction). Meanwhile, σ_disp = 12 is the average CCD read-noise per wavelength bin in the BOSS spectra (D. J. Schlegel et al., in preparation). The quantities s_{λ, i}, S_{λ, i}, and σ_disp(λ) (sky flux, calibration vector, and dispersion, respectively) are taken directly from the BOSS pipeline.

In addition, we assume that the pixel noise can be modeled as a Gaussian distribution with a variance given by Equation (5). The first, photon counting, term in the equation should formally be modeled as a Poisson distribution, but since the BOSS spectrograph always receives ≳ 30–40 counts even at the blue end of the spectrograph where the counts are the lowest, it is reasonable to use the Gaussian approximation because even in the limit of low S/N (i.e., when the spectrum is dominated by the sky flux), the moderate resolution ensures that there are at least several dozen sky photons per pixel in each exposure.

For each BOSS spectrum, we use the MCMC procedure described in the Appendix to combine the multiple exposures while simultaneously estimating the noise parameters A_j and true observed spectrum, $\mathcal {F}_{\lambda }$. With the optimal estimates of A_j and $\mathcal {F}_{\lambda }$ for a given spectrum, the estimated noise variance is then simply Equation (5).

An important advantage of the form in Equation (5) is that the object photon noise $\propto \mathcal {F}_{\lambda }$ is explicitly separated out. This facilitates the construction of a mock spectrum with the same noise characteristics as a true spectrum, but with a different spectral flux. For example, a mock spectrum of the Lyα forest will have a very different transmission field than the original data, and so the variance due to object photon counting noise can be added appropriately, in addition to contributions from the known sky, and the read noise term (Equation (5)). Our empirical determination of the parameters governing this noise model for each individual spectrum form a crucial ingredient in our forward model, which we will describe in Section 4.

Our MCMC procedure works for spectra from a single camera, either red or blue; we have not yet generalized it to combine blue and red spectra of each object. However, the spectral range of the blue camera alone (≈3600–6400 Å) covers the Lyα forest up to z  ∼  5, i.e., most practical redshifts for Lyα forest analysis. For the purposes of this paper, we restrict ourselves to spectra from the blue camera alone.

In Figures 4 and 5, we show examples of co-added BOSS quasar spectra, using both the MCMC procedure and the standard BOSS pipeline. In the upper panels, the MCMC co-adds are not noticeably different from the BOSS pipeline, although the numerical values are different. In the lower panels, we show the estimated noise from both methods—the differences are larger than in the fluxes but still difficult to distinguish by eye.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We therefore return to the statistical analysis by calculating 〈σ_side/σ_λ〉, the ratio of the pixel variance against the estimated noise from the flat 1460 Å  <  λ_rest  <  1510 Å region of BOSS quasars; this ratio, computed for our MCMC coadds, is plotted in Figure 3. With these new co-adds, we see that this ratio is within roughly ±3% of unity across the entire λ  ∼  3800–5000 Å wavelength range relevant to our subsequent analysis, with an overall bias of 1% (i.e., the noise is still underestimated by this level). Crucially, we have removed the strong wavelength dependence of 〈σ_side/σ_λ〉 that was present in the standard pipeline, and we suspect most of the scatter about unity is caused by the limited number of quasars (500 per bin) available for this estimate, which will be mitigated by the larger number of quasars spectra available in the subsequent BOSS data releases. In principle, we could correct the remaining 1% noise bias, but since our selected spectra have S/N > 6, this remaining 1% noise bias would smooth the forest transmission PDF by an amount roughly 1/25 of the average PDF bin width (ΔF = 0.05). As we shall see, there are other systematic uncertainties in our modeling that have much larger effects than this, therefore we regard our noise estimates as adequate for the subsequent transmission PDF analysis, without requiring any further correction.

In order to obtain the transmitted flux F of the Lyα forest²⁰ we first need to divide the observed flux, $\mathcal {F}_{\lambda }$, by an estimate for the quasar continuum, c. We use the version of MF-PCA continuum fitting (Lee et al. 2012) described in Lee et al. (2013). Initially, PCA fitting with eight eigenvectors is performed on each quasar spectrum redward of the Lyα line (λ_rest = 1216–1600 Å) in order to obtain a prediction for the continuum shape in the λ_rest  <  1216 Å Lyα forest region (e.g., Suzuki et al. 2005). The slope and amplitude of this initial continuum estimate is then corrected to agree with the Lyα forest mean transmission, 〈F〉_cont(z), at the corresponding absorber redshifts, using a linear correction function.

The only difference in our continuum-fitting with that in Lee et al. (2013) is that here we use the latest mean flux measurements of Becker et al. (2013) to constrain our continua. Their final result yielded the power-law redshift evolution of the effective optical depth in the unshielded Lyα forest, defined in their paper as $N_{{\rm H\,\scriptsize{I}}}\le 10^{17.2}\ \mathrm{cm^{-2}}$ (although they only removed contributions from $N_{{\rm H\,\scriptsize{I}}}\ge 10^{19}\ \mathrm{cm^{-2}}$ absorbers). This is given by

$$\begin{equation} \tau _{\rm Ly\alpha,B13}(z) \equiv {-}\ln (\langle F \rangle (z)) = \tau _0\left(\frac{1+z}{1+z_0}\right)^{\beta } + C, \end{equation} \tag{ 7 }$$

with best-fit values of [τ₀, β, C] = [0.751, 2.90, −0.132] at z₀ = 3.5.

However, the actual raw measurement made by Becker et al. (2013) is the effective total absorption within the Lyα forest region of their quasars, which also contain contributions from metals and optically thick systems:

$$\begin{equation} \tau _{\rm eff}(z) \equiv \tau _{\rm Ly\alpha,B13}(z) + \tau _{\rm {metals}} + \tau _{\rm {LLS}}(z), \end{equation} \tag{ 8 }$$

where τ_metals and τ_LLS(z) denote the IGM optical depth contributions from metals and LLSs, respectively. For the purposes of our continuum-fitting, the quantity we require is τ_eff(z), since the τ_metals and τ_LLS(z) contributions are also present in our BOSS spectra. Becker et al. (2013) did not publish their raw τ_eff(z), therefore we must now "uncorrect" the metal and LLS contributions from the published τ_{Lyα, B13}(z). The discussion below therefore attempts to retrace their footsteps and does not necessarily reflect our own beliefs regarding the actual level of these contributions.

We find τ_metals = 0.02525 by simply averaging over the Schaye et al. (2003) metal correction tabulated by Faucher-Giguère et al. (2008) (i.e., the 2.2 ⩽ z  ⩽  2.5 values in Δz = 0.1 bins from their Table 4), that were used by Becker et al. (2013) to normalize their relative mean flux measurements. Note that there is no redshift dependence on τ_metals in this context, because Becker et al. (2013) argued that the metal contribution does not vary significantly over their redshift range. Whether or not this is really true is unimportant to us at the moment, since we are merely "uncorrecting" their measurement.

The LLS contribution to the optical depth is re-introduced by integrating over $f(N_{{\rm H\,\scriptsize{I}}}, b, z)$, the column-density distribution of neutral hydrogen absorbers:

$$\begin{eqnarray} \tau _{\rm LLS}(z) &\approx & \frac{1+z}{\lambda _{\rm Ly\alpha }} \int ^{N_{{\rm max}}}_{N_{{\rm min}}} {dN}_{{\rm H\,\scriptsize{I}}}\int {d} b \nonumber \\ & & \times f(N_{{\rm H\,\scriptsize{I}}}, b, z) W_0(N_{{\rm H\,\scriptsize{I}}},b), \end{eqnarray} \tag{ 9 }$$

where b is the Doppler parameter and $W_0(N_{{\rm H\,\scriptsize{I}}},b)$ is the rest-frame equivalent width (EW; we use the analytic approximation given by Draine 2011, valid in the saturated regime).

Following Becker et al. (2013), we adopted a fixed value of b = 20 km s⁻¹ and assumed that $f(N_{{\rm H\,\scriptsize{I}}},z) = f(N_{{\rm H\,\scriptsize{I}}}) {d}n/{d}z$, where $f(N_{{\rm H\,\scriptsize{I}}})$ is given by the z = 3.7 broken power-law column density distribution of Prochaska et al. (2010) and dn/dz∝(1 + z)². Becker et al. (2013) had corrected for super-LLSs and DLAs in the column-density range [N_min, N_max] = [10¹⁹ cm⁻², 10²² cm⁻²], but as discussed above we have discarded all sightlines that include $N_{{\rm H\,\scriptsize{I}}}\ge 10^{20.3} \,\mathrm{cm^{-2}}$ DLAs; therefore, we reintroduce the optical depth contribution for super-LLSs, i.e., [N_min, N_max] = [10¹⁹ cm⁻², 10^20.3 cm⁻²]. We find τ_LLS(z) = 0.0022 × [(1 + z)/3]³. This is a small correction, giving rise to only a 0.5% change in 〈F〉 at z = 3.0.

This estimate of the raw absorption, 〈F〉_eff(z) = exp [ − τ_eff(z)], is now the constraint used to fit the continua of the BOSS quasars, i.e., we set 〈F〉_cont = 〈F〉_eff(z). Note that in our subsequent modeling of the data, we will use the same 〈F〉_cont(z) to fit the mock spectra to ensure an equal treatment between data and mocks. Since 〈F〉_cont(z) includes a contribution from $N_{{\rm H\,\scriptsize{I}}}\,{<}\,10^{20.3}\ \mathrm{cm^{-2}}$ optically thick systems, our mock spectra will need to account for these systems as we shall describe in Section 4.2.

The MF-PCA technique requires spectral coverage in the quasar rest-frame interval 1000–1600 Å. However, as noted in the previous section, we work with co-added BOSS spectra from only the blue cameras covering λ  ≲  6400 Å; this covers the full 1000–1600 Å interval required for the PCA fitting only for z  ≲  3 quasars. However, the differences in the fluxes between our MCMC co-adds and the BOSS pipeline co-adds are relatively small, and we do not expect the relative shape of the quasar spectrum to vary significantly. We can thus carry out PCA fitting on the BOSS pipeline co-adds, which cover the full observed range (3700–10000 Å), to predict the overall quasar continuum shape. This initial prediction is then used to perform mean flux regulation using the MCMC co-adds and noise estimates, to fine-tune the amplitude of the continuum fits.

The observed flux, f_λ, is divided by the continuum estimate, c, to derive the Lyα forest transmission, F = f_λ/c. For each quasar, we define the Lyα forest as the rest wavelength interval 1041–1185 Å. This wavelength range conservatively avoids the quasar's Lyβ/O vi emission line blend by Δv  ∼  3000 km s⁻¹ on the blue end, as well as the proximity zone close to the quasar redshift by staying Δv  ∼  10, 000 km s⁻¹ from the nominal quasar systemic redshift. We are now in a position to measure the transmission PDF, which is simply the histogram of pixel transmissions F ≡ exp (− τ).

Since the Lyα forest evolves as a function of redshift, we measure the BOSS Lyα forest transmission PDF in three bins with mean redshifts of 〈z〉 = 2.3, 〈z〉 = 2.6, and 〈z〉 = 3.0, and bin sizes of Δz = 0.3. These redshift bins were chosen to match the simulations outputs (Section 4.1) that we will later use to make mock spectra to compare with the observed PDF; this choice of binning leads to the gap at 2.75  <  z  <  2.85 as seen in Figure 2. In this paper, we restrict ourselves to z  ≲  3 since the primary purpose is to develop the machinery to model the BOSS spectra. In subsequent papers, we will apply these techniques to analyze the transmission PDF in the full 2  ≲  z  ≲  4 range using the larger samples of subsequent BOSS data releases (DR10; Ahn et al. 2014).

Another consideration is that the transmission PDF is strongly affected by the noise in the data. While we will model this effect in detail (Section 4), there is a large distribution of S/N within our subsample ranging from S/N = 6 pixel⁻¹ to S/N  ∼  20 pixel⁻¹. We therefore further divide the sample into three bins depending on the median S/N per pixel within the Lyα forest: 6  <  S/N  <  8, 8  <  S/N  <  10, S/N > 10. The consistency of our results across the S/N bins will act as an important check for the robustness of our noise model (Section 3.1).

We now have nine redshift and S/N bins in which we evaluate the transmission PDF from BOSS; the sample sizes are summarized in Table 1. For each bin, we have selected quasars that have at least 30 Lyα forest pixels within the required redshift range, and which occupy the quasar rest-frame interval 1041–1185 Å. The co-added spectrum is divided with its MF-PCA continuum estimate (described in the previous section) to obtain the transmitted flux, F, in the desired pixels. We then compute the transmission PDF from these pixels.

Physically, the possible values of the Lyα forest transmission range from F = 0 (full absorption) to F = 1 (no absorption). However, the noise in the BOSS Lyα forest pixels, as well as continuum fitting errors, leads to pixels with F  <  0 and F > 1. We therefore measure the transmission PDF in the range −0.2  <  F  <  1.5, in 35 bins with width Δ(F) = 0.05, and normalized such that the area under the curve is unity. The statistical errors on the transmission PDF are estimated by the following method: we concatenate all the individual Lyα forest segments that contribute to each PDF, and then carry out bootstrap resampling over Δv = 2 × 10⁴ km s⁻¹ segments with 200 iterations. This choice of Δv corresponds to ∼250–300 Å in the observed frame at z  ∼  2–3—according to Rollinde et al. (2013), this choice of Δv and number of iterations should be sufficient for the errors to converge (see also Appendix B in McDonald et al. 2000).

In Figure 6, we show the Lyα forest transmission PDF measured from the various redshift and S/N subsamples in our BOSS sample. At fixed redshift, the PDFs from the lower S/N data have a broader shape as expected from increased noise variance. With increasing redshift, there are more absorbed pixels, causing the transmission PDFs to shift toward lower F values. As discussed previously, there is a significant portion of F  >  1 pixels due to a combination of pixel noise and continuum errors, with a greater proportion of F  >  1 pixels in the lower-S/N subsamples as expected. Unlike the high-resolution transmission PDF, at 〈z〉  ≲  3 there are few pixels that reach F = 0. This effect is due to the resolution of the BOSS spectrograph, which smooths over the observed Lyα forest such that even saturated Lyα forest absorbers with $N_{{\rm H\,\scriptsize{I}}}\,{\gtrsim}\,10^{14}\hbox{--}10^{16}\, \mathrm{cm^{-2}}$ rarely reach transmission values of F  ≲  0.3. The pixels with F  ≲  0.3 are usually contributed either by blends of absorbers or optically thick LLSs (see also Pieri et al. 2014).

Zoom In Zoom Out Reset image size

An advantage of our large sample size is that it is also able to directly estimate the error covariances, C_boot, via bootstrap resampling—an example is shown in Figure 7. In contrast to the Lyα forest transmission PDF from high-resolution data which have significant off-diagonal covariances (Bolton et al. 2008), the error covariance from the BOSS transmission PDF is nearly diagonal with just some small correlations between neighboring bins, although we also see some anti-correlation between transmission bins at F  ∼  0.8 and F  ∼  1.

Zoom In Zoom Out Reset image size

It is interesting to compare the transmission PDF from our data with that measured by Desjacques et al. (2007) from SDSS DR3. This comparison is shown in Figure 8, in which the transmission PDFs calculated from SDSS DR3 Lyα forest spectra with S/N > 4 (kindly provided by Dr. V. Desjacques) are shown for two redshift bins, juxtaposed with the BOSS transmission PDFs calculated from spectra with the same redshift and S/N cuts.

Zoom In Zoom Out Reset image size

While there is some resemblance between the two PDFs, the most immediate difference is that the Desjacques et al. (2007) PDFs are shifted to lower transmission values, i.e., the mean transmission, 〈F〉, is considerably smaller than that from our BOSS data: 〈F〉(z = 2.4) = 0.73 and 〈F〉(z = 3.0) = 0.64 from their measurement, whereas the BOSS PDFs have 〈F〉(z = 2.4) = 0.80 and 〈F〉(z = 3.0) = 0.70. This difference arises because the Desjacques et al. (2007) used a power-law continuum (albeit with corrections for the weak emission lines in the quasar continuum) extrapolated from λ_rest > 1216 Å in the quasar rest-frame; this does not take into account the power-law break that appears to occur in low-redshift quasar spectra at λ_rest ≈ 1200 Å (Telfer et al. 2002; Suzuki 2006). Later in their paper, Desjacques et al. (2007) indeed conclude that this must be the case in order to be consistent with other 〈F〉(z) measurements. Our continua, in contrast, have been constrained to match existing measurements of 〈F〉(z), for which there is good agreement between different authors at z  ≲  3 (e.g., Faucher-Giguère et al. 2008; Becker et al. 2013).

Another point of interest in Figure 8 is that the error bars of the BOSS sample are considerably smaller than those of the earlier measurement. This difference is largely due to the significantly larger sample size of BOSS. The proportion of pixels with F  ≲  0 appears to be smaller in the BOSS PDFs compared with the older data set, but this is because Desjacques et al. (2007) did not remove DLAs from their data.

We next describe the creation of mock Lyα absorption spectra designed to match the properties of the BOSS data.

As the basis for our mock spectra, we use hydrodynamic simulations run with a modification of the publicly available GADGET-2 code. This code implements a simplified star formation criterion (Springel et al. 2005) that converts all gas particles that have an overdensity above 1000 and a temperature below 10⁵ K into star particles (see Viel et al. 2004). The simulations used are described in detail in Becker et al. (2011) and in Viel et al. (2013a).

The reference model that we use is a box of length 20 h⁻¹ comoving Mpc with 2 × 512³ gas and cold DM particles (with a gravitational softening length of 1.3 h⁻¹ kpc) in a flat ΛCDM universe with cosmological parameters Ω_m = 0.274, Ω_b = 0.0457, n_s = 0.968, H₀ = 70.2 km s⁻¹ Mpc⁻¹ and σ₈ = 0.816, in agreement both with nine-year WMAP (Komatsu et al. 2011) and Planck data (Planck Collaboration et al. 2014). The initial condition power spectra are generated with CAMB (Lewis et al. 2000). For the boxes considered in this work, we have verified that the transmission PDF has converged in terms of box size and resolution.

We explore the impact of different thermal histories on the Lyα forest by modifying the ultraviolet (UV) background photo-heating rates in the simulations as done in, e.g., Bolton et al. (2008). A power-law temperature–density relation, T = T₀Δ^{γ − 1}, arises in the low-density IGM (Δ < 10) as a natural consequence of the interplay between photo-heating and adiabatic cooling (Hui et al. 1997; Gnedin & Hui 1998). The value of γ within a simulation can be modified by varying a density-dependent heating term (see, e.g., Bolton et al. 2008). We consider a range of values for the temperature at mean density, T₀, and the power-law index of the temperature–density relation, γ, based on the observational measurements presented recently by Becker et al. (2011). These consist of a set of three different indices for the temperature–density relation, γ(z = 2.5)  ∼  1.0, 1.3, 1.6, that are kept roughly constant over the redshift range z = [2–6] and three different temperatures at mean density, T₀(z = 2.5)  ∼  [11000, 16000, 21500] K, which evolve with redshift, yielding a total of nine different thermal histories. Between z = 2 and z = 3 there is some temperature evolution and the IGM becomes hotter at low redshift; at z = 2.3, the models have T₀  ∼  [13000, 18000, 23000] K. We refer to the intermediate temperature model as our "reference" model, or T_REF, while the hot and cold models are referred to as T_HOT and T_COLD, respectively. The values of T₀ of our simulations at the various redshifts are summarized in Table 2.

Approximately 4000 core hours were required for each simulation run to reach z = 2. The physical properties of the Lyα forest obtained from the Tree PM/SPH code GADGET-2 are in agreement at the percent level with those inferred from the moving-mesh code AREPO (Bird et al. 2013) and with the Eulerian code ENZO (O'Shea et al. 2004).

For this study, the simulation outputs were saved at z = [2.3, 2.6, 3.0], from which we extract 5000 optical depth sightlines binned to 2048 pixels each. To convert these to transmission spectra, the optical depths were rescaled such that the skewers collectively yielded a desired mean transmission, 〈F〉_Lyα ≡ exp (− τ_Lyα). For our fiducial models, we would like to use the mean transmission values estimated by Becker et al. (2013), which we denote as 〈F〉_{Lyα, B13} ≡ exp (− τ_{Lyα, B13}). However, their estimates assume certain corrections from optically thick systems and metal absorption. We therefore add back in the corrections they made (see discussion in Section 3.2) to get their "raw" measurement for 〈F〉 that now includes all optically thick systems and metals, and then remove these contributions assuming our own LLS and metal absorption models (see below).

Later in the paper, we will argue that our PDF analysis in fact places independent constraints on 〈F〉_Lyα.

In principle, all optically thick Lyα absorbers such as LLSs and DLAs should be discarded from Lyα forest analyses, since they do not trace the underlying matter density field in the same way as the optically thin forest (Equation (1)), and require radiative transfer simulations to accurately capture their properties (e.g., McQuinn et al. 2011; Rahmati et al. 2013).

While DLAs are straightforward to identify through their saturated absorption and broad damping wings even in noisy BOSS data (see, e.g., Noterdaeme et al. 2012), the detection completeness of optically thick systems through their Lyα absorption drops rapidly at $N_{{\rm H\,\scriptsize{I}}}\,{\lesssim}\,10^{20} \ \mathrm{cm^{-2}}$. Even in high-S/N, high-resolution spectra, optically thick systems can only be reliably detected through their Lyα absorption at $N_{{\rm H\,\scriptsize{I}}}\gtrsim 10^{19}\ \mathrm{cm^{-2}}$ ("super-LLS"). Below these column densities, optically thick systems can be identified either through their rest-frame 912 Å Lyman-limit (albeit only one per spectrum) or using higher-order Lyman-series lines (e.g., Rudie et al. 2013). Neither of these approaches have been applied in previous Lyα forest transmission PDF analyses (McDonald et al. 2000; Kim et al. 2007; Calura et al. 2012; Rollinde et al. 2013), so arguably all these analyses are contaminated by LLSs.

Instead of attempting to remove LLSs from our observed spectra, we incorporate them into our mock spectra through the following procedure. For each PDF bin, we evaluate the total redshift pathlength of the contributing BOSS spectra (and corresponding mocks)—this quantity is summarized in Table 1. This is multiplied by l_LLS(z), the number of LLS per unit redshift, to give the total number of LLS expected within our sample. We used the published estimates of this quantity by Ribaudo et al. (2011)²¹ which is valid over 0.24  <  z  <  4.9:

$$\begin{equation} l_{\rm LLS}(z) = l_{z0} (1 + z)^{\gamma _{\rm LLS}}, \end{equation} \tag{ 10 }$$

where l_z0 = 0.1157 and γ_LLS = 1.83.

After estimating the total number of LLSs in our mock spectra, l_LLS(z)Δz, we add them at random points within our set of simulated optical depth skewers. We also experimented with adding LLSs such that they are correlated with regions that already have high column density (e.g., Font-Ribera & Miralda-Escudé 2012), but we found little significant changes to the transmission PDF and therefore stick to the less computationally intensive random LLSs.

For each model LLS, we then draw a column density using the published LLS column density distribution, $f(N_{{\rm H\,\scriptsize{I}}})$, from Prochaska et al. (2010). This distribution is measured at z ≈ 3.7, so we make the assumption that $f(N_{{\rm H\,\scriptsize{I}}})$ does not evolve with redshift between 2  ≲  z  ≲  3.7. For our column densities of interest, this distribution is represented by the broken power laws:

$$\begin{equation} f(N_{{\rm H\,\scriptsize{I}}}) =\left\lbrace \begin{array}{@{}l@{\quad }l@{}}k_1 N_{{\rm H\,\scriptsize{I}}}^{-0.8} & \mathrm{if}\; 10^{17.5}\,{<}\,N_{{\rm H\,\scriptsize{I}}}\,{<}\,10^{19.0} \\ k_2 N_{{\rm H\,\scriptsize{I}}}^{-1.2} & \mathrm{if}\; 10^{19.0}\,{<}\,N_{{\rm H\,\scriptsize{I}}}\,{<}\,10^{20.3}\end{array}\right.. \end{equation} \tag{ 11 }$$

For the normalizations k₁ and k₂, we demand that

$$\begin{equation} \int ^{10^{19.0}}_{10^{17.5}} k_1 N_{{\rm H\,\scriptsize{I}}}^{-0.8}\, {d}N_{{\rm H\,\scriptsize{I}}}+ \int ^{10^{20.3}}_{10^{19.0}} k_2 N_{{\rm H\,\scriptsize{I}}}^{-1.2}\, {d}N_{{\rm H\,\scriptsize{I}}}= 1, \end{equation} \tag{ 12 }$$

and require both power laws to be continuous at $N_{{\rm H\,\scriptsize{I}}}= 10^{19.0}\ \mathrm{cm^{-2}}$. These constraints produce k₁ = 10^−4.505 and k₂ = 10^3.095. After drawing a random value for the column density of each LLS, we add the corresponding Voigt profile to the optical depth in the simulated skewer.

In addition to the LLS with column densities of $10^{17.5}\ \mathrm{cm^{-2}}\,{<}\,N_{{\rm H\,\scriptsize{I}}}\,{<}\,10^{20.3}\ \mathrm{cm^{-2}}$ that are defined to have $\tau _{{\rm H\,\scriptsize{I}}} \ge 2$, there is also a population of partial Lyman-limit systems (pLLSs) that are not well-captured in our hydrodynamical simulations since they have column densities ($10^{16.5}\ \mathrm{cm^{-2}}\,{\lesssim}\,N_{{\rm H\,\scriptsize{I}}}< 10^{17.5}\ \mathrm{cm^{-2}}$) at which radiative transfer effects become significant ($\tau _{{\rm H\,\scriptsize{I}}} \gtrsim 0.1$). However, the incidence rates and column-density distribution of pLLSs are ill-constrained since they are difficult to detect in normal LLS searches. We therefore account for the pLLS by extrapolating the low end of the power-law distribution in Equation (11) down to $N_{{\rm H\,\scriptsize{I}}}= 10^{16.5}\ \mathrm{cm^{-2}}$, i.e.,

$$\begin{equation} f(10^{16.5}\ \mathrm{cm^{-2}}\,{<}\,N_{{\rm H\,\scriptsize{I}}}\,{<}\,10^{17.5}\ \mathrm{cm^{-2}}) = k_1 N_{{\rm H\,\scriptsize{I}}}^{-0.8}. \end{equation} \tag{ 13 }$$

This simple extrapolation does not take into account constraints from the mean free path of ionizing photons (e.g., Prochaska et al. 2010), which predicts a steeper slope for the pLLS distribution, but we will explore this later in Section 5.2.

Comparing the integral of this extrapolated pLLS distribution with Equation (12) leads us to conclude that

$$\begin{equation} l_{\rm pLLS}(z) = 0.197\; l_{\rm LLS}(z), \end{equation} \tag{ 14 }$$

and we proceed to randomly add pLLSs to our mock spectra in the same way as LLSs.

The other free parameter in our LLS model is their effective b-parameter distribution. However, due to the observational difficulty in identifying $N_{{\rm H\,\scriptsize{I}}}\,{\lesssim}\,18.5\ \mathrm{cm^{-2}}$ LLSs, the b-parameter distribution of this distribution has, to our knowledge, never been quantified. Due to this lack of knowledge, it is common to simply adopt a single b-value when attempting to model LLSs (e.g., Font-Ribera & Miralda-Escudé 2012; Becker et al. 2013). We therefore assume that all our pLLSs and LLSs have a b-parameter of b = 70 km s⁻¹ similar to DLAs (Prochaska & Wolfe 1997), an "effective" value meant to capture the blending of multiple Lyα components. However, the b-parameter for this population of absorbers is a highly uncertain quantity and as we shall see, it will need to be modified to provide a satisfactory fit to the data although it will turn out to not strongly affect our conclusions regarding the IGM temperature–density relationship. Figure 10 illustrates the effect of adding a LLS to a mock spectrum, and its subsequent smoothing with the spectrograph resolution kernel (next section).

Zoom In Zoom Out Reset image size

The spectral resolution of SDSS/BOSS spectra is R ≡ λ/Δλ ≈ 1500–2500 (Smee et al. 2013). The exact value varies significantly both as a function of wavelength, and across different fibers and plates depending on observing conditions (Figure 1).

For each spectrum, the BOSS pipeline provides an estimate of the 1σ wavelength dispersion at each pixel, σ_disp, in units of the co-added wavelength grid size (Δlog₁₀λ = 10⁻⁴). The spectral resolution at that pixel can then be obtained from the dispersion, through the following conversion: R ≈ (2.35 × 1 × 10⁻⁴ln 10 σ_disp)⁻¹. Figure 1 shows the pixel dispersions from 236 randomly selected BOSS quasar as a function of wavelength at the blue end of the spectrograph. Even at fixed wavelength, there is a considerable spread in the dispersion, e.g., ranging from σ_disp ≈ 0.9–1.8 at 3700 Å. The value of σ_disp typically decreases with wavelength (i.e., the resolution increases).

In their analysis of the Lyα forest 1D transmission power spectrum, Palanque-Delabrouille et al. (2013) made their own study of the BOSS spectral resolution by directly analyzing the line profiles of the mercury and cadmium arc lamps used in the wavelength calibration. They found that the pipeline underestimates the spectral resolution as a function of fiber position (i.e., CCD row) and wavelength: the discrepancy is <1% at blue wavelengths and near the CCD edges, but increases to as much as 10% at λ  ∼  6000 Å near the center of the blue CCD (compare with Figure 4 in Palanque-Delabrouille et al. 2013). Our analysis is limited to λ  ⩽  5045 Å, i.e., z  ⩽  3.15, where the discrepancy is under 4%. Nevertheless, we implement these corrections to the BOSS resolution estimate to ensure that we model the spectral resolution to an accuracy of <1%.

For each BOSS Lyα forest segment that contributes to the observed transmission PDFs discussed in Section 3, we concatenate randomly selected transmission skewers from the simulations described in the previous section. This is because the simulation box size of L = 20 h⁻¹ Mpc (Δv  ∼  2000 km s⁻¹) is significantly shorter than the path length of our redshift bins (Δz = 0.3, or Δv ≈ 27, 000 km s⁻¹). This ensures that each BOSS spectrum in our sample has a mock spectrum that is exactly matched in pathlength.

We then directly convolve the simulated skewers by a Gaussian kernel with a standard deviation that varies with wavelength, using the estimated resolution from the real spectrum, multiplied by the Palanque-Delabrouille et al. (2013) resolution corrections. The effect of smoothing on the transmission PDF is illustrated by the dashed red curve in Figure 9(b). Smoothing decreases the proportion of pixels with high transmission (F ≈ 1) and with high absorption (F ≈ 0), and increases the number of pixels with intermediate transmission values.

Metal absorption along our observed Lyα forest sightlines acts as a contaminant since their presence alters the observed statistics of the Lyα forest. In high-resolution data, this contamination is usually treated by directly identifying and masking the metal absorbers, although in the presence of line blending it is unclear how thorough this approach can be.

With the lower S/N and moderate resolution of the BOSS data, direct metal identification and masking is not a viable approach. Furthermore, most of the weak metal absorbers seen in high-resolution spectra are not resolved in the BOSS data.

Rather than removing metals from the BOSS Lyα forest spectra, we instead add metals as observed in lower-redshift quasar spectra. In other words, we add absorbers observed in the rest-frame λ_rest ≈ 1260–1390 Å region of lower-redshift quasars with 1 + z_qso ≈ (1216 Å/1300 Å)(1 + 〈z〉), such that the observed wavelengths are matched to the Lyα forest segment with average redshift 〈z〉. Figure 11 is an illustration that illustrates this concept. This method makes no assumption about the nature of the metal absorption in the Lyα forest, and includes all resolved metal absorption spanning the whole range of redshifts down to z  ∼  0. The disadvantage of this method is that it does not include metals with intrinsic wavelengths λ  ≲  1300 Å, but the relative contribution of such metal species toward the transmission PDF should be small²² since most of the metal contamination comes from low-redshift (z  ≲  2) C iv and Mg ii.

Zoom In Zoom Out Reset image size

We use a metal catalog generated by B. Lundgren et al. (in preparation; see also Lundgren et al. 2009), which lists absorbers in SDSS (Schneider et al. 2010) and BOSS quasar spectra (Pâris et al. 2012)—the SDSS spectra were included in order to increase the number of z_qso ≈ 1.9–2.0 quasars needed to introduce metals into the 〈z〉 = 2.3 Lyα forest mock spectra, which are not well sampled by the BOSS target selection (Ross et al. 2012). We emphasize that we work with the "raw" absorber catalog, i.e., the individual absorption lines have not been identified in terms of metal species or redshift. For each quasar, the catalog provides a line list with the observed wavelength, EW (W_r), FWHM, and detection S/N, $W_r/\sigma _{W_r}$. To ensure a clean catalog, we use only $W_r/\sigma _{W_r} \ge 3.5$ absorbers in the catalog that were identified from quasar spectra with S/N > 15 per angstrom redward of Lyα. The latter criterion ensures that even relatively weak lines (with EW ≳ 0.5 Å) are accounted for in our catalog. Figure 12 shows an example of the lower-redshift quasar spectra that we use for the metal modeling.

Zoom In Zoom Out Reset image size

However, we want to add a smooth model of the metal-line absorption to add to our mock spectra, rather than adding in a noisy spectrum. We therefore use a simple model as follows: For each Lyα forest segment we wish to model at redshift 〈z〉, we select an absorber line-list from a random quasar with 1 + z_qso ≈ (1216 Å/1300 Å)(1 + 〈z〉). We next assume that all resolved metals in the SDSS/BOSS spectra are saturated and thus in the flat regime of the curve-of-growth. The EW is then given by

$$\begin{equation} W_r \approx \left(\frac{2b}{c} \right) \sqrt{\ln (\tau _0/ \ln 2)}, \end{equation} \tag{ 15 }$$

where τ₀ is the optical depth at line center, b is the velocity width and c is the speed of light. In the saturated regime, W_r is mostly sensitive to changes in b while being highly insensitive to changes in τ₀. We can thus adopt τ₀ as a global constant and solve for b, given the W_r of each listed absorber in the selected "sideband" quasar. We have found that τ₀ = 3 provides a good fit for most of the absorbers.

We then add the Gaussian profile into our simulated optical depth skewers:

$$\begin{equation} \tau = \tau _0 \exp \left[-\left(\frac{c}{b}\right) \left(\frac{\Delta \lambda }{\lambda }\right)^2 \right] \end{equation} \tag{ 16 }$$

centered at the same observed wavelength, λ, as the real absorber. The red curve in Figure 12 shows our model for the observed absorbers, using just the observed wavelength, λ, and EW, W_r, from the absorber catalog.

Our method for incorporating metals is somewhat crude since one should, in principle, first deconvolve the spectrograph resolution from the input absorbers, and then add the metal absorbers into the mock spectra prior to convolving with the BOSS spectral resolution. In contrast, we fit b-parameters to the absorber catalog without spectral deconvolution, and therefore these b-parameters can be thought of as combinations of the true absorber width, b_abs, and the spectral dispersion, σ_disp, i.e., $b^2\,{\sim}\,b^2_{\rm abs} + \sigma ^2_{\rm disp}$. While technically incorrect, this seems reasonable since the template quasar spectra and forest spectra that we are attempting to model both have approximately the same resolution, and in practical terms this ad hoc approach does seem to be able to reproduce the observed metals in the lower-redshift quasar spectra (Figure 12). The other possible criticism of our approach is that it does not incorporate weak metal absorbers, although we attempted to mitigate this by setting a very high S/N threshold on the template quasars for the metals. However, we have checked that such weak metals do not significantly change the forest PDF (and indeed metals in general do not seriously affect the PDF, compare with Figure 9(c)).

We also tried adding metals with similar redshifts to—and correlated with—forest absorbers (e.g., absorption by Si ii and Si iii) measured in Pieri et al. (2010, 2014) using a method described in the Appendix of Slosar et al. (2011). We found a negligible impact on the transmission PDF owing mainly to the fact that these correlated metals contribute only ∼0.3% to the overall flux decrement, so we neglect this contribution in our subsequent analysis.

It is non-trivial to introduce the correct noise to a simulated Lyα forest spectrum: given a noise estimate from the observed spectrum, one needs to first ensure that the mock spectrum has approximately the same flux normalization as the data. This is challenging, as the Lyα forest transmission at any given pixel, which ranges from 0 to 1, will vary considerably between the simulated spectrum and the real data.

The simplest method of adding noise to a mock spectrum is simply to introduce Gaussian deviates using the pipeline noise estimate for each spectrum—this was essentially the method used by Desjacques et al. (2007) and the BOSS mocks described in Font-Ribera et al. (2012). However, with the MCMC co-addition procedure described in Section 3.1, we are in a position to model the noise in a more robust and self-consistent fashion.

Recall that the MCMC procedure returns posterior probabilities for two quantities: the true underlying spectral flux density, $\mathcal {F}_{\lambda }$, and the four free parameters A_j, which parameterize the noise in each spectrum. This estimate of the A_j from each quasar spectrum allows us to accurately model the pixel noise using Equation (5).

The MF-PCA method (Section 3.2) produces an estimate of the quasar continuum, c, providing approximately the correct flux level at each point in the spectrum. We can now multiply c with the simulated Lyα forest transmission spectra, F, which had already been smoothed to the same dispersion as its real counterpart (the estimated quasar continuum is already at approximately the correct smoothing, since it was fitted to the observed spectrum).

This procedure produces a noiseless mock spectrum with the correct flux normalization and smoothing. We can now generate noisy spectra corresponding to a given BOSS quasar, using our MCMC noise estimation described in Section 3.1. First, we substitute our mock spectrum as $\mathcal {F}_{\lambda }$ into Equation (5), and then combine the A_j noise parameters (estimated through our MCMC procedure) as well as the calibration vectors S_{λ, i} and sky estimates s_{λ, i}. This lets us generate self-consistent noise vectors corresponding to each individual exposure that make up the mock quasar spectrum, σ_λi. The noise vectors are then used to draw random Gaussian deviates that are added to the mock spectrum, on a per-pixel basis, to create the mock spectral flux density, f_λi. Finally, we combine these individual mock exposures into the optimal spectral flux density for the mock spectrum, through the expression (see the Appendix):

$$\begin{equation} f_{{\rm opt},\lambda } \equiv \frac{1}{\sigma ^2_{{\rm opt},\lambda }} \sum _{i} \frac{f_{\lambda i}}{\sigma _{\lambda i}^2}, \end{equation} \tag{ 17 }$$

where

$$\begin{equation} \frac{1}{\sigma _{{\rm opt},\lambda }^2} \equiv \sum _i \frac{1}{\sigma _{\lambda i}^2}. \end{equation} \tag{ 18 }$$

Figure 9(c) illustrates the effect of adding pixel noise to the smoothed Lyα forest transmission PDF. As expected, this scatters a significant fraction of pixels to F  >  1, and also to F  <  0 to a smaller extent.

With the noisy mock spectrum in hand (see, e.g., bottom panel of Figure 13), we can self-consistently include the effect of continuum errors into our model transmission PDFs by simply carrying out our MF-PCA continuum-fitting procedure on the individual noisy mock spectra. Dividing out the mock spectra with the new continuum fits then incorporates an estimate of the continuum errors (estimated by Lee et al. 2012 to be at the ∼4%–5% rms level) into the evaluated model transmission PDF. This estimated error includes uncertainties stemming from the estimation of the quasar continuum shape due to pixel noise, as well as the random variance in the mean Lyα forest absorption in individual lines-of-sight.

Zoom In Zoom Out Reset image size

Note that regardless of the overall mean absorption in the mock spectra (i.e., inclusive of our models for metals, LLSs, and mean forest absorption—see Section 5.4), we always use 〈F〉_cont(z), the same input mean transmission derived from Becker et al. (2013; described in Section 3.2) to fit the continua in both the data and mock spectra. While the overall absorption in our fiducial model is consistent with that from Becker et al. (2013), as we shall see later, the shape of the transmission PDF retains information on the true underlying mean transmission even if fitted with a mean flux regulated continuum with a wrong input 〈F〉(z).

The effect of continuum errors on the transmission PDF is shown in Figure 9(e): like pixel noise, it degrades the peak of the PDF, but only near F  ∼  1.

For each of our nine hydrodynamical simulations (sampling three points each in T₀ and γ), we determine the transmission PDF from the Lyα forest mock spectra that include the effects described in the previous section, for the various redshift and S/N subsamples in which we had measured the PDF in BOSS (Section 3.3). In Figure 14, we show the transmission PDFs for all our redshift and S/N subsamples in BOSS, compared with the corresponding simulated transmission PDFs from the T_REF simulation with γ = [1.0, 1.3, 1.6]. Note that the error bars shown are the diagonal elements of the covariance matrix estimated through bootstrap resampling on the data.

Zoom In Zoom Out Reset image size

At first glance, the model transmission PDFs seem to be a reasonable match for the data, especially considering we have carried out purely forward modeling without fitting for any parameters. However, when comparing the "pull," (p_{data, i} − p_{model, i})/σ_{p, i}, between the data and model (bottom panels of Figure 14), we see significant discrepancies in part due to the extremely small bootstrap error bars. Nevertheless, it is gratifying to see that the shape of the residuals is relatively consistent across the different S/N subsamples at fixed redshift and γ, since this indicates that our spectral noise model is robust.

We proceed to quantify the differences between the simulated transmission PDFs, p_model, and observed transmission PDFs, p_data, with the χ² statistic:

$$\begin{equation} \chi ^2 = \sum _{ij} (p_{\mathrm{model},i}-p_{\mathrm{data},i})^T C^{-1}_{ij}(p_{\mathrm{model},j}-p_{\mathrm{data},j}), \end{equation} \tag{ 19 }$$

where we use the bootstrap error covariance matrix, C_boot. Note that we also include a bootstrap error term that accounts for the sample variance in the model transmission PDFs, since our pipeline for generating mock spectra is too computationally expensive to include sufficiently large amounts of skewers to fully beat down the sample variance in the models.²³

We limit our model comparison to the range −0.1  ⩽  F  ⩽  1.2, i.e., 27 transmission bins with bin width Δ(F) = 0.05. This range covers pixels that have been scattered to "unphysical" values of F < 0 or F > 1 due to pixel noise, as is expected from the low-S/N of our BOSS data, and also captures >99.8% of the pixels within each of our data subsets. In particular, it is important to retain the bins with F  >  1 because the F  ∼  1 transmission bins are highly sensitive to γ (Lee 2012) and therefore we want to fully sample that region of the PDF even if it will require careful modeling of pixel noise and continuum errors.

There are two constraints on all our transmission PDFs: the normalization convention

$$\begin{equation} \int p(F)\; {dF} = 1 \end{equation} \tag{ 20 }$$

and the imposition of the same mean transmission due to the mean flux regulated continuum-fitting

$$\begin{equation} \int F\, p(F)\; {dF} = \langle F \rangle _{\rm cont} \end{equation} \tag{ 21 }$$

such that all the mock spectra have the same absorption, 〈F〉_cont(z). This is because the mock spectra have been continuum-fitted (Section 4.6) in exactly the same way as the BOSS spectra, which assumes the same mean Lyα transmission inferred from the Becker et al. (2013) measurements (Section 3.2). The "true" optically thin mean transmission, 〈F〉_Lyα, imposed on the simulation skewers is in principle a different quantity from 〈F〉_cont, since the latter includes contribution from metal contamination and optically thick LLSs.

This leaves us with ν = 27 − 1 − 2 = 24 degrees of freedom (dof) in our χ² comparison. The χ² for all the models shown in Figure 14 are shown in the corresponding figure legends.

In this initial comparison, the χ² values for the models in Figure 14 are clearly unacceptable: we find χ² ≳ 200 for 24 dof in all cases. However, it is interesting to note that the γ = 1.6 or γ = 1.3 models are preferred at all redshifts and S/N cuts. Note that the S/N = 8–10 subsamples (middle column in Figure 14) tends to have a slightly better agreement between model and data compared to the other S/N cuts at the same redshift: this simply reflects the smaller quantity of data of the subsample (compare with Table 1) and hence larger bootstrap errors.

A closer inspection of the residuals in Figure 14 indicate that there are two major sources of discrepancy between the models and data: first, at the low-transmission end, we underproduce pixels at 0.1  ≲  F  ≲  0.4 while simultaneously overproducing F  ≲  0.1 pixels, especially at 〈z〉 = 2.3 and 〈z〉 = 2.6. This seems to affect all γ models equally. Pieri et al. (2014) found that at BOSS resolution, pixels with F  ≲  0.3 come predominantly from saturated Lyα absorption from LLS. We therefore investigate possible modifications to our LLS model in Section 5.2.

The other discrepancy in the model transmission PDFs manifests at the higher-transmission end in the 〈z〉 = 2.6 and 〈z〉 = 3 subsamples, where we see a sinusoidal shape in the residuals at F  >  0.6 that appears consistent across different S/N. This portion of the transmission PDF depends on both γ and, as we shall see, on the assumed mean transmission 〈F〉(z), which we shall discuss in more detail in Section 5.4.

Finally, our transmission PDF model includes various uncertainties in the modeling of metals, LLSs, and continuum-fitting which have not yet been taken into account. In Section 5.3, we will estimate the contribution of these uncertainties, by means of a Monte Carlo method, in our error covariances.

With the moderate spectral resolution of BOSS, there are few individual pixels in the optically thin Lyα forest that reach transmission values of F  ≲  0.4. Such low-transmission pixels are typically due to either the blending of multiple absorbers (see, e.g., Figure 2 in Pieri et al. 2014), or optically thick systems (see Figure 10 in this paper).

As we have seen in Figure 14, at low-transmission values the discrepancy between data and model has a distinct shape, which is particularly clear at 〈z〉 = 2.3: the models underproduce pixels at 0.1  ≲  F  ≲  0.4 while at the same time overproducing saturated pixels with F ≈ 0.

To resolve this particular discrepancy would therefore require either drastically increasing the amount of clustering in the Lyα forest, or modifying our assumptions on the LLSs in our mock spectra. The first possibility seems rather unlikely since the Lyα forest power on relevant scales are well-constrained (Palanque-Delabrouille et al. 2013), and would in any case require new simulation suites to address—beyond the scope of this paper.

On the other hand, it is not altogether surprising that our fiducial column density distribution (Section 4.2)—which was measured at z ≈ 3.7 (Prochaska et al. 2010)—do not reproduce the BOSS data at 〈z〉 = 2.3–2.6. We therefore search for an LLS model that better describes the low-transmission end of the BOSS Lyα forest. Looking at the 〈z〉 = 2.3 PDFs in Figure 14, we see that our fiducial model overproduces pixels at F = 0, yet is deficient at slightly higher F. This suggests that our model is overproducing super-LLS ($N_{{\rm H\,\scriptsize{I}}}> 10^{19}\ \mathrm{cm^{-2}}$) that contribute large absorption troughs with F = 0, while not providing sufficient lower-column density absorbers that can individually reach minima of 0.1  ≲  F  ≲  0.4 when smoothed to BOSS resolution. In other words, our fiducial model appears to have an excessively "top-heavy" LLS column density distribution.

For a change, we will try an LLS column density distribution with a more ample bottom end, using the steepest power laws within the 1σ limits estimated by Prochaska et al. (2010):

$$\begin{equation} f(N_{{\rm H\,\scriptsize{I}}}) =\left\lbrace \begin{array}{@{}l@{\quad }l@{}}k_1 N_{{\rm H\,\scriptsize{I}}}^{-1.2} & \mathrm{if}\; 10^{17.5}\,{<}\,N_{{\rm H\,\scriptsize{I}}}\,{<}\,10^{19.0} \\ k_2 N_{{\rm H\,\scriptsize{I}}}^{-1.4} & \mathrm{if}\; 10^{19.0}\,{<}\,N_{{\rm H\,\scriptsize{I}}}\,{<}\,10^{20.3}\end{array}\right.. \end{equation} \tag{ 22 }$$

We use the same l_LLS(z) as before, and obey the integral constraints from Prochaska et al. (2010) that demand that the ratios of $\int f(N_{{\rm H\,\scriptsize{I}}})\; {d}N_{{\rm H\,\scriptsize{I}}}$ between the two column-density regimes be fixed. This gives us k₁ = 10^2.819 and k₂ = 10^7.039, although the new distribution is no longer continuous at $N_{{\rm H\,\scriptsize{I}}}= 10^{19} \ \mathrm{cm^{-2}}$. This new distribution is illustrated by the red power laws in Figure 15.

Zoom In Zoom Out Reset image size

Another change we have made is to the partial LLS model, which was possibly too conservative in the fiducial model. Instead of extrapolating from the LLS distribution, we now adopt the pLLS power-law slope of β_pLLS = −2.0 inferred from the total mean free path to ionizing photons by Prochaska et al. (2010). This dramatically increases the incidence of pLLS in our spectra relative to LLS: we now have l_pLLS = 1.8 l_LLS, where l_LLS is the same value we used previously (Equation (10)). This increase, while large, is not unreasonable in light of the large uncertainties in direct measurements on the H i column-density distribution from direct Lyα line-profile fitting (e.g., Janknecht et al. 2006; Rudie et al. 2013). Note also that even this increased pLLS incidence only amounts to, on average, less than one pLLS per quasar (Δ(z)  ∼  0.3–0.4 per quasar at our redshifts).

We found that while increasing the number of pLLS relieves the tension between data and model at 0.1  ≲  F  ≲  0.4, it does not resolve the excess at the fully absorbed F ≈ 0 pixels in the models. However, changing the b-parameter of the LLS and pLLS from our original fiducial value of b = 70 km s⁻¹ modifies the PDF in a way that improves the agreement. This is a reasonable step, since the effective b-parameter is otherwise observationally ill-constrained for the LLS and pLLS populations. This is because LLSs are typically complexes of multiple systems separated in velocity space, and while there have been analyses of the b-parameter in these individual components, the "effective" b-parameter for complete LLS systems has never been quantified to our knowledge.

We therefore search for the best-fit b-parameter with respect to the T_REF, γ = 1.3 model at 〈z〉 = 2.3, focusing primarily on the agreement in the 0  ⩽  F  ⩽  0.4 bins (Figure 16). Our choice of model for this purpose should not significantly affect our subsequent conclusions regarding the IGM temperature–density slope, since there is little sensitivity toward the latter in the relevant low-transmission bins (compare with Figure 14). However, there will be some degeneracy between the LLS b-parameter and T₀ (Figure 17) since changing the latter does somewhat change the low-transmission portion of the PDF—we will come back to this point in Section 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As shown in Figure 16, a value of b = 45 km s⁻¹ gives the best agreement with the data at 0  ⩽  F  ⩽  0.4. This yields χ² = 116 for 24 dof, which is dramatically improved over that quoted in Figure 14, but still not quite a good fit. In the subsequent results, we will adopt this steeper pLLS/LLS model and b-parameter as the fiducial model in our analysis, and will correspondingly decrease the dof in our χ² analysis to account for the fitting of b.

Note that while significantly improving the PDF fit, this new b-parameter still does not give a perfect fit to the low-transmission (F < 0.4) end. This is probably due to the simplified nature of our LLS model, which neglects the finite distribution of b-parameters and internal velocity dispersion of individual components. These properties are currently not well known, and it seems likely that an improved model would allow a better fit to the low-transmission end of the PDF.

While we have estimated the sample variance of our BOSS transmission PDFs by bootstrap resampling on the spectra, there are significant uncertainties associated with each component of our transmission PDF model as described above, e.g., the LLS incidence rate and level of continuum error. These uncertainties can be incorporated into a systematics covariance matrix, C_sys that can then be added to the bootstrap covariance, C_boot, when computing the model likelihoods. This requires assuming that C_sys and C_boot are uncorrelated, and that the errors are Gaussian distributed.

We adopt a Monte Carlo approach to estimate C_sys by generating 200 model transmission PDFs that randomly vary the systematics. We then evaluate the covariance of the transmission PDFs, p_i, relative to the fiducial model, p_{ref, i} at each transmission bin i. This allows us to construct a covariance matrix with the elements

$$\begin{equation} C_{\mathrm{sys},ij} = \langle (p_i - p_{\mathrm{ref},i})(p_j - p_{\mathrm{ref},j}) \rangle \end{equation} \tag{ 23 }$$

that encompasses the errors from the uncertainties in the LLS model, metal absorption, and continuum scatter. Note that estimation of systematic uncertainties is typically a subjective process, and for most of these contributions we can only make educated guesses as to their uncertainty.

Our Monte Carlo iterations sample the various components of our model as follows:

• 1.  
  LLS incidence. We sample the uncertainty in the power-law exponent γ_LLS of the redshift evolution in LLS incidence rate (Equation (10)), which is $\sigma _{\gamma _{\rm LLS}}\pm 0.21$ as reported by Ribaudo et al. (2011). We assume this uncertainty is Gaussian and draw l_LLS(z) accordingly. This primarily affects the low-flux regions −0.1  ≲  F  ≲  0.3 of the PDF.
• 2.  
  Partial-LLS slope. Our choice of slope for the distribution of partial LLS ($N_{{\rm H\,\scriptsize{I}}}<10^{17.5}\ \mathrm{cm^{-2}}$ absorbers is from an indirect constraint with significant uncertainty (Prochaska et al. 2010). We therefore vary the pLLS slope around the fiducial β_pLLS = −2.0 by ±0.5 assuming a flat prior in this range, which primarily alters the 0  ≲  F  ≲  0.4 portion of the PDF since pLLS typically do not saturate at BOSS resolution.
• 3.  
  LLS b-parameters. Also in the previous section, we found that a global b-parameter of b = 45 km s⁻¹ gives the best agreement with the data, but this is an ad hoc approach with significant uncertainties. In our Monte Carlo sampling we therefore adopt a conservative b = 45 km s⁻¹ ± 20 km s⁻¹ with a uniform prior. This primarily affects the PDF at −0.1  ⩽  F  ⩽  0.4 as can be seen in Figure 16.
• 4.  
  Intervening metals. Although we used an empirical method to model intervening metals (Section 4.4), we may have missed metals with rest wavelengths λ  ≲  1300 Å. Furthermore, we have a relatively small set (∼300–400) of "template" quasars from which our metal model is derived, which may contribute some sampling variance. We therefore guess at a Gaussian error of ±30% for the metal incidence rate. This modulates the extent to which metals pulls the overall PDF toward lower F-values (compare with Figure 9(c)).
• 5.  
  Continuum errors. The overall rms scatter in our continuum estimation also affect the flux PDF (Figure 9(e)). This can be varied in our model by rescaling the quantity c'(λ)/c(λ) − 1, where c is the "true" continuum used to generate the mock spectrum, while c' is the model continuum that we subsequently fit (Figure 13). For each iteration in our Monte Carlo systematics estimation, we dilate or reduce c'(λ)/c(λ) − 1 by a Gaussian deviate assuming ±20% scatter. This primarily affects the high-transmission (F  >  0.8) end of the PDF.

For these Monte Carlo iterations, we used the identical thermal model (γ = 1.6, T_REF) as well as fixed the same random number seeds used for the selection of simulation skewers and generation of noise vectors in our spectra, in order to ensure that the only variation between the different iterations are from the randomly sampled systematics. Figure 18 shows 50 of these Monte Carlo iterations on the transmission PDF for the 〈z〉 = 2.3, S/N = 8–10 subsample.

Zoom In Zoom Out Reset image size

Figure 19 shows an example of the systematic contribution to the covariance matrix. The overall amplitude of the systematic contribution is considerably higher than that estimated from the bootstrap resampling (compare with Figure 7), indicating that we are in the systematics-limited regime. We also see significant anti-correlations at almost the same level as the positive correlations, which are due mostly to correlations between transmission bins on either side of "pivot points" as the transmission PDF varies from the systematics—these anti-correlations will somewhat counteract the increased size of the diagonal components. In the subsequent analysis, we will use an error covariance matrix, C = C_boot + C_sys, in which the systematics covariance matrix estimated in this subsection is added to the bootstrap covariance matrix (described in Section 3) estimated from the BOSS transmission PDFs.

Zoom In Zoom Out Reset image size

We have at this point yet to address one more parameter that can significantly change the shape of our model transmission PDFs, namely the Lyα forest mean transmission assumed in the mock spectra, 〈F〉_Lyα. However, this is an important astrophysical parameter which we did not want to treat as a "systematic," so the next subsection will describe our treatment of 〈F〉_Lyα.

In the initial comparison of the model transmission PDFs shown in Figure 14, the models show a discrepancy with the data at higher transmission bins F ≳ 0.6. Such differences can be alleviated by varying the mean transmission of the pure Lyα forest, 〈F〉_Lyα ≡ exp (− τ_Lyα), i.e., ignoring the contribution from metals and LLS. This quantity can be varied directly in the simulation skewers (Section 4.1). When we vary 〈F〉_Lyα in the simulations, the quantity 〈F〉_cont, which is used to normalize the continuum level of the mock quasar spectrum, is always kept fixed to 〈F〉_eff(z) = exp [ − (τ_Lyα + τ_metals + τ_LLS)] as derived from Becker et al. (2013; see Section 3.2). However, since we are applying the same 〈F〉_cont to both the real and mock spectra, 〈F〉_cont can be best thought of as a normalization that does not actually need to match 〈F〉_eff. Once both the real and mock spectra have been normalized by 〈F〉_cont, the transmission PDF retains information on the respective contributions from the Lyα forest, metals, and LLSs regardless of the assumed 〈F〉_cont, because these contributions affect the shape of the PDF in different ways. In principle, it is possible to vary these all components to infer their relative contributions, but due to the crudeness of our metal and LLS models, we choose to have only 〈F〉_Lyα as a free parameter while keeping 〈F〉_metals = exp (− τ_metals) and 〈F〉_LLS = exp (− τ_LLS) fixed. The possible variation of these latter two components are instead incorporated into the systematic uncertainties determined in Section 5.3. The effect of varying 〈F〉_Lyα is illustrated in Figure 20, where we plot the same IGM model with different underlying values of 〈F〉_Lyα in the simulation skewers while keeping fixed the contribution from metals, LLSs, etc.

Zoom In Zoom Out Reset image size

We therefore explore a range of 〈F〉_Lyα around the vicinity of that estimated by Becker et al. (2013), 〈F〉_{Lyα, B13}, and at each value of 〈F〉_Lyα evaluate the χ² summed over all the S/N subsamples for each 〈z〉 and γ combination. In addition, we now adopt the updated LLS/pLLS model described in Section 5.2, while the χ² evaluation now uses the full covariance matrix including both the bootstrap and systematics (Section 5.3) uncertainties to compare with the transmission PDFs measured from the BOSS data.

The models are compared with the BOSS data as we vary 〈F〉_Lyα, and for each 〈F〉_Lyα we compute the total chi-squared summed over all three S/N subsamples, where each subsample contributes 27 − 1 − 2 = 24 dof (compare with Equations (20) and (21)) along with a further reduction of one dof since we have effectively fitted for the LLS b-parameters in Section 5.2, for a total of ν = 71 dof. The result of this exercise is shown in Figure 21 which shows the χ² values for the T_REF models with different γ—we only vary γ and not T₀ because the F ≳ 0.6 portions of the transmission PDF that change the most with 〈F〉_Lyα do not vary as much with respect to changes in T₀ (compare with Figure 17). Examples of the corresponding best-fit model PDFs in one S/N subsample are shown in Figure 22, where we see that varying 〈F〉_Lyα can indeed change the shape of the F ≳ 0.6 portion of the transmission PDF sufficiently, improving the fits in those transmission ranges compared to the fiducial models (Figure 14).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In all our redshift bins, the best-fitting models seen in Figure 21 are γ = 1.6 with χ² = [69, 67, 54] for 70 dof²⁴ at 〈z〉 = [2.3, 2.6, 3.0] (for the combined data using all S/N bins), respectively, implying probabilities of P = [52%, 59%, and 92%] of obtaining larger values.²⁵ At the higher redshifts best-fitting mean transmission for the γ = 1.6 case is pushed to significantly discrepant values with respect to the fiducial Becker et al. (2013) values (Figure 21).

The γ = 1.3 model also provides acceptable fits to the models, with χ² = [71, 73, 58] for 70 dof (P = [43%, 40%, 84%]) at 〈z〉 = [2.3, 2.6, 3.0], but at the two higher redshift bins this requires 〈F〉_Lyα values that are increasingly discrepant compared to Becker et al. (2013) (+2.3σ and +5σ respectively at 〈z〉 = [2.6, 3.0]). The isothermal γ = 1.0 models are disfavored at the two lower redshift bins, with best-fit values of χ² = [98, 97] for 70 dof (P = [2%, 2%]) at 〈z〉 = [2.3, 2.6], whereas at 〈z〉 = 3, the error bars on the PDF are sufficiently large that acceptable fits are obtainable using γ = 1.0, with χ² = 68 for 70 dof (P = 54%). However, this requires a +5σ discrepancy in 〈F〉_Lyα with respect to Becker et al. (2013). In Figure 22, one sees that fitting for 〈F〉_Lyα allows the γ = 1.0 models to be in good agreement with the data in the F  >  0.7 portion of the PDF, but gives rise to discrepancies in the 0.4  ≲  F  ≲  0.7 range, which limits the goodness-of-fit, and cannot easily be compensated by modifying the metals or LLS model.

From Figure 21, it is clear that as we move to higher redshifts, we require increasingly higher 〈F〉_Lyα relative to the fiducial Becker et al. (2013) values in order to agree with the data: at 〈z〉 = 2.3, our best-fit mean transmission for the γ = 1.6 model agrees with Becker et al. (2013), but at 〈z〉 = 3 there is a significant deviation of +2σ with respect to the Becker et al. (2013) measurement. The same trend is true for the best-fit γ = 1.3 and γ = 1.0 models, but these require even greater discrepancies with respect to the fiducial 〈F〉_Lyα.

One possible explanation for this discrepancy is the effect on the Becker et al. (2013) measurement of u-band selection bias in the SDSS quasars. This was first noted by Worseck & Prochaska (2011), who found that the color–color criteria used to select SDSS quasars preferentially selected quasars, specifically in the redshift range 3  ≲  z_qso  ≲  3.5, that have intervening Lyman-breaks at λ_rest  <  912 Å. The 3  ≲  z_qso  ≲  3.5 SDSS quasars are thus more likely to have intervening LLSs in their sightlines, yielding an additional contribution to the Lyα absorption and hence causing Becker et al. (2013) to possibly underestimate 〈F〉_Lyα when stacking the impacted quasars. Becker et al. (2013) mentioned this effect in their paper but argued that it was much smaller than their estimated errors by referencing theoretical IGM transmission curves estimated by Worseck & Prochaska (2011; Figure 17 in the latter paper).

Dr. G. Worseck has kindly provided us with these transmission curves, T_IGM(λ), which were generated for both the average IGM absorption and that extracted from SDSS quasars affected by the color–color selection bias. In Figure 23 we plot the relative difference between the biased Lyα transmission deduced from z_qso = 3.2 and z_qso = 3.4 quasars and the true mean IGM transmission, using the Worseck & Prochaska (2011) transmission curves. It is clear that at Lyα absorption redshifts of z_abs ≈ 3, the excess LLSs picked up from such quasars contribute an additional ∼1% compared to the mean IGM decrement, a discrepancy that is of the same magnitude as the error bars in the Becker et al. (2013) measurement, indicated by the dashed line.

Zoom In Zoom Out Reset image size

This could partially explain the higher 〈F〉_Lyα required to make our 〈z〉 = 3 models fit the data in Figure 21. Note that we expect this UV color selection bias to be much less significant in our BOSS data, since we have selected bright quasars in the top 5th percentile of the S/N distribution. Given that such quasars have high-S/N photometry, their colors separate much more cleanly from stellar contaminants. Furthermore, such bright quasars are much more likely to have been selected with multi-wavelength data (e.g., including near-IR and radio in addition to optical photometry see Ross et al. 2012). For both of these reasons, we expect our quasars to be much less susceptible to biases in color-selection related to the presence of an LLS. A careful accounting of this bias is beyond the scope of this paper, but from now on we will inflate by a factor of two the corresponding errors on 〈F〉_Lyα at 〈z〉 = 3 to account for this possible bias in the mean transmission measurements (dashed vertical lines in bottom panel of Figure 21).

Another possibility that could explain a bias in the 〈F〉_Lyα measured by Becker et al. (2013) is their assumption that the metal contamination of the Lyα forest does not evolve with redshift. While there are few clear constraints on the aggregate metal contamination within the forest, assuming that the metals actually decrease with increasing redshift (e.g., in the case of C iv, Cooksey et al. 2013), the assumption of an unevolving metal contribution calibrated at z ≈ 2.3 would lead to an underestimate of 〈F〉_Lyα at higher redshifts, which could explain the trend we seem to be seeing.

It is clear from the previous discussion that there is some degeneracy between γ and 〈F〉_Lyα in our transmission PDFs. However, we are primarily interested in γ, while the 〈F〉_Lyα has been extensively measured over the years, which allows the placement of strong priors. In the next section, we will therefore marginalize over 〈F〉_Lyα in order to obtain our final results.

How does this compare with other results on the thermal state of the IGM? McDonald et al. (2001) analyzed the transmission PDF from eight high-resolution, high S/N spectra and compared with them now-obsolete hydrodynamical simulations. They found the data to be consistent with a temperature–density relationship with the expected values of γ ≈ 1.5 (Hui & Gnedin 1997). More recently, Bolton et al. (2008) and Viel et al. (2009) carried out analyses of the transmission PDF measured from a larger sample (18 spectra) of Lyα forest sightlines measured by Kim et al. (2007) and found evidence for an inverted temperature–density relationship (γ  <  1). Viel et al. (2009) found that at z ≈ 3.0, the temperature–density relation was highly inverted (γ ≈ 0.5), and remained so as low as z ≈ 2.0 although at the lower redshifts the data was marginally consistent with an isothermal IGM. They suggested that the difference between their results and those of McDonald et al. (2001) was due to the now-obsolescent cosmological parameters and less-detailed treatment of intervening metals in the earlier study. However, Lee (2012) then pointed out that there is a sensitivity of the measured values of γ from the transmission PDF on continuum-fitting. Since continuum-fitting of high-resolution data generally involves manually placing the continuum at Lyα forest transmission peaks that do not necessarily reach the true continuum, it is conceivable that continuum biases combined with underestimated jacknife errors bars (e.g., Rollinde et al. 2013) could have led Bolton et al. (2008) and Viel et al. (2009) to erroneously deduce an inverted temperature–density relationship (see Bolton et al. 2014 for a detailed discussion on this point). In our analysis we have fitted our continua using an automated process that is free from the same continuum-fitting bias, although it does require an assumption on the underlying Lyα forest transmission, which we have marginalized over in our analysis.

Most recent measurements of the transmission PDF from high-resolution data have continued to favor an isothermal or inverted γ—Calura et al. (2012) analyzed the transmission PDF from a sample of z ≈ 3.3–3.8 quasars and also found an isothermal temperature–density relationship at z = 3, although combining with the Kim et al. (2007) data drove the estimated γ to inverted values at z  <  3. However, Rollinde et al. (2013) carried out a re-analysis of the transmission PDF from various high-resolution echelle data sets, which included significant overlap with the Kim et al. (2007) data. They argue that previous analyses have underestimated the error on the transmission PDF, and found the observed transmission PDF to be consistent with simulations that have γ ≈ 1.4 over 2  <  z  <  3—this discrepancy is probably also driven by a different continuum-estimation from the Kim et al. (2007) measurement.

The use of other statistics on high-resolution spectra have, however, tended to disfavor an isothermal or inverted temperature–density relationship. Rudie et al. (2012) analyzed the lower-end of the b–$N_{{\rm H\,\scriptsize{I}}}$ cutoff from individual Lyα forest absorbers measured in a set of 15 very high-S/N quasar echelle spectra, and estimated γ ≈ 1.5 at z = 2.4. Bolton et al. (2014) compared the Rudie et al. (2012) measurements to hydrodynamical simulations and corroborated their determination of the temperature–density relationship slope.

Garzilli et al. (2012) analyzed the Kim et al. (2007) sample and found that while the transmission PDF supports an isothermal or inverted temperature–density relationship, a wavelet analysis favors γ  >  1. Note, however, that the b–$N_{{\rm H\,\scriptsize{I}}}$ cutoff and the transmission PDF are sensitive to different density ranges, with the PDF probing gas densities predominantly below the mean (e.g., Bolton et al. 2014).

Our result of γ ≈ 1.6 at 〈z〉 = [2.3, 2.6, 3.0] is thus in rough agreement with measurements that do not involve the transmission PDF from high-resolution Lyα forest spectra (with the exception of Rollinde et al. 2013). Our value of γ at 〈z〉 = 3 is somewhat unexpected because one expects a flattening of the temperature–density relationship close to the He ii reionization epoch at z  ∼  3 (Furlanetto & Oh 2008; McQuinn et al. 2009; but see Gleser et al. 2005; Meiksin & Tittley 2012), although γ = 1.3 is not strongly disfavored ($\sqrt{(\Delta \chi ^2)}\,{\sim}\,2.6$).

Taken at face value, the temperature–density relationship during He ii reionization can be made steeper by a density-independent reionization and/or a lower heating rate in the IGM (Furlanetto & Oh 2008), which could be reconciled with an extended He ii event (Shull et al. 2010; Worseck et al. 2011).

Our constraints on γ appear to be in conflict with the prediction of the theories of Broderick et al. (2012) and Chang et al. (2012), who elucidated a relativistic pair-beam channel for plasma-instability heating of the IGM from TeV gamma-rays produced by a population of luminous blazars. This mechanism provides a uniform volumetric heating rate, which would cause an inverted temperature–density relationship in the IGM (Puchwein et al. 2012) since voids would experience a higher specific heating rate compared with heating by He ii reionization alone. This picture has been challenged by the recent study of Sironi & Giannios (2014), who dispute the amount of heating this mechanism could provide, since they found that the momentum dispersion of such relativistic pair beams allows ≪10% of the beam energies to be deposited into the IGM.

However, in this paper we have assumed relatively simple temperature–density relationships in which the bulk of the IGM in the density range 0.1  ≲  Δ  ≲  5 follows a relatively tight power law. We have therefore not studied more complicated T–Δ relationships, e.g., with a spread of temperatures at fixed density (e.g., Meiksin & Tittley 2012; Compostella et al. 2013) that might be caused by He ii reionization or other phenomena. It is therefore possible that such complicated temperature–density relationships could result in Lyα forest transmission PDFs that mimic the γ ≈ 1.6 power law; this is something that needs to be examined in more detail in future work.

Looking forward, the subsequent BOSS data releases will significantly enlarge our sample size, e.g., DR10 (Ahn et al. 2014) is nearly double the size of the DR9 sample used in this paper, while the final BOSS sample (DR12) should be three times as large as DR9. In particular, the newer data sets should be sufficiently large for us to analyze the transmission PDF and constrain γ during the epoch of He ii reionization at z  >  3. This would be a valuable measurement, since high-resolution spectra are particularly affected by continuum-fitting biases at these redshifts (Faucher-Giguère et al. 2008; Lee 2012).

The analysis of the optically thin Lyα forest transmission PDF from these expanded data sets will have vanishingly small sample errors, and the errors will be dominated by systematic and astrophysical uncertainties. At the high-transmission end, our uncertainties are dominated by the scatter of the continuum-fitting, which is dominated by the question of whether our quasar PCA templates, derived from low-luminosity low-redshift quasars (Suzuki et al. 2005), or high-luminosity SDSS quasars (Pâris et al. 2011), respectively, are an accurate representation of the BOSS quasars. This uncertainty should be eliminated in the near future by PCA templates derived self-consistently from the BOSS data (N. Suzuki et al. 2014, in preparation). The modeling of metal contamination could also be improved in the near future by advances in our understanding of how metals are distributed in the IGM (e.g., Zhu et al. 2014), although metals are a comparatively minor contribution to the uncertainty in our transmission PDF.

We also aim to improve on the rather ad hoc data analysis in this paper, in which we accounted for some uncertainties in our modeling by incorporating them into our error covariances (e.g., LLSs, metals, continuum errors), while 〈F〉_Lyα was marginalized over a fixed grid. In future analyses, it would make sense to carry out a full MCMC treatment of all these parameters which would rigorously account for all the uncertainties and allow straightforward marginalization over nuisance parameters.

Since this paper was initially focused on modeling the BOSS spectra, for the model comparison we used only simulations sampling a very coarse 3 × 3 grid in T₀ and γ parameter space, and were unable to take account for uncertainties in other cosmological (σ₈, n_s etc) and astrophysical (e.g., Jeans scale, Rorai et al. 2013; or galactic winds, Viel et al. 2013b) parameters in our analysis. However, methods already exist to interpolate Lyα forest statistics from hydrodynamical simulations given a set of IGM and cosmological parameters (e.g., Viel & Haehnelt 2006; Borde et al. 2014; Rorai et al. 2013). In the near future we expect to perform joint analyses using other Lyα forest statistics in conjunction with the transmission PDF, such as new measurements of the small-scale (k ≳ 0.2 s km⁻¹) 1D transmission power spectrum (M. Walther et al. 2014, in preparation), moderate-scale (0.002 s km⁻¹  ≲  k  ≲  0.2 s km⁻¹) transmission power spectrum in both 1D (e.g., Palanque-Delabrouille et al. 2013) and 3D (from ultra-dense Lyα forest surveys using high-redshift star-forming galaxies, Lee et al. 2014a, 2014b), the phase angle PDF determined from close quasar pair sightlines (Rorai et al. 2013), and others. Such efforts would require a fine grid sampling the full set of cosmological and IGM thermal parameters in order to ensure that the interpolation errors are small compared to the uncertainties in the data (see e.g., Rorai et al. 2013). Efforts are underway to utilize massively parallel adaptive-mesh refinement codes (Almgren et al. 2013) to generate such parameter grids to study the IGM (Lukić et al. 2014) However, one of the findings of this paper is the importance of correct modeling of LLSs, in particular partial LLSs ($10^{16.5}\ \mathrm{cm^{-2}}\,{\lesssim}\,N_{{\rm H\,\scriptsize{I}}}\,{\lesssim}\,10^{17.5}\ \mathrm{cm^{-2}}$), in accounting for the shape of the observed Lyα transmission PDF. Since our hydrodynamical simulations did not include radiative transfer and cannot accurately capture optically thick systems, we had to add these in an ad hoc manner based on observational constraints that are currently rather imprecise. In the near future, we would want to use hydrodynamical simulations with radiative transfer (even if only in post-processing, e.g., Altay et al. 2011; McQuinn et al. 2011; Altay et al. 2013; Rahmati et al. 2013) to self-consistently model the optically thick absorbers in the IGM. With the unprecedented statistical power of the full BOSS Lyα forest sample, this could provide the opportunity to place unique constraints on the column-density distribution function of partial LLSs.