S Peirani

S Prunet

S Colombi

C Pichon

D H Weinberg

C Laigle

G Lavaux

Y Dubois

J Devriendt

LyMAS reloaded: improving the predictions of the large-scale Lyman-α forest statistics from dark matter density and velocity fields

Keywords:

published or not. The documents may come

INTRODUCTION

Distant quasars emit light that crosses a large part of the Universe before being observed with instruments on Earth. In particular, the spectrum of each quasar presents fluctuating absorption that corresponds to the ⋆ E-mail: sebastien.peirani@oca.eu Lyman-α forest (Ly-α, Lynds (1971); Sargent et al. (1980)). The study of the Lyman-α forest has become a major focus of modern cosmology, as it is supposed to trace the neutral hydrogen density that fills most of the Universe in a way that approximately corresponds to the underlying dark matter density (Croft et al. 1999;Peeples et al. 2010). Since a single background source only provides one dimen-sional information along the corresponding line-ofsight (or "skewer"), characterizing the 3-d density of the high redshift Universe with the Ly-α forest requires large samples of quasar spectra. Successful surveys such as the (extended) Baryon Oscillation Spectroscopic Survey (BOSS/eBOSS, Dawson et al. (2013Dawson et al. ( , 2016))), of the Sloan Digital Sky Survey (SDSS-III and SDSS-IV, Blanton et al. (2017); Eisenstein et al. (2011) have measured the Lyman-α forest spectra of 160,000 quasars at redshifts 2.2 < z < 3. Thanks to this large sample, the study of the Lyman-α forest has proved to be a complementary probe to low redshift galaxy surveys. For instance, the large sample of quasar spectra have permitted accurate measurements of 3-d flux auto-correlation functions (Slosar et al. 2011) as well as the cross-correlation between the Ly-α Forest and specific tracers, namely damped-Lyα systems (DLAs) and quasars (Font-Ribera et al. 2012, 2013, 2014). Such 3-d measurements also enable measurements of the distance-redshift relation and the Hubble expansion via baryon accoustic oscillations (Busca et al. 2013;Slosar et al. 2013;Delubac et al. 2015;Bautista et al. 2017;du Mas des Bourboux et al. 2020). Moreover, BOSS spectra also permit accurate measurements of the line-of-sight power spectrum (Palanque-Delabrouille et al. 2013) and flux probability distribution function (PDF) (Lee et al. 2015). In the near future, the Dark Energy Spectroscopic Instrument (DESI, DESI Collaboration et al. ( 2016)), the William Herschel Telescope Enhanced Area Velocity Explorer (WEAVE-QSO, [START_REF] Pieri | Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics[END_REF]; [START_REF] Dalton | Society of Photo-Optical Instrumentation Engineers[END_REF][START_REF] Dalton | Society of Photo-Optical Instrumentation Engineers[END_REF]) and the Subaru Prime Focus Spectrograph (PFS, Takada et al. (2014)) will go well beyond the present surveys and will open new perspectives on the high redshift intergalactic medium probed by the Ly-α forest.

In parallel to the development of these large quasar surveys, theoretical modeling needs to reach the high level of complexity and accuracy to interpret the observational data. Nowadays, hydrodynamical cosmological simulations represent an ideal tool as they manage to model the intergalactic medium with a high degree of realism with appropriate resolution (e.g. Dubois et al. 2014;Vogelsberger et al. 2014;Schaye et al. 2015;Bolton et al. 2017). However, to properly model the 3-d correlations of the Ly-α forest, one needs to resolve the pressure-support scale (Jeans scale) of the diffuse intergalactic medium (IGM), which is typically ∼0.25 Mpc/h comoving for matter overdensity of ∼10 (Peeples et al. 2010), while considering ∼Gpc 3 simulation volumes to exploit the statistical precision achieved by the different observational surveys while avoiding box size effects. Combining such resolution and simulation volume is currently not feasible mainly because of computational limits. To tackle such an issue, several methods exist in the literature. One of the most popular is to use the so-called "Fluctuating Gunn-Peterson Approximation" (FGPA [START_REF] Weinberg | Structure and Evolution of the Intergalactic Medium from QSO Absorption Line System[END_REF]Croft et al. 1998) which links the Ly-α optical depth to the local dark matter density. This approach is relatively straightforward as it assumes a deterministic relation and only information on the density field (extracted from N-body simulations or log-normal density fields) is required. However, the FGPA approach is expected to be accurate enough only on very large scales, e.g. those of the BAO features (∼ 100 Mpc/h, e.g. Fig. 5 of Sinigaglia et al. (2021)). But 3d Ly-α forest surveys also enable precise measurements of flux correlations at much smaller scales where the FGPA might not be adequate.

Another approach is to apply relevant calibrations, derived first from small hydrodynamical simulations, to large-scale dark matter distributions extracted from pure dark matter simulations, which are much cheaper to perform. In particular, Peirani et al. (2014) (hereafter, P14) have developed the Lyman-α Mass Association Scheme (LyMAS) which follows such a philosophy. The main idea is that flux correlations on small and large scales are mainly driven by the correlations of the dark matter density field. More specifically, the flux statistics can be estimated by combining the DM density field with the conditional probability distribution P (F |ρ) of the transmitted flux F on the DM density contrast ρ. In its most sophisticated form, LyMAS creates ensemble of coherent pseudo-spectra at the BOSS resolution using the Gaussianized percentile distribution of the conditional flux, while rescaling the line-of-sight power spectrum and PDF at the last step. One of the main results of LyMAS is to improve the predictions of flux 3d correlations especially with respect to deterministic mapping (e.g. FGPA) which tends to significantly overestimate them especially when the DM density is smoothed at scale greater than 0.3 Mpc/h. Similarly, Sorini et al. (2016) have developed "Iterative Matched Statistics" (IMS) in which the PDF and the power spectrum of the realspace Ly-α flux are derived from small hydrodynamical simulations. Then, these two statistics are 1d (1D-IMS) or 3d (3D-IMS) iteratively mapped onto a pseudo-flux field of an N-body simulation from which the matter density is first Gaussian smoothed. In 3D-IMS, smoothing is followed by matching the 3D power spectrum and PDF of the flux in real space to the reference hydrodynamic simulation. With 1D-IMS, the 1d power spectrum and PDF of the flux are additionally matched. Both methods have proved to be again more accurate than the FGPA approach (which strongly relies on the DM smoothing scale) when reproducing line-of-sight observables, such as the PDF and power spectrum as well as the 3-d flux power spectrum (5-20%). Finally, Machine Learning based methods start to be considered and lead to promising results (Harrington et al. 2021;Sinigaglia et al. 2021;Chopitan et al. 2021).

Although the LyMAS full scheme is able to model the BOSS 3-d clustering quite accurately and has been already used in different analysis related to the quasar-Ly-α forest cross-correlation (Lochhaas et al. 2016), the three-point correlation functions (Tie et al. 2019) and the correlations between the Ly-α transmitted flux and the mass overdensity (Cai et al. 2016), we aim at investigating whether other sets of calibrations could still improve the theoretical predictions. To this regard, we consider a new approach based on Wiener Filtering, which has been used for 3-d map reconstruction from an ensemble of 1d lines-of-sight (e.g. Pichon et al. 2001;Caucci et al. 2008;Ozbek et al. 2016;Lee et al. 2018;Japelj et al. 2019;[START_REF] Ravoux | [END_REF]. The underlying philosophy in LyMAS2 remains unchanged. We still find that the flux correlations are driven mainly by the correlations of the DM density field, but with potential refinements from the correlations of the DM velocity field.

The paper is organised as follows. In section 2, we describe the Wiener equations as multivariate normal conditional probabilities, and we explain their application to hydrodynamical simulations. Section 3 briefly describes how we extract the flux and all relevant dark matter fields from the Horizon-noAGN simulation. We also present the potential correlations that arise between these different fields. Then in section 4 we present the statistics in the line of sight power spectrum, the PDF and the 2-point correlation function of pseudo-spectra produced when LyMAS2 is applied to different associations of dark matter fields of Horizon-noAGN. Such trends are compared to the flux statistics from the hydrodynamical simulation ("hydro flux"). In section 5 we apply LyMAS2 to N-body simulations of 1.0 and 1.5 Gpc/h comoving volumes. We summarize our results and conclusions in section 6. We also add three appendices. In appendix A, we compute the mean two point correlations functions from five different hydrodynamical simulations of lower resolution to check the robustness of the results presented in section 4. Appendix B presents the performance of specific deterministic samplings. Appendix C provides details on how estimates of density and velocities are performed on the dark matter particle distribution, relying on adaptive softening.

WIENER EQUATIONS

Multivariate conditional probabilities

Let us assume a complex Gaussian random (vector) variable x that can be separated into two sub-vectors x = (x1, x2), whose mean and covariance can be written as

µ = (µ1, µ2), Σ = Σ11 Σ12 Σ H 12 Σ22 ,
where the H superscript denotes hermitian conjugate. By using formulae for block inverses, it is possible to derive from the joint distribution of x1 and x2 the formula for the conditional distribution of x1 given x2. As expected, it is a gaussian multivariate distribution, of mean and covariance:

μ1 = µ1 + Σ12Σ -1 22 (x2 -µ2) , (1)
Σ = Σ11 -Σ12Σ -1 22 Σ H 12 . (2)
Now, consider the joint gaussian distribution of the (complex) spatial Fourier modes (f k , δ k , θ k , ω k ), where, by definition for field a, we have a k = a(x)e -ik•x d 3 x. We assume they are centered fields (zero mean), and of covariance (we drop the subscript k for lisibility):

Σ =     P ff P f δ P f θ P f ω P * f δ P δδ P δθ P δω P * f θ P * δθ P θθ P θω P * f ω P * δω P * θω Pωω     = Σ11 Σ12 Σ H 12 Σ22
,

where by definition P ab (k) ≡ a * k b k is the crossspectrum of fields a and b at wavenumber k, and we have partitioned the fields according to x1 ≡ f k and x2 ≡ (δ k , θ k , ω k ) T . Applying the equations ( 1) and ( 2) we get the conditional mean and variance of the field

f k given (δ k , θ k , ω k ): fk = P f δ , P f θ , P f ω • Σ -1 22 •   δ k θ k ω k   , (3)
Σ = P ff -P f δ , P f θ , P f ω • Σ -1 22 •   P * f δ P * f θ P * f ω   . (4)
Computing the inverse Σ -1 22 using the cofactor matrix formula, we get where "c.c" denotes the conjugate complex.

Σ -1 22 = 1 |Σ22|   A11 A12 A13 A21 A22 A23 A31 A32 A33
Note that we limit our study to a maximum of three different input fields (i.e. δ, θ and ω) to construct the field f . But obviously, this formalism can be extended to a higher number of fields, leading to more and more complex analytical solution. However, we will see that the statistical trends derived when considering two and three input fields are quite similar (when judiciously chosen), suggesting that adding more than two fields will not noticeably improve the results anymore.

Application to hydro simulations

Let us call F , ρ, v1 and v2 respectively the local Lyα transmitted flux, the local dark matter density, velocity divergence and vorticity amplitude, extracted from a given hydrodynamical simulation. Let us call now G the cumulative distribution of a standard normal N (0; 1)

distribution (G[x] =

x exp(-u 2 /2)/ √ 2πdu), and

Gρ = ρ -∞ PDF(ρ ′ )dρ ′ , (5)
Gv i = v i -∞ PDF(v ′ )dv ′ , (6)
GF = F 0 PDF(F ′ )dF ′ , (7)
the cumulative distributions of the measured DM density and velocity fields (in the hydrodynamical simulation), and G f the cumulative distribution of the measured flux, F . We can then define new fields, namely

f = G -1 (GF (F ), δ = G -1 (Gρ(ρ)), (8)
θ = G -1 (Gv 1 (v1)), ω = G -1 (Gv 2 (v2)) . ( 9
)
which should be normally distributed by construction (or "gaussianized"). Let us compute these different fields from the hydro simulations and extract from them the relevant cross-spectra, using Fourier space:

P ff = f * k f k , P δδ = δ * k δ k , P θθ = θ * k θ k , Pωω = ω * k ω k , P f δ = f * k δ k , P f θ = f * k θ k , P δθ = δ * k θ k , P θω = θ * k ω k , etc...
which are going to depend typically on a transverse and a longitudinal separation radius. If we assume that the fields f , δ, θ and ω are Gaussian random fields (GRF, not just its one point statistics is now required to be normal) then for a given measured set of ρ, v1 and v2 (correspondingly a set of δ k , θ k and ω k ) say along a set of lines of sight (LOS), one can estimate the most likely field f following equations (3):

fk = T1 • δ k + T2 • θ k + T3 • ω k , (10)
where T1, T2 and T3 are functions of cross-spectrum P ab (k). This approach can be done along a given LOS or in volume. However, in the present study we will only analyse LOS individually and independently, ignoring transverse correlations between different LOS for now. This allows us to use 1D Fast Fourier Transforms (FFT), and assuming stationarity along the LOS, the multiplication is simply done frequency by frequency since in Fourier space the covariance sub-blocks P ab are diagonal. To illustrate, if only one field is considered (i.e. the DM overdensity field in our study) we simply get:

fk = P f δ • P -1 δδ • δ k . ( 11
)
If we add additional information from a specific velocity field (e.g. velocity dispersion), the expression of f becomes:

fk = P f δ (P θθ -P δθ ) P δδ P θθ -P δθ P * δθ •δ k + P f θ (P δδ -P * δθ ) P δδ P θθ -P δθ P * δθ •θ k . (12)
Adding a second velocity field leads to a more complex expression of fk .

We can then draw samples as fk ≡ fk + ∆f k where ∆f k obeys a GRF of mean zero and variance as equation (4):

∆f k ∼ G(0, P ff -P f δ • P -1 δδ • P * f δ ) ,
when only the dark matter density is considered and

∆f k ∼ G(0, P ff - P f δ P f θ T • P δδ P δθ P * δθ P θθ -1 • P * f δ P * f θ ),
in the case of two dark matter fields. After computing the inverse Fourier transform of fk to get f , the corre-

sponding flux obeys F = G f [G -1 ( f )]
. By construction the one point statistics of f will be random normal, so that the one point PDF of its de-gaussianized transform will be that of the original field. The power spectrum P f f of f will be the same as that of f . Indeed, let us consider the case with only one input file for simplicity.

We have:

P fk fk ≡ | fk + ∆f k | 2 = 1 P δδ P f δ δ k 2 + |∆f k | 2 ,
because the expectation and the fluctuations are uncorrelated. Therefore

P fk fk = P 2 f δ P 2 δδ |δ k | 2 + P ff - P 2 f δ P δδ = P ff , since |δ k | 2 ≡ P δδ .
Recall that all equations above are valid independently for each Fourier mode k, and for each mode, all P ab terms are scalars for a given pair of fields (a, b).

FLUX AND DARK MATTER FIELDS

Throughout the present analysis, we have used the Horizon-noAGN simulation (Peirani et al. 2017) to characterize any relevant cross correlations between the transmitted flux and the different dark matter fields. Horizon-noAGN is a hydrodynamical simulation of 100 h -1 Mpc comoving box side run with the RAM-SES code (Teyssier 2002). It evolves 1024 3 dark matter particles with a mass resolution of 8.27×10 7 M⊙ while the initially uniform grid is refined in an adaptive way down to ∆x = 1 proper kpc at all times. The simulation adopts a standard ΛCDM cosmology compatible with WMAP-7 results (Komatsu et al. 2011), namely a total matter density Ωm=0.272, a dark energy density ΩΛ=0.728, an amplitude of the matter power spectrum σ8=0.81, a baryon density Ω b =0.045, a Hubble constant H0 = 70.4 km s -1 Mpc -1 and ns=0.967. Horizon-noAGN is the twin simulation of Horizon-AGN (Dubois et al. 2014). It contains all relevant physical processes such as metal-dependent cooling, photoionization and heating from a UV background, supernova feedback and metal enrichment, but does not include black hole growth and therefore AGN feedback.

The choice of using Horizon-noAGN instead of Horizon-AGN was mainly motivated by the fact that we have performed five additional but slightly lower resolution hydrodynamical simulations to estimate the accuracy and robustness of the LyMAS2 predictions. Thus, turning out the AGN feedback processes in the simulations has permitted us to limit the computational time. These results are presented in appendix A. Furthermore, we have tuned in the present study the UV background intensity in the process of generating the "noAGN" flux grid (see below) in order to get the same mean transmitted Ly-α forest flux F derived from Horizon-AGN. By doing this, the flux statistical predictions from the two simulations in the 3-d Ly-α clustering tend to be almost the same. This has been already noticed in Lochhaas et al. (2016) when studying the crosscorrelations between dark matter haloes and transmitted flux in the Ly-α forest. Note however that AGN feedback is expected to have non negligible effect on the Ly-α 3d clustering such as, for instance, the 1-d power spectrum of the Ly-α forest (e.g. Viel et al. 2013;Chabanier et al. 2020).

In the following, we describe briefly how we derived the hydro spectra field and the different DM density and velocity fields from the Horizon-noAGN simulation. Similarly to P14, we analyse the simulation outputs at redshift 2.5. Each field is sampled in a regular grid of 1024 3 pixels and the size of a single pixel is therefore ∼0.1 comoving Mpc/h or 0.04 physical Mpc.

Transmitted flux

From the Horizon-noAGN, we follow the method to generate the hydro spectra that is fully described in (Peirani et al. 2014). The optical depth of Ly-α absorption is calculated based on the neutral hydrogen density along each line of sight. Basically, the opacity at observer-frame frequency ν obs is τ (ν obs ) = cells nHIσ(ν obs )dl where the sum extends over all cells traversed by the line-of-sight, nHI is the numerical density of neutral H atoms in each cell, σ(ν obs ) is the cross section of Hydrogen to Ly-α photons, and dl is the physical size of the cell. Then each spectrum is smoothed with a 1-d Gaussian of dispersion 0.696 h -1 Mpc, equivalent to BOSS spectral resolution at z ≈ 2.5. The optical depth along each spectrum is converted to Ly-α forest flux by F = e -τ . Following common practice in Lyα forest modelling, the UV background intensity is chosen to give a mean transmitted Ly-α forest flux F = e -τ =0.795, matching the metal-corrected z = 2.5 value measured from high-resolution spectra by Faucher-Giguère et al. (2008).

3.2 Density, velocity and mean square velocity DM skewers that correspond to the "hydro" spectra are also extracted from the hydrodynamical simulation. We use the same 3-step algorithm introduced in P14 to derive both the overdensity, the velocity field and the velocity dispersion fields:

(i) adaptive interpolation of the DM particle distribution on a regular grid (Colombi et al. 2007), as detailed in Appendix C;

(ii) smoothing with a Gaussian window in Fourier space;

(iii) Extraction of the skewers from a grid of LOS aligned along the z-axis.

In step (ii), DM field is three-dimensionally smoothed using different choices of smoothing scales. In P14, we found that a smoothing scale of 0.3 Mpc/h has proved to be optimal leading to the most accurate predictions. However, we prefer a value of 0.5 Mpc/h in the present study since the predictions are very similar to those obtained with 0.3 Mpc/h (see Appendix A). Furthermore, we anticipate with the fact that it is computationally much easier to smooth a DM field to 0.5 Mpc/h than 0.3 Mpc/h when considering large volumes namely with box sides at least greater or equal to 1 Gpc/h.

Dark matter vorticity

The velocity field is projected (using Cloud-in-Cell interpolations) on a regular grid of resolution 1024 and smoothed over 0.5 Mpc/h with a gaussian filter. The vorticity Ω is then computed as being the curl of the velocity field using FFT. Slightly smoothing the input velocity field allows us to avoid Gibbs artefacts.

Dark matter velocity divergence

We have considered both the 3d velocity divergence and the 1d velocity divergence along the line of sight direction.

For the 3d case, we employed two different methods to see whether this could affect our results and trends. The first one is based on a centered finite difference approximation, namely the divergence ∇ i 3d at a pixel i is given by

∇ i 3d ≃ V i+1 x -V i-1 x 2h + V i+1 y -V i-1 y 2h + V i+1 z -V i-1 z 2h , (13)
where V i x , V i y and V i z are the velocity components at pixel i and h the size of a pixel (i.e. 100/1024 Mpc/h here). The second method uses the exact expression of the divergence in Fourier space. However, we found that the two methods lead to very similar results so we will only show results from the Fourier space method for the 3d case.

For the 1d case, we simply use the finite difference approach and the divergence ∇ i 1d at a pixel i becomes:

∇ i 1d ≃ V i+1 z -V i-1 z 2h , (14)
since we define the z-axis as the direction of the LOS.

We summarize in Table 1 the different DM fields used in this study. It is worth mentioning that we changed some of the notations that can be found in P14. We first replaced the definition of the smoothed flux at the BOSS resolution Fs in P14 to F here. We also changed the definition of the 3d smoothed DM overdensity, ρ = 1 + δ instead of ∆s = 1 + δs in P14.

In Fig. 1, we show a slice of the hydro flux (smoothed at BOSS resolution) and the corresponding DM density and velocity fields (smoothed at 0.5 Mpc/h) extracted from the Horizon-noAGN simulation at redshift 2.5. As expected, clear correlations are noticeable between the transmitted flux and the different DM fields. This trend can also be seen when studying the evolution of each field along the same LOS, and a typical example is given in Fig. 2. We note that high absorptions in the flux correspond to high density regions or high values in the vorticity or the velocity dispersion. But the relative amplitudes of peaks in the density contrast may differ from those of the the velocity dispersion/vorticity. Indeed, the density contrast and the velocicy dispersion/vorticity do not necessary put the emphasis of the same structures (e.g. walls, filaments) as suggested by Fig. 1 or, for instance, by Fig. 2 of Buehlmann & Hahn (2019). On the contrary, these high absorptions rather coincide with high negative values in the 3-d or 1-d velocity divergence. This is consistent since high density regions are associated with DM haloes in which matter tends to sink toward the center of the objects. Note also that the variations of the modulus of the vorticity and velocity dispersion are very similar.

Field cross correlations

In order to characterize the correlations that emerge from Figures 1 and2, we first plot in Fig. 3 some relevant scatter plots between the optical depth τ and the DM overdensity and velocity fields. The correlations between the optical depth in the hydro spectra and the smoothed DM overdensity (1+δ) is quite similar to the trend found in P14 using the 50 h -1 Mpc "Horizon MareNostrum" simulation. In Fig. 4, we additionally show the correlations between the different DM fields. As noticed in Fig. 2, the velocity dispersion and vorticity field are highly correlated. We do not show the correlations using the 1d velocity divergence as there are quite similar with trends found using the 3d velocity divergence.

All these plots suggest that there are more or less pronounced correlations between the different input DM fields. It is however tricky to anticipate which combinations of fields through the LyMAS2 scheme would lead to the most accurate theoretical predictions. As specified in section 2, we consider combinations with up to three different DM fields which offers 85 different possibilities (5, 20 and 60 respectively for one, two and three input fields). However, as the main philosophy of Ly-MAS is to trace the Ly-α flux from the underlying DM distribution with potential corrections from the DM veloctiy field, we will always consider the DM overdensity field in each combination reducing this number to 17. Moreover, since the velocity dispersion and the vorticity fields are highly correlated, we will also always use the velocity dispersion in the 3d case. Consequently, we limit our study to 8 different combinations presented in Table 2. Nevertheless, we have checked that combinations using only DM velocity fields do not lead to satisfactory theoretical predictions.

From each specific association of DM fields, and each Fourier mode k along a LOS, we have estimated the relevant cross-spectra P ab defined in section 2.1, where a and b refer either to the transmitted flux, the DM overdensity or a specific DM scalar field derived from the DM velocity field. For each mode k, the covariance matrix P ab is hermitian, and its linear dimension is equal to the total number of fields considered. Examples of cross-power spectra are shown in Fig. 5. We also derived the relevant 1D power spectrum P k required to computed the covariance ∆f k defined in equation ( 4). An example of P k is shown in Fig. 11.

CREATING PSEUDO-SPECTRA

In this section, we apply the LyMAS2 scheme to the DM fields extracted from the Horizon-noAGN simulation to generate grids of pseudo-spectra at BOSS resolution.

The objective is to recover the 3d Ly-α clustering statistics of the "true" hydro spectra. For a given skewer, we summarize the main steps to follow to produce a corresponding pseudo-spectrum: 1. We gaussianize the overdensity field (ρ) and eventually one or two additional DM velocity fields (e.g. σ)

2. We compute the FFT of each gaussianized field. This gives new (complex) fields, ρ k , σ k , etc.. 3. We compute (in Fourier space) the most probable flux fk = T1 • ρ k + T2 • σ k + ..., by applying the relevant filters T1, T2,... (see for instance equation ( 12) for the 2-fields case).

4. We generate a 1d Gaussian field of mean 0 and variance defined in equation ( 4) to get the covariance ∆f k . 5. After computing the inverse Fourier transform of fk + ∆f k to get f , we degaussianize to get the pseudospectrum: F = G -1 (f ).

6. One full iteration. We follow the same strategy described in P14. We first measure the 1-d flux power spectrum Pps(k) of the pseudo-spectra created in this way. Then we Fourier transform each pseudo-spectrum and multiply each of its Fourier component by the ratio [PF (k)/Pps(k)] 1/2 , inverse transform to get the same 1-d flux power spectrum than of the true hydro spectra PF (k) (Weinberg & Cole 1992). The second step of the full iteration is to compute the PDF of the pseudospectra after the 1-d Pk re-scaling and then monotonically map the flux value to match the PDF of the true hydro spectra. This full iteration can be repeated several times. However, as we will see, one or two full iterations are enough to get excellent agreement with the 1-d power spectrum up to quite high k.

In Fig. 6, we compare the same slice through the hydro flux and through different realisations of LyMAS2 using different combinations of DM fields. The clustering of each pseudo-spectrum is in fair agreement with the clustering of the hydro spectra. Comparing pseudo-spectra and hydro flux along a specific skewer, also shown in Fig. 6, confirms that LyMAS2 correctly models low and high absorptions at good locations, though amplitudes may differ. The second line of Fig. 6 compares three pseudo-spectra generated with LyMAS2 using three different field combinations but using the same seed for the random process to get the variance ∆f k . It's interesting to see that these different pseudo-spectra look also the same, which explains why the slices presented in Fig. 6 are very similar. On the contrary, the third line of Fig. 6 shows four different realisations of pseudo-spectrum from LyMAS2 using the DM overdensity and velocity dispersion fields and different seeds to get the covariance ∆f k . In this case, the amplitude of absorptions can be quite different.

In the next sections, we study in more detail the clustering statistics of each catalog of pseudo-spectra produced with LyMAS2. We aim at recovering three observationally relevant statistics of the transmitted flux: the probability density function (PDF), the line-of-sight power spectrum and the 3d clustering (through the 2point correlation function). As we will see, both LyMAS and LyMAS2 reproduce the PDF of the hydro simulations by construction and nearly reproduce the hydro simulations line-of-sight power spectrum by construction (step 6 above). The power of LyMAS is to produce accurate large-scale 3-d clustering while also reproducing these line-of-sight statistics.

3d-clustering

In order to compare the 3d clustering between the hydro and pseudo spectra, we rely on the 2-point correlation function ξ(r) defined by

ξ(r) = F (x)F (x + r) F 2 -1 , (15)
as a function of the separation r. To study the effect of redshift distortions, we also consider the 2-point correlation function averaged over bin of angle µ defined for a pair of pixels (i,j) by (rirj) /r, where r = |ri -rj| and (rirj) the component along the line of sight.

The top panel of Fig. 7 shows the full 2-point correlation functions derived from pseudo-spectra using either the first version of LyMAS (red line) or LyMAS2 using the DM overdensity and velocity dispersion field (blue line). Compared to the results of the hydro flux (green line), one can see that LyMAS2 is significantly improving the predictions which are remarkably close to the hydro spectra results. In order to estimate the precision of these reconstructions, we plot in the bottom panel of Fig. 7 the relative difference, i.e. ξ/ξ hydro -1 for different combinations of DM fields. It appears clearly that LyMAS2 leads in general to much more accurate predictions than LyMAS. Indeed, when considering the DM overdensity field only, LyMAS2(ρ) give slightly better results (magenta line) but the addition of the velocity dispersion lead to errors that are generally lower than 10% and close to 5% (e.g. blue or black lines). Also, similar trends are obtained when the vorticity is taken into account (bottom panel, green line), which is not surprising as these two fields are highly correlated. On 4) for each LyMAS2 realisation. We then note that the difference between the different pseudo-spectra is quite weak. In the bottom panel, we consider LyMAS2(ρ,σ) to produce different realisations of pseudo-spectra using this time different seeds to get the covariance. In this case, the difference between pseudo-spectra can be much more pronounced.

the contrary, the 1d and 3d velocity divergence, when associated to the DM overdensity only (cyan lines), do not seem to improve much the predictions as respect to the first version of LyMAS. Note that in linear theory, the 3-d velocity divergence is fully correlated with the density field and therefore adds no additional information. On the other hand, the vorticity and/or velocity dispersion are sourced by the non-linear evolution of the matter fields and therefore add complementary information on small scales. Finally, we also note that errors are close to 20% at r ∼1-2 Mpc/h probably due to effect of smoothing.

We now investigate how the predictions of the two point correlation functions vary when considering an angle µ. The trends are presented in Fig. 8 for three ranges of values (0.1 < µ < 0.5, 0.5 < µ < 0.8 and 0.8 < µ) following P14. The results confirm that LyMAS2 significantly improve the predictions of the Ly-α clustering.

In particular, some combinations such as (ρ,σ) still lead to errors generally lower than 10% and most of the time close to 5%. We also note that LyMAS2 is particularly efficient for reproducing the correlations along transverse separations or high angles (i.e. µ >0.8) in which the error is most of the time lower than 5%. The top panels of In particular the DM overdensity field associated to the DM velocity dispersion (blue line) or the vorticity (green line) leads to relative difference lower than 10% and close to 5% in most of the range we are interested in. On the contrary, using the velocity divergence (cyan lines) does not seem to improve much the results.

LyMAS (ρ) LyMAS2 (ρ) LyMAS2 (ρ,Ω) LyMAS2 (ρ,σ) LyMAS2 (ρ, 1d ) ∆ LyMAS2 (ρ, 3d ) ∆ LyMAS2 (ρ,σ,Ω) LyMAS2 (ρ,σ, 1d ) LyMAS2 (ρ,σ, 3d ) ∆ ∆ r (Mpc/h) ξ(r)/ξ hydro (r)-1
hydro flux and those derived from pseudo-spectra produced from LyMAS2(ρ,σ) even for high angles µ >0.8) where the previous version of LyMAS is quite inaccurate. For the two larget µ bins, the correlation functions eventually drop rapidly to zero at large r. In this regime, the fractional error in ξ are inevitably large, even though the absolute errors are small. It is evident that LyMAS2 captures the scale of these zero-crossing more accurately than LyMAS. As a first conclusion, the LyMAS2 scheme is significantly improving the predictions of the Ly-α 3d clustering especially when the DM overdensity field is associated with the velocity dispersion or the vorticity field.

1-d flux power spectrum along LOS

We also aim at producing catalogs of pseudo-spectra that look like spectra measured by a specific instrument i.e. BOSS in the present study. Therefore, the 1-d flux power spectrum of each LyMAS mock should be as close as possible to the hydro spectra 1-d Pk. In the following, we only present the results derived from LyMAS2(ρ,σ) as same trends are obtained when considering any other combination of DM fields. Fig. 9 shows the dimensionless 1-d power spectrum before power spectrum transformation (red line) and after applying the power spectrum and PDF transformation described in the text (black line). We first note that LyMAS2 without iteration reproduces the 1-d power spectrum more accurately than original LyMAS (see Fig. 13 of P14). Then as expected, the 1d-P k transformation leads to same power spectrum as the hydro simulation (blue line), by construction. The second step of the iteration is to re-scale the flux PDF, and this transformation slightly alters the 1d power spectrum. However, as illustrated in the bottom part of Fig. 9, the relative difference is close to 2% up to high values of k (i.e k ∼2 h Mpc -1 ). If one repeats a full iteration a second time, the same accuracy is reached for even higher k values (∼4 h Mpc -1 ).

1-point PDF of the Flux

Since the LyMAS scheme ends after a flux PDF rescaling (second step of a full iteration), this ensures that the 1-point PDF of the hydro flux and the pseudo spectra match exactly. To illustrate, we present in Fig. 10 the results obtained with LyMAS2(ρ,σ) before (red line) and after (blue line) a full iteration (1-d Pk and PDF re-scaling). Without transformation, the PDF of the pseudo-spectra is already close to the PDF of the hydro spectra (green line). But fractional errors on the PDF can be quite large for low or high values of F . After the 1-d Pk re-scaling, the PDF of the pseudo-spectra has slightly changed and has some non-physical values ( F < 0 and F > 1). Then, the second step of the iteration corrects the PDF to exactly match that of the hydro spectra.

LyMAS (ρ) LyMAS2 (ρ) LyMAS2 (ρ,Ω) LyMAS2 (ρ,σ) LyMAS2 (ρ, 1d) ∆ LyMAS2 (ρ, 3d) ∆ LyMAS2 (ρ,σ,Ω) LyMAS2 (ρ,σ, 1d) LyMAS2 (ρ,σ, 3d) ∆ ∆ µ 0.5 < < 0.8 µ 0.1 < < 0.5 µ 0.8 < ξ( ) r Relative Diff.
r (Mpc/h) Figure 8. Same as Fig. 7 but with a dependency to range of angles µ, as defined in the text. These plots confirm that the association of the DM overdensity and the velocity dispersion (or vorticity) lead to very accurate predictions (see blue or green lines). The corresponding relative differences shown in the bottom panel are generally within 10% and 5% over a wide range of r. It is also clear that the new LyMAS2 scheme is much more accurate for large angles where the traditional LyMAS lead to errors that grows quickly with the distance r (red lines). DM fields are again smoothed at = 0.5 Mph/h.

Comparison with the FGPA

The FGPA essentially converts DM density to optical depth using a physical model motivated by photoionization equilibrium, assuming that all gas contributing to the Ly-α lies on a temperature-density relation T ∝ (ρs/ρs) γ-1 . The predicted flux is:

F = Ae -(ρs/ρs) 2-0.7(γ-1) , ( 16
)
where 2 -0.7(γ -1) ≈ 0.6 for the values of γ expected well after reionization [START_REF] Weinberg | [END_REF]Croft et al. 1998;Peeples et al. 2010;McQuinn 2009). This relation is reasonable for modeling high-resolution spectra. However, due to existing non-linear relation between flux and optical depth, it does not automatically apply at low resolution (though it omits some physical effects in the high-resolution case). From the Horizon-noAGN DM overdensity grid smoothed at 0.5 Mpc/h, we have first generated 1024×1024 pseudospectra using equation ( 16) by estimating A so that F = 0.795. Then, we one-dimensionally smoothed each pseudo spectrum to BOSS resolution. Similarly to the LyMAS scheme, we end the process by rescaling the flux 1d-Pk and PDF. The correlation function is shown in the top panel of Fig. 7 and is considerably overes-timated as respect to the hydro flux with a relative error greater than 50% (omitted in the bottom panel for the sake of clarity). Such a trend is consistent with the results of Sorini et al. (2016) who found that typical relative errors in the 3-d power spectrum are ∼80% when a DM smoothing scale of 0.4 Mpc/h is considered.

In Appendix B, we will investigate other deterministic mapping than the FGPA. But will we see that the main conclusion remains unchanged: deterministic sampling generally tends to significantly overestimate the flux 3dcorrelation especially when the dark matter density is smoothed to scales greater than 0.3Mpc/h. Note that a similar trend is obtained when studying the correlation between the Ly-α transmitted flux and the mass overdensity (see Fig. 1 of Cai et al. (2016)).

APPLICATION TO LARGE COSMOLOGICAL DARK MATTER SIMULATIONS

In this section we apply our LyMAS2 scheme to large cosmological DM simulations to produce ensembles of BOSS pseudo-spectra. We first ran several cosmological N-body simulations using Gadget2 (Springel 2005) . Top panel: the dimensionless 1-d power spectrum of the true spectra from the Horizon-noAGN simulation (green line) at z = 2.5 and from coherent pseudo spectra using LyMAS2 considering the DM overdensity and velocity dispersion fields. We show results before (red line) and after (blue and black lines) 1d power spectrum and PDF transformations (described in the text). Bottom panel: the relative difference as respect to the hydro results (i.e. P k /P k,hydro -1). A full iteration (flux 1d-Pk and PDF rescaling) permits to recover the hydro power spectrum with an error of ∼ 2% over a wide range of k. The light and dark grey shaded areas indicate regions where the error is less than 5 and 2% respectively.

with two different box lengths, 1.0 and 1.5 Gpc/h, using the same cosmological parameters as Horizon-noAGN.

In each simulation, the adopted value of the Plummerequivalent force softening is 5% of the mean interparticle distance (24.4 kpc/h and 36.6 kpc/h for the 1.0 Gpc/h and 1.5 Gpc/h boxside respectively) and kept constant in comoving units. . Top panel: the PDF of the true spectra from the Horizon-noAGN simulation at z = 2.5 (green line) and from coherent pseudo spectra using LyMAS2 considering the DM density and velocity dispersion fields. We show results before (red line) and after (blue and black lines) 1-d power spectrum and PDF transformations (described in the text). The 1d Pk re-scaling can lead to non-physical Flux values (i.e. F < 0 or F > 1). Bottom panel: the relative difference with respect to the hydro results (i.e. PDF/PDF hydro -1). The full scheme permits to recover the hydro flux PDF exactly. The light and dark grey shaded areas indicate regions where the error is less than 10 and 5% respectively.

From each cosmological simulation, the corresponding DM density and velocity dispersion fields are computed and sampled on grids of 4096 3 pixels. This allows us to smooth each 1.0 Gpc/h field to 0.5 Mpc/h and each 1.5 Gpc/h one to 1.0 Mpc/h. According to section 3 and appendix A, the combination of the DM overdensity and 4)), derived from the Horizon-noAGN simulation at z = 2.5 and considering the overdensity and velocity dispersion in the Wiener filtering (green dashed line). The resolution of the flux and DM grids are 1024 3 . The red and blue lines are interpolations and extrapolation of the green line to construct the corresponding 1-d Pk for larger 1.0 and 1.5 Gpc/h boxes (using grids of resolution 4096 3 ). The grey shaded area indicates the common k range between the 100 Mpc/h and 1 Gpc/h boxes.

velocity dispersion fields leads to accurate and robust Ly-α clustering predictions. We therefore produce our fiducial large BOSS pseudo-spectra with LyMAS2(ρ,σ). Note that we choose the velocity dispersion field instead of the vorticity mainly for practical reasons, as the computational and memory costs to compute the latter on a large regular grid is much higher.

Once the different dark matter fields are extracted and smoothed to the appropriate scales, the last inputs we need are the relevant transfer functions T defined in equation ( 10) whose detailed expressions can be found in equation ( 11) and equation ( 12) for the 1d or 2d case respectively. We also need the corresponding 1d power spectrum P k to generate the covariance ∆f k at the considered box side. Since the calibrations are derived from the Horizon-noAGN simulation, one potential issue arises from the hydrodynamical box being much smaller than the large DM simulations, so that lower modes are not represented. For the missing modes (k 2π/100), we have extrapolated the values of T1, T2 and P k , while in the common k range, we have proceeded with interpolations. As an illustration, Fig. 11 shows the 1-d power spectrum required to compute the covariance ∆f k when considering the DM density and velocity fields extracted from the 100 Mpc/h hydrodynamical simulation as well as the resulting 1-d power spectra when considering a 1.0 or 1.5 Gpc/h box side.

Fig. 12 illustrates a reconstruction of pseudospectra from a given slice of 4096×4096 pixels through a 1 Gpc/h box simulation. It appears clearly that the 2-d clustering of the pseudo-spectra agrees very well with the clustering of the DM overdensity field. Another visual inspection of an individual skewer also shows that peaks of density match with high absorption. It is also interesting to see that the specific skewer shown in Fig. 12 has in its center a large absorption that corresponds to a large and high density region. Note that the study of groups of so called "Coherently Strong Ly-α Absorption" (CoSLA) systems imprinted in the absorption spectra of a number of quasars (from e.g. BOSS) is of particular interest, as they can potentially detect and trace high redshift proto-clusters (see, for instance Francis & Hewett 1993;Cai et al. 2016;Lee et al. 2018;Shi et al. 2021).

Fig. 13 shows the dimensionless 1-d power spectrum of the pseudo-spectra from a 1 Gpc/h simulation before and after iterations. First, in the common k-range area between the hydro and the pseudo spectra, we find similar trends to those obtained when applying LyMAS2 to the Horizon-noAGN simulation (see Fig. 9). For instance, after two full iterations, the relative difference is close to 2% even for high values of k (∼4 h Mpc -1 ), and similar results are obtained for the 1.5 Gpc/h pseudo spectra. For lower values of k (k 2π/100 ∼ 0.0628 hMpc -1 ), the power spectrum seems to have a natural and consistent extension from the hydro spectra power spectrum. Note also that the highest values of k for the 1.0 and 1.5 Gpc/h box side simulations and grid of resolution 4096 3 are respectively 12.87 and 8.58 Mpc -1 , lower than (2π/100) × 512 ∼ 32.17 hMpc -1 for the hydro box. But the power at these high values is negligible, and missing them in the calculations will not have a noticeable impact on spectra. Also, since a full iteration ends with a Flux PDF re-scaling, this ensures exactly match to the PDF of the hydro flux.

Fig. 14 shows the 2-point correlation functions derived from several large-scale pseudo-spectra. In particular, we show the predictions from the first version of LyMAS (red lines) and those obtained from LyMAS2 using the DM density field only (black lines) and with additional velocity dispersion field (blue lines). We also add the predictions derived from the 1.5 Gpc/h simulation (magenta lines), using again LyMAS2(ρ,σ). These plots confirm first that the traditional LyMAS (red lines) tends to overestimate the correlations and this trend is more pronounced when considering high angles (µ > 0.8), as already noted in section 3. The result is quite similar with the LyMAS2 scheme when considering the DM overdensity only. However, the difference from LyMAS is more and more noticeable as µ increases. These results are again consistent with those presented in section 3. The difference becomes even stronger when adding the DM velocity dispersion field. In this case, LyMAS2(ρ,σ) tends to significantly reduce the correlations and most probably lead to more reliable predictions. In the range 2 r 10 Mpc/h, the correlations are very close to those of the hydro simulation. It is also impressive that the 1.5 Gpc/h mock generated with LyMAS2(ρ,σ) leads to very similar trends (for separations r 2 Mpc/h), though the DM fields are now smoothed to 1.0 Mpc/h. This success is consistent with the results presented in the appendix A, where we compare the performance of LyMAS2 using different DM smoothing scales. This robustness is one of the key improvements accomplished with LyMAS2.

Finally, we show in Fig. 15 the 2-point correlation function averaged from five different realizations of 1 Gpc/h Ly-α pseudo-spectra obtained by applying LyMAS2(ρ,σ) to different DM cosmological simulations, at z = 2.5. The plots show clear features of BAO at r ∼105 Mpc/h and variations with respect to the angle µ, consistent with observational trends (see for instance du Mas des Bourboux et al. ( 2020)) This illustrates the ability of LyMAS to properly describe redshift distor-sions and to model realistic large BOSS Ly-α forest spectra catalogues.

CONCLUSIONS

We have introduced LyMAS2, an improved version of the LyMAS scheme (Peirani et al. 2014). In this new version, we have used the Horizon-noAGN (Peirani et al. 2017) simulation to characterize the relevant cross-correlations between the transmitted flux and the different DM fields. In particular, we have considered not only the DM overdensity but also specific DM velocity fields (i.e. velocity dispersion, vorticity, 1d and 3d divergence) and used Wiener filtering to generate the specific calibrations. LyMAS2 shares the same philosophy as LyMAS that flux correlations are mainly driven by the correlations of the underlying DM (over)density, and it uses additional information from the DM velocity correlations to refine the theoretical predictions. In a second step, we have applied LyMAS2 to DM fields extracted from the hydrodynamical or large DM-only simulations to create large ensembles of pseudo-spectra with redshift distortions, at z = 2.5 and at the BOSS resolution. Throughout the analysis, we use a DM smoothing of 0.5 Mpc/h to derive the main trends and results. Our main conclusions can be summarized as follows:

• LyMAS2 greatly improves the predictions for flux statistics of the 3-d Ly-α forest on small and large scales. More specifically, we found that the DM overdensity combined with the DM velocity dispersion (or the vorticity) recovers the 2-point correlation functions of the (reference) hydro flux within 10% and (most of the time within 5%) even when high angles are considered. This is a major improvement with respect to the original version of LyMAS, which is rather inaccurate in predicting the Ly-α correlations for large separations and high angles.

• Like LyMAS, LyMAS2 reproduces the 1-point PDF of the flux from the calibrating hydro simulation exactly, by construction. It also reproduces the 1-d (lineof-sight) power spectrum with en error of about 2% up high k values. The LyMAS2 pseudo spectra therefore have realistic observable properties on small scales while also having accurate large scale 3-d clustering when applied to a large volume DM-only simulation.

• The trends derived from five different and slightly lower resolution hydrodynamical simulations are consistent with those obtained from the fiducial Horizon-noAGN simulation. This suggests that the results presented in this study are robust. Moreover, this allows us to estimate error bars on the 2-point correlations functions, which are generally low.

• We have considered three different DM smoothing scales (0.3, 0.5 and 1.0 Mpc/h) and found similar trends in the flux clustering predictions. It is encouraging that a DM smoothing of 1.0 Mpc/h still leads to very accurate predictions, especially in the 2-point correlation functions even at high angles and large separation. Indeed, the errors are typically lower than 5%, whereas they are generally higher than 30% with the original version of LyMAS.

• LyMAS2 applied to large DM cosmological simulations of box side either 1.0 or 1.5 Gpc/h indicates that the predicted flux statistics follow the same trends obtained from the (100 Mpc/h) Horizon-noAGN DM fields. Indeed, we found again that the first version of LyMAS tends to overestimate the flux correlations at large separations and/or at high angles. On the contrary, LyMAS2 using for instance the DM overdensity and the velocity dispersion clearly reduces the 2-point correlation functions to lead to more reliable and accurate predictions. Moreover LyMAS2 adequately models large scale Ly-α absorptions systems which correspond to massive overdensity regions.

• Deterministic mappings such as the Fluctuating Gunn-Peterson Approximation tend to considerably overestimate the 3-d flux correlations especially at large separation or when high angles are considered.

LyMAS2 offers a sophisticated tool to accurately model and predict large scale Ly-α forest 3d statistics. This opens new opportunities to improve diversified studies such as Ly-α forest cross-correlation (e.g. Lochhaas et al. 2016), two-point correlations or threepoint correlations analysis (e.g. Tie et al. 2019) or BAO feature predictions. Moreover, large Ly-α catalogs produced with LyMAS2 can be used to characterize massive overdensity regions such as proto-clusters through groups of coherent large absorptions analysis (Cai et al. 2016;Shi et al. 2021;Lee et al. 2018). Compared to previous work, we recall that the main objective of LyMAS is to create large Ly-α mocks for a specific instrument (here BOSS) with 3d flux statistics as close as possible to those that would be obtained from a very large volume (but computationally intractable) hydrodynamical simulation.

The Iteratively Matched Statistics (IMS) developped by Sorini et al. (2016) does not present such predictions and limits their analysis to small simulation boxes ( 100 Mpc/h). However, when comparing the flux PDF and 1d power spectrum, LyMAS2 and 1D-IMS (see introduction) lead to similar performances: the 1D-IMS scheme perfectly reproduces these statistics, while errors of ∼2% are obtained with LyMAS2 for the flux 1d-Pk. Regarding the 3D-IMS scheme, errors are much higher, of order of 15 and 20% respectively for the flux PDF and 1d power spectrum. As far as the 3d flux statistics are concerned, at a DM smoothing of 0.4 Mpc/h, the 1D-IMS and 3D-IMS present errors of 20% and 10-20% (for a DM smoothing of 0.4 Mpc/h) respectively regarding the reconstruction of the power spectrum. In this study, LyMAS2 mainly considers a DM smoothing of 0.5 Mpc/h, which leads to errors generally lower than 5% for the 2-point correlation functions. Again, it is worth mentioning that similar (low) errors are also obtained with LyMAS2 when considering a DM smoothing of 1.0 Mpc/h. It would be then interesting to compare the performance of the 1D and 3D-IMS scheme at this specific smoothing scale in the perspective of creating large ( 1.0 Gpc/h) Ly-α mocks.

Recently, Harrington et al. (2021) have trained a convolutional neural network from hydrodynamical simulations of side 20 Mpc/h to predict both the density, the temperature and the velocity fields. This method is quite flexible and the predictions of the flux PDF and 1-d power spectrum (i.e. within ∼5% up to k∼10 Mpc/h) are promising and more accurate than the FGPA. Similarly, Sinigaglia et al. (2021) has developed a new physically-motivated supervised machine learning method (HYDRO-BAM) from a reference hydrodynamical simulation of comoving side 100 Mpc/h. The PDF, power spectrum and bi-spectra can be reconstructed with error of a few percent up to modes k= 0.9 Mpc/h. Again, it would be very interesting to see the perfor-mance of such approaches on the 3-d correlations and when considering smoothed spectra and larger boxes.

Improvements can still be done in the LyMAS scheme. For instance, one main assumption is to consider that the transverse correlations are mainly driven by the effect of DM smoothing. In the present study, we stress again that all the approach is based on creating pseudo-spectra individually and independently from each other. Because the draws of ∆f k are independent on each line-of-sight, spectra at small but nonzero transverse separations can look quite different on small scales. Since the predictions on the clustering of the flux are already very accurate with LyMAS2, we haven't considered the same approach in volume. This would take into account transverse correlations between lines of sight that have been neglected in this work: instead of predicting the flux from dark matter fields independently for each LOS, one would predict the entire cube of flux from the dark matter field cubes, using the full 3-d covariance structure. One would still use an assumption of spatial homogeneity (stationarity), so that 3D Fourier space coefficients could be computed independently, however one would need to take care of the statistical anisotropy in the LOS direction, therefore all statistics in Fourier space would depend on |k ⊥ | and k . Taking into account transverse correlations would thus further reduce the covariance of the flux conditionally to the dark matter fields, in other words reduce the noise in the predicted flux field. Among future prospects, we plan to extend this work to predict the flux clustering for other surveys such as the Dark Energy Spectroscopic Instrument (DESI, DESI Collaboration et al. 2016), the William Herschel Telescope Enhanced Area Velocity Explorer (WEAVE-QSO, [START_REF] Pieri | Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics[END_REF] or Subaru Prime Focus Spectrograph (PFS, Takada et al. 2014). They will open new vistas on the high redshift intergalactic medium probed by the Ly-α forest. It would be then interesting to estimate the level of performance of LyMAS2 when the transmitted flux has a higher resolution than BOSS spectra, which might require reducing the DM smoothing. Finally, we also intend to use Machine Learning in the process to see whether we can still improve the predicted flux statistics (Chopitan et al. 2021).

APPENDIX A: GENERAL TRENDS

In section 4.1, we have presented the predictions regarding the 2-point correlation functions for height different combinations of DM fields, with calibrations derived from Horizon-noAGN. To check the robustness of the results, the analysis of other similar hydrodynamical simulations is definitely required. To limit the computational time, we ran five additional hydrodynamical simulations with the same box side and same physics than Horizon-noAGN but with two times lower resolution (i.e. 512 3 DM particles instead of 1024 3 and a minimal cell size of ∆x=2 kpc instead of 1 kpc). The first simulation uses degraded Horizon-noAGN initial conditions while the other ones have different initial phases. For each of the five new simulations, we generated the corresponding grids of transmitted flux, DM overdensity and velocity fields and calibrations folowing the same methodology presented in section 3. We consider here flux and DM fields sampled on grids of 512×512×1024 namely 512×512 spectra of resolution 1024.

In the first step, we consider all DM fields smoothed at 0.5 Mpc/h. After checking first that the "high" and "low" resolution Horizon-noAGN simulations give consistent trends, we took an interest in the variations of then mean of the absolute relative difference (1/5)| 5 i=1 (ξi/ξ hydro,i -1)|, where we compare the 2point correlation function of the hydro spectra ξ hydro,i from a given simulation "i" to those derived from pseudo-spectra generated with LyMAS2 ξ hydro,i . In Fig. A1, we summarize the results obtained with the original LyMAS and LyMAS2 considering the same DM field combinations than in Figures 7 and8. The main conclusion is that we do find very similar trends than those obtained with Horizon-noAGN , which strongly suggest that our results are robust. In particular, the use of the velocity dispersion (σ) or the vorticity (Ω) lead to relative errors that are remarkably low i.e. in general lower than 5% even for the different ranges of angle. The plots also confirm that the 1-d and 3-d velocity divergence fields do no permit to reach the same level of accuracy.

In the next step, we present the trends obtained when the DM fields are smoothed to 0.3 or 1.0 Mpc/h. We only present in Fig. A2 the results for LyMAS, LyMAS2(ρ,σ) and LyMAS2(ρ,σ,Ω) to have a clear overview of the general trends. In P14, we found that a DM smoothing of 0.3 Mpc/h was an optimal value to reach the highest accuracy in the predictions. This is confirmed here since we get errors of 10% compared to 20% and 30% with values 0.5 and 1.0 Mpc/h respectively. As expected, LyMAS2 permits to reduces such errors that are in general much lower than 10% and most of the time lower than 5%. It is also very promising that LyMAS2 applied to DM fields smoothed at 1 Mpc/h gives such accurate predictions even for high values of µ. This is definitely not the case with the original LyMAS leading to very high errors. Note also that due the smoothing scale, the predictions are less accurate for distance lower than 2 Mpc/h but acceptable for large scale analysis.

APPENDIX B: DETERMINISTIC MAPPING

One commonly way to produce large mocks of Lyα forest, from Gaussian fields or DM distributions extracted from cosmological simulations, is to use a physically motivated deterministic relation which links the Ly-α optical depth (or transmitted flux) to the dark matter overdensity. This is the case with the so-called Fluctuating Gunn-Peterson Approximation which has been extensively used in the literature. However, the FGPA is supposed to be more suitable for modelling high-resolution spectra and can be strongly limited when the DM density field is smoothed to a scale greater than 0.1 Mpc/h (see for instance the analysis of Sorini et al. (2016)) and confirmed by our results in section 4.4. For this reason, we have derived in P14 an "optimal" deterministic relations by matching the corresponding cumulative distributions of the smoothed transmitted Flux Fs and dark matter overdensity ρs as:

Fs 0 P (F ′ s )dF ′ s = ∞ ρs P (ρ ′ s )dρ ′ s ,
where P (Fs) and P (ρs) are the 1-point PDFs of the flux and DM overdensity measure from the simulation. One advantage of choosing such deterministic relation is to recover by construction the PDF of the hydro flux. However, the 2-point correlation function of pseudo-spectra generated with this approach is still highly overestimated (see for instance Figures 10 and 19 in P14).

In this section, we consider another choices of deterministic relations. In particular, Tie et al. ( 2019) have used the conditional probability P (F |1+δ) of the transmitted flux on the DM overdensity to get the conditional mean flux:

F (1 + δ) = F.P (F |1 + δ)dF . (B1)
It can be analytically demonstrated that the 2-point correlation function of pseudo-spectra obtained from such a determinisic mapping is the same that the one obtained with the first version of LyMAS (Tie et al. 2019). Since this deterministic relation can be easily extended to several DM fields, our aim is to investigate whether the inclusion of different DM velocity fields in such deterministic mappings may improve the trends or not. As an illustration, Fig. B1 shows examples of 1ddeterministic relations constructed for different smoothing of the DM overdensity field and and one example of 2d-deterministic relation F (1 + δ, σ) using both the DM overdensity and the velocity dispersion fields (smoothed at 0.5 Mpc/h). All relations are derived from the Horizon-noAGN simulation at z = 2.5 (in redshift space). In principle, the 2d-deterministic relation is supposed to refine the results as repect to the 1d deterministic one. Indeed, let's take for instance a DM overdensity of 1 + δ = 1. This leads to an unique mean flux of F =0.75 from the 1-d deterministic relation (using a DM smoothing of 0.5 Mpc/h). The 2d-determinisic relation provides, however, a wide range of possible values of F depending this time on the velocity dispersion (see figure B2).

We have then produced grids of pseudospectra from DM fields extracted from the Horizon-noAGN simulation (see section 3), smoothed at 0.5 Mpc/h and using three different deterministic relations. The first one considers the DM overdensity field only (1-field), the second one both overdensity and velocity dispersion fields (2-fields) while the last one associates the DM density field to the velocity dispersion and vorticity fields (3-fields). To estimate the mean value of the flux from a given value of ρ or a given set of (ρ,σ) or (ρ,σ,Ω), we use respectively interpolations, bilinear interpolations and trilinear interpolations depending of the number of input DM fields. First, Fig. B3 shows the 1-d power spectrum and PDF of pseudo spectra (without iteration) for the 2-fields case. We notice that the 1-point PDF is in general not well recovered. The predictions of the 1d-P k are also not as accurate than LyMAS and things especially for the 2-fields and 3-fields cases. Indeed, the power spectra at small scales are considerably overestimated. Although a full iteration can improved these trends, the predicted 2 point correlation functions, shown in Fig. B4, present errors that are generally quite high especially when an angle µ is considered. For instance the errors are much higher than 10% when µ > 0.8. Also, one main drawback of this approach when creating large mocks (>1 Gpc/h) is to have enough statistics from the hydrodynamics simulations to cover most of the parameter space of the large cosmological simulation. This should not be an issue for the 1-field case since simple interpolations and extrapolations can be done (e.g. from Fig. B1 F ∼1 and F ∼0 for 1 + δ < -1.5 and 1 + δ > 1.5 respectively). For the 2-fields and the 3-fields cases, efficient interpolations and extrapolations could be obviously much more complicated to realized.

APPENDIX C: ADAPTIVE SMOOTHING

In this section we provide some detail on how the dark matter density field, velocity field and velocity dispersion are interpolated adaptively on a mesh from the output of a cosmological dark matter N -body simulation, before further treatment by LyMAS2, in particular additional smoothing with a Gaussian window.

Note that, while the notations below assume standard configuration space, the calculation naturally extend to redshift space, by just modifying particle coordinates accordingly with the local peculiar velocity component contribution. In LyMAS2, we use the infinitely remote observer approximation by just accounting for redshift distortion along z axis.

For a smooth phase-space distribution function f (x, v), the projected density is given by ρE

(x) = d 3 v f (x, v), (C1)
the Eulerian mean velocity field by

vE(x) = v E = 1 ρE(x) d 3 v v f (x, v), (C2)
and the local mean square velocity reads

v 2 E(x) = 1 ρE(x) d 3 v v 2 f (x, v). (C3)
Obviously, the last two equations stand for points of space where ρ(x) > 0. From equation (C3), we can derive the local velocity dispersion

σ 2 v,E ≡ v 2 E -v 2 E . (C4)
We notice that the local velocity field can be considered as a statistical average, this is why we used the

• • • E notation above, v E(x) = d 3 v v fE(x, v), (C5)
with the local density probability

fE(x, v) ≡ 1 ρE(x) f (x, v). (C6)
In this probabilistic approach, the mean square velocity is given by equation (C3), since

v 2 E(x) = d 3 v v 2 fE(x, v). (C7)
and, likewise, its local variance by equation ( C4).

What we have actually access to is not a smooth distribution function, unfortunately, but a distribution of N simulation particles of individual masses mi, positions xi and velocities vi. This means that the phasespace distribution function has the following form

f (x, v) = i mi δD(v -vi) δD(x -xi), ( C8
)
where δD is the Dirac distribution function. From eq. ( C8), one can compute the projected Eulerian density

ρE(x) = i mi δD(x -xi), (C9)
but the Eulerian velocity field vE(x) and mean square velocity are ill-defined. However, the underlying distribution of true dark matter particles is much smoother than its crude numerical representation in terms of macroparticles of the N -body simulation. While, strictly speaking, the phasespace density is still of the form (C8) at the microscopic level, it can be considered as a smooth function at the macroscopic level, at least in terms of probability density.

In order to recover a good approximation of the continuum, Colombi et al. (2007) proposed a locally adaptive smoothing algorithm, SmoothDens, inspired from smooth particle hydrodynamics (hereafter SPH), using, to compute various fiels, an interpolation window Fx, of which the shape parameters, in particular the typical size ℓ(x), depend on position x. This window function is normalized to unity, i.e. d 3 x ′ Fx(x ′ ) = 1.

In principle, for a given function h(x, v), the smoothed counterpart is given by

[Fx * h](x, v) = d 3 x ′ Fx(x -x ′ ) h(x ′ , v). (C10)
Setting h = ρE, after simple algebraic calculations exploiting the properties of the Dirac distribution function, we obtain the simple expression for the adaptively smoothed density:

ρF(x) = i mi Fx(x -xi). (C11)
From there, we can formally define the analogous of equation (C6) but with adaptive smoothing performed in the spatial position,

fF(x, v) ≡ 1 ρF(x) [Fx * f ](x, v), (C12)
from which one can derive estimates in the mean field limit of velocity related quantities: ) In the algorithm SmoothDens, the adaptive procedure is used to compute various fields for a particular set of positions x = xj on a cubical grid of size ng. Function F is a compact Monaghan & Lattanzio (1985) spline of size ℓ(x), with ℓ(x) the distance of the N th SPH closest simulation particle to position x. The value of NSPH we adopt here is NSPH = 32. As an additional recipe, adaptive smoothing is locally replaced with nearest grid point (NGP) interpolation [START_REF] Hockney | Computer simulation using particles Japelj[END_REF] when ℓ(x) is smaller or of the order of the grid cell size L/ng, where L is the simulation box size. Also, a local weight is given to each particle i so that at the end, its total contribution sums up to the particle mass mi. Note that, due to the finite extension of the spline function, some particles belonging to dense clusters or close to dense clusters may not contribute at all. In the latest implementation of the algorithm, SmoothDens5, which we use in LyMAS2, an option allows one to affect these particles to the grid with NGP interpolation in order to conserve total mass. The effect of not doing so is however generally small.

vF(x) = i mi vi Fx(x -xi) i mi Fx(x -xi) , (C13) v 2 G = i mi v 2 i Fx(x -xi) i mi Fx(x -xi) . ( C14
The outcome of SmoothDens mainly depends on two parameters, the resolution ng of the grid and the value used for the number of neighbors, NSPH. Changing both these parameters can have drastic impact on the results, especially the local velocity dispersion and the velocity derivatives estimates. Additional Gaussian smoothing performed in LyMAS2 is however expected to reduce considerably the dependence on these two parameters provided that the smoothing scale RG associated to the smoothing window is large enough compared to L/ng. Yet one has to bear in mind that the influence of NSPH cannot be negligible in underdense regions as long as it can influence scales larger than RG, which is unfortunately very likely. Despite these non trivial issues, the reason why LyMAS2 still works so accurately is that it is calibrated relying on probability distributions mappings, which naturally corrects for intrinsic biases introduced by local adaptive smoothing. k (h Mpc -1 ) Relative Diff.

Figure 13. Top panel: the dimensionless, redshift-space, 1-d power spectrum of pseudo and true spectra from the Horizon-noAGN simulation (green line) and from coherent pseudo spectra using LyMAS2 considering the DM density and the velocity dispersion fields from a 1 Gpc/h box side N-body simulation. Here again, we show results before (red line) and after (blue line) one full iteration (i.e. flux 1d power spectrum and PDF transformations). We also show the results when repeating a second iteration (black line). The grey shape defines the common k-range between hydro spectra and pseudospecra. Bottom panel: the relative difference with respect to the hydro results (i.e. P k /P k,hydro -1). The light grey and dark grey bands define regions where the error is less than 5 and 2 % respectively. and from 1 Gpc/h pseudo-spectra using the full LyMAS (red line) and LyMAS2 scheme using the DM overdensity field only (black line) or combined with the DM velocity field (blue line). We also show the 2-point correlation function derived from 1.5 Gpc/h pseudo spectra produced with LyMAS2(ρ,σ) (magenta line). In the other panels, we show the corresponding flux correlation functions averaged over bins of angle µ, as labeled. LyMAS tends to overestimate the correlation especially for large separations and high angles. LyMAS2(ρ,σ) tends to significantly reduce such correlations, which suggests more reliable predictions. For the sake of clarity, we only show results for the original LyMAS (red curves), LyMAS2(ρ,σ) (blue curves) and LyMAS2(ρ,σ,Ω) (black curves). Note that even with a DM smoothing of 1.0 Mpc/h, LyMAS2 still permits to reach a high level accuracy, even for high angles (mu > 0.8).

- . The relative difference of the two-point correlation functions of pseudo-spectra produced with deterministic sampling (and a full iteration in the scheme). It appears clearly that the different deterministic sampling do no reach the level of accuracy of LyMAS2 especially for high angles where the error are quite high (> 10%). The light and dark grey shade represent error lower that (absolute) 10 and 5% respectively.

Figure 1 .

 1 Figure1. Slices through the Flux (1-dimensionally smoothed at the BOSS resolution) as well as corresponding DM overdensity and velocity fields (3-dimensionally smoothed at 0.5 Mpc/h) in redshift space (horizontal direction). Each field has been extracted from the Horizon-noAGN simulation at z = 2.5 and has been normalized to help the visual comparison. The dashed lines correspond to the same line of sight (see Figure2).

Figure 2 .

 2 Figure 2. An example showing clear correlations (or anti-correlations) between the evolution of the hydro flux (green line) and DM density and velocity fields (blue and red lines) along the same skewer extracted from the Horizon-noAGN at z = 2.5. These skewers are extracted from the slices studied in Fig. 1.

Figure 3 .

 3 Figure 3. Correlations between the optical depth τ = -lnF in the hydro spectra (smoothed at the BOSS resolution) and dark matter quantities smoothed at 0.5 Mpc/h namely the overdensity (1 + δ), the vorticity (Ω), the velocity dispersion (σ), the 3d-velocity divergence and the 1d-velocity divergence (along the line of sight) at z = 2.5. Colors show the density of pixel using normalized values and contour line mark areas enclosing 68.27% and 95.45%. Hydro spectra and DM fields have been computed in redshift space and extracted from the Horizon-noAGN simulation.

Figure 4 .

 4 Figure 4. Some correlations between different dark matter quantities from the Horizon-noAGN simulation at z = 2.5. All field have been smoothed to 0.5 Mpc/h and have been computed in redshift space. Colors show again the density of pixel using normalized values and contour line mark areas enclosing 68.27% and 95.45%.

  Figure 5. Examples of cross-spectrum P ab (k) = a * k b k where a and b refer either to the hydro flux or a DM field, derived from the Horizon-noAGN simulation at z = 2.5. The flux field is 1-d smoothed at the BOSS resolution while the DM field are smoothed at 0.5 Mpc/h.

Figure 6 .

 6 Figure6. First line: comparison between corresponding slices through the hydro flux (top left panel), 1-dimensionally smoothed at the BOSS resolution, and through pseudo-spectra generated with LyMAS2 and different associations of DM fields (3dimensionally smoothed at 0.5 Mpc/h). The direction of redshift distorsion is horizontal. Flux and DM fields have been extracted from the Horizon-noAGN simulation at z = 2.5. Visual inspection suggests a very good agreement between the clustering of the hydro flux and those of pseudo-spectra. The two bottom panels show the evolution through a specific line of sight (located by the dashed line in the different slices) of the hydro flux (green dotted line) and pseudo-spectra generated with LyMAS2. In the top panel, we use the same seed to randomly generate the covariance from equation (4) for each LyMAS2 realisation. We then note that the difference between the different pseudo-spectra is quite weak. In the bottom panel, we consider LyMAS2(ρ,σ) to produce different realisations of pseudo-spectra using this time different seeds to get the covariance. In this case, the difference between pseudo-spectra can be much more pronounced.

  Fig. 8 indicate again a remarkably good agreement between the 2 point correlation functions of the

Figure 7 .

 7 Figure 7. Top panel: the correlation function ξ as a function of the separation r. We present results from the true hydro spectra (green line) and results from the LyMAS (red line) and the new version LyMAS2 considering the DM overdensity and the velocity dispersion fields. The predictions from FGPA (cyan line) are also shown. DM fields are all smoothed at the scale 0.5 Mpc/h. Bottom panel: the relative difference as respect to the true hydro spectra (i.e. ξ(r)/ξ hydro (r) -1) are shown for different combination of choices in the Wiener filtering. Here all DM fields have been smoothed to = 0.5 Mph/h. Compared to the traditional LyMAS scheme (red lines), the new version significantly improves the predictions.In particular the DM overdensity field associated to the DM velocity dispersion (blue line) or the vorticity (green line) leads to relative difference lower than 10% and close to 5% in most of the range we are interested in. On the contrary, using the velocity divergence (cyan lines) does not seem to improve much the results.

Figure 9

 9 Figure9. Top panel: the dimensionless 1-d power spectrum of the true spectra from the Horizon-noAGN simulation (green line) at z = 2.5 and from coherent pseudo spectra using LyMAS2 considering the DM overdensity and velocity dispersion fields. We show results before (red line) and after (blue and black lines) 1d power spectrum and PDF transformations (described in the text). Bottom panel: the relative difference as respect to the hydro results (i.e. P k /P k,hydro -1). A full iteration (flux 1d-Pk and PDF rescaling) permits to recover the hydro power spectrum with an error of ∼ 2% over a wide range of k. The light and dark grey shaded areas indicate regions where the error is less than 5 and 2% respectively.

Figure 10

 10 Figure10. Top panel: the PDF of the true spectra from the Horizon-noAGN simulation at z = 2.5 (green line) and from coherent pseudo spectra using LyMAS2 considering the DM density and velocity dispersion fields. We show results before (red line) and after (blue and black lines) 1-d power spectrum and PDF transformations (described in the text). The 1d Pk re-scaling can lead to non-physical Flux values (i.e. F < 0 or F > 1). Bottom panel: the relative difference with respect to the hydro results (i.e. PDF/PDF hydro -1). The full scheme permits to recover the hydro flux PDF exactly. The light and dark grey shaded areas indicate regions where the error is less than 10 and 5% respectively.

Figure 11 .

 11 Figure 11. An example of dimensionless 1-d Pk required to compute the covariance ∆f k (see Eq. (4)), derived from the Horizon-noAGN simulation at z = 2.5 and considering the overdensity and velocity dispersion in the Wiener filtering (green dashed line). The resolution of the flux and DM grids are 1024 3 . The red and blue lines are interpolations and extrapolation of the green line to construct the corresponding 1-d Pk for larger 1.0 and 1.5 Gpc/h boxes (using grids of resolution 4096 3 ). The grey shaded area indicates the common k range between the 100 Mpc/h and 1 Gpc/h boxes.

Figure 12 .

 12 Figure 12. Application of the full LyMAS2 scheme to a large dark matter simulation (1 Gpc/h -2048 3 particles) at z = 2.5. From the DM overdensity and velocity dispersion fields, both sampled on 4096 3 regular grids and smoothed at 0.5 Mpc/h, we derive pseudo-sepctra at the BOSS resolution using the full LyMAS2 scheme. The left part of the figure shows corresponding slices which suggest a fair agreement between the clustering of the DM density and the pseudo spectra. The right part shows an individual skewer confirming that DM density peaks are associated to high flux absorption.

Figure 14 .

 14 Figure14. Top left panel: the 2-point correlation functions of the flux derived from the 100 Mpc/h Horizon-noAGN simulation (green dotted line) and from 1 Gpc/h pseudo-spectra using the full LyMAS (red line) and LyMAS2 scheme using the DM overdensity field only (black line) or combined with the DM velocity field (blue line). We also show the 2-point correlation function derived from 1.5 Gpc/h pseudo spectra produced with LyMAS2(ρ,σ) (magenta line). In the other panels, we show the corresponding flux correlation functions averaged over bins of angle µ, as labeled. LyMAS tends to overestimate the correlation especially for large separations and high angles. LyMAS2(ρ,σ) tends to significantly reduce such correlations, which suggests more reliable predictions.

Figure 15 .Figure A1 .Figure A2 .

 15A1A2 Figure 15. The 2-point correlation function averaged from five different realizations of 1 Gpc/h Ly-α pseudo-spectra obtained by applying LyMAS2(ρ,σ) to DM cosmological simulations. The shaded areas represent the error on the mean (rms).

Figure B1 .

 B1 Figure B1. Examples of 1-d and 2-d deterministic relations derived from the Horizon-noAGN simulation at z = 2.5. Top panel: The mean flux < F > only depends on the DM overdensity 1 + δ. Three different DM smoothing have been considered. For example, the mean flux has a value of 0.75 for an overdensity of 1 (DM smoothing=0.5 Mpc/h). Bottom panel: scatter plot showing the mean flux < F > as respect to the DM overdensity and velocity dispersion field (σ). In this case, < F > can have a wide range of values for 1 + δ=1 (see also Figure B2).

Figure B2 .Figure B3 .

 B2B3 Figure B2. The evolution of the mean flux as respect to the velocity dispersion (σ) and for a given overdensity (i.e. 1 + δ=1). The DM fields are smoothed at 0.5 Mpc/h. This variation is directly derived from the 2-d deterministic sampling presented in Figure B1. Compared to the 1-d deterministic relation, a wide range of values of < F > is obtained and should refine the predictions.

Table 1 .

 1 Summary of fields and corresponding notation used in the text.

	Hydro spectra
	Flux (smoothed at BOSS res.)	F
	Optical depth	τ =-ln F
	Dark Matter Fields
	Smoothed density	ρs
	Overdensity	ρ = 1 + δ = ρs/ ρs
	Vorticity	Ω
	Velocity dispersion	σ
	1d Vel. divergence	∇ 1d
	3d Vel. divergence	∇ 3d

Table 2 .

 2 Summary of the different DM field combinations considered in the LyMAS2 scheme

	Name	Field 1	Field 2	Field 3
	LyMAS2(ρ)	Overdens.
	LyMAS2(ρ, σ)	Overdens. Vel.disp.
	LyMAS2(ρ, Ω)	Overdens. Vorticity
	LyMAS2(ρ, ∇ 1d )	Overdens. 1-d div.
	LyMAS2(ρ, ∇ 3d )	Overdens. 3-d div.
	LyMAS2(ρ, σ, Ω)	Overdens. Vel.disp.	Vorticity
	LyMAS2(ρ, σ, ∇ 1d ) Overdens. Vel.disp.	1-d div.
	LyMAS2(ρ, σ, ∇ 3d ) Overdens. Vel.disp.	3-d div.

c 0000 The Authors

MNRAS 000, 000-000 (0000)

Acknowledgements

This work was carried within the framework of the Horizon project (http://www.projet-horizon.fr). Most of the numerical modeling presented here was done on the Horizon cluster at IAP. This work was supported by the Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, cofunded by CEA and CNES. DW acknowledges support of U.S. National Foundation grant AST-2009735. We warmly thank T. Sousbie, B. Wandelt, O. Hahn, M. Buehlmann and S. Rouberol for stimulating discussions. We also thank D. Munro for freely distributing his Yorick programming language (available at http://yorick.sourceforge.net/) which was used during the course of this work.

Data availability

The data and numerical codes underlying this article were produced by the authors. They will be shared on reasonable request to the corresponding author.