We present the measurement of Baryon Acoustic Oscillations BAO from the Lymanalpha Lyalpha forest of highredshift quasars with the firstyear dataset of the Dark Energy Spectroscopic Instrument DESI. Our analysis uses over 420000 Lyalpha forest spectra and their correlation with the spatial distribution of more than 700000 quasars. An essential facet of this work is the development of a new analysis methodology on a blinded dataset. We conducted rigorous tests using synthetic data to ensure the reliability of our methodology and findings before unblinding. Additionally we conducted multiple data splits to assess the consistency of the results and scrutinized various analysis approaches to confirm their robustness. For a given value of the sound horizon rd we measure the expansion at zrm eff2.33 with 2 precision Hzrm eff  239.2 pm 4.8 147.09rm Mpc rd kmsMpc. Similarly we present a 2.4 measurement of the transverse comoving distance to the same redshift DMzrm eff  5.84 pm 0.14 rd147.09rm Mpc Gpc. Together with other DESI BAO measurements at lower redshifts these results are used in a companion paper to constrain cosmological parameters.
['Astrophysics - Cosmology and Nongalactic Astrophysics']
DESI 2024 IV Baryon Acoustic Oscillations from the Lyman Alpha Forest
{'DESI 2024 IV: Baryon Acoustic Oscillations from the Lyman Alpha Forest': "- 75 , 10 2 33 , 65 , 47 43  \n```\nM. Manera , D. Margala , P. Martini , M. Maus, P. McDonald , 2 L. Medina-Varela, 54 A. Meisner , 21 J. Mena-Fern'andez , 76 R. Miquel, 77 , 10 J. Moon, 78 S. Moore, 6 J. Moustakas , 79 E. Mueller, 30 A. Mu˜noz-Guti'errez, 37 A. D. Myers, 80 S. Nadathur , 27 L. Napolitano , 80 R. Neveux, 18 J. A. Newman , 39 N. M. Nguyen , 7 J. Nie , 81 G. Niz , 36 , 11 H. E. Noriega , 12 , 37 N. Padmanabhan, 13 E. Paillas , 67 , 69 N. Palanque-Delabrouille , 9 , 2 J. Pan , 7 S. Penmetsa, 67 W. J. Percival , 67 , 68 , 69 M. M. Pieri, 82 M. Pinon , 9 C. Poppett, 2 , 42 , 43 A. Porredon , 18 , 83 , 47 F. Prada , 50 A. P'erez-Fern'andez , 37 , 78 I. P'erez-R'afols , 84 D. Rabinowitz, 13 A. Raichoor , 2 C. Ram'ırez-P'erez, 10 S. Ramirez-Solano, 37 M. Rashkovetskyi , 45 C. Ravoux , 15 , 9 , 85 M. Rezaie , 16 J. Rich, 9 A. Rocher , 51 , 9 C. Rockosi , 40 , 72 , 86 N.A. Roe, 2 A. Rosado-Marin, 87 A. J. Ross , 33 , 65 , 47 G. Rossi, 88 R. Ruggeri , 20 , 35 V. Ruhlmann-Kleider , 9 L. Samushia , 89 , 16 , 90 E. Sanchez , 38 C. Saulder , 78 E. F. Schlafly , 91 D. Schlegel, 2 M. Schubnell, 7 H. Seo , 87 R. Sharples , 92 , 6 J. Silber , 2 F. Sinigaglia , 28 , 29 A. Slosar, 93 A. Smith , 6 D. Sprayberry, 21 T. Tan , 9 G. Tarl'e , 7 S. Trusov, 71 R. Vaisakh , 66 D. Valcin , 87 F. Valdes , 21 M. Vargas-Maga˜na , 37 L. Verde , 77 , 58 M. Walther , 61 , 62 B. Wang , 94 , 95 M. S. Wang , 18 B. A. Weaver, 21 N. Weaverdyck , 2 R. H. Wechsler , 48 , 96 , 49 D. H. Weinberg , 65 , 47 M. White , 97 , 43 J. Yu , 51 Y. Yu , 41 S. Yuan , 49 C. Y'eche , 9 E. A. Zaborowski , 33 , 46 , 47 P. Zarrouk , 71 H. Zhang , 67 , 69 C. Zhao , 95 R. Zhao , 27 , 81 R. Zhou , 2 H. Zou , 81\n```  \nAffiliations are in Appendix G  \nE-mail: spokespersons@desi.lbl.gov  \nAbstract. We present the measurement of Baryon Acoustic Oscillations (BAO) from the Lymanα (Ly α ) forest of high-redshift quasars with the first-year dataset of the Dark Energy Spectroscopic Instrument (DESI). Our analysis uses over 420 000 Ly α forest spectra and their correlation with the spatial distribution of more than 700 000 quasars. An essential facet of this work is the development of a new analysis methodology on a blinded dataset. We conducted rigorous tests using synthetic data to ensure the reliability of our methodology and findings before unblinding. Additionally, we conducted multiple data splits to assess the consistency of the results and scrutinized various analysis approaches to confirm their robustness. For a given value of the sound horizon ( r d ), we measure the expansion at z eff = 2 . 33 with 2% precision, H ( z eff ) = (239 . 2 ± 4 . 8) (147 . 09 Mpc /r d ) km/s/Mpc. Similarly, we present a 2.4% measurement of the transverse comoving distance to the same redshift, D M ( z eff ) = (5 . 84 ± 0 . 14) ( r d / 147 . 09 Mpc) Gpc. Together with other DESI BAO  \nmeasurements at lower redshifts, these results are used in a companion paper to constrain cosmological parameters.", '1 Introduction': "Baryon Acoustic Oscillations (BAO) in the distribution of matter are a unique probe of the cosmic expansion history and the geometry of the Universe [1]. By themselves, BAO measurements at different redshifts enable precise measurements of the energy density parameters. In combination with studies of the Cosmic Microwave Background (CMB), they allow us to better constrain extensions to the ΛCDM model [2, 3].  \nThe Dark Energy Spectroscopic Instrument (DESI, [4]) project aims to measure BAO with unprecedented precision over a wide range of redshifts. DESI is in the midst of a five-year campaign to obtain accurate redshifts for 40 million galaxies and quasars over 14 000 square degrees. The survey started in May 2021 and the upcoming DESI Data Release 1 (DESI DR1, [5]), covering data collected until June 14, 2022, contains about 13 million galaxies, 1.5 million quasars, and 4 million stars over an area of more than 9 500 square degrees 1 . DESI has obtained seven BAO measurements at different redshifts with this DR1 dataset.  \nThe BAO measurements from the clustering of DESI galaxies and quasars at z < 2 are presented in [6]. In this publication we present a BAO measurement at z = 2 . 33 using the auto-correlation of the Lymanα (Ly α ) forest dataset from DESI DR1 and its crosscorrelation with quasar positions. Neutral hydrogen along the line-of-sight towards highredshift quasars cause absorption features in their spectra. While gas pressure dominates the distribution of gas on scales of tens of kiloparsecs [7], on larger scales the Ly α forest is a powerful tracer of the density fluctuations [8]. The cosmological interpretation of all DESI DR1 BAO measurements is presented in [9].  \nWhile the first BAO measurements [10, 11] used the distribution of galaxies to trace the density fluctuations, the Ly α forest can also be used to extend the BAO measurements  \nto higher redshifts than those available with current galaxy surveys. The Baryon Oscillation Spectroscopic Survey (BOSS, [12]) presented the first BAO measurement using the autocorrelation of the Ly α forest [13-15] with measurements of 50 000 quasar spectra from the 9th data release of the Sloan Digital Sky Survey (SDSS DR9, [16, 17]). With the 11th data release (DR11, [18]), BOSS presented as well the first BAO measurement from the crosscorrelation of quasars and the Ly α forest [19], which doubled the amount of information available from the same dataset.  \nUpdated Ly α BAO measurements were published following subsequent SDSS data releases [20-24], including data from the Extended Baryon Oscillation Spectroscopic Survey (eBOSS, [25]). The final Ly α BAO results used 210 000 Ly α forests from BOSS and eBOSS observations in the 16th data release (SDSS DR16, [26]). These results were published in du Mas des Bourboux et al. (2020) [27] and are referred to as dMdB20 in the rest of this article. dMdB20 has been the state-of-the art in Ly α BAO measurements until this publication, and we compare the measurements and discuss the methodological differences in the next sections (notably the use of a cross-covariance between the auto-correlation of the Ly α forest and its cross-correlation with quasars). In this work we use over 420 000 Ly α forests, doubling the number of lines of sight used in dMdB20 .  \nThe structure of this paper is designed to guide readers through the main aspects of the BAO measurements with the Ly α forest using DESI's first year of data. We start in section 2 with a description of the DESI survey and the datasets used in our analysis, i.e. the Ly α forest absorption and quasar catalogues [28, 29], along with the masking of Broad Absorption Lines and Damped Lymanα systems found in the spectra. The methodology to extract the Ly α forest fluctuations from the spectra of DESI quasars is presented in [30], and it is summarised in section 2.5. We present the measurement of correlations in section 3 and how we model them in section 4. The methodology employed in both of these sections is described in more detail in [31]. In section 4 we also summarize the (minor) contamination from sky residuals and other pipeline-induced systematics, which is described in detail in a companion paper [32]. In section 5 we present the main result from this work: the most precise BAO measurement at z > 2 to date.  \nIn order to minimize confirmation (and other) biases, the DESI Collaboration used blinding strategies in the BAO measurements. The analysis methodology was entirely developed on blinded measurements, and we were only allowed to unblind the measurement once we had fulfilled a long list of validation tests. The blinding strategy used in the Ly α BAO measurement is described in appendix C, and the analysis validation is presented in section 6. In [33] we present the validation tests in more detail. These tests are based on 150 synthetic realisations of our dataset (or mocks). These mocks were generated following the methods described in [34]. In section 7 we contextualize our findings, compare them to the Ly α BAO measurement from dMdB20 , and discuss a few options to improve future Ly α BAO analyses. We conclude in section 8.", '2 Data Sample': "We use quasar spectra collected during the first year of the DESI main survey. The observations were conducted with the Mayall 4-m telescope at Kitt Peak National Observatory, in Arizona, with a new prime focus, multi-fiber spectrograph. It consists of a new corrector equipped with an atmospheric dispersion compensator [35] providing a 3.2 degree diameter field of view, and a focal plane composed of 5000 robotically actuated fibers [36] that dis-  \ntribute the light to 10 spectrographs situated below the telescope, in a temperature-controlled room. Each spectrograph is composed of three arms, blue (3600-5930 ˚ A), red (5600-7720 ˚ A) and near infrared (7470-9800 ˚ A). In this analysis we use only data from the blue arms. In each blue camera, the light from 500 fibers is dispersed and refocused, forming 500 spectral traces on a 4096x4096 pixel STA4150 CCD. Despite this dense fiber packing, the crosscontamination of the spectra from adjacent fibers was measured to be only of 0.2% or smaller after spectral extraction (for wavelengths shorter than 8900 ˚ A, see [37] section 4.1). This effect is corrected in post-processing and we expect the residual cross-talk to be negligible. The CCD pixel size is 15 µ m, corresponding to 0.6 ˚ A in wavelength. The spectrograph line spread function is 1.8 ˚ A FWHM in the blue channel, which corresponds to a resolution from 2000 to 3400 depending on the wavelength. Many more details on the instrument can be found in [38].  \nThe main survey started on May 14, 2021 after a survey validation period (see [39] and references therein) used to tune the target selection, the fiber assignment, the exposure times, and exercise the data processing pipeline that is run every day to validate the observations. DESI is running two programs for bright and dark time, switching dynamically between them with the observation conditions (moonlight, but also seeing, sky brightness and sky transparency). The choice of pointing (called a tile when combined with a specific selection of targets to observe) and the allocation of fibers to targets is performed automatically during the night [40, 41]. The exposure times are adjusted in real time in order to obtain the same signal to noise for a fiducial target in each pointing. The quasar sample is observed during dark time. The quasar targets that are spectroscopically confirmed and have a redshift > 2 . 1 are scheduled for re-observation to build up the signal to noise necessary for the Ly α forest measurements. The objective is to acquire four exposures for each quasar, each with an effective exposure time 2 of 1000 sec. We use data collected until June 14, 2022 for this analysis which will be made public in the DESI Data Release 1 (DESI DR1) [5].  \nData collected on the mountain are transferred to NERSC (National Energy Research Scientific Computing Center) within minutes and processed by the offline pipeline [37]. CCD images are preprocessed, the spectra extracted and calibrated for each exposure. It is worth noting that the extraction algorithm returns spectra with uncorrelated noise on the same wavelength grid of bin 0.8 ˚ A for all fibers and exposures, which simplifies the co-addition (or averaging) of spectra from the same target obtained in different exposures. The data processing pipeline provides 2 sets of co-added spectra, one per spectrograph and tile, combining data of several exposures and nights for the same pointing and fiber allocation, and one per HEALPix pixel [42] on the sky where the full co-added spectra (across exposures and tiles) of all the targets in a HEALPix pixel are saved. We use this later data set for this analysis. These spectra include fluxes, a wavelength array (in vacuum, and in the solar system barycenter frame), an estimate of the flux variance, a mask to flag bad pixels, and a resolution matrix that encodes the line spread function of the instrument [37].  \nIn the subsequent sections, we delineate the processing applied to our quasar sample, both at the catalog and spectral level, including the identification of Damped Lymanα (DLAs) systems and Broad Absorption Line (BAL) troughs, to yield the ultimate datasets for our analysis. While a comprehensive overview is presented in [30], here we emphasize on the small but yet significant adjustments made to enhance our methodology.  \nTable 1 . Total number of quasars in DESI DR1 that pass our selection criteria. Those with z > 1 . 77, referred to as 'tracers', contribute to the measurement of the Ly α × QSO cross-correlation. The bottom part of the table shows the statistics of number of lines of sight available in each of the three rest-frame wavelength regions (or forests) discussed in section 2.4. Some of the forests are discarded because they do not have 150 valid pixels (too short), or because they do not have a valid continuum fit (negative continuum). These cuts are described in section 2.5.", '2.1 The quasar catalog': "Studies of the Ly α forest require a pure sample of quasars with accurately determined redshifts. To meet these stringent criteria, DESI employs a multi-step classification and redshift fitting procedure for the identification of quasar spectra. This process is extensively outlined in [28] (see their figure 9 summarizing the workflow) and complementary information is also provided in [5, 29, 43]. The cornerstone of DESI's classification and redshift fitting is Redrock , a template fitting code [44]. It utilizes Principal Component Analysis (PCA) templates representing three broad object classes (stars, galaxies, and quasars) while scanning a range of redshifts. The determination of the most probable spectral class and redshift for a given spectrum relies on the lowest χ 2 fit. Additionally, DESI employs two quasar-specific classifiers, referred to as 'afterburners', to enhance the completeness of the quasar sample by ∼ 10% [28, 29]. One afterburner scrutinizes spectra not classified as quasars by Redrock for broad MgII emission; upon detection, it reclassifies the spectrum as a quasar without altering the assigned redshift. The second afterburner, known as QuasarNET [45, 46], employs a deep convolutional neural network to identify potential quasar emission features and estimate the redshift. If QuasarNET identifies a target as a quasar, there is another Redrock fitting iteration, although this time with only quasar templates and a redshift prior based on the QuasarNET result.  \nThe resulting DESI DR1 quasar catalog was produced by the three classifiers, and is slated for publication as part of first DESI Data Release (DR1) [5]. This catalog boasts an estimated completeness and redshift purity that surpasses 95% and 98%, respectively [47].  \nFrom this catalog we keep only quasars with ZWARN=0 or ZWARN=4 , rejecting quasars with any issue reported by the spectroscopic pipeline except for the low ∆ χ 2 flag of Redrock that is irrelevant for QSOs. We also use updated redshifts for this Ly α analysis. These are based on new quasar templates tailored for high redshifts, specifically 1 . 4 < z < 7 . 0. Additionally, a slightly modified version of the Redrock code was used that integrates a more recent optical depth model from [48]. These subtle adjustments were prompted by the identification of a bias for redshifts z > 2, as identified by their impact on the Ly α quasar cross-correlation (see [49] for more details). The statistical precision of the updated redshifts is estimated to be better than 150 km s -1 , with a catastrophic redshift failure rate of ∼ 2.5% ( ∼ 4%) for redshifts in error by more than 3000 km s -1 (1000 km s -1 ).  \nIn the left panel of figure 1 we show the spatial distribution of quasars in the DESI DR1 sample (green points), compared to the SDSS DR16 footprint dMdB20 (red curve), and the  \nFigure 1 . Left: Expected final DESI (blue) and SDSS-DR16 footprint (red) together with the spatial distribution of DESI DR1 observed quasars (green). For reference we also show the Galactic plane (solid black) and the Ecliptic plane (dotted black). Right: Number of observations for Ly α quasars in the DESI DR1 sample.  \n<!-- image -->  \nexpected final footprint DESI (blue curve), while the right panel shows the distribution of number of observations for Ly α quasars in DESI DR1. The number of quasars in the DESI DR1 sample is given in table 1, together with the number of quasar spectra covering each of the rest-frame wavelength 'regions' described in section 2.4.  \nThe redshift distribution of DESI DR1 and SDSS DR16 quasars are compared in the left panel of figure 2, together with the redshift distribution of Ly α pixels. In the right panel of the same figure we show the contribution of different redshifts to the measurement of the four correlation functions discussed in section 3. Integrating the curve, one finds that 95% of the measurement of the Ly α × Ly α auto-correlation comes from 1 . 96 < z < 2 . 8, while 95% of the Ly α × QSO cross-correlation comes from 1 . 96 < z < 2 . 95. As discussed in section 3, however, quasars at redshifts as low as z = 1 . 77 and as high as z = 4 . 16 can also contribute to the measurement of the cross-correlation.", '2.2 Damped Ly α systems': "The Ly α forest BAO measurement exploits the fact that typical fluctuations in the forest trace matter fluctuations of moderate over density, or under density, along the line of sight. DLA systems are produced by neutral hydrogen concentrations with column density, N HI > 2 × 10 20 cm -2 , typically arising in the interstellar medium of high redshift galaxies. Because of the damping wings, each DLA can affect a noticeable fraction of a Ly α forest spectrum (corresponding to thousands of km/s), so even though these systems are rare, it is important to mask the contaminated regions of the spectra, in order to make a robust BAO measurement with uncertainties that can be easily modeled. In DESI we have used two DLA finders, one based on a convolutional neural network (CNN), an algorithm first proposed in [50], and another one based on a Gaussian Process (GP) model [51]. Their adaptation to DESI and performance in DESI's simulated spectra is presented in [52]. The performance of both algorithms in the DESI DR1 dataset will be presented in an upcoming paper.  \nFor the purpose of this work, we construct a DLA catalog based on the combined result of the two DLA finders. We select DLAs found by the two algorithms with a threshold probability of 50%, as done in [30, 31]. However, in order to ensure the purity of the DLA catalog we also limit it to DLAs detected in spectra with signal-to-noise (SNR) larger than  \nFigure 2 . Left: Redshift distribution of quasars in DESI DR1 (orange), compared to the distribution in SDSS DR16 (green) and to the distribution of Ly α pixels (blue, divided by 200). Right: Contribution of different redshifts to the four measured correlation functions. In particular, we show the sum of weights used in equations (3.5) and (3.6) as a function of redshift, for large transverse separations (to reduce the biasing effect from the clustering of background quasars) and r ∥ = 0 (the contribution varies as a function of r ∥ , especially for quasar cross-correlations). As described in section 3, we measure Ly α correlations in two different rest-frame wavelength regions of the quasar spectra (A and B).  \n<!-- image -->  \n3 (mean SNR within the Ly α forest region). The DLAs included in the final catalog are masked in the spectra before computing the forest fluctuations (see section 2.4) and the contamination from undetected ones is taken into account in the modelling of the correlation functions (see section 4.6).", '2.3 Broad Absorption Line quasars': 'Numerous studies of large quasar samples have identified broad absorption line (BAL) troughs in 10 - 30% of quasars [e.g., 53, 54]. While these features are most commonly seen as blue shifted absorption relative to the C IV emission feature, the absorption is also associated with many other features, including several that contribute absorption in the Ly α forest region [55]. Significant BAL features associated with C IV and other strong emission lines are also expected to introduce some systematic redshift biases and increase redshift errors for a subset of BAL quasars [56]. We consequently need to identify and characterize BALs in the DESI quasar sample to mitigate their impact on the cosmological analysis.  \nWe catalog BALs in the DESI quasar sample following the same approach employed by [57] for the DESI Early Data Release (EDR). This approach fits a series of continuum components to each quasar, and iteratively identifies and masks regions that meet the Absorptivity Index [AI, 58] and Balnicity Index [BI, 59] criteria. For the DR1 BAL catalog, we use the same components developed by [60]. The catalog includes measurements of the velocity range of each trough, AI and BI values, and other quantities measured relative to the C IV and Si IV features, similar to the DESI EDR catalog. This catalog is used in our analysis to mask the corresponding locations of BAL features before computing the Ly α forest fluctuations (see section 2.4).', '2.4 Catalog of the Ly α forest': "We discuss here the subset of quasar spectra that we use to study the fluctuations in the Ly α forest. We consider two different catalogues that correspond to Ly α absorption in two  \nFigure 3 . QSO spectrum from the first year of DESI data at redshift z = 3 . 14 (TargetID = 39627581225438176). The region B is highlighted in purple. The region A is highlighted in indigo. The C IV and C III regions are highlighted in various shades of green. While there is almost no CIII absorption, the C IV absorption spans leftward of the C IV doublet, contaminating the Ly α regions A and B. The Ly α absorption extends into region B. For better visualization, we have chosen a relatively high signal-to-noise spectrum and we have smoothed it by averaging the pixels in groups of 25.  \n<!-- image -->  \ndifferent regions of the quasar spectra. Region 'A' is defined as the rest-frame wavelength range from 1040 to 1205 ˚ A 3 , while region 'B' covers the range from 920 to 1020 ˚ A. Pixels in the B region are also affected by absorption from other Lyman lines, since it extends beyond the Ly β , Ly γ and Ly δ emission lines, and therefore it require special treatment 4 . See figure 3 for an illustration of these regions in a quasar observed with DESI.  \nIn addition to the rest-frame wavelength cut, we also restrict our observed wavelength range from 3600 ˚ A to 5772 ˚ A. The lower bound is the minimum wavelength provided by the pipeline; the upper bound corresponds to the middle of the overlap region between the blue and the red arms of the spectrographs 5 . These cuts in observer-frame and rest-frame wavelength limit the redshift range of the background quasars from z = 2 . 1 to z = 4 . 4 for the region A and z = 2 . 6 to z = 5 . 1 for the region B. The number of quasars available for these catalogues is given in table 1.  \nOnce we have the initial set of Ly α forests, we apply four different masks to remove bad pixels and astrophysical contaminants. The first mask comes from the data reduction pipeline, which flags bad pixels in the spectra [37]. These are usually caused by cosmic rays or defects in the CCD cameras. As shown in [30], the fraction of pixels masked as a function of observed wavelength is generally constant and below 1%. The second mask consists of a set of narrow windows in observed wavelength: 3933 . 0 to 3935 . 8 ˚ A, 3967 . 3 to 3971 . 0 ˚ A, and 5570 . 5 to 5586 . 5 ˚ A (see [30]). This mask aims to remove galactic absorption not observed in the calibration stars (first two lines) and residual sky lines from the sky subtraction models  \n[37].  \nThe other two masks aim at removing the effects of astrophysical contaminants. In particular, we mask DLAs (see section 2.2) and BALs (see section 2.3). For DLAs, we mask the regions of the Ly α forest where the presence of a DLA decreases the transmitted flux by 20% or more. We then correct the remaining DLA wings using a Voigt profile. For BALs, we follow the procedure described in [61], which we summarize here. We mask the expected locations of all potential BAL features associated with, irrespective of whether or not the absorption is apparent. These include Ly α , N IV, C III, Si IV, and P V in region A, and O VI, OI, Ly β , Ly γ , N III and Ly δ in region B. Because of this, the mask may remove some path length that is unaffected by BALs. However, we note that there is a net increase in path length compared to completely discarding these spectra as done in previous analyses [e.g. 27]. The paper by [62] investigates the impact of other masking strategies on the BAO analysis.  \nAt this stage, we discard forests that are too short (defined as having less than 150 valid pixels corresponding to a length of 120 ˚ A). This is necessary as we later fit the unabsorbed quasar continua to extract the transmitted flux field (see section 2.5), and having forests that are too short interferes with the continuum fitting procedure. The number of forests lost because of this in each of the samples is given in table 1.  \nFinally, to check for unaccounted-for calibration residuals, we compute the mean transmitted flux fraction in a quasar rest-frame region where there is no Ly α absorption. As in [30], we use the C III region ( λ rf ∈ [1600 , 1850] ˚ A). We note that previous analysis from eBOSS used longer rest-frame wavelengths, but [30] showed that the measured calibration residuals are the same irrespective of the metal region used to measure them, and there are more spectra available in the C III region. We use this small correction (less than 5% in the relevant wavelengths range) to re-calibrate our fluxes and instrumental noise estimates.", '2.5 Continuum fitting': 'The initial step in our data analysis is to compute the transmitted flux field δ q ( λ ). In general terms, it is defined as the ratio of the observed flux to the expected flux:  \nδ q ( λ ) = f q ( λ ) F ( λ ) C q ( λ ) -1 , (2.1)  \nwhere F ( λ ) is the mean transmission and C q is the unabsorbed quasar continuum. The sub-index q indicates that these are line-of-sight (quasar) dependent. We label the process of estimating their product as continuum fitting .  \nThe continuum fitting procedure we use is described in detail in [30] and we summarize it here. The product FC q ( λ ) is taken to be a universal function of the quasar C ( λ rf ), corrected by a first degree polynomial in Λ ≡ log λ to account for quasar diversity and the redshift evolution of ¯ F ( z ):  \nF ( λ ) C q ( λ ) = C ( λ rf ) ( a q + b q Λ -Λ min Λ max -Λ min ) . (2.2)  \nHere, the universal function C ( λ rf ) is computed to be the weighted mean (inverse variance weighting) of all the lines of sight entering the analysis. The parameters a q and b q are estimated by maximizing the likelihood function  \n2 ln L = -∑ i ( f i -FC q ( λ i , a q , b q ) ) 2 σ 2 q ( λ i , a q , b q ) -∑ i ln σ 2 q ( λ i , a q , b q ) , (2.3)  \nwhere σ 2 q is the total variance of the data, including the contribution from the noise estimated from the pipeline ( σ 2 pip ,q ), modulated by a function η pip ( λ ) to account for small inaccuracies in the noise estimation, and the intrinsic variance of the Ly α forest, σ 2 LSS , multiplied by ( FC q ) 2 :  \nσ 2 q ( λ ) = η pip ( λ ) σ 2 pip ,q ( λ ) + σ 2 LSS ( λ ) ( FC q ) 2 ( λ ) . (2.4)  \nWe keep the original wavelength bin size of 0 . 8 ˚ A when computing the δ q ( λ ) following the study of [30]. We note that the continuum fitting procedure distorts the δ q field, since line-ofsight fluctuations on scales comparable to the length of a given forest will be suppressed when fitting the ( a q , b q ) parameters. This is taken into account with the definition of a projected flux transmission field in section 3.1 and in the correlation function model in section 4.3. We perform the continuum fit in each of the analysed regions independently. This means that we have up to 3 independent sets of ( a q , b q ) per quasar as each quasar can be used in the Ly α regions A and B, and in the calibration region.  \nOccasionally, the best-fit values of a q and b q , combined with the shape of C ( λ rf ) result in a negative estimated product FC q ( λ ). Whenever this happens for at least one pixel, we deem the continuum fit as problematic and discard the entire region. The number of forests lost because of this in each of the samples is given in table 1. Note that, because the a q and b q are fitted independently in each of the regions of interest, a particular spectrum might be removed from one of the regions but kept in another.', '3 Measurement of correlations': 'In section 2 we have described the three cosmological datasets that we use in our BAO measurement: the catalog of quasar positions (angles and redshifts) and the catalogs of Ly α fluctuations in the rest-frame regions A and B. There are six (3x2) possible correlations that we could use, but following dMdB20 we focus on the four correlations that have higher signal to noise: we ignore the auto-correlation of quasars, and the auto-correlation of Ly α fluctuations in the B region. For convenience, in some sections we group the correlations in two subsets: we use the term auto-correlation to refer to the combination of the auto-correlation of the Ly α fluctuations in region A, Ly α (A) × Ly α (A), and the correlation of Ly α fluctuations in region A with the Ly α fluctuations in region B, Ly α (A) × Ly α (B) 6 ; we use the term crosscorrelation to refer to the correlation of quasar positions with the Ly α fluctuations in region A, Ly α (A) × QSO, and in region B, Ly α (B) × QSO.  \nIn this section we summarise the methodology used to measure the different correlations and their covariance matrix, and we refer the reader to [31] for more details. We measure the correlations with the software picca 7 [73] that was developed originally for the Ly α analysis of BOSS and eBOSS data, and that was also used to estimate the Ly α fluctuations. The only relevant difference with respect to [31] is that we now include the cross-covariance of the different correlations, which were considered independent in previous analyses. We discuss this in section 3.3.  \nTable 2 . Parameters of the flat-ΛCDM cosmological model used in the analysis, both to compute comoving separations and in the modelling (described in section 4). The first part of the table gives the cosmological parameters, the second part gives derived quantities used in this paper.  \nWe measure the correlations in bins of comoving separation along ( r ∥ ) and across ( r ⊥ ) the line of sight, computed from the angular and redshift separations using a fiducial cosmology based on the best-fit flat ΛCDM model from Planck 2018 [2] (see table 2). We present a single measurement of the correlations averaged over a wide range of redshifts (see figure 2). This is motivated by the fact that the BAO scale in comoving coordinates does not vary much with redshift if the true cosmic expansion is not very far from the fiducial one for those redshifts 8 , and that considering several redshift bins would make the BAO peak less significant and possibly more difficult to fit [74]. As discussed in section 4, we model the correlations at an effective redshift ( z eff = 2 . 33) and report the BAO measurement at that redshift. In our fiducial model and at our effective redshift, the BAO scale corresponds to an angular separation of 1 . 46 degrees and an observed wavelength separation of 141 ˚ A.', '3.1 Measurement of the auto-correlation': 'The correlation function measurement relies on a fiducial cosmology to convert the angular separations and redshifts into comoving separations. We use the same approach as in [31] and earlier studies. We first define the longitudinal and transverse separations ( r ∥ , r ⊥ ) for a pair of measurements ( i, j ) at redshifts ( z i , z j ) and angular separation θ ij as  \nr ∥ = [ D C ( z i ) -D C ( z j )] cos ( θ ij / 2) , r ⊥ = [ D M ( z i ) + D M ( z j )] sin ( θ ij / 2) , (3.1)  \nwhere D C ( z ) is the comoving distance and D M ( z ) the angular (or transverse) comoving distance. For our fiducial cosmology with Ω k = 0, they are identical.  \nThe estimator for the correlation function of the transmitted flux field is a simple weighted average. For this purpose, we define the following weight for the flux decrement δ q ( λ ),  \nw q ( λ ) = ( 1 + z λ 1 + z 0 ) γ α -1 [ η pip ( λ ) ( σ pip ,q ( λ ) FC q ( λ ) ) 2 + η LSS σ 2 LSS ( λ ) ] -1 . (3.2)  \nThe right hand term is the inverse of the variance of the flux decrement (equation (2.4)), but with an additional scaling term, η LSS , that increases the contribution of the large scale structure variance compared to the pipeline noise term. This correction was introduced by [30] to improve the precision of the correlation function measurement 9 . Our estimator is sub-optimal because it does not take into account the correlations between neighboring flux decrement values caused by the large scale structure. The adjustment of weights with η LSS is a simple way to avoid over-weighting pixels that have a high signal to noise ratio but correlated information content. We use a value of η LSS = 7 . 5, which was found to minimise the covariance matrix of the Ly α auto-correlation given our wavelength bin size of 0.8 ˚ A [30]. The left hand term is a scaling factor to account for the variation of the amplitude of the correlation function with redshift. It is the optimal term for the measurement of the shape of the correlation function. We set the bias evolution index γ α = 2 . 9 as in previous studies, following [75] (the value of z 0 does not affect the results as it cancels in the estimator of the correlation function).  \nIn order to estimate the Ly α fluctuations, we have divided the observed flux by a model of the continuum of each quasar (see equation (2.1)). As noted in section 2.5, the parameters ( a q , b q ) of this model have been fitted using all the Ly α pixels in the spectrum and this makes our estimated fluctuation at a given pixel ( δ i ) dependent on the flux of the other pixels in the spectrum, and therefore on their true Ly α fluctuation ( δ t j ) [8, 76]. Following earlier work [21, 22], we linearise this relation and approximate the distortion 10 as δ i = η cont ij δ t j , with η cont defined as  \nη cont ij = δ K ij -w cont j ∑ k w cont k -w cont j Λ i Λ j ∑ k w cont k Λ 2 k , (3.3)  \nwhere δ K ij is the Kronecker delta function, w cont q = ( FC q ) 2 /σ 2 q are the weights used in continuum fitting (equation (2.3)) and Λ i ≡ log( λ i ) -∑ k w cont k log( λ k ) / ∑ k w cont k .  \nIt is clear from equation (3.3) that during continuum fitting the weighted mean and the linear trend over each Ly α forest are set to zero. In order to make sure that these are also zero when using the weights w defined in equation (3.2), we apply a similar linear projection to the measured fluctuations δ j and define a projected field ˜ δ i = η ij δ j , where η is equivalent to η cont but using the weights w instead of w cont in equation (3.3). One can show that, to first order in δ t , η is also the relation between the projected field and the true fluctuation 11 :  \n˜ δ i = η ij δ j = η ij η cont jk δ t k = η ij δ t j . (3.4)  \nAs a result we do not need to know the exact values of η cont ij , but only those of η ij . In section 4.3 we will use this linear projection operator η to model the correlations of the projected field ˜ δ q ( λ ).  \nThe correlation function is measured in bins of 4 h -1 Mpc ranging from 0 to 200 h -1 Mpc along both r ∥ and r ⊥ for a total of 2 500 bins. Calling M a two-dimensional bin of comoving separation, the estimator of the correlation function averaged in the separation bin M is the weighted average:  \nξ M = ∑ ( i,j ) ∈ M w i w j ˜ δ i ˜ δ j / ∑ ( i,j ) ∈ M w i w j (3.5)  \nwhere ( i, j ) ∈ M are all the pairs of projected transmitted flux measurements ˜ δ i and ˜ δ j from quasar lines of sight separated by an angle θ ij and at redshift (or wavelength) z i and z j , such that their comoving separation is in M . We only consider pairs from different quasar spectra to avoid the contribution of correlated residuals within a spectrum caused by continuum fitting errors.  \nIn figure 4 we show the measurement of the Ly α forest auto-correlation when using pixels from region A (Ly α (A) × Ly α (A), top panel) and when correlating pixels from region A with pixels in region B (Ly α (A) × Ly α (B), bottom panel). We show the measurement as a function of the total separation r = ( r 2 ∥ + r 2 ⊥ ) 1 / 2 in wedges defined by a range of the cosine between the orientation of the pair and the line-of-sight, µ = r ∥ /r . This is obtained by resampling the original data, which has the form of a rectangular 2D grid of ( r ∥ , r ⊥ ) bins, into ( r, µ ) bins. We evaluate the uncertainties in those bins using the covariance matrix, the computation of which is described in section 3.3. The resampling results in additional correlation between neighboring ( r, µ ) bins that share fractions of the original bins. We note that this resampling into wedges is performed for display purpose only and is not used for the fit, which is realized on the original rectangular grid.', '3.2 Measurement of the cross-correlation': 'With the same notations as in the previous section, we use the following estimator [76] for the quasar Lymanα cross-correlation in a separation bin M :  \nξ M = ∑ ( i,j ) ∈ M w i w Q j ˜ δ i / ∑ ( i,j ) ∈ M w i w Q j . (3.6)  \nHere ( i, j ) ∈ M stands for pairs made of a projected transmitted flux measurement ˜ δ i at a redshift z i along a quasar line of sight, and another quasar j at another redshift z j separated by an angle θ ij from the first one, such that their comoving separation is in the bin M (once the redshifts and θ ij are converted to comoving separation using equation 3.1). We use for quasars the weights  \nw Q j = [(1 + z Q ) / (1 + z 0 )] γ Q -1 (3.7)  \nwhere z Q is the quasar redshift, and we choose γ Q = 1 . 44 which follows closely the measured bias evolution of quasars [77].  \nWe use the same bin size of 4 h -1 Mpc for the cross-correlation but we differentiate positive and negative longitudinal separation, giving us a range from -200 to +200 h -1 Mpc for r ∥ because the cross-correlation is asymmetric. We define the sign of r ∥ such that r ∥ < 0 for pairs where the quasar is behind the Ly α transmitted flux measurement. We have 5 000 cross-correlation bins.  \nLy (A)×Ly (A)  \nFigure 4 . Measured Ly α auto-correlation when using pixels from region A (top, colored markers) and when correlating pixels from region A with pixels from region B (bottom), along with the best fit model (solid black curves), described in section 4. The different colors and markers correspond to different orientations with respect to the line-of-sight, with blue correlations being close to the line-ofsight 0 . 95 < µ < 1. The dotted curves show the best fit model with additive polynomial corrections (see section 6.4.4).  \n<!-- image -->  \nLy (A)×QSO  \nFigure 5 . Measured Ly α × QSO cross-correlation functions in region A (top, colored markers) and region B (bottom) along with the best fit model (solid black curves), described in section 4. The different colors and markers correspond to different orientations with respect to the line-of-sight, with blue correlations being close to the line-of-sight 0 . 95 < µ < 1. The dotted curves show the best fit model with additive polynomial corrections (see section 6.4.4).  \n<!-- image -->  \nIn figure 5 we show the measurement of the cross-correlation of quasars with Ly α pixels in region A (Ly α (A) × QSO, top panel) and with Ly α pixels in region B (Ly α (B) × QSO, bottom panel).', '3.3 Covariance matrix': 'We use the method described in more detail in dMdB20 to compute the covariances. In brief, the correlation function is first measured independently in sub-samples defined by HEALpix pixels on the sky. Each sub-sample correlation is saved with its weights W M in each bin M , which are denominators in equations (3.5) and (3.6). We do not lose or double-count pairs because each possible pair is assigned a unique sub-sample. The combined correlation function is simply the weighted mean of the sub-sample correlations, and its covariance is determined by replacing the unknown covariance of each sub-sample by the square of its difference with the mean. We ignore the cross-covariance between sub-samples, which is negligible for our scales of interest given the size of a HEALpix pixel of about (250 h -1 Mpc) 2 at z ∼ 2 . 3 for our choice of NSIDE =16. The correlation is measured in 1028 pixels, with more than 400 valid forests per pixel on average. In figure 11 we show that using NSIDE =32 instead has a negligible impact on the BAO parameters.  \nThis noisy estimate of the covariance C is then smoothed. We replace all the nondiagonal elements of the correlation matrix Corr MN ≡ C MN / ( C MM C NN ) 1 / 2 in which indices correspond to the same differences | r ∥ ( M ) -r ∥ ( N ) | and | r ⊥ ( M ) -r ⊥ ( N ) | by their average. This method has proven to be a good approximation of the covariance when compared with other methods, like a Gaussian covariance computed with the Wick expansion (see [20] and Appendix C of dMdB20 ), and it is discussed in detail in the companion paper [33].  \nWe measure with this sub-sampling technique the covariance of the full data set composed of the four correlation functions, Ly α (A) × Ly α (A), Ly α (A) × Ly α (B), Ly α (A) × QSO, and Ly α (B) × QSO. Previous works (including dMdB20 ) measured the covariances of each correlation function in isolation, ignoring the cross-covariances between them. Indeed, we show in section 6.2.2 that ignoring the cross-covariance between the correlation functions results in a 10% underestimate of uncertainties on the BAO parameters. The new correlation matrix is shown in figure 6. This full covariance estimation represents the most important change in the methodology from previous analyses of the 3D correlations in the Ly α forest.', '4 Modelling of correlations': 'In this section we describe how the correlations are modeled. The approach is similar to the one used in previous BOSS and eBOSS analyses (see dMdB20 for the latest results from eBOSS) except for few methodological changes that are summarised below, and that we describe in more detail in companion papers [31-33]. We now use the software Vega 12 , which is an improved version of the picca-fitter2 software used in these previous Ly α BAO analyses.  \nFrom the point of view of this analysis, our model has 2 important (BAO) parameters and 15 nuisance parameters that we marginalize over. We start in sections 4.1 to 4.3 with a simplified version of the model that will allow us to introduce some of the key aspects of a Ly α BAO analysis, before describing the modelling of contaminants in sections 4.4 to 4.8. The priors, best-fit values and uncertainties for the nuisance parameters are presented in appendix B.  \nFigure 6 . Global correlation matrix for the four 2pt functions included in this analysis, as measured from the scatter between correlations measured in more than 1000 HEALPix pixels, after applying the smoothing discussed in the main text. The first block of 2 500 × 2 500 in the top left corresponds to the correlation matrix of the Ly α (A) × Ly α (A) measurement, while the second block of the same size corresponds to the Ly α (A) × Ly α (B) one. The third, larger block of 5 000 × 5 000 corresponds to the Ly α (A) × QSO cross-correlation, with positive and negative values of r ∥ , and the last block (bottom right) is the correlation matrix of the Ly α (B) × QSO measurement. It is clear in that cosmic variance (in the form of off-diagonal stripes in red and blue) is relevant, even thought it is only detected at the sub-percent level. ote that we only have strong correlations ( > 10%) for bins with the same value of r ⊥ and neighbouring value of r ∥ , and that we limit the range of the color scale to 1% for visualization purposes.  \n<!-- image -->', '4.1 Isolating the BAO information': 'We start by building a model for the (anisotropic) large-scale power spectrum of fluctuations in the Ly α forest, based on linear perturbation theory. We use the linear power spectrum of matter fluctuations for our fiducial cosmology (see table 2) evaluated at our effective redshift (see below), and linear bias ( b F ) and redshift-space distortion ( β F ) parameters. When modelling the cross-correlation with quasars, we introduce another linear bias parameter ( b Q ) and we follow [78] to model the linear redshift-space distortions of quasars with β Q = f/b Q , where f is the logarithmic growth rate 13 .  \nWe then compute its inverse Fourier transform to obtain a model for the correlation function ξ ( r ⊥ , r ∥ ) that we can compare to our measurement. In our fiducial cosmology, there is an excess correlation (the BAO peak 14 ) at around 100 h -1 Mpc. We then introduce two scaling parameters, α ⊥ and α ∥ , which multiply r ⊥ and r ∥ respectively in our model, and we vary them in order to better match the BAO peak seen in the measured correlations.  \nIn order to make sure that we only extract BAO information from the fits, we decompose the model of the correlations into a peak and a smooth component following [15], and the scaling parameters ( α ⊥ , α ∥ ) are only applied to the peak component 15 . Following [15] again, we apply a Gaussian smoothing to the peak component in order to model the non-linear broadening of the BAO feature caused by non-linear growth of structure.', '4.2 Redshift evolution': 'When computing the binned correlations in equations (3.5) and (3.6), we use the same weights to compute the mean separations ( r ⊥ , r ∥ ) and to compute the mean redshift ( z ) of each bin. We evaluate the model at these coordinates. The redshift evolution of the model is captured by the linear growth of the matter power spectrum ( σ 8 ( z )), the logarithmic growth rate ( f ( z )), and the redshift evolution of the bias parameters. We define our parameters at an effective redshift ( z eff = 2 . 33, see below), and model the redshift evolution of the bias parameters with a power law, b ( z ) = b ( z eff ) [(1 + z ) / (1 + z eff )] γ . Following [27, 31], we use γ Q = 1 . 44 to describe the redshift evolution of the quasar bias, γ α = 2 . 9 for bias of the Ly α forest, and we assume that the RSD parameter of the Ly α forest ( β α ) does not vary with redshift.  \nIn order to estimate the effective redshift of our BAO measurement, we compute the mean of the redshifts of each correlation bin in the range 80 h -1 Mpc < r < 120 h -1 Mpc, weighted with their inverse variance. We do this separately for the auto-correlation ( z eff = 2 . 339) and for the cross-correlation ( z eff = 2 . 325), and we then compute a simple mean of those two values for the combined BAO measurement (rounded to two decimal values) to obtain the effective redshift of our BAO measurement ( z eff = 2 . 33).', '4.3 Distortion matrix': "As discussed in section 3, the distortions introduced with the continuum fitting have led us to the use of the projected field ˜ δ defined by equation (3.3). We must therefore fit the correlations of the projected field ˜ δ with a model that has suffered the same projection. This is performed with the 'distortion matrix' formalism introduced in [21], ξ M = ∑ N D MN ξ ' N , where ξ M refers to a ( r ⊥ , r ∥ ) bin of the distorted model and ξ ' N to a bin of the undistorted model. The matrices D MN are constructed using the same linear operators η used to compute the projected field (equation (3.3)). The elements of the matrices are given by equations (21) and (22) of dMdB20 for, respectively, the auto- and cross-correlations.  \nTwo improvements have been made to the distortion matrix treatment over that of dMdB20 . First, the undistorted model bins are now calculated on a grid of 2 h -1 Mpc (instead of the 4 h -1 Mpc bins used to measure the correlations). Second, in the r ∥ directions, we extend the modelling of the undistorted correlations to 300 h -1 Mpc rather than to 200 h -1 Mpc. This improves the accuracy of the distortion calculation at high r ∥ , but has a negligible impact on the BAO results (see section 6.4).  \nThe computation of the distortion matrix D MN is computationally intensive, and following previous work we only use a small fraction of the dataset to approximate it [21, 31]. By default we use 1% of the Ly α pixels, but in section 6.4 we show that doubling that number does not affect the BAO results.", '4.4 Metal contamination': "The Ly α forest is contaminated by absorption from atomic transitions of elements other than neutral hydrogen. The auto-correlation of those absorbers and their cross-correlation with Ly α , QSOs, or other transitions can contribute to the measured Ly α auto-correlation and its cross-correlation with QSOs, and has to be taken into account. We call them metal absorbers. Labeling δ m their contribution to the total flux decrement δ , the Ly α forest auto-correlation will have contributions of the form  \n⟨ δδ ⟩ = ⟨ δ α δ α ⟩ + ∑ m ⟨ δ m δ m ⟩ + ∑ m ⟨ δ α δ m ⟩ + ∑ m ∑ m ' = m ⟨ δ m δ m ' ⟩ . (4.1)  \n̸  \nThe second term is the contribution of the metal auto-correlations that are present for all foreground absorbers with λ m > λ min where λ min is the minimum wavelength of the forest in the quasar rest-frame. The third term is the contribution of the cross-correlation of metal absorbers with Ly α . Only transition wavelengths close to the Ly α line will have a significant contribution for the range of longitudinal separation we are studying. The fourth term is the cross-correlation of metals, which can introduce a lot of complexity in the interpretation of the measured correlation function. Fortunately, the metal absorbers are much less abundant than neutral hydrogen and hence have much smaller absorption. We include in the fit the cross-correlation of Si II and Si III absorbers, whose cross-correlation with Ly α is also observed, but we neglect the cross-correlation of other foreground absorbers. For the Ly α quasar cross-correlation, the situation is simpler, as we have only contributions from absorbers with transition wavelengths close to the Lymanα line.  \nPrevious analyses based on measurements of the Ly α 1D correlation [21], or crosscorrelation with strong absorbers [81-83], have shown that only a few transitions had to be considered for the cross-correlation terms cited above. Those are Si II lines at 1190 ˚ A, 1193 ˚ A, and 1260 ˚ A, and one Si III line at 1207 ˚ A. Other lines are present but can be neglected. While all of the above transitions will also have an auto-correlation term that we account for in the modeling, it is their cross-correlation with Ly α that will allow us to differentiate their signals and measure their biases. Other metal transitions at longer wavelengths and lower redshifts only contribute significantly with their auto-correlation. Those cannot be easily separated from the Ly α auto-correlation signal. In a companion paper [32], we show that they are dominated by the C IV absorption, and that one can measure their contribution from the auto-correlation in the side bands (at wavelengths larger than the Ly α line in the quasar rest-frame). Following a first estimate from their analysis 16 , we use a prior on the effective CIV bias of b eff CIV = -0 . 0243 ± 0 . 0015 with β CIV = 0 . 5 which combines the signal from C IV and other transitions (notably Mg II and Si IV).  \nWe use the same set of metal absorbers and priors when modelling the Ly α (A) × Ly α (B) correlations. Absorbers with wavelength close to the Ly β line like the O VI lines at 1032 ˚ A and 1038 ˚ A do not contribute in cross-correlation with Ly α to our measurements because  \ntheir wavelength is much smaller that the Ly α line. Their auto-correlation, peaking at zero separation, does not contribute either because we do not measure the auto-correlation of pixels in the B region, Ly α (B) × Ly α (B).  \nAs explained in [31], for each pair of absorbers ( m,n ) we compute a metal matrix that provides the mapping between the true co-moving separation ( r ∥ , r ⊥ ) of the two absorbers and their apparent separation when assuming both are caused by the Lymanα transition. In previous works this mapping was computed numerically for a small fraction of the total number of pairs in the sample. This estimation was not precise enough at small separation and expensive in computing time. We now compute uniquely the shift along r ∥ and ignore the few percent change in r ⊥ . This simplification allows us to measure more precisely the effect with one dimensional integrals using the sum of the weights as a function of wavelength.", '4.5 Correlated noise from the data processing': 'The data processing [37] introduces correlated noise among the spectra from fibers of the same spectrograph, which correspond to fibers from a unique petal in the focal plane 17 .  \nThis contamination is studied in a companion paper [32], where we show that the dominant contribution is the sky background model noise (as in BOSS/eBOSS, see [21] and dMdB20 ). We find that the following expression is sufficient to describe this contamination to the DESI Ly α auto-correlation function  \nξ cont ( r ∥ , r ⊥ ) = a noise δ K ( r ∥ ) f ( r ⊥ ) ≡ a noise ξ noise ( r ∥ , r ⊥ ) (4.2)  \nwith a noise an amplitude and f ( r ⊥ ) a decreasing function of r ⊥ proportional to the fraction of pairs at r ⊥ that belong to the same petal. This function is evaluated numerically in [32] assuming pairs at z = 2 . 4 for the fiducial cosmology to convert angles to co-moving separations. We have f ( r ⊥ > 110 h -1 Mpc) = 0 when the separation exceeds the size of a petal. We chose an arbitrary normalization f (0) = 1 such that a noise is close to the value of this contamination in the first ( r ∥ , r ⊥ ) bin.', '4.6 HCD contamination': "The presence of High Column Density systems (HCDs) in the quasar spectra, including Lyman Limit Systems (LLS, log N HI > 17 . 2) and Damped Lyman α systems (DLA, log N HI > 20 . 3), complicates the modeling of 3D correlations in the Ly α forest [84-86]. Like the Ly α forest itself, HCDs are tracers of the underlying matter density and on very large scales (larger than the width of their absorption profiles) their contamination is limited to a change in the linear bias parameters of the Ly α forest (see section 4.2 in [85]). However, the damping wings of DLAs can extend to fairly large line-of-sight separations, and adds an extra scale dependence to the correlation function on scales of tens of megaparsecs.  \nTo diminish the effect of the HCDs on the correlation function, we mask the highly absorbed wavelength range of identified Damped Lyman alpha systems (DLAs). As described in section 2, we have good efficiency for identifying DLAs in the spectra with high signalto-noise. However, the smoothing effect of unidentified HCDs remains, and this must be modeled.  \nWe use a scale-dependent Ly α bias of the form b ' α = b α + b HCD F HCD ( k ∥ ) and a similar form for the RSD parameter β α . As discussed in appendix A (see also [87]), the form of F HCD ( k ∥ ) and the magnitude of b HCD can be related to the column density distribution  \nof unmasked HCDs and the bias of their host halos. The k ∥ dependence is given by the Fourier transform of the HCD Voigt profiles (an absorption profile with a Gaussian core and Lorentzian tails). If the column-density distribution of the unmasked HCDs was known, it would then be possible to calculate b HCD F HCD ( k ∥ ) . Unfortunately, at present we do not know precisely the efficiency of the DLA detection in noisy spectra, and we only know approximately the bias of halos hosting DLAs [76, 88, 89].  \nBecause of this, and motivated by the fact that the Fourier transform of a Lorentzian is an exponential function, we model the contamination with F ( k ∥ ) = exp( -L HCD k ∥ ) and treat b HCD and L HCD as free parameters of the model. The parameter β HCD is poorly determined in the fits and we choose to use a prior β HCD = 0 . 5 ± 0 . 1, a value motivated by the measured bias of b DLA ∼ 2 of the DLA hosts. In the dMdB20 analysis, L HCD = 10 h -1 Mpc was fixed, but in the present analysis we vary this parameter and find L HCD = (6 . 51 ± 0 . 9) h -1 Mpc. As discussed in section 6.4, variations in the modelling of the HCD contamination have a negligible impact on the BAO results.", '4.7 Quasar redshift errors': "Similar to the impact of random peculiar velocities (or 'Fingers of God'), random errors in the estimation of quasar redshifts dilute the clustering of quasars along the line of sight. This has an impact in the cross-correlation of quasars and the Ly α forest, as first discussed in [90]. Following dMdB20 , in our main analysis we model this smoothing with a Lorentzian with free parameter σ z , but in section 6.4 we show that using a Gaussian instead has a negligible impact on the BAO results.  \nA small systematic error in the quasar redshifts would be very difficult to detect in the auto-correlation of quasars. However, such an offset would shift the Ly α × QSO crosscorrelation such that it would no longer peak at r ∥ = 0 18 . We parameterize this shift with a free parameter ∆ r ∥ .  \nThe impact of redshift errors in the Ly α × QSO cross-correlation can be seen as a nuisance in Ly α BAO studies, but as discussed in [49] it also provides a great diagnosis tool to better calibrate the redshifts of quasars.  \nWe find that the quasar redshifts are unbiased, with ∆ r ∥ = (0 . 066 ± 0 . 058) h -1 Mpc. We also find that the combination of random peculiar velocities and redshift errors are described by a Lorentzian with σ z = (3 . 67 ± 0 . 14) h -1 Mpc, significantly smaller than the value reported in the eBOSS analysis of dMdB20 , σ z = (6 . 86 ± 0 . 27) h -1 Mpc. The reduced redshift errors can be explained by the updated quasar templates used in Redrock (see table 6 of [47]).", '4.8 Quasar radiation (proximity effect)': 'Quasars are some of the brightest objects in the Universe. Therefore, we expect them to significantly ionize their surroundings, an effect sometimes referred to as the proximity effect . The Ly α forest of neighbouring quasars is therefore affected by two competing effects: the gas density is higher than average (quasars live in high-density regions), but the neutral fraction is lower than average (due to the quasar radiation).  \nFollowing [90], we use a simple model to account for the proximity effect in the Ly α × QSO cross-correlation. In particular, we use the implementation of dMdB20 that assumes isotropic radiation from the quasar, a mean-free path of UV photons of λ UV = 300 h -1 Mpc, and has a single free parameter ξ TP 0 that sets the amplitude of the contamination.', '4.9 Small-scales corrections': 'The BAO parameters only shift the peak component of the model, and as discussed in [15] this is by construction zero on scales smaller than 80 h -1 Mpc. Therefore, one could decide to limit the BAO analysis to these very large separations. However, some of the nuisance parameters described in this section can only be constrained when extending the analysis to smaller separations. This is the case for the parameters describing quasar redshift errors, or the parameter describing the Si III line at 1207 ˚ A that causes a sharp feature at ( r ⊥ ∼ 0, r ∥ ∼ 20 h -1 Mpc). Moreover, the distortion matrix discussed in section 4.3 spreads the impact of some of these small-scale effects to larger separations (in particular along the line of sight). For this reason, in our main analysis we include the measurement of correlations down to r > 10 h -1 Mpc, and in section 6.4 we show that the BAO results do not depend on the minimum separation 19 .  \nIn order to fit the Ly α auto-correlation to 10 h -1 Mpc, we follow previous Ly α BAO analyses and use a multiplicative correction to the model of the Ly α power spectrum. In particular we use the correction from [91], calibrated with hydrodynamical simulations, that models both the effect of non-linearities in the densities and velocities, but also thermal effects in the IGM (thermal broadening, pressure) 20 . In section 6.4 we show that this correction has a negligible impact on our BAO parameters.  \nAn equivalent model was proposed in [92] for the cross-correlation of the Ly α forest with halos of intermediate masses. However, the simulations used were too small to contain enough massive halos to accurately study the clustering of quasars. We decided to not include this correction in our model. However, as discussed in section 4.7 we do take into account the impact of non-linear peculiar velocities on the Ly α × QSO cross-correlation.  \nFinally, following [31] and previous Ly α BAO analyses, we take into account the finite size of our correlation function bins.', '5 Measurement of Baryon Acoustic Oscillations': 'After presenting the measured correlations in section 3 and discussing the model we used to describe them in section 4, in this section we summarise the statistical method used to fit the measurement and present the main results of this paper: the measurement of the BAO scale along and perpendicular to the line of sight.  \nWe use the Vega package both for the modelling of the correlations and for the parameter inference. We use a Gaussian likelihood, and the main results presented in this section were obtained using the Nested Sampler Polychord [93, 94]. The best-fit values are the mean of the posteriors, and the reported uncertainties are the 68% credible intervals. However, this analysis is computationally intensive, and in most of the tests in section 6 we use instead a simpler method: we use the iminuit software [95, 96] to find the maximum of the likelihood, and use the derivatives of the likelihood around the best-fit point to estimate a Gaussian posterior. As discussed in appendix D, both BAO estimates are very similar.  \nWe start with a data vector composed of 4 different 2-point functions, two of them with 2 500 data points (Ly α (A) × Ly α (A) and Ly α (A) × Ly α (B)) and two of them with 5 000 data points (Ly α (A) × QSO and Ly α (B) × QSO), for a total of 15 000 data points. While these are mostly independent, as discussed in section 3 we include their small cross-covariance, and  \nFigure 7 . Measurements of the BAO parameters along the line of sight ( α ∥ ) and across the line of sight ( α ⊥ ), with contours corresponding to the 68% and 95% confidence regions. The auto-correlation results (filled blue contours) are the combined measurement of the Ly α forest auto-correlations in the regions A and B. The cross-correlation results (dashed black) are the correlations of the forest in these two regions with quasars. The combined results (solid red) simultaneously fit all four correlations taking into account their cross-covariance, and are the main result of this publication.  \n<!-- image -->  \ntherefore we use a 15 000 × 15 000 covariance matrix. However, following dMdB20 we limit the range of separations used in the fits to 10 < r < 180 h -1 Mpc, reducing the number of data points used in the combined fit to 9540 (see table 3) 21 .  \nThe constraints on the BAO parameters are listed in table 3 and shown in figure 7, where we show constraints from the Ly α auto-correlation (in filled blue contours, including both Ly α (A) × Ly α (A) and Ly α (A) × Ly α (B)) as well as constraints from the cross-correlation with quasars (dashed black, including both Ly α (A) × QSO and Ly α (B) × QSO). Both mea-  \nTable 3 . Best fit BAO parameters (mean of the posterior), uncertainties (68% credible intervals) and correlation coefficient ρ from the three main analyses: auto-correlations (Ly α (A) × Ly α (A) and Ly α (A) × Ly α (B)), cross-correlations (Ly α (A) × QSO and Ly α (B) × QSO) and their combination. All parameters are given at z eff = 2 . 33. The p-value is only accurate for the combined analysis, because in the other analyses we fix the value of one of the nuisance parameters to the best-fit value in the combined analysis (see discussion in appendix B).  \nre consistent, and their combined constraints (shown in solid red lines) are α ⊥ = 1 . 013 ± 0 . 024 and α ∥ = 0 . 989 ± 0 . 020, with a correlation coefficient of ρ = -0 . 48.  \nAs discussed in section 4, in addition to the two BAO parameters our model has 15 nuisance parameters that we marginalise over. The number of degrees of freedom in the combined fit is 9523 (9540-17), and the χ 2 of the best-fit model is 9624.36, with a probability of having a value larger than this of 23%. The best-fit values of the nuisance parameters are discussed in appendix B. Some of these nuisance parameters only affect the auto-correlation or the cross-correlation, and are therefore ignored when fitting these correlations individually (see table 3 and table 5). Moreover, when analysing these correlations separately we are not able to break some of the degeneracies between nuisance parameters, and we use extra priors as described in table 5.  \nIn the latest Ly α BAO analyses from eBOSS, the auto-correlation provided ∼ 20% better constraints on α ∥ than the cross-correlation, while providing ∼ 10% weaker constraints on α ⊥ (see figure 12 in dMdB20 ). Redshift space distortions in the quasar dataset ( β Q ∼ 0 . 3) are milder than in the Ly α forest ( β α ∼ 1 . 7), reducing the constraining power along the lineof-sight direction. On the other hand, the Ly α BAO measurement from DESI DR1 seems to be dominated by the cross-correlation for both α ⊥ and α ∥ (see figure 7). Most Ly α quasars in DESI DR1 have only one observation (see right panel of figure 1). In future data releases their signal-to-noise will increase as we collect more observations, and the constraining power of the auto-correlation (doubly affected by noise in the spectra) will increase more than that of the cross-correlation.', '6 Analysis validation': 'Baryon Acoustic Oscillations (BAO) imprint a characteristic three-dimensional feature in the measured correlations. On the other hand, most spurious correlations from instrumental systematics and astrophysical contaminants are smooth and featureless 22 . This makes BAO  \nmeasurements particularly robust.  \nHowever, some of the analysis choices presented in section 2 cause small changes in the dataset that introduce statistical fluctuations in the BAO measurement. Examples of these are the observed wavelength range, the rest-frame limits of the Ly α regions A and B, or the masking of pixels (due to sky lines, DLAs or BAL features). Moreover, differences in the quasar redshift estimators cause pixels to fall in and out of the A and B regions, adding an extra source of statistical fluctuations that is also seen in mocks (see appendix B of dMdB20 ).', '6.1 Blinding strategy': 'In order to avoid unconscious or confirmation biases, the development and testing of the analysis framework defined in sections 2 to 5 was done using synthetic datasets and blinded measurements. We considered several blinding strategies for the DESI Ly α BAO analysis, including the possibility of blinding the data at the catalog level as done in the galaxy BAO analysis of DESI [6]. However, the presence of known sky lines in the spectra, as well as the presence of absorption lines with small restframe wavelength separations from Ly α , made it challenging to apply a robust blinding to the data at the catalog level.  \nFor this reason, we opted instead for blinding the measured correlation functions following a simple method described in appendix C. In short, we applied an additive correction to the measured correlation functions to mimic a blind shift in the BAO parameters.  \nWe defined a list of tests that we needed to pass in order to validate the analysis before unblinding the measurement. First, we validated the analysis using synthetic datasets (or mocks), as explained in detail in a companion paper [33], and as summarised in section 6.2. Second, we studied the consistency of the results under various data splits, as discussed in section 6.3. Finally, we tested the robustness of the results under variations in the analysis setup, as described in section 6.4. We report here test results applied to the unblinded data, but the same tests were first performed on the blinded data set.', '6.2 Validation using synthetic data': 'A detailed description of the procedure to generate synthetic DESI spectra (or mocks) for Ly α studies can be found in [34]. The analysis validation of the Ly α BAO measurement using mocks is presented in a companion paper [33]. Here we give a brief summary of the mocks, and we show some of the main tests validating the analysis.', '6.2.1 DESI Ly α mocks': 'We generated DESI mocks from two different sets of fast simulations: 100 realisations of LyaCoLoRe mocks [97, 98] and 50 realisations of Saclay mocks [99]. Both sets of mocks use a log-normal description of the density field, and use simplified recipes to distribute quasars and simulate the optical depth of Ly α absorption in redshift space. These recipes were calibrated in order to approximately reproduce the mean flux, the 1D power spectrum, and the large-scale biases of the Ly α forest as measured by the eBOSS Collaboration. The LyaCoLoRe simulations were also used in the final Ly α BAO analysis of eBOSS presented in dMdB20 .  \nAs described in [34], these simulations are post-processed with the script quickquasars of the desisim package 23 , where the DESI specificities are introduced, namely the footprint,  \nsignal to noise, spectrograph resolution, quasar redshift errors, etc. At the same time, astrophysical contaminants are introduced such as Damped Lyman alpha systems (DLAs), Broad Absorption Line features (BALs), and absorption from metal lines.  \nWe have made two small changes in the procedure with respect to the description in [34]. First, we have improved the way in which we imprint the footprint inhomogeneities caused by the survey strategy of DESI. Second, we have slightly modified the recipes to add metal absorption to the mocks in order to better match the amount of metal contamination seen in the data. Both of these changes are discussed in more detail in [33].  \nWe analyse these mocks using the same analysis applied to the data, with the following minor exceptions in the modelling: (i) we ignore the contamination from CIV, the transverse proximity effect and correlated sky residuals, since these are not included in the mocks; (ii) we include an extra smoothing to the model to account for pixelisation effects coming from the finite cell size of the log-normal simulations; (iii) we do not smooth the peak component of our model since we use lognormal simulations that do not capture the non-linear broadening of the BAO peak. These differences are also discussed in more detail in [33].  \nThe contaminants that we cannot study with synthetic data are studied with variations in their modelling in section 6.4 or discussed in a companion paper. For instance, the CIV contamination and the correlated signal from the sky subtraction are studied in [32] as discussed in section 4.5. Given that the fitted amplitudes of these contaminants are uncorrelated with α ∥ and α ⊥ (see figure 14 of appendix B), we are confident that they do not bias our estimate of the BAO parameters. The lack of non-linear broadening of the BAO peak in the mocks results in an underestimation of the BAO uncertainties measured from mocks (see figure 9). As discussed in section 4, when fitting the DESI DR1 measurement we model the broadening using calculations based on Lagrangian Perturbation Theory [15]. This model reproduces the shape of the BAO peak in the correlations measured from Ly α mocks constructed from N-body simulations [100].  \nFinally, the algorithm used to identify DLAs in the data requires a significant amount of computing time, and therefore we decided to not run it on the many realisations of mocks. Instead, we assume that we find all HCDs with log N HI > 20 . 3, and none below this column density, and we mask them in the analysis. Similarly, we assume that we can find all BAL features in the data, and mask the corresponding region of spectra accordingly. The impact of BAL completeness is studied in detail in [62].', '6.2.2 Validation of the covariance matrix': "In order to validate the covariance matrices of the correlation functions for the purpose of measuring the BAO scale, we measured the quantity  \nδ α ≡ ( ∂ α M T C -1 ∂ α M ) -1 ∂ α M T C -1 R (6.1)  \nfor each mock, where ∂ α M is the derivative of the best fit model from a stack of mocks with respect to a BAO scale parameter α , C the covariance matrix of the correlation functions and R the data minus the best fit model from the stack. The covariance C is determined for each mock independently with the sub-sampling and smoothing method described in section 3.3, so the noise of the covariance itself is directly comparable to that of the true data.  \nThis is not a true fit (which involves a non-linear model plus many nuisance parameters, see section 4) but a straightforward compression of the data and its covariance that is best suited to describe the fluctuations that matter for measuring α . We also measure the  \nFigure 8 . Left: BAO constraints from the 'stack' of 100 LyaColore (blue) and 50 Saclay (orange) DESI DR1 synthetic datasets, compared to the constraints from data (scaled down by a factor of 3, dashed black ellipse). Right: scatter plot of the best fit α ⊥ and α ∥ from each of the LyaColore (blue) and Saclay (orange) mocks. Note the difference of scale between the two plots.  \n<!-- image -->  \n1.06  \n1.04  \n1.02  \n1.00  \n0.98  \n0.96  \n0.94  \n0.92  \n0.94  \n0.96  \n0.98  \n1.00  \n1.02  \n1.04  \n1.06  \n1.08  \nuncertainty on this parameter from the covariance,  \nσ α ≡ ( ∂ α M T C -1 ∂ α M ) -1 / 2 . (6.2)  \nWe then measure the rms of the distribution of ( δ α /σ α ). If the covariance matrix is correct we expect this distribution to be Gaussian with a rms of 1. We looked at α ∥ and α ⊥ and linear combinations of both (major and minor axes of 2D uncertainty contour). When using the full covariance matrix (including the cross-covariance between the correlation functions), we find a scatter of 0 . 96 ± 0 . 05 (0 . 99 ± 0 . 06) for α ∥ ( α ⊥ ) for the LyaCoLoRe mocks, and 1 . 07 ± 0 . 07 (1 . 03 ± 0 . 07) for the Saclay mocks. This shows that the statistical uncertainties derived from the covariance matrix are a good estimate of the dispersion among random mock realizations. This gives us confidence in the covariance matrix at the scales of interest for the measurement of the BAO. As discussed in more detail in [33], the scatter was 10% larger for both sets of mocks when ignoring the cross-covariance between the correlation functions.", '6.2.3 Validation of BAO estimates': "We measured the correlations in 100 LyaCoLoRe mocks and 50 Saclay mocks mimicking the DESI DR1 dataset, and combine their measurements of the correlations into 'stacks of correlations', with a statistical power much larger than that of the final DESI dataset. The BAO constraints from these stacks are shown in the left panel of figure 8. We used the same cosmology to generate the mocks and to analyse them, so we should expect to recover values of ( α ⊥ , α ∥ ) consistent with unity.  \nIn order to proceed with the unblinding, we had set a requirement that the measurement on the stack of many mocks could not present a bias larger than a third of the statistical uncertainty obtained when fitting the blinded data 24 . This corresponded to a tolerance of ∼ 0 . 005 in α ∥ and ∼ 0 . 007 in α ⊥ .  \n|  \n|  \nDistribution of BAO results in mocks  \nLy Colore  \nSaclay  \nFigure 9 . Distributions of BAO uncertainties along ( σ α || , left) and across ( σ α ⊥ , right) the line of sight. The blue histogram shows the distribution from the analysis of 150 DESI DR1 mocks, while the vertical dashed lines are the uncertainties measured in the data. The solid black line shows the distribution of BAO uncertainties from Monte Carlo realisations of the data covariance matrix, when using the best-fit model. The solid red line, on the other hand, shows an equivalent distribution from Monte Carlo realisations generated around a linear model that does not include the non-linear broadening of the BAO peak. These Monte Carlo realisations are discussed in detail in [33].  \n<!-- image -->  \nThere is a small bias in the measurement of BAO from the stack of 100 LyaCoLoRe mocks (blue contours), but it is smaller than the requirement accuracy (black dashed contours). The difference between the results from the two sets of mocks also suggests that any offsets seen due to analysis problems are at the same level as systematics in the creation of the mocks. Moreover, combining these with the results from the stack of 50 Saclay mocks (orange contours) would further reduce the bias.", '6.2.4 Validation of BAO uncertainties': 'We discuss here the distribution of BAO results when fitting each of the DESI DR1 mocks individually, and use them to validate several aspects of our analysis. In the right panel of figure 8 we show the best-fit values for the BAO parameters ( α ∥ and α ⊥ ) for the 50 Saclay and the 100 LyaCoLoRe mocks, with errorbars representing the 1σ uncertainties. In the following we combine the results from these 150 mocks.  \nFigure 9 shows the distribution of uncertainties on α ∥ and α ⊥ in the mocks (blue histogram) and compares it to the uncertainty measured in DESI DR1 (vertical dashed line). One can see that the BAO uncertainties vary significantly from mock to mock, as expected, but that the uncertainties from DESI DR1 are larger than those from mocks. The BAO uncertainties from analyses of the DESI DR1 mocks are expected to be a bit smaller than in the data, since the mocks used in this analysis do not have non-linear broadening of the BAO peak. In order to study this, we generated 1 000 Monte Carlo (MC) realisations of the correlation functions using the covariance matrix from the data, adding the random fluctuations to the best-fit model from our main analysis 25 . The distribution of BAO uncertainties  \nfrom these MC realisations is shown in black in figure 9. It clearly shows that the DESI DR1 Ly α BAO result is not an outlier, and that the constraining power of the mocks is larger than that of the data. In the same figure we show (in red) the distribution of uncertainties from another set of 1 000 MC realisations, where the fluctuations have now been added to a model that ignores the non-linear broadening of the BAO peak. This distribution is in very good agreement with the distribution of uncertainties from the fits on individual mocks, and confirms the hypothesis that the non-linear broadening of the BAO peak degrades the Ly α BAO result significantly.  \nAs discussed in [33], the χ 2 value from the data is consistent with the distribution of χ 2 in the mocks. In the same publication, we also look at the distribution of BAO residuals (∆ α ∥ /σ α ∥ and ∆ α ⊥ /σ α ⊥ ). Their rms is found to be of 1 . 01 ± 0 . 07 for α ∥ and 1 . 11 ± 0 . 06 for α ⊥ , with uncertainties obtained through bootstrap. Those values which are close to one validate our error propagation. We note that this is a more comprehensive test than the one presented in section 6.2.2 because it is based on the results from the full non-linear fit that includes numerous nuisance parameters.  \nFrom the tests discussed above we conclude that our BAO estimates on mocks are unbiased at the level of precision required by the current dataset, and that the scatter of best-fit values is consistent with the reported uncertainties.', '6.3 Data splits': "The second set of tests that we use to validate the analysis are data splits, where we measure BAO using different subsets of the data. A first data split was already shown in figure 7, where we presented the consistency of BAO measurements from the auto-correlations (including both the Ly α (A) × Ly α (A) and the Ly α (A) × Ly α (B) correlations) and from the cross-correlations (including Ly α (A) × QSO and Ly α (B) × QSO). In the bottom right panel of figure 10 we group them instead in correlations that only use pixels in the A region (Ly α (A) × Ly α (A) and Ly α (A) × QSO, in green) and correlations that use pixels in the B region (Ly α (A) × Ly α (B) and Ly α (B) × QSO, in blue). The B region does not contain as much information as the A region for several reasons: only quasars at higher redshift have this region in the DESI spectrograph, so the number of Ly α lines-of-sight in the B region is smaller; the B region, as defined in section 2.4, is significantly shorter than the A region, so there are fewer pixels per line-of-sight; finally, the B region is affected by other Lyman lines that add extra variance to the fluctuations.  \nWe will now discuss the consistency of the BAO constraints when splitting the quasar catalog in three ways: by imaging survey used in the target selection; by C IV emission line equivalent width (EW); and by signal-to-noise in the spectrum. In these cases we run alternative end-to-end analyses starting from new sub-catalogs, i.e., fitting new quasar continua, measuring correlations and fitting them separately for each subset of the data 26 .  \nWe start in the bottom left panel of figure 10 by splitting the catalog based on the imaging survey that was used for target selection. Most of the DESI footprint was observed with the DECam camera on the Blanco telescope in Chile, including the entire South Galactic Cap and the southern fraction of the North Galactic Cap [101-103]. We refer to this subset of the data as 'South', while we designate the area that was imaged in the BASS and MzLS surveys at δ> 32 . 375 · as 'North'. The South sample is significantly larger, as it contains 82% of the quasars (579 166 quasars with z > 1 . 77 in the South and 130,399 in the North).  \nFigure 10 . BAO constraints from the main analysis (grey) and from data splits. Top left: low (green) vs high (blue) SNR in the quasar spectrum. Top right: low (green) vs high (blue) CIV equivalent width (EW) in the quasar spectrum. Bottom left: South (green) vs North (blue) imaging used in the quasar target selection. Bottom right: correlations from region A (green) and region B (blue); the A region shows the combined measurement from the auto-correlation of the forest measured in the A region (Ly α × Ly α ) and the cross-correlation of this region with quasars (Ly α (A) × QSO). The contours labeled region B show the combined measurement of the forest auto-correlation measured in the B region (Ly α (A) × Ly α (B)) and the cross-correlation of this region with quasars (Ly α (B) × QSO).  \n<!-- image -->  \nIn the top right panel of figure 10 we look at a second catalog split based on the C IV EW. We do this split because we expect the shape of the quasar spectral energy distribution to depend on EW, due to the well known anti-correlation between the EW of quasar emission lines and the continuum luminosity known as the Baldwin Effect [104]. We measure the C IV EW of the quasars with fastspecfit 27 , finding a median of 37.3 ˚ A for all quasars, and a median of 41.6 ˚ A for 3 σ measurements of the C IV emission line. We split the sample at 39 ˚ A, which produce a low (high) EW sample of 371 751 (337,814) quasars at z > 1 . 77. As expected from the Baldwin Effect, the low EW sample is somewhat higher luminosity and has a somewhat higher effective redshift of 2.36 compared to 2.29 for the high EW sample.  \nFinally, in the top left panel of figure 10 we present the third quasar catalog split, based on mean signal-to-noise ratio (SNR) in the spectra. Instead of splitting the quasar catalog into two subsets of equal size, we chose a SNR threshold of 2.25 such that both subsets have the same weight in the measurement of the auto-correlation function. There are different ways of estimating the SNR of quasar spectra, and we decided to use the mean value of SNR per pixel averaged over the Ly α region, as reported by the picca code at the end of the continuum fitting process of the main analysis. This results in a lower SNR catalogue with 321 767 quasars and a higher SNR catalogue with 106 636 quasars. The sum of these  \ncatalogues does not match the total size of the catalogue used in the main analysis (709,565 quasars at z > 1 . 77), since we do not detect the forest continuum at z < 2 and therefore do not have a SNR estimate.  \nWith the exception of the North vs South data split, the subsets discussed here share the same footprint and redshift range. However, cosmic variance is a very small component of the covariance matrix and to a first approximation the data splits can be considered as independent. Taking this into account, the BAO constraints on the various data splits are consistent with statistical fluctuations.", '6.4 Alternative analyses': "In this section we show a final set of validation tests: robustness of the BAO measurement under variations in the analysis configuration. We set a requirement for unblinding that none of the variations would cause a shift on the BAO parameters larger than a third of the statistical uncertainty from the main (blind) analysis. This corresponded to a tolerance of ∼ 0 . 005 in α ∥ and ∼ 0 . 007 in α ⊥ . However, some of the analysis variations result in a small change in the size of the dataset. In these cases we relax the requirement and take into account the probability of the measured shift being caused by statistical fluctuations. These are discussed in appendix E.  \nIn order to reduce the amount of computing time needed in these alternative analyses, we do not run the nested sampler algorithm and only report the maximum likelihood values and Gaussian errors computed by iminuit . As discussed in appendix D, the BAO results do not vary significantly with the sampling method used. In some of the alternative analyses iminuit has difficulty breaking internal degeneracies between nuisance parameters, particularly between L HCD and some of the bias parameters. In order to avoid this, we fix L HCD in all the alternative analyses to the best-fit value of the main analysis (6 . 51 h -1 Mpc), and we show in the variation 'vary L HCD ' that this has a negligible impact on the BAO results.", '6.4.1 Variations in the estimation of the fluctuations': "The method to estimate the Ly α fluctuations starting from the observed quasar spectra is described in section 2.5. In the first set of variations shown (in purple) in figure 11 we quantify the impact on the BAO results when we vary different aspects of the method. In particular we look at the following variations:  \n- · no calibration: We do not re-calibrate the spectra with the C III region mentioned in section 2.4 and described in [30].\n- · η pip = 1: We do not apply the re-calibration of the instrumental noise η ( λ ) mentioned in section 2.4 and described in [30].\n- · ϵ free: We include an extra term ( ϵ in dMdB20 ) in the computation of the Ly α weights to try to capture quasar diversity.\n- · η LSS = 3 . 5: We reduce the contribution from the intrinsic Ly α forest variance to the weights by a factor of two.\n- · ∆ λ = 2 . 4 ˚ A: We coadd three pixels into one before the continuum fitting step (as done in dMdB20 ). In this variation we also use σ 2 mod = 3 . 1 as suggested by [30] when coadding DESI data by three pixels.  \nFigure 11 . Shifts in the BAO parameters from alternative analyses, including variations in the method to estimate the fluctuations (purple); variations in the dataset (red); variations in the method to compute correlations and covariances (green); variations in the range of separations used (orange); and variations in the modelling (blue). The red shaded contours show the one σ uncertainty from the main analysis, while the smaller gray area shows the threshold set to these tests ( σ/ 3). Note that the two parameters are anti-correlated ( ρ = -0 . 48). Variations of the dataset (in red) are subject to statistical fluctuations as described in appendix E.  \n<!-- image -->  \nWe continue with a second set of variations (in red) in figure 11 where we look at variations that cause changes in the dataset by removing (or adding) pixels or entire quasars. As discussed in appendix E, these can cause statistical fluctuations in the BAO results. In particular we look at the following variations:  \n- · λ obs < 5500 ˚ A: We use only Ly α pixels below this observed wavelength ( λ obs < 5577 ˚ A in the main analysis).\n- · λ obs > 3650 ˚ A: We use only Ly α pixels above this observed wavelength ( λ obs > 3600 ˚ A in the main analysis).\n- · λ RF < 1200 ˚ A: We use only Ly α pixels below this rest-frame wavelength ( λ RF < 1205 ˚ A in the main analysis).\n- · z Q < 3 . 78: We use only quasars with z Q < 3 . 78 (highest redshift included in the mocks discussed in section 6.2).\n- · > 50 pixels in forest: We include lines-of-sight with more than 50 valid Ly α pixels (150 in the main analysis).\n- · original redshift estimates: We use the quasar redshifts in the original run of Redrock (slightly biased as discussed in [47]).\n- · mask-Lya redshift estimates: We use a different estimator for the quasar redshifts, rerunning Redrock with only wavelengths longer than the Ly α emission line (as done in dMdB20 ).\n- · only quasar targets: We use only quasars that were considered quasar targets.\n- · DLAs SNR > 1: We mask DLAs in spectra with SNR > 1 (SNR > 3 in the main analysis).\n- · weak BALs: We remove spectra with strong BAL features ( A I > 840, 50% percentile of strongest BALs in eBOSS [61]).\n- · no sharp lines mask: We do not mask the 4 sharp lines discussed in section 2.4.  \nSome of the variations (like 'mask-Lya redshift estimates') are slightly outside the threshold of 1 / 3 σ . As we discuss in appendix E, these shifts can be explained by statistical fluctuations caused by the addition or subtraction of Ly α pixels from the dataset.", '6.4.2 Variations in the measurement of correlations': 'We now move to the third set of variations (in green) shown in figure 11, where we look at the impact of varying the setup in the measurement of correlations, their covariances and the distortion matrices described in section 3. In particular we look at the following variations:  \n- · dmat r ∥ < 200 Mpc/h: We model the distortion matrix up to r ∥ = 200 h -1 Mpc as done in dMdB20 ( r ∥ = 300 h -1 Mpc in the main analysis).\n- · dmat 2%: We use 2% of pixels to compute the distortion matrix (1% in the main analysis).  \nTable 4 . Extra Gaussian priors added to some of the variations discussed in section 6.4. They correspond to the best-fit values and uncertainties from the main analysis as reported in table 5.  \n- · dmat model 4 Mpc/h: We model the distortion matrix using the same 4 h -1 Mpc binning as the correlation function (2 h -1 Mpc in the main analysis).\n- · ∆ λ = 3 . 2 ˚ A: We rebin the Ly α fluctuations in groups of 4 pixels (3 in the main analysis).\n- · ∆ λ = 1 . 6 ˚ A: We rebin the Ly α fluctuations in groups of 2 pixels (3 in the main analysis).\n- · nside = 32: We measure the correlations in HEALPix pixels defined by nside=32 (16 in the main analysis).\n- · ∆ r = 5 Mpc/h: We use 5 h -1 Mpc bins in the correlation functions (4 h -1 Mpc in the main analysis).\n- · no cross-covariance: We ignore the cross-covariance of the different correlation functions (as done in eBOSS, dMdB20 ).  \nWe do not see any problematic variation related to the measurement of correlations and their covariances.', '6.4.3 Variations in the parameter estimation': 'We move now to the impact of variations in the parameter estimation. We start with the four set of variations (orange) shown in figure 11 by looking at the impact of the range of separations used:  \n- · r < 200 Mpc/h: We fit separations with r < 200 h -1 Mpc (180 h -1 Mpc in the main analysis).\n- · r < 160 Mpc/h: We fit separations with r < 160 h -1 Mpc (180 h -1 Mpc in the main analysis).\n- · r > 20 Mpc/h: We fit separations with r > 20 h -1 Mpc (10 h -1 Mpc in the main analysis).\n- · r > 40 Mpc/h with priors: We fit separations with r > 40 h -1 Mpc (10 h -1 Mpc in the main analysis). Without the smaller scales we are not able to constraint several nuisance parameters, and therefore we add the informative priors described in table 4.  \nFinally, in the last set of variations (in blue) shown in figure 11 we look at the impact of different modelling choices. In particular we look at the following variations 28 :  \n- · eBOSS metals: We model the contamination by Silicon lines (metals) following the method used in eBOSS [27] instead of the new method described in section 4.4.\n- · vary L HCD : We vary the parameter L HCD in the model of the contamination by HCDs as done in the main analysis. This parameters was fixed to L HCD = 6 . 51 h -1 Mpc in the other variations to minimise the degeneracies between this and other nuisance parameters. See the related discussion in appendix A.\n- · L HCD = 10 Mpc/h: We use a fixed value of L HCD = 10 h -1 Mpc to model HCD contamination (free parameter in the main analysis), as was done in dMdB20 .\n- · L HCD = 3 Mpc/h: We use a fixed value of L HCD = 3 h -1 Mpc to model HCD contamination (free parameter in the main analysis).\n- · Gaussian redshift errors: We use a Gaussian distribution to model quasar redshift errors and quasar peculiar velocities (a Lorentzian distribution is used in the main analysis).\n- · weak CIV bias prior: We use a significantly weaker flat prior on the CIV bias parameter of -0 . 03 < b CIV < 0 (instead of -0 . 0243 ± 0 . 0015 in the main analysis).\n- · no small-scales correction: We ignore the small-scales multiplicative correction from [91] in the modelling of the Ly α auto-correlation.\n- · UV fluctuations: Following the prescription of [21], we model the impact of fluctuations in the UV background [105, 106] on the Ly α forest auto-correlation. We do not detect these fluctuations in our analysis.\n- · no sky residuals: We ignore the contamination from correlated sky residuals in the Ly α auto-correlation, discussed in section 4.5.\n- · no proximity effect: We ignore the proximity effect , the impact of quasar radiation in the cross-correlation discussed in section 4.8.  \nNone of the variations in the parameter estimation cause a significant shift in the inference of the BAO parameters. Some of the nuisance parameters, on the other hand, are more sensitive to changes in the analysis setup.', '6.4.4 Broadband polynomial corrections': 'As a final test to demonstrate the robustness of the BAO measurement, we run an alternative analysis where we introduce a flexible but smooth additive component to each of the four modelled correlations. In particular, we follow the procedure of [21, 27] and use (for each correlation) Legendre polynomials L j ( µ ) of order j = 0 , 2 , 4 and 6 to describe the angular dependence of the additive terms, divided by powers of r i with i = 0 , 1 , 2 (corresponding to a parabola in r 2 ξ ( r )). The total number of broad-band parameters is therefore 48 (12 for each of the four correlations).  \nWethen fit the baseline model with those additional parameters in the limited separation range 40 Mpc/h < r < 180 Mpc/h while adding extra priors on nuisance parameters as described in table 4 to break degeneracies between the polynomial coefficients and other parameters. The best fit model with and without those broadband corrections is compared to the data in figures 4 and 5. The shifts in the best fit BAO parameters when adding broadband terms are ∆ α ∥ = +0 . 001 and ∆ α ⊥ = -0 . 001, and the change in the uncertainties  \nis negligible ( | ∆ σ α | < 0 . 001, both for α ∥ and α ⊥ ). We also note that the value of χ 2 is reduced by only 47.5 points when adding 48 new parameters.', '6.5 Conclusion on the validation tests': 'We have presented numerous validation tests of the analysis using mocks, data splits, variations in choices made for the analysis, and adding ad-hoc broadband terms to improve the fit of the correlations around the BAO peak. All tests where the data set is left unchanged result in variations much smaller than 1/3 of the final statistical error. All data splits and other tests where the data set was altered are consistent with statistical fluctuations.  \nAdding broadband terms moved the best fit by less than a tenth of the statistical error. The largest offsets were found with mocks with average shifts of ∆ α ∥ = -0 . 003 ± 0 . 0014 and ∆ α ⊥ = +0 . 004 ± 0 . 0018 for the LyaColore mocks (see figure 8). Because we find similar discrepancies among the two sets of mocks (LyaColore and Saclay), we can not use those offsets to correct the measurements and we hence have to treat them as systematic uncertainties. Adopting a conservative systematic uncertainty of 0.5%, this results in a 3% (2%) increase of the total uncertainties on the longitudinal (transverse) BAO measurement when combined quadratically with the statistical errors. We consider that this increase in uncertainties is small enough to be ignored. As a result we only report statistical uncertainties in the following sections.', '7 Discussion': 'In section 5 we presented a measurement of the BAO parameters ( α ⊥ , α ∥ ) at z eff = 2 . 33 from DESI DR1. Thanks to the peak-smooth decomposition introduced in section 4.1, we interpret these parameters as:  \nα ∥ = D H ( z eff ) /r d [ D H ( z eff ) /r d ] fid , α ⊥ = D M ( z eff ) /r d [( D M ( z eff ) /r d ] fid , (7.1)  \nwhere D M ( z ) is the transverse comoving distance, D H ( z ) = c/H ( z ) and r d is the sound horizon at the drag epoch. Quantities with [] fid are computed with the fiducial cosmology (see table 2). In combination with the z < 2 BAO measurements from [6], our BAO measurement enables the state-of-the-art cosmological constraints presented in [9].', '7.1 Cosmological distances': "Substituting values from our fiducial cosmology, we can rewrite the ( α ∥ , α ⊥ ) constraints as the following constraints on ratios of distances:  \nD H ( z eff ) /r d = 8 . 52 ± 0 . 17 , D M ( z eff ) /r d = 39 . 71 ± 0 . 95 , (7.2)  \nwith a correlation coefficient of ρ = -0 . 48. For a given value of the sound horizon r d , these translate into a measurement of the expansion rate at z eff = 2 . 33:  \nH ( z eff ) = (239 . 2 ± 4 . 8) 147 . 09 Mpc r d km / s / Mpc (7.3)  \nand a measurement of the comoving transverse distance to z eff :  \nD M ( z eff ) = (5 . 84 ± 0 . 14) r d 147 . 09 Mpc Gpc . (7.4)  \nIt is also convenient to report the BAO information as an isotropic dilation parameter  \nD V ( z eff ) /r d = ( z eff D 2 M D H ) 1 / 3 /r d = 31 . 51 ± 0 . 44 , (7.5)  \nand an anisotropic (or Alcock-Paczy'nski [107]) parameter f AP :  \nf AP ( z eff ) = D M ( z eff ) D H ( z eff ) = 4 . 66 ± 0 . 18 . (7.6)  \nHowever, the ratio D V /r d is only the optimal definition of the isotropic BAO parameter in the absence of redshift space distortions. Every BAO measurement will have a different combination of ( D H , D M ) that will minimise the correlation with f AP and will therefore have a smaller relative uncertainty. From the posterior of the Ly α BAO measurement of DESI DR1, we find that this combination is approximately:  \nD M ( z eff ) 9 / 20 D H ( z eff ) 11 / 20 /r d = 17 . 03 ± 0 . 19 , (7.7)  \nwhich corresponds to a 1.1% measurement of the (isotropic) BAO scale at z eff = 2 . 33.", '7.2 Comparison with previous measurements from SDSS': 'Prior to this work, the BAO scale at z ∼ 2 . 3 had been measured by the BOSS and eBOSS collaborations with successively larger Ly α forest datasets [13-15, 19-24, 27]. The BAO measurements from SDSS DR11 [19, 20] showed a mild ∼ 2 . 3 σ tension with the best-fit ΛCDM model of Planck. As discussed in dMdB20 , this tension gradually disappeared with the addition of more data, and in the last measurement from eBOSS with SDSS DR16 the disagreement was only at the 1 . 5 σ level. In figure 12 we compare the BAO measurements from this work (solid blue) with those from dMdB20 using the quasar catalog of SDSS DR16 (solid gray).  \nThe differences between the analyses have been discussed in the previous sections, but we summarise them here. The main difference is the input quasar catalog. There are ∼ 150 000 quasars with z > 1 . 77 in common in both datasets, representing ∼ 45% and ∼ 25% of the quasars in SDSS DR16 and DESI DR1 respectively. DESI targeted quasars down to r < 23, compared to g < 22 in BOSS and g < 22 . 5 in eBOSS. We expect most Ly α quasars in DESI to receive four observations by the end of the survey, at which point the distribution of signal-to-noise per Angstrom in DESI and SDSS quasars will be similar, but in DESI DR1 the majority of quasars have only been observed once and are therefore on average noisier than in SDSS DR16.  \nThe spectrographs are also different: in the blue arm (relevant for Ly α science) DESI has a resolving power R in the range 2000-3000, compared to 1500-2000 in BOSS/eBOSS. While DESI pixels have a constant width of 0 . 8 ˚ A, SDSS pixels were constant in log λ with widths ranging from 0 . 83 to 1 . 27 ˚ A in the relevant range of wavelengths. The spectro-photometry is also greatly improved thanks to an atmospheric dispersion corrector, and the spectrograph optics are much more stable, allowing a finer 2D spectral extraction and improved sky subtraction [37]. The higher quality of the DESI spectra allows us to better determine quasar  \nFigure 12 . Ly α BAO measurement from DESI DR1 (blue) compared to the equivalent measurement from eBOSS using data from SDSS DR16 ( dMdB20 , gray). The red contours show the combined measurement taking into account the cross-survey covariance as discussed in appendix F.  \n<!-- image -->  \nredshifts, as captured by the σ v parameter in our fits that is almost a factor of two smaller than the one reported in dMdB20 (see section 4.7 and appendix B).  \nThe continuum fitting of quasars for both this work and dMdB20 was performed using the picca software, but the code has been re-written and re-structured in a more modular way since the analysis of dMdB20 . There have been several minor changes in the methodology as well: we now use the CIII region to re-calibrate the spectra (instead of a region redder than the MgII emission line); we extended the Ly α region to 1205 ˚ A (instead of 1200 ˚ A); we include quasars with BAL features and mask them following the prescription by [61] (instead of rejecting the entire line of sight); we use modified Ly α weights, presented in [30]. These changes were motivated by the studies of [30] using early DESI data, and the (minor) impact of these changes are discussed in the variations of section 6.4.  \nThe correlations in both analyses were measured with picca as well, and the only methodological change is the inclusion of the cross-covariance between correlations discussed in [33]. There are also some minor changes in the modelling, discussed in [31, 33]: the distortion matrix is now modelled at higher resolution and up to 300 h -1 Mpc (instead of  \n200 h -1 Mpc); the contamination of metals is modelled in a slightly different way; the model describing the contamination of HCDs now has a free parameter L HCD (fixed at 10 h -1 Mpc in dMdB20 ); the correlated sky residuals are now modelled with an improved method described in [32] that only needs a single free parameter (versus four in dMdB20 ). While these changes are a clear improvement in the modelling of the correlations, in figure 11 of section 6.4 we show that their impact on the BAO results is negligible.  \nFinally, an important difference between the analyses is the effort that we did to validate the analysis pipeline using only blinded measurements and synthetic datasets. As described in section 6, we only unblinded the measurements once we had passed a long list of robustness and consistency tests, including variations of the analysis and data splits.', '7.3 DESI - SDSS cross-survey covariance': "Even though there is a large overlap in redshift range and footprint between the Ly α samples of SDSS DR16 and DESI DR1 (see figure 2), the contribution from cosmic variance to these measurements is small, and their cross-survey covariance is smaller than one might naively think. We quantify the cross-covariance in appendix F, where we also use it to combine both results into a 'DESI + SDSS' Ly α BAO measurement:  \nα DESI+SDSS ∥ = 1 . 012 ± 0 . 016 , α DESI+SDSS ⊥ = 0 . 990 ± 0 . 019 , (7.8)  \nwith a correlation coefficient of ρ = -0 . 47. These contours are also plotted as red contours in figure 12, and can be used in cosmological inference assuming a Gaussian posterior. Given our fiducial cosmology, this gives the following constraints on ratios of distances:  \nD H ( z eff ) /r d = 8 . 72 ± 0 . 14 , D M ( z eff ) /r d = 38 . 80 ± 0 . 76 . (7.9)  \nIt is important to note that we have not redone the SDSS analysis with our analysis pipeline. We have taken the (2 × 2) Ly α BAO posterior from dMdB20 , and combined it with our own posterior after taking into account an approximate cross-survey covariance that was computed as described in appendix F. Ignoring the cross-survey covariance would lead to an underestimate of the covariance of the combined result by 10%.  \nOne could do a joint analyses of both surveys, starting by co-adding the roughly 100 000 spectra of quasars at z > 2 . 1 observed with both telescopes. However, the modelling would be complicated, since there are several differences in the contaminants of both surveys. These differences include: correlated sky residuals extend to different transverse separations; DESI observes fainter quasars, and the effective quasar bias could be different; thanks to the higher spectral resolution of DESI, quasar redshift errors are smaller than in SDSS; the efficiency of the DLA finder might also be quite different, impacting the level of HCD contamination.  \nThe combined 'DESI + SDSS' BAO measurement has an uncertainty 20% smaller than the one from DESI DR1. As the DESI survey observes more and more data, the SDSS dataset will gradually add less and less to the joint analysis.", '7.4 Future work': 'There are several known issues that we have not addressed in this publication, and that we leave for future work.  \n- · Non-linear evolution causes a small systematic shift on the BAO peak [108-110]. At low redshift, this is expected to be a sub-percent bias, and at high redshift it is expected to be even lower, below our current systematic uncertainties of 0.5% coming from the analysis of mocks (it is a conservative estimate, see section 6.5).\n- · Relative velocities between dark matter and baryons can bias the position of the BAO peak [111]. Studies using hydrodynamical simulations have shown that the bias should be small in the Ly α forest auto-correlation [112, 113], but we do not have equivalent studies for the cross-correlation with quasars.\n- · The combination of quasar redshift errors and quasar continuum errors can lead to spurious correlations in 3D analyses of the Ly α forest [114]. In a companion paper [33], we have analysed synthetic data to show that the impact of this systematic on BAO measurements is limited to 0 . 1 σ of the current statistical uncertainties. However, this could become a problem for the final analysis of DESI, or for other Ly α studies that obtain cosmological information from the full shape of the correlations [115-119].\n- · In this work we used mocks (synthetic datasets) derived from log-normal density fields. We are working on the development of more complex mocks, including mocks based on N-body simulations [100] and mocks based on Lagrangian Perturbation Theory; these mocks will include the non-linear broadening of the BAO peak discussed in figure 9.\n- · BAO analyses using galaxy samples at low redshift often use a reconstruction algorithm to undo some of the non-linear effects and improve the accuracy and precision of the BAO measurements [ ? ? ]. However, the reconstruction method is not recommended for surveys limited by shot noise (as opposed to cosmic variance), since the estimates of the displacement field would be noisy and could degrade the BAO precision. While this will be an interesting possibility for future Ly α surveys with higher quasar densities, we do not expect this method to improve the precision of the BAO measurement with the current DESI data set.', '8 Conclusions': 'We have measured the three-dimensional correlations in the Ly α forest dataset from the first year of DESI data, as well as its cross-correlation with the position of DESI quasars. Using these correlations we have measured the BAO scale parallel ( α ∥ ) and perpendicular ( α ⊥ ) to the line of sight with a precision of 2.0% and 2.4% respectively. The statistical uncertainties on the BAO parameters from just one year of DESI observations are already smaller than the ones from dMdB20 , who used 10 years of BOSS and eBOSS observations.  \nThis analysis is the first Ly α BAO measurement that was fully blinded 29 . It is also the first time that the (small) cross-covariance between the auto-correlation of Ly α fluctuations and its cross-correlation with quasars is taken into account (see section 3.3). In a companion paper [33], we validated the analysis with 150 mocks mimicking the DESI DR1 Ly α sample. We also characterized the correlated noise introduced by the data processing pipeline and the imprint of foreground absorbers on the Ly α auto-correlation in [32]. In parallel, we set a long list of robustness tests that we needed to pass before we could unblind the measurements. These are discussed in sections 6.3 and 6.4.  \nThe BAO measurement presented here can be translated into constraints on the following ratio of distances: D H ( z eff = 2 . 33) /r d = 8 . 52 ± 0 . 17 and D M ( z eff = 2 . 33) /r d = 39 . 71 ± 0 . 95, where D H = c/H ( z ), D M ( z ) is the transverse comoving distance and r d is the sound horizon at the drag epoch.  \nThis publication is part of a series of publications presenting clustering measurements of the different tracers in DESI DR1 [5]. Besides this Ly α BAO measurement at z = 2 . 33, the clustering of galaxies and quasars at z < 2 is presented in [120], the BAO measurements derived from this data set are described in [6], and the cosmological constraints from the combined galaxy, quasar and Ly α forests BAO can be found in [9]. A complementary analysis of the two points clustering statistics of galaxies and quasars using a larger range of scales beyond the BAO scale will be presented in [121] with their corresponding cosmological constraints in [122].  \nIn the near future, we will also present other cosmological analyses using the Ly α forest dataset from DESI DR1. We will extract non-BAO information from the full shape of 3D correlations on large, linear scales (see [116-119]), which should further improve the precision of these results. We will also constrain the linear matter power spectrum on small scales using the 1D power spectrum of fluctuations in the Ly α forest (see [123, 124]). DESI has completed approximately three years of observations and has collected more than half of its planned dataset. Upon completion of the five year survey, we expect the precision of the BAO results to improve by a factor of two.', 'Data Availability': 'The data used in this analysis will be made public along the Data Release 1 ( https://data. desi.lbl.gov/doc/releases/ ). The data points corresponding to the most relevant figures in this paper are available at Zenodo 30', 'Acknowledgments': "This material is based upon work supported by the U.S. Department of Energy (DOE), Office of Science, Office of High-Energy Physics, under Contract No. DE-AC02-05CH11231, and by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract. Additional support for DESI was provided by the U.S. National Science Foundation (NSF), Division of Astronomical Sciences under Contract No. AST-0950945 to the NSF's National Optical-Infrared Astronomy Research Laboratory; the Science and Technology Facilities Council of the United Kingdom; the Gordon and Betty Moore Foundation; the Heising-Simons Foundation; the French Alternative Energies and Atomic Energy Commission (CEA); the National Council of Humanities, Science and Technology of Mexico (CONAHCYT); the Ministry of Science and Innovation of Spain (MICINN), and by the DESI Member Institutions: https://www.desi.lbl.gov/ collaborating-institutions . Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the U. S. National Science Foundation, the U. S. Department of Energy, or any of the listed funding agencies.  \nThe authors are honored to be permitted to conduct scientific research on Iolkam Du'ag (Kitt Peak), a mountain with particular significance to the Tohono O'odham Nation.", 'References': '- [1] D. H. Weinberg, M. J. Mortonson, D. J. Eisenstein, C. Hirata, A. G. Riess and E. Rozo, Observational probes of cosmic acceleration , Phys. Rept. 530 (2013) 87 [ 1201.2434 ].\n- [2] Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi et al., Planck 2018 results. VI. Cosmological parameters , Astron. Astrophys. 641 (2020) A6 [ 1807.06209 ].\n- [3] S. Alam, M. Aubert, S. Avila, C. Balland, J. E. Bautista, M. A. Bershady et al., Completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: Cosmological implications from two decades of spectroscopic surveys at the Apache Point Observatory , Phys. Rev. D 103 (2021) 083533 [ 2007.08991 ].\n- [4] DESI Collaboration, A. Aghamousa, J. Aguilar, S. Ahlen, S. Alam, L. E. Allen et al., The DESI Experiment Part I: Science,Targeting, and Survey Design , arXiv e-prints (2016) arXiv:1611.00036 [ 1611.00036 ].\n- [5] DESI Collaboration, DESI 2024 I: Data Release 1 of the Dark Energy Spectroscopic Instrument , in preparation (2025) .\n- [6] DESI Collaboration, A. G. Adame, J. Aguilar, S. Ahlen, S. Alam, D. M. Alexander et al., DESI 2024 III: Baryon Acoustic Oscillations from Galaxies and Quasars , arXiv e-prints (2024) arXiv:2404.03000 [ 2404.03000 ].\n- [7] M. McQuinn, The Evolution of the Intergalactic Medium , Annual Review of Astron and Astrophys 54 (2016) 313 [ 1512.00086 ].\n- [8] A. Slosar, A. Font-Ribera, M. M. Pieri, J. Rich, J.-M. Le Goff, \' E. Aubourg et al., The Lymanα forest in three dimensions: measurements of large scale flux correlations from BOSS 1st-year data , JCAP 2011 (2011) 001 [ 1104.5244 ].\n- [9] DESI Collaboration, A. G. Adame, J. Aguilar, S. Ahlen, S. Alam, D. M. Alexander et al., DESI 2024 VI: Cosmological Constraints from the Measurements of Baryon Acoustic Oscillations , arXiv e-prints (2024) arXiv:2404.03002 [ 2404.03002 ].\n- [10] D. J. Eisenstein, I. Zehavi, D. W. Hogg, R. Scoccimarro, M. R. Blanton, R. C. Nichol et al., Detection of the Baryon Acoustic Peak in the Large-Scale Correlation Function of SDSS Luminous Red Galaxies , ApJ 633 (2005) 560 [ astro-ph/0501171 ].\n- [11] S. Cole, W. J. Percival, J. A. Peacock, P. Norberg, C. M. Baugh, C. S. Frenk et al., The 2dF Galaxy Redshift Survey: power-spectrum analysis of the final data set and cosmological implications , Mon. Not. Roy. Astron. Soc. 362 (2005) 505 [ astro-ph/0501174 ].\n- [12] K. S. Dawson, D. J. Schlegel, C. P. Ahn, S. F. Anderson, \' E. Aubourg, S. Bailey et al., The Baryon Oscillation Spectroscopic Survey of SDSS-III , AJ 145 (2013) 10 [ 1208.0022 ].\n- [13] N. G. Busca, T. Delubac, J. Rich, S. Bailey, A. Font-Ribera, D. Kirkby et al., Baryon acoustic oscillations in the Ly α forest of BOSS quasars , Astron. Astrophys. 552 (2013) A96 [ 1211.2616 ].\n- [14] A. Slosar, V. Irˇsiˇc, D. Kirkby, S. Bailey, N. G. Busca, T. Delubac et al., Measurement of baryon acoustic oscillations in the Lymanα forest fluctuations in BOSS data release 9 , JCAP 2013 (2013) 026 [ 1301.3459 ].\n- [15] D. Kirkby, D. Margala, A. Slosar, S. Bailey, N. G. Busca, T. Delubac et al., Fitting methods for baryon acoustic oscillations in the Lymanα forest fluctuations in BOSS data release 9 , JCAP 3 (2013) 024 [ 1301.3456 ].\n- [16] D. J. Eisenstein, D. H. Weinberg, E. Agol, H. Aihara, C. Allende Prieto, S. F. Anderson et al., SDSS-III: Massive Spectroscopic Surveys of the Distant Universe, the Milky Way, and Extra-Solar Planetary Systems , AJ 142 (2011) 72 [ 1101.1529 ].  \n- [17] C. P. Ahn, R. Alexandroff, C. Allende Prieto, S. F. Anderson, T. Anderton, B. H. Andrews et al., The Ninth Data Release of the Sloan Digital Sky Survey: First Spectroscopic Data from the SDSS-III Baryon Oscillation Spectroscopic Survey , ApJS 203 (2012) 21 [ 1207.7137 ].\n- [18] S. Alam, F. D. Albareti, C. Allende Prieto, F. Anders, S. F. Anderson, T. Anderton et al., The Eleventh and Twelfth Data Releases of the Sloan Digital Sky Survey: Final Data from SDSS-III , ApJS 219 (2015) 12 [ 1501.00963 ].\n- [19] A. Font-Ribera, D. Kirkby, N. Busca, J. Miralda-Escud\'e, N. P. Ross, A. Slosar et al., Quasar-Lyman α forest cross-correlation from BOSS DR11: Baryon Acoustic Oscillations , JCAP 5 (2014) 27 [ 1311.1767 ].\n- [20] T. Delubac, J. E. Bautista, N. G. Busca, J. Rich, D. Kirkby, S. Bailey et al., Baryon acoustic oscillations in the Ly α forest of BOSS DR11 quasars , Astron. Astrophys. 574 (2015) A59 [ 1404.1801 ].\n- [21] J. E. Bautista, N. G. Busca, J. Guy, J. Rich, M. Blomqvist, H. du Mas des Bourboux et al., Measurement of baryon acoustic oscillation correlations at z = 2.3 with SDSS DR12 Ly α -Forests , Astron. Astrophys. 603 (2017) A12 [ 1702.00176 ].\n- [22] H. du Mas des Bourboux, J.-M. Le Goff, M. Blomqvist, N. G. Busca, J. Guy, J. Rich et al., Baryon acoustic oscillations from the complete SDSS-III Ly α -quasar cross-correlation function at z = 2.4 , Astron. Astrophys. 608 (2017) A130 [ 1708.02225 ].\n- [23] V. de Sainte Agathe, C. Balland, H. du Mas des Bourboux, N. G. Busca, M. Blomqvist, J. Guy et al., Baryon acoustic oscillations at z = 2.34 from the correlations of Ly α absorption in eBOSS DR14 , Astron. Astrophys. 629 (2019) A85 [ 1904.03400 ].\n- [24] M. Blomqvist, H. du Mas des Bourboux, N. G. Busca, V. de Sainte Agathe, J. Rich, C. Balland et al., Baryon acoustic oscillations from the cross-correlation of Ly α absorption and quasars in eBOSS DR14 , Astron. Astrophys. 629 (2019) A86 [ 1904.03430 ].\n- [25] K. S. Dawson, J.-P. Kneib, W. J. Percival, S. Alam, F. D. Albareti, S. F. Anderson et al., The SDSS-IV Extended Baryon Oscillation Spectroscopic Survey: Overview and Early Data , AJ 151 (2016) 44 [ 1508.04473 ].\n- [26] R. Ahumada, C. Allende Prieto, A. Almeida, F. Anders, S. F. Anderson, B. H. Andrews et al., The 16th Data Release of the Sloan Digital Sky Surveys: First Release from the APOGEE-2 Southern Survey and Full Release of eBOSS Spectra , ApJS 249 (2020) 3 [ 1912.02905 ].\n- [27] H. du Mas des Bourboux, J. Rich, A. Font-Ribera, V. de Sainte Agathe, J. Farr, T. Etourneau et al., The Completed SDSS-IV Extended Baryon Oscillation Spectroscopic Survey: Baryon Acoustic Oscillations with Ly α Forests , ApJ 901 (2020) 153 [ 2007.08995 ].\n- [28] E. Chaussidon, C. Y\'eche, N. Palanque-Delabrouille, D. M. Alexander, J. Yang, S. Ahlen et al., Target Selection and Validation of DESI Quasars , ApJ 944 (2023) 107 [ 2208.08511 ].\n- [29] D. M. Alexander, T. M. Davis, E. Chaussidon, V. A. Fawcett, A. X. Gonzalez-Morales, T.-W. Lan et al., The DESI Survey Validation: Results from Visual Inspection of the Quasar Survey Spectra , AJ 165 (2023) 124 [ 2208.08517 ].\n- [30] C. Ram\'ırez-P\'erez, I. P\'erez-R\'afols, A. Font-Ribera, M. A. Karim, E. Armengaud, J. Bautista et al., The Lymanα forest catalogue from the Dark Energy Spectroscopic Instrument Early Data Release , Mon. Not. Roy. Astron. Soc. 528 (2024) 6666 [ 2306.06312 ].\n- [31] C. Gordon, A. Cuceu, J. Chaves-Montero, A. Font-Ribera, A. X. Gonz\'alez-Morales, J. Aguilar et al., 3D correlations in the Lymanα forest from early DESI data , JCAP 2023 (2023) 045 [ 2308.10950 ].\n- [32] J. Guy, S. G. A. Gontcho, E. Armengaud, A. Brodzeller, A. Cuceu, A. Font-Ribera et al., Characterization of contaminants in the Lyman-alpha forest auto-correlation with DESI , arXiv e-prints (2024) arXiv:2404.03003 [ 2404.03003 ].  \n- [33] A. Cuceu, H. K. Herrera-Alcantar, C. Gordon, P. Martini, J. Guy, A. Font-Ribera et al., Validation of the DESI 2024 Ly α forest BAO analysis using synthetic datasets , arXiv e-prints (2024) arXiv:2404.03004 [ 2404.03004 ].\n- [34] H. K. Herrera-Alcantar, A. Mu˜noz-Guti\'errez, T. Tan, A. X. Gonz\'alez-Morales, A. Font-Ribera, J. Guy et al., Synthetic spectra for Lymanα forest analysis in the Dark Energy Spectroscopic Instrument , arXiv e-prints (2023) arXiv:2401.00303 [ 2401.00303 ].\n- [35] T. N. Miller, P. Doel, G. Gutierrez, R. Besuner, D. Brooks, G. Gallo et al., The Optical Corrector for the Dark Energy Spectroscopic Instrument , arXiv e-prints (2023) arXiv:2306.06310 [ 2306.06310 ].\n- [36] J. H. Silber, P. Fagrelius, K. Fanning, M. Schubnell, J. N. Aguilar, S. Ahlen et al., The Robotic Multiobject Focal Plane System of the Dark Energy Spectroscopic Instrument (DESI) , AJ 165 (2023) 9 [ 2205.09014 ].\n- [37] J. Guy, S. Bailey, A. Kremin, S. Alam, D. M. Alexander, C. Allende Prieto et al., The Spectroscopic Data Processing Pipeline for the Dark Energy Spectroscopic Instrument , AJ 165 (2023) 144 [ 2209.14482 ].\n- [38] DESI Collaboration, B. Abareshi, J. Aguilar, S. Ahlen, S. Alam, D. M. Alexander et al., Overview of the Instrumentation for the Dark Energy Spectroscopic Instrument , AJ 164 (2022) 207 [ 2205.10939 ].\n- [39] DESI Collaboration, A. G. Adame, J. Aguilar, S. Ahlen, S. Alam, G. Aldering et al., Validation of the Scientific Program for the Dark Energy Spectroscopic Instrument , AJ 167 (2024) 62 [ 2306.06307 ].\n- [40] A. D. Myers, J. Moustakas, S. Bailey, B. A. Weaver, A. P. Cooper, J. E. Forero-Romero et al., The Target-selection Pipeline for the Dark Energy Spectroscopic Instrument , AJ 165 (2023) 50 [ 2208.08518 ].\n- [41] E. F. Schlafly, D. Kirkby, D. J. Schlegel, A. D. Myers, A. Raichoor, K. Dawson et al., Survey Operations for the Dark Energy Spectroscopic Instrument , AJ 166 (2023) 259 [ 2306.06309 ].\n- [42] K. M. G\'orski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke et al., HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere , ApJ 622 (2005) 759 [ astro-ph/0409513 ].\n- [43] C. Y\'eche, N. Palanque-Delabrouille, C.-A. Claveau, D. D. Brooks, E. Chaussidon, T. M. Davis et al., Preliminary Target Selection for the DESI Quasar (QSO) Sample , Research Notes of the American Astronomical Society 4 (2020) 179 [ 2010.11280 ].\n- [44] Bailey et al., Redrock: Spectroscopic classification and redshift fitting for the dark energy spectroscopic instrument , in preparation (2024) .\n- [45] N. Busca and C. Balland, QuasarNET: Human-level spectral classification and redshifting with Deep Neural Networks , arXiv e-prints (2018) arXiv:1808.09955 [ 1808.09955 ].\n- [46] J. Farr, A. Font-Ribera and A. Pontzen, Optimal strategies for identifying quasars in DESI , JCAP 2020 (2020) 015 [ 2007.10348 ].\n- [47] A. Brodzeller, K. Dawson, S. Bailey, J. Yu, A. J. Ross, A. Bault et al., Performance of the Quasar Spectral Templates for the Dark Energy Spectroscopic Instrument , AJ 166 (2023) 66 [ 2305.10426 ].\n- [48] V. Kamble, K. Dawson, H. du Mas des Bourboux, J. Bautista and D. P. Scheinder, Measurements of Effective Optical Depth in the Ly α Forest from the BOSS DR12 Quasar Sample , ApJ 892 (2020) 70 [ 1904.01110 ].\n- [49] A. Bault, D. Kirkby, J. Guy, A. Brodzeller, J. Aguilar, S. Ahlen et al., Impact of Systematic Redshift Errors on the Cross-correlation of the Lymanα Forest with Quasars at Small Scales Using DESI Early Data , arXiv e-prints (2024) arXiv:2402.18009 [ 2402.18009 ].  \n- [50] D. Parks, J. X. Prochaska, S. Dong and Z. Cai, Deep learning of quasar spectra to discover and characterize damped Ly α systems , Mon. Not. Roy. Astron. Soc. 476 (2018) 1151 [ 1709.04962 ].\n- [51] M.-F. Ho, S. Bird and R. Garnett, Detecting multiple DLAs per spectrum in SDSS DR12 with Gaussian processes , Mon. Not. Roy. Astron. Soc. 496 (2020) 5436 [ 2003.11036 ].\n- [52] B. Wang, J. Zou, Z. Cai, J. X. Prochaska, Z. Sun, J. Ding et al., Deep Learning of Dark Energy Spectroscopic Instrument Mock Spectra to Find Damped Ly α Systems , ApJS 259 (2022) 28.\n- [53] C. B. Foltz, F. H. Chaffee, P. C. Hewett, R. J. Weymann and S. L. Morris, On the Fraction of Optically-Selected QSOs with Broad Absorption Lines in Their Spectra , in Bulletin of the American Astronomical Society , vol. 22, p. 806, Mar., 1990.\n- [54] J. R. Trump, P. B. Hall, T. A. Reichard, G. T. Richards, D. P. Schneider, D. E. Vanden Berk et al., A Catalog of Broad Absorption Line Quasars from the Sloan Digital Sky Survey Third Data Release , ApJS 165 (2006) 1 [ astro-ph/0603070 ].\n- [55] L. Mas-Ribas and R. Mauland, The ubiquitous imprint of radiative acceleration in the mean absorption spectrum of quasar outflows , The Astrophysical Journal 886 (2019) 151.\n- [56] L. \' A. Garc\'ıa, P. Martini, A. X. Gonzalez-Morales, A. Font-Ribera, H. K. Herrera-Alcantar, J. N. Aguilar et al., Analysis of the impact of broad absorption lines on quasar redshift measurements with synthetic observations , Mon. Not. Roy. Astron. Soc. 526 (2023) 4848 [ 2304.05855 ].\n- [57] S. Filbert, P. Martini, K. Seebaluck, L. Ennesser, D. M. Alexander, A. Bault et al., Broad Absorption Line Quasars in the Dark Energy Spectroscopic Instrument Early Data Release , arXiv e-prints (2023) arXiv:2309.03434 [ 2309.03434 ].\n- [58] P. B. Hall, S. F. Anderson, M. A. Strauss, D. G. York, G. T. Richards, X. Fan et al., Unusual Broad Absorption Line Quasars from the Sloan Digital Sky Survey , ApJS 141 (2002) 267 [ astro-ph/0203252 ].\n- [59] R. J. Weymann, S. L. Morris, C. B. Foltz and P. C. Hewett, Comparisons of the emission-line and continuum properties of broad absorption line and normal quasi-stellar objects , Astrophysical Journal 373 (1991) 23.\n- [60] A. Brodzeller and K. Dawson, Modeling the Spectral Diversity of Quasars in the Sixteenth Data Release from the Sloan Digital Sky Survey , AJ 163 (2022) 110 [ 2110.07748 ].\n- [61] L. Ennesser, P. Martini, A. Font-Ribera and I. P\'erez-R\'afols, The impact and mitigation of broad-absorption-line quasars in Lyman α forest correlations , Mon. Not. Roy. Astron. Soc. 511 (2022) 3514 [ 2111.09439 ].\n- [62] P. Martini et al., Validation of the DESI 2024 Lyman Alpha Forest BAL Masking Strategy , in preparation (2024) .\n- [63] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau et al., Array programming with NumPy , Nature 585 (2020) 357.\n- [64] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau et al., SciPy 1.0: fundamental algorithms for scientific computing in Python , Nature Methods 17 (2020) 261 [ 1907.10121 ].\n- [65] Astropy Collaboration, T. P. Robitaille, E. J. Tollerud, P. Greenfield, M. Droettboom, E. Bray et al., Astropy: A community Python package for astronomy , Astron. Astrophys. 558 (2013) A33 [ 1307.6212 ].\n- [66] Astropy Collaboration, A. M. Price-Whelan, B. M. Sip"ocz, H. M. Gunther, P. L. Lim, S. M. Crawford et al., The Astropy Project: Building an Open-science Project and Status of the v2.0 Core Package , AJ 156 (2018) 123 [ 1801.02634 ].  \n- [67] Astropy Collaboration, A. M. Price-Whelan, P. L. Lim, N. Earl, N. Starkman, L. Bradley et al., The Astropy Project: Sustaining and Growing a Community-oriented Open-source Project and the Latest Major Release (v5.0) of the Core Package , ApJ 935 (2022) 167 [ 2206.14220 ].\n- [68] L. Dalcin and Y.-L. L. Fang, mpi4py: Status update after 12 years of development , Computing in Science & Engineering 23 (2021) 47.\n- [69] A. Zonca, L. Singer, D. Lenz, M. Reinecke, C. Rosset, E. Hivon et al., healpy: equal area pixelization and spherical harmonics transforms for data on the sphere in python , Journal of Open Source Software 4 (2019) 1298.\n- [70] J. D. Hunter, Matplotlib: A 2d graphics environment , Computing in Science & Engineering 9 (2007) 90.\n- [71] A. Lewis, GetDist: a Python package for analysing Monte Carlo samples , arXiv e-prints (2019) arXiv:1910.13970 [ 1910.13970 ].\n- [72] S. K. Lam, A. Pitrou and S. Seibert, Numba: A LLVM-based Python JIT Compiler , in Proc. Second Workshop on the LLVM Compiler Infrastructure in HPC , pp. 1-6, Nov., 2015, DOI.\n- [73] H. du Mas des Bourboux, J. Rich, A. Font-Ribera, V. de Sainte Agathe, J. Farr, T. Etourneau et al., \'picca: Package for Igm Cosmological-Correlations Analyses.\' Astrophysics Source Code Library, record ascl:2106.018, June, 2021.\n- [74] A. Cuceu, A. Font-Ribera and B. Joachimi, Bayesian methods for fitting Baryon Acoustic Oscillations in the Lymanα forest , JCAP 2020 (2020) 035 [ 2004.02761 ].\n- [75] P. McDonald, U. Seljak, S. Burles, D. J. Schlegel, D. H. Weinberg, R. Cen et al., The Ly α Forest Power Spectrum from the Sloan Digital Sky Survey , ApJS 163 (2006) 80 [ astro-ph/0405013 ].\n- [76] A. Font-Ribera, J. Miralda-Escud\'e, E. Arnau, B. Carithers, K.-G. Lee, P. Noterdaeme et al., The large-scale cross-correlation of Damped Lyman alpha systems with the Lyman alpha forest: first measurements from BOSS , JCAP 2012 (2012) 059 [ 1209.4596 ].\n- [77] H. du Mas des Bourboux, K. S. Dawson, N. G. Busca, M. Blomqvist, V. de Sainte Agathe, C. Balland et al., The Extended Baryon Oscillation Spectroscopic Survey: Measuring the Cross-correlation between the Mg II Flux Transmission Field and Quasars and Galaxies at z = 0.59 , ApJ 878 (2019) 47 [ 1901.01950 ].\n- [78] N. Kaiser, Clustering in real space and in redshift space , Mon. Not. Roy. Astron. Soc. 227 (1987) 1.\n- [79] P. McDonald, J. Miralda-Escud\'e, M. Rauch, W. L. W. Sargent, T. A. Barlow, R. Cen et al., The Observed Probability Distribution Function, Power Spectrum, and Correlation Function of the Transmitted Flux in the Ly α Forest , ApJ 543 (2000) 1 [ astro-ph/9911196 ].\n- [80] P. McDonald, Toward a Measurement of the Cosmological Geometry at z ˜ 2: Predicting Ly α Forest Correlation in Three Dimensions and the Potential of Future Data Sets , ApJ 585 (2003) 34 [ astro-ph/0108064 ].\n- [81] M. M. Pieri, M. J. Mortonson, S. Frank, N. Crighton, D. H. Weinberg, K.-G. Lee et al., Probing the circumgalactic medium at high-redshift using composite BOSS spectra of strong Lyman α forest absorbers , Mon. Not. Roy. Astron. Soc. 441 (2014) 1718 [ 1309.6768 ].\n- [82] L. Yang, Z. Zheng, H. du Mas des Bourboux, K. Dawson, M. M. Pieri, G. Rossi et al., Metal Lines Associated with the Ly α Forest from eBOSS Data , ApJ 935 (2022) 121 [ 2206.11385 ].\n- [83] S. Morrison, D. Som, M. M. Pieri, I. P\'erez-R\'afols and M. Blomqvist, A Strong Blend in the Morning: Studying the Circumgalactic Medium Before Cosmic Noon with Strong, Blended Lymanα Forest Systems , arXiv e-prints (2023) arXiv:2309.06813 [ 2309.06813 ].  \n- [84] M. McQuinn and M. White, On estimating Ly α forest correlations between multiple sightlines , Mon. Not. Roy. Astron. Soc. 415 (2011) 2257 [ 1102.1752 ].\n- [85] A. Font-Ribera and J. Miralda-Escud\'e, The effect of high column density systems on the measurement of the Lymanα forest correlation function , JCAP 2012 (2012) 028 [ 1205.2018 ].\n- [86] K. K. Rogers, S. Bird, H. V. Peiris, A. Pontzen, A. Font-Ribera and B. Leistedt, Correlations in the three-dimensional Lyman-alpha forest contaminated by high column density absorbers , Mon. Not. Roy. Astron. Soc. 476 (2018) 3716 [ 1711.06275 ].\n- [87] T. Tan et al., Modeling of the high column density systems in the lyman-alpha forest , in preparation (2024) .\n- [88] I. P\'erez-R\'afols, A. Font-Ribera, J. Miralda-Escud\'e, M. Blomqvist, S. Bird, N. Busca et al., The SDSS-DR12 large-scale cross-correlation of damped Lyman alpha systems with the Lyman alpha forest , Mon. Not. Roy. Astron. Soc. 473 (2018) 3019 [ 1709.00889 ].\n- [89] I. P\'erez-R\'afols, M. M. Pieri, M. Blomqvist, S. Morrison, D. Som and A. Cuceu, The cross-correlation of galaxies in absorption with the Lyman α forest , Mon. Not. Roy. Astron. Soc. 524 (2023) 1464 [ 2210.02973 ].\n- [90] A. Font-Ribera, E. Arnau, J. Miralda-Escud\'e, E. Rollinde, J. Brinkmann, J. R. Brownstein et al., The large-scale quasar-Lyman α forest cross-correlation from BOSS , JCAP 5 (2013) 018 [ 1303.1937 ].\n- [91] A. Arinyo-i-Prats, J. Miralda-Escud\'e, M. Viel and R. Cen, The non-linear power spectrum of the Lyman alpha forest , JCAP 2015 (2015) 017 [ 1506.04519 ].\n- [92] J. J. Givans, A. Font-Ribera, A. Slosar, L. Seeyave, C. Pedersen, K. K. Rogers et al., Non-linearities in the Lymanα forest and in its cross-correlation with dark matter halos , JCAP 2022 (2022) 070 [ 2205.00962 ].\n- [93] W. J. Handley, M. P. Hobson and A. N. Lasenby, polychord: nested sampling for cosmology. , Mon. Not. Roy. Astron. Soc. 450 (2015) L61 [ 1502.01856 ].\n- [94] W. J. Handley, M. P. Hobson and A. N. Lasenby, POLYCHORD: next-generation nested sampling , Mon. Not. Roy. Astron. Soc. 453 (2015) 4384 [ 1506.00171 ].\n- [95] H. Dembinski and P. O. et al., scikit-hep/iminuit , .\n- [96] F. James and M. Roos, Minuit: A System for Function Minimization and Analysis of the Parameter Errors and Correlations , Comput. Phys. Commun. 10 (1975) 343.\n- [97] J. Farr, A. Font-Ribera, H. du Mas des Bourboux, A. Mu˜noz-Guti\'errez, F. J. S\'anchez, A. Pontzen et al., LyaCoLoRe: synthetic datasets for current and future Lymanα forest BAO surveys , JCAP 2020 (2020) 068 [ 1912.02763 ].\n- [98] C. Ram\'ırez-P\'erez, J. Sanchez, D. Alonso and A. Font-Ribera, CoLoRe: fast cosmological realisations over large volumes with multiple tracers , JCAP 2022 (2022) 002 [ 2111.05069 ].\n- [99] T. Etourneau, J.-M. Le Goff, J. Rich, T. Tan, A. Cuceu, S. Ahlen et al., Mock data sets for the Eboss and DESI Lymanα forest surveys , arXiv e-prints (2023) arXiv:2310.18996 [ 2310.18996 ].\n- [100] B. Hadzhiyska et al., Planting a Lyman alpha forest on AbacusSummit , Mon. Not. Roy. Astron. Soc. 524 (2023) 1008 [ 2305.08899 ].\n- [101] A. Dey, D. J. Schlegel, D. Lang, R. Blum, K. Burleigh, X. Fan et al., Overview of the DESI Legacy Imaging Surveys , AJ 157 (2019) 168 [ 1804.08657 ].\n- [102] H. Zou, X. Zhou, X. Fan, T. Zhang, Z. Zhou, J. Nie et al., Project Overview of the Beijing-Arizona Sky Survey , PASP 129 (2017) 064101 [ 1702.03653 ].\n- [103] Schlegel et al., Desi legacy imaging surveys data release 9 , in preparation (2024) .  \n- [104] J. A. Baldwin, Luminosity Indicators in the Spectra of Quasi-Stellar Objects , ApJ 214 (1977) 679.\n- [105] A. Pontzen, Scale-dependent bias in the baryonic-acoustic-oscillation-scale intergalactic neutral hydrogen , Phys. Rev. D 89 (2014) 083010 [ 1402.0506 ].\n- [106] S. Gontcho A Gontcho, J. Miralda-Escud\'e and N. G. Busca, On the effect of the ionizing background on the Ly α forest autocorrelation function , Mon. Not. Roy. Astron. Soc. 442 (2014) 187 [ 1404.7425 ].\n- [107] C. Alcock and B. Paczynski, An evolution free test for non-zero cosmological constant , Nature 281 (1979) 358.\n- [108] M. Crocce and R. Scoccimarro, Nonlinear evolution of baryon acoustic oscillations , Phys. Rev. D 77 (2008) 023533 [ 0704.2783 ].\n- [109] R. E. Smith, R. Scoccimarro and R. K. Sheth, Motion of the acoustic peak in the correlation function , Phys. Rev. D 77 (2008) 043525 [ astro-ph/0703620 ].\n- [110] H.-J. Seo, E. R. Siegel, D. J. Eisenstein and M. White, Nonlinear Structure Formation and the Acoustic Scale , ApJ 686 (2008) 13 [ 0805.0117 ].\n- [111] D. Tseliakhovich and C. Hirata, Relative velocity of dark matter and baryonic fluids and the formation of the first structures , Phys. Rev. D 82 (2010) 083520 [ 1005.2416 ].\n- [112] C. M. Hirata, Small-scale structure and the Lymanα forest baryon acoustic oscillation feature , Mon. Not. Roy. Astron. Soc. 474 (2018) 2173 [ 1707.03358 ].\n- [113] J. J. Givans and C. M. Hirata, Redshift-space streaming velocity effects on the Lymanα forest baryon acoustic oscillation scale , Phys. Rev. D 102 (2020) 023515 [ 2002.12296 ].\n- [114] S. Youles, J. E. Bautista, A. Font-Ribera, D. Bacon, J. Rich, D. Brooks et al., The effect of quasar redshift errors on Lymanα forest correlation functions , Mon. Not. Roy. Astron. Soc. 516 (2022) 421 [ 2205.06648 ].\n- [115] A. Font-Ribera, P. McDonald and A. Slosar, How to estimate the 3D power spectrum of the Lymanα forest , JCAP 2018 (2018) 003 [ 1710.11036 ].\n- [116] A. Cuceu, A. Font-Ribera, B. Joachimi and S. Nadathur, Cosmology beyond BAO from the 3D distribution of the Lymanα forest , Mon. Not. Roy. Astron. Soc. 506 (2021) 5439 [ 2103.14075 ].\n- [117] A. Cuceu, A. Font-Ribera, S. Nadathur, B. Joachimi and P. Martini, Constraints on the Cosmic Expansion Rate at Redshift 2.3 from the Lymanα Forest , Phys. Rev. Lett. 130 (2023) 191003 [ 2209.13942 ].\n- [118] A. Cuceu, A. Font-Ribera, P. Martini, B. Joachimi, S. Nadathur, J. Rich et al., The Alcock-Paczy\'nski effect from Lymanα forest correlations: analysis validation with synthetic data , Mon. Not. Roy. Astron. Soc. 523 (2023) 3773 [ 2209.12931 ].\n- [119] F. Gerardi, A. Cuceu, A. Font-Ribera, B. Joachimi and P. Lemos, Direct cosmological inference from three-dimensional correlations of the Lyman α forest , Mon. Not. Roy. Astron. Soc. 518 (2023) 2567 [ 2209.11263 ].\n- [120] DESI Collaboration, DESI 2024 II: Sample definitions, characteristics and two-point clustering statistics , in preparation (2024) .\n- [121] DESI Collaboration, DESI 2024 V: Analysis of the full shape of two-point clustering statistics from galaxies and quasars , in preparation (2024) .\n- [122] DESI Collaboration, DESI 2024 VII: Cosmological constraints from full-shape analyses of the two-point clustering statistics measurements , in preparation (2024) .\n- [123] C. Ravoux, M. L. Abdul Karim, E. Armengaud, M. Walther, N. G. Kara¸caylı, P. Martini et al., The Dark Energy Spectroscopic Instrument: one-dimensional power spectrum from first  \nLy α forest samples with Fast Fourier Transform , Mon. Not. Roy. Astron. Soc. 526 (2023) 5118 [ 2306.06311 ].  \n[124] N. G. Kara¸caylı, P. Martini, J. Guy, C. Ravoux, M. L. A. Karim, E. Armengaud et al., Optimal 1D Ly α Forest Power Spectrum Estimation - III. DESI early data , Mon. Not. Roy. Astron. Soc. (2024) [ 2306.06316 ].  \n[125] J. X. Prochaska, P. Madau, J. M. O\'Meara and M. Fumagalli, Towards a unified description of the intergalactic medium at redshift z ≈ 2.5 , Mon. Not. Roy. Astron. Soc. 438 (2014) 476 [ 1310.0052 ].  \n[126] E. Chaussidon et al., Constraints on Primordial Non-Gaussianities for DESI 2024 , in preparation (2024) .  \n[127] S. Gratton and A. Challinor, Understanding parameter differences between analyses employing nested data subsets , Mon. Not. Roy. Astron. Soc. 499 (2020) 3410 [ 1911.07754 ].', 'A Alternative modelling of HCD contaminations': 'As discussed in section 4.6, the presence of high column density systems (HCDs) changes the flux correlation function by smearing the δ field in the radial direction [84-86]. The contamination can be described by introducing a scale-dependent component to the bias that depends on k ∥ that is added to the normal IGM-induced Ly α bias. Following [21] we write this term as b HCD F HCD ( k ∥ ) where  \nF HCD ( k ∥ , z ) = A ( z ) ∫ dN HI f ( N HI , z ) V ( N HI , k ∥ , z ) (A.1)  \nwhere V ( N HI , k ∥ ) is the Fourier transform of the Voigt profile for an HCD of column density N HI , and f ( N HI ) is the column-density distribution of HCDs. The normalization factor A ( z ) can be chosen so that F HCD ( k ∥ = 0) = 1, in which case b HCD is proportional to the product of the HCD halo bias and the mean absorption caused by HCDs (see appendix B in [84] and eq. 4.19 in [85]).  \nGiven our lack of precise knowledge of the HCD distribution f ( N HI ), following dMdB20 we model F HCD ( k ∥ ) = exp( -L HCD k ∥ ) as an exponential with unknown scale parameter L HCD that characterizes the typical size of unmasked HCDs.  \nWe compare the functional forms of F HCD ( k ∥ ) in figure 13. The solid orange line shows the computation from equation (A.1) when we use the column density distribution f ( N HI ) from [125]. The solid blue line uses the same model, but it only integrates up to log( N HI ) = 20 . 3 to mimic the effect of perfectly masking all DLAs. The dashed red and green lines show the exponential model for L HCD = 7 and L HCD = 3 h -1 Mpc, respectively, and they capture fairly well the suppression of power from equation (A.1).', 'B Nuisance parameters': 'As discussed in section 4, our model to describe the measured correlations has 17 free parameters, including the two BAO parameters ( α ⊥ , α ∥ ) and 15 nuisance parameters that we marginalize over. The best-fit values and uncertainties for all 17 parameters are shown in table 5, together with the priors used in the analysis.  \nWhile most of the parameters are well constrained by the data and the priors are uninformative, we use informative Gaussian priors for two parameters. These are the linear  \nFigure 13 . Comparison of the Voigt model for HCD contamination (solid curves) and the exponential model used in this work (dashed curves), for the HCD column-density distribution predicted by [125]. The orange solid line includes HCDs of all column densities and the blue solid line includes only those with log( N HI ) < 20 . 3, which is appropriate for perfect masking of DLAs. The dashed lines show that exponential models can reproduce the Voigt model, except at large k ∥ where the Voigt model has a longer tail.  \n<!-- image -->  \nbias of CIV absorption, b CIV(eff) , and the RSD parameter of absorption caused by HCDs, β HCD . When fitting the auto- and cross-correlation independently, we are not able to break the degeneracy between L HCD and other nuisance parameters, and we fix this value to the best-fit value of the combined analysis ( L HCD = 6 . 51 h -1 Mpc). Moreover, when fitting the cross-correlation alone we also need to add an informative Gaussian prior on the quasar bias 31 of b Q = 3 . 5 ± 0 . 1 to break the degeneracy between this parameter and the bias of the Ly α forest ( b α ).  \nComparison of the different columns in table 5 shows that the best-fit values from the auto-correlations alone (including regions A and B) are in agreement with those from the fit of the cross-correlations with quasars (also including both regions). However, we would like to add a word of caution when interpreting these nuisance parameters. While we have extensively tested the robustness of the BAO results under different data splits and analysis settings, some of the nuisance parameters do vary significantly with reasonable changes in the analysis choices. For instance, depending on how aggressively we mask DLAs, we obtain different values for the parameters that model the contamination by HCDs ( b HCD , β HCD and L HCD ), but the differences also propagate to other parameters that are degenerate with these, including b α and β α .  \nFinally, in figure 14 we show that none of the 15 nuisance parameters is correlated with any of the BAO parameters ( α ⊥ , α ∥ ).  \nTable 5 . Priors, best-fit values (mean of the posterior) and uncertainties (68% credible intervals) for the 17 free parameters in the fits. When analysing the auto-correlation or the cross-correlation alone, we fix L HCD to the best-fit value of the combination ( L HCD = 6 . 51 h -1 Mpc). This is necessary to break internal degeneracies, but it makes the p-value of these analyses difficult to interpret. When analysing the cross-correlation alone we also use an extra prior on the quasar bias parameter ( b Q = 3 . 5 ± 0 . 1) to break the degeneracy with the Ly α biases. Some parameters are not needed when fitting the auto-correlation or the cross-correlation alone.', 'C Blinding': 'Our analysis validation was performed entirely on blinded data, with clearly defined requirements that needed to be achieved before unblinding (Section 6). The main goal of our blinding strategy was to blind the BAO measurement in a way that does not impede the analysis process. We considered multiple blinding methods, including blinding at the catalog level (as done for the DESI galaxy BAO analyses, see [6]). However, the disadvantage of catalog-level blinding is that metal contamination results in very well-measured peaks along the line-ofsight due to Ly α × Metal correlations (see section 4.4). Any catalogue-level blinding would have also shifted the positions of these peaks, giving away the direction and magnitude of the blinding. We therefore instead developed a blinding method that only shifts the position  \nFigure 14 . Correlation between the BAO parameters ( α ⊥ , α ∥ ) and the 15 nuisance parameters, for the combined analysis (solid red), the cross-correlation alone (dashed black) and the auto-correlation alone (filled blue). Not all nuisance parameters are varied when fitting the auto- and cross-correlations separately. The BAO parameters are not strongly correlated with any of the nuisance parameters.  \n<!-- image -->  \nof the BAO peak at the level of the measured correlation functions.  \nOur blinding strategy starts with a correlation function model ξ t with all nuisance  \nFigure 15 . BAO parameters along ( α ∥ ) and across the line of sight ( α ⊥ ) from DESI DR1 (solid red) and from its blinded measurement (dashed red). The blinded measurement was in tension with Planck (dashed gray) and with the results from dMdB20 using SDSS DR16 (solid blue).  \n<!-- image -->  \nparameters set to their best-fit values from [27]. We then used Vega to compute the blinding template:  \nξ b = ξ t ( α || = 1 + ∆ α || , α ⊥ = 1 + ∆ α ⊥ ) -ξ t ( α || = 1 , α ⊥ = 1) , (C.1)  \nwhere ∆ α || and ∆ α ⊥ were randomly drawn from Gaussian distributions with mean zero and variance equal to two times the uncertainties on α || and α ⊥ from [27]. The random values used in the DR1 BAO analysis were ∆ α || = -0 . 011 and ∆ α ⊥ = -0 . 084, i.e., a fairly small shift along the line of sight and a very large blinding shift in the transverse BAO measurement, close to the maximum allowed by the blinding strategy. These values of ∆ α || and ∆ α ⊥ were unknown to us and never stored anywhere. We did save the template 32 and our pipeline ( picca ) automatically added this ξ b template to the correlation function at the moment of writing it to file, effectively blinding the BAO scale.  \nThe random shift applied to the blinded measurements can be seen in figure 15, where we show the DESI Ly α BAO results before (dotted red) and after (solid red) unblinding. The difference in the best fit BAO parameters before and after unblinding (∆ α || = -0 . 005,  \nFigure 16 . 68% and 95% credible regions of BAO parameters along the line-of-sight ( α ∥ ) and across the line-of-sight ( α ⊥ ). The filled blue contours are based on the full-posterior distribution computed by Polychord and the solid red contours are based on the approximate Gaussian fit computed by iminuit .  \n<!-- image -->  \n∆ α ⊥ = -0 . 098, red arrow in the figure) is consistent with the random shifts applied in the blinding (∆ α || = -0 . 011, ∆ α ⊥ = -0 . 084).  \nWe unblinded our analysis on December 8th, 2023, after passing the extensive validation tests described in Section 6. We only made two minor changes in the methodology after unblinding. First, we fixed a small bug in picca related to the masking of BAL features in the Ly α region B, with an impact on the best-fit BAO parameters smaller than 0.1%. Second, we obtained a more accurate measurement of the C IV bias parameter from [32], with no impact on the BAO results.', 'D Comparison of sampling methods': 'The main results in this publication (including figure 7) were obtained using the Nested Sampler Polychord [93, 94]. The quoted parameter values are given by the mean of the posterior distribution, and the reported uncertainties are the 68% credible regions. However, computing the full posterior distribution is computationally intensive, and in most of the tests in Section 6 we instead use a faster approximate method. This involves the use of the  \niminuit package [95] to find the maximum likelihood (minimum χ 2 ) point in parameter space. Approximate Gaussian uncertainties are then computed by taking the second derivative with respect to parameters around the best-fit point [96]. In figure 16 we show that both methods lead to very similar BAO contours. Therefore, the faster approximate method is good enough to check for shifts in the BAO position as done in Section 6.  \nIn contrast, previous Ly α forest BAO analyses used a frequentist approach to obtain their main results [see e.g. 21, 22, 27]. This involved using the Profile Likelihood method to create a two-dimensional χ 2 grid of α || and α ⊥ , and then calibrating the size of contours based on ∆ χ 2 values obtained from large sets of Monte Carlo simulations of the correlation functions. [74] found that BAO measurements obtained with this method agree well with Bayesian results based on the full posterior distribution. Therefore, in this work, we rely on the simpler and often faster method of sampling the full posterior for our main results.', 'E Significance of BAO shifts': "In section 6.4 we present an exhaustive list of robustness tests, where we look at the impact on the BAO parameters when changing different settings in the analysis. In most of these variations, the dataset is exactly the same, and therefore we expect shifts in the BAO parameters to be caused exclusively by the analysis settings. Before unblinding the measurements we required that none of these variations caused a shift in the BAO parameters larger than a third of the statistical uncertainty from the results on (blinded) data.  \nIn the second set of alternative analyses (in red) in figure 11, however, we also showed variations that caused minor changes in the dataset, by adding / removing quasars from the sample or by adding / removing Ly α pixels from the data vector. In these variations we relaxed the requirement, since changes in the dataset will cause statistical fluctuations in the measurement of the BAO parameters. This explains why the shifts shown in red in figure 11 are somewhat larger than those in the other variations.  \nTwo of the larger shifts correspond to the variations 'only quasar targets' and 'weak BALs' (see section 6.4 for details on the variations). These variations cause the datasets to be 11% and 8% less constraining than the main analysis (based on the increase in the errorbars on the measured correlations). Following [127], we estimate the statistical fluctuations for these variations to be on average 0 . 28 and 0 . 33 σ respectively. We conclude that the shifts on BAO parameters in these variations are therefore consistent with statistical fluctuations.  \nThere is another variation that has caused a shift larger than the requirement, and this is the 'mask-Lya redshift estimates' (∆ α ∥ = 0 . 001, ∆ α ⊥ = -0 . 013). In this variation, we have used an alternative redshift estimator that masks wavelengths bluer than the Ly α emission line of the quasar. The rms of the differences in redshifts above z > 2 is 443 km/s 33 , causing differences in pixel-quasar pairs of order 4 . 2 h -1 Mpc (using the fiducial cosmology to compute H ( z eff = 2 . 33)) and causing statistical fluctuations in the measurement of the crosscorrelation (see appendix B of dMdB20 ). Besides the expected effect in the cross-correlation, differences in the redshift estimates have a more subtle effect in the auto-correlation as well. The rest-frame wavelength of a given pixel changes with the quasar redshift, and we find that on average 2% of the pixels of a given Ly α forest are moved in or out of the restframe wavelength range used in the analysis, causing further statistical fluctuations.  \nFinally, we want to assess whether the shift caused by the use of a different redshift estimates is a systematic shift or if it can be explained by statistical fluctuations. For this  \nFigure 17 . Correlation matrix corresponding to the cross-covariance of SDSS DR16 and DESI DR1 correlations. The first block of 2 500 × 2 500 in the top left corresponds to the correlation matrix of the Ly α (A) × Ly α (A) measurement of SDSS, while the second block of the same size corresponds to the same measurement in DESI. The third, larger block of 5 000 × 5 000 corresponds to the Ly α (A) × QSO cross-correlation in SDSS, and the last block (bottom right) has the same measurement in DESI.  \n<!-- image -->  \npurpose, we perform a bootstrap analysis based on a random sampling with replacement of the HEALpix pixels used to compute the average correlation functions and we generate 1000 bootstrap samples for each of the two analyses. We then fit the BAO parameters for each of these samples and look at the distribution of shifts between the two analyses. We estimate with this technique that the statistical uncertainties on the shifts of α ∥ and α ⊥ are of 0.007 and 0.009, respectively. We conclude from this study that the observed shifts of ∆ α ∥ = 0 . 001 and ∆ α ⊥ = -0 . 013 are consistent with statistical fluctuations.", 'F Estimation of the DESI - SDSS covariance': 'The cross-covariance of the four correlation functions measured with the DESI DR1 Ly α forest dataset is shown in figure 6. As described in section 3.3, we compute a noisy estimate of this 15 000 × 15 000 covariance from the scatter of the measurements obtained in different HEALPix pixels. Following [27, 31], we smooth the correlation matrix so that it is invertible and we can use it to define a likelihood function.  \nHere we use the same method to compute the cross-covariance between the DESI DR1 correlations functions and those measured in dMdB20 using the SDSS DR16 Ly α dataset. A matrix describing the covariance of all eight correlations would be 30 000 × 30 000. However, given that the measurements using the B region carry a small amount of information (see the bottom right panel in figure 10) we ignore them in this appendix. In figure 17 we show the correlation matrix of SDSS DR16 and DESI DR1 measurements of Ly α (A) × Ly α (A) and Ly α (A) × QSO.  \nUsing this cross-covariance and the best-fit model from the combined analysis from table 5, we generated 1 000 Monte Carlo realisations of the two main correlation functions (Ly α (A) × Ly α (A) and Ly α (A) × QSO) of SDSS and of DESI, along with the correct crosscovariance between the surveys. We then minimised the likelihood for the SDSS and DESI correlation functions separately, and looked at the correlation between the best-fit BAO parameters ( α ⊥ , α ∥ ) to obtain the correlation matrix of the four BAO parameters:  \n     1 ρ ( α S ⊥ , α S ∥ ) ρ ( α S ⊥ , α D ⊥ ) ρ ( α S ⊥ , α D ∥ ) ρ ( α S ∥ , α S ⊥ ) 1 ρ ( α S ∥ , α D ⊥ ) ρ ( α S ∥ , α D ∥ ) ρ ( α D ⊥ , α S ⊥ ) ρ ( α D ⊥ , α S ∥ ) 1 ρ ( α D ⊥ , α D ∥ ) ρ ( α D ∥ , α S ⊥ ) ρ ( α D ∥ , α S ∥ ) ρ ( α D ∥ , α D ⊥ ) 1      =     1 -0 . 53 0 . 09 -0 . 05 -0 . 53 1 -0 . 05 0 . 10 0 . 09 -0 . 05 1 -0 . 50 -0 . 05 0 . 10 -0 . 50 1     , (F.1)  \nwhere the superscript S ( D ) refers to SDSS (DESI) measurements of BAO.  \nThe correlation of the BAO parameters from eBOSS and DESI can be computed separately from their posteriors, as discussed in the main text. For this reason, we modify the matrix above and instead use ρ ( α S ⊥ , α S ∥ ) = -0 . 45 and ρ ( α D ⊥ , α D ∥ ) = -0 . 48. We use this modified correlation matrix, and the variances of each measurement reported by each survey, to build a 4 × 4 multi-survey covariance C for the multi-survey data vector d = ( α S ⊥ , α S ∥ , α D ⊥ , α D ∥ ).  \nWe use these results to compute a combined BAO measurement d DS = ( α DS ⊥ , α DS ∥ ) with covariance C DS using linear algebra:  \nC -1 DS = S T C -1 S and C -1 DS d DS = S T C -1 d , (F.2)  \nwhere we have defined the matrix S as:  \nS =     1 0 0 1 1 0 0 1     . (F.3)  \nUsing these equations we obtain the following combined BAO measurements:  \nα DESI+SDSS ⊥ = 0 . 990 ± 0 . 019 , α DESI+SDSS ∥ = 1 . 012 ± 0 . 016 , (F.4)  \nwith correlation coefficient of ρ = -0 . 47.', 'G Author Affiliations': "1 Instituto de F'ısica Te'orica (IFT) UAM/CSIC, Universidad Aut'onoma de Madrid, Cantoblanco, E-28049, Madrid, Spain  \n- 59 Consejo Nacional de Ciencia y Tecnolog'ıa, Av. Insurgentes Sur 1582. Colonia Cr'edito Constructor, Del. Benito Ju'arez C.P. 03940, M'exico D.F. M'exico\n- 60 Centro de Investigaci'on Avanzada en F'ısica Fundamental (CIAFF), Facultad de Ciencias, Universidad Aut'onoma de Madrid, ES-28049 Madrid, Spain\n- 61 Excellence Cluster ORIGINS, Boltzmannstrasse 2, D-85748 Garching, Germany\n- 62 University Observatory, Faculty of Physics, Ludwig-Maximilians-Universitat, Scheinerstr. 1, 81677 Munchen, Germany\n- 63 Department of Astrophysical Sciences, Princeton University, Princeton NJ 08544, USA 64 Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK\n- 65 Department of Astronomy, The Ohio State University, 4055 McPherson Laboratory, 140 W18th Avenue, Columbus, OH 43210, USA\n- 66 Department of Physics, Southern Methodist University, 3215 Daniel Avenue, Dallas, TX 75275, USA\n- 67 Department of Physics and Astronomy, University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada\n- 68 Perimeter Institute for Theoretical Physics, 31 Caroline St. North, Waterloo, ON N2L 2Y5, Canada\n- 69 Waterloo Centre for Astrophysics, University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada\n- 70 Graduate Institute of Astrophysics and Department of Physics, National Taiwan University, No. 1, Sec. 4, Roosevelt Rd., Taipei 10617, Taiwan\n- 71 Sorbonne Universit'e, CNRS/IN2P3, Laboratoire de Physique Nucl'eaire et de Hautes Energies (LPNHE), FR-75005 Paris, France\n- 72 Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95065, USA\n- 73 Department of Astronomy & Astrophysics, University of Toronto, Toronto, ON M5S 3H4, Canada  \n- 84 Departament de F'ısica, EEBE, Universitat Polit'ecnica de Catalunya, c/Eduard Maristany 10, 08930 Barcelona, Spain\n- 85 Universit'e Clermont-Auvergne, CNRS, LPCA, 63000 Clermont-Ferrand, France\n- 86 University of California Observatories, 1156 High Street, Sana Cruz, CA 95065, USA\n- 87 Department of Physics & Astronomy, Ohio University, Athens, OH 45701, USA\n- 88\n- Department of Physics and Astronomy, Sejong University, Seoul, 143-747, Korea 89\n- Abastumani Astrophysical Observatory, Tbilisi, GE-0179, Georgia\n- 90 Faculty of Natural Sciences and Medicine, Ilia State University, 0194 Tbilisi, Georgia\n- 91 Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA\n- 92 Centre for Advanced Instrumentation, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK\n- 93 Physics Department, Brookhaven National Laboratory, Upton, NY 11973, USA 94 Beihang University, Beijing 100191, China  \n- 97 Department of Physics, University of California, Berkeley, 366 LeConte Hall MC 7300, Berkeley, CA 94720-7300, USA"}