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 data sets to better quantify the performance of the DLA finder and its impact on our results.

We found that the resulting $P_{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 $P_{1D}$ results are 5–15 per cent 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 per cent, and QMLE is well suited for DESI $P_{1D}$ measurements. The next five years will bring incredible accuracy, precision, and power to Ly$\alpha$ forest cosmology.

## ACKNOWLEDGEMENTS

This material is based upon work supported by the U.S. Department of Energy (DOE), Office of Science, Office of High-Energy Physics, under Contract No. DE-AC02-05CH11231, and by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract. Additional support for DESI was provided by the U.S. National Science Foundation (NSF), Division of Astronomical Sciences under Contract No. AST-0950945 to the NSF’s National Optical-Infrared Astronomy Research Laboratory; the Science and Technology Facilities Council of the United Kingdom; the Gordon and Betty Moore Foundation; the Heising-Simons Foundation; the French Alternative Energies and Atomic Energy Commission (CEA); the National Council of Science and Technology of Mexico (CONACYT); the Ministry of Science and Innovation of Spain (MICINN), and by the DESI Member Institutions: https://www.desi.lbl.gov/collaborating-institutions. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the U.S. National Science Foundation, the U.S. Department of Energy, or any of the listed funding agencies.

The authors are honoured to be permitted to conduct scientific research on Iolkam Du'ag (Kitt Peak), a mountain with particular significance to the Tohono O'odham Nation.

## DATA AVAILABILITY

All data points shown in the figures are available in simple text files on the following website: https://doi.org/10.5281/zenodo.8007370. Some underlying DESI spectra will be publicly available as early data release (EDR) in 2023. We added two months of main surveys to improve our statistics. The main survey spectra will be made publicly available as part of Year 1 DR in the future.

## REFERENCES

Alexander D. M. et al., 2023, AJ, 165, 124

Astropy Collaboration, 2013, A&A, 558, A33

Astropy Collaboration, 2018, AJ, 156, 123

Astropy Collaboration, 2022, ApJ, 935, 167

Baur J., Palanque-Delabrouille N., Yèche C., Magneville C., Viel M., 2016, J. Cosmol. Astropart. Phys., 2016, 012

Bautista J. E. et al., 2017, A&A, 603, A12

Bi H., Davidsen A. F., 1997, ApJ, 479, 523

Bi H. G., Boerner G., Chu Y., 1992, A&A, 266, 1

Boera E., Becker G. D., Bolton J. S., Nasir F., 2019, ApJ, 872, 101

Bolton A. S., Schlegel D. J., 2010, PASP, 122, 248

Boyarsky A., Lesgourgues J., Ruchayskiy O., Viel M., 2009, J. Cosmol. Astropart. Phys., 2009, 012

Busca N., Balland C., 2018, preprint (arXiv:1808.09955)

Cabayol-Garcia L., Chaves-Montero J., Font-Ribera A., Pedersen C., 2023, MNRAS, 525, 3499

Cen R., Miralda-Escudé J., Ostriker J. P., Rauch M., 1994, ApJ, 437, L9

Chabanier S. et al., 2019, J. Cosmol. Astropart. Phys., 2019, 017

Chabanier S. et al., 2023, MNRAS, 518, 3754

Chaussidon E. et al., 2023, ApJ, 944, 107

Croft R. A. C., Weinberg D. H., Katz N., Hernquist L., 1998, ApJ, 495, 44

Croft R. A. C., Hu W., Davé R., 1999, Phys. Rev. Lett., 83, 1092

DESI Collaboration, 2016a, preprint (arXiv:1611.00036)

DESI Collaboration, 2016b, preprint (arXiv:1611.00037)

DESI Collaboration, 2022, AJ, 164, 207

DESI Collaboration, 2023a, preprint (arXiv:2306.06307)

DESI Collaboration, 2023b, preprint (arXiv:2306.06308)

Dalcin L., Fang Y.-L. L., 2021, Comput. Sci. Eng., 23, 47

Dawson K. S. et al., 2016, AJ, 151, 44

Dembinski H. et al., 2022, Zenodo, p. https://doi.org/10.5281/zenodo.6998452

Dey A. et al., 2019, AJ, 157, 168

Dierckx P., 1981, IMA J. Num. Anal., 1, 267

Dierckx P., 1995, Curve and Surface Fitting with Splines. Monographs on Numerical Analysis, Oxford University Press, New York, https://books.google.com/books?id = -RIQ3SR0sZMC

Ennesser L., Martini P., Font-Ribera A., Pérez-Ràfols I., 2022, MNRAS, 511, 3514

Farr J., Font-Ribera A., Pontzen A., 2020, J. Cosmol. Astropart. Phys., 2020, 015

Filbert S. et al., 2023, preprint (arXiv:2309.03434)

Font-Ribera A., McDonald P., Slosar A., 2018, J. Cosmol. Astropart. Phys., 2018, 003

Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

Frigo M., Johnson S. G., 2005, Proc. IEEE, 93, 216
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 data sets to better quantify the performance of the DLA finder and its impact on our results.

We found that the resulting $P_{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 $P_{1D}$ results are 5–15 per cent 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 per cent, and QMLE is well suited for DESI $P_{1D}$ measurements. The next five years will bring incredible accuracy, precision, and power to Ly$\alpha$ forest cosmology.

## ACKNOWLEDGEMENTS

This material is based upon work supported by the U.S. Department of Energy (DOE), Office of Science, Office of High-Energy Physics, under Contract No. DE-AC02-05CH11231, and by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract. Additional support for DESI was provided by the U.S. National Science Foundation (NSF), Division of Astronomical Sciences under Contract No. AST-0950945 to the NSF's National Optical-Infrared Astronomy Research Laboratory; the Science and Technology Facilities Council of the United Kingdom; the Gordon and Betty Moore Foundation; the Heising-Simons Foundation; the French Alternative Energies and Atomic Energy Commission (CEA); the National Council of Science and Technology of Mexico (CONACYT); the Ministry of Science and Innovation of Spain (MICINN), and by the DESI Member Institutions: https://www.desi.lbl.gov/collaborating-institutions. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the U.S. National Science Foundation, the U.S. Department of Energy, or any of the listed funding agencies.

The authors are honoured to be permitted to conduct scientific research on Iolkam Du'ag (Kitt Peak), a mountain with particular significance to the Tohono O'odham Nation.

## DATA AVAILABILITY

All data points shown in the figures are available in simple text files on the following website: https://doi.org/10.5281/zenodo.8007370. Some underlying DESI spectra will be publicly available as early data release (EDR) in 2023. We added two months of main surveys to improve our statistics. The main survey spectra will be made publicly available as part of Year 1 DR in the future.

## REFERENCES

Alexander D. M. et al., 2023, AJ, 165, 124

Astropy Collaboration, 2013, A&A, 558, A33

Astropy Collaboration, 2018, AJ, 156, 123

Astropy Collaboration, 2022, ApJ, 935, 167

Baur J., Palanque-Delabrouille N., Yèche C., Magneville C., Viel M., 2016, J. Cosmol. Astropart. Phys., 2016, 012

Bautista J. E. et al., 2017, A&A, 603, A12

Bi H., Davidsen A. F., 1997, ApJ, 479, 523

Bi H. G., Boerner G., Chu Y., 1992, A&A, 266, 1

Boera E., Becker G. D., Bolton J. S., Nasir F., 2019, ApJ, 872, 101

Bolton A. S., Schlegel D. J., 2010, PASP, 122, 248

Boyarsky A., Lesgourgues J., Ruchayskiy O., Viel M., 2009, J. Cosmol. Astropart. Phys., 2009, 012

Busca N., Balland C., 2018, preprint (arXiv:1808.09955)

Cabayol-Garcia L., Chaves-Montero J., Font-Ribera A., Pedersen C., 2023, MNRAS, 525, 3499

Cen R., Miralda-Escudé J., Ostriker J. P., Rauch M., 1994, ApJ, 437, L9

Chabanier S. et al., 2019, J. Cosmol. Astropart. Phys., 2019, 017

Chabanier S. et al., 2023, MNRAS, 518, 3754

Chaussidon E. et al., 2023, ApJ, 944, 107

Croft R. A. C., Weinberg D. H., Katz N., Hernquist L., 1998, ApJ, 495, 44

Croft R. A. C., Hu W., Davé R., 1999, Phys. Rev. Lett., 83, 1092

DESI Collaboration, 2016a, preprint (arXiv:1611.00036)

DESI Collaboration, 2016b, preprint (arXiv:1611.00037)

DESI Collaboration, 2022, AJ, 164, 207

DESI Collaboration, 2023a, preprint (arXiv:2306.06307)

DESI Collaboration, 2023b, preprint (arXiv:2306.06308)

Dalcin L., Fang Y.-L. L., 2021, Comput. Sci. Eng., 23, 47

Dawson K. S. et al., 2016, AJ, 151, 44

Dembinski H. et al., 2022, Zenodo, p. https://doi.org/10.5281/zenodo.6998452

Dey A. et al., 2019, AJ, 157, 168

Dierckx P., 1981, IMA J. Num. Anal., 1, 267

Dierckx P., 1995, Curve and Surface Fitting with Splines. Monographs on Numerical Analysis, Oxford University Press, New York, https://books.google.com/books?id = -RIQ3SR0sZMC

Ennesser L., Martini P., Font-Ribera A., Pérez-Ràfols I., 2022, MNRAS, 511, 3514

Farr J., Font-Ribera A., Pontzen A., 2020, J. Cosmol. Astropart. Phys., 2020, 015

Filbert S. et al., 2023, preprint (arXiv:2309.03434)

Font-Ribera A., McDonald P., Slosar A., 2018, J. Cosmol. Astropart. Phys., 2018, 003

Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

Frigo M., Johnson S. G., 2005, Proc. IEEE, 93, 216
Galassi M. et al., 2021, GNU Scientific Library: Reference manual (release 2.7), [Online]. Available: https://www.gnu.org/software/gsl/

García L. Á. et al., 2023, MNRAS, 526, 4848

Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759

Gunn J. E., Peterson B. A., 1965, ApJ, 142, 1633

Guo Z., Martini P., 2019, ApJ, 879, 72

Guy J. et al., 2023, AJ, 165, 144

Hamilton A. J. S., 1997, MNRAS, 289, 285

Harris C. R. et al., 2020, Nature, 585, 357

Herrera-Alcantar H. K. et al., 2023, preprint (arXiv:2401.00303)

Ho M.-F., Bird S., Garnett R., 2021, MNRAS, 507, 704

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Iršič V. et al., 2017, MNRAS, 466, 4332

Iršič V. et al., 2017a, Phys. Rev. D, 96, 023522

Iršič V., Viel M., Haehnelt M. G., Bolton J. S., Becker G. D., 2017b, Phys. Rev. Lett., 119, 031302

Karaçaylı N. G., Font-Ribera A., Padmanabhan N., 2020, MNRAS, 497, 4742

Karaçaylı N. G. et al., 2022, MNRAS, 509, 2842

Karaçaylı N. G. et al., 2023, MNRAS, 522, 5980

Kirkby D. et al., 2021, desihub/specsim: August 2020 Release, Zenodo,

Kramida A., Ralchenko Yu., Reader J., NIST ASD Team, 2021, NIST Atomic Spectra Database (ver. 5.9), [Online]. Available: https://physics.nist.gov/asd[2022, September 1]. National Institute of Standards and Technology, Gaithersburg, MD.

Levi M. et al., 2013, preprint (arXiv:1308.0847)

Lewis A., 2019, preprint (arXiv:1910.13970)

Lukić Z., Stark C. W., Nugent P., White M., Meiksin A. A., Almgren A., 2015, MNRAS, 446, 3697

Lynds R., 1971, ApJ, 164, L73

Mas-Ribas L., Mauland R., 2019, ApJ, 886, 151

McDonald P., Seljak U., Cen R., Bode P., Ostriker J. P., 2005, MNRAS, 360, 1471

McDonald P. et al., 2006, ApJS, 163, 80

Myers A. D. et al., 2023, AJ, 165, 50

Napolitano L. et al., 2023, AJ, 166, 99

Narayanan V. K., Spergel D. N., Davé R., Ma C.-P., 2000, ApJ, 543, L103

Padmanabhan N., White M., Zhou H. H., O'Connell R., 2016, MNRAS, 460, 1567

Palanque-Delabrouille N. et al., 2013, A&A, 559, A85

Palanque-Delabrouille N. et al., 2015, J. Cosmol. Astropart. Phys., 2015, 045

Pedersen C., Font-Ribera A., Rogers K. K., McDonald P., Peiris H. V., Pontzen A., Slosar A., 2021, J. Cosmol. Astropart. Phys., 2021, 033

Ravoux C. et al., 2023, MNRAS, 526, 5118

Rogers K. K., Bird S., Peiris H. V., Pontzen A., Font-Ribera A., Leistedt B., 2018, MNRAS, 474, 3032

Schmidt M., 1965, ApJ, 141, 1295

Seljak U., 1998, ApJ, 506, 64

Seljak U., Makarov A., McDonald P., Trac H., 2006, Phys. Rev. Lett., 97, 191303

Sherman J., Morrison W. J., 1950, Ann. Math. Stat., 21, 124

Silber J. H. et al., 2023, AJ, 165, 9

du Mas des Bourboux H. et al., 2019, ApJ, 878, 47

Slosar A. et al., 2011, J. Cosmol. Astropart. Phys., 2011, 001

Slosar A. et al., 2013, J. Cosmol. Astropart. Phys., 2013, 026

Tegmark M., Taylor A. N., Heavens A. F., 1997, ApJ, 480, 22

du Mas des Bourboux H. et al., 2020, ApJ, 901, 153

Tegmark M., Hamilton A. J. S., Strauss M. A., Vogeley M. S., Szalay A. S., 1998, ApJ, 499, 555

Viel M., Weller J., Haehnelt M. G., 2004, MNRAS, 355, L23

Viel M., Becker G. D., Bolton J. S., Haehnelt M. G., 2013, Phys. Rev. D, 88, 043502

Villasenor B., Robertson B., Madau P., Schneider E., 2022, ApJ, 933, 59

Villasenor B., Robertson B., Madau P., Schneider E., 2023, Phys. Rev. D, 108, 023502

Walther M., Hennawi J. F., Hiss H., Oñorbe J., Lee K.-G., Rorai A., O'Meara J., 2017, ApJ, 852, 22

Walther M., Oñorbe J., Hennawi J. F., Lukić Z., 2019, ApJ, 872, 13

Walther M., Armengaud E., Ravoux C., Palanque-Delabrouille N., Yèche C., Lukić Z., 2021, J. Cosmol. Astropart. Phys., 2021, 059

Wang M.-Y., Croft R. A. C., Peter A. H. G., Zentner A. R., Purcell C. W., 2013, Phys. Rev. D, 88, 123515

Wang B. et al., 2022, ApJS, 259, 28

Woodbury M. A., 1950, Memorandum Rept. 42, Statistical Research Group. Princeton Univ., Princeton, p. 4

Wright E. L. et al., 2010, AJ, 140, 1868

Yèche C., Palanque-Delabrouille N., Baur J., du Mas des Bourboux H., 2017, J. Cosmol. Astropart. Phys., 2017, 047

Yèche C. et al., 2020, Res. Notes Am. Astron. Soc., 4, 179

Zhang Y., Anninos P., Norman M. L., 1995, ApJ, 453, L57

Zonca A., Singer L., Lenz D., Reinecke M., Rosset C., Hivon E., Gorski K., 2019, J. Open Source Softw., 4, 1298

## APPENDIX A: CONTINUUM ERROR MARGINALIZATION

We fit the quasar continuum in the forest region as discussed previously. This method not only fits for the quasar continuum (and the mean transmission of the IGM), but also fits the large-scale density modes in the forest. This removes the large-scale information from the estimated flux transmission fluctuations and further cause distortions in the 3D correlation function (Slosar et al. 2011; Bautista et al. 2017; du Mas des Bourboux et al. 2020). Fortunately, the modes of error are known; and the QMLE is capable of marginalizing out these errors (Slosar et al. 2013, Appendix B).

The continuum marginalization in Karaçaylı, Font-Ribera & Padmanabhan (2020) is implemented by modifying the covariance matrix to $C' = C + Nvv^T$, where N is large and v is the mode to be marginalized out. Then, one can show that $C'^{-1}(\delta'_{F} + \alpha v) \approx C^{-1}\delta'_{F}$, where the new data vector $\delta'_{F}$ is orthogonal to v, which effectively removes any information from data that is in mode v.

We updated our continuum marginalization technique for more numerical stability. Instead of adding large numbers to the covariance matrix, which could destabilize the inversion, we take advantage of the Sherman-Morrison formula (Sherman & Morrison 1950):

$$ (C + vv^T)^{-1} = C^{-1} - \frac{C^{-1}vv^T C^{-1}}{1 + v^T C^{-1} v} \quad (A1) $$

Since the covariance matrix C is symmetric, this formula can be calculated by defining an intermediate vector $y = C^{-1}v$. Moreover, the marginalization vector mode v is theoretically multiplied by a large number, such that $v^T C^{-1} v \gg 1$, and that large number cancels out by the division, hence does not need to be explicitly defined anymore. Putting these together also makes the workings of the marginalization clearer:

$$ C^{-1} = C^{-1} - \frac{yy^T}{y^T v} \quad (A2) $$

$$ C^{-1}v = C^{-1}v - \frac{yy^T v}{y^T v} = 0 \quad (A3) $$

There are two caveats to this approach. First, when we marginalize over multiple templates, we have to iterate this formula with the updated covariance matrix in each step. Fortunately, this does not introduce additional memory strain in our case since we already store
Galassi M. et al., 2021, GNU Scientific Library: Reference manual (release 2.7), [Online]. Available: https://www.gnu.org/software/gsl/

García L. Á. et al., 2023, MNRAS, 526, 4848

Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759

Gunn J. E., Peterson B. A., 1965, ApJ, 142, 1633

Guo Z., Martini P., 2019, ApJ, 879, 72

Guy J. et al., 2023, AJ, 165, 144

Hamilton A. J. S., 1997, MNRAS, 289, 285

Harris C. R. et al., 2020, Nature, 585, 357

Herrera-Alcantar H. K. et al., 2023, preprint (arXiv:2401.00303)

Ho M.-F., Bird S., Garnett R., 2021, MNRAS, 507, 704

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Iršič V. et al., 2017, MNRAS, 466, 4332

Iršič V. et al., 2017a, Phys. Rev. D, 96, 023522

Iršič V., Viel M., Haehnelt M. G., Bolton J. S., Becker G. D., 2017b, Phys. Rev. Lett., 119, 031302

Karaçaylı N. G., Font-Ribera A., Padmanabhan N., 2020, MNRAS, 497, 4742

Karaçaylı N. G. et al., 2022, MNRAS, 509, 2842

Karaçaylı N. G. et al., 2023, MNRAS, 522, 5980

Kirkby D. et al., 2021, desihub/specsim: August 2020 Release, Zenodo,

Kramida A., Ralchenko Yu., Reader J., NIST ASD Team, 2021, NIST Atomic Spectra Database (ver. 5.9), [Online]. Available: https://physics.nist.gov/asd[2022, September 1]. National Institute of Standards and Technology, Gaithersburg, MD.

Levi M. et al., 2013, preprint (arXiv:1308.0847)

Lewis A., 2019, preprint (arXiv:1910.13970)

Lukić Z., Stark C. W., Nugent P., White M., Meiksin A. A., Almgren A., 2015, MNRAS, 446, 3697

Lynds R., 1971, ApJ, 164, L73

Mas-Ribas L., Mauland R., 2019, ApJ, 886, 151

McDonald P., Seljak U., Cen R., Bode P., Ostriker J. P., 2005, MNRAS, 360, 1471

McDonald P. et al., 2006, ApJS, 163, 80

Myers A. D. et al., 2023, AJ, 165, 50

Napolitano L. et al., 2023, AJ, 166, 99

Narayanan V. K., Spergel D. N., Davé R., Ma C.-P., 2000, ApJ, 543, L103

Padmanabhan N., White M., Zhou H. H., O'Connell R., 2016, MNRAS, 460, 1567

Palanque-Delabrouille N. et al., 2013, A&A, 559, A85

Palanque-Delabrouille N. et al., 2015, J. Cosmol. Astropart. Phys., 2015, 045

Pedersen C., Font-Ribera A., Rogers K. K., McDonald P., Peiris H. V., Pontzen A., Slosar A., 2021, J. Cosmol. Astropart. Phys., 2021, 033

Ravoux C. et al., 2023, MNRAS, 526, 5118

Rogers K. K., Bird S., Peiris H. V., Pontzen A., Font-Ribera A., Leistedt B., 2018, MNRAS, 474, 3032

Schmidt M., 1965, ApJ, 141, 1295

Seljak U., 1998, ApJ, 506, 64

Seljak U., Makarov A., McDonald P., Trac H., 2006, Phys. Rev. Lett., 97, 191303

Sherman J., Morrison W. J., 1950, Ann. Math. Stat., 21, 124

Silber J. H. et al., 2023, AJ, 165, 9

du Mas des Bourboux H. et al., 2019, ApJ, 878, 47

Slosar A. et al., 2011, J. Cosmol. Astropart. Phys., 2011, 001

Slosar A. et al., 2013, J. Cosmol. Astropart. Phys., 2013, 026

Tegmark M., Taylor A. N., Heavens A. F., 1997, ApJ, 480, 22

du Mas des Bourboux H. et al., 2020, ApJ, 901, 153

Tegmark M., Hamilton A. J. S., Strauss M. A., Vogeley M. S., Szalay A. S., 1998, ApJ, 499, 555

Viel M., Weller J., Haehnelt M. G., 2004, MNRAS, 355, L23

Viel M., Becker G. D., Bolton J. S., Haehnelt M. G., 2013, Phys. Rev. D, 88, 043502

Villasenor B., Robertson B., Madau P., Schneider E., 2022, ApJ, 933, 59

Villasenor B., Robertson B., Madau P., Schneider E., 2023, Phys. Rev. D, 108, 023502

Walther M., Hennawi J. F., Hiss H., Oñorbe J., Lee K.-G., Rorai A., O'Meara J., 2017, ApJ, 852, 22

Walther M., Oñorbe J., Hennawi J. F., Lukić Z., 2019, ApJ, 872, 13

Walther M., Armengaud E., Ravoux C., Palanque-Delabrouille N., Yèche C., Lukić Z., 2021, J. Cosmol. Astropart. Phys., 2021, 059

Wang M.-Y., Croft R. A. C., Peter A. H. G., Zentner A. R., Purcell C. W., 2013, Phys. Rev. D, 88, 123515

Wang B. et al., 2022, ApJS, 259, 28

Woodbury M. A., 1950, Memorandum Rept. 42, Statistical Research Group. Princeton Univ., Princeton, p. 4

Wright E. L. et al., 2010, AJ, 140, 1868

Yèche C., Palanque-Delabrouille N., Baur J., du Mas des Bourboux H., 2017, J. Cosmol. Astropart. Phys., 2017, 047

Yèche C. et al., 2020, Res. Notes Am. Astron. Soc., 4, 179

Zhang Y., Anninos P., Norman M. L., 1995, ApJ, 453, L57

Zonca A., Singer L., Lenz D., Reinecke M., Rosset C., Hivon E., Gorski K., 2019, J. Open Source Softw., 4, 1298

## APPENDIX A: CONTINUUM ERROR MARGINALIZATION

We fit the quasar continuum in the forest region as discussed previously. This method not only fits for the quasar continuum (and the mean transmission of the IGM), but also fits the large-scale density modes in the forest. This removes the large-scale information from the estimated flux transmission fluctuations and further cause distortions in the 3D correlation function (Slosar et al. 2011; Bautista et al. 2017; du Mas des Bourboux et al. 2020). Fortunately, the modes of error are known; and the QMLE is capable of marginalizing out these errors (Slosar et al. 2013, Appendix B).

The continuum marginalization in Karaçaylı, Font-Ribera & Padmanabhan (2020) is implemented by modifying the covariance matrix to $C' = C + Nvv^T$, where N is large and v is the mode to be marginalized out. Then, one can show that $C'^{-1}(\delta'_F + \alpha v) \approx C^{-1}\delta'_F$, where the new data vector $\delta'_F$ is orthogonal to v, which effectively removes any information from data that is in mode v.

We updated our continuum marginalization technique for more numerical stability. Instead of adding large numbers to the covariance matrix, which could destabilize the inversion, we take advantage of the Sherman-Morrison formula (Sherman & Morrison 1950):

$$ (C + vv^T)^{-1} = C^{-1} - \frac{C^{-1}vv^T C^{-1}}{1 + v^T C^{-1} v} \quad (A1) $$

Since the covariance matrix C is symmetric, this formula can be calculated by defining an intermediate vector $y = C^{-1}v$. Moreover, the marginalization vector mode v is theoretically multiplied by a large number, such that $v^T C^{-1} v \gg 1$, and that large number cancels out by the division, hence does not need to be explicitly defined anymore. Putting these together also makes the workings of the marginalization clearer:

$$ C^{-1} = C^{-1} - \frac{yy^T}{y^T v} \quad (A2) $$

$$ C^{-1}v = C^{-1}v - \frac{yy^T v}{y^T v} = 0 \quad (A3) $$

There are two caveats to this approach. First, when we marginalize over multiple templates, we have to iterate this formula with the updated covariance matrix in each step. Fortunately, this does not introduce additional memory strain in our case since we already store
the inverse covariance matrix in memory, but other applications may prefer the Woodbury formula (Woodbury 1950). Second, trying to remove two templates that are approximately the same results in a division by near zero. After removing the first vector $v_1$, the second vector results in $y_2 = C^{-1}v_2 \approx 0$. To prevent numerical problems, we first store all n template vectors in a rectangular matrix and perform a singular value decomposition. This yields an orthogonal basis for the templates and vectors with small singular values can be left out of the marginalization process.

## APPENDIX B: SIGNAL-NOISE COUPLING AND SMOOTHING

Another important caveat from our preliminary analysis on mocks is the well-known signal-noise coupling in the Ly$\alpha$ forest. The pipeline noise is coupled to the observed flux by construction through Poisson statistics, which means pixels with less flux are weighted more in the quadratic estimator formalism (McDonald et al. 2006). However, this does not mean the noise estimates are wrong, the noise is non-diagonal, or there are cross-correlations between signal and noise. The simulations truthfully populate the mocks with noise using independent random numbers. The core problem is that these (correct) noise amplitudes are biased weights for the inverse variance weighted averages since they depend upon the signal (see below for an analytical description of this effect). To mitigate this problem, the DESI pipeline determines a smoothed version of the sky-subtracted spectrum of each target that is obtained with a convolution using a Gaussian kernel of $\sigma = 10$Å (Guy et al. 2023). However, our preliminary analysis showed this is not enough to fully uncouple the signal and noise, so we further smooth out noise estimates in only their contribution to the covariance matrix in QMLE.

First, we find the median $M(x)$ and the median absolute deviation $D(x) = M(|x_i - M(x)|)$ in the pipeline noise estimates while ignoring pixels with $\sigma_{\text{pipe}} > 1000$. Noise outliers are then identified by

$$ \sigma_{\text{pipe}} > M(\sigma_{\text{pipe}}) + 3.5 \times D(\sigma_{\text{pipe}}), \quad (B1) $$

and replaced with the median pipeline noise value. We pad the noise array by 25 pixels at both ends with the edge values to mitigate 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 $t_{\text{exp}}$, and has the following contributions:

(i) Source electrons $N_{\text{source}} = C_q(\lambda)F(\lambda) \propto t_{\text{exp}}$. Note $F = \bar{F}(1+\delta)$

(ii) Sky electrons $N_{\text{sky}} \propto t_{\text{exp}}$.

(iii) Dark electrons $\propto t_{\text{exp}}$. We can absorb this term into sky contribution.

(iv) Read noise electrons as Gaussian noise with constant $\sigma_G^2$ independent of $t_{\text{exp}}$.

Then, pipeline variance estimate on flux $f$ is $\sigma_{f,\text{pipe}}^2 = N_{\text{source}} + N_{\text{sky}} + \sigma_G^2$. We can write the pipeline noise estimate on $F$ as follows:

$$ \sigma_{F, \text{pipe}}^2 = \frac{\sigma_{f, \text{pipe}}^2}{C_q^2(\lambda)} = \frac{F}{C_q} + \frac{N_{\text{sky}} + \sigma_G^2}{C_q^2}. \quad (B2) $$

<figure><img src="image_20.png" /><figcaption>**Figure B1.** Coupling parameter $\Lambda$ vs quasar brightness $C_q$ 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 $\Lambda = 34$ (low coupling).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure B2.** Biased power spectrum estimates for different $\Lambda$ 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 $\Lambda = 50$ (red line) as expected from the fiducial value.</figcaption></figure>

The total variance includes large-scale structure fluctuations $\sigma_{\text{LSS}}^2$, which should be multiplied by the mean IGM flux $\bar{F}$.

$$ \sigma_F^2 = \bar{F}^2 \sigma_{\text{LSS}}^2 + \frac{F}{C_q} + \frac{N_{\text{sky}} + \sigma_G^2}{C_q^2}. \quad (B3) $$

Let us write this variance on $F$ as follows:

$$ \sigma_F^2 = \frac{1}{C_q} \left[ F + C_q \bar{F}^2 \sigma_{\text{LSS}}^2 + \frac{N_{\text{sky}} + \sigma_G^2}{C_q} \right] \quad (B4) $$

$$ \sigma_F^2 = \frac{F + \Lambda}{C_q}, \quad (B5) $$

where we defined

$$ \Lambda = C_q \bar{F}^2 \sigma_{\text{LSS}}^2 + \frac{N_{\text{sky}} + \sigma_G^2}{C_q}. \quad (B6) $$

Using inverse variance weights means $w_i \propto (F + \Lambda)^{-1}$. Since $0 < F < 1$, our weights are highly correlated with the signal if $\Lambda \lesssim 1$.
the inverse covariance matrix in memory, but other applications may prefer the Woodbury formula (Woodbury 1950). Second, trying to remove two templates that are approximately the same results in a division by near zero. After removing the first vector $v_1$, the second vector results in $y_2 = C^{-1}v_2 \approx 0$. To prevent numerical problems, we first store all n template vectors in a rectangular matrix and perform a singular value decomposition. This yields an orthogonal basis for the templates and vectors with small singular values can be left out of the marginalization process.

## APPENDIX B: SIGNAL-NOISE COUPLING AND SMOOTHING

Another important caveat from our preliminary analysis on mocks is the well-known signal-noise coupling in the Ly$\alpha$ forest. The pipeline noise is coupled to the observed flux by construction through Poisson statistics, which means pixels with less flux are weighted more in the quadratic estimator formalism (McDonald et al. 2006). However, this does not mean the noise estimates are wrong, the noise is non-diagonal, or there are cross-correlations between signal and noise. The simulations truthfully populate the mocks with noise using independent random numbers. The core problem is that these (correct) noise amplitudes are biased weights for the inverse variance weighted averages since they depend upon the signal (see below for an analytical description of this effect). To mitigate this problem, the DESI pipeline determines a smoothed version of the sky-subtracted spectrum of each target that is obtained with a convolution using a Gaussian kernel of $\sigma = 10$Å (Guy et al. 2023). However, our preliminary analysis showed this is not enough to fully uncouple the signal and noise, so we further smooth out noise estimates in only their contribution to the covariance matrix in QMLE.

First, we find the median $M(x)$ and the median absolute deviation $D(x) = M(|x_i - M(x)|)$ in the pipeline noise estimates while ignoring pixels with $\sigma_{\text{pipe}} > 1000$. Noise outliers are then identified by

$$ \sigma_{\text{pipe}} > M(\sigma_{\text{pipe}}) + 3.5 \times D(\sigma_{\text{pipe}}), \quad (B1) $$

and replaced with the median pipeline noise value. We pad the noise array by 25 pixels at both ends with the edge values to mitigate 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 $t_{\text{exp}}$, and has the following contributions:

(i) Source electrons $N_{\text{source}} = C_q(\lambda)F(\lambda) \propto t_{\text{exp}}$. Note $F = \bar{F}(1 + \delta)$

(ii) Sky electrons $N_{\text{sky}} \propto t_{\text{exp}}$.

(iii) Dark electrons $\propto t_{\text{exp}}$. We can absorb this term into sky contribution.

(iv) Read noise electrons as Gaussian noise with constant $\sigma_G^2$ independent of $t_{\text{exp}}$.

Then, pipeline variance estimate on flux $f$ is $\sigma_{f,\text{pipe}}^2 = N_{\text{source}} + N_{\text{sky}} + \sigma_G^2$. We can write the pipeline noise estimate on $F$ as follows:

$$ \sigma_{F, \text{pipe}}^2 = \frac{\sigma_{f, \text{pipe}}^2}{C_q^2(\lambda)} = \frac{F}{C_q} + \frac{N_{\text{sky}} + \sigma_G^2}{C_q^2}. \quad (B2) $$

<figure><img src="image_20.png" /><figcaption>**Figure B1.** Coupling parameter $\Lambda$ vs quasar brightness $C_q$ 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 $\Lambda = 34$ (low coupling).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure B2.** Biased power spectrum estimates for different $\Lambda$ 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 $\Lambda = 50$ (red line) as expected from the fiducial value.</figcaption></figure>

The total variance includes large-scale structure fluctuations $\sigma_{\text{LSS}}^2$, which should be multiplied by the mean IGM flux $\bar{F}$.

$$ \sigma_F^2 = \bar{F}^2 \sigma_{\text{LSS}}^2 + \frac{F}{C_q} + \frac{N_{\text{sky}} + \sigma_G^2}{C_q^2}. \quad (B3) $$

Let us write this variance on $F$ as follows:

$$ \sigma_F^2 = \frac{1}{C_q} \left[ F + C_q \bar{F}^2 \sigma_{\text{LSS}}^2 + \frac{N_{\text{sky}} + \sigma_G^2}{C_q} \right] \quad (B4) $$

$$ \sigma_F^2 = \frac{F + \Lambda}{C_q}, \quad (B5) $$

where we defined

$$ \Lambda = C_q \bar{F}^2 \sigma_{\text{LSS}}^2 + \frac{N_{\text{sky}} + \sigma_G^2}{C_q}. \quad (B6) $$

Using inverse variance weights means $w_i \propto (F + \Lambda)^{-1}$. Since $0 < F < 1$, our weights are highly correlated with the signal if $\Lambda \lesssim 1$.
Unbiased limit $\Lambda \gg 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 $\text{Ly}\alpha$ forest. As we mentioned in the main text, DESI decouples signal and noise by smoothing source contribution at $10\AA$ 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. Fig. B1 shows coupling parameter $\Lambda$ vs quasar SNR. Low SNR quasars manifest less coupling (high $\Lambda$) 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 $\Lambda = 34$.

The weighted average estimator for the mean flux is given by $\hat{F} = \sum w_i F_i / \sum w_i$, where $w_i = (F_i + \Lambda)^{-1}$, where F is signal only as before. The expected value of this estimator can be calculated by using the probability distribution function $P(F)$.

$$ \langle \hat{F} \rangle = A \int_0^1 \frac{F}{F + \Lambda} P(F) dF \quad (B7) $$

$$ A^{-1} = \int_{0}^{1} \frac{1}{F + \Lambda} P(F) dF \quad (B8) $$

We can Taylor expand these expressions for $\Lambda \gg 1$. The normalization is $A \approx \Lambda (1 + \bar{F}/\Lambda)$. The mean flux estimate is

$$ \langle \hat{F} \rangle = A \int_0^1 \frac{F}{F + \Lambda} P(F) dF \quad (B9) $$

$$ \approx \frac{A}{\Lambda} \int F \left(1 - \frac{F}{\Lambda}\right) P(F) dF \quad (B10) $$

$$ = \frac{A}{\Lambda} \left( \bar{F} - \frac{\bar{F}^2}{\Lambda} \right) = \left( 1 + \frac{\bar{F}}{\Lambda} \right) \left( \bar{F} - \frac{\bar{F}^2}{\Lambda} \right) \quad (B11) $$

Finally, the error on the mean flux estimate is

$$ \Delta \bar{F} = - \frac{\bar{F}^2 \sigma_{\text{LSS}}^2}{\Lambda} \quad (B12) $$

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 $\delta$s. It is worth noting that using biased mean flux does not scale but shifts estimated deltas.

$$ \delta_i \rightarrow \frac{\bar{F}(1 + \delta_i)}{\bar{F} + \Delta \bar{F}} - 1 = \delta_i - \frac{\Delta \bar{F}}{\bar{F}} \quad (B13) $$

This does not make a big difference and only affects the $k=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

$$ \Delta \xi_{ij} = - \frac{\bar{F}}{\Lambda} \left( \langle \delta_i^2 \delta_j \rangle + \langle \delta_i \delta_j^2 \rangle \right) \quad (B14) $$

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, Fig. B2 shows only biased $P_{1D}$ with respect to different $\Lambda$ 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 $\Lambda = 50$ (red line) as expected from the fiducial value. As we note in the main text, smoothing the noise estimates solves this problem.

<footnote n="1">¹Center for Cosmology and AstroParticle Physics, The Ohio State University, 191 West Woodruff Avenue, Columbus, OH 43210, USA</footnote>

<footnote n="2">²Department of Astronomy, The Ohio State University, 4055 McPherson Laboratory, 140 W 18th Avenue, Columbus, OH 43210, USA</footnote>

<footnote n="3">³Department of Physics, The Ohio State University, 191 West Woodruff Avenue, Columbus, OH 43210, USA</footnote>

<footnote n="4">⁴Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA</footnote>

<footnote n="5">⁵Aix Marseille Univ, CNRS/IN2P3, CPPM, F-13288 Marseille, France</footnote>

<footnote n="6">⁶IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="7">⁷Excellence Cluster ORIGINS, Boltzmannstrasse 2, D-85748 Garching, Germany</footnote>

<footnote n="8">⁸Faculty of Physics, University Observatory, Ludwig-Maximilians-Universität, Scheinerstr. 1, D-81677 München, Germany</footnote>

<footnote n="9">⁹Physics Department, Boston University, 590 Commonwealth Avenue, Boston, MA 02215, USA</footnote>

<footnote n="10">¹⁰Departamento de Física, Universidad de Guanajuato - DCI, C.P. 37150 Leon, Guanajuato, México</footnote>

<footnote n="11">¹¹Department of Physics & Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="12">¹²Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, Bellaterra, E-08193 Barcelona, Spain</footnote>

<footnote n="13">¹³Department of Physics and Astronomy, The University of Utah, 115 South 1400 East, Salt Lake City, UT 84112, USA</footnote>

<footnote n="14">¹⁴Instituto de Física, Universidad Nacional Autónoma de México, Cd. de México, CP 04510, México</footnote>

<footnote n="15">¹⁵Departamento de Física, Universidad de los Andes, Cra. 1 No. 18A-10, Edificio Ip, CP 111711, Bogotá, Colombia</footnote>

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

<footnote n="17">¹⁷Consejo Nacional de Ciencia y Tecnología, Av. Insurgentes Sur 1582. Colonia Crédito Constructor, Del. Benito Juárez, CP 03940, México</footnote>

<footnote n="18">¹⁸Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="19">¹⁹Department of Physics, The University of Texas at Dallas, Richardson, TX 75080, USA</footnote>

<footnote n="20">²⁰Department of Physics, Southern Methodist University, 3215 Daniel Avenue, Dallas, TX 75275, USA</footnote>

<footnote n="21">²¹CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies (LPNHE), Sorbonne Université, F-75005 Paris, France</footnote>

<footnote n="22">²²NSF's NOIRLab, 950 N. Cherry Avenue, Tucson, AZ 85719, USA</footnote>

<footnote n="23">²³Institució Catalana de Recerca i Estudis Avançats, Passeig de Lluís Companys, 23, E-08010 Barcelona, Spain</footnote>

<footnote n="24">²⁴Department of Physics and Astronomy, Siena College, 515 Loudon Road, Loudonville, NY 12211, USA</footnote>

<footnote n="25">²⁵Institute of Cosmology & Gravitation, University of Portsmouth, Dennis Sciama Building, Portsmouth PO1 3FX, UK</footnote>

<footnote n="26">²⁶Department of Physics & Astronomy, University of Wyoming, 1000 E. University, Dept. 3905, Laramie, WY 82071, USA</footnote>

<footnote n="27">²⁷National Astronomical Observatories, Chinese Academy of Sciences, A20 Datun Road, Chaoyang District, Beijing 100012, P.R. China</footnote>

<footnote n="28">²⁸Instituto Avanzado de Cosmología A. C. San Marcos 11 - Atenas 202, Magdalena Contreras, 10720 Ciudad de México, México</footnote>

<footnote n="29">²⁹Department of Physics and Astronomy, University of Waterloo, 200 University Avenue W, Waterloo, ON N2L 3G1, Canada</footnote>

<footnote n="30">³⁰Perimeter Institute for Theoretical Physics, 31 Caroline St. North, Waterloo, ON N2L 2Y5, Canada</footnote>

<footnote n="31">³¹Waterloo Centre for Astrophysics, University of Waterloo, 200 University Avenue W, Waterloo, ON N2L 3G1, Canada</footnote>
Unbiased limit $\Lambda \gg 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 $\text{Ly}\alpha$ forest. As we mentioned in the main text, DESI decouples signal and noise by smoothing source contribution at $10\AA$ 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. Fig. B1 shows coupling parameter $\Lambda$ vs quasar SNR. Low SNR quasars manifest less coupling (high $\Lambda$) 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 $\Lambda = 34$.

The weighted average estimator for the mean flux is given by $\hat{F} = \sum w_i F_i / \sum w_i$, where $w_i = (F_i + \Lambda)^{-1}$, where F is signal only as before. The expected value of this estimator can be calculated by using the probability distribution function $P(F)$.

$$ \langle \hat{F} \rangle = A \int_0^1 \frac{F}{F + \Lambda} P(F) dF \qquad (B7) $$

$$ A^{-1} = \int_{0}^{1} \frac{1}{F + \Lambda} P(F) dF \qquad (B8) $$

We can Taylor expand these expressions for $\Lambda \gg 1$. The normalization is $A \approx \Lambda (1 + \bar{F}/\Lambda)$. The mean flux estimate is

$$ \langle \hat{F} \rangle = A \int_0^1 \frac{F}{F + \Lambda} P(F) dF \qquad (B9) $$

$$ \approx \frac{A}{\Lambda} \int F \left(1 - \frac{F}{\Lambda}\right) P(F) dF \qquad (B10) $$

$$ = \frac{A}{\Lambda} \left( \bar{F} - \frac{\bar{F}^2}{\Lambda} \right) = \left( 1 + \frac{\bar{F}}{\Lambda} \right) \left( \bar{F} - \frac{\bar{F}^2}{\Lambda} \right) \qquad (B11) $$

Finally, the error on the mean flux estimate is

$$ \Delta \bar{F} = - \frac{\bar{F}^2 \sigma_{\text{LSS}}^2}{\Lambda} \qquad (B12) $$

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 $\delta$s. It is worth noting that using biased mean flux does not scale but shifts estimated deltas.

$$ \delta_i \rightarrow \frac{\bar{F}(1 + \delta_i)}{\bar{F} + \Delta \bar{F}} - 1 = \delta_i - \frac{\Delta \bar{F}}{\bar{F}} \qquad (B13) $$

This does not make a big difference and only affects the $k=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

$$ \Delta \xi_{ij} = - \frac{\bar{F}}{\Lambda} \left( \langle \delta_i^2 \delta_j \rangle + \langle \delta_i \delta_j^2 \rangle \right) \qquad (B14) $$

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, Fig. B2 shows only biased $P_{1D}$ with respect to different $\Lambda$ 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 $\Lambda = 50$ (red line) as expected from the fiducial value. As we note in the main text, smoothing the noise estimates solves this problem.

<footnote n="1">¹Center for Cosmology and AstroParticle Physics, The Ohio State University, 191 West Woodruff Avenue, Columbus, OH 43210, USA</footnote>

<footnote n="2">²Department of Astronomy, The Ohio State University, 4055 McPherson Laboratory, 140 W 18th Avenue, Columbus, OH 43210, USA</footnote>

<footnote n="3">³Department of Physics, The Ohio State University, 191 West Woodruff Avenue, Columbus, OH 43210, USA</footnote>

<footnote n="4">⁴Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA</footnote>

<footnote n="5">⁵Aix Marseille Univ, CNRS/IN2P3, CPPM, F-13288 Marseille, France</footnote>

<footnote n="6">⁶IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="7">⁷Excellence Cluster ORIGINS, Boltzmannstrasse 2, D-85748 Garching, Germany</footnote>

<footnote n="8">⁸Faculty of Physics, University Observatory, Ludwig-Maximilians-Universität, Scheinerstr. 1, D-81677 München, Germany</footnote>

<footnote n="9">⁹Physics Department, Boston University, 590 Commonwealth Avenue, Boston, MA 02215, USA</footnote>

<footnote n="10">¹⁰Departamento de Física, Universidad de Guanajuato - DCI, C.P. 37150 Leon, Guanajuato, México</footnote>

<footnote n="11">¹¹Department of Physics & Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="12">¹²Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, Bellaterra, E-08193 Barcelona, Spain</footnote>

<footnote n="13">¹³Department of Physics and Astronomy, The University of Utah, 115 South 1400 East, Salt Lake City, UT 84112, USA</footnote>

<footnote n="14">¹⁴Instituto de Física, Universidad Nacional Autónoma de México, Cd. de México, CP 04510, México</footnote>

<footnote n="15">¹⁵Departamento de Física, Universidad de los Andes, Cra. 1 No. 18A-10, Edificio Ip, CP 111711, Bogotá, Colombia</footnote>

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

<footnote n="17">¹⁷Consejo Nacional de Ciencia y Tecnología, Av. Insurgentes Sur 1582. Colonia Crédito Constructor, Del. Benito Juárez, CP 03940, México</footnote>

<footnote n="18">¹⁸Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="19">¹⁹Department of Physics, The University of Texas at Dallas, Richardson, TX 75080, USA</footnote>

<footnote n="20">²⁰Department of Physics, Southern Methodist University, 3215 Daniel Avenue, Dallas, TX 75275, USA</footnote>

<footnote n="21">²¹CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies (LPNHE), Sorbonne Université, F-75005 Paris, France</footnote>

<footnote n="22">²²NSF's NOIRLab, 950 N. Cherry Avenue, Tucson, AZ 85719, USA</footnote>

<footnote n="23">²³Institució Catalana de Recerca i Estudis Avançats, Passeig de Lluís Companys, 23, E-08010 Barcelona, Spain</footnote>

<footnote n="24">²⁴Department of Physics and Astronomy, Siena College, 515 Loudon Road, Loudonville, NY 12211, USA</footnote>

<footnote n="25">²⁵Institute of Cosmology & Gravitation, University of Portsmouth, Dennis Sciama Building, Portsmouth PO1 3FX, UK</footnote>

<footnote n="26">²⁶Department of Physics & Astronomy, University of Wyoming, 1000 E. University, Dept. 3905, Laramie, WY 82071, USA</footnote>

<footnote n="27">²⁷National Astronomical Observatories, Chinese Academy of Sciences, A20 Datun Road, Chaoyang District, Beijing 100012, P.R. China</footnote>

<footnote n="28">²⁸Instituto Avanzado de Cosmología A. C. San Marcos 11 - Atenas 202, Magdalena Contreras, 10720 Ciudad de México, México</footnote>

<footnote n="29">²⁹Department of Physics and Astronomy, University of Waterloo, 200 University Avenue W, Waterloo, ON N2L 3G1, Canada</footnote>

<footnote n="30">³⁰Perimeter Institute for Theoretical Physics, 31 Caroline St. North, Waterloo, ON N2L 2Y5, Canada</footnote>

<footnote n="31">³¹Waterloo Centre for Astrophysics, University of Waterloo, 200 University Avenue W, Waterloo, ON N2L 3G1, Canada</footnote>
<footnote n="32">³²Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="33">³³Space Sciences Laboratory, University of California, Berkeley, 7 Gauss Way, Berkeley, CA 94720, USA</footnote>

<footnote n="34">³⁴Instituto de Astrofísica de Andalucía (CSIC), Glorieta de la Astronomía, s/n, E-18008 Granada, Spain</footnote>

<footnote n="35">³⁵Department of Physics and Astronomy, Sejong University, Seoul 143-747, Korea</footnote>

<footnote n="36">³⁶CIEMAT, Avenida Complutense 40, E-28040 Madrid, Spain</footnote>

<footnote n="37">³⁷Department of Physics & Astronomy, Ohio University, Athens, OH 45701, USA</footnote>

<footnote n="38">³⁸Instituto de Astrofísica de Canarias, C/Vía Láctea, s/n, E-38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="39">³⁹Dept. de Astrofísica, Universidad de La Laguna, E-38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="40">⁴⁰Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA</footnote>

<footnote n="41">⁴¹Department of Astronomy, Tsinghua University, 30 Shuangqing Road, Haidian District, Beijing 100190, China</footnote>

<footnote>This paper has been typeset from a *TEX/LATEX* file prepared by the author.</footnote>
<footnote n="32">³²Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="33">³³Space Sciences Laboratory, University of California, Berkeley, 7 Gauss Way, Berkeley, CA 94720, USA</footnote>

<footnote n="34">³⁴Instituto de Astrofísica de Andalucía (CSIC), Glorieta de la Astronomía, s/n, E-18008 Granada, Spain</footnote>

<footnote n="35">³⁵Department of Physics and Astronomy, Sejong University, Seoul 143-747, Korea</footnote>

<footnote n="36">³⁶CIEMAT, Avenida Complutense 40, E-28040 Madrid, Spain</footnote>

<footnote n="37">³⁷Department of Physics & Astronomy, Ohio University, Athens, OH 45701, USA</footnote>

<footnote n="38">³⁸Instituto de Astrofísica de Canarias, C/Vía Láctea, s/n, E-38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="39">³⁹Dept. de Astrofísica, Universidad de La Laguna, E-38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="40">⁴⁰Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA</footnote>

<footnote n="41">⁴¹Department of Astronomy, Tsinghua University, 30 Shuangqing Road, Haidian District, Beijing 100190, China</footnote>

<footnote>This paper has been typeset from a *TEX/LATEX* file prepared by the author.</footnote>