# The Alcock-Paczyński effect from Lyman-$\alpha$ forest correlations:
Analysis validation with synthetic data

Andrei Cuceu,1,2,3,4 Andreu Font-Ribera,5,4 Paul Martini,2,1 Benjamin
Joachimi,4 Seshadri Nadathur,6 James Rich,7 Alma X. González-Morales,8 Hélion
du Mas des Bourboux9 and James Farr4
1Center for Cosmology and Astro-Particle Physics, The Ohio State University,
Columbus, Ohio 43210, USA
2Department of Astronomy, The Ohio State University, Columbus, Ohio 43210, USA
3Department of Physics, The Ohio State University, Columbus, Ohio 43210, USA
4Department of Physics and Astronomy, University College London, Gower Street,
London WC1E 6BT, United Kingdom
5Institut de Física d’Altes Energies, The Barcelona Institute of Science and
Technology, Campus UAB, 08193 Bellaterra (Barcelona), Spain
6Institute of Cosmology and Gravitation, University of Portsmouth, Burnaby
Road, Portsmouth, PO1 3FX, United Kingdom
7IRFU, CEA, Université Paris-Saclay, F91191 Gif-sur-Yvette, France
8Departamento de Física, DCI, Campus León, Universidad de Guanajuato, 37150,
León, Guanajuato, México
9Department of Physics and Astronomy, University of Utah, 115 S 1400 E, Salt
Lake City, UT 84112, USA E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The three-dimensional distribution of the Ly$\alpha$ 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$\alpha$ 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$\alpha$ 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$\alpha$ correlations performs well on current data sets. We
show that we are able to obtain an unbiased full-shape measurement of
$D_{M}/D_{H}(z_{\mathrm{eff}})$, where $D_{M}$ is the transverse comoving
distance, $D_{H}$ is the Hubble distance, and $z_{\mathrm{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 $r_{\mathrm{min}}=25$
$h^{-1}\text{Mpc}$. We also study and discuss the impact of the main
contaminants affecting Ly$\alpha$ forest correlations, and give
recommendations on how to perform such analysis with real data. While the
final eBOSS Ly$\alpha$ BAO analysis measured
$D_{M}/D_{H}(z_{\mathrm{eff}}=2.33)$ with $4\%$ statistical precision, a full-
shape fit of the same correlations could provide a $\sim 2\%$ measurement.

###### keywords:

large-scale structure – distance scale – methods: data analysis

††pubyear: 2022††pagerange: The Alcock-Paczyński effect from Lyman-$\alpha$
forest correlations: Analysis validation with synthetic data–D

##  1 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 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-$\alpha$ (Ly$\alpha$) forest, this type of
measurement has not yet been performed (Cuceu et al., 2021).

The Ly$\alpha$ forest is made up of absorption lines blueward of the
Ly$\alpha$ 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$\alpha$ line is progressively redshifted as the photons
travel towards us. This means the distribution of Ly$\alpha$ 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-dimensional111The one
dimension here is along the line-of-sight. 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$\alpha$ 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; 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.

The BOSS and eBOSS Ly$\alpha$ cosmological analyses have progressively
improved our understanding of the different processes that contribute to the
auto-correlation function of Ly$\alpha$ 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$\alpha$ absorption in the Ly$\beta$ region (blueward of the
Ly$\beta$ peak) was also used, first through its correlation with Ly$\alpha$
absorption in the Ly$\alpha$ region (de Sainte Agathe et al., 2019), and then
through its cross-correlation with quasars (Blomqvist et al., 2019). The model
of the Ly$\alpha$ 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$\alpha$ 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
auto- and 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$\alpha$ forest. RSD measurements themselves are also generally
used to measure the growth rate of cosmic structure, but for the Ly$\alpha$
auto-correlation this is degenerate with an unknown RSD bias parameter
(McDonald, 2003; Slosar et al., 2011; Givans & Hirata, 2020; Chen et al.,
2021). Therefore, in this work, we only focus on the measurement of the AP
effect, and choose to marginalize over RSD.222Cuceu et al. (2021) did show
that it is possible to measure the growth rate through a joint analysis of the
Ly$\alpha$ auto-correlation and its cross-correlation with quasars. However,
this would not be a significant measurement with eBOSS.

Our objective in this work is to analyse the full shape of the Ly$\alpha$
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$\alpha$ 3D correlation functions. The 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.

##  2 Framework

We begin with an overview of the framework we use to validate a full-shape
analysis of the Ly$\alpha$ 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$\alpha$ 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.

###  2.1 Synthetic data-sets

In order to test the model of the Ly$\alpha$ 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$\alpha$ eBOSS
DR16 data set. Each mock is based on a Gaussian random field generated with
the CoLoRe333https://github.com/damonge/CoLoRe 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
field444Computed using the gradient of the Newtonian potential. from each
quasar to the centre of the box.

The skewers generated by CoLoRe require significant post-processing in order
to turn them into realistic simulations of the Ly$\alpha$ forest. This is
performed by the LyaCoLoRe555https://github.com/igmhub/LyaCoLoRe package (Farr
et al., 2020). In order to create the large boxes needed to simulate an all-
sky light-cone to $z=3.7$ ($\sim 10$ $h^{-1}\text{Gpc}$ in this case), the
CoLoRe grid is limited by memory and computational constraints to a resolution
of $\sim 2.4$ $h^{-1}\text{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$\alpha$ forest (McDonald & Eisenstein, 2007), often
characterized by the one-dimensional power spectrum, or $P_{\rm 1D}$. In order
to simulate spectra with a realistic $P_{\rm 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 ($\tau$) 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 ($s$) is given by $F(s)=\exp{[-\tau(s)]}$. For more details on
this process, and how it is tuned to produce realistic forests, see Farr et
al. (2020).

Beyond the Ly$\alpha$ 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$\alpha$ 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$\beta$ 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$\alpha$ 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 $\sigma_{z}=400\text{ 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$\alpha$ 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.

###  2.2 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$\alpha$ 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 article666With the exception of those in Appendix A. come from
the same mock correlation functions computed for the eBOSS Ly$\alpha$ BAO
analysis.

The first step in the analysis process is to compute the flux transmission
fluctuation, $\delta_{q}(\lambda)$, of each quasar $q$, based on the observed
flux, $f_{q}(\lambda)$:

$\delta_{q}(\lambda)=\frac{f_{q}(\lambda)}{\bar{F}(\lambda)C_{q}(\lambda)}-1,$
(1)

where $\bar{F}(\lambda)$ is the global mean transmission, and $C_{q}(\lambda)$
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
$\bar{F}(\lambda)C_{q}(\lambda)$. This fit also necessarily includes density
modes of the size of the forest and larger, which biases $\delta_{q}$ 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, $\eta_{ij}^{q}$, such that:

$\sum_{j}\eta_{ij}^{q}\delta_{q}^{m}(\lambda_{j})=\sum_{j}\eta_{ij}^{q}\delta_{q}^{t}(\lambda_{j}),$
(2)

where $\delta_{q}^{m}$ and $\delta_{q}^{t}$ are the measured and true flux
fluctuations, $j$ indexes forest pixels before the projection, and $i$ 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).

Type | Name | Tracer 1 | Tracer 2
---|---|---|---
Auto | Ly$\alpha$(Ly$\alpha$)$\times$Ly$\alpha$(Ly$\alpha$) | Ly$\alpha$ flux in the Ly$\alpha$ region | Ly$\alpha$ flux in the Ly$\alpha$ region
Ly$\alpha$(Ly$\alpha$)$\times$Ly$\alpha$(Ly$\beta$) | Ly$\alpha$ flux in the Ly$\alpha$ region | Ly$\alpha$ flux in the Ly$\beta$ region
Cross | Ly$\alpha$(Ly$\alpha$)$\times$QSO | Ly$\alpha$ flux in the Ly$\alpha$ region | Quasars
Ly$\alpha$(Ly$\beta$)$\times$QSO | Ly$\alpha$ flux in the Ly$\beta$ region | Quasars
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$\alpha$ and Ly$\beta$ regions
are defined by the rest-frame intervals
$\lambda_{\mathrm{RF}}\in[104,120]\;$nm (situated between the Ly$\alpha$ and
Ly$\beta$ peaks) and $\lambda_{\mathrm{RF}}\in[92,102]\;$nm (situated blueward
of the Ly$\beta$ peak), respectively.

Forests that contain HCDs with column densities $\log
N_{\mathrm{H}\textsc{i}}>20.3$777This is close to the detection threshold of
the DLA detector used when analysing real data. du Mas des Bourboux et al.
(2020) also tested the DLA detector on the mock spectra in order to validate
this approach. 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
N_{\mathrm{H}\textsc{i}}<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, $i$ and $j$, at redshifts $z_{i}$ and $z_{j}$, and separated by
an angle $\Delta\theta$, the comoving separations along and across the line-
of-sight are given by (de Sainte Agathe et al., 2019):

$\displaystyle r_{||}$
$\displaystyle=[D_{\mathrm{c}}(z_{i})-D_{\mathrm{c}}(z_{j})]\cos{\frac{\Delta\theta}{2}},$
(3) $\displaystyle r_{\bot}$
$\displaystyle=[D_{\mathrm{M}}(z_{i})+D_{\mathrm{M}}(z_{j})]\sin{\frac{\Delta\theta}{2}},$
(4)

where $D_{\mathrm{M}}$ is the transverse comoving distance,
$D_{\mathrm{c}}(z)=c\int_{0}^{z}dz^{\prime}/H(z^{\prime})$ is the radial
comoving distance, with $H(z)$ as the Hubble parameter and $c$ as the speed of
light. In this work, we will also refer to the $(r,\mu)$ parametrization,
defined as $r^{2}=r_{||}^{2}+r_{\bot}^{2}$ and $\mu=r_{||}/r$. du Mas des
Bourboux et al. (2020) used a fiducial cosmological model888The fiducial
cosmology is based on fits to cosmic microwave background anisotropy data from
the Planck satellite Planck Collaboration et al. (2016); see Table 2 of du Mas
des Bourboux et al. (2020). to compute the two distances, which will result in
an extra anisotropy in the measured correlations if it differs from the true
cosmology. This is the Alcock-Paczyński effect (Alcock & Paczyński, 1979).

The correlation function is first computed independently in each
HEALPix999https://healpix.jpl.nasa.gov/index.shtml pixel. For the eBOSS
footprint, there are about $880$ pixels (nside $=16$), with each covering
$3.7\times 3.7=13.4\;\mathrm{deg}^{2}$. This corresponds to a $250\times
250\;(h^{-1}\mathrm{Mpc})^{2}$ patch at $z_{\mathrm{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. The small correlations
between subsamples are neglected.

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$\alpha$ flux
in front of a quasar, which is assigned negative $r_{||}$, and Ly$\alpha$ flux
behind a quasar, which is assigned positive $r_{||}$. 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$\alpha$ BAO analysis. The first two are the
auto-correlation of the Ly$\alpha$ flux in the Ly$\alpha$ region,
Ly$\alpha$(Ly$\alpha$)$\times$Ly$\alpha$(Ly$\alpha$), and its cross-
correlation with Ly$\alpha$ flux in the Ly$\beta$ region,
Ly$\alpha$(Ly$\alpha$)$\times$Ly$\alpha$(Ly$\beta$). The other two are cross-
correlations between Ly$\alpha$ forest flux and quasars, and include the
cross-correlation of quasars with Ly$\alpha$ flux in the Ly$\alpha$ region:
Ly$\alpha$(Ly$\alpha$)$\times$QSO, and with Ly$\alpha$ flux in the Ly$\beta$
region: Ly$\alpha$(Ly$\beta$)$\times$QSO. The Ly$\alpha$ region is found
between the Ly$\alpha$ and Ly$\beta$ peaks, and is defined by the rest-frame
interval $\lambda_{\mathrm{RF}}\in[104,120]\;$nm. The Ly$\beta$ region is
blueward of the Ly$\beta$ peak, and is defined by the rest-frame interval
$\lambda_{\mathrm{RF}}\in[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 $\mu$
bins each. Besides the BAO peak at $\sim 100$ $h^{-1}\text{Mpc}$, the other
features present 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
$r_{||}=64$ $h^{-1}\text{Mpc}$ and $r_{||}=56$ $h^{-1}\text{Mpc}$,
respectively. In the auto-correlation, the SiIII$(1207)$ peak can also be seen
at $r_{||}=21$ $h^{-1}\text{Mpc}$. Finally, a fourth metal peak due to
SiII$(1260)$ is present at $r_{||}=111$ $h^{-1}\text{Mpc}$, however, it is not
visible due to its proximity to the BAO peak.

Figure 1:  Stacked correlation functions computed from one hundred eBOSS DR16
mocks, compressed into four $\mu=r_{||}/r$ bins. The top panels show the auto-
correlation of Ly$\alpha$ forest transmitted flux within the Ly$\alpha$ region
on the left and between the Ly$\alpha$ and Ly$\beta$ regions on the right. The
bottom plots show the cross-correlation between Ly$\alpha$ forest transmitted
flux in the Ly$\alpha$ region (left) and Ly$\beta$ 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.

###  2.3 Modelling the correlations

Our model for the Ly$\alpha$ forest auto-correlation
(Ly$\alpha\times$Ly$\alpha$) and its cross-correlation with quasars
(Ly$\alpha\times$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,
$P_{\mathrm{fid}}(k)$, which is split into a peak (or wiggles) component,
$P^{\mathrm{p}}_{\mathrm{fid}}(k)$, and a smooth (or no-wiggles) component,
$P^{\mathrm{s}}_{\mathrm{fid}}(k)$, 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^{\mathrm{p}}_{\mathrm{fid}}(k)$ and $P^{\mathrm{s}}_{\mathrm{fid}}(k)$, in
parallel. We use the Vega package to perform this
analysis.101010https://github.com/andreicuceu/vega

The Ly$\alpha$ auto and cross power spectra are given by:

$\displaystyle P_{\mathrm{auto}}(k,\mu_{k},z)$ $\displaystyle=b^{\prime
2}_{\mathrm{Ly}\alpha}(1+\beta^{\prime}_{\mathrm{Ly}\alpha}\mu_{k}^{2})^{2}\;G(k,\mu_{k})F_{\mathrm{sm}}(k,\mu_{k})P_{\mathrm{fid}}(k),$
(5) $\displaystyle P_{\mathrm{cross}}(k,\mu_{k},z)$
$\displaystyle=b^{\prime}_{\mathrm{Ly}\alpha}(1+\beta^{\prime}_{\mathrm{Ly}\alpha}\mu_{k}^{2})\;G(k,\mu_{k})F_{\mathrm{sm}}(k,\mu_{k})$
$\displaystyle\times
b_{\mathrm{QSO}}(1+\beta_{\mathrm{QSO}}\mu_{k}^{2})\;F_{\mathrm{NL}}(k,\mu_{k})P_{\mathrm{fid}}(k),$
(6)

where $b^{\prime}_{\mathrm{Ly}\alpha}$ and
$\beta^{\prime}_{\mathrm{Ly}\alpha}$ are the effective Ly$\alpha$ bias and RSD
parameter, and $b_{\mathrm{QSO}}$ and $\beta_{\mathrm{QSO}}$ are the quasar
linear bias and RSD parameter.
$G(k,\mu_{k})=\mathrm{sinc}(k_{||}R_{||}/2)\;\mathrm{sinc}(k_{\bot}R_{\bot}/2)$
accounts for the binning of the correlation function, with $R_{||}$ and
$R_{\bot}$ the radial and transverse bin widths, respectively. The
$F_{\mathrm{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:

$\displaystyle F^{2}_{\mathrm{NL,Lorentz}}$
$\displaystyle=[1+(k_{||}\sigma_{v})^{2}]^{-1},$ (7) $\displaystyle
F^{2}_{\mathrm{NL,Gauss}}$
$\displaystyle=\exp\bigg{[}-\frac{1}{2}(k_{||}\sigma_{v})^{2}\bigg{]},$ (8)

where $\sigma_{v}$ 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., 2017, 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,
$F_{\mathrm{sm}}(k,\mu_{k})$, to account for the low-resolution of the CoLoRe
grid, following Farr et al. (2020). This smoothing model has two free
parameters, $\sigma_{||}$ and $\sigma_{\bot}$, which describe the smoothing
along and across the line-of-sight, respectively. They are expected to have
values of $\sim 2.4$ $h^{-1}\text{Mpc}$, corresponding to the resolution of
the CoLoRe grid.

The Ly$\alpha$ 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}\mathrm{cm}^{-2}$, which means they include
both Lyman limit systems and damped Ly$\alpha$ (DLA) systems. As described in
Section 2.2, large HCDs are detected and masked. However, those with widths
smaller than $\sim 14$ $h^{-1}\text{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 $k_{||}$ dependent term in the effective bias and RSD
parameters. This is given by:

$\displaystyle b^{\prime}_{\mathrm{Ly}\alpha}$
$\displaystyle=b_{\mathrm{Ly}\alpha}+b_{\mathrm{HCD}}F_{\mathrm{HCD}}(k_{||}),$
(9) $\displaystyle
b^{\prime}_{\mathrm{Ly}\alpha}\beta^{\prime}_{\mathrm{Ly}\alpha}$
$\displaystyle=b_{\mathrm{Ly}\alpha}\beta_{\mathrm{Ly}\alpha}+b_{\mathrm{HCD}}\beta_{\mathrm{HCD}}F_{\mathrm{HCD}}(k_{||}),$
(10)

where $b_{\mathrm{Ly}\alpha}$ and $\beta_{\mathrm{Ly}\alpha}$ are the linear
bias and RSD parameters associated with the IGM and $b_{\mathrm{HCD}}$ and
$\beta_{\mathrm{Ly}\alpha}$ are the linear bias and RSD parameters associated
with HCDs. Following Rogers et al. (2018) and de Sainte Agathe et al. (2019),
we model $F_{\mathrm{HCD}}(k_{||})$ as an exponential:
$F_{\mathrm{HCD}}=\exp(-L_{\mathrm{HCD}}k_{||})$. The parameter
$L_{\mathrm{HCD}}$ can be interpreted as the typical length scale of unmasked
HCDs (de Sainte Agathe et al., 2019), and is fixed to $L_{\mathrm{HCD}}=10$
$h^{-1}\text{Mpc}$. We test the impact of this choice in Section 3.3.

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 transform111111Performed using the FFTLog
algorithm (Hamilton, 2000). We use the version implemented in the mcfit
package (https://github.com/eelregit/mcfit). to turn each power spectrum
multipole into the corresponding correlation function multipole, and finally
the reconstruction of the two-dimensional correlation function from the
individual multipoles. Following past Ly$\alpha$ BAO analyses, we use
multipole values up to $\ell=6$. The main components of the theoretical model
are the Ly$\alpha$ transmitted flux correlations:
$\xi_{\mathrm{Ly}\alpha\times\mathrm{Ly}\alpha}$ for the auto-correlation, and
$\xi_{\mathrm{Ly}\alpha\times\mathrm{QSO}}$ for the cross-correlation with
quasars. They are computed from the two power spectra in Equations 5 and 6,
respectively.

Analysis type | Parameters | Description | Cosmological Information
---|---|---|---
Split AP measurement | $\phi_{\mathrm{p}}$ | Peak anisotropy | AP from BAO
$\phi_{\mathrm{s}}$ | Broadband anisotropy | AP from broadband
$\alpha_{\mathrm{p}}$ | Peak isotropic scale | BAO scale
$\alpha_{\mathrm{s}}$ | Broadband isotropic scale | Marginalized out
Full AP measurement | $\phi_{f}$ | Full-shape anisotropy | AP from full correlation
$\alpha_{\mathrm{p}}$ | Peak isotropic scale | BAO scale
$\alpha_{\mathrm{s}}$ | Broadband isotropic scale | Marginalized out
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.

We also model the contamination due to metal absorption, using the same models
as for the Ly$\alpha$ auto and cross-correlation. We compute models for the
correlations between each metal line and Ly$\alpha$
($\xi_{\mathrm{Ly}\alpha\times m}$), and among all metal pairs
($\xi_{m_{1}\times m_{2}}$) for the auto-correlation, and between each metal
line and quasars ($\xi_{\mathrm{QSO}\times m}$) for the cross-correlation.
Each metal line has its own bias and RSD parameter ($b_{m},\beta_{m}$). The
same Ly$\alpha$ and QSO $(b,\beta)$ parameters above are also used for the
cross-correlations between metals and Ly$\alpha$ and quasars. In this case, we
neglect the HCD effects for the Ly$\alpha$ parameters. Following past
Ly$\alpha$ BAO analyses, we fix the metal $\beta$ parameters to $0.5$
(Bautista et al., 2017; du Mas des Bourboux et al., 2020).

Because we assume all forest pixels correspond to Ly$\alpha$ absorption,
correlations involving metals are assigned to the wrong bin in comoving
coordinates due to one (for Ly$\alpha\times m$ and QSO$\times m$) or both (for
$m\times m$) tracers being assigned the wrong redshift. This is modelled
through a metal matrix which transforms the model correlations computed using
the correct separations, $\tilde{r}_{||}$ and $\tilde{r}_{\bot}$, to the
measured model correlations:

$\xi^{\mathrm{metal}}_{A}=\sum_{B}M_{AB}\xi^{\mathrm{metal}}(\tilde{r}_{||}(B),\tilde{r}_{\bot}(B)),$
(11)

where $A$ and $B$ are comoving separation bins, and the metal matrix is given
by:

$M_{AB}=\frac{1}{W_{A}}\sum_{(i,j)\in A,(i,j)\in B}w_{i}w_{j},$ (12)

where $(i,j)\in A$ refers to separations computed using the assigned redshifts
(assuming only Ly$\alpha$ absorption), and $(i,j)\in B$ refers to separations
computed using the true redshifts of each tracer. The weights, $w_{i}$, are
given by the inverse variance of $\delta_{q}$ (see du Mas des Bourboux et al.
2020), and $W_{A}=\sum_{(i,j)\in A}w_{i}w_{j}$. Note that we drop the quasar
specification from our notation, as the sums are understood to be over all
pixels in each comoving separation bin.

The full model correlations are then given by:

$\displaystyle\xi^{t}_{\mathrm{auto}}$
$\displaystyle=\xi_{\mathrm{Ly}\alpha\times\mathrm{Ly}\alpha}+\sum_{m}\xi_{\mathrm{Ly}\alpha\times
m}+\sum_{m_{1},m_{2}}\xi_{m_{1}\times m_{2}},$ (13)
$\displaystyle\xi^{t}_{\mathrm{cross}}$
$\displaystyle=\xi_{\mathrm{Ly}\alpha\times\mathrm{QSO}}+\sum_{m}\xi_{\mathrm{QSO}\times
m},$ (14)

where the sums are performed over the four metal lines introduced in Section
2.2.

In order to model the distortion due to continuum fitting, we need to apply
the projection $\eta_{ij}^{q}$ to the modelled flux fluctuation
$\delta^{t}_{q}$ (Equation 2). This is done by using $\eta_{ij}^{q}$ to
construct a distortion matrix that is applied to the theoretical model
correlations, $\xi_{\mathrm{auto}}^{t}$ and $\xi_{\mathrm{cross}}^{t}$ to give
the projected correlation function:

$\hat{\xi}_{A}=\sum_{B}D_{AB}\xi_{B}^{t},$ (15)

where $A$ is a bin of the data and $B$ is a bin of the model. The distortion
matrices for the auto and cross-correlations are given by:

$\displaystyle D_{AB}^{\mathrm{auto}}$ $\displaystyle=W_{A}^{-1}\sum_{ij\in
A}w_{i}w_{j}\sum_{i^{\prime}j^{\prime}\in
B}\eta_{ii^{\prime}}\eta_{jj^{\prime}},$ (16) $\displaystyle
D_{AB}^{\mathrm{cross}}$ $\displaystyle=W_{A}^{-1}\sum_{ij\in
A}w_{i}w_{j}\sum_{i^{\prime}j\in B}\eta_{ii^{\prime}}.$ (17)

As described at the start of this section, we perform all these steps using
both isotropic power spectrum components, $P^{\mathrm{p}}_{\mathrm{fid}}(k)$
and $P^{\mathrm{s}}_{\mathrm{fid}}(k)$, in parallel. Therefore, we compute the
projected correlation function for both the peak component,
$\hat{\xi}_{\mathrm{p}}$, and the smooth component, $\hat{\xi}_{\mathrm{s}}$.
The two model correlations are then combined into the final theoretical model,
allowing for their comoving coordinates to vary:

$\xi(r_{||},r_{\bot})=\hat{\xi}_{\mathrm{s}}(q_{||}^{\mathrm{s}}r_{||},q_{\bot}^{\mathrm{s}}r_{\bot})+\hat{\xi}_{\mathrm{p}}(q_{||}^{\mathrm{p}}r_{||},q_{\bot}^{\mathrm{p}}r_{\bot}),$
(18)

where $r_{||}$ and $r_{\bot}$ are comoving separations along and across the
line-of-sight, respectively. These separations are rescaled using the scale
parameters ($q_{||}^{\mathrm{p}}$, $q_{\bot}^{\mathrm{p}}$) for the peak, and
($q_{||}^{\mathrm{s}}$, $q_{\bot}^{\mathrm{s}}$) for the smooth component. The
two peak scale parameters are equivalent to the ($\alpha_{||}$,
$\alpha_{\bot}$) parameters used in BAO analyses.

In the case of the cross-correlation, we also allow for a systematic shift in
the $r_{||}$ coordinate of the cross-correlation through a nuisance parameter,
$\Delta r_{||}$. This means the coordinate transform is given by
$r_{||}\xrightarrow{}q_{||}(r_{||}+\Delta r_{||})$. 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$\alpha$ BAO analyses
(e.g. Blomqvist et al., 2019; du Mas des Bourboux et al., 2020).

In order to isolate the Alcock-Paczyński effect, we transform the scale
parameter system, following Cuceu et al. (2021):

$\phi(z)\equiv\frac{q_{\bot}(z)}{q_{||}(z)}\;\;\text{and}\;\;\alpha(z)\equiv\sqrt{q_{\bot}(z)q_{||}(z)},$
(19)

where $\phi$ is the anisotropic parameter, and $\alpha$ is the isotropic scale
parameter. Cuceu et al. (2021) showed that a measurement of $\phi$ is
equivalent to a measurement of the AP parameter $D_{M}(z)H(z)$. As we have two
sets of ($q_{||}$,$q_{\bot}$) parameters, we also have two sets of ($\phi$,
$\alpha$) 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, $\phi_{\mathrm{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 $\phi$ parameters for the peak and smooth components separately.
In the second one, called "full AP measurement", we fit one $\phi$ 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 $\phi_{\mathrm{s}}$ quantifies the gain in information over a standard
BAO-only analysis.

The isotropic scale parameter from the smooth component,
$\alpha_{\mathrm{s}}$, also contains some cosmological information. 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$\alpha$ bias. Furthermore, Cuceu et al. (2021) found that the constraint on
the isotropic scale from the peak, $\alpha_{\mathrm{p}}$, is much tighter than
the one from the smooth component. This is in contrast to $\phi$, where the
constraint from the broadband is the one that dominates (see Figure 2 of Cuceu
et al. 2021). Therefore, for the purposes of this analysis, we will treat
$\alpha_{\mathrm{s}}$ as a nuisance parameter, and marginalize over it.

Parameters |  Description | Prior
---|---|---
$\phi$, $\alpha$ |  Scale parameters (see Table 2) | $U(0.01,2.0)$
$b_{\mathrm{Ly}\alpha}$ |  Ly$\alpha$ linear bias | $U(-2.0,0.0)$
$b_{\mathrm{QSO}}$ |  QSO linear bias | $U(0.0,10.0)$
$\beta_{\mathrm{Ly}\alpha}$, $\beta_{\mathrm{QSO}}$ |  RSD parameters of Ly$\alpha$ and QSOs | $U(0.0,5.0)$
$b_{\mathrm{HCD}}$ |  HCD linear bias | $U(-0.2,0.0)$
$\beta_{\mathrm{HCD}}$ |  RSD parameter for HCDs | $\mathcal{N}(0.5,0.2^{2})$
$\sigma_{v}$ [ $h^{-1}\text{Mpc}$] |  Smoothing for redshift errors and QSO non-linear velocities | $U(0.0,15.0)$
$\Delta r_{||}$ [ $h^{-1}\text{Mpc}$] |  Shift due to QSO redshift errors | $U(-3.0,3.0)$
$b_{\mathrm{SiII}(1190)}$ |  Linear bias of metal absorber | $U(-0.02,0.0)$
$b_{\mathrm{SiII}(1193)}$ |  Linear bias of metal absorber | $U(-0.02,0.0)$
$b_{\mathrm{SiIII}(1207)}$ |  Linear bias of metal absorber | $U(-0.02,0.0)$
$b_{\mathrm{SiII}(1260)}$ |  Linear bias of metal absorber | $U(-0.02,0.0)$
$\sigma_{||}$ [ $h^{-1}\text{Mpc}$] |  Smoothing along the line-of-sight | $U(0.0,10.0)$
$\sigma_{\bot}$ [ $h^{-1}\text{Mpc}$] |  Smoothing across the line-of-sight | $U(0.0,10.0)$
Table 3: List of the free parameters in our model, along with their
description and priors. $U(\text{min,max})$ represents a flat prior within
that interval, while $\mathcal{N}(\mu,\sigma^{2})$ represents a Gaussian
prior.

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 $r_{\mathrm{min}}=25$
$h^{-1}\text{Mpc}$ and $r_{\mathrm{max}}=180$ $h^{-1}\text{Mpc}$ following the
analysis in Section 3.2.121212Note that we still fit for the bias parameter of
the SiIII$(1207)$ line, even 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. This
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 $\sim 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.

##  3 Results

A large part of the validation of our analysis will focus on studying full-
shape fits using a stack of the Ly$\alpha$ 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$\alpha$ 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.

###  3.1 Analysis of stacked correlations

Figure 2:  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$\alpha$ flux, Ly$\alpha$(Ly$\alpha$)$\times$Ly$\alpha$(Ly$\alpha$) and
Ly$\alpha$(Ly$\alpha$)$\times$Ly$\alpha$(Ly$\beta$), is denoted as "Auto", and
the one from the two cross-correlations with quasars,
Ly$\alpha$(Ly$\alpha$)$\times$QSO and Ly$\alpha$(Ly$\beta$)$\times$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\sigma$ 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$\alpha$ forest correlations.

As we have four different correlations, we categorise them into the two auto-
correlations of Ly$\alpha$ flux,
Ly$\alpha$(Ly$\alpha$)$\times$Ly$\alpha$(Ly$\alpha$) and
Ly$\alpha$(Ly$\alpha$)$\times$Ly$\alpha$(Ly$\beta$), and the two cross-
correlations with quasars, Ly$\alpha$(Ly$\alpha$)$\times$QSO and
Ly$\alpha$(Ly$\beta$)$\times$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$\alpha$ flux in the Ly$\beta$ region the same as the
Ly$\alpha$ flux in the Ly$\alpha$ 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$\alpha$ bias
and RSD parameters. Therefore, we fix the bias and RSD parameter of quasars
($b_{\mathrm{QSO}},\beta_{\mathrm{QSO}}$) to a fiducial value when fitting the
cross-correlation alone.131313The fiducial value is based on the joint
constraint of all four correlations. Note that due to correlations between the
AP effect and RSD, this may lead the uncertainty in $\phi$ to be
underestimated. However, we will not use these $\phi$ 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
PolyChord141414https://github.com/PolyChord/PolyChordLite 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 2, we show the results for $\phi_{s}$, $\phi_{p}$, and $\alpha_{p}$
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 $r>r_{\mathrm{min}}=25$
$h^{-1}\text{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.151515See
Appendix A for tests with different fiducial cosmologies. Therefore, we expect
to recover unity for all ($\phi,\alpha$) parameters. We do recover the
expected $\phi$ values from both the auto and cross-correlation, using both
the peak and smooth components independently. The fact that we are able to
constrain $\phi_{\mathrm{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
$\phi_{\mathrm{s}}$ and important nuisance parameters see Appendix B. Finally,
Figure 2 shows that we are indeed able to perform an unbiased analysis of the
full-shape of Ly$\alpha$ forest correlations for mocks with known
contaminants.

For the isotropic scale of the peak, $\alpha_{\mathrm{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 $\sim 0.22\%$
larger than unity. While this appears as a $\sim 2\sigma$ systematic bias in
the analysis on the stack of 100 eBOSS mocks, it is a small fraction ($\sim
0.2\sigma$) of the expected eBOSS statistical uncertainty. We will return to
discuss this deviation in Sections 3.4 and 4.

While for the BAO peak parameters, $\alpha_{\mathrm{p}}$ and
$\phi_{\mathrm{p}}$ the constraints obtained from the auto and cross-
correlations are very similar, this is not the case for $\phi_{\mathrm{s}}$,
where the measurement from the auto-correlation is $\sim 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$\alpha$ flux in
the Ly$\beta$ forest shrinks the $\phi_{\mathrm{s}}$ constraint by $16\%$ for
the auto-correlation and by $10\%$ for the cross-correlation. In Figure 2 we
only focus on the parameters of interest for cosmology. However, unlike BAO
parameters, $\phi_{\mathrm{s}}$ does have correlations with some nuisance
parameters, the most important being $b_{\mathrm{Ly}\alpha}$,
$\beta_{\mathrm{Ly}\alpha}$ and $b_{\mathrm{HCD}}$. We show these in Appendix
B.

Figure 3:  Marginalized posterior distributions of the anisotropic scale
parameter $\phi$, from a joint analysis of all four Ly$\alpha$ stacked
correlation functions. We compare $\phi$ constraints independently measured
from the peak and smooth components (denoted as $\phi_{\mathrm{p}}$ and
$\phi_{\mathrm{s}}$), with constraints measured from the full-shape of the
correlations (denoted $\phi_{f}$). All measurements are consistent with the
cosmology used to create and analyse the mocks, given by $\phi=1$.

In Figure 3, 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 $\phi$ measurements from the peak and smooth
components ($\phi_{\mathrm{p}}$ and $\phi_{\mathrm{s}}$) and the measurement
from the full shape ($\phi_{f}$). 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 $\sim 32\%$ tighter than the
constraint from the peak component. Note that we do not expect the
$\phi_{\mathrm{p}}$ and $\phi_{\mathrm{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 2.

The measurement of $\phi_{f}$ is remarkably consistent (to within $0.1\sigma$)
with the expected $\phi=1$. Overall, measuring $\phi$ from the full-shape of
the correlations leads to a $\sim 41\%$ improvement in the $1\sigma$
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.

###  3.2 Statistical and systematic errors

The results we have shown so far are already quite promising, with $\phi_{s}$
and $\phi_{f}$ measurements that are consistent with the true cosmology of the
mocks ($\phi=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,
$r_{\mathrm{min}}$.

In order to quantify the systematic bias of our $\phi$ measurements, we
compute the difference between the best fit value of $\phi$ from the stacked
correlations and the true value in the mocks, i.e.,
$\Delta\phi=\phi_{\mathrm{best}}-1$. While we do have the size of the
statistical constraints from the stack of the correlations, we want to compare
$\Delta\phi$ 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, $\sigma_{\phi}$.

Figure 4:  Comparison of the expected statistical constraint on $\phi$ from
eBOSS ($\sigma_{\phi}$) and the systematic bias ($\Delta\phi$) as a function
of minimum separation fitted, $r_{\mathrm{min}}$. The $\phi$ parameter here is
fitted from the full shapes of all four Ly$\alpha$ correlation functions
($\phi_{f}$). 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
$r_{\mathrm{min}}=25$ $h^{-1}\text{Mpc}$ in our baseline analysis.

We compare the statistical constraint and the systematic error on $\phi$ for
different values of $r_{\mathrm{min}}$ in Figure 4. We use the full AP setup
and show results for $\phi_{f}$ from the joint analysis of all four
correlation functions. The error-bars on the systematic errors are given by
the $1\sigma$ statistical constraints of each measurement, which represents
the constraining power of the stack of the one hundred mocks.

As expected, $\sigma_{\phi}$ grows significantly as we cut more data (larger
$r_{\mathrm{min}}$), with values ranging from $1.6\%$ to $2.7\%$. The
systematic bias starts high for $r_{\mathrm{min}}=10$ $h^{-1}\text{Mpc}$
($0.7\%$), but becomes consistent with zero by $r_{\mathrm{min}}=25$
$h^{-1}\text{Mpc}$. Note that we initially performed this test in
$r_{\mathrm{min}}$ steps of $10$ $h^{-1}\text{Mpc}$, but chose to compute the
$r_{\mathrm{min}}=25$ $h^{-1}\text{Mpc}$ point after seeing the initial
trends. We also performed the same analysis on $\phi_{\mathrm{s}}$, using the
split AP approach, and found similar results (systematic bias consistent with
zero for $r_{\mathrm{min}}\geq 25$ $h^{-1}\text{Mpc}$).

These results confirm that we are able to obtain unbiased measurements of
$\phi$ using the full shapes of the Ly$\alpha$ forest correlation functions in
analyses using $r_{\mathrm{min}}\geq 25$ $h^{-1}\text{Mpc}$. Therefore,
throughout the rest of this work, we will be using $r_{\mathrm{min}}=25$
$h^{-1}\text{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
$r_{\mathrm{max}}=180$ $h^{-1}\text{Mpc}$, and not impose any anisotropic
cuts.

###  3.3 Modelling of contaminants

The results we have seen so far indicate that our current modelling of the
relevant contaminants in the Ly$\alpha$ 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 $k_{\parallel})$, 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
$\phi_{\mathrm{s}}$ and $\phi_{f}$ (Figure 3), 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 $\sim 14$ $h^{-1}\text{Mpc}$
remain (de Sainte Agathe et al., 2019; du Mas des Bourboux et al., 2020). The
current model of this effect (Equations 9 and 10) follows the description
introduced by Font-Ribera & Miralda-Escudé (2012), combined with an
exponential shape for the $F_{\mathrm{HCD}}(k_{||})$ 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
N_{\mathrm{H}\textsc{i}}<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$\alpha$ and HCD auto-correlations and their cross-
correlation). However, Font-Ribera & Miralda-Escudé (2012) found that the
three-point function,
$\langle\delta_{\mathrm{Ly}\alpha}\delta_{\mathrm{HCD}}\delta_{\mathrm{Ly}\alpha}\rangle$,
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 $r_{\mathrm{min}}=25$ $h^{-1}\text{Mpc}$. However, they could
be contributing to the bias we observe for $r_{\mathrm{min}}<25$
$h^{-1}\text{Mpc}$. We have also tested the impact of marginalizing over the
$L_{\mathrm{HCD}}$ parameter (instead of fixing it to $10$
$h^{-1}\text{Mpc}$), and found that it has no impact on $\phi_{\mathrm{s}}$
measurements.

Contamination due to metal absorption results in clear extra features in the
Ly$\alpha$ correlation functions along the line-of-sight. These features
correspond to correlations at vanishing separation between metal absorption
and either Ly$\alpha$ absorption (auto) or quasars (cross). These correlations
get assigned to the wrong bins in comoving coordinates because we interpret
everything as Ly$\alpha$ absorption. This is a simple enough process that we
can replicate it in the model through the metal matrix formalism described in
Section 2.3. 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.

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 $r_{||}$ 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, $\Delta r_{||}$, 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
($\alpha_{\mathrm{p}}$, $\phi_{\mathrm{p}}$, $\phi_{\mathrm{s}}$).

Figure 5:  Posterior distributions of the anisotropic scale parameter, $\phi$,
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 $\phi_{\mathrm{s}}$, the most
significant shift appears when removing the model for HCDs.

In Figure 5, we compare the baseline results on $\phi$ 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$\alpha$ 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 $\phi_{\mathrm{s}}$, with the largest resulting
from the removal of the HCD model, leading to a $\sim 1\%$ shift. On the other
hand, for $\phi_{\mathrm{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 $\sim
0.5\sigma$ of the expected statistical constraint for eBOSS.

###  3.4 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 $r_{\mathrm{min}}=25$ $h^{-1}\text{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 analyses. Therefore, we now wish to test if our conclusions
above are also accurate when analysing individual eBOSS mocks.

Figure 6:  Histograms of the best fit values of scale parameters obtained by
minimizing the $\chi^{2}$ over all parameters in each of the 100 eBOSS mocks.
The red dashed lines represent the mean of the populations, while the dashed-
black lines show the expected values of the scale parameters (unity). The top
plots show fits of $\phi_{\mathrm{p}}$ (left) and $\phi_{\mathrm{s}}$ (right).
These show that we are able to obtain unbiased measurement of $\phi$ from the
full shape of the Ly$\alpha$ forest correlations, which is confirmed by the
bottom right plot, showing results for $\phi_{f}$. The bottom left plot shows
results from fits of the BAO isotropic scale parameter $\alpha_{\mathrm{p}}$.
The mean of the population in this case is slightly higher (by $\sim
0.25\sigma$) than the expected value.

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.161616Running the sampler for a joint analysis of
all four Ly$\alpha$ correlations with all the contaminants requires $\sim
2-3\times 10^{3}$ CPU hours, which leads to a rough estimate of $\sim
2-3\times 10^{5}$ thousand CPU hours to analyse all 100 mocks. Therefore, in
this section, we compute the best fit values of the parameters using a
minimizer171717We use the iminuit package (Dembinski et al., 2021) for
minimization.. 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
$\phi_{\mathrm{p}}$ and $\phi_{\mathrm{s}}$, and the full AP analysis with
$\phi_{f}$.

We show histograms of the best fit values of the scale parameters in Figure 6.
The top plots show results for $\phi_{\mathrm{p}}$ and $\phi_{\mathrm{s}}$,
both of which are consistent with the expected value of unity. The means of
the populations are $\mu_{\phi_{\mathrm{p}}}=1.0028$ and
$\mu_{\phi_{\mathrm{s}}}=1.0013$, corresponding to shifts from unity of $\sim
0.06\sigma$ and $\sim 0.05\sigma$, respectively.181818The $\sigma$ values here
are based on the expected eBOSS uncertainty given by the standard deviation of
these histograms. This is also consistent with the scaled uncertainty from the
stack of mocks. These 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 6 shows results for the isotropic BAO scale
($\alpha_{\mathrm{p}}$) and the full-shape anisotropic parameter ($\phi_{f}$).
The means of the populations in this case are
$\mu_{\alpha_{\mathrm{p}}}=1.0024$ and $\mu_{\phi_{f}}=1.0024$, corresponding
with shifts from unity of $\sim 0.23\sigma$ and $\sim 0.11\sigma$,
respectively. While these shifts are more significant compared to the ones
from $\phi_{\mathrm{p}}$ and $\phi_{\mathrm{s}}$, they are still relatively
small, in line with the results from Section 3.1.

##  4 Discussion

In this work, we have used a set of one hundred eBOSS Ly$\alpha$ forest mocks
to study how well we could perform a full-shape analysis of the eBOSS DR16
Ly$\alpha$ correlation functions. We found that we are able to recover
unbiased measurements of the anisotropic scale parameter, $\phi$, from a joint
analysis of the full shapes of the four Ly$\alpha$ correlations when using a
minimum separation of $r_{\mathrm{min}}=25$ $h^{-1}\text{Mpc}$. The only
significant systematic biases on $\phi$ were found when either we did not
model one of the major contaminants (Section 3.3), or we included scales
smaller than $25$ $h^{-1}\text{Mpc}$ (Section 3.2).

We did find an apparent bias of $0.22\%-0.24\%$ in the isotropic scale of the
peak component, $\alpha_{\mathrm{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 $\sim 0.21\sigma-0.25\sigma$ bias for
eBOSS.191919Note that the eBOSS uncertainty used for these values are based on
estimates from this work, not the real measurement. du 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 $\phi_{\mathrm{s}}$. However, our model allows the separate fit of BAO
($\alpha_{\mathrm{p}}$ and $\phi_{\mathrm{p}}$) and broadband
($\phi_{\mathrm{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$\alpha$ 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$\alpha$ 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.202020We do not use this model here because
the small-scale forest-related non-linearities are not included in the mocks.
This 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$\alpha$ power spectrum and the cross-spectrum of Ly$\alpha$ with massive
halos. For Ly$\alpha$, 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 small-scale power. Therefore, for the cross-correlation,
they recommend using scales larger than $30$ $h^{-1}\text{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$
$h^{-1}\text{Mpc}$, as proposed here, and of $30$ $h^{-1}\text{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$\alpha$ auto and cross-correlation 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 full-shape 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$\alpha$ 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$\alpha$ peak. This type of contamination is mostly given by the cross-
correlation between these absorptions lines and either Ly$\alpha$ 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$\alpha$ 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\sigma$ (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$\alpha$ BAO analyses.

Fluctuations of the ionizing UV background introduce a scale dependence to the
Ly$\alpha$ bias and RSD parameters (Pontzen, 2014; Gontcho A Gontcho et al.,
2014). This scale dependence was modelled in the Ly$\alpha$ 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 $\sim 2\sigma$), 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 $\phi$ 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 cross-correlation. The
other effects discussed here have been studied before, and modelled at some
level in past Ly$\alpha$ 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 $\phi_{\mathrm{s}}$. If they are
found to have an impact, it will motivate the development of more accurate
mocks in the future.

##  5 Summary

The Lyman-$\alpha$ (Ly$\alpha$) forest provides one of the best tracers of
large scale structure (LSS) at high redshifts ($2<z<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 the
potential of performing a full-shape analysis of the Ly$\alpha$ auto-
correlation, and its cross-correlation with quasars, in order to measure the
Alcock-Paczyński (AP) effect. 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 correlations using a similar approach to that used in past
Ly$\alpha$ 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, $\alpha$, and a parameter for the
anisotropy, $\phi$, 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
$\phi$ parameters ($\phi_{\mathrm{p}}$ and $\phi_{\mathrm{s}}$), while in the
second we fit only one $\phi$ parameter ($\phi_{f}$) 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 2 and 3. We find that
the results from the two Ly$\alpha$ auto-correlations are consistent with the
results from the two cross-correlations. Furthermore, we find that we are able
to recover unbiased measurements of $\phi$ 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,
$r_{\mathrm{min}}$. We run the full AP analysis using $r_{\mathrm{min}}$
values from $10$ $h^{-1}\text{Mpc}$ to $50$ $h^{-1}\text{Mpc}$, as shown in
Figure 4. We then compare the systematic bias of $\phi_{f}$ measurements,
given by the difference between the best fit value and the expected value in
the mocks ($\phi=1$), with an estimate of the expected statistical constraint
from eBOSS. We find that we are able to obtain unbiased measurements of
$\phi_{f}$ for $r_{\mathrm{min}}\geq 25$ $h^{-1}\text{Mpc}$. Therefore, we use
$r_{\mathrm{min}}=25$ $h^{-1}\text{Mpc}$ for the other parts of the analysis.

In Section 3.3, we discuss the main contaminants present in the Ly$\alpha$
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 $\phi$ measurements. We
show the relative impact of not modelling each main contaminant in Figure 5.

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 6. Our results here reinforce the earlier conclusions. The
means of these distributions are consistent with the values expected from the
mocks, $\phi=1$. This again indicates that we are able to obtain unbiased
measurements of $\phi$, 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 $\phi$ measurements for
$r_{\mathrm{min}}\geq 25$ $h^{-1}\text{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$\alpha$ BAO analysis, and were not found to have a significant impact.
However, these tests should also be replicated in the context of a full-shape
analysis. We end by noting that an eBOSS full-shape analysis, as studied here,
could improve $\phi$ constraints by as much as a factor of two compared to a
BAO only analysis.

## Data Availability

Data supporting this research is available on reasonable request from the
corresponding author.

## Acknowledgements

AC and PM acknowledge support from the United States Department of Energy,
Office of High Energy Physics under Award Number DE-SC-0011726. AFR
acknowledges support through the program Ramon y Cajal (RYC-2018-025210) of
the Spanish Ministry of Science and Innovation and from the European Union’s
Horizon Europe research and innovation programme (COSMO-LYA, grant agreement
101044612). IFAE is partially funded by the CERCA program of the Generalitat
de Catalunya. BJ acknowledges support by STFC Consolidated Grant ST/V000780/1.
SN acknowledges support from an STFC Ernest Rutherford Fellowship, grant
reference ST/T005009/2. We acknowledge the use of the GetDist (Lewis, 2019)
and numpy (Harris et al., 2020) packages. For the purpose of open access, the
authors have applied a CC BY public copyright licence to any Author Accepted
Manuscript version arising.

## References

  * Abbott et al. (2019) Abbott T. M. C., et al., 2019, MNRAS, 483, 4866
  * Abbott et al. (2022) Abbott T. M. C., et al., 2022, Phys. Rev. D , 105, 043512
  * Alam et al. (2017) Alam S., et al., 2017, MNRAS, 470, 2617
  * Alam et al. (2021) Alam S., et al., 2021, Phys. Rev. D , 103, 083533
  * Alcock & Paczyński (1979) Alcock C., Paczyński B., 1979, Nature, 281, 358
  * Arinyo-i-Prats et al. (2015) Arinyo-i-Prats A., Miralda-Escudé J., Viel M., Cen R., 2015, J. Cosmology Astropart. Phys., 2015, 017
  * Bautista et al. (2017) Bautista J. E., et al., 2017, A&A, 603, A12
  * Bautista et al. (2021) Bautista J. E., et al., 2021, MNRAS, 500, 736
  * Beutler et al. (2011) Beutler F., et al., 2011, MNRAS, 416, 3017
  * Beutler et al. (2012) Beutler F., et al., 2012, MNRAS, 423, 3430
  * Beutler et al. (2017) Beutler F., et al., 2017, MNRAS, 466, 2242
  * Bi & Davidsen (1997) Bi H., Davidsen A. F., 1997, ApJ, 479, 523
  * Blake et al. (2011) Blake C., et al., 2011, MNRAS, 415, 2876
  * Blomqvist et al. (2015) Blomqvist M., et al., 2015, J. Cosmology Astropart. Phys., 2015, 034
  * Blomqvist et al. (2019) Blomqvist M., et al., 2019, A&A, 629, A86
  * Busca et al. (2013) Busca N. G., et al., 2013, A&A, 552, A96
  * Chen et al. (2021) Chen S.-F., Vlah Z., White M., 2021, J. Cosmology Astropart. Phys., 2021, 053
  * Cole et al. (2005) Cole S., et al., 2005, MNRAS, 362, 505
  * Croft et al. (1998) Croft R. A. C., Weinberg D. H., Katz N., Hernquist L., 1998, ApJ, 495, 44
  * Cuceu et al. (2020) Cuceu A., Font-Ribera A., Joachimi B., 2020, J. Cosmology Astropart. Phys., 2020, 035
  * Cuceu et al. (2021) Cuceu A., Font-Ribera A., Joachimi B., Nadathur S., 2021, MNRAS, 506, 5439
  * Delubac et al. (2015) Delubac T., et al., 2015, A&A, 574, A59
  * Dembinski et al. (2021) Dembinski H., et al., 2021, iminuit: Jupyter-friendly Python interface for C++ MINUIT2, Astrophysics Source Code Library, record ascl:2108.024 (ascl:2108.024)
  * Eisenstein et al. (2005) Eisenstein D. J., et al., 2005, ApJ, 633, 560
  * Farr et al. (2020) Farr J., et al., 2020, J. Cosmology Astropart. Phys., 2020, 068
  * Font-Ribera & Miralda-Escudé (2012) Font-Ribera A., Miralda-Escudé J., 2012, J. Cosmology Astropart. Phys., 2012, 028
  * Font-Ribera et al. (2012) Font-Ribera A., et al., 2012, J. Cosmology Astropart. Phys., 2012, 059
  * Font-Ribera et al. (2013) Font-Ribera A., et al., 2013, J. Cosmology Astropart. Phys., 2013, 018
  * Font-Ribera et al. (2014) Font-Ribera A., et al., 2014, Journal of Cosmology and Astro-Particle Physics, 2014, 027
  * Gil-Marín et al. (2020) Gil-Marín H., et al., 2020, MNRAS, 498, 2492
  * Givans & Hirata (2020) Givans J. J., Hirata C. M., 2020, Phys. Rev. D , 102, 023515
  * Givans et al. (2022) Givans J. J., Font-Ribera A., Slosar A., Seeyave L., Pedersen C., Rogers K. K., Blas D., Iršič V., 2022, arXiv e-prints, p. arXiv:2205.00962
  * Gontcho A Gontcho et al. (2014) Gontcho A Gontcho S., Miralda-Escudé J., Busca N. G., 2014, MNRAS, 442, 187
  * Hamilton (2000) Hamilton A. J. S., 2000, MNRAS, 312, 257
  * Handley et al. (2015a) Handley W. J., Hobson M. P., Lasenby A. N., 2015a, MNRAS, 450, L61
  * Handley et al. (2015b) Handley W. J., Hobson M. P., Lasenby A. N., 2015b, MNRAS, 453, 4384
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hou et al. (2021) Hou J., et al., 2021, MNRAS, 500, 1201
  * Iršič et al. (2013) Iršič V., et al., 2013, J. Cosmology Astropart. Phys., 2013, 016
  * Kirkby et al. (2013) Kirkby D., et al., 2013, J. Cosmology Astropart. Phys., 2013, 024
  * Kirkby et al. (2016) Kirkby D., Bailey S., Guy J., Weaver B. A., 2016, Quick simulations of fiber spectrograph response v0.5, doi:10.5281/zenodo.154130
  * Lewis (2019) Lewis A., 2019, arXiv e-prints, p. arXiv:1910.13970
  * Lyke et al. (2020) Lyke B. W., et al., 2020, ApJS, 250, 8
  * Lynds (1971) Lynds R., 1971, The Astrophysical Journal, 164, L73
  * McDonald (2003) McDonald P., 2003, ApJ, 585, 34
  * McDonald & Eisenstein (2007) McDonald P., Eisenstein D. J., 2007, Phys. Rev. D , 76, 063009
  * McDonald et al. (2006) McDonald P., et al., 2006, ApJS, 163, 80
  * Neveux et al. (2020) Neveux R., et al., 2020, MNRAS, 499, 210
  * Noterdaeme et al. (2009) Noterdaeme P., Petitjean P., Ledoux C., Srianand R., 2009, A&A, 505, 1087
  * Noterdaeme et al. (2012) Noterdaeme P., et al., 2012, A&A, 547, L1
  * Palanque-Delabrouille et al. (2015) Palanque-Delabrouille N., et al., 2015, J. Cosmology Astropart. Phys., 2015, 011
  * Palanque-Delabrouille et al. (2016) Palanque-Delabrouille N., et al., 2016, A&A, 587, A41
  * Palanque-Delabrouille et al. (2020) Palanque-Delabrouille N., Yèche C., Schöneberg N., Lesgourgues J., Walther M., Chabanier S., Armengaud E., 2020, J. Cosmology Astropart. Phys., 2020, 038
  * Percival & White (2009) Percival W. J., White M., 2009, MNRAS, 393, 297
  * Pérez-Ràfols et al. (2018) Pérez-Ràfols I., et al., 2018, MNRAS, 473, 3019
  * Pieri et al. (2010) Pieri M. M., Frank S., Weinberg D. H., Mathur S., York D. G., 2010, ApJ, 724, L69
  * Pieri et al. (2014) Pieri M. M., et al., 2014, MNRAS, 441, 1718
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Pontzen (2014) Pontzen A., 2014, Phys. Rev. D , 89, 083010
  * Prochaska et al. (2010) Prochaska J. X., O’Meara J. M., Worseck G., 2010, ApJ, 718, 392
  * Ramírez-Pérez et al. (2022) Ramírez-Pérez C., Sanchez J., Alonso D., Font-Ribera A., 2022, J. Cosmology Astropart. Phys., 2022, 002
  * Rauch (1998) Rauch M., 1998, Annual Review of Astronomy and Astrophysics, 36, 267
  * Reid et al. (2012) Reid B. A., et al., 2012, MNRAS, 426, 2719
  * Rogers et al. (2018) Rogers K. K., Bird S., Peiris H. V., Pontzen A., Font-Ribera A., Leistedt B., 2018, MNRAS, 476, 3716
  * Ross et al. (2015) Ross A. J., Samushia L., Howlett C., Percival W. J., Burden A., Manera M., 2015, MNRAS, 449, 835
  * Samushia et al. (2014) Samushia L., et al., 2014, MNRAS, 439, 3504
  * Satpathy et al. (2017) Satpathy S., et al., 2017, MNRAS, 469, 1369
  * Seljak et al. (2005) Seljak U., et al., 2005, Phys. Rev. D , 71, 103515
  * Slosar et al. (2011) Slosar A., et al., 2011, J. Cosmology Astropart. Phys., 2011, 001
  * Slosar et al. (2013) Slosar A., et al., 2013, J. Cosmology Astropart. Phys., 2013, 026
  * Viel et al. (2010) Viel M., Haehnelt M. G., Springel V., 2010, J. Cosmology Astropart. Phys., 2010, 015
  * Youles et al. (2022) Youles S., et al., 2022, arXiv e-prints, p. arXiv:2205.06648
  * Zafar et al. (2013) Zafar T., Péroux C., Popping A., Milliard B., Deharveng J. M., Frank S., 2013, A&A, 556, A141
  * de Sainte Agathe et al. (2019) de Sainte Agathe V., et al., 2019, A&A, 629, A85
  * du Mas des Bourboux et al. (2017) du Mas des Bourboux H., et al., 2017, A&A, 608, A130
  * du Mas des Bourboux et al. (2020) du Mas des Bourboux H., et al., 2020, ApJ, 901, 153

##  Appendix A Impact of fiducial cosmology

Figure 7: Posterior distributions of the combination of distances we measure
for different $\Omega_{m}$ values used in the fiducial cosmology. The
combination $D_{M}/D_{H}$ corresponds to $\phi_{f}$, while
$D_{M}D_{H}/r_{d}^{2}$ corresponds to $\alpha_{\mathrm{p}}$. This shows that
we measure consistent distances even when using a very different fiducial
cosmology.

Our method requires the use of a fiducial cosmology, both for defining the
comoving coordinate grid on which the correlations are computed, and for the
template power spectrum that is used to construct the theoretical model for
the Ly$\alpha$ correlations. This fiducial cosmology is based on cosmic
microwave background anisotropy data from the Planck satellite (Planck
Collaboration et al., 2016). In this appendix, we test the impact of this
choice of fiducial cosmology.

We recompute the Ly$\alpha$ auto and cross-correlations of one of the mocks
using two different values of the matter density parameter, $\Omega_{m}$. The
two values are $\Omega_{m}=0.27$ and $\Omega_{m}=0.36$. We then perform our
standard analysis and measure $\alpha_{\mathrm{p}}$ and $\phi_{f}$, but we use
different power spectrum templates to keep the fiducial $\Omega_{m}$ value
consistent in each analysis. In order to compare the results, we compute the
relevant combinations of distances as derived parameters. Following Cuceu et
al. (2021), measurements of $\phi$ and $\alpha_{\mathrm{p}}$ correspond to:

$\displaystyle\textbf{AP: }\;\phi(z)$
$\displaystyle=\frac{D_{\mathrm{M}}(z)H(z)}{[D_{\mathrm{M}}(z)H(z)]_{\text{fid}}},$
(20) $\displaystyle\textbf{BAO: }\;\alpha_{\mathrm{p}}(z)$
$\displaystyle=\sqrt{\frac{D_{\mathrm{M}}(z)D_{\mathrm{H}}(z)/r_{\mathrm{d}}^{2}}{[D_{\mathrm{M}}(z)D_{\mathrm{H}}(z)/r_{\mathrm{d}}^{2}]_{\text{fid}}}},$
(21)

where $D_{H}=c/H$ and $r_{d}$ is the size of the sound horizon at the end of
the drag epoch. Based on these equations and the relevant fiducial distances
for each value of $\Omega_{m}$, we transform $\phi_{f}$ into $D_{M}/D_{H}$,
and $\alpha_{\mathrm{p}}$ into $D_{M}D_{H}/r_{d}^{2}$.

Figure 8: 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$\alpha$ auto-correlations, and for the
joint analysis of all four correlations. We only plot the parameters that have
correlations with $\phi_{\mathrm{s}}$, in order to show its sensitivity to the
values of these parameters.

We compare the results for the large and small values of $\Omega_{m}$ with the
one obtained from our standard analysis (using $\Omega_{m}=0.315$) in Figure
7. All three results give consistent measurements, however, they are not
identical. This can be attributed to the fact that pixel pairs are distributed
to different bins in the correlation when using different values of
$\Omega_{m}$, because the comoving distances change. This will have an impact
on the measured parameters, especially for BAO, where small changes in the
correlation measurements (consistent 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 $\Omega_{m}$, the smallest scale that
we fit will be different for the same $r_{\mathrm{min}}$. For example, the
smaller value of $\Omega_{m}$ leads to smaller scales being assigned larger
comoving separations. This is why the size of the constraints increases as
$\Omega_{m}$ increases.212121This also means that our choice of
$r_{\mathrm{min}}=25$ $h^{-1}\text{Mpc}$ is based on the Planck cosmology.
Therefore, the analysis should be re-done, and a new minimum scale chosen, if
we had reason to believe the real cosmology was very different from the one
inferred by Planck. Nevertheless, all three results are consistent with each
other even for these extreme values of $\Omega_{m}$.

##  Appendix B 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 5 what biases are introduced when not modelling certain
contaminants. In this appendix, we go into more detail on the sensitivity of
$\phi_{\mathrm{s}}$ measurements to other nuisance parameters.

In Figure 8 we show the posterior distributions of the main parameters that
are correlated with $\phi_{\mathrm{s}}$. The two posteriors are from the fit
to the two Ly$\alpha$ auto-correlations, and from the joint fit of all four
correlations. The cross-correlations between Ly$\alpha$ and quasars cannot be
included in such a comparison because the full system of parameters is
degenerate.222222See Section 3.1 above and Section 4 of Cuceu et al. (2021)
The most important correlations are between $\phi_{\mathrm{s}}$ and the
Ly$\alpha$ bias and RSD parameters, $b_{\mathrm{Ly}\alpha}$ and
$\beta_{\mathrm{Ly}\alpha}$. These two are highly correlated between
themselves and with the HCD bias, $b_{\mathrm{HCD}}$. The Ly$\alpha$ auto-
correlation provides most of the information on these parameters. On the other
hand, the information on $\beta_{\mathrm{QSO}}$ comes from the Ly$\alpha$-QSO
cross-correlation. However, a measurement of $\beta_{\mathrm{Ly}\alpha}$ (from
the auto) is necessary to constrain it due to their degeneracy in the model of
the cross-correlation. While we treat $\beta_{\mathrm{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.

##  Appendix C Other scale cuts

Figure 9: Comparison of the expected statistical constraint on $\phi$ from
eBOSS ($\sigma_{\phi}$) and the systematic bias ($\Delta\phi$) as a function
of maximum scale fitted, $r_{\mathrm{max}}$. The $\phi$ parameter here is
fitted from the full shapes of all four Ly$\alpha$ correlation functions
($\phi_{f}$). 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. Figure 10: Same as Figure 9, but as a
function of the minimum scale across the line-of-sight,
$r_{\bot}^{\mathrm{min}}$. We are progressively cutting the columns along the
line-of-sight, where the impact of the contaminants is highest. The systematic
error on $\phi$ is increasing as we cut more data, which indicates that we are
not able to accurately constrain the nuisance parameters that account for the
contaminants when ignoring the information along the line-of-sight.

In Section 3.2 we have focused on the minimum scale fitted,
$r_{\mathrm{min}}$. However, we also have to choose the largest scale we fit
for, $r_{\mathrm{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,
$r_{\mathrm{max}}$, in Figure 9. We perform the same type of comparison as we
did for $r_{\mathrm{min}}$ in Figure 4. The scales we tested range from $140$
$h^{-1}\text{Mpc}$ to $200$ $h^{-1}\text{Mpc}$. We find that measurements of
$\phi_{f}$ are unbiased for all the $r_{\mathrm{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 $r_{\mathrm{max}}=180$
$h^{-1}\text{Mpc}$ in our standard analysis.

As most of the contaminants of the Ly$\alpha$ forest correlation functions
have their largest impact along the line-of-sight, we also decided to test a
cut in $r_{\bot}^{\mathrm{min}}$. This is the smallest transverse scale that
we fit, and values $r_{\bot}^{\mathrm{min}}>0$ $h^{-1}\text{Mpc}$ mean that we
are removing the columns in the correlation that lie along the line-of-sight.
The systematic error and expected statistical constraint on $\phi_{f}$ for
different values of $r_{\bot}^{\mathrm{min}}$ are show in Figure 10. We find
that removing the first correlation column along the line-of-sight
($r_{\bot}^{\mathrm{min}}=4$ $h^{-1}\text{Mpc}$) still results in an unbiased
measurement. However, removing more data after that starts to introduce a
significant systematic bias.232323Note that we are plotting the absolute value
of $\Delta\phi$. Therefore, the trend of increasing systematic bias in Figure
10 does not mean the bias is all in the same direction. This could be
explained by the fact that removing the columns along the line-of-sight is
degrading our ability to fit the parameters that control the contaminants
(e.g. HCDs and metals). Even though these contaminants have their strongest
impact along the line-of-sight, they have a small impact on the correlation
even at large transverse separations. Therefore, it is still important to
model them accurately. Given these results, we chose to not impose any cut in
$r_{\bot}^{\mathrm{min}}$.

##  Appendix D Population constraints

Figure 11: Histograms of the $68\%$ constraints on $\alpha_{\mathrm{p}}$,
$\phi_{\mathrm{p}}$ and $\phi_{\mathrm{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.

As described in Section 3.4, when running the analysis on each of the one
hundred mocks, we are restricted to using a minimizer due to computational
constraints. This means that the reported constraints on the measured
parameters are just estimates due to two reasons. Firstly, they assume the
posterior is Gaussian around the best fitting value. Secondly, the best fit
for the parameters of interest ($\phi_{f}$, $\alpha_{\mathrm{p}}$) is found
without properly marginalizing over the nuisance parameters.242424The best fit
value of $\phi_{\mathrm{s}}$ reported by the minimizer is conditional on the
best fitting values of all the nuisance parameters. However, in the Bayesian
framework, the quantity we are interested in is the value of
$\phi_{\mathrm{s}}$ conditional on the posterior distributions of the nuisance
parameters. Only in the second case is the uncertainty due to these parameters
properly taken into account. While these approximations might be correct for
the BAO parameters due to their lack of correlations with the nuisance
parameters (Cuceu et al., 2020), $\phi_{\mathrm{s}}$ definitely has such
correlations (Appendix B). Therefore, in order to study accurate parameter
constraints from the population of mocks, we need to use a sampler.

For the purposes of this appendix, we ran the sampler on correlations computed
only using the Ly$\alpha$ region (no Ly$\beta$ region) from uncontaminated
synthetic data sets. These mocks do not include the effects of HCDs, metal
absorption or redshift errors (they are also described in du Mas des Bourboux
et al., 2020). Therefore, the parameter space is much smaller, allowing for
the computation of the full posterior distributions using a sampler. We
performed the joint analysis of the auto and cross-correlations for all one
hundred uncontaminated mocks.

In Figure 11, we show histograms of the constraints on $\alpha_{\mathrm{p}}$,
$\phi_{\mathrm{p}}$, and $\phi_{\mathrm{s}}$ from the population of one
hundred mocks. The main result is the one given by the sampler, in blue. We
also plot the standard deviation of the best fitting values for each of the
parameters over the population of mocks. This shows that the variance of the
parameter values obtained from the one hundred realizations is consistent with
the measured uncertainty in individual realizations, which means the
uncertainties are properly accounted for.

We also compare results from the sampler with the estimates from the minimizer
in Figure 11. While they roughly agree for the BAO parameters
($\alpha_{\mathrm{p}}$ and $\phi_{\mathrm{p}}$), this is not the case for
$\phi_{\mathrm{s}}$. The two give different constraints, as shown by the
difference between the red and blue histograms. This is not surprising given
the assumptions involved in the fitter constraints. However, the variance of
the $\phi_{\mathrm{s}}$ best fits is also different between the two (the red
and blue vertical lines). This means that marginalization over the nuisance
parameters is important in this case.252525In a frequentist framework, a
possible solution for this would be to use these results to compute non-
Gaussian $\Delta\chi^{2}$ values in order to set more accurate confidence
intervals on $\phi_{\mathrm{s}}$.

While this has been a rough test due to computational constraints, most of the
parameters that $\phi_{\mathrm{s}}$ is correlated with are still present in
uncontaminated mocks ($b_{\mathrm{Ly}\alpha}$, $\beta_{\mathrm{Ly}\alpha}$ and
$\beta_{\mathrm{QSO}}$). The main omission is $b_{\mathrm{HCD}}$. Therefore,
it would be interesting to perform this analysis with the effect of HCDs
modelled in the mocks. Alternatively, future studies should perform the
analysis on the fully contaminated mocks if there is access to enough
computational resources.