1 Introduction

Modern observational phenomena, such as dark energy and dark matter, are the great challenge for present cosmology, astrophysics and theoretical physics. Within the scope of standard models, a satisfactory explanation to these problems has not been offered yet. This forces the search of their solutions beyond the standard models, for example, by considering modified gravitational theories. One of the possible generalizations consists in consideration of nonlinear (with respect to the scalar curvature \(R\)) models \(f(R)\). Nonlinear models may arise either due to quantum fluctuations of matter fields including gravity [1], or as a result of compactification of extra spatial dimensions [2]. Starting from the pioneering paper [3], the nonlinear theories of gravity \(f(R)\) have attracted a great deal of interest because these models can provide a natural mechanism of the early inflation. It was also realized that these models can explain the late-time acceleration of the Universe. This fact resulted in a new wave of papers devoted to this topic (see, e.g., the reviews [4–9]).

The cosmological perturbation theory is very important for the current cosmological investigations of the large-scale structure. Thus, it would be interesting to make the corresponding cosmological perturbation analysis in nonlinear \(f(R)\) theories of gravity. In the hydrodynamical approach, such investigation was performed in a number of papers (see, e.g., Sec. 8 in the review [6] and references therein). In particular, matter density perturbations in a class of viable cosmological \(f(R)\) models were studied in [10, 11]. We consider the Universe at the late stage of its evolution when galaxies and clusters of galaxies have already formed. At scales much larger than the characteristic distance between these inhomogeneities, the Universe is well described by the homogeneous and isotropic FRW metric. These scales are approximately 190 Mpc and larger [12]. At these distances, the matter fields (e.g., cold dark matter) are well described by the hydrodynamical approach. However, at smaller scales the Universe is highly inhomogeneous, and we need to take into account the inhomogeneities in the form of galaxies, groups and clusters of galaxies. The peculiar velocities of these inhomogeneities are much less than the speed of light, and we can use the nonrelativistic approximation. This means that in equations for scalar perturbations we first neglect peculiar velocities and solve these equations with respect to scalar perturbation functions \(\Phi \) and \(\Psi \). The function \(\Phi \) represents the gravitational potential of the inhomogeneities. Then we use the explicit expression for \(\Phi \) to describe the motion of inhomogeneities. Such mechanical approach is well known in astrophysics (see e.g., [13]). We generalized it to the case of dynamical cosmological background [12, 14]. In the case of the linear model (i.e., the conventional \(\Lambda \)CDM model), we used this procedure to describe the mutual motion of galaxies and dwarf galaxies [14, 15]. Due to the great popularity of the nonlinear \(f(R)\) models, it is of interest to apply this scheme to them. However, first of all, we should show that nonlinear theories are compatible with the mechanical approach. In other words, we have to examine equations for scalar perturbations of the metrics in nonlinear \(f(R)\) gravity within the framework of the mechanical approach to show their integrability up to the required accuracy. This is the main aim of our paper.

As a result, we demonstrate that considered in our paper nonlinear theories are compatible with the mechanical approach. We get the expressions for both \(\Phi \) and \(\Psi \) in different approximations. Moreover, the exact form of the gravitational potential \(\Phi \) gives a possibility to take into account both the effects of nonlinearity of the original model and the dynamics of the cosmological background. The explicit form of this function makes it possible to carry out analytical and numerical study of mutual motion of galaxies in nonlinear models. Therefore, our formulas can be used to analyze the large-scale structure dynamics in the late Universe for nonlinear \(f(R)\) models. This is the main result of our paper.

The paper is structured as follows. In Sect. 2 we present the main background equations as well as the equations for scalar perturbations for an arbitrary \(f(R)\) model. Here the equations for scalar perturbations are written within the framework of the mechanical approach. In Sect. 3 we solve these equations in three approximations: the astrophysical approach, the large scalaron mass case and the quasi-static approximation. In all three cases we obtain the expressions for the scalar perturbation functions \(\Phi \) and \(\Psi \) up to the required accuracy. The main results are summarized in concluding Sect. 4.

2 Basic equations

In this section we reproduce some known equations of the nonlinear \(f(R)\) gravitational model that we will use hereinafter. We follow mainly the review [6] using the notation and the sign convention accepted in this paper. In \(f(R)\) gravity, the action reads

$$\begin{aligned} S=\frac{1}{2\kappa ^2}\int d^4x\sqrt{-g}f(R) + S_m, \end{aligned}$$
(2.1)

where \(f(R)\) is an arbitrary smooth function of the scalar curvature \(R\), \(S_m\) is the action of matter, \(\kappa ^2=8\pi G_N\), and \(G_N\) is the Newtonian gravitational constant. The equation of motion corresponding to this action is (see e.g., Eq. 2.4 in [6])

$$\begin{aligned}&F(R)R_{\mu \nu } - \frac{1}{2}f(R)g_{\mu \nu } - \nabla _{\mu }\nabla _{\nu }F(R) + g_{\mu \nu }\square F(R)\nonumber \\&\quad ={\kappa }^2 T_{\mu \nu }\, \quad \mu ,\nu =0,1,2,3. \end{aligned}$$
(2.2)

The trace of this equation gives

$$\begin{aligned} 3\square F(R)+F(R)R-2f(R)=\kappa ^2 T. \end{aligned}$$
(2.3)

Here, \(F(R)=f'(R)\) and \(T=g^{\mu \nu }T_{\mu \nu }\). Besides, \(\square F = (1/\sqrt{-g})\partial _{\mu }(\sqrt{-g}g^{\mu \nu }\partial _{\nu }F)\). In what follows, the prime denotes the derivative with respect to the scalar curvature \(R\). In our paper we consider a special class of \(f(R)\) models which have solutions \(R_{\mathrm {dS}}\) of the equation

$$\begin{aligned} F(R)R-2f(R)=0. \end{aligned}$$
(2.4)

As follows from Eq. (2.3), they are vacuum solutions (\(T=0\)) of this equation for which the Ricci scalar is constant (\(\square F(R)=0\)). Such solutions are called de Sitter points [6, 16, 17]. According to [18, 19], viable nonlinear models should have stable de Sitter points in the late Universe. We can expand the function \(f(R)\) in the vicinity of one of these points:

$$\begin{aligned} f(R)&= f(R_{\mathrm {dS}})+f'(R_{\mathrm {dS}})(R-R_{\mathrm {dS}})+o(R-R_{\mathrm {dS}})\nonumber \\&= -f(R_{\mathrm {dS}})+ \frac{2f(R_{\mathrm {dS}})}{R_{\mathrm {dS}}}R + o(R-R_{\mathrm {dS}}), \end{aligned}$$
(2.5)

where we used Eq. (2.4). Now, in order to have linear gravity at the late stage of the Universe evolution [6], without loss of generality we choose the parameters of the model in such a way that

$$\begin{aligned} 2\frac{f(R_{\mathrm {dS}})}{R_{\mathrm {dS}}}=1 \quad \Rightarrow \quad f(R_{\mathrm {dS}})=\frac{R_{\mathrm {dS}}}{2}. \end{aligned}$$
(2.6)

Therefore, we get

$$\begin{aligned} f(R)=-2\Lambda +R +o(R-R_{\mathrm {dS}}), \end{aligned}$$
(2.7)

where \(\Lambda \equiv R_{\mathrm {dS}}/4\). The stability of these points was discussed in [6, 20]. Obviously, these models go asymptotically to the de Sitter space when \(R\rightarrow R_{\mathrm {dS}}\ne 0\) with a cosmological constant \(\Lambda =R_{\mathrm {dS}}/4\). This happens when the matter content becomes negligible with respect to \(\Lambda \) as in the late Friedman–Robertson–Walker (FRW) cosmology. We can also consider a zero solution \(R_{\mathrm {dS}}=0\) of Eq. (2.4). It is more correct to call this point a Minkowski one. Here, \(\Lambda =0\), and such models go asymptotically to the Minkowski space. In particular, three popular models, Starobinsky [21], Hu–Sawicki [22] and MJWQ [23], have stable nonzero de Sitter points in future (approximately at the redshift \(z=-1\)) [24, 25]. The explicit search for dS points in both future and past was considered in [17, 26]. It is worth noting that in papers [18, 27] the authors point to the oscillating behavior of the parameter of the equation of state near the value \(-1\) in the future. Moreover, the number of times of such oscillations can be infinite [19].

In the case of the spatially flat background spacetime with the FRW metric

$$\begin{aligned} ds^2=g_{\mu \nu }dx^{\mu }dx^{\nu }=-dt^2+a^2(t) \left( dx^2+dy^2+dz^2\right) \end{aligned}$$
(2.8)

and matter in the form of a perfect fluid with the energy-momentum tensor components \(\bar{T}^{\mu }_{\nu }=\mathrm {diag}(-\bar{\rho },\bar{P},\bar{P},\bar{P})\), Eq. (2.2) results in the following system:

$$\begin{aligned} 3FH^2 = (FR - f)/2 - 3H\dot{F} + \kappa ^2 \bar{\rho }, \end{aligned}$$
(2.9)

and

$$\begin{aligned} -2F\dot{H}= \ddot{F} - H\dot{F} + \kappa ^2 (\bar{\rho }+ \bar{P}), \end{aligned}$$
(2.10)

where the bar denotes the homogeneous background quantities, the Hubble parameter \(H=\dot{a}/a\) (the dot everywhere denotes the derivative with respect to the synchronous time \(t\)) and the scalar curvature

$$\begin{aligned} R=6\left( 2H^2 + \dot{H}\right) . \end{aligned}$$
(2.11)

The perfect fluid satisfies the continuity equation

$$\begin{aligned} \dot{\bar{\rho }} + 3H(\bar{\rho }+ \bar{P})=0, \end{aligned}$$
(2.12)

which for nonrelativistic matter with \(P=0\) has the solution

$$\begin{aligned} \bar{\rho }= \bar{\rho }_c/a^3, \end{aligned}$$
(2.13)

where \(\bar{\rho }_c=\text{ const }\) is the rest mass density in the comoving coordinates.

Above, Eqs. (2.8)–(2.13) describe the homogeneous background. As we have written in the Sect. 1, we consider the Universe at late stages of its evolution when galaxies and clusters of galaxies have already formed and when the Universe is highly inhomogeneous inside the cell of uniformity which is approximately 190 Mpc in size [12]. Obviously, these inhomogeneities perturb the homogeneous background. At scales larger than the cell of uniformity size, the matter fields (e.g., cold dark matter) are well described by the hydrodynamical approach. However, at smaller scales the mechanical approach looks more adequate [12, 14]. In the framework of the mechanical approach galaxies, dwarf galaxies and clusters of galaxies (composed of baryonic and dark matter) can be considered as separate compact objects. Moreover, at distances much greater than their characteristic sizes they can be described well as point-like matter sources with the rest mass density

$$\begin{aligned} \rho = \frac{1}{a^3}\sum _i m_i\delta {(\mathbf {r}-\mathbf {r}}_i)\equiv \rho _c/a^3, \end{aligned}$$
(2.14)

here \(\mathbf{{r}}_i\) is the radius-vector of the \(i\)th gravitating mass in the comoving coordinates. This is the generalization of the well-known astrophysical approach [13] to the case of dynamical cosmological background. Usually, the gravitational fields of these inhomogeneities are weak and their peculiar velocities are much less than the speed of light. All these inhomogeneities/fluctuations result in scalar perturbations of the FRW metric (2.8). In the conformal Newtonian (longitudinal) gauge, such perturbed metric is [6, 28, 29]

$$\begin{aligned} ds^2= -(1+2\Phi )dt^2+a^2(1-2\Psi )\left( dx^2+dy^2+dz^2\right) , \end{aligned}$$
(2.15)

where the introduced scalar perturbations \(\Phi ,\Psi \ll 1\). These functions of all spacetime coordinates, representing deviations of metric coefficients from their average/background values, may be associated with famous Bardeen’s potentials [28] under the made gauge choice. It is worth noting that smallness of these nonrelativistic gravitational potentials \(\Phi \) and \(\Psi \) and smallness of peculiar velocities are two independent conditions (e.g., for very light relativistic masses the gravitational potential can still remain small). Therefore, similar to the astrophysical approach described in [13] (see §106), we split the investigation of galaxy dynamics in the late Universe into two steps. First, we neglect peculiar velocities and define the gravitational potential \(\Phi \). Then we use this potential to determine dynamical behavior of galaxies. It leads to a possibility to take into account both the gravitational attraction between inhomogeneities and the global cosmological expansion of the Universe. For example, for the linear model \(f(R)=R\) this procedure was performed in [15]. Our present paper is devoted to the first step in this program. In other words, we are going to define scalar perturbations \(\Phi ,\Psi \) for the \(f(R)\) gravitational models. Under our assumptions and according to [6, 30, 31], these perturbations satisfy the following system of equations:

$$\begin{aligned}&-\frac{\Delta \Psi }{a^2} + 3H\left( H\Phi +\dot{\Psi }\right) \nonumber \\&\quad =-\frac{1}{2F}\left[ \left( 3H^2+ 3 \dot{H} +\frac{\Delta }{a^2}\right) \delta F - 3H\dot{\delta F}\right. \nonumber \\&\qquad + \left. 3H\dot{F}\Phi +3\dot{F} \left( H\Phi +\dot{\Psi }\right) +\kappa ^2 \delta \rho \right] , \end{aligned}$$
(2.16)
$$\begin{aligned}&H\Phi +\dot{\Psi } = \frac{1}{2F}\left( \dot{\delta F} - H \delta {F} - \dot{F}\Phi \right) , \end{aligned}$$
(2.17)
$$\begin{aligned}&-F(\Phi -\Psi )=\delta F, \end{aligned}$$
(2.18)
$$\begin{aligned}&3\left( \dot{H}\Phi +H\dot{\Phi }+\ddot{\Psi }\right) + 6H\left( H\Phi +\dot{\Psi }\right) + 3\dot{H}\Phi +\frac{\Delta \Phi }{a^2} \nonumber \\&\quad =\frac{1}{2F}\left[ 3\ddot{\delta F}+3H\dot{\delta F}-6H^2\delta F - \frac{\Delta \delta F}{a^2} - 3\dot{F}\dot{\Phi }\right. \nonumber \\&\qquad - \left. 3\dot{F}\left( H\Phi +\dot{\Psi }\right) - \left( 3H\dot{F} + 6\ddot{F}\right) \Phi + \kappa ^2\delta \rho \right] , \end{aligned}$$
(2.19)
$$\begin{aligned}&\ddot{\delta F}+ 3H \dot{\delta F} - \frac{\Delta \delta F}{a^2}-\frac{1}{3}R\delta F\nonumber \\&\quad =\frac{1}{3}\kappa ^2(\delta \rho - 3\delta P) + \dot{F}(3H\Phi +3\dot{\Psi }+\dot{\Phi })\nonumber \\&\qquad +2\ddot{F}\Phi + 3H\dot{F}\Phi - \frac{1}{3}F\delta {R}, \end{aligned}$$
(2.20)
$$\begin{aligned}&\delta F=F'\delta R,\nonumber \\&\delta R=-2\bigg [3\left( \dot{H}\Phi +H\dot{\Phi }+\ddot{\Psi }\right) + 12H\left( H\Phi +\dot{\Psi }\right) \nonumber \\&\qquad \qquad +\frac{\Delta \Phi }{a^2}+3\dot{H}\Phi -2\frac{\Delta \Psi }{a^2}\bigg ]. \end{aligned}$$
(2.21)

In these equations the function \(F\), its derivative \(F'\) and the scalar curvature \(R\) are unperturbed background quantities. Here \(\Delta \) is the Laplacian in the comoving coordinates. As a matter source, we consider dust-like matter. Therefore, \(\delta P =0\) and

$$\begin{aligned} \delta \rho = \rho -\bar{\rho }=(\rho _c-\bar{\rho }_c)/a^3, \end{aligned}$$
(2.22)

where \(\bar{\rho }\) and \(\rho \) are defined in Eqs. (2.13) and (2.14), respectively.

It can easily be verified that in the linear case \(f(R)=R\ \Rightarrow \ F(R)=1\) this system of equations is reduced to Eqs. (2.18)–(2.20) in [14].

3 Astrophysical and cosmological approaches

3.1 Astrophysical approach

First, we consider Eqs. (2.16)–(2.21) in the astrophysical approach. This means that we neglect the time dependence of functions in these equations by setting all time derivatives equal to zero. It is supposed also that the background model is matter free, i.e., \(\bar{\rho }=0\). As we mentioned above, there are two types of vacuum background solutions of Eq. (2.3): de Sitter spacetime with \(R_{\mathrm {dS}}=12H^2=\text{ const }\ne 0\) and Minkowski spacetime with \(R=0,H=0\). However, the system of Eqs. (2.16)–(2.21) was obtained for the FRW metric (2.15) where we explicitly took into account the dependence of the scale factor \(a\) on time. Therefore, if we want to get the time independent astrophysical equations directly from Eqs. (2.16)–(2.21), we should also neglect the time dependence of \(a\), i.e., the background Hubble parameter \(H=0\). This means that the background solution is the Minkowski spacetime. This background is perturbed by dust-like matter with the rest mass density Eq. (2.14). Keeping in mind that \(\bar{\rho }=0\), we have \(\delta \rho = \rho \).

In the case of Minkowski spacetime background, dropping the time derivatives, Eqs. (2.16)–(2.21) in the astrophysical approach are reduced to the following system:

$$\begin{aligned} -\frac{\Delta }{a^2}\Psi =-\frac{1}{2F}\left( \frac{\Delta }{a^2}\delta F + \kappa ^2 \delta \rho \right) , \end{aligned}$$
(3.1)
$$\begin{aligned} -F(\Phi - \Psi )=\delta F, \end{aligned}$$
(3.2)
$$\begin{aligned} \frac{\Delta }{a^2}\Phi =\frac{1}{2F}\left( -\frac{\Delta }{a^2}\delta F + \kappa ^2 \delta \rho \right) , \end{aligned}$$
(3.3)
$$\begin{aligned} -\frac{\Delta }{a^2}\delta F=\frac{1}{3} \kappa ^2\delta \rho -\frac{1}{3}F\delta R, \end{aligned}$$
(3.4)
$$\begin{aligned} \delta F=F'\delta R,\quad \delta R=-2\left( \frac{\Delta }{a^2}\Phi - 2\frac{\Delta }{a^2}\Psi \right) . \end{aligned}$$
(3.5)

From Eqs. (3.1) and (3.3) we obtain, respectively,

$$\begin{aligned} \begin{aligned}&\Psi =\frac{1}{2F}\delta F +\frac{\varphi }{a}=\frac{F'}{2F}\delta R +\frac{\varphi }{a},\\&\Phi =-\frac{1}{2F}\delta F+\frac{\varphi }{a}=-\frac{F'}{2F}\delta R+\frac{\varphi }{a}, \end{aligned} \end{aligned}$$
(3.6)

where the function \(\varphi \) satisfies the equation

$$\begin{aligned} \Delta \varphi =\frac{1}{2F}\kappa ^2a^3\delta \rho =\frac{1}{2F}\kappa ^2\rho _c=\frac{4\pi G_N}{F}\rho _c. \end{aligned}$$
(3.7)

Here, we took into consideration that in the astrophysical approach \(\delta \rho _c=\rho _c\) where \(\rho _c\) is defined by Eq. (2.14). It is worth noting that in the Poisson equation (3.7) the Newtonian gravitational constant \(G_N\) is replaced by an effective one \(G_{\mathrm {eff}}=G_N/F\).

Equation (3.2) follows directly from Eq. (3.6) and therefore may be dropped, while from Eq. (3.4) we get the following Helmholtz equation with respect to the scalaron function \(\delta R\):

$$\begin{aligned} \Delta \delta R - \frac{a^2 F}{3 F'}\delta R = -\frac{a^2F}{3F'}\frac{\kappa ^2}{Fa^3} \delta \rho _c. \end{aligned}$$
(3.8)

On the other hand, it can easily be seen that the substitution of Eqs. (3.6) and (3.7) into Eq. (3.5) results in the same Eq. (3.8). Therefore, in the case of Minkowski background, the mass of the scalaron squared is

$$\begin{aligned} M^2 = \frac{a^2 }{3}\frac{F}{F'}. \end{aligned}$$
(3.9)

A similar formula (up to the evident substitution \(a=1\) and the zero background scalar curvature) for the mass squared can be found, e.g., in [32–35] (see also Eq. (5.2) in the review [6]).

3.2 Cosmological approach

Now we want to take into consideration cosmological evolution. This means that background functions may depend on time. In this case it is hardly possible to solve the system (2.16)–(2.21) directly. Therefore, first, we study the case of the very large mass of the scalaron. It should be noted also that we investigate the Universe filled with nonrelativistic matter with the rest mass density \(\bar{\rho }\sim 1/a^3\). Hence, we will drop all terms which decrease (with increasing \(a\)) faster than \(1/a^3\). This is the accuracy of our approach. Within this approach, \(\delta \rho \sim 1/a^3\) [14].

3.2.1 Large scalaron mass case

As we can see from Eq. (3.9), the limit of the large scalaron mass corresponds to \(F' \rightarrow 0\). Then \(\delta F\) is also negligible [see Eq. (2.21)]. Therefore, Eqs. (2.16)–(2.21) read

$$\begin{aligned}&-\frac{\Delta \Psi }{a^2} + 3H\left( H\Phi +\dot{\Psi }\right) \nonumber \\&\quad = -\frac{1}{2F}\left[ 3H\dot{F}\Phi +3\dot{F} \left( H\Phi +\dot{\Psi }\right) {+\kappa ^2\delta \rho }\right] , \end{aligned}$$
(3.10)
$$\begin{aligned} H\Phi +\dot{\Psi } = \frac{1}{2F}\left( - \dot{F}\Phi \right) , \end{aligned}$$
(3.11)
$$\begin{aligned} \Phi -\Psi =0, \end{aligned}$$
(3.12)
$$\begin{aligned}&3\left( \dot{H}\Phi +H\dot{\Phi }+\ddot{\Psi }\right) + 6H\left( H\Phi +\dot{\Psi }\right) + 3\dot{H}\Phi +\frac{\Delta \Phi }{a^2} \nonumber \\&\quad =\frac{1}{2F}\Big [- 3\dot{F}\dot{\Phi }-3\dot{F}\left( H\Phi +\dot{\Psi }\right) \nonumber \\&\qquad \qquad - \left( 3H\dot{F} + 6\ddot{F}\right) \Phi {+\kappa ^2\delta \rho }\Big ], \end{aligned}$$
(3.13)
$$\begin{aligned} 0&= \dot{F}(3H\Phi +3\dot{\Psi }+\dot{\Phi })+2\ddot{F}\Phi + 3H\dot{F}\Phi \nonumber \\&+\,\frac{1}{3}\kappa ^2\delta \rho -\frac{1}{3}F\delta R, \end{aligned}$$
(3.14)
$$\begin{aligned} -\frac{1}{2}\delta R&= 3\left( \dot{H}\Phi +H\dot{\Phi }+\ddot{\Psi }\right) + 12H\left( H\Phi +\dot{\Psi }\right) \nonumber \\&+\frac{\Delta \Phi }{a^2}+3\dot{H}\Phi - 2\frac{\Delta \Psi }{a^2}. \end{aligned}$$
(3.15)

From Eqs. (3.11) and (3.12) we get

$$\begin{aligned} \Psi =\Phi =\frac{\varphi }{a\sqrt{F}}, \end{aligned}$$
(3.16)

where the introduced function \(\varphi \) depends only on spatial coordinates. Substituting Eq. (3.16) into Eq. (3.10), we obtain

$$\begin{aligned} \frac{1}{a^3\sqrt{F}}\Delta \varphi + \frac{{3\dot{F}}^2}{4aF^2\sqrt{F}}\varphi ={\frac{1}{2F}\kappa ^2\delta \rho }. \end{aligned}$$
(3.17)

As we mentioned above, neglecting relativistic matter in the late Universe, we have \(\delta \rho \sim 1/a^3\) [14]. This approximation is getting better and better performed in the limit \(a\rightarrow +\infty \). We assume that this limit corresponds to the final stage of the Universe evolution. The similar limit with respect to the scalar curvature is \(R\rightarrow R_{\infty }\), where the value \(R_{\infty }\) is just finite. Then from Eq. (3.17) we immediately come to the condition

$$\begin{aligned} F=\mathrm {const}+o(1), \end{aligned}$$
(3.18)

where \(o(1)\) is any decreasing (with increasing \(a\)) function of \(a\). This condition holds at the considered late stage. One can also naively suppose that in the late Universe \(\dot{F}\sim 1/a+o(1/a)\). However, as we will see below, this is wrong. Obviously, without loss of generality we can suppose that \(\mathrm {const}=1\). From the condition in Eq. (3.18) we get

$$\begin{aligned} F=1+o(1) \quad \Rightarrow \quad f=-2\Lambda +R+o(R-R_{\infty }), \end{aligned}$$
(3.19)

where \(\Lambda \) is the cosmological constant. Therefore, the limit of the large scalaron mass takes place for models which possess the asymptotic form (3.19). For example, \(R_{\infty }\) may correspond to the de Sitter point \(R_{\mathrm {dS}}\) in future [see Eq. (2.7)]. As we have written in Sect. 2, all three popular models, Starobinsky [21], Hu–Sawicki [22] and MJWQ [23], have such stable de Sitter points in future (approximately at the redshift \(z=-1\)) [24, 25]. The condition of stability is \(0<R F'/F< 1\) (see, e.g., (4.80) in [6]). Since \(F\approx 1\), this condition reads \(0<R< 1/F'\) which is fulfilled for the de Sitter points in the above-mentioned models. The reason of it consists in the smallness of \(F'\).

We now return to the remaining Eqs. (3.13)–(3.15) to show that they are satisfied within the considered accuracy. First, we study Eq. (3.13) which after the substitution of Eqs. (3.16) and (3.17) and some simple algebra takes the form

$$\begin{aligned} \frac{\varphi }{a}\dot{H} -\frac{\varphi }{2aF}\left( H\dot{F}- \ddot{F}\right) =0. \end{aligned}$$
(3.20)

To estimate \(\dot{F}\) and \(\ddot{F}\), we take into account that in the limit \(R\rightarrow R_{\infty }\), \(F\approx 1\), \(H\approx \mathrm {const}\ \Rightarrow \ \dot{H} \approx 0\), and \(F'(R_{\infty })\) is some finite positive value. Then \(\dot{F}=F'\dot{R}\sim F'(R_{\infty })\dot{R}\sim \dot{T} \sim d\left( 1/a^3\right) /dt\sim H\left( 1/a^3\right) \sim 1/a^3\) and \(\ddot{F} \sim \dot{a}/a^4\sim 1/a^3\). Therefore, the left hand side of Eq. (3.20) is of the order \(o(1/a^3)\) and we can put it zero within the accuracy of our approach. Similarly, Eqs. (3.14) and (3.15) are satisfied within the considered accuracy. It can also be seen that the second term on the left hand side of Eq. (3.17) is of the order \(O(1/a^7)\) and should be eliminated.

Thus, in the case of the large enough scalaron mass we reproduce the linear cosmology from the nonlinear one, as it should be.

3.2.2 Quasi-static approximation

Now we do not want to assume a priori that the scalaron mass is large, i.e., \(F'\) can have arbitrary values. Hence, we will preserve the \(\delta F\) terms in Eqs. (2.16)–(2.21). Moreover, we should keep the time derivatives in these equations. Such a system is very complicated for the direct integration. However, we can investigate it in the quasi-static approximation [36, 37] (see also Sec. 8.1 in [6]). According to this approximation, the spatial derivatives give the main contribution to Eqs. (2.16)–(2.21). Therefore, first, we should solve “astrophysical” Eqs. (3.1)–(3.5), and then check whether the found solutions satisfy (up to the adopted accuracy) the full system of equations. In other words, in the quasi-static approximation it is naturally supposed that the gravitational potentials (the functions \(\Phi ,\Psi \)) are produced mainly by the spatial distribution of astrophysical/cosmological bodies. As we have seen, Eqs. (3.1)–(3.5) result in Eqs. (3.6)–(3.9). Now, we should keep in mind that we have the cosmological background. Moreover, we consider the late Universe which is not far from the de Sitter point \(R_{\mathrm {dS}}\) in future.Footnote 1 This means that \(\delta \rho = \rho - \bar{\rho }\) in Eq. (3.7), all background quantities (e.g., \(F,\,F'\)) are calculated roughly speaking at \(R_{\mathrm {dS}}\) and the scalaron mass squared (3.9) reads now [6, 26, 32–34]

$$\begin{aligned} M^2 = \frac{a^2 }{3}\left( \frac{F}{F'}-R_{\mathrm {dS}}\right) . \end{aligned}$$
(3.21)

Let us consider now Eq. (3.8) with the mass squared (3.21). Taking into account that now \(\delta \rho _c=\rho _c-\bar{\rho }_c\), we can rewrite this equation as follows:

$$\begin{aligned} \Delta \widetilde{\delta R} - M^2\widetilde{\delta R} + \frac{a^2}{3F'}\frac{\kappa ^2}{a^3} \sum _i m_i\delta {(\mathbf {r}-\mathbf {r}}_i)=0, \end{aligned}$$
(3.22)

where

$$\begin{aligned} \widetilde{\delta R}= \delta R + \frac{a^2}{3F'M^2}\frac{\kappa ^2}{a^3}\bar{\rho }_c=\delta R + \frac{\kappa ^2}{(F-F'R_{\mathrm {dS}})a^3}\bar{\rho }_c. \end{aligned}$$
(3.23)

Equation (3.22) demonstrates that we can apply the principle of superposition solving this Helmholtz equation for one gravitating mass \(m_i\). Then the general solution for a full system is the sum over all gravitating masses. As boundary conditions, we require for each gravitating mass the behavior \(\widetilde{\delta R}\sim 1/r\) at small distances \(r\) and \(\widetilde{\delta R}\rightarrow 0\) for \(r\rightarrow \infty \). Taking all these remarks into consideration, we obtain for a full system

$$\begin{aligned} \delta R&= \frac{\kappa ^2}{12\pi aF'}\sum _i\frac{m_i\exp \left( -M {|\mathbf {r}-\mathbf {r}}_i|\right) }{{|\mathbf {r}-\mathbf {r}}_i|}\nonumber \\&-\frac{\kappa ^2}{(F-F'R_{\mathrm {dS}})a^3}\bar{\rho }_c \, . \end{aligned}$$
(3.24)

It is worth noting that averaging over the whole comoving spatial volume \(V\) gives the zero value \(\overline{\delta R}=0\). Really, since \(\sum _im_i/V=\bar{\rho }_c\),

$$\begin{aligned} \overline{\delta R}&= \frac{1}{V}\int \limits _V\delta R\mathrm{d}V=\frac{1}{V}\frac{\kappa ^2}{12\pi aF'}\nonumber \\&\quad \times \sum _im_i\frac{4\pi }{M^2}-\frac{\kappa ^2}{(F-F'R_{\mathrm {dS}})a^3}\bar{\rho }_c=0. \end{aligned}$$
(3.25)

This result is reasonable because the rest mass density fluctuation \(\delta \rho \), representing the source of the metric and scalar curvature fluctuations \(\Phi ,\Psi \) and \(\delta R\), has a zero average value \(\overline{\delta \rho }=0\). Consequently, all enumerated quantities should also have zero average values: \(\bar{\Phi }=\bar{\Psi }=0\) and \(\overline{\delta R}=0\), in agreement with Eq. (3.25).

From Eq. (3.6) we get the scalar perturbation functions \(\Psi \) and \(\Phi \) in the following form:

$$\begin{aligned} \Psi&= \frac{F'}{2F}\left[ \frac{\kappa ^2}{12\pi aF'}\sum _i\frac{m_i\exp \left( -M {|\mathbf {r}-\mathbf {r}}_i|\right) }{{|\mathbf {r}-\mathbf {r}}_i|}\right. \nonumber \\&-\left. \frac{\kappa ^2}{(F-F'R_{\mathrm {dS}})a^3}\bar{\rho }_c\right] +\frac{\varphi }{a}, \end{aligned}$$
(3.26)
$$\begin{aligned} \Phi&= -\frac{F'}{2F}\left[ \frac{\kappa ^2}{12\pi aF'}\sum _i\frac{m_i\exp \left( -M {|\mathbf {r}-\mathbf {r}}_i|\right) }{{|\mathbf {r}-\mathbf {r}}_i|}\right. \nonumber \\&-\left. \frac{\kappa ^2}{(F-F'R_{\mathrm {dS}})a^3}\bar{\rho }_c\right] +\frac{\varphi }{a}, \end{aligned}$$
(3.27)

where \(\varphi \) satisfies Eq. (3.7) with \(\delta \rho \) in the form (2.22) (i.e., \(\bar{\rho }_c\ne 0\)). Obviously, when \(F'\rightarrow 0\), \(M\rightarrow \infty \), and we have \(\exp \left( -M {|\mathbf {r}-\mathbf {r}}_i|\right) /{|\mathbf {r}-\mathbf {r}}_i|\rightarrow 4\pi \delta (\mathbf{r}-\mathbf{r}_i)/M^2\), so the expression in the square brackets in Eqs. (3.26) and (3.27) is equal to \(\kappa ^{2} \delta \rho _{c}/\left[ (F-F'R_{\mathrm {dS}})a^3\right] \). Therefore, in the considered limit \(F'\rightarrow 0\) we reproduce the scalar perturbations \(\Phi ,\Psi \) from the previous large scalaron mass case, as it certainly should be.

Thus, neglecting for a moment the influence of the cosmological background, but without neglecting the scalaron’s contribution, we have found the scalar perturbations. They represent the mix of the standard potential \(\varphi /a\) (see the linear case [14]) and the additional Yukawa term which follows from the nonlinearity.

Now we should check that these solutions satisfy the full system (2.16)–(2.21). To do it, we substitute Eqs. (3.24), (3.26), and (3.27) into this system of equations. Obviously, the spatial derivatives disappear. Keeping in mind this fact, the system (2.16)–(2.21) is reduced to the following equations:

$$\begin{aligned} 3H\left( H\Phi +\dot{\Psi }\right)&= -\frac{1}{2F}\bigg [\left( 3H^2+ 3 \dot{H}\right) \delta F \nonumber \\&\!-\! 3H\dot{\delta F}\!+\!3H\dot{F}\Phi \!+\!3\dot{F} \left( H\Phi \!+\!\dot{\Psi }\right) \bigg ],\nonumber \\ \end{aligned}$$
(3.28)
$$\begin{aligned} H\Phi +\dot{\Psi } = \frac{1}{2F}\left( \dot{\delta F} - H \delta {F} - \dot{F}\Phi \right) , \end{aligned}$$
(3.29)
$$\begin{aligned}&3\left( \dot{H}\Phi +H\dot{\Phi }+\ddot{\Psi }\right) + 6H\left( H\Phi +\dot{\Psi }\right) + 3\dot{H}\Phi \nonumber \\&\quad =\frac{1}{2F}\left[ 3\ddot{\delta F}+3H\dot{\delta F}-6H^2\delta F - 3\dot{F}\dot{\Phi }\right. \nonumber \\&\qquad - \left. 3\dot{F}\left( H\Phi +\dot{\Psi }\right) - \left( 3H\dot{F} + 6\ddot{F}\right) \Phi \right] , \end{aligned}$$
(3.30)
$$\begin{aligned} \ddot{\delta F}+ 3H \dot{\delta F} = \dot{F}(3H\Phi +3\dot{\Psi }+\dot{\Phi })+2\ddot{F}\Phi + 3H\dot{F}\Phi ,\nonumber \\ \end{aligned}$$
(3.31)
$$\begin{aligned} \delta F&= F'\delta R,\nonumber \\ \frac{F'}{F}R_{\mathrm {dS}}\delta R&= -2\left[ 3\left( \dot{H}\Phi +H\dot{\Phi }+\ddot{\Psi }\right) \right. \nonumber \\&+\left. 12H\left( H\Phi +\dot{\Psi }\right) +3\dot{H}\Phi \right] . \end{aligned}$$
(3.32)

Here, the term \(-(R/3)\delta F\) (the term \((F'/F)R_{\mathrm {dS}}\delta R\)) on the left hand side of Eq. (3.31) ((3.32)) disappears (appears) due to the redefinition of the scalaron mass squared (3.21).

All terms in Eqs. (3.24), (3.26), and (3.27) depend on time, and therefore may contribute to Eqs. (3.28)–(3.32). As we wrote above, according to our nonrelativistic approach, we neglect the terms of the order \(o(1/a^3)\). On the other hand, exponential functions decrease faster than any power function. Moreover, we can write the exponential term in Eq. (3.24) as follows:

$$\begin{aligned} \frac{\kappa ^2}{12\pi F'}\sum _i\frac{m_i\exp \left( - \sqrt{\frac{1 }{3}\left( \frac{F}{F'}-R_{\mathrm {dS}}\right) } |\mathbf{{r}}_{\mathrm {ph}}-\mathbf{{r}}_{\mathrm {ph}i}|\right) }{|\mathbf{{r}}_{\mathrm {ph}}-\mathbf{{r}}_{\mathrm {ph}i}|},\nonumber \\ \end{aligned}$$
(3.33)

where we introduced the physical distance \(r_{\mathrm {ph}} = a r\). It is well known that astrophysical tests impose strong restrictions on the nonlinearity [38, 39] (the local gravity tests impose even stronger constraints [35, 38, 39]). According to these constraints, Eq. (3.33) should be small at the astrophysical scales. Consequently, on the cosmological scales it will be even much smaller. So, we will not take into account the exponential terms in the above equations. However, in Eqs. (3.24), (3.26), and (3.27), we have also \(1/a^3\) and \(1/a\) terms which we should examine. Before performing this, it should be recalled that we consider the late Universe which is rather close to the de Sitter point. Therefore, as we already noted in the previous subsection, \(F\approx 1\), \(H\approx \mathrm {const}\ \Rightarrow \ \dot{H} \approx 0,\, R_{\mathrm {dS}}=12 H^2\) and \(F'(R_{\mathrm {dS}})\) is some finite positive value. Additionally, \(\dot{F},\ddot{F},\dot{F}' \sim 1/a^3\). Hence, all terms of the form of \(\dot{F},\ddot{F},\dot{F}' \; \times \; \Phi ,\Psi ,\dot{\Phi }, \dot{\Psi }\) are of the order \(o(1/a^3)\) and should be dropped. In other words, the functions \(F\) and \(F'\) can be considered as time independent.

First, let us consider the terms \(\Psi =\Phi = \varphi /a\) in Eqs. (3.26) and (3.27) and substitute them into Eqs. (3.28)–(3.32). Such \(1/a\) term is absent in \(\delta R\). So, we should put \(\delta R=0\), \(\delta F=0\). Obviously, this is the linear theory case. It can easily be seen that all equations are satisfied. Indeed, the functions \(\Phi \) and \(\Psi \) are included in Eqs. (3.28)–(3.30) and (3.32) (Eq. (3.31) is satisfied identically) in combinations \(H\Phi +\dot{\Psi }\) and \(H\dot{\Phi }+\ddot{\Psi }\) which are equal to zero.

Now, we study the terms \(\sim 1/a^3\), i.e.,

$$\begin{aligned} \begin{aligned}&\delta R=-\frac{\kappa ^2}{(F-F'R_{\mathrm {dS}})}\frac{\bar{\rho }_c}{a^3},\\&\Psi =-\frac{\kappa ^2F'}{2F(F-F'R_{\mathrm {dS}})}\frac{\bar{\rho }_c}{a^3},\\&\Phi =\frac{\kappa ^2F'}{2F(F-F'R_{\mathrm {dS}})}\frac{\bar{\rho }_c}{a^3}. \end{aligned} \end{aligned}$$
(3.34)

Let us examine, for example, Eq. (3.28). Keeping in mind that \(\delta F=F'\delta R\), one can easily get

$$\begin{aligned}&12H^2\frac{\kappa ^2F'}{2F(F-F'R_{\mathrm {dS}})}\frac{\bar{\rho }_c}{a^3}\nonumber \\&\qquad = 12H^2\frac{\kappa ^2F'}{2F(F-F'R_{\mathrm {dS}})}\frac{\bar{\rho }_c}{a^3}+o(1/a^3). \end{aligned}$$
(3.35)

Therefore, the terms \(\sim 1/a^3\) exactly cancel each other, and this equation is satisfied up to the adopted accuracy \(o(1/a^3)\). One can easily show that the remaining Eqs. (3.29)–(3.32) are fulfilled with the same accuracy.

Thus, we have proved that the scalar perturbation functions \(\Psi \) and \(\Phi \) in the form (3.26) and (3.27) satisfy the system of Eqs. (2.16)–(2.21) with the required accuracy. Both of these functions contain the nonlinearity function \(F\) and the scale factor \(a\). Therefore, both the effects of nonlinearity and the dynamics of the cosmological background are taken into account. The function \(\Phi \) corresponds to the gravitational potential of the system of inhomogeneities. Hence, we can study the dynamical behavior of the inhomogeneities (e.g., galaxies and dwarf galaxies) including into consideration their gravitational attraction and cosmological expansion, and also taking into account the effects of nonlinearity. For example, the nonrelativistic Lagrange function for a test body of the mass \(m\) in the gravitational field described by the metric (2.15) has the form (see [14] for details):

$$\begin{aligned} L \approx -m\Phi + \frac{m a^2\mathbf{{v}}^2}{2},\,\, \mathbf{{v}}^2=\dot{x}^2 + \dot{y}^2 + \dot{z}^2. \end{aligned}$$
(3.36)

We can use this Lagrange function for analytical and numerical study of mutual motion of galaxies. In the case of the linear theory, such investigation was performed, e.g., in [15]. With the help of the explicit expression (3.27) we can perform now similar numerical and analytical investigations for different \(f(R)\) models.

4 Conclusion

In our paper we have studied scalar perturbations of the metric in nonlinear \(f(R)\) gravity. The Universe has been considered at the late stage of its evolution and at scales much less than the cell of uniformity size which is approximately 190 Mpc [12]. At such distances, our Universe is highly inhomogeneous, and the averaged hydrodynamic approach does not work. Here, the mechanical approach [12, 14] is more adequate. Therefore, we have used the mechanical approach to investigate the scalar perturbations in nonlinear theories. We have considered a class of viable \(f(R)\) models which have de Sitter points in future with respect to the present moment [18, 19, 27].

The main objective of this paper was to find explicit expressions for \(\Phi \) and \(\Psi \) in the framework of nonlinear \(f(R)\) models. Unfortunately, in the case of nonlinearity the system of equations for scalar perturbations is very complicated. It is hardly possible to solve it directly. Therefore, we have considered the following approximations: the astrophysical approach, the large scalaron mass case and the quasi-static approximation. In all three cases we found the explicit expressions for the scalar perturbation functions \(\Phi \) and \(\Psi \) up to the required accuracy. The latter means that, because we consider nonrelativistic matter with the averaged rest mass density \(\bar{\rho }\sim 1/a^3\), all quantities in the cosmological approximation are also calculated up to the corresponding orders of \(1/a\). It should also be noticed that in the cosmological approach our consideration is valid for nonlinear models where functions \(f(R)\) have the stable de Sitter points \(R_{\mathrm {dS}}\) in future with respect to the present time, and the closer to \(R_{\mathrm {dS}}\) we are, the more correct our approximation is. All three popular models, Starobinsky [21], Hu–Sawicki [22] and MJWQ [23] (see also [17, 26]) have such stable de Sitter points in future (approximately at the redshift \(z=-1\)) [24, 25].

The quasi-static approximation is of most interest from the point of view of the large-scale structure investigations. Here the gravitational potential \(\Phi \) of Eq. (3.27) contains both the nonlinearity function \(F\) and the scale factor \(a\). Hence, we can study the dynamical behavior of the inhomogeneities (e.g., galaxies and dwarf galaxies) including into consideration their gravitational attraction and the cosmological expansion, and also taking into account the effects of nonlinearity. All this makes it possible to carry out the numerical and analytical analysis of the large-scale structure dynamics in the late Universe for nonlinear \(f(R)\) models.