© The Authors 2025

1. Introduction

The Lyman-α forest refers to absorption lines in the spectra of high-redshift quasars resulting from Lyman-α absorption by neutral hydrogen in the intergalactic medium (IGM; for a review, see McQuinn 2016). Even though quasars can be observed at very high redshifts with relatively short exposure times, the scarcity of these sources limits their direct use for precision cosmology. Conversely, Lyman-α forest measurements from a single quasar spectrum provide information about density fluctuations over hundreds of megaparsecs along the line of sight, making this observable an excellent tracer of large-scale structure at high redshifts.

Cosmological analyses of the Lyman-α forest rely on either three-dimensional correlations of the Lyman-α transmission field (ξ_3D; e.g., Slosar et al. 2011) or correlations along the line-of-sight of each quasar (i.e., the one-dimensional flux power spectrum P_1D, Croft et al. 1998; McDonald et al. 2000). The first analyses set constraints on the expansion history of the Universe by measuring baryonic acoustic oscillations (BAO; e.g., Busca et al. 2013; Slosar et al. 2013; du Mas des Bourboux et al. 2020), for which linear or perturbation theory is accurate enough. On the other hand, P_1D analyses measure the small-scale amplitude and slope of the linear power spectrum (e.g., Croft et al. 1998; McDonald et al. 2000, 2005; Zaldarriaga et al. 2001; Viel et al. 2004), the nature of dark matter (e.g., Seljak et al. 2006a; Viel et al. 2013; Iršič et al. 2017; Palanque-Delabrouille et al. 2020; Rogers & Peiris 2021a; Iršič et al. 2024), the thermal history of the IGM (e.g., Viel & Haehnelt 2006; Bolton et al. 2008; Lee et al. 2015; Walther et al. 2019; Boera et al. 2019; Gaikwad et al. 2020, 2021) and the reionization history of the Universe (see the reviews Meiksin 2009; McQuinn 2016). In combination with cosmic microwave background constraints, P_1D analyses also set tight constraints on the sum of neutrino masses and the running of the spectral index (e.g., Spergel et al. 2003; Verde et al. 2003; Viel et al. 2004; Seljak et al. 2005, 2006b; Palanque-Delabrouille et al. 2015, 2020).

Unlike ξ_3D studies, P_1D analyses go deep into the nonlinear regime and require time-demanding hydrodynamical simulations (e.g., Cen et al. 1994; Miralda-Escudé et al. 1996; Meiksin et al. 2001; Lukić et al. 2015; Bolton et al. 2017; Walther et al. 2021; Chabanier et al. 2023; Puchwein et al. 2023; Bird et al. 2023). Naive analyses would demand running millions of hydrodynamical simulations, which is currently unfeasible. Rather, the preferred solution is constructing fast surrogate models that make precise predictions across the input parameter space using simulation measurements as the training set. The main advantage of these surrogate models, known as emulators, is reducing the number of simulations required for Bayesian inference from millions to dozens or hundreds. In the context of Lyman-α forest studies, the first P_1D emulators involved simple linear interpolation (McDonald et al. 2006) and progressively moved toward using Gaussian processes (GPs; Sacks et al. 1989; MacKay et al. 1998) and neural networks (NNs; McCulloch & Pitts 1943); for instance, Bird et al. (2019, 2023), Rogers et al. (2019), Walther et al. (2019), Pedersen et al. (2021), Takhtaganov et al. (2021), Rogers & Peiris (2021b), Fernandez et al. (2022), Molaro et al. (2023), and Cabayol-Garcia et al. (2023).

The primary purpose of this work is to provide consistent predictions for Lyman-α forest clustering statistics from linear to nonlinear scales. There are three main approaches to achieve this. The first relies on perturbation theory (e.g., Givans & Hirata 2020; Chen et al. 2021; Ivanov 2024), which delivers precise predictions on perturbative scales at the cost of marginalizing over a large number of free parameters. The second involves emulating power spectrum modes measured from a suite of cosmological hydrodynamical simulations, which provides precise predictions from quasilinear to nonlinear scales. The main limitation of this approach is that accessing the largest scales used in BAO analyses, r ≃ 300 Mpc, would require hydrodynamical simulations at least three times larger than this scale (e.g., Angulo et al. 2008), which is currently unfeasible due to the computational demands of these simulations. Nonetheless, there are recent strides in this direction using approximated methods (e.g., Jacobus et al. 2023).

Instead, we adopt a third approach: emulating the best-fitting parameters of a physically motivated Lyman-α clustering model to measurements from a suite of cosmological hydrodynamical simulations as a function of cosmology and IGM physics. In what follows, we refer to this strategy as FORESTFLOW¹. In this work, we use conditional normalizing flows (cNF; Winkler et al. 2019; Papamakarios et al. 2019) to emulate the two Lyman-α linear biases (b_δ and b_η), which completely determine the large-scale behavior of P_3D in conjunction with the linear power spectrum, along with six parameters capturing small-scale deviations of P_3D from linear theory as a function of six parameters capturing the cosmological and IGM dependence of Lyman-α clustering (Pedersen et al. 2021). We show that this strategy enables precise P_3D predictions from nonlinear scales to arbitrarily large (linear) scales even when training on a suite of simulations with moderate size (Pedersen et al. 2021; Cabayol-Garcia et al. 2023). FORESTFLOW also allows for the prediction of any statistic derived from P_3D without requiring interpolation or extrapolation. For instance, we can compute ξ_3D by taking the Fourier transform of P_3D or derive P_1D by integrating its perpendicular modes

$$\begin{array}{c}\hfill P_{1\mathrm{D}}(k_{\Vert })={(2\pi )}^{-1}{\displaystyle {\int }_0^{\infty }}\mathrm{d}{k}_{\perp }{k}_{\perp }{P}_{3\mathrm{D}}(k_{\Vert },k_{\perp }),\end{array}$$(1)

where k_∥ and k_⊥ indicate parallel and perpendicular modes, respectively.

The release of FORESTFLOW is quite timely for BAO and P_1D analyses of the ongoing Dark Energy Spectroscopic Instrument survey (DESI; DESI Collaboration 2016), which will quadruple the number of line-of-sights observed by first the Baryon Oscillation Spectroscopic Survey (BOSS; Dawson et al. 2013) and its extension (eBOSS; Dawson et al. 2016). DESI has already proven the constraining power of Lyman-α studies by measuring the isotropic BAO scale with ≃1% precision from the Data Release 1 (DESI Collaboration 2024) and P_1D at nine redshift bins with a precision of a few percent from the Early Data Release (Ravoux et al. 2023; Karaçaylı et al. 2024). In addition to being used for BAO and P_1D studies, FORESTFLOW has the potential to extend toward nonlinear scales the current full-shape analyses of ξ_3D (Cuceu et al. 2023; Gerardi et al. 2023) and P_3D (Font-Ribera et al. 2018; de Belsunce et al. 2024; Horowitz et al. 2025), and can be used to interpret alternative three-dimensional statistics (Hui et al. 1999; Font-Ribera et al. 2018; Abdul Karim et al. 2024).

The structure of this paper is as follows: we describe the strategy adopted by FORESTFLOW, the input data for the cNF, and its architecture in Sects. 2, 3, and 4, respectively. Next, we assess the performance of FORESTFLOW in Sect. 5 and highlight some novel analyses enabled by this framework in Sect. 6. Finally, we summarize the main findings and conclude in Sect. 7. Throughout this paper, all statistics and distances are in comoving units.

2. Strategy adopted by FORESTFLOW

In FORESTFLOW, we emulate the parameters of a physically motivated parametric model for Lyman-α forest clustering as a function of parameters that capture the influence of cosmology and IGM physics on this observable. This section begins with an overview of the physically motivated model, followed by an introduction to the parameters used to characterize the dependence of Lyman-α forest clustering on cosmology and IGM physics.

2.1. Parametric model for P_3D

We can express fluctuations in the Lyman-α forest flux as ${\delta }_F(\mathbf{s})={\overline{F}}^{-1}(\mathbf{s})F(\mathbf{s})-1$, where F = exp(−τ) and $\overline{F}$ are the transmitted flux fraction and its mean, respectively, τ is the optical depth to Lyman-α absorption, and s is the redshift-space coordinate. On linear scales, these fluctuations depend upon the matter field as follows (e.g., McDonald 2003)

$$\begin{array}{c}\hfill {\delta }_F=b_{\delta }\delta +b_{\eta }\eta ,\end{array}$$(2)

where δ refers to matter density fluctuations, η = −(a H)⁻¹(∂v_r/∂r) stands for the dimensionless line-of-sight gradient of radial peculiar velocities, a is the cosmological expansion factor, H is the Hubble expansion factor, v_r is the radial velocity, and r stands for the radial comoving coordinate. The linear bias coefficients b_δ and b_η capture the response of δ_F to large-scale fluctuations in the δ and η fields, respectively.

Following McDonald (2003), we decompose the three-dimensional power spectrum of δ_F into three terms

$$\begin{array}{c}\hfill P_{3\mathrm{D}}(k,\mu )={(b_{\delta }+b_{\eta }f{\mu }^2)}^2{D}_{\mathrm{NL}}(k,\mu )P_{\mathrm{lin}}(k),\end{array}$$(3)

where f = dlog G/dlog a is the logarithmic derivative of the growth factor G, μ = k⁻¹k_∥ is the cosine of the angle between the Fourier mode and the line of sight, (b_δ + b_η f μ²)² accounts for both linear biasing and large-scale redshift space distortions (Kaiser 1987; McDonald et al. 2000), P_lin is the linear matter power spectrum², and D_NL is a physically motivated parametric correction accounting for the nonlinear growth of the density field, nonlinear peculiar velocities, thermal broadening, and pressure.

The large-scale behavior of P_3D is set by the bias coefficients b_δ and b_η together with the linear power spectrum, and the latter can be computed using a Boltzmann solver (e.g., Lewis et al. 2000; Lesgourgues 2011). Therefore, the emulation of the two Lyman-α linear biases enables predicting P_3D on arbitrarily large (linear) scales³. In contrast with direct emulation of power spectrum modes, this approach only requires simulations large enough for measuring the two Lyman-α linear biases precisely.

Predicting P_3D on small scales is more challenging than on large scales due to the variety of effects influencing this statistic in the nonlinear regime (e.g., McDonald 2003). In this work, we capture these small-scale effects using the physically motivated Arinyo-i-Prats et al. (2015) parameterization

$$\begin{array}{c}\hfill D_{\mathrm{NL}}=exp\{(q_1{\mathrm{\Delta }}^2+q_2{\mathrm{\Delta }}^4)[1-{\left(\frac{k}{k_{\mathrm{v}}}\right)}^{a_{\mathrm{v}}}{\mu }^{b_{\mathrm{v}}}]-{\left(\frac{k}{k_{\mathrm{p}}}\right)}^2\},\end{array}$$(4)

where Δ²(k)≡(2π²)⁻¹k³P_lin(k) is the dimensionless linear matter power spectrum and the free parameters k_v and k_p are in  Mpc⁻¹ units throughout this work. The terms involving {q₁,  q₂}, {k_v,  a_v,  b_v}, and {k_p} account for nonlinear growth, peculiar velocities and thermal broadening, and gas pressure, respectively. The expression above does not account for a shot noise term (e.g., Iršič & McQuinn 2018). While Givans et al. (2022) successfully described P_3D and P_1D measurements down to highly nonlinear scales using this formulation, it is possible that the shot noise contribution was implicitly absorbed into their fit through free parameters representing other effects. A more detailed investigation of shot noise is deferred to future work.

In the top panel of Fig. 1, dotted lines show the ratio of measurements from the CENTRAL simulation at z = 3 and the linear power spectrum, while solid lines do so for the best-fitting model to these measurements (Eqs. (3) and (4)) and the linear power spectrum. See Sect. 3 for details about this simulation and the fitting procedure. The dashed lines depict the results for the best-fitting model when setting D_NL = 1 after carrying out the fit; in other words, the behavior of the best-fitting model on linear scales. We can readily see that nonlinear growth isotropically increases the power with growing k, while peculiar velocities and thermal broadening suppress the power of parallel modes as k increases. On even smaller scales, pressure takes over and causes an isotropic suppression. Nonlinear growth modifies the perpendicular power relative to linear theory by 10% for scales as large as k = 0.5 Mpc⁻¹, indicating that small-scale corrections are important for most of the scales sampled by our simulations. Nevertheless, in Appendix A, we show that we can measure the two Lyman-α linear biases with percent precision from these simulations. Deviations from linear theory are less pronounced down to smaller scales for modes with μ ≃ 0.5 because nonlinear growth and the combination of peculiar velocities and thermal broadening tend to cancel each other out. As we can see, the parametric model achieves an average accuracy of 2% for k > 0.5 Mpc⁻¹, supporting the validity of Eq. (4) for capturing small-scale deviations from linear theory.

Fig. 1. Accuracy of the P3D model (see Eqs. (3) and (4)) in reproducing measurements from the CENTRAL simulation at z = 3. In the top panel, dotted and solid lines show the ratio of simulation measurements and model predictions relative to the linear power spectrum, respectively. Dashed lines do so for the linear part of the best-fitting model (DNL = 1). Line colors correspond to different μ wedges, and vertical dashed lines mark the minimum scale used for computing the best-fitting model, k = 5 Mpc−1. The bottom panel displays the relative difference between the best-fitting model and simulation measurements. The overall accuracy of the model is 2% on scales in which simulation measurements are not strongly affected by cosmic variance (k > 0.5 Mpc−1; see text).

On the largest scales, we find strong variations between consecutive k-bins for the same μ-wedge. Some of these oscillations are driven by differences in the average value of μ between consecutive bins due to the limited number of modes entering each bin on large scales. To ensure an accurate comparison between simulation measurements and model predictions, we individually evaluate the P_3D model for all the modes within each k − μ bin from our simulation boxes. We then calculate the mean of the resulting distribution and assign this mean value to the bin, thereby mirroring the approach used to compute P_3D measurements from the simulations. This process is also important on small scales, where the number of modes increases rapidly with k. Throughout this work, we adopt this approach to make predictions from the P_3D model.

Using this approach to generate theoretical predictions, the best-fitting model successfully captures most of the aforementioned large-scale oscillations. However, a fluctuation at k ≃ 0.25 Mpc⁻¹ in the 0 < μ < 0.25 wedge remains unaccounted for by the model. The negligible differences between model predictions and simulation measurements in the adjacent bins suggest that this fluctuation is likely due to cosmic variance. We assess the impact of this source of uncertainty on P_3D in Appendix A, finding that it can induce fluctuations of up to 10% on scales k < 0.5 Mpc⁻¹. As a result, cosmic variance limits our ability to evaluate the model’s performance on the largest scales shown. However, this does not necessarily reflect reduced model accuracy but rather highlights the limitations of using our simulations for validating the model. Proper validation on large scales would require either a larger simulation or multiple simulations with different initial conditions.

2.2. Input and output parameters

In addition to the density and velocity fields, the Lyman-α forest depends upon the ionization and thermal state of the IGM (e.g., McDonald 2003). Following Pedersen et al. (2021), we use six parameters to describe the dependency of this observable with cosmology and IGM physics:

• Amplitude and slope of the linear matter power spectrum on small scales. We define the amplitude (Δ_p²) and slope (n_p) as

  $$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{\mathrm{p}}^2(z)& ={(2{\pi }^2)}^{-1}{k}^3{P}_{\mathrm{lin}}(k_{\mathrm{p}},z),\hfill \end{array}$$(5)

  $$\begin{array}{cc}\hfill n_{\mathrm{p}}(z)& =(\mathrm{d}log{P}_{\mathrm{lin}}/\mathrm{d}logk){\mid }_{k=k_{\mathrm{p}}},\hfill \end{array}$$(6)

  where we use k_p = 0.7 Mpc⁻¹ as the pivot scale because it is at the center of the range of interest for DESI small-scale studies. These parameters capture multiple physical effects modifying the linear power spectrum on small scales (see Pedersen et al. 2021, for a detailed discussion), including cosmological parameters such as the amplitude (A_s) and slope (n_s) of the primordial power spectrum, the Hubble parameter, and the matter density (Ω_M), or ΛCDM extensions such as curvature and massive neutrinos. The advantage of using this parameterization rather than ΛCDM parameters is twofold. First, it reduces the dimensionality of the input to the cNF (see Sect. 4), which decreases the number of simulations required for precise training. Second, the resulting emulator has the potential for making precise predictions for variations in cosmological parameters and ΛCDM extensions not considered in the training set (Pedersen et al. 2021, 2023; Cabayol-Garcia et al. 2023). Note that we do not consider cosmological parameters capturing changes in the growth rate or expansion history because the Lyman-α forest probes cosmic times during which the universe is practically Einstein de-Sitter, and both vary very little with cosmology in this regime.

• Mean transmitted flux fraction. The mean transmitted flux fraction ($\overline{F}$) depends on the intensity of the cosmic ionizing background and evolves strongly with redshift. One of the advantages of using this parameter is that it encodes the majority of the redshift dependence of the signal, serving as a proxy for cosmic time.

• Amplitude and slope of the temperature-density relation. The thermal state of the IGM can be approximated by a power law on the densities probed by the Lyman-α forest (Lukić et al. 2015): T₀Δ_b^γ − 1, where Δ_b is the baryon overdensity, T₀ is the gas temperature at mean density, and γ − 1 is the slope of the relation. These parameters influence the ionization of the IGM, which is captured by $\overline{F}$, and the thermal motion of gas particles, which causes Doppler broadening that suppresses the parallel power. Instead of using T₀ as an emulator parameter, we follow Pedersen et al. (2021) and use the thermal broadening scale in comoving units. First, we express the thermal broadening in velocity units as ${\stackrel{\sim }{\sigma }}_{\mathrm{T}}=9.1{(T_0[\mathrm{K}]/{10}^4)}^{1/2}$, and then we convert it to comoving units, ${\sigma }_{\mathrm{T}}={\stackrel{\sim }{\sigma }}_{\mathrm{T}}(1+z)H^{-1}$.

• Pressure smoothing scale. Gas pressure supports baryons on small scales, leading to a strong isotropic power suppression in this regime that depends upon the entire thermal history of the gas (Gnedin & Hui 1998). We parameterize this effect using the pressure smoothing scale in units of comoving Mpc⁻¹, k_F (see Pedersen et al. 2021, for more details).

In summary, our cNF predicts the eight free parameters of the physically motivated model for Lyman-α clustering introduced by Eqs. (3) and (4), y = {b_δ,  b_η,  q₁,  q₂,  k_v,  a_v,  b_v,  k_p}, as a function of the aforementioned six parameters capturing the cosmological and IGM dependence of the Lyman-α forest, $\mathbf{x}=\{{\mathrm{\Delta }}_{\mathrm{p}}^2,n_{\mathrm{p}},\overline{F},{\sigma }_{\mathrm{T}},\gamma ,k_{\mathrm{F}}\}$.

3. Training and testing set

In this section, we describe how we generated the training and testing data for the cNF described in Sect. 4. In Sect. 3.1, we present a suite of cosmological hydrodynamical simulations from which we generated mock Lyman-α forest measurements, and we detail our approach for extracting P_3D and P_1D measurements from these simulations in Sect. 3.2. In Sect. 3.3, we compute the best-fitting parameters of the model introduced by Eqs. (3) and (4) to measurements of these statistics, and we evaluate the performance of the fits in Sect. 3.4.

3.1. Simulations

We extracted Lyman-α forest simulated measurements from a suite of simulations run with MP-GADGET⁴ (Feng et al. 2018; Bird et al. 2019), a massively scalable version of the cosmological structure formation code GADGET-3 (last described in Springel 2005). This suite was first presented and used in Pedersen et al. (2021); we briefly describe it next. Each simulation tracked the evolution of 768³ dark matter and baryon particles from z = 99 to z = 2 inside a box of L = 67.5 Mpc on a side, producing as output 11 snapshots uniformly spaced in redshift between z = 4.5 and 2. This configuration ensures convergence for P_1D measurements down to k_∥ = 4 Mpc⁻¹ (the smallest scale used in this work) at z = 2 and less than 10% errors for this scale at z = 4 (see Bolton et al. 2017, for more details). On the other hand, this configuration may cause non-negligible biases for P_3D at high redshift (Lukić et al. 2015).

Two realizations were run for each combination of cosmological and astrophysical parameters using the “fixed-and-paired” technique (Angulo & Pontzen 2016; Pontzen et al. 2016), which significantly reduces cosmic variance for multiple observables, including the Lyman-α forest (Villaescusa-Navarro et al. 2018; Anderson et al. 2019). The initial conditions were generated using the following configuration of MP-GENIC (Bird et al. 2020): initial displacements produced using the Zel’dovich approximation and baryons and dark matter initialized on an offset grid using species-specific transfer functions. Some studies have suggested that this configuration might lead to incorrect evolution of linear modes (Bird et al. 2020). However, in a recent study, Khan et al. (2024) showed that variations in the specific setting of MP-GENIC initial conditions have a minimal impact on P_1D measurements across the range of redshifts and scales used in this work.

To increase computational efficiency, the simulations utilized a simplified prescription for star formation that turns regions of baryon overdensity Δ_b > 1000 and temperature T < 10⁵ K into collisionless stars (e.g., Viel et al. 2004), implemented a spatially uniform ultraviolet background (Haardt & Madau 2012), and did not consider active galactic nuclei (AGN) feedback (e.g., Chabanier et al. 2020). These approximations are justified because we focus on emulating the Lyman-α forest in the absence of astrophysical contaminants like AGN feedback, damped Lyman-alpha absorbers (DLAs), or metal absorbers, and we will model these before comparing our predictions with observational measurements (e.g., McDonald et al. 2005; Palanque-Delabrouille et al. 2015, 2020).

In Sect. 4, we train a cNF using data from 30 fixed-and-paired simulations from the previous suite, covering combinations of cosmological and astrophysical parameters selected via a Latin hypercube sampling method (McKay et al. 1979). Hereafter, we refer to these as TRAINING simulations. The Latin hypercube spans the parameters {Δ_p²(z = 3), n_p(z = 3), z_H,  H_A,  H_S}, where we use z = 3 because it is approximately at the center of the range of interest for DESI studies (Ravoux et al. 2023; Karaçaylı et al. 2024), z_H is the midpoint of hydrogen reionization, and the last two parameters rescale the total photoheating rate ϵ₀ as ϵ = H_AΔ_b^H_Sϵ₀ (Oñorbe et al. 2017). Cosmological parameters were generated within the ranges Δ_p²(z = 3)∈[0.25,  0.45], n_p(z = 3)∈[−2.35,   − 2.25] by exploring values of the amplitude and slope of the primordial power spectrum within the intervals A_s ∈ [1.35,  2.71]×10⁻⁹ and n_s ∈ [0.92,  1.02]. Any other ΛCDM parameter was held fixed to values approximately following Planck Collaboration VI (2020): dimensionless Hubble parameter h = 0.67, physical cold dark matter density Ω_ch² = 0.12, and physical baryon density Ω_bh² = 0.022. As for the IGM parameters, these explored the ranges z_H ∈ [5.5,  15], H_A ∈ [0.5,  1.5], and H_S ∈ [0.5,  1.5]. All simulation pairs share the same set of initial Fourier phases, making their P_3D and P_1D subject to the same large-scale noise pattern.

We evaluated different aspects of the emulation strategy using six fixed-and-paired simulations with cosmological and astrophysical parameters not considered in the TRAINING simulations:

• The CENTRAL simulation uses cosmological and astrophysical parameters at the center of the TRAINING parameter space: A_s = 2.01 × 10⁻⁹, n_s = 0.97, z_H = 10.5, H_A = 1, and H_S = 1. We use this simulation for an out-of-sample test to evaluate the performance of FORESTFLOW under optimal conditions, recognizing that the accuracy of machine-learning models generally declines near the boundaries of the convex hull defined by the training set.

• The SEED simulation uses the same parameters as the CENTRAL simulation while considering a different distribution of initial Fourier phases. Given that all TRAINING simulations use the same initial Fourier phases, SEED is useful to evaluate the impact of cosmic variance in the training set on FORESTFLOW predictions.

• The GROWTH, NEUTRINOS, and CURVED simulations adopt the same values of Δ_p²(z = 3), n_p(z = 3), physical cold dark matter and baryonic densities, and astrophysical parameters as the CENTRAL simulation. However, the GROWTH simulation uses 10% larger Hubble parameter (h = 0.74) and 18% smaller matter density (Ω_M = 0.259) while using the same value of Ω_Mh² as the TRAINING simulations, the NEUTRINOS simulation includes massive neutrinos (∑m_ν = 0.3 eV), and the CURVED simulation considers an open universe (Ω_k = 0.03). The NEUTRINOS and CURVED simulations also modify the value of the cosmological constant while holding fixed h to compensate for the increase in the matter density and the addition of curvature, respectively. We use the testing simulations to evaluate the performance of the emulation strategy for cosmologies not included in the training set.

• The REIONISATION simulation uses the same cosmological parameters as the CENTRAL simulation while implementing a distinct helium ionization history relative to the CENTRAL and TRAINING simulations (Puchwein et al. 2019). The main difference between the ionization histories of these simulations is that the one implemented in the REIONISATION simulation peaks at a later time than the others, leading to a significantly different thermal history. The REIONISATION simulation therefore tests the performance of FORESTFLOW for thermal histories not considered in the TRAINING simulations.

3.2. Simulating Lyman-α forest data

To extract Lyman-α forest measurements from each simulation, we first selected one of the simulation axes as the line of sight and displace the simulation particles from real to redshift space along this axis. Then, we computed the transmitted flux fraction along 768² uniformly distributed line of sights along this axis using FSFE⁵ (Bird 2017); these lines of sight are commonly known as skewers. Following Springel (2005), we computed pressure forces using the density-entropy formulation of smoothed particle hydrodynamics (SPH) with a cubic spline kernel and 33 neighbors⁶. We set the resolution of the skewers to 0.05 Mpc, which is enough to resolve the thermal broadening and pressure scales, and spaced these by 0.09 Mpc in the transverse direction. We checked that P_3D and P_1D measurements within the range of interest (see Sect. 3.3) do not vary by increasing the line-of-sight resolution or the transverse sampling. We repeated all the previous steps for the three simulation axes to extract further cosmological information, as each simulation axis samples the velocity field in a different direction. Finally, we scaled the effective optical depth of the skewers to 0.90, 0.95, 1.05, and 1.10 times its original value, which is equivalent to running simulations with different UV background photoionization rates (see Lukić et al. 2015).

Using this data as input, we measured P_3D by first computing the three-dimensional Fourier transform of the skewers. Then, we took the average of the square norm of all modes within 20 logarithmically spaced bins in wavenumber k from the fundamental mode of the box, k_min = 2πL⁻¹ ≃ 0.09 Mpc⁻¹, to k_max = 40 Mpc⁻¹ and 16 linearly spaced bins in the cosine of the angle between Fourier modes and the line of sight from μ = 0 to 1. We measured P_1D by first computing the one-dimensional Fourier transform of each skewer without applying any binning, and then by taking the average of the square norm of all these Fourier transforms. The impact of cosmic variance on fixed-and-paired simulations is not straightforward (Maion et al. 2022), and thus we would ideally use multiple fixed-and-paired simulations with different initial distributions of Fourier phases to estimate the precision of P_3D and P_1D measurements from our simulations. However, such simulations are not available, and we instead relied on the comparison between two simulations with same configuration but different initial conditions (see Appendix A). We found that the impact of cosmic variance on P_3D and P_1D measurements can be as large as 10 and 1% on intermediate scales, respectively.

We measured P_3D and P_1D from the 30 TRAINING and the six test simulations, ending up with 2 (opposite Fourier phases) × 3 (simulation axes) × 11 (snapshots) × 5 (mean flux rescalings) = 330 measurements of each of these statistics per simulation. To reduce cosmic variance, we computed the average of measurements from different axes and phases of fixed-and-paired simulations, which decreased the number of measurements per simulation to 55. The training and testing sets of FORESTFLOW are thus comprised of 1650 and 330 Lyman-α power spectra measurements⁷, respectively.

3.3. Fitting the parametric model

To generate training and testing data for our emulator, we computed the best-fitting parameters of Eqs. (3) and (4) to measurements from the simulations described in Sect. 3.1. We fitted the model using P_3D measurements from k = 0.09 to 5 Mpc⁻¹ and P_1D measurements from k_∥ = 0.09 to 4 Mpc⁻¹. The size of our simulation boxes determines the largest scales used, while the smallest scales measured by DESI set the maximum wavenumbers (Ravoux et al. 2023; Karaçaylı et al. 2024). It is important to note that the two Lyman-α linear biases determine the large-scale behavior of P_3D (see Eq. (3)), and thus FORESTFLOW could make accurate predictions for P_3D on arbitrarily large (linear) scales as long as the Lyman-α linear biases are measured accurately.

We computed the best-fitting value of model parameters y = {b_δ,  b_η,  q₁,  q₂,  k_v,  a_v,  b_v,  k_p} to simulation measurements by minimizing the pseudo-χ²:

$$\begin{array}{cc}\hfill {\chi }^2(\mathbf{y})& ={\displaystyle \sum _i^{M_{3\mathrm{D}}}}{w}_{3\mathrm{D}}{[P_{3\mathrm{D}}^{\mathrm{data}}(k_i,{\mu }_i)-P_{3\mathrm{D}}^{\mathrm{model}}(k_i,{\mu }_i,\mathbf{y})]}^2\hfill \\ \hfill & +{\displaystyle \sum _i^{M_{1\mathrm{D}}}}{w}_{1\mathrm{D}}{[P_{1\mathrm{D}}^{\mathrm{data}}(k_{\Vert ,i})-P_{1\mathrm{D}}^{\mathrm{model}}(k_{\Vert ,i},\mathbf{y})]}^2,\hfill \end{array}$$(7)

where M_3D = 164 and M_1D = 42 were the number of P_3D and P_1D bins employed in the fit, respectively, the superscripts data and model refer to simulation measurements and model predictions, and w_3D and w_1D weighed the fit. We used the Nelder-Mead algorithm implemented in the routine MINIMIZE of SCIPY (Virtanen et al. 2020) to carry out the minimization⁸. The results of the fits are publicly accessible⁹.

Ideally, we would have used the covariance of P_3D and P_1D measurements to weigh the previous expression. However, estimating this covariance requires multiple realizations of the same simulation with different initial distributions of Fourier phases, and we do not have these simulations available. Instead, we disregarded correlations between P_3D and P_1D and weighed these by w_3D = N_3D(k, μ)/(1 + μ²)² and w_1D = α(1 + k_∥/k₀)², where N_3D is the number of modes in each k − μ bin and k₀ = 2 Mpc⁻¹. The terms involving N_3D, μ, and k₀ attempt to ensure an unbiased fit of P_3D and P_1D across the full range of scales used. The parameter α = 8000 controls the relative weight of P_3D and P_1D in the fit, and we set this value motivated by the different impact of cosmic variance on these statistics (see Appendix A).

We expect significant correlations between the best-fitting value of the parameters to measurements from relatively small simulation boxes. As shown by Arinyo-i-Prats et al. (2015), these correlations are especially significant for the parameters accounting for nonlinear growth of structure, q₁ and q₂. Givans et al. (2022) advocated for setting q₂ = 0 since this parameter is not necessary for describing P_3D at z = 2.8. However, we found non-zero values of this parameter indispensable for describing P_3D at redshifts below z = 2.5. This is not surprising since the gravitational evolution of density perturbations becomes increasingly more nonlinear as cosmic time progresses.

3.4. Accuracy of the parametric model

In the previous section, we computed the best-fitting parameters of the P_3D model to measurements from the TRAINING simulations. Two main sources of uncertainty can affect these fits: model inaccuracies and cosmic variance. The first relates to using a parametric model without enough flexibility to describe Lyman-α clustering accurately, while the second has to do with the limited size of the TRAINING simulations. The influence of cosmic variance on the training set is amplified because all the TRAINING simulations use the same initial distribution of Fourier phases, meaning all simulations are subject to the same large-scale noise. We study this source of uncertainty in Appendix A, where we compared the best-fitting models to the CENTRAL and SEED simulations, whose only difference is in their initial distribution of Fourier phases. We proceed to study model inaccuracies next.

In Fig. 2, we show the performance of the parametric model in reproducing P_3D and P_1D measurements from the 1650 snapshots of the TRAINING simulations. As discussed in Sect. 2.1, cosmic variance limits our ability to evaluate the accuracy of the model for P_3D on scales k < 0.5 Mpc⁻¹; therefore, we quote the model accuracy from k = 0.5 Mpc⁻¹ down to the smallest scale used in the fit, k = 5 Mpc⁻¹. In contrast, since cosmic variance has a much smaller impact on P_1D, we evaluate the model performance for this statistic using all scales considered in the fit (0.09 < k_‖ [Mpc⁻¹] < 4). We adopt the same approach when evaluating the performance of FORESTFLOW in Sect. 5. Under these considerations, the overall accuracy of the parametric model is 2.4 and 0.6% for P_3D and P_1D, respectively. Given that we estimate the accuracy of the parametric model using the TRAINING simulations, the previous numbers account for both the limited flexibility of such model and cosmic variance. As discussed in Appendix A, the impact of cosmic variance on measurements of P_3D and P_1D from these simulations is 1.3 and 0.5%, respectively. If we assume that cosmic variance and errors coming from the limited flexibility of the parametric model are uncorrelated and add in quadrature, the second are responsible for 2.0 and 0.3% errors on P_3D and P_1D, respectively.

Fig. 2. Accuracy of the parametric model (see Eqs. (3) and (4)) in reproducing P3D and P1D measurements from all the TRAINING simulations. Lines and shaded areas show the mean and standard deviation of the relative difference between simulation measurements from the 1650 snapshots of the TRAINING simulations and best-fitting models to these, respectively. The accuracy of the model in recovering P3D and P1D is 2.4 and 0.6%, respectively, on scales not strongly affected by cosmic variance.

4. Emulator

In this section, we use a cNF to predict the two Lyman-α linear biases and six parameters describing small-scale deviations of P_3D from linear theory as a function of cosmology and IGM physics. We detail the architecture and implementation of this emulator in Sects. 4.1 and 4.2, respectively.

4.1. Conditional normalizing flows

Normalizing flows (NFs; Jimenez Rezende & Mohamed 2015) are a class of machine-learning generative models designed to predict complex distributions by applying a sequence of bijective mappings to simple base distributions. A natural extension to this framework is conditional NFs (cNFs; Winkler et al. 2019; Papamakarios et al. 2019), a type of NFs that condition the mapping between the base and target distributions on a series of input variables. Given an input x ∈ X and target y ∈ Y, cNFs predict the conditional distribution p_Y|X(y|x) by applying a parametric, bijective mapping f_ϕ : Y × X → Z to a base distribution p_Z(z) as follows

$$\begin{array}{c}\hfill p_{Y|X}(\mathbf{y}|\mathbf{x})=p_Z(f_{\varphi }(\mathbf{y},\mathbf{x})|\mathbf{x})\left|\frac{\partial f_{\varphi }(\mathbf{y},\mathbf{x})}{\partial \mathbf{y}}\right|,\end{array}$$(8)

where ϕ are the parameters of the mapping, while the last term of the previous equation is the Jacobian determinant of the mapping. In our cNF, the input is given by the parameters capturing the dependence of the Lyman-α forest on cosmology and IGM physics, $\mathbf{x}=\{{\mathrm{\Delta }}_{\mathrm{p}}^2,n_{\mathrm{p}},\overline{F},{\sigma }_{\mathrm{T}},\gamma ,k_{\mathrm{F}}\}$, the target by the parameters of the P_3D model, y = {b_δ,  b_η,  q₁,  q₂,  k_v,  a_v,  b_v,  k_p}, and the base distribution is an eight-dimensional Normal distribution N⁸(0, 1), where the dimension is determined by the number of P_3D model parameters.

Once trained, cNFs are a generative process from x to y. In our implementation, we start by randomly sampling from the base distribution, and then we pass this realization through a sequence of mappings conditioned on a particular combination of cosmology and IGM parameters, $\stackrel{\sim }{\mathbf{y}}=f_{\varphi }^{-1}(p_Z(\mathbf{z}),\mathbf{x})$, ending up with a prediction for the value of the P_3D model parameters. Repeating this process multiple times, the emulator yields a distribution of P_3D parameters $p_{\stackrel{\sim }{Y}|X}$ that, for a sufficiently large number of samples, approaches the target distribution p_Y|X. The breadth of this distribution captures uncertainties arising from the limited size of the training set. Finally, outside the cNF, we use each combination of P_3D parameters to evaluate Eqs. (3) and (1), obtaining predictions and uncertainties for P_3D and P_1D.

The main challenge when using cNFs is finding the mapping between the target and the base distribution, typically done using an N-layer neural network with bijective layers. This process runs in reverse relative to the generating process: we start by applying the mapping f_ϕ to the target data y conditioned on the input x, yielding z. Then, we optimize the model parameters by minimizing the loss function

$$\begin{array}{c}\hfill \mathcal{L}=\frac{1}{2}{\displaystyle \sum {\mathbf{z}}^2}-log\left|\frac{\delta f_{\varphi }(\mathbf{y},\mathbf{x})}{\delta \mathbf{y}}\right|.\end{array}$$(9)

We carried out this optimization process using stochastic gradient descent applied to minibatches, a methodology commonly employed for training neural networks.

4.2. Implementation

Neural Autoregressive Flows (Huang et al. 2018) use a series of invertible univariate operations to build a bijective transformation between a conditional distribution and a base distribution. In FORESTFLOW, we created a bijective mapping between the best-fitting parameters of the P_3D model and an eight-dimensional Normal distribution by applying N_ACB = 12 consecutive Affine-Coupling Block (ACB; Dinh et al. 2016) conditioned on cosmology and IGM physics. The transformation goes from the best-fitting parameters of the P_3D model to the base distribution when training the model, and in the opposite direction when evaluating it.

Each ACB conducts a series of operations $g_{i,{\stackrel{\sim }{\varphi }}_i}$ on its input data w_i, with i going from 1 to N_ACB and ${\stackrel{\sim }{\varphi }}_i$ standing for the parameters of the transformation. First, it splits the input data into two subsamples with approximately the same number of elements, w′_i and w″_i. Then, it applies an affine transformation to the first subsample w′_i

$$\begin{array}{c}\hfill T({\mathbf{w}}_i^{\prime })={\alpha }_i{\mathbf{w}}_i^{\prime }+{\beta }_i,\end{array}$$(10)

where α_i and β_i are neural networks with a single hidden layer of 128 neuron units. Third, the ACB merges the output from the affine transformation and the unchanged subsample, and then it applies a permutation layer to randomly rearrange these elements, obtaining ${\stackrel{\sim }{\mathbf{w}}}_i$. Fourth, the ACB applies an affine transformation to this sample, $\stackrel{\sim }{T}({\stackrel{\sim }{\mathbf{w}}}_i)$. The first and second affine transformations involve a subset of the training set and the entire training set, respectively, enabling the model to capture local and global features.

In Fig. 3, we show the architecture of the cNF. The blue arrow indicates the training direction, while the green arrow depicts the emulation direction. In the training direction, the input to the first ACB, u₁ = w₁, is a 1650-dimensional array composed of 14-dimensional vectors, where 1650 is the number of simulation snapshots in the training set. Each vector includes the eight best-fitting P_3D model parameters to each snapshot and the six parameters describing the cosmology and IGM physics of this snapshot. The input to the i ACB, u_i, is a 1650-array containing 14-dimensional vectors with the output of the i − 1 ACB and, once again, the six parameters describing the cosmology and IGM physics of each snapshot. Each ACB applies a transformation $f_{i,{\varphi }_i}=g_{i,{\stackrel{\sim }{\varphi }}_i}$, and the consecutive application of all ACBs results in the mapping between the target and the base distributions z = f_ϕ(y, x), where $f_{\varphi }={\displaystyle \prod _{i=1}^{N_{\mathrm{ACB}}}}{f}_{i,{\varphi }_i}$.

Fig. 3. Architecture of the Lyman-α forest clustering emulator. The blue arrow indicates the training direction, where the cNF optimizes a bijective mapping between the best-fitting parameters of the P3D model to measurements from the TRAINING simulations and an eight-dimensional Normal distribution. The mapping is conditioned on cosmology and IGM physics, and performed using 12 consecutive affine coupling blocks. The green arrow denotes the emulation direction, where the cNF applies the inverse of the mapping to random samples from the base distribution to predict the value of the P3D model parameters. Outside the cNF, FORESTFLOW introduces these parameters in Eq. (3) and (1) to obtain predictions for P3D and P1D, respectively.

In the emulation direction, the input to the first ACB, v₁ = w₁, is a 14-dimensional vector containing random draws from an eight-dimensional Normal distribution and the six parameters describing the cosmology and IGM physics for which we want to obtain predictions. As in the training direction, the input to each subsequent ACB relies on the output from the previous ACB, each conditioned on cosmology and IGM physics. The ACBs apply the transformations $f_{i,{\varphi }_i}^{-1}=g_{i,{\stackrel{\sim }{\varphi }}_i}$, which are the inverse of the corresponding transformations in the training direction, f_i, ϕi. The cNF makes predictions for P_3D model parameters by applying the composition of the inverse of all ACBs to random samples from the base distribution, $\stackrel{\sim }{\mathbf{y}}=f_{\varphi }^{-1}(p_Z(\mathbf{z}),\mathbf{x})$, where $f_{\varphi }^{-1}={\displaystyle \prod _{i=1}^{N_{\mathrm{ACB}}}}{f}_{i,{\varphi }_i}^{-1}$.

We implemented the emulator within the FreIA framework (Ardizzone et al. 2018-2022), which uses PyTorch (Ansel et al. 2024) in the backend. We trained it by minimizing Eq. (9) using an Adam optimizer (Kingma & Ba 2014) for 300 epochs with an initial learning rate of 10⁻³. We used the Optuna framework (Akiba et al. 2019) to select the number of ACBs and epochs, as well as the value of the learning rate. First, Optuna trains our cNF for a particular combination of these hyperparameters. Then, it computes the average value of Eq. (7) for all simulations in the training set. After that, depending on the goodness of the fit to P_3D and P_1D measurements, Optuna selects a new value of the hyperparameters. We iterated with Optuna 50 times through a hyperparameter grid, selecting the hyperparameters that yield the highest accuracy. We checked that the performance of the cNF depended weakly on small variations in the value of the hyperparameters.

5. Performance of FORESTFLOW

In Sect. 5.1, we analyze the performance of FORESTFLOW across the parameter space of the training set. Then, in Sect. 5.2, we test its accuracy using simulations with cosmologies and IGM models that are not part of the training set. All performance evaluations are conducted using the simulations described in Sect. 3.1, which were run using the same code and resolution. Before employing FORESTFLOW for cosmological inference, it will be crucial to validate it against large, high-resolution simulations produced with alternative codes. We defer this task to future work.

5.1. Cosmologies and IGM histories in the training set

In this section, we evaluate the performance of FORESTFLOW in recovering the two Lyman-α linear biases, which determine the behavior of P_3D on linear scales, as well as P_3D and P_1D measurements from simulations on the intervals 0.5 < k [Mpc⁻¹]< 5 and 0.09 < k_‖ [Mpc⁻¹] < 4, respectively. These are the ranges of scales used when fitting the parametric model in Sect. 3 that are not strongly affected by cosmic variance (see Sect. 2.1). We begin by assessing the accuracy of FORESTFLOW at the center of the training set, where machine-learning methods typically perform best, and then extend our evaluation across the entire input parameter space.

In Fig. 4, we compare measurements of P_3D and P_1D from the CENTRAL simulation at z = 3 with FORESTFLOW predictions. Dotted lines show simulation measurements, while solid lines and shaded areas display the average and 68% credible interval of FORESTFLOW predictions, respectively. We characterize the accuracy of the credible intervals in Appendix B. As we can see, FORESTFLOW captures the amplitude and scale-dependence of P_3D and P_1D precisely. In Fig. 5, we present the relative difference between P_3D and P_1D measurements from the CENTRAL simulation and FORESTFLOW predictions as a function of redshift. The model’s accuracy remains consistent for redshifts above z = 2 but shows a slight decline at this redshift. This is likely because z = 2 is the lowest redshift included in the training set and is therefore near the boundary of the convex hull defined by the training data.

Fig. 4. Accuracy of FORESTFLOW in recovering P3D and P1D measurements from the CENTRAL simulation at z = 3. Dotted lines show measurements from simulations, solid lines and shaded areas display the average and 68% credible interval of FORESTFLOW predictions, respectively, and vertical dashed lines indicate the minimum scales considered for computing the training data for the cNF. The overall performance of FORESTFLOW in recovering P3D is 2.0% on scales not strongly affected by cosmic variance and 0.6% for P1D.

Fig. 5. Accuracy of FORESTFLOW in recovering P3D and P1D measurements from the CENTRAL simulation as a function of redshift. The upper four panels show the results for P3D across different μ bins, while the bottom panel displays the results for P1D. Each color represents a different redshift. The model’s accuracy remains consistent for redshifts above z = 2 but exhibits a slight decline at this redshift.

To better characterize the performance of FORESTFLOW, we compute the average accuracy of FORESTFLOW in recovering measurements from CENTRAL across redshift. We find that it is 1.2 and 0.3% for b_δ and b_η, respectively, which translates into 1.1 and 1.2% for perpendicular and parallel P_3D modes on linear scales, and 2.6 and 0.8% for P_3D and P_1D. Note that cosmic variance hinders our ability to test the performance of the model; however, this does not necessarily indicate a decrease in the model’s accuracy for P_3D on the largest scales sampled by our simulation.

We expect the efficiency of FORESTFLOW to decrease away from the center of the input space. We could assess its performance across the parameter space using the training simulations; however, the cNF has been optimized for these points, which introduces the risk of overfitting. Overfitting could result in high precision for these specific points of the parameter space but not for others nearby. As a result, we would ideally evaluate the performance of FORESTFLOW using multiple test simulations covering the entire parameter space, but such simulations are unavailable. Instead, we conduct leave-one-out tests, which are widely used to assess the performance of an emulator when the number of training points is insufficient for out-of-sample tests (e.g., Hastie et al. 2001). In a leave-one-out test, we optimize a cNF after removing a subsample from the training set; for example, all measurements from one of the TRAINING simulations. We then check the accuracy of FORESTFLOW for the new cNF using the subsample held back. The rationale is that the new emulator should closely approximate the original emulator everywhere in the parameter space except near the excluded simulation, and more importantly, there is no risk of overfitting. By repeating this process for other subsamples, we can estimate the performance of FORESTFLOW across the parameter space. Since each cNF is trained without using the entire dataset, leave-one-out tests provide a lower bound on FORESTFLOW performance. Additionally, leave-one-out tests may require extrapolating the predictions from the cNF, and it is widely known that machine-learning methods do not extrapolate well.

In the top panels of Fig. 6, lines and shaded areas display the average and standard deviation of 30 leave-simulation-out tests. Each test requires optimizing a cNF with 29 distinct TRAINING simulations, and then using the remaining simulation as the validation. Each panel shows the results for a different redshift, and we check that the results are similar for redshifts not shown. As we can see, the large-scale noise is similar for all TRAINING simulations; this is because they use the same initial distribution of Fourier phases. The overall performance of FORESTFLOW in recovering b_δ and b_η is 1.0 and 3.1%, respectively, which translates into 2.0 and 2.9% for perpendicular and parallel P_3D modes on linear scales, and 3.4 and 1.8% for P_3D and P_1D.

Fig. 6. Accuracy of FORESTFLOW across the input parameter space estimated via leave-simulation-out (top panels) and leave-redshift-out tests (bottom panels). Top panels. Each leave-simulation-out test involves training one independent emulator with measurements from 29 distinct simulations, and then using the measurements from the remaining simulation as the validation set. Lines and shaded areas show the average and standard deviation of 30 leave-simulation-out tests, and each panel shows the results for a different redshift. Bottom panels. Leave-redshift-out tests require optimizing one emulator with all measurements but the ones at a particular redshift, and then using measurements from this redshift as validation. Each panel shows the results of a different test.

In Table 1, we gather the accuracy of FORESTFLOW at the center and across the parameter space, as well as the expected level of uncertainties due to cosmic variance and the limited flexibility of the P_3D model. Due to the limited size of our simulations, the maximum levels of accuracy we can test for P_3D and P_1D are 1.3 and 0.5% (see Appendix A), respectively. These levels would decrease by evaluating the accuracy of FORESTFLOW using bigger simulations with the same resolution. On the other hand, the combined impact of impact of cosmic variance on the training data and the limited flexibility of the P_3D model are 2.4 and 0.6% for P_3D and P_1D, respectively, which is 1.1 and 0.1% worse than the minimum accuracy we can test for these statistics. At the center of the parameter space, the accuracy of FORESTFLOW for P_3D and P_1D is only 0.2% worse than the previous levels, letting us conclude that the primary factors limiting the performance of FORESTFLOW at the center of the parameter space are the size of the training simulations and model inaccuracies.

The efficiency of FORESTFLOW across the parameter space is 1.2 and 1.0% worse than at the center for P_3D and P_1D, respectively. Consequently, its accuracy would likely improve by increasing the number of training simulations. However, leave-one-out tests significantly underestimate the performance of an emulator at the edges of the training set, especially for a small number of simulations, because it often requires extrapolating the emulator’s predictions. We can thus conclude that the quality of the training data, the accuracy of the model, and the number of training simulations have a similar impact on the performance of FORESTFLOW. Given that leave-simulation-out tests tend to provide a lower performance bound, we conclude that the overall accuracy of FORESTFLOW in predicting P_3D from linear scales to k = 5 Mpc⁻¹ is approximately 3%, and ≃1.5% for P_1D down to k_∥ = 4 Mpc⁻¹.

As discussed in Sect. 2.2, FORESTFLOW does not use as input “traditional” cosmological parameters such as Ω_m, A_s, or H₀. Instead, it uses a set of parameters measured from the outputs of individual simulation snapshots. This strategy enables training FORESTFLOW without specifying the input redshift and making predictions for redshifts not present in the training set. To test this assumption, we carry out two leave-redshift-out tests. The first involves optimizing one emulator with all TRAINING measurements but the ones at z = 2.5, and then validating it with data from this redshift. For the second, we follow the same approach but using measurements at z = 3.5. We display the results of these tests in the bottom panels of Fig. 6. The performance of FORESTFLOW is similar for leave-redshift-out and leave-simulation-out tests, validating the approach mentioned above. We find similar results for leave-redshift-out tests at other redshifts.

5.2. Other cosmologies and IGM histories

In Fig. 7, we examine the accuracy of FORESTFLOW in reproducing P_3D and P_1D measurements from simulations not included in the training set. Lines indicate the redshift average of the relative difference between model predictions and simulation measurements. The first two rows show the results for the CENTRAL and SEED simulations, whose only difference is their initial distribution of phases. Consequently, the predictions of FORESTFLOW are the same for the two simulations. As we can see, these simulations present a different large-scale pattern of fluctuations, signaling that are caused by cosmic variance. Once we ignore these, we find that the performance of FORESTFLOW is practically the same for the two simulations. We can thus conclude that FORESTFLOW predictions are largely insensitive to the impact of cosmic variance on the training set.

Fig. 7. Performance of FORESTFLOW in recovering P3D and P1D for test simulations not included in the training set. Lines and shaded areas display the average and standard deviation of the results for 11 snapshots between z = 2 and 4.5, respectively. From top to bottom, the rows show the results for the CENTRAL, SEED, GROWTH, NEUTRINOS, CURVED, and REIONISATION simulations, where the CENTRAL and SEED simulations are at the center of the input parameter space and employ the same and different initial distribution of Fourier phases as the training simulations, respectively, the GROWTH and REIONISATION simulations use a different growth and reionization history relative to those used by the TRAINING simulations, and the NEUTRINOS and CURVED simulations consider massive neutrinos and curvature. The efficiency of FORESTFLOW is approximately the same for all simulations.

In the third, fourth, and fifth rows of Fig. 7, we use the GROWTH, NEUTRINOS, and CURVED simulations to evaluate the accuracy of FORESTFLOW for three different scenarios not contemplated in the training set: different growth history, massive neutrinos, and curvature. As we can see, the performance of FORESTFLOW for all these simulations is approximately the same as for the CENTRAL simulation. These results support that using the small-scale amplitude and slope of the linear power spectrum to capture cosmological information enables setting precise constraints on growth histories and ΛCDM extensions not included in the training set (see also Pedersen et al. 2021, 2023; Cabayol-Garcia et al. 2023).

In the last row of Fig. 7, we examine the accuracy of FORESTFLOW for the REIONISATION simulation, which employs a He II reionization history significantly different from those used by the TRAINING simulations. The performance of FORESTFLOW for this and the CENTRAL simulation is similar, which is noteworthy given that the performance of P_1D emulators for the REIONISATION simulation is significantly worse than for the CENTRAL simulation (Cabayol-Garcia et al. 2023). The outstanding performance of FORESTFLOW is likely because the relationship between IGM physics and the parameters of the P_3D model is more straightforward than with P_1D variations.

6. Discussion

Cosmological analyses of the Lyman-α forest come in two flavors: one-dimensional studies focused on small, nonlinear scales and three-dimensional analyses of large, linear scales. With FORESTFLOW, we can now consistently model Lyman-α correlations from nonlinear to linear scales, enabling a variety of promising analyses that we discuss next.

6.1. Connecting large-scale biases with small-scale physics

Small-scale Lyman-α analyses use emulators to predict P_1D as a function of cosmology and IGM physics (e.g., Cabayol-Garcia et al. 2023), while large-scale analyses use linear or perturbation theory models to predict ξ_3D together with Lyman-α linear bias parameters that need to be marginalized over. FORESTFLOW provides a relationship between IGM physics and linear biases, enabling the use of P_1D studies to inform three-dimensional analyses and vice versa.

We could use FORESTFLOW to set constraints on b_δ and b_η by fitting P_1D measurements, and then use these constraints as priors in three-dimensional studies. As a result, we would break degeneracies between Lyman-α linear bias parameters and cosmology, allowing us to measure the amplitude of linear density and velocity fluctuations, σ₈(z) and fσ₈(z), rather than b_δσ₈ and b_ηfσ₈ like in traditional Lyman-α forest analyses. To illustrate this application, we proceed to compare measurements of b_δ and β ≡ b_δ⁻¹b_ηf from BAO analyses with FORESTFLOW predictions for these parameters based on small-scale P_1D analyses. The analysis of BAO in the Lyman-α forest from the first data release of DESI yields b_δ = −0.108 ± 0.005 and β = 1.74 ± 0.09 at z = 2.33 (DESI Collaboration 2024). On the other hand, FORESTFLOW predicts b_δ = −0.118 and β = 1.57 at z = 2.33 for a Planck cosmology when using as input the best-fitting constraints on IGM parameters from Table 4 of Walther et al. (2019), which were derived from high-resolution P_1D measurements. The constraints on IGM parameters were derived using a P_1D emulator trained on a suite of simulations with the same input cosmology and possibly slightly different definitions of IGM parameters relative to those used in this work. Nonetheless, FORESTFLOW predictions and DESI measurements agree at the two sigma level, encouraging this new type of study.

In the left panels of Fig. 8, we display FORESTFLOW predictions for the response of the Lyman-α linear biases and β to variations in cosmology and IGM physics. The response of b_δ to these changes is strong and has a different redshift dependence for cosmology and IGM parameters; therefore, we could use FORESTFLOW to analyze P_3D measurements from different redshifts to further break degeneracies between b_δ and σ₈. On the other hand, the response of b_η to these changes is weak, and it is thus challenging to use this approach to break degeneracies between b_η and fσ₈. Note that the response of the Lyman-α linear biases and β to A_s variations broadly agrees with measurements from simulations run while only varying σ₈ (Arinyo-i-Prats et al. 2015).

Fig. 8. Response of Lyman-α clustering to variations in cosmology and IGM physics according to FORESTFLOW. The top, middle, and bottom panels of the left column show the results for β, bδ, and bη, respectively, while those of the right column do so for the perpendicular modes of P3D, the parallel modes of P3D, and P1D. Blue, orange, and red lines show the response of the previous quantities to a 5% increase in As, ΩMh2, and σT, respectively, while green lines do so for a 1% increase in $\overline{F}$.

Similarly, we could use measurements of linear bias parameters from three-dimensional analyses (du Mas des Bourboux et al. 2020; DESI Collaboration 2024) to make predictions for IGM parameters, which could be used in P_1D studies to break degeneracies between cosmology and IGM physics. In the right panels of Fig. 8, we display FORESTFLOW predictions for the response of P_3D and P_1D to variations in cosmology and IGM physics. As we can see, the response of P_1D to A_s and $\overline{F}$ variations is largely scale-independent down to k_∥ = 1 Mpc⁻¹ where many other effects are at play, and thus these two parameters are largely degenerated. On the other hand, this is not the case for P_3D; consequently, we could use information from P_3D analyses to break degeneracies in P_1D studies. Note that the response of P_3D and P_1D to A_s, $\overline{F}$, and σ_T variations broadly agrees with measurements from simulations run varying only one of these parameters at a time (McDonald 2003; McDonald et al. 2005).

We also observe that P_3D and P_1D respond significantly to variations in Ω_Mh², suggesting that the Lyman-α clustering is highly sensitive to the expansion and growth history. However, we find that variations in A_s and n_s can absorb the changes in P_3D and P_1D to the 2% level, and completely do so for P_3D at the pivot scale of the cosmological parameters of FORESTFLOW, k_p = 0.7 Mpc⁻¹. Furthermore, Planck measured Ω_Mh² with 0.8% precision Planck Collaboration VI (2020), and A_s and n_s absorb 1% variations in Ω_Mh² to the ≃0.4% level. This result supports the approach of not considering any cosmological parameter related to variations in the expansion of growth history as input for FORESTFLOW (see also Sect. 5.2).

6.2. Alcock-Paczyński on mildly nonlinear scales

Thanks to the increasing precision of galaxy surveys, there is a growing interest in extracting cosmological information from increasingly smaller scales in three-dimensional analyses. An avenue to do so is to analyze anisotropies in the correlation function Alcock & Paczynski (AP test; 1979), first proposed in the context of the Lyman-α forest by McDonald & Miralda-Escudé (1999) and Hui et al. (1999). Recently, Cuceu et al. (2023) followed this approach to analyze Lyman-α forest measurements from the Sloan Digital Sky Survey (SDSS) data release 16 (DR16; Ahumada et al. 2020), yielding constraints on some cosmological parameters a factor of two tighter than those from BAO-only analyses.

This study modeled three-dimensional correlations using linear theory, which restricted the range of scales analyzed to those larger than 25 h⁻¹ Mpc. We could significantly extend the range of scales used in this type of analysis by modeling three-dimensional correlations using FORESTFLOW. As a result, the constraining power of AP analyses would be much larger. Furthermore, we could use FORESTFLOW to extract information from P_1D analyses to reduce degeneracies between cosmology and the parameters describing ξ_3D (see Sect. 6.1).

6.3. Extending 3D analyses to the smallest scales

The ultimate goal of FORESTFLOW is to perform a joint analysis of one- and three-dimensional measurements from small to large scales. An interesting approach to do so is to measure the Lyman-α forest cross-spectrum (P_×; e.g., Hui et al. 1999; Font-Ribera et al. 2018), which captures the correlation between one-dimensional Fourier modes from two neighboring quasars separated by a transverse separation (r_⊥). We can model this statistic by taking the inverse Fourier transform of P_3D only along the perpendicular directions

$$\begin{array}{cc}\hfill P_{\times }(k_{\Vert },r_{\perp })& \equiv \frac{1}{{(2\pi )}^2}{\displaystyle \int \mathrm{d}}{\mathit{k}}_{\perp }{e}^{i{\mathit{k}}_{\perp }·{\mathit{r}}_{\perp }}{P}_{3\mathrm{D}}(k_{\Vert },k_{\perp })\hfill \\ \hfill & =\frac{1}{2\pi }{\displaystyle {\int }_0^{\infty }}\mathrm{d}{k}_{\perp }{k}_{\perp }{J}_0(k_{\perp }{r}_{\perp })P_{3\mathrm{D}}(k_{\Vert },k_{\perp }).\hfill \end{array}$$(11)

Comparing this equation with Eq. (1), it becomes clear that P_1D is a special case of P_×, corresponding to the limit where the transverse separation is zero.

In Sect. 3.3, we optimize the P_3D model to describe measurements of P_3D and P_1D from the TRAINING simulations. Then, in Sect. 4, we use the distribution of best-fitting parameters as the training set for FORESTFLOW, which predicts the value of P_3D model parameters as a function of cosmology and IGM physics. Even though neither the best-fitting model nor FORESTFLOW use P_× for their optimization, we can make predictions of P_× for the two. To do so, we first estimate P_3D using the value of the model parameters using Eq. (3), and then we integrate it using Eq. (11). We carry out the integration using the fast Hankel transform algorithm FFTlog (Hamilton 2000) implemented in the hankl package (Karamanis & Beutler 2021).

We use P_× measurements from the simulations described in Sect. 3.1 to evaluate the accuracy of FORESTFLOW for this statistic. We first define four bins in r_⊥, the transverse separation between skewers in configuration space, with edges 0.13, 0.32, 0.80, 2, and 6 Mpc. Then, we measure P_× using all pairs of skewers with r_⊥ separation within the previous bins

$$\begin{array}{c}\hfill P_{\times }(r_{\perp },k_{\Vert })=\langle \mathfrak{R}[\stackrel{\sim }{{\delta }_{\mathrm{i}}}(k_{\Vert }){\stackrel{\sim }{{\delta }_{\mathrm{j}}}}^{\ast }(k_{\Vert })]\rangle \end{array}$$(12)

where ${\stackrel{\sim }{\delta }}_{\mathrm{i}}$ and ${\stackrel{\sim }{\delta }}_{\mathrm{j}}^{\ast }$ stand for the Fourier transform of a skewer i and the complex conjugate of its partner j, respectively, the average ⟨⟩ includes all possible pairs in the bin without repetition or permutation, and ℜ indicates that we only use the real part of the expression between brackets because the average of the imaginary part is zero. The r_⊥ on the left-hand side denotes the effective center of the bin, accounting for the skewed distribution of r_⊥ within each bin: the number of skewers separated by a small distance dr_⊥ is proportional to r_⊥, and therefore the effective center is not at the halfway point. To compute P_× at the effective center, we perform the integration using ten sub-bins within each r_⊥ bin and calculate the average of these weighed by r_⊥.

In Fig. 9, we study the performance of FORESTFLOW in reproducing P_× measurements from the CENTRAL simulation at z = 3. Dots display simulation measurements, dashed lines the best-fitting model to P_3D and P_1D measurements from this simulation, and the solid lines FORESTFLOW predictions. As we can see, P_× decreases as the r_⊥ separation increases; this is because more distant sightlines are sampling increasingly uncorrelated regions. In the middle panel, we examine the accuracy of the best-fitting model in describing simulation measurements, finding that it is better than 10% throughout all the scales shown. The performance of the model improves for smaller r_⊥ separations. This is likely because the fit’s likelihood function (Eq. (7)) considers P_1D, which is equivalent to P_× at r_⊥ = 0 separation, but not P_×. The bottom panel illustrates the performance of FORESTFLOW relative to the best-fitting model, providing an approximate assessment of its ability to reproduce the training data. FORESTFLOW achieves an accuracy better than 5.

Fig. 9. Accuracy of the parametric model and FORESTFLOW in describing P× measurements from the CENTRAL simulation at z = 3. Dots show simulation measurements, dashed lines depict predictions from the best-fitting parametric model to P3D and P1D measurements, and solid lines and shaded areas display the average and 68% credible interval of FORESTFLOW predictions. The color of the lines indicates the results for different bins in transverse separation r⊥. The middle panel shows the residual between simulation measurements and the best-fitting parametric model, while the bottom panel displays the residual between predictions from the parametric model and FORESTFLOW. The performance of FORESTFLOW in reproducing simulation measurements is similar to that of the best-fitting model.

Future studies could use FORESTFLOW for extracting constraints on cosmology and IGM physics from the analysis of P_× measurements (e.g., Abdul Karim et al. 2024). Nevertheless, as with P_1D, these analyses would also require modeling multiple systematics affecting Lyman-α measurements such as damped Lyman-α systems, metal line contamination, and AGN feedback.

7. Conclusions

We present FORESTFLOW, a novel framework for predicting Lyman-α clustering from linear and nonlinear scales as a function of cosmology and IGM physics. FORESTFLOW employs conditional normalizing flows to emulate the eight parameters of a physically motivated model for Lyman-α clustering: the two linear Lyman-α biases (b_δ and b_η) and six parameters that capture small-scale deviations of the three-dimensional flux power spectrum (P_3D) from linear theory. By combining this model with a Boltzmann solver, FORESTFLOW predicts P_3D and any derived statistics, including the two-point correlation function (ξ_3D, the primary statistic for large-scale analyses), the one-dimensional Lyman-α flux power spectrum (P_1D, central to small-scale studies), and the cross-spectrum (P_×, a promising tool for full-scale analyses).

We trained the emulator using the best-fitting parameters of the physically motivated model to measurements from a suite of 30 fixed-and-paired cosmological hydrodynamical simulations spanning 11 redshifts equally spaced between z = 2 and 4.5 (Pedersen et al. 2021). Despite the moderate size of these simulations, FORESTFLOW achieves an accuracy of 3% for P_3D from linear scales to k = 5 Mpc⁻¹ and 1.5% for P_1D down to k_∥ = 4 Mpc⁻¹. FORESTFLOW uncertainties arise from three comparable sources: the size of the training simulations, their number, and the limited flexibility of the physically motivated model. However, we only evaluate the accuracy of FORESTFLOW using simulations from our suite, and its performance may vary for higher-resolution simulations or those generated with other codes.

FORESTFLOW demonstrates comparable performance for two extensions to the ΛCDM model – massive neutrinos and curvature – and ionization histories not included in the training set. This generalization is made possible by emulating the parameters of the physically motivated model as a function of the small-scale amplitude and slope of the linear power spectrum, the mean transmitted flux fraction, the amplitude and slope of the temperature-density relation, and the pressure smoothing scale (see Pedersen et al. 2023; Cabayol-Garcia et al. 2023).

The release of FORESTFLOW is particularly timely for Lyman-α forest analyses with the ongoing Dark Energy Spectroscopic Instrument (DESI) survey. As discussed in Sect. 6, FORESTFLOW enables a range of novel analyses with DESI data, including linking large- and small-scale studies and extending three-dimensional analyses to smaller scales. However, before applying FORESTFLOW to cosmological inference, it will be necessary to model the effects of various astrophysical processes on Lyman-α clustering, including metal contamination, damped Lyman-α systems, and AGN feedback.

Data availability

FORESTFLOW and all the notebooks used to generate the plots of this paper can be found in Github¹⁰, as well as all data points shown in the published graphs. The simulations utilized for training and testing the emulator are also publicly accessible¹¹.

Acknowledgments

We thank Eric Armengaud, Roger de Belsunce, Andrei Cuceu, Vid Irsic, Ignasi Pérez-Ràfols, and Michael Walther for their useful comments and helpful discussions. We also thank the referee for his detailed comments and suggestions. JCM, LCG, ML, and AFR acknowledge support from the European Union (ERC Consolidator Grant, COSMO-LYA, grant agreement 101044612). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. AFR acknowledges financial support from the Spanish Ministry of Science and Innovation under the Ramon y Cajal program (RYC-2018-025210) and the PGC2021-123012NB-C41 project. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya. The analysis has been performed at Port d’Informació Científica (PIC); we acknowledge the support provided by PIC in granting us access to their computing resources. This material is based upon work supported by the U.S. Department of Energy (DOE), Office of Science, Office of High-Energy Physics, under Contract No. DE–AC02–05CH11231, and by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract. Additional support for DESI was provided by the U.S. National Science Foundation (NSF), Division of Astronomical Sciences under Contract No. AST-0950945 to the NSF’s National Optical-Infrared Astronomy Research Laboratory; the Science and Technology Facilities Council of the United Kingdom; the Gordon and Betty Moore Foundation; the Heising-Simons Foundation; the French Alternative Energies and Atomic Energy Commission (CEA); the National Council of Humanities, Science and Technology of Mexico (CONAHCYT); the Ministry of Science and Innovation of Spain (MICINN), and by the DESI Member Institutions (https://www.desi.lbl.gov/collaborating-institutions). Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the U. S. National Science Foundation, the U. S. Department of Energy, or any of the listed funding agencies. The authors are honored to be permitted to conduct scientific research on Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation.

References

Appendix A: Cosmic variance

Throughout this work, we trained and tested FORESTFLOW using simulations run employing the "fixed-and-paired" technique (Angulo & Pontzen 2016; Pontzen et al. 2016), which significantly reduces cosmic variance for the clustering of the Lyman-α forest (Anderson et al. 2019). We could further mitigate the impact of cosmic variance by using control variates (Kokron et al. 2022), but this is outside the scope of the current work. The impact of cosmic variance on fixed-and-paired simulations is not straightforward (Maion et al. 2022), and thus we would ideally use multiple fixed-and-paired simulations with different initial distributions of Fourier phases to estimate the precision of measurements from our simulations. However, we only have two simulations with these properties: CENTRAL and SEED. In this section, we use these two simulations to estimate the impact of cosmic variance on simulation measurements and best-fitting models. It is crucial to acknowledge that our findings are subject to significant noise because we only have access to two independent realizations.

In Fig. A.1, we show the difference between measurements from the CENTRAL and SEED simulations at z = 3, normalized by $\sqrt{2}$ times their average¹². The CENTRAL and SEED simulations differ only in their initial Fourier phase distributions, allowing their difference to isolate the effects of cosmic variance. Unlike traditional simulations, where cosmic variance for P_3D scales inversely with the square root of the number of modes, this uncertainty peaks at ≃10% around k ≃ 0.3 Mpc⁻¹ and decreases at both larger and smaller scales. This behavior can be explained as follows: at the largest scales, the fix-and-paired technique cancels cosmic variance for linear density modes, reducing variance. At intermediate scales, nonlinear evolution, particularly mode coupling, reintroduces cosmic variance, increasing the uncertainty. At smaller scales, the increasing number of modes leads to a decrease in cosmic variance, similar to trends observed in traditional simulations.

Fig. A.1. Impact of cosmic variance on P3D (top panel) and P1D (bottom panel) measurements from our simulations at z = 3. Lines show the difference between measurements from the CENTRAL and SEED simulations, which only differ on their initial distribution of Fourier phases, divided by $\sqrt{2}$ times their average. Cosmic variance induces errors as large as 10% on P3D for k ≃ 0.3 Mpc−1, while these are on the order of 1% for P1D.

Thus, cosmic variance limits our ability to accurately evaluate both the parametric model and FORESTFLOW using our simulations. To mitigate its impact, we restricted the analysis of P_3D performance to scales k > 0.5 Mpc⁻¹ in the main results. In contrast, cosmic variance has a much smaller effect on P_1D, contributing only about 1.5% uncertainty at k_∥ < 2 Mpc⁻¹. This smaller effect allowed us to include all scales in P_1D tests without concern. To more precisely quantify the impact of cosmic variance on simulation measurements, we computed the standard deviation of the results shown in Fig. A.1 across redshift. We did it within the intervals 0.5 < k[ Mpc⁻¹]< 5 for P_3D and 0.09 < k_∥[ Mpc⁻¹]< 4 for P_1D, based on the scales discussed earlier and those used for fitting the P_3D model in Sect. 3.3. We found that the average impact of cosmic variance is 1.3% for P_3D and 0.5% for P_1D.

We expect the impact of cosmic variance on the best-fitting model to P_3D and P_1D measurements to be weaker than on simulation measurements, as multiple P_3D and P_1D bins collectively constrain the eight parameters of the model. In Fig. A.2, we show the difference between the best-fitting models to the CENTRAL and SEED simulations, normalized by $\sqrt{2}$ times the best-fitting model to their average. The top panel displays results for the two Lyman-α linear biases, b_δ and b_η. The standard deviation of the differences is 0.6% and 1.8% for b_δ and b_η, respectively, demonstrating that we can measure the two Lyman-α linear biases with percent-level accuracy using our simulations. This precision is achievable because the combination of small and large scales in the fits helps break degeneracies between the two linear biases and the other six model parameters.

Fig. A.2. Impact of cosmic variance on predictions from the parametric model. Lines show the difference between the best-fitting models to P3D and P1D measurements from the CENTRAL and SEED simulations, divided by $\sqrt{2}$ times the best-fitting model to their average. The top panel shows the results for the Lyman-α linear biases (bδ and bη), while the middle and bottom panels display the results for P3D and P1D at z = 3, respectively. The impact of cosmic variance on model predictions is approximately an order of magnitude smaller than on simulation measurements (see Fig. A.1).

By propagating these uncertainties to the behavior of P_3D on linear scales, we find that the impact of cosmic variance on perpendicular and parallel modes is 1.2% and 1.8%, respectively. In the middle and bottom panels of Fig. A.2, we analyze the influence of cosmic variance on model predictions for P_3D and P_1D. The overall impact of cosmic variance on P_3D and P_1D predictions is 0.8% and 0.1%, respectively, confirming that the best-fitting model is less sensitive to cosmic variance than individual simulation measurements. Consequently, FORESTFLOW is more robust against cosmic variance than emulators that predict the power spectrum at fixed k-bins.

Appendix B: Validation of uncertainties predicted for P_3D and P_1D

Normalizing flows predict the full posterior distribution of the target data rather than only their mean like fully connected neural networks or their mean and width like Mixture Density Networks (see Ramachandra et al. 2022; Cabayol-Garcia et al. 2023, for some applications in cosmology). This is achieved through multiple sampling iterations from the target latent distribution, an eight-dimensional Gaussian in our case. In FORESTFLOW, each sampled realization of the P_3D model parameters is propagated to generate predictions for P_3D and P_1D (see Sect. 4.1), producing a covariance matrix for these statistics. In this appendix, we validate its diagonal elements. Note that well-calibrated uncertainties are critical for future uses of FORESTFLOW such as cosmology inference.

We validated the uncertainty in P_3D and P_1D predictions using the Probability Integral Transform test (PIT), which is the value of the cumulative distribution function (CDF) of a distribution evaluated at the ground-truth value z_t

$$\begin{array}{c}\hfill \mathrm{PIT}=\mathrm{CDF}[\mathrm{p},z_{\mathrm{t}}]={\displaystyle {\int }_{-\infty }^{z_{\mathrm{t}}}}p(z)\mathrm{d}\mathrm{z},\end{array}$$(B.1)

where p is in our case the distribution of FORESTFLOW predictions for P_3D or P_1D and z_t stands for measurements of these statistics from the simulations. A model that displays a well-calibrated uncertainty distribution yields PIT values that are uniformly distributed between zero and one. This indicates that the observed outcomes have an equal likelihood of falling at any point along the predicted CDF. In contrast, an excess of values close to zero or one indicates that the width of the distribution is underestimated.

In Fig. B.1, we display a PIT test produced using all the TRAINING simulations via a leave-simulation-out approach (see Sect. 5.1). This process validates average predictions and uncertainties against simulations excluded in the training process. The red and blue lines display the results for P_3D and P_1D, respectively, which were generated by combining results from different scales and redshifts. As we can see, the PIT distribution is approximately uniform for the two statistics but it presents peaks at the low and high ends, indicating underestimated uncertainties for some samples. The cause behind this feature is unclear and it demands further investigation beyond the scope of this project.

Fig. B.1. PIT distribution for P1D (blue) and P3D (red). This plot validates the uncertainties predicted by FORESTFLOW across TRAINING simulations via a leave-simulation-out approach. The PIT distribution is approximately uniform, indicating well-calibrated uncertainties for most samples, while the peaks at the edges indicate underestimated uncertainties for some samples.

All Tables

All Figures

Fig. 1. Accuracy of the P3D model (see Eqs. (3) and (4)) in reproducing measurements from the CENTRAL simulation at z = 3. In the top panel, dotted and solid lines show the ratio of simulation measurements and model predictions relative to the linear power spectrum, respectively. Dashed lines do so for the linear part of the best-fitting model (DNL = 1). Line colors correspond to different μ wedges, and vertical dashed lines mark the minimum scale used for computing the best-fitting model, k = 5 Mpc−1. The bottom panel displays the relative difference between the best-fitting model and simulation measurements. The overall accuracy of the model is 2% on scales in which simulation measurements are not strongly affected by cosmic variance (k > 0.5 Mpc−1; see text).

In the text

Fig. 2. Accuracy of the parametric model (see Eqs. (3) and (4)) in reproducing P3D and P1D measurements from all the TRAINING simulations. Lines and shaded areas show the mean and standard deviation of the relative difference between simulation measurements from the 1650 snapshots of the TRAINING simulations and best-fitting models to these, respectively. The accuracy of the model in recovering P3D and P1D is 2.4 and 0.6%, respectively, on scales not strongly affected by cosmic variance.

In the text

Fig. 3. Architecture of the Lyman-α forest clustering emulator. The blue arrow indicates the training direction, where the cNF optimizes a bijective mapping between the best-fitting parameters of the P3D model to measurements from the TRAINING simulations and an eight-dimensional Normal distribution. The mapping is conditioned on cosmology and IGM physics, and performed using 12 consecutive affine coupling blocks. The green arrow denotes the emulation direction, where the cNF applies the inverse of the mapping to random samples from the base distribution to predict the value of the P3D model parameters. Outside the cNF, FORESTFLOW introduces these parameters in Eq. (3) and (1) to obtain predictions for P3D and P1D, respectively.

In the text

Fig. 4. Accuracy of FORESTFLOW in recovering P3D and P1D measurements from the CENTRAL simulation at z = 3. Dotted lines show measurements from simulations, solid lines and shaded areas display the average and 68% credible interval of FORESTFLOW predictions, respectively, and vertical dashed lines indicate the minimum scales considered for computing the training data for the cNF. The overall performance of FORESTFLOW in recovering P3D is 2.0% on scales not strongly affected by cosmic variance and 0.6% for P1D.

In the text

Fig. 5. Accuracy of FORESTFLOW in recovering P3D and P1D measurements from the CENTRAL simulation as a function of redshift. The upper four panels show the results for P3D across different μ bins, while the bottom panel displays the results for P1D. Each color represents a different redshift. The model’s accuracy remains consistent for redshifts above z = 2 but exhibits a slight decline at this redshift.

In the text

Fig. 6. Accuracy of FORESTFLOW across the input parameter space estimated via leave-simulation-out (top panels) and leave-redshift-out tests (bottom panels). Top panels. Each leave-simulation-out test involves training one independent emulator with measurements from 29 distinct simulations, and then using the measurements from the remaining simulation as the validation set. Lines and shaded areas show the average and standard deviation of 30 leave-simulation-out tests, and each panel shows the results for a different redshift. Bottom panels. Leave-redshift-out tests require optimizing one emulator with all measurements but the ones at a particular redshift, and then using measurements from this redshift as validation. Each panel shows the results of a different test.

In the text

Fig. 7. Performance of FORESTFLOW in recovering P3D and P1D for test simulations not included in the training set. Lines and shaded areas display the average and standard deviation of the results for 11 snapshots between z = 2 and 4.5, respectively. From top to bottom, the rows show the results for the CENTRAL, SEED, GROWTH, NEUTRINOS, CURVED, and REIONISATION simulations, where the CENTRAL and SEED simulations are at the center of the input parameter space and employ the same and different initial distribution of Fourier phases as the training simulations, respectively, the GROWTH and REIONISATION simulations use a different growth and reionization history relative to those used by the TRAINING simulations, and the NEUTRINOS and CURVED simulations consider massive neutrinos and curvature. The efficiency of FORESTFLOW is approximately the same for all simulations.

In the text

Fig. 8. Response of Lyman-α clustering to variations in cosmology and IGM physics according to FORESTFLOW. The top, middle, and bottom panels of the left column show the results for β, bδ, and bη, respectively, while those of the right column do so for the perpendicular modes of P3D, the parallel modes of P3D, and P1D. Blue, orange, and red lines show the response of the previous quantities to a 5% increase in As, ΩMh2, and σT, respectively, while green lines do so for a 1% increase in $\overline{F}$.

In the text

Fig. 9. Accuracy of the parametric model and FORESTFLOW in describing P× measurements from the CENTRAL simulation at z = 3. Dots show simulation measurements, dashed lines depict predictions from the best-fitting parametric model to P3D and P1D measurements, and solid lines and shaded areas display the average and 68% credible interval of FORESTFLOW predictions. The color of the lines indicates the results for different bins in transverse separation r⊥. The middle panel shows the residual between simulation measurements and the best-fitting parametric model, while the bottom panel displays the residual between predictions from the parametric model and FORESTFLOW. The performance of FORESTFLOW in reproducing simulation measurements is similar to that of the best-fitting model.

In the text

Fig. A.1. Impact of cosmic variance on P3D (top panel) and P1D (bottom panel) measurements from our simulations at z = 3. Lines show the difference between measurements from the CENTRAL and SEED simulations, which only differ on their initial distribution of Fourier phases, divided by $\sqrt{2}$ times their average. Cosmic variance induces errors as large as 10% on P3D for k ≃ 0.3 Mpc−1, while these are on the order of 1% for P1D.

In the text

Fig. A.2. Impact of cosmic variance on predictions from the parametric model. Lines show the difference between the best-fitting models to P3D and P1D measurements from the CENTRAL and SEED simulations, divided by $\sqrt{2}$ times the best-fitting model to their average. The top panel shows the results for the Lyman-α linear biases (bδ and bη), while the middle and bottom panels display the results for P3D and P1D at z = 3, respectively. The impact of cosmic variance on model predictions is approximately an order of magnitude smaller than on simulation measurements (see Fig. A.1).

In the text

	Fig. B.1. PIT distribution for P1D (blue) and P3D (red). This plot validates the uncertainties predicted by FORESTFLOW across TRAINING simulations via a leave-simulation-out approach. The PIT distribution is approximately uniform, indicating well-calibrated uncertainties for most samples, while the peaks at the edges indicate underestimated uncertainties for some samples.
In the text