{'fraction_non_alphanumeric': 0.03560966548063046, 'fraction_numerical': 0.011917520469841814, 'mean_word_length': 4.537426808039247, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
From inflation to dark energy in scalar-tensor cosmology

The methods of dynamical systems have found wide applications in cosmology, with focus either upon inflation or upon the passage into dark energy era. In this paper, we endeavor to capture the whole history of the universe into a dynamical system by considering generic scalar tensor gravity with radiation and dust matter fluids in flat Friedmann-Lemaitre-Robertson-Walker spacetime. We construct the dynamical variables in such a way that the main stages of the cosmic evolution, viz. inflation, radiation domination, matter domination, and dark energy domination can be represented by the respective fixed points in the phase space. As the evolution of solutions is ruled by the sequence of these fixed points with appropriate properties, we can determine the conditions that the scalar potential and nonminimal coupling must satisfy for the model to deliver viable cosmic history in a generic manner. We illustrate the construction by a scalar field with quartic potential, with and without quadratic nonminimal coupling to curvature.


I. INTRODUCTION
The methods of dynamical systems offer a set of helpful tools to analyze and visualize the qualitative behavior of solutions of ordinary differential equations. These methods have been applied with great success also in cosmology, see Refs. [1][2][3] for general reviews and e.g. Refs. [4][5][6][7][8][9][10][11] for a sample of very recent papers. In a typical approach, to scrutinize the evolution for a model given before, one writes the cosmological equations as a dynamical system, finds the fixed points and from these to deduce the generic features of cosmological evolution, as the fixed points shape the course of trajectories in the phase space. In cosmology, most authors so far have focused either upon early (inflation) or late (matter to dark energy domination) stages of the cosmic history. The present work ventures beyond the usual approach in two respects.
First, we intend to describe by a dynamical system the full history of the universe, starting from the onset of inflation till the rule of dark energy. Although the idea to unify the early and late stages of the universe within a single theory is not new [12][13][14][15][16][17][18][19][20][21][22][23], the papers that attempt to tackle the problem in the context of dynamical systems are surprisingly few [24][25][26]. Yet, mapping out the global phase space should help to address some basic questions in theoretical cosmology: how common with respect to all possible initial conditions is the occurrence of inflation (early period of accelerated expansion), how graceful is the exit from inflation into a radiation domination era, how common it is for the solutions to trace through all the eras (e.g. not to skip some of them, or get stuck in one of them), etc. In our paper, we do not claim to provide comprehensive answers to these questions, however, we offer a set of ideas how a viable cosmic history should be realized, which could serve as a basis for future detailed investigations.
Second, we slightly alter the usual mode of study which takes a model and then analyzes its phase space. Instead, we first postulate that in order to reproduce the correct history of the universe, the phase space must be endowed with fixed points corresponding to all the eras, and then investigate which models belonging to a wide family of theories are able to provide these points with appropriate characteristics. The aim is to determine the conditions on the model functions and parameters that allow to realize the full cosmic history in a generic manner for many solutions (initial conditions) within the model. It turns out the requirement most difficult to satisfy is to have the early universe dynamics ruled by a saddle point with de Sitter like expansion and with a repulsive eigendirection pointing towards the fixed point exhibiting radiation domination. The fixed points of radiation and dust matter domination are quite naturally saddles, possibly with a heteroclinic orbit between them, to guide the flow of trajectories towards the latter. Finally, there must exist an attractive de Sitter fixed point responsible for dark energy era, that collects all the available trajectories in the phase space. Scalar-tensor gravity (STG), a family of generalized Jordan-Fierz-Brans-Dicke theories extend general relativity (GR) by introducing a scalar field nonminimally coupled to curvature, which mediates the gravitational interaction along with the usual tensor degree of freedom [27][28][29]. In effect, the gravitational constant becomes dynamical depending on the value of the scalar field. While astrophysical observations put quite stringent limits on the variation of the gravitational constant [30][31][32], a large class of STG models naturally accommodates it by so-called attractor mechanism, whereby the cosmological evolution from the onset of matter domination era pulls the solutions to a regime, where the scalar field relaxes around a fixed value, the effective gravitational constant stabilizes, and the theory starts to behave rather like GR [33][34][35][36][37][38][39][40].
The dynamical system tools have been successfully applied in the context of STG cosmology with a special focus on concrete classes of coupling functions and potentials [43,52,54,[64][65][66][67][68], or considering arbitrary functions [69][70][71], and including a single matter component [24,66,[72][73][74][75][76]. General geometric properties of the STG phase space were discussed in Refs. [77][78][79]. In [80], a detailed analysis was performed for STG models with both radiation and matter, but focussing on late time unverse. An interesting complete history of the Universe from inflation till the dark energy era was obtained using the dynamical systems approach in Ref. [24], although here the radiation dominated epoch was generated by nonminimal coupling only. For an interested reader we may recommend Ref. [3] as a handy introduction to the dynamical systems tools and applications in cosmology together with a comprehensive literature review.
In the present work, we focus on the cosmological dynamics of the generic class of STG theories with a general non-minimal coupling function and self interacting potential. By using the methods of dynamical systems, we derive conditions on the coupling functions and potentials in which viable cosmological evolutionary history can be achieved by a large set of solutions. In Sec. II, we introduce the key features of STG and give the cosmological equations. Then in Sec. III, we argue for the choice of dynamical variables and build the dynamical system. The ensuing set of fixed points along with their properties is described in Sec. IV. Equipped with this general knowledge, we analyze the examples of GR with a cosmological constant in Sec. V, minimally coupled scalar field with a quartic potential and cosmological constant in Sec. VI, and nonminimally coupled scalar field with the same potential in Sec. VII. Finally, Sec. VIII presents the summary and outlook. In the appendix A, we apply center manifold theory for nonhyperbolic de Sitter points in minimally coupled models, which to our knowledge has not been done before in the literature.

II. SCALAR-TENSOR COSMOLOGY
A generic action functional for scalar field nonminimally coupled to curvature and endowed with a potential can be written as [61,62] Integrating over a 4-dimensional Lorentzian curved spacetime manifold V 4 , it incorporates four arbitrary functions of the dimensionless scalar field Φ: the coupling A(Φ) to the curvature scalar R, the generic kinetic coupling of the scalar field B(Φ), the self-interaction potential of the scalar field V(Φ), and the conformal coupling e 2α(Φ) between the metric g µν and matter fields χ, while we are not considering derivative couplings or higher derivative terms. Fixing the functions {A, B, V, α} gives a specific theory in certain frame and parametrization. Indeed, one of the functions can be fixed by a conformal transformation of the metric g µν = e 2γ(Φ)ḡ µν which sets the frame. Another function of the four can be fixed by scalar field redefinition Φ = f(Φ) [61,62]. We may use a conformal transformation to choose the Jordan frame, where α ≡ 0 and matter test particles follow geodesics. Second, we may redefine the scalar field to give the kinetic term a canonical form, B ≡ 1. Finally, without a loss of generality, we may separate the curvature coupling function into a constant part (representing a bare gravitational constant) and a field dependent part, A = 1 + f (φ). This leaves us with an action Here different possible choices of f (φ) and V (φ) yield different models in the scalar-tensor family, set in the Jordan frame and this particular parametrization. If the scalar field is nondynamical, the action reduces to that of GR where a constant nonzero value of the potential plays the role of the cosmological constant Λ. In the case f (φ) ≡ 0 we have a scalar field minimally coupled to gravity. The theory is free of ghosts provided since then in the Einstein frame, where the the tensor and scalar degrees of freedom are clearly separated, the scalar kinetic term is endowed with the correct sign [51,62]. Note that we have adopted the units with c = 1, but do not immediately fix the values of the nonvariable part of the effective gravitational "constant" κ 2 and keep a positive constant parameter which has the dimension of length, e.g. the Planck length. From a "gravity" convention [κ 2 ] = 1, it follows that [S] = [ ] = L 2 and from a "particle physics" convention [S] = [ ] = 1 it follows that [κ 2 ] = L 2 . In either case the scalar field φ is dimensionless.
The modified Einstein field equation obtained by varying the action (2) with respect to the metric g µν is given by where G µν is the Einstein tensor, T µν is the usual matter energy momentum tensor and 2 ≡ ∇ µ ∇ µ . We see that in STG, the Newtonian constant is replaced by an effective gravitational "constant" which is scalar field dependent. Several astrophysical observations significantly constrain the space and time variations of the gravitational constant, for instance in cosmology we know that since the recombination era |Grec−Gnow| Gnow < 5×10 −2 [30][31][32]. From the theoretical point of view, negative G eff associates the scalar with a ghost and indicates an instability [51]. In general, diverging G eff triggers a curvature singularity which blocks the system to evolve from an attractive gravity to antigravity regime [42,70,73,77,81].
On varying the action (2) with respect to the scalar field φ, one obtains where comma ( ,φ ) was introduced to denote a derivative with respect to the scalar field φ. We may eliminate R in the last equation by substituting in the trace of Eq. (4), the result is where T is the trace of the energy momentum tensor of the matter fields and E(φ) is the combination (3). In the case of minimal coupling this is just the Klein-Gordon equation (in curved spacetime) with V ,φ on the r.h.s. and no matter contribution. In the case of nonminimal coupling two additional aspects appear. First, the scalar field dynamics is largely governed by the effective potential [62,68,82] as it is (1+f (φ)) 3 E V eff,φ which in Eq. (7) plays the same role as V ,φ in the minimal coupling case. (Note that the factors 1 + f and E are related to the inverse of the effective gravitational constant (5) and the no ghost expression (3) and should not change their sign here.) Second, the matter fields also act as a source to the scalar field via T , like they perform as a source for the tensor field in the Einstein's equations (4) via T µν . However, conformal matter like electromagnetic radiation or ultrarelativistic particles for which the trace of the energy momentum tensor vanishes, T = 0, do not source the scalar field even in the nonminimal coupling case. As we will see, this nuance will be significant for cosmology as well.
Assuming the background to be spatially flat described by the Friedmann-Lemaître-Robertson-Walker (FLRW) metric, and the matter content to be homogeneously and isotropically distributed radiation and dust matter fluids, the field equations (4), (6) become [50,51,83] 3 (1 + f (φ)) H 2 = κ 2 ρ r + κ 2 ρ m +φ for the metric,φ for the scalar field, andρ for the radiation (ρ r ) and dust matter (ρ m ) densities, respectively. Here, H =ȧ a is the Hubble parameter, a denotes a scale factor and dot (˙) indicates a derivative with respect of the cosmological time t. This set of equations reduces to the familiar cosmology in GR with a cosmological constant Λ if the evolution of the scalar field halts,φ → 0, φ → φ * , such that [63,84] ( and The conditions (15), (16) guarantee that the φ dynamics is really stopped for good and is not regenerated in Eq. (12) by the slope of the effective potential or matter sources, respectively. When the scalar field is constant then in the Einstein's equation (4) only the value of the potential V (φ * ) appears, playing the role of the cosmological constant.
A simple way to satisfy (15), (16) is to have V ,φ and f ,φ vanish at the same value φ * . If the model functions are suitable, the GR limit may be realized by the solutions as an asymptotic state (past or future). Not every theory in the STG family is endowed with such f (φ) and V (φ) that their extrema overlap, in such cases the GR limit could be obtained in the regime where either the potential or matter can be neglected. If the conditions (15), (16) are met in a particular model, it would be convenient to have as well, as in this case the effective gravitational constant (5) coincides with the bare gravitational constant κ 2 , and the value of the potential tied with the cosmological constant. We may reparametrize the scalar field by adding a constant to φ and adjust the theory constants κ 2 and so that at φ * the relations (17) hold. (We can do such manipulation only for one point, i.e. one value of φ in a given model.) We may introduce the relative energy densities for radiation, dust matter, and the scalar field, respectively In the case of minimal coupling, f (φ) = 0, the quantities (18) reduce to the familiar relative densities of radiation, dust matter, and the scalar field, while in the GR limit (15) the latter becomes the relative density of cosmological constant. However, in the nonminimal case, the physical interpretation of these quantities should be exercised with caution since with the involvement of nonminimal coupling their values are not necessarily and always bounded to be less and equal to one. Still, assuming radiation and dust matter to be normal with nonnegative energy densities ρ r ≥ 0, ρ m ≥ 0, and ruling out the effective gravitational "constant" (5) being negative, Ω r and Ω m are positive definite. Tentatively, the first two terms in the definition of Ω φ in Eq. (18) may be called the relative kinetic energy density of the scalar field, while the last term there allows the interpretation of relative potential energy density of the scalar field. Under the condition that the potential can not be negative, the relative potential energy density is always nonnegative, but the relative kinetic energy density can under certain conditions dip below zero in the nonminimal coupling case. With the definitions (18), the Friedmann equation (11) obtains a very simple form, From here we can define radiation domination as the situation where Ω r is practically one while Ω m and Ω φ are negligible. Similarly, the dust matter and scalar field domination mean that Ω m and Ω φ dominate over the other energy densities in (19), respectively. We can also speak of scalar kinetic domination and scalar potential domination when the first two terms in Ω φ or the last term in Ω φ overpowers the others. As mentioned before, the relative kinetic energy of the scalar field can become negative in certain situations whereby the value of the other terms can exceed one. Additionally, we define the effective energy density and pressure of all components combined as Then the effective equation of state pertaining to the overall universe is given by As expected, in the GR limit, under radiation domination w eff = 1 3 , under matter domination w eff = 0, and under scalar potential (cosmological constant) domination w eff = −1. For accelerated expansion of the universe (ä > 0) one requires w eff < − 1 3 , while superacceleration (Ḣ > 0) is engendered by w eff < −1. De Sitter solution is characterized by w eff = −1 whereby the Hubble parameter remains constant and proportional to the cosmological term, H = Λ 3 . A minimally coupled scalar is not capable of triggering superacceleration and is bound to show expansions with −1 ≤ w eff ≤ 1 only (unless it is a ghost by nature). However, nonminimal coupling allows expansions sporting w eff < −1 and without necessarily running into instability [52][53][54][55][56][57]. For example, a wide class of STG models is known to approach the GR limit in the potential domination regime wherby w eff exhibits damped oscillations around the −1 value [39,62].

III. CONSTRUCTING THE PHASE SPACE
In the previous section, we have presented the STG cosmological equations for flat FLRW spacetime with radiation and dust matter. Due to their nonlinear nature, it is a daunting task to try solving them analytically in full generality. However, the methods of dynamical systems offer several tools to extract qualitative information about the behavior of solutions in the system. Let us begin here by writing the equations in the form of a dynamical system, and proceed to analyze the fixed points in the next section.
First, we write the Friedmann equation in a dimensionless form as denotes a derivative with respect to e-folds N = ln a. Then, let us introduce the following dimensionless variables in order to convert the cosmological equations into an autonomous system of equations, These variables are constrained by the Friedmann equation (23) which appears as This choice of variables warrants a few comments of motivation. The variables x and y are the usual expansion normalized variables [1-3] adopted for the setting of nonminimal coupling. They are clearly well suited for our purposes. For instance, in the radiation domination situation (when the dust matter and scalar field densities can be neglected), the Friedmann equation reduces to 1 = x. In other words, in the radiation domination setting x remains constant, which is exactly what should happen at a fixed point. Note that this does not imply that the radiation density ρ r itself is constant in cosmological time t. In this regime both ρ r and H evolve, but in a manner where x is constant. Similarly, in the dust matter domination case the Friedmann constraint imposes constant y = 1, characteristic of the dust matter fixed point. Obviously, we would not like to consider the Hubble parameter as a dynamical variable, since such choice would only show Minkowski and de Sitter solutions (constant H) as fixed points.
Next, one may analogously try to introduce the remaining dynamical variables in such a way that they describe the relative kinetic and potential energy densities of the scalar field [3,85], whereas the Hubble parameter gets completely absorbed by the variables and does not occur independently in the dynamical system. While doing that we would naturally get the scalar kinetic dominated (stiff fluid) and scalar potential dominated eras as fixed points. With some luck in picking a good potential, we may also witness scaling solutions as fixed points whereby the relative energy densities of one matter component and of the scalar field scale identically as the universe expands [86,87]. Such fixed point means a situation where the shares of Ω r,m and Ω φ in the Friedmann equation stay put, while ρ r,m and φ both evolve in cosmological time. However, the drawback with the approach described above is that in order to close the equations of the dynamical system (to make the system autonomous), one has to eliminate φ by expressing it in terms of the last dynamical variable (proportional to V ), which is not always easy or even possible [3].
In this paper, we have made φ dimensionless from the beginning and opt to keep it as an independent variable. This allows us to consider any function of f (φ) or V (φ) without worrying on how to close the dynamical system. The drawback is that we will not be able to capture some interesting behaviors in the system where the scalar field evolves in cosmological time (like the scaling solutions) as fixed points. This does not mean that such behavior is not present in the system any more, it just does not manifest itself as a fixed point. Also, the physical interpretation is now a bit less straightforward, as in general, the scalar field φ itself does not carry any direct physical meaning in gravity or cosmology, only certain invariant combinations of it like the energy density or gravitational constant (5) are observable [62]. (However, when the scalar happens to couple to other fields nongravitationally like Higgs, then its value would be relevant in particle physics, of course.) Finally, since the scalar field enters with second derivatives in the system (10)- (14), we need one more variable to account for its derivative. The simplest choice might seem to treatφ itself as a dynamical variable of the system [68]. However, with hindsight z =φ H seems a better choice, since in this case the phase space will show not only regular situations withφ = 0 as fixed points, but also singular situations whereφ and H diverge but z still vanishes can emerge as fixed points (at specific values of φ). The latter turns out to be crucial in our program to associate some unstable fixed points in the phase space with inflation. In addition, we have φ = z which makes writing out the dynamical system much more straightforward and simple.
These four variables encode the dynamics of the cosmological equations (10)- (14). From the definitions (25), we compute the derivatives with respect to e-folds, and then substitute inḢ andφ from the field equations (11) and (12),ρ r andρ m from the continuity equations (13) and (14), and finally H from the Friedmann equation (26). In the latter, we have picked the positive branch of the square root for expanding universe, and assumed H > 0 in all situations of interest (no "bounce", "turnaround" or asymptotic Minkowski scenarios for the universe), implying we can divide by H without extra concerns. As a result, we get the following autonomous system where the quantities were introduced for the sake of more compact notation and the overall factor E is just the no ghost expression (3).
From the outset, we note that the system (27)-(30) has two invariant submanifolds x = 0 and y = 0, as vanishing x implies x = 0 and vanishing y implies y = 0. In physics terms, this reflects the premise that if there is no radiation or dust matter present in the universe, these components will not be generated during further evolution. (We will come back to this premise in the discussion Sec. VIII.) It follows, that we can immediately deduce that global attractors of the system must lie on the intersection of the x = 0, y = 0 surfaces, i.e. in the regime where the radiation and matter contributions vanish.
The z = 0 surface does not appear as an invariant manifold, since nonzero y as well as nonzero V eff,φ (encoded in the last two terms 12A 12 − 6A 1 B 1 of (29)) would give dynamics to the scalar field (z = 0) even when z vanishes. More precisely, in the minimally coupled case z = 0 and V ,φ = 0 are sufficient to make z = 0. In the nonminimally coupled case z = 0 and V eff,φ = 0 are not sufficient to imply z = 0, one must also include y = 0 but not necessarily x = 0. The latter is so because in Eq. (29) the free term 12A 12 − 6A 1 B 1 also appears as a factor in front of the single x, hence when the former vanishes the whole x-term vanishes. Obviously, all these details just reflect the physics encoded in the scalar field equation (7). As already discussed there, the nonminimally coupled scalar field is sourced by V eff,φ as well as different types of matter, except for radiation.
The various cosmological parameters can be written in terms of variables (25) as follows The discussion after Eq. (19) implies that radiation domination means x = 1, y = 0, matter domination means x = 0, y = 1, while the scalar field domination happens when x = 0, y = 0. Several limits of the effective barotropic index w eff were observed after Eq. (22) and apply here as well.
The Friedmann constraint (26) puts restrictions on the domain of the phase space variables x, y, z, φ. We have assumed that matter components are usual with positive definite energy densities, ρ r , ρ m ≥ 0. We also expect the scalar potential to be positive definite, V > 0, ruling out zeros of the potential, since the current accelerated expansion suggests that at least some tiny nonzero value of the potential must be there, whichever is its actual physical origin. In addition, we only consider the positive definite values of the effective gravitational constant (5) with a possibility that G eff can diverge at the boundary of the phase space. Hence, the four dimensional phase space of the system (27)-(30) is given by As the phase space (34) is unbounded in some of the directions, we may lose sight of the interesting asymptotic features of the system pertaining to the infinite values of the variables. In order to have a complete picture on the global dynamics, we can compactify the phase space by applying the following Poincaré's transformation: In the compactified phase space The points at infinity in Ψ are associated with points on unit hypersphere The equations for the compactified dynamical system in x p , y p , z p , φ p can be obtained by using the transformation (35) on the system (27)- (30). We will not give the resulting equations here explicitly due to their extremely lengthy and unwieldy expressions. However, we will refer to these equations in the discussion of asymptotic fixed points, and use them numerically in plotting the trajectories for concrete examples in the following sections.

IV. FIXED POINTS AND THEIR PROPERTIES
In order to understand the behavior of the solutions in the system (27)-(30), we will carry out the dynamical systems analysis in a standard manner [1-3]. We extract the fixed points from the equations, compute the eigenvalues of the perturbed matrix around these fixed points, and find the eigenvectors corresponding to each eigenvalue. We also discuss the fixed points at infinity. In this section, we present the main mathematical results while the next three sections discuss more physics by looking at concrete examples with increasing complexity.
The fixed points (also "equilibrium points") of a dynamical system correspond to the solutions with no dynamics in terms of the variables of the system, i.e.
For a general coupling function f (φ) and general potential V (φ), the dynamical system (27)- (30) has three classes of fixed points which entertain different cosmological eras: points corresponding to radiation domination (R) at x = 1, y = 0, z = 0, points corresponding to matter domination (M) at x = 0, y = 1, z = 0, and points corresponding to scalar field domination x = 0, y = 0, z = 0 which behave effectively like de Sitter (dS). In certain cases, the phase space may also contain some other fixed points where the variables assume other values, but these cases depend on the specifics of the functions f (φ) and V (φ) which makes their general assessment cumbersome and not very informative. In the following, let us briefly characterize each of the three classes of points mentioned above.
R The radiation dominated points (1, 0, 0, φ r ) exist for any choice of the coupling function f (φ) and potential V (φ). In fact, they form a line in the phase space since the condition (38) is satisfied in radiation domination for any value of the scalar field, φ r (such feature does also occur when one employs some other definitions of dynamical variables [34,66,76]). The effective barotropic index at these points, w eff = 1 3 , takes the usual value for relativistic matter.
The eigenvalues of the corresponding perturbed matrix are with the corresponding eigenvectors where and all the quantities must be evaluated at the value φ r under consideration. We see that all the fixed points along the fixed line maintain the same characteristic of a saddle, where two directions are repulsive to the solutions (positive eigenvalues) and one is attractive (negative eigenvalue), although the precise orientation of the eigendirections changes point by point. The eigenvector corresponding to the zero eigenvalue is aligned along the fixed line itself, as should be the case, since it is impossible for the solutions to evolve along the fixed line. The attractive eigenvector (corresponding to the negative eigenvalue) lies on the surface y = 0 where dust matter vanishes. This implies that trajectories tend to approach the radiation domination point in a regime where dust matter component is negligible. The first repulsive eigendirection does also reside on the y = 0 surface, and guides the solutions away from the radiation domination to the de Sitter attractor when dust matter is not present. In the case of minimal coupling, the second repulsive eigenvector lies on the x − y plane and points directly towards the the matter domination fixed point. In the case of nonminimal coupling this eigenvector can deviate from the x − y plane. This tells that solutions with at least some tiny amount of dust matter present will evolve away from radiation domination towards matter domination or scalar field domination, depending on the model and initial conditions.
M The dust matter dominated points (0, 1, 0, φ m ) exist when φ m satisfies It means in the minimally coupled case the dust matter dominated points form a fixed line in the φ dimension like the radiation dominated points. However, in the case of nonminimal coupling, the dust matter dominated points exist only at certain values of φ, or may not exist at all in a given model, if the coupling f (φ) has no extrema. This is a reflection of the "attractor mechanism" whereby in the regime where the dust matter dominates and the scalar potential can be neglected, the scalar field will rest only at the GR limit (16) [34,36,40]. From Eq. (33), we see that all dust matter fixed points correspond to an unaccelerated universe with w eff = 0 like in GR.
The eigenvalues of the perturbed matrix are and the corresponding eigenvectors are given by where and all quantities must be computed at the value φ = φ m under consideration.
The fixed point M has the character of a saddle with one or two repulsive eigendirections and three or two attractive eigendirections, depending on the model function f (φ). In the case of minimal coupling where the M fixed points form a line since φ m is left arbitrary, the third eigenvalue vanishes and the related eigenvector points along this line, like in the case of the R points. The fourth eigenvector that corresponds to an attractive eigenvalue lies on the x − y plane and attracts the trajectories coming from the radiation domination point.
In the case of minimal coupling, the repulsive eigendirection of R and attractive eigendirection of M align, indicating a heteroclinic orbit between them (a trajectory running from one fixed point to another). However, the eigenvector that corresponds to another attractive eigenvalue does not lie on the x−y plane which implies that also trajectories coming from scalar field domination or infinity approach matter domination without passing through radiation domination. Further, the eigenvector that corresponds to the first repulsive eigenvalue lies on the surface x = 0. This signifies that the trajectories which repel away from this point correspond to vanishing radiation solutions and hence must be attracted towards a de Sitter point.
In a system with general coupling function, but without radiation and potential, approximate analytic solutions near this fixed point can be found in Ref. [40].
dS The scalar field dominated points (0, 0, 0, φ * ) exist at the values φ * where and recall we have assumed V (φ) = 0. Again the existence of such points in a given model depends on the functions f (φ) and V (φ). For the regular values of the potential, the condition (49) is satisfied at the extrema of V eff (8), in the case of minimal coupling at the extrema of V . Moreover, by comparing with Eq. (15) we see the existence of this fixed point coincides with the GR limit. In fact, the latter is not an arbitrary coincidence, but derives from how we have set up our variables and constructed our system. In the discussions of Sec. II, we defined the GR limit occurring at the value where the scalar field feels no force, which in the absence of matter is provided by the extremum of the effective potential. On the other hand, the absence of force is just the definition of a fixed point. As far as expansion is concerned, this point corresponds to an accelerated universe with de Sitter behavior w eff = −1.
The eigenvalues of the perturbed matrix are associated with the eigenvectors (taking first '+' and then '−' in the above) where and E was defined at the no ghost condition (3) and all the quantities must be evaluated at the value φ * under consideration.
In the case of minimal coupling the situation is quite clear. The dS point is an attractor (all eigenvalues negative) when it occurs at a local minimum of the potential, V ,φφ (φ * ) > 0, and a saddle with one repulsive eigendirection when it occurs at a local maximum of the potential, V ,φφ (φ * ) < 0. In particular, it is a stable node for 0 < V ,φφ V < 4 3 and a stable focus (spiral) when The linear stability analysis is inconclusive if V ,φφ (φ * ) = 0, but in such case, the method of center manifold theory can be employed to determine the behavior near the corresponding center manifold of the point. For minimal coupling, we present such analysis in the appendix A.
In the particular subcase where the minimally coupled scalar field has a constant potential (i.e. there is only a cosmological constant), the condition (49) does not fix the value of φ * and the dS fixed points form a line in the φ dimension of the phase space (like the R points). In such a situation, the eigenvalues reduce to 0, −3, −3, −4 and the corresponding eigenvectors are Note that the first eigenvector corresponding to the zero eigenvalue points along the φ dimension and signifies no dynamics along the line of fixed points.
In the nonminimal coupling case the dS fixed points have the character of an attractor, unless where one of the eigenvalues is positive. In fact, the condition (55) can be reconstructed as with the fixed point condition (49) substituted in. This accords with our observation after Eq. (8) that in the nonminimal coupling case the scalar field dynamics is largely governed by the effective potential. The fixed point occurs at the local extremum of the effective potential, it is an attractor at the local minimum and a saddle at the local maximum of the effective potential. If V eff,φφ vanishes, this point becomes nonhyperbolic and full analysis would again require the application of the center manifold theory, as in the minimally coupled case (in appendix A). However, we will not perform this analysis here, since in the general case the formulas become rather unwieldy and not very informative.
In (51), the first eigendirection is the only one that can be repulsive, it is orthogonal to the x − y plane and guides solutions towards or from a fixed point at another value of φ. The second eigendirection is attractive and also orthogonal to the x − y plane. Therefore, when a model allows e.g. two dS fixed points, one saddle and one attractor, most likely there will be a heteroclinic orbit between them, flowing from the saddle to attractor along the plane where radiation and matter are zero. The third eigendirection is also attractive and pulls the solutions from the side of the radiation dominated point R. The fourth eigendirection is attractive and collects the solutions from the matter dominated point M . One would expect that an unstable dS saddle is very useful in engendering reasonable cosmic history, providing a natural setting for inflation. Indeed, the three attractive eigendirections would gather many solutions from the phase space into the vicinity of this point, these solutions would then undergo almost de Sitter type accelerated expansion as they pass by the fixed point, and then leave the neighborhood of this point following the first (repulsive) eigendirection, whereby the inflationary expansion terminates. For an extensive set of initial conditions, this point would provide a transient phase of accelerated expansion. However, there is a complication. The first eigendirection does not immediately point to the radiation dominated point, that according to expectations should take over the relay baton in the cosmic race. Instead, if the system has besides a saddle dS also an attractor dS, many of the solutions would flow directly from the saddle to the stable node dS, not giving much chance for radiation or matter to increase their share in the cosmic budget. Yet, both radiation and matter dominated points were endowed with attractive eigendirections, these could still pull a chunk of trajectories into their zones of influence. How this scenario precisely unfolds in the phase space depends on the relative positions of the fixed points. In the general context we can at best state that for the dS point to play a role in inflation, it has to be a saddle, i.e. the condition (55) must hold. Otherwise, the dS point is a stable node, probably responsible for dark energy. As the dS attractor lies on the intersection of invariant submanifolds, it is, therefore, representing a global attractor irrespective of the choice of initial conditions.
In a system with general nonminimal coupling and potential, but without radiation and dust matter this fixed point was described before in Refs. [69][70][71], and the analytic form of approximate solutions near it can be found in Refs. [39,62].
dS ∞ A closer consideration of the dS point condition (49) reveals, that it can also exist at the scalar field infinity, |φ| → ∞. Indeed, by inspecting the system (27)-(30) for minimal coupling, we see that the limit x = 0, y = 0, z = 0, |φ| → ∞ can imply x = 0, y = 0, z → 0, φ = 0, i.e. there is a fixed point, provided For instance, this condition is satisfied for power law potentials. In the nonminimal coupling case, a model with quadratic coupling function f = ξφ 2 and quartic potential V = V 0 φ 4 + Λ is endowed by such an asymptotic fixed point, for example. Other models possessing it can be also constructed. The eigenvalues characterizing the point need to be properly recalculated in the compact variables (35) with considerable extra effort, however, quick numerical checks of the evolution near this point in the compact phase space (36) reveal that the point is a saddle when the effective potential (8) has an asymptotic maximum there, and an attractor, when it has an asymptotic minimum.
It is interesting that in this regime we have de Sitter like accelerated expansion with w eff = −1, but neither the scalar field nor the Hubble parameter are constant. In fact, the Friedmann constraint (26) allows V → ∞, H → ∞ while z =φ H → 0. Hence strictly speaking we might better call these points "asymptotic de Sitter" or "quasi de Sitter", in contrast with the true de Sitter points where H is a constant. Such asymptotic behavior of the solutions has been noted before [68], but not in connection with a fixed point. Our choice of dynamical variables brings out this regime as a fixed point. One should also note, that contrary to the regular dS fixed point, the asymptotic dS ∞ point does not satisfy the condition (15) of GR with a cosmological constant. Nevertheless, it is very useful in the description of inflation.
dS s Finally, our variables also formally permit to satisfy the fixed point condition (49) at the limit f (φ s ) → −1 where the effective gravitational constant (5) becomes singular. At this value φ s , the effective potential (8) turns singular as well and the situation is in some ways similar to dS ∞ , i.e. de Sitter like accelerated expansion with w eff = −1, but H → ∞ while z =φ H → 0. Of course, this point is only available in the nonminimally coupled case. Numerical checks show this point behaves as a saddle, consistent with the property that the singular value of the effective potential can only be a maximum according to our assumptions. Although the variation of the effective gravitational constant is rather restricted by observations in the late universe, the possibility of this regime as a rallying point for inflation can be entertained. Like dS ∞ , the dS s point does not satisfy the condition (15) of GR limit either.
The determination of the eigenvectors associated with each critical point suggests that scalar-tensor theories are able to explain the thermal history of the universe, i.e. going from a scalar field dominated inflationary expansion to radiation dominated era to matter dominated era to the final scalar field dominated accelerated expansion (dark energy) epoch. This sequence can be realized generically in a model where the corresponding fixed points exist as three saddles and an attractor, respectively. These fixed points shape the flow of trajectories and ensure that the passing by solutions exhibit the expansion rate appropriate for the relevant era. While the radiation domination, matter domination, and dark energy domination points are generated quite easily, not every STG model provides de Sitter point suitable for inflation, i.e. a saddle where the repulsive direction guides the trajectories away, so that the early accelerated expansion can end. Still, in many models de Sitter saddle is available, either as a regular local maximum (dS), an asymptotic maximum (dS ∞ ), or a singular value (dS s ) of the effective potential.

V. EXAMPLE: GENERAL RELATIVITY WITH COSMOLOGICAL CONSTANT
To get a glimpse of how the physics is captured in our dynamical system, let us begin by considering the simple case of GR with radiation and matter, while taking φ ≡ constant, f (φ) ≡ 0, V (φ) ≡ Λ. In this case, only the variables x and y are dynamical and the system (27)-(30) reduces to Assuming the cosmological constant is positive, Λ > 0, and denoting Ω Λ = Λ 3H 2 , the Friedmann constraint implies that the two dimensional phase space is bounded, delimited by Ψ Λ = (x, y) ∈ R 2 0 ≤ x, y < 1; 0 ≤ x + y < 1 .
The same system has been presented before in the literature [3,25], and here we invoke it to guide the reader into more involved cases which follow.
From the system (58)- (59), it is easy to see that there are three fixed points situated at the boundary of the phase space, characterized by their properties as follows: R at (1, 0) is a repeller and corresponds to radiation dominated universe, Ω r = 1, w eff = 1 3 ; M at (0, 1) is a saddle and corresponds to dust matter dominated universe, Ω m = 1, w eff = 0; dS at (0, 0) is an attractor and corresponds to de Sitter universe, Ω Λ = 1, H 2 = Λ 3 , w eff = −1. The properties of the fixed points can be inferred from the eigenvalues of the linearized system at the fixed points, as a particular subcase in the discussion of the previous section, but are also directly evident from the phase diagram on Fig. 1(a). It should be stressed that the solutions sitting exactly at the fixed points are relatively uninteresting, corresponding here to universes filled with only radiation, dust matter, and cosmological constant, respectively. The importance of fixed points, however, lies in how they determine the behavior of solutions around them: repellers act like sources as nearby trajectories depart from them, attractors act like sinks as nearby trajectories arrive at them, while the saddles repel in one direction and attract in another direction. On Fig. 1(a), we see how the generic trajectories start near the radiation domination fixed point (R), evolve towards the matter domination fixed point (M), but turn away from it and end up at the de Sitter point (dS). At the boundaries, there are also three special ("heteroclinic") trajectories, one starting at R and running directly into M, another starting at R and running into dS, and third starting at M and running into dS. These correspond to universes without cosmological constant, matter, and radiation, respectively. Technically speaking, the points R and M, as well as the orbit connecting them are not part of the phase space Ψ Λ if Λ > 0, but to understand the physics it is certainly illuminating to include them. The de Sitter point dS is a global attractor as it lies in the intersection of invariant submanifolds x = 0 and y = 0.
The behavior of generic trajectories in the phase space is in qualitative agreement with the history of our observable universe. Earlier the universe was dominated by radiation, then as radiation dilutes the dust matter takes the role of the leading component for a certain period, while in the future both radiation and dust matter densities drop and the reign of the cosmological constant begins. At some point before the complete domination of the cosmological constant, the universe starts to expand in an accelerated fashion, as the trajectories enter into w eff < − 1 3 regime, marked by the shaded area on Fig. 1(a). On the figure, the late history of our observable universe is best represented by the green trajectory which comes very close to the M point. We can observe the evolution of the relative densities of radiation, matter, and cosmological constant, as well as the effective barotropic index, along this trajectory on Fig.  1(b). (Here we employ a useful convention to follow the evolution in terms of the logarithm of redshift, related to the scale factor by 1 + z = 1 a , whereby the present scale factor is normalized, a(t now ) = 1. In this way, different epochs can be conveniently and visibly captured on a single plot, as in e.g. Refs. [5,73].) The majority of other trajectories in the phase space correspond to larger values of the cosmological constant, whereby the matter component in the universe has less chance to dominate.
The system (58)-(59) gave us a correct qualitative picture of the latter stages of the evolution of the universe, owing to the properties of the fixed points. What is missing is the epoch of early accelerated expansion, aka inflation. Although there is an accelerating region in the phase space surrounding the dS fixed point, it is unable to provide early time acceleration, since all the trajectories only enter it, while none goes out from this patch. This follows from the attractive nature of the dS fixed point, which draws all trajectories into its fold and lets none to escape. For an early transient period of accelerated expansion, we need another fixed point which is of the dS type, but has the characteristics of a repeller or at least a saddle.

VI. EXAMPLE: MINIMALLY COUPLED SCALAR FIELD WITH QUARTIC POTENTIAL
As the next step, let us consider minimally coupled scalar with a quartic potential, where V 0 > 0, Λ > 0 as an example. The corresponding dynamical system, although cast in different variables was considered before in the case of absent matter fluids in Ref. [68] and with a single fluid but without cosmological constant in Ref. [88]. Compared to the previous section, including a dynamical scalar field makes the phase space four dimensional. The positivity of the potential imposes a constraint via the Friedmann equation (26), From the notes in Sec. IV, we can easily find out the fixed points: R at (1, 0, 0, φ) form a line along the φ direction, all these points are saddles in character.
M at (0, 1, 0, φ) also form a line along the φ direction, and are saddles as well. dS 0 at (0, 0, 0, 0) corresponds to de Sitter universe with H 2 = Λ 3 . Using the results of appendix A, we can determine it to be a stable node.
The behavior of one representative phase trajectory in a finite patch of (x, y, φ) subspace is depicted in Fig. 2(a), where the region of accelerated expansion is indicated by the blue contour and the fixed points shown in red. It can be seen that the trajectory evolves from the accelerating region towards a radiation dominated point R, then to a matter dominated point M , and eventually returns to the accelerated region and settles down approaching the de Sitter point dS 0 . Recall the result of Sec. IV that the attractive direction of the radiation domination point was lying in the y = 0 plane where dust matter density vanishes. On the plot 2(a) this reflects in the fact that the trajectory approaches R by being close to the y = 0 plane. When the universe evolves from inflation to radiation domination the matter density stays low, it will increase only after passing by the radiation domination point. Analogously, the repulsive direction of M lies in the x = 0 plane, thus after the matter domination, radiation density has no tendency to increase again. We just notice how such basic thermodynamic features are encoded in the phase space properties.
The time evolution of energy density parameters along with the effective equation of state for this solution can be read off from Fig. 2(b). Here, we see how the universe evolves from the scalar field dominated inflationary era (Ω de ≈ 1, w eff ≈ −1), passes on to the radiation domination epoch (Ω r ≈ 1, w eff ≈ 1 3 ), followed by the matter domination epoch (Ω m ≈ 1, w eff ≈ 0), and eventually settles in an accelerated expansion regime dominated by the scalar field as dark energy (Ω de ≈ 1, w eff ≈ −1). Although the history is presented here on a logarithmic scale, we can still notice how the fixed points are instrumental in making each of the epochs pronounced enough. Indeed, near the fixed points the flow along the trajectory (the evolution of the solution measured in time, scale factor, or redshift) "slows down" as the "force" (r.h.s. of the system (27)- (30)) becomes small and hence the derivatives of the dynamical variables become small too. On the other hand, the transition from one epoch to another occurs relatively "quickly", since while the trajectory passes from the vicinity of one fixed point to another, the "force" is bigger in the middle regions and the dynamical variables evolve "faster". The finite phase space 2(a) does not reveal the earlier path of the trajectory, where does it originate from. Thus let us check a hybrid plot 2(c), where the variables z, φ are made compact in analogy with (35), i.e. mapped to the range (−1, 1). Tracing the same trajectory backwards takes us into the vicinity of the asymptotic fixed point dS − ∞ . Therefore this point, as well as its cousin dS + ∞ are indeed responsible for inflation, because around these points the trajectories can exhibit de Sitter like expansion (w eff ≈ −1) which lasts sufficiently long. Since this point is a saddle, there is only one trajectory that exactly starts at this point. It is the one that flows out of dS ∞ along the repulsive eigendirection. All other trajectories in the phase neighbourhood come from somewhere else (other asymptotic regions), but closely approach this particular solution as they evolve, since all other eigendirections are attractive. In principle, we could follow our example trajectory further back in time, but as we already are in the asymptotics the numerical errors mount. The main feature is that inflation is a generic phenomenon for a wide range of initial conditions. For a true glimpse of the global dynamics of the system, we should compactify the full phase space by the Poincaré transformation (36). The asymptotic fixed points reside at the boundary of the compact space, dS + ∞ (0, 0, 0, 1) and dS − ∞ (0, 0, 0, −1). Fig. 2(d) depicts the x = y = 0 slice through the 4-dimensional compact phase space, where points on the circumference of the circle represent the asymptotic infinity. The gray zone is the region prohibited by the Friedmann constraint (63), and the blue shaded area stands for the region of accelerated expansion. We have also plotted a projection of the direction of the trajectories' flow, as they pass through or along this phase space slice. By closely observing the phase flow and also doing quick numerical investigations, we can recognize that the points dS + ∞ and dS − ∞ behave as saddle points, since the trajectories come close to them and then evolve away, eventually heading towards the attractor de Sitter point (dS 0 ). The logic of the phase flow also implies that the dS ± ∞ and dS 0 points are connected by a heteroclinic orbit lying exactly on the x = y = 0 plane. We do not see the R and M points here because they reside in the dimensions that are not shown on this plot. However, before reaching the global attractor dS 0 , the majority of the trajectories also wander in the x and y dimensions and pay a visit to these points as well. Putting all this information together in a global picture, we conclude that the basic cosmic history of inflation, radiation domination, matter domination, and final dark energy occurs as a rather generic property for the trajectories of the model, guaranteed by the suitable structure of fixed points.

VII. EXAMPLE: QUADRATIC NONMINIMAL COUPLING, QUARTIC POTENTIAL
As the main example, let us consider the case of quadratic nonminimal coupling and quartic potential with an additional constant term, It can be considered as a simplified version of the nonminimally coupled Higgs model [46] as the quartic piece dominates at large field values and is responsible for early universe inflation, while the constant piece dominates late universe as dark energy. The model is also interesting because it allows us to witness attractive and repulsive regular dS fixed points as well as asymptotic dS ∞ and singular dS s fixed points, described in Sec. IV. Without the radiation and dust matter components, the model (64) was analyzed before in Ref. [68] in different variables. The predictions for inflationary observables can be found in the literature [89][90][91] or quickly computed from the generic algorithm [92]. The phase space is still 4-dimensional with a few extra features and novel fixed points compared to the minimally coupled case. The positivity of the potential imposes a constraint via the Friedmann equation (26), In the case of ξ < 0 the positivity of the effective gravitational constant (5) also sets a limit while that boundary is marked by the effective potential (8) becoming singular.
Recalling the general results of Sec. IV, the fixed points R and M are independent of the model parameters, but the set and character of de Sitter points is determined by the conditions (49), (55) and depends on the value of the nonminimal coupling parameter ξ, most easily to understood in reference to the shape of the effective potential (8) plotted on Fig. 3. In summary, the fixed points are the following: R at (1, 0, 0, φ) form a line along the φ direction, all these points are saddles in character.
M at (0, 1, 0, 0) is just one point (not a line as in the minimally coupled case) matching the GR limit (16), it is a saddle.
dS 0 at (0, 0, 0, 0) is a saddle for ξ > 0 when it resides at the local maximum of the effective potential, and an attractor for ξ < 0 when it sits at the local minimum of the effective potential. In the latter case it is a stable node if − 3 16 < ξ < 0 and stable focus if ξ < − 3 16 . It corresponds to the GR limit (15).
dS ± at (0, 0, 0, ± Λξ V0 ) exist only for ξ > 0 and are attractors occupying the local minimimum of the effective potential. More precisely, their nature is stable focus for all values of ξ. They also satisfy the GR limit (15).
dS ± ∞ at (0, 0, 0, ±φ * ) are situated in the asymptotics φ * → ∞. These points have a meaning only for ξ > 0 when they coincide with the asymptotic maximum of the effective potential and correspond to saddles. (For ξ < 0 the asymptotic region is unphysical.) dS ± s at 0, 0, 0, ± 1 √ −ξ exist only in the case of negative nonminimal coupling ξ < 0. These points are saddles and correspond to a singularity of the effective potential.
Note that the emergence of dS ± happened by the virtue of adding a positive constant Λ to the potential. Without it, the system would have had only one regular de Sitter fixed point dS 0 at the minimum of the effective potential. Equipped with the general knowledge of the phase space and fixed points, let us illustrate the cosmological dynamics and follow a particular trajectory in the phase space, treating the ξ > 0 and ξ < 0 cases separately.
For positive nonminimal coupling, ξ > 0 the results are depicted on Fig. 4. Here 4(a) shows a finite patch of the (x, y, φ) subspace where a trajectory flows from the accelerated region towards to the radiation domination point R, then turns to the matter domination point M , passes by the saddle de Sitter point dS 0 , and finally ends up at de Sitter attractor dS + . The evolution of the relative energy densities and effective barotropic index is given on Fig.  4(b). A hybrid plot 4(c) has the range of the scalar field compactified to a finite interval −1 ≤ φ p ≤ 1, and we see how the trajectory actually goes back to the asymptotic fixed point dS − ∞ . The boundary of accelerating region w eff < − 1 3 , indicated by the blue wireframe, of course includes all dS points and extends further the x − y directions near these points. This can explain the "bump" in the w eff on Fig. 4(b). When the trajectory enters in the scalar field dominated era and first comes closer to dS 0 the effective barotropic index gets close to −1 as well, but while the trajectory then cruises away from dS 0 towards dS + the index bounces little bit back as the edge of the accelerating region is nearby for a while (see Fig. 4(c)). In fact, the final approach to dS + is marked by damped oscillations of w eff below and above −1, i.e. periods of minute superacceleration. The first of such periods is barely noticeable on Fig. 4(b), while in each period to follow the superaccelerating effect is an order magnitude less pronounced than at the previous one. Such damped oscillations in the dynamics are a characteristic feature of the focus nature of the fixed point dS + . On the plots 4(d) and 4(f) of the x p = y p = 0 slice of the fully compactified phase space (36), the superaccelerating region w eff < −1 is shaded yellow. The plot 4(f) shows a projection of the example trajectory onto this slice, and we see how it indeed lies on the superacceleration zone. The actual trajectory, however, does have nonvanishing radiation and matter, and it approaches dS + not exactly on the x p = y p = 0 plane, but spirals down to it from these dimensions. At last, the Fig. 4(e) shows a projection of the trajectory onto the z p = φ p = 0 slice of the compact phase space. We can compare with Fig. 1(a) and witness the correct sequence of cosmological eras to play out, viz. from a transient accelerating era, to radiation domination, to matter domination to late accelerating era of dark energy, as supplied by the fixed points dS − ∞ → R → M → dS + . Due to the symmetry of the effective potential a similar cosmological sequence can be also realized by the fixed points dS + ∞ → R → M → dS − . In the negative coupling case, ξ < 0, the story is basically the same with few modifications, as depicted in Fig.  5. There is just one regular de Sitter point dS 0 acting as the final attractor. Also, the range of the scalar field is constrained by |φ| < 1 √ −ξ , and the role of the dS saddle point is taken by dS s , occurring at the boundary of the allowed region where the effective potential diverges. Indeed, on the plots, we can see how an example trajectory exhibits the eras of inflation, radiation domination, matter domination, and dark energy, as it flows guided by the fixed points dS − s → R → M → dS 0 . The final point in this particular example is a stable node and there is no region of superaccelerated expansion for this choice of parameters. (However, for large negative values of ξ where dS 0 is a stable focus, there would also occur a superaccelerated region in the physical phase space and similar to the case where ξ > 0, one can observe damped oscillatory behavior about w eff = −1.) Due to the symmetry of the effective potential, a similar cosmological sequence can be also realized by the fixed points dS + s → R → M → dS 0 . The last extra subplot 5(g) shows how the trajectory does indeed go back close to the dS − s point and not at where would be dS − ∞ . A red curve on that subfigure indicates the limit where the effective gravitational constant diverges and the effective potential becomes singular. The unphysical region beyond that is not dynamically accessible from the physical part of the phase space. As explained before in the minimally coupled case, the saddle nature of dS s implies that only one particular phase trajectory originates from it, while all other trajectories in the vicinity converge to this particular solution. In principle, we could try to follow our example solution backwards beyond the immediate neighborhood of dS − s , but because of the singularity in the effective potential, the numerics becomes unreliable. In any case, the earlier adventures of the trajectory are not very relevant for subsequent history, as the dS − s repulsive eigendirection guides all the nearby trajectories irrespective of their origin or way of approach.
Although the two examples give the same qualitative picture of expansion history, the evolution of the effective gravitational constant (5) is radically different in them. In the positive coupling case (ξ > 0), inflation takes place at very small values of G eff , which then grows around the radiation domination era and reaches a maximal value at the matter domination, but then drops again when the solution reaches dS + . The drift of the effective gravitational constant in the matter domination era and later is rather restricted observationally, but can be suppressed if Λξ 2 V0 1, i.e. the fixed points M and dS + are sufficiently closely aligned in the φ direction of the phase space. On the contrary, in the negative coupling case (ξ < 0) inflation takes place at very large values of G eff , which then drops around the radiation domination era and reaches a minimal value at matter domination. Since the matter domination M and final dS 0 occur at the same value of φ, the evolution of the effective gravitational constant in the late universe would be marginally small for generic solutions (as they come close to the direct heteroclinic orbit from M to dS 0 ).
Because unlike matter domination, radiation domination is not governed by a single fixed point but rather a line of fixed points extending in the φ direction, there is not much control on the value of the scalar field in this era, except that the evolution of φ in time should somewhat slow down at the peak of the radiation domination when a trajectory is most close to the line of R fixed points. For an arbitrary trajectory, the scalar field value will evolve before and after the radiation domination peak. This has implications not only on the effective gravitational constant (5) but also on the relationship between the dust matter and radiation densities and the respective relative densities (18), which also evolve with φ. Therefore, it is interesting to speculate that the slightly varying scalar field value around the radiation domination era could be related to the observational puzzles of H 0 [93][94][95][96] (for investigations in this direction see e.g. [97][98][99]) or even to the EDGES 21 cm line [100].

VIII. DISCUSSION
Contemporary cosmology is proliferated by a plethora of models proposed to describe the phenomena of inflation and dark energy. Most of these models are designed to have an effect at either high or low energy scales. However, a true model of fundamental physics must deliver correct cosmic history from early to late eras, and better exhibit this history as a generic property for a large class of solutions, not only certain particular ones. In the present work, we address this issue in the context of generic scalar fields minimally or nonminimally coupled to gravity, by employing the methods of dynamical systems with a useful combination of variables, whereby the principal cosmic eras appeared as fixed points with appropriate properties. To summarize, a viable cosmological model should have the following features: 1. inflation supported by a de Sitter fixed point with one repulsive eigendirection, i.e. either a regular dS point at φ * such that (49) and (55) are satisfied, or an asymptotic dS ∞ or singular dS s described in Sec. IV, these collect a large set of solutions and guide them to an inflationary path; 2. radiation domination fixed point which is a saddle, it is always guaranteed to exist; 3. dust matter domination fixed point which is also a saddle, for nonminimal coupling, it exists provided there is a value φ m which satisfies (44); 4. final de Sitter fixed point manifesting as dark energy filled universe, when the condition (49) is satisfied but the condition (55) is not, all the eigendirections are attractive as necessary.
It turns out that the main crux in confirming a good model is to check whether it has at least two de Sitter fixed points, one with a repulsive direction and another totally attractive. These fixed points can be understood to reside at the local or asymptotic maxima and minima of the effective potential (8). If a model possesses the fixed points listed above, it has a good chance to reproduce qualitatively correct cosmic background evolution by a large set of its solutions. This list of fixed points provides the necessary barebones for a good model, but does not give an exhaustively sufficient condition, since there might be some other fixed points or features at the boundaries of the allowed phase space, which might alter or ruin the correct dynamics. However, as our variables allow to describe inflation in terms of a saddle de Sitter point, the question of how generic are the initial conditions for inflation and the eras that follow, can now be mathematically framed using the tools of dynamical systems, left for further work.
It must be also noted that in our approach we considered all the cosmic components to be immutable and without interaction or change in character. This was a very strong simplifying assumption, which does not hold in real cosmology. First, the stage of inflation should end with the process of reheating where the scalar field kinetic energy gets transferred into radiation and matter fluids. Second, the massive matter particles which later behave as collisionless nonrelativistic dust were at earlier times relativistic and behaved like radiation. These two phenomena should be included in the treatment of the full history of the universe as suitable interactions between the components, but we leave the engineering of these features for further work. When done properly the present set of fixed points should remain intact, but their basins of attraction would widen.
The aim of the current paper was to present the generic approach, while the systematic analysis of different nonminimal coupling functions and potentials enabling viable cosmic history is left for another publication. In the models with nonminimal coupling, the effective gravitational constant becomes dynamical and brings in additional aspects to pay attention to. Of course, background history is just one facet of observational viability, the other important issue is the treatment of perturbations, which remains beyond the scope of the current work. Besides cosmology, a correct model of gravity must also pass the tests coming from the Solar System, black holes, neutron stars, gravitational waves, etc. Still, for any scalar-tensor model of interest the results presented here do offer a quick viability check of basic cosmology, and allow to focus on next issues.
Finally, let us emphasize that the method presented here for scalar-tensor type of theories is completely general and can be easily adjusted to analyze other classes of theories as well. It should be rather straightforward to treat models with a scalar field nonminimally coupled to Palatini Ricci scalar [101], Gauss-Bonnet invariant [102], torsion scalar [103,104], nonmetricity scalar [105], etc., in a rather analogous manner, while in principle any more complex generalizations involving derivative couplings, higher derivative terms, multiple fields, and so on, can be subject to similar analysis as well.
The system of equations (27)- (30) can be recast in terms of these new set of variables as where g 1 , g 2 , g 3 , f are polynomials of degree greater than 2 in (X, Y, Z, Φ), with g 1 = 9 Z 2 X + 4 X 2 + 3 XY (A1) g 2 = 9 Z 2 Y + 4 XY + 3 Y 2 (A2) where V is a function of new variable Z + Φ (which is φ in terms of old variables) and D denotes derivative of V with respect to Z + Φ. The local center manifold is given by where h 1 , h 2 are approximated as respectively. The invariant property of the center manifold implies that the function h has to satisfy a quasilinear partial differential equation given by The dynamics of the reduced system is eventually determined by the equation and hence Therefore, if D [1] J. Wainwright and Ellis G. F. R., Dynamical systems in cosmology (Cambridge University Press, 1997).