$f(R)$ gravity: scalar perturbations in the late Universe

In this paper we study scalar perturbations of the metric for nonlinear $f(R)$ models. We consider the Universe at the late stage of its evolution and deep inside the cell of uniformity. We investigate the astrophysical approach in the case of Minkowski spacetime background and two cases in the cosmological approach, the large scalaron mass approximation and the quasi-static approximation, getting explicit expressions for scalar perturbations for both these cases. In the most interesting quasi-static approximation, the scalar perturbation functions depend on both the nonlinearity function $f(R)$ 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. Our investigation is valid for functions $f(R)$ which have stable de Sitter points in future with respect to the present time, that is typical for the most popular $f(R)$ models.

models. We consider the Universe at the late stage of its evolution and deep inside the cell of uniformity. We investigate the astrophysical approach in the case of Minkowski spacetime background and two cases in the cosmological approach: the large scalaron mass approximation and the quasi-static approximation, and get explicit expressions for scalar perturbations for both these cases. In the most interesting quasi-static approximation, the scalar perturbation functions depend on both the nonlinearity function f (R) and the scale factor a. Hence, we can study the dynamical behavior of the inhomogeneities (e.g., galaxies and dwarf galaxies) including into the consideration their gravitational attraction and the cosmological expansion, and also taking into account the effects of the nonlinearity. Our investigation is valid for functions f (R) which have stable de Sitter points in future with respect to the present time, that is typical for the most popular f (R) models.

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 the standard models, it has not still been offered a satisfactory explanation to these problems. 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][5][6]).
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). 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 metrics. This is approximately 190 Mpc and larger [7]. At these scales, the matter fields (e.g., cold dark matter) are well described by the hydrodynamical approach. However, at smaller scales the Universe is highly inhomogeneous. Here, the mechanical approach looks more adequate [7,8]. Therefore, it is of interest to study the compatibility of nonlinear theories of gravity with the mechanical approach. In other words, we examine equations for scalar perturbations 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 get the expressions for scalar perturbations. The scalar perturbation Φ is of special interest because it corresponds to the gravitational potential. We demonstrate that this function gives a possibility to take into account both the effects of nonlinearity of the original model and the dynamics of the cosmological background. It should be noted that the explicit form of this function makes it possible to carry out analytical and numerical study of mutual motion of galaxies in nonlinear models. In the case of the linear model, similar investigations were performed in the framework of the mechanical approach in the papers [8,9]. Therefore, our formulae 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 Sec. 2, we present the main background equations as well as 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 Sec. 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 Φ and Ψ up to the required accuracy. The main results are summarized in concluding Sec. 4.

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 where f (R) is an arbitrary smooth function of the scalar curvature R, S m is the action of matter, κ 2 = 8π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]) The trace of this equation gives In what follows, the 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 dS of the equation As it follows from Eq. (2.3), they are vacuum solutions (T = 0) of this equation for which the Ricci scalar is constant ( F (R) = 0). Such solutions are called de Sitter points [6]. We can expand the function f (R) in the vicinity of one of these points: where we used Eq. (2.4). Now, we suppose that parameters of the model can be chosen in such a way that Therefore, we get where Λ ≡ R dS /4. The stability of these points was discussed in [6,10]. Obviously, these models go asymptotically to the de Sitter space when R → R dS = 0 with a cosmological constant Λ = R dS /4. This happens when the matter content becomes negligible with respect to Λ as it is the case with a late Friedman-Robertson-Walker (FRW) cosmology. We can also consider a zero solution R dS = 0 of Eq. (2.4). It is more correct to call this point a Minkowski one. Here, Λ = 0, and such models go asymptotically to the Minkowski space.
In the case of the spatially flat background spacetime with the FRW metrics and matter in the form of a perfect fluid with the energy-momentum tensor components T µ ν = diag(−ρ,P ,P ,P ), Eq. (2.2) results in the following system: where the bar denotes the homogeneous background quantities, the Hubble parameter H = a/a (the dot everywhere denotes the derivative with respect to the synchronous time t) and the scalar curvature The perfect fluid satisfies the continuity equatioṅ which for nonrelativistic matter with P = 0 has the solution whereρ c = 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 Introduction, 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 [7]. 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 [7,8]. In 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 well described as point-like matter sources with the rest mass density where 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 [16] 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 metrics (2.8). In the conformal Newtonian (longitudinal) gauge, such perturbed metrics is [6] where scalar perturbations Φ, Ψ ≪ 1. It is worth noting that smallness of the nonrelativistic gravitational potentials Φ and Ψ 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 [16] (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 Φ. Then, we use this potential to determine dynamical behavior of galaxies. It gives us 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 [9]. Our present paper is devoted to the first step in this program. In other words, we are going to define scalar perturbations Φ, Ψ for the f (R) gravitational models. Under our assumptions and according to [6], these perturbations satisfy the following system of equations: (2.21) In these equations, the function F , its derivative F ′ and the scalar curvature R are unperturbed background quantities. Here, ∆ is a Laplacian in the comoving coordinates. As a matter source, we consider dust-like matter. Therefore, δP = 0 and

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.ρ = 0. As we mentioned above, there are two types of vacuum background solutions of Eq. (2.3): de Sitter spacetime with R dS = 12H 2 = const = 0 and Minkowski spacetime with R = 0, H = 0. However, the system of equations (2.16)-(2.21) was obtained for the FRW metrics (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 (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 (2.14). Keeping in mind thatρ = 0, we have δρ = ρ.
In the case of Minkowski spacetime background and dropping the time derivatives, Eqs. (2.16)-(2.21) in the astrophysical approach are reduced to the following system: where the function ϕ satisfies the equation Here, we took into consideration that in the astrophysical approach δρ c = ρ c where ρ c is defined by (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 eff = G N /F . Eq. (3.2) follows directly from (3.6) and, consequently, may be dropped, while from (3.4) we get the following Helmholtz equation with respect to the scalaron function δR: On the other hand, it can be easily seen that the substitution of Eqs. (3.6) and (3.7) into (3.5) results in the same Eq. (3.8). Therefore, in the case of Minkowski background, the mass squared of the scalaron is 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 [17][18][19][20] (see also Eq. (5.2) in the review [6]).

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ρ ∼ 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, δρ ∼ 1/a 3 [8].

Large scalaron mass case
As we can see from (3.9), the limit of the large scalaron mass corresponds to F ′ → 0. Then, δF is also negligible (see (2.21) where the introduced function ϕ depends only on spatial coordinates. Substituting (3.16) into (3.10), we obtain 1 As we mentioned above, neglecting relativistic matter in the late Universe, we have δρ ∼ 1/a 3 [8]. This approximation is getting better and better performed in the limit a → +∞. 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 → R ∞ , where the value R ∞ is just finite. Then, from (3.17) we immediately come to the condition F = const + o(1) , (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Ḟ ∼ 1/a+o(1/a). However, as we will see below, this is wrong. Obviously, without loss of generality we can suppose that const = 1. From the condition (3.18) we get where Λ 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 ∞ may correspond to the de Sitter point R dS in future (see Eq. (2.7)). As we have written in Sec. 2, all three popular models, Starobinsky [11], Hu-Sawicki [12] and MJWQ [13], have such stable de Sitter points in future (approximately at the redshift z = −1) [14,15]. The condition of stability is 0 < RF ′ /F < 1 (see, e.g., (4.80) in [6]). Since F ≈ 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 (3.13) which after the substitution of (3.16) and (3.17) and some simple algebra takes the form To estimateḞ andF , we take into account that in the limit R → R ∞ , F ≈ 1, H ≈ const ⇒ H ≈ 0, and F ′ (R ∞ ) is some finite positive value. Then,Ḟ = F ′Ṙ ∼ F ′ (R ∞ )Ṙ ∼Ṫ ∼ d 1/a 3 /dt ∼ H 1/a 3 ∼ 1/a 3 andF ∼ȧ/a 4 ∼ 1/a 3 . Therefore, the left hand side of Eq. Thus, in the case of the large enough scalaron mass we reproduce the linear cosmology from the nonlinear one, as it should be.

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 δ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 [21,22] (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 Φ, Ψ) are produced mainly by the spatial distribution of astrophysical/cosmological bodies. As we have seen, Eqs. (3.1)-(3.5) result in (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 dS in future 1 . This means that δρ = ρ −ρ in (3.7), all background quantities (e.g., F, F ′ ) are calculated roughly speaking at R dS and the scalaron mass squared (3.9) reads now [6,[17][18][19]] Let us consider now Eq. (3.8) with the mass squared (3.21). Taking into account that now δρ c = ρ c −ρ c , we can rewrite this equation as follows: where (3.23) Eq. (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 δR ∼ 1/r at small distances r and δR → 0 for r → ∞. Taking all these remarks into consideration, we obtain for a full system It is worth noting that averaging over the whole comoving spatial volume V gives the zero value δR = 0. Really, since i m i /V =ρ c , This result is reasonable because the rest mass density fluctuation δρ, representing the source of the metric and scalar curvature fluctuations Φ, Ψ and δR, has a zero average value δρ = 0.
Consequently, all enumerated quantities should also have zero average values:Φ =Ψ = 0 and δR = 0, in agreement with (3.25). From Eq. (3.6) we get the scalar perturbation functions Ψ and Φ in the following form: where ϕ satisfies Eq. (3.7) with δρ in the form (2.22) (i.e.ρ c = 0). Obviously, when F ′ → 0, M → ∞, and we have exp (−M |r − r i |) /|r − r i | → 4πδ(r − r i )/M 2 , so the expression in the square brackets in (3.26) and (3.27) is equal to κ 2 δρ c / (F − F ′ R dS )a 3 . Therefore, in the considered limit F ′ → 0 we reproduce the scalar perturbations Φ, Ψ from the previous large scalaron mass case, as it certainly should be. Thus, neglecting for a moment the influence of the cosmological background, but not neglecting the scalaron's contribution, we have found the scalar perturbations. They represent the mix of the standard potential ϕ/a (see the linear case [8]) 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 (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: 3 Ḣ Φ + HΦ +Ψ + 6H HΦ +Ψ + 3ḢΦ = Here, the term −(R/3)δF (the term (F ′ /F )R dS δR) in the left hand side of (3.31) ((3.32)) disappears (appears) due to the redefinition of the scalaron mass squared (3.21). It can be easily seen that all terms in (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 (3.24) as follows: where we introduced the physical distance r ph = ar. It is well known that astrophysical tests impose strong restrictions on the nonlinearity [23,24] (the local gravity tests impose even stronger constraints [20,23,24]). According to these constraints, (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 (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 ≈ 1, H ≈ const ⇒Ḣ ≈ 0, R dS = 12H 2 and F ′ (R dS ) is some finite positive value. Additionally,Ḟ ,F ,Ḟ ′ ∼ 1/a 3 . Hence, all terms of the form ofḞ ,F ,Ḟ ′ × Φ, Ψ,Φ,Ψ 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 Ψ = Φ = ϕ/a in Eqs. (3.26) and (3.27) and substitute them into Eqs. (3.28)-(3.32). Such 1/a term is absent in δR. So, we should put δR = 0, δF = 0. Obviously, this is the linear theory case. It can be easily seen that all equations are satisfied. Indeed, the functions Φ and Ψ are included in Eqs. (3.28)-(3.30) and (3.32) (Eq. (3.31) is satisfied identically) in combinations HΦ +Ψ and HΦ +Ψ which are equal to zero. Now, we study the terms ∼ 1/a 3 , i.e.
Let us examine, for example, Eq. (3.28). Keeping in mind that δF = F ′ δR, one can easily get (3.35) Therefore, the terms ∼ 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 Ψ and Φ 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 Φ 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 metrics (2.15) has the form (see [8] for details): 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 [9].

Conclusion
In our paper, we have studied scalar perturbations of the metrics 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 [7]. At such distances, our Universe is highly inhomogeneous, and the averaged hydrodynamic approach does not work here. 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 Φ and Ψ. The function Φ represents the gravitational potential of inhomogeneities. Then, we use the explicit expression for Φ to describe the motion of inhomogeneities. Such mechanical approach is well known in astrophysics (see, e.g., [16]). We generalized it to the case of dynamical cosmological background [7,8].
The main objective of this paper was to find explicit expressions for Φ and Ψ 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 Φ and Ψ up to the required accuracy. The latter means that, because we consider nonrelativistic matter with the averaged rest mass densityρ ∼ 1/a 3 , all quantities in the cosmological approximation are also calculated up to the corresponding orders of 1/a. It should be also noted that in the cosmological approach our consideration is valid for nonlinear models where functions f (R) have the stable de Sitter points R dS in future with respect to the present time, and the closer to R dS we are, the more correct our approximation is. All three popular models, Starobinsky [11], Hu-Sawicki [12] and MJWQ [13] have such stable de Sitter points in future (approximately at the redshift z = −1) [14,15].
The quasi-static approximation is of most interest from the point of view of the largescale structure investigations. Here, the gravitational potential Φ (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.