Baryon acoustic oscillations from the cross-correlation of Ly{\alpha} absorption and quasars in eBOSS DR14

We present a measurement of the baryon acoustic oscillation (BAO) scale at redshift $z=2.35$ from the three-dimensional correlation of Lyman-$\alpha$ (Ly$\alpha$) forest absorption and quasars. The study uses 266,590 quasars in the redshift range $1.77<z<3.5$ from the Sloan Digital Sky Survey (SDSS) Data Release 14 (DR14). The sample includes the first two years of observations by the SDSS-IV extended Baryon Oscillation Spectroscopic Survey (eBOSS), providing new quasars and re-observations of BOSS quasars for improved statistical precision. Statistics are further improved by including Ly$\alpha$ absorption occurring in the Ly$\beta$ wavelength band of the spectra. From the measured BAO peak position along and across the line of sight, we determine the Hubble distance $D_{H}$ and the comoving angular diameter distance $D_{M}$ relative to the sound horizon at the drag epoch $r_{d}$: $D_{H}(z=2.35)/r_{d}=9.20\pm 0.36$ and $D_{M}(z=2.35)/r_{d}=36.3\pm 1.8$. These results are consistent at $1.5\sigma$ with the prediction of the best-fit flat $\Lambda$CDM cosmological model reported for the Planck (2016) analysis of CMB anisotropies. Combined with the Ly$\alpha$ auto-correlation measurement presented in a companion paper (de Sainte Agathe et al. 2019) the BAO measurements at $z=2.34$ are within $1.7\sigma$ of the predictions of this model.


Introduction
The baryon acoustic oscillation (BAO) peak in the cosmological matter correlation function at a distance corresponding to the sound horizon, r d ∼ 100h −1 Mpc, has been seen at several redshifts using a variety of tracers.Following the original measurements (Eisenstein et al. 2005;Cole et al. 2005), the most precise results have been obtained using bright galaxies in the redshift range 0.35 < z < 0.65 (Anderson et al. 2012(Anderson et al. , 2014b,a;,a;Alam et al. 2017) from the Baryon Oscillation Spectroscopy Survey (BOSS; Dawson et al. 2013) of the Sloan Digital Sky Survey-III (SDSS-III; Eisenstein et al. 2011).Other measurements using galaxies cover the range 0.1 < z < 0.8 (Percival et al. 2007(Percival et al. , 2010;;Beutler et al. 2011;Blake et al. 2011;Padmanabhan et al. 2012;Mehta et al. 2012;Chuang & Wang 2012;Xu et al. 2013;Ross et al. 2015;Bautista et al. 2018).At higher redshift, the peak has been seen in the correlation function of quasars at a mean redshift z ∼ 1.5 (Ata et al. 2018;Gil-Marín et al. 2018;Hou et al. 2018;Zarrouk et al. 2018) and in the flux-transmission E-mail: michael.blomqvist@lam.frcorrelation function in Lyman-α (Lyα) forests at z ∼ 2.3 (Busca et al. 2013;Slosar et al. 2013;Kirkby et al. 2013;Delubac et al. 2015;Bautista et al. 2017) and in the forest cross-correlation with quasars (Font-Ribera et al. 2014;du Mas des Bourboux et al. 2017).These observations all yield measurements of comoving angular-diameter distances and Hubble distances at the corresponding redshift, D M (z)/r d and D H (z)/r d = c/(H(z)r d ), relative to the sound horizon.
BAO measurements have found an important role in testing the robustness of the flat ΛCDM model that is consistent with observed CMB anisotropies (Planck Collaboration et al. 2016).While the parameters of this model are precisely determined by the CMB data by itself, more general models are less well constrained.Most significantly, adding BAO data improves constraints on curvature (Planck Collaboration et al. 2016).The addition of BAO and type Ia supernova (SN Ia) data (Betoule et al. 2014) generalizes the "CMB" measurement of H 0 , which assumes flatness, to give an "inverse-ladder" measurement of H 0 (Aubourg et al. 2015), that can be compared with distance-ladder measurements (Riess et al. 2016(Riess et al. , 2018a,b),b).Here, the inverse-

Data sample and reduction
The quasars and forests used in this study are drawn from SDSS DR14.This release includes data from DR12 taken in the first two generations SDSS-I/II, in the BOSS program of SDSS-III and in the eBOSS pilot program SEQUELS (Myers et al. 2015).These data were used in the measurement of the quasar-forest cross-correlation of dMdB17.Here, we use in addition data from the first two years of the eBOSS program and the completed SE-QUELS.
The catalog of identified quasars, DR14Q (Pâris et al. 2018), includes 266,590 quasars2 in the redshift range 1.77 < z q < 3.5.The distribution on the sky of these quasars is shown in Fig. 1 and the redshift distribution in Fig. 2.
All spectra used for this analysis were obtained using the BOSS spectrograph (Smee et al. 2013) on the 2.5 m SDSS telescope (Gunn et al. 2006) at Apache Point Observatory (APO).The spectrograph covers observed wavelengths 360.0 λ 1040.0 nm, with a resolving power R ≡ λ/∆λ FWHM increasing from ∼ 1300 to ∼ 2600 across the wavelength range.The data were processed by the eBOSS pipeline, the same (but a marginally updated version) as that used for the cross-correlation measurement of dMdB17.The pipeline performs wavelength calibration, flux calibration and sky subtraction of the spectra.The individual exposures (typically four of 15 minutes) of a given object are combined into a coadded spectrum that is rebinned onto pixels on a uniform grid with ∆ log 10 (λ) = 10 −4 (velocity width ∆v ≈ 69 km s −1 ).The pipeline additionally provides an automatic classification into object type (galaxy, quasar or star) and a redshift estimate by fitting a model spectrum (Bolton et al. 2012).
Visual inspection of quasar spectra was an important procedure during the first three generations of SDSS to correct for misclassifications of object type and inaccurate redshift determinations by the pipeline (Schneider et al. 2010;Pâris et al. 2017).Starting in SDSS-IV, most of the objects are securely classified by the pipeline, with less than 10% of the spectra requiring visual inspection (Dawson et al. 2016).The visual-inspection redshifts, when available, are taken as the definitive quasar redshifts, while the remaining quasars have redshifts estimated by the pipeline.
The cross-correlation analysis presented here involves the selection of three quasar samples from DR14Q: tracer quasars (for which we only need the redshifts and positions on the sky), quasars providing Lyα forest absorption in the Lyα region, and quasars providing Lyα forest absorption in the Lyβ region.
The selected sample of tracer quasars contains 266,590 quasars in the range 1.77 < z q < 3.5.It includes 13,406 SDSS DR7 quasars (Schneider et al. 2010) and 18,418 broad absorption line (BAL) quasars, the latter identified as having a CIV balnicity index (Weymann et al. 1991) BI_CIV>0 in DR14Q.Quasars with redshifts less than 1.77 are excluded because they are necessarily separated from observable forest pixels (see below) by more than 200 h −1 Mpc, the maximum distance where the correlation function is measured.The upper limit of z q = 3.5 is adopted because of the low number of higher-redshift quasars that both limits their usefulness for correlation measurements and make them subject to contamination due to redshift errors of the much more numerous low-redshift quasars (Busca & Balland 2018).Such contaminations would be expected to add noise (but not signal) to the cross-correlation.
The summary of the Lyα forest data covering the Lyα or Lyβ region of the quasar spectrum is given in Table 1.Both samples exclude SDSS DR7 quasars and BAL quasars.The Lyα sample is derived from a super set consisting of 194,166 quasars in the  redshift range 2.05 < z q < 3.5, whereas the Lyβ sample is taken from a super set containing 76,650 quasars with 2.55 < z q < 3.5.The lower redshift limits are a consequence of the forests exiting the wavelength coverage of the spectrograph for quasars with z q < 2 and z q < 2.53, respectively.Spectra with the same object identification THING_ID (re-observed quasars) are coadded using inverse-variance weighting.For the selected forest samples, 17% of the quasars have duplicate spectra (less than 2% have more than one reobservation) taken with the BOSS spectrograph.
The forest spectra are prepared for analysis by discarding pixels which were flagged as problematic in the flux calibra-tion or sky subtraction by the pipeline.We mask pixels around bright sky lines using the condition 10 4 log 10 (λ/λ sky ) ≤ 1.5, where λ sky is the wavelength at the pixel center of the sky line where the pipeline sky subtraction is found to be inaccurate.Finally, we double the mask width to remove pixels around the observed CaII H&K lines arising from absorption by the interstellar medium of the Milky Way.
Forests featuring identified damped Lyα systems (DLAs) are given a special treatment using an updated (DR14) version of the catalog of Noterdaeme et al. (2012).All pixels where the absorption due to the DLA is greater than 20% are masked.The absorption in the wings is corrected using a Voigt profile following the procedure of Noterdaeme et al. (2012).
To facilitate the computation of the cross-correlation, we follow the approach in Bautista et al. (2017) to combine three adjacent pipeline pixels into wider "analysis pixels" defined as the inverse-variance-weighted flux average.Requiring a minimum of 20 analysis pixels in each spectrum discards 2447 (6155) forests for the Lyα (Lyβ) region.Lastly, 3087 (1882) forests failed the continuum-fitting procedure (see section 3) for the Lyα (Lyβ) region by having negative continua due to their low spectral signal-to-noise ratios.The final samples include 188,632 forests for the Lyα region and 68,613 forests for the Lyβ region.Figure 2 shows the redshift distributions for the tracer quasars and the Lyα absorption pixels.Our samples can be compared to those of dMdB17, which included 234,367 quasars (217,780 with 1.8 < z q < 3.5) and 168,889 forests (157,845 with 2.0 < z q < 3.5) over a wider redshift range.

The Lyα forest flux-transmission field
The transmitted flux fraction F in a pixel of the forest region of quasar q is defined as the ratio of the observed flux density f q with the continuum flux C q (the flux density that would be observed in the absence of absorption).We will be studying the transmission relative to the mean value at the observed wave-  length F(λ), and refer to this quantity as the "delta-field": We employ a similar method to the one established by previous Lyα forest BAO analyses (Busca et al. 2013;Delubac et al. 2015) in which the delta-field is derived by estimating the product C q (λ)F(λ) for each quasar.Each spectrum is modeled assuming a uniform forest spectral template which is multiplied by a quasar-dependent linear function, setting the overall amplitude and slope, to account for the diversity of quasar luminosity and spectral shape: where a q and b q are free parameters fitted to the observed flux of the quasar.The forest spectral template f (λ rf ) is derived from the data as a weighted mean normalized flux, obtained by stacking the spectra in the quasar rest-frame.The continuum fitting procedure is handled separately for the Lyα and Lyβ regions.The total variance of the delta-field is modeled as where the noise variance σ 2 noise = σ 2 pipe /(C q F) 2 .The first term represents the pipeline estimate of the flux variance, corrected  (Planck Collaboration et al. 2016).The sound horizon at the drag epoch, r d , is calculated using CAMB (Lewis et al. 2000).The Hubble distance D H and the comoving angular diameter distance D M relative to r d are given at the effective redshift of the measurement z eff .

Parameter
Value by a function η(λ) that accounts for possible misestimation.The second term gives the contribution due to the large-scale structure (LSS) and acts as a lower limit on the variance at high signal-to-noise ratio.Lastly, the third term absorbs additional variance from quasar diversity apparent at high signal-to-noise ratio.In bins of σ 2 noise and observed wavelength, we measure the variance of the delta-field and fit for the values of η, σ 2 LSS and as a function of observed wavelength.These three functions are different for the Lyα and Lyβ regions.The procedure of stacking the spectra, fitting the continua and measuring the variance of δ is iterated, until the three functions converge.We find that five iterations is sufficient.Figure 3 presents an example spectrum and the best-fit model C q (λ)F(λ) for the Lyα and Lyβ regions.
As detailed in Bautista et al. (2017), the delta-field can be redefined in two steps to make exact the biases introduced by the continuum fitting procedure.In the first step, we define where the over-bars refer to weighted averages over individual forests.Next, we transform the δq (λ) by subtracting the weighted average at each observed wavelength: (5)

The Lyα forest -quasar cross-correlation
The three-dimensional positions of the quasars and the Lyα forest delta-field are determined by their redshifts and angular positions on the sky.We transform the observed angular and redshift separations (∆θ, ∆z) of the quasar-Lyα absorption pixel pairs into Cartesian coordinates (r ⊥ , r ) assuming a spatially flat fiducial cosmology.The comoving separations along the line of sight r (parallel direction) and transverse to the line of sight r ⊥ (perpendicular direction) are calculated as where D α ≡ D c (z α ) and D q ≡ D c (z q ) are the comoving distances to the Lyα absorption pixel and the quasar, respectively.Line of sight separations r > 0 (< 0) thus correspond to background (foreground) absorption with respect to the tracer quasar position.In this paper, we will also refer to the coordinates (r, µ), where r 2 = r 2 + r 2 ⊥ and µ = r /r, the cosine of the angle of the vector r from the line of sight.The pair redshift is defined as z pair = (z α + z q )/2.A histogram of the pair redshifts is displayed in Fig. 4. We do not include pairs involving a quasar and pixels from its own forest in the cross-correlation analysis, because the correlation of such pairs vanishes due to the continuum fit and delta-field redefinition (eqn 4).
The fiducial cosmology used in the analysis is a flat ΛCDM model with parameter values taken from the Planck (2016) result for the TT+lowP combination (Planck Collaboration et al. 2016) described in Table 2.It is the same fiducial cosmology employed by dMdB17.

Cross-correlation
We estimate the cross-correlation at a separation bin A, ξ A , as the weighted mean of the delta-field in pairs of pixel i and quasar k at a separation within the bin The weights w i are defined as the inverse of the total pixel variance (see equation 3), multiplied by redshift evolution factors for the forest and quasar, so as to approximately minimize the relative error on ξA (Busca et al. 2013): where γ α = 2.9 (McDonald et al. 2006) and γ q = 1.44 (du Mas des Bourboux et al. 2019).The validity of the correlation estimator, as well as the accuracy of the distortion matrix (section 4.2) and covariance matrix estimation (section 4.3) were tested and confirmed on simulated data in dMdB17.
Our separation grid consists of 100 bins of 4 h −1 Mpc for separations r ∈ [−200, 200] h −1 Mpc in the parallel direction and 50 bins of 4 h −1 Mpc for separations r ⊥ ∈ [0, 200] h −1 Mpc in the perpendicular direction; the total number of bins is N bin = 5000.Each bin is defined by the weighted mean (r ⊥ , r ) of the quasarpixel pairs of that bin, and its redshift by the weighted mean pair redshift.The mean redshifts range from z = 2.29 to z = 2.40.The effective redshift of the cross-correlation measurement is defined to be the inverse-variance-weighted mean of the redshifts of the bins with separations in the range 80 < r < 120 h −1 Mpc around the BAO scale.Its value is z eff = 2.35.
Because the Lyβ transition is sufficiently separated in wavelength from the Lyα transition, corresponding to large physical separations > 441 h −1 Mpc for the wavelength range of the analysis, we neglect the contamination from Lyβ absorption interpreted as Lyα absorption.The total number of pairs of the crosscorrelation measurement is 9.7 × 10 9 .The Lyα absorption in the Lyβ region contributes 1.2 × 10 9 pairs (13%) and reduces the mean variance of the correlation function by 9% compared to the Lyα region-only measurement.Our cross-correlation measurement has 39% lower mean variance than the measurement of dMdB17.

Distortion matrix
The procedure used to estimate the delta-field (section 3) suppresses fluctuations of characteristic scales corresponding to the forest length, since the estimate of the product CF (equation 1) would typically erase such a fluctuation.The result is a suppression of power in the radial direction on large scales and a significant distortion of the correlation function.As first noted in Slosar et al. (2011) and further investigated in Blomqvist et al. (2015), the distortion effect can be modeled in Fourier space as a multiplicative function of the radial component k on the Lyα forest transmission power spectrum.The suppression scale was treated as a free parameter to be determined in the fit of the correlation function.
We use the method introduced by Bautista et al. ( 2017) for the Lyα auto-correlation and adapted to the cross-correlation by dMdB17 which allows one to encode the effect of this distortion on the correlation function in a "distortion matrix".This approach was extensively validated in these publications using simulated data.The measured correlation function for a separation bin A is given as a linear combination of the true correlation function for bins A : ξA The distortion matrix D AA depends only on the geometry of the survey, the lengths of the forests and the pixel weights, where the projection matrix and δ K is the Kronecker delta.The indices i and j in equation (11) refer to pixels from the same forest, k refers to a quasar, and the sums run over all pixel-quasar pairs that contribute to the separation bins A and A .The diagonal elements dominate the distortion matrix and are close to unity, D AA ≈ 0.97, whereas the off-diagonal elements are small, |D AA | 0.03.We use the distortion matrix when performing fits of the measured crosscorrelation function (see equation 16).

Covariance matrix
We estimate the covariance matrix of the cross-correlation from the data by using the subsampling technique introduced by Busca et al. (2013) and adapted to the cross-correlation by dMdB17.We divide the DR14 footprint of Figure 1 into subsamples and measure the covariance from the variability across the subsamples.Such estimates of the covariance matrix are unbiased, but the noise due to the finite number of subsamples leads to biases in the inverse of the covariance (Joachimi & Taylor 2014).As was done in dMdB17, we smooth the noise by assuming, to good approximation, that the covariance between separation bins A and B depends only on the absolute difference (∆r ).We define the subsamples through a HEALPix (Górski et al. 2005) pixelization of the sky.A quasar-absorption pixel pair is assigned to a subsample s if the forest that contains the absorption belongs to that HEALPix pixel.We use HEALPix parameter nside=32, resulting in 3262 subsamples.Using fewer but larger HEALPix pixels (nside=16, 876 subsamples) has no significant impact on the covariance matrix or the BAO peak position measurement.
The (noisy) covariance matrix is calculated as where the sum runs over all subsamples and W A is the sum of the pair weights w belonging to bin A, From the covariance, we calculate the correlation matrix: The smoothing procedure is applied to this correlation matrix by averaging as a function of (∆r , ∆r ⊥ ).The final covariance used in the fits is obtained by multiplying the smoothed correlation matrix by the diagonal elements of the original covariance matrix.Figure 5 displays the smoothed correlation matrix as a function of ∆r for the three lowest values of ∆r ⊥ .

Model of the cross-correlation
We fit the measured cross-correlation function, ξA , in the (r ⊥ , r ) bin A, to a cosmological correlation function ξ cosmo where D AA is the distortion matrix (equation 11).The broadband term, ξ bb A , is an optional function used to test for imperfections in the model and for systematic errors.The set of parameters for the model is summarized in Table 3.The model is calculated at the weighted mean (r ⊥ , r ) and redshift of each bin of the correlation function.Because of the relatively narrow redshift distribution of the bins (∆z = 0.11), most model parameters can be assumed as redshift independent to good accuracy.
The cosmological cross-correlation function is the sum of several contributions The first term represents the standard correlation between quasars, q, and Lyα absorption in the IGM.It is the most important part of the correlation function and, used by itself, would lead to an accurate determination of the BAO peak position (see results in section 6).
The remaining terms in equation ( 17) represent subdominant effects but contribute towards improving the fit of the correlation function outside the BAO peak.The second term is the sum over correlations from metal absorbers in the IGM.The third term represents Lyα absorption by high column density systems (HCDs).The fourth term is the correlation from the effect of a quasar's radiation on a neighboring forest ("transverse proximity effect").The fifth term is a relativistic correction leading to odd-multipoles in the correlation function, and the final term includes other sources of odd-multipoles (Bonvin et al. 2014).These terms will be described in detail below.

Quasar-Lyα correlation term
The quasar-Lyα cross-correlation, ξ qα , is the dominant contribution to the cosmological cross-correlation.It is assumed to be a biased version of the total matter auto-correlation of the appropriate flat ΛCDM model, separated into a smooth component and a peak component to free the position of the BAO peak: where A peak is the BAO peak amplitude.The anisotropic shift of the observed BAO peak position relative to the peak position of the fiducial cosmological model from Table 2 is described by the line-of-sight and transverse scale parameters The nominal correlation function, ξ qα (r ⊥ , r , α ⊥ = α = 1), is the Fourier transform of the quasar-Lyα cross-power spectrum: where k = (k , k ⊥ ) is the wavenumber of modulus k with components k along the line of sight and k ⊥ across, and µ k = k /k is the cosine of the angle of the wavenumber from the line of sight.As described in detail below, P QL is the (quasi) linear matter spectrum, d q and d Lyα are the standard linear-theory factors describing the tracer bias and redshift-space distortion (Kaiser 1987), V NL describes further non-linear corrections not included in P QL , and G(k) gives the effects of (r ⊥ , r ) binning on the measurement.
The first term in (20) provides for the aforementioned decoupling of the peak component (Eq.18): where the smooth component, P sm , is derived from the linear power spectrum, P L (k, z), via the side-band technique Transverse binning smoothing parameter A peak = 1 BAO peak amplitude γ α = 2.9 Lyα transmission bias evolution exponent γ m = 1 Metal transmission bias evolution exponent (Kirkby et al. 2013) and P peak = P L − P sm .The redshiftdependent linear power spectrum is obtained from CAMB (Lewis et al. 2000) with the fiducial cosmology.
The correction for non-linear broadening of the BAO peak is parameterized by Σ = (Σ , Σ ⊥ ), with Σ ⊥ = 3.26 h −1 Mpc and where f = d(ln g)/d(ln a) ≈ Ω 0.55 m (z) is the linear growth rate of structure.
The second term in (20) describes the quasar bias and redshift-space distortion Because the fit of the cross-correlation is only sensitive to the product of the quasar and Lyα biases, we set b q ≡ b q (z eff ) = 3.77 and assume a redshift dependence of the quasar bias given by (Croom et al. 2005) The quasar redshift-space distortion, assumed to be redshift independent, is Setting f = 0.969 for our fiducial cosmology yields β q = 0.257.
The third term in ( 20) is the Lyα forest factor, We assume that the transmission bias evolves with redshift as with γ α = 2.9 (McDonald et al. 2006), while β α is assumed to be redshift independent.We choose to fit for β α and the velocity gradient bias of the Lyα forest: Beyond our standard treatment of the Lyα transmission bias, we also consider the effect of fluctuations of ionizing UV radiation which lead to a scale-dependence of b α (Pontzen 2014; where W(x) = arctan(x)/x (following the parameterization of Gontcho A Gontcho et al. 2014).Our standard fit does not include the effect of UV fluctuations due to its minor contribution to the fit quality.A fit that includes the UV modeling is presented in Table A.1 for which we fix the UV photon mean free path λ UV = 300 h −1 Mpc (Rudie et al. 2013) and b a = −2/3 (Gontcho A Gontcho et al. 2014), and fit for b Γ , as was done in dMdB17.
The effect of quasar non-linear velocities and statistical redshift errors on the power spectrum is modeled as a Lorentz damping (Percival & White 2009), where σ v is a free parameter.
The last term in (20), G(k), accounts for smoothing due to the binning of the measurement of the correlation function (Bautista et al. 2017).We use where R and R ⊥ are the scales of the smoothing.In the transverse direction, this form is not exact, but we have verified that it generates a sufficiently accurate correlation function.We fix both to the bin width, R = R ⊥ = 4 h −1 Mpc.Systematic errors in the quasar redshift estimates lead to a shift of the cross-correlation along the line of sight which is accounted for in the fit using the free parameter ∆r = r ,true − r ,measured = (1 + z)∆v H(z) . (32)

Quasar-metal correlation terms
Absorption by metals in the intergalactic medium (e.g., Pieri et al. 2014) with similar rest-frame wavelengths to Lyα yields a sub-dominant contribution to the measured cross-correlation.
Assuming that these contaminant absorptions have redshifts corresponding to Lyα absorption results in an apparent shift of the quasar-metal cross-correlations along the line of sight in the where is a "metal distortion matrix" that allows us to calculate the shifted quasar-metal cross-correlation function for a given nonshifted quasar-metal cross-correlation function.The condition (i, k) ∈ A refers to pixel distances calculated using z α , but (i, k) ∈ B refers to pixel distances calculated using z m .For each metal absorption line, the (non-shifted) quasar-metal correlation is modeled using (20) with d α replaced by The metal absorption lines included in the fit are listed in Table 4.
Because the redshift-space distortion parameter of each metal is poorly determined in the fit, we fix β m = 0.5, the value derived for DLA host halos (Font-Ribera et al. 2012;Pérez-Ràfols et al. 2018).Transmission biases are assumed to evolve with redshift as a power-law with exponent γ m = 1, similar to the measured evolution of the CIV bias (Blomqvist et al. 2018), but our results are not sensitive to this choice.

Other correlation terms
Correlations due to absorption by unidentified HCD systems are modeled using ( 20) with d α replaced by (Rogers et al. 2018) where the bias b Lyα−quasar HCD and the redshift-space distortion β HCD are free parameters in the fit.Due to degeneracies, we add a Gaussian prior on β HCD of mean 0.5 and standard deviation 0.2.The parameter L HCD is the characteristic size of unidentified HCD systems.Our DLA-identification procedure requires their width (wavelength interval for absorption greater than 20% ) to be above ∼ 2.0 nm, corresponding to ∼ 14 h −1 Mpc.Following the study of Rogers et al. (2018), the corresponding unidentified HCD systems are well-modeled with L HCD = 10h −1 Mpc and we fix L HCD to this value in the fits.We have verified that varying this parameter over the range 5 < L HCD < 15h −1 Mpc does not change the fitted position of the BAO peak.
The term in (17) representing the transverse proximity effect takes the form (Font-Ribera et al. 2013): This form supposes isotropic emission from the quasars.We fix λ UV = 300 h −1 Mpc (Rudie et al. 2013) and fit for the amplitude ξ TP 0 .
In addition to accounting for asymmetries in the crosscorrelation introduced by metal absorptions, continuum-fitting distortion and systematic redshift errors, the standard fit includes modeling of relativistic effects (Bonvin et al. 2014).The relativistic correction in ( 17) is the sum of two components describing a dipole and an octupole, where L 1 and L 3 are the Legendre polynomial of degree 1 and 3 respectively, A rel1 and A rel3 are the amplitudes, and where j is the spherical Bessel function.The relativistic dipole is expected to be the dominant contribution of odd-asymmetry and our standard fit therefore neglects the relativistic octupole (A rel3 = 0).Dipole and octupole asymmetries also arise in the "standard" correlation function due to the evolution of the tracer biases and growth factor, as well as from the wide-angle correction (Bonvin et al. 2014): where Here, the two amplitudes A asy0 and A asy2 determine the dipole contribution, while A asy3 is the octupole amplitude.The ξ asy term is neglected in the standard fit, but we check the robustness of the BAO measurement with respect to the odd-multipoles in Table A.1.

Broadband function
The optional ξ bb term of ( 16) is a "broadband function" that is a slowly varying function of (r , r ⊥ ): where L j is the Legendre polynomial of degree j.Its purpose is to account for unknown physical, instrumental or analytical effects missing in the model that could potentially impact the BAO measurement.The standard fit features no broadband function.The result of adding a broadband function of the form (i min , i max , j min , j max ) = (0, 2, 0, 6) is presented in Table A.1.

Fits of the cross-correlation
Our "standard" fit of the cross-correlation function uses the 14 parameters in the first group of Table 3.The fit includes 3180 data bins in the range 10 < r < 180 h −1 Mpc.The best-fit values are presented in the column "Lyα-quasar" of Table 5. Figure 6 shows the best fit for four ranges of µ and Figure 7 for the two lowest r ⊥ bins.Constraints on the BAO parameters (α ⊥ , α ) are presented in Fig. 8.Following the method introduced and described in detail in dMdB17, we estimate the relation between ∆χ 2 = χ 2 − χ 2 min and confidence levels 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 6, indicate that the (68.27,95.45%)confidence levels for (α ⊥ , α ) correspond to ∆χ 2 = (2.51,6.67) (instead of the nominal values ∆χ 2 = (2.3,6.18)).These levels are shown as contours in Fig. 8.The best-fit values and confidence level (68.27,95.45%)ranges are: corresponding to These results are consistent at 1.5 standard deviations with the prediction of the Planck (2016) best-fit flat ΛCDM model.Using a model without the BAO peak (A peak = 0) degrades the quality of the fit by ∆χ 2 = 22.48.
Our BAO constraints can be compared with the DR12 measurement of dMdB17 at a slightly higher redshift: D M (2.40)/r d = 35.7 ± 1.7 and D H (2.40) = 9.01 ± 0.36 corresponding to α ⊥ = 0.898 ± 0.042 and α = 1.077 ± 0.042, relative to the same Planck model.The results ( 43) and (44) thus represent a movement of ∼ 0.3σ toward the Planck-inspired model through a shift in α ⊥ .As a cross-check of the results, we apply our analysis to the DR12 data set of dMdB17, without including the absorption in the Lyβ region.The best-fit values are α ⊥ = 0.889 ± 0.040 and α = 1.080 ± 0.039 (errors correspond to ∆χ 2 = 1), in good agreement with the measurement of dMdB17.This result indicates that the movement towards the fiducial model in DR14 is driven by the data.
Model predictions for D M /r d and D H /r d depend both on prerecombination physics, which determine r d , and on late-time physics, which determine D M and D H . Taking the ratio, yielding the Alcock-Paczyński parameter F AP = D M /D H (Alcock & Paczynski 1979), isolates the late-time effects which, in the ΛCDM model depend only on (Ω m , Ω Λ ).We find where the ∆χ 2 curve is shown in Fig. 9 and we have adopted that the (68.27,95.45%)confidence levels correspond to ∆χ 2 = (1.13,4.74) (instead of the nominal values ∆χ 2 = (1, 4)).This result is 1.8 standard deviations from the prediction of the Planckinspired model, F AP (z = 2.35) = 4.60.
The fit of the cross-correlation prefers a vanishing contribution from the quasar-HCD correlation term (b HCD ≈ 0).This preference is in contrast to the Lyα auto-correlation of  2019) where the HCD model is a crucial element to obtain a good fit (but does not affect the BAO peak position measurement).The best-fit radial coordinate shift ∆r is consistent with zero systematic redshift error, but the parameter is strongly correlated with the amplitude of the relativistic dipole.Setting A rel1 = 0 in the fit yields ∆r = −0.92± 0.12, in good agreement with the value reported in dMdB17.The best fit suggests marginal support for a non-zero value of A rel1 , and the combined fit increases the significance of this result.However, even with sufficient statistical significance, its correlation with ∆r (as well as other potential systematic errors on A rel1 ) prevents claims of a discovery of relativistic effects.The parameters β α , σ v and ξ TP 0 all have best-fit values in agreement with the result of dMdB17.Among the metals, only SiIII(120.7)has a bias parameter significantly different from zero (> 4σ) to show evidence for large-scale correlations with quasars.The imprints of the metal correlations are visible in the line-of-sight direction in Figure 7.
Besides the standard approach, we have also performed nonstandard analyses, described in Appendix A, to search for unexpected systematic errors in the BAO peak-position measurement.The results for the non-standard fits of the cross-correlation are summarized in Table A.1.No significant changes of the best-fit values of (α ⊥ , α ) are observed.We have also divided the data to perform fits of the cross-correlation for a low-and a high-redshift bin as described in Appendix B. These fits are summarized in Table B.1 and yield consistent best-fit BAO parameters for the two bins.

Combination with the Lyα auto-correlation
We combine our measurement of the Lyα-quasar crosscorrelation with the DR14 Lyα auto-correlation of de Sainte Agathe et al. ( 2019) by performing a combined fit of the correlation functions.For the auto-correlation, we use the combination Lyα(Lyα)xLyα(Lyα) + Lyα(Lyα)xLyα(Lyβ) described in de Sainte Agathe et al. (2019).Because the covariance between the auto-and cross-correlation is sufficiently small to be ignored, as studied in Delubac et al. (2015) and dMdB17, we treat their errors as independent.The combined fit uses the standard fit model of each analysis.In addition to the 14 free parameters in Table 3, we let free the redshift-space distortion of quasars, β q , and the auto-correlation fit introduces three additional bias parameters (b CIV(154.9) , b ), for a total of 18 free parameters.The effective redshift of the combined fit is z eff = 2.34.
The best-fit results are presented in the column "combined" of Table 5. Figure 8 displays the constraints on (α ⊥ , α ) from the combined measurement as black contours indicating the (68.27,95.45%)confidence levels (corresponding to  2019).The fits are over the range 10 < r < 180 h −1 Mpc.Errors on BAO parameters correspond to CL = 68.27%,while the other parameters have errors corresponding to ∆χ 2 = 1.The parameter β q is fixed for the cross-correlation fit.The bottom section of the table gives the minimum χ 2 , the number of data bins (N bin ) and free parameters (N param ) in the fit, the probability, the effective redshift, the correlation coefficient (ρ) for the BAO parameters, and the χ 2 for the fit with the fixed fiducial BAO peak position.

Implications for cosmological parameters
The combined-fit measurement of (D M /r d , D H /r d ) at z = 2.34 presented here is within 1.7 standard deviations of the predictions of the flat ΛCDM model favored by the measurement of CMB anisotropies (Planck Collaboration et al. 2016).This result thus does not constitute statistically significant evidence for new physics or unidentified systematic errors in the measurement.Figure 10 illustrates the agreement with the Planck prediction for the ensemble of BAO measurements.A value of H 0 can be obtained either by using the CMB measurement of r d or by using the Primordial-Nucleosynethsis value of Ω b h 2 to constrain r d .Adopting the value Ω b h 2 = 0.02156 ± 0.002 derived from the deuterium abundance measurement of Cooke et al. (2018) and assuming flat ΛCDM, we derive the constraints on (H 0 , Ω m h 2 ) shown in Figure 12 with h = 0.678 ± 0.020 (54) or h < 0.72 at 95% C.L. The limit degrades to h < 0.735 (95% C.L.) if one adopts a more conservative uncertainty on the baryon density: Ω b h 2 = 0.02156 ± 0.003.Nevertheless, as previously noted (Aubourg et al. 2015;Addison et al. 2018), the combination of BAO and nucleosynthsis provides a CMB-free confirmation of the tension with the distance-ladder determinations of H 0 (Riess et al. 2016(Riess et al. , 2018a,b),b).

Conclusions
Using the entirety of BOSS and the first two years of eBOSS observations from SDSS DR14, this paper has presented a measurement of the cross-correlation of quasars and the Lyα flux transmission at redshift 2.35.In addition to the new and reobserved quasars provided in DR14, we have improved statistics further by extending the Lyα forest to include Lyα absorption in the Lyβ region of the spectra.The position of the BAO peak is 1.5σ from the flat ΛCDM model favored by CMB anisotropy measurements (Planck Collaboration et al. 2016).We emphasize that the measured peak position shows no significant variation when adding astrophysical elements to the fit model.The basic Lyα-only model on its own provides an accurate determination of the peak position, while still yielding an acceptable fit to the data.Compared to the BAO measurement for the DR12 data set reported by dMdB17, our result represents a movement of ∼ 0.3σ toward the Planck-cosmology prediction through a shift in the transverse BAO parameter α ⊥ .This change is driven by the data and not by differences in the analyses.The inclusion of Lyα absorption in the Lyβ region has no impact on the best-fit value of α ⊥ .Combined with the Lyα-flux-transmission auto-correlation measurement presented in a companion paper (de Sainte Agathe et al. 2019), the BAO peak at z = 2.34 is 1.7σ from the expected value.
The ensemble of BAO measurements is in good agreement with the CMB-inspired flat ΛCDM model.By themselves, the BAO data provide a good confirmation of this model.The use of SNIa to measure cosmological distances (Scolnic et al. 2018) provides independent measurements of the model parameters.As can be seen in Fig. 11 they are in agreement with the BAO measurements.
The BAO measurements presented here will be improved by the upcoming DESI (DESI Collaboration et al. 2016) and WEAVE-QSO (Pieri et al. 2016) projects both by increasing the number of quasars and improving the spectral resolution.
The best-fit results and the χ 2 scans for the cross-correlation by itself and the combination with the auto-correlation are publicly available. 3 Table A.1.Results of non-standard fits.The first group presents results of successively adding complications from physical effects to the basic Lyα-only model.These complications are: metals, absorption by high-column density systems, the transverse proximity effect, and the relativistic dipole, corresponding to the standard fit from column 1 of Table 5.The second group presents fits which include fluctuations of the UV background radiation, the odd multipoles = (1, 3) or the broadband function (for this group we set ξ qHCD = 0).The last group presents fits for non-standard data samples: no absorption in the Lyβ region or no correction of DLAs in the spectra.The fit is over the range 10 < r < 180 h −1 Mpc.Errors correspond to ∆χ 2 = 1.For these fits, we set ξ qHCD = 0 to facilitate the parameter error estimation with no impact on the best fits.No significant changes of the BAO parameters are observed.The third group of fits concern non-standard data samples that either omits the correlation pairs from the Lyβ region ("no Lyβ") or leaves the DLAs uncorrected in the spectra ("keep DLAs").Even for these modified data samples the best-fit values of (α ⊥ , α ) do not deviate significantly from those of our standard analysis.

Fig. 1 .
Fig.1.Sky distribution for the sample of 266,590 tracer quasars (1.77 < z q < 3.5) from DR14Q in J2000 equatorial coordinates.The solid black curve is the Galactic plane.The high-density regions are the eBOSS and SEQUELS observations (for the northern regions of the two Galactic hemispheres) and SDSS-stripe 82 (for declination δ ∼ 0).The discontiguous small areas contain only SDSS DR7 quasars.

Fig. 3 .
Fig. 3. Example spectrum of the DR14Q quasar identified by (Plate, MJD, FiberID) = (7305, 56991, 570) at z q = 3.0.The blue line indicates the best-fit model F(z)C q (λ) for the Lyα region covering the rest-frame wavelength interval 104.0 < λ rf < 120.0 nm.The red line indicates the same for the Lyβ region over the range 97.4 < λ rf < 102.0 nm.The Lyα and Lyβ emission lines are located at λ α = 121.567nm and λ β = 102.572nm in the quasar rest-frame.The spectrum has not been rebinned into analysis pixels in this figure.

Fig. 4 .
Fig. 4. Redshift distribution of the 9.7×10 9 correlation pairs.The dashed vertical black line indicates the effective redshift of the BAO measurement, z eff = 2.35, calculated as the weighted mean of the pair redshifts for separations in the range 80 < r < 120 h −1 Mpc.

Fig. 5 .
Fig. 5.The smoothed correlation matrix from sub-sampling as a function of ∆r = |r ,A − r ,B |.The curves are for constant ∆r ⊥ = |r ⊥,A − r ⊥,B | for the three lowest values ∆r ⊥ = [0, 4, 8] h −1 Mpc.The right panel shows an expansion of the region ∆r < 140 h −1 Mpc.Table3.List of the parameters of the cross-correlation model.The 14 parameters of the standard fit are given in the first section of the table.The second section lists parameters that are fixed in the standard fit.

Fig. 6 .
Fig. 6.The r 2 -weighted cross-correlation function averaged in four ranges of µ = r /r.The red curves show the best-fit model of the standard fit obtained for the fitting range 10 < r < 180 h −1 Mpc.The curves have been extrapolated outside this range.

Fig. 7 .
Fig. 7.The cross-correlation function as a function of r for the two lowest values r ⊥ = [2, 6] h −1 Mpc.The red curves indicate the best-fit model of the standard fit obtained for the fitting range 10 < r < 180 h −1 Mpc.The curves have been extrapolated outside this range.The imprints of quasarmetal correlations are visible as peaks indicated by the dashed black lines at r ≈ −21 h −1 Mpc (SiIII(120.7)),r ≈ −53 h −1 Mpc (SiII(119.0)),r ≈ −59 h −1 Mpc (SiII(119.3)),and r ≈ +103 h −1 Mpc (SiII(126.0)).Table 5. Fit results for the cross-correlation, the auto-correlation of de Sainte Agathe et al. (2019), and the combined fit.The auto-correlation fit uses the combination Lyα(Lyα)xLyα(Lyα) + Lyα(Lyα)xLyα(Lyβ) described in de Sainte Agathe et al. (2019).The fits are over the range 10 < r < 180 h −1 Mpc.Errors on BAO parameters correspond to CL = 68.27%,while the other parameters have errors corresponding to ∆χ 2 = 1.The parameter β q is fixed for the cross-correlation fit.The bottom section of the table gives the minimum χ 2 , the number of data bins (N bin ) and free parameters (N param ) in the fit, the probability, the effective redshift, the correlation coefficient (ρ) for the BAO parameters, and the χ 2 for the fit with the fixed fiducial BAO peak position.

Fig. 9 .
Fig. 9. Constraints on the Alcock-Paczyński parameter F AP for the cross-correlation (red) and the combination with the auto-correlation (black).Confidence levels of (68.27%, 95.45%) are indicated with the horizontal dotted lines for the cross-correlation and dashed lines for the combined fit.The prediction of the Planck (2016) bestfit flat ΛCDM cosmology is indicated with the vertical dotted line at F AP (z = 2.35) = 4.60 for the cross-correlation and dashed line at F AP (z = 2.34) = 4.57 for the combined fit.

Fig. 11 .
Fig. 11.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 de Sainte Agathe 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 blue contours do not use the Lyα auto-correlation measurement of de Sainte Agathe 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.

Fig. 12 .
Fig. 12.One and two standard deviation constraints on H 0 and Ω m h 2 derived the BAO data used in Fig. 11 and from Big-Bang Nucleosynthesis.This figure assumes a flat universe and a Gaussian prior Ω b h 2 = 0.02156 ± 0.002 derived from the deuterium abundance measurement of Cooke et al. (2018).

Table 1 .
Definition of the Lyα and Lyβ regions of the quasar spectrum in which we measure Lyα forest absorption.The table shows the restand observer-frame wavelength ranges defining the regions, the range of quasar redshifts, and the number of forests available in our analysis sample.

Table 2 .
Parameters of the flat ΛCDM fiducial cosmological model

Table 4 .
The most important metal absorptions of the intergalactic medium that imprint correlations observed in the Lyα-quasar crosscorrelation for r ∈[−200, 200]h −1 Mpc.The second column lists the rest-frame wavelength of the metal line and the third column its ratio with λ α (using the shorter of the two wavelengths in the denominator).The last column gives the apparent radial distance difference between the Lyα and metal absorption, r = D c (z α ) − D c (z m ), for observed wavelength λ = 407.2nm (corresponding to Lyα absorption at z eff = 2.35).

Table 6 .
Values of ∆χ 2 corresponding to CL = (68.27,95.45%).Values are derived from 10,000 Monte Carlo simulations of the correlation function that are fit using the model containing only Lyα absorption.Confidence levels are the fractions of the generated data sets that have best fits below the ∆χ 2 limit.The uncertainties are statistical and estimated using bootstrap.