Averaging Generalized Scalar Field Cosmologies II: Locally Rotationally Symmetric Bianchi I and flat Friedmann-Lemaître-Robertson-Walker models
Genly Leon, Sebastián Cuéllar, Esteban González, Samuel Lepe, Claudio Michea, Alfredo D. Millano
AAveraging Generalized Scalar Field Cosmologies II: Locally Rotationally SymmetricBianchi I and Flat Friedmann-Lemaˆıtre-Robertson-Walker Models
Genly Leon, ∗ Sebasti´an Cu´ellar,
1, †
Esteban Gonz´alez,
2, ‡
Samuel Lepe,
3, §
Claudio Michea,
1, ¶ and Alfredo D. Millano ∗∗ Departamento de Matem´aticas, Universidad Cat´olica del Norte,Avda. Angamos 0610, Casilla 1280 Antofagasta, Chile Universidad de Santiago de Chile (USACH), Facultad de Ciencia, Departamento de F´ısica, Chile. Instituto de F´ısica, Facultad de Ciencias, Pontificia Universidad Cat´olica de Valpara´ıso, Av. Brasil 2950, Valpara´ıso, Chile (Dated: February 15, 2021)Scalar field cosmologies with a generalized harmonic potential and a background matter with a barotropicEquation of State with barotropic index γ are investigated for Locally Rotationally Symmetric Bianchi I andFlat Friedmann-Lemaˆıtre-Robertson-Walker (FLRW) metrics. Using methods from the Theory of Averaging ofNonlinear Dynamical Systems and numerical simulations are proved that the full equations of the time depen-dent system and their corresponding time-averaged versions have the same late-time dynamics. Therefore, thesimplest time-averaged system determines the future asymptotic of the full system. In particular, depending onthe metric and the values of the barotropic index γ are found generic late-time attractors of physical interestslike: flat quintessence dominated FLRW Universe and a flat Friedman Einstein-de-Sitter solution. With this ap-proach the oscillations entering the full system through Klein-Gordon equation can be controlled and smoothedout as a time-dependent perturbative parameter given by the Hubble factor H goes monotonically to zero. PACS numbers: 98.80.-k, 95.35.+d, 95.36.+xKeywords: Generalized scalar field cosmologies, Anisotropic models, Early Universe, Equilibrium-points, Harmonic oscilla-tor
Contents
1. Introduction
2. The model
3. Averaging scalar field cosmologies k =
0. 6
4. Qualitative analysis of averaged systems k =
0. 94.2.1. Late-time behaviour 10
5. Conclusions Acknowledgements A. Proof of Theorem 1 B. Numerical simulation ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] § Electronic address: [email protected] ¶ Electronic address: [email protected] ∗∗ Electronic address: [email protected]
References
1. INTRODUCTION
Mathematical methods have been widely used inCosmology. For example, in reference [1] the methodof Lie symmetries was applied to the Wheeler-De Wittequation in Bianchi Class A cosmologies for minimallycoupled scalar field gravity and Hybrid Gravity in Gen-eral Relativity (GR). Several invariant solutions weredetermined and classified according to the form of thescalar field potential by means of the first principleof the existence of symmetries. In reference [2] wasused a model-independent criterion based on first in-tegrals of motion, and in reference [3] were used dy-namical symmetries of the field equations in order toclassify the Dark Energy (DE) models in the contextof scalar field (Quintessence or Phantom) Friedmann-Lemaˆıtre-Robertson-Walker (FLRW) cosmologies. Byusing Noether symmetries the system was simplifiedand its integrability was determined in [2]. For theexponential potential, as well as, some types of hy-perbolic potentials were found extra Noether symme-tries apart of the conservation law, suggesting that thesepotentials should be preferred along the hierarchy ofscalar field potentials. In [3] the requirement that thefield equations admit dynamical symmetries resultedin two potentials one of them is the well known Uni-fied Dark Matter (UDM) potential and another hyper-bolic model. In reference [4] was used a mathemati-cal approach based on the constraints systems in order a r X i v : . [ g r- q c ] F e b to reconstruct the equation of state and the inflationarypotential for the inflaton field from the observed spec-tral indices for the density perturbations and the tensorto scalar ratio. In reference [5] was used an algorithmto generate new solutions of the scalar field equationsin homogeneous and isotropic universes. Solutions canbe found for pure scalar fields with various potentials inthe absence and presence of spatial curvature and otherperfect fluids. A series of generalizations of the Chap-lygin Gas and bulk viscous cosmological solutions forinflationary universes were found. In reference [6] wasstudied the f ( T ) cosmological scenario. In particular,analytical solutions for a isotropic and homogeneousUniverse containing a dust fluid and radiation and foran empty anisotropic Bianchi I Universe were found bymeans of the method of movable singularities of dif-ferential equations. For the isotropic Universe, the so-lutions are expressed in terms of a Laurent expansionwhile for the anisotropic Universe are found a family ofexact Kasner-like solutions in vacuum. In reference [7]was performed the symmetry classification of the KleinGordon (KG) equation in Bianchi I spacetime. It wasapplied a geometric method which relates the Lie sym-metries of the KG equation with the conformal algebraof the underlying geometry. Furthermore, by means ofLie symmetries which follow from the conformal al-gebra, which are also Noether symmetries for the KGequation, were determined all the potentials in whichthe KG equation admits Lie and Noether symmetries.The Lie admitted symmetries were used to determinethe corresponding invariant solution of the KG equa-tion for specific potentials. Also was solved the prob-lem of classification of Lie/Noether point symmetriesof the wave equation in Bianchi I spacetime and invari-ant solutions of the wave equation were determined. Inreference [8] was developed a new method in order toclassify the Bianchi I spacetimes which admit Confor-mal Killing Vectors (CKV). The method was used tostudy the Conformal Algebra of the Kasner spacetimeand other Bianchi type I matter solutions of GR.Other useful mathematical methods applied to Cos-mology are Asymptotic Methods and Averaging The-ory [9–15]. These methods were used in Cosmologyfor example in [16–23] with interest in early and late-time dynamics. In reference [24] were studied LocallyRotationally Symmetric (LRS) Bianchi type III cos-mologies with a massive scalar field by using the The-ory of Averaging of Nonlinear Dynamical Systems tocontrol the oscillations entering the system through theKG equation. In reference [25] it was presented a theo-rem about large-time behaviour for solutions of a gen-eral class spatially homogeneous cosmologies with os-cillatory behaviour based on a first order approximationof H when H is non negative and monotonic decreasingto zero. As in reference [24] the main idea is to con-struct a time-averaged version of the original system,and prove that it shares the same late-time dynamics of the original system. It turns out a physical parame-ter, say, the quantity H denoting the Hubble parameterwhich plays the role of a time dependent perturbationparameter; controlling the magnitude of the error be-tween the solutions of the full and the time-averagedproblems. The analysis of the system is therefore re-duced to study the corresponding averaged equationsusing dynamical systems techniques [26–54].Inspired in the works [16–19] we started in refer-ence [20] the program “Averaging Generalized ScalarField cosmologies” which consists in using Asymp-totic Methods and Averaging Theory to obtain rele-vant information about the solution’s space of scalarfield cosmologies in the presence of the backgroundmatter with a barotropic Equation of State (EoS) withbarotropic index γ minimally coupled to a scalar fieldwith generalized harmonic potential. This researchprogram has three steps according to the three casesof study: (I) Bianchi III and Open FLRW model, (II)Bianchi I and Flat FLRW model and (III) Kantowski-Sachs and Closed FLRW.Relevant results for the aforementioned programwere presented in reference [19], where interactingscalar field cosmologies with generalized harmonic po-tentials were studied for flat and negatively curvedFriedmann-Lemaˆıtre-Robertson-Walker (FLRW) met-rics and for Bianchi I metric. Asymptotic and Aver-aging Methods were used to obtain stability conditionsfor several solutions of interest as H → H isthe Hubble parameter; suggesting that the asymptoticbehavior of the time-averaged model is independent ofthe coupling function and the geometry.Following analogous procedures, in reference [20]was studied the case (I) Bianchi III and OpenFLRW model proving that full time dependent equa-tions and their time-averaged versions have the samelate-time dynamics. Therefore, the simplest time-averaged system determines the future asymptoticof the full system. It was presented the gen-eral metric element: ds = − dt + (cid:2) e ( t ) (cid:3) − dr + (cid:2) e ( t ) (cid:3) − (cid:104) d ϑ + k − sin ( √ k ϑ ) d ζ (cid:105) where e , e and e = √ ke / sin ( √ k ϑ ) are functions of t [55].This family of metrics contains LRS Bianchi III ( k = − k =
0) and Kantowski-Sach ( k =+
1) which satisfy e ( t ) (cid:54) = e ( t ) . FLRW metrics areobtained when e ( t ) = e ( t ) = a ( t ) − (where a ( t ) isthe scale factor): Open (negative curvature) Universe iswhen k = −
1, Flat Universe ( k =
0) and Closed (pos-itive curvature) Universe ( k = +
1) [56–61]. Hence,according to the values of spatial curvature in (3) wehave divided the study in three cases: (I) Bianchi IIIand Open FLRW model, (II) Bianchi I and Flat FLRWmodel and (III) Kantowski- Sachs and Closed FLRW.For LRS Bianchi III is was proved in paper I[20] that the late-time attractors of the full andtime-averaged systems are: a matter dominated flatFLRW universe if 0 ≤ γ ≤ (mimicking de Sitter,quintessence or zero acceleration solutions), a matter-curvature scaling solution if < γ < ≤ γ ≤
2. For FLRW metric with k = − ≤ γ ≤ (mimicking de Sit-ter, quintessence or zero acceleration solutions) and theMilne solution if < γ <
2. In all metrics, the matterdominated flat FLRW universe represents quintessencefluid if 0 < γ < .This paper is devoted to the case (II) and case (III)will be studied in a companion paper [62].The paper is organized as follows. In section 2 weintroduce the model under study. In section 3 we applyAveraging Methods to analyze the periodic solutions ofa scalar field with self-interacting potentials within theclass of generalized harmonic potentials [17]. In par-ticular, section 3.1 is devoted to LRS Bianchi I modeland section 3.2 is devoted to flat FLRW metric. In sec-tion 4 we study the resulting time-averaged systems.In particular, section 4.1 is devoted to LRS Bianchi Imodels and section 4.2 is devoted the flat FLRW met-ric. Finally, in section 5 our main results are discussed.In appendix A is given the Proof of our main theorem,and in appendix B is presented numerical evidence sup-porting the results of section 3.
2. THE MODEL
It is well-known that there is an interesting hierar-chy of Bianchi models [63–66]. In particular, the LRSBianchi I model naturally appears as a boundary sub-set of the LRS Bianchi III model. The last one isan invariant boundary of the LRS Bianchi type VIIImodel as well. Additionally, LRS Bianchi type VIIIcan be viewed as an invariant boundary of the LRSBianchi type IX models [67–72]. Bianchi spacetimescontain many important cosmological models that havebeen used to study anisotropies of primordial Universeand its evolution towards the observed isotropy of thepresent epoch [73–76]. The list includes the standardFLRW model in the limit of the isotropization.In GR the Hubble parameter H is always mono-tonic for Bianchi I and the anisotropy decays on timefor H >
0. Therefore, isotropization occurs [77, 78].The exact solution of the field equations for an ex-ponential potential have been found in some particu-lar Bianchi spacetimes [79–81]. These exact solutionslead to isotropic homogeneous spacetimes as it wasfound in reference [82, 83]. An anisotropic solutionof special interest is the Kasner spacetime [84–91], es-sential for the description of the BKL singularity [92].The action integral of interest is S = (cid:90) d x (cid:112) | g | (cid:20) R + L φ + L m (cid:21) (1) expressed in a system of units in which 8 π G = c = (cid:125) =
1. In equation (1), R is the curvature scalar, φ is thescalar field with Lagrangian density L φ = − g µν ∇ µ φ ∇ ν φ − V ( φ ) . (2) ∇ α is the covariant derivative and the potential is(14). The metric element in hyper-spherical curvature-normalized coordinates [ t , r , ϑ , ζ ] is: ds = − dt + (cid:2) e ( t ) (cid:3) − dr + (cid:2) e ( t ) (cid:3) − (cid:2) d ϑ + ϑ d ζ (cid:3) , (3)where e , e and e = e / ϑ are functions of t which are the components of the frame vectors [55]: e = ∂ t , e = e ∂ r , e = e ∂ ϑ , e = e ∂ ζ . (4)This family of metrics contains the LRS Bianchi I( e ( t ) (cid:54) = e ( t ) ) and the flat FLRW metric when e ( t ) = e ( t ) = a ( t ) − (where a ( t ) is the scale fac-tor) [56]. It can be defined the average scale factor a ( t ) and the anisotropic parameter σ + through˙ aa = H : = − ddt ln (cid:2) e ( t )( e ( t )) (cid:3) , (5) σ + = ddt ln (cid:2) e ( t )( e ( t )) − (cid:3) . (6)Taking variation of (1) for the 1-parameter family ofmetrics (3) and after some algebra we obtain [56]:˙ H = − H − σ + − ( γ − ) ρ m −
13 ˙ φ + V ( φ ) , (7)˙ σ + = − H σ + , (8)3 H = σ + + ρ m +
12 ˙ φ + V ( φ ) , (9)where it is used for the matter component a barotropicEoS p m = ( γ − ) ρ m with p m the pressure of the mattercomponent, ρ m the energy density and γ is known asbarotropic index satisfying 0 ≤ γ ≤ σ + = σ / a . The term G ( a ) = σ / a which corresponds to anisotropies inBianchi I metric does not corresponds to a fluid com-ponent in the model. However, it can be interpreted asa stiff- matter in a flat FLRW background. The term σ + dilutes very fast with the expansion in Bianchi I metric,isotropizing if H > K =( e ( t )) , ( ) R =
0. Furthermore, the evolution ofGauss curvature is˙ K = − ( σ + + H ) K , (10)while the evolution for e is given by [55]:˙ e = − ( H − σ + ) e . (11)Finally, the evolution equation for matter is˙ ρ m = − γ H ρ m , (12)and the Klein-Gordon equation is¨ φ = − H ˙ φ − dV ( φ ) d φ . (13)Scalar fields are relevant in the physical descriptionof the Universe, particularly, in the inflationary sce-nario [93–97]. Chaotic inflation is a model of cosmicinflation which takes the potential term V ( φ ) = m φ φ ,the so-called harmonic potential [94–97].In this research we assume the generalized harmonicpotential plus cosine-like corrections with small phase V ( φ ) = µ (cid:20) b f (cid:18) − cos (cid:18) φ f (cid:19)(cid:19) + φ µ (cid:21) , b > . (14)inspired in the potential of [98, 99].Potential (14) can be expressed as V ( φ ) = µ φ + f (cid:0) ω − µ (cid:1) (cid:18) − cos (cid:18) φ f (cid:19)(cid:19) (15)by imposing the condition b µ + f µ − f ω = ω ∈ R satisfy-ing ω − µ > φ =
0, wehave V ( φ ) ∼ ω φ + O (cid:0) φ (cid:1) , as φ →
0. That is, ω can be related to the mass of the scalar field near itsglobal minimum. As φ → ± ∞ the cosine- correction isbounded. Then, V ( φ ) ∼ µ φ + O ( ) as φ → ± ∞ .
3. AVERAGING SCALAR FIELD COSMOLOGIES
The term averaging here is related to the approxi-mation of initial value problems in ordinary differen-tial equations which involves perturbations (chapter 11,[15]). In [25] it was studied the large-time behaviorfor solutions of a general class of systems in standardform: (cid:18) ˙ H ˙ x (cid:19) = H (cid:18) f ( t , x ) (cid:19) + H (cid:18) f [ ] ( t , x ) (cid:19) ,where H is the Hubble parameter and is positive strictlydecreasing in t and lim t → ∞ H ( t ) =
0. This result has ap-plications for spatially homogeneous cosmologies withoscillatory behaviour as H → ω inthe amplitude-phase transformation. In this section the Averaging Methods are applied forBianchi I metrics for the generalized harmonic poten-tial (14) minimally coupled to matter.In this case the field equations are¨ φ = − H ˙ φ − V (cid:48) ( φ ) , ˙ ρ m = − γ H ρ m , ˙ a = aH , ˙ H = − (cid:0) γρ m + ˙ φ (cid:1) − σ a , H = ρ m +
12 ˙ φ + V ( φ ) + σ a . (16)Defining Ω = (cid:115) ω φ + ˙ φ H , Σ = σ √ a H , Φ = t ω − tan − (cid:18) ωφ ˙ φ (cid:19) , (17)we obtain˙ Ω = − b γ f µ Ω H sin √ H Ω sin ( t ω − Φ ) f ω − b µ √ H cos ( t ω − Φ ) sin (cid:18) √ H Ω sin ( t ω − Φ ) f ω (cid:19) + H Ω (cid:18) − ( γ − ) Σ + γ − γ µ Ω ω (cid:19) + H cos ( t ω − Φ ) (cid:18) Ω (cid:18) γ (cid:18) µ ω − (cid:19) + (cid:19) − Ω (cid:19) + (cid:0) ω − µ (cid:1) Ω sin ( t ω − Φ ) ω , (18a)˙ Σ = − b γ f µ Σ H sin √ H Ω sin ( t ω − Φ ) f ω + H (cid:32) − ( γ − ) Σ (cid:0) Σ − (cid:1) − γ µ ΣΩ sin ( t ω − Φ ) ω (cid:33) − ( γ − ) Σ H Ω cos ( t ω − Φ ) , (18b)˙ Φ = − b µ sin ( t ω − Φ ) √ H Ω sin (cid:18) √ H Ω sin ( t ω − Φ ) f ω (cid:19) − H sin ( ( t ω − Φ ))+ (cid:0) ω − µ (cid:1) sin ( t ω − Φ ) ω , (18c)˙ H = − ( + q ) H , (18d)where the deceleration parameter is given by q = − − b γ f µ H sin √ H Ω sin ( t ω − Φ ) f ω + (cid:18) − ( γ − ) Σ + γ − γ µ Ω sin ( t ω − Φ ) ω (cid:19) − ( γ − ) Ω cos ( t ω − Φ ) . (19)With f = b µ ω − µ , we have1 = Ω + Σ + ρ m H + H ( ω − µ ) U ( φ ) , (20)where the function defined by U ( φ ) = f (cid:16) − cos (cid:16) φ f (cid:17)(cid:17) − φ , satisfies U (cid:48) ( ) = U (cid:48)(cid:48) ( ) = U (cid:48)(cid:48)(cid:48) ( ) = U ( iv ) ( ) = − f , which implies φ = U ( φ ) with critical value zero. Then, for ω > µ , Ω , Σ and Ω m : = ρ m H can be greater than 1.The fractional energy density of matter Ω m : = ρ m H is parameterized by the equation Ω m = − Σ − Ω + Ω (cid:18) − µ ω (cid:19) sin ( Φ − t ω )+ b µ H ( µ − ω ) sin √ H Ω ( µ − ω ) sin ( Φ − t ω ) b µ ω = − Σ − Ω + O ( H ) . (21)Setting f = b µ ω − µ >
0, it is obtained the series expan-sion near H = ˙ x = H f ( t , x ) + O ( H ) , ˙ H = − H (cid:32) γ (cid:16) − Σ − Ω (cid:17) + Σ + Ω cos ( Φ − t ω ) (cid:33) + O ( H ) , (22) f ( t , x ) = Ω (cid:32) γ (cid:0) − Σ − Ω (cid:1) + Σ + (cid:0) Ω − (cid:1) cos ( Φ − t ω ) (cid:33) Σ (cid:32) − γ (cid:0) Σ + Ω − (cid:1) + Σ + Ω cos ( Φ − t ω ) − (cid:33) − sin ( t ω − Φ ) . (23) Replacing ˙ x = H f ( t , x ) with x = ( Ω , Σ , Φ ) T and f ( t , x ) defined by (23) by ˙ y = H ¯ f ( y ) , with y = (cid:0) ¯ Ω , ¯ Σ , ¯ Φ (cid:1) T and¯ f defined by¯ f ( · ) : = L (cid:90) L f ( · , t ) dt , L = πω (24)we obtain the averaged system:˙¯ Ω = H ¯ Ω (cid:0) γ (cid:0) − ¯ Σ − ¯ Ω (cid:1) + Σ + ¯ Ω − (cid:1) , (25)˙¯ Σ = H ¯ Σ (cid:0) γ (cid:0) − ¯ Σ − ¯ Ω (cid:1) + Σ + ¯ Ω − (cid:1) , (26)˙¯ Φ = , (27)˙ H = − H (cid:0) γ (cid:0) − ¯ Σ − ¯ Ω (cid:1) + Σ + ¯ Ω (cid:1) . (28) Theorem 1
Let be defined the functions ¯ Ω , ¯ Σ , ¯ Φ and Hsatisfying the averaged equations (25) , (26) , (27) and (28) . Then, there exists three functions g , g and g such that if Ω = Ω + Hg ( H , Ω , Σ , Φ , t ) , Σ = Σ + Hg ( H , Ω , Σ , Φ , t ) , Φ = Φ + Hg ( H , Ω , Σ , Φ , t ) , (29) the functions Ω , Σ , Φ and ¯ Ω , ¯ Σ , ¯ Φ have the samelimit as t → ∞ . Setting Σ = Σ = are derived analo-gous results for flat FLRW model. Proof.
The proof is based on the following steps:1. Proving that H is a monotonic decreasing func-tion of t , if 0 < Ω + Σ <
1. Then, can be definedrecursively the bootstrapping sequences t = t ∗ H = H ( t ∗ ) , t n + = t n + H n H n + = H ( t n + ) , (30)such that lim n → ∞ H n = n → ∞ t n = ∞ .2. Finding g i such that it can be defined Ω , Σ , Φ through equation (29) and proved that g i , i = , . . . , t ∈ [ t n , t n + ] . By continu-ity of ¯ Ω , Ω , ¯ Σ , Σ , Φ , Ω k , can be defined constants L , M , M , M which are finite and independent of H .3. Using the expansions (29), defining ∆Ω = Ω − ¯ Ω , ∆Σ = Σ − ¯ Σ , ∆Φ = Φ − ¯ Φ and taking thesame initial conditions at t = t n : Ω ( t n ) = ¯ Ω ( t n ) = Ω n , Σ ( t n ) = ¯ Σ ( t n ) = Σ n , Φ ( t n ) = ¯ Φ ( t n ) = Φ n , < Ω n < , − < Σ n < . Using Gronwall’s Lemma 4, wehave for t ∈ [ t n , t n + ] : (cid:12)(cid:12)(cid:12) ∆Ω ( t ) (cid:12)(cid:12)(cid:12) ≤ M H n ( t − t n ) (cid:104) + LH n ( t − t n ) e LH n ( t − t n ) (cid:105) (cid:46) M ( t − t n ) H n + O ( H n ) ≤ M H n + O ( H n ) , (cid:12)(cid:12)(cid:12) ∆Σ ( t ) (cid:12)(cid:12)(cid:12) ≤ M H n ( t − t n ) (cid:104) + LH n ( t − t n ) e LH n ( t − t n ) (cid:105) (cid:46) M ( t − t n ) H n + O ( H n ) ≤ M H n + O ( H n ) , | ∆Φ ( t ) | ≤ M H n , due to t − t n ≤ t n + − t n = H n .
4. Finally, taking the limit as n → ∞ we obtain H n → H n → Ω , Σ , Φ and¯ Ω , ¯ Σ , ¯ Φ have the same limit as τ → ∞ .The complete proof is given in appendix A.Theorem 1 implies that Ω , Σ , Φ evolve according tothe averaged equations(25), (26), (27) as H → k = . In this case the field equations are¨ φ = − H ˙ φ − V (cid:48) ( φ ) , (31a)˙ ρ m = − γ H ρ m , (31b)˙ H = − (cid:0) γρ m + ˙ φ (cid:1) , (31c)3 H = ρ m +
12 ˙ φ + V ( φ ) , (31d)Defining Ω = (cid:115) ω φ + ˙ φ H , Φ = t ω − tan − (cid:18) ωφ ˙ φ (cid:19) , (32)we obtain the equations˙ Ω = − b γ f µ Ω H sin √ H Ω sin ( t ω − Φ ) f ω − b µ √ H cos ( t ω − Φ ) sin (cid:18) √ H Ω sin ( t ω − Φ ) f ω (cid:19) + H cos ( t ω − Φ ) (cid:18) Ω (cid:18) γ (cid:18) µ ω − (cid:19) + (cid:19) − Ω (cid:19) + H Ω (cid:18) γ − γ µ Ω ω (cid:19) + (cid:0) ω − µ (cid:1) Ω sin ( t ω − Φ ) ω , (33a)˙ Φ = − b µ √ H Ω sin ( t ω − Φ ) sin (cid:18) √ H Ω sin ( t ω − Φ ) f ω (cid:19) − H sin ( t ω − Φ ) cos ( t ω − Φ )+ (cid:0) ω − µ (cid:1) ω sin ( t ω − Φ ) , (33b) ˙ H = − ( + q ) H , (33c) q = − + γ − b γ f µ H sin √ H Ω sin ( t ω − Φ ) f ω − γ µ Ω sin ( t ω − Φ ) ω − ( γ − ) Ω cos ( t ω − Φ ) . (33d)For flat FLRW the fractional energy density of matter Ω m : = ρ m H is parameterized by Ω m = − Ω + Ω (cid:18) − µ ω (cid:19) sin ( Φ − t ω )+ b µ H ( µ − ω ) sin √ H Ω ( µ − ω ) sin ( Φ − t ω ) b µ ω = − Ω + O ( H ) . (34)Setting f = b µ ω − µ >
0, we obtain:˙ x = H f ( t , x ) + O ( H ) , x = ( Ω , Φ ) T , ˙ H = − H (cid:18) γ (cid:0) − Ω (cid:1) + Ω cos ( t ω − Φ ) (cid:19) + O ( H ) . (35) f ( t , x )= γ (cid:0) − Ω (cid:1) + Ω (cid:0) Ω − (cid:1) cos ( t ω − Φ ) − sin ( t ω − Φ ) . (36)Replacing ˙ x = H f ( t , x ) with x = ( Ω , Φ ) T and f ( t , x ) de-fined by (36) with ˙ y = H ¯ f ( y ) , where y = (cid:0) ¯ Ω , ¯ Φ (cid:1) T withthe time averaging (24), we obtain the following time-averaged system:˙¯ Ω = − H ¯ Ω ( γ − ) (cid:0) ¯ Ω − (cid:1) , (37)˙¯ Φ = . (38)Theorem 1 applies for Bianchi I and the invariant set Σ =
4. QUALITATIVE ANALYSIS OF AVERAGEDSYSTEMS
According to Theorem 1 for Bianchi I and flatFLRW the quantity H plays the role of time dependentperturbation parameters which controls the magnitudeof error between the solutions of the full and time-averaged problems. Thus, the oscillations are viewedas perturbations. In the time-averaged system the Ray-chaudhuri equation decouples. Therefore, the analysisof the system is reduced to study the corresponding av-eraged equations. In this case the averaged system is (25), (26) and(27). Introducing the new time variable τ through d fd τ = H d fdt , we obtain the guiding system: d ¯ Ω d τ =
32 ¯ Ω (cid:0) γ (cid:0) − ¯ Σ − ¯ Ω (cid:1) + Σ + Ω − (cid:1) , (39a) d ¯ Σ d τ =
32 ¯ Σ (cid:0) γ (cid:0) − ¯ Σ − ¯ Ω (cid:1) + Σ + ¯ Ω − (cid:1) . (39b)We have ¯ Σ + ¯ Ω + ¯ Ω m =
1. Therefore, from the con-dition ¯ Ω m ≥ (cid:8) ( ¯ Ω , ¯ Ω k ) ∈ R : ¯ Ω + ¯ Ω k ≤ , ¯ Ω k ≥ (cid:9) . (40)The 2D guiding system (39a)-(39b) has the followingequilibrium points:1. T : ( ¯ Ω , ¯ Σ ) = ( , − ) with eigenvalues (cid:8) , ( − γ ) (cid:9) .i) It is a source for 0 ≤ γ < , ii) It is nonhyperbolic for γ = . In the last case T becomes a line of points L parametrized by ¯ Σ ∗ given by ¯ Ω =
0, and ¯ Σ = ¯ Σ ∗ . The eigenvalues are { , } and the eigenvectorassociated with the zero eigenvalue ( , ) is par-allel to the tangent vector of the line at the givenpoint. This means that the line is a normally hy-perbolic set of equilibrium points and the stabil-ity is determined by the nonzero eigenvalue. Bythe previous statements, the line is unstable.The line element (3) becomes ds = − dt + ( H t + ) c dr + c (cid:2) d ϑ + ϑ d ζ (cid:3) . (41)Therefore, the corresponding solution can be ex-pressed in the form of the Taub-Kasner solution( p = , p = , p =
0) [100] (Sect 6.2.2 and p193, Eq. (9.6)).2. Q : ( ¯ Ω , ¯ Σ ) = ( , ) with eigenvalues (cid:8) , ( − γ ) (cid:9) .i) It is a source for 0 ≤ γ < . ii) It is nonhyperbolic for γ = . In this case, Q is included in the line L . The line element (3) becomes ds = − dt + c − ( H t + ) − dr + c − ( H t + ) / (cid:2) d ϑ + ϑ d ζ (cid:3) . (42) Therefore, the corresponding solution can be ex-pressed in the form of the non-flat LRS Kasner( p = − , p = , p = ) Bianchi I solution([100] Sect. 6.2.2 and Sect. 9.1.1 (2)).3. F : ( ¯ Ω , ¯ Σ ) = ( , ) with eigenvalues (cid:110) ( γ − ) , ( γ − ) (cid:111) .i) It is a sink for 0 ≤ γ < < γ < γ = , . If γ = , F is included in the line L . If γ = , F isincluded in the line L . The line element (3) becomes ds = − dt + a (cid:18) γ H t + (cid:19) γ dr a (cid:18) γ H t + (cid:19) γ (cid:2) d ϑ + ϑ d ζ (cid:3) . (43)The corresponding solution is the flat matterdominated FLRW Universe with ¯ Ω m = F : ( ¯ Ω , ¯ Σ ) = ( , ) , with eigenvalues (cid:8) − , − ( γ − ) (cid:9) .i) It is a saddle 0 ≤ γ < < γ ≤ γ = F becomes a line of points L parametrized by ¯ Ω ∗ given by ¯ Ω = ¯ Ω ∗ and ¯ Σ = . The eigenvalues are {− , } . As before, thisline is a normally hyperbolic set of equilibriumpoints and the stability is determined by thenonzero eigenvalue. By the previous statements,the line is stable when γ =
1. Evaluating theKlein-Gordon equation (18d) at the equilibriumpoint F we have˙ H = b γ µ ω − µ sin √ H ( ω − µ ) sin ( Φ − t ω ) b µ ω + H (cid:16) − γ + γ µ sin ( Φ − t ω ) ω + ( γ − ) cos ( Φ − t ω ) (cid:17) . (44)Therefore,˙ H ∼ − H cos ( t ω − Φ ) , (45) F QFT - - Σ Ω γ = (a) F QFT - - Σ Ω γ = (b) F QFT - - Σ Ω γ = / (c) F QFT - - Σ Ω γ = (d) Figure 1: Phase plane for system (39a), (39b) for different choices of γ . For γ = L corresponds to the set of equilibrium points ¯ Σ = , ¯ Ω = ¯ Ω ∗ . For γ = L corresponds tothe set of equilibrium points ¯ Σ = ¯ Σ ∗ , ¯ Ω = . for large t . In average, Φ is a constant, setting Φ = H ( t ) = H ω H t ω + H sin ( t ω ) + ω , (46)where H is the current value of H ( t ) . Finally, H ( t ) ∼ t for large t . Gauss equation (10) andevolution equation (11) become˙ e = − e t , ˙ K = − K t , (47) with general solution˙ e ( t ) = c t / , K ( t ) = c t / . (48)That is, the line element (3) becomes ds = − dt + c − t / dr + c − t / (cid:2) d ϑ + ϑ d ζ (cid:3) . (49)Hence for large t and c = c the equilibriumpoint can be associated with the flat FriedmanEinstein-de-Sitter solution ([100], Sec 9.1.1 (1))with γ = γ . The guiding system (39a), (39b) when γ = d ¯ Ω d τ =
32 ¯ Σ ¯ Ω (50) d ¯ Σ d τ = −
32 ¯ Σ (cid:0) − ¯ Σ (cid:1) . (51)For the above system we obtain dd τ ln (cid:18) ¯ Σ ¯ Ω (cid:19) = − = ⇒ ¯ Σ ¯ Σ = ¯ Ω ¯ Ω e − τ (52)The orbits of system (50) and (51) are given by¯ Σ = − (cid:0) − ¯ Σ (cid:1) ¯ Ω ¯ Ω . (53) The results from the linear stability analysis com-bined with Theorem 1 lead to:
Theorem 2
The late-time attractors of full system (39) and time-averaged system (B2) for Bianchi I line ele-ment are:(i) The flat matter dominated FLRW Universe F with the line elementds = − dt + a (cid:18) γ H t + (cid:19) γ dr a (cid:18) γ H t + (cid:19) γ (cid:2) d ϑ + ϑ d ζ (cid:3) (54) if < γ < . F represents a quintessence fluidfor < γ < or a zero-acceleration model for γ = . In the limit γ = we have a ( t ) = a (cid:16) γ H t + (cid:17) γ → a e H t , i.e., the de Sitter so-lution.(ii) The scalar field dominated solution F with theline element (3) ds = − dt + c − t / dr + c − t / (cid:2) d ϑ + ϑ d ζ (cid:3) . (55) Hence for large t and c = c the equilibriumpoint can be associated with the flat FriedmanEinstein-de-Sitter solution. k = . For flat Universe k = γ (cid:54) =
1, we obtain thefollowing time-averaged system:˙¯ Ω = − H ¯ Ω (cid:0) − ¯ Ω (cid:1) ( − γ ) , ˙¯ Φ = . (56)Introducing the new time variable τ through d fd τ = H d fdt ,we obtain the guiding system d ¯ Ω d τ = −
32 ¯ Ω (cid:0) − ¯ Ω (cid:1) ( − γ ) , (57)with solution¯ Ω ( τ ) = Ω e γτ (cid:113) Ω e γτ + e τ (cid:0) − Ω (cid:1) , (58)where ¯ Ω ( ) = Ω .Equation (57) has the following equilibrium points1. F with eigenvalue − ( − γ ) . It is a sink for0 < γ < < γ < F with eigenvalue 3 ( − γ ) . It is a source for0 < γ < < γ < . In Figure 2 is presented a geometric analysis ofequation (57) for different choices of γ . For γ =
1, the time-averaged system truncated at or-der O ( H ) is d ¯ Ω d τ = H ¯ Ω (cid:0) ω − µ (cid:1) b µ ω (59a) d ¯ Φ d τ = H ¯ Ω (cid:0) ω − µ (cid:1) b µ ω − H ¯ Ω (cid:0) ω − µ (cid:1) b µ ω , (59b) dHd τ = − H − H ¯ Ω (cid:0) ω − µ (cid:1) b µ ω . (59c)From (59) we obtain H (cid:48) ( ¯ Ω ) = b µ ω Ω H ( ¯ Ω ) ( µ − ω ) − H ( ¯ Ω ) ¯ Ω . (60)0Assuming H ( ¯ Ω ) = H we obtain the solution H ( ¯ Ω ) = (cid:114) b µ ω ( ω − µ ) (cid:16) − ¯ Ω ¯ Ω (cid:17) + H Ω ¯ Ω ¯ Ω . (61)Then, equations (59) can be expressed as d ¯ Ω d η = Ω − Ω (cid:16) b µ ω + H ¯ Ω (cid:0) µ − ω (cid:1) (cid:17) b µ ω ¯ Ω , (62a) d Φ d τ = (cid:0) ω − µ (cid:1) (cid:118)(cid:117)(cid:117)(cid:116) b µ ω (cid:18) − ¯ Ω Ω (cid:19) ( ω − µ ) + H ¯ Ω ¯ Ω b µ ω − (cid:0) ω − µ (cid:1) b µ ω (cid:18) − ¯ Ω Ω (cid:19) ( ω − µ ) + H ¯ Ω ¯ Ω / b µ ω ¯ Ω . (62b)By integration we obtain ¯ Ω ( τ ) = b µ e τ / ω ¯ Ω (cid:113) b µ e τ ω + H ( − e τ ) ¯ Ω ( ω − µ ) , satisfying lim τ →− ∞ Ω ( τ ) =
0, and lim τ → + ∞ Ω ( τ ) = b µ ω ¯ Ω (cid:113) b µ ω + H ¯ Ω ( µ − ω ) .Furthermore, H ( τ ) = H e − τ (cid:113) b µ e τ ω + H ( e τ − ) ¯ Ω ( µ − ω ) b µ ω , satisfies lim τ →− ∞ H ( τ ) = ∞ and lim τ → + ∞ H ( τ ) =
0. Finally, Φ (cid:48) ( τ ) = − H e − τ ¯ Ω (cid:0) ω − µ (cid:1) (cid:16) H ¯ Ω (cid:0) ω − µ (cid:1) − b µ e τ ω (cid:17) b µ ω (cid:113) b µ e τ ω + H ( e τ − ) ¯ Ω ( µ − ω ) , is integrable, giving Φ ( τ ) . The result is quite long sowe omit the explicit expression. The results from the linear stability analysis com-bined with Theorem 1 (for Σ =
0) lead to:
Theorem 3
The late-time attractors of the full system (33) and the averaged system (57) are:(i) The flat matter dominated FLRW Universe F with line elementds = − dt + a (cid:18) γ H t + (cid:19) γ dr a (cid:18) γ H t + (cid:19) γ (cid:2) d ϑ + ϑ d ζ (cid:3) (63) for ≤ γ < . F represents a quintessencefluid if < γ < or a zero-acceleration modelfor γ = . In the limit γ = we have a ( t ) = a (cid:16) γ H t + (cid:17) γ → a e H t , i.e., the de Sitter so-lution.(ii) The scalar field dominated solution F with lineelementds = − dt + c − t / dr + c − t / (cid:2) d ϑ + ϑ d ζ (cid:3) , (64) for < γ < . Hence for large t and c = c theequilibrium point can be associated with the flatFriedman Einstein-de-Sitter solution. F F - - - - - - Ω Ω ′ γ = (a) F F Ω Ω ′ γ = / (b) Figure 2: Analysis for equation (57) for different choices of γ .
5. CONCLUSIONS
This is the second paper which follows the pro-gram “Averaging Generalized Scalar Field Cosmolo-gies” which started in reference [20]. This programconsists in using Asymptotic Methods and AveragingTheory to obtain relevant information about the solu-tion’s space of scalar field cosmologies in the presenceof a background matter with a barotropic Equation ofState (EoS) with barotropic index γ minimally coupledto a scalar field with generalized harmonic potential V ( φ ) = µ (cid:104) b f (cid:16) − cos (cid:16) φ f (cid:17)(cid:17) + φ µ (cid:105) , b > ≤ γ ≤ (mimicking de Sitter, quintessence orzero acceleration solutions), a matter-curvature scalingsolution if < γ < ≤ γ ≤
2. For FLRW metric with k = − ≤ γ ≤ (mimicking de Sitter,quintessence or zero acceleration solutions) and theMilne solution if < γ <
2. In all metrics, the matter dominated flat FLRW universe represents quintessencefluid if 0 < γ < .For LRS Bianchi I and flat FLRW metrics as wellas for LRS Bianchi III and Open FLRW, we can useTaylor expansion with respect to H near H =
0. Hence,the resulting system can be expressed in the standardform [25] (cid:18) ˙ H ˙ x (cid:19) = H (cid:18) f ( t , x ) (cid:19) + H (cid:18) f [ ] ( t , x ) (cid:19) , where H is the Hubble parameter and is positivestrictly decreasing in t and lim t → ∞ H ( t ) = F = (cid:18) f ( t , x ) (cid:19) is the first order term and F [ ] = (cid:18) f [ ] ( t , x ) (cid:19) is the second-order remainder of the Tay-lor series of the left hand side of (cid:18) ˙ H ˙ x (cid:19) = f ( t , x , H ) at H = Ω , Σ , and Φ evolve according to thetime-averaged equations (25), (26) and (27) as H → H is a time de-pendent perturbation parameter controlling the magni-tude of error between the solutions of the full and thetime-averaged problems. The analysis of the systemis therefore reduced to study the corresponding time-averaged equations.In this regard, we have formulated theorems 2 and3 concerning to the late-time behavior of our model,whose proofs are based on the previous theorem andthe linear stability analysis. For LRS Bianchi I the late-time attractors of full system (39) and averaged system(B2) are:(i) The flat matter dominated FLRW Universe F with line element ds = − dt + a (cid:18) γ H t + (cid:19) γ dr a (cid:18) γ H t + (cid:19) γ (cid:2) d ϑ + ϑ d ζ (cid:3) (65)if 0 < γ < F represents a quintessence fluid ora zero-acceleration model for γ = . In the limit γ = a ( t ) = a (cid:16) γ H t + (cid:17) γ → a e H t ,2i.e., the de Sitter solution.(ii) The scalar field dominated solution F with lineelement ds = − dt + c − t / dr + c − t / (cid:2) d ϑ + ϑ d ζ (cid:3) (66)for 1 < γ <
2. Hence for large t and c = c theequilibrium point can be associated with the flatFriedman Einstein-de-Sitter solution.For flat FLRW metric the late-time attractors of thefull system (33) and the averaged system are:(i) The flat matter dominated FLRW Universe F with line element ds = − dt + a (cid:18) γ H t + (cid:19) γ dr a (cid:18) γ H t + (cid:19) γ (cid:2) d ϑ + ϑ d ζ (cid:3) . (67)for 0 ≤ γ < F represents a quintessence fluidor a zero-acceleration model for γ = . In thelimit γ = a ( t ) = a (cid:16) γ H t + (cid:17) γ → a e H t , i.e., the de Sitter solution.(ii) The scalar field dominated solution F with lineelement ds = − dt + c − t / dr + c − t / (cid:2) d ϑ + ϑ d ζ (cid:3) , (68)for 1 < γ <
2. Hence for large t and c = c theequilibrium point can be associated with the flatFriedman Einstein-de-Sitter solution.Our analytical results were strongly supported bynumerics in appendix B. It is interesting to note thatfor LRS Bianchi I and flat FLRW cases when the back-ground fluid corresponds to a cosmological constant, H tends asymptotically to constant values dependingon the initial conditions which is consistent to de Sitterexpansion (see figures 3(a) and 7(a)). For dust γ = Ω ( τ ) tends to a constant and H ( τ ) tends to zero as τ → + ∞ . Observe in figure 7(b)that as H → Ω lying on the orange line(solution of the time-averaged system at fourth order)give an upper bound to the value Ω of the original sys-tem. Therefore, by controlling the error of the averagedhigher order system one can also control the error of theoriginal system.Summarizing, in this paper (paper II), we have ex-amined in more detail the sub-case Ω k = ≤ γ <
1, or an equi-librium solution for 1 < γ ≤ t and c = c it can be associated with the flat FriedmanEinstein-de-Sitter solution. For flat FLRW metric, thelate-time attractors of the full system and the time-averaged system are: the flat matter dominated FLRW(mimicking de Sitter, quintessence or zero accelera-tion solutions) for 0 < γ < < γ <
2. In all met-rics, the matter dominated flat FLRW Universe repre-sents quintessence fluid if 0 < γ < .We have shown that asymptotic methods and averag-ing theory are powerful tools to investigate scalar fieldcosmologies with generalized harmonic potential. Acknowledgements
This research was funded by Agencia Nacional deInvestigaci´on y Desarrollo- ANID through the pro-gram FONDECYT Iniciaci´on grant no. 11180126 andby Vicerrector´ıa de Investigaci´on y Desarrollo Tec-nol´ogico at Universidad Cat´olica del Norte. Ellen delos M. Fern´andez Flores is acknowledged for proof-reading this manuscript and for improving the English.
Appendix A: Proof of Theorem 1
Lemma 4 (Gronwall’s Lemma) Differentialform .i) Let be η ( t ) a nonnegative function, abso-lutely continuous over [ , T ] , which satis-fies almost everywhere the differential in-equality η (cid:48) ( t ) ≤ φ ( t ) η ( t ) + ψ ( t ) , (A1) where φ ( t ) and ψ ( t ) are nonnegativesummable functions over [ , T ] . Then, η ( t ) ≤ e (cid:82) t φ ( s ) ds (cid:20) η ( ) + (cid:90) t ψ ( s ) ds (cid:21) , (A2) is satisfied for all ≤ t ≤ T .ii) In particular, if η (cid:48) ( t ) ≤ φ ( t ) η ( t ) , t ∈ [ , T ] , η ( ) = , (A3) then η ≡ , t ∈ [ , T ] .2. Integral form i) Let be ξ ( t ) a nonnegative function,summable over [ , T ] which satisfies almost everywhere the integral inequality ξ ( t ) ≤ C (cid:90) t ξ ( s ) ds + C , C , C ≥ . (A4) Then ξ ( t ) ≤ C ( + C te C t ) , (A5) almost everywhere for t in ≤ t ≤ T . ii) In particular, if ξ ( t ) ≤ C (cid:90) t ξ ( s ) ds , C ≥ , (A6) almost everywhere for t in ≤ t ≤ T , then ξ ≡ , almost everywhere for t in ≤ t ≤ T .
Theorem 1 states:
Let be defined the functions ¯ Ω , ¯ Σ , ¯ Φ , and H, which satisfy the time-averaged equations (25) , (26) , (27) and (28) . Then, there exists three functions g , g and g such that if Ω = Ω + Hg ( H , Ω , Σ , Φ , t ) , Σ = Σ + Hg ( H , Ω , Σ , Φ , t ) , Φ = Φ + Hg ( H , Ω , Σ , Φ , t ) , (A7) the functions Ω , Σ , Φ and ¯ Ω , ¯ Σ , ¯ Φ have the same limit as t → ∞ . Setting Σ = Σ = are derived analogousresults for flat FLRW model. Proof . From the equation (22) it follows that H is a monotonic decreasing function of t , if 0 < Ω + Σ < t = t ∗ H = H ( t ∗ ) , t n + = t n + H n H n + = H ( t n + ) , (A8)such that lim n → ∞ H n = n → ∞ t n = ∞ .Giving the expansions (A7), we obtain˙ Ω = H Ω (cid:32) γ (cid:0) − Σ − Ω (cid:1) + Σ + (cid:0) Ω − (cid:1) cos ( Φ − t ω ) (cid:33) − ∂ g ∂ t H + O (cid:0) H (cid:1) , ˙ Σ = H Σ (cid:32) − ( γ − ) (cid:0) Σ − (cid:1) − ( γ − ) Ω + Ω cos ( ( Φ − t ω )) (cid:33) − ∂ g ∂ t H + O (cid:0) H (cid:1) , ˙ Φ = − ( ( Φ − t ω ) sin ( Φ − t ω )) H − ∂ g ∂ t H + O (cid:0) H (cid:1) . Let be defined ∆Ω = Ω − ¯ Ω , ∆Σ = Σ − ¯ Σ , ∆Φ = Φ − ¯ Φ and take the same initial conditions at t = t n , suchthat Ω ( t n ) = ¯ Ω ( t n ) = Ω n , Σ ( t n ) = ¯ Σ ( t n ) = Σ n , Φ ( t n ) = ¯ Φ ( t n ) = Φ n , < Ω n < , − < Σ n < ∂ g ∂ t = (cid:0) − ( γ − ) (cid:0) Σ − (cid:1) Ω − ( γ − ) Ω (cid:1) + Ω (cid:0) Ω − (cid:1) cos ( ( Φ − t ω ))+ Ω ¯ Ω + ¯ Σ (cid:18) ( γ − ) ¯ Ω + Ω (cid:19) + ( γ − ) ¯ Ω + (cid:18) − γ (cid:19) ¯ Ω , (A9a) ∂ g ∂ t = Σ (cid:0) γ (cid:0) − Σ − Ω (cid:1) + Σ + Ω − (cid:1) + Σ Ω cos ( ( Φ − t ω ))+ ¯ Σ (cid:18) ( γ − ) ¯ Ω − ( γ − ) (cid:19) + ( γ − ) ¯ Σ + Σ ¯ Σ + Σ ¯ Ω , (A9b) ∂ g ∂ t = −
32 sin ( ( Φ − t ω )) . (A9c)4The explicit expressions for g i are obtained by integration of (A9) as follows: g ( H , Ω , Σ , Φ , t ) = t Ω (cid:0) − γ (cid:0) Σ + Ω − (cid:1) + Σ + Ω − (cid:1) + Ω (cid:0) − Ω (cid:1) sin ( ( Φ − t ω )) ω + (cid:90) (cid:34) ¯ Σ ( t ) (cid:0) ( γ − ) ¯ Ω ( t ) + Ω (cid:1) + ¯ Ω ( t ) (cid:0) ¯ Ω ( t ) (cid:0) ( γ − ) ¯ Ω ( t ) + Ω (cid:1) − γ + (cid:1) (cid:35) dt + C ( Σ , Ω , Φ ) , (A10) g ( H , Ω , Σ , Φ , t ) = Σ t (cid:0) − γ (cid:0) Σ + Ω − (cid:1) + Σ + Ω − (cid:1) − Σ Ω sin ( ( Φ − t ω )) ω + (cid:90) (cid:34) ¯ Ω ( t ) (cid:0) ( γ − ) ¯ Σ ( t ) + Σ (cid:1) + ¯ Σ ( t ) (cid:0) ¯ Σ ( t ) (cid:0) ( γ − ) ¯ Σ ( t ) + Σ (cid:1) − γ + (cid:1) (cid:35) dt + C ( Σ , Ω , Φ ) , (A11) g ( H , Ω , Σ , Φ , t ) = ( Φ − t ω ) ω + C ( Σ , Ω , Φ ) , (A12)where C i are integration functions which are independent of H because the terms ∂ g i ∂ H ˙ H in the procedure contributewith O ( H ) -factors. The terms unbarred inside the integral are taken as constants during the integration.Then, we obtain ˙ ∆Ω = − H ∆Ω ( − ¯ Ω − ¯ Σ ) + O ( H ) , ˙ ∆Σ = − H ∆Σ ( − ¯ Ω − ¯ Σ ) + O (cid:0) H (cid:1) , ˙ ∆Φ = H ω (cid:34) (cid:0) µ − ω (cid:1) ¯ Ω b µ + ω cos ( Φ − t ω ) (cid:16) γ (cid:0) − Σ − Ω (cid:1) + Σ + (cid:0) Ω + (cid:1) cos ( ( Φ − t ω )) + Ω (cid:17)(cid:35) + O (cid:0) H (cid:1) . Let t ∈ [ t n , t n + ] such that t − t n ≤ t n + − t n = H n , we have: | ∆Ω ( t ) | = | Ω ( t ) − ¯ Ω ( t ) | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:90) tt n (cid:16) ˙ Ω ( s ) − ˙¯ Ω ( s ) (cid:17) ds (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:90) tt n (cid:20) H (cid:0) − ¯ Ω − ¯ Σ (cid:1) ∆Ω ( s ) + O ( H ) (cid:21) ds (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:90) tt n (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:0) − ¯ Ω − ¯ Σ (cid:1) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:124) (cid:123)(cid:122) (cid:125) ≤ L H (cid:124)(cid:123)(cid:122)(cid:125) ≤ H n | ∆Ω ( s ) | ds + M H n ( t − t n ) + O ( H n ) ≤ LH n (cid:90) tt n | ∆Ω ( s ) | ds + M H n ( t − t n ) + O ( H n ) . | ∆Σ ( t ) | = | Σ ( t ) − ¯ Σ ( t ) | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:90) tt n (cid:16) ˙ Σ ( s ) − ˙¯ Σ ( s ) (cid:17) ds (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:90) tt n (cid:20) H (cid:0) − ¯ Ω − ¯ Σ (cid:1) ∆Σ + O ( H ) (cid:21) ds (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:90) tt n (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:0) − ¯ Ω − ¯ Σ (cid:1) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:124) (cid:123)(cid:122) (cid:125) ≤ L H (cid:124)(cid:123)(cid:122)(cid:125) ≤ H n | ∆Σ ( s ) | ds + M H n ( t − t n ) + O ( H n ) ≤ LH n (cid:90) tt n | ∆Σ ( s ) | ds + M H n ( t − t n ) + O ( H n ) , whereby choosing g such that | ∂ g ∂ Σ | and | ∂ g ∂ Ω | and g such that | ∂ g ∂ Σ | and | ∂ g ∂ Ω | are bounded on t ∈ [ t n , t n + ] and5by continuity of ¯ Ω , ¯ Σ , Ω , Σ , Φ the following holds: L = max t ∈ [ t n , t n + ] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:0) − ¯ Ω − ¯ Σ (cid:1) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) < ∞ , M = max t ∈ [ t n , t n + ] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:32) ¯ Ω (cid:0) ( γ − ) ¯ Σ + Σ (cid:1) + (cid:0) ¯ Σ − (cid:1) (cid:0) ( γ − ) ¯ Σ + Σ (cid:1) (cid:33) ∂ g ∂ Σ + (cid:32) ¯ Σ (cid:0) ( γ − ) ¯ Ω + Ω (cid:1) + (cid:0) ¯ Ω − (cid:1) (cid:0) ( γ − ) ¯ Ω + Ω (cid:1) (cid:33) ∂ g ∂ Ω + (cid:32) − γ (cid:0) Σ + Ω − (cid:1) + Σ + (cid:0) Ω − (cid:1) cos ( ( Φ − t ω )) + Ω − (cid:33) g − ( γ − ) Σ Ω g − Ω (cid:0) Ω − (cid:1) sin ( Φ − t ω ) cos ( Φ − t ω ) ω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) < ∞ . M = max t ∈ [ t n , t n + ] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:32) ¯ Ω (cid:0) ( γ − ) ¯ Σ + Σ (cid:1) + (cid:0) ¯ Σ − (cid:1) (cid:0) ( γ − ) ¯ Σ + Σ (cid:1) (cid:33) ∂ g ∂ Σ + (cid:32) ¯ Σ (cid:0) ( γ − ) ¯ Ω + Ω (cid:1) + (cid:0) ¯ Ω − (cid:1) (cid:0) ( γ − ) ¯ Ω + Ω (cid:1) (cid:33) ∂ g ∂ Ω + Σ Ω ( − γ + cos ( ( Φ − t ω )) + ) g + (cid:32) − γ (cid:0) Σ + Ω − (cid:1) + Σ + Ω cos ( ( Φ − t ω )) + Ω − (cid:33) g − Σ Ω sin ( Φ − t ω ) cos ( Φ − t ω ) ω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) < ∞ . Using Gronwall’s Lemma 4, we have for t ∈ [ t n , t n + ] : (cid:12)(cid:12)(cid:12) ∆Ω ( t ) (cid:12)(cid:12)(cid:12) ≤ M H n ( t − t n ) (cid:104) + LH n ( t − t n ) e LH n ( t − t n ) (cid:105) (cid:46) M ( t − t n ) H n + O ( H n ) ≤ M H n + O ( H n ) . (cid:12)(cid:12)(cid:12) ∆Σ ( t ) (cid:12)(cid:12)(cid:12) ≤ M H n ( t − t n ) (cid:104) + LH n ( t − t n ) e LH n ( t − t n ) (cid:105) (cid:46) M ( t − t n ) H n + O ( H n ) ≤ M H n + O ( H n ) , due to t − t n ≤ t n + − t n = H n .Furthermore, we have | ∆Φ ( t ) | = | Φ ( t ) − ¯ Φ ( t ) | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:90) tt n (cid:16) ˙ Φ ( s ) − ˙¯ Φ ( s ) (cid:17) ds (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:90) tt n H ω (cid:34) (cid:0) µ − ω (cid:1) ¯ Ω b µ + ω cos ( Φ − s ω ) (cid:16) γ (cid:0) − Σ − Ω (cid:1) + Σ + (cid:0) Ω + (cid:1) cos ( ( Φ − s ω )) + Ω (cid:17)(cid:35) ds (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + O ( H ) ≤ M H n ( t − t n ) + O ( H n ) , where by choosing C such that | ∂ C ∂ Σ | and | ∂ C ∂ Ω | are bounded on t ∈ [ t n , t n + ] , say, C ≡ Ω , Ω , Σ , Φ the following holds: M = max t ∈ [ t n , t n + ] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ω (cid:40) (cid:0) µ − ω (cid:1) ¯ Ω b µ + ω cos ( Φ − t ω ) (cid:16) γ (cid:0) − Σ − Ω (cid:1) + Σ + (cid:0) Ω + (cid:1) cos ( ( Φ − t ω )) + Ω (cid:17)(cid:41)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) < ∞ . Hence, | ∆Φ ( t ) | ≤ M H n , due to t − t n ≤ t n + − t n = H n . Finally, taking the limitas n → ∞ we obtain H n →
0. Then, as H n →
0, itfollows the functions Ω , Σ , Φ and ¯ Ω , ¯ Σ , ¯ Φ have thesame limit as τ → ∞ . (cid:3) Appendix B: Numerical simulation
In this section we present the numerical evidencethat support the main theorem presented in the section3 by solving numerically the full and averaged systemsobtained for each metric, namely LRS Bianchi I andflat FLRW.For this purpose, it was elaborated an algorithmin the programming language
Python where the sys-tems of differential equations were solved using the solve ivp code provided by the
SciPy open-source
Python -based ecosystem. The integration method usedwas
Radau that is an implicit Runge-Kutta method ofthe Radau IIa family of order 5 with a relative and ab-solute tolerances of 10 − and 10 − , respectively. Allsystems of differential equations were integrated withrespect to τ instead of t , with an integration range of − ≤ τ ≤
10 for the original systems and − ≤ τ ≤
100 for the averaged systems. All of them partitionedin 10000 data points. Furthermore, each full and time-averaged systems were solved considering only onematter component. These are cosmological constant( γ = γ = γ = /
3) and stiff fluid ( γ = Ω = Ω m ≡ Ω m ≡
0. Finally we have considered the fixed constants µ = √ / b = √ / ω = √
2, that leads to a valueof f = b µ ω − µ = /
10 and fulfills the condition f ≥ V ( φ ) = φ + ( − cos ( φ )) . (B1)
1. LRS Bianchi I model
For the Bianchi I metric we integrate:1. The full system given by (18).2. The time-averaged system: ∂ τ ¯ Ω = ¯ Ω (cid:0) γ (cid:0) − ¯ Σ − ¯ Ω (cid:1) + Σ + Ω − (cid:1) , ∂ τ ¯ Σ = ¯ Σ (cid:0) γ (cid:0) − ¯ Σ − ¯ Ω (cid:1) + Σ + ¯ Ω − (cid:1) , ∂ τ ¯ Φ = , ∂ τ H = − H (cid:0) γ (cid:0) − ¯ Σ − ¯ Ω (cid:1) + Σ + ¯ Ω (cid:1) , ∂ τ t = / H , (B2)where as initial conditions we use the seven data setpresented in the Table I for a better comparison of bothsystems.In Figures 3(a)-6(b) are presented projections ofsome solutions of the full system (18) and time-averaged system (B2) in the ( Σ , H , Ω ) space with theirrespective projection when H = γ = γ = γ = / γ = γ = H tends asymptotically to a constant,which is consistent to de Sitter expansion.
2. Flat FLRW model
For flat FLRW model we integrate1. The full system given by (33).7 H i ii iii ivvvivii (a) Projections in the space ( Σ , H , Ω ) . The surface is given by the constraint Ω = − Σ . i ii iiiiv-v vivii i ii iiiiv-v vivii i ii iiiiv-v vivii i ii iiiiv-v vivii i ii iiiiv-v vivii i ii iiiiv-v vivii i ii iiiiv-v vivii 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.160.600.650.700.750.800.850.900.951.00 ivii ivii ivii ivii ivii ivii ivii0.3 0.4 0.5 0.6 0.7 0.80.000.020.040.060.080.100.120.140.16 ii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v vi 0.42 0.44 0.46 0.48 0.50 0.52 0.54 0.56 0.58 0.600.00750.01000.01250.01500.01750.02000.02250.0250 iv vviiv vviiv vviiv vviiv vviiv vviiv vvi (b) Projection in the space ( Σ , Ω ) . The black line represent the constraint Ω = − Σ . Figure 3: Some solutions of the full system (18) (blue) and time-averaged system (B2) (orange) for the Bianchi I metric when γ =
0. We have used for both systems the initial data sets presented in the Table I. H i ii iii ivvvivii (a) Projections in the space ( Σ , H , Ω ) . The surface is given by the constraint Ω = − Σ . i ii iiiiv-v vivii i ii iiiiv-v vivii i ii iiiiv-v vivii i ii iiiiv-v vivii i ii iiiiv-v vivii i ii iiiiv-v vivii i ii iiiiv-v vivii 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.160.600.650.700.750.800.850.900.951.00 ivii ivii ivii ivii ivii ivii ivii0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.700.000.020.040.060.080.100.120.14 ii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v vi 0.42 0.44 0.46 0.48 0.50 0.52 0.54 0.560.0080.0100.0120.0140.0160.0180.0200.0220.024 ivv viivv viivv viivv viivv viivv viivv vi (b) Projection in the space ( Σ , Ω ) . The black line represent the constraint Ω = − Σ . Figure 4: Some solutions of the full system (18) (blue) and time-averaged system (B2) (orange) for the Bianchi I metric when γ =
1. The line denoted by red diamonds corresponds to the attracting line of equilibrium points. We have used for both systemsthe initial data sets presented in the Table I. H i ii iii ivvvivii (a) Projections in the space ( Σ , H , Ω ) . The surface is given by the constraint Ω = − Σ . i ii iiiiv-v vivii i ii iiiiv-v vivii i ii iiiiv-v vivii i ii iiiiv-v vivii i ii iiiiv-v vivii i ii iiiiv-v vivii i ii iiiiv-v vivii 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.160.600.650.700.750.800.850.900.951.00 ivii ivii ivii ivii ivii ivii ivii0.3 0.4 0.5 0.6 0.7 0.80.000.020.040.060.080.100.120.140.16 ii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v vi 0.42 0.44 0.46 0.48 0.50 0.52 0.54 0.56 0.58 0.600.00750.01000.01250.01500.01750.02000.02250.02500.0275 iv vviiv vviiv vviiv vviiv vviiv vviiv vvi (b) Projection in the space ( Σ , Ω ) . The black line represent the constraint Ω = − Σ . Figure 5: Some solutions of the full system (18) (blue) and time-averaged system (B2) (orange) for the Bianchi I metric when γ = /
3. We have used for both systems the initial data sets presented in the Table I. H i ii iiiiv v vivii (a) Projections in the space ( Σ , H , Ω ) . The surface is given by the constraint Ω = − Σ . i ii iiiiv-v vivii i ii iiiiv-v vivii i ii iiiiv-v vivii i ii iiiiv-v vivii i ii iiiiv-v vivii i ii iiiiv-v vivii i ii iiiiv-v vivii 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.160.600.650.700.750.800.850.900.951.00 ivii ivii ivii ivii ivii ivii ivii0.35 0.40 0.45 0.50 0.55 0.60 0.650.000.020.040.060.080.100.120.140.16 ii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v vi 0.470 0.475 0.480 0.485 0.490 0.495 0.500 0.505 0.5100.0080.0100.0120.0140.0160.0180.0200.0220.0240.026 iv-v viiv-v viiv-v viiv-v viiv-v viiv-v viiv-v vi (b) Projection in the space ( Σ , Ω ) . The black line represent the constraint Ω = − Σ . Figure 6: Some solutions of the full system (18) (blue) and time-averaged system (B2) (orange) for the Bianchi I metric when γ =
2. We have used for both systems the initial data sets presented in the Table I. Table I: Seven initial data sets for the simulation of the full system (18) and averaged system (B2) for the Bianchi I metric. Allinitial conditions are chosen in order to fulfill the equality Σ ( ) + Ω ( ) + Ω m ( ) =1. Sol. H ( ) Σ ( ) Ω ( ) Ω m ( ) ϕ ( ) t ( ) i 0 . . . .
09 0 0ii 0 . . . .
74 0 0iii 0 . . . .
54 0 0iv 0 .
02 0 .
48 0 .
02 0 . . .
48 0 .
02 0 . . . .
01 0 .
74 0 0vii 0 . .
685 0 .
315 0 0
Table II: Six initial data sets for the simulation of the full system (33) and time-averaged system (B3) for the flat FLRW metric( k = Ω ( ) + Ω m ( ) = Sol. H ( ) Ω ( ) Ω k ( ) Ω m ( ) ϕ ( ) t ( ) i 0 . . . . . . . . . . . . . . . . . . ∂ τ ¯ Ω = − ¯ Ω ( γ − ) (cid:0) ¯ Ω − (cid:1) , ∂ τ ¯ Φ = , ∂ τ H = − H (cid:2) γ ( − Ω ) + Ω (cid:3) , ∂ τ t = / H , (B3)Finally, for dust γ = f = b µ ω − µ ≥ γ = O ( H ) ): ∂ τ ¯ Ω = H ¯ Ω ( ω − µ ) b µ ω ∂ τ ¯ Φ = H ¯ Ω ( ω − µ ) b µ ω − H ¯ Ω ( ω − µ ) b µ ω , ∂ τ H = − H − H ¯ Ω ( ω − µ ) b µ ω , ∂ τ t = / H . (B4)Taking as initial conditions we use the six data setspresented in the Table II. In Figures 7(a)-7(d) are pre- sented projections of some solutions of the full system(33) and time-averaged system (B3) for the flat FLRWmetric ( k =
0) in the ( H , Ω ) space considering in bothsystems the same initial data sets from table II. Fig-ure 7(a) shows solutions for a background fluid corre-sponding to cosmological constant ( γ = γ = γ = / γ = H tends asymptot-ically to constant values depending on the initial con-ditions which is consistent to de Sitter expansion (seefigure 7(a)).For γ =
1, we have from the qualitativeanalysis in section 4.2 that lim τ → + ∞ ¯ Ω ( τ ) = b µ ω ¯ Ω (cid:113) b µ ω + H ¯ Ω ( µ − ω ) and lim τ → + ∞ H ( τ ) = H and ¯ Ω are the initial values when τ = H → Ω lying on the orange line (solution of the time-averagedsystem at fourth order) give an upper bound to thevalue Ω of the original system. Therefore, by control-ling the error of the averaged higher order system onecan also control the error of the original system. [1] A. Paliathanasis, L. Karpathopoulos, A. Wojnar andS. Capozziello, Eur. Phys. J. C (2016) no.4, 225 [2] S. Basilakos, M. Tsamparlis and A. Paliathanasis, H iiiiiiivvvi (a) Projections in the space ( H , Ω ) for γ = H iiiiiiivvvi (b) Projections in the space ( H , Ω ) for γ = H iiiiiiivvvi (c) Projections in the space ( H , Ω ) for γ = / H iiiiiiivvvi (d) Projections in the space ( H , Ω ) for γ = Figure 7: Some solutions of the full system (33) (blue) and time-averaged system (B3) or (B4) if γ = k =
0) when the background fluid correspond to: (a) cosmological constant, (b) dust, (c) radiation and (d) stifffluid. The black line represent the constraint Ω =
1. We have used for both systems the initial data sets presented in the TableII. Phys. Rev. D (2011), 103512[3] A. Paliathanasis, M. Tsamparlis and S. Basilakos,Phys. Rev. D (2014) no.10, 103524[4] J. D. Barrow and A. Paliathanasis, Gen. Rel. Grav. (2018) no.7, 82 doi:10.1007/s10714-018-2402-4[5] J. D. Barrow and A. Paliathanasis, Phys. Rev. D (2016) no.8, 083518[6] A. Paliathanasis, J. D. Barrow and P. G. L. Leach,Phys. Rev. D (2016) no.2, 023525[7] A. Paliathanasis, M. Tsamparlis and M. T. Mustafa,Int. J. Geom. Meth. Mod. Phys. (2015) no.03,1550033[8] M. Tsamparlis, A. Paliathanasis and L. Karpathopou-los, Gen. Rel. Grav. (2015) no.2, 15[9] F. Dumortier and R. Roussarie (1995) “Canard cy-cles and center manifolds”, (Memoirs of the AmericanMathematical Society, 577).[10] N. Fenichel (1979) Geometric singular perturbationtheory for ordinary differential equations. Journal ofDifferential Equations , 53-98[11] G. Fusco and J.K. Hale, Journal of Dynamics and Dif-ferential Equations 1, 75 (1988)[12] N. Berglund and B. Gentz, Noise-Induced Phenomena in Slow-Fast Dynamical Systems, Series: Probabilityand Applications, Springer-Verlag: London, (2006)[13] M. H. Holmes (2013) “Introduction to PerturbationMethods”, (Springer Science+Business Media NewYork, ISBN 978-1-4614-5477-9)[14] Jirair Kevorkian, J.D. Cole (1981) “Perturbation Meth-ods in Applied Mathematics” (Applied Mathemati-cal Sciences Series, Volume 34, Springer-Verlag NewYork eBook ISBN 978-1-4757-4213-8[15] Ferdinand Verhulst, (2000) “Methods and Applica-tions of Singular Perturbations: Boundary Layers andMultiple Timescale Dynamics” (Springer-Verlag NewYork, ISBN 978-0-387-22966-9)[16] G. Leon and F. O. F. Silva, Class. Quant. Grav. (2020) no.24, 245005[17] G. Leon and F. O. F. Silva, [arXiv:1912.09856 [gr-qc]][18] Genly Leon and Felipe Orlando Franz Silva 2021Class. Quantum Grav. , 667 (2007).[22] A. Alho, J. Hell and C. Uggla, Class. Quant. Grav. , no. 14, 145005 (2015)[23] J. Llibre and C. Vidal, J. Math. Phys. , 012702(2012)[24] D. Fajman, G. Heißel and M. Maliborski, Class.Quant. Grav. , no.13, 135009 (2020)[25] David Fajman, Gernot Heißel, and Jin Woo JangarXiv: 2006.12844v2[26] E. A. Coddington y Levinson, N. Theory of Ordi-nary Differential Equations, New York, MacGraw-Hill, (1955)[27] J. K. Hale, Ordinary Differential Equations, New York,Wiley (1969)[28] D. K. Arrowsmith y C. M. Place, An introductionto dynamical systems, Cambridge University Press,Cambridge, England, (1990)[29] S. Wiggins. Introduction to Applied Nonlinear dynam-ical systems and Chaos. Springer (2003)[30] L. Perko, Differential equations and dynamical sys-tems, third edition (Springer-Verlag, New York, 2001).[31] V.I. Arnold, Ordinary differential equations. Cam-bridge: M.I.T. Press., 1973[32] M. W. Hirsch and S. Smale. Differential equations, dy-namical systems, and linear algebra. New York: Aca-demic Press (1974)[33] J. Hale. Ordinary differential equations. Malabar,Florida: Robert E. Krieger Publishing Co., Inc. (1980)[34] Lasalle, J. P., J. Diff. Eq., , pp. 57-65, 1968[35] B. Aulbach, Continuous and Discrete Dynamics nearManifolds of Equilibria (Lecture Notes in Mathemat-ics No. 1058, Springer, 1984)[36] R. Tavakol, Introduction to dynamical systems, ch4. Part one, pp. 84–98, Cambridge University Press,Cambridge, England, (1997)[37] S. Foster, Class. Quant. Grav. , 3485 (1998)[38] J. Miritzis, Class. Quant. Grav. , 2981 (2003)[39] R. Giambo, F. Giannoni and G. Magli, Gen. Rel. Grav. , 21 (2009).[40] G. Leon and C. R. Fadragas, Dynamical Systems: AndTheir Applications (Saarbr¨ucken: LAP Lambert Aca-demic Publishing), arXiv:1412.5701 [gr-qc][41] G. Leon, P. Silveira and C. R. Fadragas, “Phase-spaceof flat Friedmann-Robertson-Walker models with botha scalar field coupled to matter and radiation,” in Clas-sical and Quantum Gravity: Theory, Analysis and Ap-plications ed V R Frignanni (New York: Nova SciencePublisher) ch 10 [arXiv:1009.0689 [gr-qc]][42] C. R. Fadragas and G. Leon Class. Quant. Grav. ,no. 19, 195011 (2014)[43] D. Gonz´alez Morales, Y. N´apoles Alvarez, Quintaes-encia con acoplamiento no m´ınimo a la materia oscuradesde la perspectiva de los sistemas din´amicos, Bach-elor Thesis, Universidad Central Marta Abreu de LasVillas, 2008[44] G. Leon, Class. Quant. Grav. , 035008 (2009)[45] R. Giambo and J. Miritzis, Class. Quant. Grav. ,095003 (2010)[46] K. Tzanni and J. Miritzis, Phys. Rev. D , no. 10,103540 (2014) Addendum: [Phys. Rev. D , no. 12,129902 (2014)][47] R. J. van den Hoogen, A. A. Coley and D. Wands,Class. Quant. Grav. , 1843 (1999)[48] E. J. Copeland, A. R. Liddle and D. Wands, Phys. Rev. D , 4686 (1998)[49] A.A. Coley, 2003, Dynamical systems and cosmology(Kluwer Academic, Dordrecht: ISBN 1-4020-1403-1). doi:10.1007/978-94-017-0327-7 pages 7-26[50] A. A. Coley, Introduction to Dynamical Systems. Lec-ture Notes for Math 4190/5190 (1994).[51] A. A. Coley, gr-qc/9910074.[52] Bassemah Alhulaimi (2017), Einstein-Aether Cosmo-logical Scalar Field Models (Phd Thesis, DalhousieUniversity).[53] V. G. LeBlanc, D. Kerr and J. Wainwright, Class.Quant. Grav. , 513 (1995).[54] J. M. Heinzle and C. Uggla, Class. Quant. Grav. ,015009 (2010).[55] A. A. Coley, W. C. Lim and G. Leon,[arXiv:0803.0905 [gr-qc]].[56] C. R. Fadragas, G. Leon and E. N. Saridakis, Class.Quant. Grav. (2014) 07501[57] A. S. Kompaneets and A. S. Chernov, Zh. Eksp. Teor.Fiz. (J. Exptl. Theoret. Phys. (U.S.S.R.)) 47 (1964)1939 [Sov. Phys. JETP 20, 1303 (1965)][58] R. Kantowski and R. K. Sachs, J. Math. Phys. 7, 443(1966)[59] A. B. Burd and J. D. Barrow, Nucl. Phys. B 308, 929(1988)[60] J. Yearsley and J. D. Barrow, Class. Quant. Grav. 13,2693 (1996)[61] S. Byland and D. Scialom, Phys. Rev. D , 6065-6074(1998)[62] G. Leon, E. Gonz´alez, S. Lepe, C. Michea andA. D. Millano, [arXiv:2102.05551 [gr-qc]].[63] Edward Arthur Milne, “Relativity, Gravitation andWorld Structure”, Oxford University Press, 1935.[64] S. M. Carroll, “Spacetime and Geometry,” San Fran-cisco, USA: Addison-Wesley (2004) 513 p[65] V. Mukhanov, “Physical Foundations of Cosmology,”UK: Cambridge University Press (2005), 27p[66] C. W. Misner, K. S. Thorne and J. A. Wheeler, “Grav-itation,” San Francisco 1973, 1279p[67] U. Nilsson and C. Uggla, Class. Quant. Grav. 13, 1601(1996)[68] C. Uggla and H. Zur-Muhlen, Class. Quant. Grav. 71365 (1990)[69] M. Goliath, U. S. Nilsson and C. Uggla, Class. Quant.Grav. 15, 167 (1998)[70] B. J. Carr, A. A. Coley, M. Goliath, U. S. Nilsson andC. Uggla, Class. Quant. Grav. 18, 303 (2001)[71] A. Coley and M. Goliath, Class. Quant. Grav. 17, 2557(2000)[72] A. Coley and M. Goliath, Phys. Rev. D 62, 043526(2000)[73] K.C. Jacobs, Astrophys J. 153, 661 (1968)[74] C.B Collins and S.W. Hawking, Astroph. J. 180, 317(1973)[75] J.D. Barrow, Mon. Not. R. astron. Soc. 175, 359 (1976)[76] J.D. Barrow and D.H. Sonoda, Phys. Reports, 139, 1(1986)[77] M. Thorsrud, B.D. Normann and T.S. Pereira, Class.Quantum Grav. , 065015 (2020)[78] M. Heusler, Phys. Lett. B 253, 33 (1991)[79] J.M. Aguirregabiria, A. Feinstein and J. Ibanez, Phys.Rev. D 48, 4662 (1993) [80] T. Christodoulakis, Th. Grammenos, Ch. Helias andP.G. Kevrekidis, J. Math. Phys. 47, 042505 (2006)[81] M. Tsamparlis and A. Paliathanasis, Gen. Relat.Gravit. 43, 1861 (2011)[82] A.A. Coley, J. Ibanez and R.J. van den Hoogen, J.Math. Phys. 38, 5256 (1997)[83] J. Ibanez, R.J. van den Hoogen and A.A. Coley, Phys.Rev. D 51, 928 (1995)[84] K. Adhav, A. Nimkar, R. Holey, Int. J. Theor. Phys. L1 (2006)[89] T. Clifton and J.D. Barrow, Class Quantum Grav. , 347 (1981) [Adv. Ser.Astrophys. Cosmol. , 139 (1987)][94] A. D. Linde, Phys. Lett. B (1983), 177-181[95] A. D. Linde, Phys. Lett. B (1986), 395-400[96] A. D. Linde, [arXiv:hep-th/0205259 [hep-th]][97] A. H. Guth, J. Phys. A (2007), 6811-6826[98] M. Sharma, M. Shahalam, Q. Wu and A. Wang, JCAP , 003 (2018)[99] L. McAllister, E. Silverstein, A. Westphal andT. Wrase, JHEP (2014), 123[100] J. Wainwright and G. F. R. Ellis, Eds. Dynamical Sys-tems in Cosmology. Cambridge Univ. Press, Cam-bridge, 1997[101] N. Aghanim et al. [Planck Collaboration] Astron. As-trophys.641