An accurate linear model for redshift space distortions in the voidgalaxy correlation function
Abstract
Redshift space distortions within voids provide a unique method to test for environmental dependence of the growth rate of structures in low density regions, where effects of modified gravity theories might be important. We derive a linear theory model for the redshift space voidgalaxy correlation that is valid at all pair separations, including deep within the void, and use this to obtain expressions for the monopole and quadrupole contributions. Our derivation highlights terms that have previously been neglected but are important within the void interior. As a result our model differs from previous works and predicts new physical effects, including a change in the sign of the quadrupole term within the void radius. We show how the model can be generalised to include a velocity dispersion. We compare our model predictions to measurements of the correlation function using mock void and galaxy catalogues modelled after the BOSS CMASS galaxy sample using the Big MultiDark body simulation, and show that the linear model with dispersion provides an excellent fit to the data at all scales, Mpc. While the RSD model matches simulations, the linear bias approximation does not hold within voids, and care is needed in fitting for the growth rate. Assuming a constant velocity dispersion, fits to the quadrupole recover the growth rate to a precision of using the simulation volume of at .
keywords:
gravitation – largescale structure of Universe –cosmology: observations – methods: data analysis – methods: analytical1 Introduction
Galaxy redshift surveys provide a map of the largescale structure of the Universe containing anisotropic distortions of the clustering caused by gravitationallyinduced peculiar velocities that contribute to the galaxy redshifts, as first predicted by Kaiser (1987). Measurement of these redshift space distortions (RSD) (e.g. Peacock et al., 2001; Guzzo et al., 2008; Beutler et al., 2012; Samushia et al., 2012; Reid et al., 2012; Howlett et al., 2015; Beutler et al., 2017) can be used to determine the growth rate of structure, which can provide strong tests of General Relativity (GR). The theory of RSD in the galaxy clustering is however complicated by significant nonlinear contributions which are important even at quite large pair separation scales, requiring sophisticated modelling (e.g. Scoccimarro, 2004; Matsubara, 2008; Taruya et al., 2010; Reid & White, 2011; Jennings et al., 2011).
The RSD modelling could in principle be simplified in underdense regions where the dynamics is closer to linear, such as voids. Cosmic voids are large underdensities in the galaxy distribution which can be used to trace stationary points of the gravitational potential (Nadathur et al., 2017), where velocities are dominated by coherent bulk flows. Voids have been much studied recently in other contexts, including their action as weak gravitational lenses (Krause et al., 2013; Melchior et al., 2014; Clampitt & Jain, 2015; Sánchez et al., 2016), the secondary CMB anisotropies they generate through the integrated SachsWolfe effect (e.g. Granett et al., 2008; Hotchkiss et al., 2015; Nadathur & Crittenden, 2016; Cai et al., 2017; Kovács et al., 2017) and the thermal SunyaevZeldovich effect (Alonso et al., 2017), and the Baryon Acoustic Oscillation (BAO) peak in their clustering (Kitaura et al., 2016). Velocity dynamics around voids have been studied using RSD in the voidgalaxy crosscorrelation (Paz et al., 2013; Hamaus et al., 2015; Cai et al., 2016; Hamaus et al., 2017; Achitouv et al., 2017) and that of the void positions themselves (Chuang et al., 2017).
In addition to possible simplifications in the linear modelling, an important advantage of studying RSD around voids is that it presents an opportunity to study the growth of density perturbations specifically in lowdensity regions. Several popular alternatives to GR, such as chameleon gravity (Hu & Sawicki, 2007), DvaliGabadadzePorrati (DGP) (Dvali et al., 2000), and Galileon (Nicolis et al., 2009; Deffayet et al., 2009) gravity models, among others, make use of screening mechanisms in order to suppress fifth force effects in high density regions such as the Solar System; theories with such screening effects therefore predict environmentdependent differences in the growth rate , which could be probed by the RSD effects within voids.
In this paper we study the voidgalaxy crosscorrelation function in redshift space. Cai et al. (2016) have previously studied the same problem and provided a linear model for , which was subsequently also used by Hamaus et al. (2017). However, this model correctly described simulation results only for voidgalaxy pair separations greater than the void size, , i.e. outside the lowdensity region of greatest interest. By extending the model using a form of the quasilinear streaming model, Cai et al. (2016) were able to extend the region of validity slightly to , though it was still not correct in the void centres. In fact in the void centre region the model can lead to unphysical predictions of .
We revisit the derivation of the voidgalaxy RSD model from first principles and identify linearorder terms that were neglected in the expression obtained by Cai et al. (2016) but are important in the void centre regions. These terms do not have counterparts in the linear Kaiser model for RSD in the galaxy correlation, and arise because of the restriction to void regions. They have important physical effects, in particular causing a change in sign of the quadrupole term within void regions. These terms become unimportant outside the void radius, where our expression matches that of Cai et al. (2016).
We show how the model can be extended to include a velocity dispersion without changing the linear nature of the theory. This is not the same as the form of the Gaussian streaming model used in several other studies of the voidgalaxy correlation (e.g Paz et al., 2013; Hamaus et al., 2015; Achitouv et al., 2017; Achitouv, 2017). We show why it is not correct to use the standard streaming model result for the galaxy correlation in the voidgalaxy case.
We then compare our theoretical results to data from mock void and galaxy catalogues from a large body simulation and show that our linear dispersion model provides an excellent fit to the data at all scales and performs significantly better than alternatives. We highlight the fact that the approximation of a linear galaxy bias does not hold within voids. We discuss the consequences of these results for obtaining an unbiased estimate of the growth rate within void environments.
The layout of the paper is as follows. Section 2 describes the simulation data we use and the creation of galaxy and void mocks. In Section 3 we derive the linear theory model for the voidgalaxy correlation and its multipoles, and discuss differences with previous works. In Section 4 we compare theoretical predictions to simulation data and in Section 5 we discuss strategies for measurement of the growth rate based on these results. We sum up and draw conclusions in Section 6.
2 Data
2.1 Simulation and galaxy mocks
We use data from the redshift snapshot from the Big MultiDark (BigMD) body simulation (Klypin et al., 2016) from the MultiDark simulation project (Prada et al., 2012). This simulation follows the evolution of particles in a box of side Mpc using the GADGET2 (Springel, 2005) and Adaptive Refinement Tree (Kravtsov et al., 1997; Gottloeber & Klypin, 2008) codes, with cosmological parameters , , , , and . Initial conditions for the simulation were set using the Zeldovich approximation at starting redshift .
A halo catalogue was created for the given snapshot using the Bound Density Maximum algorithm (Klypin & Holtzman, 1997; Riebe et al., 2013). We populated these halos with mock galaxies using the Halo Occupation Distribution (HOD) model of Zheng et al. (2007), assigning central and satellite galaxies to halos according to a distribution based on the halo mass. Details of the algorithm and HOD model parameters used are described more fully in Nadathur et al. (2017), who used the same mock catalogue: these parameters were taken from Manera et al. (2013) and are designed to approximately reproduce the clustering and mean number density for galaxies in the Baryon Oscillation Spectroscopic Survey (BOSS) CMASS galaxy sample.
To measure dark matter (DM) densities in the simulation, we used a cloudincell interpolation scheme to determine the DM density on a grid using the full particle output of the simulation. We used this grid density field to determine the dark matter power spectrum . We then measured the galaxy power spectrum for the mocks; by fitting for the ratio at large scales, , we determine the linear bias value for the galaxy mocks, .
Using the real space galaxy positions and velocities , we determine their redshift space positions in the planeparallel approximation assuming the line of sight direction to be along the axis of the simulation box,
(1) 
2.2 Void finding
We identify voids in the real space galaxy mocks using the ZOBOV watershed voidfinding algorithm (Neyrinck, 2008). The ZOBOV algorithm uses a Voronoi tessellation field estimator (VTFE) technique to reconstruct the galaxy density field from the discrete distribution, and then identifies local minima in this field and the watershed basins around them, to form a nonoverlapping set of voids corresponding to local density depressions. As in previous works (Nadathur, 2016; Nadathur et al., 2017), we define each individual density basin as a distinct void, without any additional merging of neighbouring regions. A fuller description of the algorithm and void properties can be found in Nadathur (2016) and Nadathur et al. (2017).
We determine the centre of each void to be the centre of the largest sphere completely empty of galaxies that can be inscribed within the void (Nadathur & Hotchkiss, 2015a; Nadathur et al., 2017). Although voids are of arbitrary shape and are in general far from spherically symmetric, it is convenient to define an effective spherical radius, , where is the total volume of the void.
It is important to stress that we apply the voidfinding algorithm to the real space galaxy mocks and not to the shifted version in redshift space. As we discuss in the next section, this is crucial because none of the theoretical models for RSD in the voidgalaxy crosscorrelation discussed in this paper or elsewhere in the literature are applicable unless the real space void positions are known. In a companion publication (Nadathur, Carter & Percival, in prep.), we show how this practical difficulty can be overcome when using survey data where the real space galaxy positions are not known.
Approximately 33000 voids are identified in the simulation box. As the voidfinding algorithm is spacefilling, these voids cover the entire box volume, and undoubtedly include some spurious identifications that do not correspond to genuine matter underdensities. We also do not expect a linear model of coherent velocity outflow from a void to successfully describe the RSD around all voids. In particular it is not expected to work where the local environment of structures outside the void is important in determining the velocity field. This is expected to be the case for small voids, and indeed Cai et al. (2016) find that RSD models for the voidgalaxy correlation do not work for small voids. However, the distinction between ‘large’ and ‘small’ voids is somewhat ambiguous, and the numerical value of the cut on void size depends both on the bias and number density of the galaxies in question (Nadathur & Hotchkiss, 2015a, b) as well as on the particular features of the voidfinder.
We therefore restrict our void sample to the half of all voids with effective radius greater than the median radius. This is an easily reproducible criterion. For the mocks and voids used in this work, this means selecting Mpc, which leaves voids. Throughout the rest of this paper, we exclusively use this sample of voids. The mean effective void radius for this sample is Mpc. In some figures in later sections we show distances from the void centre in units of this mean radius, for context. However in our analysis we do not rescale distances in units of individual void sizes. Such a rescaling requires a strong implicit assumption of selfsimilarity of void profiles dependent only on void size (Nadathur et al., 2015), which is not justified (Nadathur & Hotchkiss, 2015a, b; Nadathur et al., 2017). In addition, such a rescaling would effectively weight galaxy counts in the same volume differently depending on the assigned void effective size, which complicates the error determination.
Note that restricting our sample to the largest of voids need not necessarily be the optimal choice to ensure validity of linear theory. Nadathur et al. (2017) show that environmental effects around voids are more strongly correlated with a combination of void size and density than with void size alone, so this may provide a better selection criterion. However, the median size cut implemented here has the advantage of simplicity, and as we show later, is sufficient that a purely linear RSD model already provides an excellent fit to the data.
3 Theory
3.1 The voidgalaxy crosscorrelation in redshift space
Let denote the comoving location of a void centre, and the location of a galaxy in its vicinity. The realspace separation vector for the voidgalaxy pair, , is transformed to in redshift space. Assuming that the total number of voidgalaxy pairs is unchanged by the shift to redshiftspace, we require that
(2) 
where denotes the voidgalaxy crosscorrelation (in what follows we will suppress the subscript where there is no risk of confusion), and the supersripts and denote redshift and real space respectively. may also be viewed as the galaxy number density profile around the void, and is sometimes denoted in this context.
If we further assume that the void centre is stationary, so that only the galaxy peculiar velocity is relevant to the RSD pattern, in the distantobserver approximation we may write
(3) 
where is the scale factor and is the Hubble rate.
Note that the assumptions that void centres remain stationary and that their numbers are conserved between real and redshift space are both known to fail. That is, if voids are identified using the redshift space galaxy field, different numbers and centres of voids are obtained than if they are identified using the same algorithm applied to real space galaxies (see also the discussion in Chuang et al., 2017). In a separate work, we will show that these effects are very important, producing strong contributions to the observed RSD pattern even within the scale of the mean void radius. Given the complex and nonlinear operation of voidfinding algorithms, accounting for these effects in the modelling is challenging. In the present work we restrict ourselves to developing the theory for the more manageable case where both assumptions are valid. To ensure this, we identify voids using the real space galaxy positions in the simulation, and consider the RSD pattern about the centres thus defined. In a companion paper (Nadathur, Carter & Percival, in prep.) we discuss practical measures for subtracting this additional distortion from data so that a fair comparison with theory is possible.
From Eq. 3, the lineofsight component of the separation vector is
(4) 
Symmetry arguments show that the average coherent velocity flow must be directed along the radial direction, . The determinant of the Jacobian for the coordinate transformation is therefore
(5) 
where denotes the derivative with respect to the radial distance , and is the cosine of the angle between the lineofsight direction and the separation vector,
(6) 
Eq. 2 can therefore be rewritten as
(7) 
We will assume that the peculiar velocity field is coupled to the density as in linear dynamics, and is determined by the void alone, so that
(8) 
where is the average mass density contrast within radius of the void centre,
(9) 
where is the mass density profile of the void, and , with the growth factor and the scale factor, is the linear growth rate of density perturbations.
Combining Eqs. 7 and 8, and using the fact that , we obtain
(10) 
The expansion of the second term on the RHS must be performed with care. For the voids in highly biased galaxy tracers that we consider in this paper, the mass density contrast is small enough that we can safely drop terms of second order or higher in and . However, note that this is separate from the assumption of the validity of linear dynamics, Eq. 8, and may not hold for other void populations, even if the dynamics is close to linear.
Crucially also, independent of assumptions about the size of and , the real space galaxy correlation cannot be assumed to be small. The very fact that voids are identified by selecting regions with very few or no galaxies ensures that in the void interior , and thus that terms proportional to and must be retained to linear order. The left panel of Figure 1 demonstrates this explicitly for the voids in our simulation.
Retaining these terms in the expansion, Eq. 10 reduces to
(11)  
where the difference between real space and redshift space coordinates and is
(12) 
to the same linear order in . The effect of this coordinate shift is neglected in the linear Kaiser theory for RSD in the galaxy correlation (Kaiser, 1987) because it only appears at second order in that derivation (see, e.g., Section 4.2 of Hamilton, 1998). However, in our case, the coordinate shift appears in at linear order, via the Taylor expansion
(13) 
The effect of the terms in the model above can be understood as follows. The velocity outflow from voids introduces a mapping between the real space correlation at separation and the redshift space correlation at separation , where along the line of sight, for . Since decreases towards the centre, this results in a ‘stretching’ of the void along the line of sight direction, in accordance with naive intuition for the RSD effect around underdensities. This effect is most important where the gradient of is steepest.
On the other hand, as noted and explained by Cai et al. (2016), the terms , arising from the Jacobian of the transformation, give rise instead to a ‘squashing’ along the line of sight.^{1}^{1}1In principle it is possible for these terms to also lead to a stretching for pathological choices of , but in such a case it is unlikely that a linear model based on Eq. 8 would be viable anyway. In practice, we do not observe such behaviour for any subset of the voids in our simulation. This effect is suppressed by a factor of and is therefore most important around the void edges and outside the void. The competing effects of the stretching and squashing terms can be seen in Figure 4. The combination results in a change in the sign of the quadrupole term within the void, going from negative to positive with increasing , as will be shown in Section 4.
3.2 Linear bias approximation
Eqs. 11 and 12 together provide the correct expression for the voidgalaxy crosscorrelation in redshift space assuming linear dynamics and at linear order in , . They form the base model that we will use to explain the simulation data. However, their use requires knowledge of the true mass density contrast and its integral version , which cannot be directly determined for voids in survey data. An approximation that is often made (Cai et al., 2016; Hamaus et al., 2017) is to assume that a simple linear bias relationship holds within the voids, so that and , where is the largescale linear bias factor determined from the galaxy clustering.
Under this approximation, the full model of Eqs. 11 and 12 can be rewritten as
(14)  
where
(15) 
and .
To explore the validity of this linear bias approximation within voids, we determine the linear bias factor for our mock galaxies by computing the (real space) galaxy power spectrum and the full matter power spectrum for the simulation box, and fitting at large scales, Mpc. For this value of , the residuals and compared to the linear bias approximation are shown in the right panel of Figure 1. Importantly, the deviations are (i) scale dependent, and (ii) large, with fractional differences of from the true value within the void radius.
It is important to emphasise that the deviation from the linear bias relationship seen can arise purely due to a statistical selection effect caused by the voidfinding process, independent of nonlinear or environmentdependent bias models (which we do not consider here). Where a linear bias characterises the galaxy clustering as above, the relationship between the local galaxy overdensity and the matter overdensity is still , where represents a stochastic term. If the stochastic term is unbiased, the conditional expectation value over the entire universe or simulation box is indeed
(16) 
However, voids are by construction selected as relatively rare regions with a small galaxy density. This necessarily introduces a selection effect in : regions with negative fluctuations in are more likely to be selected as voids than those with positive , thus biasing the mean relationship such that systematically overestimates the depth of the void where . Conversely, overestimates the height of the wall around the void edge where . This will be true even when , i.e. when the tracers used to identify the voids and determine are a random subset of the DM particles in the simulation.^{2}^{2}2Under certain simplifying assumptions about the nature of the stochastic term, this selection effect can be modelled analytically; see, e.g., Gruen et al. (2016), Sec. 3.1. In general, the void selection condition will introduce correlations between the stochastic terms at different radial distances, complicating the modelling. We will not pursue this further in the current work. The selection effect will however be much reduced or even absent if the tracer profile is measured using a different set of tracers to those used to select void locations, such as an overlapping galaxy sample with higher mean number density, or the set of all halos in the simulation from which the mock galaxy hosts are drawn.
For our purposes, the consequences of this failure of the linear bias relationship within voids are twofold:

the extent of the discrepancy is somewhat mitigated by the fact that the discrepancy from the linear bias value is largest in the regions where the term in Eq. 14 is small.
For the mocks we have both the matter and galaxy density fields, so we can contrast results for the real to redshift space mapping using the directly measured density field within voids as well as that inferred from the linear bias assumption, as done in Section 4. This allows us to isolate the effect of this assumption on the RSD model.
3.3 Multipole expansion
It is convenient to expand in terms of its multipoles,
(17) 
where are the Legendre polynomials of order . At linear order in and , only the monopole and quadrupole terms are nonzero. Using and , these can be calculated by direct integration for the model using either Eqs. 11 and 12 or Eqs. 14 and 15. For all numerical calculation of model multipoles presented in this paper, we use this approach.
However, approximate analytical forms can also be obtained by using the first two terms of the expansion in Eq. 13, together with Eqs. 11 and 17, to write
(18) 
and
(19) 
Corresponding versions of these equations can be obtained when also including the linear bias approximation.
A consequence that follows from these expressions is that the quadrupoletomonopole ratio does not factorise conveniently as claimed by Cai et al. (2016), and therefore unfortunately cannot be used as an estimator for the growth rate within the void interior.
3.4 Velocity dispersion and the streaming model
So far we have used a pure Kaiser model to describe : that is, we assumed that velocities around void centres exactly follow the coherent outflow described by the linear relationship in Eq. 8. In Section 4 we will show that this is a surprisingly good approximation for the mean outflow velocity on all scales, so the model of Eq. 11 provides a good qualitative description of the voidgalaxy correlation seen in simulation.
However, to enable a quantitative fit to the data, a more realistic model must account for the dispersion of galaxy velocities around this mean. To do this, we introduce a dispersion in the lineofsight galaxy velocities, such that
(20) 
where is the coherent radial component given by Eq. 8, and is a zeromean random variable with probability distribution function . This results in an integral for :
(21)  
where and are respectively distances transverse to and along the line of sight, , with and , and is as in Eq. 5. The effect of the dispersion in is primarily through shuffling the radial and transverse distances contributing to in the integral. The dispersion has a negligible effect on the Jacobian, as the contribution to the radial outflow averages to zero.
We will take the probability distribution function to have a Gaussian form,
(22) 
The dispersion is most generally a function of the radial separation scale, . The validity of the Gaussian assumption for may also vary with . We consider these questions in more detail in Section 4.
Note that Eq. 21 is not quite the same as an adaptation of the streaming model (e.g., Fisher, 1995; Scoccimarro, 2004; Reid & White, 2011) to the voidgalaxy case. In particular, as the coherence of density and velocity fields is accounted for by the inverse Jacobian, this is still fundamentally a linear model in all senses. We still assume the coherent outflow is determined by linear dynamics as before; in evaluating the Jacobian we still drop terms higher than linear order in and ; and the convolution with merely broadens the coordinate shift effect already present at linear order in Eqs. 11 and 12. In the limit of zero dispersion, where , this model reduces to Eq. 11. Eq. 21 is thus similar to the dispersion model used for the galaxy correlation (Hamilton, 1998) before the development of the full streaming model, except that the width of the velocity distribution is allowed to depend on scale.
Unlike in the case for the galaxy correlation, this linear dispersion model is already sufficient to provide an excellent fit to the data at all scales, as we show in Section 4.3. We may therefore safely leave the development of a full streaming model for the voidgalaxy case to future work.
3.5 Comparison with previous results in the literature
The basic linear model we have presented above differs in several key respects to other models for the voidgalaxy correlation in the literature (e.g. Paz et al., 2013; Cai et al., 2016; Hamaus et al., 2017; Achitouv et al., 2017; Achitouv, 2017). In Section 4 we will compare these models to measured in the simulation data and show that our model provides a significantly better description of the true voidgalaxy correlation. Before doing that, however, it is instructive to examine the reasons for the difference in the derivations.
Cai et al. (2016) (and subsequently Hamaus et al., 2017) follow a derivation similar to that described in Section 3.1, with three important differences: (1) they drop terms of order and in the expansion of Eq. 10, (2) they approximate contrary to Eq. 12, and (3) they assume the linear bias approximation (Section 3.2) holds within voids.^{3}^{3}3Cai et al. (2016) actually provide expressions assuming and , which is equivalent to using the linear bias assumption with . They therefore obtain the substantially simpler result
(23) 
which can be compared to Eqs. 11 and 14. This model gives expressions for the monopole
(24) 
and the quadrupole
(25) 
A simple way to see that these cannot be correct within the void interior is to consider the limiting case where close to the void centre, at . Here Eq. 24 predicts a redshift space monopole , which is unphysical. More generally, Eq. 23 is a poor description of everywhere that is appreciable, which in practice for the voids considered here means at least for all radial separations within the mean void radius.
In addition, Eq. 23 entirely misses the important stretching effect along the line of sight due to the mapping of voidgalaxy separations from real to redshift space. As we show in the next section, this causes a change in the sign of the quadrupole term within the void radius which is not captured by Eq. 25. In fact this change of sign was already seen in simulation data by Cai et al. (2016) (see Figures 13 in their paper), although not satisfactorily explained. This failure of the model is the primary reason why the growth rate estimator proposed in that paper fails for , leading to negative reconstructed values of the growth rate within voids (see the discussion in Cai et al., 2016).
On the other hand, in contrast to both our results and those of Cai et al. (2016), Hamaus et al. (2017) do not observe any change of sign in the quadrupole measured in BOSS galaxy data and associated galaxy mocks. This is because they use voids identified using the galaxy positions in redshift space—when this is done, neither of the fundamental assumptions used to obtain Eqs. 2 and 3 are valid, and so none of the theoretical models for discussed in this paper or elsewhere in the literature are applicable. We discuss this problem and its resolution in a separate work (Nadathur, Carter & Percival, in prep.).
Finally, a number of different authors (e.g. Paz et al., 2013; Hamaus et al., 2015, 2016; Cai et al., 2016; Achitouv et al., 2017; Achitouv, 2017) have used an analogy with the Gaussian streaming model for the galaxy autocorrelation to model :
(26) 
with . Comparison with Eq. 21 shows that this expression differs from our dispersion model. Unlike our dispersion model, it does not reduce to the linear expression, Eq. 11, in the limit of zero dispersion when .^{4}^{4}4This model does approximately match Eqs. 11 and 21 when the conditions are satisfied and the dispersion is small, (Cai et al., 2016). In practice this holds only outside the void boundary, . Also unlike our dispersion model, it provides a poor fit to the data, as we will show in Section 4.
The reason for this is simple: the streaming model in Eq. 26 has been derived for the case of linear, Gaussian fluctuations in the galaxy or matter autocorrelation (Fisher, 1995; Scoccimarro, 2004), and cannot be extended to the voidgalaxy crosscorrelation. Not only is the moment generating function (see Scoccimarro, 2004) different for the voidgalaxy case, the derivation of the Gaussian streaming model explicitly assumes that the only nonzero cumulant of the density field is the second. However, the voidgalaxy case involves constrained averages of the fields at specially selected locations rather than over all space, so that by definition the expectation value within the void. As a result the analogy is invalid and Eq. 26 does not describe the correct streaming model for the voidgalaxy crosscorrelation.
4 Simulation results
4.1 Void density profiles
We measure the density profiles for each void using the full resolution DM particle output of the simulation at the snapshot in equallyspaced radial bins out to separations of Mpc. We use the same binning for each void and do not rescale separations by the void radius; therefore, these individual profiles can be simply averaged over all voids to give the stack average profile . We measure , which is equivalent to the average real space galaxy density profile around voids, in the same way after substituting galaxy tracers for the DM particles. To determine and we interpolate and respectively and numerically evaluate the corresponding integrals.
The profiles obtained are shown in Figure 1, where for context we have shown radial separations in units of the mean radius of all voids in the stack, Mpc. The importance of the relative sizes of the different terms and the failure of the linear bias approximation have already been discussed in Section 3.
4.2 Velocity profiles and dispersion
We measure the average (stacked) galaxy velocity profile around the void centres in the simulation as
(27) 
where is the real space separation vector between the th void and th galaxy and the velocity of that galaxy, and the sum runs over all voidgalaxy pairs with in the range . As before, we use equallyspaced bins out to a maximum separation of Mpc. We also measure the dispersion in the lineofsight velocity component, , defined by
(28) 
where is the line of sight direction, which is the same for all galaxies in the planeparallel approximation and taken to be along the axis of the simulation box, is the mean line of sight velocity component at , and the sum is over all voidgalaxy pairs in the given separation bin as before.
Note that Eqs. 27 and 28 weight each voidgalaxy pair equally. An alternative way to define these quantities, which has sometimes been used in the literature, is to calculate the average and dispersion separately for each void, and then average the results over all voids. This is equivalent to weighting each void equally: this procedure leads to smaller measured values of and at small separations, since most voids do not have any galaxies in interior bins and thus contribute zero to the average. At large both methods will largely agree. However, for comparison with the voidgalaxy crosscorrelation, which also weights each voidgalaxy pair equally, our method is the more appropriate.
Finally, to control for possible statistical effects due to the void selection criterion as seen in the relationship between and , we also measure the average halo velocity profile around void centres, using the equivalent of Eq. 27 but for the full halo catalogue. This has many more objects than the mock galaxy catalogue and therefore avoids the problem of low tracer statistics close to void centres. Under the reasonable assumption that the halo velocity field is unbiased with respect to the dark matter velocity, we refer to this quantity as .
Figure 2 plots the results as a function of distance from the void centre, together with the linear velocity model for obtained from Eq. 8 (solid black line). The linear model is calculated using the fiducial value of the growth rate for the given simulation and redshift, . The first conclusion that can be drawn from this is that linear dynamics gives an extremely good description for the mean , with deviations at over the entire range of scales tested. It is also a good model for the mean galaxy radial velocity profile at distances Mpc, where and coincide.
Closer to the void centres, , starts to deviate from both the linear model and . While some scatter is expected where the errors in become large due to the small number of voidgalaxy pairs (note that the two interiormost bins have no galaxies at all), in the range Mpc the difference is statistically significant. The fact that it coincides with the ‘kink’ visible in the profile in Figure 1 and a similar, though less pronounced, kink that is also present at the same distance in , suggests that this is a physical consequence of some aspect of the void selection algorithm. We leave a fuller explanation of this effect for future work.
Recently, Achitouv (2017) proposed a semiempirical nonlinear correction term to modify the linear velocity relationship of Eq. 8, which effectively reduces the predicted outflow term within voids. Figure 2 shows that such a correction term is disfavoured in our data, since is very close to the linear prediction and never less than it. However, note that incorrectly assuming a linear bias for the density profiles, in Eq. 8 (shown by the dashed line in Figure 2) does indeed consistently overestimate the true velocity outflow. This can be understood by reference to Figure 1, which shows that assuming a linear bias overestimates the true matter deficit within the void.
The lower panel of Figure 2 shows the velocity dispersion as a function of radial distance. Errors are estimated assuming velocities follow a normal distribution. The dispersion is approximately constant outside the void radius but decreases in the interior region, with a sharp drop in the same region as the measured radial velocity starts to deviate from the linear prediction. In the interiormost bins, the small number of voidgalaxy pairs means the dispersion cannot be reliably estimated, so these points are excluded.
Where the dispersion can be measured, it is always significantly larger than the mean radial velocity. This is in agreement with Achitouv (2017) but differs from the results of Hamaus et al. (2015) and Cai et al. (2016), who report a much smaller dispersion within the mean void radius. It is possible this is due to the use of equal weighting for each void rather than each voidgalaxy pair in the latter two papers.
Figure 3 shows the measured probability distribution function for line of sight velocities in four different bins of , compared to the Gaussian distribution assumed in Eq. 22 with set to the corresponding value of . The distribution is symmetric around zero at all scales, and close to Gaussian at small scales, although at large separation scales exponential wings are seen which deviate from Gaussianity. At these large scales the impact of dispersion itself is small. Therefore the Gaussian assumption, while not strictly correct, still provides a good fit—as we show later. Note that this is quite different from the case for the distribution of pairwise velocities for galaxies, which is highly skewed at small scales and far from Gaussian at all scales (Scoccimarro, 2004).
4.3 Redshift space correlation function
We measure the voidgalaxy crosscorrelation function in the simulation data using a modified version of the CUTE correlation function code (Alonso, 2012).^{5}^{5}5http://members.ift.uamcsic.es/dmonge/CUTE.html Note that, as discussed in Section 3, in order to ensure compatibility between data and the theoretical modelling, void positions are identified using the real space galaxy positions in the simulation.
Figure 4 shows the measured redshift space crosscorrelation as a function of the transverse and lineofsight distances, compared to the prediction from our model, Eq. 11, and that of the model of Cai et al. (2016) and Hamaus et al. (2017), Eq. 23. The qualitative features visible in the data can be understood on the basis of the equations derived in Section 3.1 above. Close to the void centre, the measured shows approximate spherical symmetry. At intermediate distances, the effect of the coordinate shift discussed in Section 3.1 causes an elongation effect, stretching the contours of equal along the line of sight direction. At large distances, Mpc, this effect is reversed, with a relative squashing of the contours of along the line of sight direction.
To better understand the relative merits of the two theoretical models, Figure 5 shows the residual difference between the measured and the predicted correlation functions in the two cases. It is clear that the model of Cai et al. (2016); Hamaus et al. (2017) is a very poor fit to the data at all radial separations within the mean void radius , and particularly so very close to the centre, where Eq. 23 predicts unphysical values of . By contrast, the inclusion of the factors in Eq. 11 correctly ensures that the theory prediction matches the measured value at the void centre—where in fact RSD effects are absent, —and also produces a much better fit to the data at all distances within the void interior.
Nevertheless, Figure 5 shows that even for the improved model, when the velocity dispersion term is not included the theory residuals are still potentially significant. To perform a quantitative analysis, we instead perform a multipole expansion of the angular correlation function from the simulation data to extract the monopole and quadrupole terms using Eq. 17 for comparison with theory. The multipoles are measured in 50 radial bins in the range Mpc, and we use 100 bins in . Note that we do not rescale radial distances by the void radius, thus using the same bin sizes for each void. Such a rescaling would constitute a strong assumption of selfsimilarity in void profiles, which is not justified (Nadathur & Hotchkiss, 2015a, b; Nadathur et al., 2017).
To estimate the error in this measurement we use a jackknife resampling technique. We divide our simulation box into nonoverlapping cubic subboxes, each measuring Mpc on a side, and determine the correlation function and multipoles excluding each subbox in turn, which we combine into the data vector , for . The covariance matrix is then determined as
(29) 
Figure 6 shows the resulting matrix of crosscorrelation coefficients. For the monopole, this matrix is far from diagonal, showing strong correlation between bins even at large separations. For the quadrupole it is closer to diagonal, although correlations between neighbouring bins are still significant, especially at large separations . The values for the fit of a given theoretical model to the data are calculated as
(30) 
Figure 7 shows the measured values of the monopole and quadrupole from the simulation data, compared to theory predictions derived in the previous section evaluated using the fiducial growth rate . For visual clarity, the monopole values are shown as the difference with respect to the real space version, . Two important features are immediately apparent: the quadrupole changes sign within the void interior, with a negative dip at intermediate values of changing to a broad feature at Mpc; and the redshift space monopole differs from the real space version only at intermediate , with both at the very centre (where ) and around the mean void radius.
The dashed curves in each panel show the predictions of the basic linear model derived from Eq. 11 without accounting for velocity dispersion: these capture the main qualitative features above but do not provide a good fit to the data, as expected from visual inspection of the residuals in Figure 5. The solid curves show the prediction for the linear dispersion model of Eqs. 21 and 22, where we have used the radial dependence of the dispersion measured in the data, shown in Figure 2. This model provides an excellent fit to the data at all scales: the values for this model are for the monopole alone, for the quadrupole alone, and for both combined.
Finally, the dotted curves in Figure 7 show the predictions for and according to the model presented by Cai et al. (2016); Hamaus et al. (2017), Eqs. 24 and 25. Unsurprisingly, this model fails qualitatively and quantitatively to describe the data at any point within the void interior. Outside the mean void radius Mpc, the multipoles from this model approach the correct expression: this explains the observation of Cai et al. (2016) that the linear growth rate estimator proposed in that paper works only in the region .
Table 1 compares the goodness of fit for the various models discussed: here ‘lin.+dispersion’ refers to the linear dispersion model above, ‘lin.+dispersion+bias’ to this model with the addition of the linear bias assumption , and ‘lin. only’ to the model without dispersion, Eq. 11. The fit for the streaming model of Eq. 26 uses (without the linear bias assumption) to obtain , and measured from simulation, so can be directly compared to the dispersion model without bias. The fit for the Cai et al. (2016) model includes the bias assumption, as in Eq. 23.
Model  d.o.f  

monopole  quadrupole  combined  
lin.+dispersion  0.86  1.13  1.24 
lin.+dispersion+bias  2.14  1.58  2.34 
lin. only  5.35  4.76  6.48 
streaming model, Eq. 26  5.59  3.27  4.72 
Cai et al. (2016), Eq. 23  315.6  45.2  203.7 
5 Measuring the growth rate
We have seen that the fiducial linear dispersion model for the RSD in the voidgalaxy correlation provides an excellent description of the multipoles measured from simulation at all scales. An interesting question is how this model might be used to determine the value of growth rate within voids from fits to data.
As discussed in Section 3.3, it is not possible to construct an estimator for from the ratio of measured redshift space multipoles. Therefore to fully specify the model and calculate the multipoles for a given value of in principle requires knowledge of three functions: , , and . All three depend strongly on the voidfinding algorithm and the choice of centre for stacking, as well as potentially having a smaller cosmological dependence.
Of these, the most important input to the model is the real space correlation function , which cannot be measured directly. A plausible procedure would be to calibrate using fits to the real space void profile in simulated mocks before use in the RSD model. This is the procedure followed by Hamaus et al. (2015, 2016); Achitouv et al. (2017); Achitouv (2017). Various fitting formulae have been described in the literature, e.g. Hamaus et al. (2014); Ricciardelli et al. (2014); Nadathur et al. (2015); Nadathur & Hotchkiss (2015a, b); Barreira et al. (2015); Cautun et al. (2016). However, there is no consensus on any particular fitting form, and in any case the fit will strongly depend on the voidfinding and stacking algorithms as well as on the properties of the galaxy sample in question (Nadathur & Hotchkiss, 2015b), so the calibration needs to be performed on a casebycase basis. Alternative methods to obtain are by deprojecting the redshift space monopole (Pisani et al., 2014) or by using reconstruction techniques (Nadathur, Carter & Percival, in prep.).
In principle the matter overdensity profile can be determined independently of through void lensing measurements (Krause et al., 2013; Melchior et al., 2014; Clampitt & Jain, 2015; Sánchez et al., 2016). Alternatively, it could also be calibrated from simulations as for . Use of the linear bias assumption is however not a good approximation if the voidgalaxy correlation is to be determined for the same galaxy tracers used to identify the voids. This is because of the selection effect problem described in Section 3.2. For instance, using this approximation for our simulation data leads to errors in the predicted quadrupole of up to .
In the previous section we took the dispersion directly from the simulation data in order to demonstrate the goodnessoffit of the fiducial model. It is possible that future work on a complete streaming model for the voidgalaxy case may allow this function to be determined on the basis of and . A popular alternative (e.g. Hamaus et al., 2015; Cai et al., 2016; Achitouv, 2017) is to parametrise the dispersion using some functional form with one or more free parameters, which are to be marginalised over. The simplest possible such parametrisation is to take it to be a constant,
(31) 
Figure 2 suggests that this is a reasonable approximation separately in the void interior and exterior regions. Sensitivity to the dispersion also decreases in the void exterior, where . Thus it might be hoped that this approximation is sufficient to reconstruct .
To test this, we fit multipoles for the linear dispersion model of Eq. 21 to the simulation data with and as two free parameters. and are taken from the fits to the simulation data. Unlike in some previous works, our model works well at all scales within voids, so we use multipole data at all separations . Figure 8 shows the resulting and confidence level contours determined from fitting the quadrupole (blue) and monopole (red) separately. The fit to the quadrupole provides an unbiased reconstruction of the fiducial growth rate, with the reconstructed value at c.l. after marginalizing over . This is consistent with the fiducial value , and corresponds to a precision for our simulation volume of . The recovered constraint on is also broadly consistent with the results in Figure 2.
However, fitting to the monopole gives a biased reconstruction, , which is more than from the fiducial value. This can be understood in the following way. Since voidfinding identifies regions of low galaxy density, it effectively acts as an indirect selection cut on the real space monopole . For spherical underdensity voidfinders (e.g. Padilla et al., 2005; Cai et al., 2015) this cut is easy to understand and operates directly at a fixed distance scale, but in fact it will apply to some extent for all voidfinding algorithms, including our watershedbased algorithm. The effect of such a cut is to introduce a form of Malmquist bias: at relatively small amplitudes, the individual redshift space monopoles for voids in the sample are driven by voids that have a larger (i.e., more negative) real space monopole but scatter to smaller redshift space amplitudes. This has the potential to bias the recovered value of f, which we should then expect to be low when measured from the average redshiftspace monopole , as observed. Void selection has less of an effect on the quadrupole, since is significantly less correlated with , and thus the quadrupole gives an unbiased estimate of . An additional complication is that the monopole is more sensitive to the precise value of the dispersion, as can be seen from the stronger degeneracy between and , and that therefore the use of the simplistic assumption of a constant dispersion, while sufficient for the quadrupole, is less appropriate for the monopole.
Finally, we tested the effect on the measurement of the growth rate of using the linear bias assumption . Due to the failures described above, this leads to a strongly biased estimator of , with the reconstructed value being more than smaller than the fiducial for fits to both the quadrupole and monopole.
6 Conclusions
Measurement of redshift space distortions in the voidgalaxy correlation function is an important tool that can be used to test for possible environmental dependence of the growth rate of structures. In particular, the growth rate in the lowest density regions close to void centres might contain information about possible nonstandard theories of gravity. However, reconstructing the growth rate in voids requires a model for which can be trusted in these regions.
We have derived a configurationspace model for the voidgalaxy correlation in redshift space using linear theory, which we characterize in terms of the multipole moments of the correlation. Our model accounts for several terms that are important within voids but have previously been neglected. As a result we are able to account for important physical effects that had not been appreciated, including the change in sign of the quadrupole term within the void, indicating a turnover point between stretching and squashing of the contours of the correlation function along the line of sight direction. The model can be broadened to include a dispersion in galaxy velocities along the line of sight; the dispersion model thus obtained differs from the streaming model used in previous studies, which was based on an inappropriate application of the formula for the Gaussian streaming model derived for the galaxy autocorrelation.
Comparing our model predictions to measurements of using purposebuilt void and galaxy catalogues at redshift 0.52 in the Big MultiDark simulation shows that the linear dispersion model provides an excellent fit to the data for all values of the voidgalaxy separation. This is an important contrast to the modelling of RSD in the galaxy correlation, where nonlinear effects are very important at small pair separations and complicate the use of smallscale data. Comparison with the data shows that our linear model for the voidgalaxy correlation performs significantly better than others in the literature, especially within the low density region of interest.
A consequence of our results is that the ratio of the redshift space monopole and quadrupole cannot be used as an estimator for the growth rate within voids, contrary to previous claims (Cai et al., 2016; Hamaus et al., 2017). This is because the linear model proposed in these papers fails at . Determination of the growth rate using the correct expression we provide requires knowledge of the real space correlation function , which can be taken from fits to simulations or possibly reconstructed from redshift space data. Assuming a constant dispersion taken as a free parameter and using values of and determined from fits to the simulation data, we found that fitting to the redshift space quadrupole provides an unbiased estimate of the true growth rate with precision from our simulation volume of . This precision may improve with more optimal choices of the void selection criteria, as discussed in Section 2. The monopole data has higher constraining power than the quadrupole, but leads to a biased reconstruction of the true , due to a form of Malmquist bias arising from void selection.
Our results also highlight two very important points which have relevance to practical applications of this method to measure RSD effects within voids. Firstly, if the galaxy tracers used to measure the RSD effects are the same as those used to identify the voids, an assumption of linear galaxy bias within the void interior is not appropriate due to statistical selection effects on the void size scale. Assuming a constant linear bias leads to incorrect predictions from the model, which possibly require the addition of empirical nonlinear correction terms (Achitouv, 2017), even though in reality linear theory continues to provide a good description of the dynamics. In order to overcome this problem, in practical terms it would be better to use two differently biased galaxy tracers covering the same volume to identify void regions and to measure the RSD effects.
Secondly, we highlighted how the assumptions inherent in the derivation of the model require knowledge of the real space positions of voids, such that all RSD effects are due to motions of the galaxies only. In general if the voids are identified using galaxies in redshift space, this will not be the case (Chuang et al., 2017). Practical methods to reconstruct the real space void positions in data will be discussed in a companion paper.
Finally, a comment on the relative merits of using the voidgalaxy correlation and the galaxy autocorrelation to measure RSD is in order. It is true that unlike for the galaxy RSD case, the model presented in this work is based only on linear theory, and fits the simulation data extremely well at all pair separations. It is particularly noteworthy that dispersion effects can be easily accommodated within a linear model. These appear to be significant advantages over the galaxy RSD case, for which quasilinear modelling is required and smallscale data is excluded when performing fits. However, since the voidgalaxy correlation uses only those galaxies in the neighbourhood of voids, this method also discards some of the available data. One therefore should not necessarily expect it to outperform the standard analysis, as has sometimes been claimed (Hamaus et al., 2017). Instead an advantage of the voidgalaxy RSD measurement is specifically to constrain possible environmental dependence of the growth rate in lowdensity regions. It also provides a complementary measurement technique that is sensitive to different systematics—in particular, if the real space can be reconstructed directly from data rather than fit from simulations, this method constrains directly instead of the combination .
7 Acknowledgements
We thank Yanchuan Cai for helpful correspondence, and Davide Bianchi, Florian Beutler, Hans Winther and Rossana Ruggeri for stimulating discussions. SN acknowledges funding from the Marie SkłodowskaCurie Actions under the H2020 Framework of the European Commission, project 660053 COSMOVOID. WJP acknowledges support from the European Research Council through the Darksurvey grant 614030, and from the UK Science and Technology Facilities Council grant ST/N000668/1 and the UK Space Agency grant ST/N00180X/1. The BigMultiDark simulations were performed on the SuperMUC supercomputer at the LeibnizRechenzentrum in Munich, using resources awarded to PRACE project number 2012060963. Numerical computations were performed on the EREBOS, THEIA and GERAS clusters at the Leibniz Institut für Astrophysik (AIP).
References
 Achitouv (2017) Achitouv I., 2017, Phys.Rev.D, 96, 083506
 Achitouv et al. (2017) Achitouv I., Blake C., Carter P., Koda J., Beutler F., 2017, Phys.Rev.D, 95, 083502
 Alonso (2012) Alonso D., 2012, ArXiv eprints, arXiv:1210.1833
 Alonso et al. (2017) Alonso D., Hill J. C., Hložek R., Spergel D. N., 2017, ArXiv eprints
 Barreira et al. (2015) Barreira A., Cautun M., Li B., Baugh C. M., Pascoli S., 2015, J. Cosmol. Astropart. Phys., 8, 028
 Beutler et al. (2012) Beutler F. et al., 2012, MNRAS, 423, 3430
 Beutler et al. (2017) Beutler F. et al., 2017, MNRAS, 466, 2242
 Cai et al. (2017) Cai Y.C., Neyrinck M., Mao Q., Peacock J. A., Szapudi I., Berlind A. A., 2017, MNRAS, 466, 3364
 Cai et al. (2015) Cai Y.C., Padilla N., Li B., 2015, MNRAS, 451, 5555
 Cai et al. (2016) Cai Y.C., Taylor A., Peacock J. A., Padilla N., 2016, MNRAS, 462, 2465
 Cautun et al. (2016) Cautun M., Cai Y.C., Frenk C. S., 2016, MNRAS, 457, 2540
 Chuang et al. (2017) Chuang C.H., Kitaura F.S., Liang Y., FontRibera A., Zhao C., McDonald P., Tao C., 2017, Phys.Rev.D, 95, 063528
 Clampitt & Jain (2015) Clampitt J., Jain B., 2015, MNRAS, 454, 3357
 Deffayet et al. (2009) Deffayet C., Deser S., EspositoFarèse G., 2009, Phys.Rev.D, 80, 064015
 Dvali et al. (2000) Dvali G., Gabadadze G., Porrati M., 2000, Physics Letters B, 485, 208
 Fisher (1995) Fisher K. B., 1995, ApJ, 448, 494
 Gottloeber & Klypin (2008) Gottloeber S., Klypin A., 2008, ArXiv eprints
 Granett et al. (2008) Granett B. R., Neyrinck M. C., Szapudi I., 2008, ApJ, 683, L99
 Gruen et al. (2016) Gruen D. et al., 2016, MNRAS, 455, 3367
 Guzzo et al. (2008) Guzzo L. et al., 2008, Nature, 451, 541
 Hamaus et al. (2017) Hamaus N., Cousinou M.C., Pisani A., Aubert M., Escoffier S., Weller J., 2017, J. Cosmol. Astropart. Phys., 7, 014
 Hamaus et al. (2016) Hamaus N., Pisani A., Sutter P. M., Lavaux G., Escoffier S., Wandelt B. D., Weller J., 2016, Physical Review Letters, 117, 091302
 Hamaus et al. (2015) Hamaus N., Sutter P. M., Lavaux G., Wandelt B. D., 2015, J. Cosmol. Astropart. Phys., 11, 036
 Hamaus et al. (2014) Hamaus N., Sutter P. M., Wandelt B. D., 2014, Phys. Rev. Lett., 112, 251302
 Hamilton (1998) Hamilton A. J. S., 1998, in Hamilton D., ed., Astrophysics and Space Science Library Vol. 231, The Evolving Universe. p. 185
 Hotchkiss et al. (2015) Hotchkiss S., Nadathur S., Gottlöber S., Iliev I. T., Knebe A., Watson W. A., Yepes G., 2015, MNRAS, 446, 1321
 Howlett et al. (2015) Howlett C., Ross A. J., Samushia L., Percival W. J., Manera M., 2015, MNRAS, 449, 848
 Hu & Sawicki (2007) Hu W., Sawicki I., 2007, Phys.Rev.D, 76, 064004
 Jennings et al. (2011) Jennings E., Baugh C. M., Pascoli S., 2011, MNRAS, 410, 2081
 Kaiser (1987) Kaiser N., 1987, MNRAS, 227, 1
 Kitaura et al. (2016) Kitaura F.S. et al., 2016, Physical Review Letters, 116, 171301
 Klypin & Holtzman (1997) Klypin A., Holtzman J., 1997, ArXiv eprints
 Klypin et al. (2016) Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457, 4340
 Kovács et al. (2017) Kovács A. et al., 2017, MNRAS, 465, 4166
 Krause et al. (2013) Krause E., Chang T.C., Doré O., Umetsu K., 2013, ApJ, 762, L20
 Kravtsov et al. (1997) Kravtsov A. V., Klypin A. A., Khokhlov A. M., 1997, ApJS, 111, 73
 Manera et al. (2013) Manera M. et al., 2013, MNRAS, 428, 1036
 Matsubara (2008) Matsubara T., 2008, Phys.Rev.D, 78, 083519
 Melchior et al. (2014) Melchior P., Sutter P. M., Sheldon E. S., Krause E., Wandelt B. D., 2014, MNRAS, 440, 2922
 Nadathur (2016) Nadathur S., 2016, MNRAS, 461, 358
 Nadathur & Crittenden (2016) Nadathur S., Crittenden R., 2016, ApJ, 830, L19
 Nadathur & Hotchkiss (2015a) Nadathur S., Hotchkiss S., 2015a, MNRAS, 454, 2228
 Nadathur & Hotchkiss (2015b) Nadathur S., Hotchkiss S., 2015b, MNRAS, 454, 889
 Nadathur et al. (2017) Nadathur S., Hotchkiss S., Crittenden R., 2017, MNRAS, 467, 4067
 Nadathur et al. (2015) Nadathur S., Hotchkiss S., Diego J. M., Iliev I. T., Gottlöber S., Watson W. A., Yepes G., 2015, MNRAS, 449, 3997
 Neyrinck (2008) Neyrinck M. C., 2008, MNRAS, 386, 2101
 Nicolis et al. (2009) Nicolis A., Rattazzi R., Trincherini E., 2009, Phys.Rev.D, 79, 064036
 Padilla et al. (2005) Padilla N. D., Ceccarelli L., Lambas D. G., 2005, MNRAS, 363, 977
 Paz et al. (2013) Paz D., Lares M., Ceccarelli L., Padilla N., Lambas D. G., 2013, MNRAS, 436, 3480
 Peacock et al. (2001) Peacock J. A. et al., 2001, Nature, 410, 169
 Pisani et al. (2014) Pisani A., Lavaux G., Sutter P. M., Wandelt B. D., 2014, MNRAS, 443, 3238
 Prada et al. (2012) Prada F., Klypin A. A., Cuesta A. J., BetancortRijo J. E., Primack J., 2012, MNRAS, 423, 3018
 Reid et al. (2012) Reid B. A. et al., 2012, MNRAS, 426, 2719
 Reid & White (2011) Reid B. A., White M., 2011, MNRAS, 417, 1913
 Ricciardelli et al. (2014) Ricciardelli E., Quilis V., Varela J., 2014, MNRAS, 440, 601
 Riebe et al. (2013) Riebe K. et al., 2013, Astronomische Nachrichten, 334, 691
 Samushia et al. (2012) Samushia L., Percival W. J., Raccanelli A., 2012, MNRAS, 420, 2102
 Sánchez et al. (2016) Sánchez C. et al., 2016, ArXiv eprints
 Scoccimarro (2004) Scoccimarro R., 2004, Phys.Rev.D, 70, 083007
 Springel (2005) Springel V., 2005, MNRAS, 364, 1105
 Taruya et al. (2010) Taruya A., Nishimichi T., Saito S., 2010, Phys.Rev.D, 82, 063522
 Zheng et al. (2007) Zheng Z., Coil A. L., Zehavi I., 2007, ApJ, 667, 760