# The Completed SDSS-IV Extended Baryon Oscillation Spectroscopic Survey: Baryon Acoustic Oscillations with Lyα Forests

Helion Du Mas Des Bourboux, James Rich, Andreu Font-Ribera, Victoria de Sainte Agathe, James Farr, Thomas Etourneau, Jean-Marc Le Goff, Andrei Cuceu, Christophe Balland, Julian E. Bautista, et al.

► To cite this version:

Helion Du Mas Des Bourboux, James Rich, Andreu Font-Ribera, Victoria de Sainte Agathe, James Farr, et al.. The Completed SDSS-IV Extended Baryon Oscillation Spectroscopic Survey: Baryon Acoustic Oscillations with Lyα Forests. The Astrophysical Journal, 2020, 901 (2), pp.153. 10.3847/1538-4357/abb085. hal-02927703

HAL Id: hal-02927703

https://hal.science/hal-02927703v1

Submitted on 7 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The Completed SDSS-IV Extended Baryon Oscillation Spectroscopic Survey: Baryon Acoustic Oscillations with Ly$\alpha$ Forests

Hélion du Mas des Bourboux¹, James Rich², Andreu Font-Ribera³,⁴, Victoria de Sainte Agathe⁵, James Farr³, Thomas Etourneau², Jean-Marc Le Goff², Andrei Cuceu³, Christophe Balland⁵, Julian E. Bautista⁶, Michael Blomqvist⁷, Jonathan Brinkmann⁸, Joel R. Brownstein¹, Solène Chabanier², Edmond Chaussidon⁹, Kyle Dawson¹, Alma X. González-Morales¹⁰,¹¹, Julien Guy⁹, Brad W. Lyke¹², Axel de la Macorba¹³, Eva-Maria Mueller¹⁴, Adam D. Myers¹², Christian Nitschelm¹⁵, Andrea Muñoz Gutiérrez¹³, Nathalie Palanque-Delabrouille², James Parker⁸, Will J. Percival¹⁶,¹⁷,¹⁸, Ignasi Pérez-Ràfols⁵, Patrick Petitjean¹⁹, Matthew M. Pieri⁷, Corentin Ravoux², Graziano Rossi²⁰, Donald P. Schneider²¹,²², Hee-Jong Seo²³, Anže Slosar²⁴, Julianna Stermer⁵, M. Vivek²⁵, Christophe Yèche², and Samantha Youles⁶

¹Department of Physics and Astronomy, University of Utah, 115 S. 1400 E, Salt Lake City, UT 84112, USA; h.du.mas.des.bourboux@gmail.com

²IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France; james.rich@cea.fr

³Department of Physics and Astronomy, University College London, London, UK

⁴Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, E-08193 Bellaterra (Barcelona), Spain

⁵Sorbonne Université, Université Paris Diderot, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies, LPNHE, 4 Place Jussieu, F-75252 Paris, France

⁶Institute of Cosmology & Gravitation, University of Portsmouth, Dennis Sciama Building, Portsmouth, PO1 3FX, UK

⁷Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France

⁸Apache Point Observatory, P.O. Box 59, Sunspot, NM 88349, USA

⁹Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA

¹⁰Departamento de Física, DCI, Campus León, Universidad de Guanajuato, 37150, León, Guanajuato, México

¹¹Consejo Nacional de Ciencia y Tecnología, Av. Insurgentes Sur 1582. Colonia Crédito Constructor, Del. Benito Juárez, C.P. 03940, México D.F. México
¹²Department of Physics and Astronomy, University of Wyoming, Laramie, WY 82071, USA

¹³Instituto de Física, Universidad Nacional Autónoma de México, Ciudad de México, 04510, México

¹⁴Sub-department of Astrophysics, Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK
¹⁵Centro de Astronomía, Universidad de Antofagasta, Avenida Angamos 601, Antofagasta 1270300, Chile

¹⁶Waterloo Centre for Astrophysics, University of Waterloo, 200 University Ave. W, Waterloo, ON N2L 3G1, Canada

¹⁷Department of Physics and Astronomy, University of Waterloo, 200 University Ave. W, Waterloo, ON N2L 3G1, Canada

¹⁸Perimeter Institute for Theoretical Physics, 31 Caroline St. N, Waterloo, ON N2L 2Y5, Canada

¹⁹IAP, Université Paris 6 et CNRS, 98bis blvd. Ara go, F-75014 Paris, France

²⁰Department of Astronomy and Space Science, Sejong University, 209, Neungdong-ro, Gwangjin-gu, Seoul, Republic of Korea
²¹Department of Astronomy and Astrophysics, The Pennsylvania State University, University Park, PA 16802, USA

²²Institute for Gravitation and the Cosmos, The Pennsylvania State University, University Park, PA 16802, USA

²³Department of Physics and Astronomy, Ohio University, Clippinger Labs, Athens, OH 45701, USA

²⁴Brookhaven National Laboratory, Physics Department, Upton, NY 11973, USA

²⁵Indian Institute of Astrophysics, Koramangala, Bangalore 560034, India

Received 2020 May 12; revised 2020 August 11; accepted 2020 August 17; published 2020 October 5

## Abstract

We present a measurement of baryonic acoustic oscillations (BAOs) from Ly$\alpha$ absorption and quasars at an effective redshift $z = 2.33$ using the complete extended Baryonic Oscillation Spectroscopic Survey (eBOSS). The 16th and final eBOSS data release (SDSS DR16) contains all data from eBOSS and its predecessor, the Baryonic Oscillation Spectroscopic Survey (BOSS), providing 210,005 quasars with $z_q > 2.10$ that are used to measure Ly$\alpha$ absorption. We measure the BAO scale both in the autocorrelation of Ly$\alpha$ absorption and in its cross-correlation with 341,468 quasars with redshift $z_q > 1.77$. Apart from the statistical gain from new quasars and deeper observations, the main improvements over previous work come from more accurate modeling of physical and instrumental correlations and the use of new sets of mock data. Combining the BAO measurement from the auto- and cross-correlation yields the constraints of the two ratios $D_H(z = 2.33)/r_d = 8.99 \pm 0.19$ and $D_M(z = 2.33)/r_d = 37.5 \pm 1.1$, where the error bars are statistical. These results are within $1.5\sigma$ of the prediction of the flat-$\Lambda$CDM cosmology of Planck (2016). The analysis code, *picca*, the catalog of the flux transmission field measurements, and the $\Delta\chi^2$ surfaces are publicly available.

*Unified Astronomy Thesaurus concepts: Dark energy (351); Cosmology (343); Large-scale structure of the universe (902); Baryon acoustic oscillations (138); Lyman alpha forest (980)*

## 1. Introduction

One of the surprising characteristics of the widely accepted $\Lambda$CDM cosmological model is that the expansion of the universe is presently accelerating. The most direct evidence for this acceleration comes from the redshift dependence of distances and expansion rates. Luminosity distances to Type Ia supernovae (SNe Ia) provided the first evidence for accelerated expansion (Riess et al. 1998; Perlmutter et al. 1999). More

complete information comes from the baryonic acoustic oscillation (BAO) feature in the matter correlation function, whose position yields both distances and expansion rates normalized to the sound horizon, $r_d$. Constraints on cosmological parameters from BAOs (e.g., Aubourg et al. 2015) are consistent with those from SNe Ia (Scolnic et al. 2018; Jones et al. 2019) and with indirect probes of acceleration like
the spectrum of cosmic microwave background (CMB) anisotropies (Planck Collaboration et al. 2016, 2018).

The BAOs in the pre-recombination universe (Peebles & Yu 1970; Sunyaev & Zeldovich 1970) left their imprints as a peak in the matter correlation function at the sound horizon. The spectrum of CMB anisotropies allows one to set the comoving scale for the BAO peak: $r_d = 147.3 \pm 0.5$ Mpc (Planck Collaboration et al. 2016). This scale is used as a “comoving standard ruler,” i.e., the ruler expands with the expansion of the universe. This BAO scale is arguably simpler than the use of SNe Ia as standard candles since mean SN luminosities may depend on astrophysical conditions.

BAO surveys at a redshift $z$ yield measurements of $D_M(z)/r_d$ and $D_H(z)/r_d$, where $D_M(z) = (1+z)D_A(z)$ is the comoving angular diameter distance to $z$ and $D_H(z) = c/H(z)$ is the Hubble distance corresponding to the expansion rate $H(z)$. The first measurement of BAOs was performed using the autocorrelation determined from galaxy positions (Eisenstein et al. 2005) at $z \sim 0.35$ and the galaxy power spectrum (Cole et al. 2005) at $z \sim 0.1$. At redshifts $z \lesssim 2$ the BAO scale has been studied using discrete tracers such as galaxies (Percival et al. 2007, 2010; Beutler et al. 2011; Blake et al. 2011; Anderson et al. 2012, 2014a, 2014b; Chuang & Wang 2012; Mehta et al. 2012; Padmanabhan et al. 2012; Xu et al. 2013; Ross et al. 2015; Alam et al. 2017; Bautista et al. 2018), galaxy clusters (Hong et al. 2016), and quasars (Ata et al. 2018).

At redshifts $z \gtrsim 2$, the number density of observable discrete tracers is insufficient for high-precision clustering measurements, so BAO studies have been performed using instead opacity fluctuations from the $Ly\alpha$ transition in neutral hydrogen toward background quasars (QSOs). This transition has a rest-frame wavelength of $\lambda_{Ly\alpha} = 121.567$ nm and traces density fluctuations of neutral hydrogen in the intergalactic medium (IGM). From the ground it can be observed along the line of sight to objects with $z \gtrsim 2$, where the transition is redshifted to observed wavelengths $\lambda_{Obs.} \gtrsim 360$ nm. Even though this continuous tracer of the matter density fluctuations has a lower bias than discrete tracers, the statistical power from the large wavelength range where $Ly\alpha$ absorption can be measured allows the observation of BAOs, as suggested by McDonald (2003) and White (2003). The BAO scale was first measured in the $Ly\alpha$ autocorrelation function (Busca et al. 2013; Kirkby et al. 2013; Slosar et al. 2013; Delubac et al. 2015; Bautista et al. 2017; de Sainte Agathe et al. 2019), and then in the $Ly\alpha$–quasar cross-correlation function (Font-Ribera et al. 2014; du Mas des Bourboux et al. 2017; Blomqvist et al. 2019).

Other continuous tracers of the matter density field have been employed but have not yet provided competitive BAO measurements. At redshifts $z \lesssim 2$, Blomqvist et al. (2018) used triply ionized carbon absorption and du Mas des Bourboux et al. (2019) used singly ionized magnesium absorption correlated with the quasar and galaxy distribution to measure large-scale clustering. Both measure the biased 3D correlation of the matter density field and find a signal consistent with that of BAOs. At redshifts $z \gtrsim 2$, Laurent et al. (2016) measured the 3D autocorrelation of quasars and Pérez-Ràfols et al. (2018) the 3D cross-correlation of damped $Ly\alpha$ absorption (DLA) and $Ly\alpha$ systems, but neither of them report a measurement of BAOs.

This study presents the measurements of BAOs from the $Ly\alpha$ autocorrelation function and the $Ly\alpha$–quasar cross-

<figure><img src="image_3.png" /><figcaption>**Figure 1.** High signal-to-noise ratio eBOSS quasar spectrum as a function of observed wavelength. The two spectral regions are $Ly\alpha$ (blue), $104 < \lambda_{RF} < 120$ nm, and $Ly\beta$ (orange), $92 < \lambda_{RF} < 102$ nm. The two solid curves give the two independent best-fit models for $\bar{F}(z)C_q(\lambda_{RF})$, and the dashed curves the unabsorbed continuum $C_q(\lambda_{RF})$ assuming the $\bar{F}(z)$ of Calura et al. (2012). The quasar has a redshift $z_q = 3.058$ and is defined in the catalog DR16Q by Thing_id = 498518806. The $Ly\alpha$ ($\approx 493$ nm) and the overlapping $Ly\beta+O VI$ emission lines ($\approx 420$ nm) are marked with vertical gray dashed lines.</figcaption></figure>

correlation function updated to the 16th data release (DR16; Ahumada et al. 2020) of the fourth generation of the Sloan Digital Sky Survey (SDSS-IV; York et al. 2000; Blanton et al. 2017). This data release contains all of the clustering and $Ly\alpha$ forest data from the completed extended Baryonic Oscillation Spectroscopic Survey (eBOSS; Dawson et al. 2016). We use $Ly\alpha$ absorption in two spectral regions illustrated in the quasar spectrum of Figure 1: the “$Ly\alpha$” region between the quasar $Ly\alpha$ and $Ly\beta-O VI$ emission peaks, and the “$Ly\beta$” region between the $Ly\beta-O VI$ emission peak and the quasar rest-frame Lyman limit. We thus measure two $Ly\alpha$ autocorrelation functions:

$$1. \quad Ly\alpha(Ly\alpha) \times Ly\alpha(Ly\alpha)$$

$$2. \quad Ly\alpha(Ly\alpha) \times Ly\alpha(Ly\beta),$$

where the “transition(spectral region)” notation means, for example, that $Ly\alpha(Ly\beta)$ signifies $Ly\alpha$ absorption in the $Ly\beta$ region. Similarly, we measure two $Ly\alpha$–quasar cross-correlation functions:

$$1. \quad Ly\alpha(Ly\alpha) \times \text{quasar}$$

$$2. \quad Ly\alpha(Ly\beta) \times \text{quasar}.$$

Final BAO constraints are derived using various combinations of these four correlation functions. We do not use $Ly\beta$ absorption because its oscillator strength is $\approx 20\%$ that of $Ly\alpha$, making its contribution to the BAO measurement insignificant. Neither do we use $Ly\alpha(Ly\beta) \times Ly\alpha(Ly\beta)$ since the BAO measurement would be complicated by its superposition with the $Ly\beta(Ly\beta) \times Ly\beta(Ly\beta)$ correlation.

This work builds on previous studies of the $Ly\alpha$ auto- and $Ly\alpha$–quasar cross-correlations from SDSS. Those previous studies are the autocorrelation in SDSS DR12 (BOSS) by Bautista et al. (2017, hereafter B17) and in SDSS DR14 (eBOSS) by de Sainte Agathe et al. (2019, hereafter dSA19) and the cross-correlation in DR12 by du Mas des Bourboux et al. (2017, hereafter dMdB17) and in DR14 by Blomqvist et al. (2019, hereafter B19).

Compared to the studies of dSA19 and B19 on DR14, this analysis using DR16 adds the following improvements:
1. Compared to DR14, DR16 provides 67,606 (25%) more quasars in the redshift range $z_q > 1.77$ as tracers and 24,944 (12%) more in the range $z_q > 2.10$, with spectra allowing measurement of Ly$\alpha$ absorption (see Section 2.1). Compared to DR12, the DR16 sample has 117,940 (52%) more quasars at $z_q > 1.77$ and 44,257 (23%) more quasars at $z_q > 2.10$.

DR16 provides additional observations for greater depth, leading to spectra with improved mean signal-to-noise ratio (S/N) by comparison to the fewer observations of DR12. Furthermore, B17 and dMdB17 used only the best observation of each object, while this study, following B19 and dSA19, co-adds acceptable observations together.

The larger sample of deeper spectra leads to a decrease in the variance of both the auto- and cross-correlation functions by a factor of 0.5 relative to DR12.

2. In addition to the quasar redshift estimators used in DR12 and DR14 analyses, we develop a new estimator referred to as Z_LYAWG. This estimator is similar to the estimator of the standard pipeline (Bolton et al. 2012), but we do not use wavelengths in the vicinity of the Ly$\alpha$ emission line or at shorter wavelengths. The accuracy of the redshifts with this estimator depends less on redshift than that of previously used estimators (see Section 2.2 and Appendix B).

3. Compared to dSA19, the so-called “Ly$\beta$ spectral region” has been extended from the range [97.4,102] nm to [92,102] nm. This extends the spectral region beyond the Ly$\gamma$ absorption line ($\lambda_{Ly\gamma} = 97.2537$ nm) but not as far as the Lyman limit (91.18 nm). Doing so increases the number of Ly$\alpha$ absorption measurements in the Ly$\beta$ spectral region by a factor of two (~1.5 in weighted number; see Section 2.3).

4. We measure any spurious correlations due to the sky subtraction procedure and model their effect on the measured correlation function (see Sections 3.2 and 4). This results in a better fit of the correlation function away from the BAO peak.

5. Two new sets of 3D Gaussian random field simulations have been developed to characterize our combined measurement of the auto- and cross-correlations (see Section 5). These mocks allow us to test the steps of the data analysis, and thus study potential sources of systematic errors and test the estimation of statistical errors. Unlike our previous mock spectra, the Ly$\beta$ spectral region is also simulated and studied.

The analysis code developed by our team, picca: “Package for Igm Cosmological-Correlations Analyses,” is publicly available on GitHub.<fnref n="26" /> Modules include the algorithms to estimate the Ly$\alpha$ forest signal, to compute the auto- and cross-correlations and estimate their covariance, and to fit the measured correlations using an input matter power spectrum. The best-fit BAO results and likelihood, as well as a tutorial, are also given in this repository.<fnref n="27" /> We also publicly release the catalog of estimated Ly$\alpha$ forest signal (Appendix G). The correlation functions, their covariance matrix, and other files used in fits are available upon request.

The publication of this study is coordinated with the release of the final eBOSS measurements of BAOs and redshift-space distortions (RSDs) in the clustering of luminous red galaxies (LRGs) in the redshift range $0.6 < z < 1.0$ (Bautista et al. 2020; Gil-Marín et al. 2020), of emission-line galaxies (ELGs) in the range $0.6 < z < 1.1$ (Raichoor et al. 2020; Tamone et al. 2020; de Mattia et al. 2020), and of quasars in the range $0.8 < z < 2.2$ (Hou et al. 2020; Neveux et al. 2020). The cosmological interpretation of all DR16 eBOSS and BOSS results, in combination with other probes, is found in eBOSS Collaboration et al. (2020).<fnref n="28" />

The analysis procedure we have used is reflected in the organization of this paper and can be summarized as follows. The first step is the selection of quasar and forests samples (Section 2.1), followed by an estimation of quasar redshifts (Section 2.2). The determination of the fluxes in the selected forests is described in Section 2.3. The calculation of the fluctuations in the transmitted flux fraction, $\delta_q(\lambda)$, is presented in Section 2.4.

The calculation of the correlations of $\delta_q(\lambda)$ with itself (autocorrelation) and with quasars (cross-correlation) is described in Section 3. Correlations are a priori a function of the angular and redshift separations and of redshift. It is useful to condense these three dimensions to two dimensions by adopting a fiducial cosmology (Section 3.1). This allows the correlation functions to be calculated simply as functions of radial reparation ($r_\|\|$) and of transverse separation ($r_\perp$).

The continuum fitting procedure uses all flux measurements in a forest to determine each individual $\delta_q(\lambda)$. Thus, individual measured $\delta_q(\lambda)$ have small admixtures of all the true $\delta_q(\lambda)$ in the same forest. This effect is accounted for by the “distortion matrix,” as described in Section 3.5.

In Section 4, we summarize the physical model of the correlation functions. We include the effects of the dominant Ly$\alpha$ absorption, of absorption by high column density (HCD) systems, and of metals. The model of the cross-correlation includes effects of quasar proper motion and of the proximity effect due to quasar UV radiation. Unsuspected correlations not included in the model can be accounted for by the addition of terms that are polynomial functions of ($r_\perp$, $r_\|$).

The validation of the analysis using synthetic data is presented in Section 5. The fits of the correlation functions and the best-fit results for the BAO parameters are presented in Section 6.

We present a comparison with previous Ly$\alpha$ BAO analyses in Section 7, and we summarize our study in Section 8. The constraints on cosmological parameters derived from this analysis, from other DR16 analyses, and from other cosmological probes are presented in a companion paper (eBOSS Collaboration et al. 2020).

## 2. Catalog of Quasar and Ly$\alpha$ Tracers

In our study, we use two different tracers of matter density fluctuations: quasars and Ly$\alpha$ absorption in quasar spectra. We thus have two overlapping quasar samples. The first consists of quasars used as discrete tracers at $z_q \gtrsim 1.77$ and is referred to hereafter as “tracer quasars.” The second sample consists of

<footnote n="26">²⁶ https://github.com/igmhub/picca</footnote>

<footnote n="27">²⁷ https://github.com/igmhub/picca/tree/master/data</footnote>

<footnote n="28">²⁸ A summary of all SDSS BAO and RSD measurements with accompanying legacy figures can be found here: www.sdss.org/science/final-bao-and-rsd-measurements/. The full cosmological interpretation of these measurements can be found here: www.sdss.org/science/cosmology-results-from-eboss/.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Footprint of the eBOSS DR16 survey in a Mollweide projection. The south Galactic cap (SGC) is on the left, and the north Galactic cap (NGC) is on the right. The gray curve shows the position of the Galactic plane. The color scale gives the ratio of the weighted number of pairs for the autocorrelation between eBOSS DR16 and BOSS DR12.</figcaption></figure>

quasars whose spectra are used to measure Ly$\alpha$ absorption, referred to hereafter as “background quasars.” Background quasars have redshifts $z_q > 2.10$ and provide measurements (pixels) of Ly$\alpha$ absorption at $z > 1.96$. This section presents the catalog of both tracers and the details of the pipeline that is used to collect and process the data.

## 2.1. Tracer and Background Quasar Spectra

This analysis benefits from more than 10 years of cosmological observations from SDSS (York et al. 2000) on the 2.5 m Sloan Foundation telescope (Gunn et al. 2006) at the Apache Point Observatory. Most of the tracer quasar spectra and the entirety of the background quasar spectra were gathered during SDSS-III (Eisenstein et al. 2011) in the Baryon Oscillation Spectroscopic Survey (BOSS; Dawson et al. 2013) and in the eBOSS component of SDSS-IV. A small fraction of tracer quasars were observed during SDSS-I and SDSS-II and are found in the seventh data release (Schneider et al. 2010).

The DR16 data used in this paper (Ahumada et al. 2020) contain the complete 5 yr BOSS sample (Fall 2009–Spring 2014) and the 5 yr eBOSS sample (2014 July–2019 March), including its pilot program SEQUELS (Myers et al. 2015), which began in the final year of SDSS-III and concluded in the first year of SDSS-IV. The quasar target selection algorithms for BOSS are summarized in Ross et al. (2012), and those for eBOSS are summarized in Myers et al. (2015). An additional sample of targets was observed as part of the Time-Domain Spectroscopic Survey (TDSS; Morganson et al. 2015; Ruan et al. 2016) and the SPectroscopic IDentification of ERosita Sources survey (SPIDERS; Dwelly et al. 2017). An example high-S/N quasar spectrum is shown in Figure 1. The DR16 footprint is presented in Figure 2, where the color map reflects the statistical improvement between DR12 and DR16 for the autocorrelation of Ly$\alpha$ pixels.

The tracer and background quasars are taken from the DR16 quasar catalog (DR16Q; Lyke et al. 2020), which includes the DR7 quasar catalog (Schneider et al. 2010). Most of the DR7 quasars at $z > 2.1$ were reobserved in BOSS or eBOSS; they are all included in DR16Q. For the few quasars that were not reobserved, we use them as tracer quasars, but not as

background quasars, since their spectra were processed with another pipeline.

Each spectrum from BOSS or eBOSS was acquired using one of the two spectrographs, with wavelength coverage ranging from 360 to 1000 nm (Smee et al. 2013). Each spectrograph views 500 fibers, roughly 450 of which are dedicated to science targets and the remaining to standard F stars for flux calibration and to empty sky locations for sky background subtraction.

The data were processed by version v5_13_0 of the eBOSS pipeline. The reduction is organized in two steps. The pipeline initially extracts the 2D raw data into a 1D flux-calibrated spectrum. During this procedure, the spectra are wavelength and flux calibrated and the individual exposures of one object are co-added into a rebinned spectrum with $\Delta \log(\lambda) = 10^{-4}$. The spectra are then classified as STAR, GALAXY, or QSO, and their redshift is estimated.

DLA systems and broad emission lines (BALs) were identified in the quasar spectra. The details of these searches are presented in Lyke et al. (2020). The BAL search used a procedure similar to that used by Guo & Martini (2019). DLAs were identified using a neural network as described in Parks et al. (2018).

The current DR16 pipeline slightly differs from its last public release (DR14; Abolfathi et al. 2018). Details of the changes can be found in Ahumada et al. (2020). We focus on two relevant changes for the Ly$\alpha$ analysis here. The first concerned improvements in the background estimates on the spectral CCD images. This allowed for more accurate determination of spectral densities $f(\lambda)$. The second improves modeling of the spectra of calibration stars by using a new set of stellar templates. The set was produced for the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration et al. 2016) pipeline and was provided to eBOSS. These templates are able to reduce residuals in flux calibration by improving the modeling of spectral lines in F stars.

For the purposes of measuring Ly$\alpha$ correlations, we rebin three original pipeline $\Delta \log_{10}(\lambda) \sim 10^{-4}$ spectral pixels into one $\Delta \log_{10}(\lambda) \sim 3 \times 10^{-4}$ “analysis pixel.” Throughout this paper, the use of the word “pixel” refers to these rebinned analysis pixels, unless otherwise stated.
<table><caption>**Table 1**
Definitions and Statistics for the Matter Density Tracers</caption><thead><tr><th>Tracer</th><th>λ<sup>min</sup><sub>RF</sub><br>(nm)</th><th>λ<sup>max</sup><sub>RF</sub><br>(nm)</th><th>z<sub>q, min</sub></th><th>z<sub>q, max</sub></th><th>N<sub>q</sub></th><th>N<sub>pix</sub><br>(10<sup>6</sup>)</th><th>S/N</th></tr></thead><tbody><tr><td>Quasar</td><td></td><td></td><td>1.77</td><td>4</td><td>341,468</td><td></td><td></td></tr><tr><td>Lyα(Lyβ)</td><td>92</td><td>102</td><td>2.65</td><td>4</td><td>69,656</td><td>8.4</td><td>1.88</td></tr><tr><td>Lyα(Lyα)</td><td>104</td><td>120</td><td>2.10</td><td>4</td><td>210,005</td><td>34.3</td><td>2.56</td></tr></tbody></table>

Note. The rest-frame wavelength intervals are shown in the second and third columns, the associated redshift interval in the fourth and fifth columns, the number of background or tracer quasars in the sixth column, the total number of spectral pixels in the seventh column, and, in the last column, the mean S/N of the pixels, S/N = (δ + 1) * √w, where the weight, w, is given by the inverse of Equation (4).

## 2.2. Quasar Redshifts

The DR16Q catalog provides up to four estimates of quasar redshift based on the position of the broad emission lines. The absence of visible narrow spectral lines that would be associated with the quasar host galaxies means that individual quasar redshifts have statistical uncertainties of order 100 km s⁻¹ and the different estimators have somewhat different systematic differences. As described in Appendix B, we studied in detail these differences using synthetic quasar spectra. Based on these studies, we choose to employ an estimator that does not use any spectral information in the vicinity of the Lyα emission line or at shorter wavelengths. This choice of rest-frame wavelength coverage mitigates redshift-dependent systematic errors due, for example, to evolution in the mean Lyα opacity.

The redshift range of useful quasars is set by the goal of measuring the quasar–Lyα cross-correlations up to separations of 200 h⁻¹ Mpc. The lowest-redshift Lyα pixel has z = 1.96, meaning that quasars with redshifts z > 1.77 can be used to sample the cross-correlations up to the maximum separation. We impose a maximum redshift of z = 4, beyond which the number of tracers is insufficient for a useful correlation measurement. Our final sample is thus composed of 341,468 tracer quasars, as summarized in Table 1. The redshift distribution is shown in the left panel of Figure 3.

## 2.3. Forest Selection and Residual Calibration

As was done in DR14 (dSA19, B19), we use two different spectral regions shown in Figure 1: the “Lyα region,” $\lambda_{RF} \in [104, 120]$ nm, and the “Lyβ region,” $\lambda_{RF} \in [92, 102]$ nm. The Lyα region is defined between the Lyβ+O VI quasar emission line and the Lyα emission line. Defined as such, we ensure that this region only has Lyα absorption and weak metal absorption (see Pieri et al. 2014). This is the same definition as in previous analyses. The Lyβ+O VI and Lyα emission-line wings are avoided because of their higher quasar-to-quasar diversity.

The Lyβ spectral region covers most of the range between the Lyman limit, $\lambda_{Ly-limit} = 91.18$ nm, and the Lyβ+O VI quasar emission line. This increase in range over the [97.4, 102] nm in DR14 to [92,102] nm increases the number of spectral pixels in the region by a factor of two.

This analysis uses pixels in the observed wavelength range $\lambda \in [360,600]$ nm. The upper limit comes from the requirement $z_q < 4$, which is motivated by the decreasing number of quasars and the increasing contamination from sky emission lines. The lower limit is motivated by high atmospheric

absorption in the UV and the low CCD throughput at the lowest wavelengths. From the previously defined rest-frame wavelength range of the two spectral regions and the observed wavelength range, the minimum background quasar redshift is $z_q > 2$ for the Lyα region and $z_q > 2.53$ for the Lyβ region. This selection gives us 265,328 background quasars for the Lyα region and 103,080 background quasars for the Lyβ region. We removed BALs with BALPROB > 0.9 in the DR16Q Catalog (Lyke et al. 2020). This reduced the numbers to 249,814 and 97,124.

As in dSA19 and B19, all good spectral observations of a given background quasar are stacked to give a higher-S/N spectrum. The use of only good BOSS or eBOSS observations reduces the number of forests to 244,514 and 96,110 for the Lyα and Lyβ regions, respectively.

A spectral region is used only if it contains at least 50 pixels, since short lines of sight may be overfitted, leading to erasing of structure in a correlated manner. This requirement shifts the minimum quasar redshift to $z_q > 2.10$ for the Lyα region and $z_q > 2.65$ for the Lyβ region, leaving 214,542 and 71,602 background quasars, respectively. The fit of the continuum (Equation (2)) fails for a few low-S/N or outlier spectral regions. The final sample has 210,005 and 69,656 background quasars, respectively. We summarize the sample in Table 1.

Spectra and noise variance estimates are corrected for Galactic extinction using the dust map of Schlegel et al. (1998). This had a very small effect in the final results, since extinction corrections are very smooth and can be absorbed in the continuum fitting.

Following B17 and du Mas des Bourboux et al. (2019), we correct for small residual flux calibration errors from the eBOSS pipeline. To do so, we designate a “calibration spectral region” on the red side of the Mg II emission line, $\lambda_{RF} \in [290, 312]$ nm. We compute the mean flux as a function of observed wavelength $f_{calib.}(\lambda_i)$, correcting for the shape of the quasar continuum (see below). The new flux is then defined to be $f(\lambda_i) \rightarrow f(\lambda_i) / f_{calib.}(\lambda_i)$, and its inverse variance $ivar(\lambda_i) \rightarrow ivar(\lambda_i) \times f_{calib.}^2(\lambda_i)$. We find that the correction varies by at most 3%. This process captures residual errors from incomplete modeling of standard F stars and sky emission lines, as well as the Ca II H and K absorption of the Milky Way.

The next step is to mask out spectral intervals, in observed wavelength, where the variance increases sharply owing to unmodeled emission lines from the sky. However, a balance has to be found between pixel quantity and pixel quality. We use the previously defined calibration spectral region to compute the variance as a function of observed wavelength. Comparing to the intrinsic variance of LSS in the Lyα and Lyβ regions, we mask the following six intervals: $\lambda \in [404.30, 405.13]$, $\lambda \in [435.31, 436.51]$, $\lambda \in [545.68, 546.73]$, $\lambda \in [557.05,559.00]$, $\lambda \in [588.33, 590.33]$, and $\lambda \in [629.48, 631.68]$ nm. This mask produces less data loss compared to previous studies, in large part due to the improved sky calibration in the latest eBOSS pipeline. In addition, we mask both Ca II H and K absorption of the Milky Way, $\lambda \in [393.01, 393.82]$ and $\lambda \in [396.55, 397.38]$ nm. Even though we do not observe a significant excess of variance in these regions, masking them avoids spurious large-scale angular correlated absorption between lines of sight.
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Left: redshift distribution of Lyα(Lyα) and Lyα(Lyβ) pixels (numbers divided by 50) and of quasars. Right: normalized weighted redshift distribution of pixel–pixel pairs in the autocorrelation. The pairs are in the BAO region of the correlation functions: $r \in [80, 120] h^{-1} \text{Mpc}$. The orange histogram gives the contribution of correlations involving pixels in the Lyα region, i.e., the Lyα(Lyα) × Lyα(Lyα) correlation function. The green histogram for one pixel is in the Lyβ region, the other in the Lyα region, i.e., the Lyα(Lyα) × Lyα(Lyβ) correlation function. The blue histogram is the combination of the two. The black dotted line shows the effective redshift of the measurements, i.e., the pivot redshift for the BAO measurement (Section 6).</figcaption></figure>

## 2.4. The Flux Transmission Field $\delta_q(\lambda)$

For each spectral region in each line of sight, $q$, the flux transmission field, $\delta_q(\lambda)$, at the observed wavelength, $\lambda$, is obtained from the ratio of the observed flux, $f_q(\lambda)$, to the mean expected flux, $\bar{F}(z)C_q(\lambda)$:

$$ \delta_q(\lambda) = \frac{f_q(\lambda)}{\bar{F}(\lambda)C_q(\lambda)} - 1. \quad (1) $$

Here $C_q(\lambda)$ is the unabsorbed quasar continuum and $\bar{F}(\lambda)$ is the mean transmission. Their product is taken to be a universal function of the rest-frame wavelength, $\bar{C}(\lambda_{\text{RF}})$, corrected by a first-degree polynomial in $\log \lambda$:

$$ \bar{F}(\lambda)C_q(\lambda) = \bar{C}(\lambda_{\text{RF}})(a_q + b_q\Lambda). \quad \Lambda \equiv \log \lambda. \quad (2) $$

The coefficients $(a_q, b_q)$ are fitted separately for the Lyα and Lyβ spectral regions of each line of sight by maximizing the likelihood function

$$ 2 \ln L = -\sum_i \frac{[f_i - \bar{F}C_q(\lambda_i, a_q, b_q)]^2}{\sigma_q^2(\lambda_i)} - \ln[\sigma_q^2(\lambda_i)], \quad (3) $$

where the sum is over all forest pixels for the quasar $q$ and $\sigma_q^2(\lambda)$ is the estimated variance of the flux, $f_q$. Our estimate of the $\sigma_q(\lambda)$ depends on $(a_q, b_q)$, so we include this dependence in the likelihood function of Equation (3).

The variance $\sigma_q^2$ receives a contribution from both instrumental noise (readout and photostatistics) and large-scale structure (LSS), and we model it as the sum of three terms:

$$ \frac{\sigma_q^2(\lambda)}{(\bar{F}C_q(\lambda))^2} = \eta(\lambda)\tilde{\sigma}_{\text{pip},q}^2(\lambda) + \sigma_{\text{LSS}}^2(\lambda) + \frac{\epsilon(\lambda)}{\tilde{\sigma}_{\text{pip},q}^2(\lambda)}. \quad (4) $$

The first term is the instrumental noise, $\tilde{\sigma}_{\text{pip},q}(\lambda) = \sigma_{\text{pip},q}(\lambda)/\bar{F}C_q(\lambda)$, where $\sigma_{\text{pip},q}^2$ is the pipeline estimate of the flux variance. We multiply this by a wavelength-dependent correction, $\eta(\lambda)$, that is found to range from 1.04 at $\lambda = 360$ nm to 1.20 at $\lambda = 580$ nm. The LSS variance, $\sigma_{\text{LSS}}^2$, gives a minimum value of the variance equal to the intrinsic variance of the flux transmission field. Finally, we add an ad hoc factor

proportional to $1/\tilde{\sigma}_{\text{pip},q}^2$. It describes the observed increase of variance, at high S/N, most probably due to quasar-to-quasar spectral diversity. The three terms in Equation (4) dominate for different ranges of flux and wavelength: for our sample of quasars, the instrumental and LSS terms dominate, respectively, for wavelengths less than or greater than $\lambda \approx 450$ nm. The $\epsilon$ term is significant only for the brightest quasars, corresponding to $\approx 10^{-4}$ of the pixels.

The functions $\bar{C}$, $\eta$, $\sigma_{\text{LSS}}^2$, and $\epsilon$ are computed iteratively. One starts with an initial estimate of $\bar{C}(\lambda_{\text{RF}})$ as the mean spectrum using as weights an initial estimate of the $1/\sigma_q^2(\lambda)$. The first estimates of the quasar parameters $a_q$ and $b_q$ are then calculated. The resulting $\delta_q(\lambda)$ are calculated and their variance determined in bins of $\tilde{\sigma}_{\text{pip},q}$ and $\lambda$. The functions $\eta(\lambda)$, $\epsilon(\lambda)$, and $\sigma_{\text{LSS}}(\lambda)$ are then determined by fitting the variance of $\delta_q(\lambda)$ as a function of $\tilde{\sigma}_{\text{pip},q}$ and $\lambda$. The mean spectrum, $\bar{C}(\lambda_{\text{RF}})$, is then recalculated with the new weights. This process is repeated until stable values are obtained after about five iterations. This process is executed independently for the Lyα and Lyβ regions. Figure 4 shows histograms of $\sigma_q^{-2}$ after correction for bias evolution (Equation (7)).

This analysis does not determine separately the mean transmission, $\bar{F}(z)$, and the quasar continuum, $C_q(\lambda)$. In Figure 1, the latter is calculated assuming the $\bar{F}(z)$ measured by Calura et al. (2012) and shown purely for illustration purposes.

The $\delta_q(\lambda)$ for forests with identified DLAs are calculated masking pixels where a DLA reduces the transmission by more than 20%. The masked pixels are not used for the calculation for correlation functions. The absorption in the wings is corrected using a Voigt profile following the procedure of Noterdaeme et al. (2012).

The catalog of $\delta_q(\lambda)$ is publicly available with a data format described in Appendix G. The catalog for data used in this analysis consists of 34.3 and 8.4 million (analysis) pixels in the Lyα and Lyβ regions, respectively. The summary of the Lyα and Lyβ regions is given in Table 1. The redshift distribution of the pixels, assuming Lyα absorption, is presented in the left panel of Figure 3.
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Distribution of weights (Equation (7)) in four wavelength bands. The left and right panels are for the Lyα and Lyβ regions, respectively. For $\lambda < 400$ nm (blue), the distribution has a peak at zero weight (corresponding to noisy spectra) and a second peak at the maximum weight allowed by LSS. With increasing wavelength and the accompanying decreasing noise, the weight distribution becomes more and more concentrated near the maximum weight. Compared to the Lyα forest, the maximum weight of the Lyβ forest is reduced by a factor of $\approx 2$ because of the increased absorption fluctuations, and the wavelength distribution is shifted to lower wavelengths.</figcaption></figure>

As demonstrated in the DR12 Lyα analysis (B17, dMdB17), the quasar fitting of Equation (2) biases the mean and the spectral slope of each individual $\delta_q$ toward zero for each line of sight, resulting in a distortion of the correlation function. This distortion can be modeled if we make the biases exact by redefining the $\delta_q$, per spectral region:

$$ \delta_q(\lambda_i) \rightarrow \sum_j \eta_{ij}^q \delta_q(\lambda_j), \qquad (5) $$

where

$$ \eta_{ij}^q = \delta_{ij}^K - \frac{w_j}{\sum_k w_k} - \frac{w_j(\Lambda_i - \overline{\Lambda}_q)(\Lambda_j - \overline{\Lambda}_q)}{\sum_k w_k(\Lambda_k - \overline{\Lambda}_q)^2}, \qquad (6) $$

and $\overline{\Lambda}_q$ is the mean of $\Lambda = \log \lambda$ for spectrum $q$. In the definition of the projection (Equation (6)), the weights are the ones used in the measurement of the correlation functions:

$$ w_i = \sigma_q^{-2}(\lambda_i) \left( \frac{1 + z_i}{1 + 2.25} \right)^{\gamma_{Ly\alpha} - 1}, \qquad (7) $$

where the redshift evolution of the Lyα bias is taken into account ($\gamma_{Lyα} = 2.9$; McDonald et al. 2006). In Section 3.5, we will use Equations (5) and (6) to establish the relation between the measured correlations and the true correlations.

The mapping applied by Equation (5) changes the evolution of mean $\delta_q$ per wavelength bin slightly, allowing it to deviate from zero. We reintroduce this property for the cross-correlation (Section 3.3) by redefining the $\delta_q$ per observed wavelength bin:

$$ \delta_q(\lambda_i) \rightarrow \delta_q(\lambda_i) - \overline{\delta}(\lambda), \qquad (8) $$

where the weighted average, $\overline{\delta}(\lambda)$, is computed using the same weights as the ones used when computing the correlation function (Equation (7)). This modification ensures that the cross-correlation approaches zero at large scales, whatever the redshift distribution of the quasars.

### 3. Measurement of the Auto- and Cross-correlations

This section presents the measurements of the Lyα auto-correlation functions and of the Lyα-quasar cross-correlation

functions. In Section 2, we presented the catalog of tracer quasars and the two catalogs of Lyα absorption in pixels, covering the Lyα and Lyβ spectral regions. Since these two spectral regions have different sizes, levels of noise, shapes of the quasar continuum, and covariance matrices, we determine the measurements independently. We thus use Lyα absorption in the Lyα region to compute the autocorrelation $\text{Lyα}(\text{Lyα}) \times \text{Lyα}(\text{Lyα})$ and the cross-correlation $\text{Lyα}(\text{Lyα}) \times \text{quasar}$. We use Lyα absorption in the Lyβ region to compute two additional correlation functions $\text{Lyα}(\text{Lyα}) \times \text{Lyα}(\text{Lyβ})$ and $\text{Lyα}(\text{Lyβ}) \times \text{quasar}$.

#### 3.1. From Redshifts and Angles to Distances

In principle, the correlation functions can be computed as a function of redshift and angular separation, $\xi(\Delta z, \Delta \theta)$, as those are directly observed quantities. However, because of the redshift dependence of the comoving angular diameter distance $D_M(z) = (1+z)D_A(z)$ and of the Hubble distance $D_H(z) = c/H(z)$, this would widen the BAO peak unless the data were sorted into multiple redshift bins. To avoid degradation of the BAO feature, we convert angular separations to $r_\perp$ and redshifts separations to $r_\parallel$ by adopting a “fiducial” cosmology. This acts as an optimal data compression designed to maximize the BAO signal if the fiducial cosmology is correct. As demonstrated in Appendix A, the conversion from observed to comoving coordinates does not bias the measurements of the BAO distance scale.

The fiducial cosmology that we adopt is the $\Lambda$CDM cosmology of Planck Collaboration et al. (2016, hereafter Planck 2016). If the fiducial cosmology approximates the true cosmology, then the BAO peak will be at a constant separation, $r_{\text{BAO}} \sim 100 h^{-1} \text{ Mpc}$, for all redshifts. The cosmological parameters for this model are shown in the first part of Table 2 and the derived parameters in the second part; they are computed using the “Code for Anisotropies in the Microwave Background” (CAMB; Lewis et al. 2000). The same assumed cosmology is used to produce and analyze the mock data of Section 5. This cosmology is the same as the one used in DR12 and DR14 studies of Lyα BAOs.

The separation between two tracers is determined along and across the line of sight, $\mathbf{r} = (r_\parallel, r_\perp)$. For tracers $i$ and $j$, of redshift $z_i$ and $z_j$ and offset by an observed angle $\Delta\theta$, the
<table><caption>**Table 2**
Parameters of the Flat-ΛCDM Cosmological Model, from Planck Collaboration et al. (2016), Used for the Production and Analysis of the Mock Data and the Analysis of the Data</caption><thead><tr><th>Parameter</th><th>Planck (2016) Cosmology<br>(TT + lowP)</th></tr></thead><tbody><tr><td>Ω<sub>m</sub>h<sup>2</sup></td><td>0.14252</td></tr><tr><td>=Ω<sub>c</sub>h<sup>2</sup></td><td>0.1197</td></tr><tr><td>+Ω<sub>b</sub>h<sup>2</sup></td><td>0.02222</td></tr><tr><td>+Ω<sub>ν</sub>h<sup>2</sup></td><td>0.0006</td></tr><tr><td>h</td><td>0.6731</td></tr><tr><td>N<sub>ν</sub></td><td>3</td></tr><tr><td>σ<sub>8</sub></td><td>0.8299</td></tr><tr><td>n<sub>s</sub></td><td>0.9655</td></tr><tr><td>Ω<sub>m</sub></td><td>0.31457</td></tr><tr><td>Ω<sub>r</sub></td><td>7.975 10<sup>-5</sup></td></tr><tr><td>r<sub>d</sub> [Mpc]</td><td>147.33</td></tr><tr><td>r<sub>d</sub> [h<sup>-1</sup> Mpc]</td><td>99.17</td></tr><tr><td>D<sub>H</sub>(z = 2.334)/r<sub>d</sub></td><td>8.6011</td></tr><tr><td>D<sub>M</sub>(z = 2.334)/r<sub>d</sub></td><td>39.2035</td></tr><tr><td>f(z = 2.334)</td><td>0.9704</td></tr></tbody></table>

**Note.** The first part of the table gives the cosmological parameters; the second part gives derived quantities used in this paper. They are computed using CAMB (Lewis et al. 2000).

separation is computed as follows:

$$r_{\parallel} = [D_c(z_i) - D_c(z_j)]\cos\left(\frac{\Delta\theta}{2}\right) \quad (9)$$

and

$$r_{\perp} = [D_M(z_i) + D_M(z_j)]\sin\left(\frac{\Delta\theta}{2}\right), \quad (10)$$

where $D_c(z) = \int_0^z dz/H(z)$ is the comoving distance. We will also refer to $(r, \mu)$ in this study; they are defined as $r^2 = r_{\parallel}^2 + r_{\perp}^2$ and $\mu = r_{\parallel}/r$. The quantity $\mu$ is the cosine of the angle formed by the median line of sight of both tracers and the vector $r$.

For quasars, the redshift is defined to be $Z_{\perp}$ (see Section 2.2). We test the effects of other estimators of redshifts on the measurement of the BAO in Appendix D. For spectral pixels, the main absorber, and the one used to measure BAO, is $Ly\alpha$. We thus define the redshift from its rest-frame wavelength. For each spectral pixel $i$, of observed wavelength $\lambda_i$, the redshift is then given by $z_i = \lambda_i/\lambda_{Ly\alpha} - 1 = \lambda_i/121.567 - 1$. The consequences of the presence of other absorption, like Si II(126) or the C IV doublet, are discussed in Section 4.

As given in Table 2, the BAO peak is expected at the separation $r_d \sim 100 h^{-1} \text{ Mpc}$; furthermore, we know from linear theory that the peak has a width (FWHM) of $\approx 20 h^{-1} \text{ Mpc}$. For both these reasons, we compute the different correlations up to $\pm 200 h^{-1} \text{ Mpc}$ along and across the line of sight, i.e., twice the expected BAO scale, with a bin size of $4 h^{-1} \text{ Mpc}$. At the effective redshift of this study, $z_{\text{eff.}} = 2.33$, the BAO scale is $\Delta\theta_{\text{BAO}} \sim 1.5$ deg across the line of sight, and $\Delta z_{\text{BAO}} \sim 0.12$ (corresponding to $\approx 50$ analysis pixels) along the line of sight.

### 3.2. The $Ly\alpha$ Autocorrelation

For the estimator of the 3D autocorrelation of $Ly\alpha$ absorption in spectral pixels, we use

$$\xi_A = \frac{\sum_{(i,j) \in A} w_i w_j \delta_i \delta_j}{\sum_{(i,j) \in A} w_i w_j}. \quad (11)$$

This is the standard “covariance” estimator when the mean has been subtracted: $\langle \delta_i \rangle = 0$, by definition of the projection of Equation (5).

In this equation, $i$ and $j$ refer to two spectral pixels of the flux transmission field, $\delta_i$ and $\delta_j$, from Equation (1). The weights $w_i$ and $w_j$ are defined by Equation (7), where the factor $(1+z)^\gamma$ $\gamma_{Ly\alpha} = 2.9$ is chosen to favor high-redshift pixels where the amplitude of the correlation function is greatest (McDonald et al. 2006). We have tested that small changes in the assumed redshift evolution of the weights do not translate into noticeable differences in the BAO results presented in the next sections.

As explained in Section 3.1, the computation of the correlation is done for all possible pairs of pixels $(i, j)$, of separation $(r_{\parallel}, r_{\perp})$, within $[0, 200] h^{-1} \text{ Mpc}$ in both directions. Each bin $A$ is $4 h^{-1} \text{ Mpc}$ wide in both directions. As a result, each correlation function has $N_{\text{bin}} = 50 \times 50 = 2500$ bins. The sum runs over all possible pairs of pixels from different lines of sight. We exclude pairs of pixels from the same line of sight because of correlated continuum errors that could bias our measurement of the correlation function.

In previous measurements of the $Ly\alpha$ autocorrelation, e.g., B17 and dSA19, pairs of pixels involving the same spectrograph and $r_{\parallel} < 4 h^{-1} \text{ Mpc}$ were avoided to minimize spurious correlations due to the sky subtraction procedure. In this analysis, we keep these pairs, allowing us to model their direct effect on the lowest $r_{\parallel}$ bins and their indirect effect on other bins through distortion due to continuum fitting. This procedure is detailed in Section 4.5.

The resulting correlations have $6.8 \times 10^{11}$ pairs of pixels for the $Ly\alpha(Ly\alpha) \times Ly\alpha(Ly\alpha)$ autocorrelation and $3.8 \times 10^{11}$ for the $Ly\alpha(Ly\alpha) \times Ly\alpha(Ly\beta)$ autocorrelation. The lower number of pixels in the $Ly\beta$ region, combined with the increased absorption fluctuations and Poisson noise, translates to a variance in the $Ly\alpha(Ly\alpha) \times Ly\alpha(Ly\beta)$ correlation function that is approximately three times larger than that in the $Ly\alpha(Ly\alpha) \times Ly\alpha(Ly\alpha)$ function.

The right panel of Figure 3 shows the normalized redshift distribution of the pairs within the BAO region. The $Ly\alpha(Ly\alpha) \times Ly\alpha(Ly\alpha)$ distribution is presented in orange, $Ly\alpha(Ly\alpha) \times Ly\alpha(Ly\beta)$ in green, and the sum of the two in blue. Figure 2 presents the comparison of the weighted number of pairs between this analysis and the DR12 autocorrelation (B17) reproduced using the different improvements of this analysis. The ratio is shown over the BOSS+eBOSS survey footprint, sampled by HEALPix pixels (Górski et al. 2005).

Figure 5 presents the $Ly\alpha(Ly\alpha) \times Ly\alpha(Ly\alpha)$ autocorrelation for the data on the left and for the best-fit model on the right (Section 6). The correlation is multiplied by the separation $|r|$ for visual purposes, and the color bar is saturated and symmetric around zero. Even though each individual bin of the correlation is noisy when presented in 2D, the BAO scale is seen in the data at large $\mu$ ($r \approx r_{\parallel} \approx 100 h^{-1} \text{ Mpc}$) by the transition from blue (negative values) to white (zeros).

The measured autocorrelation, $Ly\alpha(Ly\alpha) \times Ly\alpha(Ly\alpha)$, is also presented in the top four panels of Figure 6. These panels
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Measured (left) and best-fit model (right) Lyα autocorrelation function for two pixels in the Lyα region: Lyα(Lyα) × Lyα(Lyα). The correlation is multiplied by the separation |r|, and the color bar is saturated and symmetric around zero for visualization purposes. The BAO can be observed as a quarter of a ring at r ~ 100 h⁻¹ Mpc.</figcaption></figure>

show the 2D correlation of Figure 5 reduced to a weighted 1D correlation for four different wedges of |μ| = |r_‖/r|. In the same figure, the best-fit model is shown in red and is discussed in Section 6. The BAO scale peak at r ~ 100 h⁻¹ Mpc is visible, especially for μ > 0.8. Four similar panels in Appendix F present the autocorrelation Lyα(Lyα) × Lyα(Lyβ) (Figure F1).

The autocorrelation, Lyα(Lyα) × Lyα(Lyα), is also measured in two redshift bins. This allows us to look at the evolution of the different parameters of the best-fit models, e.g., the bias of Lyα, along with allowing systematic tests. In their DR14 study, dSA19 showed that for the Lyα autocorrelation the optimal way to assign pixel pairs to the high- or low-redshift measurement is to cut on the mean of the maximum redshift of the two forests, rather than simply on the redshift of the two pixels. This minimizes the cross-covariance between the low-z and high-z autocorrelations (Section 3.4). As in dSA19, we chose zcut = 2.5, in order to approximately get the same weighted number of pairs in both redshift splits. The resulting autocorrelation function is presented in the four top panels of Figure F2 in Appendix F.

The covariance matrix of each autocorrelation function has $N_{\text{bin}}^2 = 2500 \times 2500 = 6,250,000$ elements. For two bins A and B of the correlation function $\xi$, the covariance is defined by

$$C_{AB} = \langle \xi_A \xi_B \rangle - \langle \xi_A \rangle \langle \xi_B \rangle. \quad (12)$$

In this study, the covariance matrix is estimated by dividing the sky into subsamples defined by HEALPix pixels and computing their weighted covariance. Using nside = 16 over the eBOSS footprint, we get around 880 subsamples, each covering 3.7 × 3.7 = 13.4 deg² on the sky. This solid angle is equivalent to a 250 × 250 (h⁻¹ Mpc)² patch at zeff = 2.33. Lyα pixel pairs are assigned to HEALPix pixels according to the Lyα pixel with the smallest right ascension. The covariance is then given by the following, neglecting the small correlations between subsamples:

$$C_{AB} = \frac{1}{W_A W_B} \sum_s W_A^s W_B^s [\xi_A^s \xi_B^s - \xi_A \xi_B]. \quad (13)$$

Here s is a subsample with summed weight $W_A^s$ and measured correlation $\xi^s$, and $W_A = \sum_s W_A^s$.

The covariance is dominated by the diagonal elements that are of order $C_{AA} \approx \langle \delta^2 \rangle^2/N_A$, where $\langle \delta^2 \rangle$ is the pixel variance and $N_A$ is the number of pixel pairs in the bin A. Deviations from this simple expression are due to intraforest correlations

that make the effective number of independent pairs less than $N_A$. We find

$$\text{Var}_A \approx \frac{\langle \delta^2 \rangle^2}{f N_A^{\text{pair}}}, \quad (14)$$

where $\langle \delta^2 \rangle \approx 0.13$ (0.24) and $f \approx 0.4$ (0.8) for the Lyα and Lyβ regions. For the Lyα(Lyα) × Lyα(Lyα) correlation this is equivalent to $\text{Var}_A \approx 1.5 \times 10^{-10} (100 h^{-1} \text{ Mpc}/r_\perp)$, where the $r_\perp$ dependence reflects the approximate proportionality between $N_A$ and $r_\perp$.

Off-diagonal elements of the covariance matrix are due to intraforest correlations that break the independence of pixel pairs in different bins, $A \neq B$.<fnref n="29" /> Intraforest correlations reflect both physical correlations of absorption and the effect of continuum fitting. Imperfections in the continuum fitting increase both $\langle \delta^2 \rangle$ and the intraforest correlations and are therefore reflected in the covariance matrix, as long as the imperfections themselves are not correlated between different forests.

The size of the off-diagonal elements is made clearer through the normalized covariance matrix, i.e., the correlation matrix, with elements in [-1, 1]:

$$\text{Corr}_{AB} = \frac{C_{AB}}{\sqrt{\text{Var}_A \text{Var}_B}}, \quad (15)$$

where $\text{Var}_A = C_{AA}$ is the variance. The largest values are $\text{Corr}_{AB} \approx 0.4$ for $r_\perp^A = r_\perp^B$ and $|r_\perp^A - r_\perp^B| = 4h^{-1} \text{ Mpc}$. Elements with $r_\perp^A \neq r_\perp^B$ are very small, <0.03.

The estimates of the off-diagonal elements of the correlation matrix from Equation (15) are noisy, and we smooth them by modeling them as a function of the difference of separation along and across the line of sight: $\text{Corr}_{AB} = \text{Corr}(r_\perp^A, r_\perp^B, r_\perp^A, r_\perp^B) = \text{Corr}(\Delta r_\perp, \Delta r_\perp)$, with $\Delta r_\perp = |r_\perp^A - r_\perp^B|$ and $\Delta r_\perp = |r_\perp^A - r_\perp^B|$. They are presented in the left panel of Figure 7, where they are shown to decrease rapidly as a function of $\Delta r_\perp$ and $\Delta r_\perp$.

We have also calculated the covariance matrix using other techniques presented in Appendix C. We find no significant change from the estimates of the covariance matrix that is used in the main analysis.

### 3.3. The Lyα–Quasar Cross-correlation

For the Lyα–quasar cross-correlation we use the same estimator as the one from previous studies (Font-Ribera et al. 2012b, 2013, dMdB17, B19). It is defined to be the weighted mean of the flux transmission field at a given separation from a quasar:

$$\xi_A = \frac{\sum_{(i,j) \in A} w_i w_j \delta_i}{\sum_{(i,j) \in A} w_i w_j}. \quad (16)$$

In this equation, i indexes a flux pixel and j a quasar. The weights $w_i$ are as defined in Equation (7) for the Lyα absorption fluctuations. In their study, du Mas des Bourboux et al. (2019) fitted the redshift evolution of the quasar bias and found a best-fit power law of index $\gamma_{\text{quasar}} = 1.44 \pm 0.08$. The

<footnote n="29">²⁹ The role of intraforest correlations in the covariance matrix is made explicit in the “Wick” calculation of the covariance in Appendix C.</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 6. Lyα autocorrelation function (top four panels) and Lyα-quasar cross-correlation (bottom four panels), for pixels in the Lyα region: Lyα(Lyα) × Lyα(Lyα) and Lyα(Lyα) × quasar. The correlations are multiplied by the separation r² to better see the BAO scale. The black points give the measured correlation of Figures 5 and 8, and the red curves give the best-fit models, in four wedges of |μ| = |r|/r. The dashed red lines give the best-fit models, interpolated to lower and higher separation than the fitted range: r ∈ [10, 180] h⁻¹ Mpc. Along the line of sight, |μ| ∈ [0.95, 1], the contribution from metals, Si II(119), Si II(119.3), and Si III(120.7), can be observed as extra knees or peaks at, respectively, r ~ 20 and 60 h⁻¹ Mpc.</figcaption></figure>

quasar weights are then defined to be

$$w_j = \left( \frac{1 + z_j}{1 + 2.25} \right)^{\gamma_{\text{quasar}} - 1} . \qquad (17)$$

In a similar way as for the autocorrelation, the cross-correlation is computed for all pixel-quasar pairs, though omitting pairings of pixels with their own background quasar whose mean correlation vanishes owing to the continuum fitting procedure. The correlation is computed for all pairs within r⊥ ∈ [0, 200] h⁻¹ Mpc. Unlike the autocorrelation, the cross-correlation is not symmetric by permutation of the two tracers. We thus have the opportunity to

define positive values of the separation along the line of sight, r‖, when the tracer quasar is in front of the Lyα pixel tracer, i.e., zLyα > zquasar. The line-of-sight separation then ranges over r‖ ∈ [−200, 200] h⁻¹ Mpc. With a width of 4 h⁻¹ Mpc, the correlation is computed on Nbin = 100 × 50 = 5000 bins.

The weighted distribution of the pair redshifts in the BAO region is similar for the cross-correlations to what is presented for the autocorrelation in the right panel of Figure 3. The comparison on the sky of the weighted number of pairs between DR12 (dMdB17) and this study is also similar to what is shown for the autocorrelation in Figure 2.
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Normalized covariance matrix, i.e., the correlation matrix (Equation (15)), as a function of $(\Delta r_{\perp}, \Delta r_{\parallel})$ for the autocorrelation (left) and cross-correlation (right) as estimated by subsampling (Equation (13)). The solid lines are for the $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{Ly}\alpha(\text{Ly}\alpha)$ and $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{quasar}$ functions, and the dashed lines for the $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{Ly}\alpha(\text{Ly}\beta)$ and $\text{Ly}\alpha(\text{Ly}\beta) \times \text{quasar}$ functions. The correlations have been smoothed as explained in the text.</figcaption></figure>

Figure 8 presents in 2D the measured cross-correlation, $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{quasar}$, and its best-fit model (Section 6). In such a display, the BAO scale would be seen as a half ring of radius $r \sim 100 h^{-1}$ Mpc. However, the signal is difficult to see owing to the noise in the data. We show in the bottom four panels of Figure 6 the same cross-correlation but averaged over four different wedges of $|\mu|$. There the BAO scale can be clearly observed as a dip at $r \sim 100 h^{-1}$ Mpc. We note that the correlation is reversed, i.e., negative, with respect to the matter correlation function, because the bias of $\text{Ly}\alpha$ absorption is negative and the bias of quasars is positive, giving a negative product of biases. The cross-correlation $\text{Ly}\alpha(\text{Ly}\beta) \times \text{quasar}$ is presented in Figure F1 of Appendix F.

As in B19, we alternatively compute the cross-correlation, $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{quasar}$, in two redshift bins. For this measurement, it is important that $\langle\delta\rangle = 0$ per bin of observed wavelength. To ensure this, we recompute the $\delta$ of Section 2.3, independently for the two redshift bins. We split lines of sight by their background quasar redshift, $z_q$, and select the splitting redshift, $z_{\text{cut}} = 2.57$, to have approximately the same weighted number of pairs in each binned cross-correlation. This definition of redshift bins allows us to minimize the cross-covariance between the two samples (Section 3.4). Figure F2 of Appendix F presents, in its four bottom panels, the two redshift bins of the cross-correlation, in four wedges of $\mu$.

The covariance matrix is computed using the same estimators as in Section 3.2 for the autocorrelation. The covariance matrix is dominated by its diagonal elements, i.e., the variance. The latter is approximately inversely proportional to the number of pairs:

$$ \text{Var}_A \approx \frac{\langle\delta^2\rangle}{0.7 N_A^{\text{pair}}}, \qquad (18) $$

where the factor 0.7 gives the effective loss of number of pairs, due to correlations between neighboring pixels. For the DR16 data set, the resulting variance for the $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{quasar}$ correlation is $\text{Var}_A \approx 8.6 \times 10^{-8}(100 h^{-1} \text{Mpc}/r_{\perp})$. There are a total of $\approx 1.2 \times 10^9$ quasar–pixel pairs in the BAO region, $80 < r < 120 h^{-1} \text{Mpc}$.

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Measured (left) and best-fit model (right) $\text{Ly}\alpha$–quasar cross-correlation function where the pixel is in the $\text{Ly}\alpha$ region: $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{quasar}$. The correlation is multiplied by the separation $|r|$, and the color bar is saturated and symmetric around zero for visualization purposes. The BAO can be observed as half a ring at $r \sim 100 h^{-1} \text{Mpc}$.</figcaption></figure>

The off-diagonal elements of the covariance matrix are presented in the right panel of Figure 7. As with those for the autocorrelation, they decrease rapidly as a function of $\Delta r_{\parallel}$ and $\Delta r_{\perp}$.

### 3.4. The Cross-covariance between the Correlations

We modify the estimator of the covariance matrix via subsampling defined by Equation (13) to compute the cross-covariance between the different measured correlation functions:

$$ C_{AB}^{12} = \frac{1}{W_A W_B} \sum_s W_A^s W_B^s [\xi_A^{1,s} \xi_B^{2,s} - \xi_A^1 \xi_B^2]. \qquad (19) $$

In this equation, $\xi^1$ and $\xi^2$ are two different measured correlation functions, e.g., $\xi^1$ could be the autocorrelation
Lyα(Lyα) × Lyα(Lyα) and ξ² the cross-correlation Lyα(Lyα) × quasar.

The different correlation functions are found to be marginally correlated. For example, the correlations between the Lyα (Lyα) × Lyα(Lyα) and Lyα(Lyα) × quasar functions are less than 1%, as shown in Figure 11.

### 3.5. The Distortion Matrices

The continuum fitting (Equation (2)) and the projection (Equation (5)) mix the δq within a given forest and thereby modify significantly the correlation function. The transformation (Equation (8)) also has a small effect on the correlations. Under the assumption that the δ transformations are linear, the true and distorted correlation functions are related by a “distortion matrix,” $D_{AB}$:

$$ \hat{\xi}_{\text{distorted}}(A) = \sum_B D_{AB} \xi_{\text{true}}(B). \quad (20) $$

To the extent that the projection given by Equations (5) and (6) captures the full distortion, the distortion matrices for the auto- and cross-correlations can be read off from the coefficients in Equation (6):

$$ D_{AB}^{\text{auto}} = W_A^{-1} \sum_{ij \in A} w_i w_j \left( \sum_{i'j' \in B} \eta_{ii'} \eta_{jj'} \right), \quad (21) $$

$$ D_{AA'}^{\text{cross}} = W_A^{-1} \sum_{(i,j) \in A} w_i w_j \left( \sum_{(i',j') \in A'} \eta_{ii'} \right). \quad (22) $$

In the fits of the data described in Section 6, the physical model of the correlation function (Section 4) is multiplied by the distortion matrix before comparison with the data. This procedure is validated with the analysis of the mock data sets (Section 5), where it is shown that the BAO parameters are accurately recovered in the presence of distortion.

### 4. Model of the Correlations

This section presents the theoretical model of the auto- and cross-correlations. Even though this analysis focuses on measuring the BAO scale along and across the line of sight, much effort has also been invested into better understanding the global shape of the correlation functions. The improved modeling of the global shape allows us to better understand and test possible sources of systematic errors and to better measure BAOs independently of the overall correlation function. Apart from the contribution of the sky calibration to the autocorrelation (Section 4.5), the model is the same as that used in the DR12 and DR14 analyses (B17; dMdB17; dSA19; B19).

A cosmological analysis of Lyα forest correlations is possible only because the dominant Lyα absorption has two essential characteristics: it traces the underlying matter fluctuations and, for a fixed set of cosmological parameters, there is a unique mapping between wavelength and distance. Absorption by metals also traces matter fluctuations but has a different wavelength-distance relation. Correlations due to instrumental and analysis effects do not relate at all to cosmology. Fortunately, the BAO feature due to the Lyα correlations is not degenerate with any of these secondary effects.

For the dominant Lyα–Lyα or Lyα–quasar correlations, the BAO peak in the space of angular and redshift separation appears at $\Delta\theta \sim r_d/D_H(z)$ and $\Delta z \sim r_d/D_M(z)$, where $D_H$ and

$D_M$ are the Hubble and comoving angular diameter distances calculated assuming Lyα absorption. We want to measure $r_d/D_H(z)$ and $r_d/D_M(z)$ in a way that does not depend significantly on the smooth part of the correlation function underneath the BAO peak. To do this, we follow the procedure described in Section 4.1 of separating the correlation function into two components:

$$ \xi(r_{\parallel}, r_{\perp}, \alpha_{\parallel}, \alpha_{\perp}) = \xi_{\text{sm}}(r_{\parallel}, r_{\perp}) + \xi_{\text{peak}}(\alpha_{\parallel}r_{\parallel}, \alpha_{\perp}r_{\perp}), \quad (23) $$

where $\xi_{\text{sm}}$ is the smooth correlation function, without the BAO peak, and $\xi_{\text{peak}}$ is the peak-only correlation. The BAO parameters in this equation are

$$ \alpha_{\parallel} = \frac{[D_H(z_{\text{eff}})/r_d]}{[D_H(z_{\text{eff}})/r_d]_{\text{fid}}} \quad \text{and} \quad \alpha_{\perp} = \frac{[D_M(z_{\text{eff}})/r_d]}{[D_M(z_{\text{eff}})/r_d]_{\text{fid}}}. \quad (24) $$

In the standard fits, we assume that the correlation function for $\alpha_{\perp} = \alpha_{\parallel} = 1$ is that of the fiducial cosmology of Table 2 as calculated by CAMB.<fnref n="30" /> The smooth and the peak parts of the correlation are obtained at the level of the matter power spectrum, as described below (Section 4.1).

Though the BAO parameters ($\alpha_{\parallel}$, $\alpha_{\perp}$) depend on the assumed cosmology, the measured ($D_H(z_{\text{eff}})/r_d$, $D_M(z_{\text{eff}})/r_d$) do not. This was studied in detail by Carter et al. (2020) in the context of galaxy correlations. We further illustrate this property in Appendix A, where we use a different fiducial cosmology.

One feature of the peak-smooth splitting (Equation (23)) is that unless the fit yields $\alpha_{\perp} = \alpha_{\parallel} = 1$, the best-fit correlation function does not obviously correspond to any physical cosmological model. In Appendix A we choose a physical model with values of ($h, \Omega_{k}h^{2}$) that yield $\alpha_{\perp} = \alpha_{\parallel} = 1$, thereby yielding a best fit that is a physical correlation function. As demonstrated in Appendix A, changing the model of the smooth continuum through a change in cosmological parameters produces a difference in the derived values of ($D_H(z_{\text{eff}})/r_d$, $D_M(z_{\text{eff}})/r_d$) of less than 1 part in 300, negligible compared to the statistical precision.

In addition to the primary correlations from which we measure the BAO peak, the full correlation function receives subdominant contributions from other effects that we describe here. The theoretical model of the correlation function, $\xi^t$, of Equation (25) is composed of different correlations. For the autocorrelation of Lyα, it is given by

$$ \xi^t = \xi^{\text{Ly}\alpha \times \text{Ly}\alpha} + \sum_m \xi^{\text{Ly}\alpha \times m} + \sum_{m_1, m_2} \xi^{m_1 \times m_2} + \xi^{\text{sky}}, \quad (25) $$

where $\xi^{\text{Ly}\alpha \times \text{Ly}\alpha}$ (Section 4.2) is the model of the autocorrelation of Lyα absorption, and $\xi^{\text{Ly}\alpha \times m}$ and $\xi^{m_1 \times m_2}$ (Section 4.3) give the contribution of other absorbers in the Lyα and Lyβ regions. The correlation $\xi^{\text{sky}}$ (Section 4.5) models the correlations due to the sky subtraction procedure of the eBOSS pipeline.

In a similar way, for the measured Lyα–quasar cross-correlation the theoretical model is given by

$$ \xi^t = \xi^{\text{Ly}\alpha \times \text{QSO}} + \sum_m \xi^{\text{QSO} \times m} + \xi^{\text{TP}}. \quad (26) $$

In this equation, the first term gives the cross-correlation between Lyα and quasars, and the second gives the

<footnote n="30">³⁰ https://github.com/cmbant/CAMB</footnote>
contribution of the cross-correlation between quasars and other absorbers in the Lyα and Lyβ spectral regions (Section 4.3). Finally, the third term models the “transverse proximity” (TP) effect of quasar radiation on the surrounding gas (Section 4.4).

## 4.1. The Power Spectra

The components of the correlation function that reflect the underlying matter correlations are given by the Fourier transform of the tracer biased power spectrum:

$$ \hat{P}(\mathbf{k}) = b_i b_j (1 + \beta_i \mu_k^2)(1 + \beta_j \mu_k^2) P_{\text{QL}}(\mathbf{k}) F_{\text{NL}}(\mathbf{k}) G(\mathbf{k}), \quad (27) $$

where the vector $\mathbf{k} = (k_\parallel, k_\perp) = (k, \mu_k)$ of modulus $k$ has components along and across the line of sight, $(k_\parallel, k_\perp)$, with $\mu_k = k_\parallel/k$. The bias and RSD parameters, $(b, \beta)$, are for the tracer $i$ or $j$. $P_{\text{QL}}$ is the quasi-linear power spectrum defined below, $F_{\text{NL}}$ corrects for nonlinear effects at large $\mathbf{k}$, and $G(\mathbf{k})$ is a damping term that accounts for averaging of the correlation function in individual $(r_\parallel, r_\perp)$ bins. For the autocorrelation, the main tracer is the Lyα absorption, $i = j$, and for the cross-correlation two main tracers play a role: Lyα absorption and quasars, $i \neq j$.

For the cross-correlation power spectrum we do not include relativistic effects that lead to an asymmetry in $r_\parallel$ (quasar more distant or less distant than the forest). These effects have been calculated to be sufficiently small to be negligible for this study (Lepori et al. 2020). They were included in the study of B19 and were shown to be partially degenerate with the parameter $(\Delta r_\parallel, q_{so})$ describing systematic errors in quasar redshift.

The quasi-linear power spectrum provides for the aforementioned decoupling of the peak component:

$$ P_{\text{QL}}(\mathbf{k}, z) = P_{\text{sm}}(k, z) + \exp\left[-\frac{k_\parallel^2 \Sigma_\parallel^2 + k_\perp^2 \Sigma_\perp^2}{2}\right] P_{\text{peak}}(k, z). \quad (28) $$

The smooth component, $P_{\text{sm}}$, is derived from the linear power spectrum, $P_L(k, z)$, via the sideband technique (Kirkby et al. 2013), which is implemented in picca, with the help of nbodykit<fnref n="31" /> (Hand et al. 2018). The CAMB $P_L(z_{\text{eff}})$ is Fourier transformed into the correlation function where two sidebands are defined on both side of the BAO. A smooth function is then fitted to connect the two sidebands, allowing us to produce a correlation function, without the BAO peak. The latter correlation is Fourier transformed back to a smooth power spectrum, $P_{\text{sm}}$, that does not have the BAO wiggle features. The peak power spectrum is thus defined as $P_{\text{peak}} = P_L - P_{\text{sm}}$.

The correction for nonlinear broadening of the BAO peak (Eisenstein et al. 2007) is parameterized by $(\Sigma_\parallel, \Sigma_\perp)$, with $\Sigma_\perp = 3.26 h^{-1} \text{ Mpc}$ and

$$ \frac{\Sigma_\parallel}{\Sigma_\perp} = 1 + f, \quad (29) $$

where $f = d(\ln g)/d(\ln a) \approx \Omega_m^{0.55}(z)$ is the linear growth rate of structure, resulting in $\Sigma_\parallel = 6.42 h^{-1} \text{ Mpc}$.

The function $F_{\text{NL}}(\mathbf{k})$ in Equation (27) accounts for nonlinear effects at small scales. For the autocorrelation, the most

important effects are thermal broadening, peculiar velocities, and nonlinear structure growth. We use Equation (3.6) of Arinyo-i-Prats et al. (2015) with parameter values from their Table 7 interpolated to our effective redshift $z = 2.334$. For the cross-correlation, the most important effect is that of quasar nonlinear velocities, and, following Percival & White (2009), we adopt

$$ F_{\text{NL}}^{\text{cross}}(k_\parallel) = \frac{1}{1 + (k_\parallel \sigma_v)^2}, \quad (30) $$

where $\sigma_v$ is a free parameter. This function also takes into account statistical quasar redshift errors.

The final term in Equation (27) accounts for the effect of the binning of the correlation function on the separation grid. We assume the distribution to be homogeneous on each bin<fnref n="32" /> and compute the function $G(\mathbf{k})$ as the product of the Fourier transforms of the rectangle functions that model a uniform square bin:

$$ G(\mathbf{k}) = \text{sinc}\left(\frac{k_\parallel R_\parallel}{2}\right) \text{sinc}\left(\frac{k_\perp R_\perp}{2}\right), \quad (31) $$

where $R_\parallel$ and $R_\perp$ are the radial and transverse widths of the bins, respectively.

## 4.2. Lyα and Quasar Bias Parameters

The dominant Lyα absorption can be viewed as the sum of contributions from the diffuse IGM and from HCD<fnref n="33" /> systems. 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 Section 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_\parallel$ dependence of the effective bias (Font-Ribera & Miralda-Escudé 2012):

$$ \left\{ \begin{aligned} b'_{\text{Ly}\alpha} &= b_{\text{Ly}\alpha} + b_{\text{HCD}} F_{\text{HCD}}(k_\parallel) \\ b'_{\text{Ly}\alpha} \beta'_{\text{Ly}\alpha} &= b_{\text{Ly}\alpha} \beta_{\text{Ly}\alpha} + b_{\text{HCD}} \beta_{\text{HCD}} F_{\text{HCD}}(k_\parallel), \end{aligned} \right. \quad (32) $$

where $(b_{\text{Ly}\alpha}, \beta_{\text{Ly}\alpha})$ and $(b_{\text{HCD}}, \beta_{\text{HCD}})$ are the bias parameters associated with the IGM and HCD systems, respectively, and $F_{\text{HCD}}$ is a function that depends on the number and column density distribution of HCDs. Rogers et al. (2018) numerically calculated $F_{\text{HCD}}$ using hydrodynamical simulations, and we find that a simple exponential form, $F_{\text{HCD}} = \exp(-L_{\text{HCD}}k_\parallel)$, well approximates their results. Here $L_{\text{HCD}}$ is a typical length scale for unmasked HCDs. For eBOSS spectral resolution, DLA identification is possible for DLA widths (wavelength interval for absorption greater than 20%) greater than ~2.0 nm, corresponding to ~14 $h^{-1}$ Mpc in our sample. Degeneracies between $L_{\text{HCD}}$ and other anisotropic parameters make the fitting

<footnote n="32">³² In fact, in the perpendicular direction the distribution is approximately proportional to $r_\perp$; however, assuming homogeneity produces a sufficiently accurate correlation function (B17).</footnote>

<footnote n="33">³³ In this and other eBOSS publications we use the term “HCD systems” to describe systems with a neutral hydrogen column density above $10^{17.2} \text{cm}^{-2}$, i.e., including both Lyman limit systems and DLAs.</footnote>

<footnote n="31">³¹ https://github.com/bccp/nbodykit</footnote>
somewhat unstable. We therefore impose $L_{HCD} = 10 h^{-1}$ Mpc while fitting for the bias parameters $b_{HCD}$ and $\beta_{HCD}$. 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. Cuceu et al. (2020) also showed that fixing $L_{HCD}$ has a minimal impact on BAO results when doing a Bayesian analysis of Ly$\alpha$ BAOs.

In our measured correlations functions, different ($r_\parallel$, $r_\perp$) bins have mean redshifts that vary over the range $2.32 < \bar{z} < 2.39$. Therefore, in order to fit a unique function to the data, we need to assume a redshift dependence of the bias parameters. Following McDonald et al. (2006), we assume that the product of $b_{Ly\alpha}$ and the growth factor of structures varies with redshift as $(1+z)^{\gamma_\alpha-1}$, with $\gamma_\alpha = 2.9$, while we make use of the approximation that $\beta_{Ly\alpha}$ does not depend on redshift. Because the fit of the cross-correlation is only sensitive to the product of the quasar and Ly$\alpha$ biases, we choose to adopt a value for the quasar bias. Following the analysis of du Mas des Bourboux et al. (2019) based on a compilation of quasar bias measurements (Croom et al. 2005; Shen et al. 2013; Laurent et al. 2016, 2017), we set $b_q \equiv b_q(z_{eff}) = 3.77$ with a redshift dependence given by

$$b_q(z) = 3.60 \left( \frac{1+z}{3.334} \right)^{1.44} \quad (33)$$

The quasar RSD, assumed to be redshift independent in the matter-dominated epochs explored here, is

$$\beta_q = \frac{f}{b_q}. \quad (34)$$

Setting $f = 0.9704$ for our fiducial cosmology yields $\beta_q = 0.269$, which we take to be redshift independent.

### 4.3. Absorption by Metals

We model the power spectrum for correlations involving metals with the same form we use for Ly$\alpha$–Ly$\alpha$ and Ly$\alpha$–quasar correlations (Equation (27)) except that HCD effects are neglected. Each metal species, $n$, has bias parameters ($b_n$, $\beta_n$). The Fourier transform of the power spectrum $P_{mn}(k, z)$ for the absorber pair ($m, n$) is then the model correlation function of the pair ($m, n$): $\xi_{\text{mod}}^{m-n}(\tilde{r}_{\parallel}, \tilde{r}_{\perp})$, where $(\tilde{r}_{\parallel}, \tilde{r}_{\perp})$ are the separations calculated using the correct rest-frame wavelengths, ($\lambda_m$, $\lambda_n$). For the cross-correlation, the same form is used except that the bias parameters for the metal $n$ are replaced with the quasar bias parameters and the separations are calculated using the quasar redshift.

Since we assign redshifts to pixels assuming Ly$\alpha$ absorption, the rest-frame wavelength we ascribe to a metal transition is not equal to the true rest-frame wavelength. For $m \neq n$ and $n = \text{Ly}\alpha$, this misidentification results in a shift of the model correlation function with vanishing separation of absorbers corresponding to reconstructed separations $r_\perp = 0$ and $r_\parallel \approx (1+z)D_H(z)(\lambda_m - \lambda_{Ly\alpha})/\lambda_{Ly\alpha}$. The values of $r_\parallel$ for the major metal transition are given in Table 3. For $m = n$, the reconstructed separations are scaled from the true separations by a factor $D_H(z)/D_H(z_\alpha)$ for $r_\parallel$ and by a factor $D_M(z)/D_M(z_\alpha)$ for $r_\perp$, where $z$ and $z_\alpha$ are the true and reconstructed redshifts, respectively.

<table><caption>**Table 3**
Metal Lines Included in the Synthetic Spectra</caption><thead><tr><th>Metal Line</th><th>λ<sub>m</sub> (nm)</th><th>Relative Strength (×10<sup>3</sup>)</th><th>r<sub>||</sub> (h<sup>-1</sup> Mpc)</th></tr></thead><tbody><tr><td>Si III</td><td>120.7</td><td>1.892</td><td>-21</td></tr><tr><td>Si IIa</td><td>119.0</td><td>0.642</td><td>-64</td></tr><tr><td>Si IIb</td><td>119.3</td><td>0.908</td><td>-56</td></tr><tr><td>Si IIc</td><td>126.0</td><td>0.354</td><td>+111</td></tr></tbody></table>

**Note.** For each metal line their transition wavelength and their relative strength with respect to the Ly$\alpha$ optical depth are given. The last column gives the position of maximum apparent correlation $r_\perp = 0$ and $r_\parallel \approx (1+z)D_H(z)(\lambda_m - \lambda_{Ly\alpha})/\lambda_\alpha$, corresponding to metal and Ly$\alpha$ absorption at the same physical position but reconstructed assuming only Ly$\alpha$ absorption at $z_{eff} = 2.334$.

For each pair ($m, n$) of contaminants, we compute the shifted–model correlation function with respect to the unshifted–model correlation function, $\xi_{\text{mod}}^{m-n}$, by introducing a metal matrix $M_{AB}$ (Blomqvist et al. 2018), such that

$$\xi_{\text{mod}}^{m-n}(A) \rightarrow \sum_{B} M_{AB} \xi_{\text{mod}}^{m-n}(\tilde{r}_{\parallel}(B), \tilde{r}_{\perp}(B)), \quad (35)$$

where

$$M_{AB} = \frac{1}{W_A} \sum_{(m,n) \in A, (m,n) \in B} w_m w_n, \quad (36)$$

where $W_A = \sum_{(m,n) \in A} w_m w_n$, ($m, n \in A$ refers to pixel separations computed assuming $z_\alpha$, and ($m, n \in B$ to pixel separations 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 in the computation of $w_m$ and $w_n$.

In the fits of the data, we include terms corresponding to metal–Ly$\alpha$ correlations for the four silicon transitions listed in Table 3. Since these correlations are only visible in ($r_\parallel, r_\perp$) bins corresponding to small physical separations, $b_m$ and $\beta_m$ cannot be determined separately. We therefore fix $\beta_m = 0.50$ as done in B17 based on measurements of the cross-correlation between DLAs and the Ly$\alpha$ forest (Font-Ribera et al. 2012b). Correlations between C IV and Ly$\alpha$ give a contribution outside the $r_\parallel$ range studied here, but we include effects of the C IV–C IV autocorrelation and fix $\beta_{C IV(eff)} = 0.27$ (Blomqvist et al. 2018).

### 4.4. Proximity Effect of Quasars on Ly$\alpha$ Absorption

The term in Equation (26) representing the transverse proximity effect takes the form (Font-Ribera et al. 2013)

$$\xi^{\text{TP}} = \xi_0^{\text{TP}} \left( \frac{1 h^{-1} \text{ Mpc}}{r} \right)^2 \exp(-r/\lambda_{\text{UV}}). \quad (37)$$

This form supposes isotropic emission from the quasars. We fix $\lambda_{UV} = 300 h^{-1}$ Mpc (Rudie et al. 2013) and fit for the amplitude $\xi_0^{\text{TP}}$.

### 4.5. Correlations due to Sky Subtraction

The correlation $\xi^{\text{sky}}$, of Equation (25), models the correlations induced by the sky subtraction procedure. The sky subtraction for each spectrum is done independently for each spectrograph, i.e., per half-plate, for the 500 fibers (450 science fibers and ~40 sky fibers). The Poisson fluctuations
<figure><img src="image_13.png" /><figcaption>Figure 9. Correlations due to the sky subtraction procedure in the Mg II(11) spectral region: $\lambda_{RF} \in [260, 276]$ nm. The correlations are shown for $\tilde{r}_{\perp} = 0$ as a function of $\tilde{r}_{\perp}$, where $(\tilde{r}_{\parallel}, \tilde{r}_{\perp})$ are the coordinates calculated assuming Mg II absorption. The green points show the correlations for pairs of pixels that are not taken with the same spectrograph (half-plate). Due to the small amount of Mg II(11) absorption, the correlations are consistent with zero, as expected. The blue points show the correlations for pairs of points taken with the same spectrograph ($\lambda_1 = \lambda_2$ pairs only). Here the sky subtraction procedure generates significant correlations that diminish with increasing angular separation (increasing $\tilde{r}_{\perp}$). The red points show the correlations for all pairs, and the red curve shows the fit using the form of Equation (38).</figcaption></figure>

in the sky spectra that are subtracted induce correlations in spectra obtained with the same spectrograph at the same observed wavelength, leading to an excess correlation in $r_{\parallel} = 0$ bins. This was observed in the DR12 autocorrelation measurement of B17, where they decided to reject such same-spectrograph pairs in their measurement. This effectively removes the excess correlation at $r_{\parallel} = 0$. However, because of continuum fitting, the excess correlation at $r_{\parallel} = 0$ generates a smooth distortion of the correlation function for all $r_{\parallel}$, and this was not removed by the procedure of B17. Here we do not remove the same-spectrograph pairs, which allows us to fit for its amplitude and thereby take into account the smooth component induced by continuum fitting. Apart from improving the fit to the data, this procedure is necessary because the DR16 analysis combines all measurements of each quasar, so spectra no longer correspond to a unique spectrograph.

The sky-fiber-induced correlations are easily seen at rest-frame wavelength longer than the $\text{Ly}\alpha$ quasar emission line, where correlations due to absorption are small. An example is shown in Figure 9, showing correlations in the “Mg II(11)” spectral region, $\lambda_{RF} \in [260, 276]$ nm (du Mas des Bourboux et al. 2019), where the pixels are immediately blueward of the Mg II quasar emission line, $\lambda_{RF} = 279.6$ nm. The correlation is consistent with zero for pairs of pixels taken with different spectrographs but is significant for pairs from the same spectrograph. This spurious correlation decreases rapidly with increasing angular separation. This decrease is due to the eBOSS pipeline that adds broadband functions to the sky calibration, modeling its variation across the size of the plate and thus decreasing the correlations with increasing angular separation.

Since, as the angular separation increases, the number of same-spectrograph pairs decreases and the number of different-

spectrograph pairs increases, the weighted sum of the two, shown as the red curve in the figure, can be modeled as a Gaussian function of $r_{\perp}$:

$$ \xi^{\text{sky}}(r_{\parallel}, r_{\perp}) = \begin{cases} \frac{A_{\text{sky}}}{\sigma_{\text{sky}}\sqrt{2\pi}} \exp\left(-\frac{1}{2}\left(\frac{r_{\perp}}{\sigma_{\text{sky}}}\right)^2\right), & \text{if } r_{\parallel} = 0 \\ 0, & \text{if } r_{\parallel} \neq 0 \end{cases}. \quad (38) $$

The two free parameters ($A, \sigma$)$_{\text{sky}}$ give the scale and the width of the correlation.

We fit the measured autocorrelation, Mg II(Mg II(11)) × Mg II(Mg II(11)), via the distorted model of Equation (40) and the model of the sky correlation of Equation (38). The fit is performed in the same condition as for the $\text{Ly}\alpha$ autocorrelation, hence $r \in [10, 180]$ h$^{-1}$ Mpc (Section 6), resulting in 1590 bins. We get a best fit with $\chi^2/(\text{DOF}) = 1551.94/(1590 - 2)$, corresponding to a probability $p = 0.76$ when using only the best observation. We find $\chi^2/(\text{DOF}) = 1531.49/(1590 - 2)$ and $p = 0.86$ when using all observations. In both cases, we find $(A, \sigma)_{\text{sky}} \sim (1 \times 10^{-3}, 20 \text{ h}^{-1} \text{ Mpc})$. A null test on the autocorrelation, using all observations, proves the very high significance of these two parameters: $\chi^2 = 2824.85$ and $\Delta\chi^2 = 1293.36$, for a difference of two parameters. Interestingly, a similar null test, removing the $r_{\parallel} = 0$ bins, still yields a significant measurement of the two parameters, $\Delta\chi^2 = 223.50$. This test shows that the effect of the sky subtraction residuals is mainly in the $r_{\parallel} = 0$ bins but that a nonnegligible fraction is brought to other bins by the effects of continuum fitting.

For the $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{Ly}\alpha(\text{Ly}\alpha)$ correlation function, the spurious correlation from the sky calibration in the $r_{\parallel} = 0$ bins reaches ≈20% of the physical correlation at $r_{\perp} \approx 50 \text{ h}^{-1} \text{ Mpc}$. For $r_{\parallel} \neq 0$, the distortion-induced correlation is smooth, so it does not affect our BAO measurement. However, it has a nonnegligible effect on the measurement of the $\text{Ly}\alpha$ bias parameters (Table D1), because of its $\mu$ dependence.

### 4.6. Power-law Correlations

An important test of systematic effects in the position of the BAO peak is performed by adding polynomial “broadband” terms to the correlation function. We follow the procedure and choice of broadband forms used by B17 and adopt the form

$$ B(r, \mu) = \sum_{j=0}^{j_{\text{max}}} \sum_{i=i_{\text{min}}}^{i_{\text{max}}} a_{ij} \frac{L_j(\mu)}{r^i} \quad (j \text{ even}), \quad (39) $$

where the $L_j$ are Legendre polynomials. Following B17, we fit with $(i_{\text{min}}, i_{\text{max}}) = (0, 2)$ corresponding to a parabola in $r^2\xi_{\text{smooth}}$ underneath the BAO peak. We set $j_{\text{max}} = 6$, giving four values of $j$ corresponding to approximately independent broadbands in each of the four angular ranges.

### 4.7. The Distorted Model

The expected measured correlation function, $\hat{\xi}$, is related to the true $\xi^t$, by the distortion matrix:

$$ \hat{\xi}_A = \sum_{A'} D_{AA'} \xi_{A'}^t. \quad (40) $$

The distortion matrix, $D_{AA'}$, is computed in Section 3.5 for the auto- and cross-correlations. The matrix accounts for the
correlations introduced by the continuum fitting procedure. In this equation, $A'$ is a bin of the model and $A$ is a bin of the measurement.<fnref n="34" />

# 5. Validation of the Analysis with Mocks

In this section we present a set of synthetic realizations of the eBOSS survey, and we use them to validate our BAO analysis. We start in Section 5.1 by presenting the methodology used to simulate the quasar positions and the Ly$\alpha$ forest fluctuations, and in Section 5.2 we describe how we simulate the quasar continua and instrumental artifacts. We present the BAO fits on the simulated data sets in Section 5.4.

## 5.1. Gaussian Random Field Simulations

Font-Ribera et al. (2012a) presented a method to efficiently simulate Ly$\alpha$ forest spectra for a given survey configuration. This method was used in multiple analyses of the Ly$\alpha$ autocorrelation in BOSS (Slosar et al. 2011, 2013; Busca et al. 2013; Delubac et al. 2015; Bautista et al. 2017), and it enabled multiple tests of potential systematics and the validation of the BAO analysis pipeline.

However, because the quasar positions behind the simulated forests did not correspond to density peaks, these simulated data sets did not have the correct cross-correlation between the quasars and the Ly$\alpha$ forest. This made them unfit for BAO studies using this cross-correlation. For this reason, the DR12 analysis of the cross-correlation in dMdB17 presented a new set of simulations using a different algorithm described in Le Goff et al. (2011), and that included a realistic cross-correlation. This allowed the first validation of BAOs in cross-correlation using mock data, and it showed that the auto- and cross-correlation measurements of the BAO scale have a very small correlation and can therefore be combined. These simulations, however, did not have the correct survey geometry and assumed instead that the lines of sight were parallel.

We have developed two different sets of simulations that have been generated in collaboration with the Ly$\alpha$ working group of DESI. The results presented in this section use simulations described in Farr et al. (2020), developed primarily at University College London and referred to here as the London mocks. A second set (T. Etourneau et al. 2020, in preparation) was developed primarily at CEA Saclay and is referred to here as the Saclay mocks. These mocks were completed after the London mocks and benefit from refinements inspired by experience with the London mocks. They will be used for future studies of the non-BAO component of the correlation function, including the effect of HCDs.

### 5.1.1. London Mocks

The methodology of the London mocks is described in Farr et al. (2020), and we refer the reader to that publication for further details. We use the CoLoRe package<fnref n="35" /> to generate a low-resolution Gaussian field on a very large box, of length $L \approx 10 \ h^{-1} \ \text{Gpc}$, enclosing an all-sky “light cone” to $z = 3.8$. The same package allows us to populate the peaks of the

Gaussian field with quasars, following an input quasar bias and number density. Finally, CoLoRe also provides the interpolated values of the Gaussian field on the lines of sight from the center of the box toward the quasars, as well as the interpolated radial velocity computed using the gradient of the gravitational potential.

The boxes used in this publication had 4096³ cells, resulting in a resolution of $\approx 2.4 \ h^{-1} \ \text{Mpc}$. As described in Farr et al. (2020), we use the Ly$\alpha$CoLoRe software<fnref n="36" /> to add extra small-scale fluctuations to each line of sight in order to reproduce the variance in the Ly$\alpha$ forest in the data. We then apply a lognormal transformation to the resulting Gaussian field and use the Fluctuating Gunn-Peterson Approximation (FGPA) to compute the optical depth in each cell. Finally, we use the radial velocities to include RSDs, and we compute the transmitted flux fraction for each Ly$\alpha$ spectrum.

### 5.1.2. Saclay Mocks

The methodology of the Saclay mocks is described in T. Etourneau et al. (2020, in preparation), and we refer the reader to that publication for further details. The Saclay mocks are similar to the London mocks. The main difference is that velocity gradient boxes are produced in addition to density field boxes, and RSDs are implemented by applying the lognormal transformation to the sum of the density and velocity gradient fields before applying FGPA. A nice feature of this implementation is that it allows for a prediction of the correlation function of the mocks, in a way similar to Font-Ribera et al. (2012a). The code<fnref n="37" /> does not benefit from the high level of parallelization of the CoLoRe package. The Gaussian field boxes are then smaller, and seven of them are needed to cover the eBOSS footprint.

## 5.2. Simulating eBOSS Spectra

Once we have the simulated transmitted flux fraction along each line of sight, we use these to simulate synthetic quasar spectra with the relevant astrophysical and instrumental artifacts. The methodology here is similar to the one used in the BOSS Ly$\alpha$ analyses, described in Bautista et al. (2015), and it can be summarized as follows:

1. Add contaminant absorption, including higher-order hydrogen lines, metal absorbers, and HCD systems.

2. Multiply each transmitted flux fraction by a simulated quasar continuum.

3. Convolve the simulated spectrum with the resolution of the spectrograph, pixelize it, and add Gaussian noise simulating eBOSS observations.

Figure 10 shows a simulated spectrum with different levels of complexity.

### 5.2.1. Adding Contaminants

In Section 5.1 we have described how we simulate our signal, the Ly$\alpha$ forest transmitted flux fraction. In order to study potential systematic biases caused by the presence of contaminants, we have the option to add other absorption lines from the Lyman series, as well as metal lines.

We simulate the optical depth of each contaminant as a rescaled version of the Ly$\alpha$ optical depth, mapped into a

<footnote n="34">³⁴ There are no theoretical requirements for the model to be computed using the same binning as the measurement. However, because of computing memory limitations, we decide to keep the same binning. In Appendix D we test the consequences of binning the distortion matrix on a finer grid with twice as many bins for the model as the measurement.</footnote>

<footnote n="35">³⁵ https://github.com/damonge/CoLoRe</footnote>

<footnote n="36">³⁶ https://github.com/igmhub/LyaCoLoRe</footnote>

<footnote n="37">³⁷ https://github.com/igmhub/SaclayMocks</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Example of synthetic spectra, for one of the brightest quasars in our sample. The green line shows the unabsorbed continuum; the solid blue and orange lines show the continuum multiplied by the mean transmitted flux fraction in the Lyα region and in the Lyβ region, respectively; the gray line shows the transmitted flux fraction generated with *LyaCoLoRe*, multiplied by a factor of 10 for better visualization; the black line shows the final spectrum, including instrumental noise. The vertical dashed blue line shows the position of a DLA, while its corresponding Lyβ absorption is marked with a vertical dashed orange line.</figcaption></figure>

different observed wavelength using its rest-frame wavelength. The scaling factors for the higher-order Lyman lines are computed using the oscillator strengths of each transition (see Equation (1.1) in Iršič et al. 2013). The scaling factor for Lyβ is 0.1901, and we include absorption up to the fifth Lyman line. The scaling factors for the metal lines have been tuned to approximately match the level of contamination observed in the data and are presented in Table 3. We have only included the four silicon lines that have a rest-frame wavelength similar to that of the Lyα and that are detected in the data in both the autocorrelation and the cross-correlation with quasars.

Using the method described in Farr et al. (2020), we generate a catalog of HCDs, with a column density and redshift distribution in agreement with current observations. We included systems with column densities in the range $log N_{H1} = [17.2, 22.5]$ using the software *pyigm* (Prochaska et al. 2017), itself calibrated using observations from Prochaska et al. (2014). The clustering of HCDs in the mocks has a large-scale bias consistent with the observed clustering of DLAs (Font-Ribera et al. 2012b; Pérez-Ràfols et al. 2018). A Voigt profile is then constructed for each absorber, and the absorption is included in the simulated transmitted flux fraction.

### 5.2.2. Simulating Quasar Continua

We assign a random magnitude to each simulated quasar, following the quasar luminosity function measured in Ross et al. (2013) using quasars from the ninth data release of SDSS, truncated at a maximum magnitude of $r = 21.3$. We then use the publicly available software package *simqso*<fnref n="38" /> to generate an unabsorbed continuum for each quasar. *simqso* is based on the simulations used in McGreer et al. (2013), and we refer the reader to that publication for a detailed description. In short, each continuum is constructed by adding a set of emission lines on top of a broken power law.

The distribution of quasar redshift errors is discussed in Section 2.2 and in Appendix B. Redshift errors have an important

impact in the cross-correlation of quasars and the Lyα forest, since they cause an extra smoothing of the correlations along the line of sight. In order to emulate this effect when measuring correlations in the synthetic data sets, we add a random redshift error drawn from a Gaussian distribution with $\sigma_z = 400$ km s⁻¹.

### 5.2.3. Simulating Instrumental Effects

The next step is to simulate the instrumental artifacts, including a Gaussian smoothing due to the finite spectral resolution of spectrograph, the pixelization of the spectra, and the instrumental noise.

While these effects have a dramatic impact on studies of the small-scale correlations of the Lyα forest, only the noise level has an impact on BAO measurements. This justifies our choice of *Specsim* (Kirkby et al. 2016), a simulator of the DESI spectrograph, to simulate the BOSS spectrograph. The small differences in the resolution and pixel size do not impact BAO measurements, and we have chosen an exposure time to reproduce the S/N in the eBOSS survey.

## 5.3. Comparison of Mocks and Data

As described in Farr et al. (2020), the synthetic Lyα spectra were tuned in order to reproduce large-scale measurements from the BOSS DR12 results, as well as reproducing the line-of-sight power spectrum over the relevant scales. In Appendix F we present figures comparing the measured correlations in the mocks and in the data, including for the first time comparisons of measurements of Lyα absorption in the Lyβ region.

An important role of the mocks is to verify the subsampling procedure for calculating the covariance matrix of the correlation functions. This can be done by comparing the subsampling calculation with the mock-to-mock variations of the correlation functions. We see no difference between the two at more than the 1% level. For example, in the left panel of Figure 11 we illustrate off-diagonal elements of the correlation (normalized covariance) matrix, which has been smoothed as described in Section 3.2 for the $Ly\alpha(Ly\alpha) \times Ly\alpha(Ly\alpha)$ and $Ly\alpha(Ly\alpha) \times$ quasar correlation functions. The correlations

<footnote n="38">³⁸ https://github.com/imcgreer/simqso</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Normalized mock covariances (correlations) measured by mock-to-mock variations and by subsampling and data covariances measured by subsampling. The left panel shows the $\Delta r_{\perp} = 0$ correlations vs. $\Delta r_{\parallel}$ for the $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{Ly}\alpha(\text{Ly}\alpha)$ and $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{quasar}$ correlation functions. The right panel shows the correlations between the $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{Ly}\alpha(\text{Ly}\alpha)$ and $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{quasar}$ functions. The largest correlations are those for $(\Delta r_{\perp} = 0, \Delta r_{\parallel} = 4 \text{ h}^{-1} \text{ Mpc})$ with $\text{Corr} \approx 0.4$ for the autocorrelation and $\approx 0.2$ for the cross-correlation. All correlations between the auto- and cross-correlations are at the percent level. For the mocks, the subsampling and mock-to-mock correlations agree at the percent level. The differences between mock and data correlations at $\Delta r_{\parallel} \approx 30 \text{ h}^{-1} \text{ Mpc}$ are due to different amounts of metals in the mocks and data, leading to differences in the intraforest correlation function.</figcaption></figure>

measured by subsampling and by mock-to-mock variations are
in excellent agreement, verifying our procedure.

Figure 11 shows that the covariances for mocks and data have a very similar structure, the differences being due in part to the smaller metal component in the mocks than in the data, leading to less correlation between nearby wavelengths. The mock variances are a factor of two smaller for the mocks for the autocorrelation and a factor of 1.5 for the cross-correlations. This is due to fewer low-flux spectra in the mocks compared to the data.

We ran all four analyses ($\text{Ly}\alpha$ autocorrelations and quasar cross-correlations using the $\text{Ly}\alpha$ and the $\text{Ly}\beta$ regions) on 10 realizations for each of the synthetic data sets described above. We analyze 90 extra realizations for two of the settings above (+continuum+noise) and ($+\text{HCDs}+\sigma_v$), allowing us to test possible systematic biases on the BAO measurements at a level 10 times smaller than our statistical uncertainty.

The main conclusion of this analysis is that we are able to
recover the right BAO scale even in the presence of
contaminants, and that our model is able to describe the
correlations measured on the synthetic data sets.

5.4. Analysis of BAOs in the Mocks

In Table 4 we present the results of our BAO analysis on multiple synthetic realizations of our data set. In order to study the effect of the different systematics included in the mock spectra, we run our analysis on different versions of each synthetic data set:

1. $\text{Ly}\alpha$: analysis run directly on the simulated transmitted flux fractions, including only $\text{Ly}\alpha$ absorption. This analysis does not need to do a continuum fitting and is therefore unaffected by the continuum distortions discussed in Section 3.5.

2. + continuum + noise: analysis run on simulated spectra, with quasar continua and instrumental noise.

3. + metals: analysis run on the same simulated spectra as above, but including also absorption from metal lines.

4. + HCDs + $\sigma_v$: analysis run on the same simulated spectra as above, but including also absorption from HCDs and including random redshift errors for the quasars. When analyzing the real data, we mask HCDs identified with an algorithm that is able to find the largest systems, with $\log N_{\text{H}\text{I}} > 20.3$. Similarly, when computing the delta field in the mock spectra, HCDs with $\log N_{\text{H}\text{I}} > 20.3$ are corrected using the input $N_{\text{H}\text{I}}$ and pixels with absorption larger than 20% are masked.<fnref n="39" />

5.4.1. Correlations in a BAO Measurement

Anisotropic BAO studies provide a measurement of the line-of-sight BAO scale, $\alpha_{\parallel}$, and of the transverse scale, $\alpha_{\perp}$. In BAO analysis using galaxy clustering these parameters are anticorrelated at the 40% level, i.e., they have a correlation coefficient around $\rho(\alpha_{\parallel}, \alpha_{\perp}) = -0.4$.

In Table 5 we use the BAO analyses on the 100 synthetic data sets to measure the correlation coefficient $\rho(\alpha_{\parallel}, \alpha_{\perp})$ in the four BAO measurements discussed in this publication. One can see that the results are also anticorrelated, with values ranging from -0.347 to -0.615. The correlations involving $(\alpha_{\parallel}, \alpha_{\perp})$ from distinct pairs of correlation functions are consistent with zero.

5.4.2. Covariance between Different Large-scale Correlations

The BAO results from the $\text{Ly}\alpha$ autocorrelation and from its cross-correlation with quasars were first combined in Font-Ribera et al. (2014). The authors used a Fisher matrix approach to argue that cosmic variance was subdominant in both results, and therefore the results could be considered independent.

This assumption was confirmed in du Mas des Bourboux et al. (2017), where the covariance between the measurements was measured to be very small on synthetic data sets. In the right panel of Figure 11 we present an updated study of this covariance when using the synthetic data sets described above.

<footnote n="39">³⁹ For one mock set, we have applied our DLA detector and verified that masking according to detected DLA parameters or to the input DLA parameters leads to nearly identical results.</footnote>
<table><caption>**Table 4**
Fit Results for Mock Data Sets</caption><thead><tr><th>Mock Set</th><th>α<sub>∥</sub> σ</th><th>α<sub>⊥</sub> σ</th><th>b<sub>η, Lyα</sub> σ</th><th>β<sub>Lyα</sub> σ</th><th>χ<sup>2</sup><sub>min</sub>/DOF, proba</th></tr></thead><tbody><tr><td>Lyα(Lyα) × Lyα(Lyα)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>raw</td><td>1.012 0.021</td><td>0.985 0.028</td><td>-0.2 0.001</td><td>1.568 0.021</td><td>1629.67/(1590-4) 0.27</td></tr><tr><td>+cont+noise</td><td>1.003 0.027</td><td>0.995 0.04</td><td>-0.201 0.002</td><td>1.486 0.028</td><td>1603.57/(1590-4) 0.41</td></tr><tr><td>+metals</td><td>1.012 0.029</td><td>0.987 0.05</td><td>-0.202 0.002</td><td>1.485 0.03</td><td>1600.32/(1590-8) 0.43</td></tr><tr><td>+HCD+σ<sub>v</sub></td><td>1.004 0.029</td><td>1.001 0.041</td><td>-0.205 0.003</td><td>1.48 0.051</td><td>1597.99/(1590-10) 0.42</td></tr><tr><td>Lyα(Lyα) × quasar</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>raw</td><td>1.008 0.025</td><td>0.999 0.024</td><td>-0.189 0.003</td><td>1.568 0.041</td><td>3227.31/(3180-6) 0.29</td></tr><tr><td>+cont+noise</td><td>1.008 0.029</td><td>0.992 0.033</td><td>-0.192 0.004</td><td>1.491 0.061</td><td>3203.5/(3180-6) 0.39</td></tr><tr><td>+metals</td><td>1.006 0.029</td><td>0.994 0.033</td><td>-0.193 0.004</td><td>1.51 0.063</td><td>3194.72/(3180-10) 0.4</td></tr><tr><td>+HCD+σ<sub>v</sub></td><td>1.003 0.033</td><td>0.998 0.033</td><td>-0.199 0.007</td><td>1.48 0.081</td><td>3212.8/(3180-10) 0.35</td></tr><tr><td>Lyα(Lyα) × Lyα(Lyβ)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>raw</td><td>1.005 0.025</td><td>0.996 0.034</td><td>-0.2 0.002</td><td>1.588 0.026</td><td>1609.71/(1590-4) 0.41</td></tr><tr><td>+cont+noise</td><td>1.014 0.049</td><td>0.983 0.069</td><td>-0.202 0.003</td><td>1.509 0.05</td><td>1592.39/(1590-4) 0.47</td></tr><tr><td>+metals</td><td>1.02 0.049</td><td>0.994 0.065</td><td>-0.203 0.004</td><td>1.528 0.054</td><td>1589.8/(1590-8) 0.46</td></tr><tr><td>+HCD+σ<sub>v</sub></td><td>1.009 0.054</td><td>1.019 0.087</td><td>-0.206 0.004</td><td>1.502 0.085</td><td>1609.88/(1590-10) 0.35</td></tr><tr><td>Lyα(Lyβ) × quasar</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>raw</td><td>1.028 0.042</td><td>1.009 0.044</td><td>-0.189 0.005</td><td>1.595 0.073</td><td>3189.37/(3180-6) 0.45</td></tr><tr><td>+cont+noise</td><td>1.008 0.07</td><td>1.015 0.082</td><td>-0.193 0.01</td><td>1.527 0.146</td><td>3183.49/(3180-6) 0.47</td></tr><tr><td>+metals</td><td>0.994 0.071</td><td>1.002 0.093</td><td>-0.19 0.01</td><td>1.495 0.149</td><td>3216.12/(3180-10) 0.36</td></tr><tr><td>+HCD+σ<sub>v</sub></td><td>1.011 0.08</td><td>1.013 0.099</td><td>-0.192 0.015</td><td>1.447 0.186</td><td>3195.57/(3180-10) 0.4</td></tr><tr><td>all combined</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>raw</td><td>1.009 0.012</td><td>0.995 0.014</td><td>-0.203 0.001</td><td>1.628 0.015</td><td>9948.08/(9540-7) 0.02</td></tr><tr><td>+cont+noise</td><td>1.005 0.017</td><td>0.992 0.022</td><td>-0.206 0.002</td><td>1.553 0.023</td><td>9788.16/(9540-7) 0.1</td></tr><tr><td>+metals</td><td>1.01 0.018</td><td>0.989 0.023</td><td>-0.206 0.002</td><td>1.558 0.025</td><td>9822.24/(9540-11) 0.08</td></tr><tr><td>+HCD+σ<sub>v</sub></td><td>1.005 0.019</td><td>0.998 0.023</td><td>-0.205 0.002</td><td>1.464 0.036</td><td>9625.56/(9540-13) 0.31</td></tr></tbody></table>

**Note.** Shown are the mean values and mean standard deviations ($\Delta\chi^2 = 1$) for $\alpha_\parallel$, $\alpha_\perp$, $b_\eta$, $b_{\eta, Ly\alpha}$, and $\beta_{Ly\alpha}$. Results are shown for “raw” mocks (no quasar continuum or noise added to the transmission field) and for three progressively more realistic mocks (adding continua, noise, metals, HCDs, and quasar velocity dispersion). Means are based on 100 mocks for (+cont+noise) and (HCD+σ<sub>v</sub>) and 10 otherwise.

The correlation coefficients are smaller than 2% for all combinations, probing the assumption that the measurements are independent.

## 6. Fits to the Data

This section presents the results of the fits to the data.<fnref n="40" /> We fit the two autocorrelation functions, $Ly\alpha(Ly\alpha) \times Ly\alpha(Ly\alpha)$ and $Ly\alpha(Ly\alpha) \times Ly\alpha(Ly\beta)$, and the two cross-correlation functions, $Ly\alpha(Ly\alpha) \times$ quasar and $Ly\alpha(Ly\beta) \times$ quasar. Table 6 lists the parameters used to fit the correlation functions and their best-fit values. Various combinations of the four functions have also been fit and are listed in the table.

### 6.1. The Correlation Functions

The correlation functions are fit for separations $r \in [10, 180] h^{-1} Mpc$, and for all directions, $\mu \in [0, 1]$ for the autocorrelations and $\mu \in [-1, 1]$ for the cross-correlations. This gives 1590 bins for the two autocorrelations and twice as many, 3180, for both cross-correlations. The different best-fit models differ in their number and the nature of the free parameters. However, they all share the two BAO parameters, ($\alpha_\parallel$, $\alpha_\perp$), that are the focus of this study.

The autocorrelation baseline model is composed of the Kaiser model for biased matter density tracers, the effect of

<table><caption>**Table 5**
Correlations between the BAO Parameters Measured in the Different Correlation Functions as Measured in Mocks</caption><thead><tr><th>Correlation Functions</th><th>ρ(α<sub>∥</sub>, α<sub>⊥</sub>)</th></tr></thead><tbody><tr><td>auto × auto</td><td>-0.572 ± 0.093</td></tr><tr><td>cross × cross</td><td>-0.470 ± 0.093</td></tr><tr><td>auto Lyβ × auto Lyβ</td><td>-0.615 ± 0.064</td></tr><tr><td>cross Lyβ × cross Lyβ</td><td>-0.347 ± 0.052</td></tr></tbody></table>

correlated sky residuals, the Kaiser model for metal-induced correlations, and finally the effect of HCDs. The model has 13 free parameters when fitting separately the $Ly\alpha(Ly\alpha) \times Ly\alpha(Ly\alpha)$ or $Ly\alpha(Ly\alpha) \times Ly\alpha(Ly\beta)$ functions. The simultaneous fit for the two functions has 15 free parameters because the two sky subtraction parameters are independent for the two functions, while the other parameters are common to them.

The cross-correlation baseline model is composed of the Kaiser model for quasars, $Ly\alpha$, and the contamination by metals, along with the modeling of the quasar systemic and statistical redshift errors. Because the RSD parameter of quasars, $\beta_{QSO}$, is highly correlated to the bias and beta parameters of $Ly\alpha$, we fix it, following Equations (33) and (34), at the effective redshift of the fit. We fix the HCD parameters at the values determined by the combined (auto + cross) fits since they are highly correlated with other

<footnote n="40">⁴⁰ We use the minimizer https://github.com/scikit-hep/iminuit.</footnote>
<table><caption>**Table 6**

Best-fit Parameters for Seven Different Fits to the Correlation Functions</caption><thead><tr><th>Parameter</th><th>Lyα(Lyα)<br>× Lyα(Lyα)</th><th>Lyα(Lyα)<br>×Lyα(Lyβ)</th><th>Lyα(Lyα)<br>×Quasar</th><th>Lyα(Lyβ)<br>×Quasar</th><th>Full Auto</th><th>Full Cross</th><th>All Combined</th></tr></thead><tbody><tr><td>α<sub>‖</sub></td><td>1.047 ± 0.035</td><td>1</td><td>1.059 ± 0.039</td><td>1</td><td>1.038 ± 0.032</td><td>1.056 ± 0.039</td><td>1.045 ± 0.023</td></tr><tr><td>α<sub>⊥</sub></td><td>0.980 ± 0.052</td><td>1</td><td>0.932 ± 0.047</td><td>1</td><td>0.959 ± 0.049</td><td>0.952 ± 0.042</td><td>0.956 ± 0.029</td></tr><tr><td>b<sub>η,Lyα</sub></td><td>-0.2009 ± 0.0039</td><td>-0.2045 ± 0.0065</td><td>-0.225 ± 0.010</td><td>-0.202 ± 0.024</td><td>-0.201 ± 0.0034</td><td>-0.2222 ± 0.0093</td><td>-0.2014 ± 0.0032</td></tr><tr><td>β<sub>Lyα</sub></td><td>1.657 ± 0.088</td><td>1.74 ± 0.16</td><td>1.95 ± 0.14</td><td>1.67 ± 0.30</td><td>1.669 ± 0.077</td><td>1.92 ± 0.13</td><td>1.669 ± 0.071</td></tr><tr><td>10<sup>3</sup>b<sub>η,Si II(119)</sub></td><td>-2.96 ± 0.50</td><td>-2.85 ± 0.96</td><td>-4.5 ± 1.2</td><td>3.0 ± 3.1</td><td>-2.94 ± 0.45</td><td>-3.6 ± 1.1</td><td>-2.62 ± 0.39</td></tr><tr><td>10<sup>3</sup>b<sub>η,Si II(119.3)</sub></td><td>-2.08 ± 0.50</td><td>-1.73 ± 0.93</td><td>2.1 ± 1.2</td><td>-3.6 ± 3.1</td><td>-2.02 ± 0.44</td><td>1.4 ± 1.1</td><td>-1.21 ± 0.39</td></tr><tr><td>10<sup>3</sup>b<sub>η,Si II(126)</sub></td><td>-2.2 ± 0.63</td><td>-4.1 ± 1.1</td><td>-1.81 ± 0.77</td><td>-1.7 ± 1.9</td><td>-2.68 ± 0.55</td><td>-1.8 ± 0.72</td><td>-2.29 ± 0.43</td></tr><tr><td>10<sup>3</sup>b<sub>η,Si iii(120.7)</sub></td><td>-4.54 ± 0.51</td><td>-4.34 ± 0.95</td><td>-0.98 ± 0.96</td><td>-0.5 ± 2.4</td><td>-4.5 ± 0.45</td><td>-0.92 ± 0.89</td><td>-3.72 ± 0.40</td></tr><tr><td>10<sup>3</sup>b<sub>η,C IV(ε<sup>eff</sup>)</sub></td><td>-5.2 ± 2.6</td><td>-5.1 ± 2.6</td><td>-4.8</td><td>-4.8</td><td>-5.2 ± 2.7</td><td>-4.8</td><td>-4.9 ± 2.6</td></tr><tr><td>b<sub>HCD</sub></td><td>-0.0522 ± 0.0044</td><td>-0.0556 ± 0.0078</td><td>-0.0501</td><td>-0.0501</td><td>-0.0523 ± 0.0039</td><td>-0.0501</td><td>-0.0501 ± 0.0036</td></tr><tr><td>β<sub>HCD</sub></td><td>0.610 ± 0.083</td><td>0.549 ± 0.087</td><td>0.703</td><td>0.703</td><td>0.646 ± 0.081</td><td>0.703</td><td>0.704 ± 0.080</td></tr><tr><td>β<sub>QSO</sub></td><td></td><td></td><td>0.2602</td><td>0.2602</td><td></td><td>0.2602</td><td>0.2601 ± 0.0059</td></tr><tr><td>Δr<sub>‖,QSO</sub> (h<sup>-1</sup> Mpc)</td><td></td><td></td><td>0.23 ± 0.13</td><td>0.39 ± 0.33</td><td></td><td>0.25 ± 0.12</td><td>0.10 ± 0.11</td></tr><tr><td>σ<sub>v</sub> (h<sup>-1</sup> Mpc)</td><td></td><td></td><td>7.73 ± 0.44</td><td>7.9 ± 1.1</td><td></td><td>7.77 ± 0.41</td><td>6.86 ± 0.27</td></tr><tr><td>ξ<sup>TP</sup><sub>0</sub></td><td></td><td></td><td>0.739</td><td>0.739</td><td></td><td>0.739</td><td>0.739 ± 0.092</td></tr><tr><td>10<sup>2</sup>A<sub>sky,auto</sub></td><td>0.941 ± 0.060</td><td></td><td></td><td></td><td>0.940 ± 0.058</td><td></td><td>0.930 ± 0.058</td></tr><tr><td>σ<sub>sky,auto</sub></td><td>31.4 ± 1.7</td><td></td><td></td><td></td><td>31.4 ± 1.7</td><td></td><td>31.5 ± 1.7</td></tr><tr><td>10<sup>2</sup>A<sub>sky,auto,Lyb</sub></td><td></td><td>1.32 ± 0.10</td><td></td><td></td><td>1.332 ± 0.095</td><td></td><td>1.323 ± 0.094</td></tr><tr><td>σ<sub>sky,auto,Lyb</sub></td><td></td><td>34.2 ± 2.4</td><td></td><td></td><td>34.1 ± 2.3</td><td></td><td>34.2 ± 2.3</td></tr><tr><td>N<sub>bin</sub></td><td>1590</td><td>1590</td><td>3180</td><td>3180</td><td>3180</td><td>6360</td><td>9540</td></tr><tr><td>N<sub>param</sub></td><td>13</td><td>11</td><td>10</td><td>8</td><td>15</td><td>10</td><td>19</td></tr><tr><td>λ<sup>2</sup><sub>min</sub></td><td>1604.79</td><td>1583.61</td><td>3238.96</td><td>3193.29</td><td>3190.88</td><td>6436.80</td><td>9654.56</td></tr><tr><td>probability</td><td>0.31</td><td>0.46</td><td>0.19</td><td>0.39</td><td>0.37</td><td>0.22</td><td>0.17</td></tr><tr><td>λ<sup>2</sup><sub>min,auto</sub></td><td></td><td></td><td></td><td></td><td>1603.74</td><td></td><td>1608.52</td></tr><tr><td>λ<sup>2</sup><sub>min,cross</sub></td><td></td><td></td><td></td><td></td><td></td><td>3235.31</td><td>3258.41</td></tr><tr><td>λ<sup>2</sup><sub>min,auto,Lyb</sub></td><td></td><td></td><td></td><td></td><td>1584.48</td><td></td><td>1585.48</td></tr><tr><td>λ<sup>2</sup><sub>min,cross,Lyb</sub></td><td></td><td></td><td></td><td></td><td></td><td>3196.39</td><td>3197.02</td></tr><tr><td>Δχ<sup>2</sup>(α<sub>‖</sub> = α<sub>⊥</sub> = 1)</td><td>2.00(0.81σ)</td><td>...</td><td>3.94(1.33σ)</td><td>...</td><td>1.72(0.77σ)</td><td>3.22(1.16σ)</td><td>4.62(1.55σ)</td></tr><tr><td>free A<sub>BAO</sub></td><td>1.37 ± 0.25<br>(5.48σ)</td><td>1.30 ± 0.41<br>(3.17σ)</td><td>1.08 ± 0.20<br>(5.4σ)</td><td>1.25 ± 0.58<br>(2.16σ)</td><td>1.32 ± 0.22<br>(6.0σ)</td><td>0.92 ± 0.19<br>(4.84σ)</td><td>1.10 ± 0.15<br>(7.33σ)</td></tr></tbody></table>

**Note.** The first four columns give fits to each measurement independently, while the last three give the results when running combined fits. All parameters are given at $z_{eff} = 2.334$. The first two lines give the BAO parameters, with error bars giving 68.27% of trials as estimated with fastMC, as explained in the text. Error bars of other best-fit parameters are given by minuit. Quantities without error bars are fixed in the fit. The second section gives the attributes of the best fit, along with the significance of the shift with respect to the Planck (2016) cosmology and the significance of the BAO peak.
parameters, especially the quasar velocity offset. We also fix the quasar proximity effect parameters to those determined by the combined fit. This results in a best-fit model with 10 parameters, for both cross-correlations. Since all these parameters are shared between $Ly\alpha(Ly\alpha) \times$ quasar and $Ly\alpha(Ly\beta) \times$ quasar, the combined fit is composed also of 10 free parameters.

The combined fit to the four measured correlation functions, two auto- and two cross-correlations, is composed of the different models cited in the two previous paragraphs. Running a combined fit to the auto- and cross-correlations allows us to break degeneracies and thus allows us to free $\beta_{QSO}$, to fit for HCDs in the cross-correlation and take into account the proximity effect of quasars onto $Ly\alpha$ absorption. The resulting best-fit model has 19 parameters.

When fitting simultaneously more than one of the four correlation functions, we neglect the cross-covariance between them. The subsampling covariance of Section 3.4 was at the level of 1%, and this was confirmed by the mock-to-mock variations shown in Figure 11. The expected low correlation between BAO parameters measured with different correlation functions was confirmed by the mock-to-mock variations discussed in Section 5.

We follow the prescriptions of dSA19 (their Appendix A) and compute the effective redshift of the measurement to be $z_{eff.} = 2.334$, when running the combined fit to all four correlations. This effective redshift is defined to be the pivot redshift for the two BAO parameters $\alpha_\parallel$ and $\alpha_\perp$. All parameters are assumed to be redshift independent except for the bias parameters (Section 4.2), which the fitter returns at the effective redshift.

The best-fit model of the individual fit to the autocorrelation $Ly\alpha(Ly\alpha) \times Ly\alpha(Ly\alpha)$ is given in 2D in the right panel of Figure 5 and in four wedges of $|\mu|$ in the top four panels of Figure 6. In a similar way the best-fit model of the cross-correlation $Ly\alpha(Ly\alpha) \times$ quasar is given in 2D in the right panel of Figure 8 and for four wedges of $|\mu|$ at the bottom of Figure 6. Finally, both auto- and cross-correlations, using the $Ly\beta$ region, $Ly\alpha(Ly\alpha) \times Ly\alpha(Ly\beta)$, and $Ly\alpha(Ly\beta) \times$ quasar, are given in Figure F1 of Appendix F.

The results of the seven different fits are given in Table 6. The first four columns give the results for each individual measured correlation. The last three columns give combined fits: “full auto,” “full cross,” and “all combined.”$^{41}$ In the first section, the first two lines give the BAO parameters $\alpha_\parallel$ and $\alpha_\perp$ then follows the different parameters that describe the overall shape of the correlation functions. The second section gives the attributes of the fit. The last two lines give properties of two different models. The second-to-last line is for a model where the BAO parameters are fixed to the cosmology of Planck (2016). The first number is the difference of $\chi^2$, and the second number in parentheses gives the conversion to the significance in number of trials, as computed using Monte Carlo realizations (Table 7). The last line gives the significance of the BAO peak.

In Table 6, parameters without error bars are fixed in the fit. This is the case for the BAO parameters using only the $Ly\beta$ region, where the BAO peak is not detected with high significance and we choose to fix $(\alpha_\parallel, \alpha_\perp) = (1, 1)$. In fits for the cross-correlation we fix the values of some parameters to the values found in the combined fit. For the HCD bias parameters, this is motivated by their high correlation with the quasar velocity distribution parameter $\sigma_v$. The parameters

<table><caption>**Table 7**
Values of $\Delta\chi^2$ Corresponding to CL = (68.27, 95.45%)</caption><thead><tr><th>Parameter</th><th>Δχ<sup>2</sup> (68.27%)</th><th>Δχ<sup>2</sup> (95.45%)</th></tr></thead><tbody><tr><td>Lyα(Lyα) × Lyα(Lyα)</td><td></td><td></td></tr><tr><td>α<sub>||</sub></td><td>1.06 ± 0.06</td><td>5.06 ± 0.35</td></tr><tr><td>α<sub>⊥</sub></td><td>1.24 ± 0.08</td><td>4.51 ± 0.25</td></tr><tr><td>(α<sub>||</sub>, α<sub>⊥</sub>)</td><td>2.64 ± 0.09</td><td>7.07 ± 0.33</td></tr><tr><td>Lyα(Lyα) × Lyα(Lyβ)</td><td></td><td></td></tr><tr><td>α<sub>||</sub></td><td>1.25 ± 0.09</td><td>5.1 ± 0.41</td></tr><tr><td>α<sub>⊥</sub></td><td>1.39 ± 0.07</td><td>5.03 ± 0.26</td></tr><tr><td>(α<sub>||</sub>, α<sub>⊥</sub>)</td><td>2.91 ± 0.11</td><td>7.9 ± 0.4</td></tr><tr><td>Lyα(Lyα) × quasar</td><td></td><td></td></tr><tr><td>α<sub>||</sub></td><td>1.21 ± 0.06</td><td>4.4 ± 0.22</td></tr><tr><td>α<sub>⊥</sub></td><td>1.22 ± 0.08</td><td>4.98 ± 0.26</td></tr><tr><td>(α<sub>||</sub>, α<sub>⊥</sub>)</td><td>2.59 ± 0.12</td><td>6.79 ± 0.31</td></tr><tr><td>Lyα(Lyβ) × quasar</td><td></td><td></td></tr><tr><td>α<sub>||</sub></td><td>1.43 ± 0.07</td><td>4.5 ± 0.18</td></tr><tr><td>α<sub>⊥</sub></td><td>1.41 ± 0.07</td><td>4.74 ± 0.2</td></tr><tr><td>(α<sub>||</sub>, α<sub>⊥</sub>)</td><td>2.93 ± 0.13</td><td>6.69 ± 0.22</td></tr><tr><td>auto all</td><td></td><td></td></tr><tr><td>α<sub>||</sub></td><td>1.06 ± 0.07</td><td>4.83 ± 0.46</td></tr><tr><td>α<sub>⊥</sub></td><td>1.14 ± 0.08</td><td>4.86 ± 0.29</td></tr><tr><td>(α<sub>||</sub>, α<sub>⊥</sub>)</td><td>2.42 ± 0.1</td><td>7.16 ± 0.33</td></tr><tr><td>cross all</td><td></td><td></td></tr><tr><td>α<sub>||</sub></td><td>1.23 ± 0.06</td><td>4.3 ± 0.19</td></tr><tr><td>α<sub>⊥</sub></td><td>1.11 ± 0.09</td><td>5.02 ± 0.41</td></tr><tr><td>(α<sub>||</sub>, α<sub>⊥</sub>)</td><td>2.64 ± 0.11</td><td>6.97 ± 0.36</td></tr><tr><td>combined</td><td></td><td></td></tr><tr><td>α<sub>||</sub></td><td>1.05 ± 0.05</td><td>4.07 ± 0.37</td></tr><tr><td>α<sub>⊥</sub></td><td>1.05 ± 0.08</td><td>4.37 ± 0.27</td></tr><tr><td>(α<sub>||</sub>, α<sub>⊥</sub>)</td><td>2.4 ± 0.07</td><td>6.4 ± 0.27</td></tr></tbody></table>

**Note.** Values are derived from 1000 Monte Carlo simulations of the correlation function that are fit using the model containing only $Ly\alpha$ absorption. Confidence levels are the fractions of the generated data sets that have best fits below the $\Delta\chi^2$ limit. The uncertainties are statistical and estimated using the bootstrap technique.

$b_{\eta, C_{IV}(eff)}$, $\beta_{QSO}$, and $\xi_0^{TP}$ are fixed because they are not significantly constrained by the cross-only fits.

We can notice that both widths, $\sigma_{sky}$, of the model of the correlated calibration residuals are of the same order of magnitude as that found for the Mg II(Mg II(11)) × Mg II(Mg II (11)) autocorrelation in Section 4.5. The strength, $A_{sky}$, is, however, 1 mag stronger, due to the difference in effective redshift of the correlation functions between $Ly\alpha$ and Mg II.

In another aspect we can note that the parameter describing the systematic quasar redshift error, $\Delta r_{\parallel, QSO} = 0.31 \pm 0.11 h^{-1}$ Mpc, is compatible with zero at the $3\sigma$ level. Using DR12 data, dMdB17 measured this parameter to be $-0.79 \pm 0.13 h^{-1}$ Mpc, which was $6\sigma$ from zero. This difference is linked to the change in the quasar redshift estimator (Section 2.2). Other estimators give even larger values, as shown in Table B1. Though this suggests that our new quasar redshift estimator is less biased, it is not possible to draw any definitive conclusion. Indeed, B19 showed that this parameter is highly correlated with models for relativistic effects (Lepori et al. 2020).

<footnote>$^{41}$ Best-fit values and likelihoods are also given in picca/tree/master/data/duMasdesBourbouxetal2020/fits.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 12.** The 68% and 95% confidence level contours in the ($\alpha_{\parallel}$, $\alpha_{\perp}$) plane from fits to the autocorrelation and cross-correlation functions and from the combined fit. The black dot shows the position of the flat-$\Lambda$CDM model of Planck (2016).</figcaption></figure>

## 6.2. Measurement of the BAO Parameters

Combining the two measurements of the Ly$\alpha$ autocorrelation yields

$$ \begin{cases} D_H(z = 2.334)/r_d = 8.93^{+0.28}_{-0.27}^{+0.64}_{-0.63} \\ D_M(z = 2.334)/r_d = 37.6^{+1.9}_{-1.9}^{+4.2}_{-4.0} \\ \rho(D_H(z)/r_d, D_M(z)/r_d) = -0.49 \end{cases} \quad (41) $$

Combining the two measurements of the Ly$\alpha$–quasar cross-correlation yields

$$ \begin{cases} D_H(z = 2.334)/r_d = 9.08^{+0.34}_{-0.34}^{+0.61}_{-0.60} \\ D_M(z = 2.334)/r_d = 37.3^{+1.7}_{-1.6}^{+4.2}_{-3.6} \\ \rho(D_H(z)/r_d, D_M(z)/r_d) = -0.43 \end{cases} \quad (42) $$

Finally, combining the two autocorrelations along with the two cross-correlations yields

$$ \begin{cases} D_H(z = 2.334)/r_d = 8.99^{+0.20}_{-0.19}^{+0.38}_{-0.38} \\ D_M(z = 2.334)/r_d = 37.5^{+1.2}_{-1.1}^{+2.5}_{-2.3} \\ \rho(D_H(z)/r_d, D_M(z)/r_d) = -0.45 \end{cases} \quad (43) $$

These three results are presented in Figure 12. The black point gives the Planck (2016) cosmology, while the contours give the 1 and 2 confidence levels, corresponding to (68.27%, 95.45%).

As first shown in dMdB17, the BAO parameters ($\alpha_{\parallel}$, $\alpha_{\perp}$) have non-Gaussian uncertainties. This means that, e.g., $\Delta\chi^2 = 1$ does not correspond to 68% of trials for one degree of freedom. To estimate the mapping between (68.27%, 95.45%) confidence levels to $\Delta\chi^2$ values, we generate 1000 fast Monte Carlo (fastMC) realizations of our seven fits. Each fit is reduced to the simple Kaiser model, in order to build enough statistics given the available computation time. A complete model would give similar results, as shown by dMdB17. In each realization, we generate a random noisy measurement of the best-fit model, given the covariance matrix, with BAO parameters set to 1. Then, the fit is performed in four different combinations: for the BAO parameters free, another for both of them fixed, then two more for one fixed and the other free. This allows for each of the seven results to give the error bars for $\alpha_{\parallel}$ and $\alpha_{\perp}$ independently. Table 7 gives the results of the estimation of the error bars. The BAO error bars given in Table 6 are the results of this fastMC correction.

The distribution of $\Delta\chi^2(\alpha_{\parallel} = \alpha_{\perp} = 1)$ in the 1000 fastMC allows us to compute the confidence levels in two dimensions. Furthermore, it allows us to estimate how significant the shift of our best-fit measurement is against the Planck (2016) cosmology. This result is given in Table 6 in the second-to-last line. Our final measurement, “all combined,” is $1.5\sigma$ from the cosmology of Planck (2016).

The significance of the BAO detection is estimated by leaving free the parameter $A_{\text{BAO}}$ that describes the size of the BAO peak relative to what is expected according to the Planck (2016) cosmology. We give the results for this extension to the baseline model in the last line of Table 6. Using fastMC realizations, we verify that this parameter is linear in its mapping from $\Delta\chi^2$ to confidence levels. Using 1000 fastMC, we verify this property up to $3\sigma$ and then assume that it holds true for higher confidence levels.

Our measurements of the autocorrelation, Ly$\alpha$(Ly$\alpha$) $\times$ Ly$\alpha$(Ly$\alpha$), and cross-correlation, Ly$\alpha$(Ly$\alpha$) $\times$ quasar, both have a BAO peak that has a significance of higher than $4\sigma$. It is barely $3\sigma$ for the autocorrelation Ly$\alpha$(Ly$\alpha$) $\times$ Ly$\alpha$(Ly$\beta$), and barely $2\sigma$ for the Ly$\alpha$(Ly$\beta$) $\times$ quasar. Because of this low significance of the BAO measurement in these two extra correlations, the $2\sigma$ error bars are very nonlinear, with Ly$\alpha$(Ly$\beta$) $\times$ quasar having none. We thus choose not to give the BAO measurement for these two correlations because it should not be used alone, but only in combination with other measurements. The final BAO measurement that combines all measured correlation functions has a BAO significance of more than $7\sigma$.

Our measured ($\alpha_{\parallel}$, $\alpha_{\perp}$) are marginally correlated to other parameters but significantly correlated, $\rho \sim -50\%$, with one another. The only other parameter that is correlated with ($\alpha_{\parallel}$, $\alpha_{\perp}$) is the bias parameter of the Si II(126) metal line. Indeed, the correlations of absorption by this metal and Ly$\alpha$ absorption at the same physical position generate an apparent peak in the reconstructed correlation function at ($r_{\parallel}$, $r_{\perp}$) $\sim$ (+111, 0) $h^{-1}$ Mpc (Table 3). This correlation is $\approx -22\%$ with $\alpha_{\parallel}$ and $+9\%$ with $\alpha_{\perp}$ in the case of the “full auto.” Since in the cross-correlations we have access to positive and negative distances along the line of sight, this correlation is reduced to $-14\%$ with $\alpha_{\parallel}$ and $+5\%$ with $\alpha_{\perp}$ for the “full cross.” Combining all results gives a correlation of $-19\%$ with $\alpha_{\parallel}$ and $+8\%$ with $\alpha_{\perp}$ for the “all combined” fit. All other parameters are less than $\pm 3\%$ correlated with the BAO parameters.

In Appendix D, we test for systematic errors in our measurement of the BAO parameters, along with changes in the estimation of their error bars. We either modify the best-fit model in Table D1 or change some aspects of the analysis or study data splits in Table D2. In Table D1 the strongest effect on both the recovered BAO best-fit value and error bars happens when taking the effect of metals into account, which results in a $0.5\sigma$ shift in $\alpha_{\parallel}$. As previously discussed, the Si II(126) absorption line produces a scale along the line of sight of $+111 h^{-1}$ Mpc. The change in the error bars is explained by the effect of preventing the best-fit model to build false BAO significance from this line.

Of the fits in Table D1, those that add polynomial “broadband” curves to the model are of special importance because they test the sensitivity of the BAO parameters to unidentified systematic errors in the model. We performed two fits, allowing the broadband curves different amounts of freedom in each. The first placed “physical priors” on ($b_{\text{Ly}\alpha}$, $\beta_{\text{Ly}\alpha}$, $b_{\text{HCD}}$) in the form of a Gaussian of mean and width of the fit without broadband terms. Such priors ensured
<figure><img src="image_17.png" /><figcaption>**Figure 13.** The 68% and 95% confidence level contours for the BAO parameters ($\alpha_\parallel$, $\alpha_\perp$) for the combined Ly$\alpha$ results from different data releases of SDSS: DR12 (B17, dMdB17, in blue), DR14 (dSA19, B19, in black), and DR16 (this work, in red). The black dot shows the position of the flat-ΛCDM model of Planck (2016).</figcaption></figure>

that the broadband terms were relatively small perturbations to the physical model. The second type of fit placed no priors on ($b_{\text{Ly}\alpha}$, $\beta_{\text{Ly}\alpha}$, $b_{\text{HCD}}$). The results of these fits are given in Table D1. We see that the addition of such terms does not change significantly the values of ($\alpha_\parallel$, $\alpha_\perp$) in any of the fits.

## 7. Comparison with Previous Ly$\alpha$ BAO Studies

The final Ly$\alpha$ BAO analyses from BOSS, using data from DR12, were presented in B17 and dMdB17. When combined, the measurement of the BAO scale was in mild ($\approx 2.3\sigma$) tension with the predictions from the best-fit ΛCDM cosmology from the Planck satellite (dMdB17). In the present study, combining all data from BOSS and eBOSS, the tension has now reduced to only $\approx 1.5\sigma$.

In Figure 13 we show the BAO contours for the combined Ly$\alpha$ BAO studies from different SDSS data releases: DR12 (B17, dMdB17), DR14 (dSA19, B19), and DR16 (this work). The contours have shifted toward the fiducial cosmology ($\alpha_\parallel = 1$, $\alpha_\perp = 1$), and the area of the contours has shrunk by roughly 25% between DR12 and DR16. It is important to remember, though, that the BAO uncertainty is itself a random variable, and it varies significantly from realization to realization.

In this section, we present first a comparison of the expected statistical power of the DR12 and DR16 data sets. We then discuss possible origins for the evolution of the BAO results from DR12 to DR16. The possibility that changes in the analysis pipeline induced significant changes is studied and rejected. We then investigate the effects of statistical changes, both in the enlarged data set and in the changed catalog of quasar redshifts. We will show that these two purely statistical effects are sufficient to explain the DR12-to-DR16 evolution.

We would like to have a comparison of the statistical power of the different data sets that does not suffer from this randomness. In order to do so, we generate Monte Carlo measurements of the different correlations that have the same covariance matrix as the real data, but where the measured correlations have been set to their theoretical prediction obtained using the combined best fit presented in the last column of Table 6. These mock correlations are then fitted using the same configuration as our main analyses, resulting in a fit that is clearly too good ($\chi^2 = 0$), but with parameter uncertainties that should be equivalent to those obtained in a Fisher forecast.

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 14.** Gain in statistical power on BAO parameters over that in DR12, for $\alpha_\parallel$ (top) and $\alpha_\perp$ (bottom). The blue bars show the constraining power of the Ly$\alpha$ autocorrelation in the Ly$\alpha$ region; orange bars add the cross-correlation with quasars, still using the Ly$\alpha$ region only (this result was not published for DR14); green bars show the autocorrelation using both the Ly$\alpha$ region and the Ly$\beta$ region (not computed in DR12); red bars show the combination of autocorrelation and cross-correlation with quasars, using both the Ly$\alpha$ region and the Ly$\beta$ region (not computed in DR12).</figcaption></figure>

In Figure 14 we show this gain in statistical power with respect to DR12 combined results, quantified as the ratio of uncertainties expected for the BAO parameters $\alpha_\parallel$ (top) and $\alpha_\perp$ (bottom). From this plot we can read that the statistical power in the Ly$\alpha$ autocorrelation (blue) of DR16 is 20%–25% larger than that in DR12 for both BAO parameters. These numbers are consistent with the average ratio of the error bars in the measured correlations. They are also comparable to the 27%–30% gain forecasted for eBOSS in Dawson et al. (2016), where we have used their forecasted values for a 4 yr survey with an area of 4500 square degrees, similar to the final eBOSS area. The gain is a bit smaller for the combination of the auto- and cross-correlation (orange), where the gain is roughly 15%–20% from DR12 to DR16. When we also include measurements in the Ly$\beta$ region, the gain in statistical power increases to 25%–30%. This gain corresponds to the average gain that we would see in a large ensemble of survey realizations, but each realization will have a different value. As seen in Figure 13, the
<table><caption>**Table 8**
Comparison of Lyα BAO Results Using DR12</caption><thead><tr><th>DR12 Analysis</th><th>α<sub>||</sub></th><th>α<sub>⊥</sub></th><th>χ<sup>2</sup><sub>min</sub>/DOF, proba</th></tr></thead><tbody><tr><td>Lyα(Lyα) × Lyα(Lyα):</td><td></td><td></td><td></td></tr><tr><td>Bautista et al. (2017)</td><td>1.053 ± 0.036</td><td>0.965 ± 0.055</td><td>1556.5/(1590-13), p = 0.639</td></tr><tr><td>Reanalysis (this work)</td><td>1.055 ± 0.036</td><td>0.987 ± 0.051</td><td>1582.59/(1590-13), p = 0.456</td></tr><tr><td>Lyα(Lyα) × quasar:</td><td></td><td></td><td></td></tr><tr><td>du Mas des Bourboux et al. (2017)</td><td>1.077 ± 0.038</td><td>0.898 ± 0.038</td><td>2576.3/(2504-15), p = 0.11</td></tr><tr><td>Reanalysis (this work, 10 &lt; r &lt; 160 h<sup>-1</sup> Mpc)</td><td>1.080 ± 0.037</td><td>0.896 ± 0.036</td><td>2544.58/(2504-15), p = 0.214</td></tr><tr><td>Reanalysis (this work, 10 &lt; r &lt; 180 h<sup>-1</sup> Mpc)</td><td>1.078 ± 0.037</td><td>0.893 ± 0.035</td><td>3292.87/(3180-10), p = 0.063</td></tr><tr><td>all combined:</td><td></td><td></td><td></td></tr><tr><td>du Mas des Bourboux et al. (2017)</td><td>1.069 ± 0.027</td><td>0.920 ± 0.031</td><td>3833.16/(3756-17), p = 0.14</td></tr><tr><td>Reanalysis (this work, 10 &lt; r &lt; 160 h<sup>-1</sup> Mpc)</td><td>1.074 ± 0.026</td><td>0.926 ± 0.030</td><td>3810.71/(3756-17), p = 0.203</td></tr><tr><td>Reanalysis (this work, 10 &lt; r &lt; 180 h<sup>-1</sup> Mpc)</td><td>1.066 ± 0.026</td><td>0.928 ± 0.031</td><td>4883.70/(4770-17), p = 0.091</td></tr></tbody></table>

**Note.** The results published in B17 and dMdB17 are compared to a reanalysis of the measurement using the data analysis pipeline used in this work. Top section: results for the Lyα autocorrelation; middle section: results for the cross-correlation; bottom section: combined results. Note that all analyses exclude the Lyβ region, as done in the DR12 publications. dMdB17 used a shorter range of separations in the fits (10 < r < 160 h⁻¹ Mpc), resulting in a different number of degrees of freedom. We show a reanalysis for both ranges of separations with the value of degrees of freedom in the last column specifying the range. BAO errors correspond to Δχ² = 1.

statistical gain between DR12 and DR16 in the actual realization observed is only 12%.

We now investigate the possibility that changes in the analysis pipeline could explain the evolution of the BAO results. The DR12 and DR16 results presented in Figure 13 used slightly different software to measure the correlations and different models to fit BAO. In order to address whether the shift is caused by the addition of new data, or by changes in the analysis pipeline, in Table 8 we compare the DR12 results published with our reanalysis of the DR12 data set, using the current analysis software picca, and the modeling described in the sections above. In order to better compare the analyses, we have used the exact DR12 quasar catalog and the Z<sub>VI</sub> redshifts that were used in the original DR12 analyses, as well as the DR12 DLA catalog. Similarly, we do not include Lyα information from the Lyβ region, as was done in the DR12 analyses.

Even though there are many differences in the analysis pipeline, the BAO results are in very good agreement, highlighting again the robustness of these measurements with respect to choices in the modeling (see Table D2).

We now turn to the question of whether the statistical changes due to the addition of new quasars and forests can explain the DR12-to-DR16 evolution. We estimated the expected changes in the BAO parameters using the fastMC technique described in Section 6.2. We first created 100 DR12 mock correlation functions, ξ<sub>12</sub>, Gaussian distributed about the best-fit DR16 model according to the observed DR12 covariance matrix, C<sub>12</sub>. For each DR12 correlation function, we created 100 correlation functions for eBOSS-only data using the noise-dominated covariance C<sub>16-12</sub> = (C<sub>16</sub><sup>-1</sup> - C<sub>12</sub><sup>-1</sup>)<sup>-1</sup>, where C<sub>16</sub> is the observed DR16 covariance matrix. Mock DR16 correlation functions were then created by adding the DR12 and eBOSS-only functions: ξ<sub>16</sub> = C<sub>16</sub>(C<sub>12</sub><sup>-1</sup>ξ<sub>12</sub> + C<sub>16-12</sub><sup>-1</sup>ξ<sub>16-12</sub>). The DR12-DR16 mock pairs were then fit for the BAO parameters. Approximately 30% of the mock pairs had changes in the BAO parameters that were greater than that observed for the combined fit in the data. We can then conclude that the observed changes in BAO parameters are consistent with those expected from statistical fluctuations.

A second source of statistical differences between DR12 and DR16 is our change in quasar redshift estimator from Z<sub>VI</sub> to Z<sub>LYAWG</sub>, made necessary by the lack of inspection for all DR16 quasars. This leads to random migration of quasar-pixel

pairs in r<sub>||</sub> space due to random differences between the two estimators. The statistical effect of this on the BAO parameters for the cross-correlation is studied in Appendix B, with the conclusion that it has a nonnegligible impact on DR12-DR16 evolution of BAO parameters for the cross-correlation.

## 8. Summary and Conclusions

This paper presents the measurement of BAOs with Lyα absorption and quasars from BOSS and eBOSS. It benefits from 10 yr of SDSS observations, from 2009 to 2019. Since the first measurement of the large-scale 3D Lyα autocorrelation and 3D Lyα-quasar cross-correlation, many improvements have been made in the analysis:

1. Slosar et al. (2011) measured for the first time the 3D large-scale autocorrelation of Lyα, from ~10,000 spectra, from the first few months of BOSS data. The autocorrelation was measured in the (r, μ) plane up to separations 100 h⁻¹ Mpc. The analysis also presented a set of mocks (Font-Ribera et al. 2012a) allowing tests of the measurement with the different astrophysical and observational effects: Poisson noise, quasar continuum, high absorption systems, metals. The first discussions of distortion and of the Wick expansion were presented.

2. Font-Ribera et al. (2013) measured for the first time the 3D large-scale cross-correlation of Lyα and quasars using the techniques introduced by Font-Ribera et al. (2012b) applied to ~60,000 DR9 spectra. The cross-correlation was measured in the (r<sub>||</sub>, r<sub>⊥</sub>) plane up to separations 80 h⁻¹ Mpc. The analysis estimated the bias of quasars to agree with other studies.

3. The DR9 autocorrelation measurement (Busca et al. 2013; Kirkby et al. 2013; Slosar et al. 2013) made the first detection and measurement of the BAO scale in the Lyα autocorrelation with DR9 and ~50,000 spectra. Two different studies measured this correlation either in the Cartesian coordinates (r, μ) or in the observable coordinates (Δ log λ, θ, z). Without any models for the distortion of the correlation, from the fit of the quasar continuum, the correlation was fit using broadband functions in addition to the Kaiser model.
4. The DR11 quasar cross-correlation measurement (Font-Ribera et al. 2014) was the first to observe the BAO peak in cross-correlation. Delubac et al. (2015) updated the autocorrelation measurement, and combined constraints on BAO parameters were given.

5. The DR12 autocorrelation (B17) and cross-correlation (dMdB17) measurements used the complete BOSS data set. They performed the important breakthrough of modeling the effect of the distortion of the correlation from the fit of the quasar continuum, allowing for the first physical fit of the 3D correlation. Broadband functions were kept to test for systematics in the measurement of BAO but no longer played a main role in the analysis. In addition to improvements in the data analysis, a main effort was dedicated to improving the realism of the mock spectra (Bautista et al. 2015) for metals and HCDs. The first mocks with realistic quasar-forest correlations were produced (Le Goff et al. 2011). The first combined fit between the auto- and cross-correlations was performed, breaking parameter degeneracies. Covariances were better understood through the study of Wick expansions.

6. The DR14 autocorrelation (dSA19) and cross-correlation (B19) measurements were the first to use eBOSS data. Developments were made to benefit from Lyα absorption blueward of the Lyβ+O VI emission line in ~60,000 spectra. The analysis used all observations of the same quasar, instead of the best one as was done in previous analyses. Furthermore, the analysis investigated different models for the effect of HCDs on the autocorrelation function. In addition, the analysis developed a new way of splitting the sample in order to measure the autocorrelation in two redshift bins, while limiting the cross-covariance to subpercent levels.

7. The DR16 study presented here is the first to use the complete BOSS and eBOSS data sets. It uses ~210,000 spectra for the Lyα spectral region and 70,000 for the Lyβ region. The analysis gives BAO for both the auto- and cross-correlations. This study focuses on improving the understanding of the effects of the calibration from sky and standard-star fibers, expanding the Lyβ spectral region, better understanding the effects of the quasar redshift estimator, and continuing the development of realistic 3D Gaussian random field mocks (T. Etourneau et al. 2020, in preparation; Farr et al. 2020).

With the beginning of the DESI (DESI Collaboration et al. 2016) and WEAVE (Pieri et al. 2016) projects, the measurement of BAOs with Lyα forest will continue. In order to profit from the improved statistical precision, it will be useful to continue to improve various parts of the analysis, among which are better understanding of quasar redshift estimators (Section 2.2) and improved mock spectra through improved modeling of quasar spectral diversity and absorption by metals (Section 5). Improved modeling of the correlation function is essential in light of the fact that including additional ad hoc polynomial broadband terms improves the fit, though without affecting the BAO parameters. Improvements here will allow us to profit fully from the full shape of the correlation function and yield a reliable measurement of the growth rate of cosmological structure (du Mas des Bourboux 2017).

The measurement of BAO at redshift $z \approx 2.4$ has had an important impact on cosmology. While these data do not

significantly constrain simple models when CMB (Planck Collaboration et al. 2016) and galaxy BAO data (Alam et al. 2017) are already used, the Lyα BAO data are essential for providing a measurement of ΛCDM parameters using only low-redshift data (Aubourg et al. 2015). They thus provide an independent check on the CMB-inspired flat ΛCDM model. The cosmological constraints from the DR16 data presented here are given in a companion paper (eBOSS Collaboration et al. 2020).

It is a pleasure to thank Nicolas Busca for spearheading this analysis over an 8 yr period, during which he energetically provided creative solutions to innumerable issues.

The work of H.d.M.d.B. and K.D. was supported in part by U.S. Department of Energy, Office of Science, Office of High Energy Physics, under award No. DESC0009959. A.F.-R. acknowledges support by an STFC Ernest Rutherford Fellowship, grant reference ST/N003853/1, and by FSE funds through the program Ramon y Cajal (RYC-2018-025210) of the Spanish Ministry of Science and Innovation.

This work was supported by the A\*MIDEX project (ANR-11-IDEX-0001-02) funded by the “Investissements d’Avenir” French Government program, managed by the French National Research Agency (ANR), and by ANR under contracts ANR-14-ACHN-0021 and ANR-16-CE31-0021.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS website is www.sdss.org. SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration, including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatório Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

In addition, this research relied on resources provided to the eBOSS Collaboration by the National Energy Research Scientific Computing Center (NERSC). NERSC is a U.S. Department of Energy Office of Science User Facility operated under contract No. DE-AC02-05CH11231.
<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>**Figure A1.** Autocorrelation $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{Ly}\alpha(\text{Ly}\alpha)$ (left) and cross-correlation $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{quasars}$ (right) for the wedge $0.8 < |\mu| < 0.95$ as a function of the separation in Mpc. The blue points and curves show these correlations assuming Planck (2016), Table 2, while the red points and curves show them assuming the cosmology that measures $\alpha_\parallel = \alpha_\perp = 1$ for the autocorrelation on the left and the different cosmology that measures $\alpha_\parallel = \alpha_\perp = 1$ for the cross-correlation on the right.</figcaption></figure>

<table><caption>**Table A1**
BAO Parameters for Different Assumed Cosmologies</caption><thead><tr><th>Cosmology</th><th>α<sub>||</sub></th><th>α<sub>⊥</sub></th><th>D<sub>H</sub>(z = 2.33)/r<sub>d</sub></th><th>D<sub>M</sub>(z = 2.33)/r<sub>d</sub></th><th>χ<sup>2</sup><sub>min</sub>/DOF, proba</th></tr></thead><tbody><tr><td>Lyα(Lyα) × Lyα(Lyα):</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Planck (2016, baseline)</td><td>1.043 ± 0.035</td><td>0.986 ± 0.051</td><td>8.97 ± 0.30</td><td>38.6 ± 2.0</td><td>1599.55/(1590-14), p = 0.33</td></tr><tr><td>new cosmo. auto</td><td>0.997 ± 0.035</td><td>1.003 ± 0.053</td><td>8.94 ± 0.32</td><td>38.7 ± 2.1</td><td>1631.33/(1590-14), p = 0.16</td></tr><tr><td>Lyα(Lyα) × quasar:</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Planck (2016, baseline)</td><td>1.058 ± 0.042</td><td>0.929 ± 0.053</td><td>9.10 ± 0.36</td><td>36.4 ± 2.1</td><td>3232.82/(3180-10), p = 0.21</td></tr><tr><td>new cosmo. cross</td><td>1.002 ± 0.040</td><td>1.002 ± 0.057</td><td>9.11 ± 0.36</td><td>36.5 ± 2.1</td><td>3211.82/(3180-10), p = 0.30</td></tr></tbody></table>

**Note.** The cosmologies are “Planck (2016),” the baseline of the study of this work, “new cosmo. auto,” the custom cosmology for the $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{Ly}\alpha(\text{Ly}\alpha)$ autocorrelation, and “new cosmo. cross,” the custom cosmology for the $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{quasars}$ cross-correlation. The error bars on the BAO parameters are given without fastMC corrections. The last column gives the goodness of the fit.

## Appendix A
Effect of the Assumed Cosmology

The BAO parameters ($\alpha_\parallel$, $\alpha_\perp$) are given in Table 6 with respect to the assumed cosmology of Planck (2016). For this reason they depend on the assumed cosmology. However, as given in Equation (24), the parameters ($D_H(z)/r_d$, $D_M(z)/r_d$) are independent of this assumption, as long as the real cosmology is not “too far” from the Planck (2016) cosmology. Here we demonstrate this property for the auto- and cross-correlations by replacing the fiducial cosmology of Table 2 with another model that yields ($\alpha_\parallel$, $\alpha_\perp$) = (1, 1). This test has been performed on mock galaxy catalogs (Carter et al. 2020), but this is the first test using $\text{Ly}\alpha$ correlations.

We modify the Planck (2016) model by changing the values of ($h$, $\Omega_k$) while maintaining the values of ($\Omega_c h^2$, $\Omega_b h^2$), which are precisely determined by the CMB anisotropies, independent of the curvature and dark energy model. This maintains the values of $r_d$ but changes the functions $D_M(z)$ and $D_H(z)$. For a given ($\alpha_\parallel$, $\alpha_\perp$) found using the Planck (2016) model (Table 6), we can predict the values of ($h$, $\Omega_k$) that allow us to find ($\alpha_\parallel$, $\alpha_\perp$) = (1, 1). We do this separately for the auto- and cross-correlations. As shown in Table A1, with the new cosmologies we recover the same BAO parameters ($D_H(z)/r_d$, $D_M(z)/r_d$).

The two 3D correlations, auto-correlation $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{Ly}\alpha(\text{Ly}\alpha)$ and cross-correlation $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{quasars}$, are treated independently in this section. Using the Planck (2016) cosmology of Table 2, we compute ($D_H(z = 2.33)/r_d$, $D_M(z = 2.33)/r_d$) for different values of ($h$, $\Omega_k$). We find that

($h$, $\Omega_k$) = (0.70743, -0.09837) and ($h$, $\Omega_k$) = (0.80829, -0.12035) should reproduce the auto- and cross-correlation BAO results (Equations (41) and (42)) but with ($\alpha_\parallel$, $\alpha_\perp$) = (1, 1). Both 3D correlations, with their distortion matrix and their metal matrix, are computed using the updated new cosmological parameters ($\Omega_m$, $\Omega_r$, $\Omega_k$) to map the observable coordinates ($\Delta\lambda$, $\Delta\theta$) to their respective Cartesian coordinates ($r_\parallel$, $r_\perp$). In order to fit these new correlations, we compute their respective matter power spectrum at $z = 2.33$ with CAMB; then, using *picca* and *nbodykit*, we extract the wiggle-less power spectrum of the sideband (Section 4). The resulting measured and best-fit correlation function is shown for the wedge near the line of sight in Figure A1. For these different cosmologies, the BAO scale, $r_d$, is unchanged at better than 0.01%; however, the assumed values for the Hubble parameter are very different, and thus we show the separation on the x-axis in Mpc as opposed to $h^{-1}$ Mpc as was done in Figure 6. This allows us to present the BAO scale at nearly the same position for all the different assumed cosmologies.

Table A1 summarizes the best-fit BAO results. As expected, the two customized cosmologies for the auto- and cross-correlations give ($\alpha_\parallel$, $\alpha_\perp$) = (1, 1). We also observe that the $\alpha_i$ parameters are dependent on the assumed cosmology. On the other hand, the BAO parameters ($D_H(z = 2.33)/r_d$, $D_M(z = 2.33)/r_d$) are measured to be the same regardless of the assumed cosmology; furthermore, they are measured with the same error bars.
<table><caption>Table B1
Characteristics of Six Quasar Redshift Estimators as Described in the Text</caption><thead><tr><th>Redshift<br>Estimator</th><th>Estimator - Z_LYAWG<br>Mean (rms)</th><th>Δr<sub>‖,QSO</sub><br>(h<sup>-1</sup> Mpc)</th><th>σ<sub>v</sub><br>(h<sup>-1</sup> Mpc)</th><th>χ<sup>2</sup><sub>min</sub>/DOF, proba</th></tr></thead><tbody><tr><td>Z_LYAWG</td><td>0 (0)</td><td>0.31 ± 0.11</td><td>6.88 ± 0.28</td><td>9626.45/(9540-20), p = 0.22</td></tr><tr><td>Z_PCA</td><td>-0.0022 (0.0077)</td><td>-0.84 ± 0.11</td><td>6.81 ± 0.31</td><td>9734.68/(9540-20), p = 0.061</td></tr><tr><td>Z_CIV</td><td>-0.0059 (0.0084)</td><td>-4.31 ± 0.12</td><td>8.54 ± 0.34</td><td>9746.02/(9540-20), p = 0.052</td></tr><tr><td>Z_CIII</td><td>-0.0066 (0.0113)</td><td>-4.34 ± 0.11</td><td>6.90 ± 0.33</td><td>9786.63/(9540-20), p = 0.027</td></tr><tr><td>Z</td><td>3·10<sup>-6</sup> (0.0076)</td><td>-0.783 ± 0.099</td><td>4.94 ± 0.23</td><td>9730.85/(9540-20), p = 0.064</td></tr><tr><td>Z_VI (DR12)</td><td>0.0003 (0.0087)</td><td>-0.79 ± 0.14</td><td>4.46 ± 0.30</td><td>4853.35/(4770-18), p = 0.15</td></tr></tbody></table>

Note. The second column gives the mean and standard deviation of the difference between the estimator and Z_LYAWG, the estimator used in this study. All mean and standard deviations are calculated for quasars in DR16Q with redshifts >1.77 and after eliminating outliers with $|\Delta z_q| > 0.05$. The third and fourth columns give the best-fit values of the parameters $\Delta r_{\parallel, QSO}$ and $\sigma_v$. The relatively large values of the means of Z_CIV-Z_LYAWG and Z_CIII-Z_LYAWG in Column (2) are correlated with the large values of $\Delta r_{\parallel, QSO}$ in Column (3), with the expected constant of proportionality, $D_H(z_{eff}) = 852$ Mpc/h. The results for Z_LYAWG differ slightly from those of the baseline fits (Table 6) because in this table the HCD model of B17 was used.

## Appendix B
Quasar Redshifts

Quasar redshifts enter our study in two ways. For the autocorrelation, they are used only to define the continuum model over the $Ly\alpha$ and $Ly\beta$ spectral ranges. Errors on quasar redshifts therefore only add an effective quasar spectral diversity by moving features in the spectral template randomly in wavelength. This increases the noise in the measurement of the transmission field but does not systematically shift the ($r_{\parallel}$, $r_{\perp}$) of pixel pairs. On the other hand, for the cross-correlation, random errors in quasar redshift estimates smear the correlation function in the $r_{\parallel}$ direction, while a systematic under- or overestimate of redshifts would produce an asymmetry between positive and negative $r_{\parallel}$. Both of these effects are accounted for in the fits: the smearing via the parameter $\sigma_v$ (Equation (30)) and the asymmetry via the parameter $\Delta r_{\parallel}$. As such, we do not expect quasar redshift errors to systematically bias BAO parameters.

In this appendix we quantify the effects of redshift uncertainties on several redshift estimators by using mock and real spectra. The different quasar redshift estimates that are included in the DR16Q catalog are described in Lyke et al. (2020). The important estimators are listed in Table B1 with their best-fit values of $\Delta r_{\parallel, QSO}$ and $\sigma_v$. The Z_LYAWG, Z_PCA, Z_CIV, and Z_CIII estimates are derived by the code redvsblue as described below. A composite redshift estimate, called Z, was derived primarily from the visual inspection redshift, Z_VI, and the pipeline redshift, Z_Pipe, using a prior from the neural network quasarNET,<fnref n="42" /> (Busca & Balland 2018). Z_Pipe was derived by fitting a library of spectral templates based on principal component analysis (PCA). Z_VI was estimated on all quasar spectra in the DR12 sample (Pâris et al. 2017). Visual inspections were performed for only a subset of the spectra obtained during eBOSS observations to reduce the rate of false detections (Pâris et al. 2018).

The DR12 analyses (B17, dMdB17) used Z_VI as the redshift estimator, while the DR14 analyses (dSA19, B19) used its closest descendant, Z. All eBOSS DR16 BAO analyses, including this one, use PCA-derived redshifts because the lack of systematic visual scanning makes the Z estimate inhomogeneous over the whole sample.

A known shortcoming of the pipeline redshift is the lack of redshift evolution of $Ly\alpha$ absorption in the spectral templates. To mitigate the biases resulting from incorrect models of the quasar spectrum at blue wavelengths, we developed the publicly available code “redshift versus blueshift”: redvsblue.<fnref n="43" /> As with the Z_Pipe redshifts, this code performs redshift estimates using four PCA eigenvectors and three nuisance terms represented by Legendre polynomials. However, redvsblue can correct the eigenvectors for the redshift evolution of the mean $Ly\alpha$ optical depth. The fits can also be restricted to a chosen redshift range informed by the Z estimate. The estimates Z_CIV and Z_CIII restrict the range to a given emission line. The estimate Z_LYAWG used in this study restricts the range to rest-frame wavelengths greater than 135.6 nm, i.e., near the beginning of rise of the spectrum toward the $Ly\alpha$ emission peak. The estimate Z_PCA does not make this wavelength restriction but corrects the model for the $Ly\alpha$ optical depth of Calura et al. (2012).

We assessed the accuracy and precision of the redvsblue estimates using the mock spectra of Section 5. Figure B1 shows the standard deviation (dashed curves) and mean difference (solid curves) between the input redshift and the estimated redshift for four different configurations of redvsblue. The four configurations labeled in the figure are (1) “no $Ly\alpha$ Corr,” which is nearly identical to Z_Pipe; (2) “$Ly\alpha$ Mask” (Z_LYAWG); (3) “$Ly\alpha$ Corr” (Z_PCA); and (4) “$Ly\alpha$ Corr, diff model,” also correcting for the evolution of the $Ly\alpha$ optical depth but using the model of Kamble et al. (2020).

Figure B1 shows that not correcting for the $Ly\alpha$ absorption in the spectral models leads to systematic errors and statistical errors that increase with redshift. The bias in redshift estimates is reduced to values less than 100 km s⁻¹ when using either correction for the $Ly\alpha$ absorption, albeit with slightly different signatures of the systematic errors. Of the four methods, the statistical errors are minimized when masking the spectra around the $Ly\alpha$ emission line and at shorter wavelengths. This improvement can likely be explained by variance in LSS variance and absorption from structures proximate to the quasar that lead to distortions of the $Ly\alpha$ emission line that are difficult to model. These results motivated our choice of Z_LYAWG for the DR16 BAO analysis presented in this paper.

<footnote n="42">42 https://github.com/ngbusca/QuasarNET</footnote>

<footnote n="43">43 https://github.com/londumas/redvsblue</footnote>
<figure><img src="image_22.png" /><figcaption>**Figure B1.** Velocity difference between the estimated quasar redshift and the true input redshift of mocked quasar spectra as a function of the true redshift. Solid curves give the mean of the distribution; dashed curves give the standard deviation. The four redvsblue redshift estimators are Z_pipe' (no correction for redshift evolution of Lyα optical depth), Z_LYAWG (masking of Lyα emission and forest wavelengths), Z_PCA (correcting for redshift evolution of Lyα optical depth using the model of Calura et al. 2012), and Z_PCA' using the model of Kamble et al. (2020).</figcaption></figure>

Table B1 shows the effect of using different quasar redshift estimators on the nuisance variables $\Delta r_{\parallel, QSO}$ and $\sigma_v$. The largest $r_{\parallel}$ asymmetries, reflected in $\Delta r_{\parallel, QSO}$, appear in the Z_CIV and Z_CIII estimates. The values of $\Delta r_{\parallel, QSO}$ relative to that for Z_LYAWG are correlated with the difference between the estimator and Z_LYAWG (Column (2)) with the expected factor of proportionality, $D_H(z_{\text{eff}}) = 852 h^{-1} \text{ Mpc}$.

The best-fit values of ($\alpha_{\parallel}$, $\alpha_{\perp}$) for four alternative estimators, along with those for the adopted Z_LYAWG, are listed in Table D2. We do not expect the choice of quasar redshift estimator to systematically change the BAO parameters in one way or another because any small mean shift in quasar redshifts would be absorbed into the nuisance parameter $\Delta r_{\parallel}$. Pairs of estimators do, however, have rms differences of $\Delta z_q \approx 0.008$. Use of different estimators would randomly shift quasar–forest pixel pairs in $r_{\parallel}$ space with $\Delta r_{\parallel} = 0.008 D_H \approx 7 h^{-1} \text{ Mpc}$, corresponding to $\approx 2$ bins in $r_{\parallel}$. This effect results in random changes in the cross-correlation function, resulting in random changes in the BAO parameters.

We estimated the range of expected differences in BAO parameters when using different quasar redshift estimators by using 20 Saclay mock sets (Section 5). For each mock set, a pair of new sets was created, with the only change being in the quasar redshift: for each quasar the mean of the two redshifts was equal to the original redshift, and the difference of the two redshifts was drawn randomly from the observed distribution of Z-Z_LYAWG in the DR16Q catalog (rms = 0.0076 from Table B1, Column (2)). The mock pairs were then fit for $\alpha_{\parallel}$ and $\alpha_{\perp}$. The rms differences between pairs were 0.49σ and 0.35σ for $\alpha_{\parallel}$ and $\alpha_{\perp}$, respectively. A $\chi^2$ statistic characterizing the overall change in BAO parameters was calculated. Of 20 mock pairs, 8 had changes greater than the observed difference between the use of Z and Z_LYAWG DR12 (Table D2). This indicates that the differences observed for different estimators are statistical in origin. Furthermore, four of the mock pairs had changes greater than the observed change in the Lyα(Lyα) × quasar BAO parameters going from DR12 (Table 8) to DR16 (Table 6). This indicates that the change in redshift estimator had a nonnegligible impact on the evolution from DR12 to DR16.

## Appendix C
### Tests of the Covariance Matrix

The covariance matrix for the correlation function was calculated in Section 3 by the subsampling methods. The accuracy of this technique was verified with mock spectra (Section 5) by comparing the subsampling covariance with the mock-to-mock variations of the correlation functions. Here we present two other techniques to estimate the covariance, the “Wick expansion” and “shuffling.”

The Wick method (Delubac et al. 2015, B17, dMdB17) uses the expansion of the four-point function in terms of products of two-point functions:<fnref n="44" />

$$ \frac{1}{W_A W_B} \sum_{(i,j) \in A} \sum_{(k,l) \in B} w_i w_j w_k w_l \langle \delta_i \delta_j \delta_k \delta_l \rangle \\ \approx \frac{1}{W_A W_B} \sum_{(i,j) \in A} \sum_{(k,l) \in B} \sum_{(\alpha, \beta, \mu, \nu) \in (i,j,k,l)} \\ \times w_{\alpha} w_{\beta} w_{\mu} w_{\nu} \langle \delta_{\alpha} \delta_{\beta} \rangle \langle \delta_{\mu} \delta_{\nu} \rangle, \qquad (C1) $$

where A and B are two bins of the autocorrelation and $W_{A,B} = \sum_{(i,j) \in A, B} w_i w_j$. The covariance, $C_{AB}$, is the left-hand side minus the $(i, j, k, l) = (\alpha, \beta, \mu, \nu)$ term on the right-hand side.

Calculation of all terms in Equation (C1) is computationally expensive. Fortunately, the terms involving only two forests dominate, in which case the covariance is

$$ C_{AB} = \frac{1}{W_A W_B} \sum_{ij \in A} \sum_{kl \in B} w_i w_j w_k w_l \xi_{1d}(\lambda_i / \lambda_k) \xi_{1d}(\lambda_j / \lambda_l), \quad (C2) $$

where $\xi_{1d}$ is the intraforest correlation function.

The resulting correlation matrix is presented in Figure C1, by the orange curve. There, “Diagram T123” refers to the decomposition of Equation (C1) into different configuration diagrams, as presented in Figure A.1 of Delubac et al. (2015), and where T1, T2, and T3 are two forests diagrams. The figure shows that the direct computation estimator of the correlation matrix agrees very well with the subsampling estimation. The small discrepancies come from taking into account only intraforest correlations.

The direct computation of Equation (C1) is slightly different for the cross-correlation, since there only one delta is used in each pair:

$$ C_{AB} = \frac{\sum_{(i,j) \in A} \sum_{(k,l) \in B} w_i w_j w_k w_l \langle \delta_i \delta_k \rangle}{W_A W_B}. \qquad (C3) $$

In this equation, $\delta_i$ and $\delta_k$ are in the pairs A and B, and j and l are the quasars of pairs A and B. We show in Figure C1 that the direct computation, “Diagram T1234” (see Figure A.1 of dMdB17), describes to good approximation the correlation matrix estimated via subsampling.

In the shuffling method, we shuffle the angular distribution of the lines of sight and compute for each realization the autocorrelation function of Equation (11). This allows us to have multiple measured autocorrelations, on the same footprint and with the same redshift distribution. However, since we shuffled the angular position of the line of sight, the expected value of the autocorrelation is zero, $\langle \xi_{\text{shuffle}} \rangle = 0$. The

<footnote n="44">⁴⁴ This is strictly true only for Gaussian fields, but it works well in our context because intraforest correlations are much larger than interforest correlations.</footnote>
<figure><img src="image_23.png" /><figcaption>**Figure C1.** Normalized covariance matrix, i.e., the correlation matrix (Equation (15)), for the autocorrelation, $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{Ly}\alpha(\text{Ly}\alpha)$ (top panels), and for the cross-correlation, $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{quasar}$ (bottom panels). In each panel, the blue points and curves give the subsampling correlation matrix (Equation (13)), the orange curve the direct computation (Equations (C1) and (C3)), and the green curve the shuffle estimator (Equation (13)). The two panels on the left (right) give the correlation as a function of $\Delta r_{\parallel} = |r_{\parallel}^{A} - r_{\parallel}^{B}|$ ($\Delta r_{\perp} = |r_{\perp}^{A} - r_{\perp}^{B}|$). The solid curves are for $\Delta r_{\perp} = 0 \text{ h}^{-1} \text{ Mpc}$ ($\Delta r_{\parallel} = 0 \text{ h}^{-1} \text{ Mpc}$). The dashed curves are for $\Delta r_{\perp} = 4 \text{ h}^{-1} \text{ Mpc}$ ($\Delta r_{\parallel} = 4 \text{ h}^{-1} \text{ Mpc}$).</figcaption></figure>

covariance matrix is then given by the same equation as the one defined in Equation (13), where *s* is one of the shuffle realizations and where $W_A^s = 1$. Because of the shuffling, any correlations from the 3D distribution of the lines of sight are lost; however, any correlations linked to pixels from the same line of sight are preserved.

Figure C1 shows in the top two panels, in green, the estimated correlation matrix from shuffling the lines of sight. It has more noise than other estimators, for the simple reason that we used only 100 realizations. This correlation matrix agrees very well with the one from subsampling. We observe some differences of less than 5% of correlation in the $\Delta r_{\perp}$ direction, linked to the lack of 3D correlation.

## Appendix D
### Alternative Analyses of the Data

In this appendix we present, in Tables D1 and D2, the results of alternative analyses. Table D1 and Figure D1 show the evolution of the BAO and bias parameters as components of the model are added. We note that after the addition of metals the BAO parameters are very stable. This is even true with the addition of polynomial broadband (BB) terms to simulate possible unknown components in the smooth part of the power spectrum. As can be expected, the bias parameters are less stable, changing significantly when the model includes HCDs and polynomial broadbands.

Table D2 presents the results of fits using the baseline model but under different data selection criteria or parameter choices.

With the exceptions discussed below, these alternative analyses have no significant effects on the BAO parameters. After the nominal baseline model, the table shows, first, the effect of not masking DLAs, thereby treating the intra-DLA flux–flux transmission field in the same way as the field due to the IGM. The next four entries in the table show the effect of using other quasar redshift estimators. The variations in the BAO parameters for the cross-correlation are at the level of $\le 0.5\sigma$. As discussed in Appendix B, these variations are due to the statistical fluctuations in the number of pairs in ($r_{\parallel}$, $r_{\perp}$) bins because of random fluctuations in $r_{\parallel}$ separation of quasars and pixels.

The entry labeled “dmat up 300” uses a distortion matrix calculated up to separations of 300 $\text{h}^{-1}$ Mpc, i.e., beyond the nominal maximum distance of 200 $\text{h}^{-1}$ Mpc. This allows us to show that the correlation function at large distance does not have an important impact on the distorted correlation function in the region of the fit. The entry labeled “dmat rect.” uses a distortion matrix that samples the model on an ($r_{\parallel}$, $r_{\perp}$) grid with spacings half the size of the grid used to measure the correlation functions. The next two entries change the range over which the correlation functions are fitted.

The final two entries in Table D2 give the results found by splitting the sample into two independent redshift bins. These results are the only ones in the table that differ at the level of $1\sigma$ from the baseline model, as expected because of the independence of the two samples.
<table><caption>Table D1
Evolution of the BAO and Lyα Bias Parameters with the Addition of Components of the Model, When Analyzing the Correlation Functions Measured in the Main Analysis</caption><thead><tr><th>Models</th><th>α<sub>‖</sub></th><th>α<sub>⊥</sub></th><th>b<sub>Lyα</sub></th><th>β<sub>Lyα</sub></th><th>χ<sup>2</sup><sub>min</sub>/DOF, proba</th></tr></thead><tbody><tr><td colspan="6">Lyα(Lyα) × Lyα(Lyα):</td></tr><tr><td>Lyα only</td><td>1.028 ± 0.028</td><td>1.002 ± 0.048</td><td>-0.1956 ± 0.0030</td><td>1.207 ± 0.034</td><td>2241.84/(1590 - 4), p = 0</td></tr><tr><td>+Sky</td><td>1.028 ± 0.029</td><td>1.001 ± 0.048</td><td>-0.1767 ± 0.0032</td><td>1.055 ± 0.032</td><td>1942.24/(1590-6), p = 1.3 × 10<sup>-9</sup></td></tr><tr><td>+Metals</td><td>1.043 ± 0.034</td><td>0.991 ± 0.050</td><td>-0.181 ± 0.0033</td><td>1.098 ± 0.035</td><td>1766.05/(1590 - 11), p = 6.5 × 10<sup>-4</sup></td></tr><tr><td>baseline</td><td>1.047 ± 0.034</td><td>0.98 ± 0.042</td><td>-0.2009 ± 0.0039</td><td>1.657 ± 0.088</td><td>1604.79/(1590-13), p = 0.307</td></tr><tr><td>BBphysical</td><td>1.048 ± 0.031</td><td>0.988 ± 0.04</td><td>-0.2113 ± 0.0052</td><td>1.526 ± 0.095</td><td>1488.95/(1515-25), p = 0.503</td></tr><tr><td>BB</td><td>1.046 ± 0.031</td><td>0.989 ± 0.041</td><td>-0.2382 ± 0.0074</td><td>1.669 ± 0.138</td><td>1480.7/(1515-25), p = 0.563</td></tr><tr><td colspan="6">Lyα(Lyα) × quasar:</td></tr><tr><td>Lyα only</td><td>1.052 ± 0.040</td><td>0.938 ± 0.050</td><td>-0.0942 ± 0.0051</td><td>0.604 ± 0.045</td><td>3989.12/(3180 - 5), p = 0</td></tr><tr><td>+z<sub>q</sub></td><td>1.054 ± 0.041</td><td>0.931 ± 0.053</td><td>-0.294 ± 0.012</td><td>2.51 ± 0.17</td><td>3258.81/(3180 - 6), p = 0.14</td></tr><tr><td>+Metals</td><td>1.058 ± 0.042</td><td>0.929 ± 0.053</td><td>-0.295 ± 0.012</td><td>2.48 ± 0.17</td><td>3232.82/(3180 - 10), p = 0.21</td></tr><tr><td>baseline</td><td>1.059 ± 0.032</td><td>0.932 ± 0.039</td><td>-0.2249 ± 0.0102</td><td>1.946 ± 0.142</td><td>3238.96/(3180 - 10), p = 0.193</td></tr><tr><td>BBphysical</td><td>1.058 ± 0.031</td><td>0.931 ± 0.038</td><td>-0.2306 ± 0.013</td><td>1.824 ± 0.166</td><td>3041.11/(3030 - 22), p = 0.332</td></tr><tr><td>BB</td><td>1.058 ± 0.031</td><td>0.932 ± 0.038</td><td>-0.2303 ± 0.0187</td><td>1.762 ± 0.211</td><td>3040.73/(3030 - 22), p = 0.334</td></tr><tr><td colspan="6">all combined:</td></tr><tr><td>Lyα only</td><td>1.026 ± 0.020</td><td>0.975 ± 0.031</td><td>-0.1872 ± 0.0023</td><td>1.107 ± 0.024</td><td>11666.49/(9540 - 6), p = 0</td></tr><tr><td>+Sky</td><td>1.027 ± 0.020</td><td>0.974 ± 0.031</td><td>-0.1674 ± 0.0025</td><td>0.965 ± 0.023</td><td>11113.07/(9540 - 10), p = 0</td></tr><tr><td>+z<sub>q</sub></td><td>1.027 ± 0.020</td><td>0.974 ± 0.030</td><td>-0.1892 ± 0.0027</td><td>1.190 ± 0.030</td><td>10255.42/(9540 - 11), p = 1.4 × 10<sup>-7</sup></td></tr><tr><td>+Metals</td><td>1.042 ± 0.023</td><td>0.967 ± 0.031</td><td>-0.191 ± 0.0030</td><td>1.300 ± 0.061</td><td>9999.77/(9540 - 16), p = 3.4 × 10<sup>-4</sup></td></tr><tr><td>baseline</td><td>1.045 ± 0.022</td><td>0.956 ± 0.028</td><td>-0.2014 ± 0.0032</td><td>1.669 ± 0.071</td><td>9654.56/(9540 - 19), p = 0.166</td></tr></tbody></table>

Note. Uncertainties correspond to $\Delta\chi^2 = 1$. As discussed in Section 6.1, the only significant change to the BAO parameters occurs with the addition of the metals, which, because Si(1260) absorption interferes with the BAO peak near $r_⊥ = 0$, changes $α_‖$ by $\approx 0.5σ$.
<table><caption>Table D2
Results of Nonstandard Fits as Explained in the Text</caption><thead><tr><th>Analysis</th><th>α<sub>||</sub></th><th>α<sub>⊥</sub></th><th>χ<sup>2</sup><sub>min</sub>/DOF, proba</th></tr></thead><tbody><tr><td>Lyα(Lyα) × Lyα(Lyα):</td><td></td><td></td><td></td></tr><tr><td>baseline</td><td>1.043 ± 0.035</td><td>0.986 ± 0.051</td><td>1599.55/(1590 - 14), p = 0.33</td></tr><tr><td>no SHP pairs, r<sub>||</sub> &lt; 4</td><td>1.042 ± 0.034</td><td>0.985 ± 0.051</td><td>1639.66/(1590 - 12), p = 0.14</td></tr><tr><td>no DLA mask</td><td>1.037 ± 0.034</td><td>0.978 ± 0.053</td><td>1625.09/(1590 - 14), p = 0.19</td></tr><tr><td>Z</td><td>1.049 ± 0.034</td><td>0.971 ± 0.054</td><td>1623.63/(1590 - 14), p = 0.20</td></tr><tr><td>Z_PCA(no, Lyα corr.)</td><td>1.046 ± 0.035</td><td>0.980 ± 0.055</td><td>1618.25/(1590 - 14), p = 0.22</td></tr><tr><td>Z_C IV</td><td>1.058 ± 0.036</td><td>0.967 ± 0.052</td><td>1571.92/(1590 - 14), p = 0.52</td></tr><tr><td>Z_C III</td><td>1.051 ± 0.039</td><td>0.973 ± 0.060</td><td>1628.38/(1590 - 14), p = 0.17</td></tr><tr><td>dmat up 300</td><td>1.043 ± 0.034</td><td>0.984 ± 0.051</td><td>1590.51/(1590 - 14), p = 0.39</td></tr><tr><td>dmat rect.</td><td>1.040 ± 0.035</td><td>0.986 ± 0.053</td><td>1609.73/(1590 - 14), p = 0.27</td></tr><tr><td>r<sub>min</sub> = 40</td><td>1.038 ± 0.033</td><td>0.994 ± 0.049</td><td>1490.83/(1515 - 14), p = 0.57</td></tr><tr><td>r<sub>max</sub> = 150</td><td>1.042 ± 0.035</td><td>0.988 ± 0.051</td><td>1089.86/(1097 - 14), p = 0.44</td></tr><tr><td>low z</td><td>1.024 ± 0.052</td><td>0.84 ± 0.12</td><td>1542.74/(1590 - 14), p = 0.72</td></tr><tr><td>high z</td><td>1.066 ± 0.048</td><td>0.982 ± 0.055</td><td>1624.99/(1590 - 14), p = 0.19</td></tr><tr><td>Lyα(Lyα) × quasar:</td><td></td><td></td><td></td></tr><tr><td>baseline</td><td>1.058 ± 0.042</td><td>0.929 ± 0.053</td><td>3232.82/(3180 - 10), p = 0.21</td></tr><tr><td>no DLA mask</td><td>1.052 ± 0.044</td><td>0.943 ± 0.056</td><td>3225.07/(3180 - 10), p = 0.24</td></tr><tr><td>Z</td><td>1.079 ± 0.039</td><td>0.913 ± 0.047</td><td>3187.32/(3180 - 10), p = 0.41</td></tr><tr><td>Z_PCA(no, Lyα corr.)</td><td>1.066 ± 0.046</td><td>0.930 ± 0.052</td><td>3268.92/(3180 - 10), p = 0.11</td></tr><tr><td>Z_C IV</td><td>1.062 ± 0.046</td><td>0.948 ± 0.053</td><td>3298.39/(3180 - 10), p = 0.055</td></tr><tr><td>Z_C III</td><td>1.069 ± 0.048</td><td>0.929 ± 0.064</td><td>3259.71/(3180 - 10), p = 0.13</td></tr><tr><td>dmat up 300</td><td>1.059 ± 0.042</td><td>0.927 ± 0.052</td><td>3220.68/(3180 - 10), p = 0.26</td></tr><tr><td>dmat rect.</td><td>1.059 ± 0.042</td><td>0.928 ± 0.053</td><td>3226.78/(3180 - 10), p = 0.24</td></tr><tr><td>r<sub>min</sub> = 40</td><td>1.054 ± 0.040</td><td>0.936 ± 0.052</td><td>3042.47/(3030 - 10), p = 0.38</td></tr><tr><td>r<sub>max</sub> = 150</td><td>1.054 ± 0.043</td><td>0.929 ± 0.054</td><td>2288.26/(2194 - 10), p = 0.059</td></tr><tr><td>low z</td><td>1.021 ± 0.065</td><td>0.989 ± 0.072</td><td>3191.71/(3180 - 10), p = 0.39</td></tr><tr><td>high z</td><td>1.092 ± 0.056</td><td>0.879 ± 0.064</td><td>3267.25/(3180 - 10), p = 0.11</td></tr><tr><td>all combined:</td><td></td><td></td><td></td></tr><tr><td>baseline</td><td>1.043 ± 0.022</td><td>0.960 ± 0.029</td><td>9626.45/(9540 - 20), p = 0.22</td></tr><tr><td>no DLA mask</td><td>1.038 ± 0.023</td><td>0.958 ± 0.030</td><td>9705.05/(9540 - 20), p = 0.091</td></tr><tr><td>Z</td><td>1.053 ± 0.022</td><td>0.940 ± 0.027</td><td>9730.85/(9540 - 20), p = 0.064</td></tr><tr><td>Z_PCA(no, Lyα corr.)</td><td>1.049 ± 0.024</td><td>0.954 ± 0.030</td><td>9734.68/(9540 - 20), p = 0.061</td></tr><tr><td>Z_C IV</td><td>1.047 ± 0.024</td><td>0.962 ± 0.029</td><td>9746.02/(9540-20), p = 0.052</td></tr><tr><td>Z_C III</td><td>1.048 ± 0.026</td><td>0.957 ± 0.031</td><td>9786.63/(9540-20), p = 0.027</td></tr><tr><td>dmat up 300</td><td>1.043 ± 0.022</td><td>0.958 ± 0.029</td><td>9587.65/(9540-20), p = 0.31</td></tr><tr><td>dmat rect.</td><td>1.041 ± 0.023</td><td>0.960 ± 0.030</td><td>9631.07/(9540-20), p = 0.21</td></tr><tr><td>r<sub>min</sub> = 40</td><td>1.039 ± 0.021</td><td>0.968 ± 0.028</td><td>9069.50/(9090-20), p = 0.50</td></tr><tr><td>r<sub>max</sub> = 150</td><td>1.041 ± 0.023</td><td>0.961 ± 0.030</td><td>6706.63/(6582-20), p = 0.10</td></tr><tr><td>low z</td><td>1.019 ± 0.035</td><td>0.979 ± 0.061</td><td>4751.87/(4770-18), p = 0.50</td></tr><tr><td>high z</td><td>1.076 ± 0.031</td><td>0.941 ± 0.040</td><td>4905.41/(4770-18), p = 0.059</td></tr></tbody></table>

**Note.** Uncertainties correspond to $\Delta\chi^2 = 1$. All fits use the HCD model of B17, which has one more free parameter than the HCD model used in this paper, but has no significant impact on the BAO parameters.
<figure><img src="image_24.png" /><figcaption>**Figure D1.** The $Ly\alpha(Ly\alpha) \times Ly\alpha(Ly\alpha)$ (top) and $Ly\alpha(Ly\alpha) \times$ quasar (bottom) correlation functions in four ranges of $\mu$. The different lines correspond to models with an increasing number of components as shown in Table D1: $Ly\alpha$ only (including $z_q$ smearing for $Ly\alpha(Ly\alpha) \times$ quasar), addition of metals (and the sky subtraction model for $Ly\alpha(Ly\alpha) \times Ly\alpha(Ly\alpha)$), addition of HCDs (i.e., the baseline model), and the addition of a polynomial broadband with physical priors.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Figure E1.** BAO results for the autocorrelation (left) and cross-correlation (right) computed using different methodologies: 68% and 95% CLs from the $\chi^2$ scan, computed by minimizing the $\chi^2$ over all nuisance parameters (red); 68% and 95% posterior probability from the full posterior sampling (blue); 68% and 95% CLs from the Gaussian approximation of the likelihood (green). In both frequentist measurements, the CLs were computed using the $\Delta\chi^2$ obtained from fastMC simulations (Table 7).</figcaption></figure>

## Appendix E
MCMC Sampler

The analysis presented in this work, as well as in several previous Ly$\alpha$ BAO analyses, has used a frequentist interpretation of probability. In this appendix we present an alternative analysis using a Bayesian framework and discuss the differences between the two.

The current approach, referred to here as the $\chi^2$ scan, uses the frequentist method of profiling the parameters of interest (in our case $\alpha_\parallel$ and $\alpha_\perp$). A grid is made over these parameters, and at each point in the grid the $\chi^2$ is minimized over all the other (nuisance) parameters. This is not equivalent to the Bayesian concept of marginalized likelihood, where the full likelihood is integrated over the nuisance parameters. However, in some special cases, i.e., when the likelihood is Gaussian, the profile likelihood is proportional to the marginalized one. Assuming flat priors on the BAO parameters, the $\chi^2$ scan can be used to construct a proxy for the marginalized posterior distribution of $\alpha_\parallel$ and $\alpha_\perp$, allowing us to use the $\chi^2$ scans in popular packages like CosmoMC (Lewis & Bridle 2002) or MontePython (Audren et al. 2013), where different cosmological results are combined within a Bayesian framework.

The $\chi^2$ scan can also be used to compute frequentist CLs. In order to do that, a set of Monte Carlo (MC) simulations of the correlation function is used to compute the values of $\Delta\chi^2$ that will correspond to each confidence level, as described in Appendix C of dMdB17.

In Figure E1 we compare the CLs obtained using the $\chi^2$ scan (red) with Bayesian results obtained through marginalization over the full posterior distribution (blue). This analysis assumed flat priors for all nuisance parameters and used the nested sampler Polychord (Handley et al. 2015a, 2015b).

In the same figure we also include results computed assuming that the full likelihood is Gaussian (green), where instead of computing a $\chi^2$ scan we have found the maximum of the likelihood and computed its second derivatives. The CLs were then computed using the same values of $\Delta\chi^2$ obtained from the MC simulations described above. The three methods give very similar results, which indicates that the profile likelihood works well in this case, and the $\chi^2$ scan can be interpreted as a posterior and safely used in Bayesian packages. A more detailed description and comparison of the different methods is presented in Cuceu et al. (2020).

## Appendix F
Supplementary Figures

This appendix presents figures that are referenced in the core of the analysis but that are not essential to the comprehension of our work. The two correlations presented in Figures F1 and F2 are displayed in a very similar manner to Figure 6. They display the autocorrelation in their four top panels and the cross-correlation in their four bottom panels, for four wedges of $|\mu| = |r_\parallel/r|$. Figure F1 shows these two correlation functions when using one pixel in the Ly$\beta$ spectral region, instead of all in the Ly$\alpha$ spectral region. The autocorrelation is then called $\text{Ly}\alpha(\text{Ly}\alpha) \times \text{Ly}\alpha(\text{Ly}\beta)$, and the cross-correlation is called $\text{Ly}\alpha(\text{Ly}\beta) \times \text{quasar}$. Figure F2 displays the two correlation functions, using pixels in the Ly$\alpha$ spectral region, split into two redshift bins.

In Figures F3 and F4 we present a comparison of the correlations of the data and mocks. The differences reflect the different values of the bias parameters for the data (Table 6) and mocks (Table 4).
Lyα(Lyα)×Lyα(Lyβ)

<figure><img src="image_26.png" /></figure>

Lyα(Lyβ)×QSO

<figure><img src="image_27.png" /><figcaption>**Figure F1.** Same as Figure 6, but using forest pixels in the Lyβ region: Lyα(Lyα) × Lyα(Lyβ) (top four panels) and Lyα(Lyβ) × quasar (bottom four panels). Because of the low significance of the BAO signal for this sample, the fits impose $\alpha_{\parallel} = \alpha_{\perp} = 1$.</figcaption></figure>
Lyα(Lyα)×Lyα(Lyα)

<figure><img src="image_28.png" /></figure>

Lyα(Lyα)×QSO

<figure><img src="image_29.png" /></figure>

**Figure F2.** Correlation functions in two redshift bins: Lyα(Lyα) × Lyα(Lyα) (top four panels) and Lyα(Lyα) × quasar (bottom four panels). Same as Figure 6, but with two redshift bins, split at zcut = 2.5 for the autocorrelation and at zcut = 2.57 for the cross-correlation.
Lyα(Lyα)×Lyα(Lyα)

<figure><img src="image_30.png" /></figure>

Lyα(Lyα)×QSO

<figure><img src="image_31.png" /><figcaption>**Figure F3.** Comparison of the autocorrelation and cross-correlation measured in mocks and in the data, when using the Lyα absorption in the Lyα region.</figcaption></figure>
Lyα(Lyα)×Lyα(Lyβ)

<figure><img src="image_32.png" /></figure>

Lyα(Lyβ)×QSO

<figure><img src="image_33.png" /><figcaption>**Figure F4.** Comparison of the autocorrelation and cross-correlation measured in mocks and in the data, when using the Lyα absorption in the Lyβ region.</figcaption></figure>