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

The reconstruction is based on the Cosmicflows-3 (CF3) catalog (Tully et al., 2016)¹¹1Available as a table at http://edd.ifa.hawaii.edu. 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; Springob et al., 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 catalog. Nearby, coverage is dense and distances are accurate at the level of 5%. 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.

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 modeled within our $\Lambda$CDM linear framework. Our solution is to average information over the small scale of groups. Tully (2015) provides a catalog 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 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 catalog.

Our Bayesian methodology (see Section 3) needs priors on the statistical distribution of distances. Because it is a composite catalog, 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 5 main subsamples based on the original observational surveys:

1. 1.

   6dFGSv data provides the most well defined subsample: it has a high degree of completeness up to a sharp cutoff at $z=0.054$. The subsample contains 5,777 galaxies or groups with a median redshift of $z=0.039$.

2. 2.

   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 1,546 galaxies or groups included in this subsample have a median redshift of $z=0.009$.

3. 3.

   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 1,628 galaxies or groups with a median redshift of $z=0.023$.

4. 4.

   Other TF data than in the Arecibo declination range. The subsample contains 1,569 galaxies or groups with a median redshift of $z=0.019$.

5. 5.

   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 catalog is $N=11,483$. Figure 1 shows the redshift histograms of all these subsamples. We see that 6dF data plays a major role in CF3 and has a singular behavior with a sharp cutoff 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.

The methodology relies on the assumption of a flat $\Lambda$CDM model parametrized by the set of parameters $(\Omega_{m},H_{0})$. The $\Lambda$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:

where $\Omega_{\Lambda}=1-\Omega_{m}$. Eq. 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(\textbf{{r}})$ is small $|\delta(\textbf{{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)$:

where $\delta(\textbf{{k}})$ denotes the Fourier transform of $\delta(\textbf{{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 v and density fields are related by:

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

Consequently, the statistic of v is also Gaussian. The statistic of v is described in the Appendix A.

Because of the peculiar velocity field v, the Hubble flow described in Eq. 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:

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

where $\mathcal{N}(x;m,s^{2})$ denotes the normal distribution of mean $m$ and variance $s^{2}$ over a variable $x$ and $\mathcal{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 v. Since the measurement of extragalactic distances is very noisy (most often about 20$\%$ 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. The Section 3.4 present how to describe the observations given the above model.

Our model is motivated by a rigourous 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 toward a stronger inflow onto 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 toward their true positions. For a more complete analysis of the Malmquist bias, we refer the reader to  Strauss & Willick (1995).

The appropriate modeling 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 the Eq. (6), where the luminosity distances $d$ are computed from the distance moduli using:

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 Log-normal 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\gtrapprox 20$ Mpc.

We use in this paper a different approach. Instead of analyzing 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 Eq. (6). The statistical linkage between distance moduli and luminosity distances will be explained in Section 3.4.1.

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 $\mathcal{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 $\Lambda$CDM model presented in the above section (Section 3.4.4). Third, we sample from from the posterior probability (Section 3.5).

After we have constructed the likelihood and the priors of the model, we need to sample v from the posterior distribution given by Eq. (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$:

Our parameter space is:

and we need to compute each conditional probability from the likelihood. However, we note that $h_{\text{eff}}$ is strongly correlated with 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:

1. 1.

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

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

   (20)

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

2. 2.

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

   $$\mathcal{P}\left(\sigma_{\text{NL}}|h_{\text{eff}},d,\textbf{{v}}\right)=\prod^{N}_{i}\mathcal{N}\left(v^{r}_{i};\textbf{{v}}(\textbf{{r}}_{i})\cdot\hat{r}_{i},\sigma_{z,i}^{2}(1+\bar{z}_{i})^{-2}+\sigma^{2}_{NL}\right)$$

   (21)

3. 3.

   We draw a constrained realization of v. This is done by using the Hoffman-Ribak algorithm (Hoffman & Ribak, 1991). Appendix B reminds the reader about the algorithm.

4. 4.

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

   $$\mathcal{P}\left(d|h_{\text{eff}},\sigma_{\text{NL}},\textbf{{v}}\right)\propto\mathcal{L}\times\mathcal{P}\left(d\right)$$

   (22)

5. 5.

   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 Eq. (14). The reconstructed velocity field is assumed to be the mean of the linear velocity field samples $\langle\textbf{{v}}\rangle_{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.

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$ km/s/Mpc. 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 catalog. We find that imposing uniform priors on the distances lowers the value of $h_{\text{eff}}$ by around $5\%$ (equivalent of around 3 km/s/Mpc). This suggests that the parameter $h_{\text{eff}}$ not only absorb uncertainties on the assumed Hubble Constant, but also compensate for residual homogeneous Malmquist bias that is not modeled by the simple shape of our distance priors.

The CF3 catalog 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 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% 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%. 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}}$.

The modeling 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 catalog (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.

The remaining non linearities (outside of the groups) at redshift zero are modeled by a unique Gaussian dispersion to the observed velocities. The width of this distribution, namely $\sigma_{\text{NL}}$, is fitted with the other parameters. This is an approximation: non-linear motions are correlated down to $k\gtrsim 0.1$ Mpc^-1 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 $\sigma_{\text{NL}}$ parameter at a value constrained by $\Lambda$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.

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 modeled. In our Bayesian approach, the selection effects are modeled by constructing the probability $\mathcal{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 $\mathcal{P}(d)$ does not depend on the observed distance modulus $\mu$.

We note however that selection effects would need further investigation in future works on peculiar velocities. As suggested by Hinton et al. (2017), a proper model needs the introduction of the generalized likelihood $\mathcal{L}^{\prime}$:

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 $\mathcal{P}(\text{selection}|D,\text{parameters})$ representing the probability of selecting a galaxy measurement given its overall properties and the model’s free parameters. Sampling from the extended likelihood $\mathcal{L}^{\prime}$ is consequently hard and unpractical 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 catalog 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 mentionned in Section 6.1, we tested the effect of changing the priors on the distances and the results are presented in Appendix D.