The 3D Lyman-𝛼 Forest Power Spectrum from eBOSS DR16

We measure the three-dimensional power spectrum (P3D) of the transmitted flux in the Lyman-𝛼 (Ly-𝛼 ) forest using the complete extended Baryon Oscillation Spectroscopic Survey data release 16 (eBOSS DR16). This sample consists of ∼ 205 , 000 quasar spectra in the redshift range 2 ≤ 𝑧 ≤ 4 at an effective redshift 𝑧 = 2 . 334. We propose a pair-count spectral estimator in configuration space, weighting each pair by exp ( 𝑖 k · r ) , for wave vector k and pixel pair separation r , effectively measuring the anisotropic power spectrum without the need for fast Fourier transforms. This accounts for the window matrix in a tractable way, avoiding artifacts found in Fourier-transform based power spectrum estimators due to the sparse sampling transverse to the line-of-sight of Ly-𝛼 skewers. We extensively test our pipeline on two sets of mocks: (i) idealized Gaussian random fields with a sparse sampling of Ly-𝛼 skewers, and (ii) log-normal LyaCoLoRe mocks including realistic noise levels, the eBOSS survey geometry and contaminants. On eBOSS DR16 data, the Kaiser formula with a non-linear correction term obtained from hydrodynamic simulations yields a good fit to the power spectrum data in the range ( 0 . 02 ≤ 𝑘 ≤ 0 . 35 ) ℎ Mpc − 1 at the 1 − 2 𝜎 level with a covariance matrix derived from LyaCoLoRe mocks. We demonstrate a promising new approach for full-shape cosmological analyses of Ly-𝛼 forest data from cosmological surveys such as eBOSS, the currently observing Dark Energy Spectroscopic Instrument and future surveys such as the Prime Focus Spectrograph, WEAVE-QSO and 4MOST.


INTRODUCTION
Light is absorbed by neutral hydrogen in the low-density, highly ionized intergalactic medium (IGM), producing a series of characteristic absorption features in quasar spectra, called the Lyman- (Ly-) forest.These absorption patterns are a measure of the distribution of neutral hydrogen along the line-of-sight.They represent a fundamental probe of large-scale structure at Mpc scales and below at high redshifts (2 < ∼  < ∼ 5) using ground-based observations.Given the sparseness of the Ly- forest 1 , the two-point correlation function (Slosar et al. 2011) and the one-dimensional power spectrum (P1D) along the line-of-sight (McDonald et al. 2006) have been widely adopted in cosmological analyses of Ly- data (e.g., tational cost of O ( 2 pix  3 0 ) (for number of spectral pixels  pix and maximum pair separation  0 ).
Over the past decades, Ly- surveys have improved in accuracy, size, and depth, resulting in large samples of mediumresolution spectra from the extended Baryon Oscillation Spectroscopic Survey (eBOSS; Dawson et al. 2016) and the currently observing Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration et al. 2016).Whilst the eBOSS DR16 sample consists of 205, 012 quasar spectra with absorber redshifts in the range 1.96 ≤  ≤ 3.93, DESI observes the sky with increased spectral resolution of 2000 < ∼ R < ∼ 5500 (Abareshi et al. 2022) and covers a much larger fraction of the sky (14,000 deg 2 ).Over its survey duration it will provide ∼ 840, 000 quasar spectra at  > 2.1, a more than fourfold increase compared to BOSS/eBOSS (Adame et al. 2024).To complement the picture, smaller samples of highresolution measurements from, e.g., the High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994;O'Meara et al. 2021) and the UV-Visual Echelle Spectrograph (UVES; Dekker et al. 2000;Murphy et al. 2019) allow for analyses deeper into the small-scale regime to almost 10 larger  max .A number of future surveys will also capture spectra in both high and medium resolution modes, such as the WEAVE-QSO survey (Pieri et al. 2016), the Prime Focus Spectrograph (PFS; Greene et al. 2022) and 4MOST (de Jong et al. 2019).A fast and robust estimator for the analysis of Ly- forest at all scales is essential for accessing the full potential of the combination of these data sets.
One of the key advantages of the Ly- forest is that, since the density fields are only mildly non-linear at the respective redshifts, a much wider range of scales can be used to robustly probe cosmology than with most galaxy surveys.This makes the Ly- fluctuations a particularly powerful probe of early-Universe physics when combined with tracers that are sensitive to very large scales, e.g, cosmic microwave background anisotropies.The Ly- flux density contrast,   , is used to derive constraints for cosmology and astrophysics, defined through where   () is the observed transmitted flux as a function of , the optical depth,  ( Ly- ) is the mean transmitted flux at the HI absorber redshift and () is the unabsorbed continuum of the background quasar.The three-dimensional power spectrum,  3D , of the fluctuations in the Ly- forest, given in Eq. ( 1), is principally sensitive to the amplitude of dark matter clustering  8 , the shape of the matter power spectrum Γ = Ω m ℎ, the spectral tilt   and the sum of neutrino masses   (see, e.g., Philcox & Ivanov 2022).
In this work, we present a measurement of the anisotropic 3D Ly- forest power spectrum using eBOSS DR16 Ly- forest spectra.We compare our measurement to the best-fit results from the 2PCF analysis in the range (0.019 ≤  ≤ 0.35) ℎ Mpc −1 (du Mas des Bourboux et al. 2020, hereafter dMdB20).It includes a Kaiser formula (Kaiser 1987) with a non-linear correction term obtained from hydrodynamical simulations (Arinyo-i-Prats et al. 2015) and modeling of contaminants (e.g., metals, distortion from continuum fitting, damped Ly- absorbers) specific to Ly- forest analyses.Predictions for the Ly- forest clustering have greatly improved over the last decade, using either hydrodynamical simulations (see, e.g., Arinyo-i-Prats et al. 2015;Bolton et al. 2017;Puchwein et al. 2023;Doughty et al. 2023) or perturbation theory 2 (see, e.g., Givans 2 On analyses of the connection between the physics of the Ly- forest & Hirata 2020; Garny et al. 2021;Chen et al. 2021;Givans et al. 2022;Ivanov 2023).
In the context of simulations, different P3D methods have been proposed and tested, see McDonald (2003); Arinyo-i-Prats et al. (2015); Font-Ribera et al. (2018); Horowitz et al. (2024).The latter two methods provide close-to-optimal ways of constructing a covariance matrix.Whilst we measure the anisotropic clustering in multipole space and as a function of  = √︃  2 ∥ + k 2 ⊥ , Font-Ribera et al. ( 2018) measure the P3D in the (more natural basis) { ∥ , k ⊥ }basis (which is the Fourier space equivalent to the { ∥ , r ⊥ }-basis of the configuration space analysis in, e.g., dMdB20).Whilst the information content is (in principle) the same as in the present analysis, the { ∥ , k ⊥ }-basis allows to isolate (and marginalize out) modes that are, e.g., sensitive to distortions in the continuum fitting ( ∥ = 0) (see Font-Ribera et al. 2018, for a discussion).In the present work we compute a simulation-based covariance matrix, whereas Font-Ribera et al. (2018) compute an approximate version of the global covariance matrix within the optimal quadratic estimator framework.A recent proof-of-principle presented in Karim et al. (2023) presents a small-scale measurement on simulations and eBOSS DR16 data which measures the cross-spectrum  × (, ,  ∥ ), originally developed in Hui et al. (1999);Font-Ribera et al. (2018).
The remainder of this paper is organized as follows: We present the 3D power spectrum estimator in Sec. 2, before discussing the theoretical modeling of the Ly- flux power spectrum in Sec. 3 and the forward modeling of the window function in Sec.3.1.In Sec. 4 we describe the observational eBOSS DR16 Ly- forest data set as well as the employed synthetic Ly- spectra used to model covariances.In Sec. 5 we extensively test our pipeline on mocks before presenting the main result of our analysis in Sec.6: the eBOSS DR16 three-dimensional power spectrum of the Ly- and perturbative approaches, i.e., response function approaches, see Seljak (2012); Cieplak & Slosar (2016); Iršič & McQuinn (2018).
forest.Sec.7 presents our conclusions.We make our HIPSTER-lya implementation publicly available.3

THE 3D POWER SPECTRUM ESTIMATOR
In cosmological data analysis, one typically employs data compression, both to reduce the size of the data vector to a computationally manageable level and average over stochastic fluctuations.The particular compression statistic is specific to the problem at hand; for Ly-, one typically uses the 2PCF or the one-dimensional power spectrum along the line-of-sight (e.g., Slosar et al. 2013;Busca et al. 2013;Palanque-Delabrouille et al. 2013;Ravoux et al. 2023;Karaçaylı et al. 2024).Here, we compress the large number of Ly- forest spectra instead into the 3D power spectrum,4 paving the way for future full-shape cosmological analyses in Fourier space of eBOSS data (and beyond). 5 Usually, Ly- summary statistics are computed using quadratic maximum likelihood (e.g., McDonald et al. 2006;Karaçaylı et al. 2020) or FFT-based algorithms (e.g., Palanque-Delabrouille et al. 2013;Chabanier et al. 2019); in this work, we instead make use of configuration space pair-count estimators, following the HIPSTER algorithm of (Philcox & Eisenstein 2020;Philcox 2021).These were originally introduced in the context of the small-scale galaxy power spectrum and bispectrum, and translate well to our problem, given the sparse (dense) sampling of pixels transverse to (along) the lineof-sight, which typically leads to aliasing for FFT-based methods.As we discuss below, HIPSTER estimates spectra by summing pairs of points in configuration space, weighted by exp (k • r   ), for wave vector k and pair separation by r   ≡ r  − r  .
Given some data field  (x) (encoding Ly- fluctuations or galaxy overdensities, with some weighting), the (pseudo-)power spectrum is defined as a Fourier-transform of the two-point function : where F is the Fourier transform operator, and we have assumed  to be suitably normalized.Usually, one considers angularlyaveraged power spectra, defined by 6 where we have restricted to a -bin , with volume   = ∫ k Θ  (|k|), for top-hat function Θ  ().Comparing to Eq. ( 2), we see that the k-dependence is captured through the kernel where  0 is the 0-th order spherical Bessel function.In the last step, we assumed the thin-shell limit and aligned the | k|-axis with |r  − r  |.Noting that  (x) = (x)(x)(x) for background density spectra measured on the data.Theory spectra are denoted by  ( ) which, when convolved with the window matrix, are labeled P ( ).The multipoles of the power spectra are denoted by the subscript ℓ, i.e.,  ℓ ( ).
5 For a compressed full-shape analysis in configuration space, see, e.g., Cuceu et al. (2023). 6We use the shorthand notation: , weights , and overdensity  (which is the quantity of interest), the power spectrum can be written explicitly as adding a normalization term appropriate to an ideal uniform survey of volume .In this approach, the power spectrum can be computed by explicitly evaluating Eq. ( 5), i.e. counting each pair of Ly- pixels, weighted by above kernel.
In practice, we also wish to estimate the anisotropy of the Ly- power spectrum.This is possible by a simple generalization for multipole ℓ: where ∡[k, 1 2 (r  + r  )] denotes the cosine of the angle between k and the local line-of-sight (r  + r  )/2.Eq. ( 6) is the basis of the three dimensional power spectrum estimator used herein 7 .
In practice, all pairs of Ly- pixels contribute to the measured power spectrum, which results in the naïve estimator of Eq. ( 6) being extremely expensive to compute.In practice however, wellseparated pairs have negligible contributions to the power spectrum, thus we can truncate the pair counts at separations greater than  0 ≥ 200 ℎ −1 Mpc , with minimal loss of information (Philcox & Eisenstein 2020).This is done via a smooth polynomial window function which modifies the anisotropic power spectrum estimator in Eq. ( 6) by a factor  ((r  − r  )/ 0 ).Though this necessarily induces slight distortions to the measured spectra, its effect can be straightforwardly included in the theoretical model for   ℓ , thus avoiding any potential bias.
Finally, we must consider the window function of the data, denoted by Φ(r) ≡  (r)/ (r), where  (r) = ⟨(x)(x + r)⟩, is the true 2PCF of the data.Due to the sampling of the Ly- skewers, this has a non-trivial form, and leads to our estimators returning a window-convolved power spectrum measurement.Explicitly, this effect can be included by forward modeling the window function in the theoretical model, using where  (r) is the correlation function model, and P(k) is the output window-(and pair-separation-)convolved power.This procedure will be discussed in more detail in Sec.3.1.
In discrete form, the full estimator for the Ly- power spectrum used in this work is given by a sum over each pair of pixels, , : where we drop the (arbitrary) normalization factor and additionally remove the 'self-skewer' counts; an advantage of doing the measurement in configuration space to avoid (correlated) uncertainties stemming from continuum estimation.The weights  are obtained from the continuum fitting (dMdB20), and  , are the estimated Ly- fluctuations.We use the same set of pixels to compute the window Φ (needed to forward-model geometry effects): restricting to radial bin  and Legendre multipole ℓ.The measurement of the multipoles of the power spectrum and corresponding window functions for a data set with ∼ 108 pixels (similar to the data set used in this work) takes ∼ 5 hours on the Perlmutter computer at NERSC (using one AMD Milan CPU).

THEORY POWER SPECTRUM MODELING
One of the key advantages of the Ly- forest is that it probes the Universe at redshifts 2 ≤  ≤ 5, i.e., with many more linear modes than seen in galaxy surveys.For our theory modeling, we use the same procedure to the one outlined in dMdB20 and briefly summarize it below. 8The real-space quasi-linear power spectrum9 ,  QL (, , ), is modulated by a simulation-based non-linear fitting function,  NL , and connected to the redshift-space flux power through the wellknown Kaiser formula (Kaiser 1987) where   denotes the redshift-dependent linear bias parameter,   is the redshift-space distortion (RSD) parameter (which is assumed to be redshift-independent) and  is the angle of k = { ∥ , k ⊥ } to the line-of-sight,  ≡  ∥ /.The redshift-evolution of the Ly- forest bias parameter enters the weights, , which are taken from Eq. 7 in dMdB20, and enters the power spectrum measurement through Eqs. ( 9) and (10).Following McDonald et al. (2006); dMdB20; Gordon et al. (2023b), the product of Ly- forest bias and growth factor is assumed to have a redshift dependence described by (1 + )  Ly- −1 where  Ly- = 2.9 is the redshift evolution parameter of the Ly- forest bias.Over sufficiently large scales the Ly- forest fluctuations are, to linear order, given by   =  , + ,  where  are mass fluctuations and  = −(  /  )/ is the gradient of the peculiar velocity   over the comoving line-of-sight coordinate   with the scale factor  and the Hubble constant .As such, the bias factors are the partial derivatives with respect to  and  (see, e.g., Arinyo-i-Prats et al. 2015).The redshift space distortion parameter for Ly- forest surveys is for our fiducial cosmology we obtain  = 0.9704.We can estimate   and   from simulations, by either fitting the (i) auto-power spectrum using the Kaiser formula (e.g., Arinyo-i-Prats et al. 2015); (ii) cross-power of Ly- and matter (e.g., Givans et al. 2022); or (iii) from separate Universe simulations (e.g., Cieplak & Slosar 2016).
For the Ly- forest this results in an RSD parameter greater than unity, yielding a quadrupole that is larger than the monopole (in contrast to the case for galaxies).To account for the broadening in  ∥ stemming from high column density (HCD) systems (see, e.g., Rogers et al. 2018, and references therein), the Ly- bias parameters (, ) are remapped to where  HCD ( ∥ ) = exp(− HCD  ∥ ) is fitted to hydrodynamical simulations and  HCD is the typical scale for unmasked HCDs, set to 10 ℎ −1 Mpc in dMdB20.
The non-linear (NL) correction to the power spectrum is obtained from fits to hydrodynamical simulations (Arinyo-i-Prats et al. 2015) with the usual dimensionless power spectrum Δ 2 () ≡  3  lin ()/(2 2 ).The parameter  1 is a dimensionless parameter and controls the importance of the non-linear enhancement. 10hile  1 is only (very) weakly dependent on redshift, it is, together with   and     , inversely proportional to the amplitude of the linear power spectrum.The remaining parameters model the Jeans smoothing and line-of-sight broadening stemming from non-linear peculiar velocities and thermal broadening (see Mc-Donald 2003;Arinyo-i-Prats et al. 2015).We set the parameters to  1 = 0.8558,   = 1.11454 ℎ −1 Mpc ,   = 0.5378,   = 1.607 and   = 19.47 ℎ −1 Mpc , following dMdB20.
The quasi-linear power spectrum is given by which decomposes the power spectrum into a smooth (no BAO feature) and a peak (isolating the BAO feature) component.This description includes a correction for the non-linear broadening of the BAO peak denoted by Σ ∥ /Σ ⊥ = 1 +  where  is the (linear) growth rate (see, e.g., Eisenstein et al. 2007).Analogously to dMdB20, we set the smoothing parameters to Σ ⊥ = 3.26 ℎ −1 Mpc and Σ ⊥ = 6.42 ℎ −1 Mpc for a growth rate of  ≈ 0.97.Contaminants in the Ly- forest, such as absorption by metals and correlations due to the sky-subtraction procedure of the eBOSS data pipeline, are modeled as additive terms in the Ly- auto-correlation function (see Sec. 4 in dMdB20, for a detailed discussion).The computation of comoving separations is based on the (erroneous) assumption that all absorption stems solely from the H i Ly- transition.Transitions from other elements (e.g.silicon and carbon often denoted as metals) contribute to the measured Ly- forest fluctuations as well.Thus, the measured absorption at a given wavelength is a mixture of absorption at different redshifts (from different transitions).For each pair of possible transitions the offset between the true and assumed redshift is computed.The resulting re-mapping matrix between the true and assumed comoving separation is called the metal matrix, see, e.g., Blomqvist et al. (2018), dMdB20, relating the measured to the 'true' summary statistic in configuration space, e.g. the correlation function.
To estimate the window-convolved power spectra, we will require the correlation function multipoles: these are defined as an angular integral over the full correlation function , As described below, these will then be used to compute the mode mixing introduced by the window matrix in Eq. ( 23), allowing accurate comparison of theory and data.

Forward modeling the window function
In the present analysis, we forward model the effect of the survey geometry, captured by the configuration-space window matrix Φ(r) ≡  (r)/ (r), on the theoretically expected power spectrum instead of removing the window function from the data.This factor is independent of the Ly- fluctuations,   , and can be explicitly computed by counting pairs of Ly- pixels, as in Eq. ( 10).This is analogous to the approach used in galaxy surveys e.g., Castorina & White (2018); Beutler et al. (2017); Beutler & McDonald (2021), whence one writes for (weighted) background density (x) = (x)(x), and an ideal bin volume  model (r) = 4  3 ( 3 max − 3 min ), given some set of (thin) bins in    .In the galaxy case, this is estimated using catalogs of random particles; for us, the procedure simply involves counting the unweighted Ly- pixels. 11 Expressed in terms of (even) Legendre multipoles, the paircount estimator of Eq. ( 9) measures the power spectrum convolved with the functions where  (;  0 ) is the pair-truncation function given in Eq. ( 7).
The convolved power spectrum multipoles can be obtained by first considering the distortion induced to the 2PCF (see, e.g., Beutler et al. (2017)): 11 In full, the Ly- case can be imagined as a set of (correlated) pencil-beam surveys, whose sampling is known precisely.
From Eq. ( 8), the power spectra can then be expressed as a Hankel transform: where   ℓ () is the usual spherical Bessel function  ℓ () integrated over a -bin to reduce oscillations: as in Eq. ( 3).12For a thin -bin centered at   ,   ℓ () ≈  ℓ (  ).In this work, we truncate our expansion of the window-convolved correlation function, ξℓ , at ℓ max = 4 for the 2PCF.Note that we include all corresponding window matrix multipoles for each multipole in ξℓ . 13In Appendix A, we discuss the contribution of each correlation function multipole to the final result.

LYMAN-A FOREST DATA FROM EBOSS DR16
In the present analysis, we compute the 3D power spectrum from Ly- forest spectra from the 16th data release (DR16Q; Lyke et al. 2020) of the completed extended Baryonic Oscillation Spectroscopic Survey (eBOSS; Dawson et al. 2016) of the fourth generation of the Sloan Digital Sky Survey (SDSS-IV; York et al. 2000).The primary scientific goal of eBOSS was to constrain cosmological parameters using BAO and redshift space distortions (see, e.g., Alam et al. 2021;du Mas des Bourboux et al. 2020).Briefly summarized, the DR16 data consists of the complete 5-year BOSS and 5-year eBOSS survey.The Ly- spectra were observed with a double spectrograph mounted on the 2.5m Apache Point telescope to map 10,000 deg 2 of the sky.The observations are conducted in the observed wavelength range of 3600 <  obs < 10, 000 Å with a spectral resolution of  ∼ 1500 − 2500 (Dawson et al. 2016).The quasar selection and algorithms of the pipeline are explained in detail in Myers et al. (2015) as well as technical details related to the survey itself in Dawson et al. (2016).

Data selection
The full sample consists of 205, 012 quasar spectra with absorber redshifts in the range 1.96 ≤  ≤ 3.93.The quasar sample is split into two disjoint regions on the sky; the northern (NGC) and southern (SGC) galactic caps, with 147, 392 and 57, 620 sight lines respectively.The number of spectral pixels are 34.3 • 10 6 for the entire data set with an average signal-to-noise ratio of 2.56 per pixel calculated in the Ly- forest.The forest used for cosmological analysis is defined to be 1040 ≤  rf ≤ 1200 Å.For ease of comparison to dMdB20, we use the same set of flux decrements,   , introduced in the following section.We note that BOSS/eBOSS observed spectra in spectral pixels of with Δ log 10 () ∼ 10 −4 which have been rebinned for the purpose of the final analysis onto a grid of Δ log 10 () ∼ 3 × 10 −4 .Damped Ly- absorbers (DLAs), defined as regions where the transmission is reduced by more than 20% in the flux decrements, are masked out.A Voigt profile is fitted to each DLA to correct for the absorption in the wings (Noterdaeme et al. 2012).Quasars with broad absorption lines (BALs) have been removed for the present analyses (see Lyke et al. 2020, for more details).Both are the main contaminants affecting the data selection.

Continuum fitting
To measure clustering statistics such as the P3D, we use the flux decrement given in Eq. ( 1).Extracting the unabsorbed flux, i.e., the quasar continuum, is a daunting task in general.In the present paper, we use the publicly available continuum-fitted spectra from eBOSS DR1614 .In the following, we briefly outline the continuum fitting procedure from dMdB20: to each Ly- forest a slope and an amplitude is fitted to the stacked spectra in the quasar sample using the picca package. 15The product in the denominator of Eq. ( 1) assumes a common continuum for all quasars in restframe, ( rf ), corrected by a first order polynomial for each quasar  in Λ ≡ log  which accounts for the diversity of the quasars which results in a biased mean and spectral slope of each forest flux decrement δ() =  ()/ (  +   Λ())() −1.The fluctuations are then given by with   =  ′ ( ′ ) and Λ() = log −log .We use this quantity to measure the power spectrum.Estimating the continuum directly from the forest distorts the field, requiring correction at the level of the flux decrement (see, e.g., dMdB20).The continuum fitting and projection to center the mean flux decrement at zero suppresses modes along the line-of-sight in the forest.A good approximation is to treat these distortions as linear and model them in the theoretically expected correlation function through a distortion matrix (DM; dMdB20) where  and  ′ denote two bins of the correlation function with for forest pixel pairs (, ) and ( ′ ,  ′ ) in sight lines  and  ′ .Whilst this approach marginalizes out large-scale modes, it lends itself well to measure the small-scale continuum for the one-dimensional power spectrum.In the current implementation we apply the distortion matrix to the theory correlation function in { ∥ , r ⊥ }−space. 16We denote the theory correlation function to be 'distorted' after applying the distortion matrix.)The distorted multipoles are then used to compute the window convolved theory power spectrum in Eq. ( 23).See Appendix B for the effect of the distortion matrix on theory power spectra and correlation function multipoles.

CONSISTENCY TESTS ON SIMULATIONS
Before presenting the main results of this paper, we first test our pipeline by applying the power spectrum estimator, presented in Sec. 2, to idealized Gaussian realizations, with results displayed in Sec.5.1.In Sec.5.2 we apply our pipeline to realistic eBOSS DR16 mocks including contaminants (such as high column density absorbers and metals).In particular, we discuss the effect of continuum fitting and the extraction of the unabsorbed flux, on the resulting power spectrum measurements.

Gaussian random field simulations
The first series of tests of the power spectrum estimator presented in Sec. 2 is on anisotropic Gaussian random fields (GRFs), generated with a known linear matter power spectrum  lin ().The aim is to investigate biases as well as the range of validity and robustness of the estimator to non-trivial survey geometries with a known input matter power spectrum. 17The GRFs are realized in a periodic box of size  = 1380 ℎ −1 Mpc with  = 512 cells at redshift  = 2.4 with a ΛCDM linear power spectrum as theory input (Planck Collaboration 2020).The fundamental mode is   = 2/ = 0.0045 ℎ Mpc −1 with the Nyquist frequency being  Ny = 1.17 ℎ Mpc −1 .We sample  s = 900 sight lines from the box, corresponding to a (very sparse) quasar density of ∼ 2−3 qso/deg 2 , and 460, 800 Ly- pixels. 18For the Ly- forest fluctuations, , we use the pixel value, i.e., the modes from the GRF, and model the selection function using weights of unity.For the power spectrum measurement, we do not include pixels in the same skewer, as in the main Ly- analysis.Additionally, we apply a non-trivial selection function, i.e., with edges and sparse sampling transverse to the line-of-sight to measure the effect of the window function on the measured power spectrum.To this end, we remove the periodicity of the box and add a 10% root-mean-squared error, relative to the mean, Gaussian, white (uncorrelated), noise to each of the pixels in the box, corresponding to somewhat realistic Ly- noise levels for continuum fitting and pipeline noise.
We measure the monopole, quadrupole and hexadecapole power spectra from the GRFs and compute the 2PCF contributions up to ℓ max = 4 with all corresponding window matrix multipoles.We use 30 equidistant -bins in the range (0.01 ≤  ≤ 0.50) ℎ Mpc −1 and   = 1, 000 linearly spaced bins in the range (0 ≤  ≤  0 = 200) ℎ −1 Mpc for the window matrix, respectively. 19We compare the measured power spectrum multipoles to the theory input power spectrum in Fig. 1 and recover the input spectrum well within the 1 error bars up to  < ∼ 0.5 ℎ Mpc −1 with the measured anisotropic power spectra.Note that the quadrupole is larger than the monopole due to the large RSD parameter   = 1.5, typical for Ly- surveys.Note that the pair separation window  (;  0 ) suppresses power on scales larger than the ones defined the 2PCF (and the corresponding window matrix multipole).
MNRAS 000, 1-15 ( 2024 by  = 2/ 0 which is at  = 0.031 ℎ Mpc −1 for our chosen value of  0 .Thus, we discard -values below the cutoff in our analysis pipeline.In Appendix A, we discuss the effect of the number of skewers on the multipoles of the window matrix itself as well as their convolution with the window convolved theory power spectra. 20he error bars shown in Fig. 1 capture the variation between differing GRF realizations and are obtained from the diagonal of the covariance matrix computed from  = 100 realizations of the GRFs.For Legendre multipoles ℓ 1 , ℓ 2 of the anisotropic power in interchangeably.
-bins , , this is given by with the mean of the power spectra given by where  indexes the  simulations.The correlation matrix is defined as by construction, this is unity along the diagonal.We show all combinations of the correlation matrix moments in Fig. 2: the first column displays the correlation between -bins in the monopole with itself, the quadruple and the hexadecapole (from bottom to top, respectively).Analogously, the second column displays the correlations of the quadrupole with the same power spectrum multipoles and the third column shows the correlations of the hexadecapole with said moments.The anisotropies sourced by the imposed survey geometry and intrinsic correlations from the generated anisotropies are (weakly) visible as correlations along the diagonal between the monopole and quadrupole (and the quadrupole with the hexadecapole, respectively.).We find negligible correlation between the monopole and hexadecapole (the latter being mostly dominated by noise).For the ℓ 1 = ℓ 2 correlations, we observe a correlation length of up to two -bins, resulting in an approximately bi-diagonal correlation matrix.Note that for the quadrupole and the hexadecapole, the lowest three -bins are highly correlated.This can be removed by increasing  0 in the pair separation function, though we caution that computation time scales as O ( 2 pix  3 0 ).Whilst the above tests on GRFs give us confidence that the P3D estimator recovers the input power spectrum to high precision in the presence of survey geometries typical for Ly- data analysis, we stress that these simplistic simulations ignore crucial systematic and instrumental effects as well as details of the data reduction pipeline present in real data, e.g., the co-addition of individual spectra, per spectral pixel noise estimates, sky residuals, continuum fitting, metal contamination, damped Ly- absorbers and broad absorption lines.Whilst resolution mostly affects high- modes, the key systematic arising along the line-of-sight, i.e. at large scales, is the uncertainty in the continuum measurement. 21In the following section, we address (some of) these issues using more realistic Ly- forest synthetic data.

Synthetic Ly-𝛼 spectra: LyaCoLoRe simulations
In this section, we present tests of our estimator on realistic Ly- forest simulations (for a more complete discussion of the LyaCoLoRe mocks, see Farr et al. (2020) and dMdB20).In the following we will briefly summarize the key steps to generate these simulations: Each realization is based on a GRF of length 10 ℎ −1 Gpc and 4096 3 particles yielding a resolution of ∼ 2.4 ℎ −1 Mpc generated using CoLoRe's log-normal density model (Ramírez-Pérez et al. 2022). 22he input cosmology is based on the best fit ΛCDM-parameters obtained by the Planck satellite (see column 1 of table 3 in Planck Collaboration et al. 2016).The Gaussian field is subsequently transformed to a log-normal density field which is Poisson sampled assuming a quasar density of 59 quasars per deg 2 and a functional form for the quasar bias redshift evolution (see, e.g., dMdB20).The line-of-sight skewers are computed by interpolating from the GRF onto the 'observed' pixel positions and the radial velocity field to the center of the box.The obtained skewers are then post-processed using the LyaCoLoRe package (Farr et al. 2020;Herrera-Alcantar et al. 2023). 23The sparse sampling of the Ly- forest transverse to the line-of-sight (also denoted by -sampling) results in power on small scales contributing to the error on large scales (see, e.g., McDonald & Eisenstein 2007) which is quantified by the one-dimensional power spectrum (P1D).Thus, additional power is added to each skewer by sampling from a Gaussian with variance set by the P1D.From the log-normal transformation of the final skewer the optical depth field, , is computed using the fluctuating Gunn-Peterson approximation (FGPA; Croft et al. 1998).Redshift space distortions are obtained by convolving the  field with the peculiar velocity field.The observed flux is then related to  via  = exp (−).
The LyaCoLoRe mocks mimic noise and instrumental systematics present in the eBOSS data. 24In particular, the effect of continuum fitting, instrumental noise, high column density (HCD) absorbers 25 , Ly- absorption and metals.Additionally, random redshift errors for the spectra are included by drawing from a Gaussian with mean zero and   = 400 km/s.The eBOSS analysis pipeline and spectral resolution are simulated using the specsim package (Kirkby et al. 2021). 26The co-addition of the resulting skewers with instrumental noise and a realistic quasar continuum is done using the desisim package. 27 In the following, we present tests of our power spectrum estimator on  realizations of mocks with increasing levels of systematics.We compare four different set of mocks (and use a similar nomenclature as the one adopted in table 4 of dMdB20): (i) Ly- raw mocks (eboss-raw): We estimate the P3D directly from the simulated Ly- forest fluctuations.This analysis is similar to the case of the GRFs but using the eBOSS Ly- survey geometry with the corresponding effective volume.We use  = 10 realizations.
(ii) distorted Ly- raw mocks (eboss-raw-dist): We include a variation of the 'raw' mocks of which we remove the mean of each skewer and denote them by 'raw-dist'.We use  = 10 realizations.
(iii) +continuum+noise (eboss-0.0):We add the effect of measuring quasar continua of the spectra and include instrumental noise.The performed continuum fitting, introduced in Sec.4.2, requires forward modeling the suppression of modes along the line-of-sight through a so-called distortion matrix given in Eq. ( 27).We use  = 10 realizations.
(iv) + metals + HCDs +   (eboss-0.2):We estimate P3D from realistic Ly- mock spectra including the effects stemming from analyses and have also been tested in the context of full-shape analyses of the Ly- forest with the 2PCF (see, e.g., Cuceu et al. 2023), these have not been explicitly developed for the purpose of P3D analyses. 25Following dMdB20, we denote systems with neutral hydrogen exceeding 10 17.2 atoms cm −2 as HCDs. 26Publicly available at https://github.com/desihub/specsim. 27 Publicly available at https://github.com/desihub/desisim. continuum fitting and instrumental noise (from step ii) and adding metals, high column density absorbers and random redshift errors.
Analogous to the analysis on real data, HCDs are masked based on the cut log  H i > 20.3 cm −2 .We use  = 100 realizations.
In Fig. 3 we compare the power spectra of the LyaCoLoRe mocks for the four different mock configurations, alongside their errors.As in the main analysis, we measure the power spectra up to ℓ max = 4 and compute the window function multipoles up to ℓ max = 8.Therefore, we use 24 equidistant -bins in the range (0.02 ≤  ≤ 0.35) ℎ Mpc −1 and   = 1, 000 linearly spaced bins in the range (0 ≤  ≤  0 = 200) ℎ −1 Mpc for the window matrix, respectively.We treat the mean power spectra measured from the 'raw' Ly- forest fluctuations mocks (solid colored lines in the plot), as ground truth to assess the effect of contaminants on the power spectra. 28irst, we consider the effects of removing the mean and the slope of each line of sight to mimic the effect of the distortion matrix (dashed-dot lines in Fig. 3).This isolates the effect of the distortion matrix modeling and can be seen to reduce the monopole power on all scales with most pronounced impact at large scales where it reduces the power up to a factor of two.The quadrupole is also affected at all scales, leading to an increase in power by approximately 20%.The hexadecapole is also affected at all scales, and exhibits a sign change at large scales.Note that the effect becomes increasingly small above  ≥ 0.3 ℎ Mpc −1 for all spectra.
Second, we assess the impact of the continuum and instrumental noise (dashed lines in Fig. 3), denoted by 'eboss-0.0',discussed in Sec.B. This affects all multipoles at large scales.For the monopole it removes, as expected, power at  < ∼ 0.10 ℎ Mpc −1 , since the largescale modes are projected out of the continuum.For the quadrupole (red dashed line) it enhances power at all scales with a further en- survey geometry.This shows the power spectra computed from an isotropic fiducial correlation function (without stochasticity), estimated from skewers taken from a LyaCoLoRe mock using pairs up to  0 = 200 ℎ −1 Mpc .The dots represent the measurement of the P3D estimator (black, red and green for ℓ = 0, 2, 4, respectively).The dashed lines are the fiducial correlation function,  ( ) multiplied by the window function and the pair count separation window, Φ ℓ  (;  0 ) (same color coding).For ease of comparison, we divide the monopole by a factor of 10.
hancement of 10 − 15% up to  < ∼ 0.15 ℎ Mpc −1 .The hexadecapole is strongly affected and switches sign below  < ∼ 0.20 ℎ Mpc −1 .Note that we do not expect the 'raw-dist' mocks to agree exactly with the '0.0' mocks since removing the mean from each line of sight is not exactly the same as introducing a continuum in each spectrum (the latter also affecting the weights).These differences will propagate into the distortions and yield the observed qualitative agreement with a small offset.
Thirdly (dotted lines in Fig. 3), we include metals, HCDs and random redshift errors.For the monopole (dashed to dotted) this adds ∼ 10% power at all scales, creating a small off-set.For the quadrupole this mostly affects scales up to  < ∼ 0.15 ℎ Mpc −1 .The hexadecapole, however, is barely affected.We see that the measurements are noisier and that the power is further decreased in two regions around  ∼ 0.05 ℎ Mpc −1 and  ∼ 0.15 ℎ Mpc −1 .

Anisotropies from the survey geometry
Finally, we illustrate the effects of window function-induced anisotropy by constructing a simple test built upon a known isotropic correlation function,  fid (r) = exp (−|r|/10).To remove the stochasticity of the data set, we replace the         term in our pair count estimator (given in Eq. ( 9)) with the explicit correlation function  fid , evaluated at the pair separation r   = r  − r  .We apply the resulting estimator to the LyaCoLoRe mocks as before, which tests the sensitivity of the present approach to recovering anisotropies purely sourced by the survey geometry, i.e., our measurement of the window matrix.Although the input correlation function is isotropic (with  ℓ>0 = 0), mode-mixing is expected to occur due to the window; thus Pℓ () will be non-zero in practice.To compare data and theory, we use Eqs.( 20)-( 22), as before, setting  ℓ>0 = 0.In this limit,  ℓ () =  0 ()Γ 2 ℓ(). 29n Fig. 4 we show the results using all the skewers from our 'eboss-0.2'mocks which mimic the realistic eBOSS survey geometry. 30The -fields of these mocks have HCDs that are masked out, i.e. this introduces an additional level of realism in reconstructing the window function.We recover the monopole (black), quadrupole (red) and hexadecapole (green) to better than 0.1% precision over the entire -range of (0.0 <  ≤ 1.0) ℎ Mpc −1 , indicating that our modeling of the window function is highly accurate.This gives us confidence in the present approach to measure the three-dimensional power spectrum at all scales from Ly- data.

RESULTS
In this section we use the pair-count estimator described in Sec. 2 to measure the three-dimensional power spectrum from 205,012 eBOSS DR16 quasar spectra.In Sec.6.1 and Fig. 5 we present the key result of the paper: the P3D measurement on eBOSS data in a single redshift bin 2 ≤  ≤ 4 for both data sets (NGC and SGC).
We compare the measurements to window convolved theory power spectra described in Sec. 3. Throughout this section we adopt a Planck 2016 best-fit ΛCDM cosmology: Planck Collaboration et al. 2016).For the Ly- forest, we use an effective mean redshift of  = 2.334 and for theory power spectra, given in Eq. ( 11), we use the eBOSS DR16 best-fit values from the Ly- auto-correlation presented in table 6 in dMdB20.The BAO parameters are defined as usual where   and   are the Hubble and comoving angular diameter distances, respectively, evaluated at some effective redshift,  eff .The denominator is evaluated at a fiducial ΛCDM cosmology.The bias and RSD best-fit parameters are (dMdB20) the metal biases with the associated restframe wavelength in parentheses are and the bias parameters for the HCD systems (Rogers et al. 2018) In Sec.6.2 we discuss our error bars and covariance matrix for the P3D measurement obtained from  = 100 realistic Ly- simulations.Increasing the pair separation in the pair truncation window function  (;  0 ) to  0 = 400 ℎ −1 Mpc did not improve the agreement between theory and measured power spectrum, thus we do not show the corresponding results in this work.For our main analysis, we measure the power spectra up to ℓ max = 4 in 24 equidistant -bins in the range (0.02 ≤  ≤ 0.35) ℎ Mpc −1 and the window function multipoles up to ℓ max = 8 with   = 1, 000 linearly spaced bins in the range (0 ≤  ≤  0 = 200) ℎ Mpc −1 , respectively.
0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.05 0.10 0.15 0.20 0.25 0.30 0.35 Figure 5. Power spectrum multipoles measured from eBOSS DR16 data for the NGC (left) and SGC (right), showing the monopole in black, quadrupole in red, hexadecapole in green.The dots give the measurements and the solid lines show the quasi-linear theory predictions using eBOSS DR16 best-fit values (dMdB20).The error bars are the square root of the diagonal of the covariance matrix obtained from  = 100 realistic LyaCoLoRe mocks.The window function contributions are computed up to ℓ max = 8 and convolved with the distorted theory power spectrum.All measurements use the pair-count estimator of Eq. ( 9), truncated at  0 = 200 ℎ −1 Mpc via a continuous window function  (;  0 ).

Power spectrum multipoles from eBOSS DR16
In Fig. 5 we show the three dimensional power spectrum multipoles measured on eBOSS DR16 data for the NGC (SGC) in the left (right) panel and compare them to the window convolved quasi nonlinear theory power spectrum multipoles.We show the monopole, quadrupole and hexadecapole from the data and compute the window function multipoles up to ℓ max = 8, as described in Sec. 3. The error bars are given by the square root of the diagonal of the covariance matrix, introduced in Sec.6.2.To account for the employed continuum fitting procedure that projects out large-scale modes, we apply a distortion matrix, further altering the large scale shape. 31he estimator uses a pair separation window truncating the pair count at  0 ≥ 200 ℎ −1 Mpc , as before.
For the monopole, we find excellent agreement with the theory prediction over the entire -range.The quadrupole is also in good agreement with the theory prediction, bar the last few measurement points in the  > ∼ 0.20 ℎ Mpc −1 range.The quadrupole and hexadecapole both have correlated features at large scales: We conjecture that this stems from the mode mixing introduced by the distortion matrix.In particular, the ∼ 1 − 2 discrepancy at large scales in the hexadecapole has to be dealt with cautiously when extracting cosmology from the Ly- forest.In Sec.B and in Fig. B1, we justify this by showing that the forward modeling of the DM does not capture exactly the hexadecapole on large scales: First, the distortion matrix is computed out to separations of > ∼ 400 ℎ −1 Mpc but should, in principle, be computed up to the scales of the survey.Second, it assumes a correlation function measurement in bins of size {4 ℎ −1 Mpc ×  ∥ , 4 ℎ −1 Mpc ×  ⊥ }.In the range  > ∼ 0.2 ℎ Mpc −1 , the quadrupole and hexadecapole measurements and the theory predictions are in ∼ 1 agreement with each other.Note that the error bars from SGC are larger than the ones from the NGC which stems from the smaller statistics, i.e., NGC approximately has three times the number of spectra than the SGC data set.
In Fig. 6 we show the ratio of the quadrupole and hexadecapole to the monopole to obtain an estimate for the RSD parameter,   .Note that although the measurements are quite noisy, we get good agreement for the measured power spectra to the non-linear theory predictions at scales  > ∼ 0.1 ℎ Mpc −1 .For comparison, we include the linear theory prediction (horizontal dashed lines) to guide the eye which we compute by integrating Eq. ( 17) in Fourier space.Assuming non-linear effects are negligible, the Kaiser effect gives the analytic prefactors  ℓ=0 () = 1 + 2/3 +  2 /5,  ℓ=2 () = 4/3+4 2 /7, and  ℓ=4 () = 8 2 /35 for the monopole, quadrupole and hexadecapole, respectively.In this limit, the multipoles of our theory power spectrum are given by  ℓ () =  ℓ () 2  lin (); inserting the value of  = 1.627 gives 1.39 and 0.23 for the ratios of quadrupole to monopole and hexadecapole to monopole, respectively.The errorbars are taken from the diagonals of the covariance matrices, rescaled by the monopole power, and should be taken indicatively rather than at face value.

Covariance matrix from mocks catalogs
To capture the variance between spectra and account for instrumental and systematic noise in the data, we compute a covariance matrix from  = 100 realizations of 'eboss-0.2'LyaCoLoRe simulations, introduced in Sec.5.2.The error bars in Fig. 5 are the square root of the diagonal of the covariance matrix with the corresponding correlation matrix shown in Fig. 7 for which we vary both Legendre multipoles and -bins, as in Sec.5.1.The resulting correlation matrix has a diagonal structure with non-negligible off-diagonal terms and blocks: It is interesting to note that the eBOSS DR16 survey geometry mixes modes between 'neighboring multipoles', e.g., between monopole and quadrupole as well as between quadrupole and hexadecapole.For the ℓ 1 = ℓ 2 correlations, we observe a correlation length of up to two -bins.Note that we do not include any scales with  −1 larger than the pair separation window  0 , thus (in contrast to the GRF correlation matrix) the lowest -bins are not artificially correlated.The resulting covariance matrix is positive definite and can thus be used for cosmological inference.We note that the covariance matrix of the monopole,  ℓ=0 ×  ℓ=0 , shows correlations between -bins of order 10% on larger and up to 30% on smaller scales.This may stem from the application of the distortion matrix, which effectively mixes -modes.Thus, we expect a less correlated correlation matrix when using an alternative continuum fitting method that does not project out large-scale modes at the expense of increased noise in the forest.

SUMMARY AND CONCLUSIONS
The Ly- forest is a treasure trove of cosmological information on the expansion history of the Universe and beyond.Given its high-redshift range (2 ≤  ≤ 4) and sensitivity to Mpc scales and below, it is an ideal tracer to probe a wide range of scales: early-Universe physics from combinations with a large-scale tracer, e.g., the cosmic microwave background anisotropies, from large scales or dark matter models as well as thermal properties of the ionized (cold) IGM from small scales.) have now paved the way for three-dimensional power spectrum analyses.To achieve this goal, we require accurate estimates of the 3D Ly- power spectrum: this is challenging, given the sparse sampling transverse to and dense sampling along the lineof-sight.In this paper, we attack this problem by presenting a pair count estimator to measure the three-dimensional power spectrum from the Ly- forest.Our approach is based on the pair counting estimator HIPSTER (Philcox 2021;Philcox & Eisenstein 2020) and weighs each pair by exp(k • r   ), for wave vector k and particle pair separation r   = r  − r  .This directly measures the power spectrum without the need for grid-based Fourier transforms which are affected by the sparse −function type sampling of the Ly- forest.
We have presented the first large-scale three-dimensional power spectrum measurement of Lyman- forest spectra on 205,012 medium-resolution Ly- forest spectra from the extended Baryon Oscillation Spectroscopic Survey (eBOSS) DR16 in two disjoint regions of the sky: the northern and southern galactic cap.Furthermore, we have compared our P3D measurement to the best-fit quasi-linear theory prediction from the 2PCF presented in dMdB20.We extensively test the estimator on Gaussian random fields, i.e., realizations of a linear input power spectrum, with Gaussian noise and realistic continuum error levels of 10%.We recover the input power spectrum at the ∼ (1) 3 level for (dense) sparse configurations, respectively.Further, our pipeline has been extensively tested on synthetic Ly- spectra (LyaCoLoRe) with increasing levels of realism to probe for the effects of distorted continua, instrumental noise, HCDs, metals and random redshift errors.We demonstrated that we can forward model the effects of anisotropies sourced by the survey geometry and the distortions introduced by the continuum fitting (modes along the line-of-sight,  ∥ = 0) accurately for the multipoles.
We present a covariance matrix derived from 100 realistic LyaCoLoRe simulations for further cosmological analysis of the Ly- forest data set.We obtain a fairly diagonal covariance matrix with correlations between adjacent multipoles stemming from the survey geometry and, for the auto-correlation of the monopole, contaminants and the mode-mixing distortion matrix.We present the anisotropic clustering of eBOSS DR16 Ly- spectra up to the hexadecapole and obtain good agreement to the Kaiser formula with a non-linear correction term obtained from hydrodynamical simulations (Arinyo-i-Prats et al. 2015) up to  ≤ 0.35 ℎ Mpc −1 .The best-fit theory power spectrum obtained from dMdB20 shows deviations at the ∼ 2 level on the largest scales, in particular for the hexadecapole.The quadrupole and hexadecapole mix a large range of scales when convolving the window matrix with the correlation function and Hankel transforming it to Fourier space.This effect is even more pronounced if a distortion matrix (an artifact from the present continuum fitting approach) is applied.We leave for future work to explore alternative continuum fitting methods as well as the accuracy of the presented covariance matrix.
Our main conclusion is that the novel estimator is well suited to measure clustering statistics from the Ly- forest and can deal with non-trivial survey geometries and masked data vectors.This will facilitate robust measurements of P3D for the currently observing Dark Energy Spectroscopic Instrument and future surveys such as the Prime Focus Spectrograph, WEAVE-QSO and 4MOST, opening the door to a wide variety of novel cosmological analyses.(i) sparse with ∼ 2 qso/deg 2 ; and (ii) dense with an eBOSS DR16-like sampling of 30 qso/deg 2 .The dense window function has been normalized by its number density for ease of comparison.
Here, Γ 2 denotes the window matrix obtained from the survey geometry through pair-counting without the   weights.It is interesting to note that the window pair separation function only affects very large scales below  < ∼ 0.08 ℎ Mpc −1 ; this is as expected since the  (;  0 ) function given in Eq. ( 7) was chosen to reduce aliasing (Philcox & Eisenstein 2020).For the monopole spectra (top panel), we find large contributions from the monopole and quadrupole windows at all scales, though negligible effects from the hexadecapole.In the middle panel we show the power spectrum quadrupole with contributions up to ℓ max = 6 where the hexadecapole contributions are strongest and higher order terms have vanishing power.In the bottom panel we show the theory convolved hexadecapole with window matrix contributions up to ℓ max = 8.Here it is important to note that all multipoles of the window matrix contribute visibly to the final window convolved power spectrum.This stems from the integrand of the Hankel transform from the window convolved correlation function to the power spectrum multipoles.Note that the integrands of the quadrupole and hexadecapole mix a wide range of scales; this is discussed in the context of the eBOSS DR16 data in Appendix B. and hexadecapole (green lines), respectively.The solid lines are the best-fit correlation functions obtained from the mean correlation function measurements (dMdB20) which agree well with the mean measured power spectrum for each mock configuration in Fig. 3: solid lines correspond to the ones from 'raw' mocks, dashed-dotted lines from 'raw-dist' mocks including the mean subtraction displaying the effect of the distortion matrix, dashed lines from the 'eboss-0.0'mocks which include the effect of noise and continuum fitting through the distortion matrix and the dotted lines are the correlation function multipoles of the 'eboss-0.2'mocks which include contaminants and effects from continuum fitting.
introduced in Sec.4.2 suppresses modes along the line-of-sight, effectively mixing  fields resulting in distortions in the measured power spectra (or correlation functions).The resulting distortion matrix is computed from the 'eboss-0.0'mocks using Eq. ( 27).When comparing the measured power spectra of the 'raw' to the 'eboss-0.0'mocks we find a suppression of the monopole on all scales which is most pronounced at  < ∼ 0.1 ℎ Mpc −1 by up to ∼ 20%.The quadrupole (red) is enhanced by ∼ 20% on all scales.The hexadecapole (green) switches sign at the largest scales and is strongly damped.We compare the measured spectra to the lineartheory power spectra (gray lines) which are only valid on large scales (since the true power of these mocks is not known).Applying the DM to the multipoles, we find that the forward modeling of the distortion is accurate at the ∼ 5% level.For the quadrupole and the hexadecapole we observe that the discrepancy between the distorted theory and the measurement for the 'eboss-0.0'mocks is reduced at the largest scales by increasing the computed maximum separation of the distortion matrix from 200 ℎ −1 Mpc to 400 ℎ −1 Mpc in 4 ℎ −1 Mpc bins in { ∥ , r ⊥ }.We additionally tested that including higher order contributions of the theory correlation function,  ℓ () up to ℓ max = 6, and including all the window matrix multipoles in Eqs. ( 20)-( 22), did have a negligible effect on the theory prediction of the multipoles.This additional level of complexity can be removed from the analysis pipeline by using a continuum fitting method that does not introduce distortions, e.g., a PCA-based continuum fitting (de Belsunce et al. 2024) at the expense of obtaining noisier continuum estimates.
For comparison we show the effect of the distortion matrix and contaminants on the multipoles of the theory correlation function in Fig. B2 for the four LyaCoLoRe mock configurations.It is interesting to note that the application of the DM introduces a correlation between modes that are much further apart than a skewer length of a few hundred ℎ −1 Mpc 's, i.e. in Fig. B1 quadrupole modes at  ∼ 0.1 ℎ Mpc −1 and  ∼ 0.35 ℎ Mpc −1 are both affected by the DM.We conjecture that for a given mode  the power spectrum quadrupole is sensitive to larger scales in the correlation function  23) from the multipole of the correlation function to the corresponding multipole of the power spectrum (monopole top, quadrupole center, hexadecapole bottom).The pair truncation window  (;  0 = 200 ℎ −1 Mpc ) is applied to the theory correlation functions.The solid (dashed) lines are for  = 0.1 ℎ Mpc −1 ( = 0.3 ℎ Mpc −1 ); red (black) lines are with (with out) the distortion matrix.
quadrupole than for the monopole.Therefore, we show in Fig. B3 the integrand of the Hankel transform of the window-convolved correlation function to the power spectrum multipoles for both wavenumbers  ∼ 0.1 ℎ Mpc −1 as a black solid and  ∼ 0.35 ℎ Mpc −1 as a dashed line, respectively.The application of the distortion matrix (black to red) mixes a wide range of scales for quadrupole and hexadecapole.

Figure 1 .Figure 2 .
Figure 1.Consistency test of the power spectrum estimator, comparing the window-convolved linear power spectrum at  = 2.4 (monopole in black, quadrupole in red, hexadecapole in green, respectively) to the mean power spectrum of  = 100 Gaussian random field simulations (dots), computed using pair counts truncated  0 = 200 ℎ −1 Mpc .The shaded regions give the square root of the diagonal of the covariance matrix and illustrate the variance between the  = 100 realizations; this is also shown by the errorbars, which include the factor of 1/ √  = 1/10.The 2PCF contributions are computed up to ℓ max = 4 with all corresponding window matrix multipoles, as described in Sec. 3.

Figure 3 .
Figure3.Comparison of Ly- power spectra measured from LyaCoLoRe mocks for four different configurations: 'raw' are the raw Ly- forest fluctuations (solid lines); 'raw-dist' are the distorted fluctuations where the mean and the slope along each sight line have been removed (dashed-dotted lines); '0.0' include continua and instrumental noise (dashed lines); '0.2' include metals, HCDs and random redshift errors (dotted lines).The colored lines (black for the monopole, red for the quadrupole and green for the hexadecapole, respectively) are the mean power spectra and the error bars are the variance between the different realizations.

Figure 4 .
Figure 4. Consistency test for anisotropies induced by the (eBOSS-like) survey geometry.This shows the power spectra computed from an isotropic fiducial correlation function (without stochasticity), estimated from skewers taken from a LyaCoLoRe mock using pairs up to  0 = 200 ℎ −1 Mpc .The dots represent the measurement of the P3D estimator (black, red and green for ℓ = 0, 2, 4, respectively).The dashed lines are the fiducial correlation function,  ( ) multiplied by the window function and the pair count separation window, Φ ℓ  (;  0 ) (same color coding).For ease of comparison, we divide the monopole by a factor of 10.

Figure 6 .Figure 7 .
Figure6.Ratio of measured quadrupole (black dots) and hexadecapole (black inverted triangles) power spectra to monopole power for eBOSS DR16 data (left panel: NGC, right panel: SGC ), with linear expectation for   = 1.627 given by the horizontal dashed lines and ratio of the (non-linear) theory power spectra for the ratio of the quadrupole (red line) and hexadecapole (green line) to the theory monopole.The error bars are taken from the power spectrum, i.e., the square root of the diagonal of the covariance matrix obtained from LyaCoLoRe simulations, and are rescaled by the monopole power.

Figure A1 .
Figure A1.Window function of the data, Φ ℓ ( ), for two different configurations discussed in Sec.5.1: (i) sparse with ∼ 2 qso/deg 2 ; and (ii) dense with an eBOSS DR16-like sampling of 30 qso/deg 2 .The dense window function has been normalized by its number density for ease of comparison.

Figure B2 .
Figure B2.Effect of the distortion matrix and contaminants on the theory correlation function multipoles: monopole (black), quadrupole (red) and hexadecapole (green lines), respectively.The solid lines are the best-fit correlation functions obtained from the mean correlation function measurements (dMdB20) which agree well with the mean measured power spectrum for each mock configuration in Fig.3: solid lines correspond to the ones from 'raw' mocks, dashed-dotted lines from 'raw-dist' mocks including the mean subtraction displaying the effect of the distortion matrix, dashed lines from the 'eboss-0.0'mocks which include the effect of noise and continuum fitting through the distortion matrix and the dotted lines are the correlation function multipoles of the 'eboss-0.2'mocks which include contaminants and effects from continuum fitting.

Figure B3 .
Figure B3.Integrand of the Hankel transform given in Eq. (23) from the multipole of the correlation function to the corresponding multipole of the power spectrum (monopole top, quadrupole center, hexadecapole bottom).The pair truncation window  (;  0 = 200 ℎ −1 Mpc ) is applied to the theory correlation functions.The solid (dashed) lines are for  = 0.1 ℎ Mpc −1 ( = 0.3 ℎ Mpc −1 ); red (black) lines are with (with out) the distortion matrix.