Issue |
A&A
Volume 698, May 2025
|
|
---|---|---|
Article Number | A273 | |
Number of page(s) | 21 | |
Section | Cosmology (including clusters of galaxies) | |
DOI | https://doi.org/10.1051/0004-6361/202554319 | |
Published online | 20 June 2025 |
Generalized framework for likelihood-based field-level inference of growth rate from velocity and density fields
1
Université Clermont-Auvergne, CNRS, LPCA, 63000 Clermont-Ferrand, France
2
Aix Marseille Université, CNRS/IN2P3, CPPM, Marseille, France
3
Department of Physics, Duke University, Durham, NC 27708, USA
4
Korea Astronomy and Space Science Institute, 776 Daedeok-daero, Yuseong-gu, Daejeon 34055, South Korea
5
Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA
⋆ Corresponding author: corentin.ravoux@clermont.in2p3.fr
Received:
28
February
2025
Accepted:
24
April
2025
Measuring the growth rate of large-scale structures (f) as a function of redshift has the potential to break degeneracies between modified gravity and dark energy models, when combined with expansion-rate probes. Direct estimates of peculiar velocities of galaxies have attracted interest as a means of estimating fσ8. In particular, field-level methods can be used to fit the field nuisance parameter along with cosmological parameters simultaneously. This article aims to provide the community with a unified framework for the theoretical modeling of the likelihood-based field-level inference by performing fast field covariance calculations for velocity and density fields. Our purpose is to lay the foundations for a nonlinear extension of the likelihood-based method at the field level. We have developed a generalized framework, implemented in the dedicated software flip to perform a likelihood-based inference of fσ8. We derived a new field covariance model, which includes wide-angle corrections. We also included the models previously described in the literature inside our framework. We compared their performance against ours, and we validated our model by comparing it with the two-point statistics of a recent N-body simulation. The tests we performed have allowed us to validate our software and determine the appropriate wavenumber range to integrate our covariance model and its validity in terms of separation. Our framework allows for a wider wavenumber coverage to be used in our calculations than in previous works, which is particularly interesting for nonlinear model extensions. Finally, our generalized framework allows us to efficiently perform a survey geometry-dependent Fisher forecast of the fσ8 parameter. We show that the Fisher forecast method we developed gives an error bar that is 30% closer to a full likelihood-based estimation than a standard volume Fisher forecast.
Key words: methods: numerical / cosmology: observations / large-scale structure of Universe
© The Authors 2025
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
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σ8 (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 N2, 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 flip1, 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 Cab. 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 r1 and r2, respectively. Their separation vector is r=r1−r2 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 r1 and r2 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 r1, r2, 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:
We computed the covariance of the fields a and b, respectively evaluated at positions r1 and r2. 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:
In addition we define the individual line-of-sight angle cosine where
and
with k and ri are the norms of the vectors k and ri, respectively. We also define the cosine angle
. The power spectrum is defined by
where δD is the Dirac delta function. The corresponding correlation is then defined as the Fourier transform of the power spectrum as:
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
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 Fab,n(k,μ1,μ2) 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 wab,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 wab,n when maximizing the likelihood. Furthermore, integrals of the geometrical part (Fab,n(k,μ1,μ2)) 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
where Cab,n(r1,r2) can be expressed as an integral over the spherical coordinates of the vector k, defined in Eq. (A.1), as
Using the Legendre polynomial expansion (Eq. A.3) on Fab,n(k,μ1,μ2) with respect to μ1 and μ2 we get
where the Lℓ are the Legendre polynomials. For conciseness, we define a wavenumber term that contains the analytical integration of the Fab,n(k,μ1,μ2) term over μ1 and μ2 such that
Combining Eqs. (7), (8) and using the plane-wave expansion (Eq. A.5) to decompose the eik·r term, we obtain
Using Eq. (A.9) to express the three Legendre polynomial product, we obtain
For conciseness, we write the coefficient containing the linear combination of spherical harmonics such that
We can explicitly express the dependency of the spherical harmonic terms inside the terms considering the definition of angles in Figure 1:
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
Finally, the individual covariance terms can be written:
This equation is in an algorithmically optimized form for a fast computation of the field covariance matrix. The integration of the and
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 wab,n can depend on the redshift of the two considered points. Furthermore, the power spectra, , 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 . The covariance of the Eq. (7) is now written as
Using the Legendre polynomial expansion (Eq. A.3), we obtain
Similarly to the wide-angle case, we define the coefficients as
In addition, we use the plane-wave expansion (Eq. A.5) such that the Eq. (17) becomes
Using the result of Eq. (A.12) for the angular integral in Eq. (20), considering as a reminder that , we get
We define the geometrical term , and simplify it using the spherical harmonic addition theorem (Eq. A.6):
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
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 , where the field x, the field covariance Cxx, and the observational covariance Cobs 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 Cab,n(r1,r2), 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:
where the a is the scale factor, H the Hubble parameter, f the growth rate of structure, σ8 the amplitude of the matter perturbations in spheres of 8 h−1 Mpc comoving radius, Pθθ is the velocity dispersion power spectrum normalized by fiducial value, and Du 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:
A simple sympy2 (Meurer et al. 2017) computation of the and
terms (in Equations 10 and 13) states that only following terms are nonzero:
Adding up the different terms to compute the field covariance matrix, we obtain
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:
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−1 Mpc, it means that the model is only valid for low values of k (typically k<0.1 h Mpc−1). It is limiting, especially for an extension to nonlinear models which needs small-scale modes for integration stability, directly at the power spectrum level () or in the full power spectrum model (Pab(k,μ1,μ2)). The authors of the L22 model decided to separate small-scale modeling (0.2<k<1.0 h Mpc−1) 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−1), 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 and
terms (respectively, Nab,ℓ(ϕ) and
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 (Nab,ℓ(ϕ)) are computed for all the field pairs considered. In a second step, the term
is computed numerically by first expressing
(or
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 Cgv and deduce Cvg. Adams & Blake (2017) (Appendix A.3) showed that a minus sign appears between Cvg and Cgv 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 Cvg as the transpose of Cgv.
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 or
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:
where is the logarithmic distance ratio of the galaxy i, zobs,i and zcos,i are, respectively, the observed and cosmological redshifts, and D is the comoving radial distance expressed in h−1 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
where p are the parameters needed to estimate the velocities. The currently implemented estimators are given in Table 1.
Velocity redshift dependence implemented in flip to convert Hubble diagram residuals into peculiar velocities.
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):
where mB,i is the apparent magnitude of SN Ia i in the B-band, M0 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 x1,i and color ci parameters. For this SN Ia implementation, the uncertainty on Hubble residuals is computed as:
where the Clcfit,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 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 pypower3 software. The galaxy density field in a specific cell i is defined by where Ncell,i is the number of galaxies associated with the cell, i, and Nexp,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, Ncell,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
, considering that Nexp,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:
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:
where the covariance matrix used is C=Cxx+Cobs where Cxx is the theoretical field covariance of the scalar field x, and Cobs 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
and the field covariance is organized as:
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 iminuit4 (James & Roos 1975; Dembinski et al. 2020), or by performing a Markov-chain Monte Carlo (MCMC) sampling with the emcee5 (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:
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 CLASS6 (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 Cvv 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 pypairV7 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),
where the vectors are those depicted in Figure 1. Focusing on the integration part, the parallel and perpendicular components are defined as
and
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 and
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 r1 to a specific value which becomes a reference for expressing the strength of wide-angle correction. The impact of the r1 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:
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 r1,⊥ = 0 h−1 Mpc and r1,∥ = 20 h−1 Mpc to highlight the effect of the wide-angle model. To perform the covariance matrix calculation, the three power spectrum terms Pmm, Pmθ, 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−1. 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 Pmm by Γ(k,L)2 from Equation (37) (n = 1), and Pmθ by Γ(k,L). We take a grid size value of L = 10 h−1 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 wab,n parameters of all models with the example values fσ8 = 0.4, bσ8 = 0.8, and βf=f/b = 0.5. We use an RSD FoG parameter σg = 1.0 h−1 Mpc and for the velocity side, we choose a value of σu = 15 h−1 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−1 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 r1 reference to a very large value (e.g., r1∼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 Cgg and Cgv contractions in Figure 7. The σg parameter is varied by a 1.0 h−1 Mpc step. The σg value shown in the figure are the same as in Figure 6 (i.e. 1.0 h−1 Mpc), and the first values for which the two models start to disagree (3.0 h−1 Mpc for gg and 7.0 h−1 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−1. 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−1 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 kmax = 0.2 h Mpc−1 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 kmax = 0.2 h Mpc−1, 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 kmax = 1.0 h Mpc−1, the L22 and RC25 models give the same gg contraction that is stable for both models. However, reaching kmax = 3.0 h Mpc−1, L22 becomes unstable due to the approximation used to derive this model, and these instabilities increase when kmax increases. Our model is stable with respect to maximal wavenumber and gives the same results for all values kmax>1.0 h Mpc−1 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 badd 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)〉=〈Ncen(M)〉+〈Nsat(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 abacusutils8 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):
where DD(r) = 〈ng(r1)ng(r2)〉 and RR=〈nr(r1)nr(r2)〉 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) = 〈pg(r1)pg(r2)〉 and cross momentum-density correlation VD(r) = 〈pg(r1)ng(r2)〉. In the flat-sky approximation, the Peebles analog for peculiar velocities is given by
The estimation of two-point correlation functions was performed using the vegacorr9 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 Gpc3 h−3), we scale uncertainties up to match smaller volumes typically measured in peculiar velocity surveys, i.e. 1.12 Gpc3 h−3 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 kmax = 100 h Mpc−1 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σ8, bσ8, σ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 χ2 function used in the fit is then the sum of χ2 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 rmin = 20 h−1 Mpc for the ℓ = 0 and ℓ = 1 multipoles, and rmin = 40 h−1 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−1 Mpc for monopoles (ξvv,0 and ξgg,0), r⪆25 h−1 Mpc for ξvg,1, and r⪆45 h−1 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−1 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σ8 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
where Ωsky is the solid angle from the sky coverage of the survey, rmin and rmax are the lower and upper limits of the considered survey in term of comoving distance, kmin and kmax defines the wavenumber of integration. The covariance matrix CH is defined in Fourier space for velocity-only Fisher forecast as
where a shot-noise term is added with the typical peculiar velocity uncertainty σobs(r), and the number density for the considered velocity tracer . 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 version10 of the code given in Howlett et al. (2017), and use it to derive the expected uncertainties on fσ8.
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:
with the matrices defined by Ai=C−1C,i and 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 Mij = 0.
In the formalism of flip, considering the parameters of interest wab,n, we can compute the covariance derivatives from Equation (6) such that:
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 snsim11 package. We use the peculiar velocities extracted from the ZTF simulations to compare the estimated error bar on fσ8 between flip Fisher, a standard volume Fisher, and a full likelihood-based study.
Figure 10 shows the results for the percentage error on fσ8 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σ8.
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σ8 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.
flip: Field Level Inference Package https://github.com/corentinravoux/flip
References
- Adams, C., & Blake, C. 2017, MNRAS, 471, 839 [NASA ADS] [CrossRef] [Google Scholar]
- Adams, C., & Blake, C. 2020, MNRAS, 494, 3275 [NASA ADS] [CrossRef] [Google Scholar]
- Bel, J., Pezzotta, A., Carbone, C., Sefusatti, E., & Guzzo, L. 2019, A&A, 622, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2019, PASP, 131, 018002 [Google Scholar]
- Beutler, F., Castorina, E., & Zhang, P. 2019, JCAP, 2019, 040 [CrossRef] [Google Scholar]
- Blake, C., & Turner, R. J. 2023, ArXiv e-prints [arXiv:2308.15735] [Google Scholar]
- Blas, D., Lesgourgues, J., & Tram, T. 2011, JCAP, 2011, 034 [Google Scholar]
- Boruah, S. S., Hudson, M. J., & Lavaux, G. 2020, MNRAS, 498, 2703 [Google Scholar]
- Boruah, S. S., Lavaux, G., & Hudson, M. J. 2021, ArXiv e-prints [arXiv:2111.15535] [Google Scholar]
- Boubel, P., Colless, M., Said, K., & Staveley-Smith, L. 2024, MNRAS, 531, 84 [Google Scholar]
- Carrasco, J. J. M., Foreman, S., Green, D., & Senatore, L. 2014, JCAP, 07, 057 [CrossRef] [Google Scholar]
- Carreres, B., Bautista, J. E., Feinstein, F., et al. 2023, A&A, 674, A197 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Carrick, J., Turnbull, S. J., Lavaux, G., & Hudson, M. J. 2015, MNRAS, 450, 317 [Google Scholar]
- Castorina, E., & White, M. 2020, MNRAS, 499, 893 [Google Scholar]
- Dam, L., Bolejko, K., & Lewis, G. F. 2021, JCAP, 2021, 018 [CrossRef] [Google Scholar]
- D’Amico, G., Senatore, L., & Zhang, P. 2021, JCAP, 01, 006 [Google Scholar]
- Davis, T. M., Hui, L., Frieman, J. A., et al. 2011, ApJ, 741, 67 [NASA ADS] [CrossRef] [Google Scholar]
- de Jong, R. S., Agertz, O., Berbel, A. A., et al. 2019, The Messenger, 175, 3 [NASA ADS] [Google Scholar]
- Dembinski, H., Ongmongkolkul, P., Deil, C., et al. 2020, https://doi.org/10.5281/zenodo.3949207 [Google Scholar]
- DESI Collaboration (Aghamousa, A., et al.) 2016a, ArXiv e-prints [arXiv:1611.00036] [Google Scholar]
- DESI Collaboration (Aghamousa, A., et al.) 2016b, ArXiv e-prints [arXiv:1611.00037] [Google Scholar]
- Djorgovski, S., & Davis, M. 1987, ApJ, 313, 59 [Google Scholar]
- Dupuy, A., Courtois, H. M., & Kubik, B. 2019, MNRAS, 486, 440 [NASA ADS] [CrossRef] [Google Scholar]
- Eisenstein, D. J., & Hu, W. 1998, ApJ, 496, 605 [Google Scholar]
- Ferreira, P. G., Juszkiewicz, R., Feldman, H. A., Davis, M., & Jaffe, A. H. 1999, ApJ, 515, L1 [NASA ADS] [CrossRef] [Google Scholar]
- Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
- Garrison, L. H., Eisenstein, D. J., Ferrer, D., Maksimova, N. A., & Pinto, P. A. 2021, MNRAS, 508, 575 [NASA ADS] [CrossRef] [Google Scholar]
- Gordon, C., Land, K., & Slosar, A. 2007, Phys. Rev. Lett., 99, 081301 [Google Scholar]
- Gorski, K. 1988, ApJ, 332, L7 [Google Scholar]
- Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [NASA ADS] [CrossRef] [Google Scholar]
- Howlett, C. 2019, MNRAS, 487, 5209 [Google Scholar]
- Howlett, C., Robotham, A. S. G., Lagos, C. D. P., & Kim, A. G. 2017, ApJ, 847, 128 [NASA ADS] [CrossRef] [Google Scholar]
- Howlett, C., Said, K., Lucey, J. R., et al. 2022, MNRAS, 515, 953 [NASA ADS] [CrossRef] [Google Scholar]
- Hui, L., & Greene, P. B. 2006, Phys. Rev. D, 73, 123526 [NASA ADS] [CrossRef] [Google Scholar]
- Huterer, D., Shafer, D. L., Scolnic, D., & Schmidt, F. 2017, JCAP, 2017, 015 [CrossRef] [Google Scholar]
- Ivanov, M. M. 2023, Effective Field Theory for Large-Scale Structure (Springer Nature) [Google Scholar]
- Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2018, ArXiv e-prints [arXiv:0805.2366] [Google Scholar]
- James, F., & Roos, M. 1975, Comput. Phys. Commun., 10, 343 [Google Scholar]
- Johnson, A., Blake, C., Koda, J., et al. 2014, MNRAS, 444, 3926 [NASA ADS] [CrossRef] [Google Scholar]
- Karamanis, M., & Beutler, F. 2021, ArXiv e-prints [arXiv:2106.06331] [Google Scholar]
- Kim, A. G., & Linder, E. V. 2020, Phys. Rev. D, 101, 023516 [NASA ADS] [CrossRef] [Google Scholar]
- Koda, J., Blake, C., Davis, T., et al. 2014, MNRAS, 445, 4267 [NASA ADS] [CrossRef] [Google Scholar]
- Kourkchi, E., Tully, R. B., Courtois, H. M., Dupuy, A., & Guinet, D. 2022, MNRAS, 511, 6160 [NASA ADS] [CrossRef] [Google Scholar]
- Lai, Y., Howlett, C., & Davis, T. M. 2022, MNRAS, 518, 1840 [CrossRef] [Google Scholar]
- Maksimova, N. A., Garrison, L. H., Eisenstein, D. J., et al. 2021, MNRAS, 508, 4017 [NASA ADS] [CrossRef] [Google Scholar]
- Martini, P., Bailey, S., Besuner, R. W., et al. 2018, Proc. SPIE, 10702, 107021F [Google Scholar]
- Meurer, A., Smith, C. P., Paprocki, M., et al. 2017, PeerJ Comput. Sci., 3, e103 [CrossRef] [Google Scholar]
- Peebles, P. J. E. 1974, ApJ, 189, L51 [NASA ADS] [CrossRef] [Google Scholar]
- Perko, A., Senatore, L., Jennings, E., & Wechsler, R. H. 2016, ArXiv e-prints [arXiv:1610.09321] [Google Scholar]
- Pfeifer, S., Valade, A., Gottlöber, S., et al. 2023, MNRAS, 523, 5985 [NASA ADS] [CrossRef] [Google Scholar]
- Prada, F., Ereza, J., Smith, A., et al. 2025, A&A, 698, A170 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Qin, F., Parkinson, D., Hong, S. E., & Sabiu, C. G. 2023, JCAP, 06, 062 [Google Scholar]
- Qin, F., Howlett, C., & Parkinson, D. 2025, ApJ, 978, 7 [Google Scholar]
- Richard, J., Kneib, J. -P., Blake, C., et al. 2019, The Messenger, 175, 50 [NASA ADS] [Google Scholar]
- Rigault, M., Smith, M., Goobar, A., et al. 2025, A&A, 694, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Said, K., Colless, M., Magoulas, C., Lucey, J. R., & Hudson, M. J. 2020, MNRAS, 497, 1275 [NASA ADS] [CrossRef] [Google Scholar]
- Said, K., Howlett, C., Davis, T., et al. 2024, ArXiv e-prints [arXiv:2408.13842] [Google Scholar]
- Saulder, C., Howlett, C., Douglass, K. A., et al. 2023, MNRAS, 525, 1106 [Google Scholar]
- Scolnic, D., Perlmutter, S., Aldering, G., et al. 2019, ArXiv e-prints [arXiv:1903.05128] [Google Scholar]
- Springob, C. M., Magoulas, C., Colless, M., et al. 2014, MNRAS, 445, 2677 [Google Scholar]
- Talman, J. D. 1978, J. Comput. Phys., 29, 35 [NASA ADS] [CrossRef] [Google Scholar]
- Taruya, A., Nishimichi, T., & Saito, S. 2010, Phys. Rev. D, 82, 063522 [NASA ADS] [CrossRef] [Google Scholar]
- Taruya, A., Bernardeau, F., Nishimichi, T., & Codis, S. 2012, Phys. Rev. D, 86, 103528 [NASA ADS] [CrossRef] [Google Scholar]
- Taylor, E. N., Cluver, M., Bell, E., et al. 2023, The Messenger, 190, 46 [NASA ADS] [Google Scholar]
- Tegmark, M., Taylor, A., & Heavens, A. 1997, ApJ, 480, 22 [NASA ADS] [CrossRef] [Google Scholar]
- The mpmath development team 2023, mpmath: A Python Library for Arbitrary-precision Floating-point Arithmetic (version 1.3.0), http://mpmath.org [Google Scholar]
- Tripp, R. 1998, A&A, 331, 815 [NASA ADS] [Google Scholar]
- Tully, R. B., & Fisher, J. R. 1977, A&A, 54, 661 [NASA ADS] [Google Scholar]
- Turner, R. J. 2024, ArXiv e-prints [arXiv:2411.19484] [Google Scholar]
- Turner, R. J., Blake, C., & Ruggeri, R. 2022, MNRAS, 518, 2436 [NASA ADS] [CrossRef] [Google Scholar]
- Valade, A., Hoffman, Y., Libeskind, N. I., & Graziani, R. 2022a, MNRAS, 513, 5148 [NASA ADS] [CrossRef] [Google Scholar]
- Valade, A., Libeskind, N. I., Hoffman, Y., & Pfeifer, S. 2022b, MNRAS, 519, 2981 [Google Scholar]
- Yuan, S., Garrison, L. H., Hadzhiyska, B., Bose, S., & Eisenstein, D. J. 2022, MNRAS, 510, 3301 [NASA ADS] [CrossRef] [Google Scholar]
- Zheng, Z., Coil, A. L., & Zehavi, I. 2007, ApJ, 667, 760 [Google Scholar]
Appendix A: Spherical harmonic properties
A.1. Spherical coordinates
It is useful to integrate over k using the spherical coordinates as
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
with the coefficient aℓ such that
One useful Legendre expansion is the plane-wave expansion:
A.3. Product of Legendre polynomials
A useful theorem to express the product of Legendre polynomials is the spherical harmonic addition theorem, which gives
A.3.1. Product of three Legendre polynomials
From Eq. A.6 we obtains that
The Gaunt coefficients are defined such that
where we used the Wigner 3-j symbols.
We then can express the integral of the product of three Legendre polynomials over a solid angle as
A.3.2. Product of two Legendre polynomials
From Eq. A.6 we obtains that
Spherical harmonics are orthonormal, this property can be expressed as
where δK is the Kronecker's delta.
We then obtain
Appendix B: Extension of covariance calculation with redshift dependency
Both wide-angle and plane-parallel power spectrum models can be extended to directly account for the redshift dependency of the parameters to fit. We consider that this redshift evolution can be written at the power spectrum stage such that:
where z1 and z2 are the redshifts of the considered pair. This decomposition is always valid when considering linear models for velocities and densities, both with and without wide-angle.
The covariance of one pair of points can then be written:
and the full covariance can be easily computed using the outer product wa,n and wb,n for all redshifts considered.
The redshift dependencies are optionally implemented in the flip software. This extension is used to account for the redshift dependence of the growth rate in velocity terms, often parameterized with the growth index γ such that
Appendix C: Vector field covariance in the flip framework
Let's consider a vector field x in the basis x,y,z, such as that for the velocity vector field. The vector correlation tensor is defined by
where the indexes i and j vary for all axes x,y,z. By projecting the (r⊥,r∥) basis into x,y,z we can obtain the following expression for the tensor components:
For any directional unit vector and
the correlation of the projected vector field is then defined by:
Similarly to the formalism developed in section 2, we can express the field covariance matrix of projected vectors by performing additional Legendre expansions for the directional vectors. An expression for the generalized field covariance of two vector fields a and b projected over the direction and
is then given by:
where r=r2−r1, ,
,
and
. We make note that this expansion is not yet implemented in the flip framework.
Appendix D: Decomposition of previous models in the flip framework
The detail decomposition of AB17 Adams & Blake (2017), AB20 Adams & Blake (2020), L22 (Lai et al. 2022), and RC25 (This work) models in the flip framework detailed in section 2 is given in table D.1
Decomposition, of all the power spectra models used in Adams & Blake (2017, 2020), Lai et al. (2022), Carreres et al. (2023), and the one developed in this study.
Appendix E: Numerical considerations
For some models in table D.1, we need special care when computing the (or
for plane-parallel) terms that are used in the computation of the field covariance. In particular, when an exponential term is included in the Fab,n part of the power spectrum model, for example for AB20 or RC25,
and
terms correspond, respectively, to a linear combination of very large values when the product kσg reaches low values. For the AB20 model, in Appendix A of Adams & Blake (2020), the Kgg,n,l terms (corresponding to the
product in the flip formalism) have terms that are proportional to (kσg)−9. For the wide-angle model RC25, it can be proportional to (kσg)−10 in the worst case. To illustrate this, we use the expression of
in our model RC25 for the third term (n = 2) of the density-density (gg) correlation, for the l1 = 0 and l2 = 2 case:
In figure E.1, we show the shape of this function for the value σg = 2.0 h−1·Mpc. At low wavenumber, this function is a linear combination of very large floats (up to 1020) that cancel out to give a small value of the order of unity or less. It leads to what is sometimes called catastrophic cancellation, i.e., the small difference will be largely misestimated by numpy12 (Harris et al. 2020), leading to large numerical instabilities. To solve this issue, we use the mpmath
13 (The mpmath development team 2023) Python package, which allows us to make floating-point arithmetic calculations with arbitrary precision. We use this package for all the
, which integrates exponential RSD FoG contaminants terms, as shown in the figure E.1 for the aforementioned term. We check that the mpmath solution is enough to cover all reasonable values of σg (typically higher than 1 h−1·Mpc) without instabilities.
![]() |
Fig. E.1. Wavenumber term |
A second numerical issue encountered by our formalism is the instabilities introduced by the use of Hankel transforms in the numerical integration of the equation 16 (or 23 in the plane-parallel case). When performing a Hankel transform, the wavenumber range should cover all the separations (r=|r|) inside the covariance matrix element. For the case of small separations (k∼1.0 h·Mpc−1), we choose to switch to standard Simpson's numerical integration, which is costly but only for a small number of terms in the full covariance. Note that the linear modeling theories we are using are insufficient for those small separations, and we plan to improve nonlinear corrections in future works.
For the special case when we are using a simulation, if the minimal wavenumber of the power spectrum employed to perform covariance calculation is lower than the maximal size of the box, the Hankel transform accounts for a too large number of modes in the integration. It adds additional terms to the field covariance matrix which can make it non-positive defined. We solved this issue by computing the numerical integration with standard Simpson's numerical integration between the lowest wavenumber of the power spectrum and a wavenumber kmin, corresponding to the lowest mode of the considered simulation. We then subtract this additional integral from the Hankel transform to remove those terms. This regularization should only be used for simulations, not for data containing the Universe's large-scale wavenumber modes.
Appendix F: Wide-angle impact on covariance contraction
This section shows complementary tests conducted on the impact of wide-angle on covariance matrix calculation.
Figure F.1 gives the effect of varying the reference position r1 on the gg covariance contraction computed with the RC25 model. This impact is only shown as a function of perpendicular separation r⊥ because the wide-angle effects in the covariance calculation is not sensitive to parallel separation when the perpendicular separation is null. Both for large r1,⊥ and r1,∥ values, the contraction gives the shape of the plane-parallel limit shown in figure 6. At low values of reference position, the impact of wide-angle is significant and gives a non-symmetric shape with respect to the separation.
![]() |
Fig. F.1. Galaxy-galaxy covariance contraction for the RC25 model with the same parameterization than 6 (and σg = 3.0 h−1·Mpc) while varying the reference position r1 in the parallel r1,∥ (left) and perpendicular r1,⊥ (right) direction. Those contractions are expressed in one-dimension as a function of perpendicular separation r⊥ when r∥ = 0 h−1·Mpc. |
Figure F.2 shows the same covariance contraction but varying the definition of the line-of-sight reference vector d in equation 1. It is computed for r1,⊥ = 10 h−1·Mpc and r1,∥ = 10 h−1·Mpc in order to obtain large wide-angle effects. The endpoint definition is unsuitable for d as it gives non-symmetric contraction as a function of parallel separation r∥ for which the covariance should not be sensitive. Furthermore, it does not exhibit the expected asymmetric behavior over the perpendicular separation axis. Endpoint is the less-used definition for wide-angle correction, as it requires arbitrarily selecting one point as a reference over the other. The midpoint and bisector definition gives a similar contraction profile, but the midpoint definition shows a slightly jagged line pattern near 25 h−1·Mpc separation. This test convinces us to adopt the bisector definition, which is already the most widely used for wide-angle correlation computation.
![]() |
Fig. F.2. Galaxy-galaxy covariance contraction for the RC25 model with the same parameterization than F.1 while varying the reference position the definition for the line-of-sight reference vector d between bisector, midpoint, and endpoint. The reference vector taken is r1,⊥ = 10 h−1·Mpc and r1,∥ = 10 h−1·Mpc. |
All Tables
Velocity redshift dependence implemented in flip to convert Hubble diagram residuals into peculiar velocities.
Decomposition, of all the power spectra models used in Adams & Blake (2017, 2020), Lai et al. (2022), Carreres et al. (2023), and the one developed in this study.
All Figures
![]() |
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 the text |
![]() |
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. |
In the text |
![]() |
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. |
In the text |
![]() |
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. |
In the text |
![]() |
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. |
In the text |
![]() |
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. |
In the text |
![]() |
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). |
In the text |
![]() |
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. |
In the text |
![]() |
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. |
In the text |
![]() |
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. |
In the text |
![]() |
Fig. E.1. Wavenumber term |
In the text |
![]() |
Fig. F.1. Galaxy-galaxy covariance contraction for the RC25 model with the same parameterization than 6 (and σg = 3.0 h−1·Mpc) while varying the reference position r1 in the parallel r1,∥ (left) and perpendicular r1,⊥ (right) direction. Those contractions are expressed in one-dimension as a function of perpendicular separation r⊥ when r∥ = 0 h−1·Mpc. |
In the text |
![]() |
Fig. F.2. Galaxy-galaxy covariance contraction for the RC25 model with the same parameterization than F.1 while varying the reference position the definition for the line-of-sight reference vector d between bisector, midpoint, and endpoint. The reference vector taken is r1,⊥ = 10 h−1·Mpc and r1,∥ = 10 h−1·Mpc. |
In the text |
Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.