LyAl-Net: A high-efficiency Lyman-$α$ forest simulation with a
  neural network
The inference of cosmological quantities requires accurate and large
hydrodynamical cosmological simulations. Unfortunately, their computational
time can take millions of CPU hours for a modest coverage in cosmological
scales ($\approx (100 {h^{-1}}\,\text{Mpc})^3)$). The possibility to generate
large quantities of mock Lyman-$\alpha$ observations opens up the possibility
of much better control on covariance matrices estimate for cosmological
parameters inference, and on the impact of systematics due to baryonic effects.
We present a machine learning approach to emulate the hydrodynamical simulation
of intergalactic medium physics for the Lyman-$\alpha$ forest called LyAl-Net.
The main goal of this work is to provide highly efficient and cheap simulations
retaining interpretation abilities about the gas field level, and as a tool for
other cosmological exploration. We use a neural network based on the U-net
architecture, a variant of convolutional neural networks, to predict the
neutral hydrogen physical properties, density, and temperature. We train the
LyAl-Net model with the Horizon-noAGN simulation, though using only 9% of the
volume. We also explore the resilience of the model through tests of a transfer
learning framework using cosmological simulations containing different baryonic
feedback. We test our results by analysing one and two-point statistics of
emulated fields in different scenarios, as well as their stochastic properties.
The ensemble average of the emulated Lyman-$\alpha$ forest absorption as a
function of redshift lies within 2.5% of one derived from the full
hydrodynamical simulation. The computation of individual fields from the dark
matter density agrees well with regular physical regimes of cosmological
fields. The results tested on IllustrisTNG100 showed a drastic improvement in
the Lyman-$\alpha$ forest flux without arbitrary rescaling.
Chotipan Boonkongkird, Guilhem Lavaux, Sebastien Peirani, Yohan Dubois, Natalia Porqueres, Eleni Tsaprazi
# LyAI-Net: A high-efficiency Lyman-\(\alpha\) forest simulation with a neural network

###### Abstract

Context:The inference of cosmological quantities requires accurate and large hydrodynamical cosmological simulations. Unfortunately, their computational time can take millions of CPU hours for a modest coverage in cosmological scales (\(\approx(100h^{-1}\,\mathrm{Mpc})^{3}\)). The possibility to generate large quantities of mock Lyman-\(\alpha\) observations opens up the possibility of much better control on covariance matrices estimate for cosmological parameters inference, and on the impact of systematics due to baryonic effects.

Aims:We present a machine learning approach to emulate the hydrodynamical simulation of intergalactic medium physics for the Lyman-\(\alpha\) forest called _LyAI-Net_. The main goal of this work is to provide highly efficient and cheap simulations retaining interpretation abilities about the gas field level, and as a tool for other cosmological exploration.

Methods:We use a neural network based on the U-net architecture, a variant of convolutional neural networks, to predict the neutral hydrogen physical properties, density, and temperature. We train the _LyAI-Net_ model with the Horizon-noAGN simulation, though using only 9% of the volume. We also explore the resilience of the model through tests of a transfer learning framework using cosmological simulations containing different baryonic feedback. We test our results by analysing one and two-point statistics of emulated fields in different scenarios, as well as their stochastic properties.

Results:The ensemble average of the emulated Lyman-\(\alpha\) forest absorption as a function of redshift lies within 2.5% of one derived from the full hydrodynamical simulation. The computation of individual fields from the dark matter density agrees well with regular physical regimes of cosmological fields. The results tested on Illustris TNG100 showed a drastic improvement in the Lyman-\(\alpha\) forest flux without arbitrary rescaling.

Conclusions:Such emulators could be critical for the exploitation of upcoming surveys like WEAVE-QSO. The transfer learning technique is showing promises to alleviate sensitivity to the learned physical model. This technique could have a decisive impact on the results derived from present experiments, such as QSO surveys derived from SDSS3, which has a resolution power \(R=1500\), and SDSS4 (\(R=2000\)) data, which are the resolution ranges available with our emulator.

## 1 Introduction

The Lyman-\(\alpha\) forest is a collection of sharp absorption features, generally observed from distant Quasi-Stellar object (QSO) spectra, first observed by Lynds (1971). This phenomenon arises from the photons emitted from background sources, with frequencies higher or equal to Lyman-\(\alpha\) frequency, undergoing redshift due to the Hubble expansion and absorbed by the intervening intergalactic medium. Since the absorptions occur at different distances from an observer in a line-of-sight, they are also redshifted. This mechanism creates a forest-like ensemble of features on the observed QSO broadband Lyman-\(\alpha\) emission. While the process only happens with neutral hydrogen (HI) atoms, which represent a small fraction of the intergalactic medium (IGM), one can use the properties of the Lyman-\(\alpha\) forest to infer the spatial structure of the baryonic matter in the universe. For the concordance \(\Lambda\)CDM model, the majority of matter in the universe is under the form of dark matter (DM) (Jarosik et al., 2011; Planck Collaboration VI, 2020).

Thus the baryonic spatial distribution is mainly dictated by the gravitational potential imprinted by the dark matter density. Therefore, Lyman-\(\alpha\) forest observations are good tracers of the total matter in the context of this cosmological model (Petitjean et al., 1995; Hui & Gnedin, 1997; Croft et al., 1998). Moreover, it is crucially important for cosmology, as at the relevant redshifts probed by Lyman-\(\alpha\), we do not have as many easily observable galaxies as at low redshift (\(z\lesssim 1\)). The Lyman-\(\alpha\) forest is thus a complementary probe of the matter clustering at \(z\simeq 2-3\), also known as the cosmic noon era.

The Lyman-\(\alpha\) forest can provide much insight into the cosmological model. Past work on the subject includes, for example, the recovering of the properties of the intergalactic medium using the Bayesian inversion method (Pichon et al., 2001), the constraining of the neutrino masses (Palanque-Delabrouille et al., 2015; Yeche et al., 2017), the inference of the three-dimensional matter distribution from 1D Lyman-\(\alpha\) absorption using Bayesian
Forward modelling (Porquares et al., 2019; Horowitz et al., 2019; Porquares et al., 2020; Kraljic et al., 2022), the constraints on the thermal history of IGM (Villasenor et al., 2022), and the measurement of the Baryonic Acoustic Oscillation (BAO), which provides a tight constraint on the expansion history of the universe (Font-Ribera et al., 2014). Several large cosmological surveys include nowadays a QSO survey with spectrum measurement (BOSS, eBOSS, DESI Lee et al., 2013; du Mas des Bourboux et al., 2020; Levi et al., 2019).

To analyse the Lyman-\(\alpha\) forest observations, we require hydrodynamical simulations of the IGM, which has to be computed alongside a dark matter \(N\)-body simulation. Though the IGM is coupled on large scales with the gravitational potential induced by dark matter overdensity, on small scales, it is dominated by baryonic shocks, cooling, and feedback, for which we do not have analytical solutions. Therefore, the IGM state has to be solved numerically. With upcoming cosmological surveys such as WEAVE-QSO (Pieri et al., 2016), the resolution and volume become higher and larger. They consequentially increase the required computational time drastically. Therefore using full hydrodynamical simulation to match the same volume is becoming intractable. Alternatively, \(N\)-body simulations of pure dark matter are much faster than hydrodynamical simulations. Thus, one could avoid solving for the hydrodynamic altogether and trying to correct phenomenologically the effect at small scales (\(\sim 100h^{-1}\) kpc, the Jeans' scale where Lyman-\(\alpha\) happens).

The fluctuating Gunn-Peterson approximation (FGPA, Gunn & Peterson, 1965; Weinberg et al., 1997) is a main framework to estimate a mock Lyman-\(\alpha\) from dark matter overdensity. However, this approximation does come with a few limitations. For example, it fails to capture the baryonic feedbacks at small scales, the power-law breaks down in the regions with higher density and strongly heated gas (Lukic et al., 2014; Kooistra et al., 2022), and the power law starts to decouple where \(z\gtrapproprox\). In most cases, the dark matter overdensity is smoothed before being used for emulation. For example, the baryonic pressure is emulated by convolving the matter density with a Gaussian kernel (Hui & Gnedin, 1997). Other work (Sorini et al., 2016) introduced a numerical simulation called Iteratively Matched Statistics (IMS) with an effort to improve the accuracy higher than Gaussian smoothing. Finally, the LyMAS-2 technique (Peirani et al., 2014; Peirani et al., 2022) relies on a conditional probability to map the flux from a smoothed dark matter field and introduces a Wiener filter to represent the coherence along line-of-sight better.

In recent years, the astrophysical community growingly adopted machine learning models because of their ability to provide a universal fitting function and its prediction speed through accelerators such as Graphics Processing Units (GPUs), which allows us to reach both a highly accurate and fast model, provided training data exist. Inspired by the work of Peirani et al. (2014), we propose to swap the conditional PDF with a neural network. In this work, we will test a _U-Net_ architecture, focusing on emulating the most essential fields to describe the IGM for Lyman-\(\alpha\). We introduce such an emulator, called _LyAl-Net_, trained to map a dark-matter density field from an \(N\)-body simulation to IGM fields (neutral hydrogen density and temperature) derived from a sibling hydrodynamical simulation. In this work, the two simulations are Horizon-DM and Horizon-noAGN at \(z=2.43\)(from the Horizon-AGN suite of simulations Dubois et al., 2014; Peirani et al., 2017). We then derive the Lyman-\(\alpha\) forest absorption features from the emulated fields. Previous work followed a similar initial strategy to derive emission and absorption features of the IGM (Villaescusa-Navarro et al., 2013; Harrington et al., 2021). In this work, we improve on their work by a larger set of tests of the accuracy of the models and allowing for flexibility to allow the models to be used in an observational context.

We also explored the method to re-calibrate the equation-of-state to improve the emulated Lyman-\(\alpha\) absorption for different gas physics, allowing the existing _LyAl-Net_ to be used in vastly different conditions, only paying with a small additional training set. This allows us to directly affect the quantity instead of simply adjusting the mean transmitted flux. This is a powerful tool for a fast simulation of Lyman-\(\alpha\) forest, with the potential of generalising a machine learning model to emulate different gas physics beyond the training set. This work aims to provide a framework and a proof of concept for generating more accurate emulators of the gas hydrodynamics and the Lyman-\(\alpha\) forest in particular.

This paper is organised as follows. We discuss the physical process and the model we use for the Lyman-\(\alpha\) absorption system in Section 2. In Section 3, we motivate and describe the architecture of the neural network used for training, including the data transformations and the hyperparameters for training. We present the prediction results and accuracy benchmarks using Horizon-noAGN and Horizon-DM in Section 4. In Section 5, we discuss the results using IllustrisTNG (Nelson et al., 2019) as a test set and the transfer learning method to fine-tune the equation of state. We discuss the results and conclude in Section 6.

## 2 Physical model of the Lyman-\(\alpha\) absorption

In this section, we give a brief reminder of the fundamentals of the physical mechanisms that produce the observed Lyman-\(\alpha\) forest. The details of the atomic physics underlying this derivation are given in earlier work by Meiksin (2009). We focus here on the relevant physical quantities we seek to model with our framework.

We can break down the mechanisms behind the observed Lyman-\(\alpha\) absorption features into a few components: the intrinsic linewidth of the resonance of the Lyman-\(\alpha\) transition, the thermal broadening owing to the small-scale random motions of the atoms, and the Doppler shift due to the bulk motions of the hydrogen clouds.

The absorption strength of the broadband Lyman-\(\alpha\) emission from one QSO allows us to infer the density of the neutral hydrogen at different distances along a line-of-sight. This absorption rate depends on the number of interactions. We can refer to it as opacity which is a function of the observed frequency \(\nu_{\rm obs}\) as

\[\tau(D_{\rm QSO},\hat{n},\nu_{\rm obs})=\int_{0}^{\nu_{\rm QSO}(D_{\rm QSO} )}n_{\rm HI}(s,\hat{n})\sigma_{\rm HI}(\nu_{\rm obs},s,\hat{n})\ {\rm d}s\,, \tag{1}\]

where \(n_{\rm HI}\) is the number density of neutral atomic hydrogen, \(\sigma_{\rm HI}\) is the interaction cross-section of photons with the HI atoms. We note that \(s\) is a physical proper distance. As such, \(s_{\rm QSO}(D_{\rm QSO})\) is the physical proper distance corresponding to the comoving distance \(D_{\rm QSO}\) of the QSO on the considered line-of-sight.

To derive the expression for all these quantities, we first consider the situation in the rest frame of some atomic hydrogen in space. The atomic cross-section can be described by the Lorentz profile with the following expression

\[\sigma_{\rm HI}(\nu)=\frac{\pi e^{2}}{m_{e}c}f_{\rm IL}L(\nu)=\frac{\pi e^{2}}{ m_{e}c}f_{\rm IL}\frac{\Gamma_{nl}/(4\pi^{2})}{(\nu-\nu_{\rm in})^{2}+( \Gamma_{nl}/(4\pi)^{2})}\,, \tag{2}\]

where \(L(\nu)\) is the Lorentz profile, \(\nu_{\rm in}\) is the frequency of the transition from lower (\(l\)) to upper (\(u\)) energy level, \(\Gamma_{nl}\) is the upper
energy level damping width, \(f_{\rm{J}u}\) is the oscillator strength, \(m_{e}\) and \(e\) are electron mass and its charge, respectively. We note that the transition levels for the Lyman-\(\alpha\) system are \(l=1\) and \(u=2\), which gives \(f_{12}=0.4162\).

In addition to the probabilities of the fundamental transition, we have to consider two other components. The HI atoms have finite positive temperatures that induce their thermal motions, which for non-relativistic gas is well described by a Maxwell distribution. This effect is called Doppler broadening, and we denote the distribution by \(G(\nu,T)\). This type of broadening depends on the temperature of the gas, which is itself position-dependent. As HI clouds are in motion with respect to the photons emitted by the QSO, we have to apply a further Doppler boost to put the incoming photon in the rest frame of the cloud.

By taking the redshift of the photon emitter into account, this gives the relation between the frequency in the two frames to be

\[\nu=\nu_{\rm{obs}}(1+z)\left(1+\frac{\nu_{z}}{c}\right) \tag{3}\]

where \(\nu\) is the central frequency, \(\nu_{\rm{obs}}\) the observed frequency,and \(\nu_{z}\) is the line-of-sight velocity. Therefore the cross-section in this scenario is the convolution of the Lorentz and Doppler profile,

\[\sigma_{\rm{HI}}(\nu,T)=\frac{\pi e^{2}}{m_{e}c}f_{12}(L*G)(\nu,T)=\frac{\pi e ^{2}}{m_{e}c}f_{12}V(\nu,T). \tag{4}\]

The convolution of these two profiles is called a Voigt profile, which we write as \(V(\nu,T)\):

\[V(\nu,T)=\int_{\nu=-\infty}^{\infty}{\rm d}\nu^{\prime}\frac{ \Gamma_{21}/(4\pi^{2})}{(\nu^{\prime}-\nu_{12})^{2}+(\Gamma_{21}/(4\pi)^{2})} \\ \frac{1}{\sqrt{\pi}\Delta\nu_{D}(T)}\exp\left[-\left(\frac{\nu- \nu^{\prime}}{\Delta\nu_{D}(T)}\right)^{2}\right]\,, \tag{5}\]

where

\[\Delta\nu_{D}=\frac{\nu_{12}}{c}\sqrt{\frac{2k_{B}T}{m_{H}}} \tag{6}\]

is the Doppler width, which is the frequency broadening of the Lyman-\(\alpha\) frequency due to the Doppler effect, and \(\nu_{12}\) is the Lyman-\(\alpha\) frequency.

We may rearrange the terms to simplify the expression by introducing \(x=(\nu-\nu^{\prime})/\Delta\nu_{D}(T)\), which is the frequency offset in the unit of \(\Delta\nu_{D}\). We set \(a=\Gamma_{21}/(4\pi\Delta\nu_{D})\) to be a ratio of the Lyman-\(\alpha\) line width to the Doppler frequency width \(\Delta\nu_{D}\). We can now rewrite the equation (5) to be:

\[V(\nu,T)=\] \[\int_{\nu^{\prime}}{\rm d}\nu^{\prime}\left(\frac{1}{\pi\Delta \nu_{D}}\right)\left(\frac{a}{(x(\nu)-x^{\prime})^{2}+a^{2}}\right)\frac{1}{ \sqrt{\pi}\Delta\nu_{D}}\exp\left(-x^{\prime 2}\right)\,. \tag{7}\]

The Voigt profile then reads

\[V(a,x)=\frac{a}{\pi^{3/2}\Delta\nu_{D}}\int_{-\infty}^{\infty}\frac{e^{-y^{2} }}{(x-y)^{2}+a^{2}}dy\,. \tag{8}\]

The normalisation of the Voigt profile is unity. We further transform the expression to express the Voigt profile \(V(a,x)\) in terms of the Voigt function \(H\) as

\[V(a,x)=\frac{1}{\sqrt{\pi}\Delta\nu_{D}}H(a,x)\,. \tag{9}\]

In the next section, we come back to the reason for using this function to increase computational efficiency. This convention allows us to use the amplitude of the absorption, which yields

\[H(a,x)=\frac{a}{\pi}\int_{-\infty}^{\infty}\frac{e^{-y^{2}}dy}{a^{2}+(x-y)^{2 }}\,. \tag{10}\]

Therefore the cross-section of the Lyman-\(\alpha\) absorption in the frame of the observer is,

\[\sigma_{\rm{HI}}(\nu_{\rm{obs}})=\frac{\pi e^{2}}{m_{e}c}f_{12}\frac{H(a,x)} {\sqrt{\pi}\Delta\nu}\,. \tag{11}\]

With the cross-section in Equation (11), we then compute the normalised absorption rate for Lyman-\(\alpha\), also named the Lyman-\(\alpha\) forest, as

\[F(\lambda)=\exp\left[-\tau\left(\nu_{\rm{obs}}\right)\right]. \tag{12}\]

Thus, we have an explicit algorithm for computing the absorption features encoded by \(F(\lambda)\). We detail the numerical implementation in Section 3.3. We note that we have not discussed the influence of the UV background, which directly affects the amount of neutral hydrogen (Becker & Bolton 2013), the relation between the spin temperature and the gas temperature (Liszt 2001), and possible relativistic effects (Irsic et al. 2016). Also, by considering the function \(F(\lambda)\), we assume that the Lyman-\(\alpha\) profile of the QSO is fully known since this is not the subject of this work. We only focus on how to emulate the quickest and the broadest possible effect of hydrodynamical simulation to obtain the function \(F(\lambda)\).

### Benchmark Model : modified-FGPA

To assess the performance of our model to generate Lyman-\(\alpha\) absorption, we employ a similar approach to the fluctuating Gunn-Peterson approximation (FGPA, Gunn & Peterson 1965; Weinberg et al. 1997) which is an analytical framework to estimate a mock Lyman-\(\alpha\) from dark matter overdensity based on the idea that dark matter gravitational potential attracts baryonic matter (Bi & Davidsen 1997). The approximation expresses as

\[\tau(z,{\bf r})\propto(1+\delta({\bf r}))^{g}, \tag{13}\]

where \(\tau\) is the IGM optical depth at redshift \(z\) at a comoving position \({\bf r}\), \(\delta\) is the dark matter overdensity, and \(\beta\) is derived from the slope of temperature-density power law relation, which relies on the equilibrium between HI collisional recombination and photoionisation.

We use a similar approximation for each individual fields, which yields the following relation for the neutral hydrogen density:

\[n_{\rm{HI}}{}^{\rm{mock}}=\bar{n}_{\rm{HI}}\left(\frac{\rho_{\rm{DM}}}{\bar{ \rho}_{\rm{DM}}}\right)^{\alpha_{\pi}}\,\,, \tag{14}\]

where \(\bar{\rho}_{\rm{DM}}\) is dark matter average density, \(n_{\rm{HI}}^{\rm{mock}}\) is an estimated neutral hydrogen number density, and \(\bar{n}_{\rm{HI}}\) is the coefficient of the power law.

And from the temperature-density power law relation, we can estimate a mock hydrogen temperature by using

\[T_{\rm{HI}}^{\rm{mock}}=\bar{T}_{\rm{HI}}\left(\frac{\rho_{\rm{DM}}}{\bar{ \rho}_{\rm{DM}}}\right)^{\alpha_{\cal{T}}} \tag{15}\]

where \(T^{\rm{mock}}\) is an estimated hydrogen temperature and \(\bar{T}_{\rm{HI}}\) is the coefficient of the power law. The parameters of the dark matter
from Horizon-DM fitting with Horizon-noAGN hydrogen density and temperature shows in Table 1. From here, we refer to the benchmark model as modified-FGPA (mod-FGPA). Furthermore, to compare modified-FGPA to _LyAl-Net_, we apply the estimated density and temperature into the Equation (1), which also includes the Voigt profile, and we assume that the line-of-sight velocity of dark matter and gas are approximately the same.

We also note that this model is not strictly speaking the FGPA, but it employs the same power-law approximation for the relation between the gas-state fields and the dark matter. We expect the result of this test to be more optimistic than classical FGPA but worse than the one obtained with _LyAl-Net_ because it is less flexible and it ignores the stochastic components of the gas.

## 3 Method

The model we use to calculate the absorption at different frequencies of the Lyman-\(\alpha\) broad emission line needs three ingredients that characterise the physical state of the neutral hydrogen: the temperature, the gas density, and the line-of-sight gas velocity. Therefore, the main objective of this work is to construct the neural network which can produce these parameters for the computation of the model of the Lyman-\(\alpha\) forest. We train _LyAl-Net_ for each gas parameter separately to accommodate the optimisation and the prediction interpretation. This section discuses the details of the cosmological simulation used as the training set, the numerical approach for generating the mock Lyman-\(\alpha\) absorption, and the data pre-processing procedure.

### The simulation datasets

The simulation that we used is the Horizon-AGN simulation. It is a cosmological hydrodynamic simulation based on the adaptive mesh refinement code Ramses (Teyssier, 2002). The simulation contains \(1024^{3}\) dark matter particles with a resolution of \(M_{\rm DM_{\rm{less}}}=8\times 10^{7}\)  M\({}_{\odot}\) and the volume span a cube of \(100h^{-1}\,\rm Mpc\), which took \(4\) million CPU hours to complete with a redshift down to \(z=1.2\). Horizon-AGN has been computed assuming a standard \(\Lambda\)CDM cosmology and adopting cosmological parameters compatible with the data of the WMAP mission (Komatsu et al., 2011): \(\Omega_{m}=0.272\), \(\Omega_{\Lambda}=0.728\), \(\sigma_{8}=0.81\), \(\Omega_{b}=0.045\), \(H_{0}=70.4\) \(\,{\rm{km}}^{-1}\,\rm Mpc^{-1}\), and \(n_{s}=0.967\). Furthermore, the simulation models the gas as ideal monatomic gas with adiabatic index \(\gamma=5/3\). The gas cools down via a metal-dependent model of hydrogen and helium photon emission (Sutherland & Dopita, 1993). This allows the temperature of the gas to cool down to \(10^{4}\) K. At the same time, the gas is heated via a uniform UV background following a model by Haardt & Madau (1996) after redshift \(z_{\rm{eion}}=10\). In this work, we used the Horizon-noAGN as a training set. It is the sibling simulation of Horizon-AGN with one major difference: the noAGN variant does not have AGN feedback due to the lack of black hole growth. The reason for using noAGN is that the computational time for the simulation was reduced by turning off the AGN feedback and emulating their effect by fine-tuning the UV background radiation to achieve the same gas properties (Peirani et al., 2022). Another upside of using Horizon-noAGN is that we can study the effects and sensitivities of the Lyman-\(\alpha\) forest model. We also used TNG100 simulation as a validation set to properly cross-check the results in Section 5.4.

### Dataset generation for training and validation

In this section, we describe the steps involved in creating the different datasets required for the training and validating the Machine Learning model.

The phase space distribution of dark matter is sampled with particles in the cosmological simulations that we considered. For our emulation strategy, we require the densities to be on a regular grid mesh to facilitate the computations. We have assigned the particles to a \(1024^{3}\) grid mesh, leading to a resolution of \(\sim 98h^{-1}\,\rm kpc\) for each voxel. Several mass assignment schemes exist, among which the most used is the Cloud-In-Cell assignment (Hockney & Eastwood, 1988). Unfortunately, this may lead to mesh cells completely devoid of mass and prone to very poor sampling in cosmic voids, which is exactly where the signal of the Lyman-\(\alpha\) forest is the most relevant (Porquerers et al., 2020). We rely on the particle assignment developed by Colombi et al. (2007) to avoid these problems. This assignment relies on an adaptive filter derived from Smooth Particle Hydrodynamics (SPH) kernel (Monaghan, 1992).

We briefly remind the reader of the procedure here. The complete classical description of a particle system is given by the

\begin{table}
\begin{tabular}{l c c} \hline \hline
**Gas Parameter (\(X\))** & \(\mathbf{\widehat{X}}\) & \(\alpha_{\mathbf{X}}\) \\ \hline Temperature & 3.98 & 0.58 \\ Density & -10.43 & 1.27 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Fitted parameters of the power-law relation of dark matter and IGM parameters in Equation (14) and (15). We use Horizon-DM dark matter overdensity and Horizon-noAGN gas density and temperature.

Figure 1: The flowchart illustrates a schematic of the _LyAl-Net_ pipeline for the Lyman Alpha forest simulation. The flux absorption model (Equation (12)) requires three gas physical quantities, namely temperature, density, and the line-of-sight velocity. LyNet emulates only the gas temperature and the gas density of neutral atomic hydrogen. We use dark matter line-of-sight velocity of the same _N_-body simulation as an approximated gas line-of-sight velocity.
phase-space distribution, \(f({\bf r},{\bf v}_{L})\). For our computations, we typically require the moment of that distribution, which is for the velocity:

\[{\bf v}({\bf r})=\frac{1}{\rho}\int{\rm d}^{3}{\bf v}_{L}\,{\bf v}_{L}f({\bf r}, {\bf v}_{L}), \tag{16}\]

where \(f({\bf r},{\bf v})\) is the phase space density of the dark matter particles. The mass density \(\rho({\bf r})\) is obtained similarly from the integral

\[\rho({\bf r})=\int{\rm d}^{3}{\bf v}_{L}f({\bf r},{\bf v}_{L}). \tag{17}\]

In practice, we compute Equations (16) and (17) into a mesh made of small cubic patch \(\Delta r^{3}\), and approximate the resulting fields piecewise. We do these computations using an approach similar to smooth particle hydrodynamics (Monaghan 1992) to reduce shot noise. We can represent each particle as a smooth cloud with a finite size based on their local density, which depends on the distances between neighbouring particles. Generally, we define the cloud size as \(2R_{\rm SPH}\) and the number of neighbouring particles as \(N_{\rm SPH}\). The algorithm we use preserves the conservation of mass and momentum by an adequate choice of mesh node weighing, as detailed below.

To interpolate particles into the grid site, we use the following equation

\[\tilde{A}(i,j,k)=\frac{1}{\left[R_{\rm SPH}(i,j,k)\right]^{3}}\left[\sum_{l=1 }^{N_{\rm S}-1}A_{l}W_{l}{\cal S}\left(\frac{d_{l}}{R_{\rm SPH}(i,j,k)}\right) \right]. \tag{18}\]

This equation provides a weighted sum of a quantity \(A_{l}\) carried by a particle \(l\). Each neighbouring particle is multiplied by \({\cal S}(x)\), an SPH kernel where \(x\) is the relative distance by calculating the ratio of the distance \(d_{l}\) of a particle \(l\) from the grid mesh node \((i,j,k)\) and \(R_{\rm SPH}(i,j,k)\). They are weighted using the \(W_{l}\) factor such that the total sum for each grid site is equal to unity. This constraint yields the following identity for the weight of each particle:

\[W_{l}=1/S_{l}=\left[\sum_{i,j,k}\frac{1}{\left[R_{\rm SPH}(i,j,k)\right]^{3}} {\cal S}\left(\frac{d_{l}}{R_{\rm SPH}(i,j,k)}\right)\right]^{-1}. \tag{19}\]

For the particular case of assigning the mass to the grid, we end up with a mass per mesh node of the grid. We divide by the volume corresponding to this mesh node to obtain the mass density. This yields:

\[\tilde{\rho}(i,j,k)=\tilde{m}(i,j,k)/\Delta r^{3}. \tag{20}\]

Similarly, the interpolated velocity is obtained following the equation (16) which gives the mesh the following identity:

\[\tilde{v}(i,j,k)=\tilde{p}(i,j,k)/\tilde{m}(i,j,k)\, \tag{21}\]

where \(\tilde{p}\) is the interpolated momentum of the cell.

Unlike dark matter, the gas simulation from Ramses has an adaptive mesh refinement structure. Each raw gas cells have a different comoving size. The largest cells have an extent of \(0.0977h^{-1}\,{\rm Mpc}\), while the smallest cells have 1/32th of this length. We homogenised the resolution by assigning the gas cell based on its coordinate to the nearest grid, similar to the nearest grid point scheme. Some of the mesh cells may contain multiple sub-cell, e.g. two \(L_{\rm largest}/4\) cells and one \(L_{\rm largest}/2\) cells. We used the same weight regardless of the length to average the quantities, which is a sufficient approximation given the scale we are working on with the Lyman-\(\alpha\) system. The algorithm can easily be expressed as

\[X_{cell}(i,j,k)=\frac{1}{N}\sum_{i=1}^{N}x_{i}, \tag{22}\]

where \(X_{cell}\) is the target parameter at coordinate \((i,j,k)\) and \(x_{i}\) is the sub-cell value of the target parameter.

### Mock Lyman-\(\alpha\) absorption generation

We have discussed the physics of the Lyman-\(\alpha\) forest earlier in Section 2. Let us now move to the numerical approach to calculate the normalised flux from the neutral hydrogen parameters. The opacity of the gas as a function of the observed frequency of a given line-of-sight defined in Equation (1) is in the continuous form. However, since the dataset is in a grid mesh structure, we have replaced Equation (1) with a Riemann sum to calculate the integral. The opacity now reads as

\[\tau_{a,k}=\sum_{j=0}^{j=N}n_{\rm HI,a,j}\sigma_{\rm HI}(y_{\rm obs,j},\mathbf{a},j)\times\delta l\, \tag{23}\]

where \(j\) is the index of the position with respect to the line-of-sight, \(l_{j}=j\times\delta l\) is the proper distance from an observer to a source, \(\mathbf{a}\) is the unit vector which is perpendicular to a source plane, \(\delta l\) is the physical length of each cell width, which in our case is \(\delta l=0.098h^{-1}\,{\rm Mpc}\). The Voigt function is expensive to compute numerically. We opted to rely on the Faddeeva function, which is precomputed as part of SciPy (Jones et al. 2001). We refer the reader to Appendix A for the relation between this function and the actual computation we do. We then compute the normalised absorption flux \(F(\nu)\) by following Equation (12). We do not consider the exact shape of the broadband emission Lyman-\(\alpha\) line of the mock QSO itself, which has to be estimated as usual in the observational context.

Converting the spatial resolution of the simulation to spectrometer resolution gives a resolution of \(R=\frac{\langle\Delta l\rangle}{\langle\Delta l\rangle}\simeq 30000\). This resolution is much higher than what is possible to achieve with current cosmological surveys. For example, WEAVE-QSO has \(R=5000\) and \(R=20000\) in low and high-resolution mode (Dalton et al. 2014), and SDSS IV (BOSS and eBOSS) has \(R=2000\). This allows us to have an error budget in the simulation compared to the real surveys. Of course, a high resolution does not mean the physics of the simulation is accurate, nor does it contains all of the baryonic physics at that resolution. We will consider the reliability of our model by comparing results obtained with other simulations in Section 5.

### Pre-transformation of the fields before training

The dark matter density field is usually described in terms of the density contrast defined as

\[\delta_{\rm DM}=\frac{\rho_{\rm DM}}{\tilde{\rho}_{\rm DM}}-1\, \tag{24}\]

where \(\tilde{\rho}_{\rm DM}\) is dark matter average density. The 1D distribution of \(\delta_{\rm DM}\) is heavily skewed with a range from -1 to infinity by construction. This is not necessarily the most convenient for emulators to manipulate. Furthermore, given enough flexibility, the skewed dataset does not affect the neural network. However, this transformation benefits several things: it compresses the dynamic range, allowing a scale-invariant process, and should help
the training speed and reduce the required sample size to achieve a given accuracy.

Therefore, we will provide the transformed density \(A_{\rho}\) derived from the density contrast as

\[A_{\rho}(\rho_{\rm DM})=\log_{10}\left(\frac{\rho_{\rm DM}}{\bar{\rho}_{\rm DM}} \right)=\log_{10}(1+\delta_{\rm DM}). \tag{25}\]

Similarly, we applied a decimal logarithmic transformation for the density and temperature of neutral hydrogen.

\[A_{\rm n,HI}=\log_{10}\frac{n_{\rm HI}}{1\ {\rm g\ cm^{-3}}}\, \tag{26}\]

and

\[A_{\rm T,HI}=\log_{10}\frac{T_{\rm HI}}{1\ {\rm K}}. \tag{27}\]

The choice of this transformation for density and temperature is not only for the coherence of the scale used between the input and prediction. This transformation should facilitate our interpretation of the differences in the dark matter distribution in the prediction analysis. Most importantly, following the suggestion that the relation between density and temperature is a power law in a diffuse IGM regime, which is typically located in the void regions (Hui & Gnedin, 1997; Schaye et al., 1999; Meiksin, 2009), this temperature-density relation expresses as

\[T=T_{0}\left(\frac{\rho_{b}}{\bar{\rho}_{b}}\right)^{\gamma-1}, \tag{28}\]

where \(\gamma\) is the adiabatic coefficient, \(T_{0}\) is the temperature at the mean density, \(\rho_{b}\) is the density of baryonic matter, and \(\bar{\rho}_{b}\) is the mean baryonic density. Since we are interested in neutral hydrogen, for simplicity, we assume \(\frac{\rho_{b}}{\bar{\rho}_{b}}\approx\frac{n_{\rm HI}}{n_{\rm HI}}\). Hence, a log transformation of this relation becomes a linear relation expressed as

\[A_{\rm T,HI}=A_{\rm T,HI}+\left(\gamma-1\right)\left(A_{\rm n,HI}-A_{\rm n,HI} \right). \tag{29}\]

We show a density-temperature plot in Figure 2, and we highlight the diffuse IGM region using the threshold of the temperature is \(T>10^{5}\) K and for the density is \(n_{\rm H}<10^{-4}(1+z)\ {\rm cm^{-3}}\)(Martizzi et al., 2019). This gas phase contains the majority of baryonic density, corresponding to \(~{}99.5~{}\%\) of the total volume and \(~{}76.7~{}\%\) of the total hydrogen mass. Atomic hydrogen density is also embedded within this field, therefore, carrying the most important information for a Lyman-\(\alpha\) forest. For dark matter and gas velocities, we do not re-scale for the neural network training since they are similar in magnitude. The hope is that a neural network can construct a phenomenological model that can relate corrections to this mean relation.

### LyAi-Net Architecture

This section discusses the neural network architecture used for this project and its advantage. Then we will discuss the training process, including the size of the training sample, the choice of the loss function, and the details of the architecture of _LyAi-Net_. For this work, we chose a U-net architecture for the neural network, firstly introduced by Ronneberger et al. (2015) for biomedical image segmentation. A U-net is based on the convolutional neural network (CNN, O'Shea & Nash, 2015), which can extract and create feature mapping similar to how the human visual cortex works. It is well known for its pattern recognition ability, spatial invariance, and scale invariance. This architecture is much faster than the traditional artificial neural network (ANNs, O'Shea & Nash, 2015), such as multilayer perceptron. It has become an invaluable technique for solving various problems, especially image-focus problems. The significant difference between CNNs and ANNs consists of the convolutional layers. These layers allow the machine to perform convolutions with kernels and leverage on three crucial points: sparse interactions, parameter sharing, and equivariant representations (Goodfellow et al., 2016). The U-net architecture allows us to extract and summarise the information on a small scale while retaining large-scale information thanks to network skipping. With the knowledge that hydrogen gas is the tracer of dark matter and vice versa, the U-net becomes a practical approach to generating the gas fields in a short amount of time compared to the typical hydrodynamic simulation, which also has been used by Wadekar et al. (2020) and Bernardini et al. (2021).

The schematic of the _LyAi-Net_ is shown in Figure 3. It comprises three major steps: the convolution side (on the left), the bottleneck (at the bottom), and the up-sampling side (on the right). We chose the \(3\times 3\times 3\) convolution kernel for each convolution step while maintaining the dimension using 'SAME' padding and adopting a dropout layer with a rate value of \(0.2\).1. The \(3\times 3\times 3\) max-pooling is also used to reduce the dimension for each contracting step.2 After the bottleneck phase, we used the \(3\times 3\times 3\) up-sampling layer followed by the same configuration of convolution. For each convolution layer of the up-sampling phase, we used 16 filters applied on the input concatenation and the bypass from the contracting side.

Footnote 1: A convolution kernel with ’SAME’ padding adds zeros around the input, such as the output as the same dimensions as the input.

Footnote 2: Max-pooling operates by extracting the maximum value within the \(3\times 3\times 3\) kernel, scanning in the tensor.

The overall dimension of the input is \(81\times 81\times 81\) pixels, and the output is \(27\times 27\times 27\) pixels, which is cropped to prevent the possible edge effect caused by the convolution. Thus, we need to tile the output into the larger cube to obtain the full simulation from the output, which will be explained later in Section 3.5.3.

Figure 2: A Density-Temperature (\(n-T\)) diagram of hydrogen density vs temperature in the Horizon-noAGN simulation at \(z\approx 2.4\). The black dashed line illustrates the fitted strong power-law relation within the diffused IGM regime, which is the main contribution of Lyman-\(\alpha\) forest, corresponding to \(~{}99.5~{}\%\) of the total volume and \(~{}76.7~{}\%\) of the total hydrogen mass. Light grey dashed line separates the different baryonic regimes. We provide the density plot in this phase diagram by colouring each point with a log scale on the right-hand side.
We note that this architecture is not optimal, and a better architecture can exist to provide an equivalent or better result with the smaller training parameters.

#### 3.5.1 Loss Function

We choose a mean squared error (MSE) as the loss function. The main advantage of using this loss function is that it is more stable compared to GAN loss (Nguyen et al., 2021). It also allows a better understanding of what the neural network learns, therefore, better interpretability. It takes a form as

\[\text{MSE}\left(Y_{i},\hat{Y}_{i}\right)=\frac{1}{n}\sum_{i=1}^{n}(Y_{i}-\hat{ Y}_{i})^{2}\, \tag{30}\]

where \(n\) is the batch size, \(Y_{i}\) is the prediction, and \(\hat{Y}_{i}\) is the true value.

We note that this loss is applied to the transformed variables (\(A_{X}\)), and it only calculates the loss of the central part with \(27\times 27\times 27\) pixels. For the velocity field, this stays linear. For the other quantities, as \(T_{\text{HI}}\) and \(n_{\text{HI}}\), owing to the log transformation, the effective loss function is the mean squared log error (MSLE) which reads

\[\text{MSLE}(Y_{i},\hat{Y}_{i})=\frac{1}{n}\sum_{i=1}^{n}\left(\log\frac{Y_{i} }{\hat{Y}_{i}}\right)^{2}. \tag{31}\]

This loss weighs high and low values equally and eases the training process. The MSLE loss punishes the underestimation more than the overestimation. It is a two-sided choice, both carrying a risk and an advantage. Indeed it can lead to bias in the mean response, as we will investigate later. Doing so also avoids being overly sensitive to the high-value excursion, which is irrelevant for the Lyman-\(\alpha\) forest. As long as the prediction is in a sensible window of error and the baryonic observation's resolution, this loss function should not be the primary concern. We will discuss prediction bias and error later in Section 4.

#### 3.5.2 Training Process

We train _LyAl-Net_ to predict each gas parameter separately because we can tweak the training hyper-parameters and interpret the prediction results independently from the other fields. We used Horizon-noAGN simulation as the training set, which contains both dark matter and gas components. It spans a volume of \((100h^{-1}\,\text{Mpc})^{3}\), with a resolution of 1024 pixels on each side. This snapshot has the redshift range from \(z=2.44\) to \(z=2.33\). Since we were limited to a single full simulation, we randomly selected 5120 sub-boxes, with each sub-box having a size of \(81^{3}\) voxels or \(\sim(7.9\ h^{-1}\,\text{Mpc})^{3}\) in volume, which is equivalent to \(\sim 9\%\) of the total simulation volume. The data sampling also allows us to use repeated samples and overlapping of the training set, which encourages the network to learn translational invariance. The size of the sub-box is intended to fit the GPU memory and encourage the neural network to learn the information at a small scale. We note that a physically larger kernel size could improve the solution by better capturing the transfer of power from large to small scales. We will discuss this further in Section 6.

Figure 3: We show here the schematic of the _LyAl-Net_ architecture. The input of the network is the transformed dark matter field, \(\log_{10}\left(\delta_{DM}+1\right)\), with the size of \(81\times 81\times 81\) voxels. The output is the targeted gas parameter (temperature or gas density) with the size of \(27\times 27\times 27\) pixels. We crop the volume at its centre from the last convolutional layer. The arrow represents the direction of the data where each colour represents a different operation with \(3\times 3\times 3\) kernel size. We note that all of the convolution layers extract 16 feature maps. The model is trained to predict each gas parameter separately for the simplicity of optimisation and interpretation.
We also perform rotational augmentation 3 in all three axes, except the line-of-sight velocity which is rotated solely around the z-axis. The rotation is required for the trained neural network to obtain the rotational invariance (He et al., 2019; Kaushal et al., 2022; Ramanah et al., 2019), reflecting the physical process involved in the prediction. We expect this process to make the emulator robust. This property will be useful for transfer learning for a different dark matter simulation.

Footnote 3: To increase the diversity of the dataset from the existing dataset

Table 2 shows summarised details of the parameters of the training loop that we used, which includes the number of epochs and the learning rate. We used a batch size of 64 with Adam Optimiser (Kingma and Ba, 2014). We note that the hyper-parameters of the _LyAl-Net_ were obtained via numerical experimentation to ensure a repeatable and quick convergence.

#### 3.5.3 Tiling Algorithm

As mentioned, the simulation size is discretised into \(1024^{3}\) voxels. Since _LyAl-Net_ produces emulated fields of only \(27^{3}\) voxels from an \(81^{3}\) voxels input, we have to tile emulated elementary cubes to achieve the desired volume before assessing the accuracy and quality of the prediction model. We implemented a tiling algorithm by applying a 27 pixels sliding window, horizontally scanning from the XY-plane 38 times per row, column, and layer (Z direction). We describe with pseudo-code the algorithm in Algorithm 1. The total total volume consists of \(38\times 38\times 38\) elementary cubes (equivalent to \(1026^{3}\) voxels), which is then trimmed to obtain the \(1024^{3}\) voxels. The construction of the complete volume took \(\sim 8\) minutes on Nvidia Tesla V100 and Intel(R) Xeon(R) Gold 6230, 40 cores. This cost can be further reduced by re-arranging some operations directly on the GPU to alleviate GPU-to-CPU data transfer. We postpone this optimisation for later, which could be useful for a fully differentiable model of the Lyman-\(\alpha\) forest.

## 4 Results on Horizon-noAGN

In the following, we present the metrics we use to assess the quality of the emulated fields in Section 4.1. Then we consider in turn each of the physical fields in the following section: the number density of atomic neutral hydrogen in Section 4.2, the temperature in Section 4.3 and the gas bulk velocity in Section 4.4. We emphasise that we train _LyAl-Net_ with Horizon-noAGN, and use Horizon-DM as a validation set to see how the model responds to different baryonic feedbacks. We also refer to the term "predicted fields" and "emulated fields" interchangeably.

### Prediction Quality Assessment

We consider several metrics to assess the degree of trust that we have on _LyAl-Net_ to predict physical fields. We will first consider the one-point statistics, i.e. the relation between the true and the predicted field at a given point in space. The second statistic that provide insight into the reproduction of the mean spatial structure of the two respective fields.

#### 4.1.1 One-point Statistics

We consider several metrics to assess the degree of trust we may have on _LyAl-Net_ to predict physical fields. We will first consider the one-point statistics, i.e. the relation between the true and the predicted field at a given point in space. The second statistic that we will consider is the two-point correlation function which will provide insight into the reproduction of the mean spatial structure of the two respective fields.

#### 4.1.2 One-point Statistics

We quantify the emulator's performance by considering the conditional probability distribution function of the predicted value \(X_{\text{Pred}}\) of the field \(X\) given the true value \(X_{\text{true}}\) of that same \(X\) from the simulation. This conditional probability is obtained from Bayes' theorem directly, which is expressed as

\[P(X_{\text{Pred}}|X_{\text{True}})=\frac{P(X_{\text{Pred}},X_{\text{True}}) }{P(X_{\text{True}})}, \tag{32}\]

where \(P(X_{\text{Pred}}|X_{\text{True}})\) is the probability of the predicted value of \(X\) given the true value of \(X=X_{\text{Pred}}\), \(P(X_{\text{Pred}},X_{\text{True}})\) is the joint probability of \(X=X_{\text{Pred}}\) and \(X=X_{\text{True}}\), and \(P(X_{\text{True}})\) the probability of \(X_{\text{True}}\).

We have estimated the joint probability distribution function for the predicted and true values \(P(X_{\text{Pred}},X_{\text{True}})\) and \(P(X_{\text{True}})\) by using the kernel density estimation (KDE, Silverman, 1986; Scott, 2015). It is a technique to infer a continuous probability density function (PDF) from a finite sample by smoothing the distribution using a weighted kernel function. Owing to the computational cost of this estimator, we built them on a reduced number of voxels. We, therefore, use approximately 50% of the total volume to calculate the KDE estimator. The full details are indicated in Appendix B.

#### 4.1.3 Two-point Statistics

The two-point correlation function \(\xi(r)\) is defined as

\[\xi(|\mathbf{r}|)=\langle\delta_{A}(\mathbf{r}^{\prime})\delta_{B}(\mathbf{r} ^{\prime}+\mathbf{r})\rangle. \tag{33}\]

Compared to a homogeneous distribution, this measures the excess probability of finding two-point objects, field fluctuations, separated by a distance \(|\mathbf{r}|\). These objects or fields may be,for example, galaxies or gas quantities. The power spectrum is the image of the correlation function in the spatial frequency domain. It is defined as

\[P(|\mathbf{k}|)=\int\mathrm{d}^{3}\mathbf{r}\ \xi(r)e^{-i\mathbf{k}\cdot \mathbf{r}}. \tag{34}\]

The metrics that we used to assess both magnitudes and spatial information of the emulated fields are the transfer function \(T_{X}(k)\) and the cross-correlation function \(r_{X}(k)\).

The transfer function \(T_{X}(k)\), \(k\) being the spatial comoving wave numbers, is a convenient way to represent the departure of the predicted field to the original correlation structure of the true field \(X\) for each \(k\) scale. This function is commonly used in literature as a benchmark of the performance of an emulator or

\begin{table}
\begin{tabular}{l c c} \hline \hline
**Target HI Parameter** & **Learning Rate** & **Number of Epochs** \\ \hline Density & 0.002 & 1000 \\ Temperature & 0.002 & 1000 \\ Velocity & 0.003 & 100 \\ \hline \hline \end{tabular}
\end{table}
Table 2: The table summarises the hyper-parameters used for _LyAl-Net_ training with Horizon-noAGN.
an approximation (e.g. Bardeen et al., 1986; Leclercq et al., 2013; Vlah et al., 2016; He et al., 2019; Dai & Seljak, 2021). The transfer function is defined as

\[T_{X}(k)=\sqrt{\frac{P_{\rm pred,{X}}(k)}{P_{\rm true,{X}}(k)}}\, \tag{35}\]

where \(P_{\rm pred}(k)\) is the power spectrum of the emulated field, and \(P_{\rm true}(k)\) is the reference or true field.

Finally, we consider the correlation rate \(r_{X}(k)\) to check the linear correspondence

\[r_{X}(k)=\frac{P_{\rm true,{X}\cup{\rm pred,{X}}(k)}}{\sqrt{P_{\rm true,{X}}( k)P_{\rm pred,{X}}(k)}}. \tag{36}\]

To compute the power spectrum of the different scalar fields considered in this work, we rely on the use of nbodykit(Hand et al., 2018). For the prediction quality assessment of Horizon-noAGN, we decide to include the trained portions of the dataset in the power spectra calculation because only \(\leq 9\%\)4 of the simulation volume was used to train _LyAI-Net_. The main reason is to avoid an irregular volume for power spectra computation, and the volume should be small enough compared to the rest of the validation volume.

Footnote 4: The data sampling algorithm for a training set allows repeated sampling.

### Neutral hydrogen Density (\(n_{\rm HI}\))

Among the different quantities required to compute the absorption rate in the Lyman-\(\alpha\) forest, the atomic hydrogen number density is most critical to obtain correctly since the neutral hydrogen clouds are the absorber in the Lyman-\(\alpha\) forest process. In the top row, Figure 4 shows a slice of the atomic hydrogen number density \(n_{\rm HI}\) from the Horizon-noAGN prediction from dark matter overdensity and the (logarithmic) difference between the ground truth and prediction. Visually, the structures of both slices are in good agreement. We note no apparent defects caused by _LyAI-Net_ nor edge effects induced by our tiling procedure to cover the entire volume except for spots where a high relative error is present, typically correlated with high-density regions and where Lyman-\(\alpha\) absorptions are saturated.

We now follow the model evaluation procedure in Section 4.1 to compute the true number density's joint and conditional probability distribution from Horizon-noAGN and the emulated number density. The estimated probability distribution function \(P(\log_{10}n_{\rm HI}^{\rm Prod}\log_{10}n_{\rm HI}^{\rm True})\) is shown in Figure 5 with contour lines representing the different levels of probability. The green dashed line, referred to as \(\mu\left(\log_{10}n_{\rm HI}{}^{\rm Pred}|n_{\rm HI}^{\rm True}\right)\) is the mean predicted neutral hydrogen density from the conditional probability, along with the \(68\%\) confidence interval represented in the green band. This is modelled after the fact that the conditional probability is quite close to a Gaussian distribution. The red dashed line is given as a reference where \(n_{\rm HI}^{\rm Pred}=n_{\rm HI}^{\rm True}\). We further added a red band representing the upper and lower boundary of \(10\%\) error from the mid-point as a qualitative reference to detect the mean bias of the emulator, which we call the prediction bias. The contour levels and the mean imply that the prediction behaves well within the main region of the dataset, given that the \(0\)th and \(99\)th percentiles are \(4.5\times 10^{-9}\) and \(3.8\times 10^{-13}\,{\rm cm}^{-3}\) respectively. The bias only starts to occur when \(n_{\rm HI}<10^{-12}\,{\rm cm}^{-3}\) and \(n_{\rm HI}>3\times 10^{-9}\,{\rm cm}^{-3}\) which might be caused by the longer

Figure 4: Comparisons of sample slices for different hydrodynamic quantities in the decimal logarithm scale. Left column: fields directly extracted from Horizon-noAGN. Middle column: fields obtained from the application of _LyAI-Net_ to the dark matter field. Right column: the difference between the leftmost panel and the middle panel. The atomic hydrogen number density (temperature, respectively) is presented on the top row (bottom row, respectively).
distribution tail and sample variance affecting the KDE itself. To get a better picture of the prediction bias, we define the bias as a function of true value as

\[\text{Bias}=\left[10^{n\left(\log_{10}\lambda_{\text{HI}}^{\text{Mean}}\log_{10} \lambda_{\text{HI}}^{\text{True}}\right)}-\lambda_{\text{HI}}^{\text{True}} \right]\frac{1}{\lambda_{\text{HI}}^{\text{True}}}, \tag{37}\]

where \(X_{\text{HI}}\) is a gas parameter. We visualise this function for the emulated \(n_{\text{HI}}\) in Figure 6, and it shows the prediction bias where \(10^{-12}\,\text{cm}^{-3}<n_{\text{HI}}<10^{-10}\,\text{cm}^{-3}\) stay within 10% fiducial bias. By comparing _LyAl-Net_ to the benchmark, the prediction from modified-FGPA is unbiased only around the density \(n_{\text{HI}}\)\(\approx 10^{-11}\,\text{cm}^{-3}\), and fail to capture the broader range of the density.

While one-point statistics show that _LyAl-Net_ performs well, we extend our analysis to the two-point statistics of the emulated fields. Before calculating the power spectra, we mask the extreme values for emulated and true \(n_{\text{HI}}\) to zero. We based this on the fact that the bias is starting to occur in the tails of the distribution due to the fewer data points. Furthermore, masking has to be done because power spectra are very sensitive to extreme values (see Figure C.1). Therefore, we implement three different masking ranges and refer to them as _fidelity range_: low, medium, and much. Table 3 summarises the threshold values used in this process. These masking values are arbitrary, so we can identify the regions in which _LyAl-Net_ performs best. We emphasise that the maximum number of masked voxels is less than 1%. Therefore, this should not artificially increase the correlation rate \(r(k)\).

We compare the transfer function and cross-correlation using different fidelity ranges illustrated in Figure 7. The transfer function is approximately stable for all fidelity ranges while \(T_{\text{HI}}(k)\) of high fidelity performs the best and stays around 0.9 up to \(k\approx 11h\,\text{ Mpc}^{-1}\). The cross-correlation functions \(r_{\text{shift}}\) show the convergence of all fidelity ranges with the amplitude staying above 0.9 up to \(k\approx 3h\,\text{ Mpc}^{-1}\). Based on two-point statistics, this assessment shows that the best performance of the _LyAl-Net_ is up to \(n_{\text{HI}}\leq 1.10\times 10^{-9}\,\text{ cm}^{3}\), which is around the estimated diffuse IGM limit. 5

Footnote 5: We estimated the diffuse IGM limit of neutral hydrogen density by scaling the Martzzi et al. (2019)\(n_{\text{H}}\) limit; \(n_{\text{HI}}=\alpha n_{\text{H}}=10^{-4}(1+z)\,\text{ cm}^{-3}\) where \(\alpha=\text{mean}(n_{\text{HI}}/n_{\text{H}})\), obtained directly from Horizon-noAGN simulation at \(z=2.43\).

### HI gas temperature (\(T_{\text{HI}}\))

We perform the same analysis for the HI gas temperature as the density. We show the bottom panel in Figure 4 a slice of through the temperature field in the Horizon-noAGN simulation, the emulated temperature from dark matter density and the logarithmic difference between those two fields. The emulated and true tem

Figure 5: The conditional probability of the predicted decimal logarithm of density \(n_{\text{HI}}\) given the true density of the Horizon-noAGN and their corresponding marginal distribution where the black solid contour lines represent different levels of probabilities (see Equation (32)). The diagonal red-dashed line indicates the unbiased relation between the predicted and true density alongside a 10% fiducial error budget. The green dashed green line indicates the mean of the conditional probability distribution alongside the \(67\%\) confidence interval. We note a good agreement over several orders of magnitude of the predicted gas density from \(1\times 10^{-12}\,\text{cm}^{-3}\) and \(6\times 10^{-10}\,\text{cm}^{-3}\). A significant saturation occurs for very low density at \(n_{\text{HI}}\leq 10^{-12}\,\text{cm}^{-3}\).

\begin{table}
\begin{tabular}{c c c} \hline \hline \multirow{2}{*}{**Fidelity Range**} & \multicolumn{2}{c}{**Upper Limit Masking Value**} \\ \cline{2-3}  & **Density**\([\text{cm}^{-3}]\) & **Temperature**\([\text{K}]\) \\ \hline low & \(5.16\times 10^{-7}\) & \(3.08\times 10^{6}\) \\ medium & \(1.79\times 10^{-8}\) & \(5.92\times 10^{5}\) \\ high & \(1.10\times 10^{-9}\) & \(8.03\times 10^{4}\) \\ \hline \hline \end{tabular}
\end{table}
Table 3: The different density ranges used in the masking process to evaluate the _LyAl-Net_ performance where the upper-density boundaries are based on 99.95th, 99.99th,and 99.5th percentile for low,medium,and much fidelity range respectively. The lower boundary of neutral hydrogen density is \(4.36\times 10^{-15}\text{cm}^{-3}\) and for temperature is 600 K. The maximum number of masked voxels is less than 0.5%, therefore, this should not artificially increase the cross-correlation \(r(k)\).

Figure 6: The emulated \(n_{\text{HI}}\) bias as a function of the true \(n_{\text{HI}}\) of Horizon-noAGN at \(z=2.44\) of _LyAl-Net_ (an orange solid line) and the modified-FGPA (a blue solid line). The bias of modified-FGPA is only unbiased at \(n_{\text{HI}}\approx 10^{-11}\,\text{cm}^{-3}\). The _LyAl-Net_’s bias shows a better prediction overall. It is mostly stable and stays inside the 10 % fiducial bias levels for density range between 1st and 90th percentile, highlighted in a light blue strip which corresponds to \(1\times 10^{-12}\) and \(6\times 10^{-11}\,\text{cm}^{-3}\) respectively. The drastic behaviour of the bias outside the blue strip is expected since the extremities of the density on both sides rarely occur.
perature also agree without any apparent defects caused by the _LyAl-Net_. The logarithmic difference map of temperature shows that the map is, on average better than density, with a lower relative error for dense regions.

As for the density section, we also compute the emulated temperature's joint and conditional probability distribution given the true temperature, \(P(\log_{10}T_{\rm HI}^{\rm pred}|\log_{10}T_{\rm HI}^{\rm True})\) illustrated in Figure 8, and we follow the same colour scheme as \(n_{\rm HI}\). The contour lines represent different probability levels, the green dashed line is the mean predicted temperature, and the green band is the 68% confidence interval taken from its maximum. The red dashed line is given as a reference where \(T_{\rm HI}^{\rm pred}=T_{\rm HI}^{\rm True}\) with a red band representing the 10% fiducial error from the midpoint as a qualitative reference. Visually, the mean values of the emulated temperature appear to have a small bias across the temperature range, which becomes more prominent in the temperature regime. The details of the emulated temperature bias as a function of the true values are illustrated in Figure 9 in a solid orange line. Thus, we consider our emulator to be nearly unbiased and stable in the range of \(2\times 10^{3}~{}{\rm K}<T_{\rm HI}<10^{4}~{}{\rm K}\). The overprediction only appears when for the extreme temperature above \(10^{4}~{}{\rm K}\). In addition, _LyAl-Net_ performs impressively better than modified-FGPA, which suffers from prediction bias in all ranges of the temperature value.

To estimate the fidelity range of _LyAl-Net_, we also picked three masking ranges, summarised in Table 3. From Figure 10, the transfer function indicates the amplitude stays well greater than 0.90 up to \(k\approx 10~{}h~{}{\rm Mpc}^{-1}\) for the most fidelity ranges, with a small drop at \(k<0.3h~{}{\rm Mpc}^{-1}\). The spatial information measured by the cross-correlation function indicates \(r>0.9\) up to \(k=2h~{}{\rm Mpc}^{-1}\) for the low and med fidelity range and up to \(k=3h~{}{\rm Mpc}^{-1}\) in the high fidelity range. This implies the highest performance of the _LyAl-Net_ lies in the scale of \(14h^{-1}~{}{\rm Mpc}\) to \(63h^{-1}~{}{\rm Mpc}\) scale for the temperature range \(1.92\times 10^{3}<T_{\rm HI}<8.03\times 10^{4}~{}{\rm K}\). It is worth pointing out that the mean prediction bias is stable when the values are within the 1st and 90th percentile on both \(n_{\rm HI}\) and \(T_{\rm HI}\). This range is empirical evidence that _LyAl-Net_ performs best in the main distribution of both fields, while the high bias only occurs in the rare and extreme environment in the tails of the distributions. Thus, this implies an intrinsic limitation of the neural network weighing the temper

Figure 8: The conditional probability of the predicted decimal logarithm of temperature given the true temperature of the Horizon-noAGN where the contour lines represent different probabilities. The diagonal red-dashed line indicates the unbiased relation between those two temperatures alongside a fiducial 10% error budget in shaded red. The dashed green line indicates the mean from the conditional distribution alongside the 67% confidence interval (in shaded green).

Figure 7: Horizon-noAGN prediction quality assessment using two-point correlations for a different choice of fidelity range for \(n_{\rm HI}\). **Top panel:** We show the transfer functions \(T_{n_{\rm HI}}(k)\) for the different for low (blue solid line), and (orange solid line), and (green solid line). The descriptions are given in Table 3. **Bottom panel:** We show the correlation rate \(r_{n_{\rm HI}}(k)\) for the same fidelity ranges as in the top panel, with the same colour convention.

Figure 9: The emulated \(T_{\rm HI}\) bias as a function of the true \(T_{\rm HI}\) of Horizon-noAGN at \(z=2.44\) of _LyAl-Net_ (in an orange solid line) and the modified-FGPA (a blue solid line). The modified-FGPA bias lies outside the 10% fiducial bias levels for all temperature range. While the bias of _LyAl-Net_ is mostly stable and stays inside the 10 % fiducial bias levels for temperature values between 1st and 90th percentile highlighted in a light blue strip which corresponds to \(2\times 10^{3}~{}{\rm K}\) to \(2\times 10^{4}~{}{\rm K}\).
ature and density in the tails of distribution less due to the lack of data points. Fortunately, in practice, the Lyman-\(\alpha\) absorption is not sensitive to the highly dense region because the saturation occurs at a density lower than where the discrepancy happens (see Equation 12).

### The line-of-sight velocity of the HI gas (\(v_{z,\mathrm{HI}}\))

The gas line-of-sight velocity, \(v_{z,\mathrm{HI}}\), plays a significant role in the Lyman-\(\alpha\) absorption by shifting the frequency of an incoming photon due to the Doppler effect (for a non-relativistic case). With inaccurate values of the velocity field, the emulated flux will result in a change in absorption rate and a shift of Lyman-\(\alpha\) absorption in some frequencies. To estimate \(v_{z,\mathrm{HI}}\), we can use the line-of-sight velocity of dark matter \(v_{z,\mathrm{DM}}\) directly, because their velocity fields are nearly identical since the gravity from dark matter influences the IGM. This estimation is effective on large scales, but it is inevitable for the IGM velocity to decouple at small scales due to various shocks and feedback (Pando et al., 2004). Therefore, naively using \(v_{z,\mathrm{DM}}\) could potentially impact the accuracy of absorption at small scales, which means there is room for improvements to emulate \(v_{\mathrm{HI}}\). To show the decoupling effect empirically, we compare velocity divergences of \(V_{\mathrm{DM}}\) to \(V_{\mathrm{HI}}\) using the transfer function and cross-correlation in Figure 11 in a solid blue line. The transfer function \(T_{\nabla_{\nu}}(k)\) of dark matter shows a steady amplitude \(\approx 1.1\) where \(k<2h\) Mpc\({}^{-1}\), which means the dark matter velocity fluctuation is higher than the gas. The cross-correlation rate \(r_{\nu,\nu}(k)\) on the bottom panel shows the value drops below 0.95 when \(k>2h\) Mpc\({}^{-1}\), meaning the IGM velocity starts to decouple with below \(3\,h^{-1}\) Mpc. Based on the idea that _LyAl-Net_ should be able to correct the small-scale effect, we train _LyAl-Net_ to predict a HI line-of-sight velocity using dark matter line-of-sight velocity as an input. We check whether _LyAl-Net_ improves the results instead of simply adding more systematic effects to the dark matter velocity field and illustrate in a solid orange line in Figure 11. Both methods show similar performance in the spatial distribution of cross-correlation of velocity divergences. However, the transfer function of dark matter shows better reliability in the high-\(k\) regime, while _LyAl-Net_ introduces a loss of amplitude over the same range. Despite this, we use \(v_{z,\mathrm{DM}}\) as a proxy for \(v_{z,\mathrm{HI}}\) and postpone the improvement for velocity prediction in future work.

### The effect from an absence of the baryonic feedback

With the additional complexity of the equations and the large increase in the number of timesteps, cosmological hydrodynamic simulations are largely more computationally expensive by at least an order of magnitude (Periani et al., 2017) than an \(N\)-body simulation. While it is cheaper to produce dark matter-only fields, neglecting the baryonic feedback can reduce the accuracy of the emulated gas fields on smaller scales. Since the main goal is to utilise _LyAl-Net_ on pure dark matter \(N\)-body simulations for time efficiency, in this section, we explore how a pure dark matter \(N\)-body simulation can affect the quality of emulated fields. We then look at the impact of this approximation on the Lyman-\(\alpha\) forest absorption in Section 4.6. For this test, we emulate HI density and temperature from Horizon-DM. This simulation is a sibling simulation of Horizon-AGN, which considers only pure dark matter dynamics. We use the same snapshot, i.e. the same redshift range, as the one used for the training of _LyAl-Net_. Since the neural network model is trained on Horizon-noAGN with the dark matter density transformation described by Equation (25), it allows _LyAl-Net_ to be immune to the difference in dark matter particle mass of pure \(N\)-body and a full hydrodynamical simulation. We can use these emulated fields as benchmarks for the case of emulated fields from Horizon-DM directly.

Figure 11: Horizon-noAGN prediction quality assessment using two-point correlations of velocity divergence from Horizon-noAGN dark matter velocity and emulated Horizon-noAGN HI velocity, labelled as **Dark Matter** in blue and _LyAl-Net_ in orange respectively. **Top Panel:** Transfer functions comparison shows both have a similar amplitude. However, Horizon-noAGN dark matter velocity divergence outperforms where \(k\gtrsim 6h\) Mpc\({}^{-1}\). **Bottom Panel:** Cross-correlation functions comparison shows both are identical, which implies that _LyAl-Net_ does not improve the spatial distribution.

Figure 10: Horizon-noAGN prediction quality assessment using two-point correlations for different choices of fidelity range for emulated \(T_{\mathrm{HI}}\). **Top Panel:** The transfer functions comparison for low, med, and high fidelity range in blue, orange and green respectively. The temperature range is indicated in Table 3. **Bottom Panel:** A comparison of correlation rates of the same fidelity ranges with the same colour scheme.
Let us introduce a new metric, which is the ratio of transfer functions,

\[\frac{T_{x,\mathrm{H-DM}}(k)}{T_{x,\mathrm{H-noAGN}}(k)}=\sqrt{\frac{P_{x,\mathrm{ H-DM}}(k)}{P_{x,\mathrm{H-noAGN}}(k)}}, \tag{38}\]

where the a subscripts \(x,\mathrm{H-noAGN}\) and \(x,\mathrm{H-DM}\) are denoted as a gas field \(x\) emulated using Horizon-noAGN and Horizon-DM respectively. We also use a high fidelity masking range introduced in Table 3 for both fields.

The top panel of Figure 12 shows the ratio of the transfer functions for the emulated density and the temperature. The ratio is relatively stable across the whole range in the considered scales for the neutral hydrogen density. However, this stable value stays at \(\approx 0.85\), which means the absence of baryonic feedback impacts the fluctuation of \(n_{\mathrm{HI}}\). The cross-correlation function of \(n_{\mathrm{HI}}\) (middle panel) shows a slight increase in its accuracy at \(k\approx 5h\) Mpc\({}^{-1}\). This is slightly better than the Horizon-noAGN and contradicts the expectation that the emulated fields from Horizon-DM should be worse due to the lack of feedbacks. On the other hand, the emulated temperature from Horizon-DM performs much closer to using Horizon-noAGN because the ratio stays relatively stable at \(0.95\). Its cross-correlation function shows that using Horizon-DM slightly impacts the spatial distribution of the emulated field on a large scale. It is worth pointing out that when \(k>10h\) Mpc\({}^{-1}\), the transfer function ratio starts to deviate for both fields, which implies that the lack of baryonic feedback dominantly affects scales above \(10h\) Mpc\({}^{-1}\), i.e. \(0.62h^{-1}\) Mpc. The accuracy reduction of emulated gas fields by using a dark matter simulation with different settings of baryonic effect is expected, especially when the shocks and feedbacks are absent. In Section 5.4, we will discuss this in more detail and the possible treatment/improvements using IllustrisTNG where the difference of the simulation parameters is much more complicated.

### Lyman Alpha forest Absorption (\(F\))

In this section, we discuss the Lyman alpha forest prediction quality using emulated HI density and temperature using Horizon-noAGN and especially what happens when the baryonic feedback is absent using Horizon-DM.

We follow the Equation (23) and calculate the normalised Lyman-\(\alpha\) absorption. Figure 13 shows Lyman-\(\alpha\) absorption derived from the emulator for a sample skewer, and the flux computed from the gas parameters from Horizon-noAGN as a ground truth. The bottom panel shows the residual, defined as

\[\mathrm{Residual}\equiv\exp(-\tau_{\mathrm{True}})-\exp(-\tau_{\mathrm{Pred}}) =F_{\mathrm{True}}-F_{\mathrm{Pred}}. \tag{39}\]

We also include the flux residual estimated by FGPA as a benchmark model6. A visual inspection shows that the emulated spectra agree with the overall residual staying within the \(\pm 0.1\) fiducial error. Some prominent errors exist in the highly absorbed region where _LyAl-Net_ shows limited performance in extreme-density environments. However, it shows the overall prediction performance better than the FGPA model. Nonetheless, Individual-skewer residuals do not provide enough information to assess the emulated field quality due to a few reasons. The absorption requires three different parameters which can affect different behaviours of the emulated flux:

Footnote 6: We include the Voigt profile described in Section 2.1

1. the density \(n_{\mathrm{HI}}\) directly affects the amount of absorption;
2. the temperature \(T_{\mathrm{HI}}\) affects the absorption width due to thermal broadening;
3. the line-of-sight velocity affects the absorption frequency, which will appear as the a shift on spectra

Therefore, to suppress the small-scale fluctuation and detect the systematic error of the emulation, we perform an ensemble average of the transmitted flux along the XY-plane to obtain the mean transmitted flux as a function of redshift (\(\langle F(z)\rangle\)) and illustrate in Figure 14. The top panel of the figure shows that the fluctuations of the _LyAl-Net_ flux versus the true flux (referred to as "ground truth") agree on all the redshift ranges. The only main difference is the amount of absorption from emulated fields. The bottom panel shows that the relative error of _LyAl-Net_ using Horizon-noAGN and Horizon-DM are approximately 2.5% and 3.5%, respectively. This is directly affected by the underprediction affecting the column density in the high matter density region of _LyAl-Net_. The relative error also shows the prominent fluctuations around redshift \(z\approx 2.35\), 2.40, and 2.43. The stochastic components can explain this fluctuation from the various feedbacks that _LyAl-Net_ might ignore. To alleviate these effects would require designing a new loss function that would mitigate the impact of large deviations during the training. The relative error of modified-FGPA is approximately 1%. However, the 1D distribution of Lyman-\(\alpha\) optical depth \(\tau\) shows that _LyAl-Net_ outperforms the benchmark model, particularly in the low optical depth region where the modified-FGPA fails to predict.

For the two-point summary statistics illustrated in Figure 15, the transfer function using Horizon-noAGN (in blue) shows that

Figure 12: Impact of the alteration of dark matter fields from baryonic physics on the quality of the emulation. **Top Panel** shows the ratio of the transfer function of emulated fields using Horizon-noAGN and Horizon-DM showing that the temperature does not suffer from the absence of baryonic feedback. However, the density of neutral hydrogen loses power, meaning the density is smoother. **Middle Panel** and **Bottom Panel** show comparisons of high-fidelity range cross-correlation functions of emulated density and temperature from Horizon-noAGN (illustrated previously in Figure 7 and 10) and Horizon-DM.
the amplitude stays above \(0.90\) up to \(k\approx 10h\) Mpc\({}^{-1}\). We note that \(T_{F}(k)\) the initial drop followed by the slow rise in large scales (\(k\approx 0.1h\) Mpc\({}^{-1}\)) implies the predicted flux is smoother than the ground truth. When turning baryonic feedback off, the flux emulated from Horizon-mD behaves the same way as Horizon-noAGN with a smaller amplitude. This reduction in the amplitude is the direct impact of the \(n_{\rm HI}\) underprediction discussed in Section 4.5. Meanwhile, the modified-FGPA model shows a higher amplitude in larger scales which also reflects on the ensemble average in Figure 14 being closer to Horizon-noAGN. The transfer function, however, noticeably diverges when \(k>3\), which implies that _LyAl-Net_ successfully outperforms the benchmark model on smaller scales. The cross-correlation function on the bottom panel shows a remarkable agreement between the emulated fluxes and the ground truth up to \(k\approx 6h\) Mpc\({}^{-1}\) for using Horizon-noAGN and \(k\approx 9h\) Mpc\({}^{-1}\) for Horizon-mM. the modified-FGPA also shows a similar level of agreement as _LyAl-Net_. However, the spatial distribution is better with Horizon-DM. We can think of two possible reasons for this: the cross-correlation of Horizon-DM is slightly better (see Figure 12), or its velocity field is coincidentally closer to the velocity field of the gas. We will investigate this in future work.

We now consider if the inaccuracies of the model would be problematic to model observations. For a practical example, Bautista et al. (2015) modelled mock quasar spectra as

\[f(\lambda)=\left\{[F(\lambda)\cdot C(\lambda)]*\tilde{W}\left(\lambda,R_{p},R _{w}\right)+N(\lambda)\right\}\cdot M(\lambda)+\delta f_{\rm sky}(\lambda). \tag{40}\]

Where \(f(\lambda)\) is a quasar flux, \(F(\lambda)\) the Lyman-\(\alpha\) absorption, \(C(\lambda)\) is a quasar continuum, \(\tilde{W}\left(\lambda,R_{p},R_{w}\right)\) is a BOSS pixelisation kernel, \(N(\lambda)\) is a noise from Bautista et al. (2015), \(M(\lambda)\) is a correction linear function, and \(\delta f_{\rm sky}\left(\lambda\right)\) is the added sky subtraction residual.

Let us consider a Lyman-\(\alpha\) absorption from full hydrodynamical simulation to be \(F(z)\). the output from _LyAl-Net_ can be expressed as

\[F(z)_{LyAl-Net}=F(z)B(z) \tag{41}\]

where \(B(z)\) is an arbitrary intrinsic bias from _LyAl-Net_ responding to the \(N\)-body simulation input. Therefore, one has to take into account a systematic bias from the prediction model. Therefore, the Equation (40), \(F(\lambda)\) now has an extra term of \(B(\lambda)\), which has to be treated.

If we assume as a minimum requirement that the QSO's spectra signal-to-noise ratio be \({\rm SNR}=10\) and further take the bias to have a linear scaling with \(2.5\%\) bias. Hence, the requirement for SNR will be changed to \(10.3\).

## 5 Transfer Learning with IllustrisTNG

In the previous sections, we test the performance of the direct application of _LyAl-Net_ on simulations. The real universe, however, is not likely to follow any of those. Therefore, it is essential to check the portability of this model. To guarantee this, we need a framework that generalises to other physical feedback and other updates of the gas physics. In this section, we explore the efficiency of the porting of the trained _LyAl-Net_, which is trained on Horizon-noAGN, to different gas physics. We are using IllustrisTNG as a benchmark, and a validation set, to test the accuracy of the predicted density and temperature fields. We first test this by applying _LyAl-Net_ straightaway (out of the box), and then we

Figure 13: **Top Panel:**A comparison of Lyman-\(\alpha\) forest normalised fluxes from the same sample skewer of the ground truth and emulated flux, labelled as Ground Truth and _LyAl-Net_, respectively. The fluxes are calculated from Horizon-noAGN and emulated HI parameters from _LyAl-Net_ using Horizon-noAGN and Horizon-DM dark matter. **Bottom Panel:** Residuals of emulated flux and true fluxes along with fiducial error plotted in dashed line at \(-0.1\), and \(0.1\).Both _LyAl-Net_ and modified-FGPA show that the values mostly stay within the fiducial region while the residuals and peaks within the highly absorbed region with _LyAl-Net_ being more stable.
will discuss the possibility of improving the prediction results by using a technique inspired by transfer learning.

### TNG100 specification

IllustrisTNG is a set of simulations with different settings and resolutions. We choose TNG100 since it has the volume closest to the Horizon-noAGN simulation. TNG100 is a hydrodynamic cosmological simulation where its volume spans a cube with a dimension of \(75h^{-1}\,\mathrm{Mpc}\) on each side. The simulation relies on AREPO, which solves the dynamics equation on a moving unstructured mesh defined by the Voronoi tessellation of a set of discrete points (Springel 2010). We picked snapshot 28, which corresponds to the redshift \(z=2.58\). It is the closest snapshot to the trained _LyAl-Net_ with Horizon-noAGN \(z=2.5\).

### Precomputation

As mentioned above, the data structure of TNG100 is slightly different from Horizon-noAGN, as it is based on the Voronoi tessellation. To derive the dataset of the hydrodynamic quantities for our test, we deploy the same technique using the same adaptive filter to efficiently assign the dark matter and gas fields into a mesh structure with the same resolution as the Horizon-noAGN training set (\(\sim 9.78h^{-1}\,\mathrm{Mpc}\)) which yields \(768\times 768\times 768\) voxels.

We derived the temperature of the gas from the internal energy, \(u\), of the TNG100 simulation by using the relation 7

Footnote 7: We followed the guide from the official The IllustrisTNG Project ([https://www.tng-project.org/data/docs/faq/gen5](https://www.tng-project.org/data/docs/faq/gen5)).

\[T=\mu(\gamma-1)\frac{u}{k_{B}}\, \tag{42}\]

where \(T\) is the temperature of the gas in each cell, \(\gamma=5/3\) is the adiabatic index for an ideal monoatomic gas, \(k_{B}\) is the Boltzmann constant, and \(\mu\) is the mean molecular weight. The mean molecular weight is derived from the electron abundance, a ratio of the number density of free electrons and a total hydrogen number density (Kim et al. 2022): This relation is expressed by,

\[\mu=\frac{4m_{p}}{1+3X_{H}+4X_{H}x_{e}}\, \tag{43}\]

where \(m_{p}\) is proton mass and \(x_{e}=n_{e}/n_{H}\) is the electron abundance which is a fraction of free electrons with respect to the number density of total hydrogen, \(X_{H}\) is the hydrogen mass fraction which we assume to be constant and uniform.

The database of the IllustrisTNG simulations does not publicly provide the neutral hydrogen density. We instead choose

Figure 14: One-point statistics comparisons of Lyman-\(\alpha\) absorption between Horizon-noAGN (Ground Truth) and emulated absorptions, using _LyAl-Net_ and the benchmark model (mod-FGPA) **Top Panel**: A mean transmitted flux as a function of redshift, \((F(z))\), shows that _LyAl-Net_ tracks well with the ground truth. **Middle Panel**: A comparison of relative error of mean transmitted flux as a function of redshift between flux emulated. Both _LyAl-Net_’s relative errors exhibit the same shape, while Horizon-noAGN has the mean closer to the ground truth. The benchmark has a closer mean transmitted flux. **Bottom Panel**: A comparison of optical depths 1D distribution shows _LyAl-Net_ flux using Horizon-noAGN and Horizon-DM perform better than the benchmark model.

Figure 15: The comparisons of transfer functions of and cross-correlation function (Top panel) and cross-correlation function (Bottom panel) of the emulated using _LyAl-Net_ on Horizon-noAGN, Horizon-DM, and modified-FGPA. The transfer function of modified-FGPA noticeably diverges when \(k>3\) which implies that _LyAl-Net_ performs better and more stable in the smaller scales.
to estimate it from the total hydrogen density, which is available. Therefore, we estimated the density of neutral hydrogen, \(n_{\rm HIMmock}\), by scaling it from its mean relative abundance in the Horizon-noAGN simulation. This can be expressed as

\[n_{\rm HI}{}_{\rm mock}^{\rm TNG100}=\left(\frac{n_{\rm HI}{}^{\rm HnoAGN}}{n_{ \rm H}^{\rm HnoAGN}}\right)n_{\rm H}^{\rm TNG100}. \tag{44}\]

The above equation is an important simplification of the complicated physics of the hydrogen ionisation state. In future work, we intend to better model this effect for the training set through adequate computation of this equilibrium and matching it to Horizon-noAGN. For a performance comparison, we train the _LyAl-Net_ to predict the total hydrogen density, \(n_{H}^{\rm HnoAGN}\), with the same Horizon-noAGN, allowing us to scale the density using Equation (44) to obtain the TNG100 mock neutral hydrogen density.

### Assessment of LyAl-Net capabilities

We now discuss the performance of the emulated neutral hydrogen density and temperature using _LyAl-Net_ on dark matter overdensity on TNG100.

#### 5.3.1 Hydrogen Density

We perform the same KDE analysis for one-point statistics on the 400 slices (\(768\times 768\times 400\) voxels). Figure 16 shows that _LyAl-Net_'s prediction bias remains relatively stable and stays around the 10% fiducial bias in the first and ninetieth percentile, as highlighted in blue. The overall prediction is more prominent compared to Horizon-noAGN counterpart, and this is expected to occur. Figure 17 illustrates a transfer function (top panel) and correlation rate (bottom panel), which stay above 90% up to \(k\approx 7h\) Mpc\({}^{-1}\) and \(6h\) Mpc\({}^{-1}\) respectively. In addition, the accuracy drops from \(k\approx 10\)\(h\) Mpc\({}^{-1}\) compared to Horizon-noAGN.

#### 5.3.2 Temperature

Figure 18 compares the 1D distribution between true and predicted temperature. The predicted temperature appears to have a higher median with a shorter distribution tail in the high temperature where \(T>10^{4}\) K compared to TNG100. This overall underprediction also appears on the bias in Figure 19. This should be mainly the result of the lack of the AGN feedback. Two-point statistics in Figure 20 show that the transfer function of _LyAl-Net_ stays under 50%. This lack of power shows that the fluctuation of the predicted temperature on _LyAl-Net_ is lower than the ground truth, which is also implied from the bias plot. The cross-correlation function on the bottom panel, on the other hand, shows that _LyAl-Net_ performs the best up to \(k\approx 1h\) Mpc\({}^{-1}\) while the _LyAl-Net_ with Horizon-noAGN performs best up to \(k\approx 3h\) Mpc\({}^{-1}\).

Figure 16: The comparisons of the prediction bias before and after applying the transfer learning on a TNG100 snapshot 28 (Sec 5.4).The light blue band represents the density within the 1st and 90th percentile. It shows the general improvement in the high-density regime after the 0th percentile while maintaining a similar bias across the density range compared to the out-of-the-box _LyAl-Net_.

Figure 17: Top panel: transfer functions and Bottom panel: cross-correlation function of the neutral hydrogen temperature emulated using _LyAl-Net_. The blue line represents the performance of _LyAl-Net_ out-of-the-box, and the orange line represents the performance of _LyAl-Net_ with transfer layers. We also used the (high fidelity) 99.5th percentile range masking corresponds to \(1.13\times 10^{-7}\) to \(1.26\times 10^{-4}\) cm\({}^{-3}\)

Figure 18: The decimal logarithm Temperature 1D distribution comparison Between _LyAl-Net_ prediction and the true value of TNG100
### Transfer learning

_LyAl-Net_ mainly translates the dark matter overdensity to HI density and temperature. As the previous sections show, it relies on universal spatial correlations between the matter density and gas parameters. Even though dark matter generally traces gas clouds, there is a spatial limit that _LyAl-Net_ can predict. Results presented in the previous sections show that applying _LyAl-Net_ onto other dark matter simulations, in this case, IllustrisTNG may result in inaccurate emulated gas fields. This is expected since, at minimum, baryonic feedbacks are not identical to Horizon-noAGN. The main goal of applying _LyAl-Net_ is to have an emulator framework that can generalise and mimic the gas physics for a different equation of states, such as cooling rate and UV index, with the possibility of generalising to different cosmology. Achieving this portability of _LyAl-Net_ would ultimately unlock its application for cosmological analysis and the computation of other gas quantities in the IGM physics. In particular, cosmological inference requires lots of forward simulation to properly marginalise systematic effects and uncertainty on gas physics. A generic and cheap way of generating maps of the IGM would be crucial for future surveys, such as Lyman-\(\alpha\)-forest ones. In this section, we explore the possibility of slightly augmenting the _LyAl-Net_ to obtain a hydrodynamic simulation which contains a desirable cosmology. _LyAl-Net_ can be considered as a complex non-linear transformation, where the output gas field (\(y\)) is the response of the input dark matter (\(x\)), which can be expressed as

\[y=f(x)\,. \tag{45}\]

The emulator, _LyAl-Net_, is only able to achieve adequate transformations of dark matter to gas properties on a restricted range of scales and for a given physics. A _LyAl-Net_ trained on Horizon-noAGN has good emulation capabilities of the gas equation of state. The first thing we should try to preserve is exactly that physical property, equivalent an nth-order expansion of an arbitrary analytical solution of the gas physics. By adjusting the input and output _LyAl-Net_, without retraining the entire network, we may achieve higher-order accuracy to the emulated gas parameters. Based on the Equation (45), the equation yields

\[\ddot{y}=g[f[h(x)]]\,, \tag{46}\]

where \(h(x)\) and \(g(y)\) are arbitrary transformation functions of the input dark matter overdensity and the output gas fields, respectively. To test this idea, we introduce a custom transformation function which takes the form

\[L(x)=Ax+B, \tag{47}\]

where \(A\) and \(B\) are two free parameters. This function performs a linear scaling and shifting in the log space for _LyAl-Net_. We assume that the physical feedbacks can be summarised and encoded into a few parameters, reducing the bias, especially the diffuse region of the gas, and ultimately improving the Lyman-\(\alpha\) flux. As indicated in the introduction, to first order, this is expected to improve the emulation of the gas equation of state and density. We note that a better transformation may exist, which will be investigated in future works. We implemented this function into _LyAl-Net_ as a layer without an activation function, which we will refer to as a _transfer layer_. The transfer layers are placed before the dark matter input and after the prediction from _LyAl-Net_, effectively encasing it. By implementing the transfer layers in this way, we lose the linearity of the function

\[y^{\prime}=L(y;A,B)=L(LyA|I(x);C,D];A,B)\,. \tag{48}\]

In Figure 21 we provide schematics on implementing transfer layers. The weights of _LyAl-Net_ are frozen. As for the original training of the network, the mean square error was used to optimise these parameters. This allows us to leverage a GPU and TensorFlow to optimise the coefficients of transfer layers, namely A, B, C, and D, based on the \(\chi^{2}\) loss function of the ground truth and predicted values of IllustrisTNG-100. Only 512 sub-boxes were sampled, \(27^{3}\) voxels each, equivalent to \(\sim 1\%\) of the total volume. Table 4 summarises the parameters obtained through the described optimisation. We then re-computed the emulated total hydrogen density \(n_{\rm HI}\) and temperature \(T_{\rm HI}\).

Figure 19: The comparisons of the prediction bias before and after applying the transfer learning on the snapshot 28 of TNG100 redshift \(z=2.58\). We use the same colour and style conventions as in Figure 16. It shows the general improvement in the high-temperature regime. We note that we nearly reach the requirement of the diffuse IGM limit with this modification.

Figure 20: The comparisons of transfer functions of and cross-correlation function (Top panel) and cross-correlation function (Bottom panel) of the neutral hydrogen temperature emulated using _LyAl-Net_, and we also used the nbm fidelity range masking. The blue solid line shows the results obtained with the native _LyAl-Net_, whereas the orange solid line is the result of the transfer after optimising for the four additional parameters. We note a dramatic improvement of the performances of the network to different physics.
#### 5.4.1 Hydrogen density

For the hydrogen gas density, the bias function in orange in Figure 16 shows that the overall shift in the blue region (1st-90th percentile) is slightly outside the fiducial bias region on the one hand. On the other hand, the bias of extremely high density is suppressed. For the two-point statistics, Figure 17 shows an overall improvement in the transfer function amplitude, though at very small scales at \(k=10h\) Mpc\({}^{-1}\) the predicted density is overboosted by 10% in the process. The cross-correlation, on the other hand, shows an impressive improvement across visible spatial scales. The wave number corresponding to a correlation rate of 90% has increased by \(\sim 40\%\), from \(6.5h\) Mpc\({}^{-1}\) to \(9h\) Mpc\({}^{-1}\).

#### 5.4.2 Temperature

After applying transfer layers, the bias function of temperature in Figure 19 shows an improvement for the temperature range within the 1st and 90th percentile, where the emulator predicts mean values stable within the 10% fiducial error. From Figure 20, the transfer function shows a drastic improvement of the amplitude being nearly double compared to the original _LyAl-Net_. The cross-correlation also shows that the spatial distribution accuracy now increases to \(k\sim 3h\) Mpc\({}^{-1}\), which performs similarly to the Horizon-noAGN counterpart.

### Lyman-\(\alpha\) forest absorption rate

We illustrated a sample skewer of a normalised Lyman-\(\alpha\) flux of TNG100 in Figure 22. The top panel directly compares ground truth in blue, a calculated flux from out-of-the-box _LyAl-Net_ in orange, and _LyAl-Net_ equipped with transfer layers in green. Fluxes from both configurations of _LyAl-Net_ trace the true flux well. The residual plot on the bottom panel shows a similar performance to this skewer. However, the original _LyAl-Net_ performs poorly in the highly-absorbed region. After applying the transfer learning layer, this region has been improved significantly. This is mainly the result of temperature field correction via applying transfer layers, which tends to smooth out and broaden the absorption. The mean transmitted flux comparison in Figure 23 shows _LyAl-Net_ has a relative error of \(\approx 4\%\). This error drops drastically with the transfer layers. We can see that transfer layers improve the absorption accuracy, illustrated on the transfer function from the top panel Figure 24. The amplitude increases to 0.9 up to \(k\approx 10h\) Mpc\({}^{-1}\). The spatial accuracy measured by the cross-correlation function at the bottom shows only a tiny improvement. This implies that the current transfer layers improve mainly the amount of absorption. Moreover, the overall performance is nearly identical to Horizon-noAGN flux prediction, as shown in Figure 15.

This finding shows that _LyAl-Net_ with the experimental transfer layer model shows the framework's portability. It can also be extended to extrapolate the different configurations of the baryonic feedback for future cosmological emulation. Moreover, the eight parameters (four for gas density and four for temperature) may be inferred jointly with the density field in a Bayesian hierarchical model of Lyman-\(\alpha\) forest observations, for example, in Porqueres et al. (2020).

## 6 Discussion and conclusion

We have presented _LyAl-Net_, a neural network architecture and training procedure to derive hydrogen gas field quantities from dark matter density and velocity fields at a cosmological scale. We have also shown a number of benchmarks on its resilience, portability, and capabilities for producing mock Lyman-\(\alpha\) forest observables. One of the most relevant points of this network is its portability to other gas models. In the existing literature, the mean transmitted flux from simulations has been accounted

\begin{table}
\begin{tabular}{l c c c c} \hline \hline \multicolumn{1}{c}{**Gas**} & \multicolumn{3}{c}{**Transfer Layer Parameters**} \\ \cline{2-5}
**Parameter** & **A** & **B** & **C** & **D** \\ \hline Density & 1.180 & 0.657 & 1.133 & 0.084 \\ Temperature & 0.898 & -0.862 & 1.34 & -0.466 \\ \hline \hline \end{tabular}
\end{table}
Table 4: The parameters of the transfer layers for IllustrisTNG-100 total hydrogen density and temperature at the redshift \(z=2.58\). All parameters are optimised by the MSE function using a sample size of 512 sub-boxes, 27\({}^{3}\) voxels each.

Figure 21: A simple schematic of the modified emulator to implement the transfer learning procedure. The weights of the core network, _LyAl-Net_, are now frozen, meaning they are not further trained. We introduce transfer layers between the dark matter overdensity and this network to tune the input to the new problem. We only optimise \(A\), \(B\), \(C\), and \(D\) to fine-tune the prediction results, which can also be considered high-order corrections.
for using a simplistic re-scaling factor to match the observations (Regan et al., 2007; Lukic et al., 2014; Borde et al., 2014). This method is equivalent to re-scaling the intensity of the UV background and the photoionisation rate from observed \(\tau\)(Meiksin, 2009). However, this re-scaling does not provide a deeper insight into what happens at the IGM level. In contrast, recent work by Borde et al. (2014) has shown that the UV background is dominated by the UV background and the photoionisation rate.

Figure 24: A performance assessment using two-point correlation function of transmitted fluxes emulated by _LyAl-Net_ (labelled in blue) and _LyAl-Net_ with transfer layers (labelled in orange) from IllustrisTNG-100 at \(z=2.58\) dark matter overdensity, where the top panel is a transfer function comparison, and the bottom panel is a cross-correlation function comparison.

Figure 23: **Top Panel:** A comparison of the IllustrisTNG-100 mean transmitted flux as a function of redshift. A true flux referred to as _Ground Truth_ in solid blue, and two emulated fluxes calculated from emulated IGM from _LyAl-Net_; 1.) the trained _LyAl-Net_ in orange and 2.) trained _LyAl-Net_ with transfer layers in green. **Bottom Panel:** A comparison of relative error plots of emulated mean transmitted fluxes using trained _LyAl-Net_ and _LyAl-Net_ with transfer layers. It shows the relative error drastically reduces after applying transfer layers.

Figure 22: **Top Panel:**A comparison of Lyman-\(\alpha\) forest normalised fluxes from the same sample skewer of the ground truth and emulated fluxes, labelled as Ground Truth,_LyAl-Net_, and _LyAl-Net_ with Transfer Layers respectively from IllustrisTNG100 snapshot 28, redshift \(z=2.58\). The _LyAl-Net_ with transfer layers shows an improvement, especially within the highly absorbed region. **Bottom Panel:** Residuals of the fluxes along with 10% fiducial error plotted in dashed line shows that the residual values of both configurations mostly stay within the fiducial region
works focused on improving the accuracy of the field, which yields a better emulation of the flux absorption. These works also utilise a neural network approach to generate hydrodynamical simulations from \(N\)-body simulations of dark matter. For example, Horowitz et al. (2022) used a fully convolutional variational auto-encoder similar to U-Net architecture. A follow-up paper, Wadekar et al. (2020) used a special _U-Net_ architecture to predict the neutral hydrogen density from dark matter simulation. However, to capture the physics of underdense and overdense regions, they had to combine two separate neural networks, whose output is selected depending on the amplitude of the local dark matter density. Sinigaglia et al. (2022) used a hierarchical domain-specific machine learning to obtain Lyman-\(\alpha\) forest flux from dark matter overdensity (nicknamed BAM), but also requires a reconstruction of HI and HII.

In this work, we improved on several aspects of existing scientific literature: a characterisation of the stochastic properties of the emulated fields with respect to the simulated fields for different implementations of the hydrodynamic physics; a flexible layer to transpose our results to unknown new physics; a simplification of the overall U-net architecture for Lyman-\(\alpha\); and a mitigation of particle shot-noise effect in voids using an adaptive filter for the dark matter as used in LyMAS (Peirani et al., 2014; Peirani et al., 2022). The first point is critical for using our model as part of a Bayesian analysis. We need to know what are the limits of the model, at least at the level of the 1-point statistics, to be able to fit it to observations, e.g. as a part of a field-level analysis (e.g. BORG in Porqueres et al.2020), a simulation-based approach (e.g. SELFI, Leclercq et al.2019), or a standard correlation function analysis (Slosar et al., 2011).

The second point in the above list is the most important, as we have no warranty that the parameters (e.g. star formation, AGN activity) governing the simulated fields (gas density, temperature) agree with observations. A portable model allows us to probe this directly and note the model's limits from observations. A similar procedure was followed in Lee et al. (2015); Palanque-Delabrouille et al. (2015), but limited to matching 2-point statistics for the FGPA model. Once _LyAl-Net_ model is calibrated on observations, it also indirectly allows us to utilise these emulated fields as a by-product. The emulator predicts the hydrodynamics simulation at the field level, which allows _LyAl-Net_ to predict the IGM properties from any \(N\)-body simulation code. We made several improvements with respect to LyMAS-2 (Peirani et al., 2014). _LyAl-Net_ does not require a smoothing filter, but the emulated fields' resolution is fixed, depending on the training set. Therefore, the framework can be applied to any desired resolution. To check the portability of _LyAl-Net_, we first assessed by applying a pre-trained network naively onto IllustrisTNG-100 dark matter. The physics of the baryonic feedbacks are present, whereas it was absent in the training set, which originates from Horizon-noAGN. The results, given in Figure 5.5, showed that the relative error of \(\langle F(z)\rangle\) increases, which was directly impacted by the underprediction of the hydrogen density as expected. We solved this by a simple numerical treatment using transfer layers, drastically improving emulated flux and matching the performance emulated Horizon-noAGN. The transfer layers we introduced may be simply interpreted as a tuning of the equation-of-state by shifting the median of \(n_{H}\) and \(T\) closer to the true values of TNG100. Of course, the form of the equation we used has only a few parameters, only allowing one-dimensional remapping of the value, which are not expected to be able to emulate the higher-order corrections. Simply put, we can improve the mean transmission without arbitrarily scaling the flux, and we apply transfer layers that fine-tune the equation of state of IGM. This approach provides extra information and estimations of IGM for a given baryonic physics while obtaining the individual fields of hydrodynamical simulations as by-products.

We remind the reader about the mean transmitted flux as a function of redshift in Figure 14. The relative error when using Horizon-noAGN and Horizon-DM behaved in the same way, while Horizon-DM has a higher relative error. From this empirical observation, we can safely claim that _LyAl-Net_ can transform the unit of dark matter into IGM properties based on overdensity fluctuation. However, it has some limitations in emulating the small-scale effects due to the stochastic nature of the baryonic feedbacks in which _LyAl-Net_ might need to be correctly captured or provide only a mean response, possibly at a low order of accuracy. The approach followed in this work still has some limitations. Notably, we have tested the emulation at a fixed redshift, in term of conformal time wise, and we have no warranty of the portability of _LyAl-Net_ to other simulation time. We have not fully explored the possibility of reducing the complexity of the U-net, and a smaller version may allow reaching the same level of accuracy. We seek to improve also the interpretability in future work. For example, while staying compact, the transfer layer may take a more optimal form than Equation 38. We also seek to reduce the free parameters required to emulate a different environment and feedback sufficiently enough to match the scale and resolution for cosmological surveys.

We developed a fast simulator for Lyman-\(\alpha\) forest simulation using a neural network to emulate hydrodynamical simulations of the intergalactic medium. We use Horizon-noAGN as a training set. The current model works well with simulation at \(z\sim 2.5\) while treating the entire simulation independently from redshift. We tested the sensitivity of the emulator using dark matter overdensity from IllustrisTNG-100 at \(z\sim 2.4\). This simulation contains different baryonic feedback from the training set, and the prediction bias occurred as expected. To post-process the emulated fields, we introduced a method called _linear layer_ to calibrate an IGM equation of state. We found that it improves the Lyman-\(\alpha\) absorption flux closer to the full hydrodynamical simulation. This improvement has shown that one can perform a correction from the trained _LyAl-Net_ by a fixed cosmology and baryonic feedbacks, in the form of the linear layer, and the corresponding parameters still need to contain interpretability.

This version of _LyAl-Net_ assumes redshift independence for a given training set, which means it is limited to a single redshift that covers a specific volume. While the transfer layers can fine-tune the IGM equation of state, even with a cosmological simulation with a different redshift, as shown for the IllustrisTNG test. However, the IGM density-temperature relation the IGM density-temperature relation becomes less reliable in higher redshift because the dark matter distribution is less evolved. Moreover, with the observations of higher redshift Lyman-\(\alpha\) forest spectra becoming increasingly available, future work will look into redshift dependency for _LyAl-Net_ so that it can be more resilient and robust.

Another development for _LyAl-Net_ is to integrate into inference models such as DELFI (Alsing et al., 2019), BOLFI (Leclercq, 2018), and SELFI (Leclercq et al., 2019). Such models also require the emulation speed, which is the main priority in this work since such models require a massive amount of data for a usable analysis. _LyAl-Net_ also allows an improvement on classical cross-correlation analysis of large cosmological surveys (Seljak et al., 2006; Amendola et al., 2013; Ivezic et al., 2019). Further, _LyAl-Net_ opens the way for field-level large-scale structure inferences through a combination of Lyman-\(\alpha\) tracers and galaxy clustering at high redshifts with next-generation galaxy
surveys (Jasche & Wandelt 2013; Porquetres et al. 2019; Ivezic et al. 2019; Tsaprazi et al. 2023). Contrary to galaxy clustering, the Lyman-\(\alpha\) forest is sensitive to underdensities and probes the large-scale structure at a higher resolution. Therefore, combining the two is expected to improve the constraining power of field-level inferences of structure formation while providing the ability to emulate the physics of the intergalactic medium without the large computational cost of hydrodynamical simulations.

###### Acknowledgements.

 We thank valuable discussions with Francisco Villaescusus-Navarro, Shy Genel, Simon Prunet, Benjamin D. Wandelt. This work was supported by the ANR BIG4 project, grant ANR-16-CE23-0002 of the French Agence Nationale de la Recherche. This work was supported by the Simons Collaboration on "Learning the Universe". CB acknowledges the financial support from the Sorbonne Center for Artificial Intelligence (SCA). This work was enabled by the research project grant "Understanding the Dynamic Universe" funded by the Knut and Alice Wallenberg Foundation under Drk AAW 2018.0067. This work has made use of the Horizon cluster hosted by the Institut \(A\) Astrophysics de Paris in which the cosmological simulations were post-processed. We thank Stephane Rouberl for running smoothly his cluster for us. This work is conducted within the Aquila Consortium8. We acknowledge the use of the following packages in this analysis: Numpy (Harris et al. 2020), TensorFlow (Aaddi et al. 2015), JAX (Bradbury et al. 2018), IPython (Preze & Granger 2007), Matplotlib (Hunter 2007), Numba (Lam et al. 2015), SciPy (Jones et al. 2001).

Footnote 8: [https://www.aquila-consortium.org/](https://www.aquila-consortium.org/)
## Appendix A Voigt Profile calculation

The cross-section of Lyman-\(\alpha\) is modelled by a Voigt profile. This function is very expensive to compute. We rely on the function scipy.special.wofz to efficiently obtain the results. This software implements the Faddeeva algorithm to compute it (Johnson, 2012), which is a complex scaled complementary error function (Faddeeva and Terent'ev, 1961). We note that this is the most computationally expensive part of the entire pipeline.

## Appendix B Kernel Density Estimation implementation

We estimated probability density functions using scipy.stats.gaussian_kde, which applies a Gaussian kernel with Scott's Rule for bandwidth selection (Scott, 2015). However, the computational time does not increases linearly with the sample size. Therefore, the practical approach is to use a limited number of voxels. For our test, a full simulation volume of Horizon-noAGN has \(1024\times 1024\times 1024\) voxels. We use \(1024\times 1024\times 500\) voxels which is approximately 50% of the total volume for the analysis. For this reduced number of points, the required wall time to generate the density estimate costs a total of approximately 38 hours on AMD EPYC 7702 (128 cores).

## Appendix C Power spectra masking process

The \(n_{\mathrm{HI}}\) prediction suffers much more in the high-density regime where the emulator incorrectly produces invalid predictions. We presume that the bias from the high density is the culprit for the misbehaviour of the diagnostics using two-point statistics documented in the main text. To check this hypothesis, we have selected the density pixels below three different fidelity ranges (see Table 3) and masked the density outside this range to zero. Figure 14 shows the comparisons of density power spectra in the high-fidelity range, which now agrees well with the Horizon-noAGN and confirms the idea mentioned earlier.
