# 3D Correlations in the Lyman-α Forest from Early DESI Data

Calum Gordon, Andrei Cuceu, Jonás Chaves-Montero, A Font-Ribera, A.X González-Morales, J Aguilar, S Ahlen, E Armengaud, S Bailey, A Bault, et al.

► To cite this version:

Calum Gordon, Andrei Cuceu, Jonás Chaves-Montero, A Font-Ribera, A.X González-Morales, et al. 3D Correlations in the Lyman-α Forest from Early DESI Data. JCAP, 2023, 11, pp.045. 10.1088/1475-7516/2023/11/045. [hal-04190251]

HAL Id: hal-04190251

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

Submitted on 31 Jul 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
PAPER

# 3D correlations in the Lyman-α forest from early DESI data

To cite this article: C. Gordon et al. JCAP11(2023)045

View the [article online](#) for updates and enhancements.

## You may also like

- Performance of the CMS high-level trigger during LHC Run 2

A. Hayrapetyan, A. Tumasyan, W. Adam et al.

- Identification of low-momentum muons in the CMS detector using multivariate techniques in proton-proton collisions at sqrt(s) = 13.6 TeV

V. Chekhovsky, A. Hayrapetyan, V. Makarenko et al.

- DESI 2024 III: baryon acoustic oscillations from galaxies and quasars

A.G. Adame, J. Aguilar, S. Ahlen et al.
<figure><img src="image_2.png" /></figure>

JCAP ANNIVERSARY
SPECIAL ISSUE

3D correlations in the Lyman-α forest from early
DESI data

C. Gordon,<sup>a</sup> A. Cuceu,<sup>b,1</sup> J. Chaves-Montero,<sup>a</sup> A. Font-Ribera,<sup>a,c</sup>

A.X. González-Morales,<sup>d,e</sup> J. Aguilar,<sup>f</sup> S. Ahlen,<sup>g</sup> E. Armengaud,<sup>h</sup> S. Bailey,<sup>f</sup> A. Bault,<sup>i</sup>

A. Brodzeller,<sup>j</sup> D. Brooks,<sup>c</sup> T. Claybaugh,<sup>f</sup> R. de la Cruz,<sup>e</sup> K. Dawson,<sup>j</sup> P. Doel,<sup>c</sup>

J.E. Forero-Romero,<sup>k,l</sup> S. Gontcho A Gontcho,<sup>f</sup> J. Guy,<sup>f</sup> H.K. Herrera-Alcantar,<sup>e</sup>

V. Iršič,<sup>m</sup> N. G. Karaçaylı,<sup>b</sup> D. Kirkby,<sup>i</sup> M. Landriau,<sup>f</sup> L. Le Guillou,<sup>n</sup> M.E. Levi,<sup>f</sup>

A. de la Macorra,<sup>o</sup> M. Manera,<sup>a,p</sup> P. Martini,<sup>b</sup> A. Meisner,<sup>q</sup> R. Miquel,<sup>a,r</sup>

P. Montero-Camacho,<sup>s</sup> A. Muñoz-Gutiérrez,<sup>o</sup> L. Napolitano,<sup>t</sup> J. Nie,<sup>u</sup> G. Niz,<sup>e,v</sup>

N. Palanque-Delabrouille,<sup>f,h</sup> W.J. Percival,<sup>w,x</sup> M. Pieri,<sup>y</sup> C. Poppett,<sup>f,z</sup> F. Prada,<sup>aa</sup>

I. Pérez-Ràfols,<sup>ab</sup> C. Ramírez-Pérez,<sup>a</sup> C. Ravoux,<sup>h,ac</sup> M. Rezaie,<sup>ad</sup> A.J. Ross,<sup>b</sup> G. Rossi,<sup>ae</sup>

E. Sanchez,<sup>af</sup> D. Schlegel,<sup>f</sup> M. Schubnell,<sup>ag</sup> H. Seo,<sup>ah</sup> F. Sinigaglia,<sup>ai</sup> T. Tan,<sup>n</sup>

G. Tarlé,<sup>aj</sup> M. Walther,<sup>ak</sup> B.A. Weaver,<sup>p</sup> C. Yèche,<sup>h</sup> Z. Zhou<sup>u</sup> and H. Zou<sup>u</sup>

<sup>a</sup>Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, 08193 Bellaterra (Barcelona), Spain

<sup>b</sup>Center for Cosmology and Astro-Particle Physics, The Ohio State University, Columbus, Ohio 43210, U.S.A.

<sup>c</sup>Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, U.K.

<sup>d</sup>Consejo Nacional de Ciencia y Tecnología, Av. Insurgentes Sur 1582. Colonia Crédito Constructor, Del. Benito Juárez C.P. 03940, México D.F., Mexico

<sup>e</sup>Departamento de Física, Universidad de Guanajuato — DCI, C.P. 37150, Leon, Guanajuato, Mexico

<sup>f</sup>Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, U.S.A.

<sup>g</sup>Physics Dept., Boston University, 590 Commonwealth Avenue, Boston, MA 02215, U.S.A.

<sup>h</sup>IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France

<sup>i</sup>Department of Physics and Astronomy, University of California, Irvine, 92697, U.S.A.

<sup>j</sup>Department of Physics and Astronomy, The University of Utah, 115 South 1400 East, Salt Lake City, UT 84112, U.S.A.

<footnote n="1">¹NASA Einstein Fellow.</footnote>
<sup>k</sup>Observatorio Astronómico, Universidad de los Andes,
Cra. 1 No. 18A-10, Edificio H, CP 111711 Bogotá, Colombia

<sup>l</sup>Departamento de Física, Universidad de los Andes,
Cra. 1 No. 18A-10, Edificio Ip, CP 111711, Bogotá, Colombia

<sup>m</sup>Kavli Institute for Cosmology, University of Cambridge,
Madingley Road, Cambridge CB3 0HA, U.K.

<sup>n</sup>Sorbonne Université, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies
(LPNHE), FR-75005 Paris, France

<sup>o</sup>Instituto de Física, Universidad Nacional Autónoma de México, Cd. de México C.P. 04510, México

<sup>p</sup>Departament de Física, Serra Húnter, Universitat Autònoma de Barcelona,
08193 Bellaterra (Barcelona), Spain

<sup>q</sup>NSF's NOIRLab, 950 N. Cherry Ave., Tucson, AZ 85719, U.S.A.

<sup>r</sup>Institució Catalana de Recerca i Estudis Avançats,
Passeig de Lluís Companys, 23, 08010 Barcelona, Spain

<sup>s</sup>Department of Astronomy, Tsinghua University,
30 Shuangqing Road, Haidian District, Beijing, China, 100190

<sup>t</sup>Department of Physics & Astronomy, University of Wyoming,
1000 E. University, Dept. 3905, Laramie, WY 82071, U.S.A.

<sup>u</sup>National Astronomical Observatories, Chinese Academy of Sciences,
A20 Datun Rd., Chaoyang District, Beijing, 100012, P.R. China

<sup>v</sup>Instituto Avanzado de Cosmología A. C.,
San Marcos 11 — Atenas 202. Magdalena Contreras, 10720. Ciudad de México, México

<sup>w</sup>Perimeter Institute for Theoretical Physics,
31 Caroline St. North, Waterloo, ON N2L 2Y5, Canada

<sup>x</sup>Waterloo Centre for Astrophysics, University of Waterloo,
200 University Ave W, Waterloo, ON N2L 3G1, Canada

<sup>y</sup>Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France

<sup>z</sup>Space Sciences Laboratory, University of California,
Berkeley, 7 Gauss Way, Berkeley, CA 94720, U.S.A.

<sup>aa</sup>Instituto de Astrofísica de Andalucía (CSIC),
Glorieta de la Astronomía, s/n, E-18008 Granada, Spain

<sup>ab</sup>Departament de Física Quàntica i Astrofísica, Universitat de Barcelona,
Martí i Franquès 1, E08028 Barcelona, Spain

<sup>ac</sup>Aix Marseille Univ, CNRS/IN2P3, CPPM, Marseille, France

<sup>ad</sup>Department of Physics, Kansas State University, 116 Cardwell Hall, Manhattan, KS 66506, U.S.A.

<sup>ae</sup>Department of Physics and Astronomy, Sejong University, Seoul, 143-747, Korea

<sup>af</sup>CIEMAT, Avenida Complutense 40, E-28040 Madrid, Spain

<sup>ag</sup>Department of Physics, University of Michigan, Ann Arbor, MI 48109, U.S.A.

<sup>ah</sup>Department of Physics & Astronomy, Ohio University, Athens, OH 45701, U.S.A.

<sup>ai</sup>Instituto de Astrofísica de Canarias, C/ Vía Láctea, s/n, E-38205 La Laguna, Tenerife, Spain

<sup>aj</sup>University of Michigan, Ann Arbor, MI 48109, U.S.A.

<sup>ak</sup>University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität,
Scheinerstr. 1, 81677 München, Germany

E-mail: cgordon@ifae.es
**ABSTRACT:** We present the first measurements of Lyman-$\alpha$ (Ly$\alpha$) forest correlations using early data from the Dark Energy Spectroscopic Instrument (DESI). We measure the auto-correlation of Ly$\alpha$ absorption using 88 509 quasars at $z > 2$, and its cross-correlation with quasars using a further 147 899 tracer quasars at $z \gtrsim 1.77$. Then, we fit these correlations using a 13-parameter model based on linear perturbation theory and find that it provides a good description of the data across a broad range of scales. We detect the BAO peak with a signal-to-noise ratio of 3.8$\sigma$, and show that our measurements of the auto- and cross-correlations are fully-consistent with previous measurements by the Extended Baryon Oscillation Spectroscopic Survey (eBOSS). Even though we only use here a small fraction of the final DESI dataset, our uncertainties are only a factor of 1.7 larger than those from the final eBOSS measurement. We validate the existing analysis methods of Ly$\alpha$ correlations in preparation for making a robust measurement of the BAO scale with the first year of DESI data.

**KEYWORDS:** Lyman alpha forest, baryon acoustic oscillations, dark energy experiments

ARXIV ePRINT: 2308.10950
# Contents

<table><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Data</td><td>3</td></tr><tr><td>2.1 Quasar spectra</td><td>3</td></tr><tr><td>2.2 Pixel masking</td><td>5</td></tr><tr><td>2.3 The flux transmission field</td><td>6</td></tr><tr><td>3 Measuring the correlation functions</td><td>8</td></tr><tr><td>3.1 Lyα weights</td><td>8</td></tr><tr><td>3.2 Continuum distortion</td><td>8</td></tr><tr><td>3.3 Auto-correlation</td><td>9</td></tr><tr><td>3.4 Cross-correlation</td><td>10</td></tr><tr><td>3.5 Covariance matrix</td><td>11</td></tr><tr><td>4 Correlation function models</td><td>13</td></tr><tr><td>4.1 Power spectra</td><td>13</td></tr><tr><td>4.2 Non-linear corrections</td><td>14</td></tr><tr><td>4.3 Modeling contaminants</td><td>14</td></tr><tr><td>4.3.1 Quasar radiation effects</td><td>15</td></tr><tr><td>4.3.2 Instrumental systematics</td><td>15</td></tr><tr><td>4.3.3 Absorption by high column density systems</td><td>16</td></tr><tr><td>4.3.4 Absorption by metal lines</td><td>17</td></tr><tr><td>4.4 Modelling the continuum distortion</td><td>18</td></tr><tr><td>5 Fits to the data</td><td>18</td></tr><tr><td>5.1 Fit probability</td><td>20</td></tr><tr><td>5.2 Comparison with eBOSS DR16</td><td>22</td></tr><tr><td>6 Conclusion</td><td>24</td></tr><tr><td>A Complete fit results</td><td>25</td></tr></table>

## 1 Introduction

The Lyman-α (Lyα) forest is a dense series of absorption lines found in the spectra of distant quasars, that we use to study the large-scale structure of the universe. These absorption features arise as a result of diffuse neutral hydrogen in the intergalactic galactic medium [1], which readily absorbs light at the Lyα wavelength $\lambda_{Lyα} = 1215.67$ Å, and can be used as a biased tracer of the underlying matter density [2]. Due to the low rest-frame wavelength of the Lyα transition, its absorption is only detected by optical spectrographs for quasars
with redshift $z \ge 2$. The Ly$\alpha$ forest allows us to measure the expansion of the Universe at redshifts larger than those accessible by spectroscopic galaxy surveys ($z < 1.5$, [3]).

The environment that gives rise to the Ly$\alpha$ forest has been well studied using hydrodynamical simulations [4, 5], where a combination of photo-ionisation heating and adiabatic cooling leads to a tight relation between temperature and gas density $T \propto (\rho/\bar{\rho})^{\gamma-1}$ [6]. The transmitted flux fraction is related to the optical depth of the Ly$\alpha$ forest by $F = e^{-\tau}$. At sufficiently large scales ($> 10h^{-1}\text{Mpc}$), the power spectrum of transmitted flux $\delta_F = F/\bar{F} - 1$, is found to trace linearly the power spectrum of $\delta(x) = \rho(x)/\bar{\rho} - 1$ [7]. However, due to non-linear growth, thermal broadening, and Jeans smoothing, this relation becomes more complex and non-linear at smaller scales [8–10].

The first studies of clustering in the Ly$\alpha$ forest looked at one-dimensional (1D) correlations along the line-of-sight towards a handful of quasars [11, 12]. These 1D measurements were later repeated using larger data sets from the Sloan Digital Sky Survey (SDSS; [13]), from the Baryon Acoustic Oscillation Survey (BOSS; [14]) and from its extension (eBOSS; [15]).

This paper, on the other hand, focuses on three-dimensional (3D) correlations in the Ly$\alpha$ forest. Using early data from the BOSS survey, [7] presented the first measurement of 3D correlations in the Ly$\alpha$ forest. These correlations were also computed from larger datasets to measure the scale of Baryon Acoustic Oscillations (BAO) around $z = 2.3$, in BOSS and eBOSS [16–22]. The cross-correlation of quasars and the Ly$\alpha$ forest was also measured in BOSS [23], and it was used to improve the BAO measurements from the Ly$\alpha$ datasets [22, 24–26]. The final eBOSS analysis [22], using data from SDSS DR16, included 210 005 Ly$\alpha$ forest from quasars at $z > 2.1$ that were used in the auto-correlation, and 341 468 at $z > 1.77$ that were also used in the cross-correlation. More recently, the full shape of the 3D Ly$\alpha$ correlation has been used to measure the Alcock-Paczynski effect in eBOSS data [27], which combined with Ly$\alpha$ BAO gives the best measurement of the expansion rate of the Universe above $z = 1$.

We use early data from the Dark Energy Spectroscopic Instrument (DESI) [28, 29] to measure the auto-correlation of the Ly$\alpha$ forest and its cross-correlation with quasars, and compare our results with those from the eBOSS collaboration [22]. We also show in section 5 that we detect the BAO peak to high-confidence in our dataset. Given the early stage and relatively low statistical power of our data, we choose not to present measurements of the BAO scale parameters. Throughout the latter half of the paper, we make use of the best fit $\Lambda$CDM cosmological parameters from Planck 2018 [30],<fnref n="1" /> where $\Omega_M = 0.3153$, to compute correlations in terms of co-moving separation. This allows us to combine a wide range of redshifts into a single measurement while preserving the BAO feature.

The structure of the paper follows roughly the steps required to perform the end-to-end analysis. We start by describing the DESI instrument and the early data release, and outline how we measure the Ly$\alpha$ transmission catalogue from high-redshift quasar spectra in section 2. We then compute the auto- and cross-correlation and covariance matrices in section 3. In section 4, we discuss the physical model of Ly$\alpha$ forest correlations and contaminants, and in section 5 we present the results of fitting this model to our data and discuss the consistency of these results with the final eBOSS analysis [22]. Finally, we summarise our results in section 6.

<footnote n="1">¹See table 2 of [30], where we use TT,TE,EE+lowE+lensing results.</footnote>
The Lyα forest dataset used here is described in detail in [31]. Our analysis is also accompanied by a paper analysis on early DESI mock datasets [32] which utilises the same 13-parameter model to describe the auto- and cross-correlation. We further highlight two papers where systematic effects, including quasar redshift errors [33] and instrumental effects [34], are studied on the same early DESI data set.

## 2 Data

DESI is a multi-object spectrograph installed at the Mayall 4-meter telescope at Kitt Peak National Observatory [35]. DESI is equipped with 5000 robotically controlled fibres across 10 independent spectrographs, that have an approximately 3-degree field of view in each pointing [36–38]. Each measured spectrum covers an observed wavelength range between 3600 and 9800 Å with a pixel width of 0.8 Å, which is split between 3 independent spectrograph arms ranging between 3600–5930, 5600–7720, and 7470–9800 Å respectively [35].

In May 2021, DESI started a 5-year campaign which will observe 3 million quasar [39] and 36 million galaxy targets [40–42], over an area of 14 000 deg<fnref n="2" />. These targets are selected from the DESI Legacy Imaging Survey [43–45]. These observations are made thanks to extensive supporting software packages: one to select targets for spectroscopic observation [46], one to assign fibres to targets [47], one to optimise survey operations [48], an exposure-time calculator [49], and a spectroscopic pipeline to reduce the data and calibrate the spectra [50].

We analyse data from the Early Data Release (EDR) of DESI [51], which is publicly available<fnref n="2" /> and contains 68 750 quasars. EDR contains data from a pre-survey designed to optimise to optimise target selection algorithms [39, 52], refine the observational procedure, and generally evaluate the quality of spectra and the redshift accuracy [53, 54]. To increase the statistical power of the measurements, we also include data from the first two months of the main survey (M2) containing 249 941 quasars. We use “EDR+M2” to refer to the combination of both data sets.

The sample footprint is shown in figure 1, where we highlight the each survey phase and the eBOSS Data Release 16 (DR16) survey footprint [22].

Even though we analyse all sources from EDR+M2 together, the typical signal-to-noise (SNR) of M2 and EDR spectra varies. Targets observed in M2 are limited to one observing pass,<fnref n="3" /> which will increase to 4 by the end of DESI observations. Quasars identified as Lyα quasars ($z > 2$) receive more passes in general than other targets to reach the required SNR per pixel across the Lyα forest. On the other hand, roughly half of EDR spectra present 4 passes and half 1.3 [39].

### 2.1 Quasar spectra

The DESI pipeline uses a robust spectroscopic reduction pipeline [50], in combination with a template-fitting algorithm Redrock (RR) [55] to give object classifications and estimate their redshifts. Templates used by RR are formed by a linear combination of bases resulting

<footnote n="2">²https://data.desi.lbl.gov/public/edr/.</footnote>

<footnote n="3">³Equivalent to nominal 1000s of effective exposure time, which is exposure time accounting for observing conditions and fibre effects.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Footprint of quasar targets in DESI Early Data Release (purple), the first two months of main survey (red), and eBOSS DR16 (blue) [22]. The EDR and M2 surveys cover an area of 250.1 deg² and 1290.9 deg².</figcaption></figure>

from the principal component analysis (PCA) decomposition of Sloan Digital Sky Survey spectra [55] (see [54] for comparison when using a new set of quasar templates).

To quantify the ability of the pipeline to classify quasars we use two terms, purity and completeness. Purity is the percentage of classified quasars that are true quasars, which is > 99% in our pipeline. Completeness is the percentage of true quasars identified as quasars, which we find to be ~86%. To increase this, we run two algorithms on the RR output. The first attempts to identify the MgII broad emission line — a feature only found in quasars — in the spectra of sources classified as galaxies [39], and is especially effective at 0.5 < z < 1.5. The other is QuasarNet (QN) [56, 57], a convolutional neural network trained using a large number of visually inspected BOSS spectra, to identify faint quasars missed by both RR and the first afterburner, and is more useful at z > 2. These algorithms increase the spectral completeness from ~86% to ~94% [53], and only cause a ~1–2% drop in spectral purity at z < 1 [39].

Redrock is very reliable in estimating the redshifts of high-quality spectra [53]. The number of catastrophic redshifts, where the difference between visually-inspected and RR redshifts is $\Delta z > 0.01$, is only ~1.5%. For sources presenting a disagreement between RR and QN redshifts, the pipeline computes the final redshift by re-running RR using the QN redshift as prior and only quasar templates. This reduces the fraction of catastrophic redshifts from 1.5 to 1% [53].

We obtain an overall density of 210 quasars deg⁻² (60 for z > 2.1) [39], which surpasses the DESI requirements of 170 quasars deg⁻² (50 at z > 2.1) [29]. Throughout this work, we will only analyse spectroscopically identified quasars that were also selected as quasar targets (i.e. we ignore the small fraction of quasars found in other target classes) [46], to keep the purity of the catalogue as high as possible. These have a density of around 54 deg⁻² at z > 2.1.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Redshift distribution of high-redshift quasars in DESI EDR+M2 and eBOSS DR16 [22], with 147 899 and 341 468 in each sample respectively. Only quasars between $z = 1.77$ and $z = 3.75$ feature in our analysis.</figcaption></figure>

In figure 2 we show the redshift distributions of our EDR+M2 $Ly\alpha$ quasars and include for comparison the catalogue used in the $Ly\alpha$ analysis of eBOSS DR16 ([22], hereby dMdB20).

Figure 3 shows one such quasar with $z = 2.44$. Purple labels indicate broad, distinctive emission lines often observed in high-redshift quasars, including the MgII line. We also highlight the $Ly\alpha$ forest in blue, which exhibits a high density of absorption lines. Due to the rest-frame wavelength range of the $Ly\alpha$ forest and the observed redshift distribution of quasars, most of our data falls within the blue arm of the spectrographs.

## 2.2 Pixel masking

The $Ly\alpha$ catalogue used to obtain the results in this paper is described in [31]. We present here a summary of how we mask contaminants in the data, and refer the reader to [31] for a detailed description of the process.

In this section and while computing the flux transmission field, we use $Ly\alpha$ pixels of width 0.8 Å in the observed frame — the resolution provided by the spectroscopic pipeline [50]. When computing the correlation functions (section 3) and the metal matrices (section 4.3.4), we combine these into 2.4 Å analysis pixels to reduce the computation time, without affecting the measured correlations.

The two primary contaminants of our $Ly\alpha$ forests are broad absorption line features (BALs) and damped $Ly\alpha$ absorbers (DLAs). BALs are broad spectral troughs that arise from high-velocity gas outflows from quasars. They often fall into the $Ly\alpha$ forest region and are difficult to model. These have been identified in the past using convolutional neural networks [58], and a similar approach has been used in DESI [59]. Identified BALs are then masked in the affected spectra depending on the width of the trough, following the procedure described in [60]. Masking avoids discarding whole $Ly\alpha$ forests where there are often valuable pixels, as was done in previous $Ly\alpha$ forest BAO studies, like dMdB20.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** An example of a DESI EDR spectra at z = 2.44. Black dashed lines indicate the central wavelength of the strongest quasar emission lines, including the Lyα used in our analysis. The shaded blue region shows the Lyα forest between the Lyα and Lyβ lines, the latter only partially observed in this specific example. Also of note is the MgII line, used in the pipeline to re-classify galaxies as quasars and the region between the CIV and CIII lines, used for the re-calibration step described in section 2.3.</figcaption></figure>

DLAs are large absorption systems with column densities of neutral hydrogen $N_{HI} > 10^{20.3}cm^{-2}$. Like the Lyα forest, DLAs are a tracer of the underlying matter distribution but present a higher bias than the Lyα forest [61–63]. Their wide absorption profiles, however, increase the noise in the clustering measurements. Following dMdB20, we mask the regions where more than 20% of flux is absorbed DLAs and correct their damping wings in the unmasked regions using a Voigt profile. We use a DLA catalog created using two different techniques, a convolutional neural network (CNN) and a Gaussian process (GP) finder, used in previous Lyα analyses [64, 65] and re-developed for DESI [66, 67]. We decide to mask only those DLAs that are detected by both algorithms with a confidence higher than 50%, resulting in 7% (or ~8500) of the forests in our analysis containing at least one DLA.

To further remove noise from our data, we perform two additional steps. The first is to mask areas affected by sharp spectral features, like the bright [OI] 5577.3 Å sky line and the calcium H and calcium K lines from the interstellar medium. We mask the wavelength ranges between [5570.5,5586.7] Å, [3967.3,3971.0] Å and [3933.0,3935.8] Å respectively to avoid these [31]. Other features caused by systematic flux calibration errors are corrected using the CIII] spectral region (between CIV and CIII] in figure 3), see [31] for more details.

## 2.3 The flux transmission field

Here we define the limits on our Lyα forests, then describe the measurement of the flux transmission field used in our correlation functions. We only discuss the most important steps; for a complete description of these processes and the DESI EDR Lyα catalogue, see [31].

In the observed wavelength regime we restrict our forests between $\lambda_{min} = 3600$ Å and $\lambda_{max} = 5772$ Å. The lower limit is motivated by the limit of the DESI spectrograph, which for
our Lyα catalogue translates to only including quasars above z~2. We set the upper limit to maximise the signal-to-noise of our correlation functions. In the rest-frame, our wavelength range is also reduced from the full extent of the Lyα forest to between $\lambda_{RF,min} = 1040$ Å and $\lambda_{RF,max} = 1205$ Å, to avoid the wings of the Lyα and Lyβ emission lines, which can vary from quasar to quasar and thus add error to the continuum fitting. Compared to dMdB20, we have a higher upper wavelength limit in both the observed and rest-frame cases. In [31], they found that $\lambda_{RF,max} = 1205$ Å was the limit corresponding to the lowest possible variance in the Lyα correlation function monopole, achieving a balance between the number of Lyα forest pixels and avoiding the wings of the Lyα emission line.

The limits above result in 109 900 Lyα forests, with a maximum Lyα absorption pixel redshift of $z = 3.75$. In practice we also place an upper limit on the quasar redshift of $z = 3.75$ to be consistent with our mock catalogs [32], but this is done when computing the correlation functions and removes only a negligible number of forests. For the Lyα-quasar cross-correlation, we use a tracer catalog of quasars with a minimum redshift of $z > 1.77$, containing 147 899 quasars. This limit corresponds to a maximum possible separation of $300h^{-1}$Mpc in the cross-correlation. In practice, we only measure correlations up to $200h^{-1}$Mpc but scales up to $300h^{-1}$Mpc are needed to accurately estimate the distortion (section 4.4). We place final restrictions on our sample by rejecting forests with fewer than 50 valid Lyα forest pixels, and forests with a failed continuum fit [31]. The resulting dataset contains 88 509 Lyα forests.

Fluctuations around the mean transmitted flux fraction ($\delta_q(\lambda)$) are computed as a function of wavelength as:

$$ \delta_q(\lambda) = \frac{f_q(\lambda)}{\bar{F}(z)C_q(\lambda)} - 1, \qquad (2.1) $$

where $f_q(\lambda)$ is the observed flux of quasar q, $\bar{F}(z)$ is the mean transmission of the IGM at redshift $z$ ($z(\lambda) = \lambda/\lambda_{Ly\alpha} - 1$) and $C_q(\lambda)$ is the unabsorbed quasar continuum. The product of the two terms in the denominator is equivalent to the mean expected flux for quasar q at observed wavelength $\lambda$.

To measure the flux transmission field of each forest, we first estimate the product of the denominator terms. The first step is to make the approximation:

$$ \bar{F}(z)C_q(\lambda) = \bar{C}(\lambda_{RF}) \left( a_q + b_q \frac{\Lambda - \Lambda_{min}}{\Lambda_{max} - \Lambda_{min}} \right) \quad \Lambda = \log \lambda, \qquad (2.2) $$

where $\bar{C}(\lambda_{RF})$ is an estimate of the mean rest-frame spectra of all quasars, with the two quasar parameters ($a_q, b_q$) adjusting for quasar spectral diversity and the redshift evolution of the mean flux $\bar{F}(z)$.

Using a first estimate of the continuum of each quasar, we characterise the flux variance ($\sigma_q^2$) in the pixels, later used to calculate the pixel weights (equation (3.1)), as the sum of two terms:

$$ \sigma_q^2 = [\bar{F}C_q(\lambda)]^2 \sigma_{LSS}^2(\lambda) + \eta(\lambda) \sigma_{pip,q}^2. \qquad (2.3) $$

Both $\eta(\lambda)$ and $\sigma_{LSS}^2(\lambda)$ functions are fitted iteratively with the mean quasar continuum $\bar{C}(\lambda_{RF})$ and the per-quasar parameters ($a_q, b_q$). We use the publicly available code Picca<fnref n="4" /> to compute

<footnote n="4">⁴https://github.com/igmhub/picca/.</footnote>
these fits, and they typically converge in five iterations. The first term captures the intrinsic variance of the Lyα forest, and the second term the instrumental noise variance. The noise variance estimated by the pipeline ($\sigma_{pip,q}^2$) is normalized with the expected flux and is further multiplied by a function $\eta(\lambda)$ to correct for possible mis-calibrations of the instrumental noise. In the dMdB20 analysis, $\eta(\lambda)$ ranged from 1.05 and 1.2 between 3600 and 5500 Å respectively. Over the same range in observed wavelength of the EDR+M2 dataset, this correction ranges from 1 to 1.01 [31], demonstrating that the DESI pipeline gives a better noise estimate.

# 3 Measuring the correlation functions

In this section we present the measurement of the correlation functions. We first describe the weights entering in the estimator of the correlation function in section 3.1, and discuss a small correction applied to the forest to model the distortion introduced by the continuum fitting (section 3.2). We then proceed to measure the auto-correlation of the Lyα forest (section 3.3) and its cross-correlation with quasars (section 3.4), as well as the uncertainty in these measurements (section 3.5).

## 3.1 Lyα weights

Because of varying quasar brightness and exposure time, the signal-to-noise in our Lyα forest is quite diverse. Therefore we weight the Lyα pixels when measuring correlations. These weights are used to compute the auto-correlation in equation (3.7) and the cross-correlation in equation (3.8). An optimal analysis would weight all **pairs** of pixels according to their inverse covariance matrix [17, 68]. However, this makes the analysis significantly more complex, so we follow previous analyses of eBOSS data (dMdB20) and weight pixels individually.

Instead of directly using the inverse of the pixel variance (equation (2.3)), we use the modified Lyα weights introduced in [31]:

$$w_i = \frac{[(1 + z_i)/(1 + z_{\text{fid}})]^{\gamma_{\text{Ly}\alpha}-1}}{\sigma_{\text{LSS}}^2(\lambda)\sigma_{\text{mod}}^2 + \eta(\lambda)\sigma_{\text{pip,q}}^2 / [\overline{FC_q}(\lambda)]^2}, \quad (3.1)$$

where $z_{\text{fid}} = 2.25$ (although this factor cancels in equation (3.7)). The extra parameter $\sigma_{\text{mod}}^2$ was introduced in [31] to modulate the relative importance of instrumental noise and intrinsic Lyα fluctuations in the weights. The value of $\sigma_{\text{mod}}^2 = 7.5$ was found to be optimal for this dataset [31], and it reduces the uncertainties on the auto- and cross-correlation by 20% and 10% respectively.

Finally, the numerator of equation (3.1) takes into account the redshift evolution of the Lyα bias ($\gamma_{\text{Lyα}}$) and of the growth factor, and up-weights higher redshift pixels to favour absorption with higher signal. Following dMdB20, we use a value of $\gamma_{\text{Lyα}} = 2.9$, motivated by the observation that the amplitude of the Lyα power spectrum with redshift evolves as $(1 + z)^{3.8}$ [13]. This correction ranges from $\approx 0.85$ to 2.3 between $z = 2$ and 4.

## 3.2 Continuum distortion

As discussed in section 2.3, while fitting the quasar continua we fit two parameters per quasar ($a_q, b_q$) that modulate the amplitude and slope of the continuum in each forest. However, this
biases the mean and slope of each $\delta_q(\lambda)$ toward 0 within each line of sight and distorts the measured correlations, often referred to as the continuum distortion [7, 20]. Following [20], and to simplify the modelling of this distortion in section 4.4, we modify the fluctuations $\delta_q$ and make this distortion explicit:

$$ \hat{\delta}_i = \sum_i \eta_{ij} \delta_j, \qquad (3.2) $$

where the projection matrix $\eta_{ij}$ is given by

$$ \eta_{ij}^{q} = \delta_{ij}^{K} - \frac{w_{j}}{\sum_{k} w_{k}} - \frac{w_{j} \kappa_{i} \kappa_{j}}{\sum_{k} w_{k} \kappa_{k}^{2}} \quad \kappa_{k} = \log \lambda_{k} - \overline{\log \lambda_{q}}, \qquad (3.3) $$

where $\delta_{ij}^K$ is the Kronecker delta and weights $w$ are given in equation (3.1).

Finally, the transformation in equation (3.2) means that the average $\delta_q$ at wavelength $\lambda$ is biased towards 0. Following dMdB20, we subtract $\overline{\delta(\lambda)}$ when calculating the cross-correlation (section 3.4) such that:

$$ \tilde{\delta}_q(\lambda) = \hat{\delta}_q(\lambda) - \overline{\delta(\lambda)}. \qquad (3.4) $$

This correction guarantees that the cross-correlation (section 3.4) tends to 0 at large separations regardless of the quasar redshift distribution.

## 3.3 Auto-correlation

To measure the correlation functions, redshift and angular ($\Delta z, \Delta \theta$) separations are first transformed to co-moving separations along ($r_\|$) and perpendicular ($r_\perp$) to the line of sight. We calculate ($r_\|, r_\perp$) as:

$$ r_{\|} = [D_{\mathrm{C}}(z_{\mathrm{i}}) - D_{\mathrm{C}}(z_{\mathrm{j}})]\cos\left(\frac{\Delta\theta}{2}\right), \qquad (3.5) $$

and

$$ r_{\perp} = [D_{\mathrm{M}}(z_{\mathrm{i}}) + D_{\mathrm{M}}(z_{\mathrm{j}})]\sin\left(\frac{\Delta\theta}{2}\right), \qquad (3.6) $$

where (i,j) represent pixel-pixel or pixel-quasar pairs and $\Delta\theta$ the angle between them. The redshift of a pixel i is calculated assuming Ly$\alpha$ absorption, such that $z_i = \lambda_{\mathrm{OBS}}/\lambda_{\mathrm{Ly}\alpha} - 1$, where $\lambda_{\mathrm{Ly}\alpha} = 1215.67$ Å. The co-moving distance is given by $D_{\mathrm{C}} = c/100 \int_{0}^{z} \frac{dz'}{E(z')}$ in units of $h^{-1}\mathrm{Mpc}$, where $E(z) = \sqrt{\Omega_{\mathrm{M}}(1+z)^3 + \Omega_{\Lambda}}$. $D_{\mathrm{M}}$ is the transverse co-moving distance (also known as angular co-moving distance), and $D_{\mathrm{M}} = D_{\mathrm{C}}$ [69], as we assume $\Omega_{\mathrm{k}} = 0$.

For the auto-correlation measurement we use a weighted covariance estimator following previous analyses [7, 22]:

$$ \xi_{\mathrm{A}} = \frac{\sum_{i,j \in \mathrm{A}} w_i w_j \delta_i \delta_j}{\sum_{i,j \in \mathrm{A}} w_i w_j}, \qquad (3.7) $$

where A is a bin in ($r_\|, r_\perp$) space with width $4h^{-1}\mathrm{Mpc}$ and the weights $w_i$ are described in section 3.1. We sum over all pixels (i,j) across different lines of sight, but do not include pixels in the same line of sight because of continuum fitting errors that affect the entire
forest, leading to spurious correlations which bias our measurement. In each bin A, our model correlation $\xi^{\text{mod}}$ is evaluated at the weighted (by the number of pixel pairs) mean separation ($r_{\parallel}$, $r_{\perp}$) of the Ly$\alpha$ pixels (equation (3.1)) in our data.

We measure the auto-correlation from [0, 200]$h^{-1}$Mpc using 50 bins of $4h^{-1}$Mpc along and perpendicular to the line of sight, giving us 2500 bins in total. The given limits on separation mean that the Ly$\alpha$ auto-correlation has $3.75 \times 10^{11}$ pixel pairs. In figure 4 we present our Ly$\alpha$ auto-correlation measurement as a function of ($\mu$, $r$) where $\mu = r_{\parallel}/r$ and $r^2 = r_{\parallel}^2 + r_{\perp}^2$, for EDR+M2 and eBOSS DR16. We show the auto-correlation in 4 wedges, computed by averaging the 2D correlation in different selections of $\mu$, ranging from closest ($\mu \in [0.95, 1]$) to furthest ($\mu \in [0, 0.5]$) from the line of sight. In general, our DESI EDR+M2 measurement is consistent with dMdB20, and across all bins ($r_{\parallel}$, $r_{\perp}$), the uncertainties are only ~1.9 times larger. It should be noted that there is overlap between the two datasets, since DESI is re-observing quasars observed in dMdB20. We present a constraint on the BAO peak amplitude from our best-fit model in section 5. For $\mu \in [0.95, 1]$, we can also see peaks caused by other atomic transitions from elements like Si, referred to as metal absorption or contamination in the rest of the paper. For example, we see the SiIII(1207) peak at $r \sim 20h^{-1}$Mpc and the SiII(1190)/SiII(1193) peak at $r \sim 60h^{-1}$Mpc. At $r \sim 100h^{-1}$Mpc we observe the BAO peak, which appears more prominently in the $\mu \in [0.95, 1]$ plot, due to its overlap with the SiII(1260)xLy$\alpha$ peak at $111h^{-1}$Mpc. In section 4.3.4, we discuss the origin of the metal peaks and how we model them.

## 3.4 Cross-correlation

Similarly to the auto-correlation, we define an estimator for the cross-correlation as [22, 61]:

$$ \xi_A = \frac{\sum_{i,j \in A} w_i w_j \delta_i}{\sum_{i,j \in A} w_i w_j}, \qquad (3.8) $$

for a pixel of Ly$\alpha$ absorption i and quasar j. In this case, the weights $w_j$ are corrected for the quasar bias evolution and are given by:

$$ w_j = \left( \frac{1 + z_j}{1 + z_{\text{fid}}} \right)^{\gamma_q - 1}, \qquad (3.9) $$

where $\gamma_q = 1.44 \pm 0.08$ [70]. The $(1 + z_{\text{fid}})$ term, with $z_{\text{fid}} = 2.25$, cancels again in equation (3.8). The sum in equation (3.8) applies over all quasar-pixel pairs except pixels from their background quasar. Again, the correlation is binned in terms of line of sight and transverse separation ($r_{\parallel}$, $r_{\perp}$), with $r_{\perp} \in [0, 200]$h$^{-1}$Mpc but now $r_{\parallel} \in [-200, 200]$h$^{-1}$Mpc. The cross-correlation is not symmetric under permutation of the two tracers, primarily because peaks from metal contamination appear at either $r_{\parallel} > 0$ or $r_{\parallel} < 0$, and because of systematic redshift errors (equation (4.9)) and because the bias of each tracer evolves differently with redshift. Therefore we choose to define -ve separations for the case when the quasar is behind the Ly$\alpha$ pixel ($z_q > z_{Ly\alpha}$) and likewise +ve separations when the quasar is between the observer and the Ly$\alpha$ pixel ($z_{Ly\alpha} > z_q$). This asymmetry allows us to study systematics like quasar redshift errors [33]. The given range of separation in both directions and the $4h^{-1}$Mpc bin width now result in $N = 100 \times 50 = 5000$ bins.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** The DESI EDR+M2 (blue points) and eBOSS DR16 [22] (shaded green) Lyα auto-correlation compressed into weighted averages of $\mu = r_\parallel/r$, where $r = \sqrt{r_\parallel^2 + r_\perp^2}$. We also include the best-fitting model to EDR+M2 described in section 4 (blue dashed). We have multiplied the correlation by $r^2$ to visualise the BAO peak, which is visible in both data sets. From these plots, we can see the consistency between our measurements and those in eBOSS DR16 — a validation of the quality of DESI data at an early survey stage. Note that the presence of three other bumps in the line of sight plot (bottom right) at 20, 60 and 111h$^{-1}$Mpc is due to correlations between the Lyα and SiII(1207), SiII(1190)/SiII(1193) and SiII(1260) lines respectively. We model these contributions to the correlation in section 4.3.4.</figcaption></figure>

In figure 5 we again show the cross-correlation as a function of $r$ and $\mu$, this time averaging over $\mu \in [-1, 1]$, since the cross-correlation has -ve values of $r_\parallel$. Here, the correlation is reversed with respect to the matter correlation function. This behaviour arises from the fact that the bias of the Lyα forest is negative. By convention we trace the flux-transmission field, which is higher in more under-dense regions — and the quasar bias is positive, giving a negative product. The uncertainties in our DESI EDR+M2 cross-correlation measurement across all ($r_\parallel, r_\perp$) bins are ~1.7 times larger than in eBOSS DR16. In the following sections, we discuss the covariance of the measurements of both the auto- and cross-correlations.

## 3.5 Covariance matrix

We estimate the covariances of the auto- and cross-correlations using a sub-sampling method, following [19]. In this method, we define our sub-samples in 440 HEALPix [71] pixels using Healpy[72], equivalent to an nside = 16, each of resolution 3.7 deg or solid angle $3.7 \times 3.7 = 13.4$ deg$^2$. We then calculate the weighted covariance $C_{AB}$ between samples as:

$$C_{AB} = \frac{1}{W_A W_B} \sum_{i,j \in A} \sum_{k,l \in B} W_A^s W_B^s [\hat{\xi}_A^s \hat{\xi}_B^s - \hat{\xi}_A \hat{\xi}_B], \quad (3.10)$$
<figure><img src="image_7.png" /><figcaption>**Figure 5.** The 3D EDR+M2 (blue points) and eBOSS DR16 [22] (green shaded regions) Lyα-quasar cross-correlation, and the baseline fit (blue dashed) to EDR+M2 described in section 4. Because we have negative values of $r_\parallel$ (when $z_q > z_{Ly\alpha}$), we have negative values of $\mu = r_\parallel/r$ and therefore average over $\mu \in [-1, 1]$. The cross-correlation is expectedly noisier than the auto-correlation, but still we see a good level of consistency between eBOSS and DESI at this early stage.</figcaption></figure>

where $\hat{\xi}_A^s$ is the measured correlation and $W_A^s$ is the summed weight in sub-sample s, and $W_A = \sum_s W_A^s$. We assign each Lyα pixel pair, or quasar-Lyα pair, to a given HEALPix pixel based on the object with the highest right ascension. Our effective redshift, $z_{eff} = 2.376$, is calculated by taking the average redshift of each bin of the correlation function weighted by the number of pixel pairs. The area of one HEALPix pixel at this redshift is $\sim 255 \times 255$ $(h^{-1}Mpc)^2$. The sub-sampling method of estimating covariance has also been validated in studies with mock catalogues in dMdB20.

Our auto- and cross-covariance matrices have $N = 2500 \times 2500 = 6.25$ million and $N = 5000 \times 5000 = 25$ million entries, with the diagonal entries dominating.

The normalised covariance, or correlation matrix, is defined as:

$$ Corr_{AB} = \frac{C_{AB}}{\sqrt{Var_A Var_B}}, \qquad (3.11) $$

where $Var_A = C_{AA}$ is the variance. The off-diagonal elements of equation (3.10) are noisy, and we need to smooth them to make our covariance matrices invertible. Following dMdB20, we model the off-diagonal elements of the correlation matrices as a function of line of sight and transverse separation, such that $Corr_{AB} = Corr(\Delta r_\parallel, \Delta r_\perp) = Corr(|r_\parallel^A - r_\parallel^B|, |r_\perp^A - r_\perp^B|)$. In figure 7 of dMdB20 we see that the correlation matrix in this model decreases rapidly with $(\Delta r_\parallel, \Delta r_\perp)$.
# 4 Correlation function models

In this section, we discuss the theoretical model for our correlation functions, which includes the effect of all the major contaminants. For the most part, our model follows that in dMdB20, except the instrumental systematic effects model, which we adapt to account for the differences between the DESI and eBOSS instruments. The final fit, which we describe in section 5, consists of 13 free parameters shown in table 1.

## 4.1 Power spectra

For our correlation model, we start with an isotropic linear power-spectrum template, decomposed into a peak and a smooth component. We then add the Lyα/quasar Kaiser terms, some non-linear corrections, and contaminants, and convert the resulting anisotropic power spectrum into a correlation function. Finally, we multiply by the distortion matrix and recombine the peak and smooth components to obtain the final correlation model.

The auto- and cross-correlation functions are derived from the tracer biased power-spectrum, given by:

$$P_{ij}(k, \mu_k, z) = b_i(z)b_j(z)(1 + \beta_i\mu_k^2)(1 + \beta_j\mu_k^2)P_{QL}(k, \mu_k, z)F_{NL,ij}(k, \mu_k)G(k, \mu_k), \quad (4.1)$$

where the vector $(k, \mu_k)$ is defined such that $k^2 = k_{\parallel}^2 + k_{\perp}^2$ and $\mu_k = k_{\parallel}/k$, G is a term that accounts for binning on an $(r_{\parallel}, r_{\perp})$ grid (equation (4.5)) and $P_{QL}$ is a quasi linear power spectrum that is multiplied by a non-linear correction term $F_{NL,ij}(k, \mu_k)$ (section 4.2). The parameters $b_i$ and $\beta_i$ are the linear bias and RSD of tracer i.

To fit a unique function across each $(r_{\parallel}, r_{\perp})$ bin where the effective redshift varies from that of the full dataset ($z_{eff} = 2.376$), we need to evolve the Lyα and quasar biases using:

$$b_{Ly\alpha}(z) = b_{Ly\alpha}(z_{eff}) \left( \frac{1+z}{1+z_{eff}} \right)^{\gamma_{Ly\alpha}}, \quad (4.2)$$

$$b_{q}(z) = b_{q}(z_{eff}) \left( \frac{1+z}{1+z_{eff}} \right)^{\gamma_{q}}, \quad (4.3)$$

where $\gamma_{Ly\alpha} = 2.9$ [13] and $\gamma_{q} = 1.44$ [70].

For the auto-correlation (i = j) the only tracer is Lyα absorption and for the cross (i≠j) the tracers are Lyα absorption and quasars. Following past BAO analyses (dMdB20), we take the linear matter power spectrum from Planck 2018⁵ [30] and separate it into a smooth component and peak component with non-linear broadening corrections:

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

where $P_{peak}$ contains BAO and $P_{smooth}$ does not. We make this decomposition using the side band technique from [18] such that $P_{peak} = P_{lin} - P_{smooth}$, and allows us to isolate the model for the BAO peak from the rest of the correlation and fit its position and amplitude. For the non-linear broadening correction of the BAO peak, we use $\Sigma_{\parallel} = 6.36h^{-1}\text{Mpc}$ and $\Sigma_{\perp} = 3.24h^{-1}\text{Mpc}$ as in dMdB20 [22, 73], assuming that our small difference in effective redshift has negligible effect.

<footnote n="5">⁵Derived using CAMB https://github.com/cmbant/CAMB.</footnote>
The last term in equation (4.1), G, corrects for the fact that we bin correlations in a grid in $(r_\parallel, r_\perp)$ and is the product of the Fourier transform of two top-hat functions:

$$G(k_\parallel, k_\perp) = \text{sinc}\left(\frac{k_\parallel \Delta A_\parallel}{2}\right) \text{sinc}\left(\frac{k_\perp \Delta A_\perp}{2}\right), \qquad (4.5)$$

where $\Delta A_\parallel = 4h^{-1}\text{Mpc}$ and $\Delta A_\perp = 4h^{-1}\text{Mpc}$ are the bin widths in the line-of-sight and transverse directions respectively. To compute equation (4.5), we make the approximation that our correlation is distributed homogeneously within each bin. In reality, it scales with $r_\perp$ in the transverse direction, but [20] have shown that the approximation made here produces an accurate correlation function.

## 4.2 Non-linear corrections

We multiply the quasi-linear power spectrum $P_{\text{QL,ij}}(k, \mu_k)$ of equation (4.1) by a non-linear correction $F_{\text{NL,ij}}(k, \mu_k)$. In the auto-correlation, the dominant effects encompassed by this correction are non-linear structure growth, peculiar motion, and broadening due to the thermal motion of particles, modelled in [9].

Statistical redshift errors and galaxy peculiar velocities of quasars smear the cross-correlation along the line of sight $r_\parallel$, which following [74] we model as Lorentzian damping:

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

where $\sigma_v$ is one of the free parameters in our model, given in table 1. Statistical errors also affect the auto-correlation by adding random shifts to quasar spectral templates that in turn introduce systematic error in the mean continuum estimate $\bar{F}(z)C_q(\lambda)$ (equation (2.2)). The impact of statistical redshift errors on BAO measurements from both the auto- and cross-correlation is described in [75].

## 4.3 Modeling contaminants

Both the auto- and cross-correlation functions receive small contributions from other effects, which for the auto-correlation we model as:

$$\xi_{\text{auto}} = \xi_{\text{Ly}\alpha\times\text{Ly}\alpha} + \sum_{\text{m}} \xi_{\text{Ly}\alpha\times\text{m}} + \sum_{\text{m}_i, \text{m}_j} \xi_{\text{m}_i\times\text{m}_j} + \xi_{\text{inst}}, \qquad (4.7)$$

where $\xi_{\text{Ly}\alpha\times\text{Ly}\alpha}$ is transformed from $P_{\text{Ly}\alpha}(k, \mu_k)$ (equation (4.1)) via a Fast Hankel transform. The Ly$\alpha$ bias parameter is also modelled to account for contamination by high column density (HCD, see section 4.3.3) systems. $\xi_{\text{Ly}\alpha\times\text{m}}$ is the cross-correlation between Ly$\alpha$ absorption and metal absorption and $\xi_{\text{m}_i\times\text{m}_j}$ is the metal-metal auto-correlation discussed in section 4.3.4. Finally, $\xi_{\text{inst}}$ is the contribution from DESI instrumental systematics, including sky subtraction residuals and spectro-photometric calibrations, described in full in [34] and modelled in section 4.3.2. Likewise, for the cross-correlation, we write:

$$\xi_{\text{cross}} = \xi_{\text{Ly}\alpha\times\text{q}} + \sum_{\text{m}} \xi_{\text{q}\times\text{m}} + \xi^{\text{TP}}, \qquad (4.8)$$

where $\xi_{\text{Ly}\alpha\times\text{q}}$ is transformed from Fourier space in the same way as the auto-correlation, $\xi_{\text{q}\times\text{m}}$ is the quasar-metal line cross-correlation and $\xi^{\text{TP}}$ is the contribution from the effect of quasar radiation given in section 4.3.1.
Finally we introduce another free parameter to our model to account for any systematic shift in line-of-sight separation between quasars and Lyα absorption:

$$ \Delta r_{\parallel,q} = r_{\parallel,\text{True}} - r_{\parallel,\text{Measured}}. \qquad (4.9) $$

This effect arises due to systematic quasar redshift errors and thus produces an asymmetry over positive and negative $r_\parallel$ which can be studied using the cross-correlation of quasars and Lyα absorption, as is done in [33]. In figure 6, we show the effect of separately excluding the quasar-metal cross-correlation model (left hand panel) and the quasar radiation and systematic redshift error models (right hand panel).

### 4.3.1 Quasar radiation effects

At small scales, quasar radiation can strongly affect its surrounding gas or a nearby Lyα forest. This radiation dominates over the UV background [25], so the ionisation fraction, and therefore the number of passing Lyα photons, increases. This is modelled following [22, 23], assuming isotropic quasar emission:

$$ \xi^{\text{TP}} = \xi_0^{\text{TP}} \left( \frac{1h^{-1}\text{Mpc}}{r} \right)^2 \exp \left( \frac{-r}{\lambda_{\text{UV}}} \right), \qquad (4.10) $$

where we fit for the scale parameter $\xi_0^{\text{TP}}$, and $\lambda_{\text{UV}}$ is the UV photon mean free path set to $300h^{-1}\text{Mpc}$ following [76].

### 4.3.2 Instrumental systematics

The focal plane of DESI is divided in 10 “petals”, with 500 fibers each that are fed to one of the 10 spectrographs. The calibration of DESI spectra (sky subtraction using sky fibers and flux calibration using standard stars) are performed separately for each of the 10 petals [50]. This introduces a possible source of systematic correlations for pairs of spectra obtained with the same spectrograph. The impact of correlated sky residuals and flux calibration errors in the Lyα correlations in DESI is described in detail in [34]. Poisson fluctuations at the observed wavelength of the sky lines within each spectrograph, and flux calibration residuals, produce an excess correlation at $r_\parallel = 0$. We model the contributions of both of these, which contribute roughly equally, in the auto-correlation with an empirical functional form:

$$ \xi_{\text{inst}} = \begin{cases} A_{\text{inst}} \left( \frac{r_{\perp}}{80} - 1 \right)^2 & r_{\perp} < 80h^{-1}\text{Mpc}, r_{\parallel} = 0 \\ 0 & r_{\perp} > 80h^{-1}\text{Mpc} \end{cases} \qquad (4.11) $$

where $A_{\text{inst}}$ is fit and given in section 5. The limit of $80h^{-1}\text{Mpc}$ roughly corresponds to the angular size of a DESI petal on the sky at the effective redshift of our data. Note that we do not see this effect in the cross-correlation since it requires two correlating subtracted sky lines in two Lyα forests. In dMdB20 the sky subtraction residuals were modelled with a Gaussian with two free parameters instead of one. The differences in the DESI and eBOSS instruments — such as eBOSS having only two spectrographs with 500 fibres each — lead to a difference in the shape and amplitude of the instrumental effects model. In section 5, we compare the amplitude of both models at $r_\parallel = r_\perp = 0$.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** The Lyα-quasar cross-correlation as a function of $r_\parallel$, for $r_\perp \in [0, 12]h^{-1}\text{Mpc}$. On the right, we over-plot our baseline model (red solid) reported in table 1 and a model with no metal line component (blue dashed). On the left, we show two additional lines to highlight the effect of setting $\Delta r_{\parallel,q}$ to 0 and ignoring quasar radiation effects (green dashed). We can see the asymmetry of the cross-correlation through metals and systematic redshift errors (parameterised by $\Delta r_{\parallel,q}$) in the data here.</figcaption></figure>

### 4.3.3 Absorption by high column density systems

The Lyα forest contains absorption from neutral hydrogen in the diffuse IGM and systems of neutral hydrogen with significantly higher column density.<fnref n="6" /> In section 2.2, we defined DLAs as having column density $N_{HI} > 10^{20.3}\text{cm}^{-2}$ and discussed how we identify and account for these using pixel masks and a Voigt profile to model the wings. Assuming we identify and treat all of the DLAs, the noise in our correlation measurements should be significantly reduced. We still expect to be contaminated by some objects of column density $N_{HI} < 10^{20.3}\text{cm}^{-2}$ (sometimes referred to as sub-DLAs) which have a similar absorption profile, but lesser effect on our data. Ideally, we would also mask HCDs below the DLA limit, but these are difficult to identify with our current pipeline and resolution. Following [61, 77] we model the contamination from these systems as a modification of the large-scale bias parameters:

$$ \tilde{b}_{\text{Ly}\alpha} = b_{\text{Ly}\alpha} + b_{\text{HCD}} F_{\text{HCD}}(k_{\parallel}), \qquad (4.12) $$

$$ \tilde{b}_{\text{Ly}\alpha} \tilde{\beta}_{\text{Ly}\alpha} = b_{\text{Ly}\alpha} \beta_{\text{Ly}\alpha} + b_{\text{HCD}} \beta_{\text{HCD}} F_{\text{HCD}}(k_{\parallel}), \qquad (4.13) $$

where $b_{\text{HCD}}$ and $\beta_{\text{HCD}}$ are the bias and RSD parameter of the HCD systems and $F_{\text{HCD}}(k_\parallel)$ is a factor depending on the distribution of HCDs, modelled following dMdB20 as:

$$ F_{\text{HCD}}(k_{\parallel}) = \exp(-L_{\text{HCD}} k_{\parallel}). \qquad (4.14) $$

<footnote n="6">⁶High column density systems is a term used in past eBOSS publications to refer to neutral hydrogen with column density $N_{HI} > 10^{17.2}\text{cm}^{-2}$, of which DLAs are the highest column density subset and the most impactful contaminants on our data.</footnote>
Here $L_{\text{HCD}}$ is the typical HCD length scale, and is set to $10h^{-1}\text{Mpc}$ following dMdB20 because of its degeneracies with other parameters in the fit. Note that it was verified in [78, 79] that varying $L_{\text{HCD}}$ between ~7–13$h^{-1}\text{Mpc}$ does not significantly affect BAO measurements.

### 4.3.4 Absorption by metal lines

As stated earlier and shown in equations (4.7) and (4.8), our measured correlations are not only the result of Lyα absorption in the IGM or HCDs but also heavier elements in the IGM due to galactic gas outflows [80]. Therefore, we introduce a cross power spectrum $P_{\text{mn}}(\mathbf{k}, z)$ for each pair of absorbers (m,n) and their configuration-space counterparts $\xi^{\text{m}\times\text{n}}$.

Because we measure our flux-transmission field and correlations assuming all absorption in the Lyα is due to the Lyα transition, we assign to each absorption pixel a redshift $z = \lambda_{\text{obs}}/1215.67\AA - 1$. However, because of the metal lines from heavier elements, an excess correlation appears at the separation between the true and assigned redshifts $z_{\text{m}}$ and $z_{\text{Lyα}}$:

$$r_{\parallel} = d_{\text{m}} - d_{\text{Ly}\alpha} = c \int_{z_{\text{Ly}\alpha}}^{z_{\text{m}}} \frac{dz'}{H(\bar{z})} \approx \frac{c(1+z) \lambda_{\text{m}} - \lambda_{\text{Ly}\alpha}}{H(z) \bar{\lambda}}, \quad (4.15)$$

where $\bar{\lambda}$ and $\bar{z}$ are the mean values of the Lyα and metal absorption. The excess correlation is more clearly observed for vanishing physical separations between two absorbers, when $r_{\perp} \sim 0$. We represent the correlation shift for each pair of absorbers (m,n) using the metal matrix formalism of [81]. In this method, the shifted correlation function $\xi_1^{\text{m}\times\text{n}}$ for each absorber pair (m,n) is given in terms of the un-shifted correlation function $\xi_0^{\text{m}\times\text{n}}$ as:

$$\xi_1^{\text{m}\times\text{n}}(A) = \sum_{\text{B}} M_{\text{AB}}, \xi_0^{\text{m}\times\text{n}}(r_{\parallel}(\text{B}), r_{\perp}(\text{B})), \quad (4.16)$$

where $(\text{m,n}) \in \text{B}$ is the pixel separation using the rest wavelength of each absorber m and n, $(\text{m,n}) \in \text{A}$ refers to the pixel separation assuming both pixels have rest wavelength $\lambda_{\text{Lyα}}$, and:

$$M_{\text{AB}} = \frac{\sum_{(\text{m,n}) \in \text{A}, (\text{m,n}) \in \text{B}} w_{\text{m}} w_{\text{n}}}{\sum_{(\text{m,n}) \in \text{A}} w_{\text{m}} w_{\text{n}}}. \quad (4.17)$$

For the cross-correlation, index m is a quasar with redshift $z_q$ used to calculate the projected separation. To compute $M_{\text{AB}}$ in a reasonable time, we only use 1% of pixel pairs. We have verified that this does not affect the measured metal biases.

Different rest-frame wavelengths peaks at different line of sight separations. Given the range of separations that we measure, we model 4 dominant metal lines — SiII(1190), SiII(1193), SiIII(1207), and SiII(1260) [20]. Each metal line has its own bias and RSD parameters ($b_{\text{m}}$, $\beta_{\text{m}}$) and corresponding power spectrum $P_{\text{mn}}$ as mentioned above. We cannot determine these two parameters separately since the correlations $\xi_1^{\text{m}\times\text{n}}$ only have a significant impact at small $r_{\perp}$, and therefore we set $\beta_{\text{m}} = 0.5$ following [20]. The values of $b_{\text{m}}$ for each metal line are constrained in the fits discussed in section 5. The CIV metal bias has also been typically fit in previous Lyα forest correlation studies [81, 82], but we decide not to model it in this study because it was weakly detected and had a negligible effect on the dMdB20 analysis, and we report the same for EDR+M2 data.
## 4.4 Modelling the continuum distortion

Due to distortions created in the quasar continuum fitting process, we apply two corrections to our measured delta fields in section 2.3. Doing this, each measured delta $\tilde{\delta}_q(\lambda)$ becomes a linear combination of all other pixels in the same forest, distorting the correlation functions measured from these deltas. To correct for the distortion we relate the measured $(\hat{\xi})$ and true $(\xi^t)$ correlations following [20, 22, 25] with:

$$ \hat{\xi}_A = \sum_{A'} D_{AA'} \xi_{A'}^t, \qquad (4.18) $$

where

$$ D_{AA',\text{auto}} = \frac{\sum_{(i,j) \in A} w_i w_j \left( \sum_{(i',j') \in A'} \eta_{ii'} \eta_{jj'} \right)}{\sum_{(i,j) \in A} w_i w_j}, \qquad (4.19) $$

for the auto-correlation, where $\eta$ is the projection matrix given in equation (3.3). We also have

$$ D_{AA',\text{cross}} = \frac{\sum_{(i,k) \in A} w_i \sum_{(j,k) \in A'} \eta_{ij}}{\sum_{(i,k) \in A} w_i}, \qquad (4.20) $$

for the cross-correlation, where i, j are pixels in the same forest and k is a quasar. Like the metal matrices, we compute the distortion matrices with only 1% of total pixel pairs to reduce computing hours. In our analysis, we model the effect of distortion out to a separation of $300h^{-1}$Mpc, since distortions at this scale are known to affect the correlation function at scales relevant to our analysis. The two expressions above equate to multiplying our physical model of the correlation by the distortion matrix before any comparison with data. The effects of distortion have also been validated in previous studies using mocks [20].

## 5 Fits to the data

We now present the best-fitting solution of the model introduced in section 4 to the $\Lambda$CDM auto-correlation and $\Lambda$CDM quasar cross-correlation described in section 3. For this study, we have established a baseline model with 13 free parameters based on Kaiser models for biased tracers with additional contaminants, described in detail in section 4.

In the anisotropic correlation model of equation (4.1), it is evident that the cross-correlation we measure is only sensitive to the product of the $\Lambda$CDM and quasar biases. As a result, we present only the auto and auto+cross combined fits. In the latter, we break the bias degeneracy to measure the quasar bias $b_q$ distinctly. Then, using the fact that for point tracers $\beta_q = f/b_q$, where $f = 0.97$ [30] is the growth-rate for our fiducial cosmology at $z = 2.376$, we infer a value of $\beta_q$. The other advantage of combining both correlations is improving constraints on key parameters. For example, our constraints on $A_{\text{BAO}}$ improve by 1.6x going from auto-correlation fit alone to the combined fit.

In table 1, mean values and 68% credible regions of the posterior distributions for each free parameter in the auto and combined auto+cross correlations are given.<fnref n="7" /> To perform

<footnote n="7">⁷Computed using the plotting tool `Getdist` [83].</footnote>
<table><caption>**Table 1.** Table of early DESI fits on the auto-correlation and the combined auto+cross correlations. We do not perform a cross-correlation-only fit due to degeneracies between parameters. In the first part we show the mean and 68% credible region for each free parameter in our model computed from the posterior distributions from the nested sampler. We evaluate each parameter at the effective redshift of our data $z_{eff} = 2.376$, where redshift is relevant. Where a parameter does not feature in the auto-correlation fit, we use the “-” symbol. In the second part, we give the characteristics of the fit, including the number of bins, free parameters, and the $\chi^2_{min}$ probability as estimated with a minimiser.</caption><thead><tr><th>Parameter</th><th>Auto</th><th>Combined</th></tr></thead><tbody><tr><td>b<sub>Lyα</sub></td><td>-0.129±0.010</td><td>-0.134±0.009</td></tr><tr><td>β<sub>Lyα</sub></td><td>1.49<sup>+0.14</sup><sub>-0.19</sub></td><td>1.41<sup>+0.12</sup><sub>-0.15</sub></td></tr><tr><td>b<sub>HCD</sub></td><td>-0.045±0.009</td><td>-0.039±0.009</td></tr><tr><td>10<sup>3</sup>b<sub>η,SiII(1190)</sub></td><td>-3.3±1.0</td><td>-2.2±0.8</td></tr><tr><td>10<sup>3</sup>b<sub>η,SiII(1193)</sub></td><td>-1.7<sup>+1.1</sup><sub>-0.8</sub></td><td>-0.9<sup>+0.8</sup><sub>-0.3</sub></td></tr><tr><td>10<sup>3</sup>b<sub>η,SiII(1260)</sub></td><td>-2.9±1.1</td><td>-2.6±0.9</td></tr><tr><td>10<sup>3</sup>b<sub>η,SiIII(1207)</sub></td><td>-3.9±1.0</td><td>-3.4±0.9</td></tr><tr><td>b<sub>q</sub></td><td>-</td><td>3.41±0.16</td></tr><tr><td>Δr<sub>‖,q</sub>(h<sup>-1</sup>Mpc)</td><td>-</td><td>-2.21±0.18</td></tr><tr><td>σ<sub>v</sub>(h<sup>-1</sup>Mpc)</td><td>-</td><td>5.2±0.5</td></tr><tr><td>ξ<sup>TP</sup><sub>0</sub></td><td>-</td><td>0.68±0.18</td></tr><tr><td>10<sup>4</sup>A<sub>inst</sub></td><td>2.6<sup>+0.4</sup><sub>-0.8</sub></td><td>2.4<sup>+0.3</sup><sub>-0.5</sub></td></tr><tr><td>A<sub>BAO</sub></td><td>1.04±0.45</td><td>1.21±0.32</td></tr><tr><td>N<sub>bin</sub></td><td>1590</td><td>4770</td></tr><tr><td>N<sub>params</sub></td><td>9</td><td>13</td></tr><tr><td>χ<sup>2</sup><sub>min</sub></td><td>1711</td><td>4945</td></tr><tr><td>probability</td><td>0.01</td><td>0.03</td></tr></tbody></table>

each fit, we use the Python package Vega,<fnref n="8" /> which employs the nested sampler Polychord [84] to sample each parameter. We find that the auto-only results are very consistent with the auto+cross combined. In each column we put a “-” symbol in entries where parameters are only relevant to the combined correlation and not to the auto-correlation alone.

The second part of table 1 gives the characteristics of each fit, including the number of bins in each measurement, the number of free parameters, the minimum $\chi^2$ value, and the probability (p-value). Note that to obtain the latter two values, we use a minimiser instead of the nested sampler used to obtain parameter means and confidence regions. The minimiser makes use of the Python package Iminuit<fnref n="9" /> We perform the fit for $r \in [10, 180]h^{-1}Mpc$ due to the increasing complexity of modelling the correlation function at small scales (non-linearities). The $4h^{-1}Mpc$ bin width of the correlation functions translates to 1590 and 4770 bins for the auto and combined fits.

<footnote n="8">⁸https://github.com/andreicuceu/vega.</footnote>

<footnote n="9">⁹https://github.com/scikit-hep/iminuit.</footnote>
To measure the full correlation function we re-combine the peak and smooth parts of the correlation (outlined in section 4.1):

$$ \xi(r_{\parallel}, r_{\perp}, \alpha_{\parallel}, \alpha_{\perp}) = \xi_{\text{smooth}}(r_{\parallel}, r_{\perp}) + A_{\text{BAO}}\xi_{\text{peak}}(\alpha_{\parallel}r_{\parallel}, \alpha_{\perp}r_{\perp}), \quad (5.1) $$

where $\alpha_{\parallel}$ and $\alpha_{\perp}$ are the BAO parameters (see dMdB20 for recent constraints on these). Here $A_{\text{BAO}}$ is the amplitude of BAO, which we detect at $3.8\sigma$ confidence (table 1). We fix both $\alpha_{\parallel}$ and $\alpha_{\perp}$ to 1, or the fiducial cosmology, and focus on preparing for a robust measurement with DESI year 1 data.

The best-fitting model is shown on four auto-correlation wedges in figure 4 and the same four cross-correlation wedges in figure 5. In figure 7, we also show the auto-correlation along the line of sight, averaged over $0 < r_{\perp} < 12h^{-1}\text{Mpc}$ and overlaid with the baseline model from the minimiser fit. We also include a fit that does not model the effect of metal contamination, and a fit that does not model metal contamination or BAO. We multiply by $r$ to see that the bumps caused by the metal lines SiIII(1207) and SiII(1190)/SiII(1193) at $r \sim 20h^{-1}\text{Mpc}$ and $r \sim 60h^{-1}\text{Mpc}$ respectively are well-fit by our model, and poorly-fit by the no metal model. The SiII(1260) line gives a peak at $r \sim 111h^{-1}\text{Mpc}$, which along the line of sight dominates over the BAO peak with which it overlaps.

## 5.1 Fit probability

The combined fit presented in table 1 has a $\chi^2$ probability of 0.03. Because this is a relatively low value, we discuss several caveats here.

In figure 8, the residuals are spread evenly across $(r_{\parallel}, r_{\perp})$. There are some groups of bins with higher residuals closer to the line-of-sight ($r_{\perp} = 0h^{-1}\text{Mpc}$) where continuum distortions are strongest, and at larger separations $r \sim 150h^{-1}\text{Mpc}$, but these are likely due to the high degree of correlation between data points which we ignore in this plot (we use the diagonal of the covariance matrix here).

If we change the lower limit on our fit of both correlations to $r = 20h^{-1}\text{Mpc}$ rather than $10h^{-1}\text{Mpc}$, our fit probability increases to 0.08 from 0.03. It's possible that our constraints on the baseline model are slightly biased when fitting down to $10h^{-1}\text{Mpc}$, but we will leave this mock catalogues studies in DESI year 1. Recent results using Ly$\alpha$ mocks on AbacusSummit [85] also suggest that the Kaiser model that we employ here does not adequately model non-linear effects in the Ly$\alpha$ quasar cross-correlation for $r < 30h^{-1}\text{Mpc}$. If we set this as a lower limit on the cross-correlation while keeping $10h^{-1}\text{Mpc}$ for the auto-correlation ([85] find the Kaiser model is a good approximation in this case), we increase our probability to 0.05, without biasing our measurements.

Changing the upper limit to $r = 150h^{-1}\text{Mpc}$ also increases the fit probability slightly to 0.05, and has a negligible effect on measured parameter values. For all of these tests, the most important thing for future analyses will be the effect on the BAO scale parameters, which we do not measure here.<fnref n="10" /> Therefore, we keep the baseline analysis between $r \in [10, 180]h^{-1}\text{Mpc}$ as in dMdB20.

<footnote n="10">¹⁰The Ly$\alpha$ working group key paper on DESI year 1 data will be centred around measuring $(\alpha_{\parallel}, \alpha_{\perp})$, but studies of the Alcock-Paczyński (AP) effect using both the smooth and peak part of the Ly$\alpha$ forest correlations will happen.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** The best fitting solution (red solid) of our model to the DESI EDR+M2 Lyα auto-correlation (black points) function along the line of sight, where we average over $0 < r_{\perp} < 12h^{-1}\text{Mpc}$. We also include a fit which does not model metal contamination (blue dashed) and a fit which does not model metal contamination or BAO (green dotted). We multiply the y-axis by $r$ for visualisation. The bumps in correlation caused by SiIII(1207) and SiII(1190)/SiII(1193) at $r_{\parallel} \sim 20h^{-1}\text{Mpc}$ and $r_{\parallel} \sim 60h^{-1}\text{Mpc}$ are visible. The most prominent peak is driven by SiII(1260) at $r_{\parallel} \sim 111h^{-1}\text{Mpc}$, rather than the overlapping BAO peak at $r_{\parallel} = 100h^{-1}\text{Mpc}$.</figcaption></figure>

In section 2 we outlined how we compute the Lyα flux transmission catalogue using the full resolution of the DESI spectrograph (0.8Å). Using the optimal weighting of equation (3.1), by setting $\sigma_{\text{mod}}^2 = 7.5$, we reach the same precision on the correlation function as using 2.4Å pixels with $\sigma_{\text{mod}}^2 = 3$. However, when computing the fits from the latter weighting, we find a $\chi^2$ probability of $p = 0.13$ for the combined fit. Despite this higher fit probability, we prefer to use the full resolution of the DESI spectrograph in this analysis and leave this effect for study in the DESI year 1 analysis. Reducing the maximum rest-frame wavelength of the Lyα forest to 1200 Å from 1205 Å (chosen to minimise correlation function uncertainties [31]) increases the fit probability to $p = 0.18$, but we refrain from making this change until there is a further study on mock data sets with greater statistical power. Several sources of contamination like DLAs, metals, and other effects such as quasar radiation are modelled in the same way as in dMdB20 and require an in-depth study for DESI.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Residuals of our fit to the DESI EDR+M2 auto-correlation, defined as the difference between data and model correlation divided by the uncertainties of each data point (square root of the diagonal of the covariance matrix). The area between the black dashed lines indicates the BAO region between 80–120h$^{-1}$Mpc. We see that overall there are no strong areas of discrepancy.</figcaption></figure>

## 5.2 Comparison with eBOSS DR16

To compare with dMdB20, we fit their measured correlation functions with the settings described in this paper. Except for the instrumental systematics effect, every other part of the model is the same. In figure 9 we present the marginalized 2D posterior of $\beta_{Ly\alpha}$ and $b_{Ly\alpha}$. Note that the effective redshift of the eBOSS combined correlation is $z = 2.334$, and thus we use the Ly$\alpha$ bias-redshift relation in equation (4.2) to shift our bias measurement to that effective redshift.

We summarise the differences in the mean value, 68% confidence intervals, and the consistency between $b_{Ly\alpha}$, $\beta_{Ly\alpha}$ and a selection of parameters between DESI EDR+M2 and the eBOSS DR16 fits in table 2. In the case of $b_q$, we also shift the DESI EDR+M2 measurement to the effective redshift of the DR16 data set, this time using the relation $b_q(z) \propto (1+z)^\gamma$ where $\gamma \sim 1.44$ [70]. For the other two parameters in table 2, we are only interested in the mean value and the relative values of the uncertainties between surveys. The $\sigma_v$ constraint tells us that the DESI pipeline more precisely estimates quasar redshifts than was done by eBOSS DR16.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** 68% and 95% credible regions on bias and redshift-space distortions from the auto- and cross-correlation of the Lyα forest for DESI EDR+M2 data at the effective redshift of eBOSS DR16 (green) and eBOSS DR16 data. We performed fits on both datasets with the same model (presented in section 4). To evolve the value of the DESI EDR+M2 bias with redshift we use the relation $b_{\text{Ly}\alpha} \propto (1+z)^{\gamma-1}$ [13] where $\gamma = 2.9$ [22].</figcaption></figure>

The final tracers that we constrain in our model are HCDs and metals. We measure the bias of HCDs to be $b_{\text{HCD}} = -0.039 \pm 0.009$ and $b_{\text{HCD}} = -0.053 \pm 0.0043$ in early DESI and eBOSS (data) respectively. A smaller value of $b_{\text{HCD}}$ indicates weaker contamination by DLAs in our data, from which we conclude that the DESI pipeline [67] masks a higher proportion of DLAs, than the dMdB20 pipeline. The metals are all detected with high confidence, each with a non-zero bias to a minimum of $2.25\sigma$ except SiII(1193), and are in good agreement with dMdB20. The effect of metal correlations is visible in figure 7.

The systematic error in the quasar redshift estimations is parameterised by $\Delta r_{\parallel,q}$, which is measured through an asymmetry in the cross-correlation as described in the previous section. We report a value of $\Delta r_{\parallel,q} = -2.21 \pm 0.18h^{-1}\text{Mpc}$, a clear detection of systematic redshift errors. This offset could be related to an inconsistency in the treatment of the mean flux transmission $\bar{F}(z)$ in the quasar templates used by Redrock at $z \gtrsim 2$. In order to avoid this type of bias, dMdB20 used a redshift estimator which did not use wavelengths around the Lyα emission line or lower. They measured a value of $\Delta r_{\parallel,q} = 0.10 \pm 0.11h^{-1}\text{Mpc}$, consistent with zero. The impact of redshift errors in DESI quasars is studied in more detail in [33], and mitigation strategies are currently being considered in DESI.

At $r_{\perp} = r_{\parallel} = 0$, the spurious correlation produced by instrumental effects is $A_{\text{inst}} = 2.4_{-0.5}^{+0.3} \times 10^4$. In dMdB20, they used a Gaussian instrumental effects model, in which the spurious correlation at $r_{\perp} = r_{\parallel} = 0$ is $A_{\text{sky}}/\sigma_{\text{sky}}\sqrt{2\pi} = 1 \times 10^4$. We expect to see discrepancies
<table><caption>**Table 2.** Table highlighting the main sampled posteriors of the auto+cross correlations for DESI EDR+M2 and eBOSS DR16 data, their relative uncertainties (68% credible region), and their consistency with each other, defined as $\Delta\sigma = |\mu_{eBOSS} - \mu_{DESI}|/\sqrt{\sigma_{eBOSS}^2 + \sigma_{DESI}^2}$, where $\mu$ is the mean parameter value. We use the same model (outlined in section 4) for the fits in both cases, except for instrumental systematics contribution that is slightly different because of variations in the eBOSS and DESI instruments [34].</caption><thead><tr><th>Parameter</th><th>DESI EDR+M2</th><th>eBOSS DR16</th><th>σ<sub>DESI</sub>/σ<sub>eBOSS</sub></th><th>Δσ</th></tr></thead><tbody><tr><td>b<sub>Lyα</sub> (z = z<sub>eff,DR16</sub>)</td><td>-0.131±0.009</td><td>-0.117±0.004</td><td>2.25</td><td>1.39</td></tr><tr><td>β<sub>Lyα</sub></td><td>1.41<sup>+0.12</sup><sub>-0.15</sub></td><td>1.71±0.09</td><td>1.5</td><td>1.88</td></tr><tr><td>b<sub>q</sub> (z = z<sub>eff,DR16</sub>)</td><td>3.36±0.16</td><td>3.70±0.10</td><td>1.6</td><td>1.8</td></tr><tr><td>σ<sub>v</sub>[h<sup>-1</sup>Mpc]</td><td>5.2±0.5</td><td>6.77±0.29</td><td>1.72</td><td>-</td></tr><tr><td>A<sub>BAO</sub></td><td>1.21±0.32</td><td>1.16±0.16</td><td>2</td><td>-</td></tr></tbody></table>

in amplitudes here, given how distinct the two instruments are (section 4). To summarise the differences between early DESI and eBOSS DR16 data discussed above, we include a 13-parameter corner plot for both data sets in appendix A, figure 10.

## 6 Conclusion

We present the first study of the Lyα auto-correlation and its cross-correlation with quasars from Dark Energy Spectroscopic Instrument (DESI) data. Our data sample, EDR+M2, consists of 318 691 quasar target spectra from the DESI survey validation [52] phase and the first two months of the DESI main survey, resulting in 88 509 Lyα forests with redshift $z > 2$.

We use the catalogue of Lyα fluctuations from [31], which made several adaptations to previous analysis pipelines [22], leading to a 20% and 10% improvement in the auto- and cross-correlations. Our 2-point correlation functions are estimated following the eBOSS DR16 analysis [22], with which we report high consistency and an average of 1.7x larger uncertainties across all bins in co-moving separation ($r_\|, r_\perp$). We do a combined auto+cross correlation fit with our 13 free parameter model, based on linear perturbation theory, and make constraints on key parameters like the Lyα forest bias, quasar bias, and the amplitude of the BAO peak. We detect Baryon Acoustic Oscillations to $3.8\sigma$ confidence, showing the constraining power of DESI data at an early survey stage. The model used is again similar to that in [22], except the instrumental systematics model [34] that is modified to account for the differences between both instruments.

At this stage of the DESI survey, we have highlighted the quality of our data and validated existing analysis methods of 3D correlations using the Lyα forest. With relatively low statistical power, we do not report constraints on the scale parameters of BAO and instead leave that to the DESI year 1 analysis. This future work will require an in-depth study of systematics and contaminants like DLAs and metals before performing a proper cosmological analysis of BAO.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Baseline 13 free parameter fit on DESI EDR and eBOSS DR16 auto+cross correlations. We model the instrumental systematics effects of eBOSS with a slightly different model with two parameters (as opposed to just $A_{\text{DESI}}^{\text{sky}}$ for DESI), which we have not included here. For $b_{\text{Ly}\alpha}$ and $b_q$ we present posteriors for both surveys at the effective redshift of the eBOSS DR16 data set.</figcaption></figure>

A Complete fit results

Here we show the triangle plot of the full posterior distribution — made using the Polychord nested sampler [84] and the plotting tool `Getdist` [83] — of all 13 parameters in our model of the combined auto+cross correlations, for both the DESI EDR+M2 and eBOSS DR16 data set. We measure both $b_{\text{Ly}\alpha}$ and $b_q$ at the effective redshift of eBOSS DR16 by using the redshift evolution of each parameter on the EDR+M2 result. For the Ly$\alpha$ and quasars this is $b \propto (1+z)^\gamma$, where $\gamma_{\text{Ly}\alpha} = 2.9$ [13, 22] and $\gamma_{\text{quasar}} = 1.44$ [70].
We can see several parameters in the model which are correlated. The bias and RSD parameters of the Lyα forest are strongly correlated with each other and with the bias factor of HCDs. The correlations that we show between the quasar bias and the bias of HCDs, and the quasar radiation strength are why we do not fit the cross-correlation independently of the auto-correlation, as mentioned in section 5.

Finally, there is an apparent correlation between the instrumental systematic effects parameter $A_{\text{inst}}$ and the BAO peak amplitude $A_{\text{BAO}}$. This is somewhat unexpected because the instrumental systematics model should only contribute at small scales ($r_{\parallel} = 0, r_{\perp} < 40h^{-1}\text{Mpc}$) and needs to be studied further for future analyses of DESI data.

For each parameter, we see the EDR+M2 result is generally within $1-2\sigma$ of the DR16 results, except $\Delta r_{\parallel,q}$, which has a much larger absolute value for our data set. As mentioned in section 5, this is likely because of systematic redshift errors that arise from poorly performing quasar templates. Improvements in the quasar templates used for redshift estimation will be made [54] for DESI year 1 to reduce this effect, but for now, it is a well-constrained nuisance parameter in our model. Finally, we see that for EDR+M2, the posterior of $\text{SiII}(1193)$ is not Gaussian (it is hitting the upper prior bound of 0) because of difficulties constraining the parameter from our data. From the DR16 fits, we can also see that this is one of the weaker detected metals (although it is still non-zero with $3\sigma$ confidence). Thus we expect that, with lower statistical power, it is difficult to detect. DESI year 1 data will present an exciting opportunity to make the most precise constraints on the presence of metals in the Lyα forest, including lines that are not included in our model here but were in the past (CIV(eff) [22]).

## Acknowledgments

CG is partially supported by the Spanish Ministry of Science and Innovation (MICINN) under grants PGC-2018-094773-B-C31 and SEV-2016-0588. AC acknowledges support provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51526.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. Both JCM and AFR from the European Union's Horizon Europe research and innovation programme (COSMO-LYA, grant agreement 101044612). AFR acknowledges financial support from the Spanish Ministry of Science and Innovation under the Ramon y Cajal program (RYC-2018-025210) and the PGC2021-123012NB-C41 project. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya. AXGM acknowledges support from Dirección de Apoyo a la Investigación y al Posgrado, Universidad de Guanajuato, research Grant No. 179/2023 and CONACYT México under Grants No. 286897 and A1-S-17899.

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 Science and Technology of Mexico (CONACYT); 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.

The 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.

**Data points.** The data points corresponding to each figure in this paper can be accessed in the Zenodo repository at https://doi.org/10.5281/zenodo.8244702.

## References

[1] R. Lynds, *The Absorption-Line Spectrum of $^4c$ 05.34*, Astrophys. J. 164 (1971) L73.

[2] D.H. Weinberg, R. Dave, N. Katz and J.A. Kollmeier, *The Lyman-$\alpha$ forest as a cosmological tool*, AIP Conf. Proc. 666 (2003) 157 [astro-ph/0301186] [INSPIRE].

[3] EBOSS collaboration, *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 [arXiv:2007.08991] [INSPIRE].

[4] R.-Y. Cen, J. Miralda-Escudé, J.P. Ostriker and M. Rauch, *Gravitational collapse of small scale structure as the origin of the Lyman $\alpha$ forest*, Astrophys. J. Lett. 437 (1994) L9 [astro-ph/9409017] [INSPIRE].

[5] J. Miralda-Escudé, R.-Y. Cen, J.P. Ostriker and M. Rauch, *The Lyman $\alpha$ forest from gravitational collapse in the CDM + Lambda Model*, Astrophys. J. 471 (1996) 582 [astro-ph/9511013] [INSPIRE].

[6] L. Hui and N.Y. Gnedin, *Equation of state of the photoionized intergalactic medium*, Mon. Not. Roy. Astron. Soc. 292 (1997) 27 [astro-ph/9612232] [INSPIRE].

[7] A. Slosar et al., *The Lyman-$\alpha$ forest in three dimensions: measurements of large scale flux correlations from BOSS 1st-year data*, JCAP 09 (2011) 001 [arXiv:1104.5244] [INSPIRE].

[8] P. McDonald, *Toward a Measurement of the Cosmological Geometry at $z \sim 2$: Predicting Lyα Forest Correlation in Three Dimensions and the Potential of Future Data Sets*, Astrophys. J. 585 (2003) 34 [astro-ph/0108064] [INSPIRE].

[9] A. Arinyo-i-Prats, J. Miralda-Escudé, M. Viel and R. Cen, *The Non-Linear Power Spectrum of the Lyman Alpha Forest*, JCAP 12 (2015) 017 [arXiv:1506.04519] [INSPIRE].

[10] J.J. Givans et al., *Non-linearities in the Lyman-$\alpha$ forest and in its cross-correlation with dark matter halos*, JCAP 09 (2022) 070 [arXiv:2205.00962] [INSPIRE].

[11] R.A.C. Croft, D.H. Weinberg, N. Katz and L. Hernquist, *Recovery of the Power Spectrum of Mass Fluctuations from Observations of the Lyα Forest*, Astrophys. J. 495 (1998) 44 [astro-ph/9708018] [INSPIRE].

[12] P. McDonald et al., *The observed probability distribution function, power spectrum, and correlation function of the transmitted flux in the Lyα forest*, Astrophys. J. 543 (2000) 1 [astro-ph/9911196] [INSPIRE].
[13] SDSS collaboration, *The Ly$\alpha$ forest power spectrum from the sloan digital sky survey*, Astrophys. J. Suppl. 163 (2006) 80 [astro-ph/0405013] [INSPIRE].

[14] N. Palanque-Delabrouille et al., *The one-dimensional Ly$\alpha$ forest power spectrum from BOSS*, Astron. Astrophys. 559 (2013) A85 [arXiv:1306.5896] [INSPIRE].

[15] S. Chabanier et al., *The one-dimensional power spectrum from the SDSS DR14 Ly$\alpha$ forests*, JCAP 07 (2019) 017 [arXiv:1812.03554] [INSPIRE].

[16] N.G. Busca et al., *Baryon Acoustic Oscillations in the Ly$\alpha$ forest of BOSS quasars*, Astron. Astrophys. 552 (2013) A96 [arXiv:1211.2616] [INSPIRE].

[17] A. Slosar et al., *Measurement of Baryon Acoustic Oscillations in the Lyman-$\alpha$ Forest Fluctuations in BOSS Data Release 9*, JCAP 04 (2013) 026 [arXiv:1301.3459] [INSPIRE].

[18] D. Kirkby et al., *Fitting Methods for Baryon Acoustic Oscillations in the Lyman-$\alpha$ Forest Fluctuations in BOSS Data Release 9*, JCAP 03 (2013) 024 [arXiv:1301.3456] [INSPIRE].

[19] BOSS collaboration, *Baryon acoustic oscillations in the Ly$\alpha$ forest of BOSS DR11 quasars*, Astron. Astrophys. 574 (2015) A59 [arXiv:1404.1801] [INSPIRE].

[20] J.E. Bautista et al., *Measurement of baryon acoustic oscillation correlations at $z=2.3$ with SDSS DR12 Ly$\alpha$-Forests*, Astron. Astrophys. 603 (2017) A12 [arXiv:1702.00176] [INSPIRE].

[21] V. de Sainte Agathe et al., *Baryon acoustic oscillations at $z=2.34$ from the correlations of Ly$\alpha$ absorption in eBOSS DR14*, Astron. Astrophys. 629 (2019) A85 [arXiv:1904.03400] [INSPIRE].

[22] H. du Mas des Bourboux et al., *The Completed SDSS-IV Extended Baryon Oscillation Spectroscopic Survey: Baryon Acoustic Oscillations with Ly$\alpha$ Forests*, Astrophys. J. 901 (2020) 153 [arXiv:2007.08995] [INSPIRE].

[23] A. Font-Ribera et al., *The large-scale Quasar-Lyman $\alpha$ Forest Cross-Correlation from BOSS*, JCAP 05 (2013) 018 [arXiv:1303.1937] [INSPIRE].

[24] BOSS collaboration, *Quasar-Lyman $\alpha$ Forest Cross-Correlation from BOSS DR11: Baryon Acoustic Oscillations*, JCAP 05 (2014) 027 [arXiv:1311.1767] [INSPIRE].

[25] H. du Mas des Bourboux et al., *Baryon acoustic oscillations from the complete SDSS-III Ly$\alpha$-quasar cross-correlation function at $z=2.4$, Astron. Astrophys. 608 (2017) A130 [arXiv:1708.02225] [INSPIRE].

[26] M. Blomqvist et al., *Baryon acoustic oscillations from the cross-correlation of Ly$\alpha$ absorption and quasars in eBOSS DR14*, Astron. Astrophys. 629 (2019) A86 [arXiv:1904.03430] [INSPIRE].

[27] A. Cuceu et al., *Constraints on the Cosmic Expansion Rate at Redshift 2.3 from the Lyman-$\alpha$ Forest*, Phys. Rev. Lett. 130 (2023) 191003 [arXiv:2209.13942] [INSPIRE].

[28] DESI collaboration, *The DESI Experiment, a whitepaper for Snowmass 2013*, arXiv:1308.0847 [INSPIRE].

[29] DESI collaboration et al., *The DESI Experiment Part I: Science, Targeting, and Survey Design*, arXiv:1611.00036.

[30] PLANCK collaboration, *Planck 2018 results. VI. Cosmological parameters*, Astron. Astrophys. 641 (2020) A6 [Erratum ibid. 652 (2021) C4] [arXiv:1807.06209] [INSPIRE].

[31] DESI collaboration, *The Lyman-$\alpha$ forest catalog from the Dark Energy Spectroscopic Instrument Early Data Release*, arXiv:2306.06312 [INSPIRE].

[32] H. Herrera et al., in preparation (2023).

[33] A. Bault et al., in preparation (2023).
[34] S. Gontcho A Gontcho et al., in preparation (2023).

[35] DESI collaboration, Overview of the Instrumentation for the Dark Energy Spectroscopic Instrument, Astron. J. 164 (2022) 207 [arXiv:2205.10939] [INSPIRE].

[36] DESI collaboration, The DESI Experiment Part II: Instrument Design, arXiv:1611.00037 [INSPIRE].

[37] DESI collaboration, The Robotic Multiobject Focal Plane System of the Dark Energy Spectroscopic Instrument (DESI), Astron. J. 165 (2023) 9 [arXiv:2205.09014] [INSPIRE].

[38] DESI collaboration, The Optical Corrector for the Dark Energy Spectroscopic Instrument, arXiv:2306.06310 [INSPIRE].

[39] E. Chaussidon et al., Target Selection and Validation of DESI Quasars, Astrophys. J. 944 (2023) 107 [arXiv:2208.08511] [INSPIRE].

[40] C.H. Hahn et al., The DESI Bright Galaxy Survey: Final Target Selection, Design, and Validation, Astron. J. 165 (2023) 253 [arXiv:2208.08512] [INSPIRE].

[41] DESI collaboration, Target Selection and Validation of DESI Luminous Red Galaxies, Astron. J. 165 (2023) 58 [arXiv:2208.08515] [INSPIRE].

[42] A. Raichoor et al., Target Selection and Validation of DESI Emission Line Galaxies, Astron. J. 165 (2023) 126 [arXiv:2208.08513] [INSPIRE].

[43] DESI collaboration, Overview of the DESI Legacy Imaging Surveys, Astron. J. 157 (2019) 168 [arXiv:1804.08657] [INSPIRE].

[44] H. Zou et al., Project Overview of the Beijing-Arizona Sky Survey, Publ. Astron. Soc. Pac. 129 (2017) 064101 [arXiv:1702.03653].

[45] D. Schlegel et al., in preparation (2023).

[46] A.D. Myers et al., The Target-selection Pipeline for the Dark Energy Spectroscopic Instrument, Astron. J. 165 (2023) 50 [arXiv:2208.08518] [INSPIRE].

[47] A. Raichoor et al., in preparation (2023).

[48] DESI collaboration, Survey Operations for the Dark Energy Spectroscopic Instrument, arXiv:2306.06309 [INSPIRE].

[49] D. Kirkby et al., in preparation (2023).

[50] J. Guy et al., The Spectroscopic Data Processing Pipeline for the Dark Energy Spectroscopic Instrument, arXiv:2209.14482.

[51] DESI collaboration, The Early Data Release of the Dark Energy Spectroscopic Instrument, arXiv:2306.06308 [DOI:10.5281/zenodo.7964161] [INSPIRE].

[52] DESI collaboration, Validation of the Scientific Program for the Dark Energy Spectroscopic Instrument, arXiv:2306.06307 [DOI:10.5281/zenodo.7858207] [INSPIRE].

[53] D.M. Alexander et al., The DESI Survey Validation: Results from Visual Inspection of the Quasar Survey Spectra, arXiv:2208.08517.

[54] DESI collaboration, Performance of the Quasar Spectral Templates for the Dark Energy Spectroscopic Instrument, Astron. J. 166 (2023) 66 [arXiv:2305.10426] [INSPIRE].

[55] S. Bailey et al., in preparation (2023).

[56] N. Busca and C. Balland, QuasarNET: Human-level spectral classification and redshifting with Deep Neural Networks, arXiv:1808.09955 [INSPIRE].
[57] J. Farr, A. Font-Ribera and A. Pontzen, *Optimal strategies for identifying quasars in DESI*, JCAP 11 (2020) 015 [arXiv:2007.10348] [INSPIRE].

[58] Z. Guo and P. Martini, *Classification of Broad Absorption Line Quasars with a Convolutional Neural Network*, Astrophys. J. 879 (2019) 72 [arXiv:1901.04506] [INSPIRE].

[59] S.M. Filbert et al., in preparation (2023).

[60] L. Ennesser, P. Martini, A. Font-Ribera and I. Pérez-Ràfols, *The impact and mitigation of broad-absorption-line quasars in Lyman α forest correlations*, Mon. Not. Roy. Astron. Soc. 511 (2022) 3514 [arXiv:2111.09439] [INSPIRE].

[61] A. Font-Ribera et al., *The large-scale cross-correlation of Damped Lyman Alpha Systems with the Lyman Alpha Forest: First Measurements from BOSS*, JCAP 11 (2012) 059 [arXiv:1209.4596] [INSPIRE].

[62] I. Pérez-Ràfols 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 [arXiv:1709.00889] [INSPIRE].

[63] I. Pérez-Ràfols et al., *The cross-correlation of galaxies in absorption with the Lyman α forest*, Mon. Not. Roy. Astron. Soc. 524 (2023) 1464 [arXiv:2210.02973] [INSPIRE].

[64] 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.

[65] M.-F. Ho, S. Bird and R. Garnett, *Damped Lyman-α absorbers from Sloan digital sky survey DR16Q with Gaussian processes*, Mon. Not. Roy. Astron. Soc. 507 (2021) 704.

[66] B. Wang et al., *Deep Learning of Dark Energy Spectroscopic Instrument Mock Spectra to Find Damped Lyα Systems*, Astrophys. J. Supp. 259 (2022) 28 [arXiv:2201.00827] [INSPIRE].

[67] J. Zou et al., in preparation (2023).

[68] A. Font-Ribera, P. McDonald and A. Slosar, *How to estimate the 3D power spectrum of the Lyman-α forest*, JCAP 01 (2018) 003 [arXiv:1710.11036] [INSPIRE].

[69] D.W. Hogg, *Distance measures in cosmology*, astro-ph/9905116 [INSPIRE].

[70] H. du Mas des Bourboux et al., *The Extended Baryon Oscillation Spectroscopic Survey: Measuring the Cross-correlation between the MgII Flux Transmission Field and Quasars and Galaxies at z = 0.59*, Astrophys. J. 878 (2019) 47 [arXiv:1901.01950] [INSPIRE].

[71] K.M. Górski et al., *HEALPix — A Framework for high resolution discretization, and fast analysis of data distributed on the sphere*, Astrophys. J. 622 (2005) 759 [astro-ph/0409513] [INSPIRE].

[72] A. Zonca et al., *healpy: equal area pixelization and spherical harmonics transforms for data on the sphere in Python*, J. Open Source Softw. 4 (2019) 1298 [INSPIRE].

[73] D.J. Eisenstein, H.-J. Seo and M.J. White, *On the Robustness of the Acoustic Scale in the Low-Redshift Clustering of Matter*, Astrophys. J. 664 (2007) 660 [astro-ph/0604361] [INSPIRE].

[74] W.J. Percival and M. White, *Testing cosmological structure formation using redshift-space distortions*, Mon. Not. Roy. Astron. Soc. 393 (2009) 297 [arXiv:0808.0003] [INSPIRE].

[75] S. Youles et al., *The effect of quasar redshift errors on Lyman-α forest correlation functions*, Mon. Not. Roy. Astron. Soc. 516 (2022) 421 [arXiv:2205.06648] [INSPIRE].

[76] G.C. Rudie, C.C. Steidel, A.E. Shapley and M. Pettini, *The Column Density Distribution and Continuum Opacity of the Intergalactic and Circumgalactic Medium at Redshift $\langle z \rangle \ge 2.4$*, Astrophys. J. 769 (2013) 146 [arXiv:1304.6719] [INSPIRE].
[77] K.K. Rogers et al., *Correlations in the three-dimensional Lyman-alpha forest contaminated by high column density absorbers*, Mon. Not. Roy. Astron. Soc. **476** (2018) 3716 [arXiv:1711.06275] [INSPIRE].

[78] A. Cuceu et al., *The Alcock-Paczyński effect from Lyman-α forest correlations: analysis validation with synthetic data*, Mon. Not. Roy. Astron. Soc. **523** (2023) 3773 [arXiv:2209.12931] [INSPIRE].

[79] A. Cuceu, A. Font-Ribera and B. Joachimi, *Bayesian methods for fitting Baryon Acoustic Oscillations in the Lyman-α forest*, JCAP **07** (2020) 035 [arXiv:2004.02761] [INSPIRE].

[80] L. Yang et al., *Metal Lines Associated with the Lyα Forest from eBOSS Data*, Astrophys. J. **935** (2022) 121 [arXiv:2206.11385] [INSPIRE].

[81] M. Blomqvist et al., *The triply-ionized carbon forest from eBOSS: cosmological correlations with quasars in SDSS-IV DR14*, JCAP **05** (2018) 029 [arXiv:1801.01852] [INSPIRE].

[82] S. Gontcho A Gontcho et al., *Quasar — CIV forest cross-correlation with SDSS DR12*, Mon. Not. Roy. Astron. Soc. **480** (2018) 610 [arXiv:1712.09886] [INSPIRE].

[83] A. Lewis, *GetDist: a Python package for analysing Monte Carlo samples*, arXiv:1910.13970 [INSPIRE].

[84] W.J. Handley, M.P. Hobson and A.N. Lasenby, *PolyChord: nested sampling for cosmology*, Mon. Not. Roy. Astron. Soc. **450** (2015) L61 [arXiv:1502.01856] [INSPIRE].

[85] B. Hadzhiyska et al., *Planting a Lyman alpha forest on AbacusSummit*, Mon. Not. Roy. Astron. Soc. **524** (2023) 1008 [arXiv:2305.08899] [INSPIRE].