Baryon acoustic oscillations at z = 2.34 from the correlations of Ly$\alpha$ absorption in eBOSS DR14

We measure the imprint of primordial baryon acoustic oscillations (BAO) in the correlation function of Ly$\alpha$ absorption in quasar spectra from the Baryon Oscillation Spectroscopic Survey (BOSS) and the extended BOSS (eBOSS) in Data Release 14 (DR14) of the Sloan Digital Sky Survey (SDSS)-IV. In addition to 179,965 spectra with absorption in the Lyman-$\alpha$ (Ly$\alpha$) region, we use, for the first time, Ly$\alpha$ absorption in the Lyman-$\beta$ region of 56,154 spectra. We measure the Hubble distance, $D_H$, and the comoving angular diameter distance, $D_M$, relative to the sound horizon at the drag epoch $r_d$ at an effective redshift $z=2.34$. Using a physical model of the correlation function outside the BAO peak, we find $D_H(2.34)/r_d=8.86\pm 0.29$ and $D_M(2.34)/r_d=37.41\pm 1.86$, within 1$\sigma$ from the flat-$\Lambda$CDM model consistent with CMB anisotropy measurements. With the addition of polynomial"broadband"terms, the results remain within one standard deviation of the CMB-inspired model. Combined with the quasar-Ly$\alpha$ cross-correlation measurement presented in a companion paper Blomqvist19, the BAO measurements at $z=2.35$ are within 1.7$\sigma$ of the predictions of this model.


Introduction
Since the first observations of the imprint of primordial baryonic acoustic oscillations (BAO) as a peak in the galaxy correlation function (Eisenstein et al. 2005) or as a periodic modulation in the corresponding power spectrum (Cole et al. 2005), the BAO signal has led to significant constraints on cosmological parameters.The BAO peak in the radial direction at a redshift z yields D H (z)/r d = c/(r d H(z)), where H(z) is the Hubble parameter and r d is the sound horizon at the drag epoch (Eisenstein & Hu 1998).The transverse measurement constrains the quantity D M (z)/r d = (1 + z)D A (z)/r d , where D A (z) is the angular diameter distance.Because of its sensitivity to both the distance and the expansion rate, the ensemble of BAO measurements yield tight constraints on ΛCDM parameters (Aubourg et al. 2015) even without the use of Cosmic Microwave Background (CMB) data.
Contact e-mail: victoria.de.sainte.agathe@lpnhe.in2p3.frMost BAO measurements have employed discrete objects like galaxies (Percival et al. 2010;Reid et al. 2010;Beutler et al. 2011;Blake et al. 2011;Anderson et al. 2012Anderson et al. , 2014a,b;,b;Ross et al. 2015;Alam et al. 2017;Bautista et al. 2018) or quasars (Ata et al. 2018;Gil-Marín et al. 2018;Hou et al. 2018;Zarrouk et al. 2018).An alternative tracer of the density is the intergalactic medium (IGM), itself traced by Lyα absorption in quasar spectra.Such measurements at z ≈ 2.4 were suggested by McDonald (2003) and McDonald & Eisenstein (2007).The first detections of a BAO peak in the Lyα auto-correlation function (Busca et al. 2013;Slosar et al. 2013) used data from the Baryon Oscillation Spectroscopic Survey (BOSS) in the Sloan Digital Sky Survey (SDSS) data-release 9 (DR9), at an effective redshift of z = 2.3.Delubac et al. (2015), using BOSS in SDSS-DR11, confirmed the detection of a BAO acoustic peak in the Lyα autocorrelation function at the 5σ level.Most recently, Bautista et al. (2017) (B17 hereafter) used Lyα forests from BOSS DR12 data and provided a measurement of D H /r d at 3.4% precision level (or of the optimal combination D 0.7 H D 0.3 M /r d at the 2.5% level).2014) and the auto-correlation of Delubac et al. (2015).
In the present paper, we use quasar spectra from the BOSS survey and from its extended version eBOSS in the SDSS DR14 to study BAO in the Lyα auto-correlation function.The quasar-Lyα cross-correlation is studied in a companion paper (Blomqvist et al. 2019) As in previous measurements, we use Lyα absorption in the "Lyα region" of quasar spectra, i.e., quasar rest-frame wavelengths in the range 104 < λ RF < 120 nm.We call the auto-correlation function using only this region the Lyα(Lyα) × Lyα(Lyα) correlation1 .To increase the statistical power we also include Lyα absorption in the Lyβ regions of quasars, 97.4 < λ RF < 102 nm, correlated with the Lyα absorption in Lyα regions, i.e., the Lyα(Lyα)xLyα(Lyβ) correlation function.The Lyβ region was previously used by Iršič et al. (2013) to investigate the flux transmission power spectrum within individual spectra.
Besides the use of the Lyβ region, the analysis presented here differs in several ways from that of Bautista et al. (2017) based on DR12 data.For each quasar, we now use all observations instead of just the best one.We analyze ∼ 15% more Lyα regions.We have refined the modelling of the weights ( §3.2), the way we take into account the effect of unmasked High Column Density (HCD) systems, and the modeling of nonlinearities in the power spectrum ( §4.1).We have not developed new mock spectra beyond those used in the DR12 analysis though this is being done for the final eBOSS analysis.
Fig. 2. The Lyα and Lyβ spectral regions defined in Table 1.
The layout of the paper is the following.In §2, we present the Lyα and Lyβ spectral region samples used in the present study.We compute correlation function of Lyα absorption for the DR14 data in §3 and present our physical model for this function in §4.The results of fitting the data are presented and discussed in §5.We draw cosmological conclusions in §6 and summarize our results in §7.
The computations of the correlation functions presented in this paper have been performed using a dedicated software package, Package for Igm Cosmological-Correlations Analyses (picca), developed by our team2 .

Data sample and reduction
The extended Baryon Oscillation Spectroscopic Survey (eBOSS; Dawson et al. 2016) is the extension of the BOSS experiment (Dawson et al. 2013) which aims to measure cosmology with BAO using optical spectra from quasars, emission line galaxies and luminous red galaxies.It is one of the four projects of the fourth stage of the Sloan Digital Sky Survey (SDSS-IV; Blanton et al. 2017).
The optical spectra are collected from 1000 fibers, attached to the focal plane of a 2.5 m telescope in Apache Point Observatory (Gunn et al. 2006), by two spectrographs in the wavelength range [360,1000] nm (Smee et al. 2013).The spectral resolution of the spectrographs is ≈ 2000.
In this paper, we use the forests of the high-redshift quasar sample from the SDSS Data Release 14 (DR14; Abolfathi et al. 2017).This sample contains the first two years of eBOSS data and the five years of BOSS observations reprocessed using the eBOSS pipeline.It also includes data from the ancillary programs Time-Domain Spectroscopic Survey (TDSS) and SPectroscopic IDentification of ERosita Sources (SPIDERS).The quasar target selection is presented in Myers et al. (2015).Note that eBOSS also targets quasars at low redshifts (where the Lyα region is not observable) to be used in other programs (Ata et al. 2018;Wang et al. 2018;Blomqvist et al. 2018;Zhao et al. 2019;du Mas des Bourboux et al. 2019).
The automated data reduction is organized in two steps (Dawson et al. 2016).The pipeline initially extracts the twodimensional raw data into a one-dimensional flux-calibrated spectrum.During this procedure, the spectra are wavelength and flux calibrated and the individual exposures of one object are coadded into a rebinned spectrum with ∆ log(λ) = 10 −4 .The spectra are then classified as STAR or GALAXY or QSO, and their redshift is estimated.Objects that cannot be automatically classified are visually inspected (Pâris et al. 2017) and a quasar catalog is produced, which contains 526,356 quasar spectra with redshift 0 < z < 7.Among these objects, 144,046 were not in DR12.The coverage footprint of DR14 quasars is presented in Fig. 1 In this study, we examine both Lyα and Lyβ regions (see Fig. 2).The Lyα region in the quasar spectrum lies between the Lyα and the Lyβ emission peaks.We limit its coverage to the rest-frame wavelength range [104,120] nm in order not to Table 1.Definitions of the Lyα and Lyβ regions in terms of restframe wavelength range.Also listed are the allowed observer frame wavelength ranges, the corresponding quasar redshift ranges and the number of forests available in our sample.
In the DR14 quasar catalog, selecting quasar redshifts in the range [2.0, 3.5] yields to 216,162 spectra containing, at least partially, the Lyα region, and selecting quasar redshifts in the range z q ∈ [2.53, 3.5] yields 86,245 spectra containing the Lyβ region.We choose z q = 3.5 as an upper limit, as beyond this redshift the quasar density is insufficient to measure correlations and the rate of redshift misidentification is large (Busca & Balland 2018).The requirement that the observed wavelength must be greater than 360 nm is due to the low CCD response and atmospheric transmission in the UV region.
In order to mask damped Lyα systems (DLA), we use the updated DR14 DLA catalog of Noterdaeme et al. (2009Noterdaeme et al. ( , 2012)), which contains 34,541 DLA in 27,212 forests.The absorption of the identified DLAs are modeled with a Voigt profile and the regions with more than 20% of absorbed flux are masked.For the Lyβ regions, we apply this procedure both for Lyα and Lyβ strong absorbers.We also mask the sky emission and absorption lines listed on the SDSS website 3 .The Broad Absorption Line (BAL) quasars are automatically identified (Pâris et al. 2017) and excluded from the data, leaving a sample of 201,286 objects for the Lyα regions and 80,443 for the Lyβ regions.
For the determination of the correlation function, we divide spectra into "analysis pixels" that are the inverse-varianceweighted flux average over three adjacent pipeline pixels.Throughout the rest of this paper, "pixel" refers to analysis pixels unless otherwise stated.Spectral regions with less than 50 such pixels in regions or which have failed the continuum-fitting procedure (Sec.3.1) are discarded.These selection criteria produce 179,965 Lyα regions (compared to 157,783 in B17) and 56,154 Lyβ regions (see Table 1).
The analysis procedure described in the next section assigns redshifts to the observed pixel wavelengths by assuming that flux decrements, in both the Lyα and Lyβ regions, are due to Lyα absorption.The effect of non-Lyα absorption is taken into account in the correlation-function model presented in §4.The weighted distribution of the redshifts of pairs of pixels used to measure Lyα(Lyα) × Lyα(Lyα) and Lyα(Lyα) × Lyα(Lyβ) correlations are presented in Fig. 3.The mean redshift of the combined set of pixel pairs is 2.34.3).The computation of the correlation function requires an estimation of the transmission field along the line-of-sight (LOS) towards surveyed quasars.This field arises due to the presence along the LOS of intergalactic gas.More precisely, for the correlation calculation, we only need to know the flux fluctuations around the average transmitted flux spectrum in the forests of quasar q at wavelength λ.We thus define the field δ q (λ), for each quasar q under investigation, as:

Computing the Lyα correlation function from the data
where f q (λ) denotes the observed flux of quasar q at observed wavelength λ, C q (λ) is the continuum flux and F(z) the mean transmission at the absorber redshift z.
We estimate the quantity C q F(z) from the average of the transmitted flux of all forest spectra in the sample: q w q (λ RF ) f q (λ RF ) q w q (λ RF ) . ( where λ RF is the rest-frame wavelength and w q is a weight (see § 3.2).For each quasar, f (λ RF ) is then multiplied by a linear polynomial function of Λ ≡ log(λ) to account for the diversity of quasar luminosity and spectral shape: There are thus two adjustable parameters per quasar, a q and b q .
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).
The fitting procedure to determine (a q , b q ) forces to zero the mean and spectral slope of δ q (λ) for each quasar, thus introducing spurious correlations in the measured field.To make it easier to deal with this distortion in the analysis, we follow B17 by transforming the measured δ q (λ) to δq : where where δ K i j denotes the Kronecker symbol.The advantage of this transformation is that it makes the distortion of the true field introduced by the continuum fit procedure explicit, and, as a consequence, simplifies the link between the true correlation function and the measured, distorted one (see § 4.3).
The statistics of δq (λ) within individual forests are described (in part) by the so-called one-dimensional correlation function,  3 lists the important observed transition pairs.(See also Pieri et al. (2014)).

Pixel weights
The pixel weights are proportional to the inverse of the variance of δ q (λ).Following Blomqvist et al. (2018), the variance is modeled as the sum of three terms: The noise pixel variance is σ 2 noise = σ 2 pip /(C q F) 2 where σ 2 pip is the pipeline estimate of the pixel variance.The intrinsic, redshift dependent, contribution of the density fluctuations underlying Lyα regions is σ 2 LSS .The third term, (λ)/σ 2 noise , takes into and distance separation r i j .The radial separation r ,i j is the projection of r i j on the median LOS and the transverse separation r ⊥,i j is the LOS perpendicular component of r i j , assuming the flat Pl2015 model (Table 2).
account differences between the fitted quasar spectrum and the individual spectrum of quasar q (these differences appear at high signal-to-noise).The functions η(λ) and (λ) correct for imperfections of the pipeline estimates and differences between the average and individual spectra, respectively.Following Busca et al. (2013), the weights are corrected to take into account the expected redshift dependence of the correlation function amplitude: where the Lyα bias redshift-evolution parameter, γ α = 2.9 (Mc-Donald et al. 2006) and λ α is the Lyα restframe wavelength.
In practice, one starts with an initial estimate of the weights, allowing a first estimate of the mean spectrum f (λ RF ) (eqn. 2) and the quasar parameters a q and b q (eqn.3).The functions η(λ), (λ) and σ LSS (λ) are then fit and the mean spectrum is then recalculated with the new weights.This process is repeated until stable values are obtained after about five iterations.

The correlation function
To compute the correlation function, we correlate absorption at an observed wavelength λ i in the LOS of a given quasar q, with absorption at an observed wavelength λ j in the LOS of another quasar q .Assuming the absorption is due to the Lyα transition, one can compute, from the values of λ i and λ j , the redshifts z i and z j of the matter absorbing these lines.Each pair of absorbers (z, q) entering the computation defines a "pixel" in real space and we call r i j the physical separation between two such pixels i and j (see Fig. 5).This distance is calculated assuming the Pl2015 cosmology (Table 2).The distance r i j can be projected on the radial and the transverse directions, leading to two components r ,i j and r ⊥,i j .These components can be expressed in terms of the comoving distances D(z i ) and D(z j ) from us to absorbers i and j and the subtended angle between the two LOS, θ i j , as: We then define bins of (r ,i j , r ⊥,i j ) on a 2D grid.In practice, the grid uses 2500 bins of dimensions 4h −1 Mpc × 4h −1 Mpc over 0 < r ⊥ < 200h −1 Mpc and 0 < r < 200h −1 Mpc.For a given bin in this grid, A, we consider each pair of pixels (i, j) whose r and r ⊥ coordinates fall on this bin.The measured correlation function in bin A reads: with w k ≡ w q k (λ k ) and δk ≡ δq k (λ k ).
We discard from the computation all pixel pairs belonging to the same LOS, since two pixels belonging to the same quasar spectrum are affected in a correlated way by the fitting procedure described in §3.1.Likewise, pixels belonging to the same half plate at the same wavelength are excluded, to avoid unphysical correlations induced by the extraction pipeline.

The covariance matrix
The covariance between two bins A and B is defined as: where .... denotes an ensemble average.Following Delubac et al. (2015) and B17, we estimate equation ( 10) by dividing the eBOSS footprint in N h = 876 sky pixels, using the HEALPix tessellation scheme (see Górski et al. 2005), and by equating the ensemble averages of equation ( 10) with the weighted mean over these sky pixels: and with W h A the sum of the weights of pairs in sky pixels h contributing to bin A. Similarly, ξ h A is the correlation function of pairs in sky pixels h that contribute to bin A.
In practice, for the computation of the correlation function, a pair (i, j) is attributed to the sky pixel of the first quasar of the pair, and the pair ( j, i) is never considered, insuring that a pair is not counted twice in the calculation.
In this approximation, we assume that each sky pixel provides an independent realization of the δ field.This statement is not exactly true as correlations do exist between pairs in different sky pixels, but these correlations are small (e.g., Delubac et al. 2015).
We thus compute the covariance matrix defined in equation (10) using the following expression: where ξA is given by ( 9).Due to the finite number of sky pixels, the estimate ( 13) is noisy and must be smoothed before it can be used in fits.We perform the smoothing by approximating the correlation,   As a check of the subsampling method, the covariance can also be estimated by neglecting inter-forest correlations, in which case the four-point function vanishes unless the four pixels are drawn from just two spectra: where ξ 1d is the intra-forest correlation function shown in Fig. 4.
The sum can then be estimated from a random sample of forest pairs.Because neighboring forests are nearly parallel, the sum necessarily produces C AB = 0 unless r A ⊥ ∼ r B ⊥ .Because the Lyα and Lyβ forests have different ξ 1d , we expect differences between the covariances for Lyα(Lyα) × Lyα(Lyα) and Lyα(Lyα) × Lyα(Lyβ) correlations.These differences are illustrated in Fig. 6 showing, for the two lowest values of ∆r ⊥ , the correlation, Corr AB .For ∆r ⊥ = 0, there is good agreement between the subsampling (13) and independent-forest (14) calculations.

Modeling the correlation function
This section details the model of the Lyα auto-correlation function that will be fitted against the estimator ξA of equation ( 9).Table 4 lists the different parameters of the model.The model includes two components: one is the Lyα-only correlation function computed from Lyα absorption only; the other component incorporates the contribution to the correlation function of absorption by metals, for which the nominal separation for the pixel is not the true separation (see §4.2).We thus write: The following subsections describe these two components.Section 4.3 explains how the model is "distorted" to fit the data.

Lyα−Lyα mod
We start from the CAMB linear power spectrum (Lewis et al. 2000) which is decomposed into a smooth and a peak components, following the side-band technique described by Kirkby et al. (2013).This allows one to constrain the position of the BAO peak independently of the correlation function at scales much smaller or much larger than the BAO distance scale.We thus model the matter power spectrum as the sum of two terms corresponding to the smooth and the peak terms in the correlation function.Moreover, in order to incorporate the effects of the non-linear growth of matter that lead to broadening of the BAO peak, the peak term is corrected by a Gaussian factor (Eisenstein et al. 2007).The "quasi-linear" power spectrum hence reads: where P smooth (k, z) and P peak (k, z) are the power spectra of the smooth and peak components, and Σ and Σ ⊥ represent the RMS displacements in the parallel and transverse directions, respectively.We adopt the values of Kirkby et al. (2013) for these parameters: Σ = 6.41h −1 Mpc and Σ ⊥ = 3.26h −1 Mpc.The power spectrum is obtained from P QL as where d Lyα is the Kaiser factor (Kaiser 1987) for the Lyα absorption and D NL (k) takes into account non-linear effects.The function G(k) models the effect of binning of the correlation function on the separation grid.
The Kaiser factor can be written as: where (b Lyα , β Lyα ) and (b HCD , β HCD ) are the bias parameters associated with the IGM and HCD systems and F HCD is a function defined below.
Following McDonald et al. (2006), we assume that the product of b Lyα and the growth factor of structures varies with redshift as (1 + z) γ α −1 , with γ α = 2.9, while we make use of the approximation that β Lyα does not depend on redshift.
HCD absorbers are expected to trace the underlying density field and their effect on the flux-transmission 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.This broadening effect introduces a k dependence of the effective bias (Font-Ribera & Miralda-Escudé 2012).Following the study of Rogers et al. (2018), we adopt a simple exponential form, F HCD = exp(−L HCD k ), where L HCD is a typical length scale for these systems.DLA identification is possible if their width (wavelength interval for absorption greater than 20% ) is above ∼ 2.0 nm, corresponding to ∼ 14h −1 Mpc in our sample.Based on results from Rogers et al. ( 2018) results, we impose L HCD ∼ 10h −1 Mpc while fitting for the bias parameters b HCD and β HCD .Fixing L HCD is necessary because otherwise the model becomes too unconstrained.We have verified that setting L HCD in the range 7 < L HCD < 13 h −1 Mpc does not significantly change the inferred BAO peak position.
We focus on the minimal model able to reproduce the data, designated as "baseline model".This baseline model does not include the correction of the UV background fluctuations (Pontzen et al. 2014;Gontcho A Gontcho et al. 2014) used in B17.We discuss the improvement of the fit when this UV correction is added, in section 5.
The function D NL (k) accounts for non-linear effects such as thermal broadening, peculiar velocities and non-linear structure growth.A fitting formula for D NL is given by equation ( 21) of McDonald (2003) and has been extensively used in previous studies.More recently, Arinyo-i-Prats et al. (2015) proposed a new fitting formula involving 6 free parameters given by their equation (3.6).Besides reducing number of free parameters with respect to McDonald (2003), it has the correct behavior at small wavenumber k and an explicit dependence on P QL (k), whereas this dependence is only implicit in the McDonald (2003) formula.In practice, the two approaches yield similar results but for the above reasons, we adopt the formula of Arinyo-i-Prats et al. (2015) in the present work and linearly interpolate the parameter values from their Table 7 at the effective redshift z = 2.34.
To account for the effect of the binning of the correlation function on the separation grid, we assume the distribution to be homogeneous on each bin4 and compute the function G(k) as the product of the Fourier transforms of the rectangle functions that model a uniform square bin: where R and R ⊥ are the radial and transverse widths of the bins, respectively.
The two terms in P QL (k, z) (eqn.16) are Fourier transformed to the smooth and peak components of the correlation function: The amplitude of the peak, A peak , is fixed to unity in the standard fit.In the peak component we have introduced the parameters (α , α ⊥ ) which allow us to fit for the peak position independently of the smooth component: where z is the effective redshift of the measurement and the suffix 'fid' denotes the Pl2015 cosmology from Table 2.

The contamination by metals
The second term in the model correlation function (eq.15) accounts for absorption by metals along the quasar LOS.Such absorption is correlated with Lyα absorption (Pieri et al. 2014) and can be used as a tracer of the density field (Blomqvist et al. 2018;du Mas des Bourboux et al. 2019).Here, it is a complicating factor in the analysis because the redshifts of pixels are calculated assuming Lyα absorption.
The important metals can be seen in the 1D correlation function, ξ 1d (λ 1 /λ 2 ), shown in Fig. 4. Column 2 of Table 3 lists the wavelength ratios for the main metal/metal and metal/Lyα absorption correlations, relevant for the Lyα autocorrelation function computation.The corresponding apparent radial separation at vanishing physical separation is where z is the mean redshift of the pair.Values are given in Table 3 for z = 2.34.We model the power spectrum of each pair of absorbers, (m, n), with the same form as that for Lyα -Lyα absorption (17) except that HCD effects are neglected: Since the b m and β m are mostly determined near (r ⊥ , r ) ∼ (0, r ap ), they cannot be determined separately.We therefore fix β CIV(eff) = 0.27 (Blomqvist et al. 2018).For the other metal species we keep the value β = 0.50 used in Bautista et al. (2017) which comes from DLA measurements (Font-Ribera & Miralda-Escudé 2012).
The Fourier transform of P mn (k, z) is then the model correlation function of the pair (m, n): ξ m−n mod (r , r⊥ ), where (r , r⊥ ) are the separations calculated using the correct restframe wavelengths, (λ m , λ n ).
Since we assign a redshift, z α , assuming Lyα absorption, the rest-frame wavelength we ascribe to a metal transition m observed at wavelength λ i is not equal to the true rest-frame wavelength λ i /(1 + z m ), where z m is the true redshift of the metal absorber.This misidentification results in a shift of the model contaminant correlation function.For each pair (m, n) of contaminants, we compute the shifted model correlation function with respect to the unshifted model metal correlation function ξ m−n mod by introducing a metal matrix M AB (Blomqvist et al. 2018), such that: where: and (m, n) ∈ A refers to pixel separation computed assuming z α , and (m, n) ∈ B to pixel separation computed using the redshifts of the m and n absorbers, z m and z n .We take into account the redshift dependence of the weights, equation ( 7), in the computation of w m and w n .The total metal contaminant correlation function, ξ metals mod , is the sum of all the ξ m−n mod contributions, where (m, n) runs over all the involved transition pairs for the Lyα auto-correlation function, see Table 3:

The distorted model
The model correlation function, ξ mod of eqn.( 15), cannot be fit directly to the estimated correlation function ( 9) because the measured δ(λ) are only related to the true δ(λ) through the transformation (4).Following B17, we can account for this effect in the fit by using a distorted model: where D AB is the distortion matrix which, following equation ( 5), is given by The accuracy of this method of accounting for the distortion of the correlation function was tested with mock data sets by Bautista et al. (2017).
In practice, to avoid prohibitive computational time, the distortion matrix is computed using only a random 5% portion of the total number of pairs.

Fitting the BAO peak position
Table 5 presents the best-fit parameters for the Lyα(Lyα) × Lyα(Lyα) correlation function alone and those including the Lyβ region, i.e. the Lyα(Lyα) × Lyα(Lyα+Lyβ) correlation function.Figure 8 displays data for the latter in four ranges of µ along with the best fits.The BAO peak is apparent for µ > 0.8 and is suggested for 0.5 < µ < 0.8.
The BAO parameters for the fit using both Lyα and Lyβ regions are Using the D/r d values for the Pl2015 cosmology in These results can be compared with those of B17, who found α = 1.053 ± 0.36 and α ⊥ = 0.965 ± 0.55 at z = 2.33 using only the Lyα region.These values are very near the present results using only the Lyα region: α = 1.047 ± 0.35 and α ⊥ = 0.960 ± 0.41.Our use of the Lyβ region produces consistent results, given the increase in the data set; the main improvement is that on the precision on D M /r d by ∼ 25%.Constraints on the BAO parameters (α , α ⊥ ) are presented in Fig. 9. Following the method introduced and described in detail in du Mas des Bourboux et al. (2017), we estimate the relation between ∆χ 2 = χ 2 − χ 2 min and confidence levels (CLs) for the BAO parameters using a large number of simulated correlation functions generated from the best-fit model and the covariance matrix measured with the data.The results of the study, summarized in Table D.1, indicate that the (68.27,95.45%)confidence levels for (α , α ⊥ ) correspond to ∆χ 2 = (2.74,7.41) (instead of the nominal values ∆χ 2 = (2.3,6.18)).These levels are shown as the red contours in Fig. 9 for the Lyα(Lyα) × Lyα(Lyα+Lyβ) fit.The best fit is within one standard deviation of the Pl2015 model.
In addition to the baseline fits, we performed a variety of non-standard fits to verify that our BAO results are robust.The results of this exercise are given in Table 6.The first group of fits starts with the simple Kaiser redshift-space model and then includes progressively metals, HCD absorption, and UV background fluctuations.Including HCD absorption is necessary to obtain a good χ 2 but adding UV fluctuations such as those characterized in Gontcho A Gontcho et al. (2014) does not improve the fit, justifying our choice of ignoring the UV issue in the baseline fit.
An important test of systematic effects in the position of the BAO peak is performed by adding polynomial "broadband" terms to the correlation function (before distortion).We follow the procedure and choice of broadband forms used by B17 and adopt the form where the L j are Legendre polynomials.
We want to ensure that the power-law terms model variations of the slowly-varying part of the correlation function under the BAO peak.We therefore perform these fits only over the restricted range 40 < r < 180h −1 Mpc, avoiding introducing undue influence of the 10 < r < 40h −1 Mpc range on the amplitudes of the power laws.Following B17 we fit with (i min , i max ) = (0, 2) corresponding to a parabola in r 2 ξ smooth underneath the BAO peak.We set j max = 6, giving four values of j corresponding to approximately independent broadbands in each of the four angular ranges in Fig. 8.
We performed the broadband fits in two ways.The first placed "physical priors" on (b Lyα , β Lyα , b HCD ) in the form of a Gaussian of mean and width of the fit without broadband terms.Such priors ensured that the broadband terms were relatively small perturbations to the physical model.The second type of fit placed no priors on (b Lyα , β Lyα , b HCD ).
The results of these fits are given in Table 6 and Fig. 9.We see that the addition of such terms does not change significantly The curves show the standard fit and the two fits with broadband terms defined by eqn.32 with (i min , i max , j max ) = (0, 2, 6) with and without additional priors, as described in the text.
α but does shift α ⊥ by 0.5σ or 0.7σ for fits with and without physical priors.This effect was already seen in B17 but, with less significance.Figure 9 shows that in all cases the BAO peak position is within one standard deviation of the prediction of the Pl2015 model.
The fits described above of the Lyα(Lyα) × Lyα(Lyα) and Lyα(Lyα) × Lyα(Lyβ) correlation functions are the primary results of this paper.We also performed fits with two redshift bins, as described in Appendix B. Each of the two redshifts yielded values of (α , α ⊥ ) that are within 1.2σ of the Pl2015 model.(Fig. B.2).We also fit the Lyα(Lyα) × Lyβ(Lyβ) correlation as described in Appendix C. Adding the Lyβ absorption data does not add a significant signal to the BAO peak, but it does allow us to measure the Lyβ bias parameters.
Finally, we combine the measurement of Lyα autocorrelation function of the present analysis with the Lyα -quasar cross-correlation measurement of Blomqvist et al. (2019) by performing a joint fit of the two correlation functions.We use the baseline models of the two analyses and consider the errors to be independent.The joint fit has 18 free parameters and the effective redshift is z = 2.34.The results are given in the column four of Table 7 and the constraints on (α ⊥ , α ) in the right panel of Fig. 9. From this combined fit, we obtain: corresponding to: The value of χ 2 for (α = 1, α ⊥ = 1) is 4.99 greater than the best fit.Using the confidence levels of Table D.1, we conclude that the results of the combined fit are 1.7σ from the predictions of the Pl2015 model (Planck Collaboration et al. 2016).
Fig. 11 plots this value along with other measurements.The data are consistent with the expected behavior of deceleration at high redshift followed by acceleration at low redshift.Independent of CMB data and without assuming flatness, the BAO data by themselves constrain the parameters   2019) is the blue dot.The other points are computed using galaxy measurements (Beutler et al. (2011), Ross et al. (2015), Alam et al. (2017)).The points at z = 0.106 (Beutler et al. 2011) and z = 0.15 (Ross et al. 2015) are converted from D V to H(z) using the SNIa measurement of q 0 given by Betoule et al. (2014).Solid black line shows the Pl2015 values (Planck Collaboration et al. 2016).

Conclusions
We have used Lyα and Lyβ spectral regions from the BOSS and eBOSS DR14 data sample to study BAO.Following B17, we have built a model for the Lyα auto-correlation function that we have then fit to the data.Our model incorporates the effects of redshift space distortions, the non-linear growth of matter, the contamination by metals and the modeling of high column density systems along lines-of-sight to quasars.Including UV fluctuations has only a minor impact on the fit results.We measure the ratios D H /r d and D M /r d at the average redshift of pixel pairs, z = 2.34.We have also performed a measurement of these ratios from the Lyα auto-correlation function in two redshift bins, at z = 2.19 and z = 2.49.
The D H /r d ratio is measured with a precision of ∼ 3.3%, a slight improvement over the precision obtained by B17 for this ratio.The D M /r d ratio is measured with a precision of ∼ 4.4%, which represents an improvement of about 25% with respect to B17.The cosmological measurements obtained in this analysis are in agreement with the predictions of the flat ΛCDM model (Pl2015) favored by the measurement of CMB anisotropies by Planck.
We have also combined the measurements of the present analysis with the ones obtained from the cross-correlation of Lyα absorption and quasars by Blomqvist et al. (2019).The latter alone favors a value of the D H /r d ratio ∼ 3% higher than the one favored by the Lyα auto-correlation.As a result, the best-fit value of D H /r d for the combined fit is shifted towards a higher value than the best-fit from the Lyα auto-correlation alone.Combining the measurement of Lyα auto-correlation (this paper) with the quasar -Lyα cross-correlation of Blomqvist et al. (2019), the BAO measurements at z = 2.34 are within 1.7σ of the predictions of the Pl2015 model.
The ensemble of BAO measurements is in good agreement with the Pl2015 model (Planck Collaboration et al. 2016).They provide an independent way of determining cosmological parameters that is based only on low-redshift measurements.As illustrated in Fig. 12, the BAO results are also consistent with the recent Pantheon SNIa results (Scolnic et al. 2018).
The present measurements will be much improved by the greater statistical power of the upcoming DESI (DESI Collaboration et al. 2016) and WEAVE-QSO (Pieri et al. 2016) projects.The challenge will be to improve the physical modeling of the correlation function in order to fully profit from the improved data.

Appendix A: Effective redshift of the fitted parameters
In this section we present a method to determine the region in (r ⊥ , r , z) space that is most constraining for the various parameters in the fits of the correlation function.We can expect that the parameters (α ⊥ , α ) are mostly determined by (r ⊥ , r ) bins near the BAO peak and at a redshift near the mean redshift of the pixel pairs used in the BAO region.In fact, previous studies (B17, Busca et al. (2013); Delubac et al. (2015)) defined the effective redshift of the BAO measurement in this way.Here, we make this intuitive conclusion more precise by using a Fisher matrix analysis.We use the Fisher matrix formalism as follows: given a parameter p varying linearly with redshift, we define the effective redshift z 0 at which it is measured by: where p 0 is the value given by the fit at z = z 0 .The covariance matrix C p between two parameters p 0 and p 1 is given by: where σ 2 i is the variance of the parameter p i and ρ is the correlation coefficient between p 0 and p 1 .By definition, C p is the inverse of the Fisher matrix F p : with m i the model at bin i.In the case of the linear redshift dependency (A.1), the Fisher matrix F p at redshift z, computed using the set of all the fitted parameter values, {λ 0 }, is given by : with z i the mean redshift of the pairs in bin i.We represent the quantities ∂m i ∂p for 4 of the 12 fitted parameters of the Lyα auto-correlation function in Fig. .1.The covariance matrix C p (z) then reads : Since M i j is symmetric, the determinant of the Fisher matrix, |F p |, does not depend on redshift and is given by: The variance of p 0 at redshift z becomes: The effective redshift z 0 is the value which minimizes the error on p 0 : M i j (z i − z 0 ) = 0, (A.9) i.e., z 0 = i j M i j z i i j M i j . (A.10) In the case of a combined fit, we compute one matrix M d for each correlation function entering the fit:  The graphs show which pixels contribute the most to the constraints on the considered parameter.BAO parameters α and α ⊥ are constrained by the bins around the location of the BAO peak, while the Lyα bias is mostly constrained by the bins at ∼ zero separation.We also note that the Si iii(1207) bias is mostly constrained by the bins r ⊥ ∼ 0, r ∼ 21h −1 Mpc, in agreement with the r ap apparent separation given in Table 3.
to assign forest pairs to the high or low-redshift sample by cutting on the mean of the maximum z of the two forests.
We thus evaluate, for all pairs of forests (i, j), the following quantity: where 1 + z k max,abs = max(λ k obs )/λ abs , with max λ k obs the last pixel of forest k, and λ abs the rest-frame wavelength of the considered transition.The condition z i j < z cut defines the low redshift bin, while the opposite condition defines the high redshift one.
The value of z cut is tuned so that the sum of the weights of all absorber pairs in the high redshift correlation function equals the one of all pairs in the low redshift correlation function.This ensures the two bins have a comparable statistical power.This process leads us to select z cut = 2.5.The redshift distribution of absorber pairs obtained in this way are shown on  Table B.1 shows the associated best-fit parameters.From the table, we note that β Lyα is notably different at low and high redshift.When fitting the full sample, we assumed, following Kirkby et al. (2013), β Lyα to be constant.A redshift dependent β Lyα could thus be an improvement in future analyses.
From the lower right panel of Figure B.3, we see that the amplitude of the high-redshift correlation function is higher than the amplitude of the low-redshift one.This is expected, as b Lyα increases with redshift (Kirkby et al. 2013).

Fig. 1 .
Fig. 1.Sky distribution of the 216,163 quasars with redshift in the [2.0,3.5]range in the DR14 footprint of the BOSS and eBOSS surveys.The high-density regions are the eBOSS and SEQUELS observations (for the highest declinations in the two Galactic caps) and SDSS-stripe 82 (on the celestial equator in the south galactic cap).

Fig. 4 .
Fig. 4. The one-dimensional correlation functions, ξ 1d , in the Lyα (red curve) and Lyβ (blue curve) regions as a function of the ratio of transition wavelengths.Peaks are due to absorption by the two labeled elements at zero physical separation (Table3).
Fig. 4 presents this function for the Lyα and Lyβ forests.The peaks are due to absorption by different transitions at the same physical position.Table

Fig. 5 .
Fig.5.Definition of the coordinates of pixels used in the computation of the correlation function.Absorbers i and j have angular separation θ i j and distance separation r i j .The radial separation r ,i j is the projection of r i j on the median LOS and the transverse separation r ⊥,i j is the LOS perpendicular component of r i j , assuming the flat Pl2015 model (Table2).
as a function of ∆r = |r A − r B | and ∆r ⊥ = |r A ⊥ − r B ⊥ | only, ignoring the small dependence on r and r ⊥ .

Table 4 .
Figure 7 displays the Corr AB between the Lyα(Lyα) × Lyα(Lyα) and Lyα(Lyα) × Lyα(Lyβ) correlation functions.The Corr AB are less than a percent, and will be ignored in the fits of the correlation functions.
where b Lyα is the effective bias of Lyα absorbers with respect to the underlying matter density field, β Lyα is the effective redshift space distortion (RSD) parameter, and µ k = k /k.The two effective parameters (b Lyα and b Lyα β Lyα ) combine Lyα absorption in the IGM and in unmasked high-column density (HCD) systems, i.e., HI absorbers with column densities N HI > 10 17.2 cm −2 :

Fig. 8 .
Fig. 8. Weighted combination between measured Lyα(Lyα)xLyα(Lyα) and Lyα(Lyα)xLyα(Lyβ) correlation functions along with the model best fits in four ranges of µ = r /r.The curves show the standard fit and the two fits with broadband terms defined by eqn.32 with (i min , i max , j max ) = (0, 2, 6) with and without additional priors, as described in the text.

Fig. 9 .
Fig. 9.The left panel shows the 68% and 95% confidence level contours in the (α ,α ⊥ ) plane from the Lyα(Lyα) × Lyα(Lyα+Lyβ) auto-correlation function for the standard fit and for fits with polynomial broadband (BB) terms with and without additional priors, as described in the text.The right panel shows the contours for Lyα(Lyα) × Lyα(Lyα+Lyβ) auto-correlation standard fit and those from the combined fit of the auto-correlation and the quasar-Lyα cross-correlation of Blomqvist et al. (2019).In both panels the value for the Pl2015 model (Planck Collaboration et al. 2016) is shown as a black point.

Fig. 12 .
Fig. 12.One and two standard deviation constraints on (Ω m , Ω Λ ).The red contours use BAO measurements of D M /r d and D H /r d of this work, of Blomqvist et al. (2019) and Alam et al. (2017), and the measurements of D V /r d of Beutler et al. (2011), Ross et al. (2015), Ata et al. (2018) and Bautista et al. (2018).The gray contours do not use the Lyα -quasar cross-correlation measurement of Blomqvist et al. (2019).The green contours show the constraints from SN Ia Pantheon sample (Scolnic et al. 2018).The black point indicates the values for the Planck (2016) best-fit flat ΛCDM cosmology.
model for the correlation function d at bin i.In this case, z 0 reads: the effective redshifts at which the α and α ⊥ parameters are measured for the different correlation functions computed in this paper.The effective redshift values differ by less than 0.5% for α and α ⊥ .Figure .1 shows the quantities ∂m/∂p in the (r ⊥ , r ) plane for the fitted parameters p ∈ [α , α ⊥ , b Lyα , b SiIII(1207) ]. m is the baseline model for the Lyα auto-correlation.

Fig. . 1 .
Fig. .1.The quantity ∂m/∂p in the (r ⊥ , r ) plane for the fitted parameters p ∈ [α , α ⊥ , b Lyα , b SiIII(1207) ]. m is the baseline model of the Lyα autocorrelation function.The graphs show which pixels contribute the most to the constraints on the considered parameter.BAO parameters α and α ⊥ are constrained by the bins around the location of the BAO peak, while the Lyα bias is mostly constrained by the bins at ∼ zero separation.We also note that the Si iii(1207) bias is mostly constrained by the bins r ⊥ ∼ 0, r ∼ 21h −1 Mpc, in agreement with the r ap apparent separation given in Table3.
Figure B.1.The average pair redshift is z = 2.19 and z = 2.49 for the low and high redshift bin, respectively.
Figure B.3 presents the result of fitting the Lyα(Lyα) × Lyα(Lyα) correlation baseline model to the data in the low (blue points) and high (red points) redshift bins in the usual four µ wedges.

Figure B. 2
Figure B.2 presents the constraints obtained, in the (α ⊥ ,α ) parameter space, from fitting the Lyα auto-correlation Lyα(Lyα)xLyα(Lyα) in the low (blue contours) and high (red contours) redshift bins.The values of ∆χ 2 corresponding to a given confidence level were taken from Table D.1.Both the high and low redshift measurements are within 1.2σ of the Pl1015 model.

Fig. B. 1 .
Fig. B.1.Pixel pair redshift distribution of the subsamples used in the present analysis: full Lyα auto-correlation function (gray), low-redshift Lyα auto-correlation function (blue), high redshift Lyα auto-correlation function (red).The two latter subsamples are used to produce a measurement of H(z) at z = 2.19 and z = 2.49.

Fig
Fig. B.2.The 68% and 95% confidence level contours in the (α ,α ⊥ ) plane from the Lyα(Lyα) × Lyα(Lyα) computed with the low and high redshift bins.The ∆χ 2 values corresponding to confidence levels are taken from Table D.1.The black dot corresponds to the Pl2015 model.

Table 2 .
Parameters of the "Pl2015 model", i.e. the flat ΛCDM model of Planck Collaboration et al. (2016) that we use here to transform redshifts and angular separations into radial and transverse separations.

Table 3 .
The Lyα/metal and metal/metal pairs contributing to the flux correlation function.The table shows the ratio of transition wavelengths and the corresponding apparent separation, r ap ,

Table 6 .
Best fit values of (α , α ⊥ ) for the Lyα(Lyα) × Lyα(Lyα+Lyβ) correlation function fit with various models.The first group includes physical models starting with the basic Kaiser redshift-space model and then including, progressively, metals, HCD, and UV corrections.Fits in the second group include polynomial broadband terms, as described in the text.

Table 7 .
Best Errors on BAO parameters correspond to CL = 68.27%,while the other parameters have errors corresponding to ∆χ 2 = 1.

Table A .
1. Effective redshifts at which the α and α ⊥ parameters are measured.The average redshift of pairs is also given.Lyα(Lyα) × Lyα(Lyα) (low z) and (high z) are introduced in Appendix B.