EINSTEIN'S SIGNATURE IN COSMOLOGICAL LARGE-SCALE STRUCTURE

Einstein's general relativity provides a coherent, causal framework in which to describe classical cosmological dynamics. ΛCDM cosmology is a remarkably successful model for our observed universe, based on a spatially flat Friedmann–Lemaître–Robertson–Walker (FLRW) spacetime containing nonrelativistic, collisionless (cold) dark matter (CDM) and a cosmological constant (Λ). While the homogeneous and isotropic FLRW background can be studied analytically (using relativistic or Newtonian theory) the fully nonlinear evolution of the inhomogeneous matter distribution in ΛCDM cosmology is usually studied using Newtonian N-body simulations. These are used to make detailed predictions for comparison against large-scale galaxy surveys. As the scale and accuracy of these surveys, and hence that required from numerical simulations, continues to improve, there has been growing scrutiny of the reliability of results derived from Newtonian gravity (Wands & Slosar 2009; Chisari & Zaldarriaga 2011; Green & Wald 2012; Bruni et al. 2014b). In this Letter, we examine the characteristic signature of general relativity in the large-scale matter density and hence the galaxy distribution.

Within ΛCDM cosmology, the distribution of matter is described by a pressureless fluid. The relativistic and Newtonian descriptions of the fluid are very similar if the appropriate variables are used. The fluid is characterized by its density, ρ, and its motion, represented by kinematical variables related to its velocity at every point; Θ describes the expansion and σ describes its anisotropic deformation, or shear. We neglect vorticity to be consistent with the predictions of inflationary cosmology at early times. If we work in terms of the matter density seen by comoving observers, ρ ≡ T_μνu^μu^ν, and the expansion of the matter four-velocity, Θ ≡ ∇^μu_μ, then the relativistic and Newtonian evolution equations are formally exactly the same (Ellis 1971). The continuity equation for the matter density is

$$\begin{equation} \dot{\rho } + \Theta \rho = 0, \end{equation} \tag{ 1 }$$

while the Raychaudhuri equation for the expansion is

$$\begin{equation} \dot{\Theta } + \frac{1}{3} \Theta ^2 + 2\sigma ^2 + 4 \pi G \rho - \Lambda = 0, \end{equation} \tag{ 2 }$$

where a dot denotes derivatives with respect to the proper time of the comoving observers, corresponding to a Lagrangian time derivative in the Newtonian description.

The difference between Newton and Einstein formalisms becomes evident in the constraint equations. At the heart of Newtonian gravity is the Poisson equation,

$$\begin{equation} \nabla ^2_{\rm p} \phi = 4\pi G \rho , \end{equation} \tag{ 3 }$$

where $\nabla ^2_{\rm p}$ is the spatial Laplacian in physical coordinates. This gives a linear relation between the gravitational field, ϕ, and the matter density, ρ. In general relativity, the density and expansion are related to the intrinsic curvature of the three-dimensional space orthogonal to u^μ, denoted by ⁽³⁾R. This is the energy constraint equation from Einstein's equations:

$$\begin{equation} \frac{2}{3} \Theta ^2 - 2\sigma ^2 + {}^{(3)}R = 16\pi G \rho + 2 \Lambda . \end{equation} \tag{ 4 }$$

For the homogeneous and isotropic (σ = 0) background, we have $\rho =\bar{\rho }(t)$ and Θ = 3H(t), where H is the Hubble expansion. The evolution equations (Equations (1) and (2)) then become the familiar FLRW equations,

$$\begin{eqnarray} \dot{\bar{\rho }}=&{-}3H\bar{\rho }, \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} \dot{H}=-H^2& {-} \frac{4\pi G}{3}\bar{\rho } +\frac{\Lambda }{3}, \end{eqnarray} \tag{ 6 }$$

while the energy constraint (Equation (4)) reduces to the Friedmann constraint,

$$\begin{equation} H^2=\frac{8\pi G}{3}\bar{\rho } + \frac{\Lambda }{3}, \end{equation} \tag{ 7 }$$

with ⁽³⁾R = 0 for a spatially flat cosmology. We characterize this background model by the present-day value of the dimensionless density parameter $\Omega _{m}\equiv {8\pi G \bar{\rho }_0}/{3 H_0^2}$.

Considering inhomogeneities about the FLRW cosmology, we have

$$\begin{eqnarray} \Theta (t,x^i) &=& 3H(t) + \theta (t,x^i), \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} \rho (t,x^i) &=& \bar{\rho }(t) [ 1+\delta (t,x^i)], \end{eqnarray} \tag{ 9 }$$

and the inhomogeneous metric can be written in comoving-synchronous coordinates as

$$\begin{equation} ds^2=-dt^2 +a^2(t) e^{2\zeta (t,x^i)}\, \gamma _{kj} dx^k dx^j , \end{equation} \tag{ 10 }$$

where a(t) is the cosmological scale factor and γ_ij(t, xⁱ) has unit determinant.

The specific initial conditions for ΛCDM cosmology are set by a period of inflation in the very early universe. In particular, inflation produces an almost scale-invariant distribution for the primordial metric perturbation ζ in Equation (10) (Lyth & Liddle 2009). This allows us to consider small initial inhomogeneities on large scales, and to perform a gradient expansion (or long-wavelength approximation; Lifshitz & Khalatnikov 1963; Tomita 1975; Lyth 1985; Salopek & Bond 1990; Deruelle & Langlois 1995; Bruni & Sopuerta 2003; Rampf & Rigopoulos 2013), keeping only leading order terms, i.e., terms that are, at most, second order in spatial gradients of the metric, in particular, ζ. In this approximation, we have³ (Bruni et al. 2014a)

$$\begin{equation} \delta \sim \theta \sim \sigma \sim {}^{(3)}R \sim \nabla ^2, \end{equation} \tag{ 11 }$$

where ∇ is the spatial gradient in comoving coordinates. We emphasize that δ, θ, and ⁽³⁾R contain all orders in a conventional perturbative expansion. They are leading order quantities only in terms of spatial gradients.

With this proviso, δ and θ satisfy the simple evolution equations, from Equations (1) and (2):

$$\begin{eqnarray} \dot{\delta }+ \theta &=& {\cal O}(\nabla ^4), \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} \dot{\theta }+2H\theta +4\pi G \bar{\rho }\delta &=& {\cal O}(\nabla ^4). \end{eqnarray} \tag{ 13 }$$

These quantities are subject to the energy constraint, from Equation (4),

$$\begin{equation} \frac{{}^{(3)}R}{4} +H\theta = 4\pi G \bar{\rho }\delta + {\cal O}(\nabla ^4). \end{equation} \tag{ 14 }$$

Taking the time derivative of this equation and using the evolution equations for θ and δ, we generalize the well-known first-order result that the conformal curvature,

$$\begin{equation} R \equiv {}^{(3)}R a^2, \end{equation} \tag{ 15 }$$

remains constant in this large-scale limit (Lukash 1980; Lyth 1985; Bruni et al. 1992), i.e., R is a first integral of Equations (12) and (13).

Equations (12)–(14) are well known in perturbation theory: they are the same linear differential equations and constraints that can be derived at first order in a conventional perturbative expansion in the synchronous-comoving gauge (Bruni et al. 2014a), or in a covariant gauge-invariant fashion for corresponding quantities (Bruni et al. 1992). Their solution is, therefore, formally the same as in the first-order perturbation theory. The two independent solutions of these linear differential equations are a decaying and a growing mode (Peebles 1980). Thanks to inflation, the decaying mode is negligible; thus, we focus on the growing-mode solution. In the large-scale limit, we thus have (Bruni et al. 2014a)

$$\begin{equation} \delta = C(x^i) D_+(t) , \quad \theta = -C(x^i) \dot{D}_+(t), \end{equation} \tag{ 16 }$$

where the growth factor, D₊(t), is proportional to the scale factor, a(t), in an Einstein–de Sitter (Ω_m = 1) cosmology (Bernardeau et al. 2002). The growing-mode amplitude, C(xⁱ), is related to the conformal curvature on large scales through the energy constraint equation (Equation (14)) evaluated at an initial time t_IN early in the matter-dominated era,

$$\begin{equation} C(x) = \frac{R}{10a_{{\rm IN}}^2H_{{\rm IN}}^2 D_{+{\rm IN}}} . \end{equation} \tag{ 17 }$$

The growing-mode solution for δ and θ on large scales, Equation (16), has the same time dependence as the first-order perturbative solution; thus, it may be referred to as the linearly growing mode. However, we remark again that, in our nonlinear case, R is only conserved at leading order in our gradient expansion;⁴δ, θ, R, and C here contain the large-scale part of all perturbative orders.

In single-field, slow-roll inflation, the primordial metric perturbation, ζ(t, xⁱ), is predicted to have an almost Gaussian distribution (Maldacena 2003; Acquaviva et al. 2003). Crucially, the conformal curvature R, Equation (15), is a nonlinear function of the spatial metric in Equation (10). Considering only the scalar part of the initial metric perturbation at leading order on large scales, the spatial metric can be taken to have a simplified form,⁵γ_ij ≃ δ_ij, and the conformal curvature R is then a nonlinear function of only the perturbation ζ in Equation (10). With γ_ij ≃ δ_ij, the function a(t)exp [ζ(x)] effectively acts as a local scale factor in the so-called "separate universe" picture corresponding to our gradient expansion. R then represents the corresponding local spatial curvature and takes a beautifully simple and exact form (Wald 1984):

$$\begin{eqnarray} R &\simeq & \exp {(-2\zeta)} [ -4 \nabla ^2 \zeta - 2 \left(\nabla \zeta \right)^2 ] \nonumber \\ &\simeq & {-} 4 \left[ \nabla ^2\zeta + \frac{1}{2} \left(\nabla \zeta \right)^2 -2 \zeta \nabla ^2 \zeta\right. \nonumber\\ &&\left.-\zeta \left(\nabla \zeta \right)^2 + 2 \zeta ^2 \nabla ^2 \zeta +\cdots \, \right] . \end{eqnarray} \tag{ 18 }$$

This expression is second order in spatial gradients, consistent with Equation (11), but nonlinear in terms of the metric perturbation, ζ. Consequently, even if ζ is described by a Gaussian distribution, its nonlinear relation to the curvature R leads to a nonGaussian distribution (K. Koyama 2014, in preparation) for the comoving density contrast, δ, determined by the amplitude (Equation (17)) of the growing mode (Equation (16)).

At first order, in a perturbative expansion, we have from Equation (18)

$$\begin{equation} R_1 = -4\nabla ^2 \zeta _1 . \end{equation} \tag{ 19 }$$

Substituting this in the constraint equation (Equation (14)), we recover the Poisson equation (Equation (3)), where we identify the Newtonian potential in terms of the first-order Ricci curvature and the inhomogeneous expansion in the comoving-synchronous gauge:

$$\begin{eqnarray} \nabla ^{2} \phi _1 &=& a^2 \left[ \frac{1}{4} {}^{(3)}R_1 + H\theta _1 \right] \nonumber \\ &=& - \nabla ^{2} \zeta _1 + a^2 H \theta _1 . \end{eqnarray} \tag{ 20 }$$

Using the full nonlinear expression in Equation (18), we can write the conformal curvature in terms of ζ as an infinite series:

$$\begin{eqnarray} &&R \simeq - 4 \nabla ^2\zeta + \sum _{m = 0}^{\infty } \frac{(-2)^{m+1}}{m!} [ (m+1) (\nabla \zeta)^2 -4 \zeta \nabla ^2 \zeta] \zeta ^m . \nonumber\\ \end{eqnarray} \tag{ 21 }$$

It is this nonlinear curvature that determines the nonlinear amplitude (Equation (17)) of the growing-mode density perturbation (Equation (16)).

We wish to relate this density contrast to the observed distribution of galaxies revealed by astronomical surveys. Although a full description requires complex, nonlinear astrophysics, we can assume that in ΛCDM cosmology, galaxies form in virialized dark matter halos, which are biased tracers of the underlying matter distribution on large scales (Peacock 1999). In the simplest model of spherical collapse in Einstein–de Sitter, written in comoving-synchronous coordinates, there is an exact parametric solution (Peebles 1980),

$$\begin{eqnarray} \delta & = & \frac{9(\psi -\sin \psi)^2}{2(1-\cos \psi)^3} - 1 , \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} t & = & \frac{6^{3/2}}{2R^{3/2}} \left(\psi - \sin \psi \right) , \end{eqnarray} \tag{ 23 }$$

which can be expanded term by term as

$$\begin{equation} \delta = C D_+ + \frac{38}{21} (C D_+)^2 + \cdots \ \,, \end{equation} \tag{ 24 }$$

where the linearly growing mode (Equation (16)), with Equation (17), is given by CD₊ = Ra/10 for Einstein–de Sitter in both Newtonian theory and general relativity (Wands & Slosar 2009). Halos collapse when ψ = 2π and the linearly evolved density contrast reaches a critical value δ_* = 1.686. Thus, we can predict the number of collapsed halos (of a given mass) at a given time in terms of the number of peaks of the initial growing mode of the comoving density contrast (smoothed on a given mass scale) above a critical value (Press & Schechter 1974). Going beyond the spherical collapse, this is the barrier crossing approach, in which halos form where the linearly growing mode exceeds a critical value.

Note that it is the nonlinear amplitude, C, of the linearly evolved growing mode (Equation (16)) that determines the halo density, and this is given by the full nonlinear conformal curvature, R, in Equation (18). In a general-relativistic description of spherical collapse (Wands & Slosar 2009), it is thus the initial density contrast in the local comoving matter, δ, that predicts the distribution of halos (Bruni et al. 2012) and, as we have seen, this is nonlinearly related to the primordial metric perturbation ζ.

To understand the effect of this nonlinearity on structure formation, we consider a peak-background split (Peacock 1999), where one decomposes a field into shorter-wavelength modes that generate local peaks and much longer-wavelength modes that modulate the number density of peaks.⁶ Note that, since we have already made a gradient expansion in the above (wavenumbers k < k_max), spatial gradients of our "shorter-wavelength" modes should still be small (k_split < k < k_max), and we will now completely drop all gradients of the very long-wavelength modes (k < k_split).

For simplicity, from now on, we shall assume the simplest inflationary scenario where ζ is Gaussian, focusing on the specific general-relativistic nonGaussianity introduced by the nonlinearity of Equation (18) in the constraint (Equation (14)). We can split

$$\begin{equation} \zeta \equiv \zeta _\ell + \zeta _s , \end{equation} \tag{ 25 }$$

where the longer- and shorter-wavelength modes are independent for an initially Gaussian metric perturbation. Substituting Equation (25) into Equation (18) we obtain

$$\begin{equation} R \simeq \exp (-2\zeta _\ell) R_s + 4 \exp (-2\zeta _\ell -2\zeta _s) \nabla \zeta _\ell \nabla \zeta _s + \exp (-2\zeta _s) R_\ell , \end{equation} \tag{ 26 }$$

where

$$\begin{equation} R_s = \exp {(-2\zeta _s)} [ -4 \nabla ^2 \zeta _s - 2 (\nabla \zeta _s)^2] , \end{equation} \tag{ 27 }$$

and similarly for R_ℓ. Dropping all spatial gradients of long-wavelength modes, ∇ζ_ℓ, i.e., taking these modes to define a locally homogeneous background, we find that the spatial curvature due to short-wavelength modes is modulated such that R ≃ exp (− 2ζ_ℓ)R_s. This is consistent with the interpretation that the long-wavelength metric perturbation is a rescaling of the local background scale factor (Maldacena 2003; Creminelli & Zaldarriaga 2004; Bartolo et al. 2005b; Creminelli et al. 2011, 2013):

$$\begin{equation} a \rightarrow a_\ell = \exp (\zeta _\ell) a . \end{equation} \tag{ 28 }$$

Hence the local amplitude of the growing mode of the density contrast is also modulated (see Equations (16) and (17)):

$$\begin{equation} \delta = \exp (-2\zeta _\ell) \delta _s + \mathcal {O}(\nabla \zeta _\ell). \end{equation} \tag{ 29 }$$

The nonlinear effect of a long-wavelength overdensity, ζ_ℓ > 0, suppresses the amplitude of shorter-wavelength modes since ζ_ℓ > 0 increases the local effective scale factor, suppressing spatial curvature and thus the density contrast.

We can compare Equation (29) with local-type primordial nonGaussianity (Wands 2010) in a Newtonian approach where the amplitude of the linearly growing mode of the density is determined by the Newtonian potential,

$$\begin{eqnarray} &&\phi = \phi _1 + f_{\rm NL} \left(\phi _1^2 - \left\langle \phi _1^2 \right\rangle \right) + g_{\rm NL} \phi _1^3 + h_{\rm NL} \left(\phi _1^4 - \left\langle \phi _1^4 \right\rangle \right) +\cdots .\nonumber \\ \end{eqnarray} \tag{ 30 }$$

If we split the first-order Newtonian potential into longer- and shorter-wavelength modes ϕ₁ = ϕ_ℓ + ϕ_s and drop the gradients of ϕ_ℓ, we find

$$\begin{equation} \delta = \left(1 + 2f_{\rm NL} \phi _\ell + 3g_{\rm NL} \phi _\ell ^2 + 4h_{\rm NL} \phi _\ell ^3 +\cdots \ \right) \delta _s + \cdots. \end{equation} \tag{ 31 }$$

The modulation of the amplitude of smaller-scale density fluctuations δ_s by the long-wavelength potential, ϕ_ℓ, modifies the halo density giving rise to a strong modulation of the halo power spectrum on sufficiently large scales, where ϕ_ℓ remains finite even though the long-wavelength density contrast is suppressed, δ ∼ ∇². This leads to a scale-dependent bias in the distribution of galaxies on large scales (Dalal et al. 2008; Matarrese & Verde 2008).

If we impose the same linear relation between ζ and the Newtonian potential, ϕ = (3/5)ζ, which is valid for first-order perturbations in the matter-dominated era, then in single field, slow-roll inflation f_NL and all higher-order coefficients are suppressed. This results in an effectively Gaussian distribution for the Newtonian potential and hence (in Newtonian theory) the density field. However, expanding the exponential in Equation (29) and comparing term by term with the equivalent Newtonian expression (Equation (31)), we can identify the effective nonGaussianity on large scales in general relativity:

$$\begin{equation} f_{\rm NL}^{\rm GR}= -\frac{5}{3} ,\quad g_{\rm NL}^{\rm GR}= \frac{50}{27} , \quad h_{\rm NL}^{\rm GR}= -\frac{125}{81} , \cdots. \end{equation} \tag{ 32 }$$

More generally, we find

$$\begin{equation} f_{\rm NL}^{(n)\ {\rm GR}}= \frac{1}{n!} \left(-\frac{10}{3} \right)^{n-1} , \end{equation} \tag{ 33 }$$

where we write the local expansion (Equation (30)) as

$$\begin{equation} \phi = \phi _1 + \sum _{n=2}^\infty f_{\rm NL}^{(n)} \left(\phi _1^n - \left\langle \phi _1^n \right\rangle \right) , \end{equation} \tag{ 34 }$$

extending the previous result at second order for $f_{\rm NL}^{\rm GR}$ (Bartolo et al. 2005a, 2010; Verde & Matarrese 2009; Hidalgo et al. 2013; Bruni et al. 2014a; Uggla & Wainwright 2014) to higher orders.

Traditionally, primordial nonGaussianity is described in terms of the Newtonian gravitational potential, ϕ (for example, Equation (34)), linearly related to the density field through the Poisson equation (Equation (3)). On the other hand, inflationary predictions are expressed in terms of the primordial metric perturbation, ζ in Equation (10). Our results, valid in full nonlinearity and at large scales, show how the essential nonlinearity of general relativity produces an intrinsic nonGaussianity in the matter density field and hence the galaxy distribution on large scales, even starting from purely Gaussian primordial metric perturbations, generalizing previous results in second-order perturbation theory (Tomita 2005; Bartolo et al. 2005a, 2010; Verde & Matarrese 2009; Hidalgo et al. 2013; Bruni et al. 2014a; Uggla & Wainwright 2014). We also need a detailed modeling of observational surveys, including all geometrical and relativistic effects, to fully disentangle effects of primordial nonGaussianity from intrinsic nonlinearity in general relativity (Bruni et al. 2012; Raccanelli et al. 2014; K. Koyama 2014, in preparation). Most studies of general-relativistic effects on observations of large-scale structure have been restricted to the linear perturbation theory (Yoo et al. 2009; Bonvin & Durrer 2011; Challinor & Lewis 2011), but there have been recent attempts to include nonlinear effects; see, e.g., Thomas et al. (2014), Bertacca et al. (2014), Yoo & Zaldarriaga (2014), Di Dio et al. (2014), and Jeong & Schmidt (2014).

Alternative gravity theories may impose different constraints between the primordial metric perturbation, ζ, and the comoving density contrast, δ, and hence could, in principle, be distinguished by a different galaxy distribution on large scales. This could be an interesting approach to testing gravity on cosmological scales, complementary to existing work, which probes gravity through the growth of cosmic structure at late times (Zhao et al. 2012).

Even within the context of general relativity, the constraint equation (Equation (4)) could include additional contributions due to other fields such as dark energy/quintessence. Fields that have a negligible effect in the background could still contribute to the inhomogeneous perturbations, e.g., magnetic fields or gravitational waves. In particular, we have considered only the growing mode of scalar perturbations at early times. Tensor metric perturbations are decoupled from scalar density perturbations at first order, but do contribute to the Ricci curvature at second order, even in the large-scale limit, and hence could contribute to the nonlinear density perturbation (Matarrese et al. 1998; Dai et al. 2013), though this is expected to be sub-dominant.

In summary, in this Letter, we have obtained for the first time the fully nonlinear general-relativistic initial distribution of primordial density perturbations in ΛCDM on large scales:

$$\begin{equation} \delta = \frac{ \exp {(-2\zeta)} [ -4 \nabla ^2 \zeta - 2 (\nabla \zeta )^2 ] }{10a_{{\rm IN}}^2H_{{\rm IN}}^2 D_{+{\rm IN}}} D_+(t), \end{equation} \tag{ 35 }$$

an expression including all perturbative orders. This fully nonlinear relation between δ and ζ clearly shows that, even for a Gaussian-distributed ζ, the corresponding matter density field is nonGaussian. Assuming simple inflationary Gaussian initial conditions in ζ, and using a peak-background split, we have derived the corresponding specific general-relativistic effective nonGaussianity parameters, Equations (32) and (33), that results when a Newtonian treatment is used, i.e., a Poisson equation as the relation between δ and the Newtonian potential ϕ, and a linear relation is assumed between ζ and ϕ. Although Newtonian simulations are commonly used to study the nonlinear evolution of the matter density contrast δ, a relativistic approach is essential to properly determine the initial conditions set by a period of inflation in the very early universe. Thus, we have shown how Einstein's gravity imprints a characteristic signature in the large-scale structure of our universe.

The authors are grateful to Rob Crittenden, Roy Maartens, and Gianmassimo Tasinato for useful discussions. This work was supported by STFC grants ST/K00090X/1 and ST/L005573/1, and by PAPIIT-UNAM grants IN103413-3 and IA101414-1.