1 + 3 covariant perturbations in power-law f(R) gravity

We applied the 1+3 covariant approach around the Friedmann–Lemaître–Robertson–Walker (FLRW) background, together with the equivalence between f(R) gravity and scalar-tensor theory to study cosmological perturbations. We defined the gradient variables in the 1 + 3 covariant approach which we used to derive a set of evolution equations. Harmonic decomposition was applied to partial differential equations to obtain ordinary differential equations used to analyse the behavior of the perturbation quantities. We focused on dust dominated area and the perturbation equations were applied to background solution of αR+βRn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha R+\beta R^{n}$$\end{document} model, n being a positive constant. The transformation of the perturbation equations into redshift dependence was done. After numerical solutions, it was found that the evolution of energy-density perturbations in a dust-dominated universe for different values of n decays with increasing redshift.


Introduction
The standard model of cosmology, ΛCDM, is by far the most successful model in explaining the origin of cosmic microwave background (CMB) radiations, the observations of Type I a supernovae, the formation and distribution of large-scale structures, the synthesis of light elements in the universe and the expansion of the universe [1][2][3][4].
Although ΛCDM fits most of current cosmological data, it does not describe for instance, the horizon problem the flatness problem and the structure (homogeneity/smoothness) problem. It also does not describe the mechanism for dark energy production and the current acceleration of the universe. Some of these cosmological problems emerge in the early universe and others in the late universe [5][6][7][8].
The greatest challenge to the standard model of cosmology came when observation of type I a supernovae showed that we live in an accelerated expanding universe where spacea e-mail: murorunkwereb@gmail.com (corresponding author) time between galaxies is continuously growing. This epoch in the evolution of the universe can not be explained by the presence of baryonic matter. The observational data showed that 4% of the energy content of the universe is baryonic matter, 23% is dark matter and 72% is dark energy that is believed to be the root cause of the cosmic expansion of the universe [9,10].
Different theories have been constructed in attempt to explain the above mentioned problems such as exotic scalar fields, cosmological constant and modified theories of gravity [11,12]. There have been several proposed modified theories of gravity such as theories of gravity with extra fields, higher-dimensional theories of gravity and higher-order theories of gravity [13]. One of the classes of higher order theories that explains the shortcomings of general relativity (GR) are fourth-order gravity models and they are generally obtained by including higher-order curvature invariant in the Einstein-Hilbert action. Some of the higher order theories of gravity are f (R) theories of gravity in which the Hilbert-Einstein action of general relativity is modified where R represents the Ricci scalar and f (R) is the arbitrary function of Ricci scalar [14]. This way of modification was first instituted by Buchdahl in 1970s, see his original work in [15]. The f (R) theories are the most explored alternatives to general relativity GR in the context of late-time accelerated expansion of the universe and they became popular due to their ability of producing late-time acceleration of the universe without introduction of cosmological constant [16,17].
Evolution of small density perturbations lead to large scale structures in the universe. The 1 + 3 covariant approach is chosen to study the evolution of the scalar perturbation equations and it has been introduced by Ellis and Bruni in 1980s [18]. One of the importance of this approach is that, it allows one to define gradient variables that describe the inhomogeneity and anisotropy of the universe. The 1 + 3 covariant approach has been applied to study scalar perturbations in f (R) gravity [9,[19][20][21]. It is also studied in [22], where the approach was used to explore the equivalence between scalar-tensor theory and f (R) gravity but they used a different definition of the scalar field.
Linear perturbations have severally been discussed in f (R) models in different studies using different approach [19,[21][22][23][24][25][26]. In these studies, the investigation between dark energy and f (R) models was done by analysing the behavior of scalar perturbations around FLRW background. They found that for all values of positive n of R n model, there is always a growing mode for energy density perturbations. In particular, the α R + β R n model has been analysed at the level of background and perturbation using many different approaches [19,20,24,27,28].
Different types of scalar fields have been used in modified theories of gravity for different studies, for instance the study of perturbations in scalar-tensor theory, focusing on the Brans-Dicke theory of gravity as subclass of scalar-tensortheory of gravity was done using the scalar field φ = f − 1 [22], the same authors studied f (R) theory in scalar-tensor theory of gravity limiting only to inflationary parameters using the scalar field φ = f − 1 [29].
The investigation of the relationship between scalar-tensor theory and f (R) theories of gravity for palatini formalism using the scalar field φ = ln f was done [30] and the study of how to solve higher field equations in higher curvature gravity theory using the scalar field φ = ln f was done [31]. However, the study of 1+3 covariant perturbations using the scalar field φ = ln f in α R + β R n model has not been explored.Therefore, there is a need to study if the scalar field φ = ln f at perturbation level produce also significant results.
In this work, we study the evolution of the matter energydensity fluctuations in the universe using the 1 + 3 covariant approach and their effect on structures formation (galaxies, clusters and superclusters) by considering φ = ln f . Using the gradient variables defined in 1 + 3 covariant approach, we will derive the perturbation equations and subject them to α R + β R n model. We will study the dynamics of the matter energy-density perturbations both in the short-and longwavelength limits of α R + β R n model of gravity. This paper is organised as follows: In Sect. 2, we describe the 1+3 covariant approach. In Sect. 3, f (R) theory and scalar-tensor theory and their equivalence are reviewed. In Sect. 4, gradient variables are defined. In Sect. 5, evolution equations are generated. Sect. 6 is devoted for harmonic decomposition. In Sect. 7, we look for the solutions in α R + β R n theory of gravity, In Sect. 8, we provide discussion. Section 9 summarises our conclusion. The adopted spacetime signature is (− + ++) and unless stated otherwise. The symbol ∇ refer to covariant derivative, ∂ is partial differentiation and the over dot shows differentiation with respect to cosmic time.

The 1 + 3 covariant approach
The 1 + 3 covariant approach is a framework for studying the linear evolution of the cosmological perturbations. In this approach, a fundamental observer divides space-time into hyper-surfaces and a perpendicular 4-velocity field vector where 1 + 3 indicate the number of dimensions involved in each slice. One of the importance of the 1 + 3 covariant approach is to identity a set of covariant variables which describe the inhomogeneity and anisotropy of the universe .
The 4-velocity field vector u a is defined as [22,[32][33][34] where τ is the proper time such that u a u a = −1. From the 4-velocity u a , we construct a projection tensor that projects perpendicular to u a and is given by where g ab is the metric of the spacetime. The projection tensor that projects parallel to u a is given by The covariant time derivative along the fundamental wordlines is defined aṡ The kinematic quantities that define the dynamics of spacetime, expansion, shear and vorticity are obtained by splitting ∇ a u b into irreducible parts, as follows [35] ∇ a u b =∇ a u b − u aub where w ab = ∇[au b ] is the vorticity tensor which satisfies u a w ab = 0 which describes the rotation of the matter relative to a non-rotating frame, σ ab is the shear tensor which describes the rate of distortion of the matter flow, θ = ∇ a u a = 3H is the rate of expansion of the fluid where H = θ 3 =ȧ a , H is the hubble parameter, andu a is the acceleration.
The volume element for 3-restspaces is given by where η abcd is the 4-dimensional volume element. The matter stress energy tensor that describes the matterenergy content of the universe is decomposed irreducible with respect to u a as [35,36] T ab = ρu a u b + q a u b + q b u a + Ph ab + π ab , where ρ is the energy-density measured by the observer moving with 4-velocity u a , q a is the relativistic energy flux (momentum density), P is the anisotropic pressure of the matter and π ab is the trace-free anisotropic pressure.
In the FLRW limit , q a and π ab vanish since the stress energy tensor is restricted to perfect fluid form. The tensor T ab then reduces to where the equation of state for perfect fluid is p = p(ρ).

The f (R) theory of gravity
The action of f (R) theory of gravity is given by the following equation [30,37] where k = 8π G, G is the gravitational constant, R is the Ricci scalar, L m is the Lagrangian that contains matter (normal matter, CDM and radiation), g μν is the metric, ψ is the matter fields, g is the determinant of the metric and f is the general differentiable function of the Ricci scalar. The arbitrary function introduced in the Lagrangian gives us more freedom to explain the observed cosmic acceleration and formation of large-scale structures without inclusion of the cosmological constant. The variation of the action given by Eq. (9) with respect to the metric, gives where , R μν is the Ricci tensor and ∇ c ∇ c = is the D'Alembert operator.
The term T m μν represent the stress energy tensor of standard matter and the Eq. (10) reduces to the standard Einstein field equation when f (R) = R.
The general action that represents the scalar-tensor theory is given by [30] where U (φ) is the potential of the scalar field φ, y(φ) and φ , the scalar-tensor theory reduces to Brans-Dicke theory of gravity as For vanishing coupling parameter, W = 0, the Brans-Dicke theory reduces to For a given f (R) lagrangian, one can define an auxiliarly field χ such that it is a function of the scalar field φ as χ(φ) so that the action of f (R) gravity given by Eq. (9) becomes where variation with respect to Setting the potential as the action of f (R) gravity given by Eq. (14) takes the form The action given by Eq. (16) is equal to the action of Brans-Dicke theory for vanishing coupling parameter given by Eq. (13), hence f (R) is a special case of the scalar-tensor theory.

The gradient variables
The following are gradient variables defined in 1 + 3 covariant approach. These variables are used in analysis of the cosmological perturbation and they are defined as follow [19,36,38]: The gradient variable responsible for perturbation in scalar field is given by where a is the FLRW cosmological scale factor that shows the dynamics of the universe with time and∇ is 3-spatial gradient.
The gradient variable that characterises perturbation due to the momentum of the scalar field is given by The comoving volume expansion gradient variable is given by The comoving fractional gradient in the energy density is given by The quantities with index m refer to matter (dust and radiation) and these gradient variables contains both scalar and vector variables.
The scalar variables are believed to be the ones that lead to the structure formation that we see in the universe such as galaxies, clusters and super-clusters and they are obtained by taking the divergence of the vector quantities.
The local decomposition method is defined as [19] where X ab describes shear and X [ab] describes vorticity and they describe the rotation of the fluid relative to a non-rotating frame and distortion of the fluid flow respectively. Using this decomposition method, the scalar variables are defined as: The background under consideration is FLRW spacetime; and the background quantities: energy-density and isotropic pressure are given by [22]: For a scalar field, the energy-density and isotropic pressure are given as The four velocity is given bẏ These background quantities together with the gradient variables are used to derive 1+3 covariant perturbation equations for the universe.

The 1 + 3 covariant perturbation equations
The derivation of the perturbation equations around an homogeneous and isotropic background is done by derivating gradient variables with respect to cosmic time. Using Eqs. (26)(27)(28)(29)(30), one can derive the evolution equations of the vector variables (17)(18)(19)(20) as follow: The linear evolution equation for the gradient variable responsible for the perturbation in the energy-density is given bẏ The Eq. (31) can be found in the existing literatures [9,19,20]. The linear evolution equation for the gradient variable responsible for the perturbation in the volume expansion is given bẏ The linear evolution equation for the gradient variable responsible for perturbation in the scalar field is given bẏ The linear evolution equation for the gradient variable that characterises the perturbation due to momentum of the scalar field is given bẏ The system of Eqs. (31)(32)(33)(34) describe the evolution equations of the vector gradient variables and they determine the evolution of every single perturbation variable. The scalar evolution equations that characterise the evolution of the spherical symmetry of the gradient variables are given by taking time derivative of Eqs. (22)(23)(24)(25) and they are given bẏ From the scalar evolution equations, we get their second order evolution equations by derivating them with respect to cosmic time.
The second order evolution equation for the gradient variable responsible for perturbation in the scalar field is given by the time derivative of Eq. (36) as The second order evolution equation for the gradient variable responsible for the perturbation in the energy-density is given by time derivative of Eq. (35) as These second order differential Eqs. (39) and (40) are to be analyzed after getting their harmonic decomposed equations by considering different scales.

Harmonic decomposition
Harmonic decomposition analysis is a way of reducing partial differential equations into ordinary differential equations which are easier to be solved [22,36]. The evolution Eqs. (39) and (40) can be thought as a coupled system of harmonic oscillators of the form [36,39]Ẍ where the A represents the friction (damping) term, the B represents the restoring force while C represents the source forcing term. Harmonic analysis applies the separation of variables techniques such that and write where Q k are the eigenfunctions of the covariantly defined Laplace-Beltrami operator on FLRW spacetime that is given bỹ The wave-number given by k = 2πa λ and λ is the wavelength of the mode. In this way the evolution Eqs. (35)(36)(37)(38)(39)(40) can be converted into ordinary differential equationṡ We note that Eqs. (47-49) can be found in the existing literature [22]. The Eqs. (51) and (52) are the ones responsible for evolution of linear perturbations of single content. The harmonic decomposed quantities are associated with the wave-number k with the corresponding wavelength, so we can make analysis by considering different scale such as short-wavelength and long-wavelength limit.

The case of f (R) = α R + β R n
Consider the case where f (R) theory of gravity is defined as where α and β are some constants parameter, n is a positive constant and H 0 is the present Hubble parameter. This f (R) theory is characterised by the following action This type of f (R) theory of gravity is the generalization of both GR and R n action. If α = 1 and β = 0, this reduces to GR theory of gravity and for α = 0 and β different from 0, it reduces to R n gravity. This model has been analyzed using the dynamical system approach [40] and proved that like R n gravity, it has unstable fixed point associated with the Friedmann like solution given by the following equation.
where w is the equation of state parameter, m is a positive constant and n = m . This phase was argued to be suitable for the structure formation to take place and it is for this reason we used Eq. (55) as background. The expression for the volume expansion, the Ricci scalar and matter energy density in this model are respectively given by [41] The f (R) = α R + β H 2 0 R H 2 0 n model has been used in different studies using different approach [40,42] and it is a fourth order gravity theory with additional scale for the gravitational interaction so we are going to apply perturbation equations to this background and check the behavior of the solutions.

Long-wavelength limit
In the analysis of the evolution of scalar perturbation equations in long-wavelength limit, the wave-number k is considered to be small so that the wavelength λ = 2πa k associated with it is much larger than the Hubble radius [9]. From Eq. (46), it implies that all laplacian term can be neglected and the spatial dependence of the perturbation variable can be factored out. So In this limit, we study the evolution of energydensity perturbation equations by considering scales which have wavelength much larger than the Jean's wavelength that is to say λ λ J .

Dust-dominated universe
After the epoch of radiation, the Universe was dominated by pressureless matter (dust) fluid. In such epoch, the density of dust is much greater than the density of the radiation (ρ d >> ρ r ). Due to the fact that the matter content in the universe is dominated by the dust, the equation of state parameter w = 0, since dust is pressureless. Consequently Eqs. (51) and (52) read: To transform any time derivative function f and H into redshift dependence, it follows aṡ where N = ln a and a is the scale factor that is related to cosmological redshift as We apply the Eqs. (63) and (64) to perturbation equations to get them in redshift dependence. Then we have perturbation Eqs. (59) and (60) in redshift space as For the case of f (R) = R − 2Λ, known as Lagrangian for general relativity with a cosmological constant, the Eq. (66) reduces to [1]: where the density parameter of cosmological constant fluid is given by The Eq. (67) admits the following solution When Λ = 0, this reduces to the original field equations of general relativity, then Eq. (67) reduces to [43]: and it admits the following solution .
(71) C 1 , C 2 , C 3 and C 4 are constants of intergration and they are obtained using some initial conditions. A graphic representation of the behavior of the Eq. (66) as n changes for f (R) theory of gravity is given by Fig. 1. We plot the normalized dust density on the vertical axis δ(z) = versus redshift where z 0 = 2000 is the initial redshift.

Short-wavelength limit
In small scale limit, we study the evolution of the shortwavelength modes, that is to say large values of the wavenumber k using the equations for scalar field fluid-dust sys- tem. The quasi-static approximation for the matter perturbations is used. In this approximation all time derivative of gradient variable responsible for perturbation in scalar field and momentum of the scalar field are discarded that is to saẏ Φ(z) = 0 andΨ (z) = 0 and only those including energydensity perturbations are kept. This approximation method has been widely used in different literatures such as [44,45].
In this limit, we study the evolution of energy-density perturbation equations by considering scales which have wavelength much smaller than the Jean's wavelength that is to say λ λ J .

Dust-dominated universe
If, we assume that the universe is dominated by dust as matter fluid, then the barotropic factor of dust w d = 0. Applying quasi-static approximation method whereΦ(z) = 0 anḋ Ψ (z) = 0 and the equation of state parameter w d = 0 to the Eq. (51) giveṡ Substituting the value of Φ(z) given by Eq.(72) into the energy-density perturbation equation given by Eq. (52), the following expression is obtained In redshift space, the Eq. (73) becomes The plots in Fig. 2 show the dynamics of energy density fluctuations of dust obtained by solving numerically the full systems Eq. (69) for ΛCDM case and (74) for α R + β R n model. We plot the normalized energy-density of dust on the vertical axis δ(z) = versus redshift.

Discussions
The evolution of energy-density perturbations is studied by solving numerically different perturbation equations for both the long-wavelength and short-wavelength limits. And dust dominated universe was taken into account. In the long-wavelength limit, the wave-number k is so small that the k 2 terms are neglected because they are too small compared to other terms. The wavelength λ = 2πa k λ J where λ J is the Jean's length given by [46] λ J = C s π Gρ , C s is the speed of sound, G is the gravitational constant and ρ is the energy-density. Thus, we considered the scales of the fluid inhomogeneity which are larger than the Jean's length. In the short-wavelength limit, we considered the scales which are less than the Jean's length that is λ = 2πa k λ J and quasi-static approximation has been used in this limit. In the small scales limits, we have fixed the wavelenght λ to be less than the Jean's wavelength λ J as λ = 1.5 Mpc.
We obtained the numerical solutions of energy-density perturbations equations for different values of the parameter n such as n = 1.1, n = 1.2, n = 1.3 for α R +β R n model and n = 1 for ΛCDM case. The choice of n was made according to [9,20].
We solved the perturbation equations setting the initial values to beΔ(z 0 ) = 0, Δ(z 0 ) = 0.00001,Φ(z 0 ) = 0 and Φ(z 0 ) = 0.000001 where z 0 is the initial redshift when the radiation and dust were equal in the universe, that is approximately z 0 = 2000. We computed the numerical solutions of the energy-density perturbation equations for α R + β R n model using α = 0.01 in the long-wavelength mode and α = 0.001 in the short-wavelength model for different values of the constant β. The density parameter Ω d = 0.32 and the hubble parameter H 0 = 73 km s −1 M pc −1 were used to get numerical solutions.
The scalar field φ = ln f has been used as extra degree of freedom to develop perturbation equations. The prime indicates derivative of f (R) with respect to R, knowing that in GR case, we do not have scalar field since f = 1.
The evolution of the energy-density perturbation Eq. (66) in the long-wavelength limit and the Eq. (74) in the shortwavelength limit for α R + β R n model in dust-dominated universe together with Eq. (69) of Λ CDM case is plotted in Figs. 1 and 2 respectively. From these figures, we found that the perturbations of the energy density decrease with increasing redshift for different values of constant n. The figures shows that the evolution of density perturbations for f (R) models is above the one of GR, which indicates that the modified gravity exhibits the overdensity referring to the one produced by GR. It appears that we have an increase in the energy-density perturbations as we move from early universe to late universe and this implies that there are more structure formations (galaxies, clusters and super-clusters) and more distanced galaxies in the present universe. This is in agreement with the existing literature about f (R) models [9,19,22,35]. The formation of structures in the universe is one of the characteristics of cosmic expansion, because of this, one could observe the advantage of having more degree of freedom in considering f (R) theory than in GR.

Conclusions
In this work, we presented a detailed analysis of 1+3 covariant approach where the universe is described by dust-fluid system . We used the equivalence between metric f (R) gravity and scalar-tensor theory with the definition of the scalarfield φ = ln f . We developed 1 + 3 covariant perturbation equations based on the gradient variables. We applied Harmonic decomposition method to partial differential equations to obtain ordinary differential equations for easier analysis. The background solutions of α R + β R n model was taken into account for fluid system dominated by dust and checked for solutions in short and long-wavelength limit. We used quasi-static approximation in the short-wavelength limit. We transformed time dependent perturbation equations into redshift dependence and we obtained numerical solutions using Maple software. Both in the short-and long-wavelength limits, we found that perturbations of the energy-density in dusts decay exponentially with increasing redshift for different values of n of α R + β R n model.
As an extension to this work however, further research is required on multi-fluid systems and other forms of f (R) gravity models.
The Equations (103-106) can be found in the existing literature [22].
To deal with the evolution of projected gradients we use the identity equation