Accurate Determination of halo velocity bias in simulations and its cosmological implications
Abstract
A long standing issue in peculiar velocity cosmology is whether the halo/galaxy velocity bias at large scale. The resolution of this important issue must resort to high precision cosmological simulations. However, this is hampered by another long standing “sampling artifact” problem in volume weighted velocity measurement. We circumvent this problem with a hybrid approach. We first measure statistics free of sampling artifact, then link them to volume weighted statistics in theory, and finally solve for the velocity bias. determined by our method is not only free of sampling artifact, but also free of cosmic variance. We apply this method to a CDM Nbody simulation of simulation particles and box size. For the first time we determine the halo velocity bias of various mass and redshift to  accuracy. Major findings are as follows. (1) at . The deviation from unity () increases with . Depending on halo mass and redshift, it may reach at and at . The discovered has statistically significant impact on structure growth rate measurement by spectroscopic redshift surveys, including DESI, Euclid and SKA. (2) Both the sign and the amplitude of depend on mass and redshift. These results disagree with the peak model prediction on protohalos in that has much weaker deviation from unity, varies with redshift, and can be bigger than unity. (3) Most of the mass and redshift dependences can be compressed into a single dependence on the halo density bias. Based on this finding, we provide an approximate twoparameter fitting formula.
Subject headings:
cosmology: observations: largescale structure of universe: dark matter: dark energy1. introduction
A crucial (but often ignored) assumption in cosmology based on large scale peculiar velocity is that the galaxy/halo velocity bias equals unity (). The usual argument is based on the weak equivalence principle. Since galaxies/halos are not the dominant sources of gravity over Mpc scale, they can be treated as test particles. Therefore their motions should faithfully follow dark matter and therefore at Mpc scale.
However, this argument overlooks the fact that halos/galaxies only reside at (local) density peaks. Along with the fact that the density gradient is tightly correlated with the velocity field, the seminal BBKS (Bardeen et al., 1986) paper predicted . This result was derived using onepoint Gaussian statistics. Desjacques & Sheth (2010) (hereafter DS10) extended the peak model to 2point statistics. They derived an elegant expression, . The prefactor depends on the halo mass , but not redshift. DS10 predicted significant and redshift independent deviation of from unity, even at scales Mpc. For example, at for protohalos (e.g. Elia et al. (2012)). Later theoretical and numerical works (Elia et al., 2012; Chan et al., 2012; Baldauf et al., 2014; Chan, 2015) investigated and verified the DS10 prediction. Nonetheless, these works are all for protohalos defined in the linear and Gaussian initial conditions, instead of real halos which host galaxies in observations. Due to stochastic relation between protohalos and halos, and complexities in halo velocity evolution (e.g. Colberg et al. (2000)), ambiguities exist to extrapolate these works to peculiar velocities of real halos/galaxies.
Therefore, despite of decades of effort, the velocity bias remains an unresolved issue. Even worse, it will become a limiting factor of peculiar velocity cosmology in the near future (e.g. Howlett et al. (2017)). In particular, stage IV dark energy surveys such as DESI, PFS, Euclid, SKA and WFIRST (e.g. Schlegel et al. (2011); Spergel et al. (2015); Abdalla et al. (2015); DESI Collaboration et al. (2016); Amendola et al. (2016)) aim to measure the structure growth rate to or higher accuracy, through redshift space distortion (RSD). However, what RSD actually measures is the galaxy peculiar velocity and therefore the combination . Systematic bias in the understanding of then induces a systematic error in ,
(1) 
is in principle dependent of scale . So the induced systematic error depends on the small scale cut (namely we only use the velocity information at ). Various cuts have been adopted in peculiar velocity cosmology. DESI Collaboration et al. (2016) adopts in DESI forecast. The eBOSS collaboration adopts in the quasar power spectrum analysis, and in the quasar correlation function analysis (GilMarín et al., 2018; Ruggeri et al., 2018; Zhao et al., 2018; Hou et al., 2018; Zarrouk et al., 2018). Ambitious stage V dark energy projects (Dodelson et al., 2016) can measure RSD to smaller scales, and will further significantly improve cosmological constraints. To match the survey capability, we have to theoretically understand at  level accuracy over the range .
Set ID  mass range  

37.7  7.1  1.36  
29.8  5.5  1.89  
23.7  3.5  2.68  
17.7  0.88  4.98  
110  2.7  54.4  0.81  
110  2.6  52.5  1.04  
110  2.5  46.9  1.48  
110  2.3  28.6  2.64  
0.51  0.70  52.9  0.70  
0.51  0.70  54.0  0.86  
0.51  0.69  52.4  1.15  
0.51  0.69  39.7  1.97  
0.51  0.68  21.8  3.15  
0.51  0.70  52.9  0.70  
17  2.4  51.5  0.80  
710  8.2  3.3  1.02  
37.7  7.1  1.36  
710  8.2  3.3  1.02  
1.24.0  2.1  34.1  1.02  
0.310.35  0.33  20.3  1.01 
Since halos are highly nonlinear objects, we have to resort to high precision cosmological simulations to accurately measure their velocity statistics. State of art simulations are already able to reliably simulate the formation and evolution of halos hosting target galaxies in cosmological surveys, and generate accurate phasespace distribution of halos.
Nevertheless, translating the accurately simulated halo phasespace distribution into the volume weighted velocity statistics^{1}^{1}1Unlike the density weighted velocity statistics, the volume weighted velocity does not depend on the galaxy density bias. So it is preferred in constraining cosmology with peculiar velocity. Furthermore, the volume weighted velocity statistics is preferred in many RSD models (e.g. Kaiser (1987); Scoccimarro (2004); Taruya et al. (2010); Zhang et al. (2013)). But there are exceptions such as the distribution function approach (Seljak & McDonald, 2011) and the streaming model (Peebles, 1980; White et al., 2015). is highly nontrivial, due to a long standing problem of “sampling artifact”. This problem exists for galaxies, halos and simulation particles. One way to demonstrate its existence is to randomly select a fraction of these objects and then measure the velocity power spectrum of this subsample. The measured velocity power spectrum should be independent of . However, both theoretical and numerical investigations show significant dependence, and therefore the existence of sampling artifact (Zhang et al., 2015; Zheng et al., 2015b). This sampling artifact problem arises from the fact that we only know the velocities where there are objects (halos, galaxies, simulation particles). Their distributions are not only inhomogeneous, but also correlated with their velocity fields. So the sampling of their velocity fields is biased, leading to biased measurement of volume weighted velocity statistics (e.g. Bernardeau & van de Weygaert (1996); Bernardeau et al. (1997); Schaap & van de Weygaert (2000); Pueblas & Scoccimarro (2009); Zheng et al. (2013); Zhang et al. (2015); Zheng et al. (2015b); Jennings et al. (2015)).
Set ID  

For the dark matter velocity statistics, this is essentially an issue of mass resolution. If on the average there are simulation particles in a volume , the velocity field above the scale is well sampled and the sampling artifact is negligible. Increasing the mass resolution increases the number density of simulation particles and pushes the reliable measurement to smaller scales. Unfortunately, the sampling artifact in the halo velocity statistics is much worse and can not be reduced by increasing the mass resolution. First, halos are much more sparse than simulation particles, causing much severer sampling artifact. Second, increasing the mass resolution can not alleviate the sampling artifact problem, since the halo number density is fixed at given redshift and mass. This sampling artifact prohibits accurate determination of halo velocity bias in simulations, with existing methods.
We have tried various approaches to measure the volume weighted velocity statistics in simulations. We have designed new velocity assignment methods, including the NP method (Zheng et al., 2013) and the Kriging method (Yu et al., 2015, 2017). We have built theoretical model of the sampling artifact (Zhang et al., 2015), verified it in simulated dark matter velocity field (Zheng et al., 2015b) and then applied it to correct the sampling artifact in the halo velocity field (Zheng et al., 2015a). Despite these efforts, we have not yet succeeded in measuring with accuracy at . The measurement at larger is even more challenging.
The current paper presents an exact approach to determine from simulations. It circumvents the problem of sampling artifact by a hybrid method. For the first time, we determine to  at and , for various halo mass bins. In §2, we present our method of determining . We leave further technical details in the appendix. In §3 we describe the simulation used for the velocity bias measurement. In §4 we present the determined for various halo mass and redshift. In §5 we discuss its impact on cosmological surveys. We discuss and conclude in §6.
2. The method
The method is hybrid in the sense that it is composed of two steps, direct measurements and theoretical interpretation.

Direct measurements are for three quantities, all with negligible sampling artifact. (1) One is the halo momentum . It is the density weighted velocity, free of sampling artifact. (2) One is the matter velocity . In principle it contains sampling artifact. Fortunately, it has been suppressed to a negligible level in our simulation. Our simulation has on the average simulation particles per volume. Therefore the dark matter velocity field at scales of interest () is well sampled and is essentially free of sampling artifact (Zheng et al., 2013). (3) One is the halo number overdensity , for which the sampling artifact is irrelevant.

These direct measurements are then exactly linked to the volume weighted statistics in theory, where the only unknown parameter is . We then solve for . The determination of is then free of sampling artifact.
Step 1. We first measure the correlation function
(2) 
and its Fourier counterpart . These measurements are then linked to the following correlation functions
(3) 
and their power spectra (Fourier components)
(4) 
Here denote the volume/ensemble averaging.
Step 2. We then solve Eq. 4 for the velocity bias , defined through
(5) 
A key point to pay attention is that is the only unknown quantity in Eq. 4.^{2}^{2}2In Fourier space, we can decompose , where (). Therefore does not contribute to . Furthermore it does not contribute to . The direction of is uncorrelated with due to the statistical isotropy of the Universe. It is also uncorrelated with , by definition. Averaging over its direction, we have . Therefore does not contribute to the right hand side of Eq. 4. Since we know and from the same simulation, is all we need to fix the right hand side of Eq. 4. Then comparing the left and right hand sides drawn from the same simulation, we determine . For this reason, the determined is essentially free of cosmic variance. In the appendix, we present the maximum likelihood approach to solve Eq. 4 for .
3. The simulation
The simulation we analyze (J6620) adopts the standard CDM cosmology, with , , , , and . It has box size , particle number and the mass resolution . J6620 is run with a particleparticleparticlemesh () code, detailed in Jing et al. (2007). The halo catalogue is constructed by the FriendsofFriends (FOF) method. The linking length is . In the catalogue gravitationally unbound “halos” have been excluded. The halo center is defined as the mass weighted center and the halo velocity is defined as the velocity averaged over all member particles. We adopt various mass and redshift bins to calculate the mass and redshift dependence of velocity bias. Table 1 shows details of these bins.
The number density and momentum fields of both halos and dark matter are measured using the NGP method. Namely . The sum () is over all particles in the given cell. is the mean number of halos in each cell. We adopt grid points. The grid cell size is . Each cell has 216 simulation particles on the average. Therefore we have excellent sampling of the dark matter velocity field above such scale. This allows us to robustly measure the dark matter velocity field through , with negligible sampling artifact. The aliasing effect is also negligible, since we are interested in the scales of , much smaller than the Neyquist wavenumber (Jing, 2005).
Fig. 1 shows the directly measured for and halos at . For comparison, we also show . The three are almost identical to each other until . These terms have two contributions. The contribution from dominates at . The contribution from becomes significant at and becomes dominant at . These results already imply at . The exact determination of requires to solve Eq. 4 using the method described in the appendix.
4. The velocity bias
We obtain the bestfit value and the associate error of over the ranges of , , , , , , , . As explained early, the determined is essentially free of cosmic variance, since it is obtained by comparing the halo and dark matter fields in the same simulation box. The only source of noise is shot noise in the halo distribution. The large number of halos () then enables us to determine with  statistical error. Such accuracy also enables us to detect any significant deviation of from unity.
Fig. 2, 3 & 4 show the redshift dependence of velocity bias for three mass bins (halo set A1, A2, A3 respectively). Fig. 5 shows the mass dependence of velocity bias at (set B). Furthermore, table 2 lists the velocity bias at selected ranges of .
We detect statistically significant deviation of from unity at . This invalidates the assumption of commonly adopted in peculiar velocity cosmology. It will significantly impact RSD cosmology of stage IV dark energy projects. The deviation shows rich behavior in , halo mass and redshift . Nonetheless, it shows significant difference to the peak model prediction and poses new question to the halo peculiar velocity theory. Major findings are as follows.
4.1. The dependence
can have either negative or positive sign. This challenges the peak mode, which predits a negative sign. The sign of does not vary with . Furthermore, increases with , and roughly speaking . (1) At , the deviation is very weak. Over all the halo mass and redshift investigated, the deviation is or less and it is statistically insignificant. is orders of magnitude weaker than the peak model prediction on protohalos. It means that we can not simply extrapolate the predictions on protohalo velocity to real halo velocity. (2) At , may show statistically significant deviation from unity. Depending on the type of halos, the deviation may reach . As discussed later, despite its weakness, it will become a significant source of systematic error for DESI RSD cosmology. (3) At , some halo samples show deviation from unity. One task of RSD cosmology is to extract cosmological information deep into the nonlinear regime of (e.g., the cosmic vision dark energy report (Dodelson et al., 2016)). The existence of significant deviation of from unity at this regime is a challenge to this task.
4.2. The mass and redshift dependence
increases with increasing redshift (Fig. 2, 3 & 4) and decreasing mass (Fig. 5). As a consequence, the sign of depends on halo mass and redshift. For example, is always negative for halos of at all redshifts (Fig. 2). But it changes from positive at to negative at for less massive halos (Fig. 3 & 4). Another consequence is that has strong dependence on the halo mass and redshift. For example, for halos at and , . But for halos of , at . Fig. 5 compares of various mass at . Now the biggest deviation from unity happens for the least massive halos.
To translate the above results into impact on peculiar velocity cosmology, we need specifications of galaxy surveys, since different surveys probe different galaxies in different halos and different redshifts. Here we just present a qualitative description on halos that may be probed by various surveys. Later will quantify the impact of velocity bias to some of these surveys in §5. (1) halos at may be probed by LRGs in DESI (e.g. Guo et al. (2015); DESI Collaboration et al. (2016)). Galaxies in the TAIPAN redshift and peculiar velocity survey are also expected to reside in these halos, but at (Howlett et al., 2017). (2) For smaller halos (), 21cm surveys may be capable of detecting them. SKA are capable of detecting billions of 21cm emitting galaxies residing in these halos, given its sensitivity in HI mass (Abdalla & Rawlings, 2005; Yang & Zhang, 2011) and the observationally constrained HI masshalo mass relation (e.g. Guo et al. (2017)). SKA is also able to indirectly detect them through the intensity mapping. (3) For halos of , a large fraction of ELGs in surveys such as DESI and PFS may reside in these halos (Favole et al., 2016; GonzalezPerez et al., 2017). DESI can probe them at while PFS can probe them to higher redshift (). 21cm intensity mapping surveys such as CHIME (Bandura et al., 2014) and Tianlai (Xu et al., 2015) are also sensitive to these halos, although they may not be able to detect individual galaxies.
4.3. The dependence on the halo density bias
An interesting observation is that the mass and redshift dependence of may be absorbed into a single dependence on the halo density bias . This can be seen by first checking in the table 1 and then comparing of halos with similar . More explicitly, halos of set C (Table 1 ) at different redshifts are chosen to have identical density bias (). Table. 2 shows that they have roughly the same velocity bias. This motivates us to propose the following fitting formula,
(6) 
Here . This is basically the Taylor expansion of around . The isotropy of the velocity bias () requires that terms of odd power in such as and () vanish in the Taylor expansion. Therefore the leading order term is . We find that and (Fig. 6). Small error in shows that the dependence on is statistically significant. We caution that this fitting formular is only approximate, since it ignores dependence beyond and ignores higher order dependence (e.g. ). Nevertheless, it is sufficient to describe the over dependence of the velocity bias on halo property and scale (Fig. 6).
5. Implications for peculiar velocity surveys
We discuss two implications of velocity bias on cosmology. The first is that it may bias the structure growth rate measurement in spectroscopic redshift surveys. The second is that it may open a window of testing the equivalence principle at cosmological scales.
5.1. Impact on structure growth rate constraint
A major task of cosmological surveys is to constrain the structure growth rate through peculiar velocity. The velocity bias, if ignored or modelled inappropriately, will become a source of systematic error. Whether it is of statistical significance depends on surveys and galaxy types. The low redshift TAIPAN survey aims to measure the peculiar velocities of galaxies. It will constrain at with accuracy, using information at (Howlett et al., 2017). The target galaxies (with ) have at (Fig. 2 and Table 2). Therefore the usual assumption of results into negligible systematic error, and can be adopted safely.
On the other hand, the spectroscopic redshift survey DESI can measure to level over a number of redshift bins, and resulting into an overall statistical error (DESI Collaboration et al., 2016). Fig. 7 plots the predicted of various galaxies in DESI. We estimate using the fitting formula of Eq. 6. The density biases of various galaxy types are adopted as , and (DESI Collaboration et al., 2016). Here is the linear density growth rate and it is normalized as unity at . For DESI, we are no longer able to approximate , otherwise systematic error at level can be induced.
For PFS ELGs, the predicted at is similar to that of DESI ELGs. So for clarity we neglect them in Fig. 7. We only show the results at , where we adopt (the PFS SSP proposal). If only using information at ,the systematic error in caused by ignoring is by less than . Due to a factor of smaller sky coverage than DESI, the overall statistical error of constrained by PFS RSD to is expect to be . Therefore the impact of velocity bias on PFS RSD is subdominant and we can simply approximate . However, due to higher number density of PFS galaxies, it can measure RSD to smaller scales and has the potential to further reduce the statistical error in . Fig. 7 shows that, if we push to , we may no longer adopt .
The proposed SKA HI survey has the capability of detecting 21cm emitting galaxies to over deg (Abdalla et al., 2015). The statistical error in is for each redshift bins over . If assuming , the induced systematic error will overwhelm the statistical error. This will also be true for the proposed stage V survey of measuring a billion spectra of LSST galaxies (Dodelson et al., 2016). The situation for Euclid may fall between DESI and SKA.
5.2. A cosmological test of the equivalence principle
Velocity bias also provides a new test of CDM cosmology. Observationally we are not able to measure the velocity bias (with respect to dark matter) directly. However, we are able to measure the relative velocity bias between two tracers of the large scale structure (LSS). Furthermore, if the two tracers overlap in space, the measured ratio will be free of cosmic variance (McDonald & Seljak, 2009). Our result predicts that in CDM, the velocity ratio of two tracers is
(7) 
This weak deviation from unity is a genuine consequence of the equivalence principle (EP). Therefore if or large deviation at is detected, it may be a smoking gun of EP violation and therefore modifications of general relativity (e.g. Hui et al. (2009)). We will further investigate this issue in future works.
6. Discussions and conclusions
We invent a novel method to determine the volume weighted halo velocity bias . This method is free of the long standing sampling artifact problem, which has hindered accurate determination of velocity bias. We apply it to a particle simulation and measure to with better than accuracy. Our findings confront both the standard assumption in peculiar velocity data analysis, and the peak model prediction. (1) There exists statistically significant deviation of from unity at . Depending on halo mass, redshift, may reach at and at . If ignored, this velocity bias will become a significant source of systematic error in RSD cosmology of DESI. Its impacts on SKA HI galaxy survey and Euclid are even stronger. (2) However, is a factor of smaller than the prediction of peak model. Furthermore, its mass and redshift dependence do not agree with the peak model prediction. varies with redshift, while the peak model predicts the opposite. of less massive halos can be bigger than unity, while the peak model always predicts . The peak model is based on protohalo statistics. Therefore we have to consider the mismatch between protohalos and real halos, and the displacement of halos from the corresponding protohalos, to improve the theoretical understanding of velocity bias. We also expect that the environment that halos reside (e.g. filaments or clusters) may play a role in the halo velocity bias. For example, the infall velocity within filaments may be responsible or partly responsible to the behavior of less massive halos.
There are many issues for further investigations. For example, since the velocity bias depends on the density bias, would it also depend on the halo formation time? Or more general, besides the density bias, what could affect the velocity bias? How does it depend on parameters within the standard cosmology? How does it behave in modified gravity cosmology? Also, to robustly predict its impact on RSD cosmology, we need to generate realistic mocks for target surveys and measure the velocity bias of LRGs, ELGs, 21cm emitting galaxies, etc.
Acknowledgments
This work was supported by the National Science Foundation of China (11621303, 11433001, 11653003, 11320101002, 11533006), National Basic Research Program of China (2015CB85701, 2015CB857003).
Appendix A The algorithm to solve for the scale dependent velocity bias
In §2 we prove that is the only unknown quantity in Eq. 4, where all other quantities are provided by the same simulation. This allows us to determine uniquely. The only complexity in determining is the nonlocal dependence of the right hand side of Eq. 4 on . It is caused by , in which also contributes. Here we present the maximum likelihood solution to .
We bin the unknown into a number of bins, . if , and zero otherwise. is the averaged value of in the range . The power spectrum . Here is in which we replace with . The calculation is done with several FFTs. First we obtain from the simulated . We then inverse Fourier transform and denote it as . Finally we Fourier transform , multiply it by , and obtain . The estimated/modelled power spectrum is then
(A1) 
Here both and are measured from the same simulation, and the only set of unknown parameters are . Both sides are drawn from the same simulation, therefore the determined will be essentially free of cosmic variance.
The only relevant statistical error in determining then arises from shot noise in the halo distribution. This allows us to write down the likelihood function straightforwardly.
(A2) 
is the r.m.s. fluctuation of the data caused by the finite number of halos. We estimate it by splitting halos in a given mass bin into nonoverlapping subsamples by randomly selecting among these halos. We measure the dispersion between these subsamples, divide it by and obtain . Since cosmic variance is irrelevant, the errors are uncorrelated over different . The bestfit is given by the following equation,
(A3) 
The error covariance matrix is given by
(A4) 
References
 Abdalla & Rawlings (2005) Abdalla, F. B., & Rawlings, S. 2005, MNRAS, 360, 27
 Abdalla et al. (2015) Abdalla, F. B., Bull, P., Camera, S., et al. 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), 17
 Amendola et al. (2016) Amendola, L., Appleby, S., Avgoustidis, A., et al. 2016, ArXiv eprints, arXiv:1606.00180
 Baldauf et al. (2014) Baldauf, T., Desjacques, V., & Seljak, U. 2014, ArXiv eprints, arXiv:1405.5885
 Bandura et al. (2014) Bandura, K., Addison, G. E., Amiri, M., et al. 2014, in Proc. SPIE, Vol. 9145, Groundbased and Airborne Telescopes V, 914522
 Bardeen et al. (1986) Bardeen, J. M., Bond, J. R., Kaiser, N., & Szalay, A. S. 1986, ApJ, 304, 15
 Bernardeau & van de Weygaert (1996) Bernardeau, F., & van de Weygaert, R. 1996, MNRAS, 279, 693
 Bernardeau et al. (1997) Bernardeau, F., van de Weygaert, R., Hivon, E., & Bouchet, F. R. 1997, MNRAS, 290, 566
 Chan (2015) Chan, K. C. 2015, Phys. Rev. D, 92, 123525
 Chan et al. (2012) Chan, K. C., Scoccimarro, R., & Sheth, R. K. 2012, Phys. Rev. D, 85, 083509
 Colberg et al. (2000) Colberg, J. M., White, S. D. M., MacFarland, T. J., et al. 2000, MNRAS, 313, 229
 DESI Collaboration et al. (2016) DESI Collaboration, Aghamousa, A., Aguilar, J., et al. 2016, ArXiv eprints, arXiv:1611.00036
 Desjacques & Sheth (2010) Desjacques, V., & Sheth, R. K. 2010, Phys. Rev. D, 81, 023526
 Dodelson et al. (2016) Dodelson, S., Heitmann, K., Hirata, C., et al. 2016, ArXiv eprints, arXiv:1604.07626
 Elia et al. (2012) Elia, A., Ludlow, A. D., & Porciani, C. 2012, MNRAS, 421, 3472
 Favole et al. (2016) Favole, G., Comparat, J., Prada, F., et al. 2016, MNRAS, 461, 3421
 GilMarín et al. (2018) GilMarín, H., Guy, J., Zarrouk, P., et al. 2018, ArXiv eprints, arXiv:1801.02689
 GonzalezPerez et al. (2017) GonzalezPerez, V., Comparat, J., Norberg, P., et al. 2017, ArXiv eprints, arXiv:1708.07628
 Guo et al. (2017) Guo, H., Li, C., Zheng, Z., et al. 2017, ArXiv eprints, arXiv:1707.01999
 Guo et al. (2015) Guo, H., Zheng, Z., Zehavi, I., et al. 2015, MNRAS, 446, 578
 Hou et al. (2018) Hou, J., Sánchez, A. G., Scoccimarro, R., et al. 2018, ArXiv eprints, arXiv:1801.02656
 Howlett et al. (2017) Howlett, C., StaveleySmith, L., & Blake, C. 2017, MNRAS, 464, 2517
 Hui et al. (2009) Hui, L., Nicolis, A., & Stubbs, C. W. 2009, Phys. Rev. D, 80, 104002
 Jennings et al. (2015) Jennings, E., Baugh, C. M., & Hatt, D. 2015, MNRAS, 446, 793
 Jing (2005) Jing, Y. P. 2005, ApJ, 620, 559
 Jing et al. (2007) Jing, Y. P., Suto, Y., & Mo, H. J. 2007, ApJ, 657, 664
 Kaiser (1987) Kaiser, N. 1987, MNRAS, 227, 1
 McDonald & Seljak (2009) McDonald, P., & Seljak, U. 2009, J. Cosmology Astropart. Phys, 10, 7
 Peebles (1980) Peebles, P. J. E. 1980, The largescale structure of the universe, ed. Peebles, P. J. E.
 Pueblas & Scoccimarro (2009) Pueblas, S., & Scoccimarro, R. 2009, Phys. Rev. D, 80, 043504
 Ruggeri et al. (2018) Ruggeri, R., Percival, W. J., Gil Marin, H., et al. 2018, ArXiv eprints, arXiv:1801.02891
 Schaap & van de Weygaert (2000) Schaap, W. E., & van de Weygaert, R. 2000, A&A, 363, L29
 Schlegel et al. (2011) Schlegel, D., Abdalla, F., Abraham, T., et al. 2011, ArXiv eprints, arXiv:1106.1706
 Scoccimarro (2004) Scoccimarro, R. 2004, Phys. Rev. D, 70, 083007
 Seljak & McDonald (2011) Seljak, U., & McDonald, P. 2011, J. Cosmology Astropart. Phys, 11, 39
 Spergel et al. (2015) Spergel, D., Gehrels, N., Baltay, C., et al. 2015, ArXiv eprints, arXiv:1503.03757
 Taruya et al. (2010) Taruya, A., Nishimichi, T., & Saito, S. 2010, Phys. Rev. D, 82, 063522
 White et al. (2015) White, M., Reid, B., Chuang, C.H., et al. 2015, MNRAS, 447, 234
 Xu et al. (2015) Xu, Y., Wang, X., & Chen, X. 2015, ApJ, 798, 40
 Yang & Zhang (2011) Yang, X., & Zhang, P. 2011, MNRAS, 415, 3485
 Yu et al. (2015) Yu, Y., Zhang, J., Jing, Y., & Zhang, P. 2015, Phys. Rev. D, 92, 083527
 Yu et al. (2017) —. 2017, Phys. Rev. D, 95, 043536
 Zarrouk et al. (2018) Zarrouk, P., Burtin, E., GilMarin, H., et al. 2018, ArXiv eprints, arXiv:1801.03062
 Zhang et al. (2013) Zhang, P., Pan, J., & Zheng, Y. 2013, Phys. Rev. D, 87, 063526
 Zhang et al. (2015) Zhang, P., Zheng, Y., & Jing, Y. 2015, Phys. Rev. D, 91, 043522
 Zhao et al. (2018) Zhao, G.B., Wang, Y., Saito, S., et al. 2018, ArXiv eprints, arXiv:1801.03043
 Zheng et al. (2015a) Zheng, Y., Zhang, P., & Jing, Y. 2015a, Phys. Rev. D, 91, 123512
 Zheng et al. (2015b) —. 2015b, Phys. Rev. D, 91, 043523
 Zheng et al. (2013) Zheng, Y., Zhang, P., Jing, Y., Lin, W., & Pan, J. 2013, Phys. Rev. D, 88, 103510