© The Authors 2026

1. Introduction

The motion of galaxies inspired cosmologists for at least a century when Hubble noted that the redshifting of absorption lines was proportional to distance, thereby discovering the expansion of the Universe (Hubble 1929). The Hubble law, predicted by Friedman some years prior, is simple: cz = H₀d, where H₀ is (the Hubble) constant, z is the redshift, c is the speed of light, and d is the distance. The redshift Hubble used, however, is now known to be a combination of the velocity due to cosmic expansion and the line-of-sight component of a galaxy’s peculiar velocity, namely the velocity caused by the net gravitational pull exerted by the full mass distribution of the Universe (Kaiser 1987; Dekel et al. 1993). At first order, the equation should read cz = H₀d + v_r, with (v_r) the radial component of the peculiar velocity. Thus, the peculiar radial velocity of a galaxy can be measured based on a joint measurement of redshift and distance,

$$\begin{array}{c}\hfill v_r=cz-H_{0}d.\end{array}$$(1)

Catalogs combining redshifts and distances have demonstrated their relevance for cosmology. They enabled the discovery of the accelerated expansion of the Universe (Riess et al. 1998; Perlmutter et al. 1999), yielded a local value of H₀ (see Di Valentino et al. 2025, for a review on the H₀ tension), permitted us to constrain cosmological parameters such as the growth factor fσ₈ (e.g., Carreres et al. 2023; Rosselli et al. 2025) or the matter density parameter Ω_M (Brout et al. 2022), tested nonstandard models of (evolving) dark energy with cold dark matter (CDM) such as w₀w_aCDM (Brout et al. 2022), explored anisotropies in the expansion of the Universe (Giani et al. 2024; Boubel et al. 2025; Kalbouneh et al. 2025b,a; Stiskalek et al. 2026), or investigated the misalignment between local motion and the Cosmic Microwave Background (CMB) dipole (e.g., Sorrenti et al. 2023). Locally, peculiar velocities have been used to map the large-scale structure (e.g., Courtois et al. 2013; Valade et al. 2024) and the discovery of local dynamical structures such as the Great Attractor (Dressler et al. 1987), the Laniakea basin of attraction (Tully et al. 2014), or the dipole repeller (Hoffman et al. 2017). Last but not least, measurements of peculiar velocities can capture the so-called moments of the velocity field (Regos & Szalay 1989, namely the dipole, quadrupole, etc) which can be easily compared to theoretical predictions, providing a valuable test for cosmological models. This is the path we explore in this work.

While systematic redshift surveys with sufficient precision became available in the mid 1980s (de Lapparent et al. 1986), it remains difficult to measure the physical distances to galaxies. All distance measurements rely on essentially the same principle: the comparison between the observed magnitude or size with an assumed absolute magnitude or size. In the local Universe, where individual stars can be resolved, this is typically more accurate since (e.g.) the magnitude of stars at the tip of the red giant branch (Lee et al. 1993) is a well-established standard candle. Beyond this region, absolute magnitudes need to be derived from supernovae and scaling relations such as the Tully-Fisher (Tully & Fisher 1977) or the fundamental plane relations (Dressler et al. 1987). Supernovae, although observable to great distances, are rare, unfortunately.

For observational reasons, the vast majority of distance measurements are provided in distance moduli μ ∝ log(d). Issues arise when errors on μ (assumed to be normally distributed) are directly placed on the distance, since these then transform into log-normal distance errors, leading to biases (Strauss & Willick 1995). This breathing mode motivated the development of various distance estimators (e.g., Watkins & Feldman 2015; Sorce et al. 2016; Hoffman et al. 2021; Watkins et al. 2023; Sorce et al. 2024), which need to be applied to these data prior to any analysis, such as the moments of the velocity field we explore here. While carefully derived, these methods remain very spontaneous and result in the designing of heterogeneous analysis pipelines.

Methods designed to reconstruct the bulk flow from corrected velocity data (i.e., data from which this log-normal bias was removed) can be separated into three families. First, the maximum likelihood method proposed by Kaiser (1988) and applied to various datasets by, for example, Ma et al. (2011), Ma & Pan (2014), Qin et al. (2018), Qin et al. (2021), Howlett et al. (2022), and Whitford et al. (2023). Second, the minimum-variance method developed by Watkins et al. (2009), later applied by, for example, Feldman et al. (2010), Watkins & Feldman (2015), Watkins et al. (2023), Scrimgeour et al. (2016), Peery et al. (2018), and Whitford et al. (2023). Third, the Wiener filter framework, augmented by constrained realizations to sample uncertainties (Hoffman & Ribak 1992; Zaroubi et al. 1995), which has been applied in successive releases of the Cosmicflows catalogs (CF2, CF3, CF4) (Tully & Courtois 2012; Tully et al. 2016, 2023) and associated analyses (Hoffman et al. 2015, 2021, 2024). The ability to reconstruct the full 3D velocity field of this last approach, assuming a ΛCDM prior and the CMB-derived matter power spectrum, is among its strengths. Since it enforces a theoretical ΛCDM velocity-velocity correlation, this approach cannot truly probe the ΛCDM model, especially in regions beyond a few dozen megaparsec, where the noise dominates the signal and data are scarce.

An alternative approach, known as forward modeling, consists of deriving the full posterior probability of the model parameters of interest (e.g., radial flow, bulk flow, and H₀) given the data using a Bayesian framework. The posterior probability is then explored with a Monte Carlo algorithm. Although computationally expensive, forward modeling surpasses more traditional methods as it can describe the complexity of the distance observations with a fidelity higher than that of conventional distance estimators.

The application of forward modeling to reconstructing the velocity field was first proposed by (Lavaux 2016) and was notably applied by Graziani et al. (2019) to the Cosmicflows-3 catalog. The method has been further improved, as described by Boruah et al. (2021) and Valade et al. (2022), and was applied to the Cosmicflows-4 data notably by Courtois et al. (2023) and Valade et al. (2024). While the forward-modeling approach has shown a greater ability to reconstruct small scales than the traditional Wiener filter approach (Valade et al. 2023), its ability to recover the tension with ΛCDM on the large scales is even more limited than that of the Wiener filter. In addition to assuming an underlying ΛCDM prior, the reconstructions are performed in periodic boxes (e.g., 1 Gpc/h in Valade et al. (2024)), which is known to further dampen the recovered large-scale motion due to the missing large-scale tidal modes. On the one hand, strong priors tend to hide the tension in the data with respect to the model. On the other hand, in the absence of priors, parameters of interest may be less well constrained and the scientific statement is more modest.

A forward-modeling reconstruction of the first spherical modes of the velocity field has been presented by Boubel et al. (2024), who also provided an estimate of the Hubble constant H₀. With a similar model, the anisotropy of the Hubble constant has been discussed by Boubel et al. (2025). The approach developed in these works, however, focused on distances obtained from the Tully-Fisher relation, thus limiting the application of their method to the Cosmicflows-4 Tully-Fisher sample (Kourkchi et al. 2022). Although it is affected by complex magnitude-selection effects, this dataset is relatively isotropic.

We present a forward modeling of the first two spherical modes of the velocity field (radial flow and bulk flow) from peculiar velocity data that makes minimum ΛCDM assumptions. As it relies on the Bayesian framework, the method can directly treat the observational data without relying on an intermediary distance estimator such as Watkins & Feldman (2015). Limiting the reconstruction to the first modes of the velocity enables us to depart from strong priors on the power spectrum of the velocity field, such as in Hoffman et al. (2024), minimizing the effect of the assumed cosmology on the resulting radial and bulk flows. Finally, compared to Boubel et al. (2024), our more generalist, but less realistic, treatment of the data allows us to apply our method to the entire Cosmicflows-4 catalog of peculiar velocities. This leads to results on the velocity flow at greater depth while raising issues of anisotropy coupled to the selection in magnitude (Andersen et al. 2016).

2. Observations

As mentioned in the preceding section, we analyzed the radial peculiar velocity field. For this purpose, we employed the Cosmicflows-4 (CF4) catalog (Tully et al. 2023), which constitutes the most comprehensive collection of galaxy peculiar velocities to date, containing about 56 000 galaxies gathered in about 38 000 groups of galaxies¹ so as to mitigate the effects for nonlinear motion within bound structures. This work is based on the grouped catalog, and all observed redshifts are considered in the CMB frame.

Rather than a single survey, the CF4 catalog is an intercalibrated compilation of several subcatalogs found in the literature: a catalog of catalogs. The fundamental-plane scaling relation (Djorgovski & Davis 1987) provides about three-quarters of the measured distances, with an uncertainty on the distance modulus of σ_μ ∼ 0.4. Many of these data are gathered from the Sloan Digital Sky Survey Peculiar Velocity (Howlett et al. 2022, SDSS-PV) catalog, which extends to z = 0.1 (d ∼ 300 Mpc/h), but is restricted to the relatively narrow footprint; another significant component comes from the 6 degree field galaxy redshift survey (Magoulas et al. 2016, 6dFGRSv), which covers the entire southern celestial hemisphere up to z = 0.055 (d ∼ 160 Mpc/h).

The Tully-Fisher scaling relation (Tully & Fisher 1977) is used to derive the majority of the remaining entries, with an uncertainty similar to that of the fundamental plane. These distances are mostly provided by the CF4-TF catalog (Kourkchi et al. 2022), which is an analysis of SDSS HI fluxes, line widths, and redshifts. CF4-TF covers roughly the entire sky, with a slight preference for the northern celestial hemisphere. The redshift distribution of the Tully-Fisher sample peaks at z = 0.02 (d ∼ 60 Mpc/h) and slowly decreases until z = 0.055.

Cosmicflows-4 contains about 1 000 type 1a supernovae extracted from the Pantheon+ and SH0ES catalogs (Riess et al. 2022; Scolnic et al. 2022, respecitvely), found mostly below z = 0.03 (d ∼ 120 Mpc/h) and reaching z ∼ 0.1. Finally, various higher-precision distance estimators (σ_μ ≤ 0.1) are also present in the CF4 compilation, but are limited to z < 0.02.

The complex footprint of CF4 can be appreciated in Fig. 1, where we show the angular distribution of the catalog, colored by the number of groups in each pixel on the sky in three different radial shells. We considered three different regimes: dense and isotropic below 80 Mpc/h, relatively dense and relatively isotropic between 80 − 160 Mpc/h, and last, strongly anisotropic between 160 − 240 Mpc/h. The last shell was cropped to 240 Mpc/h so as to keep shells of constant thickness. The Milky Way shadows an opening angle of about 10° around the Galactic plane; a region of the sky dubbed the zone of avoidance, where it is relatively impossible to obtain distance modulus measurements due to foreground obscuration.

Fig. 1. Visualization of the CF4 data, split by redshift distance dz. Some of the local clusters within the given range are marked with a black cross with their names. Top: Number of grouped galaxies within the HEALPix pixel. Bottom: Median of radial peculiar velocity for each direction. The pixels containing fewer than 10 data points are shaded in gray.

3. Method

In this section, we detail the model and the method we used to reconstruct the moments of the velocity field from measurements of peculiar velocities. The approach can be summarized as follows. The first two moments of the velocity field (i.e., a radial in or out flow plus a net flow in some direction) were chosen at random. Assuming that the errors in the redshift measurements are negligible, we reconstructed the distance to each galaxy from (a more accurate version of) Eq. (1) (see Eq. (4)). The inferred distances, transformed into distance moduli, were then compared with the Cosmicflows-4 observations. Finally, the moments of the velocity field were updated with a Monte Carlo step. The result is a large collection of probable values of radial and bulk flows compatible with the Cosmicflows-4 data. This approach, dubbed forward modeling, is detailed in the following section.

3.1. Forward modeling

Forward modeling heavily relies on Bayes’s theorem, which encodes the posterior probability of the model parameters π given that the data 𝒟 were obtained as

$$\begin{array}{c}\hfill P(\pi \mid \mathcal{D})\propto P(\mathcal{D}\mid \pi )P(\pi ),\end{array}$$(2)

where the likelihood P(𝒟 ∣ π) and the prior P(π) encode the observational and the underlying physical models, respectively. As discussed in Section 6, we assumed no prior distribution on our parameters, that is, P(π) = 1.

We reconstructed the zeroth and first modes of the velocity field (Regos & Szalay 1989). These are often called the monopole and the dipole, respectively. However, since the term “dipole” has come to be understood in common parlance as the motion of an observer relative to background, we refer with it to the mean motion of galaxies flowing either radially (radial flow) or in some direction (bulk flow).

Throughout this paper, the radial flow is referred to as V_r. Although the radial flow has the dimensions of speed (km/s) and one of two directions (inward or outward), it is not a vector. The bulk flow is a 3D vector in the traditional sense and is referred to here as V = (V_X, V_Y, V_Z). Our parameter space thus consists of four parameters, namely π = (V_r, V_X, V_Y, V_Z). As detailed below, these quantities fit by the forward model represent the velocity field in concentric shells, while the associated quantities in spheres were constructed as detailed in Section 3.3.

We considered 𝒟 = {μ^obs_i}_i < n, where n is the number of observations (namely the number of entries in the CosmicFlows-4 catalog). The redshift does not feature in 𝒟 since the error on it was neglected and therefore did not need to be inferred. Moreover, distance moduli are not per se observables, but were derived from other observables through relations such as the Tully-Fisher relation or the fundamental plane.

3.2. Likelihood

As mentioned above, the goal was to make a suggestion for each galaxy’s (or galaxy group’s) distance modulus, for the purpose of comparison with the measured value. Thus, the distance modulus for each entry was constructed as follows. First, the peculiar velocity of the constraint i was computed as

$$\begin{array}{c}\hfill v_r^i=V_r+\mathit{V}·{\widehat{\mathit{r}}}^i,\end{array}$$(3)

where $\widehat{\mathit{r}}$ is a unit vector pointing toward the constraint. Assuming that the full redshift z is known with negligible uncertainty, we computed the proposed cosmological redshift (namely the redshift due to the Hubble expansion) from the proposed peculiar velocity,

$$\begin{array}{c}\hfill z_{\cos }=\frac{z-v_r/c}{1+v_r/c},\end{array}$$(4)

which is an accurate version of Eq. (1). Then, the proposed physical distance d was obtained by numerically integrating

$$\begin{array}{c}\hfill d=\frac{c}{H_0}{\displaystyle {\int }_0^{z_{\cos }}}\frac{\text{d}z}{\sqrt{{(1+z)}^3{\mathrm{\Omega }}_{\mathrm{M}}+(1-{\mathrm{\Omega }}_{\mathrm{M}})}}.\end{array}$$(5)

Last, the distance modulus was obtained,

$$\begin{array}{c}\hfill \mu =5{log}_{10}[(1+z_{\cos })d]+25.\end{array}$$(6)

The Hubble parameter H₀ and the density parameter Ω_M were not explored by our forward model. Instead, their respective effects on the results are discussed in Section 4 and Appendix B. We note that Eq. (5) makes assumptions on the background cosmology, which are further discussed in Section 6.

After these calculations were carried out, the likelihood function was constructed, namely, the ensemble of proposed values of the distance moduli were compared with the observed values and their errors. The distance modulus errors were assumed to be normally distributed around the mean value and independent of each other, which is the standard approach (e.g., Hoffman et al. 2021; Boruah et al. 2021; Valade et al. 2022; Tully et al. 2023). The likelihood over all the reconstructed distance moduli is expressed as the exponential product of the difference between the observed distance modulus and the proposed one, divided by the error on the distance modulus,

$$\begin{array}{cc}\hfill P(\{{{\mu }^{\mathrm{obs}}}_i\}\mid \pi )& ={\displaystyle \prod _i}P({{\mu }^{\mathrm{obs}}}_i\mid \pi )\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill & \propto {\displaystyle \prod _i}exp(-\frac{{({{\mu }^{\mathrm{obs}}}_i-{\mu }_i(\pi ))}^2}{2{\sigma }_{\mu ,i}^2}).\hfill \end{array}$$(8)

We did not attempt to model the volume or selection effects (i.e., P(d) = constant). This is loosely consistent with the distance distribution of the entire CF4, which is roughly flat, as shown in Fig. 2.

Fig. 2. 3D spatial distribution of the CF4 and mock catalog of each subtype shown as solid and dashed lines, respectively. The shade in a lighter color is the scatter among the mocks.

The posterior distribution was explored by means of a Hamiltonian Monte Carlo (HMC) technique (Neal 2011). An HMC is an efficient Monte Carlo technique that enables the exploration of the full permissible parameter space. Each Monte Carlo chain comprises 1000 states, preceded by a 100-step burn-in phase. Given the simplicity of our problem, any reasonable choice of hyperparameters led to an optimal convergence of the algorithm. The HMC implementation made use of TENSORFLOW-PROBABILITY (Dillon et al. 2017) with a JAX backend (Bradbury et al. 2018).

3.3. Inference of the velocity field in concentric shells

Since distances are measured via their observed distance moduli, the error on the peculiar velocity (σ_v) for each entry of the catalog is proportional to the measured error on the distance modulus (σ_μ) and its (true) distance. For the vast majority of data points, the redshift distance is more precise than the physical distance inferred from the distance ladder,

$$\begin{array}{c}\hfill {\sigma }_v\approx \frac{log(10)}{5}{\sigma }_{\mu }\times cz.\end{array}$$(9)

Owing to distance-dependent errors, most of the constraining power of the data is found near the observer. This implies that running the inference method on the entire catalog at once would lead to a bulk flow estimated at an effective radius

$$\begin{array}{c}\hfill r_{\mathrm{eff}}=\frac{c}{H_0}\frac{{\sum }_i{z}_i/{\sigma }_{v,i}}{{\sum }_{i}1/{\sigma }_{v,i}}\approx 100\mathrm{Mpc}(\approx 75\mathrm{Mpc}/h),\end{array}$$(10)

which is much shallower than that of the data (whose median distance is 155 Mpc/h). This is due to (1) the decrease in the number of constraints with distance and (2) the increase in individual error on the constraints. Additionally, it is to be noted that the CF4 catalog becomes very anisotropic beyond 165 Mpc/h (z = 0.055).

We designed an approach to overcome this limitation. Essentially, the Monte Carlo fits were first run independently in concentric shells, and then, the values in concentric spheres were reconstructed by averaging the contributions of the shells by their respective volume. More precisely, the catalog was split (according to each galaxy redshift distance r_i ≈ cz_i/H₁₀₀ Mpc/h) into concentric shells of thickness Δr = 20 Mpc/h, in each of which the Monte Carlo algorithm was run separately. The value of the sphere of radius (R) was then computed as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill r_j& =j\times \mathrm{\Delta }r,\hfill \\ \hfill {\mathcal{V}}_j& =\frac{4}{3}\pi (r_j^3-r_{j-1}^3),\hfill \\ \hfill Q^{(s)}(R)& =\frac{{\sum }_j^{r_j<R}{Q}_j^{(s)}{\mathcal{V}}_j}{{\sum }_j^{r_j<R}{\mathcal{V}}_j},\hfill \end{array}\end{array}$$(11)

where Q_j^(s) is the sth Monte Carlo step of the quantity Q in the jth shell (i.e., Q_j). Equation (11) was applied independently to the radial flow (V_r) and the components of the bulk flow (V_X, V_Y, V_Z).

Q^(s)(R) only constitute a pseudo-Monte Carlo chain and not an actual radius-dependent reconstruction. Namely, correlations between fitted values at different radii cannot be interpreted.

3.4. Mock of CF4

Before the forward-model method described above was applied to the CF4 data, we determined its precision and accuracy, namely how well it measures the two modes on the velocity field in which we were interested (the radial infall and the bulk flow). It was thus applied to mock observations of the CF4 grouped catalog. In this way, the estimated velocity field can be compared directly to the ground truth. Specifically, due to its complex footprint (anisotropy, sparse sampling, large errors, and luminosity limitation), the construction of reliable mocks of the CF catalog is nontrivial, and several methods have been developed (e.g., Doumler et al. 2013; Sorce et al. 2016; Qin et al. 2021; Valade et al. 2023; Whitford et al. 2023; Boubel et al. 2024). We note that applying the method to mock catalogs drawn from simulations necessitates choosing a cosmology. In this way, the accuracy of the method loses its prior-free nature. However, we recall that the assumption of a cosmology only affects the interpretation of the results: the method itself remains largely cosmology-free.

The BIG MULTIDARKN-body simulation (Riebe et al. 2013), which covers a volume of L_box = (2.5 Gpc/h)³ with N = 3840³ particles, was used for this purpose. The cosmology of the simulation is similar to that of Planck XIII (h = 0.6777, Ω_M = 0.307; Planck Collaboration XIII 2016). Halos with more than 20 particles were considered, namely a lower-mass resolution of 4.7 × 10¹¹ M_⊙/h.

The cosmological box was sliced into 4 × 4 × 4 = 64 cubic subboxes of side length (625 Mpc/h)³, each of which was sufficiently large to enclose the full CF4 catalog. The observer was set at the center of each subbox. We acknowledge that this simplistic approach implies that the observers are statistically found in a cosmological underdense region, which affects results on the radial flow as discussed in Section 4.1. The box coordinates were (arbitrarily) associated with supergalactic coordinates (SGX, SGY, and SGZ).

To grasp the complexity of the CF4 dataset, we separately constructed three of its subcatalogs: SDSS-PV, 6dFGRSv, and “other”. This enabled us to model the anisotropic footprints of the scaling-relation surveys contributing to CF4, since the catalog anisotropy is dominated by the two wide-angle sky surveys (SDSS-PV and 6dFGRSv).

Our method also focused on modeling a (simple) magnitude-limitation². The large number of instruments and distance-derivation methods on which the CF4 catalog is based makes a magnitude analysis very difficult. Not only are CF4 entries extracted from surveys with different magnitude limitations (among other method-dependent selection effects), but their magnitudes are observed in different bands. Furthermore, the simulation from which the halos are taken is dark-matter only, meaning that neither are halos populated by galaxies, nor are their luminosities simulated. We thus opted for an empirical approach detailed in Appendix A. A more accurate procedure is left to further works.

For each subbox, CF4-mock halos were extracted from the simulation according to the following recipe.

1. An angular mask was applied that mimicked the sky coverage of the surveys: the footprint of the SDSS-PV, the southern celestial hemisphere (dec < 0°) for 6dFGRSv, and the full sky for other. We also included an additional cut of 10° around the Galactic plane so as to model the zone of avoidance. Halos were included or excluded based on whether they were inside or outside of the survey footprints and the zone of avoidance.

2. Each halo was given an observed redshift and distance modulus by applying Eqs. (4) to (6) to the radial peculiar velocity and distance of each halo.

3. Galaxies were painted onto halos by computing an absolute magnitude from the halo mass. For each subcatalog, the empirical luminosity-mass relation was fit alongside an upper (apparent) magnitude cut so as to reconstruct the radial distribution of the observations. The exact procedure is detailed in Appendix A.

4. Halos were then randomly selected out of the magnitude-limited sample.

5. An observational scatter depending on the subcatalog was added on the distance moduli following Table 1.

The mock CF4 catalog was constructed by merging the mock SDSS-PV, 6dFGRSz, and the other data.

In Fig. 2 we show the angular and redshift distributions of the mock observations and the CF4 catalog. The spatial distribution of the CF4 catalog is reproduced with great fidelity, demonstrating the ability of our mock procedure to reproduce the complexity of the CF4 compilation well.

While the distributions of SDSS-PV and 6dFGRSv agree well in Fig. 2, some disparities can be seen in the other component. Although the distribution in d_z is well captured, the angular distributions (SGB and SGL) have minor differences. This is due to the crude modeling of the other component, which in reality contains a number of different anisotropic surveys.

4. Results: Mock data

4.1. Radial flow

The reconstruction of the radial flow of the 64 mocks in spherical volumes out to 300 Mpc/h is presented in the top panel of Fig. 3. The recovered radial flow closely follows the ground truth of the simulations, demonstrating that our method reliably reconstructs the radial flow within the expected statistical uncertainty. The error in the method due to the sampling increases with distance, with a remarkable change in behavior after 160 Mpc/h. As expected, this effect stems from the observational limitations modeled in our mock data, namely scarcity of the data, large observational uncertainties, and luminosity limitations.

Fig. 3. Comparison between the average radial peculiar velocity from the simulation and the fit results. In the top panel, the black line represents the median of the true value across 64 mocks, with the shaded gray regions indicating 68%–95%–99.7% percentile intervals, equivalent to 1, 2, and 3σ. The purple error bars represent the median of the fit results, with 68% error bars. The middle and bottom rows show the difference Δ and ratio R between the estimated and simulation values, respectively. For each of these parameters, the median across the 64 mocks is shown by the central light purple line, with error bars indicating the percentile interval covering the central 68% of the distribution.

In the bottom panels of Fig. 3, the under- or overestimate of the radial flow as a function of distance is shown in two ways. We focus on the misestimate in km/s represented by $\mathrm{\Delta }=V_r-V_r^{\mathrm{sim}}$, namely the estimated infall velocity at a given distance minus the corresponding value from the simulation. We also examined the ratio $R=|V_r|/|V_r^{\mathrm{sim}}|$. The method underestimates the slight outflow found below 50 Mpc/h seen in the simulations. The rising of the median ⟨V_r⟩ within this region reflects the observer’s underdense position. The underestimation of this outflow is due to the (poor) simulation resolution (in the exact region where the observations can see the faintest objects) and implies a disparity between the CF4 catalog and our mocks in the innermost region (also seen in Fig. 2). In the absence of constraints, the mean tends to zero while the scatter increases, but due to the absence of a ΛCDM prior, this convergence is expected to be weaker than it is for field-level inference methods such as Hoffman et al. (2024).

On larger scales, beyond ∼200 Mpc/h, the method slightly underestimates the true infall velocity by about 20 km/s, which is small with respect to the expected ΛCDM cosmic variance (more than 100 km/s at 200 Mpc/h). More prominently, observational limitations amplify the signal starting 160 Mpc/h, reaching a factor of two at the 200 Mpc/h. The reconstructed velocity and true velocity begin to diverge at a distance of about 160 Mpc/h. This result is expected, as CF4 becomes strongly anisotropic beyond this distance, as it is then limited to the SDSS-PV sample, which covers only a small fraction of the sky.

4.2. Bulk flow

How well the forward model proposed here works is shown in the bottom two rows of Fig. 4. Out to ∼140 Mpc/h, where the CF4 catalog is relatively isotropic, the method is able to recover each component with a variance comparable to that of ΛCDM. No significant over- or underestimation is found in any of the components. Beyond 160 Mpc/h, where only the SDSS-PV dominates, the method amplifies the underlying signal, leading to errors that exceed cosmic variance and reach a factor of 5 at the edge of the catalog. We note that the R factor increases faster on the Y component than on X and Z, which we ascribe to the geometry of the survey: the SDSS-PV is concentrated around the Y-axis.

Fig. 4. Same as Fig. 3 but for the three Cartesian components of the bulk flow (VX/Y/Z, first three columns) and its magnitude (right-most panel).

The magnitude of the bulk flow (namely the square root of the sum of the squares of each component) is presented in the right-most panel of Fig. 4. We show the direct application of Eq. (11) to the Monte Carlo computed magnitude Q^(s) = |V|^(s) (in purple). The estimated magnitude of the bulk flow is overestimated at all distances: the greater the distance, the greater the overestimation. We distinguish three regimes: (1) up to 60 Mpc/h, the reconstructed amplitude is accurate, (2) between 60 Mpc/h and 160 Mpc/h, the bulk flow decreases more slowly than the ΛCDM predictions; it is characterized by a slow, distance-increasing amplification, and (3) beyond 160 Mpc/h, the amplification becomes stronger and the reconstructed bulk flow begins to increase artificially, which is at odds with the ΛCDM prediction.

The salient points to be gleaned from Fig. 4 are thus that

• The components and the amplitude of the bulk flow are well estimated within ∼60 Mpc/h;

• A modest amplification is found on all components, mostly Y, between 60 Mpc/h and 160 Mpc/h, resulting in a noticeable overestimation of the amplitude of the bulk flow in this distance range;

• A strong amplification of each of the components and thus of the amplitude of the bulk flow beyond 160 Mpc/h appears to be symptomatic of the strong anisotropy of CF4 in this region.

This analysis shows that the CF4 observational limits affect the reconstructed bulk flow and lead to an artificial amplification of the estimated bulk velocity. Thus, a naive interpretation of the reconstructed bulk flow would result in an important but artificial tension with the cosmological model.

4.3. Simulation-based correction

The tests on mock catalogs conducted in the preceding section not only give us confidence in the regime where we expect the method to perform well, but also indicate how we can interpret (and thus correct for) the regions in which the method overestimates the radial or bulk flow. In each shell of each component of the bulk flow and separately for the radial flow, we can use the value of R and Δ to correct the measured value and thus recover the true velocity field moments. This was done by

1. Removing the mean over- or underestimation, namely subtracting the value of Δ, (the Δ correction);

2. Dividing by a factor of the amplification R (R correction).

The correction we applied was based on a set of mocks, and there is an obvious variance across these mocks and their corrections. In other words, in some of our mocks, the method might perform well, while in others, it might perform worse. In general, we took the median correction, but we also included this variance in the error bar of our corrected results.

To account for the variance stemming from the correction, we write

$$\begin{array}{cc}\hfill P_{\mathrm{\Delta }}(Q_{\mathrm{corr}}\mid D)& ={\displaystyle \int P}(Q_{\mathrm{fit}}=Q_{\mathrm{corr}}-\mathrm{\Delta }\mid D)P(\mathrm{\Delta })\mathrm{d}\mathrm{\Delta },\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill P_R(Q_{\mathrm{corr}}\mid D)& ={\displaystyle \int P}(Q_{\mathrm{fit}}=Q_{\mathrm{corr}}\times R\mid D)P(R)\mathrm{d},R\hfill \end{array}$$(13)

where P(Q_fit ∣ D) is the result of the forward model described above, and P(Δ) and P(R) are the probability distributions of the correction factors. In practice, the uncorrected output of the forward model was redressed separately by each of the 64 R or Δ corrections derived from the 64 mock universes, so as to sample the uncertainty in the correction arising from cosmic variance. The final results consist thus of 64 corrected profiles, on which the median and quantiles were computed. We note that this method might overestimate the uncertainty.

As shown in Figs. 5 and 6, the R and Δ corrections behave differently depending on whether the bulk flow or the radial flow is examined: the Δ correction is better suited to radial measurements, while the R correction works ideally for the bulk flow. We highlight that the correction on the bulk flow was performed on the components, not on the amplitude.

Fig. 5. Application of the corrections to the radial flow fit of the mocks. The median of the radial flow in the simulation, expected to be minimal, among the mocks is shown as a black line, with gray shades indicating the 68%–95%–99.7% intervals demonstrating the cosmic variance. Corrections based on differences and ratios are shown as dashed red and blue lines, respectively, with error bars indicating the 68% interval across the 64 corrected radial flows.

Fig. 6. Same as Fig. 5 but for the bulk flow components and magnitude. Note that corrections were applied to the components, not to the magnitude of the vector.

5. Results: Cosmicflows-4 data

5.1. Hubble parameter and radial flow

Catalogs of calibrated distances, such as CF4, are the best way to obtain a local estimate of the Hubble parameter H₀. A simple regression of the cz = H₀d scatter plot might be run, but this is prone to several issues. First, the fact that most of the information is concentrated near the observer (where the observations are dense and uncertainties are small) leads to a very small effective radius of the fitting. Second, this approach overlooks the large-scale radial modes of the velocity. Furthermore, in the presence of an anisotropic sampling of the sky, such as CF4, the bulk flow might pollute the estimate of H₀ and so might higher modes of the velocity field, although observations tend to show that they have significantly less power (Kalbouneh et al. 2025a).

We introduce a more elaborate approach. For several values of H₀, the radial component of the velocity flow was fit to the data. The corresponding profiles are presented in Fig. 7, where the uncorrected (top) and Δ corrected (bottom) mean radial flow is presented. Two regimes appear. Up to 200 Mpc/h, the velocity profiles are straight, with a slope depending on the value of H₀. This is the expected behavior in the presence of an offset in the value of H₀. Beyond 200 Mpc/h, all the profiles plummet. This can be explained by the combination of (1) the limitation in the sky of the coverage in the SDSS-PV volume, and (2) the overwhelming effect of the Sloan Great Wall (Gott et al. 2005; Valade et al. 2024) in this region, whose distance coincides with the elbow in the radial flow (i.e., about 240 Mpc/h). In essence, all entries in the CF4 catalog beyond 200 Mpc/h fall onto this massive object, creating a local positive radial flow between 200 Mpc/h and 240 Mpc/h and a negative radial flow beyond this distance. A correction of the radial flow changes the local inflow, as indicated in Fig. 3, and the method underestimates the mean radial flow below 100 Mpc/h.

Fig. 7. Radial flow fit of the CF4 data. The top panel shows the result directly from the forward modeling as solid lines, and the bottom panel shows the Δ-corrected result as dashed lines. Each line, corresponding to a different value of H0 assumed in the model, is distinguished by color. The error bars are the 68% interval of the Monte Carlo steps.

In the absence of prior information on the radial flow, all the values of H₀ are equally likely. In other words, without knowing on which scale no radial flow is expected, it is impossible to further constrain the expansion rate. This scale is a prediction of any cosmological model (e.g., a consequence of the power spectrum of initial perturbations). In the ΛCDM model, ⟨V_r⟩ is driven to zero with increasing distance according to the Copernican principle. Symptomatic of the so-called cosmic variance, the progressive convergence to homogeneity can be described by a radius-dependent probability law on ⟨V_r⟩. This allowed us to constrain H₀ at each distance probabilistically. Our method essentially reconstructs the probability of H₀ at a certain distance by marginalizing the fitted radial flows over the values predicted by cosmic variance.

The expected ΛCDM ⟨V_r⟩ is a Gaussian distribution, namely P(⟨V_r⟩) = 𝒩(0, σ_⟨Vr⟩). That this is centered on zero reflects the cosmological principle of homogeneity, while σ_⟨Vr⟩ is the scatter in the radial flow due to cosmic variance. σ_⟨Vr⟩ can be analytically derived assuming a power spectrum or empirically sampled from a ΛCDM simulation. In the scope of the scope of the linear theory of structure formation, ⟨V_r⟩ and σ_⟨Vr⟩ both decrease as one over distance³.

The probability of a given H₀ given the measured ⟨V_r⟩ and the ΛCDM prior assumption of the scatter σ_⟨Vr⟩ is

$$\begin{array}{cc}\hfill P(H_0\mid \mathrm{CF}4,{\sigma }_{⟨V_r⟩})& ={\displaystyle \int P}(H_0\mid ⟨V_r⟩,\mathrm{CF}4)P(⟨V_r⟩\mid {\sigma }_{⟨V_r⟩})\mathrm{d}⟨V_r⟩\hfill \end{array}$$(14)

$$\begin{array}{cc}& ={\displaystyle \int P}(⟨V_r⟩\mid H_0,\mathrm{CF}4)\mathcal{N}(0,{\sigma }_{⟨V_r⟩})\mathrm{d}⟨V_r⟩.\hfill \end{array}$$(15)

We considered that P(H₀ ∣ ⟨V_r⟩, CF4) was obtained from the forward model of Section 3.2 and does not depend on the cosmological model⁴. The equality P(⟨V_r⟩ ∣ H₀, CF4) = P(H₀ ∣ ⟨V_r⟩, CF4) is the application of the Bayes theorem in the absence of a prior on either H₀ or ⟨V_r⟩, which is consistent with the absence of a prior in the forward model. The quantity P(⟨V_r⟩ ∣ H₀, CF4) is effectively the posterior of the forward model described above, whose results are shown in Fig. 7 at various distances. The integral Eq. (15) was computed numerically.

Figure 8 shows the resulting distribution of H₀ for the corrected and uncorrected radial flow estimator at several distances. The distributions appear to be close to normal, and the interpretation is therefore straightforward. The corrected and uncorrected values agree except at the smallest distance (60 Mpc/h). As they are computed in concentric spheres, the values of H₀ at different distances are not independent. To summarize our results, we chose 160 Mpc/h as our fiducial distance. This is the largest distance we trusted – it lies at the edge of the 6dF data, outside the volume dominated by the Sloan Great Wall, and appears insensitive to the correction (see Fig. 5). The final result thus reads H₀ = 75.9 ± 1 km/s/Mpc. We highlight that the uncertainty might increase strongly with a proper treatment of further systematics errors in the data and in the approach⁵.

Fig. 8. Probability distribution of the Hubble constant given the CF4 data and ΛCDM radial flow variance at different sphere radii shown in colors, with μ ± 1σ of each curve. The vertical lines emphasize the curve mean. The marginalized likelihood function is proportional to the (un-) corrected radial flow distribution, whose resulting distributions are shown in the (top) bottom panel.

5.2. Bulk flow

The V_Y and V_Z components of the bulk flow are generally well behaved and fully consistent with the ΛCDM predictions, with and without the R correction. This is seen since the magenta and blue curves are within the 3σΛCDM expectation for essentially all distances, with a tension below 1σ before 160 Mpc/h. Beyond 160 Mpc/h, the R-correction brings the bulk flow components in line with ΛCDM expectations and alleviates the moderate tension found in the direct results of the forward model.

This is markedly not the case for V_X, which exits the 3σ range of ΛCDM at a distance of ∼140 Mpc/h, reaching a maximal tension with the cosmological paradigm at more 4σ level at 240 Mpc/h for the corrected and about 8σ at 240 Mpc/h for the uncorrected curves. Although the uncorrected case looks fully inconsistent with the ΛCDM model in the region 140 − 300 Mpc/h, the corrected V_X is marginally consistent (i.e., at the 3σ level) with ΛCDM between 140 − 220 Mpc/h.

We now analyze the amplitude of the bulk flow, which we show in Fig. 9. The direct application of the method and the corrected bulk flow decrease with increasing radius for the innermost ∼60 Mpc/h. This behavior is typical of ΛCDM. Beyond this distance, however, the CF4 data display a rise which, in the case of the corrected flow, peaks at 350 km/s at ∼140 Mpc/h. This value is well within the 3σ scatter of the ΛCDM model. We assume that this behavior is due to the discrepant V_X component. The (uncorrected) bulk flow, on the other hand, exploded and reached a peak of more than 600 km/s at 240 Mpc/h.

Fig. 9. Bulk flow of the CF4 catalog, given H0 = 75.9 km/s/Mpc. The purple solid line shows the direct application of our method with its 1σ uncertainty band, and the dashed blue line shows the R-corrected results presented in Section 4.3. The color shaded region corresponds to the uncertainty stemming from the forward model while the dotted blue lines delimit the corridor of the uncertainty coming from the R correction described in Eq. (13). Other reported values of the bulk flow are denoted by numbers referring to Table 2, with values derived from the Cosmicflows series shown as circles and those from other surveys shown as squares. The radius for some literature values is approximate for illustrative purposes; the exact values are provided in Table 2.

We summarize the results of the bulk flow by noting the following: while a rise in one component (V_X) and in the amplitude of the bulk flow is rare in ΛCDM, the result presented here is consistent with all other published analyses of CF4, and indeed with a few other studies that were based on other methods and data. Figure 9 presents the values of bulk flow found in the literature at various distances, using various methods. Our approach yields results that are consistent with the literature and the ΛCDM model for the regime out to 60 Mpc/h, where the direct reconstruction of the bulk flow appears to be reliable, and beyond 60 Mpc/h, where the simulation-based correction is effective. The rise in the amplitude of the bulk flow is unusual but appears to be real, or is at least a feature of the currently available data. Some methods even measured even higher bulk flows than presented here, specifically when compared to the corrected bulk flow.

6. Discussion

We have explored a forward-modeling approach to jointly reconstruct the first two modes of the velocity field from measurements of peculiar velocities: the radial inflow, and the bulk flow. The forward model was fit independently on the data, which were split into concentric shells before they were reassembled in spherical volumes, each shell being weighed by its volume. This enabled us to reconstruct the radial and bulk flow as a function of the distance. The method was applied to the Cosmicflows-4 (CF4) dataset after being tested on a simulated mock data. From the analysis of the radial flow in a sphere of 160 Mpc/h, augmented by a simple assumption on the shape of the velocity field on large scales, we provided an estimate of the Hubble parameter H₀.

The application of our method to a series of 64 realistic mock observations of the CF4 catalog yielded two results. First, we quantified the selection effects in space and magnitude on the radial and bulk flow. Although no systematic shift appears, the radial flow and the Cartesian components of the bulk flow are amplified with respect to the underlying signal, as previously discussed in Andersen et al. (2016). This results in an overestimated amplitude of the bulk flow, with a qualitative degradation of the reconstructed amplitude beyond 140 Mpc/h. Second, a simple simulation-based correction was extracted from the application of our method to the 64 mock catalogs. This approach enabled a corrected estimation of the radial and bulk flow of the CF4 data at all distances up to the edge of the data.

The uncorrected bulk flow differs significantly from ΛCDM beyond 140 Mpc/h in the −SGX direction and on the magnitude of the bulk flow. The tension is extreme: the bulk flow measured naively rules out ΛCDM. The increase in the bulk flow around 100 − 200 Mpc/h has been consistently described in the literature, with various amplitudes. Beyond this distance, the amplitude predicted by the uncorrected approach stagnates, while the one predicted by the corrected estimator decreases slowly. We note that our uncorrected estimation of the bulk flow of CF4 is consistent with that of Whitford et al. (2023) and Watkins et al. (2023), which might indicate that their methods are also sensitive to spurious amplification issues. However, our studies of the effect of the CF4 survey footprint (more specifically, the angular and radial selection function) have demonstrated that the method overestimates the true velocity field. When this is taken into account, the results are more consistent with ΛCDM. This does not need to be the case: our ΛCDM based correction does not force agreement between the method and the model. It simply indicates that if we lived in a ΛCDM model, the bulk flow estimated by our forward model would overestimate the true one by a quantifiable amount.

The radial flow in spheres is degenerate with the value of H₀ chosen when running the fitting procedure. Instead of exploring the tension with ΛCDM using the H₀ provided by Tully et al. (2023), we chose to estimate H₀ in the light of the CF4 data in a probabilistic way, assigning each value of H₀ a probability based on the likelihood of that corresponding mean inflow in ΛCDM. This analysis resulted in a value of H₀ = 75.9 ± 1 km/s/Mpc, which is at the higher end of the range of local values of H₀ found by other authors. We note that systematic errors on the value should be treated more carefully to obtain a more realistic uncertainty range. We also highlight that this estimate differs from a crude regression of the cz = H₀d on the data, as in Tully et al. (2023). In the presence of an anisotropic dataset, the radial and bulk flow are correlated, and the bulk flow thus biases the estimation of H₀. Our procedure extracts the radial component from the signal, resulting in a value of H₀ that is not polluted by the directional bulk flow.

Throughout this work, we tried to control the effect of the assumed cosmology on the reported results.

• To reconstruct the radial and bulk flows, an isotropic H₀ was assumed, as well as the existence of a velocity field. No prior on the power spectrum of this field was assumed, however, as opposed to previous works (e.g., Courtois et al. 2023; Valade et al. 2024; Hoffman et al. 2024).

• To derive the distance-redshift relation in Eq. (5), a flat ΛCDM Universe was considered. Any cosmology assuming an FLRW metric (e.g., Robertson 1933) yields a similar relation, however.

• Our series of synthetic observations of the CF4 catalog and the resulting correction were based on ΛCDM simulations. The departure of the reconstructed velocity profile with respect to the ground truth is essentially due to the selection function of the data and is thus a priori unrelated to the simulation cosmology. A series of tests that we do not present here demonstrated, however, that the correction scheme depends on cosmological parameters such as fσ₈. We would like to highlight that assuming ΛCDM in the correction model is consistent with our aim to test ΛCDM as a cosmological model. Furthermore, the correction works in favor of ΛCDM and might only dampen the tension in the data.

• In order to give an estimate of H₀, we assumed the cosmological principle, i.e., that the radial velocities of spheres of increasing radii asymptotically converge to zero. Furthermore, a single-parameter scatter on the radial flow of a sphere of 160 Mpc/h was introduced to model the local cosmic variance. The numerical value of this scatter was extracted from our ΛCDM simulations. We highlight, however, that (1) ΛCDM agrees with linear theory on these scales, making our chosen value more general than ΛCDM; (2) that this scatter only affects the uncertainty on H₀ and not its value, and (3) that this result can easily be extended to any other cosmology, given that the cosmic variance on the radial velocity of sphere of 160 Mpc/h can be derived.

• It appears that the bulk flow reconstructed by our algorithm is mildly sensitive to the value of H₀. Considering, for instance, Planck’s H₀ = 67.77 km/s/Mpc value slightly decreases the magnitude of the bulk flow, but simultaneously induces extreme outflows (on the order of thousands of km/s) at the edge of the data. This only marginally mitigates the tension on the amplitude of the bulk flow while creating a major (unphysical) radial outflow of matter.

This method differs from traditional bulk flow estimators, such as those presented in Whitford et al. (2023), from several perspectives. First, the forward-modeling approach enables the writing of complex, nonlinear, but more realistic observational models (e.g., Boubel et al. 2024), thus bypassing external Gaussianization methods (e.g., Watkins & Feldman 2015; Qin et al. 2021; Hoffman et al. 2021). The fitting of the forward model in shells that are combined in the postprocessing into spheres constitutes a new weighting technique. While a simple cut in the data at the high-redshift end causes all the profiles to plateau to their values at small distances because the information is concentrated near the observer, this approach enables a better estimation of the evolution of the fitted values with distance. Compared to field-level inference approaches (e.g., Courtois et al. 2023; Hoffman et al. 2024), our method does not rely on the linear theory of fields, nor does it assume a correlation function of the velocity field, which is known to dampen the tension between the data and the underlying model. Last but not least, our approach is similar to that of Boubel et al. (2024, 2025). However, we treated the entire CF4 catalog with a more general, but less realistic, model than did Boubel et al. (2024), who focused on the Tully-Fisher relation and were thus limited to about 10 000 constraints out of the 38 000 of the entire CF4 grouped catalog. While interpreting the dipolar anisotropic expansion as bulk flow, Kalbouneh et al. (2025a) reported similar bulk flow results using the ungrouped CF4 data, complementing the findings of Stiskalek et al. (2026) using the redshift-limited (z ≲ 0.05) CF4 subsamples, where no local H₀ anisotropy was detected, and the results remain consistent with ΛCDM at these shallower scales.

In future work, we suggest improving the model to include possible systematics in the observational data, notably issues of selection function or intercalibration, which have been proven to lead to significant spurious flows (Whitford et al. 2023, Fig. 9). The correction framework could be improved by considering more elaborated methods. We also look forward to more peculiar velocity data, specifically, to data covering the heavens more spherically.

Acknowledgments

CD appreciates the DPST scholarship. AV thanks A. Whitford for sharing the literature data used in Fig. 9. AV has been supported by the Agence Nationale de la Recherche of the French government through the program ANR-21-CE31-0016-03. NIL acknowledges funding from the European Union Horizon Europe research and innovation program (EXCOSM, grant No.101159513)/The MultiDark Database used in this paper and the web application providing online access to it were constructed as part of the activities of the German Astrophysical Virtual Observatory as a result of a collaboration between the Leibniz-Institute for Astrophysics Potsdam (AIP) and the Spanish MultiDark Consolider Project CSD2009-00064. The Bolshoi and MultiDark simulations were run on the NASA’s Pleiades supercomputer at the NASA Ames Research Center. The MultiDark-Planck (MDPL) and the BigMD simulation suite have been performed in the Supermuc supercomputer at LRZ using time granted by PRACE. Plots are made with Python libraries: Healpy and Matplotlib. HMC module is made available by TensorFlow Probability, coupled with NumPy, Scipy, and Astropy for handy mathematical constants and expressions.

References

Appendix A: Modeling the magnitude-limited mocks

The mocks presented in this work are based on halos issued from a dark matter-only simulation. In order to properly model the magnitude cut, the halos have to be affected a luminosity. This is done by introducing a simplistic mass-luminosity relationship

$$\begin{array}{c}\hfill \frac{L}{L_0}=\alpha {\left(\frac{M}{M_0}\right)}^{\beta }\end{array}$$(A.1)

where L is the constructed luminosity, M is the tabulated mass of the halos (using the 200c convention) and L₀ = 10¹⁰ and M₀ = 10¹² are normalizing constants. Given the halo mass function P(log(M)) predicted by Tinker et al. (2008), the expected number of halos as function of distance in the presence of a magnitude cut at m_cut is given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill P^{\mathrm{mod}}(r)=& r^2\times {\delta }^K(r<r_{\mathrm{cut}})\hfill \\ \hfill & \times {\displaystyle {\int }_{-\infty }^{m_{\mathrm{cut}}}}P(log(M)={\beta }^{-1}(m-5log(r)-25))\mathrm{d}m\hfill \end{array}\end{array}$$(A.2)

where δ^K is the Kronecker delta function and δ^K(r < r_cut) is a crude modeling of the cut in redshift space present in the 6dF and SDSS catalogs. Note that we have set α = 1 as this parameter is perfectly degenerated with the parameter m_cut.

The parameters β and m_cut are fitted to each one of the major subcomponents of the CF4 catalog. The optimized quantity is

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _i}{(P^{\mathrm{obs}}(r_i)-P^{\mathrm{mod}}(r_i))}^2\end{array}$$(A.3)

where r_i = i × Δr with Δr = 5 Mpc/h and P^obs(r) is the histogram of the CF4 redshift-derived distances.

The results of the fitting procedure are presented in table 1 and in fig. 2 (right panel). The mass-distance relationship, displaying the simulated magnitude cut, is presented in fig. A.1 in the appendix. A visual evaluation of the distance histogram shows a good convergence of the approach. The variations in the CF4 observations appear in agreement with the cosmic variance over our 64 mocks, represented by the shaded corridor around the mean histogram. Even though the simplicity of the model prevents a physical interpretation for the fitted values of m_cut and β, it is noted that all the subcatalogs first display a volume-limited regime before being limited in magnitude, at 60 Mpc/h for the other component, 90 Mpc/h for 6dF and 230 Mpc/h for SDSS.

Fig. A.1. An example of the halo mass as a function of redshift distance of the selected halo in the mock, given the criteria in Table 1. The minimum halo mass cut (equivalent to the minimum magnitude cut) is in a black dashed line.

Appendix B: Effect of Ω_M

Equation (5) introduces a dependence of our model on H₀ and Ω_M. As demonstrated in fig. B.1, the effect of Ω_M is of second order compared to that of H₀. A strong variation of Ω_M (±3σ the uncertainty of Brout et al. (2022)) results in a variation of less than 100 km/s in the radial flow at 240 Mpc/h, comparable to a variation of H₀ of about ≈0.4km/s/Mpc. In future works, this should be added in the systematic errors budget on H₀.

Fig. B.1. Radial flow fit of CF4 data with different matter density parameter, given the same H0.

All Tables

All Figures

	Fig. 1. Visualization of the CF4 data, split by redshift distance dz. Some of the local clusters within the given range are marked with a black cross with their names. Top: Number of grouped galaxies within the HEALPix pixel. Bottom: Median of radial peculiar velocity for each direction. The pixels containing fewer than 10 data points are shaded in gray.
In the text

	Fig. 2. 3D spatial distribution of the CF4 and mock catalog of each subtype shown as solid and dashed lines, respectively. The shade in a lighter color is the scatter among the mocks.
In the text

Fig. 3. Comparison between the average radial peculiar velocity from the simulation and the fit results. In the top panel, the black line represents the median of the true value across 64 mocks, with the shaded gray regions indicating 68%–95%–99.7% percentile intervals, equivalent to 1, 2, and 3σ. The purple error bars represent the median of the fit results, with 68% error bars. The middle and bottom rows show the difference Δ and ratio R between the estimated and simulation values, respectively. For each of these parameters, the median across the 64 mocks is shown by the central light purple line, with error bars indicating the percentile interval covering the central 68% of the distribution.

In the text

	Fig. 4. Same as Fig. 3 but for the three Cartesian components of the bulk flow (VX/Y/Z, first three columns) and its magnitude (right-most panel).
In the text

Fig. 5. Application of the corrections to the radial flow fit of the mocks. The median of the radial flow in the simulation, expected to be minimal, among the mocks is shown as a black line, with gray shades indicating the 68%–95%–99.7% intervals demonstrating the cosmic variance. Corrections based on differences and ratios are shown as dashed red and blue lines, respectively, with error bars indicating the 68% interval across the 64 corrected radial flows.

In the text

	Fig. 6. Same as Fig. 5 but for the bulk flow components and magnitude. Note that corrections were applied to the components, not to the magnitude of the vector.
In the text

	Fig. 7. Radial flow fit of the CF4 data. The top panel shows the result directly from the forward modeling as solid lines, and the bottom panel shows the Δ-corrected result as dashed lines. Each line, corresponding to a different value of H0 assumed in the model, is distinguished by color. The error bars are the 68% interval of the Monte Carlo steps.
In the text

	Fig. 8. Probability distribution of the Hubble constant given the CF4 data and ΛCDM radial flow variance at different sphere radii shown in colors, with μ ± 1σ of each curve. The vertical lines emphasize the curve mean. The marginalized likelihood function is proportional to the (un-) corrected radial flow distribution, whose resulting distributions are shown in the (top) bottom panel.
In the text

Fig. 9. Bulk flow of the CF4 catalog, given H0 = 75.9 km/s/Mpc. The purple solid line shows the direct application of our method with its 1σ uncertainty band, and the dashed blue line shows the R-corrected results presented in Section 4.3. The color shaded region corresponds to the uncertainty stemming from the forward model while the dotted blue lines delimit the corridor of the uncertainty coming from the R correction described in Eq. (13). Other reported values of the bulk flow are denoted by numbers referring to Table 2, with values derived from the Cosmicflows series shown as circles and those from other surveys shown as squares. The radius for some literature values is approximate for illustrative purposes; the exact values are provided in Table 2.

In the text

	Fig. A.1. An example of the halo mass as a function of redshift distance of the selected halo in the mock, given the criteria in Table 1. The minimum halo mass cut (equivalent to the minimum magnitude cut) is in a black dashed line.
In the text

	Fig. B.1. Radial flow fit of CF4 data with different matter density parameter, given the same H0.
In the text