Dark effects in $f(R,P)$ gravity

In the present paper a new cosmological model is proposed by extending the Einstein--Hilbert lagrangian with a generic functional $f(R,P)$, which depends on the scalar curvature $R$ and a term $P$ which encodes a possible influence from specific cubic contractions of the Riemann tensor. After proposing the corresponding action, the associated modified Friedmann relations are deduced, in the case where the generic functional has a power law decomposition, $f(R,P)=f_0 R^n+g_0 P^m$, with $f_0, g_0, n, m$ constant parameters. In this case the specific method of dynamical system analysis is employed, revealing the fundamental properties of the phase space structure, discussing the dynamical consequences for the cosmological solutions obtained. It is revealed that the cosmological solutions associated to the critical points can explain various dynamical eras, with a high sensitivity to the values of the $n$ and $m$ parameters, encoding various effects due to the geometrical nature of the specific couplings.


I. INTRODUCTION
The accelerated expansion of the Universe represents an important fundamental question in modern physics. The solution to this question can offer new insights for the gravitational theory, leading to a golden era for the cosmological context, opening various directions in modern physics. The phenomenon behind the accelerated expansion of the Universe represents an enigma, having different implications to the development of science and technology. This phenomenon has been studied intensively in the past two decades, being supported by various different observational analyses [1][2][3][4][5][6][7][8].
The modified gravity theories [9][10][11][12][13] represent an attempt of correcting the basic ΛCDM model, a theoretical approach which aims towards a more complete framework capable of solving different inconsistencies [14][15][16][17][18][19] of the present cosmological scenarios. These theories are based on specific modifications or replacements of the Einstein-Hilbert action, leading to a dynamical evolution of the dark energy equation of state. The main aim of these theories is to explain the known evolution of the Universe from early times to the present days. Within these theories a particular direction has been established recently which include higher order terms in the corresponding action, a specific model which further replaces or complements the basic gravitational theory based on the Einstein-Hilbert action [20].
The higher order gravities [20] represent an alternative approach in the modified gravity theories, capable of explaining various aspects and effects for the gravitational interaction. In this framework the Einsteinian cubic gravity was proposed [21], a specific theory based on specific contractions of the Riemann tensor in the third order. After this direction has emerged, various authors have investigated different properties of the Ein- * mihai.marciu@drd.unibuc.ro steinian cubic gravity or in similar approaches, analyzing the black-holes solutions [22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39]. Moreover, for this scenario, the wormholes properties have been addressed in [40,41]. In general the cubic gravity term is expected to have a considerable influence especially at early times during the inflationary epoch [42]. The effects due to the cubic term in the inflationary epoch have been studied intensively in the last years [42][43][44][45][46]. In Ref. [47] the authors have investigated the inflationary epoch by including a scalar field in a cubic gravity theory, showing the viability of such an epoch in this scenario. The basic aspects for thermodynamics in the case of a generic cubicquartic gravity have been investigated recently [48]. The fundamental properties of the anisotropic instabilities in the cubic gravity have been also analyzed [49]. In this regard, it has been shown that in the Einsteinian cubic extension various specific pathological instabilities might emerge. Recently, in Ref. [50] various pathological aspects have been analyzed in the case of Einsteinian cubic gravity and generalised quasi-topological gravity models.
The extended cubic gravity based on the third order contractions of the Riemann tensor has been proposed recently [51], a theory which further corrects the Einstein-Hilbert action with a generic functional f (P ) which encodes specific effects due to the topological invariant P . In this case the topological invariant P is based on some specific contractions of the Riemann tensor [51]. The dynamical features for the generalized extended cubic gravity have been addressed in [52], by taking into account two possible configurations for the cubic term which corrects the Einstein-Hilbert action. An alternative proposal has been analyzed in [53] for the Einsteinian cubic gravity, a specific theory which also includes a cosmological constant. The inclusion of the cubic invariant into possible theories which are embedding scalar fields has been investigated in [54], considering a specific coupling for the quintessence or phantom models. Furthermore, different studies have considered a holographic approach for the Einsteinian cubic gravity [55,56]. Recently, the Starobinsky's model of inflation has been corrected by arXiv:2103.08420v2 [gr-qc] 4 Dec 2021 adding a cubic component [57]. The authors have investigated the slow-roll regime and discussed the possibility of validation, taking into consideration the scalar and tensor perturbations.
In the present paper we shall further generalise the cubic extension of the Einstein-Hilbert action [51] by including viable effects from the curvature. Hence, in this approach we shall add to the Einstein-Hilbert action the functionalf (R, P ), a generalization which takes into account possible geometrical effects due to the curvature and the inclusion of the cubic invariant. In this case the physical features corresponding to the present model shall be investigated by considering specific methods associated to the dynamical system analysis [58]. This methods have been considered in many cosmological scenarios, representing important analytical tools in the study of physical systems [58].
The plan of the present paper is the following. In Sec. II we shall describe the action for the current cosmological model, taking into account a simple decomposition for the generic extensionf (R, P ). After writing the action, we shall present the associated modified Friedmann relations, obtained by the variation of the action with respect to the inverse metric. Then, in Sec. III we adopt the specific power-law parameterizationf (R, P ) = f 0 R n + g 0 P m with f 0 , g 0 and n, m as constant parameters. For this specific model we propose the special form of the auxiliary variables, considering the dynamical system analysis in the case of the powerlaw parameterization. In Sec. IV we analyze the exponential decompositionf (R, P ) = f 0 e nR + g 0 e mP , with f 0 , g 0 , n, m constant parameters. After we describe the main features of the phase space for the two models we present a short summary and the final concluding remarks in Sec. V.

II. THE ACTION AND THE FIELD EQUATIONS
In the present study we shall propose a new cosmological model which is described by the following action: We note that the Einstein-Hilbert action is further extended by adding a generic functionalf (R, P ) which depends on the scalar curvature R and an additional term P which encodes specific cubic contractions of the Riemann tensor [21]. In this expression,g represents the determinant of the metric, while S m is the action corresponding to the matter sector. Notice that in our action (1) the radiation component is neglected since we are interested in late-time dynamics. The additional term P embedded here has the following decomposition [21,51]: encoding specific effects due to the cubic contractions of the Riemann tensor. In this formula the β i (i = 1, 8) components are constant parameters. In what follows we shall consider the case of a Robertson-Walker metric: where a(t) represents the corresponding scale factor, and t the specific cosmic time. In this case the Hubble parameter will be denoted with H = 1 a da dt . Note that in our approach we have set the spatial curvature index to zero, a specific value which is compatible with various astrophysical observations. Next, we shall consider the following specific relations between the constant parameters β i , i = 1, 8 [21,51], β = (−β 1 + 4β 2 + 2β 3 + 8β 4 ).
In this specific case the cubic term is equal to the following expression [51], describing a topological invariant in the four dimensional space-time. Furthermore, we shall implement the following decomposition for the specific extension of the Einstein-Hilbert action, by considering In the previous formula we shall encode the effects due to the curvature couplings in the f (R) part [59], while the cubic couplings are embedded into the behavior of the g(P ) functional. The scalar curvature in the case of the Robertson-Walker metric (3) is equal to R = 6(2H 2 +Ḣ).
For the matter component, which is assimilated into the dark matter sector, the energy-momentum tensor is defined as: If we take into account the Robertson-Walker metric (3), we have the following representation for the matter energy-momentum tensor: with ρ m the density and p m the pressure, connected through a barotropic equation of state where w m is a constant parameter which describes the properties of the dark matter component. For simplicity we shall consider the dust case where w m = 0. This implies that the the matter component can be regarded as a pressure-less gas which is not too dense.
If we consider the variation of the action (1) with respect to the inverse metric g µν , for the previous men-tioned decomposition (9), we obtain the following modified Friedmann relations: The densities ρ f , ρ g and pressures p f , p g have the following expressions [51]: For the present cosmological model the dark energy is encoded into the curvature and cubic couplings, having a geometrical nature. From the definitions of the energy densities and pressures, we note that the current scenario satisfies the continuity equation for the matter and dark energy sector. In the above equations ∂ t and the dotḋ efines the partial derivative with respect to the cosmic time, while ∂ 2 t represents the double time derivative. In the case where f (R) = 0 the present action (1) describes the so-called cubic extension of gravity, proposed in Ref. [51] and studied in the recent years. Furthermore, by considering g(P ) = 0, we obtain a particular extension of the Einstein-Hilbert action which encodes specific effects due to the curvature, the modified f (R) theory of gravity [59,60]. This framework has been studied into the past years and represents a viable extension, one of the most studied directions in the scalar tensor theories of gravity [60][61][62][63].

III. THE POWER LAW DECOMPOSITION
In this section we shall analyze the physical features of the present cosmological model by adopting the dynamical system analysis, an important method specific to physical systems. We shall consider that thef (R, P ) functional is decomposed into a direct sum based on a power-law behavior, i.e.f (R, P ) = f 0 R n + g 0 P m . For this specific model we shall assume that f 0 , g 0 and n, m are constant parameters. In what follows we shall introduce the following auxiliary dimension-less variables: At this point we note that the choice of the dimensionless variables is not unique. This specific choice is motivated by the form of the first modified Friedmann relation, the constraint equation. If we take into account that f (R) = f 0 R n with f 0 , n constant parameters, then we obtain an additional relation between x 2 and x 1 variable, x 2 = nx 1 . In this case the dynamical system becomes a six dimensional system with the associated variables [s, x 1 , x 3 , z, y 1 , y 2 ]. Here the first auxiliary variable s is associated to the matter component, acting as an effective density parameter influenced only by the Hubble parameter and the first variation of the scalar curvature coupling, embedded into the f (R) function. Taking into account the existence conditions specific to f (R) gravity theories, we have: s ≥ 0. The second auxiliary variable which is independent, the x 1 term, is associated with the specific form of the scalar curvature coupling, embedded into the f (R) function. Next, the z variable is connected to the specific value of the scalar curvature R, balanced by the square of the Hubble parameter. Finally, the y 1 variable is associated to the effects due to the cubic component, while in y 2 we notice the influence from the variation of the P invariant term.
We can define the effective equation of state for the present cosmological system, The effective or total equation of state is connected to the value of the z variable sincė This relation appears due to the value of the scalar curvature, valid for the present metric. Taking into account the previous definitions for the auxiliary variables we can rewrite the Friedmann constraint equations as the following expression, reducing the dimensionality of the dynamical system by determining the s variable, The next step for the dynamical system analysis method assumes the transformation from the cosmic time t to N , where N = log(a). In this case we shall obtain an autonomous system of ordinary differential equations described in the following relations: In the dynamical system determined by the equations (30)- (34) we notice that we have two components which have to be determined,R H 4 andP H 8 . These components are determined by considering the acceleration equation (15) with the specific expressions for the scalar curvature and the P invariant term, expressed into the relations (10) and (8). Hence, the acceleration equation can be rewritten in terms of the auxiliary variables in the following way: Another relation betweenR andP is determined by considering the specific expressions for these invariants in the case of the Robertson-Walker metric. Taking into account the previous mentioned considerations, we obtain: At this step we can note that the dynamical system of differential equations (30)-(34) is completely autonomous and closed, ready for the dynamical system analysis. For the present model we have obtained two critical points, determined by setting the right hand side of the autonomous equations (30)-(34) to zero.
The first cosmological solution represents a critical line located at the following coordinates: A = x 1 , x 3 = 0, z = 12, For this critical line we note that the n, m parameters are influencing the location in the phase space structure. This critical line is characterized by an indefinite value In this case the value of the s variable is zero, describing a possible late time stage of the Universe. Hence, this solution can be regarded as a geometrical de-Sitter epoch, with the physical effects encoded into the specific form of the f (R) and g(P ) functions which are correcting the Einstein-Hilbert action. From a dynamical perspective this solution is always saddle, due to the specific form of the obtained eigenvalues: with the following expressions: For this solution we have displayed in Figs. 1-2 the specific real values of the fourth and fifth eigenvalues, by setting the value of the x 1 variable. At this point, we note that the n and m parameters are only influencing the spiral behavior in the phase space. In this case a spiral trajectory is obtained if the corresponding eigenvalues contain imaginary values. For the A solution a spiral behavior is obtained if the C 2 component represents a complex number.
The second cosmological solution is located at the following coordinates: , . (43) The effective equations of state presents a sensitivity to the values of the n and m parameters, In Fig. 3 we have displayed the graphical representation of the value of the effective equation of state as a function of the n and m parameters. These parameters are encoding the effects due to the curvature and the cubic invariant term, respectively. At this critical point the value of the s variable is the following: s = 8m −3m 2 8n 2 − 7n − 6 + mn 8n 2 + 5n − 30 − 4(n − 2)n 2 3n (m 3 (64n 2 − 112n + 25) + 2m 2 (76n − 5) − 16mn(7n + 2) + 32n 2 ) .
The values of the s variable for the second cosmological solution B is displayed in Figs. 4-5. We note that in this figure we have displayed only the [0, 1] interval. As previously stated, the s component is associated to an effective matter density variable which has to be positive from a physical point of view. Hence, we have displayed in Figs. 4-5 some possible non-exclusive intervals where the B solution is physically viable.
Furthermore, for the dynamical analysis we have obtained the following eigenvalues:

IV. THE EXPONENTIAL DECOMPOSITION
In this section we shall adopt a second parameterization which involves a specific exponential decomposition wheref (R, P ) = f (R) + g(P ) = f 0 e nR + g 0 e mP , with f 0 , g 0 , n, m constant parameters. In this case we shall introduce the following auxiliary variables by analyzing the Friedmann constraint equation: If we take into account the exponential decompositioñ f (R, P ) = f 0 e nR + g 0 e mP the x 2 variable becomes a dependent component, while the remaining independent variables are the following: [s, the Friedmann constraint equation (14) has the following form: reducing the dimension of the corresponding dynamical system with one unit. If we take into account a pressureless matter fluid (w m = 0), then the second modified Friedmann relation can be written as: In this case the final dynamical system obtained is described by the following differential equations: (62) In order to close the dynamical system an additional relation betweenR andP is obtained, by differentiating the values specific to R and P components with respect to time. Considering the auxiliary variables, we obtain the following relation: At this point the second system of ordinary differential equations (57)-(62) is completely autonomous and can be analyzed using a dynamical system approach. In this case we have identified two critical points by analysing the right hand side of the (57)-(62) equations.
The first class of critical points is located at the following coordinates: , x 3 = 0, z = 12, y 2 = 0 .
From a dynamical point of view the D solution is always saddle with one positive eigenvalue and one negative eigenvalue. We also note that the first two eigenvalues are equal to zero, describing a non-hyperbolic solution. Hence, the values of the n and m coefficients are affecting the behavior in the phase space for the spiral properties of the corresponding trajectory. A non-exclusive region where the trajectory in the phase space is spiral due to the complex behavior of the last two eigenvalues in eq. (65) is displayed in Fig. 9.
Lastly, the second class of critical points is located at the following coordinates: The cosmological epoch associated to this class represents also a de-Sitter era (w ef f = −1) where the geometrical dark energy component completely dominates in terms of density parameters (s = 0). For this solution we note that the auxiliary variables x 1 is independent. From a physical point of view, the F solution represents a geometrical de-Sitter epoch where the values of n, m and β parameters are affecting the location in the phase space structure through the value of the y 1 component. The eigenvalues of this solution are the following: For this solution the specific form of the last two eigenvalues Ξ 5 , Ξ 6 is not displayed due to the complicated form of the obtained relations. From a dynamical point of view, the E solution has a similar behavior to the previous one, describing a non-hyperbolic epoch always saddle where the de-Sitter dynamics appears from the curvature and cubic extensions. As in the previous point, the values of the n and m coefficients are affecting only the spiral properties of the phase space. In the Fig. 10 we have displayed a possible region characterized by a spiral behavior due to the imaginary values of the last eigenvalues Ξ 5 and Ξ 6 .

V. CONCLUSIONS
In this paper we have extended the Einstein-Hilbert action by adding a generic termf (R, P )) which depends on the scalar curvature R and the specific invariant P , based on different cubic contractions of the Riemann tensor. This action can be regarded as an attempt which corrects the ΛCDM model by encoding particular geometrical invariants, leading to possible interesting dynamical effects. In our approach we have first considered that the generic termf (R, P ) can be written as a direct sum between the geometrical constituents, taking into account the power law decomposition, i.e.f (R, P ) → f 0 R n + g 0 P m .
After we have obtained the specific modified Friedmann equations, the constraint and the acceleration equation, we have investigated the physical characteristics of the present cosmological model by considering an analytical approach based on the dynamical system analysis. Hence, we have introduced specific dimensionless variables, approximating the evolution of the cosmological model as an autonomous system of ordinary differential equations, applying the analytical methods associated to dynamical systems. In this regard, we have obtained the critical points specific to the present cosmological model, analyzing the location in the physical phase space and the viability of these solutions. For each critical point we have determined the associated eigenvalues, detecting the dynamical characteristics. In the case of the power law decomposition the present dynamical system has two critical points.
The first critical point represents a de-Sitter epoch, where the effective equation of state corresponds to a cosmological constant-like solution. Analyzing the eigenvalues corresponding to this critical point, we have noticed the this point cannot be stable or unstable, always representing a saddle solution, due to the existence of one eigenvalue with positive real part, and one eigenvalue with a negative real part. Moreover, due to the presence of one zero eigenvalue, we also note that this solution represents a non-hyperbolic equilibrium in the physical phase space.
The second critical point represents a cosmological era where the dynamical features depend on the specific values of the geometric couplings, the n and m parameters, encoding effects due to the curvature and the cubic gravity type parameterization. In this case the effective equation of state is sensitive to the values of the n and m parameters and can describe various epochs. Hence this solution can also be associated to a matter or radiation era for specific values of the n and m constant parameters. Moreover, it can describe also a quintessence or phantom-like epoch, explaining the super-acceleration scenario. Furthermore, it can also be associated to some cosmological solutions, stiff or super-stiff dynamics. For this solution the form of the eigenvalues depend on the specific values of the n and m parameters, describing a hyperbolic equilibrium. We have analyzed the second equilibrium point from a dynamical point of view, showing that in certain cases this solution can be stable, saddle, or unstable. Lastly, some specific solutions have been analyzed, determining the evolution and the viability of the analytical solutions.
In our study we have also considered an exponential decompositionf (R, P ) → f 0 e nR + g 0 e mP , with f 0 , g 0 , n, m constant parameters in Sec. IV. For this specific case we have analyzed the structure and properties of the corresponding phase space, revealing the cosmological solutions obtained. For the exponential model we have identified two cosmological de-Sitter epochs where the values of the n and m parameters are affecting the location in the phase space and the dynamical properties of the corresponding trajectories.
The present paper offers a generalization to the f (R) gravity theories, by embedding an invariant component denoted g(P ), based on specific contractions of the Riemann tensor in the third order. In this way the f (R) gravity theory is extended, by including specific geometrical manifestations from the cubic component, offering a generalization to the basic Einstein-Hilbert action to-wards a more complete gravity theory. In the current manuscript the investigation is based on the usage of dynamical system analysis, an important tool in cosmology [59]. We mention here that in order to better discriminate between the effects due to the scalar curvature part and the cubic component, various reconstruction methods can in principle be applied [11] to the present model, obtaining possible constraints due to different dynamical behaviors. Another important aspect in the dynamical system analysis is related to the viability of the corresponding singularities at finite time. For the f (R) gravity theory the properties of the resulting singularities have been investigated in Ref. [63]. In principle, such an analysis can be applied for specific models of f (R, P ), obtaining possible constraints from the physical properties of the corresponding models.
We note that the present cosmological model can be also extended in various applications, by analyzing different specific models and solutions to the gravitational field equations, by considering effects due to the curvature and the cubic term. One particular extension is related to the inclusion of the Gauss-Bonnet topological invariant, leading to a more generic theory of gravity. A different approach is related to the study of the inflationary era, an epoch where this theory can have visible effects in the early times. Furthermore, one can consider different specific models for thef (R, P ) gravity by taking into account possible parameterizations, as an ex-ponential decomposition or more complex models which can include a cosmological constant. All of the previous mentioned questions and directions are open and left as possible future investigations.
Another important aspect in modified gravity theories is related to the study of various instabilities which can manifest. In the case of Einsteinian cubic gravity various studies [49,50] have indicated that some pathological instabilities might emerge, leaving the corresponding theory unhealthy from a physical point of view. The gravitational theory studied at the background level in the present paper represents a generalised attempt based on non-linear cubic and curvature extensions for the Einstein-Hilbert action. To this regard, the current model have to be considered only as an effective theory and needs to be analyzed from this point of view. An expanded discussion on this issue is presented in a recent publication [57].

VI. ACKNOWLEDGEMENTS
For this project we have considered different analytical computations in W olf ram M athematica [64] and xAct [65]. The computational part of this work was performed utilizing the computer stations provided by CN-FIS through the project CNFIS-FDI-2020-0355.