Measurements of effective optical depth in the Ly$\alpha$ forest from the BOSS DR12 Quasar sample

TL;DR

This study measures the effective optical depth in the Lyα forest using a large quasar sample, providing highly precise estimates of its redshift evolution and examining implications for cosmic reionization.

Contribution

It presents the most precise measurements of the optical depth parameters in the Lyα forest from BOSS DR12 data, with a detailed analysis of systematic uncertainties.

Findings

Power-law model fits the optical depth evolution with $ au_0$ and $eta$ parameters.

Systematic errors dominate the uncertainty in the redshift evolution parameter.

No significant evidence for He II reionization signals was found.

Abstract

We measure the effective optical depth in the Lyman-alpha (Ly$\alpha$) forest using 40,035 quasar spectra from the Twelfth Data Release (DR12) of the Baryon Oscillation Spectroscopic Survey (BOSS) of SDSS-IV. A rigorous selection based on spectral index and the equivalent width of the C IV emission line is applied to choose seven uniform samples with minimal intrinsic variations across redshifts. Modeling the redshift evolution of the effective optical depth with a power-law, $\tau_\mathrm{eff}=\tau_0 (1 + z)^\gamma$ produces $\tau_0=(5.54 \pm 0.64)\times 10^{-3}$ and $\gamma=3.182 \pm 0.074$. The 2.2% precision estimate on $\gamma$ is dominated by systematic errors, likely arising from the bias and uncertainties in spectral index estimates. Even after incorporating the systematic errors, this work provides the most precise estimates of optical depth parameters to date. Finally, using the reconstructed Ly$\alpha$ forest continuum to directly measure the transmitted flux ratio as a function of redshift, we find deviations of less than 2.5% from the predictions from the global model and no convincing evidence for signal associated with He II reionization.

Figures

Tables

Table 1. Table 1 The Number of Quasars Remaining in the Sample After Each Selection Criterion is Applied

Selection	Remaining	Percent
Criterion	quasars	remaining
All	\num297301	–
Objects with ZWARNING = 0	\num283405	95
Remove quasars with BALs	\num257138	86
Remove quasars with DLAs	\num231785	78
$z_q>1.6$	\num163263	55
$S/N>5$	\num58062	19
Cuts in parameter space	\num40035	13

Table 2. Table 2 Best fit values for the optical depth parameters with statistical errors

bin	$\#$ quasars	$\ln {\tau }_0$	$\gamma$	p-value (stat.)	p-value (stat. + sys.)
bin 1	4625	$-5.54\pm 0.15$	$3.42\pm 0.10$	$5.8\times {10}^{-7}$	0.28
bin 2	4203	$-5.33\pm 0.21$	$3.28\pm 0.13$	$1.4\times {10}^{-2}$	0.71
bin 3	7477	$-5.04\pm 0.13$	$3.10\pm 0.08$	$4.5\times {10}^{-6}$	0.69
bin 4	7715	$-4.91\pm 0.15$	$3.00\pm 0.09$	$1.0\times {10}^{-1}$	0.51
bin 5	3855	$-5.58\pm 0.20$	$3.40\pm 0.12$	$9.8\times {10}^{-2}$	0.73
bin 6	8824	$-4.97\pm 0.13$	$3.00\pm 0.08$	$2.3\times {10}^{-8}$	0.36
bin 7	3336	$-6.06\pm 0.36$	$3.59\pm 0.23$	$3.1\times {10}^{-25}$	0.06
Coadd	40,035	$-5.22\pm 0.06$	$3.19\pm 0.04$
Coadd (w/o bin 7)	36,699	$-5.13\pm 0.06$	$3.14\pm 0.04$
Coadd with systematic errors	40,035	$-5.31\pm 0.14$	$3.23\pm 0.09$
Coadd with systematic errors (w/o bin 7)	36,699	$-5.20\pm 0.11$	$3.18\pm 0.07$

Table 3. Table 3 Percent change in continuum in the forest due to changes in spectral index on the redside

Bin	$\overline{\alpha }$
Bin 1	-2.35	-3.02	-3.87	-2.91	-1.29	-2.19	-2.86
Bin 2	-2.37	-1.81	-3.47	-0.69	0.36	-0.82	NA
Bin 3	-1.75	-1.21	-1.40	1.62	4.55	2.52	1.23
Bin 4	-1.78	-0.08	0.22	3.08	4.72	1.64	2.69
Bin 5	-1.81	0.44	1.98	4.14	3.97	3.60	NA
Bin 6	-1.21	0.30	2.06	5.01	6.35	6.41	5.13
Bin 7	-1.25	0.00	3.00	6.94	8.58	6.21	NA

Table 4. Table 4 Previous measurements of the effective optical depth parameters

Work	# quasars	${\tau }_0$	$\gamma$
Kirkman+2005111Calculated using their errorbars on data, but using our model	24	$0.0049\pm 0.0011$	$3.03\pm 0.17$
Faucher+2008	86	0.0018	3.92
Bernardi+2003	1061	$0.0024\pm 0.0014$	$3.79\pm 0.18$
Becker+2013	6065	$0.0097\pm 0.0021$	$2.90\pm 0.12$
This Work	40,035	$0.0055\pm 0.0006$	$3.18\pm 0.07$

Equations

$$\bar{f}(z,{\lambda }_{\mathrm{rf}})=C({\lambda }_{\mathrm{rf}})\overline{T(z)},$$

$$1+z=\frac{{\lambda }_{\mathrm{obs}}}{1215.67\mathring{\mathrm{A}}}.$$

$${\lambda }_{\mathrm{rf}}(1+z_q)={\lambda }_{\mathrm{obs}},$$

$${\tau }_{\mathrm{eff}}(z)={\tau }_0(1+z{)}^{\gamma }.$$

$${\sigma }_T^2(z)=A{\left(\frac{1+z}{1+z_r}\right)}^B{\overline{T(z)}}^2,$$

$${\sigma }_f^2(z,{\lambda }_{\text{rf}})=A{\left(\frac{1+z}{1+z_r}\right)}^B{\bar{f}}^2(z,{\lambda }_{\text{rf}}).$$

$$x_0$$

$$x_1$$

$$\begin{gathered} \ln {\tau }_{\mathrm{eff}}=(-0.8563x_0+0.5165x_1-5.27)+ \\ (0.5165x_0+0.8563x_1+3.21)\ln (1+z). \end{gathered}$$

$$C^{\prime }=\left[\begin{array}{cc}{s}_0& 0\\ 0& s_1\end{array}\right]C\left[\begin{array}{cc}{s}_0& 0\\ 0& s_1\end{array}\right].$$

$$C_{\mathrm{sys}}[\ln {\tau }_0,\gamma ]=\left[\begin{array}{cc}0.106& -0.057\\ -0.057& 0.032\end{array}\right].$$

$${\tau }_0=(5.01\pm 0.76)\times 10^{-3},$$

$$\gamma =3.231\pm 0.086,$$

$$corr=\mathrm{0.99.}$$

$$C_{\mathrm{sys}}[\ln {\tau }_0,\gamma ]=\left[\begin{array}{cc}0.044& -0.029\\ -0.029& 0.020\end{array}\right].$$

$${\tau }_0=(5.54\pm 0.64)\times 10^{-3},$$

$$\gamma =3.182\pm 0.074,$$

$$corr=\mathrm{0.98.}$$

$$DM({\lambda }_{\mathrm{rf}})=0.01564$$

$$\displaystyle DM(\lambda_{\mathrm{obs}})=0.01686\$$

$$C_{\mathrm{sys}}[\ln {\tau }_0,\gamma ]=\left[\begin{array}{cc}0.114& -0.067\\ -0.067& 0.040\end{array}\right].$$

$$C_{\mathrm{sys}}[\ln {\tau }_0,\gamma ]=\left[\begin{array}{cc}1.39& -0.85\\ -0.85& 0.52\end{array}\right].$$

$$\Delta {\alpha }_i=\log ⟨\frac{f_i}{f_0}⟩/\log \left(\frac{⟨{\lambda }_{\mathrm{rf}}⟩}{1450\text{A˚}}\right),$$

$$m=\frac{d\log \alpha }{dz}.$$

$${\tau }_{\mathrm{eff}}(z)={\tau }_0{\left(\frac{1+z}{1+z_0}\right)}^{\gamma }+C$$

$${\tau }_{\mathrm{eff}}={\tau }_0(1+z{)}^{\gamma }+A\exp \left[\frac{-(z-z_{\mathrm{cen}}{)}^2}{2{\sigma }^2}\right];$$

Full text

\SIthousandsep

,

**affiliationtext: vikrant.kamble@utah.edu

1. Introduction

Following the formation of the first stars and quasars, the neutral Hydrogen gas present in the Intergalactic Medium (IGM) is ionized in a process known as the epoch of reionization (Zaroubi, 2013). The ionization state of the IGM is maintained in equilibrium through a balance of photoionization from background UV radiation and the adiabatic cooling of gas (Haardt & Madau, 2012). This ionization state evolves with redshift owing to the expansion of the Universe and the change in the rate of photoionization radiation from quasars and massive stars.

The strong absorption features imprinted in the spectra of high-redshift quasars (Gunn & Peterson, 1965) provide an important means to study the density, temperature and redshift evolution of the IGM following reionization (Haardt & Madau, 1996; Madau et al., 1999). These absorption features, arising from the resonant scattering of the background light from neutral Hydrogen gas with column densities $N_{\text{HI}}\sim 10^{14}\mathrm{cm}^{-2}$, are known as the Ly$\alpha$ forest (Lynds, 1971). Simulations that model the H I gas following a given temperature-density relation, a given background ionization intensity and perturbations following those of matter have been shown to reproduce the Ly$\alpha$ forest along with various properties remarkably well. One of the most basic quantities that can be obtained from Ly$\alpha$ forest absorption is the effective optical depth of the IGM to Ly$\alpha$ photons. The effective optical depth directly constrains the intensity of the background ionizing flux (McDonald & Miralda-Escudé, 2001; Bolton et al., 2005) and provides a link between the matter power spectrum and the flux power spectrum measured from Ly$\alpha$ forest (Croft et al., 1998; Seljak et al., 2003; Tytler et al., 2004; Palanque-Delabrouille et al., 2011; Chabanier et al., 2018).

Numerous studies using quasar spectra have found a smooth increase in the effective optical depth of the IGM with increasing redshift that can be described by a power-law evolution (Schneider et al., 1991). However, the estimates of the mean opacity at $z=0$, $\tau_{0}$, and the exponent of the power-law, $\gamma$, vary between measurements by more than the reported uncertainties. A possible source of these systematic errors is the estimation of the unabsorbed quasar continuum. High-resolution spectroscopic studies perform a direct fitting using the peaks in the Ly$\alpha$ forest to estimate the unabsorbed continuum (Schaye et al., 2003; Kim et al., 2007; Faucher-Giguère et al., 2008). These regions of each spectrum likely lead to underestimates of the continuum level, due to finite optical depth even in the most underdense regions. These systematic underestimates of the continuum level are expected to increase with redshift due to the increase in matter density.

Studies involving large samples of low signal-to-noise ratio ($S/N$) spectra have used composite spectra to obtain high $S/N$ representations of the average absorbed continuum. For this approach to work, it is important that the spectra being averaged have the same underlying unabsorbed continuum over all redshifts. This requirement introduces new complexities to the analysis as one must rely on sample selection based on features at wavelengths greater than the restframe Ly$\alpha$ emission and assume that these trends can be extrapolated into the Ly$\alpha$ forest. Moreover, since one does not know the underlying continuum, this method provides a direct estimate of only the relative optical depth. Previous studies have followed different approaches to model absolute optical depth measurements from low resolution spectra. Pâris et al. (2011) used Principal Components Analysis (PCA) to model the continuum at wavelengths longer than Ly$\alpha$ emission to predict the continuum in the Ly$\alpha$ forest. One can also perform a joint modeling of a parameterized continuum and optical depth evolution using the flux in the Ly$\alpha$ forest (Bernardi et al., 2003; Prochaska et al., 2009). Another approach is to perform measurements of the relative optical depth and to fit those measurements with a model for the evolving optical depth (Becker et al., 2013). Our modeling of the Ly$\alpha$ forest flux measurements is similar to those of Becker et al. (2013), except that we do not use composite spectra.

Bernardi et al. (2003) found a sharp bump at $z\sim 3.2$ of width $\Delta z\sim 0.4$ in the optical depth evolution. They attribute this feature to the increase in temperature of the IGM following He II reionization. A careful study using high $S/N$, high resolution quasar spectra by Faucher-Giguère et al. (2008) found a similar feature. However, other studies have not detected the bump (Kirkman et al., 2005; Pâris et al., 2011). Recent measurements using composite spectra created from a sample of 6065 quasar spectra (Becker et al., 2013) also reveal no such feature.

Given the size of the dataset available from the final Baryon Oscillation Spectroscopic Survey (BOSS; Dawson et al., 2013), we can explore the systematic errors that arise from modeling assumptions while also improving the statistical constraints on the redshift evolution of the mean Ly$\alpha$ transmission. Following a recent study on quasar spectral diversity by Jensen et al. (2016), we control the sample to have similar physical properties by dividing the spectra into seven bins based on observable parameters: spectral index and Carbon IV (C IV) equivalent width.

Assuming a smooth evolution of the effective optical depth modeled as a power-law, we constrain its parameters from each of these seven subsamples. The major differences from previous work are:

We bin by quasar properties and make seven independent measurements of the effective optical depth. This process allows a test of assumptions about uniformity of quasar continuum in the Ly$\alpha$ forest over the full redshift range. 2. 2.

We account for correlations between rest-frame wavelengths primarily due to cosmological fluctuations on small scales. 3. 3.

We quantify the systematic errors in the analysis by introducing systematic covariance matrices.

These refined measurements can supplement other areas of research that involve Ly$\alpha$ forest such as the measurement of Baryon Acoustic Oscillations (BAO) feature (Bautista et al., 2017; du Mas des Bourboux et al., 2017) and the one-dimensional power spectrum (Palanque-Delabrouille et al., 2013; Chabanier et al., 2018).

The structure of this paper is as follows: we describe our data and sample selection in Section 2. We also detail how the observable parameters were estimated and the methods used to correct for flux miscalibrations and measurement uncertainty. Section 3 describes the method used to obtain optical depth estimates from the raw flux data. Possible sources of systematic errors in our measurements are investigated in Section 4. The reconstructed continuum for each bin and their interpretation are presented in Section 5. Comparison of our results to optical depth measurements from previous studies are presented in Section LABEL:other_results, along with a discussion on the evidence of a He II feature. We summarize the analysis in Section 7.

2. Data and sample selection

The quasar spectra used to measure the effective optical depth were obtained from BOSS (Dawson et al., 2013), a part of the third generation of the Sloan Digital Sky Survey (SDSS-III; Eisenstein et al., 2011). The primary goal of BOSS was to extract cosmological constraints from BAO that behave like a standard ruler (Alam et al., 2017; Bautista et al., 2017; du Mas des Bourboux et al., 2017). This section describes how we identify an appropriate spectroscopic quasar sample for this study. We then characterize the spectral diversity of the sample and bin by common features. Finally, we perform corrections to the flux calibrations and uncertainties in the flux estimates.

2.1. Spectroscopic data

BOSS uses a pair of double spectrographs (Smee et al., 2013) mounted on the Apache Point 2.5 m Telescope (Gunn et al., 2006). The quasar selection involves a combination of algorithms (Kirkpatrick et al., 2011; Palanque-Delabrouille et al., 2011; Bovy et al., 2011) that are detailed in Ross et al. (2012). The observations are conducted using aluminum plates; each plate subtends an angle of $3^{o}$ on the sky and contains 1000 holes of $2^{\prime\prime}$ diameter drilled at locations corresponding to spectroscopic targets. Optical fibers are manually inserted into each plate to feed the pair of spectrographs. After observations of roughly one hour per plate, the raw data are processed into one-dimensional spectra and classified (Stoughton et al., 2002; Bolton et al., 2012). The reduced spectra used in this analysis correspond to version v5_10_0 of the pipeline presented in the 14th Data Release (DR14; Abolfathi et al., 2018). This version of the data reduction pipeline corrects for Atmospheric Differential Refraction (Margala et al., 2016; Jensen et al., 2016). This dataset is the first spectroscopic sample released publicly from eBOSS (Dawson et al., 2016), a component of SDSS-IV (Blanton et al., 2017).

We use the quasar classifications from the quasar catalog (DR12Q; Pâris et al., 2017) released with the 12th Data Release (DR12; Alam et al., 2015). We include the best spectrum of each quasar with ZWARNING=0. The redshift estimates obtained from visual inspection (Z_VI) were used as the systemic redshift of the quasars. Quasars with broad absorption lines (BALs) were removed using the BAL_FLAG_VI attribute specified in DR12Q. Quasars containing absorption from Damped Lyman-Alpha systems (DLAs) identified in the updated catalog from Noterdaeme et al. (2012) were also excluded. We select quasars with $z_{q}>1.6$ and a median $S/N$ per pixel greater than five computed over the bandpasses $1280<\lambda_{\mathrm{rf}}<1290$, $1320<\lambda_{\mathrm{rf}}<1330$, $1345<\lambda_{\mathrm{rf}}<1360$ and $1440<\lambda_{\mathrm{rf}}<1480$ Å in the quasar restframe. This cut on $S/N$ was made to obtain a reasonable precision in the placement of quasars into their appropriate bins based on the observable properties. The sample size after application of each criterion is reported in Table 1. The $S/N$ cut produces the largest fractional change in the sample size.

To correct for the effects of Galactic extinction, we adopt the Fitzpatrick model (Fitzpatrick, 1999) and the Galactic dust extinction map from Schlegel et al. (1998). Possible contamination from atmospheric emission were removed using a skyline mask from Delubac et al. (2015). Measured flux at $\lambda_{\mathrm{obs}}<3700$ Å was excluded because of lower $S/N$ and increased uncertainty in flux calibration. Each spectrum was shifted into its restframe wavelength solution in logarithmic units with the same pixel width as that of BOSS. The shifts were performed to the nearest pixel to avoid resampling.

2.2. Sample selection

It is crucial to minimize spectral diversity in the quasar sample so the redshift evolution of the observed flux levels in the Ly$\alpha$ forest can be attributed to the IGM and not to the quasars themselves. One prominent source of diversity originates from the Baldwin Effect, wherein luminosity is found to be anti-correlated with equivalent widths (EW) of broad emission lines (Baldwin, 1977). Jensen et al. (2016) found the spectroscopic signature associated with luminosity to be dominated by emission lines across the spectrum, a trend generally consistent with the Baldwin effect. It is likely that the variations in line strengths are correlated with variations in intrinsic quasar properties rather than simply luminosity. For example, Richards et al. (2011) argue that the C IV emission varies with quasar winds and the spectrum of the ionizing continuum. The color of the quasar defined by the spectral index is another feature of spectral diversity. E.g., Ivashchenko et al. (2014) report variations in spectral index lead to systematic errors in redshift estimates.

We introduce a new technique of measuring optical depth from low resolution spectroscopy compared previous studies (Bernardi et al., 2003; Becker et al., 2013) by exploring spectral diversity in the sample. This refinement is enabled by the sheer number of quasar spectra in the BOSS sample.

Motivated by the observations above, we assume that a significant fraction of the variations in the Ly$\alpha$ forest continuum arise from variations in spectral index ($\alpha_{\lambda}$) and C IV equivalent width $W_{\lambda}(\mathrm{CIV})$. Spectral index captures the continuum level and its shape in the Ly$\alpha$ forest relative to the shape of the continuum at wavelengths longer than 1216 Å. Variations in emission line features in the forest, particularly the strong O VI feature, are likely to be related to variations in C IV equivalent width. Hence to reduce the diversity within each sample, we first measure $\alpha_{\lambda}$ and $W_{\lambda}$(C IV) for each quasar so that the quasars can be binned in a grid spanned by those parameters.

To estimate the spectral index, $\alpha_{\lambda}$, a pseudo-continuum of the form $f=b\ \lambda^{\alpha_{\lambda}}$ is fit to each spectrum over the range $1280-1290,1320-1330,1345-1360$ and $1440-1480$ Å. These regions are relatively free of emission lines after visually inspecting the high $S/N$ composite spectrum from Harris et al. (2016). To ensure reliable estimates of spectral index, only spectra that contain more than 20 good pixels in each of the wavelength intervals were used. The fit was iterated three times, each time clipping points that were three standard deviations below the median estimate. This outlier rejection was performed to mitigate the bias arising from narrow metal absorption lines. We do not use restframe windows above C IV (1550 Å) emission to avoid contamination by iron emission line complexes. Unfortunately, the values of the spectral indices depend strongly on the restframe ranges used for the calculation. However, since our aim is to create a uniform sample across redshift for each bin, potential error in spectral index determination should not bias the results as long as the measurements do not vary with redshift. We explore this assumption further in Section 4.

For calculating the equivalent width of C IV emission, the line is modeled as a sum of two Gaussians. A simple estimation of the underlying continuum was performed using a linear fit over the range $1450-1465$ Å and $1685-1700$ Å. The continuum-subtracted spectrum was fit by a double Gaussian over the wavelength range $1500-1580$ Å. It is widely known that the broad and narrow components of the line arise from gas with different kinematics (Marziani et al., 2010); hence the parameters (location, scale and amplitude) of the two components are varied independently. We excluded pixels at wavelengths larger than 1580 Å to avoid contamination from He II (1640 Å) emission. A first stage of rejection of possible absorption lines was performed by smoothing with a box kernel of size 20 pixels and then removing pixels from the raw spectrum whose flux values were three standard deviations below the smoothed spectrum. Iterative clipping against the two component Gaussian model was then performed, rejecting negative outliers at the three standard deviation level.

Figure 1 shows the fit of power-law continuum and two component Gaussian for a sample quasar. The distribution of all 58,062 quasars that meet our $S/N$ criterion in the $\alpha_{\lambda}-$ $W_{\lambda}$(C IV) plane is shown in the top panel of Figure 2.

To remove the contribution from varying luminosity, the spectra were normalized using the median value calculated in the restframe wavelength window $1450-1460$ Å. The objects were then divided into seven bins in the $\alpha_{\lambda}-W_{\lambda}(\mathrm{CIV})$ plane. The bins were created using the central 84% of objects in the projected $\alpha_{\lambda}$ distribution and the central 84% of objects in the projected $W_{\lambda}$(C IV) distribution. Three bins spaced evenly in $\alpha_{\lambda}$ were created from the remaining sample. The first and third $\alpha_{\lambda}$ bin were divided in half according to $W_{\lambda}$(C IV), while the middle bin was divided into three equal parts. This binning scheme and the numbers assigned to each bin are illustrated in the top panel of Figure 2. The bottom panel of Figure 2 presents the redshift distribution of quasars for each bin. The seven bins contain a total of 40,035 quasars.

Objects in each bin are then warped by a power-law to reduce broadband continuum variations that may appear in the Ly$\alpha$ forest. The common spectral index to which they are warped was chosen to be the mean spectral index of all the objects in a given bin.

We then require that the $\alpha_{\lambda}-$ $W_{\lambda}$(C IV) distribution of objects be identical across quasar redshifts. For this purpose, we first discretized the quasar redshifts into intervals of width 0.2, with the first interval being $2.1<z_{q}<2.3$. Weights are assigned to each object as a function of $\alpha_{\lambda}$ and $W_{\lambda}$(C IV), such that the weighted probability distribution is identical across all redshift intervals. These weights are multiplied by the inverse variance vector for each object, thus reducing biases owing to selection effects. We test our assumptions of mitigating spectral diversity in each bin in Section 5.

2.3. Calibration

BOSS spectra have been shown to suffer from systematic errors in the estimates of measurement uncertainty assigned by the reduction pipeline (Palanque-Delabrouille et al., 2013; Delubac et al., 2015). To recalibrate the uncertainty estimates, we follow a similar analysis as was done in Lee et al. (2013). We identify a bandpass of width 10 Å in the restframe with relatively few emission lines or gradient in the continuum level. For each quasar spectrum, we perform a $\chi^{2}$ fit against an average flux level and rescale the pipeline uncertainty until we achieve a $\chi^{2}$ equal to the number of degrees-of-freedom. The sideband ranges chosen for this process are $1350-1360$ Å and $1470-1480$ Å. This ratio was averaged in the observer frame over all quasars with $1.6<z<4.0$. The fractional error in the pipeline uncertainty estimate, $\eta(\lambda_{\mathrm{obs}})$, as a function of observed wavelength, is displayed in the top panel of Figure 3. The pipeline overestimates the measurement error for $\lambda_{\mathrm{obs}}>4000$ Å  with a maximum difference of $20\%$ at $\lambda_{\mathrm{obs}}\approx 5800$ Å. We find an overall shape that is roughly consistent with that measured by (Palanque-Delabrouille et al., 2013) but with a positive offset. The difference could be caused by the new reduction algorithms that were introduced in DR14. Earlier studies also used a larger sideband range of 50 Å, potentially leading to additional dispersion caused by spectral features in the quasars. To examine this assumption, we repeat our analysis using the same 50 Å sidebands and find a suppression in $\eta(\lambda_{\mathrm{obs}})$, similar to previous studies. In this work, we apply corrections to the estimates of the flux uncertainty using the 10 Å sideband results.

Flux calibration in BOSS relies on theoretical models for F-type stars used as spectroscopic standards. Incorrect modeling of stellar features or Galactic absorption can distort estimates of the Ly$\alpha$ transmission in a redshift-dependent fashion. We assume that all quasars in a bin share the same spectral properties at wavelengths longer than 1216 Å and measure variations in the restframe spectrum as a function of observed wavelength. We attribute systematic difference at any observed wavelength to flux calibration errors. In correcting these errors, we use the same ranges that were employed to measure the spectral indices. We model the flux as a function of observed wavelength over each restframe wavelength pixel, normalized with the average flux measured over the observer wavelength range $4600-4640$ Å. The results are presented in the bottom panel of Figure 3. The flux calibration errors are in agreement with those presented in (Bautista et al., 2017; Lan et al., 2018). The flux calibration appears to deviate by few percent at wavelengths below $3700$ Å. In the subsequent analysis, we only use flux measurements at observed wavelengths greater than this value. It is important to note that we cannot remove large scale flux calibration errors, as they are degenerate with the power-law distortion applied to each quasar spectrum. We revisit the effects of absolute flux calibrations on the estimation of optical depth parameters in Section 4.

While we apply these corrections to measurement uncertainties and flux calibration, we are confident that they have no significant effect on the analysis. Hence we do not investigate their effects in detail.

3. Optical depth measurements

Optical depth studies have historically used high-resolution spectra with direct models for the continuum or relied on composite spectra built from large samples of low-resolution data (e.g. Bernardi et al., 2003; Becker et al., 2013). We constructed a large sample of low-resolution spectra but make use of the raw flux values at each pixel for each quasar rather than composite spectra. This allows one to model the variance contribution from Large Scale Structure into the likelihood definition. With the seven quasar bins available, we perform seven independent measurements of the redshift evolution of optical depth. Since the quasars in each bin illuminate the same global matter density field, any differences in the derived optical depth can be attributed to systematic errors associated with residual quasar diversity or with the parameter estimation.

3.1. General framework

The Ly$\alpha$ forest region used for this analysis covers the restframe wavelength range $1070\leq\lambda_{\mathrm{rf}}\leq 1160$ Å. In this region of each quasar spectrum, the average observed flux, $\bar{f}(z,\lambda_{\mathrm{rf}})$, at a given pixel is given by:

[TABLE]

where $C(\lambda_{\mathrm{rf}})$ represents the unabsorbed quasar continuum and $\overline{T(z)}$ is the mean transmission of the neutral Hydrogen in the IGM that gives rise to the Ly$\alpha$ forest. The mean transmission is related to the effective optical depth according to the relation $\overline{T(z)}=e^{-\tau_{\mathrm{eff}}(z)}$. Since the Ly$\alpha$ absorption takes place at the restframe of the absorber, the absorber redshift $z$ is related to the observed wavelength $\lambda_{\mathrm{obs}}$ as:

[TABLE]

However, since we are working in the restframe of the quasar, the restframe wavelength $\lambda_{\mathrm{rf}}$ is also related to the observed wavelength as:

[TABLE]

where $z_{q}$ is the quasar redshift.

The mean transmission is the product of Ly$\alpha$ and metal absorptions at different redshifts whose contributions are degenerate with those from neutral Hydrogen. In Subsection 3.3, we argue that metals do not significantly affect the measured redshift evolution of neutral Hydrogen density, but primarily change the estimation of the quasar continuum.

As in previous studies (Press et al., 1993; Bernardi et al., 2003; Kim et al., 2007), we parameterize the effective optical depth by a power-law:

[TABLE]

Each flux value, $f(z,\lambda_{\mathrm{rf}})$, from each quasar, is assumed to be drawn from a Gaussian distribution with mean $\bar{f}(z,\lambda_{\mathrm{rf}})$ and variance $\sigma^{2}=\sigma_{f}(z)^{2}+e_{i}^{2}$. The measurement error, $e_{i}$, is assigned after taking into consideration read noise, photon noise, and other processes and corrected as per Subsection 2.3. The redshift-dependent contribution from large-scale structure (LSS), $\sigma_{f}(z)$, is approximated to be Gaussian and hence added in quadrature to the measurement errors.

Following Lee et al. (2015), the variance in the transmission from Large Scale Structure is modeled as:

[TABLE]

where $z_{r}=2.25$, $A=0.065$, and $B=3.8$ as given in McDonald et al. (2000). However, since we are using measurements of flux, the variance in transmission is translated to variance in flux according to the relation:

[TABLE]

A joint likelihood fit was performed for all quasars in a given bin at a given restframe wavelength. Our model thus consists of three free parameters: the unabsorbed continuum $C(\lambda_{\text{rf}})$ and the two optical depth parameters $\tau_{0}$ and $\gamma$. The posterior distribution was marginalized over the parameter $C(\lambda_{\mathrm{rf}})$ to create the likelihood contours over $\tau_{0}$ and $\gamma$. This process was performed over all the restframe wavelength pixels in the forest range $1070\leq\lambda_{\mathrm{rf}}\leq 1160\ \mathrm{\AA}$. We use a more conservative range than BOSS BAO studies (Bautista et al., 2017; du Mas des Bourboux et al., 2017) to mitigate contamination from Ly$\alpha$ emission, O VI emission, and extrapolation of the power-law correction to short wavelengths. We also performed fitting allowing $A$ and $B$ to be free parameters. The best-fit values were consistent with McDonald et al. (2000) and did not significantly change the $\tau_{0}-\gamma$ contour, hence they were fixed as described above.

The final likelihood surface over the optical depth parameters for each bin of quasars was obtained by combining the contours over all 351 restframe wavelengths. Since there is a strong correlation between $\ln\tau_{0}$ and $\gamma$, it was computationally efficient to sum the likelihood in a rotated basis where the transformed parameters are approximately orthogonal. The two IGM parameters in this new basis are represented as ($x_{0},x_{1}$) which are related to ($\ln\tau_{0},\gamma$) as:

[TABLE]

The effective optical depth in terms of $x_{0}$ and $x_{1}$ is given as:

[TABLE]

A high degree of correlation exists between measurements of both $x_{0}$ and $x_{1}$ amongst neighboring pixels. The neighboring pixels are separated by less than one Mpc along the line of sight. Hence, we expect LSS correlations to appear in the statistics of the optical depth estimates from one restframe wavelength to another. We also attribute this correlation to the resampling of pixels and the resolution of the instrument.

To account for the effects of these processes on the statistical significance of the IGM parameter estimates, we computed the mean correlation function as a function of pixel separation for each bin. To isolate this component and ensure a positive-definite covariance matrix, each correlation function was modeled as an exponential fit to the first five data points. This exponential function is a decent fit to each mean correlation function with a typical scale length $r_{0}=2.24$ pixels along the $x_{0}$ direction and $r_{0}=3.4$ pixels along the $x_{1}$ direction. We converted the correlation function into a correlation matrix by assigning each band diagonal element the corresponding value from the analytic fit to the correlation function. The correlation matrix was transformed into a covariance matrix by weighting each element by the variance estimates defined earlier.

Incorporating the effects of covariance into the final likelihood surface is complicated by the fact that the likelihood surfaces for each restframe wavelength are not Gaussian. When ignoring the correlation between pixels and approximating the likelihood surface for each restframe wavelength as an uncorrelated two-dimensional Gaussian, we find that the confidence intervals are preserved in the combined likelihood surface when compared to a fit using the original likelihoods. However, biased estimates of the central values of $x_{0}$ and $x_{1}$ are produced. To account for the correlations and the fact that the contours for each restframe wavelength are not Gaussian, we adopted the following procedure: we first defined a stretch factor as the amount by which the uncertainty estimates change when taking the small scale correlations into account, as compared to a simple weighted estimate without correlations. This calculation was done independently along the $x_{0}$ and $x_{1}$ direction, yielding stretch factors $s_{0}$ and $s1$, respectively. We incorporate the effect of small scale correlations into the covariance matrix of the full likelihood estimate of $x_{0}$ and $x_{1}$ without correlations by rescaling as follows:

[TABLE]

The addition of cosmological correlations to the analysis has the impact of increasing the uncertainty in the $x_{0}$ estimates by an average factor of $s_{0}=2.04$. The uncertainty in $x_{1}$ estimates is increased by an average factor of $x_{1}=2.23$. The correlation between neighboring pixels effectively reduces the number of independent restframe wavelengths by 75%.

3.2. Combined estimates

Figure 4 presents the optical depth parameter estimates for all seven bins. The left panel displays the two-dimensional contours for each quasar bin in the orthogonal basis. The contours for the combined statistical likelihood from these seven measurements are shown in magenta in the right panel. The number of quasars in each bin, the best-fit estimates of $\ln\tau_{0}$ and $\gamma$, and their statistical errors are listed in Table 2. The constraints on the parameters arising from statistical errors for each bin are similar, ranging from 2.5% to 5.9% on $\ln\tau_{0}$ and 2.5% to 6.3% on $\gamma$. However, the central values deviate by more than what would be expected from the reported statistical errors, as seen in the left panel of Figure 4. The fifth column of Table 2 reports the probability that each individual measurement would produce the central value resulting from the combined likelihood. Four of the seven bins have a p-value less than one percent.

We attribute the extra scatter in the optical depth parameters in each bin to systematic errors not accounted for in the analysis. We model this systematic component using a nuisance covariance matrix that is added to the statistical covariance matrix of each data point in the modified basis. The seven measurements are fit with a model that includes additional free parameters describing the systematic error on $x_{0}$, the systematic error on $x_{1}$, and their correlation.

After rotating back to the $\tau_{0}-\gamma$ basis, the most likely values of the systematic covariance matrix assigned to each bin are estimated to be:

[TABLE]

We marginalize over these nuisance parameters to obtain the maximum likelihood estimates for the optical depth parameters. The net effect of the systematic errors is to increase the uncertainty on $\ln\tau_{0}$ by a factor of 2.5 and on $\gamma$ by a factor of 2.3. These estimates, transformed back to the $\tau_{0}-\gamma$ basis are:

[TABLE]

We then recompute the p-value for each of the seven measurements using the statistical+systematic covariance matrices. These p-values are listed in the last column of Table 2. Given the sample size none of the measurements are significantly deviant.

The 68% and 95% contours for the full two-dimensional likelihood are presented as dashed black curves in the left panel of Figure 4. These contours overlap the 68% statistical contours for each of bin 1 through bin 6, but do not overlap the 95% contours in bin 7. The central value of bin 7 appears to be driving the systematic error on $x_{1}$ toward larger values. To determine the extent of this effect, we repeat the analysis excluding the data from bin 7. In this case, the systematic covariance matrix assigned to each bin is:

[TABLE]

The variance in the systematic contribution to $\ln\tau_{0}$ is decreased by a factor of 2.4 after removing bin 7 and the variance in the systematic contribution to $\gamma$ is decreased by a factor of 1.6. The reduction in the size of the systematic errors alone is not sufficient evidence to remove bin 7 from the analysis; however, we provide additional evidence in Section 4 and Section 5 in support of removing this bin from the analysis. Excluding bin 7, the estimates of the optical depth parameters become:

[TABLE]

The uncertainty estimates on both $\tau_{0}$ and $\gamma$ are reduced by approximately 15% after removing bin 7. The central values shift slightly towards shallower redshift evolution and higher optical depth at $z=0$. The likelihood contours after modeling the systematic errors and excluding bin 7 are shown in solid black in the $x_{0}-x_{1}$ basis and the $\ln\tau_{0}-\gamma$ basis in Figure 4.

3.3. Metal contamination

The mean transmission is the product of absorption from metals and Ly$\alpha$. This effect will be a function of both observer wavelength owing to the evolution of optical depth and restframe wavelength as more metal lines contribute to the absorption at shorter wavelengths. Given the low resolution of the BOSS spectra, it is impossible to correct metal absorption in each spectrum individually; instead, we assess its statistical contribution. We adopt the corrections obtained by Kirkman et al. (2005) determined from 52 high resolution quasar spectra published in Sargent et al. (1988). They model the metal absorption as $DM=1-T=1-e^{-\tau_{\mathrm{eff,m}}}$:

[TABLE]

For the evolution of metal transmission over the observed wavelength range, we evaluate the range of transmission values at a restframe wavelength $\lambda_{\mathrm{rf}}=1115$ Å. At this restframe wavelength $DM(\lambda_{\mathrm{obs}}=3700\ \text{\AA})=0.028$ and $DM(\lambda_{\mathrm{obs}}=7000\ \text{\AA})=0.023$, corresponding to the extremes in the redshift range of our study. The evolution only amounts to 0.5%, sufficiently small to be neglected. This contribution to absorption from metals will lead to an underestimation of the unabsorbed continuum, $C(\lambda_{\mathrm{rf}})$ by $2.5\%$. When reconstructing the continuum in each bin using the optical depth parameter estimates, we correct for this contribution from metals by scaling the estimated value of the unabsorbed continuum by the inverse of the average transmission of the metals.

4. Systematic errors

We make two core assumptions in our analysis: that spectral index and equivalent width sufficiently capture the spectral diversity and that the continuum can be standardized across the entire sample of quasars in each bin. Here we investigate these assumptions.

4.1. Different basis

Since our measurements of optical depth parameters are dominated by systematic errors, we explored other ways of binning the sample. We choose two different parameter spaces to bin the quasar population and assess the size of systematic errors compared to those found in Equation 10.

We first binned on spectral index and C IV FWHM because the latter was shown to be highly correlated with quasar diversity in Jensen et al. (2016). The total sample of 42,615 quasars covers seven unique bins divided by the same percentiles in spectral index and C IV FWHM as was done in Section 2. The FWHM was determined from the best-fit double Gaussian model as described in Section 2.2. Constructing composite spectra revealed a redshift evolution in the equivalent widths that is likely due to Malmquist bias. The equivalent widths of emission lines in the forest are small (Harris et al., 2016), so this variation is not expected to be significant. We follow the procedure presented in Section 3.1 and Section 3.2 to estimate $\tau_{0}$ and $\gamma$ and their associated systematic errors. The best-fit optical depth parameters lie 2.8 standard deviations from those found in Equation 3.2, indicating a shallower evolution in $\tau_{\mathrm{eff}}$. The systematic error covariance matrix is:

[TABLE]

The variance between estimates of $\gamma$ due to systematics errors in this basis are roughly 25% higher than in the $\alpha_{\lambda}-$ $W_{\lambda}$(C IV) basis used for the main results of this work.

We next binned on Eddington ratio and black-hole mass based on single-epoch spectroscopic scaling relationships. We used the relationship between black-hole mass, luminosity and C IV FWHM presented in Shen et al. (2011). Using this basis is an attempt to divide the sample by the physical parameters of the quasars. Although we do not bin by spectral index, we still perform a power-law correction to the continuum as in Section 2.2. The total sample of 42,439 quasars covers seven unique bins covering approximately the same fraction of the parameter space as in Section 2. We again follow the procedure presented in Section 3.1 and Section 3.2 to estimate $\tau_{0}$ and $\gamma$ and their associated systematic errors. The best-fit optical depth parameter $\tau_{0}$ lies 2.5 standard deviations below and $\gamma$ lies 3.9 standard deviations above from those found in Equation 3.2. The systematic error covariance matrix is:

[TABLE]

The variance between measurements due to systematics errors in this basis are so much larger than in the $\alpha_{\lambda}-$ $W_{\lambda}$(C IV) basis, that we infer a flaw in the assumptions of diversity based on estimates of Eddington ratio and black hole mass. Another indication for the problem in this basis is the low number of quasars above redshift $z=3$ in several of the bins.

This exercise highlights that care is needed in identifying subsamples to have accurate optical depth measurements. The test using Eddington ratio and black hole mass to isolate diversity demonstrates that poor coverage at high redshifts leads to a model that prefers higher values of $\gamma$, at least in the case of the BOSS spectroscopic sample. Indeed, a similar trend can be seen in bins 1, 2, 5 and 7 in the main analysis. These four bins are each roughly half the sample size of the other bins. There is also a slightly lower representation at higher redshifts; 10% of the quasar sample in bin 3 lies at $z>3$, while only 5% of bins 2 and 7 lies at these higher redshifts.

4.2. Bias in spectral index measurements

One possible source of error neglected in our analysis is the uncertainty in the values of spectral indices. Variations in the spectral index imply variations amongst the pixels in the Ly$\alpha$ forest that are correlated for a given quasar. Incorporating this correlated noise is complicated because data from individual quasars are modeled independently across restframe wavelengths. Here, we use a MonteCarlo approach to assess the contribution of statistical and systematic errors on spectral index to the systematic error matrix of our measurements of optical depth parameters.

The typical uncertainty in spectral index measurement for a given quasar is $\sigma_{\alpha_{\lambda}}\sim 0.1$. Extrapolating this uncertainty from the region used to measure the spectral index to shorter wavelengths results in an uncertainty of $\sim 3\%$ in the continuum level in the Ly$\alpha$ forest. This quantity is small compared to the dispersion introduced by the LSS term and the typical measurement uncertainty at each pixel; however, this term doesn’t decrease as $\sqrt{n}$ when combining the measurements across restframe wavelengths due to its correlated nature.

To estimate the effect of correlated continuum errors across the whole sample, we use simulated power-law quasar spectra created with uncertainty in the spectral index measurements. To each spectral index, we add an extra term that is sampled from a normal distribution with zero mean and standard deviation of 0.1. We then measure the best-fit $\ln\tau_{0}$ and $\gamma$ values for a sample of 7000 simulated spectra. This process is repeated on 100 unique realizations. The scatter in the central values of $\ln\tau_{0}$ and $\gamma$ captures the covariance associated with continuum errors. The covariances account for only 6% of the systematic covariance matrix quoted in Equation 10.

We now investigate the bias in spectral index measurement as a possible source of systematic error. In Section 2.2, we chose the relatively line-free restframe wavelength region between 1280 and 1480 Å to estimate the spectral index. After warping all spectra to follow the same continuum over this wavelength range, the composite spectra divided by redshift in each quasar bin, shown in Figure 5, reveal residual spectral diversity at restframe wavelengths longer than 1600 Å.

To gauge the possible effect of this diversity in the Ly$\alpha$ forest continuum, we model the variations at longer wavelength according to a systematic error in the spectral index estimation. The difference in spectral index for a composite spectrum in a given redshift interval for a given bin with respect to the first redshift interval is computed as:

[TABLE]

where $f_{0}$ and $f_{i}$ are the flux vectors for the first redshift interval and the $i^{th}$ redshift interval, respectively. The mean flux ratio, $\left<\frac{f_{i}}{f_{0}}\right>$, is computed using these flux vectors over the restframe wavelength range $1600-1800$ Å and $\left<\lambda_{\mathrm{rf}}\right>=1700\ \text{\AA}$. The relative change in the Ly$\alpha$ forest continuum compared to the composite spectrum in the first bin is approximated by evaluating the quantity in Equation 16 at restframe wavelength of 1100 Å.

Table 3 presents the percent change in the Ly$\alpha$ forest continuum values if these redshift-dependent variations at 1700 Å were present at 1100 Å. To predict the systematic errors on $\tau_{0}$ and $\gamma$, we simulate composite spectra following the inverse variance and composite redshifts of bin 3, using our central values of $\ln\tau_{0}$ and $\gamma$. The continuum was perturbed at each redshift interval following the measured changes quoted in the table for bin 3. The mean optical depth parameters across all the restframe wavelengths were then measured on these simulated composite spectra. The deviations of the measured values from the input values are much larger than the systematic errors given in Equation 10. These results indicate that the variations observed at 1700 Å cannot be used to predict the absolute variations in the Ly$\alpha$ forest continuum. However, this result does not imply the extrapolated continuum is free of redshift dependent biases.

To quantify the sensitivity of the optical depth model to redshift dependent errors in the spectral index measurements, we investigate the effect under a set of simple assumptions. We assume that the continuum for the first redshift bin ($z_{q}\sim 2.1$) is correctly modeled, and assume the systematic offset in spectral index measurements to be linearly increasing with redshift. This variation is parameterized by a slope parameter, $m$, that represents the fractional change in spectral index per unit redshift:

[TABLE]

As before, we measure the mean optical depth parameters across all the restframe wavelengths on these simulated composite spectra. The offsets in $\tau_{0}$ and $\gamma$ are compared to the input values as a monotonic function of $m$. As the slope increases, the central value of $\tau_{0}$ decreases and vice versa for $\gamma$. We focus on bin 7 because it is the largest outlier in optical depth parameters in Table 2 and in continuum offsets in Table 3. The central value of $\ln\tau_{0}$ measured in bin 7 differs from the global central value by 0.75 and the central value of $\gamma$ measured in bin 7 differs from the global central value by 0.36. A fractional change in the spectral index per unit redshift as large as 7.6% is required to explain this shift in the optical depth parameters.

5. Interpretation

In this section, we consider the measurements of the optical depth parameters presented in Section 3 in the context of predictions for the unabsorbed continuum. We examine the diversity of features in the Ly$\alpha$ forest of the quasar spectra. We then present the measurements of the mean transmission from the observer frame to test how accurately we predict the underlying continuum.

5.1. Reconstructed continuum

We reconstruct the Ly$\alpha$ forest continuum for each quasar spectrum using the best-fit optical depth parameters from its respective bin in Table 2; then apply the same spectral warping to each spectrum as was done in Section 3.1. A composite spectrum for each bin is created from these unabsorbed spectra. We correct for the contribution of metals as discussed in Section 3.3. The resulting composite spectra covering the forest range $1050-1170$ Å are shown in Figure 6.

The top panel shows the reconstructed continuum normalized by a power-law extrapolation using the mean spectral index for each bin. If the estimates of spectral index were unbiased and the same power-law describes the continuum in the Ly$\alpha$ forest, then the continuum level for all composite spectra would be unity. Instead, using 1100 Å as a proxy for the line-free continuum, the reconstructed spectra deviate over a range of -10% to +5%. The deviations from bins 1, 2, 5 and 7 indicate an underestimate of the continuum level. As stated in Section 4.1, these bins have less representation at higher redshifts and are best described with a steeper evolution of optical depth. If the suppression of estimated continuum level is consistent across all redshifts, then no systematic bias in optical depth parameters should arise in these bins. Conversely, it is possible that the variations in predicted continuum levels are an artifact of redshift or signal-to-noise ratio dependent errors in the spectral index. Further investigation with larger, deeper samples is required to disentangle these possible causes.

The bottom panel of Figure 6 presents spectral diversity amongst the bins using bin 4 as a reference. Doing so highlights the diversity in line features. The reconstructed composite spectrum from bin 4 is warped using a powerlaw and fit to each of the other six bins. The reconstructed continuum in each bin is shown in color following the same scheme as the top panel, relative to the model based on bin 4 in black. The warped model matches remarkably well with the reconstructed continuum for the other bins. The maximum difference occurs near the Fe III blend centered around 1123 Å. For bins 3 and 6, the differences are only 2.2% and 2.7%, respectively, indicating a mild negative correlation with the strength of the C IV line. These results indicate that the assumptions of a uniform model for the features in the Ly$\alpha$ forest continuum made in BOSS cosmology studies (Bautista et al., 2017; du Mas des Bourboux et al., 2017) are reasonable.

5.2. Mean Transmission

We now revisit the analysis from the perspective of the observer frame. Our goal is to assess how well the direct estimates of mean transmission follow the power-law model for the evolution of optical depth and how the data from each bin compare to predictions from the global model.

The model in Section 3.1 and 3.2 leads to degeneracy between the optical depth parameters $\tau_{0}$ and $\gamma$, and the estimate of the unabsorbed continuum at each restframe wavelength. The work presented in Section 5.1, however, uses the joint estimates of $\tau_{0}$ and $\gamma$ from all the restframe wavelengths to provide a high-precision estimate of the reconstructed continuum for each bin. Using this continuum as a model, we measure the transmission directly as a function of redshift for each quasar in the Ly$\alpha$ forest region over the observed wavelengths. The left panels of Figure 7 show the mean of these transmission values as a function of redshift, in bins of size $\Delta z=0.05$.

The models using the best-fit parameters for each bin are indicated in red in Figure 7. Using 10,000 bootstrap resamplings to assess the covariance between data-points, we compare each model to its respective data sample. We find a typical $\chi^{2}\sim 140$ for a total of 40 data-points, indicating a discrepancy between the data and the power-law model. The break from a power-law is evident in bins 1 and 2 in the form of a wiggle with a maximum excursion around $z=2.5$. However, the signature of the excursion changes across all the bins. The variable nature of the excursions indicate that the high $\chi^{2}$ values are not likely due to flux calibration errors or a consistent failure in the power-law model to describe the data. We were unable to identify the true source of these deviations but hypothesize that residual quasar diversity across the $2<z<4$ redshift range is responsible. For example, the inconsistencies at restframe wavelength 1700 Å across quasar redshifts, as seen in Figure 5, could produce such features in the measured transmission if indeed they are present at wavelengths shorter than 1200 Å.

The right panels in Figure 7 present the difference between the measured transmission and the global best fit model as given by the parameters in Equation 3.2. For bins 1 through 6, the average deviations range from 0.2% to 2.4%, compared to a typical uncertainty of 1.8% on the modeled transmission. Bin 7 deviates by an average of 11.5%; indicating that the continuum is systematically underestimated when using the optical depth parameters derived from that sample of quasars ($\ln\tau_{0}=-6.06,\gamma=3.59$). Bin 7 was also the largest outlier in optical depth parameter estimates (Figure 4) and showed the largest trend in the continuum residuals at restframe wavelengths around 1700 Å (Figure 5). This overall outlier behavior of the quasar spectra in Bin 7 is likely a result of systematic errors in the estimates of the spectral index or uncontrolled diversity in the sample. Our final constraints on the optical depth parameters henceforth (and in the abstract) are reported using only bins 1 through 6.

6. Comparison to Previous results

This section compares our results to those from previous studies. We include works from both high-resolution and low-resolution spectroscopy in our comparison. We first compare the smooth evolution of optical depth parameterized by a power-law, then discuss whether the BOSS data provide evidence for the previously-reported He II feature.

6.1. Optical depth evolution

Measurements of the effective Ly$\alpha$ forest optical depth in the literature fall into two camps. High-resolution, high signal-to-noise ratio spectra allow for a direct but approximate continuum fitting using the peaks in the forest region. This method suffers from a potential underestimation of the continuum, especially at higher redshifts. However, the approach allows one to directly measure the transmission probability distribution function and the associated statistics of the transmission field. Studies such as this work that rely on low-resolution, low $S/N$ spectra typically contain more objects and therefore a potential for higher precision statistical estimates. We present a summary of our optical depth measurements relative to four other works in Figure 8 and Table 4.

Kirkman et al. (2005) obtained a sample of 24 high resolution spectra from HIRES spectrograph (Vogt et al., 1994) on the Keck Telescope. They measured the evolution of transmission over the redshift range $2.2<z<3.2$, after removing contributions from Lyman Limit systems (LLS) and metal lines, and attempted to remove biases due to continuum fitting by using artificial quasar spectra. Their model differs from that used in this work, in that they assume a power-law in the quantity $DA=1-e^{-\tau_{\mathrm{eff}}}$. Table 4 presents constraints on the optical depth parameters produced by the model used in this paper, assuming their measurements shown in Figure 8 to be independent.

A detailed analysis using a larger sample of high-resolution quasar spectra was performed by Faucher-Giguère et al. (2008). They measured the effective optical depth over the redshift range $2\leq z\leq 4$ from a sample of 86 high-resolution, high-$S/N$ quasar spectra obtained with the ESI (Sheinis et al., 2002) and HIRES spectrographs on Keck and with the MIKE spectrograph (Bernstein et al., 2003) on Magellan. Their study contained a detailed analysis of systematic errors from continuum fitting and metal corrections. Fitting the data with a two parameter power-law produced a $\chi^{2}=40.1$ for 21 degrees-of-freedom. Adding more free parameters in the form of a three-component Gaussian, they improve the fit to a $\chi^{2}=26.8$, providing positive evidence for the additional free parameters according to the Bayesian Inference Criterion (BIC). The evidence for an expanded functional form could be attributed to He II reionization, a general deviation from a power-law, or systematic errors in the approach.

Bernardi et al. (2003) employed 1061 low-resolution SDSS quasar spectra to map the redshift evolution of optical depth over the Ly$\alpha$ redshift range $2.5\leq z\leq 4$. They constructed composite spectra in bins of width $\Delta z_{q}=0.2$ after normalizing each spectrum by its flux in the restframe wavelength range $1450-1470$ Å. They adopted a parametric form to describe the continuum in the forest consisting of a power-law and Gaussian features to account for emission lines at 1073 Å, 1123 Å and Ly$\alpha$ 1215.67 Å. In addition to a smooth power-law evolution, they reported a clear ‘bump’ at $z\sim 3.2$ that they attribute to He II reionization. Their best-fit model, including the modeled reionization feature, is shown as the brown dashed line in Figure 8.

A study that closely resembles our work is that of Becker et al. (2013), who used 6065 low resolution quasar spectra over the quasar redshift range $2\leq z\leq 5$ from the SDSS DR7 quasar catalog (Schneider et al., 2010). They constructed composite spectra in bins of quasar redshifts with a typical $\Delta z_{q}=0.1$, after correcting for differences in the spectral indices. Using the normalized transmitted flux measurements in bins of $\Delta z=0.1$, they fit the optical depth parameters simultaneously across all restframe wavelengths. To obtain absolute measurements of transmission, they scale their results to those from Faucher-Giguère et al. (2008) at $z\sim 2.35$ using an additional free parameter. This approach leads to an additional offset parameter in the effective optical depth parameterization given as:

[TABLE]

They estimate the value of this offset to be $C=-0.132$. For comparison, our most discrepant measurement in bin 7 would require a value of $C=0.097$ to bring the measured optical depth into alignment with the global model. Their best fit three-component model for the effective optical depth is indicated in black with 1$\sigma$ confidence intervals in Figure 8. They reported no evidence for a He II reionization feature.

Our best-fit model for the effective optical depth including the $\sigma$ confidence interval is shown in red in Figure 8. Our model predicts a larger typical opacity due to Ly$\alpha$ at $z<3.3$ than the measurements of all the other works. Over this redshift range, the average discrepancy ranges from $\Delta\tau_{\mathrm{eff}}=0.01$ in case of Bernardi et al. (2003) to $\Delta\tau_{\mathrm{eff}}=0.09$ in the case of Kirkman et al. (2005). We predict lower opacity than Bernardi et al. (2003) by an average of $\Delta\tau_{\mathrm{eff}}=0.05$ over the redshift range $3.3<z<4.2$. At these higher redshifts, we agree with the measurements of Faucher-Giguère et al. (2008) and Becker et al. (2013) at a level better than $\Delta\tau_{\mathrm{eff}}=0.04$.

A more quantitative comparison of the best-fit models is presented in Table 4. At face value, our constraints are similar to those of Kirkman et al. (2005), in that the $\tau_{0}$ and $\gamma$ estimates are consistent to better than one standard deviation. Likewise, our measurements lie within two standard deviations of the Becker et al. (2013) $\tau_{0}$ and $\gamma$ estimates. However, simply comparing the marginalized estimates of each parameter independently neglects the high degree of correlation between the parameters, as shown in the right panel of Figure 8. There is a large disagreement between all of the studies and our own. The measurement of Becker et al. (2013) and Bernardi et al. (2003) are in closest agreement, although they lie right at the edge of our 8$\sigma$ contour. The prediction of steeper evolution compared to our work as found in Faucher-Giguère et al. (2008) and Bernardi et al. (2003) indicates that their estimates lie in the extremes of our $\tau_{0}-\gamma$ degeneracy curve.

6.2. He II feature

A few previous studies (Bernardi et al., 2003; Faucher-Giguère et al., 2008) report detection of a narrow feature superposed on a smooth power-law evolution of the effective optical depth. This feature appears as a decrement with width $\Delta z\approx 0.1$ at $z\approx 3.2$. This feature is commonly attributed to He II reionization. Other studies do not detect such a feature, leaving doubt as to whether it arises from astrophysical sources or systematic errors.

To test for the presence of this He II feature, we use the mean transmission estimates as a function of redshift, and their covariances, from Section 5.2. As in previous studies, the effective optical depth is fit with a model of the form:

[TABLE]

to the measured transmission. We searched for evidence of a Gaussian feature against a simple power-law model using the the BIC, with 40 data-points. There are two free parameters in the smooth model and five free parameters in the model with Gaussian departure.

Figure 9 shows the $\Delta\chi^{2}$ surface for each of the seven bins, as a function of the location of the He II Gaussian bump. This $\Delta\chi^{2}$ is reported as the difference in $\chi^{2}$ between a five-component model and a two-component model. The $\Delta\mathrm{BIC}$ values shown are computed at the best-fit redshift of the Gaussian feature. Bins 1 and 2 produce positive $\Delta$BIC estimates at a redshift $z_{\mathrm{cen}}\approx 2.4$. None of the other bins produce a positive $\Delta$BIC nor do they indicate a local minimum at this redshift. The bottom panel displays the joint model to the 280 data-points that span all seven bins. We do see a global minimum in the $\chi^{2}$ surface at $z_{\mathrm{cen}}\approx 2.92$; however, the negative value of the $\Delta$BIC suggests that no meaningful information is provided by the additional free parameters.

7. Conclusions

This work has produced the tightest constraints to date on the redshift evolution of the mean effective optical depth due to Ly$\alpha$ absorption by neutral hydrogen. The evolution is fully-described by a powerlaw with no convincing evidence for He II reionization in the redshift interval $2.0<z<3.5$. The measurement of the powerlaw exponents differ by 0.59 between the extreme values over seven independent analyses and are discrepant by 3.3$\sigma$ in the two measurements that show the largest tension. The final measurements on the optical depth parameters yield $\tau_{0}=0.00554\pm 0.00064$ and $\gamma=3.182\pm 0.074$ after excluding a single deviant measurement and combining the results produced by the other six datasets.

One can compare the high-redshift results of this work to those at low redshifts based on Hubble Space Telescope (HST) observations of active galactic nuclei (AGN). Danforth et al. (2016) measured the average Ly$\alpha$ flux decrement at redshifts $0<z<0.4$ using a sample of 82 medium resolution spectra of UV-bright AGN obtained from the Cosmic Origins Spectrograph. They report best-fit optical depth model with $\tau_{0}=0.014\pm 0.001$ and $\gamma=2.2\pm 0.2$. The shallower evolution is consistent with our model if we allow a break in the power-law at a redshift around $z=1.6$. Future measurements of optical depth evolution over the redshift range $1<z<2$ will allow a direct test for such a break.

Systematic errors in the analysis appear to dominate the final measurement uncertainty. The uncertainty on the powerlaw exponent in the combined analysis is diluted from a statistical uncertainty of 1.2% to 2.3% when incorporating these systematic errors. Even after including these systematic errors, the measurement produces a precision that is a factor of 1.6 better than the previous measurement using the SDSS quasar sample (Becker et al., 2013). The improvement is enabled by the larger sample of quasar spectra from the BOSS program. Beyond the increased sample size, the analysis is improved over previous work by the inclusion of covariances due to Large Scale Structure and the division by quasar diversity. The former of these effects has been neglected in prior work but reduces the effective number of independent SDSS/BOSS wavelength bins by a factor of four. The latter allows a characterization of the source of systematic errors in the analysis.

Investigation of the systematic errors indicates that the measurements of spectral index used to normalize the Ly$\alpha$ forest continuum are subject to scatter and biased estimates. The algorithm in this analysis relies on a powerlaw fit to only 182 pixels over a wavelength range $1280<\lambda_{\rm rf}<1480$ Å. Redshift-dependent variation in N V emission (1240 Å), the O IV/Si IV emission line complex (1400 Å), or C IV (1549 Å) emission may lead to contamination in the region used to estimate the continuum and thus explain part of the bias. Likewise, redshift dependence on the signal-to-noise ratio or some other affect of small statistical size of the continuum region may be biasing the powerlaw estimates. A future analysis may be able to mitigate these errors by taking a more comprehensive approach to continuum estimation and normalization. For example, archetype spectra can be identified based on a Set Cover Problem (e.g. Zhu, 2016), and a controlled sample of quasars can be identified around each archetype based on a $\chi^{2}$ nearest neighbor determination in the unabsorbed continuum at wavelengths longer than 1216 Å. The spectral warping that we apply based on the simple continuum estimates could instead be performed using the entire unabsorbed wavelength region, thus producing a higher precision estimate that naturally incorporates emission line diversity.

Finally, we do not assess systematic errors due to flux calibration in this study, as such an analysis is non-trivial given the lack of an independent broadband reference. The median residual when comparing eBOSS spectrophotometric fluxes to imaging fluxes in ($g,r,i,r-i$) have been shown to be ($-0.001,0.004,-0.022,0.032$) magnitudes, respectively (Jensen et al., 2016). These systematic biases are at the level of the final precision in our analysis of the redshift evolving mean transmission. In addition, the lack of any strong deviation in the mean transmission from a powerlaw optical depth model provides evidence against significant spectroscopic calibration errors. The flux calibration errors would have to follow a powerlaw to avoid detection in that analysis.

The final sample from eBOSS exceeds 200,000 high-redshift quasar spectra. Such a sample can improve our results or extend the study by incorporating the Ly$\beta$ forest covering the restframe wavelength range $978-1014$ Å (e.g. Iršič et al., 2013). Further improvement with the eBOSS spectra would require a reduction in the systematic errors to make use of the larger sample. An even larger sample of quasar spectra, with more stable flux calibration and higher $S/N$ ratio will be produced by the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration et al., 2016a, b) over the time period 2020–2025. DESI will cover 14,000 square degrees and produce a sample of 700,000 quasar spectra at redshifts $z>2.1$. Observations of these quasars will have an effective exposure time roughly four times higher than that of eBOSS, thus leading to spectra that should be less susceptible to biases in continuum estimation arising from low $S/N$ ratio data. This DESI sample is designed for 1% precision constraints on the Hubble parameter from BAO and will be very well suited for a new measurement of the mean transmission in the Ly$\alpha$ forest.

Acknowledgments: This paper represents an effort by both the SDSS-III and SDSS-IV collaborations. Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III web site is http://www.sdss3.org/.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org.

SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

The work of VK, KD and HdMdB was supported in part by U.S. Department of Energy, Office of Science, Office of High Energy Physics, under Award Number DESC0009959.