Uneven flows: On cosmic bulk flows, local observers, and gravity

# Uneven flows: On cosmic bulk flows, local observers, and gravity

## Abstract

Using N-body simulations we study the impact of various systematic effects on the low-order moments of the cosmic velocity field: the bulk flow (BF) and the Cosmic Mach Number (CMN). We consider two types of systematics: those related to survey properties and those induced by observer’s location in the Universe. In the former category we model sparse sampling, velocity errors, and survey incompleteness (radial and geometrical). In the latter, we consider Local Group (LG) analogue observers, placed in a specific location within the Cosmic Web, satisfying various observational criteria. We differentiate such LG observers from Copernican ones, who are at random locations. We report strong systematic effects on the measured BF and CMN induced by sparse sampling, velocity errors and radial incompleteness. For BF most of these effects exceed 10% for scales . For CMN some of these systematics can be catastrophically large (i.e. ) also on bigger scales. Moreover, we find that the position of the observer in the Cosmic Web significantly affects the locally measured BF (CMN), with effects as large as ( at for a LG-like observer as compared to a random one. This effect is comparable to the sample variance at the same scales. Such location-dependent effects have not been considered previously in BF and CMN studies and here we report their magnitude and scale for the first time. To highlight the importance of these systematics, we additionally study a model of modified gravity with enhanced growth rate (compared to general relativity). We found that the systematic effects can mimic the modified gravity signal. The worst-case scenario is realized for a case of a LG-like observer, when the effects induced by local structures are degenerate with the enhanced growth rate fostered by modified gravity. Our results indicate that dedicated constrained simulations and realistic mock galaxy catalogs will be absolutely necessary to fully benefit from the statistical power of the forthcoming peculiar velocity data from surveys such as TAIPAN, WALLABY, Cosmic Flows-4 and SKA.

## I Introduction

The standard model of cosmology – Lambda Cold Dark Matter (LCDM) - is extremely successful in explaining a plethora of observations. These include the features of the Cosmic Microwave Background (e.g. Hinshaw et al., 2013; Planck Collaboration et al., 2016), the primordial nucleosynthesis and light element abundances (Yang et al., 1984; Walker et al., 1991), the growth of primordial density perturbations into the present-day large-scale structure (LSS) (Percival et al., 2001; Tegmark et al., 2004; Beutler et al., 2017), as well as the late-time accelerated expansion of the Universe (Riess et al., 1998; Perlmutter et al., 1999; Percival et al., 2010; Weinberg et al., 2013). However, since LCDM is mostly phenomenological in its nature, this spectacular success comes at a price of accepting that the main contributors to the cosmic energy budget are dark matter (DM) and dark energy (DE), which have not been directly detected in any experiments so far. Therefore, it is desirable to look for other probes of the cosmological model, especially those which do not share at least some of the systematics of the aforementioned measurements.

In this context, the peculiar motions of galaxies – i.e. deviations from the uniform Hubble flow – are considered as particularly valuable Strauss and Willick (1995); Koda et al. (2014); Howlett et al. (2017). Induced by gravity only, they are not affected by such systematics as galaxy bias, which plagues for instance the measurements of galaxy clustering. Peculiar velocities can be therefore used, at least in principle, to obtain constraints on various cosmological parameters such as the mean matter density or the growth of structure Nusser and Davis (2011); Hudson and Turnbull (2012), independently of other methods.

Arguably the most popular statistics of the velocity field is the bulk flow (BF), i.e. the net peculiar motion of galaxies contained in a given volume. BF probes large-scale fluctuations of matter distribution, and should generally diminish with increased volume. Over the decades, BF measurements have often been subject to various controversies. An example from early studies is by Lauer and Postman (1994), who measured a net motion of Abell clusters amounting to within a radius of 15,000 , which was however not confirmed by subsequent analyses (e.g. Giovanelli et al., 1996; Dale et al., 1999) (but see Hudson et al. (1999)). More recently, Watkins et al. (2009) claimed significant BF () on scales of from a combined sample of galaxies and clusters, which also is not supported by several other studies (e.g. Nusser and Davis, 2011; Hoffman et al., 2015) (see however Watkins and Feldman (2015)). Even more controversial are the claims of the very large scale () ‘dark flow’ by Kashlinsky et al. (2008), which again is not corroborated by related analyses Osborne et al. (2011); Planck Collaboration et al. (2014a). Thanks to the ever growing amount of observational data, there is continued interest in measuring the BF and, if these discrepancies could be resolved, using it as a cosmological probe; for some more recent results see Colin et al. (2011); Nusser et al. (2011); Weyant et al. (2011); Branchini et al. (2012); Mody and Hajian (2012); Turnbull et al. (2012); Feindt et al. (2013); Lavaux et al. (2013); Ma and Scott (2013); Rathaus et al. (2013); Feix et al. (2014); Hong et al. (2014); Ma and Pan (2014); Appleby et al. (2015); Huterer et al. (2015); Scrimgeour et al. (2016); Springob et al. (2016).

Part of the BF ‘controversy’ (or more precisely, inconsistency between some measurements) is due to the fact that many of the BF assessments are not directly comparable due to different estimators used, with specific sensitivity to various scales and systematics Nusser (2014); Andersen et al. (2016); Nusser (2016). The quality and volume of the velocity data is another important issue here. We note that some of the developed estimators do not use peculiar velocities at all to estimate the BF (e.g. Kashlinsky et al., 2008; Nusser et al., 2011; Branchini et al., 2012; Mody and Hajian, 2012; Lavaux et al., 2013; Feix et al., 2014), they are thus not sensitive to the related biases, although this of course does not make them immune to other, often major, systematic effects.

The BF continues to be regarded as a promising probe of cosmology especially taking into account that larger, denser, and more accurate samples of peculiar velocities are expected to appear in the coming years from such surveys as Taipan da Cunha et al. (2017), WALLABY Duffy et al. (2012), or CosmicFlows-4 Tully et al. (2016). However, agreement is gradually building up that in order to take full advantage of these future datasets for BF and other velocity-based measurements, the control of systematic effects and biases is crucial for proper data interpretation. Recent developments of e.g.Andersen et al. (2016); Hellwing et al. (2017a); Li et al. (2012a) highlight the importance of selection and observer-driven effects for peculiar velocity studies. Ref. Andersen et al. (2016) considered the impact of purely geometrical selection effects on the inferred bulk flows, including the partial sky coverage. In addition, Ref. Li et al. (2012a) investigated mainly the effect of different radial selection functions and the corresponding galaxy/halo weighting. Both works report the importance of these two systematic effects that can bias the data, but the effect they studied referred to a hypothetical Copernican observer. The results of Hellwing et al. (2017a) however underline that for relatively shallow and sparse velocity data, the specific location of the observer within the cosmic web affects in a non-trivial way the cosmic variance of the velocity observables. Inspired by these previous results, in this work we will readdress this issue by looking closely at the impact of the observer location (i.e. importance of the local cosmic structures) on the inferred BF and related statistics. We will show that the BF itself is very sensitive to such effects, which must be therefore properly accounted for when measuring it from the current and forthcoming datasets.

A statistics related to the BF, which uses additionally the third moment of the peculiar velocity field, is the cosmic Mach number (CMN) defined as the ratio of the BF to the peculiar velocity dispersion in the same volume Ostriker and Suto (1990). In the original proposal, CMN was regarded as a ‘critical test for current theories’, and more recently quoted as ‘a sensitive probe for the growth of structure’ Ma et al. (2012). For other theoretical and observational studies of CMN and its importance for cosmology, see Suto et al. (1992); Strauss et al. (1993); Nagamine et al. (2001); Atrio-Barandela et al. (2004); Agarwal and Feldman (2013). In this paper we will examine the sensitivity of the CMN to the same systematics as those studies for the BF. Similar conclusions regarding the importance of such effects for CMN as in the case of BF will apply.

The cosmic velocity field reflects a continuous action of gravity integrated over the history of large-scale-structure growth. Thus it offers, in principle, a very sensitive probe of the very nature of gravity itself. Here, even small possible deviations from a general relativity (GR)-like force law provide minute galaxy acceleration changes that are amplified when integrated over time. This has been shown by other authors for a range of velocity field statistics and viable modified gravity (MG) models (see e.g.Li et al. (2013a); Hellwing et al. (2014); Ivarsen et al. (2016)). Thus, if one is able to control various systematic effects, and in the case of known (assumed) cosmological parameters like and (taken for example from CMB observations), then the galaxy velocity field (and its low-order moments) provides a potentially powerful way of constraining non-GR models. Such constraints would foster an independent, thus complimentary, way of testing GR and measuring the local value of the growth rate, Koda et al. (2014); Howlett et al. (2017). In order to be able to use the velocity data for testing gravity one needs to recognize and control all important systematic effects. Consequently, in this paper we also consider a modified gravity model (deviating by in the growth rate from GR) and compare its signal with the magnitude of various systematics in the GR case.

As briefly indicated above, various systematic and statistical effects that disturb the velocity data were a subject of careful study in the past. However, except for the early work of Ref. Tormen et al. (1993), analyses of the impact of a specific location of the observer within the large-scale structure were not conducted. Ref. Tormen et al. (1993) studied only the 2-point velocity statistics and they did not require the presence of any nearby Cosmic Web structures such as the Virgo cluster. Here, we will conduct a joint study of various systematic effects, starting from sparse sampling and radial selection, up to the impact of the observer location. We will identify scales and magnitudes of various effects and compare them against expected statistical fluctuations in a systematic fashion. In this way we will obtain insight into scales, magnitudes and the interdependence of various systematical effects troubling BF and CMN measurements. This will constitute another important step for peculiar velocity studies towards the precision cosmology era.

The paper is organized as follows: in §II we describe in details computer simulations used in this study; in §III some theoretical preliminaries and relevant considerations are given; §IV contains description of mock catalogs and various observational effects that we model; in §V we discuss the effects induced by systematics independent from a specific observer’s position, while in §VI the focus is given to signals measured by Local Group-analogue observers; §VII compares signal from a modified gravity model with known GR systematic effects. Finally §VIII summarizes our findings, this is followed by §IX where discussion and conclusions are given. Some additional tests and discussion about the influence of the simulation box size are given in the Appendix.

## Ii Simulations

To study cosmic flows we employ a set of large N-body simulations conducted with the use of the ECOSMOG code Li et al. (2012b). Time evolution of cosmic structures is here followed with respect to a background spatially flat Universe described by cosmological parameters consistent with the 2013 results from the Planck mission Planck Collaboration et al. (2014b). We imposed the following values: , , , . The growth of density fluctuations is modeled by assuming that all non-relativistic matter is collisionless, i.e. we treat the baryonic component as DM. Ignoring baryonic physics will not introduce any significant biases as long as we are not interested in internal properties of individual halos but only in their spatial distribution and peculiar velocities (e.g.Hellwing et al. (2016)).Thus in our simulation we place DM particles in a cubic box of comoving size of . This particular set-up fixes the mass resolution at . ECOSMOG is an extension of the RAMSES code Teyssier (2002) and uses adaptive mesh refinement (AMR) and dynamical grid relaxation methods to compute the gravitational potential and forces. Thus our simulations are not characterized by a single fixed force resolution, but to gauge our dynamical spatial resolution we can use the cell size of the most refined AMR grid. In all our runs such a grid had a rank of resulting in the finest force resolution of .

In this paper we aim to study various systematic effects that affect the lowest moments of the cosmic velocity field. For that reason we will be mostly concerned with the fiducial cosmological model. However, in Sec.VII we will compare the magnitude of various systematics with the predicted amplitude of a non-GR signature expected in the case of a modified gravity model. As a representative guinea-pig we chose the so-called normal branch of Dvali-Gabadadze-Poratti (henceforth nDGP) model Dvali et al. (2000); Schmidt (2009a). For a more detailed description of that model and its implementation in simulations, see the relevant Section.

Real astronomical observations measure the radial component of the peculiar velocity of a galaxy rather than of its host halo. Our simulations do not attempt to model assembly of galaxies in any way, but we can safely use the bulk velocities of DM halos found in our simulations as faithful proxies for real galaxy peculiar velocities. This is the case since the studies of other authors (e.g. Hellwing et al., 2016; Ye et al., 2017) have shown that for central galaxies residing in massive halos their relative velocities with respect to their host halos are very small () compared to the bulk-flow magnitude we will study here. In addition we do not expect that any non-zero galaxy velocity in relation to its host halo would be correlated to the large-scale matter distribution which induces bulk flows. Thus, due to global isotropy these velocities should average-out to zero for scales much larger than a given halo radius.

To identify halos and subhalos in the simulations we employ ROCKSTAR Behroozi et al. (2013), a phase-space friends-of-friends halo finder. To define a halo edge we use the virial radius , defined as the radius within which the enclosed density is , where is the critical cosmic density. For further analysis we keep all gravitationally bound halos that contain at least 20 DM particles each. This sets our minimal halo mass to . Based on the initial halo catalogs, we build our test halo populations by distinguishing the central halos from satellites (subhalos). For further analysis we keep only the centrals, which we will treat as rough mocks for population of central galaxies. To obtain additional halo samples with lower number densities we perform random subsampling. Our main catalog includes all central halos resolved in our simulation at and has a number density of Mpc. To obtain sparser samples we consequently dilute this main sample by randomly (and spatially uniformly) keeping only every -th halo. Thus we also obtain the following samples: Mpc, Mpcand Mpc.

Finally, since peculiar velocity catalogs are rather shallow, rarely reaching at present deeper than (see e.g. Springob et al. (2007, 2014, 2016); Tully et al. (2016)) we constrain all our analysis only to the snapshot of our simulations and to scales up to . This being said, it is also imperative to comment on the convergence of the simulation results at large scales. The velocity field is much more sensitive to contributions from perturbation modes much larger than a given scale one considers. In other words, we can expect that the finite-volume effects will be more pronounced here than in the case of the density field. To check what scales we can trust, we have run additional tests involving three more simplified simulations with a varying box size. The details and analysis of these are given in the Appendix. The results of these tests indicate that on scales the amplitude of our BF is systematically biased down by or more. However, the size of various systematic effects expressed as a relative BF magnitude difference appears to be only weakly affected by the box size up to . This supports our choice of the maximal scale we consider in this paper.

## Iii Theoretical preliminaries

Throughout our work we assume the homogeneous and isotropic cosmological model, in which the background obeys Friedman-Lemâitre equations with a scale factor . All the quantities will be expressed in comoving coordinates i.e. . For background density and density contrast , the Poisson equation linking the peculiar gravitational potential with density perturbations is

 ∇2ϕ(→x,t)=4πGρb(t)a2δ(→x,t). (1)

By integration, we obtain the expression for peculiar accelerations Peebles (1980)

 →g(→x)=−∇ϕa=Gaρ0∫δ(→x′)(→x′−→x)|→x′−→x|3d→x′. (2)

Peculiar velocities , defined as deviations from the Hubble flow, are coupled to the density field via the continuity equation:

 ∂δ∂t+1a∇⋅[(1+δ)→v]=0. (3)

### iii.1 Linear theory predictions

We can model the cosmic velocity field by performing a decomposition of the full three-dimensional (3D) field into a sum of longitudinal (non-rotational) and transverse (rotational) components:

 →v=→vL+→vT,where: (4) ∇×→vL=0and∇⋅→vT=0. (5)

In the linear regime, the velocity field is curl-free, thus and the field is purely potential. Henceforth it can be expressed as a gradient of a scalar function (called the velocity potential)

 →v=−∇Ψv/a. (6)

Now considering the continuity equation (3) it can be shown that the velocity potential obeys

 ∇2Ψv=Hfa2δ, (7)

where we have used the definition of the dimensionless growth rate . The growth rate only very weakly depends on Lahav et al. (1991) and for a flat LCDM universe Linder (2005). However, in general for some alternative cosmologies (like coupled DE or modified gravity) it can take a different value and also be a scale-dependent function.

Finally, in the linear regime we have and (where is the peculiar gravitational acceleration), and in particular at one has

 →v=H0f4πGρ0→g=2f3H0Ωm→g. (8)

In the linear regime we can also express the relation between the power spectrum of density fluctuations, , and the dimensionless expansion scalar, , which is the scaled velocity divergence (also called the expansion scalar)

 θ≡−∇→vaH0,and→vk=aH0i→kk2θk,soPθθ(k)=⟨θkθ∗k⟩. (9)

In the linear regime for a potential flow it follows from the continuity equation (3) that

 P(k)=f−2Pθθ(k). (10)

The above relation is often used in the literature to approximate velocity power spectrum by linear velocity divergence, thus neglecting dispersion and vorticity (see e.g.Strauss and Willick (1995)). Such approximation however holds only on sufficiently large scales; those scales are generally larger (i.e. ) than in relevant analyses of the density field (see e.g.Chodorowski and Ciecielag (2002); Ciecielg et al. (2003); Pueblas and Scoccimarro (2009)).

### iii.2 Bulk flow, velocity dispersion, and cosmic Mach number

The bulk flow (BF) is the dipole (second) moment of the peculiar velocity field, , in a given region of space (volume). Non-zero BF reflects a net streaming motion towards a particular direction in space. Thus in the continuous limit of the field , for a spherical region with a radius , it will be

 →B(R)=34πR3∫R0→v(→x)d3x. (11)

Throughout this paper we will interchangeably use BF and to denote the bulk flow amplitude.

When the velocity field is sampled by a set of discrete tracers (e.g. galaxies) then the above integral becomes a finite sum. If each individual galaxy is assigned a weight , then the 3D bulk flow vector will be

 →B(R)=∑Ni=1wi→vi∑Ni=1wi, (12)

where is the peculiar velocity of the -th galaxy. The corresponding dispersion (VD) of the peculiar velocities with respect to the averaged bulk flow is

 →D(R)=∑Ni=1[wi→vi−→B(R)]2∑Ni=1wi−1, (13)

where the sum of the weights needs to be , so the denominator does not take a zero value. If the density fluctuations are a random Gaussian field, then in the linear theory (i.e. on sufficiently large scales) the corresponding velocity field will also be a random variable (for each vector component separately) with a zero mean and the variance given by the velocity power spectrum , where and we already assumed global isotropy. Thus the predicted root mean square value of the bulk flow amplitude is

 B2(R)=12π2∫dkk2Pvv(k)|^W(kR)|2. (14)

Here is the Fourier image of the window function. Usually one takes to be spherical top-hat, which implies , but some authors consider also the so-called all-sky Gaussian selection function with .

Now, if there is no velocity bias and the velocity field is curl-free, then , and equation (14) becomes

 B2(R)=H202π2∫dkPθθ(k)|^W(kR)|2. (15)

In the regime where the velocity vorticity is negligible and the Eqn. (10) holds, one finally obtains

 B2(R)=H20f22π2∫dkP(k)|^W(kR)|2. (16)

The above equation is commonly used as the linear theory prediction for the bulk-flow amplitude in a universe described by a particular choice of and . Consequently, the corresponding variance of the residual velocity field (after the BF was subtracted) for that case takes the form

 D2(R)=H20f22π2∫dkP(k)(1−|^W(kR)|2). (17)

Now to obtain predictions for the bulk flow amplitude and some significance intervals, a model distribution function for peculiar velocities is needed. This is obtained by noticing that for sufficiently large smoothing scales, the distribution for a single velocity component approaches a Gaussian, thus the distribution for the bulk flow magnitude becomes Maxwellian (see Bahcall et al. (1994); Coles and Lucchin (2002)). Hence for a velocity field with rms velocity of , this is given by

 p(v)dv=√2π(3B2)3/2v2exp(−3v22B2)dv. (18)

Considering gives in the limit the most likely value (MLV) and the expected value (EV) . MLV and EV are widely used as common linear theory (LT) predictions for the BF amplitude, and in the reminder of this manuscript we shall adopt the same strategy whenever we will be invoking LT formulas. We caution however, that in this context it is important to bear in mind that such predictions only hold if the distribution of the components of is Gaussian. The validity of this assumption depends on scales which one considers. Although in general it was established that for most scales dealt with in modern velocity analysis (i.e. ) this assumption generally holds Andersen et al. (2016), results shown by other authors imply that caution should be taken (see also Chodorowski and Ciecielag, 2002; Ciecielg et al., 2003).

A separate note should be made here about the limits of the integrals used to calculate and from Eqns. (14)-(17). To obtain predictions for the physical Universe one should consider the obvious limits from to . However, when we want to compare LT predictions with numerical simulations that used some finite computation box, we should account for the fact that the modeled velocity field will miss the contribution from the modes larger than the box length . Also due to discretization of both mass and volume there is some characteristic minimal scale that is still resolved by the simulation, usually taken to be the force resolution . In such a case, the corresponding integration limits are then . Whenever we will be comparing LT predictions with the simulation results we will employ the above integration limits.

Some authors Ostriker and Suto (1990); Ma et al. (2012); Agarwal and Feldman (2013) advocated also another type of statistics, namely the cosmic Mach number (CMN, or ), that we can define now as

 M(R)≡B(R)D(R), (19)

which in the linear regime should be only a function of the shape of the matter power spectrum (or the effective slope of around ) Suto et al. (1992); Strauss et al. (1993); Watkins and Feldman (2007).

The above considerations suggest that the linear theory prediction for the bulk flow and associated statistics should strongly depend on two parameters of the underlying cosmological model, namely the growth rate and the amplitude of , which can be evaluated by the parameter. These dependencies have motivated many authors to advocate the use of the low-order velocity field statistics as cosmological probes Strauss and Willick (1995); Nusser et al. (2011); Branchini et al. (2012); Ma et al. (2012); Hong et al. (2014); Feix et al. (2014); Koda et al. (2014); Hoffman et al. (2015); Scrimgeour et al. (2016); Howlett et al. (2017).

To gauge the magnitude of variations and their co-dependence on and we have considered a number of power spectra variants. The fiducial case is i) the Planck13 cosmology (Planck Collaboration et al. (2014b); the same as used in our simulations) and we also examined four cases: ii) high (); iii) low (); iv) high (=0.9) and v) low (=0.75). Here, for each case i)–v) we kept fixed all the remaining parameters, imposing and , and varied only the value of a given matter density or power spectrum normalization. By changing we probe different values of the growth rate (by around the fiducial case) and by varying we sample different power spectrum amplitudes. For all the cases we have used the CAMB software package Lewis et al. (2000) to obtain high-accuracy linear matter power spectra and then applied the halofit model Smith et al. (2003) to evolve the spectra to the non-linear regime. In addition, we also considered one more case, where we used the fully non-linear estimated from our simulation. The non-linear velocity divergence power spectrum was used only for , where it deviates by more than from the non-linear ; for smaller it was substituted by the CAMB-provided , rescaled by . We checked the effect of the non-linear divergence spectra, since the scales at which the velocity field is curl-free and the scales at which are not necessary the same Ciecielg et al. (2003); Pueblas and Scoccimarro (2009).

In Fig. 1 we compare all the examined power spectra with our fiducial Planck13 case i). The velocity divergence power spectrum was scaled by the corresponding factor. We can observe that for the cases where is varied, the corresponding changes in are limited to large scales, . Small deviations seen above reflect the different length of the matter-dominated epochs in low and high- universes and so different degree of non-linearity in the density field. However, this appears at scales too small to be relevant for the large-scale velocity field. As expected the high-(low-) case is characterized by a smaller (larger) amplitude of the power spectrum than the fiducial case at these scales. For both cases the changes in the large-scale amplitudes are quite dramatic. Variations in alone affect the spectrum on all scales, but the overall effect is much smaller (typically within ). Here we can also note that the small-scale variance of is strongly suppressed compared to the matter . This is expected, once one considers that in the non-linear regime, while the collapsed objects increase the density field variance, the corresponding velocity field around and inside those objects attains a high-degree of vorticity and dispersion due to shell crossing and virialization Pichon and Bernardeau (1999); Pueblas and Scoccimarro (2009); Libeskind et al. (2013, 2014).

Figure 2 illustrates the changes, imposed due to variations in shape and amplitude, in the corresponding estimated  (symbols) and  (lines). The previously seen dramatic differences in amplitudes are translated to rather mild impact on the resulting linear-theory bulk flows. Here, for the most cases, the changes are within , thus of the same magnitude as our variations in both and . We can also notice the known degeneracy, where the effect of increasing one parameter can be to a large extent compensated by the decrease of the other. The effect of using the non-linear to predict   is minimal for . In contrast, the use of the non-linear velocity divergence power spectrum results in a much more dramatic effects onto the   estimator. This suggests that modeling of non-linearities in the density and velocity field is not that important for   predictors, but might be crucial for the prediction of the expected Mach number. The latter fact was already emphasized to some extent by Agarwal and Feldman (2013), who noticed that in order to obtain more accurate predictions for the Mach number some non-linear corrections for   have to be applied. This reflects the fact that the velocity dispersion is intrinsically a local quantity, and non-linear effects such as virialization and shell-crossing have a significant effect (see e.g.Hahn et al. (2013, 2015); Rampf and Frisch (2017)).

## Iv The velocity mocks and non-linear observables

To move beyond the linear theory we employ the set of -body simulations described in Sec. II. To study various systematics, non-linear effects and biases, and to get a closer connection with real astronomical observations, we construct a set of mock catalogs and observables from our simulations. As an input for all our analysis we consider halo and subhalo catalogs saved at .

Generally, when considering various observational errors and systematics (like survey geometry, selection function, radial distribution, etc.) one can apply their modeling to the simulation data and then analyze the mock catalog by computing various statistics from it. We adopt this routine approach by calculating various data points weights, which characterize different modeled effects in separate mock catalogs.

We consider the following “observational effects” on the data:

• observer location – all the relevant quantities, such as distances and angles, depend on a specific observer location, whether it would be a random or pre-selected observer; computations are done in the CMB rest frame;

• radial selection – we model the following radial selections: 1) full completeness (i.e. no radial selection nor distance limit); 2) CosmicFlows-3-like Tully et al. (2016) selection functions (see below);

• geometry/Zone of Avoidance – since all our catalogs are observer-dependent, it is natural to also include the effect of the so-called Zone of Avoidance (ZoA) caused by obscuration of the far-away objects by the Galactic disk. This is done by removing galaxies from the appropriate part of the volume. See more details below. In our analysis we do not model the importance of particular structures hidden behind the ZoA, such as the Norma Cluster Kraan-Korteweg et al. (1996) or recently discovered Vela Supercluster Kraan-Korteweg et al. (2017), as this would require detailed constrained simulations. We postpone such studies for future work;

• radial velocity error – to model peculiar velocity errors associated with the uncertainties of galaxy scaling relations that are used to infer galaxy velocities from redshifts (see more below).

In our analysis we are concerned with lower-order velocity statistics that are estimated from specific observer-dependent mock catalogs. Therefore, all our results (unless clearly emphasized otherwise) are computed as ensemble averages over all mock observers in a given sample. Refs.Li et al. (2012a) and Andersen et al. (2016) have shown that the distribution of bulk flows amplitudes inferred from simulations deviates from a Gaussian. We have checked that this is the case for all our samples, both for the bulk flows as well as for the velocity dispersions. For that reason, a simple averaged mean and associated variance might not be a faithful characterization of the underlying ensemble. Thus we decided to use medians and associated and percentiles to characterize all our results.

Observer location. All the observables we discuss later in the paper were estimated for a fixed given number of observer locations. By construction, all our observers must sit in a DM-halo. We consider two types of observer locations: random and pre-selected. Random observers are chosen randomly from all halo positions in a given catalog, while the pre-selected are contained in a closed list of locations predefined by some user provided criteria. In this paper we consider various criteria of a hypothetical Local Group (LG)-like observer. See more in Sec. VI.

Radial selection. Generally, to obtain the desired radial selection, we would have to select multiple times from an input data set according to probability that is inversely proportional to the defined shape of the selection function, keeping finally the data product with mock radial selection that is closest to the imposed one. Such a procedure for large samples as ours is however very unpractical. We decided to use a simple data weighting scheme instead, where each galaxy is given a weight exactly as defined by the input selection function. For a large number of galaxies, the results of both procedures give comparable results. Therefore, since we do not compare our results to any particular galaxy survey, but rather aim at providing general observational data modeling, we are satisfied with the much faster data weighting method. We opt to use weighting scheme that follows the radial selection of the Cosmic-Flows 3 catalog for the sake of simplicity and generality. CF3 is currently the largest peculiar velocity catalog, thus by studying CF3-like radial selection our model will be close to the best-case data scenario. In the case 1) listed above, all the halos have equal unit weights, as in Eq. (12). When modeling the CF3-like radial selection 2), we impose Hellwing et al. (2017a):

 wh={1,if r≤rw(r/rw)−m,otherwise. (20)

Here, is the characteristic radial depth of the catalog (in ). For our CF3-like catalogs we consider and two values for the exponent ;

Geometry / Zone of Avoidance. Most extragalactic observations, including those of peculiar velocities, do not have access to low Galactic latitudes due to the obscuration by dust, gas, and stars in the Milky Way – the ‘Zone of Avoidance’. To model it, we consider a small opening angle Kraan-Korteweg and Lahav (2000) chosen with respect to a fixed observer-dependent local plane. Galaxies falling inside are removed.

Radial velocity error. Galaxy peculiar velocity surveys rely on redshift-independent distance-indicators relations (DIs) to extract the cosmological and peculiar components from a galaxy redshift. The most commonly used methods are based on galaxy scaling relations, such as Tully-Fisher Tully and Fisher (1977) or Fundamental Plane Djorgovski and Davis (1987). Such methods are unavoidably plagued with significant relative errors on estimated velocities stemming from intrinsic scatter in used relations and various systematic (usually non-linear) biases. The peculiar velocity errors are a source of a serious worry and their magnitude sets a fundamental limit on cosmic velocity data usability. A constant relative error in distance determination translates here to a velocity uncertainty that grows linearly with galaxy redshift. We attempt to model this by a simple relation of the from:

 σv=ΔvH0Dz. (21)

Here is the Hubble parameter, is the galaxy co-moving distance and models the typical scatter of the logarithmic distance ratio error. The ratio is used to estimate the peculiar velocity. Here is the co-moving distance to a galaxy inferred via DIs (see more in e.g. Courtois et al. (2013); Scrimgeour et al. (2016); Tully et al. (2016)) and the spectroscopic galaxy redshift . We choose , which is a conservative value when compared with smaller scatter typically found in modern velocity data Tully et al. (2016). We assume that the above velocity error is Gaussian with zero mean and dispersion . In reality such an assumption is often broken for various velocity estimators, but we adopt it for simplicity, as non-Gaussian contributions to velocity errors depend strongly on particular galaxy catalog specifics.

Once parameters for mock galaxy catalogs are chosen, we compute the bulk flow and the dispersion of the residual velocity field by assigning specific halo/galaxy weights and using Eqn.(12)-(13). We sum separately over the three Cartesian velocity vector components in concentric spheres of radius around a fixed observer location. This procedure yields us specific weighted bulk flow components, i.e. , and . The bulk flow amplitude is then

 ~B(R)=(3∑iBi(R)2)1/2. (22)

Here the sum runs over three Cartesian components of a 3D velocity vector field and the procedure for the residual velocity dispersion is analogous.

In reality, the above procedure cannot be applied to real data, since except for a very few cases, we do not have full 3D peculiar velocity information. What is directly accessible is only the line-of-sight (l.o.s.) velocity component. Thus for observational data one usually adopts an estimator of the BF that is based on the radial velocity component. For example, in the most popular Maximum Likelihood (ML) method, the BF components are obtained via

 ~Bi=N∑nwi,nVn, (23)

where again indicates one of the three Cartesian indexes, is an th l.o.s. velocity measurement. Here, is an associated weight of a given velocity measurement, which usually is taken to be

 wi,n=3∑jA−1ij^rn,jσ2n+σ2∗, (24)

where is a unit vector from the observer to a given galaxy , is the uncertainty of a given velocity measurement, and describes 1D velocity dispersion due local non-linear virial motions. The matrix describes geometric moments of the whole sample of tracers, and is given by

 Aij=N∑n^rn,i^rn,jσ2n+σ2∗. (25)

The above estimator is based of inverse variance weighting method of Ref. Kaiser (1988).

We do not choose to implement the above estimator for various reasons. First, it is uniquely defined for a given astronomical data set, with its specific radial and geometrical selections and errors of velocity estimates. To keep our discussion as general as possible we opt to use a much simpler estimator of Eqn. (22) instead. This is justified since all our mock catalogs are isotropic and spatially uniform. For such a case the geometric matrix is uniform and approximates a product of a constant factor and a unit matrix. In addition, since we only use central halos, contributions from any non-linear virial motions are strongly suppressed. The last statement does not hold for non-relaxed systems, but those constitute a marginal fraction of our halo catalog. Thus we also opt to drop the non-linear velocity dispersion contribution, , from our modeling.

Finally, taking into account above considerations and for the sake of simplicity, we choose to use a maximally simplified ML estimator, which only includes individual velocity errors in the data weights drawn from a Gaussian distribution independently for each velocity component according to the prescription of Eqn. (21).

## V Observer-independent systematics

Here we will present the results of our analysis of the BF and CMN inferred from mock catalogs where the observer location was kept random and unspecified, i.e. it corresponds statistically (after averaging) to a Copernican observer (see more in Hellwing et al. (2017a)). By adopting this approach we will be able to study various systematic effects that are, in principle, independent from the location. By doing this we can assess how much the various systematics can affect the measurements in an idealized survey.

### v.1 Bulk flow

We begin by investigating how the sampling rate or number density of tracers used for the measurement affects the resulting BF. In Fig. 3 we show the median bulk flow measured for the full sample which is characterized by Mpc, and for three catalogs with lower number density of tracers, namely , , and Mpc, respectively. We also plot two LT predictions for the MLV and EV. The lower panel of Fig. 3 illustrates the relative differences for the various samples, taken always with respect to the fiducial full one, which includes all the central halos. For the scales all the samples agree with the fiducial one down to . However, at smaller scales we can notice a clear departure of the BF in the lower number density samples from the fiducial case. The scale at which such deviations start to be noticeable, as well as the magnitude of the effect itself, depend on the number density of objects in the sample. The most diluted sample of Mpc is at characterized by median BF amplitude already higher by than for the full one, and this grows dramatically to at . This discrepancy gets the less dramatic the larger number density we consider. For a sample of Mpc, the scale at which the measured BF departs significantly from the fiducial result shrinks to , but the magnitude still can attain quite remarkable difference at the smallest scales we consider (i.e. ). The sub-sample of one order of magnitude larger number density also deviates from the fiducial case, but only at very small scales , and the relative difference reaches only for the smallest considered radius.

There is no physical reason for sparser samples to be characterized by larger bulk flow magnitudes. In particular, we expect that all samples trace the same large-scale regions of a simulated universe. The increase of the amplitude we observe is a purely statistical effect. Since the BF distribution is not Gaussian, for sparser samples the shot noise enlarges the width of the BF distribution. This effect combined with overweighted contribution of the outliers results in the observed artificial increase of the measured BF amplitude. Still, despite the fact that all the differences between the samples are contained within the 16-th and 84-th percentile variation from the median of the fiducial one, they are of a systematic nature and if ignored could be a source of a significant BF bias, especially at small scales, where in real astronomical surveys the target selection is rather non-uniform. We will discuss the implications of these systematic effects in the discussion Section IX.

Separately, we note that the LT MLV is a reasonably good prediction for the true BF at nearly all scales probed. This indicates that choosing suitable integral limits for the LT predictors (as discussed earlier) allows to properly account for the missing large-scale power.

We now consider the effects induced on our measured BF by applying various weighting schemes. The galaxy weighting prescriptions from Sec. IV are meant to roughly mimic various systematic effects present in real data. Again, we will be gauging the measured bulk flow amplitude with respect to our full sample, which constitutes an idealized fiducial case with the best sampling rates and no systematics present. For each effect we consider, we apply the specific weighting and data transformation separately from all the other effects, every time taking the fiducial full sample as a starting point. The situation as presented in our Fig. 4 looks quite the opposite to what was shown in the previous plot 3. Here, we observe that the systematic effects (if present) start to matter at large scales and grow in magnitude with scale. The effect related to the observational error modeling as in the ML method is quite easy to understand, as the error on the velocity grows linearly with scale. This taken together with the Malmquist bias (Lynden-Bell et al., 1988a, b) produces a systematic overestimation of the measured BF in relation to the full sample Nusser (2014). The scale dependence of the velocity error makes it actually quite easy to model: for the scales of , this weighting overestimates the BF by less than . At larger scales , it saturates the departure that is rather flat, as no clear scale dependence can be seen. At even larger scales the effect grows up to reaching the maximal expected effect related to the scatter of the intrinsic galaxy relation we use of .

The situation is significantly more complicated for the case of a radial selection function that is characteristic of a CosmicFlows-3 like data set. Here, there is a clear trend that grows systematically with scale, and is related to the effective depth of the sample. At scales that are above this characteristic depth, , which for our case is , the BF is grossly overestimated. At such a radial selection already biases the measurement by and this quickly grows to values of and larger for .

When we look at the geometrical selection effect of the ZoA as modeled by us, our results confirm the findings of other authors. Namely we find it to have a negligible effect on the measurements, as expected for the case of a symmetric data masking. We re-emphasize however that this is valid under the assumption of no significant nearby structures present in the ZoA, an effect that we do not investigate in the current paper.

### v.2 Cosmic Mach number

In this paper we analyze also the cosmic Mach number (CMN or , interchangeably), which, as mentioned earlier, is the ratio of the BF and peculiar velocity dispersion. The in a given sphere centered on the observer is not directly observable, however there have been some indirect methods proposed to measure the CMN Ostriker and Suto (1990); Strauss et al. (1993); Atrio-Barandela et al. (2004). Thus, we will not present and separately discuss the above mentioned sampling and weighting effects for the VD alone, but rather for the sake of brevity we show the combined effects on the actual CMN itself. This is presented in Fig. 5. Again as the reference line we take the fiducial measurement from the full sample.

The first observation to make is that the magnitude of all visible systematic effects is significantly larger for the   than it was for the BF. This is not surprising and stems from two facts. First, the VD is a much more non-linear quantity than the BF, as the former strongly depends on short-wavelengths modes; and second, the CMN is a ratio of two quantities and thus the overall effect of systematic biases and uncertainties is boosted. Moving towards more specific cases, we note that a sparse sample of Mpc leads to a strongly biased   estimate for scales . Here, the deviation from the fiducial case increases with diminishing scale, from up to more than bias in a sphere of radius . We were not able to probe the CMN for that sample on smaller scales, since the shot noise from small number counts in such a sparse sample dominates there. Even for we should be careful with interpreting our result, as the mean number of objects in such a volume is then .

For the case of modeled velocity errors, the estimator clearly provides too low a  . We have checked that this is a combination of two effects. Namely, as previously shown, the velocity errors lead to overestimation of the BF, which enters the denominator in the CMN formula. At the same time the velocity errors naturally lead to underestimation of the local . These two combined effects make such a CMN estimator, mimicking real data properties, significantly biased at all probed scales.

The situation becomes even more severe for the case of CF3-like radial selection functions. Here, both examined selections offer highly biased   estimators for all scales larger than the characteristic survey depth , and the systematic effects quickly become catastrophically large. At the estimated CMN is very close to the fiducial case; this is not a surprise, as here the radial selection is still complete (i.e. is equal to unity). Finally, it is also important to note that for the case of the  , the LT predictor does not offer a reliable estimator. This is clearly shown in Fig. 5: LT significantly underestimates the CMN for all the considered scales. As we have already assessed that the LT offers a reasonably good prediction for the BF, we then conclude that it must be the VD which is underestimated. Indeed, this is clearly the case, as was already hinted at by the results shown in Fig. 2.

## Vi Biases for Local Group observers

In this section we will test and quantify potential biases that arise in the measurements of the lowest moments of the peculiar velocity field if one neglects the fact that the related observations available to us come from a specific location in the Universe. In other words, we will compare ensemble medians of the low-order moments of the galaxy velocity field measured by unspecified observers, whom we will call random observers (RNDO), and different observers placed at specific locations which fulfill various criteria we consider to be related to the position of a Local Group (LG) analog observer. Work of Ref. Hellwing et al. (2017a) have shown that such LG observers (LGO) can exhibit highly biased local velocity correlation measurements.

To stay consistent with this previous study we will consider exactly the same selection criteria used to define a set of LG-analog observers. For clarity we give here all the essential information, referring the reader looking for more specific details or discussion to the original work. The LG is a gravitationally bound system of a dozen major galaxies with the Milky Way (MW) and its neighboring M31 as the major gravitational players. The region of 5 Mpc distance from the LG barycenter is characterized by moderate density (see e.g. Tully and Fisher, 1987, 1988; Hudson, 1993; Tully et al., 2008; Courtois et al., 2013) and a quiet flow (Sandage et al., 1972; Schlegel et al., 1994; Karachentsev et al., 2002, 2003). Located at a distance of is the Virgo cluster, whose gravitational effects extend to tens of Mpcs around us, as evident from the corresponding infall flow pattern of galaxies (Tully and Shaya, 1984; Tammann and Sandage, 1985; Lu et al., 1994; Gudehus, 1995; Karachentsev et al., 2014; Libeskind et al., 2015). The presence of such a large non-linear mass aggregation can and does have substantial impact on peculiar velocity field of the local galaxies.

To find locations of prospective LG-like observers we use the following criteria:

1. the observer is located in a MW-like host halo of mass (Busha et al., 2011; Phelps et al., 2013; Cautun et al., 2014; Guo et al., 2015),

2. the bulk velocity (of smoothed DM velocity field) within a sphere of centered on the observer is (Kogut et al., 1993),

3. the mean density contrast within the same sphere is in the range of (Karachentsev et al., 2012; Elyiv et al., 2013; Tully et al., 2014),

4. a Virgo-like cluster of mass is present at a distance from the observer (Tammann and Sandage, 1985; Mei et al., 2007).

To examine the role of individual criteria we also study results for sets of observers selected without imposing all the above constraints. The sets of observers we consider are:

LGO0

a set of the most constrained 2294 observers, each satisfying all the selection criteria 1 through 4;

LGO1

consists of 5051 candidate observers without the velocity constraint 2, but satisfying the remaining criteria 1, 3 & 4;

LGO2

includes 4978 candidates without the density contrast condition 3, but with 1, 2 & 4;

LGO3

of 4840 candidates with the conditions 2 & 3 relaxed simultaneously, i.e. meeting 1 & 4;

LGO4

a set of 6245 observers without imposing the constraint on the host halo mass 1, but with all the other criteria 24 fulfilled;

LGO5

contains 288424 candidate observers satisfying the conditions 13 but not the proximity to a Virgo-like cluster 4;

RNDO

is a list of observers with randomly selected positions in the simulation box. This set is used as a benchmark for comparison.

Since the number of prospective candidates in each set is large, to keep the sampling noise at the same level and also to speed up the calculations we will only consider 125 observers from each set. Since positions of observers are not independent of each other, we subsample the candidates by placing a grid in the simulation box and keeping only one unique observer location within each grid cell. All the results shown in this section were obtained by taking the median of the distribution for all the 125 observers in each set.

### vi.1 Bulk flow

Figure 6 illustrates the systematic effects on the median BF as measured by various observers. As the reference we take the Copernican observer of an unspecified location. In other words, we expect that the RNDO observers measure the expected cosmic mean values. Indeed, the results shown in the previous section agree with this assumption, as the BF measured for the random observers agrees well with the LT prediction (Fig. 3). The shaded region in Fig. 6 again illustrates the width of the distribution of measured bulk flows between the 16-th and the 84-th percentiles.

A quick look at the results for different non-random observers already allows us to find a striking feature: there is only one criterion really discriminatory for the results. Namely, what matters here is the proximity of a Virgo-like cluster to the observer. All LGO analogues who fulfill the latter requirement measure a BF that is systematically smaller than the cosmic mean for . This effect is around at and grows to even for scales smaller than . Additionally, we see that the LG position requirements considered without the proximity of a Virgo-like analogue also have an effect on the measured BF. Interestingly, this seems to work in the opposite direction than the other joint criteria, and an LG-analogue but no-Virgo observer would measure actually a systematically larger BF than a random one. This means that the effect of the Virgo-like object proximity is actually stronger than shown by our LG-analogues. We have used a small set of observers with just the Virgo-criterion to check that this is indeed the case.

We propose the following interpretation of these findings. The criterion that an observer should be located nearby a massive structure of a Virgo-like mass induces a constraint on the local density (hence also velocity) field when compared to a fully random observer. Such a constraint naturally lowers the scatter among observers (Hoffman and Ribak, 1992; van de Weygaert and Bertschinger, 1996), thus also the BF magnitude.

### vi.2 Velocity Dispersion and Cosmic Mach number

We now turn to the importance of observer location for the VD and   statistics. In Fig. 7 we plot the comparison of median velocity dispersions obtained for the different observers we consider. Here, we notice that the effects imposed by a Virgo-like proximity are contained to somewhat smaller ( scales than in the BF case. All our LG-analogues with a nearby cluster measure much higher (up to ) at small scales. This clearly indicates that the effect is purely driven by the presence of a massive non-linear structure of the cluster. Interestingly however, all the measurements converge to the random value at .

The effects of the observer location for the CMN statistics are illustrated in Fig. 8. Not surprisingly, it is clear that the overall LGO effect is driven mostly by the presence or absence of a nearby Virgo-analogue cluster. This amounts to LGO   bias of the order of at , which reduces to at . Thus, in the case of CMN one is concerned with an even stronger observer bias than in the BF case. This should be remembered and accounted for before any cosmological analysis of this statistics is performed.

## Vii Gravity and growth rate

In the previous sections we have observed that various systematics can significantly change the bulk flow amplitude on wide range of scales. This fact has fundamental consequences for all applications that hope to use the BF and related statistics such as   to search for non-GR signature. As an illustration, here we will compare a velocity signal from a modified gravity model against the various observational effects present in the GR case.

For our guinea pig MG model we choose the normal branch of Dvali-Gabadadze-Poratti (henceforth nDGP) model Sahni and Shtanov (2003); Lue and Starkman (2004); Schmidt (2009b), which implements the non-linear fifth-force screening (the Vainshtein mechanism) Vainshtein (1972) and can be characterized at large-scales by a nearly constant (i.e. scale-independent) enhancement of the growth rate of structures (see also Barreira et al. (2013)). Specifically we choose to take the value of the so-called crossing-over scale to be . This value represents the scale at which gravity becomes 5-dimensional in this model. The smaller this scale, the stronger deviations from GR-dynamics (due to the fifth-force) can be expected. Our choice of gives moderate modifications to GR that are characterized by a linear growth rate (the logarithmic derivative of linear density growing mode) Koyama and Maartens (2006); Li et al. (2013b); Bose et al. (2017). Except for the modified dynamics induced by the scalar field present in the nDGP model, our MG simulation shares exactly the same set-up and parameters as the fiducial GR-case. For the sake of speeding-up the numerical computations we have employed the Truncated DGP method described in details in Barreira et al. (2015). The speed-up is obtained at the expense of the resolution of the scalar-field spatial fluctuations, solving of which was truncated beyond the 4-th mesh refinement level. This sets the resolution of the scalar force at , which is still considerably smaller than the smallest halos we consider. As we use the same initial conditions for both GR and nDGP, the large-scale cosmic variance effects should be of the same magnitude in both runs (see also Hellwing et al. (2017b)), and the observed discrepancies should reflect the differences in the underlying gravitational dynamics.

Figure 9 compares the BF measured by two Copernican observers, one in GR and one in the nDGP model (marked as MG) versus the amplitudes expected in the GR case with different systematic effects. For the sake of brevity, we choose to compare with only the strongest systematics elucidated in the previous Section. In particular, we show the LGO0 and LGO5 signals, as well as RNDO observers with sparse sampling of hMpc. For , the MG bulk flow is enhanced by , as one can expect from the LT prediction of Eqn. (16). This potentially observable effect can be easily obscured by various systematics that have larger magnitudes on the same scales. Specifically, we see that realistic modeling of the Local Group analogue observers, which includes the effects of the Virgo cluster proximity, gives opposite sign to the MG enhancement. Thus, in the worst case scenario, we could have a conspiracy, where a BF signal for a LGO0 observer in an MG universe would look like a BF expected for a RNDO observer in the GR universe. On the other hand, the signal expected for a LGO observer modeled without a Virgo-like cluster presence can mimic the scale-dependence and amplitude of a RNDO MG signal. For a very sparse sample, these two observations would be dwarfed by a systematic effect that on small scales () can be by a factor of a few times larger than what we can expect for a reasonably mild MG model enhancement.

The main merit of our work here is to systematically study potential biases of low-order velocity measurements, but it is illustrative to compare the scales and amplitudes of the effects we report with some measurements reported in the literature. We have selected arbitrarily seven such measurements and marked them in Fig. 9. We show results from Branchini et al.(Branchini et al., 2012), Hoffman et al.(Hoffman et al., 2015), Hong et al.(Hong et al., 2014), Scrimgeour et al.(Scrimgeour et al., 2016), Lavaux et al.(Lavaux et al., 2013), Nusser&Davis(Nusser and Davis, 2011) and Feldman et al.(Feldman et al., 2010). The methods and datasets used in these references vary significantly, so this collection is a fair representation of approaches and data used currently in peculiar velocity studies. Except for Refs. (Hoffman et al., 2015) and (Feldman et al., 2010), all the results are consistent within -th and th percentiles with median of random and LGO observers and even with the MG model. If we took the size of the error bars reported by those authors at their face-values, some of our results (such as the MG model) would be marginally inconsistent with that data. However, we clearly see that the variance added by the systematic effects will boost the reported error bars significantly.

The results shown here can have potentially profound repercussions, as it would seem that the lower-order velocity statistics are plagued by potentially overwhelming systematic effects that can completely obscure even relatively strong () deviations from the GR case. We shall discuss the implication of these findings in the next section.

## Viii Summary

In this paper our main aim was to methodically check various possible systematic effects that could affect the measured values of the bulk flow, peculiar velocity dispersion, and cosmic Mach number. Peculiar velocities of galaxies strongly depend on the underlying cosmic parameters, such as the logarithmic growth rate () and the non-relativistic matter energy density (). Velocity data are prone to large uncertainties stemming from the intrinsic scatter of various empirical galaxy scaling relations used to measure redshift-independent distances, and consequently, infer peculiar velocities. The latter are additionally affected by such issues as non-Gaussian errors and non-linear (i.e. virial) contributions. Many methods have been proposed and implemented to deal with these issues. However, once the velocity data had been corrected for intrinsic errors and Malmquist biases, it was commonly assumed that the relevant statistics could be directly related to the underlying cosmology using theoretical modeling (such as linear perturbation theory). This assumed advantage was one of the main arguments for using the galaxy velocities as alternative cosmological probes. In our analysis we have revisited this assumption, and our results indicate that there are many systematic effects that need to be accurately modeled and accounted for, in order to infer cosmological parameters from low-order velocity statistics in an unbiased manner.

Below we summarize and comment on all our important findings and their implications.

• Perturbation theory estimators:

1. The results encapsulated in Fig. 2 show that strong modulations of the density power spectrum amplitude lead to only mild variations in BF and VD. More precisely, we have found that the changes in both statistics are roughly proportional to changes in and .

2. Using the non-linear velocity divergence power spectrum instead of the non-linear density has a strong effect on the predicted , and therefore also on the  . This is because the magnitudes of the effects induced in the non-linear regime are opposite in those two spectra. Namely, at small scales the non-linear takes smaller values than the linear theory prediction, while the non-linear has actually a boosted amplitude w.r.t. the linear theory. At small scales in the non-linear regime the motions of galaxies lose the character of a potential flow. This reflects significant growth of vorticity and velocity dispersion due to shell-crossing (Pichon and Bernardeau, 1999; Pueblas and Scoccimarro, 2009; Libeskind et al., 2014; Cusin et al., 2017). For that reason it is important to take into account and model properly these non-linear effects, in order to get a more realistic perturbation theory prediction for the CMN.

• Sampling rate effects:

1. At small scales (i.e. ) we found a significant effect on BF magnitude from under-sampling. For increasingly diluted samples the inferred BF is biased towards higher values when compared to our fiducial full sample, and the effect is the larger, the smaller sample density. This is a statistical effect induced by the non-Gaussian distribution of BF magnitudes and increased shot-noise contribution due to sparse sampling. For the sparsest sample with Mpc, the bias could attain or more. This indicates that one should carefully account for this effect when low-order velocity statistics are measured from very sparse samples like supernovae Weyant et al. (2011); Turnbull et al. (2012); Feindt et al. (2013); Huterer et al. (2015).

2. The down-sampling of the data we have performed is essentially equivalent to weighting the data by the local value of the density correlation function, since galaxies are biased tracers and form preferentially in higher-density regions. Such regions are characterized by higher values of local bulk flows.

3. Our results show that the MLV estimator from linear perturbation theory is in good agreement with the -body data for a random-location BF observer.

• Selection effects:

1. Weighting halo velocities according to velocity errors (like in the ML method) leads to overestimation of the BF. This effect grows with distance from the observer, but it saturates to a maximum value that is close to the considered typical scatter of the distance indicator error .

2. For a radial selection function of a CF3-like survey form, it is evident that the limited depth of the catalog is reflected in the measured BF. The BF is formally an integral over a sphere, but for a realistic survey of discrete galaxies, it becomes a sum over concentric spherical shells of growing radius. The radial selection function is effectively down-weighting the outer shells, and thus the BF value derived from inferior spheres is spuriously propagated to larger scales. This is a strong effect, which indicates that BF measurements from scales comparable or larger than the characteristic depth of a given catalog should be interpreted with care.

3. A symmetric sky-map angular incompleteness with an opening angle of – such as the Zone of Avoidance – has a negligible effect on the inferred BF.

4. All systematic effects related to galaxy selection are much more pronounced in the CMN statistic. In particular, both sparse sampling and velocity errors induce a significant   bias for scales .

5. A radial selection function of the CF3-like form induces a catastrophically large CMN over-prediction for scales larger than the given survey characteristic depth.

• Observer location effects:

1. For all our low-order statistics the most important LGO-analogue criterion is the proximity of a Virgo-like cluster.

2. Local Group-analogue observers measure systematically smaller BF amplitudes than the cosmic mean (i.e. a random observer) on scales up to . This systematic attains at and grows to for .

3. A no-Virgo observer (i.e. LGO5) at the same time exhibits an opposite bias, inferring a BF that is larger by on similar scales than measured by the RNDO observers.

4. For the VD the LG-observer bias is contained to somewhat smaller scales , but its magnitude reaches a quite dramatic value of up to . This effect is purely driven by the proximity of a Virgo-analogue, as a no-Virgo observer measures VD values compatible with RNDO observers. This indicates that an infall region around a massive cluster is significantly heating up the local velocity field.

5. The effects observed for BF and VD combine into a biased   for LGO observers, which manifests itself as under-evaluation at , and still takes too small a value at .

• Growth rate and gravity

1. The effect of an increased growth rate observed in a representative MG model is degenerate with the specific bias induced by a LGO observer, and the Virgo-like cluster proximity in particular.

2. Similarly, a non-Virgo LGO-like BF signal appears very similar to the MG signal for a RNDO observer.

3. Finally, we notice that the BF magnitude increase observed in GR for a sparse sample of Mpc is stronger than the non-GR effect of our MG model.

## Ix Discussion

The above summary of all our important findings regarding various systematic effects that impact low-order moments of the galaxy velocity field, underlines a number of crucial observations. The linear theory predictions (for MLV) obtained using Eqns.(14-18) render quite accurate values of BF for a “cosmic mean” (i.e. Copernican) observer in the case of high-density clean data. Thus, they can be used as a first order prediction for the case when all systematic effects can be ignored, or when there are no computer simulations to be compared with. However, we caution that whenever one wants to compare such LT predictions to the real data, one needs to remember that the local galaxy velocity field is biased w.r.t. the LT prediction at small scales and for sparse or radially incomplete samples.

Currently, the velocity data is sparse and noisy; however, in the near future they will increase significantly in volume. There is also hope for better modeling and understanding of galaxy intrinsic scaling relations, which can lead to further suppression of individual velocity errors. The surveys such as TAIPAN da Cunha et al. (2017) or WALLABY Duffy et al. (2012) will lead to an increase both in richness as well as depth of galaxy peculiar velocity catalogs. In addition, the possibility to obtain transverse velocities from surveys like GAIA (de Bruijne, 2012) or LSST (Telescope, ) could result in additional largely unbiased velocity measurements for the nearby galaxies (for a dedicated discussion see Nusser et al. (2012)). In this new era of velocity data, the linear theory predictions for and related statistics will be too inaccurate to be used for model testing and data analysis, except for the large scales (i.e. ), where the precision and data abundance will continue to be poor.

The various observer-independent systematic effects surfacing strongly in our analysis suggest that the bulk flow amplitude and related measurements at small distances should be carefully reanalyzed and compared with predictions based on galaxy-mock catalogs. The higher BF values reported by Refs. (Branchini et al., 2012; Hong et al., 2014; Kashlinsky et al., 2008; Ma and Scott, 2013; Ma and Pan, 2014; Hoffman et al., 2015) might be a signature of biases induced by sparse sampling and radial selection.

The importance of proper modeling of such non-linear effects is of paramount importance for the cosmic Mach number predictions. This was already emphasized by Ref. Agarwal and Feldman (2013) for the case of improving the LT predictions by using the non-linear matter density power spectrum rather than the linear one. Here, our analysis adds further that sparse sampling induces a very strong effect on the VD and thus on the resulting CMN. In addition, other systematics such as radial selection and velocity errors affect   to a much stronger extent than BF. This suggests in particular that one should aim at using the richest possible galaxy samples when considering CMN measurements.

The proximity of a Virgo-like cluster to a Local-Group-like observer is equally significant and needs to be considered as an additional important contribution to the local bulk flow. If this effect is not properly accounted for in the BF analysis, it will result in an additional non-Gaussian systematic for the BF measured on scales .

Finally, the combination of all the aforementioned systematic effects, if not accounted for carefully, can lead to strong degeneracies of the cosmological signals encoded in galaxy velocities and in their low-order moments. We have clearly demonstrated that the signal of a non-GR cosmological model, such as the nDGP modified gravity we considered, that employs a moderately strong modification to the cosmic growth rate of structures, can be easily absorbed by the non-trivial systematics effects we studied. In light of this evidence, analyses such as for example Ref. Seiler and Parkinson (2016), where the BF deviations induced by modified gravity were studied, should be definitely revisited.

All this is a source of potentially major concern, as the cosmic velocity data offers, at least in principle, a model-independent way to constrain growth rate and gravity on cosmic scales (Hellwing et al., 2014; Feix et al., 2014; Hoffman et al., 2015; Nusser, 2016; Hoffman et al., 2016). Other authors have also shown that current and future data show promise to become competitive cosmological probes (Koda et al., 2014; Howlett et al., 2017). In principle, this can still be achieved. However, the results presented here clearly indicate that all the various systematic effects need to be carefully addressed and accounted for, before any high-accuracy cosmological analysis can be performed. This is especially important that in the near future, the amount of peculiar velocity data is expected to significantly increase, and therefore systematics will likely dominate over statistical errors in relevant studies. In this context, using dedicated computer simulations employing constrained local density (and velocity) realizations (such as (Gottloeber et al., 2010; Courtois and Tully, 2012; Heß et al., 2013; Sorce et al., 2014; Leclercq et al., 2015; Sawala et al., 2016; Sorce et al., 2016; Leclercq et al., 2017; Desmond et al., 2018)) look very promising.

## Acknowledgments

The authors are very grateful to Baojiu Li for inspiring discussions and for providing the ECOSMOG code that was used to run the simulations in this paper. Enzo Branchini, Martin Feix and Adi Nusser are also acknowledged for stimulating discussions and for careful reading of the manuscript. WAH is supported by an Individual Fellowship of the Marie Skłodowska-Curie Actions and therefore acknowledges that this project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 748525. MB is supported by the Netherlands Organization for Scientific Research, NWO, through grant number 614.001.451, and by the Polish National Science Center under contract UMO-2012/07/D/ST9/02785. This project has also benefited from numerical computations performed at the Interdisciplinary Centre for Mathematical and Computational Modelling (ICM) University of Warsaw under grants no GA67-17 and GA65-30. This work used the COSMA Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk. This equipment was funded by a BIS National E-infrastructure capital grant ST/K00042X/1, DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure.

## Appendix

Here we will investigate how and on what scales the limited simulation volume affects our measurements. As mentioned in the main text, the cosmic velocity field is characterized by a large correlation length. This means that the convergence of velocity moments is slower than in the case of the density field. For that reason we can expect that scales which are normally considered as converged will be still affected by missing large-scale modes in our simulations. To assess this we have conducted a series of auxiliary approximated simulations varying the box size. For this we employ the COmoving Lagrangian Accelerator (COLA) method (Tassev et al., 2013). The parallel implementation of the COLA algorithm (called PI-COLA, see (Howlett et al., 2015)) allows to run large simulations at a reduced computational cost, with the trade-off of limited spatial and temporal resolution. We are however here interested in the effect of the missing large-scale modes, thus the scales which we will study, namely , are large enough to be fully resolved by the PI-COLA method. In particular, we have used the publicly available optimized branch of the COLA family, the MG-COLA, introduced by Ref. (Winther et al., 2017).

We ran three simulations, all containing volume elements, with three boxes: and on a side. The simulations are set to use the same cosmology as our main N-body runs used in this work. We process the simulation outputs at in the same manner as our full N-body simulations. We use the final ROCKSTAR halo catalogs as our input data. The COLA method is known to bias weakly the resulting halo velocities. However, this bias is small (up to for our case) and concerns mostly the small-scale halo velocity field (see more in Ref. (Munari et al., 2017)). Since our primary concern here is to study the effects of missing large-scale power, we are confident that halo catalogs obtained via the simplified COLA method are suitable for our purpose.

In FIG. 10 we show probability distribution functions of  magnitudes for our COLA runs computed and binned for spheres of three radii: in the left-hand panel, for the middle one, and for the right-hand panel. The PDF for the fiducial run of box is illustrated by filled purple boxes, while two consecutively smaller simulations are depicted by open green () and blue () boxes. For comparison we also plot the relevant PDFs for the full -body simulation used in the paper. The immediate impression is that all the results for the smallest box are significantly biased w.r.t. the fiducial case. Both the median of the distribution (), as well as its spread () 1 are visibly smaller for all three radii. The corresponding relative differences taken w.r.t the fiducial case here (i.e. the largest box) are () and () for the median and spread of the PDF at () and and at . This is a clear manifestation of the lacking large-scale power in the simulation box. Thus we can, not surprisingly, conclude that the halo velocity field is not converged at those scales in the box. The situation for the medium box is much better though. Here, the medians are off by only () at (), while the corresponding distribution widths are smaller by (). However, at the largest radius of the biases grow to for the median and for the width. Our full -body simulations use a box, and it is reassuring to find biases of their adequate distributions to lie between COLA and boxes. At the median is biased by more than already. This shows that all the results in this paper for spheres larger than are noticeably affected by missing large-scale power.

Obtaining a reliable and accurate absolute bulk flow (and corresponding dispersion) magnitude is of paramount importance when one wants to compare it with astronomical data and use such a comparison for parameter constraints. In this paper however, we are more interested in some specific effects that affect the velocity field statistics in a systematic way. Thus, to study to what scales we can trust our results we present FIG. 11. Here, we compare size of relative differences (taken always w.r.t. the fiducial unweighted sample) for three systematic effects: ML velocity error weights (left panel), CF-3-like (with ) radial selection function (middle panel) and the effect of the Zone of Avoidance (right panel). Reassuringly, we denote that in all three cases the scale-dependence as well as  magnitude are very similar (within the sampling error) for the three COLA run and our full -body. The largest differences appear for the individual velocity error weights. Here, the -body results for are consistently below the COLA line. This indicates that for the case of this systematic effect and at those scales our results render only the lower-bound, and amore realistic modeling will probably foster larger effects.

Finally, we can report that the   distributions (both from COLA and -body) at all probed scales are deviating significantly from a Gaussian, with typical skewness () and kurtosis () taking values .

### Footnotes

1. For our purpose here, for a measure of the distribution spread (i.e. ) we use the half-width between the 16-th and 84-th percentiles.

### References

1. G. Hinshaw, D. Larson, E. Komatsu, D. N. Spergel, C. L. Bennett, J. Dunkley, M. R. Nolta, M. Halpern, R. S. Hill, N. Odegard, L. Page, K. M. Smith, J. L. Weiland, B. Gold, N. Jarosik, A. Kogut, M. Limon, S. S. Meyer, G. S. Tucker, E. Wollack,  and E. L. Wright, ApJS 208, 19 (2013)arXiv:1212.5226 .
2. Planck Collaboration, P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, J. G. Bartlett, N. Bartolo, E. Battaner, R. Battye, K. Benabed, A. Benoît, A. Benoit-Lévy, J.-P. Bernard, M. Bersanelli, P. Bielewicz, J. J. Bock,  and et al., A&A 594, A13 (2016)arXiv:1502.01589 .
3. J. Yang, M. S. Turner, D. N. Schramm, G. Steigman,  and K. A. Olive, ApJ 281, 493 (1984).
4. T. P. Walker, G. Steigman, H.-S. Kang, D. M. Schramm,  and K. A. Olive, ApJ 376, 51 (1991).
5. W. J. Percival, C. M. Baugh, J. Bland-Hawthorn, T. Bridges, R. Cannon, S. Cole, M. Colless, C. Collins, W. Couch, G. Dalton, R. De Propris, S. P. Driver, G. Efstathiou, R. S. Ellis, C. S. Frenk, K. Glazebrook, C. Jackson, O. Lahav, I. Lewis, S. Lumsden, S. Maddox, S. Moody, P. Norberg, J. A. Peacock, B. A. Peterson, W. Sutherland,  and K. Taylor, MNRAS 327, 1297 (2001)astro-ph/0105252 .
6. M. Tegmark, M. R. Blanton, M. A. Strauss, F. Hoyle, D. Schlegel, R. Scoccimarro, M. S. Vogeley, D. H. Weinberg, I. Zehavi, A. Berlind, T. Budavari, A. Connolly, D. J. Eisenstein, D. Finkbeiner, J. A. Frieman, J. E. Gunn, A. J. S. Hamilton, L. Hui, B. Jain, D. Johnston, S. Kent, H. Lin, R. Nakajima, R. C. Nichol, J. P. Ostriker, A. Pope, R. Scranton, U. Seljak, R. K. Sheth, A. Stebbins, A. S. Szalay, I. Szapudi, L. Verde, Y. Xu, J. Annis, N. A. Bahcall, J. Brinkmann, S. Burles, F. J. Castander, I. Csabai, J. Loveday, M. Doi, M. Fukugita, J. R. Gott, III, G. Hennessy, D. W. Hogg, Ž. Ivezić, G. R. Knapp, D. Q. Lamb, B. C. Lee, R. H. Lupton, T. A. McKay, P. Kunszt, J. A. Munn, L. O’Connell, J. Peoples, J. R. Pier, M. Richmond, C. Rockosi, D. P. Schneider, C. Stoughton, D. L. Tucker, D. E. Vanden Berk, B. Yanny, D. G. York,  and SDSS Collaboration, ApJ 606, 702 (2004)astro-ph/0310725 .
7. F. Beutler, H.-J. Seo, S. Saito, C.-H. Chuang, A. J. Cuesta, D. J. Eisenstein, H. Gil-Marín, J. N. Grieb, N. Hand, F.-S. Kitaura, C. Modi, R. C. Nichol, M. D. Olmstead, W. J. Percival, F. Prada, A. G. Sánchez, S. Rodriguez-Torres, A. J. Ross, N. P. Ross, D. P. Schneider, J. Tinker, R. Tojeiro,  and M. Vargas-Magaña, MNRAS 466, 2242 (2017)arXiv:1607.03150 .
8. A. G. Riess, A. V. Filippenko, P. Challis, A. Clocchiatti, A. Diercks, P. M. Garnavich, R. L. Gilliland, C. J. Hogan, S. Jha, R. P. Kirshner, B. Leibundgut, M. M. Phillips, D. Reiss, B. P. Schmidt, R. A. Schommer, R. C. Smith, J. Spyromilio, C. Stubbs, N. B. Suntzeff,  and J. Tonry, AJ 116, 1009 (1998)arXiv:astro-ph/9805201 .
9. S. Perlmutter, G. Aldering, G. Goldhaber, R. A. Knop, P. Nugent, P. G. Castro, S. Deustua, S. Fabbro, A. Goobar, D. E. Groom, I. M. Hook, A. G. Kim, M. Y. Kim, J. C. Lee, N. J. Nunes, R. Pain, C. R. Pennypacker, R. Quimby, C. Lidman, R. S. Ellis, M. Irwin, R. G. McMahon, P. Ruiz-Lapuente, N. Walton, B. Schaefer, B. J. Boyle, A. V. Filippenko, T. Matheson, A. S. Fruchter, N. Panagia, H. J. M. Newberg, W. J. Couch,  and The Supernova Cosmology Project, ApJ 517, 565 (1999)arXiv:astro-ph/9812133 .
10. W. J. Percival, B. A. Reid, D. J. Eisenstein, N. A. Bahcall, T. Budavari, J. A. Frieman, M. Fukugita, J. E. Gunn, Ž. Ivezić, G. R. Knapp, R. G. Kron, J. Loveday, R. H. Lupton, T. A. McKay, A. Meiksin, R. C. Nichol, A. C. Pope, D. J. Schlegel, D. P. Schneider, D. N. Spergel, C. Stoughton, M. A. Strauss, A. S. Szalay, M. Tegmark, M. S. Vogeley, D. H. Weinberg, D. G. York,  and I. Zehavi, MNRAS 401, 2148 (2010)arXiv:0907.1660 [astro-ph.CO] .
11. D. H. Weinberg, M. J. Mortonson, D. J. Eisenstein, C. Hirata, A. G. Riess,  and E. Rozo, Phys. Rep. 530, 87 (2013)arXiv:1201.2434 .
12. M. A. Strauss and J. A. Willick, Phys. Rep. 261, 271 (1995)astro-ph/9502079 .
13. J. Koda, C. Blake, T. Davis, C. Magoulas, C. M. Springob, M. Scrimgeour, A. Johnson, G. B. Poole,  and L. Staveley-Smith, MNRAS 445, 4267 (2014)arXiv:1312.1022 .
14. C. Howlett, L. Staveley-Smith,  and C. Blake, MNRAS 464, 2517 (2017)arXiv:1609.08247 .
15. A. Nusser and M. Davis, ApJ 736, 93 (2011)arXiv:1101.1650 .
16. M. J. Hudson and S. J. Turnbull, ApJ 751, L30 (2012)arXiv:1203.4814 .
17. T. R. Lauer and M. Postman, ApJ 425, 418 (1994).
18. R. Giovanelli, M. P. Haynes, G. Wegner, L. N. da Costa, W. Freudling,  and J. J. Salzer, ApJ 464, L99 (1996)astro-ph/9610130 .
19. D. A. Dale, R. Giovanelli, M. P. Haynes, L. E. Campusano, E. Hardy,  and S. Borgani, ApJ 510, L11 (1999)astro-ph/9810467 .
20. M. J. Hudson, R. J. Smith, J. R. Lucey, D. J. Schlegel,  and R. L. Davies, ApJ 512, L79 (1999)astro-ph/9901001 .
21. R. Watkins, H. A. Feldman,  and M. J. Hudson, MNRAS 392, 743 (2009)arXiv:0809.4041 .
22. Y. Hoffman, H. M. Courtois,  and R. B. Tully, MNRAS 449, 4494 (2015)arXiv:1503.05422 .
23. R. Watkins and H. A. Feldman, MNRAS 447, 132 (2015)arXiv:1407.6940 .
24. A. Kashlinsky, F. Atrio-Barandela, D. Kocevski,  and H. Ebeling, ApJ 686, L49 (2008)arXiv:0809.3734 .
25. S. J. Osborne, D. S. Y. Mak, S. E. Church,  and E. Pierpaoli, ApJ 737, 98 (2011)arXiv:1011.2781 .
26. Planck Collaboration, P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. Balbi, A. J. Banday, R. B. Barreiro, E. Battaner, K. Benabed, A. Benoit-Lévy, J.-P. Bernard, M. Bersanelli, P. Bielewicz, I. Bikmaev, J. Bobin, J. J. Bock, A. Bonaldi, J. R. Bond, J. Borrill, F. R. Bouchet, C. Burigana, R. C. Butler, P. Cabella, J.-F. Cardoso, A. Catalano, A. Chamballu, L.-Y. Chiang, G. Chon, P. R. Christensen, D. L. Clements, S. Colombi, L. P. L. Colombo, B. P. Crill, F. Cuttaia, A. Da Silva, H. Dahle, R. D. Davies, R. J. Davis, P. de Bernardis, G. de Gasperis, G. de Zotti, J. Delabrouille, J. Démoclès, J. M. Diego, K. Dolag, H. Dole, S. Donzelli, O. Doré, U. Dörl, M. Douspis, X. Dupac, T. A. Enßlin, F. Finelli, I. Flores-Cacho, O. Forni, M. Frailis, M. Frommert, S. Galeotta, K. Ganga, R. T. Génova-Santos, M. Giard, G. Giardino, J. Gonzáalez-Nuevo, A. Gregorio, A. Gruppuso, F. K. Hansen, D. Harrison, C. Hernández-Monteagudo, D. Herranz, S. R. Hildebrandt, E. Hivon, W. A. Holmes, W. Hovest, K. M. Huffenberger, G. Hurier, T. R. Jaffe, A. H. Jaffe, J. Jasche, W. C. Jones, M. Juvela, E. Keihánen, R. Keskitalo, I. Khamitov, T. S. Kisner, J. Knoche, M. Kunz, H. Kurki-Suonio, G. Lagache, A. Lähteenmäki, J.-M. Lamarre, A. Lasenby, C. R. Lawrence, M. Le Jeune, R. Leonardi, P. B. Lilje, M. Linden-Vørnle, M. López-Caniego, J. F. Macías-Pérez, D. Maino, D. S. Y. Mak, N. Mandolesi, M. Maris, F. Marleau, E. Martínez-González, S. Masi, S. Matarrese, P. Mazzotta, A. Melchiorri, J.-B. Melin, L. Mendes, A. Mennella, M. Migliaccio, S. Mitra, M.-A. Miville-Deschênes, A. Moneti, L. Montier, G. Morgante, D. Mortlock, A. Moss, D. Munshi, J. A. Murphy, P. Naselsky, F. Nati, P. Natoli, C. B. Netterfield, H. U. Nørgaard-Nielsen, F. Noviello, D. Novikov, I. Novikov, S. Osborne, L. Pagano, D. Paoletti, O. Perdereau, F. Perrotta, F. Piacentini, M. Piat, E. Pierpaoli, D. Pietrobon, S. Plaszczynski, E. Pointecouteau, G. Polenta, L. Popa, T. Poutanen, G. W. Pratt, S. Prunet, J.-L. Puget, S. Puisieux, J. P. Rachen, R. Rebolo, M. Reinecke, M. Remazeilles, C. Renault, S. Ricciardi, M. Roman, J. A. Rubiño-Martín, B. Rusholme, M. Sandri, G. Savini, D. Scott, L. Spencer, R. Sunyaev, D. Sutton, A.-S. Suur-Uski, J.-F. Sygnet, J. A. Tauber, L. Terenzi, L. Toffolatti, M. Tomasi, M. Tristram, M. Tucci, L. Valenziano, J. Valiviita, B. Van Tent, P. Vielva, F. Villa, N. Vittorio, L. A. Wade, N. Welikala, D. Yvon, A. Zacchei, J. P. Zibin,  and A. Zonca, A&A 561, A97 (2014a)arXiv:1303.5090 .
27. J. Colin, R. Mohayaee, S. Sarkar,  and A. Shafieloo, MNRAS 414, 264 (2011)arXiv:1011.6292 .
28. A. Nusser, E. Branchini,  and M. Davis, ApJ 735, 77 (2011)arXiv:1102.4189 .
29. A. Weyant, M. Wood-Vasey, L. Wasserman,  and P. Freeman, ApJ 732, 65 (2011)arXiv:1103.1603 .
30. E. Branchini, M. Davis,  and A. Nusser, MNRAS 424, 472 (2012)arXiv:1202.5206 .
31. K. Mody and A. Hajian, ApJ 758, 4 (2012)arXiv:1202.1339 [astro-ph.CO] .
32. S. J. Turnbull, M. J. Hudson, H. A. Feldman, M. Hicken, R. P. Kirshner,  and R. Watkins, MNRAS 420, 447 (2012)arXiv:1111.0631 .
33. U. Feindt, M. Kerschhaggl, M. Kowalski, G. Aldering, P. Antilogus, C. Aragon, S. Bailey, C. Baltay, S. Bongard, C. Buton, A. Canto, F. Cellier-Holzem, M. Childress, N. Chotard, Y. Copin, H. K. Fakhouri, E. Gangler, J. Guy, A. Kim, P. Nugent, J. Nordin, K. Paech, R. Pain, E. Pecontal, R. Pereira, S. Perlmutter, D. Rabinowitz, M. Rigault, K. Runge, C. Saunders, R. Scalzo, G. Smadja, C. Tao, R. C. Thomas, B. A. Weaver,  and C. Wu, A&A 560, A90 (2013)arXiv:1310.4184 .
34. G. Lavaux, N. Afshordi,  and M. J. Hudson, MNRAS 430, 1617 (2013)arXiv:1207.1721 .
35. Y.-Z. Ma and D. Scott, MNRAS 428, 2017 (2013)arXiv:1208.2028 .
36. B. Rathaus, E. D. Kovetz,  and N. Itzhaki, MNRAS 431, 3678 (2013)arXiv:1301.7710 .
37. M. Feix, A. Nusser,  and E. Branchini, J. Cosmology Astropart. Phys 9, 019 (2014)arXiv:1405.6710 .
38. T. Hong, C. M. Springob, L. Staveley-Smith, M. I. Scrimgeour, K. L. Masters, L. M. Macri, B. S. Koribalski, D. H. Jones,  and T. H. Jarrett, MNRAS 445, 402 (2014)arXiv:1409.0287 .
39. Y.-Z. Ma and J. Pan, MNRAS 437, 1996 (2014)arXiv:1311.6888 .
40. S. Appleby, A. Shafieloo,  and A. Johnson, ApJ 801, 76 (2015)arXiv:1410.5562 .
41. D. Huterer, D. L. Shafer,  and F. Schmidt, J. Cosmology Astropart. Phys 12, 033 (2015)arXiv:1509.04708 .
42. M. I. Scrimgeour, T. M. Davis, C. Blake, L. Staveley-Smith, C. Magoulas, C. M. Springob, F. Beutler, M. Colless, A. Johnson, D. H. Jones, J. Koda, J. R. Lucey, Y.-Z. Ma, J. Mould,  and G. B. Poole, MNRAS 455, 386 (2016)arXiv:1511.06930 .
43. C. M. Springob, T. Hong, L. Staveley-Smith, K. L. Masters, L. M. Macri, B. S. Koribalski, D. H. Jones, T. H. Jarrett, C. Magoulas,  and P. Erdoğdu, MNRAS 456, 1886 (2016)arXiv:1511.04849 .
44. A. Nusser, ApJ 795, 3 (2014)arXiv:1405.6271 .
45. P. Andersen, T. M. Davis,  and C. Howlett, MNRAS 463, 4083 (2016)arXiv:1609.04022 .
46. A. Nusser, MNRAS 455, 178 (2016)arXiv:1508.03860 .
47. E. da Cunha, A. M. Hopkins, M. Colless, E. N. Taylor, C. Blake, C. Howlett, C. Magoulas, J. R. Lucey, C. Lagos, K. Kuehn, Y. Gordon, D. Barat, F. Bian, C. Wolf, M. J. Cowley, M. White, I. Achitouv, M. Bilicki, J. Bland-Hawthorn, K. Bolejko, M. J. I. Brown, R. Brown, J. Bryant, S. Croom, T. M. Davis, S. P. Driver, M. D. Filipovic, S. R. Hinton, M. Johnston-Hollitt, D. H. Jones, B. Koribalski, D. Kleiner, J. Lawrence, N. Lorente, J. Mould, M. S. Owers, K. Pimbblet, C. G. Tinney, N. F. H. Tothill,  and F. Watson, ArXiv e-prints  (2017), arXiv:1706.01246 .
48. A. R. Duffy, M. J. Meyer, L. Staveley-Smith, M. Bernyk, D. J. Croton, B. S. Koribalski, D. Gerstmann,  and S. Westerlund, MNRAS 426, 3385 (2012)arXiv:1208.5592 .
49. R. B. Tully, H. M. Courtois,  and J. G. Sorce, AJ 152, 50 (2016)arXiv:1605.01765 .
50. W. A. Hellwing, A. Nusser, M. Feix,  and M. Bilicki, MNRAS 467, 2787 (2017a)arXiv:1609.07120 .
51. M. Li, J. Pan, L. Gao, Y. Jing, X. Yang, X. Chi, L. Feng, X. Kang, W. Lin, G. Shan, L. Wang, D. Zhao,  and P. Zhang, ApJ 761, 151 (2012a)arXiv:1207.5338 .
52. J. P. Ostriker and Y. Suto, ApJ 348, 378 (1990).
53. Y.-Z. Ma, J. P. Ostriker,  and G.-B. Zhao, J. Cosmology Astropart. Phys 6, 026 (2012)arXiv:1106.3327 [astro-ph.CO] .
54. Y. Suto, R. Cen,  and J. P. Ostriker, ApJ 395, 1 (1992).
55. M. A. Strauss, R. Cen,  and J. P. Ostriker, ApJ 408, 389 (1993).
56. K. Nagamine, J. P. Ostriker,  and R. Cen, ApJ 553, 513 (2001)astro-ph/0010253 .
57. F. Atrio-Barandela, A. Kashlinsky,  and J. P. Mücket, ApJ 601, L111 (2004)astro-ph/0312394 .
58. S. Agarwal and H. A. Feldman, MNRAS 432, 307 (2013)arXiv:1301.1039 .
59. B. Li, W. A. Hellwing, K. Koyama, G.-B. Zhao, E. Jennings,  and C. M. Baugh, MNRAS 428, 743 (2013a)arXiv:1206.4317 [astro-ph.CO] .
60. W. A. Hellwing, A. Barreira, C. S. Frenk, B. Li,  and S. Cole, Physical Review Letters 112, 221102 (2014)arXiv:1401.0706 .
61. M. F. Ivarsen, P. Bull, C. Llinares,  and D. Mota, A&A 595, A40 (2016)arXiv:1603.03072 .
62. G. Tormen, L. Moscardini, F. Lucchin,  and S. Matarrese, ApJ 411, 16 (1993).
63. B. Li, G.-B. Zhao, R. Teyssier,  and K. Koyama, J. Cosmology Astropart. Phys 1, 51 (2012b)arXiv:1110.1379 [astro-ph.CO] .
64. Planck Collaboration, P. A. R. Ade, N. Aghanim, C. Armitage-Caplan, M. Arnaud, M. Ashdown, F. Atrio-Barandela, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, J. G. Bartlett, E. Battaner, K. Benabed, A. Benoît, A. Benoit-Lévy, J.-P. Bernard, M. Bersanelli, P. Bielewicz, J. Bobin,  and et al., A&A 571, A16 (2014b)arXiv:1303.5076 .
65. W. A. Hellwing, M. Schaller, C. S. Frenk, T. Theuns, J. Schaye, R. G. Bower,  and R. A. Crain, MNRAS 461, L11 (2016)arXiv:1603.03328 .
66. R. Teyssier, A&A 385, 337 (2002)arXiv:astro-ph/0111367 .
67. G. Dvali, G. Gabadadze,  and M. Porrati, Physics Letters B 485, 208 (2000)hep-th/0005016 .
68. J.-N. Ye, H. Guo, Z. Zheng,  and I. Zehavi, ApJ 841, 45 (2017)arXiv:1705.02071 .
69. P. S. Behroozi, R. H. Wechsler,  and H.-Y. Wu, ApJ 762, 109 (2013)arXiv:1110.4372 [astro-ph.CO] .
70. C. M. Springob, K. L. Masters, M. P. Haynes, R. Giovanelli,  and C. Marinoni, ApJS 172, 599 (2007).
71. C. M. Springob, C. Magoulas, M. Colless, J. Mould, P. Erdoğdu, D. H. Jones, J. R. Lucey, L. Campbell,  and C. J. Fluke, MNRAS 445, 2677 (2014)arXiv:1409.6161 .
72. P. J. E. Peebles, The large-scale structure of the universe (Research supported by the National Science Foundation. Princeton, N.J., Princeton University Press, 1980. 435 p., 1980).
73. O. Lahav, P. B. Lilje, J. R. Primack,  and M. J. Rees, MNRAS 251, 128 (1991).
74. E. V. Linder, Phys. Rev. D 72, 043529 (2005)astro-ph/0507263 .
75. M. J. Chodorowski and P. Ciecielag, MNRAS 331, 133 (2002)astro-ph/0109291 .
76. P. Ciecielg, M. J. Chodorowski, M. Kiraga, M. A. Strauss, A. Kudlicki,  and F. R. Bouchet, MNRAS 339, 641 (2003)astro-ph/0010364 .
77. S. Pueblas and R. Scoccimarro, Phys. Rev. D 80, 043504 (2009)arXiv:0809.4606 .
78. N. A. Bahcall, M. Gramann,  and R. Cen, ApJ 436, 23 (1994)astro-ph/9410061 .
79. P. Coles and F. Lucchin, Cosmology: The Origin and Evolution of Cosmic Structure (Cosmology: The Origin and Evolution of Cosmic Structure, Second Edition, by Peter Coles, Francesco Lucchin, pp. 512. ISBN 0-471-48909-3. Wiley-VCH , July 2002., 2002) p. 512.
80. R. Watkins and H. A. Feldman, MNRAS 379, 343 (2007)astro-ph/0702751 .
81. A. Lewis, A. Challinor,  and A. Lasenby, ApJ 538, 473 (2000)astro-ph/9911177 .
82. R. E. Smith, J. A. Peacock, A. Jenkins, S. D. M. White, C. S. Frenk, F. R. Pearce, P. A. Thomas, G. Efstathiou,  and H. M. P. Couchman, MNRAS 341, 1311 (2003)astro-ph/0207664 .
83. C. Pichon and F. Bernardeau, A&A 343, 663 (1999), astro-ph/9902142 .
84. N. I. Libeskind, Y. Hoffman, M. Steinmetz, S. Gottlöber, A. Knebe,  and S. Hess, ApJ 766, L15 (2013)arXiv:1212.1454 .
85. N. I. Libeskind, Y. Hoffman,  and S. Gottlöber, MNRAS 441, 1974 (2014)arXiv:1310.5706 .
86. O. Hahn, T. Abel,  and R. Kaehler, MNRAS 434, 1171 (2013)arXiv:1210.6652 .
87. O. Hahn, R. E. Angulo,  and T. Abel, MNRAS 454, 3920 (2015)arXiv:1404.2280 .
88. C. Rampf and U. Frisch, MNRAS 471, 671 (2017)arXiv:1705.08456 .
89. R. C. Kraan-Korteweg, P. A. Woudt, V. Cayatte, A. P. Fairall, C. Balkowski,  and P. A. Henning, Nature 379, 519 (1996).
90. R. C. Kraan-Korteweg, M. E. Cluver, M. Bilicki, T. H. Jarrett, M. Colless, A. Elagali, H. Böhringer,  and G. Chon, MNRAS 466, L29 (2017)arXiv:1611.04615 .
91. R. C. Kraan-Korteweg and O. Lahav, A&A Rev. 10, 211 (2000)astro-ph/0005501 .
92. R. B. Tully and J. R. Fisher, A&A 54, 661 (1977).
93. S. Djorgovski and M. Davis, ApJ 313, 59 (1987).
94. H. M. Courtois, D. Pomarède, R. B. Tully, Y. Hoffman,  and D. Courtois, AJ 146, 69 (2013)arXiv:1306.0091 .
95. N. Kaiser, MNRAS 231, 149 (1988).
96. D. Lynden-Bell, D. Burstein, R. L. Davies, A. Dressler,  and S. M. Faber, in The Extragalactic Distance Scale, Astronomical Society of the Pacific Conference Series, Vol. 4, edited by S. van den Bergh and C. J. Pritchet (1988) pp. 307–316.
97. D. Lynden-Bell, S. M. Faber, D. Burstein, R. L. Davies, A. Dressler, R. J. Terlevich,  and G. Wegner, ApJ 326, 19 (1988b).
98. R. B. Tully and J. R. Fisher, Annales de Geophysique (1987).
99. R. B. Tully and J. R. Fisher, Catalog of Nearby Galaxies (1988) p. 224.
100. M. J. Hudson, MNRAS 265, 43 (1993).
101. R. B. Tully, E. J. Shaya, I. D. Karachentsev, H. M. Courtois, D. D. Kocevski, L. Rizzi,  and A. Peel, ApJ 676, 184-205 (2008)arXiv:0705.4139 .
102. A. Sandage, G. A. Tammann,  and E. Hardy, ApJ 172, 253 (1972).
103. D. Schlegel, M. Davis, F. Summers,  and J. A. Holtzman, ApJ 427, 527 (1994)astro-ph/9310040 .
104. I. D. Karachentsev, M. E. Sharina, D. I. Makarov, A. E. Dolphin, E. K. Grebel, D. Geisler, P. Guhathakurta, P. W. Hodge, V. E. Karachentseva, A. Sarajedini,  and P. Seitzer, A&A 389, 812 (2002)astro-ph/0204507 .
105. I. D. Karachentsev, D. I. Makarov, M. E. Sharina, A. E. Dolphin, E. K. Grebel, D. Geisler, P. Guhathakurta, P. W. Hodge, V. E. Karachentseva, A. Sarajedini,  and P. Seitzer, A&A 398, 479 (2003)astro-ph/0211011 .
106. R. B. Tully and E. J. Shaya, ApJ 281, 31 (1984).
107. G. A. Tammann and A. Sandage, ApJ 294, 81 (1985).
108. N. Y. Lu, E. E. Salpeter,  and G. L. Hoffman, ApJ 426, 473 (1994).
109. D. Gudehus, A&A 302, 21 (1995).
110. I. D. Karachentsev, R. B. Tully, P.-F. Wu, E. J. Shaya,  and A. E. Dolphin, ApJ 782, 4 (2014)arXiv:1312.6769 .
111. N. I. Libeskind, Y. Hoffman, R. B. Tully, H. M. Courtois, D. Pomarède, S. Gottlöber,  and M. Steinmetz, MNRAS 452, 1052 (2015)arXiv:1503.05915 .
112. M. T. Busha, R. H. Wechsler, P. S. Behroozi, B. F. Gerke, A. A. Klypin,  and J. R. Primack, ApJ 743, 117 (2011)arXiv:1011.6373 .
113. S. Phelps, A. Nusser,  and V. Desjacques, ApJ 775, 102 (2013)arXiv:1306.4013 [astro-ph.CO] .
114. M. Cautun, C. S. Frenk, R. van de Weygaert, W. A. Hellwing,  and B. J. T. Jones, MNRAS 445, 2049 (2014)arXiv:1405.7697 .
115. Q. Guo, A. P. Cooper, C. Frenk, J. Helly,  and W. A. Hellwing, MNRAS 454, 550 (2015)arXiv:1503.08508 .
116. A. Kogut, C. Lineweaver, G. F. Smoot, C. L. Bennett, A. Banday, N. W. Boggess, E. S. Cheng, G. de Amici, D. J. Fixsen, G. Hinshaw, P. D. Jackson, M. Janssen, P. Keegstra, K. Loewenstein, P. Lubin, J. C. Mather, L. Tenorio, R. Weiss, D. T. Wilkinson,  and E. L. Wright, ApJ 419, 1 (1993)astro-ph/9312056 .
117. I. D. Karachentsev, V. E. Karachentseva, O. V. Melnyk, A. A. Elyiv,  and D. I. Makarov, Astrophysical Bulletin 67, 353 (2012)arXiv:1210.6571 .
118. A. A. Elyiv, I. D. Karachentsev, V. E. Karachentseva, O. V. Melnyk,  and D. I. Makarov, Astrophysical Bulletin 68, 1 (2013)arXiv:1302.2369 .
119. R. B. Tully, H. Courtois, Y. Hoffman,  and D. Pomarède, Nature 513, 71 (2014)arXiv:1409.0880 .
120. S. Mei, J. P. Blakeslee, P. Côté, J. L. Tonry, M. J. West, L. Ferrarese, A. Jordán, E. W. Peng, A. Anthony,  and D. Merritt, ApJ 655, 144 (2007)astro-ph/0702510 .
121. Y. Hoffman and E. Ribak, ApJ 384, 448 (1992).
122. R. van de Weygaert and E. Bertschinger, MNRAS 281, 84 (1996)astro-ph/9507024 .
123. H. A. Feldman, R. Watkins,  and M. J. Hudson, MNRAS 407, 2328 (2010)arXiv:0911.5516 [astro-ph.CO] .
124. V. Sahni and Y. Shtanov, J. Cosmology Astropart. Phys 11, 014 (2003)astro-ph/0202346 .
125. A. Lue and G. D. Starkman, Phys. Rev. D 70, 101501 (2004)astro-ph/0408246 .
126. A. Vainshtein, Phys.Lett. B39, 393 (1972).
127. A. Barreira, B. Li, W. A. Hellwing, C. M. Baugh,  and S. Pascoli, J. Cosmology Astropart. Phys 10, 027 (2013)arXiv:1306.3219 .
128. K. Koyama and R. Maartens, J. Cosmology Astropart. Phys 1, 016 (2006)astro-ph/0511634 .
129. B. Li, G.-B. Zhao,  and K. Koyama, J. Cosmology Astropart. Phys 5, 023 (2013b)arXiv:1303.0008 [astro-ph.CO] .
130. B. Bose, K. Koyama, W. A. Hellwing,  and G.-B. Zhao, ArXiv e-prints  (2017), arXiv:1702.02348 .
131. A. Barreira, S. Bose,  and B. Li, J. Cosmology Astropart. Phys 12, 059 (2015)arXiv:1511.08200 .
132. W. A. Hellwing, K. Koyama, B. Bose,  and G.-B. Zhao, Phys. Rev. D 96, 023515 (2017b)arXiv:1703.03395 .
133. G. Cusin, V. Tansella,  and R. Durrer, Phys. Rev. D 95, 063527 (2017)arXiv:1612.00783 .
134. J. H. J. de Bruijne, Ap&SS 341, 31 (2012)arXiv:1201.3238 [astro-ph.IM] .
135. L. S. S. Telescope,  .
136. A. Nusser, E. Branchini,  and M. Davis, ApJ 755, 58 (2012)arXiv:1202.4138 [astro-ph.CO] .
137. J. Seiler and D. Parkinson, MNRAS 462, 75 (2016)arXiv:1607.01837 .
138. Y. Hoffman, A. Nusser, H. M. Courtois,  and R. B. Tully, MNRAS 461, 4176 (2016)arXiv:1605.02285 .
139. S. Gottloeber, Y. Hoffman,  and G. Yepes, ArXiv e-prints  (2010), arXiv:1005.2687 .
140. H. M. Courtois and R. B. Tully, Astronomische Nachrichten 333, 436 (2012)arXiv:1205.4627 [astro-ph.CO] .
141. S. Heß, F.-S. Kitaura,  and S. Gottlöber, MNRAS 435, 2065 (2013)arXiv:1304.6565 .
142. J. G. Sorce, H. M. Courtois, S. Gottlöber, Y. Hoffman,  and R. B. Tully, MNRAS 437, 3586 (2014)arXiv:1311.2253 .
143. F. Leclercq, J. Jasche,  and B. Wandelt, J. Cosmology Astropart. Phys 6, 015 (2015)arXiv:1502.02690 .
144. T. Sawala, C. S. Frenk, A. Fattahi, J. F. Navarro, R. G. Bower, R. A. Crain, C. Dalla Vecchia, M. Furlong, J. C. Helly, A. Jenkins, K. A. Oman, M. Schaller, J. Schaye, T. Theuns, J. Trayford,  and S. D. M. White, MNRAS 457, 1931 (2016)arXiv:1511.01098 .
145. J. G. Sorce, S. Gottlöber, G. Yepes, Y. Hoffman, H. M. Courtois, M. Steinmetz, R. B. Tully, D. Pomarède,  and E. Carlesi, MNRAS 455, 2078 (2016)arXiv:1510.04900 .
146. F. Leclercq, J. Jasche, G. Lavaux, B. Wandelt,  and W. Percival, J. Cosmology Astropart. Phys 6, 049 (2017)arXiv:1601.00093 .
147. H. Desmond, P. G. Ferreira, G. Lavaux,  and J. Jasche, MNRAS 474, 3152 (2018)arXiv:1705.02420 .
148. S. Tassev, M. Zaldarriaga,  and D. J. Eisenstein, J. Cosmology Astropart. Phys 6, 036 (2013)arXiv:1301.0322 [astro-ph.CO] .
149. C. Howlett, M. Manera,  and W. J. Percival, Astronomy and Computing 12, 109 (2015)arXiv:1506.03737 .
150. H. A. Winther, K. Koyama, M. Manera, B. S. Wright,  and G.-B. Zhao, J. Cosmology Astropart. Phys 8, 006 (2017)arXiv:1703.00879 .
151. E. Munari, P. Monaco, J. Koda, F.-S. Kitaura, E. Sefusatti,  and S. Borgani, J. Cosmology Astropart. Phys 7, 050 (2017)arXiv:1704.00920 .
152. For our purpose here, for a measure of the distribution spread (i.e. ) we use the half-width between the 16-th and 84-th percentiles.