Momentum power spectrum of SDSS galaxies by massE cosmic ruler: 2.1$\times$ improvement in measure of growth rate

The distance measurement of the SDSS MGS elliptical sample with the massE cosmic ruler is detailed in Shi et al. (2022). Briefly, the massE ruler estimates the angular diameter distance of an elliptical galaxy with stellar mass well above 10⁸ $M_{\odot}$ as:

where $R_{\rm e,as}$ is the apparent effective radius in arcsec, $\sigma_{\rm e}$ is the velocity dispersion within $R_{\rm e,as}$, $M_{\rm star,1Mpc}$ is the galaxy stellar mass placed at a distance of 1 Mpc and $z_{\rm cos}$ is the cosmic redshift of the galaxy. Besides the above four observables, the two nuisance parameters in Equation 3 are $\beta$ and $S_{D_{0}}$, where $\beta$ is the slope of the massE relationship, and $S_{D_{0}}$ is the intercept of the relation and is related to the absolute calibration of observables and is degenerate with the local Hubble constant $H_{0}$. By fitting to the low-redshift sample compiled in Shi et al. (2021), Shi et al. (2022) estimated:

and

The angular diameter distance shown in the above equation can be converted the luminosity distance through $D_{\rm L}$=$D_{\rm A}$(1+z)². For a flat cosmology, the angular distance is related to the comoving distance through $D_{\rm c}$=$D_{\rm A}$(1+z).

In practice, the observed redshift $z_{\rm obs}$ is used to approximate $z_{\rm cos}$. Although the effect on the distance is negligible, the effect on the PV is noticeable as shown in Equation 10. By assuming $z_{\rm cos}$=$z_{\rm obs}$, the shift in the comoving distance can be derived as:

in which we define $\kappa$=$(\frac{2\beta+0.25}{2\beta-0.25})\frac{1}{c}$, and assume that $v_{\rm p}$ is a small number as compared to the speed of light.

As shown later, the PV is not sensitive to the accurate calibration of the local massE relationship, i.e., $\beta$ and $S_{D_{0}}$ can vary significantly while the derived PV and momentum power spectrum remains almost unchanged.

The procedure to construct the PV catalog is shown in Table 1 that lists the number of objects for each step. We start with the SDSS Large-Scale-Structure (LSS) safe-0 catalog (Blanton et al., 2005) whose completeness and selection function have been quantified ¹¹1http://sdss.physics.nyu.edu/vagc/lss.html, http://sdss.physics.nyu.edu/lss/dr72/. The redshift in the catalog is converted from the heliocentric frame to the CMB frame following Fixsen et al. (1996).

To define the elliptical sub-sample, we first cross-match with the MPA-JHU catalog with a search radius of 1.5 arcsec to query the stellar mass, Petrosion half-light radius in the $r$ band (PETROR50_R), velocity dispersion (V_DISP) and the fraction of de Vaucouleurs component in the $r$ band (FRACDEV_R) from the MPA-JHU catalog (Kauffmann et al., 2003; Brinchmann et al., 2004). The fraction of the successful match is 0.938, which will be included in the galaxy weight when calculating the momentum power spectrum. The SDSS MGS-elliptical sample is then defined as FRACDEV_R $>$ 0.80, PETROR50_R > 0, along with V_DISP > 0 and 0.02 $<$ $z_{\rm obs}$ $<$ 0.2. The lower redshift limit is adopted so that even in the lowest redshift bin of [0.02, 0.021], 95% of galaxies still have stellar masses larger than 3$\times$10⁸ M_⊙ to ensure the validity of the distance calculation with Equation 3. The upper redshift limit is set to ensure a reasonable high galaxy number density. The footprint of the MGS-elliptical sample is shown in Figure 1(a).

To measure PV for the MGS-elliptical sample, we first derive quantities as needed by the distance calculator given by Equation 3. As detailed in Shi et al. (2022), the apparent size $R_{\rm e,as}$ is derived from the Petrosian half-light radius PETROR50_R by adopting the cosmological model of the MPA-JHU catalog, and the flux-like $M_{\rm star,1Mpc}$ is converted from the stellar mass in the MPA-JHU catalog using the same cosmological model. The velocity dispersion within the fiber aperture from that catalog is corrected to the effective radius to give $\sigma_{\rm e}$.

In practice it is not recommended to estimate $v_{\rm p}$ directly from Equation 2, because the log-normal distribution of the distance errors would make errors of $v_{\rm p}$ non-Gaussian, thus complicating the error analysis of momentum power spectra and cosmological parameter constraints. In previous studies (Adams & Blake, 2017), it has been demonstrated that $v_{\rm p}$ can be estimated from the logarithmic ratio of the comoving distance at $z$=$z_{\rm obs}$ and that at $z$=$z_{\rm cos}$, i.e., $\eta_{\rm gal}$ = ln($\frac{D_{\rm c}(z_{\rm obs})}{D_{\rm c}(z_{\rm cos})}$):

where

with requirements of $v_{\rm p}$/c $\ll$ 1 and $\alpha{v_{\rm p}}$ $\ll$ 1. The typical $v_{\rm p}$ is about 300 km s^-1  so that the first condition is satisfied. The second condition requires $z_{\rm obs}$ to be reasonably larger than 0, e.g. $\alpha{v_{\rm p}}$ $\sim$ 0.05 at $z_{\rm obs}$=0.02 and $v_{\rm p}$=300 km s^-1.

In the above equation, the true comoving distance of a galaxy, i.e., $D_{\rm c}(z_{\rm cos})$, is measured directly from the massE cosmic ruler. In practice we replace $z_{\rm cos}$ with $z_{\rm obs}$ in the massE equation, which causes a systematic shift in ln($D_{\rm c}(z_{\rm cos})$) as derived in Equation 6. To account for this shift, an updated equation to derive PV for the massE cosmic ruler is given as

where $\kappa$ (defined in Equation 6) dominates over $\alpha$ at increasing redshift and reduces PV errors accordingly compared to Equation 8.

The remaining key is to obtain unbiased estimate of the comoving distance at the observed redshift of the galaxy, i.e., $D_{\rm c}(z_{\rm obs})$. Our overall strategy is to divide the interested redshift range into small bins and estimate $D_{\rm c}(z_{\rm obs})$ therein. The redshift bin size needs to balance two facts. On the one hand, galaxies in each redshift bin are composed of those located at slightly lower true distances with $v_{\rm p}$ > 0 and those located at slightly larger distances with $v_{\rm p}$ < 0. As a result, an increasing redshift bin size contains a larger number of galaxies to converge their mean true distance ${\langle}{\rm ln}(D_{\rm c}(z_{\rm cos})){\rangle}$ to the distance at $z_{\rm obs}$, i.e., ln($D_{\rm c}(z_{\rm obs})$). On the other hand, the redshift bin cannot be too large so that $D_{\rm c}(z_{\rm obs})$ is no longer a constant over the bin.

Our fiducial measurements adopt a redshift bin size of $d{z_{\rm obs}}$=5$\times$10^-3. In each redshift bin, as shown in Figure 2(a) for an example, the ln($D_{\rm c}(z_{\rm cos})$) distribution is more or less normal. But because a small fraction of the sample has failed measurements in the velocity dispersion or galaxy sizes, there are outliers. We remove outliers in the following ways: for those within three times the standard deviation of the best-fit Gaussian, we include all of them; for those outside, we randomly include part of them so that the number of objects follows the Gaussian expectation at its ln($D_{\rm c,z_{\rm obs}}$) bin. The new distribution is shown in Figure 2 (b). The exclusion of these outliers decreases the sample completeness and increases the error of the momentum power spectrum. The $\eta_{\rm bin}$ is then defined as the minus of the offset of a galaxy from the Gaussian mean in that redshift bin. Because the redshift bin size is not infinitely small, we correct additional small offset between the observed redshift and the mean redshift of galaxies within the bin:

where ${\langle}z{\rangle}$ is the mean redshift of galaxies within the given redshift bin. Besides the fiducial redshift bin size, we show that if varying $dz_{\rm obs}$ from 1$\times$10^-3 to 1$\times$10^-2, the result remains unchanged as discussed later.

Our final PV catalog, listed in Table 1 and illustrated in Figure 3, contains 229,890  galaxies, representing the largest number of objects in any PV catalog thus far. It is seven times larger than the previously largest one (Howlett et al., 2022) and even exceeds the expected number of PV objects generated by applying the FP plus Tully-Fisher methods to the DESI survey (Saulder et al., 2023). At $z$ $<$ 0.1 we include the whole SDSS MGS early-type galaxies with few selection, while the SDSS-FP sample has additional selections to have higher PV precision for individual galaxies, such as a brighter flux limit, galaxy inclination angle, H$\alpha$ equivalent width etc. As a result, our sample at $z<$0.1 contains a factor of 2.5  more galaxies than the SDSS-FP one. Another important difference is because of our ability to apply the massE method up to a redshift of 0.2, which is well beyond the limit ($\sim$0.1) of the Tully-Fisher and FP methods. The redshift limit of the latter two methods is not a rigid boundary based on the first principle, but rather derived from empirical usage. The Tully-Fisher method requires spatially-resolved kinematic measurements, limiting its redshift range. Although the FP can rely on single fiber measurements, its zero point is a function of galaxy properties (Magoulas et al., 2012; Howlett et al., 2022). While adjustments can be made to the FP zero point based on galaxy properties, even a small remaining systematic shift may lead to a large offset in the PV estimation at relatively high redshifts. On the hand, the massE contains solely two free parameters that correspond to the intercept and slope of the massE relationship, respectively. Consequently, we divide the redshift range into narrow individual bins, whereby PV is measured by the offset from the Gaussian mean of galaxies in narrow redshift bins, largely removing the necessity of accurately measuring the zero point throughout the entire redshift range. We will present the investigation using mock data-sets to support the above statement in § 4.4. We do not see any apparent limit for massE to extend beyond a redshift of 0.2. However, given the rapid drop in the galaxy number density of the SDSS MGS toward higher redshift, in this study, we define our sample to be below a redshift of 0.2.

For the momentum power spectrum measurement, we select those that have completeness fgotten $\geq$ 0.9 and lie within the contiguous region of the North Galactic cap. The corresponding footprint area is 6725 deg² as illustrated in Fig. 1(b).

In the above approach to calculate PV, we use the Gaussian mean of the distance distribution within a small redshift bin as the velocity zero point. This assumption is correct as along as the large scale structure within the redshift bin is not coherent so that no bulk motion exists. The best way to assure this is an all-sky survey that covers both northern and southern hemisphere. However the SDSS is not large enough especially for galaxies below $z$=0.1. Fortunately, the SDSS MGS have additional coverage along the equatorial plane besides a continuous area in the northern hemisphere. As a result, the Gaussian mean is measured only using galaxies in two great circles of the sky sphere, i.e., those along the equatorial plane plus a stripe along longitude as enclosed by the dashed boxes in Figure 1 (a). We demonstrate this with mock data in § 4.3.

The SDSS MGS is a flux-limited sample. Due to the luminosity function, sources that are intrinsically faint and below the flux limit have a higher likelihood of scattering to be above the limit, compared to sources that are intrinsically bright, which have a lower chance of falling below the limit. As a result, some fraction of sources that are above the flux limit should be intrinsically faint, thus overestimating its stellar mass and underestimating its distance.

The log distribution of $r$-band magnitude errors as a function of the magnitude is shown for the entire photometric catalog from $m_{r}$=14 to 24 in Figure 5(a). We fit a double Gaussian function to the log distribution of magnitude errors in magnitude bins with $\Delta$$m_{r}$=0.4 mag. An example of the fitting is illustrated in Figure 5(b) for $m_{r}$=17.8-18.2. The SDSS MGS sample for the large scale structure has a flux limit around $m_{\rm r}$=17.6. Below the redshift of 0.1, this flux limit still probes those below the characteristic magnitude as shown in Figure 5(c).

We randomly draw 3$\times$10⁶ absolute magnitudes from the luminosity function for a given redshift. The distribution ranges from $M_{\rm r}$=$-$24 to $M_{\rm r}$=23.25-DM, where DM is the distance module, with 23.25 corresponding to the SDSS image 3-$\sigma$ depth. Next, we assign the flux error to each absolute magnitude following Figure 5(a), and calculate the difference in the magnitude (${\Delta}M_{r}$) between the intrinsic and observed ones. Since the mass-to-light ratio is a function of color but SDSS MGS is a purely flux-limited sample, this magnitude difference corresponds to a difference in distance such that ${\Delta}{\rm ln}{D_{\rm c}}$ = $\frac{{\rm ln}(10)}{2.5}$$\frac{\beta}{(2\beta-0.25)}$${\Delta}M_{r}$.

We then obtain the ${\Delta}{\rm ln}{D_{\rm c}}$ distribution for two separate samples: one with intrinsic apparent magnitudes below the SDSS MGS bright limit and the other with observed apparent magnitudes below the limit, as shown in Figure 5(d) for $z$=0.2. The difference between these two distributions arises from the fact that the latter is affected by the flux limit cut.

In Figure 5(e), we show the fractional difference (the denominator is the observed distribution) in individual bins of ${\Delta}{\rm ln}{D_{\rm c}}$/ std(${\rm ln}{D_{\rm c}}$). Here, std(${\rm ln}{D_{\rm c}}$) represents the standard deviation of ${\Delta}{\rm ln}{D_{\rm c}}$, and the offset is expressed in terms of standard deviations. This approach ensures that the above trend remains unchanged if ${\rm ln}{D_{\rm c}}$ experiences increasing Gaussian errors as in the real case. After conducting 10 simulations for a given redshift, we fit a linear function to the above trend with zero intercept. We obtain the slope of the trend for $z$=0.05, 0.1, 0.15 and 0.2, and then fit a power law to interpolate at any given redshift. To correct the bias introduced by the flux limit for each redshift bin in real case, we select a portion of galaxies according to its offset from the Gaussian mean and the linear function at given redshift, and reverse its sign of the offset.

The line-of-sight mass-weighted momentum field is described as

where $\delta(\textbf{{r}})$ is the galaxy over-density at the location vector r. For a catalog with discrete galaxies, it can be written as (Feldman et al., 1994; Hand et al., 2017; Howlett, 2019):

where $w(\textbf{{r}})$ is the galaxy total weight that is $\geq$ 1, $n_{\rm g}(\textbf{{r}})$ is the observed galaxy number density, and $A$ is the normalization given in Equation 18 below.

The momentum power spectrum $P^{P}(\textbf{{k}})$=${\langle}$$F^{P}(\textbf{{k}})F^{P}(\textbf{{k}}^{\prime})$${\rangle}$, the Fourier analog of the two-point correlation function ${\langle}F^{P}(\textbf{{r}})F^{P}(\textbf{{r}}^{\prime}){\rangle}$ where k=$\frac{2{\pi}(\hat{\textbf{{r}}}-\hat{\textbf{{r}}}^{\prime})}{|\textbf{{r}}-\textbf{{r}}^{\prime}|}$ ²²2The hat symbol denotes the unit vector., can expand through Legendre polynomials:

where $\mu=\hat{\textbf{{k}}}\cdot\hat{\textbf{{r}}}_{h}$ and $\hat{\textbf{{r}}}_{h}=(\hat{\textbf{{r}}}+\hat{\textbf{{r}}}^{\prime})/2$. The estimator for the $P^{P}_{\ell}(k)$ is given by (Feldman et al., 1994; Hand et al., 2017):

where $\Omega_{k}$ is the solid angle in Fourier space. By assuming the local plane-parallel approximation of ($\hat{\textbf{{k}}}\cdot\hat{\textbf{{r}}}_{h}$)$\approx$($\hat{\textbf{{k}}}\cdot\hat{\textbf{{r}}}$)$\approx$($\hat{\textbf{{k}}}\cdot\hat{\textbf{{r}}}^{\prime}$), the integral over $F^{P}(\textbf{{r}})$ and $F^{P}(\textbf{{r}}^{\prime})$ in the above equation can be separated into independent integrals (Yamamoto et al., 2006). By further decomposing the $\hat{\textbf{{k}}}{\cdot}\hat{\textbf{{r}}}$ into a product of spherical harmonics, the calculation of the above equation can be done with a small number of Fast Fourier Transforms (FFTs), which significantly saves the computing time (Hand et al., 2017). In this study, we are only interested in $\ell$=0, i.e. the angle-averaged momentum power spectrum $P^{P}_{\ell=0}(k)$, so that the local plane-parallel is exactly valid.

The python code nbodykit (Hand et al., 2018) implements the above calculation for the density power spectrum. We make some minor modification of their python codes FKPCatalog and ConvolvedFFTPower to produce FKPMomentumCatalog and MomentumConvolvedFFTPower for the momentum power spectrum calculation. Here we list the key inputs in order to run the code:

(1) $w_{\rm mpa-jhu}$ – the completeness weight due to the failure match of the LSS catalog with the MPA-JHU catalog: the reciprocal of 0.938  as derived in Table 1 is taken as $w_{\rm mpa-jhu}$ for all galaxies.

(2) $w_{v_{\rm p}}$ – the completeness weight for the $v_{\rm p}$ measurement that is due to the failure in measurements of $\sigma_{\rm e}$: we divide the RA, DEC and $z$-range of our sample into 10 bins, respectively. In each 3-D cell, the completeness is calculated as the ratio of the number of objects with available $v_{\rm p}$ and the total number. The completeness is found to be almost a constant with the median value of 0.977  and a standard deviation of 0.029. As a result we assign $w_{v_{\rm p}}$ as a reciprocal of the median completeness for all galaxies.

(3) $w_{\rm sel}$ – the weight for the selection function of the SDSS main galaxy LSS sample. This number is available in the catalog, which is 1/fgotten.

(4) The observed number density $\bar{n}(z)$ as a function of redshift: this function is shown in Figure 4(a). We interpolate the function to get the number density at the redshift of each galaxy.

(5) The $\langle{v_{\rm p}^{2}}\rangle_{z}$ as a function of redshift: this function is shown in Fig. 4(b). Similarly we interpolate it to obtain the velocity variance at the redshift of each galaxy.

(6) The FKP weight

where $P^{\rm P}$ is set to be $10^{9}$ $h^{-3}$ Mpc³ km² s^-2, which is roughly the momentum power spectrum at our median spatial scale of 0.1 $h$/Mpc.

(7) The range of spatial scales as bracketed by $k_{\rm min}$ and $k_{\rm max}$: $k_{\rm min}$ is crucial for the measurement of structure growth that is sensitive to large spatial scales. We set it to be 0.02 $h$/Mpc. Our test with mock data shows that the momentum power spectrum is recovered well above $k_{\rm min}$. For $k_{\rm max}$, we test several values and found that the growth rate is insensitive to it.

(8) With the combined completeness weight $w_{\rm complete}$=$w_{\rm mpa-jhu}$$w_{v_{\rm p}}$$w_{\rm sel}$ and FKP weight $w_{\rm FKP}$, we ran FKPMomentumCatalog and MomentumConvolvedFFTPower implemented in nbodykit to get the momentum power spectrum which is further subtracted by the shot-noise contribution:

where $A$ is the normalization as given by

Given $w_{\rm tot}$=$w_{\rm complete}$$w_{\rm FKP}$, the effective redshift of the sample is given by

where $i$ and $j$ are members of all galaxy pairs in the sample. For our whole sample, $z_{\rm eff}$ is found to be 0.0959.

We use 32 Uchuu-SDSS mock galaxy catalogs that are constructed from the Uchuu 2.1 trillion N-body simulation (Dong-Páez et al., 2022). The catalogs match the depth, footprint, selection function etc of the SDSS MGS in the North Cap region. To construct the mock catalog, the $r$-band absolute magnitude $M_{r}$ was assigned to a halo through subhalo abundance matching (SHAM) techniques. The $g$-$r$ color was further assigned randomly through the probability distribution of the color-magnitude diagram. With $M_{r}$, $g$-$r$ color and observed redshift, k-d tree was used to find the closest real galaxy for each mock galaxy, after which properties of real galaxies were attached to mock galaxies. For our use, we need one additional parameter – FRACDEV, which is obtained with k-d tree match also based on $M_{r}$, $g$-$r$ color and observed redshift.

We define elliptical galaxies of each mock catalog following the same criteria as the observed one: (1) FRACDEV $\geq$ 0.8 and (2) 0.02 $<$ $z_{\rm obs}$ $<$ 0.2. With true peculiar velocities obtained through Equation 2, the true momentum power spectrum is then calculated using nbodykit. The median of intrinsic momentum power spectra of mocks is presented as solid lines in Figure 8.

To produce the PV for the mock galaxy catalog, we carry out the following simulation as outlined in Figure 6:

(1) The apparent half-light radius in arcsec ($R_{\rm e,as}^{\rm mock}$) for each mock galaxy: we first plot the observed physical half-light radius as a function of $M_{r}$ of the SDSS MGS-elliptical sample in Figure 7(a). Here we include ellipticals with $z_{\rm obs}$=0.01-0.2 to extend the dynamic range. As shown in the figure, we fit the median and standard deviations of half-light radii at given $M_{r}$ with second-order polynomial functions, respectively. We then assign the intrinsic physical size to each mock galaxy based on the above median trend along with its mock $M_{r}$. The associated errors are assigned based on standard deviations of trends (dashed lines in Figure 6 (a)). The mock physical size is then further converted to the angular size $R_{\rm e,as}^{\rm mock}$ based on the mock true angular diameter distance.

(2) The stellar mass normalized at 1 Mpc ($M_{\rm star,1Mpc}^{\rm mock}$) for each mock galaxy: we first plot the observed mass-to-light ratio in terms of log$M_{\rm star}$+$M_{\rm r}$/2.5 as a function of $g$-$r$ color based on the SDSS MGS in Figure 7(b). Here we include all types of galaxies in order to extend the dynamic range. Again, second-order polynomial functions are fitted to the median and standard deviations of the relationship, respectively. We then assign an intrinsic stellar mass to each mock galaxy based on its mock $g$-$r$ color and $M_{\rm r}$ using the median trend, while the observed errors are assigned based on the fitting to the standard deviations of the relationship (dashed lines in Figure 6 (b)). The stellar mass is then converted to that at 1 Mpc ($M_{\rm star,1Mpc}^{\rm mock}$) based on the mock true luminosity distance.

(3) The velocity dispersion ($\sigma_{\rm e}^{\rm mock}$) for each mock galaxy: we first use Equation 3 to calculate an intrinsic velocity dispersion from the intrinsic $M_{\rm star,1Mpc}^{\rm mock}$ and $R_{\rm e,as}^{\rm mock}$ along with the mock true distance, for which the slope $\beta$ of the massE relationship is fixed to the best-fit value in Equation 4. About 10⁴ objects of the PV sample have duplicate measurements of the velocity dispersion in the MPA-JHU catalog. As shown in Figure 6 (c), the distribution of the velocity dispersion errors can be described by two Gaussian profiles. We thus assign each galaxies following the error distribution in the figure.

(4) The mock distance for each mock galaxy: based on the mock $R_{\rm e,as}^{\rm mock}$, $M_{\rm star,1Mpc}^{\rm mock}$ and $\sigma_{\rm e}^{\rm mock}$ with errors, we calculate the mock distance through Equation 3. Here we perturb the $\beta$ and the corresponding $D_{0}$ with their errors, respectively.

(5) The mock PV for each mock galaxy: we then follow the same procedure as the observed one to obtain the PV for each mock galaxy. For each of 31 mocks from mock-0 to mock-30, we run 60 simulations so that in total we have 1860 simulations for each observed set of PV measurement. Note that the last mock (mock-31) is used to test the zero velocity point. Each simulation includes random perturbation of the intrinsic $R_{\rm e,as}^{\rm mock}$, $M_{\rm star,1Mpc}^{\rm mock}$ and $\sigma_{\rm e}^{\rm mock}$ with their error distributions, and one perturbation of $\beta$ and $D_{0}$, respectively.

By construction, the mock catalog does not have failure measurements, so we do not remove any outliers. We confirm that the velocity errors produced by mocks are similar to the observed ones at different redshift.

In this section, we examine the impact of the velocity zero point on the momentum power spectrum. In our method to measure the PV, the zero velocity represents the Gaussian mean of the distance distribution within a narrow redshift bin. We first present the case where each mock uses its own distance distribution to calculate the Gaussian mean. As a result, the zero point is derived from galaxies in an about 7000 deg² area in the northern hemisphere. As shown in the first row of Figure 8, in the low redshift range of $\Delta$z=0.02-0.10 the derived momentum power spectrum is systematically lower than the intrinsic one, while the one in $\Delta$z=0.10-0.20 overall recovers the intrinsic one very well except for small k that is below 0.02 $h$/Mpc. This demonstrates that because of the limited volume at low redshift a small bulk motion still exists and affects our methodology.

The Uchuu-SDSS provides 32 independent mocks to simulate the SDSS MGS in the northern hemisphere. In the second case, we use the last mock (mock-31) to emulate a SDSS-MGS-like survey in the southern hemisphere. For each of the remaining mocks, we combine their galaxies with those from mock 31 to measure the Gaussian mean. This mean is then used to derive PV for the galaxies in that particular mock. Note that galaxies from mock 31 are solely used for determining the Gaussian mean and are not used in the calculation of the momentum power spectra. As shown in the middle panel of Figure 8, the intrinsic momentum power spectra for three different redshift ranges are well recovered. This indicates, even at redshifts below $z$=0.1, the effect of the velocity zero point on the momentum power spectrum is negligible if large enough survey is conducted.

Unfortunately, the SDSS MGS lacks a comparable survey in the southern hemisphere. But it does contain galaxies along the equatorial plane. In the third case, we thus combine part of galaxies from mock 31 with part of galaxies from each of the remaining mocks to calculate the Gaussian mean for that mock. For the part from mock 31, we select galaxies with a declination below 5 degrees, which roughly yields a similar sky area to the part of the SDSS MGS that lies in the equatorial plane but not in the continuous northern regions. This sky area is represented by the two orange boxes along the equatorial plane in Figure 1. For each of the remaining mocks, we select their galaxies along the longitude plane that is similar to the stripe in Figure 1. The above two parts of galaxies lie on two perpendicular great circles of the sky sphere, and thus effectively reduce the coherent bulk motion. As shown in the lower panel of Figure 8, in both the high and low redshift bins, the intrinsic momentum power spectra are recovered reasonably well. We thus use this strategy for both the observation data and mock data.

To assess the validity of our methodology, we conduct two types of comparison between the massE-based measurements and the true ones in mock data-sets. The first one is the PV of individual galaxies, while the second focuses on the momentum power spectrum, including both the monopole and quadrupole. Figure 9 presents a comparison between the derived PV and the true velocity in three redshift bins for galaxies, where the derived PV is the mean value after 2000 simulations. The two velocities in all three redshift bins exhibit an almost one-to-one correlation. This suggests that, on a statistical level, the massE effectively recovers the true velocity. A small offset is observed in the zero velocity, which represents a random error on the zero velocity point. This offset adds to the shot noise (Howlett, 2019), resulting in a slight increase in the momentum power spectrum above 0.1 $h$/Mpc, as seen in Figure 10. However, this should not affect the constraint on the growth rate, as it is not sensitive to the high $k$ regime.

To further evaluate the velocity field, Figure 10 presents a comparison of the momentum power spectra between the massE-based measurements and true ones. The solid line indicates the mean value from 31 mocks, while the symbols represents the mean values of massE-based measurements for 31 mocks with each mock consisting of 60 simulations. The left panel of the figure displays the monopole of the momentum power spectrum for three different redshift ranges. The massE-based measurement well recovers the true momentum power spectrum except for a slight increment above $k\sim$ 0.1 $h$/Mpc as stated above. Although the survey volume of mock data is not large enough to measure the quadrupole of the momentum power spectrum, by examining the average of simulations it can give insights into whether the velocity field is recovered by the massE. Note that the quadrupole power spectrum is affected by the wide angle effect but for the SDSS MGS survey volume it is minimal compared to errors in the measurement (e.g. Castorina & White, 2018). Additionally, when comparing the true quadrupole to the one derived from massE, both measurements should be affected by the wide-angle effect to the same extent. As shown in the right panel of the figure, across the whole redshift range ($\Delta$$z$=0.02-0.20), the massE-based measurement follows the true quadrupole momentum with no systematic shift. For two redshift sub-ranges, only at high $k$ values exceeding 0.2 h/Mpc, the massE-based measurements appear to be higher than the true values, but note that the error is large too.

We perform additional tests, taking into account that the intrinsic massE (used to derive the velocity dispersion) does not show any redshift evolution but we mistakenly incorporate redshift evolution in the massE when calculating the PV (see Fig. 6 for two relations). The first case assumes the intercept in the second massE has a strong redshift evolution with a form of $S_{D_{0}}$(z)=$S_{D_{0}}$(1+z)⁵. As shown in Figure 11, the best fit to PVs of individual galaxies in three redshift bins still show good one-to-one correlations with slightly increase in errors. The upper panel of Figure 13 shows that both monopole and quadrupole of the momentum power spectrum over $\Delta{z}$=0.02-0.2 is well recovered too. The second case assumes the slope in the second massE has a redshift evolution in a form of $\beta$(z)=$\beta$(1+z)^0.4. At $z$=0.2 the evolution results in an approximately 20-$\sigma$ offset from our best estimate of $\beta$ as listed in Equation 4. Figure 12 shows that best fits to PVs of individual galaxies in three redshift bins also exhibit one-to-one correlations, but with larger errors. The lower panel of Figure 13 indicates the mean values of both the monopole and quadrupole momentum power spectra derived from massE still closely follow the intrinsic one.

Figure 14 (a) shows the result of the momentum power spectrum of our massE-based PV sample over a redshift range of $\Delta$$z_{\rm obs}$=0.02-0.2 where PV is measured with a redshift bin of $dz_{\rm obs}$=0.005. The covariance matrix is shown in Figure 15, which is based on the result from the mock. As compared to the true median of 31 Uchuu-SDSS mocks, our observed momentum power spectrum is slightly higher at low $k$ but the deviation becomes larger at $k$ $\gtrsim$ 0.07 $h$/Mpc. The difference could be the result that on small spatial scales non-linear effects affect the power spectrum, while the Uchuu-SDSS simulation is purely N-body (Ishiyama et al., 2021). In the figure, we overlay the result of TNG-300 hydrodynamics simulation that include baryons (Nelson et al., 2019). For this simulation, we use the snapshot at $z$=0 and include all sub-halos with stellar masses larger than 10¹⁰ $M_{\odot}$. The calculation of the momentum power spectrum is obtained through FFTPower implemented in nbodykit. Given a specific simulation box, we can have three lines of sight, each giving three independent measurements, as shown in the figure. The disparity between TNG and our observation above 0.07 $h$/Mpc is now much smaller. It remains unclear what causes the difference between our observations and Uchuu-SDSS mock data. The momentum power spectrum at high $k$ is dominated by the product of the density and velocity (e.g., see Fig. 13 in Howlett (2019)). The Uchuu-SDSS mock is designed to have the same galaxy bias as the observed one, a fact we confirm by comparing the density power spectra of our observations and the Uchuu-SDSS mock. Accordingly the difference at high-$k$ should not be attributed to different galaxy biases between the two. It is probable that the N-body simulation, such as Uchuu, does not account adequately for non-linear effects as compared with TNG hydrodynamic simulation. Fortunately, the constraint on the growth rate is not particularly sensitive to the maximum wavenumber $k_{\rm max}$.

Figure 14(b) and (c) presents two additional momentum power spectra for two sub-redshift ranges of ${\Delta}z_{\rm obs}$=0.10-0.20 and 0.02-0.10, respectively. Similar to the case of the full redshift range, the momentum power spectra within two sub-ranges exhibit a slightly elevation when compared to the Uchuu-SDSS result at low $k$. However, they show considerably larger disparities above $k$=0.07 $h$/Mpc, but tend to match more or less the result of the TNG hydrodynamic simulation at high $k$.

In this section, we discuss the result if varying some parameters that we use to derive the momentum power spectrum including the massE relation itself and the redshift bin $dz_{\rm obs}$.

In § 4.4, we confirm with mock data that the massE-based PV is insensitive to the accurate calibration of the massE relationship. Here we further demonstrate this with the observational data. In Figure 16(a), we overlay two momentum power spectra with the fiducial one: one is the case that the normalization $S_{D_{0}}$ of the massE relation evolves strongly with redshift as $S_{D_{0}}$ $\propto$ (1+$z$)⁵ , and another is the one that the slope $\beta$ of the relationship evolves strongly with the redshift as $\beta$ $\propto$ (1+$z$)^0.4. For the $\beta$, such a redshift evolution gives a value at $z$=0.2 that deviates from the observed one by 20-$\sigma$. In the first case, any alteration in the value of $S_{D_{0}}$ causes all distances to shift uniformly with the same amount, thereby maintaining the offset of a distance from the Gaussian mean. As a result, the momentum power spectrum remains the same regardless of the magnitude of the $S_{D_{0}}$ shift. In the second case, the momentum power spectrum changes a little bit but well within the error bar, which still benefits from our method. For galaxies within a narrow redshift bin, their physical properties including sizes, velocity dispersions and stellar masses exhibit a narrow range. If there is a change in the slope $\beta$, it only introduces a small amount of scatter in the distribution of distances. Therefore Figure 16(a) demonstrates that the momentum power spectrum with the massE-based PV is not sensitive to the accurate calibration of the massE relation. Figure 16(b) presents the momentum power spectra whose PV is measured with redshift bin sizes of $d{z_{\rm obs}}$=0.001 and $d{z_{\rm obs}}$=0.01, as compared to the fiducial one with $d{z_{\rm obs}}$=0.005. As shown in the figure, the three spectra are almost the same, indicating that the PV measurements is not sensitive to $d{z_{\rm obs}}$ which can vary by a factor of 10.

Before detailed modeling, we first roughly evaluate the statistical power of our sample in constraining a cosmological parameter, given that the signal to noise of the parameter constrained at a given $k$ is more or less proportional to the square root of the survey effective volume at that $k$ (e.g. Koda et al., 2014). Similar to the definition of the density power spectrum (Feldman et al., 1994; Tegmark, 1997), an effective volume for the momentum power spectrum can be written as

where $\bar{n}(z)$ and $\langle{v_{\rm p}^{2}}\rangle_{z}$ are the galaxy number density and PV variance as a function of the redshift, respectively. The above equation indicates that for high galaxy number densities $\bar{n}(z)$ $\gg$ $\langle{v_{\rm p}^{2}}\rangle_{z}$/$P^{\rm P}(k)$, the effective volume is just the geometrical volume and thus increases with the cube of the redshift. On the other hand, for low galaxy number densities, because $\langle{v_{\rm p}^{2}}\rangle_{z}$ roughly scales with the square of the redshift at $z$ < 0.2, the effective volume only scales linearly with the redshift.

To calculate the effective volume for our sample, we adopt the median of 31 Uchuu-SDSS mocks for $P^{\rm P}(k)$. The $\bar{n}(z)$ and $\langle{v_{\rm p}^{2}}\rangle_{z}$ are shown in Figure 4 for our sample. We derived $V_{\rm eff}(k)$ for the whole redshift range as well as two redshift sub-ranges as shown in Figure 17, along with 6dFGSv and SDSS-FP surveys. The 6dFGSv data has a spatial coverage of 17 000 deg² and a redshift range of 0.0001 to 0.0534 (Scrimgeour et al., 2016). We use $(cz/(1+z){\sigma_{\eta}})^{2}$ to get its velocity variance by setting $\sigma_{\eta}$=0.324 (Scrimgeour et al., 2016). The galaxy number density is almost a constant around 1.2$\times$10^-3 $h^{3}$ Mpc^-3. We use the same $P^{\rm P}(k)$ as for our sample. The SDSS FP-based PV sample covers a redshift range of 0.0033 to 0.1 and an area of 7016 deg² (Howlett et al., 2022). We use the above equation for the velocity variance with $\sigma_{\eta}$=0.23 (Howlett et al., 2022), measure $\bar{n}(z)$ from their catalog and also adopt the median of Uchuu-SDSS mocks for $P^{\rm P}(k)$.

The integral quantity $N_{\rm k,modes}$=$\int_{0.02}^{0.07}$$k^{2}V_{\rm eff}dk$ gives the total number of $k$ modes on linear scales where the growth rate is sensitive to. As compared to the 6dFGSv and SDSS FP-based PV catalogs, our sample offers a factor of about 2.4  and 6.7  increase in the above quantity, respectively. As shown in the figure, even for $z$<0.1, our effective volume is still larger than the SDSS-FP catalog, with about a factor of 1.7 larger in terms of $N_{\rm k,modes}$. This is because of our larger sample even at $z$ $<$ 0.1 as presented in § 2.3. Note that, although our sample extends to $z_{\rm obs}$=0.2, the galaxy number density above $z_{\rm obs}$=0.1 of the SDSS MGS drops rapidly.

We use the model that is based on the distribution function and Eulerian perturbation theory as developed by Vlah et al. (2012, 2013); Okumura et al. (2014) and Saito et al. (2014). The analytic formula of the model is compiled in the appendix of Howlett (2019). Here we give the basic one:

where

The term $P_{11}$, $P_{12}$, $P_{13}$ and $P_{22}$ can be calculated from the linear power spectrum while containing free parameters such as the growth rate $f$, the linear-scale bias $b_{\rm 1}$, the non-linear bias $b_{\rm 2}$, the non-linear velocity dispersion $\sigma_{v}$ etc. Through these parameters the model accounts for the effect of small-scale non-linear motion on the power spectrum.

Before running the MCMC fitting, we first quantify the effect of the window function for the momentum power spectrum measurement. Due to the survey volume as well as the sample selection as a function of both spatial positions and redshift, the power spectrum suffers from the effect of the window function. We evaluate this effect using the python code pypower.CatalogSmoothWindow³³3https://github.com/cosmodesi/pypower with following steps: (1) a random elliptical catalog: we start with a random catalog available for the SDSS LSS sample⁴⁴4http://sdss.physics.nyu.edu/vagc/lss.html, http://sdss.physics.nyu.edu/lss/dr72/ (Blanton et al., 2005). We first follow the observed fraction of ellipticals in the LSS sample to randomly pick objects from the above catalog to define the random elliptical catalog. By fitting the observed redshift distribution of the MGS-elliptical sample with multiple Gaussian profiles, we assign a random redshift to those in the random elliptical catalog. After that, we remain those in the North Cap region and limit the redshift range to 0.02-0.2. (2) the window function: based on the above random catalog, we then use pypower.CatalogSmoothWindow to derive multipoles of the window function up to $\ell$=9 with the wide angle up to the order of $n$=1. The result in the configuration space is shown in Figure 18(a) for ${\Delta}z$=0.02-0.20. (3) the theory power spectrum: we then use Equation 21 to produce theoretical momentum power spectrum at $\ell$=0, 2 & 4 with $k$ from 0.001 to 10 with a resolution of 0.005. (4) the power spectrum with window function and wide angle: we then convolve the above theoretical momentum power spectrum with the window function and wide-angle, and show its ratio to the pure theoretical spectrum in Figure 18(b). It is found that the effect of the window function and wide-angle is small but a difference of a few to ten percent is seen for $k$ values between 0.02 and 0.1 $h$/Mpc. We thus carry out the MCMC fitting by incorporating correction factors into the theoretical power spectrum. The correction factors for two redshift sub-ranges of ${\Delta}z$=0.02-0.10 and ${\Delta}z$=0.10-0.20 are also produced following the same procedure, and their differences from the above one are small.

As discussed in Qin et al. (2019), the distribution of the momentum power spectrum at given $k$ is not always normal as shown in Fig. 19(a), and it is thus suggested to first Box-Cox transform (Box & Cox, 1964) the momentum power spectrum to $Z$ for which the MCMC fitting is carried out:

Here, because of large errors in $P(k)$ some P($k$) is negative so that we add a $P_{\rm shift}$ which is set to be the minimum negative value of the distribution. The code scipy.stats.boxcox is used to find the best $\lambda$ for each transformation.

We set the background cosmological parameter to the Planck 2018 result with $\sigma_{8}$($z$=0)=0.811 (Planck Collaboration et al., 2020), based on which the linear power spectrum at $z$=0 is calculated as an input of the theoretical power spectrum through nbodykit using the CLASS transfer function. In this case, if we multiply the derived $f$ with $\sigma_{8}$($z$=0), it gives $f\sigma_{8}$ at $z=z_{\rm eff}$. Following Qin et al. (2019), we set the following four free parameters: $f$, $b_{\rm 1}$, $b_{\rm 2}$ and $\sigma_{v}$. The MCMC fitting has been carried out through Python code emcee (Foreman-Mackey et al., 2013). Priors of four parameters are listed in Table 2.

We set $k_{\rm min}=0.02$ $h$/Mpc and conduct the MCMC fitting for $k_{\rm max}$ values of 0.15, 0.20 and 0.25 $h$/Mpc, respectively, for galaxies within $\Delta$z=0.02-0.2. The corner plots of three fitting results are shown in Figure 20. The growth rate $f\sigma_{8}$ remains unaffected by the non-linear parameters $b_{2}$ and $\sigma_{v}$. While there are some correlations between $f\sigma_{8}$ and $b_{1}$ at a given $k_{\rm max}$, the best-fit values of $f\sigma_{8}$ exhibit no sensitivity to $k_{\rm max}$. This is consistent with the expected advantage of the momentum power spectrum that is affected little by the bias and non-linear effects. Figure 22 offers a further comparison in $f\sigma_{8}$ between observations and Uchuu-SDSS mocks, revealing that our observations and the mock data conform to the predictions of General Relativity (GR) and display no trends with respect to $k_{\rm max}$. The errors from observation are slightly larger than those from mocks. This is reasonable because the median momentum power spectrum of mocks at each $k$ almost precisely lies on the theoretical prediction. As shown in Figure 23, the best fit is slightly below the median of data points, which is caused by the fact that the covariance matrix contains significant non-diagonal parts. This is also the reason that the best-fitted $f\sigma_{8}$ of mocks are not exactly but slightly smaller than the GR’s prediction.

As shown in Figure 20 and Figure 21, there are noticeable discrepancies in the best-fit $b_{1}$ between observations and mocks. This discrepancy is indicative of a disparity in the momentum power spectrum at high $k$ regime between the two as shown in Figure 14. While both observations and mocks have a similar galaxy bias ($\sim$ 1.5) as seen by their density power spectrum, the difference in the above result is likely due to the velocity field on small scales. As stated above, the comparison between Uchuu-SDSS and TNG indicates that Uchuu-SDSS as a N-body simulation may not fully account for non-linear effects. In addition, the observations exhibit higher values than TNG above $k$= 0.1 $h$/Mpc by about 20-30%. There may be two reasons for this: one is that there is some random zero-point offsets in velocity similar to those seen in mock simulations (see Figure 9), which adds a constant to the shot noise (Howlett, 2019). This addition predominantly affects the high-$k$ regime but becomes negligible at low-$k$; another is that since the mocks underestimate the power at high-$k$, the corresponding error at high-$k$ is also underestimated. Nevertheless, these factors should not affect our conclusion regarding $f\sigma_{8}$, as it primarily relies on the low $k$ regime.

We also carry out fit to the momentum power spectra of two sub-redshift ranges of ${\Delta}z_{\rm obs}$=0.02-0.1 and 0.1-0.2, which gives $f\sigma_{8}$ of 0.416${}^{+0.074}_{-0.076}$  and 0.526${}^{+0.133}_{-0.148}$, respectively, for $k_{\rm max}$=0.20 $h$/Mpc. As shown in Figure 23, they are consistent with the result of the entire redshift range, which yields a $f\sigma_{8}$ value of 0.459${}^{+0.068}_{-0.069}$, well within the 1-$\sigma$ confidence interval. The inverse variance of $f\sigma_{8}$ for the whole redshift range is almost the same as the sum of those of two redshift sub-ranges too. These results offer additional evidence for the robustness of our results.

Table 2 lists $\chi^{2}$/d.o.f. values for three redshift ranges, which are between 2 and 3. We suspect that this is due to the fact that the mock data under-estimate the momentum power spectra and associated errors at high $k$ values as stated in § 5.1. To test this for $\Delta{z}$=0.02-0.2, we increase the error by a factor of 1.5 for $k$ above 0.1 $h$/Mpc and corresponding non-diagonal parts in the covariance matrix. By carrying out the MCMC fitting again, it is found that the best-fit $f\sigma_{8}$ remains unchanged but now $\chi^{2}$/d.o.f. is 36/32. This further demonstrates that the growth rate measured from momentum power spectrum is insensitive to momentum power spectra on small scales.

Qin et al. (2019) adopted the same theoretical model with the same set of free parameters for the 6dFGSv and 2MTF PV data. They obtained $f\sigma_{8}$=0.404${}^{+0.082}_{-0.081}$ with the combined density power spectrum and momentum power spectrum of two data-sets together.This constraint is similar to that by Turner et al. (2023) using the combined velocity and density data of the 6dFGSv survey. If adopting the momentum power spectrum of the 6dFGSv survey alone, Qin et al. (2019) estimated 68% confidence range to be 0.226 (see their Table 2). Our measurement from ${\Delta}z$=$0.02-0.20$ offers a factor of 1.7 higher accuracy. This improvement in accuracy is somewhat less than that expected by the square root of the increase in the effective volume ($\times$6.7), but are consistent within 30%. This is reasonable given different mocks or fitting techniques (Eisenstein et al., 2005). Lai et al. (2023) measured the growth rate using the maximum-likelihood fields method through both density and PV data of the SDSS-FP survey (Howlett et al., 2022). The derived growth rate is 0.405${}^{+0.076}_{-0.071}$(stat)$\pm$0.009 (sys), whose confidence range is consistent with our low-$z$ one (0.416${}^{+0.074}_{-0.076}$  for ${\Delta}z$=0.02-0.1). They did not report the result when only using PV data. However, as mentioned earlier, although both our catalog and the SDSS-FP catalog are based on the same parent sample, our massE-based PV catalog includes 2.5  more objects than the SDSS-FP catalog at $z$ < 0.1. This could account for our PV-only result exhibiting a similar level of error as theirs of combined PV and density data. Saulder et al. (2023) have provided the predicted number density and PV error for the DESI survey, from which they estimated a 68% confidence range for $f\sigma_{8}$, which is 0.126 and 0.092 for $k_{\rm max}$ values of 0.1 $h$/Mpc and 0.2 $h$/Mpc, respectively. We measure their effective volume and find that their number of $k$-modes between 0.02 and 0.07 $h$/Mpc is 1.7 times more than our entire redshift range. This corresponds to a 1.3-fold enhancement in their precision, consistent with our constraint of $f\sigma_{8}$ at a value of 0.459${}^{+0.068}_{-0.069}$  at $k_{\rm max}$=0.2 $h$/Mpc.

Howlett et al. (2015) measured the redshift space distortion of the SDSS MGS sample through two-point correlation function, resulting in a derived $f\sigma_{8}$ of 0.49${}^{+0.15}_{-0.14}$. Our massE-based PV measurement of the same sample thus offers a factor of 2.1  improvement in the error of $f\sigma_{8}$. In Figure 24, we overlay our measurement on most recent collection of $f\sigma_{8}$ that has been measured with redshift space distortion (Alam et al., 2021).

Our work has demonstrated that the massE can serve as a new cosmic ruler with only two nuisance parameters, allowing for the probing of the PV beyond the current redshift limit of 0.1. For future perspective, we envision that there is no apparent limit for the PV measurement with the massE to higher redshift, as it is insensitive to the accurate calibration of the massE relation itself. In addition, because of $\kappa$ term in Equation 10 that is used to convert the distance measurement to PV, the PV variance scales with the redshift slower than $z/(1+z)$ by a few times, especially above $z$=0.2, as shown in Figure 25(a). As a result, although measuring the momentum power spectrum still requires a high galaxy number density above $z$=0.2, it has a potential to provide strong constraints on the growth rate in the late universe where dark energy dominates. To illustrate this quantitatively, we use the Fisher-matrix forecast code⁵⁵5https://github.com/CullanHowlett/PV_fisher (Howlett et al., 2017; da Cunha et al., 2017) with a survey area of 14 k square degree, a distance error of $\sigma_{\eta}$=20%, a galaxy number density of 8$\times$10^-3 $h^{3}$/Mpc³ and Equation 10 for the PV error. As shown in Figure 25(a), PV-only data yields stronger constraints compared to the RSD only for this specific data-set where $k_{\rm max}$ is set to be 0.2 $h$/Mpc. On the other hand, the momentum power spectrum mainly relies on linear scales to extract cosmological information, complementing the RSD. The combination of these two offers significant improvement in the growth rate, potentially exceeding the DESI’s constraints by a factor of 1.7 to 4.7 in terms of errors.