TESTING GRAVITY USING THE GROWTH OF LARGE-SCALE STRUCTURE IN THE UNIVERSE

In the framework of GR, the growth of a density fluctuation, $\delta \equiv (\rho (x,t)-\bar{\rho }_{\rm m})/\bar{\rho }_{\rm m}$, where $\bar{\rho }_{ \rm m}$ is the average matter density, depends only on the expansion history, H(a). Using the perturbed equations of motion, within GR, the growth of perturbations follows:

$$\begin{eqnarray} &&\ddot{\delta } +2H\dot{\delta } -4\pi G_N \rho _{\rm m} \delta = 0, \end{eqnarray} \tag{ 1 }$$

where G_N is the present gravitational constant found in laboratory experiments and a dot denotes a time derivative. The growth rate is f ≡ dlnδ/dlna, where δ(a) is the growing mode solution to Equation (1). Changing variables to g ≡ δ/a and allowing the gravitational constant to vary in time, denoted by $\tilde{G}$, gives (Linder 2005)

$$\begin{eqnarray} \frac{{d}^2g}{{d}a^2} &+& \left(5 + \frac{1}{2} \frac{{d\rm { ln}}H^2}{{d\rm { ln}} a}\right)\frac{1}{a}\frac{{d}g}{{d}a} \\ \nonumber &+& \left(3 + \frac{1}{2}\frac{{d\rm { ln}}H^2}{{d\rm { ln}} a} -\frac{3}{2}\frac{\tilde{G}(a)}{G_N}\Omega _{\rm {m}}(a) \right)g = 0 \,, \end{eqnarray} \tag{ 2 }$$

where $\Omega _{\rm {m}}(a)$ is the matter density parameter. Equation (2) shows that in GR, $\tilde{G}(a)/G_N =1$ and the growth of perturbations depends only on the expansion history, H(a). In modified gravity theories, however, the growth of perturbations depends on both H(a) and $\tilde{G}(a)$.

Modifications to GR provide an alternative explanation to dark energy for the accelerating cosmic expansion. Modified gravity theories can generally be divided into models which introduce a new scalar degree of freedom to Einstein's equations, e.g., scalar tensor or f(R) theories, and those which change dimensionality of space, e.g., braneworld gravity. In many such models, the time variation of fundamental constants, such as Newton's gravitational constant, G_N, is naturally present.

Self-consistent scalar tensor theories are viable alternatives to GR and give rise to an accelerating expansion at late epochs. We refer to these as "extended quintessence"   models. Calculations which follow spatial variations in the scalar field have shown that, in practice, a broad range of these models can be effectively described with a time-varying Newton's constant (Pettorino & Baccigalupi 2008; Li et al. 2010).

The variation of G_N is constrained by various observations, such as the lifespan of stars (Teller 1948), the age of globular clusters (degl'Innocenti et al. 1996), the mass of neutron stars (Thorsett 1996), and the synthesis of light nuclei (Umezu et al. 2005; Clifton et al. 2005). A time-varying G_N would also modify the temperature fluctuations in the cosmic microwave background (CMB), shifting the peaks to larger (smaller) scales on increasing (decreasing) G_N. This leads to a constraint on the variation of G, $\dot{G}/G = (-9.6\mbox{--}8.1) \times 10^{-12}$ yr⁻¹ (Chan & Chu 2007).

Here, we consider a simple model for $\tilde{G}$ (Zahn & Zaldarriaga 2003; Umezu et al. 2005; Chan & Chu 2007):

$$\begin{eqnarray} \tilde{G} &=& \mu ^2 G_N \,, \end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray} \mu ^2 = \left\lbrace \begin{array}{ll}\mu _0^2 & \mbox{if} \, \, a<a_* \\ 1 - \frac{a_s-a}{a_s-a_*}(1-\mu _0^2) &\mbox{if} \, \, a_*\le a \le a_s \\ 1 & \mbox{if} \, \, a>a_s \,. \end{array} \right. \end{eqnarray} \tag{ 4 }$$

This parameterization describes a smoothly varying $\tilde{G}$ which converges slowly to its present value, G_N, and is more physical than those based on step functions (e.g., Cui et al. 2010). The parameter, a_*, denotes the scale factor of photon decoupling, and the parameters μ₀ and a_s quantify the deviation of $\tilde{G}$ from the laboratory measured value, G_N, and the scale factor at which $\tilde{G}$ and G_N are equal, respectively. The background evolution is given by

$$\begin{eqnarray} &&H^2 = H^2_0 \frac{\tilde{G}}{G_N} \left(\frac{ \Omega _{\rm m}}{a^3} + \Omega _{{\rm DE}} e^{3\int _a^1 {d\rm { ln}} a^{\prime } [1+w(a^{\prime })]}\right). \end{eqnarray} \tag{ 5 }$$

Note that we assume an equation of state w(a) = −1 in the modified gravity model to match ΛCDM. In Equation (5), Ω_DE is the ratio of the dark energy density to the critical density today. In the left panel of Figure 1, we plot the ratio of the Hubble rate for two different cosmological models with varying $\tilde{G}$, to the Hubble rate in a ΛCDM cosmology as a function of redshift. We chose to simulate the model with the maximum deviation of $\tilde{G}$ from G_N which is still compatible with CMB measurements and solar system constraints ($\tilde{G} \rightarrow G$ as a → 1), which occurs for a stabilization redshift corresponding to a_s = 1(i.e., the green dot-dashed line in Figure 1).

Zoom In Zoom Out Reset image size

We use large volume N-body simulations to carry out the first direct test of the idea that dark energy and modified gravity cosmologies which, by construction, have exactly the same expansion history, can be distinguished by a measurement of the rate at which structure grows. The modified gravity model we simulate has the maximum deviation from Newton's constant that is compatible with observational constraints, as discussed above. We construct a quintessence model by fitting the expansion history to match the varying $\tilde{G}$ model within 0.25% over 0 ⩽ z ⩽ 200. This model is consistent with constraints on dynamical dark energy (Komatsu et al. 2009; Sánchez et al. 2009).

The simulations were carried out using a memory efficient version of the TreePM code Gadget-2, called L-Gadget-2 (Springel 2005). The simulation used N = 1024³ ∼ 1 × 10⁹ particles in a box of comoving length 1500 h⁻¹ Mpc. The comoving softening length was = 50 h⁻¹ kpc and the present day linear rms fluctuation in spheres of radius 8 h⁻¹ Mpc is σ₈ = 0.8. Simulations of extended quintessence cosmologies need to account for both the gravitational correction due to a varying $\tilde{G}$ in the Poisson equation and a modified expansion history (see Pettorino & Baccigalupi 2008). In the modified gravity simulation, both the long- and short-range TreePM algorithm force computations are modified to include a time-dependent gravitational constant. In both the dark energy and modified gravity simulations the Hubble parameter computed by the code was also changed as in Jennings et al. (2010a).

The linear theory power spectrum used to generate the initial conditions was obtained using CAMB (Lewis & Bridle 2002). We adopt a ΛCDM linear theory power spectrum at z = 0 and use consistent linear growth factors in each cosmology to obtain the power spectrum amplitude at z = 200. In principle, as the quintessence cosmology could be classed as an early dark energy model, the linear theory spectrum should be modified in shape. However, as we have shown, such a change has a negligible impact on the nonlinear spectrum and on the ratio of the quadrupole to monopole moments (Jennings et al. 2010a, 2010b).

To obtain errors on our measurements we ran 10 lower resolution simulations with 512³ particles, also in a box of comoving length 1500 h⁻¹ Mpc, with different realizations of the density field. The power spectrum was computed using the cloud in cell assignment scheme and performing a fast Fourier transform. For the initial conditions the linear growth rate for each model and ΛCDM was obtained by solving Equation (2) numerically and is plotted in the right-hand panel of Figure 1 as a function of redshift. For all the models we used the following cosmological parameters: Ω_m = 0.26, $\Omega _{\rm {DE}}=0.74$, Ω_b = 0.044, h₀ = H₀/100 km s⁻¹ Mpc⁻¹ = 0.715, and a spectral index of n_s = 0.96 (Sánchez et al. 2009). We have verified that our modifications to Gadget-2 are accurate by checking that the growth of the fundamental mode in the simulations agrees with the linear theory predictions.

The matter power spectrum in redshift space can be decomposed into multipole moments using Legendre polynomials. The ratio of the quadrupole and monopole moments of the matter power spectrum is plotted in Figure 2. We model redshift space distortions in the distant observer approximation by perturbing the particle positions down one of the Cartesian axes, using the suitably scaled component of the peculiar velocity. The simulation results show that this ratio has a strong dependence on wavenumber. This can be contrasted with the linear perturbation theory prediction (Cole et al. 1994),

$$\begin{eqnarray} \frac{P_2(k)}{P_0(k)} &=& \frac{4\beta /3 +4\beta ^2/7}{1 + 2\beta /3 + \beta ^2/5}, \end{eqnarray} \tag{ 6 }$$

where β = f/b and b is the linear bias, which is unity for dark matter; Equation (6) is independent of scale (horizontal lines in Figure 2). We note that, by considering redshift space distortions in the clustering of the dark matter, we are testing theoretical models against the simplest possible case. The distortions will inevitably be more complicated for dark matter halos and galaxies, for which the bias factor b can have scale dependence (e.g., Angulo et al. 2008). In Figure 2, the quadrupole to monopole ratio increases in amplitude with redshift, due to the evolution in the matter density parameter. At z = 0 there is a 2.5% difference between the linear theory growth rates in the two models. However, at this level, the measured ratios P₂/P₀ are indistinguishable on the very largest scales k < 0.02 h Mpc⁻¹(green dotted and blue dashed horizontal lines). At z = 0.5 and z = 1, the linear theory predictions for the growth rates in the two models differ by 4% and 6%, respectively. The error on this ratio measured from the lower resolution simulations is shown by the shaded region in Figure 2.

Zoom In Zoom Out Reset image size

In addition to the linear theory model we consider two variants. The first is the Gaussian model (Peacock & Dodds 1994),

$$\begin{eqnarray} &&P^s(k,\mu) = P^r(k) (1+\beta \mu ^2)^2 e^{(-k^2 \mu ^2 \sigma _p^2)}, \end{eqnarray} \tag{ 7 }$$

where σ_p is the pairwise velocity dispersion along the line of sight, which is treated as a parameter to be fitted. We refer to Equation (7) as the "linear theory plus damping" model. The damping introduces a scale dependence into the ratio P₂/P₀. The second variant model takes into account departures from linear theory as well as including small-scale damping (Scoccimarro 2004):

$$\begin{eqnarray} P^s(k,\mu)= &&(P_{\delta \delta }(k) + 2 f\mu ^2 P_{\delta \theta }(k) \\ && + f^2\mu ^4P_{\theta \theta }(k)) \times e^{-(f k \mu \sigma _v)^2}, \nonumber \end{eqnarray} \tag{ 8 }$$

where σ_v is the one-dimensional linear velocity dispersion, and P_θθ and P_δθ are the velocity divergence auto and cross-power spectrum, respectively, measured from the simulations (see also Jennings et al. 2010b). We refer to Equation (8) as the "quasi-linear plus damping" model. We note that P₂/P₀ is more sensitive to changes in f than the ratio of the monopole moment of the redshift space to real space P(k), see Figure 3, and, as a result, the 1σ errors for f are smaller when fitting to P₂/P₀.

Zoom In Zoom Out Reset image size

We now apply the above models to the simulation results. In Figure 4, we plot the measured ratio P₂/P₀, for the modified gravity simulation at z = 0.5, together with the theoretical predictions. In the left panel, the correct value of f for this cosmology together with the best-fit value for σ_p and σ_v in the range 0.01 ⩽ k( h Mpc⁻¹) ⩽ 0.25 was used in the linear theory plus damping and quasi-linear plus damping models, respectively. In the right panel, the best-fit value for f obtained by fitting over the same range of wavenumbers is used for all models plotted. The value of f obtained for the linear theory model is sensitive to the maximum value of k used in the fit. It is clear that both the linear theory and the linear theory plus damping models fail to predict the correct value for f, with the best-fitting values differing by ∼40% and ∼6%, respectively, from the true value, see Figure 5. All the models plotted in the right panel in Figure 4 use the value of f recovered when k_max = 0.25 h Mpc⁻¹. The quasi-linear plus damping model recovers the correct value of f over this wavenumber range to within ∼0.64%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To test these models for the redshift space power spectrum further we vary the maximum wavenumber, k_max, used in the fit and plot the recovered growth rate as a function of k_max in Figure 5. With an accurate model we would recover the correct value for the growth rate f and the answer would be independent of the value of k_max adopted, with the only sensitivity to k_max being in the error on the growth rate. Figure 5 shows that the quasi-linear plus damping model comes closest to meeting this ideal. Even this model breaks down beyond k_max ∼ 0.3h Mpc⁻¹, which suggests that the modeling of the small-scale velocity dispersion can be improved. Most importantly, this model recovers the correct value for f and can distinguish between the two cosmologies. The models based on linear theory perform less well. In fact, the answer depends strongly on the maximum wavenumber used in the fit. In Figure 5, filled symbols are plotted for scales over which the model is a good description of the measured ratio (i.e., χ²/ν ∼ 1, where ν is the number of degrees of freedom).