© The Authors 2025

1. Introduction

In recent years, a large number of cosmological analyses have used peculiar velocities of galaxies to infer the growth rate of large-scale structures, denoted as fσ₈ (see e.g., Turner 2024 for a recent review). These velocities can be determined using a variety of methods, all of which involve determining redshift-independent distance indicators and spectroscopic redshifts.

Two of the most widely used distance indicators are the Tully-Fisher (TF) relation (Tully & Fisher 1977) and the Fundamental plane (FP) relation (Djorgovski & Davis 1987). Two large recent samples of peculiar velocities determined with such methods are the Cosmicflows project for the TF relation (Kourkchi et al. 2022) and SDSS-PV for the FP relation (Howlett et al. 2022). Type Ia supernovae (SNe Ia) have also been considered as another distance indicator to determine the peculiar velocities Howlett et al. (2017), Huterer et al. (2017), Scolnic et al. (2019), Kim & Linder (2020), Carreres et al. (2023).

Peculiar velocities can be combined with measurements of redshift-space distortion (RSD) of galaxies in the same volume, allowing us to reach stronger constraints on the growth rate. In particular, this can be done with compressed two-point statistics such as the density and momentum correlation function, power spectra, or the average pair-wise velocities (Ferreira et al. 1999; Dupuy et al. 2019; Howlett 2019; Turner et al. 2022; Qin et al. 2025). An alternative class of methodologies infer cosmological parameters directly from the velocity and density fields, without compression. The density-velocity comparison method (Springob et al. 2014; Carrick et al. 2015; Boruah et al. 2020; Said et al. 2020; Qin et al. 2023; Boubel et al. 2024) compares the observed velocity field to a predicted model based on an observed galaxy density field, typically using reconstruction techniques. A more complex method often referred to as forward modeling, or simulation-based inference (Boruah et al. 2021; Valade et al. 2022a, b; Pfeifer et al. 2023), consists of evolving initial conditions of density and velocity fields using some theoretical model, generally through perturbation theory or simplified simulations, to fit observations.

The focus of this work is on another type of commonly used method: the likelihood-based field-level estimator (Johnson et al. 2014; Howlett et al. 2017; Adams & Blake 2017, 2020; Huterer et al. 2017; Lai et al. 2022; Carreres et al. 2023). It consists of calculating the theoretical correlations for every pair of positions in the considered field (density or velocities), as a function of cosmological and nuisance parameters. Those parameters are varied to maximize the likelihood of the model given the data, which is assumed to be drawn from a multivariate Gaussian distribution. This maximization can be performed in a statistical or Bayesian way. Note that in this method, correlations are not compressed into binned statistics such as the correlation function or the power spectrum, and the theoretical covariance matrix corresponds to the field covariance.

Compared to traditional two-point statistics, the likelihood-based estimator allows for simultaneous fitting of the analysis parameters used to derive the fields, along with cosmological and nuisance parameters, while catching potential degeneracies between parameters. To give a more practical example, when considering SN Ia peculiar velocities, the velocities are estimated from Hubble diagram residuals, which depend on standardization parameters. The parameter of interest is the growth rate of structures, but some nuisance parameters are needed to perform the fit. In addition, the likelihood-based method allows one to capture two-point correlations in the field without compression, maximizing the potential extracted cosmological information. This method also allows one to add elements of the model predictions not described by Gaussian statistics, such as non-Gaussian likelihood corrections. The likelihood-based method fits a small number of parameters compared, for instance, to forward-modeling techniques in which the field data points considered are also free parameters, increasing the computational cost and complexity of the inference. However, one disadvantage of the likelihood-based method is its expensiveness in computation time and memory, as the likelihood calculation generally requires a covariance matrix inversion, which scales as N², where N is the number of objects or mesh cells. The latter can scale up to approximately 10 000 elements compared to compressed statistics, where the data vector is typically smaller than O(100) elements. Furthermore, in previous studies, the likelihood-based estimator generally computes the field theoretical correlations with linear models, and its expansion to a nonlinear model yields a higher mathematical complexity than compressed statistics and is generally computationally expensive.

To address the aforementioned issues of the likelihood-based field-level estimator, we developed a generalized framework that aims to reproduce the past works models (Adams & Blake 2017, 2020; Lai et al. 2022; Carreres et al. 2023), with or without wide-angle modeling, and to extend them to more complex field covariance models. The mathematical foundations of this framework allows one to treat all those models in a consistent way, with adapted algorithmic optimization. This framework has been implemented in the Python package flip¹, which is a generalized extension of an early version used for likelihood-based inference of SNe Ia peculiar velocities in Carreres et al. (2023). The purpose of this paper is to present the mathematical formalism used in this package and all the applications currently included.

This work paves the way for likelihood-based field-level method improvements, such as nonlinear power spectra models, the direct inclusion in the analysis of observational effects at the field level, or the unified treatment of several velocity fields (FP, TF, SNe Ia) with galaxy density field estimates.

This article is organized as follows. In Section 2, we detail our generalized framework for field-level covariance computation, in a wide-angle or plane-parallel configuration. Section 3 gives an overview of the likelihood-based field-level method, along with previous and new covariance models developed in this study. Section 4 shows the different validation tests performed on the flip algorithm. In Section 5, we apply our covariance calculation to create a survey-dependent Fisher forecasting tool. Finally, Section 6 discusses future improvements of the flip package.

2. Generalized theoretical covariance framework

When performing a likelihood-based inference from two fields, denoted as a and b, which can be for example radial velocities, v, galaxy number density, δ, velocity divergence, θ, or logarithmic distance ratios, η, we assume that the fields can be described by perturbations on a Gaussian random field that depends on a theoretical covariance denoted as C_ab. In the method presented here, this covariance matrix is estimated in an analytical way by computing the theoretical correlation between each considered field. Note that the covariance terms here refer to the covariance at the field level, and is different from the covariance matrix used to fit correlation functions and power spectra. The main objective of this section is to derive a framework for fast computation of theoretical correlations of velocities and densities.

2.1. Coordinate definition

The coordinates used to develop our framework are shown in Fig. 1. We consider two points of a given field (e.g., density or velocity) with distances to the observer represented by the vectors r₁ and r₂, respectively. Their separation vector is r=r₁−r₂ and their angular separation is represented by the angle α. We also define the α-angle bisector d and the angle ϕ between d and r.

Fig. 1. Schematic representation of the two field elements, a1 and b2, for which we want to compute the theoretical correlation. Those two fields can be the peculiar velocity of a considered galaxy or group of galaxies, or the galaxy density field itself. The definition of the vector, d, and consequently the angle, ϕ, depends on the chosen wide-angle definition.

In this coordinate definition, we are considering the general case where the individual line-of-sight r₁ and r₂ are not parallel. Calculating the correlation in the non-parallel case is often referred to as wide-angle modeling. The plane-parallel approximation is when the vectors r₁, r₂, and d are all considered to be parallel. We note that in the wide-angle case, there exist several potential definitions of the d vector, considered as the line-of-sight reference for correlation calculation (Beutler et al. 2019). Changing the definition of this vector can have an impact on the modeling of correlation calculation (see Castorina & White 2020). The flip package implements the bisector (the one used in this article), midpoint and endpoint definitions of d:

$$\begin{array}{cc}\hfill {\text{d}}_{\mathrm{bisector}}& =\frac{r_1{r}_2}{r_1+r_2}\left(\frac{{\text{r}}_1}{r_1}+\frac{{\text{r}}_2}{r_2}\right)\hfill \\ \hfill {\text{d}}_{\mathrm{midpoint}}& =\frac{{\text{r}}_1+{\text{r}}_2}{2}\hfill \\ \hfill {\text{d}}_{\mathrm{endpoint}}& ={\text{r}}_1.\hfill \end{array}$$(1)

We computed the covariance of the fields a and b, respectively evaluated at positions r₁ and r₂. Note that this formalism is kept general on purpose and can be extended to any considered correlated field described by some power spectrum model. We used the following Fourier transform convention:

$$\mathrm{a}(\text{r})=\frac{1}{(2\pi {)}^3}{\displaystyle {\int }_{\text{k}}}{d}^3\text{k}\hspace{0.5em}\tilde{\mathrm{a}}(\text{k})e^{i\text{k}·\text{r}}.$$(2)

In addition we define the individual line-of-sight angle cosine ${\mu }_i=\hat{\text{k}}·{\hat{\text{r}}}_i$ where $\hat{\text{k}}=\text{k}/k$ and ${\hat{\text{r}}}_i={\text{r}}_i/r_i$ with k and r_i are the norms of the vectors k and r_i, respectively. We also define the cosine angle $\mu =\hat{\text{k}}·\hat{\text{d}}=cos(\varphi )$. The power spectrum is defined by

$$\langle \tilde{\mathrm{a}}(\text{k}){\tilde{\mathrm{b}}}^{*}({\text{k}}^{\prime })\rangle =(2\pi {)}^3{P}_{\mathrm{ab}}(k,{\mu }_1,{\mu }_2){\delta }_D(\text{k}-{\text{k}}^{\prime }),$$(3)

where δ_D is the Dirac delta function. The corresponding correlation is then defined as the Fourier transform of the power spectrum as:

$$C_{\mathrm{ab}}({\text{r}}_1,{\text{r}}_2)=\frac{1}{(2\pi {)}^3}{\displaystyle {\int }_{\text{k}}}{d}^3\text{k}{P}_{\mathrm{ab}}(k,{\mu }_1,{\mu }_2)e^{i\text{k}·\text{r}}.$$(4)

2.2. Wide-angle framework

In the general wide-angle formalism, we consider that for every pair of fields a and b, we can decompose the power spectrum model as

$$P_{\mathrm{ab}}(k,{\mu }_1,{\mu }_2)={\displaystyle \sum _n}{\mathit{w}}_{\mathrm{ab},n}{F}_{\mathrm{ab},n}(k,{\mu }_1,{\mu }_2){\mathcal{P}}_{\mathrm{ab},n}(k),$$(5)

where n stands for the index of the term in the power spectrum model.

The logic of this decomposition is to have a linear combination of geometrical terms F_ab,n(k,μ₁,μ₂) that contain the angular information to account for wide-angle effects, and isotropic power spectrum terms 𝒫_ab,n(k) that can be computed by a Boltzmann solver software. The coefficient of the linear decomposition w_ab,n are the parameters we want to fit in the likelihood-based framework.

The main interest of this decomposition is to compute separately each covariance term just once and only vary w_ab,n when maximizing the likelihood. Furthermore, integrals of the geometrical part (F_ab,n(k,μ₁,μ₂)) can be performed analytically while the power spectrum part (𝒫_ab,n(k)) can be integrated numerically in an algorithmically optimized way. Note that this decomposition might not be directly applicable for more complex nonlinear models, but this is beyond the scope of this paper.

The corresponding covariance matrix given in Equation (4) can then be written

$$C_{\mathrm{ab}}({\text{r}}_1,{\text{r}}_2)={\displaystyle \sum _n}{\mathit{w}}_{\mathrm{ab},n}{C}_{\mathrm{ab},n}({\text{r}}_1,{\text{r}}_2),$$(6)

where C_ab,n(r₁,r₂) can be expressed as an integral over the spherical coordinates of the vector k, defined in Eq. (A.1), as

$$\begin{array}{cc}\hfill C_{\mathrm{ab},n}({\text{r}}_1,{\text{r}}_2)=& \frac{1}{(2\pi {)}^3}{\displaystyle {\int }_{\text{k}}}{d}^3\text{k}{F}_{\mathrm{ab},n}(k,{\mu }_1,{\mu }_2){\mathcal{P}}_{\mathrm{ab},n}(k)e^{i\text{k}·\text{r}}\hfill \\ \hfill =& {\displaystyle {\int }_0^{\infty }}\frac{k^2\mathit{dk}}{(2\pi {)}^3}{\mathcal{P}}_{\mathrm{ab},n}(k)\hfill \\ \hfill & \times {\displaystyle {\int }_{\mathrm{\Omega }}}d\mathrm{\Omega }{F}_{\mathrm{ab},n}(k,{\mu }_1,{\mu }_2)e^{i\text{k}·\text{r}}.\hfill \end{array}$$(7)

Using the Legendre polynomial expansion (Eq. A.3) on F_ab,n(k,μ₁,μ₂) with respect to μ₁ and μ₂ we get

$$F_{\mathrm{ab},n}(k,{\mu }_1,{\mu }_2)={\displaystyle \sum _{{\ell }_1,{\ell }_2}}\frac{(2{\ell }_1+1)(2{\ell }_2+1)}{4}{L}_{{\ell }_1}({\mu }_1)L_{{\ell }_2}({\mu }_2)$$(8)

$$\times {\displaystyle {\int }_{-1}^1}{\displaystyle {\int }_{-1}^1}d{\mu }_1^{\prime }d{\mu }_2^{\prime }{F}_{\mathrm{ab},n}(k,{\mu }_1^{\prime },{\mu }_2^{\prime })L_{{\ell }_1}({\mu }_1^{\prime })L_{{\ell }_2}({\mu }_2^{\prime }),$$(9)

where the L_ℓ are the Legendre polynomials. For conciseness, we define a wavenumber term $M_{\mathrm{ab},n}^{{\ell }_1,{\ell }_2}(k)$ that contains the analytical integration of the F_ab,n(k,μ₁,μ₂) term over μ₁ and μ₂ such that

$$M_{\mathrm{ab},n}^{{\ell }_1,{\ell }_2}(k)=\frac{1}{4}{\displaystyle {\int }_{-1}^1}{\displaystyle {\int }_{-1}^1}d{\mu }_1^{\prime }d{\mu }_2^{\prime }{F}_{\mathrm{ab},n}(k,{\mu }_1^{\prime },{\mu }_2^{\prime })L_{{\ell }_1}({\mu }_1^{\prime })L_{{\ell }_2}({\mu }_2^{\prime }).$$(10)

Combining Eqs. (7), (8) and using the plane-wave expansion (Eq. A.5) to decompose the e^ik·r term, we obtain

$$\begin{array}{cc}\hfill C_{\mathrm{ab},n}({\text{r}}_1,{\text{r}}_2)=& {\displaystyle {\int }_0^{\infty }}\frac{k^2\mathit{dk}}{(2\pi {)}^3}{\mathcal{P}}_{\mathrm{ab},n}(k)\hfill \\ \hfill & \times {\displaystyle \sum _{\ell ,{\ell }_1,{\ell }_2}}{i}^{\ell }{j}_{\ell }(\mathit{kr})(2\ell +1)(2{\ell }_1+1)(2{\ell }_2+1)M_{\mathrm{ab},n}^{{\ell }_1,{\ell }_2}(k)\hfill \\ \hfill & \times {\displaystyle {\int }_{\mathrm{\Omega }}}d\mathrm{\Omega }{L}_{{\ell }_{}}(\hat{\text{k}}·\hat{\text{r}})L_{{\ell }_1}({\mu }_1)L_{{\ell }_2}({\mu }_2).\hfill \end{array}$$(11)

Using Eq. (A.9) to express the three Legendre polynomial product, we obtain

$$\begin{array}{cc}\hfill C_{\mathrm{ab},n}({\text{r}}_1,{\text{r}}_2)=& {\displaystyle \sum _{\ell ,{\ell }_1,{\ell }_2}}(4\pi {)}^2{\displaystyle \sum _{m,m_1,m_2}}{G}_{\ell ,{\ell }_1,{\ell }_2}^{m,m_1,m_2}{Y}_{{\ell }_{}^{}{m}_{}^{}}^{*}(\hat{\text{r}})Y_{{\ell }_1^{}{m}_1^{}}^{*}(\widehat{{\text{r}}_1})Y_{{\ell }_2^{}{m}_2^{}}^{*}(\widehat{{\text{r}}_2})\hfill \\ \hfill & \times i^{\ell }{\displaystyle {\int }_0^{\infty }}\frac{k^2\mathit{dk}}{2{\pi }^2}{\mathcal{P}}_{\mathrm{ab},n}(k)M_{\mathrm{ab},n}^{{\ell }_1,{\ell }_2}(k)j_{\ell }(\mathit{kr}).\hfill \end{array}$$(12)

For conciseness, we write the coefficient $N_{\mathrm{ab},\ell }^{{\ell }_1,{\ell }_2}$ containing the linear combination of spherical harmonics such that

$$N_{\mathrm{ab},\ell }^{{\ell }_1,{\ell }_2}(\hat{\text{r}},{\hat{\text{r}}}_1,{\hat{\text{r}}}_2)=(4\pi {)}^2{\displaystyle \sum _{m,m_1,m_2}}{G}_{\ell ,{\ell }_1,{\ell }_2}^{m,m_1,m_2}{Y}_{{\ell }_{}^{}{m}_{}^{}}^{*}(\hat{\text{r}})Y_{{\ell }_1^{}{m}_1^{}}^{*}(\widehat{{\text{r}}_1})Y_{{\ell }_2^{}{m}_2^{}}^{*}(\widehat{{\text{r}}_2}).$$(13)

We can explicitly express the dependency of the spherical harmonic terms inside the $N_{\mathrm{ab},\ell }^{{\ell }_1,{\ell }_2}$ terms considering the definition of angles in Figure 1:

$$\begin{array}{cc}\hfill N_{\mathrm{ab},\ell }^{{\ell }_1,{\ell }_2}(\hat{\text{r}},{\hat{\text{r}}}_1,{\hat{\text{r}}}_2)=(4\pi {)}^2& {\displaystyle \sum _{m,m_1,m_2}}{G}_{\ell ,{\ell }_1,{\ell }_2}^{m,m_1,m_2}{Y}_{{\ell }_{}^{}{m}_{}^{}}^{*}(\pi -\varphi ,0)\hfill \\ \hfill & \times Y_{{\ell }_1^{}{m}_1^{}}^{*}(\alpha /2,0)Y_{{\ell }_2^{}{m}_2^{}}^{*}(\alpha /2,\pi ).\hfill \end{array}$$(14)

The choice regarding the π terms in the orientation of the angles depends entirely on the definition of r. Furthermore, the value of the ϕ angle depends on the definition of the d vector, and can have an impact on field covariance modeling.

We define the zeroth-order Hankel transform ℋ_ℓ in a more cosmologically relevant form following Karamanis & Beutler (2021) as

$${\mathcal{H}}_{\ell }\left[f(k)\right](r)=i^{\ell }{\displaystyle {\int }_0^{\infty }}\frac{k^2\mathit{dk}}{2{\pi }^2}{j}_{\ell }(\mathit{kr})f(k).$$(15)

Finally, the individual covariance terms can be written:

$$C_{\mathrm{ab},n}({\text{r}}_1,{\text{r}}_2)={\displaystyle \sum _{\ell ,{\ell }_1,{\ell }_2}}{N}_{\mathrm{ab},\ell }^{{\ell }_1,{\ell }_2}(\varphi ,\alpha ){\mathcal{H}}_{\ell }\left[{\mathcal{P}}_{\mathrm{ab},n}(k)M_{\mathrm{ab},n}^{{\ell }_1,{\ell }_2}(k)\right](r).$$(16)

This equation is in an algorithmically optimized form for a fast computation of the field covariance matrix. The integration of the $N_{\mathrm{ab},\ell }^{{\ell }_1,{\ell }_2}$ and $M_{\mathrm{ab},n}^{{\ell }_1,{\ell }_2}$ terms can be performed analytically, and a numerical integration is performed over the wavenumber k. The latter can easily be accelerated as it is in the form of a Hankel transform.

In practice, for N field points, the integral of Eq. (4) has to be computed N(N+1)/2 times if we include the diagonal term. This process is computationally intensive, particularly in the context of the new generation of surveys (e.g., the Legacy Survey of Space and Time Rubin-LSST; Ivezić et al. 2018 or the Dark Energy Spectroscopic Instrument DESI; DESI Collaboration 2016a, b; Martini et al. 2018). The flip software, with the formalism developed here, proposes an efficient way to compute the covariance matrix using Hankel transforms and parallelized processes.

In our derivation, we did not consider the redshift at which the model was calculated. For small surveys concentrated on a limited redshift range, accounting for this might not be necessary, but this is not the case for the new aforementioned surveys. In particular, the parameters w_ab,n can depend on the redshift of the two considered points. Furthermore, the power spectra, ${\mathcal{P}}_{\mathrm{ab},n}$, computed with Boltzmann solvers also depend on redshift. Assuming that the power spectrum redshift dependency can be factorized, we developed an option to account for different redshift between fields a and b and provide details in Appendix B.

We have only considered a covariance calculation between scalar fields, either between components of a vector field or a scalar field itself. Our formalism can be extended to correlations between vector fields, which means computing covariances between all components of each vector, yielding a covariance tensor. This derivation is given in Appendix C.

2.3. Plane-parallel approximation

In the case in which the angles between the two points of the considered field are small, we can safely use the plane-parallel approximation and simplify the modeling of the field covariance matrix. Here, we show how this approximation is implemented in flip. The principal purpose of this implementation is to have reference models and to compare them with previous implementations (Adams & Blake 2017, 2020).

In the plane-parallel approximation it is assumed that the two galaxies are on the same plane at the same distance, d. Then we have ${\mu }_1={\mu }_2=\mu =\hat{\text{k}}·\hat{\text{d}}$. The covariance of the Eq. (7) is now written as

$$\begin{array}{cc}\hfill C_{\mathrm{ab},n}({\text{r}}_1,{\text{r}}_2)=& {\displaystyle {\int }_0^{\infty }}\frac{k^2\mathit{dk}}{(2\pi {)}^3}{\mathcal{P}}_{\mathrm{ab},n}(k)\hfill \\ \hfill & \times {\displaystyle {\int }_{\mathrm{\Omega }}}d\mathrm{\Omega }{F}_{\mathrm{ab},n}(k,\mu )e^{i\text{k}·\text{r}}.\hfill \end{array}$$(17)

Using the Legendre polynomial expansion (Eq. A.3), we obtain

$$F_{\mathrm{ab},n}(k,\mu )={\displaystyle \sum _{{\ell }^{\prime }}}\frac{(2{\ell }^{\prime }+1)}{2}{L}_{{\ell }^{\prime }}(\mu ){\displaystyle {\int }_{-1}^1}d{\mu }^{\prime }{F}_{\mathrm{ab},n}(k,{\mu }^{\prime })L_{{\ell }^{\prime }}({\mu }^{\prime }).$$(18)

Similarly to the wide-angle case, we define the coefficients $M_{\mathrm{ab},n}^{{\ell }^{\prime }}(k)$ as

$$M_{\mathrm{ab},n}^{{\ell }^{\prime }}(k)=\frac{1}{2}{\displaystyle {\int }_{-1}^1}d{\mu }^{\prime }{F}_{\mathrm{ab},n}(k,{\mu }^{\prime })L_{{\ell }^{\prime }}({\mu }^{\prime }).$$(19)

In addition, we use the plane-wave expansion (Eq. A.5) such that the Eq. (17) becomes

$$\begin{array}{cc}\hfill C_{\mathrm{ab},n}({\text{r}}_1,{\text{r}}_2)=& {\displaystyle {\int }_0^{\infty }}\frac{k^2\mathit{dk}}{(2\pi {)}^3}{\mathcal{P}}_{\mathrm{ab},n}(k)\hfill \\ \hfill & \times {\displaystyle \sum _{\ell ,{\ell }^{\prime }}}{i}^{\ell }{j}_{\ell }(\mathit{kr})(2\ell +1)(2{\ell }^{\prime }+1)M_{\mathrm{ab},n}^{{\ell }^{\prime }}(k)\hfill \\ \hfill & \times {\displaystyle {\int }_{\mathrm{\Omega }}}d\mathrm{\Omega }{L}_{{\ell }_{}}(\hat{\text{k}}·\hat{\text{r}})L_{{\ell }^{\prime }}(\mu ).\hfill \end{array}$$(20)

Using the result of Eq. (A.12) for the angular integral in Eq. (20), considering as a reminder that $\mu =\hat{\text{k}}·\hat{\text{d}}$, we get

$$\begin{array}{cc}\hfill C_{\mathrm{ab},n}({\text{r}}_1,{\text{r}}_2)& ={\displaystyle \sum _{\ell }}4\pi {\displaystyle \sum _{m=-\ell }^{\ell }}{Y}_{{\ell }_{}^{}{m}_{}^{}}(\hat{\text{d}})Y_{{\ell }_{}^{}{m}_{}^{}}^{*}(\hat{\text{r}})\hfill \\ \hfill & \times {\displaystyle {\int }_0^{\infty }}{i}^{\ell }\frac{k^2\mathit{dk}}{2{\pi }^2}{\mathcal{P}}_{\mathrm{ab},n}(k)M_{\mathrm{ab},n}^{\ell }(k)j_{\ell }(\mathit{kr}).\hfill \end{array}$$(21)

We define the geometrical term $N_{\mathrm{ab},i}^{\ell }(\hat{\text{r}},\hat{\text{d}})$, and simplify it using the spherical harmonic addition theorem (Eq. A.6):

$$N_{\mathrm{ab},\ell }(\hat{\text{r}},\hat{\text{d}})=4\pi {\displaystyle \sum _{m=-\ell }^{\ell }}{Y}_{{\ell }_{}^{}{m}_{}^{}}(\hat{\text{d}})Y_{{\ell }_{}^{}{m}_{}^{}}^{*}(\hat{\text{r}})=(2\ell +1)L_{\ell }(cos\left(\pi -\varphi \right)).$$(22)

In this plane-parallel case, those geometric terms are dependent only on the ϕ angle, and thus also on the definition of the referential line-of-sight d. The final expression for the covariance term in the plane-parallel model is given by

$$C_{\mathrm{ab},n}({\text{r}}_1,{\text{r}}_2)={\displaystyle \sum _{\ell }}{N}_{\mathrm{ab},\ell }(\varphi ){\mathcal{H}}_{\ell }\left[{\mathcal{P}}_{\mathrm{ab},n}(k)M_{\mathrm{ab},n}^{\ell }(k)\right](r).$$(23)

This equation corresponds to a similar form to the wide-angle case (Equation 16) with simpler analytical integrations. The properties concerning the algorithmic efficiency are still valid. Our framework gives the advantage to include wide-angle effects, which can be very computationally intensive and mathematically complex, in the same formalism.

3. Likelihood-based field-level inference

The first direct application of the previously derived framework is the inference of cosmological parameters with a likelihood maximization. This method allows one to make an inference on a set of parameters that we note Θ by maximizing the likelihood whose general form is $\mathcal{L}\left[x(\mathrm{\Theta }),C_{\mathit{xx}}(\mathrm{\Theta }),C_{\mathrm{obs}}(\mathrm{\Theta })\right]$, where the field x, the field covariance C_xx, and the observational covariance C_obs depend on a set of parameters Θ. For example, the latter matrix can contain the field error bar or be a full covariance comprising additional field correlations caused by spurious instrumental effects.

The schematic implementation of this method in the flip package is represented in Fig. 2. We use the input data to generate the field covariance components C_ab,n(r₁,r₂), following the derivation from the Section 2, and to create a data vector x formatted in a class which also contains the observational covariance. The likelihood is built from the theoretical covariance and the data vector. The parameters are fitted from this likelihood either by best-fit minimizer or with an MCMC sampling.

Fig. 2. Schematical implementation of the different modules in the flip package. The arrow represents data flows between two modules. The link in gray between vector construction and the covariance calculation represents an alternative way to compute the covariance matrix directly from the data vector object.

3.1. Covariance model developed for inference

To better assess the performance of our generalized field covariance framework, we first aimed at reproducing the latest models available in the literature (Adams & Blake 2017, 2020; Lai et al. 2022; Carreres et al. 2023). In the second step, we used lessons learned to develop a new covariance model for densities and velocities, accounting for wide-angle effects and with superior performance.

As an illustrative example, we start by decomposing a model used in literature for velocity correlations with wide-angle modeling in Adams & Blake (2017), Lai et al. (2022), Carreres et al. (2023). Their velocity power spectrum model contains one term:

$$P_{\mathrm{vv}}=(\mathit{aHf}{\sigma }_8{)}^2\frac{{\mu }_1{\mu }_2}{k^2}{P}_{\theta \theta }(k)D_u^2(k,{\sigma }_u),$$(24)

where the a is the scale factor, H the Hubble parameter, f the growth rate of structure, σ₈ the amplitude of the matter perturbations in spheres of 8 h⁻¹ Mpc comoving radius, P_θθ is the velocity dispersion power spectrum normalized by fiducial ${\sigma }_8^2$ value, and D_u is an empirical damping function often used in peculiar velocity studies to model the effect of RSD on the position of velocities itself (see e.g., Koda et al. 2014). In the flip framework of Carreres et al. (2023) model, we decompose the Equation (5) over one covariance term:

$${\mathit{w}}_{\mathrm{vv},0}=(f{\sigma }_8{)}^2,$$(25)

$$F_{\mathrm{vv},0}(k,{\mu }_1,{\mu }_2)=(\mathit{aH}{)}^2\frac{{\mu }_1{\mu }_2}{k^2},$$(26)

$${\mathcal{P}}_{\mathrm{vv},0}=P_{\theta \theta }(k)D_u^2(k,{\sigma }_u).$$(27)

A simple sympy² (Meurer et al. 2017) computation of the $N_{\mathrm{ab},\ell }^{{\ell }_1,{\ell }_2}$ and $M_{\mathrm{ab},n}^{{\ell }_1,{\ell }_2}$ terms (in Equations 10 and 13) states that only following terms are nonzero:

$$M_{\mathrm{vv},n=0}^{{\ell }_1=1,{\ell }_2=1}=\frac{(\mathit{aH}{)}^2}{9k^2},$$(28)

$$N_{\mathrm{vv},l=0}^{{\ell }_1=1,{\ell }_2=1}=3cos(\alpha ),$$(29)

$$N_{\mathrm{vv},l=2}^{{\ell }_1=1,{\ell }_2=1}=\frac{(9cos2\varphi +3cos\alpha )}{2k^2}·$$(30)

Adding up the different terms to compute the field covariance matrix, we obtain

$$\begin{array}{cc}\hfill C_{\mathrm{vv}}({\text{r}}_1,{\text{r}}_2)=& \frac{(\mathit{aHf}{\sigma }_8{)}^2}{2{\pi }^2}{\displaystyle {\int }_k}\mathit{dk}{P}_{\theta \theta }{D}_u^2(k,{\sigma }_u)\hfill \\ \hfill & \left[j_0(\mathit{kr})\frac{cos(\alpha )}{3}-j_2(\mathit{kr})\frac{(3cos2\varphi +cos\alpha )}{6}\right].\hfill \end{array}$$(31)

Following the demonstration of Lai et al. (2022) in Appendix B.3, we obtain the same formula in previous articles (Adams & Blake 2017; Lai et al. 2022; Carreres et al. 2023). Note that the demonstration uses the bisector theorem, meaning that the mathematical equivalency is only valid when the line-of-sight reference is the bisector. The flip covariance framework is unnecessarily complex for simple linear velocity modeling. However, this is not the case for multiple-term power spectrum models with complex forms, which is the case for density models, especially considering nonlinear corrections.

Table D.1 lists the terms used in the flip framework to reproduce all the previously designed models and a new model we developed for this study. Adams & Blake (2017) developed a model (AB17 hereafter) aiming at linking density and velocity information with a density model without RSD. Thanks to this simplification, they were able to account for wide-angles. The first inclusion of RSD in the density terms was realized in the model presented in Adams & Blake (2020), named AB20. As a simplification, they choose to perform the integration for a plane parallel model, i.e., taking the RSD anisotropic power spectra for vv, gv, and vv terms. This assumption breaks for large-area surveys at low redshift, and it is thus insufficient to properly model the clustering for the current and future generations of survey. However, this plane-parallel model is adapted to test the wide-angle models in the limit of very high redshift. For the case of Carreres et al. (2023), named C23, which only contains a velocity modeling vv, the terms correspond to the wide-angle vv of AB17. The authors in Lai et al. (2022) propose a wide-angle model, which we call L22, by performing a Taylor expansion on the RSD finger-of-god (FoG) small-scale term:

$$e^{-\frac{{\left(k\mu {\sigma }_{\mathrm{g}}\right)}^2}{2}}={\displaystyle \sum _{p=0}^{\infty }}\frac{(-1{)}^p{\left(k{\sigma }_{\mathrm{g}}\right)}^{2p}}{2^{p}p!}{\mu }^{2p}.$$(32)

Since this decomposition is performed two times for the gg terms, one time for the gv terms, and calculated to the p = 4 order, it yields many covariance terms. In Table D.1, we have decomposed this model in the flip formalism by summing the density-density (gg) terms which have the same value for the sum m=p+q. For the gv cross-correlation terms, we have the same formalism as L22, renaming the p index by m for consistency. This model allows ont to directly fit the FoG parameter (σ_g) and manages to express all the covariance computations in terms of Hankel transforms, thus speeding up the calculation. The main weaknesses of this model are that the Taylor expansion is computed up to a certain fixed order and that it is based on the assumption that the product kμσ_g is small compared to unity. Knowing that σ_g values are generally considered ranging [1−10] h⁻¹ Mpc, it means that the model is only valid for low values of k (typically k<0.1 h Mpc⁻¹). It is limiting, especially for an extension to nonlinear models which needs small-scale modes for integration stability, directly at the power spectrum level (${\mathcal{P}}_{\mathrm{ab},n}$) or in the full power spectrum model (P_ab(k,μ₁,μ₂)). The authors of the L22 model decided to separate small-scale modeling (0.2<k<1.0 h Mpc⁻¹) with the rest of the scales considered and introduce a dedicated nuisance parameter for the amplitude of the covariance in the small scales.

In this study, we created the new covariance model, denoted as RC25, which frees itself from the Taylor expansion approximation while keeping the wide-angle modeling. Removing the Taylor expansion allows one to extend modeling to larger wavenumbers (k>0.1 h Mpc⁻¹), which is necessary for the use of nonlinear power spectrum models. As shown in Table D.1, we obtain this model by reverting the L22 model Taylor-expansion, or equivalently by including wide-angle terms in the AB20 model. Contrary to the L22 model, RC25 cannot account for σ_g nuisance parameter directly in the linear decomposition of the power spectrum model. We overcome this issue by pre-computing covariances for several values of σ_g and interpolating during the inference. This interpolation procedure will be detailed in Section 3.2.

The flip.covariance module of the Figure 2 contains for all the models previously described, the terms $N_{\mathrm{ab},\ell }^{{\ell }_1,{\ell }_2}(\varphi ,\alpha )$ and $M_{\mathrm{ab},n}^{{\ell }_1,{\ell }_2}(k)$ terms (respectively, N_ab,ℓ(ϕ) and $M_{\mathrm{ab},n}^{\ell }(k)$ for plane-parallel models). The integration inside those terms is performed analytically with the symbolic library sympy. The derived functions are stored inside the flip.covariance module by symbolic code generation.

When a covariance model is generated following the Equation (16) (or 23 in the plane-parallel case), the geometrical terms $N_{\mathrm{ab},\ell }^{{\ell }_1,{\ell }_2}(\varphi ,\alpha )$ (N_ab,ℓ(ϕ)) are computed for all the field pairs considered. In a second step, the term ${\mathcal{H}}_{\ell }\left[{\mathcal{P}}_{\mathrm{ab},n}(k)M_{\mathrm{ab},n}^{{\ell }_1,{\ell }_2}(k)\right]$ is computed numerically by first expressing $M_{\mathrm{ab},n}^{{\ell }_1,{\ell }_2}(k)$ (or $M_{\mathrm{ab},n}^{\ell }(k)$ in the plane-parallel case) for the considered pair and by computing the Hankel transform ℋ_ℓ. This last step is a numerical integration performed over the wavenumber k and is accelerated by employing the FFTLog algorithm (Talman 1978). For a given data field vector x of size N, we must compute N(N+1)/2 times this equation. In order to speed up this process, we distribute the computation of the covariance with multi-threading by separating the covariance into arrays.

When the correlation is computed between two different fields, for example the radial velocity, v, and galaxy density, δ_g, scalar field, we only compute C_gv and deduce C_vg. Adams & Blake (2017) (Appendix A.3) showed that a minus sign appears between C_vg and C_gv when employing their definitions. However, the equation for cross terms is not symmetric because only odd Legendre polynomials are nonzero, so changing the r vector between gv and vg correlations introduces another minus sign. Consequently, we can simply take C_vg as the transpose of C_gv.

Some numerical considerations must be considered to compute all the covariance models shown in Table D.1. For models including RSD FoG terms without Taylor expansion such as RC25 and AB20, the $M_{\mathrm{ab},n}^{{\ell }_1,{\ell }_2}(k)$ or $M_{\mathrm{ab},n}^{\ell }(k)$ terms can lead to numerical instabilities due to linear combination of large floats. A numerical correction, detailed in the Appendix E, allows us to solve this issue. Additionally, the use of the Hankel transform can introduce numerical instabilities. The latter can be caused by considering the smaller separations between the pair for which we want to compute the covariance. In the special case of simulations, other instabilities can occur at large scales when the numerical integration over wavenumber is performed on a range not covered by the considered simulation. We solve this issue by introducing a regularization term for low wavenumbers, also detailed in the Appendix E.

3.2. Vector and Likelihood estimators

As shown in Fig. 2, the flip package includes modules to perform cosmological parameter fits. The input data are handled by different Python classes implemented in the flip.data_vector module. This module includes various Python classes to handle density and velocity survey data.

For peculiar velocities, several estimators are implemented to cover different types of surveys, such as TF, FP, and SNe Ia. For simulation, where true peculiar velocity is accessible, a data vector mode with the direct velocities is implemented, with an option for grouping them when they are near each other. For TF and FP studies, the flip software contains a commonly used velocity estimator for a galaxy i based on the logarithmic distance ratio:

$$v_i=ln(10)\frac{D\left(z_{\mathrm{obs},i}\right)H\left(z_{\mathrm{obs},i}\right)}{\left(1+z_{\mathrm{obs},i}\right)}{\eta }_i,$$(33)

where ${\eta }_i=log\left[D\left(z_{\mathrm{obs},i}\right)/D\left(z_{\cos ,i}\right)\right]$ is the logarithmic distance ratio of the galaxy i, z_obs,i and z_cos,i are, respectively, the observed and cosmological redshifts, and D is the comoving radial distance expressed in h⁻¹ Mpc.

The flip software also contains a method to go from Hubble diagram residuals Δμ_i to peculiar velocities, which can be used for any distance indicator (TF, FP, SNe Ia). Several methods for transforming Hubble diagram residuals to peculiar velocities are implemented; they take the form

$$v_i=-\frac{cln(10)}{5}J(p_i)\mathrm{\Delta }{\mu }_i,$$(34)

where p are the parameters needed to estimate the velocities. The currently implemented estimators are given in Table 1.

The uncertainty of the radial peculiar velocity estimate σ_v,i is computed by propagating the estimated error in the observable used. If velocities are directly used, flip needs an estimate of the error on those velocities as an input. For the logarithmic distance (respectively, Hubble residuals), the error bar on velocities is computed by replacing η_i (resp. Δμ_i) on Equation (33) (resp. 34) by the error on logarithmic distance σ_η,i (resp. Hubble residual σ_Δμ,i).

For the special case of SNe Ia, we express explicitly the Hubble residuals as a function of light-curve fitting parameters. Those residuals are expressed using the Tripp relation (Tripp 1998) and given in greater detail in Carreres et al. (2023):

$$\begin{array}{cc}\hfill \mathrm{\Delta }{\mu }_i({\mathrm{\Theta }}_{\mathrm{HD}})& ={\mu }_{\mathrm{obs},i}({\mathrm{\Theta }}_{\mathrm{HD}})-{\mu }_{\mathrm{model},i}\hfill \\ \hfill & =m_{B,i}-M_0+\alpha x_{1,i}-\beta c_i\hfill \\ \hfill & -\left(5log\left[(1+z_{\mathrm{obs},i})D(z_{\mathrm{obs},i})\right]+25\right),\hfill \end{array}$$(35)

where m_B,i is the apparent magnitude of SN Ia i in the B-band, M₀ is the absolute magnitude of SNe Ia which can be seen as an offset of the SN Ia Hubble diagram, α and β correct for correlations of Hubble residuals with SN Ia stretch x_1,i and color c_i parameters. For this SN Ia implementation, the uncertainty on Hubble residuals is computed as:

$${\sigma }_{\mathrm{\Delta }\mu ,i}^2({\mathrm{\Theta }}_{\mathrm{HD}})=\left(1\hspace{0.5em}\hspace{0.5em}\alpha \hspace{0.5em}\hspace{0.5em}-\beta \right)·C_{\mathrm{lcfit},i}·\left({1\hspace{0.5em}\hspace{0.5em}\alpha \hspace{0.5em}\hspace{0.5em}-\beta }^T\right)+{\sigma }_M^2,$$(36)

where the C_lcfit,i is the covariance of the light-curve fit for SN Ia i, and σ_M is the SN Ia intrinsic scatter (assumed to be color-independent). This explicit parameterization allows us to fit simultaneously for SN Ia Hubble diagram's nuisance parameters ${\mathrm{\Theta }}_{\mathrm{HD}}=\left\{M_0,\alpha ,\beta ,{\sigma }_M\right\}$ along with the cosmological parameters of interest. This type of parameterization can be adapted for TF and FP empirical relations, and we plan to develop it in future studies.

The velocity estimator presented here or in Equation (33) can have an average not necessarily null. It can rise from an actual bulk flow in the cosmic web. In this model, the SN Ia intrinsic magnitude dispersion is also added as a free parameter to the covariance matrix at the fitting stage. We added an option to all velocity estimators to fit for a velocity zero point offset.

The galaxy density field, δ_g, is computed from a galaxy catalog using mesh assignment schemes available in pypower³ software. The galaxy density field in a specific cell i is defined by ${\delta }_{\mathrm{g},i}=\left(N_{\mathrm{cell},i}/N_{\exp ,i}\right)-1$ where N_cell,i is the number of galaxies associated with the cell, i, and N_exp,i is the expected number. The latter is computed using a random catalog normalized to the same number of galaxies as the main catalog. The number of assigned galaxies, N_cell,i, can vary depending on the adopted resampling scheme: nearest grid point (NGP), cloud-in-cell (CIC), triangular-shaped cloud (TSC), or piecewise cubic spline (PCS). This resampling method changes which neighboring cells are associated with a galaxy, effectively changing the smoothing of the density field. The uncertainty associated with the density field is estimated as ${\sigma }_{\delta ,\mathrm{g},i}=1/\sqrt{N_{\exp ,i}}$, considering that N_exp,i≫1. As the assignment scheme suppresses power at small scales, the input power spectra are corrected by multiplying them to the following grid window function:

$$\mathrm{\Gamma }(k,L)={\langle \frac{8}{L^3}{\left(\frac{sin\left(k_x\frac{L}{2}\right)}{k_x}\frac{sin\left(k_y\frac{L}{2}\right)}{k_y}\frac{sin\left(k_z\frac{L}{2}\right)}{k_z}\right)}^n\rangle }_{\text{k}},$$(37)

where n depends on the resampling scheme chosen (1 for NGP, 2 for CIC, 3 for TSC, and 4 for PCS).

The likelihood implemented in flip to perform the likelihood-based field-level inference is a multivariate Gaussian function:

$$\mathcal{L}[x(\mathrm{\Theta }),C(\mathrm{\Theta })]=(2\pi {)}^{-N_x/2}|C{|}^{-\frac{1}{2}}exp\left[-\frac{1}{2}{x}^T{C}^{-1}x\right],$$(38)

where the covariance matrix used is C=C_xx+C_obs where C_xx is the theoretical field covariance of the scalar field x, and C_obs is an observational covariance matrix. The latter can be a diagonal matrix containing the error bar associated with the field and scaled by a field-dependent parameter σ_x or the observational covariance matrix provided by a given survey.

Furthermore, the flip package contains an option to interpolate the field covariance matrix when parameters to fit do not factor into a coefficient of the linear decomposition of the model power spectrum in Equation (5). For peculiar velocities, we interpolate the covariance matrix depending on the velocity position FoG parameter σ_u and include that in the fit. For density models that do not include the FoG parameter σ_g in the model power spectrum linear decomposition, such as AB20 or RC25, we also interpolate with respect to this parameter. The interpolation can also be performed for two parameters, for example, when the fit is performed jointly for galaxy density δ_g and velocity v field. For that case, the data vector is given by

$$x=\left(\begin{array}{c}\hfill {\delta }_{\mathrm{g}}\hfill \\ \hfill v\hfill \end{array}\right),$$(39)

and the field covariance is organized as:

$$C_{\mathit{xx}}=\left[\begin{array}{cc}{C}_{\delta \delta }\hfill & C_{\delta v}\hfill \\ C_{v\delta }\hfill & C_{\mathit{vv}}\hfill \end{array}\right].$$(40)

The shape of the likelihood implies that the field x is considered as Gaussian. This approximation holds for noisy fields and large cosmological scales but generally breaks down on small scales. We plan to extend the likelihood of the flip software to account for the observational effects changing the field's Gaussianity, such as selection effects. Finally, the likelihoods are minimized in flip either with the best-fit minimizer iminuit⁴ (James & Roos 1975; Dembinski et al. 2020), or by performing a Markov-chain Monte Carlo (MCMC) sampling with the emcee⁵ (Foreman-Mackey et al. 2013) software. When performing an MCMC sampling, the flip software contains implementations for positive, uniform, and Gaussian priors P(Θ). The posterior distribution P(Θ∣x) which is sampled is then given by the Bayes theorem:

$$P(\mathrm{\Theta }\mid x)=\frac{\mathcal{L}[x(\mathrm{\Theta }),C(\mathrm{\Theta })]P(\mathrm{\Theta })}{P(x)},$$(41)

where the evidence P(x) is not varied and used for normalizing the posterior distribution to unity.

The flip software also contains a module for isotropic power spectrum generation for velocity divergence and matter. Those power spectrum can be generated using the CLASS⁶ (Blas et al. 2011) Boltzmann solver, and velocity dispersion models are generated following the Bel et al. (2019) models.

4. Validation of the flip software

The main objective of this section is to show the additional applications of the flip software while validating the field covariance models we implemented. In particular, we aim at comparing flip to previous codes, comparing the covariance models between each other for a fixed regular grid of coordinates, and validating those models with an estimation of the two-point correlation function on an N-body simulation.

4.1. Comparison with previous code

In previous SN Ia studies, the velocity covariance matrix C_vv was generally calculated using the code pairV (Hui & Greene 2006). This code was created to show that the coherent motion caused by large-scale structure is in fact important for parameter estimation in SN Ia cosmology, and as such calculates the coherent-motion-induced magnitude covariance between arbitrary points in space, as demonstrated by Davis et al. (2011). The pairV code can therefore be used as the standard to compare with flip, numerically and performance-wise. To demonstrate this, we generate a mock set of positions mimicking a SN Ia survey at a range of distances and angular positions and calculate the covariance using each code. The sky positions were chosen to have uniform random angular separations from 0° to 180°, and the redshifts were normally distributed with a mean z = 0.025 and a standard deviation of z = 0.025 (with a minimum distance cut of 10 Mpc) to approximately reproduce a low-z survey.

The ratio of these matrices is shown in Figure 3. We used the same cosmology (that of Carreres et al. 2023) for both and the C23 model (equivalent to RC25 velocity) in flip, and converted them to be in the same units (arbitrarily, flip to magnitude-space). We used the native power spectrum for each code, which we chose to be fully linear for this comparison. For display purposes, the data are ordered by redshift to show the dependence on separation. For ease of comparison between the codes, we created a Python version of pairV called pypairV⁷ in which we have only updated the original code from FORTRAN77 to Python and validated that they give precisely the same results.

Fig. 3. Ratio of Cvv as calculated using the C23 model (equivalent to RC25 velocity) and the pairV for a mock set of 3D positions using the same cosmology. The agreement is on average very close, but there are regions with instabilities.

The agreement between the codes is excellent, but the different choices in each code cause some structured residuals. The two main ones are the power spectrum generation and the covariance integration method. The pairV code opted to approximate the power spectrum with the transfer function developed by Eisenstein & Hu (1998) that very closely reproduces the full Boltzmann code calculation, by for example CLASS, used in flip. Both are shown in Figure 4. The percent-level differences between the power spectra also cause percent-level gradients in the covariance ratio.

Fig. 4. Comparison of the linear power spectra generated using flip and pairV. Whereas flip makes use of the CLASS Boltzmann code, pairV uses a fitting formula from Eisenstein & Hu (1998). The differences at the Pθθ level cause similar differences at the covariance level.

The most striking feature, however, is the numerical instability in the residual plot (Figure 3). This is caused primarily by the method pairV uses to decompose the integral over the power spectrum. Neglecting prefactors and focusing on the power spectrum integration, the velocity covariance can be written as the addition of a parallel component to the line-of-sight, Π(r), and perpendicular, Σ(r), as (Gorski 1988; Gordon et al. 2007; Davis et al. 2011),

$$C_{\mathit{vv}}^{\mathit{ij}}\propto ({\hat{\text{r}}}_1·\hat{\text{r}})({\hat{\text{r}}}_2·\hat{\text{r}})\mathrm{\Pi }(r)+\left({\hat{\text{r}}}_1·{\hat{\text{r}}}_2-({\hat{\text{r}}}_1·\hat{\text{r}})({\hat{\text{r}}}_2·\hat{\text{r}})\right)\mathrm{\Sigma }(r),$$(42)

where the vectors are those depicted in Figure 1. Focusing on the integration part, the parallel and perpendicular components are defined as

$$\mathrm{\Pi }(r)\propto {\displaystyle {\int }_k}\frac{\mathit{dk}}{2{\pi }^2}{P}_{\theta \theta }\left(j_0(\mathit{kr})-\frac{2j_1(\mathit{kr})}{\mathit{kr}}\right),$$(43)

and

$$\mathrm{\Sigma }(r)\propto {\displaystyle {\int }_k}\frac{\mathit{dk}}{2{\pi }^2}{P}_{\theta \theta }{j}_1(\mathit{kr}).$$(44)

The components are pre-computed and then interpolated within the pairV code, which enables the speed-up compared to the direct, observer-centric approach (for the proof of equivalence of both forms, see Davis et al. 2011). The separation-centric form is used due to its computation speed; while pairV also contains the observer-centric form, it becomes unfeasible to use beyond even a couple of hundred SNe Ia because of a large increase in the CPU time requirement. However, the same calculation can be performed by flip on the order of seconds, rather than pairV's many hours.

Figure 5 shows the decomposition components as a function of separation, where it can be seen that Π(r) becomes numerically unstable at large separations and also turns negative at around 150 Mpc for the chosen cosmology. As a direct consequence, this is the same scale above which the numerical instabilities appear in Figure 3; the curves of high instability delineate the regions where separations are always larger than 150 Mpc (upper left and lower right) and always below 150 Mpc (lower left). In addition to this source of instability, there is also the fact that the covariance at large separation is so small that tiny differences cause large ratios. At low separations (lower left region) the cause of instability is not as clear, although it improves notably when using identical power spectra between codes. We conclude that flip calculates an almost identical velocity field covariance matrix as the standard pairV code on the range of scales where the integration is performed similarly, and any differences are ultimately negligible for cosmology.

Fig. 5. Parallel (Π, solid blue line) and perpendicular (Σ, solid orange line) correlation components as a function of physical separation. Since Π(r) turns negative at around 150 Mpc, we plot −Π(r) with a dashed blue line.

4.2. Covariance contraction

All the field covariance models in Table D.1 are integrated consistently in the flip framework, making it possible to perform a model comparison. By expressing the field covariance matrix on a fixed regular coordinate grid, a step that we call contraction, we can express a model giving a theoretical estimator of the two-point correlation function.

We link the coordinates defined in Figure 1 to a basis commonly used in the two-point correlation modeling, i.e., line-of-sight and transverse separation (r_⊥,r_∥) or (r,μ). We define this basis with respect to the vector d so that the $r=\hat{\text{r}}$ and $\mu =\hat{\text{k}}·\hat{\text{d}}$ definition are equivalent to the one we previously defined. The coordinates r_⊥ and r_∥ are defined as transverse and parallel separation such that r_∥=rcos(ϕ) = rμ and r_⊥=rsin(ϕ).

In order to compute the field covariance model, we need to express the (r,ϕ,α) parameters in the (r_⊥,r_∥) or (r,μ) basis. Considering the previous basis definitions, only the α parameter cannot be directly expressed. Consequently, for wide-angle models, which explicitly depend on α, we need to provide additional information to define this angle. It can be done by fixing r₁ to a specific value which becomes a reference for expressing the strength of wide-angle correction. The impact of the r₁ reference position is shown in more detail in Appendix F. Fixing this reference, we can express the α angle depending on the definition of the line-of-sight reference d (see Equation 1), such that:

$${\alpha }_{\mathrm{bisector}}=2arcsin\left(\frac{rsin(\varphi )}{\sqrt{r_{1,\parallel }^2+r_{1,\perp }^2}+\sqrt{{\left(r_{1,\parallel }+r_{\parallel }\right)}^2+{\left(r_{1,\perp }+r_{\perp }\right)}^2}}\right)$$(45)

$$\begin{array}{cc}\hfill {\alpha }_{\mathrm{midpoint}}=& arcsin\left(\frac{rsin(\varphi )}{2\sqrt{r_{1,\parallel }^2+r_{1,\perp }^2}}\right)\hfill \\ \hfill & +arcsin\left(\frac{rsin(\varphi )}{2\sqrt{{\left(r_{1,\parallel }+r_{\parallel }\right)}^2+{\left(r_{1,\perp }+r_{\perp }\right)}^2}}\right)\hfill \end{array}$$(46)

$${\alpha }_{\mathrm{endpoint}}=arctan\left(\frac{r_{\perp }}{\sqrt{r_{1,\parallel }^2+r_{1,\perp }^2}+r_{\perp }}\right).$$(47)

The Appendix F is also showing the impact of the d vector definition on the contraction. For the plane-parallel case, which does not depend on the α angle, there is no need to provide a reference, and the definition of the line-of-sight reference does not change the result. Once the link between the basis is explicitly expressed, we can compute the contraction by computing the field covariance in the same way detailed in Section 3.1, but for a fixed grid of coordinates.

4.3. Comparing covariance models

Figure 6 shows the contraction for the models of Table D.1 on (r_⊥,r_∥) basis, considering as an illustration a reference r_1,⊥ = 0 h⁻¹ Mpc and r_1,∥ = 20 h⁻¹ Mpc to highlight the effect of the wide-angle model. To perform the covariance matrix calculation, the three power spectrum terms P_mm, P_mθ, P_θθ are computed using the modeling from Bel et al. (2019) and the CLASS Boltzmann solver with nonlinear halofit correction. The wavenumbers used are logarithmically spaced in the range k∈[0.0005,1.0] h Mpc⁻¹. We chose a rather large value of maximal wavenumber because this article aims to test our model and the previous ones and to extend their applicability to larger wavenumber ranges. The cosmological parameters used are following AbacusSummit (Maksimova et al. 2021; Garrison et al. 2021) baseline simulation. We also apply a grid window correction to the power spectrum for the simplest grid model (NGP), i.e. multiplying P_mm by Γ(k,L)² from Equation (37) (n = 1), and P_mθ by Γ(k,L). We take a grid size value of L = 10 h⁻¹ Mpc, which is smaller than typical values used but it minimizes the impact of this mesh assignment correction. To compute the contraction, we express the w_ab,n parameters of all models with the example values fσ₈ = 0.4, bσ₈ = 0.8, and β_f=f/b = 0.5. We use an RSD FoG parameter σ_g = 1.0 h⁻¹ Mpc and for the velocity side, we choose a value of σ_u = 15 h⁻¹ Mpc, which is a commonly used value for this nuisance parameter (Koda et al. 2014; Carreres et al. 2023).

Fig. 6. Expression of the field covariance for a uniform coordinate grid, i.e. covariance contraction, for all the models implemented in flip. From top to bottom the AB17, AB20, L22, and RC25 models are represented (see Table D.1 for reference). The C23 model is fully equivalent to vv contractions in the AB17, L22, and RC25 models. From left to right are shown the density-density (gg), density-velocity (gv) and velocity-velocity (vv) covariance contractions as a function of transverse (r⊥) and line-of-sight (r∥) separations. The gg contraction is multiplied by the squared separation to highlight the BAO pattern.

For all gg contractions, we see the effect of baryon acoustic oscillations (BAO) at r∼100 h⁻¹ Mpc contained in the input power spectrum. The AB17 model is expressed in a wide-angle framework but does not contain the RSD modeling in the gg correlations, giving a fully isotropic correlation function. The AB20 model contains RSD modeling without the wide-angle definition; thus, the gv and vv contractions are simpler (it is symmetric in r_∥ but insufficient to model low-redshift field covariance. Finally, for the value of σ_g used in the contraction, the L22 model and our RC25 model are qualitatively and quantitatively equivalent. They accurately describe the RSD term in gg while fully modeling the wide-angle for all the considered contractions. We verified that fixing the r₁ reference to a very large value (e.g., r₁∼10 Gpc) gives equal wide-angle and plane-parallel models with equivalent power spectrum models (e.g., L22 and RC25 compared to AB20). This means that our wide-angle models are correctly expressed within the limit of parallel lines-of-sight.

To estimate the validity range of L22 and RC25 models, we computed their contractions for different values of σ_g and report illustrative cases for one-dimensional C_gg and C_gv contractions in Figure 7. The σ_g parameter is varied by a 1.0 h⁻¹ Mpc step. The σ_g value shown in the figure are the same as in Figure 6 (i.e. 1.0 h⁻¹ Mpc), and the first values for which the two models start to disagree (3.0 h⁻¹ Mpc for gg and 7.0 h⁻¹ Mpc for gv). Our model is able to cover a larger σ_g range when integrating in the wavenumber range k∈[0.0005,1.0] h Mpc⁻¹. For too large values of σ_g, the L22 model shows spurious oscillations that tend to increase in amplitude. The gg contraction is the worst case and gives instabilities for σ_g values as low as 3.0 h⁻¹ Mpc. We note that those instabilities are more important for radial separations r_∥. As in Lai et al. (2022), we verified that reducing the maximal wavenumber to k_max = 0.2 h Mpc⁻¹ improves the stability of L22 model with respect to the σ_g parameter. This behavior from the L22 is expected since it is based on the assumption that kμσ_g is small compared to unity, i.e., reducing k or μ stabilizes the model.

Fig. 7. Covariance contraction with the same parameterization than Fig. 6 while varying the σg parameter for L22 and RC25 models. Only Cgg (top) and Cgv (bottom) are represented, as velocity-velocity correlations do not depend on σg. Those contractions are expressed in one-dimension as a function of r⊥ when r∥ = 0 h−1 Mpc (left) and as a function of r∥ when r⊥ = 0 h−1 Mpc (right).

Figure 8 shows the one-dimensional field covariance contraction for the galaxy-galaxy correlation for different maximal wavenumbers. For low maximal wavenumber k_max = 0.2 h Mpc⁻¹, since the galaxy-galaxy input power spectrum is still high at this typical wavenumber value, the cut-off of wavenumber is creating some unwanted oscillation pattern due to Hankel transform. This pattern is absent for vv and vg for which the input power spectrum is damped at large wavenumber. It indicates that using Hankel transform on an undamped power spectrum cut at a wavenumber that is too low is not advised. For higher maximal wavenumber k_max = 1.0 h Mpc⁻¹, the L22 and RC25 models give the same gg contraction that is stable for both models. However, reaching k_max = 3.0 h Mpc⁻¹, L22 becomes unstable due to the approximation used to derive this model, and these instabilities increase when k_max increases. Our model is stable with respect to maximal wavenumber and gives the same results for all values k_max>1.0 h Mpc⁻¹ and for all field covariance contractions.

Fig. 8. Covariance contraction for RC25 (top) and L22 (bottom) models with the same parameterization than Fig. 7, and σg = 3.0 h−1 Mpc, while varying the maximal wavenumber used in the integration and denoted here as kmax. Those contractions are expressed in one dimension as a function of r⊥ when r∥ = 0 h−1 Mpc. For RC25 the kmax = 1.0 h Mpc−1 and kmax = 50.0 h Mpc−1 lines are superimposed.

We conclude that our model is more stable for a wider range of wavenumber, orientation, and σ_g nuisance parameter. Furthermore, our model does not need to introduce an additional nuisance bias parameter b_add to integrate the small scales as is done in Adams & Blake (2020) and Lai et al. (2022).

4.4. Validation on N-body simulation

We want to show whether our new model can reproduce the clustering of an N-body simulation and at which scale it breaks down. It will notably give insight into the minimal separation that can be used with the current models. This study does not aim to develop a fitter for two-point correlation functions, as the flip formalism is not adapted for this, but rather to validate our field covariance models on a more compact representation.

We use the publicly available AbacusSummit N-body simulation (Maksimova et al. 2021; Garrison et al. 2021), and in particular, the halo catalogs of the 25 baseline ΛCDM cosmo000 (Planck 2018 cosmology) simulations. We assign galaxies to the dark matter halos of the simulations by using the halo occupation distribution (HOD) framework introduced in Zheng et al. (2007). We use a vanilla HOD with five parameters, which assumes that the number of galaxies N in a given dark matter halo follows a probabilistic distribution that only depends on the halo's mass M. Moreover, the occupation is decomposed into the contribution of central and satellite galaxies 〈N(M)〉=〈N_cen(M)〉+〈N_sat(M)〉. The number of central galaxy occupations follows a Bernoulli distribution, while the satellite galaxy occupation follows a Poisson distribution. In particular, we use the AbacusHOD implementation (Yuan et al. 2022), available in the publicly available abacusutils⁸ package. We chose the HOD five parameter ranges following the DESI Bright Galaxy Survey (Prada et al. 2025), a typical sample for this kind of study. To reduce shot noise in our clustering measurements, we draw five independent realizations of galaxies over each halo catalog and using the same HOD model.

This study aims to measure the true clustering of velocity and density in the cosmic web. Therefore, we place ourselves in a flat-sky approximation, do not consider observational effects, and use the true velocities and densities of the galaxy catalog. For the density correlations, because we use a periodic box, we can make use of the simple Peebles estimator (Peebles 1974):

$${\xi }_{\mathit{gg}}(\text{r})=\frac{\mathit{DD}(\text{r})}{\mathit{RR}(\text{r})}-1,$$(48)

where DD(r) = 〈n_g(r₁)n_g(r₂)〉 and RR=〈n_r(r₁)n_r(r₂)〉 the normalized auto pair counts of the data and random catalogs. We consider the full galaxy catalog for computing velocity correlations, which means that we can use the same random catalog. Estimators for velocity autocorrelation and cross-correlation with galaxy density can be derived in a similar way as the Peebles estimator by defining the normalized auto galaxy momentum correlation VV(r) = 〈p_g(r₁)p_g(r₂)〉 and cross momentum-density correlation VD(r) = 〈p_g(r₁)n_g(r₂)〉. In the flat-sky approximation, the Peebles analog for peculiar velocities is given by

$${\xi }_{\mathit{vv}}(\text{r})=\frac{\mathit{VV}(\text{r})}{\mathit{RR}(\text{r})},{\xi }_{\mathit{vg}}(\text{r})=\frac{\mathit{VD}(\text{r})}{\mathit{RR}(\text{r})}·$$(49)

The estimation of two-point correlation functions was performed using the vegacorr⁹ Python package. Additionally, ξ_gg, ξ_gv, and ξ_vv were decomposed into multipoles ξ_ab,ℓ following Equation (A.4). We kept the monopole (ℓ = 0) and quadrupole (ℓ = 2) of the vv and gg terms, and the dipole (ℓ = 1) of the vg term, as they are the only non-null terms. We computed the error bar associated with each multipole as the standard deviation over the 25 realizations. Given the large volume of the simulations used for two-point correlation estimation (8 Gpc³ h⁻³), we scale uncertainties up to match smaller volumes typically measured in peculiar velocity surveys, i.e. 1.12 Gpc³ h⁻³ for a sphere up to z<0.1.

Figure 9 shows the multipoles computed on the AbacusSummit simulation, as well as the best-fit flip contracted model. Those models are generated similarly to Section 4.3 but without the grid window correction since gridding is not used for building the two-point correlation functions. We set k_max = 100 h Mpc⁻¹ to avoid aliasing in the integration without the grid window function. The fit was performed in contracted space directly with iminuit by varying the parameters β_f, fσ₈, bσ₈, σ_u, and σ_g. Since the errors are extremely underestimated at small scales, we used in the fitting procedure a constant uncertainty for each multipole ξ_ab,ℓ equal to the maximal error bar. The χ² function used in the fit is then the sum of χ² for each multipole. Since our model is not able to reproduce the very small scales of the two-point correlation functions, very sensitive to the HOD model considered, we take a minimal fitting value of r_min = 20 h⁻¹ Mpc for the ℓ = 0 and ℓ = 1 multipoles, and r_min = 40 h⁻¹ Mpc for the quadrupoles. To match our measurements, we do not include wide-angle effects in this fit. Consequently, we use our model RC25 expressed at large reference separation, equivalent to the AB20 model.

Fig. 9. Representation of the vv monopole (top left), quadrupole (top right), vg dip (middle left), gg monopole (bottom left) and quadrupole (bottom right) of the two-point correlation functions estimated on the AbacusSummit baseline halo catalog. The orange points are the average of the measurement over 25 simulations, and their associated error bar corresponds the standard deviation over those simulations, rescaled in volume to a sphere of maximal redshift z = 0.1. The plain (RC25) and dashed (AB17) blue lines are the best-fit models determined with flip covariance contraction. The lower panel of each monopole shows the absolute difference between RC25 model and the AbacusSummit multipoles, normalized by the multipole error bars.

For comparison, we also show the best-fit AB17 model, which does not model RSD in the gg and gv correlations. We clearly see that AB17 fails to model the gg quadrupole (by construction) as well as the damping of the BAO peak caused by the FoG RSD effect.

The RC25 model gives a reasonable agreement with the N-body simulation for scales r⪆15 h⁻¹ Mpc for monopoles (ξ_vv,0 and ξ_gg,0), r⪆25 h⁻¹ Mpc for ξ_vg,1, and r⪆45 h⁻¹ Mpc for quadrupoles (ξ_vv,2 and ξ_gg,2). This model is suitable for large linear scales but should be taken with care in the intermediate range scale 15<r<45 h⁻¹ Mpc and ultimately breaks down for the smallest scales. Those discrepancies highlight the potential room for improvement of our model. In particular, we expect nonlinear extensions of the model to improve small scales. On the gg correlation, another potential improvement would be to include HOD modeling. On the vv side, both monopole and quadrupole are controlled only by two parameters (fσ₈ and σ_u), which does not seem to be sufficient for modeling all scales. Our study exposes the need for better nonlinear models of velocity clustering, such as Dam et al. (2021). Finally, we note that this fit is not optimal. We can improve it with a full covariance matrix treatment and better data uncertainties accounting for HOD variability, but this is out of the scope of this study. Furthermore, this analysis represents an extremely accurate noise-free measurement, far beyond current or future observational surveys. The suggested enhancements in the modeling may not be needed when treating current datasets.

5. Survey-dependent Fisher forecast

The framework we detail in Section 2 allows for the fast generation of field covariance matrices following the geometry of a given data sample. The likelihood-based method is time-consuming mainly due to the inversion of the covariance matrix at each likelihood calculation, which can reach many occurrences for minimization or MCMC sampling. Flip implements a faster way to estimate the uncertainties of parameters of interest by computing the Fisher information from the field covariance.

A standard Fisher forecast (Koda et al. 2014; Howlett et al. 2017) is based on the information contained in the peculiar velocity power spectrum. In this method, the Fisher information matrix is computed using

$$\begin{array}{cc}\hfill & F_{\mathit{ij}}=\frac{{\mathrm{\Omega }}_{\mathrm{sky}}}{4{\pi }^2}{\displaystyle {\int }_{r_{min}}^{r_{max}}}{r}^2\mathit{dr}{\displaystyle {\int }_{k_{min}}^{k_{max}}}{k}^2\mathit{dk}{\displaystyle {\int }_0^1}d\mu \mathrm{Tr}\left[{{\text{C}}^{\text{H}}}^{-1}{\text{C}}_{,\text{i}}^{\text{H}}{{\text{C}}^{\text{H}}}^{-1}{\text{C}}_{,\text{j}}^{\text{H}}\right],\hfill \end{array}$$(50)

where Ω_sky is the solid angle from the sky coverage of the survey, r_min and r_max are the lower and upper limits of the considered survey in term of comoving distance, k_min and k_max defines the wavenumber of integration. The covariance matrix C^H is defined in Fourier space for velocity-only Fisher forecast as

$${\text{C}}^{\text{H}}=\left[\begin{array}{cc}\hfill P_{\mathit{vv}}\left(k,\mu \right)+\frac{{\sigma }_{\mathrm{obs}}^2(r)}{{\overline{n}}_v(r)}\hfill & \hfill \hfill \end{array}\right],$$(51)

where a shot-noise term is added with the typical peculiar velocity uncertainty σ_obs(r), and the number density for the considered velocity tracer ${\overline{n}}_v(r)$. The positions of each velocity tracer is not taken into account. Only the redshift distribution of the considered sample and the survey area probed are considered. We created a Python version¹⁰ of the code given in Howlett et al. (2017), and use it to derive the expected uncertainties on fσ₈.

In contrast, our Fisher forecast is computed from the precise position of each data point in the considered field and can account for geometrical effects and velocity/density distribution in a specific survey configuration. Additionally, our forecast can account for various systematic effects thanks to the error vector construction and the potential addition of an observational field covariance matrix. We computed the Fisher information matrix following Tegmark et al. (1997), in the Gaussian case:

$$F_{\mathit{ij}}=\langle {\mathcal{L}}_{,\mathit{ij}}\rangle =\frac{1}{2}\mathrm{Tr}\left[A_i{A}_j+C^{-1}{M}_{\mathit{ij}}\right],$$(52)

with the matrices defined by A_i=C⁻¹C_,i and $M_{\mathit{ij}}={\mu }_{,i}{\mu }_{,j}^t+{\mu }_{,j}{\mu }_{,i}^t$ where the μ is the average of the considered field vector x. This expression can be simplified in our case because we are considering only vectors with null averages (density, velocity, and logarithmic distance), i.e., we have M_ij = 0.

In the formalism of flip, considering the parameters of interest w_ab,n, we can compute the covariance derivatives from Equation (6) such that:

$$C_{\mathrm{ab}}({\text{r}}_1,{\text{r}}_2{)}_{,i}={\displaystyle \sum _n}\frac{\partial {\mathit{w}}_{\mathrm{ab},n}}{\partial {\mathit{w}}_{\mathrm{ab},i}}{C}_{\mathrm{ab},n}({\text{r}}_1,{\text{r}}_2).$$(53)

Using the covariance derivatives, we build the Fisher information matrix given in Eq. (52), and by inverting the Fisher matrix, we obtain the expected uncertainties on the parameters of interest.

We test our Fisher forecast formalism on a ZTF-like simulation of SN Ia peculiar velocities. Specifically, we analyze 27 realizations of a 6-year ZTF survey simulated by Carreres et al. (2023) using the snsim¹¹ package. We use the peculiar velocities extracted from the ZTF simulations to compare the estimated error bar on fσ₈ between flip Fisher, a standard volume Fisher, and a full likelihood-based study.

Figure 10 shows the results for the percentage error on fσ₈ of our flip Fisher forecast, compared to the standard volume Fisher and the likelihood-based method. The likelihood-based results are taken from Figure 13 in Carreres et al. (2023). Our Fisher method achieves errors that are more comparable to the likelihood-based analysis with respect to the standard volume Fisher forecast, especially for larger redshift ranges. We interpret the inaccuracy of the volume Fisher forecast to be caused by the data geometry, and the nuisance parameters not being taken into account. On average, over all the redshift bins, our Fisher method gives uncertainties 30% closer to the likelihood-based estimate in comparison to the standard volume one.

Fig. 10. Percentage error on fσ8 for different redshift bins. The orange points show the results for our flip Fisher forecast, the green points the results for the standard volume Fisher, and the blue points the ones for the likelihood-based method. The results are obtained for different redshift bins with the same minimum redshift, zmin = 0.02, and increasing maximum redshift. The results are plotted at the maximum redshift of each bin and are the average over 27 realizations of a ZTF 6-year survey.

6. Conclusion and future prospects

In this work, we have developed a generalized framework for inferring the growth rate of large-scale structures with the likelihood-based field-level method. The generality of this framework also makes it work with a large variety of peculiar velocity probes (TF, FP, and SNe Ia) and galaxy density estimators. We implemented most of previously used models in the literature and we introduced a new density and velocity covariance model that accounts for wide-angle effects, redshift-space distortions and FoG with better numerical convergence. The covariance generation, the likelihood-based inference, and the wide variety of vector estimators are all implemented consistently in the flip software.

This generalized framework includes several applications that allow us to validate it. First, we compared the field covariance generation to a previous code publicly available. We derived inside our framework the field covariance models derived previously in the literature (Adams & Blake 2017, 2020; Lai et al. 2022; Carreres et al. 2023) and compared them to our new model by expressing the covariance in regular gridding. This contraction of the covariance matrix allows us to test the validity of the models with respect to wavenumber, fitted parameters, and orientation ranges. In particular, we concluded that the model we developed in this work is the most suitable for covering a larger range of wavenumber, orientation, and RSD FoG nuisance parameters than previous models while keeping full wide-angle modeling. By comparing this contraction to two-point correlation functions measured on N-body simulations, we obtain insights into the minimal separation for which our model is still valid.

Secondly, as our framework allows for an efficient computation of field covariance, we developed a Fisher forecast method that accounts for the geometry of the considered survey. We validated this method on an N-body simulation by comparing its performance to a standard volume Fisher forecast software and a full likelihood-based inference. This comparison was performed on a ZTF-like simulation of SNe Ia peculiar velocities. We concluded that since our Fisher forecast method accounts for the entire geometry of the survey, i.e., the position of all the objects used, it gives error bars 30% closer to the full likelihood-based estimate on average in comparison to the standard volume method.

We identify several avenues for extending this work. The most-straightforward improvement is the inclusion of nonlinear models such as the TNS model (Taruya et al. 2010), the regularized PT model at two-loop order model (Taruya et al. 2012), or the effective field theory of large scale structures (EFTofLSS, see e.g., Carrasco et al. 2014; Perko et al. 2016; D’Amico et al. 2021; Ivanov 2023) in the likelihood-based framework. Such models will increase the accuracy of fits on small scales, enabling smaller mesh cell sizes to be used for the density field. The latter can help reducing the uncertainty on fσ₈.

The generalized field covariance framework we developed can also be extended to the covariance of two-point compressed statistics such as a two-point correlation function or power spectrum, allowing one to integrate complex power spectrum models. It can be done by adapting the equations presented in Blake & Turner (2023) to our generalized framework.

The likelihoods and vector estimators developed in our framework can be improved. In particular, we can develop more complex non-Gaussian likelihoods to include observational effects such as selection effects due to magnitude cuts of the considered tracer. Our framework was mainly built to deal with SN Ia peculiar velocities, particularly to simultaneously fit Hubble diagram parameters along with cosmological ones. We plan to adapt this simultaneous fitting to the case of FP and TF peculiar velocity estimators. Those estimators contain normalization parameters characteristic of both empirical relations, and adapting our method will allow us to account for all parameter correlations.

Our survey geometry-dependent Fisher forecast method and the likelihood-based method can be applied to a larger variety of surveys. With the Fisher forecast, we plan to forecast the cosmic variance for the geometry of recent and future surveys, i.e., for which number of data densities the improvement on the fσ₈ error no longer evolves. On the likelihood-based method side, we plan to apply it to the Zwicky Transient Facility (ZTF) (Bellm et al. 2019; Rigault et al. 2025) and the Legacy Survey of Space and Time Rubin-(LSST) (Ivezić et al. 2018) for SNe Ia. In addition, the flip software will allow us to combine SN Ia surveys with galaxy surveys such as the Dark Energy Spectroscopic Instrument (DESI) (DESI Collaboration 2016a, b; Martini et al. 2018) with the peculiar velocity studies (Saulder et al. 2023; Said et al. 2024) or the 4 meter Multi-Object Spectroscopic Telescope (4MOST) (de Jong et al. 2019) with the combination of the Cosmology Redshift Survey (4MOST-CRS) (Richard et al. 2019) with the 4MOST Hemisphere Survey of the Nearby Universe (4MOST-4HS) (Taylor et al. 2023). Finally, the framework we developed will be necessary for consistently combining the fit of all the different samples of peculiar velocities. We plan to extend our framework so that it will be able to treat simultaneously SN Ia peculiar velocities (ZTF, LSST), FP/TF peculiar velocities from galaxy surveys (DESI, 4MOST-4HS), and the galaxy density field (DESI, 4MOST-CRS, 4MOST-4HST).

Acknowledgments

The project leading to this publication has received funding from Excellence Initiative of Aix-Marseille University – A*MIDEX, a French “Investissements d’Avenir” program (AMX-20-CE-02 – DARKUNI). AGK was supported by the U.S. Department of Energy (DOE), Office of Science, Office of High-Energy Physics, under Contract No. DE-AC02-05CH11231.

References

Appendix A: Spherical harmonic properties

A.1. Spherical coordinates

It is useful to integrate over k using the spherical coordinates as

$${\displaystyle {\int }_{\text{k}}}{d}^3\text{k}={\displaystyle {\int }_0^{\infty }}{\displaystyle {\int }_{{\alpha }^{\prime }=0}^{\pi }}{\displaystyle {\int }_{{\varphi }^{\prime }=0}^{2\pi }}{k}^2\mathit{dk}sin{\alpha }^{\prime }d{\alpha }^{\prime }d{\varphi }^{\prime }$$(A.1)

$$={\displaystyle {\int }_k}{k}^2\mathit{dk}{\displaystyle {\int }_{\mathrm{\Omega }}}d\mathrm{\Omega }$$(A.2)

A.2. Legendre polynomial expansion

In the flip formalism we take advantage of the Legendre polynomial expansion. This expansion gives that a function f(x) with x∈[−1,1] can be expanded as

$$f(x)={\displaystyle \sum _{l=0}^{+\infty }}{a}_{\ell }{L}_{\ell }(x),$$(A.3)

with the coefficient a_ℓ such that

$$a_{\ell }=\frac{2\ell +1}{2}{\displaystyle {\int }_{-1}^1}{\mathit{dx}}^{\prime }f(x^{\prime })L_{\ell }(x^{\prime }).$$(A.4)

One useful Legendre expansion is the plane-wave expansion:

$$e^{i\text{k}·\text{r}}={\displaystyle \sum _{l=0}^{+\infty }}(2\ell +1)i^{\ell }{j}_{\ell }(\mathit{kr})L_{\ell }(\hat{\text{k}}·\hat{\text{r}}).$$(A.5)

A.3. Product of Legendre polynomials

A useful theorem to express the product of Legendre polynomials is the spherical harmonic addition theorem, which gives

$$L_{\ell }(\hat{\text{k}}.\hat{\text{r}})=\frac{4\pi }{(2\ell +1)}{\displaystyle \sum _{m=-\ell }^{\ell }}{Y}_{\ell m}(\hat{\text{k}})Y_{\ell m}(\hat{\text{r}}{)}^{*}.$$(A.6)

A.3.1. Product of three Legendre polynomials

From Eq. A.6 we obtains that

$$\begin{array}{cc}\hfill & L_{{\ell }_{}}(\hat{\text{k}}·\hat{\text{r}})L_{{\ell }_1}(\hat{\text{k}}·{\hat{\text{r}}}_1)L_{{\ell }_2}(\hat{\text{k}}·{\hat{\text{r}}}_2)=\frac{(4\pi {)}^3}{(2\ell +1)(2{\ell }_1+1)(2{\ell }_2+1)}\hfill \\ \hfill & \times {\displaystyle \sum _{m=-\ell }^{\ell }}{\displaystyle \sum _{m_1=-{\ell }_1}^{{\ell }_1}}{\displaystyle \sum _{m=-{\ell }_2}^{{\ell }_2}}{Y}_{{\ell }_{}^{}{m}_{}^{}}(\hat{\text{k}})Y_{{\ell }_1^{}{m}_1^{}}(\hat{\text{k}})Y_{{\ell }_2^{}{m}_2^{}}(\hat{\text{k}})\hfill \\ \hfill & \times Y_{{\ell }_{}^{}{m}_{}^{}}^{*}(\hat{\text{r}})Y_{{\ell }_1^{}{m}_1^{}}^{*}({\hat{\text{r}}}_1)Y_{{\ell }_2^{}{m}_2^{}}^{*}({\hat{\text{r}}}_2).\hfill \end{array}$$(A.7)

The Gaunt coefficients are defined such that

$$\begin{array}{cc}\hfill G_{\ell ,{\ell }_1,{\ell }_2}^{m,m_1,m_2}& ={\displaystyle {\int }_{\mathrm{\Omega }}}d\mathrm{\Omega }{Y}_{{\ell }_{}^{}{m}_{}^{}}(\hat{\text{k}})Y_{{\ell }_1^{}{m}_1^{}}(\hat{\text{k}})Y_{{\ell }_2^{}{m}_2^{}}(\hat{\text{k}})\hfill \\ \hfill & =\sqrt{\frac{(2\ell +1)(2{\ell }_1+1)(2{\ell }_2+1)}{4\pi }}\hfill \\ \hfill & \times \left(\begin{array}{ccc}\ell \hfill & {\ell }_1\hfill & {\ell }_2\hfill \\ 0\hfill & 0\hfill & 0\hfill \end{array}\right)\left(\begin{array}{ccc}\ell \hfill & {\ell }_1\hfill & {\ell }_2\hfill \\ m\hfill & m_1\hfill & m_2\hfill \end{array}\right),\hfill \end{array}$$(A.8)