Cosmic Flow Measurement and Mock Sampling Algorithm of Cosmicflows-4 Tully−Fisher Catalog

Producing realistic mock surveys that model nonlinear gravitational interactions requires N-body simulations. However, such full N-body simulations, which model the gravitational forces and motions down to a fine resolution, can require an extensive HPC infrastructure and a large amount of CPU time. To make a large ensemble of mock surveys quickly and efficiently, we use the COmoving Lagrangian Acceleration (COLA) approach (Tassev et al. 2013) implemented in the l-picola code (Howlett et al. 2015a, 2015b).

We use l-picola to generate 250 dark matter simulations. Though approximate, this method and code have been demonstrated to reproduce the clustering of dark matter particles and halos well on all scales of interest for galaxy and peculiar velocity clustering (Howlett et al. 2015a, 2015b; Blot et al. 2019), at a ≳95% accuracy for both the power spectrum and bispectrum at k = 0.3h Mpc⁻¹. Halos are identified in these simulations using the 3D friends-of-friends algorithm in the VELOCIraptor code (Elahi 2009). This code has also been demonstrated to recover realistic halos, subhalos, and halo merger trees (although we only use the first of these here), in agreement with a variety of other halo finders (Onions et al. 2012; Knebe et al. 2013).

The fiducial cosmology is given by Ω_m = 0.3121, Ω_b =0.0488, σ₈ = 0.815, and h = 0.6751. The redshift of the simulations is at z = 0. The box size is 1800 Mpc h⁻¹. The number of particles in each simulation is 2560³. The minimum halo mass is ∼5 × 10¹¹ M_☉ h⁻¹ (20 particles per halo).

We have 250 simulation boxes in total; each simulation provides us with the masses m_hl of halos and the Cartesian positions [x, y, z] and velocities [v_x , v_y , v_z ] of halos. We divide each simulation box into eight identically sized cubes. In each cube, the observer is placed at a galaxy close to the center of that cube. Therefore, we have generated 2000 mock CF4TF catalogs in total.

The halos in the simulations are parent halos; it is required to generate subhalos and galaxies from these parent halos. The starting point is the halo concentration, defined as

$$\begin{eqnarray}&&{c}_{v}\equiv \displaystyle \frac{{r}_{v}}{{r}_{s}}\end{eqnarray} \tag{ 2 }$$

which characterizes the mass distribution in a halo. The term r_v is the virial radius of a halo, r_s is the break radius between an inner and outer density profile in a halo, and c_v is computed purely as a function of parent halo mass m_v using the fitted relationship from Prada et al. (2012), calibrated on high-resolution N-body simulations. The virial radius is also computed using the parent halo mass and assuming the halo is both spherical and has a density 200× the critical density of the universe. Given c_v and r_v , one can generate the subhalos and galaxies using the Navarro–Frenk–White (NFW) profile (Navarro et al. 1997).

The subhalo generation algorithm is based on Vale & Ostriker (2004), Conroy et al. (2006), Howlett et al. (2015c, 2017) and Qin et al. (2019a) and is detailed below. First, we assume a power-law distribution for the mass ratio between a subhalo and its parent halo. From this power-law distribution, we can compute the expected number of subhalos as (Springel et al. 2008; Giocoli et al. 2008; Elahi et al. 2018)

$$\begin{eqnarray}&&{\lambda }_{\mathrm{Poi}}={\int }_{{f}_{\min }}^{1}{{Af}}_{M}^{-{\alpha }_{h}}{{df}}_{M}\end{eqnarray} \tag{ 3 }$$

where f_M is the mass ratio between the subhalo and its parent. Ratio f_min is the minimum mass ratio we consider, which we set, based on the minimum halo mass of the simulation (see Section 3.1), to 20 times the dark matter particle mass divided by the parent halo mass (i.e., a fixed minimum mass, such that the minimum mass ratio varies based on the parent halo). The free parameters A and α_h will be fitted by matching the density power spectrum of mocks to real data.

To add realism to the sampling process, the actual number of subhalos in each parent halo N_sub is generated assuming a Poisson distribution with mean λ_Poi. The mass ratios of the N_sub subhalos are generated by drawing uniform random numbers in the interval [0, 1], and inverting the quantile distribution for Equation (3). This results in mass ratios in the interval $[{f}_{\min },\,1]$. The subhalo masses are then obtained by multiplying the set of f_M by the mass of the parent halo. As the largest value of f_M generated for the N_sub subhalos is 1, the masses of the subhalos will not be larger than their parent halos.

The positions of subhalos are generated from the NFW profile following the arguments in Robotham & Howlett (2018), and rely on the concentration and virial radius for each parent halo computed as described previously . We generate N_sub random numbers in the interval R ∈ [0, 1], then calculate the position parameter (Robotham & Howlett 2018)

$$\begin{eqnarray}&&p=R\left(\mathrm{ln}(1+{c}_{v})-\displaystyle \frac{{c}_{v}}{1+{c}_{v}}\right)\end{eqnarray} \tag{ 4 }$$

and the radius

$$\begin{eqnarray}&&r=-\displaystyle \frac{{r}_{v}}{{c}_{v}}\left(1+\displaystyle \frac{1}{{W}_{0}(-{e}^{-p-1})}\right)\end{eqnarray} \tag{ 5 }$$

which gives the radial position of the subhalos relative to the center of their parent halo. Here, W₀ is the Lambert Function. ³ N_sub random numbers in the intervals [−π, π] and [0, 2π] are generated as the polar coordinates of the subhalos, respectively. Then we convert r into Cartesian coordinates using the polar coordinates. Finally we add them to the Cartesian positions of their parent halo and subtract the observer's position to obtain the final positions of subhalos.

The velocities of subhalos are calculated from (Navarro et al. 1997; Howlett et al. 2015c)

$$\begin{eqnarray}&&v=\sqrt{\displaystyle \frac{{{Gm}}_{v}s}{{r}_{v}}}\end{eqnarray} \tag{ 6 }$$

where G is the Newtonian Gravitational Constant and s is the ratio between the velocities of subhalos and the circular velocity of their parent halo, given by (Navarro et al. 1997, Howlett et al. 2015c)

$$\begin{eqnarray}&&s=\displaystyle \frac{1}{q}\displaystyle \frac{\mathrm{ln}(1+{c}_{v}q)-\tfrac{{c}_{v}q}{1+{c}_{v}q}}{\mathrm{ln}(1+{c}_{v})-\tfrac{{c}_{v}}{1+{c}_{v}}}.\end{eqnarray} \tag{ 7 }$$

Then we randomly draw the Cartesian components of the velocities v_x , v_y , and v_z for the subhalos from the Gaussian function with the mean equal to the velocity of their parent halo and the standard deviation (std) equal to $v/\sqrt{3}$ (Howlett et al. 2015c).

The i-band absolute magnitude is generated using the luminosity function (see Appendix A for more discussion about the i-band luminosity function). We re-order these magnitudes in descending order, then assign them to the halos (both parent and subhalos) in descending order of mass. We add additional scatter in this one-to-one assignment via a standard deviation parameter ${\sigma }_{\mathrm{log}M}$ to account for the expected scatter between halo/subhalo mass and galaxy luminosity. This is used to draw a "proxy" luminosity for the pairwise-matching of halo mass and luminosity using a Gaussian centered on the actual simulated luminosity. Note that this "proxy" is used only for the matching; the simulated luminosity is the one actually stored and used for later applications of the mocks. Then, we obtain mock galaxies. The normalization factor ϕ_⋆ of the luminosity function will be fitted by matching the density power spectrum of mocks to real data.

We calculate the sky completeness of the CF4TF survey by splitting into five distinct patches (see Appendix B) and comparing to the 2M++ redshift survey (Lavaux et al. 2010). To do this, the CF4TF and 2M++ are gridded using HEALPix (Górski et al. 2005; Zonca et al. 2019). For each HEALPix pixel, we then calculate the ratio of galaxies in CF4TF and 2M++, treating this as the completeness. However there are two caveats to this. First, 2M++ is not quite uniform; it is 1 K-band magnitude deeper in the Sloan Digital Sky Survey (SDSS) and 6dF regions, which we correct by normalizing the sum of the values in the HEALPix pixels in these regions to be the same. Second, we would not expect every galaxy with a redshift in 2M++ to be capable of providing a Tully−Fisher (TF) distance. However, as we are only interested in measuring the completeness of the CF4TF measurements in one region of the sky relative to other regions, the completeness ratio is then normalized such that the mean completeness of the two patches containing ALFALFA data is one. If the completeness of a pixel is greater than one, it is set to be one. We note that this does not give an absolute measure of completeness (i.e., given all of the galaxies that could have TF measurements in a region of the real universe, how many actually have measurements in CF4TF), but it does provide a relative measure, which is all that is needed for the mocks given the following steps. Importantly, this procedure removes any signatures of large-scale structure from the completeness mask, as these structures should be present in both CF4TF and 2M++. For each mock, the simulated galaxies are subsampled to match the smoothed redshift distribution of CF4TF data in each of the five distinct sky patches separately. Given the relative sky completeness above, this downsampling ensures the mocks have comparable numbers of objects to those we actually have TF measurements for, and that the completeness in both the angular and radial coordinates is representative of the real universe.

In a given redshift bin of the real CF4TF data, we use the Gaussian kernel distribution function ⁴ (KDE) to smooth the distribution of (which denotes the measurement error of log-distance ratio η) of that redshift bin. Then we calculate the spline cumulative distribution function (CDF) from the Gaussian KDE for of that redshift bin. The spline CDF is used to generate measurements errors of log-distance ratios for mock galaxies of that redshift bin. We repeat this step for all of the other redshift bins to obtain the measurements errors of log-distance ratio _mock for all mock galaxies. The measured log-distance ratio for each mock galaxy is generated using a Gaussian function centered on the true log-distance ratio η_t with standard deviation _mock. ⁵ As the measured log-distance ratios are centered on the true values, our mocks do not reproduce possible systematics such as Malmquist bias that may affect the real data. To do so would require instead simulating the observed quantities of the Tully−Fisher relationship, and running these through the same pipeline/fitting procedure as the real data. However, such a procedure was already performed when obtaining the CF4TF data in Kourkchi et al. (2020b), so here we assume the Malmquist bias in the data has been corrected for appropriately, and the mocks only need to reproduce the remaining effects of cosmic variance and measurement errors.

Finally, we can measure the redshift space density power spectra of the 2000 mocks using the method ⁶ in Howlett (2019) and Qin et al. (2019a). This is shown in Figure 3, which compares the average of these measured power spectra (blue curve) to the density power spectrum of the real CF4TF data (black filled squares) to find the best values for the parameters A, α_h , ${\sigma }_{\mathrm{log}M}$ and ϕ_⋆, as presented in Table 1. The χ²/d.o.f. = 39/(38 − 4) indicates that the mocks are in excellent agreement with the data.

Zoom In Zoom Out Reset image size

Table 1. The Best-fit Values of the Parameters used in the Mock Sampling Algorithm

A	αh	${\sigma }_{\mathrm{log}M}$	ϕ⋆
2.882	0.102	2.121	0.005

Download table as:  ASCIITypeset image

Figure 4 shows the comparison between the real CF4TF data and the mocks. The top left-hand-side panel shows the i-band absolute magnitude distribution. The black curve is the average of 2000 mocks, and the shaded areas indicate the 1σ, 2σ, and 3σ regions. Here, σ denotes the standard deviation of the magnitude distributions of the 2000 mocks. The dashed blue curve is for the real CF4TF data. The top right-hand-side panel shows the distribution of the log-distance ratio for the data and mock average. The bottom left-hand-side panel shows the distribution of the error of log-distance ratio for the data and mock average. The bottom right-hand-side panel shows the redshift distribution of the data and mock average. In Figure 5, the top panel shows the sky coverage of an example mock in equatorial coordinates. For comparison, the bottom panel shows the sky coverage of the real CF4TF data in equatorial coordinates. The colors of the dots indicate the redshifts based on the color bars.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The mocks are in good agreement with the data. They will be continuously updated based on the updates of the real CF4TF catalog and L-PICOLA simulations.

Using the Taylor series expansion, the line-of-sight peculiar velocity field, v(d_h ), can be expanded to the first order (Jaffe & Kaiser 1995; Parnovsky et al. 2001; Feldman & Watkins 2008; Feldman et al. 2010; Qin et al. 2019b)

$$\begin{eqnarray}&&v({d}_{h})=\sum _{p=1}^{9}{U}_{p}{g}_{p}({d}_{h})\,,\end{eqnarray} \tag{ 8 }$$

where U_p are the nine moment components

$$\begin{eqnarray}&&{U}_{p}=\{{B}_{x},{B}_{y},{B}_{z},{Q}_{{xx}},{Q}_{{yy}},{Q}_{{zz}},{Q}_{{xy}},{Q}_{{xz}},{Q}_{{yz}}\}\,,\end{eqnarray} \tag{ 9 }$$

and where the zeroth-order vector [B_x , B_y , B_z ] is the so-called bulk flow velocity. The first-order symmetric tensor Q_ij , (i, j =x, y, z) is the so-called shear moment. The mode functions are given by

$$\begin{eqnarray}&&\begin{array}{l}{g}_{p}({d}_{h})=\{{\hat{{\boldsymbol{r}}}}_{x},{\hat{{\boldsymbol{r}}}}_{y},{\hat{{\boldsymbol{r}}}}_{z},{d}_{h}{\hat{{\boldsymbol{r}}}}_{x}^{2},{d}_{h}{\hat{{\boldsymbol{r}}}}_{y}^{2},{d}_{h}{\hat{{\boldsymbol{r}}}}_{z}^{2},\\ \quad 2{d}_{h}{\hat{{\boldsymbol{r}}}}_{x}{\hat{{\boldsymbol{r}}}}_{y},2{d}_{h}{\hat{{\boldsymbol{r}}}}_{x}{\hat{{\boldsymbol{r}}}}_{z},2{d}_{h}{\hat{{\boldsymbol{r}}}}_{y}{\hat{{\boldsymbol{r}}}}_{z}\}\end{array}\end{eqnarray} \tag{ 10 }$$

where the true comoving distance d_h of a galaxy is given by

$$\begin{eqnarray}&&{d}_{h}({z}_{h})=\displaystyle \frac{c}{{H}_{0}}{\int }_{0}^{{z}_{h}}\displaystyle \frac{{dz}^{\prime} }{\sqrt{{{\rm{\Omega }}}_{m}{\left(1+{z}_{h}\right)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}},\end{eqnarray} \tag{ 11 }$$

and where z_h is the Hubble recessional redshift of that galaxy. H₀ is the Hubble constant of the present-day universe. The terms Ω_m and Ω_Λ are the matter and dark energy densities of the present-day universe, respectively.

Measuring cosmic flows can provide us with an intuitive understanding of the amplitude and direction of the measured (not the modeled or reconstructed) velocity field. In addition, comparing the measured U_p to the ΛCDM prediction enables us to test whether ΛCDM accurately describe the motions of galaxies in the nearby universe. To avoid the non-Gaussianity of the peculiar velocities, we use the so-called ηMLE (Nusser & Davis 1995, 2011; Qin et al. 2018, 2019b) to estimate the nine moment components U_p from the full-sky galaxy catalog CF4TF.

Assuming the log-distance ratio of the nth galaxy, η_n has Gaussian errors _n , then the Gaussian likelihood function of a set of N galaxies can be written as (Nusser & Davis 1995, 2011; Qin et al. 2018, 2019b)

$$\begin{eqnarray}&&\begin{array}{l}P({U}_{p},{\epsilon }_{\star })=\prod _{n=1}^{N}\displaystyle \frac{1}{\sqrt{2\pi \left({\epsilon }_{n}^{2}+{\epsilon }_{\star ,n}^{2}\right)}}\\ \qquad \qquad \qquad \exp \left(-\displaystyle \frac{1}{2}\displaystyle \frac{{\left({\eta }_{n}-{\tilde{\eta }}_{n}({U}_{p})\right)}^{2}}{{\epsilon }_{n}^{2}+{\epsilon }_{\star ,n}^{2}}\right),\end{array}\end{eqnarray} \tag{ 12 }$$

where the intrinsic scatter of the log-distance ratio _⋆ arises from the nonlinear motions of the galaxies. We set it as a free parameter to fit.

Following the steps of the second paragraph in Section 4.1 of Qin et al. (2019b), the modeled log-distance ratio ${\tilde{\eta }}_{n}({U}_{p})$ of the nth galaxy is converted by first substituting Equation (8) into the usual peculiar velocity estimator (Colless et al. 2001; Hui & Greene 2006; Davis & Scrimgeour 2014; Scrimgeour et al. 2016; Qin et al. 2018, 2019b)

$$\begin{eqnarray}&&v=c\left(\displaystyle \frac{z-{z}_{h}}{1+{z}_{h}}\right)\,,\end{eqnarray} \tag{ 13 }$$

to replace v to calculate d_h (U_p ). Then we can calculate ${\tilde{\eta }}_{n}({U}_{p})$ from d_h (U_p ) and the observed redshift z of that galaxy.

We choose flat priors for the 10 independent parameters. For the three bulk flow components, we use flat priors in the interval B_i ∈ [−1200, 1200] km s⁻¹. For the six shear components, we use flat priors in the interval Q_ij ∈ [−100, 100] h km s⁻¹ Mpc⁻¹. We use flat priors in the interval _⋆ ∈ [−1000, 1000] h km s⁻¹ Mpc⁻¹ for _⋆. Combining these flat priors with the likelihood in Equation (12) to obtain the posterior, then we use the Metropolis−Hastings Markov Chain Monte Carlo ⁷ (MCMC) algorithm to estimate the parameters.

The measurement error of the bulk flow component, ${e}_{{B}_{i}}$ (i = x, y, z) is calculated from the standard deviation of the MCMC samples of the corresponding MCMC chain. They can be converted to the error of the bulk flow amplitude, e_B , using (Scrimgeour et al. 2016; Qin et al. 2018):

$$\begin{eqnarray}&&{e}_{B}^{2}={{JC}}_{{ij}}{J}^{T}\,,(i=1,2,3),\end{eqnarray} \tag{ 14 }$$

where J is the Jacobian ∂B/∂B_i . The covariance of the bulk flow components C_ij is computed using the MCMC samples.

To explore how well the estimator recovers the true moments, we test it using mocks (a comparison between data and cosmological models is presented in Section 5) . The true moment U_p,t of the mocks is defined as the weighted average of the true velocities of the galaxies

$$\begin{eqnarray}&&{U}_{p,t}=\sum _{n=1}^{N}{w}_{p,n}{{\boldsymbol{v}}}_{n,t}\cdot {\hat{{\boldsymbol{r}}}}_{n}\,,\end{eqnarray} \tag{ 15 }$$

where v _n,t is the true velocity of the nth galaxy and is known from the simulation, and ${\hat{{\boldsymbol{r}}}}_{n}$ is the unit vector point to that galaxy. The weights are calculated using the mode function Equation (10), given by

$$\begin{eqnarray}&&{w}_{p,n}=\sum _{q=1}^{9}{A}_{{pq}}^{-1}\displaystyle \frac{{g}_{q,n}}{{\alpha }_{n}^{2}+{\alpha }_{\star }^{2}}\,,\,\,{A}_{{pq}}=\sum _{n=1}^{N}\displaystyle \frac{{g}_{p,n}{g}_{q,n}}{{\alpha }_{n}^{2}+{\alpha }_{\star }^{2}}\,,\end{eqnarray} \tag{ 16 }$$

where α_n is given by Hui & Greene (2006), Johnson et al. (2014), and Adams & Blake (2017):

$$\begin{eqnarray}&&{\alpha }_{n}=\displaystyle \frac{\mathrm{ln}(10){{cz}}_{n}}{1+{z}_{n}}{\epsilon }_{n},\end{eqnarray} \tag{ 17 }$$

which is the measurement error of peculiar velocity. The term α_⋆ is converted from _⋆ using a similar expression.

Figure 6 shows the bulk (top panels) and shear moments (middle and bottom panels) measured from the mocks under equatorial coordinates. The error bars are the measurement errors calculated using the MCMC samples, and hence include the effect of peculiar velocity measurement errors, but not cosmic variance. The impact of cosmic variance can instead be inferred from the spread in true bulk flows between mock catalogs. When comparing our results to cosmological models in Section 5, we include both measurement errors and cosmic variance. A total of 1000 example mocks are shown here. The black dashed lines are the expected one-to-one relation. The colored solid lines are the best fits to the co-responding colored points (moments). The best-fit lines are almost consistent with the one-to-one relation, indicating the ηMLE can recover the true moments. The slopes of the dashed lines are slightly different from the one-to-one relation. The reason for this is most likely due to the fact that we assume that the cosmic flow field is represented simply as bulk and shear moments, without higher order components. We leave the testing of this hypothesis for future work since accurate measurements of higher moments need both larger numbers of galaxies and a more isotropic and homogeneous sky coverage of galaxies.

Zoom In Zoom Out Reset image size

The reduced χ² between the measured moments and true moments is given by

$$\begin{eqnarray}&&{\chi }_{\mathrm{red}}^{2}=\displaystyle \frac{1}{9\times 1000-1}({{\boldsymbol{U}}}_{m}-{{\boldsymbol{U}}}_{t}){{\mathsf{C}}}^{-1}{\left({{\boldsymbol{U}}}_{m}-{{\boldsymbol{U}}}_{t}\right)}^{T},\end{eqnarray} \tag{ 18 }$$

where the measured moments U _m and the true moments U _t contain 9000 elements (1000 mocks times 9 moments), respectively. C is the covariance matrix. Equation (8) gives the nine moments U_p as the components of the cosmic flow field. Since these moments are all being measured from a single observation location, and with incomplete sky coverage, uncertainty of the measurement of one of these moments will be correlated with those of another, and this is true for all moments. C is a 9000 × 9000 matrix with 1000 9 × 9 diagonal blocks and zeros elsewhere. The 1000 diagonal blocks are calculated from the MCMC chains. The ${\chi }_{\mathrm{red}}^{2}=1.806$. This slightly larger value of ${\chi }_{\mathrm{red}}^{2}$ is due to the intrinsic scatter between the measured moments and true moments, which arises from the intrinsic scatter of the peculiar velocities α_⋆ (or _⋆ of log-distance ratio). Parameter α_⋆ (or _⋆) is to encapsulate the nonlinear peculiar motions of galaxies. The peculiar velocity estimator Equation (13) is not good enough to account for the nonlinear peculiar motions. Therefore a more robust peculiar velocity estimator needs to be developed in future work.

Applying ηMLE to the real CF4TF data, we obtain the measured bulk and shear moments in the local universe, as presented in Table 2. The moments are measured in Galactic coordinates. Figure 7 shows the 2D contours and histograms of the MCMC samples of the moments. The shaded areas of the histograms indicate the 1σ errors of the moments. To compare to theory, the characteristic scale of cosmic flow measurement is defined as (Scrimgeour et al. 2016)

$$\begin{eqnarray}&&{d}_{\mathrm{MLE}}=\displaystyle \frac{\sum | {{\boldsymbol{d}}}_{h,n}| {W}_{n}}{\sum {W}_{n}}\,,\end{eqnarray} \tag{ 19 }$$

where the weight factors ${W}_{n}=1/({\alpha }_{n}^{2}+{\alpha }_{\star }^{2})$.

Zoom In Zoom Out Reset image size

Table 2. The Bulk and Shear Moments Measured from the CF4TF Data

	ηMLE	CRMS
Bx (km s−1)	175.9 ± 23.5	±173.4
By (km s−1)	−330.2 ± 21.9	±175.1
Bz (km s−1)	−38.9 ± 17.1	±183.5
Qxx (h km s−1 Mpc−1)	2.53 ± 0.74	±2.62
Qxy (h km s−1 Mpc−1)	−3.58 ± 0.58	±1.70
Qxz (h km s−1 Mpc−1)	1.52 ± 0.43	±1.43
Qyy (h km s−1 Mpc−1)	2.14 ± 0.79	±3.65
Qyz (h km s−1 Mpc−1)	0.47 ± 0.51	±1.77
Qzz (h km s−1 Mpc−1)	−1.16 ± 0.48	±2.72
χ2	10.4
p-value	0.319
Direction of bulk flow	l = 2984 ± 34, b = −59 ± 27
dMLE (h−1 Mpc)	35

Note. To compare to the ΛCDM prediction, we also list the CRMS in the last column. The number of degrees of freedom is 9.

Download table as:  ASCIITypeset image

The ΛCDM predicted moments should have a zero mean and "cosmic root mean square" (CRMS) variation (Feldman et al. 2010). The CRMS is calculated from the following equation (Feldman et al. 2010; Ma et al. 2011; Johnson et al. 2014)

$$\begin{eqnarray}&&{R}_{{pq}}^{v}=\displaystyle \frac{{{\rm{\Omega }}}_{m}^{1.1}{H}_{0}^{2}}{2{\pi }^{2}}\int {{ \mathcal W }}_{{pq}}^{2}(k){ \mathcal P }(k){dk},\end{eqnarray} \tag{ 20 }$$

where the matter density power spectrum ${ \mathcal P }(k)$ is generated using the CAMB package (Lewis et al. 2000). The window function is given by (Feldman et al. 2010; Ma et al. 2011; Johnson et al. 2014):

$$\begin{eqnarray}&&{{ \mathcal W }}_{{pq}}^{2}(k)=\sum _{m,n}^{N}{w}_{p,m}{w}_{q,n}{f}_{{mn}}(k),\end{eqnarray} \tag{ 21 }$$

where f_mn (k) is given by Ma et al. (2011) and Johnson et al. (2014),

$$\begin{eqnarray}&&\begin{array}{l}{f}_{{mn}}(k)=\displaystyle \frac{1}{3}\left[{j}_{0}({{kD}}_{{mn}})-2{j}_{2}({{kD}}_{{mn}})\right]{\hat{{\boldsymbol{r}}}}_{m}\cdot {\hat{{\boldsymbol{r}}}}_{n}\\ \quad \times \displaystyle \frac{1}{{D}_{{mn}}^{2}}{j}_{2}({{kD}}_{{mn}}){r}_{m}{r}_{n}{\sin }^{2}({\alpha }_{{mn}}),\end{array}\end{eqnarray} \tag{ 22 }$$

and where D_mn ≡ ∣ r _m − r _n ∣ and ${\alpha }_{{mn}}={\cos }^{-1}({\hat{{\boldsymbol{r}}}}_{m}\cdot {\hat{{\boldsymbol{r}}}}_{n})$. Variable r _n is the position vector of the nth galaxy, and j₀(x) and j₂(x) are the zeroth- and second-order spherical Bessel functions of the first kind. The weight factors are given by Equation (16).

Figure 8 shows the p = q components of the window function W² _pq of the CF4TF data in Galactic coordinates. In Table 2, we list the CRMS predicted by ΛCDM. As shown in the top panel of Figure 8, the amplitudes of window functions for the three bulk flow components are similar; therefore, the co-responding CRMS of the three bulk flow components are similar in Table 2. As shown in the middle and bottom panels of Figure 8, the window functions for the y direction related components [Q_yy , Q_xy , Q_yz ] have larger amplitudes, indicating a larger cosmic variance in this direction. This is due to the an-isotropic sky coverage caused by the denser region contributed by ALFALFA galaxies. Therefore, in Table 2, the CRMS for the y direction related components, especially Q_yy , is larger than other shear moment components.

Zoom In Zoom Out Reset image size

The ΛCDM predicted moments should have a zero mean. Therefore, the χ² between the nine measured moments and the ΛCDM prediction is given by:

$$\begin{eqnarray}&&{\chi }^{2}={U}_{p}{\left({C}_{{pq}}+{R}_{{pq}}^{v}\right)}^{-1}{U}_{q}^{T}\,,\end{eqnarray} \tag{ 23 }$$

where C_pq is the covariance of the measurement errors and is calculated from the MCMC samples. The χ² = 10.4, and the corresponding p-value is 0.319, indicating the measurements are consistent with the ΛCDM prediction.

To visualize the comparison of different bulk flow measurements from other surveys on a single figure, we need to standardize the window function due to the differing geometries and depths of the surveys. In this work, we use the spherical top-hat window function:

$$\begin{eqnarray}&&{ \mathcal W }(k)=3{\left(\sin {kR}-{kR}\cos {kR})/({kR}\right)}^{3}.\end{eqnarray} \tag{ 24 }$$

The other surveys we will compare to are: W09—Watkins et al. (2009); C11—Colin et al. (2011); D11—Dai et al. (2011); N11—Nusser & Davis (2011); T12—Turnbull et al. (2012); M13—Ma & Scott (2013); H14—Hong et al. (2014); S16—Scrimgeour et al. (2016); Q18—Qin et al. (2018); Q19—Qin et al. (2019b); B20—Boruah et al. (2020); S21—Stahl et al. (2021).

In Figure 9, the pink curve is the ΛCDM prediction calculated from the spherical top-hat window function. The shaded areas indicates the 1σ and 2σ cosmic variance. The yellow stars are the measurements from other surveys. The red dot is our measurement using CF4TF. Following the arguments in Scrimgeour et al. (2016), to be comparable to the top-hat window function prediction, we plot the W09 and T12 measurements at twice their quoted depth since they have Gaussian windows. S16 measures the bulk flow from the 6dFGSv data (Springob et al. 2014), which is a hemispherical top-hat. Following the arguments in S16, we plot their measurement at effective radius ${R}_{\mathrm{eff}}\,={\left({R}^{3}/2\right)}^{1/3}$. Q18 revised the bulk flow measurement of 6dFGSv by calibrating a bias in the Malmquist bias correction of 6dFGSv distance measurements; the revised result is shown as "Q18-6dFGSv." Generally, most of the measurements are consistent with the ΛCDM prediction, although the W09 result does not agree with the ΛCDM prediction.

Zoom In Zoom Out Reset image size

Comparing to other measurements (yellow stars), our measurement (red dot) has a smaller error. B20 has the smallest error bar; they measure the bulk flow by comparing the reconstructed velocity field of 2M++ (Carrick et al. 2015) to the measured velocity field of 465 supernovae. However, as presented in Section 5.1 of Carrick et al. (2015), the reconstructed field has many potential systematic effects that are hard to test. Q19 measures the bulk flow using the Cosmicflows-3 catalog (Tully et al. 2016), which contains 391 Type Ia supernovae and, therefore, has a smaller error bar too. However, the major part of the published Cosmicflows-3 catalog is the 6dFGSv, which does not have a perfect Malmquist bias correction (Qin et al. 2018); therefore, this measurement has potential systematic errors.

Figure 10 shows directions of bulk flows measured from different surveys in Galactic coordinates. The S16 measurements is significantly biased from other surveys. The revised result is shown as "Q18-6dFGSv."

Zoom In Zoom Out Reset image size