Issue 
A&A
Volume 603, July 2017



Article Number  A12  
Number of page(s)  23  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201730533  
Published online  30 June 2017 
Measurement of baryon acoustic oscillation correlations at z = 2.3 with SDSS DR12 LyαForests
^{1} Department of Physics and Astronomy, University of Utah, 115 S 1400 E, Salt Lake City, UT 84112, USA
^{2} APC, Université Paris DiderotParis 7, CNRS/IN2P3, CEA, Observatoire de Paris, 10 rue A. Domon & L. Duquet, 75013 Paris, France
^{3} LPNHE, CNRS/IN2P3, Université Pierre et Marie Curie Paris 6, Université Denis Diderot Paris 7, 4 place Jussieu, 75252 Paris CEDEX, France
^{4} IRFU, CEA, Université ParisSaclay, 91191 GifsurYvette, France
email: rich@hep.saclay.cea.fr
^{5} Aix Marseille Université, CNRS, LAM (Laboratoire d’Astrophysique de Marseille) UMR 7326, 13388 Marseille, France
^{6} Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA
^{7} Department of Physics and Astronomy, University College London, Gower Street, London, UK
^{8} Laboratoire d’Astrophysique, École Polytechnique Fedérale de Lausanne, 1015 Lausanne, Switzerland
^{9} Department of Physics and Astronomy, University of California, Irvine, CA 92697, USA
^{10} Brookhaven National Laboratory, 2 Center Road, Upton, NY 11973, USA
^{11} HarvardSmithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138, USA
^{12} Institució Catalana de Recerca i Estudis Avançats, 08010 Barcelona, Catalonia, Spain
^{13} Institució de Ciéncies del Cosmos, Universitat de Barcelona (UBIEEC), 08007 Barcelona, Catalonia, Spain
^{14} Université Paris 6 et CNRS, Institut d’Astrophysique de Paris, 98bis boulevard Arago, 75014 Paris, France
^{15} Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh, EH9 3HJ, UK
^{16} Department of Astronomy and Astrophysics, The Pennsylvania State University, University Park, PA 16802, USA
^{17} Institute for Gravitation and the Cosmos, The Pennsylvania State University, University Park, PA 16802, USA
^{18} Department of Astronomy, Ohio State University, 140 West 18th Avenue, Columbus, OH 43210, USA
Received: 31 January 2017
Accepted: 20 March 2017
We have used fluxtransmission correlations in Lyα forests to measure the imprint of baryon acoustic oscillations (BAO). The study uses spectra of 157 783 quasars in the redshift range 2.1 ≤ z ≤ 3.5 from the Sloan Digital Sky Survey (SDSS) data release 12 (DR12). Besides the statistical improvements on our previous studies using SDSS DR9 and DR11, we have implemented numerous improvements in the analysis procedure, allowing us to construct a physical model of the correlation function and to investigate potential systematic errors in the determination of the BAO peak position. The Hubble distance, D_{H} = c/H(z), relative to the sound horizon is D_{H}(z = 2.33) /r_{d} = 9.07 ± 0.31. The bestdetermined combination of comoving angulardiameter distance, D_{M}, and the Hubble distance is found to be D_{H}^{0.7}D_{M}^{0.3} /r_{d} = 13.94 ± 0.35. This value is 1.028 ± 0.026 times the prediction of the flatΛCDM model consistent with the cosmic microwave background (CMB) anisotropy spectrum. The errors include marginalization over the effects of unidentified highdensity absorption systems and fluctuations in ultraviolet ionizing radiation. Independently of the CMB measurements, the combination of our results and other BAO observations determine the openΛCDM density parameters to be Ω_{M} = 0.296 ± 0.029, Ω_{Λ} = 0.699 ± 0.100 and Ω_{k} = −0.002 ± 0.119.
Key words: cosmological parameters / dark energy
© ESO, 2017
1. Introduction
The sound waves that propagated in the prerecombination Universe produced a pronounced peak in the twopoint correlation function of the cosmologicaldensity field (Peebles & Yu 1970; Sunyaev & Zeldovich 1970). This “baryonacoustic oscillation” (BAO) peak, first observed by Eisenstein et al. (2005) and Cole et al. (2005), is centered on a comoving distance equal to the sound horizon at the drag epoch, r_{d}. Observed at a redshift z, the BAO peak position in the transverse (angular) direction determines the ratio D_{M}(z) /r_{d}, where D_{M}(z) = (1 + z)D_{A}(z) is the “comoving” angulardiameter distance (D_{A} is the “traditional” angulardiameter distance). In the radial (redshift) direction, the peak position determines D_{H}(z) /r_{d}, where D_{H}(z) = c/H(z) is the Hubble distance. Since D_{M}, D_{H} and r_{d} all have simple dependencies on the cosmological parameters, observation of the BAO feature constrains those parameters, especially when combined with cosmic microwave background (CMB) data (Planck Collaboration XVI 2014; Planck Collaboration XIII 2016; Aubourg et al. 2015). In particular, one can constrain models beyond the flatΛCDM model that describes CMB data, deriving constraints on cosmological curvature and the darkenergy equation of state. Furthermore, the shape of the spectrum of CMB anisotropies can be used to calculate the value of r_{d} to percentlevel precision, r_{d} ~ 147 Mpc and the use of this value allows one to derive D_{M}(z) and D_{H}(z) from BAO measurements. These absolute distances can be combined with relative distances determined with type Ia supernovae (Betoule et al. 2014), to extrapolate to z = 0, yielding a “topdown” measurement of H_{0} (Aubourg et al. 2015).
Most studies of the BAO peak have used galaxies at redshifts z< 0.8 as tracers of the density. The first observations used z ~ 0.35 data from the Sloan Digital Sky Survey (SDSS; Eisenstein et al. 2005) and z ~ 0.2 data from the TwoDegree Field Redshift Survey (Cole et al. 2005), and the combination of these sets (Percival et al. 2007, 2010). Since then, results of increasingly higher precision have been obtained, most significantly in the redshift range 0.35 <z< 0.65 from the Baryon Oscillation Spectroscopy Survey (BOSS) of SDSSIII (Anderson et al. 2012, 2014b,a) with the results of the complete survey being summarized in Alam et al. (2016). Results at other redshifts have been obtained by 6dFGRS at z ~ 0.11 (Beutler et al. 2011), WiggleZ at 0.4 <z< 0.8 (Blake et al. 2011), SDSSI at z ~ 0.35 (Padmanabhan et al. 2012; Mehta et al. 2012; Chuang & Wang 2012; Xu et al. 2013) and SDSSI at z ~ 0.15 (Ross et al. 2015). There is an impressive agreement of the results of these studies with the expectations of flatΛCDM models based on CMB data, as emphasized by Planck Collaboration XIII (2016).
At higher redshifts, BAO correlations can be seen using quasars and their Lymanα (Lyα) forests as mass tracers (McDonald & Eisenstein 2007). The correlations in the Lyαforest fluxtransmission field of BOSS quasars were first studied in Slosar et al. (2011) and the BAO peak was seen in SDSS DR9 (Busca et al. 2013; Slosar et al. 2013; Kirkby et al. 2013) and DR11 (Delubac et al. 2015). Crosscorrelations of the Lyα absorption with the distribution of quasars were detected in DR9 (FontRibera et al. 2013), and the first BAO detection was presented in FontRibera et al. (2014).
In this paper, we study the autocorrelation function, ξ, of the Lyα fluxtransmission field using the SDSS DR12 and we update the cosmological constraints reported in Delubac et al. (2015). Our study of the quasarforest crosscorrelation function will be presented in a future publication (du Mas des Bourboux et al., in prep.). In addition to using the 15% increase of survey area of DR12 over DR11, the following improvements over the analysis of Delubac et al. (2015) have been implemented:

Spectroscopic pipeline improvements (Sect. 2) that includes anew algorithm for the extraction of spectra from the CCD imagesthat results in a more linear flux response. This modificationallows us to correct for the mean distortion of the fluxtransmission field due to imperfect spectral modeling of standardstars. We also correct for the differential positioning of quasar andstellar fibers in the focal plane due to optimization at differentwavelengths (Margala et al. 2016).

The use of mock spectra with improved modeling of metal absorbers (Sect. 3), including both Lyαmetal and metalmetal correlations.

Modeling of the distortions of the correlation function due to quasarcontinuum fitting (Sect. 4). This allows us to construct a physical model of the correlation function over the range 10 <r< 200 h^{1} Mpc.

Modeling of spurious correlations introduced by the pipeline (Sect. 5), calculation of their effect on the correlation function, and searches for unidentified spurious correlations using the CarbonIV (CIV) forest, 142.0 <λ_{RF}< 152.0 nm.

Fits of the data (Sect. 6) that marginalize over the contributions to the correlation function of metals, unidentified highcolumndensity systems, and fluctuations of ionizing UV flux.
This paper is organized as follows. Section 2 describes the DR12 data used in this analysis. Section 3 gives a brief description of the mock spectra used to test the analysis procedure, with a more detailed description being found in Bautista et al. (2015). Section 4 presents our method of estimating the fluxtransmission field, its correlation function, and the associated covariance matrix. Section 5 studies spurious correlations induced by the pipeline. In Sect. 6 we fit the mocks and data to derive the BAO peak position parameters, D_{M}(z = 2.33) /r_{d} and D_{H}(z = 2.33) /r_{d}. Section 7 investigates possible systematic errors in the measurement. In Sect. 8 we compare our measured peak position with the predictions of ΛCDM models and derive constraints on the associated cosmological parameters. Section 9 presents a brief summary.
Fig. 1
SDSS DR12 footprint (in J2000 equatorial coordinates) used in this work. The survey covers one quarter of the sky (10^{4}deg^{2}). The light blue regions are those added beyond the area covered by our previous study (Delubac et al. 2015). The dotted line is the Galactic plane. 
2. The BOSS quasar sample and data reduction
The quasar sample was gathered over a fiveyear period by the SDSSIII Collaboration (Eisenstein et al. 2011; Gunn et al. 1998, 2006; ?). We use the data from the twelfth data release (DR12) of SDSS as presented in Alam et al. (2015). The associated quasar catalog is described in Pâris et al. (2017). Most of the quasar spectra were obtained by the Baryon Oscillation Spectroscopic Survey, BOSS (Dawson et al. 2013), but DR12 also includes six months of data from the SEQUELS^{1} program. The DR12 celestial footprint covering ~π sr ~ 10^{4} deg^{2} is displayed in Fig. 1. An example of a quasar spectrum in the Lyαforest region is shown in Fig. 2. The redshift distribution of measurement pairs in the forest is shown in Fig. 3.
The quasar target selection used in BOSS, summarized in Ross et al. (2012), combines different targeting methods described in Yèche et al. (2010), Kirkpatrick et al. (2011), and Bovy et al. (2011). The selection algorithms use SDSS photometry and, when available, data from the GALEX survey (Martin et al. 2005) in the UV; the UKIDSS survey (Lawrence et al. 2007) in the NIR, and the FIRST survey (Becker et al. 1995) in the radio.
The DR12 data were processed using a new software package that differs from the standard DR12 SDSSIII pipeline (Bolton et al. 2012) and which has become the standard pipeline for SDSS DR13 (Albareti et al. 2016). For each object, both pipelines provide a flux calibrated spectrum, f(λ), errors, and an object classification (galaxy, quasar, star). A model spectrum is fit to f(λ) providing a redshift estimate. For this study, we use the “coadded” spectra constructed from typically four exposures of 15 min resampled at wavelength pixels of width Δlog _{10}λ = 10^{4} (cΔλ/λ ~ 69 km s^{1}). For the small number of quasars with repeated observations, the coadded spectra can include exposures widely separated in time.
An important difference with respect to the DR12 pipeline is that pixels on the CCD image are combined to give a flux f(λ) with pixelweights determined only by the CCD readout noise. While this method is suboptimal because it ignores photoelectron Poisson noise, compared to the DR12 method it yields an unbiased estimate of f(λ) since the weights do not depend on the observed CCD counts which are needed to estimate Poisson noise. A more detailed description of the changes to the extraction pipeline is given in Appendix A.
The ratio of observed flux to model flux, averaged over all spectra, is an important diagnostic of pipeline systematic errors. Figure 4 shows the average ratio, R(λ), as a function of observed wavelength for quasar spectra on the red side of Lyα emission. Since model imperfections for individual quasar spectra are averaged over in this figure (because of the range of quasar redshifts), the deviations from unity reflect imperfections in the spectrograph flux calibration. The figure reveals percentlevel deviations that are mostly due to imperfect modeling of photospectroscopic standard stars. The calcium H and K lines from interstellar absorption are also visible. In the analysis to be presented here, the f(λ) are divided by R(λ) to correct on average for these artifacts. This procedure is effective only with the fluxes from the DR13 pipeline since the nonlinearities present in the DR12 pipeline make the correction fluxdependent. We show in Sect. 5 that after this global correction, remaining calibration artifacts (due to their timedependence) are sufficiently small to have a negligible effect on the measurement of the correlation function.
The spectra of all quasar targets were visually inspected (Pâris et al. 2012, 2014, 2017) to correct for misidentifications, to flag broad absorption lines (BALs), and to establish the definitive quasar redshift. Damped Lyα troughs (DLAs) were visually flagged, but also identified and characterized automatically (Noterdaeme et al. 2012). The visual inspection of DR12 confirmed 297 301 quasars, of which 181 719 are in the redshift range appropriate for this study, 2.1 ≤ z_{q} ≤ 3.5. We discarded quasars with visually identified BALs (to avoid the necessity of modeling their profiles in the forest) leaving 160 868 quasars. A further cut requiring a minimum number of unmasked forest pixels (50 “analysis pixels”; see below) yielded a sample of 157 922 quasars. Finally, 139 spectra failed the continuumfitting procedure (Sect. 4.1), leaving 157 783 spectra compared to 137 562 in the Delubac et al. (2015) investigation.
Fig. 2
Example of a BOSS quasar spectrum of redshift 2.91 (smoothed to the width of analysis pixels). The red and blue lines cover the forest region used in our analysis, 104.0 <λ_{RF}< 120.0 nm. This region is sandwiched between the quasar’s Lyβ and Lyα emission lines, respectively at 400.9 and 475.4 nm (restframe 102.572 and 121.567 nm). The blue line is the model of the continuum, C_{q}(λ); the red line is the product of the continuum and the mean transmission, , as calculated by the method described in Sect. 4. 
Fig. 3
Weighted redshift distribution of pairs of Lyα forest pixels. The mean is ⟨ z ⟩ = 2.33. Included in the distribution are the ~5 × 10^{10} pairs within 20 h^{1} Mpc of the center of the BAO peak. 
For the measurement of the flux transmission, we use the restframe wavelength interval (1)As illustrated in Fig. 2, this range is bracketed by the emission lines λ_{Lyβ} = 102.572 nm and λ_{Lyα} = 121.567 nm. This region was chosen as the maximum range that avoids the large pixel variances on the slopes of the two lines due to quasartoquasar diversity of lineemission strength and profile. The absorber redshift, z = λ/λ_{Lyα}−1, is required to lie in the range 1.96 <z< 3.44. The lower limit is set by the requirement that the observed wavelength be greater than 360 nm, below which the system throughput is less than 10% of its peak value. The upper limit is produced by the maximum quasar redshift of 3.5, beyond which the BOSS surface density of quasars is not sufficient to be useful for this study. The weighted distribution of redshifts of absorber pairs near the BAO peak position is shown in Fig. 3. The distribution has a mean of .
For the determination of the correlation function, we use analysis pixels that are the inversevarianceweighted flux average over three adjacent pipeline pixels. Throughout the rest of this paper, “pixel” refers to analysis pixels unless otherwise stated. The width of these pixels is 207 km s^{1} corresponding at z ~ 2.33 to an observedwavelength width ~0.28 nm and a comoving radial distance of ~2.0 h^{1} Mpc. The total sample of 157 783 spectra thus provides ~3 × 10^{7} measurements of Lyα absorption over an effective volume of ~50 Gpc^{3}.
Fig. 4
Mean ratio, R(λ), of observed flux to pipelinemodel flux as a function of observed wavelength for quasar spectra to the red of the Lyα emission line (λ_{RF}> 130 nm). (The mean is calculated by weighting each measurement by the inverse of the pipeline variance.) In this mostly unabsorbed region of quasar spectra, the percentlevel wavelengthdependent deviations from unity are due to imperfect modeling of calibration stars and to the calcium H and K lines (393.4 and 396.9 nm) due to Galactic absorption. 
3. Mock quasar spectra
In order to test the analysis procedure and investigate statistical and possible systematic errors, we created 100 sets of mock spectra that reproduce the essential physical and instrumental characteristics of the BOSS spectra. The basic method for the production of the mocks with Lyα absorption is described in FontRibera et al. (2012). Except for the implementation of absorption due to metals and high columndensity systems (HCDs, i.e., damped Lyαsystems and Lymanlimit systems), the mocks used in this study are identical to the DR11 mocks (Bautista et al. 2015) that were used in Delubac et al. (2015). The DR11 mocks therefore cover ~15% less solid angle than the DR12 data.
For each set of spectra, the background quasars were assigned the angular positions and redshifts of the DR11 quasars. The unabsorbed spectra (continua) of the quasars were generated using the Principal Component Analysis eigenspectra of Suzuki et al. (2005). The amplitudes for each eigenspectrum were randomly drawn from Gaussian distributions with a dispersion equal to that of the corresponding eigenvalues in Table 1 of Suzuki (2006). The overall normalization was chosen by fitting the mock spectrum to the corresponding observed spectrum.
The Lyα absorption field was generated using the technique described in FontRibera et al. (2012) in which a Gaussian random field, δ_{G}, is first defined at the positions of the ~10^{7} observed forest pixels. Flux transmissions were defined by the nonlinear transformation The two functions, a(z) and b(z) are chosen so that the resulting mean transmission and variance are near those observed in the data. The correlations of the field δ_{G} are chosen so that the correlations of F(δ_{G}) follow the linear correlations of the “mock” cosmology of Table 1 modified by nonlinear effects (McDonald 2003). For pixels randomly distributed in space, this procedure would involve inverting a 10^{7} × 10^{7} matrix. To reduce the problem to a manageable size, use was made of the fact that the forest pixels are nearly parallel, allowing a separate treatment of radial and transverse coordinates.
In the final step, the spectra were modified to include the effects of the BOSS spectrograph point spread function (PSF), readout noise, photon noise, and flux systematic errors.
For each of the 100 mock data sets, four types of spectra were produced and analyzed. The first type consists simply of Lyα fluxtransmission fraction, F_{Lyα}(λ), modified for the wavelength resolution but without multiplication by a quasar continuum spectrum, C(λ). Analysis of this mock type allowed us to study the recovery of the BAO peak position under the most favorable conditions. With the introduction of the quasar continuum, the second type consists of more realistic spectra, F_{Lyα}(λ)C(λ). Analysis of this type tests our ability to fit the quasar continuum and to model the resulting distortion of the correlation function.
The final two types of spectra add to F_{Lyα}(λ) absorption that is due to HCDs and to metals. Following FontRibera & MiraldaEscudé (2012), HCDs are placed at randomly chosen pixels where the optical depth was above a chosen threshold. The neutralhydrogen column densities were drawn randomly with N_{HI}> 10^{17.2} cm^{2}, assuming an intrinsic powerlaw distribution corrected for selfshielding and normalized to match the observations of Prochaska et al. (2005). Because the HCDs are placed in redshift space, the resulting HCDs have correlations that trace the underlying density field but with a redshiftspace distortion parameter that is not necessarily equal to that of physical HCDs.
Fig. 5
Onedimensional fluxcorrelation function, ξ_{1d}, for BOSS quasars showing correlations of δ_{q}(λ) within the same forest. The correlation function is shown as a function of wavelength ratio for the data and for the mocks (procedure Met1). Prominent peaks due to Lyαmetal and metalmetal correlations are indicated. The peak at λ_{1}/λ_{2} ~ 1.051, which is seen in the data but not in the mocks, is due to CII(133.5)SiIV(140.277) at z ~ 1.85, outside the redshift range covered by the mocks. 
Absorption due to metal transitions was also simulated by adding absorption in proportion to the Lyα absorption at the same redshift. The important transitions can be seen in the “one dimensional” correlation function, ξ_{1d}(λ_{1}/λ_{2}) for fluxes at two wavelengths within the same forest, as shown in Fig. 5. In this figure, the peaks correspond to absorption by two different transitions by material at the same physical position but a different wavelength. The most prominent features are due to pairs comprising Lyα and one metal transition, but there are also features due to metalmetal pairs. More detailed information can be obtained by stacking spectra around strong absorbers in the forest (Pieri et al. 2010, 2014). Table 2 lists the transitions included in the mock spectra.
We have tested two procedures for adding metal absorption to the primary Lyα absorption. The methods are defined by the assumed relation between Lyα absorption and metal absorption at the same physical position. For each wavelength, λ, with Lyα absorption A_{Lyα}(λ) = 1−F_{Lyα}(λ), an appropriate metal absorption A_{m}(λ + Δλ_{m}) for each transition, m, is chosen. A simple procedure would be to choose A_{m} ∝ A_{Lyα} with a proportionality constant chosen to reproduce the Lyαmetal features in the observed ξ_{1d}. However, this procedure does not produce significant features corresponding to metalmetal correlations. To do this it is necessary to add, for a randomly chosen small fraction of wavelengths, a much larger metal absorption.
Parameters of the flat ΛCDM cosmological model used for the production and analysis of the mock spectra and the CMBbased flat ΛCDM model from Planck Collaboration XIII (2016) used for the analysis of the data.
The first procedure (hereafter Met1) aims to match the observed lineofsight correlation function ξ_{1d}, in particular the amplitudes of the peaks associated to Lyαmetal and metalmetal correlations. We first transform the pure mock Lyα flux into “optical depth” τ_{Lyα} and define the metal absorption of transition m as τ_{m} = a_{m}τ_{Lyα}. A quadratic term in τ_{Lyα} is added for a fraction of the strong Lyα absorbers to simulate strong metal absorption. The parameters a_{m}, the quadratic terms and rates of strong absorption are set to match the observed amplitudes in the data ξ_{1d}.
The second procedure (hereafter Met2) aims to match the observed stack of high signaltonoise ratio Lyα absorbers, while maintaining consistency with the ξ_{1d} in the resulting mock ξ_{1d}. Following Pieri et al. (2014), stacks were produced using the DR12 Lyα forest sample. As a function of the Lyα transmission, we measured amplitudes for the metal absorption features. The obtained “flux decrements” caused by metals were treated as target average flux decrements. Metals are implemented in Met2 mocks as a mix of weak and strong metal absorption, generating a mock ξ_{1d} consistent with the observed ξ_{1d} and consistent with the absorption frame measurements of Pieri et al. (2014).
The metal absorption added by these methods is due to the presence of metals at the same physical position as HI. Because our mocks provide the HI density only for redshifts z> 1.9, we cannot generate absorption at lower redshifts that nevertheless appears in the Lyα forest because of transition wavelengths much greater than that of Lyα. An important example is CIV doublet ( nm) where the absorption at an observed wavelength of 400 nm is due to material at z ~ 1.6. Fortunately this absorption has little effect on the correlation function, as we see in Sect. 6.
The statistical properties of metal absorption in the mocks are determined by the underlying density field. However, the analysis procedure interprets absorption at a given wavelength as absorption due to the Lyα transition. Because of this behavior, the metal contribution to the measured correlation function is shifted and deformed in (r_{⊥},r_{∥}) space. In particular, the large correlation due to HI and metals at the same physical position is seen at (r_{⊥} ~ 0,r_{∥}) with r_{∥} = (1 + z)D_{H}(z)Δλ/λ where Δλ/λ is the relative wavelength separation of the metal feature with respect to Lyα. This leads to a large correlation at this r_{∥} (and r_{⊥} ~ 0), so if the amplitude is significant, it can be measured. Table 2 lists the apparent r_{∥} corresponding to vanishing Lyαmetal separation.
Transitions contributing to absorption by metals in the mock spectra.
4. Measurement of fluxtransmission field and its correlation function
In this section we describe the measurement of the correlation function of the transmitted flux fraction: (2)Here, f_{q}(λ) is the observed flux density for quasar q at observed wavelength λ, C_{q}(λ) is the unabsorbed flux density (the socalled “continuum”) and is the mean transmitted fraction at the absorber redshift, z(λ) = λ/λ_{Lyα}−1. Measurement of the fluxtransmission field δ_{q}(λ) requires estimates of the product for each quasar. An example is shown in Fig. 2. The estimation procedures, described in this section, differ slightly from those of our previous studies (Busca et al. 2013; Delubac et al. 2015). One important modification is that we now calculate the distortion of the correlation function due to continuumfitting (Sect. 4.3).
4.1. Estimation of δ_{q}(λ)
As in our previous investigations, we assume the quasar continuum, C_{q}(λ), is the product of a universal function of the restframe wavelength, λ_{RF} = λ/ (1 + z_{q}), and a linear function of log λ, included to account for quasar spectral diversity: (3)with C(λ_{RF}) being normalized so that its integral over the forest is equal to unity. The (a_{q},b_{q}) and C(λ_{RF}) are determined by maximizing the likelihood function given by (4)Here P(f_{q}(λ)  C_{q}(λ)) is the probability to observe a flux f_{q}(λ) for a given continuum found by convolving the intrinsic probability, D(F = f_{q}(λ) /C_{q}(λ),z), with the observational resolution assumed to be Gaussian: (5)Here, σ_{q}(λ)^{2} is the variance due to readout noise and photon statistics.
Fig. 6
Mean transmission, , calculated using the DR12 pipeline of our previous investigations, and using the DR13 pipeline applied to DR12 data of this analysis. The DR13 pipeline allows us to divide forest spectra by the correction factor (Fig. 4) derived from the spectra on the red side of Lyα emission. This procedure eliminates most of the small scale structure found with the old pipeline that is due to artifacts of stellar modeling, Galactic calcium H and K absorption, and sky lines. 
A continuumdetermination method is defined by the assumed form of D(F,z). For this work, it is taken to be the lognormal model of absorption used to generate the mock data, corresponding to “method 2” of our previous studies (Delubac et al. 2015). As a check, we also use “method 1”, which is equivalent to using a narrow Gaussian for D(F,z), thereby producing only the product for each forest.
In practice, we maximize the likelihood iteratively by assuming a C(λ_{RF}) to determine the (a_{q},b_{q}). The mean absorption is then calculated by an appropriately weighted average of f_{q}(λ) /C_{q}(λ) (for fixed λ) after which C(λ_{RF}) is recalculated as an average of (for fixed λ_{RF}). The procedure is stopped after ten iterations, at which point a stability of ~10^{4} is reached for C(λ_{RF}) and . Figure 2 shows a spectrum with its C_{q}(λ) and .
The function to be used in Eq. (2) to calculate δ_{q}(λ) is the mean of f_{q}(λ) /C_{q}(λ) and can thus differ slightly from used in the model to estimate the C_{q}(λ). The use of the mean of f_{q}(λ) /C_{q}(λ) ensures that the mean of δ is zero at each redshift (or wavelength).
Figure 6 displays the calculated , both for the DR13 pipeline used here and the DR12 pipeline used previously. The use of the new pipeline removes most of the artifacts in the old analysis. We emphasize, however, that the derived is dependent on the assumed form of D(F,z) and should therefore not be considered as a measurement of the mean absorption. For example, the flattening of for z< 2.1 (λ< 377 nm) suggests that we have slightly underestimated at these redshifts. Since, by construction, the mean δ vanishes at each redshift, this implies that our procedure makes a compensating overestimate of the C_{q}(λ). Since it is the product that determines δ_{q}(λ), the measured correlation function is therefore not strongly affected.
Those forests with identified DLAs are given a special treatment. All pixels where the absorption due to the DLA is higher than 20% are not used. The absorption in the wings is corrected using a Voigt profile following the procedure of Noterdaeme et al. (2012).
We denote as the estimate of δ_{q}(λ) using the relation (2). Because forest data is used to fit for C_{q}(λ) and the measured is not equal to the original δ_{q}(λ). We can identify two effects. First, the use of is equivalent to the transformation (6)where the overbar refers to the average over forests at fixed λ. Second, the fitting of (a_{q},b_{q}) with the data biases toward zero the mean and mean . To simplify this effect and facilitate its correction we explicitly subtract the mean and first moment for each forest (7)where the overbars refer to averages within a given forest using the weights defined in the next section. The accompanying distortion of the correlation function is nonnegligible, as we demonstrate in Sect. 6.2 with the mock spectra.
The transformation (7) has two interesting effects on the measured fluxtransmission field. Most importantly, it makes it simple to calculate the distortion of the correlation function (Sect. 4.3) and thus simplifies the relation between the underlying physical model and the measured correlation function. Second, it nearly eliminates the difference between the correlations functions calculated with the two continuum fitting methods used in Delubac et al. (2015) with the rms difference in the two ξ(r_{⊥},r_{∥}) being 0.056 of the rms uncertainty per (r_{⊥},r_{∥}) bin.
4.2. Estimation of the correlation functions
For the estimator of the flux autocorrelation function, we adopt a simple weighted sum of products of the : (8)where the w_{i} are weights (see below) and each i (or j) indexes a measurement on a quasar q_{i} at wavelength λ_{i}. The sum over (i,j) is understood to run over all pairs of pixels within a bin A in the space of pixel separations, (A) → (r_{⊥},r_{∥}). We exclude pairs of pixels from the same quasar to avoid the correlated errors in and arising from the estimate of C_{q}(λ) for the spectrum of the quasar. The bins A are defined by a range of width 4 h^{1} Mpc of the components perpendicular and parallel to the line of sight, r_{⊥} and r_{∥}. We use 50 bins in each component, spanning the range from 0 to 200 h^{1} Mpc; the total number of bins used for evaluating the correlation function is therefore 2500. Separations in observational pixel coordinates (RA, Dec, λ) are transformed to (r_{⊥},r_{∥}) in units of h^{1} Mpc by assuming that absorption is due to the Lyα transition and using the cosmological parameters from Table 1 (Planck cosmology for the data and the mock cosmology for the mocks).
As described in Delubac et al. (2015), the weights, w_{i}, are chosen so as to account for both Poisson noise in the flux measurement and for the intrinsic fluctuations in due to cosmological largescale structure. The weights are set to zero for pixels flagged by the pipeline as having problems due, for example to sky emission lines or cosmic rays. To reduce the pipeline systematics discussed in Sect. 5, we also do not use pairs of pixels that have nearly the same wavelength (r_{∥}< 4 h^{1} Mpc) and that were taken on the same exposures.
4.3. The distortion matrix
The transformations (6) and (7) mix pixels so that the correlation is equal to the original ⟨ δ_{q}(λ)δ_{q′}(λ′) ⟩ plus a linear combination of the correlations of other pixelpairs ⟨ δ_{q′′}(λ′′)δ_{q′′′}(λ^{′′′}) ⟩. This approach means that the measured correlation function is a “distorted” version of the true correlation function ξ. Since the transformations (6) and (7) are linear, the relation between measured and true correlation functions is given by a distortion matrix D_{AB}: (9)where A and B refer to bins in pixel separation space. Writing produces (10)where w_{A} = ∑ _{i,j ∈ A}w_{i}w_{j}. We ignore the small effect of transformation (6), in which case η_{ii′} = 0 unless the pixels i and i′ are in the same forest: (11)where the two sums over k include only pixels in the same forest as that of i and i′ and the overbars refer to averages in that forest. The matrix D_{AB} thus depends only on the geometry and weights of the survey. We see its effect on the mock correlation function in Sect. 6.2 (Fig. 11).
Previous analyses (Busca et al. 2013; Slosar et al. 2013; Delubac et al. 2015; Blomqvist et al. 2015) have dealt with the distortions introduced by continuumfitting in different ways. Busca et al. (2013) and Delubac et al. (2015) model it as an additive power law in r and μ with 12 free parameters. Blomqvist et al. (2015) assumed that continuumfitting reduces the observed amplitude of longwavelength modes in the direction parallel to the line of sight. They then model it as a multiplicative function of k_{∥} that tends to zero at large scales (with k_{∥} ≲ 2π/L, where L ~ 380 h^{1} Mpc is the typical length of a forest) and to one at small scales (k_{∥} ≫ 2π/L). They tune the shape of the function using simulated data to ultimately reduce the number of free parameters to one.
In our approach we do not introduce free parameters to account for the effects of continuumfitting. Instead, we follow the assumption, first proposed by Slosar et al. (2013), that at each line of sight the continuumfit delta field differs from the true delta field by a linear function in Λ. Slosar et al. (2013) then deweight these “linear modes” in their covariance matrix. Alternatively, we use the transformation from δ to to remove the “linear modes” from the model via the distortion matrix.
Fig. 7
Correlation, Corr, plotted vs. for the two lowest intervals of . The correlation is averaged over . The left panel is the Metal (Met1) mocks and the right panel is for the data. As labeled, the correlation is calculated by subsampling (Eq. (13)), by assuming “independent forests” (Eq. (14), offset by 0.5 h^{1} Mpc for clarity), or by the mocktomock variation (Eq. (15)). Good agreement between the methods is seen, though the independentforest method necessarily underestimates the correlation for Δr_{⊥} ≠ 0. The differences between the mocks and the data reflect the differences in ξ_{1d} (Fig. 5). 
4.4. The covariance matrix
The covariance matrix associated with is: (12)where W_{A} = ∑ _{ij ∈ A}w_{i}w_{j}. Following Delubac et al. (2015), we evaluate C_{AB} by dividing the BOSS footprint into subsamples and measuring and in each subsample s. Neglecting the small correlations between subsamples, and replacing the fourpoint function by the measured product of correlations in subsamples, the covariance is given by (13)where the are the sums of weights in the subsample s. As in Delubac et al. (2015), the SDSS plates define the subsamples.
As a check of the subsampling method, the sum (12) can also be estimated by neglecting interforest correlations, in which case the fourpoint function vanishes unless the four pixels are drawn from just two spectra: (14)The sum can then be estimated from a random sample of forest pairs. Because neighboring forests are nearly parallel, the sum necessarily gives C_{AB} = 0 unless .
Finally, for the analysis of the mock data, the covariance matrix can also by calculated from the mocktomock variations of ξ: (15)where the overbar refers to averages over the set of 100 mocks. Given the approximations used in the calculation of C_{AB} via Eqs. (13) or (14), it is of great importance that the mocktomock variations confirm the accuracy of the other two methods when applied to mock data. This comparison is shown in Fig. 7.
The 2500 × 2500 element matrix, C_{AB}, has a relatively simple structure. By far the most important elements are the diagonal elements which are, to good approximation, inversely proportional to the number of pixel pairs used in the calculation of the correlation function; the number of pairs is roughly proportional to r_{⊥}:(16)The variance is about twice as large as what one would calculate assuming all (analysis) pixels used to calculate ξ(r_{⊥},r_{∥}) are independent. This decrease in the effective number of pixels is due to the physical and instrumental correlations between neighboring pixels in a given forest (Eq. (14)).
The offdiagonal elements of the covariance matrix also have a simple structure. As previously noted, the covariance is mostly due to the twoforest part of the fourpoint function which, because neighboring forests are nearly parallel, only contribute to the covariance matrix elements with . This behavior is illustrated in Fig. 7, which displays the mean values of the correlation matrix elements as a function of for the smallest values of .
The statistical precision of the subsampling calculation is ~0.02 for individual elements of the correlation matrix. Figure 7 reveals that only correlations with Δr_{⊥} ~ 0 and Δr_{∥}< 20 h^{1} Mpc are greater than the statistical precision and therefore sufficiently large for individual matrix elements to be measured accurately by subsampling. As in Delubac et al. (2015), we therefore average the subsampling correlation matrix over and use the resulting covariance matrix that is a function of .
5. Correlations introduced by the optics and data pipeline
Spurious correlations in δ_{q}(λ) are introduced by the telescope and spectrometer optics, and by the pipeline reductions of the data. These correlations are superimposed on the physical ξ(r_{⊥},r_{∥}) that we wish to measure; in this section we estimate the various contributions.
5.1. Optical crosstalk
At the optical level, correlations are introduced by signals from one object “scattering” into the spectra of other objects, that is through optical crosstalk. There is negligible crosstalk introduced in the telescope focal plane by photons from one object entering into the fiber of another object. However, there can be measurable crosstalk between neighboring fibers downstream of the wavelength dispersion where they are focused on the CCD. This contamination arises through imperfect modeling of the pointspread function when transforming the twodimensional CCD image into a series of onedimensional spectra, f(λ). We have measured the crosstalk as a function of fiber separation by fitting the signal in sky fibers to a sky model and a weighted sum of the spectra of the objects (quasars or galaxies) in neighboring fibers. Consistent with the results of Croft et al. (2016), we find the crosstalk for the neighboring fibers is ~0.2% and ~0.05% for fibers separated by two or three rows respectively, and is consistent with zero for larger separations.
The crosstalk directly introduces correlations between pixels at the same wavelength but these pixel pairs are, at any rate, not used in the analysis. The crosstalk introduced between pixels λ_{1} and λ_{2} of two quasars is proportional to the product of the crosstalk amplitude and to ξ_{1d}(λ_{1}/λ_{2}). We have verified that these correlations are insignificant compared to the measured ξ. This fact is in part due to the fiberassignment strategy which avoided placing two quasar candidates on neighboring fibers.
5.2. Pipelineinduced correlations
The pipeline treatment of the spectrometer data transforms flatfield corrected CCD counts, , to fluxes, f_{q}(λ). This process requires subtracting a sky contribution, B_{q}(λ), and multiplying by a calibration vector, A_{q}(λ), that corrects for the wavelengthdependent throughput of the system: (17)Both A_{q}(λ) and B_{q}(λ) are determined from spectra taken simultaneously with the data: spectra of spectrophotometric standard stars for A_{q}, and spectra of “sky fibers” (fibers pointing to empty sky regions) for B_{q}. The spectra from the 1000 fibers of a given exposure are treated independently for the two 500fiber spectrographs because of possible differential variations of instrumental properties such as throughput, PSF, and scattered light. (See Smee et al. 2013, for the relevant details of the spectrometer construction.) Hereafter, the collection of fibers of a given plate assigned to one of the two spectrographs will be referred to as a “halfplate.” Inaccuracies in the determinations of A_{q} and B_{q} will lead directly to correlated inaccuracies in the f_{q}(λ) of a given halfplate.
Errors in the sky subtraction, B_{q}(λ), are simple to model because they are mostly due to wellunderstood Poisson fluctuations of photoelectron counts in the ~50 sky fibers per halfplate. The pipeline interpolates these measurements in the focal plane to produce the sky background to be subtracted from a given quasar. The limited number of skyfibers results in correlations between f_{q}(λ) and f_{q′}(λ′) for q ≠ q′ and λ = λ′ since the sky subtractions where calculated using the same noisy skyfibers. These correlations for λ = λ′ generate small correlations for λ ≠ λ′ through continuum fitting, as discussed in Sect. 4.2. Because the correlations are primarily for λ = λ′ they contribute a spurious ξ(r_{⊥},r_{∥} = 0) for pairs of pixels on the same halfplate.
Errors on the calibration vector A_{q} are also Poissonian in that they are due to the fluctuations in the imperfect modeling of ~10 calibration stars per halfplate. The mean over all spectra of the imperfections are visible in the mean “transmission” of the unabsorbed parts of quasar spectra shown in Fig. 4. As mentioned in Sect. 2, the mean imperfections are removed by dividing forest spectra by R(λ) from Fig. 4. However, fluctuations from exposure to exposure remain because of the spectral variations of the small number of calibration stars used. The statistical properties of these variations have been estimated by comparing the A_{q} measured with random subsamples of stars and, independently, by comparing the A_{q} on different halfplates for the same exposure. Unlike the correlations due to skysubtraction, which are confined to λ = λ′, the errors in the calibration vectors are strongly correlated over the characteristic wavelength range of stellar spectral features, Δλ ~ 1 nm.
We constructed a model for pipelineinduced correlations that uses the measured statistical properties of the A_{q}(λ) and B_{q}(λ). The contribution of the sky noise was estimated using a fit of random realizations of the signal in the sky fibers (using a degree 2 polynomial fit as a function of the fiber number). The realizations were based on the statistical uncertainty estimation from the pipeline (which we have found to be accurate at the 5% level, see Fig. A.1, bottom panel). A normalization factor was applied to account for the average difference of calibration between the sky and target fibers (a ~10% correction). Performing a fit allowed us to capture accurately the resulting correlations as a function of fiber separation.
The contribution of the calibration was estimated using the fact that because the flux calibration of the two spectrographs is determined independently, the difference of those two calibration solutions for each exposure provides us with an estimate of their statistical fluctuation (up to a normalization factor ) while retaining the correlation of those fluctuations as a function of wavelength. Those differences for all DR12 plates was computed. Their average per observation run was subtracted in order to account for the actual difference of throughput of the two spectrographs and their evolution during the course of the survey. This data set provided us with a library of calibration uncertainties that was used to calculate calibrationinduced correlations between spectra.
This model was used to calculate the expected correlation between pixels in the Lyα forest of one quasar and pixels in the CIV forest (142 <λ_{RF}< 152 nm) of another quasar. The physical correlation of the two forests is due mostly to the autocorrelation of the weak CIV absorption in the two forests, but any pipelineinduced correlations should be present at full strength. Since this crosscorrelation is designed to isolate pipelineinduced effects, we denote this correlation function as ξ^{pl}(λ_{2}−λ_{1},θ) where (λ_{1},λ_{2}) are the wavelengths in the two forests and θ is the angular separation. We also distinguish between ξ^{pl} and correlations between flux pairs on the same halfplate, ξ^{pl}(hp = hp′), and correlations between flux pairs on different halfplates, ξ^{pl}(hp ≠ hp′). While these correlations are naturally given as a function of wavelength and angular separations, for convenience, the ξ^{pl} will be given as functions of pseudoseparations calculated using the Lyα restframe wavelength to define redshifts in both spectra. With this approach, corresponds to absorption at the same observed wavelength.
This correlation for samehalfplate pairs (measured using the techniques of Sect. 4.2) is shown with the red points in Fig. 8 for the first bin (). Superimposed on the data is the prediction of the model of the pipelineinduced correlations. For the first bin, the correlation is dominated by those induced by to the skysubtraction model. There is good agreement between this simple model and the observed correlations.
The pipelineinduced correlations that we have considered do not contribute to correlations between pixels observed on different halfplates. The differenthalfplates correlation functions for the bin is shown by the blue points in Fig. 8. The correlations for different halfplates (blue points) are clearly far less than those for same halfplates (red points). For , χ^{2} = 28.6 (N_{d.o.f.} = 15) for the nocorrelation hypothesis.
For the bins, the model predicts much smaller correlations. In particular, the skysubtraction model noise induces nonzero correlations only because of the continuum fit which distorts the original correlation function as described in Sect. 4.2. This behavior is illustrated in Fig. 9 showing ξ^{pl} in four ranges of μ. Because of the low level of absorption fluctuations redward of Lyα emission, the variance of the points is a factor ~7 smaller than the corresponding variance for the forestforest correlation function.
Fig. 8
Correlation between pixels in the Lyαforest and pixels in the CIVforest, ξ^{pl}(λ_{2}−λ_{1},θ) for small wavelength differences. The wavelength difference and angular separation have been transformed to pseudoseparations assuming Lyα absorption in both forests. The red points show the correlation for pairs on the same halfplate . The much smaller correlations for different halfplates are the blue points. The red line represents the prediction for our model of calibration and sky noise. 
The model for pipelineinduced correlations assumes that the skysubtraction errors are entirely due to Poisson statistics of photoelectron counts, neglecting possible wavelength and positiondependent systematic misestimates of the sky flux. While the changes applied to the pipeline (Appendix A) considerably reduced the systematic sky residuals, significant residuals remain on bright sky lines due to an imperfect PSF model, PSF variations and displacement of spectral traces (due to changes of temperature and variations of the gravity load of the spectrographs which are at the Cassegrain focus of the telescope). We tested the effect of those sky residuals by computing the correlation function of fake Lyα forests consisting of the residual signal in the nearest sky fiber in the CCD divided by the quasar continuum model (we made sure to use only once each sky fiber in this process to avoid introducing an artificial correlation). The measured correlation signal could be entirely explained by the Poissonian sky model noise described above. We hence conclude that the systematic sky residuals induce a negligible contamination to the Lyα correlation function.
This analysis presented above evaluates the contamination due to any additive signal to the Lyα forest that leads to a correlated signal in the CCD. It includes the systematic sky residuals, their fluctuations from plate to plate, but also the potential effect of scattered light in the spectrograph.
Fig. 9
Null test for pipelineinduced correlations. As in Fig. 8, the correlation between the Lyαforest pixels and the CIVforest pixels, is shown, but now for four angular ranges as labeled and offset by 0.1 successively for clarity. Pixel pairs with are excluded. There is no evidence for pipelineinduced systematics with the χ^{2} for vanishing correlation for is (32.5, 31.2, 22.7, 40.0) for the four ranges, each with 38 data points. 
We have used our model of pipelineinduced correlations to calculate their effect on the determination of the BAO peak position. Since the model contains no scale near the BAO scale, it is not surprising that it predicts no measurable influence, as reported in Table 7 of Sect. 7. Furthermore, to facilitate the comparison of the measured correlation function with the physical model developed in the next section, we do not use pairs of pixels on the same exposure that would contribute to the r_{⊥}< 4 h^{1} Mpc bins of ξ(r_{⊥},r_{∥}).
Fig. 10
Twodimensional representation of r^{2}ξ(r_{⊥},r_{∥}) in units of (h^{1} Mpc)^{2}. The right panel shows the measurement and the left panel the bestfit model, ξ^{ph}, modified by the distortion matrix via Eq. (18). The BAO feature is at r ~ 100 h^{1} Mpc. The effects of metalLyα correlations are seen in the lowest r_{⊥} bin, in particular the peak at 50 <r_{∥}< 70 h^{1} Mpc due to SiIIa and SiIIb. 
6. Fits for the BAO peak position
To determine the position of the BAO peak in the transverse and radial directions, we fit the measured ξ(r_{⊥},r_{∥}) to functions that describe the underlying largescalestructure correlations. These correlations are primarily due to Lyα absorption in the intergalactic medium (IGM), but we also include absorption by metals in the IGM and neutral hydrogen in highcolumndensity systems (HCDs). The physical correlations are corrected for distortions that are introduced by the procedure for determining the quasar continuum (Sect. 4.3). The fitting routine we use gives results consistent with those found using the publicly available baofit^{2} package (Kirkby et al. 2013; Blomqvist et al. 2015) modified to include the effects described in this section. The parameters of the fits are described below and in Table 3. The bestfit model correlation function is shown in Fig. 10.
Parameters of the physical model for the correlation function and their bestfit values for the data fit over the range 10 < r < 180 h^{1} Mpc.
6.1. The model of the correlation function
The general form for the model correlation function ξ_{A} of the (r_{⊥},r_{∥}) bin A is a distorted sum of the cosmological, or “physical” correlation function, , and a slowly varying function, B_{A}, used to test for systematics: (18)Here, D_{AA′} is the distortion matrix (Eq. (10)) that models the effects of continuum fitting.
The physical component of the model is dominated by the autocorrelation due to Lyα absorption. It is assumed to be a biased version of the total matter autocorrelation of the appropriate flatΛCDM model (Table 1) modified to free the position of the BAO peak (Kirkby et al. 2013): (19)where the BAO peakposition parameters to be fit are (20)and where the subscript “fid” refers to the fiducial cosmological model from Table 1 used to transform angle differences and redshift differences to (r_{⊥},r_{∥}). The nominal correlation function, ξ_{Lyα}(r_{⊥},r_{∥},α_{⊥} = α_{∥} = 1), is derived from its Fourier transform (21)where P_{QL} is the (quasi) linear power spectrum decomposed into a smooth component and a peak component corrected for nonlinear broadening of the BAO peak: (22)The smooth component is derived from the CAMBcalculated linear power spectrum, P_{L}, via the sideband technique (Kirkby et al. 2013) and P_{peak} = P_{L}−P_{smooth}. The correction for nonlinear broadening of the BAO peak is parameterized by . The nominal values used are Σ_{∥} = 6.41 h^{1} Mpc and Σ_{⊥} = 3.26 h^{1} Mpc (Kirkby et al. 2013).
In Eq. (21), the bias, b_{Lyα}, is assumed to have a redshift dependence b_{n} ∝ (1 + z)^{γ} with γ = 2.9 (so that P(k,z)b_{Lyα}(z)^{2} ∝ (1 + z)^{3.8}), and β_{Lyα} is assumed redshift independent. The function F_{NL} corrects for nonlinear effects at large k due to the isotropic enhancement of power due to nonlinear growth, the isotropic suppression of power due to gas pressure, and the suppression of power due to lineofsight nonlinear peculiar velocity and thermal broadening. We use the form given by Eq. (21) and Table 1 of McDonald (2003). The forms proposed by ArinyoiPrats et al. (2015) produce nearly identical results for the range of (r_{⊥},r_{∥}) used in this study.
Fig. 11
Correlation function for the metalfree mocks in four ranges of μ. The black points and curves correspond to mocks with Lyα absorption but without the addition of a quasar continuum. The red points and curves correspond to mocks with the addition of a continuum. The points correspond to stacks of 100 mocks and the light curves to individual mocks. The heavy curves correspond to the input model of Table 1 (after distortion by the matrix D_{AA′} (Eq. (10)) for the red curve). 
The last term in Eq. (21), G(k) = G_{∥}(k_{∥})G_{⊥}(k_{⊥}), accounts for the binning in (r_{⊥},r_{∥}) space which effectively averages the correlation function over a bin. (The large width of these bins renders unnecessary a term that accounts for the spectrometer resolution.) If the distribution of observed (r_{⊥},r_{∥}) were uniform in a bin, G(k) would be the Fourier transform of the bin. The distribution is approximately uniform in the radial direction, which implies G_{∥} = sinc^{2}(R_{∥}k_{∥}/ 2) where R_{∥} = 4 h^{1} Mpc (the bin width). In the perpendicular direction, the distribution is approximately proportional to r_{⊥} so G_{⊥} = [(J_{1}(k_{⊥}r_{+})−J_{1}(k_{⊥}r_{−})) /k_{⊥}Δr] ^{2} where r_{+} and r_{−} are the extrema of a given r_{⊥} bin and Δr = r_{+}−r_{−} is the bin width. The expression is bindependent and performing a Fourier transform on each bin would be too timeconsuming. We therefore replace G_{⊥} by sinc^{2}(k_{⊥}Δr) and evaluate the correlation function at the average r_{⊥} of the pairs in a bin. This procedure produces a sufficiently accurate correlation function.
We allow for fluctuations of ionizing UV radiation (Pontzen 2014; Gontcho A Gontcho et al. 2014) which lead to a scaledependent bias of Lyα absorption given by Eq. (12) of Gontcho A Gontcho et al. (2014). The effect of the fluctuations is to increase b_{Lyα} from its nominal value at small scale to a different value at large scale. The transition scale is determined by the UV photon meanfreepath which we set to a comoving value of 300 h^{1} Mpc (Rudie et al. 2013). We then fit for one parameter, b_{Γ} corresponding to the b_{Γ}(b_{s}−b_{a}) of Gontcho A Gontcho et al. (2014); it determines the change in b_{Lyα} between large and small scales. A second bias, , that determines the precise dependence of the bias on scale, is set to the nominal value of − 2/3 used by Gontcho A Gontcho et al. (2014).
In our previous investigations, ξ_{Lyα} was assumed to be a sufficiently accurate approximation for the correlation function. In this study, we include absorption by unidentified HighColumnDensity Systems (HCDs) and by metals. The total physical correlation function is then the sum of the auto and crosscorrelations of the various components: (23)Each component, n or m, has its own bias and redshiftspace distortion parameter.
The correlation function for HCDs and their crosscorrelation with Lyα absorption is derived from a power spectrum as in Eq. (21) but with different bias parameters (FontRibera & MiraldaEscudé 2012). While these absorbers are expected to trace the underlying density field their effect on the fluxtransmission field depends on whether they are identified and given the special treatment described in Sect. 2. If they are correctly identified with the total absorption region masked and the wings correctly modeled, they can be expected to have no significant effect on the field. Conversely, if they are not identified, the measured correlation function will be modified because their absorption is spread along the radial direction. We model this by multiplying the bias for unidentified HCDs, b_{HCD}, by a function of k_{∥}, F_{HCD}(k_{∥}), that reflects the typical column densities of the unidentified HCDs. Following the studies with the mock spectra described in Sect. 3, we adopted the form F_{HCD}(k_{∥}) = sinc(L_{HCD}k_{∥}), where L_{HCD} is a free parameter. The parameter β_{HCD} is poorly determined by the fit and we impose a Gaussian prior β_{HCD} = 0.5 ± 0.2.
Because there is little absorption by metals the treatment of metal components is simplified without the separation into peak and smooth components. The fiducial correlation function is directly used to calculate in real space both metalmetal and metalLyα correlations. The parameters (b,β) for each species (Table 3) are included by using a multipole expansion.
Absorption by metals is complicated by the fact that in the data analysis the Lyα transition is naturally assumed when transforming wavelength differences to position differences. Because of this procedure, for the crosscorrelation between absorption by transition n and transition m, the nominal (r_{⊥},r_{∥}) of the bin A is not equal to the true separation coordinates . The model correlation, for the bin A is thus evaluated at the true, rather than nominal, separation.
Because amplitudes for SiII and SiIII are mostly determined by the excess correlation at (r_{⊥} ~ 0,r_{∥} ≠ 0), the β for each metal is poorly determined. We therefore fix their value to β = 0.5 corresponding to host halos with bias of 2, the value found for DLAs (FontRibera & MiraldaEscudé 2012) which is also typical of starforming galaxies.
The bias, b_{CIV}, of foreground (z ~ 1.65) CIV absorption is poorly determined by the fit of the Lyαforest auto correlation. This is because the large wavelength difference between the Lyα and CIV transitions places the peak corresponding to vanishing metalLyα separation outside the forest. We have therefore measured b_{CIV} by correlating the flux in the CIV forest with quasars at z ~ 1.65. We compared these correlations with the analogous ones between the δ_{q}(λ) in the Lyα forest with z ~ 2.3 quasars. The resulting μaveraged correlation ratio is (24)Because the quasarCIV flux correlation is linear in b_{CIV} while the contribution to the flux autocorrelation is proportional to , this result implies a negligible (<1%) contamination of the Lyαforest autocorrelation by CIV autocorrelation.
Fig. 12
Correlation function for the stack of the 100 mocks in four ranges of μ. The red points represent the measured correlation function of the mocks with metals and the red curves shows the best fit. The black points indicate the correlation function after including, in addition to metals, HCDs unmasked for N_{HI}< 10^{20} cm^{2} and the black curve the best fit. The peak at r ~ 60 h^{1} Mpc due to SiIIa and SiIIb is apparent in the range μ> 0.95 but not in the range 0.8 <μ< 0.95. 
Weighted mean of fit parameters and mean uncertainties for six sets of 100 mocks of increasing realism: Lyα only; including a quasar continuum; including metal absorption (for Met1 or Met2); including HCDs (and Met1) and masking those with N_{HI}> 10^{20} cm^{2} or > 10^{21} cm^{2}.
To transform the ratio (24) into a constraint on b_{CIV} that is needed for the fits, we equate it to the appropriate function of bias parameters: (25)where ξ_{L} is the linear correlation function at 10 h^{1} Mpc. To determine b_{CIV}(1.65) we use the values (β_{Lyα},β_{q},β_{CIV}) = (1.4, 0.27, 0.5), a growth factor ξ_{L}(1.65) /ξ_{L}(2.3) = 1.5, and the redshift evolution of the quasar bias from Croom et al. (2005), b_{q}(1.65) ] /b_{q}(2.3) = 0.70. The resulting value of b_{CIV}(1.65) /b_{Lyα}(2.3) is multiplied by the measured value of b_{Lyα} (Table 5) and then used as a Gaussian prior in the fits: b_{CIV}(2.3) = −0.0224 ± 0.0112 (and b_{CIV}< 0), where we have conservatively doubled the uncertainty in (24).
Fig. 13
Measured α_{⊥} and α_{∥} (left) and β_{Lyα} and b_{Lyα}(1 + β_{Lyα}) (right) for the 100 mock catalogs including HCDs with N_{HI}> 10^{20} cm^{2} masked. There are four outliers on the left plot and two on the right. The horizontal and vertical blue lines show the weighted means of the distributions, also given in Table 4. The distribution and mean values of (α_{⊥},α_{∥}) indicates no significant bias in the reconstruction of the BAO peakposition parameters. 
To test for systematic errors, the model in Eq. (18) includes the broadband term B_{A}, a smoothly varying function of (r_{⊥},r_{∥}): (26)where the L_{j} are Legendre polynomials. In our previous studies, B_{A} had a central role because we did not attempt to model distortions due to continuum fitting as we now do with D_{AA′}. We include B_{A} now as an option to study how the fit position of the BAO peak depends on the assumed form of the smooth component of the correlation function. By including this function, we can assure that the measurement of the BAO parameters is due only to the peak position, and is not significantly influenced by the assumed form of the smooth component.
6.2. Fits with the mock data sets
The 100 mock sets were used to test the fitting procedure and to search for possible biases in the fit BAO parameters. Studies were made on individual mock sets and for the stack of the 100 sets. The former test our ability to reliably extract (α_{⊥},α_{∥}) in data sets similar to the BOSS observations. The latter test the accuracy with which we model the analysis procedure including the effects of weighting, (r_{⊥},r_{∥}) binning, and continuum fitting.
Four types of mocks were produced and analyzed. The results are summarized in Table 4 and Figs. 11 and 12. The four types and the analysis procedure are as follows:

1.
Only Lyα absorption: in the fits the distortionmatrix is set to the identity matrix and only the parameters(α_{⊥},α_{∥},b_{Lyα},β_{Lyα}) are fit.

2.
A quasar continuum is added: the fits use the distortion matrix and the parameters (α_{⊥},α_{∥},b_{Lyα},β_{Lyα}) are fit.

3.
Metal absorption is added: the additional parameters are the (b_{i},β_{i}) for each metal transition, i.

4.
High column density systems are added: three additional parameters (b_{HCD},β_{HCD},L_{HCD}) are fit.
Figure 11 shows the correlation function in four ranges of μ. The black points are the stack of the measured correlation function for type 1 mocks (only Lyα absorption). The red points are the stacked measurements for type 2 mocks (quasar continuum added). The difference between the black and red points shows the effect of continuumfitting. The heavy curves are the best fit model for the cosmology used to generate the mocks. The good agreement between the red curve and points indicates the accuracy of the modeling of the continuumfitting distortion.
Fig. 14
Measured correlation function in four ranges of μ. The most radial range (topleft μ> 0.95) has, in addition to the BAO peak at ~100 h^{1} Mpc, a peak at ~60 h^{1} Mpc due to correlated absorption by Lyα and SiIII(119.0,119.3) at the same physical position. In the next radial bin (topright, 0.8 <μ< 0.95) only the BAO peak is visible. The lines show fits including successively Lyα and metal absorption, unidentified HCDs absorption, UV flux fluctuations, and a (i_{min},i_{max},j_{max}) = (0,2,6) broadband (BB). 
Fits with type 3 mocks (including metals) were used to search for possible influences of metal absorption on derived BAO parameters. Metal absorption is mostly confined to the small r_{⊥} region and this is seen in the two top panels of Fig. 12. The most radial sector, μ> 0.95, displays the effects of metal absorption, most prominently the peak due to SiIIa,b absorption around r ~ 60 h^{1} Mpc. The peak is not visible in the next sector, 0.8 <μ< 0.95.
Two sets of fits for type 4 mocks (HCDs) were also performed. In one, HCDs with N_{HI}< 10^{20} cm^{2} were not given the DLA treatment described in Sect. 4.1. This cut approximates the threshold for DLA detection either automatically or visually. In the second set, only HCDs with N_{HI}> 10^{21} cm^{2} were treated specially, testing our sensitivity to an underestimation of the efficiency to flag DLAs.
Table 4 lists the mean fit values for the 100 mocks and for the 4 types. Figure 13 presents the corresponding distributions for type 4 mocks. For all types, the mean of the bestfit values of (α_{⊥},α_{∥}) is consistent with unity at the subpercent level. Since the precision of our measurement with the data is of order three percent, this means that the mocks do not indicate any significant bias in the measurement of the BAO peak position.
6.3. Fits of the data with physical correlation function
Results of fits of the data assuming increasingly complicated forms for ξ^{ph}(r_{⊥},r_{∥}): Lyα absorption only; including metal absorption; including unidentified HCDs; and including fluctuations in the ionizing UV flux (i.e. the complete physical model whose complete set of parameters is given in Table 3).
The data were fit over the range 10 <r< 180 h^{1} Mpc under various hypotheses of increasing complexity concerning the form of ξ^{ph}(r_{⊥},r_{∥}). The results are summarized in Table 5. The first fit assumes only Lyα absorption, the second includes absorption by metals, and the third absorption by unidentified HCDs. The fourth fit also assumes fluctuation in the flux of ionizing UV radiation. Inclusion of each effect yields a significantly better χ^{2}; however, the values of (α_{⊥},α_{∥}) do not change significantly with each step. The data and fit correlation functions are presented for the four angular ranges in Fig. 14.
In what follows, we focus our attention on the fourth fit that will be referred to as the “complete physical model”. The BAO peakposition parameters for this model are where we have used the Planck values from Table 1 that were used in the (angle, redshift) to (r_{⊥},r_{∥}) transformation. The one, two and three standard deviation contours (corresponding to Δχ^{2} = 2.3,6.18,11.83) are shown as the red contours in Fig. 15. Our result is within one standard deviation of the prediction of flatΛCDM model of Table 1 that is in agreement with the Planck data.
In the next section we study the robustness of the determination of the BAOpeak position parameters by adding powerlaw broadband terms to the complete physical model. We can also estimate the significance of the BAO peak by modifying this model by setting A_{peak} = 0. This degrades the quality of the fit by Δχ^{2} = 27.7 corresponding to a BAO peak detection of 5.2σ. Adding broadband terms reduces this to Δχ^{2} = 21.8 (4.7σ).
Other than (α_{⊥},α_{∥}), the parameters of the model listed in Table 3 have no direct cosmological implications. However, the fact that the bestfit values are physically reasonable reinforces our confidence in the model and, therefore, on the cosmological conclusions. Most importantly, the bestfit values of (b_{Lyα},β_{Lyα}) in Table 5 are both relatively stable and near the expected values. For example, the model of ArinyoiPrats et al. (2015) predicts β_{Lyα} ~ 1.4 and  b_{Lyα}  (1 + β_{Lyα}) in the range 0.25 to 0.33. In Appendix B we discuss in more detail the comparison of our measurement of (b_{Lyα},β_{Lyα}) with theory and with previous measurements.
The bestfit values of most of the other parameters in Table 3 are near those found in the mock fits. For the metal biases, this only means that we have placed a reasonable amount of metal absorption in the mock spectra. Of more significance is the value of L_{HCD} = (24.3 ± 1.13) h^{1} Mpc which is, as expected, roughly the width of a DLA of N_{HI} = 10^{20} cm^{2}, corresponding to our estimated threshold for good efficiency in identifying DLAs. In the mocks, bestfit values of L_{HCD} when masking only those DLAs with N_{HI}> 10^{20} cm^{2} cluster near L_{HCD} ~ 25 in about half the fits. For the other half, the fits find L_{HCD} ~ 0, corresponding to the addition of a second Lyα component but with a different β. This indicates that more than one way exists to model sufficiently well the correlation function with unidentified HCDs. The bias parameter found in the data, b_{HCD} = −0.0288 ± 0.0043, is about a factor of 7 higher than that found in those mocks that find L_{HCD} ~ 25. This indicates that there are unidentified HCDs in the data with N_{HI}> 10^{20} cm^{2} and or that the HCDs that we have placed in the mocks are less biased than in reality.
We have not generated mocks with fluctuations in the flux of ionizing radiation. However, we note that the bias parameter b_{Γ} = 0.1125 ± 0.0511 is near the value of 0.13 suggested by Gontcho A Gontcho et al. (2014).
6.4. Fits of the data with powerlaw broadbands added
In addition to statistical uncertainties, it is important to determine to what extent the derived values of (α_{⊥},α_{∥}) are influenced by the assumed form of ξ^{ph}. The amplitude and form of the metal components are well understood and constrained by the lowr_{⊥} bins. However, the HCD and UV parameters have a more diffuse effect on the correlation function and are not strongly constrained by independent data. For example, the UV background fluctuation model we have used includes only fluctuations that follow the underlying density field. Additional correlations (Gontcho A Gontcho et al. 2014) due to the discrete nature of sources of UV radiation would modify the absorption in neighboring forests since they are influenced by the same discrete sources. This “shotnoise” term would give an additive, isotropic and positive contribution to the Lyα autocorrelation function. Its shape and amplitude as a function of comoving separation are however poorly known as they depend on the quasar luminosity distribution and variability.
To study the effects of possible uncertainties in the form of the nonBAO component, , we have performed fits adding to the complete physical model broadband functions in the form of power laws (26) defined by (i_{min},i_{max},j_{max}). With the addition of new components of , the connection between the BAO peak amplitude and the smooth components is lost so we also performed fits releasing the constraint A_{peak} = 1.
The parameters of slowly varying broadband forms are necessarily somewhat degenerate with the parameters of the slowly varying part of the physical model: (b,β) and the parameters describing HCDs and UV fluctuations. As such, fits with powerlaw broadband terms can produce unphysical bestfit values of these parameters. Conversely, we expect that a small number of powerlaw terms will not significantly affect the BAOpeak parameters. Fits with the mock spectra indicate that this is the case if the powerlaw broadbands are restricted to the ranges − 2 ≤ i_{max} ≤ 3, (i_{max}−i_{min}) ≤ 3 and j_{max} ≤ 6. Of course in the case of the mocks we know precisely the form of the physical correlation function while for the data we do not. Shifts of (α_{⊥},α_{∥}) in the data when using powerlaw terms could indicate that the underlying physical correlation function has not been modeled with sufficient precision.
Bestfit values of (α_{⊥},α_{∥}) for fits including metals, HCDs and UV fluctuations, with and without powerlaw broadbands (BBs) of the form (i_{min},i_{max},j_{max}) = (0,2,6).
We want to ensure that the powerlaw terms model variations of the slowlyvarying part of the correlation function under the BAO peak. We therefore perform these fits only over the restricted range 40 <r< 180 h^{1} Mpc, avoiding undue influence of the 10 <r< 40 h^{1} Mpc range on the amplitudes of the power laws. Using this range, the fit with the complete physical model without powerlaw broadbands yields BAO parameters indistinguishable from those found over the full range 10 <r< 180 h^{1} Mpc (Table 5).
Fig. 15
Contours for (α_{⊥},α_{∥}) at the (68.3, 95.45, 99.73)% CL (solid, dashed, dotted). The red contours are for the physical model with HCDs and UV fluctuations from Tables 3 and 5. The blue (green) contours are for the (i_{min},i_{max},j_{max}) = (0,2,6) powerlaw broadband with (without) priors on the parameters of the physical correlation function. 
To choose appropriate values of (i_{min},i_{max},j_{max}), we note that apart from the BAO peak, r^{2}ξ does not vary wildly over the range 40 <r< 180 h^{1} Mpc. A reasonable choice is therefore (i_{min},i_{max}) = (0,2) corresponding to a parabola and which amounts to freeing the value and first two derivatives of r^{2}ξ_{smooth} underneath the BAO peak. A reasonable choice for j_{max} is j_{max} = 6 corresponding to approximately independent broadbands in each of the four angular ranges in Fig. 14. The bestfit correlation function for (i_{min},i_{max},j_{max}) = (0,2,6) and A_{peak} free is shown in the figure and the BAOpeak parameters on the third line of Table 6. The bestfit value of α_{∥} is ~1% (i.e. ~0.3σ) lower than the value for the complete physical model and uncertainty is unchanged. The value of α_{⊥} is about 3% higher than the standard value with the error increasing slightly from ~0.055 to ~0.075. The one and twostandard deviation contours are shown in green in Fig. 15. We emphasize that these results are obtained without placing any priors on the physical parameters, other than those already used in the standard fits.
We have experimented with other combinations of (i_{min},i_{max},j_{max}) in the range j_{max} ≤ 6, − 2 <i_{min}< 0, − 1 <i_{max}< 3. While most combinations yield values of (α_{⊥},α_{∥}) that are close to the nominal values, there are cases where there are significant changes. These cases always correspond to fits yielding unphysical values for (b,β) of one or more components. In all cases, the large changes in (α_{⊥},α_{∥}) can be avoided by placing a “physical priors” on the (b,β) that requires them to take values within two standard deviations of those found in the standard fits over the range 10 <r< 180 h^{1} Mpc and imposing A_{peak} = 1. The results for the (0,2,6) model are listed oneline two of Table 6 and the contours shown in blue in Fig. 15.
6.5. Summary of data fits
Table 6 and Fig. 15 summarizes the BAO parameters found for fits with our most complete physical model of the correlation function, with and without powerlaw broadbands. The results of the fits with physical priors imply that the uncertainty in the form of does not translate into a significant uncertainty in (α_{⊥},α_{∥}). Because of the priors placed on the physical parameters, this conclusion follows only if our basic physical model is not too far removed from reality. Without such priors, the uncertainty of α_{⊥} can be significantly increased (line 3 of Table 6).
The limits on α_{∥} are relatively stable with the inclusion of broadbands. A quantity that is more robust than either α_{∥} or α_{⊥} is with the exponent γ chosen to minimize the variance. For the observed correlation and variances of (α_{⊥},α_{∥}), one calculates γ ~ 0.7 with (29)where the uncertainty is chosen to cover the onestandarddeviations of the three fits in Table 6. Using the Planck model values from Table 1, we deduce (30)This result is largely independent of the addition of powerlaw broadbands. It contains most of the cosmological constraints of our measurement if used in combination with constraints from the forestquasar crosscorrelation, since the latter has more constraining power on D_{M}/r_{d} than the autocorrelation result.
6.6. Comparison with previous results
We can compare the results (27) and (28) with those of Delubac et al. (2015): D_{H}(2.34) /r_{d} = 9.18 ± 0.28, D_{M}/r_{d} = 37.67 ± 2.17 and . The change in D_{H}/r_{d} is ~ 0.5σ. To determine whether this change is due to the change in analysis procedure or in the data set we performed fits using only data used in the DR11 analysis. The results for a HCDS,UV analysis are D_{H}(2.33) /r_{d} = 9.14 ± 0.28 and D_{M}(2.34) /r_{d} = 37.54 ± 2.17, within 0.15σ and 0.06σ of the results of Delubac et al. (2015). This result implies that most of the change is due to the change in the data set.
The 0.5σ change in D_{H}/r_{d} is typical of the variations we can expect when increasing the sample size from DR11 to DR12. We confirmed this expectations by using jackknife samples of DR12 to simulate DR11 data sets and observing changes in D_{H}/r_{d}. Roughly 30% of the jackknife samples had a change as large as that observed for DR11 and we therefore conclude that the DR12DR11 difference is primarily statistical.
7. Systematic Errors in the BAOpeak position
Systematic errors in (α_{⊥},α_{∥}) can result if the measured correlations function is not sufficiently well modeled by the fitting function (18). Such effects can be introduced by the pipeline and analysis procedures, by improper modeling for the fluxtransmission field, and by astrophysical absorption not related to cosmological largescale structure.
Spurious correlations introduced by the pipeline were discussed in Sect. 5. The most important correlations were confined to the lowest r_{∥} bin with negligible effects on the fitting of the correlation function. No evidence was found for any additional instrumental correlations that could have an influence of the measurement of (α_{⊥},α_{∥}). Since our model of pipelineinduced correlations introduces no scales near the BAO scale, it is not surprising that it predicts negligible effects on the BAO scale. Upper limits on the pipelineinduced shifts in (α_{⊥},α_{∥}) are given in Table 7.
Biases could also be introduced at the analysis level, for example through the small correlations between weights and the measured δ_{q}(λ). These correlations are estimated to have a negligible impact on the BAO parameters. If they did not, they would have produced biased estimates of (α_{⊥},α_{∥}) in the mock data sets, while no such biases are seen.
The second type of systematic error would result from incorrect modeling of the various contributions to the physical correlations. The physical correlations are parameterized by the (b_{i},β_{i}) for each transition. These parameters are marginalized over and their uncertainties therefore contribute to the reported statistical errors on (α_{⊥},α_{∥}).
Modeling of the contribution of unidentified HCDs and of UV fluctuations also contain freeparameters that are marginalized over. We have verified that the bestfit BAO parameters are insensitive to assumption about nonlinearities and redshift evolution of the bias parameters. We have also investigated possible influence of foreground absorbers in addition to CIV. Finally, we have included powerlaw broadband parameters to place limits on our sensitivity to unmodeled slowlyvarying contributions.
In our physical model of the correlation function we have ignored two effects that could potentially bias our results at some level. First, we have assumed that the lines of sight sample random positions of the Universe, while in a real survey we only have information in pixels in front of a quasar. This makes our measured correlation sensitive to the threepoint function of quasarlyalya, but given that this is averaged over pixel pairs at different separations from the quasars we expect this contamination to be very smooth on BAO scales. The second effect is that the transmitted flux fraction could be affected by the relative velocity between baryons and dark matter in the early Universe (Tseliakhovich & Hirata 2010). These relative velocities are quite small in the redshift range of interest, but it is possible that the different velocities at early times left an imprint in the distribution of neutral hydrogen even at low redshift, similar to what has been suggested in the case of galaxy clustering (Dalal et al. 2010; Yoo et al. 2011; Blazek et al. 2016). In order to study this effect in the Lyman alpha forest we would require detailed hydrodynamical simulations, and this is beyond the scope of this paper. Given that the density of quasars should have a different dependency on the relative velocity, differences in the BAO scale measured from the autocorrelation and from the crosscorrelation with quasars (FontRibera et al. 2014) could be used to constraint the amplitude of this effect.
The third type of systematic error would be due to correlations of noncosmological origin. One possibility is atomic or molecular absorption in the Galactic interstellar medium (ISM). The average ISM absorption is included in the fit mean flux transmission, or the subsequent subtraction of the average δ_{q}(λ). However, residuals correlated across the sky might still be present. Current studies show that the ISM absorption is highly correlated with dust column density. We hence use the estimated MilkyWay dust extinction as a proxy for potential ISM absorption and fit for the correlation of δ_{q}(λ) with this dust extinction parameter. The resulting correlation is shown in Fig. 16. The Ca H and K lines are the only known ISM lines to appear significantly. Therefore ISM absorption has a negligible impact on the Lyα correlation function. We have verified this conclusion by computing the autocorrelation function of the quantity E(B−V) × dδ_{q}(λ)/dE(B−V).
Fig. 16
Derivative of the normalized Lyαforest fluxdecrement, δ_{LyaF}, with respect to the Milky Way dust column density normalized to its dispersion. The dashed curves represent the 0, 1, 2 and 3σ uncertainties. The red vertical lines mark the positions of the calcium H and K lines. 
A final source of correlations that we have considered would originate in our use of a unique quasarcontinuum template which we adapt to individual quasars by the (a_{q},b_{q}) parameters of Eq. (3). Templates could be improved by making them depend on a nonforest variable, like quasar redshift or CIV emissionline strength. These improved templates could change the derived flux correlation function if the nonforest variables of the quasars are themselves correlated at large distances. This is not expected to be the case of CIV line strength since this strength is believed to be due to local conditions and therefore not correlated at large distances. We have tested this by using two templates, one for strong CIV emitters (relative to sidebands) and one for weak emitters. This slightly improves the pixel variance (by 3%) but makes no significant change in the measured correlation function and less than 0.1% difference in the derived (α_{⊥},α_{∥}).
A continuum that evolves with redshift could be more problematic because the correlation function measurement necessarily requires quasar pairs that are near in redshift. We searched for the effect of such evolution by using a quasar continuum template that evolves linearly in redshift. The results do not differ significantly from those for the redshiftindependent template.
We have searched for unsuspected systematics by dividing the data into two redshift bins, z< 2.295 and z> 2.295. The 95%CL ranges on α_{∥} in the high and low redshift ranges are 0.90 <α_{∥}< 1.094 and 1.02 <α_{∥}< 1.15. For α_{⊥} the ranges are 0.74 <α_{⊥}< 1.5 and 0.86 <α_{⊥}< 1.01. As with the DR11 results of Delubac et al. (2015), the results are consistent with the expected absence of evolution at the twostandarddeviation level.
8. Cosmological implications
Our measurements of D_{M}/r_{d} and D_{H}/r_{d} at z = 2.33 are in good agreement with the values predicted by the flatΛCDM model consistent with the CMB data. The precision is, however, less than that of the BOSS galaxy BAO data (Alam et al. 2016), and for simple cosmological models, parameter constraints based on combinations of CMB and BAO data are not strongly changed by our results. Our constraints however, can provide tests of more general models with complicated expansion dynamics (Aubourg et al. 2015).
When combined with other BAO data, our measurement provides interesting results using only lowredshift data. For example, the lowredshift expansion history can be directly measured through the measurements of D_{H}(z) /r_{d}. The Hubble distance (27) can be written as (31)Figure 17 displays this result along with those of Alam et al. (2016), Beutler et al. (2011), Ross et al. (2015) and FontRibera et al. (2014). The expected phases of acceleration at z< 0.7 and deceleration at z> 0.7 are striking.
Fig. 17
Comoving expansion rate, r_{d}H(z)/(1 + z), measured with BAO, with r_{d} in units of 147.3 Mpc. The red square is the present measurement at z = 2.33. Also shown are measurements using the Lyαquasar cross correlation at z = 2.36 (FontRibera et al. 2014), and galaxy correlations at 0.35 <z< 0.65 (Alam et al. 2016), at z = 0.15 (Ross et al. 2015), and at z = 0.11 (Beutler et al. 2011). The points at z = 0.11 and z = 0.15 use the SNIa measurement of q_{0} (Betoule et al. 2014) to convert the measured to D_{H}(z). The black line is the prediction of the Planck flat ΛCDM model (Table 1). 
Combining D_{H}/r_{d} and D_{M}/r_{d} measurements, the cosmological parameters can be measured independently of CMB data. Nonflat ΛCDM models (oΛCDM ) have three parameters (Ω_{M},Ω_{Λ},H_{0}r_{d}) that determine D_{M}/r_{d} and D_{H}/r_{d} and constraints on (Ω_{M},Ω_{Λ}) can be derived by marginalizing over H_{0}r_{d}. This procedure is equivalent to the use of SNIa to measure (Ω_{M},Ω_{Λ}) by marginalizing over the parameters that determine supernova luminosities. Figure 18 shows the constraints derived using various combinations of BAO results. By themselves, the lowredshift BAO data provide no strong evidence for Ω_{Λ}> 0. The inclusion of the Lyαforest data results in ~10% precision on (Ω_{M},Ω_{Λ}): (32)consistent with spatial flatness: (33)The fit has χ^{2} = 12.27 for 12 data points and three parameters. The values of (Ω_{M},Ω_{Λ}) agree well with the flatΛCDM values found using CMB anisotropies (Ω_{M},Ω_{Λ}) = (0.315,0.685) (Planck Collaboration XIII 2016).
9. Conclusions
In this paper, we present an analysis of the complete SDSS3 (DR12) set of Lyα forests. The position of the BAO peak in the flux correlation function is in good agreement with that expected in flatΛCDM models favored by CMB anisotropies.
Our analysis represents a significant improvement over our previous work, both in the understanding of instrumental effects on the flux transmission and in the data analysis. The modeling of the distortion of the correlation function due to continuum fitting allows fitting the entire correlation function using a physical model without resorting to arbitrary “broadband terms”. In the future, this approach may allow alternative ways to study the cosmological parameters, for instance by eliminating the need to separate the correlation function into “smooth” and “peak” components.
Now that we have direct access to the full correlation function, further improvements in the analysis will profit from better modeling of absorption mechanisms. In particular, it will be important to fully use external constraints on HCD and metal absorption.
Fig. 18
One and twostandard deviation constraints on the oΛCDM parameters (Ω_{M},Ω_{Λ}) using BAO measurements of D_{M}(z) /r_{d} and D_{H}(z) /r_{d} without CMB measurements, i.e. without imposing the CMB value of r_{d}. The red contours use the work presented here (line 1 of Table 6), the z< 0.8 galaxy data (Alam et al. 2016; Beutler et al. 2011; Ross et al. 2015) and the quasarforest crosscorrelation (FontRibera et al. 2014). The gray contours exclude the quasarforest crosscorrelation. The blue contours use only the z< 0.8 galaxy data and the green contours show the constraints derived from the SNIa Hubble diagram (Betoule et al. 2014). The black dot is the position of the flatΛCDM model that describes CMB data (Planck Collaboration XIII 2016). 
The improvements in modeling and data analysis that we have implemented here will clearly be important for future projects with reduced statistical errors. For example, the DESI project (Aghamousa et al. 2016) will have up to 15 times the number of forestpixel pairs resulting in a significant decrease of statistical errors of the correlation function. The understanding the physics of the broadband correlation function may be the limiting factor in the extraction of BAO parameters.
The cosmological information from the Lyαforest autocorrelation function studied here is complemented by the study of the crosscorrelation between the forest and neighboring quasars. Results and implications of a new measurement of this function will be presented in a forthcoming publication (du Mas des Bourboux et al., in prep.).
Acknowledgments
Funding for SDSSIII 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 SDSSIII web site is http://www.sdss3.org/. SDSSIII is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSSIII Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University. The French Participation Group of SDSSIII was supported by the Agence Nationale de la Recherche under contracts ANR08BLAN0222 and ANR12BS05001501. NB and JG acknowledge support from the French Programme National Cosmologie et Galaxies. M.B., M.M.P. and I.P. were supported by the A*MIDEX project (ANR 11IDEX000102) funded by the “Investissements d’Avenir” French Government program, managed by the French National Research Agency (ANR), and by ANR under contract ANR14ACHN0021.The work of J.B. and K.D. was supported in part by the U.S. Department of Energy, Office of Science, Office of High Energy Physics, under Award Number DESC0009959. N.P.R acknowledges support from the S.T.F.C. and the Ernest Rutherford Fellowship scheme.
References
 Aghamousa, A., Aguilar, J., et al. (DESI Collaboration) 2016, ArXiv eprints [arXiv:1611.00036] [Google Scholar]
 Alam, S., Albareti, F. D., Allende Prieto, C., et al. 2015, ApJS, 219, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Alam, S., Ata, M., Bailey, S., et al. 2016, MNRAS, submitted [arXiv:1607.03155] [Google Scholar]
 Albareti, F. D., Allende Prieto, C., et al. (SDSS Collaboration) 2016, ApJS, submitted [arXiv:1608.02013] [Google Scholar]
 Anderson, L., Aubourg, E., Bailey, S., et al. 2012, MNRAS, 427, 3435 [NASA ADS] [CrossRef] [Google Scholar]
 Anderson, L., Aubourg, É., Bailey, S., et al. 2014a, MNRAS, 441, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Anderson, L., Aubourg, E., Bailey, S., et al. 2014b, MNRAS, 439, 83 [NASA ADS] [CrossRef] [Google Scholar]
 ArinyoiPrats, A., MiraldaEscudé, J., Viel, M., & Cen, R. 2015, J. Cosmol. Astropart. Phys., 12, 017 [NASA ADS] [CrossRef] [Google Scholar]
 Aubourg, É., Bailey, S., Bautista, J. E., et al. 2015, Phys. Rev. D, 92, 123516 [NASA ADS] [CrossRef] [Google Scholar]
 Bautista, J. E., Bailey, S., FontRibera, A., et al. 2015, J. Cosmol. Astropart. Phys., 5, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Becker, R. H., White, R. L., & Helfand, D. J. 1995, ApJ, 450, 559 [NASA ADS] [CrossRef] [Google Scholar]
 Betoule, M., Kessler, R., Guy, J., et al. 2014, A&A, 568, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beutler, F., Blake, C., Colless, M., et al. 2011, MNRAS, 416, 3017 [NASA ADS] [CrossRef] [Google Scholar]
 Blake, C., Davis, T., Poole, G. B., et al. 2011, MNRAS, 415, 2892 [NASA ADS] [CrossRef] [Google Scholar]
 Blazek, J. A., McEwen, J. E., & Hirata, C. M. 2016, Phys. Rev. Lett., 116, 121303 [NASA ADS] [CrossRef] [Google Scholar]
 Blomqvist, M., Kirkby, D., Bautista, J. E., et al. 2015, J. Cosmol. Astropart. Phys., 11, 034 [NASA ADS] [CrossRef] [Google Scholar]
 Bolton, A. S., Schlegel, D. J., Aubourg, É., et al. 2012, AJ, 144, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Bovy, J., Hennawi, J. F., Hogg, D. W., et al. 2011, ApJ, 729, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Busca, N. G., Delubac, T., Rich, J., et al. 2013, A&A, 552, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chuang, C.H., & Wang, Y. 2012, MNRAS, 426, 226 [NASA ADS] [CrossRef] [Google Scholar]
 Cole, S., Percival, W. J., Peacock, J. A., et al. 2005, MNRAS, 362, 505 [NASA ADS] [CrossRef] [Google Scholar]
 Croft, R. A. C., MiraldaEscudé, J., Zheng, Z., et al. 2016, MNRAS, 457, 3541 [NASA ADS] [CrossRef] [Google Scholar]
 Croom, S. M., Boyle, B. J., Shanks, T., et al. 2005, MNRAS, 356, 415 [NASA ADS] [CrossRef] [Google Scholar]
 Dalal, N., Pen, U.L., & Seljak, U. 2010, J. Cosmol. Astropart. Phys., 11, 007 [NASA ADS] [CrossRef] [Google Scholar]
 Dawson, K. S., Schlegel, D. J., Ahn, C. P., et al. 2013, AJ, 145, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Delubac, T., Bautista, J. E., Busca, N. G., et al. 2015, A&A, 574, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eisenstein, D. J., Zehavi, I., Hogg, D. W., et al. 2005, ApJ, 633, 560 [NASA ADS] [CrossRef] [Google Scholar]
 Eisenstein, D. J., Weinberg, D. H., Agol, E., et al. 2011, AJ, 142, 72 [NASA ADS] [CrossRef] [Google Scholar]
 FontRibera, A., & MiraldaEscudé, J. 2012, J. Cosmol. Astropart. Phys., 7, 28 [NASA ADS] [CrossRef] [Google Scholar]
 FontRibera, A., McDonald, P., & MiraldaEscudé, J. 2012, J. Cosmol. Astropart. Phys., 1, 1 [NASA ADS] [CrossRef] [Google Scholar]
 FontRibera, A., Arnau, E., MiraldaEscudé, J., et al. 2013, J. Cosmol. Astropart. Phys., 5, 18 [NASA ADS] [CrossRef] [Google Scholar]
 FontRibera, A., Kirkby, D., Busca, N., et al. 2014, J. Cosmol. Astropart. Phys., 5, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Gontcho A Gontcho, S., MiraldaEscudé, J., & Busca, N. G. 2014, MNRAS, 442, 187 [NASA ADS] [CrossRef] [Google Scholar]
 Gunn, J. E., Carr, M., Rockosi, C., et al. 1998, AJ, 116, 3040 [NASA ADS] [CrossRef] [Google Scholar]
 Gunn, J. E., Siegmund, W. A., Mannery, E. J., et al. 2006, AJ, 131, 2332 [NASA ADS] [CrossRef] [Google Scholar]
 Kirkby, D., Margala, D., Slosar, A., et al. 2013, J. Cosmol. Astropart. Phys., 3, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Kirkpatrick, J. A., Schlegel, D. J., Ross, N. P., et al. 2011, ApJ, 743, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Lawrence, A., Warren, S. J., Almaini, O., et al. 2007, MNRAS, 379, 1599 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473 [NASA ADS] [CrossRef] [Google Scholar]
 Margala, D., Kirkby, D., Dawson, K., et al. 2016, ApJ, 831, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Martin, D. C., Fanson, J., Schiminovich, D., et al. 2005, ApJ, 619, L1 [NASA ADS] [CrossRef] [Google Scholar]
 McDonald, P. 2003, ApJ, 585, 34 [NASA ADS] [CrossRef] [Google Scholar]
 McDonald, P., & Eisenstein, D. J. 2007, Phys. Rev. D, 76, 063009 [NASA ADS] [CrossRef] [Google Scholar]
 Mehta, K. T., Cuesta, A. J., Xu, X., Eisenstein, D. J., & Padmanabhan, N. 2012, MNRAS, 427, 2168 [NASA ADS] [CrossRef] [Google Scholar]
 Noterdaeme, P., Petitjean, P., Carithers, W. C., et al. 2012, A&A, 547, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Padmanabhan, N., Xu, X., Eisenstein, D. J., et al. 2012, MNRAS, 427, 2132 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Pâris, I., Petitjean, P., Aubourg, É., et al. 2012, A&A, 548, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pâris, I., Petitjean, P., Aubourg, É., et al. 2014, A&A, 563, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pâris, I., Petitjean, P., Ross, N. P., et al. 2017, A&A, 597, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Peebles, P. J. E., & Yu, J. T. 1970, ApJ, 162, 815 [NASA ADS] [CrossRef] [Google Scholar]
 Percival, W. J., Cole, S., Eisenstein, D. J., et al. 2007, MNRAS, 381, 1053 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Percival, W. J., Reid, B. A., Eisenstein, D. J., et al. 2010, MNRAS, 401, 2148 [NASA ADS] [CrossRef] [Google Scholar]
 Pieri, M. M., Frank, S., Weinberg, D. H., Mathur, S., & York, D. G. 2010, ApJ, 724, L69 [NASA ADS] [CrossRef] [Google Scholar]
 Pieri, M. M., Mortonson, M. J., Frank, S., et al. 2014, MNRAS, 441, 1718 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration. XVI. 2014, A&A, 571, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration. XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pontzen, A. 2014, Phys. Rev. D, 89, 083010 [NASA ADS] [CrossRef] [Google Scholar]
 Prochaska, J. X., HerbertFort, S., & Wolfe, A. M. 2005, ApJ, 635, 123 [NASA ADS] [CrossRef] [Google Scholar]
 Ross, N. P., Myers, A. D., Sheldon, E. S., et al. 2012, ApJS, 199, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Ross, A. J., Samushia, L., Howlett, C., et al. 2015, MNRAS, 449, 835 [NASA ADS] [CrossRef] [Google Scholar]
 Rudie, G. C., Steidel, C. C., Shapley, A. E., & Pettini, M. 2013, ApJ, 769, 146 [NASA ADS] [CrossRef] [Google Scholar]
 Slosar, A., FontRibera, A., Pieri, M. M., et al. 2011, J. Cosmol. Astropart. Phys., 9, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Slosar, A., Iršič, V., Kirkby, D., et al. 2013, J. Cosmol. Astropart. Phys., 4, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Smee, S. A., Gunn, J. E., Uomoto, A., et al. 2013, AJ, 146, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Sunyaev, R. A., & Zeldovich, Y. B. 1970, Ap&SS, 7, 3 [NASA ADS] [Google Scholar]
 Suzuki, N. 2006, ApJS, 163, 110 [NASA ADS] [CrossRef] [Google Scholar]
 Suzuki, N., Tytler, D., Kirkman, D., O’Meara, J. M., & Lubin, D. 2005, ApJ, 618, 592 [NASA ADS] [CrossRef] [Google Scholar]
 Tseliakhovich, D., & Hirata, C. 2010, Phys. Rev. D, 82, 083520 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, X., Cuesta, A. J., Padmanabhan, N., Eisenstein, D. J., & McBride, C. K. 2013, MNRAS, 431, 2834 [NASA ADS] [CrossRef] [Google Scholar]
 Yèche, C., Petitjean, P., Rich, J., et al. 2010, A&A, 523, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yoo, J., Dalal, N., & Seljak, U. 2011, J. Cosmol. Astropart. Phys., 7, 018 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Modifications applied to the BOSS spectral extraction pipeline
Fig. A.1
Comparison of flux measurement uncertainties between the DR12 data reduction (in blue) and this new extraction (in red). This analysis is based on the data from camera b1, plate 7339, where many observations of the same targets were performed. The top panel shows the empirical rms (circles) among the various observations, and the corresponding predicted uncertainty (lines) of the number of electrons per CCD row (N_{e}) as a function of N_{e}. The second panel is the same quantity normalized by the predicted uncertainty of DR12. The bottom panel displays the fraction of 3σ outliers that were discarded in the two other panels. As expected, there is a larger dispersion with the new, less optimal reduction at high flux. The larger dispersion at low flux is not easy to understand; it is probably a consequence of the biased estimator in DR12, where negative statistical fluctuations are given a high weight. 
Each of the two BOSS spectrograph slitheads consists of 500 fibers organized in 25 aligned blocks (or bundles) of 20 fibers (see Smee et al. 2013, for a description of the BOSS spectrographs). Fibers in a bundle are separated by 260 μm whereas the separation between fibers of adjacent bundles is 624 μm. This design allows one to treat independently each fiber bundle in the data reduction: whereas the overlap of light coming from fibers of the same bundle must be accounted for, one can safely ignore interbundle contamination.
The spectral processing consists in a rowbyrow extraction from each CCD, where the signal of a row of pixels illuminated by a fiber bundle is fit independently from the rest of the data (wavelength are dispersed approximately along CCD columns). The extraction consists in minimizing the following quantity
where D are the pixel values (after preprocessing, i.e. after the subtraction of the bias, conversion from digital units to electrons, and flatfielding), w the weights assigned to each pixel, P the true crossdispersion PSF profile of a fiber and the profile used (a Gaussian). The quantity is the number of collected electrons from a fiber, to be estimated. The last term is a loworder polynomial contribution used to fit any residual light in the row; it compensates both for scattered light and the imperfect modeling of the crossdispersion profile tails.
For the extraction to be optimal, one must have , where var_{i} is the variance of D_{i}. This pixel variance is, however, not known a priori. In the BOSS DR12 pipeline, the pixel data were directly used as an estimate of the expected signal: var_{i} = var_{read} + D_{i} if D_{i}> 0 and var = var_{read} if not, where var_{read} is the readout noise variance, and D_{i} is here an estimate of the Poisson noise.
This use of the data to determine the variance results in a biased estimate of the flux. As an illustration, if one considers a single fiber extraction, without the background polynomial contribution, the χ^{2} minimization yields
To first order, assuming D_{i} = NP_{i} + δ_{i}> 0, and noting the expectation value of w_{i}, one has and (A.1)The first term produces an unbiased estimate of N if , which is not strictly the case. This bias is in principle corrected for by the calibration, but since the weights depend on the amount of signal, the calibration obtained with bright spectra does not perfectly correct the bias of faint spectra. The second term is a systematic underestimation of the signal; it can be interpreted as an effective number of pixels contributing to the measurement (if the profile is normalized, with ∑ _{i}P_{i} = 1).
The fractional bias is more important at low fluxes, and as a consequence i) this bias has a non trivial wavelength dependence correlated with bright sky emission lines and ii) its amplitude relative to the calibrated flux depends on the throughput and hence the observation conditions. Those two effects could lead to correlated biases in the Lyα forest.
We have rerun the data processing using w_{i} = 1/var_{read} (this weight is in practice not a constant quantity per CCD amplifier because of the flatfield correction). We also improved the outlier rejection procedure that discards unmasked pixels affected by cosmic rays.
This new approach to extraction is unbiased but less optimal. Expressing this procedure using a matrix representation, if we define X as the parameters of the model (fluxes and polynomial corrections) such that we can interpret the data as D = HX + δ (where δ is the noise), the optimal extraction results in a covariance of the parameters , whereas the suboptimal extraction produces a larger covariance . We use this latter expression in the subsequent steps of the pipeline, using the estimated flux (and not directly the pixel values) to evaluate the Poisson noise contribution to var_{i}.
The ratio of uncertainties does not exceed 5% over the dynamic range of BOSS observations of quasars. It reaches a higher value of 10% for the brightest measurements (see Fig. A.1).
Subsequent steps of the data processing (fiber flatfield, sky spectrum model, photometric calibration) also involve the weighted average of the data. Similarly to the extraction, using directly the flux inverse variances as weights results in a bias because of the fluxinversevariance correlation. We hence took care of replacing the inverse variances of spectra by their average over many wavelength (using a spline fit) to mitigate this bias and to simultaneously minimize the loss of optimality. Those changes have been incorporated in the pipeline used for the SDSS DR13 (Albareti et al. 2016).
Appendix B: The bias parameters of the Lyα forest
The physical model of the Lyα forest autocorrelation function involves two independent parameters for the linear biasing: the density bias factor b_{Lyα}, and the redshift distortion factor β_{Lyα}. The values for our best fit with errors are given in Table 5. As usual, we give the value of b_{Lyα}(1 + β_{Lyα}) because its error is smaller than the error for b_{Lyα} and less correlated with β_{Lyα} (see Slosar et al. 2011). Some variations of these parameters occur as the model is improved to include the effects of metal lines, and possible effects of HCDs and ionizing intensity fluctuations. As expected, there is a degeneracy in the values of various parameters of the model that affect broadband terms. In particular, the value of β_{Lyα} shows a substantial change from 1.2 to 1.7 depending on the model that is used. The first model, assuming that only Lyα forest correlation is present with a constant ionizing intensity, does not yield a good fit. Taking into account the effect of metal lines is the main improvement of the model that allows for an acceptable fit, and does not change appreciably the value of β_{Lyα} although errors are substantially increased. The addition of the effects of HCDs and UV fluctuations further improves the model and increases substantially the best fit value of β_{Lyα}.
To compare our values of b and β_{Lyα} with the previous determination from the DR11 data from Blomqvist et al. (2015), we first translate their result to the cosmological model we use in this paper, taking into account that the measurement of the Lyα autocorrelation depends only on , where is the mean redshift of our Lyα forest correlation measurement. Blomqvist et al. (2015) measured β_{Lyα} = 1.39 ± 0.11, and b_{Lyα}(1 + β_{Lyα}) = −0.374 ± 0.007, using a model with Ω_{m} = 0.27 and σ_{8} = 0.7877, implying . The model used in this paper has Ω_{m} = 0.3147 and σ_{8} = 0.8298, implying . Therefore the Lyα bias of Blomqvist et al. (2015) translated to our model is b_{Lyα}(1 + β_{Lyα}) = −0.367 ± 0.007. This should be compared with our result on the first row of Table 5, because Blomqvist et al. (2015) did not consider the effects from metals, HCDs or UV fluctuations.
The main reason for the difference with the values of b_{Lyα}(1 + β_{Lyα}) and β_{Lyα} we obtain here for our Lyα forest model applied to DR12 given in Table 5 is the different radial range for the fit: Blomqvist et al. (2015) restricted their fit to r > 40 h^{1} Mpc, whereas we use r > 10 h^{1} Mpc. There are other details that contribute to the difference, e.g., the continuum distortion and bin smoothing modeling, but the radial range of the fit is most important. As mentioned earlier, the χ^{2} value shows that our Lyα fit is not good, reflecting the fact that the shape of the correlation cannot be properly fit over our broad radial range. This is why the values of ∥ b_{Lyα} ∥ (1 + β_{Lyα}) and β_{Lyα} decrease when extending the fit to smaller separations. The better fits obtained for the more complete models in Table 5 show that b_{Lyα}(1 + β_{Lyα}) is relatively stable, but the value of β_{Lyα} can be more model dependent.
These results can also be compared to theoretical predictions from hydrodynamic simulations of the intergalactic medium, which were compared to the observational measurement of Blomqvist et al. (2015) in Sect. 7 of ArinyoiPrats et al. (2015). There, it was noticed that while the predicted value of β_{Lyα} is roughly in agreement with the observational determination, the value of ∥ b_{Lyα} ∥ (1 + β_{Lyα}) was predicted to be lower than observed. Now, the lower value of ∥ b_{Lyα} ∥ (1 + β_{Lyα}) we measure from DR12, correcting for the effect of metals and extending the fit to smaller separations, is in much better agreement with the range of the expected theoretical values. Despite the considerable uncertainties that are still present both in the predictions from hydrodynamic simulations of the Lyα forest (as discussed in ArinyoiPrats et al. 2015), and in the observational determinations that are affected by the modeling of metals, HCDs and ionizing intensity fluctuations, we think that the rough agreement we now find of two independently predicted linear bias factors of the Lyα forest is fairly remarkable. In the future, we expect to be able to better measure these linear factors and their redshift evolution.
All Tables
Parameters of the flat ΛCDM cosmological model used for the production and analysis of the mock spectra and the CMBbased flat ΛCDM model from Planck Collaboration XIII (2016) used for the analysis of the data.
Parameters of the physical model for the correlation function and their bestfit values for the data fit over the range 10 < r < 180 h^{1} Mpc.
Weighted mean of fit parameters and mean uncertainties for six sets of 100 mocks of increasing realism: Lyα only; including a quasar continuum; including metal absorption (for Met1 or Met2); including HCDs (and Met1) and masking those with N_{HI}> 10^{20} cm^{2} or > 10^{21} cm^{2}.
Results of fits of the data assuming increasingly complicated forms for ξ^{ph}(r_{⊥},r_{∥}): Lyα absorption only; including metal absorption; including unidentified HCDs; and including fluctuations in the ionizing UV flux (i.e. the complete physical model whose complete set of parameters is given in Table 3).
Bestfit values of (α_{⊥},α_{∥}) for fits including metals, HCDs and UV fluctuations, with and without powerlaw broadbands (BBs) of the form (i_{min},i_{max},j_{max}) = (0,2,6).
All Figures
Fig. 1
SDSS DR12 footprint (in J2000 equatorial coordinates) used in this work. The survey covers one quarter of the sky (10^{4}deg^{2}). The light blue regions are those added beyond the area covered by our previous study (Delubac et al. 2015). The dotted line is the Galactic plane. 

In the text 
Fig. 2
Example of a BOSS quasar spectrum of redshift 2.91 (smoothed to the width of analysis pixels). The red and blue lines cover the forest region used in our analysis, 104.0 <λ_{RF}< 120.0 nm. This region is sandwiched between the quasar’s Lyβ and Lyα emission lines, respectively at 400.9 and 475.4 nm (restframe 102.572 and 121.567 nm). The blue line is the model of the continuum, C_{q}(λ); the red line is the product of the continuum and the mean transmission, , as calculated by the method described in Sect. 4. 

In the text 
Fig. 3
Weighted redshift distribution of pairs of Lyα forest pixels. The mean is ⟨ z ⟩ = 2.33. Included in the distribution are the ~5 × 10^{10} pairs within 20 h^{1} Mpc of the center of the BAO peak. 

In the text 
Fig. 4
Mean ratio, R(λ), of observed flux to pipelinemodel flux as a function of observed wavelength for quasar spectra to the red of the Lyα emission line (λ_{RF}> 130 nm). (The mean is calculated by weighting each measurement by the inverse of the pipeline variance.) In this mostly unabsorbed region of quasar spectra, the percentlevel wavelengthdependent deviations from unity are due to imperfect modeling of calibration stars and to the calcium H and K lines (393.4 and 396.9 nm) due to Galactic absorption. 

In the text 
Fig. 5
Onedimensional fluxcorrelation function, ξ_{1d}, for BOSS quasars showing correlations of δ_{q}(λ) within the same forest. The correlation function is shown as a function of wavelength ratio for the data and for the mocks (procedure Met1). Prominent peaks due to Lyαmetal and metalmetal correlations are indicated. The peak at λ_{1}/λ_{2} ~ 1.051, which is seen in the data but not in the mocks, is due to CII(133.5)SiIV(140.277) at z ~ 1.85, outside the redshift range covered by the mocks. 

In the text 
Fig. 6
Mean transmission, , calculated using the DR12 pipeline of our previous investigations, and using the DR13 pipeline applied to DR12 data of this analysis. The DR13 pipeline allows us to divide forest spectra by the correction factor (Fig. 4) derived from the spectra on the red side of Lyα emission. This procedure eliminates most of the small scale structure found with the old pipeline that is due to artifacts of stellar modeling, Galactic calcium H and K absorption, and sky lines. 

In the text 
Fig. 7
Correlation, Corr, plotted vs. for the two lowest intervals of . The correlation is averaged over . The left panel is the Metal (Met1) mocks and the right panel is for the data. As labeled, the correlation is calculated by subsampling (Eq. (13)), by assuming “independent forests” (Eq. (14), offset by 0.5 h^{1} Mpc for clarity), or by the mocktomock variation (Eq. (15)). Good agreement between the methods is seen, though the independentforest method necessarily underestimates the correlation for Δr_{⊥} ≠ 0. The differences between the mocks and the data reflect the differences in ξ_{1d} (Fig. 5). 

In the text 
Fig. 8
Correlation between pixels in the Lyαforest and pixels in the CIVforest, ξ^{pl}(λ_{2}−λ_{1},θ) for small wavelength differences. The wavelength difference and angular separation have been transformed to pseudoseparations assuming Lyα absorption in both forests. The red points show the correlation for pairs on the same halfplate . The much smaller correlations for different halfplates are the blue points. The red line represents the prediction for our model of calibration and sky noise. 

In the text 
Fig. 9
Null test for pipelineinduced correlations. As in Fig. 8, the correlation between the Lyαforest pixels and the CIVforest pixels, is shown, but now for four angular ranges as labeled and offset by 0.1 successively for clarity. Pixel pairs with are excluded. There is no evidence for pipelineinduced systematics with the χ^{2} for vanishing correlation for is (32.5, 31.2, 22.7, 40.0) for the four ranges, each with 38 data points. 

In the text 
Fig. 10
Twodimensional representation of r^{2}ξ(r_{⊥},r_{∥}) in units of (h^{1} Mpc)^{2}. The right panel shows the measurement and the left panel the bestfit model, ξ^{ph}, modified by the distortion matrix via Eq. (18). The BAO feature is at r ~ 100 h^{1} Mpc. The effects of metalLyα correlations are seen in the lowest r_{⊥} bin, in particular the peak at 50 <r_{∥}< 70 h^{1} Mpc due to SiIIa and SiIIb. 

In the text 
Fig. 11
Correlation function for the metalfree mocks in four ranges of μ. The black points and curves correspond to mocks with Lyα absorption but without the addition of a quasar continuum. The red points and curves correspond to mocks with the addition of a continuum. The points correspond to stacks of 100 mocks and the light curves to individual mocks. The heavy curves correspond to the input model of Table 1 (after distortion by the matrix D_{AA′} (Eq. (10)) for the red curve). 

In the text 
Fig. 12
Correlation function for the stack of the 100 mocks in four ranges of μ. The red points represent the measured correlation function of the mocks with metals and the red curves shows the best fit. The black points indicate the correlation function after including, in addition to metals, HCDs unmasked for N_{HI}< 10^{20} cm^{2} and the black curve the best fit. The peak at r ~ 60 h^{1} Mpc due to SiIIa and SiIIb is apparent in the range μ> 0.95 but not in the range 0.8 <μ< 0.95. 

In the text 
Fig. 13
Measured α_{⊥} and α_{∥} (left) and β_{Lyα} and b_{Lyα}(1 + β_{Lyα}) (right) for the 100 mock catalogs including HCDs with N_{HI}> 10^{20} cm^{2} masked. There are four outliers on the left plot and two on the right. The horizontal and vertical blue lines show the weighted means of the distributions, also given in Table 4. The distribution and mean values of (α_{⊥},α_{∥}) indicates no significant bias in the reconstruction of the BAO peakposition parameters. 

In the text 
Fig. 14
Measured correlation function in four ranges of μ. The most radial range (topleft μ> 0.95) has, in addition to the BAO peak at ~100 h^{1} Mpc, a peak at ~60 h^{1} Mpc due to correlated absorption by Lyα and SiIII(119.0,119.3) at the same physical position. In the next radial bin (topright, 0.8 <μ< 0.95) only the BAO peak is visible. The lines show fits including successively Lyα and metal absorption, unidentified HCDs absorption, UV flux fluctuations, and a (i_{min},i_{max},j_{max}) = (0,2,6) broadband (BB). 

In the text 
Fig. 15
Contours for (α_{⊥},α_{∥}) at the (68.3, 95.45, 99.73)% CL (solid, dashed, dotted). The red contours are for the physical model with HCDs and UV fluctuations from Tables 3 and 5. The blue (green) contours are for the (i_{min},i_{max},j_{max}) = (0,2,6) powerlaw broadband with (without) priors on the parameters of the physical correlation function. 

In the text 
Fig. 16
Derivative of the normalized Lyαforest fluxdecrement, δ_{LyaF}, with respect to the Milky Way dust column density normalized to its dispersion. The dashed curves represent the 0, 1, 2 and 3σ uncertainties. The red vertical lines mark the positions of the calcium H and K lines. 

In the text 
Fig. 17
Comoving expansion rate, r_{d}H(z)/(1 + z), measured with BAO, with r_{d} in units of 147.3 Mpc. The red square is the present measurement at z = 2.33. Also shown are measurements using the Lyαquasar cross correlation at z = 2.36 (FontRibera et al. 2014), and galaxy correlations at 0.35 <z< 0.65 (Alam et al. 2016), at z = 0.15 (Ross et al. 2015), and at z = 0.11 (Beutler et al. 2011). The points at z = 0.11 and z = 0.15 use the SNIa measurement of q_{0} (Betoule et al. 2014) to convert the measured to D_{H}(z). The black line is the prediction of the Planck flat ΛCDM model (Table 1). 

In the text 
Fig. 18
One and twostandard deviation constraints on the oΛCDM parameters (Ω_{M},Ω_{Λ}) using BAO measurements of D_{M}(z) /r_{d} and D_{H}(z) /r_{d} without CMB measurements, i.e. without imposing the CMB value of r_{d}. The red contours use the work presented here (line 1 of Table 6), the z< 0.8 galaxy data (Alam et al. 2016; Beutler et al. 2011; Ross et al. 2015) and the quasarforest crosscorrelation (FontRibera et al. 2014). The gray contours exclude the quasarforest crosscorrelation. The blue contours use only the z< 0.8 galaxy data and the green contours show the constraints derived from the SNIa Hubble diagram (Betoule et al. 2014). The black dot is the position of the flatΛCDM model that describes CMB data (Planck Collaboration XIII 2016). 

In the text 
Fig. A.1
Comparison of flux measurement uncertainties between the DR12 data reduction (in blue) and this new extraction (in red). This analysis is based on the data from camera b1, plate 7339, where many observations of the same targets were performed. The top panel shows the empirical rms (circles) among the various observations, and the corresponding predicted uncertainty (lines) of the number of electrons per CCD row (N_{e}) as a function of N_{e}. The second panel is the same quantity normalized by the predicted uncertainty of DR12. The bottom panel displays the fraction of 3σ outliers that were discarded in the two other panels. As expected, there is a larger dispersion with the new, less optimal reduction at high flux. The larger dispersion at low flux is not easy to understand; it is probably a consequence of the biased estimator in DR12, where negative statistical fluctuations are given a high weight. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.