# The peculiar velocity field up to $z \sim 0.05$ by forward-modelling Cosmicflows-3 data

R. Graziani, H.M. Courtois, G Lavaux, Y. Hoffman, R.B. Tully, Y. Copin, D.

Pomarède

► To cite this version:

R. Graziani, H.M. Courtois, G Lavaux, Y. Hoffman, R.B. Tully, et al.. The peculiar velocity field up to $z \sim 0.05$ by forward-modelling Cosmicflows-3 data. Monthly Notices of the Royal Astronomical Society, 2019, 488 (4), pp.5438-5451. 10.1093/mnras/stz078. hal-01990691

HAL Id: hal-01990691

https://hal.science/hal-01990691

Submitted on 31 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The peculiar velocity field up to $z \sim 0.05$ by forward-modelling Cosmicflows-3 data

R. Graziani<sup>1,2*</sup> H. M. Courtois,<sup>1*</sup> G. Lavaux<sup>1,3</sup> Y. Hoffman,<sup>4</sup> R. B. Tully,<sup>5</sup> Y. Copin<sup>1</sup> and D. Pomarède<sup>6</sup>

Accepted 2019 January 7. Received 2019 January 7; in original form 2018 September 6

## ABSTRACT

A hierarchical Bayesian model is applied to the Cosmicflows-3 catalogue of galaxy distances in order to derive the peculiar velocity field and distribution of matter within $z \sim 0.054$. The model assumes the ΛCDM model within the linear regime and includes the fit of the galaxy distances together with the underlying density field. By forward modelling the data, the method is able to mitigate biases inherent to peculiar velocity analyses, such as the Homogeneous Malmquist bias or the lognormal distribution of peculiar velocities. The statistical uncertainty on the recovered velocity field is about $150 \text{ km s}^{-1}$ depending on the location, and we study systematics coming from the selection function and calibration of distance indicators. The resulting velocity field and related density fields recover the cosmography of the Local Universe which is presented in an unprecedented volume of our Universe 10 times larger than previously reached. This methodology opens the doors to reconstruction of initial conditions for larger and more accurate constrained cosmological simulations. This work is also preparatory to larger peculiar velocity data sets coming from Wallaby, TAIPAN, or LSST.

**Key words:** methods: data analysis – galaxies: distances and redshifts – dark matter – large-scale structure of Universe – cosmology: observations.

## 1 INTRODUCTION

Peculiar motions of galaxies are due to the gravitational interaction with the underlying density field of matter. Thus peculiar velocities of galaxies are a powerful and unbiased tool to study the dynamics and structure of the Local Universe. Velocities have been used as probes of cosmological parameters (e.g. Zaroubi et al. 1997; Nusser & Davis 2011; Feix, Branchini & Nusser 2017; Howlett et al. 2017; Nusser 2017; Wang et al. 2018), for cosmography studies (Dekel et al. 1999; Tully et al. 2014; Courtois et al. 2017; Hoffman et al. 2017), and to set initial conditions for constrained simulations (Gottloeber, Hoffman & Yepes 2010; Sorce et al. 2016). In the past, the Local Universe peculiar velocity field has been reconstructed from the expected response to the observed redshift distribution of galaxies taken as tracers of the mass distribution (Hudson et al.

2004; Nusser & Davis 2011; Davis & Scrimgeour 2014; Hong et al. 2014; Scrimgeour et al. 2016).

An alternate approach has been followed by the Cosmicflows program. Peculiar velocities are inferred from departures of measured distances from the expectations of uniform cosmic expansion. Cosmicflows catalogues (Tully et al. 2008, 2013) have been analysed through the Wiener Filter/Constrained Realizations methodology (WF/CR; Zaroubi, Hoffman & Dekel 1999; Courtois et al. 2012). The assumption is made with the WF/CR methodology that the measured velocities have a Gaussian noise and that their two-point correlations are given by the Lambda cold dark matter (ΛCDM) model.

It is critical that steps be taken to mitigate the Malmquist bias that arises from errors in distance. Objects are preferentially misplaced from regions of higher sampling density to lower sampling density (Strauss & Willick 1995). As an example, it can be anticipated that galaxies assigned the greatest distances probably have large positive errors. These galaxies will be attributed with large negative peculiar velocities. The main objective of this paper is to explore

<footnote n="1"><sup>1</sup>Université de Lyon, F-69622, Lyon, France; Université de Lyon 1, Villeurbanne; CNRS/IN2P3, Institut de Physique Nucleaire de Lyon, France</footnote>

<footnote n="2"><sup>2</sup>Université Clermont Auvergne, CNRS/IN2P3, Laboratoire de Physique de Clermont, France</footnote>

<footnote n="3"><sup>3</sup>Institut d'Astrophysique de Paris, 75014, Paris, Sorbonne Université/CNRS, France</footnote>

<footnote n="4"><sup>4</sup>Racah Institute of Physics, Hebrew University, Jerusalem 91904, Israel</footnote>

<footnote n="5"><sup>5</sup>Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA</footnote>

<footnote n="6"><sup>6</sup>Institut de Recherche sur les Lois Fondamentales de l'Univers, CEA Université Paris-Saclay, F-91191, Gif-sur-Yvette, France</footnote>

<footnote>* E-mail romain.graziani@clermont.in2p3.fr (RG); h.courtois@ipnl.in2p3.fr (HMC)</footnote>
a more rigorous solution than the previous approach of Hoffman, Courtois & Tully (2015) in order to overcome this issue.

As a framework, Lavaux (2016) has developed a fully Bayesian algorithm that incorporates the constrained realizations technique within a statistical model accounting for the uncertainty on the location of tracers. We use here a similar method to reconstruct the 3D linear velocity field from Cosmicflows-3 (CF3) data up to redshift $z \sim 0.054$.

The paper is organized as follows: the CF3 data is briefly described in Section 2 then the method is detailed in Section 3. The results include an analysis on the reconstruction of the linear velocity field in Section 4 and an overlook on the resulting cosmography is given in Section 5. There are comments on outstanding issues in Section 6.

## 2 DATA

### 2.1 Compilation of distance moduli

The reconstruction is based on the CF3 catalogue (Tully, Courtois & Sorce 2016).<fnref n="1" /> CF3 provides a compilation of almost 18 000 galaxy distances, with the high redshift ones computed for the most part from three methodologies: the luminosity–linewidth relation of spiral galaxies, TF (Tully & Fisher 1977), the Fundamental Plane of early-type galaxies (Djorgovski & Davis 1987; Dressler et al. 1987), and Type Ia supernovae (Phillips 1993). The absolute scale of the global distance ladder of these methodologies is given by galaxies that overlap with Cepheid variables (Leavitt & Pickering 1912) or tip of the red giant branch (Da Costa & Armandroff 1990). The Surface Brightness Fluctuation (Tonry & Schneider 1988) method helps providing a bridge between the near and far field. The CF3 compilation is heterogeneous, unlike concurrent single methodology samples (Springob et al. 2007, 2014; Hong et al. 2014). Indeed, when available, CF3 incorporates the major literature contributions. For inclusion in CF3, a source must usefully complement other components while overlapping sufficiently to assure consistency of scale. Each linkage has associated uncertainties. However, what is lost in the ambiguities of linkages is surely more than compensated by the dynamic range of the CF3 catalogue. Nearby, coverage is dense and distances are accurate at the level of 5 per cent. Farther away, Fundamental Plane contributions emphasize coverage of major clusters. TF samples preferentially provide distances to galaxies in the field. SNIa hosts are scattered serendipitously. The methodologies converge in groups where there can be multiple contributions. A group, or an individual galaxy where there is a convergence, has a unique distance. Averaged over all such cases, methodologies should agree.

### 2.2 Groups

Our goal is to derive the linear velocity field from the peculiar velocities of galaxies. However, galaxies in groups or clusters are affected by non-linear motions which are not modelled within our ΛCDM linear framework. Our solution is to average information over the small scale of groups. Tully (2015) provides a catalogue of groups built from the 2MASS redshift survey complete to $K_s = 11.75$ (Huchra et al. 2012). Candidate galaxies are either directly linked to these groups as members of the 2MASS sample or indirectly linked by close spatial and velocity association. A group

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Histogram of the observed redshifts for the five subsamples described in Section 2.2.</figcaption></figure>

is assigned a velocity that averages over all known members and a distance that is the weighted average over those constituents with the necessary measurements. Uncertainties with *N* measures are reduced roughly as $\sqrt{N}$ (depending on the details of the contributing methodologies), so groups are particularly high value entries in the CF3 catalogue.

### 2.3 Selection function

Our Bayesian methodology (see Section 3) needs priors on the statistical distribution of distances. Because it is a composite catalogue, CF3 is inhomogeneous both in distance and angular coverage. Consequently it does not admit a unique and simple selection function. Still, it is possible to identify five main subsamples based on the original observational surveys:

(i) 6dFGSv data provides the most well-defined subsample: it has a high degree of completeness up to a sharp cut-off at $z = 0.054$. The subsample contains 5777 galaxies or groups with a median redshift of $z = 0.039$.

(ii) Another reasonably well-defined subsample is based on the TF method with near-infrared photometry from Spitzer Space Telescope. This constituent gives particular emphasis to coverage at low Galactic latitudes. The 1546 galaxies or groups included in this subsample have a median redshift of $z = 0.009$.

(iii) TF data is particularly deep within the region covered by Arecibo Telescope which observes galaxies of declination $\delta \in [0, 38]^{\circ}$. The subsample contains 1628 galaxies or groups with a median redshift of $z = 0.023$.

(iv) Other TF data than in the Arecibo declination range. The subsample contains 1569 galaxies or groups with a median redshift of $z = 0.019$.

(v) Other data coming from heterogeneous methodologies. Cepheids, tip of the red giant branch and surface brightness fluctuation contributions are local and of high accuracy. Groups which have more than ten members with known redshifts are included in this subsample. Group distances are mostly averaged over Fundamental Plane and TF measures with occasional SNIa contributions. Isolated SNIa lie over a wide range of redshifts up to $z = 0.1$. This last subsample does not have a simple selection function. This subsample contains 963 galaxies or groups with a median redshift of $z = 0.015$.

The total number of tracers in the catalogue is *N* = 11 483. Fig. 1 shows the redshift histograms of all these subsamples. We see that

<footnote n="1">¹Available as a table at http://edd.ifa.hawaii.edu</footnote>
6dF data plays a major role in CF3 and has a singular behaviour with a sharp cut-off at redshift $z = 0.054$. Except for the last subsample of heterogeneous inputs, the other subsamples’ redshift distributions behave as expected: a growth at small distance due to the volume effect and then a decrease due to selection effects. We will show in Section 3.4.4 how we model these distributions.

# 3 METHODOLOGY

This section presents the methodology developed to reconstruct the linear velocity field underlying CF3 data. Our framework is assuming a ΛCDM model and a Gaussian error model on the observations, the distance moduli {$\mu_i$}, and redshifts {$z_i$}. As it will be described in the following sections, the parameter space $\Theta$ includes the most probable luminosity distance of galaxies {$d_i$}, an effective Hubble constant $h_{\text{eff}}$, a non-linear velocity dispersion $\sigma_{\text{NL}}$, and the linear velocity field itself $\mathbf{v}$:

$$ \Theta = \{d_i, h_{\text{eff}}, \sigma_{\text{NL}}, \mathbf{v}\}. \qquad (1) $$

The general idea is to derive from basic assumptions the posterior probability of these parameters given the data $P(\Theta|\{\mu_i, z_i\})$ and then sample from it. The sampling will result in the posterior distribution of the linear velocity field $\mathbf{v}$ from which we will extract the mean and standard deviation. The first subsection reminds the reader about the ΛCDM model of peculiar velocity statistics while the next subsections will detail how the posterior probability is constructed and the procedure for sampling from this complex distribution.

## 3.1 Definitions and notations

The methodology relies on the assumption of a flat ΛCDM model parametrized by the set of parameters ($\Omega_m, H_0$). The ΛCDM model assumes homogeneity and isotropy so that the cosmological redshift $\bar{z}$ of a galaxy is related to its luminosity distance $d$ through the Hubble law:

$$ d(\bar{z}) = c \frac{1 + \bar{z}}{H_0} \int_0^{\bar{z}} \frac{dz}{\sqrt{\Omega_{\Lambda} + \Omega_m (1 + z)^3}}, \qquad (2) $$

where $\Omega_{\Lambda} = 1 - \Omega_m$. Equation (2) can be numerically inverted, giving $\bar{z}(d)$. In the following we make the dependence implicit and note $\bar{z}$ instead of $\bar{z}(d)$. The analysis assumes the linear regime where the overdensity field $\delta(\mathbf{r})$ is small $|\delta(\mathbf{r})| \ll 1$. In this case the linear theory of perturbations predicts that the linear overdensity field is Gaussian and is described by its power spectrum $P(k)$:

$$ \langle \delta(\mathbf{k}) \delta^*(\mathbf{k}') \rangle = (2\pi)^3 \delta_D(\mathbf{k} - \mathbf{k}') P(k), \qquad (3) $$

where $\delta(\mathbf{k})$ denotes the Fourier transform of $\delta(\mathbf{r})$ and $\delta_D$ is the Dirac delta distribution. Gravitational dynamics in an expending universe dictates that the rotational component of the velocity decays early after the onset of the instability and to linear order the velocity $\mathbf{v}$ and density fields are related by:

$$ \nabla \cdot \mathbf{v} = -H_0 f \delta(\mathbf{r}), \qquad (4) $$

where $f$ is the growth rate of structure and depends on the assumed cosmological parameters. Equation (4) is linear as we can see by transposing it into the Fourier space:

$$ i\mathbf{k} \cdot \mathbf{v}(\mathbf{k}) = -H_0 f \delta(\mathbf{k}). \qquad (5) $$

Consequently, the statistic of $\mathbf{v}$ is also Gaussian. The statistic of $\mathbf{v}$ is described in the Appendix A.

Because of the peculiar velocity field $\mathbf{v}$, the Hubble flow described in equation (2) is distorted, and the observed redshift $z$ of a galaxy is a composition of the cosmological redshift and the Doppler effect coming from the radial peculiar velocity:

$$ 1 + z = (1 + \bar{z}(d)) \left( 1 + \frac{\mathbf{v}(\mathbf{r}) \cdot \hat{\mathbf{r}}}{c} \right), \qquad (6) $$

where $\hat{\mathbf{r}} \triangleq \mathbf{r}/r$. The velocity field $\mathbf{v}$ can be divided into a linear and a non-linear part $\mathbf{v}_{\text{tot}} = \mathbf{v} + \mathbf{v}_{\text{NL}}$ and the non-linear part is approximated here by an isotropic Gaussian probability distribution function:

$$ P(\mathbf{v}_{\text{NL}}|\sigma_{\text{NL}}) \sim \mathcal{N}(\mathbf{v}_{\text{NL}}; \mathbf{0}, \sigma_{\text{NL}}^2), \qquad (7) $$

where $\mathcal{N}(x; m, s^2)$ denotes the normal distribution of mean $m$ and variance $s^2$ over a variable $x$ and $P(x|y)$ the conditional probability of $x$ given $y$. From the measurements of galaxy distance moduli {$\mu_i$} and redshifts {$z_i$}, our goal is to infer the underlying peculiar velocity field $\mathbf{v}$. Since the measurement of extragalactic distances is very noisy (most often about 20 per cent of relative error), we must take into account the correlations between the observed redshifts to extract information on the velocity field. Also, we need to take into account the possible biases appearing when dealing with peculiar velocity field reconstruction, which is done in Section 3.2. Section 3.4 presents how to describe the observations given the above model.

## 3.2 Malmquist bias

Our model is motivated by a rigorous treatment of homogeneous and inhomogeneous Malmquist bias. The homogeneous Malmquist bias is a statistical bias resulting from a combination of the volume effect nearby and the selection effects far away. The observational uncertainties on the measured distances scatter the observed galaxies along the radial direction. Closeby, the number of galaxies grows with the distance, so that overall it is more likely to underestimate the galaxy distances, and as a consequence to assign erroneously positive peculiar velocities (Strauss & Willick 1995). At large distances, the effect is the opposite: it is more likely that the galaxies we observe are scattered away from us, and so in average we overestimate their distances, and assign erroneously negative peculiar velocities. Consequently, neglecting the homogeneous Malmquist bias would create a fake outflow in the central region and a fake inflow on the edges. To overcome this bias, our statistical model fits the underlying most probable luminosity distances of galaxies, together with the velocity field. The bias is statistically handled by attributing a prior function to these distances. The shape of the prior function will be detailed in Section 3.4.4, and will be chosen so that both the volume effect and selection effects are taken into account.

The related inhomogeneous Malmquist bias arises from structure in the observed volume. In the vicinity of a dense region, galaxies are more likely to be scattered from denser to less dense regions. Consequently, the inferred peculiar velocities are biased towards a stronger inflow on to the structures, and thus bias the reconstructed velocity field. By introducing the luminosity distances as parameters, our method is able to reduce this bias. Distances are inferred by both the observations and the reconstructed velocity field, and at each step the galaxies are relocated with respect to the velocity field at their positions. Consequently, an inflow on a dense region will shift the distances towards their true positions. For a more complete analysis of the Malmquist bias, we refer the reader to Strauss & Willick (1995).
### 3.3 Observed radial peculiar velocities

The appropriate modelling of the distribution of errors on the observed peculiar velocities is an another important concern. The peculiar velocities are not directly measured but are usually derived from the distance moduli and redshifts measurements through equation (6), where the luminosity distances *d* are computed from the distance moduli using:

$$ \mu = 5 \log_{10} \frac{d}{10 \text{ pc}}. \qquad (8) $$

Since the errors on the distance moduli are supposed to be Gaussian, the resulting distribution of peculiar velocities will not be Gaussian distributed but rather lognormal distributed (Tully et al. 2016). An example of treatment of this effect is given by Watkins & Feldman (2015) who suggest the use of an unbiased estimator of peculiar velocities with Gaussian distributed errors and which is valid at distances *d* ≳ 20 Mpc.

We use in this paper a different approach. Instead of analysing the observed peculiar velocities, we rather choose to model the distance moduli observations directly. To do so, the introduction of the luminosity distances as free parameters allow us to take into account the Gaussian distribution of distance moduli errors including the relativistic effects of equation (6). The statistical linkage between distance moduli and luminosity distances will be explained in Section 3.4.1.

### 3.4 Statistical model

We aim at recovering the linear peculiar velocity field from the CF3 observations. We work in a Bayesian framework and try to estimate the posterior probability of *v* given the model described above and the observations. We proceed in three steps. First, we detail how the likelihood of our observations *L* is constructed (Sections 3.4.1, 3.4.2, 3.4.3); second we impose priors on the fitted parameters that come from the ΛCDM model presented in the above section (Section 3.4.4). Third, we sample from the posterior probability (Section 3.5).

#### 3.4.1 Distance moduli

Distance indicators used in CF3 are expressed as distance moduli rather than luminosity distances. The link between the two is given by equation (8). Our primary interest in this study concerns deviations from cosmic expansion which are determined independently of the absolute scaling of the extragalactic distance ladder. In addition, our analysis is insensitive to a potential monopole term that might reflect that we live in an overall under or overdense part of the Universe. Tully et al. (2016) argued that the value of the Hubble constant consistent with the 18 000 distances in CF3 is $H_0 = 75$ km s⁻¹ Mpc⁻¹ to within ~3 per cent, discounting the irrelevant absolute scaling. There is agreement at this level in the determination of $H_0$ between the sample within $z = 0.1$ and samples of SNIa at $z \gg 0.1$, limiting concern of a substantial local monopole to flows.

Our conclusion regarding the matter of the Hubble constant is that the value is relatively well determined but uncertainties at the few per cent level remain that are relevant for our analysis. For this reason, we introduce a dimensionless free parameter $h_{\text{eff}}$ to model this uncertainty in the constant:

$$ \mu = 5 \log_{10} \frac{d}{10 \text{ pc}} + 5 \log_{10} h_{\text{eff}}. \qquad (9) $$

If the absolute scale is compatible with the assumed $H_0$, then $h_{\text{eff}}$ should be unity.

Because the distance moduli are very noisy, we need to model the observations by assuming that the error is Gaussian (noted $\sigma_\mu$) and choose to fit for the underlying true luminosity distance *d*:

$$ \mathcal{P}(\mu|d, \sigma_\mu) = \mathcal{N} \left( \mu; 5 \log_{10} \frac{h_{\text{eff}} d}{10 \text{ pc}}, \sigma_\mu^2 \right). \qquad (10) $$

The parameter $h_{\text{eff}}$ is correlated to the reconstructed velocity field and is prone to systematics coming from shifts between zero-point scales of different methodologies, non-linear effects, selection functions, and possible external bulk flows. For this reason, one needs to take this parameter with caution and not as direct measure of $H_0$. For the location of a galaxy in space, the angular position is measured essentially without error. Hence, from the distance *d* and the angular position, we can compute the spatial position of a galaxy, which we denote *r*.

#### 3.4.2 Redshifts

Equation (6) gives the relation between the observed redshift *z* and the cosmological redshift $\bar{z}$, which can be computed from the luminosity distance *d* through equation (2). From these equations, we can compute the radial peculiar velocity of a galaxy:

$$ v^r(z, d) = c \frac{z - \bar{z}(d)}{1 + \bar{z}(d)}. \qquad (11) $$

In CF3, the errors on individual redshifts are not provided, and we chose to assume that they are measured with a Gaussian error of $c\sigma_z = \sigma_{cz} = 50$ km s⁻¹, the typical error for a spectroscopic redshift measurement.

To model the redshifts measurements, we introduce the underlying linear velocity field $v \triangleq v(r_j)_{j \in [0, M^3]}$ sampled on a grid of size $M^3$ and volume $L^3$. Since the sampled velocity field is linear, we need to model the departure of the observed velocities from the linearity. We do so by introducing a Gaussian dispersion $\sigma_{\text{NL}}$ around the linear field which is to be evaluated. The introduction of a unique parameter $\sigma_{\text{NL}}$ hence models the departure of the overall velocity field from linearity and does not model high dispersion inherent to non-linear environment such as clusters of galaxies. For this reason the use of the grouped CF3 catalogue is mandatory. Applying this model on the non-grouped catalogue would underestimate the redshift errors and lead to unphysical results near high-density regions. The probability of observing the redshift *z* knowing the luminosity distance *d* and the velocity field *v* is

$$ \mathcal{P}(z|\mathbf{r}, d, \sigma_z, v) = \mathcal{N} \left( v^r(z, d); v(\mathbf{r}) \cdot \hat{\mathbf{r}}, \sigma_{cz}^2(1 + \bar{z})^{-2} + \sigma_{\text{NL}}^2 \right), \quad (12) $$

where $v^r(z, d)$ is defined by equation (11). Note that $v^r$ is here a function of the parameter *d* and the observation *z*, and does not correspond to a peculiar velocity measurement. The probability distribution above describes the statistical link between the redshifts and the model's parameters. In practice, the linear velocity *v* in equation (12) is computed on a finite-size grid {$r_j$}. To evaluate it at any position (such as the position of a galaxy), we use trilinear interpolation between grid cells. We note from the error distribution of equation (12) that $\sigma_{cz}$ and $\sigma_{\text{NL}}$ are strongly correlated, the two dispersions being different by only the relativistic factor $1 + \bar{z} \simeq 1$. As a consequence, an over (or under) estimation of $\sigma_{cz}$ will be transferred to a under (or over) estimation of $\sigma_{\text{NL}}$, and our results will be insensitive to the choice of $\sigma_{cz}$.
### 3.4.3 Likelihood

The likelihood gives the conditional probability of our observations, namely the observed redshifts $z$ and distance moduli $\mu$, given the model and its parameters, namely the true luminosity distance $d$, the linear velocity field sampled on a grid $v$, $h_{\text{eff}}$, and $\sigma_{\text{NL}}$. Since the conditional probabilities on $\mu$ and $z$ are independent, the likelihood is given by the combination of equations (10) and (12):

$$ 
\begin{aligned}
\mathcal{L} &= \mathcal{P} (\{\mu_i, z_i\} | \{d_i, \sigma_{\mu,i}, \sigma_{z,i}\}, h_{\text{eff}}, \sigma_{\text{NL}}, \mathbf{v}) \\
&= \prod_i \mathcal{P} (\mu_i | d_i, \sigma_{\mu,i}, h_{\text{eff}}) \times \mathcal{P} (z_i | \mathbf{r}_i, d_i, \sigma_z, \mathbf{v}) \\
&= \frac{1}{2\pi} \prod_i \frac{1}{\sqrt{\sigma_{\mu,i}^2}} \exp \left( -\frac{1}{2} \frac{(\mu_i - 5 \log_{10} \frac{h_{\text{eff}} d_i}{10 \text{pc}})^2}{\sigma_{\mu,i}^2} \right) \\
&\quad \times \frac{1}{\sqrt{\sigma_{\varepsilon z}^2 (1 + \bar{z}_i)^{-2} + \sigma_{\text{NL}}^2}} \exp \\
&\quad \times \left( -\frac{1}{2} \frac{(v^r(z_i, d_i) - v(\mathbf{r}_i) \cdot \hat{\mathbf{r}}_i)^2}{\sigma_{\varepsilon z}^2 (1 + \bar{z}_i)^{-2} + \sigma_{\text{NL}}^2} \right),
\end{aligned} 
\quad (13) $$

where the index $i$ denotes the $i$th galaxy or group. The likelihood defined by equation (13) is similar to the one in Lavaux (2016) but simplified because we do not aim to fit the power spectrum properties such as the shape or the normalization. Another difference is that we directly sample the velocity field rather than the overdensity field. We do so to avoid the use of Fourier Transform to compute equation (4) and hence not be subject to periodic boundary conditions. Also, we use a trilinear interpolation to compute the linear velocity field at any point in space, while Lavaux (2016) developed a Fourier-Taylor algorithm to interpolate using only series of Fast Fourier Transforms. We note that our way of proceeding does not enforce the curl-free properties of the linear velocity field. We a posteriori check that the curl part of the reconstructed field is negligible.

From the likelihood (13), we can estimate the probability of a given velocity field $v$ (and associated parameters) from the measurements of distance moduli and redshifts. This probability is given by the Bayes theorem:

$$ \mathcal{P}(\{d_i\}, h_{\text{eff}}, \sigma_{\text{NL}}, v | \{\mu_i, z_i\}) \propto \mathcal{L} \times \mathcal{P}(\{d_i\}) \mathcal{P}(h_{\text{eff}}) \mathcal{P}(\sigma_{\text{NL}}) \mathcal{P}(v), \quad (14) $$

where $\mathcal{P}(\theta)$ denotes the prior on the parameter $\theta$, $\theta \in [(d_i), h_{\text{eff}}, \sigma_{\text{NL}}, v]$. The velocity field reconstruction is then obtained by sampling $v$ from this probability distribution. Before explaining how the sampling is done, we turn now our attention to the definitions of priors.

### 3.4.4 Priors

We assume uniform priors on $h_{\text{eff}}$ and $\sigma_{\text{NL}}$ (see Table 1). Following the model described in Section 3.1, we take into account the correlations between the peculiar velocities by adopting the following prior for $\delta$:

$$ \mathcal{P}(\delta) = \prod_k \frac{1}{\sqrt{2\pi P(k)}} \exp \left( -\frac{|\hat{\delta}(k)|^2}{2P(k)} \right), \quad (15) $$

where $P(k)$ is the power spectrum and has been defined in equation (3). The corresponding prior on the linear velocity field is

$$ \mathcal{P}(\mathbf{v}) = |2\pi\Psi_{\alpha,\beta}|^{-1/2} \exp(-v_{\alpha}\Psi_{\alpha,\beta}^{-1}v_{\beta}), \quad (16) $$

where $(\alpha, \beta)$ denoted Cartesian components and $\Psi_{\alpha,\beta}$ is the velocity-velocity correlation tensor and is defined in the Appendix A. Because of the aforementioned Malmquist biases, the prior on the distances can play a significant role in the analysis. We use two types of empirical priors. The first one is a piecewise normal distribution defined by

$$ \mathcal{P}^{(1)}(d_i | a, b, c) = \frac{1}{\sqrt{2\pi}(b+c)} \begin{cases} \exp\left(-\frac{1}{2} \frac{(d_i-a)^2}{b^2}\right) & \text{if } d_i \le a \\ \exp\left(-\frac{1}{2} \frac{(d_i-a)^2}{c^2}\right) & \text{otherwise,} \end{cases} \quad (17) $$

where $(a, b, c)$ are the shape parameters of the function.

The second one is a power law with an exponential cut-off, as proposed by Lavaux (2016):

$$ \mathcal{P}^{(2)}(d_i | a, b, c) = \frac{1}{N(a, b, c)} (d_i)^a \exp \left( -(d_i/b)^c \right). \quad (18) $$

The normalization factor $N(a, b, c)$ is non-analytical and is computed numerically.

These two priors have the properties that we expect for distance priors: they are bell-shaped curves allowing an asymmetry and with an exponential cut-off at large distance. The shape parameters $(a, b, c)$ allow for some flexibility of these priors and determine the mean value, standard deviation, and skewness of the distributions. Since we are not able to establish the selection function(s) of CF3, we use an empirical approach and choose to fit the shape parameters together with the other parameters of the current model. Leaving these parameters free allows to take into account the volume and selection effects while not imposing strong constraints on derived luminosity distances. In practice, we attribute the prior $\mathcal{P}^{(2)}$ to distances of subsamples (ii), (iii), and (iv) described in Section 2.3, and $\mathcal{P}^{(1)}$ to distances of the 6dF subsample. Distances of subsample $(v)$ are given a uniform prior. We will see in Section 4 that they describe correctly the posterior distribution of distances and in Section 6.1 the effect of changing the prior functions. The model parameters are summarized in Table 1.

## 3.5 Sampling

After we have constructed the likelihood and the priors of the model, we need to sample $v$ from the posterior distribution given by equation (14). In this section, we briefly explain how this sampling is done by a Markov Chain Monte Carlo (MCMC) method. For more technical details, the reader is referred to Lavaux (2016). The sampling method is the partially collapsed blocked Gibbs sampling algorithm. Gibbs sampling is a MCMC method where each parameter is drawn from its conditional probability given the other parameters. Schematically, if we want to sample $n$ parameters $\{x_i\}$, the sampling will be done using the following scheme down to the Markov step $j$:

$$ x_0^1 \leftarrow \mathcal{P} (x_0 | \{x_i^0\}_{i \in \{1:n\}}) $$

$$ x_i^1 \leftarrow \mathcal{P} (x_i | \{x_i^1\}_{i \in \{1:i-1\}}, \{x_i^0\}_{i \in \{i+1:n\}}) $$

$$ x_i^j \leftarrow \mathcal{P} (x_i | \{x_i^j\}_{i \in \{1:i-1\}}, \{x_i^{j-1}\}_{i \in \{i+1:n\}}) $$

Our parameter space is

$$ \Theta = \{h_{\text{eff}}, v, \sigma_{\text{NL}}, \{d_i\}, \{a, b, c\}\} \quad (19) $$

and we need to compute each conditional probability from the likelihood. However, we note that $h_{\text{eff}}$ is strongly correlated with
<table><caption>**Table 1.** Observations and parameters used in this model. For the observations, we specify the number and errors. For the parameters, we specify the corresponding priors.</caption><thead><tr><th>Fixed parameters</th><th>Description</th><th>Number</th><th>Notes and priors</th></tr></thead><tbody><tr><td>N</td><td>Number of galaxy and groups</td><td>1</td><td>N = 11483</td></tr><tr><td>L</td><td>Length of the box side</td><td>1</td><td>L = 800 Mpc h<sup>-1</sup>h<sub>75</sub></td></tr><tr><td>M</td><td>One dimensional size of the grid</td><td>1</td><td>M = 128</td></tr><tr><td>μ</td><td>Distance moduli</td><td>N</td><td>Normally distributed with standard deviation σ<sub>μ</sub>.</td></tr><tr><td>σ<sub>μ</sub></td><td>Errors on distance moduli</td><td>N</td><td>Normally distributed with standard deviation σ<sub>μ</sub></td></tr><tr><td>z</td><td>Observed redshift</td><td>N</td><td>Normally distributed with standard deviation σ<sub>z</sub></td></tr><tr><td>σ<sub>z</sub></td><td>Error on the observed redshifts</td><td>1</td><td>σ<sub>z</sub> = σ<sub>z,c</sub> = 50 km s<sup>-1</sup>c<sup>-1</sup></td></tr><tr><td colspan="4">Free parameters</td></tr><tr><td>v</td><td>The linear velocity field sampled on a grid</td><td>M<sup>3</sup></td><td>ΔCDM prior defined by equation (16)</td></tr><tr><td>d</td><td>Luminosity distances</td><td>N</td><td>Empirical priors defined by equations (17) and (18) and/or uniform, depending on the membership in the subsamples defined in Section 2.2.</td></tr><tr><td>h<sub>eff</sub></td><td>Effective shift of the distance moduli scale</td><td>1</td><td>Uniform prior within [0.5,1.5]</td></tr><tr><td>σ<sub>NL</sub></td><td>Gaussian standard deviation modelling the departure from linearity</td><td>1</td><td>Uniform prior within [50,1500] km s<sup>-1</sup></td></tr><tr><td>(a, b, c)</td><td>Hyperparameters defining the distance priors shapes</td><td>3 × 4</td><td>Uniform priors depending on the prior function</td></tr></tbody></table>

the velocity field, and we draw this parameter from its conditional probability marginalized over the velocity field to make the sampling more efficient. This is called collapsed Gibbs sampling. At the end, our sampling is the following procedure:

(i) We first sample $h_{\text{eff}}$ from the probability distribution marginalized over the velocity field:

$$ \mathcal{P}(h_{\text{eff}}|d, \sigma_{\text{NL}}) = \mathcal{N}(\mathbf{v}^r(h_{\text{eff}}); \mathbf{0}, \mathbf{C}(h_{\text{eff}})), \quad (20) $$

where $\mathbf{v}^r$ is the vector of the galaxies radial peculiar velocities and $\mathbf{C}$ is the velocity autocorrelation matrix defined in the Appendix A. Those two quantities implicitly depend on the parameters $\sigma_{\text{NL}}$ and $d$.

(ii) We draw $\sigma_{\text{NL}}$ from the conditional probability:

$$ \mathcal{P}(\sigma_{\text{NL}}|h_{\text{eff}}, d, \mathbf{v}) = \prod_{i}^{N} \mathcal{N}(\mathbf{v}_{i}^{r}; \mathbf{v}(\mathbf{r}_{i}) \cdot \hat{\mathbf{r}}_{i}, \sigma_{z,i}^{2}(1 + \bar{z}_{i})^{-2} + \sigma_{\text{NL}}^{2}). \quad (21) $$

(iii) We draw a constrained realization of $\mathbf{v}$. This is done by using the Hoffman-Ribak algorithm (Hoffman & Ribak 1991). Appendix B reminds the reader about the algorithm.

(iv) From the sampled constrained realization, we generate a new set of luminosity distances $d$ with the following probability:

$$ \mathcal{P}(d|h_{\text{eff}}, \sigma_{\text{NL}}, \mathbf{v}) \propto \mathcal{L} \times \mathcal{P}(d). \quad (22) $$

(v) Eventually we fit the hyperparameters $(a, b, c)$ of the distance prior functions over the generated distances.

This procedure is carried until convergence. At the end we have a number of realizations of all parameters following the posterior probability defined by equation (14). The reconstructed velocity field is assumed to be the mean of the linear velocity field samples $\langle\mathbf{v}\rangle_{\text{MCMC}}$ and the error is the standard deviation. The method has been tested on different mocks by Lavaux (2016) and we test our specific implementation as described in Appendix C.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** The MCMC for three parameters: (top) the effective reduced Hubble constant $h_{\text{eff}}$; (mid) the non-linear dispersion $\sigma_{\text{NL}}$; (bottom) the reconstructed overdensity field near Virgo, at coordinates (SGX, SGY, SGZ) = (-3.6, 15.6, -0.7) Mpc h<sup>-1</sup>h<sub>75</sub>. The three black lines represent the 15.9th (dashed), 50th (plain), and 84.1th (dashed) percentiles.</figcaption></figure>

## 4 THE COSMICFLOWS-3 PECULIAR VELOCITY FIELD

In this section we present the results of the method applied on CF3 catalogue. We chose to assume ($\Omega_m, \Omega_\Lambda, H_0$) = (0.3, 0.7, 75.0). We fit the distance priors by the empirical function defined by equation (18) for the subsamples (ii, iii, iv) defined in Section 2.3. The distances from the subsample (i) corresponding to 6dF data are fitted assuming the function described by equation (17). We do so to ensure the possibility for the prior function to model a sharp cut-off in distances.

We fit for a global effective zero-point shift $h_{\text{eff}}$ and assume the intercalibration of distances determined by Tully et al. (2016) and the CMB frame transformations given by Tully et al. (2008). Fig. 2 shows the resulting MCMC chain for the two parameters $\sigma_{\text{NL}}$, $h_{\text{eff}}$ and for the overdensity field reconstructed near the Virgo Cluster at location (SGX, SGY, SGZ) = (-3.6, 15.6, -0.7) Mpc h<sup>-1</sup>h<sub>75</sub> and the corresponding histograms are given in Fig. 3. We see
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The posterior distribution for three parameters: (left) the effective reduced Hubble constant $h_{\text{eff}}$; (mid) the non-linear dispersion $\sigma_{\text{NL}}$; (right) the reconstructed overdensity at Virgo, of coordinates $(\text{SGX}, \text{SGY}, \text{SGZ}) = (-3.6, 15.6, -0.7)$ Mpc. The value of $m$ given in the legend corresponds to the mean value of the posterior, while $s$ corresponds to its standard deviation. The three black lines represent the 15.9th (dashed), 50th (plain), and 84.1th (dashed) percentiles.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** The normalized autocorrelation of the MCMC for three parameters: the effective reduced Hubble constant $h_{\text{eff}}$, the non-linear dispersion $\sigma_{\text{NL}}$, and the reconstructed overdensity at Virgo, of coordinates $(\text{SGX}, \text{SGY}, \text{SGZ}) = (-3.6, 15.6, -0.7)$ Mpc. The black solid line corresponds to a null autocorrelation.</figcaption></figure>

from these two figures that the chain has globally converges and result in approximately Gaussian posterior distributions. To further evaluate the convergence of our Markov Chain, we plot in Fig. 4 the normalized autocorrelation of the chains for the two parameters $h_{\text{eff}}$ and $\sigma_{\text{NL}}$. The autocorrelation of a parameter $f$ for a correlation length $\tau$ is defined by

$$c_j(\tau) = \frac{1}{N - \tau} \sum_{i=1}^{N-\tau} (f_i - \bar{f})(f_{i+\tau} - \bar{f}), \qquad (23)$$

where $\bar{f}$ is the mean value of $f$ computed on the $N$ samples. The first decorrelation for a chain corresponds to the intersection with zero. We can see that the $h_{\text{eff}}$ has a decorrelation length of about ~80 and the non-linear dispersion around ~20 steps. This suggests that over our 1400 MC steps, there is some 20 independent samples. For each parameter, the error on the mean of the distribution decreases, as $1/\sqrt{N_{\text{sample}}}$ after convergence and running the chain further will reduce this statistical error. We estimate here that the sampling error on the mean overdensity field (typically ~0.1) is sufficient considering the computation cost of the sampling.

We find $h_{\text{eff}} = 1.02 \pm 0.01$, which suggests that the calibration of CF3 data is compatible with the assumed fiducial Hubble constant $H_0 = 75 \text{ km s}^{-1}\text{Mpc}^{-1}$ within 2 per cent. The non-linear dispersion parameter $\sigma_{\text{NL}}$, which models our lack of knowledge about the non-linear part of the velocity field, is found at $\sigma_{\text{NL}} = 280 \pm 35 \text{ km s}^{-1}$. The fitted value of around $300 \text{ km s}^{-1}$ appears to be high compared to the typical value of $100-200 \text{ km s}^{-1}$. This can be due to underestimation of distance modulus errors, or redshifts errors which were set here at $\sigma_{\text{cz}} = 50 \text{ km s}^{-1}$ for every galaxy. Also, this value depends on the efficiency of the grouping described Section 2.2 at removing entirely the non-linearities in groups of galaxies. Another possibility is that the trilinear interpolation used to compute the reconstructed radial velocity at the position of each galaxy artificially increases the departure from linearity modelled by $\sigma_{\text{NL}}$. Overall this parameter absorbs uncertainties of our model.

Fig. 5 shows the fitted priors on the distances. Overall the shape of the priors is in agreement with the underlying distance distribution. The case of 6dF data might look surprising because there is a discrepancy between the fitted distances and the original ones. However, by considering the distribution in redshifts shown in Fig. 1, we can see that the tail of the measured distance distribution is only due to observational errors and results from the convolution of the real distance distribution with the Gaussian of errors. This particular case illustrates the importance of imposing priors to model selection effects.

Eventually, the $\text{SGZ} = 0 \text{ Mpc } h^{-1}_{75}$ slice of the reconstruction is shown in Fig. 6. In this figure, the colourmap corresponds to the reconstructed overdensity field and the black arrows to the tridimensional linear velocity field. We stress that the reconstructed overdensity and velocity fields are not computed one from another using equation (5), but rather sampled simultaneously.<fnref n="2" /> The original field was computed in a box of width $800 \text{ Mpc } h^{-1}_{75}$, and we show here the central $500 \text{ Mpc } h^{-1}_{75}$ where the signal over noise is non-zero. We can see that the overdensity goes toward zero near the edges, because there are no observed galaxies at large distances. By eye, we can identify several cosmic structures that we will shortly list in Section 5. We also qualitatively see the match between the

<footnote n="2">²For each MC step, we draw both a constrained realization of the velocity field and overdensity field from a common random realization with equations (B2) and (B4).</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Histogram of measured (blue) and fitted (orange) distances for the five subsamples described in Section 2. The green lines are the fitted priors. From top to bottom: (1) 6dF sample, mixed Gaussian prior distribution; (2) Spitzer galaxies, empirical prior distribution; (3) TF by Arecibo Telescope, empirical prior distribution; (4) TF not covered by Arecibo Telescope, empirical prior distribution. The empirical prior distribution is described in equation (18).</figcaption></figure>

overdense structures and the velocity infalls. These fields come with statistical errors that can be computed from the resulting Markov Chain. The standard deviation of the reconstructed fields $\delta$ and $v'$ are plotted in Fig. 7. At large distance, we recover the $\Lambda$CDM standard deviation of the overdensity and velocity field. In particular, we recover the value of $300\,\text{km s}^{-1}$ deviation for the velocity field. This is due to the absence of data points beyond $z = 0.054$, a limit which is illustrated by the black circle. We can see that the radial peculiar velocity field seems less noisy than the overdensity field. The reason is that the velocity field is more correlated than the overdensity as can be seen with equation (A3). Consequently, the root mean square value of the posterior distribution does not capture entirely the correlated errors of the velocity field, and these errors correlate on larger scales. In the next section, we study the reconstructed overdensity field shown in the left-hand panel of Fig. 6 and compare it with the actual distribution of galaxies to check the consistency with other observables.

# 5 COSMOGRAPHY OVERVIEW

In this section we show partial results on the reconstructed overdensity field. A full description will be made in an upcoming

cosmography paper. Main results are presented in Fig. 8 where we plot the overdensity field in three different slices corresponding to coordinates $SGZ = 0\,\text{Mpc h}^{-1}h_{75}$ (twice), $SGX = 0\,\text{Mpc h}^{-1}h_{75}$, and $SGY = -93\,\text{Mpc h}^{-1}h_{75}$ going from top left to bottom right. We first look at the two top panels. On the left, the coloured dots represent CF3 data within a slab of $10\,\text{Mpc h}^{-1}h_{75}$ thickness, while on the right the dark dots are all available galaxies of the LEDA data base located at their redshift positions within the same depth. The coloured galaxies show the anisotropy of CF3 catalogue. The Southern hemisphere is mainly covered by 6dF galaxies (in red), while there are fewer data points in the northern hemisphere. Looking at the right-hand panel, we can observe a good agreement between the reconstructed overdensity field and the location of galaxies. At large distances, we notice small discrepancies (for example around $SGY = -210\,\text{Mpc h}^{-1}h_{75}$, top right-hand panel). The reconstructed structures seem to be shifted compared to the redshift positions. This issue is tightly linked to the recovered value of $h_{\text{eff}}$, which is prone to systematics. This might be a hint of a shift between zero-point calibrations depending on the region covered. We discuss this issue in Section 6.1. We add to the plots the name of some known structures, such as Coma, Shapley, Perseus-Pisces, Apus, and Pisces-Cetus. The two bottom panels correspond to other slices of the Local Universe reconstruction. Again, there is agreement between the reconstructed overdense regions and the galaxy distribution. The bottom right-hand panel exhibits some distant structures which are less known, such as the Southern Wall, Telescopium, or Lepus.

# 6 DISCUSSION

In this paper, we described our linear peculiar velocity field reconstruction method. We then applied it to the CF3 catalogue, considering generic modelling of the prior distance distributions and potential shift on the zero-point calibration of all distances. We saw in Section 5 how this reconstruction can be used to study our local environment and the unbiased distribution of dark matter in the $\Lambda$CDM framework. We now turn to the discussion of the limits and possible improvements of the presented method.

## 6.1 The Hubble constant

We introduced an effective reduced Hubble constant $h_{\text{eff}}$ to absorb uncertainty on the calibration of the distance indicators and assumed a Hubble constant of $H_0 = 75\,\text{km s}^{-1}\text{Mpc}^{-1}$. It is worth noticing that the value of the effective Hubble constant is correlated with our choice of distance priors. This systematic is studied in Appendix D on a subsample of the CF3 catalogue. We find that imposing uniform priors on the distances lowers the value of $h_{\text{eff}}$ by around 5 per cent (equivalent of around 3 km s$^{-1}$Mpc$^{-1}$). This suggests that the parameter $h_{\text{eff}}$ not only absorbs uncertainties on the assumed Hubble constant, but also compensate for residual homogeneous Malmquist bias that is not modelled by the simple shape of our distance priors.

The CF3 catalogue contains distances from six distinct methodologies. With each of these methodologies there are contributions from multiple sources. It is fundamentally important that the distances from the diverse inputs be consistent in zero-point and free of systematics with redshift. An important feature in the construction of CF3 was large overlaps between sources. The complex interlacing is discussed by Tully et al. (2016). The grounding assumption was all measurements of the distance to a galaxy or a cluster of galaxies should give the same value on average. The scale is
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Central part of the CF3 velocity field reconstruction in the SGZ = 0 Mpc h⁻¹₇₅ slice. The colour corresponds to the value of the overdensity field while the black arrows represent the tridimensional reconstructed linear peculiar velocity field. The dotted black circle illustrates the edge of the data at z = 0.054.</figcaption></figure>

bootstrapped from fundamental Cepheid and RR Lyrae calibrators
in our Galaxy. Consequently the CF3 product is distances in Mpc
derived independently from knowledge of velocities.

There is one relatively weak link. The 6dF sample is a major
component of CF3 that explores a domain that is poorly covered by
other samples. The scale linkage is established through 84 individual
galaxies with alternate measures and 381 groups with distances to
members from alternate methods. Given the importance of the 6dF
contribution, in Appendix E we consider the possibility of a scale
mismatch by introducing an independent $h_{\text{eff}}$ parameter associated
with the 6dF contribution. A small difference is found between
the optimal values for the 6dF component and the rest of the CF3
contributions, namely the effective Hubble constant is reduced by
3 per cent for the 6DF sample. We find that the reconstructed radial
velocity field is affected by the introduction of this additional
parameter, and that the non-linear dispersion parameter $\sigma_{\text{NL}}$ is
reduced by 20 per cent. It is not clear, though, that the addition of
a parameter has improved the model. The additional $h_{\text{eff}}$ parameter
for the 6dF component could mask velocity streaming specific to
the sector uniquely sampled by 6dF. Our primary interest is to map

the velocity field so for that purpose the conservative assumption
is to not add an additional parameter. Rather, we assume that the
6dF distances were properly linked to the other CF3 components,
and model departures from the fiducial $H_0$ are described by a single
value of $h_{\text{eff}}$.

## 6.2 Non-linearities

The modelling of non-linearities can also be a source of systematics
of our reconstruction. Our treatment involves the grouping described
in Section 2.2. The merging of data within groups has clear
advantages. While what we directly observe are individual galaxies,
the test particles we want to follow in a linear reconstruction are
collapsed haloes. If a halo contains multiple galaxies then we do
best to average over the constituent properties. The groups are
formulated from a large redshift catalogue (Tully 2015) and we
average over the redshifts of all associated members. Then we find
weighted averages of the distances from those group members in
CF3. Distances of groups with many contributions can have low
formal errors.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Statistical standard deviation for the sampled overdensity (left) and velocity (right) fields in the SGZ = 0 Mpc h⁻¹₇₅ slice over the full reconstructed box. The standard deviation of the velocity field takes into account all three Cartesian components. The black circle represents the edge of the data at z = 0.054.</figcaption></figure>

The remaining non-linearities (outside of the groups) at redshift zero are modelled by a unique Gaussian dispersion to the observed velocities. The width of this distribution, namely σNL, is fitted with the other parameters. This is an approximation: non-linear motions are correlated down to k ≳ 0.1 Mpc⁻¹ scale and so are not randomly independent. However, including non-linear correlations within the actual method is a technical challenge since the correlations in the non-linear regime are not isotropic anymore. An interesting alternative is proposed by Hoffman et al. (2018) who try to use constrained simulations in order to reconstruct the non-linear part of the peculiar velocity field. Other options could also be tried like fixing the σNL parameter at a value constrained by ΛCDM non-linear simulation, or adding a discrete parameter for each galaxies to model their association to non-linear regions, such as in Lavaux (2016). The advantage of the latter method is to model both non-linearities and possible outliers.

## 6.3 Selection effects

Selection effects on distances play a considerable role in peculiar velocity analyses even though less critical than in galaxy number counts analysis. We described in Section 3.2 the Malmquist biases appearing when the selection effects are not properly modelled. In our Bayesian approach, the selection effects are modelled by constructing the probability P(d) of a galaxy having a luminosity distance d knowing that we observed it. Our model suppose that this probability only depends on the actual distance d and fitted hyperparameters (a, b, c). In particular, individual CF3 distances are unbiased (Tully et al. 2016), and the probability P(d) does not depend on the observed distance modulus μ.

We note however that selection effects would need further investigation in future works on peculiar velocities. As suggested by Hinton, Kim & Davis (2017), a proper model needs the introduction of the generalized likelihood L':

$$ \mathcal{L}' = \mathcal{L} \times \frac{\mathcal{P}(\text{selection}|\text{data, parameters})}{\int \mathcal{P}(\text{selection}|\mathcal{D}, \text{parameters})\mathcal{P}(\mathcal{D}|\text{parameters})d\mathcal{D}}, \quad (24) $$

where the integration is over all possible observational data D. Since the selection function of CF3 is not analytical it is challenging to write the term P选选选|D, parameters) representing the probability of selecting a galaxy measurement given its overall properties and the model's free parameters. Sampling from the extended likelihood L' is consequently hard and impractical for the current analysis.

Our approach, while it is approximate, is robust and secure since we fit the priors directly on the reconstructed distances. Doing so, we value CF3 distances over eventual prior information. Also, it is worth noticing that the CF3 catalogue benefits from the multiplicity of methods included in it. The selection effects being different from one methodology to another, the overall reconstructed velocity field should not be subject to individual method specificity. As mentioned in Section 6.1, we tested the effect of changing the priors on the distances and the results are presented in Appendix D.

## 7 CONCLUSIONS

This article presents an algorithm to reconstruct the linear peculiar velocity field up to z ~ 0.054 from the CF3 catalogue. We have been able to reconstruct for the first time the Local Universe velocity field from CF3 data, and showed some results in terms of cosmography. We also associated a corresponding map of statistical errors on both the overdensity and velocity field. The reconstructed field will be used for both cosmological analysis and cosmography. We have stressed the limits of the current method, especially about the model of non-linearities and selection effects. We showed that these effects were not completely negligible and should be more precisely modelled in future analyses.

Above all, this article highlights the ability of peculiar velocities to probe the matter distribution. With only about 11 000 tracers, we were able to map and identify overdense and underdense regions in the Local Universe, and showed the good agreement with redshift surveys containing more than 150 000 galaxies. As suggested by Lavaux (2016), the model could be extended to estimate cosmological parameters such as the Hubble constant H₀ or the growth rate fσ₈. With upcoming large distance data sets, coming from TAIPAN (da Cunha et al. 2017, ~50 000 FP distances),
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Four slices of the reconstructed overdensity field of CF3. The coloured dots represent galaxies from CF3 catalogue. The colour depends on the appartenance of galaxies to subsamples defined in Section 2.3 and are coded the same way as in Fig. 1: (red) 6dF galaxies, subsample (i); (purple) Spitzer TF data, subsample (ii); (green) Arecibo TF data, subsample (iii); (orange) Other TF data, subsample (iv); (brown) Remaining individual and groups data, subsample (v). On other slices, the dark dot represents galaxies from LEDA located at their redshift position. Coordinates of slices: top left and top right SGZ = 0 Mpc $h^{-1}_{75}$; Bottom left SGX = 0 Mpc $h^{-1}_{75}$; Bottom right: SGY = -93 Mpc $h^{-1}_{75}$. We denominate major structures on the different shown slices.</figcaption></figure>

WALLABY (Duffy et al. 2012, ~60 000 TF distances), and LSST (Lochner et al. 2018, ~100 000 SNIa distances), peculiar velocity analyses will need an accurate model to avoid systematics in the determination of cosmological parameters. The method presented here is to be considered as a baseline of such model.

has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 759194 - USNAC, Understanding Type Ia SuperNovae for Accurate Cosmology).

## ACKNOWLEDGEMENTS

RG thanks Mickael Rigault for fruitful discussions. We are grateful to the colleagues maintaining LEDA Lyon extragalactic data base Dmitry Makarov and Philippe Prugniel. This research was supported by Institut Universitaire de France and CNES. This project

## REFERENCES

Courtois H. M., Hoffman Y., Tully R. B., Gottlöber S., 2012, ApJ, 744, 43

Courtois H. M., Tully R. B., Hoffman Y., Pomarède D., Graziani R., Dupuy A., 2017, ApJ, 847, L6

Da Costa G. S., Armandroff T. E., 1990, AJ, 100, 162

da Cunha E. et al., 2017, PASA, 34, e047
Davis T. M., Scrimgeour M. I., 2014, MNRAS, 442, 1117

Dekel A., Eldar A., Kolatt T., Yahil A., Willick J. A., Faber S. M., Courteau S., Burstein D., 1999, ApJ, 522, 1

Djorgovski S., Davis M., 1987, ApJ, 313, 59

Dressler A., Lynden-Bell D., Burstein D., Davies R. L., Faber S. M., Terlevich R., Wegner G., 1987, ApJ, 313, 42

Duffy A. R., Meyer M. J., Staveley-Smith L., Bernyk M., Croton D. J., Koribalski B. S., Gerstmann D., Westerlund S., 2012, MNRAS, 426, 3385

Feix M., Branchini E., Nusser A., 2017, MNRAS, 468, 1420

Gorski K. M., Davis M., Strauss M. A., White S. D. M., Yahil A., 1989, ApJ, 344, 1

Gottlieber S., Hoffman Y., Yepes G., 2010, preprint (arXiv:1005.2687)

Hinton S. R., Kim A., Davis T. M., 2017, preprint (arXiv:1706.03856)

Hoffman Y. et al., 2018, Nature Astron., 2, 680

Hoffman Y., Ribak E., 1991, ApJ, 380, L5

Hoffman Y., Courtois H. M., Tully R. B., 2015, MNRAS, 449, 4494

Hoffman Y., Pomarède D., Tully R. B., Courtois H. M., 2017, Nature Astron., 1, 0036

Hong T. et al., 2014, MNRAS, 445, 402

Howlett C. et al., 2017, MNRAS, 471, 3135

Huchra J. P. et al., 2012, ApJS, 199, 26

Hudson M. J., Smith R. J., Lucey J. R., Branchini E., 2004, MNRAS, 352, 61

Lavaux G., 2016, MNRAS, 457, 172

Leavitt H. S., Pickering E. C., 1912, Harv. Coll. Obs. Circ., 173, 1

Lochner M. et al., 2018, preprint (arXiv:1812.00515)

Nusser A., 2017, MNRAS, 470, 445

Nusser A., Davis M., 2011, ApJ, 736, 93

Phillips M. M., 1993, ApJ, 413, L105

Planck Collaboration XIII, 2015, A&A, 594, A13

Scrimgeour M. I. et al., 2016, MNRAS, 455, 386

Sorce J. G. et al., 2016, MNRAS, 455, 2078

Springob C. M. et al., 2014, MNRAS, 445, 2677

Springob C. M., Masters K. L., Haynes M. P., Giovanelli R., Marinoni C., 2007, ApJS, 172, 599

Strauss M. A., Willick J. A., 1995, Phys. Rep., 261, 271

Tonry J., Schneider D. P., 1988, AJ, 96, 807

Tully R. B. et al., 2013, AJ, 146, 86

Tully R. B., 2015, AJ, 149, 171

Tully R. B., Fisher J. R., 1977, A&A, 54, 661

Tully R. B., Shaya E. J., Karachentsev I. D., Courtois H. M., Kocevski D. D., Rizzi L., Peel A., 2008, ApJ, 676, 184

Tully R. B., Courtois H., Hoffman Y., Pomarède D., 2014, Nature, 513, 71

Tully R. B., Courtois H. M., Sorce J. G., 2016, AJ, 152, 50

Wang Y., Rooney C., Feldman H. A., Watkins R., 2018, MNRAS, 480, 5332

Watkins R., Feldman H. A., 2015, MNRAS, 450, 1868

Zaroubi S., Zehavi I., Dekel A., Hoffman Y., Kolatt T., 1997, ApJ, 486, 21

Zaroubi S., Hoffman Y., Dekel A., 1999, ApJ, 520, 413

## APPENDIX A: LINEAR PECULIAR VELOCITY DISTRIBUTION IN THE λCDM MODEL

Our framework is the linear theory of structure formation within a flat ΛCDM model. The initial perturbations of such universe are characterized by their Gaussian statistics as observed from the CMB (Planck Collaboration XIII 2015). In the linear theory, the Fourier modes of these perturbations has grown independently and the overdensity field of our Local Universe can still be described by its power spectrum $P(k)$:

$$ \langle \delta(\mathbf{k}) \delta^*(\mathbf{k}') \rangle = (2\pi)^3 \delta_D(\mathbf{k} - \mathbf{k}') P(k), \quad (A1) $$

where $\delta(\mathbf{k})$ is the Fourier transform of $\delta(\mathbf{r})$. Such a linear approximation is only valid at large scales, i.e. $k \lesssim 0.2 h \, \text{Mpc}^{-1}$, and consequently this analysis aims at recovering only the very large-scale structures, down to scales of few tens of Mpc. Given the

cosmological growth rate $f$ ( depends on the adopted cosmology), one can compute the velocity field $v(\mathbf{r})$ by

$$ \nabla \cdot \mathbf{v} = -H_0 f \delta, \quad (A2) $$

which can be written in Fourier space as

$$ \mathbf{v}(\mathbf{k}) = i H_0 f \frac{\mathbf{k}}{k^2} \delta(\mathbf{k}). \quad (A3) $$

Consequently, the velocity–velocity two-point correlation tensor is, in configuration space:

$$ \Psi_{\alpha,\beta}(\mathbf{r}) \triangleq \langle v(\mathbf{r}')v(\mathbf{r}'+\mathbf{r}) \rangle_{\alpha\beta} = \frac{(H_0 f)^2}{(2\pi)^3} \int_0^\infty \frac{k_\alpha k_\beta}{k^4} P(k) e^{-ik \cdot r} dk. \quad (A4) $$

In practice, $\Psi_{\alpha,\beta}(\mathbf{r})$ can be expressed using the radial and transverse correlations functions $\psi_R$ and $\psi_T$ (Gorski et al. 1989):

$$ \Psi_{\alpha,\beta}(\mathbf{r}) = \psi_T(\mathbf{r}) \delta_{\alpha,\beta} + (\psi_R(\mathbf{r}) - \psi_T(\mathbf{r})) \hat{r}_\alpha \hat{r}_\beta \quad (A5) $$

$$ \psi_R(x) \triangleq \frac{1}{2\pi^2} \int_0^\infty \left( j_0(kx) - \frac{2j_1(kx)}{kx} \right) P(k) dk \quad (A6) $$

$$ \psi_T(x) \triangleq \frac{1}{2\pi^2} \int_0^\infty \frac{j_1(kx)}{kx} P(k) dk \quad (A7) $$

where $j_0$ and $j_1$ are the zeroth and first-order spherical Bessel functions and $\delta^k$ is the Kroenecker delta. In practice, we use tabulated $\psi_R$ and $\psi_T$, and use linear interpolation between the sampled positions. We also define the covariance matrix $\mathbf{C}$ of a set of radial peculiar velocities by:

$$ [\mathbf{C}]_{i,j} \triangleq \langle v_i^r v_j^r \rangle + \langle \epsilon_i \epsilon_j \rangle = \sum_{\alpha,\beta} \hat{r}_{i,\alpha} \hat{r}_{j,\beta} \Psi_{\alpha,\beta} + \langle \epsilon_i \epsilon_j \rangle. \quad (A8) $$

Usually, the matrix of errors $\langle \epsilon_i \epsilon_j \rangle$ is taken as the sum of the error on the redshift measurement plus a dispersion due to non-linearities at $z \sim 0$, $\sigma_{\text{NL}}$:

$$ \langle \epsilon_i \epsilon_j \rangle = \delta_{ij}^K (\sigma_{cz}^2 + \sigma_{\text{NL}}^2). \quad (A9) $$

## APPENDIX B: THE HOFFMAN–RIBAK ALGORITHM

The marginalized probability density for the overdensity field $P(\{\delta(\mathbf{r}_j)\}|\{d_i\}, \sigma_{\text{NL}}, h_{\text{eff}})$ is

$$ \begin{aligned} P(\delta|\{d_i\}, \sigma_{\text{NL}}, h_{\text{eff}}) = & \propto \prod_i \frac{1}{\sqrt{\sigma_{cz}^2(1 + \bar{z}_i)^{-2} + \sigma_{\text{NL}}^2}} \exp \\ & \times \left( -\frac{1}{2} \frac{(v^r(z_i, \bar{z}_i, h_{\text{eff}}d_i)) - v \cdot \hat{r}_i}{\sigma_{cz}^2(1 + \bar{z}_i)^{-2} + \sigma_{\text{NL}}^2} \right) \\ & \times \prod_j \frac{1}{\sqrt{2\pi P(k_j)}} \exp \left( -\frac{|\delta(\mathbf{k}_j)|^2}{2P(k_j)} \right). \end{aligned} \quad (B1) $$

To sample $\delta$ from this probability density function Hoffman & Ribak (1991) proposes the following. From the power spectrum $P(k)$ is generated a random realization $\delta^{\text{RR}}$. Then, a constrained realization is computed with the previous random part and a correlated one:

$$ \delta^{\text{CR}} = \delta^{\text{RR}} + \langle \delta c_i \rangle \langle c_i c_j \rangle^{-1} (c_i - c_i^{\text{RR}}), \quad (B2) $$
where $c_i$ is a constraint on the sampled field, here the radial peculiar velocities $v^r$. The matrix $\langle c_i c_j \rangle = C$ is defined in Appendix A, and the correlation between the overdensity field and the radial peculiar velocity is

$$ \langle \delta c_i \rangle = \frac{H_0 f}{(2\pi)^3} \int \frac{\mathbf{k} \cdot \hat{\mathbf{r}}}{k^2} P(k) e^{-ik \cdot \mathbf{r}} d\mathbf{k}. \qquad (B3) $$

It is also possible to directly sample the velocity field using:

$$ v_{\alpha}^{\text{CR}} = v_{\alpha}^{\text{RR}} + \langle v_{\alpha} c_i \rangle \langle c_i c_j \rangle^{-1} (c_i - c_i^{\text{RR}}), \qquad (B4) $$

where the index $\alpha$ corresponds to the Cartesian component. The correlation $\langle v_{\alpha} c_i \rangle$ is given by

$$ \langle v_{\alpha} c_i \rangle = \sum_{\beta} \Psi_{\alpha, \beta} r_{\beta}. \qquad (B5) $$

In this paper, we sample both the velocity and density field using the same random realization. We thus do not need the periodic boundary conditions implied by the use of fast Fourier transform to evaluate equation (4).

## APPENDIX C: TEST ON MOCK

We test the implementation of the algorithm on a mock catalogue of 4000 tracers generated as follows:

(i) The angular positions are drawn from an uniform distribution,

(ii) The distances are drawn from a truncated normal distribution within $d \in [0, 200]$ Mpc,

(iii) We mimic the observations by computing the corresponding distance moduli and scatter them following a normal distribution of standard deviation $\sigma_{\mu} = 0.2$,

<figure><img src="image_11.png" /><figcaption>Figure C1. MCMC chain for the mock. The upper panel shows the $h_{\text{eff}}$ parameter while the bottom one shows the non-linear dispersion $\sigma_{\text{NL}}$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure C2. Comparison between the original radial peculiar velocities and the reconstructed ones. The left-hand panel shows the residuals $\frac{v^r - v_{\text{RR}}}{\sigma_{v^r}}$ to the original field in the $z=0$ slice. The right-hand panel shows the histogram of this quantity over the whole box. The black solid line is the unit Gaussian.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure C3. Same as Fig. C2 for the reconstruction using WF/CR technique alone.</figcaption></figure>

(iv) We input a shift in the distance moduli scale of $h_{\text{eff}} = 1.0$.

(v) From a random realization generated from a power spectrum truncated at 0.1 Mpc (to avoid non-linearities), we draw the peculiar velocities. We add a non-linear dispersion with $\sigma_{\text{NL}} = 150 \text{ km s}^{-1}$,

(vi) From the original distances and peculiar velocities, we compute the measured redshifts and add a Gaussian dispersion of $\sigma_{cZ} = 50 \text{ km s}^{-1}$,

From the simulated distance moduli and redshifts, we reconstruct the velocity field following the procedure described in Section 3.5. We compute it on a grid of size 64³ and box of 500 Mpc width. Fig. C1 shows the resulting distribution for the two parameters $h_{\text{eff}}$ and $\sigma_{\text{NL}}$, which both agree with their fiducial values. The comparison between the original Gaussian random field and the reconstructed radial velocity field is shown in Fig. C2. The left-hand panel shows the pull distribution of the reconstructed radial velocity field in the $Z=0$ slice. On the left plot, we show the histogram of these values within the white circle corresponding to the data limit. Overall the distribution is close to a unit Gaussian, showing the unbiased aspect of the reconstruction. Also, on the left-hand panel, there is no clear sign of reconstructed structures over more than $3\sigma$. For comparison, the same plot is done using the WF/CR technique alone, i.e. fixing the parameters $d$, $\sigma_{\text{NL}}$, and $h_{\text{eff}}$, is shown in Fig. C3. We can see the strong improvement from the method. Homogeneous Malmquist bias affects the distances so that there is excessive outflow in the center of the box and inflow outside, as it is described in Section 3.2.

## APPENDIX D: CHANGING THE DISTANCE PRIORS

We test the effect of changing the assumed priors on distances on a subsample of the CF3 catalogue. The test catalogue is built by

<figure><img src="image_14.png" /><figcaption>Figure D1. The posterior distribution for three parameters in the case of a reconstruction considering uniform priors en distances (except for 6dF data, see text): (left) the effective reduced Hubble constant $h_{\text{eff}}$; (mid) the non-linear dispersion $\sigma_{\text{NL}}$; (right) the reconstructed radial velocity at Virgo, of coordinates $(SGX, SGY, SGZ) = (-3.6, 15.6, -0.7)$ Mpc.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure D2.** Same as Fig. 6 in the case of three different reconstructions from the half CF3 catalogue (see text). Left: Nominal reconstruction; mid: reconstruction with uniform priors except for 6dF subsample; right: reconstruction with an independent zero-point shift for 6dF subsample.</figcaption></figure>

randomly taking half of the galaxies in each subsample defined in Section 2.3. The reconstruction is done on a grid of size 64³. The reconstruction was carried in the case of uniform priors imposed on all distances except the 6dF one. We treat 6dF as separate because the sharp cut-off on their redshift makes the analysis unrealistic with uniform prior. On 6dF, we fit an empirical prior defined by equation (18). The calculation was carried on 800 MCMC steps and the first 300 are considered as the warm-up phase. The Fig. D1 shows the resulting histograms for the $h_{\text{eff}}$ and $\sigma_{\text{NL}}$ parameters. One can see that the effective Hubble constant has been minored compared to the reconstruction shown in Fig. 3 by more than one sigma. This shows the strong dependance of any Hubble constant determination with this method on the prior distance distributions. Fig. D2 shows the comparison between the reconstructed fields in the case of a nominal reconstruction (empirical priors and unique shift in the zero-point) (left) and the case of uniform priors (middle). The differences are small, but one can notice the inflow on the $SGX > 0 \text{ Mpc } h^{-1}_{75}$ part (data which is mainly not covered by 6dF data) has been increased. Because we imposed no priors on the data, the distances are more likely to be overestimated far away, and consequently the radial velocities underestimated, biasing the field towards negative values. This results shows that the prior distribution has an impact in our Bayesian analysis and underlines the fact that selection effects have to be properly modelled to extract cosmological parameters from the velocity field reconstruction.

## APPENDIX E: CONSIDERING 6DF AS A SEPARATE SUBSAMPLE

We try to fit two different reduced effective Hubble constants $h_{\text{eff}}$, one for 6dF data and the other for the rest of the data. The distance priors are kept to empirical prior functions defined by equation (18). The resulting histograms are shown in Fig. E1. We see that the values recovered for the two $h_{\text{eff}}$ parameters are different while still compatible. We can see that the non-linear dispersion has decreased compared to the main and uniform priors reconstructions ($\sim 225 \text{ km s}^{-1}$ compared to $\sim 270 \text{ km s}^{-1}$). This suggests a better fit of the underlying velocity field. However, it could also be due to the

<figure><img src="image_16.png" /><figcaption>**Figure E1.** The posterior distribution for three parameters in the case of a reconstruction with two different effective Hubble constants (see text): (left) the effective reduced Hubble constant $h_{\text{eff}}$ for 6dF data; (mid) $h_{\text{eff}}$ for the other galaxies; (right) the non-linear dispersion $\sigma_{\text{NL}}$.</figcaption></figure>

artificial reduction of a real flow between 6dF covered region and the rest, which is why the main result of this article assumes a unique zero-point shift. Again, Fig. D2 shows the comparison between the reconstructed fields in the case of a fiducial reconstruction (empirical priors and unique shift in the zero-point) (left) and the case of two different $h_{\text{eff}}$ (right). The differences are more pronounced: The inflows in the regions $SGX > 100 \text{ Mpc } h^{-1}_{75}$ and $SGX < 200 \text{ Mpc } h^{-1}_{75}$ have been reduced. Also, we can see from the overdensity field that the overall inflows on structures have changed between the left-hand and right-hand panels, increasing in the case of a unique zero-point. We stress that the difference of flows can be artificial and the conservative approach taken in the main result of this paper (empirical priors and unique zero-point shift) should be preferred. However, this test shows that small shifts in distance indicators calibrations change the mean inflow on structures, and could have an impact on the determination of cosmological parameters, such as the growth rate $f\sigma_8$ or the Hubble constant $H_0$. This impact is however beyond the work presented in this paper. Extending the model to the determination of cosmological papers from peculiar velocities is to be considered as a possible development of the current methodology.

This paper has been typeset from a TeX/Latex file prepared by the author.