Second-order Cosmological Perturbations Engendered by Point-like Masses

TL;DR

This paper develops a second-order cosmological perturbation theory within the $$CDM model, accounting for discrete point-like masses and nonlinear density contrasts, suitable for relativistic simulations and backreaction studies.

Contribution

It introduces a new second-order perturbation framework for nonrelativistic matter inhomogeneities modeled as point masses, valid at arbitrary distances and nonlinear contrasts.

Findings

Derived equations for second-order metric corrections including scalar, vector, and tensor modes.

Verified Einstein equations and consistency in the weak field limit.

Outlined potential for relativistic simulations of cosmological backreaction.

Abstract

In the $\Lambda$CDM framework, presenting nonrelativistic matter inhomogeneities as discrete massive particles, we develop the second-order cosmological perturbation theory. Our approach relies on the weak gravitational field limit. The derived equations for the second-order scalar, vector and tensor metric corrections are suitable at arbitrary distances including regions with nonlinear contrasts of the matter density. We thoroughly verify fulfilment of all Einstein equations as well as self-consistency of order assignments. In addition, we achieve logical positive results in the Minkowski background limit. Feasible investigations of the cosmological backreaction manifestations by means of relativistic simulations are also outlined.

Equations

$$d{s}^2=a^2\left(d{\eta }^2-{\delta }_{\alpha \beta }d{x}^{\alpha }d{x}^{\beta }\right),\alpha ,\beta =1,2,3,$$

$$\frac{3{\mathcal{H}}^2}{a^2}=\kappa \overline{\epsilon }+\Lambda ,\frac{2{\mathcal{H}}^{\prime }+{\mathcal{H}}^2}{a^2}=\Lambda .$$

$$d{s}^2=a^2\left[\left(1+2\Phi \right)d{\eta }^2+2B_{\alpha }d{x}^{\alpha }d\eta -\left(1-2\Phi \right){\delta }_{\alpha \beta }d{x}^{\alpha }d{x}^{\beta }\right].$$

$$\nabla \mathbf{B}\equiv {\delta }^{\alpha \beta }\frac{\partial B_{\alpha }}{\partial x^{\beta }}\equiv 0.$$

$$\Phi$$

$$\mathbf{B}$$

$$\rho (\eta ,\mathbf{r})=\sum _n{m}_n\delta (\mathbf{r}-{\mathbf{r}}_n)=\sum _n{\rho }_n,{\rho }_n(\eta ,\mathbf{r})\equiv m_n\delta (\mathbf{r}-{\mathbf{r}}_n).$$

$$△\Phi -\frac{3\kappa \overline{\rho }{c}^2}{2a}\Phi =\frac{\kappa c^2}{2a}\delta \rho -\frac{3\kappa c^2}{2a}\mathcal{H}\Xi ,$$

$$△\mathbf{B}-\frac{2\kappa \overline{\rho }{c}^2}{a}\mathbf{B}=-\frac{2\kappa c^2}{a}\left(\sum _n{\rho }_n{\tilde{\mathbf{v}}}_n-\nabla \Xi \right),$$

$$\Xi =\frac{1}{4\pi }\sum _n{m}_n\frac{(\mathbf{r}-{\mathbf{r}}_n){\tilde{\mathbf{v}}}_n}{|\mathbf{r}-{\mathbf{r}}_n{|}^3}.$$

$${\Phi }^{\prime }+\mathcal{H}\Phi =-\frac{\kappa c^2}{2a}\Xi ,{\Phi }^{\prime \prime }+3\mathcal{H}{\Phi }^{\prime }+\left(2{\mathcal{H}}^{\prime }+{\mathcal{H}}^2\right)\Phi =0,{\mathbf{B}}^{\prime }+2\mathcal{H}\mathbf{B}=0.$$

$$△\Xi =\nabla \sum _n{\rho }_n{\tilde{\mathbf{v}}}_n=-\sum _n{\rho }_n^{\prime },$$

$$d{s}^2$$

$$\nabla {\mathbf{B}}^{(2)}\equiv {\delta }^{\alpha \beta }\frac{\partial B_{\alpha }^{(2)}}{\partial x^{\beta }}\equiv 0.$$

$${\delta }^{\alpha \beta }\frac{\partial h_{\alpha \gamma }}{\partial x^{\beta }}\equiv 0,{\delta }^{\alpha \beta }{h}_{\alpha \beta }\equiv 0.$$

$$G_i^k=\kappa T_i^k+\Lambda {\delta }_i^k,i,k=0,1,2,3,$$

$$G_i^k={\left(G_i^k\right)}^{(0)}+{\left(G_i^k\right)}^{(1)}+{\left(G_i^k\right)}^{(2)}+{\left(G_i^k\right)}^{(\mathrm{II})}.$$

$$G_{\beta }^0={\left(G_{\beta }^0\right)}^{(0)}+{\left(G_{\beta }^0\right)}^{(1)}+{\left(G_{\beta }^0\right)}^{(2)}+{\left(G_{\beta }^0\right)}^{(\mathrm{II})},\beta =1,2,3,$$

$${\left(G_{\beta }^0\right)}^{(0)}=0,$$

$${\left(G_{\beta }^0\right)}^{(1)}=\frac{1}{2a^2}△B_{\beta }+\frac{2}{a^2}\mathcal{H}\frac{\partial \Phi }{\partial x^{\beta }}+\frac{2}{a^2}\frac{\partial {\Phi }^{\prime }}{\partial x^{\beta }},$$

$${\left(G_{\beta }^0\right)}^{(2)}=\frac{1}{2a^2}△B_{\beta }^{(2)}+\frac{2}{a^2}\mathcal{H}\frac{\partial {\Phi }^{(2)}}{\partial x^{\beta }}+\frac{2}{a^2}\frac{\partial {\left({\Psi }^{(2)}\right)}^{\prime }}{\partial x^{\beta }},$$

$${\left(G_{\beta }^0\right)}^{(\mathrm{II})}$$

$${\left(G_0^0\right)}^{(0)}=\frac{3{\mathcal{H}}^2}{a^2},{\left(G_{\beta }^{\alpha }\right)}^{(0)}=\left(\frac{2{\mathcal{H}}^{\prime }+{\mathcal{H}}^2}{a^2}\right){\delta }_{\alpha \beta },$$

$${\left(G_0^0\right)}^{(1)}=\frac{2}{a^2}△\Phi -\frac{6}{a^2}{\mathcal{H}}^2\Phi -\frac{6}{a^2}\mathcal{H}{\Phi }^{\prime },$$

$${\left(G_{\beta }^{\alpha }\right)}^{(1)}$$

$${\left(G_0^0\right)}^{(2)}=\frac{2}{a^2}△{\Psi }^{(2)}-\frac{6}{a^2}{\mathcal{H}}^2{\Phi }^{(2)}-\frac{6}{a^2}\mathcal{H}{\left({\Psi }^{(2)}\right)}^{\prime },$$

$${\left(G_{\beta }^{\alpha }\right)}^{(2)}=\frac{1}{2a^2}{h}_{\alpha \beta }^{\prime \prime }+\frac{1}{a^2}\mathcal{H}{h}_{\alpha \beta }^{\prime }-\frac{1}{2a^2}△h_{\alpha \beta }$$

$$T_0^0=\frac{\rho c^2}{\sqrt{-g}}\cdot \frac{g_{00}+{\tilde{v}}^{\gamma }{g}_{0\gamma }}{\sqrt{g_{00}+2g_{0\alpha }{\tilde{v}}^{\alpha }+g_{\alpha \beta }{\tilde{v}}^{\alpha }{\tilde{v}}^{\beta }}},$$

$$T_{\alpha }^0=\frac{\rho c^2}{\sqrt{-g}}\cdot \frac{g_{0\alpha }+{\tilde{v}}^{\beta }{g}_{\alpha \beta }}{\sqrt{g_{00}+2g_{0\mu }{\tilde{v}}^{\mu }+g_{\mu \nu }{\tilde{v}}^{\mu }{\tilde{v}}^{\nu }}},$$

$$T_{\beta }^{\alpha }=\frac{\rho c^2}{\sqrt{-g}}\cdot \frac{{\tilde{v}}^{\alpha }\left(g_{0\beta }+{\tilde{v}}^{\gamma }{g}_{\gamma \beta }\right)}{\sqrt{g_{00}+2g_{0\mu }{\tilde{v}}^{\mu }+g_{\mu \nu }{\tilde{v}}^{\mu }{\tilde{v}}^{\nu }}},$$

Full text

cosmological parameters — cosmology: theory — dark energy — dark matter — gravitation — large-scale structure of universe

1 INTRODUCTION

The conventional Lambda Cold Dark Matter ($\Lambda$CDM) model conforms with the observational data (see, in particular, Ade et al. (2016)) and embodies the mainstream of modern cosmology despite the distressing fact that the nature of dark ingredients of the Universe still remains uncertain. The key assumption, being typical for this cosmological model as well as its numerous alternatives, is the existence of the homogeneous and isotropic Friedmann-Lemaître-Robertson-Walker (FLRW) background, which is only slightly perturbed by inhomogeneities inherent in the distribution of the world’s filling material. The following quite natural question arises: can galaxies and their accumulations lead to considerable metric corrections and deeply affect the average cosmic expansion? The affirmative answer would bring to general recognition of the so-called backreaction (see the reviews by Räsänen (2011); Clarkson et al. (2011); Buchert & Räsänen (2012); Bolejko & Korzyński (2017) and Refs. therein; also for a recent spirited debate about the magnitude of backreaction effects see, in particular, Green & Wald (2013); Buchert et al. (2015) and Refs. therein) and give at least a glimmer of hope to explain the apparent acceleration of expansion without mysterious dark energy.

Even if the answer is negative, it is extremely important to estimate deviations from the FLRW description and confront the theoretical predictions with the promising outcomes of such future space missions as Euclid (Scaramella et al., 2015; Amendola et al., 2016). For this purpose a reliable cosmological perturbation theory should be developed in the general relativity (GR) framework. Appropriateness at all cosmic scales and non-perturbative treatment of the matter density are those basic requirements which such a theory must meet. The first-order scheme complying with these reasonable demands and being suitable for relativistic $N$-body simulations has been successfully constructed for a system of discrete massive particles with nonrelativistic velocities by Eingorn (2016). Some predecessors and their drawbacks are reviewed in (Eingorn, 2016) as well.

The current paper is devoted to the generalization of the afore-mentioned first-order approach to the second order with respect to deviations of the metric coefficients from the corresponding background quantities. This generalization is particularly motivated by the fact that the first order is obviously insufficient for the comprehensive analysis of the possible backreaction manifestations and the trustworthy prediction of their magnitude. The appeal to the second order with the purpose of revealing the corresponding observable features is also very promising in the era of precision cosmology (see Ben-Dayan et al. (2013); Bonvin et al. (2015a, b) and Refs. therein).

The narration is organized as follows: after reviewing the basic results of Eingorn (2016) in Section 2, we switch over to derivation and verification of equations for the second-order scalar, vector and tensor perturbations in Section 3. Then, in Section 4, we focus attention on identification of the effective average energy density and pressure, and propose the research program aimed at the perturbative computation of the cosmological backreaction effects. Our achievements are laconically summarized in concluding Section 5.

2 DISCRETE PICTURE OF THE FIRST-ORDER COSMOLOGICAL PERTURBATIONS AT ALL SCALES

Let us start with an overview of the newly formulated first-order perturbation theory covering all cosmological spatial scales and permitting of nonlinear contrasts of the matter density (Eingorn, 2016). We confine ourselves to the conventional $\Lambda$CDM model (with zero spatial curvature) and concentrate on those stages of the Universe evolution when cold matter (dark and baryonic) and dark energy (being represented by the cosmological constant $\Lambda$) dominate while radiation or relativistic cosmic neutrinos make negligible contributions (see, however, Eingorn & Brilenkov (2015); Eingorn et al. (2016) for a broad generalization to the multicomponent case). Then the homogeneous and isotropic cosmological background is described by the unperturbed FLRW metric

[TABLE]

where $\eta$ and $x^{\alpha}$, $\alpha=1,2,3$, stand for the conformal time and comoving coordinates, respectively, and the corresponding Friedmann equations for the scale factor $a(\eta)$:

[TABLE]

Here ${\mathcal{H}}(\eta)\equiv a^{\prime}/a$ (with prime denoting the derivative with respect to $\eta$) and $\kappa\equiv 8\pi G_{N}/c^{4}$ (with $c$ and $G_{N}$ representing the speed of light and Newtonian gravitational constant, respectively). Further, $\varepsilon$ is the energy density of the nonrelativistic pressureless matter, and the overline indicates averaging.

In the first-order approximation the real inhomogeneous Universe is usually assumed to be described well by the perturbed metric (Bardeen, 1980; Mukhanov, 2005; Gorbunov & Rubakov, 2011; Durrer, 2008)

[TABLE]

Here $\Phi(\eta,{\bf r})$ is the scalar perturbation while the spatial vector ${\bf B}(\eta,{\bf r})\equiv\left(B_{1},B_{2},B_{3}\right)\equiv\left(B_{x},B_{y},B_{z}\right)$ stands for the vector perturbation and satisfies the prevalent gauge condition

[TABLE]

Similarly to Clarkson & Umeh (2011); Baumann et al. (2012), we have chosen the popular Poisson gauge but have not yet taken account of the tensor perturbations treating them as second-order quantities. As pointed out by Eingorn (2016), one can in principle allow for the first-order tensor perturbations associated with freely propagating gravitational waves (with no generator). They have not been explicitly included in (2.3) since below we totally neglect their possible contributions to the sources of second-order metric corrections. At the same time, similarly to Eingorn (2016) but in contrast to Clarkson & Umeh (2011); Baumann et al. (2012), the vector perturbation ${\bf B}$ has been already included in (2.3) as a first-order quantity since it has a definite nonzero generator (see the right-hand side (rhs) of the equation (2.9) for ${\bf B}$ below).

In (Eingorn, 2016) the role of the inhomogeneous gravitational field source belongs to a system of separate nonrelativistic point-like particles with masses $m_{n}$, comoving radius-vectors ${\bf r}_{n}(\eta)$ and comoving peculiar velocities ${\bf\tilde{v}}_{n}(\eta)\equiv d{\bf r}_{n}/d\eta\equiv(\tilde{v}_{n}^{1},\tilde{v}_{n}^{2},\tilde{v}_{n}^{3})$, and the following expressions for $\Phi$, ${\bf B}$ are derived:

[TABLE]

[TABLE]

Here $q_{n}(\eta,{\bf r})\equiv a|{\bf r}-{\bf r}_{n}|/\lambda$, with $\lambda(\eta)\equiv\left[2a^{3}/\left(3\kappa\overline{\rho}c^{2}\right)\right]^{1/2}\sim a^{3/2}$ defining a finite range of Yukawa interaction. Being armed with the observed values of the Hubble constant $H_{0}\approx 68\,\mathrm{km\,s^{-1}\,Mpc^{-1}}$ and the parameter $\Omega_{\mathrm{M}}\equiv\kappa\overline{\rho}c^{4}/\left(3H_{0}^{2}a_{0}^{3}\right)\approx 0.31$ (Ade et al., 2016), where $a_{0}$ denotes the today’s scale factor, one can easily estimate the current value of the introduced characteristic cutoff scale: $\lambda_{0}\approx 3.7\,\mathrm{Gpc}$ (Eingorn, 2016). An evident relationship $\overline{\varepsilon}=\overline{\rho}c^{2}/a^{3}$ establishes linkage between $\overline{\varepsilon}(\eta)\sim a^{-3}$ and the constant average rest mass density $\overline{\rho}$. As regards the corresponding non-averaged quantity, the rest mass density of the analyzed particle system in the comoving coordinates has the form

[TABLE]

The explicit analytical expressions (2.5) and (2.6) and their noteworthy features along with the underlying perturbative approach and its numerous physical implications and advantages are analyzed in detail by Eingorn (2016), and it makes no sense to manifest hair-splitting by repeating all deserving results here. Nevertheless, let us briefly enumerate those facts which are crucial for the clear purpose of the next section.

First of all, we stubbornly adhere to the following well-grounded argumentation (see also Baumann et al. (2012) for the similar reasoning): the metric corrections and peculiar particle velocities are assumed to be small at arbitrary distances, however, the smallness of the density contrast $\delta\equiv\delta\rho/\overline{\rho}$ (where $\delta\rho\equiv\rho-\overline{\rho}$) is not demanded. For instance, $\rho\gg\overline{\rho}$ within galaxies, but even for huge density contrasts the metric is still approximated well by (2.3). Therefore, one can linearize the Einstein equations in the first-order metric perturbations $\Phi$; ${\bf B}$ and velocities ${\bf\tilde{v}}_{n}$ without resorting to the unnecessary restrictive inequality $|\delta|\ll 1$. Thus, the nonlinear deviation of $\rho$ from its average value $\overline{\rho}$, actually occurring at sufficiently small scales, is absolutely unforbidden and fully taken into consideration, as opposed to the standard first-order relativistic perturbation theory (Mukhanov, 2005; Durrer, 2008; Gorbunov & Rubakov, 2011) and its generalization to the second order (see, for instance, Bartolo et al. (2006)).

Secondly, the functions (2.5) and (2.6) are found as exact solutions of the Helmholtz equations (see Eingorn (2016) for the corresponding Fourier transforms):

[TABLE]

[TABLE]

where $\triangle\equiv\delta^{\alpha\beta}\cfrac{\partial^{2}}{\partial x^{\alpha}\partial x^{\beta}}$ is the Laplace operator in the comoving coordinates, and the auxiliary function $\Xi(\eta,{\bf r})$ has been introduced:

[TABLE]

Evidently, the presence of the term $\sim\Phi$ in the equation (2.8) for the scalar perturbation $\Phi$ is the reason for the Yukawa-type cutoff in the solution (2.5). The afore-mentioned screening length $\lambda$ is naturally related to the factor in this term: $3\kappa\overline{\rho}c^{2}/(2a)\equiv a^{2}/\lambda^{2}$. The contribution $\sim\Phi$ itself as well as the contribution $\sim{\bf B}$ in the equation (2.9) arise in view of the fact that the energy-momentum fluctuations, generating the metric corrections, contain these corrections themselves (see the formulas (2.13) and (2.14) in (Eingorn, 2016) as well as the expressions (3.22)–(3.24) below).

The vector perturbation ${\bf B}$ (2.6) obeys the gauge condition (2.4). In addition, within the adopted accuracy the perturbations (2.5) and (2.6) satisfy all remaining linearized Einstein equations, which are reduced to the triplet containing temporal derivatives:

[TABLE]

Finally, the auxiliary function (2.10) is found as the exact solution of the Poisson equation:

[TABLE]

where the continuity equation $\rho^{\prime}_{n}+\nabla(\rho_{n}{\bf\tilde{v}}_{n})=0$, which is satisfied for any $n$-th particle identically, has been employed.

It presents no difficulty to show that the usual theory of hydrodynamical fluctuations emerges in the continuum limit. Indeed, let us momentarily regard $\rho(\eta,{\bf r})$ as a continuous mass density field and ${\bf\tilde{v}}(\eta,{\bf r})=\nabla\tilde{v}^{(\parallel)}+{\bf\tilde{v}}^{(\perp)}$ as a continuous velocity field (with $\tilde{v}^{(\parallel)}(\eta,{\bf r})$ and ${\bf\tilde{v}}^{(\perp)}(\eta,{\bf r})$ denoting the scalar and vector parts, respectively). Further, we replace $\Xi$ by $\overline{\rho}\tilde{v}^{(\parallel)}$ and $\left(\sum\limits_{n}\rho_{n}{\bf\tilde{v}}_{n}-\nabla\Xi\right)$ by $\overline{\rho}{\bf\tilde{v}}^{(\perp)}$ in the case of linear density fluctuations. Then we make use of the relationship $\varepsilon=\rho c^{2}/a^{3}+\left(3\overline{\rho}c^{2}/a^{3}\right)\Phi$ in order to rewrite equations in terms of the energy density field $\varepsilon(\eta,{\bf r})$ (instead of $\rho$). Finally, we introduce covariant spatial components of the 4-velocity $u_{\beta}=aB_{\beta}-a\tilde{v}^{\beta}$, $\beta=1,2,3$, and a spatial vector with these components ${\bf u}(\eta,{\bf r})\equiv\left(u_{1},u_{2},u_{3}\right)$ (instead of ${\bf\tilde{v}}$). We also single out the corresponding scalar and vector parts: ${\bf u}=\nabla u^{(\parallel)}+{\bf u}^{(\perp)}$, where $u^{(\parallel)}(\eta,{\bf r})=-a\tilde{v}^{(\parallel)}$ and ${\bf u}^{(\perp)}(\eta,{\bf r})=a{\bf B}-a{\bf\tilde{v}}^{(\perp)}$. As a result, the rewritten equations (2.8), (2.9), (2.11) become equivalent to those describing hydrodynamical perturbations (see Mukhanov (2005), section 7.3, in the pressureless matter case).

Since the first-order metric corrections inevitably contribute to the second-order ones due to nonlinearity of the Einstein equations, we make the most out of the enumerated equations for $\Phi$ and ${\bf B}$ as well as $\Xi$ in the very next section.

3 SPLENDORS OF THE SECOND-ORDER THEORY

3.1 Sources of Perturbations

Switching over to the second-order approximation, we present the metric as

[TABLE]

where $\Phi^{(2)}(\eta,{\bf r})$ and $\Psi^{(2)}(\eta,{\bf r})$ are the second-order scalar perturbations while the spatial vector ${\bf B}^{(2)}(\eta,{\bf r})\equiv\left(B_{1}^{(2)},B_{2}^{(2)},B_{3}^{(2)}\right)\equiv\left(B_{x}^{(2)},B_{y}^{(2)},B_{z}^{(2)}\right)$ stands for the second-order vector perturbation and satisfies the same gauge condition as (2.4) for ${\bf B}$:

[TABLE]

Similarly to Clarkson & Umeh (2011); Baumann et al. (2012), we have included the second-order tensor perturbations $h_{\alpha\beta}$, $\alpha,\beta=1,2,3$, in (3.1). They obey the standard “transverse-traceless” gauge conditions:

[TABLE]

It should be noted that we decompose the second-order metric corrections into scalar, vector and tensor modes with respect to the unperturbed background. Thus, for example, the covariant divergence of the total vector ${\bf B}+{\bf B}^{(2)}$, defined via the perturbed spatial metric, does not vanish under the made gauge choice.

In order to elaborate on the Einstein equations

[TABLE]

for the sought-for functions $\Phi^{(2)}$, $\Psi^{(2)}$, ${\bf B}^{(2)}$ and $h_{\alpha\beta}$, first of all it is necessary to calculate the mixed components of the Einstein tensor $G_{i}^{k}$ and matter energy-momentum tensor $T_{i}^{k}$ up to the second order in metric corrections and their sources. Let us start with presenting the left-hand side (lhs) of (3.4) as follows (see also Baumann et al. (2012) for a similar decomposition):

[TABLE]

Here $\left(G_{i}^{k}\right)^{(0)}$ corresponds to the unperturbed cosmological background metric (2.1). The expressions $\left(G_{i}^{k}\right)^{(1)}$ (or, by analogy, $\left(G_{i}^{k}\right)^{(2)}$) are constructed from the terms being linear in the first-order (or second-order) perturbations $\Phi$ and ${\bf B}$ (or $\Phi^{(2)}$, $\Psi^{(2)}$, ${\bf B}^{(2)}$ and $h_{\alpha\beta}$) and their spatio-temporal derivatives. Finally, $\left(G_{i}^{k}\right)^{(\mathrm{II})}$ represent the second-order quantities containing products of $\Phi$, ${\bf B}$ and their derivatives. We express the temporal derivatives $\Phi^{\prime}$, $\Phi^{\prime\prime}$ and ${\bf B}^{\prime}$ in $\left(G_{i}^{k}\right)^{(\mathrm{II})}$ by means of the functions $\Phi$ and ${\bf B}$ themselves with the help of the triplet (2.11). Thus, for example,

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here ${\bf B}^{2}\equiv B_{x}^{2}+B_{y}^{2}+B_{z}^{2}$. The rest of the expressions $\left(G_{i}^{k}\right)^{(0)}$ and $\left(G_{i}^{k}\right)^{(1)}$ are well-known (Mukhanov, 2005; Durrer, 2008; Gorbunov & Rubakov, 2011):

[TABLE]

[TABLE]

[TABLE]

The expressions $\left(G_{i}^{k}\right)^{(2)}$ bear a strong resemblance to $\left(G_{i}^{k}\right)^{(1)}$ with due regard for the inequality $\Phi^{(2)}\neq\Psi^{(2)}$ as well as additional contributions from the tensor perturbations:

[TABLE]

[TABLE]

Finally, the quantities $\left(G_{0}^{0}\right)^{(\mathrm{II})}$ and $\left(G_{\beta}^{\alpha}\right)^{(\mathrm{II})}$ are quite cumbersome, and it makes no sense to reproduce them here.

As regards the mixed energy-momentum tensor components in the rhs of (3.4), we resort to the well-known formulas for the analyzed system of point-like particles (Landau & Lifshitz, 2000; Eingorn & Zhuk, 2014):

[TABLE]

[TABLE]

[TABLE]

where $g\equiv\mathrm{det}(g_{ik})$. In addition, $\tilde{v}^{\alpha}$ coincides with $\tilde{v}_{n}^{\alpha}$ in the term with the factor $\delta({\bf r}-{\bf r}_{n})$. For instance,

[TABLE]

With the help of (3.1) we immediately write down the metric coefficients $g_{ik}$ up to the second order:

[TABLE]

Consequently, with the same accuracy

[TABLE]

and the formulas (3.16), (3.17) and (3.18) are reduced to the formulas

[TABLE]

[TABLE]

[TABLE]

respectively. Here $\tilde{v}^{2}\equiv\delta_{\alpha\beta}\tilde{v}^{\alpha}\tilde{v}^{\beta}$. By analogy with Eingorn (2016), we consider the quantities $\delta\rho$ and ${\bf\tilde{v}}_{n}$ as “importing” the first order of smallness in the Einstein equations. In other words, metric corrections, which are generated by these sources themselves, are assigned the first order. Indeed, in the above-mentioned equations (2.8) and (2.9) the quantities $\delta\rho$ and ${\bf\tilde{v}}_{n}$ (as well as $\Xi$ containing ${\bf\tilde{v}}_{n}$) in the rhs play the role of sources generating the first-order metric corrections $\Phi$ and ${\bf B}$ forming the lhs. Therefore, we have omitted such terms as, for example, $\sim\delta\rho\Phi^{2}$ in (3.22) since this term is much smaller than the summand $\sim\delta\rho\Phi$ in the same parentheses at all cosmological scales and would import the third order of smallness (see also Chisari & Zaldarriaga (2011) for the similar reasoning). The established thorough separation of the first- and second-order summands in the Einstein equations is strongly corroborated in subsection 3.4 below.

Once again we employ a helpful decomposition

[TABLE]

The only nonzero component with the superscript “(0)” is $\left(T_{0}^{0}\right)^{(0)}=\overline{\rho}c^{2}/a^{3}$. The components $\left(T_{i}^{k}\right)^{(1)}$ (or, by analogy, $\left(T_{i}^{k}\right)^{(2)}$) are constructed from the terms being linear in the quantities $\delta\rho$, ${\bf\tilde{v}}$, $\Phi$, ${\bf B}$ (or $\Psi^{(2)}$, ${\bf B}^{(2)}$):

[TABLE]

[TABLE]

Finally, the components with the superscript “(II)” contain products of $\delta\rho$, ${\bf\tilde{v}}$, $\Phi$ or ${\bf B}$:

[TABLE]

It is worth mentioning that owing to the zero value of $\left(T^{\alpha}_{\beta}\right)^{(1)}$ the anisotropic stress vanishes in the first-order approximation. This is the cogent reason for using the same designation $\Phi\equiv\Phi^{(1)}=\Psi^{(1)}$ for the equal first-order scalar perturbations $\Phi^{(1)}$ and $\Psi^{(1)}$ from the very outset. Nevertheless, owing to the nonzero values of $\left(G^{\alpha}_{\beta}\right)^{(\mathrm{II})}$ and $\left(T^{\alpha}_{\beta}\right)^{(\mathrm{II})}$ (or, more precisely, their combinations $Q_{\alpha\beta}$ introduced below) the anisotropic stress does not vanish in the second-order approximation. Thus, generally speaking, the second-order scalar perturbations $\Phi^{(2)}$ and $\Psi^{(2)}$ are unequal: $\Phi^{(2)}\neq\Psi^{(2)}$.

Let us conclude this subsection by introducing the promised helpful combinations

[TABLE]

and presenting their explicit expressions without concealing anything, even despite the quite cumbersome form of some of them:

[TABLE]

[TABLE]

[TABLE]

$Q_{22}$ and $Q_{33}$ are analogous,

[TABLE]

where $\{x\leftrightarrow y\}$ stands for exactly the same terms with the occurring everywhere replacement of $x$ by $y$ and vice versa, $Q_{13}$ and $Q_{23}$ are analogous. The enumerated formulas have been derived, in particular, through the instrumentality of the formulas (3.10), (3.1) and equations from Section 2 including the following direct consequence of (2.2): $\mathcal{H}^{2}-\mathcal{H}^{\prime}=\kappa a^{2}\overline{\varepsilon}/2=\kappa\overline{\rho}c^{2}/(2a)$. From (3.31) and (3.32) we get

[TABLE]

[TABLE]

respectively. Hereinafter summation over repeated Greek subscripts is implied without superfluous decoding. It should be also noted that naturally $Q_{ik}=Q_{ki}$, in complete agreement with the symmetry inherent in the Einstein equations, which we write down in the very next subsection.

3.2 Scalar, Vector and Tensor Sectors

Substituting (3.5) and (3.25) into (3.4) with due account taken of (3.29), we immediately get:

$\bullet$  $00$-component:

[TABLE]

$\bullet$  $0\beta$-components:

[TABLE]

$\bullet$  $11$-component:

[TABLE]

$22$- and $33$-components are similar;

$\bullet$  $12$-component:

[TABLE]

$13$- and $23$-components are similar.

All these equations clearly demonstrate that the second-order scalar, vector and tensor perturbations (represented by $\Phi^{(2)}$, $\Psi^{(2)}$; ${\bf B}^{(2)}$; $h_{\alpha\beta}$, respectively) do not mix (Baumann et al., 2012) and are generated, in particular, by the quadratic combinations of the first-order scalar and vector perturbations $\Phi$; ${\bf B}$.

Now it is just the right time for the standard “scalar–vector–tensor” decomposition of $Q_{\alpha\beta}$:

[TABLE]

where $Q^{(0)}$ and $Q^{(\mathrm{S})}$ describe the scalar sector while ${\bf Q}^{(\mathrm{V})}$ and $Q^{(\mathrm{T})}_{\alpha\beta}$ describe the vector and tensor sectors, respectively, and satisfy the corresponding conditions:

[TABLE]

According to (3.40) and (3.41), the introduced functions $Q^{(0)}$, $Q^{(\mathrm{S})}$ and ${\bf Q}^{(\mathrm{V})}$ can be determined as solutions of the corresponding equations

[TABLE]

[TABLE]

then the remaining unknown functions $Q^{(\mathrm{T})}_{\alpha\beta}$ can be easily found from (3.40). They act as the sole sources of gravitational waves (see the rhs of the equation (3.47) for $h_{\alpha\beta}$ below).

Let us synchronously perform the standard “scalar–vector” decomposition of $Q_{0\alpha}$:

[TABLE]

where $Q^{(\parallel)}$ and ${\bf Q}^{(\perp)}$ denote the scalar and vector contributions, respectively. According to (3.44), they can be determined as solutions of the corresponding equations

[TABLE]

Finally, in order to abridge the notation, let us introduce the handy mixed-order quantities

[TABLE]

Then the Einstein equations (3.36)–(3.39) can be eventually rewritten as follows:

$\bullet$  tensor sector:

[TABLE]

$\bullet$  vector sector:

[TABLE]

[TABLE]

$\bullet$  scalar sector:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Substituting (3.52) into (3.51), we get

[TABLE]

Recalling (3.46) along with (2.8) and (2.9), we reduce (3.48), (3.50) and (3.54) to the equations

[TABLE]

[TABLE]

Thus, we have derived the “master” equations (3.47), (3.55) and (3.56) for the sought-for second-order cosmological perturbations. In the next subsection we show that the remaining “non-master” equations (3.49), (3.52) and (3.53) are satisfied automatically provided that one takes advantage of the equations of motion governing the particle dynamics.

3.3 Verification of Equations

Without going into detail, let us outline the proof that the scalar sector equations (3.52) and (3.53) containing temporal derivatives are really satisfied. In the first place, one finds a derivative of (3.54) with respect to $\eta$ and further expresses $\left(\Psi^{(12)}\right)^{\prime}$ from (3.52). As a result,

[TABLE]

In the second place, one substitutes the expression for $\left(\Psi^{(12)}\right)^{\prime}$ from (3.52) into (3.53). As a result,

[TABLE]

Therefore, it is enough to show that the rhs of (3.57) is really equal to the rhs of (3.58) since the lhs of these equations is the same, and then to prove either equation. We have successfully coped with both these onerous tasks. Equating the right-hand sides, after lengthy calculation one eventually arrives at an identity. The following auxiliary formulas should be used on the way:

[TABLE]

in the first-order approximation (Eingorn, 2016);

[TABLE]

within the adopted accuracy. The underlying equations of motion of the $n$-th particle have the form ${\bf\tilde{v}}^{\prime}_{n}=-\mathcal{H}{\bf\tilde{v}}_{n}-\nabla\Phi-\mathcal{H}{\bf B}$ (Eingorn, 2016).

When the desired identity is achieved, it is enough to prove, for instance, the correctness of (3.58). Now the accuracy of (3.59) is insufficient, and it is necessary to take advantage of the spacetime interval for the $n$-th particle:

[TABLE]

where the metric corrections are computed at the point ${\bf r}={\bf r}_{n}$ and, as usual, do not include the divergent contributions from the considered particle itself. For the sake of simplicity we can confine ourselves to those terms in (3.58) which are not quadratic in particle velocities, then the Lagrange equations of motion have the form

[TABLE]

Multiplication of (3.62) by $\rho_{n}$ with subsequent summation over $n$ gives

[TABLE]

where the last equation of the triplet (2.11) has been used to replace the summand $\delta\rho{\bf B}^{\prime}$ by $-2\mathcal{H}\delta\rho{\bf B}$. We have also dropped all terms which would import the third order of smallness in the Einstein equations. If one additionally omits the terms importing the second order, then (3.63) is reduced exactly to (3.59).

Being armed with (3.63), after exhausting calculation one turns (3.58) into an identity. Thus, both initial non-master scalar sector equations (3.52) and (3.53) are satisfied. The same applies to (3.49). Indeed, suffice it to demonstrate that

[TABLE]

Recalling (3.48), one can reduce (3.64) to the following equation:

[TABLE]

In the framework of the above-mentioned simplification, that is without products of velocities, substitution of (3.58) and (3.63) into (3.65) eventually leads to the desired identity. Thus, the initial non-master vector sector equation (3.49) is satisfied as well. Obviously, the same applies to the gauge conditions (3.2) and (3.3) since exactly the same gauge conditions hold true for the corresponding right-hand sides of (3.55) and (3.47).

3.4 Self-consistent Separation of Summands

In the previous subsection we have demonstrated that the functions $\Phi^{(2)}$, $\Psi^{(2)}$; ${\bf B}^{(2)}$; $h_{\alpha\beta}$ determined as solutions of (3.56), (3.55) and (3.47), respectively, satisfy all Einstein equations in the second-order approximation. The following relevant question arises: is the undertaken separation of the first- and second-order terms well-grounded and self-consistent? In other words, do we correctly and logically assign orders to different summands?

Of course, the answer is affirmative. As an illustrative example, let us single out two types of terms in (3.51), namely, those which are either present, or absent in the corresponding equation in the framework of the first-order approximation

[TABLE]

which is equivalent to (2.8) in view of the gauge condition (2.4) and the first equation of the triplet (2.11). The vanishing last term in the lhs of (3.66) is momentarily reinstated since it can be considered as being initially present in the corresponding $00$-component of the Einstein equations as a part of $\left(G_{0}^{0}\right)^{(1)}$ before applying the gauge condition (2.4).

We designate the first derivatives with respect to each comoving spatial coordinate and conformal time as $1/L$ and $1/\Upsilon$, respectively, as well as ascribe the orders of smallness $\epsilon$ and $\epsilon^{2}$ to the first- and second-order metric corrections. Then, for instance, $\triangle\Phi\sim\epsilon/L^{2}$ while $\left(\Psi^{(2)}\right)^{\prime}\sim\epsilon^{2}/\Upsilon$. As a result, taking into account the explicit expression (3.30) for $Q_{00}$, we have $6$ terms of the first type (present in (3.51) as well as in (3.66)), namely,

[TABLE]

and $9$ terms of the second type (present in (3.51), but absent in (3.66)):

[TABLE]

We distinguish between the coefficients $\kappa\overline{\rho}c^{2}/a$ and $\mathcal{H}^{2}$: they evolve synchronously during the matter-dominated stage of the Universe evolution, but asynchronously during the $\Lambda$-dominated stage. The essence of the perturbative computation lies in the fact that for each term of the second type in (3.68) there must exist a counterpart of the first type in (3.67), such that their ratio is of the order of smallness $\epsilon$. This is what we intend to confirm right now.

It can be easily seen that first five terms in (3.68), divided by the corresponding first five terms in (3.67), give precisely the order of smallness $\epsilon$. Obviously, the same applies to the sixth terms. Indeed, $\rho\tilde{v}^{2}\ll|\delta\rho|$ at arbitrary distances (Chisari & Zaldarriaga, 2011), and the helpful estimate $\tilde{v}^{2}\sim\Phi\delta\rho/\rho$ (see Baumann et al. (2012)) holds true. Hence, $\rho\tilde{v}^{2}/\delta\rho\sim\Phi\sim\epsilon$. Further, since $\kappa c^{2}\Xi/a\sim\Phi^{\prime}+\mathcal{H}\Phi$ (2.11) and $\kappa c^{2}\rho{\bf\tilde{v}}/a\sim\triangle{\bf B}-2\kappa\overline{\rho}c^{2}{\bf B}/a$ (2.9), the seventh term in (3.68) is reduced to a combination of $(\mathcal{H}/\Upsilon)\epsilon^{2}$ and $\mathcal{H}^{2}\epsilon^{2}$ while the eighth term is reduced to a combination of $\epsilon^{2}/L^{2}$ and $\left(\kappa\overline{\rho}c^{2}/a\right)\epsilon^{2}$. This quartet is already present in (3.68), hence, the seventh and eighth summands add nothing new. Similarly, the last term $\kappa^{2}c^{4}\Xi^{2}/a^{2}\sim(\Phi^{\prime}+\mathcal{H}\Phi)^{2}=\Phi^{\prime 2}+2\mathcal{H}\Phi\Phi^{\prime}+\mathcal{H}^{2}\Phi^{2}$. This is a combination of $\epsilon^{2}/\Upsilon^{2}$, $(\mathcal{H}/\Upsilon)\epsilon^{2}$ and $\mathcal{H}^{2}\epsilon^{2}$. Further, $\Phi^{\prime\prime}=-3\mathcal{H}\Phi^{\prime}-\left(2\mathcal{H}^{\prime}+\mathcal{H}^{2}\right)\Phi$ (2.11) and $\mathcal{H}^{\prime}=\mathcal{H}^{2}-\kappa\overline{\rho}c^{2}/(2a)$ (2.2), hence, in its turn, $\epsilon^{2}/\Upsilon^{2}$ may be treated as a combination of $(\mathcal{H}/\Upsilon)\epsilon^{2}$, $\mathcal{H}^{2}\epsilon^{2}$ and $\left(\kappa\overline{\rho}c^{2}/a\right)\epsilon^{2}$. Consequently, the last summand in (3.68) also adds nothing new to those terms which are already available in the collection.

Thus, we have shown that the elaborated perturbative scheme is valid. This scheme elegantly resolves the formidable challenge briefly discussed in the introductory part of (Clarkson & Umeh, 2011): at any cosmological scale for each summand in the equations for the second-order metric corrections there exists a much larger counterpart in the corresponding equations for the first-order metric corrections. Therefore, in particular, the situation when magnitudes of $\Phi^{(2)}$ and $\Phi$ are comparable is really improbable. Quite the contrary, the inequality $\left|\Phi^{(2)}\right|\ll|\Phi|$ may be expected to occur everywhere, as it certainly should be in the framework of a self-consistent perturbation theory.

3.5 Minkowski Background Limit

In this subsection, again for the sake of simplicity, we momentarily ignore all terms being quadratic in particle velocities and concentrate on the Minkowski background limit: the scale factor $a$ is now just a constant, $\mathcal{H}=0$, $\overline{\rho}=0$. Then, according to Eingorn (2016),

[TABLE]

[TABLE]

The sum of Newtonian potentials (3.69) is a solution of the standard Poisson equation

[TABLE]

At the same time, from the second equation in (3.56) and (3.30) we get

[TABLE]

where an evident relationship $2\left(\nabla\Phi\right)^{2}=\triangle\left(\Phi^{2}\right)-2\Phi\triangle\Phi$ has been used along with (3.71). Hence,

[TABLE]

After lengthy calculation, being based on (3.32), (3.33), (3.35) and (3.42), one also finds

[TABLE]

Substitution of (3.73) and (3.74) into the first equation in (3.56) gives

[TABLE]

As usual, the gravitational field produced by the $n$-th particle is excluded from the factors $\left.\Phi\right|_{{\bf r}={\bf r}_{n}}$ and $\left.(\nabla\Phi)\right|_{{\bf r}={\bf r}_{n}}$.

Let us compare the solutions (3.70) and (3.75) with the corresponding adapted expressions

[TABLE]

[TABLE]

which are equivalent to those from the textbook by Landau & Lifshitz (2000) (see the formulas (106.15) and (106.13) therein). Here we still ignore velocities squared as we arranged before. Of course, neither (3.70) coincides with (3.76), nor (3.75) coincides with (3.77). As pointed out by Eingorn (2016), the reason lies in the fact that our gauge conditions differ from those applied by Landau & Lifshitz (2000). Therefore, in order to reach agreement with this textbook, suffice it to find such a transformation of coordinates that would establish desired linkage. Apparently, it is enough to transform only the temporal coordinate: $\eta\mapsto\eta-A(\eta,{\bf r})$, then $\Phi^{(2)}\mapsto\Phi^{(2)}+A^{\prime}$ and ${\bf B}\mapsto{\bf B}+\nabla A$. Demanding that

[TABLE]

with the help of (3.70), (3.75), (3.76) and (3.77) we get

[TABLE]

[TABLE]

Action of $\nabla$ on both sides of (3.79) gives

[TABLE]

and exactly the same result follows also from (3.80). This incontestable fact ensures existence of the function $A(\eta,{\bf r})$ and, consequently, of the above-mentioned coordinate transformation. Thus, agreement with Landau & Lifshitz (2000) has been reached.

4 AVERAGING INITIATIVES ON THE EVE OF COSMOLOGICAL BACKREACTION ESTIMATION

In view of the predictably zero average values of the first-order metric corrections (Eingorn, 2016), the computation of the cosmological backreaction effects should be based on the second-order perturbation theory. Without pretending to an exhaustive study, let us perform the Euclidean averaging, or smoothing (Clarkson et al., 2011), of the $00$-component of Einstein equations (3.36), multiplied by $a^{2}/2$, and the sum of $11$-, $22$- and $33$-components (see (3.38)), multiplied by $\left(-a^{2}/6\right)$. We gather all terms containing $\overline{\Psi^{(2)}}$, $\overline{\Phi^{(2)}}$ and their temporal derivatives in the lhs, while the other averaged contributions are gathered in the rhs:

[TABLE]

[TABLE]

Here the overline indicates integrating over a comoving volume $\mathcal{V}$ and dividing by this volume in the limit of the infinite integration domain ($\mathcal{V}\rightarrow+\infty$). In addition, we have introduced the effective average energy density $\overline{\varepsilon}^{(\mathrm{II})}(\eta)$ and pressure $\overline{p}^{(\mathrm{II})}(\eta)$:

[TABLE]

[TABLE]

where the explicit expressions (3.30) for $Q_{00}$ and (3.35) for $Q_{\alpha\alpha}$ have been used. Replacing $\overline{(\nabla\Phi)^{2}}$ by $-\overline{\Phi\triangle\Phi}$ and expressing $\triangle\Phi$ and $\triangle{\bf B}$ from (2.8) and (2.9), respectively, we rewrite (4.3) and (4.4) in the more compact form:

[TABLE]

[TABLE]

Expressing $\overline{\varepsilon}^{(\mathrm{II})}$ and $\overline{p}^{(\mathrm{II})}$ from (4.1) and (4.2), one can easily verify that these functions satisfy the standard conservation equation

[TABLE]

as it certainly should be. Hence, the expressions (4.5) and (4.6) for the same functions must automatically satisfy this equation as well. This can be verified through the instrumentality of the formulas (3.59) and (3.60) as well as equations from Section 2.

It is worth mentioning that if one keeps in the rhs of (4.5) and (4.6) only first two terms, which dominate at sufficiently small scales, and makes use of the relationship $\overline{\rho\Phi}=-2\overline{\rho\tilde{v}^{2}}$, which holds true for the virialized regions, then $\overline{p}^{(\mathrm{II})}\rightarrow 0$ while $\overline{\varepsilon}^{(\mathrm{II})}\rightarrow-\left[c^{2}/\left(2a^{3}\right)\right]\overline{\rho\tilde{v}^{2}}\sim 1/a^{3}$. Thus, at virialized scales the effective pressure $\overline{p}^{(\mathrm{II})}$ vanishes while the non-vanishing effective average energy density $\overline{\varepsilon}^{(\mathrm{II})}$ brings to a small time-independent renormalization of the corresponding background quantity $\overline{\varepsilon}$, in full accord with Baumann et al. (2012) (see additionally Wetterich (2003) for earlier theoretical efforts and a “cosmic virial theorem”). The interpretation of the simulation outputs by Adamek et al. (2015) also suggests that “stable clustering” (implying virialized nonlinear structures) razes backreaction from the cosmological battlefield. The underlying perturbative scheme advocated by Adamek et al. (2015) is compared with ours by Eingorn (2016). It is necessary to mention that this purely numerical scheme is characterized by the first order accuracy for large enough distances and the second order accuracy for sufficiently small distances, while the approach advocated in the current paper is characterized by the second order accuracy everywhere and is fully analytical at least with respect to the first-order cosmological perturbations (2.5); (2.6) and the sources (3.30)–(3.33) of the second-order ones.

It is noteworthy as well that the velocity-dependent summands can be easily distinguished from the velocity-independent ones in (4.5) and (4.6), and there are only two types of contributions, which do not contain particle velocities: $\sim\overline{\rho\Phi_{0}}$ and $\sim\overline{\Phi_{0}^{2}}$. Here $\Phi_{0}$ denotes the velocity-independent part of the first-order scalar perturbation $\Phi$ (2.5), that is the sum of Yukawa potentials with the same interaction range $\lambda$ (up to an additive constant $1/3$):

[TABLE]

For illustration purposes, we compute both average quantities $\overline{\rho\Phi_{0}}$ and $\overline{\Phi_{0}^{2}}$ analytically:

[TABLE]

[TABLE]

It presents no difficulty to receive evidence that both these expressions tend to zero in the homogeneous mass distribution limit ($\sum\sum m_{n}m_{k}\rightarrow\overline{\rho}^{2}\iint d{\bf r}_{n}d{\bf r}_{k}$), as it certainly should be since $\Phi_{0}=0$ at any point in this test limit (Eingorn, 2016).

Reverting to (4.5) and (4.6), we emphasize that collections of terms in the right-hand sides may assist in the cosmological backreaction estimation. We formulate the following quite feasible two-stage plan:

$\bullet$  the launch of a new generation of cosmological $N$-body simulations based on the formalism developed by Eingorn (2016) (see the equations of motion (3.6) therein);

$\bullet$  the use of outputs of these simulations for the estimation of the effective average energy density $\overline{\varepsilon}^{(\mathrm{II})}$ and pressure $\overline{p}^{(\mathrm{II})}$, and the subsequent comparison with the background quantity $\overline{\varepsilon}$.

If the underlying inequalities $\left|\overline{\varepsilon}^{(\mathrm{II})}\right|\ll\overline{\varepsilon}$ and $\left|\overline{p}^{(\mathrm{II})}\right|\ll\overline{\varepsilon}$ become doubtful at any moment during the matter-dominated or $\Lambda$-dominated stages of the Universe evolution, then this fact may serve as a sure sign of backreaction significance and inappropriateness of the FLRW metric (2.1) and Friedmann equations (2.2), with inevitable grave consequences. At the same time, if the inequalities being tested seem always unquestionable, this result by itself does not necessarily mean that backreaction is insignificant, for the simple reason that we still rely on the initial assumption of the FLRW background existence and actually have no predictive power beyond. Nevertheless, if this key assumption is valid and $\overline{\varepsilon}$ is really much greater than $\left|\overline{\varepsilon}^{(\mathrm{II})}\right|$ and $\left|\overline{p}^{(\mathrm{II})}\right|$, then we can formally celebrate the preliminary success of the elaborated perturbative approach and add a third stage relying on (3.55) and (3.56):

$\bullet$  the estimation of $\Phi^{(2)}$, $\Psi^{(2)}$ and ${\bf B}^{(2)}$, and the subsequent comparison with $\Phi$ and ${\bf B}$, respectively.

The proposed plan exploits Yukawa gravity resulting from GR (Eingorn, 2016, 2017) and therefore possesses a definite advantage over generally accepted simulations exploiting Newtonian gravity (Peebles, 1980; Springel, 2005). This advantage is clearly explained by Räsänen (2010, 2011): in the framework of the so-called Newtonian cosmology, as opposed to GR, the backreaction effects are reduced to boundary terms vanishing in the case of the standard periodic boundary conditions. It is stressed by Räsänen (2010, 2011) that Newtonian gravity does not represent the weak field limit of GR. We emphasize that the formulation of this crucial limit is no longer an open issue: the corresponding cosmological perturbation theory incorporating nonlinear density contrasts has been developed by Eingorn (2016) and extended in the current paper. As pointed out by Eingorn (2016), Yukawa screening of the interparticle attraction may be treated as a relativistic effect. This is especially important in view of the fact that backreaction significance is inseparably linked with non-Newtonian gravitational physics (Räsänen, 2010; Buchert & Räsänen, 2012).

Of course, the first-order scheme elaborated by Eingorn (2016) is insufficient for $N$-body simulations taking into consideration the second-order metric corrections in the equations of motion, and the formalism developed in the current paper should be employed instead. This is an indispensable complication of the aforesaid simplified plan from the point of view of modern precision cosmology. There are formidable obstacles on the way of solving the derived equations (3.55) and (3.56) analytically, however, owing to their linearity with respect to the sought-for functions, it presents no difficulty in principle to solve these equations numerically. Furthermore, it would be very interesting to study propagation of gravitational waves in the inhomogeneous Universe, governed by (3.47). Having at our disposal the determined first- and second-order cosmological perturbations, we can in principle introduce and calculate the “averaged Hubble rate” as well as the kinematical and dynamical backreactions (Behrend et al., 2008; Brown et al., 2009, 2013). Moreover, the thorough analysis of the “effective equation of state” (that is the relationship between $\overline{\varepsilon}$, $\overline{\varepsilon}^{(\mathrm{II})}$ and $\overline{p}^{(\mathrm{II})}$) for different cosmological epochs becomes feasible (see additionally Li & Schwarz (2007, 2008); Schwarz (2010); Bose & Majumdar (2013) for related theoretical efforts). It is also noteworthy that all formulas, where the discreteness of matter is not manifestly specified (for example, by means of the delta-functions), are valid for continuous medium as well. Thus, the transition to the hydrodynamical description is simple.

Implementation of both aforesaid “simplified” and “complicated” plans is extremely promising, but oversteps the limits of our prolonged narration.

5 CONCLUSION

Summarizing, we recollect the main results achieved in the current paper within the conventional $\Lambda$CDM model:

$\bullet$  the equations (3.56), (3.55) and (3.47) for the second-order scalar ($\Phi^{(2)}$, $\Psi^{(2)}$), vector (${\bf B}^{(2)}$) and tensor ($h_{\alpha\beta}$) cosmological perturbations have been derived. These equations are suitable at all spatial scales (naturally, except for the regions of strong gravitational fields in immediate proximity to such their generators as black holes or neutron stars) and permit of nonlinear density contrasts;

$\bullet$  the Helmholtz equations (3.55) and (3.56) incorporate exactly the same Yukawa interaction ranges as their corresponding counterparts (2.9) and (2.8) determining the first-order metric corrections;

$\bullet$  the constructed scheme passes three important corroborative tests:

1. we have verified that all Einstein equations are satisfied within the adopted accuracy along with the gauge conditions (subsection 3.3);

2. we have confirmed the self-consistency of order assignments and the related expectation that the first-order metric corrections dominate over the second-order ones everywhere (subsection 3.4);

3. in the Minkowski background limit the linkage with the textbook material has been established (subsection 3.5);

$\bullet$  the highway to investigate the cosmological backreaction effects beyond Newtonian gravitational physics has been outlined. The proposed relativistic simulations of the cosmic structure growth accompanied by the investigations of propagation of light and gravitational waves in the inhomogeneous Universe would definitely assist in deepening and testing our knowledge of the spacetime and world’s filling material including dark ingredients.

ACKNOWLEDGEMENTS

We are grateful to the anonymous Referee for the valuable comments. The work of R. Brilenkov was partially supported by the EMJMD Student Scholarship from the Erasmus + : Erasmus Mundus Joint Master Degree programme AstroMundus in Astrophysics. The work of M. Eingorn was partially supported by NSF CREST award HRD-1345219 and NASA grant NNX09AV07A.