Non-linear Cosmological Perturbations for Coupled Dark Energy

Testing gravity at cosmological scales [1, 2, 3] is finally becoming possible thanks to several large-scale surveys like BOSS [4], DES [5], DESI [6], Euclid [7] and, in the near future [8], LSST [9] and SKA [10, 11].

A common characteristic of most alternatives to Einstein’s gravity is the introduction of a non-minimal coupling of a scalar field to matter or to the metric, thereby adding at least a new parameter to the cosmological set. In this paper we focus on what could be classified as the simplest model of modified gravity: a constant coupling between a scalar field and dark matter. The scalar field can drive the cosmic acceleration, and therefore be a form of dark energy, but this is not a necessary condition. This dark-dark coupling bypasses the stringent conditions of local gravity since ordinary matter (i.e. baryons) are left uncoupled and feel only standard gravity. We refer to this model simply as coupled dark energy (CDE) [12, 13, 14, 15, 16, 17].

Previous work on the observable effects of this minimal model has focused on the linear regime (e.g. [18, 19, 20, 21, 22, 23]), corresponding to very large scales, typically larger than 50 Mpc$/h$. It is however now clear that the mildly non-linear regime, reaching down to roughly 10 Mpc$/h$, contains a wealth of additional information, allowing, in particular, to break some degeneracies among cosmological parameters.

The extension of cosmological data analysis to the non-linear regime can be performed in two ways: either resorting to $N$-body simulations, or by going higher in perturbation theory. The first avenue is of course more powerful and is clearly preferable if one deals with a restricted class of models that can be reliably and efficiently simulated. In the present context of uncertainty about what could be the preferred class of models beyond $\Lambda$CDM, however, the perturbation approach has still some appeal, since it can straightforwardly encompass large classes of models.

In this paper we consider therefore the next-to-leading order perturbation for a minimal dark-dark model, characterized by a coupling constant $\beta$ and a generic scalar field potential. We consider two potentials: an exponential potential, as suggested by theoretical considerations [24]; and a linear approximation that is suitable when the scalar field does not change much during the relevant period (i.e., within the range of observations). The linear potential has the advantage that an analytical solution for the background can be obtained. In the limit in which the potential reduces to a cosmological constant, this coupled model introduces a single constant beyond standard, and represents therefore a minimal modified gravity model.

In higher-order perturbation theory, the perturbation variables are expressed as convolutions of the first order variables. The main difficulty consists therefore in deriving the kernels of the convolutions. In this paper we obtain the kernels for a generic cosmology following the approach of [25] and then we specialize our results to the CDE model, providing analytical and/or numerical forms of the kernels as a function of time, of the coupling parameter, and of the potential slope. These forms are ready to be used in future work to forecast the performance of cosmological surveys and to analyse real data.

In our model, the dark energy scalar field (subscript $\phi$) is coupled to the non-relativistic matter (subscript $m$) field in a flat FLRW background. Therefore the individual energy-momentum tensors are not conserved due to interaction terms on the right-hand side:

Here $\beta$ is a dimensionless constant coupling term that quantifies the strength of the dark-dark interaction, while $T_{m}$ refers to the trace of the energy-momentum tensor for pressureless matter. Baryonic matter is assumed to be uncoupled so that the local gravity constraints are automatically satisfied. In CDE, the gravitational equations remain the same as in standard gravity. Other formulations of such interaction terms can be found in the literature [26, 27, 28, 29, 30].

The coupling constant $\beta$ has been constrained, mainly by CMB, to be smaller than 0.04 roughly [18]. However, this upper limit assumes a constant value across the entire cosmic evolution. A variable coupling might in general take larger values. Here, therefore, we will explore also values up to $\beta=0.2$.

The time component ($\nu=0$) of Eq. (2.1) provides evolution equations for dark matter, dark energy and radiation (subscript $r$) respectively:

Here, the prime is the derivative with respect to $\eta=\log(a)$. The radiation Eq. (2.4) is independent of the coupling since its energy-momentum tensor is traceless.

The Friedmann equation reads

The conservation equation for $\phi$ can also be written in the Klein-Gordon form

We now rewrite the background equations in terms of dimensionless parameters. These parameters can be defined as follows:

The quantities $\Omega_{K}=x^{2}$, $\Omega_{P}=y^{2}$, $\Omega_{r}=z^{2}$, and $\Omega_{b}=v^{2}$ represent the fractions of total energy associated with field kinetic energy, field potential energy, radiation, and baryons, respectively. Additionally, the field energy fraction is given by $\Omega_{\phi}=x^{2}+y^{2}$. The effective state parameter, defined as the total pressure-to-density ratio, is given in terms of the parameters in (2.7) as $w_{\rm eff}=x^{2}-y^{2}+\frac{z^{2}}{3}$; while the field state parameter, $w_{\phi}=\frac{P_{\phi}}{\rho_{\phi}}$, is expressed as $w_{\phi}=\frac{x^{2}-y^{2}}{x^{2}+y^{2}}$. Since the total energy fractions must sum to one, the cold dark matter energy density fraction is given by $\Omega_{c}=1-x^{2}-y^{2}-z^{2}-v^{2}$.

Using Eq. (2.7) we obtain an autonomous dynamical system:

where $\mu=\frac{1}{V}\frac{\partial V}{\partial\phi}$ is the potential slope parameter. To gain insight into the dynamics of the system, one can identify the critical points of the phase space and evaluate their stability conditions. This has been already investigated in detail and will not be repeated here (see e.g. the review [31]). However, for the scope of our work, we will specifically address one critical point, namely the existence of a matter-dominated era in which the scalar field is not negligible. This critical point, known as ${\phi}$MDE, will be detailed in Section 4.2.

We now move to higher orders in perturbations. We adopt the perturbation scheme of Eulerian standard perturbation theory. The general form of the conservation equations in the Newtonian limit and in Fourier space are (see e.g. the review in [33])

where we defined

Here, $(\tilde{\delta}=\frac{\rho}{\rho_{0}}-1)$ is the density contrast, $\tilde{\theta}=\frac{ik_{i}v^{i}}{Ha}$ is the peculiar velocity divergence field, and the time parameter is $\eta=\log(a)$. We consider the matter distribution to act as a pressureless fluid without vorticity. Additionally, we assume that the peculiar velocities are sufficiently small to be treated as non-relativistic and that we are at deep sub-horizon scales, $k\gg aH$, where the Newtonian approximations hold.

To solve these equations, we need to apply perturbation theory, by expressing the solution as a functional of the linear matter density, $\delta^{(1)}$. This approach may break down at very small scales due to additional stochastic effects from short-wavelength modes [34, 35, 36]; however, in the weakly nonlinear regime considered here, this effect remain negligible. Consequently, we can expand the solution perturbatively [33], as shown below:

By inserting the expression (3.5) into the conservation equations (3.1),(3.2) and solving for each order separately, one can obtain the following general form for nth order:

The $\delta_{\mathbf{q}_{n}}(\eta)$ terms on the right-hand side represent linear density contrast , which evolve from the initial density contrast via the growth function $D(\eta)$ as $\delta_{\mathbf{q}_{n}}(\eta)=D(\eta)\delta_{\mathbf{q}_{n}}$. The functions $F_{n}$ and $G_{n}$, known as kernels, arise from nonlinear terms in the conservation equations, and they describe the coupling between different modes. For instance, in the Einstein-de Sitter Universe (EdS) model, the kernels can be derived recursively, starting from the linear order [37, 33]. This form of the solution, as given in equations (3.6) and (3.7), can also be applied to non-EdS models. However, determining the kernel structure for these models can be quite challenging. Therefore, we need to explore alternative methods.

In this work, we employ the perturbation theory kernel structure studied in the paper [25], which was obtained by using the symmetries of the conservation equations. This particular symmetry is extended Galilean invariance, where the equation of motion remains invariant under spatial translations with an arbitrary time parameter [38, 39]. This approach is particularly advantageous for scenarios beyond $\Lambda$CDM. For instance, following the procedure outlined in [25], one can express a generic second-order kernel as:

The $F_{2}$ matter kernel (3.8) is constructed from two basis functions, $\gamma$ and $\tilde{\beta}$, which depend on the external momenta $\mathbf{q}_{1}$ and $\mathbf{q}_{2}$. These functions are homogeneous, meaning $F(\lambda\mathbf{q}_{1},\lambda\mathbf{q}_{2})=F(\mathbf{q}_{1},\mathbf{q}_{2})$, and symmetric under the exchange of $\mathbf{q}_{1}$ and $\mathbf{q}_{2}$, so that $F(\mathbf{q}_{1},\mathbf{q}_{2})=F(\mathbf{q}_{2},\mathbf{q}_{1})$. Due to the isotropy of the Universe, the basis functions are also rotationally invariant. For external momenta $\mathbf{p}$ and $\mathbf{q}$, using rotational invariance and homogeneity, one can construct basis functions as follows [25]

Note that $\gamma$ and $\tilde{\beta}$ are symmetric, while $\alpha_{\alpha}$ is an anti-symmetric function. Here, $a_{0}^{(2)}$, $a_{1}^{(2)}$, and $a_{2}^{(2)}$ are time-dependent coefficients, where $a_{0}^{(2)}$ and $a_{2}^{(2)}$ can be determined by applying constraints derived from the symmetries of the conservation equations. The remaining undetermined coefficient, $a_{1}^{(2)}$, can be fixed by our cosmological model. The most general forms of the matter kernels (i.e density contrast kernel $F_{i}$ and the velocity divergence kernel $G_{i}$), after applying the constraints, are presented below [25]:

and

The $h$ function is defined as:

and $\mathbf{q}_{ij}=\mathbf{q}_{i}+\mathbf{q}_{j}$ is the sum of momenta. At second order, $F_{2}$ depends on the single undetermined coefficient $a^{(2)}_{1}$, while $G_{2}$ depends on $d^{(2)}_{1}$. Moving to third order, $F_{3}$ includes three undetermined coefficients: $a^{(2)}_{1}$, $a^{(3)}_{5}$, and $a^{(3)}_{10}$, whereas $G_{3}$ includes $d^{(2)}_{1}$, $d^{(3)}_{5}$, and $d^{(3)}_{10}$.

We are finally in place to express the kernels for our coupled dark energy model. After recalling the standard result for EdS, we examine explicitly three cases: 1) the $\phi$MDE epoch; 2) a linear potential; 3) an exponential potential. In the first case we can provide analytical solutions, while in the others we need to perform numerical integrations.