# Don’t miss the forest for the trees: the Lyman alpha forest power spectrum
in effective field theory

Mikhail M. Ivanov11footnotetext<EMAIL_ADDRESS>
###### Abstract

We derive an effective field theory (EFT) for cosmological Lyman alpha forest
fluctuations valid for the power spectrum at the one-loop order. The “bottom-
up” EFT expansion at the level of the transmitted flux is identical to the
line-of-sight dependent bias model first derived by Desjacques et al. We
confirm this result by a “top-down” derivation based on the exponential map of
the optical depth field. Specifically, we show that the combination of the
exponential map and conditions of renormalizability generates the same EFT
expansion as the “bottom-up” approach. In passing, we point out
inconsistencies of the tree-level perturbative expansion of the exponential
map without counterterms. To facilitate practical applications, we generalize
the FFTLog method for efficient calculations of one-loop integrals from new
line-of-sight dependent operators. Finally, we compare the one-loop EFT model
against data from the Sherwood hydrodynamic simulation. The theory fits the
data with sub percent accuracy up to $k=$ 3 $h$Mpc-1 at $z=2.8$ for both 3D
and 1D correlations. Our model can be readily used for cosmological full-shape
analyses of Lyman alpha forest data.

MIT-CTP/5607

##  1 Introduction

The Lyman alpha forest is the key source of information on the matter power
spectrum at small scales ($\sim$ Mpc and smaller) and high redshifts,
$2\lesssim z\lesssim 5$ [1, 2, 3, 4, 5, 6]. It has been extensively used to
constrain the $\Lambda$CDM model [7, 8, 9, 10, 11, 12, 13, 14, 15], as well as
neutrino masses [12, 16, 17, 18, 19], primordial black hole dark matter [20],
dark-matter baryon interactions [21], and other non-minimal cosmological
models [22, 23, 24]. The forest is also a unique probe of warm dark matter
[25, 26, 27, 28, 29, 30, 31, 32, 33, 34] and fuzzy dark matter models [35, 36,
37] that is not affected by ambiguities of the galaxy formation physics. In
addition, the large-scale Lyman alpha forest is a valuable source of high-
redshift baryon acoustic oscillation (BAO) information, see e.g. [38, 39, 40,
41]. In addition to the BAO, the broadband shape of the large-scale 3D Lyman
alpha power spectrum may also be a powerful probe of cosmology [42, 43, 44,
45].

The Lyman alpha forest has been one of the key targets for past and ongoing
large-scale structure surveys such as BOSS [46, 15], eBOSS [41], XQ-100 [14],
MIKE/HIRES [28], and DESI [47, 48, 45]. Entering the regime of high-precision
Lyman alpha cosmology with DESI requires a robust theoretical understanding of
the forest. The standard approach to the forest modeling has been
hydrodynamical simulations, see e.g. [49, 50, 51, 52, 53, 54, 55, 56]. In this
work we explore an alternative method and study to what extent the Lyman alpha
fluctuations can be modeled perturbatively [57, 22, 19, 58, 59, 60]. On the
one hand, perturbative approaches break down at short scales, where the power
spectrum is fully nonlinear, which ultimately limits their utility. On the
other hand, they provide a high level of accuracy and flexibility on mildly
non-linear scales. Indeed, as long as the EFT formally applies,222The
conditions of the applicability of the EFT will be specified shortly. it gives
a systematic program of consecutive approximations that can be executed to
arbitrary accuracy. The high flexibility of the EFT is especially important
for the efficient exploration of beyond-$\Lambda$CDM models. For example, it
is currently unfeasible to run a full simulation-based Monte Carlo analysis of
a large grid of $\Lambda$CDM extensions similar to the analysis by the Planck
collaboration [61].333 Although recently there has been significant progress
in cosmological analyses of Lyman alpha forest beyond $\Lambda$CDM based on
approximating, interpolating [17, 62, 63] and emulating [64, 65, 66] a grid of
simulations. Such an analysis, is, however, perfectly possible with
perturbation theory-based pipelines analogous to the ones recently applied to
galaxy clustering data, see e.g. [67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77,
78]. In addition, perturbative descriptions are based on different assumptions
than simulations. In principle, they provide alternative first-principle
models that are agnostic about the physics of intergalactic medium.444See [79,
80] for relationship between the underlying physics of the Lyman alpha forest
and perturbation theory parameters. Therefore, they represent a valuable
addition to the Lyman alpha cosmology toolbox.

The aim of this work is to develop an effective field theory (EFT) description
for the Lyman alpha forest, analogous to that of other large-scale structure
tracers [81, 82, 83, 84, 85]. The main philosophy of the EFT is to describe
dynamics on large scales using only relevant symmetries. These ideas have been
explored before both in the context of the Lyman alpha forest [22, 58, 59] and
for general line-of-sight dependent bias tracers [86]. In particular, our work
is a direct extension of Refs. [86, 59]. Desjacques et al. [86] derived a
general perturbative model for a tracer of matter that depends on the line-of-
sight selection effects. Although this description has been developed
primarily in the context of galaxies [87], one can argue that it applies,
practically without modifications, to the Lyman alpha forest. This is quite
natural, as the EFT bias expansion [83] is based on symmetries, and the
galaxies in the presence of selection effects enjoy the same key symmetries as
the Lyman alpha forest: the $SO(2)$ rotations around the line-of-sight and the
equivalence principle.

Chen et al. [59] pursued an alternative route and built a perturbative
expansion for the Lyman alpha forest flux $F$ starting from the exponential
map of the effective optical depth field $\tau$, $F=\exp(-\tau)$. This field
is often assumed to be a selection-independent tracer of dark matter, whose
density is conserved during transformations from real to redshift space, see
e.g. [79, 80]. Chen et al. then showed how the combination of dynamics and the
exponential map produces terms in the anisotropic bias expansion for the Lyman
alpha forest derived in [83]. In this regard, the main theoretical goal of
this work is to establish a direct connection between the two approaches. We
will explicitly show that the combination of the exponential map and the EFT
renormalizability principle explicitly reproduces the complete set of EFT
operators. In order to emphasize the role of renormalization, we also carry
out an alternative derivation in which we ignore counterterms, i.e. treat the
perturbative expansion à la standard perturbation theory (SPT) [88]. In that
case the exponential map produces an incomplete set of line-of-sight dependent
bias operators at the tree-level, with strong constraints on their bias
coefficients. We call this naive SPT-like expansion “the tree-level Tau
model.” We show explicitly that this model is mathematically inconsistent.
Once the model is properly renormalized, its tree-level constraints get
violated by loops. We argue that within this particular derivation the loops
should also generate the rest of the EFT operators allowed by symmetries. Our
calculations thus suggest that the “renormalized” SPT-like approach is
equivalent to the full EFT555At least at the level of the one loop power
spectrum.  once the loop corrections are taken into account.

It is also important to note that Refs. [89, 90, 53] found that the optical
depth is actually a selection-dependent tracer, i.e. in contrast to idealized
galaxies666From now on when discussing galaxies, we will always assume an
idealized situation when the line-of-sight selection effects [87, 86] are
negligible, although this assumptions may not always be the case, see e.g.
[91, 92].  its fluctuations trace the line-of-sight velocity gradients. This
further strengthens the motivation for a bottom-up approach such as the EFT.

On a more practical side, we extend a fast logarithmic Fourier transform
(FFTLog) algorithm of Ref. [93, 94] to quickly compute the full EFT one-loop
power spectrum of the 3D Lyman alpha forest. Having done this calculation, we
compare our EFT power spectrum model to the Lyman alpha forest data from the
Sherwood simulations [52]. We find this model to be sub percent level accurate
up to $k_{\rm max}=3~{}h\text{Mpc}^{-1}$ for $z=2.8$, and up to $k_{\rm
max}=5~{}h\text{Mpc}^{-1}$ for $z=3.2$. We have also detected strong
deviations from the predictions of the tree-level Tau model. This may be
considered an “experimental” evidence in favor of the full EFT expansion and
validity of the renormalization program. Finally, we show how the EFT approach
works out at the level of the one-dimensional flux power spectrum. The methods
that we have presented are ready to be applied to data. This paves the way for
systematic and efficient EFT-based cosmological analyses of the Lyman alpha
forest.

Our paper is structured as follows. In Sec. 2 we recap the Lyman alpha forest
physics and discuss its aspects relevant for the EFT: scales and power
counting. Section 3 outlines the EFT model first derived by [83] and discusses
its features relevant for the Lyman alpha. There we also discuss the FFTLog
implementation of our one-loop calculation. Section 4 establishes the
relationship between the EFT model and the exponential map. There we also
introduce the tree-level Tau model, which we use as a case study to underline
the importance of loop corrections and renormalization. In Sec. 5 we compare
the full one-loop EFT model with Sherwood simulations. Sec. 6 is devoted to
the one-dimensional flux power spectrum modeling. Finally, we outline our main
results and draw conclusions in Sec. 7. Some technical material is collected
in appendices.

##  2 Preliminaries

###  2.1 Basics of Lyman alpha forest

Our Universe is filled with optically thin neutral hydrogen (HI) clouds at
redshifts $2\lesssim z\lesssim 5$. Background quasars emit UV radiation that
propagates towards us through these clouds along given lines of sights. The
atoms in the could absorb the background radiation if the radiation’s rest-
frame frequency happens to match that of a transition between hydrogen levels.
Notably, if the relevant wavelength is 121.6 nm, a hydrogen atom undergoes a
Lyman alpha transition. Since the quasar spectrum is continuous in a wide
frequency range, and there are many clouds along the line of sight at
different redshifts, the observed quasar spectrum features a dense “comb” of
absorption lines, called the Lyman alpha forest. The fraction of transmitted
flux for each quasar spectrum is not homogeneous, and its fluctuations trace
cosmological mass fluctuations along the line of sight. Specifically, under
the assumption of photoionization-recombination equilibrium, the optical depth
of neutral hydrogen clouds is proportional the Lyman alpha absorption cross
section and the neutral hydrogen density $n_{HI}$ along the line of sight. The
latter scales with the fractional overdensity of gas (baryons) $\delta_{b}$,
gas temperature $T$, and the ionizing background amplitude $J$ as

$n_{HI}\propto\frac{(1+\delta_{b})^{2}}{T^{0.7}J}\,.$ (2.1)

Assumptions of the adiabatic expansion and the photoionization equilibrium
lead to the tight relation between temperature and density,
$T=T_{0}(1+\delta_{b})^{\gamma-1}$, where $\gamma\simeq 1.6$. Assuming that
$T_{0}$, $J$ and $\gamma$ are constants (i.e. do not have spatial
fluctuations), we arrive at the well-known conclusion (see e.g. [95, 35, 96])
that the fluctuations of $n_{HI}$ observed through the forest trace the
underlying fluctuations of gas density, which in turn, must reflect the matter
over-density $\delta$.777Recall that in standard linear cosmological
perturbation theory $\delta_{b}=\delta$ sufficiently deep into the matter
domination regime. Note that this condition is not always correct beyond
$\Lambda$CDM, see e.g. [97, 98]. In this simplified model, known as the
fluctuating Gunn-Peterson approximation [99, 4], the optical depth $\tau$ is
proportional to $n_{HI}$, while the transmitted flux we measure is simply
given by $F=\exp(-\tau)$.

One way to study the Lyman alpha fluctuations is to model the optical depth
field. In our work, however, we employ a different approach. Since we are
interested in fluctuations of the flux, we focus directly on this observable.
Instead of modeling explicitly underlying physics that affects the flux, we
will build a perturbative (gradient) expansion that is based on symmetries of
the problem and a minimal set of additional assumptions. In particular, we
will not explicitly employ relations such as Eq. (2.1), but simply assume that
the large-scale flux contrast is a certain unknown function of only a few
variables (degrees of freedom): the tidal field, velocity gradients, and
stochastic noise. This function is then represented as a general Taylor
expansion that contains all possible dependencies allowed by symmetries. This
approach in general is known as the EFT of large-scale structure [81, 82, 84,
85].

In the EFT of large-scale structure one builds a perturbative expansion of
relevant observables in terms of the matter over field, which is linear on
large scales. This linear field $\delta^{(1)}$ is the seed function of all
perturbative EFT calculations. Its statistical properties are fully fixed it
terms of the linear matter power spectrum,

$\langle\delta^{(1)}({\bf k})\delta^{(1)}({\bf
k}^{\prime})\rangle=(2\pi)^{3}\delta_{D}^{(3)}({\bf k}+{\bf k}^{\prime})P_{\rm
lin}(k,z)\,.$ (2.2)

In what follows we will often drop the explicit redshift-dependence and assume
that all quantities are evaluated at a given redshift of the Lyman alpha
forest, which we take to be $z=2.8$ to match the simulation data that we use
here. This redshift is similar to the redshifts of Lyman alpha forest data
from BOSS [15], eBOSS [100] and DESI [47].

###  2.2 Relevant scales

The EFT is a gradient expansion based on the assumption of scale separation.
Let us discuss how well this assumption is satisfied by the Lyman alpha forest
fluctuations.

The first relevant scale in the Lyman alpha forest physics is the non-linear
scale where the density field becomes fully non-linear. In momentum space, it
can be estimated as a wavenumber for which the amplitude of the dimensionless
linear matter power spectrum becomes unity,

$\frac{k^{3}_{\rm NL}}{2\pi^{2}}P_{\rm lin}(k_{\rm
NL},z)=1\,,~{}~{}~{}\Longrightarrow\quad k_{\rm NL}\simeq
5~{}h\text{Mpc}^{-1}\,,\quad\text{for}\quad z\simeq 3\,.$ (2.3)

A perturbative description of the Lyman alpha forest is possible only if $k\ll
k_{\rm NL}$. The second relevant parameter is the smoothing scale. The forest
is smoothed by two effects888It is worth mentioning that the actual data
features certain observational effects that effectively act as an additional
source of smoothing, e.g. an imperfect resolution [17]. These effects can also
be absorbed, to some extent, into the EFT expansion.  : gas pressure (3D
effect) [101] and thermal broadening (1D effect) [102, 103, 8]. These
smoothing effects are naturally incorporated within the EFT as part of the
gradient expansion. The gas perturbations are smoothed by pressure on the
Jeans scale $k_{J}\sim 18~{}h\text{Mpc}^{-1}$ [101, 19]. The absorption lines
are also subject to thermal broadening, whose characteristic scale $k_{S}\sim
11~{}h\text{Mpc}^{-1}$ in the 1D flux power spectrum is comparable to that of
the Jeans smoothing [19].999We stress that both scales are not known exactly.
The values quited in [19] are benchmark values that are sufficient for our
order-of magnitude estimates.  Both smoothing length scales are typically
shorter than the non-linear scale at the lower end of redshifts relevant for
the forest, $z\simeq 3$. This gives us the following hierarchy of scales in
the EFT for Lyman alpha:

$k\ll k_{\rm NL}\lesssim k_{J},~{}k_{s}\,.$ (2.4)

High redshift Lyman alpha forest is affected by potentially significant
spatial fluctuations of the ionizing background, and temperature fluctuations
due to inhomogeneous (patchy) reionization.101010In our notation this can be
thought of a spatial variation of $T_{0}$ in Eq. (2.1). The first effect
modulates the power spectrum shape on a broad range of scales, up to 100
$h^{-1}\text{Mpc}$ [104, 105]. These modulations are, however, quite smooth,
which suggest that they can be captured perturbatively with an appropriately
modified bias expansion. As far as the patchy reionisation is concerned, Ref.
[106] suggests that the gas temperature fluctuates on a scale of $\sim
25~{}h^{-1}\text{Mpc}$ at $z\sim 5$. These results are consistent with
simulations of Ref. [107], which additionally imply that patchy reionisation
may impact the 3D Lyman alpha power spectrum on wavenumber as large as
$0.2~{}h\text{Mpc}^{-1}$ at a few percent level at low redshifts. This
suggests that the typical scales associated with UV background and temperature
fluctuations are significantly longer than the non-linear scale, and hence
their systematic description within a gradient expansion approach, such as
EFT, may present a theoretical challenge111111See however Ref. [108] for an
approach to resumm the large gradients due to the radiation transfer effects..
In this work we build the EFT for Lyman alpha forest starting with an
idealized scenario where the photoionization and temperature fluctuations are
absent. We stress though that our description is still realistic at low
redshifts ($z\sim 3$) and small scales ($k\gtrsim 0.1~{}h\text{Mpc}^{-1}$),
where the effects in question are suppressed [105, 107].

Within our approximation, the most relevant scale for the EFT expansion is the
non-linear scale, just like in the usual EFT of LSS. In what follows we will
construct a perturbative expansion of the Lyman alpha forest fluctuations
around linear theory, which will be implicitly controlled by a small parameter
$k/k_{\rm NL}$. Let us sketch the structure of this expansion.

###  2.3 Estimates of perturbative terms in the scaling Universe

Figure 1: The dimensionless linear matter power spectrum in $\Lambda$CDM and a
scaling Universe from Eq. (2.5). Here we assume $z=2.8$, a redshift of the
Sherwood simulation snapshot that we use in this work.

Let us estimate the size of various terms in perturbation theory. It is
convenient to approximate the actual $\Lambda$CDM power spectrum with a power-
law $P_{\rm lin}\propto k^{n}$. For the range of scales relevant for our
analysis, and for $z\approx 3$, $k\approx 3~{}h\text{Mpc}^{-1}$, the linear
matter power spectrum can be well approximated as

$P_{\rm lin}(k)=2\pi^{2}\frac{k^{n}}{k^{n+3}_{\rm NL}}\,,\quad\text{with}\quad
n\approx-2.5\,,\quad k_{\rm NL}\approx 5.3~{}h\text{Mpc}^{-1}\,,$ (2.5)

see Fig. 1. From the dimensional analysis, the $L$’th-loop correction to this
power spectrum must scale as [88, 109]

$\frac{k^{3}}{2\pi^{2}}P_{L-\text{loop}}(k)=\frac{\alpha_{L}(2\pi^{2})}{k^{3}_{\rm
NL}}\left(\frac{k}{k_{\rm NL}}\right)^{(n+3)(L+1)}\,,$ (2.6)

where $\alpha_{L}$ is an order one number. These terms will be refereed to as
“mode-coupling corrections” in what follows. The higher derivative
(counterterm) corrections scale as $k^{2}P_{\rm lin}\propto k^{n+2}$. Thus,
the total power spectrum in the EFT reads, schematically

$\begin{split}P_{\rm NL}(k)&=2\pi^{2}\frac{k^{n}}{k^{n+3}_{\rm
NL}}\left(\underbrace{1}_{\text{linear}}+\underbrace{\alpha_{1}\left(\frac{k}{k_{\rm
NL}}\right)^{n+3}}_{\text{1-loop}}+\underbrace{\alpha_{2}\left(\frac{k}{k_{\rm
NL}}\right)^{2(n+3)}}_{\text{2-loop}}+\underbrace{\alpha_{k^{2}}\left(\frac{k}{k_{\rm
NL}}\right)^{2}}_{\text{higher deriv.}}+...\right)\\\
&=2\pi^{2}\frac{k^{-2.5}}{k^{0.5}_{\rm
NL}}\left(\underbrace{1}_{\text{linear}}+\underbrace{\alpha_{1}\left(\frac{k}{k_{\rm
NL}}\right)^{0.5}}_{\text{1-loop}}+\underbrace{\alpha_{2}\left(\frac{k}{k_{\rm
NL}}\right)^{1}}_{\text{2-loop}}+\underbrace{\alpha_{k^{2}}\left(\frac{k}{k_{\rm
NL}}\right)^{2}}_{\text{higher deriv.}}+...\right)\,.\end{split}$ (2.7)

where we explicitly included only the mode-coupling and $k^{2}$-corrections.
We see that in our toy model scaling Universe the one-loop corrections due to
mode coupling are the leading contributions in the regime $k\ll k_{\rm NL}$.
The $k^{2}$ corrections become important only at the three loop
order.121212See [110, 111] for the three loop order matter calculations.

Our discussion so far has considered non-linearity in the context of the
matter density field. In the EFT, additional non-linearities appear due to
biasing and redshift-space distortions. These effects are conceptually similar
to the dark matter non-linearities, and hence our scaling universe estimates
are expected to apply to the more realistic case without significant
modifications.131313It is important to note that the velocity field is more
non-linear than the density one, see e.g. [112, 113], i.e. the wavenumber
associated with the velocity non-linear scale is smaller. One aspect worth
mentioning is the presence of constant stochastic shot noise contribution,
which is, however, negligibly small in the scaling Universe where it scales as
$P_{\rm shot}\sim k_{\rm NL}^{-3}$.

From Eq. (2.7) we see that in contrast to galaxies (see e.g. [114, 115]), the
leading order corrections to linear theory in the Lyman alpha forest case come
from the one-loop mode coupling effects.141414Here we distinguish between the
one-loop mode-coupling integrals and higher-derivative counterterms, both of
which are parts of the one-loop power spectrum. Recall that for the scales
relevant to galaxy clustering the effective power law is less steep,
$n\approx-1.5$, so that the $k^{2}-$type higher derivative terms become
important already at the one-loop level. Our $\Lambda$CDM Universe is, of
course, more complicated than the toy model we consider here. Nevertheless,
the scaling universe model gives us a good sense of the importance of various
terms in the perturbation expansion.

##  3 Effective field theory for Lyman alpha fluctuations

###  3.1 Overview of galaxies in redshift space

It is instructive to start our discussion of the bias model for the Lyman
alpha with a recap of the bias model for galaxies. In real space, and on large
scales, the galaxy overdensity field $\delta_{g}$ can be expressed through a
perturbative expansion over the velocity gradients and tidal fields [116, 117,
118, 119, 83]. This expansion is, in general, non-local in time, but at the
cubic in density order, relevant for the one-loop power spectrum calculation,
it can be expressed through local in time operators [83]. Specifically, at
this order, we have the following expression for the deterministic part of the
galaxy bias relation,

$\begin{split}\delta_{g}({\bf x})=b^{g}_{1}\delta({\bf
x})+\frac{b^{g}_{2}}{2}(\delta({\bf x})^{2}-\langle\delta({\bf
x})^{2}\rangle)+b^{g}_{\mathcal{G}_{2}}\mathcal{G}_{2}({\bf
x})+b^{g}_{\Gamma_{3}}\Gamma_{3}({\bf
x})+b^{g}_{\nabla^{2}\delta}R_{*}^{2}\nabla^{2}\delta\,,\end{split}$ (3.1)

where $\delta$ is the (nonlinear) matter overdensity field, and we also
defined

$\begin{split}&\mathcal{G}_{2}({\bf
x})=(\partial_{i}\partial_{j}\Phi)^{2}-(\partial^{2}\Phi)^{2}\,,\\\
&\Gamma_{3}({\bf
x})=\mathcal{G}_{2}[\Phi]-\mathcal{G}_{2}[\Phi_{v}]\,,\end{split}$ (3.2)

where we introduced velocity and density potentials $\Phi$ and $\Phi_{v}$,
respectively, which satisfy

$\partial^{2}\Phi=\delta\,,\quad\partial_{i}\Phi_{v}=v^{i}\,,$ (3.3)

where $v^{i}$ is the peculiar velocity field. The bias coefficients
$b^{g}_{1},b^{g}_{2},b^{g}_{\mathcal{G}_{2}},b^{g}_{\Gamma_{3}},b^{g}_{\nabla^{2}\delta}$
are free parameters that characterize a given galaxy selection; $R_{*}$ is the
typical length scale associated with galaxies, which we inserted to make
$b^{g}_{\nabla^{2}\delta}$ dimensionless.

There are several aspects of the galaxy bias model in Eq. (3.1) that are worth
emphasizing. It is the most general local in time perturbative expansion in
the matter density in real space, modulo terms that do not contribute at the
level of the one-loop galaxy power spectrum, e.g. $\delta^{3}$. This expansion
involves all possible operators allowed by symmetries of the problem, the
spherical symmetry and the equivalence principle. In our context this means
the relevant operators are scalars under $SO(3)$ rotations, and there is no
velocity bias on large scales, i.e. galaxies undergo the same acceleration as
dark matter. Using the general symmetry-motivated expansions like (3.1) is the
key principle of effective field theory. In this philosophy, the distribution
of galaxies on large scales is expressed in terms of underlying basic long-
wavelength degrees of freedom, which in our context are matter density,
velocity gradients, and tidal fields, plus stochastic terms that will be
discussed later. The tree-level bias coefficients that appear in (3.1) are,
strictly speaking, not well defined since they receive formally infinite loop
corrections. In the EFT jargon these are called “bare” parameters, or Wilson
coefficients. Only their finite “renormalized” parts make physical sense and
can be matched to the data that describes physical observables.

Since the galaxies are observed in redshift space, one needs to transform the
galaxy bias expansion (3.1) using the following velocity-dependent mapping
(see e.g. [88]),151515Using this expansion, we assume a flat sky approximation
that is accurate for small scales relevant for our work.

$\delta^{(s)}_{g}({\bf k})=\delta_{g}({\bf k})+\int d^{3}x~{}e^{-i{\bf k}{\bf
x}}\left(e^{-ik_{z}v_{gz}({\bf x})/(aH)}-1\right)(1+\delta_{g}({\bf x}))\,,$
(3.4)

where we have switched to the Fourier-space representation, and used $a$ and
$H$ to denote the metric scale factor and the Hubble parameter, respectively,
${\cal H}=aH$ is the conformal Hubble parameter, whilst the subscript $z$
stands for the projection onto the line-of-sight, which we will describe with
the unit vector $\hat{\bf z}$. $v_{g}$ above is the tracer’s velocity, which
is equal to the matter peculiar velocity $v$ at the zeroth order in the
derivative expansion due to the equivalence principle.

In what follows it will be convenient to use the following notation for the
normalized velocity gradient of matter along the line-of-sight,

$\eta\equiv\frac{\partial_{z}v_{z}}{aH}\,.$ (3.5)

Note that in linear theory in Fourier space we have

$\eta\approx-f\mu^{2}\delta^{(1)}_{\bf k}\,,\quad\mu=\frac{(\hat{\bf
z}\cdot{\bf k})}{k}\,,$ (3.6)

where we introduced $f\equiv d\mathop{\rm ln}\nolimits D_{+}/d\mathop{\rm
ln}\nolimits a$ ($D_{+}$ is the growth factor).

Assuming that the velocity field is perturbative, one Taylor expands the
exponent in (3.4) and finds, at linear order,

$\delta_{g}^{(s)}=b_{1}^{g}\delta-\eta\,.$ (3.7)

Crucially, the second term, which depends on the line-of-sight, does not have
any free parameter in front of it. This embodies the absence of selection
effects, and follows from the conservation of the galaxy number density under
the coordinate transformation from the galaxy rest frame to observed redshift
space.

Going to higher order in the Taylor expansion of (3.4), one obtains
corrections to the observed galaxy overdensity due to redshift-space
distortions (RSD). Note that this expansion generates extra UV sensitivity,
which must be removed by appropriate counterterms [120, 121], which start at
order $k^{2}\delta$. Apart from these higher-derivative contributions, RSD
mapping (3.4) does not generate new free parameters for the deterministic part
of the galaxy power spectrum in redshift space. The redshift space galaxy
overdensity reads161616We use the notation $\int_{\bf
q}\equiv\int\frac{d^{3}q}{(2\pi)^{3}}$.

$\begin{split}\delta_{g}^{(s)}({\bf k})=&Z_{1}({\bf k})\delta^{(1)}_{\bf
k}+\int_{{\bf q}_{1}}\int_{{\bf q}_{2}}(2\pi)^{3}\delta_{D}({\bf k}-{\bf
q}_{12})Z_{2}({\bf q}_{1},{\bf q}_{2})\delta^{(1)}_{{\bf
q}_{1}}\delta^{(1)}_{{\bf q}_{2}}+\\\ &+\int_{{\bf q}_{1}}\int_{{\bf
q}_{2}}\int_{{\bf q}_{3}}(2\pi)^{3}\delta_{D}({\bf k}-{\bf q}_{123})Z_{3}({\bf
q}_{1},{\bf q}_{2},{\bf q}_{3})\delta^{(1)}_{{\bf q}_{1}}\delta^{(1)}_{{\bf
q}_{2}}\delta^{(1)}_{{\bf q}_{3}}+...\,,\end{split}$ (3.8)

where the dots denote higher loop corrections and higher derivate
counterterms, $\delta^{(1)}_{\bf k}$ is the linear matter density field, while
the redshift-space kernels $Z_{1,2,3}$ are given by

$\displaystyle Z_{1}({\bf k})=b^{g}_{1}+f\mu^{2}\,,$ (3.9a) $\displaystyle
Z_{2}({\bf k}_{1},{\bf
k}_{2})=\frac{b^{g}_{2}}{2}+b^{g}_{\mathcal{G}_{2}}\left(\frac{({\bf
k}_{1}\cdot{\bf k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}-1\right)+b^{g}_{1}F_{2}({\bf
k}_{1},{\bf k}_{2})+f\mu^{2}G_{2}({\bf k}_{1},{\bf k}_{2})$
$\displaystyle\qquad\qquad\quad~{}~{}+\frac{f\mu
k}{2}\left(\frac{\mu_{1}}{k_{1}}(b^{g}_{1}+f\mu_{2}^{2})+\frac{\mu_{2}}{k_{2}}(b^{g}_{1}+f\mu_{1}^{2})\right)\,,$
(3.9b) $\displaystyle Z_{3}({\bf k}_{1},{\bf k}_{2},{\bf
k}_{3})=2b^{g}_{\Gamma_{3}}\left[\frac{({\bf k}_{1}\cdot({\bf k}_{2}+{\bf
k}_{3}))^{2}}{k_{1}^{2}({\bf k}_{2}+{\bf
k}_{3})^{2}}-1\right]\big{[}F_{2}({\bf k}_{2},{\bf k}_{3})-G_{2}({\bf
k}_{2},{\bf k}_{3})\big{]}$ $\displaystyle\quad+b^{g}_{1}F_{3}({\bf
k}_{1},{\bf k}_{2},{\bf k}_{3})+f\mu^{2}G_{3}({\bf k}_{1},{\bf k}_{2},{\bf
k}_{3})+\frac{(f\mu
k)^{2}}{2}(b^{g}_{1}+f\mu_{1}^{2})\frac{\mu_{2}}{k_{2}}\frac{\mu_{3}}{k_{3}}$
$\displaystyle\quad+f\mu k\frac{\mu_{3}}{k_{3}}\left[b^{g}_{1}F_{2}({\bf
k}_{1},{\bf k}_{2})+f\mu^{2}_{12}G_{2}({\bf k}_{1},{\bf k}_{2})\right]+f\mu
k(b^{g}_{1}+f\mu^{2}_{1})\frac{\mu_{23}}{k_{23}}G_{2}({\bf k}_{2},{\bf
k}_{3})$ $\displaystyle\quad+b^{g}_{2}F_{2}({\bf k}_{1},{\bf
k}_{2})+2b_{\mathcal{G}_{2}}\left[\frac{({\bf k}_{1}\cdot({\bf k}_{2}+{\bf
k}_{3}))^{2}}{k_{1}^{2}({\bf k}_{2}+{\bf k}_{3})^{2}}-1\right]F_{2}({\bf
k}_{2},{\bf k}_{3})+\frac{b^{g}_{2}f\mu k}{2}\frac{\mu_{1}}{k_{1}}$
$\displaystyle\quad+b^{g}_{\mathcal{G}_{2}}f\mu
k\frac{\mu_{1}}{k_{1}}\left[\frac{({\bf k}_{2}\cdot{\bf
k}_{3})^{2}}{k_{2}^{2}k_{3}^{2}}-1\right]\,,$ (3.9c)

where the density and velocity SPT kernels $F_{n}$ and $G_{n}$ can be found,
e.g. in [88]. Note that $Z_{3}$ must be symmetrized over its arguments. We
also introduced cosines between the line-of-sight and momentum vectors,

$\mu_{i}\equiv\frac{(\hat{\bf z}\cdot{\bf k}_{i})}{k_{i}}\,,$ (3.10)

and use ${\bf k}$ to denote the sum of all arguments, i.e. for $Z_{3}$ we have
${\bf k}\equiv{\bf k}_{1}+{\bf k}_{2}+{\bf k}_{3}$.

Note that the density field variance $\langle\delta({\bf x})^{2}\rangle$ is
used in the Eq. (3.1) in order to enforce the vanishing of the overdensity
v.e.v., $\langle\delta\rangle=0$, which a properly defined fluctuating field
must satisfy. The ${\bf k}=0$ correction to (3.9) produced by the density
variance is not explicitly shown above, but it is needed to ensure that the
line-of-sight dependent bias operators are not generated by loop corrections
for galaxies.

###  3.2 General bias expansion for the Lyman alpha forest

The Lyman alpha forest is different from galaxies in several important
aspects. First, unlike galaxies, the forest has a preferred direction - the
line-of-sight. Indeed, in contrast to galaxies, the flux fluctuations
intrinsically exist only in redshift space. This means that new operators are
allowed in the Lyman alpha forest bias expansion. Physically, this happens
because the absorption probability for the Lyman-$\alpha$ line depends on the
tidal field along the line of sight. Mathematically, this property may traced
to the nonlinear mapping between flux and the optical depth,
$F=\exp{(-\tau)}$. If one assumes that the $\tau$ field is a matter tracer
whose number density is conserved, just like for galaxies, there would be no
selection effects, and the bias expansion would take the form (3.7) at the
lowest order. The nonlinear mapping, however, breaks the conservation of
tracer number density [79] and forces us to introduce line-of-sight selection
effects. At linear order, and in the rest frame of a neutral hydrogen cloud,
we should write

$\frac{F-\bar{F}}{\bar{F}}\equiv\delta_{F}=b_{1}\delta+b_{\eta}\hat{z}^{i}\hat{z}^{j}\partial_{i}\partial_{j}\Phi\,,$
(3.11)

where $\bar{F}\equiv\langle F\rangle$ is the mean transmitted flux. In linear
theory the tidal field is related to the velocity gradient $\eta$ defined in
Eq. (3.5), so that one can re-write the above expression as

$\delta_{F}=b_{1}\delta+b_{\eta}\eta\,,$ (3.12)

The new selection-dependent parameter $b_{\eta}$ will be referred to as
“velocity gradient bias” in what follows [80]. The linear model (3.12)
produces the well-known tree-level result

$P_{\rm tree}(k,\mu)=(b_{1}-b_{\eta}f\mu^{2})^{2}P_{\rm lin}(k)\,,$ (3.13)

which is a simple generalization of the Kaiser formula for galaxies [122]. We
stress that Eq. (3.11) is not the most general expression. As in the context
of galaxies, a consistent bias expansion should depend on the history of the
relevant operators along their fluid trajectories ${\bf x}_{\rm fl}({\bf
x},\tau)$ (${\bf x}$ is the last point of the trajectory) [123, 116, 118, 83,
85],171717In analogy with the smoothing scale of the Lyman alpha fluctuations
[101], one may say that the Lyman alpha bias parameters in general are
expected to depend on the whole reionization history of the Universe.

$\delta_{F}({\bf
x},\tau)\supset\int^{\tau}T(\tau,\tau^{\prime})\mathcal{O}(\tau^{\prime},{\bf
x}_{\rm fl}({\bf x},\tau^{\prime}))\,,$ (3.14)

where $\tau$ is conformal time. Using the matter equations of motion, the
dependence on the fluid trajectory and time can be Taylor expanded and summed
up into a finite number of linearly independent terms, which eventually
reproduce the local-in-time bias expansion similar to Eq. (3.12). All in all,
it is important to keep in mind that the bias expansion is intrinsically non-
local in time, but for the purposes of the one-loop power spectrum calculation
this non-locality is irrelevant, so we will proceed with the local in time
expansion.

It is straightforward to generalize the bias model (3.12) to operators that
are higher order in the linear density field. Specifically, Ref. [86] showed
that a general selection-dependent biased tracer has the following
deterministic bias model up to cubic order:

$\begin{split}\delta_{F}=&b_{1}\delta+\frac{b_{2}}{2}\delta^{2}+b_{\mathcal{G}_{2}}\mathcal{G}_{2}+b_{\Gamma_{3}}\Gamma_{3}\\\
&+b_{\eta}\eta+b_{(KK)_{\parallel}}K_{ij}K_{jl}\hat{z}^{i}\hat{z}^{l}+b_{\delta\eta}\eta\delta+b_{\eta^{2}}\eta^{2}+b_{\Pi^{[2]}_{\parallel}}\Pi^{[2]}_{ij}\hat{z}^{i}\hat{z}^{j}\\\
&+b_{\delta\Pi^{[2]}_{\parallel}}\delta\Pi^{[2]}_{ij}\hat{z}^{i}\hat{z}^{j}+b_{\eta\Pi^{[2]}_{\parallel}}\eta\Pi^{[2]}_{ij}\hat{z}^{i}\hat{z}^{j}+b_{(K\Pi^{[2]})_{\parallel}}K_{ij}\Pi^{[2]}_{jl}\hat{z}^{i}\hat{z}^{l}+b_{\Pi^{[3]}_{\parallel}}\Pi^{[3]}_{ij}\hat{z}^{i}\hat{z}^{j}\,,\end{split}$
(3.15)

where we kept only the terms that contribute non-trivially to the one-loop
power spectrum, and introduced the new operators

$\begin{split}&K_{ij}=\frac{\partial_{i}\partial_{j}\delta}{\Delta}-\frac{1}{3}\delta_{ij}\delta\\\
&\Pi^{[1]}_{ij}=\partial_{i}\partial_{j}\Phi=\frac{\partial_{i}\partial_{j}\delta}{\Delta}\,,\\\
&\Pi^{[n]}_{ij}=\frac{1}{(n-1)!}\left[(aHf)^{-1}\frac{D}{D\tau}\Pi^{[n-1]}_{ij}-\Pi^{[n-1]}_{ij}\right],\quad
n>1\,,\end{split}$ (3.16)

where we used the convective derivative
$\frac{D}{D\tau}=\partial_{\tau}+v^{i}\partial_{i}$. The explicit expressions
for the tensors $\Pi^{[2]}_{ij}$ and $\Pi^{[3]}_{ij}$ in perturbation theory
can be found in Ref. [86]. In total, we have 9 new selection bias coefficients
in addition to the 4 selection free ones.181818Our discussion so far does not
include higher derivative counterterms and stochastic bias contributions,
which will generate additional free parameters.

The expression (3.15) is written in the rest-frame of the tracer. In order to
describe the actual observations, this has to be transferred to the observer’s
frame using the standard RSD mapping (3.4) [83, 86]. Note that the physical
meaning of the RSD mapping here is somewhat different from that in the
galaxies’ case. Neglecting selection effects, the mapping simply acts on the
real-space bias expansion that captures the actual number density in tracer’s
rest frame. If the probability to detect a galaxy depends on line-of-sight
properties [87, 92], the line-of-sight selection effects should be treated as
part of the bias expansion that describes the apparent number density of the
tracer in its rest frame. The apparent density can be defined as the observed
galaxy distribution formally inversely transformed from redshift space to real
space. Obviously, it is different from the actual, rest-frame (real-space)
galaxy number density, which, in principle, can be reconstructed with better
selection criteria or better measurements [87]. Contrasting the actual number
density and the observed one in redshift-space, we may say that the RSD
mapping does not conserve the galaxy number density in the presence of
selection effects. In the Lyman alpha case, the line-of-sight dependent bias
expansion also captures the apparent number density. This number density is by
definition subject to the RSD mapping that conserves it. The conceptual
difference with galaxies is that the Lyman alpha fluctuations intrinsically
live along lines on sight, so that a notion of the actual rest frame density
does not make sense in this case. From the technical point of view, however,
the treatment of the mapping is similar: on has to apply the RSD mapping (3.4)
to the selection-dependent bias expansion (3.15).

When expanded, the RSD mapping will produce additional velocity-dependent
contributions into (3.15), most of which will be degenerate with the selection
terms that are already present. This will lead to an unobservable re-
definition of the selection bias coefficients. There will be, however, a few
extra operators that are not present in (3.15). These terms correspond to pure
projection effects, i.e. they explicitly depend on locally unobservable
tracer’s velocity. The coefficients in front of these terms are protected by
the equivalence principle and do not get renormalized by non-linear effects.
With these terms included, we have the net expression [86]:

$\begin{split}\delta^{(s)}_{F}=\delta_{F}[\text{from
Eq.~{}(\ref{eq:F_full})}]+\frac{v_{z}}{aH}\hat{z}^{i}\partial_{i}\left[b_{1}\delta+b_{\eta}\eta+b_{\Pi^{[2]}_{\parallel}}\Pi^{[2]}_{kl}\hat{z}^{k}\hat{z}^{l}\right]\,.\end{split}$
(3.17)

The final mode-coupling three-dimensional one-loop power spectrum that comes
from the expansion (3.17) is given by

$\begin{split}\Delta P^{\text{1-loop}}(k,\mu)\equiv&2\int_{\bf
q}K_{2}^{2}({\bf q},{\bf k}-{\bf q})P_{\text{lin}}(|{\bf k}-{\bf
q}|)P_{\text{lin}}(q)\\\ &+6K_{1}({\bf k})P_{\text{lin}}(k)\int_{\bf
q}K_{3}({\bf k},-{\bf q},{\bf q})P_{\text{lin}}(q)\,,\end{split}$ (3.18)

where we introduced the new selection-dependent redshift space kernels,

$\begin{split}K_{1}({\bf k})=&b_{1}-b_{\eta}f\mu^{2}\,,\\\ K_{2}({\bf
k}_{1},{\bf k}_{2})=&\frac{b_{2}}{2}+b_{\mathcal{G}_{2}}\left(\frac{({\bf
k}_{1}\cdot{\bf k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}-1\right)+b_{1}F_{2}({\bf
k}_{1},{\bf k}_{2})-b_{\eta}f\mu^{2}G_{2}({\bf k}_{1},{\bf k}_{2})\\\
&-fb_{\delta\eta}\frac{\mu_{2}^{2}+\mu_{1}^{2}}{2}+b_{\eta^{2}}f^{2}\mu_{1}^{2}\mu_{2}^{2}\\\
&+b_{1}f\frac{\mu_{1}\mu_{2}}{2}\left(\frac{k_{2}}{k_{1}}+\frac{k_{1}}{k_{2}}\right)-b_{\eta}f^{2}\frac{\mu_{1}\mu_{2}}{2}\left(\frac{k_{2}}{k_{1}}\mu_{2}^{2}+\frac{k_{1}}{k_{2}}\mu_{1}^{2}\right)\\\
&+b_{(KK)_{\parallel}}\left(\mu_{1}\mu_{2}\frac{({\bf k}_{1}\cdot{\bf
k}_{2})}{k_{1}k_{2}}-\frac{\mu_{1}^{2}+\mu_{2}^{2}}{3}+\frac{1}{9}\right)\\\
&+b_{\Pi^{(2)}_{\parallel}}\left(\mu_{1}\mu_{2}\frac{({\bf k}_{1}\cdot{\bf
k}_{2})}{k_{1}k_{2}}+\frac{5}{7}\mu^{2}\left(1-\frac{({\bf k}_{1}\cdot{\bf
k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}\right)\right)\,,\end{split}$ (3.19)
$\begin{split}&\int_{{\bf q}}K_{3}({\bf k},{\bf q},-{\bf q})P_{\rm lin}(q)\\\
=\,&b_{1}\int_{{\bf q}}F_{3}({\bf k},{\bf q},-{\bf
q})P_{\text{lin}}(q)-fb_{\eta}\mu^{2}\int_{{\bf q}}G_{3}({\bf q},-{\bf q},{\bf
k})P_{\text{lin}}(q)\\\ &\>+\int_{{\bf
q}}\left[1-\left(\hat{k}\cdot\hat{q}\right)^{2}\right]P_{\rm lin}(q)\\\
&\>\times\Bigg{\\{}\frac{4}{21}(5b_{\mathcal{G}_{2}}+2b_{\Gamma_{3}})\left[\left(\frac{({\bf
k}-{\bf q})\cdot{\bf q}}{|{\bf k}-{\bf
q}|q}\right)^{2}-1\right]-\frac{2}{21}fb_{\delta\eta}\left[\frac{3(k_{\parallel}-q_{\parallel})^{2}}{|{\bf
k}-{\bf q}|^{2}}+\frac{5q_{\parallel}^{2}}{q^{2}}\right]\\\
&\quad\>+\frac{4}{7}f^{2}b_{\eta^{2}}\frac{q_{\parallel}^{2}}{q^{2}}\frac{(k_{\parallel}-q_{\parallel})^{2}}{|{\bf
k}-{\bf q}|^{2}}+\frac{20}{21}b_{(KK)_{\parallel}}\left[\frac{({\bf
k}\cdot{\bf q}-q^{2})(k_{\parallel}-q_{\parallel})q_{\parallel}}{|{\bf k}-{\bf
q}|^{2}q^{2}}-\frac{1}{3}\frac{(k_{\parallel}-q_{\parallel})^{2}}{|{\bf
k}-{\bf
q}|^{2}}-\frac{1}{3}\frac{q_{\parallel}^{2}}{q^{2}}+\frac{1}{9}\right]\\\
&\quad\>+\frac{10}{21}b_{\Pi_{\parallel}^{[2]}}\frac{({\bf k}\cdot{\bf
q}-q^{2})}{|{\bf k}-{\bf
q}|^{2}}\frac{(k_{\parallel}-q_{\parallel})^{2}}{q^{2}}+\frac{10}{21}\left[b_{\delta\Pi_{\parallel}^{[2]}}-\frac{1}{3}b_{(K\Pi^{[2]})_{\parallel}}-fb_{\eta\Pi_{\parallel}^{[2]}}\frac{q_{\parallel}^{2}}{q^{2}}\right]\frac{(k_{\parallel}-q_{\parallel})^{2}}{|{\bf
k}-{\bf q}|^{2}}\\\
&\quad\>+\frac{10}{21}b_{(K\Pi^{[2]})_{\parallel}}\frac{({\bf q}\cdot{\bf
k}-q^{2})}{q|{\bf k}-{\bf
q}|}\frac{q_{\parallel}(k_{\parallel}-q_{\parallel})}{q|{\bf k}-{\bf
q}|}+\frac{10}{21}fb_{\Pi^{[2]}_{\parallel}}\frac{q_{\parallel}(k_{\parallel}-q_{\parallel})^{3}}{q^{2}|{\bf
k}-{\bf q}|^{2}}\\\
&\quad\>+(b_{\Pi_{\parallel}^{[3]}}+2b_{\Pi_{\parallel}^{[2]}})\Bigg{[}\frac{13}{21}\frac{{\bf
k}\cdot{\bf q}-q^{2}}{|{\bf k}-{\bf
q}|^{2}}\frac{q_{\parallel}(k_{\parallel}-q_{\parallel})}{q^{2}}-\frac{5\mu^{2}}{9}\left[\left(\frac{({\bf
k}-{\bf q})\cdot{\bf q}}{|{\bf k}-{\bf
q}|q}\right)^{2}-\frac{1}{3}\right]\Bigg{]}\\\
&\quad\>+\frac{2}{21}fb_{1}\left[5\frac{q_{\parallel}(k_{\parallel}-q_{\parallel})}{q^{2}}+3\frac{q_{\parallel}(k_{\parallel}-q_{\parallel})}{|{\bf
k}-{\bf
q}|^{2}}\right]-\frac{2}{7}f^{2}b_{\eta}\frac{q_{\parallel}(k_{\parallel}-q_{\parallel})}{q^{2}|{\bf
k}-{\bf
q}|^{2}}\left[(k_{\parallel}-q_{\parallel})^{2}+q_{\parallel}^{2}\right]\Bigg{\\}}\,.\end{split}$
(3.20)

Note that for convenience, we presented the $K_{3}$ integrand instead of the
kernel, as in this case the expression greatly simplifies.

###  3.3 Higher derivative counterterms and stochastic contributions

The line-of-sight dependent higher derivative and stochastic contributions
relevant for the one-loop power spectrum of a generic selection-dependent
tracer are also known [86],

$\begin{split}\delta_{F}^{(s)}\Big{|}_{\text{stoch+h.d.}}=\epsilon+b^{g}_{\nabla^{2}\delta}R_{*}^{2}\nabla^{2}\delta+b_{\eta}R_{*}^{2}\nabla^{2}(\beta_{\nabla^{2}\textbf{v}}\nabla^{2}\eta+\beta_{\partial^{2}_{||}\textbf{v}}\nabla^{2}\eta)\,,\end{split}$
(3.21)

where $\epsilon$ is the stochastic field uncorrelated with $\delta$. Note that
we do not explicitly add the redshift-space EFT counterterms $\sim
k^{2}\mu^{2n}\delta$ [120] in Eq. (3.21) as they are fully degenerate with the
higher-derivative operators that we already have. On the one hand, our scaling
Universe analysis suggests that the higher-derivative contributions scale as
three-loop order contributions, and hence may be ignored at the one-loop order
we are interested in. On the other hand, the precision of the Sherwood
simulation data is sufficient to nominally detect these terms even when their
Wilson coefficients are very small, of the order
$10^{-3}~{}[h^{-1}\text{Mpc}]^{2}$. The inclusion of these corrections
noticeably improves the fit at $k\sim 3~{}h\text{Mpc}^{-1}$, which is why we
prefer to keep them. The theoretical rationale behind this choice is that
these corrections can be treated as a proxy for higher loop contributions,
which can be partially absorbed into the counterterms (3.21). Thus, their
inclusion may indeed improve the fit even if two-loop and three-loop
contributions are not explicitly included.

As far as the stochastic corrections are concerned, their expression up to
$(k/k_{\rm NL})^{2}$ order is given by [124, 86]

$\langle\epsilon({\bf k})\epsilon({\bf
k}^{\prime})\rangle=(2\pi)^{3}\delta_{D}^{(3)}({\bf k}+{\bf
k}^{\prime})\left(P_{\text{shot}}+a_{0}\frac{k^{2}}{k_{\rm
NL}^{2}}+a_{2}\frac{k^{2}\mu^{2}}{k_{\rm NL}^{2}}\right)\,,$ (3.22)

where $P_{\rm shot}$, $a_{0,2}$ are dimension-full stochastic counterterms.
The shot noise corrections can be ignored for the 3D correlators, but are
important for the 1D power spectrum, see Section 6 for more detail.

###  3.4 IR resummation

Baryon acoustic oscillations are affected by large IR effects, which may be
resummed along the lines of Refs. [125, 126, 127, 128, 129, 130, 131]. These
effects are fully controlled by the displacement and velocity fields, which
are both locally unobservable. In the context of the Eulerian fluid
description, the structure of these IR sensitive terms is dictated by the
equivalence principle, which demands that the IR-enhanced contributions cancel
in the formal limit $\\{q_{i}\\}\to 0$, where $\\{q_{i}\\}$ denote a set of
loop momenta. The cancellation is inexact if some $q_{i}$’s are greater than
the BAO scale $r_{\rm BAO}\sim 110~{}h\text{Mpc}^{-1}$. In that case (IR)
resummation is needed. Since IR resummation is controlled by the IR divergent
terms whose form is completely fixed by the equivalence principle, it would be
unaffected by the line-of-sight selection bias operators, which are finite in
the IR limit by construction. This implies that IR resummation for a generic
line-of-sight dependent tracer would be identical to that of redshift space
galaxies [130], as was first pointed out in [86]. Indeed, one can explicitly
check that the new line-of-sight dependent operators do not produce any new IR
singularities in the soft limit, see e.g. (3.20).

Even though IR resummation can be straightforwardly included, we defer its
explicit implementation for future work. First, we main goal of this paper is
to study the performance of the EFT on small scales, $k\sim
3~{}h\text{Mpc}^{-1}$, where the BAO feature is completely washed out. Second,
the hydrodynamical simulations that we use here do not have enough large-scale
resolution for a significant detection of the BAO in the data. Given these
reasons, we believe that the non-linear effects on the BAO deserve a separate
dedicated analysis that goes beyond the scope of this paper.

###  3.5 Calculation of loop integrals with FFTLog

The one-loop expression (3.18) can be transformed to the form suitable for the
evaluation with the FFTLog method [93] (see also [132, 133]). The basic idea
of this method is to represent the linear matter power spectrum, which is an
input in all loop calculations, as a sum of the power-law functions. Then the
loop integrals for each power law can be done analytically, and a loop
calculation reduces to a matrix multiplication problem.

In redshift space, following Ref. [94] (see also [134]), one has to expand the
loop integrands over a basis of master operators involving contractions of the
line of sight vector $\hat{\bf z}$ and loop momenta. One technical difficulty
w.r.t. to the case of galaxies is that one needs such master expression up to
eight momenta instead of four in the galaxy case. This is quite
straightforward albeit the expressions are rather lengthly. The eventual
expressions for the relevant master integrals after taking all the necessary
contractions of the basis operators are given by

$\begin{split}&\int_{\bf q}\frac{(\hat{\bf z}\cdot{\bf
q})^{5}}{q^{2\nu_{1}}|{\bf q}-{\bf k}|^{2\nu_{2}}}=k^{3-2\nu_{12}}\cdot
k^{5}\mu(A_{5}+\mu^{2}B_{5}+\mu^{4}C_{5})\,,\\\ &\int_{\bf q}\frac{(\hat{\bf
z}\cdot{\bf q})^{6}}{q^{2\nu_{1}}|{\bf q}-{\bf
k}|^{2\nu_{2}}}=k^{3-2\nu_{12}}\cdot
k^{6}(A_{6}+\mu^{2}B_{6}+\mu^{4}C_{6}+\mu^{6}D_{6})\,,\\\ &\int_{\bf
q}\frac{(\hat{\bf z}\cdot{\bf q})^{7}}{q^{2\nu_{1}}|{\bf q}-{\bf
k}|^{2\nu_{2}}}=k^{3-2\nu_{12}}\cdot
k^{7}\mu(A_{7}+\mu^{2}B_{7}+\mu^{4}C_{7}+\mu^{6}D_{7})\,,\\\ &\int_{\bf
q}\frac{(\hat{\bf z}\cdot{\bf q})^{8}}{q^{2\nu_{1}}|{\bf q}-{\bf
k}|^{2\nu_{2}}}=k^{3-2\nu_{12}}\cdot
k^{8}(A_{8}+\mu^{2}B_{8}+\mu^{4}C_{8}+\mu^{6}D_{8}+\mu^{8}E_{8})\,,\end{split}$
(3.23)

where $\nu_{12}\equiv\nu_{1}+\nu_{2}$ and the master functions
$A_{n},B_{n},C_{n},D_{n},E_{n}$ for $n=5,6,7,8$ are given in Appendix A. The
rest of the master integrals are presented in Appendix A of [94]. With the
help of master integrals like (3.23), the total mode-coupling power spectrum
may be written as

$\Delta P^{\rm
1-loop}=\sum_{n=0}^{4}\mu^{2n}P^{(22)}_{n}+(b_{1}-b_{\eta}f\mu^{2})\sum_{n=0}^{2}\mu^{2n}P^{(13)}_{n}\,,$
(3.24)

where $P^{(22)}_{n}$ and $P^{(13)}_{n}$ depend on the biases that appear in
Eq. (3.17).

An important aspect of the FFTLog method is the choice of the so-called FFTLog
bias parameter $\nu$ which appears in the FFTLog power law approximation as
$P_{\rm lin}\propto k^{\nu}$, and hence it controls the convergence of the
loop integrals. Most of the 22-type integrals that we encounter converge for
$-1<\nu<1/2$. The lower end of this range corresponds to spurious IR
contributions that cancel in the power spectrum due to the equivalence
principle [135, 128]. If this IR contribution is taken into account, one can
choose an FFTLog bias in the range $-3<\nu<-1$. Certain bias integrals, like
the ones proportional to $b_{2}^{2}$, converge only for $-3<\nu<-3/2$. For
$\nu>-3/2$ they have formally UV divergent constant parts that behave as [117,
93]

$I_{\delta^{2}\delta^{2}}(0)\equiv\frac{1}{2}\int_{{\bf q}}P^{2}_{\rm
lin}(q)\,.$ (3.25)

In the case of galaxies the standard practice is to subtract these terms from
the final answer as they simply renormalize the constant shot noise
contribution. For the Lyman alpha, however, the intrinsic shot noise
contribution produced by tracer discreteness is negligibly small, while the
constant piece (3.25) is quite large for the physical $\Lambda$CDM power
spectrum,

$I_{\delta^{2}\delta^{2}}(0)\Big{|}_{z=2.8}\approx
28~{}[h^{-1}\text{Mpc}]^{3}\,.$ (3.26)

Note that this integral is saturated in the IR, so it would have the same
value even if we cut it off at $\Lambda\sim k_{\rm NL}$. Hence, Eq. (3.26)
represents a calculable, physical constant power spectrum contribution that is
deterministic by origin, and must be present in the model. Let us estimate its
magnitude. When discussing the Tau model below, we will get a natural estimate
$b_{2}\sim 0.3$, so that

$P_{\rm shot}\sim\frac{b_{2}^{2}}{2}\int_{{\bf q},|q|\leq k_{\rm
NL}}P^{2}_{\rm lin}(q)\approx b_{2}^{2}\times
I_{\delta^{2}\delta^{2}}(0)\approx 2.5~{}[h^{-1}\text{Mpc}]^{3}\,.$ (3.27)

This estimate is two orders of magnitude larger than the naive scaling
Universe result $P_{\rm shot}\sim k_{\rm NL}^{-3}\sim
10^{-2}~{}[h^{-1}\text{Mpc}]^{3}$, and also many orders of magnitudes larger
than the actual discreteness shot noise.

The constant terms like (3.27) are fully calculable. They are produced mostly
by perturbative modes and they do not require introducing new parameters as
counterterms, as their UV-sensitivity is very mild for $\Lambda$CDM. The
deterministic constant pieces like (3.27) should be contrasted with the
stochastic constant shot noise contributions due to discreteness effects. The
latter are extremely small for the Lyman alpha forest due to high column
densities, and hence can be neglected for the precision we are interested in.
Note that unlike the stochastic constant pieces, the deterministic constant
pieces do not require free parameters – their contributions are fully fixed in
terms of non-trivial bias operators that survive in the limit $K_{2}({\bf
q},-{\bf q})$, like $b_{2}$. In order to retain these terms, we choose the
FFTLog bias for $P_{22}$-type integrals to be $\nu=-1.7$

As far as the 13-type corrections are concerned, there is no single choice of
$\nu$ that would make them converge both in the IR and UV. Following the
discussion above, we use $\nu=-1.7$ that leads to the convergence in the UV,
but misses the spurious IR-sensitive contributions. Together with our choice
$\nu=-1.7$ for the $P_{22}$-type integrals, all spurious IR corrections are
consistently absent in our loop integrals.

The convergent UV contributions of one-loop integrals are proportional to the
displacement variance,

$\sigma_{v}^{2}=\frac{1}{3}\int_{\bf q}\frac{P_{\rm
lin}(q)}{q^{2}}\Big{|}_{\Lambda\text{CDM},z=2.8}=4.0~{}[h^{-1}\text{Mpc}]^{2}\,.$
(3.28)

With our choice of the FFTLog bias, this correction is computed over the
entire range of scales. The physical calculable part, however, should only
involve modes with $k\lesssim k_{\rm NL}$. This means that our choice of the
FFTLog bias will introduce an error of

$\frac{1}{6\pi^{2}}\int_{k_{\rm NL}}^{\infty}dq~{}q^{2}\frac{P_{\rm
lin}(q)}{q^{2}}\approx 4\cdot 10^{-3}~{}[h^{-1}\text{Mpc}]^{2}~{}\,.$ (3.29)

Even if this error is very small in absolute terms, we will find that the
Sherwood data is actually sensitive to these corrections. This is another
reason to keep the higher-derivative counterterms in the fit: they compensate
for errors in loop calculations.

##  4 Top-down derivation from the exponential map

The symmetry-based EFT expansions (3.15) and (3.17) are complete and general.
It is interesting to see directly how these expressions arise from the non-
linear transformation of the optical depth field. The optical depth itself is
assumed to be a selection-free biased tracer of the underlying matter density,
just like the galaxies. In that case all selection-dependent contributions in
the flux bias model should stem from the exponential map. The main objective
of this section is to show how this happens. Since we aim here at reproducing
the known results (3.17), the content of this Section would be somewhat
academic in nature. We will argue that it still important to the general
understanding of the perturbative Lyman alpha modeling.

The question how the exponential map generates selection-dependent bias
contributions was previously studied by Ref. [59]. Specifically, this work
showed how the exponential map and the SPT kernels generate certain terms in
the EFT expansion. In this chapter we build on results of [59] and study the
relationship between the exponential map and the EFT bias expansion in more
detail. The main novelty of our work is that we show explicitly that the
mapping and the requirement of renormalization generate a full set of the
line-of-sight dependent EFT operators (3.17). If we formally drop the
requirement of renormalization, the exponential map will produce only a
particular incomplete set of EFT operators with fixed Wilson coefficients, at
the tree level. We call the corresponding perturbative model “the tree-level
Tau-model.” Naively, this model is extremely predictive, as it requires only
one extra parameter (velocity gradient bias, or the mean optical depth) in
addition to the standard set of selection-free bias parameters. We show,
however, that this model receives infinitely large loop contributions, which
renders it inconsistent. Once the loop corrections are appropriately
renormalized, their finite pieces generate order one corrections to all bias
parameters, which will break the tree-level constraints imposed by the
exponential map. All in all, this means that even if we ignore counterterms in
the perturbative expansion of the Lyman alpha forest flux at the beginning,
the loop corrections will generate them and make the exponential map expansion
equivalent to the full EFT.

###  4.1 Renormalized exponential map expansion

We start with a top down derivation involving the renormalization of the flux
contrast field. It is a non linear transformation of the optical depth field
$\tau$, which we assume to be a selection-independent tracer of matter. The
exponential map reads

$F=e^{-\tau_{0}(1+\epsilon\delta_{\tau})}\,,$ (4.1)

where $\tau_{0}$ is the mean optical depth. In this spirit of EFT, this is a
“bare” parameter that has to be properly renormalized. We will address this
issue soon. $\epsilon$ is our order counting parameter that keeps track of the
order of our perturbative solution. It is to be set to 1 in the final result.
Note that $\tau$ is a function of the redshift space space coordinate [5],
i.e. the full RSD mapping (3.4) has already been applied to it. Taylor
expanding Eq. (3.4) up to cubic order we obtain,191919We use the notation
$[g]_{\bf k}=\int d^{3}x~{}e^{-i{\bf k}{\bf x}}g({\bf x})~{}$ for an arbitrary
function $g({\bf x})$.

$\begin{split}\delta^{(s)}_{\tau}({\bf k})=&\delta^{(r)}_{\tau}({\bf
k})-\frac{ik_{z}v_{\tau,z}({\bf k})}{{\cal H}}-\epsilon\frac{ik_{z}}{{\cal
H}}[\delta^{(r)}_{\tau}v_{\tau,z}]_{\bf
k}+\epsilon\frac{i^{2}k_{z}^{2}}{2{\cal H}^{2}}[v^{2}_{\tau,z}]_{\bf k}\\\
&+\epsilon^{2}\frac{i^{2}k_{z}^{2}}{2{\cal
H}^{2}}[v^{2}_{\tau,z}\delta^{(r)}_{\tau}]_{\bf
k}-\epsilon^{2}\frac{i^{3}k_{z}^{3}}{6{\cal H}^{3}}[v^{2}_{\tau,z}]_{\bf
k}\,,\end{split}$ (4.2)

where $\delta_{\tau}^{(r)}$ is the optical depth fluctuation field in real
space, and we have used our bookkeeping parameter $\epsilon$ to trace
perturbative order. In what follows we will drop the $(s)$ superscript and use
$\delta_{\tau}$ for the optical depth contrast in redshift space. As a next
step, we Taylor expand the expression (4.1),

$F=e^{-\tau_{0}}\left(1-\epsilon\tau_{0}\delta_{\tau}+\frac{(\epsilon\tau_{0}\delta_{\tau})^{2}}{2}-\frac{(\epsilon\tau_{0}\delta_{\tau})^{3}}{6}+O(\epsilon^{4})\right)\,,$
(4.3)

We are interested in one-loop order so the series goes up to
$O(\epsilon^{4})$. At this order, assuming $\langle\delta_{\tau}\rangle=0$,
the mean flux is given by

$\langle
F\rangle=e^{-\tau_{0}}\left(1+\epsilon^{2}\frac{\tau_{0}^{2}}{2}\langle\delta^{2}_{\tau}\rangle+...\right)\,.$
(4.4)

The mean flux is a function of the “bare coupling” $\tau_{0}$ that receives
corrections order by order in perturbation theory. These corrections are
unobservable since they are absorbed into $\tau_{0}$ in order to produce the
physically observable mean flux $\langle F\rangle$ after renormalization. The
situation here is similar to the renormalization of the “bare” mean number
density of galaxies in the context of the perturbative bias expansion [83].
The flux fluctuations are given by

$\begin{split}\delta_{F}=-\tau_{0}\epsilon\delta_{\tau}+\frac{\tau_{0}^{2}}{2}\epsilon^{2}(\delta^{2}_{\tau}-\langle\delta^{2}_{\tau}\rangle)-\frac{(\epsilon\tau_{0}\delta_{\tau})^{3}}{6}+{O}(\epsilon^{4})\,.\end{split}$
(4.5)

Setting $\epsilon=1$, and going into Fourier space we get,

$\begin{split}\delta_{F}({\bf k})=-\tau_{0}\delta_{\tau}({\bf
k})+\frac{\tau_{0}^{2}}{2}[\delta^{2}_{\tau}]_{\bf
k}-\frac{(\tau_{0})^{3}}{6}[\delta^{3}_{\tau}]_{\bf k}\,,\end{split}$ (4.6)

where $\delta_{\tau}$ is given in Eq. (4.2). Note that so far we have not used
the selection-independent bias model for $\delta^{(r)}_{\tau}$. Similar to Eq.
(4.2), the Taylor expansion above produces operators evaluated at the same
point in space that are highly UV sensitive. In the EFT approach, we need to
make sure that $\delta_{F}$ receives contributions only from modes up to
$k_{\rm NL}$ that are under perturbative control. That means all functions in
the r.h.s. of Eq. (4.3), e.g. $\delta_{\tau}$ as well as the composite local
operators [$\delta_{\tau}^{2}$] and $[\delta_{\tau}^{3}]$ need to be
renormalzied.

Let us start with the linear term $[\delta_{\tau}]$, see Eq. (4.2). First, one
can plug the perturbative bias model for the real space $\delta_{\tau}$,
which, according to our assumptions, does not have line-of-sight selection
effects, see Eq. (3.1). This model is closed under renormalization so we do
not need to add new counterterms. The velocity terms in Eq. (4.2) are
renormalized in the standard way, see Ref. [120]. Hence, this term also does
not require additional counterterms. Let us move to the operators
[$\delta_{\tau}^{2}$] and $[\delta_{\tau}^{3}]$. Their renormalized versions
must include all possible counterterms built out of the long-wavelength fields
and allowed by symmetries. In our case these are the equivalence principle and
$SO(2)$ rotations around the line-of-sight. The equivalence principle means
the counterterms can only be functions of $\partial_{i}\partial_{j}\Phi$ and
$\partial_{i}v_{j}$ and time derivatives along the flow. Up to cubic order,
these operators have already been introduced in Eq. (3.15). Thus, we write

$\begin{split}&[\delta^{2}_{\tau}]_{\bf k}\Big{|}_{\rm
ren}=[\delta^{2}_{\tau}]_{\bf
k}+b^{\prime}_{1}\delta+b^{\prime}_{\eta}\hat{z}^{i}\hat{z}^{j}\partial_{i}\partial_{j}\Phi+\sum
b^{\prime}_{\mathcal{O}}\mathcal{O}\,,\\\ &[\delta^{3}_{\tau}]_{\bf
k}\Big{|}_{\rm ren}=[\delta^{3}_{\tau}]_{\bf
k}+b^{{}^{\prime\prime}}_{1}\delta+b^{{}^{\prime\prime}}_{\eta}\hat{z}^{i}\hat{z}^{j}\partial_{i}\partial_{j}\Phi+\sum
b^{{}^{\prime\prime}}_{\mathcal{O}}\mathcal{O}\,,\end{split}$ (4.7)

where $\sum b^{{}^{\prime},^{\prime\prime}}_{\mathcal{O}}\mathcal{O}$ stand
for the quadratic and cubic operators in Eq. (3.15). This way we generate all
the terms in the expansion that we already had before in the “bottom-up” case.
The bias coefficients $b^{{}^{\prime}}_{\mathcal{O}}$ and
$b^{{}^{\prime\prime}}_{\mathcal{O}}$, along with the biases present in
$[\delta^{n}_{\tau}]_{\bf k}$, will sum up into the bias coefficients of Eq.
(3.15). Note that velocity-dependent contributions dictated by the equivalence
principle are contained inside the expressions $\delta^{2}_{\tau}$ and
$\delta^{3}_{\tau}$. They are unaffected by the counterterms, as the EFT
cannot generate terms that depend on the velocity itself, only on its
gradients. All in all, this way we have arrived at the same EFT expansion as
Eq. (3.15). We see that even if we start directly from the exponential map,
the smoothing and renormalization will generate all the necessary selection-
dependent bias operators.

To sum up, the renormalization procedure requires that coarse-grained
composite operators stemming from the exponential map contain all possible
line-of-sight dependent bias operators as counterterms,

$[\delta^{n}_{\tau}]_{\bf k}\Big{|}_{\rm ren}=[\delta^{n}_{\tau}]_{\bf
k}+\text{counterterms}\,.$ (4.8)

This is similar to usual renormalization conditions for local operators in
quantum field theory. If we were to ignore the counterterms, the perturbative
expansion generated by the exponential map would only involve simple
contractions of $\delta^{n}_{\tau}$ whose bias coefficients would be
completely fixed by the mapping. It is instructive to study the structure of
this expansion and contrast it with the full EFT. Since the counterterms are
typically considered along the loop corrections, with some abuse of language
we call such an expansion “the tree-level Tau model.” The goal of this
discussion is to give a pedagogical example of how the counterterms in Eq.
(4.8) will reappear from loops even if we ignore them to begin with.

###  4.2 Tree-level Tau model

Let us study in more detail the “unrenormalized” part of the exponential map
expansion. A similar derivation was recently done by [59]. We reproduce parts
of this derivation here and extend the results of [59] to non-linear
selection-dependent bias terms. It is interesting to re-derive them again in
the context of the standard perturbation theory (SPT) without renormalization
[88]. Our starting point would be the exponential map (4.1) that we will
Taylor expand and arrive at Eq. (4.5) for the fluctuations. The difference now
is that we assume that fluctuations of $\tau$ are perturbative in the matter
density, and they satisfy the standard bias expansion similar to
(3.1),202020In the SPT context this assumption is wrong because it is
precisely coarse-graining and renormalization that guarantee the adequacy of
the perturbative regime.

$\delta_{\tau}=\delta^{(1)}_{\tau}+\epsilon\delta^{(2)}_{\tau}+\epsilon^{2}\delta^{(3)}_{\tau}\,,$

where $\delta^{(n)}_{\tau}\sim Z_{n}[\delta^{(1)}]^{n}$. Plugging this into
(4.5) we get:

$\begin{split}\delta_{F}=&-\tau_{0}\epsilon\delta^{(1)}_{\tau}-\tau_{0}\epsilon^{2}\delta^{(2)}_{\tau}-\tau_{0}\epsilon^{3}\delta^{(3)}_{\tau}+\\\
&+\frac{\tau_{0}^{2}\epsilon^{2}}{2}((\delta^{(1)}_{\tau})^{2}+2\epsilon\delta^{(2)}_{\tau}\delta^{(1)}_{\tau}-\langle\delta^{2}_{\tau}\rangle)-\frac{(\epsilon\tau_{0}\delta^{(1)}_{\tau})^{3}}{6}+O(\epsilon^{4})\,.\end{split}$
(4.9)

Using the explicit expression for the standard redshift-space kernels (3.8),
the expansion above can be rewritten in Fourier space in the standard
perturbation theory form as

$\begin{split}\delta_{g}^{(s)}({\bf k})=&\tilde{K}_{1}({\bf
k})\delta^{(1)}_{\bf k}+\int_{{\bf q}_{1}}\int_{{\bf
q}_{2}}(2\pi)^{3}\delta_{D}({\bf k}-{\bf q}_{12})\tilde{K}_{2}({\bf
q}_{1},{\bf q}_{2})\delta^{(1)}_{{\bf q}_{1}}\delta^{(1)}_{{\bf q}_{2}}+\\\
&+\int_{{\bf q}_{1}}\int_{{\bf q}_{2}}\int_{{\bf
q}_{3}}(2\pi)^{3}\delta_{D}({\bf k}-{\bf q}_{123})\tilde{K}_{3}({\bf
q}_{1},{\bf q}_{2},{\bf q}_{3})\delta^{(1)}_{{\bf q}_{1}}\delta^{(1)}_{{\bf
q}_{2}}\delta^{(1)}_{{\bf q}_{3}}+...\end{split}$ (4.10)

where

$\begin{split}&\tilde{K}_{1}({\bf
k})=(-\tau_{0}b^{\tau}_{1})+(-\tau_{0})f\mu^{2}\,,\\\ &\tilde{K}_{2}({\bf
k}_{1},{\bf
k}_{2})=\frac{(-\tau_{0}b^{\tau}_{2})}{2}+(-\tau_{0}b^{\tau}_{\mathcal{G}_{2}})\left(\frac{({\bf
k}_{1}\cdot{\bf
k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}-1\right)+(-\tau_{0}b^{\tau}_{1})F_{2}({\bf
k}_{1},{\bf k}_{2})-\tau_{0}f\mu^{2}G_{2}({\bf k}_{1},{\bf k}_{2})\\\
&-\tau_{0}\frac{f\mu
k}{2}\left(\frac{\mu_{1}}{k_{1}}(b^{\tau}_{1}+f\mu_{2}^{2})+\frac{\mu_{2}}{k_{2}}(b^{\tau}_{1}+f\mu_{1}^{2})\right)+\frac{\tau_{0}^{2}}{2}((b^{\tau}_{1})^{2}+b^{\tau}_{1}f(\mu_{1}^{2}+\mu_{2}^{2})+f^{2}\mu_{1}^{2}\mu_{2}^{2})\,,\\\
&\tilde{K}_{3}({\bf k}_{1},{\bf k}_{2},{\bf k}_{3})=-\tau_{0}Z_{3}({\bf
q}_{1},{\bf q}_{2},{\bf q}_{3})+\tau_{0}^{2}[Z_{2}({\bf q}_{1},{\bf
q}_{2})Z_{1}({\bf q}_{3})]_{\rm symm.}-\frac{\tau_{0}^{3}}{3!}Z_{1}({\bf
q}_{1})Z_{1}({\bf q}_{2})Z_{1}({\bf q}_{3})\,.\end{split}$ (4.11)

Comparing the linear kernel $\tilde{K}_{1}$ with the EFT expression (3.19) we
identify

$\begin{split}&b_{1}=-\tau_{0}b^{\tau}_{1}\,,\\\
&b_{\eta}=\tau_{0}\,.\end{split}$ (4.12)

As anticipated, the exponential map generates the velocity gradient bias
already at the linear level [79]. Let us focus now on the quadratic kernel
$\tilde{K}_{2}$. Rewriting an expression in second line of the $\tilde{K}_{2}$
expansion above as

$\begin{split}&-\tau_{0}\frac{f\mu
k}{2}\left(\frac{\mu_{1}}{k_{1}}(b^{\tau}_{1}+f\mu_{2}^{2})+\frac{\mu_{2}}{k_{2}}(b^{\tau}_{1}+f\mu_{1}^{2})\right)+\frac{\tau_{0}^{2}}{2}(b^{\tau}_{1}f(\mu_{1}^{2}+\mu_{2}^{2})+f^{2}\mu_{1}^{2}\mu_{2}^{2})\\\
&=-f(\tau_{0}b^{\tau}_{1}-\tau_{0}^{2}b^{\tau}_{1})\frac{\mu_{2}^{2}+\mu_{1}^{2}}{2}+(-\tau_{0}+\frac{\tau_{0}^{2}}{2})f^{2}\mu_{1}^{2}\mu_{2}^{2}\\\
&+(-\tau_{0}b^{\tau}_{1})f\frac{\mu_{1}\mu_{2}}{2}\left(\frac{k_{2}}{k_{1}}+\frac{k_{1}}{k_{2}}\right)-\tau_{0}f^{2}\frac{\mu_{1}\mu_{2}}{2}\left(\frac{k_{2}}{k_{1}}\mu_{2}^{2}+\frac{k_{1}}{k_{2}}\mu_{1}^{2}\right)\,,\\\
\end{split}$ (4.13)

and matching the entire kernel with (3.19) we obtain:

$\begin{split}&b_{2}=-\tau_{0}b^{\tau}_{2}+\tau_{0}^{2}(b^{\tau}_{1})^{2}\,,\quad
b_{\mathcal{G}_{2}}=-\tau_{0}b^{\tau}_{\mathcal{G}_{2}}\,,\quad
b_{\delta\eta}=\tau_{0}b^{\tau}_{1}-\tau_{0}^{2}b^{\tau}_{1}\,,\\\
&b_{\eta^{2}}=-\tau_{0}+\frac{\tau_{0}^{2}}{2}\,,\quad~{}~{}~{}~{}~{}~{}~{}~{}b_{(KK)_{\parallel}}=0\,,\quad~{}~{}~{}~{}~{}b_{\Pi^{[2]}_{\parallel}}=0\,.\end{split}$
(4.14)

Assuming that the bias parameters of the optical depth field are order one
numbers, just like those of galaxies, $b_{1}^{\tau}\sim b_{2}^{\tau}\sim 1$,
we get an estimate $b_{1}\sim b_{2}\sim\tau_{0}\approx 0.3$ [53].

As far as the $\tilde{K}_{3}$ are concerned, we see that in addition to the
standard term $\sim Z_{3}$, there are new contributions proportional to
$Z_{1}^{3}$ and $Z_{2}Z_{1}$ kernels, which are generated by the exponential
map. However, one can easily check that only the $Z_{2}Z_{1}$ term contributes
non-trivially to the one-loop power spectrum. More specifically, only the term
containing $G_{2}$ has a finite contribution that consistency reproduces the
cubic order terms that must stem from the non-trivial $\delta\eta$, $\eta^{2}$
and $\eta$ operators. The rest of the additional contributions are
proportional to the UV-diverging mass variance

$\sigma^{2}_{\Lambda}=\int_{0}^{\Lambda}\frac{dq}{2\pi^{2}}~{}q^{2}P_{\text{lin}}(q)\,,$
(4.15)

where $\Lambda$ is a UV cutoff. We will discuss the renormalization of these
13-type terms shortly. At the kinematic configurations relevant for the one-
loop power spectrum we thus have

$\begin{split}&\int_{{\bf q}}\tilde{K}_{3}({\bf k},{\bf q},-{\bf
q})P_{\text{lin}}(q)=-\tau_{0}b^{\tau}_{1}\int_{{\bf q}}F_{3}({\bf k},{\bf
q},-{\bf q})P_{\text{lin}}(q)-f\tau_{0}\mu^{2}\int_{{\bf q}}G_{3}({\bf
q},-{\bf q},{\bf k})P_{\text{lin}}(q)\\\ &\>+\int_{{\bf
q}}\left[1-\left(\hat{\bf k}\cdot\hat{\bf
q}\right)^{2}\right]P_{\text{lin}}(q)\\\
&\>\times\Bigg{\\{}\frac{4}{21}(5b_{\mathcal{G}_{2}}+2b_{\Gamma_{3}})\left[\left(\frac{({\bf
k}-{\bf q})\cdot{\bf q}}{|{\bf k}-{\bf
q}|q}\right)^{2}-1\right]-\frac{2}{21}f\left(-\tau_{0}+\frac{\tau_{0}^{2}}{2}\right)\left[\frac{3(k_{\parallel}-q_{\parallel})^{2}}{|{\bf
k}-{\bf q}|^{2}}+\frac{5q_{\parallel}^{2}}{q^{2}}\right]\\\
&\quad\>+\frac{4}{7}f^{2}\left(\tau_{0}b_{1}^{\tau}-b_{1}^{\tau}\tau_{0}^{2}\right)\frac{q_{\parallel}^{2}}{q^{2}}\frac{(k_{\parallel}-q_{\parallel})^{2}}{|{\bf
k}-{\bf q}|^{2}}\\\
&\quad\>+\frac{2}{21}f(-\tau_{0}b^{\tau}_{1})\left[5\frac{q_{\parallel}(k_{\parallel}-q_{\parallel})}{q^{2}}+3\frac{q_{\parallel}(k_{\parallel}-q_{\parallel})}{|{\bf
k}-{\bf
q}|^{2}}\right]-\frac{2}{7}f^{2}\tau_{0}\frac{q_{\parallel}(k_{\parallel}-q_{\parallel})}{q^{2}|{\bf
k}-{\bf
q}|^{2}}\left[(k_{\parallel}-q_{\parallel})^{2}+q_{\parallel}^{2}\right]\Bigg{\\}}\,.\end{split}$
(4.16)

Comparing this with the full EFT expression, we observe that the Tau model
predicts that all EFT cubic biases are zero except
$b_{\Gamma_{3}}=-\tau_{0}b_{\Gamma_{3}}^{\tau}$.

We call Eq. (4.11) the tree-level Tau model. This model involves only one
extra parameter, $\tau_{0}$. It also suggests that most of the one-loop
selection-dependent EFT operators vanish,

$\begin{split}b_{(KK)_{\parallel}}=b_{\Pi^{[2]}_{\parallel}}=b_{\Pi^{[3]}_{\parallel}}=b_{\delta\Pi^{[2]}_{\parallel}}=b_{\eta\Pi^{[2]}_{\parallel}}=b_{(K\Pi^{[2]})_{\parallel}}=0\,,\end{split}$
(4.17)

and that all non-trivial quadratic selection-dependent biases in the EFT
expansion can be expressed through linear biases as

$\begin{split}&b_{\delta\eta}=-b_{1}(1-b_{\eta})\,,\quad
b_{\eta^{2}}=-b_{\eta}+\frac{b_{\eta}^{2}}{2}\,.\end{split}$ (4.18)

Our discussion so far has not included the loop corrections, which we discuss
in detail now.

###  4.3 The effect of loop corrections

The inconsistencies of the tree-level Tau model become apparent once we take
loops into account. In particular, they will generate formally infinite
corrections to the tree-level bias parameters of the linear density and
projected velocity gradient fields. These biases get renormalized by the
following $P_{13}$-type contribution:

$\begin{split}&\delta
P_{13}=(b_{1}-b_{\eta}f\mu^{2})P_{\text{lin}}(k)\int_{\bf
q}P_{\text{lin}}(q)\left[(-\tau_{0})^{3}(b^{\tau}_{1}+f\mu^{2})Z^{2}_{1}({\bf
q})+4\tau_{0}^{2}Z_{2}({\bf k},{\bf q})Z_{1}({\bf q})\right]\,,\end{split}$
(4.19)

which produces the leading order in $k$ UV contributions. The above integral
can be done exactly, except for the term that contains $\sim G_{2}$ in the
rightmost term. The exact part of this expression is simply proportional to
the short-scale mass variance (4.15), see Appendix B. The non-trivial part of
the $Z_{2}Z_{1}$ contribution has the same UV-limit, proportional to
$\sigma^{2}_{\Lambda}$. Thus, for the purposes of complete renormalization of
the UV behavior, we need to add this part. This way we obtain

$\delta P_{13}\Big{|}_{\rm UV}=2(b_{1}-b_{\eta}f\mu^{2})P_{\rm
lin}(k)(A_{0}-A_{2}f\mu^{2})\sigma_{\Lambda}^{2}\,,$ (4.20)

where the oder-one constants $A_{0}$ and $A_{2}$ can be found in Appendix B.
Note that $b_{1}$ also receives contributions from the usual bias operators,
but they can be treated in the standard way [117], and we will ignore these
terms in our discussion for clarity.

It is convenient now to recover our bookkeeping perturbation theory parameter
$\epsilon$. Adding the linear result, we obtain the following expression for
the tree-level plus the UV part of the one-loop power spectrum in the Tau
model,

$\begin{split}P_{\rm tree}(k)+P_{\rm 1-loop}(k)\Big{|}_{\rm
UV}=&(b_{1}-b_{\eta}f\mu^{2})^{2}P_{\rm
lin}(k)+2(b_{1}-b_{\eta}f\mu^{2})P_{\rm
lin}(A_{0}-A_{2}f\mu^{2})\epsilon^{2}\sigma_{\Lambda}^{2}\\\ =&P_{\rm
lin}(k)(b_{1}+A_{0}\sigma_{\Lambda}^{2}\epsilon^{2}-(\tau_{0}+A_{2}\epsilon^{2}\sigma_{\Lambda}^{2})f\mu^{2})^{2}+O(\epsilon^{4})\,,\end{split}$
(4.21)

where $O(\epsilon^{4})$ are formally two-loop corrections that can be ignored
at the precision level we are interested in. We have seen that the leading
order loop corrections in the Tau model produce infinite contributions to the
linear bias parameters, including the velocity gradient bias. This merely
implies that the SPT expansion for the Lyman alpha forest is inconsistent, and
one has to use the full EFT model that includes all operators allowed by
symmetries.

In the EFT approach one must treat $b_{1}$ and $b_{\eta}$ as formally infinite
“bare” parameters that absorb the leading order UV behavior from the loops.
Then their finite, physically observable, renormalized values of $b_{1}$ and
$b_{\eta}$ are given by the sum of the “bare” parameters and the UV-sensitive
loop contribution,

$b_{\eta}^{\rm ren}=\tau_{0}+A_{2}\sigma_{\Lambda}^{2}\,,\quad b^{\rm
ren}_{1}=-\tau_{0}b^{\tau}_{1}+A_{0}\sigma_{\Lambda}^{2}\epsilon^{2}\,.$
(4.22)

Let us estimate the finite part of $b_{\eta}^{\rm ren}$ now. To that end we
can use $\Lambda\sim k_{\rm NL}$, so that $\sigma^{2}_{\rm NL}\sim 1$ by
definition. In the absence of fine tuning, it should be of the order of

$b_{\eta}^{\rm ren}\sim A_{2}\sigma^{2}_{\rm NL}\sim 1\,,\quad b_{1}^{\rm
ren}\sim A_{0}\sigma^{2}_{\rm NL}\sim 1\,.$ (4.23)

Hence, even if we take the tree-level Tau model at face value, the loops would
generate an order-one breaking212121From the structure of the renormalization
condition Eq. (4.22) one can expect certain suppressions of the loop
corrections due to the smallness of the mean optical depth $\tau_{0}\approx
0.3$. We will not entertain this possibility in what follows. of the
relationship between the EFT operators following from the tree-level Tau
model. Specifically, Eq. (4.23) implies that $b_{\eta}^{\rm ren}$ and
$b_{1}^{\rm ren}$ are order one different from $-\tau_{0}b_{1}^{\tau}$ and
$\tau_{0}$, and hence the tree-level exponential mapping constraints (4.18) do
not hold once the loops are taken into account. This also suggests that the
other constraint (4.17), i.e. the vanishing of specific selection-dependent
EFT operators, would not hold true as well because of the loops. This is quite
natural to expect, as there is no symmetry that would guarantee the vanishing
of these operators in the presence of loops. In the EFT jargon one says that
operators like $\Pi^{[2]}_{\parallel}$ are generated by loops even if they are
not present in the exponential mapping at the tree-level.

##  5 Comparison with simulations

###  5.1 Fitting details

In this Section we compare the one-loop EFT power spectrum model for the Lyman
alpha flux power spectrum against simulation data. To that end we use the
state-of-the-art suite of Sherwood simulations [52]. These are large-scale,
high resolution hydrodynamic simulations that reproduce the evolution of
intergalactic medium with a large number of particles (up to $1.7\times
10^{10}$). The fiducial cosmology of these simulations is a flat $\Lambda$CDM
with massless neutrinos and $\Omega_{m}=0.308$, $\Omega_{b}=0.0482$,
$\sigma_{8}=0.829$, $n_{s}=0.961$, $h=0.678$. The details of the Sherwood
simulations can be found in Refs. [52, 60].

In this work, we fit the 3D Lyman alpha power spectrum at $z=2.8$, extracted
from the simulation box L160_N2048, which followed 20483 gas and cold dark
matter particles in a cub with side length L=160$~{}h^{-1}\text{Mpc}$. The
simulations use a homogeneous ionizing background model, and assume that the
gas is in ionization equilibrium and optically thin [52]. The final 3D power
spectrum measurements are publicly available222222
https://github.com/andreufont/sherwood_p3d. They are presented as functions of
the Fourier mode wavenumber $k$ and the cosine $\mu$ of the angle between the
corresponding Fourier vector and the line of sight. In practice, the $\mu$
space is sampled by 16 bins, while the $k$ space is log-sampled in the range
[$k_{F},k_{\rm Ny}$], where $k_{F}=0.039~{}h\text{Mpc}^{-1}$ is the
fundamental mode, and $k_{\rm Ny}=20~{}h\text{Mpc}^{-1}$ is the Nyquist
frequency.

Following [60], we fit the power spectrum by Monte-Carlo Markov Chain sampling
the following pseudo-$\chi^{2}$ function,

$\chi^{2}=\sum_{i}\frac{[P^{\rm data}_{i}-P^{\rm
model}(k_{i},\mu_{i})]^{2}}{2(P^{\rm data}_{i})^{2}/N_{i}}\,,$ (5.1)

where $P^{\rm data}_{i}$ are measurements power spectrum bins, and $N_{i}$ is
the number of modes in the bin. Note that in Eq. (5.1) we have explicitly
assumed a Gaussian diagonal covariance for the power spectrum. While this
assumption is definitely correct on large scales, it has not been
systematically validated on small scales, relevant for our analysis. Studies
of the galaxy power spectrum covariance in Refs. [136, 137, 138] showed that
the Gaussian approximation is very accurate on mildly non-linear scales
because in this regime the effective covariance is dominated by the
theoretical error due to the marginalization over nuisance parameters [139,
140, 141]. In what follows we assume that the same argument is true for the
Lyman alpha power spectrum, and we proceed with the Gaussian diagonal
covariance. We warn, however, that the interpretation of our results is,
strictly speaking, contingent upon our covariance matrix assumptions.

Note that unlike [60], we do no add the ad hoc noise floor correction to the
covariance, as our goal here is to fit the data given its actual precision.
This allows us to illustrate the main advantage of the EFT over
phenomenological approaches: its high precision on the scales where the EFT is
applicable.

The model vector $P^{\rm model}(k_{i},\mu_{i})$ is a one-loop EFT model that
includes the necessary stochastic contributions and counterterms,

$\begin{split}P^{\rm
1-loop}(k,\mu)&=(b_{1}-b_{\eta}f\mu^{2})^{2}P_{\text{lin}}(k)\\\ &+2\int_{\bf
q}K_{2}^{2}({\bf q},{\bf k}-{\bf q})P_{\text{lin}}(|{\bf k}-{\bf
q}|)P_{\text{lin}}(q)\\\ &+6K_{1}({\bf k})P_{\text{lin}}(k)\int_{\bf
q}K_{3}({\bf k},-{\bf q},{\bf q})P_{\text{lin}}(q)\\\
&-2(c_{0}+c_{2}\mu^{2}+c_{4}\mu^{4})k^{2}P_{\text{lin}}(k)\,.\end{split}$
(5.2)

The $k^{2}P_{\rm lin}$ corrections are added here as proxy for the two-loop
contributions. For the purposes of the data analysis, we switched to a
simplified version of the higher derivative contribution that is motivated by
the EFT-based galaxy power spectrum analyses, e.g. [70]. We have checked that
adding an extra parameter $\propto k^{2}\mu^{6}P_{\text{lin}}(k)$ has no
impact on the fit. In addition, we found that the effect of the stochastic
counterterms (see Eq. (3.22)) is negligible for our data, and hence we set
them to zero. This is consistent with our theoretical arguments.

At face value, the one-loop EFT model depends on 16 free parameters: 2 linear
biases, 11 non-linear biases, and 3 higher derivative operators. We found,
however, that the quality of our data does not allow us to break the known
degeneracy between $b_{\mathcal{G}_{2}}$ and $b_{\Gamma_{3}}$ that exists at
the power spectrum level [67, 142, 72, 78]. We thus set $b_{\Gamma_{3}}=0$
following [67, 143]. For the remaining parameters, we use the following
priors:

$\begin{split}&b_{1}\in[-2,2]\,,\quad b_{\eta}\in[-2,2]\,,\quad
b_{2}\sim\mathcal{N}(0,2^{2})\,,\quad
b_{\mathcal{G}_{2}}\sim\mathcal{N}(0,2^{2})\,,\\\
&b_{(KK)_{\parallel}}\sim\mathcal{N}(0,2^{2})\,,\quad
b_{\Pi^{[2]}_{\parallel}}\sim\mathcal{N}(0,2^{2})\,,\quad
b_{\Pi^{[3]}_{\parallel}}\sim\mathcal{N}(0,2^{2})\,,\quad\\\
&b_{\delta\eta}\sim\mathcal{N}(0,2^{2})\,,\quad
b_{\eta^{2}}\sim\mathcal{N}(0,2^{2})\,,\quad\frac{c_{0,2,4}}{[h^{-1}\text{Mpc}]^{2}}\sim\mathcal{N}(0,0.1^{2})\,,\\\
&b_{(K\Pi^{[2]})_{\parallel}}\sim\mathcal{N}(0,2^{2})\,,\quad
b_{\delta\Pi^{[2]}_{\parallel}}\sim\mathcal{N}(0,2^{2})\,,\quad
b_{\eta\Pi^{[2]}_{\parallel}}\sim\mathcal{N}(0,2^{2})\,,\end{split}$ (5.3)

where $\mathcal{N}(\mu,\sigma^{2})$ stands for a Gaussian distribution with
mean $\mu$ and r.m.s. $\sigma$. Let us discuss in detail the motivation for
this choice of priors. First, the priors for $b_{1}$ and $b_{\eta}$ are chosen
to be flat and uninformative. Second, for all the quadratic and cubic
selection-dependent bias parameters we choose large uninformative Gaussian
priors with zero mean and variance square roots of 2. In general, the EFT
naturalness arguments dictate that they should be $O(1)$ numbers, i.e. we
expect them to roughly be in the range $1\div 10$. We note, however, that some
additional suppression may come from the mapping, i.e. the flux bias
parameters are $O(1)$ bias parameters of the $\tau$ field that have got
multiplied by $\tau_{0}\simeq 0.3$. Given this reason, we reduced the r.m.s.
of our prior to 2 instead, e.g. a conservative choice of 10 that could be made
without any a priori knowledge.

As far as the counterterms are concerned, as we discussed earlier, we include
them mainly in order to absorb some two-loop corrections that are not
explicitly computed. Indeed, we have found that they actually improve the fit,
which justifies their presence a posteriori. The typical size of these
corrections expected in the EFT is $k_{\rm NL}^{-2}\sim
0.05~{}[h^{-1}\text{Mpc}]^{2}$. We additionally multiply this by 2 in order to
be conservative. Indeed, one might expect the velocity field to be more
nonlinear than the density one, and hence the $k^{2}$-contributions may be
enhanced, as it is the case for galaxies [67, 144]. Jumping ahead, let us note
that we will actually find the higher derivative corrections whose values are
consistent with $k_{\rm NL}^{-2}$, with errorbars that are much tighter than
the priors.

Only the 15 nuisance parameters (5.3) are assumed to vary in our MCMC chains.
We leave the fitting of cosmological parameters for future work. We sample the
likelihood analytically marginalized over all parameters that enter the
likelihood quadratically, i.e.

$\\{P_{\rm
shot},c_{0},c_{2},c_{4},b_{\Pi^{[3]}_{\parallel}},b_{(K\Pi^{[2]})_{\parallel}},b_{\delta\Pi^{[2]}_{\parallel}},b_{\eta\Pi^{[2]}_{\parallel}}\\}\,,$

as in Refs. [138, 141]. The approximate posterior distributions and best-fit
values for these parameters are later recovered from the MCMC chains for the
marginalized likelihood. The rest of the EFT parameters in Eq. (5.3) are
explicitly varied in our MCMC chains.

The MCMC chains are run with the Montepython sampler [145, 146] and analyzed
with getdist [147].

###  5.2 Results

Figure 2: Triangle plot and marginalized projections for bias parameters of the EFT model for the Lyman alpha flux power spectrum of the Sherwood simulations at $z=2.8$. We show results for four choices of $k_{\rm max}$: 2, 3, 4, and 5 $h\text{Mpc}^{-1}$ (gray, green, blue, red, respectively). Param | best-fit | mean$\pm\sigma$ | 95% lower | 95% upper
---|---|---|---|---
$b_{1}$ | $-0.220$ | $-0.217_{-0.01}^{+0.01}$ | $-0.235$ | $-0.198$
$b_{\eta}$ | $0.322$ | $0.324_{-0.024}^{+0.024}$ | $0.276$ | $0.372$
$b_{2}$ | $-0.29$ | $-0.24_{-0.28}^{+0.24}$ | $-0.72$ | $0.27$
$b_{\mathcal{G}_{2}}$ | $-0.13$ | $-0.09_{-0.21}^{+0.17}$ | $-0.45$ | $0.30$
$b_{\eta^{2}}$ | $0.003$ | $0.075_{-0.2}^{+0.1}$ | $-0.28$ | $0.47$
$b_{\delta\eta}$ | $-0.85$ | $-0.74_{-0.27}^{+0.16}$ | $-1.13$ | $-0.27$
$b_{(KK)_{\parallel}}$ | $0.38$ | $0.49_{-0.62}^{+0.6}$ | $-0.69$ | $1.67$
$b_{\Pi^{[2]}_{\parallel}}$ | $-0.38$ | $-0.41_{-0.16}^{+0.17}$ | $-0.73$ | $-0.07$
$10^{3}c_{0}/[h^{-1}\text{Mpc}]^{2}$ | $3.59$ | $3.59_{-1.1}^{+1.1}$ | $1.47$ | $5.71$
$10^{3}c_{2}/[h^{-1}\text{Mpc}]^{2}$ | $2.64$ | $2.64_{-0.62}^{+0.62}$ | $1.39$ | $3.88$
$10^{3}c_{4}/[h^{-1}\text{Mpc}]^{2}$ | $-0.83$ | $-0.83_{-0.75}^{+0.75}$ | $-2.32$ | $0.67$
$b_{\Pi^{[3]}_{\parallel}}$ | $0.82$ | $0.82_{-0.05}^{+0.05}$ | $0.72$ | $0.92$
$b_{\delta\Pi^{[2]}_{\parallel}}$ | $-0.81$ | $-0.81_{-0.05}^{+0.05}$ | $-0.71$ | $0.91$
$b_{(K\Pi^{[2]})_{\parallel}}$ | $-1.52$ | $-1.52_{-0.07}^{+0.07}$ | $-1.68$ | $-1.38$
$b_{\eta\Pi^{[2]}_{\parallel}}$ | $1.95$ | $1.95_{-0.1}^{+0.1}$ | $1.76$ | $2.14$
Table 1: One-dimensional marginalized constraints on nuisance parameters of
the one-loop EFT model fit to the Lyman alpha forest flux power spectrum of
the Sherwood simulation at $z=2.8$. The upper group of parameters were
directly sampled in our MCMC chains. The lower parameters are analytically
marginalized over in the likelihood. Their posteriors are recovered from the
chains a posteriori.

In Fig. 2 we present 1d and 2d marginalized posterior distributions for the
bias parameters of the one-loop EFT model that were directly sampled in our
chains. We show results for four choices of $k_{\rm max}=2,3,4$ and $5$
$h\text{Mpc}^{-1}$. We see that the posteriors for $k_{\rm
max}=2~{}h\text{Mpc}^{-1}$ and $k_{\rm max}=3~{}h\text{Mpc}^{-1}$ are fully
consistent with each other. The contour for the $k_{\rm
max}=4~{}h\text{Mpc}^{-1}$ case is, however, significantly shifted w.r.t. the
$2$ and $3$ $h\text{Mpc}^{-1}$ posteriors. Although the $k_{\rm
max}=5~{}h\text{Mpc}^{-1}$ contour shifts back inside the $k_{\rm
max}=2~{}h\text{Mpc}^{-1}$ posterior, it is still inconsistent with the
$k_{\rm max}=3~{}h\text{Mpc}^{-1}$ contour. These dramatic shifts and
reduction of the posterior volume suggest that the fit is biased for $k_{\rm
max}=4~{}h\text{Mpc}^{-1}$ and $k_{\rm max}=5~{}h\text{Mpc}^{-1}$ due to two-
loop corrections that are not included in our model. As a frequentest
confirmation of our scale cuts, we also see a significant deterioration of the
best-fit $\chi^{2}$ statistic for $k_{\rm max}>3~{}h\text{Mpc}^{-1}$. This
suggests to choose the $k_{\rm max}=3~{}h\text{Mpc}^{-1}$ case as a baseline.
Note that a similar behavior of the posteriors when sliding $k_{\rm max}$ was
previously observed in the case of redshift-space galaxy mocks [143], where
the same criterion of the “stability” and consistency of posteriors w.r.t.
variations of $k_{\rm max}$ was applied.

The best-fit values, 68% and 95% confidence limits for the fit parameters at
$k_{\rm max}=3~{}h\text{Mpc}^{-1}$ are presented in Table 1. First, we see
that the linear bias parameters $b_{1}$ and $b_{\eta}$ are measured to $\sim
5\%$ precision. The second important observation is that many Lyman alpha bias
parameters are detected at high significance:
$b_{\delta\eta},b_{\Pi^{[2]}_{\parallel}}$, $b_{\Pi^{[3]}_{\parallel}}$,
$b_{(K\Pi^{[2]})_{\parallel}}$,$b_{\delta\Pi^{[2]}_{\parallel}}$,
$b_{\eta\Pi^{[2]}_{\parallel}}$. Detections of these operators in the data
clearly indicate an inconsistency of the tree-level Tau model, and the power
of the EFT approach. The best-fit value of $b_{\delta\eta}$ is 4$\sigma$ away
from the prediction of the tree-level Tau model
$\tilde{b}_{\delta\eta}=-b_{1}(1-b_{\eta})\approx 0.15$, evaluated with Eq.
(4.18) using the best-fit value of $b_{\eta}$. The apparent tension between
the tree-level Tau model and the data is an “experimental” proof of the
inconsistency of the tree-level calculations for the Lyman alpha forest. As
far as higher derivative counterterms are concerned, we have a significant
detection of $c_{0}$, and $c_{2}$.

Figure 3: Best-fit EFT models against the simulated power spectra (left
panel), and the residuals between the model and the data (right panel).

Figure 4: Left panel: the magnitude of one-loop corrections relative to the
linear theory answer. Right panel: the fit for the transverse modes.

The best-fit models for four angular bins are shown against the data in Fig.
3. The nominal $\chi^{2}$ statistics across the 161 data points is 188, which
indicates a satisfactory fit for 15 free parameters. It is important to keep
in mind, however, that our EFT parameters are quite degenerate and hence the
counting of degrees of freedom is not straightforward in our case. In
addition, our Gaussian covariance assumptions are likely not very accurate at
the scales of interest. This is important to keep in mind when interpreting
the $\chi^{2}$ statistics.

An important consistency check is that the one-loop corrections are smaller
than the tree level result. To verify that this is the case, in Fig. 4 we plot
the one-loop contributions divided by the tree level model (3.13). We see that
the magnitude of the perturbative corrections does not exceed the linear
theory power spectrum, which is a necessary condition of the applicability of
perturbation theory on the scales of interest. The size of the one-loop
corrections, however, is almost as large as the tree-level result for $k\simeq
3~{}h\text{Mpc}^{-1}$, which suggests that higher loop corrections may not be
negligible. Given that no biases are observed at the level of the parameter
estimation, it is likely that the effect of the two-loop corrections was
partly absorbed by the one-loop nuisance parameters and counterterms. A
separate dedicated analysis is needed in order to quantify the effect of the
two loop corrections more accurately. We note however, that even if two-loop
corrections are not negligible at $k=3~{}h\text{Mpc}^{-1}$ in the strict
sense, the EFT can still be applied on these scales as a phenomenological
model that is capable of fitting the data with sub-percent accuracy.

Finally, the right panel of Fig. 4 illustrates that for the transverse
($\mu\approx 0$) modes the EFT provides a good fit even for $k_{\rm
max}=4~{}h\text{Mpc}^{-1}$. This suggests that the breakdown of the one-loop
model happens because of the velocity field that dominates the signal along
the line of sight. Indeed, in general, the velocity field is more nonlinear
than the density field for dark matter and galaxies [112, 148, 141, 144, 73].
Fig. 4 shows, however, that the modes along the line of sight and transverse
to the line of sight are roughly equally nonlinear for the Lyman alpha forest
at $z=2.8$. Overall, the quality of the $z=2.8$ data does not allow us to make
a definitive conclusion about the role of the velocity non-linearities. More
insight may be obtained from the analysis of the $z=3.2$ data, which we
discuss below.

We present the analysis of the Lyman alpha forest power spectrum at $z=3.2$ in
Appendix C. At this redshift, the non-linear wavenumber becomes significantly
large, $k_{\rm NL}=10~{}h\text{Mpc}^{-1}$, suggesting that the EFT expansion
should converge better and one can push the one-loop fit to smaller scales
where the errorbars are smaller. We found that this is indeed the case. The
EFT provides a good fit to the data up to $k_{\rm max}=5~{}h\text{Mpc}^{-1}$.
The optimal values of the Lyman alpha bias parameters are larger at $z=3.2$,
and they are detected with higher significance, which is consistent with the
expectations that tracers become more biased with redshift for uniform
selection criteria. In particular, we have a significant detection of the RSD
counterterms, whose impact increases towards $\mu=1$, e.g. the $c_{4}$
counterterm is detected at almost $20\sigma$. This supports the above argument
that the velocity field is primarily responsible for the breakdown of
perturbation theory at small scales. On the practical side, this suggests that
a more optimal strategy for parameter inference would be to use a
$\mu$-dependent $k_{\rm max}$ cutoff. We will explore this option in future
work.

##  6 One dimensional power spectrum

We have seen that the EFT can accurately describe the full 3D power spectrum
of the Lyman alpha forest. Let us show now how the EFT technique can be used
to calculate the one-dimensional (1D) power spectrum of the Lyman alpha
forests. Imagine that the flux fluctuations have an underlying power spectrum
$P_{\rm 3D}(k,k_{\parallel})$ in 3D, where $k_{\parallel}$ is the wavevector
projection along the line-of-sight. Then the power spectrum of these
fluctuations seen only along the line-of-sight is given by [149]

$P_{\rm{1D}}(k_{\parallel})=\frac{1}{2\pi}\int_{k_{\parallel}}^{\infty}dk~{}kP_{\rm
3D}(k,k_{\parallel})\,.$ (6.1)

A conceptual difficulty of this observable is that the 1D power spectrum is an
integral over the 3D power spectrum that involves UV modes that are not under
robust analytic control. This is precisely the issue that the EFT is aimed to
resolve. In fact, Eq. (6.1) has the form of a loop integral, so we can treat
this UV sensitivity along the lines of the one-loop EFT renormalization for
the matter power spectrum [85]. To that end we assume that the integral in Eq.
(6.1) is cut off at a scale $\Lambda\lesssim k_{\rm NL}$. A first relevant
observation is that in perturbation theory $P_{\rm 3D}(k,k_{\parallel})$ has a
simple polynomial dependence on $k_{\parallel}=k\cdot\mu$. For instance, at
one loop order we have

$P_{\rm
3D}(k,k_{\parallel})=\sum_{n=0}^{4}\left(\frac{k_{\parallel}}{k}\right)^{2n}P_{2n}(k)\,.$
(6.2)

A calculation of $P_{\rm 3D}(k,k_{\parallel})$ at a given order in
perturbation theory is valid only up to a certain scale that we call $k_{\rm
trust}$. Then Eq. (6.1) can be equivalently recast as

$P_{\rm{1D}}(k_{\parallel})=\frac{1}{2\pi}\int_{k_{\parallel}}^{k_{\rm
trust}}dk~{}kP_{\rm 3D}(k,k_{\parallel})+\frac{1}{2\pi}\int_{k_{\rm
trust}}^{\Lambda}dk~{}kP_{\rm 3D}(k,k_{\parallel})\,.$ (6.3)

The rightmost integral is given by simple polynomials of $k_{\parallel}$ times
some function of the cutoff,

$\begin{split}&\frac{1}{2\pi}\int_{k_{\rm trust}}^{\Lambda}dk~{}kP_{\rm
3D}(k,k_{\parallel})=\sum_{n=0}^{4}k^{2n}_{\parallel}c_{n}^{\Lambda}\,,\quad
c_{n}^{\Lambda}\equiv\frac{1}{2\pi}\int_{k_{\rm
trust}}^{\Lambda}dk~{}k^{1-2n}P_{2n}(k)\,.\end{split}$ (6.4)

These UV-sensitive contributions are naturally absorbed into the EFT
stochastic counterterms, which have the scale dependence that exactly match
Eq. (6.4). Indeed, Eq. (3.22) implies that

$P^{\rm stoch}_{\rm 3D}(k,k_{\parallel})=P_{\rm shot}+a_{0}\frac{k^{2}}{k_{\rm
NL}^{2}}+a_{2}\frac{k^{2}_{\parallel}}{k_{\rm NL}^{2}}+\cdots\,,$ (6.5)

where “$\cdots$” stand for terms higher order in $k/k_{\rm NL}$. $P_{\rm
shot},a_{0},a_{2}$ above are $\Lambda$-dependent Wilson coefficients that
cancel the UV sensitivity from the integrals (6.4). Plugging Eq. (6.5) into
Eq. (6.3) we get:

$\begin{split}&\frac{1}{2\pi}\int_{k_{\parallel}}^{k_{\rm trust}}dk~{}kP^{\rm
stoch}_{\rm 3D}(k,k_{\parallel})={P_{\rm shot}}\frac{k_{\rm
trust}^{2}-k_{\parallel}^{2}}{4\pi}+{a_{0}}\frac{k_{\rm
trust}^{4}-k_{\parallel}^{4}}{8\pi k_{\rm
NL}^{2}}+{a_{2}}k_{\parallel}^{2}\frac{k_{\rm
trust}^{2}-k_{\parallel}^{2}}{4\pi k_{\rm NL}^{2}}\\\ &=\frac{2P_{\rm
shot}k_{\rm trust}^{2}+a_{0}k_{\rm trust}^{4}/k_{\rm
NL}^{2}}{8\pi}+k_{\parallel}^{2}\left(a_{2}\frac{k_{\rm trust}^{2}}{4\pi
k_{\rm NL}^{2}}-\frac{P_{\rm shot}}{4\pi}\right)+\frac{k_{\parallel}^{4}}{4\pi
k_{\rm NL}^{2}}\left(-\frac{a_{0}}{2}-{a_{2}}\right)\,,\\\ &\equiv
C^{\prime}_{0}+k_{\parallel}^{2}C^{\prime}_{1}+k_{\parallel}^{4}C^{\prime}_{2}\,.\end{split}$
(6.6)

Adding this to (6.4) we obtain a finite cutoff-independent expression,

$\begin{split}&\frac{1}{2\pi}\int_{k_{\parallel}}^{k_{\rm trust}}dk~{}kP^{\rm
shoch}_{\rm 3D}+\frac{1}{2\pi}\int_{k_{\rm trust}}^{\Lambda}dk~{}kP_{\rm
3D}\\\
&=\left(C_{0}^{\Lambda}+C^{\prime}_{0}\right)+k_{\parallel}^{2}\left(C_{1}^{\Lambda}+C^{\prime}_{1}\right)+k_{\parallel}^{4}\left(C_{2}^{\Lambda}+C^{\prime}_{2}\right)+\cdots\\\
&=C_{0}^{\rm finite}+C_{1}^{\rm finite}k_{\parallel}^{2}+C_{2}^{\rm
finite}k_{\parallel}^{4}+\cdots\,,\end{split}$ (6.7)

where $c_{n}^{\rm finite}$ are finite parts of the EFT coefficients that have
to be matched to data. To sum up, in the EFT the 1D power spectrum calculation
amounts to doing the integral over the deterministic part of the 3D spectrum
up to $k_{\rm trust}$, which is natural to identify with $k_{\rm max}$ from
our 3D fits, and supplementing it with simple power-law functions of
$k_{\parallel}^{2}$ whose free parameters $c_{n}^{\rm finite}$ need to be fit
from the data,

$P^{\rm
EFT}_{\rm{1D}}(k_{\parallel})=\frac{1}{2\pi}\int_{k_{\parallel}}^{k_{\rm
max}}dk~{}k[P_{\rm 3D}-P_{\rm stoch}]+C_{0}^{\rm finite}+C_{1}^{\rm
finite}k_{\parallel}^{2}+C_{2}^{\rm finite}k_{\parallel}^{4}~{}\,.$ (6.8)

The part that contains “$-P_{\rm stoch}$” can be re-absorbed into $c_{0}^{\rm
finite}$, yielding

$P^{\rm
EFT}_{\rm{1D}}(k_{\parallel})=\frac{1}{2\pi}\int_{k_{\parallel}}^{k_{\rm
max}}dk~{}kP_{\rm
3D}(k,k_{\parallel})+\mathcal{C}_{0}+\mathcal{C}_{1}k_{\parallel}^{2}+\mathcal{C}_{2}k_{\parallel}^{4}~{}\,.$
(6.9)

Note that counterterms similar to $\mathcal{C}_{n}$ were first discussed in
the context of regularization of the 1D Lyman alpha power spectrum integrals
in [22]. In particular, our $\mathcal{C}_{0}$ matches the counterterm
$\bar{I}_{0}$ from [22].

Figure 5: Left panel: dimensionless 1D flux power spectrum of the Sherwood
simulations against the best-fit theory model from Eq. (6.9), where all the
EFT parameters except $\mathcal{C}_{n}$’s are kept fixed to their 3D best-fit
values. Right panel: residuals between the theory and the data.

Let us estimate calculable parts of the counter terms $\mathcal{C}_{n}$ in our
Universe for $z=2.8$. To that end we can approximate them assuming $P_{\rm
3D}$ from linear theory up to $k_{\rm max}=3~{}h\text{Mpc}^{-1}$ and $k_{\rm
NL}=5~{}h\text{Mpc}^{-1}$,

$\mathcal{C}_{n}\simeq\Bigg{\\{}\begin{array}[]{lr}&\int^{k_{\rm NL}}_{k_{\rm
max}}\frac{dk}{2\pi}~{}kb_{1}^{2}P_{\rm
lin}(k)=0.02~{}[h^{-1}\text{Mpc}],\quad n=0\,,\\\ &\int^{k_{\rm NL}}_{k_{\rm
max}}\frac{dk}{2\pi}~{}k^{-1}(-2b_{1}b_{\eta})P_{\rm lin}(k)=4\cdot
10^{-3}~{}[h^{-1}\text{Mpc}]^{3},\quad n=1\,,\\\ &\int^{k_{\rm NL}}_{k_{\rm
max}}\frac{dk}{2\pi}~{}k^{-3}b^{2}_{\eta}P_{\rm lin}(k)=3\cdot
10^{-4}~{}[h^{-1}\text{Mpc}]^{5},\quad n=2\,,\end{array}$ (6.10)

where we used the best-fit values for $b_{1}$ and $b_{\eta}$ from Tab. 1. We
see that the power law contributions in the r.h.s. of Eq. (6.9) are suppressed
for small $k_{\parallel}$, so they represent a convergent gradient expansion,
which justifies dropping higher order terms in Eq. (6.9). In a power law
Universe these terms naturally scale as powers of $k_{\parallel}/k_{\rm NL}$.
Using Eq. (2.5) we get the following order-of-magnitude estimates for the 1D
power spectrum

$\begin{split}&\frac{k_{\parallel}}{2\pi}P_{\rm
1D}(k_{\parallel})\sim{k_{\parallel}}\int_{k_{\parallel}}^{k_{\rm
NL}}dk~{}\Bigg{[}\underbrace{\frac{k^{n+1}}{k_{\rm NL}^{n+3}}}_{\rm
linear}+\underbrace{\frac{k^{2n+4}}{k_{\rm NL}^{2n+6}}}_{\rm
1-loop}+\underbrace{\frac{k}{k_{\rm NL}^{3}}+\frac{k^{3}}{k_{\rm
NL}^{5}}+\frac{k{k_{\parallel}}^{2}}{k_{\rm NL}^{5}}}_{\rm
stochastic}\Bigg{]}\\\ &=\underbrace{\frac{{k_{\parallel}}}{k_{\rm
NL}}-\frac{k^{n+3}_{\parallel}}{k_{\rm NL}^{n+3}}}_{\rm
linear}+\underbrace{\frac{{k_{\parallel}}}{k_{\rm
NL}}-\frac{k^{2n+6}_{\parallel}}{k_{\rm NL}^{2n+6}}}_{\rm
1-loop}+\underbrace{\frac{{k_{\parallel}}}{k_{\rm
NL}}-\frac{k^{3}_{\parallel}}{k_{\rm
NL}^{3}}+\frac{{k_{\parallel}}^{3}}{k_{\rm
NL}^{3}}-\frac{{k_{\parallel}}^{5}}{k_{\rm NL}^{5}}}_{\rm
stochastic}\,.\end{split}$ (6.11)

First, we see that each term in this expansion, even the deterministic ones,
produce stochastic-type contributions $\propto{k_{\parallel}}/k_{\rm NL}$ from
the upper limit of the integration. They contribute to the calculable part of
the constant shot noise contribution in 1D. The calculable parts of the
higher-derivative stochastic contributions in 3D scale as
$({k_{\parallel}}/k_{\rm NL})^{3}$ and $({k_{\parallel}}/k_{\rm NL})^{5}$, and
hence are clearly suppressed in the limit ${k_{\parallel}}/k_{\rm NL}\to 0$.
The deterministic parts of the linear and one-loop corrections dominate over
the stochastic contributions, but not so strongly as they did in 3D. In
particular, the proper one loop part scales as $({k_{\parallel}}/k_{\rm
NL})^{1}$, the same way as the shot noise $({k_{\parallel}}/k_{\rm NL})$
terms, in contrast with $({k_{\parallel}}/k_{\rm NL})$ against
$({k_{\parallel}}/k_{\rm NL})^{3}$ in the 3D case, see Eq. (2.7). This
explains why the stochastic counterterms $P_{\rm shot}$, $a_{0}$ and $a_{2}$
were not needed in the 3D fits, but appeared important for the 1D
correlations.

To illustrate how does the EFT model for the 1D power spectrum Eq. (6.9) work
in practice, it is natural to evaluate this expression for the best-fit model
to the 3D power spectrum from the previous section at $k_{\rm
max}=3~{}h\text{Mpc}^{-1}$. The 3D power spectrum was actually fit only up to
$k_{\rm max}\approx 2.5~{}h\text{Mpc}^{-1}$ since there were no data points in
the range $2.5<k/(h\text{Mpc}^{-1})<3$. Extrapolating our 3D best-fit into
this range leads to an additional error, which affects the 1D modes with
$k_{\parallel}\simeq k_{\rm trust}$. We found, however, that this error can be
largely canceled by introducing a higher order counterterm
$\mathcal{C}_{6}{k_{\parallel}}^{6}$. Thus, in addition to the usual 3D EFT
counterterms, we need to introduce 4 free parameters in order to describe the
1D power spectrum at the percent precision. We fit all these extra
counterterms $\mathcal{C}_{0,2,4,6}$ to the 1D power spectrum data extracted
from the same simulation box. The best-fit counterterm piece is given by

$\begin{split}\frac{P_{\rm 1D}^{\rm stoch}}{h^{-1}\text{Mpc}}=&0.085+1.8\cdot
10^{-3}[h^{-1}\text{Mpc}]^{2}k_{\parallel}^{2}\\\ &+3.1\cdot
10^{-4}[h^{-1}\text{Mpc}]^{4}k_{\parallel}^{4}-4.2\cdot
10^{-5}[h^{-1}\text{Mpc}]^{6}k_{\parallel}^{6}~{}\,.\end{split}$ (6.12)

Note that these best fit values are in perfect agreement with our estimates
(6.10). The results are shown in Fig. 5. We observe a perfect match between
the theory and the data within the percent statistical scatter characteristic
to the 1D power spectra of the Sherwood simulations, cf. Fig. 10 of [60]. Note
however, that this good fit is obtained at the price of four extra parameters
(6.12). Given that one can probe only line-of-sight modes, a large number of
nuisance parameters may represent a problem in the context of actual
cosmological analyses of the 1D Lyman-$\alpha$ flux power spectrum. We believe
that the best way to proceed is to determine these parameters from high-
fidelity Ly-$\alpha$ simulations and use these measurements as priors when
analyzing the actual data. Similar ideas have been recently discussed in the
context of galaxy clustering, see e.g. [150, 151, 114, 152]. We leave this
line of research for future work.

##  7 Discussion and conclusions

We have developed a one-loop EFT model for the 3D and 1D power spectra of the
Lyman alpha forest. As we discussed in Sec. 3, this model is identical to that
presented by Desjacques et al. [86] in the context of line-of-sight dependent
selection effects in galaxy bias. We showed that the relevant loop integrals
can be used to efficiently evaluated with the FFTLog method, which allows for
an opportunity to carry out global full-shape analyses of the Lyman alpha data
that include variations of cosmological parameters. This will enable one to
explore the information content of the Lyman alpha forest along the lines of
the recent EFT-based full-shape analyses of the galaxy power spectra [67, 68].
This extension of our work is currently underway.

From the theoretical point of view, it is not immediately obvious that the
line-of-sight dependent selection effect bias model of Desjacques et al. [86]
would exactly match a perturbative expansion of the exponential map of the
optical depth field in the context of the Lyman alpha forest transmission. The
fact that the two should be in direct relationship was first pointed out by
Chen et al. [59]. We have carried out an explicit deviations that shows how
the full EFT model emerges from the exponential map of the optical depth. From
the technical point of view, the EFT expansion gets generated by the
renormalization of the contact operators that stem down from the exponential
map. To underline the importance of renormalization, we also carried out a
naive expansion of the exponential map as one does in SPT. In that case an
incomplete set of line-of-sight dependent operators is generated. Their bias
coefficients are fixed by the optical depth contrast biases and the background
optical depth value. We call this expansion the tree-level Tau model. We
explicitly show this model is inconsistent as the loop corrections break
constraints between the bias coefficients. Our calculation suggests that even
if we ignored the counterterms to start with, the loop corrections would
generate all line-of-sight dependent operators that are missing in the tree-
level Tau model.232323Although this exercise is of purely academic interest,
the cleanest way to explicitly check this would be a calculation of the one-
loop bispectrum of the Lyman alpha forest, which can be done along the lines
of the recent one-loop redshift space galaxy bispectrum calculations [153,
154].

On the data-related side, we have carried out a precision comparison of the
one-loop EFT model and the highly accurate 3D Lyman alpha power spectra
extracted from the Sherwood hydrodynamical simulations. We have found that the
EFT model can successfully describe the broadband shape of 3D and 1D power
spectra, with $\lesssim 1\%$ residuals at $k_{\rm max}=3~{}h\text{Mpc}^{-1}$
for $z=2.8$. Note that at face value, these residuals are smaller than the
ones of the phenomenological models used to fit the same the data in Ref.
[60]. The results of this work, however, cannot be directly compared to ours,
as in contrast to us, Ref. [60] assumed an additional $5\%$ noise floor in
their covariance.

As far as the measured values of the Lyman alpha biases are concerned, we have
found strong deviations from the naive tree-level Tau model. First, we have
observed a breakdown of the tight relationships between quadratic and linear
bias parameters. In addition, we have detected quadratic and cubic line-of-
sight operators missing in the tree-level Tau model. We warn, however, that
our measurements could be affected by degeneracies between the EFT nuisance
parameters, covariance matrix assumptions, and simulations’ systematics. It
would be interesting to carry out a more detailed search for line-of-sight-
dependent biases, especially at the bispectrum level. This line of research is
also motivated by measurements of the Lyman alpha 3-point correlations in
simulations [155] and the real data [7, 156]. We leave this for future work.

Finally, we have also described how to consistently extend the EFT framework
to the 1D power spectrum. Thus, our calculations can be applied to the one
dimensional flux power spectrum data from BOSS [15].

Going forward, it would be important to explicitly implement IR resummation
and use it to quantify the systematics of the BAO measurements. In addition,
it would be interesting to study cross-correlations of galaxies and quasars
with the Lyman alpha forest, similar to Ref. [157, 60]. Another important line
of work is to convert our formalism into configuration space, either using the
FFTLog technique similar to [67], or developing the Lagrangian space EFT,
along the lines of [158, 159, 160, 161, 134, 162]. Eventually, we will have to
account for realistic effects present in the actual data, such as the
photoionization fluctuations and patchy reionization. We leave all these
research directions for future exploration.

#### Acknowledgments.

We would like to thank Stephen Chen, Giovanni Cabass and Matias Zaldarriaga
for useful discussions. We are grateful to Diego Blas, Stephen Chen, Vincent
Desjacques, Andreu Font-Ribera, Mathias Garny, Donghui Jeong, Fabian Schmidt,
and Martin White for their valuable comments on the draft.

##  Appendix A Master Integrals

In this Appendix we present master integrals that are necessary to compute the
selection-dependent one-loop corrections. We will use the notation that
($\nu_{12}\equiv\nu_{1}+\nu_{2}$)

$J[\nu_{1},\nu_{2}]\equiv\frac{1}{8\pi^{3/2}}\frac{\Gamma\left(\frac{3}{2}-\nu_{1}\right)\Gamma\left(\frac{3}{2}-\nu_{2}\right)\Gamma\left(\nu_{12}-\frac{3}{2}\right)}{\Gamma(\nu_{1})\Gamma(\nu_{2})\Gamma(3-\nu_{12})}\,.$
(A.1)

We have:

$\displaystyle
A_{5}=\frac{15}{256}\Bigg{(}J[-5+\nu_{1},\nu_{2}]-5J[-4+\nu_{1},-1+\nu_{2}]-3J[-4+\nu_{1},\nu_{2}]$
$\displaystyle+10J[-3+\nu_{1},-2+\nu_{2}]+4J[-3+\nu_{1},-1+\nu_{2}]+2J[-3+\nu_{1},\nu_{2}]$
$\displaystyle-10J[-2+\nu_{1},-3+\nu_{2}]+6J[-2+\nu_{1},-2+\nu_{2}]$
$\displaystyle+2J[-2+\nu_{1},-1+\nu_{2}]+2J[-2+\nu_{1},\nu_{2}]+5J[-1+\nu_{1},-4+\nu_{2}]$
$\displaystyle-12J[-1+\nu_{1},-3+\nu_{2}]+6J[-1+\nu_{1},-2+\nu_{2}]+4J[-1+\nu_{1},-1+\nu_{2}]$
$\displaystyle-3J[-1+\nu_{1},\nu_{2}]-J[\nu_{1},-5+\nu_{2}]+5J[\nu_{1},-4+\nu_{2}]-10J[\nu_{1},-3+\nu_{2}]$
$\displaystyle+10J[\nu_{1},-2+\nu_{2}]-5J[\nu_{1},-1+\nu_{2}]+J[\nu_{1},\nu_{2}]\Bigg{)}$

$\displaystyle B_{5}$
$\displaystyle=-\frac{35}{128}J[-5+\nu_{1},\nu_{2}]+\frac{175}{128}J[-4+\nu_{1},-1+\nu_{2}]+\frac{25}{128}J[-4+\nu_{1},\nu_{2}]$
$\displaystyle\quad-\frac{175}{64}J[-3+\nu_{1},-2+\nu_{2}]+\frac{25}{32}J[-3+\nu_{1},-1+\nu_{2}]+\frac{5}{64}J[-3+\nu_{1},\nu_{2}]$
$\displaystyle\quad+\frac{175}{64}J[-2+\nu_{1},-3+\nu_{2}]-\frac{225}{64}J[-2+\nu_{1},-2+\nu_{2}]+\frac{45}{64}J[-2+\nu_{1},-1+\nu_{2}]$
$\displaystyle\quad+\frac{5}{64}J[-2+\nu_{1},\nu_{2}]-\frac{175}{128}J[-1+\nu_{1},-4+\nu_{2}]+\frac{125}{32}J[-1+\nu_{1},-3+\nu_{2}]$
$\displaystyle\quad-\frac{225}{64}J[-1+\nu_{1},-2+\nu_{2}]+\frac{25}{32}J[-1+\nu_{1},-1+\nu_{2}]+\frac{25}{128}J[-1+\nu_{1},\nu_{2}]$
$\displaystyle\quad+\frac{35}{128}J[\nu_{1},-5+\nu_{2}]-\frac{175}{128}J[\nu_{1},-4+\nu_{2}]+\frac{175}{64}J[\nu_{1},-3+\nu_{2}]$
$\displaystyle\quad-\frac{175}{64}J[\nu_{1},-2+\nu_{2}]+\frac{175}{128}J[\nu_{1},-1+\nu_{2}]-\frac{35}{128}J[\nu_{1},\nu_{2}],$

$\displaystyle C_{5}$
$\displaystyle=\frac{63}{256}J[-5+\nu_{1},\nu_{2}]-\frac{315}{256}J[-4+\nu_{1},-1+\nu_{2}]+\frac{35}{256}J[-4+\nu_{1},\nu_{2}]$
$\displaystyle\quad+\frac{315}{128}J[-3+\nu_{1},-2+\nu_{2}]-\frac{105}{64}J[-3+\nu_{1},-1+\nu_{2}]+\frac{15}{128}J[-3+\nu_{1},\nu_{2}]$
$\displaystyle\quad-\frac{315}{128}J[-2+\nu_{1},-3+\nu_{2}]+\frac{525}{128}J[-2+\nu_{1},-2+\nu_{2}]-\frac{225}{128}J[-2+\nu_{1},-1+\nu_{2}]$
$\displaystyle\quad+\frac{15}{128}J[-2+\nu_{1},\nu_{2}]+\frac{315}{256}J[-1+\nu_{1},-4+\nu_{2}]-\frac{245}{64}J[-1+\nu_{1},-3+\nu_{2}]$
$\displaystyle\quad+\frac{525}{128}J[-1+\nu_{1},-2+\nu_{2}]-\frac{105}{64}J[-1+\nu_{1},-1+\nu_{2}]+\frac{35}{256}J[-1+\nu_{1},\nu_{2}]$
$\displaystyle\quad-\frac{63}{256}J[\nu_{1},-5+\nu_{2}]+\frac{315}{256}J[\nu_{1},-4+\nu_{2}]-\frac{315}{128}J[\nu_{1},-3+\nu_{2}]$
$\displaystyle\quad+\frac{315}{128}J[\nu_{1},-2+\nu_{2}]-\frac{315}{256}J[\nu_{1},-1+\nu_{2}]+\frac{63}{256}J[\nu_{1},\nu_{2}].$

$\displaystyle A_{6}=-\frac{1}{1024}\times
5\Big{(}J[-6+\nu_{1},\nu_{2}]-6J[-5+\nu_{1},-1+\nu_{2}]-6J[-5+\nu_{1},\nu_{2}]$
$\displaystyle+15J[-4+\nu_{1},-2+\nu_{2}]+18J[-4+\nu_{1},-1+\nu_{2}]+15J[-4+\nu_{1},\nu_{2}]$
$\displaystyle-20J[-3+\nu_{1},-3+\nu_{2}]-12J[-3+\nu_{1},-2+\nu_{2}]-12J[-3+\nu_{1},-1+\nu_{2}]$
$\displaystyle-20J[-3+\nu_{1},\nu_{2}]+15J[-2+\nu_{1},-4+\nu_{2}]-12J[-2+\nu_{1},-3+\nu_{2}]$
$\displaystyle-6J[-2+\nu_{1},-2+\nu_{2}]-12J[-2+\nu_{1},-1+\nu_{2}]+15J[-2+\nu_{1},\nu_{2}]$
$\displaystyle-6J[-1+\nu_{1},-5+\nu_{2}]+18J[-1+\nu_{1},-4+\nu_{2}]-12J[-1+\nu_{1},-3+\nu_{2}]$
$\displaystyle-12J[-1+\nu_{1},-2+\nu_{2}]+18J[-1+\nu_{1},-1+\nu_{2}]-6J[-1+\nu_{1},\nu_{2}]$
$\displaystyle+J[\nu_{1},-6+\nu_{2}]-6J[\nu_{1},-5+\nu_{2}]+15J[\nu_{1},-4+\nu_{2}]$
$\displaystyle-20J[\nu_{1},-3+\nu_{2}]+15J[\nu_{1},-2+\nu_{2}]-6J[\nu_{1},-1+\nu_{2}]+J[\nu_{1},\nu_{2}]\Big{)}$

$\displaystyle B_{6}$
$\displaystyle=\frac{105J[-6+\nu_{1},\nu_{2}]}{1024}-\frac{315}{512}J[-5+\nu_{1},-1+\nu_{2}]-\frac{135}{512}J[-5+\nu_{1},\nu_{2}]+$
$\displaystyle\begin{aligned}
&\frac{1575J[-4+\nu_{1},-2+\nu_{2}]}{1024}+\frac{225}{512}J[-4+\nu_{1},-1+\nu_{2}]+\frac{135J[-4+\nu_{1},\nu_{2}]}{1024}-\\\
&\frac{525}{256}J[-3+\nu_{1},-3+\nu_{2}]+\frac{225}{256}J[-3+\nu_{1},-2+\nu_{2}]+\frac{45}{256}J[-3+\nu_{1},-1+\nu_{2}]+\\\
&\frac{15}{256}J[-3+\nu_{1},\nu_{2}]+\frac{1575J[-2+\nu_{1},-4+\nu_{2}]}{1024}-\frac{675}{256}J[-2+\nu_{1},-3+\nu_{2}]+\\\
&\frac{405}{512}J[-2+\nu_{1},-2+\nu_{2}]+\frac{45}{256}J[-2+\nu_{1},-1+\nu_{2}]+\frac{135J[-2+\nu_{1},\nu_{2}]}{1024}-\\\
&\frac{315}{512}J[-1+\nu_{1},-5+\nu_{2}]-\frac{675}{512}J[-1+\nu_{1},-4+\nu_{2}]+\frac{405}{512}J[-1+\nu_{1},-3+\nu_{2}]+\\\
&\frac{225}{512}J[-1+\nu_{1},-2+\nu_{2}]-\frac{135}{512}J[-1+\nu_{1},-1+\nu_{2}]+\frac{105J[-1+\nu_{1},\nu_{2}]}{1024}+\\\
&\frac{105J[\nu_{1},-6+\nu_{2}]}{1024}-\frac{315}{512}J[\nu_{1},-5+\nu_{2}]+\frac{225}{512}J[\nu_{1},-4+\nu_{2}]-\\\
&\frac{675}{512}J[\nu_{1},-3+\nu_{2}]+\frac{1575J[\nu_{1},-2+\nu_{2}]}{1024}-\frac{315}{512}J[\nu_{1},-1+\nu_{2}]+\frac{105J[\nu_{1},\nu_{2}]}{1024}\end{aligned}$

$\displaystyle
C_{6}=-\frac{315J[-6+\nu_{1},\nu_{2}]}{1024}+\frac{945}{512}J[-5+\nu_{1},-1+\nu_{2}]+\frac{105}{512}J[-5+\nu_{1},\nu_{2}]$
$\displaystyle-\frac{4725J[-4+\nu_{1},-2+\nu_{2}]}{1024}+\frac{525}{512}J[-4+\nu_{1},-1+\nu_{2}]+\frac{75J[-4+\nu_{1},\nu_{2}]}{1024}$
$\displaystyle+\frac{1575}{256}J[-3+\nu_{1},-3+\nu_{2}]-\frac{1575}{256}J[-3+\nu_{1},-2+\nu_{2}]+\frac{225}{256}J[-3+\nu_{1},-1+\nu_{2}]$
$\displaystyle+\frac{15}{256}J[-3+\nu_{1},\nu_{2}]-\frac{4725J[-2+\nu_{1},-4+\nu_{2}]}{1024}+\frac{2625}{256}J[-2+\nu_{1},-3+\nu_{2}]$
$\displaystyle-\frac{3375}{512}J[-2+\nu_{1},-2+\nu_{2}]+\frac{225}{256}J[-2+\nu_{1},-1+\nu_{2}]+\frac{75J[-2+\nu_{1},\nu_{2}]}{1024}$
$\displaystyle+\frac{945}{512}J[-1+\nu_{1},-5+\nu_{2}]-\frac{3675}{512}J[-1+\nu_{1},-4+\nu_{2}]+\frac{2625}{256}J[-1+\nu_{1},-3+\nu_{2}]$
$\displaystyle-\frac{1575}{256}J[-1+\nu_{1},-2+\nu_{2}]+\frac{525}{512}J[-1+\nu_{1},-1+\nu_{2}]+\frac{105}{512}J[-1+\nu_{1},\nu_{2}]$
$\displaystyle-\frac{315J[\nu_{1},-6+\nu_{2}]}{1024}+\frac{945}{512}J[\nu_{1},-5+\nu_{2}]-\frac{4725J[\nu_{1},-4+\nu_{2}]}{1024}$
$\displaystyle+\frac{1575}{256}J[\nu_{1},-3+\nu_{2}]-\frac{4725J[\nu_{1},-2+\nu_{2}]}{1024}+\frac{945}{512}J[\nu_{1},-1+\nu_{2}]$
$\displaystyle-\frac{315J[\nu_{1},\nu_{2}]}{1024};$

$\displaystyle
D_{6}=\frac{231J[-6+\nu_{1},\nu_{2}]}{1024}-\frac{693}{512}J[-5+\nu_{1},-1+\nu_{2}]+\frac{63}{512}J[-5+\nu_{1},\nu_{2}]$
$\displaystyle+\frac{3465J[-4+\nu_{1},-2+\nu_{2}]}{1024}-\frac{945}{512}J[-4+\nu_{1},-1+\nu_{2}]+\frac{105J[-4+\nu_{1},\nu_{2}]}{1024}$
$\displaystyle-\frac{1155}{256}J[-3+\nu_{1},-3+\nu_{2}]+\frac{1575}{256}J[-3+\nu_{1},-2+\nu_{2}]-\frac{525}{256}J[-3+\nu_{1},-1+\nu_{2}]$
$\displaystyle+\frac{25}{256}J[-3+\nu_{1},\nu_{2}]+\frac{3465J[-2+\nu_{1},-4+\nu_{2}]}{1024}-\frac{2205}{256}J[-2+\nu_{1},-3+\nu_{2}]$
$\displaystyle+\frac{3675}{512}J[-2+\nu_{1},-2+\nu_{2}]-\frac{525}{256}J[-2+\nu_{1},-1+\nu_{2}]+\frac{105J[-2+\nu_{1},\nu_{2}]}{1024}$
$\displaystyle-\frac{693}{512}J[-1+\nu_{1},-5+\nu_{2}]+\frac{2835}{512}J[-1+\nu_{1},-4+\nu_{2}]-\frac{2205}{256}J[-1+\nu_{1},-3+\nu_{2}]$
$\displaystyle+\frac{1575}{256}J[-1+\nu_{1},-2+\nu_{2}]-\frac{945}{512}J[-1+\nu_{1},-1+\nu_{2}]+\frac{63}{512}J[-1+\nu_{1},\nu_{2}]$
$\displaystyle+\frac{231J[\nu_{1},-6+\nu_{2}]}{1024}-\frac{693}{512}J[\nu_{1},-5+\nu_{2}]+\frac{3465J[\nu_{1},-4+\nu_{2}]}{1024}$
$\displaystyle-\frac{1155}{256}J[\nu_{1},-3+\nu_{2}]+\frac{3465J[\nu_{1},-2+\nu_{2}]}{1024}-\frac{693}{512}J[\nu_{1},-1+\nu_{2}]+\frac{231J[\nu_{1},\nu_{2}]}{1024};$

$\displaystyle
A_{7}=-\frac{1}{2048}\bigg{(}35J[-7+\nu_{1},\nu_{2}]-7J[-6+\nu_{1},-1+\nu_{2}]-5J[-6+\nu_{1},\nu_{2}]$
$\displaystyle+21J[-5+\nu_{1},-2+\nu_{2}]+18J[-5+\nu_{1},-1+\nu_{2}]+9J[-5+\nu_{1},\nu_{2}]$
$\displaystyle-35J[-4+\nu_{1},-3+\nu_{2}]-15J[-4+\nu_{1},-2+\nu_{2}]-9J[-4+\nu_{1},-1+\nu_{2}]$
$\displaystyle-5J[-4+\nu_{1},\nu_{2}]+35J[-3+\nu_{1},-4+\nu_{2}]-20J[-3+\nu_{1},-3+\nu_{2}]$
$\displaystyle-6J[-3+\nu_{1},-2+\nu_{2}]-4J[-3+\nu_{1},-1+\nu_{2}]-5J[-3+\nu_{1},\nu_{2}]$
$\displaystyle-21J[-2+\nu_{1},-5+\nu_{2}]+45J[-2+\nu_{1},-4+\nu_{2}]-18J[-2+\nu_{1},-3+\nu_{2}]$
$\displaystyle-6J[-2+\nu_{1},-2+\nu_{2}]-9J[-2+\nu_{1},-1+\nu_{2}]+9J[-2+\nu_{1},\nu_{2}]$
$\displaystyle+7J[-1+\nu_{1},-6+\nu_{2}]-30J[-1+\nu_{1},-5+\nu_{2}]+45J[-1+\nu_{1},-4+\nu_{2}]$
$\displaystyle-20J[-1+\nu_{1},-3+\nu_{2}]-15J[-1+\nu_{1},-2+\nu_{2}]+18J[-1+\nu_{1},-1+\nu_{2}]$
$\displaystyle-5J[-1+\nu_{1},\nu_{2}]-J[\nu_{1},-7+\nu_{2}]+7J[\nu_{1},-6+\nu_{2}]$
$\displaystyle-21J[\nu_{1},-5+\nu_{2}]+35J[\nu_{1},-4+\nu_{2}]-35J[\nu_{1},-3+\nu_{2}]$
$\displaystyle+21J[\nu_{1},-2+\nu_{2}]-7J[\nu_{1},-1+\nu_{2}]+J[\nu_{1},\nu_{2}]\bigg{)};$

$\displaystyle
B_{7}=\frac{1}{2048}\bigg{(}105(3J[-7+\nu_{1},\nu_{2}]-21J[-6+\nu_{1},-1+\nu_{2}]-7J[-6+\nu_{1},\nu_{2}]$
$\displaystyle+63J[-5+\nu_{1},-2+\nu_{2}]+14J[-5+\nu_{1},-1+\nu_{2}]+3J[-5+\nu_{1},\nu_{2}]$
$\displaystyle-105J[-4+\nu_{1},-3+\nu_{2}]+35J[-4+\nu_{1},-2+\nu_{2}]+5J[-4+\nu_{1},-1+\nu_{2}]$
$\displaystyle+J[-4+\nu_{1},\nu_{2}]+105J[-3+\nu_{1},-4+\nu_{2}]-140J[-3+\nu_{1},-3+\nu_{2}]$
$\displaystyle+30J[-3+\nu_{1},-2+\nu_{2}]+4J[-3+\nu_{1},-1+\nu_{2}]+J[-3+\nu_{1},\nu_{2}]$
$\displaystyle-63J[-2+\nu_{1},-5+\nu_{2}]+175J[-2+\nu_{1},-4+\nu_{2}]-150J[-2+\nu_{1},-3+\nu_{2}]$
$\displaystyle+30J[-2+\nu_{1},-2+\nu_{2}]+5J[-2+\nu_{1},-1+\nu_{2}]+3J[-2+\nu_{1},\nu_{2}]$
$\displaystyle+21J[-1+\nu_{1},-6+\nu_{2}]-98J[-1+\nu_{1},-5+\nu_{2}]+175J[-1+\nu_{1},-4+\nu_{2}]$
$\displaystyle-140J[-1+\nu_{1},-3+\nu_{2}]+35J[-1+\nu_{1},-2+\nu_{2}]+14J[-1+\nu_{1},-1+\nu_{2}]$
$\displaystyle-7J[-1+\nu_{1},\nu_{2}]-3J[\nu_{1},-7+\nu_{2}]+21J[\nu_{1},-6+\nu_{2}]$
$\displaystyle-63J[\nu_{1},-5+\nu_{2}]+105J[\nu_{1},-4+\nu_{2}]-105J[\nu_{1},-3+\nu_{2}]$
$\displaystyle+63J[\nu_{1},-2+\nu_{2}]-21J[\nu_{1},-1+\nu_{2}]+3J[\nu_{1},\nu_{2}]\bigg{)};$

$\displaystyle
C_{7}[\nu_{1},\nu_{2}]=\frac{1}{2048}\bigg{(}21(33J[-7+\nu_{1},\nu_{2}]-231J[-6+\nu_{1},-1+\nu_{2}]-21J[-6+\nu_{1},\nu_{2}]$
$\displaystyle+693J[-5+\nu_{1},-2+\nu_{2}]-126J[-5+\nu_{1},-1+\nu_{2}]-7J[-5+\nu_{1},\nu_{2}]$
$\displaystyle-1155J[-4+\nu_{1},-3+\nu_{2}]+945J[-4+\nu_{1},-2+\nu_{2}]-105J[-4+\nu_{1},-1+\nu_{2}]$
$\displaystyle-5J[-4+\nu_{1},\nu_{2}]+1155J[-3+\nu_{1},-4+\nu_{2}]-2100J[-3+\nu_{1},-3+\nu_{2}]$
$\displaystyle+1050J[-3+\nu_{1},-2+\nu_{2}]-100J[-3+\nu_{1},-1+\nu_{2}]-5J[-3+\nu_{1},\nu_{2}]$
$\displaystyle-693J[-2+\nu_{1},-5+\nu_{2}]+2205J[-2+\nu_{1},-4+\nu_{2}]-2450J[-2+\nu_{1},-3+\nu_{2}]$
$\displaystyle+1050J[-2+\nu_{1},-2+\nu_{2}]-105J[-2+\nu_{1},-1+\nu_{2}]-7J[-2+\nu_{1},\nu_{2}]$
$\displaystyle+231J[-1+\nu_{1},-6+\nu_{2}]-1134J[-1+\nu_{1},-5+\nu_{2}]+2205J[-1+\nu_{1},-4+\nu_{2}]$
$\displaystyle-2100J[-1+\nu_{1},-3+\nu_{2}]+945J[-1+\nu_{1},-2+\nu_{2}]-126J[-1+\nu_{1},-1+\nu_{2}]$
$\displaystyle-21J[-1+\nu_{1},\nu_{2}]-33J[\nu_{1},-7+\nu_{2}]+231J[\nu_{1},-6+\nu_{2}]$
$\displaystyle-693J[\nu_{1},-5+\nu_{2}]+1155J[\nu_{1},-4+\nu_{2}]-1155J[\nu_{1},-3+\nu_{2}]$
$\displaystyle+693J[\nu_{1},-2+\nu_{2}]-231J[\nu_{1},-1+\nu_{2}]+33J[\nu_{1},\nu_{2}]\bigg{)};$

$\displaystyle
D_{7}=\frac{1}{2048}\bigg{(}429J[-7+\nu_{1},\nu_{2}]-3003J[-6+\nu_{1},-1+\nu_{2}]+231J[-6+\nu_{1},\nu_{2}]$
$\displaystyle+9009J[-5+\nu_{1},-2+\nu_{2}]-4158J[-5+\nu_{1},-1+\nu_{2}]+189J[-5+\nu_{1},\nu_{2}]$
$\displaystyle-15015J[-4+\nu_{1},-3+\nu_{2}]+17325J[-4+\nu_{1},-2+\nu_{2}]-4725J[-4+\nu_{1},-1+\nu_{2}]$
$\displaystyle+175J[-4+\nu_{1},\nu_{2}]+15015J[-3+\nu_{1},-4+\nu_{2}]-32340J[-3+\nu_{1},-3+\nu_{2}]$
$\displaystyle+22050J[-3+\nu_{1},-2+\nu_{2}]-4900J[-3+\nu_{1},-1+\nu_{2}]+175J[-3+\nu_{1},\nu_{2}]$
$\displaystyle-9009J[-2+\nu_{1},-5+\nu_{2}]+31185J[-2+\nu_{1},-4+\nu_{2}]-39690J[-2+\nu_{1},-3+\nu_{2}]$
$\displaystyle+22050J[-2+\nu_{1},-2+\nu_{2}]-4725J[-2+\nu_{1},-1+\nu_{2}]+189J[-2+\nu_{1},\nu_{2}]$
$\displaystyle+3003J[-1+\nu_{1},-6+\nu_{2}]-15246J[-1+\nu_{1},-5+\nu_{2}]+31185J[-1+\nu_{1},-4+\nu_{2}]$
$\displaystyle-32340J[-1+\nu_{1},-3+\nu_{2}]+17325J[-1+\nu_{1},-2+\nu_{2}]-4158J[-1+\nu_{1},-1+\nu_{2}]$
$\displaystyle+231J[-1+\nu_{1},\nu_{2}]-429J[\nu_{1},-7+\nu_{2}]+3003J[\nu_{1},-6+\nu_{2}]$
$\displaystyle-9009J[\nu_{1},-5+\nu_{2}]+15015J[\nu_{1},-4+\nu_{2}]-15015J[\nu_{1},-3+\nu_{2}]$
$\displaystyle+9009J[\nu_{1},-2+\nu_{2}]-3003J[\nu_{1},-1+\nu_{2}]+429J[\nu_{1},\nu_{2}]\bigg{)};$

$\displaystyle
A_{8}=\frac{1}{32768}\bigg{(}35J[-8+\nu_{1},\nu_{2}]-8J[-7+\nu_{1},-1+\nu_{2}]-8J[-7+\nu_{1},\nu_{2}]$
$\displaystyle+28J[-6+\nu_{1},-2+\nu_{2}]+40J[-6+\nu_{1},-1+\nu_{2}]+28J[-6+\nu_{1},\nu_{2}]$
$\displaystyle-56J[-5+\nu_{1},-3+\nu_{2}]-72J[-5+\nu_{1},-2+\nu_{2}]-72J[-5+\nu_{1},-1+\nu_{2}]$
$\displaystyle-56J[-5+\nu_{1},\nu_{2}]+70J[-4+\nu_{1},-4+\nu_{2}]+40J[-4+\nu_{1},-3+\nu_{2}]$
$\displaystyle+36J[-4+\nu_{1},-2+\nu_{2}]+40J[-4+\nu_{1},-1+\nu_{2}]+70J[-4+\nu_{1},\nu_{2}]$
$\displaystyle-56J[-3+\nu_{1},-5+\nu_{2}]+40J[-3+\nu_{1},-4+\nu_{2}]+16J[-3+\nu_{1},-3+\nu_{2}]$
$\displaystyle+16J[-3+\nu_{1},-2+\nu_{2}]+40J[-3+\nu_{1},-1+\nu_{2}]-56J[-3+\nu_{1},\nu_{2}]$
$\displaystyle+28J[-2+\nu_{1},-6+\nu_{2}]-72J[-2+\nu_{1},-5+\nu_{2}]+36J[-2+\nu_{1},-4+\nu_{2}]$
$\displaystyle+16J[-2+\nu_{1},-3+\nu_{2}]+36J[-2+\nu_{1},-2+\nu_{2}]-72J[-2+\nu_{1},-1+\nu_{2}]$
$\displaystyle+28J[-2+\nu_{1},\nu_{2}]-8J[-1+\nu_{1},-7+\nu_{2}]+40J[-1+\nu_{1},-6+\nu_{2}]$
$\displaystyle-72J[-1+\nu_{1},-5+\nu_{2}]+40J[-1+\nu_{1},-4+\nu_{2}]+40J[-1+\nu_{1},-3+\nu_{2}]$
$\displaystyle-72J[-1+\nu_{1},-2+\nu_{2}]+40J[-1+\nu_{1},-1+\nu_{2}]-8J[-1+\nu_{1},\nu_{2}]$
$\displaystyle+J[\nu_{1},-8+\nu_{2}]-8J[\nu_{1},-7+\nu_{2}]+28J[\nu_{1},-6+\nu_{2}]$
$\displaystyle-56J[\nu_{1},-5+\nu_{2}]+70J[\nu_{1},-4+\nu_{2}]-56J[\nu_{1},-3+\nu_{2}]$
$\displaystyle+28J[\nu_{1},-2+\nu_{2}]-8J[\nu_{1},-1+\nu_{2}]+J[\nu_{1},\nu_{2}]\bigg{)};$

$\displaystyle B_{8}[\nu_{1},\nu_{2}]=-\frac{1}{8192}\bigg{(}35\cdot
9J[-8+\nu_{1},\nu_{2}]-72J[-7+\nu_{1},-1+\nu_{2}]-40J[-7+\nu_{1},\nu_{2}]$
$\displaystyle+252J[-6+\nu_{1},-2+\nu_{2}]+168J[-6+\nu_{1},-1+\nu_{2}]+60J[-6+\nu_{1},\nu_{2}]$
$\displaystyle-504J[-5+\nu_{1},-3+\nu_{2}]-168J[-5+\nu_{1},-2+\nu_{2}]-72J[-5+\nu_{1},-1+\nu_{2}]$
$\displaystyle-24J[-5+\nu_{1},\nu_{2}]+630J[-4+\nu_{1},-4+\nu_{2}]-280J[-4+\nu_{1},-3+\nu_{2}]$
$\displaystyle-60J[-4+\nu_{1},-2+\nu_{2}]-24J[-4+\nu_{1},-1+\nu_{2}]-10J[-4+\nu_{1},\nu_{2}]$
$\displaystyle-504J[-3+\nu_{1},-5+\nu_{2}]+840J[-3+\nu_{1},-4+\nu_{2}]-240J[-3+\nu_{1},-3+\nu_{2}]$
$\displaystyle-48J[-3+\nu_{1},-2+\nu_{2}]-24J[-3+\nu_{1},-1+\nu_{2}]-24J[-3+\nu_{1},\nu_{2}]$
$\displaystyle+252J[-2+\nu_{1},-6+\nu_{2}]-840J[-2+\nu_{1},-5+\nu_{2}]+900J[-2+\nu_{1},-4+\nu_{2}]$
$\displaystyle-240J[-2+\nu_{1},-3+\nu_{2}]-60J[-2+\nu_{1},-2+\nu_{2}]-72J[-2+\nu_{1},-1+\nu_{2}]$
$\displaystyle+60J[-2+\nu_{1},\nu_{2}]-72J[-1+\nu_{1},-7+\nu_{2}]+392J[-1+\nu_{1},-6+\nu_{2}]$
$\displaystyle-840J[-1+\nu_{1},-5+\nu_{2}]+840J[-1+\nu_{1},-4+\nu_{2}]-280J[-1+\nu_{1},-3+\nu_{2}]$
$\displaystyle-168J[-1+\nu_{1},-2+\nu_{2}]+168J[-1+\nu_{1},-1+\nu_{2}]-40J[-1+\nu_{1},\nu_{2}]$
$\displaystyle+9J[\nu_{1},-8+\nu_{2}]-72J[\nu_{1},-7+\nu_{2}]+252J[\nu_{1},-6+\nu_{2}]$
$\displaystyle-504J[\nu_{1},-5+\nu_{2}]+630J[\nu_{1},-4+\nu_{2}]-504J[\nu_{1},-3+\nu_{2}]$
$\displaystyle+252J[\nu_{1},-2+\nu_{2}]-72J[\nu_{1},-1+\nu_{2}]+9J[\nu_{1},\nu_{2}]\bigg{)};$

$\displaystyle C_{8}=\frac{1}{16384}\bigg{(}105\cdot
33J[-8+\nu_{1},\nu_{2}]-264J[-7+\nu_{1},-1+\nu_{2}]-72J[-7+\nu_{1},\nu_{2}]$
$\displaystyle+924J[-6+\nu_{1},-2+\nu_{2}]+168J[-6+\nu_{1},-1+\nu_{2}]+28J[-6+\nu_{1},\nu_{2}]$
$\displaystyle-1848J[-5+\nu_{1},-3+\nu_{2}]+504J[-5+\nu_{1},-2+\nu_{2}]+56J[-5+\nu_{1},-1+\nu_{2}]$
$\displaystyle+8J[-5+\nu_{1},\nu_{2}]+2310J[-4+\nu_{1},-4+\nu_{2}]-2520J[-4+\nu_{1},-3+\nu_{2}]$
$\displaystyle+420J[-4+\nu_{1},-2+\nu_{2}]+40J[-4+\nu_{1},-1+\nu_{2}]+6J[-4+\nu_{1},\nu_{2}]$
$\displaystyle-1848J[-3+\nu_{1},-5+\nu_{2}]+4200J[-3+\nu_{1},-4+\nu_{2}]-2800J[-3+\nu_{1},-3+\nu_{2}]$
$\displaystyle+400J[-3+\nu_{1},-2+\nu_{2}]+40J[-3+\nu_{1},-1+\nu_{2}]+8J[-3+\nu_{1},\nu_{2}]$
$\displaystyle+924J[-2+\nu_{1},-6+\nu_{2}]-3528J[-2+\nu_{1},-5+\nu_{2}]+4900J[-2+\nu_{1},-4+\nu_{2}]$
$\displaystyle-2800J[-2+\nu_{1},-3+\nu_{2}]+420J[-2+\nu_{1},-2+\nu_{2}]+56J[-2+\nu_{1},-1+\nu_{2}]$
$\displaystyle+28J[-2+\nu_{1},\nu_{2}]-264J[-1+\nu_{1},-7+\nu_{2}]+1512J[-1+\nu_{1},-6+\nu_{2}]$
$\displaystyle-3528J[-1+\nu_{1},-5+\nu_{2}]+4200J[-1+\nu_{1},-4+\nu_{2}]-2520J[-1+\nu_{1},-3+\nu_{2}]$
$\displaystyle+504J[-1+\nu_{1},-2+\nu_{2}]+168J[-1+\nu_{1},-1+\nu_{2}]-72J[-1+\nu_{1},\nu_{2}]$
$\displaystyle+33J[\nu_{1},-8+\nu_{2}]-264J[\nu_{1},-7+\nu_{2}]+924J[\nu_{1},-6+\nu_{2}]$
$\displaystyle-1848J[\nu_{1},-5+\nu_{2}]+2310J[\nu_{1},-4+\nu_{2}]-1848J[\nu_{1},-3+\nu_{2}]$
$\displaystyle+924J[\nu_{1},-2+\nu_{2}]-264J[\nu_{1},-1+\nu_{2}]+33J[\nu_{1},\nu_{2}]\bigg{)};$

$\displaystyle D_{8}=-\frac{1}{8192}\bigg{(}7\cdot
429J[-8+\nu_{1},\nu_{2}]-3432J[-7+\nu_{1},-1+\nu_{2}]-264J[-7+\nu_{1},\nu_{2}]$
$\displaystyle+12012J[-6+\nu_{1},-2+\nu_{2}]-1848J[-6+\nu_{1},-1+\nu_{2}]-84J[-6+\nu_{1},\nu_{2}]$
$\displaystyle-24024J[-5+\nu_{1},-3+\nu_{2}]+16632J[-5+\nu_{1},-2+\nu_{2}]-1512J[-5+\nu_{1},-1+\nu_{2}]$
$\displaystyle-56J[-5+\nu_{1},\nu_{2}]+30030J[-4+\nu_{1},-4+\nu_{2}]-46200J[-4+\nu_{1},-3+\nu_{2}]$
$\displaystyle+18900J[-4+\nu_{1},-2+\nu_{2}]-1400J[-4+\nu_{1},-1+\nu_{2}]-50J[-4+\nu_{1},\nu_{2}]$
$\displaystyle-24024J[-3+\nu_{1},-5+\nu_{2}]+64680J[-3+\nu_{1},-4+\nu_{2}]-58800J[-3+\nu_{1},-3+\nu_{2}]$
$\displaystyle+19600J[-3+\nu_{1},-2+\nu_{2}]-1400J[-3+\nu_{1},-1+\nu_{2}]-56J[-3+\nu_{1},\nu_{2}]$
$\displaystyle+12012J[-2+\nu_{1},-6+\nu_{2}]-49896J[-2+\nu_{1},-5+\nu_{2}]+79380J[-2+\nu_{1},-4+\nu_{2}]$
$\displaystyle-58800J[-2+\nu_{1},-3+\nu_{2}]+18900J[-2+\nu_{1},-2+\nu_{2}]-1512J[-2+\nu_{1},-1+\nu_{2}]$
$\displaystyle-84J[-2+\nu_{1},\nu_{2}]-3432J[-1+\nu_{1},-7+\nu_{2}]+20328J[-1+\nu_{1},-6+\nu_{2}]$
$\displaystyle-49896J[-1+\nu_{1},-5+\nu_{2}]+64680J[-1+\nu_{1},-4+\nu_{2}]-46200J[-1+\nu_{1},-3+\nu_{2}]$
$\displaystyle+16632J[-1+\nu_{1},-2+\nu_{2}]-1848J[-1+\nu_{1},-1+\nu_{2}]-264J[-1+\nu_{1},\nu_{2}]$
$\displaystyle+429J[\nu_{1},-8+\nu_{2}]-3432J[\nu_{1},-7+\nu_{2}]+12012J[\nu_{1},-6+\nu_{2}]$
$\displaystyle-24024J[\nu_{1},-5+\nu_{2}]+30030J[\nu_{1},-4+\nu_{2}]-24024J[\nu_{1},-3+\nu_{2}]$
$\displaystyle+12012J[\nu_{1},-2+\nu_{2}]-3432J[\nu_{1},-1+\nu_{2}]+429J[\nu_{1},\nu_{2}]\bigg{)};$

$\displaystyle
E_{8}=\frac{1}{32768}(6435J[-8+\nu_{1},\nu_{2}]-51480J[-7+\nu_{1},-1+\nu_{2}]$
$\displaystyle+3432J[-7+\nu_{1},\nu_{2}]+180180J[-6+\nu_{1},-2+\nu_{2}]-72072J[-6+\nu_{1},-1+\nu_{2}]$
$\displaystyle+2772J[-6+\nu_{1},\nu_{2}]-360360J[-5+\nu_{1},-3+\nu_{2}]+360360J[-5+\nu_{1},-2+\nu_{2}]$
$\displaystyle-83160J[-5+\nu_{1},-1+\nu_{2}]+2520J[-5+\nu_{1},\nu_{2}]+450450J[-4+\nu_{1},-4+\nu_{2}]$
$\displaystyle-840840J[-4+\nu_{1},-3+\nu_{2}]+485100J[-4+\nu_{1},-2+\nu_{2}]-88200J[-4+\nu_{1},-1+\nu_{2}]$
$\displaystyle+2450J[-4+\nu_{1},\nu_{2}]-360360J[-3+\nu_{1},-5+\nu_{2}]+1081080J[-3+\nu_{1},-4+\nu_{2}]$
$\displaystyle-1164240J[-3+\nu_{1},-3+\nu_{2}]+529200J[-3+\nu_{1},-2+\nu_{2}]-88200J[-3+\nu_{1},-1+\nu_{2}]$
$\displaystyle+2520J[-3+\nu_{1},\nu_{2}]+180180J[-2+\nu_{1},-6+\nu_{2}]-792792J[-2+\nu_{1},-5+\nu_{2}]$
$\displaystyle+1372140J[-2+\nu_{1},-4+\nu_{2}]-1164240J[-2+\nu_{1},-3+\nu_{2}]+485100J[-2+\nu_{1},-2+\nu_{2}]$