The Redshift-Space Momentum Power Spectrum I: Optimal Estimation From Peculiar Velocity Surveys

TL;DR

This paper develops a comprehensive framework for extracting the redshift-space momentum power spectrum from galaxy peculiar velocity surveys, improving cosmological parameter constraints by including momentum information.

Contribution

It introduces a new estimator for the redshift-space momentum power spectrum, accounting for shot noise, survey geometry, and measurement errors, with a proof-of-concept demonstration.

Findings

Successfully measured and fitted galaxy density and momentum power spectra from mock surveys.

Including momentum power spectrum improves constraints on the growth rate of structure.

Provided theoretical models for non-linear redshift-space power spectrum multipoles.

Abstract

Low redshift surveys of galaxy peculiar velocities provide a wealth of cosmological information. We revisit the idea of extracting this information by directly measuring the redshift-space momentum power spectrum from such surveys. We provide a comprehensive theoretical and practical framework for estimating and fitting this from data, analogous to well understood techniques used to measure the galaxy density power spectrum from redshift surveys. We formally derive a new estimator, which includes the effects of shot noise and survey geometry; we evaluate the variance of the estimator in the Gaussian regime; we compute the optimal weights for the estimator; we demonstrate that the measurements are Gaussian distributed, allowing for easy extraction of cosmological parameters; and we explore the effects of peculiar velocity measurement errors. We finish with a proof-of-concept using realistic mock galaxy catalogues, which demonstrates that we can measure and fit both the redshift-space galaxy density and momentum power spectra from peculiar velocity surveys and that including the latter substantially improves our constraints on the growth rate of structure. We also provide theoretical descriptions for modelling the non-linear redshift-space density and momentum power spectrum multipoles, and forecasting the constraints on cosmological parameters using the Fisher information contained in these measurements for arbitrary weights. These may be useful for measurements of the galaxy density power spectrum even in the absence of peculiar velocities.

Figures

Tables

Table 1. Table 1: Constraints and forecasts on the growth rate of structure from fitting various multipoles of the density and momentum power spectra from mock galaxy catalogues. We use the 16 th superscript 16 th 16^{\mathrm{th}} , 50 th superscript 50 th 50^{\mathrm{th}} and 84 th superscript 84 th 84^{\mathrm{th}} percentiles to quote the median and errors on the marginalised f ​ σ 𝑓 𝜎 f\sigma measurements. Given the input cosmology of the mocks, the expected value is f ​ σ 8 ≈ 0.432 𝑓 subscript 𝜎 8 0.432 f\sigma_{8}\approx 0.432 . We also show forecasts using the Fisher matrix method in Appendix  B which quantifies the lower error bound using the information in the velocity or density fields for different weighting schemes. We provide forecasts for the optimal weighting scheme derived in Section  5 , and where we use the common technique of using fixed values for the power spectrum in these weights. The latter case does reduce the growth rate information in both the density and momentum power spectra, but the latter is affected more greatly due to its 1 / k 2 1 superscript 𝑘 2 1/k^{2} dependence on linear scales.

Measurements	$f{\sigma }_8$ (Percentiles)	$\sigma \left(f{\sigma }_8\right)$ Predictions
		Optimum weights	Constant $P_{FKP}$ weights
Momentum $\ell =0$	${0.411}_{-0.064}^{+0.056}$	$0.032$	$0.043$
Momentum $\ell =0,2$	${0.410}_{-0.059}^{+0.053}$	$0.031$	$0.039$
Density $\ell =0,2$	${0.371}_{-0.070}^{+0.081}$	$0.049$	$0.051$
Density $\ell =0,2,4$	${0.392}_{-0.063}^{+0.062}$	$0.048$	$0.051$
Density $\ell =0$ + Momentum $\ell =0$	${0.430}_{-0.054}^{+0.056}$	$0.031$	$0.041$
Density $\ell =0,2$ + Momentum $\ell =0$	${0.383}_{-0.038}^{+0.038}$	$0.028$	$0.035$
Density $\ell =0,2,4$ + Momentum $\ell =0$	${0.390}_{-0.039}^{+0.039}$	$0.028$	$0.035$
Density $\ell =0$ + Momentum $\ell =0,2$	${0.345}_{-0.044}^{+0.040}$	$0.030$	$0.037$
Density $\ell =0,2$ + Momentum $\ell =0,2$	${0.364}_{-0.035}^{+0.034}$	$0.027$	$0.033$
Density $\ell =0,2,4$ + Momentum $\ell =0,2$	${0.367}_{-0.034}^{+0.032}$	$0.027$	$0.033$

Equations

$$\mathbf{\rho }(\mathbf{r})=(1+{\delta }_g(\mathbf{r}))\mathbf{v}(\mathbf{r}),$$

$$(2\pi {)}^3{\delta }^D(\mathbf{k}-{\mathbf{k}}^{\prime })P^p(\mathbf{k})$$

$$=⟨u(\mathbf{k})u({\mathbf{k}}^{\prime })⟩+⟨u(\mathbf{k}){\delta }_g({\mathbf{k}}^{\prime })u({\mathbf{k}}^{\prime })⟩+⟨{\delta }_g(\mathbf{k})u(\mathbf{k})u({\mathbf{k}}^{\prime })⟩+⟨{\delta }_g(\mathbf{k})u(\mathbf{k}){\delta }_g({\mathbf{k}}^{\prime })u({\mathbf{k}}^{\prime })⟩,$$

$$\xi (\mathbf{r})=\xi (r,\mu )=\sum _{\ell }{\xi }_{\ell }(r)L_{\ell }(\mu ),$$

$$F^p(\mathbf{r})=\frac{w(\mathbf{r})n_g(\mathbf{r})u(\mathbf{r})}{A},$$

$$F^{\delta }(\mathbf{r})=w(\mathbf{r})[n_g(\mathbf{r})-\alpha n_s(\mathbf{r})]/A,$$

$$⟨F^p(\mathbf{r})F^p({\mathbf{r}}^{\prime })⟩=\frac{1}{A^2}w(\mathbf{r})w({\mathbf{r}}^{\prime })⟨n_g(\mathbf{r})n_g({\mathbf{r}}^{\prime })u(\mathbf{r})u({\mathbf{r}}^{\prime })⟩.$$

$$⟨n_g(\mathbf{r})n_g({\mathbf{r}}^{\prime })u(\mathbf{r})u({\mathbf{r}}^{\prime })⟩=\bar{n}(\mathbf{r})\bar{n}({\mathbf{r}}^{\prime }){\xi }^p(\mathbf{r}-{\mathbf{r}}^{\prime })+\bar{n}(\mathbf{r})⟨v^2(\mathbf{r})⟩{\delta }^D(\mathbf{r}-{\mathbf{r}}^{\prime }),$$

$$\begin{gathered} ⟨F^p(\mathbf{r})F^p({\mathbf{r}}^{\prime })⟩=\frac{1}{A^2}[w(\mathbf{r})w({\mathbf{r}}^{\prime })\bar{n}(\mathbf{r})\bar{n}({\mathbf{r}}^{\prime }){\xi }^p(\mathbf{r}-{\mathbf{r}}^{\prime }) \\ +w^2(\mathbf{r})\bar{n}(\mathbf{r})⟨v^2(\mathbf{r})⟩{\delta }^D(\mathbf{r}-{\mathbf{r}}^{\prime })], \end{gathered}$$

$$P^p(\mathbf{k})=\sum _{\ell }{P}_{\ell }^p(k)L_{\ell }(\mu ),$$

$$|F_{\ell }^p(k){|}^2=\frac{2\ell +1}{V}\int \frac{d{\Omega }_k}{4\pi }\int d^{3}r\int d^3{r}^{\prime }{F}^p(\mathbf{r})F^p({\mathbf{r}}^{\prime })L_{\ell }(\hat{\mathbf{k}}\cdot {\hat{\mathbf{r}}}^{\prime })e^{i\mathbf{k}\cdot (\mathbf{r}-{\mathbf{r}}^{\prime })},$$

$$⟨|F_{\ell }^p(k){|}^2⟩$$

$$\displaystyle=\frac{2\ell+1}{A^{2}}\int\frac{d\Omega_{k}}{4\pi}\biggl{[}\frac{1}{V}\int d^{3}r\int d^{3}r^{\prime}w(\boldsymbol{r})w(\boldsymbol{r}^{\prime})\bar{n}(\boldsymbol{r})\bar{n}(\boldsymbol{r}^{\prime})\int\frac{d^{3}k^{\prime}}{(2\pi)^{3}}\,P^{p}(\boldsymbol{k}^{\prime},\boldsymbol{r}^{\prime})L_{\ell}(\hat{\boldsymbol{k}}\cdot\hat{\boldsymbol{r}}^{\prime})e^{i(\boldsymbol{k}-\boldsymbol{k}^{\prime})\cdot(\boldsymbol{r}-\boldsymbol{r}^{\prime})}+N^{p}_{\ell}(\boldsymbol{k})\biggl{]},$$

$$N_{\ell }^p(\mathbf{k})=\int d^{3}r{w}^2(\mathbf{r})\bar{n}(\mathbf{r})⟨v^2(\mathbf{r})⟩L_{\ell }(\hat{\mathbf{k}}\cdot \hat{\mathbf{r}}).$$

$$\langle|F^{p}_{\ell}(k)|^{2}\rangle=\frac{(2\ell+1)}{A^{2}}\int\frac{d\Omega_{k}}{4\pi}\biggl{[}\sum_{\ell^{\prime}}\int\frac{d^{3}k^{\prime}}{(2\pi)^{3}}P^{p}_{\ell^{\prime}}(k^{\prime})G(\boldsymbol{k}-\boldsymbol{k}^{\prime})S_{\ell,\ell^{\prime}}(\boldsymbol{k},\boldsymbol{k}^{\prime})+N^{p}_{\ell}(\boldsymbol{k})\biggl{]},$$

$$G(\mathbf{k}-{\mathbf{k}}^{\prime })$$

$$S_{\ell ,{\ell }^{\prime }}^{\ast }(\mathbf{k},{\mathbf{k}}^{\prime })$$

$$\widehat{P_{\ell }^{\delta }}(k)=|F_{\ell }^{\delta }(k){|}^2-N_{\ell }^{\delta }(k),$$

$$\widehat{P_{\ell }^p}(k)=|F_{\ell }^p(k){|}^2-N_{\ell }^p(k).$$

$$⟨v^2(\mathbf{r})⟩=\frac{1}{N}\sum _N{u}^2(\mathbf{r}).$$

$$P_L^p(k,\mu )\approx {\mu }^2{f}^2(aH{)}^2{k}^{-2}{P}_L(k)$$

$$P_{L,0}^p(k)=\frac{f^2(aH{)}^2}{3k^2}{P}_L(k),P_{L,2}^p(k)=\frac{2f^2(aH{)}^2}{3k^2}{P}_L(k),P_{L,\ell }^p(k)=0\mathrm{for}\ell \ne 0,2.$$

$$C_{\ell {\ell }^{\prime }}^P(k,k^{\prime })=⟨|F_{\ell }^p(\mathbf{k}){|}^2|F_{{\ell }^{\prime }}^p({\mathbf{k}}^{\prime }){|}^2⟩-⟨|F_{\ell }^p(\mathbf{k}){|}^2⟩⟨|F_{{\ell }^{\prime }}^p({\mathbf{k}}^{\prime }){|}^2⟩.$$

$$\begin{gathered} C_{\ell {\ell }^{\prime }}^P(k,k^{\prime })=\frac{(2\ell +1)(2{\ell }^{\prime }+1)}{A^4{V}^2}\int \frac{d{\Omega }_k}{4\pi }\int \frac{d{\Omega }_{k^{\prime }}}{4\pi }\int d^3{r}_1\int d^3{r}_2\int d^3{r}_3\int d^3{r}_{4}w({\mathbf{r}}_1)w({\mathbf{r}}_2)w({\mathbf{r}}_3)w({\mathbf{r}}_4)L_{\ell }(\hat{\mathbf{k}}\cdot \widehat{{\mathbf{r}}_2})L_{{\ell }^{\prime }}(\widehat{{\mathbf{k}}^{\prime }}\cdot \widehat{{\mathbf{r}}_4}) \\ [⟨n_g({\mathbf{r}}_1)n_g({\mathbf{r}}_2)n_g({\mathbf{r}}_3)n_g({\mathbf{r}}_4)u({\mathbf{r}}_1)u({\mathbf{r}}_2)u({\mathbf{r}}_3)u({\mathbf{r}}_4)⟩-⟨n_g({\mathbf{r}}_1)n_g({\mathbf{r}}_2)u({\mathbf{r}}_1)u({\mathbf{r}}_2)⟩⟨n_g({\mathbf{r}}_3)n_g({\mathbf{r}}_4)u({\mathbf{r}}_3)u({\mathbf{r}}_4)⟩]e^{i\mathbf{k}\cdot ({\mathbf{r}}_1-{\mathbf{r}}_2)}{e}^{i\mathbf{k}\cdot ({\mathbf{r}}_3-{\mathbf{r}}_4)}. \end{gathered}$$

$$⟨n_g({\mathbf{r}}_1)n_g({\mathbf{r}}_2)n_g({\mathbf{r}}_3)n_g({\mathbf{r}}_4)u({\mathbf{r}}_1)u({\mathbf{r}}_2)u({\mathbf{r}}_3)u({\mathbf{r}}_4)⟩=⟨n_g({\mathbf{r}}_1)n_g({\mathbf{r}}_2)u({\mathbf{r}}_1)u({\mathbf{r}}_2)⟩⟨n_g({\mathbf{r}}_3)n_g({\mathbf{r}}_4)u({\mathbf{r}}_3)u({\mathbf{r}}_4)⟩+2cyc.,$$

$$C^{P}_{\ell\ell^{\prime}}(k,k^{\prime})=\frac{(2\ell+1)(2\ell^{\prime}+1)}{A^{4}V^{2}}\int\frac{d\Omega_{k}}{4\pi}\int\frac{d\Omega_{k^{\prime}}}{4\pi}\biggl{[}Q_{\ell\ell^{\prime}}(\boldsymbol{k},\boldsymbol{k}^{\prime})Q^{*}_{00}(\boldsymbol{k},\boldsymbol{k}^{\prime})+Q_{0\ell^{\prime}}(\boldsymbol{k},\boldsymbol{k}^{\prime})Q^{*}_{\ell 0}(\boldsymbol{k},\boldsymbol{k}^{\prime})\biggl{]},$$

$$Q_{\ell {\ell }^{\prime }}(\mathbf{k},{\mathbf{k}}^{\prime })$$

$$=\int d^{3}r\int d^3{r}^{\prime }w(\mathbf{r})w({\mathbf{r}}^{\prime })\bar{n}(\mathbf{r})\bar{n}({\mathbf{r}}^{\prime })\int \frac{d^3{k}^{\prime \prime }}{(2\pi {)}^3}{P}^p({\mathbf{k}}^{\prime \prime },{\mathbf{r}}^{\prime })L_{\ell }(\hat{\mathbf{k}}\cdot \hat{\mathbf{r}})L_{{\ell }^{\prime }}({\hat{\mathbf{k}}}^{\prime }\cdot \widehat{{\mathbf{r}}^{\prime }})e^{i(\mathbf{k}\cdot \mathbf{r}-{\mathbf{k}}^{\prime }\cdot {\mathbf{r}}^{\prime }+{\mathbf{k}}^{\prime \prime }\cdot (\mathbf{r}-{\mathbf{r}}^{\prime }))}+M_{\ell {\ell }^{\prime }}^p(\mathbf{k},{\mathbf{k}}^{\prime }),$$

$$M_{\ell {\ell }^{\prime }}^p(\mathbf{k},{\mathbf{k}}^{\prime })=\int d^3\mathbf{r}{w}^2(\mathbf{r})\bar{n}(\mathbf{r})⟨v^2(\mathbf{r})⟩L_{\ell }(\hat{\mathbf{k}}\cdot \hat{\mathbf{r}})L_{{\ell }^{\prime }}({\hat{\mathbf{k}}}^{\prime }\cdot \widehat{{\mathbf{r}}^{\prime }})e^{i(\mathbf{k}-{\mathbf{k}}^{\prime })\cdot \mathbf{r}}.$$

$$P^{p}_{\mathrm{eff}}(\boldsymbol{k},\boldsymbol{k}^{\prime},\boldsymbol{r})=\frac{1}{2}\biggl{[}P^{p}(\boldsymbol{k},\boldsymbol{r})+P^{p}(\boldsymbol{k}^{\prime},\boldsymbol{r})\biggl{]}=\frac{1}{2}\sum_{\ell^{\prime\prime}}\biggl{[}P^{p}_{\ell^{\prime\prime}}(k)L_{\ell^{\prime\prime}}(\hat{\boldsymbol{k}}\cdot\hat{\boldsymbol{r}})+P^{p}_{\ell^{\prime\prime}}(k^{\prime})L_{\ell^{\prime\prime}}(\hat{\boldsymbol{k}}^{\prime}\cdot\hat{\boldsymbol{r}})\biggl{]}.$$

$$Q_{\ell {\ell }^{\prime }}(\mathbf{k},{\mathbf{k}}^{\prime })=\int d^{3}r{w}^2(\mathbf{r}){\bar{n}}^2(\mathbf{r})P_{\mathrm{eff}}^p(\mathbf{k},{\mathbf{k}}^{\prime },\mathbf{r})L_{\ell }(\hat{\mathbf{k}}\cdot \hat{\mathbf{r}})L_{{\ell }^{\prime }}({\hat{\mathbf{k}}}^{\prime }\cdot \hat{\mathbf{r}})e^{i(\mathbf{k}-{\mathbf{k}}^{\prime })\cdot \mathbf{r}}+M_{\ell {\ell }^{\prime }}^p(\mathbf{k},{\mathbf{k}}^{\prime }).$$

$$C^{p}_{\ell\ell^{\prime}}(k,k^{\prime})=\frac{(2\pi)^{3}}{V_{k}V}(2\ell+1)(2\ell^{\prime}+1)\delta^{D}(k-k^{\prime})\int_{-1}^{1}d\mu\,L_{\ell}(\mu)L_{\ell^{\prime}}(\mu)\biggl{[}P^{p}(k,\mu)+\frac{\langle v^{2}\rangle}{\bar{n}}\biggl{]}^{2},$$

Full text

††pagerange: The Redshift-Space Momentum Power Spectrum I: Optimal Estimation From Peculiar Velocity Surveys–LABEL:lastpage††pubyear: 2019

1 Introduction

The nature of dark matter and dark energy is one of the greatest mysteries in cosmology. Modification to our current theory of gravity, General Relativity (GR; Einstein 1916) offers a way to solve this. Large galaxy surveys, through the use of Redshift Space Distortion (RSD; Kaiser 1984), give us a way to test different theories of gravity on scales dominated by dark matter and dark energy. RSD induces an anisotropy in the clustering of galaxies proportional to the rate at which structure grows, commonly parameterised on linear scales by the parameter $f$. GR provides stringent predictions on this ‘growth rate of structure’, namely that it is scale-independent and evolves with redshift as $f(z)\approx\Omega_{m}(z)^{0.55}$.

Many galaxy surveys have been used to place constraints on the growth rate of structure (e.g., Blake et al. 2011; Beutler et al. 2012; de la Torre et al. 2013; Howlett et al. 2015b; Alam et al. 2017; Gil-Marín et al. 2018, to name a few), and these will tighten significantly with the advent of next generation surveys such as DESI (Levi et al., 2013), 4MOST (de Jong et al., 2012) and Euclid (Laureijs et al., 2011). In addition, there exists a long history of works (e.g., Strauss 1989; Kaiser et al. 1991; Nusser & Davis 1994; Hudson et al. 1995; Park 2000; Silberman et al. 2001; Pike & Hudson 2005; Carrick et al. 2015 or see Strauss & Willick (1995) for an excellent and comprehensive review) which demonstrate that combining peculiar velocity (PV) surveys with galaxy redshifts offers a way to test our cosmological model. More recent studies (Burkey & Taylor, 2004; Koda et al., 2014; Adams & Blake, 2017; Howlett et al., 2017a, c) have emphasised that growth rate measurements from future redshift surveys will be fundamentally limited by cosmic variance - the volume of the universe, and hence the number of independent ‘modes’ with which we can measure cosmology, is too small - but have also repeatedly demonstrated that the addition of PV data is one of the most promising routes to overcome this fundamental limitation. The use of PV data acts in the same way as the multi-tracer method (McDonald & Seljak, 2009), where one tracer has a known large scale bias of one (Desjacques & Sheth, 2010; Jennings et al., 2015; Chen et al., 2018). Multi-tracer analyses allow us to improve our constraints on the growth rate by including information from the relative clustering amplitudes of multiple tracers of the same underlying density field (Abramo & Leonard, 2013). In the case of PV surveys, this is the density field traced separately by the distribution of luminous matter, and the motions of galaxies induced through gravity.

Realising the potential of PV surveys for constraining gravity requires development of techniques to measure the correlations in the velocity field and between the density and velocity fields. Kolatt & Dekel (1997) estimated the power spectrum by first reconstructing the density field on a grid from a set of PV measurements. Many studies (Zaroubi et al., 1997; Macaulay et al., 2012; Johnson et al., 2014; Howlett et al., 2017b; Huterer et al., 2017) have also successfully used a likelihood based method to estimate the velocity power spectrum from existing datasets. This was extended in Adams & Blake (2017) to include measurements of the galaxy density power spectrum, and the cross-power spectrum between the density and velocity fields. This latter approach, in theory, fully captures the information in the two fields. However, it requires a computation of the integral of the power spectrum for each individual pair of galaxies, and each choice of model. Next generation datasets such as the Taipan Galaxy Survey (da Cunha et al., 2017), WALLABY survey (Koribalski, 2012), or supernovae (Howlett et al., 2017c; Kim et al., 2019) will contain up to an order of magnitude as many PV measurements as all current surveys combined (Tully et al., 2016), and require more complex non-linear modelling. This may make the above approaches prohibitively expensive. An alternative is measuring the velocity correlation function (Gorski et al., 1989; Groth et al., 1989; Juszkiewicz et al., 2000) or the density-velocity cross-correlation (Nusser, 2017) directly, but recent results (Wang et al., 2018) suggest the presence of non-Gaussian posteriors in such measurements makes extracting cosmological information difficult.

In this work, we present an alternative; estimation of the momentum (or mass-weighted velocity) power spectrum from a real survey. Much like the galaxy density power spectrum and redshift surveys, this presents an efficient compression of the information contained in a PV survey. This is not the first work to use this; measurements of the momentum power spectrum are common when analysing simulations (Scoccimarro, 2004; Pueblas & Scoccimarro, 2009; Jennings et al., 2011; Chen et al., 2018) although often the underlying goal is analysis of the non-linear velocity field, which cannot be easily extracted from such a measurement. This has also been applied to data. Two key references that present this as a viable method for extracting cosmology from PV surveys are Park (2000) and Park & Park (2006). However, to the best of the author’s knowledge this has not been revisited over the last decade, whilst techniques to measure cosmology from RSD have improved tremendously. Recently, PV surveys have also seen a resurgence, with significant improvements in data quality and volume that will only continue over the coming years. As such, we feel it is time to revisit this, and we completely re-interpret the formalism in the framework of modern galaxy surveys, bringing this technique up to the same standard as is used for modern RSD measurements. As we will show, this approach overcomes some of the short comings with the above methods, although does come with its own challenges (in particular estimation of the covariance matrix, as with the galaxy density power spectrum, is theoretically difficult and computationally demanding).

With this work, we aim to present the momentum power spectrum as a viable method to extract the full cosmological information from PV surveys and provide a complete description for how this can be measured in current or future surveys. This paper has the following structure: In Sections 2 and 3 we introduce the momentum field and present a new estimator for the redshift-space momentum power spectrum based on the well-known Feldman et al. 1994 (FKP) and Yamamoto et al. (2006) estimators for the galaxy density power spectrum. We then measure this in simulated halo catalogues and identify the similarities and differences between the two power spectra; i.e., the similar normalisations and window function convolutions. In Sections 4 and 5, we derive the covariance of our estimator in the Gaussian regime, compare this to the same halo simulations and derive the optimal weights for our estimator. Throughout these three Sections there is, by design, considerable analogy between the momentum and galaxy density power spectra. In Section 6, we then turn to realistic mock PV catalogues to demonstrate how our estimator works in reality, including several quirks of PV surveys that do not affect galaxy redshift surveys such as PV measurement errors. In particular, in Section 7 we bring everything together and demonstrate that we can use combined measurements of the galaxy density power spectrum and momentum power spectrum multipoles to improve our constraints on the growth rate of structure. We also compare our results to forecasts for the same mock data to show that we are recovering the expected amount of information on the growth rate. We conclude in Section 8. Finally, Appendices A and B provide theoretical descriptions of how to compute the non-linear redshift-space density and momentum power spectrum, and forecast the constraints on cosmological parameters using the Fisher information contained in the multipoles of the galaxy density and velocity/momentum power spectrum for arbitrary weights. These may be useful for measurements of the galaxy density power spectrum even in the absence of peculiar velocities.

This is just the first in a planned series of works using the momentum power spectrum. In a follow up paper in this series (Qin et. al, in preparation; Paper II) we apply this to the combined 2MTF (Hong et al. 2014; Hong et al., in preparation) and 6dFGSv (Campbell et al., 2014) data and show that we can recover accurate and robust constraints on the growth rate of structure at low redshift. We also plan to investigate how to extract the additional information in the cross-power spectrum of the density and momentum fields (which as shown in Section 7 is expected to improve the growth rate constraints by at least a further $\sim 10\%$), and to investigate potential constraints on alternative models of gravity that could be achieved with our method.

2 The Momentum Field

We start by defining the momentum field, following Park (2000), as the velocity field traced by discrete galaxies. Mathematically, we write the momentum field, $\boldsymbol{\rho}$ as

[TABLE]

where $\delta_{g}(\boldsymbol{r})$ and $\boldsymbol{v}(\boldsymbol{r})$ are the galaxy overdensity and the velocity field, respectively, at some location $\boldsymbol{r}$. In a real survey, we are only able to measure the radial peculiar velocities of galaxies $u(\boldsymbol{r})$, which means that we are also only able to construct the radial momentum field $\rho_{||}(\boldsymbol{r})$. Additionally, both the measured density and velocity fields contain non-linear contributions and are potentially biased tracers of the true underlying density and velocity fields (although there is good reason to believe the latter is unbiased on large scales, Desjacques & Sheth 2010; Jennings et al. 2015; Chen et al. 2018). Nonetheless, the measured radial momentum field can still be used as a tracer of cosmology.

As with any field, we can consider measuring the radial momentum field on a grid and taking the Fourier transform. Given the nature of the radial momentum field, we can simply do this by averaging the peculiar velocities of galaxies in cells using some interpolation scheme, which can be corrected for after the Fourier transform has been taken. Note that, unlike when trying to estimate the velocity field directly, the momentum field does not suffer from ‘zero-valued’ regions, i.e., parts of the gridded space where the velocity field is artificially set to zero due to the absence of tracers. In the absence of any galaxies due to a survey selection function, the momentum field is both measured and expected to be zero (because the overdensity is expected to be zero given our knowledge of the survey selection function), and so we are able to obtain results without complex corrections or gridding schemes (such as those used in Jennings et al. 2011; Jennings 2012 or Koda et al. 2014). However, we do have to modify our expectations for what we have measured; we have not measured the velocity field on the grid, and so we are not able to measure the velocity power spectrum directly.

We can demonstrate what we have measured by correlating the radial momentum field in Fourier space. We can write the radial momentum power spectrum as

[TABLE]

which highlights the fact that the power spectrum of the radial momentum field is related to the velocity power spectrum (it is the first term in the above equation), but is not fully equivalent. On linear scales, the velocity power spectrum dominates due to the $1/k^{2}$ scaling of the velocity field relative to the density field. However, on non-linear scales higher order terms arise from the convolution of the density and velocity fields at the same location. Although at first glance these appear to be three- and four-point functions, in practice they still only contain two distinct locations, and so can be computed relatively simply using perturbation theory (Vlah et al. 2012, 2013; Okumura et al. 2014, see also Appendix A). The equivalence of the velocity and momentum power spectra on linear scales makes this measurement an excellent tracer of gravity, whilst the density field dependence on small scales can be used to partially break the degeneracy between galaxy bias and the growth rate of structure seen in measurements of the galaxy density power spectrum.

Applying this to a real survey introduces some complexity. Firstly, in order to assign the velocities to a grid, we need an estimate of their position. This can be done one of three ways: using the measured redshifts and assuming a cosmological model to put the galaxies at their redshift-space distance; using the redshifts and peculiar velocities to estimate the true distance to each galaxy; or using their positions inferred from reconstructions of the density field using linear theory. In this work, we advocate the former. This is in contrast to the work of Park (2000) and Park & Park (2006), but as we will show, the use of redshift-space positions (which effectively replaces $\delta_{g}$ is Eq. 2 with the redshift space overdensity $\delta^{s}_{g}$ and can also be modelled with perturbation theory) only changes the measured momentum power spectrum on small, noise-dominated scales. Using instead the true distances (from the peculiar velocities or via reconstruction) could introduce considerable measurement or systematic error which may be harder to theoretically model.

Secondly, to apply this to a real survey, we need to formally derive an estimator that accounts for the discrete sampling of the fields and the survey geometry. This will be the focus of the next Section. It is also worth noting that extracting constraints from these measurements does require accurate knowledge of the survey geometry and selection function. However, ensuring a well-defined selection function for future peculiar velocity surveys will be extremely advantageous for a number of reasons beyond the requirements of the method developed herein.

3 Estimating the Redshift Space Power Spectrum

Galaxies typically fall towards dense regions of the Universe and have 3-dimensional velocities, however the redshift-space and real-space positions only differ due to the line-of-sight velocity. The net result of this on large scales is that an otherwise spherical distribution of galaxies will appear squashed or extended along the line-of-sight; the transverse positions are unaffected, whilst the line-of-sight positions are moved closer to the overdensities influencing the large scale motion of the galaxy. In terms of the clustering of galaxies, whereas in an isotropic and homogeneous universe we might expect the correlation function or power spectrum to depend only on the separation between galaxies, the use of redshift-space positions introduces a measurable difference along and transverse to the line-of-sight. The strength of this anisotropy is proportional to the typical velocity of objects along of the line-of-sight and so can be used to determine the growth rate of structure.

Although the spherical symmetry of two-point clustering measurements is broken in redshift-space, they do retain symmetry about the line-of-sight. Hence, a common technique is to decompose the redshift-space correlation function $\xi(\boldsymbol{r})$ or power spectrum $P(\boldsymbol{k})$ into multipoles using the Legendre polynomials $L(\mu)$ e.g.,

[TABLE]

where $r$ is the separation and $\mu$ is the cosine of the angle to the line-of-sight. Fig. 1 shows a diagrammatic representation of this. As the angle $\mu$ arises from the bisection of the vector between pairs of objects and the line-of-sight it will vary across a survey. For sufficiently distant galaxies or small survey areas all possible values of $\mu$ are approximately equal and we only need to consider a single fixed line-of-sight. This is commonly referred to as the “distant-observer” or “global plane-parallel” approximation. Given the wide-area and low redshift of most peculiar velocity surveys we do not wish to make such an approximation here a priori, although it will be used in some cases later.

Allowing for a varying line-of-sight, measuring the anisotropic correlation function can be achieved by counting the number of pairs of galaxies as both a function of their separation and $\mu$. The power spectrum on the other hand cannot be calculated so easily because the varying line-of-sight does not align with a Cartesian grid, which in turn makes Fourier transforms difficult. One alternative is to decompose the redshift-space fields into spherical harmonics (Heavens & Taylor, 1995), but a more popular method is to define a different $\mu$ for different subsamples of objects (Cole et al., 1994). This technique was extended to its logical conclusion by defining a separate line-of-sight for each pair of galaxies (Yamamoto et al., 2006) such that each pair is considered ‘locally parallel’, but the line-of-sight can still vary across the survey. As such this assumption is typically called the “local plane-parallel” approximation (Samushia et al., 2012; Beutler et al., 2014; Yoo & Seljak, 2015). It was recently demonstrated (Bianchi et al., 2015; Scoccimarro, 2015) that the monopole, quadrupole and hexadecapole ($l=0$,$2$,$4$ moments) of the galaxy power spectrum can be computed with only 1+6+15 Fast Fourier transformations respectively under this approximation (although there is a subtly different convention used to calculate the line-of-sight, see Beutler et al. 2019) and this technique has been quickly adopted into the standard practice for estimating the redshift-space clustering of galaxies from large galaxy redshift surveys (e.g., Blake et al. 2011; Oka et al. 2014; Beutler et al. 2014; Gil-Marín et al. 2016; Beutler et al. 2017; Gil-Marín et al. 2018). In this work, we will demonstrate that a similar technique can be developed for measuring the redshift-space clustering of the momentum field.

It is worth noting here that neither the Yamamoto et al. (2006) estimator, nor it’s implementation using Fast Fourier Transformations (Bianchi et al., 2015; Scoccimarro, 2015), are immune from ‘wide-angle effects’ on large scales, which if left unaccounted may bias constraints obtained from anistropic power spectrum measurements. There is a long history of works attempting to account for these effects (see i.e., Hamilton & Culhane (1996); Zaroubi & Hoffman (1996) for some early examples, or Castorina & White (2018); Beutler et al. (2019) and references therein for more recent studies). We do not consider such effects in this work and did not find any associated bias in the results when the methods developed in this paper were applied to current data in Paper II. However, given the low redshift, large sky area and increased precision of upcoming and planned peculiar velocity surveys similar techniques should be explored to ensure we can obtain unbiased results in the future.

3.1 Momentum Power Spectrum

In this section, we present a new estimator for the multipoles of the redshift-space line-of-sight momentum power spectrum (which from now on will simply be called the momentum power spectrum), as measured from a PV survey with arbitrary number density and weighting. This formalism effectively combines the methods of Park & Park (2006), Yamamoto et al. (2006) and Blake et al. (2018) resulting in an estimator that looks and behaves similar to the multipoles of the galaxy density power spectrum.

We begin by defining the weighted radial momentum field

[TABLE]

where $n_{g}(\boldsymbol{r})$ and $u(\boldsymbol{r})$ are the number of galaxies and line-of-sight peculiar velocity at location $\boldsymbol{r}$ respectively. $w(\boldsymbol{r})$ and $A$ are some convenient weights and normalisation that we will define later111For comparison, the weighted density field used to define the galaxy power spectrum estimator in Feldman et al. (1994); Yamamoto et al. (2006) is

$F^{\delta}(\boldsymbol{r})=w(\boldsymbol{r})[n_{g}(\boldsymbol{r})-\alpha n_{s}(\boldsymbol{r})]/A,$

(5)

where a synthetic random catalogue containing $1/\alpha$ times as many objects as there are galaxies is used to calculate $n_{s}$. The weights for these two fields may not be the same..

If we then take the expectation value of the correlation of the weighted momentum field at different locations we obtain

[TABLE]

Using the classic trick of dividing the survey volume V into infinitesimally small cells with occupancy [math] or $1$, Park & Park (2006) showed that we can write

[TABLE]

such that

[TABLE]

and where $\bar{n}(\boldsymbol{r})$ is the expected number of galaxies at $\boldsymbol{r}$, $\xi^{p}(\boldsymbol{r}-\boldsymbol{r}^{\prime})$ is the radial momentum correlation function, and $\langle v^{2}(\boldsymbol{r})\rangle$ is the value of this at zero-lag. Switching now to Fourier space, we define the multipole moments of the momentum power spectrum in the usual manner

[TABLE]

and, under the ‘local plane-parallel’ approximation (such that $L_{\ell}(\hat{\boldsymbol{k}}\cdot\hat{\boldsymbol{r}})=L_{\ell}(\hat{\boldsymbol{k}}\cdot\hat{\boldsymbol{r}}^{\prime})$), write the momentum version of the Yamamoto et al. (2006) estimator

[TABLE]

where $L_{\ell}(\hat{\boldsymbol{k}}\cdot\hat{\boldsymbol{r}}^{\prime})$ are the Legendre polynomials, and we average over shells in $\boldsymbol{k}$-space as represented by the integral over $d\Omega_{k}$. The volume factor, $1/V$ arises due to our Fourier transformation convention. Taking the expectation value of $|F^{p}_{\ell}(\boldsymbol{k})|^{2}$, and making use of Eq. 8

[TABLE]

where we have substituted the momentum correlation function for its Fourier counterpart, the momentum power spectrum $P^{p}(\boldsymbol{k},\boldsymbol{r})$, and compressed the $\boldsymbol{k}$-dependent shot-noise term into

[TABLE]

From here, we follow the method of Blake et al. (2018) (Eqs. 13 and 14 therein) and substitute the multipole expansion in Eq. 9 into the above expression

[TABLE]

where

[TABLE]

Finally, in the absence of any window function, the weights and number density must be constant as a function of $\boldsymbol{r}$. In this case $G(\boldsymbol{k}-\boldsymbol{k}^{\prime})=\bar{n}wV\delta^{D}(\boldsymbol{k}-\boldsymbol{k}^{\prime})$ and the orthogonality condition of the Legendre polynomials means $S^{*}_{\ell,\ell^{\prime}}(\boldsymbol{k},\boldsymbol{k}^{\prime})=\bar{n}w\delta^{D}(\ell-\ell^{\prime})\delta^{D}(\boldsymbol{k}-\boldsymbol{k}^{\prime})/(2\ell+1)$ where $\delta^{D}$ is the Dirac delta. Hence, if we set the normalisation constant to the same value used when estimating the galaxy density power spectrum $A^{2}=\int d^{3}r\,w^{2}(\boldsymbol{r})\bar{n}^{2}(\boldsymbol{r})$, then in the absence of any window function $\langle|F^{p}_{\ell}(k)|^{2}\rangle=P^{p}_{\ell}(k)+N^{p}_{\ell}(k)$ and we can formally express our estimator for the multipoles of the redshift-space momentum power spectrum in the familiar form222Again, for comparison, the same derivation for the galaxy density field in Yamamoto et al. (2006) results in

$\widehat{P^{\delta}_{\ell}}(k)=|F^{\delta}_{\ell}(k)|^{2}-N^{\delta}_{\ell}(k),$

(17)

where $N^{\delta}_{\ell}(k)=\int d^{3}r\,w^{2}(\boldsymbol{r})\bar{n}(\boldsymbol{r})L_{\ell}(\hat{\boldsymbol{k}}\cdot\hat{\boldsymbol{r}})$.

[TABLE]

This completes the formal definition of our estimator for the multipoles of the momentum power spectrum. In practice, when estimating this from a dataset we use the Fourier-based method of Bianchi et al. (2015) (see also Scoccimarro 2015), replacing the weighted density field $F^{\delta}(\boldsymbol{r})$ used in their work with $F^{p}(\boldsymbol{r})$ from Eq. 5. For the shot-noise terms, we approximate the momentum correlation function at zero-lag as a sum over $N$ measurements

[TABLE]

Looking at the shot-noise term in Eq. 13, the above derivation demonstrates that the measured momentum power spectrum depends on both the use of discrete tracers of the velocity field, and the variance in the velocity field itself. This can be seen more clearly if we take the case for $\ell=0$ without any window function $N^{p}_{0}=\langle v^{2}\rangle/\bar{n}$, which can be compared to the well-known limit for the shot-noise on the density power spectrum in this regime $N^{\delta}_{0}=1/\bar{n}$. The variance in the velocity field contains contributions from the intrinsic variance and the measurement errors, and so remains even for a field where the peculiar velocities are perfectly known. In the absence of any errors on the PV measurements the expectation of the velocity field is independent of position, and so we can use all the measurements to evaluate the shot noise term. In practice, the presence of errors on the PV measurements adds a radial dependence to $\langle v^{2}(\boldsymbol{r})\rangle$. A suitable value for this can be obtained by fitting the errors on the data as a function of redshift or distance (i.e., as in Howlett et al. 2017b), or by averaging over measurements in radial bins. We will further examine the contribution from measurement errors, and the fact that the mean velocity itself may not be zero, in Section 6.

The obvious similarities between the estimator for the multipoles of the density and momentum power spectra highlight its potential as a tool for extracting cosmological information from current and future PV surveys:

• •

Codes to estimate the galaxy density power spectrum can be trivially modified to also measure the momentum power spectrum, allowing for efficient measurement of both of these from a survey.

• •

In a real survey, the measurements of the momentum power spectrum are affected by the survey window function through the convolution with the $G$ and $S_{\ell,\ell}$ terms in Eqs. 14-16. However, direct comparison between these equations and the convolution between the galaxy density power spectrum and the window function in Feldman et al. (1994), Yamamoto (2003) or Blake et al. (2018), shows that the convolution of the the momentum power spectrum with the survey window function is the same. Hence, it can be evaluated using the same random catalogue that mimics the survey geometry.

• •

The covariance matrix of the multipoles of the momentum power spectrum estimated from a PV survey can be estimated both theoretically (under some approximations, see Section 4) or using an ensemble of realistic simulations of the survey, again in an identical manner to that used for measurements of the galaxy density power spectrum.

• •

Given all of the above, and a model for the redshift-space momentum power spectrum (e.g., Appendix A) extracting cosmological constraints from the momentum power spectrum can be performed using the same methods as used for the galaxy density power spectrum, e.g., using a Monte-Carlo Markov Chain approach to compute the posterior probabilities.

We demonstrate our estimator of the momentum power spectrum using simulated halos in two different mass bins at redshift $z=0$. The simulation consists of $2560^{3}$ particles in a box of length $1800\,h^{-1}\,{\rm Mpc}$ and halo catalogues were generated using the implementation of the 3D Friend-of-Friends algorithm in the VELOCIraptor code (Elahi et al., 2011). As a first step towards creating conditions representative of real observations, we calculate the momentum power spectrum by placing an observer in the centre of the simulation box and computing line-of-sight velocities for each halo. These are then assigned to a $512^{3}$ grid using nearest neighbour interpolation. The corresponding Nyquist frequency is $k_{Ny}=0.89\,h\,{\rm Mpc}^{-1}$, much larger than the scales we are interested in. We compare the momentum power spectra in both real- and redshift-space, where for the latter we use the line-of-sight velocity to perturb the positions of the halos.

Our measurements are shown in Fig. 2 and highlight some interesting features of the momentum power spectrum:

• •

Unlike the galaxy density power spectrum, the real-space momentum power spectrum measured using line-of-sight velocities is not isotropic. The use of peculiar velocities adds a direction dependence to the momentum power spectrum even in real-space.

• •

Using the redshift-space positions of the halos does not change the momentum power spectrum on linear scales. The lower panels of Fig. 2 show the ratio of the redshift- and real-space spectra. The only effect of using the redshifts to infer the positions of objects when computing the momentum power spectrum is a damping on non-linear scales. The equivalence of the real and redshift-space momentum power spectrum on linear scales is demonstrated theoretically in Vlah et al. (2013) and Okumura et al. (2014) and arises due to the fact that the velocity power spectrum contribution to the momentum power spectrum does not change if we use redshift-space positions.

• •

In the linear regime, all odd multipoles and multipoles beyond the quadrupole are expected to vanish. This can be demonstrated using the fact that the large scale radial momentum power spectrum is equivalent to the radial velocity power spectrum (Park, 2000), which in turn means

[TABLE]

where $\mu$ is cosine of the angle to the line-of-sight, $f$ is the growth rate of structure, $a$ is the scale factor, $H$ is the Hubble parameter and $P^{p}_{L}(k)$, $P_{L}(k)$ are the linear momentum and matter power spectra. Multiplying by the relevant Legendre polynomials and integrating over the line-of-sight gives

[TABLE]

We plot the linear theory predictions for the monopole and quadrupole of the momentum power spectrum in Fig. 2. The momentum power spectrum is a pure tracer of the growth rate of structure on large scales, independent of galaxy bias. This property makes it competitive as a probe of the growth rate, even given the typically small numbers of galaxies with direct PV measurements and the intrinsic error in each velocity measurement.

• •

On non-linear scales the momentum power spectrum becomes dominated by the power spectrum of the convolution between density and velocity. That is, the power spectrum arising from correlations between $\delta(\boldsymbol{r})v(\boldsymbol{r})$ and $\delta(\boldsymbol{r}^{\prime})v(\boldsymbol{r}^{\prime})$. This causes an increase in the momentum power spectrum above that expected from the velocity power spectrum and proportional to the halo/galaxy bias. This can be seen in Fig. 2 as the excess in power for $k\gtrsim 0.1\,h\,{\rm Mpc}^{-1}$, which is more evident for larger mass halos. As shown in Section 8 and Appendix A this can be well modelled using perturbation theory.

The final feature of note in Fig. 2 is the departure from linear theory seen on large scales. This arises due to the window function imposed on the simulation when placing an observer. As with the galaxy density power spectrum, this window function convolves the observed power, increasing the correlation between modes and reducing the large scale power. Given the relatively shallow depth of PV surveys, we expect this to cause a significant reduction in the measure power spectra from real data and will need to be well modelled. Fortunately, the effect of a window function on the momentum power spectrum has already been derived in Eqs. 14-16. In Section 8 we include this convolution using the same (matrix-multiplication) method as is commonly used for the galaxy density power spectrum.

4 Covariance matrix of the momentum power spectrum

Given the above estimator for the momentum power spectrum we can compute the covariance matrix of the momentum power spectrum as

[TABLE]

Substituting Eqs. 5 and 10 into this expression, we find

[TABLE]

Although this expression is daunting, we can make some progress by rewriting $n_{g}(\boldsymbol{r}_{i})u(\boldsymbol{r}_{i})=\bar{n}(\boldsymbol{r}_{i})(1+\delta_{g}(\boldsymbol{r}_{i}))u(\boldsymbol{r}_{i})=\bar{n}(\boldsymbol{r}_{i})\rho(\boldsymbol{r}_{i})$ (where $\delta_{g}(\boldsymbol{r})$ is the galaxy overdensity) and assuming that the radial momentum field $\rho(\boldsymbol{r}_{i})$ is Gaussian. Although this assumption is certainly not true (even if the density and velocity fields are both Gaussian, their product would not be), on large scales the radial momentum field is simply proportional to the radial velocity. The large scale radial velocity is approximately Gaussian in nature333Note, this is not to be confused with the pairwise velocity between objects, which is non-Gaussian (Scoccimarro, 2004), or the distribution of velocity amplitudes, which has a Maxwell-Boltzmann distribution (Coles & Lucchin, 1995) and so this gives credence to the assumption that the radial momentum field is too.

In this case, we can then use Wick’s Theorem for a Gaussian random field to write

[TABLE]

and substituting this into Eq. 23,

[TABLE]

where

[TABLE]

and

[TABLE]

Finally, we substitute in Eq. 7 and follow the steps in Eqs. (33) and (34) of Blake et al. (2018), neglecting the convolution with $P(\boldsymbol{k}^{\prime\prime},\boldsymbol{r}^{\prime})$ which allows us to take this outside the integral over $\boldsymbol{k}^{\prime\prime}$, and defining an effective momentum power spectrum

[TABLE]

Following these steps,

[TABLE]

From here, the covariance matrix of the momentum power spectrum in the Gaussian limit has the same form as the galaxy density power spectrum in Blake et al. (2018). Hence, this can be computed for a given survey using the same decomposition into spherical harmonics, and a model for the momentum power spectrum and associated shot noise. Given that the theory and procedure is the same as that in Blake et al. (2018) it will not be repeated here. It is interesting however, to consider the case without a survey window function and constant variance in the velocity field (such that $\langle v^{2}(\boldsymbol{r})\rangle$ is also independent of radius). In this limit, the weights, galaxy number density and velocity field variance must be independent of position and can be moved outside of the integral in Eqs 28 and 30. Following some simplification, reverting to the distant-observer approximation (as opposed to the local plane-parallel approximation used to derive the previous expressions), and rewriting the equation in terms of $\mu$ instead of $\hat{\boldsymbol{k}}\cdot\hat{\boldsymbol{r}}$, we find that the covariance matrix of the momentum power spectrum can be expressed

[TABLE]

i.e., in the same form as the galaxy density power spectrum covariance matrix (Yamamoto et al., 2006; Taruya et al., 2010; Grieb et al., 2016), where $V_{k}=4\pi k^{2}\Delta k$ and $\Delta k$ is the bin width.

We can demonstrate the validity of our analytic formulae by comparing to cubic simulations. We generated $256$ approximate dark matter simulations using the l-picola simulation code (Howlett et al., 2015a) and the same simulation parameters as in Section 3. Halos are also identified in these simulations using the same procedure. From these simulations we computed the redshift space density and momentum power spectrum for halos in the mass range $12.4h^{-1}\mathrm{M_{\odot}}<M<12.6h^{-1}\mathrm{M_{\odot}}$ using the z-axis of the simulation as the line-of-sight. The full covariance matrix for $P^{p}_{0}(k)$ and $P^{p}_{2}(k)$ as well as the galaxy density power spectrum multipoles was then computed from all 256 measurements. The use of a fixed line-of-sight for the whole simulation differentiates this from the results shown in Figure 2, and was chosen to remove the effect of the window function, which is expected to change the diagonal elements of the covariance matrix of both the density and momentum power spectrum away from the Gaussian prediction (see Howlett & Percival 2017 for a detailed discussion of how a window function affects the galaxy density power spectrum covariance matrix).

We calculated the Gaussian prediction using a model for the anisotropic power spectrum (both density and momentum) based on the averaged measurements of the multipoles from the simulations, $P(k,\mu)=P_{0}(k)+L_{2}(\mu)P_{2}(k)$. This was then used as input to Eq. 31, and the similar equation for the galaxy density power spectrum multipoles. Although the galaxy density power spectrum should contain higher order multipoles, even in the linear regime, these measurements are typically noise dominated and are not expected to contribute significantly to the Gaussian covariance matrix. Fig. 3 shows the measured diagonal elements of the covariance matrices of the density and momentum power spectra, as well as the cross-covariance between the monopole and quadrupole alongside these Gaussian predictions.

In all cases we see good agreement between the measurements and theory on the largest scales where the density and momentum fields are approximately Gaussian. However, the variance is underestimated as we go to smaller scales, where non-Gaussian terms enter the covariance matrices and the disagreement is much more severe for the momentum power spectrum. This is due to the contribution to the momentum power spectrum from the correlation between the density and velocity fields, which even in the case of Gaussian random fields would produce a non-Gaussian signature.

Using the same simulations, in Fig. 4, we highlight portions of the correlation matrix $R_{ij}=C_{ij}/\sqrt{C_{ii}C_{jj}}$ for the density and momentum multipoles. We can see that the the cubic simulations show negligible off-diagonal covariance for the galaxy density power spectrum, as to be expected for halos in the absence of a window function and given the agreement between the diagonal elements and the Gaussian prediction. However, the momentum power spectrum shows significant correlation between modes on quasi-linear scales driven by the correlated nature of the density and velocity fields. This is particularly striking for the monopole, but reduced for the quadrupole moments. The cross-covariance between the two multipoles is also much more significant away from $k_{i}=k_{j}$ for the momentum power spectra compared to the galaxy density power spectrum.

In Fig. 5 we plot slices through the joint correlation matrix of the density and momentum power spectra alongside the Gaussian prediction of Eq. 31. From this we see that the high degree of cross-covariance between bins for the momentum power spectrum is a result of non-Gaussian contributions to the covariance matrix of the momentum power spectra. However, we also find that the density and momentum power spectra are largely uncorrelated except when looking at measurements in the same $k$-bin. Given that the momentum field contains information from the density field, this is perhaps unexpected. This will be explored further in the next section, where we also demonstrate how we can compute a Gaussian prediction cross-covariance matrix of the density and momentum power spectra.

4.1 Cross-covariance between density and momentum power spectra

In order to combine measurements of the density and momentum power spectrum we need to account for potential cross-covariance between these two sets of measurements. Given a set of mock catalogues, these can be computed by brute force. However, we can also compute these in the Gaussian limit using the same method as for the covariance of the momentum power spectrum.

Starting from the weighted density (Feldman et al., 1994) and momentum fields, we express the cross-covariance between density and momentum power spectra as

[TABLE]

The large number of terms compared to the previous derivation arises due to the use of a random field in the calculation of the weighted density field which is not used for the weighted velocity field. Although this seems tedious to evaluate, we can again make headway using Wick’s theorem, .e.g.,

[TABLE]

and similarly for the other terms. Substituting these into Eq. 33 results in a large cancellation of terms and a number of terms involving the correlation between the momentum field and random field at two different locations (which vanishes), such that all that remains is

[TABLE]

Finally, we introduce the cross-correlation function between the density and momentum fields $\xi^{p\delta}(\boldsymbol{r}-\boldsymbol{r}^{\prime})$, which has Fourier counterpart $P^{p\delta}(\boldsymbol{k},\boldsymbol{r}^{\prime})$, using

[TABLE]

where we have neglected shot-noise in the cross-correlation function. Then, following the same steps as in Section 4, we finally arrive at

[TABLE]

where

[TABLE]

In the limit of no window function and following the same procedure used to derive Eq. 31, this reduces to

[TABLE]

i.e., the cross-covariance between the density and momentum power spectra is a function of the cross power spectrum between the two fields. This mimics the cross-covariance between two different measurements of the galaxy density power spectrum from different galaxy samples often considered in a multi-tracer analysis (i.e., Smith 2009; White et al. 2009).

As before, we can use simulations to test this theory. We use the same halos as Section 4 and measure the multipoles of the cross-power spectrum using the same procedure before computing the analytic cross covariance. The results are shown in Fig. 6, where we plot the diagonal elements of the cross-covariance of the density and momentum power spectra alongside the correlation matrix. Again, the theoretical result agrees well with the measurements on large-scales, and then underestimates the cross-covariance when non-Gaussian terms become important. This transition occurs on smaller scales than for the momentum power spectrum covariance matrix. In Gaussian theory, the cross-covariance between the monopole and quadrupole of the density and momentum fields respectively should be the same as the cross-covariance between the monopole and quadrupole of the momentum and density fields respectively, i.e., Eq. 39 is invariant under the exchange of $\ell$ and $\ell^{\prime}$. The measurements from the simulations agree with this prediction within noise. In Fig. 6, we also show the cross-correlation matrix between the multipoles of the density and momentum fields. The full set of measurements is only correlated on large scales when $k_{i}=k_{j}$; on small scales, the large intrinsic variance in the velocity field (on the order of $150^{2}-300^{2}\mathrm{km^{2}\,s^{-2}}$; Turnbull et al. 2012) decorrelates the density and momentum power spectra. An alternative way of viewing this is that the $\langle v^{2}(\boldsymbol{r})\rangle$ component of the shot-noise adds such a large component to the variance of the non-linear momentum power spectrum that the cross correlation coefficient, which is proportional to the inverse of this variance, becomes negligible.

As a final note, it is worth stating that we do not explore further the cross-power spectrum between the density and momentum fields in this work, using it only to test the covariance matrix of the auto-correlation of the density and momentum fields. Our forecasts for the mock catalogues used in Section 7 suggest this can improve the growth rate constraints by $\sim 10\%$ or more over the already considerable improvement obtained from including the momentum power spectrum multipoles in addition to the galaxy density power spectrum multipoles. As such this is an interesting measurement in its own right and an important step toward fully realising the information in the combined velocity and density fields measured from PV surveys. However, measuring this in the presence of observational effects is significantly more involved than in the cubic simulations used in this sections and requires in-depth study. As such we elect to study this in future work and focus the remainder of this paper on the auto-power spectrum of the momentum field. It is important to realise that we do include the cross-covariance between the density and momentum power spectra in Section 7; it is only the measurement of the cross-spectrum and the resulting cross-covariance between this measurement and the density and momentum power spectra that we defer to later work.

5 Optimal weights

In the above derivation of our estimator for the momentum power spectrum and its covariance matrix we have purposely included a free weight parameter $w(\boldsymbol{r})$. Following the standard set by Feldman et al. (1994) and Yamamoto (2003); Yamamoto et al. (2006), we can compute the value for this weight that minimises the variance in our estimation of the momentum power spectrum multipoles. To make an analytic solution for this weight tractable, we start by applying the distant observer approximation to Eq. 25 such that we can combine the $Q_{\ell,\ell}$ terms into a single integral over $\boldsymbol{r}$,

[TABLE]

The optimal weights are then found by requiring this variance to be stationary with respect to a small change in weights $w(\boldsymbol{r})\rightarrow w_{0}(\boldsymbol{r})+\delta w(\boldsymbol{r})$, .i.e, $\delta\langle\Delta P^{p}_{\ell}(k)^{2}\rangle/\delta w(\boldsymbol{r})=0$. This can be solved as in Feldman et al. (1994); Yamamoto (2003); Yamamoto et al. (2006) using functional derivatives (see also Smith & Marian 2015 for a much more in depth discussion of how this is accomplished). Alternatively, we can save considerable time by direct comparison with these previous results, finding

[TABLE]

as the weight factor that minimises the variance444We did verify that performing the full derivation using functional derivatives gives the same solution!. As with the FKP weights, this weighting scheme balances contributions from shot-noise, which in this case depends on both the number density of galaxies and the variance in the velocity field, and the power spectrum itself. In practice, we set the weights by fixing $\langle v^{2}(\boldsymbol{r})\rangle$ based on the variance in the measured PVs in the data and choosing a constant value $P^{p}(\boldsymbol{k},\boldsymbol{r})=P^{P}_{FKP}$, as is done for the usual FKP weights. We explore different choices for $P^{p}_{FKP}$ in Section 6, where we test our estimator on realistic mocks with varying number density.

The above weight factor minimises the variance of the estimator in the Gaussian limit. This is not necessarily the optimum choice of weights when including non-Gaussian contributions (which as shown in Section 4 are significant for the momentum power spectrum), nor when optimising cosmological parameter constraints. In particular, one of the main uses of the momentum power spectrum will be constraining the growth rate of structure in combination with galaxy density power spectrum measurements. In this case, optimal weights could be obtained numerically following the method of Pearson et al. (2016). In Section 7 we use such a method to identify the optimum choices for the constants $P_{FKP}$ (used for the density power spectrum measurement) and $P^{p}_{FKP}$ (used for the momentum power spectrum measurement) by maximising the Fisher information (or equivalently minimising the ideal error) on $f\sigma_{8}$.

6 Tests on realistic mocks

In the remainder of this work, we focus on applying the new estimator for the momentum power spectrum to a set of realistic simulated peculiar velocity surveys (mocks) that have a sky coverage and redshift distribution representative of what could be achieved with next generation surveys. For our test case we consider a sample of Type IA supernovae measured with LSST (Ivezic et al., 2008) matching the predictions of Howlett et al. (2017c); Kim et al. (2019). We also incorporate errors on the PVs for each mock galaxy and look at the effects of a non-zero Bulk Flow on the measured power spectra.

The mocks were generated using the same approximate l-picola halo catalogues as in Section 4. Galaxies were added into these halo catalogues using a hybrid Halo Occupation Distribution (HOD, e.g., Zheng et al. 2007) and Subhalo Abundance Matching (SHAM, e.g., Conroy et al. 2006) approach. Full details of the construction of the mock catalogues is detailed in a companion paper (Qin et. al., in preparation, herein Paper II), where we use the same technique to produce mock galaxy catalogues as part of the application of this work to real data.

In brief, whilst the approximate nature of the l-picola simulations allows us to reproduce halos reasonably well (see Howlett et al. (2015b) for a comparison with fully non-linear N-Body simulations) it makes them unsuitable for identification of subhalos. Instead, motivated by measurements of the subhalo mass function (Springel et al., 2008; Giocoli et al., 2008; Han et al., 2016; Elahi et al., 2018), we add subhalos to halos by drawing from a power law. Subhalos are placed within the halos assuming an NFW profile; drawing from the inverse CDF for the positions (Robotham & Howlett, 2018) and using the virial theorem to then determine the satellite velocity as a function of position (Lagos et al., 2018). This prescription contains two free parameters, which are tuned to the galaxy density power spectrum in the data. For the mock catalogues in this paper, we use the parameters listed under the “2MTF” sample in Paper II, which is chosen to roughly account for the idea that Type IA supernovae are preferentially found in star-forming galaxies (Smith et al., 2012; Andersen & Hjorth, 2018), which are also those used in the construction of the 2MTF sample.

Given a catalogue now containing central and satellite galaxies, we assign luminosities to these following Howlett et al. (2017b) and Qin et al. (2018), drawing from the K-band luminosity function of Kochanek et al. (2001). These are assigned to the halos and subhalos by rank ordering against their mass. We artificially introduce scatter into this matching to mimic the observed scatter between halo mass and luminosity by drawing the ‘matching’ log luminosity from a Gaussian distribution where the standard deviation is an additional free parameter, also tuned to the data (i.e., in the same spirit as Avila et al. 2018). Note that this ‘matching’ log luminosity is only used to assign the luminosities to the galaxies, and not for each galaxy’s actual luminosity.

Next, we place 8 maximally separated observers within each simulation, compute apparent magnitudes for the galaxies and apply a magnitude cut of $K\leq 14.75$, a redshift cut of $z\leq 0.2$ and keep only objects with declination $-80.0<\delta<15.0$, which roughly matches what could be observed with LSST. This returns a total number of objects $(\sim 150,000)$ and a redshift distribution similar to that proposed in Kim et al. (2019). Finally, errors are assigned to the velocities of each galaxy using $0.08H_{0}r(z)$ where $r(z)$ is the comoving distance to the galaxy (Howlett et al., 2017a, c).

Optimal measurements of the momentum power spectra require an estimate of $\langle v^{2}(\boldsymbol{r})\rangle$. As mentioned in Section 3, for real PV measurements this will be dominated by the measurement error on the individual galaxies, which typically increases as a function of distance from the observer. As such, we can obtain an estimate for this by fitting the errors in the data as a function of redshift or distance. For the mocks used herein we have assigned errors using the distance itself and so we simply adopt $\langle v^{2}(\boldsymbol{r})\rangle=(0.08H_{0}r(z))^{2}+300^{2}\mathrm{km^{2}\,s^{-2}}$, where the additional $300^{2}\mathrm{km^{2}\,s^{-2}}$ models the contribution from the intrinsic variance in the velocity field. This value is used for calculating the optimal weights and the window function of the data. For computing the shot-noise we use a sum over all galaxies (and their velocities) as is done for the shot-noise in the galaxy density power spectrum.

6.1 Effects of PV measurement errors

We begin by looking at the effect of PV measurement errors on the measurements of the momentum power spectrum. For our mocks we find we are able to produce robust measurements of both the monopole and quadrupole, although this may not be true for current or future surveys depending on the cosmological volume, shot noise or velocity errors. We do not explore higher order multipoles as these are expected to be zero from linear theory and were found to be completely noise-dominated in our mocks. Figures 7-10 show a series of measurements without the inclusion of measurement errors and also when including errors and four different weighting schemes; unweighted, and $P^{\,p}_{FKP}=[2\times 10^{8},1\times 10^{9},5\times 10^{9}]\,{\rm km^{2}\,s^{-2}}\,h^{-3}\,{\rm Mpc^{3}}$.

Turning first to Fig. 7 we see that different choices for $P^{\,p}_{FKP}$ change the effective depth of the measurement, in the same way as the standard FKP weights. However, the FKP weights typically (relatively) upweight more distant regions of the survey where the number density of objects is lower, increasing the effective cosmological volume and decreasing the variance in the power spectrum. Instead, the optimal weights for the momentum power spectrum have to also account for the fact that more distant objects have larger errors, and hence larger $\langle v^{2}(\boldsymbol{r})\rangle$. This means that decreasing the variance is not as simple as upweighting more distant objects.

Next, looking at the monopole and quadrupole power spectra for the different cases (Fig 8), we see that the inclusion of PV errors does not significantly change the mean power spectrum measured from the mocks. For all choices of weights the measurements with errors are almost indistinguishable from those without for the monopole. There are some differences for the quadrupole depending on the choice of $P^{p}_{FKP}$, however on large scales we attribute this to the fact that the distribution of quadrupole measurements can become non-Gaussian especially for the case of $P^{p}_{FKP}=2\times 10^{8}\,{\rm km^{2}\,s^{-2}}\,h^{-3}\,{\rm Mpc^{3}}$ given the emphasis this places on lower redshift measurements (a point which will be explored later and is a discussed at length in Paper II when measurements are made from real data). Overall, the measurements match what we would expect; little change in the power spectrum in the case of Gaussian distributed peculiar velocity errors, which, whilst included in our mocks by design, are also obtainable for modern PV surveys555Modern PV surveys such as 2MTF (Hong et al., 2014) or 6dFGSv (Campbell et al., 2014; Springob et al., 2014) typically contain measurements of the log-distance ratio, which is (for the most part) Gaussian distributed. Whilst the velocities are then generally log-normally distributed, we can obtain Gaussian velocities if we make use of the Watkins & Feldman (2015) estimator.. A subtle point of note is the change in the amplitude of the non-linear momentum power spectrum for different weightings, which arises due to different weights changing the effective galaxy bias of the sample, which in turn contributes to the portion of the momentum power spectrum arising from the convolution between the density and velocity fields.

Whilst we believe that the presence of Gaussian distributed PV errors should not bias the power spectrum measurements, what they do do is significantly increase the variance of the power spectrum, acting as additional shot-noise. The introduction of PV errors also adds random non-linear noise to the momentum power spectrum, significantly reducing the correlation between the monopole and quadrupole and between different scales. Although not highlighted here, the cross-correlation between the momentum and galaxy density power spectra is also reduced. These effects can be seen in Figs 9 and 10, where regardless of the weighting scheme the errors on the monopole and quadrupole are increased by up to an order of magnitude due to the PV errors whilst the cross-correlation coefficient between the monopole and quadrupole at fixed scale is reduced. The reduction in the correlation between different scales in the monopole and quadrupole can be best seen by comparing Fig. 10 with the right-hand panel of Fig, 4. It’s worth noting that whilst the correlation is reduced between Figs 10 and 4, the overall amplitude of the off-diagonal covariance is not; the window function actually increases the off-diagonal covariance in the mocks compared to the cubic simulations used in previous sections, but the inclusion of uncorrelated PV errors then adds so substantially to the diagonal covariance that the net effect is a large reduction in correlation.

Finally, although the effect of PV errors matches our expectations, different choices for $P^{\,p}_{FKP}$ in the presence of these errors change the relative variance of linear and non-linear scales and the cross-correlations in a complex way. The overall take-away from this is that whilst we have derived the optimal weights for the momentum power spectrum, the fact that it is common to choose a fixed $P^{p}(k)=P^{\,p}_{FKP}$ in evaluating these weights means some choices for $P^{\,p}_{FKP}$ will be better than others depending on how the momentum power spectrum is to be used. In particular, for measuring the growth rate of structure it may be more beneficial to minimize the errors on linear scales, whilst sacrificing the constraining power on small scales where galaxy bias becomes a nuisance parameter. We explore this further in Section 7 and Appendix B where we identify the choice of $P^{\,p}_{FKP}$ that maximises the information on $f\sigma_{8}$ obtainable from our mocks.

6.2 Effects of non-zero Bulk Flow

Another feature of real PV surveys is the presence of a non-zero Bulk Flow. Physically, this arises due to our own motion within the cosmic web, and is equivalent in our notation to a non-zero mean velocity, $\langle v(\boldsymbol{r})\rangle\neq 0$. A measured Bulk Flow (even if the true Bulk Flow is zero) can also arise due to non-uniform survey geometry and sampling of the velocity field (Andersen et al., 2016). It is interesting to ask what effect this has on the measured momentum power spectrum. To test this, we measure the momentum power spectrum monopole in our mock catalogues for a single weighting scheme ($P^{\,p}_{FKP}=2.4\times 10^{9}\,{\rm km^{2}\,s^{-2}}\,h^{-3}\,{\rm Mpc^{3}}$) after removing the observers Bulk Flow from each mock galaxy. Although modern, sophisticated methods for this exist (Nusser & Davis, 2011; Peery et al., 2018; Qin et al., 2018), we do this simply by first estimating the 3-dimensional Bulk Flow $\boldsymbol{B}$ using the standard Maximum Likelihood estimator Kaiser (1987) (which is valid given how our measured PVs are largely Gaussian) and then subtracting this component from the radial velocity of each galaxy using $u(\boldsymbol{r})\rightarrow u(\boldsymbol{r})-\boldsymbol{B}\cdot\hat{\boldsymbol{r}}$.

The results for this test are shown in Figure 11. We find that removing the Bulk Flow reduces the large scale monopole of the momentum power spectrum and slightly increases the quadrupole (although the latter remains unchanged within the errors). The variance and covariance of both of these is reduced on large ($k<0.05\,h\,{\rm Mpc}^{-1}$) scales. In subtracting the Bulk Flow, we are artificially setting $\langle v(\boldsymbol{r})\rangle=0$, which forces the spherically averaged power spectrum to go to zero as $k=0$. This is related to the well-known integral constraint in the galaxy density power spectrum arising from the assumption that the average overdensity over the observed survey volume is zero (which we are also assuming here too). The variance of the power spectrum is reduced as we are removing a portion of the cosmic variance associated the observer’s own motion, which even in $\Lambda$CDM can have a significant range of values. Although not shown, we also found that the presence of a Bulk Flow has negligible effect on the correlation between different scales in the density/momentum power spectrum multipoles.

When fitting the mocks in Section 7 we tested the fits with and without removing the Bulk Flow, finding that the errors on the growth rate were marginally improved after its removal, consistent with the reduction in the variance of the momentum power spectrum. The best-fit value of the growth rate did not change. During the fitting we included the standard correction for the integral constraint (i.e., following Ross et al. 2013; Beutler et al. 2014 for the galaxy density power spectrum), but do not explicitly include further correction for artificially setting the mean velocity to zero. Our mocks advocate that there is some small benefit in removing the bulk flow. In their original work, Park & Park (2006) removed the bulk flow and our mocks seem to validate their approach. However, our results could be strongly influenced by the fact that the mocks used here extend over quite a large cosmological volume and we only fit for $k\geq 0.02\,h\,{\rm Mpc}^{-1}$; our mock Bulk Flows are likely smaller than in current datasets with lower redshift limits or that cover smaller sky areas and their effects do not extend much past our lower $k$-limit. Without further theoretical work towards understanding the effect that removing the Bulk Flow has on the momentum power spectrum and its relation to the well-known integral constraint, we conclude that this should be tested on a case-by-case basis.

In addition to the Bulk Flow, peculiar velocity surveys often suffer from an unknown ‘zero-point’ and the monopole of the velocity field is artificially set to zero. Howlett et al. (2017a) demonstrate that the impact of a constant addition to the velocity field would be to add as shot-noise to the velocity power spectrum. The same would be true for the momentum power spectrum. However, even for the optimistic sample we consider in this work, the contribution to the shot-noise from a zero-point offset would remain subdominant compared to the shot-noise arising from the distance errors. Hence we believe that any systematic error associated with the monopole of the velocity field will be small compared to the typical statistical error on the momentum power spectrum.

6.3 Gaussianity of power spectrum measurements

One final thing that we wish to test before applying this work to real data is the Gaussianity of the distribution of momentum power spectrum measurements. This affects how well we can describe the data using only the measurements themselves and an estimate of their covariance matrix; and how easily we are able to recover cosmological constraints from these. We use our mock catalogues as a test of the Gaussianity of the moments of the momentum power spectrum with and without PV errors by evaluating the difference between the distribution of the $P^{\,p}_{\ell}(k)$ values and the mean at fixed $k$. Our results are shown in Fig. 12, where we compute the skewness and kurtosis of this distribution compared to a true Gaussian with the same number of measurements666More precisely, we draw a number of ‘fake’ momentum power spectra from a multivariate Gaussian distribution equal to the number of mocks, with mean and covariance equal to the measured mean and covariance. We then compute the skewness and kurtosis. We repeat this 1000 times, obtaining 1000 values for the skewness and kurtosis in each $k$-bin. We then plot the $16^{th}$ and $84^{th}$ percentiles of the skewness and kurtosis values.. We also plot the distribution for three individual $k$ bins.

For the measurements without errors, we find the monopole and quadrupoole of the momentum power spectrum to be closer to log-normally distributed rather than Gaussian (similar results are found for the galaxy density power spectrum e.g., Ross et al. 2015) resulting in considerable skewness and kurtosis. If we instead perform the experiment above on the logarithm of the monopole, the distribution becomes close to Gaussian; there is some residual skewness and kurtosis, but this is not particularly evident in the individual distributions and we do not believe it will have much impact on parameter constraints if the logarithm is fit. However, such a transformation cannot be naively applied to the quadrupole due to the presence of negative $P^{\,p}_{2}(k)$ values. Once PV measurement errors are included, the distribution of monopole and quadrupole power spectrum measurements becomes much more Gaussian given that the PV errors are also Gaussian (seen by comparing the left and right panels of Fig 12). On large scales the distribution still contains significant skewness and kurtosis which, again, we cannot get round using the logarithm of the measurements due to the presence of negative $P^{\,p}(k)$ values in both the monopole and quadrupole.

Fortunately, this is not unmanageable. When fitting the mocks in Section 7 we find that this does not overly bias the growth rate recovered from fitting the momentum power spectrum. In Paper II we also found no bias in the growth rate constraints due to non-Gaussianity in measurements from the 6dFGSv survey. For the 2MTF data used in Paper II we did find that the small volume of the survey leads to a high degree of non-Gaussianity which biases the growth rate measurements777As highlighted in Kalus et al. (2016), for a field whose real and imaginary parts are close to individually Gaussian distributed (as might be expected for our Universe) the power spectrum is actually Rayleigh distributed. The distribution of measured power spectra only approaches Gaussian if there are sufficient modes that the central limit theorem holds. For the small volume of the 2MTF survey, the number of modes is sufficiently low that the central-limit theorem does not hold even for quite non-linear scales.. However, in this same work, we developed and verified that performing a simple transformation to the measurements from the mocks and data prior to fitting is enough to fully account for this. Hence even though we do not apply them to our mocks, the techniques therein could be applied to the momentum power spectrum measurements from any dataset. We conclude that, given an estimated covariance matrix, we can easily extract cosmological information from measurements of the momentum power spectrum from real data. The presence of Gaussian PV errors (even small ones used here) is enough to ensure that whilst the intrinsic measurement in roughly log-normal, the measurement from data is Gaussian enough for robust constraints to be obtained.

7 Final proof of concept

As a final proof of concept that measurements of the momentum power spectrum can be used to access the cosmological information in the velocity field from current and next-generation surveys, we bring everything discussed so far in this work together and measure the growth rate of structure from our mock PV surveys introduced in Section 6.

We take measurements of the monopole, quadrupole and hexadecapole ($\ell=0$,$2$ and $4$ moments) of the galaxy density field of the mock galaxies using an FKP weight with constant $P_{FKP}=8000\,h^{-3}\,{\rm Mpc^{3}}$ and measurements of the monopole and quadrupole of the momentum power spectrum using weights with constant $P^{\,p}_{FKP}=2.4\times 10^{9}\,{\rm km^{2}\,s^{-2}}\,h^{-3}\,{\rm Mpc^{3}}$. These weights were chosen to maximise the information on the growth rate based on the procedure given in Appendix B, where we show how to estimate the Fisher information on the growth rate contained in different multipoles of the density and momentum power spectrum using different weights. The same calculation is used to produce forecasts for the constraints on the growth rate for a scenario exactly matching the measurements we make here.

The measurements are then combined with the redshift-space model given in Appendix A. We correct this model for the effects of the window function using a modified version of the code provided with Blake et al. (2018), which computes a convolution matrix that multiplies the fine-binned model multipoles and returns a convolved model with the same binning as the measurements. We correct for the integral constraint using the same code, and remove the Bulk Flow before making the measurements, consistent with our findings in Section 6.2. We then fit the average measurements from the mocks for $k=0.02-0.20\,h\,{\rm Mpc}^{-1}$ using a Gaussian likelihood and covariance matrix also estimated from the mocks.

From Appendix A, there are five free parameters for the density and momentum power spectra: $b_{1}\sigma_{8}$, $b_{2}\sigma_{8}$, $f\sigma_{8}$, and the velocity dispersions $\sigma_{v,1}$ and $\sigma_{v,2}$. This models is found to recover unbiased constraints when fitting either the density field or momentum field multipoles separately. However, when fitting all five multipoles simultaneously, we find that we need to use separate $\sigma_{v,1}$ parameters for the density and momentum measurements. We attribute this to the fact that we are weighting the same sample differently for the two measurements, which changes the effective depth and bias of the sample. However, as the bias and velocity dispersion terms all enter the momentum power spectrum measurements on non-linear scales, we are able to account for this with a single extra parameter (rather than also introducing a second set of bias parameters). Hence when fitting both the density and momentum power spectra simultaneously we have a total of 6 free parameters.

We start looking at how different components of the theoretical redshift-space momentum power spectrum go together to fit our data. Figure 13 shows best-fitting model for the monopole and quadrupole of the momentum power spectrum alongside the measurements and errors. We identify different components contributing to the best-fit model; the velocity power spectrum (the terms in $P_{11}$ in Appendix A independent of galaxy bias), the power spectrum of the convolution between the density and velocity field (the $P_{11}$ terms containing galaxy bias), and non-linear redshift-space distortions ($P_{12}$, $P_{13}$ and $P_{22}$). As expected from Fig. 2, for the monopole the velocity power spectrum dominates on large scales, and the density-velocity convolution spectrum on non-linear scales. The non-linear RSD terms stay subdominant on all the scales we fit. For the quadrupole, the story is a little different. On the largest scales we consider, the velocity power spectrum remains dominant, but non-linear redshift-space distortions become important on relatively large scales compared to the monopole. The contributions from the convolution between the density and velocity fields and non-linear RSD are similar for $k>0.1\,h\,{\rm Mpc}^{-1}$. Given that both the monopole and quadrupole are proportional to the growth rate on large scales, and have different dependencies on bias and non-linear RSD on small scales, we might expect there to be some improvement when adding the momentum power spectrum quadrupole. However, the effects of galaxy bias and non-linear RSD combined with the noise on non-linear scales means these are likely very degenerate and that the real benefit comes from having an independent measurement of the galaxy bias from the density power spectrum.

To explore this and the information contained in other combinations of measurements of the galaxy density and momentum multipoles we turn next to the growth rate constraints themselves. These are listed in Table 1 for various combinations of moments of the combined density and momentum power spectra. For each case we also include forecasts based on the method in Appendix B for the growth rate measurements assuming truly optimal weights and using weights with constant $P(k)$ as given above. We also show the $b_{1}\sigma_{8}-f\sigma_{8}$ contours for the density field measurements and then when the momentum power spectrum is added in Fig. 14. We found negligible information in the $\ell>4$ and $\ell>2$ moments of the density and momentum power spectra respectively, which is corroborated by our forecasting code.

To begin, the growth rate constraints are largely consistent with the expected value for the mocks although with a small bias likely due to non-linear modelling inaccuracies or residual non-Gaussianity in the measurements (the latter of which is solved in Paper II, but not a concern for our proof of concept). The first point of note is that the measurements using only the monopole of the momentum power spectrum are tighter than using all three moments of the galaxy density power spectrum for our sample. This is perhaps an unfair comparison as our sample contains equal numbers of redshifts and PV, whereas in reality a substantially higher number of redshifts could be obtained for the same ‘effort’ as the 150,000 PVs we use here. Nonetheless, the fact is that we would expect the information in the density field alone to saturate relatively quickly as the cosmic variance limit is reached, and the momentum power spectrum is proportional to the growth rate on linear scales. This means we get tighter constraints on the growth rate than with the galaxy density power spectrum and it is far easier to model and extract that information from the monopole of the momentum power spectrum than the full set of galaxy density multipoles.

The real benefit comes from combining these two sets of measurements. The combined constraints from the full set of 5 power spectrum multipoles is a factor of $\sim 1.5$-$2$ times better than the constraints from either the density or momentum field alone. We find that most of the constraints come from the first two moments of the galaxy density power spectrum and the monopole of the momentum power spectrum. Including higher order moments in the combination does provide some advantage but in practice this may not be worth the additional modelling complexity or systematic errors this could introduce.

It is also interesting to ask how much of the total information in the velocity field is recovered using the momentum power spectrum (which is ‘contaminated’ by bias on small scales). We can answer these by looking at the forecasts in Table 1. To provide a benchmark, we first compare the forecasts for the optimal weights and using the fixed $P_{FKP}$ weighting scheme. The use of fixed values to make the measurements reduces the information by about $10\%$ and $30\%$ for the density and momentum power spectrum respectively. We additionally remind the reader that the fixed values were chosen as the weights that minimize this difference. Hence we conclude that the use of a fixed weight impacts the constraints somewhat, and this commonly used technique could be improved upon. This is true for the galaxy density power spectrum (and so likely true for a number of published works), but is more evident in the momentum power spectrum due to the fact that this has a strong $1/k^{2}$ scale dependence. Future improvements to this work could develop a method to implement weights that vary as a function of $k$ which would benefit even measurements of the galaxy density power spectrum from redshift-only surveys.

Looking at the case for the momentum power spectrum alone, we see the measurements to be quite a bit worse than the forecasts even for non-optimum weights. Some portion of this can be attributed to the overoptimistic nature of the Fisher matrix forecasts as the measurements from the density field alone are also larger than the forecasts by $\sim 25\%$, however this will also be due to the fact that the momentum power spectrum is not the velocity power spectrum and so does not contain much growth rate information on non-linear scales where the density-velocity convolution starts to dominate ($k\approx 0.1\,h\,{\rm Mpc}^{-1}$ as seen in Fig. 13). However, when combining the density and momentum power spectra the story is different. The constraints we obtain are in-line with the forecasts for our weighting scheme; simultaneously fitting the density and momentum measurements serves to break the degeneracy between galaxy bias and the growth rate in both of these. This proves that our method for measuring and modelling the combined density and velocity field information using the momentum power spectrum does an excellent job of recovering the information on the growth rate from the PV data. This serves as a proof-of-concept that the work in this paper can be used to extract accurate constraints from real data.

A final question we can ask is how much information we may be missing by only considering the power spectra of the density and momentum field and their cross-covariance and neglecting the cross-spectrum between these two fields. This was mentioned in Section 4.1, where we argued that measuring this introduces enough complexity that we leave it for future work. Although we also leave the full details of the calculation to the next paper in this series, we can include the multipoles of the cross-spectrum, its covariance matrix, and its cross-covariance with the multipoles of the galaxy density and momentum power spectra in the Fisher matrix method of Appendix B. Doing so, we forecast errors on $f\sigma_{8}$ of 0.026 and 0.030 for the optimum and fixed weights used herein. This is an improvement of $\sim 10\%$ over the case where we do not include this observable, which is somewhat less than the results of Adams & Blake (2017). However, as they pointed out in their comparisons with simulations, a larger improvement may be seen in data with lower number density or worse velocity errors, where the growth rate and galaxy bias are less well constrained. A larger improvement may be also be seen in the case of real data where including the cross-spectrum could help resolve unknown systematics.

8 Conclusions

In this work we have formally derived an estimator for the redshift-space momentum (or mass-weighted velocity) power spectrum from a peculiar velocity (PV) survey. The method holds for any tracer of the velocity field, including galaxy motions from Tully-Fisher or Fundamental Plane surveys, or for Type-IA supernovae. This derivation mimics the well-known Feldman et al. (1994) (FKP) and Yamamoto et al. (2006) estimators for the galaxy density power spectrum, and has a number of useful properties. In particular:

• •

It is an extremely efficient compression of information compared to other techniques for extracting the velocity power spectrum (Johnson et al., 2014; Adams & Blake, 2017; Howlett et al., 2017b; Huterer et al., 2017).

• •

The normalisation of the galaxy and momentum power spectrum is the same.

• •

The effect of a survey window function introduces a convolution with the momentum power spectrum identical to the galaxy density power spectrum. This can be best seen by comparing Eqs. 14-16 with the expressions for the convolution of the galaxy density power spectrum in Feldman et al. (1994), Yamamoto et al. (2006) or Blake et al. (2018).

• •

The momentum power spectrum is also subject to shot-noise dependent on the galaxy number density, however there is an additional contribution from the variance in the velocity field, which will contain both an intrinsic component, and a component due to measurement error.

• •

The covariance matrix of the momentum power spectrum, and the cross-covariance with the galaxy density power spectrum can both be calculated analytically in the Gaussian regime. We find good agreement between the Gaussian prediction and measurements from simulations on large scales, but the scale at which this breaks down is much larger than for the galaxy density power spectrum due to the convolution of density and velocity fields which enters the momentum power spectrum.

• •

We have computed the optimal weights for the momentum power spectrum. Much like the FKP weights, these depend on the galaxy number density and power spectrum itself. We also find dependence on the intrinsic variance of the velocity field.

We have also investigated several properties unique to PV surveys that affect the momentum power spectrum measurements. We find that the presence of errors on individual PV measurements adds as shot-noise to the momentum power spectrum, such that our estimator remains unbiased, but the variance increases significantly. Removing the bulk motion of the observer before estimating the momentum power spectrum removes power on large scales, acting as an integral constraint, but also reduces the variance accordingly. In the mocks used here, we find it beneficial to remove the bulk flow when fitting the growth rate of structure as it slightly reduces the uncertainty on the constraints but without adding bias. However we caution that this result could arise from of the relatively small bulk flow we might expect in our mocks compared to current datasets and that more theoretical work is needed to understand how this changes the measured momentum power spectrum. Finally, we show that the momentum power spectrum measured from a galaxy survey is Gaussian distributed on scales sufficiently smaller than the survey volume, allowing for extraction of cosmological constraints using only an estimate of the covariance matrix. This point is further investigated in Paper II where we validate a method to recover unbiased cosmological constrains even in the presence of greater levels of non-Gaussianity than shown here.

Finally, we have provided a complete proof of concept demonstrating that our estimator of the momentum power spectrum makes the information in the velocity field, in particular on the growth rate of structure, readily obtainable. This will be a key tool in the analysis of the next generation Taipan Galaxy survey (da Cunha et al., 2017), WALLABY survey (Koribalski, 2012) or of Type-IA supernovae (Howlett et al., 2017c; Kim et al., 2019) and will allow these surveys to overcome the limitations of cosmic variance in the low redshift universe, leading to some of the most precise test of gravity to date. However, this work is only the tip of the iceberg; in Paper II (Qin et. al., in preparation), we present the first application of all the techniques in this work to the combined 2MASS Tully-Fisher (Hong et al. 2014; Hong et. al, in preparation) and 6-degree Field Galaxy Redshift Survey peculiar velocity samples (Campbell et al., 2014). In subsequent papers in this series we will extend this to look at the remaining information contained within in the density-momentum cross-power spectrum and the prospects of constraining modified gravity theories.

Acknowledgements

We thank Fei Qin, David Parkinson, Lister Staveley-Smith and Chris Blake for many useful discussions and providing comments on this manuscript. We especially thank the latter for providing their code to measure and model the multipoles of the galaxy density power spectrum.

This research was conducted by the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), through project number CE110001020. This work was performed on the Pleiades HPC cluster at the International Centre for Radio Astronomy Research at the University of Western Australia and on the OzSTAR and swinSTAR supercomputers at Swinburne University of Technology.

This research has made use of NASA’s Astrophysics Data System Bibliographic Services and the astro-ph pre-print archive at https://arxiv.org/. All plots in this paper were made using the matplotlib plotting library (Hunter, 2007).

Appendix A Theoretical model for the redshift-space density and momentum power spectrum

In this work we model the redshift-space density and momentum power spectrum for biased tracers using the combined distribution function and Eulerian perturbation theory approach of Vlah et al. (2012, 2013); Okumura et al. (2014); Saito et al. (2014), with some small modification to higher order terms for easier computation. Full details of the model can be found therein, but requires some work piecing together various components; here we provide a more condensed overview.

We begin by defining the multipoles of the power spectra as integrals with respect to the angle, $\mu$

[TABLE]

where the quantities $P^{\delta}(k,\mu)$ and $P^{p}(k,\mu)$ are the full redshift-space anisotropic galaxy density and momentum power spectra respectively. We model the two anisotropic power spectra (dropping the explicit dependencies on $z$, $k$ and $\mu$) as

[TABLE]

where

[TABLE]

In the above equations, $D(z)$ is the linear growth factor, $f(z)$ is the growth rate of structure, $H(z)$ is the Hubble parameter, $a(z)$ is the scale factor, $P_{L}(k)$ is the $z=0$ linear matter power spectrum and $\sigma_{v}$ is the velocity dispersion. We have used the galaxy bias expansion of McDonald & Roy (2009) which encapsulates local and non-local bias up to third order in the four terms $b_{1}$, $b_{2}$, $b_{s}$ and $b_{3nl}$. The remaining terms $I_{mn}$, $J_{mn}$, $K_{mn}$, $\sigma^{2}_{3}$ and $\sigma^{2}_{4}$ are integrals over the linear power spectrum.

$I_{mn}$ and $J_{mn}$ are neatly presented in Appendix D of Vlah et al. (2012) and so we do not repeat them here. We can write the terms $K^{(s)}_{mn}$ in the same form

[TABLE]

where $r=q/k$, $x=\frac{\boldsymbol{k}\cdot\boldsymbol{q}}{kq}$ and

[TABLE]

We normalise $K_{01}\rightarrow K_{01}-\int\frac{d^{3}q}{(2\pi)^{3}}P^{2}_{L}(q)$, $K^{s}_{01}\rightarrow K^{s}_{01}-\frac{4}{9}\int\frac{d^{3}q}{(2\pi)^{3}}P^{2}_{L}(q)$ and $K^{s}_{02}\rightarrow K^{s}_{02}-\frac{2}{3}\int\frac{d^{3}q}{(2\pi)^{3}}P^{2}_{L}(q)$, which ensures that the higher order bias terms go to zero as $k\rightarrow 0$ (Saito et al., 2014).

The final pieces we need are

[TABLE]

A.1 Choice of free parameters

In the above equations there are a number of free parameters. In addition to the cosmological parameters that affect the shape of the linear power spectrum, these include the growth rate of structure $f$ (our main parameter of interest), galaxy bias parameters $b_{1}$, $b_{2}$, $b_{s}$ and $b_{3,nl}$, the normalisation of the linear power spectrum commonly parameterised as $\sigma_{8}$ and the velocity dispersion $\sigma_{v}$. However, measurements of the density and momentum power spectrum are too noisy to break the degeneracy between all these parameters. Instead, when fitting our mock data we fix the linear power spectrum (and hence the $I_{mn}$, $J_{mn}$ and $K_{mn}$) to the input cosmology of the mocks with $\sigma_{8}=0.815$ and adopt the parameter combinations $b_{1}\sigma_{8}$, $b_{2}\sigma_{8}$, $f\sigma_{8}$. We then write the remaining galaxy bias terms as a function of the linear bias $b_{s}=-4/7(b_{1}-1)$ and $b_{3,nl}=32/315(b_{1}-1)$ (Saito et al., 2014). A similar technique was used in Beutler et al. (2017). This allows for constraints to be put on $f\sigma_{8}$ from the overall amplitude of the momentum power spectrum (which still allows for tests of different gravitational models; Song & Percival 2009) whilst the systematic effects of fixing $\sigma_{8}$ can be absorbed into the nuisance terms causing negligible change in the recovered values for the growth rate.

The velocity dispersion can be evaluated on linear scales as

[TABLE]

however as demonstrated in Vlah et al. (2012) this will not fully capture the non-linear motions of galaxies within their host halos. In applying the above model to our mock galaxy catalogues in Section 7 we found that even leaving $\sigma_{v}$ to vary leads to significant biases in the recovered values of $f\sigma_{8}$. Instead, motivated by their Table 1 we make the following redefinitions

[TABLE]

We are hence swapping a single $\sigma_{v}$ (free or otherwise) for 2 free parameters $\sigma_{v,1}$ and $\sigma_{v,2}$ that enable the model to more accurately reproduce the non-linearities in the motions of galaxies. Overall we fit for 5 parameters [$b_{1}\sigma_{8}$, $b_{2}\sigma_{8}$, $f\sigma_{8}$, $\sigma_{v,1}$, $\sigma_{v,2}$].

Appendix B Fisher information in the combined velocity and density field multipoles

In this appendix we develop a simple Fisher matrix calculation for the cosmological information contained in measurements of the galaxy density or velocity power spectrum multipoles for arbitrary weights. This is based on a combination of Taruya et al. (2011), Koda et al. (2014) and Howlett et al. (2017a). Taking results from these works, the Fisher information for a set of parameters $\theta$ contained within a series of power spectrum multipoles $P_{\ell}(k)$ can be calculated as

[TABLE]

where $\Omega$ is the sky area of the survey (in steradians) and the $\mathsf{C}_{\ell\ell^{\prime}}(k,w)$ is the covariance matrix between the different measured multipoles which depends on the weights $w$ assigned to the data. In this work we consider multipoles of the galaxy density and momentum power spectra (with some consideration also given to the cross-spectrum). The covariance of these various measurements has been derived in various Sections of this work and can be pieced together to create $C_{\ell\ell^{\prime}}(k,w)$ depending on the exact combination of observed multipoles. Assuming that the power spectra depend only on $k$ and $\mu$ and that the weights, number density and velocity errors of galaxies vary only radially, we can write these as a modified version of Eq. 40,

[TABLE]

where we have carefully maintained the distinction between the weights and number density used in the estimation of the density and momentum power spectra (denoted via the $\delta$ and $p$ superscripts respectively). $r_{max}$ is the maximum comoving distance of the survey; in solving the above equations it is often useful to perform a change of variable from comoving distance to redshift, such that $r_{max}\rightarrow z_{max}$.

The other piece we need to compute the Fisher matrix is the derivatives of the power spectrum multipoles as a function of cosmological parameters. For our Fisher forecasts in Table 1 we are interested in the information in the velocity $u$, power spectrum (rather than the momentum $p$, power spectrum, which may contain less) and so instead of the models in Appendix A, we use models for the density, velocity and cross power spectrum from Koda et al. (2014); Howlett et al. (2017a)

[TABLE]

where the damping terms $D_{g}(k,\mu)=[1+0.5(k\mu\sigma_{v,1})^{2}]^{-1/2}$ and $D_{u}(k)=\mathrm{sinc}(k\sigma_{v,2})$ account for non-linear redshift-space distortions (although not to the same extent as the model in Appendix A). The models in Eq. 68 are then integrated with respect to the Legendre polynomials to obtain the multipoles. The derivatives with respect to the parameters (in this case galaxy bias, growth rate and two velocity dispersion parameters for 4 in total) can be performed before the integration. The above equations for the covariances still hold when we swap out the momentum power spectrum for the velocity power spectrum. Given these derivatives and the above equations we can compute Eq. 64 relatively simply for arbitrary combinations of multipoles. Finally, the Fisher matrix can be inverted to get the covariance matrix for the parameters, where the diagonal element corresponding to the growth rate is the predicted variance after marginalising over the other parameters in the matrix. In applying this to our mock catalogues we use the number density as a function of redshift averaged over all our mocks and marginalise over $b_{1}\sigma_{8}$, $\sigma_{v,1}$ and $\sigma_{v,2}$ when computing the growth rate errors.

The above method also allows us to numerically minimise the predicted error on the growth rate as a function of the two numbers $P^{p}_{FKP}$ and $P_{FKP}$ used as inputs to Eq. 41 and the standard FKP weight. When doing this we are interested in actually modelling the momentum power spectrum (not just the velocity component, as the optimal constant values should take into account the fact that the density-velocity convolution term enters on non-linear scales) and so we do use the models of Appendix A. We compute the derivatives using these models and then substitute the weights for each $P^{p}_{FKP}$ and $P_{FKP}$ into the above covariance matrices. Performing this optimisation for the mock catalogues in Section 7 returns rounded values $P_{FKP}=6000\,h^{-3}\,{\rm Mpc^{3}}$ and $P^{p}_{FKP}=2.4\times 10^{9}\,{\rm km^{2}\,s^{-2}}\,h^{-3}\,{\rm Mpc^{3}}$. Alternatively, the Fisher matrix could be numerically optimised for a particular choice of $k$ and $\mu$, such that the model anisotropic density and momentum power spectra are first calculated for this $k$ and $\mu$, and then these two values are used to compute the fixed weights.