# Mapping Lyman-alpha forest three-dimensional large scale structure in real
and redshift space

Francesco Sinigaglia Department of Physics and Astronomy, Università degli
Studi di Padova, Vicolo dell’Osservatorio 3, I-35122, Padova, Italy  INAF -
Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, I-35122,
Padova, Italy  Instituto de Astrofísica de Canarias, s/n, E-38205, La Laguna,
Tenerife, Spain  Departamento de Astrofísica, Universidad de La Laguna,
E-38206, La Laguna, Tenerife, Spain  Francisco-Shu Kitaura Instituto de
Astrofísica de Canarias, s/n, E-38205, La Laguna, Tenerife, Spain
Departamento de Astrofísica, Universidad de La Laguna, E-38206, La Laguna,
Tenerife, Spain  Andrés Balaguera-Antolínez Instituto de Astrofísica de
Canarias, s/n, E-38205, La Laguna, Tenerife, Spain  Departamento de
Astrofísica, Universidad de La Laguna, E-38206, La Laguna, Tenerife, Spain
Ikkoh Shimizu  Shikoku Gakuin University, 3-2-1 Bunkyocho, Zentsuji, Kagawa,
765-8505, Japan  Kentaro Nagamine Theoretical Astrophysics, Department of
Earth and Space Science, Graduate School of Science, Osaka University,
1-1 Machikaneyama, Toyonaka, Osaka 560-0043, Japan  Kavli-IPMU (WPI),
University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba, 277-8583, Japan
Department of Physics & Astronomy, University of Nevada, Las Vegas, 4505 S.
Maryland Pkwy, Las Vegas, NV 89154-4002, USA  Manuel Sánchez-Benavente
Departamento de Astrofísica, Universidad de La Laguna, E-38206, La Laguna,
Tenerife, Spain  Metin Ata Kavli-IPMU (WPI), University of Tokyo, 5-1-5
Kashiwanoha, Kashiwa, Chiba, 277-8583, Japan

(Accepted February 1, 2022)

###### Abstract

This work presents a new physically-motivated supervised machine learning
method, hydro-bam, to reproduce the three-dimensional Lyman-$\alpha$ forest
field in real and in redshift space learning from a reference hydrodynamic
simulation, thereby saving about 7 orders of magnitude in computing time. We
show that our method is accurate up to $k\sim 1\,h\,\rm{Mpc}^{-1}$ in the one-
(probability distribution function), two- (power-spectra) and three-point (bi-
spectra) statistics of the reconstructed fields. When compared to the
reference simulation including redshift space distortions, our method achieves
deviations of $\lesssim 2\%$ up to $k=0.6\,h\,\rm{Mpc}^{-1}$ in the monopole,
$\lesssim 5\%$ up to $k=0.9\,h\,\rm{Mpc}^{-1}$ in the quadrupole. The bi-
spectrum is well reproduced for triangle configurations with sides up to
$k=0.8\,h\,\rm{Mpc}^{-1}$. In contrast, the commonly-adopted Fluctuating Gunn-
Peterson approximation shows significant deviations already neglecting
peculiar motions at configurations with sides of $k=0.2-0.4\,h\,\rm{Mpc}^{-1}$
in the bi-spectrum, being also significantly less accurate in the power-
spectrum (within 5% up to $k=0.7\,h\,\rm{Mpc}^{-1}$). We conclude that an
accurate analysis of the Lyman-$\alpha$ forest requires considering the
complex baryonic thermodynamical large-scale structure relations. Our
hierarchical domain specific machine learning method can efficiently exploit
this and is ready to generate accurate Lyman-$\alpha$ forest mock catalogues
covering large volumes required by surveys such as DESI and WEAVE.

cosmology, large scale structure — Lyman-$\alpha$ forest — hydrodynamic
simulations — mock catalogs — intergalactic medium

##  1 Introduction

One of the most outstanding discoveries in modern Cosmology has been the
acceleration of the expansion of the Universe (Riess et al., 1998; Perlmutter
et al., 1999). This finding supported the so-called $\Lambda$ cold dark matter
($\Lambda$CDM) model and revived the existence of an accelerating force as
opposed to gravity, called dark energy (see e.g. Peebles & Ratra, 2003, and
references therein), initially postulated by Einstein’s General Relativity.
Since then, $\Lambda$CDM and other alternative theories have been put under
test to unveil the nature of the accelerating expansion and probe the geometry
of our Universe. To this end, the Large Scale Structure of the Universe (LSS
hereafter, see, e.g. Springel et al., 2006) offers invaluable insights into
these profound questions, carrying an imprint of the underlying fundamental
physics which determines the evolution of the Universe and ultimately its
fate.

While most of the efforts have been focused on mapping the Universe on large
scales with galaxies at $z\lesssim 2$, other complementary probes open the
possibility to trace the LSS at higher redshifts.

In particular, the Lyman-$\alpha$ (Ly$\alpha$ hereafter) absorption features
(dubbed Ly$\alpha$ forest), originated in quasar spectra by interactions of
rest-frame UV photons with HI atoms in gaseous clouds along the sightline,
have been increasingly employed to investigate the high-redshift LSS.
Differently from galaxies, which identify discrete points in redshift space,
the Ly$\alpha$ forest is a valuable proxy for the continuous matter density
field along the line-of-sight at $2\lesssim z\lesssim 4$, allowing to use
multiple quasar sightlines to reconstruct the three-dimensional matter density
field (Nusser & Haehnelt, 1999; Pichon et al., 2001; Gallerani et al., 2011;
Kitaura et al., 2012b; Horowitz et al., 2019; Huang et al., 2020; Porqueres et
al., 2020). Detailed cosmological simulations (e.g. Cen et al., 1994;
Hernquist et al., 1996; Miralda-Escudé et al., 1996; Zhang et al., 1998;
Gnedin & Hui, 1998; Davé et al., 1999; Machacek et al., 1999, 2000; Regan et
al., 2007; Norman et al., 2009; Lukić et al., 2015; Nagamine et al., 2021)
indicate that the forest traces preferentially the underdense and mildly
overdense ($\delta\lesssim 1-10$) regions, whereas the absorption lines
rapidly saturate in corresponding high-density gaseous clouds, probing
therefore a complementary density regime to that traced by galaxies.

The Ly$\alpha$ forest has become a valuable cosmological probe, exploiting the
Baryon Oscillation Spectroscopic Survey (BOSS, Eisenstein et al., 2011; Dawson
et al., 2013) $\sim 160,000$ quasars dataset to measure the baryon acoustic
oscillations (BAO) (Slosar et al., 2011; Busca et al., 2013; Slosar et al.,
2013; Kirkby et al., 2013; Font-Ribera et al., 2013, 2014; Delubac et al.,
2015; Bautista et al., 2017; du Mas des Bourboux et al., 2017, 2020) and
recently even proposed to constrain the growth rate (Cuceu et al., 2021). As
surveys, such as the Dark Energy Spectroscopic Instrument (DESI) survey (Levi
et al., 2019), get ready to improve on cosmological measurements, the
modelling of the Ly$\alpha$ forest needs to improve as well.

While the Ly$\alpha$ forest is a linear tracer of the underlying density field
on large cosmological scales ($k<0.1\,h\,\rm{Mpc}^{-1}$), on smaller scales
non-linear corrections become unavoidably necessary (e.g. Chen et al., 2021),
such that modelling the absorption field at high-resolution turns out to be a
highly non-trivial task. Additional effects such as HI self-shielding (Font-
Ribera & Miralda-Escudé, 2012), redshift space distorsions (RSD hereafter, see
Kaiser (1987); Hamilton (1998) and in particular for the Ly$\alpha$ forest
Seljak (2012)), and thermal broadening (Cieplak & Slosar, 2016), among others,
complicate even further the overall analytical framework. The large degree of
complexity accounted for by these features traditionally requires resorting to
full hydrodynamic simulations, at the expense of huge time and computational
resources, making it practically unfeasible to produce the large amount of
independent cosmic volumes necessary to perform cosmological analysis and
exploit in full power the unprecedent incoming datasets from future surveys.

To overcome this problem, a customary approach has been to generate low-
resolution Ly$\alpha$ forest cosmological volumes using scaling relations
(Fluctuating Gunn-Peterson Approximation, FGPA hereafter, Rauch, 1998; Croft
et al., 1998; Weinberg et al., 1999) applied to independent dark matter (DM)
fields, obtained from DM-only N-body simulations (Meiksin & White, 2001; Viel
et al., 2002; Slosar et al., 2009; White et al., 2010; Rorai et al., 2013; Lee
et al., 2014), approximated gravity solvers (Horowitz et al., 2019) (e.g.
FastPM (Feng et al., 2016), COLA (Tassev et al., 2013)), Gaussian random
fields (Le Goff et al., 2011; Font-Ribera et al., 2012; Bautista et al., 2015;
Farr et al., 2020) and rank-ordering transformations from linear theory to
non-linear distributions obtained from a hydrodynamic simulation (Viel et al.,
2002; Greig et al., 2011). Although the FGPA succeeds in roughly accounting
for the spatial correlations on large scales, in the case of $N$-body
simulations it is found to lose accuracy at $k>0.1\,h\,\rm{Mpc^{-1}}$ (as will
be reported in more detail in §2 and 3).

Approaches relying on Gaussian random fields need to be complemented with an
input power-spectrum to impose a priori the two-point target summary
statistics. However, they do not guarantee to accurately reproduce the
realistic three-dimensional cosmic-web matter distribution, as would be
obtained from a simulation, or a detailed analytical structure formation
model. In particular, Font-Ribera et al. (2012) assumes a lognormal model for
the Ly$\alpha$ forest optical depth $\tau$ distribution, obtained applying a
non-linear transform to the Gaussian random field and fitting the free
parameters to match the mean transmitted flux to observations. This
prescription holds in first approximation, although it does not ensure that
the flux distribution accurately matches the observed one. A similar approach
adopted by Le Goff et al. (2011), combining the lognormal model with the FGPA,
features substantial deviations in the $1$D power spectrum at $k\gtrsim
0.4\,h\,\rm{Mpc}^{-1}$ (Le Goff et al., 2011).

The accuracy of the Ly$\alpha$ forest modelling methods is expected to improve
when using (fast) $N$-body codes rather than approximated gravity solvers,
which fail to model the evolution of cosmic structures beyond the mildly non-
linear regime. However, given the high redshift of the data and when the
studies are restricted to large scales of say, above $\sim 5\,h^{-1}\rm{Mpc}$,
fast methods correcting for shell-crossing may be accurate enough (Kitaura &
Hess, 2013; Tosone et al., 2021).

A number of efforts have been dedicated to develop alternative techniques
based on matching just the Ly$\alpha$ forest probability distribution function
(PDF, hereafter) (e.g. the LyMAS code, Peirani et al., 2014) and/or the power-
spectrum (as e.g. in the Iteratively Matched Statistics, IMS, Sorini et al.,
2016). However, in such cases the achieved accuracy in the one- and two-point
summary statistics (IMS: $4\%$ mean deviations in the PDF and $5\%-7\%$ in the
power-spectrum) goes in detriment of the bias relation and bi-spectrum (for a
similar situation modelling the halo distribution, see Kitaura et al., 2015).
An accurate modelling of the two- and three-point statistics is mandatory to
obtain proper covariance matrices that leads to unbiased likelihood analysis
of the data (see e.g., Baumgarten & Chuang, 2018).

Two recent companion papers have proposed to use convolutional neural networks
to synthetize hydrodynamic fields, including the Lyman-$\alpha$ forest,
starting from the dark matter field, using a deterministic (Harrington et al.,
2021) or a probabilistic (HyPhy, Horowitz et al., 2021, and references
therein) deep learning approach. Such methods outperform previous efforts,
although they still feature significant biases already at the level of
reproducing the Ly$\alpha$ forest PDF.

In this work we present the latest application of the Bias Assignment Method
(Balaguera-Antolínez et al., 2018, BAM hereafter) to map the Ly$\alpha$ forest
field in three dimensions from a reference cosmological hydrodynamic
simulation in order to establish a pipeline and yield an accurate (and fast)
alternative to the aforementioned methods for the generations of hydro-
simulations.

In particular, we show that accurate reconstructions of the Ly$\alpha$ forest
cannot be obtained directly from the dark matter field in our framework, but
require using a hierarchical domain specific machine learning approach which
involves the reconstruction of ionized (HII) and neutral (HI) gas density
fields first, as proposed in Sinigaglia et al. (2021) (HydroBAM-I hereafter,
equivalent to BAM-V). Subsequently, we use such fields to model for the first
time the non-trivial mapping of Ly$\alpha$ fluxes, thereby including the
dependence on baryon effects, such as cooling and feedback. We introduce here
also a novel framework for the treatment of redshift-space distortion, not
included in HydroBAM-I. In summary, while HydroBAM-I introduces the general
hierarchical bias showing reconstructions of HII, HI and $T$ fields in real
space, in this work we extend our framework to the complex treatment of
Ly$\alpha$ forest, both in real and in redshift space. In particular, we start
by considering the problem of predicting the Ly$\alpha$ forest field from the
dark matter field and argue that a direct mapping between these two quantities
is not sufficient to achieve the desired accuracy in the power spectrum and
higher-order summary statistics. Then, we show how the bias hierarchy
introduced in HydroBAM-I helps in accomplishing our goal and provide physical
insights to explain this performance. Furthermore, we extend the whole
framework to redshift space, devising a suitable model for redshift space
distortions.

The performance of our method is assessed by means of the deviation of
relevant summary statistics (power- and bi-spectrum in real and redshift
space) with respect to the reference simulation. The results suggest that the
implementation of our method with hydrodynamic simulations on cosmological
volumes will be in position to produce a set of accurate Ly$\alpha$ forest
mock realizations in a fast accurate way, with direct applications to upcoming
surveys such as DESI and WEAVE (Dalton et al., 2012).

The paper is organized as follows. In §2 we briefly introduce the BAM method,
describe the reference cosmological hydrodynamic simulation and summarize the
FGPA. In §3 we present our pipeline, which is applied to Ly$\alpha$ forest
reconstructions in real and redshift space in §4 and §5, respectively.
Possible limitations of our method are discussed in §6. We end with
conclusions in §7.

##  2 Hydro-BAM

In this section, we recap first the machine learning method we use in our
framework. Secondly, we describe the training data set used to learn the
modelling of the Ly$\alpha$ forest, in particular a large hydrodynamic
simulation.

###  2.1 The machine learning method

We build on the BAM algorithm (Balaguera-Antolínez et al., 2018, 2019, BAM-I
and BAM-II hereafter). As stated in HydroBAM-I, this is a special domain
specific machine learning approach which learns the values of a binned
isotropic kernel in Fourier-space, minimizing the cost function defined by the
Mahalanobis distance of the reconstructed power-spectrum and the ground truth
given by the reference simulation. The minimization is achieved through an
iterative rejection Metropolis-Hastings algorithm. The kernel preserves the
dimensionality of the problem, meaning that it is equivalent to a quadratic
matrix in configuration space. Hence, non-local dependencies need to be
explicitly provided within the BAM algorithm, otherwise they will only be
accounted for isotropically through the kernel and the anisotropic
contributions will be lost111We know that the relevant non-local dependencies
are not isotropic (see e.g. Kitaura et al., 2020, and references therein).. In
this sense, BAM represents a physically motivated supervised learning
algorithm. The advantage is that knowing the non-local bias contributions and
hierarchical bias relations, we can speed up the learning process using only
one (sufficiently large) or a few reference simulations with same primordial
power spectrum, setup, and cosmological parameters, e.g. fixed-amplitude
paired simulations (Angulo & Pontzen, 2016, which have shown that just two
fixed-amplitude paired simulations can suppress cosmic variance and achieve
the same accuracy in the power spectrum as tens of traditional simulations),
thereby increasing the accuracy of the reconstruction, as we will show below.
In the latter case, using different simulations is mainly meant at increasing
the effective volume covered by the training set, hence making the bias
relations less prone to be affected by cosmic variance. In this sense, it is
crucial that such different simulations share the same numerical setup,
parameters and parametrizations, as mentioned above, not to introduce
incoherencies in the training set.

The resulting kernels and bias relations (the latter understood as the
probability $\mathcal{P}(\eta|{\Theta})_{\partial V}$222We will omit the
symbol $\partial V$ hereafter. of having a given number of dark matter tracers
$\eta$ in a cell of volume $\partial V$ conditional to a set $\Theta$ of
properties of the dark matter density field, DMDF hereafter in this section,
in the same cell can be used to generate independent realizations of the
tracer spatial distribution, when mapped over independent DMDF (evolved from a
set of initial conditions generated with the cosmological model and parameters
used for the reference simulation). BAM has been shown to be able to generate
ensembles of mock catalogs with an accuracy of $1-10\%$ in the two- and three-
point statistics of dark matter halos (see e.g. BAM-II), with higher accuracy
that previous bias-mapping methods (e.g. Pellejero-Ibañez et al., 2020, BAM-
III) and currently aims at generating the realistic set of galaxy mock
catalogs by including intrinsic properties of tracers such as halo masses,
velocity concentrations, spins (Balaguera-Antolínez & Kitaura, in
preparation).

Our classification of BAM as a machine learning algorithm relies on the fact
that it attempts to solve the problem of assessing the optimal link (bias
relation) between the dark matter density field and some tracer field. By
optimal we mean the relation yielding the most accurate predictions of one-
and two-point statistics of the target tracers field when mapping it from the
dark matter field. In particular, BAM learns how to modify, through a kernel,
the bias initially measured from a reference simulation acting as training
set, applying it, in a second stage, to an independent and approximated DMDF,
in order to predict mock catalogs of the desired target tracers field. As
anticipated, in BAM the cost function is defined as the Mahalanobis distance
between the reconstructed power-spectrum and the ground truth given by the
reference simulation. In this picture, the Metropolis-Hastings algorithm
employed in the iterative determination of the kernel, is the optimization
step usually embedded in the training phase of a machine learning algorithm,
as e.g. gradient descent in a regression problem. However, since we are in the
presence of a stochastic mapping between the dark matter field and its tracers
field, we cannot rely on a deterministic minimization algorithm based e.g. on
gradient computation. Instead, we need to employ an alternative algorithm (the
Metropolis-Hastings in our case) which performs optimization by means of
Monte-Carlo methods.

In HydroBAM-I, the approach of BAM was extended to the context of hydrodynamic
simulations333This application is dubbed Hydro-BAM, but we will sometimes
refer to it throughout the paper as just BAM for shortness., exploring the
link between gas properties and the underlying DMDF, as well as the scaling
relations between different gas properties. According to our findings, the
spatial distribution of HII represents an intermediate stage between the dark
matter density field and other baryon quantities, tracing both the large-scale
matter distribution, as well as small-scale baryonic processes and the
underlying thermodynamic relations.

The BAM procedure can be briefly summarized as follows:

  * •

Start with a DMDF $\delta_{\rm dm}(\vec{r})$ interpolated on a mesh with
resolution $\partial V$.

  * •

For a given dark matter tracer $\eta$, obtain the interpolation on the same
mesh as the DMDF.

  * •

Measure the bias $\mathcal{P}(\eta|\\{\Theta\\})$, where $\\{\Theta\\}$
denotes a set of properties of the underlying dark matter density field
$\delta_{\rm dm}(\vec{r})$ retrieved from the reference simulation. This bias
is obtained from the assessment of the joint distribution
$\mathcal{P}(\eta,\\{\Theta\\})$, which is represented by a multi-dimensional
histogram consisting on $N_{\eta}$ bins of the target tracer property and
$N_{\Theta}$ properties of the DMDF.

  * •

Sample a new variable $\tilde{\eta}$ at each spatial cell (with DM properties
$\\{\tilde{\Theta}\\}$) as

$\tilde{\eta}\curvearrowleft\mathcal{P}(\tilde{\eta}=\eta|\\{\Theta\\}=\\{\tilde{\Theta}\\})$
(1)

  * •

Measure the power-spectrum of this quantity and compare with the same
statistics from the reference.

In general, the first sampling of the tracer field performed according to the
bias relation learnt from the simulation results in biased reconstructions of
the power-spectrum $P_{\tilde{\eta}}(k)$ of the variable $\tilde{\eta}$ with
respect to the reference simulation statistics $P_{\eta}(k)$. Such deviation
is accounted for and corrected iteratively by means of a convolution of the
dark matter density field with an isotropic kernel $\mathcal{K}(k=|\vec{k}|)$,
selected with a Metropolis-Hastings algorithm and updated throughout the
iterations. Briefly, the process to obtain the kernel starts with the
definition of the ratio $T_{0}(k)\equiv
P^{i=0}_{\tilde{\eta}}(k)/P_{\eta}(k)$, where $P^{i=0}_{\tilde{\eta}}(k)$ is
the power-spectrum obtained from the sampling described by Eq. (1). At each
iteration $i$ (producing a new estimate of $P^{i}_{\tilde{\eta}}(k)$), a value
of the ratio $T_{i}(k)$ is computed. At each wavenumber $k$, it is subject to
a Metropolis-Hasting (MH) rejection algorithm444This is done by computing a
transition probability ${\rm min}(1,{\rm
exp}(\mathcal{H}^{2}_{0j}-\mathcal{H}^{2}_{1j}))$ with
$\mathcal{H}_{ij}=(P_{\rm ref}-P^{i}_{\tilde{\eta}})/\sigma$,
$\sigma=\sigma(k)$ the (Gaussian) variance associated to the reference power-
spectrum.. The outcome of the MH sampling process defines a sets of weights
$w_{i}(k)=T_{i}(k)$ if the current value of $T_{i}(k)$ is accepted, or
$w_{i}(k)=1$ otherwise. The BAM kernel is then defined as
$\mathcal{K}(k)\equiv\prod_{j=0}^{j=i}w_{j}(k)$ such that, under successive
convolutions with the DMDF, the limit
$\lim_{n\to\infty}P^{(n)}_{\tilde{\eta}}(k)/P_{\eta}(k)\to 1$ (where $n$ is
the number of iterations) is induced.

The aim of this operation is to transform the underlying DMDF such that, once
the sampling of the tracer-bias over the DMDF is performed, we obtain a tracer
field with the reference power-spectrum. This is achieved in BAM with
typically $\lesssim 1\%$ deviations555Throughout this work, deviations
(sometimes referred to as residuals) between two quantities $A(x)$ and $B(x)$
are computed as $\sum_{i}|A(x_{i})/B(x_{i})-1|/N_{x}$ where the sum is done
over all probed (bins of) values of the variable $x$ and $N_{x}$ is the number
of $x$ values. after $\sim 200$ iterations. The convergence is nevertheless
sensitive to i) the cross-correlation between the DMDF and the tracer field,
ii) the properties of the dark matter density field included in the bias
model, and iii) the way these properties are explored (i.e., binning
strategy). A weak fulfillment of any of these conditions may potentially
undermine the success of the whole procedure (see also HydroBAM-I).

Furthermore, it is desirable that $\mathcal{K}(k)<1~{}\forall k$, since
otherwise it is equivalent to a deconvolution in a similar way to correcting
for aliasing introduced by a mass assignment scheme (Jing, 2005). It is known
that deconvolving a signal in the presence of noise can enhance the noise
contribution666Let us suppose that we aim at recovering a signal $s$ from some
simple data model expressed as: $d={\rm R}s+\epsilon$, where R is a response
function and $\epsilon$ some noise contribution. If R is modelled as a
convolution, then reconstructing the signal with $s_{\rm rec}={\rm
R}^{-1}d=s+{\rm R}^{-1}\epsilon$ will yield reasonable results as long as the
additive noise term remains small. (see e.g. Kitaura & Enßlin, 2008, and
references therein). In our case, the bias relations are subject to the
particular realisation and include stochastic components, which are also
affected by the mathematical representation such as the binning, hence,
introducing a noise component. In the absence of a prior optimal non-linear
transformation which regularises the kernel to values $\mathcal{K}<1$, we
allow the kernel to go beyond $1$, which is theoretically sub-optimal as might
in principle enhance the noise contribution. As we show in §4, however, from a
practical point of view it yields already very good results. A further
improvement on this is left for future work when smaller scales will be
investigated in more detail.

Since BAM is based on a direct measurement of the bias relation (described as
a multi-dimensional histogram777See BAM-II for the details.), devising an
efficient way to represent the different properties is a key ingredient in our
analysis. This is done by means of a suitable scaling and binning, aiming at
efficiently extracting the information encoded in the joint probability
distribution describing the correlations between physical quantities. The
total number of bins used to model the properties in $\\{\Theta\\}$ is $N_{\rm
bins,tot}=\prod_{k=1}^{n}\theta_{k}$, where $\theta_{k}$ is the number of bins
used to model property $k\in\\{\Theta\\}$, and therefore rapidly grows with
$n$. For this reason, we need a sufficiently large volume or several
simulations to avoid over-fitting.

We refer the reader to §4.4 of HydroBAM-I for a preliminary assessment of
overfitting in the context of Hydro-BAM when learning from small data sets. We
further discuss this subject in forthcoming sections.

In comparison to other machine learning approaches, adopting a deep learning
neutral networks perspective (e.g., Harrington et al., 2021; Horowitz et al.,
2021), Hydro-BAM relies on a simpler architecture with no hidden layers. As
mentioned above in this Section, the crucial physical dependencies to map the
desired dark matter tracer field need to be explicitly specified in the bias
model implemented in the code. The kernel, indeed, modifies iteratively the
dark matter field in order to try to compensate the lack of fundamental
physical information in the bias formulation. In particular, previous
applications (BAM-I to BAM-V) have clearly highlighted that including non-
local anisotropic contributions in the bias accounts for significant
improvement of the tracer mapping.

###  2.2 The training data set

The training data set is determined by a reference cosmological hydrodynamic
simulation. This simulation has been obtained using the cosmological smoothed-
particle hydrodynamics (SPH) code GADGET3-OSAKA (Aoyama et al., 2018; Shimizu
et al., 2019; Nagamine et al., 2021)888GADGET3-OSAKA is a modified version of
GADGET-3, an updated version of the popular $N$-body/SPH code GADGET-2
(Springel, 2005), with comoving volume $V=(100\,h^{-1}\text{Mpc})^{3}$,
$N=2\times 512^{3}$ particles of mass $m_{\rm DM}=5.38\times
10^{8}\,h^{-1}\text{M}_{\odot}$ for DM particles and $m_{\rm gas}=1.0\times
10^{8}\,h^{-1}\text{M}_{\odot}$ for gas particles, gravitational softening
length $\epsilon_{g}=7.8\,h^{-1}$ kpc (comoving) and allowing the baryonic
smoothing length to become as small as $0.1\,\epsilon_{g}$. As a result, the
minimum baryonic smoothing at $z=2$ is about physical $260\,h^{-1}$ pc, which
is sufficient to resolve the HI distribution in the circumgalactic medium.
Some basic statistics of Ly$\alpha$ forest have already been discussed in
Nagamine et al. (2021). The simulation includes detailed models for supernova
and stellar feedback (Shimizu et al., 2019), photo-heating and photo-
ionization under a uniform UV background Haardt & Madau (2012), radiative
cooling, chemical evolution of atomic and molecular species by Grackle (Smith
et al., 2016), among others. The simulation is evolved from initial conditions
generated at $z=99$ using MUSIC (Hahn & Abel, 2011), adopting Planck
Collaboration et al. (2016) cosmological parameters.

For the purpose of this work, we consider the output at $z=2$ of the
simulation, interpolating dark matter density, HII density, HI number density
and HI optical depth onto a $128^{3}$ cubic mesh using a CIC mass assignment
scheme (MAS hereafter). In this setup, the physical cell-volume corresponds to
$\partial V\sim(0.78\,h^{-1}\rm{Mpc})^{3}$ and the Nyquist frequency is
$k_{N}=4.02\,h\,\rm{Mpc}^{-1}$.

Figure 1: Phase diagrams showing the scaling relations of the transmitted flux
$F/F_{c}$ as a function of different properties of the gas and dark matter,
from the reference simulation (top row) and the BAM reconstructions (bottom
row). The color intensity indicates the height of the joint probability
distribution between each pair of properties. White regions correspond to
empty two-dimensional bins, i.e. to zero-probability regions. The smoothing is
achieved through a bilinear interpolation over a $100\times 100$ bins grid.

The measurements of power- and bi-spectra are performed within spherical
averages in Fourier space of width equal to the fundamental mode, with proper
corrections for the MAS implemented (through its Fourier transform, see e.g.
Hockney & Eastwood, 1981) and for their (Poisson) shot-noise.

We refer the reader to §2 of HydroBAM-I (and references therein, e.g.,
Nagamine et al., 2021) for a detailed description of the reference simulation
and the properties of the dark matter and baryon fields taken into
consideration.

A major novel aspect of this work with respect to HydroBAM-I is the treatment
of HI optical depth $\tau$ and the corresponding Ly$\alpha$ forest field
999The notation $F=\exp(-\tau)$ for the Ly$\alpha$ forest flux field actually
refers to the transmitted flux $F/F_{\rm c}$, i.e. to the quasar spectrum
normalized to its continuum $F_{\rm c}$. We will omit hereafter the
normalization and denote the transmitted flux just as $F$, unless specified
otherwise.$F=\exp(-\tau)$. In this context, we do not adopt the FGPA, but
rather obtain the HI opacity by means of a line-of-sight integration. The
optical depth is computed as (Nagamine et al., 2021):

$\tau=\frac{\pi e^{2}}{m_{e}c}\sum_{j}f\,\phi(x-x_{j})n_{\rm HI}(x_{j})\Delta
l,$ (2)

where $e$, $m_{e}$, $c$, $n_{\rm HI}$, $f$, $x_{j}$, $\Delta l$ denote
respectively the electron charge, electron mass, speed of light in vacuum, HI
number density, the absorption oscillator strength, the line-of-sight
coordinate of the $j$-th cell and the physical cell size. The Voigt-line
profile $\phi(x)$ in Eq. (2) is provided by the fitting formula of Tasitsiomi
(2006). Where necessary, relevant quantities (e.g. HI number density) are
previously interpolated on the mesh according to the SPH kernel of the
simulation. Coordinates $x_{j}$ of cells along the line-of-sight refer to the
outcome of interpolation of particles based either on their positions in real
space $r_{j}$ or redshift-space $s_{j}=r_{j}+v^{\rm los}_{j}/aH$, where
$v^{\rm los}_{j}$, $a$ and $H$ are the $j$-th particle velocity component
along the line-of-sight, the scale factor and the Hubble parameter at $z=2$,
respectively.

Figure 2: Ly$\alpha$ forest slices through the simulation box in real (top
row) and redshift space (bottom row). From left to right: a slice from
reference simulation, from BAM reconstruction using HII and HI previously
calibrated within BAM, and from FGPA box. These panels evidence how BAM is
optimal at reproducing both large and small-scale structures, in contrast to
the reconstructions performed with the FGPA, where small-scale structures are
evidently suppressed.

###  2.3 Fluctuating Gunn-Peterson Approximation

Hydrodynamic simulations (see e.g., Lukić et al., 2015) show that a
considerable fraction ($>90\%$ in volume and $>50\%$ in mass) of the gas
probed by the Ly$\alpha$ forest, found in regions with mildly-linear
($\delta_{\rm dm}<10$) regime, is in a diffuse state not subject to shock
processes. In these conditions, the gas density $\rho_{\rm gas}$ and
temperature $T_{\rm gas}$ are found to display a tight scaling relation
characterized by a power-law relation $T_{\rm gas}=T_{0}(\rho_{\rm
gas}/\bar{\rho}_{\rm gas})^{\alpha}$ (with $\bar{\rho}_{\rm gas}$ being the
average gas density), whose parameters (amplitude and slope) depend on the
reionization history and on the spectral slope of the UV background. These
vary commonly within the ranges $4000\,\rm{K}\lesssim T_{0}\lesssim
10^{3}\,\rm{K}$ and $0.3\lesssim\alpha\lesssim 0.6$ (see e.g., Hui & Gnedin,
1997). Assuming photo-ionization equilibrium, $\tau\propto n_{\rm HI}$ and
$n_{\rm HI}\propto\rho^{2}T_{\rm gas}^{-0.7}/\Gamma_{\rm HI}$, where $n_{\rm
HI}$ and $\Gamma_{\rm HI}$ are the HI number density photo-ionization rate and
$\tau$ denotes the optical depth, respectively, it is possible to express the
latter as

$\tau=A(1+\delta_{\rm gas})^{\beta},$ (3)

which represents the FGPA. In this expression, $\beta=2-0.7,\alpha\sim 1.6$,
and $A$ is a constant which depends on redshift and on the details of the
hydrodynamics (see e.g. Weinberg et al., 1999). Provided that dark matter and
gas densities are highly correlated in the cool low-density regions (where the
Ly$\alpha$ absorption mostly takes place) $\delta_{\rm gas}$ can be replaced
with $\delta_{\rm dm}$ in Eq. (3). Furthermore, in some cases an additional
Gaussian smoothing is applied to transform the dark matter field into a
pseudo-baryon density field, accounting for the fact that the gas density
field tends to be smoother than dark matter density in correspondence of low-
density dark matter peaks (Gnedin & Hui, 1998; Bryan et al., 1999).

In this work we compute the FGPA based on our reference simulation and compare
it against the results obtained with BAM. To this end, we use the dark matter
field at our fiducial resolution ($0.78\,h^{-1}\,\rm{Mpc}$ with a mesh of
$128^{3}$ cells), comparable to that used for cosmological studies, but
coarser than typical smoothing lengths employed in detailed studies based on
hydrodynamic simulations (e.g., $228\,\rm{kpc}$ in Sorini et al., 2016). We do
not apply a Gaussian smoothing. Following relevant works (e.g., Weinberg et
al., 1999; Viel et al., 2002; Seljak, 2012; Rorai et al., 2013; Lukić et al.,
2015; Cieplak & Slosar, 2016; Horowitz et al., 2019), we use the value
$\beta=1.6$, and $A\sim 0.44$, the latter obtained from a linear fit in the
$\log_{10}(\tau)-\log_{10}(\delta_{\rm dm})$ plane. This ensures that the
resulting scaling relation is consistent with that observed in the reference
simulation.

##  3 The calibration pipeline

We introduce here a change of paradigm with respect to previous studies of the
Ly$\alpha$ forest in which the relation of the optical depth is expressed as a
function of the dark matter field only (e.g., McDonald et al., 2000; McDonald,
2003; Seljak, 2012; Wang et al., 2015; Arinyo-i-Prats et al., 2015; Cieplak &
Slosar, 2016; Givans & Hirata, 2020; Chen et al., 2021). Instead, we establish
a hierarchy of bias relations starting from a dark matter field, continuing
over the ionised and the neutral gas, to the flux directly. We circumvent
hereby the problem of translating optical depths into fluxes with coarse
resolutions. These operations are known to introduce biases, especially taking
redshift space distortions into account. In particular, the commutation of
non-linear operations (as are down-sampling from the native resolution of the
simulation to the resolution adopted in this work, and the mapping from
optical depth to Ly$\alpha$ forest fluxes $F=\exp(-\tau)$) are known to
introduce velocity bias contributions in redshift-space Ly$\alpha$ forest
representations on a mesh (Seljak, 2012). We have explicitly checked with
numerical experiments that at the volume and resolution we are working at,
this bias is negligible and hence commuting the two aforementioned operations
does not imply a problem. This fact offers us the opportunity of applying the
$F=\exp(-\tau)$ transform at $128^{3}$ resolution, thereby yielding directly a
reconstruction $F^{\prime}$ of the final target quantity $F$ and ensuring that
our procedure does not introduce further biases. This may have not been the
case if we had had a reconstruction $\tau^{\prime}$ of the optical depth
$\tau$ as final step of our hierarchical method, and then computed fluxes as
$F^{\prime}=\exp(-\tau^{\prime})$. In fact, in that case our procedure would
not ensure that the summary statistics of $F^{\prime}$ reproduces the
reference one within the required few percent accuracy.

We construct first a proxy of the gas distribution which is fast to compute.
To this end we take the dark matter field from the reference simulation down-
sampled to the chosen resolution. A large number of fast gravity solvers are
now available achieving that accuracy at those high redshifts ($z=2$) on Mpc
scales (Kitaura & Hess, 2013; Tassev et al., 2013; Feng et al., 2016; Tosone
et al., 2021). Then we make a suitable non-linear transform of that field
which improves the bias mapping relation to the targeted tracer distribution.
In particular, we use $\log_{10}(1+X)$ (see Neyrinck et al., 2009; Falck et
al., 2012; Kitaura & Angulo, 2012; Kitaura et al., 2012a, for a motivation),
with $X$ being $\delta_{\rm dm}$, $\delta_{\rm H{\small{II}}{}}$ and
$\delta_{\rm H{\small{I}}{}}$ (although other kind of transforms can be
considered, Kitaura et al., 2020; Sinigaglia et al., 2021).

The key aspect of starting from the dark matter field is that gravitational
evolution has shifted the matter perturbations forming a cosmic web which
correlates with the gas distribution. In fact, as we showed in HydroBAM-I the
correlation between the ionised gas (HII) and the dark matter field is very
high. Despite of this, the bias relation is highly non-linear and non-local.
Hence, in a second step, we apply our machine learning algorithm to extract
the corresponding kernel and multidimensional bias relation.

We must stress that the BAM algorithm learns the statistical relations between
two fields which can be applied to any spatial configuration. One could, in
principle, learn the different relations between the dark matter and the
baryonic quantities, such as HII or HI, independently. Although, the
statistical relations would be individually correct, the thermodynamical
equilibrium relations among the various baryonic quantities would not be
satisfied, as they would have been obtained ignoring each other.

For this reason, the step in which the first baryonic quantity is sampled
fixes the choice of the overall gas distribution. The hierarchy of bias
relations continues with sampling HI from HII. We could continue with the
temperature of the gas, but our studies have shown that it is not necessary to
obtain the Ly$\alpha$ forest fluxes. In particular, introducing the dependence
on the temperature does not yield substantial improvements in the flux power
spectrum calibration, with respect to our fiducial model based just on HII and
HI. Therefore, we believe the temperature does not account for crucial missing
physical dependencies at our mesh resolution. We will consider potential
extensions to the temperature of our model in future works, investigating the
accuracy of Ly$\alpha$ forest reconstruction towards smaller scales than the
ones probed in this paper.

Figure 3: Transmitted flux $F/F_{\rm c}$ skewers along the line-of-sight
extracted from the simulation box at $2$ random positions. The Ly$\alpha$
forest pseudo-spectra reconstructed by BAM using HII and HI fields sampled by
BAM (yellow dotted) are found to qualitatively reproduce the reference skewers
(red solid) more accurately than the FGPA (purple dash-dotted). The spectral
resolution is $l\sim 0.78\,h^{-1}\rm{Mpc}$.

The first part of our procedure is entirely performed in real space, whereas
the Ly$\alpha$ forest mocks are obtained in redshift space. To accomplish this
goal we apply BAM in real space to obtain HII and HI densities, and then map
these reconstructed quantities in redshift space. This last step requires an
adequate modelling of RSD, which will be discussed in §5.1.

The full calibration strategy can be summarized in the following five steps:

  1. 1.

$\tilde{\delta}_{\rm HII}(\vec{r})\curvearrowleft\mathcal{P}\left(\delta_{\rm
HII}(\vec{r})|\\{\Theta(\delta_{\rm
dm}(\vec{r})\otimes\mathcal{K}_{1})\\}\right)\equiv\mathcal{B}_{1}(\vec{r})$

  2. 2.

$\tilde{\delta}_{\rm HI}(\vec{r})\curvearrowleft\mathcal{P}(\delta_{\rm
HI}(\vec{r})|\\{\Theta(\tilde{\delta}_{\rm
HII}(\vec{r})\otimes\mathcal{K}_{2})\\})\equiv\mathcal{B}_{2}(\vec{r})$

  3. 3.

$\tilde{\delta}_{\rm HII}(\vec{s})=\mathcal{S}\left[\tilde{\delta}_{\rm
HII}(\vec{r}),\vec{v}(\vec{r})\right]$

  4. 4.

$\tilde{\delta}_{\rm HI}(\vec{s})=\mathcal{S}\left[\tilde{\delta}_{\rm
HI}(\vec{r}),\vec{v}(\vec{r})\right]$

  5. 5.

$\tilde{F}(\vec{s})\curvearrowleft\mathcal{P}(F(\vec{s})|\\{\Theta(\tilde{\delta}_{\rm
HII}(\vec{s})\otimes\mathcal{K}_{3},\tilde{\delta}_{\rm
HI}(\vec{s}))\\})=\mathcal{B}_{3}(\vec{s})$

where $\tilde{X}$ denotes the reconstruction of the property $X$ and
$\mathcal{S}$ denotes the mapping from real to redshift space. That is, in
steps 1, 2 and 5 the field $\tilde{X}$ on the left-hand side is randomly-
sampled following the stochastic bias relation $\mathcal{B}$ on the right-hand
side, expressing the probability of having a certain value of $X$ in a cell
conditional to the set of properties $\Theta$, after convolving with the
kernel $\mathcal{K}$ the field upon which $\Theta$ is built. As anticipated,
steps 1-2, presented in details in HydroBAM-I, are carried out in real space,
whereas step 5 is performed in redshift space. $\vec{r}$ and $\vec{s}$ stand
for Eulerian coordinates in real and redshift space, respectively. In steps
3-4, $\mathcal{S}$ denotes the mapping from real to redshift space, which will
be presented in details in §5.1. In summary, HII and HI densities are obtained
in real space starting from the dark matter field in steps 1-2, are moved to
redshift space in steps 3-4, and the resulting redshift-space HII and HI
fields are used to sample fluxes in step 5.

Figure 4: Panel (a): distributions of $\log_{10}(1+\delta_{\rm
H{\small{II}}{}})$ in real space from the simulation and as reconstructed by
BAM . Panel (b): distributions of $\log_{10}(1+\delta_{\rm H{\small{I}}{}})$
in real space from the simulation and as reconstructed by BAM. Panel (c):
distribution of Ly$\alpha$ forest fluxes in real space from the reference
simulation, from the FGPA and as reconstructed by BAM based on HII and HI from
the simulation and from reconstructions shown in panels (a) and (b). Panel
(d): distribution of Ly$\alpha$ forest fluxes in redshift space from the
reference simulation, as reconstructed by BAM based on HII and HI from the
simulation and from reconstructions shown in panels (a) and (b) and mapped in
redshift space as explained in §5.1.

The output of the calibration consists therefore in the sets of multi-
dimensional bias relations $\\{\mathcal{B}_{\rm
i}(\vec{x}=\vec{r},\vec{s})\\}$ and kernels $\\{\mathcal{K}_{\rm i}\\}$,
$\rm{i}=1,2,3$. These are then applied to independent dark matter fields to
obtain the corresponding Ly$\alpha$ forest mock realization, following all the
steps depicted above.

In the following sections, we explore the Ly$\alpha$ forest reconstruction in
real and redshift space based on baryon density fields and present the results
of dedicated BAM runs.

##  4 Reconstruction of the Lyman-$\alpha$ forest: real space

In this section we show the performance of the Hydro-BAM approach in real
space.

###  4.1 The bias model

We start analysing the reconstruction procedure in real space, i.e., excluding
steps 3 and 4 in the hierarchical sampling process introduced in §3. This
enables us to assess the accuracy of the method and compare it to the FGPA
before additional physical effects (i.e. RSD) are taken into account.

The top panels of Fig. 1 show the phase diagrams (or bias relations) linking
the transmitted flux $F$ to other gas properties as well as to the dark matter
density in the reference simulation. From this figure it is evident how the
relation between the flux and dark matter density shows to be less
concentrated (large scatter around the mean relation), in contrast to the
scaling relations between the transmitted flux and the HII and HI density
fields. This is an indication that the bias model $\mathcal{P}(F|\,\delta_{\rm
dm}\otimes\mathcal{K})$ can be sub-optimal for our reconstruction procedure.
Indeed, when calibrating the bias relation, the resulting power-spectrum of
the reconstructed flux field displays $>10\%$ deviations throughout all the
probed Fourier modes. This result suggests that an arbitrarily complex (i.e,
non-linear) $\delta_{\rm dm}-F$ scaling relation is not yet sufficient to
capture the full physical correlation between these quantities. Non-local
contributions have to be accounted for, which however cannot rely solely on
the dark matter field.

The reference hydrodynamic simulation includes a series of baryonic processes
(e.g. cooling and feedback by supernovae and active galactic nuclei) relevant
to determine the spatial distribution of the gas, especially of its cold
neutral phase.

We therefore extend the bias model to include baryons, replacing the dark
matter density with gas density fields, i.e, $\mathcal{P}(F|\\{\Theta\\})$,
where
$\\{\Theta\\}=\\{\delta_{\rm{HII}}\otimes\mathcal{K},\delta_{\rm{HI}}\\}$.

In the case in which HII and HI densities are extracted from the reference
simulation, the aforementioned bias model generates reconstructions with
percent accuracy in the two-point statistics of the Ly$\alpha$ forest, as will
be shown in forthcoming sections. The situation in which this bias model uses
the HII and HI fields reconstructed with BAM is more subtle. In the latter
scenario, we obtain $5\%-10\%$ deviations in the flux power-spectrum on small
and intermediate scales ($k\sim 1.0-2.0\,h\,\rm{Mpc}^{-1}$). To reduce these
deviations, we add the dependence on the dark matter density with the model
$\\{\Theta^{\prime}\\}=\\{\delta_{\rm{HII}}\otimes\mathcal{K},\delta_{\rm{HI}},\delta_{\rm{dm}}\\}$.
This does not further improve the BAM calibration with respect to
$\\{\Theta\\}=\\{\delta_{\rm{HII}}\otimes\mathcal{K},\delta_{\rm{HI}}\\}$ when
using baryons from the simulation. However, in the case of reconstructed HII
and HI the explicit dependency with the dark matter density helps to ensure
the large-scale clustering signal, allowing for an overall improvement of the
reconstruction. We therefore use $\\{\Theta\\}$ as our fiducial model when HII
and HI are read from the reference simulation, while we use
$\\{\Theta^{\prime}\\}$ when these fields are reconstructed within BAM.

Adopting two different fiducial models in the two situations poses no conflict
if one considers that $\\{\Theta^{\prime}\\}$ is just an extension of
$\\{\Theta\\}$, alleviating the slight loss of accuracy discussed above and/or
the effect of random noise in the spatial assignment of tracers for
reconstructed fields. In the ideal case where HII and HI are extracted from
the simulation, we regard $\\{\Theta\\}$ and $\\{\Theta^{\prime}\\}$ to carry
an equivalent amount of information of Ly$\alpha$ forest clustering, as
demonstrated by the numerical tests we mentioned.

Figure 5: Comparison between Ly$\alpha$ forest power-spectra in real space
(panels (a) and (c)) and redshift space (panels (b) and (d)) of reference
simulation, BAM reconstruction using density fields from the reference
simulation and BAM reconstruction using BAM-reconstructed density fields. Gray
shaded area in panels (c) and (d) indicate $1\%,2\%$ and $5\%$ deviations from
the reference simulation.

###  4.2 Results

A first visual assessment of the performance of our method in real space is
achieved by comparing slices through the simulation box from the reference,
the reconstruction by BAM and the FGPA Ly$\alpha$ forest fields. This is shown
in the top row of Fig. 2, from where we observe that BAM replicates the flux
structures, both on large and small scales, in good agreement with the
reference. The FGPA, on the other hand, reproduces the large-scale spatial
distribution, whereas it tends to suppress small-scale structures and to
overestimate the high density regions. This can be also observed in Fig. 3,
showing transmitted flux skewers along the line-of-sight extracted at two
random positions from the reference simulation, the BAM reconstructions using
baryon fields and the FGPA. The origin of the sizable differences between
these three cases can be due to the fact that the FGPA accounts only for the
$F-\delta_{\rm dm}$ mean relation and for local dependencies on the dark
matter field, whereas BAM captures the full non-linear non-local stochastic
nature of the tracer-bias.

---
Figure 6: Top plots set: real space Ly$\alpha$ forest reference bi-spectrum
(red solid), bi-spectrum of BAM reconstruction obtained with density fields
from the reference simulation (green dashed), FGPA bi-spectrum (purple dash-
dotted). Bottom plots set: redshift space Ly$\alpha$ forest reference bi-
spectrum (red solid), bi-spectrum of BAM reconstruction obtained with density
fields from the reference simulation mapped in redshift space (blue dashed)
and bi-spectrum of BAM reconstruction obtained with BAM-reconstructed density
fields (brown dash-dotted). BAM well-reproduces the reference bi-spectrum up
to $k\sim 0.8\,h\,\rm{Mpc}^{-1}$ both in real and in redshift space, whereas
the FGPA considerably underestimate the reference bi-spectrum at
$k>0.1\,h\,\rm{Mpc}^{-1}$.

The quality of the HI and HII reconstructions performed by BAM can be firstly
assessed through the distribution function of each property in these two
scenarios. This is shown in panels (a) and (b) of Fig. 4, where we see that
BAM reproduces the bi-modality presented in the HI distribution. On the other
hand, the qualitative differences observed in Fig. 2 between FGPA and BAM are
represented in Fig. 4 panel (c) by an excess of abundance in the low-density
($F\sim 1$, i.e., full transmission) regime of the flux distribution function,
with a correspondent drop of abundance in high-density regions ($F\sim 0$,
full absorption). Conversely, BAM reproduces with very high accuracy the flux
PDF (panel (d)), as it is designed to match the tracer PDF by construction
learning from the reference simulation.

The panel (a) of Fig. 5 and the top set of plots in Fig. 6 show the BAM-
reconstructed (green dashed based on reference HII and HI, yellow dash-dotted
based on reconstructed HII and HI) and the FGPA (purple dash-dotted) power-
spectrum and bi-spectrum101010For the measurements of bi-spectrum we use the
code https://github.com/cheng-zhao/bispec made available by Cheng Zhao.,
compared to the corresponding reference flux field summary statistics (red
solid). As can be seen from panel (c) of Fig. 5, reconstructed power-spectrum
based on reconstructed baryon densities (reference baryon densities) average
and maximum deviations are respectively $<2\%$ and $\sim 3\%$ up to
$k=2.0\,h\,\rm{Mpc}^{-1}$ (up to the Nyquist frequency) with respect to the
reference power-spectrum. The FGPA power-spectrum oscillates around the
reference one with deviations of $\sim 5\%$ at $k<1.0\,h\,\rm{Mpc}^{-1}$ and
is affected by an abrupt drop of power at $k>1.0\,h\,\rm{Mpc}^{-1}$. Regarding
the bi-spectrum, BAM replicates well the reference bi-spectrum for all the
probed configurations, whereas the FGPA bi-spectrum reproduces the qualitative
trend but severely underestimates it already at configurations probing scales
of $k>0.1\,h\,\rm{Mpc}^{-1}$.

Further experiments considering the HII and HI densities separately, i.e.
$\\{\Theta\\}=\\{\delta_{\rm{HII}}\otimes\mathcal{K}\\}$ or
$\\{\Theta\\}=\\{\delta_{\rm{HI}}\otimes\mathcal{K}\\}$, do not yield the same
level of power-spectrum convergence as the case in which the two are used in
combination. This reinforces our argument that a comprehensive modelling of
both large-scale and small-scale physics is required here.

Our results can be understood in more detail by comparing panels (g)-(h) and
(l)-(m) in Fig. 1, which show respectively the relations
$\mathcal{P}(\delta_{\rm HII}|\,\delta_{\rm dm})$ and $\mathcal{P}(\delta_{\rm
HI}|\,\delta_{\rm HII})$ as seen from the simulation and as reconstructed by
BAM. Although BAM reproduces the reference simulation remarkably well and
manages to replicate the two-phases (hot and cool) thermodynamic structure of
the gas over nearly 12 orders of magnitude, the resulting scaling relations
turn out to have slightly more scatter than the reference ones, implying that
they pose a weaker constraint. This translates into broader bias relations
also shown in panels (d) and (f), which need therefore to be complemented with
additional information supplied by the dark matter spatial distribution.

We should bare in mind that the reference hydrodynamic simulation has been run
at a considerably higher resolution (5123 vs 1283). Hence, we cannot expect
these relations to be equal, but to be the optimal ones reproducing the
summary statistics of the various quantities of interest at a given (lower)
resolution (than the reference one). The computing time gain at this
resolution is already of 7 orders of magnitude as compared to the hydrodynamic
simulation. This difference will, however, dramatically increase when going to
larger volumes and/or lower mesh resolutions than the ones studied here. In
particular, assuming to run Hydro-BAM with fixed number of cores and bias
model, the code speed scales approximately as $\mathcal{O}(N_{\rm c})$, with
$N_{\rm c}$ the number of cells used in the grid. Therefore, passing from
$128^{3}$ to $64^{3}$ mesh resolution would imply a gain in computing time of
a factor $\sim 8$ in the case of the reference simulation used in this paper.
Similarly, adopting a simulation with comoving volume
$V=(500\,h^{-1}\,\rm{Mpc})^{3}$, $N_{\rm p}=2048^{3}$ particles, interpolated
on a $N_{\rm c}=128^{3}$ cells regular mesh (corresponding to approximately
$l\sim 3.9\,h^{-1}\,\rm{Mpc}$ comoving cell size) and run with a TreePM N-body
code scaling as $\mathcal{O}(N_{\rm p}\log(N_{\rm p}))$, Hydro-BAM would keep
its performance unchanged whereas the simulation would take approximately
$\sim 78$ times the CPU time required by our reference simulation to run up to
$z=2$. Lastly, BAM performance does not depend on redshift, whereas it takes a
much longer CPU time to evolve the initial conditions to lower redshift, e.g.
up to $z\sim 1$. In general, there are hence no drawbacks in using larger
volume simulations, as long as a reasonable trade-off between number of
particles and mesh resolution is guaranteed. In some unfavourable cases, there
may not be any further gain in computing time with respect to the reference
case investigated in this paper, being anyway left with the 7-orders-of-
magnitude computing time gain quoted above, i.e. HydroBAM is not penalized.

When the HII and HI fields are available from the simulation, the bias is
modelled with just two conditional properties (HII and HI), so that $N_{\rm
bins,tot}\ll N_{\rm cell}$ (or, equivalently, the average number of cells per
$\Theta$-bin is $N_{\Theta}=N_{\rm cell}/N_{\rm bins,tot}\gg 1$) and the
results are therefore not likely to be affected by overfitting. On the other
hand, the case in which BAM learns from reconstructed HII and HI fields, and
includes also information on dark matter in the bias model, is more delicate
and deserves a more careful look. Because the majority of the bins used to
model $\\{\Theta\\}$ in the bias histogram are actually devoid of cells,
following §4.4 of HydroBAM-I we account only for the bins which are filled
with cells instead of considering all the bins used to describe $\Theta$, and
define the effective average number of cells per $\Theta$-bin as
$N_{\Theta}^{\rm eff}=N_{\rm cell}/N_{\rm bins,filled}$. The purpose of this
is to shift the focus only onto the $\Theta$-bins which are filled with cells
and for which is therefore meaningful to define a non-trivial one-dimensional
tracer probability distribution, according to which the subsequent random
sampling is performed. In other words, considering also the $\Theta$-bins
which are devoid of cells would be misleading, as no random assignment of
tracers is performed in such cases.111111In the context of BAM, all the
$\Theta$-bins used to describe regions of zero probability in non-pathological
bias relations are purposeless and can be effectively considered as one single
bin. Increasing the binning resolution in such regions helps just in gaining
accuracy when defining the transition between zero and non-zero probability
regions of the distribution. This computation yields
$N_{\Theta}^{\rm{eff}}\sim 40-50$ depending on the considered iteration within
the calibration procedure, which is reasonably high number to regard that we
do not incur into overfitting.

From a perturbation theory perspective (see e.g. McDonald et al., 2000;
McDonald, 2003; Seljak, 2012), this would correspond to a Taylor expansion of
the flux field at resolution $\Delta l$, which to linear order is expressed as

$\displaystyle F_{\Delta l}/F_{\Delta l}(\delta_{{\rm HII},\Delta
l}=\delta_{{\rm HI},\Delta l}=0)=$ $\displaystyle 1+b_{{\rm HII},\Delta
l}\delta_{{\rm HII},\Delta l}({\rm DM})+b_{{\rm HI},\Delta l}\delta_{{\rm
HI},\Delta l}({\rm HII}({\rm DM}))\,,$

where convenient expansions for $\delta_{\rm HII}$ and $\delta_{\rm HI}$ can
be substituted into the equation to eventually relate it to the dark matter
field. In fact, BAM is able to perform linear combinations of the non-linear
perturbative Taylor expansions of the properties included in
$\\{\Theta\\}$(see Kitaura et al., 2020, for a more detailed discussion on the
connection between BAM and perturbation theory).

The novel aspects to the standard perturbative framework are that our bias
treatment (i) is explicitly based on baryon rather than on dark matter
densities, (ii) that the dependencies are formulated in a multi-variate
fashion, and (iii) that the relations are iteratively learnt from a reference
simulation. Separating the HII and HI components, which could be in principle
unified into one single baryon density field term, allows to treat them
independently, give them different weights, and leave the freedom to build
distinct higher order local and non-local terms, to model different baryonic
phenomena down to highly non-linear scales.

##  5 Reconstruction of the Lyman$-\alpha$ forest: redshift space

In this section we show the performance of the Hydro-BAM approach including
redshift space distortions.

###  5.1 RSD modelling

As anticipated in §3, the second part of the calibration pipeline requires the
transition from real to redshift space. Given that the HII and HI density
fields are reconstructed in real space (as in HydroBAM-I), an explicit
modelling of RSD for such quantities is required at this stage. In order to do
so, we devise a procedure to apply RSD to density fields already interpolated
on a mesh. The procedure we adopt can be summarized in the following steps:

  1. 1.

Assign to each cell $i$ characterized by a local density $\rho_{i}$, a number
$N=16$ pseudo-particles121212In general, the performance of our method is
found to improve with increasing $N$, though rapidly saturate already at
$N=32-64$. The choice $N=16$ represents the optimal trade-off between accuracy
and code speed. with mass $M=\rho_{i}V_{i}/N$, and positions $\vec{r}_{j}$
coincident to the center of the cell;

  2. 2.

Consider the dark matter velocity field $\vec{v}^{\,\rm sim}_{{\rm dm},i}$
from the simulation, interpolated on a mesh of equal size and resolution as
the density fields. Apply RSD and transform real-space $\vec{r}_{j}$ to
redshift space $\vec{s}_{j}$ coordinates by means of the mapping (see e.g.
Kaiser, 1987; Hamilton, 1998)

$\vec{s}_{j}=\vec{r}_{j}+[b_{\rm v}\,(\vec{v}_{{\rm
dm},j}\cdot\hat{r}_{j})\,\hat{r}_{j}]/(aH),$ (5)

where $\vec{v}_{{\rm dm},j}=\vec{v}^{\,\rm coh}_{{\rm dm},j}+\vec{v}^{\,\rm
disp}_{{\rm dm},j}$, with $\vec{v}^{\,\rm coh}_{{\rm dm},j}=\vec{v}^{\,\rm
sim}_{{\rm dm},i}$ is the coherent component of the velocity field, while
$\vec{v}^{\,\rm disp}_{{\rm dm},j}$ accounts for the small-scale quasi-
virialized motions. This last term is originally smoothed-out from the
distribution of gas by the interpolation of particle velocities on the mesh at
physical cell-resolution. Therefore, we restore it assuming that its
components follow a Gaussian distribution of zero mean and variance
$A(1+\delta_{i})^{\alpha}$ with $A$ and $\alpha$ as free parameters (Heß et
al., 2013; Kitaura et al., 2014). The term $b_{\rm v}$ in Eq. (5) is a
velocity bias factor, while $a$ is the scale factor, $H$ is the Hubble
parameter and $\hat{r}=\vec{r}/|\vec{r}|$. We have adopted the distant
observer approximation and translated to redshift space along the z-component.

  3. 3.

Interpolate again the redshift-space pseudo-particles on a grid using a CIC
MAS, with $M_{j}$ as weight.

This description of RSD implies a set of three free parameters $\\{b_{\rm
v},A,\alpha\\}$, which are determined by comparing the signal of the multipole
decomposition of the $3D$ flux power-spectrum $P_{\ell}(k)$, in particular the
quadrupole ($\ell=2$) and hexadecapole ($\ell=4$) 131313Multipoles have been
measured using nbodykit (https://nbodykit.readthedocs.io, Hand et al., 2018)
in Fourier space, leading to $\\{b_{\rm v},A,\alpha\\}\sim\\{0.93,1,1.03\\}$.
In general, the parameters depend on the grid resolution, and hence these
figures might not be the optimal choice in other cases. We will further
address this aspect in future works, in which we will apply our method to
higher-resolution mesh representations of the same reference simulation.

Figure 7: Multipole decomposition of the Ly$\alpha$ forest in Fourier space.
Panels (a) and (b) show respectively the quadrupole $P_{2}(k)$ and the
hexadecapole $P_{4}(k)$ measured from the reference simulation, BAM
reconstruction using density fields from the reference simulation and the BAM
reconstruction using BAM-reconstructed density fields. The bottom panels (c)
and (d) show the ratios between the multipoles form the two BAM
reconstructions and the reference simulation. In these panels, gray shaded
area in indicate the $1\%,2\%$ and $5\%$ deviations from the reference
simulation.

###  5.2 Reconstruction

The final phase of our pipeline consists of reconstructing the Ly$\alpha$
forest field in redshift space, with the bias characterized by the set of
underlying properties
$\\{\Theta\\}=\\{\delta_{\rm{HII}}(\vec{s})\otimes\mathcal{K},\delta_{\rm{HI}}(\vec{s})\\}$.

We first apply BAM to the reference HII and HI density fields (mapped in
redshift space through $\mathcal{S}$), to address the performance of our RSD
model in the bias. Subsequently, we repeat the same using BAM-reconstructed
HII and HI fields, to mimic a realistic mock production context.

Similarly to the case of real space, BAM is found to yield accurate results.
Dedicated experiments reveal that in the latter the outcome suffers larger
deviations than expected. This is due to the progressive loss of accuracy at
$k>1.0\,h\,\rm{Mpc}^{-1}$ affecting the reconstructed HII and HI fields.
Therefore, we re-introduce the dependence in the bias on the dark matter field
($\\{\Theta^{\prime}(\vec{s})\\}=\\{\delta^{z}_{\rm{HII}}(\vec{s})\otimes\mathcal{K},\delta_{\rm{HI}}(\vec{s}),\delta_{\rm{dm}}(\vec{s})\\}$,
as previously discussed in§3), conveniently mapped in redshift space as well.

The preliminary qualitative assessment based on the visual inspection of
slices through the simulation box (Fig. 2, bottom row) highlights that, the
same as in real space, BAM (central panel) succeeds in reproducing the small-
and large-scale structure of the reference slice (left panel) with a high
level of agreement.

A more quantitative assessment can be performed through the analysis of the
reconstructed summary statistics. Panel (d) in Fig. 4 shows that both the
Ly$\alpha$ forest field reconstructions built on the simulation-based (blue
dashed) and BAM-based (brown dash-dotted) HII and HI fields, are found to have
extremely accurate PDF if compared to the reference (red solid).

Panel (a) in Fig. 5 shows the power-spectrum reconstruction of the flux field
performed by BAM, based on HII and HI density fields from the simulation (blue
dashed) and from previous BAM runs (brown dash-dotted) respectively, compared
to the reference flux field power-spectrum (red solid). The panel (b) of the
same figure shows the ratio $P_{\rm mock}(k)/P_{\rm ref}(k)$, indicating that
in this setup BAM is able to learn the bias relation from the simulation
reproducing the reference power-spectrum with average deviations $\lesssim
1.7\%$ ($<1\%$) and maximum deviation $\sim 3\%$ ($\sim 2\%$) up to the
Nyquist frequency (up to $k=1.0\,h\,\rm{Mpc}^{-1}$). Learning the bias
relation from baryon fields previously reconstructed by BAM yields a moderate
loss of accuracy at $k>1.0\,h\,\rm{Mpc}^{-1}$, albeit preserving the accuracy
on large-scales, accounting for average and maximum deviations $\lesssim 2\%$
and $\sim 3.5\%$ respectively at $k\lesssim 1.0\,h\,\rm{Mpc}^{-1}$. As
mentioned, the average mild lack of power trend at $k\gtrsim
1.0\,h\,\rm{Mpc}^{-1}$, dominated by superimposed random noise oscillations,
is due to the progressive loss of accuracy in HII and HI reconstructions,
which unavoidably propagates to the Ly$\alpha$ forest as well. In fact,
modelling the power spectrum at $k>1.0\,h\,\rm{Mpc}^{-1}$ (where AGN and
supernova feedback have a significant effect) in redshift space is a highly
non-trivial task. Hence, to some extent it is not surprising that a loss of
accuracy in BAM reconstructions of HII and HI has a negative impact also on
the subsequent reconstruction of the Ly$\alpha$ forest, based on the previous
two (see also §6).

To further assess the accuracy of our reconstruction and the performance of
RSD modelling, we compute the multipole expansions of the power-spectrum. In
particular, the quadrupole and the hexadecapole of Ly$\alpha$ forest
reconstructions are shown in panels (a) and (b) of Fig. 7, respectively. While
the power-spectrum is constrained in BAM through the iterative procedure, the
quadrupole and hexadecapole are not explicitly constrained. Therefore, the
accuracy in the reproduction of such multipoles is primarily a consequence of
the approach we have adopted for RSD. The reference quadrupole is shown to be
recovered with $<5\%$ maximum deviation at $k\lesssim 1.0\,h\,\rm{Mpc}^{-1}$,
for reconstructions based both on the reference simulation (green dashed) and
on BAM baryon density fields (purple dash-dotted). Average deviations of the
hexadecapole with respect to the reference (red solid) are $\sim 10\%$ and
$\sim 15\%$ at $k\lesssim 0.5\,h\,\rm{Mpc}^{-1}$ when learning from the
simulation gas (blue dashed) and from the BAM gas (brown dash-dotted),
respectively.

The bottom plots set in Fig. 6 shows the reconstructed flux bi-spectrum at
different triangular configurations in Fourier space. BAM credibly reproduces
the reference bi-spectrum (red solid) at $k\lesssim 0.8\,h\,\rm{Mpc}^{-1}$
when building on the HII and HI fields both from the simulation (blue dashed)
and from reconstructions (brown dash-dotted), starting to lose accuracy in
$\theta_{12}\approx[0,1.3]$ only at $k>0.8\,h\,\rm{Mpc}^{-1}$.

As in the case of real space, an assessment of overfitting is needed only in
the case of reconstruction based on HII and HI fields obtained by BAM, in
which the dark matter field is included as well in the bias. Again,
$N_{\Theta}^{\rm{eff}}\sim 40-50$ and we can rule out the possibility that
overfitting is affecting our results.

We consider the RSD modelling acceptable at this stage. A number of
modifications can be introduced to improve the accuracy in the multipoles at
the expense of either introducing more parameters, or another machine learning
step.

For the first strategy, we could consider separate $b_{\rm v}$ velocity bias
for HII, HI, and the DM. Also the corresponding dispersion terms could adopt
different $A$ and $\alpha$ values. These improvements may allow to account for
intrinsic velocity bias differences among these fields. We do not expect these
refinements to introduce dramatic changes in the RSD free parameters we have
estimated in this work. However, leaving the algorithm the freedom of choosing
different values for the parameters related to different fields might improve
even further the RSD modelling, with a consequent potential projected
improvement in the accuracy of our reconstructions.

In a second strategy, we could introduce BAM-like additional machine learning
steps to extract the degree of anisotropy as a function of scale and density
required to minimize a cost function determined by, e.g., the quadrupole or
hexadecapole. In particular, the parameters of our RSD model would be
automatically determined in a machine learning fashion by explicitly requiring
that deviations of higher order multipoles from the reference statistics are
minimized in the cost function, currently not implemented in our algorithm.
Concretely, in this case the new machine learning step would consist in adding
a further optimization stage in the BAM flowchart sketched in §3, repeating
iteratively steps 3, 4, and 5, to determine the values of $b_{v}$, $A$ and
$\alpha$ which optimize the quadrupole and hexadecapole. Extracting such
parameters as a function of scale and density might allow to unveil non-
obvious dependencies and, thus, improve the model.

We leave these series of studies for later work.

##  6 Limitations of the method

In what follows, we summarize the main potential drawbacks and limitations of
our method.

  * •

Availability of a suitable reference simulation: since BAM is in principle
able to self-consistently learn the kernel from just one, or few, reference
simulations, the main limitation to its application is the availability of
such a suitable simulation. In machine learning language, such limitation
translates into the availability of a suitable, sufficiently large, training
dataset. While this condition is unavoidable by construction and represents a
cause of likely failure if not fulfilled, as already mentioned Hydro-BAM
alleviates the need for a rather large training set (i.e. tens of
simulations), commonly shared by other machine learning algorithms through an
explicit modelling of the physical dependencies underlying the bias relation.
In this sense, as already mentioned in §2.1, one or few simulations spanning a
sufficiently large total volume, are enough to guarantee that the learning
processes is successful and not prone to cosmic variance effects. A volume in
our context is sufficiently large when:

    1. 1.

the mode coupling on large scales is correctly included (Crocce & Scoccimarro,
2006; Klypin & Prada, 2018), which agrees with a flat kernel towards low $k$
values (Balaguera-Antolínez et al., 2019);

    2. 2.

the number of cells exceeds the number of parameters of the bias model (i.e.,
the total number of bins) to ensure good enough statistics for each bias
dependency.

  * •

Need for an effective and feasible prescription for the bias: while the
architecture of Hydro-BAM — relying on an explicit modelling of the key
physical dependencies in the bias relation — has been proven to be a major
strength of our method, it can become a drawback if one has little clue of
which dependencies are the leading ones. This is particularly true for non-
local bias dependencies. While a common ML method is able to explore any kind
of non-local dependence without understanding the underlying physics, the
Hydro-BAM approach needs some explicit information such as the tidal field
tensor. An extension of this problem is that Hydro-BAM needs in principle all
relevant physical relations given to reach a certain accuracy. This drawback
could become especially problematic when going to smaller scales, when we need
to specify not only the tidal field tensor, but also the density fluctuation
curvature, or even the velocity shear. A ML algorithm can potentially
automatically explore all dependencies and even blindly handle them with
multiple layers. For the considered scales in this work, we consider
nonetheless that Hydro-BAM is very competitive and could inspire ML
algorithms. Even though all the applications shown thus far in BAM papers are
supported by a long-standing theoretical discussion in the literature,
offering valuable insights into the topic, this might not be the case for
other applications not yet explored. Moreover, a similar issue arises if the
bias relation requires the modelling of quantities beyond the dark matter
density field, which cannot be reliably obtained within BAM. In this case, one
can still extract the needed fields from the reference simulations to carry
out a purely academic study, but will not be able to exploit such information
to generate mock catalogs. E.g., in this work, we make use of the BAM-
reconstructed HII and HI fields to map the Ly$\alpha$ forest. If we had found
other gas properties to play a crucial role in the Ly$\alpha$ forest mapping,
we would have had to go beyond the hierarchical approach presented in
HydroBAM-I and devise new strategies.

  * •

Loss of accuracy in the various BAM hierarchical reconstruction steps: we have
shown in this paper that the application of the whole Hydro-BAM pipeline
results in a very accurate reconstruction in redshift space of the power
spectrum up to $k\sim 1.0\,h^{-1}\,\rm{Mpc}$, and of the bi-spectrum up to
$k\sim 0.8\,h^{-1}\,\rm{Mpc}$. The loss of accuracy at $k\gtrsim
1.0\,h^{-1}\,\rm{Mpc}$ is a natural consequence of the fact that baryon
processes have an important impact in determining the Ly$\alpha$ forest power
spectrum at such frequencies. Therefore, achieving high accuracy in such a
highly-nonlinear regime is very non-trivial, and a loss of accuracy in the
various reconstruction steps progressively forward-propagate until the last
step.

  * •

Dependence on cosmological and astrophysical models of the reference
simulation: as BAM extracts the relations linking different quantities from
the reference simulations, where such quantities have been obtained through
the full $N$-body/SPH calculation, the resulting BAM mock catalogs will
forcefully represent cosmic realizations bound to the input parameters and
models adopted to run the reference training simulation. This implies that, to
perform predictions with different cosmological models (e.g. modified
gravity), or different baryon feedback implementations, it will be necessary
to dispose of a suite of different training simulations, one for each
different sets of models and parameters to be investigated.

##  7 Conclusions

In this paper, we have presented a hierarchical domain specific machine
learning approach to efficiently map the three-dimensional Ly$\alpha$ forest
field starting from the dark matter density field. We have shown that a direct
relation cannot be accurate if we ignore the HII and HI distributions.
Therefore, the machine learning approach has to be first applied to obtain the
HII field, which shows the highest correlation with the dark matter field, and
then in turn that baryonic field has to be used to get the HI field, as
described in Sinigaglia et al. (2021). These fields are combined in a non-
linear, non-local multivariate fashion to obtain the Ly$\alpha$ forest
absorption flux field. In this way, a few percent accuracy is achieved in the
summary statistics.

Our results, first obtained in real space, show maximum and average deviations
in the power-spectrum up to the Nyquist frequency of $\sim 3\%$ and $<1.7\%$,
respectively; and good agreement in the bi-spectra probing scales down to
$k\sim 1.0\,h\,\rm{Mpc}^{-1}$. In contrast, the commonly-adopted FGPA
reproduces the Ly$\alpha$ forest power-spectrum with oscillations of $\sim
5\%$ around the reference up to $k\sim 1.0\,h\,\rm{Mpc}^{-1}$, failing to
match the one-point PDF and is found to progressively lose accuracy in the bi-
spectrum at configurations already with $k>0.1\,h\,\rm{Mpc}^{-1}$.

In order to represent the Ly$\alpha$ forest in redshift space (as probed by
galaxy surveys), we map the HII and HI fields from real to redshift space
through an appropriate modelling of RSD, not included in HydroBAM-I. In
particular, we explicitly account for both large-scale coherent flows and
small-scale quasi-virialized motions contributions.

Having these consistent anisotropic fields we apply the same machine learning
scheme as in real space to reconstruct the Ly$\alpha$ forest field in redshift
space.

The reconstructed Ly$\alpha$ forest power-spectrum presents average deviations
of $\lesssim 1.7\%$ ($\lesssim 2\%$) and maximum deviations of $\sim 2\%$
($\sim 3\%$) up to $k\sim 1.0\,h\,\rm{Mpc}^{-1}$ (up to the Nyquist frequency)
when the reconstruction of fluxes is based on HII and HI density fields
previously reconstructed with BAM following HydroBAM-I (from the simulation).
The reconstructed quadrupole deviates $<5\%$ with respect to that from the
reference at $k<0.9\,h\,\rm{Mpc}^{-1}$ when using BAM-reconstructed baryon
densities and up to $k\sim 2.0\,h\,\rm{Mpc}^{-1}$ when relying on baryon
densities from the reference simulation. The hexadecapole reproduces the
reference with reasonable accuracy considering that it is highly affected by
cosmic variance. Finally, the bi-spectrum is well-reproduced for Fourier
triangular configurations with side lenghts up to $k\sim
0.8\,h\,\rm{Mpc}^{-1}$.

We stress that our hierarchical bias mapping is applied for the first time in
this paper also using in the calibration HII and HI density fields self-
consistently reconstructed within BAM, whereas in HydroBAM-I we showcased only
the ideal situation in which such fields can be read from the simulation.

In forthcoming publications we will extend the technique presented in this
work to a larger volume simulation, aiming at producing mock Ly$\alpha$ forest
realizations to cover the cosmic volumes probed by on-going and future
observational campaigns. In parallel, we will address the bias modelling
towards smaller scales, to improve the accuracy of reconstructed quantities in
Hydro-BAM on scales where relevant astrophysical phenomena, such as feedback
and cooling, plays a major role in determining the spatial distribution of
baryons. As mentioned, in this context we also plan to extend our study to
higher-resolution mesh representation of the same reference simulations, to
get Hydro-BAM ready to explore new avenues and study phenomena which can
impact the Ly$\alpha$ forest summary statistics on very small scales ($k\sim
5-10\,h\,\rm{Mpc}^{-1}$), such as the mass of warm dark matter particles (Viel
et al., 2005, 2013), the mass of neutrinos (Palanque-Delabrouille et al.,
2015) and the thermal state of the intergalactic medium (Garzilli et al.,
2017, and references therein), among others.

Furthermore, our method can be applied in a similar fashion as in this paper
to obtain independent HI density fields in redshift space which, if properly
complemented with instrumental effects, can pave the way to the generation of
accurate and detailed 21cm-line mock catalogs for HI intensity mapping.

We stress that BAM as a particularly efficient machine learning technique has
still plenty of potential applications to be explored in the fields of large-
scale structure and galaxy formation and evolution, limited a priori just by
the availability of a suitable reference simulation.

## Acknowledgments

We acknowledge the anonymous referee for their comments, which have helped to
improve the quality of the manuscript. We warmly thank Marc Huertas-Company
for valuable comments and insights on machine learning. FS acknowledges the
support of the doctoral grant funded by the University of Padova and by the
Italian Ministry of Education, University and Research (MIUR). FS is also
grateful to the _Instituto de Astrofísica de Canarias_ (IAC) for hospitality
and computing resources. FSK and ABA acknowledge the IAC facilities and the
Spanish Ministry of Economy and Competitiveness (MINECO) under the Severo
Ochoa program SEV-2015-0548, AYA2017-89891-P and CEX2019-000920-S grants. FSK
also thanks the RYC2015-18693 grant. K.N. is grateful to Volker Springel for
providing the original version of GADGET-3, on which the GADGET3-Osaka code is
based on. Our numerical simulations and analyses were carried out on the XC50
systems at the Center for Computational Astrophysics (CfCA) of the National
Astronomical Observatory of Japan (NAOJ), the XC40 and the Yukawa-21 at YITP
in Kyoto University, Octopus at the Cybermedia Center, Osaka University, and
Oakforest-PACS at the University of Tokyo as part of the HPCI system Research
Project (hp200041, hp210090). This work is supported in part by the JSPS
KAKENHI Grant Number JP17H01111, 19H05810, 20H00180. K.N. acknowledges the
travel support from the Kavli IPMU, World Premier Research Center Initiative
(WPI), where part of this work was conducted. MA thanks for the support of the
Kavli IPMU fellowship and was supported by JSPS KAKENHI Grant Number
JP21K13911. We acknowledge Cheng Zhao for making available his code to measure
the bi-spectrum.

## References

  * Angulo & Pontzen (2016) Angulo, R. E., & Pontzen, A. 2016, MNRAS, 462, L1, doi: 10.1093/mnrasl/slw098
  * Aoyama et al. (2018) Aoyama, S., Hou, K.-C., Hirashita, H., Nagamine, K., & Shimizu, I. 2018, Monthly Notices of the Royal Astronomical Society, 478, 4905, doi: 10.1093/mnras/sty1431
  * Arinyo-i-Prats et al. (2015) Arinyo-i-Prats, A., Miralda-Escudé, J., Viel, M., & Cen, R. 2015, J. Cosmology Astropart. Phys, 2015, 017, doi: 10.1088/1475-7516/2015/12/017
  * Balaguera-Antolínez et al. (2018) Balaguera-Antolínez, A., Kitaura, F.-S., Pellejero-Ibáñez, M., Zhao, C., & Abel, T. 2018, Monthly Notices of the Royal Astronomical Society: Letters, 483, L58, doi: 10.1093/mnrasl/sly220
  * Balaguera-Antolínez et al. (2019) Balaguera-Antolínez, A., Kitaura, F.-S., Pellejero-Ibáñez, M., et al. 2019, Monthly Notices of the Royal Astronomical Society, 491, 2565, doi: 10.1093/mnras/stz3206
  * Baumgarten & Chuang (2018) Baumgarten, F., & Chuang, C.-H. 2018, MNRAS, 480, 2535, doi: 10.1093/mnras/sty1971
  * Bautista et al. (2015) Bautista, J. E., Bailey, S., Font-Ribera, A., et al. 2015, J. Cosmology Astropart. Phys, 2015, 060, doi: 10.1088/1475-7516/2015/05/060
  * Bautista et al. (2017) Bautista, J. E., Busca, N. G., Guy, J., et al. 2017, A&A, 603, A12, doi: 10.1051/0004-6361/201730533
  * Bryan et al. (1999) Bryan, G. L., Machacek, M., Anninos, P., & Norman, M. L. 1999, ApJ, 517, 13, doi: 10.1086/307173
  * Busca et al. (2013) Busca, N. G., Delubac, T., Rich, J., et al. 2013, A&A, 552, A96, doi: 10.1051/0004-6361/201220724
  * Cen et al. (1994) Cen, R., Miralda-Escudé, J., Ostriker, J. P., & Rauch, M. 1994, ApJ, 437, L9, doi: 10.1086/187670
  * Chen et al. (2021) Chen, S.-F., Vlah, Z., & White, M. 2021, arXiv e-prints, arXiv:2103.13498.  https://arxiv.org/abs/2103.13498
  * Cieplak & Slosar (2016) Cieplak, A. M., & Slosar, A. 2016, J. Cosmology Astropart. Phys, 2016, 016, doi: 10.1088/1475-7516/2016/03/016
  * Crocce & Scoccimarro (2006) Crocce, M., & Scoccimarro, R. 2006, Phys. Rev. D, 73, 063519, doi: 10.1103/PhysRevD.73.063519
  * Croft et al. (1998) Croft, R. A. C., Weinberg, D. H., Katz, N., & Hernquist, L. 1998, in Large Scale Structure: Tracks and Traces, ed. V. Mueller, S. Gottloeber, J. P. Muecket, & J. Wambsganss, 69–75.  https://arxiv.org/abs/astro-ph/9801255
  * Cuceu et al. (2021) Cuceu, A., Font-Ribera, A., Joachimi, B., & Nadathur, S. 2021, MNRAS, 506, 5439, doi: 10.1093/mnras/stab1999
  * Dalton et al. (2012) Dalton, G., Trager, S. C., Abrams, D. C., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8446, Ground-based and Airborne Instrumentation for Astronomy IV, ed. I. S. McLean, S. K. Ramsay, & H. Takami, 84460P, doi: 10.1117/12.925950
  * Davé et al. (1999) Davé, R., Hernquist, L., Katz, N., & Weinberg, D. H. 1999, ApJ, 511, 521, doi: 10.1086/306722
  * Dawson et al. (2013) Dawson, K. S., Schlegel, D. J., Ahn, C. P., et al. 2013, AJ, 145, 10, doi: 10.1088/0004-6256/145/1/10
  * Delubac et al. (2015) Delubac, T., Bautista, J. E., Busca, N. G., et al. 2015, A&A, 574, A59, doi: 10.1051/0004-6361/201423969
  * du Mas des Bourboux et al. (2017) du Mas des Bourboux, H., Le Goff, J.-M., Blomqvist, M., et al. 2017, A&A, 608, A130, doi: 10.1051/0004-6361/201731731
  * du Mas des Bourboux et al. (2020) du Mas des Bourboux, H., Rich, J., Font-Ribera, A., et al. 2020, ApJ, 901, 153, doi: 10.3847/1538-4357/abb085
  * Eisenstein et al. (2011) Eisenstein, D. J., Weinberg, D. H., Agol, E., et al. 2011, AJ, 142, 72, doi: 10.1088/0004-6256/142/3/72
  * Falck et al. (2012) Falck, B. L., Neyrinck, M. C., Aragon-Calvo, M. A., Lavaux, G., & Szalay, A. S. 2012, ApJ, 745, 17, doi: 10.1088/0004-637X/745/1/17
  * Farr et al. (2020) Farr, J., Font-Ribera, A., du Mas des Bourboux, H., et al. 2020, J. Cosmology Astropart. Phys, 2020, 068, doi: 10.1088/1475-7516/2020/03/068
  * Feng et al. (2016) Feng, Y., Chu, M.-Y., Seljak, U., & McDonald, P. 2016, MNRAS, 463, 2273, doi: 10.1093/mnras/stw2123
  * Font-Ribera et al. (2012) Font-Ribera, A., McDonald, P., & Miralda-Escudé, J. 2012, J. Cosmology Astropart. Phys, 2012, 001, doi: 10.1088/1475-7516/2012/01/001
  * Font-Ribera & Miralda-Escudé (2012) Font-Ribera, A., & Miralda-Escudé, J. 2012, J. Cosmology Astropart. Phys, 2012, 028, doi: 10.1088/1475-7516/2012/07/028
  * Font-Ribera et al. (2013) Font-Ribera, A., Arnau, E., Miralda-Escudé, J., et al. 2013, J. Cosmology Astropart. Phys, 2013, 018, doi: 10.1088/1475-7516/2013/05/018
  * Font-Ribera et al. (2014) Font-Ribera, A., Kirkby, D., Busca, N., et al. 2014, J. Cosmology Astropart. Phys, 2014, 027, doi: 10.1088/1475-7516/2014/05/027
  * Gallerani et al. (2011) Gallerani, S., Kitaura, F. S., & Ferrara, A. 2011, MNRAS, 413, L6, doi: 10.1111/j.1745-3933.2011.01020.x
  * Garzilli et al. (2017) Garzilli, A., Boyarsky, A., & Ruchayskiy, O. 2017, Physics Letters B, 773, 258, doi: 10.1016/j.physletb.2017.08.022
  * Givans & Hirata (2020) Givans, J. J., & Hirata, C. M. 2020, Phys. Rev. D, 102, 023515, doi: 10.1103/PhysRevD.102.023515
  * Gnedin & Hui (1998) Gnedin, N. Y., & Hui, L. 1998, MNRAS, 296, 44, doi: 10.1046/j.1365-8711.1998.01249.x
  * Greig et al. (2011) Greig, B., Bolton, J. S., & Wyithe, J. S. B. 2011, MNRAS, 418, 1980, doi: 10.1111/j.1365-2966.2011.19617.x
  * Haardt & Madau (2012) Haardt, F., & Madau, P. 2012, The Astrophysical Journal, 746, 125, doi: 10.1088/0004-637x/746/2/125
  * Hahn & Abel (2011) Hahn, O., & Abel, T. 2011, Monthly Notices of the Royal Astronomical Society, 415, 2101, doi: 10.1111/j.1365-2966.2011.18820.x
  * Hamilton (1998) Hamilton, A. J. S. 1998, Linear Redshift Distortions: a Review, ed. D. Hamilton, Vol. 231, 185, doi: 10.1007/978-94-011-4960-0_17
  * Hand et al. (2018) Hand, N., Feng, Y., Beutler, F., et al. 2018, AJ, 156, 160, doi: 10.3847/1538-3881/aadae0
  * Harrington et al. (2021) Harrington, P., Mustafa, M., Dornfest, M., Horowitz, B., & Lukić, Z. 2021, arXiv e-prints, arXiv:2106.12662.  https://arxiv.org/abs/2106.12662
  * Hernquist et al. (1996) Hernquist, L., Katz, N., Weinberg, D. H., & Miralda-Escudé, J. 1996, ApJ, 457, L51, doi: 10.1086/309899
  * Heß et al. (2013) Heß, S., Kitaura, F.-S., & Gottlöber, S. 2013, MNRAS, 435, 2065, doi: 10.1093/mnras/stt1428
  * Hockney & Eastwood (1981) Hockney, R. W., & Eastwood, J. W. 1981, Computer Simulation Using Particles
  * Horowitz et al. (2021) Horowitz, B., Dornfest, M., Lukić, Z., & Harrington, P. 2021, arXiv e-prints, arXiv:2106.12675.  https://arxiv.org/abs/2106.12675
  * Horowitz et al. (2019) Horowitz, B., Lee, K.-G., White, M., Krolewski, A., & Ata, M. 2019, ApJ, 887, 61, doi: 10.3847/1538-4357/ab4d4c
  * Huang et al. (2020) Huang, L., Croft, R. A. C., & Arora, H. 2020, arXiv e-prints, arXiv:2009.10673.  https://arxiv.org/abs/2009.10673
  * Hui & Gnedin (1997) Hui, L., & Gnedin, N. Y. 1997, MNRAS, 292, 27, doi: 10.1093/mnras/292.1.27
  * Jing (2005) Jing, Y. P. 2005, ApJ, 620, 559, doi: 10.1086/427087
  * Kaiser (1987) Kaiser, N. 1987, MNRAS, 227, 1, doi: 10.1093/mnras/227.1.1
  * Kirkby et al. (2013) Kirkby, D., Margala, D., Slosar, A., et al. 2013, J. Cosmology Astropart. Phys, 2013, 024, doi: 10.1088/1475-7516/2013/03/024
  * Kitaura & Angulo (2012) Kitaura, F.-S., & Angulo, R. E. 2012, MNRAS, 425, 2443, doi: 10.1111/j.1365-2966.2012.21614.x
  * Kitaura et al. (2012a) Kitaura, F.-S., Angulo, R. E., Hoffman, Y., & Gottlöber, S. 2012a, MNRAS, 425, 2422, doi: 10.1111/j.1365-2966.2012.21589.x
  * Kitaura et al. (2020) Kitaura, F.-S., Balaguera-Antolínez, A., Sinigaglia, F., & Pellejero-Ibáñez, M. 2020, The cosmic web connection to the dark matter halo distribution through gravity.  https://arxiv.org/abs/2005.11598
  * Kitaura & Enßlin (2008) Kitaura, F. S., & Enßlin, T. A. 2008, MNRAS, 389, 497, doi: 10.1111/j.1365-2966.2008.13341.x
  * Kitaura et al. (2012b) Kitaura, F.-S., Gallerani, S., & Ferrara, A. 2012b, MNRAS, 420, 61, doi: 10.1111/j.1365-2966.2011.19997.x
  * Kitaura et al. (2015) Kitaura, F.-S., Gil-Marín, H., Scóccola, C. G., et al. 2015, MNRAS, 450, 1836, doi: 10.1093/mnras/stv645
  * Kitaura & Hess (2013) Kitaura, F. S., & Hess, S. 2013, MNRAS, 435, L78, doi: 10.1093/mnrasl/slt101
  * Kitaura et al. (2014) Kitaura, F.-S., Yepes, G., & Prada, F. 2014, MNRAS, 439, L21, doi: 10.1093/mnrasl/slt172
  * Klypin & Prada (2018) Klypin, A., & Prada, F. 2018, MNRAS, 478, 4602, doi: 10.1093/mnras/sty1340
  * Le Goff et al. (2011) Le Goff, J. M., Magneville, C., Rollinde, E., et al. 2011, A&A, 534, A135, doi: 10.1051/0004-6361/201117736
  * Lee et al. (2014) Lee, K.-G., Hennawi, J. F., White, M., Croft, R. A. C., & Ozbek, M. 2014, ApJ, 788, 49, doi: 10.1088/0004-637X/788/1/49
  * Levi et al. (2019) Levi, M. E., Allen, L. E., Raichoor, A., et al. 2019, Ground-based and Airborne Instrumentation for Astronomy V, 9147, 91470S.  https://arxiv.org/abs/1907.10688v1http://arxiv.org/abs/1907.10688
  * Lukić et al. (2015) Lukić, Z., Stark, C. W., Nugent, P., et al. 2015, MNRAS, 446, 3697, doi: 10.1093/mnras/stu2377
  * Machacek et al. (1999) Machacek, M., Bryan, G. L., Anninos, P., et al. 1999, in American Institute of Physics Conference Series, Vol. 470, After the Dark Ages: When Galaxies were Young (the Universe at 2 < Z < 5), ed. S. Holt & E. Smith, 106–109, doi: 10.1063/1.58592
  * Machacek et al. (2000) Machacek, M. E., Bryan, G. L., Meiksin, A., et al. 2000, ApJ, 532, 118, doi: 10.1086/308551
  * McDonald (2003) McDonald, P. 2003, ApJ, 585, 34, doi: 10.1086/345945
  * McDonald et al. (2000) McDonald, P., Miralda-Escudé, J., Rauch, M., et al. 2000, ApJ, 543, 1, doi: 10.1086/317079
  * Meiksin & White (2001) Meiksin, A., & White, M. 2001, MNRAS, 324, 141, doi: 10.1046/j.1365-8711.2001.04264.x
  * Miralda-Escudé et al. (1996) Miralda-Escudé, J., Cen, R., Ostriker, J. P., & Rauch, M. 1996, ApJ, 471, 582, doi: 10.1086/177992
  * Nagamine et al. (2021) Nagamine, K., Shimizu, I., Fujita, K., et al. 2021, ApJ, 914, 66, doi: 10.3847/1538-4357/abfa16
  * Neyrinck et al. (2009) Neyrinck, M. C., Szapudi, I., & Szalay, A. S. 2009, ApJ, 698, L90, doi: 10.1088/0004-637X/698/2/L90
  * Norman et al. (2009) Norman, M. L., Paschos, P., & Harkness, R. 2009, in Journal of Physics Conference Series, Vol. 180, Journal of Physics Conference Series, 012021, doi: 10.1088/1742-6596/180/1/012021
  * Nusser & Haehnelt (1999) Nusser, A., & Haehnelt, M. 1999, MNRAS, 303, 179, doi: 10.1046/j.1365-8711.1999.02272.x
  * Palanque-Delabrouille et al. (2015) Palanque-Delabrouille, N., Yèche, C., Baur, J., et al. 2015, J. Cosmology Astropart. Phys, 2015, 011, doi: 10.1088/1475-7516/2015/11/011
  * Peebles & Ratra (2003) Peebles, P. J., & Ratra, B. 2003, Reviews of Modern Physics, 75, 559, doi: 10.1103/RevModPhys.75.559
  * Peirani et al. (2014) Peirani, S., Weinberg, D. H., Colombi, S., et al. 2014, ApJ, 784, 11, doi: 10.1088/0004-637X/784/1/11
  * Pellejero-Ibañez et al. (2020) Pellejero-Ibañez, M., Balaguera-Antolínez, A., Kitaura, F.-S., et al. 2020, Monthly Notices of the Royal Astronomical Society, 493, 586–593, doi: 10.1093/mnras/staa270
  * Perlmutter et al. (1999) Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, ApJ, 517, 565, doi: 10.1086/307221
  * Pichon et al. (2001) Pichon, C., Vergely, J. L., Rollinde, E., Colombi, S., & Petitjean, P. 2001, MNRAS, 326, 597, doi: 10.1046/j.1365-8711.2001.04595.x
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13, doi: 10.1051/0004-6361/201525830
  * Porqueres et al. (2020) Porqueres, N., Hahn, O., Jasche, J., & Lavaux, G. 2020, A&A, 642, A139, doi: 10.1051/0004-6361/202038482
  * Rauch (1998) Rauch, M. 1998, ARA&A, 36, 267, doi: 10.1146/annurev.astro.36.1.267
  * Regan et al. (2007) Regan, J. A., Haehnelt, M. G., & Viel, M. 2007, MNRAS, 374, 196, doi: 10.1111/j.1365-2966.2006.11132.x
  * Riess et al. (1998) Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009, doi: 10.1086/300499
  * Rorai et al. (2013) Rorai, A., Hennawi, J. F., & White, M. 2013, ApJ, 775, 81, doi: 10.1088/0004-637X/775/2/81
  * Seljak (2012) Seljak, U. 2012, J. Cosmology Astropart. Phys, 2012, 004, doi: 10.1088/1475-7516/2012/03/004
  * Shimizu et al. (2019) Shimizu, I., Todoroki, K., Yajima, H., & Nagamine, K. 2019, Monthly Notices of the Royal Astronomical Society, 484, 2632, doi: 10.1093/mnras/stz098
  * Sinigaglia et al. (2021) Sinigaglia, F., Kitaura, F.-S., Balaguera-Antolínez, A., et al. 2021, ApJ, 921, 66, doi: 10.3847/1538-4357/ac158b
  * Slosar et al. (2009) Slosar, A., Ho, S., White, M., & Louis, T. 2009, J. Cosmology Astropart. Phys, 2009, 019, doi: 10.1088/1475-7516/2009/10/019
  * Slosar et al. (2011) Slosar, A., Font-Ribera, A., Pieri, M. M., et al. 2011, J. Cosmology Astropart. Phys, 2011, 001, doi: 10.1088/1475-7516/2011/09/001
  * Slosar et al. (2013) Slosar, A., Iršič, V., Kirkby, D., et al. 2013, J. Cosmology Astropart. Phys, 2013, 026, doi: 10.1088/1475-7516/2013/04/026
  * Smith et al. (2016) Smith, B. D., Bryan, G. L., Glover, S. C. O., et al. 2016, Monthly Notices of the Royal Astronomical Society, 466, 2217, doi: 10.1093/mnras/stw3291
  * Sorini et al. (2016) Sorini, D., Oñorbe, J., Lukić, Z., & Hennawi, J. F. 2016, ApJ, 827, 97, doi: 10.3847/0004-637X/827/2/97
  * Springel (2005) Springel, V. 2005, Monthly Notices of the Royal Astronomical Society, 364, 1105, doi: 10.1111/j.1365-2966.2005.09655.x
  * Springel et al. (2006) Springel, V., Frenk, C. S., & White, S. D. M. 2006, Nature, 440, 1137, doi: 10.1038/nature04805
  * Tasitsiomi (2006) Tasitsiomi, A. 2006, ApJ, 645, 792, doi: 10.1086/504460
  * Tassev et al. (2013) Tassev, S., Zaldarriaga, M., & Eisenstein, D. J. 2013, Journal of Cosmology and Astroparticle Physics, 2013, doi: 10.1088/1475-7516/2013/06/036
  * Tosone et al. (2021) Tosone, F., Neyrinck, M. C., Granett, B. R., Guzzo, L., & Vittorio, N. 2021, MNRAS, 505, 2999, doi: 10.1093/mnras/stab1517
  * Viel et al. (2013) Viel, M., Becker, G. D., Bolton, J. S., & Haehnelt, M. G. 2013, Phys. Rev. D, 88, 043502, doi: 10.1103/PhysRevD.88.043502
  * Viel et al. (2005) Viel, M., Lesgourgues, J., Haehnelt, M. G., Matarrese, S., & Riotto, A. 2005, Phys. Rev. D, 71, 063534, doi: 10.1103/PhysRevD.71.063534
  * Viel et al. (2002) Viel, M., Matarrese, S., Mo, H. J., Theuns, T., & Haehnelt, M. G. 2002, MNRAS, 336, 685, doi: 10.1046/j.1365-8711.2002.05803.x
  * Wang et al. (2015) Wang, X., Font-Ribera, A., & Seljak, U. 2015, J. Cosmology Astropart. Phys, 2015, 009, doi: 10.1088/1475-7516/2015/04/009
  * Weinberg et al. (1999) Weinberg, D., Burles, S., Croft, R. A. C., et al. 1999, in Evolution of Large Scale Structure : From Recombination to Garching, ed. A. J. Banday, R. K. Sheth, & L. N. da Costa, 346.  https://arxiv.org/abs/astro-ph/9810142
  * White et al. (2010) White, M., Pope, A., Carlson, J., et al. 2010, ApJ, 713, 383, doi: 10.1088/0004-637X/713/1/383
  * Zhang et al. (1998) Zhang, Y., Meiksin, A., Anninos, P., & Norman, M. L. 1998, ApJ, 495, 63, doi: 10.1086/305260