Issue 
A&A
Volume 683, March 2024



Article Number  A215  
Number of page(s)  6  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/202345876  
Published online  20 March 2024 
The cosmic web from perturbation theory
^{1}
Instituto de Astrofísica de Canarias, s/n, 38205 La Laguna, Tenerife, Spain
email: fkitaura@iac.es
^{2}
Departamento de Astrofísica, Universidad de La Laguna, 38206 La Laguna, Tenerife, Spain
^{3}
Department of Physics and Astronomy, Università degli Studi di Padova, Vicolo dell’Osservatorio 3, 35122 Padova, Italy
^{4}
INAF–Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, 35122 Padova, Italy
Received:
9
January
2023
Accepted:
2
February
2024
Context. Analysing the largescale structure (LSS) in the Universe with galaxy surveys demands accurate structure formation models. Such models should ideally be fast and have a clear theoretical framework in order to rapidly scan a variety of cosmological parameter spaces without requiring large training data sets.
Aims. This study aims to extend Lagrangian perturbation theory (LPT), including viscosity and vorticity, to reproduce the cosmic evolution from dark matter Nbody calculations at the field level.
Methods. We extend LPT to a Eulerian framework, which we dub eALPT. An ultraviolet regularisation through the spherical collapse model provided by Augmented LPT turns out to be crucial at low redshifts. This iterative method enables modelling of the stress tensor and introduces vorticity. The eALPT model has two free parameters apart from the choice of cosmology, redshift snapshots, cosmic volume, and the number of particles.
Results. We find that compared to Nbody solvers, the crosscorrelation of the dark matter distribution increases at k = 1 h Mpc^{−1} and z = 0 from ∼55% with the Zel’dovich approximation (∼70% with ALPT), to ∼95% with the threetimestep eALPT, and the power spectra show percentage accuracy up to k ≃ 0.3 h Mpc^{−1}.
Key words: methods: analytical / cosmology: theory / dark matter / largescale structure of Universe
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
The cosmic web is the largescale 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 threedimensional distribution of matter in the Universe through galaxy surveys such as DESI (Levi et al. 2013), Euclid (Amendola et al. 2016), JPAS (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 nonGaussianities (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, 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 (BalagueraAntolínez et al. 2019, 2020; Kitaura et al. 2022). These techniques can also accurately map the Lymanα forest (see, e.g., Sinigaglia et al. 2021, 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 particlemesh based Nbody solvers (see Tassev et al. 2013; Feng et al. 2016). Despite these developments, Nbody codes are costly when aiming to massproduce 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 nonlocal bias description (see the application of BAM to galaxy catalogues, BalagueraAntolínez et al. 2023).
Machine learning techniques applied to large training data sets (of thousands of Nbody 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.
2. 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, v, t), with position r and velocity v, such that the probability of finding a dark matter particle in the phasespace volume drdv centred on r, v at time t is given by f(r, v, t) drdv (see, e.g., Mo et al. 2010). The total number of particles is then given by integrating over the whole phasespace volume: N = ∫∫drdv f(r, v, t). Given the functional form of the distribution function, total and continuous changes to the phasespace of the system of particles can be expressed with product and chain derivative rules: , with g being the gravityinduced acceleration. Due to probability conservation, the generalised continuity equation in phasespace must be fulfilled (Liouville theorem): . Hence, . 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(v)^{0} and integrating over dv: , with ρ being the density and u the velocity field. Hence,
From the first moment, the Euler equation multiplying Boltzmann’s equation with m = m(v)^{1} and integrating over dv:
with w ≡ u − v being the random velocity.
The process of virialisation is expressed by a stress tensor T term: ⟨ww⟩=black−P1+T. The pressure term black−P1 can be neglected for collisionless dark matter. From the combination of the Euler and continuity equations, we get: . Introducing comoving 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 x ≡ aẋ instead of proper velocity u = ṙ = ȧ(t)x + v, yielding:
using the Poisson equation with density contrast and average density .
One usually assumes curlfree 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 and the divergence of the peculiar velocity field (see Bernardeau et al. 2002, and references therein). Alternatively, one can consider Eq. (4) and make an expansion in Lagrangian coordinates considering the total derivative , yielding Lagrangian perturbation theory (LPT) solutions (Buchert & Ehlers 1993; Bouchet et al. 1995; Catelan 1995). (black For nth order LPT, see Schmidt 2021).
The stress tensor is commonly represented by the linear elastic model 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}v with a single viscosity parameter μ′ (see Shandarin & Zeldovich 1989). However, the generation of vorticity (∇ × v) has been studied in simulations as an important component when going to the nonlinear regime (see, e.g., Pueblas & Scoccimarro 2009; JelicCizmek et al. 2018).
As an alternative to expensive Nbody 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. 2015, 2016; Foreman et al. 2016). (See also other perturbative, McDonald 2011; Pietroni et al. 2012; Rampf 2012; Cusin et al. 2017; and nonperturbative approaches, Buchert & Domínguez 2005 to modelling the curl.)
In this work, we propose modelling the viscosity induced by a stress tensor black 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.
3. Method
The method is based on a smallscale 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 singlestep mapping through a displacement Ψ (with ): x = q + Ψ(q).
3.1. 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 shortrange component Ψ_{S} (see Kitaura & Hess 2013): Ψ(q, z)=Ψ_{L}(q, z)+Ψ_{S}(q, z).
The longrange displacement field is obtained from the convolution of a Gaussian kernel with an LPT solution: Ψ_{L}(q, z)=𝒦(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 , where i = 1 is the linear primordial density field used as the initial conditions, and i = 2 is determined by δ^{(2)} = ∑_{i, j < i}(∂_{ii}Φ^{(1)}∂_{jj}Φ^{(1)}−(∂_{ij}Φ^{(1)})^{2}).
The shortrange displacement is given by Ψ_{S}(q, z)=(1−𝒦(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) .
The shortrange force has two regimes, one given by the spherical collapse model in the milddensity regime and one given by the perfect collapse to a point (0 = ∇_{q} ⋅ x = ∇_{q} ⋅ (q + Ψ)=3 + ∇_{q} ⋅ Ψ) in the highdensity regime when the EPT approach yields imaginary solutions. The combination of both long and shortrange forces obtained through a kernel with a transition scale parameter r_{s}) smears out the perfect collapse, partially emulating the result from quasivirialisation. A second convolution with is applied to the primordial density field to determine the collapsing regions avoiding voidsinclouds and is restricted to the highdensity regime (Sheth & van de Weygaert 2004; Neyrinck 2016). Hence, the model has, so far, two free parameters, r_{s} and , which are adjusted to improve the clustering towards small scales.
3.2. 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 EulerianALPT, 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} = ρ(x_{l})dx_{l}. Hence, , 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. black 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 curlfree fields in each step. Hence, ∇_{xm} × Ψ(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, for i = 1, 2. 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 curlfree fields initially.
3.3. Vorticity corrections
In particular, in each timestep we neglect the stress tensor term T and the divergencefree 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 curlfree component. The displacement field can be decomposed according to the Helmholtz theorem into an irrotational or longitudinal (curlfree) and a solenoidal or transversal (divergencefree) component: Ψ = Ψ_{∇ × Ψ = 0} + Ψ_{∇ ⋅ Ψ = 0}.
To get the solenoidal component, we exploit the fact that the divergence of the curl of a vector field vanishes:
where Ψ_{∇ × Ψ = 0} ≡ ∇∇^{−2}∇⋅Ψ(q, z_{l}).
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 𝒯(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:
black We explore the cases with and without α boost contributions in the next section.
4. Numerical results
We consider in this study cubical boxes with a side length of 200 h^{−1} Mpc and 256^{3} particles using LambdaCDM 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 black two reference sets of 25 Nbody calculations, black each at z = 0 black 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 secondorder 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 twoblack, three, and five steps and eALPTvc with three, five, and six steps black(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.
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 1sigma contours have been computed from 25 FastPM simulations. 
Fig. 2. Cosmic density and displacement field components as obtained with different methods. Upper panels: logdensity 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. 
Fig. 3. Power and cross power spectra at redshift z = 0 (averaging in Δkbins for fields A and B: ) for the various approximate gravity solvers compared to fast particlemesh based Nbody calculations (FastPM). The Zel’dovich, 2LPT, and ALPT cases require only a onestep 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 Nbody simulations, which require iterations on the order of 10^{3}. Also shown are the 1sigma contours from 25 fields for FastPM and the eALPTlimiting cases with the lowest and highest crosscorrelations, 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 black k ≃ 0.3 and 4 h Mpc^{−1}(the Nyquist frequency) as well as correlations of black 95% and 98% at k ≃ 1 h Mpc^{−1}, respectively. 
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, singlestep 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 . We note that the exact choice of is not critical, and results within percentage accuracy are obtained for different values in that range. The timestep size has to be accordingly corrected: Δz(z_{l + 1}, z_{l})=D^{−1}(ΔD(z_{l + 1}, z_{l})−(b(z_{l + 1}, z_{l})−1)). 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 cloudincell or triangular shaped cloud (TSC) massassignment schemes (Hockney & Eastwood 1988) with and without tetrahedral tesselations (THT) (Abel et al. 2012). Finally, we used TSC with THT.
blackFrom our numerical tests, we draw following conclusions (see Figs. 2–4):

The ideal transition scale parameter r_{s} between long and shortrange 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 shellcrossing 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.

The scale is restricted to the range 0 ≤ ≤ 2 dL, where dL is the cell side length. Neither larger smoothing scales nor additional convolutions have any impact on the results.

The fraction of additional vorticity is constrained to the range black α ∈ [1, 10] [%] and requires crosscorrelations 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 threestep 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 crosscorrelation and are already worse in terms of power spectrum for k > 0.4 h Mpc^{−1}. It is only by applying more timesteps that we can benefit from the vorticity boost.

The threestep 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 % crosscorrelation at k = 1 h Mpc^{−1}.

The threestep 2LPT calculation considerably improves the clustering towards small scales compared to the onestep 2LPT, but it lacks an ultraviolet regularisation, which significantly limits its accuracy towards high kvalues, especially at low redshifts (compare results from Figs. 3 and 4).

The z = 1 calculations show that the iterative application of LPT at subsequent times considerably improves the results compared to the onestep calculations. We also find that the ultraviolet regularisation becomes less critical.

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 blackfor the eALPT runs are the following: For two steps, {r_{s} = 5, = 2, α = 0}; for three steps, black{r_{s} = 6, = 2, α = 0} and {r_{s} = 3.1, = 1.2, α = 7}; for five steps, {r_{s} = 7, = 2, α = 0} and {r_{s} = 3, = 2, α = 3.5}; and for six steps, {r_{s} = 2.75, = 1.4, α = 3}. The numerical results demonstrate a higher accuracy when applying vorticity correctionsblack but at the expense of depending on reference simulations, as the results need to be crosschecked with reference fullvolume density fields. We also find that vorticity partially substitutes for the shortrange 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 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.
5. 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 crosscorrelations 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 Nbody calculations. Moreover, this model is only about three to four times more expensive than setting up initial conditions for an Nbody 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 multitracer 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.
Acknowledgments
FSK, FS, ABA, and GF acknowledge the Spanish Ministry of Economy and Competitiveness (MINECO) for financing the Big Data of the Cosmic Web project: PID2020120612GBI00 and the IAC for continuous support to the Cosmology with LSS probes project. FS thanks financial support from the University of Padova, ABA thanks support from the SEV20150548 grant, and GF thanks support from the IJC2020044343I grant.
References
 Abel, T., Hahn, O., & Kaehler, R. 2012, MNRAS, 427, 61 [Google Scholar]
 Amendola, L., Appleby, S., Avgoustidis, A., et al. 2016, Liv. Rev. Rel., 2018, 21 [Google Scholar]
 Angulo, R. E., & Hahn, O. 2022, Liv. Rev. Comput. Astrophys., 8, 1 [CrossRef] [Google Scholar]
 Angulo, R. E., Springel, V., White, S. D. M., et al. 2012, MNRAS, 426, 2046 [NASA ADS] [CrossRef] [Google Scholar]
 Angulo, R. E., Foreman, S., Schmittfull, M., & Senatore, L. 2015, JCAP, 2015, 039 [CrossRef] [Google Scholar]
 BalagueraAntolínez, A., Kitaura, F.S., PellejeroIbáñez, M., Zhao, C., & Abel, T. 2019, MNRAS, 483, L58 [CrossRef] [Google Scholar]
 BalagueraAntolínez, A., Kitaura, F.S., PellejeroIbáñez, M., et al. 2020, MNRAS, 491, 2565 [Google Scholar]
 BalagueraAntolínez, A., Kitaura, F.S., Alam, S., et al. 2023, A&A, 673, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baldauf, T., Mercolli, L., & Zaldarriaga, M. 2015, Phys. Rev. D, 92, 123007 [NASA ADS] [CrossRef] [Google Scholar]
 Baldauf, T., Schaan, E., & Zaldarriaga, M. 2016, JCAP, 2016, 017 [CrossRef] [Google Scholar]
 Baumann, D., Nicolis, A., Senatore, L., & Zaldarriaga, M. 2012, JCAP, 2012, 051 [Google Scholar]
 Benitez, N., Dupke, R., Moles, M., et al. 2014, arXiv eprints [arXiv:1403.5237] [Google Scholar]
 Bernardeau, F. 1994, ApJ, 427, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Bernardeau, F., Colombi, S., Gaztañaga, E., & Scoccimarro, R. 2002, Phys. Rep., 367, 1 [Google Scholar]
 Blas, D., Lesgourgues, J., & Tram, T. 2011, JCAP, 2011, 034 [Google Scholar]
 Bouchet, F. R., Colombi, S., Hivon, E., & Juszkiewicz, R. 1995, A&A, 296, 575 [NASA ADS] [Google Scholar]
 Buchert, T., & Domínguez, A. 2005, A&A, 438, 443 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Buchert, T., & Ehlers, J. 1993, MNRAS, 264, 375 [Google Scholar]
 Carrasco, J. J. M., Hertzberg, M. P., & Senatore, L. 2012, J. High Energy Phys., 2012, 82 [CrossRef] [Google Scholar]
 Catelan, P. 1995, MNRAS, 276, 115 [NASA ADS] [Google Scholar]
 Chuang, C.H., Kitaura, F.S., Prada, F., Zhao, C., & Yepes, G. 2015, MNRAS, 446, 2621 [NASA ADS] [CrossRef] [Google Scholar]
 Chuang, C.H., Yepes, G., Kitaura, F.S., et al. 2019, MNRAS, 487, 48 [NASA ADS] [CrossRef] [Google Scholar]
 Chudaykin, A., & Ivanov, M. M. 2019, JCAP, 2019, 034 [Google Scholar]
 Cusin, G., Tansella, V., & Durrer, R. 2017, Phys. Rev. D, 95 [CrossRef] [Google Scholar]
 Dai, B., & Seljak, U. 2021, Proc. Nat. Acad. Sci., 118, e2020324118 [NASA ADS] [CrossRef] [Google Scholar]
 DESI Collaboration (Aghamousa, A., et al.) 2016, arXiv eprints [arXiv:1611.00036] [Google Scholar]
 Feng, Y., Chu, M.Y., Seljak, U., & McDonald, P. 2016, MNRAS, 463, 2273 [NASA ADS] [CrossRef] [Google Scholar]
 Foreman, S., Perrier, H., & Senatore, L. 2016, JCAP, 2016, 027 [CrossRef] [Google Scholar]
 Fosalba, P., Crocce, M., Gaztañaga, E., & Castander, F. J. 2015, MNRAS, 448, 2987 [NASA ADS] [CrossRef] [Google Scholar]
 He, S., Li, Y., Feng, Y., et al. 2019, Proc. Nat. Acad. Sci., 116, 13825 [Google Scholar]
 Heath, D. J. 1977, MNRAS, 179, 351 [NASA ADS] [CrossRef] [Google Scholar]
 Hockney, R. W., & Eastwood, J. W. 1988, Computer Simulation Using Particles (Bristol: Hilger) [Google Scholar]
 Jamieson, D., Li, Y., Alves de Oliveira, R., et al. 2023, ApJ, 952, 145 [NASA ADS] [CrossRef] [Google Scholar]
 JelicCizmek, G., Lepori, F., Adamek, J., & Durrer, R. 2018, JCAP, 2018, 006 [CrossRef] [Google Scholar]
 Kitaura, F.S., & Angulo, R. E. 2012, MNRAS, 425, 2443 [Google Scholar]
 Kitaura, F.S., & Hess, S. 2013, MNRAS, 435, L78 [NASA ADS] [CrossRef] [Google Scholar]
 Kitaura, F.S., Yepes, G., & Prada, F. 2014, MNRAS, 439, L21 [NASA ADS] [CrossRef] [Google Scholar]
 Kitaura, F.S., GilMarín, H., Scóccola, C. G., et al. 2015, MNRAS, 450, 1836 [NASA ADS] [CrossRef] [Google Scholar]
 Kitaura, F.S., RodríguezTorres, S., Chuang, C.H., et al. 2016, MNRAS, 456, 4156 [NASA ADS] [CrossRef] [Google Scholar]
 Kitaura, F.S., Ata, M., RodríguezTorres, S. A., et al. 2021, MNRAS, 502, 3456 [NASA ADS] [CrossRef] [Google Scholar]
 Kitaura, F.S., BalagueraAntolínez, A., Sinigaglia, F., & PellejeroIbáñez, M. 2022, MNRAS, 512, 2245 [NASA ADS] [CrossRef] [Google Scholar]
 Lesgourgues, J. 2011, arXiv eprints [arXiv:1104.2932] [Google Scholar]
 Levi, M., Bebek, C., Beers, T., et al. 2013, arXiv eprints [arXiv:1308.0847] [Google Scholar]
 McDonald, P. 2011, JCAP, 2011, 032 [CrossRef] [Google Scholar]
 Meerburg, P. D., Green, D., Flauger, R., et al. 2019, BAAS, 51, 107 [NASA ADS] [Google Scholar]
 Mercolli, L., & Pajer, E. 2014, JCAP, 2014, 006 [CrossRef] [Google Scholar]
 Mo, H., van den Bosch, F. C., & White, S. 2010, Galaxy Formation and Evolution (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Mohayaee, R., Mathis, H., Colombi, S., & Silk, J. 2006, MNRAS, 365, 939 [CrossRef] [Google Scholar]
 Neyrinck, M. C. 2013, MNRAS, 428, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Neyrinck, M. C. 2016, MNRAS, 455, L11 [Google Scholar]
 Pajer, E., & Zaldarriaga, M. 2013, JCAP, 2013, 037 [CrossRef] [Google Scholar]
 Pietroni, M., Mangano, G., Saviano, N., & Viel, M. 2012, JCAP, 2012, 019 [CrossRef] [Google Scholar]
 Porto, R. A., Senatore, L., & Zaldarriaga, M. 2014, JCAP, 2014, 022 [CrossRef] [Google Scholar]
 Pueblas, S., & Scoccimarro, R. 2009, Phys. Rev. D, 80, 043504 [NASA ADS] [CrossRef] [Google Scholar]
 Rampf, C. 2012, JCAP, 2012, 004 [Google Scholar]
 Schmidt, F. 2021, JCAP, 2021, 033 [CrossRef] [Google Scholar]
 Shandarin, S. F., & Zeldovich, Y. B. 1989, Rev. Mod. Phys., 61, 185 [NASA ADS] [CrossRef] [Google Scholar]
 Sheth, R. K., & van de Weygaert, R. 2004, MNRAS, 350, 517 [NASA ADS] [CrossRef] [Google Scholar]
 Sinigaglia, F., Kitaura, F.S., BalagueraAntolínez, A., et al. 2021, ApJ, 921, 66 [NASA ADS] [CrossRef] [Google Scholar]
 Sinigaglia, F., Kitaura, F.S., BalagueraAntolínez, A., et al. 2022, ApJ, 927, 230 [NASA ADS] [CrossRef] [Google Scholar]
 Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv eprints [arXiv:1503.03757] [Google Scholar]
 Tassev, S., Zaldarriaga, M., & Eisenstein, D. J. 2013, JCAP, 6, 036 [CrossRef] [Google Scholar]
 Tosone, F., Neyrinck, M. C., Granett, B. R., Guzzo, L., & Vittorio, N. 2020, MNRAS, 498, 2663 [NASA ADS] [CrossRef] [Google Scholar]
 Zel’dovich, Y. B. 1970, A&A, 5, 84 [NASA ADS] [Google Scholar]
 Zhao, C., Chuang, C.H., Bautista, J., et al. 2021, MNRAS, 503, 1149 [CrossRef] [Google Scholar]
All Figures
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 1sigma contours have been computed from 25 FastPM simulations. 

In the text 
Fig. 2. Cosmic density and displacement field components as obtained with different methods. Upper panels: logdensity 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. 

In the text 
Fig. 3. Power and cross power spectra at redshift z = 0 (averaging in Δkbins for fields A and B: ) for the various approximate gravity solvers compared to fast particlemesh based Nbody calculations (FastPM). The Zel’dovich, 2LPT, and ALPT cases require only a onestep 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 Nbody simulations, which require iterations on the order of 10^{3}. Also shown are the 1sigma contours from 25 fields for FastPM and the eALPTlimiting cases with the lowest and highest crosscorrelations, 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 black k ≃ 0.3 and 4 h Mpc^{−1}(the Nyquist frequency) as well as correlations of black 95% and 98% at k ≃ 1 h Mpc^{−1}, respectively. 

In the text 
Fig. 4. Power and cross power spectra at redshift z = 1 analogously to Fig. 3. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.