DESI Mock Challenge: Constructing DESI galaxy catalogues based on FastPM simulations

Together with larger spectroscopic surveys such as the Dark Energy Spectroscopic Instrument (DESI), the precision of large scale structure studies and thus the constraints on the cosmological parameters are rapidly improving. Therefore, one must build realistic simulations and robust covariance matrices. We build galaxy catalogues by applying a Halo Occupation Distribution (HOD) model upon the \textsc{FastPM} simulations, such that the resulting galaxy clustering reproduces high resolution $N$-body simulations. While the resolution and halo finder are different from the reference simulations, we reproduce the reference galaxy two-point clustering measurements -- monopole and quadrupole -- to a precision required by the DESI Year 1 Emission Line Galaxy sample down to non-linear scales, i.e. $k<0.5\,h\mathrm{Mpc}$ or $s>10\,\mathrm{Mpc}/h$. Furthermore, we compute covariance matrices based on the resulting \textsc{FastPM} galaxy clustering -- monopole and quadrupole. We study for the first time the effect of fitting on Fourier conjugate [e.g. power spectrum] on the covariance matrix of the Fourier counterpart [e.g. correlation function]. We estimate the uncertainties of the two parameters of a simple clustering model and observe a maximum variation of 20 per cent for the different covariance matrices. Nevertheless, for most studied scales the scatter is between two to ten per cent Consequently, using the current pipeline we can precisely reproduce the clustering of $N$-body simulations and the resulting covariance matrices provide robust uncertainty estimations against HOD fitting scenarios. We expect our methodology will be useful for the coming DESI data analyses and their extension for other studies.


INTRODUCTION
The study of Large Scale Structure of the Universe has significantly improved in the last two decades leading to Baryon Oscillation Spectroscopic Survey (BOSS; Alam et al. 2017) and extended-BOSS (eBOSS; Alam et al. 2021a) surveys.They have published the largest 3D map of over 2 millions galaxies and quasars (Alam et al. 2021a).This has allowed the measurement of cosmological parameters to a percent-level precision studying Baryonic Acoustic Oscillations (BAO) and Redshift Space Distortions (RSD).
Currently, the Dark Energy Spectroscopic Instrument (DESI; Levi et al. 2013;DESI Collaboration et al. 2022) is a five years long spectroscopic survey that will outperform previous surveys by a an order of magnitude (DESI Collaboration et al. 2016a), aiming to constrain the cosmological parameters with precision at a sub-percent level.
With its 5000 robotically controlled optical fibres (Silber et al. 2023;Miller et al. 2023;DESI Collaboration et al. 2016b), DESI will scan a third of the sky to map 40 millions galaxies (Lan et al. 2023) and quasars (Alexander et al. 2023).Only after the five-month Survey Validation (DESI Collaboration et al. 2023a), DESI has measured the spectra of more than one million galaxies leading to the recent Early Data Release (EDR; DESI Collaboration et al. 2023b).
The sub-percent precision measurements expected from ongoing and future surveys require careful analyses of the systematic effects.To this end, the DESI Mock Challenge was launched as a series of studies and projects to build and validate the methodology for the cosmological analysis.In particular, one must find the most robust way to estimate the uncertainty of the measurements (Chuang et al. 2023).To achieve this goal, one needs to create multiple realistic simulations of the large-scale structure, which is required to lower the noise on covariance matrix and to describe accurately the non-linear scales.
In this study, we investigate the possibility to tune FastPM catalogues to reproduce the clustering of SLICS reference with the final goal of estimating the covariance matrix.In contrast to the other fast methods, FastPM uses accelerated particle-mesh solvers to evolve the dark-matter field, that should provide a higher accuracy of the large scale structure.The additional accuracy provided by FastPM can be important given the unprecedented statistical power of the DESI survey.Therefore, the final FastPM covariance matrix together with other methods -i.e.BAM, EZmock, Jackknife (Zhang et al. 2023), analytical models (Xu et al. 2013;Wadekar & Scoccimarro 2020;Wadekar et al. 2020) -are compared to the SLICS reference covariance matrix, in a parallel DESI Mock Challenge paper (Chuang et al. 2023).
Fundamentally similar to full -body simulations, FastPM evolves the dark matter field into the cosmic web, the skeleton of the large scale structure in the Universe (e.g.Mo et al. 2010;Wechsler & Tinker 2018).After the dark matter haloes are selected, one must implement galaxy-halo connection models (Wechsler & Tinker 2018) to assign galaxies.There are more empirically inspired models such as the Halo Occupation Distribution (HOD; e.g.Benson et al. 2000;Seljak 2000;Peacock & Smith 2000;White et al. 2001;Berlind & Weinberg 2002;Cooray & Sheth 2002) and Sub-Halo Abundance Matching (SHAM; e.g.Kravtsov et al. 2004;Tasitsiomi et al. 2004;Vale & Ostriker 2004) and more physically inspired ones such as full hydro-dynamical simulations (e.g.Schaye et al. 2010Schaye et al. , 2015;;Dubois et al. 2014;McCarthy et al. 2017;Pillepich et al. 2018;Davé et al. 2019) or Semi Analytical Models (SAMs; e.g.Guo et al. 2011;Gonzalez-Perez et al. 2014).In this case, we adopt a HOD model as it is one the most efficient ways to create mock galaxy catalogues.
Comparisons between the matter power spectra of FastPM and full -body references have shown deviations by ≈ 0.5 per cent at  = 0.3 ℎ/Mpc (Feng et al. 2016) and by ≈ 20 per cent at  = 1 ℎ/Mpc (Grove et al. 2022).In addition, Feng et al. (2016) have observed that the FastPM halo power spectrum can deviate up to 5 per cent around  ≃ 0.5 ℎ/Mpc.Lastly, they have noticed that the FastPM less massive haloes are not as well matched to the reference as the higher mass haloes.As a consequence, the purpose of the current paper is to show that using a HOD model to assign galaxies on FastPM halo catalogues, one can overcome the differences between the FastPM and the full -body simulations and thus reproduce the reference SLICS galaxy clustering.We note that the vanilla HOD model used in this work may not be sufficient to capture the physics of galaxy formation, as shown in (Alam et al. 2023;Chaves-Montero et al. 2023).Nevertheless, given the aim of this article, we are only concerned with the ability to produce the covariance matrix for twopoint galaxy clustering using FastPM.We defer the study of other observables such as galaxy-galaxy lensing, count-in-cell statistics, and the Voronoi volume function (Paranjape & Alam 2020) to future work.
Furthermore, we compare the impact of different clustering statistics and examine the effects of various scales on the HOD fitting.Finally, we calculate covariance matrices for all the studied scenarios and perform a comparison to understand the influence of the HOD modelling on the parameter uncertainty estimation.
In Section 2, we present the SLICS and FastPM simulations.The methodology that we follow is detailed in Section 3. We describe our results on the HOD fitting performance and the covariance matrix comparison in Section 4. In the end, Section 5 concludes the article.
The dark matter haloes have been selected by applying a spherical over-density halo-finder (Harnois-Déraps et al. 2013).Their mass function follows precisely the Sheth et al. (2001) fitting function, as shown in Figure 2 of Harnois-Déraps et al. (2018).The redshift of the halo catalogues included in this study is  = 1.041.Lastly, given that some halo catalogues have been corrupted at the run time, we are limited to only 139 independent mocks.This study, together with the BAM (Balaguera-Antolínez et al. 2022) and the DESI covariance matrix comparison paper (Chuang et al. 2023) are focused on the DESI Emission Line Galaxies (ELGs) sample.Therefore, the SLICS haloes have been populated using the HMQ HOD model presented in Alam et al. (2020), capable of describing the ELG galaxies -consult (Alam et al. 2021b) for a comparison between different HOD models for ELGs.In addition, the parameters of the HOD model have been set to match the expected linear bias of the DESI ELG sample.The final product is a set of 139 galaxy catalogues that are used as reference in all the studies mentioned before.More details about the SLICS galaxy catalogues production can be found in the DESI covariance matrix comparison paper (Chuang et al. 2023).

Fast Particle-Mesh
Accelerated Particle-Mesh (PM) solvers -such as the FastPM software (Feng et al. 2016) -are able to produce accurate halo populations with respect to the full -body simulations.Thus, they are suitable to accurately simulate large volumes.
FastPM makes use of a pencil domain-decomposition Poisson solver and Fourier-space four-point differential kernel to compute the force.Additionally, the vanilla leap-frog scheme for the time integration is adjusted to account for the acceleration of velocity during a step, allowing for the accurate tracking of the linear growth of large-scale modes regardless of the number of time steps.
For the current analysis.we have run FastPM with two resolutions, resulting in one set of 778 Low Resolution boxes (LR; 12963 particles) and one set of 141 High Resolution (HR; 1536 3 particles) catalogues.Both sets output snapshots at the same redshift ( = 1.041), and have the same box side length ( box = 505 Mpc/ℎ) and cosmology as the SLICS simulations.The particle mass of the HR simulations is 2.86444 × 10 9  ⊙ /ℎ, while the one of LR is 4.77 × 10 9  ⊙ /ℎ.The resolution of the force mesh is boosted by a factor of  = 2 compared to the number of particles per side, for both LR and HR.Lastly, 40 linear steps have been used to evolve the density field from  = 0.05 to  = 0.96.
Due to the small number of SLICS galaxy realisations, for 123 runs of the FastPM (LR and HR likewise), we use the SLICS initial conditions.This plays an important role to reduce the effect of the cosmic variance in the clustering statistics and thus in the HOD fitting.SLICS initial conditions have been built on a 1536 3 regular grid using the Zel'dovich approximation (Zel'dovich 1970) to displace the particles from the grid positions.Lastly, the initial conditions have been downgraded to the LR by cutting in Fourier space the high frequency modes larger than the Nyquist frequency corresponding to the LR field.
The halos have been selected from the dark matter field with the Friends-of-Friends halo finder in nbodykit (Hand et al. 2018).During the galaxy assignment process -Section 3.3 -we only make use of halos with a minimum mass of 5.72 × 10 10  ⊙ /ℎ.
Finally, in Section 3, when we mention FastPM, we imply for simplicity both HR and LR.We only make the distinction in the results section, i.e.Section 4.

Two point correlation function
Mathematically, the two-point correlation function (2PCF) is a continuous function that can describe the clustering of galaxies.However, given the discrete nature of the galaxy distribution in the Universe, the 2PCF is measured using discrete estimators.In the case of cubic mocks, one can implement the natural estimator (Peebles & Hauser 1974): where  (, ) and (, ) are the data and the random pair counts, respectively, as functions of the radial distance and the cosine of the angle between s and the line-of-sight (3) In the previous equations,  ⊥ and  ∥ are the perpendicular (⊥) and parallel (∥) to the line-of-sight components of s , respectively.While the  term is evaluated directly on the data catalogue,  is calculated theoretically.

Power spectrum
From the mathematical point of view, the power spectrum (k) is the Fourier Transform of the 2PCF.However, the limited volume of a survey or a simulation creates mode coupling and makes the two clustering measurements not completely equivalent.Consequently, we compute the power spectrum multipoles  ℓ () starting from the density field in Fourier space (k), as follows (de Mattia et al. 2021): where the unit vector k represents the direction of k and η is the global line-of-sight unit vector, chosen as the  axis of the cubic simulations.Finally, given the finite galaxy number density ng , the monopole shot-noise is computed as: while for the higher-order multipoles, the shot-noise is zero.
In practice, we harness the versatility of POWSPEC 3 described in Zhao et al. (2021) through its python wrapper4 to calculate the power spectra and their multipoles starting from the galaxy catalogues.We estimate the density field on a grid of size 512 3 , by applying the Cloud-In-Cell (CIC; Sefusatti et al. 2016) particle assignment scheme on the catalogues of galaxies.Lastly, we exploit the grid interlacing technique (Sefusatti et al. 2016) to reduce the alias effects at small scales.

Bi-spectrum
The power spectrum and the 2PCF are two-point clustering statistics, but higher order statistics are necessary to characterise more precisely the galaxy distributions.In this study, we also look at the threepoint clustering statistics, namely the bi-spectrum (k 1 , k 2 , k 3 ), the Fourier pair of the three-point correlation function (e.g.Bernardeau et al. 2002): The three vectors k 1 , k 2 , k 3 are chosen to form a triangle whose two of the three sides are fixed ( 1 = 0.1 ± 0.05 and  2 = 0.2 ± 0.05), but the angle  12 between k 1 and k 2 is varied from 0 to .In practice, we compute the monopole of the bi-spectrum on a grid of size 512 3 .

FastPM HOD model
The galaxy population and its associated clustering covariance matrix can potentially be influenced by halo properties beyond just mass, as shown in Alam et al. (2023).Nonetheless, such effects are expected to be small for large volume surveys such as DESI and hence we plan to address them in future work.Additionally, the FastPM haloes are less accurate than the ones from a -body simulation, thus we do not expect that the final HOD model and parameters maintain the same physical interpretation.
As a consequence, we can adopt the simple five-parameter HOD model described in Zheng et al. (2005) to assign galaxies to the FastPM halo catalogues, as long as the resulting clustering and covariance matrix match the reference.Nevertheless, in future work one can study more complex and more adapted models for the studied ELG sample.
The current model assumes that each halo can host at most one central galaxy with a probability B (1) = ⟨ cen ⟩( h ) dependent on the halo mass  h , where B () denotes the Bernoulli distribution and: with erf the error function: log  min is the halo mass at which the probability to host a central galaxy is one half and  log  controls the steepness of the transition from a probability of one to zero.Lastly, the positions and velocities of the central galaxies are precisely the values of their parent haloes.
In contrast, the number of satellite galaxies  sat per halo is sampled from a Poisson distribution P ( sat |⟨ sat ⟩( h )) with the mean: where  0 is a minimum halo mass threshold below which haloes cannot host satellite galaxies;  0 and  1 indicate the halo mass at which one halo hosts on average one satellite galaxy, and  is the power-law index.Furthermore, the positions and velocities of the satellite galaxies follow the Navarro-Frenk-White (Navarro et al. 1996, NFW) density profile.
In the interest of adjusting the smaller scales and the quadrupole, we introduce a velocity dispersion factor ( disp ) for the velocity parallel (∥) to the line-of-sight (i.e o in the current case) of the satellite galaxies, in addition to the five HOD parameters: where  halo ∥ is the velocity parallel to the line-of-sight of the satellites' parent halo.Finally, the six free parameters are fitted so that the resulting FastPM clustering matches the SLICS one.

HOD fitting
We would like to draw the attention of the reader to Eq. ( 22): The covariance matrix used to compute the  2  , Eq.( 21).It is not used for fitting.
Table 1.A summary of some of the most important and possibly confusing notations and their meaning.
MultiNest is a sampler based on Bayes' theorem that provides the maximum likelihood (best-fitting) parameters, as well as the posterior probability distribution of parameters alongside the Bayesian evidence.Bayes' theorem combines prior knowledge about the Θ parameters of a model  with information from the data  to calculate the posterior probability density of the Θ parameters: where (Θ|) is the prior distribution of Θ of the model , (|Θ, ) is the likelihood, and (|) is a normalizing factor called Bayesian evidence.
The uniform prior distributions that we impose on all six parameters are shown in Table 2. Furthermore, we approximate the likelihood by a multivariate Gaussian: with the chi-squared: where v is the difference between the data and model vectors v =  data −  model (Θ), and C is the covariance matrix.
The purpose of a covariance matrix C is to estimate the noise in the data, in the context of a noise-free model.Nevertheless, the peculiarity of this study is that both the model ( model (Θ), FastPM) and the data ( data , SLICS) are affected by noise.Due to the small volume of the SLICS and FastPM boxes, the cosmic variance component of the noise would be larger than the expected precision of Table 2.The limits of the uniform prior distributions included in the HOD fitting.Note that  0 from Eq. ( 10) is  0 ≡  ×  min . ⊙ denotes the solar mass.
ongoing surveys such as DESI.Therefore, the simulations have been run with matching initial conditions.In this case, the relevant noise component is no longer the cosmic variance but rather the accumulated noise due to gravitational evolution while starting with exactly the same initial conditions.Thus, the mock covariance estimated by SLICS or FastPM substantially over-estimates the error for our fittings, see Figure 2. Mathematically, since one computes the difference vector v =  data −  model (Θ) to estimate the  2 , one needs to estimate the covariance matrix of v, i.e.C(v): where C × represents the cross covariance matrix (Gubner 2006) between the data and the model vectors.In the more common case of independent data and model vectors and noise-free model, one obtains C(v) = C( data ).Nonetheless, in our case, the FastPM clustering model has positively correlated noise with the SLICS clustering -given the matching initial conditions -hence C( model ) ≠ 0 and C × ≠ 0.
Consequently, in order to more appropriately estimate the noisei.e. an estimation of C(v) -we perform a two-step HOD fitting as schematically shown in Figure 1: (i) we fit the monopole and quadrupole of the 2PCF [ 0 ,  2 ] and the galaxy number density  g using a diagonal covariance matrix (Σ diag ) and thus obtain an initial-guess (IG) best-fitting FastPM galaxy catalogues (IG-FastPM), see Section 3.3.1; (ii) we compute the differences [Δ 0 , Δ 2 ] between the clustering (monopole, quadrupole) of the IG best-fitting FastPM and the SLICS galaxy catalogues; we use these differences to calculate a new covariance matrix (Σ diff ) with which we perform again the fitting, see Section 3.3.2.
In both cases, we use 20 FastPM (F) and 20 SLICS (S) halo boxes ( fit mocks = 20) -sharing the same initial conditions -for the purpose of additionally decreasing the noise.Nonetheless, the average nS g is computed using 139 realisations, while the average nF g is calculated using the 20 realisations included in the HOD fitting.There are three main reasons behind this discrepancy: first, it quickly becomes expensive to apply galaxies using HOD to more than 20 FastPM simulations; second, the number of SLICS reference simulations has to be the same as for FastPM, so that the cosmic variance is reduced in the clustering by the shared initial conditions; third, the noise in the galaxy number density is not reduced by the shared initial conditions, thus one needs more realisations to estimate a (practically) noiseless SLICS reference galaxy number density.The galaxy number density is an important constraint as it governs the shot-noise which has a significant role in the covariance matrix.

MultiNest fitting
Final FastPM galaxy boxes, The two-step HOD fitting process that is detailed in Section 3.3.

The First Step
Initially, we perform the HOD fitting on the monopole and the quadrupole of the 2PCF, together with the galaxy number density.Hence, the data vector  data is formed by concatenating their respective averages for the SLICS (S) mocks: Considering that the computing time of clustering measurements scales with the maximum separation, we need a large enough upperlimit to constrain relevant parameters, but small enough to keep a reasonable execution time for model evaluation during the HOD fitting.Additionally, since we are interested in capturing the nonlinear effects, the lower-limit is set to 0. Consequently, the monopole and the quadrupole of the 2PCF are evaluated for  ∈ [0, 50] Mpc/ℎ, with a bin size of 5 Mpc/ℎ.Thus,  is an array containing 10 elements ( 1 , . . .,  10 ).
As previously argued, in the first step, there is no appropriate noise estimation.Therefore, we can use an approximate covariance matrix that enables us to proceed to the second step and calculate a more suitable one.In this regard, we create a diagonal covariance matrix: where the first 20 elements are defined as follows: This selection of the diagonal covariance matrix is based on an examination of the  2  SLICS () values, where  SLICS () represents the standard deviation of the SLICS 2PCF.Notably, the highest value is approximately three; hence, we initially approximate all values as three for simplicity.The last element   g is computed as the standard deviation of 139 SLICS galaxy number densities, divided by √ 139, so that it estimates the uncertainty corresponding to the average of 139 realisations.The strong constraint on the   improves the fitting time, as HODOR initially evaluates the goodness-of-fit based only on the nF g and nS g , and does not compute the clustering if nF g is 10 away from the reference.Additionally, the lack of covariance terms in the covariance matrix should, as well, decrease the convergence time.
Finally, we apply the best-fitting HOD model to all  cov mocks = 123 FastPM halo boxes that share the initial conditions with the SLICS mocks to obtain the IG-FastPM.

The Second Step
To examine the influence of smaller scales on the HOD fitting, we compute the following for both SLICS and FastPM: (i) the power spectrum for  ∈ [0.02,  max ] ℎ/Mpc, with a bin size of 0.02 ℎ/Mpc, (ii) the 2PCF for  ∈ [ min , 50] Mpc/ℎ, with a bin size of 5 Mpc/ℎ, where the values of  max and  min are presented in Table 3.Consequently, we create the data and model vectors as follows: In order to estimate the noise in the context of shared initial conditions between SLICS and FastPM, we use the  cov mocks galaxy boxes of both SLICS and IG-FastPM, along with their corresponding clustering measurements (power spectrum or 2PCF).Furthermore, we introduce Δ  ℓ,IG =  F ℓ,IG () −  S ℓ () and Δ  ℓ,IG =  F ℓ,IG () −  S ℓ () as well as the generic vector Δ IG () = [Δ 0,IG , Δ 2,IG ] to express the difference between the SLICS and the IG-FastPM galaxy clustering that share the initial conditions.Here, the variable  represents either  or .
Taking advantage of the previous definitions, we further define a matrix M with the following elements: where Δ IG  denotes the vector corresponding to the −th (SLICS, IG-FastPM) pair, ΔIG represents the mean vector over all (SLICS, IG-FastPM) pairs and [ min ,  max ] defines the interval of points involved in the fitting, see Table 3. Starting from this matrix and its transpose, we calculate the sample covariance matrix C  as follows: Lastly, we calculate the  ′  g as the standard deviation of 139 SLICS galaxy number densities, divided by √︃  fit mocks -so that it estimates the uncertainty corresponding to the average of  fit mocks realisations -and we attach it to the C  to obtain the final covariance matrix used in the HOD fitting: Note that while the error estimate for the clustering is based on the difference in clustering due to matched initial condition, the error of the number density is directly computed from the SLICS realisations, as we aim to constrain the absolute number density, which has strong effect on the final clustering covariance.

Goodness-of-fit
In this section, we define a reduced  2 - 2  -that expresses the goodness-of-fit for the average of  fit mocks FastPM galaxy clustering realisations with respect the SLICS reference, i.e. the  g is not included: where  denotes the difference between FastPM and SLICS clustering -monopole and quadrupole -and  =  bins −  params , with (i)  params = 6 -the number of free parameters; (ii)  bins = 2× ℓ bins -the length of the Δ IG () vector, see Table 3.
The Σ −1  is the unbiased estimate of the inverse covariance matrix (Hartlap et al. 2007): where C  is defined in Eq. ( 19).Sellentin & Heavens (2016); Percival et al. (2022) have shown that this correction may not be the optimal choice for accurately determining the uncertainty of the parameters.However, since our main focus is on obtaining the best-fitting clustering and assessing its goodness-of-fit, it remains a reasonable correction.
Finally, as we fit the average of  fit mocks realisations, we must scale the covariance matrix C  by a factor of 1/ fit mocks .As a consequence, the  fit mocks factor appears in Eq.( 21).

Covariance matrix comparison
Given that the main goal is to have a robust estimation of the uncertainty on the cosmological parameters, we want to compare the constraining power of the covariance matrices.To this end, we fit the 123 individual SLICS clustering (monopole and quadrupole) with the following models: and where Pℓ 123,SLICS () and ξℓ 123,SLICS () are averages of the 123 realisations and  ℓ denotes the two free parameters.
Moreover, the covariance matrices are computed similarly to the Eq. ( 22), but using 778 LR FastPM realisations.The fitting is performed using PyMultiNest, for different fitting ranges ( ∈ [0.02, K] ℎ/Mpc and  ∈ [S, 200] Mpc/ℎ, see Table 4) for the purpose of comparing the effect of the covariance matrices at different scales.The largest fitting intervals are chosen so that they cover the nominal scales included in the BAO and RSD analyses, i.e.K ≈ 0.2 ℎ/Mpc and S ≈ 20 Mpc/ℎ (e.g.Tamone et al. 2020;de Mattia et al. 2021).Finally, the shown values are the average ( ℓ ) and standard deviation (  ℓ ) of the marginalised posterior ( ℓ ) and covariance (R [ 0 ,  2 ]) of the posterior distribution of  0 and  2 , ( 0 ,  2 ).By construction, the values of  ℓ should be one.
The main reason why we perform such a simplified test is to avoid the systematic errors that can arise due to the modelling.Consequently, the comparison between the quoted   ℓ and R [ 0 ,  2 ] should be directly related to the differences in FastPM covariance matrices.We, nevertheless, reckon that these comparisons do not show how the errors on the parameters of a realistic BAO/RSD model would behave.

RESULTS
One of the challenges of HOD fitting is addressing the high precision imposed by large volume surveys such as DESI because it requires prohibitively many large volume simulations.Figure 2 illustrates this issue as a comparison between  20,SLICS 6 the noise corresponding to the average of  fit mocks = 20 SLICS clustering realisations and the expected DESI Y5 7 and Y1 8 errors of the ELG sample.It is obvious that  fit mocks SLICS realisations do not reach the required precision 9 .In order to overcome this issue, we employ the novel matched 6 This would be the noise level in a hypothetical case where SLICS and FastPM would not share the initial conditions. 7The DESI Year 5 error is estimated by rescaling  20,SLICS to match the Y5 ELG sample volume, which is assumed to be 24 Gpc 3 ℎ −3 . 8The DESI Year 1 error is estimated by rescaling  20,SLICS to match the Y1 ELG sample volume, which is assumed to be one third of the Y5 volume. 9A simple calculation reveals that one would need 192 SLICS realisations to meet the DESI Y5 precision requirements.initial conditions simulations (SLICS and FastPM).In this case, the effect of the cosmic variance on the clustering difference is mostly removed.Therefore, as discussed in Section 3.3, the relevant error estimate is given by the covariance matrix of the clustering difference between the two simulations.Given the fact that we use  fit mocks pairs to perform the HOD fitting, the covariance matrix must be rescaled by  fit mocks .The square root of the diagonal of the resulting covariance matrix is illustrated with an dashed orange line in Figure 2. One can observe that the matched initial conditions significantly reduce the noise to values below  20,SLICS .
Furthermore, we would like to highlight that the precision depicted by the dashed orange line is either better than or equal to DESI Y1 precision up to  ≈ 0.25 ℎ/Mpc.Consequently, the results presented in this paper are precise enough with respect to the requirements of further DESI Y1 analyses.Nonetheless, it might be necessary to readdress this study for the full DESI sample, to account for even lower noise levels.For this, one could use the 1800 AbacusSummit (Maksimova et al. 2021) -body 0.5 Gpc/ℎ cubic boxes.
In addition, Figure 2 illustrates the comparison between  DIFF and the square root of the diagonal elements of Σ diff .In this context,  DIFF represents the standard deviation of the differences between the best-fitting FastPM (obtained from the second HOD fitting step) and SLICS clustering, further divided by

√︃
fit mocks .Ideally, an iterative HOD fitting process should be performed to ensure a robust Σ diff , but the close agreement between  DIFF and the diagonal elements of Σ diff suggests that Σ diff has approximately converged after a single iteration.A more detailed argument in support of the convergence of Σ diff is presented in Section A.
As pointed out in Section 3.3, it is important that the FastPM galaxy catalogues reproduce the SLICS shot-noise.Examining the FastPM galaxy number densities of all HOD fitting cases, we observed that the largest deviation, | nS g − nF g |/ ′  g , is approximately 0.5 , but most values are below 0.2 .This strongly supports that the galaxy number density is well constrained and that it is safe to define a  2  without including  g -see Eq.( 21).Furthermore, the values of the  2  are subject to uncertainties due to the finite number of realisations used to estimate the covariance matrix and the limited number of HOD realisations per halo catalogue.The most significant uncertainty, ≈ 27 per cent, arises from the limited number of HOD realisations.The remaining values are below 20 per cent, see Section B for more details.The  2  is simply used as a metric to evaluate the goodness-of-fit.For this reason it is important to consider that it is affected by a large uncertainty when comparing its magnitude to the expected value of one.
The primary focus of this paper is to investigate the limits of the FastPM capabilities to model the non-linear scales captured by -body simulations.Furthermore, we study the effect of fitting to successively more non-linear scales and either Fourier or configuration space statistics on the FastPM covariance matrix.

Power spectrum fitting
Figure 3 shows the results of the HOD fitting performed on the power spectrum for three different  intervals, defined in Table 3.The second, third and fifth rows display the difference in the clustering scaled by the difference error.We remind the reader that this error is smaller than the expected one for the given volume, due to the matched initial conditions between the two simulations, see Figure 2.
The best fitting monopoles and quadrupoles are within ±1 for most scales.Moreover, the results for the HR FastPM -presented with dashed line -are only marginally better than the ones for LR FastPM.Given the modest difference between the performances of the two resolutions, we believe that the LR FastPM is precise enough to describe the two-point clustering to non-linear scales for a DESI ELG-like galaxy sample, within the estimated DESI Y1 error bars, see Figure 2.
Considering that we only fit the first two even multipoles, there is no guarantee that the third one would match the reference.Nevertheless, the fifth row of Figure 3 illustrates that fitting the monopole and quadrupole to smaller scales improves the agreement of the hexadecapole.For instance, fitting on the Large interval pushes the ℓ = 4 multipole within ±2 for  < 0.4 ℎ/Mpc, whereas for Medium and Small intervals, the hexadecapole is placed within ±2 only for  < 0.3 ℎ/Mpc or  < 0.2 ℎ/Mpc, respectively.
Due to the fact that the power spectrum is affected by the window function, it is not obvious that a good matching in Fourier space translates as a good matching in Configuration space.Thus, we compute and display the corresponding 2PCF in the right-hand side of Figure 3.Most monopoles and quadrupoles agree within ±2 with SLICS for separations larger than 20 Mpc/ℎ.This suggests that it is possible to obtain a reasonable 2PCF above a certain minimum separation, even when performing HOD fitting on the power spectrum.However, fitting on the Medium and Large intervals, the 2 matching goes down to a separation of 10 Mpc/ℎ.
In contrast, for separations smaller than 5 Mpc/ℎ, the non-linear effects become dominant, making it difficult to replicate the velocity field.This is why increasing the fitting range up to  max = 0.5 can improve the monopole but not the quadrupole.Lastly, the 2PCF hexadecapole exhibits a bias of over 3 for  < 50 Mpc/ℎ in all six cases.
After a more qualitative description of the results, we present the  2  values in the upper panels of Figure 4. Generally, the HR FastPM produces lower  2  values than the LR, as expected from Figure 3.However,  2  [(0.02, max )] ≃ 1, which reiterates that by fitting the monopole and quadrupole of the power spectrum up to the three  max values, one can achieve a good match with the SLICS reference, within the DESI Y1 precision even with LR.In addition,  2  [ (20, 50)] ≃ 2 for the small fitting interval of the LR power spectrum, reinforcing the fact that one can get a reasonable 2PCF above a certain minimum separtion threshold when the fitting is performed on the power spectrum.
Additionally, we can observe the behaviour of  2  when it is estimated on different intervals than those used for the fitting.When the fitting is performed on the Large interval, the  2  ≃ 1 for all smaller intervals, regardless of the resolution.However, fitting on the Medium interval shows that the difference between HR and LR becomes more significant for  > 0.4 ℎ/Mpc (see also Figure 3): the  2  ≃ 2 for LR, while for HR, it is close to one.These findings imply that fitting up to  ≤ 0.4 ℎ/Mpc is satisfactory for HR FastPM, whereas smaller scales play a more significant role in LR.
Furthermore, fitting on the Small interval shows that although  2  [(0.02,0.3)] ≃ 1, it is much larger for  > 0.3 ℎ/Mpc, indicating strong clustering divergence beyond that value (see Figure 3).Therefore, both LR and HR benefit from considering the clustering information contained in smaller scales  > 0.3 ℎ/Mpc.

2PCF fitting
When the HOD fitting is performed on the power spectrum, the minimum 2PCF  2  is  2  [ (10, 50)] ≈ 2. While this translates to a 2 agreement down to the separation of 10 Mpc/ℎ between FastPM and SLICS 2PCF, we test whether fitting directly the 2PCF can improve the results.Therefore, in this section, we analyse the outcomes of the HOD fitting performed on the 2PCF monopole and quadrupole, for  ∈ [ min , 50] Mpc/ℎ, see Table 3.
Figure 5 presents the monopole, quadrupole and hexadecapole of the 2PCF computed for  ∈ [0, 200] Mpc/ℎ as well as the tensions between the FastPM and SLICS.The FastPM clustering typically falls within 2 of the reference for scales larger than 50 Mpc/ℎ and is largely unaffected by the fitting scenario.However, the HR monopoles are consistently closer to the reference than LR monopoles by approximately 0.5 at scales larger than ≈ 150 Mpc/ℎ.
Including the smallest scales (Large interval) in the HOD fitting, we observe a 1 to 2 agreement with the reference for  < 10 Mpc/ℎ in both the monopole and quadrupole.However, at intermediate scales  ∈ [10, 50] Mpc/ℎ, the monopole is significantly biased, exhibiting a deviation of 3.In contrast, for the Medium and Small scenarios, we notice that the tensions for the monopole and quadrupole at intermediate scales drop to 1, while the smallest scales can get biased by more than 3.Nevertheless, they match better the reference than the power spectrum HOD fitting case.Lastly, the hexadecapole does not depend on the resolution nor the fitting range and is strongly biased for  < 60 Mpc/ℎ, showing no improvement compared to the power spectrum fitting.
As in the previous subsection, we test the clustering statistics of the best-fitting FastPM boxes that were not included in the HOD fitting, i.e. the power spectrum in the Figure 5.The first observation is that these FastPM power spectra do not fit as well the reference as the ones from Figure 3. On one hand, for the HR case and Medium and Small fitting intervals a ±1 matching is possible up to  = 0.4 ℎ/Mpc  2) between them: left -power spectrum and right -2PCF.FastPM mocks share the white-noise through the initial conditions with the SLICS ones.The fitting has been performed: 1) on the monopole and quadrupole of the power spectrum; 2) for three different fitting ranges, see Table 3; 3) using HR (dashed) and LR (continuous) FastPM realisations.The O axis of the 2PCF panels has a linear scale from 0 to 50 Mpc/ℎ and a logarithmic scale above this limit.
and  = 0.3 ℎ/Mpc, respectively.On the other hand, the LR FastPM allows a good matching up to  ≈ 0.2 ℎ/Mpc for the same fitting intervals.While the  min = 0 case has a good matching quadrupole up to  ≈ 0.4 ℎ/Mpc, its monopole follows similar trend to the 2PCF monopole, i.e. the intermediate scales  ∈ [0.25, 0.4] ℎ/Mpc are biased and the rest are mostly within 2 deviation.Lastly, the hexadecapole is within ±2 up to  ≈ 0.3 ℎ/Mpc for the Large fitting interval and up to  ≈ 0.2 ℎ/Mpc for the other cases.
A quantitative evidence that directly fitting the 2PCF yields superior matching of the 2PCF compared to fitting the power spectrum is displayed in Figure 4.The majority of the  2  values in the lower-right panel are lower compared to those in the upper-right panel.Furthermore, fitting on the Small interval ( min = 10), the  2  ≈ 1, indicating that the 2PCF is in good agreement with the SLICS reference above a certain minimum separation.The almost constant  2  for the Large fitting interval in the lower-right panel of Figure 4 is explained by the discrepancy at the intermediate scales of the monopole for the Large fitting interval in Figure 5. Lastly, as in the previous fitting scenario, the HR FastPM generally provides a lower  2  than the LR.In con-trast, only the HR simulations can provide a  2  < 2 to both the 2PCF and the power spectra, and only when fitting with the Medium and Small intervals to 2PCF.Although not shown in the aforementioned figure, it is important to note that  2  [(0.02,0.2)] = 2.4 for the 2PCF Small interval LR case.

Bi-spectrum comparison
The two-point clustering covariance matrix is directly related to the tri-spectrum, i.e. the four-point clustering, however this is difficult to tune.Nonetheless, Baumgarten & Chuang (2018) have shown that for similar two-point clustering, the corresponding covariance matrices -at scales of  < 40 Mpc/ℎ -are sensitive to changes in the bi-spectrum.Consequently, even though we do not include the bi-spectrum into the HOD fitting, we aim to understand its behaviour when incorporating various scales of the two-point clustering in the HOD fitting.
Figure 6 compares the average bi-spectrum of the 20 best-fitting FastPM boxes with the one computed on the corresponding SLICS  as defined in Section 3.3.3.We compute  2  : 1) for different intervals (see O and O axes) of the clustering statistics (left panels -power spectrum; right panels -2PCF); 2) for different fitted clustering (upper panels -power spectrum, see Section 4.1; lower panels -2PCF, see Section 4.2); 3) for different fitting ranges (see Table 3).boxes, for 2  1 =  2 = 0.2 ℎ/Mpc configuration as done by Baumgarten & Chuang (2018).These scales are chosen because they are sensitive to BAO and RSD, the primary scientific case of DESI.It is evident that by increasing the fitting range to include smaller scales, the FastPM bi-spectrum changes to the extent that for  max = 0.5, the tension ranges from 1 to 2.In contrast, when fitting the 2PCF the resulting bi-spectrum is more biased, i.e. the lowest deviation is ≈ 5, for  min = 0 case.In addition, we have checked multiple configurations of the bi-spectrum with   ∈ [0.01, 0.21] ℎ/Mpc and  = 1 or 2 and observed that for  max = 0.5 and  min = 0, the overall agreement between the FastPM and SLICS is around 1 and 5 per cent, respectively.Lastly, there is no significant improvement in terms of the goodness-of-fit between the HR and LR.
In the previous sections, we compare the HR and LR FastPM with SLICS using the two-point clustering of the 20 cubic mocks included in the HOD fitting.The HR simulations perform better than LR to model the extremely non-linear scales, such as  ≈ 0.5 ℎ/Mpc,  ≈ 0 Mpc/ℎ.In contrast, at mildly non-linear scales ( ≈ 0.3 ℎ/Mpc,  ≈ 10 Mpc/ℎ) that are more relevant to BAO and RSD analyses (e.g.Tamone et al. 2020;de Mattia et al. 2021), LR and HR show similar performance.Moreover, Figure 6 suggests that the bi-spectrum does not depend strongly on the resolution.Nevertheless, the computing cost of HR is significantly higher than for LR and given the small difference in precision, we argue it is optimal to use LR FastPM for further analyses.
Furthermore, in Figure 7, we compare the average bi-spectracomputed from 778 LR FastPM realisations -corresponding to the six HOD fitting cases, see Table 3.In this and the next figures, we choose the  max = 0.5 case as reference because: (i) Figure 4 shows that the best-fitting power spectrum provides  2  ≈ 1; (ii) Figure 6 implies that the corresponding bi-spectrum is the closest to the SLICS reference.
One can notice that  min = 0 bi-spectrum is at most 5 per cent different than the reference, while the rest can reach 15 per cent discrepancies.The  max = 0.3 and  max = 0.   3. The reference for the bi-spectra ratios is the average bi-spectrum computed for the  max = 0.5 case.The bispectra are computed for  1 = 0.1 ± 0.05 ℎ/Mpc and  2 = 0.2 ± 0.05 ℎ/Mpc, with the  angle between k 1 and k 2 varying from 0 to .
per cent different from each other and similarly for  min = 5 and  min = 10.

Covariance matrix comparison
Having studied the behaviour of the bi-spectra with respect to the scales introduced into the HOD fitting, we now want to understand the resulting impact on the two-point clustering covariance matrices.Firstly, Figure 8 compares the  max = 0.5 and SLICS correlation matrices and standard deviations.Given the limited number of realisations used to estimate the correlation matrices, it is difficult to assess the importance of the differences.Nevertheless, there seems to be more significant differences at very small scales for the 2PCF and important differences at almost all scales for the power spectrum.Numerically, the standard deviations deviate by more than 7 to 10 per cent from the SLICS case for  < 10 Mpc/ℎ and  > 0.35 ℎ/Mpc.A more detailed comparison is performed by Chuang et al. (2023).Secondly, we compare the correlation matrices and standard deviations obtained from 778 FastPM realisations, of the six HOD fitting cases, where  max = 0.5 is the reference.

Power spectrum covariance matrix
Figure 9 presents the correlation matrices and the corresponding standard deviations  ℓ for the monopole and quadrupole of the power spectrum.The following pairs ( max = 0.4,  max = 0.3), ( min = 5,  min = 10) and ( min = 0,  max = 0.5) have very similar correlation matrices, thus we only show three cases.However, we introduce all of them in Appendix C.
The similarity to the reference correlation matrix diminishes in the following order:  min = 0,  max = 0.4,  max = 0.3,  min = 5 and  min = 10.However, for the largest scales of the quadrupole ( < 0.15 ℎ/Mpc), the correlation coefficients are practically the same for all cases.
The standard deviations in the lowest panels show similar trends.The  min = 0 case is within two per cent of the reference case.The In the first row, the upper triangular matrices display the correlation matrices of the FastPM  max = 0.5 case, while the lower triangular ones show the SLICS correlation matrices.The second row of panels contains the differences between the SLICS and FastPM correlation matrices.The third row of panels illustrate the ratios between the standard deviations.The shaded regions denote two and five per cent limits.
max = 0.4,  max = 0.3 cases overestimate the  ℓ () by approximately two per cent for  < 0.27 ℎ/Mpc and by ≈ 5 per cent for smaller scales.Nevertheless, these two cases are consistent with each other within one to two per cent.In contrast,  min = 5 and  min = 10 can overestimate the  ℓ () by ≈ 2 to 5 per cent for  < 0.27 ℎ/Mpc and by 10 to 20 per cent for smaller scales.These two cases are also consistent with each other for most scales, except for the quadrupole  > 0.3 ℎ/Mpc.These findings are in agreement with the trends observed in the bi-spectrum comparison in Figure 7.
In order to quantify the differences between the covariance matrices we adopt the method described in Section 3.4 and thus obtain the results displayed in Figures 10 and 11.The first one reveals that none of the six covariance matrices bias the two fitting parameters, regardless of the fitting range.We only present here the results of one fitting range, however all cases can be found in Appendix C.
Examining the uncertainty on  0 in Figure 11, we observe that including the smaller scales the discrepancy between the error estimates of the six covariance matrices increases, as we expect from Figure 9, reaching a maximum of ≈ 20 per cent larger error estimation for the  min = 10 covariance at K = 0.25 ℎ/Mpc.Moreover, each of the following pairs ( max = 0.4,  max = 0.3), ( min = 5,  min = 10) and ( min = 0,  max = 0.5) provide coherent estimations of the uncertainty, which is consistent with the observations on the correlation matrices and standard deviations.Lastly, a five per cent  (right).Given the similarity between some correlation matrices, only three out of six different HOD fitting cases -see Table 3 -are presented here.However, all cases can be found in Appendix C. The reference case corresponds to  max = 0.5.The upper triangular matrices display the correlation matrices, while the lower triangular ones show the differences between the correlation matrices and the reference one.The bottom panels illustrate the ratios between the standard deviations.The shaded regions denote two and five per cent limits.
consensus between all six covariance matrices is achieved when we fit the power spectra on the  ∈ [0.02, 0.1] ℎ/Mpc.
The agreement between covariance matrices on   2 is much better than   0 .Given the error bars, the six methods estimate the uncertainty with a two per cent tolerance with each other for all K values.
Finally, all six covariance matrices provide values of R [ 0 ,  2 ] that are consistent at the level of 5 per cent, given the error bars and up to K = 0.2 ℎ/Mpc.For K = 0.25 ℎ/Mpc, the largest discrepancy is shown by  min = 10 case which underestimates the value

2PCF covariance matrix
Comparing the correlation matrices obtained from 778 2PCF in Figure 9, one can observe that the largest differences occur at the smallest scales  < 30 Mpc/ℎ.Similarly to the power spectrum correlation matrices, the same pairs of cases show resembling behaviours at all scales.Equivalent qualitative comments can be made about the ratios of the standard deviations.Nonetheless, all cases are within ≈ 2 per cent from each other for  > 30 Mpc/ℎ, while at smaller scales, the differences can get larger than ≈ 20 per cent.
Following the method described in Section 3.4, we obtain the results shown in Figures 10 and 12.The first figure proves that all six covariance matrices provide unbiased measurements of  ℓ parameters.
Resembling the power spectrum fitting case, the six estimations of   0 in Figure 12   or the 2PCF can be different on highly non-linear scales, therefore an analysis using such scales needs to be more careful with the covariance matrix estimation.In contrast, in the quasi-linear regime ( ≳ 20 Mpc/ℎ;  ≲ 0.2 ℎ/Mpc) the differences are not statistically significant (see e.g. min = 0 and  max = 0.5 cases) , thus one can use either space to obtain the covariance matrices and apply them to both analyses.

CONCLUSIONS
We have implemented an HOD model to assign galaxies on the FastPM halo cubic mocks, such that the resulting clusteringmonopole and quadrupole -matches the SLICS reference one.In order to remove the cosmic variance, we have used 20 SLICS galaxy catalogues and 20 halo FastPM mocks (low resolution or high resolution) that share the initial conditions with the SLICS simulations.
Given the shared white noise, the standard covariance matrix is obsolete, thus we have performed a two-steps HOD fitting: (i) use a simple diagonal covariance matrix to get Initial-Guess best-fitting FastPM galaxy mocks; (ii) compute the covariance matrix of the 123 realisations of the difference between the IG-FastPM and the SLICS clustering, and use it to perform the final HOD fitting.
On one hand, the HR FastPM generally performs better than the LR at modelling the SLICS clustering.On the other hand, LR is also able to provide a  2  ≈ 1 for  max = 0.5,  max = 0.4,  max = 0.3 and  min = 10.The  max = 0.5 case is one of the most valuable as it additionally offers 2 matching: (i) power spectrum hexadecapole for  < 0.4 ℎ/Mpc; (ii) 2PCF monopole and quadrupole for  > 10 Mpc/ℎ; (iii) bi-spectrum.
Nevertheless, fitting the 2PCF with  min = 10, produce a 1 matching power spectrum monopole and quadrupole for  ≲ 0.2, but a strongly biased bi-spectrum.In a similar way as the power spectrum, one must include the smallest scales to better reproduce the SLICS bi-spectrum, i.e. for  min = 0 the bi-spectrum tension drops from 20 to 5.As a general remark, the power spectrum hexadecapoles can be slightly tuned by changing the values of  max or  min , but the 2PCF hexadecapole is practically independent on the fitting range.
The fact that there is no HOD fitting case that provides both power spectrum and 2PCF 1 matching with the reference might raise the question whether one needs two sets of catalogues for the DESI analysis.Therefore, it could be interesting for future studies to perform a joint fitting of both Fourier and Configuration clustering statistics to test for possible improvements in modelling non-linear scales.Future studies could test also the necessary level of agreement between the FastPM and the full -body reference in order to have a similar behaviour when observational systematic effects are included.Nonetheless, for the purpose of this article, in the second part of the study, we have exposed and compared the resulting covariance matrices.
Firstly, we have briefly shown with the 123 FastPM realisations -that share the initial conditions with SLICS -that the standard deviations of the  max = 0.5 case agree within five per cent with the SLICS reference case for  > 10 Mpc/ℎ and  < 0.35 ℎ/Mpc.Secondly, we have focused on the 778 LR FastPM realisations corresponding to the six best-fitting cases, where  max = 0.5 is considered the reference.
Since Baumgarten & Chuang (2018) have shown that the bispectrum computed with a configuration that includes the BAO and RSD effects -i.e. 2 ×  1 =  2 = 0.2 ℎ/Mpc -affects the covariance matrix, we have compared the six bi-spectra and checked how this comparison reflects into the two-point clustering covariance matrices.Lastly, we have examined the constraining power of these covariance matrices using a simplified clustering model with two scaling parameters, i.e.  0 and  2 for the monopole and quadrupole.We have focused on fitting intervals similar to the ones used in standard BAO and RSD analyses i.e.K ≲ 0.20 ℎ/Mpc and S ≳ 20 Mpc/ℎ (e.g.Tamone et al. 2020;de Mattia et al. 2021).
The  min = 0 bi-spectrum is at most five per cent different than the  max = 0.5, while the other cases can reach a discrepancy of 15 per cent.However, each of these pairs ( max = 0.4,  max = 0.3), ( min = 5,  min = 10) yield similar bi-spectra.These observations are in a good agreement with a qualitative description of the shown correlation matrices and the standard deviations.
Quantitatively, the power spectrum standard deviations of  min = 0 and ( max = 0.4,  max = 0.3) are within two percent from the reference for  < 0.5 ℎ/Mpc and  < 0.27 ℎ/Mpc, respectively.Furthermore, the 2PCF standard deviations of all cases are within two percent from the reference for  > 30 Mpc/ℎ.
Using the simplified clustering model,  0 and  2 are measured accurately for both power spectrum and 2PCF using all six covariance matrices.The six estimations of   0 from the power spectrum fitting up to K = 0.20 ℎ/Mpc are scattered within at most 20 per cent from the reference, whereas the values of   2 are within two per cent agreement, given the error bars.Lastly, the covariances between  0 and  2 are scattered within 5 per cent from the reference.
In contrast, the estimations of   0 from the 2PCF fitting down to S = 20 Mpc/ℎ are found within five per cent from each other.Similarly to the power spectrum case, the   2 values agree at the level of two per cent.Given the error bars, the covariances between  0 and  2 are consistent at the level of five per cent.
Lastly, analysing individual cases in the following pairs of HOD fitting cases: ( min = 0,  max = 0.5), ( max = 0.4,  max = 0.3), ( min = 5,  min = 10), one can observe that each pair provides similar uncertainty estimation   ℓ , irrespective of the level of matching between the FastPM and SLICS two-point clustering.In order words, it seems that the estimations of the uncertainties are closer when the bi-spectra -computed for 2 ×  1 =  2 = 0.2 ℎ/Mpc -are more similar.The most striking case is the 2PCF for  min = 0 which is more than 2 different than the  max = 0.5 2PCF, for  < 40 Mpc/ℎ, but the 2PCF standard deviations and the   ℓ estimations are practically identical.
In conclusion, one can use an HOD model on the low resolution FastPM halo catalogues to tune the galaxy clustering such that it matches the SLICS reference down to certain minimum scales.Additionally, the HOD fitting intervals can have an impact on the final FastPM based covariances.This influence is observed as a scatter in the uncertainty estimation of up to 20 per cent for power spectrum and five per cent for 2PCF at the scales interesting for BAO and RSD analyses.Nevertheless, more accurate analyses could be performed in the future using actual BAO and RSD models and a larger sample of mocks, such as AbacusSummit.3.While the left panels include monopole and quadrupole of the 2PCF, the right ones display the power spectrum. as defined in Section 3.3.3,but using different covariance matrices.We compute  2  : 1) for the six best-fitting FastPM cases, three cases for the power spectrum in the left panel (see Section 4.1), and three cases for the 2PCF in the right panel (see Section 4.2); 2) with the six difference covariance matrices obtained after the second HOD fitting step (the coloured dots).The black points show the best fitting  2  for the six cases that also appear in Figure 4.
Additionally, the HOD fitting is computationally expensive (≈ 6000 CPUhours), thus we are not able to perform hundreds of HOD fittings with different covariance matrices.Consequently, after obtaining one set of best-fitting HOD parameters, we computed the  2  with the same best-fitting FastPM clustering, but with  cov mocks different covariance matrices.The covariance matrices -Σ   , with  from 1 to  cov mocks -are estimated using Eq. ( 22), but with only  cov mocks − 1 clustering realisations.Furthermore, we compute  2, ,JK for each Σ   , as defined in Eq. ( 21) and we calculate the mean χ2  and their uncertainties introduced by the covariance matrix (Eqs.(B1) and (B2)) and the finite number of HOD realisations per FastPM halo catalogue (Eqs.(B3) and (B4)).The estimations have been performed on the LR (1296 3 ) FastPM galaxy catalogues and on both the 2PCF and the power spectrum for the three specific fitting ranges defined in Table 3.

B2 HOD induced uncertainty
During the HOD fitting, for each FastPM halo catalogue we create a single galaxy realisation, in order to reduce the optimisation time.As a consequence, we introduce additional noise in the HOD fitting process, that is not considered in the covariance matrix.
With the aim of estimating the effect of this noise on the  2  , we compute 100 galaxy realisations for a given set of best-fitting HOD parameters and per FastPM halo catalogue.Furthermore, using the 20 galaxy realisations corresponding to the halo catalogues used in the HOD fitting process and the same covariance matrix, we compute  2, ,HOD as in Eq. ( 21), where  = 1, ..., 100.Finally, we calculate the mean and the standard deviation of the 100  2, ,HOD values:

Figure 2 .
Figure 2. Monopole and quadrupole error bars: left panels -2PCF; right panels -power spectrum.Black -estimated uncertainty for the entire DESI survey (Y5); Dotted black -estimated uncertainty for the Year 1 (Y1) DESI survey; Blue -square root of the Σ diag 's terms; Dashed orange -square root of the Σ diff 's diagonal terms, divided by √ 20,  fit mocks = 20; Dashed redstandard deviation of the differences between the best-fitting FastPM clustering (from the second HOD fitting step, one HOD fitting scenario) and SLICS ( fit mocks realisations), further divided by √︃  fit mocks ; Grey -standard deviation of  fit mocks SLICS clustering realisations, further divided by √︃  fit mocks .We emphasise that the Y5 and Y1 errors are estimated by scaling the SLICS covariance to the effective volumes of 24 Gpc 3 ℎ −3 and 8 Gpc 3 ℎ −3 , respectively.

Figure 3 .
Figure3.The average of 20 SLICS (reference-black) and 20 FastPM (model-colours) clustering realisations and the tension ( DIFF is shown in Figure2) between them: left -power spectrum and right -2PCF.FastPM mocks share the white-noise through the initial conditions with the SLICS ones.The fitting has been performed: 1) on the monopole and quadrupole of the power spectrum; 2) for three different fitting ranges, see Table3; 3) using HR (dashed) and LR (continuous) FastPM realisations.The O axis of the 2PCF panels has a linear scale from 0 to 50 Mpc/ℎ and a logarithmic scale above this limit.

Figure 5 .Figure 6 .
Figure 5. Same as Figure 3, but the fitting is done on the monopole and quadrupole of the 2PCF.

Figure 7 .
Figure 7. Average bi-spectra computed using 778 LR FastPM realisations for different HOD fitting cases, see Table3.The reference for the bi-spectra ratios is the average bi-spectrum computed for the  max = 0.5 case.The bispectra are computed for  1 = 0.1 ± 0.05 ℎ/Mpc and  2 = 0.2 ± 0.05 ℎ/Mpc, with the  angle between k 1 and k 2 varying from 0 to .

Figure 8 .
Figure8.The correlation matrices and the standard deviations computed using 123 realisations of the monopole and quadrupole of the power spectrum (left) and 2PCF (right).In the first row, the upper triangular matrices display the correlation matrices of the FastPM  max = 0.5 case, while the lower triangular ones show the SLICS correlation matrices.The second row of panels contains the differences between the SLICS and FastPM correlation matrices.The third row of panels illustrate the ratios between the standard deviations.The shaded regions denote two and five per cent limits.

Figure 9 .
Figure9.The correlation matrices and the standard deviations computed using 778 realisations of the monopole and quadrupole of the power spectrum (left) and 2PCF (right).Given the similarity between some correlation matrices, only three out of six different HOD fitting cases -see Table3-are presented here.However, all cases can be found in Appendix C. The reference case corresponds to  max = 0.5.The upper triangular matrices display the correlation matrices, while the lower triangular ones show the differences between the correlation matrices and the reference one.The bottom panels illustrate the ratios between the standard deviations.The shaded regions denote two and five per cent limits.

Figure 10 .
Figure 10.The average of the 123 fitting parameters ( ℓ ) obtained from 123 SLICS clustering realisations, as described in Section 3.4.The measurements performed on the power spectra with  ∈ [0.02,K ] ℎ/Mpc are show on the left, while those based on the 2PCF with  ∈ [ S, 200] Mpc/ℎ are depicted on the right.The error bars are computed as the average of 123   ℓ , divided by √ 123.Here,   ℓ represents the standard deviation of the  ℓ posterior distribution.The colours correspond to the different FastPM covariance matrices illustrated in Figure 9. Due to the similar results, only one value for K and S are shown here.However, all tested values are presented in Appendix C

Figure 11 .
Figure 11.The averages of 123   ℓ and 123 R [ 0 ,  2 ] -the standard deviation of the  ℓ posterior distribution and the covariance between  0 and  2 , respectively, detailed in Section 3.4 -obtained from 123 SLICS power spectra fitted on  ∈ [0.02,K ] ℎ/Mpc.In order to estimate the error bars, we split the 778 FastPM realisations in six distinct sets of 123 realisations and compute for each set  a covariance matrix Σ  123,FastPM with which we fit the 123 SLICS clustering realisations.Having obtained 123 values of    ℓ per set, we compute their average σ .Finally, the error bars are the standard deviation of the six

Figure A1 .
Figure A1.The ratios between the standard deviations computed on the differences of  cov mocks = 123 (LR FastPM, SLICS) clustering pairs and the square root of the diagonal of Σ diff , i.e.Σ  diff .The colours denote the HOD fitting scenarios in the second HOD fitting step, see Table3.While the left panels include monopole and quadrupole of the 2PCF, the right ones display the power spectrum.

Figure A2 .
Figure A2.The  2 as defined in Section 3.3.3,but using different covariance matrices.We compute  2  : 1) for the six best-fitting FastPM cases, three cases for the power spectrum in the left panel (see Section 4.1), and three cases for the 2PCF in the right panel (see Section 4.2); 2) with the six difference covariance matrices obtained after the second HOD fitting step (the coloured dots).The black points show the best fitting  2  for the six cases that also appear in Figure4.
and C2 one can observe that the ( min = 5,  min = 10) and ( max = 0.3,  max = 0.4) pairs have very similar correlation matrices.Consequently, we only show  max = 0.3 and  min = 10 in the main text.FiguresC3 and C4show the results of fitting the clustering with the simplified model detailed in Section 3.4.Since most results are consistent with the expected value of one, we only display the values for K = 0.25 ℎ/Mpc and S = 20 Mpc/ℎ in the main text.This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure C1 .Figure C3 .Figure C4 .
Figure C1.Upper triangular matrices: correlation matrices of the power spectrum monopole and quadrupole, for the fitting cases defined in Table3.Lower triangular matrices: the difference between the shown correlation matrix and the reference one, i.e.  max = 0.5.

Table 1 .
It contains a summary of important symbols related to the HOD fitting.With the aim of finding the best-fitting FastPM clustering, we run