The Opacity of the Intergalactic Medium Measured along Quasar Sightlines at z ∼ 6

Determining when and how the epoch of reionization proceeded is one of the major goals of observational cosmology today. During this early evolutionary phase of our universe, the cosmic "dark ages" following recombination ended, and the intergalactic medium (IGM) transitioned from a neutral state into the ionized medium that we observe today due to the ultraviolet radiation of the first stars, galaxies, and quasars. The details of the reionization process not only reflect the nature of these primordial objects, but also the formation of large-scale structure and are therefore a subject of major interest. Despite much progress in the last decade, there are still crucial yet unanswered questions regarding the exact timing and morphology of reionization.

The most compelling constraints to date on the end of the reionization epoch come from the evolution of the Lyman-alpha (Lyα) forest opacity, observed in the spectra of z ≳ 6 quasars (Fan et al. 2006; Becker et al. 2015; Bosman et al. 2018). The detection of transmitted flux spikes in the high redshift quasar spectra and the absence of large Gunn–Peterson (GP) troughs (Gunn & Peterson 1965) below z ≲ 5.0 indicate that the epoch of reionization must be completed by that time. Both the steep rise in the observed opacity around z ≳ 5.5 as well as the increased scatter of the measurements suggest a qualitative change in the ionization state of the IGM (Becker et al. 2015), provoked by a decrease in the ionizing ultraviolet background (UVB) radiation (Calverley et al. 2011; Wyithe & Bolton 2011).

The inferred rapid decline of the UVB radiation has been interpreted as an indication for the end stages of reionization (e.g., Fan et al. 2006; Bolton & Haehnelt 2007; Calverley et al. 2011). Fan et al. (2006) argue further that the increased scatter in the opacity measurements around z ≳ 5.5 could be explained by strong variations in the UVB radiation field as expected in patchy and inhomogeneous reionization scenarios.

However, Lidz et al. (2006) have argued that large-scale density fluctuations alone could explain the significant variations between sightlines. They calculated the scatter arising solely from density fluctuations while assuming a uniform UVB, which gives results comparable to the observations by Fan et al. (2006). If this were correct, the evidence for patchy reionization based on the observations of the mean opacity would be significantly weakened.

Measurements of the Lyα forest opacity along additional quasar sightlines by Becker et al. (2015) finally showed that the observed scatter in the optical depth measurements significantly exceeds not only fluctuations expected from the density field alone, but also fluctuating UVB models with a spatially uniform mean free path of ionizing photons. They posited that, if a fluctuating UVB was in fact the source of the large scatter in optical depth, the mean free path must be spatially variable, supporting the interpretation of probing the end stages of an inhomogeneous hydrogen reionization period. Subsequently, Davies & Furlanetto (2016) modeled the UVB with a spatially varying mean free path and found that the additional fluctuations were sufficient to explain the extra scatter in the optical depth measurements (see also Chardin et al. 2017).

An alternative explanation for the observed scatter in the mean opacity was presented by D'Aloisio et al. (2015), who showed that residual spatial fluctuations in the temperature field could result in extended opacity variations. The spatially varying temperature field is a natural consequence of an extended an inhomogeneous reionization process, wherein regions that reionized early have had time to cool, while regions that reionized late are still hot. The amplitude of the resulting opacity variations then depends directly on the timing and duration of the reionization process.

The largest sample of quasar sightlines to date used for IGM opacity measurements was recently presented by Bosman et al. (2018). They compared their findings to IGM models that included either a fluctuating UVB or temperature fluctuations and concluded that neither fully captures the observed scatter in IGM opacity.

In this paper we present a new data set of 34 high redshift quasar spectra at 5.77 ≤ z_em ≤ 6.54, which we make publicly available, and which presents an ideal laboratory for studying the epoch of reionization and setting constraints on the timing and the morphology of the reionization process. For a subset of 23 quasar sightlines we present new measurements of the evolution of the mean opacity of the IGM within the Lyα forest between $4.8\lesssim z\lesssim 6.3$. We carefully mask all spectral regions that could be biasing our measurements—the region in the immediate vicinity of the quasar, its so-called proximity zone where the transmitted flux is enhanced due to the radiation from the quasar itself, as well as all patches exhibiting additional absorption due to intervening absorption systems such as damped Lyα absorbers (DLAs), along the line of sight to the quasars. Additionally we correct for possible offsets in the zero level of the quasar spectra due to potential systematics in the data reduction procedure.

We then compare our opacity measurements to predictions from a hydrodynamical simulation for three different cases—assuming a uniform UVB radiation field, a fluctuating UVB field, and a fluctuating temperature field—in order to describe the observed evolution of the IGM opacity, and to assess the excess of inhomogeneities in the density, radiation or temperature field that would be required to explain our measurements.

This paper is organized as follows: in Section 2 we describe our data set and its properties. In this section we also present the data release of this quasar sample within the igmspec⁴ database. The methods we use to continuum normalize the quasar spectra and measure the IGM opacity are described in Section 3. The results of the opacity measurements and their evolution with redshift are presented in Section 4. We compare our measurements to different outputs from a hydrodynamical simulation, which are described in Section 6. The implications for the epoch of reionization are discussed in Section 7, and we conclude in Section 8 with a summary of the main results. Throughout the paper we assume a cosmology of h = 0.685, Ω_m = 0.3 and Ω_Λ = 0.7, which is consistent within the 1σ errorbars in Planck Collaboration et al. (2016).

In this section we describe the properties of the data sample of quasar spectra that we use for our analysis of the IGM opacity. This data set has been previously introduced in Eilers et al. (2017a), in which we conducted a detailed analysis of the proximity zones of the quasars. We briefly summarize the details of the observations in Section 2.1, and the data reduction procedure in Section 2.2. The properties of this quasar sample are described in Section 2.3, before presenting the data release at the end of Section 2.4.

2.1. Properties of the Data Set

2.2. Data Reduction

A detailed description of the data reduction can be found in Eilers et al. (2017a) and will be summarized here only briefly. All spectra were reduced uniformly using the ESIRedux pipeline⁶ developed as part of the XIDL⁷ suite of astronomical routines in the Interactive Data Language. This pipeline employs standard data reduction techniques comprising the following: images are overscan subtracted, flat fielded using a normalized flat field image, and then wavelength calibrated by means of a wavelength image constructed from afternoon arc lamp calibration images. After identifying the objects in the science frames, we subtracted the background using B-spline fits (Kelson 2003; Bernstein et al. 2015) to the object-free regions of the slit. The profiles of the science objects were also fit with B-splines, and an optimal extraction was performed on the sky-subtracted frames. We combined 1D spectra of overlapping echelle orders to produce a spectrum for each exposure, and co-added individual exposures into a final 1D spectrum. For a more detailed description of the applied algorithms, see Bochanski et al. (2009). We further optimized the XIDL ESI pipeline to improve the data reduction for quasars at high redshift by differentiating two images (ideally taken during the same observing run) with similar exposure times, analogous to the standard difference imaging techniques performed for near-infrared observations, in order to improve the sky subtraction especially in the reddest echelle orders, which are affected by fringing. This procedure requires dithered exposures for which the trace of the science object lands at different spatial locations on the slit. Dithered exposures have the additional advantage that different parts of the fringing pattern are being sampled, and hence a combination of different exposures further reduces fringing issues in the data. However, since not every observer dithered their object along the slit, it was not possible for us to apply this procedure to ≈10% of the exposures that we took from the archive.

We co-added exposures from different observing runs taken by various PIs to maximize the S/N. To combine the data from different observing runs, we weighted each 1D spectrum by its squared signal-to-noise ratio (S/N²), determined in the quasar continuum region of each spectrum, i.e., at wavelengths longer than the Lyα emission line. This ensures that spectral regions with low or no transmitted flux, which are common in high redshift quasar spectra, obtain the same weight as regions with more transmitted flux.

All final reduced quasar spectra are shown in Figure 1, sorted by their emission redshift.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.3. Quasar Properties

The properties of all quasars in our data sample, such as the emission redshift z_em, the absolute magnitude M₁₄₅₀ at λ_rest = 1450 Å in the rest frame, the size of the proximity zone R_p of the quasars, and the S/N of their spectra, are presented in Table 1.

Table 1.  Overview of Our Data Sample

Object	RAhms	DECdms	zem	M1450	Rp [pMpc]	S/Na	Opacity?
SDSS J0002+2550	00h02m3939	+25°50'3496	5.82	−27.31	5.43 ± 1.49	59	yes
SDSS J0005−0006	00h05m5234	−00°06'5580	5.844	−25.73	2.87 ± 0.40	13	yes
CFHQS J0050+3445	00h55m0291	+34°45'2165	6.253	−26.7	4.09 ± 0.37	18	yes
SDSS J0100+2802	01h00m1302	+28°02'2592	6.3258	−29.14	7.12 ± 0.13	39	yes
ULAS J0148+0600	01h48m3764	+06°00'2006	5.98	−27.39	6.03 ± 0.39	27	yes
ULAS J0203+0012	02h03m3238	+00°12'2927	5.72	−26.26	⋯	9	no
CFHQS J0210−0456	02h10m1319	−04°56'2090	6.4323	−24.53	1.32 ± 0.13	1	no
PSO J036+03	02h26m0187	+03°02'5942	6.5412	−27.33	3.64 ± 0.13	11	yes
CFHQS J0227−0605	02h27m4329	−06°05'3020	6.20	−25.28	1.60 ± 1.37	3	no
SDSS J0303−0019	03h03m3140	−00°19'1290	6.078	−25.56	2.21 ± 0.38	2	no
SDSS J0353+0104	03h53m4973	+01°04'0466	6.072	−26.43	⋯	13	no
PSO J060+25	04h02m1269	+24°51'2443	6.18	−26.95	4.17 ± 1.38	12	yes
SDSS J0818+1722	08h18m2740	+17°22'5201	6.02	−27.52	5.89 ± 1.42	7	yes
SDSS J0836+0054	08h36m4386	+00°54'5326	5.81	−27.75	5.06 ± 0.40	112	yes
SDSS J0840+5624	08h40m3530	+56°24'2022	5.8441	−27.24	7.39 ± 0.30b	26	yes
SDSS J0842+1218	08h42m2943	+12°18'5058	6.069	−26.91	6.47 ± 0.38	10	yes
SDSS J0927+2001	09h27m2182	+20°01'2364	5.7722	−26.76	4.68 ± 0.15	7	no
SDSS J1030+0524	10h30m2711	+05°24'5506	6.309	−26.99	5.93 ± 0.36	24	yes
SDSS J1048+4637	10h48m4507	+46°37'1855	6.2284	−27.24	⋯	43	no
SDSS J1137+3549	11h37m1773	+35°49'5685	6.03	−27.36	6.98 ± 1.42	22	yes
SDSS J1148+5251	11h48m1665	+52°51'5039	6.4189	−27.62	4.58 ± 0.13	35	yes
SDSS J1250+3130	12h50m5193	+31°30'2190	6.15	−26.53	6.59 ± 1.38	8	yes
SDSS J1306+0356	13h06m0827	+03°59'2636	6.016	−26.81	5.39 ± 0.38	37	yes
ULAS J1319+0950	13h19m1130	+09°50'5152	6.133	−27.05	3.84 ± 0.14	10	yes
SDSS J1335+3533	13h35m5081	+35°33'1582	5.9012	−26.67	0.78 ± 0.15	6	no
SDSS J1411+1217	14h11m1129	+12°17'3728	5.904	−26.69	4.60 ± 0.39	29	yes
SDSS J1602+4228	16h02m5398	+42°28'2494	6.09	−26.94	7.11 ± 1.40	21	yes
SDSS J1623+3112	16h23m3181	+31°12'0053	6.2572	−26.55	5.05 ± 0.14	9	yes
SDSS J1630+4012	16h30m3390	+40°12'0969	6.065	−26.19	4.80 ± 0.38	15	yes
CFHQS J1641+3755	16h41m2173	+37°55'2015	6.047	−25.67	3.98 ± 0.38	3	no
SDSS J2054−0005	20h54m0649	−00°05'1480	6.0391	−26.21	3.17 ± 0.14	15	yes
CFHQS J2229+1457	22h29m0165	+14°57'0900	6.1517	−24.78	0.45 ± 0.14	2	no
SDSS J2315−0023	23h15m4657	−00°23'5810	6.117	−25.66	3.70 ± 1.39	14	yes
CFHQS J2329−0301	23h29m0828	−03°01'5880	6.417	−25.25	2.45 ± 0.35	2	no

Notes. The columns show the object name, the coordinates of the quasar given in RAhms and DECdms, the emission redshift and the quasar's magnitude M₁₄₅₀, the measurements for their proximity zones, and the S/N of the quasar spectrum. The last column states whether we used the quasar sightline for the IGM opacity analysis in this paper.

^aMedian S/N per 10 km s⁻¹ pixel; estimated between 1250 Å ≤ λ_rest ≤ 1280 Å. ^bEstimate of the proximity zone conservatively assuming Δv = 5000 km s⁻¹, because the proximity zone is prematurely truncated due to associated absorption systems (see Appendix A in Eilers et al. 2017a).

Download table as:  ASCIITypeset image

The proximity zone measurements are taken from Eilers et al. (2017a). The uncertainties of these measurements arise solely from uncertainties in the redshift estimate, since these errors provide the largest source of uncertainty for the proximity zone measurements. For one object, SDSS J0840+5624 we do not report a measurement of its proximity zone, because associated absorption systems located in the immediate vicinity of the quasar prematurely truncate its proximity zone (see Appendix A in Eilers et al. 2017a). Thus, we assume conservatively a region of Δv = 5000 km s⁻¹ to be within the influence of the quasar's radiation, resulting in an effective R_p = 7.39 ± 0.30.

For three objects (ULAS J0203+0012, SDSS J0353+0104, and SDSS J1048+4637) no estimates of their proximity zones have been determined. These quasars show broad absorption line (BAL) features, which make a precise and unbiased estimate of their proximity zones unfeasible. Due to the additional absorption features in their spectra we exclude these objects from the analysis of the IGM opacity.

2.4. Data Release

In order to measure the IGM opacity the quasar spectra must be normalized to their unabsorbed continua. In this section we explain our method for continuum normalizing the quasar spectra (Section 3.1) and present the details of the IGM opacity measurements (Section 3.2). In the last part of this section (Section 3.3) we describe a procedure for correcting small offsets in the zero level of the quasar spectra.

3.1. Quasar Continuum Normalization

We fit the quasar continua and normalize the spectra in a similar manner as that previously conducted by Eilers et al. (2017a). We summarize the main steps here briefly, but refer the reader to the previous paper for more details.

All spectra were normalized to unity at λ_rest = 1280 Å in a spectral region that is free of emission lines. The quasar continuum was then estimated with principal component spectra (PCS) from a principal component analysis (PCA) of lower redshift quasar spectra (Pâris et al. 2011). The concept of PCA is to represent each continuum spectrum $| {q}_{\lambda }\rangle$ at wavelength λ by a reconstructed spectrum comprising a mean spectrum $| {\mu }_{\lambda }\rangle$ and a sum of m weighted PCS $| {\xi }_{\lambda }\rangle$, i.e.,

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

where the index i refers to the ith PCS and α_i denotes its weight.

Since the quasars in our data sample are all at z_em ∼ 6, they experience substantial absorption due to intervening neutral hydrogen within the IGM blueward of the Lyα emission line. Thus, the continuum estimate was performed solely on the unabsorbed continuum spectrum redward of the Lyα emission line with a set of PCS from Pâris et al. (2011) covering the wavelengths 1215.67 Å ≤ λ_rest ≤ 2000 Å. We take the model that minimizes χ² using the noise vector from the spectra as the best estimate.

In order to obtain coefficients ${\boldsymbol{\alpha }}$ for a set of PCS that cover the entire spectral region between 1020 Å ≤ λ_rest ≤ 2000 Å, we use a projection matrix ${\boldsymbol{P}}$ to project the estimated coefficients for the PCS redward of Lyα onto this new set of coefficients for the entire spectrum. The projection matrix ${\boldsymbol{P}}$ has been computed by Pâris et al. (2011) using the set of PCS for both the red wavelength side only and the whole spectral region covering wavelengths blueward and redward of Lyα. Hence,

$$\begin{eqnarray}&&{\boldsymbol{\alpha }}={\boldsymbol{P}}\cdot {{\boldsymbol{\alpha }}}_{\mathrm{red}}.\end{eqnarray} \tag{ 2 }$$

This new set of coefficients together with Equation (1) provides a continuum model for each quasar covering all wavelengths 1020 Å ≤ λ_rest ≤ 2000 Å. Pâris et al. (2011) estimate that the median uncertainty on the transmitted flux in the Lyα forest to be $| {\rm{\Delta }}{F}_{\mathrm{forest}}| \approx 5 \%$. However, since we do not take all PCA components into account and do not have the full spectral coverage up to λ_rest = 2000 Å to estimate the continua, the uncertainty on the continua in our quasar spectra is most likely higher, i.e., $| {\rm{\Delta }}{F}_{\mathrm{forest}}| \approx 10 \% \mbox{--}20 \%$.

For an estimate of the continua at lower wavelength we take the composite quasar spectrum provided by Shull et al. (2012), constructed from 22 low redshift quasars observed with the Cosmic Origins Spectrograph on the Hubble Space Telescope (HST) that extends from 550–1750 Å in the rest frame, and rescale the composite spectrum to match the PCA-constructed continuum model at λ_rest = 1020 Å. We augment the continuum model by simply appending the composite spectrum at wavelengths λ_rest < 1020 Å. Note, however, that we do not use the spectrum at λ_rest < 1020 Å for the analysis of the Lyα forest opacity. A few example quasar spectra from our data set and its continuum model are shown in Figure 2, and all remaining quasar spectra that we analyzed and their estimated continua are shown in Figure 12 in Appendix C.

Zoom In Zoom Out Reset image size

Note that in some cases the N v line at λ_rest = 1240 Å is not very well represented by the continuum fit. This behavior occurs when the N v line is slightly blue- or redshifted compared to the systemic redshift of the quasar and those shifts are not accounted for in the PCA basis. Additionally, the continuum for quasar spectra with very weak emission lines, such as SDSS J0100+2802 or SDSS J1148+5251, is not very well captured by the PCA.

3.2. Measuring the Optical Depth of the IGM

We exclude BAL quasars (ULAS J0203+0012, SDSS J0353+0104, and SDSS J1048+4637) from our IGM optical depth sample to avoid potential contamination by broad non-IGM absorption in the Lyα forest, and quasars with only very low S/N data, i.e., S/N < 7, whose spectra are more subject to systematic errors. Our final IGM optical depth sample then consists of 23 quasar spectra out of the original 34.

We estimate the mean opacity of the IGM by means of the effective optical depth, which is defined as

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

where F is the continuum normalized flux. The effective optical depth is measured in discrete spectral bins along the line of sight of each quasar. We chose fixed comoving bins of size 50 comoving Mpc h⁻¹ ($\mathrm{cMpc}\,{h}^{-1}$) (Becker et al. 2015), which contains a similar path length as the bins of size Δz = 0.15 at z ∼ 5–6 previously applied by Fan et al. (2006).

In order to avoid biases in the measurement of the opacity of the IGM, we mask the spectral region around each quasar that is strongly influenced and ionized by the quasar's own radiation. We use the measurements for the proximity zones R_p as an estimate the influenced region. However, the proximity zone is defined such that it does not completely reach out to the ionization front expanding from the quasar. Thus, the influence of the quasar's radiation is expected to be still present outside of its measured proximity zone (Eilers et al. 2017a), since the radiation of the quasar still dominates the UVB radiation, i.e., Γ_QSO ≫ Γ_UVB at R_p. Hence, we mask an additional 2.5 proper Mpc (pMpc) around each quasar, i.e., the masked region measures R_p+2.5 pMpc, in order to eliminate all enhanced transmission due to the quasar's radiation.

Thus, we choose the maximum wavelength that we consider for opacity measurements to lie just blueward of this masked region. The minimum wavelength we consider for measurements within the Lyα forest is 1030 Å in the rest frame, in order to account for possible redshift uncertainties.

Another possible contamination in the measurement of optical depths are intervening absorbers along the line of sight, such as DLAs or other low-ion metal absorbers, which are likely associated with relatively high H i column density (N_{H i} ≳ 10¹⁹ cm⁻²). We mask the regions in the quasar spectra around these absorbers, as they reflect an absorption signature that is not typically resolved in IGM simulations. To this end, we searched for low-ion metal absorption lines, such as, e.g., Al ii, Fe ii, and O i, associated with absorbers in the continuum spectra redward of the Lyα emission line that are located at the same redshift as a spectral region showing complete absorption in the forest of the spectrum. Additionally, we searched the literature for DLAs and low-ion metal absorbers along the quasar sightlines in our sample. For each absorber, we conservatively masked the spectral region around the absorption system within ±30 Å in the observed wavelength frame in the Lyα forest at the corresponding wavelength. Spectral bins containing absorbers were then excluded from the IGM opacity measurements. Table 2 shows a compilation of all identified absorbers along the line of sight to the quasars in our sample. Note that most of the identified absorbers have already been found by other authors, since most quasars in our data set were previously known and observed.

Table 2.  Intervening Low-ion Absorption Systems along the Line of Sight

Object	zem	zabs	Referencesa
SDSS J1148+5251	6.4189	6.0115	1
		6.1312	1
		6.1988	1
		6.2575	1
SDSS J2054−0005	6.0391	5.9776	1
SDSS J2315−0023	6.117	5.7529	1
SDSS J1630+4012	6.065	5.8865	1
SDSS J1137+3549	6.03	5.0124	1
SDSS J1623+3112	6.2572	5.8415	1
SDSS J0840+5624	5.8441	5.5940	2
SDSS J0002+2550	5.82	4.914	3
SDSS J0100+2802	6.3258	6.1437	3
SDSS J0818+1722	6.02	5.7911	1
		5.8765	1

Note. The columns show the name of the object and its emission redshift, the absorption redshift of the intervening low-ion absorber, and the reference therefor.

^aReference for low-ion absorbers. 1: Becker et al. (2011), 2: Ryan-Weber et al. (2009), 3: this work.

Download table as:  ASCIITypeset image

Additionally we mask all spuriously high pixels within the Lyα forest of the quasar spectra by checking for single pixels showing F > 1 in the continuum normalized spectra. Because sky-subtraction systematics occasionally result in large negative sky-subtraction residuals, we also mask all negative flux pixels with the 2.5% lowest S/N to avoid biases due to large uncertainties in pixels that fall onto sky lines and have large negative residuals.

After masking all low-ion metal absorption systems, proximity zones, and spuriously high and negative pixels, the combined usable path length for the opacity measurements is 6350 $\mathrm{cMpc}\,{h}^{-1}$ for the 23 quasar sightlines in our data sample. The spectral regions in which we measure the mean flux and calculate its effective optical depth are shown for each quasar as the dark and light blue colored bars in Figure 3. Masked regions are shown in white. The gray regions show path lengths that are in principle usable but are not used, because the remaining unmaksed region would be smaller than our chosen bin size of 50 $\mathrm{cMpc}\,{h}^{-1}$.

Zoom In Zoom Out Reset image size

3.3. Correcting for Offsets in the Zero Level of the Spectra

The noise for pixels with no intrinsic flux should be symmetrically distributed around zero, since pixels with zero flux have equal probability to be scattered into positive or negative values. This idea was applied by McGreer et al. (2011), for example, to estimate the number of so-called "dark pixels" that have a flux value consistent with zero. However, a detailed inspection of the quasar spectra in our data set reveals that the zero level in some spectral regions can be slightly biased, i.e., we do not observe a symmetric distribution around zero in the noise, possibly caused by sky-subtraction systematics present in a small fraction of exposures. But even tiny offsets in the zero level can cause large differences in the opacity estimates, especially in regions with very little transmitted flux. Because these are the regions we are particularly interested in, we correct for small offsets in the zero level of the spectra. These offsets are calculated and applied to each 50 $\mathrm{cMpc}\,{h}^{-1}$ spectral bin individually in order to avoid correlations between the optical depth measurements.

A detailed description of this correction procedure can be found in Appendix B. We estimate that the systematic error in the mean flux $\langle {F}^{\mathrm{Ly}\alpha }\rangle$ due to corrections of the zero-level offset is

$$\begin{eqnarray*}&&{\sigma }_{\langle {F}^{\mathrm{Ly}\alpha }\rangle }\approx 0.0067,\end{eqnarray*}$$

i.e., less than <1%, and thus constitutes only a minor correction to the optical depth measurements.

We compute the effective optical depth τ_eff from the measurements of the observed mean flux $\langle {F}^{\mathrm{obs}}\rangle$ in bins of 50 cMpc h⁻¹ using Equation (3). We list all measurements of the mean flux within the Lyα forest for each bin along all 23 quasar sightlines in our data sample in Table 5 in Appendix D. All spectral bins indicating the respective measurements of $\langle {F}^{\mathrm{obs}}\rangle$ and ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\alpha }$ are shown in Figure 13 in Appendix E. If the mean flux is detected with less than 2σ significance, or if we measure a negative mean flux, we adopt a lower limit on the optical depth at ${\tau }_{\mathrm{eff}}=-\mathrm{ln}(2{\sigma }_{\langle {F}^{\mathrm{obs}}\rangle })$ consistent with previous works. Note that we do not include the systematic uncertainty of the mean flux (∼10%–20%, see Section 3.1) arising from the PCA continuum estimate. This uncertainty of the mean flux would lead to an additional uncertainty of ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\alpha }$ of ∼2%–5%, when most of the flux is absorbed in the quasar spectra (see also Figure 7 in Becker et al. 2015). The results of the optical depth measurements within the Lyα forest, plotted as a function of redshift, are shown in Figure 4.

Zoom In Zoom Out Reset image size

In Figure 4 our new measurements are shown in dark blue in both panels. For the majority of quasar spectra that we use in our analysis the IGM opacity has been analyzed before. However, we have co-added the data from multiple observation runs (see Section 2.2 for details) in order to achieve higher S/N. However, for one object, ULAS J0148+0600, the data obtained by Becker et al. (2015) with VLT/X-Shooter in a 10 hr observation, has a higher S/N than our spectrum. This sightline exhibits a particularly deep GP trough (Becker et al. 2015), and hence the enhanced S/N results in more stringent opacity limits, representing the strongest fluctuations in the IGM opacity at this redshift. In order to model the IGM fluctuations correctly, it is important to include these outliers (Chardin et al. 2015; D'Aloisio et al. 2015; Davies & Furlanetto 2016). Thus, we construct a master compilation of opacity measurements, presented in the left panel of Figure 4, and replace our optical depth measurements within the Lyα forest along just the sightline of ULAS J0148+0600 with the more precise measurements obtained by Becker et al. (2015) (orange data points). This mainly adds the two most stringent limits at z = 5.634 and z = 5.796 to our analysis, since the better data quality results in higher sensitivity in the GP troughs. The lower redshift ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\alpha }$ measurements for this object are consistent with our measurements. This master compilation is shown in the left panel of Figure 4. Additionally, we also show the lower redshift ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\alpha }$ measurements from Becker et al. (2015) (green data points).

We present the average opacity evolution by calculating the mean flux and the bootstrapped error on the mean in bins of Δz = 0.25 within the Lyα forest from the master compilation set, shown in the left panel of Figure 4. We then compute the binned opacity values via ${\tau }_{\mathrm{eff}}=-\mathrm{ln}\langle F\rangle$, where $\langle F\rangle$ is the mean flux computed in these bins. The uncertainties on the opacity values with uncertainties also determined via bootstrapping are shown as the red data points and tabulated in Table 3. Similar to the individual ${\tau }_{\mathrm{eff}}$ measurements, we adopt a limit if the mean flux in the bin is measured with less than 2σ significance (where σ is here the bootstrap errors on the mean flux).

Table 3 .  Measurements of the Average Flux and Optical Depth within the Lyα Forest of Our Master Compilation Sample

zabs	$\langle F\rangle$	${\sigma }_{\langle F\rangle }$	$\langle {\tau }_{\mathrm{eff}}^{\mathrm{Ly}\alpha }\rangle$	${\sigma }_{\langle {\tau }_{\mathrm{eff}}^{\mathrm{Ly}\alpha }\rangle }$
4.0	0.4046	0.0151	0.9049	0.0372
4.25	0.3595	0.0112	1.0230	0.0311
4.5	0.2927	0.0190	1.2286	0.0651
4.75	0.1944	0.0150	1.6381	0.0770
5.0	0.1247	0.0132	2.0818	0.1060
5.25	0.0795	0.0078	2.5321	0.0984
5.5	0.0531	0.0058	2.9357	0.1090
5.75	0.0182	0.0045	4.0057	0.2469
6.0	0.0052	0.0043	>4.7595	⋯
6.25a	−0.0025	0.0007	>6.5843	⋯

Note. The columns show the mean redshift z_abs of the redshift bins of size Δz = 0.25, the averaged flux $\langle F\rangle$ and its uncertainty ${\sigma }_{\langle F\rangle }$ determined via bootstrapping, and the mean optical depth $\langle {\tau }_{\mathrm{eff}}^{\mathrm{Ly}\alpha }\rangle$ in that redshift bin and its error ${\sigma }_{\langle {\tau }_{\mathrm{eff}}^{\mathrm{Ly}\alpha }\rangle }$, also determined via bootstrapping.

^aNote that this redshfit bin only contains two measurements.

Download table as:  ASCIITypeset image

Note that we do not take any systematic errors on the mean flux measurements into account that could, for instance, result from uncertainties in the continuum estimation. The dark gray regions give an estimate of the additional scatter expected due to continuum uncertainties of ∼20%, which are negligible at high redshift, where the transmitted flux is low and the scatter is dominated by fluctuations along different sightlines (Becker et al. 2015; Eilers et al. 2017b).

The right panel of Figure 4 compares our data set to opacity measurements from additional sightlines from the literature that are not in our data sample. The additional data points come from the sightlines of SDSS J0144–0125 and SDSS J1436+5007 (Fan et al. 2006), CFHQS J1509–1749 (Willott et al. 2007), ULAS J1120+0641 (Barnett et al. 2017), PSO J006.1240+39.221 (Tang et al. 2017), and J0323–4701, J0330–4025, J0410–4414, J0454–4448, J0810+5105, J1257+6349, J1609+3041, J1621+5155, and J2310+1855 (Bosman et al. 2018). In most of these analyses the bins were chosen to be Δz = 0.15, following Fan et al. (2006). This bin size covers roughly the same spectral region as the chosen bin size of 50 $\mathrm{cMpc}\,{h}^{-1}$ in our analysis and the one by Becker et al. (2015) at z ∼ 6, but in the redshift interval of 5 ≲ z ≲ 7, the bin size changes quite significantly. Overall the agreement between the various measurements with our new analysis is good, but we chose not to add these measurements to the master compilation, because of the different path lengths used to construct the measurements, very low S/N data, or the variety of different instruments and data reduction pipelines used to obtain the spectra, which enlarges the systematic uncertainties on these measurements (see Section 5).

For several quasar sightlines in our data sample the optical depth has been measured previously by Fan et al. (2006) and Becker et al. (2015), and more recently by Bosman et al. (2018). However, the quality of the data and the methods to analyze the data differ. Here, we carry out a detailed comparison of our methods and measurements to previous work, and discuss potential systematic uncertainties (Section 5.1) and resulting discrepancies in the cumulative distribution functions (CDFs) of the optical depth (Section 5.2).

5.1. Estimating Systematic Uncertainties

We compare the measurements of the IGM opacity for the 16 quasar sightlines that are both part of our analysis and the data set of Fan et al. (2006) and are not BAL quasars. The spectra from Fan et al. (2006) partially overlap with our data set, but six quasars were observed with a different telescope and instrument (MMT Red Channel, Hobby-Eberly Telescope, and KP 4 m MARS) and eight quasars have additional Keck/ESI data. We have reduced and stacked all of the Keck/ESI observations from the archive and thus the spectra analyzed in this paper have an improved quality due to their longer exposure time.

Additionally, our methods to analyze the data differ. For instance, Fan et al. (2006) applied a power law to the red side of the quasar spectra to estimate the quasar continua, whereas we estimated the quasar continua by a PCA (see Section 3.1). In our analysis we mask all spectral regions containing low-ion metal absorption systems, while in previous work it has been argued that those have a negligible influence and can thus be ignored (Fan et al. 2006; Becker et al. 2015).

All these differences contribute to the systematic error of the opacity measurements. We attempted to assess these systematic uncertainties by comparing the results from our analysis to those of Fan et al. (2006) along the sightlines that are part of both data sets. To this end we measure the mean flux in the same redshift bins as Fan et al. (2006) with a fixed bin size of Δz = 0.15, and compare our measurements to those of Fan et al. (2006) in Figure 5. We observe a large scatter in the distribution and a systematic offset toward lower mean flux values in our measurements, much larger than the formal measurement uncertainties. The negative offset is strongest at lower redshifts with higher mean flux values, i.e., lower optical depths.

Zoom In Zoom Out Reset image size

We estimate the systematic error arising due to different observations, different data reduction pipelines, and different analyses by the median of the distribution of measured flux differences ${\rm{\Delta }}\langle F\rangle =\langle {F}_{\mathrm{Eilers}\mathrm{et}\mathrm{al}.2018}\rangle -\langle {F}_{\mathrm{Fan}\mathrm{et}\mathrm{al}.2006}\rangle$. We find a median systematic error of

$$\begin{eqnarray*}&&{\sigma }_{{\rm{\Delta }}\langle F\rangle }\approx -0.023,\end{eqnarray*}$$

with a large scatter of ≈0.026 determined from the mean of the 16th and 84 percentiles of the distribution. A detailed investigation of a few of the largest outliers in this distribution suggests that differences in the spectra itself, due to the different instruments with which they were observed and potentially due to differences in the data reduction, cause the largest discrepancies.

Recently, a similar analysis of the Lyα optical depths measured from a comparable quasar sample at z_em > 5.7 was presented by Bosman et al. (2018). Of the 62 sightlines they analyzed, 22 satisfy our quality criteria, namely, that they are non-BAL quasars with a S/N > 7.⁹ Out of these, 17 objects overlap with our sample. Although the Bosman et al. (2018) sample is comparable to ours in size and partially overlapping, their methods differ in a variety of important aspects from ours. As in Fan et al. (2006) different data reduction pipelines have been used to reduce the data, the quasar continuum estimation methods differ, and while both our study and that of Bosman et al. (2018) mask the proximity zone regions, we adapt the excluded region dependent on the actual measured proximity zone size R_p (see Section 3.2), whereas their analysis excludes a fixed spectral range until λ_rest = 1178 Å, which corresponds to ΔR_p = 13.3 pMpc at z = 6. Finally, we have masked strong absorbers and account for small zero-level offsets, whereas they do not.

5.2. Comparison of the CDFs

In Figure 6 we compare the CDF from our measurements shown in blue to the CDF from previous studies by Fan et al. (2006) and Becker et al. (2015), which are shown as the gray curves, and by Bosman et al. (2018) shown in yellow, in different redshift bins centered around 5.0 ≤ z ≤ 6.0. We show the so-called GOLD sample from Bosman et al. (2018), including 33 quasar spectra for which they applied a data quality cut of S/N > 11.2 per 60 km s⁻¹ pixel to their sample, which would imply a quality cut on our sample of S/N > 4.6 per 10 km s⁻¹ pixel.

Zoom In Zoom Out Reset image size

While previous work noticed an increased scatter in the opacity measurements only at z ≳ 5.5, we also see evidence for increased scatter at 5.0 < z < 5.5. We see systematic differences toward higher optical depths in our work compared to others, most strikingly in the 5.3 < z < 5.7 bins whose excess fluctuations have been the focus of several works. However, in most redshift bins the measurements agree within 1σ uncertainties (shown as the shaded regions in Figure 6), which we determined via bootstrap resampling, the only exception being the redshift bin at z = 5.4 and z = 5.6, where our results are slightly more discrepant with previous studies.

A discrepancy between the Fan et al. (2006) and Becker et al. (2015) measurements in this bin was previously noted by Chardin et al. (2017), and in particular it seems that our (higher) ${\tau }_{\mathrm{eff}}$ measurements are more consistent with the data from Becker et al. (2015) than those from Fan et al. (2006).

We would like to compare our measurements of the IGM opacity to expectations from simulations. For this purpose we use a hydrodynamical simulation, which we briefly describe in section Section 6.1. In this simulation we use a uniform UVB radiation field and thus this simulations provides a good approximation for opacity fluctuations in the IGM long after the epoch of reionization, when we expect to have a uniform UVB. We use two more sophisticated models to compare our measurements with more realistic conditions in the post-reionization IGM. To this end, we use two semi-numerical models with fluctuating UVB and temperature field, which we describe in Section 6.2. In Section 6.3 we explain how we compute the Lyα optical depth from the skewers through the simulation box.

6.1. Nyx Hydrodynamical Simulation

6.2. Semi-numerical Models with Fluctuating UVB and Temperature Fields

6.3. Calculating the ${Ly}\alpha$ Optical Depth from Simulated Skewers

We then extract 1000 skewers through the various simulation boxes to compute the optical depths and compare the results to our measurements. Because the exact strength of the UVB radiation Γ_UVB is unknown, we have to rescale the optical depth in each skewer at each pixel i, i.e., ${\tau }_{i}^{\mathrm{Ly}\alpha ,\mathrm{unscaled}}$, to match the mean optical depth corresponding to the observed mean flux value $\langle {F}^{\mathrm{obs}}\rangle$ of our measurements, which in turn depends on the exact value of Γ_UVB. Hence, at each redshift we determine a scaling factor A₀ that solves the following equation:

$$\begin{eqnarray}\langle \exp [-{\tau }_{i}^{\mathrm{Ly}\alpha }]\rangle & = & \langle \exp [-{A}_{0}\times {\tau }_{i}^{\mathrm{Ly}\alpha ,\mathrm{unscaled}}]\rangle \\ & = & \ \langle {F}^{\mathrm{obs}}\rangle .\end{eqnarray} \tag{ 4 }$$

We then average the rescaled flux at each pixel $\langle \exp [-{\tau }_{i}^{\mathrm{Ly}\alpha }]\rangle$ within each skewer of size 50 $\mathrm{cMpc}\,{h}^{-1}$, and determine the 68% and 95% of the distribution of ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\alpha }$. This gives an estimate of the expected scatter within the Lyα optical depth measurements given a uniform UVB and IGM thermal state, that is, arising from density fluctuations alone.

In order to assess the implications of our opacity measurements for the epoch of reionization, we compare our measurements to the opacity distributions from the Nyx hydrodynamical simulation (Section 6.1). The Nyx simulation was computed with a uniform UVB, which we would expect long after the epoch of reionization or when assuming no signatures of an inhomogeneous reionization process. The evolution of the optical depth distributions from these simulations, with the mean fluxes matched to our measurements, are shown as the gray region in Figure 4. The width of the light (68th percentile) and medium (95th percentile) gray regions corresponds to the scatter in τ_eff expected due to fluctuations in the underlying density field alone. The dark gray regions indicate the additional scatter expected from ∼20% uncertainties in the quasar continuum estimation. These have been calculated by dividing the mean flux in each spectral bin by (1 + ΔC), where the continuum uncertainty ΔC was drawn randomly from a normal distribution with σ = 0.2 and μ = 0, corresponding to ∼20% uncertainties in the continuum estimate. As expected continuum uncertainties matter only very little at higher redshifts, when the mean transmitted flux is low and fluctuations between different quasar sightlines dominate, and the scatter at high redshift can thus not be explained by continuum uncertainties.

The measurements show a steep rise in ${\tau }_{\mathrm{eff}}$ for z ≳ 5 and an increased scatter in the distribution of measurements. At lower redshifts (z ≲ 5) the scatter in the observed ${\tau }_{\mathrm{eff}}$ decreases rapidly and becomes consistent with the expectations from density fluctuations alone. It is evident, however, that at high redshifts the scatter in the optical depths measurements significantly exceeds the scatter expected from density fluctuations alone, i.e., the scatter represented by our hydrodynamical simulation with uniform UVB. The tiny rare flux spikes observed in the Lyα forests of SDSS J0100+2802 and SDSS J1148+5251 at redshifts 5.8 ≲ z ≲ 5.9 are in strong contrast to the abundant transmitted flux along the sightlines toward SDSS J1306+0356 or SDSS J2054–0005 at similar redshifts, for example. We show the respective spectral regions exhibiting very high (upper panels), average (middle panels), and very low (lower panels) optical depths at similar redshifts of the aforementioned sightlines in Figure 7.

The discrepancy between our measured opacity distribution and the expectation from simulations with a uniform UVB becomes even more evident in Figure 6, where we show the cumulative distributions of ${\tau }_{\mathrm{eff}}^{\mathrm{Ly}\alpha }$ from our master compilation in different redshift bins. The CDFs of ${\tau }_{\mathrm{eff}}$ from our hydrodynamical simulation + uniform UVB are shown as the dashed red curves, where we have rescaled the pixel optical depths (Section 6.3) to match the mean ${\tau }_{\mathrm{eff}}$ in each redshift bin. This model with a uniform UVB is clearly not a good match to the observations. While they are more consistent with the measurements at lower redshift (z ∼ 5), there are large discrepancies at higher redshifts (z ≳ 5.6) between the simulated and the observed CDF, a point previously noted by Fan et al. (2006) and Becker et al. (2015). While it may seem that the model provides a better fit in the highest redshift bin at z ∼ 6.0, this apparent agreement is misleading, and arises due to the fact that we show limits on the optical depth in the same way as measurements, and the bin at z ∼ 6 contains several limits.

7.1. Comparison to Patchy Reionization Models

Multiple scenarios have been proposed to explain the increased scatter in the optical depth relative to the fluctuations one would expect from the density field of the IGM alone. One possible explanation is that the UVB is strongly fluctuating, either due to coupled variations in the mean free path of ionizing photons (Davies & Furlanetto 2016; D'Aloisio et al. 2018) or a rare source population, such as quasars (Chardin et al. 2015, 2017). Another possibility is that the thermal state of the IGM is highly inhomogeneous (Lidz & Malloy 2014; D'Aloisio et al. 2015). Such fluctuations can arise as a result of an extended and patchy reionization process, where different regions in the universe were reionized (and simultaneously photoheated) at different redshifts z_reion. The regions that reionized earlier would have had time to cool down, while regions that were reionized at a later time would still be relatively hot.

In Figure 8 we compare our measurements at 5.6 < z < 5.8 to the semi-numerical models with a fluctuating UVB and fluctuating temperature field (see Section 6.2), where we have rescaled the opacities in the Lyα forest skewers from that work to match the mean ${\tau }_{\mathrm{eff}}$ we have measured in this bin. Note that while previously the CDF of optical depths at this redshift bin containing the strong outliers in the opacity measurements in the GP trough along the sightline of ULAS J0148+0600 (Becker et al. 2015) was very challenging to reproduce in simulations (D'Aloisio et al. 2015; Davies & Furlanetto 2016; Chardin et al. 2017; Davies et al. 2018a; Keating et al. 2018), these outliers are now easier to explain because the mean ${\tau }_{\mathrm{eff}}$ in our measurements is substantially higher than the Fan et al. (2006)+Becker et al. (2015) compilation, and hence these data points do not represent such strong deviations from the mean of the distribution anymore. This first comparison of our measurements to the two semi-numerical models with a fluctuating UVB and a fluctuating thermal state of the IGM shows that both models can reproduce the observations. A more detailed comparison to these models will be part of future work.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We present a new data set of 34 quasar spectra at 5.77 ≤ z_em ≤ 6.54 that we make publicly available via the igmspec database. The spectra have all been observed with ESI on the Keck telescopes, and exposures from different observing runs have been co-added, resulting in a very rich and homogeneous data set, with a total of ∼180 hr of telescope time.

For a subsample of 23 quasar spectra, that do not show BAL features and have good quality data (i.e., S/N > 7), we measure the IGM opacity by means of the effective optical depth of the Lyα forest in bins of 50 $\mathrm{cMpc}\,{h}^{-1}$, covering a redshift range of 4.0 ≲ z ≲ 6.5. Our results are in qualitative agreement with previous studies (Fan et al. 2006; Becker et al. 2015; Bosman et al. 2018), showing a steep rise in opacity and increased scatter within the measurements at high redshift. However, while previous work noticed an increased scatter at z ≳ 5.5, we also see evidence for increased scatter at 5.0 < z < 5.5. A detailed comparison in the optical depth in several redshift bins, shows systematic differences toward higher optical depths in our work compared to others, most strikingly at 5.3 < z < 5.7. The discrepancies, however, between our measurements and previous work are mostly within the ∼1σ uncertainties, which we determined via bootstrap resampling.

Our work improves upon previous studies in several aspects. We considered possible contamination due to intervening low-ion metal absorption systems such as DLAs that have previously been ignored and carefully masked all regions in the Lyα forest that are affected by these high H i column density absorption systems. We also corrected for small offsets in the zero level of the quasar spectra, introduced presumably by improper sky subtraction of a few individual exposures. Finally and most importantly, we considered a very homogeneously reduced data sample that minimizes systematic effects due to the use of different telescopes and detectors, or data reduction pipelines. We present a master compilation sample, including mainly our newly analyzed sample with the exception of the sightline of ULAS J0148+0600 taken from Becker et al. (2015), which has a larger sensitivity in the prominent GP trough due to the higher S/N spectrum observed with VLT/X-Shooter.

We compare our measurements to a large-volume hydrodynamical simulation with a uniform UVB. As noted previously by Fan et al. (2006) and Becker et al. (2015), we find that the spread in observed ${\tau }_{\mathrm{eff}}$ cannot be explained by fluctuations of the underlying density field alone, and thus our results support an inhomogeneous reionization scenario. Whether temperature fluctuations in the IGM, a fluctuating UVB, or a combination of both can best explain this increased scatter in opacity, remains an open question. A preliminary comparison of our measurements to semi-numerical simulations of UVB and IGM temperature fluctuations shows good agreement for both scenarios.

This work presents a crucial ingredient in constraining the end stages of the epoch of reionization at 5.0 ≲ z ≲ 6.0, when the physical conditions of the post-reionization IGM can be directly measured via absorption spectroscopy of high redshift quasars. The past several years have seen an impressive fivefold increase in the number of z > 6 quasars from deep wide-field optical and infrared surveys, which are enabling precise measurements of the Lyα forest absorption at 5 < z < 6.5 (Becker et al. 2015; Gnedin et al. 2017; Davies et al. 2018a). The requirement that reionization models reproduce these high-precision measurements provides an important low redshift anchor point that all models must reproduce, and can dramatically narrow the family of viable reionization models. Statistical anaylses of the Lyα forest, such as measurements of the power spectrum (Oñorbe et al. 2017; D'Aloisio et al. 2018), or the probability distribution function (PDF) of the IGM opacity (Davies et al. 2018b), set further constraints on the reionization process, allowing us to develop accurate models about the early evolutionary stages of our universe.

We would like to thank the anonymous referee for insightful and constructive feedback. We thank José Oñorbe for helpful discussions regarding the hydrodynamical simulations, Zarija Lukić for providing the Nyx simulation, and Vikram Khaire for valuable feedback. Furthermore, we appreciate the help of Xiaohui Fan for making his quasar spectra available to us for a data comparison.

The data presented in this paper were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation.

This research has made use of the Keck Observatory Archive, which is operated by the W.M. Keck Observatory and the NASA Exoplanet Science Institute, under contract with the National Aeronautics and Space Administration.

The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have had the opportunity to conduct observations from this mountain.

Software: ESIRedux (http://www2.keck.hawaii.edu/inst/esi/ESIRedux/), XIDL (http://www.ucolick.org/xavier/IDL/), igmspec (http://specdb.readthedocs.io/en/latest/igmspec.html), Nyx (Almgren et al. 2013), emcee (Foreman-Mackey et al. 2013), numpy (van der Walt et al. 2011), scipy (Jones et al. 2001), matplotlib (Hunter 2007), astropy (The Astropy Collaboration et al. 2018), and seaborn (Waskom et al. 2017).

In Table 4 we present the catalog with the properties of our data set that will be available together with the final co-added spectra, their noise vectors, and their continuum estimates via the igmspec database.¹⁰ igmspec is a database of publicly available ultraviolet, optical, and near-infrared spectra that probe the IGM. It provides ∼500,000 spectra from various data sets, including the Sloan Digital Sky Survey (SDSS), 2dF QSO Redshift Survey (2QZ), and data from the HST, the Keck Telescopes, the Very Large Telescope, and more. The database is part of the specdb repository¹¹ within the specdb Github organization¹² , which provides software developed in Python for accessing and interacting with the quasar spectra.

In order to correct for possible offsets in the zero level of the quasar spectra, we examine the negative pixels in each spectral bin, which should appear to be a truncated Gaussian distribution, providing an estimate of the noise level in the spectra. To this end, we take all flux pixels F below zero in each spectral bin of 50 $\mathrm{cMpc}\,{h}^{-1}$, and calculate the CDF, ignoring correlations between neighboring pixels.

In the case of no offset in the zero level within a spectral bin, the estimated mean of the CDF μ_CDF should be equal to zero, as shown by the toy model example in the lower middle panel of Figure 9. The upper panels show the PDF of the same respective toy model case. However, if the zero level is slightly under- or overestimated (left and right panels in Figure 9, respectively), the estimated μ_CDF tracking the true zero level will likewise be below or above zero. It is clear that we can obtain a handle on these systematic offsets by examining the purely negative pixels and fitting a truncated CDF model.

Zoom In Zoom Out Reset image size

We use a Markov chain Monte Carlo (MCMC) approach making use of the implementation of the affine-invariant ensemble sampler emcee¹³ (Foreman-Mackey et al. 2013) to estimate the mean μ_CDF of the best fitting model CDF

$$\begin{eqnarray}&&{\mathrm{CDF}}_{\mathrm{model}}(F)=\displaystyle \frac{{A}_{\mathrm{CDF}}}{2}\left[1+\mathrm{erf}\left(\displaystyle \frac{F-{\mu }_{\mathrm{CDF}}}{\sqrt{2}{\sigma }_{\mathrm{CDF}}}\right)\right],\end{eqnarray} \tag{ 5 }$$

while marginalizing over the width of the distribution σ_CDF and the scaling factor A_CDF. The likelihood function just maximizes the least squares between the CDF model and the measured CDF, i.e.,

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=-0.5{({\mathrm{CDF}}_{\mathrm{model}}-{\mathrm{CDF}}_{\mathrm{data}})}^{2}.\end{eqnarray} \tag{ 6 }$$

We then take the median of the resulting posterior PDF as the best estimate for μ_CDF.

The three free parameters of the CDF_model, μ_CDF, σ_CDF, and A_CDF are highly degenerate with each other because we are fitting only a small part of the CDF when taking solely flux pixels below zero into account, i.e., F < 0. Thus, we have to apply strict priors to break this degeneracy. The priors we chose are flat priors within the intervals ${\mu }_{\mathrm{CDF}}\in [-0.05,0.05]$, since we expect the total offset of the zero level to be less than ±5%, ${\sigma }_{\mathrm{CDF}}\in [0.75{\sigma }_{\mathrm{eff}},1.25{\sigma }_{\mathrm{eff}}]$, which takes into account the noise vector at each pixel i of each quasar spectrum to estimate ${\sigma }_{\mathrm{eff}}=\sqrt{{\sum }_{i}{\sigma }_{i}^{2}/N}$, and ${A}_{\mathrm{CDF}}\in [0.45\,{N}_{F\lt 0},0.55\,{N}_{F\lt 0}]$. The upper and lower boundaries for A_CDF result from the fact that the number of pixels with flux below zero, i.e., N_F<0, in an unbiased case should be exactly half of the pixels, i.e., 0.5 N_F<0. In the presence of the possible offsets in the zero level of the spectra we allow A_CDF to deviate from the unbiased case.

Figure 10 shows two examples of the procedure. Both panels show the CDF for two spectral bins along the sightline of SDDS J0840+5624. We show the actually measured CDF of all negative flux pixels (blue curves) and overplot the best fitted CDF (red dashed curves, with mean μ_CDF indicated by the red dashed-dotted lines) and the ideal CDF with no zero-level offset, i.e., μ_CDF = 0, and the variance given by the noise of the data, i.e., σ_CDF = σ_eff. We can see that we obtain small zero-level offsets of about ${\rm{\Delta }}\langle {F}^{\mathrm{Ly}\alpha }\rangle \sim {\mu }_{\mathrm{CDF}}\approx 0.3 \% \mbox{--}0.4 \%$.

Zoom In Zoom Out Reset image size

In the end, we offset all pixels i within each spectral bin by the respective best estimate for μ_CDF, i.e., F_i,new = F_i − μ_CDF, and calculate the mean flux and the opacity from the offsetted pixels F_i,new.

Figure 11 shows the difference in mean flux estimates due to corrections in the zero level. We estimate this systematic uncertainty in the mean flux measurements from calculating the 84th percentile and 16th percentile of this distribution and taking their average, which results in

$$\begin{eqnarray*}&&{\sigma }_{\langle {F}^{\mathrm{Ly}\alpha }\rangle }=0.0067.\end{eqnarray*}$$

Zoom In Zoom Out Reset image size

In Figure 12 we show the estimates of the quasar continua for all objects, which we use for measurements of the IGM opacity that are not already shown in Figure 2.

Zoom In Zoom Out Reset image size

Our measurements of the mean flux in the Lyα forest are shown in Table 5 for spectral bins of size 50 $\mathrm{cMpc}\,{h}^{-1}$ along all 23 quasar sightlines in our data sample.

Table 5.  Mean Flux Measurements in the Lyα Forest

Object	zstart	zabs	zend	$\langle F\rangle$	${\sigma }_{\langle F\rangle }$
SDSS J0002+2550	5.54	5.62	5.70	0.0472	0.0004
	5.39	5.47	5.54	0.0315	0.0005
	5.25	5.32	5.39	0.0241	0.0004
	5.10	5.17	5.25	0.1170	0.0005
	4.97	5.04	5.10	0.1239	0.0006
SDSS J0005−0006	5.60	5.68	5.76	0.0507	0.0016
	5.45	5.53	5.60	0.0384	0.0028
	5.30	5.38	5.45	0.0484	0.0028
	5.16	5.23	5.30	0.0460	0.0023
	5.02	5.09	5.16	0.1128	0.0031
	4.89	4.95	5.02	0.0732	0.0028
CFHQS J0050+3445	5.96	6.05	6.14	0.0045	0.0018
	5.80	5.88	5.96	0.0217	0.0031
	5.64	5.72	5.80	0.0269	0.0022
	5.49	5.56	5.64	0.0570	0.0020
	5.34	5.41	5.49	0.1336	0.0020
	5.19	5.26	5.34	0.1009	0.0014
SDSS J0100+2802	5.81	5.89	5.98	0.0013	0.0005
	5.65	5.73	5.81	0.0038	0.0003
	5.50	5.57	5.65	0.0416	0.0003
	5.35	5.42	5.50	0.0778	0.0003
ULAS J0148+0600a	5.68	5.76	5.84	0.0056	0.0018
	5.53	5.60	5.68	−0.0020	0.0016
	5.38	5.45	5.53	0.0533	0.0018
	5.23	5.30	5.38	0.0826	0.0014
	5.09	5.16	5.23	0.1428	0.0014
	4.95	5.02	5.09	0.1340	0.0018
PSO J036+03	6.24	6.33	6.42	−0.0018	0.0015
	6.06	6.15	6.24	−0.0032	0.0011
	5.89	5.98	6.06	−0.0021	0.0009
	5.73	5.81	5.89	0.0001	0.0008
	5.57	5.65	5.73	0.0065	0.0006
	5.42	5.50	5.57	0.0632	0.0012
PSO J060+25	5.90	5.98	6.06	0.0308	0.0037
	5.73	5.81	5.90	0.0154	0.0054
	5.57	5.65	5.73	0.0176	0.0020
	5.42	5.50	5.57	0.0256	0.0023
	5.28	5.35	5.42	0.0155	0.0017
	5.13	5.20	5.28	0.0646	0.0016
SDSS J0818+1722	5.56	5.64	5.72	0.0118	0.0034
	5.41	5.49	5.56	0.0194	0.0042
	5.26	5.34	5.41	0.0462	0.0032
	5.12	5.19	5.26	0.0457	0.0028
	4.99	5.05	5.12	0.0768	0.0041
SDSS J0836+0054	5.54	5.62	5.70	0.1061	0.0003
	5.39	5.46	5.54	0.0997	0.0004
	5.24	5.31	5.39	0.0303	0.0003
	5.10	5.17	5.24	0.1104	0.0003
	4.97	5.03	5.10	0.1036	0.0004
	4.83	4.90	4.97	0.1673	0.0003
SDSS J0840+5624	5.39	5.46	5.54	0.0279	0.0012
	5.24	5.31	5.39	0.0778	0.0009
	5.10	5.17	5.24	0.1010	0.0011
	4.96	5.03	5.10	0.0937	0.0014
	4.83	4.90	4.96	0.1487	0.0008
SDSS J0842+1218	5.76	5.84	5.92	−0.0177	0.0055
	5.60	5.68	5.76	0.0160	0.0037
	5.44	5.52	5.60	0.0285	0.0052
	5.30	5.37	5.44	0.0817	0.0048
	5.15	5.22	5.30	0.1302	0.0041
	5.02	5.08	5.15	0.0311	0.0046
SDSS J1030+0524	5.98	6.07	6.16	0.0055	0.0013
	5.82	5.90	5.98	0.0140	0.0020
	5.66	5.74	5.82	0.0114	0.0011
	5.50	5.58	5.66	0.0444	0.0010
	5.35	5.43	5.50	0.1219	0.0012
	5.21	5.28	5.35	0.0561	0.0008
SDSS J1137+3549	5.71	5.79	5.87	0.0056	0.0025
	5.56	5.63	5.71	0.0807	0.0017
	5.40	5.48	5.56	0.1422	0.0020
	5.26	5.33	5.40	0.0913	0.0017
	5.12	5.19	5.26	0.1193	0.0017
SDSS J1148+5251	5.77	5.86	5.94	0.0022	0.0003
	5.61	5.69	5.77	0.0109	0.0001
	5.46	5.54	5.61	0.0123	0.0002
	5.31	5.39	5.46	0.0422	0.0002
SDSS J1250+3130	5.83	5.91	5.99	−0.0212	0.0078
	5.67	5.75	5.83	−0.0091	0.0054
	5.51	5.59	5.67	−0.0099	0.0034
	5.36	5.44	5.51	0.0037	0.0042
	5.22	5.29	5.36	0.0276	0.0030
	5.08	5.15	5.22	0.0644	0.0027
SDSS J1306+0356	5.72	5.80	5.89	0.0934	0.0013
	5.57	5.64	5.72	0.0495	0.0007
	5.42	5.49	5.57	0.0491	0.0008
	5.27	5.34	5.42	0.0436	0.0007
	5.13	5.20	5.27	0.0630	0.0006
	4.99	5.06	5.13	0.0785	0.0008
ULAS J1319+0950	5.86	5.94	6.02	0.0019	0.0058
	5.70	5.78	5.86	0.0002	0.0049
	5.54	5.62	5.70	0.0249	0.0029
	5.39	5.46	5.54	0.0697	0.0036
	5.24	5.31	5.39	0.0406	0.0027
	5.10	5.17	5.24	0.1396	0.0029
SDSS J1411+1217	5.63	5.71	5.79	0.0170	0.0008
	5.48	5.56	5.63	0.0355	0.0010
	5.33	5.40	5.48	0.0543	0.0010
	5.19	5.26	5.33	0.0384	0.0008
	5.05	5.12	5.19	0.1211	0.0010
	4.91	4.98	5.05	0.0429	0.0009
SDSS J1602+4228	5.76	5.85	5.93	0.0242	0.0028
	5.61	5.68	5.76	0.0273	0.0013
	5.45	5.53	5.61	0.0381	0.0014
	5.30	5.38	5.45	0.0504	0.0012
	5.16	5.23	5.30	0.0609	0.0009
	5.02	5.09	5.16	0.0660	0.0011
SDSS J1623+3112	5.95	6.04	6.12	0.0136	0.0031
	5.63	5.71	5.79	0.0276	0.0030
	5.47	5.55	5.63	0.0383	0.0034
	5.32	5.40	5.47	0.0198	0.0029
	5.18	5.25	5.32	0.0225	0.0018
SDSS J1630+4012	5.62	5.70	5.78	−0.0062	0.0023
	5.47	5.54	5.62	0.0193	0.0024
	5.32	5.39	5.47	0.1088	0.0028
	5.17	5.24	5.32	0.0586	0.0020
	5.04	5.10	5.17	0.1274	0.0028
SDSS J2054−0005	5.78	5.86	5.95	0.0287	0.0039
	5.62	5.70	5.78	0.0150	0.0026
	5.47	5.54	5.62	0.1131	0.0031
	5.32	5.39	5.47	0.0652	0.0032
	5.18	5.25	5.32	0.1196	0.0025
	5.04	5.11	5.18	0.2436	0.0037
SDSS J2315−0023	5.84	5.93	6.01	−0.0127	0.0041
	5.53	5.60	5.68	0.0477	0.0037
	5.38	5.45	5.53	0.0348	0.0042
	5.23	5.30	5.38	0.1012	0.0031
	5.09	5.16	5.23	0.1468	0.0031

Note. The different columns show the name of the object, the beginning of each redshift bin z_start, the mean redshift of each bin z_abs and the end of z_end the redshift bin, and the mean flux of the continuum normalized spectrum with its uncertainty.

^aNote that the measurements along this quasar sightline have been replaced by the ones from Becker et al. (2015) in our master compilation.

Download table as:  ASCIITypeset images: 1 2

In Figure 13 we show all spectral bins of 50 $\mathrm{cMpc}\,{h}^{-1}$ along all 23 quasar sightlines in our data sample, indicating the respective optical depth measurements.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size