© The Authors 2026

1. Introduction

One of the most profound scientific pursuits in cosmology is to unravel the evolutionary processes that have shaped the large-scale structures of our Universe. Galaxies exhibit peculiar motions – superimposed on their Hubble recession velocities – due to gravitational fluctuations arising from local variations in the matter density field. These peculiar velocities serve as a crucial observational probe, offering deep insights into the underlying dark matter distribution and the dynamic history of cosmic expansion and structure formation. Consequently, they have garnered growing interest as an essential diagnostic tool for assessing the accuracy and predictive power of cosmological models.

An approach to assess the cosmological models involves measuring the low-order kinematic moments of the cosmic flow field based on the radial components of galaxies’ peculiar velocities, and subsequently comparing these estimates with the predictions of the cosmological models (Kaiser 1988; Strauss & Willick 1995; Feldman et al. 2010; Scrimgeour et al. 2016; Zhang et al. 2017; Qin et al. 2018, 2019b; Qin 2021; Qin et al. 2021; Whitford et al. 2023). However, due to the significant cosmic variance inherent in low-redshift velocity surveys, the constraints imposed on cosmological models through this method are typically not tight.

An alternative approach involves estimating cosmological parameters and subsequently comparing them with predictions of cosmological models. In the field of cosmology, researchers commonly aim to investigate the rate at which large-scale structures evolve, as well as the ratio between galaxy density and dark matter density. The former quantity is associated with the structure linear growth rate, f, while the latter relates to the galaxy biasing parameter, b. Furthermore, in practical applications, normalized parameters fσ₈ and bσ₈ have attracted attention as measurable quantities. σ₈ refers to the root mean square of mass density fluctuations within spheres of 8 Mpc h⁻¹. Various methodologies have been proposed in prior studies to extract these parameters from galaxy survey data (Ishak 2019; Ishak et al. 2025).

One technique involves the reconstruction of the cosmological density and velocity fields (Nusser et al. 1991; Kitaura et al. 2012; Wang et al. 2012), or the cosmography of the local Universe (Courtois et al. 2013; Dupuy & Courtois 2023). The parameters fσ₈ and bσ₈ can be estimated by comparing the observed peculiar velocities of galaxies with the reconstructed velocity field (Carrick et al. 2015; Said et al. 2020; Boruah et al. 2022). However, systematic errors associated with this method are difficult to quantify, and previous studies have indicated that the errors tend to be underestimated (Turner & Blake 2023; Blake & Turner 2024). Recently, artificial intelligence techniques have increasingly been applied to field reconstruction (Hong et al. 2021; Qin et al. 2023a; Wu et al. 2023), offering novel opportunities in this research area. Moreover, since the covariance matrix of the density and velocity fields can be formulated in terms of fσ₈ and bσ₈, these parameters can alternatively be estimated by maximizing the likelihood function under the assumption of a Gaussian distribution (Johnson et al. 2014; Adams & Blake 2017; Lai et al. 2023; Carreres et al. 2023; Rocher et al. 2023; Ravoux et al. 2025). However, this method remains computationally demanding. Another commonly used approach is the two-point statistic method, which has played a significant role in cosmology over the past decades. Based on this method, fσ₈ and bσ₈ can be derived from the galaxy two-point correlation functions, the velocity-correlation functions, and the galaxy-velocity cross-correlation functions (Gorski et al. 1989; Howlett et al. 2015; Dupuy et al. 2019; Wang et al. 2018; Turner et al. 2021, 2023; Qin et al. 2022, 2023b; Lyall et al. 2023). Nevertheless, these approaches tend to be computationally intensive too.

Instead of studying the two-point correlations in configuration space, we can also analyse them in Fourier space, i.e. analysing the so-called power spectrum. The concept of the cosmological power spectrum was first introduced by Kaiser (1987). The first measurement of the (galaxy) density power spectrum was carried out by Kaiser & Peacock (1991). The methodology for estimating the density power spectrum was later refined by Feldman et al. (1994), forming the standard approach widely used today. In addition to these foundational contributions, Yamamoto et al. (2006) and Bianchi et al. (2015) further developed estimators for the multipoles of the density power spectrum.

However, the linear growth rate of the structure, fσ₈ is only primarily sensitive to higher-order components of the density power spectrum, resulting in relatively weak constraints when inferred from this method alone. To improve the estimation of fσ₈, the momentum power spectrum – initially explored in Park (2000) – has recently gained widespread use in cosmology. Howlett (2019, hereafter Paper I) adapted the momentum power spectrum method to be implemented in a manner analogous to the conventional density power spectrum approach. Qin et al. (2019a, hereafter Paper II) successfully applied the momentum power spectrum to constrain fσ₈ using data from the combined 6dFGS peculiar velocity survey (6dFGSv; Campbell et al. 2014) and 2MASS Tully-Fisher (2MTF; Hong et al. 2019) survey. Building upon their research, Qin et al. (2025, hereafter Paper III) advanced the existing methodology by extending it to the cross-power spectrum of density and momentum fields – both derived from a unified survey catalogue. Their work achieved a successful estimation of fσ₈ through the peculiar velocity data from the Sloan Digital Sky Survey peculiar velocity catalogue (SDSSv; Howlett et al. 2022). This paper marks the fourth in the ongoing series. Herein, we further refine and expand our analytical framework to model and measure the cross-power spectrum in more complex cases, where the density and momentum fields are drawn from distinct galaxy survey catalogues, each exhibiting unique survey geometries, selection criteria, and sample characteristics. Additional studies employing the momentum power spectrum include those by Park & Park (2006), Appleby et al. (2023) and Shi et al. (2024).

We constrained fσ₈ using Data Release 1 (DR1) from the Bright Galaxy Survey (BGS; Hahn et al. 2023) and Peculiar Velocity Survey (DESI-PV; Carr et al. 2025; Ross et al. 2025; Douglass et al. 2025) of the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration 2025a,b). The DESI-PV DR1 stands as the largest peculiar velocity survey ever assembled to date. This paper forms part of a comprehensive series of results from DESI-PV DR1 (Turner et al. 2025; Lai et al. 2025; Ross et al. 2025; Douglass et al. 2025; Bautista et al. 2025; Carr et al. 2025; Nguyen et al. 2025).

This paper is organized as follows. Section 2 presents the data and mocks employed for measuring the power spectrum. Section 3 introduces the estimators utilized to compute the power spectrum. Section 4 delves into the theoretical models of the power spectrum. In Sect. 5 we outline the fitting methodologies and validate our estimators and models through extensive testing on mocks. Section 6 showcases the final parameter fitting results derived from the survey data. In Sect. 7 we discuss the parameter constraints in light of predictions and fits from other data sets. Lastly, Sect. 8 offers a comprehensive conclusion summarizing the key findings of the study.

In this paper, we adopt the AbacusSummit (Maksimova et al. 2021) cosmological parameters¹n_s = 0.96490, A_s = 2.101198 × 10⁻⁹, Ω_m = 0.315192, Ω_bh² = 0.02237, Ω_ch² = 0.12000, σ₈ = 0.8113545, and H₀ = 100 h km s⁻¹ Mpc⁻¹, where h = 0.67360. We adopt a flat Λ cold dark matter (ΛCDM) fiducial cosmological model.

2. Data and mocks

2.1. The Dark Energy Spectroscopic Instrument

Located at the Kitt Peak National Observatory in Arizona, USA, the Dark Energy Spectroscopic Instrument is a multi-fibre spectrograph mounted on the 4-meter Mayall Telescope, designed to explore the large-scale structure of the Universe and the expansion history of the Universe (Levi et al. 2013; DESI Collaboration 2016a,b, 2022; Miller et al. 2024; Poppett et al. 2024). The DESI focal plane spans 8 square degrees and is equipped with 5000 optical fibres, each of which can be positioned onto individual targets within an observational field using robotic positioners (Silber et al. 2023). The spectrographs span the near-ultraviolet to near-infrared spectrum (3600–9800 Å), delivering a spectral resolution that ascends from 2000 in the blue camera to 5000 in the red camera. The observing strategy is detailed in Guy et al. (2023) and Schlafly et al. (2023). Following its Early Data Release (EDR; DESI Collaboration 2024a,b), the survey has now progressed to DR1 (DESI Collaboration 2025a,c). For this study, the galaxy sample was drawn from the DESI DR1 catalogue. We used the DESI DR1 large-scale structure catalogue² (DESI Collaboration 2025e).

2.1.1. The BGS data

The Bright Galaxy Survey (BGS; Hahn et al. 2023) Bright sample has undergone refinement by Turner et al. (2025), who introduced tailored selection criteria to enhance its suitability for analyses involving peculiar velocities. The final BGS data catalogue employed in our study encompasses a total of 415 523 galaxies, with an r-band absolute magnitude threshold set at −17.7. The sky distribution of these galaxies is illustrated by the black dots in Fig. 1. Furthermore, the redshifts of these galaxies, measured in the CMB frame, span the range from z = 0.01 to z = 0.1, as is depicted by the black bars in Fig. 2. The galaxy mean number density, $\overline{n}(\mathbf{r})$, is shown in the top left panel of Fig. 3. The calculation of $\overline{n}(\mathbf{r})$ is detailed in Bautista et al. (2025). In this paper, the density field was derived from the BGS catalogue.

2.1.2. The DESI-PV data

Due to the peculiar motions of galaxies, the apparent distance of a galaxy, denoted as d_z, deviates from its true comoving distance, d_h. This discrepancy is measurable and can be quantified through the logarithmic distance ratio, defined as $\eta \equiv {log}_{10}\frac{d_z}{d_h}$. It can be converted to the line-of-sight peculiar velocity using the velocity estimator Eq. (A.1). Within the DESI-PV survey, two distinct populations of galaxies are observed: spiral galaxies, also known as late-type galaxies, and elliptical galaxies, referred to as early-type galaxies. The η for late-type galaxies is derived from the Tully-Fisher (TF; Tully & Fisher 1977) relation, whereas for early-type galaxies it is inferred from the fundamental plane (FP; Djorgovski & Davis 1987, Dressler et al. 1987) relation. In DESI-PV, the TF sample (Douglass et al. 2025) comprises a total of 6,806 late-type galaxies. The sky distribution of these galaxies is illustrated by the blue dots in Fig. 1. Furthermore, the redshifts of these galaxies span the range from z = 0.015 to z = 0.1, as depicted by the blue bars in Fig. 2. In contrast, the FP sample (Ross et al. 2025) consists of 73 822 early-type galaxies. Their sky distribution is marked by grey dots in Fig. 1, and their redshifts span the interval from z = 0.01 to z = 0.1, as is shown in the grey bars of Fig. 2.

The DESI-PV catalogue is a combination of the TF sample and the FP sample. The zero-point calibration bridging these two datasets has been rigorously investigated by Carr et al. (2025), ensuring a physical alignment. The redshift distribution of this combined sample is illustrated by the yellow bars in Fig. 2. Meanwhile, the galaxy mean number density, $\overline{n}(\mathbf{r})$, is depicted in the top right panel of Fig. 3. In this study, the momentum field was constructed from the DESI-PV catalogue. In accordance with the reasoning articulated in Turner et al. (2025), we imposed a redshift cut of z_cut = 0.05 on the TF galaxies when measuring the power spectrum – excluding those galaxies lying beyond this threshold – to mitigate systematic biases inherent in the TF velocity measurements. This identical redshift cut was consistently applied across both mocks and random catalogues, ensuring coherence and comparability throughout the analysis. Following Turner et al. (2025), a 4σ-clipping was applied to the measured η values of the DESI-PV data too.

2.2. Mocks

A considerable number of mocks is indispensable for validating the algorithm and accurately estimating the uncertainties in the power spectrum and associated parameters. Bautista et al. (2025) offers a comprehensive set of 675 BGS mocks, meticulously designed to mirror the characteristics of the BGS data. In parallel, 675 FP mocks have been constructed to emulate the FP data, alongside 675 TF mocks tailored to align with the TF data. Each set of BGS, FP, and TF mocks shares the same observer, resulting in a total of 675 mock sets. These mocks were derived from the AbacusSummit simulation, and the fiducial value of fσ₈ for these mocks is 0.466 at an effective redshift of z = 0.2.

The FP and TF mocks were generated with the specific aim of replicating the selection function and survey geometry, the FP relation, the TF relation, and the galaxy clustering observed in the FP and TF data, respectively. The DESI-PV mocks were formed by combining the TF and FP mocks. The galaxy mean number density, $\overline{n}(\mathbf{r})$, of DESI-PV mocks is illustrated in the bottom right panel of Fig. 3. The galaxy mean number density, $\overline{n}(\mathbf{r})$, of BGS mocks is illustrated in the bottom left panel of Fig. 3.

2.3. The random catalogues

The random catalogues for both the data and mocks are also supplied by Bautista et al. (2025). In order to precisely extract cosmological parameters from the power spectrum, it is imperative to appropriately account for the survey geometry of both BGS and DESI-PV. Employing fast Fourier transforms to measure the power spectrum in the context of a non-periodic and incomplete survey geometry leads to a convolution of the true power spectrum with the survey’s window function. To address this, the random catalogue was employed to estimate the window function, thereby enabling the application of an equivalent convolution to our theoretical model during the data fitting process. The random catalogues for FP and TF were generated independently, with the objective of more faithfully replicating the distinct observational features of the TF and FP surveys, respectively. The random catalogue for DESI-PV was constructed by combining the random catalogues of TF and FP.

3. The estimation of the power spectrum

3.1. The field functions

In this paper, we investigate the redshift-space power spectrum of galaxies, encompassing both the auto-power spectrum and cross-power spectrum of the galaxy density and momentum fields. The density field was derived from the BGS catalogue, whereas the momentum field was obtained from the DESI-PV catalogue. Within our paper series, the density (contrast) field is defined by

$$\begin{array}{c}\hfill \delta (\mathbf{r})\equiv \frac{\rho (\mathbf{r})-\overline{\rho }}{\overline{\rho }},\end{array}$$(1)

where ρ(r) is the mass density at position r = [r_x, r_y, r_z], and $\overline{\rho }$ represents the average mass density of the Universe. The line-of-sight momentum field is defined by

$$\begin{array}{c}\hfill p(\mathbf{r})\equiv [1+\delta (\mathbf{r})]v(\mathbf{r}),\end{array}$$(2)

with v(r) denoting the line-of-sight peculiar velocity at position r. The auto-density power spectrum, P^δ, the auto-momentum power spectrum, P^p, and the cross-power spectrum of the density and momentum fields, P^δp, were originally formulated as two-point correlations of δ and/or p in Fourier space, as expressed in

$$\begin{array}{c}\hfill {(2\pi )}^3{\delta }^D(\mathbf{k}-\mathbf{k}\prime )P^{\delta }(\mathbf{k})\equiv ⟨\delta (\mathbf{k}){\delta }^{\ast }(\mathbf{k}\prime )⟩\end{array}$$(3)

$$\begin{array}{c}\hfill {(2\pi )}^3{\delta }^D(\mathbf{k}-\mathbf{k}\prime )P^p(\mathbf{k})\equiv ⟨p(\mathbf{k})p^{\ast }(\mathbf{k}\prime )⟩\end{array}$$(4)

$$\begin{array}{c}\hfill {(2\pi )}^3{\delta }^D(\mathbf{k}-\mathbf{k}\prime )P^{\delta p}(\mathbf{k})\equiv ⟨\delta (\mathbf{k})p^{\ast }(\mathbf{k}\prime )⟩,\end{array}$$(5)

respectively, where δ^D(k − k′) indicates the Dirac δ function and ‘*’ denotes the complex conjugate. k is the wave vector.

We favour the momentum power spectrum over the velocity power spectrum because, in principle, the velocity field is a continuous entity defined everywhere in space, even in regions devoid of galaxies. However, our measurements are limited to locations where galaxies (and thus mass) are present, effectively yielding a mass-weighted or momentum-based power spectrum. In regions devoid of galaxies, we record zero velocity not because the velocity is truly zero, but due to the absence of tracers. Additionally, the momentum power spectrum retains the same fundamental information as the velocity power spectrum, and with additional non-linear contributions introduced by the (1 + δ) term that are incorporated into the theoretical modelling.

Following the methodology outlined in Feldman et al. (1994), which constitutes the standard approach widely adopted today, the estimator for the density field is presented in

$$\begin{array}{c}\hfill F^{\delta }(\mathbf{r})\equiv \frac{w_{\delta }(\mathbf{r})[n_{\delta }(\mathbf{r})-\alpha n_s(\mathbf{r})]}{A_{\delta }}.\end{array}$$(6)

In this formulation w_δ(r), n_δ(r) and n_s(r) represent the weight factor, the number density of galaxies, and the number density of random points at position r. The number of random points is α times greater than the number of galaxies. A_δ serves as the normalization factor, defined by

$$\begin{array}{c}\hfill A_{\delta }^2={\displaystyle \int w_{\delta }^2}(\mathbf{r}){\overline{n}}_{\delta }^2(\mathbf{r})d^{3}r,\end{array}$$(7)

which ensures that the amplitude of the measured power aligns with it in a Universe unaffected by survey selection effects. The computation of the mean galaxy number density, ${\overline{n}}_{\delta }(\mathbf{r})$, at position r is detailed in Bautista et al. (2025) or see Fig. 3. As is indicated in Eq. (6), to estimate the density field it is necessary to subtract the random catalogue. This is because we measure the density contrast field, δ, which can only be accomplished if we possess knowledge of, and are able to subtract, the mean number of galaxies within each grid cell, estimated from randoms of that cell.

Similarly, following the derivation in Paper I, the estimator for the momentum field is given by

$$\begin{array}{c}\hfill F^p(\mathbf{r})\equiv \frac{w_p(\mathbf{r})n_p(\mathbf{r})v(\mathbf{r})}{A_p},\end{array}$$(8)

where w_p(r) and n_p(r) denote the weight factor and the number density of galaxies at position r, respectively. A_p functions as the normalization factor, defined by

$$\begin{array}{c}\hfill A_p^2={\displaystyle \int w_p^2}(\mathbf{r}){\overline{n}}_p^2(\mathbf{r})d^{3}r,\end{array}$$(9)

which again adjusts the amplitude of the measured power to conform with it in an idealized Universe without survey selection biases. The computation of the mean galaxy number density, ${\overline{n}}_p(\mathbf{r})$, at position r is also detailed in Bautista et al. (2025) (or see Fig. 3). Notably, no random catalogue is subtracted when estimating the momentum field. This is because that the mean velocity at any spatial point is inherently zero; consequently, subtracting a random catalogue during the construction of the momentum field would be equivalent to subtracting zero.

In this paper, we define the Fourier transform of the field functions presented in Eqs. (6) and (8) through the following formulation:

$$\begin{array}{c}\hfill F(\mathbf{k})\equiv \frac{1}{V}{\displaystyle \int F}(\mathbf{r})e^{i\mathbf{k}·\mathbf{r}}{d}^{3}r,\end{array}$$(10)

and inversely,

$$\begin{array}{c}\hfill F(\mathbf{r})\equiv \frac{1}{{(2\pi )}^3}{\displaystyle \int F}(\mathbf{k})e^{-i\mathbf{k}·\mathbf{r}}{d}^{3}k.\end{array}$$(11)

3.2. The estimators of the power spectrum multipoles

The observed redshift of a galaxy comprises two dominant contributions: the Hubble recessional redshift resulting from the cosmic expansion, and the peculiar velocity redshift arising from localized gravitational influences. Consequently, the inferred position of a galaxy based on its observed redshift does not represent its true comoving position; instead, it reveals a distorted position that is referred to as the ‘redshift-space position’. This phenomenon is widely recognized as redshift-space distortion (RSD).

Redshift-space distortions disrupt the spherical symmetry of the power spectrum with respect to the line of sight. To capture this anisotropic behavior, a prevalent methodology involves decomposing the redshift-space power spectrum, P(k), into a series of Legendre polynomials, L_ℓ(μ), i.e.

$$\begin{array}{c}\hfill P(\mathbf{k})={\displaystyle \sum _{\ell }}{P}_{\ell }(k)L_{\ell }(\mu ),(\ell =0,1,2,3,4,...).\end{array}$$(12)

In this decomposition, the angular dependence of the power spectrum, P(k), is encapsulated within the Legendre polynomials, L_ℓ(μ), while the amplitude information of P(k) are encoded in the multipole moments, P_ℓ(k), commonly referred to as power spectrum multipoles. Here, $\mu =\widehat{\mathbf{r}}·\widehat{\mathbf{k}}=\cos \theta$ denotes the cosine of the angle between the unit wave vector, $\widehat{\mathbf{k}}$, and the line-of-sight unit vector, $\widehat{\mathbf{r}}$.

Following the formalism established by Yamamoto et al. (2006) and Paper I, the estimators for the auto-density and auto-momentum power spectrum multipoles are defined by

$$\begin{array}{c}\hfill P_{\ell }^{\delta }(k)\equiv |F^{\delta }(k)F_{\ell }^{\delta \ast }(k)|-{\mathcal{N}}_{\ell }^{\delta }\end{array}$$(13)

and

$$\begin{array}{c}\hfill P_{\ell }^p(k)\equiv |F^p(k)F_{\ell }^{p\ast }(k)|-{\mathcal{N}}_{\ell }^p,\end{array}$$(14)

respectively. The indices ‘δ’ and ‘p’ corresponds to the auto-density and auto-momentum power spectrum, respectively. The associated shot-noise contributions are described by

$$\begin{array}{c}\hfill {\mathcal{N}}_{\ell }^{\delta }=\frac{(1+\alpha )(2\ell +1)}{A_{\delta }^{2}V}{\displaystyle \int w_{\delta }^2}(\mathbf{r}){\overline{n}}_{\delta }(\mathbf{r})L_{\ell }(\mu )d^{3}r\end{array}$$(15)

$$\begin{array}{c}\hfill {\mathcal{N}}_{\ell }^p=\frac{2\ell +1}{A_p^{2}V}{\displaystyle \int w_p^2}(\mathbf{r}){\overline{n}}_p(\mathbf{r})⟨v^2(\mathbf{r})⟩L_{\ell }(\mu )d^{3}r,\end{array}$$(16)

respectively.

As is outlined in Appendix H.1 (or see Paper III), the estimator for the cross-power spectrum multipoles is defined through³

$$\begin{array}{c}\hfill P_{\ell }^{\delta p}(k)\equiv \frac{1}{2}\mathrm{Im}\{|F^p(k)F_{\ell }^{\delta \ast }(k)|-|F^{\delta }(k)F_{\ell }^{p\ast }(k)|\}-{\mathcal{N}}_{\ell }^{\delta p},\end{array}$$(17)

with the corresponding shot-noise term detailed in

$$\begin{array}{c}\hfill {\mathcal{N}}_{\ell }^{\delta p}=\frac{2\ell +1}{A_{\delta }{A}_{p}V}{\displaystyle \int w_{\delta }}(\mathbf{r})w_p(\mathbf{r})min\{{\overline{n}}_{\delta }(\mathbf{r}),{\overline{n}}_p(\mathbf{r})\}⟨v(\mathbf{r})⟩L_{\ell }(\mu )d^{3}r.\end{array}$$(18)

In this context, ‘δp’ specifically refers to the density-momentum cross-power spectrum, while ‘Im’ denotes the imaginary part of a complex number. Given that both the time evolution of gravitational interactions and the initial conditions of the Universe remain invariant under the transformation δ → −δ and v → −v, the cross-power spectrum between the density and momentum fields must also exhibit invariance under such transformations. This leads directly to the identity $P_{\ell }^{\delta p}=-P_{\ell }^{\delta p\ast }$, which holds exclusively for purely imaginary functions. See Appendix B for more details.

The shot noise observed in the auto-power spectrum originates from the discrete nature of galaxies as tracers of the underlying density and momentum fields. In the case of the cross-power spectrum, shot noise arises due to the potential overlap of galaxies used in the estimation of both the density and momentum fields. This form of noise is proportional to the sample containing fewer objects (Smith 2009).

The expressions for the estimators of auto-power spectrum multipoles have been thoroughly examined by Bianchi et al. (2015). Nevertheless, our formulations of the field function definitions and Fourier transform definition deviate slightly from that of Bianchi et al. (2015). Consequently, below we restate the corresponding expressions for clarity and consistency. For convenience, we adopt the methodological framework proposed by Bianchi et al. (2015) and introduce the following function:

$$\begin{array}{c}\hfill T_{\ell }(\mathbf{k})\equiv {\displaystyle \int {(\widehat{\mathbf{k}}·\widehat{\mathbf{r}})}^{\ell }}F(\mathbf{r})e^{i\mathbf{k}·\mathbf{r}}{d}^{3}r,\end{array}$$(19)

where ℓ = 0, 1, 2, 3, 4, .... Consequently, based on Eq. (13) or Eq. (14), the corresponding even-multipoles of the auto-power spectrum are expressed in (see Appendix H.3 for more discussion)

$$\begin{array}{c}\hfill P_0(k)=\frac{1}{V}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }}[F(\mathbf{k})T_0^{\ast }(\mathbf{k})]-{\mathcal{N}}_0\end{array}$$(20)

$$\begin{array}{c}\hfill P_2(k)=\frac{5}{2V}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }}F(\mathbf{k})[3T_2^{\ast }(\mathbf{k})-T_0^{\ast }(\mathbf{k})]\end{array}$$(21)

$$\begin{array}{c}\hfill P_4(k)=\frac{9}{8V}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }}F(\mathbf{k})[35T_4^{\ast }(\mathbf{k})-30T_2^{\ast }(\mathbf{k})+3T_0^{\ast }(\mathbf{k})],\end{array}$$(22)

where dΩ_k denotes the solid angle. See Appendix D for the expressions of T₀, T₂, and T₄. It is important to note that, unlike in Bianchi et al. (2015), the normalization factors A_δ and A_p are inherently incorporated within the field functions in our formulation, and our definition of Fourier transformation diverges from them too. For multipoles with ℓ > 0, the shot noise becomes negligible and can be safely ignored (Blake et al. 2010; Howlett 2019). If ignoring the wide-angle or general relativistic effects, only the even-multipoles of the auto-power spectrum are non-zero (Bianchi et al. 2015; Beutler & Di Dio 2020; Castorina & White 2020).

In this paper, building upon Eq. (17) and Eq. (19), we derive the associated odd-multipoles for the cross-power spectrum, as expressed in

$$\begin{array}{c}\hfill P_1(k)=\frac{3}{V}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }}\frac{\mathrm{Im}\{F^p(\mathbf{k})T_1^{\delta \ast }(\mathbf{k})-F^{\delta }(\mathbf{k})T_1^{p\ast }(\mathbf{k})\}}{2},\end{array}$$(23)

$$\begin{array}{cc}\hfill P_3(k)& =\frac{7}{2V}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }}[5\times \frac{\mathrm{Im}\{F^p(\mathbf{k})T_3^{\delta \ast }(\mathbf{k})-F^{\delta }(\mathbf{k})T_3^{p\ast }(\mathbf{k})\}}{2}\hfill \\ \hfill & -3\times \frac{\mathrm{Im}\{F^p(\mathbf{k})T_1^{\delta \ast }(\mathbf{k})-F^{\delta }(\mathbf{k})T_1^{p\ast }(\mathbf{k})\}}{2}].\hfill \end{array}$$(24)

See Appendix D for the expressions of T₁ and T₃. For similar reasons, only the odd-multipoles exhibit non-vanishing contributions to the cross-power spectrum, and the influence of shot noise diminishes to a negligible extent for ℓ > 0.

3.3. The FKP weights

In Feldman et al. (1994), the authors introduced a weighting scheme for the density power spectrum aimed at minimizing the fractional variance of its estimation. This method is widely known as the Feldman-Kaiser-Peacock (FKP) weighting scheme and is given by

$$\begin{array}{c}\hfill w_{\delta }(\mathbf{r},k)=\frac{1}{1+{\overline{n}}_{\delta }(\mathbf{r})P_{\mathit{FKP}}^{\delta }(k)}.\end{array}$$(25)

Building upon this work, Paper I proposed an analogous weighting strategy tailored specifically for the momentum field, given by

$$\begin{array}{c}\hfill w_p(\mathbf{r},k)=\frac{1}{⟨v^2(\mathbf{r})⟩+{\overline{n}}_p(\mathbf{r})P_{\mathit{FKP}}^p(k)},\end{array}$$(26)

where ⟨v²(r)⟩ accounts for both intrinsic scatter and measurement errors associated with peculiar velocity estimates and can be calculated using Eq. (A.3).

The FKP weights were derived under the assumption that the power spectrum follows a Gaussian distribution. However, FKP weights may not represent the optimal choice when the power spectrum is not Gaussian. As is demonstrated in Sect. 6 of Paper I, varying the value of P_FKP(k) alters the effective survey depth of the galaxy survey. Different choices of FKP weights influence the trade-off between upweighting higher redshift data – thereby increasing cosmic volume and reducing power spectrum variance – and upweighting nearby objects where velocity measurement errors are smaller. Following Paper I, Paper II, and Paper III, we set $P_{\mathrm{FKP}}^{\delta }=1600h^{-3}{\mathrm{Mpc}}^3$ and $P_{\mathrm{FKP}}^p=5\times {10}^9{h}^{-3}{\mathrm{Mpc}}^3{\mathrm{km}}^2{\mathrm{s}}^{-2}$ for the density and momentum fields, respectively, to achieve optimal measurement accuracy.

3.4. Grid correction and Nyquist frequency

To compute the field functions of Eq. (6) and (8) from galaxies, it is necessary to assign the galaxies to a three-dimensional grid. This gridding process effectively convolves the field values with a top-hat window function in each spatial direction, resulting in each Fourier mode being multiplied by a Sinc function. To mitigate this effect, following the methodology outlined in Abate et al. (2008), Johnson et al. (2014), and Adams & Blake (2017), we applied a correction factor by dividing each voxel corresponding to a wave vector of k = [k_x, k_y, k_z] by the following Sinc functions for the Fourier transformed field functions (i.e. F(k) and T_ℓ(k)):

$$\begin{array}{c}\hfill \mathrm{\Gamma }(\mathbf{k})\equiv \frac{8}{L_x{L}_y{L}_z}\frac{sin(\frac{k_x{L}_x}{2})}{k_x}\frac{sin(\frac{k_y{L}_y}{2})}{k_y}\frac{sin(\frac{k_z{L}_z}{2})}{k_z},\end{array}$$(27)

where L_x, L_y, and L_z denote the size of the voxel in the x, y, and z directions, respectively. In this study, we constructed the galaxy grid using a voxel size of L_x = L_y = L_z = 2 Mpc h⁻¹ within a cubic box of size L = 800 Mpc h⁻¹ centred on the observer. This configuration provides sufficient spatial coverage to include the galaxies observed in the DESI survey.

Furthermore, we consider only those Fourier modes corresponding to wave numbers, k, that exceed the Nyquist frequency, defined by

$$\begin{array}{c}\hfill K_{\mathit{Ny}}=min\{\frac{\pi }{L_x},\frac{\pi }{L_y},\frac{\pi }{L_z}\}.\end{array}$$(28)

This ensures that the sampling rate meets the minimum requirement for an undistorted estimation of the power spectrum.

4. The theoretical model of the power spectrum

4.1. The models of power spectrum

Kaiser (1987) developed a linear model for the density power spectrum, commonly referred to as the Kaiser formula (see Appendix C), which has been extensively utilized over the past few decades. Building upon this foundation, Koda et al. (2014) extended the analysis to the power spectrum of peculiar velocities, and incorporated empirical damping terms into the linear models in order to more accurately account for non-linear motions of galaxies. However, these linear and quasi-linear models are effective only for k < 0.15 h Mpc⁻¹, where measurement errors tend to be larger.

To more effectively capture non-linear effects and extend the analysis into regions with smaller measurement errors, we employed one-loop perturbation theory models for the power spectrum in this paper. These redshift-space power spectrum models were initially developed by Vlah et al. (2012, 2013), Okumura et al. (2014) and summarized in Appendix A of Paper I and updated in Sect. 4 of Paper III. Following the theoretical framework outlined in Vlah et al. (2012, 2013) and Paper I, the density power spectrum model is described by

$$\begin{array}{cc}\hfill P^{\delta }(k,\mu )& =P_{00}+{\mu }^2(2P_{01}+P_{02}+P_{11})\hfill \\ \hfill & +{\mu }^4(P_{03}+P_{04}+P_{12}+P_{13}+\frac{1}{4}{P}_{22}).\hfill \end{array}$$(29)

Based on the theoretical formulation presented in Okumura et al. (2014) and Paper I, the momentum power spectrum model is given by

$$\begin{array}{c}\hfill P^p(k,\mu )=\frac{{[a(z)H(z)]}^2}{k^2}[P_{11}+{\mu }^2(2P_{12}+3P_{13}+P_{22})].\end{array}$$(30)

Additionally, the cross-power spectrum model was also derived in Okumura et al. (2014). However, according to Chen et al. (2025, Eq. (D.2) therein), there exist minor inaccuracies in the derivation of the cross-power spectrum expression in Okumura et al. (2014, Eq. (2.19) therein), and consequently in Paper III. In this paper, we derived the corrected version, which is given in

$$\begin{array}{cc}\hfill P^{\delta p}(k,\mu )& =\frac{a(z)H(z)}{k}\mu [P_{01}+P_{02}+P_{11}\hfill \\ \hfill & +{\mu }^2(\frac{3}{2}{P}_{03}+2P_{04}+\frac{3}{2}{P}_{12}+2P_{13}+\frac{1}{2}{P}_{22})].\hfill \end{array}$$(31)

Notably, this is the imaginary part of the cross-power spectrum model. The new formulations for P_mn, (m, n = 0, 1, 2, 3, 4) are presented in Sect. 4.2. In the aforementioned equations, the scale factor was determined from

$$\begin{array}{c}\hfill a(z)\equiv \frac{1}{1+z}.\end{array}$$(32)

Assuming ΛCDM, the Hubble parameter is given by

$$\begin{array}{c}\hfill \frac{H(z)}{H_0}=E(z)\equiv \sqrt{\frac{{\mathrm{\Omega }}_m}{a^3}+\frac{1-{\mathrm{\Omega }}_m-{\mathrm{\Omega }}_{\mathrm{\Lambda }}}{a^2}+{\mathrm{\Omega }}_{\mathrm{\Lambda }}},\end{array}$$(33)

where H₀, Ω_m, and Ω_Λ are the Hubble constant, matter density parameter, and dark energy density parameter in the present-day Universe at effective redshift z_eff = 0.

The models described in Eq. (29), (30), and (31) represent highly general frameworks, with the sole underlying assumption being the local plane-parallel approximation. These models do not presuppose the ΛCDM paradigm. To derive the ΛCDM-specific models, one simply needs to employ ΛCDM theory in the computation of E(z) (in Eq. (33)), and the growth factor D(a) (in Eq. (44)), as well as the linear matter power spectrum, P_L (as was mentioned in Sect. 4.2). However, one is entirely at liberty to utilize alternative theoretical approaches to calculate P_L, D(a) and E(z), thereby generating power spectrum models corresponding to other cosmological models.

To obtain the theoretical models for the power spectrum multipoles, one can substitute Eq. (29), (30) and (31) into the following:

$$\begin{array}{c}\hfill P_l(k)=\frac{(2\ell +1)}{2}{\displaystyle {\int }_{-1}^1}P(k,\mu )L_l(\mu )d\mu ,\end{array}$$(34)

which essentially performs a Legendre transformation corresponding to Eq. (12). The aforementioned theoretical models for the power spectrum can be readily converted into models for the galaxy two-point correlation function, the velocity correlation function, and the galaxy-velocity cross-correlation function using the HANKL package⁴ (Karamanis & Beutler 2021); see Turner et al. (2025) and Appendix C for more discussion. The momentum correlation is equivalent to the velocity correlation in linear scales, as demonstrated by Eq. (2) of Paper I; therefore, the model of the velocity correlation function can be converted from the model of the momentum power spectrum approximately.

4.2. The new loop terms for the cross power spectrum

The loop terms, P_mn, (m, n = 0, 1, 2, 3, 4), in Eq. (29), 30, and 31 were initially formulated by Vlah et al. (2012, 2013), Okumura et al. (2014), Paper I and Paper III. Their models, however, were specifically tailored for situations in which both density and momentum fields originate from the same galaxy survey dataset.

In this paper, we generalize P_mn of cross-power spectrum Eq. (31) to accommodate cases in which the density and momentum fields stem from different surveys with differing observational properties. This generalization implies that the cross-power spectrum Eq. (31) will incorporate multiple biasing parameters. The corresponding new formulations for P_mn (of Eq. (31)) are detailed in the subsequent equations:

$$\begin{array}{cc}\hfill P_{00}& =b_1^{\delta }{b}_1^p{D}^2{P}_L+b_1^{\delta }{D}^4(b_2^p{K}_{00}+b_s^p{K}_{00}^s+2b_1^p(I_{00}+3J_{00}{k}^2{P}_L))\hfill \\ \hfill & +\frac{1}{2}{D}^4(b_2^{\delta }{b}_2^p{K}_{01}+2b_1^p(b_2^{\delta }{K}_{00}+b_s^{\delta }{K}_{00}^s)+b_s^{\delta }{b}_s^p{K}_{01}^s\hfill \\ \hfill & +b_2^p{b}_s^{\delta }{K}_{02}^s+b_2^{\delta }{b}_s^p{K}_{02}^s+4b_{3nl}^p{P}_L{\sigma }_3^2),\hfill \end{array}$$(35)

$$\begin{array}{cc}\hfill P_{01}& =f{D}^2(b_1^{\delta }(P_L+2D^2(I_{01}+b_1^p{I}_{10}+3(J_{01}+b_1^p{J}_{10})k^2{P}_L))\hfill \\ \hfill & -D^2(b_2^{\delta }(K_{10}+b_1^p{K}_{11})+b_s^{\delta }{K}_{10}^s+b_1^p{b}_s^{\delta }{K}_{11}^s+b_{3nl}^p{P}_L{\sigma }_3^2)),\hfill \end{array}$$(36)

$$\begin{array}{cc}\hfill P_{11}& =f^2{D}^2({\mu }^2{P}_L+D^2(b_1^{\delta }{b}_1^p{I}_{31}+(2I_{11}+b_1^{\delta }{b}_1^p{I}_{13}+2b_1^{\delta }{I}_{22}\hfill \\ \hfill & +2b_1^p{I}_{22}+6((b_1^{\delta }+b_1^p)J_{10}+J_{11})k^2{P}_L){\mu }^2)),\hfill \end{array}$$(37)

$$\begin{array}{cc}\hfill P_{13}& =-f^4{k}^2{D}^2({\sigma }_{\delta ,vS}^2{\mu }^2(P_L+D^2(2I_{11}+2(b_1^{\delta }+b_1^p)I_{22}\hfill \\ \hfill & +6k^2{P}_L(J_{11}+(b_1^{\delta }+b_1^p)J_{10})))+{\sigma }_{\delta ,vT}^2{b}_1^{\delta }{b}_1^p{D}^2({\mu }^2{I}_{13}+I_{31})),\hfill \end{array}$$(38)

$$\begin{array}{cc}\hfill P_{04}& =-\frac{1}{2}{f}^4{b}_1^{\delta }{k}^2{\sigma }_{\delta ,vT}^2{D}^4(I_{02}+{\mu }^2{I}_{20}+2k^2{P}_L(J_{02}+{\mu }^2{J}_{20}))\hfill \\ \hfill & +\frac{1}{4}{f}^4{k}^4{b}_1^{\delta }{b}_1^p{P}_{00}({\sigma }_{\delta ,vT}^4+{\sigma }_4^2),\hfill \end{array}$$(39)

where the indices δ and p denote the parameters for density and momentum fields, respectively. The terms P₀₀ in the expression of P₀₄ should be calculated from Eq. (35) of this paper. Following the derivation outlined in Paper III, we arrive at

$$\begin{array}{cc}\hfill P_{02}& =f^2{b}_1^{\delta }{D}^4(I_{02}+{\mu }^2{I}_{20}+2k^2{P}_L(J_{02}+{\mu }^2{J}_{20}))-f^2{k}^2{\sigma }_{\delta ,vT}^2{P}_{00}\hfill \\ \hfill & +f^2{D}^4(b_2^{\delta }(K_{20}+{\mu }^2{K}_{30})+b_s^{\delta }(K_{20}^s+{\mu }^2{K}_{30}^s)),\hfill \end{array}$$(40)

$$\begin{array}{cc}\hfill P_{12}& =f^3{D}^4(I_{12}+{\mu }^2{I}_{21}-b_1^{\delta }(I_{03}+{\mu }^2{I}_{30})+2k^2{P}_L(J_{02}+{\mu }^2{J}_{20}))\hfill \\ \hfill & -f^2{k}^2{\sigma }_{\delta ,vT}^2{P}_{01}+2f^3{D}^4{k}^2{\sigma }_{\delta ,vT}^2(I_{01}+I_{10}+3k^2{P}_L(J_{01}+J_{10})),\hfill \end{array}$$(41)

where P₀₀ and P₀₁ in the above equations should be obtained from Eqs. (35) and (36) of this paper, separately. By following the derivation presented in Paper I, we also obtain

$$\begin{array}{c}\hfill P_{03}=-f^2{k}^2{\sigma }_{\delta ,vS}^2{P}_{01},\end{array}$$(42)

$$\begin{array}{cc}\hfill P_{22}& =\frac{1}{4}{f}^4{D}^4(I_{23}+2{\mu }^2{I}_{32}+{\mu }^4{I}_{33})+f^4{k}^4{\sigma }_{\delta ,vT}^4{P}_{00}-f^2{k}^2{\sigma }_{\delta ,vT}^2\hfill \\ \hfill & \times (2P_{02}-f^2{D}^4(b_2^{\delta }(K_{20}+{\mu }^2{K}_{30})+b_s^{\delta }(K_{20}^s+{\mu }^2{K}_{30}^s))).\hfill \end{array}$$(43)

In the above, the P₀₀, P₀₁, and P₀₂ should be derived from Eq. (35), (36), and (40) of this paper, respectively. The term D represents the linear growth factor. Under the assumption of general relativity (GR) and ΛCDM, its analytical expression is given by (Heath 1977; Howlett et al. 2015)

$$\begin{array}{c}\hfill D(a)=\frac{D_{\mathit{gr}}(a)}{D_{gr,0}}\end{array}$$(44)

where

$$\begin{array}{c}\hfill D_{\mathit{gr}}(a)=E(a){\displaystyle {\int }_0^a}\frac{da\prime }{a\prime {}^{3}H{(a\prime )}^3},\end{array}$$(45)

and where D_gr, 0 symbolizes the D_gr(a) of the present-day Universe at redshift z_eff = 0. E(a) is given in Eq. (33). The terms I_mn, K_mn, K_mn^s, σ₃, and σ₄ were derived from integrations of the linear matter power spectrum P_L, following the approach detailed in Vlah et al. (2012) and Paper I. In this study, we utilized the CAMB⁵ package (Lewis et al. 2000) to calculate P_L(k) assuming ΛCDM, with the resulting P_L(k) illustrated in Fig. C.1.

The parameters of the power spectrum models are summarized in Table 1. Within the context of the cross-power spectrum formalism Eq. (31), the fully parameterized model comprises ten parameters. On the other hand, by equating the biasing parameters and velocity dispersion parameters of the density field with those of the momentum field, Eqs. (35) to (43) simplify to the expressions detailed in Appendix A of Paper I. As a result, we obtain the P_mn for the auto-density and auto-momentum power spectrum models (i.e. Eqs. (29) and (30)) conform to the formulations proposed by Paper I, Paper II and Paper III. As was argued in McDonald & Roy (2009), Saito et al. (2014) and Paper I, two higher-order galaxy biasing parameters can be determined from the linear bias using local Lagrangian relations:

$$\begin{array}{c}\hfill b_s=-\frac{4}{7}(b_1-1),b_{3nl}=\frac{32}{315}(b_1-1).\end{array}$$(46)

They can also be set as free parameters when fitting the models to the measurements to gain more flexibility from the models.

Furthermore, setting all of the loop terms, I_mn, K_mn, K_mn^s, σ₃, and σ₄, to zero, in conjunction with setting the higher-order biasing and velocity dispersion parameters to zero, leads to a reduction in the power spectrum models in Eqs. (29), (30), and (31) to the Kaiser formulas, as is demonstrated in Fig. C.2.

4.3. The window function convolution

In actual astronomical observations, our view is confined to galaxies residing within a limited spatial volume and with a particular sky completeness function. Consequently, the power spectrum measured from galaxy surveys reflects only the characteristics of this restricted volume, rather than representing the entire Universe. However, theoretical formulations of the power spectrum are typically based on the infinitely extended Universe. Therefore, in order to reliably infer cosmological parameters from galaxy surveys, it is crucial to incorporate the impacts of survey geometry and observational completeness when comparing theoretical models with measured power spectrum. This necessitates the application of a convolution operation involving the window function, as formally described by

$$\begin{array}{c}\hfill {\mathbf{P}}^c(k)=\mathbf{W}·\mathbf{P}(k\prime ),\end{array}$$(47)

where the model power spectrum multipoles P_ℓ′(k′) (ℓ′ = 0, 1, 2, 3, 4, ...) are concatenated into a single vector, P(k′) = [P₀(k′), P₁(k′), P₂(k′), P₃(k′), P₄(k′), ...]. The term P^c(k) denotes the convolved model power spectrum multipoles, which serves as the counterpart for direct comparison with the measured power spectrum multipole. The window function convolution matrix, W, is characterized by dimensions (N_ℓ × N_k)×(N_ℓ′ × N_k′). Specifically, N_k denotes the number of k bins associated with each measured power spectrum multipole, while N_ℓ represents the total count of measured power spectrum multipoles. Furthermore, N_k′ indicates the number of k′ values employed in computing each model power spectrum multipole, and N_ℓ′ refers to the overall number of theoretical power spectrum multipoles considered.

The window function convolution matrix, W, was constructed using random catalogues, following the approach detailed in Appendix H.2. Each entry of the matrix, W, was defined according to

$$\begin{array}{c}\hfill {\mathbf{W}}_{\ell }(k,k\prime )=\frac{4\pi }{A^2}{\displaystyle \int d}{\mathrm{\Omega }}_k{\displaystyle \sum _{m=-\ell }^{\ell }}{Y}_m^{\ell }(\widehat{\mathbf{k}}){\displaystyle \sum _{\ell \prime }}{\displaystyle \sum _{m\prime =-\ell \prime }^{\ell \prime }}\frac{I_{mm\prime }^{\ell \ell \prime }(\widehat{\mathbf{k}},\widehat{\mathbf{k}}\prime )}{(2\ell \prime +1)V},\end{array}$$(48)

where

$$\begin{array}{c}\hfill I_{m{m}^{\prime }}^{\ell {\ell }^{\prime }}(\widehat{\mathbf{k}},{\widehat{\mathbf{k}}}^{\prime })=\frac{1}{{(2\pi )}^3}{\displaystyle \int G}(\mathbf{k}-{\mathbf{k}}^{\prime }){\stackrel{\sim }{S}}_{m{m}^{\prime }}^{\ell {\ell }^{\prime }\ast }(\mathbf{k}-{\mathbf{k}}^{\prime })Y_{m^{\prime }}^{{\ell }^{\prime }\ast }({\widehat{\mathbf{k}}}^{\prime })d^3{k}^{\prime },\end{array}$$(49)

$$\begin{array}{c}\hfill {\stackrel{\sim }{S}}_{m{m}^{\prime }}^{\ell {\ell }^{\prime }}(\mathbf{k}-{\mathbf{k}}^{\prime })=\frac{1}{V}{\displaystyle \int Y_m^{\ell }}(\widehat{\mathbf{r}})Y_{m^{\prime }}^{{\ell }^{\prime }\ast }(\widehat{\mathbf{r}})\overline{n}(\mathbf{r})w(\mathbf{r})e^{i\mathbf{r}·(\mathbf{k}-{\mathbf{k}}^{\prime })}{d}^{3}r,\end{array}$$(50)

$$\begin{array}{c}\hfill G(\mathbf{k}-{\mathbf{k}}^{\prime })={\displaystyle \int w}(\mathbf{r})\overline{n}(\mathbf{r})e^{i\mathbf{r}·(\mathbf{k}-{\mathbf{k}}^{\prime })}{d}^{3}r,\end{array}$$(51)

where $Y_m^{\ell }(\widehat{\mathbf{r}})$ corresponds to the spherical harmonics, w(r) symbolizes the FKP weight introduced in Sect. 3.3, and $\overline{n}$ stands for the mean galaxy number density outlined in Sects. 2 and 3.1.

For the auto-power spectrum, w(r) and $\overline{n}(\mathbf{r})$ should be computed individually for both the density and momentum fields, and the A² should be computed from Eqs. (7) or (9), respectively. However, when considering the cross-power spectrum, and in accordance with the arguments articulated in Appendix H.2, the following substitutions must be made to the normalization factor in Eq. (48):

$$\begin{array}{c}\hfill A^2\to A_p{A}_{\delta },\end{array}$$(52)

along with the following corresponding modifications to Eq. (49),

$$\begin{array}{cc}& G(\widehat{\mathbf{k}}-{\widehat{\mathbf{k}}}^{\prime }){\stackrel{\sim }{S}}_{m{m}^{\prime }}^{\ell {\ell }^{\prime }\ast }(\widehat{\mathbf{k}}-{\widehat{\mathbf{k}}}^{\prime })\to \hfill \\ \hfill & \frac{1}{2}[G^p(\widehat{\mathbf{k}}-\widehat{\mathbf{k}}\prime ){\stackrel{\sim }{S}}_{m{m}^{\prime }}^{\delta ,\ell {\ell }^{\prime }\ast }(\widehat{\mathbf{k}}-{\widehat{\mathbf{k}}}^{\prime })+G^{\delta }(\widehat{\mathbf{k}}-{\widehat{\mathbf{k}}}^{\prime }){\stackrel{\sim }{S}}_{m{m}^{\prime }}^{p,\ell {\ell }^{\prime }\ast }(\widehat{\mathbf{k}}-\widehat{\mathbf{k}}\prime ),]\hfill \end{array}$$(53)

where G^δ and G^p are derived from Eq. (51), applied separately to the density and momentum fields, whereas the terms ${\stackrel{\sim }{S}}_{mm\prime }^{\delta ,\ell \ell \prime *}$ and ${\stackrel{\sim }{S}}_{mm\prime }^{p,\ell \ell \prime }$ are obtained through Eq. (50), also evaluated independently for the density and momentum fields.

The window function convolution matrix was computed by assigning the random points to a three-dimensional grid. In this paper, we set n_x = n_y = n_z = 140 for a cubic box size of L = 700 Mpc h⁻¹, in order to get a grid without empty voxels in the region where random points are fulfilled.

Figure 4 illustrates the window function convolution matrix across different types of power spectrum, including the density power spectrum of the BGS data ( left panel), the momentum power spectrum of the DESI-PV data (middle panel), and their cross-power spectrum ( right panel). Notably, as is revealed in Fig. 4, the off-diagonal blocks of the W exhibit non-zero values, signifying that the convolution of any given multipole receives contributions from other multipoles. This interdependence implies that higher-order multipoles must be incorporated into the modelling framework, even if they are not directly involved in the final fitting procedure. In this work, we include all multipoles up to ℓ_max = 4 to ensure comprehensive and accurate representation. As is shown in Fig. 5, the differences arising from varying ℓ_max values are negligible, indicating that ℓ = 4 is sufficiently large to ensure convergence. We employed k′∈[0, 0.4075] h Mpc⁻¹ (with 300 bins) in computing each model’s power spectrum multipole, which should be slightly larger than the k ∈ [0, 0.3] h Mpc⁻¹ that is associated with each measured power spectrum multipole, in order to fully cover the diagonal information of each block of the convolution matrix.

5. Tests on mocks

5.1. Fitting method

In this study, we concatenated the density, momentum, and cross power spectrum to effectively constrain the growth rate fσ₈. By concatenating the models of these three spectrum, i.e. Eqs. (29), (30), and (31), we derived a comprehensive framework comprising 13 parameters in total. These encompass the growth rate, fσ₈, the biasing parameters [b₁σ₈, b₂σ₈, b_sσ₈, b_3nlσ₈] for both the density and momentum fields, as well as the velocity dispersion parameters [σ_vS, σ_vT] for the density and momentum fields, respectively. However, in practical applications, it is not necessary to treat all 13 parameters as independent variables during model fitting with measurements, because these models are not equally sensitive to every parameter. Through the following testing using mocks, we discovered that a streamlined set of only five key parameters [ fσ₈,  b₁σ₈,  b₂σ₈,  σ_vT,  σ_vS ] suffices to maintain robustness and precision in the fitting process (see Appendix E for more discussion). This led us to adopt a unified linear biasing parameter, b₁σ₈, for both the density and momentum fields, and a shared second-order biasing parameter, b₂σ₈, for both fields. Furthermore, we employed shared velocity dispersion parameters [σ_vT,  σ_vS] across both the density and momentum fields.

In our study, the covariance matrix of the power spectrum is derived not from the actual Universe, but from a finite number of mocks. This limitation introduces an uncertainty into the inverse of the covariance matrix. This issue is referred to as the Hartlap effect (Hartlap et al. 2007). To address this challenge, we adopted the methodology proposed by Sellentin & Heavens (2016), which is different from the conventional approach of directly minimizing the χ². Instead, we employed a likelihood function based on the following modified t distribution to derive more robust estimates of the parameters θ = [fσ₈, b₁σ₈, b₂σ₈, σ_vT², σ_vS²]:

$$\begin{array}{c}\hfill \mathcal{L}(P|\theta )\propto {|\mathsf{C}|}^{-\frac{1}{2}}{[1+\frac{{\chi }^2(P|\theta )}{N-1}]}^{-\frac{N}{2}},\end{array}$$(54)

where N = 675 is the number of mocks. Within this framework, the χ² is formulated according to ${\chi }^2({\mathit{P}}_o|\theta )=[{\mathbf{P}}_o-{\mathbf{P}}_m^c(\theta )]{\mathsf{C}}^{-1}{[{\mathbf{P}}_o-{\mathbf{P}}_m^c(\theta )]}^T$, where P_o signifies the measured power spectrum, ${\mathbf{P}}_m^c$ denotes the model power spectrum convolved with the survey window function, and C represents the covariance matrix.

In this research, cosmological parameters were extracted through a combination of the auto-density power spectrum monopole, P₀^δ, and quadrupole, P₂^δ, the auto-momentum power spectrum monopole, P₀^p, and the cross-power spectrum dipole, P₁^δp. Specifically, both P_o = [P₀^δ, P₂^δ, P₀^p, P₁^δp] and P_m^c = [P₀^δc, P₂^δc, P₀^pc, P₁^δpc] are represented as 4N_k vectors, while the covariance matrix, C, spans a dimensionality of 4N_k × 4N_k. Additional measured multipole components of the power spectrum of the survey data are found to be close to zero and exhibit high noise levels; therefore, they are excluded from the parameter estimation process in this work.

Parameter inference was conducted using the Metropolis-Hastings Markov chain Monte Carlo (MCMC) methodology, under uniformly flat priors defined over the interval fσ₈ ∈ (0, 1.5], b₁σ₈ ∈ [0, 3], b₂σ₈ ∈ [ − 5, 5], σ_vT² ∈ (0, 350²], and σ_vS² ∈ (0, 350²]. Our growth rate constraints were obtained through the measurement of power spectrum in bins of k ∈ [0.025,  0.3], with bin widths of 0.009h Mpc⁻¹.

5.2. Fitting fσ₈ of the mocks

We began by fitting the power spectrum of the mock average. As is illustrated in Fig. 6, the filled circles represent the average of the power spectrum measured from 675 mocks, while the solid curves depict the best-fit theoretical models. The density power spectrum was measured from the BGS mocks, whereas the momentum power spectrum was obtained from the DESI-PV mocks. Each cross-power spectrum was computed using a BGS mock and its corresponding DESI-PV mock (which has the same observer as that BGS mock).

The MCMC resulting parameter estimates are summarized in Table 2 and illustrated in the right panel of Fig. 6. Notably, the estimated value of $f{\sigma }_8=0.{468}_{-0.043}^{+0.061}$ closely aligns with the mock fiducial value of fσ_8,fid = 0.466 at effective redshift z_eff = 0.2. This agreement demonstrates that our fitting methodology effectively recovers the true growth rate embedded in the mocks, and that the power spectrum models perform robustly up to the non-linear scale of k_max = 0.3 h Mpc⁻¹ under current conditions. In this analysis, we excluded the momentum power spectrum quadrupole, P₂^p, as it is heavily dominated by noise in survey data. To further assess the adaptability of our model, we also conducted fits incorporating P₂^p using the mock average, where P₂^p exhibits sufficient smoothness to yield a reliable fit. For a more detailed exploration, please refer to Appendix F.

We proceeded to fit fσ₈ of 15 randomly selected DESI-PV mocks and their corresponding BGS mock counterparts. As is illustrated in the top panel of Fig. 7, the filled squares represent the MCMC estimates of fσ₈ derived from each of the 15 mock sets, all of which align closely with the fiducial fσ₈ value (indicated by the dashed yellow line). The shaded blue band displays the fσ₈ constraint extracted from Table 2 (i.e. mock average), providing an independent reference for comparison. The reduced χ² values for all 15 measurements are found to be near unity, underscoring the robustness of the model and confirming that it provides an excellent statistical fit to the mock data. Figure 8 presents the growth rate constraints derived from fitting the power spectrum using different wave number cut-offs k_max, demonstrating that our models yield stable fits and reliable for k_max < 0.35 h Mpc⁻¹.

6. Results

6.1. The systematic error of the FP data

In Fig. 9, we compare the measured momentum power spectrum of the FP survey data with the average momentum power spectrum derived from 675 FP mocks. As the plot illustrates, the momentum power spectrum from the FP data begins to deviate. Specifically, it becomes elevated relative to the mocks starting at k = 0.05 h Mpc⁻¹. This discrepancy becomes even more pronounced at k = 0.1 h Mpc⁻¹, where the data exhibits significantly higher values than the mock predictions. This bias originates from a systematic bias in the measured η values of the FP data. In particular, within larger galaxy groups, the measured η tends to be more negative, leading to this anomalous behavior in the momentum power spectrum. This bias was initially identified by Howlett et al. (2022) through an analysis of the SDSSv data and has since been confirmed in DESI FP data, as reported by Ross et al. (2025). However, the underlying physical cause of this bias remains elusive and not yet fully understood. In future endeavours, a comprehensive galaxy group catalogue of DESI will be meticulously constructed, serving as a pivotal foundation for the calibration of this bias.

This systematic deviation also affects the estimation of fσ₈, as the DESI-PV catalogue is predominantly composed of FP galaxies rather than TF galaxies. As is shown in the top panel of Fig. 10, the filled blue squares represent the estimated values of the growth rate, fσ₈ (derived from the power spectrum of the BGS and DESI-PV data), as a function of the cut-off wave number, k_max^p, of the momentum power spectrum. We fixed k_max = 0.3 h Mpc⁻¹ for the density and cross power spectrum. Notably, at $k_{\text{max}}^p=0.125h{\text{Mpc}}^{-1}$, the measured fσ₈ abruptly rises from approximately 0.4 to around 0.6.

To mitigate this bias, we have to adopt a conservative approach by restricting the momentum power spectrum analysis to $k_{\text{max}}^p=0.1h{\text{Mpc}}^{-1}$ when fitting fσ₈ from survey data, while retaining k_max = 0.3 h Mpc⁻¹ for the density and cross power spectrum. The information encoded in the momentum power spectrum is predominantly concentrated on the larger cosmological scales, specifically at wave numbers below k ≤ 0.1 h Mpc⁻¹ (Howlett 2019; Bautista et al. 2025). Consequently, adopting a conservative cut-off at $k_{\text{max}}^p=0.1h{\text{Mpc}}^{-1}$ enables a robust analysis of the data without sacrificing significant signal content and does not impact our main conclusion.

It is important to emphasize that $k_{\text{max}}^p=0.1h{\text{Mpc}}^{-1}$ reflects the reliability limit of the measured momentum power spectrum, rather than a breakdown of the theoretical model. Indeed, the model of momentum power spectrum remains robust and accurate up to k_max = 0.3 h Mpc⁻¹ as demonstrated in Sect. 5.2. In future studies, with access to an unbiased velocity survey, we anticipate being able to extend the momentum power spectrum analysis to k_max = 0.3 h Mpc⁻¹, thereby enhancing the precision and constraining power of our cosmological inferences.

6.2. Result of fσ₈

As is depicted in Fig. 11, the power spectrum fitting results for the BGS and DESI-PV data are presented. In the left-hand panels, the top, middle, and bottom subplots display the measured density monopole and quadrupole, momentum monopole, and cross dipole power spectrum, respectively, with filled circles representing the measurements of the data. The curves illustrate the theoretical model spectrum that have been fitted to the measured data points. The goodness-of-fit is quantified by a reduced chi-squared value of χ²/d.o.f. = 116.563/(98 − 5) = 1.2534, indicating a satisfactory agreement between the model and the data. The derived cosmological parameter constraints are summarized in Table 3 and illustrated in the right panel of Fig. 11. Notably, the estimate for the growth rate is $f{\sigma }_8(z_{\mathrm{eff}}=0.07)=0.{440}_{-0.096}^{+0.080}$ at the effective redshift z_eff = 0.07.

To compute the theoretically predicted value of fσ₈, we followed the formalism outlined in Howlett et al. (2015), wherein the growth rate as a function of redshift is expressed as

$$\begin{array}{c}\hfill f(a){\sigma }_8(a)={\mathrm{\Omega }}_m{(a)}^{\gamma }{\sigma }_8\frac{D_{\mathit{gr}}(a\ast )}{D_{gr,0}}\frac{D_{\gamma }(a)}{D_{\gamma }(a\ast )},\end{array}$$(55)

with

$$\begin{array}{c}\hfill \frac{D_{\gamma }(a)}{D_{\gamma }(a\ast )}=\exp ({\displaystyle {\int }_{a\ast }^a}{\mathrm{\Omega }}_{\mathrm{m}}{(a\prime )}^{\gamma }d\ln a\prime ),{\mathrm{\Omega }}_{\mathrm{m}}(a)=\frac{{\mathrm{\Omega }}_{\mathrm{m}}}{a^{3}E{(a)}^2},\end{array}$$(56)

and D_gr(a) obtained from Eq. (45). Under the framework of GR, the growth index, γ, was fixed at 0.55 (Linder 2005). Consequently, with the Planck Collaboration VI (2020) cosmology, the GR predicts a growth rate of fσ₈ = 0.446 at z_eff = 0.07, which aligns remarkably well with our measured value of $f{\sigma }_8=0.{440}_{-0.096}^{+0.080}$ at the 68% confidence level, underscoring the consistency between theory and observation.

The refined estimate of the growth rate, derived by synthesizing the fitting outcomes from the correlation function (Turner et al. 2025), maximum likelihood estimation (MLE, Lai et al. 2025), and the power spectrum of this paper, yields a consensus value of $f{\sigma }_8(z_{\text{eff}}=0.07)=0.{450}_{-0.055}^{+0.055}$, in concordance with theoretical predictions too. For a deeper exploration of data synthesis, see Bautista et al. (2025).

6.3. Result of γ

We present constraints on Linder’s growth index, γ (see Eq. (55)), assuming a ΛCDM background, obtained with cobaya (Torrado & Lewis 2021) interfaced to MGCAMB (Wang et al. 2023; Hojjati et al. 2011) using MG_flag = 2. In this setting, γ is mapped into the modified scale and redshift-dependent Poisson factor μ(a, k). In our analysis, we used DESI DR1 ShapeFit measurements in six redshift bins, distance-scale information from the post-reconstruction correlation function, and the baryon acoustic oscillation (BAO) Lyα likelihood used in DR1 (DESI Collaboration 2025d)⁶. We sampled the background parameters in addition to γ and report σ₈ as a derived quantity. To mitigate projection effects, we applied a Big-Bang nucleosynthesis Gaussian prior, ω_b ∼ 𝒩(0.02218, 0.00055), and a broad Gaussian prior on n_s centred on the Planck fit value with standard deviation inflated by a factor of ten, π(n_s)∼𝒩(0.9649, 0.042) (Planck Collaboration VI 2020). The addition of PV fσ₈ measurements with DESI DR1 ShapeFit breaks the σ₈–γ degeneracy and sharpens the constraints on all parameters.

Figure 12 shows the fitting results for (γ, Ω_m, σ₈) obtained from DESI DR1 ShapeFit (SF) + BAO, combined with our PV consensus, $f{\sigma }_8=0.{450}_{-0.055}^{+0.055}$ (orange contours). We find $\gamma =0.{580}_{-0.110}^{+0.110}$, ${\mathrm{\Omega }}_{\mathrm{m}}=0.{301}_{-0.011}^{+0.011}$, and ${\sigma }_8=0.{834}_{-0.032}^{+0.032}$. Our result is consistent with GR. The PV fσ₈ information reduced the errors of γ, Ω_m, and σ₈ by ≈32%, ≈33%, and ≈20% respectively, relative to SF+BAO alone (see blue contours or Table 4). Our measurement of γ falls notably below the values reported by Nguyen et al. (2023) and Calabrese et al. (2025), who obtain higher estimate values of γ ∼ 0.65.

7. Discussion

Figure 13 displays the evolution of the growth rate, fσ₈, as a function of redshift, z. The blue curve was derived from Eq. (55), under the assumption of GR and a ΛCDM cosmological model based on Planck Collaboration VI (2020). The filled green circle represents the measurement derived from this paper. In comparison with the results presented in Turner et al. (2025) and Lai et al. (2025), our measurement exhibits slightly broader error margins, a consequence of adopting a lower wave number cut-off, $k_{\text{max}}^p=0.1h{\text{Mpc}}^{-1}$, in the analysis of the momentum power spectrum, relative to the higher thresholds employed in the other two studies. The filled red diamond denotes the consensus outcome of these three measurements. The other measurements are: H17: Howlett et al. (2017); J14: Johnson et al. (2014); A20: Adams & Blake (2020); Q19: Qin et al. (2019a); Q25: Qin et al. (2025); B12: Beutler et al. (2012); As23: Appleby et al. (2023); T23: Turner et al. (2023); D19: Dupuy et al. (2019); W18: Wang et al. (2018); C15: Carrick et al. (2015); S20: Said et al. (2020); Bp24: Boubel et al. (2024); Bs20: Boruah et al. (2020); DESI: DESI Collaboration (2025c); SDSS: Alam et al. (2021); B11 WiggleZ: Blake et al. (2011).

Although these measurements largely align with the predictions of GR, a subtle discrepancy becomes apparent when they are analysed in combination. To further investigate this, we integrated the aforementioned measured fσ₈ values (i.e. all the data points in Fig. 13; we only used the filled red diamond for DESI-PV) with the Planck chain⁷ and performed a fitting analysis for γ. The resulting fit value of γ = 0.63 ± 0.02 was then transformed into a corresponding range of fσ₈ values across different redshifts and is presented as the orange-coloured curve in Fig. 13. This result indicates a moderate tendency in the data towards a higher γ value, implying a potentially weaker gravitational model compared to the predictions of GR (blue curve in Fig. 13). However, it is important to note that the fit value and associated uncertainty for this γ should be interpreted with caution, as there exists considerable overlap – and thus covariance – among many of the measurements in Fig. 13, which has not been fully incorporated into the analysis. Nonetheless, as is shown in Fig. 13, a minor tension remains between current observational data and the theoretical predictions of GR, a discrepancy that may be further clarified or resolved through future research and surveys.

8. Conclusions

In this paper, we build upon the foundational research presented in Paper I, Paper II, and Paper III by conducting a comprehensive joint analysis of the monopole and quadrupole moments of the auto-density power spectrum, the monopole of the auto-momentum power spectrum, and the dipole component of the cross-power spectrum to measure the growth rate, fσ₈. The density field was derived from the BGS, whereas the momentum field was extracted from the DESI-PV. Our investigation yields the following key findings:

1. We have systematically presented the power spectrum estimators and the window function convolution matrix.

2. We have refined the theoretical model of the cross-power spectrum and extended its applicability to scenarios in which the density and momentum fields originate from different survey catalogues.

3. Utilizing mocks, we demonstrate that the power spectrum models exhibit robust performance up to the non-linear scale of k_max = 0.3 h Mpc⁻¹ when fitting fσ₈ using the set [P₀^δ, P₂^δ, P₀^p, P₁^δp] and treating [fσ₈, b₁σ₈, b₂σ₈, σ_vT², σ_vS²] as free parameters.

4. Furthermore, our analysis of mocks reveals that the models remain reliable up to the non-linear scale of k_max = 0.25 h Mpc⁻¹ when fitting fσ₈ using [P₀^δ, P₂^δ, P₀^p, P₂^p, P₁^δp] and allowing $[f{\sigma }_8,b_1^{\delta }{\sigma }_8,b_2{\sigma }_8,b_1^p{\sigma }_8,{\sigma }_{\delta ,vT}^2,{\sigma }_{\mathit{vS}}^2,{\sigma }_{p,vT}^2]$ to vary freely.

5. Our results indicate that the systematic bias of FP data enhances the measured momentum power spectrum at k > 0.1 h Mpc⁻¹. Consequently, we restrict our fitting range of the momentum power spectrum to k ≤ 0.1 h Mpc⁻¹ to mitigate this effect.

6. We obtain an estimate of the growth rate as $f{\sigma }_8=0.{440}_{-0.096}^{+0.080}$ at z_eff = 0.07. The refined estimate of the growth rate, derived by synthesizing the fitting outcomes from Turner et al. (2025), Lai et al. (2025), and this paper, yields a consensus value of fσ₈ = 0.450 ± 0.055. These are in agreement with the prediction of the ΛCDM and GR models.

7. We obtain an estimate of γ as $\gamma =0.{580}_{-0.110}^{+0.110}$ (corresponding to the consensus fσ₈), which is in agreement with the prediction of GR. Correspondingly, we find ${\mathrm{\Omega }}_{\mathrm{m}}=0.{301}_{-0.011}^{+0.011}$ and ${\sigma }_8=0.{834}_{-0.032}^{+0.032}$.

There are several promising avenues for future research building upon this work. (1) Our methodology can be applied to Data Release 2 of the DESI peculiar velocity survey, which is significantly more homogeneous and extensive than Data Release 1, thereby enabling a more precise constraint on fσ₈. (2) While the current analysis assumes the local plane-parallel approximation, future efforts will focus on incorporating wide-angle effects into both our models and estimators. (3) As is highlighted in Appendix G of this paper, the momentum power spectrum and cross power spectrum lack Galilean invariance, making them susceptible to potential systematic biases in the data. To address this limitation, we aim to transition toward measuring and modelling the pair-wise momentum power spectrum, as has been proposed by Chen et al. (2025). (4) We have not yet accounted for the gravitational redshift of galaxies. In forthcoming studies, we intend to extend our framework to measure and model this relativistic effect, following the approach outlined by Beutler & Di Dio (2020).

Data availability

The data used to generate the figures in this paper are available at https://zenodo.org/uploads/17672675

Acknowledgments

FQ and JB is supported by the funding from Excellence Initiative of Aix-Marseille University – A*MIDEX, a French “Investissements d’Avenir” program (AMX-20-CE-02 – DARKUNI). This material is based upon work supported by the U.S. Department of Energy (DOE), Office of Science, Office of High-Energy Physics, under Contract No. DE–AC02–05CH11231, and by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract. Additional support for DESI was provided by the U.S. National Science Foundation (NSF), Division of Astronomical Sciences under Contract No. AST-0950945 to the NSF’s National Optical-Infrared Astronomy Research Laboratory; the Science and Technology Facilities Council of the United Kingdom; the Gordon and Betty Moore Foundation; the Heising-Simons Foundation; the French Alternative Energies and Atomic Energy Commission (CEA); the National Council of Humanities, Science and Technology of Mexico (CONAHCYT); the Ministry of Science, Innovation and Universities of Spain (MICIU/AEI/10.13039/501100011033), and by the DESI Member Institutions: https://www.desi.lbl.gov/collaborating-institutions. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the U.S. National Science Foundation, the U.S. Department of Energy, or any of the listed funding agencies. The authors are honored to be permitted to conduct scientific research on I’oligam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation.

References

---

Appendix A: Peculiar velocity estimator

By employing the low-redshift approximation of the log-distance ratio, and under the assumption that the true peculiar velocities of galaxies are significantly smaller than their observed redshifts, the line-of-sight peculiar velocities can then be inferred from η through (Johnson et al. 2014; Watkins & Feldman 2015; Adams & Blake 2017; Howlett et al. 2017; Qin et al. 2018; Carreres et al. 2023)

$$\begin{array}{c}\hfill v=\frac{c{z}_{\mathrm{mod}}ln10}{1+z_{\mathrm{mod}}}\eta ,\end{array}$$(A.1)

where z_mod is given by (Davis & Scrimgeour 2014; Watkins & Feldman 2015)

$$\begin{array}{c}\hfill z_{\mathrm{mod}}=z[1+\frac{1}{2}(1-q_0)z-\frac{1}{6}(1-q_0-3q_0^2+1)z^2],\end{array}$$(A.2)

and where z is the observed redshift of a galaxy. The acceleration parameter is q₀ = 0.5(Ω_m − 2Ω_Λ).

Following the methodology outlined in Watkins & Feldman (2015) and Paper II, we adopt the following for Eq. 26

$$\begin{array}{c}\hfill ⟨v^2(\mathbf{r})⟩={\left(\frac{ln(10)cz}{1+z}{\sigma }_{\eta }\right)}^2+{300}^2{\mathrm{km}}^2{\mathrm{s}}^{-2},\end{array}$$(A.3)

where σ_η denotes the measurement error in the log-distance ratio η. 300² km² s⁻² represents the intrinsic scatter in peculiar velocities caused by non-linear galactic motions. As discussed in Appendix A of Paper II, variations in the assumed value of this intrinsic scatter have minimal impact on our results. The ⟨v²(r)⟩ term for the DESI-PV random catalogue is computed by Bautista et al. (2025).

Appendix B: Verifying $P_{\ell }^{\delta p}=-P_{\ell }^{\delta p\ast }$ for cross-power spectrum

In Sect. 3.2, we indicated that the estimator of the cross-power spectrum given in Eq. 17 should satisfy the condition $P_{\ell }^{\delta p}=-P_{\ell }^{\delta p\ast }$. This property will be verified in the present section. To begin with, let i denote the imaginary unit, such that i² = −1. Upon disregarding the shot-noise term and referring to Eq. 17, we obtain

$$\begin{array}{cc}\hfill P_{\ell }^{\delta p}\to & F^p{F}_{\ell }^{\delta \ast }-F^{\delta }{F}_{\ell }^{p\ast }\hfill \\ \hfill =& (\mathrm{Re}\{F^p\}+\mathrm{Im}\{F^p\}i)\times (\mathrm{Re}\{F_{\ell }^{\delta }\}-\mathrm{Im}\{F_{\ell }^{\delta }\}i)\hfill \\ \hfill & -(\mathrm{Re}\{F^{\delta }\}+\mathrm{Im}\{F^{\delta }\}i)\times (\mathrm{Re}\{F_{\ell }^p\}-\mathrm{Im}\{F_{\ell }^p\}i)\hfill \\ \hfill =& \mathrm{Re}\{F^p\}\mathrm{Re}\{F_{\ell }^{\delta }\}-\mathrm{Re}\{F^p\}\mathrm{Im}\{F_{\ell }^{\delta }\}i\hfill \\ \hfill & +\mathrm{Im}\{F^p\}\mathrm{Re}\{F_{\ell }^{\delta }\}i+\mathrm{Im}\{F^p\}\mathrm{Im}\{F_{\ell }^{\delta }\}\hfill \\ \hfill & -\mathrm{Re}\{F^{\delta }\}\mathrm{Re}\{F_{\ell }^p\}+\mathrm{Re}\{F^{\delta }\}\mathrm{Im}\{F_{\ell }^p\}i\hfill \\ \hfill & -\mathrm{Im}\{F^{\delta }\}\mathrm{Re}\{F_{\ell }^p\}i-\mathrm{Im}\{F^{\delta }\}\mathrm{Im}\{F_{\ell }^p\}\hfill \\ \hfill =& (\mathrm{Re}\{F^p\}\mathrm{Re}\{F_{\ell }^{\delta }\}+\mathrm{Im}\{F^p\}\mathrm{Im}\{F_{\ell }^{\delta }\}\hfill \\ \hfill & -\mathrm{Re}\{F^{\delta }\}\mathrm{Re}\{F_{\ell }^p\}-\mathrm{Im}\{F^{\delta }\}\mathrm{Im}\{F_{\ell }^p\})\hfill \\ \hfill & +(\mathrm{Im}\{F^p\}\mathrm{Re}\{F_{\ell }^{\delta }\}+\mathrm{Re}\{F^{\delta }\}\mathrm{Im}\{F_{\ell }^p\}\hfill \\ \hfill & -\mathrm{Re}\{F^p\}\mathrm{Im}\{F_{\ell }^{\delta }\}-\mathrm{Im}\{F^{\delta }\}\mathrm{Re}\{F_{\ell }^p\})i\hfill \end{array}$$(B.1)

which subsequently implies

$$\begin{array}{cc}\hfill P_{\ell }^{\delta p}& \propto \mathrm{Im}\{F^p{F}_{\ell }^{\delta \ast }-F^{\delta }{F}_{\ell }^{p\ast }\}\hfill \\ \hfill & =\mathrm{Im}\{F^p\}\mathrm{Re}\{F_{\ell }^{\delta }\}+\mathrm{Re}\{F^{\delta }\}\mathrm{Im}\{F_{\ell }^p\}\hfill \\ \hfill & -\mathrm{Re}\{F^p\}\mathrm{Im}\{F_{\ell }^{\delta }\}-\mathrm{Im}\{F^{\delta }\}\mathrm{Re}\{F_{\ell }^p\}\hfill \end{array}$$(B.2)

On the other hand, we have

$$\begin{array}{cc}\hfill P_{\ell }^{\delta p\ast }\to & F^{p\ast }{F}_{\ell }^{\delta }-F^{\delta \ast }{F}_{\ell }^p\hfill \\ \hfill =& (\mathrm{Re}\{F^p\}-\mathrm{Im}\{F^p\}i)\times (\mathrm{Re}\{F_{\ell }^{\delta }\}+\mathrm{Im}\{F_{\ell }^{\delta }\}i)\hfill \\ \hfill & -(\mathrm{Re}\{F^{\delta }\}-\mathrm{Im}\{F^{\delta }\}i)\times (\mathrm{Re}\{F_{\ell }^p\}+\mathrm{Im}\{F_{\ell }^p\}i)\hfill \\ \hfill =& \mathrm{Re}\{F^p\}\mathrm{Re}\{F_{\ell }^{\delta }\}+\mathrm{Re}\{F^p\}\mathrm{Im}\{F_{\ell }^{\delta }\}i\hfill \\ \hfill & -\mathrm{Im}\{F^p\}\mathrm{Re}\{F_{\ell }^{\delta }\}i+\mathrm{Im}\{F^p\}\mathrm{Im}\{F_{\ell }^{\delta }\}\hfill \\ \hfill & -\mathrm{Re}\{F^{\delta }\}\mathrm{Re}\{F_{\ell }^p\}-\mathrm{Re}\{F^{\delta }\}\mathrm{Im}\{F_{\ell }^p\}i\hfill \\ \hfill & +\mathrm{Im}\{F^{\delta }\}\mathrm{Re}\{F_{\ell }^p\}i-\mathrm{Im}\{F^{\delta }\}\mathrm{Im}\{F_{\ell }^p\}\hfill \end{array}$$(B.3)

which subsequently implies

$$\begin{array}{cc}\hfill P_{\ell }^{\delta p\ast }& \propto \mathrm{Im}\{F^{p\ast }{F}_{\ell }^{\delta }-F^{\delta \ast }{F}_{\ell }^p\}\hfill \\ \hfill & =\mathrm{Re}\{F^p\}\mathrm{Im}\{F_{\ell }^{\delta }\}-\mathrm{Im}\{F^p\}\mathrm{Re}\{F_{\ell }^{\delta }\}\hfill \\ \hfill & -\mathrm{Re}\{F^{\delta }\}\mathrm{Im}\{F_{\ell }^p\}+\mathrm{Im}\{F^{\delta }\}\mathrm{Re}\{F_{\ell }^p\}\hfill \\ \hfill & =-(\mathrm{Im}\{F^p\}\mathrm{Re}\{F_{\ell }^{\delta }\}+\mathrm{Re}\{F^{\delta }\}\mathrm{Im}\{F_{\ell }^p\}\hfill \\ \hfill & -\mathrm{Re}\{F^p\}\mathrm{Im}\{F_{\ell }^{\delta }\}-\mathrm{Im}\{F^{\delta }\}\mathrm{Re}\{F_{\ell }^p\}).\hfill \end{array}$$(B.4)

Comparing the above equation to Eq. B.2, one can logically infer that $P_{\ell }^{\delta p}=-P_{\ell }^{\delta p\ast }$.

Appendix C: The linear models of power spectrum and the models for the correlation functions

Kaiser (1987) developed a linear model for the density power spectrum, commonly referred to as the Kaiser Formula, given by

$$\begin{array}{c}\hfill P^{\delta }(k,\mu )={(b_1+f{\mu }^2)}^2{P}_L(k),\end{array}$$(C.1)

which has been extensively utilized over the past decades. Building upon this foundation, Koda et al. (2014) extended the analysis to the power spectrum of peculiar velocities as well as cross-power, they modelling the velocity power spectrum and density-velocity cross power spectrum using

$$\begin{array}{cc}& P^p(k,\mu )=\frac{{(aHf\mu )}^2}{k^2}{P}_L(k)\hfill \end{array}$$(C.2)

$$\begin{array}{cc}& P^{\delta p}(k,\mu )=(b_1+f{\mu }^2)\frac{iaHf\mu }{k}{P}_L(k).\hfill \end{array}$$(C.3)

Furthermore, they incorporated damping terms into the above linear models in order to more accurately account for non-linear motions of galaxies, given by

$$\begin{array}{cc}& P^{\delta }(k,\mu )={(b_1+f{\mu }^2)}^2{P}_L(k)D_g^2\hfill \end{array}$$(C.4)

$$\begin{array}{cc}& P^p(k,\mu )=\frac{{(aHf\mu )}^2}{k^2}{P}_L(k)D_p^2\hfill \end{array}$$(C.5)

$$\begin{array}{cc}& P^{\delta p}(k,\mu )=(b_1+f{\mu }^2)\frac{iaHf\mu }{k}{P}_L(k)D_g{D}_p\hfill \end{array}$$(C.6)

where the damping terms are defined by

$$\begin{array}{c}\hfill D_g=\frac{1}{\sqrt{1+\frac{1}{2}{(k\mu {\sigma }_g)}^2}},D_p=\frac{sin(k{\sigma }_v)}{k{\sigma }_v}\end{array}$$(C.7)

to account for non-linear motions of galaxies. The σ_g and σ_v are the non-linear velocity dispersion parameters for the density and velocity fields, respectively.

We initialize all loop terms–I_mn, K_mn, K_mn^s, σ₃ and σ₄–to zero, and similarly set the higher-order biasing and velocity dispersion parameters to zeros. In essence, we retain only fσ₈ and b₁σ₈. Consequently, Eq. 29, 30 and 31 reduce to the classic Kaiser formulae expressed in Eq. C.1,C.2 and C.3, respectively, as illustrated in Fig.C.2. For simplicity, we manually set fσ₈ = 0.437, b₁σ₈ = 0.564 and effective redshift z = 0 to generate this plot. In the middel panel of Fig. C.2, it becomes evident that the momentum power spectrum quadrupole P₂^p (represented by yellow dots and the yellow curve) surpasses the monopole P₀^p (depicted by blue dots and its curve). This arises due to the absence of window function convolution effects in our current consideration. Illustrated in Fig.C.1 is the linear matter power spectrum utilized in our analysis.

The mathematical framework for converting power spectrum to correlation functions primarily draws upon Turner et al. (2025). In Fig.C.3, we present the model correlation functions converted from the power spectrum models presented in Sect. 4 with parameters from Table 2.

Appendix D: The multipoles of the field function

The Eq. 19 can be further expressed in a more refined form as

$$\begin{array}{c}\hfill T_{\ell }(\mathbf{k})={\displaystyle \int (}\frac{k_x{r}_x+k_y{r}_y+k_z{r}_z}{\mathit{kr}}{)}^{\ell }F(\mathbf{r})e^{i\mathbf{k}·\mathbf{r}}{d}^{3}r\end{array}$$(D.1)

Thus, by directly expanding the term $(\frac{k_x{r}_x+k_y{r}_y+k_z{r}_z}{\mathit{kr}}{)}^{\ell }$ for ℓ = 0,  2 and 4, the even-multipole components of the Fourier-transformed field function are explicitly presented in

$$\begin{array}{c}\hfill T_0(\mathbf{k})={\displaystyle \int F}(\mathbf{r})e^{i\mathbf{k}·\mathbf{r}}{d}^{3}r=V\times F(\mathbf{k})\end{array}$$(D.2)

$$\begin{array}{cc}\hfill T_2(\mathbf{k})=& \frac{1}{k^2}(k_x^2{U}_{\mathit{xx}}+k_y^2{U}_{\mathit{yy}}+k_z^2{U}_{\mathit{zz}}\hfill \\ \hfill & +2(k_x{k}_y{U}_{\mathit{xy}}+k_x{k}_z{U}_{\mathit{xz}}+k_y{k}_z{U}_{\mathit{yz}}))\hfill \end{array}$$(D.3)

and

$$\begin{array}{cc}\hfill T_4(\mathbf{k})=& \frac{1}{k^4}(k_x^4{B}_{\mathit{xxx}}+k_y^4{B}_{\mathit{yyy}}+k_z^4{B}_{\mathit{zzz}}+4(k_x^3{k}_y{B}_{\mathit{xxy}}+k_x^3{k}_z{B}_{\mathit{xxz}}\hfill \\ \hfill & +k_y^3{k}_x{B}_{\mathit{yyx}}+k_y^3{k}_z{B}_{\mathit{yyz}}+k_z^3{k}_x{B}_{\mathit{zzx}}+k_z^3{k}_y{B}_{\mathit{zzy}})\hfill \\ \hfill & +6(k_x^2{k}_y^2{B}_{\mathit{xyy}}+k_x^2{k}_z^2{B}_{\mathit{xzz}}+k_y^2{k}_z^2{B}_{\mathit{yzz}})\hfill \\ \hfill & +12k_x{k}_y{k}_z(k_x{B}_{\mathit{xyz}}+k_y{B}_{\mathit{yxz}}+k_z{B}_{\mathit{zxy}})),\hfill \end{array}$$(D.4)

respectively, and where

$$\begin{array}{c}\hfill U_{\mathit{ij}}={\displaystyle \int \frac{r_i{r}_j}{r^2}}F(\mathbf{r})e^{i\mathbf{k}·\mathbf{r}}{d}^{3}r,B_{\mathit{ijn}}={\displaystyle \int \frac{r_i^2{r}_j{r}_n}{r^4}}F(\mathbf{r})e^{i\mathbf{k}·\mathbf{r}}{d}^{3}r,\end{array}$$(D.5)

where i, j, n = x, y, z.

Furthermore, guided by the theoretical framework proposed in Paper III, based on Eq. D.1, the odd-multipole components of the Fourier-transformed field function are formulated as described in

$$\begin{array}{cc}\hfill T_1(\mathbf{k})=& \frac{1}{k}(k_x{E}_x+k_y{E}_y+k_z{E}_z)\hfill \end{array}$$(D.6)

$$\begin{array}{cc}\hfill T_3(\mathbf{k})=& \frac{1}{k^3}(k_x^3{M}_{\mathit{xxx}}+k_y^3{M}_{\mathit{yyy}}+k_z^3{M}_{\mathit{zzz}}\hfill \\ \hfill & +3k_x^2{k}_y{M}_{\mathit{xxy}}+3k_x^2{k}_z{M}_{\mathit{xxz}}\hfill \\ \hfill & +3k_y^2{k}_x{M}_{\mathit{yyx}}+3k_y^2{k}_z{M}_{\mathit{yyz}}\hfill \\ \hfill & +3k_z^2{k}_x{M}_{\mathit{zzx}}+3k_z^2{k}_y{M}_{\mathit{zzy}}+6k_x{k}_y{k}_z{M}_{\mathit{xyz}}),\hfill \end{array}$$(D.7)

and where

$$\begin{array}{c}\hfill E_i={\displaystyle \int \frac{r_i}{r}}F(\mathbf{r})e^{i\mathbf{k}·\mathbf{r}}{d}^{3}r,M_{\mathit{ijn}}={\displaystyle \int \frac{r_i{r}_j{r}_n}{r^3}}F(\mathbf{r})e^{i\mathbf{k}·\mathbf{r}}{d}^{3}r\end{array}$$(D.8)

where i, j, n = x, y, z.

Appendix E: More tests using mocks

In Sect. 5.2 we demonstrate the adaptability of the power spectrum models by performing fits on [P₀^δ, P₂^δ, P₀^p, P₁^δp] using five free parameters [ fσ₈,  b₁σ₈,  b₂σ₈,  σ_vT,  σ_vS ], which suffices to maintain robustness and precision in the fitting process. Alternatively, we can separate the biasing parameters between the density and momentum fields–i.e., adopt seven free parameters $[f{\sigma }_8,b_1^{\delta }{\sigma }_8,b_2^{\delta }{\sigma }_8,b_1^p{\sigma }_8,b_2^p{\sigma }_8,{\sigma }_{\mathit{vT}},{\sigma }_{\mathit{vS}}]$–without significantly altering the resulting fσ₈; see Fig.E.1. Consequently, to maintain model simplicity and reduce parameters, we adopt shared biasing parameters for both fields.

Appendix F: Fitting fσ₈ of mocks from $P_{0,2}^{\delta }+P_{0,2}^p+P_1^{\delta p}$

In Sect. 5.2 we demonstrate the adaptability of the power spectrum models by performing fits on [P₀^δ, P₂^δ, P₀^p, P₁^δp]. Given that the momentum power spectrum quadrupole P₂^p derived from the FP mock average displays a sufficiently smooth behavior to support a robust fitting procedure, we extend our analysis to include P₂^p in the parameter estimation, thereby enabling a more comprehensive determination of fσ₈ based on [P₀^δ, P₂^δ, P₀^p, P₂^p, P₁^δp].

As we incorporate additional data into the fitting process, it becomes essential to enhance the flexibility of the power spectrum models accordingly. To achieve this, we introduce two additional free parameters. Specifically, we distinguish the linear biasing parameters between the density and momentum fields, and similarly differentiate the velocity dispersion parameters across these two fields. Consequently, the set of free parameters is defined as $\theta =[f{\sigma }_8,b_1^{\delta }{\sigma }_8,b_2{\sigma }_8,b_1^p{\sigma }_8,{\sigma }_{\delta ,vT}^2,{\sigma }_{\mathit{vS}}^2,{\sigma }_{p,vT}^2]$, resulting in a total of seven independent parameters. All other parameters remain shared between the density and momentum fields. Furthermore, to ensure model stability and prevent breakdown at nonlinear scales, we adopt a reduced k_max value of 0.25 h Mpc⁻¹, in contrast to that used in Sect. 5.2.

As illustrated in Fig. F.1, the filled circles represent the average of the power spectrum measured from 675 mocks, while the solid curves depict the theoretical models fit to the measurements. The MCMC resulting parameter estimates are show in Fig.F.2 and summarized in Table F.1. The fit value of $f{\sigma }_8=0.{427}_{-0.048}^{+0.056}$, which well aligns with the fiducial simulation value of fσ_8, fid = 0.466 at effective redshift z_eff = 0.2. This agreement demonstrates that our fitting methodology effectively recovers the true growth rate embedded in the mocks, and that the power spectrum models perform robustly up to the non-linear scale of k_max = 0.25 h Mpc⁻¹ under current conditions.

Appendix G: The Galilean transformation of power spectrum

As elaborated in Paper III, the momentum power spectrum P^p(k) and cross-power spectrum P^δp(k) are not invariant under the Galilean transformation; that is, their values depend on the bulk motion of the sample relative to the observational frame. Specifically, for the momentum power spectrum, if a constant bulk velocity shift ϵ is introduced to the line-of-sight peculiar velocities, i.e., by substituting v(r)→v(r)+ϵ, the momentum correlation yields

$$\begin{array}{cc}& ⟨p_G(\mathbf{r})p_G({\mathbf{r}}^{\prime })⟩\hfill \\ \hfill \equiv & ⟨(1+\delta (\mathbf{r}))(v(\mathbf{r})+ϵ)(1+\delta ({\mathbf{r}}^{\prime }))(v({\mathbf{r}}^{\prime })+ϵ)⟩\hfill \\ \hfill =& ⟨(1+\delta (\mathbf{r}))v(\mathbf{r})(1+\delta ({\mathbf{r}}^{\prime }))v({\mathbf{r}}^{\prime })⟩\hfill \\ \hfill & +⟨(1+\delta (\mathbf{r}))v(\mathbf{r})⟩ϵ+⟨(1+\delta (\mathbf{r}))v(\mathbf{r})\delta ({\mathbf{r}}^{\prime })⟩ϵ\hfill \\ \hfill & +⟨(1+\delta ({\mathbf{r}}^{\prime }))v({\mathbf{r}}^{\prime })⟩ϵ+⟨(1+\delta ({\mathbf{r}}^{\prime }))v({\mathbf{r}}^{\prime })\delta (\mathbf{r})⟩ϵ\hfill \\ \hfill & +{ϵ}^2+⟨\delta (\mathbf{r})⟩{ϵ}^2+⟨\delta ({\mathbf{r}}^{\prime })⟩{ϵ}^2+⟨\delta (\mathbf{r})\delta ({\mathbf{r}}^{\prime })⟩{ϵ}^2\hfill \\ \hfill =& ⟨p(\mathbf{r})p({\mathbf{r}}^{\prime })⟩+⟨p(\mathbf{r})⟩ϵ+⟨\delta ({\mathbf{r}}^{\prime })p(\mathbf{r})⟩ϵ+⟨p({\mathbf{r}}^{\prime })⟩ϵ+⟨\delta (\mathbf{r})p({\mathbf{r}}^{\prime })⟩ϵ\hfill \\ \hfill & +⟨\delta (\mathbf{r})⟩{ϵ}^2+⟨\delta ({\mathbf{r}}^{\prime })⟩{ϵ}^2+⟨\delta (\mathbf{r})\delta ({\mathbf{r}}^{\prime })⟩{ϵ}^2+{ϵ}^2\hfill \end{array}$$(G.1)

In this expression, the average of the density contrast field, i.e. ⟨δ(r)⟩ and ⟨δ(r′)⟩ are zeros. The average of the momentum field, i.e. ⟨p(r)⟩ and ⟨p(r′)⟩ are zeros too. Consequently, the above equation simplifies to

$$\begin{array}{cc}\hfill ⟨p_G(\mathbf{r})p_G({\mathbf{r}}^{\prime })⟩& =⟨p(\mathbf{r})p({\mathbf{r}}^{\prime })⟩+⟨\delta ({\mathbf{r}}^{\prime })p(\mathbf{r})⟩ϵ+⟨\delta (\mathbf{r})p({\mathbf{r}}^{\prime })⟩ϵ\hfill \\ \hfill & +⟨\delta (\mathbf{r})\delta ({\mathbf{r}}^{\prime })⟩{ϵ}^2+{ϵ}^2\hfill \end{array}$$(G.2)

Applying the Fourier transform to the aforementioned equation yields

$$\begin{array}{cc}& \frac{1}{{(2\pi )}^6}{\displaystyle \int ⟨p_G(\mathbf{k})p_G^{\ast }({\mathbf{k}}^{\prime })⟩}{e}^{-i\mathbf{k}·\mathbf{r}}{e}^{i{\mathbf{k}}^{\prime }·{\mathbf{r}}^{\prime }}{d}^{3}k{d}^3{k}^{\prime }\hfill \\ \hfill =& \frac{1}{{(2\pi )}^6}{\displaystyle \int [}⟨p(\mathbf{k})p^{\ast }({\mathbf{k}}^{\prime })⟩+⟨{\delta }^{\ast }({\mathbf{k}}^{\prime })p(\mathbf{k})⟩ϵ+⟨\delta (\mathbf{k})p^{\ast }({\mathbf{k}}^{\prime })⟩ϵ\hfill \\ \hfill & +⟨\delta (\mathbf{k}){\delta }^{\ast }({\mathbf{k}}^{\prime })⟩{ϵ}^2]e^{-i\mathbf{k}·\mathbf{r}}{e}^{i{\mathbf{k}}^{\prime }·{\mathbf{r}}^{\prime }}{d}^{3}k{d}^3{k}^{\prime }+{ϵ}^2{\delta }^D(\mathbf{k}){\delta }^D({\mathbf{k}}^{\prime })\hfill \end{array}$$(G.3)

in the final term, the δ^D(k) and δ^D(k′) cannot equal unity simultaneously, as we are exclusively considering the case k ≠ k′ in the context of two-point statistic. Consequently, δ^D(k)δ^D(k′) = 0, and the above equation indicates

$$\begin{array}{cc}\hfill ⟨p_G(\mathbf{k})p_G({\mathbf{k}}^{\prime })⟩& =⟨p(\mathbf{k})p^{\ast }({\mathbf{k}}^{\prime })⟩+⟨{\delta }^{\ast }({\mathbf{k}}^{\prime })p(\mathbf{k})⟩ϵ\hfill \\ \hfill & +⟨\delta (\mathbf{k})p^{\ast }({\mathbf{k}}^{\prime })⟩ϵ+⟨\delta (\mathbf{k}){\delta }^{\ast }({\mathbf{k}}^{\prime })⟩{ϵ}^2\hfill \end{array}$$(G.4)

Substituting the definitions of the power spectrum provided in Eq. 3, 4 and 5 into the above equation leads to

$$\begin{array}{c}\hfill P_G^p(\mathbf{k})=P^p(\mathbf{k})+P^{\delta p\ast }(\mathbf{k})ϵ+P^{\delta p}(\mathbf{k})ϵ+{ϵ}^2{P}^{\delta }(\mathbf{k})\end{array}$$(G.5)

Given that P^δp = −P^δp* (see Appendix B), the above equation further reduces to

$$\begin{array}{c}\hfill P_G^p(\mathbf{k})=P^p(\mathbf{k})+{ϵ}^2{P}^{\delta }(\mathbf{k})\end{array}$$(G.6)

which represents the Galilean-transformed form of the momentum power spectrum. Analogously, the corresponding Galilean-transformed cross-power spectrum is expressed as

$$\begin{array}{c}\hfill P_G^{\delta p}(\mathbf{k})=P^{\delta p}(\mathbf{k})+ϵP^{\delta }(\mathbf{k})\end{array}$$(G.7)

The bulk velocity ϵ comprises two possible components: one stems from coherent bulk motion of galaxies driven by local gravitational fluctuations, and the other originates from potential systematic errors induced by observational or measurement inaccuracies.

We investigate whether the fitting outcomes exhibit sensitivity to the parameter ϵ. By substituting the power spectrum models described in Sect. 4 into Eqs. G.6 and G.7, we obtain the Galilean-transformed versions of these models. These transformed models are then fitted to the measured power spectrum of BGS and DESI-PV data. In this analysis, ϵ is treated as a free parameter with a flat prior over the broad range of ϵ ∈ [ − 1000, 1000] km s⁻¹. The resulting fits are visualized through pink-colored histograms and 2D contour plots in Fig.G.1. For comparative purposes, the blue-colored histograms and 2D contours represent the fitting results obtained without applying the Galilean transformation (identical to those shown in Fig.11). The results from both approaches are in close agreement, and the fit value for ϵ is determined to be $ϵ=-15.{94}_{-35.95}^{+35.95}$ km s⁻¹, which is statistically consistent with zero. Consequently, the influence of the Galilean transformation is found to be negligible within the context of our analysis.

Appendix H: Deriving the estimator of power spectrum and the window function convolution matrix

Corresponding to Sects. 3.2 and 4.3, in this section, we present a detailed derivation of the estimator for the power spectrum and its associated window function convolution matrix by taking the cross-power spectrum as a representative example. The mathematical framework primarily draws upon Sect. 3.1 of Paper I, Sect. 3.2 of Paper III, and Sect. 2.2 of Blake et al. (2018).

H.1. Power spectrum estimator

We begin by introducing the Legendre transformation, denoted as (Blake et al. 2018)

$$\begin{array}{c}\hfill F_{\ell }(\mathbf{r})=(2\ell +1){\displaystyle \int F}(\mathbf{r})L_{\ell }(\widehat{\mathbf{k}}·\widehat{\mathbf{r}})\frac{d{\mathrm{\Omega }}_k}{4\pi }\end{array}$$(H.1)

Following the methodology outlined in Yamamoto et al. (2006) and Sect. 3.2 of Paper III and employing the Fourier transform definition given in Eq. 10, we formulate our estimator for the cross-power spectrum under the "local plane-parallel approximation", i.e. $L_{\ell }(\widehat{\mathbf{k}}·\widehat{\mathbf{r}})=L_{\ell }(\widehat{\mathbf{k}}·\widehat{\mathbf{r}\prime })$, as

$$\begin{array}{cc}\hfill |F^{\delta }(k)F_{\ell }^{p\ast }(k)|=& \frac{1}{V}{\displaystyle \int F^{\delta }}(\mathbf{r})e^{i\mathbf{k}·\mathbf{r}}{d}^{3}r\frac{1}{V^{\prime }}{\displaystyle \int [}(2\ell +1){\displaystyle \int F^{p\ast }}({\mathbf{r}}^{\prime })L_{\ell }^{\ast }(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})\frac{d{\mathrm{\Omega }}_k}{4\pi }]e^{-i\mathbf{k}·{\mathbf{r}}^{\prime }}{d}^3{r}^{\prime }\hfill \\ \hfill =& \frac{2\ell +1}{V{V}^{\prime }}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }\int d^{3}r\int d^3}{r}^{\prime }{F}^{\delta }(\mathbf{r})F^{p\ast }({\mathbf{r}}^{\prime })L_{\ell }^{\ast }(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})e^{i\mathbf{k}·(\mathbf{r}-{\mathbf{r}}^{\prime })}.\hfill \end{array}$$(H.2)

Subsequently, the cross-correlation between the density and momentum fields is expressed as

$$\begin{array}{c}\hfill ⟨|F^{\delta }(k)F_{\ell }^{p\ast }(k)|⟩=\frac{2\ell +1}{V{V}^{\prime }}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }\int d^3}r{\displaystyle \int d^3}{r}^{\prime }⟨F^{\delta }(\mathbf{r})F^{p\ast }({\mathbf{r}}^{\prime })⟩L_{\ell }^{\ast }(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})e^{i\mathbf{k}·(\mathbf{r}-{\mathbf{r}}^{\prime })}\end{array}$$(H.3)

Building upon the methodology outlined in Sect. 3.2 of Paper III, and employing Eq. 6 and 8, we derive the cross-correlation between the density and momentum fields as

$$\begin{array}{c}\hfill ⟨F^{\delta }(\mathbf{r})F^{p\ast }({\mathbf{r}}^{\prime })⟩=\langle \frac{w_{\delta }(\mathbf{r})[n_{\delta }(\mathbf{r})-\alpha n_s(\mathbf{r})]}{A_{\delta }}\frac{w_p({\mathbf{r}}^{\prime })n_p({\mathbf{r}}^{\prime })v({\mathbf{r}}^{\prime })}{A_p}\rangle =\frac{1}{A_{\delta }{A}_p}{w}_{\delta }(\mathbf{r})w_p({\mathbf{r}}^{\prime })(⟨n_{\delta }(\mathbf{r})n_p({\mathbf{r}}^{\prime })v({\mathbf{r}}^{\prime })⟩-\alpha ⟨n_s(\mathbf{r})n_p({\mathbf{r}}^{\prime })v({\mathbf{r}}^{\prime })⟩)\end{array}$$(H.4)

Drawing from the approach described in Park & Park 2006 and Paper I, the first term within the brackets of Eq. H.4 can be rewritten as $⟨n_{\delta }(\mathbf{r})n_p(\mathbf{r}\prime )v(\mathbf{r}\prime )⟩={\overline{n}}_{\delta }(\mathbf{r}){\overline{n}}_p(\mathbf{r}\prime ){\xi }_{\delta p}(|\mathbf{r}-\mathbf{r}\prime |)+min\{{\overline{n}}_{\delta }(\mathbf{r}),{\overline{n}}_p(\mathbf{r})\}⟨v(\mathbf{r})⟩{\delta }^D(|\mathbf{r}-\mathbf{r}\prime |)$. Given that the cross-correlation between galaxies and random points is zero, the second term within the brackets of Eq. H.4 becomes zero, i.e. α⟨n_s(r)n_p(r′)v(r′)⟩ = 0. Consequently, we obtain

$$\begin{array}{c}\hfill ⟨F^{\delta }(\mathbf{r})F^{p\ast }({\mathbf{r}}^{\prime })⟩=\frac{1}{A_{\delta }{A}_p}{w}_{\delta }(\mathbf{r})w_p({\mathbf{r}}^{\prime })[{\overline{n}}_{\delta }(\mathbf{r}){\overline{n}}_p({\mathbf{r}}^{\prime }){\xi }_{\delta p}(|\mathbf{r}-{\mathbf{r}}^{\prime }|)+min\{{\overline{n}}_{\delta }(\mathbf{r}),{\overline{n}}_p(\mathbf{r})\}⟨v({\mathbf{r}}^{\prime })⟩{\delta }^D(|\mathbf{r}-{\mathbf{r}}^{\prime }|)]\end{array}$$(H.5)

By substituting the above equation into Eq. H.3 in place of ⟨F^δ(r)F^p*(r′)⟩, we arrive at

$$\begin{array}{cc}\hfill ⟨|F^{\delta }(k)F_{\ell }^{p\ast }(k)|⟩=& \frac{1}{A_{\delta }{A}_p}\frac{2\ell +1}{V}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }\int d^3}r[\frac{1}{V^{\prime }}{\displaystyle \int d^3}{r}^{\prime }{w}_{\delta }(\mathbf{r})w_p({\mathbf{r}}^{\prime }){\overline{n}}_{\delta }(\mathbf{r}){\overline{n}}_p({\mathbf{r}}^{\prime }){\xi }_{\delta p}(|\mathbf{r}-{\mathbf{r}}^{\prime }|)L_{\ell }^{\ast }(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})e^{i\mathbf{k}·(\mathbf{r}-{\mathbf{r}}^{\prime })}\hfill \\ \hfill & +\frac{1}{V^{\prime }}{\displaystyle \int d^3}{r}^{\prime }{w}_{\delta }(\mathbf{r})w_p(\mathbf{r})min\{{\overline{n}}_{\delta }(\mathbf{r}),{\overline{n}}_p(\mathbf{r})\}⟨v(\mathbf{r})⟩{\delta }^D(\mathbf{r}-{\mathbf{r}}^{\prime })L_{\ell }^{\ast }(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})e^{i\mathbf{k}·(\mathbf{r}-{\mathbf{r}}^{\prime })}]\hfill \end{array}$$(H.6)

Upon integrating the second term over the Dirac delta function about r′ using the relation $f(\mathbf{r})=\frac{1}{V}{\displaystyle \int f}(\mathbf{r}){\delta }^D(\mathbf{r}-\mathbf{r}\prime )d^{3}r\prime$, we obtain

$$\begin{array}{cc}\hfill ⟨|F^{\delta }(k)F_{\ell }^{p\ast }(k)|⟩=& \frac{1}{A_{\delta }{A}_p}\frac{2\ell +1}{V}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }\int d^3}r[\frac{1}{V^{\prime }}{\displaystyle \int d^3}{r}^{\prime }{w}_{\delta }(\mathbf{r})w_p({\mathbf{r}}^{\prime }){\overline{n}}_{\delta }(\mathbf{r}){\overline{n}}_p({\mathbf{r}}^{\prime }){\xi }_{\delta p}(|\mathbf{r}-{\mathbf{r}}^{\prime }|)L_{\ell }^{\ast }(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})e^{i\mathbf{k}·(\mathbf{r}-{\mathbf{r}}^{\prime })}\hfill \\ \hfill & +w_{\delta }(\mathbf{r})w_p(\mathbf{r})min\{{\overline{n}}_{\delta }(\mathbf{r}),{\overline{n}}_p(\mathbf{r})\}⟨v(\mathbf{r})⟩L_{\ell }^{\ast }(\widehat{\mathbf{k}}·\widehat{\mathbf{r}})]\hfill \\ \hfill =& \frac{1}{A_{\delta }{A}_p}\frac{2\ell +1}{V}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }}[{\displaystyle \int d^3}r\frac{1}{V^{\prime }}{\displaystyle \int d^3}{r}^{\prime }{w}_{\delta }(\mathbf{r})w_p({\mathbf{r}}^{\prime }){\overline{n}}_{\delta }(\mathbf{r}){\overline{n}}_p({\mathbf{r}}^{\prime }){\xi }_{\delta p}(|\mathbf{r}-{\mathbf{r}}^{\prime }|)L_{\ell }^{\ast }(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})e^{i\mathbf{k}·(\mathbf{r}-{\mathbf{r}}^{\prime })}\hfill \\ \hfill & +{\displaystyle \int w_{\delta }}(\mathbf{r})w_p(\mathbf{r})min\{{\overline{n}}_{\delta }(\mathbf{r}),{\overline{n}}_p(\mathbf{r})\}⟨v(\mathbf{r})⟩L_{\ell }^{\ast }(\widehat{\mathbf{k}}·\widehat{\mathbf{r}})d^{3}r]\hfill \end{array}$$(H.7)

We define the function

$$\begin{array}{c}\hfill N^{\delta p}\equiv {\displaystyle \int w_{\delta }}(\mathbf{r})w_p(\mathbf{r})min\{{\overline{n}}_{\delta }(\mathbf{r}),{\overline{n}}_p(\mathbf{r})\}⟨v(\mathbf{r})⟩L_{\ell }^{\ast }(\widehat{\mathbf{k}}·\widehat{\mathbf{r}})d^{3}r=\frac{A_{\delta }{A}_{p}V}{2\ell +1}\times {\mathcal{N}}_{\ell }^{\delta p}\end{array}$$(H.8)

where the shot-noise term ${\mathcal{N}}_{\ell }^{\delta p}$ is formally introduced as Eq. 18 of this paper, allowing us to simplify Eq. H.7 into

$$\begin{array}{c}\hfill ⟨|F^{\delta }(k)F_{\ell }^{p\ast }(k)|⟩=\frac{1}{A_{\delta }{A}_p}\frac{2\ell +1}{V}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }}[{\displaystyle \int d^3}r\frac{1}{V^{\prime }}{\displaystyle \int d^3}{r}^{\prime }{w}_{\delta }(\mathbf{r})w_p({\mathbf{r}}^{\prime }){\overline{n}}_{\delta }(\mathbf{r}){\overline{n}}_p({\mathbf{r}}^{\prime }){\xi }_{\delta p}(|\mathbf{r}-{\mathbf{r}}^{\prime }|)L_{\ell }^{\ast }(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})e^{i\mathbf{k}·(\mathbf{r}-{\mathbf{r}}^{\prime })}+N^{\delta p}]\end{array}$$(H.9)

The two-point correlation function maintains an intrinsic relationship with the power spectrum through

$$\begin{array}{c}\hfill {\xi }_{\delta p}(|\mathbf{r}-{\mathbf{r}}^{\prime }|)\equiv \frac{1}{{(2\pi )}^3}{\displaystyle \int P^{\delta p}}(\mathbf{k})e^{-i\mathbf{k}·(\mathbf{r}-{\mathbf{r}}^{\prime })}{d}^{3}k\end{array}$$(H.10)

Incorporating this relationship into Eq. H.9 to substitute ξ_δp leads to

$$\begin{array}{cc}& ⟨|F^{\delta }(k)F_{\ell }^{p\ast }(k)|⟩\hfill \\ \hfill =& \frac{1}{A_{\delta }{A}_p}\frac{2\ell +1}{V}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }}[{\displaystyle \int d^3}r\frac{1}{V^{\prime }}{\displaystyle \int d^3}{r}^{\prime }{w}_{\delta }(\mathbf{r})w_p({\mathbf{r}}^{\prime }){\overline{n}}_{\delta }(\mathbf{r}){\overline{n}}_p({\mathbf{r}}^{\prime })({\displaystyle \int P^{\delta p}}({\mathbf{k}}^{\prime })e^{-i{\mathbf{k}}^{\prime }·(\mathbf{r}-{\mathbf{r}}^{\prime })}\frac{d^3{k}^{\prime }}{{(2\pi )}^3})L_{\ell }^{\ast }(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})e^{i\mathbf{k}·(\mathbf{r}-{\mathbf{r}}^{\prime })}+N^{\delta p}]\hfill \\ \hfill =& \frac{1}{A_{\delta }{A}_p}\frac{2\ell +1}{V}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }}[({\displaystyle \int w_{\delta }}(\mathbf{r}){\overline{n}}_{\delta }(\mathbf{r})e^{i(\mathbf{k}-{\mathbf{k}}^{\prime })·\mathbf{r}}{d}^{3}r)\frac{1}{V^{\prime }}{\displaystyle \int d^3}{r}^{\prime }{\displaystyle \int \frac{d^3{k}^{\prime }}{{(2\pi )}^3}}{P}^{\delta p}({\mathbf{k}}^{\prime })L_{\ell }^{\ast }(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})e^{-i(\mathbf{k}-{\mathbf{k}}^{\prime })·{\mathbf{r}}^{\prime }}{w}_p({\mathbf{r}}^{\prime }){\overline{n}}_p({\mathbf{r}}^{\prime })+N^{\delta p}]\hfill \end{array}$$(H.11)

We define the function

$$\begin{array}{c}\hfill G^{\delta }(\mathbf{k}-{\mathbf{k}}^{\prime })\equiv {\displaystyle \int w_{\delta }}(\mathbf{r}){\overline{n}}_{\delta }(\mathbf{r})e^{i(\mathbf{k}-{\mathbf{k}}^{\prime })·\mathbf{r}}{d}^{3}r\end{array}$$(H.12)

as proposed in Eq. 51, Eq. H.11 can be reduced to

$$\begin{array}{c}\hfill ⟨|F^{\delta }(k)F_{\ell }^{p\ast }(k)|⟩=\frac{1}{A_{\delta }{A}_p}\frac{2\ell +1}{V}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }}[\frac{1}{V^{\prime }}{\displaystyle \int d^3}{r}^{\prime }{\displaystyle \int \frac{d^3{k}^{\prime }}{{(2\pi )}^3}}{P}^{\delta p}({\mathbf{k}}^{\prime })L_{\ell }^{\ast }(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})e^{-i(\mathbf{k}-{\mathbf{k}}^{\prime })·{\mathbf{r}}^{\prime }}{w}_p({\mathbf{r}}^{\prime }){\overline{n}}_p({\mathbf{r}}^{\prime })G^{\delta }(\mathbf{k}-{\mathbf{k}}^{\prime })+N^{\delta p}]\end{array}$$(H.13)

Expanding the power spectrum P^δp(k′) in terms of Legendre polynomials using Eq. 12 of this paper, we derive

$$\begin{array}{cc}& ⟨|F^{\delta }(k)F_{\ell }^{p\ast }(k)|⟩\hfill \\ \hfill =& \frac{1}{A_{\delta }{A}_p}\frac{2\ell +1}{V}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }}[\frac{1}{V^{\prime }}{\displaystyle \int d^3}{r}^{\prime }{\displaystyle \int \frac{d^3{k}^{\prime }}{{(2\pi )}^3}}({\displaystyle \sum _{{\ell }^{\prime }}}{P}_{{\ell }^{\prime }}^{\delta p}(k^{\prime })L_{{\ell }^{\prime }}({\widehat{\mathbf{k}}}^{\prime }·\widehat{{\mathbf{r}}^{\prime }}))L_{\ell }^{\ast }(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})e^{-i(\mathbf{k}-{\mathbf{k}}^{\prime })·{\mathbf{r}}^{\prime }}{w}_p({\mathbf{r}}^{\prime }){\overline{n}}_p({\mathbf{r}}^{\prime })G^{\delta }(\mathbf{k}-{\mathbf{k}}^{\prime })+N^{\delta p}]\hfill \\ \hfill =& \frac{2\ell +1}{A_{\delta }{A}_{p}V}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }}[{\displaystyle \sum _{{\ell }^{\prime }}}{\displaystyle \int \frac{d^3{k}^{\prime }}{{(2\pi )}^3}}{P}_{{\ell }^{\prime }}^{\delta p}(k^{\prime })G^{\delta }(\mathbf{k}-{\mathbf{k}}^{\prime })(\frac{1}{V^{\prime }}{\displaystyle \int d^3}{r}^{\prime }{w}_p({\mathbf{r}}^{\prime }){\overline{n}}_p({\mathbf{r}}^{\prime })L_{{\ell }^{\prime }}({\widehat{\mathbf{k}}}^{\prime }·\widehat{{\mathbf{r}}^{\prime }})L_{\ell }^{\ast }(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})e^{-i(\mathbf{k}-{\mathbf{k}}^{\prime })·{\mathbf{r}}^{\prime }})+N^{\delta p}]\hfill \end{array}$$(H.14)

Further defining the function

$$\begin{array}{c}\hfill S_{\ell {\ell }^{\prime }}^p(\mathbf{k},{\mathbf{k}}^{\prime })\equiv \frac{1}{V}{\displaystyle \int w_p}(\mathbf{r}){\overline{n}}_p(\mathbf{r})L_{\ell }(\widehat{\mathbf{k}}·\widehat{\mathbf{r}})L_{{\ell }^{\prime }}^{\ast }({\widehat{\mathbf{k}}}^{\prime }·\widehat{\mathbf{r}})e^{i(\mathbf{k}-{\mathbf{k}}^{\prime })·\mathbf{r}}{d}^{3}r\end{array}$$(H.15)

the conjugate form of the above equation enables us to reformulate Eq. H.14 as

$$\begin{array}{c}\hfill ⟨|F^{\delta }(k)F_{\ell }^{p\ast }(k)|⟩=\frac{2\ell +1}{A_{\delta }{A}_{p}V}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }}[{\displaystyle \sum _{{\ell }^{\prime }}}{\displaystyle \int \frac{d^3{k}^{\prime }}{{(2\pi )}^3}}{P}_{{\ell }^{\prime }}^{\delta p}(k^{\prime })G^{\delta }(\mathbf{k}-{\mathbf{k}}^{\prime })S_{\ell {\ell }^{\prime }}^{p\ast }(\mathbf{k},{\mathbf{k}}^{\prime })+N^{\delta p}]\end{array}$$(H.16)

which is similar to Equation 15 of Paper III⁸. We define the window function as

$$\begin{array}{c}\hfill \mathbf{W}\equiv \frac{2\ell +1}{A_{\delta }{A}_{p}V}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }}[{\displaystyle \sum _{{\ell }^{\prime }}}{\displaystyle \int \frac{d^3{k}^{\prime }}{{(2\pi )}^3}}{G}^{\delta }(\mathbf{k}-{\mathbf{k}}^{\prime })S_{\ell {\ell }^{\prime }}^{p\ast }(\mathbf{k},{\mathbf{k}}^{\prime })]\end{array}$$(H.17)

in the absence of this component, Eq. H.16 collapses to

$$\begin{array}{c}\hfill ⟨|F^{\delta }(k)F_{\ell }^{p\ast }(k)|⟩=P_{\ell }^{\delta p}+\frac{(2\ell +1)N^{\delta p}}{A_{\delta }{A}_{p}V}=P_{\ell }^{\delta p}+{\mathcal{N}}_{\ell }^{\delta p}\end{array}$$(H.18)

i.e. we can formally express our estimator of the density-momentum cross power spectrum as $P_{\ell }^{\delta p}=|F^{\delta }(k)F_{\ell }^{p\ast }(k)|-{\mathcal{N}}_{\ell }^{\delta p}$ with $|F^{\delta }(k)F_{\ell }^{p\ast }(k)|$ given by Eq. H.2. Consequently, we re-formulate our estimator as

$$\begin{array}{c}\hfill P_{\ell }^{\delta p}=\frac{1}{2}\mathrm{Im}\{|F^p(k)F_{\ell }^{\delta \ast }(k)|-|F^{\delta }(k)F_{\ell }^{p\ast }(k)|\}-{\mathcal{N}}_{\ell }^{\delta p}\end{array}$$(H.19)

ensuring its symmetry under $P_{\ell }^{\delta p}=-P_{\ell }^{\delta p\ast }$ as established in Appendix B.

H.2. Window function

We can expand the Legendre functions in terms of spherical harmonics following the approach outlined in Blake et al. (2018)

$$\begin{array}{c}\hfill L_{\ell }(\widehat{\mathbf{k}}·\widehat{\mathbf{r}})=\frac{4\pi }{2\ell +1}{\displaystyle \sum _{m=-\ell }^{\ell }}{Y}_m^{\ell \ast }(\widehat{\mathbf{k}})Y_m^{\ell }(\widehat{\mathbf{r}})\end{array}$$(H.20)

Substituting this expansion and its conjugate into Eq. H.15 yields

$$\begin{array}{c}\hfill S_{\ell {\ell }^{\prime }}^p(\mathbf{k},{\mathbf{k}}^{\prime })=\frac{4\pi }{2\ell +1}\frac{4\pi }{2{\ell }^{\prime }+1}{\displaystyle \sum _{m=-\ell }^{\ell }}{Y}_m^{\ell \ast }(\widehat{\mathbf{k}}){\displaystyle \sum _{m^{\prime }=-{\ell }^{\prime }}^{{\ell }^{\prime }}}{Y}_{m^{\prime }}^{{\ell }^{\prime }}({\widehat{\mathbf{k}}}^{\prime }){\stackrel{\sim }{S}}_{m{m}^{\prime }}^{p,\ell {\ell }^{\prime }}(\mathbf{k})\end{array}$$(H.21)

where ${\stackrel{\sim }{S}}_{mm\prime }^{\ell \ell \prime }(\mathbf{k})$ is defined by Eq. 50 of this paper. Substituting this result back into Eq. H.17 leads to

$$\begin{array}{c}\hfill \mathbf{W}=\frac{1}{A_{\delta }{A}_p}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }\sum _{m=-\ell }^{\ell }}{Y}_m^{\ell }(\widehat{\mathbf{k}}){\displaystyle \sum _{{\ell }^{\prime }}}\frac{{(4\pi )}^2}{2{\ell }^{\prime }+1}{\displaystyle \sum _{m^{\prime }=-{\ell }^{\prime }}^{{\ell }^{\prime }}}[\frac{1}{V}{\displaystyle \int \frac{d^3{k}^{\prime }}{{(2\pi )}^3}}{Y}_{m^{\prime }}^{{\ell }^{\prime }\ast }({\widehat{\mathbf{k}}}^{\prime })G^{\delta }(\mathbf{k}-{\mathbf{k}}^{\prime }){\stackrel{\sim }{S}}_{m{m}^{\prime }}^{p,\ell {\ell }^{\prime }\ast }(\mathbf{k})]\end{array}$$(H.22)

which corresponds to Eq. 48 of this paper, and serves as the foundation for computing the window function convolution matrix. The aforementioned expression of W corresponds to the term $|F^{\delta }(k)F_{\ell }^{p\ast }(k)|$. However, as suggested by Blake et al. (2018), in order to derive the corresponding W that conforms to the format specified in Eq. H.19, we just need to substitute the term $G^{\delta }(\mathbf{k}-\mathbf{k}\prime ){\stackrel{\sim }{S}}_{mm\prime }^{p,\ell \ell \prime *}(\mathbf{k})$ in the above equation with Eq. 53. Notably, the factor V in Equation 19 of Blake et al. (2018) (and Equation 38 of Paper III) is differ from the above Eq. H.22 due to their definition of Fourier transformation and G^δ diverges from our definition presented in Eq. 10 and 51.

H.3. The expression of multipoles

To illustrate the methodology for computing Eq. 20 to 24, we consider the ℓ = 4 auto-power spectrum as a representative example for computing the estimators defined in Eq. 22. In this specific case, Eq. H.2 simplifies to

$$\begin{array}{cc}& |F(k)F_4^{\ast }(k)|\hfill \\ \hfill =& \frac{9}{V}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }}{\displaystyle \int d^3}r\frac{1}{V^{\prime }}{\displaystyle \int d^3}{r}^{\prime }F(\mathbf{r})F^{\ast }({\mathbf{r}}^{\prime })\times \frac{1}{8}(35{(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})}^4-30{(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})}^2+3{(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})}^0)\times e^{i\mathbf{k}·(\mathbf{r}-{\mathbf{r}}^{\prime })}\hfill \\ \hfill =& \frac{9}{8}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }}(\frac{1}{V}{\displaystyle \int F}(\mathbf{r})e^{i\mathbf{k}·\mathbf{r}}{d}^{3}r)\frac{1}{V^{\prime }}{\displaystyle \int d^3}{r}^{\prime }{F}^{\ast }({\mathbf{r}}^{\prime })(35{(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})}^4-30{(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})}^2+3{(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})}^0)e^{-i\mathbf{k}·{\mathbf{r}}^{\prime }}\hfill \end{array}$$(H.23)

By applying Eq. 10, the resulting expression becomes

$$\begin{array}{cc}& |F(k)F_4^{\ast }(k)|\hfill \\ \hfill =& \frac{9}{8}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }}F(\mathbf{k})\frac{1}{V^{\prime }}{\displaystyle \int d^3}{r}^{\prime }{F}^{\ast }({\mathbf{r}}^{\prime })(35{(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})}^4-30{(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})}^2+3{(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})}^0)e^{-i\mathbf{k}·{\mathbf{r}}^{\prime }}\hfill \\ \hfill =& \frac{9}{8V^{\prime }}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }}F(\mathbf{k})[35{\displaystyle \int {(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})}^4}{F}^{\ast }({\mathbf{r}}^{\prime })e^{-i\mathbf{k}·{\mathbf{r}}^{\prime }}{d}^3{r}^{\prime }-30{\displaystyle \int {(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})}^2}{F}^{\ast }({\mathbf{r}}^{\prime })e^{-i\mathbf{k}·{\mathbf{r}}^{\prime }}{d}^3{r}^{\prime }+3{\displaystyle \int {(\widehat{\mathbf{k}}·\widehat{{\mathbf{r}}^{\prime }})}^0}{F}^{\ast }({\mathbf{r}}^{\prime })e^{-i\mathbf{k}·{\mathbf{r}}^{\prime }}{d}^3{r}^{\prime }]\hfill \end{array}$$(H.24)

Further incorporating the conjugate form of Eq. 19 for the cases where ℓ = 0, 2, 4, this expression reduces to

$$\begin{array}{c}\hfill |F(k)F_4^{\ast }(k)|=\frac{9}{8V^{\prime }}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }}F(\mathbf{k})[35T_4^{\ast }(\mathbf{k})-30T_2^{\ast }(\mathbf{k})+3T_0^{\ast }(\mathbf{k})]\end{array}$$(H.25)

Substituting this expansion into Eq. 13 or Eq. 14 yields

$$\begin{array}{c}\hfill P_4(k)=\frac{9}{8V^{\prime }}{\displaystyle \int \frac{d{\mathrm{\Omega }}_k}{4\pi }}F(\mathbf{k})[35T_4^{\ast }(\mathbf{k})-30T_2^{\ast }(\mathbf{k})+3T_0^{\ast }(\mathbf{k})]-{\mathcal{N}}_4\end{array}$$(H.26)

with 𝒩₄ ≈ 0, the above equation collapses precisely to Eq. 22.

All Tables

All Figures