The cosmic web from perturbation theory

,


Introduction
The cosmic web is the large-scale structure (LSS) pattern that emerges as a manifestation of the action of gravity in an expanding background universe.It is shaped according to the cosmological information content set in the initial conditions of cosmic times.
The scientific community is putting special effort into mapping the three-dimensional distribution of matter in the Universe through galaxy surveys such as DESI (Levi et al. 2013), Euclid (Amendola et al. 2016), J-PAS (Benitez et al. 2014), and the Nancy Grace Roman Space telescope (Spergel et al. 2015).The clustering of the LSS yields powerful constraints on the standard cosmological model, such as the nature of dark energy (e.g.DESI Collaboration 2016), primordial non-Gaussianities (see, e.g., Meerburg et al. 2019), and neutrino masses (see, e.g., Chudaykin & Ivanov 2019).
Since the observable Universe is unique, one needs mock catalogues of it in order to perform robust analyses.Such catalogues permit the computation of covariance matrices and the study of systematics in the observational data.However, the computation of such a mock catalogue is usually costly, and only a few can be done (see, e.g., Angulo et al. 2012;Fosalba et al. 2015;Chuang et al. 2019).In this context, bias mapping techniques at the field level have emerged as a solution to save enormous computational resources while maintaining high accuracy.Such examples include PATCHY (Kitaura et al. 2014(Kitaura et al. , 2015)), applied to the BOSS data (Kitaura et al. 2016), and EZmocks (Chuang et al. 2015), which was applied to the eBOSS data (Zhao et al. 2021).More recently, the BAM code has been designed to learn the complex bias relation from reference simulations (Balaguera-Antolínez et al. 2019, 2020;Kitaura et al. 2022).These techniques can also accurately map the Lyman-α forest (see, e.g., Sinigaglia et al. 2021Sinigaglia et al. , 2022)).All the methods mentioned above require a dark matter field defined on a mesh.The accuracy of that matter distribution determines the precision of the mock catalogues.
A series of ideas have been implemented to accelerate particle-mesh based N-body solvers (see Tassev et al. 2013;Feng et al. 2016).Despite these developments, N-body codes are costly when aiming to mass-produce mock catalogues covering large cosmic volumes.Therefore, approximate gravity solvers are still commonly used.While EZmocks relies on the Zel'dovich approximation (Zel'dovich 1970), PATCHY and BAM rely on ALPT, including tidal field corrections (Kitaura & Hess 2013).It has been shown that such methods can, to a great extent, correct the bias introduced by approximate gravity solvers in the nonlinear and non-local bias description (see the application of BAM to galaxy catalogues, Balaguera-Antolínez et al. 2023).
Machine learning techniques applied to large training data sets (of thousands of N-body simulations for a given cosmology) have emerged as an alternative to approximate gravity solvers (He et al. 2019;Dai & Seljak 2021;Jamieson et al. 2023).(For a review of numerical methods in LSS modelling, see Angulo & Hahn 2022) This work presents a novel approach to modelling the cosmic web by extending LPT approaches to an Eulerian framework.In Sect.2, we briefly review the theoretical background.In Sect.3, we present the methodology, emphasising the ultraviolet regularisation, the Eulerian extension of LPT and the vorticity corrections.Section 4 discusses the numerical results.Finally, in Sect.5, we present our conclusions.

Theoretical background
The Universe is considered a closed Hamiltonian system where energy is conserved.The dark matter content of the Universe can be described by a distribution function f (r, u, t), with position r and velocity u, such that the probability of finding a dark matter particle in the phase-space volume drdu centred on r, u at time t is given by f (r, u, t) drdu (see, e.g., Mo et al. 2010).The total number of particles is then given by integrating over the whole phase-space volume: N = drdu f (r, u, t).Given the functional form of the distribution function, total and continuous changes to the phase-space of the system of particles can be expressed with product and chain derivative rules: , with g being the gravity-induced acceleration.Due to probability conservation, the generalised continuity equation in phase-space must be fulfilled (Liouville theorem): Consequently, to get the Vlasov or collisionless Boltzmann equation, the sum of terms multiplied by f must vanish: which is fulfilled by inserting Hamiltonian equations of motion, yielding: Building moments of the Boltzmann equation, one gets the following: From the zeroth moment, the continuity equation multiplying Boltzmann's equation with the mass of the particles m = m(u) 0 and integrating over du: d 3 v f mu = 0, with ρ being the density and u the velocity field.Hence, From the first moment, the Euler equation multiplying Boltzmann's equation with m = m(u) 1 and integrating over du: with w ≡ u − u being the random velocity.
The process of virialisation is expressed by a stress tensor T term: ww =−P1+T.The pressure term −P1 can be neglected for collisionless dark matter.From the combination of the Euler and continuity equations, we get: Introducing co-moving coordinates x: r = a(t)x, with the scale factor a encoding the expansion of the Universe and conformal time τ determined by dt = a dτ, we can rewrite the previous equation in terms of peculiar motions u ≡ a ẋ instead of proper velocity u = ṙ = ȧ(t)x + u, yielding: using the Poisson equation ∇ 2 Φ = 4πG ρδa 2 with density contrast δ ≡ ρ/ ρ − 1 and average density ρ ≡ ρ .
One usually assumes curl-free velocity fields, neglecting the stress tensor T. Based on this, one can perform a Eulerian perturbative expansion of both the continuity and the Euler equations around the density contrast δ(x, τ) = ∞ n=1 δ (n) (x, τ) and the divergence of the peculiar velocity field θ(x, τ) ≡ Bernardeau et al. 2002, and references therein).Alternatively, one can consider Eq. ( 4) and make an expansion in Lagrangian coordinates considering the total derivative d dτ u = ∂ ∂τ u + u • ∇u, yielding Lagrangian perturbation theory (LPT) solutions (Buchert & Ehlers 1993;Bouchet et al. 1995;Catelan 1995).(For nth order LPT, see Schmidt 2021).
The stress tensor is commonly represented by the linear elas- with viscosity parameters µ and η (see Bernardeau et al. 2002), which for irrotational fields (when the velocity is the gradient of a potential field) simplifies to the adhesion model ∇ • T ∼ µ ∇ 2 u with a single viscosity parameter µ (see Shandarin & Zeldovich 1989).However, the generation of vorticity (∇ × u) has been studied in simulations as an important component when going to the non-linear regime (see, e.g., Pueblas & Scoccimarro 2009;Jelic-Cizmek et al. 2018).
As an alternative to expensive N-body simulations, effective field theories have emerged in LSS, including a modelling of the stress tensor to compute summary statistics (see, e.g., Carrasco et al. 2012;Baumann et al. 2012;Pajer & Zaldarriaga 2013;Porto et al. 2014;Mercolli & Pajer 2014;Angulo et al. 2015;Baldauf et al. 2015Baldauf et al. , 2016;;Foreman et al. 2016).(See also other perturbative, McDonald 2011; Pietroni et al. 2012;Rampf 2012;Cusin et al. 2017;and non-perturbative approaches, Buchert & Domínguez 2005 to modelling the curl.) In this work, we propose modelling the viscosity induced by a stress tensor by iteratively applying LPT within a Eulerian framework with an ultraviolet regularisation given by the spherical collapse model.With this, the gravitational potentials become increasingly deeper and change their shape.In this way, the vorticity of the displacement field emerges naturally.

Method
The method is based on a small-scale or ultraviolet regularisation obtained with Augmented Lagrangian perturbation theory, which is applied iteratively within a Eulerian framework.Optionally, one may introduce vorticity corrections, as we discuss below.
We start by considering the Lagrangian q to Eulerian x coordinate single-step mapping through a displacement Ψ (with u = d dτ Ψ): x = q + Ψ(q).

Regularisation of the small scales
The displacement is obtained according to augmented LPT (ALPT), which separates the total displacement Ψ into a longrange Ψ L and a short-range component Ψ S (see Kitaura & Hess 2013): Ψ(q, z) = Ψ L (q, z) + Ψ S (q, z).
The long-range displacement field is obtained from the convolution of a Gaussian kernel with an LPT solution: Ψ L (q, z) = K(q, z, r s )•Ψ LPT (q, z).We restricted the present study to secondorder LPT: Ψ 2LPT (q, z) = −D(z)∇ q Φ (1) (q) + D (2) (z)∇ q Φ (2) (q), where D(z) is the growth factor (see, e.g., Heath 1977), and (Bouchet et al. 1995).The normalised potentials Φ i (q) are the solutions of the Poisson equations ∇ 2 q Φ (i) = δ (i) , where i = 1 is the linear primordial density A215, page 2 of 6 Kitaura, F.-S., et al.: A&A, 683, A215 (2024) field used as the initial conditions, and i = 2 is determined by The short-range displacement is given by Ψ S (q, z) = (1 − K(q, r s )) • Ψ SC (q, z), where the displacement Ψ SC (q, z) is derived within the spherical collapse (SC) approximation (see Bernardeau 1994), ψ SC (q, z) ≡ ∇ • Ψ SC (q, z), where ψ SC (q, z) is the solution to the Poisson like equation (see, e.g., Mohayaee et al. 2006;Neyrinck 2013 − 1 .The short-range force has two regimes, one given by the spherical collapse model in the mild-density regime and one given by the perfect collapse to a point (0 ) in the high-density regime when the EPT approach yields imaginary solutions.The combination of both long-and short-range forces obtained through a kernel with a transition scale parameter r s ) smears out the perfect collapse, partially emulating the result from quasi-virialisation.A second convolution with r s is applied to the primordial density field to determine the collapsing regions avoiding voids-in-clouds and is restricted to the high-density regime (Sheth & van de Weygaert 2004;Neyrinck 2016).Hence, the model has, so far, two free parameters, r s and r s , which are adjusted to improve the clustering towards small scales.

Eulerian augmented Lagrangian perturbation theory
In this section, we generalise the approach to arbitrary Eulerian coordinates, substituting q by x l for timestep l so that x l+1 = x l +Ψ(x l , z l+1 ) (see also Kitaura & Angulo 2012).The initial dark matter particle positions in LPT are regularly distributed on a mesh (q).In Eulerian-ALPT, however, the initial positions are given by the final positions of the previous step (x l ).
After each step, we make a mass assignment of the dark matter particles onto a mesh to get the density contrast at each snapshot δ(x l+1 , z l+1 ).This field is used instead of the Gaussian field to compute the displacement field for the next step.
Due to mass conservation, we can write ρ(x l+1 )dx l+1 = , where x l = x l−1 +Ψ(x l−1 , z l ) is not a homogeneous distribution after the first step.Nonetheless, we can perturbatively expand ρ(x l+1 ) ρ(x l )(1 + δ(x l+1 )), which enables us to apply LPT within a Eulerian framework.This implies that one should, in principle, use small step sizes between redshift calculations.Only the first (A)LPT step is conveniently selected to be exceptionally long, as the initial conditions can be considered homogenous for this case.
If we consider two subsequent steps, we can write x l+1 = x l−1 + Ψ(x l−1 , z l ) + Ψ(x l , z l+1 ).The total displacement to redshift z l+1 is then given by Ψ(q, z l+1 ) ≡ x l+1 − q = Ψ(q, z l−1 ) + Ψ(x l−1 , z l ) + Ψ(x l , z l+1 ).We assume curl-free fields in each step.Hence, ∇ x m × Ψ(x m , z m+1 ) = 0 for m = 1, . . ., n, with n being the total number of redshift snapshots.The shape of the gravitational potentials changes with redshift, Therefore, ∇ q ×Ψ(q, z l+1 ) 0 for l > 0, and z 1 is the first redshift snapshot.Hence, already after two steps, vorticity at Lagrangian coordinates emerges naturally.For this reason, it is wrong to assume that there are no curl sources after one step, even if one assumes curl-free fields initially.

Vorticity corrections
In particular, in each timestep we neglect the stress tensor term T and the divergence-free velocity component in Eq. ( 4).We can compensate for this vorticity leakage by adding in the subsequent step l + 1 a fraction of the vorticity generated in the previous step l, which is computed as the difference between the total displacement to redshift z l and the corresponding curl-free component.The displacement field can be decomposed according to the Helmholtz theorem into an irrotational or longitudinal (curl-free) and a solenoidal or transversal (divergence-free) component: To get the solenoidal component, we exploit the fact that the divergence of the curl of a vector field vanishes: where From the third step on (l > 1), we naturally get a transversal component that is not properly evolved in each step according to Eq. ( 4).A potential way of improving this could be achieved by introducing a transfer function T (q, z l ) in the following way (see also Tosone et al. 2020): We consider a reduced version of the transfer function to a single scalar value, a third free parameter α, in order to modulate the amplitude of the transversal component: We explore the cases with and without α boost contributions in the next section.

Numerical results
We consider in this study cubical boxes with a side length of 200 h −1 Mpc and 256 3 particles using Lambda-CDM cosmology with an initial power spectrum obtained from CLASS (Lesgourgues 2011;Blas et al. 2011) with the cosmological parameters {Ω M = 0.307, Ω Λ = 0.693, Ω b = 0.048, σ 8 = 0.823, w = −1, n s = 0.95} and a Hubble constant (H 0 = 100 h km s −1 Mpc −1 ) given by h = 0.68.We constructed two reference sets of 25 N-body calculations, each at z = 0 and z = 1, relying on FastPM (Feng et al. 2016) with the same volume and number of particles, the force resolution ∼0.39 h −1 Mpc, and 50 timesteps.Initial conditions were generated via second-order LPT at z = 99.Using the same Gaussian initial density field, we ran eALPT with (α > 0) and without (α = 0) vorticity corrections (vc), which we dubbed eALPT and eALPTvc, respectively.In particular, we studied the cases of eALPT with two, three, and five steps and eALPTvc with three, five, and six steps (see Figs. 2-4).
Figure 1 shows the Helmholtz decomposition for the FastPM displacement fields at z = 0. We find that the difference in power between the full and the longitudinal components is below 10% for k < 1 h Mpc −1 .Hence, focusing on the range of scales below k ∼ 1 h Mpc −1 , the alpha vorticity boost parameter cannot be larger than 10% in each timestep and should decrease further with an increasing number of timesteps.
The size of each eALPT (or eALPTvc) timestep is computed based on the difference between the linear growth factors from subsequent steps: For a given cosmology, we determine the corresponding redshift: ∆z(z l+1 , z l ) = D −1 (∆D(z l+1 , z l )), inverting the D = D(z) relation through an interpolation based on a previously computed table .We then use ∆z for D(∆z) and D (2) (∆z) to compute the displacement field to go to redshift z l+1 starting at redshift z l in Eulerian coordinates.
Because of numerical uncertainties arising from the massassignment scheme and from using significant steps between Eulerian coordinates at subsequent redshifts, the power at large scales of the resulting density field does not correspond to the one expected from linear theory.Hence, the step size has to be adjusted to get a precise final result.To this end, single-step ALPT calculations can be applied as a reference to determine excess or lack of power following Kitaura et al. (2021) using Gaussian convolutions of the density field to determine the variance and thus the resulting bias: with r S ∈ [30, 70] h −1 Mpc.We note that the exact choice of r S is not critical, and results within percentage accuracy are obtained for different values in that range.The timestep size has to be accordingly corrected: ).These corrections must be computed only once for a particular setup and can then be applied for different seed perturbations.
We find that the resulting density fields are very stable, relying on either cloud-in-cell or triangular shaped cloud (TSC) mass-assignment schemes (Hockney & Eastwood 1988) with and without tetrahedral tesselations (THT) (Abel et al. 2012).Finally, we used TSC with THT.
From our numerical tests, we draw following conclusions (see Figs. 2-4): 1.The ideal transition scale parameter r s between long-and short-range forces has a lower r s > r s−min and an upper limit of r s < r s−max .While r s−min has to be greater than zero to suppress shell-crossing in knots of the cosmic web, the accuracy in the tidal field tensor limits r s−max .In practice, we find 2 < r s < 10 h −1 Mpc. 2. The scale r s is restricted to the range 0 ≤ r s ≤ 2 dL, where dL is the cell side length.Neither larger smoothing scales nor additional convolutions have any impact on the results.Mpc side length with a mesh of 256 3 cells.Lower panels: Same as the upper panels but for the divergence of the irrotational (∇ q • Ψ ∇×Ψ=0 (q)) and the absolute value of the solenoidal component of the displacement field, Ψ ∇•Ψ=0 (q) (with power spectra within percentage agreement up to k 0.4 h Mpc −1 ).The resulting maps have been saturated to the same colour scale after averaging over five slices corresponding to one simulation.
3. The fraction of additional vorticity is constrained to the range α ∈ [1, 10] [%] and requires cross-correlations with the reference simulation, as a different combination of parameters can yield similar power spectra.In particular, the vorticity boost can partially compensate for the ultraviolet regularisation (compare numerical values of r S with and without α below).We find for the three-step case at z = 0 that the vorticity boost very moderately improves the results with α = 0 up to k 1.5 h Mpc −1 in terms of cross-correlation and are already worse in terms of power spectrum for k > 0.4 h Mpc −1 .It is for the various approximate gravity solvers compared to fast particlemesh based N-body calculations (FastPM).The Zel'dovich, 2LPT, and ALPT cases require only a one-step calculation.For the rest, the number of steps is indicated in parentheses.FastPM needs the first step for the ICs relying on 2LPT and then a minimal number of iterations compared to standard N-body simulations, which require iterations on the order of 10 3 .Also shown are the 1-sigma contours from 25 fields for FastPM and the eALPT-limiting cases with the lowest and highest cross-correlations, two (w/o vc) and six (with vc) steps, respectively.The eALPT (3s) and the eALPTvc (6s) cases yield power spectra with precision within a few per cent up to k 0.3 and 4 h Mpc −1 (the Nyquist frequency) as well as correlations of 95% and 98% at k 1 h Mpc −1 , respectively.
only by applying more timesteps that we can benefit from the vorticity boost.4. The three-step eALPT calculations neglecting the vorticity boost do not require any reference simulation and achieve percentage accuracy in the power spectrum up to k 0.3 h Mpc −1 and ∼95 % cross-correlation at k = 1 h Mpc −1 . 5. The three-step 2LPT calculation considerably improves the clustering towards small scales compared to the one-step 2LPT, but it lacks an ultraviolet regularisation, which significantly limits its accuracy towards high k-values, especially at low redshifts (compare results from Figs. 3 and 4).6.The z = 1 calculations show that the iterative application of LPT at subsequent times considerably improves the results compared to the one-step calculations.We also find that the ultraviolet regularisation becomes less critical.7. The optimal timestep choice consists of first making a long step with ALPT and significantly smaller subsequent steps.However, the timesteps need to be large enough to avoid numerical uncertainties arising from being too short (see Eq. ( 8)).The last timestep can be adjusted to ensure a vanishing bias, that is, b → 1 in Eq. ( 9) (in our calculation, we demand a precision of better than 0.5% in b 1).This greatly simplifies the calibration process for each set of free parameters.The parameters are chosen in order to obtain unbiased power spectra starting from the lowest modes and continuing towards the largest ones.
The concrete values for the eALPT runs are the following: For two steps, {r s = 5, r s = 2, α = 0}; for three steps, {r s = 6, r s = 2, α = 0} and {r s = 3.1, r s = 1.2, α = 7}; for five steps, {r s = 7, r s = 2, α = 0} and {r s = 3, r s = 2, α = 3.5}; and for six steps, {r s = 2.75, r s = 1.4,α = 3}.The numerical results demonstrate a higher accuracy when applying vorticity correctionsbut at the expense of depending on reference simulations, as the results need to be cross-checked with reference full-volume density fields.We also find that vorticity partially substitutes for the short-range component, as the transition scale parameter r S becomes smaller for α > 0. As expected from our previous discussion, the optimal α parameter remains below 10%.In fact, the first time we can introduce α is in a third timestep calculation, since vorticity emerges only after two timesteps, and we found α = 7% and even lower values for more timesteps.We checked the power spectra of the solenoidal component for both the reference FastPM runs and the eALPTvc(6s) calculations, and we found an excellent agreement, within percentage accuracy, up to k ∼ 0.5 h Mpc −1 .However, the A215, page 5 of 6 Kitaura, F.-S., et al.: A&A, 683, A215 (2024) deviations for that component can become larger than 10% towards higher frequencies.We also checked that we obtain an equivalent level of accuracy with different cosmological parameters.

Conclusions
Accurate calculations of the cosmic web dark matter distribution can be made very quickly with effective models by iteratively applying LPT within a Eulerian framework with an ultraviolet regularisation given by the spherical collapse model, which we dub eALPT.We showed that within three steps we can already find cross-correlations at the level of ∼95% at k = 1 h Mpc −1 and power spectra within a percentage accuracy up to k 0.3 h Mpc −1 as compared to N-body calculations.Moreover, this model is only about three to four times more expensive than setting up initial conditions for an N-body simulation.We also explored a variation of the eALPT method by adding a small fraction of the vorticity in each subsequent timestep, finding some improvement; however, it is at the expense of requiring calibration with a reference simulation.
We plan to use the method studied in this work to massproduce lightcone dark matter fields covering the entire redshift range up to z ∼ 4, populating them with bright galaxies, luminous red galaxies, emission line galaxies, quasars, Lymanα forests, and lensing maps in order to conduct a multi-tracer analysis from galaxy surveys.The developments presented in this work have applications ranging from setting up initial conditions, producing mock catalogues for galaxy surveys, and performing Bayesian inference analysis to providing simulations for emulators.Nevertheless, further investigation needs to be done to explore the accuracy within this framework in modelling the peculiar velocity field and to study different resolutions and redshifts.

Fig. 1 .
Fig. 1.Power spectra of the displacement fields from the Helmoltz decomposition for the FastPM simulations at z = 0 (mean and standard deviation) comparing the total component |Ψ| to the longitudinal one |Ψ ∇×Ψ=0 |.The 1-sigma contours have been computed from 25 FastPM simulations.

Fig. 2 .
Fig.2.Cosmic density and displacement field components as obtained with different methods.Upper panels: log-density slices at redshift z = 0 (log(1 + δ)) corresponding to different gravity solvers through the threedimensional simulation cubical box of 200 h −1 Mpc side length with a mesh of 256 3 cells.Lower panels: Same as the upper panels but for the divergence of the irrotational (∇ q • Ψ ∇×Ψ=0 (q)) and the absolute value of the solenoidal component of the displacement field, Ψ ∇•Ψ=0 (q) (with power spectra within percentage agreement up to k 0.4 h Mpc −1 ).The resulting maps have been saturated to the same colour scale after averaging over five slices corresponding to one simulation.