Perturbations of the Lambda-CDM model in a dynamical systems perspective

TL;DR

This paper introduces dynamical systems methods to analyze linear perturbations in the $\Lambda$CDM cosmological model, providing pedagogical insights, studying asymptotic behaviors, and proposing a more accurate approximation for the growth rate of structures.

Contribution

It applies dynamical systems techniques to $\Lambda$CDM perturbations, offering new analytical tools, asymptotic analysis, and an improved growth rate approximation.

Findings

Dynamical systems effectively describe $\Lambda$CDM perturbations.

New approximation for the growth rate $f(z)$ is more accurate.

Analysis of shear and Weyl tensors' asymptotic properties.

Abstract

The observational success and simplicity of the $\Lambda$CDM model, and the explicit analytic perturbations thereof, set the standard for any alternative cosmology. It therefore serves as a comparison ground and as a test case for methods which can be extended and applied to other cosmological models. In this paper we introduce dynamical systems and methods to describe linear scalar and tensor perturbations of the $\Lambda$CDM model, which serve as pedagogical examples that show the global illustrative powers of dynamical systems in the context of cosmological perturbations. We also study the asymptotic properties of the shear and Weyl tensors and discuss the validity of the perturbations as approximations to the Einstein field equations. Furthermore, we give a new approximation for the linear growth rate, $f(z) = \frac{d\ln \delta}{d\ln a} = \Omega^{\frac{6}{11}}_m - \frac{1}{70}(1-\Omega_m)^{\frac{5}{2}}$, where $z$ is the cosmological redshift, $\Omega_m=\Omega_m(z)$, while $a$ is the background scale factor, and show that it is much more accurate than the previous ones in the literature.

Equations

$$d{s}^2=-d{t}^2+a^2{\gamma }_{ij}d{x}^{i}d{x}^j=a^2\left(-d{\eta }^2+{\gamma }_{ij}d{x}^{i}d{x}^j\right)=-H^{-2}d{N}^2+a^2{\gamma }_{ij}d{x}^{i}d{x}^j,$$

$$N=\ln x,x=\frac{a}{a_0},$$

$$x=\frac{1}{1+z},$$

$${\mathbf{D}}^2\to -k^2,$$

$$u_1^{\prime }$$

$$u_2^{\prime }$$

$$\begin{array}{rl}{\theta }^{\prime }& =d{\cos }^2\theta +(e-b)\sin \theta \cos \theta -c{\sin }^2\theta \\ & =\frac{1}{2}[(d-c)+(d+c)\cos 2\theta +(e-b)\sin 2\theta ].\end{array}$$

$$y(N,k^2)=y:=\tan \theta =u_2(N,k^2)/u_1(N,k^2),$$

$$y^{\prime }=d+(e-b)y-c{y}^2,$$

$$T^{\prime }=G(T),\text{where}G(T)=F^{\prime }(F^{-1}(T)).$$

$$T=\frac{h(x)}{1+h(x)},$$

$${\theta }^{\prime }$$

$$T^{\prime }$$

$$y^{\prime }$$

$$T^{\prime }$$

$$d{s}^2=a^2\left(-(1+2\varphi )d{\eta }^2+{\mathbf{D}}_{i}Bd\eta d{x}^i+[(1-2\psi ){\gamma }_{ij}+2{\mathbf{D}}_i{\mathbf{D}}_{j}C]d{x}^{i}d{x}^j\right),$$

$${\delta }_m=\frac{{}^{(1)}{\rho }_m}{{}^{(0)}{\rho }_m},$$

$${}^{(1)}{u}_i=a{\mathbf{D}}_{i}V.$$

$$(1+q){\partial }_N((1+q{)}^{-1}{\varphi }_{\mathrm{c}})$$

$${\partial }_N(a^2{B}_{\mathrm{c}})$$

$$q^{\prime }=-(1+q)(1+3c_s^2-2q).$$

$$1+q=\frac{3}{2}(1+w).$$

$$c_s^2=0,1+w={\Omega }_m.$$

$${\varphi }_{\mathrm{c}}^{\prime }$$

$$(\mathcal{H}{B}_{\mathrm{c}}{)}^{\prime }$$

$${\Omega }_m={\Omega }_{m0}{x}^{-1}{\left(\frac{{\mathcal{H}}_0}{\mathcal{H}}\right)}^2,{\Omega }_{\Lambda }={\Omega }_{\Lambda 0}{x}^2{\left(\frac{{\mathcal{H}}_0}{\mathcal{H}}\right)}^2,$$

$$\frac{{\Omega }_{\Lambda }}{{\Omega }_m}={\lambda }_m{x}^3,\text{where}{\lambda }_m:=\frac{{\Omega }_{\Lambda 0}}{{\Omega }_{m0}}=\frac{\Lambda }{{\rho }_{m0}}.$$

$${\Omega }_m=\frac{1}{1+{\lambda }_m{x}^3}.$$

$${\varphi }_{\mathrm{c}}^{\prime }$$

$$(\mathcal{H}{B}_{\mathrm{c}}{)}^{\prime }$$

$$y_{\mathrm{c}}=\frac{-{\varphi }_{\mathrm{c}}}{\mathcal{H}{B}_{\mathrm{c}}}.$$

$$y_{\mathrm{c}}^{\prime }=y_{\mathrm{c}}\left(\frac{5}{2}-\frac{9}{2}(1-{\Omega }_m)-y_{\mathrm{c}}\right).$$

$$1-{\Omega }_m=\frac{{\lambda }_m{x}^3}{1+{\lambda }_m{x}^3}={\Omega }_{\Lambda },$$

$$T=\frac{{\lambda }_m{x}^3}{1+{\lambda }_m{x}^3}=1-\frac{1}{1+{\lambda }_m{x}^3}={\Omega }_{\Lambda }.$$

Full text

1 Introduction

This paper is the first in a series of papers dealing with cosmological perturbations by means of dynamical systems formulations and methods. We will show how a wide variety of increasingly complex problems can be formulated as dynamical systems and how powerful dynamical systems methods can be applied to yield insights about cosmological perturbations. In particular, this will make it possible to apply approximation techniques from the theory of dynamical systems and to obtain illustrative pictures as well as mathematically rigorous results about the global structure of the various models solution spaces, thereby also providing a context for especially physically interesting solutions. Our aim is thus to provide a useful complement to traditional approaches to cosmological perturbations. Step by step we will introduce increasingly sophisticated models and methods. In this paper we focus on first order scalar and tensor perturbations of the spatially flat $\Lambda$CDM model with dust and a positive cosmological constant $\Lambda$. Since this model is mathematically simple and is compatible with a wide range of observations it is a natural choice as a first example to illustrate the most simple aspects of our dynamical systems approach to perturbative cosmology. Moreover, due to their observational success they provide a comparative testing ground for any observational contender.

Dynamical systems have been used before to analyze cosmological linear scalar perturbations in general relativity (GR) with an open Robertson-Walker (RW) geometry as background, see [24] and references therein. These models turn out to be somewhat easier to handle than models with a spatially flat background,111In contrast to the spatially flat RW models, ${\cal H}^{-2}=(aH)^{-2}$ is always bounded for the spatially open models if the energy density is non-negative. This is due to how the spatial curvature term appears in the background Gauss constraint, which also can be used to solve for ${\cal H}^{-2}$. Since ${\cal H}^{-2}$ appears in connection with the spatial derivatives of the perturbed equations, this significantly simplifies the analysis of the perturbed field equations in the open case. which is what we will focus on. However, there has also been some previous work on dynamical systems in this area in GR, notably [4] which treated dust and radiation as a single fluid.

The main foundation for standard cosmology is the spatially flat RW geometry, characterized by a line element that can be written as

[TABLE]

where $a(t)$ is the background scale factor, $\gamma_{ij}$ is the flat spatial 3-metric, which in Cartesian coordinates is given by $\delta_{ij}$, and $H=a^{-1}da/dt$ is the background Hubble variable. The different time coordinates above are the clock time $t$, the conformal time $\eta$, and the $e$-fold time

[TABLE]

where $N$ describes the number of background $e$-foldings with respect to some reference epoch $a=a_{0}$ at which $x=1$ (a negative $N$ describes the number of $e$-folds before the reference time). If this reference epoch is the present time, then

[TABLE]

where $z$ is the cosmological background redshift. Much of the work in cosmological perturbation theory uses the clock time $t$ or the conformal time $\eta$, but we find it more convenient to use the $e$-fold time $N$ as the starting point for our work.

The outline of the paper is as follows. In the next section we describe some aspects of the framework of our new dynamical systems approach to cosmological perturbations. In Section 3 we consider linear scalar perturbations of the $\Lambda$CDM models from a dynamical systems perspective. We also present and discuss previously known analytic results in the present context. In addition we describe and discuss various asymptotic approximation methods and give a new more accurate approximation for the so-called linear growth rate. In Section 4 we consider the linear tensor perturbations from a dynamical systems perspective. In Section 5 we relate the present work to the full state space of GR and discuss the validity of the perturbations as approximations to the Einstein field equations. This is done by comparing asymptotic perturbative results with asymptotic results in GR for the Hubble-normalized comoving shear and Weyl tensor. The paper is concluded with some final remarks in Section 6.

2 Dynamical systems approach to cosmological perturbations

Our aim is to analyze cosmological perturbations by formulating the governing equations as regular dynamical systems on compact state spaces. In this paper we will consider first order scalar and tensor perturbations of the spatially flat $\Lambda$CDM models with dust and a positive cosmological constant $\Lambda$. However, since similar methods apply to models whose matter content consists of a perfect fluid with a barotropic equation of state and with, or without, a cosmological constant it is useful to point out some common features these models exhibit. The governing equation for scalar perturbations for this class of models can be a single second order partial differential equation (e.g. the Bardeen equation) or a system of two coupled first order partial differential equations (e.g. the Kodama-Sasaki equations), depending on the choice of gauge and the choice of variables, see e.g. [21] and references therein. For an arbitrary equation of state these governing equations will contain the spatial Laplacian ${\bf D}^{2}$. This is also the case for the tensor perturbations which obey a linear second order differential master equation for a single variable.

In order to obtain ordinary differential equations (ODEs) we make a spatial Fourier decomposition of the perturbation variables, which involves replacing the perturbation variables by their Fourier coefficients and making the transition

[TABLE]

where $k$ is the wave number. At this stage if we have a second order ODE as a governing equation we would replace it by a system of two coupled first order ODEs by using the first order time derivative as an independent variable. In this way the Einstein field equations for first order scalar and tensor perturbations of the models under consideration can be written as a system of linear ODEs of the form

[TABLE]

for two linear perturbation variables $u_{1}=u_{1}(N,k^{2})$ and $u_{2}=u_{2}(N,k^{2})$, where a ′ denotes differentiation with respect to the $e$-fold time $N$. Due to the spatial Fourier decomposition one obtains complex variables, but thanks to that the differential equations are linear, real and complex parts satisfy the same equations, and therefore we can without loss of generality consider $u_{1}$ and $u_{2}$ to be real, as are $b(N,k^{2}),c(N,k^{2}),d(N,k^{2}),e(N,k^{2})$, whose specific form is determined by the background model. This is exemplified in detail in section 4 where tensor perturbations are treated.

The next step toward obtaining a regular dynamical system on a bounded state space is to introduce polar coordinates $u_{1}=r\cos\theta$, $u_{2}=r\sin\theta$, which thanks to the linearity of the equations lead to a decoupling of an equation for $r$ and a reduction of the system of two linear ODEs (5) to one nonlinear ODE for $\theta(N,k^{2})$:

[TABLE]

However, locally it is more convenient to replace $\theta$ with

[TABLE]

which results in a Riccati ODE given by

[TABLE]

which describes the essential “reduced” dynamics of the problem since, e.g., $u_{1}$ can be obtained as a quadrature from the decoupled equation $u_{1}^{\prime}=u_{1}(b+cy)$ once the equation for $y$ has been solved.

To obtain a dynamical system, i.e. a system of first order autonomous ODEs, that incorporates the dynamics described by the non-autonomous ODE (8) or (6), we introduce a new dependent variable $T=F(N)$, where $F$ is a non-negative, bounded and increasing function. It follows that $T$ satisfies an autonomous ODE of the form

[TABLE]

Note that $G(T)$ is obtained by differentiating $F(N)$ and then expressing $N$ in terms $T$ using the inverse function $N=F^{-1}(T)$. This equation is adjoined to (8) or (6), thereby yielding a 2-dimensional dynamical system for $(y,T)$ or $(\theta,T)$.

In practice, however, we have found it convenient to use the relation $N=\ln(x)$ to express $T=F(N)$ as a function of $x$ rather than $N$. In particular we write $T$ in terms of a subsidiary function $h(x)$ according to

[TABLE]

where $h(x)$ is a non-negative, increasing, explicitly invertible, and suitably differentiable function which satisfies $h(0)=0$ and $h(x)\rightarrow\infty$ when $x\rightarrow\infty$. This ensures that $T$ is a non-negative, bounded and increasing function that satisfies $T(0)=0$ and $\lim_{x\rightarrow\infty}T(x)=1$. To find the function $G(T)$ in (9), differentiate (10) with respect to $N$ using $dx/dN=x$ and then express $x$ in terms of $T$ using (10).

Augmenting the ODE (6) for $\theta$ with the ODE (9) for $T$ yields the dynamical system

[TABLE]

where the functions $b,c,d,e$ are expressed as functions of $T$. The state space for this system is a finite cylinder in which all orbits (i.e., solution trajectories), begin at the boundary $T=0$ and end at the boundary $T=1$. We thus extend the state space to include these boundaries, which we refer to as the (extended) compactified state space cylinder $[0,1]\times S^{1}$. Note that $T$ can be regarded as a bounded time variable for which $T=0$ and $T=1$ correspond to $x=0$ and $x\rightarrow\infty$, respectively.

As mentioned, for local analysis it is more convenient to use $y$ instead of $\theta$, which results in the system

[TABLE]

In this representation traversing the infinite strip defined by $-\infty<y<\infty$ and $0\leq T\leq 1$ twice corresponds to making one revolution on the state space cylinder $[0,1]\times S^{1}$ with the global coordinates $T$ and $\theta$.222The reason for traversing the strip twice is that $y=u_{2}/u_{1}=-u_{2}/(-u_{1})$, i.e., the mapping is two-to-one, although note that the right hand side of (6) has a periodicity $\pi$.

The key to implementing the above procedure is to make an appropriate choice for the function $h(x)$ that determines the new dependent variable $T$ through equation (10). In a particular application the choice is motivated by the form of the coefficient functions $b,c,d,e$ in the initial system (5). We will use this procedure in the following sections to describe the scalar and tensor perturbations of the $\Lambda$CDM model.

3 Dynamical systems approach to scalar perturbations of $\Lambda$CDM

Linear scalar perturbations of a spatially flat RW background geometry can be described by

[TABLE]

where $\phi,B,\psi,C$ are the metric scalar perturbation variables (see e.g. Uggla and Wainwright (2018) [21].) The scalar matter perturbations for $\Lambda$CDM are given by the first order fractional matter density perturbation

[TABLE]

where the superscripts denote the order of the perturbation. The scalar velocity perturbation ${V}$ is defined by the first order perturbation of the spatial covariant 4-velocity components according to the relation

[TABLE]

Setting $C=0$ fixes the spatial gauge and covers most of the familiar gauges, although it excludes the synchronous gauge (for a recent work using the synchronous gauge, see e.g. [8], and also [23] and [21] for further discussions and references). The temporal gauge can be fixed in a number of ways, e.g., by setting to zero one of the variables $B$, $\psi$, $V$, ${\delta_{m}}$. We use the following terminology and subscripts to label the gauges:

• i)

Poisson (Newtonian, longitudinal) gauge, subscript p, defined by $B_{\mathrm{p}}=0$,

• ii)

uniform curvature gauge, subscript c, defined by $\psi_{\mathrm{c}}=0$,

• iii)

total matter (comoving) gauge, subscript v, defined by $V_{\mathrm{v}}=0$.

As shown in e.g. [21] each choice leads to a different system of governing equations. In particular, for linear perturbations the uniform (flat) curvature gauge leads to a simple system of two first order partial differential equations, which form a natural starting point for a dynamical systems analysis. We obtain this system of equations from [21], where it is given in the following form:333See equations (10), (54a) and (54b) in [21]. We have dropped superscripts (1) on the perturbation variables, and have set $\Gamma=0$ since we are considering a barotropic fluid.

[TABLE]

with

[TABLE]

For the $\Lambda$CDM model we have444See, for example, [21], Appendix B.

[TABLE]

As basic perturbation variables for the $\Lambda$CDM model we choose $(u_{1},u_{2})=({\cal H}B_{\mathrm{c}},-\phi_{\mathrm{c}})$, where we have scaled $B_{\mathrm{c}}$ with $\cal H$ to obtain a dimensionless quantity, as discussed in [21] and [23]. We now specialize the governing equations (16) to the $\Lambda$CDM model, using equations (17) and (18). After expanding the $\partial_{N}$ derivatives and replacing $\partial_{N}$ by ′ we obtain

[TABLE]

Since the Laplacian ${\bf D}^{2}$ has dropped out, there is no need for a spatial Fourier decomposition. The system (19), which is of the form (5), constitutes the starting point for transforming the governing equations into a dynamical system.

We now digress to describe the background dynamics of the $\Lambda$CDM model in order to obtain the $x$-dependence of $\Omega_{m}$. The density parameters are given by

[TABLE]

which implies that

[TABLE]

Since $\Omega_{m}+\Omega_{\Lambda}=1$ it follows from (21) that

[TABLE]

3.1 Derivation of the dynamical systems

The metric perturbations in the uniform curvature gauge

We have shown that perturbations of the $\Lambda$CDM model are described by the system of ODEs (19) which we repeat here:

[TABLE]

with variables $(u_{1},u_{2})=({\cal H}B_{\mathrm{c}},-\phi_{\mathrm{c}})$. This problem is explicitly solvable, and we will later give the solution. Here, however, we will use it as a first illustration of our dynamical systems approach. We therefore introduce $y=u_{2}/u_{1}$, as in (7), but give $y$ a subscript, $y_{\mathrm{c}}$, as a reminder that we are using the uniform curvature gauge, thus

[TABLE]

It follows from (23) that

[TABLE]

To complete the process of constructing a dynamical system we now have to make a choice for the additional independent variable $T$. It follows from equation (22) that

[TABLE]

which suggests that we choose555The second expression is useful for the calculation that follows.

[TABLE]

In other words the function $h(x)$ in (10) is given by $h(x)=\lambda_{m}x^{3}$, which has the desired properties. Differentiating (27) with respect to $N$ yields

[TABLE]

after expressing $x$ in terms of $T$ using (27) again. This is the desired autonomous ODE for $T$. Note that $0\leq T\leq 1$ and that $T(x)$ is an increasing function, as required. We finally use that $1-\Omega_{m}=T$ in equation (25). The resulting equation, with (28), comprises the desired (analytic) dynamical system as follows:

[TABLE]

Since $y_{\mathrm{c}}$ is defined by the ratio of two metric coefficients, one orbit of the above system corresponds to a one-parameter family of solutions. Furthermore, the solutions for the dynamical system (29) correspond to solutions with all possible (non-zero) values of $\lambda_{m}$ since $\lambda_{m}$ is incorporated in the definition of $T$. In order to solve for $y_{\mathrm{c}}(N,x^{i})$ one has to impose an initial condition at $N=0$, which corresponds to the initial epoch $a=a_{0}$, of the form $y_{\mathrm{c}}(0,x^{i})=f(x^{i}),$ where $f(x^{i})$ is an arbitrary spatial function. Note that the initial value of $T(N)$, where $T(N)$ was chosen in order to arrive at a simple form for (29) and is given by (27), is $T(0)=1-\Omega_{m0}$ where $\Omega_{m0}=\Omega_{m}(0)$ is the initial value of $\Omega_{m}(N)$.

Finally we note that the full dynamical system in the uniform curvature gauge consists of the state space vector $({\cal H}B_{\mathrm{c}},\phi_{\mathrm{c}},T)$ governed by the equations (23) and (29b). This state space can then be covered by using polar coordinates ${\cal H}B_{\mathrm{c}}=r_{\mathrm{c}}\cos\theta_{\mathrm{c}},-\phi_{\mathrm{c}}=r_{\mathrm{c}}\sin\theta_{\mathrm{c}}$, which leads to a decoupling of $r_{\mathrm{c}}$ from a dynamical system with a reduced state space vector $(\theta_{\mathrm{c}},T)$. Locally it is, however, more convenient to use $y_{\mathrm{c}}(N,x^{i})$ and ${\cal H}B_{\mathrm{c}}(N,x^{i})$ as the decoupled variable, which obeys

[TABLE]

This equation is easily solved as a quadrature once a solution $y_{\mathrm{c}}(N,x^{i}),T(N)$ has been obtained, thereby also yielding the second constant of the motion which together with the constant of the motion for the reduced system for $(y_{\mathrm{c}},T)$ characterizes the various scalar perturbations.

The fractional density perturbation in the total matter gauge

The density contrast $\delta$ is defined by

[TABLE]

where ${}^{(0)}\!\rho_{m}$ is the background matter density and ${}^{(1)}\!\rho_{m}$ is the linear density perturbation. The comoving density contrast $\delta_{\mathrm{v}}$, the gauge invariant that equals $\delta$ in the the total matter/comoving gauge, plays a central role in observational cosmology. Since $\delta_{\mathrm{v}}$ satisfies a second order evolution equation we can apply our method to analyze its behaviour from a dynamical systems perspective. In the $\Lambda$CDM universe this evolution equation has the following form when using $e$-fold time $N$:666This evolution equation has been given in a general context in Uggla and Wainwright (2018) [21], equation (71a.b). When specialized to $\Lambda$CDM and to first order perturbations, equation (71a) reads ${\cal L}_{D}\delta_{\mathrm{v}}=0$, where the differential operator ${\cal L}_{D}$ in (71b) reduces to ${\cal L}_{D}=\partial_{N}^{2}+(1-q)\partial_{N}-(1+q).$ Since $1+q={\textstyle{3\over 2}}(1+w)={\textstyle{3\over 2}}\Omega_{m}$ we obtain our equation (31). To avoid confusion we note that the usual density contrast $\delta$ is defined by normalizing with ${}^{(0)}\!\rho_{m}$, while $\delta$ in [21] is defined by normalizing with ${}^{(0)}\!\rho+{}^{(0)}\!p$. In the $\Lambda$CDM universe these two normalizations are the same (see [21], Appendix B).

[TABLE]

In order to formulate this differential equation as a dynamical system we use $(u_{1},u_{2})=(\delta_{\mathrm{v}}^{\prime},\delta_{\mathrm{v}})$ as basic variables, and define

[TABLE]

in accordance with equation (7).

We obtain the ODE for $y_{\mathrm{v}}$ by differentiating (33) using (32), and we make the same choice (27) of $T$ as before. This results in the following dynamical system

[TABLE]

This formulation of the evolution equations which is based on $y_{\mathrm{v}}=\delta_{\mathrm{v}}^{\prime}/\delta_{\mathrm{v}}$ differs from the system that is based on $y_{\mathrm{c}}=-\phi_{\mathrm{c}}/{\cal H}{B}_{\mathrm{c}}$. However, $y_{\mathrm{c}}$ and $y_{\mathrm{v}}$ are closely related, as will be shown later.

3.2 Dynamical systems analysis

In this section we use dynamical systems methods to obtain qualitative information about the family of all solutions of the perturbed field equations, including asymptotic descriptions of the solutions at early and late times.

The uniform curvature gauge

When using the uniform curvature gauge the differential equations that define the reduced dynamical system are given by equations (29), which we here repeat for the reader’s convenience:

[TABLE]

The associated state space is the infinite strip:

[TABLE]

with boundaries $T=0$ and $T=1$. Alternatively we can use the angular variable $\theta_{\mathrm{c}}$ defined by $y_{\mathrm{c}}=\tan\theta_{\mathrm{c}}$, which yields the reduced regular global state space which is a finite section of a cylinder $[0,1]\times S^{1}$, as discussed in section 2.

The analysis of this dynamical system is straightforward. Equation (35b) shows that all the fixed points of the system lie on the boundaries $T=0$ and $T=1$, and that each orbit is past asymptotic to a fixed point on $T=0$ and future asymptotic to a fixed point on $T=1$. The local stability of the fixed points can be determined in the usual way by linearizing the differential equations.

When specialized to the $T=0$ and $T=1$ boundaries the differential equation (35a) results in the following equations:

[TABLE]

It follows that on the $T=0$ boundary there are two fixed points (four in the case of $\theta_{\mathrm{c}}$, although they are connected by the discrete symmetry $\theta\rightarrow\theta+\pi$ and thereby given by the two fixed points for $y_{\mathrm{c}}$):777The first subscript on a fixed point $\mathrm{P}$ denotes the value of $T$ while the + (-) denotes the larger (smaller) value of $y_{\mathrm{c}}$.

[TABLE]

Linearization shows that the fixed point $\mathrm{P_{0+}}$ is a hyperbolic saddle and that a single orbit originates from $\mathrm{P_{0+}}$. To obtain an analytical approximation for this special orbit, we make a series expansion for $y_{\mathrm{c}}$ in powers of $T$ and use the system (29) to solve for the coefficients,888This can be mathematically justified be relating this series expansion to a so-called Picard expansion, see e.g. [17] and references therein. which leads to

[TABLE]

for $T$ and $x$ close to zero. Next, linearization shows that the fixed point $\mathrm{P_{0-}}$ is a hyperbolic source, and hence a one-parameter family of orbits originates from $\mathrm{P_{0-}}$. A series expansion in $T$ results in

[TABLE]

where $C_{0}$ is a spatial function that parameterizes the orbits, depending on the spatial position.

Inserting the series expansion (40) into equation (30) for ${\cal H}B_{\mathrm{c}}$ shows that ${\cal H}B_{\mathrm{c}}\rightarrow\infty$ toward the past for all orbits that originate from $\mathrm{P_{0-}}$ (to leading order it suffices to insert $y_{c}=0$ into (30)). In contrast inserting (39) into (30) shows that ${\cal H}B_{\mathrm{c}}$ is finite asymptotically for the orbit that originates from $\mathrm{P_{0+}}$. Since ${\cal H}B_{\mathrm{c}}=-\psi_{\mathrm{p}}$ (see equation (61) below), it follows that the orbit that originates from $\mathrm{P_{0+}}$ is the only orbit along which the perturbations remain finite into the past i.e. as $T\rightarrow 0$. Below we will identify this orbit with the so-called growing mode solution toward the future. Along all other orbits into the past, i.e., solutions that originate from $\mathrm{P_{0-}}$, (some of) the perturbations increase without bound and hence will not asymptotically approximate solutions of the full Einstein equations. This is related to the fact that a general solution of the perturbed Einstein equations for scalar $\Lambda$CDM perturbations is a linear combination of a so-called growing mode and a decaying mode, where the latter has the property that it becomes unbounded into the past.

On the $T=1$ boundary it follows from (37b) that the fixed points are given by

[TABLE]

Linearization shows that the fixed point $\mathrm{P_{1-}}$ is a hyperbolic saddle that attracts a single orbit that is past asymptotic to $T=0$. A series expansion in $1-T$ yields:

[TABLE]

Finally, linearization shows that the fixed point $\mathrm{P_{1+}}$ is a hyperbolic sink, which implies that a one-parameter family of orbits is asymptotic to $\mathrm{P_{1+}}$. A series expansion in $1-T$ results in

[TABLE]

where $C_{1}$ is a spatial function that parameterizes the orbits, depending on the spatial position.

The preceding analysis of the stability of the fixed points enables one to predict the qualitative form of the state space (36) in figure 1 which shows the fixed points, the special heteroclinic orbits (i.e. solution trajectories that originate and end at two distinct fixed points), $\mathrm{P_{0-}}\rightarrow\mathrm{P_{1-}}$, $\mathrm{P_{0-}}\rightarrow\mathrm{P_{1+}}$ (note that equation (35a) shows that this latter orbit is the invariant set $y_{\mathrm{c}}=0$, which in turn corresponds to the invariant set $\phi_{\mathrm{c}}=0$ of equation (23a)), and $\mathrm{P_{0+}}\rightarrow\mathrm{P_{1+}}$, together with some typical orbits of the dynamical system (35). The most significant aspect of the state space in figure 1 is the heteroclinic orbit $\mathrm{P_{0+}}\rightarrow\mathrm{P_{1+}}$, the growing mode solution, which represents the most physically important solution of the perturbation equations (actually a one parameter family of solutions), and which thereby is the primary focus in cosmological perturbation theory. Note also that due to its observational success, any observational contender must presumably result in a solution trajectory in $y_{\mathrm{c}}$ and $T$ that is quite similar to that of the growing mode solution for the observational redshift range.

As a final side remark we note the following about the heteroclinic growing mode orbit $\mathrm{P_{0+}}\rightarrow\mathrm{P_{1+}}$. Since $\mathrm{P_{0+}}$ is a saddle point and $\mathrm{P_{1+}}$ is a local sink, this orbit acts as an “attractor solution” toward the future of nearby orbits, reminiscent of “attractor solutions” in inflationary cosmology. The latter, however, are associated with a non-hyperbolic saddle point, whose unstable manifold is a center manifold, rather than a hyperbolic saddle point, which allows the inflationary regime to last longer than it would if the saddle was hyperbolic, see e.g. [2] for a discussion of attractor solutions in inflationary cosmology.

The comoving density perturbation

The dynamical system based on the comoving density contrast $y_{\mathrm{v}}=\delta_{\mathrm{v}}^{\prime}/\delta_{\mathrm{v}}$ is given by (34), which we repeat here:

[TABLE]

The structure of the orbits of this system is very similar to that of the system (35), with the fixed points lying on the boundaries $T=0$ and $T=1$.

When specialized to these boundaries the differential equation (44a) results in the following equations:

[TABLE]

It follows that on the $T=0$ ($x=0$) boundary the fixed points are given by

[TABLE]

Linearization shows that the fixed point $\mathrm{P_{0-}}$ is a hyperbolic source, while the hyperbolic saddle $\mathrm{P}_{0+}$ has a single orbit originating from it into the interior, which, as we will see, describes the growing mode solution. A series expansion for this orbit yields

[TABLE]

On the $T=1$ boundary the fixed points are given by:

[TABLE]

The fixed point $\mathrm{P_{1+}}$ is a hyperbolic sink, while the fixed point $\mathrm{P_{1-}}$ is a hyperbolic saddle which attracts a single interior orbit.

There is a one-to-one correspondence as regards fixed points and stability properties between the dynamical system (44) and the dynamical system (35) for the state space $(y_{\mathrm{c}},T)$. This correspondence is reflected in the similarity between the form of the orbits in figure (1a) and figure (2a).999Just before submitting the present paper, Basilakos et al. published a paper [3] on the archive with a diagram that corresponds to figure (2a), but with $\Omega_{m}$ instead of $T=1-\Omega_{m}$. They also gave the explicit solution (their equation (31)) for $y_{\mathrm{v}}$, which in their notation was called $U_{m}$, but instead of using the parameters $C_{\pm}$ in the next subsection they used $U_{m0}$ and $\Omega_{m0}$, which are related to $C_{\pm}$ according to $U_{m0}=\frac{C_{+}T_{0}^{\frac{5}{6}}(1-T_{0})^{\frac{2}{3}}}{C_{+}I_{0}+C_{-}}-\frac{3}{2}(1-T_{0})$ where $T_{0}=1-\Omega_{m0}$ and $I_{0}=\frac{1}{3}\int_{0}^{T_{0}}T^{-\frac{1}{6}}(1-T)^{-\frac{1}{3}}dT$. Our parametrization is adapted to the solution structure so that e.g. the physically important growing mode solution is easily obtained by setting $C_{-}=0$. In particular the heteroclinic orbit $P_{0+}\rightarrow P_{1+}$ again represents solutions that only contain the growing mode. Moreover, since $y_{\mathrm{v}}$ is positive on this orbit and $y_{\mathrm{v}}=\delta_{\mathrm{v}}^{\prime}/\delta_{\mathrm{v}}$ it follows that when $\delta_{\mathrm{v}}>0$ then $\delta_{\mathrm{v}}$ is growing throughout its evolution (hence the name, growing mode), although $\delta_{\mathrm{v}}$ approaches a constant value toward the future since $y_{\mathrm{v}}=0$ at $\mathrm{P_{1+}}$.101010An analogue of the orbit $y_{\mathrm{c}}=0$ in figure (1a) also occurs in figure (2a). It is the heteroclinic orbit $P_{0-}\rightarrow P_{1+}$ and it is given by $y_{\mathrm{v}}=-(3/2)(1-T)$. Finally, the comoving density contrast $y_{\mathrm{v}}$ for the growing mode solution is often referred to as the linear growth rate, which we will denote as $f(z)$ when expressed in terms of the redshift $z$.

In the above discussion of $y_{\mathrm{v}}$ we have refrained from giving details about asymptotics, except for the growing mode solution. The reason for this is that $y_{\mathrm{v}}$ and $y_{\mathrm{c}}$ can be related to each other, as shown in the next section about explicit solutions. This will also establish the identification of the growing mode orbits in the $y_{\mathrm{c}}-T$ and $y_{\mathrm{v}}-T$ state spaces.

3.3 Explicit solutions

Here we derive and discuss the explicit solutions for $y_{\mathrm{c}}$ and $y_{\mathrm{v}}$ as functions of the time variable $T$. For the variable $y_{\mathrm{c}}=-\phi_{\mathrm{c}}/({\cal H}B_{\mathrm{c}})$ we begin with the governing equations (16) in the uniform curvature gauge, which when applied to the $\Lambda$CDM universe simplify to

[TABLE]

where $\partial_{a}$ refers to the partial derivative with respect to the background scale factor $a$. We can successively solve the equations to obtain

[TABLE]

where $C_{+},\,C_{-}$ are arbitrary spatial functions and where the perturbation evolution function $g(a)$ is given by111111See section 7 in [22] for properties and a discussion about the perturbation evolution function $g(a)$.

[TABLE]

In order to obtain $y_{\mathrm{c}}$ as a function of $T$ we need to express $\phi_{\mathrm{c}}$ and ${\cal H}B_{\mathrm{c}}$ as functions of $T$, using equations (20)-(22), (27) and (29b). The results are as follows:121212As intermediate steps we obtain ${\cal H}/a^{2}=CT^{-5/6}(1-T)^{-1/3}$, $1+q={\textstyle{3\over 2}}(1-T)$, and

$da/a=dT/(3T(1-T))$. Here $C$ is a constant that does not appear in the final result. The constants $C_{\pm}$ in (50) have been redefined in obtaining (52).

[TABLE]

where131313$I(T)$ is related to $g$, when expressed in $T$, according to $g(T)=\frac{3}{2}T^{-\frac{5}{6}}(1-T)^{\frac{1}{3}}I(T)$.

[TABLE]

It follows that the variable $y_{\mathrm{c}}=-\phi_{\mathrm{c}}/({\cal H}B_{\mathrm{c}})$ is given by

[TABLE]

The function $I(T)$ is well-defined at the end points,

[TABLE]

and the leading order behaviour of $I(T)$ is given by the following limits:

[TABLE]

We can now relate the explicit solution (54) to the orbits of the dynamical system in figure 1. It follows from (54) that there are two special values of $C_{-}$, namely $C_{-}=0$ and $C_{-}=-C_{+}I_{1}$ that affect the limit of $y_{\mathrm{c}}$ as $T\rightarrow 0,\,1$, as follows:

[TABLE]

Referring to figure 1 we conclude that the orbit with $C_{-}=0$ is the growing mode orbit, i.e. the orbit that is past asymptotic to the fixed point $P_{0+}$, while the orbits with $C_{-}\neq 0$ are those that are past asymptotic to the local source $P_{0-}$. Further, the orbit with $C_{-}=-C_{+}I_{1}$ is the special orbit that is future asymptotic to the fixed point $P_{1-}$, while the orbits with $C_{-}\neq-C_{+}I_{1}$ are those that are future asymptotic to the local sink $P_{1+}$.

In section 3.2 we derived the leading terms of a series expansions of $y_{\mathrm{c}}$ for each of the above four classes of orbits. We now derive a full series for $y_{\mathrm{c}}$ as given by (54), by giving a series expansion for $I(T)$, first in powers of $T$ and then in powers of $1-T$. It follows from (53) that (see e.g. [1])

[TABLE]

where ${}_{2}F_{1}$ is the Gaussian hypergeometric function, and $(p)_{n}$ the Pochhammer symbol, with $(p)_{0}=1$, and $(p)_{n}=p(p+1)...(p+n-1)$, $n\in\mathbb{N}$. A truncated version of (58) when substituted in (54) with $C_{-}=0$ yields141414Making this transition is made somewhat complicated by the fact the series for $I(T)$ is in the denominator of $y_{\mathrm{c}}$. Truncate the series after three terms if $C_{-}=0$ and after one term otherwise. the leading term expression (39), while if $C_{-}\neq 0$ we obtain the leading term expression (40) with the arbitrary function $C_{0}$ given by $C_{0}=C_{+}/C_{-}.$

Similarly, we can expand $I(T)$ in powers of $1-T$ obtaining (see e.g. [1])

[TABLE]

A truncated version of (59) in (54) with $C_{+}I_{1}+C_{-}=0$ yields the expression (42), while if $C_{+}I_{1}+C_{-}\neq 0$ we obtain the expression (43) with the arbitrary function $C_{1}$ given by $C_{1}=C_{+}/(C_{+}I_{1}+C_{-}).$

We can also use the solution for $y_{\mathrm{c}}$ to obtain an explicit expression for $y_{\mathrm{v}}$, as follows. The GR Poisson equation and the conservation of momentum equation are given by

[TABLE]

respectively, see [21]. In addition we have the relations151515The first two are standard change of gauge formulas, see e.g. [23], section 3, and the third is the velocity constraint in the uniform curvature gauge, see  [21], equation (54c).

[TABLE]

which lead to the following result:

[TABLE]

We substitute the explicit solution (52) into (62) to obtain

[TABLE]

which for the growing mode ($C_{-}=0$) reduces to

[TABLE]

It follows from (54) with $C_{-}=0$ that for this special orbit we have the simple relation

[TABLE]

This relation then establishes that the heteroclinic orbit $P_{0+}\rightarrow P_{1+}$ in both systems represents the same solution, i.e., the growing mode solution. Making a Fourier decomposition of ${\cal H}B_{\mathrm{c}}$ and $\phi_{c}$ results in that equation (62) also reduces to (65). This explains why making the above variable transformation transforms the system (34) to the system (35).

We conclude by remarking that the explicit solutions makes it possible to give an analytic example of how the growing mode solution acts as an “attractor solution” by defining

[TABLE]

as follows from equation (54) where $y_{\mathrm{c}}[G]$ stands for the growing mode solution with $C_{-}=0$. This is an analytic description of the deviation of solutions from the growing mode solution, which can be given in terms of the redshift $z$ since $T=T(z)$. Finally, note that it is possible for the spatial function $C_{-}$ to be zero at one, two or three spatial coordinates. This gives rise to a so-called permanent spike, which is an asymptotic spatial discontinuity on a surface, line, or point, respectively, a situation that can be compared with the features of the special explicit solutions of the Einstein field equations given in [13, 12].

3.4 Approximation methods

We now turn to approximation techniques and make comparisons with the exact results. The motivation for this is to develop increasingly accurate approximation schemes that work for problems that are not explicitly solvable. From comparisons with the explicit solution we find that the various series expansions obtained by dynamical systems methods give correct asymptotic approximations for the solution. In the present context it is of particular interest to obtain, preferably globally, or at least for all observable redshifts, accurate approximations for the growing mode solution, especially with methods that can be applied to other problems and dynamical systems.

We can improve the accuracy of the previous truncated series approximations by using them to derive Padé approximants.161616For a discussion of Padé approximants in a cosmological setting and further references, see e.g. [2]. For example, for the growing mode solution originating from $\mathrm{P}_{0+}$ in the $y_{\mathrm{v}}-T$ state space we get the following Padé approximation:

[TABLE]

A completely different global approximation for the growing mode solution can be obtained by observing that it is a slowly varying heteroclinic orbit in $T=\Omega_{\Lambda}$, or, equivalently $\Omega_{m}$, as can be seen from figure 2. More precisely, it is a trajectory that bends slightly to the right in figure 2 with respect to the straight line $y_{\mathrm{v}}=1-T$ that goes through the fixed points $\mathrm{P}_{0+}$ and $\mathrm{P}_{1+}$. This motivates the following variable transformation from $y_{\mathrm{v}}$ to a new function $\gamma(T)$, defined by

[TABLE]

where $\gamma(T)$ is called the growth (rate) index function (for a historical background, see the discussion below). Expressing equation (34) as a first order differential equation for $\gamma$ results in

[TABLE]

where a power series expansion of $\gamma(T)$ in $T$ yields

[TABLE]

Background models with different $\Omega_{m0}$ and $\Omega_{\Lambda 0}$, i.e., different $\lambda_{m}=\Omega_{\Lambda 0}/\Omega_{m0}$, all have the same trajectories in the dynamical systems pictures. However, when e.g. $y_{\mathrm{v}}$ is plotted against the redshift $z$, a given solution in the dynamical systems picture, e.g. the growing mode solution, results in a one-parameter set of solutions parameterized by $\lambda_{m}$ since

[TABLE]

The growing mode solution together with the $y_{\mathrm{v}}=[1,1]_{y_{\mathrm{v}}}$ Padé approximant and the approximation $y_{\mathrm{v}}=(1-T)^{\frac{6}{11}}=\Omega_{m}^{\frac{6}{11}}$ are depicted in a redshift diagram in figure 4 for $\lambda_{m}=7/3$, i.e., $\Omega_{m0}=0.3$, $\Omega_{\Lambda 0}=0.7$ (chosen for simplicity and in agreement with recent observational data, see e.g. [6]), and $\lambda_{m}=1/3$, i.e., $\Omega_{m0}=0.75$, $\Omega_{\Lambda 0}=0.25$.

Historically the type of analytical approximation given in equation (70) can be traced back to Peebles [15], who gave it in the form $y_{\mathrm{v}0}=\Omega^{0.6}_{m0}$ as an approximation at the present time. For models containing negative curvature $\Omega_{k}$, a similar approximation was given by Lightman and Schester [11], with $y_{\mathrm{v}0}=\Omega^{\frac{4}{7}}_{m0}$, see also [5]. Lahav et al. [10] realized that this type of approximations could be extended to general redshift $z$, and they further refined it according to

[TABLE]

where

[TABLE]

Later Wang and Steinhardt [25], see also references therein, clarified that the values $0.6$ and $\frac{4}{7}$ are approximations to $\frac{6}{11}$, obtained by the series expansion given in (70); see also [16] for the correct next order term in the exponent for $\Lambda$CDM, which is given in equation (70). The approximations in equation (72), $y_{\mathrm{v}}=(1-T)^{\frac{6}{11}}=\Omega_{m}^{\frac{6}{11}}$, and $y_{\mathrm{v}}=\Omega_{m}^{\frac{6}{11}+\frac{3\cdot 5}{11^{2}\cdot 17}(1-\Omega_{m})}$ are quite good when compared with the explicit solution (63), as shown by plotting the errors $\Delta=y_{\mathrm{v}}[\mathrm{Approx.}]-y_{\mathrm{v}}[\mathrm{Explicit}]$, i.e. the difference between an approximation and the explicit solution (63), in figure 5(a). For another discussion about approximations, see [9]. Let us now introduce the following simple correction to $y_{\mathrm{v}}=\Omega^{\frac{6}{11}}_{m}$, which compensates for the errors for the intermediate evolution,171717Most of the approximations in this subsection are algorithmic in nature, but $y_{\mathrm{v}}=\Omega_{m}^{\frac{4}{7}}$, $y_{\mathrm{v}}=\Omega_{m}^{0.6}+\frac{1}{70}\!\left(1-\frac{1}{2}\Omega_{m}(1+\Omega_{m})\right)$ and our new expression $y_{\mathrm{v}}=\Omega_{m}^{\frac{6}{11}}-\frac{1}{70}(1-\Omega_{m})^{\frac{5}{2}}$ are not. In the dynamical systems setting they correspond to curve fitting, but all aim at approximating the full analytical growing mode solution of $y_{\mathrm{v}}$.

[TABLE]

where the range of $z$ is determined by $T<T_{0}=\Omega_{\Lambda 0}\approx 0.7$. As seen in figure 5(b), equation (74) is an approximation which is more accurate than $y_{\mathrm{v}}=\Omega_{m}^{\frac{6}{11}+\frac{3\cdot 5}{11^{2}\cdot 17}(1-\Omega_{m})}$ by several orders of magnitude for all observational $z$. It is possible to improve the accuracy further, but we have not been able to do so significantly with an approximation that is as simple or simpler than the present one. It should be noted that all approximations are quite good for large $z$, i.e., small $T$; the differences reside in values for $z$ that are relevant for large scale structure formation at comparatively late times.

4 Tensor perturbations

Linear tensor perturbations of the spatially flat RW background geometry are characterized by a perturbed metric of the form

[TABLE]

where $h_{ij}$ is a gauge invariant that satisfies $\gamma^{ij}h_{ij}=0$, ${\bf D}^{i}h_{ij}=0$. In the absence of anisotropic stresses the perturbed Einstein equations assume the following form (e.g. Malik and Wands (2009) [14], equation (8.6)):

[TABLE]

This equation is now going to be used to illustrate the general discussion in section 2 in more detail. In order to formulate this equation as a dynamical system we first introduce $e$-fold time, which yields

[TABLE]

Using spatial Cartesian coordinates we apply the Fourier transform to this partial differential equation and write the transform of $h_{ij}$ as a linear combination of time-independent polarization tensors $e_{ij}^{+}$ and $e_{ij}^{\times}$:181818See, e.g., Weinberg (2008) [26], page 232 for more detail.

[TABLE]

where $k$ is the wave number, and $h^{+},h^{\times}$ are complex-valued functions. The outcome is that each of the functions $h^{+},h^{\times}$ satisfy the following ordinary differential equation:

[TABLE]

In the rest of this section the complex-valued function $h(N,k^{2})$ will denote either $h^{+}$ or $h^{\times}$ and we will usually not indicate the dependence on the wave number $k$ explicitly. We finally specialize this differential equation to the $\Lambda$CDM universe using (17b), (18) and, (20) to obtain

[TABLE]

where we have introduced a scaled wave number according to $k^{2}_{0}=k^{2}/(\lambda_{m}^{1/3}\Omega_{m0}{\cal H}_{0}^{2})$.

The final step in constructing a dynamical system (a system of first order autonomous differential equations) is to follow the approach used for the density perturbation and introduce a (real-valued) quotient variable analogous to $y_{\mathrm{v}}$. However, since $h$ in (80) is complex and stands for one of the two functions $h^{+},h^{\times}$, we must write $h^{+}=h^{+}_{1}+ih^{+}_{2}$, and $h^{\times}=h^{\times}_{1}+ih^{\times}_{2}$ where $h^{+}_{1},h^{+}_{2},h^{\times}_{1},h^{\times}_{2}$ are four real-valued functions which independently satisfy (80). Then for any one of these four functions, which we simply denote by $h=h(N,k^{2})$, we define

[TABLE]

To complete the process we have to choose a suitable time function $T$. Considerations of the temporally $x$-dependent functions in the system of perturbative ODEs in order to obtain a regular dynamical system result in a different $T$ than for scalar perturbations, namely191919Compare with (10) and (27).

[TABLE]

We now calculate $y_{\mathrm{t}}^{\prime}$ by differentiating (81) and using (80) and $T^{\prime}$ by differentiating (82). After expressing the coefficients in terms of $T$ using (20) and (82) we obtain

[TABLE]

Thus equation (83) describes a one-parameter family of real-valued analytic dynamical systems labelled by the parameter $k_{0}^{2}$, which yields the long wavelength approximation when $k_{0}^{2}=0$. The state space is again the infinite strip defined by $-\infty<y_{\mathrm{t}}<\infty,\,0\leq T\leq 1$, which is to be traversed twice in order to describe the global state space of $\theta_{\mathrm{t}}$ and $T$.

The structure of the orbits of this system is very similar to that of the system (35), with $T$ a monotonically increasing function and with the fixed points lying on the boundaries $T=0$ and $T=1$. We note that the fixed points do not depend on the arbitrary parameter $k_{0}^{2}$, since on the boundaries the function $F(T)$ equals zero. When specialized to these boundaries the differential equation (44a) results in the following equations:

[TABLE]

It follows that on the $T=0$ boundary the fixed points are given by

[TABLE]

Linearization shows that the fixed point $\mathrm{P_{0+}}$ is a hyperbolic saddle and that a single orbit originates from $\mathrm{P_{0+}}$ and is future asymptotic to the $T=1$ boundary. This special orbit is approximated by

[TABLE]

Next, linearization shows that the fixed point $\mathrm{P_{0-}}$ is a hyperbolic source, and hence a one-parameter family of orbits originates from $\mathrm{P_{0-}}$. A series expansion in $T$ results in

[TABLE]

where $C_{0}$ parameterizes the different orbits.

On the $T=1$ boundary the fixed points are as follows:

[TABLE]

Linearization shows that the fixed point $\mathrm{P_{1-}}$ is a hyperbolic saddle that attracts a single orbit that is past asymptotic to $T=0$. A series expansion in $1-T$ yields:

[TABLE]

Finally, linearization shows that the fixed point $\mathrm{P_{1+}}$ is a hyperbolic sink, which implies that a one-parameter family of orbits is asymptotic to $\mathrm{P_{1+}}$. A series expansion in $1-T$ results in

[TABLE]

where $C_{1}$ parameterizes the different orbits.

The orbit structure for tensor perturbations for $\Lambda$CDM cosmology for a variety of values of $k_{0}$ is illustrated in figure 6. Note that for large $k_{0}$, i.e., $k_{0}\gg 0$, orbits start to circulate the state space at an intermediate stage of the evolution. This regime is approximately described by the short wavelength limit for dust, but eventually the cosmological constant starts to dominate and the future asymptotic limits for the orbits $y_{\mathrm{t}}(N,k_{0}^{2})$ are the same for all values of $k_{0}$, since the fixed points at $T=1$ are independent of $k_{0}$.

The long wavelength limit corresponds to $k_{0}=0$ and yields an explicit solution for $y_{\mathrm{t}}$, which is most conveniently expressed by using

[TABLE]

which results in

[TABLE]

In this case there are two special orbits, one with $C_{+}=0$ going from $\mathrm{P}_{0+}$ to $\mathrm{P}_{1+}$ and one with $C_{-}=-C_{+}$, for which $y_{\mathrm{t}}=-(1+\Omega_{\Lambda}^{\frac{1}{2}})$, going from $\mathrm{P}_{0-}$ to $\mathrm{P}_{1-}$. The structure of the orbits in the long wavelength limit is illustrated in figure 7.

5 Measures of anisotropy: the shear and Weyl tensors

In this paper we have given a global analysis of the evolution of linear scalar and tensor perturbations of a $\Lambda$CDM universe by formulating the perturbed Einstein equations as dynamical systems. In order to place this analysis in perspective we now consider the family of solutions of the Einstein equations whose matter content is dust and a cosmological constant. These solutions model a class of universes that generalize the $\Lambda$CDM universe which is the unique model in this class that describes a completely isotropic universe with flat spatial geometry. We will thus refer to these universes as generalized $\Lambda$CDM universes. Complete isotropy is characterized by the requirement that the shear tensor $\sigma^{a}\!_{b}$ of the fluid congruence and the Weyl curvature tensor $C^{ab}\!_{cd}$ are zero, see e.g. [24] and references therein. From an observational point of view one is interested in universes in which the anisotropy is small, by which one means that the shear tensor and the Weyl tensor are small relative to the overall expansion of the universe, described by the Hubble scalar $H$. We thus form dimensionless scalars by normalizing the contracted shear and Weyl tensors202020One can decompose the Weyl tensor into an electric and magnetic part relative to a given timelike congruence, with spatial components denoted $C^{\eta i}\!_{\eta j}$ and $C^{\eta i}\!_{lm}\varepsilon^{lm}\!_{j}$, where $\varepsilon_{ijk}$ is the three dimensional alternating symbol. One can form space-time scalars using $C^{ab}\!_{cd}$ as below, or spatial scalars using the electric and magnetic parts. with an appropriate power of the Hubble scalar $H$:

[TABLE]

There exists a number of results about the evolution of generalized $\Lambda$CDM universes. First, as regards late times it has been shown that universes in this family that expand indefinitely approach the de Sitter model for an open set of initial conditions,212121This class of cosmologies is labelled by eight arbitrary spatial functions, Lim et al. (2004) [13], page 8. in the sense that

[TABLE]

Second it has been shown that there is a subset of models which approximate the flat RW (Friedmann-Lemaître) model on approach to the singularity,222222This class of cosmologies is labelled by three arbitrary spatial functions [13], page 11.

[TABLE]

The singularity in these models is referred to as an isotropic singularity (see Goode and Wainwright (1985) [7]). On the other hand a typical model undergoes a more complicated evolution described by BKL oscillations and possible so-called spike oscillations on approach to the singularity (see e.g. Uggla (2013) [18]). Nevertheless, the Hubble-normalized anisotropy scalars $\Sigma^{2}$ and $C^{2}$ remain bounded during this process, a result that we will use later in this section.

In this paper we have illustrated the global solution space of the linearly perturbed $\Lambda$CDM models using dynamical systems that describe scalar and tensor perturbations separately. We will now use the shear and Weyl tensors to compare the asymptotic behaviour of the perturbations at early and late times with the full state space picture of Einstein’s field equations in the Hubble-normalized state space description given in [13], and briefly described above. The purpose with this is to shed light on the important issue of assessing the validity of cosmological linear pertubations as approximations to solutions of the Einstein field equations.

We refer to [19] for expressions for the perturbed shear and Weyl tensors.232323See equations (B35a) and (B41c), which we specialize as follows. We assume zero anisotropic stress, which implies $\Psi=\Phi$, a flat background ($K=0$) and for simplicity we exclude the vector mode (${\bf B}_{i}=0$). Note that the present $h_{ij}={\textstyle{1\over 2}}{\bf C}_{ij}$. For the linear scalar and tensor perturbations the perturbed shear tensor is given by

[TABLE]

5.1 Scalar perturbations

Using the relations (61)

[TABLE]

and the scalar mode extraction operator ${\cal S}^{ij}$ we form the Hubble-normalized shear and electric Weyl scalars for the scalar perturbations:252525${\cal S}^{ij}={\textstyle{3\over 2}}({\bf D}^{-2})^{2}{\bf D}^{ij}$, where ${\bf D}^{-2}$ is the inverse spatial Laplacian, ${\bf D}_{ij}:={\bf D}_{(i}{\bf D}_{j)}-{\textstyle{1\over 3}}\delta_{ij}{\bf D}^{2}$, and where spatial indices are raised with $\delta^{ij}$, see Appendix A in [20].

[TABLE]

In the dynamical systems formulation for the scalar perturbations we have chosen a variable $T$ according to equation (27) which yields $1+q=\frac{3}{2}(1-T)$ and

[TABLE]

One can express $\Sigma_{s}$ and $W_{s}$ in the non-reduced state space $({\cal H}B_{\mathrm{c}},y_{c},T)$, which yields

[TABLE]

Using the previously obtained asymptotic expressions for $y_{\mathrm{c}}$, and inserting them into equation (30) for ${\cal H}B_{\mathrm{c}}$ to obtain asymptotic expressions for ${\cal H}B_{\mathrm{c}}$, subsequently results in asymptotic expressions when $T\rightarrow 0$ and $T\rightarrow 1$ for $\Sigma_{\mathrm{s}}$ and $W_{\mathrm{s}}$. However, since we in the present case have explicit solutions for $\phi_{\mathrm{c}}$ and ${\cal H}B_{\mathrm{c}}$, given in equation (52), we can give explicit expressions for $\Sigma_{\mathrm{s}}$ and $W_{\mathrm{s}}$ as functions of $T$:

[TABLE]

where we recall that $T=\lambda_{m}x^{3}/(1+\lambda_{m}x^{3})$. It follows that

[TABLE]

where the second result follows from $I\propto T^{\frac{5}{6}}$ when $T\rightarrow 0$. The Weyl scalar has the same limits. One can also infer the asymptotic rate of growth/decay of ${\Sigma}_{\mathrm{s}}$ and $W_{\mathrm{s}}$ as $x\rightarrow 0$ from equation (101):

[TABLE]

The limits (102) suggest that the growing mode (the orbit past asymptotic to the fixed point $P_{0+}$, given by $C_{-}=0$), approximates an exact solution with an isotropic singularity. The asymptotic decay rate in (103) agrees with the asymptotic results of Lim et al. (2004) [13] (see equations (4.19) and (5.9)). On the other hand, the unboundedness of the Hubble-normalized shear and Weyl perturbation when $C_{-}\neq 0$ indicates that generic perturbations become physically unviable, i.e. they do not approximate exact solutions of the Einstein field equations, when $T$ is sufficiently close to zero (recall that investigations of the Einstein field equations indicate that the Hubble-normalized shear and Weyl tensors are expected to be bounded generically).

One can also use (101) to determine the asymptotic behaviour of ${\Sigma}_{\mathrm{s}}$ and $W_{\mathrm{s}}$ at late times ($T\rightarrow 1,x\rightarrow\infty$). Since $\lim_{x\rightarrow\infty}I$ is finite it follows that

[TABLE]

for all $C_{+},C_{-}$. Equation (101) also gives the rates of decay along generic orbits asymptotic to the fixed point $P_{1+}(C_{+}I_{1}+C_{-}\neq 0)$:

[TABLE]

The result (104) for perturbations of $\Lambda$CDM is compatible with the future asymptotic behaviour of generalized $\Lambda$CDM universes toward a future de Sitter state, as described by equation (94), where equation (105) agrees with the asymptotic results of Lim et al (2004) [13] (see equations (3.24) and (5.8)).

5.2 Tensor perturbations

To analyze the tensor perturbations we make the transition (78) to Fourier space and use one of the four real-valued functions $h$ to represent the perturbation. Thus according to (96) the shear will be represented by ${\textstyle{1\over 2}}h^{\prime}$ and the electric Weyl tensor will be represented by ${\textstyle{1\over 2}}(h^{\prime}+k^{2}{\cal H}^{-2}h)$, which motivates defining the shear scalar and electric Weyl scalar for tensor perturbations according to

[TABLE]

where

[TABLE]

In order to analyze the behaviour of these scalars along the orbits we first rewrite the defining equation $y_{\mathrm{t}}=h^{\prime}/h$ in the form $d(\ln h)/dT=y_{\mathrm{t}}dN/dT=y_{\mathrm{t}}/(T(1-T))$ and then integrate to obtain

[TABLE]

It then follows from (106) that $\Sigma_{\mathrm{t}}$ and $W_{\mathrm{t}}$ are given as functions of $T$ along the orbits by

[TABLE]

We can now use the asymptotic expansions for the four fixed points in the previous section to determine the asymptotic form of $h(t)$ and hence of $\Sigma_{\mathrm{t}}$. We first consider the orbits that are past asymptotic to $P_{0+}$ and $P_{0-}$ as $T\rightarrow 0\,(x\rightarrow 0)$, referring to equations (86) and (87):

[TABLE]

One can further obtain the leading temporal $x$-dependence, as follows:

[TABLE]

Since ${\cal H}^{-2}(0)=0$ it follows from (109) that $W_{\mathrm{t}}$ has the same asymptotic behaviour as $\Sigma_{\mathrm{t}}$.

The above results describe the asymptotic behaviour of the perturbed shear and electric Weyl scalars near the initial singularity. The results are the same as for scalar perturbations. For example, only the single orbit $\mathrm{P}_{0+}\rightarrow\mathrm{P}_{1+}$ in figure 6 and figure 7 is compatible with an isotropic singularity and this orbit is thus analogous to the growing mode orbit $\mathrm{P}_{0+}\rightarrow\mathrm{P}_{1+}$ for scalar perturbations in figure 1. Also, analogously with the scalar perturbations, generically tensor perturbations result in unbounded Hubble-normalized shear and Weyl tensors, suggesting that the perturbations asymptotically are no longer approximations to the Einstein field equations when approaching $\mathrm{P}_{0-}$ toward the past.

We now consider the orbits that are future asymptotic to the fixed points $P_{1+}$ and $P_{1-}$ as $T\rightarrow 1\,(x\rightarrow\infty)$, referring to equations (90) and (89):

[TABLE]

One can further obtain the leading $x$-dependence, as follows:

[TABLE]

Again, $W_{\mathrm{t}}$ has the same asymptotic behaviour as $\Sigma_{\mathrm{t}}$. Here there is one difference between the scalar and tensor perturbations: both $\Sigma_{\mathrm{s}}$ and $\Sigma_{\mathrm{t}}$ are ${\cal O}(x^{-2})$ as $x\rightarrow\infty$, while $W_{\mathrm{s}}$ is ${\cal O}(x^{-2})$ and $W_{\mathrm{t}}$ is ${\cal O}(x^{-3})$. Since the asymptotic perturbations agree with the asymptotic results for the full Einstein equations in (94) one expects that the perturbations will approximate exact solutions at late times toward the future asymptotic de Sitter state.262626As regards fluids with rotation, vector perturbations $B_{i}$ are given by $B_{i}=b_{i}x^{-2}$, where $b_{i}$ depends on the spatial coordinates only, as follows from equation (58a) in [19]. Equations (58b), (42a), (66d) in the same reference yields $\tilde{v}_{i}$, which together with $B_{i}$ when inserted into equation (B.41c) for the vector mode of ${}^{(1)}\!\sigma^{i}\!_{j}/{}^{(0)}\!H$ results in that this quantity is $\propto x^{-\frac{3}{2}}$ when $x\rightarrow 0$. Hence the vector mode has to be set to zero in the case of an isotropic singularity. Equation (B.41d) in [19] for the fluid rotation also results in an initial blow up, unless the vector mode is set to zero.

6 Concluding remarks

The purpose of this paper has been to develop a new approach to using dynamical systems methods to analyze linear perturbations on a spatially flat RW background. We decided to use the $\Lambda$CDM model to illustrate the method because of its importance in cosmology and because of its relative mathematical simplicity. In our approach the state space ${\cal S}$ of the dynamical system has a product structure

[TABLE]

where ${\cal B}$ is the background state space which describes the dynamics of the flat RW background, and ${\cal P}$ is the perturbation state space, which contains the gauge invariant perturbation variables.

The Einstein equations in the RW background give a system of autonomous differential equations on ${\cal B}$ and the linearly perturbed Einstein equations give a system of autonomous differential equations for the perturbation variables in ${\cal P}$, which involve the background variables in ${\cal B}$. In this way the dynamics in the background determine the dynamics of the perturbations. The advantage of the product structure ${\cal S}={\cal B}\times{\cal P}$ is that when an orbit on ${\cal B}\times{\cal P}$ is projected onto the background ${\cal B}$ it coincides with an orbit on ${\cal B}$.

A key step is to choose bounded variables so that the state space is compactified and the system of autonomous differential equations is regular. It is also desirable to take advantage of the fact that the Einstein equations make it possible to decouple some of the variables, leaving a reduced state space to describe the essential dynamics.

The mathematical simplicity of the $\Lambda$CDM model is reflected in the fact that it is possible to use only one background variable and one perturbation variable to describe the essential dynamics. In other words both ${\cal B}$ and ${\cal P}$ are one dimensional spaces. For example, in the case of scalar perturbations in the uniform curvature gauge we represented ${\cal B}$ as the unit line segment $L$ with $T=\Omega_{\Lambda}$ as the background variable, and ${\cal P}$ as the circle ${S}^{1}$ described by the angular variable $\theta_{\mathrm{c}}$ . The state space is thus ${\cal S}={L}\times{S}^{1},$ which is a finite segment of a cylinder. Because the state space is bounded we were able to give a global description of the dynamics, in particular the behaviour at early and late times and the evolution at intermediate stages that may be of physical interest. In addition the differential equations, using $e-$fold time $N$ are well-suited for performing numerical simulations.

In future papers we will show how to obtain reduced and compactified product state spaces ${\cal B}\times{\cal P}$ with systems of regular differential equations for scalar field models and models with multiple sources. As a first step we will consider the simplest case, namely a minimally coupled scalar field with exponential potential, for which ${\cal B}$ is two dimensional and ${\cal P}$ is one dimensional. This case will provide the basis for the generalization to the case of a scalar field with more general potentials, which requires that ${\cal B}$ is three dimensional. This will also illustrate how one can organize perturbation theory into hierarchical structures where simpler models act as building blocks for more complicated ones.

Acknowledgments

AA is funded by the FCT grant SFRH/BPD/85194/2012, and supported by the project (GPSEinstein) PTDC/MAT-ANA/1275/2014, and CAMGSD, Instituto Superior Técnico by FCT/Portugal through UID/MAT/04459/2013. Furthermore, AA thanks the warm hospitality of Karlstad University. CU would like to thank the CAMGSD, Instituto Superior Técnico in Lisbon and the University of Waterloo, Canada, for kind hospitality.