Cosmological Parameter Constraints from the SDSS Density and Momentum Power Spectra

We extract the galaxy density and momentum power spectra from a subset of early-type galaxies in the Sloan Digital Sky Survey (SDSS) DR7 main galaxy catalog. Using galaxy distance information inferred from the improved fundamental plane described in Yoon and Park, we reconstruct the peculiar velocities of the galaxies and generate number density and density-weighted velocity fields, from which we extract the galaxy density and momentum power spectra. We compare the measured values to the theoretical expectation of the same statistics, assuming an input Λ cold dark matter (CDM) model and using a third-order perturbative expansion. After validating our analysis pipeline with a series of mock data sets, we apply our methodology to the SDSS data and arrive at constraints fσ8=0.471−0.080+0.077 and b1σ8=0.920−0.070+0.070 at a mean redshift z¯=0.04 . Our result is consistent with the Planck cosmological best-fit parameters for the ΛCDM model. The momentum power spectrum is found to be strongly contaminated by small-scale velocity dispersion, which suppresses power by ∼(30%) on intermediate scales k ∼ 0.05 h Mpc−1.


INTRODUCTION
Extracting information from the spatial distribution of galaxies is a perpetual occupation within the cosmological community.From angular positions and redshifts, one can reconstruct the galaxy number density field and hence measure the galaxy N -point statistics in position and Fourier space.Typically the two-point correlation function and power spectrum are utilised, owing to their relative simplicity and considerable constraining power (Blake et al. 2011;Anderson et al. 2014;Oka et al. 2014;Okumura et al. 2016;de la Torre et al. 2017).Previously intractable higher point statistics are increasingly studied due to rapid advances in numerical statistical modelling and larger data volumes (Gil-Marín et al. 2017;Yankelevich & Porciani 2019;Philcox 2022).Alternative approaches to the standard N -point functions are motonari.tonegawa@apctp.orgalso progressively being explored (Pisani et al. 2015;Pan et al. 2020;Uhlemann et al. 2020;Villaescusa-Navarro et al. 2021;Appleby et al. 2022;Qin et al. 2023).
Galaxy properties are measured in redshift space, in which cosmological redshifts are contaminated by local peculiar velocities parallel to the line of sight.By modelling this redshift-distortion effect we can constrain not only the parameters governing the shape and amplitude of the power spectrum, but also the rate at which structures are forming.The reason is that peculiar velocities trace in-fall into gravitational potentials, a phenomenon that occurs on all scales (Kaiser 1987).Unfortunately, the effect of redshift space distortion is strongly degenerate with the overall amplitude of the power spectrum -at the level of linearized perturbations the two are exactly degenerate.However, if we complement the galaxy positional information by also measuring their velocities, then we can additionally extract the velocity power spectrum and break this degeneracy, placing simultaneous constraints on the amplitude of the galaxy power spectrum and the growth rate of structure.By modelling the quasi-linear scales using higher order perturbation theory, one can obtain additional information on galaxy bias and improved constraints on cosmological parameters (d' Amico et al. 2020).On small scales, stochastic velocities arising from bound structures dominate the cosmological signal, a phenomenon known as the Fingerof-God effect (Jackson 1972;Park et al. 1994;Fisher 1995;Juszkiewicz et al. 1998;Hikage & Yamamoto 2013;Tonegawa et al. 2020;Scoccimarro 2004;Jennings et al. 2011Jennings et al. , 2010;;Okumura & Jing 2010;Kwan et al. 2012a;White et al. 2014).
Distance measurements to galaxies are difficult to infer, as they require prerequisites for a physical scaling relation that can be used to convert dimensionless quantities such as redshift into distance.The complicated dynamics and evolution of galaxies means that there is no simple universal scale associated with their morphology, although certain subsets are known to be reasonably well approximated as virialised systems.For such subsets, we are able to directly measure the distance to recover the velocity field of the matter distribution.
The velocity field is problematic to reconstruct since it is expected to be non-zero even in empty spaces such as voids1 .To evade this problem, the density weighted velocity field -or momentum field -was first proposed as a more amenable statistic in Park (2000), since it naturally approaches zero in regions where the density is low.This statistic was used in Park & Park (2006) to place cosmological parameter constraints with some nascent large scale structure catalogs, but then the idea lay dormant for a time, before being picked up in a series of recent works (Howlett 2019;Qin et al. 2019).Significant improvements towards understanding the quasilinear ensemble average of the momentum field has been made in the intervening period (Vlah et al. 2012;Saito et al. 2014).
In this work we take the redshift and velocity information of a subset of SDSS galaxies measured in Yoon & Park (2020), extract the galaxy density and momentum power spectra and fit the standard cosmological model to the statistics, inferring a constraint on the galaxy power spectrum amplitude b 1 σ 8 and growth rate f σ 8 .While working towards this goal, we encounter a series of difficulties that we document in the following sections, pertaining to the convolution of the mask, residual parameter degeneracies and the effect of non-linear velocity dispersion on the shape of the momentum power spectrum.We highlight the assumptions and approximations required to arrive at our parameter constraints throughout.
In Section 2 we review the data used in our analysis, and the method that we use to extract the power spectra from the point distribution.We present the ensemble expectation values of the statistics, to which we compare our measurements in Section 3 including the effect of the mask.Some preliminary measurements from mock galaxy catalogs are presented in Section 4. Our main results can be found in Section 5, where we provide measurements of the power spectra from the data, and the resulting cosmological parameter constraints are in Section 5.2.We discuss our findings in Section 6 and review the assumptions made in arriving at our results.

DATA
A subset of the SDSS Data Release 7 (DR7; Abazajian et al. 2009) classified as early-type galaxies (ETGs) in the KIAS value-added galaxy catalog (VAGC; Choi et al. (2010)) is used to measure the galaxy density and momentum power spectra.It is based on the New York University Value-Added Galaxy Catalog (Blanton et al. 2005), where missing redshifts are supplemented from various other redshift catalogs to improve the completeness at higher redshift.The classification between early and late type is based on u − r colour, g − i colour gradient and inverse concentration index in the i-band (Park & Choi 2005;Choi et al. 2007).The result of the automated classification is corrected by visual inspection; we direct the reader to Park & Choi (2005) for details.We use a sub-sample obtained by applying redshift and magnitude cuts.ETGs are selected in the redshift range 0.025 ≤ z spec < 0.055, with the lower limit corresponding to a distance ∼ 75 h −1 Mpc to mitigate large peculiar velocity effects in the local Universe.A de Vaucouleur absolute magnitude cut of M dV r ≤ −19.5 is also applied.There are a total of 16, 283 galaxies in the sample, although we make further cuts below.In the top left panel of Figure 1 we present the angular distribution of the ETGs on the sky.We have removed the stripes from the data, focusing only on the largest contiguous SDSS footprint region.
ETGs are selected since they are bulge-dominated systems with velocity dispersions well described using the virial theorem.As a result, ETGs lie on the fundamental plane (FP) in the space of three-dimensional variables In the lower panels we present measured galaxy velocity v against redshift (lower left panel) and velocity measurement uncertainty ∆v against redshift (lower right panel).The black points in the lower panels represent the galaxies used in our analysis, and the yellow points those that are cut by our choice of velocity uncertainty limit.The uncertainty in the velocity measurement ∆v is negatively correlated with the measured velocity v, hence the ∆v cut removes the large negative velocity galaxies (cf.lower left panel).
log 10 R e = a log 10 σ 0 + bµ e + c (1) where σ 0 is the central velocity dispersion, and µ e is the mean surface brightness within half-light radius R e .
In principle, a and b can be deduced from the virial theorem, but practically a, b, c are empirically determined parameters.In Yoon & Park (2020) it was noted that a subsample of ETG galaxies exhibit a significantly smaller scatter on the FP plane and were selected to generate a catalog with smaller intrinsic scatter in the velocity reconstruction.Specifically, the FP of old ETGs with age≳ 9 Gyr has a smaller scatter of ∼ 0.06 dex (∼ 14% in the linear scale) than that of relatively young ETGs with age ≲ 6 Gyr that exhibits a larger scatter of ∼ 0.075 dex (∼ 17%).For the subsample of young ETGs, less compact ETGs have a smaller scatter on the FP (∼ 0.065 dex; ∼ 15%) than more compact ones (∼ 0.10 dex; ∼ 23%).By contrast, the scatter on the FP of old ETGs does not depend on the compactness of galaxy structure.The use of the FPs with smaller scatters allows for more precise distance measurements,2 which in turn facilitates a more accurate determination of the peculiar velocities of ETGs.
From the sample, we exclude galaxies with a large velocity uncertainty ∆v.In the lower right panel of Figure 5 we present the velocity uncertainty of each galaxy in the SDSS FP catalog as a function of redshift.We observe a systematic increase in velocity uncertainty as a function of z, which is typical in velocity measurements.The galaxies with extreme velocity uncertainties will generate a large, non-Gaussian noise contribution to the momentum power spectrum that is undesirable for our purposes.For this reason, we generate a Constrained B-spline fit in the ∆v -z plane (Ng & Maechler 2015;Kawinwanichakij et al. 2017;Braga et al. 2018;Li et al. 2019), and then cut all galaxies in the sample above the 95% percentile of the velocity uncertainty.The yellow points are those that are cut using this procedure, which constitute a total of 734 objects.This is approximately 5% of the 14, 180 FP galaxies that lie in the range 0.025 ≤ z ≤ 0.055 and within the main SDSS footprint.In the bottom left panel of Figure 5 we present the velocities of each galaxy as a function of redshift.The gold points are those cut by our procedure; we observe that those objects with a large negative velocity are removed.The velocity distribution is skewed towards negative values because galaxies are more likely to be scattered into the redshift range that out.Our choice largely removes this systematic artifact.Note that one could instead choose to weight each galaxy with some function their velocity uncertainty when constructing the momentum field (Qin et al. 2019).This achieves practically the same goal.
Making a velocity selection may introduce an artificial gradient in the number density, and this must be checked.In Figure 1 we present the number density of the SDSS FP catalog used in this work, as a function of comoving distance (top right panel).The black line is the number density without a cut, and the red line is the sample with the velocity uncertainty cut applied.The blue points/error bars are the sample mean and standard deviation of the number density of a set of mock samples, which are described in Section 4. Our selection does not introduce any significant redshift gradient, only acting to reduce the overall number density by ∼ 5%.However, the effect is considerably less than the 1σ density variation in the mocks that we use in our analysis (cf.blue points/error bars).
We also observe a large fluctuation in the number density of the full sample, with a higher number density of ETGs in the first three bins (cf.black line).We expect that this is a statistical fluctuation arising from the modest volume of the sample, because the number density of the full sample does not show a systematic gradient with distance.Rather, the first three bins show a ∼ 1σ high number density.We do not sub-sample the galaxies to remove this overdensity, on the grounds that statistical fluctuations should not be removed by pre-processing the data.This local overdensity is very likely due to the presence of the CfA great wall, spanning the distance of approximately 60 to 110 h −1 Mpc and consisting of rich galaxy clusters such as Coma, Leo, and many in the Hercules Superclusters.Regardless, the variation of the density is a ∼ 1σ effect, and is not unusual compared to individual mock samples.In Appendix C we present some n-z curves from a random sample of mocks, finding similar variation in number density with redshift.The mean of the mock data presents no systematic evolution in number density with redshift, as expected.

Power Spectrum Estimator
To estimate the power spectrum, we first bin the galaxy catalog into a regular lattice that we then Fourier transform.With this aim in mind, we generate a box of size L = 512 h −1 Mpc to enclose the data, and create a regular lattice of grid points with resolution ∆ = 512/256 = 2 h −1 Mpc per side.The observer is placed at the center of the box.All galaxies are assigned to the pixels according to the cloud-in-cell (CIC) algorithm, and for both the galaxy density and momentum fields each galaxy is inversely weighted by its KIAS VAGC angular selection function weight3 .We use the redshift of each galaxy corrected to the CMB frame, and assume a flat ΛCDM cosmology with parameters Ω m = 0.3, w de = −1 to infer the comoving distance.Once all galaxies have been assigned to the grid, we apply the SDSS angular selection function as a binary HEALPix4 mask of resolution N side = 512 (Gorski et al. 2005), zeroing any pixel if its angular weight is less than w < 0.9.We also apply a radial cut, and zero all pixels that lie r ≤ 80h −1 Mpc and r ≥ 160h −1 Mpc relative to the observer.The ℓ = 0 mode of the number den-sity, or momentum, power spectrum is then given by (Yamamoto et al. 2006), where F (r) is constructed from the galaxy density n(r), its mean n, and the line-of-sight velocity u(r) (Howlett 2019), and L ℓ (x) are the Legendre polynomials.We perform the integrals in Equation (2) as a double sum over unmasked pixels, which is a tractable procedure due to the limited volume of the data.Here, r and r ′ are vectors pointing to the pixels in the double sum.Using the Rayleigh expansion of a plane wave - where j ℓ is the spherical Bessel function of ℓ-th order, the integral over Ω k reduces to where y = r − r ′ .We are only extracting the L = 0 mode from the data, so take L L=0 (ŷ • r′ ) = 1.The contributions to the Ω k integral from higher ℓ-modes in the Rayleigh expansion are negligible, consistent with noise.
The FP peculiar velocity uncertainty of each galaxy is the dominant contaminant in the reconstruction of the momentum field, and in what follows we will treat the velocity uncertainty as a Gaussian white noise contribution to the momentum power spectrum.This is discussed further in Section 5.
We write the density of the matter distribution as ρ(x) and the corresponding cosmological velocity vector field as v(x).The density can be decomposed into a time dependent average ρ(t) and fluctuations δ according to ρ/ρ = 1 + δ(x).The momentum field is defined as the density-weighted velocity p(x) = [1 + δ(x)]v(x) (Park 2000).We can only measure the radial component of the galaxy velocities, p ∥ = (1 + δ)v ∥ where v ∥ = v • ê∥ , and ê∥ is the unit vector pointing along the line of sight.
We measure galaxy positions in redshift space.If we denote real and redshift space comoving distances as r and s respectively, they are related according to the relation where a and H are the scale factor and Hubble parameter, respectively.The density field in real (δ) and redshift ( δ) space are correspondingly related according to Kaiser (1987) δ where f = d log D/d log a is the linear growth rate, and D is the growth factor.The quantity in the square bracket on the right hand side is the radial redshift space distortion operator.This relation is valid to linear order in the perturbations δ and v.
The plane parallel approximation is often used in redshift space analysis.It is an approximation in which a constant, common line of sight vector is assigned to every galaxy.This neglects the radial nature of redshift space distortion, and is appropriate in the limit in which the galaxy sample is far from the observer and localised to a patch in the sky.With this approximation, the redshift space density field can be written in simplified form in Fourier space as where µ = k ∥ /k, and k ∥ is the component of the Fourier modes aligned with the line of sight vector, which is constant in the plane parallel limit 5 .This relation is valid in the large scale limit where the linear perturbation theory can be applied.The power spectrum of δ is Throughout this work we use an accent `to distinguish theoretical power spectra from the same quantities measured from galaxy distributions.Similarly, in the same limits (plane parallel, linear perturbations), the momentum and velocity fields are equivalent and the k ∥ component of the velocity power spectrum can be approximated as Linear perturbation theory is not sufficient for the extraction of information on quasi-linear scales.Also, the relations above are valid for the dark matter field and the galaxy bias -both linear and at higher order in the perturbations -must be accounted for.In a series of works (Vlah et al. 2012;Saito et al. 2014), the galaxy density and momentum power spectra to third order in redshift space have been constructed.The form is neatly summarized in Appendix A of Howlett ( 2019), and we direct the reader to that work and Vlah et al. (2012) for the complete derivation.The power spectra can be written as P δ (µ, k) = P 00 + µ 2 (2P 01 + P 02 + P 11 ) + µ 4 (P 03 + P 04 + P 12 + P 13 + P 22 /4) , We do not repeat the expressions for P ij in this work; they are combinations of convolutions of the linear matter power spectrum.Their exact functional form can be found in Appendix A of Howlett (2019) and Appendix D of Vlah et al. (2012).
In redshift space the power spectrum is no longer only a function of the magnitude of Fourier mode, due to the anisotropy generated along the line of sight 6 .The standard approach is to decompose the power spectra into multipoles; The non-linear power spectra (Eqs.12,13) are predicated on the global plane parallel assumption -there is 5 The density field in equation ( 9) cannot be physically realised without some window function convolution, as the condition implied by the plane parallel assumption -localisation to a patch on the sky -is not consistent with periodicity. 6If we additionally drop the plane parallel approximation, then translational invariance is lost and the power spectrum becomes a function of the angle between line of sight and separation vector between tracers.
a constant line of sight vector against which the angle µ between it and Fourier modes can be defined.In contrast, the power spectrum extracted from the data utilises the 'local plane parallel approximation', for which , where r and r′ are unit vectors.One can expect discrepancies between theoretical expectation and data due to these two different interpretations of the plane parallel approximation.However, higher order perturbation theory modelling almost exclusively utilises the global plane parallel limit7 .Since we expect breakdown of the plane parallel limit to only be relevant on the largest scales, where the statistical error is large, we do not anticipate any significant bias will occur by fitting the global plane parallel approximated power spectra to the data.However, this is a consequence of the modest data volume available; as the data quality improves (DESI Collaboration et al. 2016;Amendola et al. 2018;Doré et al. 2014) this issue must be addressed.One possibility is a hybrid splitting of The non-linear theoretical power spectra P δ (left panel) and P p∥ (right panel) for different values σv = (5, 6, 7)h −1 Mpc (pink, black and tan lines).The black dashed lines are the corresponding linear power spectra, with σv = 0.
small and large scale modes, as suggested in Castorina & White (2018).
The theoretical power spectra are sensitive to a large number of parameters -both cosmological and those pertaining to galaxy bias.In terms of cosmology, the shape of the galaxy density power spectrum is sensitive to Ω b h 2 , Ω c h 2 , n s and the amplitude to b 1 σ 8 , where b 1 is the linear galaxy bias and σ 8 is the rms fluctuation of the density.In addition, in redshift space both density and momentum power spectra are sensitive to the combination f σ 8 .Since we measure the statistics via biased tracers, at third order in perturbation theory we have four additional parameters -b 1 , b 2 , b 3,nl , b s .The definitions of the higher order bias terms can be found in Mc-Donald & Roy (2009); Howlett (2019).To linear order in the perturbations, the linear galaxy bias b 1 is completely degenerate with σ 8 , but this is not the case when including higher order contributions, and we must simultaneously vary b 1 and b 1 σ 8 .The quantities P ij used in equations (12,13) contain a velocity dispersion (labelled σ v , defined in units of h −1 Mpc).We treat σ v as an additional free parameter used to model the non-linear Finger-of-God effect.It has been argued elsewhere that multiple velocity dispersion parameters should be used to accurately model the Finger of God effect (Okumura et al. 2015), but for the limited range of scales that we measure, a single parameter is sufficient.Some concessions must be made -we cannot constrain all parameters from the available data, but fortunately, the power spectra are practically insensitive to a subset of them over the scales probed in this work.For this reason, we fix Ω b = 0.048, n s = 0.96, b s = −4(b 1 − 1)/7, and b 3,nl = 32(b 1 − 1)/315 (Saito et al. 2014).Practically, we have found that the non-linear power spectrum contributions pertaining to b 1 are the most significant.We fix Ω b as the baryon fraction has been measured accurately and independently using astrophysical sources and Big Bang Nucleosynthesis (BBN), which is sufficiently robust for our purposes.The primordial tilt n s will affect the shape of the matter power spectrum, but it is measured so accurately by the CMB that we can fix this parameter.We have found that order ∼ 10% variations of n s will not impact our final constraints, whereas the parameter is constrained to ∼ 0.4% accuracy from the CMB temperature data (Aghanim et al. 2020).Although it would be preferable to simultaneously constrain n s , Ω b , Ω m , and σ 8 , this is simply not feasible with the quality of data currently available.In the literature, multiple parameters pertaining to the shape of the galaxy power spectrum are fixed by convention.This highlights the fact that current velocity catalogues are not competitive cosmological probes compared to the CMB or expansion rate data such as type Ia supernova.However, improvements in both quality and volume of large scale structure data will allow us to further refine our analysis and perform a more complete cosmological parameter estimation in the future Saulder et al. (2023).
For similar reasons, the second order bias b 2 is fixed as b 2 = −0.05.We have tried varying this parameter and found our results to be insensitive to its value over the prior range −1 ≤ b 2 ≤ 1. Hence we fix b 2 = −0.05,which is the best fit value inferred from the mock galaxies (see Appendix I A).The power spectra are only very weakly sensitive to b s and b 3,nl over the scales probed, and our choice of these quantities does not affect our conclusions.We adopt the values quoted in (Saito et al. 2014) without any impact to our results.The final list of parameters varied in this work is Ω m , b 1 σ 8 , σ v , b 1 .The purpose of this paper is to infer a constraint on b 1 σ 8 and f σ 8 , so we treat f σ 8 = Ω 6/11 m (b 1 σ 8 )/b 1 as a derived parameter.
In Figure 2, we present the non-linear theoretical galaxy density (left panel) and momentum (right panel) power spectra, using fiducial parameters Ω m = 0.3, (Saito et al. 2014).We allow the velocity dispersion to take values σ v = 5, 6, 7h −1 Mpc (solid pink, black, tan lines) .We also present the linear power spectra as black dashed lines.The momentum power spectrum significantly departs from its linear limit even on intermediate scales k ∼ 0.05 h Mpc −1 , and the velocity dispersion σ v suppresses P p∥ and P δ , with the suppression in P p∥ entering on larger scales than P δ .
In the literature, the quantity f σ 8 is normally treated as a free parameter to be constrained by the data.In this work, we take a different approach by varying the standard cosmological and bias parameters to infer f σ 8 .Our reasoning is that the theoretical ensemble averages to which we compare the data to are entirely derived using General Relativity and the ΛCDM model, and for this model, f is fixed in terms of other cosmological parameters.Furthermore, beyond linear theory, the redshift-space momentum power spectrum depends on f due to higher order perturbative terms, thus only changing f σ 8 cannot fully accommodate the dependence of the theory on f .Varying Ω m allows us to predict the theoretical templates as accurately as possible.Using our approach, departures from the standard gravitational model would be detectable, via cosmological parameter posteriors that differ significantly from those obtained from other data sets such as the CMB.To map between parameters, we use the approximation f ≃ Ω γ m with γ = 6/11, which is accurate to sub-percent level for the parameter ranges considered in this work (Linder 2005).In non-standard gravity models, the exponent γ can ac-quire both redshift and scale dependence (Appleby & Weller 2010).
The premise of utilising the velocity field in cosmological parameter estimation is that the galaxy power spectrum amplitude is sensitive to both b 1 σ 8 and f σ 8 , whereas the amplitude of the velocity power spectrum only to f σ 8 , which breaks the degeneracy.However, the actual picture is complicated by two issues.First, the momentum power spectrum is a linear combination of various two-point functions, the dominant two being ⟨v ∥ (x)v ∥ (x ′ )⟩ and ⟨δ g (x)v ∥ (x)δ g (x ′ )v ∥ (x ′ )⟩ on large and small scales respectively.Most of the statistical constraining power is at small scales where the ⟨δ g (x)v ∥ (x)δ g (x ′ )v ∥ (x ′ )⟩ term dominates.As a result, adding the momentum power spectrum to the galaxy density power spectrum does not completely break the degeneracy between cosmological parameters, since f σ 8 and b 1 σ 8 increase the amplitude of both the momentum power spectrum on small scales and the galaxy density power spectrum on all scales.
Second, we are varying Ω m , and then inferring the growth rate f ≃ Ω 6/11 m .However, Ω m also changes the shape of the matter power spectrum, shifting the peak by changing the matter/radiation equality scale.Since we only measure the power spectra over a limited range of scales 0.02 h Mpc −1 ≤ k ≤ 0.2 h Mpc −1 , which does not include the peak, a shift in the peak position can be confused with an increase in amplitude.This introduces an additional source of correlation -two parameter sets can admit the same values of b 1 σ 8 and f σ 8 but present seemingly different power spectrum amplitudes over the range of scales probed.This introduces a three-way correlation between b 1 σ 8 , f σ 8 and b 1 .
In the Figure we plot the ratio of power spectra P δ (Ω m , σ 8 , b 1 )/ P δ fid and P p∥ (Ω m , σ 8 , b 1 )/ P p∥ fid , where the denominator is the power spectrum assuming the fiducial parameter set.We see that the parameters b 1 and σ 8 do not have a degenerate effect, as σ 8 changes both the overall amplitude and shape of P p∥ whereas b 1 only changes the shape.The sensitivity of the power spectra to Ω m arises from both f ∼ Ω γ m and the matter/radiation equality scale.We observe that variation of Ω m does not correspond only to an amplitude shift for either momentum or galaxy density power spectrum -the shape is sensitive to this parameter and we are ex- tracting information from both the amplitude and shape of these statistics.
We note that many of these issues will be ameliorated by reducing the statistical error associated with the large scale modes, by increasing the volume of data.Given the current data limitations, a three-way correlation between b 1 , f σ 8 (or Ω m ) and b 1 σ 8 is unavoidable.
Cosmological parameter dependence enters into our analysis in another way.When measuring the observed power spectra, we bin galaxies into a three-dimensional grid, using redshift and angular coordinates to generate comoving distances.This procedure is sensitive to Ω m , but only very weakly because the sample occupies a low redshift z ≤ 0.055.Hence when generating the field we fix Ω m = 0.3.The parameter h is absorbed into distance units throughout, but we must adopt a value of h when generating the theoretical matter power spectrum.We take h = h pl = 0.674 (Aghanim et al. 2020) for our analysis of the SDSS data, and h = h HR4 = 0.72 when applying our methodology to mock data, since this is the value used in the simulation.

Mask Convolution
The mask acts as a multiplicative window function on the field in configuration space, and hence a convolution in Fourier space.The power spectra extracted from the data is implicitly convolved, so to compare the measurement with the theoretical expectation value one must either de-convolve the data and mask, or convolve the theoretical power spectrum with the mask8 .We take the latter approach, which is more common in power spectrum analyses.One complication is the fact that the mask will couple the ℓ modes of the power spectrum due to the convolution.To proceed, we start with the ℓmodes of the unmasked power spectra P ℓ δ and P ℓ p∥ , and perform a Hankel transform to obtain the corresponding real space correlation function ℓ-modes - We then take the product of ζ ℓ δ, p∥ (r) with the mask ℓ-modes (Wilson et al. 2017) - where the c superscript denotes that the quantity has been convolved with the mask and Q ℓ (r) are the ℓmodes of the real space mask correlation function, computed from a random point distribution chosen to match the angular and radial selection of the survey volume.
We only extract the ℓ = 0 mode of the power spectrum from the data, since the data volume is not sufficient to obtain an accurate measurement of the higher multipoles.The ℓ-modes of the mask, Q ℓ (r), are obtained by generating a random point distribution of N = 10 6 points which encompass the radial and angular domain of the data, then binning pairs of points as a function of radial separation and angle µ and finally integrating out the µ dependence using the orthogonality of the Legendre polynomials.The function Q ℓ is normalised such that Q 0 (0) = 1.We present the Q 0,2,4 (r) modes of the mask correlation function in Figure 4.
To perform the masking procedure, in (15) we use the power spectra P ℓ δ and P ℓ p∥ measured from a mock galaxy catalog in a z = 0 snapshot box, using a simulation described in the following section.After performing the mask convolution and arriving at c P 0 δ, p∥ (k), we define the ratios r δ = c P 0 δ /P 0 δ and r p∥ = c P 0 p∥ /P 0 p∥ .To adjust the theoretical power spectra (12,13) to account for the effect of the mask, we multiply them by these precomputed ratios.We adopt this approach, rather than applying the Hankel transforms directly to the theoretical expectation values (12,13), because the perturbative expansion breaks down on small scales and this could generate artificial numerical artifacts when performing the k-space integrals.Practically, we find no significant change to the convolution when we include the Q ℓ=2 , Q ℓ=4 modes of the mask, compared to just using the Q ℓ=0 component in the convolution.We include them for completeness only.
Since we only extract the ℓ = 0 mode from the data, to simplify our notation in what follows we drop the ℓ = 0 subscript and simply denote the measured(theoretical) monopole power spectra as P δ, p∥ ( P δ, p∥ ), and the corresponding 'convolved' theoretical power spectra are c P δ = r δ P δ and c P p∥ = r p∥ P p∥ .

MOCK DATA AND COVARIANCE MATRICES
We initially extract the galaxy density and momentum power spectra from the Horizon Run 4 simulation z = 0 snapshot box (HR4; Kim et al. 2015), and then use this mock data to generate a covariance matrix for the SDSS data.We briefly describe the simulation used and present power spectrum and covariance matrix measurements.
HR4 is a cosmological scale cold dark matter simulation in which 6300 3 particles were gravitationally evolved in a V = (3150 h −1 Mpc) 3 box using a modified version of GOTPM code9 .Details of the simulation can be found in Kim et al. (2015).
Dark Matter (DM) halos are constructed using the friends-of-friends (FoF) algorithm with a linking length 0.2 times the mean DM particle separation for each snapshot.Only those with more than 30 DM particles are identified as halos, making the minimum halo mass 2.7 × 10 11 h −1 M ⊙ .Then, the merger tree of HR4 DM halos is constructed for 75 snapshots between z = 12 and 0 with a timestep of ∼ 0.1Gyr.During the construction of the merger tree, the most-bound halo particle (MBP) of each halo is traced, even when its hosting halo is merged into a more massive halo.
The HR4 mock galaxy catalog is constructed by following Hong et al. ( 2016) using an MBP-galaxy correspondence method.If an MBP is from an isolated halo, its position and velocity are regarded as those of the central galaxy of the given halo.On the other hand, if an MBP is from a halo already merged into a more massive halo, it is regarded as a candidate for a satellite galaxy.Then, the survival time of the satellite candidate after merger event t merge is calculated using a modified model of Jiang et al. (2008), where t dyn is the orbital period of a virialized object; M host and M sat are the masses of host and satellite halos, and epsilon is the circularity of the satellite's orbit.Only the satellite candidates whose merger events occurred no earlier than t merge are successfully identified as satellite galaxies.
Typically, dark matter only simulations underestimate the population of satellites in dark matter halos and the two-point correlation function imperfectly matches observations as the one-halo term is underestimated.In our scheme the timescale for satellites to merge with the central is taken into account by tracking the mostbound halo particles at the time of their dark matter halo entrance regardless of survival of the dark matter subhalo.By adjusting the timescale we can match the correlation function of the simulated galaxies with that of the SDSS Main Galaxy Sample (Zehavi et al. 2002).We direct the reader to Figure 1 of Hong et al. (2016) for a complete description.With our choice of algorithm, we do not have the freedom to vary the small-scale distribution of satellite galaxies; this is fixed by following the MBP dynamics and selecting the timescale to match the two-point clustering of the data.
We select mock galaxies by applying two independent mass cuts -one to central galaxies and the second to satellites, such that the total number of galaxies in the snapshot box is N = 1.66 × 10 8 and the satellite fraction is f sat = 0.4.These values were chosen to match the number density of the SDSS FP catalog n = 5.3 × 10 −3 h 3 Mpc −3 and roughly match the expected satellite fraction of ETGs (Mandelbaum et al. 2006).The central/satellite mass cuts for the HR4 mock galaxies are M cen = 9.3 × 10 11 M ⊙ /h and M sat = 6.3 × 10 11 M ⊙ /h, respectively.
We generate a regular grid of resolution ∆ = 3150/768 = 4.1 h −1 Mpc in each dimension, covering the snapshot box and assign each galaxy to a pixel according to a CIC scheme.For the scales being considered in this work -0.02 h Mpc −1 ≤ k ≤ 0.2 h Mpc −1 , we have checked that there are no resolution effects due to our choice of binning.
We generate real and redshift-space fields from the data.In real space, galaxy positions are given by their comoving position within the box.To generate a red-shift space field, the position of each galaxy is perturbed according to equation ( 7), where we take ê∥ = êz so that the line of sight vector is aligned with the z-Cartesian component of the box.Because we are using the redshift zero snapshot box and absorbing h into distance units, we take (aH) −1 = 10 −2 h −1 as the numerical factor between velocity and distance.We denote the real and redshift space galaxy density fields as δ and δ respectively, where δ ijk = (n ijk − n)/n, n is the average galaxy per pixel, and n ijk is the number of particles in the (i, j, k) pixel in a cubic lattice.
We also generate momentum fields by assigning galaxy velocities to the same grid as the number density.The momentum field studied in this section is given by p ijk is the sum of all galaxy velocities assigned to the (i, j, k) pixel.We can extract p (g) ijk from the snapshot using the real and redshift space positions of the galaxies, which we denote with/without a tilde respectively.The galaxy momentum field in redshift space -p(g) -is the quantity that we will measure from the data.
In Figure 5, we present the galaxy density power spectra (left panel) and momentum power spectra (right panel).In both panels, the dotted grey curve is the linear, dark matter density and momentum expectation value, and the blue dashed, red solid lines are the galaxy density/momentum power spectra in real, redshift space respectively.
The matter power spectrum presents familiar resultsthe mock galaxies have a linear bias b 1 ≃ 1.3 and the effect of redshift space distortion amplifies the power spectrum amplitude on large scales by an additional factor of 1+2β/3+β 2 /5, with β = f /b 1 .On small scales, starting at k ≳ 0.1 h Mpc −1 , the redshift space power spectrum exhibits a suppression in power due to galaxy stochastic velocities within their host halos (cf.red curve, left panel).
The momentum field presents very different behaviour.On the largest scales measurable by the simulation k ∼ 10 −2 h Mpc −1 , the real and redshift space momentum power spectra agree with the linear theory expectation value with no velocity bias (Zheng et al. 2015;Chen et al. 2018).However, on large scales 0.05 h Mpc −1 < k < 0.1 h Mpc −1 the galaxy redshift space power spectra departs from the real space measurement.It is clear that the shape of the momentum power spectrum is significantly modified by nonlinear effects.To linear order in the fields, the real and redshift space power spectra should be indistinguishable and so the deviation between the red and blue curves in the right panel are due to a combination of higher order per- turbative, and fully non-perturbative, non-linear effects.
In Appendix I B we consider the impact of non-linear velocity contributions on the momentum power spectrum further.
Next, we generate a covariance matrix for the SDSS data, using a set of N real = 512 non-overlapping mock SDSS catalogs from the snapshot box, placing mock observers in the box then applying the SDSS angular footprint and radial selection function relative to the observer position.When generating the mocks we take radial velocities relative to the observer and galaxy positions are corrected according to equation ( 7), where we take ê∥ = êr which is the radial unit vector pointing between observer and galaxy.The observers are placed in the rest frame of the snapshot box and the radial velocities of the galaxies are used to generate the momentum field.
For each mock catalog, we measure the galaxy density and momentum power spectra using the same methodol-ogy as for the SDSS data, then define a set of covariance matrices as where 1 ≤ i, j ≤ N k denote the Fourier mode bin k i , k j , (m, n) = ( δ, δ), ( δ, p∥ ), (p ∥ , δ), (p ∥ , p∥ ), y q,i = ln[P (m) q,i ] and ȳ(m) i is the sample mean of y (m) q,i in the i th Fourier bin.We use N k = 19 Fourier bins, linearly equi-spaced over the range k = [0.02,0.2] h Mpc −1 .Explicitly, P (m) q,i is the measured value of the m = δ, p∥ power spectrum from the q th realisation in the i th Fourier mode bin.We use y (m) q,i -the logarithm of the power spectra -as the observable, as the mocks indicate that the P (m) q,i measurements are not Gaussian distributed within the k i bins.The logarithm is sufficiently Gaussianized for our purposes.In Qin et al. (2019) a different transformation was used -the Box-Cox transform -to the same effect.
Finally, the total covariance matrix is the 2N k × 2N k symmetric, block matrix constructed from these three matrices; where Σ ( p∥ , δ) = Σ ( δ, p∥ ) .The inverse of the covariance matrix Λ −1 is used for parameter estimation.
In Figure 6, we present the normalised correlation matrices where σ is the standard deviation of the mocks in the i th Fourier bin (we are not using the Einstein summation convention).The left/middle/right panels are the (m, n) = ( δ, δ), (p ∥ , p∥ ), ( δ, p∥ ) matrices, respectively.We note the strong correlation between all Fourier bins of the momentum power spectrum (middle panel).The third panel shows an important property of the momentum field -it is positively correlated with the galaxy density field (Park 2000).By constructing statistics from the ratio of these fields, one may expect a reduction in cosmic variance and improved constraints on cosmological parameters.This possibility will be considered in the future.

RESULTS
We now present measurements of the galaxy density and momentum power spectra from the SDSS FP data, and then fit the standard model to the data to infer parameter constraints.

Power Spectra Measurements
In Figure 7, we present the galaxy density (left panel) and momentum (right panel) power spectra extracted from the SDSS velocity catalog (gold lines).We also plot the median and sample 68% range of the N real = 512 mock catalogs as blue points/error bars.The power spectrum of the FP data (cf.gold solid line, left panel) exhibits large fluctuations as a result of the modest volume of data -0.025≤ z ≤ 0.055.In particular, on intermediate scales 0.05 h Mpc −1 ≤ k ≤ 0.15 h Mpc −1 there is less power in the SDSS galaxies relative to the mocks.This is likely a statistical fluctuation.
In the right panel, the blue points are again the median and 68% range of the mocks and the dashed gold dashed line is momentum power spectrum extracted from the SDSS data.We observe that the SDSS power spectrum is dominated on small scales by the statistical uncertainty associated with the FP velocity measurements.This contribution is larger than the cosmological signal on scales k ≥ 0.1 hMpc −1 .To eliminate this contribution, we treat this velocity component as Gaussian noise, uncorrelated with the cosmological signal.Then, from the SDSS FP catalog we generate realisations by drawing a random velocity for each galaxy, using the velocity uncertainty as the width of a Gaussian distribution, then reconstructing the momentum power spectrum of these random velocities.We denote the resulting power spectrum as P fpn (k i ) (FP noise).We repeat this process N = 1000 times, and construct the average Table 1.The parameters varied in this work, and the prior range imposed.The quantity σv has units h −1 Mpc.
Pfpn (k i ) and standard deviation ∆P fpn (k i ).The pink dotted line is the resulting FP noise power spectrum Pfpn , and the solid gold line is the SDSS momentum power spectrum with this noise component subtracted.
To account for the uncertainty in the FP noise removal, the quantity ∆P fpn is added to the diagonal elements of the covariance matrix Σ Since we assume the FP uncertainty is uncorrelated with the actual velocity, has mean zero and is Gaussian distributed, this component will not contribute to Σ ( δ, δ) or Σ ( p∥ , δ) .The dimensionless quantity ∆P fpn To confirm our assumption of Gaussianity, in Figure 14 of Appendix C we bin the SDSS galaxies into nine comoving distance bins and generate histograms of the velocity uncertainties.The yellow bars are the entire sample, and the dark bars correspond to the galaxies used in our analysis, after a ∆v cut has been applied.We see that the cut acts to Gaussianize the velocity uncertainty, allowing us to treat this component as Gaussian white noise.
We take the measured galaxy density and momentum power spectra (cf.solid gold lines, Figure 7), and fit the mask-convolved theoretical model described in Section 3, arriving at constraints on a set of cosmological and galaxy bias parameters.We also apply our analysis pipeline to the mean of the mock realisations to confirm that we generate unbiased parameter constraints.

Parameter Estimation
We first use the mean of the N = 512 mock realisations to test our analysis pipeline.We define the data vector as , where y δ i = ln P δ i , y p∥ j = ln P p∥ j are the logarithm of the mean of the measured mock galaxy density and momentum power spectra in the 1 ≤ i, j ≤ N k Fourier bins, hence ⃗ d has dimension 2N k .
P δ i , P p∥ j are the theoretical expectation values derived in Section 3, and they have been multiplied by the pre-computed effect of the mask convo-lution; the ratios r δ and r p∥ respectively.We minimize the Gaussian likelihood L ∝ e −χ 2 /2 , with We vary the parameter set (Ω m , b 1 σ 8 , b 1 , σ v ), then subsequently infer the derived parameter f σ 8 ≃ Ω 6/11 m (b 1 σ 8 )/b 1 .The matter density Ω m enters the analysis in two ways -changing the amplitude of the power spectra via f σ 8 and modifying the shape of the linear matter power spectrum via the equality scale k eq .
The prior ranges for our parameters are provided in Table 1.Reasonable variation of these ranges does not significantly modify our results.We exclude σ v = 0 from consideration, since this quantity describes the Finger of God effect and must be non-zero.In addition, we apply a weak and physically motivated prior on the galaxy bias, such that b 1 = 1 +∞ −0.15 .With this choice, we penalise values of the bias b 1 < 1 with an additional Gaussian contribution to the χ 2 , of width σ b1 = 0.15.We do so because the SDSS galaxies are known to have bias b 1 ≥ 1, and the Gaussian width σ b1 = 0.15 is selected as this is the expected statistical constraint on the bias that can be obtained from the volume being probed in this work.In contrast, all values of b 1 larger than unity are selected with uniform prior, indicated with the +∞ upper bound.
We apply the same methodology to the SDSS extracted power spectra, taking the data vector ⃗ d = (y δ i − ln [ c P δ i ], y p∥ j − ln [ c P p∥ j ]) with y δ i = ln P δ i , y p∥ j = ln [P p∥ j − Pfpn,j ], and P δ i , P p∥ j are the measured power spectra from the SDSS catalog.We use the same parameter set and ranges as for the mocks shown in Table 1, covariance matrix Λ to generate χ 2 for the likelihood and the prior b 1 = 1 +∞ −0.15 .When generating the theoretical linear matter power spectrum, used in the theoretical power spectra P δ and P p∥ , we adopt a value of h = h HR4 = 0.72 for the mocks and h = h pl = 0.67 for the data.
The most relevant two-dimensional marginalised contours can be found in Figure 8, where we present the parameter pairs (f The colour scheme follows Figure 7 -blue empty contours represent the posteriors from mean of the mocks and the solid gold contours are our actual results; the cosmological parameters inferred from the SDSS FP data.The pale blue points are the parameter values used to generate the simulation -Ω m = 0.26, f = Ω γ m ≃ 0.48 and σ 8 = 1/1.26,b 1 = 1.3, which the blue contours successfully reproduce.The data favours a lower value of b 1 σ 8 , b 1 and σ v compared to the mocks, but a larger value of f σ 8 .The three-way degeneracy between b 1 , b 1 σ 8 and f σ 8 is apparent from the figures, and is consistent between the mocks and the data.The correlation matrices -the normalised covariance matrix between the parameters (b 1 σ 8 , f σ 8 , b 1 ) -are The one-dimensional posteriors for the parameters b 1 σ 8 and f σ 8 are presented in Figure 9.The correct input cosmology is recovered from the mean of the mocks (blue vertical dashed lines).The pink dashed line in the right panel is the Planck ΛCDM best fit f σ 8 = 0.43.We note that the peak of the probability distribution function (PDF) of f σ 8 is not exactly the expectation value, since the PDF is skewed towards larger parameter values.This is due to the condition that f ≃ Ω 6/11 m ≥ Ω 6/11 b ≃ 0.19, since we implicitly assume that Ω m cannot be zero due to the presence of baryons -we fix Ω b = 0.048.Low values of f σ 8 could be realised by allowing σ 8 to be arbitrarily small, and the amplitude of the galaxy power spectrum could in turn be compensated by allowing the linear bias b 1 to become arbitrarily large, as the matrix indicates, but this is not favoured by the data.The upshot is that f σ 8 , which is a derived parameter that is non-linearly related to Ω m , will not generically admit a symmetric posterior.
In Table 2, we present the expectation values and 68/95% limits on the parameters varied, for the mean of the mocks and the data.The data prefers a value of f σ 8 and Ω m that is in agreement with the Planck cosmology.The data also very marginally prefers a lower value of σ v , but this parameter is not well constrained by the data.The value of b 1 σ 8 is lower than the mocks, Figure 9. One-dimensional marginalised probability distributions (not normalised) for the parameters b1σ8 (left panel) and f σ8 (right panel).The color scheme is the same as in Figure 8.The pale blue dashed lines are the parameters used in the Horizon Run 4 simulation, and the pink dashed line is the Planck ΛCDM best fit f σ8 = 0.43, using Ωm = 0.3 and σ8 = 0.811 (Aghanim et al. 2020).
which is likely due to the drop in the galaxy density power spectrum on scales k ∼ 0.1 h Mpc −1 .
Finally in Figure 10, we present the measured galaxy density (left panel) and momentum (right panel) power spectra as solid gold lines, and the best fit theoretical power spectra as green dashed lines.The tan filled region is the 68% confidence region determined from the mocks.The theoretical curve for the galaxy density power spectrum does not fit the large scale modes well (cf.left panel, green curve); the excess of power on large scales relative to the best fit could be mitigated by a lower value of Ω m .We note that an excess of power at scales k ∼ 0.05 h Mpc −1 was also observed in independent datasets in Qin et al. (2019).However, we should be careful about performing a chi-by-eye fit to the data since the bins are correlated.The best fit is reasonable; χ 2 = 37.2 for the number of degrees of freedom N dof = 34 (38 data points, 4 free parameters).

DISCUSSION
In this work we extracted the galaxy density and momentum power spectrum from an ETG subset of the SDSS main galaxy sample, using FP information to infer velocities.We compared the measurements to power spectra derived using perturbation theory up to third order, which is formulated using the ΛCDM model.After testing our analysis on mock galaxy catalogs, we arrive at a constraint of b 1 σ 8 = 0.920 +0.070 −0.070 and f σ 8 = 0.471 +0.077 −0.080 at a mean redshift z = 0.04.Our analysis is consistent with other measurements of the same parameters in the literature (Fisher et al. 1994;Peacock et al. 2001;Blake et al. 2011;Beutler et al. 2012;de la Torre et al. 2013;Alam et al. 2017;Gil-Marin et al. 2020;Reid et al. 2012;Bautista et al. 2020;Qin et al. 2019).The momentum power spectrum is a sum of various two point functions, and hence sensitive to a combination of b 1 , σ 8 , f , Ω m , not simply to f σ 8 .For this reason, we tried to simultaneously fit for b 1 , b 1 σ 8 and Ω m , relying on the well-developed theo- retical framework of the ΛCDM and treated f σ 8 as a derived parameter.It is common practice in the literature to treat f σ 8 as a free parameter and then compare it to theoretical predictions.However, to extract this quantity from the galaxy density and momentum fields, we restricted our analysis entirely to the confines of the ΛCDM model.This is because we are extracting information from both the amplitude and shape of the power spectra.We expect that any modified gravity/cosmological model would not only change f σ 8 , but also the kernels used in perturbation theory, and introduce additional scales -such as scalar field masses -that will further modify the shape of the power spectrum.In this work we did not build the theoretical templates for such a general class of gravity models.Still, any deviation from the ΛCDM model could be detected by comparing our derived f σ 8 to values inferred from other probes.
Furthermore, we found that f σ 8 does not yield a symmetric marginalised posterior when fitting the ΛCDM model to the data, since f is a nonlinear function of Ω m that is in turn bounded by the condition Ω m ≥ Ω b .In this case, taking f σ 8 as a derived parameter and using the ΛCDM expectation f ≃ Ω 6/11 m is our preferred choice.By allowing f σ 8 to be an independent, free parameter with uniform prior, we might be making a subtly different assumption that could impact our interpretation of the result.The skewness of the f σ 8 posterior will be particularly pronounced at lower best fit values.Although combinations of measurements hint at a growth rate that departs from the ΛCDM model (Nguyen et al. 2023), it may be difficult to exactly ascertain the significance of any anomalous measurements of this parameter due to the non-Gaussian nature of the posterior.
We also comment on the non-linear velocity dispersion and its effect on the power spectrum.The Fingers of God produce a very pronounced effect on the momentum power spectrum, on surprisingly large scales.This is discussed further in the appendix, but we find an order ∼ O(20 − 30%) effect on scales k < 0.1h Mpc −1 .Nonlinear stochastic velocities generate a large decrease in power on large scales, followed by an increase in power on smaller scales.Such behaviour has already been noted in the literature (Koda et al. 2014;Dam et al. 2021).This phenomena will depend on the galaxy sample under consideration -central galaxies will be less affected by non-linear velocities.Understanding the Finger of God effect on the velocity and momentum statistics remains an interesting topic of future study.In this work, we used a single free parameter to model this effect.Given the quality and volume of data, this was sufficient.Future analysis will require more consideration.
The low redshift SDSS FP velocity data is consistent with other cosmological probes such as the CMB.However, the volume of data is modest and statistical uncertainties are large; the current data does not have strong model discriminatory power.In addition, we find that non-linear velocity dispersion strongly affects the momentum power spectrum, which forces us to introduce phenomenological prescriptions to model the effect.This philosophy is counter to the standard cosmology orthodoxy, which is to perform perturbation theory on an FLRW background.Certainly, there is room for alternative spacetime metrics to fit the data with equal precision, and this idea will be pursued in the future.The SDSS FP data is at or below the homogeneity scale assumed within the standard ΛCDM model, and searching for alternative prescriptions of the low redshift Universe is an interesting direction of future study.when generating redshift space fields.The galaxy velocity is always used as the observable when constructing the momentum field p ∥ and p∥ .
We use the full snapshot box in this section, and take e ∥ = e z .The quantities P δ g and P p∥ g are directly measurable from actual data, although by using the galaxy velocities one could also reconstruct the real space fields P δ and P p ∥ (Park 2000;Park & Choi 2005), albeit with large scatter due to the FP uncertainty.The halo redshift space quantities P δ h and P p∥ h are mock constructs that allow us to disentangle non-linear stochastic peculiar velocities between galaxies occupying the same host halo and the larger scale velocities that contain cosmological information.
In Figure 12 we present the galaxy density (top left) and momentum (top right) power spectra measured from these multiple fields.In the both panels, the grey, orange, cyan lines are the power spectra in real, halo-sourced and galaxy-sourced redshift space.The lower left/right panels contain the ratios P δ h /P δ , P δ g /P δ and P p∥ h /P p ∥ , P p∥ g /P p ∥ (cyan/orange lines) respectively.
The galaxy density power spectrum behaves as expected.The effect of redshift space distortion is to increase the amplitude of the power spectrum on all scales, and there is an additional suppression of power on small scales due to the stochastic velocities of bound structures (cf.cyan curve, bottom panel).The orange curve, which is a hypothetical redshift space using the host halo velocity to correct galaxy positions, almost entirely removes the small scale suppression due to the Fingers of God.
The effect of redshift space distortion on the momentum power spectrum is different.On large scales, all three power spectra approach a common value -although there is a ∼ 5% discrepancy between real and redshift space momentum power spectra even on scales k ∼ 0.01h Mpc −1 .The finger of God effect has a particularly large impact on P p∥ (cf.cyan line, right panels).It generates a suppression of power as large as ∼ 30% on scales k ∼ 0.05h Mpc −1 , with a subsequent increase in power on small scales.This behaviour has been observed in previous works (Koda et al. 2014), and is typically modelled using a series of phenomenological prescriptions such as exponential damping and a shot noise contribution.However, even when using the halo velocities to generate redshift space fields there is a suppression in power (cf.orange lines, right panels), which suggests that both nonlinear perturbative effects and the Finger of God Galaxy density (top left) and momentum (top right) power spectra extracted from the Horizon Run 4, z = 0 snapshot box.The grey/orange/cyan lines are the statistics in real space, halo-source redshift space and galaxy-sourced redshift space.In the lower panels we present the ratio of the two redshift space power spectra and the real space counterpart.contribution act to suppress the momentum power spectrum.This was also noted in Dam et al. (2021), who studied the velocity two-point statistics.

C. ANCILLARY INFORMATION
In this appendix, we present some ancilliary information regarding the mocks.In Figure 13 we present some n vs z curves from random mock samples, overlaid with the SDSS data.We see that the variation of the density within the SDSS data is consistent with typical variations observed within the mocks.
In Figure 14 we decompose the SDSS FP galaxies into nine redshift bins and generate histograms of the velocity uncertainty in each bin.The gold bars represent the entire sample, and the black bars the remaining galaxies after making a 5% ∆z cut as described in the main body of the text.The blue curves are Gaussian distributions fitted .Some random n vs z curves obtained from mock galaxy samples, with the SDSS data (solid black curve).We observe that the variation of n within the SDSS data is typical.
Bin dcm (h to each black histogram.We observe that the ∆v cut largely Gaussianizes the velocity uncertainty, which justifies our treatment of the FP noise as a Gaussian white noise component.In Table 4 we present the sample mean, standard deviation, skewness and kurtosis of the data in the nine redshift bins, before (after) the cut is made.The cut significantly reduces the non-Gaussian skewness and kurtosis of the samples.
Finally, we check that the mock data can approximately mimic the SDSS catalog at the level of the field and subsequent power spectrum extraction, by introducing a velocity uncertainty to the mock galaxies.For each of the 512 mock realisations, to each galaxy we add a velocity uncertainty ∆v.This is drawn from a Gaussian, of zero mean and width equal to the velocity uncertainty of the SDSS galaxy closest to the mock in comoving distance10 .In Figure 15 we present the (RA, dec) of a typical mock (upper left panel), the velocity against comoving distance (lower left panel) and the velocity error against comoving distance (lower right panel).In the upper right panel we present the median and 68% limits of the momentum power spectra extracted from these mocks, along with the SDSS measurement.We see that the mocks qualitatively reproduce the power spectrum of the data.Also, the v-z and ∆v-z plots (lower panels) are similar to those in Figure 1 from the SDSS data with the cut (black points)

Figure 1 .
Figure1.From the SDSS FP catalog, the angular distribution of ETGs used in this work (top left), the number density as a function of comoving distance with (red) and without (black) a velocity uncertainty cut (top right).The blue points/errorbars in the top right panel are the mean and standard deviation of the mock catalogs, discussed further in Section 4. In the lower panels we present measured galaxy velocity v against redshift (lower left panel) and velocity measurement uncertainty ∆v against redshift (lower right panel).The black points in the lower panels represent the galaxies used in our analysis, and the yellow points those that are cut by our choice of velocity uncertainty limit.The uncertainty in the velocity measurement ∆v is negatively correlated with the measured velocity v, hence the ∆v cut removes the large negative velocity galaxies (cf.lower left panel).

Figure 4 .
Figure 4.The ℓ = 0, 2, 4 modes of the real space mask correlation function, obtained by pair counting a random point distribution that matches the radial and angular selection functions of the data.The ℓ = 2 mode is negative, so we plot |Q2|.

Figure 5 .
Figure5.[Left panel] The galaxy density power spectrum extracted from the HR4 z = 0 snapshot box in real/redshift space (blue dashed/red solid).The grey dotted line is the linear dark matter (DM) power spectrum.[Right panel] The momentum power spectrum, with the same colour scheme as in the left panel.

Figure 7 .
Figure 7. [Left panel] Galaxy density power spectrum extracted from the SDSS FP data (solid gold line) and the median and 68% limits from N = 512 mock catalogs (blue points/error bars).[Right panel] The momentum power spectrum extracted from the SDSS FP data (gold dashed line), the contribution from the FP velocity uncertainty (pink dashed line) and the FP noise-subtracted momentum power spectrum (solid gold line).The blue points/error bars are the median and 68% range inferred from the mock catalogs.

Figure 8 .
Figure 8. Two-dimensional 68/95% marginalised contours on the parameters (f σ8, b1, b1σ8, σv).The empty blue contours are obtained from the mean of the N = 512 mock samples, and the solid gold are from the SDSS FP data.The pale blue points are the parameter input values used in the HR4 simulation.

Figure 10 .
Figure10.The measured galaxy density (left panel) and momentum (right panel) power spectra (gold lines), and the best fit theoretical power spectra to the SDSS data (green dashed lines).The filled region is the ±1σ uncertainty from the mock samples.
Figure12.Galaxy density (top left) and momentum (top right) power spectra extracted from the Horizon Run 4, z = 0 snapshot box.The grey/orange/cyan lines are the statistics in real space, halo-source redshift space and galaxy-sourced redshift space.In the lower panels we present the ratio of the two redshift space power spectra and the real space counterpart.
Figure13.Some random n vs z curves obtained from mock galaxy samples, with the SDSS data (solid black curve).We observe that the variation of n within the SDSS data is typical.

Figure 14 .
Figure 14.After binning the velocity uncertainties into nine redshift bins, we fit a Gaussian distribution to the cut data (blue curves).The dark bars are the cut SDSS data and yellow bars the full sample.The choice of cut approximately Gaussianizes the velocity uncertainty in each redshift bin.

Figure 15 .
Figure15.Angular distribution of a typical mock catalog used in this work (upper left), the velocity as a function of comoving distance after a velocity uncertainty has been added to each galaxy (lower left) and the velocity uncertainty as a function of comoving distance (lower right panel).In the top right panel we present the median and 68% range of the momentum power spectra extracted from the 512 mock catalogs, and the solid yellow line is the power spectrum extracted from the SDSS.

Table 2 .
The best fit and 68/95% ranges for the most relevant parameters studied in this work.

Table 3 .
The list of the power spectra considered.Different columns mean different spaces (real space, redshift space sourced by host halo velocities, and redshift space sourced by galaxy velocities).Different rows mean different quantities from which the power spectrum is calculated.

Table 4 .
The sample mean, standard deviation, skewness and kurtosis of ∆v in nine redshift bins, before and after the ∆v cut.The numbers in the brackets indicate the values after the cut has been applied.The sample is Gaussianized by our choice.