Assessing the accuracy of cosmological parameters estimated from velocity -- density comparisons via simulations
MMNRAS , 1–10 (2020) Preprint 13 January 2021 Compiled using MNRAS L A TEX style file v3.0
Assessing the accuracy of cosmological parameters estimated fromvelocity – density comparisons via simulations
Amber M. Hollinger , and Michael J. Hudson , , Department of Physics and Astronomy, University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada Waterloo Centre for Astrophysics, University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada Perimeter Institute for Theoretical Physics, 31 Caroline St. North, Waterloo, ON N2L 2Y5, Canada
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
A promising method for measuring the cosmological parameter combination 𝑓 𝜎 is tocompare observed peculiar velocities with peculiar velocities predicted from a galaxy densityfield using perturbation theory. We use 𝑁 -body simulations and semi-analytic galaxy formationmodels to quantify the accuracy and precision of this method. Specifically, we examine anumber of technical aspects, including the optimal smoothing length applied to the densityfield, the use of dark matter halos or galaxies as tracers of the density field, the effect of noisein the halo mass estimates or in the stellar-to-halo mass relation, and the effect of finite surveyvolumes. We find that for a Gaussian smoothing of 4 ℎ − Mpc, the method has only smallsystematic biases at the level of 5%. Cosmic variance affects current measurements at the 5%level due to the volume of current redshift data sets.
Key words:
Galaxy: kinematics and dynamics – galaxies: statistics – large-scale structure ofUniverse – cosmology: observations
Peculiar velocities are the only practical way of measuring the un-derlying distribution of dark matter (DM) on large scales (Willicket al. 1997), in the nearby (low redshift) Universe. Specifically, theycan be used to measure the cosmological parameter combination 𝑓 ( Ω 𝑚 ) 𝜎 , where the first term is the logarithmic growth rate offluctuations, with 𝑓 ( Ω 𝑚 ) = Ω . 𝑚 in the Lambda Cold Dark Mat-ter ( Λ CDM) model, and 𝜎 is the root mean square matter densityfluctuation in a sphere of radius 8 ℎ − Mpc. Although peculiar ve-locities have been used for similar purposes since the early 90s (seereviews by Dekel (1994) and Strauss & Willick (1995)), there hasbeen a recent revival of interest in measuring 𝑓 𝜎 (Pike & Hudson2005; Davis et al. 2011; Turnbull et al. 2012; Hudson & Turnbull2012; Carrick et al. 2015; Huterer et al. 2017; Dupuy et al. 2019;Qin et al. 2019; Adams & Blake 2020; Said et al. 2020; Boruahet al. 2020).The use of peculiar velocities as a probe of 𝑓 𝜎 has takenon renewed importance in light of the 3.2 𝜎 conflict between Cos-mic Microwave Background (Planck Collaboration 2018) and weaklensing measurements (Asgari et al. 2019) of the parameter com-bination 𝑆 ≡ Ω . 𝑚 𝜎 , with the latter giving lower values. Thiscombination is very similar to 𝑓 𝜎 , differing by only 0.05 in theexponent of Ω 𝑚 . As shown by Boruah et al. (2020), some recentpeculiar velocity measurements also yield lower values of 𝑆 thanPlanck (as do most other probes of this parameter combination).One method for measuring 𝑓 𝜎 involves regressing the ob-served peculiar velocities (from e.g. supernovae, Tully-Fisher or Fundamental Plane standard candles/rulers) on their predicted pe-culiar velocities from the density field of galaxies, obtained from aredshift survey. Specifically, the slope of this regression, 𝛽 𝑔 , is thencombined with a measurement of the fluctuations of galaxies, 𝜎 ,𝑔 to yield an estimate of 𝑓 𝜎 = 𝛽 𝑔 𝜎 ,𝑔 , (1)as first suggested by Pike & Hudson (2005). The background theoryand assumptions implicit in this method are discussed in more detailin Section 2 below.The goal of this paper is to test the accuracy and precisionof this method using large cosmological 𝑁 -Body simulations at aredshift of zero of DM halos as a proxy for galaxies, as well assemi-analytical models of galaxy formation. Our approach is notto simulate all the observational properties of the surveys simul-taneously (as in, for example, Nusser et al. (2014)), but rather toconsider one by one the different physical effects which may biasor add uncertainty to our results. In this sense the paper is similarin spirit to Berlind et al. (2000).This paper is organized as follows. Section 2 describes the the-oretical framework of the relationship between peculiar velocitiesand the density fluctuation field from linear perturbation theory, andthe relevant cosmological parameters. Section 3 presents the simu-lation data and semi-analytic models that are used for the analysesperformed in this paper. Section 4 describes the prescription for howwe predict peculiar velocities. In section 5 we investigate how thechoice of smoothing kernel impacts our estimates of 𝛽 / 𝑓 and the © a r X i v : . [ a s t r o - ph . C O ] J a n Hollinger and Hudson amount of scatter generated in velocity – density cross-correlations.Section 6 investigates how uncertainties corresponding to 0.1 and0.2 dex in halo mass measurements influence the predictions of 𝛽 / 𝑓 and 𝜎 . In section 7 we instead explore how using either the stellarto halo mass relation or galaxy observables to weight the densityfield impacts these cosmological estimates. Section 8 focuses onhow these estimates are affected by volume limited surveys. Finally,Section 9 presents our conclusions. In linear perturbation theory, it is possible to relate the density fieldto the peculiar velocities of the galaxies at low redshift using 𝑣 ( 𝒓 ) = 𝐻 𝑓 ( Ω 𝑚 ) 𝜋 ∫ 𝛿 ( 𝒓 (cid:48) ) ( 𝒓 (cid:48) − 𝒓 )| 𝒓 (cid:48) − 𝒓 | 𝑑 𝒓 (cid:48) , (2)where 𝑣 ( 𝒓 ) is the peculiar velocity field and 𝛿 ( 𝒓 ) is the matterdensity fluctuation field given by 𝛿 ( 𝒓 ) = 𝜌 ( 𝒓 ) − ¯ 𝜌 ¯ 𝜌 , (3)where 𝜌 ( 𝒓 ) is the matter density field and ¯ 𝜌 is its cosmic average.This calculation is only valid in the linear regime, where 𝛿 (cid:46) Λ CDM, the rms matter density in spheresincreases with decreasing sphere radius, so we expect linear theoryto break down on small scales. In practice then, to apply equation 4one needs to smooth 𝛿 ( 𝑟 ) .In the linear regime, the density modes in Fourier space growindependently of one another. As a result it is easier to write equation2 in Fourier space as follows 𝒗 𝒌 = 𝑖𝐻 𝑓 𝒌 | 𝒌 | 𝛿 𝑘 , (4)where 𝐻 is the Hubble constant ( 𝐻 = ℎ km s − Mpc − ). Thesmoothing is also simpler in Fourier space, because a convolutionis a multiplication in Fourier space.The density fluctuation field used in the previous equations isthat of the underlying matter density field. Because it is dominatedby dark matter, this cannot not be measured empirically. Instead, anassumption must be made as to how the observed galaxies trace theunderlying total matter. If one assumes linear biasing, the relationis 𝛿 𝑔 = 𝑏 𝑔 𝛿, (5)where 𝑏 𝑔 is the linear bias and 𝛿 𝑔 is the density fluctuation field ofthe galaxies . Under this assumption equation 2 can be written as 𝑣 ( 𝒓 ) = 𝐻 𝑓 ( Ω 𝑚 ) 𝜋𝑏 𝑔 ∫ 𝛿 𝑔 ( 𝒓 (cid:48) ) ( 𝒓 (cid:48) − 𝒓 )| 𝒓 (cid:48) − 𝒓 | 𝑑 𝒓 (cid:48) . (6)Note that if we express distances 𝒓 in units of ℎ − Mpc (or km/s),as are naturally obtained from redshift surveys, then one must set 𝐻 =
100 km s − Mpc − (or 1, respectively) in the above expression.Thus when applying this equation to density fields derived fromredshift surveys, the predictions are independent of the true valueof 𝐻 . The other two values outside the integral can be compactedinto the parameter combination 𝛽 𝑔 ≡ 𝑓𝑏 𝑔 . (7) If linear biasing holds, then 𝜎 ,𝑔 = 𝑏 𝑔 𝜎 . Putting this into equa-tion 7, we find that the product of the observables to be 𝛽 𝑔 𝜎 ,𝑔 ,as in equation 1, and so we can set constraints on cosmologicalparameters.In reality, the assumptions of linearity, both in the context ofperturbation theory and in the context of biasing, will not holdexactly. Fortunately, peculiar velocities are primarily generated bylarge scale waves, where linearity will be a good approximation.Nevertheless, they also have a contribution from smaller, less linearscales. For these reasons, simulations are needed to calibrate anybiases arising from non-linearities. We use two publicly available simulations: Bolshoi and MultiDarkPlanck 2 (MDPL2) (Klypin et al. 2011, 2016), along with two sim-ulated semi-analytic catalogues, SAG (Cora et al. 2018) and SAGE(Croton et al. 2016) which populate the DM haloes of MDPL2. Weuse the snapshots at 𝑧 =
0. The halo and galaxy catalogues wereobtained from the COSMOSIM database .The high-resolution Bolshoi simulation (Klypin et al. 2011)follows 2048 particles in a comoving, periodic cube of length 250 ℎ − Mpc from 𝑧 =
80 to today. It has a mass and force resolution of1 . × ℎ − 𝑀 (cid:12) and 1 ℎ − kpc respectively, and the DM haloesrange from the masses of MW satellites (10 𝑀 (cid:12) ) to the largestof clusters (10 𝑀 (cid:12) ). It was run as a collisionless DM simulationwith the Adaptive Refinement Tree Code (ART; Kravtsov et al.(1997)) and assumes a flat, WMAP5 cosmology with parameters Ω 𝑚 = . Ω Λ = . ℎ = .
7, (linear) 𝜎 = .
82 and 𝑛 𝑠 = . 𝜎 ,𝑚 , measured from the particles is 0.897.The MultiDark project consists of a suite of cosmologi-cal hydrodynamical simulations (Klypin et al. 2016), all assumea flat Λ CDM cosmology with cosmological parameters: Ω 𝑀 = . Ω Λ = . ℎ = . 𝜎 = . 𝑛 𝑠 = .
96, which is consistent with Planck results (Planck Collab-oration 2018). We focus on the MDPL2 simulation which has aperiodic box of length 1000 h − Mpc evolved from a redshift of 120to 0 with a varying physical force resolution level from 13-5 ℎ − kpc and various implemented physics. The simulation uses 3840 dark matter particles of mass 1 . × ℎ − 𝑀 (cid:12) , and has identifiedmore that 10 haloes using Rockstar (Behroozi et al. 2013a), withmerger trees that were generated using ConsistentTrees (Behrooziet al. 2013b). The (non-linear) 𝜎 ,𝑚 , measured from the particles is0.95.The SAG (Cora et al. 2018) and SAGE (Croton et al. 2016)semi-analytic models include the most relevant physical processesin galaxy formation and evolution, such as radiative cooling, starformation, chemical enrichment, supernova feedback and winds,disc instabilities, starbursts, and galaxy mergers. These model werecalibrated to generate galaxy catalogues using the MDPL2 simula-tion. A comprehensive review of the models can be found in Knebeet al. (2017). Our goal is to test equation (6) using data from N-body simulations.Specifically, we will calculate the predicted peculiar velocities by MNRAS000
96, which is consistent with Planck results (Planck Collab-oration 2018). We focus on the MDPL2 simulation which has aperiodic box of length 1000 h − Mpc evolved from a redshift of 120to 0 with a varying physical force resolution level from 13-5 ℎ − kpc and various implemented physics. The simulation uses 3840 dark matter particles of mass 1 . × ℎ − 𝑀 (cid:12) , and has identifiedmore that 10 haloes using Rockstar (Behroozi et al. 2013a), withmerger trees that were generated using ConsistentTrees (Behrooziet al. 2013b). The (non-linear) 𝜎 ,𝑚 , measured from the particles is0.95.The SAG (Cora et al. 2018) and SAGE (Croton et al. 2016)semi-analytic models include the most relevant physical processesin galaxy formation and evolution, such as radiative cooling, starformation, chemical enrichment, supernova feedback and winds,disc instabilities, starbursts, and galaxy mergers. These model werecalibrated to generate galaxy catalogues using the MDPL2 simula-tion. A comprehensive review of the models can be found in Knebeet al. (2017). Our goal is to test equation (6) using data from N-body simulations.Specifically, we will calculate the predicted peculiar velocities by MNRAS000 , 1–10 (2020) esting Accuracy of Peculiar Velocity Comparisons integrating over a smoothed density tracer field, denoted by thesubscript “t”, used in place of 𝛿 𝑔 ( 𝒓 ) in the integral, 𝒗 pred ( 𝒓 ) = 𝐻 𝜋 ∫ 𝛿 𝑡 ( 𝒓 (cid:48) ) ( 𝒓 (cid:48) − 𝒓 )| 𝒓 (cid:48) − 𝒓 | 𝑑 𝒓 (cid:48) , (8)and compare this to the “observed” (unsmoothed) velocity tracers,which are used in place of 𝒗 ( 𝒓 ) on the left hand side of equation(6). Note that equation (8) omits the 𝛽 𝑡 term, which we fit to theN-body data by performing a linear regression 𝒗 𝑡 = 𝛽 𝑡 𝒗 pred , (9)where 𝒗 𝑡 is the measured N-body velocity of a tracer of the velocityfield which, in principle, need not be the tracer as used for the den-sity field). This is a cross-correlation between two samples, and wewill refer to it generically as a velocity – density cross-correlation,although the actual comparison is made between observed and pre-dicted peculiar velocities.The procedure used to obtain the density tracer field is thesame regardless of whether the tracer data are particles, DM haloesor galaxies, as each simulation provides Cartesian positions andvelocities for all tracers. The latter two also provide additional in-formation such as halo mass, stellar mass and luminosity measure-ments in different filter bands. In this paper, we will construct thedensity fluctuation field using all of the aforementioned tracers. Inthe following section however, we will focus solely on using theparticles for the density field.Likewise, one also has a choice of which velocity tracers to usein the comparison. Using the particles as the velocity tracers is themost straightforward, so we will begin with this case in Section 5.1below. Generally, however observers usually combine the peculiarvelocities from multiple galaxies in the same group or cluster to ob-tain the mean peculiar velocity of the group. The 𝑁 -body analoguefor the group-averaged peculiar velocity is the peculiar velocity ofthe host (or central) halo.In order to calculate the density field, the density tracers areplaced on a 3D cubic grid at the nearest grid point and the contri-bution of all tracers at the same grid point is summed.In the case where grid spacing is non-negligible, the griddingacts like smoothing, and it adds in quadrature with the applied toyield a total effective smoothing (Boruah et al. 2020), 𝜎 = 𝜎 + 𝜎 , (10)where 𝜎 = Δ where Δ grid is the grid spacing in ℎ − Mpc.Using a fine grid allows us to preserve some of the detail andignore the effects of grid size smoothing. All comparisons in thispaper are done assuming a grid spacing of 0.36 and 0.98 h − Mpcfor Bolshoi and MDPL2 respectively, which produces a negligibleeffect on the effective smoothing scales that will be used in thispaper.We choose to smooth the density fluctuation field using a Gaus-sian smoothing kernel, which, in configuration and Fourier space,are given by 𝑊 ( 𝑟 ) = √ 𝜋𝑅 𝐺 exp (cid:32) − 𝑟 𝑅 𝐺 (cid:33) (11) 𝑊 ( 𝑘 ) = exp (cid:32) − 𝑘 𝑅 𝐺 (cid:33) , (12)respectively. We then calculate the peculiar velocities on the 3Dgrid using equation (4). To predict the peculiar velocity of a given velocity tracer, we linearly interpolate from the velocity grid to thetracer’s location.As discussed above, we then regress the tracer velocity againstits predicted velocity, with the intercept fixed to zero, and where thefitted slope gives an estimate for 𝛽 𝑡 . The fits also yield the rms (1D)velocity scatter, 𝜎 𝑣 , of the difference between the 𝑁 -body velocitiesand the linear theory predictions at the best fit 𝛽 𝑡 . As discussed above, the density tracer field must be smoothed inorder to predict the peculiar velocities in linear perturbation theory.In this section, we assess the impact of different Gaussian smoothinglengths on the recovered value of 𝛽 𝑡 , and the 1D 𝜎 𝑣 around the lineartheory predictions. We focus on smoothing scales between 1-6 ℎ − Mpc.Specifically, we aim to confirm at what smoothing length theslope of linear regression is unbiased, which previous work hassuggested lies in the range 4 – 5 ℎ − Mpc (Carrick et al. 2015).We first consider particles as tracers of the density field, and thenconsider haloes.
We now consider the case in which we predict the peculiar velocitiesof particles using the particle density fluctuation field. The valuesof the best-fit slope divided by the value of 𝑓 in the simulation,i.e. 𝛽 / 𝑓 , as a function of smoothing scale, are shown in Figure 1.This quantity should be unity if the method is unbiased. We findthat predictions are unbiased for a Gaussian smoothing kernel 𝑅 𝐺 between 4 and 5 ℎ − Mpc. Of the two simulations, MDPL2 should bemore accurate since its minimum wavenumber 𝑘 is 2 𝜋 × − ℎ / Mpc,and whereas that of Bolshoi is a factor 4 larger due to its smaller boxlength. Thus Bolshoi fails to capture the long-wavelength modes thatgenerate a significant fraction of the rms peculiar velocity. The 1D 𝜎 𝑣 around the best fit slope, however, is minimized at a smoothinglength that is 1-1.5 ℎ − Mpc smaller. The 𝜎 𝑣 of particle peculiarvelocities is high: between 225 and 275 km/s. This is because theparticle’s peculiar velocity includes its motion within the halo aswell as the peculiar velocity of the halo itself. Only the latter iswell-predicted by linear perturbation theory.If we predict the peculiar velocities of host (or central) haloes(i.e. excluding subhaloes) using the particle density field, we findsimilar results, with unbiased results and minimum 𝜎 𝑣 occurringfor Gaussian smoothing lengths that are ∼ . ℎ − Mpc smaller thanwhen using particles as velocity tracers. The 𝜎 𝑣 as a function ofsmoothing length is quite flat; it is not much higher at the fidu-cial 4 ℎ − Mpc. However the amplitude of 𝜎 𝑣 for host haloes issignificantly lower ( ∼
150 km/s) than for particles ( ∼
250 km/s),as expected given that the velocities of the particles include theirmotion with respect to the host halo.Finally, it is interesting to explore the question of whetherhaloes of different masses have different velocities relative to thepredictions of linear theory, one scenario of “velocity bias.” Wefind that imposing a minimum mass threshold on the haloes usedto sample the velocity field of 10 𝑀 (cid:12) has little effect on themeasured 𝛽 / 𝑓 for smoothing lengths greater that 1.5 ℎ − Mpc forMDPL2 and 3 ℎ − Mpc for Bolshoi, as shown in Fig. 2. This is alsotrue for a minimum mass of 10 𝑀 (cid:12) . A small velocity bias appearsonly when considering cluster mass haloes ( > 𝑀 (cid:12) ) as peculiar MNRAS , 1–10 (2020)
Hollinger and Hudson S l o p e o f L i n e a r R e g r e ss i o n ( / f ) simulation : tracerMDPL2 : particlesBolshoi : particlesMDPL2 : haloesBolshoi : haloes1 2 3 4 5 6 Smoothing Length h Mpc v [ k m / s ] Figure 1.
Top Panel: Slope of linear regression for the predicted and N-Bodyparticle and central halo velocities (excluding subhaloes) for MDPL2 andBolshoi simulations. Because we correct for the 𝑓 parameter, a value of 𝛽 / 𝑓 of unity indicates for the smoothing length for which the velocity–densitycross-correlation is unbiased. The predicted velocities were calculated as-suming an underlying particle density field and linear theory for the Gaussiansmoothing radius shown on the horizontal axis . Bottom Panel: The scatterbetween the measured 𝑁 -body and predicted peculiar velocity associatedwith each of the linear regression slopes. A circle has been placed at thesmoothing length where the standard deviation was minimized. velocity tracers. For clusters, there is also an increase in the 1D 𝜎 𝑣 with respect to the linear theory predictions. We now consider a scenario that is closer to the observational one,where the density field is provided by DM haloes, weighted by theirmass. Whereas the particle density field is unbiased, this field willbe biased. Therefore, we no longer expect 𝛽 h / 𝑓 =
1. As discussedabove, in the linear regime, this bias can be calculated by measuringthe rms density fluctuations of the halo-weighted density field in 8 ℎ − Mpc spheres. The halo and particle 𝜎 measurements are relatedby 𝑏 h = 𝜎 𝜎 (13)which is the same 𝑏 in the 𝛽 term defined previously. We measure 𝛽 / 𝑓 from the slope of linear regression as before, but multiplyingthis by the correction factor 𝑏 h . This should return a value of unity if S l o p e o f L i n e a r R e g r e ss i o n ( / f ) MDPL2all haloesM h >10 h MM h >10 h MM h >10 h M1 2 3 4 5 6
Smoothing Length h Mpc v [ k m / s ] Figure 2.
As in Fig. 1 where again the predicted peculiar velocities arebased on the particle density field, but here these are compared to the 𝑁 -body peculiar velocities of haloes (but excluding subhaloes). For clarity thesolid line labeled ’all halos’ represents all simulation objects classified ascentrals, while the dashed lines represent centrals with masses greater thanvarious minimum masses as indicated in the legend. the method unbiased. This parameter combination will be referredto as 𝛽 ∗ = 𝛽𝑓 𝜎 𝜎 .To measure 𝜎 we placing non-overlapping spheres of 8 ℎ − Mpc covering the entire the density field, and measure thestandard deviation in the values. Doing this, we find values of 𝜎 = . ± .
04 and 1 . ± .
02 for the halo masses for Bol-shoi and MDPL2 respectively, much larger than 𝜎 , which is themeasured (non-linear) 𝜎 of the particle density field.The fitted values of 𝛽 h / 𝑓 , however, are lower than unity, asexpected. The results for 𝛽 ∗ are plotted in Figure 3, showing thatthe mass-weighted determination is nearly unbiased: at the fiducial 𝑅 𝐺 = ℎ − Mpc, the MDPL2 field has 𝛽 ∗ = .
05, whereas forBolshoi 𝛽 ∗ = .
02. This suggests that linear biasing correctionworks well, even for fields with 𝑏 h ∼ .
7. The 1D 𝜎 𝑣 ( ∼ It is not obvious why a Gaussian smoothing, with 𝑅 𝐺 ∼ − ℎ − Mpc, should be the “correct” smoothing length. We can gain someinsight by considering the problem in Fourier space, specifically
MNRAS000
MNRAS000 , 1–10 (2020) esting Accuracy of Peculiar Velocity Comparisons ( / f ) * ( , h / , M ) Minimum Halo MassMDPL2Bolshoiall haloesM h >10 h M1 2 3 4 5 6
Smoothing Length h Mpc S c a tt e r [ k m / s ] Figure 3.
Similar to Fig. 2, except here the mass-weighted halo density fieldis used to predict peculiar velocities of haloes. The vertical axis in the toppanel is now 𝛽 / 𝑓 ( 𝜎 ,ℎ / 𝜎 ,𝑚 ) . equation 4 which relates the velocity modes to the density modesin linear perturbation theory. This should be exact on large scales,but we expect it to break down at high 𝑘 .For simplicity, first assume that all Fourier modes are Gaus-sian and independent and that at a given 𝑘 the joint distribution ofvelocity and density mode amplitudes is given by a bivariate Gaus-sian. We expect the correlation coefficient, 𝜌 , to approach 1 on largescales and zero at high 𝑘 . Generally if one has a bivariate normaldistribution of two variables 𝑥 and 𝑦 , and one regresses 𝑦 on 𝑥 theslope of the relation is 𝜌 𝜎 𝑦 𝜎 𝑥 , or since the covariance 𝐶 𝑥𝑦 = 𝜌𝜎 𝑥 𝜎 𝑦 ,the slope can also be written 𝐶 𝑥𝑦 / 𝜎 𝑥 . Therefore, the predicted ve-locity mode amplitude is given by the density mode amplitude timesthe slope 𝐶 𝑣 𝛿 / 𝜎 𝛿 .The quantities 𝐶 𝑣 𝛿 and 𝜎 𝛿 can be measured in 𝑁 -body simu-lations as a function of wavenumber 𝑘 . The latter is just the powerspectrum of matter density fluctuations, 𝑃 𝛿 𝛿 , but the covarianceis less well studied. Zheng et al. (2013) have measured the cross-spectrum of the closely related quantity 𝜃 = ∇ · 𝑣 with the density 𝛿 , and its ratio with the density power spectrum. They define thenormalized window function (cid:101) 𝑊 𝑘 = 𝑓 𝑃 𝜃 𝛿 𝑃 𝛿 𝛿 , (14)which has the property that it asymptotes to 1 at low 𝑘 and goesto zero at higher 𝑘 . This function is plotted in Figure 4. Also over- k/h Mpc W ( k ) e ( kR G ) /2 , R G = 2.88 e ( kR G ) /2 , R G = 4.0 Figure 4.
The normalized window function from Zheng et al. (2013), cal-culated from the J1200 simulation, based on the correlation of the velocityand density fields (see text for details) is shown, with Window functionsfor Gaussian smoothing kernels with scales of 2.85 and 4.0 ℎ − Mpc. Bothfunctions exhibit the expected characteristics of being unity as 𝑘 → 𝑘 → ∞ . (cid:102) 𝑊 𝑘 and the Gaussian kernel of 𝑅 𝐺 = . ℎ − Mpc cross0.5 near 𝑘 = .
39 and 0 . ℎ / Mpc, respectively. plotted for comparison are two Gaussian smoothing filters with 𝑅 𝐺 = .
88 and 4.0 ℎ − Mpc. The estimated 𝛽 / 𝑓 from the par-ticle weighted density field smoothed using either (cid:101) 𝑊 𝑘 as a 𝑘 -space smoothing function or with a Gaussian kernel of 𝑅 𝐺 = . ℎ − Mpc both produce values of 0.88. There is also no significantdifference in 𝜎 𝑣 .The Gaussian smoothing with 𝑅 𝐺 = . ℎ − Mpc is a bettermatch to (cid:101) 𝑊 𝑘 at low 𝑘 .In reality, particularly at high 𝑘 , the assumption of bivariateGaussianity will no longer be correct, so the above simple argumentwill break down. In principle, with detailed knowledge of the cor-relations, it should be possible to design the optimal 𝑘 -space filter.For our purposes in this paper, we retain the simplicity of Gaussiansmoothing, and adopt a Gaussian smoothing filter with 𝑅 𝐺 = ℎ − Mpc for consistency with previous work.
In the previous sections, the masses of the haloes were assumed to beknown exactly, with no uncertainty in the measurements. Howeverin real survey data there is a some uncertainty in the true mass ofany given halo, depending on the method used to estimate it, andthis may be in the range 0.1 – 0.2 dex (see Section 7). In this section,we explore how scattering the halo mass impacts the predictions of 𝛽 / 𝑓 and 𝜎 ,ℎ . This is accomplished by varying the halo massesby a lognormal Gaussian random variable with standard deviationscorresponding to 0.1 and 0.2 dex, and calculating 𝜎 and 𝛽 for each MNRAS , 1–10 (2020)
Hollinger and Hudson / f * ( , h / , M ) Minimum Halo MassBolshoi:all haloesMDPL2:all haloesBolshoi:10 h MMDPL2:10 h M0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 scatter [dex] v [ k m / s ] Figure 5.
Effect of additional scatter to the mass of haloes as a function ofthe logarithmic scatter (0, 0.1 and 0.2 dex). For 0 dex this is the measuredvalue of the simulation, for 0.1 and 0.2 dex a random value correspondingto a lognormal Gaussian with a 𝜎 of ∼ ± 𝜎 range (shaded band) of these measurements are shown for 𝛽 / 𝑓 ( 𝜎 ,ℎ / 𝜎 ,𝑚 ) at a smoothing radius of 4 ℎ − Mpc. Bottom panel :scatter around the best fit of 𝜎 𝑣 for the same smoothing. realization. A total of 500 realizations were performed on both theBolshoi and MDPL2 halo catalogues.We find that, from realization to realization, both the measured 𝛽 and 𝜎 deviate from the no-scatter values, and these deviationsare anti-correlated. This can be understood as follows. In Section5.3, we discussed how the slope ( 𝛽 ) could be expressed as the ratioof the covariance between density and velocity, and the densitypower spectrum, 𝐶 𝑣 𝛿 / 𝜎 𝛿 𝛿 . If noise is added to the density field,then the denominator increases, but the covariance is unaffected.Consequently, the fitted 𝛽 drops. On the other hand, 𝜎 increasesbecause of the additional noise. This anti-correlation leads to somecancellation of the additional noise from realization to realizationin the product 𝛽 ∗ = 𝛽 h 𝜎 .The net result of increasing the halo mass scatter and its impacton 𝛽 ∗ are shown in the top panel of Figure 5, with a scatter of zerodex corresponding to the case where the halo masses are preciselymeasured. We find that when introducing a lognormal scatter,the value of 𝛽 ∗ increases only slightly with increasing scatterfor both Bolshoi and MDPL2. We note however that the standarddeviation of 𝛽 ∗ is higher for Bolshoi due to the smaller volumeof the simulation, which leads to greater variation from realization l o g ( M * [ h M ]) SAGESAG10 11 12 13 14 15 log(M h [ h M ]]) l o g ( M * / M * ) Figure 6.
Top Panel: The stellar-to-halo mass relation (SHMR) for the totalstellar mass of a halo (including stellar mass in subhaloes) as a functionof total halo mass. The dark curves show the means of the SHMR and thelighter bands indicate their measured standard deviation. Both SAG andSAGE semi-analytic models are shown. Bottom Panel: the stellar to halomass ratio as a function of halo mass. to realization in 𝛽 h , 𝜎 h and their product 𝛽 ∗ . For MDPL2, the neteffect of scatter in the 𝑁 independent halo masses is reduced by afactor 1 /√ 𝑁 due to the larger volume containing a larger number 𝑁 of haloes. This highlights the importance of the volume of thesample, a topic we shall discuss in greater detail in Section 8.The bottom panel shows how the 1D 𝜎 𝑣 is impacted by intro-ducing stochasticity in the halo masses. We find that the 0.1 dex casegenerates a negligible change in the measured 𝜎 𝑣 when comparedto the original value. For halo mass uncertainties of 0.2 dex, 𝜎 𝑣 increases, but the effect is still small, corresponding to ∼ Galaxies trace the underlying dark matter distribution on largescales, but observable quantities, such as a galaxy’s luminosityand stellar mass, do not necessarily have an exact relationship withthe mass of the halo in which a galaxy lies. To explore how usingthese quantities as a proxy for mass density impact 𝛽 / 𝑓 and 𝜎 , weuse two galaxy semi-analytic galaxy formation models available forthe MDPL2 simulation SAG (Cora 2006) and SAGE (Croton et al.2016).Many DM haloes host more than one galaxy, which are typ- MNRAS000
Top Panel: The stellar-to-halo mass relation (SHMR) for the totalstellar mass of a halo (including stellar mass in subhaloes) as a functionof total halo mass. The dark curves show the means of the SHMR and thelighter bands indicate their measured standard deviation. Both SAG andSAGE semi-analytic models are shown. Bottom Panel: the stellar to halomass ratio as a function of halo mass. to realization in 𝛽 h , 𝜎 h and their product 𝛽 ∗ . For MDPL2, the neteffect of scatter in the 𝑁 independent halo masses is reduced by afactor 1 /√ 𝑁 due to the larger volume containing a larger number 𝑁 of haloes. This highlights the importance of the volume of thesample, a topic we shall discuss in greater detail in Section 8.The bottom panel shows how the 1D 𝜎 𝑣 is impacted by intro-ducing stochasticity in the halo masses. We find that the 0.1 dex casegenerates a negligible change in the measured 𝜎 𝑣 when comparedto the original value. For halo mass uncertainties of 0.2 dex, 𝜎 𝑣 increases, but the effect is still small, corresponding to ∼ Galaxies trace the underlying dark matter distribution on largescales, but observable quantities, such as a galaxy’s luminosityand stellar mass, do not necessarily have an exact relationship withthe mass of the halo in which a galaxy lies. To explore how usingthese quantities as a proxy for mass density impact 𝛽 / 𝑓 and 𝜎 , weuse two galaxy semi-analytic galaxy formation models available forthe MDPL2 simulation SAG (Cora 2006) and SAGE (Croton et al.2016).Many DM haloes host more than one galaxy, which are typ- MNRAS000 , 1–10 (2020) esting Accuracy of Peculiar Velocity Comparisons / f * ( , g / , M ) SAGESAGall haloesM h >10 h M1 2 3 4 5 6
Smoothing Length h Mpc v [ k m / s ] Figure 7.
Similar to Fig. 3, except here a fitted SHMR is used to predictthe halo masses from stellar masses, and then constructing a mass-weightedhalo density field to predict the peculiar velocities of central galaxies. ically divided into centrals (the galaxy identified with the main orhost DM halo) and satellites (associated with DM subhaloes). Whenpredicting host halo mass, one can use only the stellar mass or lu-minosity of the central galaxy, or one can use the total stellar mass(or total luminosity) of all galaxies. We adopt the latter approachhere when calculating the density fields. As before, velocity com-parisons will be done solely on galaxies classified as centrals. Forconsistency with the previous work done in this paper, any cutsimposed on the data will be done using the mass of the host halo.The stellar to halo mass relation (SHMR) is different for boththe SAG and SAGE semi-analytic models in MDPL2, see 6. SAGhas a tighter SHMR with less scatter in stellar mass at a given halomass: for haloes with masses between 10 and 10 ℎ − 𝑀 (cid:12) , ithas an average scatter of 0.15 dex compared to 0.39 dex in SAGE.The average SHMR is comparable for both SAG and SAGE for halomasses greater than the characteristic pivot point at ∼ ℎ − 𝑀 (cid:12) . To construct a proxy mass density field, we can obtain proxy halomasses using the SHMR. The K-band luminosity is a good proxyfor stellar mass, however the galaxy luminosities in SAG are limitedto the 𝑢, 𝑔, 𝑟, 𝑖, 𝑧 -bands. So in this section, we use the stellar massesprovided for both SAG and SAGE. We fit the SHMR profiles of bothSAG and SAGE semi-analytic models with a broken power law, and use these to convert the stellar masses into halo masses. Then thesame analysis described in a previous section is performed. Thisprocedure will remove the nonlinearity of the SHMR relation, butcannot remove the scatter in it, since the mean SHMR is used foreach halo.From Figure 7, we find for the SAGE galaxies, 𝛽 ∗ is unbiasedfor a Gaussian smoothing of 4.1 ℎ − Mpc but generates the lowest 𝜎 𝑣 at 𝑅 𝐺 = ℎ − Mpc. For SAG these are at smoothing lengths of3.6 and 3.5 ℎ − Mpc, respectively. Comparing these results to Fig. 3,we find that at 𝑅 𝐺 = 4 ℎ − Mpc, SAG produces results that are quitesimilar (within a few percent) to what was found using MDPL2’shalo masses directly. SAGE, however, generates values of 𝛽 ∗ thatare 5% smaller at the same smoothing length than the haloes. Wefind that the generated SAG haloes produce similar estimates forthe 1D scatter in velocity predictions, contrarily the SAGE haloesgenerate ∼
20 km/s less scatter.Comparing these results to those in the next section we notethat while the halo mass does tend to be unbiased at smaller 𝑅 𝐺 theconversion of stellar to halo mass introduces 10-20 km/s of addi-tional scatter than using galaxy observables to predicted velocities.We attribute the results to the differences in the models’ stellarto halo mass ratios as a function of halo mass. SAG has a flatterstellar to halo mass ratio than SAGE, and therefore it is close tothe simple halo mass weighted case. SAGE has a steeper ratio athigh halo masses, hence SAGE puts less weight on massive clustersleading to estimates of 𝛽 ∗ slightly less than unity. A density field can also be constructed without approximating themass of individual haloes via the SHMR. In this section, we in-vestigate how weighting directly by the galaxy observables impactsthe cosmological estimates. Such a procedure is closer to what wasdone by Carrick et al. (2015), who used the 𝐾 -band luminosity-weighted density field from the 2M++ catalogue. Figure 8 showsthe summary of results discussed below.We find that, for the SAG model, weighting the density field bystellar mass produces results that yield a higher 𝛽 ∗ than luminosityweighting. We can weight the density field using luminosities in anyof the five bands provided by SAG. In the remainder of this paper wefocus on the 𝑟 -band luminosity. We note however that the longestwavelength 𝑧 -band produces the highest values of 𝛽 / 𝑓 howeverafter applying the correction factor ( 𝜎 / 𝜎 ) there is virtuallyno difference between bands for 𝛽 ∗ , provided a low minimum massthreshold.If the minimum mass threshold of the haloes is low, we findthat weighting using the stellar masses provided by SAG closelyresembles the case where the density field is halo-mass weightedand is unbiased at an 𝑅 𝐺 ∼ . ℎ − Mpc. The SAGE stellar massand SAG 𝑟 -band predict 𝛽 ∗ values that are comparable at 𝑅 𝐺 > ℎ − Mpc, but only produce unbiased estimated of 𝛽 ∗ for a Gaussiansmoothing radius of ∼ ℎ − Mpc.In the case where a minimum mass threshold of 10 ℎ − 𝑀 (cid:12) isapplied to the data, we find that 𝛽 ∗ for the three cases are comparablefor 3 ℎ − Mpc (cid:46) 𝑅 𝐺 (cid:46) ℎ − Mpc. With all the cases estimatingunbiased 𝛽 ∗ values at 3.5-3.7 ℎ − Mpc. The 𝜎 𝑣 is also comparablefor these cases with the SAGE stellar mass producing only ∼ 𝑟 -band. MNRAS , 1–10 (2020)
Hollinger and Hudson
Table 1.
Summary of 𝛽 ∗ values taken for 𝑅 𝐺 = ℎ − Mpc for the various MDPL2 tracers which weight the density field and from which the peculiar velocitiesare compared. Density Tracer all M 𝑡 𝑀 𝑡 > ℎ − 𝑀 (cid:12) Particles 0.97 – Figure 1Haloes 1.05 1.08 Figure 3SAGE: SHMR 1.00 1.02 Figure 7SAG: SHMR 1.04 1.06 Figure 7SAGE: stellar mass 0.98 1.01 Figure 8SAG: stellar mass 1.01 1.03 Figure 8SAG: 𝑟 -band luminosity 0.97 1.02 Figure 8 / f * ( , g / , M ) SAG: rbandSAG: stellar massSAGE: stellar massall haloesM h >10 h M1 2 3 4 5 6
Smoothing Length h Mpc v [ k m / s ] Figure 8.
Similar to Fig. 7 except now the density field is constructed usinggalaxy observables (i.e. weighted by stellar mass or luminosity) to predictand measure centrals’ peculiar velocities.
In real data sets both the data used to obtain the density field, basedon redshift surveys, and the peculiar velocity samples, cover limitedvolumes. This has multiple consequences. Some of these are relatedto sample variance: the local volume will have local values of themean density and of 𝜎 that are different from the global values.Moreover, the predicted peculiar velocities will be less accurate forgalaxies close to the edge of the density field than for those far fromthe edges, say at the center of the volume.In this section, we investigate these finite volume and edgeeffects by limiting the data, to spheres of radii with 𝑅 max = Figure 9.
To demonstrate the dependence of our calculated parameters onthe survey volume, we restrict simulation data to (non-overlapping) spheresof radius ranging from 50 to 300 ℎ − Mpc. The components of 𝛽 ∗ arecalculated for each test sample assuming a Gaussian smoothing kernel of 4 ℎ − Mpc. The upper panel and lower panel show the results for 𝛽 ∗ and in 𝜎 measurements, respectively, as a function of sphere radius. In both panels,the light and dark coloured bands represent the ± 𝜎 standard deviationfrom sphere to sphere and the standard error in the mean of the test cases(respectively) for a given sample spherical radius, with the mean beingshown by the solid colour lines. The colour horizontal dashed lines showthe values from the full simulation box. ℎ − Mpc, and its impact on cosmologicalestimates. In particular, we are interested in simulations of data setswith 𝑅 max ∼ ℎ − Mpc, which is comparable to the survey sizeof 2M++, the depth of which varies from 125 ℎ − Mpc to 200 ℎ − Mpc, with an effective spherically-averaged radius of 175 ℎ − Mpc.
MNRAS000
MNRAS000 , 1–10 (2020) esting Accuracy of Peculiar Velocity Comparisons / f * ( , g / , M ) SAGE: stellar massSAG: stellar massSAG: r-band1 2 3 4 5 6
Smoothing Length h Mpc v [ k m / s ] Figure 10.
Similar to Fig. 8 expect we now use non-overlapping sphere cutsof radius 150 ℎ − Mpc and the components of 𝛽 ∗ are calculated for a rangeof Gaussian kernel lengths between 1-6 ℎ − Mpc. The bands have the samemeaning as in Fig. 9.
Given the size of the MDPL2 simulation, it is possible to generatemultiple independent (non-overlapping) finite volume realizations.For each sphere, we ignore any galaxy that exists outside of thesphere, and assign 𝛿 = 𝛿 using the mean density withinthe sphere instead of the simulation box average. For each sphere,we then calculate 𝛽 / 𝑓 and 𝜎 using only that sphere’s galaxies,although 𝜎 ,𝑚 , which appears in the denominator of 𝛽 ∗ , continuesto be calculated from the full simulation. As in real analyses, toaccount for the missing contribution from structures beyond thesphere’s edge, in addition to fitting the 𝛽 / 𝑓 , we also fit a residualbulk flow term.Figure 9 demonstrates how 𝛽 ∗ and 𝜎 depend on 𝑅 max forthe case with no mass cut on the galaxy catalogue. As expected, asthe size of the sphere increases, results converge to the values foundfor the full simulation for 𝑅 max ≥ ℎ − Mpc.There is a significant amount of scatter in the measured 𝜎 values for small 𝑅 max , varying by as much as 10% of the meanvalue, and with a mean value that tends to be between 4-8 % lowerthan that of the full simulation. The bias in the mean may be due torenormalizing the density fluctuation, 𝛿 = ( 𝜌 − ¯ 𝜌 )/ ¯ 𝜌 with the local ¯ 𝜌 measured on the scale of the sphere, 𝑅 max , thus effectively filteringout the contribution to 𝜎 from large-scale waves. The scatter arises from sample variance effects, which is a combination of cosmicvariance in the underlying DM structures and stochasticity in theSHMR. This decreases as 𝑅 max increases.Likewise, the locally-measured value of 𝛽 / 𝑓 also varies fromsphere to sphere. Again, there are several reasons for this. First, thelocal ¯ 𝜌 averaged over the sphere will differ from the true average.This is important because, as noted above, it leads to a renormal-ization of 𝛿 . Second, there is also stochasticity in the SHMR whichmay affect the locally-measured 𝛽 / 𝑓 . Of these two, the first ef-fect is the dominant one: we find that the total scatter in 𝛽 / 𝑓 for 𝑅 max = ℎ − Mpc at 𝑅 𝐺 = ℎ − Mpc is typically 7.0% and ismostly due to cosmic variance in ¯ 𝜌 . The scatter in 𝛽 / 𝑓 dominatesthe scatter in 𝜎 (2.5%) when they are combined in 𝛽 ∗ .Figure 10 demonstrates how the cut samples ( 𝑅 max = ℎ − Mpc) depend on the Gaussian smoothing radius 𝑅 𝐺 . In particular,when compared to Figure 8 we find that the average 𝛽 ∗ for thesevolume-limited realizations are unbiased between 3.75-4.75 ℎ − Mpc. The value of 𝛽 ∗ at 𝑅 𝐺 = ℎ − Mpc ranges between 0.98-1.01.More importantly, the standard deviation in 𝛽 ∗ from sphere to sphereis 0.077 at 𝑅 max = ℎ − Mpc, although this declines significantlyto 0.032 at 200 ℎ − Mpc and 0.020 at 300 ℎ − Mpc. This suggeststhat sample variance effects are the dominant uncertainty in 𝑓 𝜎 for current data sets.The 1D velocity scatter as a function of 𝑅 𝐺 has a similarshape as found previously, but is 40-50 km/s larger than for thefull simulation box. This increase in the scatter is primarily dueto the degradation of the predicted velocities as one approaches 𝑅 max . When we compare the 1D 𝜎 𝑣 for the subsample of galaxieslocated within the inner half of the sphere’s volume with those inthe outer half at an 𝑅 𝐺 = ℎ − Mpc, we find scatters of ∼
165 km/sand ∼
210 km/s, respectively. This increase is attributed to the factthat the velocity tracers near the outer edge have poorer predictedpeculiar velocities because of unaccounted-for structures outside ofthe survey limits.
The key results of this paper are as follows: • The velocities of DM haloes are well predicted by linear theoryfrom the true density field with a Gaussian smoothing 𝑅 𝐺 = ℎ − Mpc with a velocity scatter of 154 km/s. This is in agreement withAppendix A of Carrick et al. (2015). • This can be understood because, in Fourier space, a Gaussianfilter with 𝑅 𝐺 = ℎ − Mpc is a good match to the cross-correlation function of the density and velocity fields. • The accuracy and precision of the linear theory predictions donot depend on the mass of the velocity tracer; there is no "velocitybias", except for clusters with 𝑀 h > ℎ − 𝑀 (cid:12) . • If DM haloes are used as tracers of the density field, and onecalculates 𝜎 of the halo-mass-weighted density field, then 𝛽 h 𝜎 is a good estimator of 𝑓 𝜎 . • If noise is added to the DM halo masses, then 𝛽 ∗ is biased highby only a percent, for a 0.1 dex noise level. • When galaxy luminosity or stellar mass are used for the densityfield, the values of 𝛽 ∗ indicate that the method is unbiased to within5%, depending on the semi-analytic galaxy formation model. • When the density field is restricted to a finite volume, there isadditional uncertainty due to cosmic variance, at the level of 7% fora 150 ℎ − Mpc sphere.
MNRAS , 1–10 (2020) Hollinger and Hudson
The results for 𝛽 ∗ calculated using the same tracers for thevelocity and density field, are summarised in Table 1. Overall wefind that the method has ∼
5% systematic uncertainties. This canbe improved with semi-analytic galaxy formation models that morecarefully match the real SHMR and its scatter. There is also uncer-tainty due to finite volumes and cosmic variance.Previous work has neglected the additional scatter due to thecoupled effects of stochasticity in the galaxy-mass relation and thecosmic variance effect of finite volumes. For example, the 2M++catalogue (Lavaux & Hudson 2011) has an effective volume of175 ℎ − Mpc. The SAGE model, weighted by stellar mass, showsthat, for 150 ℎ − Mpc spheres, the uncertainty on 𝛽 g 𝜎 ,𝑔 = 𝑓 𝜎 is 7.7%, while the expected (interpolated) value for a survey of2M++’s size is 5.2%. This is slightly higher than the uncertaintyestimated by Carrick et al. (2015), Boruah et al. (2020) and Saidet al. (2020) who adopted a 4% sampling variance uncertainty, plusobservational errors in 𝛽 due to uncertainties in peculiar velocitymeasurements, which are subdominant. For precise quantificationof the biases and systematic uncertainties in 𝑓 𝜎 derived from aspecific survey, e.g. 2M++, the best approach to minimizing thesystematic errors will be to create mock catalogues that mimic thegeometry and selection of that particular survey.The cosmic variance uncertainty can be reduced in the futurewith deeper, all-sky redshift surveys. For example, a survey extend-ing to a redshift of 0.2 (600 ℎ − Mpc) would have an uncertainty ofonly 0.4% in the mean mass density and hence 0.5 – 0.6% in theluminosity (or stellar mass) density. In the North, the DESI BrightGalaxy Survey (DESI Collaboration et al. 2016), will observe 10million nearby galaxies. In the South, 4MOST (de Jong et al. 2019)has the capability to survey large volumes in the nearby Universe.The future looks bright.
ACKNOWLEDGEMENTS
DATA AVAILABILITY STATEMENT
The data underlying this article are publicly available from theCOSMOSIM database , with theirrespective publications cited in section 3.
REFERENCES
Adams C., Blake C., 2020, MNRAS, Asgari M., et al., 2019, KiDS+VIKING-450 and DES-Y1 com-bined: Mitigating baryon feedback uncertainty with COSEBIs( arXiv:1910.05336 )Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013a, ApJ, 762, 109Behroozi P. S., Wechsler R. H., Wu H.-Y., Busha M. T., Klypin A. A.,Primack J. R., 2013b, ApJ, 763, 18Berlind A. A., Narayanan V. K., Weinberg D. H., 2000, ApJ, 537, 537Boruah S. S., Hudson M. J., Lavaux G., 2020, MNRAS, 498, 2703Carrick J., Turnbull S. J., Lavaux G., Hudson M. J., 2015, Monthly Noticesof the Royal Astronomical Society, 450, 317Cora S. A., 2006, MNRAS, 368, 1540Cora S. A., et al., 2018, Monthly Notices of the Royal Astronomical Society,479, 2–24Croton D. J., et al., 2016, The Astrophysical Journal Supplement Series,222, 22DESI Collaboration et al., 2016, arXiv e-prints, p. arXiv:1611.00036Davis M., Nusser A., Masters K. L., Springob C., Huchra J. P., Lemson G.,2011, MNRAS, 413, 2906Dekel A., 1994, ARA&A, 32, 371Dupuy A., Courtois H. M., Kubik B., 2019, MNRAS, 486, 440Hudson M. J., Turnbull S. J., 2012, The Astrophysical Journal Letters, 751,L30Huterer D., Shafer D. L., Scolnic D. M., Schmidt F., 2017, J. CosmologyAstropart. Phys., 2017, 015Klypin A. A., Trujillo-Gomez S., Primack J., 2011, ApJ, 740, 102Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457,4340Knebe A., et al., 2017, Monthly Notices of the Royal Astronomical Society,474, 5206–5231Kravtsov A. V., Klypin A. A., Khokhlov A. M., 1997, The AstrophysicalJournal Supplement Series, 111, 73–94Lavaux G., Hudson M. J., 2011, Mon. Not. R. Astron. Soc., 416, 2840Nusser A., Davis M., Branchini E., 2014, ApJ, 788, 157Peebles P. J. E., 1993, Principles of physical cosmology. Princeton UniversityPressPike R., Hudson M. J., 2005, The Astrophysical Journal, 635, 11Planck Collaboration 2018, arXiv e-prints, p. arXiv:1807.06209Qin F., Howlett C., Staveley-Smith L., 2019, MNRAS, 487, 5235Said K., Colless M., Magoulas C., Lucey J. R., Hudson M. J., 2020, MNRAS,497, 1275Strauss M. A., Willick J. A., 1995, Phys. Rep., 261, 271Turnbull S. J., Hudson M. J., Feldman H. A., Hicken M., Kirshner R. P.,Watkins R., 2012, MNRAS, 420, 447Willick J. A., Courteau S., Faber S., Burstein D., Dekel A., Strauss M. A.,1997, The Astrophysical Journal Supplement Series, 109, 333Zheng Y., Zhang P., Jing Y., Lin W., Pan J., 2013, Physical Review D, 88de Jong R. S., et al., 2019, The Messenger, 175, 3This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000