CosmoMIA: cosmic web-based redshift space halo distribution

Modern galaxy surveys demand extensive survey volumes and resolutions surpassing current dark matter-only simulations' capabilities. To address this, many methods employ effective bias models on the dark matter field to approximate object counts on a grid. However, realistic catalogs necessitate specific coordinates and velocities for a comprehensive understanding of the Universe. In this research, we explore sub-grid modeling to create accurate catalogs, beginning with coarse grid number counts at resolutions of approximately 5.5 h -1 Mpc per side. These resolutions strike a balance between modeling nonlinear damping of baryon acoustic oscillations and facilitating large-volume simulations. Augmented Lagrangian Perturbation Theory (ALPT) is utilized to model the dark matter field and motions, replicating the clustering of a halo catalog derived from a massive simulation at z = 1.1. Our approach involves four key stages: Tracer Assignment: Allocating dark matter particles to tracers based on grid cell counts, generating additional particles to address discrepancies. Attractor Identification: Defining attractors based on particle cosmic web environments, acting as gravitational focal points. Tracer Collapse: Guiding tracers towards attractors, simulating structure collapse. Redshift Space Distortions: Introducing redshift space distortions to simulated catalogs using ALPT and a random dispersion term. Results demonstrate accurate reproduction of monopoles and quadrupoles up to wave numbers of approximately k = 0.6 h Mpc-1. This method holds significant promise for galaxy surveys like DESI, EUCLID, and LSST, enhancing our understanding of the cosmos across scales.


Introduction
Large Scale Structure (LSS) analyses provide information on the content and evolution of the Universe.Over the last 20 years, various experiments have pushed the boundaries of observations exponentially.The pioneering Two Degree Field Galaxy Redshift Survey (2dFGRS) captured around 390,000 galaxy redshifts within a sky area of 2 × 103 deg2 and a volume of 0.12 Gpc 3 /h 3 [1].Following suit, the Baryon Oscillation Spectroscopic Survey (BOSS) under the Sloan Digital Sky Survey (SDSS) probed 1.5 million galaxies sprawled over an area of ∼ 104 deg 2 [2,3], corresponding to a volume of around 2.5 Gpc 3 /h 3 .The extended BOSS (eBOSS) survey extended the depth of BOSS, which increased observed volume by ∼ 1.5 Gpc 3 /h 3 .The Dark Energy Spectroscopic Survey (DESI) is expected to observe over 30 million redshifts in 14 × 10 3 deg 2 [4][5][6] which amounts to a volume of ∼ 20 Gpc 3 /h 3 .Simultaneously, the Euclid satellite 1 is observing approximately 35 million galaxy redshifts and 1.5 billion shapes in order to combine galaxy clustering and weak lensing measurements [7].
Additionally, future-generation surveys are on the horizon, including 4MOST 2 (4-metre Multi-Object Spectroscopic Telescope) [8], HETDEX 3 [9] (Hobby-Eberly Telescope Dark Energy Experiment), PFS 4 (Subaru Prime Focus Spectrograph) [10] and the Roman Space Tele-scope 5 .These undertakings aim to decipher the properties of dark energy (DE) through meticulous measurements of the Baryon Acoustic Oscillations (BAO) scale and the universe's growth rate via Redshift Space Distortions (RSD).In essence, the forthcoming years are poised to witness a substantial amplification in the observable volume of the universe.
Uncertainties in the cosmological measurement are traditionally estimated using a sample of simulated universes that emulate the observed data.This requires, in particular, that these mocks have a volume at least as large as the observations, in order to properly estimate the statistical uncertainties in the data.However, current surveys such as DESI are already too big for an N -body simulation to be able to capture the total observed volume.The AbacusSummit [11] simulation suite contains boxes of only 8 Gpc 3 /h 3 which are not enough to contain the observed light cone.To encompass all tracers, including quasars, boxes with a volume of approximately 1000 Gpc 3 /h 3 would be necessary.In addition, in order to estimate accurate covariances for cosmological measurements, it is necessary to perform thousands of simulations.Moreover, the number of simulations which are required increases as the data vector size does, so the advent of multi-tracer and alternative analyses [12][13][14] is expected increase the number of mocks to be produced.
To counter the huge computational burden of precise N -body simulations, fast and approximate simulations are performed.These simulations sacrifice both mass resolution and gravitational evolution accuracy in order to decrease computing time.Over the years, many such methods have emerged.Some of these use variations of Lagrangian perturbation theory in order to estimate the evolved dark matter (DM) density field [15][16][17][18], whereas some others use particle mesh (PM) methods [19][20][21].Regardless of the method, if the volume requirements of the survey are large enough, the resolution of the mocks must decrease in order to maintain the simulations computationally tractable, which means that only large scales will be appropriately modelled by the simulations.In order to take full advantage of observations, especially of small-scale information, our mocks and models should be able to accurately capture them, which low-resolution simulations alone cannot do.
In this work we introduce the Cosmological Multiscale Infall Algorithm (CosmoMIA) 6 ; a model that corrects the positions of particles given by a low resolution approximate simulation in order to accurately emulate the subgrid clustering from a target N -body simulation, that is, the clustering on scales smaller than the spatial resolution of the approximate simulation.In section 2 we briefly introduce the data used to showcase the capabilities of our model as well as including a short description the approximate gravity solver used in this work, then in section 3 we fully describe our subgrid method.Section 4 introduces the metrics used to evaluate and optimize the model parameters.In section 5 we show our results and finally we discuss and conclude in section 6.

Data
In order to apply our model, we need two main ingredients: a reference and a fast approximate simulation.For the former, we use an N -body simulation to extract the relevant target clustering on the scales desired, plus it allows to minimize the effects of cosmic variance through the use of its initial conditions in the approximate solver.On the other hand, the approximate simulation is the basis upon which we want to improve, correcting their smallscale clustering.

N -body reference simulations
To introduce our model in this particular case, we use halo catalogs from the publicly available AbacusSummit suite of simulations [11] which are based on the Abacus code [22,23].The simulation suite consists of 97 different cosmologies around the central "base" Planck cosmology [24], each of which was run with a mass resolution of 2 × 10 9 M ⊙ /h (6912 3 particles in a (2 Gpc/h) 3 volume) in order to meet the requirements of the DESI [4] tracers.In this work, we use a subsample of one realization of the Abacus halo catalogs in the base cosmology, at a redshift of z = 1.1.This downsampling is done such that the most massive halos are used until the number density of halos in the box is 3×10 −3 (h/Mpc) 3 .This corresponds to 8×10 6 halos.
Furthermore, we use the corresponding initial conditions in order to run an approximate low-resolution simulation using the WebON code [25].The Abacus initial density fields, initially on a grid of 576 3 cells, were downsampled with a sharp Fourier space filter to a mesh of N low = 360 3 cells, which corresponds to a cell side size of ∆x ≈ 5.5 Mpc/h.In terms of mass resolution, our approximate simulation has ∼ 10 −4 × less resolution than the reference.
We must emphasize the multiple reasons behind conducting this study at the chosen resolution.Firstly, a minimum resolution of between 5 and 10 Mpc is necessary to adequately resolve the displacement and peculiar velocity fields for accurately modeling the BAO and the RSD [see, e.g., 26].Additionally, we aim at exploring the lowest possible resolution to showcase the effectiveness of our method.
Such low resolutions facilitate the efficient generation of large-volume mock datasets.In fact, we conducted a series of 200 large-volume light-cone dark matter simulations using the WebON code within cubical volumes of 1000 Gpc 3 /h 3 [25], precisely with the resolution explored in this study.

Approximate simulations
A very precise description of the evolution of the ensemble of simulation particles under gravity can in general be computed with N -body codes such as Abacus [22], Gadget3 [27] or PkdGrav3 [28], given that enough particles are used i.e. there is a fine enough mass resolution.On top of the high resolution requirements, current and future surveys impose large volume requirements that are practically intractable to an N -body simulation.Even when these requirements are relaxed, these simulations remain too expensive to be used to estimate cosmological measurement covariances and thus estimate uncertainties in cosmological parameters.In order to estimate these, cosmological analyses have used fast simulations.These are simulations that are relatively cheaper in terms of compute time, but sacrifice mass resolution and accuracy in the gravitational evolution, usually because they use an approximate gravity solver, which models the real displacement field analytically or with a reduced-resolution particle-mesh (PM) solution.In this study, we use Augmented Lagrangian Perturbation Theory (ALPT) [16], an analytical approximation to the nonlinear displacement field that improves upon the second order LPT (2LPT) by modelling the small-scale displacement with spherical collapse (SC).This gravity solver was used in the Patchy mocks [15] used by BOSS.
ALPT combines 2LPT for the long-range displacement and SC for small scales through a Gaussian kernel G(q, r s ) evaluated at the Lagrangian particle positions q with a width of r s = 6 Mpc/h (in this study).The ALPT displacement at redshift z is then given by: where " * " denotes convolution.The displacement field encodes all the necessary information for the application of our model.The final (Eulerian) DM particle positions are given by x = q + Ψ ALPT , whereas the Lagrangian coherent velocities are computed as and the Eulerian velocities can be computed from these through v coh (x, z) = v coh (q + Ψ ALPT (q, z), z).The growth rates f 1,2 are given by [29] f Finally, we rely on the cosmological environment of the particles, not only in the use of the local dark matter density, but also in the form of a cosmic web type, which can be computed within the δ-web prescription introduced in the companion papers [25,30,31].This new cosmic web classification scheme identifies cells based on the eigenvalues λ δ of the Hessian of the dark matter field Γ ij ≡ ∂ i ∂ j δ.Alternatively, the ϕ-web prescription introduced by [32] may also be used.This classification relies on the eigenvalues λ ϕ of the tidal field tensor T ij ≡ ∂ i ∂ j ϕ, i.e. the Hessian of the potential.Both cosmic-web classification schemes require a threshold λ th [33] and define regions as The present study, however, only uses the δ-web classification.The inclusion of a hierarchical cosmic web classification (see [31]) in our subgrid modelling is left for future work.Table 1 summarizes the input information that our model requires.
In this work, we assume that a perfect bias model is available and our approximate number counts will be given by the real Abacus halo catalog painted into a mesh of N low cells using Nearest Grid Point (NGP) mass assignment.Velocity kernel used to increase the power of the approximate coherent velocity field.Left: Fourier-space representation of the kernel for different values of the gain parameter, Γ -including our fiducial choice Γ = 10 -and a fixed scale parameter q = 0.6.Right: Effect of the velocity kernel on the monopole and quadrupole of the power spectrum of a sample with random subgrid assignment.We show the Abacus power spectrum as reference.

Velocity Kernel
To accurately fit the halo velocities from the Abacus simulation, correcting the peculiar velocity field becomes imperative when approximate gravity solvers are utilized, as they frequently struggle to capture the nonlinear regime of structure formation effectively.
In the nonlinear regime, particles begin to undergo virialization, and the velocity field exhibits an additional dispersed component alongside the coherent one.Handling the dispersed velocity component stochastically grants us a degree of freedom in adjusting the accuracy of the coherent component.The concept involves empowering the coherent velocity component with sufficient strength at smaller scales.
To achieve this, we convolve the velocity field with an isotropic kernel designed to preserve power on large scales and inject more in the small scales.The enhanced coherent velocity field is given by (2.4) In principle, the gain and scale parameters of the kernel (Γ, q respectively) can be free.However, we find that the scale parameter does not significantly affect the two-point clustering when varying in the range of q ∈ [0.2, 0.8] and thus fix it to q = 0.6.The gain parameter, Γ on the other hand, does significantly affect the clustering, specially the quadrupole, as can be seen in fig. 1.We find that a gain of Γ = 10 gives enough power on large scales for this particular application.These may vary from reference to reference.
Following this, at a later stage (see next section), we introduce a dispersed velocity component to systematically diminish this strength, thus enabling us to adjust the overall peculiar velocity on a global scale.

Method
The method comprises four stages.Initially, we compute the dark matter density field on a mesh using approximate gravity solvers (see section 2.2).We allocate dark matter particle positions to the tracers with a novel adaptive perturbative prescription (see section 3.1).
Additionally, we compute the coherent peculiar velocity field using a nonlinear enhancement prescription (see section 2.2.1).The final and most novel aspect of the method, as outlined in section 3.2, entails assigning halo positions and peculiar velocities exploiting the phase-space and cosmic web information.

Tracer Assignment
Given the tracer number-counts field, to obtain mocks appropriate for LSS clustering it is necessary to obtain a discrete particle distribution.In order to do this, the appropriate number of particles are often generated randomly around the center of each cell, which ignores the existence of any substructure below the scale of the grid.As an alternative, the DM particles themselves can be regarded as a proxy for the tracers, which preserves some information on the sub-grid physics.However, there are often not enough particles to satisfy the requirements of the target tracer field, which usually implies that a larger (higher-resolution) simulation is required.
Our approach is to use a hybrid method.Let N t be the number of tracers required in each cell and N DM the number of dark matter particles located in the same cell.Ideally, there are enough DM particles to assign to all tracers (N t ≤ N DM ).However, there may be cases in which this condition is not satisfied.If N t > N DM > 0, then the position of the (N DM + l)−th tracer is given by x N DM +l = x i≤N DM + G, that is, the new position is assumed to be the position (x i≤N DM ) of another particle in the cell plus a Gaussian-distributed random perturbation G ∼ N (µ, σ), where µ = 0, σ = 0.1∆x and ∆x is the cell size.Finally, for low resolution simulations, the case where N DM = 0 becomes more likely; in which case we sample the particle position randomly around the center of the cell, c.The new particle position will then be given by x N DM +l = c + 0.5∆x sign(U)(1 − |U|) [17], where U ∼ U(−1, 1) is a random vector sampled from a uniform distribution between -1 and 1.
As a baseline, we will compare our method with catalogs generated from the same number-count field but using only randomly distributed particles as described in the previous paragraph.We refer to this catalog as the random sample.
In addition, during this step we assign coherent velocities v ′ coh to the new particles via Cloud-In-Cell (CIC) interpolation from the enhanced eulerian velocity fields (see section 3.2).We have also tested assigning the cell's velocity (i.e., NGP mass assignment) but found that the velocity distribution in such case would not be centered around zero.Moreover, using the same CIC interpolation scheme we assign dark matter overdensity δ DM .Finally, we use the NGP mass assignment scheme to assign a cosmic web type to each particle.These will be necessary in order to identify which particles will behave as attractors in a given collapse iteration (see section 3.2).

Tracer Collapse
Figure 2 shows (top-left panel) that randomly sampling particles within the cell will cause a loss of power with respect to the reference of 10% at k ≈ 0.1 h/Mpc and 50% already by k ≈ 0.2 h/Mpc.In order to remedy this we perform a two-step particle collapse.In principle, the number of collapse steps is arbitrary, however we find that two steps are sufficient, preserve some computational efficiency, and are physically motivated.
This collapse consists on identifying pairs of particles, one of which will be the "central" whereas the other will be the "satellite".In the collapse, the satellite will be moved radially towards the central by a factor ϵ c .We do not collapse all particles towards all other particles, instead we use some criteria that allow us to decide whether a pair of particles is to be collapsed.These criteria correspond first to whether it the particle is a real DM particle from the simulation or has been sampled randomly in the cell; and second on their local cosmic web environment.Using these criteria we select the subset of particles that 1. are real DM and 2. do not reside in cosmic voids.These are denominated attractors, and reflect the fact that we rely more on the DM particle positions than the randomly sampled ones and that we expect gravity to cause particles to collapse towards peaks in the DM field.
The first collapse step is done among these attractors exclusively, i.e., both central and satellite populations are the attractors and the selection criteria for collapsing pairs is simply given by the closest minimum distance.This means that we collapse attractors to their closest attractor.The second step is done taking the central population to be the attractors and the satellite population to be the rest.In this case, we modify the selection criteria to also include information on the local DM density.In particular, we will collapse satellites either to the closest central or to the central located in the highest density environment.

Redshift Space Distortions
The collapse step is isotropic and completely independent on tracer velocities.In order to generate accurate mocks, we must provide accurate redshift-space clustering, which requires us to tune the tracer velocities.In section 2.2.1 we introduced the isotropic kernel that we employ to enhance the coherent flow that the approximate gravity solver provides.This kernel has two parameters that depend on the reference but need not be fine tuned, that is, we only require the quadrupole (of the redshift-space clustering) to have enough power on small scales such that it can later be fine-tuned.This subsection explains how this fine-tuning is performed.
Once the quadrupole of the redshift-space clustering has enough power (i.e. more than the reference) on all scales, we add a dispersive velocity component to the satellite particles in order to account for tracer virialisation.The final form of the velocity dispersion is where G ′ is a Standard-Gaussian-distributed random vector, δDM is the DM density of the central particle, clipped to the [0, ∞) range and γ is the free velocity dispersion parameter.

Summary of the model
Let us now summarize the collapse step, i: 1. Neighbour identification: we limit the radius r c,i of influence of each collapse, which allows us to fine-tune small scales and be more computationally efficient.These searches are efficiently implemented in Julia using the CellListMap.jlpackage [34].
2. Pair identification: we define which particles will be collapsed ("satellites") and towards which other particle ("central") they will move.We do so based on their cosmic web environment, their radial distance and whether they were gravitationally evolved during the simulation or were randomly generated within the cell.We stress that the set of centrals/satellites is not constant and different subsets of particles can be used in each collapse iteration as explained before.
3. Collapse: this steps effectively multiplies the radial distance between the satellites and the corresponding central s by a factor ϵ c,i , that is, for ϵ c,i < 1 the pair collapses.
A summary of the free parameters of our collapse model can be found on table 2.

Evaluation metrics: Clustering statistics
In this section, we evaluate the performance of the method using several summary statistics.We begin by examining the two-point correlation function in real and redshift space, both in Fourier and configuration space, including multipole expansions (see section 4.1).Following that, we present our results on the three-point statistics in Fourier space, specifically focusing on the bispectrum (see section 4.2).

Two-point clustering
The two-point functions contain all the information of a Gaussian field.They are widely used in cosmology in order to extract information from the large-scale structure of the Universe [35,36] and are the main observables of modern spectroscopic instruments such as DESI.
Our main objective is to provide two-point statistics that are accurate in scales smaller than the approximate simulation's grid would allow.Given a catalog, we use the Fourier-space two-point function (i.e. the power spectrum).
where δ D is a Dirac delta distribution, k and k ′ are wavevectors, δ(k) is the Fourier-space overdensity field and ⟨ • ⟩ denotes averaging.Through this work we estimate power spectra using the CosmoCorr.jl7and pypowspec8 codes, both based on the C powspec9 code.However, spectroscopic observations cannot disentangle the cosmological redshift from the redshift due to the peculiar velocity of the tracers. in order to compute the redshift space clustering, we add the contribution of the peculiar velocity to the position of each tracer in the plane-parallel approximation, assuming a line of sight along the Z (third) axis: where Z z denotes the redshift space position along the line of sight, v is the halo's peculiar velocity along the same direction, and H(z) is the Hubble parameter at redshift z.
The peculiar velocity contributions induce anisotropy in the power spectrum which can be easily detected in its multipoles, P ℓ (k).These are the projections of the 2-dimensional power spectrum P (k, µ) onto a Legendre polynomial basis L ℓ (µ), where µ is the cosine of the angle to the line of sight, that is Alternatively, in configuration space we compute the two-point correlation function ξ(s).In practice, we estimate it using the natural estimator [37], which allows us to compute very small scale clustering that is impractical to compute in Fourier space: Equivalently, the RR(s, µ) term is the number of such pairs in a random catalogue.Given that our tests are performed on periodic boxes, the RR factor is computed analytically.The multipoles of the correlation function are analogously defined as Through this paper we use the correlation function computation implementations in the CosmoCorr.jland pyfcfc10 , the latter of which is based on the C FCFC11 [18,38] code.

Bispectrum
While the two-point clustering is a good descriptor of the matter field, higher order statistics such as the bispectrum are being increasingly used in order to constrain cosmology [39][40][41] and break degeneracies present in the two-point analyses.Moreover, these higher order statistics are instrumental in constraining alternative cosmological models including modified gravity [42,43] and primordial non-Gaussianity [44][45][46].Current and next generation mocks should then be able to properly emulate three-point functions and their covariances.The bispectrum is defined as In this work we evaluate the bispectrum in a configuration with k 1 = 0.1 and k 2 = nk 1 with n = 2, 3. k 3 is given by the closure relation enforced by the Dirac distribution δ D .For this bispectrum projection, we show it as a function of the angle θ 12 between the wavevectors k 1 and k 2 , given by k 2 3 = k 2 2 sin 2 θ 12 + (k 2 cos θ 12 + k 1 ) 2 .In practice, we use the GPU-enabled bispectrum implementation available in the jax-powspec12 package, based on Pylians313 and a Julia wrapper to the C library libbispec.
Moreover, recent analyses [41,47,48], have used the Sugiyama estimator [49] of the bispectrum in order to constrain cosmology.This estimator computes the bispectra in different k 1 , k 2 configurations for different multipoles ℓ 1 , ℓ 2 , L. A common configuation is the diagonal k 1 = k 2 ≡ k configuration.[49] shows that the ℓ 1 = ℓ 2 = L = 0 and ℓ 1 = ℓ 2 = 2, L = 0 have a high signal to noise ratio and are thus important for cosmological analyses.In practice we use the Triumvirate package [50] to compute the bispectra.

Free Parameters & Optimization
Our model has just 5 free parameters, θ, that can be optimized by trial and error to a certain extent.Our approach however is to use a minimizer to compute steps towards a reasonablelooking minimum.To do so, we define the loss function as the weighted sum of the mean absolute errors, compared to the target statistics x i .Our total loss is then where L P is the MAE of the real-space power spectrum monopole, L P,ℓ is the redshift space power spectrum ℓ−multipole and L ξ,0 is the redshift space correlation function monopole.The weighting factors a i not only help to regularize the dynamic range of the different statistics, but modifying them during the training of the model will help to speed up convergence.In the first stages of training we focus on the power spectrum in redshift space, making sure that a 2 > a 3 > a 1 and a 4 = 0. Once the larger scales have been fixed, we introduce the correlation function term with a 4 = 50 in order to help the model better fit the smallest scales.In practice we use the Optim.jl14package.

Results
Our fitting procedure can be run for an arbitrary amount of time, however after a couple of hours the precision is already good enough compared to the reference.This usually means that the redshift-space power spectrum has reached 1 − 2% accuracy down to scales of 0.4 − 0.5 h/Mpc.Notice that given that we share the same initial conditions, we do not take cosmic variance into account for the fit.The resulting values of our fit are also shown in table 2. The values of the collapse fractions are interesting since the first collapse seems to be reducing the pairwise distances to about 70% of the original value, this will generate excess power at scales under the first cutoff radius.The subsequent "collapse" step actually separates particles from the chosen attractors, diluting power in the very small scales, under the second collapse radius.The combination of these two accurately emulates the small-scale clustering.We comfortably recover the 50% loss in power at large scales and inject the necessary power to scales significantly smaller than the minimum requirements for BAO and RSD analyses.
Figure 2 shows the Fourier-space two-point clustering of the reference Abacus halo catalog compared to the baseline random subgrid particle assignment and our collapsed model for the subgrid clustering.The leftmost column shows that the real-space clustering is within 1% of the target up to scales of k ≈ 0.6 h/Mpc.This is a substantial improvement over the usual particle distribution technique where at these scales the power is already negligible.However, our main target is the two-point redshift space clustering, which we upweight in the loss function.The monopole is consistent with the reference to 1% precision up to k ≈ 0.5 h/Mpc and 2% up to k ≈ 0.7 h/Mpc.The quadrupole shows also a significant improvement, as it is consistent within 5% up to scales of k ≈ 0.4 h/Mpc.
The late addition of the configuration space monopole to the loss, allows for fine-tuning the smallest scales in order to properly reproduce the redshift space monopole.fig. 3 shows that the term used in the loss (redshift space monopole) is accurate within 2% down to scales of s ≈ 5 Mpc/h and 5% down to s ≈ 1 Mpc/h.On the other hand, the real-space case does not fit as well, specially on these small scales.The feature at s ≈ 0.5 Mpc/h, likely due to halo exclusion, makes the fit to these scales considerably more difficult, and the model was not explicitly trained on this particular measurement.Notice how the redshift space mapping, largely mitigates these effects and the model is able to capture the small scales better.We expect that this feature will not be present in galaxy clustering, thus improving the overall fit to these scales.Nonetheless, we find a 5% agreement with the reference down to s ≈ 4 Mpc/h in the real space two-point function.Similarly the quadrupole is also not explicitly trained on thus we do not capture the feature at s ≈ 2 Mpc/h.However, we are able to recover the agreement with the reference even at s ≈ 5 Mpc/h.This is also a significant improvement over the baseline, which in this case does not contain a dispersion velocity term.Finally, fig. 4 shows that the random assignment of particles performed in the baseline, dampens the strength of the BAO peak significantly.This was already seen for the EZmocks [17], where the artificial strengthening of the BAO peak in the initial conditions was necessary.Evidently, our model not only allows to fix small-scale clustering but also corrects smaller effects on relevant larger scales such as the sound horizon.
In order to model the covariance correctly, it is necessary to properly model the threepoint clustering of the sample [51].In fig. 5 we show the real and redshift space for one large-scale k 1 = 0.1 h/Mpc, k 2 = 2k 1 and one smaller scale k 1 = 0.1 h/Mpc, k 2 = 3k 1 B(θ 12 ) configurations.Evidently, the baseline particle assignment does not generate accurate bispectra, showing discrepancies of around 50 and 100% for the large and small scale distributions respectively.Given that we are using ideal number counts, this also shows that properly fitting the global PDF of number counts with a bias model is not enough for an accurate bispectrum and that an accurate subgrid description is necessary.On the other hand, our collapsed sample shows a significantly better agreement with the reference, with per cent differences of 5 and 10% for the large and small scale bispectrum configurations in redshift space, while the agreement for the real-space clustering is slightly worse.This is not surprising given that we focused on fitting the redshift space measurements.Furthermore, we want to check whether our model can reproduce some of the common bispectrum projections used for cosmological measurements.In order to do this, we use the Sugiyama estimator [49] and choose two configurations with a high signal-to-noise ratio, thus the most relevant for cosmological parameter estimation.We use the ℓ 1 = ℓ 2 = 0 and ℓ 1 = ℓ 2 = 2, L = 0 configurations.For the former we show the real and redshift space and for the latter we show the redshift space only measurements in fig.6.We observe that the random subgrid assignment causes the bispectrum multipoles shown to have a divergence of 50% compared to the reference at scales as large as k = 0.2 h/Mpc.On the contrary, our collapse method can match the Abacus 3-point clustering within 5% down to scales of k ≈ 0.6 h/Mpc.This kind of agreement may open the door for more precise cosmological analyses provided that theoretical models can also reach these scales.

Discussion and Conclusions
The present work presents a subgrid model designed to emulate the small-scale clustering of an accurate (N -body) reference using an approximate, fast simulation with a mass resolution of 10 −4 times the original one.We have explored not only the technical details of the modelling of the small-scale clustering but have also identified theoretically motivated solutions to offset the limitations of low resolution gravity solvers.We compare our model with a baseline created using a common (random) particle assignment scheme within each cell.
We show that even under the assumption that a perfect bias model exists, such that the target number counts field is exactly the one obtained from the simulation, the random particle assignment is not enough to obtain two and three point statistics that properly emulate the reference on the relevant scales.Within our model, we employ a dark matterrandom hybrid particle position assignment within each low resolution simulation cell that potentially mitigates the dampening of the baryon acoustic peak that is observed in other methods such as EZmocks [17].In addition, we developed a two-step collapse and disperse velocity model with 5 free parameters that takes into account cosmic-web and environmental information and can automatically be tuned to reproduce the real and redshift space two point clustering of the reference down to scales of k ≈ 0.6 h/Mpc or s ≈ 1 Mpc/h with an accuracy of 1% in the monopole and of 5% in the quadrupole.This accuracy is a large improvement compared to previous approaches to approximate mock generation as seen in [52], where only the N -body based methods such as COLA yield comparable results.This agreement satisfies the requirements from current generation large spectroscopic instruments.
We employ a multi-stage training scheme that allows for faster optimization by avoiding the computation of sub-dominant terms in the objective function such as the small scale configuration space two-point clustering.This term will only be computed in later stages of the training, when large scales have already been fit by minimizing the power spectrum terms of the loss.
Despite not being trained to do so, our model is also capable of fitting various bispectra configurations B(θ 12 ) with accuracies of 5 to 10%, a factor of 10 improvement over the baseline random particle assignment.This hints at the capability of our model to properly model the two-point covariance matrix, as shown by [51].Moreover, we test different cosmologically relevant bispectrum projections using the Sugiyama estimator.We find a 5% agreement with the reference, which is also a factor of 10 improvement over the naive random assignment of particles.The combination of these shows that this novel technique can be applied to the generation of mocks for higher order statistics analyses, which will be a major focus of spectroscopic surveys in the coming years.
There are several avenues for enhancing our method further.While in this study we have employed ALPT at extremely low resolutions (four orders of magnitude less than the reference simulation), alternative approaches such as eALPT [53] operating at higher resolutions or machine learning-based solutions [54] offer promising alternatives.These methods could provide an improved baseline for determining the positions and velocities of dark matter particles.The particle collapse step could potentially be enhanced using a machine learning approach.However, in this work, we aimed to demonstrate how we can leverage the information contained in the phase-space distribution of particles from approximate gravity solvers through straightforward prescriptions without requiring large training data sets.
Our model can be regarded as a post-processing step that is completely independent of the method used to model the tracer field, thus various bias models such as [16,17,55,56] can equivalently be used.In particular, this allows us to overcome the data-volume limitations of machine learning (ML)-enabled models or models that are in general limited by (GPU) memory such as [57], given that a lower resolution mesh is required.In addition, this step can be applied to existing mocks such as the MD-Patchy mock suite used in BOSS in order to mitigate the small-scale discrepancies with respect to N -body (e.g., [58]).
Finally, we stress that a great advantage the model is its flexibility in the fact that modifications can easily be made to include more cosmic-web related information in the definition of attractors or even the number of collapse steps showing its potential paths for further improvement.This method will potentially become crucial in the analysis of current and upcoming galaxy surveys.

Figure 1 .
Figure1.Velocity kernel used to increase the power of the approximate coherent velocity field.Left: Fourier-space representation of the kernel for different values of the gain parameter, Γ -including our fiducial choice Γ = 10 -and a fixed scale parameter q = 0.6.Right: Effect of the velocity kernel on the monopole and quadrupole of the power spectrum of a sample with random subgrid assignment.We show the Abacus power spectrum as reference.

. 4 )
DD(s, µ) is the number of pairs separated by a distance s in the data catalogue, normalised by the total number of pairs N D (N D −1), where N D is the number of tracers in the catalogue.

Figure 2 .
Figure 2. Power spectra of the reference halo catalog, a randomly sampled population of halos and our collapsed sample.Left: Real space monopole, : Redshift space monopole.Right: Redshift space quadrupole.Top row: comparison of the Abacus reference to a random sample.Bottom: Comparison of Abacus against our collapsed sample.Monopoles agree to 1% (dark grey band) up to k ∼ 0.6h/Mpc.The quadrupole agrees within 5% (light grey band) on the same scale range.

Figure 3 .
Figure 3. Small-scale correlation function of the Abacus halo catalog, a randomly sampled population of halos and our collapsed sample.Left: Real space monopole, : Redshift space monopole.Right: Redshift space quadrupole.Top row: comparison of the Abacus reference to a random sample.Bottom: Comparison of Abacus against our collapsed sample.Monopoles agree to 1% (grey band) down to to s ∼ 3Mpc/h.The quadrupole agrees within 5% down to to s ∼ 10Mpc/h.The zerocrossing of the quadrupoles makes the ratio diverge.

Figure 4 .
Figure 4. Large-scale correlation function of the Abacus halo catalog, a randomly sampled population of halos and our collapsed sample.

Figure 5 .
Figure 5. Bispectrum of the reference Abacus catalog, a randomly sampled population of halos and our collapsed sample for two different triangle configurations.Left two columns: Real and redshift space curves for k 1 = 0.1 h/Mpc; k 2 = 2k 1 .Right columns: Real and redshift space curves for k 1 = 0.1 h/Mpc; k 2 = 3k 1 .Top row: comparison of the Abacus reference to a random sample.Bottom: Comparison of Abacus against our collapsed sample.Gray areas show 5% discrepancy.The k 2 = 2k 1 agrees with the reference within 5% while the k 2 = 3k 1 configurations agree within 10%.

Figure 6 .
Figure 6.Diagonal B ℓ1ℓ2L (k, k) of the bispectrum in different (ℓ 1 , ℓ 2 , L) configurations of the Sugiyama estimator as a function of k.We show two of the highest signal-to-noise configurations.Our collapsed catalog matches the reference within 5% up to k = 0.5 h/Mpc in both real and redshift space for both configurations shown.

Table 1 .
Summary of the input required for our method.The approximate quantities are expected to be generated by an approximate gravity solver.

Table 2 .
Summary of the collapse model parameters, their description and their final values for our test on Abacus halos at z = 1.1.Velocity dispersion: we modify the velocity of the satellite particles by adding a stochastic or dispersive term such that their final velocity is v sat