Averaging Generalized Scalar Field Cosmologies I: Locally Rotationally Symmetric Bianchi III and open Friedmann-Lemaître-Robertson-Walker models
Genly Leon, Esteban González, Samuel Lepe, Claudio Michea, Alfredo D. Millano
AAveraging Generalized Scalar Field Cosmologies I: Locally Rotationally SymmetricBianchi III and Open Friedmann-Lemaˆıtre-Robertson-Walker Models
Genly Leon, ∗ Esteban Gonz´alez,
2, †
Samuel Lepe,
3, ‡
Claudio Michea,
1, § and Alfredo D. Millano
1, ¶ 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 given by a barotropicEquation of State with barotropic index γ are investigated for Locally Rotationally Symmetric Bianchi III met-ric and for Open Friedmann-Lemaˆıtre-Robertson-Walker (FLRW) metric. Using methods from the Theory ofAveraging of Nonlinear Dynamical Systems and numerical simulations, it is proved that the full equations ofthe time-dependent system and their corresponding time-averaged versions have the same late-time dynamics.Therefore, the simplest time-averaged system determines the future asymptotic of the full system. In particular,depending on values of the barotropic index γ are found for LRS Bianchi III metric the generic late-time attrac-tors of physical interests given by Bianchi III flat spacetime, flat matter dominated FLRW universe (mimickingde Sitter, quintessence or zero acceleration solutions) and a matter-curvature scaling solution. For Open FLRWmetric the late-time attractors are the flat matter dominated FLRW universe (mimicking de Sitter, quintessenceor zero acceleration solutions) and the Milne solution. With this approach, the oscillations entering the fullsystem through Klein-Gordon equation can be controlled and smoothed out as the time-dependent perturbativeparameter given by the Hubble factor H tends 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 oscillator
Contents
1. Introduction
2. Spatial homogeneous and anisotropic scalarfield cosmologies
3. Averaging scalar field cosmologies k = −
1. 9
4. Qualitative analysis of averaged systems k = −
1. 154.2.1. Late-time behaviour 18
5. Conclusions Acknowledgements A. Main equations ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] § Electronic address: [email protected] ¶ Electronic address: [email protected]
B. Center Manifold Calculations F for γ = D for γ ≥ C. Proof of Theorem 7 D. Numerical simulation k = − References
1. INTRODUCTION
Scalar fields have played important roles in the phys-ical description of the Universe. For example, thescalar field known as inflaton provides the main mech-anism to explain the early accelerated era of the Uni-verse known as inflation [1]. In addition, the mostsimplest models to explain the late time acceleratedera of the Universe are a Cosmological Constant or aQuintessence scalar field [2–5]. Phantom scalar fieldcan also explain the late time accelerated expansion ofthe Universe but, as a difference with Quintessence,it has a non-canonical negative kinetic energy termwhich allows to the Equation of State (EoS) parameterof the cosmological fluid crossing below the boundaryof the Cosmological Constant (phantom divide line).However, phantom models suffer ghosts instabilities[6]. Quintom scalar field models comprise both scalarfields: a Quintessence field and a Phantom field, where a r X i v : . [ g r- q c ] F e b the two scalar fields are minimally coupled to gravityand interact between them. EoS of the quintom modelcrosses the phantom divide line without the presenceof ghosts [7–24]. On the other hand, Chiral cosmol-ogy is defined in terms of two scalar fields where thetwo kinetic energies of the scalar fields define a two-dimensional manifold of constant nonzero curvatureand the scalar fields EoS parameter is bounded in theinterval [ − , ] [25, 26]. However, quantum transitionsduring the early Universe allow the crossing of phan-tom divide line [24]. These models are used to de-scribe the early time accelerated expansion of the Uni-verse. They can be used as well as unified scenarioswhich describe both Dark Energy (DE) and Dark Mat-ter (DM) in a single fluid [27]. Multi-scalar field mod-els can provide a simple mechanism to describe vari-ous epochs of the cosmological history [28–30]. Jor-dan theory [31], Brans-Dicke theory [32], Horndeskitheories [33], extended Quintessence, Modified Grav-ity, Horava-Lifschitz, the Galileons, Phantom fields, in-teracting models, etc., have been exhaustively studied,e.g., in [34–93].Several studies provide the analysis of scalar fieldcosmologies with arbitrary potentials and arbitrarycouplings [94–109]. By means of qualitative tech-niques of dynamical systems [110–126], the stabilityanalysis of solutions can be examined. In addition, pre-cise schemes to find analytical approximations of solu-tions, as well as exact solutions or solutions in quadra-ture by choosing various approaches are provided. Forexample, Perturbation Methods and Averaging Meth-ods were used in [127, 128] with interest in early andlate-time dynamics. Slow-Fast Methods were used forexample in theories based on a Generalized Uncer-tainty Principle (GUP), say in [129]. The Amplitude-Angle Transformation was used in [130, 131] to studyoscillations of the scalar field driven by generalizedharmonic potentials. In reference [132] interactingscalar field cosmologies with generalized harmonicpotentials for flat and negatively curved Friedmann-Lemaˆıtre-Robertson-Walker (FLRW) metrics and forBianchi I metrics were studied. Asymptotic and Aver-aging Methods were used to obtain stability conditionsfor several solutions of interest as H →
0, where H is the Hubble parameter. These results suggested thatthe asymptotic behavior of the time-averaged model isindependent of the coupling function and the geome-try. In reference [133] was used the Averaging The-ory to study the periodic orbits of Hamiltonian sys-tems describing a universe filled with a scalar field. Inreference [134] was studied the future asymptotics ofLocally Rotationally Symmetric (LRS) Bianchi typeIII cosmologies with a massive scalar field by usingthe Theory of Averaging of Nonlinear Dynamical Sys-tems to control the oscillations entering the full systemthrough the Klein Gordon equation. In [135] was pre-sented a theorem about the large-time behaviour for so- lutions of a general class of systems in spatially homo-geneous cosmology with oscillatory behaviour writtenin “standard form” as the first order Taylor expansionwith respect to H centered at H =
0, when H is nonnegative and monotonic decreasing to zero.In references [131, 136] were studied scalar fieldcosmologies with arbitrary potential (and with an ar-bitrary coupling to matter). In particular, generalizedharmonic potentials and exponential couplings to mat-ter in the sense of [100–103] were examined. This pa-per is a sequel of this, where we use Asymptotic Meth-ods and Averaging Theory [137–143] to obtain relevantinformation about the solution’s space of scalar fieldcosmologies: (i) in the vacuum, (ii) in the presence ofmatter. For the scalar field potential we assume the har-monic potential plus cosine-like corrections with smallphase V ( φ ) = µ (cid:104) φ µ + b f (cid:16) − cos (cid:16) φ f (cid:17)(cid:17)(cid:105) , b > , in-spired in the potential of [144, 145]. In section 2.2 wemake some comments about this potential related to itsphysical and mathematical properties. We provide fur-ther elements to justify this choice. As in [134], weconstruct an averaged version of the original systemwhere the oscillations of the solutions are smoothedout. Then, the analysis of the system is reduced tostudy the corresponding averaged system.This research program – named “Averaging Gener-alized Scalar Field Cosmologies”– has three steps ac-cording to the three cases of study: (I) Bianchi IIIand Open FLRW model, (II) Bianchi I and Flat FLRWmodel and (III) Kantowski-Sachs and Closed FLRW.This paper is devoted to the case I and cases II and IIIwill be studied in two companion papers [146, 147].The paper is organized as follows. In section 2 we in-troduce the model. In section 3 we apply AveragingMethods to analyze the periodic solutions of a scalarfield with self-interacting potentials within the class ofgeneralized harmonic potentials [130]. In particular,in section 3.1 LRS Bianchi III metric is studied. Insection 3.2 is investigated FLRW metrics with nega-tive curvature. In section 4 we study the resulting av-eraged systems. Bianchi III metric is studied in section4.1 and section 4.2 is devoted to FLRW metric withnegative curvature (Open FLRW). Finally, in section5 are discussed our main results. Our main equationsare presented in appendix A. In appendix B are pre-sented center manifold calculations for nonhyperbolicequilibrium points. In appendix C is proved our maintheorem. In appendix D numerical evidence supportingthe results of section 3 is presented.
2. SPATIAL HOMOGENEOUS AND ANISOTROPICSCALAR FIELD COSMOLOGIES
The spatial homogeneous but anisotropic spacetimesare known as either Bianchi or Kantowski-Sachs cos-mologies. In Bianchi models, the spacetime manifoldis foliated along the time axis with three dimensionalhomogeneous hypersurfaces. Bianchi spacetimes con-tain many important cosmological models that havebeen used to study anisotropies of primordial Universeand its evolution towards the observed isotropy of thepresent epoch [148–151]. The list includes the stan-dard FLRW model in the limit of the isotropizationwith Bianchi III isotropizing to Open FLRW modelsand Bianchi I isotropizing to Flat FLRW models. InGeneral Relativity (GR) the Hubble parameter H isalways monotonic for Bianchi I and Bianchi III. ForBianchi I the anisotropy decays on time for H > The action integral of interest is (cid:90) d x (cid:112) | g | (cid:20) R − g µν ∇ µ φ ∇ ν φ − V ( φ ) + L m (cid:21) . (1)It is expressed in a system of units in which 8 π G = c = (cid:125) =
1, where R is the curvature scalar, φ is the scalar field, ∇ α is the covariant derivative and the potential is V ( φ ) = µ (cid:20) b f (cid:18) − cos (cid:18) φ f (cid:19)(cid:19) + φ µ (cid:21) , b > . (2)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:104) d ϑ + k − sin ( √ k ϑ ) d ζ (cid:105) , (3)where e , e and e = √ ke / sin ( √ k ϑ ) are func-tions of t which are the components of the frame vec-tors [177]: e = ∂ t , e = e ∂ r , e = e ∂ ϑ , e = e ∂ ζ .According to the values of k we will divide the studyin three cases: (I) LRS Bianchi III and Open FLRWmodel ( k = − k = k =+ e ( t ) (cid:54) = e ( t ) andFLRW metrics satisfy e ( t ) = e ( t ) = a ( t ) − (where a ( t ) is the scale factor) [61]. The average scale factor a ( t ) and the anisotropic parameter σ + ( t ) are defined by˙ aa = H : = − ddt ln (cid:2) e ( t )( e ( t )) (cid:3) , (4) σ + = ddt ln (cid:2) e ( t )( e ( t )) − (cid:3) . (5)Taking variation of (1) for the 1-parameter family ofmetrics (3) leads to [61]:3 H + kK = σ + + ρ m +
12 ˙ φ + V ( φ ) , (6) − ( σ + + H ) − σ + − H − kK = ( γ − ) ρ m +
12 ˙ φ − V ( φ ) , (7) − σ + + σ + H − H + ˙ σ + − H = ( γ − ) ρ m +
12 ˙ φ − V ( φ ) , (8)where is used for the matter component a barotropicEoS p m = ( γ − ) ρ m with p m the pressure of the mat-ter component, ρ m the energy density and γ ∈ [ , ] isknown as barotropic index.The Gauss curvature and 3-curvature scalar are K = ( e ( t )) , ( ) R = kK . (9)Furthermore, the evolution of Gauss curvature is˙ K = − ( σ + + H ) K , (10)while the evolution for e is given by [177]:˙ e = − ( H − σ + ) e . (11)From Eqs. (7), (8) is obtained the shear equation:˙ σ + = − H σ + − kK . (12)Eqs. (6), (7), (8), (12) give the Raychaudhuri equation:˙ H = − H − σ + − ( γ − ) ρ m −
13 ˙ φ + V ( φ ) . (13)Finally, the matter and Klein-Gordon equations are:˙ ρ m = − γ H ρ m , (14)¨ φ = − H ˙ φ − dV ( φ ) d φ . (15) Chaotic inflation is a model of cosmic inflationwhich takes the potential term V of the hypotheticalinflaton field φ to be V ( φ ) = m φ φ , the so-called har-monic potential ( φ -interaction) [178–181]. Whereasother inflationary models assume a monotonously de-creasing potential with φ ; assuming in an ad hoc waythat the inflaton field has a large amplitude “at the BigBang” to then slowly “roll down” the potential. Theidea of [178] is that instead of the inflaton rolls downand sits on its potential minimum at equilibrium, quan-tum fluctuations stochastically (“chaotically”) drive itout of its minimum back and forward. Wherever thishappens cosmic inflation sets in and blows up the re-gion of the ambient spacetime in which the inflatonhappened to fluctuate out of its equilibrium.Relevant experimental results disfavoring the φ -interaction are due to [182, 183]. These results statethat chaotic inflation generically predicts large valuesof the tensor-to-scalar ratio r . In contrast to recentmeasurements which show low upper bounds on r .Notwithstanding, we investigate variations of the φ -potential and we do not refer to tensor-to-scalar ratioissue for the potential (19), which is out of the scope ofthe present research.The scalar field potential V ( φ ) of interest in this re-search is given by (2). It is related but not equal to themonodromy potential of [144, 145] in the context ofloop-quantum gravity, which is a particular case of thegeneral monodromy potential [145]: V ( φ ) = µ − p φ p + ∑ i Λ i cos (cid:18) φ f i + δ i (cid:19) . (16)In references [130–132] were proved that the potentialof [144, 145] for p = V ( φ ) = µ − p φ p + ∑ i Λ i (cid:20) − cos (cid:18) φ f i + δ i (cid:19)(cid:21) , (17)given by V ( φ ) = φ + f (cid:20) − cos (cid:18) φ f (cid:19)(cid:21) . (18)The family of potentials (17) provides non-negative lo-cal minimums which can be related to a late-time ac-celerated Universe. The potential (18) is obtained bysetting µ = √ and b µ = χ = χ e λφ − γ where λ is a constantand the barotropic index satisfies 0 ≤ γ ≤ , γ (cid:54) = for the FLRW metrics with k = − , φ = φ ∗ whenever φ ∗ is a localnon zero minimum of V ( φ ) . For FLRW metrics, theglobal minimum is unstable to curvature perturbationsfor γ > . Therefore, the result in [107] is confirmed,that for γ > / V ( ) = V ( φ ) = Λ (cid:104) − cos (cid:16) φ f + δ (cid:17)(cid:105) which have been used in the con-text of axion models [184]. Axionic dark matter modelwith a modified periodic potential for the pseudoscalarfield V ( φ , Φ ∗ ) = m A Φ ∗ π (cid:104) − cos (cid:16) πφ Φ ∗ (cid:17)(cid:105) in the frame-work of the axionic extension of the Einstein-aethertheory has been studied in [185]. This periodic poten-tial has minima at φ = n Φ ∗ , n ∈ Z , whereas maximaare found when n → m + . Near the minimum, when φ = n Φ ∗ + ψ and | ψ | is small, V → m A ψ where m A theaxion rests mass.The previous statements justify the study of the po-tential (2), which can be expressed as V ( φ ) = µ φ + f (cid:0) ω − µ (cid:1) (cid:18) − cos (cid:18) φ f (cid:19)(cid:19) , (19)by imposing the condition b µ + f µ − f ω = ω ∈ R and satis-fies ω − µ > b > V is a real-valued smooth function V ∈ C ∞ ( R ) with lim φ →± ∞ V ( φ ) = + ∞ . - - f V H f L f f f + H - cos H f LL (a) V ( φ ) defined by (19) for µ = √ , ω = (cid:114)(cid:12)(cid:12)(cid:12) f − f (cid:12)(cid:12)(cid:12) i , f = . - - f V H f L f f f + H - cos H f LL (b) V ( φ ) defined by (19) for µ = √ , ω = √ , f = Figure 1: Generalized harmonic potentials. Comparison with φ -squared potentials. V is an even function V ( φ ) = V ( − φ ) .3. V ( φ ) has always a local minimum at φ = V ( ) = , V (cid:48) ( ) = , V (cid:48)(cid:48) ( ) = ω > φ c (cid:54) = µ φ c + f (cid:0) ω − µ (cid:1) sin (cid:16) φ c f (cid:17) =
0, whichare extreme points of V ( φ ) . They are local maxi-mums or local minimums depending on whether V (cid:48)(cid:48) ( φ c ) : = µ + (cid:0) ω − µ (cid:1) cos (cid:16) φ c f (cid:17) < V (cid:48)(cid:48) ( φ c ) >
0. For | φ c | > f ( ω − µ ) µ = φ ∗ , this setis empty.5. There exist V max = max φ ∈ [ − φ ∗ , φ ∗ ] V ( φ ) and V min = min φ ∈ [ − φ ∗ , φ ∗ ] V ( φ ) =
0. The function V has no upper bound but it has a lower boundequal to zero.Finally, the asymptotic features of the potential (19) arethe following. Near the global minimum φ =
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- correctionis bounded, then, V ( φ ) ∼ µ φ + O ( ) as φ → ± ∞ . This makes it suitable to describe oscillatory behaviorin cosmology.Setting µ = √ , ω = √
2, we have V ( φ ) = φ + f (cid:20) − cos (cid:18) φ f (cid:19)(cid:21) . (20)Setting µ = √ , ω = (cid:113) f − f , the potential of [130–132]is recovered: V ( φ ) = φ + f (cid:20) − cos (cid:18) φ f (cid:19)(cid:21) . (21)Although (21) can be derived as a particular case of ourstudy, the cases of interest f (cid:28) ω = (cid:114)(cid:12)(cid:12)(cid:12) f − f (cid:12)(cid:12)(cid:12) i , in contradiction to ω ∈ R .Therefore, the potential (19) will not contains the po-tential studied in [130–132], unless we set f > ω = (cid:113) f − f . The potentials (20) and (21) are presentedin Figure 1.
3. AVERAGING SCALAR FIELD COSMOLOGIES
Given the differential equation ˙ x = f ( t , x , ε ) , with f periodic in t . One approximation scheme which can beused to solve the full problem is first solving the unper-turbed problem ˙ x = f ( t , x , ) by setting ε = ε . Next, we comment the case whenthe perturbations are measured by a time-dependingperturbation parameter ε ( t ) . In the approximation ofa vector function x ε ( t ) for t ≥ y ε ( t ) ; we say that y ε ( t ) is an O ( ε m ) -approximation of x ε ( t ) in the time scale 1 / ε n if for t ≥ x ε ( t ) − y ε ( t ) = O ( ε m ) , ≤ ε n t ≤ C follows, where theconstants m , n and C are independent of ε .For example, take this simple equation¨ φ + φ = ε ( − φ ) (22)with φ ( ) and ˙ φ ( ) given. The unperturbed prob-lem ¨ φ + φ = φ ( t ) = r cos ( t − ϕ ) , φ ( t ) = r sin ( t − ϕ ) , where r and ϕ are con-stants depending on the initial conditions. Let bedefined the Amplitude-Phase Transformation ([143],chapter 11): ˙ φ ( t ) = r ( t ) cos ( t − ϕ ( t )) , φ ( t ) = r ( t ) sin ( t − ϕ ( t )) , (23) such that r = (cid:113) ˙ φ ( t ) + φ ( t ) , ϕ = t − tan − (cid:18) φ ( t ) ˙ φ ( t ) (cid:19) . (24) Then, under the coordinate transformation ( ˙ φ , φ ) → ( r , ϕ ) given by (24), equation (22) becomes,˙ r = − r ε cos ( t − ϕ ) , ˙ ϕ = − ϕε sin ( ( t − ϕ )) . (25)These equations mean that r and ϕ are slowly varyingwith time, and system takes the form ˙ y = ε f ( y ) . Theidea is consider only the nonzero average of the right-hand-sides keeping r and ϕ fixed and leaving out theterms with average zero and ignoring the slow-varyingdependence of r and ϕ on t in the averaging process.Now, we average by replacing r and ϕ by their aver-aged approximations ¯ r and ¯ ϕ to obtain the system˙¯ r = − ε π (cid:90) π r cos ( t − ϕ ) dt = − ε ¯ r , (26)˙¯ ϕ = − ε π (cid:90) π sin ( ( t − ϕ )) dt = . (27) Solving this with ¯ r ( ) = r and ¯ ϕ ( ) = ϕ , the approx-imation takes the form ¯ φ = r e − ε t sin ( t − ϕ ) whichcoincides with the result that would be obtained usingthe two-timing expansion procedure. These two pro-cedures alleviate the failure of the regular asymptoticexpansion that would yield spurious secular terms inthe asymptotic expansion. Say, on the regular asymp-totic expansion x ( t , ε ) = sin t − ε t sin t + O ( ε ) , the“next to the leading term” ε t sin t dominates over scales ε t = O ( ) .It is important to know the nature of the averagedtransformations. Therefore, we enumerate the follow-ing theorems. Theorem 1 (Theorem 11.1 of [143])
Let be the n- di-mensional initial value problem in standard form: ˙ x = ε f ( t , x ) + ε g ( t , x , ε ) , x ( ) = x , t ≥ . (28) Supposing that f ( t , x ) is T -periodic in t, with T > a constant independent of ε . Performing the averag-ing process ¯ f ( y ) = T (cid:82) T f ( t , y ) dt, where y is consid-ered as a parameter that is kept constant during in-tegration. Let be the associated initial value problem ˙ y = ε ¯ f ( y ) , y ( ) = x . Then, we have y ( t ) = x ( t )+ O ( ε ) on the time scale / ε under fairly general conditions:1. The vector functions f and g are continuously dif-ferentiable in a bounded n-dimensional domainD with x as interior point on the time scale / ε .2. y ( t ) remains interior to the domain D on the timescale / ε to avoid boundary effects. Theorem 2 (Theorem 11.2 of [143])
Let be the n- di-mensional initial value problem in standard form: ˙ x = ε N ∑ n = f n ( t , x ) , x ( ) = x , (29) where the vector functions f n ( t , x ) have N independentperiods T n , n = , . . . , N. Averaging over the respectiveperiods, we obtain the associated system: ˙ y = ε N ∑ n = T n (cid:90) T n f n ( t , y ) dt , y ( ) = x . (30) As before, y ( t ) = x ( t ) + O ( ε ) on the time scale / ε . There are two important restrictions: N must be finiteand must be possible to write down the right-hand-sideof the equation as a sum of vector functions where eachof them is periodic with one period.Now, supposing that the slowly varying system ˙ x = ε f ( t , x ) is such that f ( t , x ) is not periodic, nor a finitesum of periodic vector fields as before we have the fol-lowing result. Theorem 3 (Theorem 11.3 of [143])
Let be the n- di-mensional initial value problem in standard form: ˙ x = ε f ( t , x ) + ε g ( t , x , ε ) , x ( ) = x , t ≥ . (31) Supposing that f ( t , x ) can be averaged over t in thesense that the limit ¯ f ( y ) = lim T → ∞ T (cid:90) T f ( t , y ) dt (32) exists where y is again considered as a parameter thatis kept constant during integration. Let be the associ-ated initial value problem ˙ y = ε ¯ f ( y ) , y ( ) = x . (33) Then, we have y ( t ) = x ( t ) + O ( δ ( ε )) (34) on the timescale / ε under fairly general conditions:1. The vector functions f and g are continuously dif-ferentiable in a bounded n-dimensional domainD with x an interior point on the timescale / ε .2. y ( t ) remains interior to the domain D on thetimescale / ε to avoid boundary effects.For the error δ ( ε ) , we have the explicit estimate δ ( ε ) = sup x ∈ D sup ≤ ε t ≤ C ε (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:90) t ( f ( s , x ) − ¯ f ( x )) ds (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (35) where C is a constant independent of ε . For periodic solutions we have the following.
Theorem 4 (Theorem 11.4 of [143])
Consider the n-dimensional system ˙ x = ε f ( t , x ) with f ( t , x ) T -periodic.Assume the time-averaged system ˙ y = ε ¯ f ( y ) , where ¯ f ( y ) = T (cid:90) T f ( t , y ) dt (36) admits the stationary solution y ( ¯ f ( y ) = ). If1. f ( t , x ) is a smooth vector field and2. for the Jacobian in y we have (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∂ ¯ f ∂ y (cid:12)(cid:12)(cid:12) y = y (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:54) = . (37) Then, there exists a T -periodic solution of the equation ˙ x = ε f ( t , x ) in an ε - neighborhood of x = y . We can establish the stability of a periodic solution asit matches exactly the stability of a stationary solutionof the time-averaged system. This reduces the stability problem of a periodic solution to determine the eigen-values of a matrix.In [135] it was presented a general theorem about thelarge-time behavior for solutions of a general class ofsystems in standard form: (cid:18) ˙ H ˙ x (cid:19) = H (cid:18) f ( t , x ) (cid:19) + H (cid:18) f [ ] ( t , x ) (cid:19) , (38)where H is the Hubble parameter and it is positivestrictly decreasing in t and lim t → ∞ H ( t ) =
0, where 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 = ω in an amplitude-phase transformation.The following theorem gives local-in-time asymp-totics for system (38). Let the norm (cid:107) · (cid:107) denotes thestandard discrete (cid:96) - norm (cid:107) u (cid:107) : = ∑ ni | u i | for u ∈ R n .Let also L ∞ t , x denotes the standard L ∞ space in both t and x variables with the corresponding norm definedas (cid:107) f (cid:107) L ∞ t , x : = sup t , x | f ( t , x ) | . Theorem 5 (Theorem 2.1 of [135])
Suppose H ( t ) > is strictly decreasing in t and lim t → ∞ H ( t ) = . Fix any ε = with ε < H ( ) and define t ∗ > such that ε = H ( t ∗ ) . Suppose that (cid:107) f (cid:107) L ∞ t , x , (cid:107) f [ ] (cid:107) L ∞ t , x < ∞ , f ( t , x ) isLipschitz continuous and f [ ] is continuous with respectto x for all t ≥ t ∗ . Also, assume that f and f [ ] areT -periodic for some T > . Then for all t > t ∗ witht = t ∗ + O ( H ( t ∗ ) − δ ) for some δ ∈ ( , ) we have x ( t ) − z ( t ) = O ( H ( t ∗ ) min { , − δ } ) , where x is the solution ofsystem (38) with initial condition x ( ) = x and z ( t ) isthe solution of the time-averaged system ˙ z = H ( t ∗ ) ¯ f ( z ) , for t > t ∗ with initial condition z ( t ∗ ) = x ( t ∗ ) where the time-averaged vector ¯ f is defined as ¯ f ( z ) = T (cid:90) t ∗ + Tt ∗ f ( z , s ) ds . The following theorem gives global-in-time asymp-totics for system (38).
Theorem 6 (Theorem 2.2 of [135])
Under the hy-potheses of Theorem 5, lim τ → ∞ (cid:107) x ( τ ) − z ( τ ) (cid:107) = . Therefore, H can be used as a time-dependent pertur-bation parameter. This result has applications for spa-tially homogeneous cosmologies with oscillatory be-haviour as H → Firstly, we consider LRS Bianchi III metrics for thegeneralized harmonic potential (2) minimally coupledto matter with field equations:¨ φ = − H ˙ φ − V (cid:48) ( φ ) , (39)˙ ρ m = − γ H ρ m , (40)˙ K = − ( σ + + H ) K , (41)˙ H = − H − σ + − ( γ − ) ρ m −
13 ˙ φ + V ( φ ) , (42)˙ σ + = − H σ + + K , (43)3 H = σ + + ρ m +
12 ˙ φ + V ( φ ) + K . (44) Defining Ω = (cid:115) ω φ + ˙ φ H , Σ = σ + H , Ω k = K H , Φ = t ω − tan − (cid:18) ωφ ˙ φ (cid:19) , (45)we obtain the full system (A2).Assuming ω > µ and b µ + f µ − f ω = f = b µ ω − µ (which is equivalent to tunethe angular frequency ω ); are eliminated the undesiredterms evolving as ∝ H in the series expansion around H = x = H f ( t , x ) + O ( H ) , x = ( Ω , Σ , Ω k , Φ ) T , (46)˙ H = − H (cid:18) γ ( − Σ − Ω k − Ω ) + Σ + Ω k + Ω cos ( t ω − Φ ) (cid:19) + O ( H ) , (47)where f ( t , x ) = Ω (cid:16) − ( γ − ) Σ + ( − γ ) Ω k + (cid:0) Ω − (cid:1) ( − γ + ( t ω − Φ )) (cid:17) (cid:16) Ω k (( − γ ) Σ + ) + Σ (cid:32) − ( γ − ) Σ + γ + Ω ( − γ + ( t ω − Φ )) − (cid:17)(cid:33) Ω k (cid:16) − γ (cid:0) Σ + Ω + Ω k − (cid:1) + Σ − Σ + Ω cos ( t ω − Φ ) + Ω k − (cid:17) − sin ( t ω − Φ ) . (48)Replacing ˙ x = H f ( t , x ) , x = ( Ω , Σ , Ω k , Φ ) T where f ( t , x ) is defined by (48) with ˙ y = H ¯ f ( y ) , y = (cid:0) ¯ Ω , ¯ Σ , ¯ Ω k , ¯ Φ (cid:1) T and¯ f ( y ) given by time-averaging ¯ f ( · ) : = L (cid:90) L f ( · , t ) dt , L = πω (49)we obtain: ˙¯ Ω = H ¯ Ω (cid:0) − γ (cid:0) ¯ Σ + ¯ Ω + ¯ Ω k − (cid:1) + Σ + Ω + Ω k − (cid:1) , (50)˙¯ Σ = H (cid:0) ¯ Σ (cid:0) − γ (cid:0) ¯ Σ + ¯ Ω + ¯ Ω k − (cid:1) + Σ + Ω + Ω k − (cid:1) + Ω k (cid:1) , (51)˙¯ Ω k = − H ¯ Ω k (cid:0) γ (cid:0) ¯ Σ + ¯ Ω + ¯ Ω k − (cid:1) − Σ + Σ − Ω − Ω k + (cid:1) , (52)˙¯ Φ = , (53)˙ H = − H (cid:16) γ (cid:0) − ¯ Σ − ¯ Ω − ¯ Ω k (cid:1) + Σ + Ω + Ω k (cid:17) . (54) Theorem 7
Let be defined the functions ¯ Ω , ¯ Σ , ¯ Ω k , ¯ Φ and H satisfying the time-averaged system (50) , (51) , (52) , (53) , (54) . Then, there exist functions g , g , g and g , such that if Ω = Ω + Hg ( H , Ω , Σ , Ω k , Φ , t ) , Σ = Σ + Hg ( H , Ω , Σ , Ω k , Φ , t ) , Ω k = Ω k + Hg ( H , Ω , Σ , Ω k , Φ , t ) , Φ = Φ + Hg ( H , Ω , Σ , Ω k , Φ , t ) , (55) the functions Ω , Σ , Ω k , Φ and ¯ Ω , ¯ Σ , ¯ Ω k , ¯ Φ have thesame limit as t → ∞ . Setting Σ = Σ = are derived theanalogous results for negatively curved 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 < Ω + Σ + Ω k <
1. Then, it can bedefined recursively the bootstrapping sequences t = t ∗ H = H ( t ∗ ) , t n + = t n + H n H n + = H ( t n + ) , (56)such that lim n → ∞ H n = n → ∞ t n = ∞ .2. Finding g i such that it can be defined Ω , Σ , Ω k , Φ through equation (55) andproving that g i , i = , . . . , t ∈ [ t n , t n + ] . By continuityof ¯ Ω , Ω , ¯ Σ , Σ , Φ , ¯ Ω k , Ω k , it can be definedconstants L , M , M , M , M which are finite andindependent of H .3. Using the expansions (55), defining ∆Ω = Ω − ¯ Ω , ∆Σ = Σ − ¯ Σ , ∆Ω k = Ω k − ¯ Ω k , ∆Φ = Φ − ¯ Φ and taking the initial conditions Ω ( t n ) = ¯ Ω ( t n ) = Ω n , Σ ( t n ) = ¯ Σ ( t n ) = Σ n , Ω k ( t n ) = ¯ Ω k ( t n ) = Ω kn , Φ ( t n ) = ¯ Φ ( t n ) = Φ n , with 0 < Ω ( t n ) + Σ ( t n ) + Ω k ( t n ) < 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 ) , (cid:12)(cid:12)(cid:12) ∆Ω k ( 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 .
5. Finally, taking the limit as n → ∞ , we obtain H n →
0. Then, as H n →
0, it follows the func-tions Ω , Σ , Ω k , Φ and ¯ Ω , ¯ Σ , ¯ Ω k , ¯ Φ have thesame limit as t → ∞ .The complete proof is given in appendix C. Theorem7 implies that Ω , Σ , Ω k and Φ evolve according to thetime-averaged system (50), (51), (52), (53) as H → k = − . Secondly, we consider FLRW metric with k = − φ = − H ˙ φ − V (cid:48) ( φ ) , (57a)˙ ρ m = − γ H ρ m , (57b)˙ a = aH , (57c)˙ H = − (cid:0) γρ m + ˙ φ (cid:1) + ka , (57d)3 H = ρ m +
12 ˙ φ + V ( φ ) − ka , (57e)Defining Ω = (cid:115) ω φ + ˙ φ H , Ω k = − ka H , Φ = t ω − tan − (cid:18) ωφ ˙ φ (cid:19) , (58)we obtain the equations (A5).Setting f = b µ ω − µ >
0, we obtain:˙ x = H f ( t , x ) + O ( H ) , x = ( Ω , Ω k , Φ ) T , (59)˙ H = − (cid:0) γ (cid:0) − Ω − Ω k (cid:1) + Ω k (cid:1) − Ω cos ( t ω − Φ ) + O ( H ) , (60)0where f ( t , x ) = Ω (cid:0) γ − γ (cid:0) Ω + Ω k (cid:1) + Ω k (cid:1) + Ω (cid:0) Ω − (cid:1) cos ( t ω − Φ ) − Ω k (cid:0) γ Ω + ( γ − )( Ω k − ) (cid:1) + Ω Ω k cos ( t ω − Φ ) − sin ( t ω − Φ ) . (61)Replacing ˙ x = H f ( t , x ) with f ( t , x ) defined by (61) with˙ y = H ¯ f ( y ) , y = (cid:0) ¯ Ω , ¯ Ω k , ¯ Φ (cid:1) T and using the time aver-aging (49) we obtain the time-averaged system:˙¯ Ω = − H ¯ Ω (cid:0) ( γ − ) (cid:0) ¯ Ω − (cid:1) + ( γ − ) ¯ Ω k (cid:1) , (62)˙¯ Ω k = − H ¯ Ω k (cid:0) ( γ − ) ¯ Ω − γ + ( γ − ) ¯ Ω k + (cid:1) , (63)˙¯ Φ = . (64)Theorem 7 applies to Bianchi III and the invariant set Σ =
4. QUALITATIVE ANALYSIS OF AVERAGEDSYSTEMS
According to Theorem 7, for Bianchi III and neg-atively curved FLRW metrics H plays the role of atime-dependent perturbation parameter controlling themagnitude of the error between the solutions and theoscillations are viewed as perturbations. In the time-averaged system the Raychaudhuri equation decouples.The analysis of the system is therefore reduced to thestudy the corresponding time-averaged system. In this case the averaged system is (50), (51), (52)and (53). Introducing the new time variable τ through d fd τ = H d fdt , we obtain the guiding system: d ¯ Ω d τ =
12 ¯ Ω (cid:0) − γ (cid:0) ¯ Σ + ¯ Ω + ¯ Ω k − (cid:1) + Σ + Ω + Ω k − (cid:1) , (65) d ¯ Σ d τ = (cid:0) ¯ Σ (cid:0) − γ (cid:0) ¯ Σ + ¯ Ω + ¯ Ω k − (cid:1) + Σ + Ω + Ω k − (cid:1) + Ω k (cid:1) , (66) d ¯ Ω k d τ = − ¯ Ω k (cid:0) γ (cid:0) ¯ Σ + ¯ Ω + ¯ Ω k − (cid:1) − Σ + Σ − Ω − Ω k + (cid:1) . (67) We have ¯ Σ + ¯ Ω + ¯ Ω k + ¯ Ω m =
1. Therefore, from thecondition ¯ Ω m ≥ (cid:8) ( ¯ Ω , ¯ Σ , ¯ Ω k ) ∈ R : ¯ Σ + ¯ Ω + ¯ Ω k ≤ , ¯ Ω k ≥ (cid:9) . (68)The equilibrium points of the guiding system(65), (66), (67) are:1. T : ( ¯ Ω , ¯ Σ , ¯ Ω k ) = ( , − , ) with eigenvalues (cid:8) , , − γ (cid:9) .i) It is a source for 0 ≤ γ < . ii) It is nonhyperbolic for γ = . Defining a ( t ) a = (cid:2) e ( t )( e ( t )) (cid:3) − , τ = ln (cid:18) a ( t ) a (cid:19) , (69)such that H = ˙ aa . (70)We denote by convention t = (cid:16) a ( ) a (cid:17) = e ( )( e ( )) = τ ( ) = T we obtain: ˙ H = − H ˙ a = aH = ⇒ H ( t ) = H H t + a ( t ) = a √ H t + . (71) Σ = − σ + = − H = − H H t + . This im-plies that the Gauss curvature is constant. In-deed, from ˙ K = − ( σ + + H ) K , (72)it follows K = ( e ( t )) = c − . Substituting backin equation (11) is obtained:˙ e = − H H t + e , e ( ) = c . (73)Hence, e ( t ) = c H t + . (74)1That is, the line element (3) becomes ds = − dt + ( H t + ) c dr + c (cid:2) d ϑ + ϑ d ζ (cid:3) . (75)Therefore, the corresponding solution can be ex-pressed in the form of the Taub-Kasner solution( p = , p = , p = Q : ( ¯ Ω , ¯ Σ , ¯ Ω k ) = ( , , ) with eigenvalues (cid:8) , , − γ (cid:9) . i) It is a source for 0 ≤ γ < . ii) It is nonhyperbolic for γ = . Evaluating Raychaudhuri equation (A2e) at D and using the previous argument we obtain H ( t ) = H H t + . (76) Σ = σ + = H = H H t + . Hence, Gaussequation (10) and evolution equation (11) be-comes ˙ K = − H K H t + , K ( ) = c − (77)and ˙ e = H H t + e , e ( ) = c . (78)Then, by integration, e ( t ) = c (cid:112) H t + , (79) K ( t ) = c ( H t + ) / . (80)That is, the line element (3) becomes ds = − dt + c − ( H t + ) − dr + c − ( H t + ) / (cid:2) d ϑ + ϑ d ζ (cid:3) . (81)Therefore, the corresponding solution can be ex-pressed in the form of the non-flat LRS Kasner( p = − , p = , p = ) Bianchi I solution([186] Sect. 6.2.2 and Sect. 9.1.1 (2)).3. D : ( ¯ Ω , ¯ Σ , ¯ Ω k ) = ( , , ) with eigenvalues (cid:8) − , , − γ (cid:9) . It is a nonhyperbolic saddle for0 ≤ γ < . i) For γ = {− , , } ii) For γ > D we obtain ˙ H = − H ˙ a = aH = ⇒ H ( t ) = H H t + a ( t ) = a (cid:16) H t + (cid:17) / . (82) Σ = / σ + = H = H H t + . Hence,Gauss equation (10) and evolution equation (11)become˙ e = , e ( ) = c , (83)˙ K = − H K H t + , K ( ) = c . (84)Hence, e = c , K = c ( H t + ) . (85)That is, the line element (3) becomes ds = − dt + c − dr + ( H t + ) c (cid:2) d ϑ + sinh ( ϑ ) d ζ (cid:3) . (86)Therefore, the corresponding solution can be ex-pressed in the form of the Bianchi III form of flatspacetime ([186] p 193, Eq. (9.7)).4. F : ( ¯ Ω , ¯ Σ , ¯ Ω k ) = ( , , ) with eigenvalues (cid:8) − , , − γ (cid:9) . This point is always a saddlebecause it has a negative and a positive eigen-value. For γ = ψ ( t ) = t ω − Φ ( t ) and evaluating Ray-chaudhuri equation (A2e) at D we obtain˙ H = b γ µ (cid:18) cos (cid:18) √ H ( µ − ω ) sin ( ψ ) b µ ω (cid:19) − (cid:19) ( µ − ω )+ H (cid:0) γ µ sin ( ψ ) + ω (cid:0) ( γ − ) cos ( ψ ) − γ (cid:1)(cid:1) ω . (87)Therefore,˙ H ∼ − H cos ( t ω − Φ ) , (88)for large t . In average, Φ is constant, setting forsimplicity Φ = H ( t ) = H ω H t ω + H sin ( t ω ) + ω , (89)where H is the current value of H ( t ) . Finally, H ( t ) ∼ t for large t .Gauss equation (10) and evolution equation (11)2become ˙ e = − e t , ˙ K = − K t , (90)with general solution˙ e ( t ) = c t / , K ( t ) = c t / . (91)That is, the line element (3) becomes ds = − dt + c − t / dr + c − t / (cid:2) d ϑ + ϑ d ζ (cid:3) . (92)Hence for large t and c = c the equilibriumpoint can be associated with the flat FriedmanEinstein-de-Sitter solution ([186], Sec 9.1.1 (1))with γ = F : ( ¯ Ω , ¯ Σ , ¯ Ω k ) = ( , , ) with eigenvalues (cid:110) ( γ − ) , ( γ − ) , γ − (cid:111) . i) It is a sink for 0 ≤ γ < . ii) It is a saddle for < γ < < γ < . iii) It is nonhyperbolic for γ = or γ = γ = . Evaluating Raychaudhuri equation (A2e) at F we obtain ˙ H = − γ H ˙ a = aH = ⇒ H ( t ) = H γ H t + a ( t ) = a (cid:16) γ H t + (cid:17) γ . (93)That is, 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) . (94)The corresponding solution is a flat matter dom-inated FLRW universe with ¯ Ω m = MC : ( ¯ Ω , ¯ Σ , ¯ Ω k ) = ( , γ − , − γ + γ − ) with eigenvalues (cid:110) ( γ − ) , (cid:16) γ + √ − γ (cid:112) γ ( γ − ) + − (cid:17) , (cid:16) γ − √ − γ (cid:112) γ ( γ − ) + − (cid:17)(cid:111) . Bydefinition Ω k ≥
0, therefore, we impose therestriction ≤ γ ≤ < γ < . ii) It is a saddle for 1 < γ < . iii) It is nonhyperbolic for γ = or γ = γ = . The corresponding solution is a matter-curvaturescaling solution with ¯ Ω m = ( − γ ) . We obtainthe same expressions in (93), for a , H and for ds we obtain the expression: ds = − dt + a (cid:18) γ H t + (cid:19) γ dr a (cid:18) γ H t + (cid:19) γ (cid:2) d ϑ + sinh ( ϑ ) d ζ (cid:3) . (95)In figure 2(a) are presented some orbits in the phasespace of the guiding system (65), (66), (67) for γ = F where the scalar field mimics a Cosmolog-ical Constant. The equilibrium point D is a saddle.In figure 2(b) are presented some orbits of the phasespace of the guiding system (65), (66), (67) for γ = .The point F coincides with MC ; it is asymptoticallystable as proved in appendix 1 by means of CenterManifold theory. The equilibrium point D is a saddle.In figures 3(a) and 3(b) are presented some orbits inthe phase space of the guiding system (65), (66), (67)for γ = . γ = . MC is a stable node and D is a saddle.It is worth to notice that for γ = ( ¯ Ω , ¯ Σ , ¯ Ω k ) = ( ¯ Ω ∗ , , ) and D ∗ : = ( ¯ Ω , ¯ Σ , ¯ Ω k ) = ( ¯ Ω ∗ , , ) , where ¯ Ω ∗ is an ar-bitrary number which satisfies ¯ Ω ∗ ∈ [ , ] . Therefore,the Bianchi III flat spacetime D , as well as F are notisolated fixed points anymore. Additionally, MC coin-cides with D . In figure 4 are presented some orbits inthe phase space of the guiding system (65), (66), (67)for γ = Ω k = F and F .According to the center manifold analysis in ap-pendix 2 and supported by Figure 10 for γ = D is unstable (saddle type) for ( Σ + Ω k − ) (cid:54) =
0. However, if we restrict the analysis to ( Σ + Ω k − ) < D is asymptotically stable and behaves asa local attractor.In figure 5(a) are presented some orbits in the phasespace of guiding system (65), (66), (67) for γ = cor-responding to radiation. In figure 5(a) are presentedsome orbits in the phase space of guiding system (65),(66), (67) for γ = Ω k = F . For γ > D is locally asymptotically stable ac-cording to the center manifold analysis in 2. For γ = T , F , Q is invariant and unstable.3 k T Q DFF0 (a) γ = k T Q DFF0 (b) γ = . Figure 2: Phase space of the guiding system (65), (66), (67) for γ = , . k T Q DFF0 MC (a) γ = . k T Q DFF0 MC (b) γ = . Figure 3: Phase space of the guiding system (65), (66), (67) for some values of γ = . , . k T Q DFF0
Figure 4: Phase space of the guiding system (65), (66), (67) for γ = The results from the linear stability analysis, theCenter Manifold calculations in appendix B and com-bined with Theorem 7 lead to:
Theorem 8
The late time attractors of full system (A2) and averaged system (D2) for Bianchi III line elementare(i) The flat matter dominated FLRW universe F withds = − dt + a (cid:18) γ H t + (cid:19) γ dr a (cid:18) γ H t + (cid:19) γ (cid:2) d ϑ + ϑ d ζ (cid:3) (96) 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 matter-curvature scaling solution MC with ¯ Ω m = ( − γ ) andds = − dt + a (cid:18) γ H t + (cid:19) γ dr a (cid:18) γ H t + (cid:19) γ (cid:2) d ϑ + sinh ( ϑ ) d ζ (cid:3) (97) if < γ < . (iii) The Bianchi III flat spacetime D withds = − dt + c − t / dr + c − t / (cid:2) d ϑ + ϑ d ζ (cid:3) (98) if ≤ γ ≤ . k = − . In this case the time-averaged system is (62), (63),(64). Introducing the new time variable τ through d fd τ = H d fdt , we obtain the guiding system d ¯ Ω d τ = −
12 ¯ Ω (cid:16) ( γ − ) (cid:16) ¯ Ω − (cid:17) + ( γ − ) ¯ Ω k (cid:17) , (99) d ¯ Ω k d τ = − ¯ Ω k (cid:16) ( γ − ) ¯ Ω − γ + ( γ − ) ¯ Ω k + (cid:17) . (100) k T Q DFF0 (a) γ = . k T Q DFF0 (b) γ = Figure 5: Phase space of the guiding system (65), (66), (67) for γ = , F CF Ω k Ω γ = Figure 6: Phase plane of system (102) for γ = We have ¯ Ω + ¯ Ω k + ¯ Ω m =
1. Therefore, from the con-dition ¯ Ω m ≥ (cid:8) ( ¯ Ω , ¯ Ω k ) ∈ R : ¯ Ω + ¯ Ω k ≤ , ¯ Ω k ≥ (cid:9) . (101)For γ =
1, the guiding system (99)-(100) is reduced to d ¯ Ω d τ = −
12 ¯ Ω ¯ Ω k , d ¯ Ω k d τ = ¯ Ω k (cid:0) − ¯ Ω k (cid:1) . (102)The solution is¯ Ω ( τ ) = Ω √ e τ Ω k − Ω k + , ¯ Ω k ( τ ) = e τ Ω k e τ Ω k − Ω k + . (103)where ¯ Ω ( ) = Ω , ¯ Ω k ( ) = Ω k . C : ( ¯ Ω , ¯ Ω k ) = ( , ) with eigenvalues {− , − } isthe late-time attractor.The line of equilibrium points ¯ Ω ( Ω ) = Ω , ¯ Ω k ( Ω ) = Ω ∈ R with eigenvalues { , } is normally hyperbolic. Therefore, by consid-ering the eigenvalues with non zero real parts, it isfound that the line is unstable. This line contains thepoints F , F . In Figure 6 is presented the phase planeof system (102) for γ = Ω k = C .For γ (cid:54) =
1, the 2D guiding system (99), (100) has thefollowing equilibrium points:1. F : ( ¯ Ω , ¯ Ω k ) = ( , ) with eigenvalues (cid:110) ( γ − ) , γ − (cid:111) .i) It is a sink for 0 < γ < . ii) It is nonhyperbolic for γ = .iii) It is a saddle for < γ < < γ < ds = − dt + a (cid:18) γ H t + (cid:19) γ dr a (cid:18) γ H t + (cid:19) γ (cid:2) d ϑ + ϑ d ζ (cid:3) . (104)The corresponding solution is a flat matter dom-inated FLRW universe with ¯ Ω m = F : ( ¯ Ω , ¯ Ω k ) = ( , ) with eigenvalues { , − ( γ − ) } .i) It is a source for 0 < γ < < γ < ds = − dt + c − t / dr + c − t / (cid:2) d ϑ + ϑ d ζ (cid:3) . (105)Hence for large t the equilibrium point can beassociated with the flat Friedman Einstein-de-Sitter solution.3. C : ( ¯ Ω , ¯ Ω k ) = ( , ) with eigenvalues (cid:8) − , − γ (cid:9) .i) It is a saddle for 0 < γ < .ii) It is nonhyperbolic for γ = .iii) It is a sink for < γ < C we have q =
0. Then, ˙ H = − H ˙ a = aH = ⇒ H ( t ) = H H t + a ( t ) = a ( H t + ) . (106)The line element (3) becomes ds = − dt + a ( H t + ) dr a ( H t + ) (cid:2) d ϑ + sinh ( ϑ ) d ζ (cid:3) . (107)That is, C represents the Milne universe ( Ω = Ω m = , Ω k = , k = −
1) [188–190]. The met-ric for the Milne universe can be expressed withhyperspherical coordinates (3) also as ds = − dt + t ( d χ + sinh χ d Ω ) where d Ω = d ϑ + sin ϑ d ζ is the metric for a two-sphere8 F CF Ω k Ω γ = (a) F CF Ω k Ω γ = / (b) F CF Ω k Ω γ = (c) F CF Ω k Ω γ = / (d) Figure 7: Phase plane for system (99), (100) for different choices of γ . and χ = sinh − r is the curvature-corrected ra-dial component for negatively curved space thatvaries between 0 and + ∞ .In Figure 7 is portrayed the phase plane of the system(99), (100) for γ = , , . γ = . The results from the linear stability analysis com-bined with Theorem 7 (for Σ =
0) lead to:
Theorem 9
The late time attractors of the full system (A5) and the averaged system are:(i) The flat matter dominated FLRW universe F withds = − dt + a (cid:18) γ H t + (cid:19) γ dr a (cid:18) γ H t + (cid:19) γ (cid:2) d ϑ + ϑ d ζ (cid:3) (108) for < γ ≤ . F represents a quintessence fluid or a zero-acceleration model for γ = . In thelimit γ = we have a ( t ) = a (cid:16) γ H t + (cid:17) γ → a e H t , i.e., the de Sitter solution.(ii) The Milne solution C with ¯ Ω k = , k = − withds = − dt + a ( H t + ) dr a ( H t + ) (cid:2) d ϑ + sinh ( ϑ ) d ζ (cid:3) (109) for < γ < .
5. CONCLUSIONS
In this paper we have used Asymptotic methods andAveraging Theory to explore the solution’s space ofscalar field cosmologies with generalized harmonic po-tential V ( φ ) = µ (cid:104) b f (cid:16) − cos (cid:16) φ f (cid:17)(cid:17) + φ µ (cid:105) , b > (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 ) = (cid:18) f ( t , x ) (cid:19) is the first order term and (cid:18) f [ ] ( t , x ) (cid:19) isthe second-order remainder of the Taylor series of theleft hand side of (cid:18) ˙ H ˙ x (cid:19) = f ( t , x , H ) at H = Z = − Σ − Ω − Ω k we have found˙ Z = − HZ (cid:16) ( γ − ) Σ + ( γ − ) Ω + ( γ − ) Ω k − Ω cos ( ( Φ − t ω )) (cid:17) + O (cid:0) H (cid:1) which implies that the sign of 1 − Σ − Ω − Ω k is in-variant as H →
0. Then, from the equation (47) it wasproved that H is a monotonic decreasing function of t if 0 < Ω + Σ + Ω k <
1. This implies lim t → ∞ H ( t ) = t .Hence, from Theorem 7 was deduced that Ω , Σ , Ω k , and Φ evolve according to the time-averaged system (50),(51), (52), (53) as H →
0. Then can be established thestability of a periodic solution as it matches exactly thestability of the stationary solution of the time-averagedsystem. We have formulated theorems 8 and 9 aboutthe late-time behavior of our model, whose proofs arebased on the previous theorem, center manifold calcu-lations and linear stability analysis.In particular, in LRS Bianchi III the late time attrac-tors of full system (A2) and averaged system (D2) forBianchi III line element are: (i) The flat matter dominated FLRW universe F with ds = − dt + a (cid:18) γ H t + (cid:19) γ dr a (cid:18) γ H t + (cid:19) γ (cid:2) d ϑ + ϑ d ζ (cid:3) 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 ,i.e., the de Sitter solution.(ii) The matter-curvature scaling solution CS with¯ Ω m = ( − γ ) and ds = − dt + a (cid:18) γ H t + (cid:19) γ dr a (cid:18) γ H t + (cid:19) γ (cid:2) d ϑ + sinh ( ϑ ) d ζ (cid:3) if < γ < . (iii) The Bianchi III flat spacetime D with ds = − dt + c − t / dr + c − t / (cid:2) d ϑ + ϑ d ζ (cid:3) if 1 < γ ≤ k = −
1, the late time attrac-tors of the full system (A5) and the averaged systemare:(i) The flat matter dominated FLRW universe F with ds = − dt + a (cid:18) γ H t + (cid:19) γ dr a (cid:18) γ H t + (cid:19) γ (cid:104) d ϑ + k − sin ( √ k ϑ ) d ζ (cid:105) 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 Milne solution C with ds = − dt + a ( H t + ) dr a ( H t + ) (cid:2) d ϑ + sinh ( ϑ ) d ζ (cid:3) for < γ < ≤ γ ≤ (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 matterdominated flat FLRW universe represents quintessencefluid if 0 < γ < .Continuing the Program “Averaging GeneralizedScalar Field Cosmologies”, the cases: (II) Bianchi Iand Flat FLRW model and (III) Kantowski-Sachs andClosed FLRW are respectively studied in two compan-ion papers [146, 147]. In [146] it is used the sameapproach used here. That is, the oscillations enteringthe full system through Klein-Gordon equation can becontrolled and smoothed out when the Hubble factor H is used as a time-dependent parameter since it tendsmonotonically to zero and preserves it sign during theevolution. However, in the case (III) this approach isnot valid given that H is not necessarily monotonicallydecreasing to zero, and it can change its sign. There-fore, we have developed an alternative procedure in[147].Finally, our analytical results were strongly sup-ported by numerics in appendix D. We showed thatasymptotic methods and averaging theory are powerfultools to investigate scalar field cosmologies with gen-eralized harmonic potential; which have evident advan-tages, e.g., to determine the stability of the full oscilla-tion it is not needed to analyze the full dynamics, butonly the late-time behavior of the time-averaged (sim-pler) system has to be analyzed. Interestingly, for LRSBianchi III and Open FLRW model, when the back-ground fluid corresponds to a cosmological constant, H tends asymptotically to constant values depending onthe initial conditions which is consistent to de Sitter ex-pansion (see Figures 13(a) and 21(a)). In addition, forOpen FLRW and any γ < and Ω k > Ω k →
0. Onthe other hand, when γ > and Ω k > Ω k → 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 improving the English.
Appendix A: Main equations
In this appendix we use an amplitude-phase transfor-mation and some dimensionless variables to provide atime-averaged system in the sense of section 3.For Bianchi III metric are defined Ω = (cid:115) ω φ + ˙ φ H , Σ = σ + H , Ω k = K H , Φ = t ω − tan − (cid:18) ωφ ˙ φ (cid:19) . (A1)Open FLRW metric corresponds to Σ = Ω = − b µ √ H cos ( t ω − Φ ) sin (cid:32) √ H sin ( t ω − Φ ) Ω f ω (cid:33) − b f γ Ω µ H sin (cid:113) H sin ( t ω − Φ ) Ω f ω + (cid:0) ω − µ (cid:1) Ω ω sin ( t ω − Φ ) + H (cid:18)(cid:18) γ (cid:18) µ ω − (cid:19) + (cid:19) Ω − Ω (cid:19) cos ( t ω − Φ )+ H Ω (cid:16) − ( γ − ) Σ − γ µ Ω ω + γ + ( − γ ) Ω k (cid:17) , (A2a)˙ Σ = − b f γ Σ µ H sin (cid:113) H sin ( t ω − Φ ) Ω f ω − ( γ − ) H ΣΩ cos ( t ω − Φ )+ H (cid:32) Ω k − ( γ − ) Σ + (cid:16) − γ µ Ω ω sin ( t ω − Φ ) + γ − γ Ω k + Ω k − (cid:17) Σ (cid:33) , (A2b)1˙ Ω k = b f γ Ω k µ H (cid:32) − + cos (cid:32) √ H sin ( t ω − Φ ) Ω f ω (cid:33)(cid:33) − ( γ − ) H Ω Ω k cos ( t ω − Φ )+ H (cid:16) − ( γ − ) Σ − Σ − ( γ − )( Ω k − ) − γ µ Ω ω sin ( t ω − Φ ) (cid:17) Ω k , (A2c)˙ Φ = − b µ √ H Ω sin ( t ω − Φ ) sin (cid:32) √ H sin ( t ω − Φ ) Ω f ω (cid:33) + (cid:0) ω − µ (cid:1) ω sin ( t ω − Φ ) − ( t ω − Φ ) H sin ( t ω − Φ ) , (A2d)˙ H = − ( + q ) H , (A2e)where the deceleration parameter is given by q = − b f γ µ H sin (cid:113) H sin ( t ω − Φ ) Ω f ω − ( γ − ) cos ( t ω − Φ ) Ω − ( γ − ) Σ − γ µ Ω ω sin ( t ω − Φ ) − ( γ − )( Ω k − ) . (A3)With f = b µ ω − µ > = Ω + Ω k + Σ + ρ m H + H ( ω − µ ) U ( φ ) , where the function defined by U ( φ ) = f (cid:16) − cos (cid:16) φ f (cid:17)(cid:17) − φ . It satisfies U (cid:48) ( ) = U (cid:48)(cid:48) ( ) = U (cid:48)(cid:48)(cid:48) ( ) = U ( iv ) ( ) = − f , which implies φ = U ( φ ) . Therefore, U ( φ ) ≤
0. This implies that for ω > µ , Ω , Σ and Ω m : = ρ m H can be greater than 1. Ω m : = ρ m H is parameterized by the equation Ω m = − Ω k − Σ − Ω + Ω (cid:18) − µ ω (cid:19) sin ( Φ − t ω ) + b µ sin (cid:32) (cid:113) H Ω ( µ − ω ) sin ( Φ − t ω ) b µ ω (cid:33) H ( µ − ω )= − Ω k − Σ − Ω + O ( H ) . (A4)For open FLRW metric, where Σ ≡
0, the full system is˙ Ω = − b γ f µ Ω H sin (cid:113) H Ω sin ( t ω − Φ ) f ω − b µ √ H cos ( t ω − Φ ) sin (cid:32) √ H Ω sin ( t ω − Φ ) f ω (cid:33) + H cos ( t ω − Φ ) (cid:18) Ω (cid:18) γ (cid:18) µ ω − (cid:19) + (cid:19) − Ω (cid:19) + H Ω (cid:18) γ − γ µ Ω ω + ( − γ ) Ω k (cid:19) + (cid:0) ω − µ (cid:1) Ω sin ( t ω − Φ ) ω , (A5a)˙ Ω k = b γ f µ Ω k H (cid:32) − + cos (cid:32) √ H Ω sin ( t ω − Φ ) f ω (cid:33)(cid:33) − H Ω k (cid:18) γ µ Ω sin ( t ω − Φ ) ω + ( γ − )( Ω k − ) (cid:19) − ( γ − ) H Ω Ω k cos ( t ω − Φ ) , (A5b)˙ Φ = − H sin ( t ω − Φ ) cos ( t ω − Φ ) + (cid:0) ω − µ (cid:1) ω sin ( t ω − Φ ) − b µ √ H Ω sin ( t ω − Φ ) sin (cid:18) √ H Ω sin ( t ω − Φ ) f ω (cid:19) , (A5c)˙ H = − ( + q ) H , (A5d)2with deceleration parameter q = − + γ − γ µ Ω sin ( t ω − Φ ) ω − ( γ − ) Ω cos ( t ω − Φ ) − γ Ω k + Ω k − b γ f µ H sin √ H Ω sin ( t ω − Φ ) f ω . (A5e)For negatively curved FLRW models the fractional energy density Ω m : = ρ m H is parameterized by Ω m = − Ω k − Ω + Ω (cid:18) − µ ω (cid:19) sin ( Φ − t ω ) + b µ H ( µ − ω ) sin √ H Ω ( µ − ω ) sin ( Φ − t ω ) b µ ω = − Ω k − Ω + O ( H ) . (A6) Appendix B: Center Manifold Calculations1. Center manifold of F for γ = Let J be the linearization matrix of the guiding sys-tem (65), (66), (67). Letting γ = , the Jordan Decom-position of J in the critical point F = ( , , ) is J ( D ) = s − diag {− , − / , } s , s = . Defining the new variables x = ¯ Σ − ¯ Ω k , y = ¯ Ω , z = ¯ Ω k , (B1)to obtain the new equations x (cid:48) = Ax + f ( x , y , z ) , (B2) y (cid:48) = By + f ( x , y , z ) , (B3) z (cid:48) = Cz + f ( x , y , z ) , (B4)with A = − , B = − , C = , (B5) f ( x , y , z ) = x + x z + xy − xz + xz − y z − z + z , (B6) f ( x , y , z ) = x y + xyz + y + yz , (B7) f ( x , y , z ) = x z + xz − xz + y z + z − z . (B8) F - - - zz Figure 8: One dimensional flow for (B16) for z ∈ [ − , ] theorigin is stable if z ≥
0. Note that z ≥ Ω k ≥ Note that if J is the linearization matrix of system(B2), (B3), (B4) then the eigensystem of J ( , , ) is (cid:18) − − { , , } { , , } { , , } (cid:19) . This implies that the local center manifold of the originfor (B2), (B3), (B4) is given by the graph W cloc ( ) = (cid:40) ( x , y , z ) ∈ R : x = h ( z ) , y = h ( z ) , h ( ) = h ( ) = , h (cid:48) ( ) = h (cid:48) ( ) = , | z | < δ (cid:41) (B9)for some δ >
0. Therefore, we can use Taylor series todefine h ( z ) = a z + a z + a z + a z + a z + O (cid:0) z (cid:1) , (B10) h ( z ) = b z + b z + b z + b z + b z + O (cid:0) z (cid:1) . (B11)3The following quasi-linear differential equations N ( h ( z )) ≡ h (cid:48) ( z )( Cz + f ( h ( z ) , h ( z ) , z )) − Ah ( z ) − f ( h ( z ) , h ( z ) , z ) , (B12) N ( h ( z )) ≡ h (cid:48) ( z )( Cz + f ( h ( z ) , h ( z ) , z )) − Bh ( z ) − f ( h ( z ) , h ( z ) , z ) , (B13)must be solved for the a i , b i , i = , . . . b i , i = , . . . a = , a = , a = , a = .Hence, we obtain h ( z ) = z + z + z + z + z + O ( z ) , (B14) h ( z ) = O ( z ) , (B15)and the dynamics in the center manifold is given by: z (cid:48) = − z + z + z + z + z + O (cid:0) z (cid:1) . (B16)Then, F is locally asymptotically stable for z ≥
0. Notethat z ≥ Ω k ≥ z = ¯ Ω k . In Figure 8 is represented the one dimensional flowfor (B16) for z ∈ [ − , ] . The origin is locally asymp-totically stable if z ≥
2. Center manifold of D for γ ≥ Let J be the linearization matrix of the guiding sys-tem (65), (66), (67). Letting γ / ∈ (cid:8) , (cid:9) , the JordanDecomposition of J in the critical point D is J ( D ) = s − diag {− / , , − γ } s , s = − γ − γ − . Defining the new variables x = γ (
12 ¯ Σ − Ω k − ) −
16 ¯ Σ + Ω k + − γ , (B17) y = ( γ − )( Σ + Ω k − ) ( γ − ) , (B18) z = ¯ Ω , (B19)to obtain the new equations x (cid:48) = Ax + f ( x , y , z , γ ) , (B20) y (cid:48) = By + f ( x , y , z , γ ) , (B21) z (cid:48) = Cz + f ( x , y , z , γ ) , (B22) with A = − , B = ( − γ ) , C = , (B23) f ( x , y , z , γ ) = ( γ − )( γ − ) x − γ + x (cid:18) − γ γ − + ( γ − )( γ − ) y ( γ − )( γ − ) (cid:19) + x (cid:32) γ ( γ − ) y ( − γ ) ( γ − )+ (cid:18) − γ + − γ + γ − + (cid:19) y (cid:33) + z (cid:32) ( γ − ) γ − + ( γ − )( γ − ) x − γ + ( γ − )( γ − ) y − γ (cid:33) − ( γ − ) y ( − γ ) ( γ − ) − ( γ ( γ ( γ − ) + ) − )( γ − ) y ( − γ ) ( γ − ) , (B24) f ( x , y , z , γ ) = ( γ − )( γ − ) x − γ + x (cid:32) (cid:18) − γ + γ − + (cid:19) + ( γ − )( γ − ) y − γ (cid:33) + x (cid:18) ( γ − )( γ − ) y ( γ − )( γ − ) + ( − γ ) y γ − (cid:19) + z (cid:32) ( γ − )( γ − ) − γ + ( γ − )( γ − ) x − γ + ( γ − )( γ − ) y − γ (cid:33) − ( γ − ) ( γ − ) y ( − γ ) ( γ − )+ (cid:18) − γ + − γ + γ − + (cid:19) y − γ y , (B25) f ( x , y , z , γ ) = z (cid:32) (cid:18) − γ (cid:19) x + x (cid:18) ( γ − ) y γ − − (cid:19) − ( γ − ) y ( − γ ) + (cid:18) − γ + − γ + (cid:19) y (cid:33) − ( γ − ) z . (B26)Note that if J is the linearization matrix of system(B20), (B21), (B22) then the eigensystem of J ( , , ) is (cid:18) − − γ ( , , ) ( , , ) ( , , ) (cid:19) . W cloc ( ) = (cid:40) ( x , y , z ) ∈ R : x = h ( z ) , y = h ( z ) , h ( ) = h ( ) = , h (cid:48) ( ) = h (cid:48) ( ) = , | z | < δ (cid:41) (B27)for some δ > h ( z ) = a z + a z + a z + a z + a z + O ( z ) , (B28) h ( z ) = b z + b z + b b + b z + b z + O ( z ) . (B29)The following quasilinear differential equations N ( h ( z )) ≡ h (cid:48) ( z )( Cz + f ( h ( z ) , h ( z ) , z , γ )) − Ah ( z ) − f ( h ( z ) , h ( z ) , z , γ ) , (B30) N ( h ( z )) ≡ h (cid:48) ( z )( Cz + f ( h ( z ) , h ( z ) , z , γ )) − Bh ( z ) − f ( h ( z ) , h ( z ) , z , γ ) , (B31)must be solved for the a i , b i up to order six to ap-proximate the center manifold. We notice that theeven terms are zero and a = γ − γ − , a = γ − γ − , a = γ − γ − , b = − γ ( γ − ) , b = − γ ( γ − ) , b = − γ ( γ − ) .Hence, we obtain h ( z ) = ( γ − ) z γ − + ( γ − ) z γ − + ( γ − ) z γ − , (B32) h ( z ) = ( − γ ) z ( γ − ) + ( − γ ) z ( γ − ) + ( − γ ) z ( γ − ) . (B33)Letting γ = the Jordan Decomposition of J in thecritical point D is J ( D ) = s − diag {− / , − , } s , s = − . Defining the new variables x = ¯ Ω k − , y = ¯ Σ + ¯ Ω k − , z = ¯ Ω , (B34) to obtain the new equations x (cid:48) = Ax + f ( x , y , z ) , (B35) y (cid:48) = By + f ( x , y , z ) , (B36) z (cid:48) = Cz + f ( x , y , z ) , (B37)where A = − , B = − , C = , (B38) f ( x , y , z ) = x − x y − x + xy − xy − xz + y − z , (B39) f ( x , y , z ) = x − x y + x − xy − xy − xz + y + y − yz − z , (B40) f ( x , y , z ) = x z − xyz − xz + y z + yz − z . (B41)Note that if J is the linearization matrix of system(B35), (B36), (B37) then the eigensystem of J ( , , ) is (cid:18) − − { , , } { , , } { , , } (cid:19) . Then, using Taylor series, the local center manifold ofthe origin for (B35), (B36), (B37) is given by the graph(B27) where h ( z ) = − z − z − z , (B42) h ( z ) = − z − z − z . (B43)Letting γ = the Jordan Decomposition of J in thecritical point D is J ( D ) = s − − −
00 0 0 s , s = − −
01 0 0 . Defining the new variables x = ¯ Ω k − , (B44) y = − ( Σ + Ω k − ) , (B45) z = ¯ Ω , (B46)5to obtain the new equations x (cid:48) = A x + A y + f ( x , y , z ) , (B47) y (cid:48) = By + f ( x , y , z ) , (B48) z (cid:48) = Cz + f ( x , y , z ) , (B49)where A = − , A = , B = − , C = , (B50) f ( x , y , z ) = x + x y − x + xy + xy − xz + y − z , (B51) f ( x , y , z ) = − x − x y − x + xy − xy + xz + y − y − yz + z , (B52) f ( x , y , z ) = x z + xyz − xz + y z − yz − z . (B53)Note that if J is the linearization matrix of system(B47), (B48), (B49) then the eigensystem of J ( , , ) is (cid:18) − − { , , } { , , } { , , } (cid:19) . Then, using Taylor series, the local center manifold ofthe origin is(B47), (B48), (B49) is given by the graph(B27) where h ( z ) = − z − z − z , (B54) h ( z ) = z + z + z . (B55)Finally, in these three cases (which corresponds to γ > z (cid:48) = − z + z + O ( z ) . (B56)This is a gradient-like equation z (cid:48) = − ∇ U ( z ) , U ( z ) = z − z , (B57)for which the origin is a degenerate minimum of sec-ond order, i.e., U (cid:48) ( ) = U (cid:48)(cid:48) ( ) = U (cid:48)(cid:48)(cid:48) ( ) = , U ( iv ) ( ) = >
0. Therefore, the origin is locally asymptoticallystable. Indeed, performing a geometrical analysis ofthe flow over (B56) we note that the origin is locallyasymptotically stable as shown in Figure 9. D - - - - zz Figure 9: One dimensional flow for (B56) for z ∈ [ − , ] . Theorigin is stable. Letting γ =
1, we deduce that the center manifoldof D is 2-dimensional. On the other hand, the sys-tem admits the lines of equilibrium points ( ¯ Ω , ¯ Σ , ¯ Ω k ) =( ¯ Ω ∗ , , ) and ( ¯ Ω , ¯ Σ , ¯ Ω k ) = ( ¯ Ω ∗ , , ) where ¯ Ω ∗ is anarbitrary number which satisfies ¯ Ω ∗ ∈ [ , ] . Therefore, D as well as F are not isolated fixed points anymore.To analyze the stability of D : ( ¯ Ω , ¯ Σ , ¯ Ω k ) = ( , , ) first notice that the Jordan decomposition of the Jaco-bian matrix evaluated at D is J ( D ) = s − − s , s = − . Defining the new variables x = ( Σ − Ω k + ) , y = ( Σ + Ω k − ) , z = ¯ Ω , (B58)we obtain the new equations x (cid:48) = Ax + f ( x , y , z ) , (B59) y (cid:48) = C y + f ( x , y , z ) , (B60) z (cid:48) = C z + f ( x , y , z ) , (B61)where A = − , C = , C = , (B62) f ( x , y , z ) = x + x (cid:18) y + (cid:19) + x (cid:18) y − y (cid:19) − y − y , (B63)6 D - - - - y z Figure 10: Two dimensional flow for (B71), (B72). The phys-ical region of the phase space is y : = ( Σ + Ω k − ) < z : = ¯ Ω > f ( x , y , z ) = − x + x (cid:18) y + (cid:19) + x (cid:18) y + y (cid:19) + y + y , (B64) f ( x , y , z ) = z (cid:18) x + x ( y + ) + y + y (cid:19) . (B65)The local center manifold of the origin for (B59),(B60), (B61) is given by the graph W cloc ( ) = (cid:40) ( x , y , z ) ∈ R : x = h ( y , z ) , h ( , ) = , ∂ h ∂ y ( , ) = ∂ h ∂ z ( , ) = , | ( y , z ) T | < δ (cid:41) , (B66) for some δ > h ( y , z ) satisfies the quasilinear partial differentialequation (cid:32) h − ( y + ) h − y ( y + ) h − y ( y + ) (cid:33) ∂ h ∂ y + (cid:18) − zh − ( y + ) zh − y ( y + ) z (cid:19) ∂ h ∂ z + h + ( y + ) h + ( y ( y − ) − ) h − y ( y + ) = . (B67) We propose the Taylor expansion h ( y , z ) = a y + a yz + a z + b y + b y z + b yz + b z + O ( ) , (B68)where O ( ) denotes terms of fourth order in the vectornorm. Therefore, equation (B67) can be expressed as14 ( − a − ) y + y z (cid:18) − a − b (cid:19) − a yz − yz ( a + b ) − a z − b z = O ( ) . (B69)Equating the terms of the same power in y , z we havea solution a = − , a = , a = , b arbitrary , b = , b = , b =
0. Then, we obtain h ( y , z ) = − y + b y + O ( ) . (B70)The dynamics at the center manifold is given by y (cid:48) = y + y , (B71) z (cid:48) = yz + y z . (B72)In Figure 10 is presented a two dimensional flow for(B71) and (B72) where is shown that the origin is un-stable (saddle type) for y (cid:54) =
0. However, if we restrictthe analysis to y < D is asymptotically stable andbehaves as a local attractor.To analyze the stability of an arbitrary point D ∗ : ( ¯ Ω , ¯ Σ , ¯ Ω k ) = ( ¯ Ω ∗ , , ) with ¯ Ω ∗ (cid:54) = D ∗ is J ( D ∗ ) = s − − s , s = Ω ∗ − Ω ∗ Ω ∗ . Defining the new variables x = ( Σ − Ω k + ) (B73) y = ( ¯ Ω ∗ ( − Σ + Ω k − ) + Ω ) (B74) z =
18 ¯ Ω ∗ ( Σ + Ω k − ) , (B75)7 D * - - - - Y Z Figure 11: Two dimensional flow for (B90), (B91). The ori-gin is unstable. we obtain the new equations x (cid:48) = Ax + f ( x , y , z ) , (B76) y (cid:48) = C y + B x + B z + f ( x , y , z ) , (B77) z (cid:48) = C z + f ( x , y , z ) , (B78)where A = − , B = Ω ∗ , B = , C = , C = , (B79) f ( x , y , z ) = x + z (cid:18) x Ω ∗ − x ¯ Ω ∗ (cid:19) + x + z (cid:18) x Ω ∗ −
14 ¯ Ω ∗ (cid:19) − z Ω ∗ , (B80) f ( x , y , z ) = − x ¯ Ω ∗ + x (cid:18) y + Ω ∗ (cid:19) + z (cid:18) − x + x (cid:18) y ¯ Ω ∗ + (cid:19) + y ¯ Ω ∗ (cid:19) + z (cid:18) x ¯ Ω ∗ + y Ω ∗ +
116 ¯ Ω ∗ (cid:19) + xy + z ¯ Ω ∗ , (B81) f ( x , y , z ) = − x ¯ Ω ∗ + x ¯ Ω ∗ + (cid:18) x + x (cid:19) z + z (cid:18) x Ω ∗ + Ω ∗ (cid:19) + z Ω ∗ . (B82)The local center manifold of the origin for (B76), (B77), (B78) is given by the graph W cloc ( ) = (cid:40) ( x , y , z ) ∈ R : x = h ( y , z ) , h ( , ) = , ∂ h ∂ y ( , ) = ∂ h ∂ z ( , ) = , | ( y , z ) T | < δ (cid:41) , (B83) for some δ > h ( y , z ) satisfies the quasilinear partialdifferential equation ∂ h ∂ z (cid:32)
34 ¯ Ω ∗ h − (cid:0) Ω ∗ + z (cid:1) h − z (cid:0)
22 ¯ Ω ∗ + z (cid:1) h Ω ∗ − z (cid:0) Ω ∗ + z (cid:1) Ω ∗ (cid:33) + ∂ h ∂ y (cid:32) ¯ Ω ∗ h + (cid:18) − Ω − y + z (cid:19) h − (cid:0) y ¯ Ω ∗ +
12 ¯ Ω ∗ + z ¯ Ω ∗ + yz + z (cid:1) h Ω − z (cid:0) y ¯ Ω ∗ + Ω ∗ + z ¯ Ω ∗ + yz + z (cid:1) Ω ∗ (cid:33) + (cid:0) ¯ Ω ∗ + z (cid:1) h Ω ∗ + (cid:32) z (cid:0) z − Ω (cid:1) ¯ Ω ∗ − (cid:33) h − z (cid:0) ¯ Ω ∗ + z (cid:1) Ω ∗ + h = . (B84)We propose the Taylor expansion h ( y , z ) = a y + a yz + a z + b y + b y z + b yz + b z + O ( ) , (B85)where O ( ) denotes terms of fourth order in the vectornorm. Therefore, equation (B84) can be expressed as yz (cid:18) − Ω ∗ (cid:16) a a + a (cid:17) − a + a Ω ∗ − b − b (cid:19) + y z (cid:18) − a a ¯ Ω ∗ − a ¯ Ω ∗ − b − b (cid:19) + yz (cid:18) − a − a (cid:19) − a y − z (cid:0) Ω ∗ (cid:0) a (cid:0) a ¯ Ω ∗ + (cid:1) + a + b ¯ Ω ∗ + b ¯ Ω ∗ (cid:1) + (cid:1)
12 ¯ Ω ∗ + z (cid:18) − a − a − Ω (cid:19) = O ( ) . (B86) Equating the terms of the same power in y , z we havea solution a = , a = , a = −
16 ¯ Ω ∗ , b = b , b = − b , b = (cid:16) Ω ∗ − b (cid:17) . For simplicity, we set8 b =
112 ¯ Ω ∗ , b =
0. Then, we obtain h ( y , z ) = y
32 ¯ Ω ∗ − y z
16 ¯ Ω ∗ + yz
12 ¯ Ω ∗ − z Ω ∗ + O ( ) . (B87)The dynamics at the center manifold is given by y (cid:48) = z + yz ¯ Ω ∗ + z Ω ∗ + y Ω ∗ − y z Ω ∗ + yz Ω ∗ + z
18 ¯ Ω ∗ , (B88) z (cid:48) = z Ω ∗ + z ¯ Ω ∗ . (B89)Using the re-scaling ( T , Y , Z ) = ( ¯ Ω ∗ τ , y ¯ Ω ∗ , z ¯ Ω ∗ ) , ¯ Ω ∗ > dYdT = Z + Y Z + Z + Y − Y Z + Y Z + Z , (B90) dZdT = Z + Z Z (cid:54) = Appendix C: Proof of Theorem 7
Lemma 10 (Gronwall’s Lemma) Differentialform .i) Let be η ( t ) a nonnegative function, absolutelycontinuous over [ , T ] which satisfies almost every-where the differential inequality η (cid:48) ( t ) ≤ φ ( t ) η ( t ) + ψ ( t ) , where φ ( t ) and ψ ( t ) are nonnegative summable func-tions over [ , T ] . Then, η ( t ) ≤ e (cid:82) t φ ( s ) ds (cid:20) η ( ) + (cid:90) t ψ ( s ) ds (cid:21) , is satisfied for all ≤ t ≤ T .ii) In particular, if η (cid:48) ( t ) ≤ φ ( t ) η ( t ) , t ∈ [ , T ] , η ( ) = , then η ≡ , t ∈ [ , T ] . 2. Integral form i) Let be ξ ( t ) a nonnegative function, summable over [ , T ] which satisfies almost everywhere the integral in-equality ξ ( t ) ≤ C (cid:90) t ξ ( s ) ds + C , C , C ≥ . Then, ξ ( t ) ≤ C ( + C te C t ) , almost everywhere for t in ≤ t ≤ T . ii) In particular,if ξ ( t ) ≤ C (cid:90) t ξ ( s ) ds , C ≥ , almost everywhere for t in ≤ t ≤ T , then ξ ≡ , al-most everywhere for t in ≤ t ≤ T .
Theorem 7 states:
Let be defined the functions ¯ Ω , ¯ Σ , ¯ Ω k , ¯ Φ , and H which satisfy the time-averaged sys-tem (50) , (51) , (52) , (53) , (54) . Then, there exists func-tions g , g , g and g , such that if Ω = Ω + Hg ( H , Ω , Σ , Ω k , Φ , t ) , Σ = Σ + Hg ( H , Ω , Σ , Ω k , Φ , t ) , Ω k = Ω k + Hg ( H , Ω , Σ , Ω k , Φ , t ) , Φ = Φ + Hg ( H , Ω , Σ , Ω k , Φ , t ) , (C1) the functions Ω , Σ , Ω k , Φ and ¯ Ω , ¯ Σ , ¯ Ω k , ¯ Φ have thesame limit as t → ∞ . Setting Σ = Σ = are derivedanalogous results for negatively curved FLRW model. Proof . Defining Z = − Σ − Ω − Ω k it followsfrom˙ Z = − HZ (cid:16) ( γ − ) Σ + ( γ − ) Ω + ( γ − ) Ω k − Ω cos ( ( Φ − t ω )) (cid:17) + O (cid:0) H (cid:1) that the sign of 1 − Σ − Ω − Ω k is invariant as H → H is a monotonicdecreasing function of t if 0 < Ω + Σ + Ω k <
1. Theseallow to define recursively the bootstrapping sequences t = t ∗ H = H ( t ∗ ) , t n + = t n + H n H n + = H ( t n + ) , (C2)such that lim n → ∞ H n = n → ∞ t n = ∞ .9Giving the expansions (C1) we obtain˙ Ω = (cid:32) Ω (cid:16) Σ + Ω + Ω k − γ (cid:0) Σ + Ω + Ω k − (cid:1) + (cid:0) Ω − (cid:1) cos ( ( t ω − Φ )) − (cid:17)(cid:33) H − ∂ g ∂ t H + O (cid:0) H (cid:1) , ˙ Σ = (cid:32) Σ cos ( ( t ω − Φ )) Ω + Ω k + Σ (cid:16) Σ + Ω + Ω k − − γ (cid:0) Σ + Ω + Ω k − (cid:1) (cid:17)(cid:33) H − ∂ g ∂ t H + O (cid:0) H (cid:1) , ˙ Ω k = (cid:32) Ω k (cid:16) Σ − Σ + Ω + Ω k − γ (cid:0) Σ + Ω + Ω k − (cid:1) + Ω cos ( ( t ω − Φ )) − (cid:17)(cid:33) H − ∂ g ∂ t H + O (cid:0) H (cid:1) , ˙ Φ = − (cid:18)
32 sin ( ( t ω − Φ )) + ∂ g ∂ t (cid:19) H + O (cid:0) H (cid:1) . Let be defined ∆Ω = Ω − ¯ Ω , ∆Σ = Σ − ¯ Σ , ∆Ω k = Ω k − ¯ Ω k , ∆Φ = Φ − ¯ Φ , such that Ω ( t n ) = ¯ Ω ( t n ) = Ω n , Σ ( t n ) = ¯ Σ ( t n ) = Σ n , Ω k ( t n ) = ¯ Ω k ( t n ) = Ω kn , Φ ( t n ) = ¯ Φ ( t n ) = Φ n , with 0 < Ω ( t n ) + Σ ( t n ) + Ω k ( t n ) < g i such that ∂ g ∂ t = Ω (cid:32) − γ (cid:0) Σ + Ω + Ω k − (cid:1) + Σ + Ω + Ω k − (cid:33) + Ω (cid:0) Ω − (cid:1) cos ( ( Φ − t ω ))+ ¯ Ω (cid:18) ( γ − ) ¯ Σ + ( γ − ) ¯ Ω k (cid:19) + ( γ − ) (cid:0) ¯ Ω − ¯ Ω (cid:1) + Ω ¯ Σ + Ω ¯ Ω + Ω ¯ Ω k , (C3a) ∂ g ∂ t = (cid:32) Σ (cid:16) − γ (cid:0) Σ + Ω + Ω k − (cid:1) + Σ + Ω + Ω k − (cid:17) + Ω k (cid:33) + Σ Ω cos ( ( Φ − t ω ))+ ¯ Ω (cid:18) ( γ − ) ¯ Σ + Σ (cid:19) + ¯ Ω k (cid:18) ( γ − ) ¯ Σ + Σ − (cid:19) + ( γ − ) ¯ Σ − ( γ − ) ¯ Σ + Σ ¯ Σ , (C3b) ∂ g ∂ t = Ω k (cid:32) − γ (cid:0) Σ + Ω + Ω k − (cid:1) + Σ ( Σ − ) + Ω + Ω k − (cid:33) + Ω Ω k cos ( ( Φ − t ω ))+ ¯ Ω k (cid:32) (cid:18) γ − (cid:19) ¯ Σ + Σ + ( − γ + Ω k + ) (cid:33) + ¯ Ω (cid:18)(cid:18) γ − (cid:19) ¯ Ω k + Ω k (cid:19) + (cid:18) γ − (cid:19) ¯ Ω k + Ω k ¯ Σ , (C3c) ∂ g ∂ t = −
32 sin ( ( Φ − t ω )) , (C3d)0where we have substituted the explicit dependence of t of the averaged functions ¯ Ω , ¯ Σ , ¯ Ω k .By integration of (C3) we obtain the explicit expressions of the g i , i = , . . . , g ( H , Ω , Σ , Ω k , Φ , t ) = Ω (cid:32) − γ (cid:0) Σ + Ω + Ω k − (cid:1) + Σ + Ω + Ω k − (cid:33) t + Ω (cid:0) − Ω (cid:1) sin ( ( Φ − t ω )) ω + (cid:90) (cid:34) ¯ Ω ( t ) (cid:18) ( γ − ) ¯ Σ ( t ) + ( γ − ) ¯ Ω k ( t ) (cid:19) + ( γ − ) (cid:0) ¯ Ω ( t ) − ¯ Ω ( t ) (cid:1) + Ω ¯ Σ ( t ) + Ω ¯ Ω ( t ) + Ω ¯ Ω k ( t ) (cid:35) dt + C ( Ω , Σ , Ω k , Φ ) , (C4) g ( H , Ω , Σ , Ω k , Φ , t ) = (cid:32) Σ (cid:16) − γ (cid:0) Σ + Ω + Ω k − (cid:1) + Σ + Ω + Ω k − (cid:17) + Ω k (cid:33) t − Σ Ω sin ( ( Φ − t ω )) ω + (cid:90) (cid:34) ¯ Ω ( t ) (cid:18) ( γ − ) ¯ Σ ( t ) + Σ (cid:19) + ¯ Ω k ( t ) (cid:18) ( γ − ) ¯ Σ ( t ) + Σ − (cid:19) + ( γ − ) ¯ Σ ( t ) − ( γ − ) ¯ Σ ( t ) + Σ ¯ Σ ( t ) (cid:35) dt + C ( Ω , Σ , Ω k , Φ ) , (C5) g ( H , Ω , Σ , Ω k , Φ , t ) = Ω k (cid:32) − γ (cid:0) Σ + Ω + Ω k − (cid:1) + Σ ( Σ − ) + Ω + Ω k − (cid:33) t − Ω Ω k sin ( ( Φ − t ω )) ω + (cid:90) (cid:34) ¯ Ω k ( t ) (cid:32) (cid:18) γ − (cid:19) ¯ Σ ( t ) + Σ ( t ) + ( − γ + Ω k + ) (cid:33) + ¯ Ω ( t ) (cid:18)(cid:18) γ − (cid:19) ¯ Ω k ( t ) + Ω k (cid:19) + (cid:18) γ − (cid:19) ¯ Ω k ( t ) + Ω k ¯ Σ ( t ) (cid:35) dt + C ( Ω , Σ , Ω k , Φ ) , (C6) g ( H , Ω , Σ , Ω k , Φ , t ) = ( Φ − t ω ) ω + C ( Ω , Σ , Ω k , Φ ) , (C7)where C i , i = , , , H becausethe terms ∂ g i ∂ H ˙ H in the procedure contribute to O ( H ) -factors. The terms unbarred inside the integral are taken asconstants during the integration.Hence, we obtain ˙ ∆Ω = H ∆Ω (cid:0) − ¯ Σ − ¯ Ω − ¯ Ω k (cid:1) + O (cid:0) H (cid:1) , ˙ ∆Σ = H ∆Σ (cid:0) − ¯ Σ − ¯ Ω − ¯ Ω k (cid:1) + O (cid:0) H (cid:1) , ˙ ∆Ω k = H ∆Ω k (cid:0) − ¯ Σ − ¯ Ω − ¯ Ω k (cid:1) + O (cid:0) H (cid:1) , ∆Φ = H ω (cid:34) Ω (cid:0) µ − ω (cid:1) b µ + ω cos ( Φ − t ω ) × (cid:32) − γ (cid:0) Σ + Ω + Ω k − (cid:1) + Σ + (cid:0) Ω + (cid:1) cos ( ( Φ − t ω )) + Ω + Ω k (cid:33)(cid:35) + O (cid:0) H (cid:1) . Continuing the proof, let t ∈ [ t n , t n + ] such that t − t n ≤ t n + − t n = H n , we have the following: | ∆Ω ( 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:90) tt n (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:16) − ¯ Σ − ¯ Ω − ¯ Ω k (cid:17)(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 ) , | ∆Σ | = | Σ − ¯ Σ | = (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:90) tt n (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:16) − ¯ Σ − ¯ Ω − ¯ Ω k (cid:17)(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 ) , | ∆Ω k | = | Ω k − ¯ Ω k | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:90) tt n (cid:16) ˙ Ω k ( s ) − ˙¯ Ω k ( s ) (cid:17) 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:16) − ¯ Σ − ¯ Ω − ¯ Ω k (cid:17)(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 | ∆Ω k ( 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 H (cid:40) Ω (cid:0) µ − ω (cid:1) b µ + ω cos ( Φ − t ω ) × (cid:32) γ (cid:0) − Σ − Ω − Ω k (cid:1) + Σ + (cid:0) Ω + (cid:1) cos ( ( Φ − t ω )) + Ω + Ω k (cid:33)(cid:41) ds (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + O ( H ) ≤ M H n ( t − t n ) + O ( H n ) , whereby assuming C i ( Ω , Σ , Ω k , Φ ) , i = , . . . , g i , i = , . . . , t ∈ [ t n , t n + ] . By continuity of¯ Ω , Ω , ¯ Σ , Σ , Φ , ¯ Ω k , Ω k , the following constants L , M , M , M , M are finite: L = max t ∈ [ t n , t n + ] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:16) − ¯ Σ − ¯ Ω − ¯ Ω k (cid:17)(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:16) Σ (cid:0) ( γ − ) ¯ Ω k + Ω k (cid:1) + Ω (cid:0) ( γ − ) ¯ Ω k + Ω k (cid:1) + (cid:0) ¯ Ω k − (cid:1)(cid:0) ( γ − ) ¯ Ω k + Ω k (cid:1) + Σ ¯ Ω k (cid:17) ∂ g ∂ Ω k + (cid:16) Ω (cid:0) ( γ − ) ¯ Σ + Σ (cid:1) + ¯ Ω k (cid:0) ( γ − ) ¯ Σ + Σ − (cid:1) + (cid:16) ¯ Σ − (cid:17)(cid:0) ( γ − ) ¯ Σ + Σ (cid:1)(cid:17) ∂ g ∂ Σ + (cid:16) Σ (cid:0) ( γ − ) ¯ Ω + Ω (cid:1) + ¯ Ω k (cid:0) ( γ − ) ¯ Ω + Ω (cid:1) + (cid:16) ¯ Ω − (cid:17)(cid:0) ( γ − ) ¯ Ω + Ω (cid:1)(cid:17) ∂ g ∂ Ω + g (cid:16) − γ (cid:16) Σ + Ω + Ω k − (cid:17) + Σ + (cid:16) Ω − (cid:17) cos ( ( Φ − t ω )) + Ω + Ω k − (cid:17) − ( γ − ) Σ Ω g + (cid:18) Ω − γ Ω (cid:19) 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:16) Σ (cid:0) ( γ − ) ¯ Ω k + Ω k (cid:1) + Ω (cid:0) ( γ − ) ¯ Ω k + Ω k (cid:1) + (cid:0) ¯ Ω k − (cid:1)(cid:0) ( γ − ) ¯ Ω k + Ω k (cid:1) + Σ ¯ Ω k (cid:17) ∂ g ∂ Ω k + (cid:16) Ω (cid:0) ( γ − ) ¯ Σ + Σ (cid:1) + ¯ Ω k (cid:0) ( γ − ) ¯ Σ + Σ − (cid:1) + (cid:16) ¯ Σ − (cid:17)(cid:0) ( γ − ) ¯ Σ + Σ (cid:1)(cid:17) ∂ g ∂ Σ + (cid:16) Σ (cid:0) ( γ − ) ¯ Ω + Ω (cid:1) + ¯ Ω k (cid:0) ( γ − ) ¯ Ω + Ω (cid:1) + (cid:16) ¯ Ω − (cid:17)(cid:0) ( γ − ) ¯ Ω + Ω (cid:1)(cid:17) ∂ g ∂ Ω + Σ Ω ( − γ + cos ( ( Φ − t ω )) + ) g + (cid:16) − γ (cid:16) Σ + Ω + Ω k − (cid:17) + Σ + Ω cos ( ( Φ − t ω )) + Ω + Ω k − (cid:17) g + (cid:18) − γ Σ + Σ + (cid:19) g − Σ Ω 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:16) Σ (cid:0) ( γ − ) ¯ Ω k + Ω k (cid:1) + Ω (cid:0) ( γ − ) ¯ Ω k + Ω k (cid:1) + (cid:0) ¯ Ω k − (cid:1)(cid:0) ( γ − ) ¯ Ω k + Ω k (cid:1) + Σ ¯ Ω k (cid:17) ∂ g ∂ Ω k + (cid:32) Ω (cid:0) ( γ − ) ¯ Σ + Σ (cid:1) + ¯ Ω k (cid:0) ( γ − ) ¯ Σ + Σ − (cid:1) + (cid:16) ¯ Σ − (cid:17)(cid:0) ( γ − ) ¯ Σ + Σ (cid:1)(cid:33) ∂ g ∂ Σ + (cid:16) Σ (cid:0) ( γ − ) ¯ Ω + Ω (cid:1) + ¯ Ω k (cid:0) ( γ − ) ¯ Ω + Ω (cid:1) + (cid:16) ¯ Ω − (cid:17)(cid:0) ( γ − ) ¯ Ω + Ω (cid:1)(cid:17) ∂ g ∂ Ω + Ω Ω k ( − γ + cos ( ( Φ − t ω )) + ) g − Ω k ( ( γ − ) Σ + ) g + (cid:32) − γ (cid:16) (cid:16) Σ + Ω − (cid:17) + Ω k (cid:17) + Σ ( Σ − ) + Ω cos ( ( Φ − t ω )) + Ω + Ω k − (cid:33) g − Ω Ω k 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) b µ + ω cos ( Φ − t ω ) × (cid:32) γ (cid:16) − Σ − Ω − Ω k (cid:17) + Σ + (cid:16) Ω + (cid:17) cos ( ( Φ − t ω )) + Ω + Ω k (cid:33)(cid:33)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) < ∞ . Using Gronwall’s Lemma 10, 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 ) , (cid:12)(cid:12)(cid:12) ∆Ω k ( 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 . Finally, taking the limit as n → ∞ , we obtain H n →
0. Then, as H n →
0, the functions Ω , Σ , Ω k , Φ and¯ Ω , ¯ Σ , ¯ Ω k , ¯ Φ have the same limit as τ → ∞ . (cid:3) Appendix D: Numerical simulation
In this section we present the numerical evidencethat support the main theorem of section 3 by solv-ing numerically the full and time-averaged systems ob-tained for each metric, namely LRS Bianchi III andopen FLRW. For this purpose was elaborated an al-gorithm in the programming language
Python wherethe systems of differential equations were solved usingthe 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 the range of integra-tion − ≤ τ ≤
10 for the original systems and − ≤ τ ≤
100 for the averaged systems. All of them parti-tioned in 10000 data points. Furthermore, each full andtime-averaged systems were solved considering onlyone matter component, these are cosmological constant( γ = γ = γ = /
3) and stiff fluid ( γ = Ω = Ω m ≡ Ω m ≡
0. Finally we have considered the fixed con-stants µ = √ / b = √ / ω = √
2, that lead toa value of f = b µ ω − µ = /
10, that fulfills the condi-tion f ≥
0. With this values a generalized harmonicpotential of the form V ( φ ) = φ + ( − cos ( φ )) (D1)is obtained.
1. LRS Bianchi III
For the LRS Bianchi III metric we integrate:1. The full system given by (A2).2. The time-averaged system: ∂ τ ¯ Ω = ¯ Ω (cid:16) γ (cid:0) − ¯ Σ − ¯ Ω − ¯ Ω k (cid:1) + Σ + Ω + Ω k − (cid:17) , ∂ τ ¯ Σ = (cid:32) ¯ Σ (cid:16) γ (cid:0) − ¯ Σ − ¯ Ω − ¯ Ω k (cid:1) + Σ + Ω + Ω k − (cid:17) + Ω k (cid:33) , ∂ τ ¯ Ω k = ¯ Ω k (cid:16) γ (cid:0) − ¯ Σ − ¯ Ω − ¯ Ω k (cid:1) + Σ − Σ + Ω + Ω k − (cid:17) , ∂ τ ¯ Φ = , ∂ τ H = − H (cid:16) γ (cid:0) − ¯ Σ − ¯ Ω − ¯ Ω k (cid:1) + Σ + Ω + Ω k (cid:17) , ∂ τ t = / H , (D2)where as initial conditions we use the seven data setpresented in Table I for a better comparison of bothsystems. It is important to mention that the firstsix initial conditions correspond to the initial condi-tions presented in Table 2 of [134] and the additionaldata set vii is obtained considering current values of Ω m ( ) = .
315 and Ω k ( ) = .
001 according to [191].Even more, the model presented in [134] is containedin our model when γ = b = f = ω = µ = Ω m = Ω k = − Σ − Ω (then γ does not appearsin the model presented in [134]) with the identification Ω (cid:55)→ Ω . As can be seen in Figures 12(a) and 12(b)where some solutions of the full system (A2) and time-averaged system (D2) are presented; showing that our results are in complete agreement with the results pre-sented in [134] for the limiting case.In Figures 13(a)-20(b) are presented projections ofsome solutions of the full system (A2) and time-averaged system (D2) in the ( Σ , H , Ω ) and ( Ω k , H , Ω ) space with their respective projection when H = γ = γ = γ = ). Figures 19(a)-20(b)show solutions for a background fluid corresponding4 Table I: Seven initial data sets for the simulation of the full system (A2) and time-averaged system (D2). All initial conditionsare chosen in order to fulfill the equality Σ ( ) + Ω ( ) + Ω k ( ) + Ω m ( ) = Sol. H ( ) Σ ( ) Ω ( ) Ω k ( ) Ω m ( ) ϕ ( ) t ( ) i 0 . . . .
09 0 0 0ii 0 . . . .
74 0 0 0iii 0 . . . .
54 0 0 0iv 0 .
02 0 .
48 0 .
02 0 . . .
48 0 .
02 0 . . . .
01 0 .
74 0 0 0vii 0 . .
684 0 .
001 0 .
315 0 0
Table II: Seven initial data sets for the simulation of the full system (A5) and time-averaged system (D3) for the FLRW metricwith negative curvature ( k = − Ω ( ) + Ω k ( ) + Ω m ( ) = Sol. H ( ) Ω ( ) Ω k ( ) Ω m ( ) ϕ ( ) t ( ) i 0 . . .
09 0 .
01 0 0ii 0 . . .
74 0 .
16 0 0iii 0 . . .
54 0 .
36 0 0iv 0 .
02 0 .
02 0 . . . .
02 0 . . . .
01 0 . .
69 0 0vii 0 . .
684 0 .
001 0 .
315 0 0to a stiff fluid ( γ =
2. FLRW k = − For the FLRW metric with negative curvature ( k = − Ω k >
0) we integrate: 1. The full system given by (A5).2. The time-averaged system: ∂ τ ¯ Ω = − ¯ Ω (cid:0) ( γ − ) (cid:0) ¯ Ω − (cid:1) + ( γ − ) ¯ Ω k (cid:1) , ∂ τ ¯ Ω k = − ¯ Ω k (cid:0) ( γ − ) ¯ Ω − γ + ( γ − ) ¯ Ω k + (cid:1) , ∂ τ ¯ Φ = , ∂ τ H = − H (cid:2) γ ( − Ω k − Ω ) + Ω + Ω k (cid:3) , ∂ τ t = / H , (D3)where as initial conditions we use the seven data setpresented in the Table II for a better comparison ofboth system and where the data set vii is obtainedconsidering the current values of Ω m ( ) = .
315 and Ω k ( ) = .
001 according to [191].In figures 21(a)-24(b) are presented projections ofsome solutions of the full system (A5) and time- averaged system (D3) for the FLRW metric with neg-ative curvature ( k = −
1) in the ( Ω k , H , Ω ) space withtheir respective projection when H = γ = γ = γ = / γ = H tends asymptotically to constant values depending on the initial conditions which is consistent to de Sit-ter expansion. In addition, for any γ < and Ω k > Ω k →
0. On the other hand, when γ > and Ω k > Ω k → [1] A. H. Guth, Phys. Rev. D , 347 (1981) [Adv. Ser.Astrophys. Cosmol. , 139 (1987)][2] B. Ratra and P.J.E. Peebles, Phys. Rev D. ,103 (2007)[8] G. Leon, A. Paliathanasis and J. L. Morales-Mart´ınez,Eur. Phys. J. C , no.9, 753 (2018)[9] Y.F. Cai, E.N. Saridakis, M.R. Setare and J.-Q. Xia,Phys. Rep. , 1 (2010)[10] Z. K. Guo, Y. S. Piao, X. M. Zhang and Y. Z. Zhang,Phys. Lett. B , 177-182 (2005)[11] B. Feng, M. Li, Y. S. Piao and X. Zhang, Phys. Lett. B , 101-105 (2006)[12] H. Wei, R. G. Cai and D. F. Zeng, Class. Quant. Grav. , 3189-3202 (2005)[13] X. F. Zhang, H. Li, Y. S. Piao and X. M. Zhang, Mod.Phys. Lett. A , 231-242 (2006)[14] X. Zhang, Commun. Theor. Phys. , 762-768 (2005)[15] R. Lazkoz and G. Leon, Phys. Lett. B , 303-309(2006)[16] M. Setare and E. Saridakis, Phys. Lett. B , 177-181(2008)[17] M. Setare and E. Saridakis, Phys. Lett. B , 331-338(2009)[18] G. Leon, R. Cardenas and J. L. Morales,[arXiv:0812.0830 [gr-qc]].[19] G. Leon, Y. Leyva and J. Socorro, Phys. Lett. B ,285-297 (2014)[20] G. Le´on Torres, “Qualitative analysis and char-acterization of two cosmologies including scalarfields,” PhD thesis. Universidad Central de Las Villas[arXiv:1412.5665 [gr-qc]][21] S. Mishra and S. Chakraborty, Eur. Phys. J. C ,no.11, 917 (2018)[22] M. Marciu, Phys. Rev. D , no.4, 043508 (2019)[23] M. Marciu, Eur. Phys. J. C (2020) no.9, 894[24] N. Dimakis and A. Paliathanasis, [arXiv:2001.09687[gr-qc]][25] S.V. Chervon, Quantum Matter , 71 (2013)[26] I.V. Fomin, J. Phys.: Conf. Ser. 918, 012009 (2017)[27] A. Paliathanasis, Class. Quant. Grav. (2020) no.19,195014 [28] E. Elizalde, S. Nojiri and S. D. Odintsov, Phys. Rev. D (2004) 043539[29] E. Elizalde, S. Nojiri, S. D. Odintsov, D. Saez-Gomezand V. Faraoni, Phys. Rev. D (2008) 106005[30] A. Paliathanasis, G. Leon and S. Pan, Gen. Rel. Grav. , no. 9, 106 (2019)[31] Jordan, P. Z. Physik 157, 112–121 (1959)[32] C. Brans and R. H. Dicke, Phys. Rev. , 925 (1961)[33] G. W. Horndeski, Int. J. Theor. Phys. , 363 (1974)[34] E. J. Copeland, E. W. Kolb, A. R. Liddle and J. E. Lid-sey, Phys. Rev. D (1993) 252[35] J. E. Lidsey, A. R. Liddle, E. W. Kolb, E. J. Copeland,T. Barreiro and M. Abney, Rev. Mod. Phys. (1997)37[36] J. Ibanez, R. J. van den Hoogen and A. A. Coley, Phys.Rev. D (1995) 928[37] A. A. Coley, J. Ibanez and R. J. van den Hoogen, J.Math. Phys. (1997) 5256[38] E. J. Copeland, I. J. Grivell, E. W. Kolb and A. R. Lid-dle, Phys. Rev. D (1998) 043002[39] A. A. Coley and R. J. van den Hoogen, Phys. Rev. D (2000) 023517[40] A. P. Billyard and A. A. Coley Phys. Rev. D (2000)083503[41] A. Coley and M. Goliath, Class. Quant. Grav. (2000) 2557[42] A. Coley and M. Goliath, Phys. Rev. D (2000)043526[43] A. Coley and Y. J. He, Gen. Rel. Grav. (2003) 707[44] R. Curbelo, T. Gonzalez, G. Leon and I. Quiros, Class.Quant. Grav. (2006) 1585[45] T. Gonzalez, G. Leon and I. Quiros, astro-ph/0502383[46] S. Capozziello, S. Nojiri and S. D. Odintsov, Phys.Lett. B (2006) 597[47] T. Gonzalez, G. Leon and I. Quiros, Class. Quant.Grav. (2006) 3165[48] T. Gonzalez and I. Quiros, Class. Quant. Grav. (2008) 175019[49] O. Hrycyna and M. Szydlowski, Phys. Rev. D (2007) 123510[50] G. Leon and E. N. Saridakis, Phys. Lett. B (2010)1[51] G. Leon and E. N. Saridakis, JCAP (2009) 006[52] G. Leon, Y. Leyva, E. N. Saridakis, O. Martin andR. Cardenas, Falsifying field-based dark energy mod-els, in Dark Energy: Theories, Developments, andImplications (New York: Nova Science Publishers)arXiv:0912.0542 [gr-qc][53] G. Leon and E. N. Saridakis, Class. Quant. Grav. (2011) 065008 [54] J. Miritzis, J. Phys. Conf. Ser. (2011) 012024.[55] S. Basilakos, M. Tsamparlis and A. Paliathanasis,Phys. Rev. D (2011) 103512[56] C. Xu, E. N. Saridakis and G. Leon, JCAP (2012)005[57] M. Jamil, D. Momeni and R. Myrzakulov, Eur. Phys.J. C (2012) 207[58] G. Leon and E. N. Saridakis JCAP (2013) 025[59] G. Leon, J. Saavedra and E. N. Saridakis Class. Quant.Grav. (2013) 135001[60] M. A. Skugoreva, S. V. Sushkov and A. V. Toporensky,Phys. Rev. D (2013) 083539 Erratum: [Phys. Rev.D (2013) no.10, 109906][61] C. R. Fadragas, G. Leon and E. N. Saridakis, Class.Quant. Grav. (2014) 07501[62] O. Minazzoli and A. Hees, Phys. Rev. D (2014)023017[63] G. Kofinas, G. Leon and E. N. Saridakis, Class. Quant.Grav. (2014) 175011[64] A. Paliathanasis and M. Tsamparlis, Phys. Rev. D (2014) no.4, 043529[65] G. Leon and E. N. Saridakis, JCAP (2015) 031[66] A. Paliathanasis, M. Tsamparlis, S. Basilakos andJ. D. Barrow, Phys. Rev. D (2015) no.12, 123535[67] G. Leon and E. N. Saridakis, JCAP (2015) 009[68] T. Harko, F. S. N. Lobo, J. P. Mimoso and D. Pav´on,Eur. Phys. J. C (2015) 386[69] A. R. Solomon, doi:10.1007/978-3-319-46621-7arXiv:1508.06859 [gr-qc].[70] R. De Arcia, T. Gonzalez, G. Leon, U. Nucamendi andI. Quiros, Class. Quant. Grav. (2016) no.12, 125036[71] J. D. Barrow and A. Paliathanasis, Phys. Rev. D (2016) no.8, 083518[72] J. D. Barrow and A. Paliathanasis, Gen. Rel. Grav. (2018) no.7, 82[73] N. Dimakis, A. Giacomini, S. Jamal, G. Leon andA. Paliathanasis, Phys. Rev. D (2017) no.6, 064031[74] M. Cruz, A. Ganguly, R. Gannouji, G. Leon andE. N. Saridakis, Class. Quant. Grav. (2017) no.12,125014[75] J. Matsumoto and S. V. Sushkov, JCAP (2018)040[76] A. Giacomini, S. Jamal, G. Leon, A. Paliathanasis andJ. Saavedra, Phys. Rev. D (2017) no.12, 124060[77] B. Alhulaimi, R. J. Van Den Hoogen and A. A. Coley,JCAP (2017) 045[78] L. Karpathopoulos, S. Basilakos, G. Leon,A. Paliathanasis and M. Tsamparlis, Gen. Rel.Grav. (2018) no.7, 79[79] A. Paliathanasis, Mod. Phys. Lett. A (2017) no.37,1750206[80] R. De Arcia, T. Gonzalez, F. A. Horta-Rangel,G. Leon, U. Nucamendi and I. Quiros, Class. Quant.Grav. (2018) no.14, 145001[81] M. Tsamparlis and A. Paliathanasis, Symmetry (2018) no.7, 233[82] J. D. Barrow and A. Paliathanasis, Eur. Phys. J. C (2018) no.9, 767[83] R. J. Van Den Hoogen, A. A. Coley, B. Alhulaimi,S. Mohandas, E. Knighton and S. O’Neil, JCAP (2018) 017[84] G. Leon, A. Paliathanasis and L. Velazquez Abab, Gen. Rel. Grav. (2020) 71[85] F. Humieja and M. Szydłowski, Eur. Phys. J. C (2019) no.9, 794[86] I. Quiros, Int. J. Mod. Phys. D (2019) no.07,1930012[87] G. Leon and A. Paliathanasis, Eur. Phys. J. C (2019)no.9, 746[88] A. Paliathanasis and G. Leon, ZnA, 75, 523, 2020[89] S. Basilakos, G. Leon, G. Papagiannopoulos andE. N. Saridakis, Phys. Rev. D (2019) no.4, 043524[90] M. Shahalam, R. Myrzakulov and M. Y. Khlopov, Gen.Rel. Grav. (2019) no.9, 125[91] A. Paliathanasis, G. Papagiannopoulos, S. Basilakosand J. D. Barrow, Eur. Phys. J. C (2019) no.8, 723[92] G. Leon, A. Coley and A. Paliathanasis, Annals Phys. (2020) 168002[93] S. Nojiri, S. D. Odintsov and V. K. Oikonomou, An-nals Phys. (2020) 168186[94] S. Foster, Class. Quant. Grav. , 3485 (1998)[95] J. Miritzis, Class. Quant. Grav. , 2981 (2003)[96] R. Giambo, F. Giannoni and G. Magli, Gen. Rel. Grav. , 21 (2009).[97] G. Leon and C. R. Fadragas, Dynamical Systems: AndTheir Applications (Saarbr¨ucken: LAP Lambert Aca-demic Publishing), arXiv:1412.5701 [gr-qc][98] 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]][99] C. R. Fadragas and G. Leon Class. Quant. Grav. ,no. 19, 195011 (2014)[100] 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[101] G. Leon, Class. Quant. Grav. , 035008 (2009)[102] R. Giambo and J. Miritzis, Class. Quant. Grav. ,095003 (2010)[103] K. Tzanni and J. Miritzis, Phys. Rev. D , no. 10,103540 (2014) Addendum: [Phys. Rev. D , no. 12,129902 (2014)][104] R. J. van den Hoogen, A. A. Coley and D. Wands,Class. Quant. Grav. , 1843 (1999)[105] A. Albrecht and C. Skordis, Phys. Rev. Lett. , 2076(2000)[106] E. J. Copeland, A. R. Liddle and D. Wands, Phys. Rev.D , 4686 (1998)[107] R. Giamb`o, J. Miritzis and A. Pezzola, Eur. Phys. J.Plus , no.4, 367 (2020)[108] A. Cid, F. Izaurieta, G. Leon, P. Medina and D. Nar-bona, JCAP , 041 (2018)[109] A. Alho and C. Uggla, J. Math. Phys. , no. 1, 012502(2015)[110] E. A. Coddington y Levinson, N. Theory of Ordi-nary Differential Equations, New York, MacGraw-Hill, (1955)[111] J. K. Hale, Ordinary Differential Equations, New York,Wiley (1969)[112] D. K. Arrowsmith y C. M. Place, An introduction to dynamical systems, Cambridge University Press,Cambridge, England, (1990)[113] S. Wiggins. Introduction to Applied Nonlinear dynam-ical systems and Chaos. Springer (2003)[114] L. Perko, Differential equations and dynamical sys-tems, third edition (Springer-Verlag, New York, 2001).[115] V.I. Arnold, Ordinary differential equations. Cam-bridge: M.I.T. Press., 1973[116] M. W. Hirsch and S. Smale. Differential equations, dy-namical systems, and linear algebra. New York: Aca-demic Press (1974)[117] J. Hale. Ordinary differential equations. Malabar,Florida: Robert E. Krieger Publishing Co., Inc. (1980)[118] Lasalle, J. P., J. Diff. Eq., , pp. 57-65, 1968[119] B. Aulbach, Continuous and Discrete Dynamics nearManifolds of Equilibria (Lecture Notes in Mathemat-ics No. 1058, Springer, 1984)[120] R. Tavakol, Introduction to dynamical systems, ch4. Part one, pp. 84–98, Cambridge University Press,Cambridge, England, (1997)[121] 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[122] A. A. Coley, Introduction to Dynamical Systems. Lec-ture Notes for Math 4190/5190 (1994).[123] A. A. Coley, gr-qc/9910074.[124] Bassemah Alhulaimi (2017), Einstein-Aether Cosmo-logical Scalar Field Models (Phd Thesis, DalhousieUniversity).[125] V. G. LeBlanc, D. Kerr and J. Wainwright, Class.Quant. Grav. , 513 (1995).[126] J. M. Heinzle and C. Uggla, Class. Quant. Grav. ,015009 (2010).[127] A. D. Rendall, Class. Quant. Grav. , 667 (2007).[128] A. Alho, J. Hell and C. Uggla, Class. Quant. Grav. ,no. 14, 145005 (2015)[129] A. Paliathanasis, S. Pan and S. Pramanik, Class.Quant. Grav. , no. 24, 245006 (2015)[130] G. Leon and F. O. F. Silva, [arXiv:1912.09856 [gr-qc]][131] Genly Leon and Felipe Orlando Franz Silva 2021Class. Quantum Grav. , 012702(2012)[134] D. Fajman, G. Heißel and M. Maliborski, Class.Quant. Grav. , no.13, 135009 (2020)[135] David Fajman, Gernot Heißel, and Jin Woo JangarXiv: 2006.12844v2[136] G. Leon and F. O. F. Silva, Class. Quant. Grav. (2020) no.24, 245005[137] F. Dumortier and R. Roussarie (1995) “Canard cy-cles and center manifolds”, (Memoirs of the AmericanMathematical Society, 577).[138] N. Fenichel (1979) Geometric singular perturbationtheory for ordinary differential equations. Journal ofDifferential Equations , 53-98[139] G. Fusco and J.K. Hale, Journal of Dynamics and Dif-ferential Equations 1, 75 (1988)[140] N. Berglund and B. Gentz, Noise-Induced Phenomenain Slow-Fast Dynamical Systems, Series: Probabilityand Applications, Springer-Verlag: London, (2006)[141] M. H. Holmes (2013) “Introduction to Perturbation Methods”, (Springer Science+Business Media NewYork, ISBN 978-1-4614-5477-9)[142] 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[143] Ferdinand Verhulst, (2000) “Methods and Applica-tions of Singular Perturbations: Boundary Layers andMultiple Timescale Dynamics” (Springer-Verlag NewYork, ISBN 978-0-387-22966-9)[144] M. Sharma, M. Shahalam, Q. Wu and A. Wang, JCAP , 003 (2018)[145] L. McAllister, E. Silverstein, A. Westphal andT. Wrase, JHEP (2014), 123[146] G. Leon, S. Cu´ellar, E. Gonz´alez, S. Lepe, C. Micheaand A. D. Millano, [arXiv:2102.05495 [gr-qc]].[147] G. Leon, E. Gonz´alez, S. Lepe, C. Michea andA. D. Millano, [arXiv:2102.05551 [gr-qc]].[148] K.C. Jacobs, Astrophys J. 153, 661 (1968)[149] C.B Collins and S.W. Hawking, Astroph. J. 180, 317(1973)[150] J.D. Barrow, Mon. Not. R. astron. Soc. 175, 359 (1976)[151] J.D. Barrow and D.H. Sonoda, Phys. Reports, 139, 1(1986)[152] M. Thorsrud, B.D. Normann and T.S. Pereira, Class.Quantum Grav. , 065015 (2020)[153] S. Byland and D. Scialom, Phys. Rev. D , 6065-6074(1998)[154] U. Nilsson and C. Uggla, Class. Quant. Grav. 13, 1601(1996)[155] C. Uggla and H. Zur-Muhlen, Class. Quant. Grav. 71365 (1990)[156] M. Goliath, U. S. Nilsson and C. Uggla, Class. Quant.Grav. 15, 167 (1998)[157] B. J. Carr, A. A. Coley, M. Goliath, U. S. Nilsson andC. Uggla, Class. Quant. Grav. 18, 303 (2001)[158] A. Coley and M. Goliath, Class. Quant. Grav. 17, 2557(2000)[159] A. Coley and M. Goliath, Phys. Rev. D 62, 043526(2000)[160] M. Heusler, Phys. Lett. B 253, 33 (1991)[161] J.M. Aguirregabiria, A. Feinstein and J. Ibanez, Phys.Rev. D 48, 4662 (1993)[162] T. Christodoulakis, Th. Grammenos, Ch. Helias andP.G. Kevrekidis, J. Math. Phys. 47, 042505 (2006)[163] M. Tsamparlis and A. Paliathanasis, Gen. Relat.Gravit. 43, 1861 (2011)[164] A.A. Coley, J. Ibanez and R.J. van den Hoogen, J.Math. Phys. 38, 5256 (1997)[165] J. Ibanez, R.J. van den Hoogen and A.A. Coley, Phys.Rev. D 51, 928 (1995)[166] V.A. Belinskii, E.M. Lifhitz and I.M. Khalatnikov,JETP 33, 1061 (1971)[167] K. Adhav, A. Nimkar, R. Holey, Int. J. Theor. Phys. L1 (2006) [172] T. Clifton and J.D. Barrow, Class Quantum Grav. , no.39, 1950326 (2019)[176] A. Paliathanasis, L. Karpathopoulos, A. Wojnar andS. Capozziello, Eur. Phys. J. C , no.4, 225 (2016)[177] A. A. Coley, W. C. Lim and G. Leon,[arXiv:0803.0905 [gr-qc]].[178] A. D. Linde, Phys. Lett. B (1983), 177-181[179] A. D. Linde, Phys. Lett. B (1986), 395-400[180] A. D. Linde, [arXiv:hep-th/0205259 [hep-th]][181] A. H. Guth, J. Phys. A (2007), 6811-6826[182] P. A. R. Ade et al. [Planck], Astron. Astrophys. (2014), A22[183] P. A. R. Ade et al. [BICEP2 and Planck], Phys. Rev. Lett. (2015), 101301[184] G. D’Amico, T. Hamill and N. Kaloper, Phys. Rev. D (2016) no.10, 103526[185] A. B. Balakin and A. F. Shakirzyanov, Universe (2020) no.11, 192[186] J. Wainwright and G. F. R. Ellis, Eds. Dynamical Sys-tems in Cosmology. Cambridge Univ. Press, Cam-bridge, 1997[187] Edward Arthur Milne, “Relativity, Gravitation andWorld Structure”, Oxford University Press, 1935.[188] S. M. Carroll, “Spacetime and Geometry,” San Fran-cisco, USA: Addison-Wesley (2004) 513 p[189] V. Mukhanov, “Physical Foundations of Cosmology,”UK: Cambridge University Press (2005), 27p[190] C. W. Misner, K. S. Thorne and J. A. Wheeler, “Grav-itation,” San Francisco 1973, 1279p[191] N. Aghanim et al. [Planck Collaboration] Astron. As-trophys. , A6 (2020). (a) Projections in the space ( Σ , H , Ω ) . The surface is given by the constraint Ω = − Σ . iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii 0.050 0.025 0.000 0.025 0.050 0.075 0.100 0.125 0.1500.50.60.70.80.91.0 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.14 ii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v vi 0.475 0.480 0.485 0.490 0.495 0.500 0.5050.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 12: Some solutions of the full system (A2) (blue) and time-averaged system (D2) (orange), corresponding to LRSBianchi III metric, when γ = b = f = ω = µ = Ω m = Ω k = − Σ − Ω , with the identification Ω (cid:55)→ Ω , forwhich the results of [134] are recovered. We have used 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 Ω = − Σ . iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.140.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.14 ii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v vi 0.475 0.480 0.485 0.490 0.495 0.500 0.5050.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 13: Some solutions of the full system (A2) (blue) and time-averaged system (D2) (orange) for the LRS Bianchi IIImetric when γ =
0, in the projection Ω k =
0. We have used for both systems the initial data sets presented in the Table I. k H i iiiii ivv-vivii (a) Projections in the space ( Ω k , H , Ω ) . The surface is given by the constraint Ω = − Ω k . k i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 k ivii ivii ivii ivii ivii ivii ivii0.5 0.6 0.7 0.8 0.9 1.0 k iiiii iv-v viiiiii iv-v viiiiii iv-v viiiiii iv-v viiiiii iv-v viiiiii iv-v viiiiii iv-v vi 0.735 0.740 0.745 0.750 0.755 0.760 k iv vvi iv vvi iv vvi iv vvi iv vvi iv vvi iv vvi (b) Projection in the space ( Ω k , Ω ) . The black line represent the constraint Ω = − Ω k . Figure 14: Some solutions of the full system (A2) (blue) and time-averaged system (D2) (orange) for the LRS Bianchi IIImetric when γ =
0, in the projection Σ =
0. 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 Ω = − Σ . iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.140.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.14 ii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v vi 0.475 0.480 0.485 0.490 0.495 0.500 0.5050.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 15: Some solutions of the full system (A2) (blue) and time-averaged system (D2) (orange) for the LRS Bianchi IIImetric when γ =
1, in the projection Ω k =
0. We have used for both systems the initial data sets presented in the Table I. k H i iiiii ivv-vivii (a) Projections in the space ( Ω k , H , Ω ) . The surface is given by the constraint Ω = − Ω k . k i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 k ivii ivii ivii ivii ivii ivii ivii0.5 0.6 0.7 0.8 0.9 1.0 k iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi 0.735 0.740 0.745 0.750 0.755 0.760 k iv vvi iv vvi iv vvi iv vvi iv vvi iv vvi iv vvi (b) Projection in the space ( Ω k , Ω ) . The black line represent the constraint Ω = − Ω k . Figure 16: Some solutions of the full system (A2) (blue) and time-averaged system (D2) (orange) for the LRS Bianchi IIImetric when γ =
1, in the projection Σ =
0. 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 Ω = − Σ . iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.140.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.14 ii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v vi 0.475 0.480 0.485 0.490 0.495 0.500 0.5050.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 17: Some solutions of the full system (A2) (blue) and time-averaged system (D2) (orange) for the LRS Bianchi IIImetric when γ = /
3, in the projection Ω k =
0. We have used for both systems the initial data sets presented in the Table I. k H i iiiii ivv-vivii (a) Projections in the space ( Ω k , H , Ω ) . The surface is given by the constraint Ω = − Ω k . k i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 k ivii ivii ivii ivii ivii ivii ivii0.5 0.6 0.7 0.8 0.9 1.0 k iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi 0.735 0.740 0.745 0.750 0.755 0.760 k iv vvi iv vvi iv vvi iv vvi iv vvi iv vvi iv vvi (b) Projection in the space ( Ω k , Ω ) . The black line represent the constraint Ω = − Ω k . Figure 18: Some solutions of the full system (A2) (blue) and time-averaged system (D2) (orange) for the LRS Bianchi IIImetric when γ = /
3, in the projection Σ =
0. 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 Ω = − Σ . iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii iii iiiiv-v vivii 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.140.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.14 ii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v viii iiiiv-v vi 0.475 0.480 0.485 0.490 0.495 0.500 0.5050.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 19: Some solutions of the full system (A2) (blue) and time-averaged system (D2) (orange) for the LRS Bianchi IIImetric when γ =
2, in the projection Ω k =
0. We have used for both systems the initial data sets presented in the Table I. k H i iiiii ivv-vivii (a) Projections in the space ( Ω k , H , Ω ) . The surface is given by the constraint Ω = − Ω k . k i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 k ivii ivii ivii ivii ivii ivii ivii0.5 0.6 0.7 0.8 0.9 1.0 k iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi 0.735 0.740 0.745 0.750 0.755 0.760 k iv vvi iv vvi iv vvi iv vvi iv vvi iv vvi iv vvi (b) Projection in the space ( Ω k , Ω ) . The black line represent the constraint Ω = − Ω k . Figure 20: Some solutions of the full system (A2) (blue) and time-averaged system (D2) (orange) for the LRS Bianchi IIImetric when γ =
2, in the projection Σ =
0. We have used for both systems the initial data sets presented in the Table I. k H i iiiii ivvvivii (a) Projections in the space ( Ω k , H , Ω ) . The surface is given by the constraint Ω = − Ω k . k i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii 0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 k ivii ivii ivii ivii ivii ivii ivii0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 k iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi 0.70 0.75 0.80 0.85 0.90 0.95 1.00 k iv viv viv viv viv viv viv v (b) Projection in the space ( Ω k , Ω ) . The black line represent the constraint Ω = − Ω k . Figure 21: Some solutions of the full system (A5) (blue) and time-averaged system (D3) (orange) for the FLRW metric withnegative curvature ( k = −
1) when γ ==
1) when γ ==
0. We have used for both systems the initial data sets presented in the Table II. k H i iiiiiiv vvivii (a) Projections in the space ( Ω k , H , Ω ) . The surface is given by the constraint Ω = − Ω k . k i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 k ivii ivii ivii ivii ivii ivii ivii0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 k iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi 0.70 0.72 0.74 0.76 0.78 0.80 0.82 0.84 k iv viv viv viv viv viv viv v (b) Projection in the space ( Ω k , Ω ) . The black line represent the constraint Ω = − Ω k . Figure 22: Some solutions of the full system (A5) (blue) and time-averaged system (D3) (orange) for the FLRW metric withnegative curvature ( k = −
1) when γ ==
1) when γ ==
1. We have used for both systems the initial data sets presented in the Table II. k H i iiiiiiv vvivii (a) Projections in the space ( Ω k , H , Ω ) . The surface is given by the constraint Ω = − Ω k . k i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 k ivii ivii ivii ivii ivii ivii ivii0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 k iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi 0.70 0.72 0.74 0.76 0.78 0.80 0.82 0.84 k iv viv viv viv viv viv viv v (b) Projection in the space ( Ω k , Ω ) . The black line represent the constraint Ω = − Ω k . Figure 23: Some solutions of the full system (A5) (blue) and time-averaged system (D3) (orange) for the FLRW metric withnegative curvature ( k = −
1) when γ ==
1) when γ == /
3. We have used for both systems the initial data sets presented in the Table II. k H i iiiiiiv vvivii (a) Projections in the space ( Ω k , H , Ω ) . The surface is given by the constraint Ω = − Ω k . k i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii i iiiii iv-vvivii 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 k ivii ivii ivii ivii ivii ivii ivii0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 k iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi iiiii iv-vvi 0.70 0.72 0.74 0.76 0.78 0.80 0.82 0.84 0.86 k iv viv viv viv viv viv viv v (b) Projection in the space ( Ω k , Ω ) . The black line represent the constraint Ω = − Ω k . Figure 24: Some solutions of the full system (A5) (blue) and time-averaged system (D3) (orange) for the FLRW metric withnegative curvature ( k = −
1) when γ ==