Analytic expressions for the second-order scalar perturbations in the $\Lambda$CDM Universe within the cosmic screening approach

TL;DR

This paper derives analytic expressions for second-order scalar perturbations in the $$CDM universe using the cosmic screening approach, highlighting the Yukawa screening effect and its implications for precision cosmology.

Contribution

It provides the first analytic formulas for second-order scalar perturbations in $$CDM within the cosmic screening framework, including velocity-dependent effects.

Findings

Yukawa screening effect appears in second-order perturbations.

Formulas reduce to post-Newtonian expressions at small distances.

Second-order corrections impact interpretation of cosmological observations.

Abstract

We study the second-order scalar perturbations in the conventional $\Lambda$CDM Universe within the cosmic screening approach. The analytic expressions for both the velocity-independent and velocity-dependent second-order scalar perturbations are derived. We demonstrate how the Yukawa screening effect, which is inherent in the first-order metric corrections, manifests itself in the second-order ones. It is shown that the obtained formulas for the second-order perturbations are reduced to the known post-Newtonian expressions at distances much smaller than the Yukawa screening length. In the era of precision cosmology, these analytic formulas play an important role since the second-order metric corrections may affect the interpretation of observational data (e.g., the luminosity-redshift relation, gravitational lensing, baryon acoustic oscillations).

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 ,$$

$${\mathcal{H}}^{\prime }-{\mathcal{H}}^2=-\frac{\kappa \overline{\rho }{c}^2}{2a}.$$

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

$$\delta \rho =\rho -\overline{\rho }.$$

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

$$\Delta \Phi -\frac{a^2}{{\lambda }^2}\Phi =\frac{\kappa c^2}{2a}\delta \rho ,$$

$$\lambda \equiv {\left[\frac{3\kappa \overline{\rho }{c}^2}{2a^3}\right]}^{-1/2}$$

$$\Phi (\eta ,\mathbf{r})=\frac{1}{3}-\frac{\kappa c^2}{8\pi a}\sum _n\frac{m_n}{|\mathbf{r}-{\mathbf{r}}_n|}{e}^{-\mu |\mathbf{r}-{\mathbf{r}}_n|},$$

$$\mu =\frac{a}{\lambda }={\left[\frac{3\kappa \overline{\rho }{c}^2}{2a}\right]}^{1/2}.$$

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

$$\Delta {\Psi }^{(2)}-\frac{a^2}{{\lambda }^2}{\Psi }^{(2)}$$

$${\Phi }^{(2)}-{\Psi }^{(2)}$$

$$Q_{00}$$

$$\Delta Q^{(\parallel )}$$

$$\Delta \Delta Q^{(S)}$$

$$Q_{\alpha \alpha }$$

$$Q_{\alpha \beta }$$

$$Q_{\beta \beta }$$

$$Q^{(\parallel )}=\frac{5}{a^2}\mathcal{H}{\Phi }^2.$$

$$Q_{00}=-\left(\frac{9\kappa \bar{\rho }{c}^2}{2a^3}+\frac{15}{a^2}{\mathcal{H}}^2\right){\Phi }^2-\frac{\kappa c^2}{a^3}\Phi \delta \rho -\frac{3}{a^2}(\nabla \Phi {)}^2.$$

$$\Delta {\Psi }^{(2)}-\frac{3\kappa \bar{\rho }{c}^2}{2a}{\Psi }^{(2)}=-\frac{9\kappa \bar{\rho }{c}^2}{4a}{\Phi }^2-\frac{\kappa c^2}{2a}\Phi \delta \rho -\frac{3}{2}(\nabla \Phi {)}^2$$

$$(\nabla \Phi {)}^2=\frac{1}{2}\Delta \left({\Phi }^2\right)-\Phi \Delta \Phi .$$

$$\chi ={\Psi }^{(2)}+\frac{3}{4}{\Phi }^2,$$

$$\Delta \chi -\frac{3\kappa \bar{\rho }{c}^2}{2a}\chi =-\frac{9\kappa \bar{\rho }{c}^2}{8a}{\Phi }^2+\frac{\kappa c^2}{4a}\Phi \delta \rho .$$

$$\Delta \chi -\frac{3\kappa \bar{\rho }{c}^2}{2a}\chi =\frac{\kappa c^2}{4a}\Phi \delta \rho .$$

$$G_{\mathrm{H}}(\text{r},{\text{r}}^{\prime })=-\frac{1}{4\pi }\frac{e^{-a|\text{r}-{\text{r}}^{\prime }|/\lambda }}{|\text{r}-{\text{r}}^{\prime }|}.$$

$$\chi$$

$${\Psi }^{(2)}$$

$$Q_{\alpha \alpha }=-\frac{15}{a^2}{\mathcal{H}}^2{\Phi }^2-\frac{4\kappa c^2}{a^3}\Phi \delta \rho -\frac{7}{a^2}(\nabla \Phi {)}^2,$$

$$\frac{{\partial }^2{Q}_{\alpha \beta }}{\partial x^{\alpha }\partial x^{\beta }}=-\frac{\kappa c^2}{a^3}\nabla \left(\delta \rho \nabla \Phi \right)-\frac{5}{a^2}{\mathcal{H}}^{\prime }\Delta ({\Phi }^2).$$

$$\Delta \Delta Q^{(S)}=-\frac{1}{2}\Delta Q_{\alpha \alpha }-\frac{15}{2a^2}{\mathcal{H}}^{\prime }\Delta ({\Phi }^2)+\frac{3\kappa \bar{\rho }{c}^2}{2a^3}\Delta \Phi -\frac{3\kappa c^2}{2a^3}\nabla \left(\rho \nabla \Phi \right).$$

$$f=\frac{1}{4\pi }\sum _k{m}_k\frac{1}{|\text{r}-{\text{r}}_k{|}^3}(\text{r}-{\text{r}}_k)\nabla \Phi (\eta ,{\text{r}}_k),$$

Full text

1 Introduction

According to the cosmological principle [1, 2], our Universe is isotropic and homogeneous at large enough scales. This follows from the natural assumption that the laws of physics should be the same wherever in the Universe we are. Starting from a certain scale, the distribution of inhomogeneities (e.g., galaxies and groups of galaxies) should be statistically homogeneous. As a result, such statistically homogeneous Universe has the Friedmann-Lema$\mathrm{\hat{\i}}$tre-Robertson-Walker (FLRW) metric, and its dynamics is described by the Friedmann equations. This is the zero-order/background approach. Obviously, inhomogeneities disturb this background resulting in first- and higher-order perturbations. They play a crucial role in the investigation of the large scale structure formation. Evidently (by definition!), the average values of the first-order perturbations should be equal to zero (see, e.g., [3]). On the other hand, the average values of the second-order perturbations are different from zero and can affect the background spacetime and matter. This effect is called backreaction (see, e.g., [4, 5, 6, 7, 8, 9] and references therein). It is important to determine how strong the backreaction is and to what extent we may use the standard FLRW metric as a background one. As an example, the backreaction may affect the baryon acoustic oscillations [5, 10]. The second-order perturbations also contribute to the luminosity-redshift relation [11] and gravitational lensing [12, 13].

Within the cosmic screening approach, the theory of the first-order perturbations was developed in the papers [14, 15, 16, 17, 18]. In its framework, inhomogeneities in the Universe (e.g., galaxies and their groups) are considered as point-like gravitating masses. These masses disturb the background spacetime and matter. For example, the energy density fluctuation reads $\delta\varepsilon\approx(c^{2}/a^{3})\delta\rho+(3\bar{\rho}c^{2}/a^{3})\Phi$, where $c$ denotes the speed of light, $a$ is the scale factor of the Universe, and we singled out the gravitational potential $\Phi$. Consequently, the 00-component of the Einstein equation is reduced to the Helmholtz-type equation (rather than the Poisson one) for the gravitational potential, where the comoving mass density fluctuation $\delta\rho=\rho-\bar{\rho}$ is the source (see Eq. (2.7) below). As a result, the first-order scalar perturbation $\Phi$ is characterized by exponential Yukawa screening at large cosmological scales. A similar effect takes place for the first-order vector perturbation ${\bf B}$. The cosmic background, namely, the nonzero value of the average mass density $\bar{\rho}$, is responsible for this effect. It is worth mentioning that in our approach we assume that peculiar velocities of gravitating masses are much less than the speed of light and the gravitational field produced by inhomogeneities is weak (i.e. we work in the weak field approximation). However, we do not demand the smallness of the mass density contrast $\delta\rho/\overline{\rho}$. Therefore, the given formalism can be used at all cosmic scales (below and above the Hubble scale).

In the framework of this formalism, the theory of the second-order perturbations was developed in [19]. It is worth noting that the alternative perturbative method was considered, e.g., in [20]. There, the authors came to the conclusion that the backreaction effects are at a few percent level at most. In their analysis, in particular, they disregard the term $\sim\Phi^{2}$. In our present paper, to get the analytic solutions for the second-order perturbations, we follow the similar method, i.e. we also drop such terms as $\sim\Phi^{2}$ from the sources of the analyzed perturbations. However, this is the only similarity with [20]. In other respects our approach is completely different. In [20] the authors actually mix the first- and second-order perturbations, so their equations are nonlinear. In contrast, in our approach there is no mixture of the orders, so the equations are linear with respect to the corresponding metric corrections. Thus, in general, we can arrive at different conclusions regarding the backreaction (at least to the extent that may be important for precision cosmology). Moreover, the interpretation of our results is clearer since the orders are not mixed and it is much easier to deal with linear equations. As is generally known, from the point of view of numerical calculus, linear equations possess an unquestionable advantage over nonlinear ones, which consists in higher computational accuracy and saving of computing time.

In [19] the system of equations for the second-order scalar, vector and tensor perturbations was derived. This system has rather complicated form. However, it is remarkable that these perturbations do not mix, and we can study them separately. In the present paper we study the second-order scalar perturbations. Even in this case the corresponding equations are still very complicated. So, we restrict ourselves to the analysis of the second-order perturbations generated by those sources which are dominant at sufficiently small spatial scales, where the mass density contrast is typically large (see also [20]). This method enables us to solve the reduced equations exactly. In the large-scale spatial regions the linear relativistic perturbation theory works very well. Hence, the above mentioned simplification is not significant. As a result, we find analytic expressions for the velocity-independent as well as velocity-dependent scalar perturbations.

The paper is structured as follows. In section 2 we briefly describe the background model and the first-order scalar perturbations. Then, we present the equations for the second-order scalar perturbations with velocity-independent sources. These equations are solved in section 3. The solutions for the second-order velocity-dependent scalar perturbations are derived in section 4. In section 5 we study the Newtonian and post-Newtonian approximations for the found analytic solutions. Here we demonstrate that both the velocity-independent and velocity-dependent expressions are reduced to the known post-Newtonian formulas. In the concluding section 6 we summarize and discuss the obtained results. In Appendix A we collect formulas which we use to solve the equations for the second-order scalar perturbations.

2 From the background model to the second-order scalar perturbations: basic equations

2.1 Background model

We start with the unperturbed FLRW metric

[TABLE]

where $a(\eta)$ is the scale factor, $\eta$ is the conformal time, and $x^{\alpha}$, $\alpha=1,2,3$, represent the comoving coordinates. It is supposed that the spatial curvature is zero. Let us write down the corresponding Friedmann equations in the framework of the $\Lambda$CDM model:

[TABLE]

where $\mathcal{H}\equiv a^{\prime}/a\equiv(da/d\eta)/a$, with the prime standing for the derivative with respect to $\eta$, while $\kappa\equiv 8\pi G_{N}/c^{4}$ (with $c$ denoting the speed of light and $G_{N}$ being the Newtonian gravitational constant). In addition, $\varepsilon$ stands for the energy density of nonrelativistic pressureless matter, $\Lambda$ denotes the cosmological constant, and the overline indicates the average. Obviously, the average energy density is defined by the constant average comoving mass density $\overline{\rho}$ as follows: $\overline{\varepsilon}=\overline{\rho}c^{2}/a^{3}$. From Eqs. (2.2) we can easily get a useful auxiliary equation

[TABLE]

2.2 First-order scalar perturbations

The described above background Universe is perturbed by inhomogeneities in the form of discrete point-like masses with mass density

[TABLE]

where $\mathbf{r}\equiv\left(x^{1},x^{2},x^{3}\right)$ is the comoving radius-vector of the observation point, and $\mathbf{r}_{n}\equiv\left(x_{n}^{1},x_{n}^{2},x_{n}^{3}\right)$ is the comoving radius-vector of the $n$-th particle. These massive particles may represent galaxies and their groups. The mass density fluctuation is

[TABLE]

It is important to note that we do not assume the smallness of the mass density contrast, i.e. $\delta\rho/\overline{\rho}$ can be much larger than unity. Hence, our scheme is valid at both superhorizon and subhorizon scales. The inhomogeneities result in scalar perturbations of the metric (2.1). In the conformal Newtonian gauge and in the first-order approximation, the perturbed metric is [21, 22, 23, 24]

[TABLE]

where the first-order scalar perturbation $\Phi(\eta,\mathbf{r})$ satisfies the inequality $|\Phi|\ll 1$. This means that we work in the weak gravitational field limit. It is well known that, e.g., in the vicinity of galaxies the mass density contrast is much larger than unity, however the gravitational field is weak. This fact allows expanding the components of the Einstein equation into series in metric corrections. Additionally, we assume that the particle peculiar velocities are much less than the speed of light: $|\mathbf{v}_{n}|=|cd\mathbf{r}_{n}/d\eta|\ll c$. For example, the today’s typical values are $(250\div 500)$ km/s. As demonstrated in [14], for such nonrelativistic velocities the contribution of the velocity-dependent part into the total expression for the first-order scalar perturbation $\Phi$ is negligible. In this case the gravitational potential $\Phi$ satisfies the following Helmholtz-type equation [14]:

[TABLE]

where $\Delta\equiv\delta^{\alpha\beta}\frac{\partial^{2}}{\partial x^{\alpha}\partial x^{\beta}}$ is the Laplace operator in comoving coordinates. The time-dependent parameter

[TABLE]

defines the characteristic Yukawa screening length of the gravitational interaction.

The solution of Eq. (2.7) is [14]

[TABLE]

where

[TABLE]

2.3 Second-order velocity-independent scalar perturbations

Let us turn now to the second-order scalar perturbations $\Phi^{(2)}(\eta,\mathbf{r})$ and $\Psi^{(2)}(\eta,\mathbf{r})$. The corresponding metric reads

[TABLE]

The main aim of the present paper consists in determination of these perturbations. According to [19], the functions $\Phi^{(2)}(\eta,\mathbf{r})$ and $\Psi^{(2)}(\eta,\mathbf{r})$ satisfy the following system of master equations:

[TABLE]

where we momentarily keep only the velocity-independent parts of the sources:

[TABLE]

[TABLE]

where $(\nabla\Phi)^{2}\equiv\nabla\Phi\nabla\Phi=\delta^{\alpha\beta}\frac{\partial\Phi}{\partial x^{\alpha}}\frac{\partial\Phi}{\partial x^{\beta}}$. These expressions follow from Eqs. (3.30), (3.45), (3.34), (3.42), (3.35), (3.33) and (3.32) in [19] where we disregard the peculiar velocities and vector perturbations. Now, we should solve these equations.

3 Analytic solutions for the second-order velocity-independent potentials $\Phi^{(2)}$ and $\Psi^{(2)}$

3.1 Potential $\Psi^{(2)}(\eta,\mathbf{r})$

The potential $\Psi^{(2)}$ is the solution of Eq. (2.12). To solve this equation, we should determine its right-hand side. First, from Eq. (2.15) we obtain the function $Q^{(\parallel)}$:

[TABLE]

Taking into account Eq. (2.7) for $\Delta\Phi$, the function $Q_{00}$ (2.14) can be rewritten in the form

[TABLE]

Then, Eq. (2.12) takes the form

[TABLE]

where we used Eq. (2.7) and the auxiliary equality

[TABLE]

It makes sense to define a new function

[TABLE]

which satisfies the equation

[TABLE]

To solve this equation analytically, we resort to a supplementary simplification. Namely, we concentrate on those sources of the second-order perturbations, which dominate at sufficiently small distances where the mass density contrast is typically large. In the case of Eq. (3.6) this means dropping the term $\sim\Phi^{2}$ while keeping the term $\sim\Phi\delta\rho$ in the right-hand side. Such a simplification implies failing to take into account all sources at large enough distances, but this failure is insignificant since the disregarded sources (such as the term $\sim\Phi^{2}$) are much less than the corresponding first-order ones in the considered large-scale spatial region, where the linear relativistic perturbation theory works very well (see the argumentation in [20] and Refs. therein). Returning to Eq. (3.6), we have

[TABLE]

This is the Helmholtz-type equation with the Green’s function

[TABLE]

Therefore, we look for its solution in the form

[TABLE]

where the integrals $I_{1}$ and $I_{2}$ are given by the formulas (A.1) and (A.2), respectively. The prime over the double sum indicates that the summation indices must not coincide. Substituting these integrals and taking into account Eq. (3.5), we finally get

[TABLE]

In spite of the presence of the first-order term in the right-hand side, this function in total is of the second order. We clearly demonstrate it for the case of the small-scale limit in section 5 (see the formula (5.8)).

3.2 Potential $\Phi^{(2)}(\eta,\mathbf{r})$

According to Eq. (2.13), to get the potential $\Phi^{(2)}$, we need to determine the function $Q^{(S)}$. This function satisfies Eq. (2.16) where $Q_{\alpha\alpha}$ and $Q_{\alpha\beta}$ are defined by Eqs. (2.17)-(2.19). The function $Q_{\alpha\alpha}$ can be rewritten as follows:

[TABLE]

where we used Eq. (2.7). After lengthy calculations one can also derive

[TABLE]

Therefore, the function $Q^{(S)}$ satisfies the following equation:

[TABLE]

To solve this equation, we introduce a new function

[TABLE]

which satisfies the equation

[TABLE]

Now, applying the inverse Laplace operator $\Delta^{-1}$ to Eq. (3.13), we get

[TABLE]

where the function $Q_{\alpha\alpha}$ is given by (3.11). With the help of Eqs. (2.7) and (3.4) as well as the formula $1/{|\text{{r}}-\text{{r}}_{k}|}=\frac{1}{2}\Delta|\text{{r}}-\text{{r}}_{k}|$ (see, e.g., $\S$ 106 in [25]) we finally arrive at the equation

[TABLE]

Using the same reasoning as for Eq. (3.6), we drop the term $\sim\Phi^{2}$ in the right-hand side. The solution of the resulting equation can be found with the help of the Laplace operator Green’s function

[TABLE]

Let us introduce two new functions $F_{1}$ and $F_{2}$:

[TABLE]

Therefore, for $Q^{(S)}$ we get

[TABLE]

where

[TABLE]

and

[TABLE]

The function $F_{2}$ can be also expressed as follows:

[TABLE]

where we used the integral (A.3).

Finally, the potential $\Phi^{(2)}$ reads:

[TABLE]

where $\Psi^{(2)}$ and $Q^{(S)}$ are given by Eqs. (3.10) and (3.21), respectively.

4 Second-order velocity-dependent potentials $\Phi^{(2)}_{(v)}$ and $\Psi^{(2)}_{(v)}$

Let us turn now to the analytic expressions for the second-order velocity-dependent potentials. In the framework of the cosmic screening approach, the velocity-independent and velocity-dependent sources of the potentials enter into the corresponding equations additively. For example, the right-hand sides of Eqs. (2.12) and (2.13) (Eqs. (3.56) in [19]) clearly demonstrate it. Therefore, one can find the respective solutions separately. The velocity-independent parts of the second-order potentials were obtained in the previous section. To derive these parts, we focused on the gravitational field sources prevailing in small enough regions with considerable density contrasts. For these regions the velocity-dependent sources of the second-order perturbations are substantially simplified as well (see section 3.4 in [19], devoted to comparison of quantities of different orders, and also [10, 20] for the similar reasoning). For instance, the term $\sim\rho{\bf\tilde{v}}{\bf B}$ (the last one in the total expression for $Q_{00}$, given by Eq. (3.30) in [19]) is negligible as compared to the term $\sim\rho\tilde{v}^{2}$ (the first one in the same formula) in view of the fact that ${\bf B}$ is of the order of ${\bf\tilde{v}}\Phi$ at small scales. Consequently, the ratio $\left(\rho{\bf\tilde{v}}{\bf B}\right)/\left(\rho\tilde{v}^{2}\right)$ is of the order of $|\Phi|\ll 1$. Adhering to the considered supplementary simplification, from Eqs. (3.30)-(3.35) in [19] we get

[TABLE]

where $\tilde{v}^{\beta}\equiv dx^{\beta}/d\eta$ is connected with the physical peculiar velocity $v^{\beta}\equiv adx^{\beta}/dt$ as follows: $v^{\beta}=c\tilde{v}^{\beta}$ (in view of the relationship $ad\eta=cdt$ defining the synchronous time $t$). In addition, $\tilde{v}^{2}\equiv\delta_{\alpha\beta}\tilde{v}^{\alpha}\tilde{v}^{\beta}$. Here and in what follows, both the superscript “$(v)$” and the subscript “$(v)$” indicate that the corresponding quantity depends on velocities of massive particles. The function $Q^{(\parallel)}$ is defined by Eq. (3.45) in [19]. Taking into account (4.4), we get

[TABLE]

The function $Q^{(S)}$ satisfies Eq. (2.16) where $Q_{\alpha\beta}$ and $Q_{\alpha\alpha}$ are given now by (4.3) and (4.5), respectively, with the mass density (2.4). Then, the solution of (2.16) can be written in the form

[TABLE]

The inverse Laplacians in this equation are easily solvable:

[TABLE]

Therefore,

[TABLE]

Now we can derive expressions for the velocity-dependent second-order potentials. The potential $\Psi^{(2)}$ satisfies Eq. (2.12) where we should substitute the functions (4.1) and (4.6). Then, this equation reads

[TABLE]

It has the following solution:

[TABLE]

where $\mu$ is given by (2.10). The potential $\Phi^{(2)}$ satisfies Eq. (2.13) where we should substitute (4.9) and (4.11):

[TABLE]

5 Newtonian and post-Newtonian cosmological approximations

In the present section we consider the derived above formulas at distances much smaller than the screening length: $\mu r=ar/\lambda=r_{\mathrm{ph}}/\lambda\ll 1$, where $r_{\mathrm{ph}}=ar$ is the physical distance. Thereby we analyze the Newtonian and post-Newtonian approximations. We call them cosmological since the obtained expressions depend on the scale factor $a$. To study these limits, we consider an auxiliary model. In this model a sphere of comoving radius $R$ contains $N$ discrete particles. Outside this sphere, the rest of the Universe is uniformly filled with matter with the constant comoving mass density $\overline{\rho}$. For such a geometrical configuration, we will get the first- and second-order perturbations in a point with radius-vector ${\bf r}$ inside the sphere: $r<R$. Additionally, we assume that the physical radius of the sphere is much less than the screening length: $aR/\lambda=\mu R\ll 1$.

5.1 Gravitational potential

It is well known (see, e.g., [25]) that the first-order perturbation $\Phi(\eta,\mathbf{r})$ corresponds to the gravitational potential. In the cosmic screening approach this function is given by Eq. (2.9). In the case of the described above model, the function (2.9) takes the form

[TABLE]

where the sum is taken over all discrete masses inside the sphere (i.e. $r_{n}<R$). Since we consider the point inside the sphere ($r<R$), the integral in (5.1) coincides with $\mathcal{I}_{1}(r<R)$ given by (A.5). Therefore,

[TABLE]

In the limit $\mu R\to 0\ \Rightarrow\ \mu r\to 0,\ \mu r_{n}\to 0$ we obtain

[TABLE]

where $\Phi_{N}$ is the Newtonian potential at the position $\mathbf{r}$ inside the sphere, produced by all $N$ discrete masses. It is worth reminding that $\kappa c^{2}/(8\pi)=G_{N}/c^{2}$ and the physical distance $r_{\mathrm{ph}}=ar$. The formula (5.3) clearly demonstrates that for the considered model and in the given approximation the gravitational potential is determined by the particles from the nearest environment, and the term $1/3$ is exactly compensated by the contribution of an infinite number of remote particles. As a result, the discussed expression is truly of the first order of smallness.

5.2 Velocity-independent second-order perturbations

Let us study now the second-order perturbation $\Psi^{(2)}(\eta,\mathbf{r})$ given by the formula (3.10). According to the previous subsection, in the Newtonian limit $\Phi^{2}\to\Phi_{N}^{2}$. Therefore, we only need to investigate the function $\chi$. For the considered model with $N$ discrete particles inside the sphere and uniformly distributed matter outside the sphere, the function $\chi$ takes the form

[TABLE]

Here the introduced term $S_{N}$ incorporates pure sums. The integrals $\mathcal{I}_{1}(r<R)$, $\mathcal{I}_{2}(r<R)$ and $\mathcal{I}_{3}(r<R)$ are given by Eqs. (A.5), (A.7) and (A.9), respectively. In the last line of (5.4) $\mathbf{r}_{k}$ are the radius-vectors of discrete particles inside the sphere. To evaluate the expression in this line, we suppose that the discrete masses are concentrated in the central part of the sphere, i.e. $r_{k}\ll R\,\Rightarrow\,r_{k}\ll r^{\prime}$. Hence,

[TABLE]

where the integrals $\mathcal{I}_{4}$ and $\mathcal{I}_{5}(r<R)$ are given by Eqs. (A.10) and (A.12), respectively. Substituting this expression into (5.4) and taking into account the integrals (A.5), (A.7), (A.9), (A.10) and (A.12), we get

[TABLE]

In the limit $\mu R\to 0\ \Rightarrow\ \mu r\to 0,\ \mu r_{k}\to 0$ we obtain

[TABLE]

Therefore, taking into account Eqs. (3.10) and (5.3), the potential $\Psi^{(2)}(\eta,\mathbf{r})$ tends to

[TABLE]

in full agreement with the formula (3.73) in [19]. This expression clearly demonstrates that $\Psi^{(2)}$ is really of the second order of smallness.

To investigate the same limit of the second-order potential $\Phi^{(2)}(\eta,\mathbf{r})$, we need to substitute (5.8) as well as the limit of $Q^{(S)}$ into Eq. (3.25). One can show that

[TABLE]

Consequently, we find

[TABLE]

The expressions (5.9) and (5.10) agree with the formulas (3.74) and (3.75) in [19], respectively.

5.3 Velocity-dependent second-order perturbations

As is easily seen from Eqs. (4.11) and (4.12), the velocity-dependent second-order potentials in the post-Newtonian approximation read, respectively,

[TABLE]

and

[TABLE]

In the same approximation, the metric corrections $h_{00}$ and $h_{0\alpha}$ are presented, e.g., in the textbook [25] (see Eqs. (106.13) and (106.15), respectively). Obviously, the sum of our expressions (5.10) and (5.12) should be in agreement with (106.13). The direct comparison shows that one-to-one coincidence is absent. Indeed, the first two terms in the right-hand side of (5.10) and the first term in the right-hand side of (5.12) coincide (up to evident redefinitions of spatial coordinates and peculiar velocities: $a{\bf r}\mapsto{\bf r}$, $a{\bf r}_{n}\mapsto{\bf r}_{n}$ and ${\bf\tilde{v}}_{n}\mapsto{\bf v}_{n}/c$) with the corresponding terms in (106.13), divided by $2$, as it should be. However, each of Eqs. (5.10) and (5.12) contain additional terms which are absent in (106.13). As was emphasized in [14] and [19], the reason for this noncoincidence with the Landau & Lifshitz formulas is the choice of different gauges here and in [25]. To prove that the results are in accordance, let us consider an appropriate coordinate transformation connecting our formulas and (106.13). Such a transformation was proposed in [19]: $\eta\mapsto\eta-A(\eta,\mathbf{r})$. Then $\Phi^{(2)}\mapsto\Phi^{(2)}+A^{\prime}$. In this case the first-order vector perturbation is also transformed: $\mathbf{B}\mapsto\mathbf{B}+\nabla A$. Hence, the problem is to determine the function $A$ which provides the transition to the Landau & Lifshitz formulas.

In the small-scale limit the first-order vector perturbation is presented in [14] and [19]. Its comparison with Eq. (106.15) in [25] gives (see Eq. (3.80) in [19]):

[TABLE]

The time derivative of this expression reads

[TABLE]

The first line of this formula is the time derivative of (5.13) with fixed positions $\mathbf{r}_{n}$ (i.e. when the time derivatives of $\mathbf{r}_{n}$ are not taken into account). This line was determined previously in [19]. The second and third lines represent the time derivative of (5.13) with fixed velocities $\mathbf{\tilde{v}}_{n}$. These two lines are novel.

On the other hand, the direct comparison of the formula (106.13) with Eqs. (5.10) and (5.12) shows that

[TABLE]

It is not difficult to verify that the gradient of (5.15) exactly coincides with (5.14). Therefore, we have proved that the sought-for function $A$ exists and the agreement with the Landau & Lifshitz formulas is achieved.

6 Conclusion

In the present paper we have studied the second-order scalar perturbations for the $\Lambda$CDM cosmological model within the cosmic screening approach. We have found the analytic expressions for both the velocity-independent and velocity-dependent perturbations (see Eqs. (3.10), (3.25), (4.11), (4.12)). We have demonstrated that the Yukawa screening effect, being inherent in the first-order metric corrections (see, e.g., [14, 16, 18]), affects the second-order ones as well. In addition, it has been shown that the obtained expressions for the second-order perturbations are reduced to the known post-Newtonian formulas [19, 25] at distances much smaller than the screening length. With respect to the velocity-independent perturbations, such correspondence was already observed in [19]. However, in the case of the velocity-dependent perturbations this is a novel result.

In the era of precision cosmology, the derived analytic formulas play an important role since they enable revealing the contributions of the second-order metric corrections to the observed physical quantities and effects (including the luminosity-redshift relation, gravitational lensing, baryon acoustic oscillations). The analytic expressions for the second-order perturbations give an opportunity to estimate the backreaction and determine how strong the backreaction is and to what extent we may use the standard FLRW metric as a background one. Obviously, the perturbative approach is robust if the second-order corrections are much smaller than the first-order ones. If this is the case, it is usually enough to be limited to the first order. However, how can we know it from numerical simulations? Performing a numerical simulation, which takes into account the second-order perturbations, is not an easy task. Instead of such a complicated procedure, we suggest the following test. The simulation can be performed on the basis of the first-order approach. After that one can calculate the first-order perturbation $\Phi$ at a number of points. On the other hand, using our analytic expressions, one can calculate the second-order perturbations $\Phi^{(2)}$ and $\Psi^{(2)}$ at the same points. Then, if $\Phi\gg\Phi^{(2)},\Psi^{(2)}$, the perturbative scheme is robust, and the backreaction is apparently negligible. Otherwise the backreaction should be certainly taken into account.

Acknowledgements

We would like to thank the Referee for the valuable comments which have considerably improved the presentation and discussion of the obtained results.

Appendix A Integrals

In this appendix we present the integrals used for our calculations. First, we list three integrals where integration over the radial coordinate runs from zero to infinity:

[TABLE]

Now we consider the integrals which correspond to the model described in section 5. In this case integration over the radial coordinate runs from the radius $R$ of the sphere to infinity:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]