# Modeling the environmental dependence of the growth rate of cosmic structure

Ixandra Achitouv, Yan-Chuan Cai

► To cite this version:

Ixandra Achitouv, Yan-Chuan Cai. Modeling the environmental dependence of the growth rate of cosmic structure. Physical Review D, 2018, 98 (10), pp.103502. 10.1103/PhysRevD.98.103502. hal-01827997

HAL Id: hal-01827997

https://hal.science/hal-01827997

Submitted on 18 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Modeling the environmental dependence of the growth rate
of cosmic structure

Ixandra Achitouv¹,² and Yan-Chuan Cai³,\*

¹Laboratoire Univers et Théories (LUTh), UMR 8102 CNRS, Observatoire de Paris,
Université Paris Diderot, 5 Place Jules Janssen, 92190 Meudon, France

²APC, Univ Paris Diderot, CNRS/IN2P3, CEA/Irfu, Obs de Paris, Sorbonne Paris Cit, France

³Institute for Astronomy, University of Edinburgh, Royal Observatory,
Blackford Hill, Edinburgh, EH9 3HJ, United Kingdom

<figure><img src="image_2.png" /></figure>

(Received 12 June 2018; published 2 November 2018)

The growth rate of cosmic structure is a powerful cosmological probe for extracting information on the gravitational interactions and dark energy. In the late-time Universe, the growth rate becomes nonlinear and is usually probed by measuring the two-point statistics of galaxy clustering in redshift space up to a limited scale, retaining the constraint on the linear growth rate *f*. In this paper, we present an alternative method to analyze the growth of structure in terms of local densities, i.e., *f*(*Δ*). Using *N*-body simulations, we measure the function of *f*(*Δ*) and show that structure grows faster in high-density regions and slower in low-density regions. We demonstrate that *f*(*Δ*) can be modeled using a log-normal Monte Carlo random walk approach, which provides a means to extract cosmological information from *f*(*Δ*). We discuss prospects for applying this approach to galaxy surveys.

DOI: 10.1103/PhysRevD.98.103502

The growth rate of cosmic structure contains important information on the matter-energy content of the Universe and the gravitational interactions that shape the cosmic web. A powerful way to extract this information is to use redshift-space distortions (RSDs) in galaxy clustering (e.g., Refs. [1–4]), or in the cross-correlation between clusters/voids and galaxies [5–8]. However, when using RSDs, among other cosmological probes, we are limited by the accuracy of our model to reproduce complex patterns in the galaxy clustering on small scales. Hence, we are often forced to throw away data in the nonlinear regime in order to extract unbiased cosmological information—in this case, the linear growth rate (e.g., Refs. [8–12]). One way to overcome this issue is to use perturbative approaches to model the global clustering in the quasin nonlinear regime down to a certain small scale where models break down. While nonlinear modeling allows us to extract an unbiased value of the linear growth rate, in principle, two-point statistics such as the correlation function are sensitive to the variance of the field. Applying them to a nonlinear field will not allow us to extract all of the information. This is because a nonlinear density field is usually non-Gaussian, and cannot be fully characterized by its variance. One can use higher-order statistics such as three-point or four-point correlation functions to regain the information beyond the variance, but this is currently computationally expensive.

In this study, we propose a different approach to the same problem: instead of measuring the globally averaged linear growth rate *f* at different scales by forward modeling the nonlinear growth of the matter power spectrum/correlation function, we accept that the growth of structure depends on local densities and aim to model this dependency, i.e., *f*(*Δ*), where *Δ* = *ρ* / ρ̄ − 1 is the local density contrast. To do this, we analyze the growth rate using numerical simulations in and around overdense and underdense regions and show how it can be predicted as a function of local density and for a given cosmology. This prediction relies on log-normal Monte Carlo random walks, a method introduced in Ref. [13]. We find that our model is successful in tracking the evolution of the growth rate in different local density environments. This, in principle, provides an independent method to extract cosmological information from the quasilinear and nonlinear regimes. Our method of understanding the nonlinear growth is in the same spirit as modeling the distribution of densities within spheres [14–16], density split statistics [17,18], the position-dependence power spectrum [19,20], and the modeling of the nonlinear aspect of the baryon acoustic oscillations [21,22]. A more complete study will be presented in a companion paper.

We perform our analysis using *N*-body simulations
from the DEUS Consortium. These are described in
Refs. [23–25] and are publicly available. These simulations
are run in a ΛCDM model with the WMAP 5-year
cosmology [26] with (w = −1; Ω_m = 0.26; σ_8 = 0.79).

<footnote>ixandra.achitouv@obspm.fr</footnote>
They have box lengths of 648 h⁻¹ Mpc with 1024³ particles. They are generated using the RAMSES code [27]; halos are found using a friend-of-friends finder with link length b = 0.2, [28] and cover a range of masses M ~ [10¹²–10¹⁵] h⁻¹ M☉.

We first identify regions of different density contrast Δ(R) (i.e., environment) in the simulations, where R is the radius of the region. We follow the method presented in Ref. [13] to identify low-density regions, i.e., voids. This algorithm imposes density thresholds at the radius of our choice, therefore allowing the flexibility to represent a large variety of void profiles. Here we choose R_v = 20 h⁻¹ Mpc (motivated by the resolution of the simulation) and the same criteria for the voids as those used in Refs. [8,12,13]: δ(R_1) < -0.9, δ(R_2) < -0.7, and δ(R_v = R_v + dR) > δ(R_v), where R_1 = 1.5 h⁻¹ Mpc, R_2 = 3 h⁻¹ Mpc, and δ is the density contrast of the halo field at R; see Ref. [13] for more details. Note that the choice of the criteria for the voids is not important for the outcome of this analysis as long as they are able to sample a wide range of Δ, covering Δ ~ -1. This also applies to our selection for overdense regions.

We run this void finder on the halo catalog and find ~2300 void centers. We measure the dark matter density profiles around our selected void centers to avoid complications due to the halo bias. We select the overdense regions by randomly sampling positions of dark matter particles belonging to halos above the mass resolution at z = 0, until we reach the same number of overdensities as the number of voids, to make sure that these two samples have similar noise properties. Keeping the same comoving coordinates for the under-/overdense regions fixed (identified from z = 0), we measure the evolution of the density profiles at redshifts z = {0.00; 0.05; 0.11; 0.67, 1.50} (corresponding to scale factors a = {1.00, 0.95, 0.90, 0.60, 0.40}, respectively).

In Fig. 1 we show the mean matter density profiles of these over-/underdense patches (dots) at different redshifts. From these profiles we measure numerically the growth rate within the radius R at a = 0.95 using three consecutive snapshots at a = {1.00, 0.95, 0.90} by computing

$$f(R) = \frac{d \ln \Delta(R)}{d \ln a}, \qquad (1)$$

where Δ(R) is the cumulative density contrast. The comoving coordinates of the centers of our under-/overdense regions are kept unchanged at different epochs.

The resulting growth rates are shown in Fig. 2 by the black data points. We bin the f values according to their local density Δ to show the values of fσ₈ as a function of the local density, where σ₈ = 0.79 is a constant. Note that the density contrasts of different scales R may end up in the same bin of Δ. In this sense, the behavior of f(Δ) is no longer an explicit function of the scale R, but rather depends solely on the local density Δ, which could be

<figure><img src="image_3.png" /><figcaption>FIG. 1. Cumulative matter density profiles around overdense (red) and underdense (blue) regions, measured from ΛCDM N-body simulations at different redshifts indicated by the legend.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 2. Growth rate parameter fσ₈ measured in a range of regions characterized by their density contrasts Δ from a ΛCDM N-body simulation (dots with black errors). The blue error bars correspond to the expectation of the statistical errors for the upcoming TAIPAN survey covering a volume of 1.3 (h⁻¹ Gpc)³ [29]. The errors are expected to shrink by another factor of ~4 and ~6 for the DESI LRG and ELG surveys, respectively [30]. The red curve shows our model prediction from the log-normal MCRW approach. The horizontal dashed line shows the linear expectation. The vertical dashed line indicates the mean density of the Universe.</figcaption></figure>
contributed by perturbations of different scales. The error bars correspond to the standard deviation computed from the mean measurements of 64 subcubes of length 162 h⁻¹ Mpc. The horizontal dashed line corresponds to the linear growth rate. Although Eq. (1) has a logarithmic divergence for Δ → 0, we can see how the growth rate varies compared to the linear one, indicated by the horizontal dashed line in Fig. 2. The growth of structure slows down in low-density regions and speeds up in high-density regions. While f(R) is expected to reach a linear value on large scales, the growth rate in terms of Δ only crosses the linear value when |Δ| is small. It is therefore important to go beyond a single value of the linear f by modeling the whole spectrum of f(Δ) to effectively extract cosmological information.

In general, we expect the overall averaged growth rate on small scales to be higher than in linear theory. This is because the amplitudes of the late-time matter power spectrum/correlation function tend to be higher than the linear version on small scales. These higher amplitudes must arise from a higher growth rate. This suggests that the larger/smaller growth rate in over-/underdense regions seen in Fig. 2 does not exactly cancel out for the global average on these scales. In fact, the matter power spectrum/ correlation function on small scales is dominated by high-density regions. Therefore, the branch of the curve with Δ > 0 shown in Fig. 2 contributes more to the global average growth rate than the Δ < 0 branch does. This is also consistent with the general trend of the inferred value of the growth rate from redshift-space distortions (e.g., Refs. [8–10]). The models used to analyze the redshift-space distortion measurements are considered within a fitting range that excludes the small-scale clustering. For instance, in Ref. [9] the authors inferred a linear growth rate using galaxy-galaxy redshift-space distortions with a cutting scale along the line of sight > 10 h⁻¹ Mpc.

While the small-scale information with the expected higher growth rate is usually disregarded due to the limitations of models, the main idea of our study is to provide a description for the growth rate on these nonlinear scales. To develop this model, we could try to reproduce the density profiles shown in Fig. 1, for instance, by using the well-known Zel'dovich approximation [31], which links the initial density profiles Δ(a ini ) to a later-time Δ( a ) assuming no shell crossing and mass conservation (e.g., Ref. [32]). These approximations as well as the spherical evolution (e.g., Refs. [33,34]) have been investigated in the literature, and recently the authors of Ref. [32] found that both Zel'dovich and spherical evolution lead to a similar evolution of an initially spherical density perturbation, which is in very good agreement with N -body simulations in some special cases [e.g., voids that are compensated, Δ( R = R$_{v}$ ) > 0, where R$_{v}$ is the radius of a void]. However, these two methods that describe the nonlinear evolution have one main disadvantage: they require as an

input the initial density perturbation Δ( R , a ini ). Because the evolution of this initial density profile becomes nonlinear at late times, a small modification in the initial input can lead to very different predictions of Δ( R , a = 1). This makes it very difficult, from an observational point of view, to precisely probe the initial densities and connect them to cosmologies, although recent developments have been made through probing projected void density profiles (e.g., Ref. [18]).

In this study, we adopt an approach that has the advantage of not requiring the initial conditions of density profiles. Instead of modeling the global nonlinear evolution of densities in terms of scales as done in perturbation theories, we generalize the nonlinear evolution of the growth rate f as a function of the local density, which is equivalent to having a model for f (Δ), where Δ( R ) is the value of the density contrast within the radius R . Our approach is referred to as log-normal Monte Carlo random walks (MCRWs) and was developed in Ref. [13]. It relies on the empirical observation that the late-time probability density function (PDF) of the galaxies (and, hence, the dark matter density fluctuations) is well described by a log-normal PDF (e.g., Refs. [35–37]). This has been confirmed by several studies using N -body simulations (e.g., Refs. [37–39]) even in the highly nonlinear regime (down to R ∼ 2 h⁻¹ Mpc for ΛCDM [40]). Using this log-normal (LN) assumption, we have generated a set of log-normal Monte Carlo random walks [13]. These walks are ensem- bles of density contrast vectors Δ$_{LN}$( R ) that are numeri- cally generated from a log-normal distribution, and aim to describe the density contrasts around random positions in the late-time Universe. The starting point of this method uses the framework of excursion set theory [41]: for Gaussian initial density perturbations, the evolution of the density contrast, smoothed on a scale R and at a random position (e.g., x = 0), is

$$
\frac{\partial \Delta(R, \mathbf{x} = 0)}{\partial R} = \int \frac{d^3 k}{2\pi^3} \tilde{\delta}_k \frac{\partial \tilde{W}(k, R)}{\partial R}, \quad (2)
$$

where ˜ δ$_{k}$ and ˜ W ( k, R ) are the Fourier transforms of the density fluctuation, and the filter function (top-hat in real space), respectively. For Gaussian initial conditions, ˜ δ$_{k}$ satisfies 〈 ˜ δ$_{k}$ ˜ δ ′ k 〉 ≡ δ$_{D}$ ( k − k $^{′}$) P$_{lin}$ ( k ), where P$_{lin}$ ( k ) is the linear matter power spectrum. For each initial realization of the density fluctuations ˜ δ$_{k}$ , the stochastic differential Eq. (2) can be solved numerically assuming that Δ( R → ∞ ) = 0 (e.g., Ref. [41]). Hence, we have a discrete set of values { Δ( R$_{1}$ ) , Δ( R$_{2}$ ) , ..Δ( R$_{N}$ ) } at each smoothing scale { R$_{1}$, R$_{2}$, ...R$_{N}$ } that is by definition one random walk. Repeating this process for a large number of initial density fluctuations allows us to generate Gaussian random walks. In order to describe the later-time nonlinear density fluctuation, we follow Ref. [13] and take the log-normal transformation of each Gaussian random walk using
$$ \Delta_{\text{LN}} + 1 = \frac{1}{\sqrt{1 + \sigma_{\text{NL}}^2(R)}} \times \exp\left(\frac{\Delta}{\sigma_{\text{lin}}(R)}\sqrt{\ln(1 + \sigma_{\text{NL}}^2(R))}\right), \quad (3) $$

with

$$ \sigma_{\text{lin}}^2(R) = \frac{1}{2\pi^2} \int P_{\text{lin}}(k) \tilde{W}^2(k, R) k^2 dk, \quad (4a) $$

$$ \sigma_{\text{NL}}^2(R) = \frac{1}{2\pi^2} \int P_{\text{NL}}(k) \tilde{W}^2(k, R) k^2 dk, \quad (4b) $$

where $P_{\text{NL}}$ is the nonlinear power spectrum. Hence, to generate these random walks, we need estimates of both $P_{\text{lin}}$ and $P_{\text{NL}}$, which we obtain using CAMB [42] with the fiducial cosmology of the DEUS N-body simulations (ΔCDM).

To compare the nonlinear growth rate obtained from the MCRWs with the one obtained from N-body simulations, we proceed as follows. We start by generating, at $a = 1$, 100 000 log-normal random walks that have “physical” properties: for the overdense regions we require that $\Delta_{LN} > \Delta$ for $\Delta_{LN} > 0$ and for the underdense regions if $\Delta_{LN} < 0$ then $\Delta < 0$. To obtain the profiles at higher redshift, we do not recompute all of the walks at different redshifts, but rather we keep the values of all of the linear trajectories at $a = 1$, $\Delta^i(a = 1)$, where $i$ is the label of one selected random walk. We can therefore directly compute $\Delta^i(a) = \Delta^i(a = 1)D^+(a)/D^+(a = 1)$ [where $D^+(a)$ is the linear growth factor at $a$] and hence $\Delta_{LN}^i(a)$ using Eq. (3). From these profiles we compute the growth rate parameter $f$ using Eq. (1) and bin the $f$ values according to $\Delta$, as we did for the simulation. Figure 2 shows the comparison of $f(\Delta)$ between the model and the simulation. Remarkably, even if the MCRW density profiles are not required to match those measured in the N-body simulation, the evolution of the nonlinear growth rates as a function of the local density matches well between the model and simulations. Note that due to the logarithmic divergence of Eq. (1), the values of $f(\Delta)$ cannot connect smoothly at $\Delta = 0$. The good agreement between our prediction and the N-body simulation measurement suggests that it is possible to extract cosmological information from these nonlinear regions.

This is a key result that shows how the nonlinear growth rate can be described by its local density. One can again draw an analogy with the island universe picture, where each region has its own growth rate depending on the mean density of the island. However, when the size of the island is small, the coupling between small and large modes becomes complex. Hence, the log-normal Monte Carlo random walks offer an alternative way to model the environmental growth rate to extract cosmological information from these nonlinear regions. Alternative methods such as those in Refs. [14–16], including the separated

universe approach [19,43], may also be helpful in improving the accuracy for the model prediction.

To summarize, we have proposed an alternative approach to extract cosmological information from the nonlinear regime. Instead of modeling “out” the nonlinear evolution of the growth rate down to a certain scale in the two-point correlation function or power spectrum, in order to recover the linear growth rate, we generalized $f$ in terms of local densities. This allowed us to map the entire spectrum of the growth rate to its underlying cosmology. We have also shown as a proof of concept that the log-normal Monte Carlo random walk approach [13] describes the function of $f(\Delta)$ reasonably well. This in principle will allow us to extract cosmological information from measurements of $f(\Delta)$.

Furthermore, because our approach goes beyond Gaussian statistics (as conventional RSD analysis uses two-point statistics), we may expect to recover more information. We expect our approach to be particularly useful for testing theories of gravity that predict a nonstandard environmental dependence for structure growth. For example, in the $F(R)$ model, due to the chameleon screening mechanism, the strength of gravity is different for different local densities [44,45]. This may alter structure formation in an environment-dependent manner, which may be better captured by measuring $f(\Delta)$. Finally, the fact that the growth rate is lower/higher in voids/clusters than its linear version indicates that one needs to employ nonlinear modeling in these low-/high-density regions (e.g., Refs. [6,10–13,46]) in order to have unbiased results.

The next question to ask is how to implement our method when analyzing data from galaxy surveys. The key is to be able to measure $f(\Delta)$ from data. We outline two possible approaches to do this. First, with the combination of a galaxy redshift survey and a lensing survey, one can use the redshift survey data to define patches of over-/underdense regions in terms of galaxy number densities $\Delta_g$ with a top-hat smoothing window. We then cross correlate these top-hat regions of different $\Delta_g$ with the lensing survey to measure their corresponding matter densities $\Delta$ and, importantly, at different tomographic bins. This is similar to measuring the lensing signal around galaxies (i.e., galaxy-galaxy lensing) except that galaxies will be replaced by top-hat regions. This will allow us to compute the numerical time derivatives of $\Delta$'s and measure $f(\Delta)$ using Eq. (1). The recent work of Refs. [17,18] has demonstrated the feasibility of this approach, where the matter densities of two-dimensional projected $\Delta_g$'s were measured using the DES survey. A challenge for this method is the requirement of having a three-dimensional galaxy redshift survey overlapping with the lensing survey of sufficient depth, which is necessary for defining $\Delta_g$ at different tomographic bins.

Second—similarly to the first approach but now with only a galaxy redshift survey—one can split the galaxy density field into top-hat regions of different $\Delta_g$'s. Different
regions would be cross correlated with the entire galaxy sample, and a RSD analysis would be performed for these $\Delta_g$-galaxy correlation functions. The $\Delta_g$-galaxy correlation is a generalized version of void-galaxy or cluster-galaxy correlations, and the latter have been demonstrated to be able to constrain the linear growth rate [5–8,10,11] with the knowledge of galaxy bias. Keeping the separated universe analogy [19], using a simple multipole decomposition and taking their ratios for those cross-correlation functions should allow us to estimate the nonlinear growth rates around those patches. One challenge for this method is the accuracy of the bias model, which is likely to be nonlinear. Also, the possible complex environmental dependence for the properties of galaxies may affect the selection of galaxies and their biases in a nontrivial way [47,48]. We

will investigate the implementation of our method in simulations and observations in more detail in future work.

## ACKNOWLEDGMENTS

We thank John Peacock, Masahiro Takada, Xin Wang, and Pengjie Zhang for useful discussions. The research leading to these results has received funding from the European Research Council under the European Community Seventh Framework Programme (FP7/2007-2013 Grant Agreement no. 279954) RC-StG EDECS. Y. C. was supported by supported by the European Research Council under Grant No. 670193. Y. C. acknowledges the support of the Royal Society through the award of a University Research Fellowship and an Enhancement Award.

[1] J. A. Peacock et al., Nature (London) **410**, 169, 2001.

[2] M. Tegmark et al., Phys. Rev. D **74**, 123507 (2006).

[3] B. A. Reid et al., Mon. Not. R. Astron. Soc. **426**, 2719 (2012).

[4] S. de la Torre et al., Astron. Astrophys. **557**, A54 (2013).

[5] Y. Zu and D. H. Weinberg, Mon. Not. R. Astron. Soc. **431**, 3319 (2013).

[6] N. Hamaus, P. M. Sutter, G. Lavaux, and B. D. Wandelt, J. Cosmol. Astropart. Phys. **11** (2015) 036.

[7] A. J. Hawken et al., Astron. Astrophys. **607**, A54 (2017).

[8] I. Achitouv, C. Blake, P. Carter, J. Koda, and F. Beutler, Phys. Rev. D **95**, 083502 (2017).

[9] F. Beutler, C. Blake, M. Colless, D. H. Jones, L. Staveley-Smith, G. B. Poole, L. Campbell, Q. Parker, W. Saunders, and F. Watson, Mon. Not. R. Astron. Soc. **423**, 3430 (2012).

[10] Y.-C. Cai, A. Taylor, J. A. Peacock, and N. Padilla, Mon. Not. R. Astron. Soc. **462**, 2465 (2016).

[11] S. Nadathur and W. J. Percival, arXiv:1712.07575.

[12] I. Achitouv, Phys. Rev. D **96**, 083506 (2017).

[13] I. Achitouv, Phys. Rev. D **94**, 103524 (2016).

[14] F. Bernardeau, S. Codis, and C. Pichon, Mon. Not. R. Astron. Soc. **449**, L105 (2015).

[15] C. Uhlemann, S. Codis, C. Pichon, F. Bernardeau, and P. Reimberg, Mon. Not. R. Astron. Soc. **460**, 1529 (2016).

[16] S. Codis, C. Pichon, F. Bernardeau, C. Uhlemann, and S. Prunet, Mon. Not. R. Astron. Soc. **460**, 1549 (2016).

[17] O. Friedrich et al., Phys. Rev. D **98**, 023508 (2018).

[18] D. Gruen et al., Phys. Rev. D **98**, 023507 (2018).

[19] C.-T. Chiang, C. Wagner, F. Schmidt, and E. Komatsu, J. Cosmol. Astropart. Phys. **05** (2014) 048.

[20] C.-T. Chiang, C. Wagner, A. G. Sánchez, F. Schmidt, and E. Komatsu, J. Cosmol. Astropart. Phys. **09** (2015) 028.

[21] I. Achitouv and C. Blake, Phys. Rev. D **92**, 083523 (2015).

[22] M. C. Neyrinck, I. Szapudi, N. McCullagh, A. S. Szalay, B. Falck, and J. Wang, Mon. Not. R. Astron. Soc. **478**, 2495 (2018).

[23] J.-M. Alimi, A. Füzfa, V. Boucher, Y. Rasera, J. Courtin, and P.-S. Corasaniti, Mon. Not. R. Astron. Soc. **401**, 775 (2010).

[24] J. Courtin, Y. Rasera, J.-M. Alimi, P.-S. Corasaniti, V. Boucher, and A. Füzfa, Mon. Not. R. Astron. Soc. **410**, 1911 (2011).

[25] Y. Rasera, J.-M. Alimi, J. Courtin, F. Roy, P.-S. Corasaniti, A. Füzfa, and V. Boucher, AIP Conf. Proc. **1241**, 1134 (2010).

[26] E. Komatsu et al., Astrophys. J. **180**, 330 (2009).

[27] R. Teyssier, Astron. Astrophys. **385**, 337 (2002).

[28] M. Davis, G. Efstathiou, C. S. Frenk, and S. D. M. White, Astrophys. J. **292**, 371 (1985).

[29] E. da Cunha et al., Pub. Astron. Soc. Aust. **34**, e047 (2017).

[30] A. Aghamousa et al. (DESI Collaboration), arXiv:1611.00036.

[31] Y. B. Zel'dovich, Astron. Astrophys. **5**, 84 (1970).

[32] P. de Fromont and J.-M. Alimi, Mon. Not. R. Astron. Soc. **473**, 5177 (2018).

[33] J. E. Gunn and J. R. Gott, III, Astrophys. J. **176**, 1 (1972).

[34] F. Bernardeau, Astrophys. J. **433**, 1 (1994).

[35] A. J. S. Hamilton, Astrophys. J. Lett. **292**, L35 (1985).

[36] F. R. Bouchet, M. A. Strauss, M. Davis, K. B. Fisher, A. Yahil, and J. P. Huchra, Astrophys. J. **417**, 36 (1993).

[37] L. Kofman, E. Bertschinger, J. M. Gelb, A. Nusser, and A. Dekel, Astrophys. J. **420**, 44 (1994).

[38] P. Coles and B. Jones, Mon. Not. R. Astron. Soc. **248**, 1 (1991).

[39] A. N. Taylor and P. I. R. Watts, Mon. Not. R. Astron. Soc. **314**, 92 (2000).

[40] I. Kayo, A. Taruya, and Y. Suto, Astrophys. J. **561**, 22 (2001).

[41] J. R. Bond, S. Cole, G. Efstathiou, and N. Kaiser, Astrophys. J. **379**, 440 (1991).

[42] A. Lewis, A. Challinor, and A. Lasenby, Astrophys. J. **538**, 473 (2000).

[43] C. Wagner, F. Schmidt, C.-T. Chiang, and E. Komatsu, Mon. Not. R. Astron. Soc. **448**, L11 (2015).
[44] J. Khoury and A. Weltman, Phys. Rev. D **69**, 044026 (2004).

[45] S. M. Carroll, A. de Felice, V. Duvvuri, D. A. Easson, M. Trodden, and M. S. Turner, Phys. Rev. D **71**, 063513 (2005).

[46] N. Hamaus, M.-C. Cousinou, A. Pisani, M. Aubert, S. Escoffier, and J. Weller, J. Cosmol. Astropart. Phys. **07** (2017) 014.

[47] B. W. O'Shea, J. H. Wise, H. Xu, and M. L. Norman, Astrophys. J. Lett. **807**, L12 (2015).

[48] E. Eardley, J. A. Peacock, T. McNaught-Roberts, C. Heymans, P. Norberg, M. Alpaslan, I. Baldry, J. Bland-Hawthorn, S. Brough, M. E. Cluver, S. P. Driver, D. J. Farrow, J. Liske, J. Loveday, and A. S. G. Robotham, Mon. Not. R. Astron. Soc. **448**, 3665 (2015).