Cosmological evolution of two-scalar fields cosmology in the Jordan frame

In the present article we study the cosmological evolution of a two-scalar field gravitational theory defined in the Jordan frame. Specifically, we assume one of the scalar fields to be minimally coupled to gravity, while the second field which is the Brans-Dicke scalar field is nonminimally coupled to gravity and also coupled to the other scalar field. In the Einstein frame this theory reduces to a two-scalar field theory where the two fields can interact only in the potential term, which means that the quintom theory is recovered. The cosmological evolution is studied by analyzing the equilibrium points of the field equations in the Jordan frame. We find that the theory can describe the cosmological evolution in large scales, while inflationary solutions are also provided.


Introduction
Observational data since last 20 years [1][2][3][4][5][6][7] have consistently reported that our universe is currently passing through a phase of accelerating expansion. This has been one of the mysterious discoveries of modern cosmology that thrilled the entire scientific community with a big question mark. Within the paradigm of standard cosmology, it is impossible to fit such accelerating phase since it demands something more in terms of some exotic component having negative pressure enabling the expansion of the universe in an accelerating manner. Usually two distinct but well known approaches are used to introduce such exotic components. Both the approaches essentially need the introduction of new terms in terms of the new degrees of freedom into the Eina e-mail: alexgiacomini@uach.cl b e-mail: genly.leon@ucn.cl c e-mail: anpaliat@phys.uoa.gr (corresponding author) d e-mail: supriya.maths@presiuniv.ac.in stein's field equations. The simplest approach to introduce such new terms is the modifications of the matter distribution of the universe within the context of Einstein gravity, and these new exotic components are usually dubbed as dark energy (DE) fluids [8][9][10][11][12][13][14][15][16][17]. Alternatively, one could introduce such new additional terms through the modifications of the Einstein-Hilbert action [18][19][20][21][22][23][24][25][26][27][28][29]. The resulting exotic terms coming from the gravity modifications are usually designated as the geometrical dark energy (GDE) fluids. According to latest observational estimations [7], the percentage of such exotic dark energy or geometrical dark energy fluids whatever the universe has within it, is around 68% of its total energy budget. Hence, the investigations towards this direction is an essential part of modern cosmological research and we still are looking for an actual cosmological theory fitting the observational results.
The theory of scalar fields have played a significant role in the two different approaches proposed by the cosmologists. When scalar fields are introduced in the gravitational action integral, they provide new degrees of freedom which drives the dynamics of the field equations in a way to explain the observable phenomena and the cosmological history [30]. Scalar fields can be minimally or nonminimally coupled to gravity. When the scalar fields are introduced in such a way where they are minimally coupled to gravity, then we are in Einstein's General Relativity, where the contribution of the scalar fields is on the energy momentum tensor. The most common known minimally coupled scalar fields are, for instance, the quintessence, phantom, quintom and Chiral cosmological models [8,[31][32][33][34][35][36][37][38][39][40].
Alternatively, scalar fields can be introduced in the gravitational action such that they nonminimally couple to gravity. In these models, the contribution of the scalar field in the field equations includes the terms which correspond to the energy momentum tensor, and in addition, dynamical terms of geo-metric origin. In these models, usually, the gravitational constant k, is not a constant anywhere, rather it becomes variable. A nonminimally coupled scalar field was introduced by Carl H. Brans and Robert H. Dicke to provide a gravitational theory which satisfies the Machian Principle [20]. Brans-Dicke theory is a special case of the so called scalartensor theory [41]. A main characteristic of the Brans-Dicke field is the parameter ω known as the Brans-Dicke parameter. Indeed, for small values of ω the contribution of the scalar field in the dynamics of the field equations is significant, while on the other hand, when ω is large, the main contribution in the dynamics is followed by the tensor part. Someone would expect that as ω reaches infinity, the theory of General Relativity will be recovered, however, that is not true, which means that General Relativity and Brans-Dicke gravity are fundamentally different theories [42]. Other wellknown nonminimally coupled models are the Galileon and the Hordenski theories [43][44][45][46][47][48].
The nonminimally coupled scalar fields can describe modified theories of gravity by attributing new degrees of freedom in terms of geometric invariants which are introduced to modify the Einstein-Hilbert gravitational action. Indeed, the fourth-order f (R)-gravity is dynamical and physically equivalent to the Brans-Dicke theory in the limit where ω is zero, such that Brans-Dicke theory reduces to O'Hanlon gravity [49]. For more details in this direction we refer the reader to [50][51][52][53][54] and the references therein.
As we discussed, Brans-Dicke theory is a different theory compared to Einstein's General Relativity, however, there is a mathematical connection between these two theories. Indeed, the two theories are conformally equivalent [41,55,56]. The latter equivalency is a mathematical equivalence between the Jordan and the Einstein frame in the solution space of the field equations while in general the physical properties of the solutions do not remain invariant under the conformal transformation. For instance, a singular universe in the Einstein frame, under a conformal transformation can become a non-singular universe in the Jordan frame which is described by a scalar tensor theory. There is a plethora of studies where the physical quantities have been studied between the Einstein and the Jordan frame for various exact solutions [56][57][58][59]. In addition, the question of which is the preferred frame is still unanswered [60].
In this work we considered a two-scalar field cosmological model where one of the fields is the Brans-Dicke field which is nonminimally coupled to gravity, and the second scalar field is minimally coupled to gravity. The model of our consideration has the Action Integral of the Brans-Dicke Action where the Brans-Dicke field is coupled to the Lagrangian of the quintessence scalar field. Consequently, the two fields interact in the kinetic and in the potential parts. Such a model was investigated before in [61]. Specifically, in [61] the integrability of the cosmological model was investigated and some exact solutions were determined. While our model is defined in the Jordan frame, when we perform a conformal transformation in order to pass to the Einstein frame, we end up with a two-scalar field theory which includes the quintom cosmological model.
Here we study the cosmological evolution of this specific cosmological model in the Jordan frame by investigating the stability of the equilibrium points. For the background space we assume this to be described by the spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) universe. Under this minimal assumption we shall be able to understand the physical evolution of the universe as it is provided by the field equations and also we will be able to discuss the viability of the cosmological model. Such an analysis has been applied in various cosmological theories with many interesting results [62][63][64][65][66][67][68][69][70][71][72][73][74][75][76]. Any equilibrium point, also known as stationary or critical point, describes a specific phase, i.e., an exact solution, in the evolution of the generic solution. The stability of the critical point indicates the stability of the solution, from where we can draw conclusions for the cosmological evolution.
The plan of the paper is follows. In Sect. 2, we present the cosmological model of our consideration. We derive the field equations and we discuss the mathematical relation of the theory between the Jordan and the Einstein frame, where we show that our model is equivalent with a twoscalar field theory, such is the quintom cosmology, in the Einstein frame. In Sect. 3, we write the field equations in dimensionless variables by using the H -normalization. We end with a five-dimensional algebraic-differential dynamical system. Section 4, includes the main analysis of our work, where we investigate the critical points and their stabilities for the model under consideration in the vacuum case, and in presence of an additional dust fluid source. The scalar field potentials have been selected in order to reduce the number of equations of the dimensionless system under study. However, such selection can also be seen as the limit for other potential forms too. Finally, in Sect. 5, we summarize our results with the main findings.

Cosmological model
The gravitation model of our consideration is described by the following Action Integral which actually describes a two-scalar field cosmology defined in the Jordan frame. More specifically, the scalar field φ (x μ ) is the Brans-Dicke field with ω = 2 3 , which is nonminimally coupled to gravity and to the second scalar field ψ (x μ ) as well, while the second scalar field ψ (x μ ) is minimally coupled to gravity. The limit ω = 2 3 corresponds to a special case where the Brans-Dicke field does not introduce any real degree of freedom and can be neglected by a scale transformation in the metric g μν . Therefore, in the following we shall assume that ω = 2 3 . An equivalent way to write the Action Integral (1) is by using the Lagrange function for the quintessence, that means, where L Q R, ψ, ψ ;μ is the Lagrangian function for the minimally coupled scalar field theory, that is, and . Consequently, when φ (x μ ) = φ 0 , the gravitational Action Integral (1) reduces to that of a minimally coupled scalar field model. From (1), one can obtain the gravitational field equations by varying this action with respect to the metric tensor g μν as follows [61] while the variation of (1) with respect to φ (x μ ), and ψ (x μ ) respectively provide with the equations of motion of the two fields as follows, It is important to mention here that while someone will expect that for φ ;κ = 0, the limit of General Relativity with a minimally coupled scalar field will be recovered, but that is not true, because of the constraint equation (5) which becomes R − 2V ,φ 0 = 0. In order to find an equivalent theory to (1) in the Einstein frame, we consider the conformal transformation g μν = φḡ μν , that is, for a four-dimensional manifold, the Ricci scalar R becomes [77] Replacing (7) and integrating by pats we end up with the gravitational Action Integral where the new scalar field, (x μ ), is a functional form of φ (x μ ) , that is, = (φ (x μ )), and specifically, ( in the Einstein frame, the two scalar fields interact. Therefore, when U (φ, ψ) = V (ψ), the Action Integral describes a two-scalar field minimally coupled cosmological model, and such model has been studied earlier in the literature [78].
Moreover, if in the Jordan frame, we select V (φ, ψ) = φ 2 F (ψ) + Z (φ), then in the Einstein frame it follows that the two scalar fields are not interacting neither in the kinetic terms nor in the potential terms. In addition, in the latter scenario, by doing the change ψ (x μ ) = iψ (x μ ), the Action integral (8) takes the form of the quintom theory [68,79].
For the background geometry of our universe we assume that in the large scales it is well described by the spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) spacetime characterized by the following line element where N (t) is the lapse function and a (t) is the expansion scale factor of the universe, i.e. the Hubble function of this universe is defined as H (t) = 1 Nȧ a . In the following we assume the co-moving observer u μ = 1 N δ μ t , u μ u μ = −1. For the line element (9), the gravitational field equations (4), (5) and (6) can be expressed as [61] 6φ 2φḢ where without any loss of generality we selected N (t) = 1, while additionally, we assumed that the two scalar fields φ (x μ ) , ψ (x ν ) inherit the symmetries of the spacetime, that means, φ (x μ ) = φ (t) and ψ (x μ ) = ψ (t).
Let us note that Eq. (10) is the constraint equation, while Eqs. (11), (12) and (13) are the second-order differential equations describing the dynamics of the universe for the two scalar fields φ (t) , ψ (t) and the scale factor a (t). Now, in presence of an additional matter source, and more specifically, in the presence of a dust fluid term minimally coupled to the Brans-Dicke field φ (x μ ) in the gravitational model, the field equations (11), (12), (13) will remain same, except the constraint equation (10) which becomes while the continuity equation for the dust fluid isρ m + 3Hρ m = 0, from where it follows that, ρ m = ρ m0 a −3 , where ρ m0 is an integration constant physically which represents the energy density of the dust fluid at current epoch. We investigate first the vacuum case ( m = 0), and after that we investigate the matter case. The gravitational field equations for this particular cosmological model are described by a point-like Lagrangian, however, such an analysis is not necessary in the present work. The purpose of this work is to study the dynamical evolution for this cosmological model. In [61], the specific forms of the scalar field potential were determined such that the field equations admit linear and quadratic conservation laws in the momentum and they are integrable. In addition, exact solutions of special interests, such that the singular solution a (t) = a 0 t p and the de Sitter solution a (t) = a 0 e H 0 t , were determined in [61].
With the analysis of this work, we will be able to understand the general behavior of the cosmological solution as well as we will be also able to identify the attractors of the cosmological model and to study the stability of the particular solutions. At the same time, we will be able to investigate the similarities and differences of two-scalar field gravitational theories defined in the Einstein and Jordan frames.
Making the selection of the potential for the two scalar field V (φ, ψ) = V (ψ) W (φ), we solve for the higher derivatives (in presence of an additional matter source and more specifically in presence of a dust fluid which is minimally coupled to the Brans-Dicke field) and obtain the field equations: with first integral equation (14).

Dimensionless system
In this section we proceed by writing the field equations with the use of dimensionless variables. In particular, we make use of the H −normalization [80], which has been applied also in Brans-Dicke theory [81][82][83][84] and in multi scalar field theories [85][86][87][88][89][90]. As far as the scalar field potential is concerned, we make the selection The new dimensionless variables are or equivalently, In the new variables, the field equations can be expressed as the following algebraic-differential first-order system in which the new variables λ and μ are defined as functions of φ while the constraint equation (14) becomes where m is the dimensionless density parameter of the dust fluid source defined as m = ρ m 3φ H 2 , in which m ∈ [0, 1]. We also have an extra evolution equation As far as the variation of the new variables μ and λ is concerned, we find that they should satisfy the following two first-order ordinary differential equations, namely, in which Consequently, the algebraic-differential system has dimension five, consisting of the first-order differential equations (21), (22), (23), (28) and (27), while for specific forms of the potentials V (ψ) and W (φ), it follows that, λ = const and μ = const, which mean that the dimension of the dynamical system can be reduced by one or by two dimensions. In addition, for the vacuum case, i.e., m = 0, the algebraic equation (25) can also be used to reduce the dimension of the dynamical system by one dimension.
Any critical point describes a specific universe where the total equation of state is w tot = −1 − 2 3Ḣ H 2 , which can be be expressed in terms of the dimensionless parameters as When w tot = −1, the scale factor at the critical point is a power-law function, while for w tot = −1, the critical point describes a de Sitter universe.

Cosmological evolution
In the following, we set V (ψ) = V 0 e λφ and W (φ) = W 0 φ μ , and we study the equilibrium points and analyze the stability when m = 0 (vacuum case) and when 0 < m ≤ 1 (matter case). For the specific functions of the two potentials, variables μ and λ are constants, i.e. the right hand side of equations (27), (28) are zero.

Vacuum case
In the vacuum case, namely for m = 0, the algebraic equation (25) gives an additional restriction that can be used to remove one phase-space variable, say z. Hence, we obtain the reduced system: The total equation of state parameter now becomes, For the choices V (ψ) = V 0 e λφ and W (φ) = W 0 φ μ , with μ and λ being constants, we study the Equilibrium points and investigate the stability of the two-dimensional system (32), (33).
The (lines of) equilibrium points are the following: The line corresponding to the massless fields, exists for The eigenvalues are: Fig.   1 we show the stability conditions of A + for (a) ω = 1, and for (b) ω = 0. For parameters in the red region (upper and lower left graphs of Fig. 1) A + is a source and for parameters in the blue region (upper and lower right graphs of Fig. 1) it is a sink. Similarly, in Fig. 2 we present the stability conditions of A − for (a) ω = 1, and for (b) ω = 0. For parameters in the red region (upper and lower left graphs of Fig. 2) A − is a source and for parameters in the blue region (upper and lower right graphs of Fig. 2) it is a sink. The total EoS parameter and the fractional matter density at the equilibrium points A ± are the following: . For In particular, it mimics a dust solution (w tot (A ± ) = 0) for . Eigenvalues are: . The critical point is a stable node for .
The critical point an unstable node for Otherwise the critical point is nonhyperbolic. The total EoS parameter and the fractional matter density at the equilibrium point B are the following: . Therefore, they represent solutions dominated by the scalar field and they are accelerating (w tot (B) < − 1 3 ) for .
Eigenvalues are: Nonhyperbolic critical point with a 1D stable manifold for: , or Nonhyperbolic critical point with a 1D unstable manifold for: , or Eigenvalues are: Nonhyperbolic critical point with a 1D stable manifold for: Nonhyperbolic with a 1D unstable manifold for The total EoS parameter and the fractional matter density at the equilibrium point C ± are the following: . They represent solutions dominated by the scalar field and they are accelerating (i.e., 3 2 < ω < μ+1 2 .
In particular, for ω = μ 6 + 7 6 , it represents a de Sitter solution. It is is decelerating for For ω = 11 6 − μ 6 it mimics a dust fluid.

Center manifold calculations
Now, we proceed to the center manifold calculation of C + , and C − .
a. The center manifold of C + : For the calculation of the center manifold of C + , we introduce the new variables with inverse Then, we obtain the dynamical system Hence, the center manifold is given by the graph (u, v) = (u, h(u)), with h(0) = 0, h (0) = 0, and satisfying the equation Using the Taylor series expansion we propose the expansion Substituting this expansion into equation (39), and equating the coefficients of the same powers of u, we obtain Fig. 3 Qualitative evolution in the phase plane (u, v) given by (38) for ω = 1. The solid red lines denote the implicit function x(u) x(u)ω + √ 6 + y(u, v) 2 + 1 = 0, with x(u) and y(u, v) defined by (37) for ω = 1, whereas the dashed black lines represent the approximated center manifold of C + . Observe the accuracy for small u a = The evolution equation on the center manifold is reduced to u = O(u 6 ), which means that u is approximately constant on the center manifold. In Fig. 3 we show the qualitative evolution in the phase plane (u, v) given by (38) for ω = 1. The solid red line in each plot denotes the implicit function x(u) x(u) + √ 6 + y(u, v) 2 + 1 = 0, with x(u) and y(u, v) defined by (63) for ω = 1, whereas the dashed black lines in each plot represent the approximated center manifold of C + . In Fig. 4 we show the qualitative evolution (u, v) given by (38) for ω = 0. The solid red lines in each plot denote the implicit function √ 6x(u) + y(u, v) 2 + 1 = 0, with x(u) and y(u, v) defined by (37) for ω = 0, whereas the dashed black lines shown in each plot represent the approximated center manifold of C + . Observe the accuracy for small u. Fig. 4 Qualitative evolution in the phase plane (u, v) given by (38) for ω = 0. The solid red lines denote the implicit function x(u) x(u)ω + √ 6 + y(u, v) 2 + 1 = 0, with x(u) and y(u, v) defined by (37) for ω = 0, whereas the dashed black lines represent the approximated center manifold of C + . Observe the accuracy for small u The numerical solutions suggest that the exact center manifold of C + is the line x xω + √ 6 + y 2 + 1 = 0, where x, y are given by the inverse relations of (37).
Indeed, the two branches of this implicit equations are Both are solutions of the equation (39). But the only one that satisfies the tangential conditions at the origin, say h(0) = 0, h (0) = 0, is (42b). The exact equation that dictates the dynamics at the center manifold is u = 0. It means that u is constant at the center manifold.
b. The center manifold of C − : For the calculation of the center manifold of C − , we introduce the new variables with inverse Then, we obtain the dynamical system Hence, the center manifold is given by the graph (u, v) = (u, h(u)), with h(0) = 0, h (0) = 0, and satisfying the equation Using the Taylor series we propose the expansion Substituting this expansion in the equation (46), and equating the coefficients of the same powers of u, to obtain The evolution equation on the center manifold is reduced to In Fig. 5, we show the qualitative evolution in the phase plane (u, v) given by (45) for ω = 1. The solid red lines denote the implicit function x(u) x(u) + √ 6 + y(u, v) 2 + 1 = 0, with x(u) and y(u, v) defined by (44) for ω = 1. The dashed black lines in Fig. 5 represent the approximated center manifold of C − . In Fig. 6, we present the qualitative evolution of (u, v) given by (45) for ω = 0. The solid red lines denote the implicit function √ 6x(u) + y(u, v) 2 + 1 = 0, with x(u) and y(u, v) defined by (44) for ω = 0, whereas the dashed black lines represent the approximated center manifold of C − . Let us observe the accuracy for small u.  (u, v) given by (45) for ω = 1. The solid red line denotes the implicit function x(u) x(u)ω + √ 6 + y(u, v) 2 + 1 = 0, with x(u) and y(u, v) defined by (44) for ω = 1, whereas the dashed black lines represents the approximated center manifold of C + . Observe the accuracy for small u There are two branches of C: Fig. 6 Qualitative evolution in the phase plane (u, v) given by (45) for ω = 0. The solid red line denotes the implicit function x(u) x(u)ω + √ 6 + y(u, v) 2 + 1 = 0, with x(u) and y(u, v) defined by (44) for ω = 0, whereas the dashed black lines represents the approximated center manifold of C + . Let us observe the accuracy for small u Both are solutions of the equation (46). But the only one that satisfies the tangential conditions at the origin, say h(0) = 0, h (0) = 0, is (50b). In this case the exact equation that dictates the dynamics over the center manifold of C − is u = 0. It means that u is constant at the center manifold.

Matter case
In this section we determine the equilibrium points of the dimensionless dynamical system (21), (22), (23). For each critical point P we calculate the physical quantities of the solution, which are the equation of state parameter for the effective fluid, i.e. w tot (P), the energy density of the dust fluid source, m (P). As far as concerns the quantity = x xω + √ 6 + y 2 + z, this describes the effective energy density for the two scalar fields. Recall that there is an interaction between the two scalar fields and the interaction in the dynamical term is included in the variable z. Hence, there is no unique definition of φ and ψ .
In order to avoid any selection of φ and ψ which can lead to false conclusion of the nature of the two scalar fields at the equilibrium points, we choose to see the two-scalar fields as an effective fluid. The only cases where a conclusion for the nature of the scalar field can be done is, if at the critical point P, one of the coordinates is zero, that means, x (P) = 0 , or y (P) = 0, or z (P) = 0. When x (P) = 0, by definition φ = const which means that the Brans-Dicke field does not contribute to the physical state of the solution. Moreover, when y (P) = 0, it follows ψ = const , where now the scalar field ψ does not contribute in the evolution of the universe. Finally, when z (P) = 0, then only the kinetic parts of the two scalar fields, namely, φ and ψ contribute in the cosmological evolution at the specific point.
Since m ≥ 0, therefore, from the restriction we have that the phase space is given by In this example we study the equilibrium points and their stability of the dimensional system (21), (22), (23).
As far as the equation of state is concerned for the effective fluid derived to be w tot (A 3 ) = − 1 μ − λ 2 3μ(μ−1) . From the latter it follows that −1 ≤ w tot (A 3 ) < 1 when In Fig. 8 we present the contour-plots for the effective equation of state parameters w tot (A 2 ) and w tot (A 3 ) in the region space of the variables {λ, μ}.
We continue our analysis by studying the stability conditions of the critical points. As we shall see, one can easily conclude about the stability of the equilibrium points A 2 , A 3 , while for points A 1± the center manifold theorem should be applied.
The eigenvalues of the linearized system around point A 2 are found to be from where we infer that the point is stable when As far as the critical point A 3 is concerned, the eigenvalues of the linearized system are found to be where Because of the nonlinear terms one can solve the stability conditions Re (e 1 ) < 0, Re (e 2± ) < 0 numerically. The regions of values {λ, μ} where the point is an attractor are presented in Fig. 9. In addition in Fig. 9 we show the regions in the space {λ, μ} where point A 2 is stable. For point A 3 we show also the regions where it is stable and satisfies the condition 0 ≤ m ≤ 1 (physically acceptable). Combining both conditions we obtain that A 3 "exists" (0 ≤ m ≤ 1) and it is stable for a. Center manifold theorem for A 1± : The parametric equation for A 1± can be derived from from which we can define The linearization matrix evaluated at A 1± is then with eigenvalues 0, √ 3 cos(θ ), √ 3λ sin(θ ) Due to the dynamics on the sets A 1± , which is essentially two-dimensional, to study the dynamics over A 1± , we project the flow on the plane (θ, z), where the dynamics is given by For the calculation of the center manifold of A 1+ we introduce the new variables with inverse where we have used θ = θ c to parametrize the system at an specific point of A 1± with coordinates Then, we obtain the dynamical system where In the coordinates u, v, the eigensystem of the linearization is -as expected - Hence, the center manifold is given by the graph (u, v) = (u, h(u)), with h(0) = 0, h (0) = 0, and satisfying the equation For which we have the trivial solution. Now, imposing the compatibility conditions we can seek for another, non-trivial solution using Taylor series with results on a = 1 4 . In the two cases the evolution equation on the center manifold is reduced to that means that u is constant at the invariant manifold. In Figs. 10 and 11 we present some orbits of the phase portraits of (65) for different values of λ, μ ad θ c .

Equilibrium points for ω = 0
The case ω = 0 is of special interest because the Brans-Dicke field reduces to that of O'Hanlon theory [49]. The latter theory is dynamically equivalent with the fourth-order theory known as f (R)-gravity [50]. The dynamically equivalent scenario can be found easily with the use of a Lagrange multiplier, for more details we refer the reader to [50] and references therein.
For ω = 0, we evaluated the cosmological observables and discuss the stability conditions of the equilibrium points of the dimensionless dynamical system (21), (22), (23), as follows.
Nonhyperbolic with a 2D unstable manifold for λ ∈ R, − √ 2 < y c < √ 2, μ < −y 2 c + √ 6λy c +5 y 2 c +1 . The cosmological observables are m (B 1 ) = 0, and w tot (B 1 ) = 1 3 1 − 2y 2 c . That means the solution corresponds to a scalar field dominated universe, that is accel- Fig. 11 Phase portrait of (65) for different values of λ, μ ad θ c erating for y 2 c > 1. In particular, it is a de Sitter solution for y c = ± √ 2. It is decelerating for y 2 c < 1. In particular, it mimics a dust solution for y c = ± In particular, it mimics a dust solution for μ = 1 3 3−λ 2 , λ λ 2 − 12 = 0. 6x c , similar expression with that of A 1± . In the case where √ 6x c = −1 where only the field φ contributes to the evolution of the universe, we find that w tot (B 1± ) | √ 6x c →−1 = 1 3 , where the critical point in which the f (R)-gravity describes the radiation era is recovered.

Center manifold of B 1
It is more convenient to use the parametrization B 1 : (x, y, z) = − y 2 c +1 √ 6 , y c , 0 , which is solved globally for x, rather than considering two branches B 1± , which provide with only a local description of the model. We choose y 2 c = 2, in which the center manifold is one dimensional. Next, we introduce the new variables with inverse z = w. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .