Optimal 1D Ly 𝛼 Forest Power Spectrum Estimation – III. DESI early data

The one-dimensional power spectrum 𝑃 1D of the Ly 𝛼 forest provides important information about cosmological and astrophysical parameters, including constraints on warm dark matter models, the sum of the masses of the three neutrino species, and the thermal state of the intergalactic medium. We present the first measurement of 𝑃 1D with the quadratic maximum likelihood estimator (QMLE) from the Dark Energy Spectroscopic Instrument (DESI) survey early data sample. This early sample of 54 600 quasars is already comparable in size to the largest previous studies, and we conduct a thorough investigation of numerous instrumental and analysis systematic errors to evaluate their impact on DESI data with QMLE. We demonstrate the excellent performance of the spectroscopic pipeline noise estimation and the impressive accuracy of the spectrograph resolution matrix with two-dimensional image simulations of raw DESI images that we processed with the DESI spectroscopic pipeline. We also study metal line contamination and noise calibration systematics with quasar spectra on the red side of the Ly 𝛼 emission line. In a companion paper, we present a similar analysis based on the Fast Fourier Transform estimate of the power spectrum. We conclude with a comparison of these two approaches and implications for the upcoming DESI Year 1 analysis.


INTRODUCTION
Neutral hydrogen gas between us and distant quasars forms absorption lines at wavelengths shorter than the Ly emission line in the quasar spectrum through absorption and scattering.These absorption lines are collectively called the Ly forest; and they trace the underlying matter distribution in the intergalactic medium (IGM) and the circumgalactic medium (CGM).The Ly forest is consequently a powerful tool to map vast volumes at redshifts 2 ≲  ≲ 5 and probing scales from hundreds of Mpc to below 1 Mpc.Gunn & Peterson (1965) first estimated the density of neutral hydrogen in the IGM.They realised that the measurement of some continuum flux of 3C 9 below the Ly emission line by Schmidt (1965) implied the IGM was mostly ionized.Later work by Lynds (1971) showed that the IGM absorption was in the form of discrete features.In the 1990s, work by many investigators (Bi et al. 1992;Cen et al. 1994;Zhang et al. 1995;Bi & Davidsen 1997) clearly established that this Ly forest originates from smooth IGM fluctuations.Based on this smooth density fluctuations picture, the 1D power spectrum ( 1D ) has emerged as an important quantity to measure in high-resolution, high-signal-to-noise (SNR) spectra (Croft et al. 1998;Iršič et al. 2017;Walther et al. 2017;Karaçaylı et al. 2022), as well as medium-resolution, medium-SNR spectra (McDonald et al. 2006;Palanque-Delabrouille et al. 2013;Chabanier et al. 2019). 1D is valuable because it is sensitive to smaller scales than are accessible in high-redshift galaxy surveys, and consequently to particular physical quantities.Applications of the Ly  1D include investigations of the thermal state of the IGM (Boera et al. 2019;Walther et al. 2019;Villasenor et al. 2022), inference of the primordial power spectrum (Viel et al. 2004), constraints on the sum of neutrino masses (Croft et al. 1999;Palanque-Delabrouille et al. 2015;Yèche et al. 2017), and explorations of the nature of dark matter (Narayanan et al. 2000;Seljak et al. 2006;Wang et al. 2013;Iršič et al. 2017b), with warm dark matter receiving particular attention (Boyarsky et al. 2009;Viel et al. 2013;Baur et al. 2016;Iršič et al. 2017a;Villasenor et al. 2023).
Even though  1D is a summary statistic for cosmological analysis, it is very sensitive to several sources of systematic errors.The fiveyear data from the Dark Energy Spectroscopic Instrument (DESI, Levi et al. 2013) will provide approximately 700 000 Ly quasar spectra with medium resolution ( ≈ 3 000), medium signal-tonoise (SNR) (≈ 2 per Å) (DESI Collaboration et al. 2016a,b), which will constitute a data set that is four times larger than the Extended Baryon Oscillation Spectroscopic Survey (eBOSS, Dawson et al. 2016).DESI will consequently substantially expand the statistical power of Ly forest measurements relative to previous work.To fully exploit this great increase in statistical power requires comprehensive studies of  1D systematics.These include systematics related to the theoretical interpretation (e.g., Lukić et al. 2015;Walther et al. 2021;Chabanier et al. 2023a), instrumental effects, and other spectroscopic extraction details.We address the latter two topics in this paper with early DESI observations.
We analyze two distinct datasets in this paper.The first set of spectra was collected between December 2020 and May 2021 during the DESI Survey Validation (SV, DESI Collaboration et al. 2023a) phase.The purpose of this phase was to perform various tests to verify the pipeline for target selection, spectral extraction, classifier performance and clustering analysis.The spectra collected during this period will be publicly available as early data release (DESI Collaboration et al. 2023b).The second set was obtained during the first two months of the DESI main survey, which began in May 2021.Together, these data span a wide range of signal-to-noise ratios (SNR).We use them to measure  1D and characterize the noise, flux calibration, and spectrograph resolution calculated by the DESI spectroscopic pipeline.
The two main methods to estimate  1D are the maximum likelihood estimator and the Fast Fourier transform (FFT).The maximum likelihood estimator is typically considered to be statistically optimal, although it is slower than FFT-based algorithms.The maximum likelihood estimator can be implemented in two different ways.A direct implementation finds the maximum likelihood solution by sampling the likelihood with respect to  1D (Palanque-Delabrouille et al. 2013).This implementation has slower convergence properties and is more prone to numerical instabilities.The second implementation takes advantage of the Newton-Raphson method and achieves a faster and more stable performance.We call this estimator the quadratic maximum likelihood estimator (QMLE;McDonald et al. 2006;Font-Ribera et al. 2018;Karaçaylı et al. 2020) and the application of QMLE to DESI data is the main focus of this paper.In a companion paper by Ravoux et al. (2023) we present the application of the FFT-based estimator to early DESI data.That paper adapts the FFT approach previously used for eBOSS (Chabanier et al. 2019).
A major virtue of QMLE is that it is robust against challenges such as strong sky emission lines, high-column density (HCD) systems, and bad CCD pixels.Pixels affected by these features must be masked  (TargetID 39627871806818826).The Ly forest is defined to be the spectral region between a quasar's Ly and Ly  emission lines.Absorption features redward of the quasar's Ly emission line may be due to metal systems.The regions from Ly to Si iv and from Si iv to C iv are called the "side bands" (SB).We call the Ly -Si iv region SB 1 and the Si iv-C iv region SB 2, and use these regions to quantify metal contamination, noise and flux calibrations.
to avoid contamination from unrelated physical effects and imperfections in instrumentation.This masking introduces a bias that must be corrected in FFT estimates; and these corrections in turn introduce uncertainties to the measurement (Chabanier et al. 2019).A major advantage of QMLE is that it can handle masked, uneven spectra without further corrections by construction.Relatedly, QMLE is capable of weighting individual pixels by the inverse pipeline noise, and hence diminishes the impact of variations in instrument noise and other noisy spectral regions such as certain sky lines.In addition, the QMLE implementation of Karaçaylı et al. (2020) interpolates pixel pairs into two redshift bins to account for the redshift evolution within the Ly forest.These properties, among others discussed later in the text, make the QMLE an excellent tool for DESI Ly  1D estimation.
For a medium-resolution, medium-SNR survey such as DESI, the potential systematics due to the pipeline noise estimation and the spectrograph resolution require the most attention.Previous experiments suffered from spectroscopic pipeline noise miscalibration levels of 15%, which necessitated separate calculations and recalibrations of the pipeline noise (McDonald et al. 2006;Palanque-Delabrouille et al. 2013).DESI was meticulously designed to abate such miscalibrations (DESI Collaboration et al. 2022;Guy et al. 2023).Yet even though the pipeline is significantly improved, the statistical power of even the early data demands ever-stringent precision.Another consideration is that spectral extraction for DESI is based on the spectro-perfectionism algorithm, which can handle arbitrarily complicated (i.e.not solely separable) 2D point-spread functions (PSF) (Bolton & Schlegel 2010).This extraction preserves the full native resolution of the 2D spectrograph without degradation in the 1D spectrum and yields an independent resolution matrix for each 1D spectrum that is based on the spectrograph resolution and the noise in each spectrum (Guy et al. 2023).QMLE can naturally incorporate this novel resolution matrix, and in this paper, we validate the spectro-perfectionism and its synergy with the QMLE by simulating CCD images and extracting spectra with the DESI spectroscopic pipeline.
The outline of this paper is as follows.First, we describe the DESI survey, target selection, the creation of quasar catalogues, the identification of damped Ly (DLA) systems and broad absorption lines (BAL), and the properties of the early spectra in Section 2. We outline the continuum fitting algorithm and detail the QMLE and various updates in Section 3. Synthetic spectra are central in our validation to make robust statistical claims.In Section 4, we validate the continuum fitting algorithm, DLA masking and damping wing corrections with extensive sets of 1D mock spectra, and validate the resolution matrix derived by the pipeline CCD image simulations that we analyze with the same spectroscopic pipeline that we use with real DESI observations.We perform various tests for systematics and present our  1D measurement from data in Section 5. Finally, we compare DESI  1D measurements from the QMLE and FFT estimators to each other and to eBOSS in Section 6.We summarize our results in Section 7. As noted before, a companion paper by Ravoux et al. (2023) presents the FFT-based results.

DATA
The DESI collaboration began a five-year survey of 40 million galaxies and quasars in May 2021.The main goal of this survey is to measure distances with the Baryon Acoustic Oscillation (BAO) method from the local universe to beyond  > 3.5 and use these data to explore the nature of dark energy.DESI will also employ Redshift Space Distortions to measure the growth of cosmic structures and test potential modifications to general relativity, measure the sum of neutrino masses, and investigate primordial density fluctuations from the inflationary epoch.The collaboration is conducting this survey with a new, high-throughput, fiber-fed spectrograph on the 4 m Mayall telescope that can obtain 5000 spectra in each observation (DESI Collaboration et al. 2016b;Silber et al. 2023).The light from each fiber is directed into one of ten, identical, bench-mounted spectrographs that record the light from 360 − 980 nm in three wavelength channels.The blue channel is optimized for Ly forest studies and extends from 360 − 593 nm with a resolution that ranges from 2000 − 3500.These spectrographs are in a climate-controlled enclosure that provides very stable calibrations and minimizes systematic errors due to instrumental effects.The instrumentation is described in detail in DESI Collaboration et al. (2022) and the spectroscopic pipeline in Guy et al. (2023).DESI targets were selected with , ,  photometry from the Legacy Imaging Surveys (Dey et al. 2019) and 1, 2 photometry from the Wide-field Infrared Explorer (WISE, Wright et al. 2010).The target selection process is described in detail in Myers et al. (2023).The targets include quasars at 0.9 <  < 2.1 that are used to trace large-scale structure and at  > 2.1 that are used to trace the matter distribution with the Ly forest (Yèche et al. 2020).The collaboration refined the target selection algorithms during the Survey Validation (SV, DESI Collaboration et al. 2023a) period in early 2022 with a significant visual inspection effort (Alexander et al. 2023).The final quasar target selection is based on a random forest algorithm and selects quasars in the magnitude range 16.5 <  < 23 (Chaussidon et al. 2023).We use the One-Percent Survey (SV3) spectra that are part of early data release (EDR, DESI Collaboration et al. 2023b), and further include two months of main survey (DESI-M2) to increase the statistical precision in our analysis.We call this combined data set EDR+.The Target Selection Validation (SV1) spectra are the deepest observations in EDR, but their pipeline noise estimates differ from the other two data sets (Ravoux et al. 2023).Therefore, we rely on these spectra only for the DLA identification and not for  1D estimation since the pipeline noise estimates do not affect DLA identification as they affect  1D .Figure 1 shows a quasar at  = 2.94 from this DESI early data.
DESI employs three classification algorithms to identify quasars.Most targets are correctly classified with Redrock1 (Bailey et al. 2023, in preparation).This algorithm performs a  2 analysis for a range of spectral templates as a function of redshift and identifies the best redshift and spectral template for each target.Our visual inspection process demonstrated that Redrock misses some quasars, so we employ QuasarNET (Busca & Balland 2018;Farr et al. 2020) and an Mg ii afterburner (Napolitano et al. 2023) to help identify additional quasars.QuasarNET is a machine learning algorithm that uses convolutional neural networks for classification and the Mg ii afterburner searches for broad Mg ii emission at the Redrock redshift in the spectral of quasar targets classified as galaxies.Chaussidon et al. (2023) describe this process in more detail.We limit ourselves to objects that are targeted as quasars in the afterburner catalogue.

Quasars with broad absorption lines
Broad absorption line (BAL) features are present in about 15% of all quasar spectra and can contaminate the Ly forest as well as impact quasar redshift errors and classifications.The vast majority of BAL quasars exhibit blueshifted absorption associated with the C iv emission feature and the BAL identification algorithm searches this region in every quasar spectrum where this spectral region is visible (1.57<  < 5).This algorithm is similar to the one presented by Guo & Martini (2019), except that it does not use the Convolutional Neural Network (CNN) classifier.Filbert et al. (2023) describe the BAL identification and characterization for the early DESI quasar catalogues in detail, including the catalogue completeness and purity, and the impact of BAL features on redshift errors (García et al. 2023).We use the measured velocity range of the BAL features associated with C iv to mask this ion and also mask the wavelengths that correspond to the same velocities associated with the S iv, P v, C iii, Ly, N v, and Si iv.All of these features may be present in BALs (Mas-Ribas & Mauland 2019), and all but S iv can contaminate the Ly forest (e.g., Ennesser et al. 2022).

Damped Ly𝛼 systems
DLAs are identified using both Convolutional Neural Network (CNN) and Gaussian process (GP) finders, then their results are combined into a concordance catalogue while adopting GP results over CNN if both detect the same DLA (Ho et al. 2021;Wang et al. 2022).We pick unique DLA identifications while combining three separate DLA catalogues for SV1, SV3 and DESI-M2 since the same quasars and DLAs can be present in different catalogues.If two DLAs are within a threshold redshift separation Δ  that corresponds to a DLA's observed redshift size, we pick the highest confidence identification, where We select systems based on average signal-to-noise ratio SNR between 1420-1480 Å in quasar's rest frame.For DLAs that are identified by CNN, we keep them in the catalogue if the host quasar spectrum has SNR > 3, but remove them from the catalogue if the confidence level is less than 0.3 in quasars with SNR < 3. We keep all systems GP identifies.There are 41 946 DLAs with  HI > 20.3 in the combined catalogue.Sub-DLA detections contain many false positives, so we do not mask them.Our selection criteria and duplicate removal reduce this number to 30 131.Introducing a minor confidence threshold of 0.2 for high SNR and 0.9 for GP systems removes 567 DLAs.We believe masking possible DLAs in this small sample is more valuable than missing them.We note that not all DLA sightlines end up in our final sample since some host quasars are left out due to quality cuts.

Redshift distribution
Figure 2 shows the quasar redshift distribution of our sample on the top panel.The quasar distribution  qso () drops off rapidly at higher redshift, as expected from the selection function.There are 67 241 quasars in our final sample.On the bottom panel, we show the SNR distribution in the forest as a function of redshift with bin size Δ = 0.1.We define SNR based on the propagated error (), where SNR=1/() and the propagated error () on the weighted mean as follows: where  −1  =  2 LSS, +  2 pipe, and summation is done over all pixels that fall into the redshift bin.This quantity is equivalent to pixel SNR values after coadding all quasar spectra in the forest region into coarse Δ = 0.1 pixels.Large-scale structure variance  2 LSS is calculated during the continuum fitting process as described in Section 3.Even though we keep BAL quasars while continuum fitting, we remove them from final  1D estimates.Removing these BAL quasars leaves us with 54 600 spectra and reduces our SNR as shown in blue in Figure 2.
We explain the details regarding DLA and BAL masking in Section 5.

Continuum fitting
The continuum fitting algorithm we use was developed over the last few years and has been applied to both 3D analyses (Bautista et al. 2017;du Mas des Bourboux et al. 2019, 2020) and  1D measurements (Chabanier et al. 2019).This algorithm is part of the software Package for Igm Cosmological-Correlations Analyses (picca) and is publicly available2 .We summarize the algorithm below.
One important aspect of the algorithm is that the definition of the quasar continuum absorbs the mean transmission  () of the IGM.Specifically, we model every quasar "continuum"   ( RF ) by a global mean continuum  ( RF ) and two quasar "diversity" parameters, amplitude   and slope   : where  RF is the wavelength in quasar's rest frame and (1,2) RF are the minimum and maximum wavelengths considered for calculation.We assume that the global mean continuum  ( RF ) does not depend on redshift, and therefore our model only adjusts  (), as well as solves for the   and   parameters for each quasar.In other words, the amplitude and slope parameters do not only fit for intrinsic quasar diversity such as brightness, but also for the IGM mean transmission.Given these definitions, transmitted flux fluctuations are given by where  = (1 +   ) RF is the observed wavelength and   () is the observed flux.The effect of spectrograph resolution has been ignored for simplicity as noted in Slosar et al. (2013), since the affected scales are small for 3D analysis.The features in the continuum are also wider than the spectrograph resolution, so this assumption should also hold for  1D .
Our continuum fitting procedure calculates   and   for each quasar, and three global functions: the mean quasar continuum  ( RF ), the large-scale Ly fluctuations  2 LSS (), and the pipeline noise correction term ().We do not assume a functional form for these three functions; instead, we construct linear interpolations based on binned estimates.Specifically,  ( RF ) is calculated between rest-frame wavelengths (1)

RF and 𝜆
(2) RF in bins of size Δ RF .The other two parameters () and  2 LSS () are calculated in the observed frame in  obs bins linearly spaced between  (1) and  (2) .These binning parameters are tuned for each analysis depending on the available statistics.Before we start our fitting process, we coadd the three spectrograph arms using the pipeline inverse variance as weights.Our fitting procedure is iterative.Each iteration  is as follows: • Fit each spectrum for   and   while keeping other parameters fixed.
• Fit for variance parameters  and  2 LSS (defined below) for each bin.
For each quasar, we find the   and   values that minimize the following cost function while keeping all other parameters fixed: where the summation  is over all pixels in the forest region and   is the observed wavelength.The major complication comes from  2 ,  , which must take into account the intrinsic large-scale Ly fluctuations  2 LSS : After every quasar is fit, we stack all continua in the rest frame and update the global mean continuum .As described above, parameters  and  2 LSS are calculated at discrete wavelength bins.For each bin, we rebin the  values with respect to the pipeline noise estimates  pipe and calculate the scatter in these  pipe bins to measure the  2  −  2 pipe relation from the data.Lastly, we fit equation 7 to this relation to find  and  2 LSS values for every wavelength bin.

Quadratic estimator
We measure  1D using the quadratic maximum likelihood estimator (QMLE), which was extensively studied in the 90s in the context of cosmic microwave background radiation, galaxy surveys and weak lensing (Hamilton 1997;Tegmark et al. 1997Tegmark et al. , 1998;;Seljak 1998), and later also applied to the Ly forest (McDonald et al. 2006;Karaçaylı et al. 2020Karaçaylı et al. , 2022)).The QMLE works in real space (instead of Fourier space) to estimate the power spectrum, and therefore allows weighting by the pipeline noise, accounts for intrinsic Ly large-scale structure correlations, and most importantly is not biased by gaps in the spectra.We refer the reader to Karaçaylı et al. (2020) and Karaçaylı et al. (2022) for our development process and application to high-resolution spectra.In this section, we provide a short summary of QMLE and then describe the important steps for the resolution matrix and shifting Nyquist frequency implementations.Details regarding the continuum error marginalization are in Appendix A and signal-noise coupling correction is in Appendix B.
One motivation for the development of QMLE is that the power spectrum is typically estimated on discrete wavenumbers  as band powers, since it cannot be estimated continuously on , and this discretization inevitably averages the underlying power over these bands.Our QMLE implementation alleviates this effect by estimating deviations from a fiducial power spectrum such that (, ) =  fid (, ) + ,  () (, ) () , where we adopt tophat  bands with   as bin edges and linear interpolation for  bins with   as bin centres.This fiducial power spectrum further improves the weighting by including large-scale Ly correlations, does not have to exactly match the true power spectrum, and can be approximated in an unbiased way if no safe guess is available (as shown by Karaçaylı et al. 2020).We use the following fitting function: where  0 = 0.009 s km −1 and  0 = 3.0, and stress that this is sufficient for a baseline estimate, which in turn can be used to weight pixels, but does not capture all of the scientific information in  1D .Given a collection of pixels representing normalized flux fluctuations   , the quadratic estimator is formulated as follows: where  is the iteration number and The covariance matrix C ≡ ⟨     ⟩ is the sum of signal and noise, C = S fid +  Q    + N, Q  = C/  and the estimate of the Fisher matrix is The covariance matrices on the right hand side of equation ( 9) are computed using parameters from the previous iteration  () .Assuming different quasar spectra are uncorrelated, the Fisher matrix   ′ and the expression in parentheses in equation ( 9) can be computed for each quasar, then accumulated, i.e.F =  F  etc.
We convert wavelength to velocity using logarithmic spacing, following the convention in cosmology: where  Ly = 1215.67Å.We assume the noise is uncorrelated at different wavelengths, which results in a diagonal noise matrix with   =  2  , where   is the continuum-normalized pipeline noise.

Resolution matrix
In our previous applications of QMLE, we approximated the spectrograph resolution effects by a continuous window function  (,  ′ ) such that the smoothed flux fluctuations   were given by Even though it is a valid and prevalent approximation, this formalism unfortunately fails to capture wavelength dependent resolution of the spectrographs.However, for DESI the spectral extraction is built on the improved spectro-perfectionism algorithm (Bolton & Schlegel 2010;Guy et al. 2023).Spectro-perfectionism produces a resolution matrix R associated with each spectrum that is based on the spectrograph resolution as well as the noise properties of each spectrum, and captures the wavelength dependent resolution on the same discrete wavelength bins as the spectrum.The observed signal becomes a matrix-vector multiplication: This redefinition is natural to incorporate into the QMLE formalism.
We achieve this by replacing the integral equations for the signal S and derivative matrices Q by the following expressions: where the subscript  denotes the smoothed matrices, and matrices without a subscript are evaluated as integrals but now without a resolution window function. where The derivative matrix for redshift bin  and wavenumber bin  is where   () is the interpolation kernel.This is 1 when  =   and 0 when  =  ±1 .However, there are more subtleties regarding the resolution matrix.First, these matrix multiplications require that all pixels are present, so we mark masked pixels with large noise estimates instead of eliminating them from the spectrum.Second, the resolution matrix does not capture the resolution outside the spectral range (by construction).This is a potential problem at the largest scales, so we implement an option in QMLE that pads the resolution matrix by mirroring its columns at the edges.Third, both synthetic spectra and the actual DESI pipeline produce this matrix on the same wavelength grid as the spectrum with the same spacing.This is natural in the spectro-perfectionism formalism in data, and we test its accuracy in Section 4.2; however, it yields an undercorrection at small scales in the mock analysis.Our solution to this problem in mocks is to oversample every row of the resolution matrix (Appendix D, Guy et al. 2023).One could model the resolution matrix at each row (i.e.wavelength) as a convolution of Gaussian and top-hat window functions, and fit for one or two free parameters for this model.One then evaluates each row of the oversampled resolution matrix using the best-fitting parameters at smaller wavelength steps.However, spectro-perfectionist resolution matrix carries negative elements and evidently does not follow this simple description.Therefore, achieving a stable oversampling method requires a nuanced procedure.We decided to use an unassuming description by interpolating the intermediate values.To correctly capture the rapid change in resolution matrix elements, we interpolate using their natural logarithms with a cubic spline.To obtain a valid natural logarithm, we shift every element to a small positive value in each row.This small positive value is the smallest absolute value in that row (using an arbitrary number breaks down in subsequent steps).We then apply a cubic spline to the natural logarithm of these elements, oversample at a desired factor (usually three), and finally trace back these changes to obtain the new resolution matrix.

Shifting Nyquist frequency on a linear wavelength grid
Another update to QMLE concerns the Fisher matrix and DESI's wavelength binning.The DESI pipeline extracts spectra on a linear wavelength grid of Δ = 0.8 Å, which results in an increasing Nyquist frequency with wavelength in velocity space  Ny = /, where  = Δ/.In other words, we can measure higher  modes at higher redshifts.However, forcing the code to measure the same  bins at lower redshifts results in numerically unstable Fisher matrix elements that could contaminate all scales when inverted.Hence, these modes should be removed from the analysis.
We decide each spectral segment's Nyquist frequency using their median , then set  >  Ny /2 modes in the Fisher matrix and the power spectrum to zero.Since this procedure results in a "singular" matrix, we replace zeros on the diagonal with one while inverting.Note this replacement does not contaminate lower  modes, because it constitutes a block diagonal matrix.This process stabilises the Fisher matrix.

Nominal estimator settings
Throughout this paper, we use 20 linear bins with Δ lin = 0.5 × 10 −3 s km −1 and 13 log-linear bins with Δ log = 0.05.We use redshift bins of size Δ = 0.2 from  = 2.0 to  = 3.8 included.To reduce computation time and help continuum marginalisation, we split the spectra into two segments if they have more than 500 pixels, and we ignore segments having less than 20 remaining pixels.We interpolate the signal and derivative matrices using 3601 points in velocity with 10 km s −1 spacing and 400 points in redshift.

Software
Our quadratic estimator 3 is written in c++.It depends on cblas and lapacke routines for matrix/vector operations, GSL4 for certain scientific calculations (Galassi et al. 2021), FFTW35 for deconvolution when needed (Frigo & Johnson 2005); and uses the Message Passing Interface (MPI) standard6,7,8 to parallelize tasks.The DESI spectra are organized using HEALPix (Górski et al. 2005) scheme on the sky.We use the following commonly-used software in python analysis: astropy9 a community-developed core python package for Astronomy (Astropy Collaboration et al. 2013, 2018, 2022), numpy10 an open source project aiming to enable numerical computing with python (Harris et al. 2020), scipy 11 an open source project with algorithms for scientific computing, healpy to interface with HEALPix in python (Zonca et al. 2019), numba 12 an open source just-in-time (JIT) compiler that translates a subset of python and numpy code into fast machine code, mpi4py13 which provides python bindings for the MPI standard (Dalcin & Fang 2021).Finally, we make plots using matplotlib14 a comprehensive library for creating static, animated, and interactive visualizations in python (Hunter 2007).

VALIDATION
Synthetic data are crucial to verify that the measurements are unbiased, and the errors are correctly captured.Our mock generation procedure consists of the generation of transmission files with forest fluctuations, diverse quasar spectra, and simulation of the DESI instrument.The lognormal mock transmission files are generated using the procedure in Karaçaylı et al. (2020).We generate them on a linear wavelength grid of 0.2 Å spacing without any resolution and noise effects.
We develop two methods to simulate and validate the DESI analysis pipeline.The first set of mocks is produced using quickquasars, which is part of the desisim package15 and uses specsim16 (Kirkby et al. 2021) for quick simulations of fiber spectrograph response (see Herrera-Alcantar et al. (2023) for a detailed description of quickquasars mocks).This program generates random quasar continua, simulates sky and instrumental noise, and incorporates wavelengthdependent camera resolution, but does not validate the computationally expensive spectral extraction.Hence, we cannot validate the spectro-perfectionist resolution matrix with these mocks.In order to apply and validate the spectro-perfectionism algorithm in the Ly forest, we create a second set of mocks called "CCD image simulations" that project mock quasar spectra onto two-dimensional images that simulate DESI raw data at the CCD pixel level with the desisim  comparison for different  cuts and continuum marginalisation polynomials on mocks.We find  2  increases for all settings as we include higher  values, which is unfortunate but expected since our correction to the quickquasars resolution matrix is not exact.The true continuum analysis results (blue triangles) stay within 1.5 of  2  = 1.Lower rows correspond to larger small-scale confidence regions.From left-most column to the right, we remove large-scale modes.When continuum errors are not marginalised (orange squares), throwing out these large-scale modes brings  2  down to 1 within error bars.We also find that first (green circles) and second (red triangle) order marginalizations remove the contamination from continuum errors at all scales.
package.These CCD image simulations are then processed in a similar manner to actual data with the algorithms that comprise the DESI spectroscopic reduction pipeline (Guy et al. 2023).This approach is more computationally expensive than one-dimensional mocks, so we only employ it on a smaller number of mock spectra.

Quickquasars mocks
For these mocks, the quasar diversity, DESI instrument and the sky are simulated through a program called quickquasars in the desisim package.This program randomly generates quasar continua from a broken power law with emission lines, convolves with the wavelength-dependent camera resolution for each arm, adds noise for a given exposure time and observation program, and finally resamples onto the output DESI wavelength grid of Δ DESI = 0.8 Å per pixel.We smooth out the source contribution to noise with a Gaussian kernel of  = 10 Å to imitate the DESI pipeline (Guy et al. 2023).
All unique targets that are identified as quasars are simulated in our mocks.However, in real data analysis, we remove certain surveys, programs and low SNR targets.We generate the transmission files with the exact redshift distribution of DESI quasars in our sample and assume a constant 4000  exposure time for all spectra.
As extensive and realistic as quickquasars is, it does not fully reproduce the spectral extraction pipeline output since it does not generate 2D CCD images.As an important consequence, the output resolution matrix does not follow the spectro-perfectionism formalism and instead it is approximated as a box-car average over rows and columns of the finely sampled camera resolution matrix.Unfortunately, this approximation is not correct as it smoothes the resolution matrix twice, once over rows and once over columns.To correct that implementation, we deconvolve a top-hat window function and oversample each row of this matrix by a factor of three in the power spectrum estimation.This yields adequately unbiased power spectrum results but is not a precise enough solution to strongly rely on  2 criteria.We also perform CCD image simulations to understand the behaviour of the resolution matrix in data.
As noted, we generate a mock spectrum for each unique target in our sample, which yields 92 780 quasars in our final sample.We define the forest to be between 1050-1180 Å in the quasar rest frame, use the analytically calculated true power spectrum as our fiducial and perform a single iteration using the QMLE (Karaçaylı et al. 2020).We define our small-scale confidence range with respect to effective velocity spacing   = Δ DESI /(1 + ) Ly of each redshift bin, where Δ DESI = 0.8 Å.

True continuum, no systematics
We start validating our analysis without any continuum fitting complications or other systematics.We obtain flux transmission fluctu- ations   using the true continuum (which is provided by quickquasars) and true mean transmission.We estimate the mean transmission from pure transmission files, confirm that this estimate is correct using the analytical mean transmission expression (Karaçaylı et al. 2020), and use the analytical expression to remove the measurement noise.
We find that the estimated power spectrum agrees with the true underlying power albeit the problems at small scales due to inaccuracy in the resolution as mentioned above.We calculate the reduced chi-square  2  =  2 /, where the number of degrees of freedom  is equal to the number of (, ) points in the range of concern, and  2 = ( −  true )  C −1 ( −  true ) where  is the measurement,  true is the underlying true power spectrum and C is the covariance matrix from QMLE. Figure 3 shows the reduced  2  from the true continuum analysis results in blue triangles. 2  values increase as we include higher  values (going from the top to the bottom row), which is unfortunate but expected since our correction to the quickquasars resolution matrix is not exact.The   < 0.9 range is firmly validated with  2  ≈ 1.The   < 1 range deteriorates the agreement between power spectra by 1.5, and finally, the agreement breaks down in the   < 1.2 range.

Continuum fitting, no systematics
We now turn to validating our continuum fitting procedure, since we do not have access to the true quasar continua.The quickquasars code generates quasar continua with broken power laws and emission lines, so our continuum fitting model with a single global mean and two diversity parameters is not exact.Therefore, our mock continuum is not tailored toward our fitting model, and the test results we present here also capture model mismatches.
There are 92, 780 quasars in our mock data set.We find that fitting for  LSS is not valid for observed wavelength  > 6000 Å due to the small number of high-redshift quasars with forest data at these wavelengths (only 883), so we limit our continuum fitting region to 3600 − 6000 Å.This sets  Ly = 3.8 as our largest redshift bin.We measure the global mean continuum  ( RF ) in 2.5Å steps.We fix  = 1, and measure  2 LSS in 20 wavelength bins in the observed frame in equation ( 7).We do not apply a SNR cut in order to keep all spectra and perform five iterations.
In Figure 4, we compare the mean continua from the true continuum analysis to the one from continuum fitting.Continuum fitting accentuates peaks and valleys in the mean continuum compared to when the true continuum is known.These deviations are interesting and merit further investigation, but our main objective is to obtain unbiased  1D results.As we discuss below, these deviations do not impede that objective.The bottom panel of Figure 4 shows  2 LSS estimated from the true and continuum fitting analyses.We find fitting the variance leads to the correct  2 LSS values.We note that  2 LSS is not only a function of  1D , but also depends on spectrograph resolution and wavelength spacing.We investigate reduced  2  values for various settings to judge the accuracy of the  1D estimate.In addition to the true continuum results, Figure 3 shows results for no continuum marginalization (orange squares), first order ln  polynomial (green circles), and second order polynomical (red triangles). 2  from no the marginalization case is not visible in the left-most column, but produces reasonable values when large-scale modes are removed in the middle and right columns.This result illustrates the importance of continuum marginalization, especially in that we can retain even the largest scales.We note that this analysis does not account for metals or DLA systematics, which dominate at these scales.Furthermore, we estimate the power spectrum produced by the remaining continuum errors by calculating   ≡  est / true − 1 and running it through QMLE.We do not subtract noise and fiducial terms in this case, but keep everything else the same.We find this continuum error power spectrum is a factor of 10 −5 smaller than the signal at most scales and redshifts as shown in Figure 5.With these results, we consider our continuum fitting and marginalization validated for this work.In future work, we will test our analysis pipeline on multiple (ideally 100) realizations and directly study the  2 distribution.

Masking high-column density systems
We finally tested masking the high-column density (HCD) systems both in continuum fitting and in the  1D estimate.We generate mocks with randomly placed HCDs and build a truth catalogue using their redshifts and column densities.There are 17 273 HCDs with  H i > 19.5 on 15 097 sightlines, which corresponds to 16.2% of all sightlines.
We calculate a DLA transmission profile based on the column density for each system.The damping wings extend to large wavelength separations from the central wavelength, such that an aggressive masking strategy would remove many data points.Therefore, we mask pixels where the model profile is below a transmission threshold and correct the damping wings at larger transmission values based on the same model profile.A higher transmission threshold results in smaller corrections, but eliminates more data points.We tested two thresholds: a nominal threshold with the DLA absorption greater than 20% and a conservative cut of greater than 10%.We find both options yield similar  2  ∼ 1 within error bars.When unmasked, these systems add power to the large scales, as shown in Figure 6.Furthermore, this extra power depends on whether these systems are correlated with the underlying matter field (Mc-Donald et al. 2005), and it has different amplitudes and shapes for different column densities (Rogers et al. 2018).Accurate simulation of these systems embedded in the Ly forest remains challenging, yet they cannot be completely removed from the measurement either.For instance, the catalogue produced by the DLA finder is approximately 90% efficient and pure (Wang et al. 2022).As we discuss in Section 5.3, we report the systematics associated with the DLA finder inefficiency based on a simple scaling of this ratio.However, the effect of uncorrected and undetected DLAs still needs to be modelled and marginalized over in cosmological inferences (Palanque-Delabrouille et al. 2013;Chabanier et al. 2019).

CCD image simulations
We validated the resolution matrix implementation in the DESI pipeline and the QMLE with CCD image simulations.The image simulations draw on extensive infrastructure in the desisim package that was built by the DESI team to develop and validate the spectroscopic data processing pipeline (Guy et al. 2023) in advance of first light.This package produces realistic two-dimensional spectroscopic image simulations that include the bias, readnoise, and gain for each amplifier for each of the 30 CCDs, models the throughput based on engineering data for each channel of each spectrograph, the model PSF and trace of each fiber as a function of wavelength and position on the detector based on the Zemax optical design models, sky emission, and applies noise appropriate for the flux of each object at the desired exposure time.
The typical input to the simulation code is a library of the model spectra for a single DESI observation and a file that describes how these objects are distributed into the fibers.This mapping of targets to fiber positions for one observation corresponds to a single DESI tile.Normally, a tile would include all of the DESI target classes, along with a selection of flux calibration stars and designated 'empty' fibers that are used to measure the night sky spectrum.Since we are just interested in validating the performance with the Ly forest and the time to construct one simulation tile of 30 CCD images is not insignificant, each of our tile designs has only quasars, standard stars, and sky fibers.There are approximately 4500 quasars in each tile, and the quasars are all at 2.6 <  < 3.6 so we observe the entire forest region for each quasar in the blue channel of the spectrographs.We generate ten tiles, or approximately 45,000 total Ly forest spectra, with apparent magnitudes representative of DESI quasar targets and noise that is representative of a single, 1000 s exposure in nominal dark conditions.We also generate a mock set of arc calibration lamp exposures, which are used by the pipeline to measure the PSF and for wavelength calibration.
We processed this simulated dataset with the DESI pipeline, which is described in detail by Guy et al. (2023).Briefly, the pipeline preprocesses the raw images to remove the bias level, fits the arc calibration lamps to measure the nightly PSF, spectral trace, and wavelength calibration, refits the nightly PSF to the PSF of each exposure, extracts the spectra, applies a flat field, subtracts the sky, fits the flux calibration stars, applies the flux calibration, and measures the redshift and identifies the best spectral type for each target.The pipeline fits the PSF with an empirical model that consists of a linear combination of Gauss-Hermite functions.The pipeline model is an excellent but not exact match to true PSF, which is well-approximated by the (non-parametric) PSF predicted by the optical design, and we generated the simulated arc and observation files with the optical model PSF to include any systematic error associated with the pipeline's PSF model fit in our analysis.To isolate the impact of fitting the PSF model from the remainder of the spectral extraction and calibration steps, we also processed the simulations with the correct PSF model.
Figure 7 (top) shows the results for the case where the pipeline starts with the correct PSF model.This panel shows the ratio of the measured power spectrum  to the input power spectrum  true based on the resolution matrix provided "as is" by the pipeline (blue squares) for a subset of the quasars at  = 2.8.The ratio is consistent with unity over the entire range in , with the exception of the largest scales (smallest ), which are in any case not usable due to continuum fitting errors.We also explored oversampling the resolution matrix to better match the input model (see Section 3.2.1 for details), and padding the resolution matrix to remove edge effects.None of these modifications is superior to using the resolution matrix provided by the pipeline.
We also analyzed the performance of the pipeline with an empirical PSF measured from the arc calibration lamps, rather than with the input PSF model as in the previous case.The bottom panel of Figure 7 shows the difference in the 1D power spectrum measurement starting with the arc calibration lamps  arcfit and starting with the correct model PSF  model for seven redshift bins from 2.2 to 3.4.The fractional difference relative to the true power spectrum is at the level of 1-2% up to  = 0.01 s km −1 after which the errors grow exponentially as expected.This behaviour is consistent with 1% precision on the resolution itself.This is our best estimate of the systematic error contribution of the resolution matrix to the measurement.
A potential limitation of these simulations is that real DESI observations have mostly galaxies, rather than quasars.While our simulations have typical numbers of standard stars and sky fibers, they have no galaxies and consequently do not simulate potential cross-talk between galaxy and quasar targets.Guy et al. (2023) carefully studied cross talk between adjacent fibers and found that this is minimal, even for bright calibration lines, so we do not anticipate this will be important for much fainter the continuum and emission lines present in galaxies and quasars.

RESULTS FROM DATA
After our comprehensive validation tests on the synthetic spectra, we now analyze DESI early data.The Ly forest is measured in the 1050 − 1180 Å rest-frame region of each quasar; and the global mean continuum is calculated using Δ RF = 0.8 Å coarse rest-frame binning pixels in this range.We limit our analysis to the observed wavelength range of 3600−6000 Å.We also remove forests with mean SNR less than 0.25 in the forest region in order to minimize possible impurities in the quasar catalogue, where SNR ≡ / pipe17 .The locations of sky lines are naturally down-weighted by the pipeline, but we mask certain particularly strong lines because they are difficult to reliably model18 .DLAs are difficult to simulate accurately and complicate cosmological inference from  1D .We thus mask DLAs at the wavelengths of their strongest absorption  < 0.8 and correct for the damping wings (due to both Ly and Ly  transitions) above this threshold.Furthermore, BAL features can contaminate the forest, and hence these features are also masked.We do not estimate the pipeline noise calibration errors simultaneously (i.e.we fix  = 1 in eqn.7) and only measure  2 LSS in 20 observed frame wavelength bins.As we show below, the noise and flux reported by the pipeline have calibration errors, but due to heavy absorption and correlations between pixels, the Ly forest region is not stable to calibrate for these errors.Instead, we carry out a meticulous study of statistics in the side band regions to calibrate our final reduction.We limit continuum fitting to five iterations where we update the global mean continuum and  2 LSS .We estimate  1D using the following fiducial power parameters in equation ( 8):  = 0.066,  = −2.685, = −0.22, = 3.59,  = −0.16 and  1 = 0.053 s km −1 .We neither oversample nor pad the resolution matrix, and use it as provided by the pipeline.We perform a single iteration to measure the power spectrum.Further iterations mostly refine Fisher matrix estimates (Karaçaylı et al. 2020), which we replace with a regularised bootstrap estimate as described below.Even though we keep BAL quasars with masked features while fitting the continuum, since Ennesser et al. (2022) showed masking BAL lines yields an uncontaminated estimate of the mean continuum, we remove them from the  1D estimation to be conservative in our approach 19 .
Wavelengths larger than the Ly line in the spectrum are free from neutral hydrogen absorption, so they can be used to statistically estimate metal contamination and other systematics (McDonald et al. 2006;Palanque-Delabrouille et al. 2013;Chabanier et al. 2019).The regions between strong emission lines at these wavelengths are called the "side bands" (SB).Table 1 lists the wavelength ranges for the side bands and the number of quasars in each region.Subtracting the SB 1 power spectrum statistically removes all power due to metals with  RF > 1380 Å, but some metal contamination remains.For example, the Si iii-Ly cross-correlation imprints oscillatory features (Mc-Donald et al. 2006;Palanque-Delabrouille et al. 2013).We provide Figure 8. Final Ly forest  1D results.We remove BAL quasars from our sample, mask DLAs and major sky lines, and correct for pipeline noise and flux miscalibrations.Metal power is subtracted using side bands as described in Section 5.2.Error bars are from 200 000 bootstrap realizations of 256 subsamples.Our  1D results are slightly larger than eBOSS measurements (Chabanier et al. 2019), which is most visible at  = 2.2 and 2.4 bins.the details regarding the metal power estimations and other studies on systematics in the following subsections.

Wavelength range [Å] # All quasars
The results after we subtract the metal power are shown in Figure 8.These  1D results are 5 − 15% larger than the eBOSS measurements (Chabanier et al. 2019), which corresponds to 1.5 − 3 tension.This tension is most visible in the  = 2.2 and 2.4 bins, but is present in all redshift bins.We present possible explanations for the origins of this discrepancy in Section 6.Furthermore, accurate noise estimates are crucial to our final  1D results.Figure 9 shows the ratio of the noise power spectrum (equation ( 11)) to the noise-subtracted Ly power spectrum, which ranges between 25-100% and is larger at higher  values.The features in this figure can mostly be attributed to the inverse of  1D , but to be exact, QMLE's noise power spectrum is an inverse covariance weighted average and therefore manifests features based on the fiducial Ly power spectrum, continuum marginalization, and characteristics of pipeline noise.We stress that this noise power strictly comes from the randomness of observed flux values, and is not related to metal absorption, continuum fitting errors, DLA masking, or noise correlations between pixels.Even small miscalibrations can directly propagate to our final  1D estimates.The SB power subtraction can balance out some miscalibration, but not all of it.We find that the SB noise power estimate is 10 − 25% of the Ly region as shown in Figure 10.Fortunately, the pipeline noise estimates are accurate at the percent level and the remaining errors can be corrected by investigating the side bands.

Bootstrap error estimates
The Fisher matrix given by the QMLE assumes Gaussianity, and hence may not be representative of the statistical errors in the data either due to non-linearities in the Ly forest at small scales or other effects in quasar selection, DLA masking etc.For that reason, we calculate the bootstrap covariance matrix for a more reliable error estimate as follows.First, we save QMLE's power spectrum and  Fisher matrix estimates in 256 sub-samples 20 .We then estimate the bootstrap covariance using these sub-samples over 200 000 realisations.As we noted in Karaçaylı et al. (2022), the bootstrap covariance is noisy (especially off-diagonal terms) and needs regularisation.We take advantage of the sparsity pattern of the covariance matrix (Padmanabhan et al. 2016), and regularise the bootstrap covariance as follows: (i) We apply a sparsity pattern on the bootstrap covariance matrix using the Gaussian covariance matrix from QMLE such that  QMLE  We repeat these steps until convergence or for a maximum of 500 iterations.We choose  min = 0.01 for Ly and  min = 0.001 for SB 1, because SB 1 values are more strongly correlated.
Figure 11 shows that the bootstrap error estimates are mostly larger than the Gaussian estimates except between 0.003 s km −1 ≲  ≲ 0.1 s km −1 for  ≳ 3.0.QMLE consequently somewhat underestimates the errors on most scales for most redshifts.

Side bands
We have investigated wavelengths larger than the Ly line in the spectrum to statistically estimate metal contamination and other systematics as they are free from neutral hydrogen absorption (McDonald et al. 2006;Palanque-Delabrouille et al. 2013;Chabanier et al. 2019).As mentioned previously, we use SB 1 to estimate the metal power in the Ly forest.In this section, we provide details for the SB 1 power spectrum measurement and further make use of the SB 2 and SB 3 regions as diagnostics of the metal power and other systematics such as noise calibration.As before, we mask BAL features on all continuum fitting reductions, then ignore these quasars in further analysis.
We first fit the continuum in the SB 1 and SB 2 regions while fixing  = 1.We find the power in SB 1 is larger than in SB 2 as expected, except at  = 2.0, where  SB 2 >  SB 1 for  ≲ 0.003 s km −1 .This likely points to some remaining continuum errors in the side bands.The C iv doublet feature is clearly visible in our estimates.We refer the reader to our companion paper on modelling the doublets in the side band power spectrum (Karaçaylı et al. 2023).
The accuracy of noise calibration and its dependence on SNR can also be studied using these side bands.We divide the spectra into high (SNR > 2) and low (SNR < 2) signal-to-noise ratio samples.This corresponds to approximately a 30/70% split in terms of the number Figure 12.Ratio of metal (SB 1) power spectrum to metals-subtracted Ly power spectrum.The metal power is well below 20% for  ≳ 0.001 s km −1 , but it is potentially a significant source of systematic error at larger scales.
of quasars for both side bands.We find that the low SNR sample sometimes has larger power, as was the case in SDSS (McDonald et al. 2006), but we note that the low SNR sample yields significantly noisier  1D estimates, which hinders strong conclusions, so we further explore this dependence using variance statistics below.Unfortunately, the power difference between the low and high SNR samples is potentially due to the SNR dependence of the pipeline noise estimates.Therefore, we remove the metal power from the Ly forest estimates using our best estimates after we recalibrate the pipeline noise, which is different from what was done in Mc-Donald et al. (2006).Finally, Figure 12 shows the ratio of the SB 1 power spectrum to the metals-subtracted Ly power spectrum.The metal power is potentially a large source of systematic error at  ≲ 0.001 s km −1 , but its effect is well below 20% at higher  values.

Noise calibration error
As we have alluded to before, the pipeline noise estimates can suffer from miscalibrations, which then directly propagate to the final Ly  1D estimates.Fortunately, smooth quasar continua in relatively absorption-free regions (i.e. the side bands) provide near-ideal data to investigate the estimated pipeline variance vs. observed variance in the data.These regions can furthermore be observed in numerous lower redshift quasars, so they provide robust statistical measurements.Our continuum fitting algorithm quantifies this noise calibration error through the  parameter, which is measured by comparing the scatter in   to the reported pipeline  pipe values.The pipeline noise is underestimated for  > 1 and overestimated for  < 1.
First, we fit the continuum on all three side bands while keeping  = 1 fixed.Then we calculate the multipoles ⟨  ⟩, ⟨ 2  ⟩ and ⟨ 4  ⟩ in logarithmic  pipe bins.We find that use of

√︃
(⟨ 4  ⟩ − ⟨ 2  ⟩ 2 )/ as an error estimate yields biased results, so we instead estimate the error on the observed scatter ⟨ 2  ⟩ with the delete-one Jackknife method over sub-samples.Finally, we fit for  and  LSS using equation ( 7), where the continuum is already taken into account in   .
We find that all three side bands yield similar  values to those shown in Figure 13.Our calculated  values fluctuate by a few Figure 13.Pipeline noise correction term  on all three side bands.These regions are relatively absorption free and can be observed in lower redshift quasars, so they provide robust statistics.(Top) All three side bands show the same  trend.We find that the pipeline noise estimates are correct at the percent level.The sharp feature at 4,800 Å occurs at the boundary between CCD amplifiers.(Bottom) Average  over three side bands.The SNR > 2 sample (orange circles) has higher  than lower SNR < 2 sample (blue triangles).This difference is 1.1% on average.We correct the pipeline noise estimates by  of all spectra (black squares) and assign 1.1% as our noise systematic error budget.percent around one, which means the pipeline noise calibration is correct at the percent level.The boundary between CCD amplifiers is responsible for the sharp feature at 4,800 Å.The average noise calibration correction  is shown in black squares in the lower panel.We correct the pipeline noise estimates using this mean  as a function of wavelength and perform a final calibrated continuum fitting, which further includes the flux calibration as discussed below.We again split the data into high (SNR > 2) and low (SNR < 2) signal sub-samples and measure the  parameter in each sample.These are shown with orange circles and blue triangles in the same lower panel of Figure 13.We observe a clear SNR dependence in the noise calibration error, as has been indicated by the power spectrum estimates on these same low and high signal sub-samples.We assign the average difference   = 1.1% to the systematic error budget.
Figure 14.Stacked normalized flux from all three side bands in observed wavelength.We smooth the normalized flux with a 4.8 Å moving box-car average to suppress spurious fluctuations.Residual errors peak at most at 3% at the Balmer and Ca ii H&K doublet lines.We correct the pipeline flux and noise estimates using the average of all three side bands (black line) and perform the final calibrated continuum fitting (which also includes the  correction).This calibration removes significant features from SB 1 (red line), such that the remaining fluctuations are 0.2% on average.
Finally, we calculate the average pipeline inverse variance as a function of flux transmission fluctuations   and wavelength.The pipeline inverse variance overweights   = −1 pixels at  obs = 5200 Å which we confirm does not happen in our mock analysis.We do not observe such a peak in SB 1, so we speculate that it is a bias in noise calibration due to substantial Ly absorption.We confirm that using  LSS in the continuum fitting and  fid in the QMLE removes this feature from the weights.Therefore, we conclude that it does not bias our results.

Flux calibration error
Following du Mas des Bourboux et al. ( 2020), we tested possible flux calibration errors by stacking all normalized quasar spectra in the observed frame using all three side bands.This stacked normalized flux is shown in Figure 14.The residual errors are at most 3%, and we find that the largest errors are at the locations of the Balmer series and Ca ii H&K doublet lines.
We average the stacked normalized flux in all three side bands and smooth with a 4.8 Å moving box-car average.We then divide the flux and noise estimates by this correction factor.This calibration removes significant features from SB 1 (red line), such that the remaining fluctuations are 0.2% on average (smoothed on the same 4.8 Å scale).

Systematic error budget
We identified four sources of systematic errors.First, given the large contribution of the noise power, our power spectrum estimates are susceptible to pipeline noise misestimates.Second, the DLA finder is not perfect and can miss some DLAs or introduce false positives, which would add power to large scales.Third, the spectrograph resolution limits the smallest scale we can measure, and its uncertainties must be propagated to the remaining scales.Finally, our continuum marginalization does not remove all modes of error, so we estimate the systematic error due to the remaining fluctuations.
Our systematic error analysis relies on various scalings of the underlying signal.Directly using the power spectrum estimates introduces spurious fluctuations, so we instead use a smooth power spectrum  smooth (, ) that we calculate as follows: We apply Smooth-BivariateSpline to the logarithms of 1 + ,  and  1D with corresponding statistical weights  stat / 1D and a smoothing factor that is five times the number of data points.This function in scipy is based on the algorithm presented in Dierckx (1981) and Dierckx (1995).
• Noise: As we discussed above, we use the average difference of  values in the side bands between high and low signal sub-samples to quantify the systematic error budget shown in Figure 13.
where   = 1.1% as discussed previously.In this case, we do not smooth the noise power spectrum   since it does not suffer from random fluctuations.Possible origins of this error include correlated CCD readout noise and errors in sky estimates per fiber.We revisit this point at the end of Section 6.
• Resolution: Based on our analysis in Section 4.2, we assign 1% precision to the pipeline resolution estimates that is approximately consistent with the redshift average values in Figure 7.
where   ≡ Δ DESI /(1 + ) Ly as defined in Section 4.1.We do not directly use the values in Figure 7 as it could introduce double counting of noise calibration errors as spectro-perfectionism couples noise and resolution.
• Incomplete DLA removal: Wang et al. (2022) reports over 90% efficiency and purity for the DLA finder for a range of SNR and column density ranges.To be conservative, we pick the worst performance number at SNR= 1 − 2 and assign an average 15% inefficiency ratio and multiply the smooth power estimate with the redshift average of the absolute ratio of  wDLA / true − 1 in Figure 6.
• Continuum errors: We use the absolute ratio of  cont / true in Figure 5, which is on the order of 10 −5 .This error also scales with the smooth power.
Note that continuum errors themselves are larger than 10 −5 as shown in Figure 4.However, our marginalization prevents the bulk of these errors from contaminating the  1D measurements and the remaining continuum errors are fortunately extremely small.Furthermore, the modes that are most affected are not in our conservative range.
Figure 15 shows the systematic errors relative to the statistical errors.DLA systematics heavily affect the large scales, whereas the resolution systematics become relevant at  ≳ 0.01 s km −1 (small scales).The strength of noise systematics decreases with redshift as a consequence of decreasing number of quasars.Change from linear to log-linear binning at  = 0.01 s km −1 increases the bin size and causes the jump at this wavenumber for the noise systematics ratio.

DISCUSSION
In this section, we first provide a qualitative comparison of our results to FFT estimates on the same DESI EDR+ data, and compare both to eBOSS measurements.We then discuss the effects of each systematics and their contribution to the covariance matrix.A full cosmological comparison with respect to physical parameters will be the subject of future analysis with DESI Year 1 data, and will require significant development to generate hydrodynamic simulations and train emulators (Pedersen et al. 2021;Cabayol-Garcia et al. 2023;Chabanier et al. 2023b).We instead provide a comparison using a simplified version of equation (8).

Consistency between methods and literature
Figure 16 compares  1D from DESI EDR+ using FFT and QMLE methods.Both methods agree with each other within error bars until half the Nyquist frequency.We note that the FFT sample is limited to SNR>1, whereas the QMLE sample is extended to SNR>0.25.Since QMLE applies a weighted average based on SNR, this choice has minimal effect, which we confirm by comparing the QMLE estimates on the SNR>1 sample.The agreement between the two estimators is noteworthy given the built-in different methods of handling noise, metals, masking, etc.Since we find  = 2 bin is highly contaminated by various systematics, we recommend not using our QMLE results at  = 2 for any cosmological inference.Specifically, we find larger power than eBOSS on large scales, and we put substantial effort into identifying the root cause.The continuum fitting method slightly differs between eBOSS and DESI analyses, which could affect these large scales.However, our tests on mocks indicate that the effect of continuum errors on the power spectrum is small for DESI.The prime candidate then is systematics related to DLAs.We find that eliminating sub-DLAs (log  ≲ 20.3) removes power from large scales, but improvement in lower redshift bins comes with deterioration at higher redshifts.Furthermore, the DLA finder yields many false-positive, low-column density systems and is not reliable to identify such systems.The second candidate is a possible error in the noise power subtraction.However, our noise calibration correction does not fix this issue.We will study this disagreement further after additional development and testing of the DLA finder and how DLAs are included in our mock datasets, as well as with additional tests on the continuum fitting method.

Off-diagonal systematics in the covariance matrix
Proper use of our systematics budget requires special guidance.Chabanier et al. (2019) and Karaçaylı et al. (2022) added the systematics errors in quadrature to the diagonal of the covariance matrix.This may be true if systematic uncertainties are uncorrelated between bins, but is inconsistent with our characterisation.Simply put, we characterise each systematic error by an unknown scalar  multiplied with some () such that the uncertainty is in , not in ().For a given true signal (), the observed data is Both methods agree with each other overall, but disagree with eBOSS mostly at large scales at  ≲ 2.6.We suspect the efficiency of the DLA finder and the noise calibration could cause this discrepancy, and we will explore this further in future work.
given by  () = () +  (), and the covariance matrix becomes  ≡ ⟨ () ( ′ )⟩ = ⟨()( ′ )⟩ + ⟨ 2 ⟩()( ′ ).In terms of the systematic error budget we presented in Section 5.3, the covariance matrix should be modified as C total = C stat +  ∈ {syst}      for each systematic error mode   .This notation also highlights the mathematical relation with the continuum marginalisation procedure.Unlike the continuum marginalization procedure which adds infinitely large error modes   to fully remove contamination, systematic error modes are down-weighted by finite numbers.
Let us now discuss how these off-diagonal systematic error terms affect the parameters of interest.As we have noted, a full cosmological comparison is out of the scope of this work, and we instead use a fitting function based on equation (8).In our preliminary analysis, we found that not all parameters in equation ( 8) can be fitted properly and securely, so we simplify it by first ignoring the redshift evolution parameters  and .We instead fit each redshift bin separately and deduce the redshift evolution of each fitting parameter from these independent fits, though this redshift evolution itself is not pertinent to our discussion.We also found that the data is not sensitive enough to constrain the  1 parameter and it is highly degenerate with others.These two aspects of the  1 parameter destabilise the fitting, so we also eliminate it from our fitting function.We therefore instead use the following simplified fitting function: where  0 = 0.009 s km −1 as before.This formulation has simplistic yet useful characterisations such as parameter  as the amplitude and  as the slope of the power spectrum.Even though this procedure does not completely describe  1D , it is adequate to compare different measurements.
Another consideration in the fitting is that the Si ii and Si iii ions absorb at wavelengths that are close to the Ly transition line (Mc-Donald et al. 2006).The more visible oscillations in the power spectrum estimates are due to the Si iii line at  RF = 1206.5Å, which corresponds to a velocity spacing of Δ = 2270 km s −1 .Rather than complicating our fitting function with numerous Si ii lines (at 1190 Å, 1193 Å and 1194.5 Å) (Kramida et al. 2021), we take a single, average velocity spacing of Δ = 5770 km s −1 to model its pattern in the power spectrum.We multiply the base fitting function in equation ( 26) with 1 +  2  + 2  cos(Δ) for each Si ion , where   is the relative bias with respect to neutral hydrogen (McDonald et al. 2006;Palanque-Delabrouille et al. 2013).We find the least  2 solution using iminuit21 (Dembinski et al. 2022), and then sample around the minimum using the emcee22 sampler (Foreman-Mackey et al. 2013).
We apply this procedure to our metal-power subtracted power spectrum results.The covariance matrices are obtained with the bootstrap method.Figure 17 shows a corner plot at redshift  = 2.2 at 68% and 95% confidence levels generated using the getdist package (Lewis 2019).Results from adding the systematic error budget to the diagonals of the covariance matrix are shown in blue contours.Our recipe for the off-diagonal contribution (shown in orange) not only enlarges the confidence area (i.e.increases the error bars), but also shifts the best-fitting parameter value, which is particularly pronounced in the amplitude parameter .This has potentially important implications for cosmological inference.These off-diagonal covariance matrix contributions will be important to incorporate into future analyses.

Effects on parameters of each systematic
The systematic error budget diminishes the benefit of having a large and statistically powerful data set.As the data size increases, the statistical error bars decrease, but the systematic error budget stays the same and eventually saturates the constraining power of any analysis.Using the same off-diagonal recipe for the covariance matrix,  we break down how each systematic influences the error bar of inferred parameters.First starting with only the statistical covariance matrix from bootstraps, we add each systematic individually to the covariance matrix, find the least  2 solution, and sample around this minimum.Figure 18 shows the confidence regions using the statistical covariance in black dashed lines, where the noise systematics are in blue, the resolution systematics are in orange, and the DLA systematics are in green at  = 2.8.The noise and resolution systematics visibly enlarge the contours, and resolution systematics further shift the best-fitting values.The DLA systematics do not seem to increase the error bars significantly compared to the other two.Based on our analysis with the minimizer, we find that the error in the  parameter increases by 50% due to noise systematics, 34% due to resolution systematics, and 74% when all systematics are included.However, the slope parameter  is affected more unevenly: its error increases only 7% due to noise systematics, but 89% due to resolution systematics.Table 2 and Table 3  fully exploit its statistical power, we will have to prioritise mitigation of the noise calibration and spectrograph resolution systematics.As mentioned in the previous section, correlated CCD readout noise and sky subtraction errors can source noise systematics.The small number of quasars in the EDR+ sample limited our ability to perform ambitious data splits.The first-year data will have about a million quasars, which will enable a granular investigation.To remedy the noise systematics, we will first split the data into ten spectrograph subsamples, and then further subdivide the data within each spectrograph into two to isolate different CCD amplifiers.This division results in about 50,000 quasars in each region, approximating the sample size in this analysis, and so should provide enough statistics to study  1D in each subsample precisely.CCD image simulations could also be studied for noise recalibration parameters as well as refining the pipeline resolution estimates.

SUMMARY
The 1D Ly forest power spectrum  1D quantifies the clustering of diffuse neutral hydrogen gas in the intergalactic medium. 1D has been measured from various data sets, and has been used to constrain the thermal state of the IGM, the sum of neutrino masses, and various dark matter models.DESI will collect over 700 000  > 2 quasars during the five-year survey and will provide enormous statistical power for future  1D cosmology.It will be able to measure  1D from  = 2 to  = 5 and from approximate scales of 60 Mpc down to 1 Mpc.This power must be matched with rigorous studies of data and estimation methods.In this work, we employed the quadratic maximum likelihood estimator (QMLE) to measure  1D .QMLE is built to be statistically optimal and robust against Ly forest-specific challenges such as gaps in spectra and errors due to continuum fitting.Additionally, QMLE benefits from the resolution matrix output of spectro-perfectionism that preserves the full 2D resolution of the spectrograph.
In order to test the pipeline at various stages, DESI collected thousands of spectra during its Survey Validation phase.We used these early spectra to determine if the DESI  1D pipeline (such as noise calibration, resolution matrix, continuum fitting and  1D estimator itself) is accurate at the desired level and measure the initial  1D from DESI.
The quasar continuum estimation is a potential source of largescale uncertainty and bias.We described each quasar with two free fitting parameters (amplitude and slope) that multiply a mean continuum, and fit each continuum in the forest region.The two-parameter description is simplistic in terms of quasar continuum diversity; and furthermore fitting the continuum in the forest region removes the large-scale density information.Using synthetic spectra, we found that the estimated large-scale variance matches the input, but the estimated mean continuum deviated from the truth.Even though the estimated mean continuum was different, we showed that QMLE successfully marginalized out large-scale biases, and the residual power due to unmarginalized continuum error modes was significantly smaller than the signal.
The small-scale structure in the Ly forest is extremely important for the many science applications of  1D , and accurate knowledge of the spectrograph resolution is important to make the best use of the smallest scales.We created CCD image simulations with approximately 45 000  > 2.1 quasars (ten simulated observations of just quasars) and extracted the 1D spectra of these quasars using the DESI pipeline.We used these simulations to demonstrate that the resolution matrix produced by the pipeline was valid and did not require any modifications.We also showed that the PSF is well approximated by the pipeline's PSF model and that systematic errors due to PSF mismatch between the true PSF and the PSF fit by the pipeline were consistent with 1% precision on the resolution.
The noise power subtraction is not an insignificant part of the  1D estimation for DESI's medium-SNR spectra.Hence, any miscalibration of the pipeline noise is directly transmitted to the final results.We showed that the side band power subtraction that removes the metal power also eliminates some of this miscalibration, although not perfectly.Therefore, we examined the variance statistics of the transmission fluctuations in the relatively absorption-free side band regions to quantify the noise and flux miscalibration of the pipeline.We showed that these calibration errors are small (around a few percent) and corrected them in our analysis.We also found that the pipeline noise estimates are still coupled to the signal, and the noise calibration error depends on SNR.We do not try to correct for this SNR dependence and instead include this effect in the systematic error budget.
In order to have accurate error estimates, we relied on bootstrap realizations, which yielded larger statistical errors than their Gaussian counterparts from QMLE on almost all scales and redshifts.We further identified and quantified four major systematic error sources: noise, resolution, incomplete DLA removal, and continuum errors.Off-diagonal terms in the covariance matrix due to these systematic errors enlarge the error bars of inferred parameters and can also shift the best-fitting values, hence biasing cosmological interpretations.We showed that the noise calibration and resolution systematics weaken the statistical power of the current DESI early data, and will be priorities for additional study and mitigation for the Year 1 analysis.For the latter, we plan to conduct data split studies to quantify and mitigate noise systematics, create more extensive image simulations to quantify resolution systematics, and produce more synthetic datasets to better quantify the performance of the DLA finder and its impact on our results.
We found that the resulting  1D measurement from DESI early data and two months of main survey was in agreement between QMLE and FFT results, which is remarkable since the two estimators have different approaches for systematics.Furthermore, the two methods apply different SNR thresholds to the DESI sample.We confirmed that QMLE results are not affected by this threshold, which is expected since QMLE applies a weighted average based on SNR.These DESI QMLE and FFT  1D results are 5 − 15% larger than previous measurements from eBOSS and in 1.5 − 3 tension with those results.We investigated: 1) DLA finder efficiency and column density cuts; 2) noise calibration errors to explain this disagreement.However, none of these investigations offered a satisfactory explanation for the disagreement with eBOSS.This feature is worth further attention in future work.
To conclude, the DESI spectral pipeline works exceptionally well.The major analysis pipeline errors are no larger than few percent, and QMLE is well suited for DESI  1D measurements.The next five years will bring incredible accuracy, precision, and power to Ly forest cosmology.any boundary effects.We convolve the resulting noise estimates by a hybrid Gaussian box-car window function, where the smoothing kernel has a size of 51 pixels and a Gaussian sigma of 20 pixels.After applying this hybrid boxcar-Gaussian smoothing kernel, we return outlier values to their original positions.This smoothing broadly captures the spectrograph behaviour while still down-weighting the masked or high-variance pixels.
Analytical expressions: The pipeline noise estimate of a CCD pixel depends on the exposure time  exp , and has the following contributions: • Source electrons  source =   () () ∝  exp .Note  =  (1 + ) • Sky electrons  sky ∝  exp .
• Dark electrons ∝  exp .We can absorb this term into sky contribution.
• Read noise electrons as Gaussian noise with constant  2  independent of  exp .
Then, pipeline variance estimate on flux  is  2  ,pipe =  source +  sky +  2  .We can write the pipeline noise estimate on  as follows: Using inverse variance weights means   ∝ ( + Λ) −1 .Since 0 <  < 1, our weights are highly correlated with the signal if Λ ≲ 1. Unbiased limit Λ ≫ 1 is satisfied when the sky dominates the signal (noisy spectra) or the quasar is so bright that the LSS term dominates the variance.
There are two important points to note.First, strictly speaking, pipeline noise estimates are correct.The problem is not an error in the pipeline but in the nature of Ly forest.As we mentioned in the main text, DESI decouples signal and noise by smoothing source contribution at 10Å scale (Guy et al. 2023).Second, noise is still diagonal, and there are no cross correlations between signal and noise.Noise on each pixel is another independently generated random number.Therefore, even though its amplitude depends on the signal, this dependence does not introduce auto or cross correlations to noise.DESI already has expected values for these quantities.Figure B1 shows coupling parameter Λ vs quasar SNR.Low SNR quasars manifest less coupling (high Λ) since they are dominated by sky and read noise.High SNR quasars also manifest less coupling as the LSS term dominates.On average, DESI expects the mean quasar SNR to be 0.45, which corresponds to Λ = 34.The weighted average estimator for the mean flux is given by F =     /   , where   = (  + Λ) −1 , where  is signal only as before.The expected value of this estimator can be calculated by using the probability distribution function P ().
Finally, the error on the mean flux estimate is The same calculation can be done for two-point statistics.Gauss-Hermite quadrature is still the best way to numerically calculate the expected value.We can either use the true mean flux or use the biased mean flux estimates from weighted averages to calculate s.It is worth noting that using biased mean flux does not scale but shifts estimated deltas.This does not make a big difference and only affects the  = 0 mode of the power spectrum.However, the correlation function asymptotically approaches a constant at large scales.Assuming true deltas are used with biased weights, the error on the two-point function estimate is In other words, estimated two-point statistics are contaminated by three-point statistics and distorted at all scales through correlations between weights.
We have numerically calculated both expressions.For brevity, Figure B2 shows only biased  1D with respect to different Λ values.Using true mean flux to extract deltas (dashed lines) shows slightly less bias than using biased delta extraction (solid lines).For comparison, our estimates from mocks without smoothing fall close to Λ = 50 (red line) as expected from the fiducial value.As we note in the main text, smoothing the noise estimates solves this problem.
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.Quasar at  = 2.94 observed during DESI survey validation (TargetID 39627871806818826).The Ly forest is defined to be the spectral region between a quasar's Ly and Ly  emission lines.Absorption features redward of the quasar's Ly emission line may be due to metal systems.The regions from Ly to Si iv and from Si iv to C iv are called the "side bands" (SB).We call the Ly -Si iv region SB 1 and the Si iv-C iv region SB 2, and use these regions to quantify metal contamination, noise and flux calibrations.

Figure 2 .
Figure 2. (Top) Quasar redshift distribution of our sample.(Bottom) Weighted mean SNR distribution in the forest as a function of redshift, where we define SNR to be 1/ () and  () is the propagated error on the weighted mean of pixel size of Δ = 0.1.Having BAL quasars (orange) improves SNR, but it also comes with possible biases in  1D (see Section 2).

Figure 3 .
Figure 3. Reduced  2 comparison for different  cuts and continuum marginalisation polynomials on mocks.We find  2  increases for all settings as we include higher  values, which is unfortunate but expected since our correction to the quickquasars resolution matrix is not exact.The true continuum analysis results (blue triangles) stay within 1.5 of  2  = 1.Lower rows correspond to larger small-scale confidence regions.From left-most column to the right, we remove large-scale modes.When continuum errors are not marginalised (orange squares), throwing out these large-scale modes brings  2  down to 1 within error bars.We also find that first (green circles) and second (red triangle) order marginalizations remove the contamination from continuum errors at all scales.

Figure 4 .
Figure 4. (Top) Mean global continuum from the true continuum analysis vs. from continuum fitting on mocks.The continuum fitting accentuates features in the mean continuum.(Bottom)  2 LSS values from the true continuum analysis vs. continuum fitting.Our fitting algorithm finds the correct  2 LSS

Figure 5 .
Figure 5. Power spectrum of the remaining continuum errors after marginalization, divided by the true underlying mock Ly power spectrum.We find the remaining continuum errors are a factor of 10 −5 smaller than the signal at most scales and redshifts.The maximum is 10 −3 at the largest scales.

Figure 6 .
Figure 6.Power spectrum when DLAs and sub-DLAs are not masked divided by the true underlying mock power spectrum.These systems add power to large scales.

Figure 7 .
Figure 7. (Top) Power spectrum estimates for various resolution matrix treatments at  = 2.8 based on CCD image simulations.The default pipeline output (blue squares) performs best, whereas oversampling (green line) deviates from the truth.Padding the edges of the resolution matrix (orange triangles) improves the agreement at the largest scales, but those modes are lost to continuum errors in any case.(Bottom) The difference between power spectrum estimates using the fitted PSF and the true model PSF, divided by the true underlying signal.The fitted PSF introduces a wavelength dependent resolution error.

Figure 9 .
Figure 9. Ratio of noise power to noise-subtracted Ly power from data.The noise power spectrum is not negligible even at large scales.Results at  = 2.0 are most sensitive to the noise power estimates.The features in this ratio mostly come from the inverse of  1D .

Figure 10 .
Figure 10.Ratio of noise power in SB 1 to Ly region from data.It is redshift and scale dependent.Noise power is smaller in side band regions, which means pipeline noise miscalibrations cannot be fully removed by side band subtraction.
>  min and    ≡    / √︁      . 20Number of sub-samples is based on the MPI tasks used, which is a current code limitation.

Figure 11 .
Figure11.Ratio of the regularized bootstrap error estimates to Gaussian (QMLE) errors for Ly .We estimate the bootstrapped covariance matrix from 200 000 realizations over 256 subsamples.This shows that QMLE underestimates the errors on most scales and redshifts.

Figure 15 .
Figure 15.Systematic error budget divided by the statistical errors estimated by bootstrap analysis.Systematics due to the noise calibration is persistent on all scales and prevalent over the statistical error at lower redshifts.Resolution systematics predominantly influence small scales (high ).

Figure 16 .
Figure16.Comparison of  1D between FFT and QMLE methods on DESI EDR+ and to eBOSS.Both methods agree with each other overall, but disagree with eBOSS mostly at large scales at  ≲ 2.6.We suspect the efficiency of the DLA finder and the noise calibration could cause this discrepancy, and we will explore this further in future work.

Figure 17 .
Figure17.Best-fitting parameters and 60% and 95% confidence contours when the covariance matrix has systematic error contributions in its offdiagonal elements (orange) and in its diagonal only (blue) at  = 2.2.Offdiagonal terms shift the best-fitting values and increase the error bars, which has potentially important implications for cosmological inference.

Figure 18 .
Figure18.The 60% and 95% confidence regions using only the statistical covariance (black dashed lines), and individually adding the noise systematics (blue), resolution systematics (orange) and DLA systematics (green) at  = 2.8.DLA systematics do not seem to increase the error bars as much as the other sources of systematic error.The noise and resolution systematics visibly enlarge the contours, and the resolution systematics further shifts the best-fitting values.

Figure B1 .
Figure B1.Coupling parameter Λ vs quasar brightness   as electrons per Angstrom for DESI fiducial values at 3800Å.The upper axis is the human-readable SNR of the quasar.DESI expects mean SNR=0.45, which corresponds to Λ = 34 (low coupling).

Figure B2 .
Figure B2.Biased power spectrum estimates for different Λ values.Using true mean flux to extract deltas (dashed lines) shows slightly less bias than using biased delta extraction (solid lines).For comparison, our estimates from mocks without smoothing fall close to Λ = 50 (red line) as expected from the fiducial value.

Table 1 .
Rest-frame wavelength ranges and number of quasars in DESI early data.Quasars with BAL features are ignored.

Table 2 .
Percentage increase in error given by the minimizer for each systematics at  = 2.8.The precision of the amplitude  is nearly equally affected by noise and resolution systematics, whereas for , it is thoroughly affected by resolution systematics.

Table 3 .
Best-fitting values and error estimates given by the minimizer for each systematics at  = 2.8.
detail these numbers.Year 1 DESI spectra will provide even more statistical power.To The total variance includes large-scale structure fluctuations  2 LSS , which should be multiplied by the mean IGM flux .