The Alcock-Paczyński effect from Lyman-𝛼 forest correlations: Analysis validation with synthetic data

The three-dimensional distribution of the Ly 𝛼 forest has been extensively used to constrain cosmology through measurements of the baryon acoustic oscillations (BAO) scale. However, more cosmological information could be extracted from the full shapes of the Ly 𝛼 forest correlations through the Alcock-Paczyński (AP) effect. In this work, we prepare for a cosmological analysis of the full shape of the Ly 𝛼 forest correlations by studying synthetic data of the extended Baryon Oscillation Spectroscopic Survey (eBOSS). We use a set of one hundred eBOSS synthetic data sets in order to validate such an analysis. These mocks undergo the same analysis process as the real data. We perform a full-shape analysis on the mean of the correlation functions measured from the one hundred eBOSS realizations, and find that our model of the Ly 𝛼 correlations performs well on current data sets. We show that we are able to obtain an unbiased full-shape measurement of 𝐷 𝑀 / 𝐷 𝐻 ( 𝑧 eff ) , where 𝐷 𝑀 is the transverse comoving distance, 𝐷 𝐻 is the Hubble distance, and 𝑧 eff is the effective redshift of the measurement. We test the fit over a range of scales, and decide to use a minimum separation of 𝑟 min = 25 ℎ − 1 Mpc. We also study and discuss the impact of the main contaminants affecting Ly 𝛼 forest correlations, and give recommendations on how to perform such analysis with real data. While the final eBOSS Ly 𝛼 BAO analysis measured 𝐷 𝑀 / 𝐷 𝐻 ( 𝑧 eff = 2 . 33 ) with 4% statistical precision, a full-shape fit of the same correlations could provide a ∼ 2% measurement.


INTRODUCTION
Over the last few decades, probes of the large-scale structure (LSS) of the Universe have become one of the main tools for studying its expansion history and the properties of its constituents.Measuring the scale of the baryon acoustic oscillations (BAO) feature is currently the most widely used method for extracting cosmological information from LSS (e.g.Eisenstein et al. 2005;Cole et al. 2005;Beutler et al. 2011;Ross et al. 2015;Alam et al. 2017;Abbott et al. 2019;Alam et al. 2021;Abbott et al. 2022).BAO produces distinct features in the two-point statistics of LSS probes, which allows us to disentangle its signal from the other cosmological, astrophysical, and instrumental effects present.In BAO analyses, other cosmological information is effectively ignored by marginalizing over the broadband of the two-point statistics.
A common way of measuring this extra cosmological information ★ E-mail: cuceu.1@osu.edu is to instead model the full shape of the correlation function or power spectrum in order to measure the Alcock-Paczyński (AP) effect and redshift space distortions (RSD; e.g.Blake et al. 2011;Reid et al. 2012;Beutler et al. 2012;Samushia et al. 2014).However, such analyses are more difficult than BAO analyses, because accurately modelling these effects requires a good understanding of the astrophysical and instrumental contaminants present in the broadband of two-point statistics.Furthermore, linear theory approaches that may work well on large scales start to break down due to non-linear growth on smaller scales.These difficulties have been successfully addressed in the case of two-point statistics of discrete tracers, leading to numerous full-shape analyses in the literature (e.g.Blake et al. 2011;Reid et al. 2012;Beutler et al. 2012;Samushia et al. 2014;Satpathy et al. 2017;Beutler et al. 2017;Gil-Marín et al. 2020;Bautista et al. 2021;Hou et al. 2021;Neveux et al. 2020).On the other hand, in the case of the three-dimensional (3D) distribution of the Lyman- (Ly) forest, this type of measurement has not yet been performed (Cuceu et al. 2021).
The Ly forest is made up of absorption lines blueward of the Ly emission peak in spectra of high-redshift quasars (e.g.Lynds 1971;Rauch 1998).This absorption is caused by neutral hydrogen in the intergalactic medium between the quasar and us.Due to the expansion of the Universe, the Ly line is progressively redshifted as the photons travel towards us.This means the distribution of Ly absorption can be used to study the structure and evolution of the Universe.The forest has long been used as a tool for cosmology through its one-dimensional1 two-point statistics (e.g.McDonald et al. 2006;Seljak et al. 2005;Viel et al. 2010;Palanque-Delabrouille et al. 2015, 2020).With the advent of large spectroscopic surveys such as the Baryon Oscillation Spectroscopic Survey (BOSS) and its successor, extended BOSS (eBOSS), the 3D distribution of the Ly forest has also become a useful source of cosmological information (Busca et al. 2013;Slosar et al. 2013;Kirkby et al. 2013;Font-Ribera et al. 2014;Delubac et al. 2015;Bautista et al. 2017 The BOSS and eBOSS Ly cosmological analyses have progressively improved our understanding of the different processes that contribute to the auto-correlation function of Ly transmitted flux, and its cross-correlation with the quasar distribution.The first physical model of the two correlations was introduced with the BOSS data release (DR) 12 analyses (Bautista et al. 2017; du Mas des Bourboux et al. 2017).The main effects that had to be understood were the distortion caused by fitting the quasar continuum (see Bautista et al. 2017;Pérez-Ràfols et al. 2018 for the modern treatment and Slosar et al. 2011;Font-Ribera et al. 2012;Blomqvist et al. 2015 for further discussion of the effect), the contamination due to metal absorption in the forest (Pieri et al. 2010(Pieri et al. , 2014)), and the contamination due to high column density (HCD) absorbers (Font-Ribera et al. 2012;Font-Ribera & Miralda-Escudé 2012).Starting with the eBOSS DR14 analyses, Ly absorption in the Ly region (blueward of the Ly peak) was also used, first through its correlation with Ly absorption in the Ly region (de Sainte Agathe et al. 2019), and then through its cross-correlation with quasars (Blomqvist et al. 2019).The model of the Ly correlation functions has been progressively improved through these analyses, and we now may be in a position where we can perform a full-shape analysis.The goal of this article is to test this.In particular, we want to study the performance of this model when it comes to fitting the full-shape of the eBOSS DR16 Ly forest correlation functions.
As mentioned above, a common way of extracting cosmological information from the full shape of two-point statistics is to measure the AP effect and RSD.The AP effect arises due to the fiducial cosmology used to transform redshifts and angles into comoving coordinates.Any differences between this fiducial cosmology and the true cosmology will introduce an extra anisotropy in the autoand cross-correlations (Alcock & Paczyński 1979).Assuming we can accurately model other sources of anisotropy (e.g.RSD, HCDs, metals, distortion due to continuum fitting), the remaining anisotropy can be used to constrain the background cosmological model.This method is already used in analyses of the BAO, but these only use a small range of scales around the BAO peak.Cuceu et al. (2021) have shown that measuring the AP effect using the full shape is a promising avenue for extracting more cosmological information from the Ly forest.RSD measurements themselves are also generally used to measure the growth rate of cosmic structure, but for the Ly autocorrelation this is degenerate with an unknown RSD bias parameter (McDonald 2003;Slosar et al. 2011;Givans & Hirata 2020;Chen et al. 2021).Cuceu et al. (2021) did show that it is possible to measure the growth rate from the quasar RSD parameter in a joint analysis of the Ly auto-correlation and its cross-correlation with quasars.However, this would not be a significant measurement with eBOSS.Therefore, in this work we only focus on the AP information, and leave RSD measurements for future work using larger surveys such as the Dark Energy Spectroscopic Instrument (DESI Collaboration et al. 2016).
Our objective in this work is to analyse the full shape of the Ly forest correlation functions in simulated data in order to understand how well current models perform relative to the expected precision of the eBOSS DR16 data.This analysis is meant as a first step towards a full-shape analysis of the eBOSS Ly 3D correlation functions.The full-shape measurement using real eBOSS data is presented in a separate article (Cuceu et al. 2023).
This article is structured as follows.In Section 2, we introduce our framework, which includes a description of the mock data sets, the analysis process for computing the 3D correlation functions, and the modelling of the correlations.We show the results of our analysis in Section 3, and discuss their implications for a full-shape analysis using real data in Section 4. We summarize and conclude in Section 5.

FRAMEWORK
We begin with an overview of the framework we use to validate a full-shape analysis of the Ly forest 3D correlation functions.We use the set of one hundred eBOSS mock data sets produced by du Mas des Bourboux et al. (2020), which were created using the method introduced by Farr et al. (2020).We describe these mock data sets in Section 2.1, and the measurement of the transmitted flux correlations in Section 2.2.The analysis methodology we use is similar to the one used in past BAO analyses of Ly BOSS and eBOSS data.The only differences appear at the step of fitting the correlation functions.Instead of only fitting for the BAO peak position, we fit the full shape of the correlations in order to extract the AP parameter.This is based on the method introduced in Cuceu et al. (2021).We describe our modelling framework in Section 2.3.

Synthetic data-sets
In order to test the model of the Ly correlation functions, we use the set of one hundred mock data sets introduced in du Mas des Bourboux et al. (2020).These mocks were produced for the BAO analysis of the Ly eBOSS DR16 data set.Each mock is based on a Gaussian random field generated with the CoLoRe2 package (Ramírez-Pérez et al. 2022).The Gaussian field is transformed into a log-normal density field, which is Poisson sampled based on an input number density and bias in order to obtain a set of quasar sources.CoLoRe then computes line-of-sight skewers by interpolating the initial Gaussian field and the radial velocity field3 from each quasar to the centre of the box.
The skewers generated by CoLoRe require significant postprocessing in order to turn them into realistic simulations of the Ly forest.This is performed by the LyaCoLoRe4 package (Farr et al. 2020).In order to create the large boxes needed to simulate an all-sky light-cone to  = 3.7 (∼ 10 ℎ −1 Gpc in this case), the CoLoRe grid is limited by memory and computational constraints to a resolution of ∼ 2.4 ℎ −1 Mpc (for the 4906 3 grid used here).However, one of the contributions to the covariance of 3D correlations (and consequently of BAO uncertainties) is related to the amount of small-scale fluctuations in the Ly forest (McDonald & Eisenstein 2007), often characterized by the one-dimensional power spectrum, or  1D .In order to simulate spectra with a realistic  1D , LyaCoLoRe adds an extra 1D Gaussian field to each line-of-sight.After that, the field undergoes a log-normal transformation, and the output is used to compute the optical depth field () using the fluctuating Gunn-Peterson approximation (FGPA; Bi & Davidsen 1997;Croft et al. 1998).Redshift-space distortions (RSD) are introduced by convolving the real-space optical depth field with a kernel based on the peculiar velocity field.Finally, the transmitted flux in redshift-space () is given by  () = exp [− (𝑠)].For more details on this process, and how it is tuned to produce realistic forests, see Farr et al. (2020).
Beyond the Ly forest transmitted flux fraction, LyaCoLoRe also simulates high column density (HCD) absorbers, as well as absorption by other Lyman lines and metals.HCDs are sampled from the Gaussian density field using an input bias and number density based on literature results.After that, each HCD is allocated a column density given a column density distribution from the literature.In contrast, the other absorption lines are produced using a rescaled version of the Ly optical depth, which is then mapped to a different observed wavelength based on the rest-frame wavelength of the absorber.In the case of higher Lyman lines, the scaling factors are tuned based on the oscillator strengths of each transition (Iršič et al. 2013;Farr et al. 2020), with the scaling factor of Ly being 0.1901.For metal absorbers, the scaling factors are tuned in order to reproduce the level of contamination in the data (Farr et al. 2020;du Mas des Bourboux et al. 2020).A list of the main metal absorbers along with their relative strength can be found in Table 2 of Farr et al. (2020).These contaminants (HCDs, metal lines, and higher-order hydrogen lines) are then added to the simulated Ly forest transmitted flux fraction.
Each quasar is assigned a random magnitude using an input quasar luminosity function based on Palanque-Delabrouille et al. (2016).This is used to generate an unabsorbed continuum for each quasar by adding a set of emission lines on top of a broken power law (du Mas des Bourboux et al. 2020).Random redshift errors are drawn from a Gaussian distribution with a standard deviation   = 400 km s −1 .The specsim package (Kirkby et al. 2016) is then used to simulate the eBOSS spectral resolution, the pixelization of the spectra, and instrumental noise (du Mas des Bourboux et al. 2020).
While these synthetic data sets were created for BAO studies, in this article we use them to study full-shape analyses.As described above, the major contaminants that are known to affect Ly forest correlations are being modelled in these mocks.This will allow us to test the performance of our model, and how full-shape analyses are affected by these contaminants.However, one of the main concerns when attempting to use the full-shape of correlation functions for cosmology is the ability to accurately model the non-linear effects on small scales.While we will be able to test some of these using the mocks described above, some ingredients are missing.These include the fact that the small-scale quasar clustering is overestimated (Youles et al. 2022), and missing IGM effects such as pressure smoothing and the impact of UV radiation on the ionization fraction.We will discuss these issues and the applicability of our results in more detail in Section 4.

From spectra to correlation functions
Once the simulated eBOSS spectra are produced, they undergo the same analysis process as the real data, in order to measure the 3D auto-correlation of Ly transmitted flux, and its cross-correlation with the quasar distribution.This is described in detail in du Mas des Bourboux et al. (2020).Here, we only give a brief overview of this process, as all the results in this article5 come from the same mock correlation functions computed for the eBOSS Ly BAO analysis.
The first step in the analysis process is to compute the flux transmission fluctuation,   (), of each quasar , based on the observed flux,   (): where F () is the global mean transmission, and   () is the unabsorbed quasar continuum.In general, the true quasar continuum is unknown, and therefore, the data is used to jointly fit for the product F ()  ().This fit also necessarily includes density modes of the size of the forest and larger, which biases   towards zero for each forest.Even though pixels from the same forest are not cross-correlated, this will still bias pixel correlations in nearby forests, resulting in a distortion of the correlation functions.
The distortion of the correlation function can be removed by building a projection,     , such that: where    and    are the measured and true flux fluctuations,  indexes forest pixels before the projection, and  indexes the projected forest pixels.Therefore, the left-hand side of Equation ( 2) is applied to the measured flux fluctuation field, while the right-hand side is forward modelled in the correlations, as described in Section 2.3.A detailed description of this projection and how it is built can be found in Bautista et al. (2017), and Appendix B of Pérez-Ràfols et al. (2018).
Forests that contain HCDs with column densities log  Hi > 20.36 are given special treatment.Firstly, the regions where the HCD reduces the transmission by more than 20% are masked.Secondly, the absorption wings are corrected using a Voigt profile (du Mas des Bourboux et al. 2020;Noterdaeme et al. 2012).HCDs with column densities log  Hi < 20.3 remain in the data, and their effect has to be included in the model of the correlations (Section 2.3).
The correlation functions were computed on a grid in comoving coordinates.For two tracers,  and , at redshifts   and   , and separated by an angle Δ, the comoving separations along and across the line-of-sight are given by (de Sainte Agathe et al. 2019): where  M is the transverse comoving distance,  c () =  ∫  0  ′ / ( ′ ) is the radial comoving distance, with  () as the Hubble parameter and  as the speed of light.In this work, we will also refer to the (, ) parametrization, defined as The correlation function is first computed independently in each HEALPix8 pixel (Górski et al. 2005).For the eBOSS footprint, there are about 880 pixels (nside = 16), with each covering 3.7 × 3.7 = 13.4 deg 2 .This corresponds to a 250×250 (ℎ −1 Mpc) 2 patch at  eff = 2.33.The population of correlations can then be used to compute the mean and covariance of the correlation function of the entire survey.When computing the correlation function in one HEALPix pixel, pairs with forests in neighbouring pixels are still counted.However, given that we are most sensitive to small scales when measuring AP, we assume that the correlation function measurements in each HEALPix pixel are independent for the purposes of computing the covariance matrix.This method of computing the covariance matrix of Ly forest correlation functions has been validated against other methods by du Mas des Bourboux et al. (2020).
The process for computing the cross-correlation with the quasar distribution, and its covariance matrix, is similar to the one used for the auto-correlation.However, in this case we also distinguish between Ly flux in front of a quasar, which is assigned negative  | | , and Ly flux behind a quasar, which is assigned positive  | | .This is because the cross-correlation is not symmetric under permutations of the two tracers.
Four types of correlation functions were computed by du Mas des Bourboux et al. (2020) for the eBOSS DR16 Ly BAO analysis.The first two are the auto-correlation of the Ly flux in the Ly region, Ly(Ly)×Ly(Ly), and its cross-correlation with Ly flux in the Ly region, Ly(Ly)×Ly(Ly).The other two are cross-correlations between Ly forest flux and quasars, and include the cross-correlation of quasars with Ly flux in the Ly region: Ly(Ly)×QSO, and with Ly flux in the Ly region: Ly(Ly)×QSO.The Ly region is found between the Ly and Ly peaks, and is defined by the rest-frame interval  RF ∈ [104, 120] nm.The Ly region is blueward of the Ly peak, and is defined by the rest-frame interval  RF ∈ [92,102] nm.This information is summarized in Table 1.
In this work, we will be extensively using the mean and covariance of the one hundred mock eBOSS correlation functions.We will also refer to this mean as the stacked correlation function.In order to compute these, we first collect all individual correlation subsamples (in HEALPix pixels) from each of the one hundred mocks.The weighted mean of all correlation subsamples over all one hundred mocks gives us the stacked correlation, while its covariance is given by the covariance of all the subsamples.This stacked correlation was computed for each of the four correlation types.
We show the four stacked correlations in Figure 1, compressed into four  bins each.Besides the BAO peak at ∼ 100 ℎ −1 Mpc, the other features are due to metal contamination.This is why they are present in the wedges closer to the line-of-sight and absent from the ones across the line-of-sight.The most prominent of these is the blended metal peak due to SiII(1190) and SiII (1193), which are at comoving separations of  | | = 64 ℎ −1 Mpc and  | | = 56 ℎ −1 Mpc, respectively.In the auto-correlation, the SiIII(1207) peak can also be seen at  | | = 21 ℎ −1 Mpc.Finally, a fourth metal peak due to SiII( 1260) is present at  | | = 111 ℎ −1 Mpc, however, it is not visible due to its proximity to the BAO peak.

Modelling the correlations
Our model for the Ly forest auto-correlation (Ly×Ly) and its cross-correlation with quasars (Ly×QSO) closely follows that used by the eBOSS collaboration for the BAO analysis of the SDSS DR16 data (du Mas des Bourboux et al. 2020).The main difference is that we use scale parameters to model the full shape of the correlation, instead of restricting their application to scales around the BAO peak.Our model is based on an isotropic template power spectrum,  fid (), which is split into a peak (or wiggles) component,  p fid (), and a smooth (or no-wiggles) component,  s fid (), following Kirkby et al. (2013).As shown in Cuceu et al. (2021), this allows us to separate the information we obtain from the BAO peak from that obtained from the rest of the correlation.All the following steps in the analysis are done with both  p fid () and  s fid (), in parallel.We use the Vega package to perform this analysis. 9he Ly auto and cross power spectra are given by: where  ′ Ly and  ′ Ly are the effective Ly bias and RSD parameter, and  QSO and  QSO are the quasar linear bias and RSD parameter.The Ly absorption is given by a combination of contributions from the diffuse IGM and from high column density (HCD) systems.Both trace the underlying density field, but HCDs correspond to regions with significantly higher clustering.Here, we define HCDs as systems with neutral hydrogen column density above 10 17.2 cm −2 , which means they include both Lyman limit systems and damped Ly (DLA) systems.As described in Section 2.2, large HCDs are detected and masked.However, those with widths smaller than ∼ 14 ℎ −1 Mpc cannot be detected and remain in the data.These remaining HCDs lead to a broadening effect along the line of sight.Font-Ribera & Miralda-Escudé (2012) showed that this can be modelled through an extra  | | dependent term in the effective bias and RSD parameters.This is given by: where  Ly and  Ly are the linear bias and RSD parameters associated with the IGM and  HCD and  HCD are the linear bias and RSD parameters associated with HCDs.The  NL term models the quasar non-linear velocities and statistical quasar redshift error.Following Percival & White (2009), we test two different models for this term, one that introduces a Lorentzian damping, and one with Gaussian damping.They are given by: where   is a free parameter.Given that the synthetic data is created with Gaussian redshift errors, the Gaussian damping model is more appropriate for modelling the mocks.However, when modelling the real data, the Lorentz damping model is used (du Mas des Bourboux et al. 2017Bourboux et al. , 2020)).This is because the quasar redshift errors generally have long tails that are not well modelled by a simple Gaussian (e.g.Lyke et al. 2020).Therefore, for the purposes of this analysis, we will test both models.We also introduce Gaussian anisotropic smoothing,  sm (,   ), to account for the low-resolution of the CoLoRe grid, following Farr et al. (2020).This smoothing model has two free parameters,  | | and  ⊥ , which describe the smoothing along and across the lineof-sight, respectively.They are expected to have values of ∼ 2.4 ℎ −1 Mpc, corresponding to the resolution of the CoLoRe grid.
In order to turn model power spectra into correlation functions, we follow the process described in Kirkby et al. (2013).This involves a multipole decomposition, followed by a Hankel transform10 to turn each power spectrum multipole into the corresponding correlation function multipole, and finally the reconstruction of the twodimensional correlation function from the individual multipoles.Following past Ly BAO analyses, we use multipole values up to ℓ = 6.The main components of the theoretical model are the Ly transmitted flux correlations:  Ly×Ly for the auto-correlation, and  Ly×QSO for the cross-correlation with quasars.They are computed from the two power spectra in Equations ( 5) and ( 6), respectively.
We also model the contamination due to metal absorption, using the same models as for the Ly auto and cross-correlation.We compute models for the correlations between each metal line and Ly ( Ly× ), and among all metal pairs (  1 × 2 ) for the autocorrelation, and between each metal line and quasars ( QSO× ) for the cross-correlation.This is described in more detail in Appendix A. Each metal line has its own bias and RSD parameter (  ,   ).The same Ly and QSO (, ) parameters above are also used for the cross-correlations between metals and Ly and quasars.In this case, we neglect the HCD effects for the Ly parameters.Following past Ly BAO analyses, we fix the metal  parameters to 0.5 (Bautista et al. 2017;du Mas des Bourboux et al. 2020).
The full model correlations are then given by: where the sums are performed over the four metal lines introduced in Section 2.2.The distortion due to quasar continuum errors is modelled using distortion matrices that multiply the results of Equations ( 11) and ( 12).See Appendix A and du Mas des Bourboux et al.
As described at the start of this section, we perform all these steps using both isotropic power spectrum components,  p fid () and  s fid (), in parallel.Therefore, we compute the projected correlation function for both the peak component, ξp , and the smooth component, ξs .The two model correlations are then combined into the final theoretical model, allowing for their comoving coordinates to vary: In the case of the cross-correlation, we also allow for a systematic shift in the  | | coordinate of the cross-correlation through a nuisance parameter, Δ | | .This means the coordinate transform is given by We introduce this parameter in order to find if it has any impact on the scale parameters of interest, which could lead to a potential bias when performing the analysis on real data, as it was found to be non-zero in past Ly BAO analyses (e.g.This shows that the two effects (AP and RSD) lead to different changes in the correlation as a function of , which explains why they are not degenerate.

Modelling the Alcock-Paczyński effect
In order to isolate the Alcock-Paczyński effect, we transform the scale parameter system, following Cuceu et al. (2021): where  is the anisotropic parameter, and  is the isotropic scale parameter.Cuceu et al. (2021) showed that a measurement of  is equivalent to a measurement of the AP parameter   () ().As we have two sets of ( | | , ⊥ ) parameters, we also have two sets of (, ) parameters, for the peak and the smooth component, respectively.The two peak scale parameters have already been tested and measured in BAO analyses.Therefore, our goal here is to use the simulated data sets described above to study measurements of the anisotropic parameter of the smooth component,  s .
The two main types of analyses we will perform are described in Table 2.In the first one, which we will refer to as the "split AP measurement", we fit two distinct  parameters for the peak and smooth components separately.In the second one, called "full AP measurement", we fit one  parameter to the full shape of the correlations.When performing this analysis on real data, we want to use the full AP measurement, because it coherently extracts the desired information.The split AP measurement can be used as a consistency check of the method or data.However, it is also useful in the current work, because  s quantifies the gain in information over a standard BAO-only analysis.
In the top panel of Figure 2, we show the impact of changing  on the correlation function compressed in a shell in isotropic separation and plotted as a function of .We choose the smallest separation shell, 25 <  < 45 ℎ −1 Mpc, in order to best illustrate the anisotropy pattern we aim to recover when fitting the AP effect.The bottom panel of Figure 2 shows the impact of changing the Ly forest RSD parameter  Ly , which also affects the anisotropy of the correlation.This figure shows that AP and RSD produce different anisotropy patterns as a function of , which means a 2D measurement of the correlation (i.e., in  | | ,  ⊥ or , ) can be used to disentangle the two effects.
The isotropic scale parameter from the smooth component,  s , also contains some cosmological information due to the scale of matterradiation equality.However, as described in Cuceu et al. (2021), it is not clear how to extract this information given that the effect is very similar to that produced by the Ly bias, and likely to be heavily impacted by contaminants.Furthermore, Cuceu et al. (2021) found that the constraint on the isotropic scale from the peak,  p , is much tighter than the one from the smooth component.This is in contrast to , where the constraint from the broadband is the one that dominates (see Figure 2 of Cuceu et al. 2021).Finally, Gerardi et al. (2023) recently performed simplified forecasts of a direct cosmological analysis of the Ly forest correlation functions and found that the information extracted can be traced back to the Alcock-Paczyński effect and redshift space distortions.Therefore, for the purposes of this analysis, we will treat  s as a nuisance parameter, and marginalize over it.
In Table 3, we give a description of all the free parameters in our standard analysis, along with their priors.We follow past BAO analyses, and assign a Gaussian prior to the RSD parameter of HCDs, based on measurements of DLA clustering (Pérez-Ràfols et al. 2018).
We plot the joint best fit model using the split AP analysis configuration in Figure 1.We fit the correlation functions between  min = 25 ℎ −1 Mpc and  max = 180 ℎ −1 Mpc following the analysis in Section 3.2.11This model works surprisingly well given the extreme circumstances of this test.While the fit statistics clearly indicate the model is not good enough for fitting the stack of 100 eBOSS mocks (the fit probability is ∼ 10 −11 ), the contrast with the eBOSS uncertainties, shown as shaded bands, visually indicates that our model could be good enough for the statistical precision achievable with current data sets.We will show this using the stack of correlations in Section 3.1 and by analysing individual eBOSS mocks in Section 3.4.

RESULTS
A large part of the validation of our analysis will focus on studying full-shape fits using a stack of the Ly forest correlation functions from the eBOSS mocks.Fitting the stack of the mock correlations is much less computationally expensive compared to separately fitting all one hundred mock correlations.Therefore, this method is advantageous when we want to test the performance of our model or study different versions of the analysis.This is because it allows us to study the trends necessary without as much of an emphasis on the effect of noise.Our main goal in this section is to test if the current modelling approach for the Ly forest correlation functions is appropriate for full-shape analyses of the eBOSS data.In order to achieve this, we want to gain an understanding of the relation between the statistical precision we expect using the eBOSS data and the systematic bias.We also want to investigate modelling choices such as using one versus two anisotropic parameters, the relative contributions of the auto and cross-correlations, and how our modelling of the contaminants affects our measurements.
though that metal peak is outside our fitted range.This is because its contamination is spread all along the line-of-sight by the distortion matrix, and therefore we have to marginalize over it.

Analysis of stacked correlations
As we have four different correlations, we categorise them into the two auto-correlations of Ly flux, Ly(Ly)×Ly(Ly) and Ly(Ly)×Ly(Ly), and the two cross-correlations with quasars, Ly(Ly)×QSO and Ly(Ly)×QSO.This is because the modelling of the two auto-correlations is identical (and similarly for the two cross-correlations) as we treat the Ly flux in the Ly region the same as the Ly flux in the Ly region.In the case of the two cross-correlations, the full system of parameters is degenerate.This degeneracy can only be broken when performing a joint analysis of the auto and cross-correlations, because the auto-correlation helps constrain the Ly bias and RSD parameters.Therefore, we fix the bias and RSD parameter of quasars ( QSO ,  QSO ) to a fiducial value when fitting the cross-correlation alone.12Note that due to correlations between the AP effect and RSD, this may lead the uncertainty in  to be underestimated.However, we will not use these  constraints from the cross-correlation to make any decision or recommendation in this article.We only compute and show them in order to check consistency and for completeness.Also note that, in general, the parameter constraints from the stack of the correlations will not be the same as the mean of the constraints from individual mocks.This is due to the non-linear dependence of the likelihood on the parameters.However, the stacked correlations are useful for testing the performance of our model.
We assume there is no cross-covariance between any of the four correlations, as du Mas des Bourboux et al. (2020) showed that these are negligible.We use a Gaussian likelihood, and compute posterior distributions using the Vega package for modelling and the PolyChord13 package for sampling (Handley et al. 2015a,b).When running PolyChord, we initialize a number of live points equal to 25 times the number of parameters, and set the length of the slice sampling chains to be of order of the number of parameters.
In Figure 3, we show the results for   ,   , and   from the joint fit of the two auto-correlations alone, the two-cross correlations alone, and the combination of all four correlations, using the split AP analysis.These fits use a minimum separation  >  min = 25 ℎ −1 Mpc following the analysis in Section 3.2.We find good agreement between the scale parameters measured from the auto-and cross-correlations.
We use the same cosmological model that was used to construct the mocks to also compute comoving coordinates and the power spectrum template.14Therefore, we expect to recover unity for all (, ) parameters.We do recover the expected  values from both the auto and cross-correlation, using both the peak and smooth components independently.The fact that we are able to constrain  s means it is not degenerate with any of the other effects that introduce their own anisotropies.For a discussion of these effects see Section 3.3, and for plots of correlations between  s and important nuisance parameters see Appendix E. Finally, Figure 3 shows that we are indeed able to perform an unbiased analysis of the full-shape of Ly forest correlations for mocks with known contaminants.
For the isotropic scale of the peak,  p , the result from the cross is consistent with unity.However, the result from the auto-correlation is higher, leading to the combined constraint being ∼ 0.22% larger than unity.While this appears as a ∼ 2 systematic bias in the analysis on the stack of 100 eBOSS mocks, it is a small fraction (∼ 0.2) of  , is denoted as "Auto", and the one from the two cross-correlations with quasars, Ly(Ly)×QSO and Ly(Ly)×QSO, is denoted as "Cross".All scale parameters are expected to be unity, given that the same cosmology was used to construct the mocks and to analyse the data.We also show the expected eBOSS constraint for comparison, but note that the 1 2D contours are all larger than the bounds of the plots here.
We find that we are able to recover an unbiased measurement of the AP parameter from the full shape of the Ly forest correlations.
the expected eBOSS statistical uncertainty.We will return to discuss this deviation in Sections 3.4 and 4.
While for the BAO peak parameters,  p and  p the constraints obtained from the auto and cross-correlations are very similar, this is not the case for  s , where the measurement from the auto-correlation is ∼ 30% tighter than that from the cross-correlation.This is in line with the forecasts by Cuceu et al. (2021).We also find that adding the correlations of Ly flux in the Ly forest shrinks the  s constraint by 16% for the auto-correlation and by 10% for the cross-correlation.In Figure 3 we only focus on the parameters of interest for cosmology.However, unlike BAO parameters,  s does have correlations with some nuisance parameters, the most important being  Ly ,  Ly and  HCD .We show these in Appendix E.
In Figure 4, we compare the split AP measurement with the full AP measurement, based on the joint analysis of all four stacked correlations.We show 1D marginalised posteriors of  measurements from the peak and smooth components ( p and  s ) and the measurement from the full shape (  ).We find that the measurements from the peak and smooth components are consistent with each other, with the smooth component providing a constraint that is ∼ 32% tighter than the constraint from the peak component.Note that we do not expect the  p and  s measurements to be the same, even though they both   measure the same effect.This is because they use different parts of the data, resulting in two mostly independent measurements.This is confirmed by the lack of correlation between the two parameters, as shown by their 2D posteriors in Figure 3.The measurement of   is remarkably consistent (to within 0.1) with the expected  = 1.Overall, measuring  from the full-shape of the correlations leads to a ∼ 41% improvement in the 1 constraints compared to the BAO-only analysis.These numbers will depend on the range of scales we include in our fits, so we next turn our attention to studying the effects of the minimum scale cut we employ.

Statistical and systematic errors
The results we have shown so far are already quite promising, with   and   measurements that are consistent with the true cosmology of the mocks ( = 1).The next step is to quantify this statement.
To achieve this, we aim to study how the statistical precision and possible systematic biases depend on the minimum scale included in our analysis.Including more data (especially at small separations) can significantly improve our constraining power (Cuceu et al. 2021), but it may also introduce systematic biases due to imperfect modelling (of e.g.non-linear effects, HCDs, etc.).Therefore, our goal is to find a balance between these two.The key variable in this analysis is the minimum separation scale included in the fit,  min .
In order to quantify the systematic bias of our  measurements, we compute the difference between the best fit value of  from the stacked correlations and the true value in the mocks, i.e., Δ =  best − 1.While we do have the size of the statistical constraints from the stack of the correlations, we want to compare Δ with the expected constraints from one eBOSS realization.Therefore, we rescale the statistical uncertainty in the fit of the stack by a factor of 10 (given the 100 mocks used to compute the stack) to obtain the expected eBOSS constraint,   .
We compare the statistical constraint and the systematic error on  for different values of  min in Figure 5.We use the full AP setup and show results for   from the joint analysis of all four correlation functions.The error-bars on the systematic errors are given by the 1 statistical constraints of each measurement, which represents the constraining power of the stack of the one hundred mocks.
As expected,   grows significantly as we cut more data (larger  min ), with values ranging from 1.6% to 2.7%.The systematic bias starts high for  min = 10 ℎ −1 Mpc (0.7%), but becomes consistent with zero by  min = 25 ℎ −1 Mpc.Note that we initially performed this test in  min steps of 10 ℎ −1 Mpc, but chose to compute the  min = 25 ℎ −1 Mpc point after seeing the initial trends.We also performed the same analysis on  s , using the split AP approach, and found similar results (systematic bias consistent with zero for  min ≥ 25 ℎ −1 Mpc).
These results confirm that we are able to obtain unbiased measurements of  using the full shapes of the Ly forest correlation functions in analyses using  min ≥ 25 ℎ −1 Mpc.Therefore, throughout the rest of this work, we will be using  min = 25 ℎ −1 Mpc as our baseline analysis.We also studied the impact of other types of scale cuts, as described in Appendix C. Based on that, we decided to follow BAO analyses and use a maximum scale of  max = 180 ℎ −1 Mpc, and not impose any anisotropic cuts.

Modelling of contaminants
The results we have seen so far indicate that our current modelling of the relevant contaminants in the Ly forest correlations is good enough to recover unbiased full-shape measurements of the AP effect with eBOSS.Before moving to the study of analyses on individual eBOSS mocks, we wish to take a closer look at the individual components of our model.In particular, we want to understand how our measurements are affected by these individual models, and how sensitive we are to different analyses choices.Of particular interest for AP are those effects that introduce their own anisotropies.
The most important effect we have to account for is that of the distortion due to continuum fitting.This effect is modelled by introducing a projection (Equation 2; Bautista et al. 2017) which removes the affected modes in both our data and our models.These distortions are limited to large scales in the power spectrum (low values of  ∥ ), but affect the correlation function on all scales (Blomqvist et al. 2015;Bautista et al. 2017).Given the comparison in Figure 1 between the stacked correlations and the best fit model, combined with the fact that we can recover unbiased values of  s and   (Figure 4), we conclude that our modelling of this distortion is accurate enough for current datasets.
While large high column density (HCD) absorbers are masked in the input spectra, those with a typical width smaller than ∼ 14 ℎ −1 Mpc remain (de Sainte Agathe et al. 2019;du Mas des Bourboux et al. 2020).The current model of this effect (Equations 7 and 8) follows the description introduced by Font-Ribera & Miralda-Escudé (2012), combined with an exponential shape for the  HCD ( | | ) function (based on the analysis by Rogers et al. 2018).Note that this shape assumes a model for the column density distribution function of HCDs based on observations (Noterdaeme et al. 2009;Prochaska et al. 2010;Noterdaeme et al. 2012;Zafar et al. 2013), which are very limited for column densities log  Hi < 20.This model takes into account the joint contribution of IGM and HCD absorption at the level of two-point functions (which includes the Ly and HCD auto-correlations and their cross-correlation).However, Font-Ribera & Miralda-Escudé (2012) found that the three-point function, ⟨ Ly  HCD  Ly ⟩, also has a significant contribution.Furthermore, the process of masking the larger HCDs could also bias our correlation measurements, because we are preferentially masking regions with high clustering.These two potential sources of systematic bias have been long known (Slosar et al. 2011;Font-Ribera & Miralda-Escudé 2012).However, so far they have not had a significant impact.This appears to still be the case with full-shape analyses of eBOSS, taking  min = 25 ℎ −1 Mpc.However, they could be contributing to the bias we observe for  min < 25 ℎ −1 Mpc.We have also tested the impact of marginalizing over the  HCD parameter (instead of fixing it to 10 ℎ −1 Mpc), and found that it has no impact on  s measurements.
Contamination due to metal absorption results in clear extra features in the Ly correlation functions along the line-of-sight.These features correspond to correlations at vanishing separation between metal absorption and either Ly absorption (auto) or quasars (cross).These correlations get assigned to the wrong bins in comoving coordinates because we interpret everything as Ly absorption.This is a simple enough process that we can replicate it in the model through the metal matrix formalism described in Appendix A. However, it means that correlations that used to be at small separations are now spread across the correlation along the line-of-sight.This could result in potential systematic biases, as imperfect modelling of the small scales can now have an impact even on large scales.The results shown above do not indicate this is a problem for eBOSS.
We also note that unlike HCDs, metals were added to these mocks using a re-scaled version of the Ly optical depth.In reality, metal absorption is associated with the circum-and intra-galactic medium more than the IGM (Pérez-Ràfols et al. 2022).Therefore, these mocks only provide a rough approximation for the clustering of metals.This is especially relevant in the case of metal RSD, as the model has metal RSD parameters fixed to  m = 0.5, following du Mas des Bourboux et al. (2020).Synthetic data sets with more realistic metal prescriptions are needed to fully understand their impact on AP measurements.
Quasar redshift errors can be a significant source of uncertainty and also potentially bias both BAO and AP measurements (Youles et al. 2022).For the auto-correlation, quasar redshifts are only used to define the rest-frame wavelength range when computing the continuum.Therefore, redshift errors result in extra spectral diversity, which increases the noise in the measurement of the flux transmission field (du Mas des Bourboux et al. 2020).On the other hand, the cross-correlation is smeared along the line-of-sight, and a potential systematic offset between negative and positive  | | is introduced.We model the first effect using either a Gaussian or a Lorentzian model with a free parameter, as described in Section 2.3.The second effect is modelled through an extra free parameter, Δ | | , to account for a potential offset.Furthermore, the Gaussian or Lorentzian models are also meant to model the non-linear quasar velocities which give rise to the Finger-of-God effect.The mocks used here do contain statistical redshift errors, modelled as a Gaussian, and the Finger-of-God effect.As analyses on real data have preferred the Lorentzian model, we tested both of them.However, we found no significant shift in any of the scale parameters of interest ( p ,  p ,  s ).
In Figure 6, we compare the baseline results on  with results where we remove the model for each of these contaminants in turn.This means that the contaminants are still present in the measured correlations, but we do not model them.These results are computed from a joint analysis of all four Ly stacked correlation functions.We also show the baseline result rescaled to correspond to one eBOSS realization.We study three cases: one where we remove the HCD model, one where we do not model metals, and one where we do not model quasar redshift errors.All three introduce some systematic bias in the measurement of  s , with the largest resulting from the removal of the HCD model, leading to a ∼ 1% shift.On the other hand, for  p , only the removal of the metal modelling results in a significant bias.This figure shows the relative importance of including each of these contaminants in the model.However, note that even the largest bias (with no HCD modelling at all) represents a shift of only ∼ 0.5 of the expected statistical constraint for eBOSS.

Analysis of individual mock realizations
We now turn our attention to analysing individual eBOSS mocks.The analysis of the stacked correlations has been useful for understanding how well our model performs in the context of a full-shape analysis.We found that, using a minimum separation  min = 25 ℎ −1 Mpc, we are able to obtain unbiased measurements of the anisotropic scale parameter from the full shape of the correlations.However, as mentioned above, parameter constraints from the stack of mocks are not, in general, the same as mean constraints from individual mock Posterior distributions of the anisotropic scale parameter, , from the peak and smooth components, respectively.We compare the baseline result, where we model all contaminants following eBOSS (du Mas des Bourboux et al. 2020), with results where we do not model each one of the contaminants in turn.The three results here are obtained by either removing the HCD modelling, the modelling of the metals or the model for quasar non-linear velocities.In each of these cases, not modelling the relevant contaminant introduces a systematic bias in our results.For  s , the most significant shift appears when removing the model for HCDs.
analyses.Therefore, we now wish to test if our conclusions above are also accurate when analysing individual eBOSS mocks.The results shown so far have been obtained by using a sampler to compute full posterior distributions.However, this becomes too computationally expensive for a full set of 100 mocks. 15Therefore, in this section, we compute the best fit values of the parameters using a minimizer 16 .We estimate the statistical uncertainty in this result using the second derivative in parameter space around the best fit point, assuming a Gaussian likelihood.This allows us to obtain rough estimates of population statistics from the set of mocks.However, note that both the best fits and the constraints obtained using this method are ultimately just an approximation, because we are not properly marginalizing over the nuisance parameters (see Cuceu et al. 2020).We only use this method here because a more complex and accurate method (i.e.sampling) is not practical.We discuss this in more detail in Appendix D.
Our modelling of individual eBOSS mock correlations closely follows the model we used for the stack of correlations, as introduced in Section 2.3.We run the iminuit minimizer on each of the one hundred mocks.This provides us with both the best fit value of each parameter and with its Gaussian uncertainty.We compute two versions of the results, using both the split AP analysis with  p and  s , and the full AP analysis with   .
We show histograms of the best fit values of the scale parameters 15 Running the sampler for a joint analysis of all four Ly correlations with all the contaminants requires ∼ 2 − 3 × 10 3 CPU hours, which leads to a rough estimate of ∼ 2 − 3 × 10 5 CPU hours to analyse all 100 mocks. 16We use the iminuit package (Dembinski et al. 2021) for minimization.The top plots show fits of  p (left) and  s (right).These show that we are able to obtain unbiased measurement of  from the full shape of the Ly forest correlations, which is confirmed by the bottom right plot, showing results for   .The bottom left plot shows results from fits of the BAO isotropic scale parameter  p .The mean of the population in this case is slightly higher (by ∼ 0.25) than the expected value.
in Figure 7.The top plots show results for  p and  s , both of which are consistent with the expected value of unity.The means of the populations are   p = 1.0028 and   s = 1.0013, corresponding to shifts from unity of ∼ 0.06 and ∼ 0.05, respectively.17These results reinforce the conclusions from Section 3.1: we are indeed able to recover unbiased measurements of the AP effect using the full shapes of the correlations.The bottom row of Figure 7 shows results for the isotropic BAO scale ( p ) and the full-shape anisotropic parameter (  ).The means of the populations in this case are   p = 1.0024 and    = 1.0024, corresponding with shifts from unity of ∼ 0.23 and ∼ 0.11, respectively.While these shifts are more significant compared to the ones from  p and  s , they are still relatively small, in line with the results from Section 3.1.

DISCUSSION
In this work, we have used a set of one hundred eBOSS Ly forest mocks to study how well we could perform a full-shape analysis of the eBOSS DR16 Ly correlation functions.We found that we are able to recover unbiased measurements of the anisotropic scale parameter, , from a joint analysis of the full shapes of the four Ly correlations when using a minimum separation of  min = 25 ℎ −1 Mpc.The only significant systematic biases on  were found when either we did not model one of the major contaminants (Section 3.3), or we included scales smaller than 25 ℎ −1 Mpc (Section 3.2).
We did find an apparent bias of 0.22% − 0.24% in the isotropic scale of the peak component,  p , both when analysing the stack of one hundred mocks (Section 3.1) and when analysing them individually (Section 3.4).This corresponds to a ∼ 0.21 − 0.25 bias for eBOSS.18du Mas des Bourboux et al. ( 2020) also found a systematic shift of similar magnitude for the BAO parameters, using the same correlation function measurements with a different parametrization.In the same publication, du Mas des Bourboux et al. ( 2020) measured BAO using real data, and tested the robustness of the result by adding broadband polynomials meant to marginalize over any unaccounted broadband systematic effects.They found no significant shift in the BAO parameters when performing this test.Therefore, they concluded the BAO results are robust.Such a test cannot be replicated in the context of a full-shape analysis, because the broadband polynomials would erase the information in  s .However, our model allows the separate fit of BAO ( p and  p ) and broadband ( s ) information.Therefore, when analysing real data, we can benchmark our BAO result using the results of the eBOSS BAO analysis.
In Section 3.3, we have discussed the impact of the major contaminants that are modelled in the mock data and shown the importance of modelling them accurately.However, we have only touched on the effects that we model and observe in the mock data.While these represent the most important known contaminants for the Ly forest 3D correlations (which is why they are modelled in the mock data), they are not the only possible sources of systematic uncertainty.Here we wish to discuss some other possible sources of contamination, how they have been treated in the past, and also to give recommendations on which of them warrant attention when performing this analysis on real data.
The mock data used in this work is based on a Gaussian density field with quasars randomly sampled using its log-normal transformation.While this is a good enough approximation on large scales, it does not correctly reproduce the quasar clustering on small scales (Farr et al. 2020).This makes it difficult to accurately test the model on small scales, especially for the cross-correlation.The problem is also exacerbated by the fact that both the distortion due to continuum fitting and the metal contamination take effects that are at small scales and spread them to larger scales along the line-of-sight.Incorrect modelling of the small scales can thus have an impact on all scales, so this analysis would benefit from more realistic mocks to test this effect.However, as these are not currently available, we have to leave such a study for future work.
For the analysis of the Ly auto-correlation from real data, an empirical model based on the work by Arinyo-i-Prats et al. ( 2015) is used to fit small-scale non-linearities. 19This model has five free parameters whose values were measured using hydro-simulations.However, no equivalent model exists for the cross-correlation.A possible approximation would be to use the same model for the cross-correlation as well, but it is not clear if the same parameter values would still work, and this would not account for quasar non-linearities (FoG).Recently, Givans et al. (2022) used high-resolution hydrodynamic simulations to study how well current models perform on small scales in both the Ly power spectrum and the cross-spectrum of Ly with massive halos.For Ly, they found that the model by Arinyo-i-Prats et al. (2015) performs very well up to wavenumbers much larger than those included here.On the other hand, for the cross-spectrum, they found that the standard Lorentzian damping combined with the Arinyo-i- Prats et al. (2015) model cannot accurately fit the smallscale power.Therefore, for the cross-correlation, they recommend using scales larger than 30 ℎ −1 Mpc when fitting it with linear theory.However, they did not include the effect of redshift errors, which suppresses power on these small scales, making deviations from linear theory less important.Therefore, in an analysis of real data, it may be worth testing scale cuts for the cross-correlation of both 25 ℎ −1 Mpc, as proposed here, and of 30 ℎ −1 Mpc, as proposed by Givans et al. (2022).
Comparing results from the auto and cross-correlations could be a good way of testing for systematic biases due to imperfect modelling of the small scales.This is because the two are driven by different effects, with the cross-correlation affected more by redshift errors and quasars, which cluster very strongly and have large non-linear peculiar velocities, and the auto-correlation driven by the IGM absorption (gas pressure effects, thermal broadening, etc.).Furthermore, for the cross-correlation, testing both the Gaussian and Lorentzian damping models would be a good consistency check for the analysis on real data.
While preparing this manuscript, a new study by Youles et al. (2022) found that large redshift errors can introduce spurious correlations along the line-of-sight in both the Ly auto and crosscorrelation with quasars.Such an effect could lead to systematic biases when performing a full-shape analysis such as the one studied here.Therefore, future work is needed to study its impact on fullshape analyses, and to potentially find ways to account for it (e.g. by forward modelling the effect).
An effect that is not included in the mock data is that of BAO broadening due to non-linear growth.However, this effect is relatively well understood and modelled in BAO analyses (see Kirkby et al. 2013 for a detailed discussion).Therefore, as long as it is correctly modelled when fitting the real data, we do not consider it a potential issue for a full-shape analysis.The effect of quasar radiation, also known as the transverse proximity effect, is another important source of contamination for the cross-correlation, because it increases the ionization fraction in the gas surrounding the quasar, leading to a decrease in Ly absorption.However, this effect has been successfully modelled (Font-Ribera et al. 2013), and this model has been included in previous BAO analyses.As this only affects the cross-correlation, comparing results from the auto-and cross-correlation should also be useful here.
The metal contamination introduced in the synthetic data used here only includes the effect of the lines that are close in wavelength to the Ly peak.This type of contamination is mostly given by the crosscorrelation between these absorptions lines and either Ly absorption (for the auto-correlation) or quasars (for the cross-correlation).However, in BOSS and eBOSS BAO analyses, the contamination due to CIV absorption was also modelled.At the relevant comoving coordinates used here, this is dominated by the auto-correlation of CIV absorption, and therefore only affects the Ly auto-correlation.It is generally modelled in the same way as the other metal lines, as described in Section 2.3.It has its own free bias parameter.However, this parameter has only been marginally detected, with a significance less than 3 (du Mas des Bourboux et al. 2020).Therefore, the impact of this contamination is minimal.Nevertheless, we recommend that it be modelled and marginalized over in a full-shape analysis on real data, similarly to Ly BAO analyses.
Fluctuations of the ionizing UV background introduce a scale dependence to the Ly bias and RSD parameters (Pontzen 2014;Gontcho A Gontcho et al. 2014).This scale dependence was modelled in the Ly BAO analyses using BOSS DR12 and eBOSS DR14 (Bautista et al. 2017;de Sainte Agathe et al. 2019;Blomqvist et al. 2019), following the framework introduced by Gontcho A Gontcho et al. (2014).However, the effect was not detected at a significant level (only ∼ 2), and modelling it did not have an impact on the BAO fits.This indicates that it is probably not of significant concern for a full-shape analysis, but its impact on  measurements should be tested when performing this analysis on real data.
In conclusion, the two effects that likely require more attention are the spurious correlations due to redshift errors (identified by Youles et al. 2022), and the small scale non-linear effects in the crosscorrelation.The other effects discussed here have been studied before, and modelled at some level in past Ly BAO analysis.They were not prioritized by the eBOSS collaboration when creating mocks because they were not found to play an important role in past BAO analyses.However, this may not be true for full-shape analyses.Therefore, an important part of the first full-shape analysis on real data will be to determine if the current models for these secondary effects have an impact on measurements of  s .If they are found to have an impact, it will motivate the development of more accurate mocks in the future.

SUMMARY
The Lyman- (Ly) forest provides one of the best tracers of large scale structure (LSS) at high redshifts (2 <  < 4).However, 3D correlations of the forest have so far only been used to measure the baryon acoustic oscillations (BAO) scale.Recently, Cuceu et al. (2021) showed that a measurement of the Alcock-Paczyński (AP) effect from the full shapes of Ly forest correlation functions can lead to significant improvements in cosmological constraints compared to BAO-only measurements.In this work, we use synthetic data of the extended Baryon Oscillation Spectroscopic Survey (eBOSS) in order to test such an analysis.
We model the Ly auto-correlation, and its cross-correlation with quasars, using a similar approach to that used in past Ly BAO analyses (Kirkby et al. 2013;Bautista et al. 2017;du Mas des Bourboux et al. 2017, 2020), as described in Section 2.3.The only difference appears from the fact that we use scale parameters to fit the full shapes of the correlations.Our parametrization of the scale parameters splits them into an isotropic scale parameter, , and a parameter for the anisotropy, , corresponding to a measurement of the AP effect (Cuceu et al. 2021).As we decompose the template power spectrum into a peak component and a smooth component, we have the option of using separate scale parameters for each.We use two setups for the analysis, as described in Table 2.In the first, we decouple the peak and smooth components by fitting two distinct  parameters ( p and  s ), while in the second we fit only one  parameter (  ) for the full shapes of the correlations.
We begin our study by analysing stacked correlation functions from the one hundred mock data sets (see Figure 1), in Section 3.1.We perform a joint analysis of all four correlations, as shown in Figures 3  and 4. We find that the results from the two Ly auto-correlations are consistent with the results from the two cross-correlations.Furthermore, we find that we are able to recover unbiased measurements of  from the full shapes of the correlations, for both types of analysis described above.
In Section 3.2, we study joint analyses of all four stacked correlations using a range of scale cuts.In particular, we focus on the minimum scale used,  min .We run the full AP analysis using  min values from 10 ℎ −1 Mpc to 50 ℎ −1 Mpc, as shown in Figure 5.We then compare the systematic bias of   measurements, given by the difference between the best fit value and the expected value in the mocks ( = 1), with an estimate of the expected statistical constraint from eBOSS.We find that we are able to obtain unbiased measurements of   for  min ≥ 25 ℎ −1 Mpc.Therefore, we use  min = 25 ℎ −1 Mpc for the other parts of the analysis.
In Section 3.3, we discuss the main contaminants present in the Ly forest correlation functions.We also test a few different setups for the analysis, but find that as long as the main contaminants are modelled correctly, there is no significant bias in full-shape  measurements.We show the relative impact of not modelling each main contaminant in Figure 6.
We study full-shape analyses on the individual correlations in Section 3.4.We used approximate results from a minimizer instead of running a sampler on each of the one hundred mocks, due to computational constraints (see Appendix D for a comparison between the two).We show histograms of the best fit parameter values in Figure 7.Our results here reinforce the earlier conclusions.The means of these distributions are consistent with the values expected from the mocks,  = 1.This again indicates that we are able to obtain unbiased measurements of , using the full shapes of the correlations.
Finally, in Section 4 we discuss the implications of our results and give recommendations for how such an analysis should be performed on real data from eBOSS.As we have shown, the most important contaminants are well modelled and do not have a significant impact on  measurements for  min ≥ 25 ℎ −1 Mpc.Therefore, we focus on discussing the effects that are not modelled in the synthetic data, and what tests can be done on the real data in order to gauge their impact.Most of these are quite well understood, have already been modelled and tested in past Ly BAO analyses, and were not found to have a significant impact.However, these tests should also be replicated in the context of a full-shape analysis.
An eBOSS full-shape analysis, as studied here, could improve AP constraints by as much as a factor of two compared to a BAO-only analysis.We have performed this measurement in a follow-up work presented in Cuceu et al. (2023).While our focus here has been on the information extracted through the AP effect, our validation of full-shape measurements could also prove important to other sources of information present in the broadband of Ly forest correlations.These include growth of structure measurements using redshift space distortions (Cuceu et al. 2021), and fluctuations in the ionising UV background (e.g.Long & Hirata 2023).

DATA AVAILABILITY
Data supporting this research is available on reasonable request from the corresponding author.with the expected noise) can produce large changes in the shape of the posterior (e.g.Cuceu et al. 2020).Furthermore, when changing the fiducial Ω  , the smallest scale that we fit will be different for the same  min .For example, the smaller value of Ω  leads to smaller scales being assigned larger comoving separations.This is why the size of the constraints increases as Ω  increases. 20Nevertheless, all three results are consistent with each other even for these extreme values of Ω  .

APPENDIX C: OTHER SCALE CUTS
In Section 3.2 we have focused on the minimum scale fitted,  min .However, we also have to choose the largest scale we fit for,  max .Furthermore, as our measurement is two-dimensional, we could also impose different scale cuts along or across the line-of-sight.We tested different versions of this analysis, which we summarize here.
We show the impact of changing the largest separation included in our fit,  max , in Figure C1.We perform the same type of comparison as we did for  min in Figure 5.The scales we tested range from 140 ℎ −1 Mpc to 200 ℎ −1 Mpc.We find that measurements of   are unbiased for all the  max values we have tried.Furthermore, the improvement in the constraint as we include more data is very weak.Therefore, following eBOSS, we decided to choose  max = 180 ℎ −1 Mpc in our standard analysis.
As most of the contaminants of the Ly forest correlation functions have their largest impact along the line-of-sight, we also decided to Histograms of the 68% constraints on  p ,  p and  s when fitting the correlations measured from one hundred uncontaminated eBOSS mocks.We compare results from a fitter which assume a Gaussian posterior, with results from a sampler where we compute the full posterior distributions.The vertical lines represent the standard deviation of the best fitting results for each parameter over the one hundred mocks.We show that when using the sampler, the variance of the population best fits is consistent with the individual constraints for all three parameters.

APPENDIX E: PARAMETER CORRELATIONS
The advantage of BAO analyses is that they measure a well-defined feature in the correlation function.BAO scale parameters rarely have any significant correlations with the nuisance parameters, which means they are much more robust to different analysis choices when it comes to modelling contaminants.On the other hand, full-shape analyses rely on the ability to accurately model correlation functions over all the scales of interest.We have already shown in Figure 6 what biases are introduced when not modelling certain contaminants.In this appendix, we go into more detail on the sensitivity of  s measurements to other nuisance parameters.
In Figure E1 we show the posterior distributions of the main parameters that are correlated with  s .The two posteriors are from the fit to the two Ly auto-correlations, and from the joint fit of all four correlations.The cross-correlations between Ly and quasars cannot be included in such a comparison because the full system of parameters is degenerate. 24The most important correlations are between  s and the Ly bias and RSD parameters,  Ly and  Ly .These two are highly correlated between themselves and with the HCD bias,  HCD .The auto-correlation provides most of the information on these parameters.On the other hand, the information on  QSO comes from the Ly-QSO cross-correlation.However, a measurement of  Ly (from the auto) is necessary to constrain it due to their degeneracy in the model of the cross-correlation.While we treat  QSO as a nuisance parameter in this work, future studies could also attempt to extract a measurement of the growth rate of structure from it.Posterior distributions obtained from fitting the mean of correlation functions computed from one hundred eBOSS mocks.We show results for the joint analysis of the two Ly auto-correlations, and for the joint analysis of all four correlations.We only plot the parameters that have correlations with  s , in order to show its sensitivity to the values of these parameters.
accounts for the binning of the correlation function, with  | | and  ⊥ the radial and transverse bin widths, respectively.

Figure 1 .
Figure 1.Stacked correlation functions computed from one hundred eBOSS DR16 mocks, compressed into four  =  || / bins.The top panels show the auto-correlation of Ly forest transmitted flux within the Ly region on the left and between the Ly and Ly regions on the right.The bottom plots show the cross-correlation between Ly forest transmitted flux in the Ly region (left) and Ly region (right) and the quasar distribution.The lines represent the joint best fit model of all four correlations.The error-bars are plotted for all points, but except at large separations, they are too small to see.The shaded bands represent the uncertainties from one eBOSS realization.
Following Rogers et al. (2018) and de Sainte Agathe et al. (2019), we model  HCD ( | | ) as an exponential:  HCD = exp(− HCD  | | ).The parameter  HCD can be interpreted as the typical length scale of unmasked HCDs (de Sainte Agathe et al. 2019), and is fixed to  HCD = 10 ℎ −1 Mpc.We test the impact of this choice in Section 3.3.
where  | | and  ⊥ are comoving separations along and across the line-of-sight, respectively.These separations are rescaled using the scale parameters ( p | | ,  p ⊥ ) for the peak, and ( s | | ,  s ⊥ ) for the smooth component.The two peak scale parameters are equivalent to the ( | | ,  ⊥ ) parameters used in BAO analyses.

Figure 2 .
Figure 2. Models of the Ly forest auto-correlation function in a shell in isotropic separation , shown as a function of the cosine of the line-ofsight angle, .The top panel shows the effect of changing the value of , the parameter measuring the Alcock-Paczyński effect, while the bottom plot shows the effect of changing the value of the Ly RSD parameter,  Ly .This shows that the two effects (AP and RSD) lead to different changes in the correlation as a function of , which explains why they are not degenerate.

Figure 3 .
Figure 3. Contour plots of the constraints on scale parameters from the stack of 100 eBOSS mock correlations.The result from the two auto-correlations of Ly flux, Ly(Ly)×Ly(Ly) and Ly(Ly)×Ly(Ly), is denoted as "Auto", and the one from the two cross-correlations with quasars, Ly(Ly)×QSO and Ly(Ly)×QSO, is denoted as "Cross".All scale parameters are expected to be unity, given that the same cosmology was used to construct the mocks and to analyse the data.We also show the expected eBOSS constraint for comparison, but note that the 1 2D contours are all larger than the bounds of the plots here.We find that we are able to recover an unbiased measurement of the AP parameter from the full shape of the Ly forest correlations.

Figure 5 .
Figure5.Comparison of the expected statistical constraint on  from eBOSS (  ) and the systematic bias (Δ) as a function of minimum separation fitted,  min .The  parameter here is fitted from the full shapes of all four Ly correlation functions (  ).The systematic bias is obtained from the difference between the best fit value on the stack of 100 mocks and the true value in the mocks.The statistical constraint is obtained by rescaling the constraint from the stack of mocks to one eBOSS realization.Based on this figure, we choose  min = 25 ℎ −1 Mpc in our baseline analysis.

Figure 7 .
Figure7.Histograms of the best fit values of scale parameters obtained by minimizing the  2 over all parameters in each of the 100 eBOSS mocks.The red dashed lines represent the mean of the populations, while the dashedblack lines show the expected values of the scale parameters (unity).The top plots show fits of  p (left) and  s (right).These show that we are able to obtain unbiased measurement of  from the full shape of the Ly forest correlations, which is confirmed by the bottom right plot, showing results for   .The bottom left plot shows results from fits of the BAO isotropic scale parameter  p .The mean of the population in this case is slightly higher (by ∼ 0.25) than the expected value.

Figure B1 .
Figure B1.Posterior distributions of the combination of distances we measure for different Ω  values used in the fiducial cosmology.The combination   /  corresponds to   , while     / 2  corresponds to  p .This shows that we measure consistent distances even when using a very different fiducial cosmology.
Figure C3.Histograms of the 68% constraints on  p ,  p and  s when fitting the correlations measured from one hundred uncontaminated eBOSS mocks.We compare results from a fitter which assume a Gaussian posterior, with results from a sampler where we compute the full posterior distributions.The vertical lines represent the standard deviation of the best fitting results for each parameter over the one hundred mocks.We show that when using the sampler, the variance of the population best fits is consistent with the individual constraints for all three parameters.
Figure E1.Posterior distributions obtained from fitting the mean of correlation functions computed from one hundred eBOSS mocks.We show results for the joint analysis of the two Ly auto-correlations, and for the joint analysis of all four correlations.We only plot the parameters that have correlations with  s , in order to show its sensitivity to the values of these parameters.
; du Mas des Bourboux et al. 2017; de Sainte Agathe et al. 2019; Blomqvist et al. 2019; du Mas des Bourboux et al. 2020).However, as mentioned above, it has only been used to measure BAO.

Table 1 .
The types and names of the four correlation functions we use in this work, along with the tracers used in each of them.Here, we use "flux" to refer to the transmitted flux fraction.The Ly and Ly regions are defined by the rest-frame intervals  RF ∈ [104, 120] nm (situated between the Ly and Ly peaks) and  RF ∈ [92, 102] nm (situated blueward of the Ly peak), respectively.

Table 2 .
The two types of analysis we perform in this work, along with the scale parameters we measure in each of them and their description.

Table 3 .
List of the free parameters in our model, along with their description and priors. (min,max) represents a flat prior within that interval, while N ( ,  2 ) represents a Gaussian prior.Marginalized posterior distributions of the anisotropic scale parameter , from a joint analysis of all four Ly stacked correlation functions.We compare  constraints independently measured from the peak and smooth components (denoted as  p and  s ), with constraints measured from the full-shape of the correlations (denoted   ).All measurements are consistent with the cosmology used to create and analyse the mocks, given by  = 1.