DESI DR1 Ly$\alpha$ 1D power spectrum: The Fast Fourier Transform estimator measurement

The Dark Energy Spectroscopic Instrument (DESI) is a multi-object spectrograph whose purpose is to measure the spectra of more than $40$ million galactic (stars) and extra-galactic objects (galaxies and quasars). The main science goal of DESI is to infer the properties of dark energy using Baryon Acoustic Oscillations (BAO) and Redshift Space Distortions (RSD) measured in galaxy, quasars, and Ly$\alpha$ correlation functions. Additionally, the significant statistics and variety of observations allows DESI to cover broader scientific objectives, including the one presented in this article. The focal plane of DESI is equipped with $5,000$ robotic fiber positioners [25, 46, 47, 48, 49] to quickly reconfigure the observation pattern and point towards pre-selected targets [50]. A full description of the DESI instrument is given in [26]. We use the first data release noted DESI-DR1 (or DR1 for conciseness), which results from the first year of DESI observation [39].

Quasars candidates are pre-selected with a dedicated targeting to obtain a denser sample of high-redshift targets [50, 51, 52]. The DESI-DR1 quasar sample is constructed by automated classifiers. A first template fitting software called redrock \faGithub¹¹1https://github.com/desihub/redrock [53] is used to classify quasars, galaxies, and stars, and to provide a robust redshift measurement. A dedicated neural network software called QuasarNet \faGithub²²2https://github.com/desihub/QuasarNP [54, 55] picks up quasars missed by redrock. As QuasarNet is designed explicitly for classification, the redshift of retrieved quasars is determined by running redrock again with an informed prior. Finally, a Mgii line fitter algorithm catches low-redshift quasars ($1.5<z<2.0$) whose redshifts were correctly measured, but which have been classified as objects of another type. The full description of the pipeline used to generate the quasar catalog is given in [52]. The histogram of the DR1 quasar redshifts is represented in figure 1. It contains $521,509$ quasars with a redshift in the Ly$\alpha$ forest range of interest ($2.0<z<5.0$).

Several modifications are made to this quasar catalog. First, spectra from secondary target programs dedicated to the search for high-redshift quasars are added. A visual inspection campaign of those targets is carried out (see the companion paper K25 for more details). Second, we remove quasars whose measured individual power spectrum appears as an outlier. The full DR1 sample cannot be entirely visually inspected and still contains outliers passing the automated classifier criteria. We carry out a pre-study in which the individual power spectrum of each Ly$\alpha$ forest is measured. Taking the full DR1 catalog, we remove quasars for which the average of individual power spectrum over large scales ($k<0.8$ Å^-1) is more than $20\sigma$ higher than the average for all quasars, where $\sigma$ is the standard deviation over all quasars. This analysis is carried out in the Ly$\alpha$ and side-band SB1 ($1,270<\lambda_{\mathrm{rf}}<1,380$ Å) regions used later in this paper. A visual inspection showed that the removed spectra are indeed quasar spectra but exhibiting spikes caused by instrumental artifacts. A small number of $7$ and $18$ quasars are removed by conducting this study in the Ly$\alpha$ and SB1 region, respectively. This removal eliminates spurious oscillations in the measurement of the one-dimensional power spectrum.

The computation of one-dimensional power spectra necessitates accounting for several effects that causes changes in the spectra (absorptions or emissions), and that are not directly related to IGM absorptions. Our choice concerning these contaminants is to remove the spectra from the data sample or mask the impacted spectral region.

First, the Broad Absorption Lines (BAL) quasars present specific blueshifted absorptions caused by the quasar torus, thus directly related to its surroundings. These specific quasars are found using the baltools \faGithub³³3https://github.com/paulmartini/baltools software [56], which is a $\chi^{2}$ minimizer that measures blueshifted Civ or Siiv absorptions compared to an unabsorbed quasar continuum. The baltools software computes for each quasar spectrum the Balnicity Index (BI) and Absorption Index (AI) [56] defined as

Specific absorptions called High-Column Density (HCD) systems are present in the spectrum of some quasars. They are caused by the passage of the quasar light in the vicinity of a galaxy, called a circumgalactic medium. HCD objects contaminate the Ly$\alpha$ forest by imprinting a saturated absorption that spans a much larger spectral range than the size of the galaxy itself and adds metal absorption lines present in the circumgalactic medium [57]. The HCD objects are detected using the combination of a convolutional neural network (CNN) algorithm desi-dlas \faGithub⁴⁴4https://github.com/cosmodesi/desi-dlas [58] and a Gaussian process (GP) algorithm \faGithub⁵⁵5https://github.com/jibanCat/gpy_dla_detection [59]. The GP algorithm provides a better estimate of the HCD column density $N_{\mathrm{H{\textsc{i}}}}$, and the CNN a better detection completeness. We only consider the Damped Ly$\alpha$(DLA) objects, defined by $N_{\mathrm{H{\textsc{i}}}}>10^{20.3}~\mathrm{cm}^{-2}$. We take the same criteria as R23 to define our sample: valid DLA detection by CNN with a confidence level higher than $0.2$ for quasars continuum-to-noise ratio $\overline{\mathrm{CNR}}>3$ and a confidence level higher than $0.3$ when $\overline{\mathrm{CNR}}<3$. For the GP algorithm, we take a minimal $0.9$ confidence level. The final catalog results from the merging of DLA detected by either algorithm. When both algorithms detect the DLA, the GP column density is used. The core region of a detected DLA, defined by an induced absorption larger than $20~\%$, is masked at the spectrum level. The remaining DLA damping wings are corrected at the spectrum stage using a Voigt profile (see e.g. [27] for details regarding DLA masking). Approximately $16$ % of the quasars in the DR1 catalog have at least one detected DLA in their spectrum. The redshift distribution of quasars presenting at least one detected DLA and that of quasars passing the $AI>0$ or $BI>0$ BAL quasar criteria are shown in figure 1.

Finally, absorption lines caused by passing through the Milky Way and atmospheric emission lines cause spurious spectra modifications even after being corrected by the DESI pipeline. We mask the lines in the line catalog adapted to DESI spectral resolution and developed in R23(\faGithub⁶⁶6https://github.com/corentinravoux/p1desi/tree/main/etc/skylines/list_mask_p1d_DESI_EDR.txt). The transmitted flux fraction of masked pixels (DLA or lines) is fixed to its average value ($F=\overline{F}$), which corresponds to a zero flux contrast. This masking induces a bias on $P_{1\mathrm{D},\alpha}$ that will be corrected in section 7.

The continuum fitting procedure estimates the flux contrast $\delta_{F}$, i.e., the normalization of the flux variation caused by all absorptions in the Ly$\alpha$ forest, from the measured flux $f$ at a specific wavelength $\lambda$:

The noise associated to each flux contrast is directly given by the estimate of the DESI pipeline noise, noted $\sigma_{\mathrm{pip,q}}(\lambda)$ (see [61] and R23), normalized by the $C_{\mathrm{q}}(\lambda,z_{\mathrm{q}})\overline{F}(\lambda)$ product. We then define the mean signal-to-noise ratio of a quasar $q$ as

Following R23, we only keep quasars with $\overline{\mathrm{SNR}}>1$. The fitting procedure is performed for an observed wavelength range of $3,600<\lambda<7,600$ Åand a rest-frame range of $1,050<\lambda_{\mathrm{rf}}<1,180$ Å. The quasar spectra are linearly binned in observed wavelength with $\Delta\lambda_{\mathrm{pix}}=0.8$ Å. We adopt the same rebinning scheme of the common continuum $C$ in the rest-frame basis with a bin size equal to $\Delta\lambda_{\mathrm{pix,rf}}=2.67$ Å. Finally, the stack of all flux contrasts is set to zero at the end of the procedure.

The one-dimensional Ly$\alpha$ power spectrum, noted $P_{1\mathrm{D},\alpha}$, is defined as the power spectrum of the Ly$\alpha$ absorption contrast $\delta_{\mathrm{Ly\alpha}}$ along the quasar line-of-sight. Following R23, we include in the $P_{1\mathrm{D},\alpha}$ definition the cross-correlations between the Ly$\alpha$ absorption line and other IGM elements whose rest-frame absorption is very close to the Ly$\alpha$ line. In particular, we include absorptions from Siii and Siiii ($\lambda_{\mathrm{Si{\textsc{ii}}}}=1,190$ and $1,193$ Å, and $\lambda_{\mathrm{Si{\textsc{iii}}}}=1,206.50$ Å) whose cross-correlation with Ly$\alpha$ absorption is modeled at the cosmological interpretation stage. Neglecting the Siii and Siiii auto-correlations as they correspond to minor metal lines, we mathematically define $P_{1\mathrm{D},\alpha}$ as

To compute $P_{1\mathrm{D},\alpha}$, we first separate the flux contrast $\delta_{F}$ from each spectrum into three non-overlapping sub-forests of equal size, $L_{\mathrm{sub}}=43.3$ Å in the rest frame, with a flux contrast noted $\delta_{F,s}$. Each of the $3N_{\mathrm{q}}$ sub-forests (with $N_{\mathrm{q}}$ the total number of quasars) is referred to with the index $s$ in the following equations. This splitting reduces correlations between redshift bins and increases the smallest accessible wavenumber to $k_{\mathrm{min}}=2\pi/(L_{\mathrm{sub}}(1+z_{\mathrm{min}}))=0.0468$ Å^-1. Furthermore, the resulting sub-forest covers at most $\Delta z=0.2$, and we use this value to define the redshift binning for $P_{1\mathrm{D},\alpha}$. The redshift associated with a sub-forest is determined by the center of its observed wavelength range, which assigns it to a redshift bin. Therefore, each redshift bin contains information from adjacent redshift bins, as sub-forests with a central redshift near a bin edge have half their spectrum extending beyond their redshift bin. It implies the existence of an inter-redshift correlation, which is not accounted for in this work but will be in future studies. The sub-forest flux contrasts are then transformed to Fourier space by picca using a multiprocessed Fast Fourier Transform (FFT) algorithm.

The $P_{1\mathrm{D},\alpha}$ FFT estimator is built under the main assumption that each sub-forest in our dataset provides an independent realization of the IGM Ly$\alpha$ absorptions at the redshift at which they occur. In other words, we assume that the average of all sub-forests gives an unbiased estimator of the one-dimensional power spectrum, and we neglect correlations between sub-forests. The FFT estimator is derived from the relation between the observed flux contrast $\delta_{F,s}$ obtained after continuum fitting and the Ly$\alpha$ contrast $\delta_{\mathrm{Ly\alpha},s}$. This relation is derived by considering all the instrumental (noise, resolution) and astrophysical (metals) effects that imprint absorptions in the Ly$\alpha$ region. Following the decomposition detailed in R23, the FFT $P_{1\mathrm{D},\alpha}$ estimator for the wavenumber bin $A$ and redshift bin $z$ is defined by

The picca software is used to compute this FFT estimator. In agreement with R23, we fix our minimal redshift to $z_{\mathrm{min}}=2.1$ to avoid the $\lambda\lesssim 3,700$ Å region of the spectra where atmospheric absorptions significantly increase the flux noise. As we apply several masks (DLA, BAL, and atmospheric emission line catalogs), we remove sub-forests shorter than 75 wavelength pixels or with more than 120 masked pixels. We compute the raw power spectrum for each sub-forest directly from the flux contrast as $P_{\mathrm{raw},s}=\left|\delta_{F,s}(k)^{2}\right|$. The noise power spectrum $P_{\mathrm{noise},s}$ is computed using the DESI pipeline noise. The software generates a large number ($N_{\mathrm{G}}=2500$) of Gaussian signals $\delta_{\mathrm{noise},s}$ with standard deviation equal to $\sigma_{\mathrm{pip,s}}(\lambda)$ for each wavelength, and the noise power spectrum is obtained by averaging those signals in Fourier space, i.e. $P_{\mathrm{noise},s}(k)=\left\langle\left|\delta_{\mathrm{noise},s}\right|^{2}\right\rangle_{N_{\mathrm{G}}}$.

The picca software averages the individual power spectra with a $\overline{\mathrm{SNR}}$ weighting scheme. In each bin $(A,z)$, we distribute the power spectra of individual sub-forests in bins in the signal-to-noise ratio of the sub-forests $\overline{\mathrm{SNR}}_{s}$. We compute the variance of the power spectra in each $\overline{\mathrm{SNR}}_{s}$ bin. We then fit these variances with a simple least-squares minimization as a function of $\overline{\mathrm{SNR}}_{s}$ according to:

The fitted terms $a_{A,z}$ and $b_{A,z}$ are used to define the inverse-variance weights of each sub-forest as $w_{A,s}=1/\sigma_{(A,z)}^{2}\left(\overline{\mathrm{SNR}}_{s}\right)$. By construction, the weights tend to zero when $\overline{\mathrm{SNR}}$ tends to 1, in agreement with the applied $\overline{\mathrm{SNR}}>1$ cut. Finally, the weighted average in equation 3.6 is explicitly computed as:

Some methods applied to the analysis of EDR+M2 DESI datasets in R23 are applied in this paper to the DR1 data set without any modifications. We review the new characterization of those instrumental features with the DR1 data set. The updated figures are located in appendix A.

The metal power spectrum, $P_{\mathrm{metals}}$ in equation 3.6, is estimated using the side-band region SB1 ($1,270<\lambda_{\mathrm{rf}}<1,380$ Å). We show the SB1 power spectrum in figure 16, which is measured using the averaging of equation 3.6 in the SB1 rest-frame range. Here, we use a rest-frame range distinct from the Ly$\alpha$ forest one, but with the same observed wavelength range. The quasars selected for SB1 are then located at lower redshifts than those selected for Ly$\alpha$. Consequently, our correction assumes that the continuum fitting process can catch the quasar emission variability as a function of redshift. As in R23, we average the side-band power spectrum over redshift bin on a common rest-frame wavenumber binning $k_{\mathrm{rest}}=(1+z)\times k_{\mathrm{obs}}$. This average is fitted with a least-square minimization by a physically motivated model considering Siiv and Civ doublets oscillations:

The evolution of the SB1 power spectrum as a function of redshift is accounted for in a second step, by multiplying and fitting a linear function independently for each redshift, $P_{\mathrm{SB1,m,z}}(k)=(a_{z}k+b_{z})P_{\mathrm{SB1,m}}(k)$. The resulting fitted functions are used as the metal power spectrum correction in equation 3.6. The fitted parameters for SB1 power spectrum are similar to the EDR+M2 ones in R23. The DR1 power spectrum contains a larger number of sub-forests (by at least a factor $5$ for all redshifts) and consequently gives lower error bars.

The R23 method is applied to measure the resolution damping with DR1 dataset. The resolution damping term in equation 3.6 is the average of the Fourier transform of the resolution matrix. The average resolution damping is shown in Fig. 15 and is almost identical to the one measured for EDR+M2 data set in R23. Consequently, we choose the same maximal wavenumber value $k_{\mathrm{max}}=2.0$ Å^-1, corresponding to a factor 5 correction.

We measure the noise power spectrum with the method detailed in section 3.2 directly from the DESI pipeline output. As highlighted in R23, characterizing the noise sources in DESI is challenging, and the noise level can be underestimated. To have an alternative way to estimate the noise level, we use the smallest scales of $P_{1\mathrm{D},\alpha}$, for which the resolution damping completely suppresses the signal, resulting in an expected equality between the raw and noise power spectra. Consequently, the measurement of $P_{\mathrm{raw},s}$ at the largest wavenumber is a suitable noise level estimator. The noise and raw power spectra are compared in figure 2 for the DR1 dataset.

We derive an additive noise correction by averaging the $P_{\mathrm{raw},s}-P_{\mathrm{noise},s}$ difference for wavenumber where the resolution damping is higher than 98 % ($k_{\mathrm{res},98}=3.1$ Å^-1), which gives $\alpha=\left\langle P_{\mathrm{raw},s}(k)\right\rangle_{\forall z,k>k_{\mathrm{res},98}}-\left\langle P_{\mathrm{noise},s}(k)\right\rangle_{\forall z,k>k_{\mathrm{res},98}}$. The measured $\alpha=6.9\times 10^{-4}$ Å for DR1 corresponds to a 2.8 % noise power spectrum underestimation. This is lower than the value $\alpha=1.09\times 10^{-3}$ obtained for DESI-M2 in R23, a dataset with similar properties, which means we have better control over the noise characterization thanks to the dataset quality. In contrast to R23, the noise power spectrum is not flat for $k\lesssim 2.0$ Å^-1. This is due to the $\overline{\mathrm{SNR}}$ weighting used for DR1 and not for R23. Indeed, the noise power spectrum is flat for all individual power spectra, but since we are considering different weights in the $k\in[0.0-2.0]$ Å^-1 range, the $\overline{\mathrm{SNR}}$ weights distribution impacts the averaging of the noise power spectrum in a wavenumber-dependent way.

To control the noise power spectrum even better, we developed an alternative power spectrum estimator called cross-exposure, which does not need to model or measure the noise power spectrum explicitly. The DESI standard pipeline coadds the multiple exposures of a given quasar to create less noisy spectra [61]. In contrast, the cross-exposure estimator computes the one-dimensional individual cross-power spectrum between all distinct exposures of a given quasar. This provides an estimator free from noise influence if the noise realizations are independent between exposures.

Therefore, to minimize the impact of correlated noise sources, we select exposures to get only one per quasar and per night. It removes the noise sources resulting from exposures using the same calibration images, such as bias and flat calibrations, but not necessarily that resulting from dark calibration, which is not performed each night (see e.g. [61] or appendix C of R23 for an extensive description of noise sources).

For a quasar with $N_{\mathrm{exp}}$ exposures in DESI, we define the following cross-exposure operator:

We replace the simple raw power spectrum coadd estimator given in section 3.2 by the cross-exposure estimator $P_{\mathrm{raw},s}(k)=\mathcal{P}_{X}\left[\vec{\delta}_{F,s}(k),\vec{\delta}_{F,s}(k)\right]$, where $\vec{\delta}_{F,s}$ is the set of sub-forest flux contrast for the $N_{\mathrm{exp}}$ exposures.

Noise correction and resolution damping must be accounted for to build a one-dimensional power spectrum estimator. Following R23, the flux contrast can be decomposed in Fourier space as $\vec{\delta}_{F,s}(k)=\left(\vec{\delta}_{\mathrm{Ly\alpha},s}(k)+\vec{\delta}_{\mathrm{noise},s}(k)\right)\times\vec{R}_{s}(k)+\vec{\delta}_{\mathrm{metal}}(k)$. Following the demonstration in R23 to build the $P_{1\mathrm{D},\alpha}$ estimator, the cross-exposure estimator $P_{1\mathrm{D},\alpha,\mathrm{cross}}$ yields:

In this decomposition, since the noise contrasts are purely random Gaussian signals, they are all independent from each other and from Ly$\alpha$ contrasts. Consequently, the last two terms of the decomposition are null. We numerically verified this assumption on the DR1 data set by assessing that the average of the last two terms is consistent with zero. Consequently, the cross-exposure estimator is given by:

We use a new implementation in picca to compute this cross-exposure estimator by performing the continuum fitting on each exposure of all the DR1 quasar catalog. Since distinct exposures of the same quasar have different $\overline{\mathrm{SNR}}$, we perform the continuum fitting independently, as if each exposure were an independent quasar. We tried alternative continuum fitting methods and concluded that using the continuum fitted on the coadd of all exposures or on the exposure with the highest $\overline{\mathrm{SNR}}$ introduces biases and spurious oscillations in $P_{1\mathrm{D},\alpha,\mathrm{cross}}$.

Figure 3 compares the one-dimensional power spectrum measured on DR1 dataset before noise correction with $P_{1\mathrm{D},\alpha}$ after noise correction, and with the cross-exposure estimation, $P_{1\mathrm{D},\alpha,\mathrm{cross}}$. The two bottom panels are showing the ratios between those three analysis variations. At large wavenumber ($k\gtrsim 1.5$ Å^-1), the cross-exposure estimator has the same trend as the noise-corrected measurement, indicating that both noise estimations match. The uncorrected $P_{1\mathrm{D},\alpha}$ have a discrepancy at small redshifts for those scales, even if the difference is still within the statistical error bar. Near $k\sim 0.6$ Å^-1, both uncorrected and corrected measurements show a slight discrepancy with the cross-exposure estimator. This bump is caused by a pattern in a specific faulty DESI amplifier and is studied in detail in section 5.4. The cross-exposure estimator can remove part of this contamination at low redshifts.

The cross-exposure estimator only uses between $35$ and $39$ % (depending on the redshift bin) of the quasar catalog because of the need for several exposures and the fact that we select exposures from different nights. Consequently, we do not directly use the cross-exposure estimator for high redshifts, where the statistical errors are large. As discussed in 7, to better control the noise estimate, we use the cross-exposure estimator as the baseline for the redshift bins where we see small-scale and amplifier discrepancies, i.e. $z=2.2,$ $2.4$, and $2.6$. For the other bins, we use the noise-corrected measurement as it agrees with the cross-exposure but results in smaller error bars due to the larger number of sub-forests and exposures available. Additionally, considering the extensive studies we performed here and in R23, we assert here that we have sufficient control over the noise level to drop the associated systematic error bar, which was already overly conservative in R23. We assume that the increase in the statistical error at low redshift due to using the cross-exposure estimator makes the systematic error due to noise negligible.

The DLA catalog used for masking $P_{1\mathrm{D},\alpha}$ was created with the specific CNN and GP confidence cuts in section 3. Variations of those parameters provide insight into DLA completeness. The left panel of figure 4 shows the two variations related to the DLA catalog. We start by not masking DLA (NOMASK) and comparing it to the baseline. In agreement with the previous measurement in R23 or on simulation in [63], not masking DLAs greatly increases the power spectrum at the largest scales considered. This test is used in section 6.2 to derive a systematic uncertainty associated with the DLA completeness.

Conjointly with the companion paper K25, we modify the confidence cut to obtain a higher purity version of the DLA catalog. This alternative cutting uses the same confidence cut ($CL>0.9$) for GP, and a continuum-to-noise ratio ($\overline{\mathrm{CNR}}$)-dependent cut for the CNN: $CL>0.2$ for $\overline{\mathrm{CNR}}>3$, $CL>0.4$ for $2<\overline{\mathrm{CNR}}<3$, and $CL>0.5$ for $\overline{\mathrm{CNR}}<1$. The resulting catalog contains approximately $31$ % less DLA than the high-completeness version used for the baseline. When using this high-purity catalog, the obtained power spectrum is compared to the baseline in the last panel of Figure 4 (PURITY). The $p$-values reported show a statistically significant correlation with respect to wavenumber, coming from small wavenumbers and small redshifts. We note, however, that this difference is small (around 2$\%$) and is well below the level of the systematic uncertainties that we derive in section 6.2 for DLA completeness. Thus, employing a purer or more complete sample of DLA does not largely impact our measurement. Furthermore, even if the two catalogs give different DLA numbers, the impact of masking is not significantly changed. We choose to keep the high-completeness catalog as the baseline.

We have a larger number of sub-forests than in R23, sufficient to study the effect of the completeness of BAL catalog, as illustrated in the right panels of figure 4. The baseline measurement (BASE) uses a catalog in which quasars tagged as BAL with the $BI>0$ criterion have been removed. The upper panel of the figure compares it to an analysis variation (AICUT) with the $AI>0$ cut, which removes $22$ % of quasars, including the $4$ % of $BI>0$ quasars. We observe a near $5$ % difference at small scales, well above the level of statistical uncertainties, giving a $0.0$ Correlation $p$-value. However, removing all $AI$ BAL quasars would induce a large drop in statistics, and we aim to match the power spectrum obtained with this severe cut while keeping most of the quasars.

Two additional variations are then performed: masking at the spectrum level all $AI$ BAL absorptions while keeping all quasars in the catalog (AIMASK) and masking all $AI$ BAL pixels while removing $BI$ BAL quasars from the catalog (BIC_AIM). AIMASK is a drastic change relative to baseline as no quasars are removed from the input catalog, while BIC_AIM can be considered as the baseline while additionally masking for BAL regions that are passing $AI>0$ criterion. As shown on the lower two right panels of the same figure, the impact of masking BAL is very similar for AIMASK and BIC_AIM: it decreases the $P_{1\mathrm{D},\alpha}$ level for all scales due to the effect of pixel masking while accounting for the small scale difference due to the unaccounted $AI$ BAL that are masked. Since the impact of masking is way more prominent when the $BI$ BAL quasars are kept (AIMASK), we choose BIC_AIM as our new baseline.

To verify that the impact of $AI$ BAL quasars is correctly accounted for, we compare our new baseline to the case where all $AI$ BAL quasars are removed from the catalog (AICUT). The new baseline is corrected for the impact of BAL masking using the correction derived in KR25. The number of $AI$ BAL quasars is larger in the DR1 data than in the mocks used in KR25. We normalize the absolute level of the BAL masking correction to account for the larger fraction of BAL in the data. The resulting new baseline (BIC_AIM with corrections) is compared to the $AI$ cutting case in figure 5. Compared to the first left panel of figure 4, the new baseline properly accounts for the $AI$ BAL, as shown by the $p=0.05$ value. This p-value is the lower value for six bins in redshift, so $p=0.05$ indicates that the test is passed for all redshift bins. We will use this study of the impact of $AI$ BAL to derive systematic uncertainties linked to BAL completeness and the BAL masking correction.

We vary some parameters chosen as fixed in section 3 to blindly search for potential sources of analysis systematic error. First, the parameters used for continuum fitting in section 3.1 are varied, and the ratio with baseline is shown in the left panels of figure 6. We carry out the continuum fitting while allowing noise-dependent weights in the likelihood used for fitting (VARW). In equation 3.3, the weights $\sigma_{q}=\overline{F}(\lambda_{i})C_{\mathrm{q}}\left(\lambda_{i},z_{\mathrm{q}},a_{\mathrm{q}},b_{\mathrm{q}}\right)$ in the likelihood are replaced by:

The four variations considered do not induce a significant decrease in the number of sub-forests: only a nearly $5$ % decrease for LSMALL due to the creation of shorter sub-forests that do not pass the size cuts. The second test (ORDER2) has a very small effect on $P_{1\mathrm{D},\alpha}$ as indicated by the reported $p$-values. All the three other tests (VARW, LSMALL, SFCUT) shows some statistically significant correlation between the ratio and the wavenumber. The VARW test has a small but significant effect on the power spectrum. However, including $\sigma^{2}_{\rm LSS}$ in the weights is biasing the evaluation of the power spectrum, so this effect is not surprising. The LSMALL and SFCUT variations show more significant fluctuations. We interpret those spikes as caused by adding smaller sub-forests that appear as outliers in the averaging, and that are potentially not accounted by the outlier finder detailed in section 2. Those two variations should be considered in a cosmological interpretation study as they could potentially yield different results.

In a separate study, we test the average $P_{1\mathrm{D},\alpha}$ computation by varying the $\overline{\mathrm{SNR}}$ cut and the weighting scheme. The right panels of figure 6 show the variation performed here: using a constant weighting scheme with the same $\overline{\mathrm{SNR}}>1$ cut (NOW), with a $\overline{\mathrm{SNR}}>3$ cut (NOWSNR3) which is close to what was done in eBOSS [45], or using the same $\overline{\mathrm{SNR}}$ weighting scheme as baseline with the $\overline{\mathrm{SNR}}>3$ cut (SNR3). The $\overline{\mathrm{SNR}}>3$ cut implies keeping only between $32$ to $40$ % of the DR1 sub-forests depending on the redshift bin. It corresponds to a drastic cut in statistics, but can be considered a robust sample for small scales, for which noisy sub-forests can significantly impact the measurement. Indeed, noisy sub-forests have a large noise power spectrum and a slight noise misestimation causes a large small-scale variation. The upper plot exhibits a large discrepancy between the baseline and the case without weighting.

For the second plot when a $\overline{\mathrm{SNR}}>3$ cut is applied to the unweighted analysis, the ”Distribution” p-value is improved, but there still exists a wavenumber correlation. However, we note that the ratio is very close to unity with a maximal value near $2$ %. The NOW discrepancy can therefore be ascribed to noisy spectra with $1<\overline{\mathrm{SNR}}<3$, and probably mostly to spectra with $\overline{\mathrm{SNR}}\gtrsim 1$. These spectra have a very small weight and a negligible contribution to the baseline power spectrum, which, therefore, is hardly affected by removing $\overline{\mathrm{SNR}}<3$ spectra, as illustrated by the lower plot, which shows a small (near $1$ %) difference, which is well below the level of small-scale systematics (e.g., resolution) associated to the final measurement in section 6.2. Even if the impact of biasing caused by low $\overline{\mathrm{SNR}}$ sub-forest is mostly accounted for by our $\overline{\mathrm{SNR}}$-weighting, the $\overline{\mathrm{SNR}}>3$ can be considered a potential variation for cosmological interpretations.

Each quasar is planned to be observed several times by the DESI instrument, and since DR1 is an intermediate data set, it contains spectra observed with very different numbers of exposures. The coadding of exposures can potentially introduce a bias, especially when it involves a large number of exposures with very different integration times and observing conditions. A data split involving quasars observed over only one (N=1), two or three (N=2$|$3), or four or more (N$>$=4) exposures is shown in the left panels of figure 7. The data splits contain, on average over redshift bins, $39$ %, $36$ %, and $25$ % of the total number of sub-forests, respectively, thus resulting in larger error bars. The ratios do not show strong variations from unity; however, only the (N=2 $|$3) and (N $>$=4) are passing the $\chi^{2}$ test. The low number of exposures (N=1) can potentially bias our measurement, especially at small scales, because a quasar with a low number of exposures tends to be noisier, thus giving different weights when averaging. Removing this sample implies a too-large decrease in statistics, but this variation should be considered for cosmological interpretation and future measurement with larger data sets and a larger number of exposures.

The focal plane of the DESI instrument is composed of ten modules called petals, each containing 500 robotic fiber positioners. The light measured by those positioners is linked to ten spectrographs by optical fibers. Each spectrograph receives the fiber associated with one petal and scatters the received light using volume phase holographic gratings. Three CCD cameras per spectrograph read the spectra, covering different wavelength ranges. Finally, four amplifiers per CCD are positioned in squares to perform the readout. In the $x$ direction, a CCD image contains different spectra chunks, while $y$ is the direction of wavelengths. A given quasar spectrum is then measured over six amplifiers (2 per CCD camera). For one spectrograph, the fibers directed to the left and right parts of the CCD cameras constitute two sets, measured with completely distinct amplifiers. At the spectrum level, the data can be split between 20 amplifier sets, according to the left and right parts of the ten spectrographs.

We measured $P_{1\mathrm{D},\alpha}$ over those 20 amplifier set splits independently, which constitutes a significant drop in terms of statistics. The right panels of figure 7 show the ratio between baseline and the power spectrum using only spectra from the amplifier set number 18 with constant-weight averaging (18NOW), with $\overline{\mathrm{SNR}}$-weight averaging (18W), and with the $\overline{\mathrm{SNR}}$-weight averaging but using the continuum from the complete DR1 data set (18WCONT). The 18^th amplifier set was chosen because it shows a considerable difference for low redshift near wavenumber $k\sim 0.65$ Å^-1. Among the 20 splits of the amplifier set, only two exhibit this behavior, with the 16^th data split showing a smaller spike amplitude. Those spikes originate from oscillatory patterns of period $\Delta\lambda=9.6$ Å(corresponding to $k\sim 0.65$ Å^-1) present on DR1 calibrations images used for bias corrections.

The source of those patterns was not identified and was highly complicated to correct due to their partial coverage of the amplifier and sporadic occurrence. In contrast to the companion paper K25, the impact of this feature drastically decreases when using the $\overline{\mathrm{SNR}}$ weighting scheme. Furthermore, comparing 18W and 18WCONT indicates that a subsequent part of the spike is caused by the continuum fitting procedure itself when using only the 18^th amplifier set. The spike is more prominent at small redshifts. As pointed out in section 4.2, the cross-exposure estimator accounts for this spike because it uses the individual cross-power spectrum between different fibers and, consequently, different amplifier sets. As we use the cross-exposure estimator for the lower redshift ($z<2.7$) where the $0.65$ Å^-1 spike is the largest ($2.3\sigma$ for $z=2.2$ and $3.5\sigma$ for $z=2.4$ and dropping to $1\sigma$ or below for larger redshifts), we consider that this feature does not impact our measurement, and we choose not to add a dedicated systematic uncertainty. We interpret the low $p$-value for the last panel to be mostly caused by fluctuations and not by the $0.65$ Å^-1 spike itself.

Different imaging surveys are used for the quasar target selection of the DESI survey [50, 51, 52]. Inhomogeneities in the photometric target pre-selection can significantly impact galaxy clustering studies, since they affect the target density and subsequently the estimation of the matter density. We do not expect an impact on the measurement of $P_{1\mathrm{D},\alpha}$ since each spectrum pixel is a proxy for the IGM density independently of the quasar density. However, as a sanity check, we still perform a data split between the North Galactic Cap (NGC) and the South Galactic Cap (SGC), as shown on the left panels of figure 8. The NGC dataset contains $65$ % of the sub-forests, and the SGC $35$ %. As expected, the $P_{1\mathrm{D},\alpha}$ measured on both data splits show no visible difference with the baseline, as also indicated by the reported $p$-values.

The last data split performed follows the BAO analysis [29, 31] and is related to the Civ equivalent width (EW) in angstrom, i.e. the spectral area associated with the quasar Civ emission line. The EW variable is a good quasar population discriminant because of the Baldwin effect, which causes EW to be anti-correlated to the quasar continuum luminosity [64]. Consequently, the high EW quasars will tend to be noisier. We use the EW measurement performed with fastspecfit\faGithub⁹⁹9https://github.com/desihub/fastspecfit in [29, 31] to create two populations equal in terms of quasar number: EW$>$39 and EW$<$39. The result of $P_{1\mathrm{D},\alpha}$ computations for those two EW populations is shown in the right panels of figure 8. Since we apply a $\overline{\mathrm{SNR}}>1$ quality cut, the noisier EW$>$39 measurement contains only $35$ % of the total sub-forests. The EW$>$39 data set shows a $5$ % discrepancy for small redshift and $k<0.2$ Å^-1. This discrepancy can be interpreted as a proximity effect impacting the larger scales of the quasar spectra. The EW$<$39 sample is also showing a discrepancy, shown by the $p$-value, but with a smaller level (near $2$ % at large scales). Removing the EW$>$39 sample would imply a statistical drop that is too large. For this analysis, we keep the two populations of quasars in the $P_{1\mathrm{D},\alpha}$ baseline measurement. However, the EW$<$39 data sample should also be considered as a variation for cosmological analysis, and this EW test should be investigated more in future studies. For now, the level of the $\sim 1\%$ induced bias compared to EW$<$39 is well below the level of large-scale systematics such as DLA completeness detailed in section 6.2.

We performed an extensive data split study to find potential biases in our $P_{1\mathrm{D},\alpha}$ measurement. Those data splits allow us to validate our measurement and improve the characterization of systematic uncertainties. Notably, this study indicates that the impact of $AI$ BAL on the $P_{1\mathrm{D},\alpha}$ measurement cannot be ignored, and we changed the baseline measurement to mask the associated absorptions. This analysis shows that some effects are slightly changing the global trend of $P_{1\mathrm{D},\alpha}$, and some imprint fluctuations sufficient not to pass our statistical tests. We can also note that the ”Correlation $p$-value” test is especially severe. Removing the sub-samples that do not pass the tests would imply a significant statistical drop; we then keep them for this paper. First, those effects (continuum fitting variations, number of exposures, amplifier splits, and Civ EW quasar population) should be considered as baseline variations for any cosmological interpretation. Secondly, those sub-samples should potentially be removed for future DESI releases with even larger data sets available.

As the aim of this DR1 measurement is to be interpreted in terms of cosmological parameters, we estimate the covariance matrix associated with $P_{1\mathrm{D},\alpha}$. Furthermore, we improve the study conducted in R23 concerning the estimation of statistical and systematic uncertainties. The full derivation of the variance and covariance estimators used here is presented in Appendix B.

The expression of the power spectrum $P_{1\mathrm{D},\alpha}(A,z)$ involves a weighted average of $P_{s}(k)$ (Eq. 3.6). In Appendix B, we derive an estimator for the variance of a weighted average $\sum w_{i}A_{i}/\sum w_{i}$ and for the covariance between two such weighted averages. However, this derivation requires that the $A_{i}$ be uncorrelated. When a sub-forest $s$ has several wavenumbers $k$ in a given $A$ bin, these wavenumbers are correlated, so we cannot use our variance estimator for the weighted average $\langle P_{s}(k)\rangle_{s\in z,k\in A}$.

We therefore introduce $A_{s}$ as the ensemble of wavenumbers associated with sub-forest $s$ that fall into bin $A$, $N_{A,s}$ the number of wavenumbers in $A_{s}$, and $P_{A,s}=\sum_{k\in A_{s}}P_{s}(k)/N_{A,s}$ the average of individual power spectra in bin $A$. We can neglect the correlation between different sub-forests as accounting for them would give a three-dimensional estimator. Therefore, we use our variance estimator for the $P_{A,s}$. We introduce the weights $w_{A,s}^{\prime}=N_{A,s}w_{A,s}$ to rewrite $\langle P_{s}(k)\rangle_{s\in z,k\in A}$ as a weighted average of the $P_{A,s}$:

We can then use our estimator for the variance

We note that our estimator fluctuates significantly and can even result in a negative variance for a small dataset size (at high redshift in our case). This issue is solved by smoothing the variance with a Savitzky-Golay filter implemented in scipy [65]. Similarly, the covariance matrix is smoothed along its two dimensions using a Savitzky-Golay filter implemented in sgolay2 \faGithub¹⁰¹⁰10https://github.com/espdev/sgolay2. The covariance and variance estimators have been validated in the companion KR25 paper by comparing the variation between $20$ DR1-like uncontaminated mocks with the average of those mocks. This validation showed that our estimator slightly overestimates the covariance by $10$ %. We correct our covariance matrix by renormalizing by this factor, independently of redshift.

Figure 9 shows the standard deviation calculated with equation 6.2 for the DR1 power spectrum. The three smallest redshift bins ($2.2$, $2.4$, and $2.6$) have a different profile because the cross-exposure estimator (see section 4.2) is used for these bins. Due to the reduced number of sub-forests available for the cross-exposure estimator, the error bars increase, especially at large wavenumber. Compared to R23 results, the uncertainties are smoother due to the applied smoothing but have a similar profile. The level of uncertainties is much lower than for R23: between $2.6$ and $3$ times decrease depending on redshift. This improvement is even higher than the five times increase in quasar statistics, which should decrease the error bar only by $\sim 2.23$ times. Considering that $\overline{\mathrm{SNR}}$ weighting is used in both measurements, we interpret this significant improvement as caused by decreasing noise in the input data set.

Figure 10 shows the correlation matrix for the DR1 data set. The correlation profile is very similar to the one obtained for eBOSS in [45], with a significant correlation (near $20$ %) at small wavenumber and decreasing as a function of redshift. There are some instabilities for high redshift caused by the low number of available sub-forests. The cosmological interpretation of this measurement should account for those instabilities.

We re-evaluated the systematic uncertainty budget of R23 to account for new sources of systematics and improvements developed in this article and in the validation companion paper KR25. More specifically, pixel masking and resolution corrections were reevaluated to match DR1 statistics. The continuum fitting correction, see below, was also derived in the same paper.

All systematic uncertainties included in the total error budget are shown in figure 11. As mentioned in section 4.2, with the development of a cross-exposure estimator that measures the noise in a fully data-driven way, we consider that the systematic error associated with the noise level estimation is negligible and can safely be removed from the systematic error budget. Other systematic uncertainties associated to a bin $(A,z)$ and noted $\sigma_{\mathrm{syst,X}}$, are defined similarly to R23 in the following way:

• •

  Resolution damping: We ascribe an uncertainty to the stability of the Point Spread Function (PSF) and its impact on the DESI resolution matrix. Following R23, the average resolution damping in figure 15 is fitted with a simplified resolution damping model $\exp\left(-0.5(k\Delta\lambda)^{2}\right)\cdot\operatorname{sinc}\left(0.5k\Delta\lambda_{\text{pix }}\right)$ to measure the effective spectral resolution $\Delta\lambda$ in Å. The $P_{1\mathrm{D},\alpha}$ uncertainty is given by $\sigma_{\mathrm{syst,res}}(A,z)=2k_{A}^{2}\Delta\lambda\sigma_{\Delta\lambda}\cdot P_{1\mathrm{D},\alpha}(A,z)$, where $\sigma_{\Delta\lambda}$ is the error on $\Delta\lambda$ for which we take a conservative absolute error of $1\%$ following the PSF stability measurement in [61].

• •

  Resolution correction: A correction for the resolution modeling $A_{\mathrm{res}}(A,z)$ is derived in the companion paper KR25. We apply this multiplicative correction to the main $P_{1\mathrm{D},\alpha}$ measurement and ascribe a systematic error that is $30$ % of the correction, i.e., $\sigma_{\mathrm{syst,res}}(A,z)=0.3\left|A_{\mathrm{res}}(A,z)-1\right|P_{1\mathrm{D},\alpha}(A,z)$.

• •

  Metal power spectrum: The metal power spectrum is measured in section 4.1 in the SB1 side-band wavelength range and shown in figure 16. The $P_{1\mathrm{D},\alpha}$ measurement is directly corrected by a physically motivated fit $P_{\mathrm{SB1,m,z}}$. We add a systematic uncertainty equal to the statistical error bar of the side-band power spectrum.

• •

  Continuum fitting: A multiplicative correction accounting for the bias introduced by the continuum fitting is derived from DR1-like mocks in KR25. This correction, noted $A_{\mathrm{cont}}(A,z)$, directly multiplies the power spectrum measurement, and a $30$ % relative systematic uncertainty is added to the total budget.

• •

  Pixel masking corrections: At the continuum fitting stage, some spectrum pixels are masked due to DLA, BAL passing the $AI>0$ criterion, and atmospheric emission lines, as discussed in sections 2 and 5.1. Multiplicative corrections, respectively noted $A_{\mathrm{dla}}(z)$, $A_{\mathrm{bal}}(A,z)$, and $A_{\mathrm{line}}(A,z)$ are computed in KR25. As for the previous multiplicative corrections, we correct the power spectrum measurement and associate a $30$ % relative systematic uncertainty. For the special case of BAL masking, a dedicated study in KR25 shows that a 30 $\%$ value is over-estimating the systematic uncertainty, and we take a 6 $\%$ value.

• •

  Catalog completeness: The DLA and BAL finding algorithms detailed in section 2 result in incomplete detections, especially for the low $\overline{\mathrm{SNR}}$ spectra. The data splits presented in section 5.1 provide insight into the impact of undetected DLA and BAL. In contrast with R23, we derive the impact of DLA and BAL from the data splits instead of mocks, making this analysis independent of the distribution of DLA and BAL in mocks. For DLA, we compare the power spectrum obtained when DLA are masked and the power spectrum corrected by $A_{\mathrm{dla}}(z)$ to the case where no DLA are masked. The resulting ratio corresponds to the top left panel of figure 4 with the additional DLA masking correction. We fit this ratio with a power-law function provided by [63]. Following R23 which uses the same DLA finding algorithms, we ascribe a conservative $20$ % incompleteness to the DLA finders. The associated systematic uncertainty is obtained by multiplying this incompleteness factor by the product of the fitted function times $P_{1\mathrm{D},\alpha}$. In a very similar way, we estimate the impact of unmasked $AI$ BAL absorptions by comparing our new baseline for which we masked $AI$ BAL, corrected by $A_{\mathrm{bal}}(A,z)$, to the previous baseline measurement where those absortions where left unmasked. We use the same fitting function and associate a systematic uncertainty considering a $15$ % incompleteness for the finder, following [56].

For multiplicative correction, the $30$ % relative systematic uncertainty is motivated by considering a random shift of the correction between $0$ and $100$ %. Considering a uniform distribution for this shift gives a $30$ % standard deviation. It is a straightforward but very conservative way of generating an error bar associated with the multiplicative corrections because it gives overestimated error bars related to larger corrections, such as atmospheric lines. In contrast, taking the standard deviation of the correction for different mock realizations gives too small and non-conservative error bars. We chose not to use this method based on mocks for systematic uncertainties.

All systematic uncertainties are added in quadrature to compute the total systematic error, which is shown in the last panels of figure 11. The correlations between wavenumber bins induced by these systematics are complex to model, and their contribution to the covariance matrix will be discussed in the upcoming inference of this work. Given the improvement in statistics compared to R23, the statistical error bars have decreased by between 260 to 300 % depending on redshift, and now the $P_{1\mathrm{D},\alpha}$ measurement is dominated by systematic errors for most redshifts and wavenumbers.