Averaging generalized scalar field cosmologies III: Kantowski--Sachs and closed Friedmann--Lema\^itre--Robertson--Walker models

Scalar field cosmologies with a generalized harmonic potential and matter with energy density $\rho_m$, pressure $p_m$, and barotropic equation of state (EoS) $p_m=(\gamma-1)\rho_m, \; \gamma\in[0,2]$ in Kantowski-Sachs (KS) and closed Friedmann--Lema\^itre--Robertson--Walker (FLRW) metrics are investigated. We use methods from non--linear dynamical systems theory and averaging theory considering a time--dependent perturbation function $D$. We define a regular dynamical system over a compact phase space, obtaining global results. That is, for KS metric the global late--time attractors of full and time--averaged systems are two anisotropic contracting solutions, which are non--flat locally rotationally symmetric (LRS) Kasner and Taub (flat LRS Kasner) for $0\leq \gamma \leq 2$, and flat FLRW matter--dominated universe if $0\leq \gamma \leq \frac{2}{3}$. For closed FLRW metric late--time attractors of full and averaged systems are a flat matter--dominated FLRW universe for $0\leq \gamma \leq \frac{2}{3}$ as in KS and Einstein-de Sitter solution for $0\leq\gamma<1$. Therefore, time--averaged system determines future asymptotics of full system. Also, oscillations entering the system through Klein-Gordon (KG) equation can be controlled and smoothed out when $D$ goes monotonically to zero, and incidentally for the whole $D$-range for KS and for closed FLRW (if $0\leq \gamma<1$) too. However, for $\gamma\geq 1$ closed FLRW solutions of the full system depart from the solutions of the averaged system as $D$ is large. Our results are supported by numerical simulations.


Introduction
Scalar fields have played important roles in the physical description of the universe , particularly, in the inflationary scenario. For example, chaotic inflation is a model of cosmic inflation in which the potential term takes the form of the harmonic potential V (φ ) = m 2 φ φ 2 2 [18][19][20][21]. Scalar field models can be examined by means of qualitative techniques of dynamical systems [22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38], which allow a stability analysis of the solutions. Complementary, asymptotic methods and averaging theory [39][40][41][42][43][44][45] are helpful to obtain relevant information about the solutions space of scalar field cosmologies: (i) in the vacuum and (ii) in the presence of matter [46,47]. In this process one idea is to construct a time-averaged version of the original system. By solving this version the oscillations of the original system are smoothed out [48]. This can be achieved for Bianchi I, flat FLRW, Bianchi III, and negatively curved FLRW metrics where the Hubble parameter H plays the role of a time-dependent perturbation function which controls the magnitude of the error between the solutions of full and time-averaged problems as H → 0 [46,47]. The conformal algebra of Bianchi III and Bianchi V spacetimes, which admits a proper conformal Killing vector, was studied in [49]. In [50] the method of Lie symmetries was applied for Wheeler-De Witt equation in Bianchi Class A cosmologies for minimally coupled scalar field gravity and Hybrid Gravity in General Relativity (GR). Several invariant solutions were deter-mined and classified according to the form of the scalar field potential by means of these symmetries. Based on references [51][52][53][54] we started the "Averaging generalized scalar field cosmologies" program [46]. The idea is to use asymptotic methods and averaging theory to obtain relevant information about the solutions space of scalar field cosmologies in the presence of matter with energy density ρ m and pressure p m with a barotropic EoS p m = (γ − 1)ρ m (with barotropic index γ ∈ [0, 2]) minimally coupled to a scalar field with generalized harmonic potential (13). This research program has three steps according to three cases of study: (I) Bianchi III and open FLRW model [46], (II) Bianchi I and flat FLRW model [47], and (III) KS and closed FLRW. In reference [46] Bianchi III metrics were studied and written conveniently as where g H 2 = dϑ 2 + sinh 2 (ϑ )dζ 2 denotes the 2-metric of negative constant curvature on hyperbolic 2-space and the lapse function N was set one. Moreover, functions A(t) and B(t) are interpreted as scale factors.
There we calculate the characteristic scale factor to obtain The formal definition of is given in eq. (43). For FLRW, it is the scale factor of the universe. By convention we set t = 0 as the present time.
In reference [46] late-time attractors of original and time-averaged systems for LRS Bianchi III were found to be the following: 1. A solution with the asymptotic metric ds 2 = −dt 2 + 2 0 3γH 0 t 2 + 1 4 3γ for 0 ≤ γ ≤ 2 3 , where 0 is the current value of the characteristic scale factor and H 0 is the current value of the Hubble factor. It represents a matterdominated flat FLRW universe.

A solution with the asymptotic metric
for 2 3 < γ < 1, which represents a matter-curvature scaling solution. 3. A solution with asymptotic metric for 1 ≤ γ ≤ 2, where c 1 is a constant. It corresponds to the Bianchi III form of flat spacetime ( [55] p 193, eq. (9.7)).
In reference [46] the open FLRW model was studied, whose metric is given by where dΩ 2 = dϑ 2 + sin 2 ϑ dζ 2 and A(t) is the scale factor of the isotropic and homogeneous universe. The late-time attractors are the following: 1. A solution with asymptotic metric for 0 ≤ γ ≤ 2 3 , corresponding to a flat matter-dominated FLRW universe.

A solution with asymptotic metric
for 2 3 < γ < 2, corresponding to a curvature dominated Milne solution ( [59]; [55] Sect. 9.1.6, eq. (9.8), [60][61][62]). In all metrics the matter-dominated flat FLRW universe represents quintessence fluid if 0 < γ < 2 3 . The chosen barotropic equation of state can mimic one of several fluids of interest in early-time and late-time cosmology. Typical values are γ = 2 corresponding to stiff matter, γ = 4 3 corresponding to radiation, γ = 1 corresponding to cold dark matter (CDM), γ = 2 3 corresponding to Dirac-Milne universe, and γ = 0 corresponding to cosmological constant (CC). According to our stability analysis, the ranges 0 ≤ γ < 2, 0 ≤ γ ≤ 2 3 , 0 < γ < 1, and 1 < γ < 2 are found. Special cases γ = 1 and γ = 2, corresponding to bifurcations parameters where the stability changes, are treated separately. It is important to mention that stiff matter is a component present in a very early evolution, and had a role before the last scattering epoch. The last scattering epoch is an important cornerstone in the cosmological history since after that, the cosmic microwave background (CMB) photons freely travelled through the universe providing a photographic picture of the universe at that epoch. Also, radiation is an early relevant cosmic component, although even today we have (tiny) traces of it. Moreover, CDM is an important component in current cosmology. The range 0 < γ < 2 3 corresponds to a quintessence and γ = 0 represents the CC ("omnipresent" in cosmic evolution). In the current state of cosmology the dark components satisfy 0 ≤ γ < 2 3 (dark energy) and γ = 1 (CDM). The Dirac-Milne universe is characterized by γ = 2 3 .
For FLRW metric, the characteristic length scale coincides with the scale factor of the universe. Thus, Friedmann's usual scheme leads to where z is the redshift and we use k to indicate if the model is a closed (k = 1), flat (k = 0), or open (k = −1). If γ = 2 3 and using 1+z = 0 / , we have from the above equations Hence, if k = 0, For vacuum ρ m = 0 and k = −1 (open case) Thus, we obtain a Dirac-Milne universe. The behavior H (z) ≈ (1 + z) is also satisfied in presence of matter (ρ m = 0) and for k = 0. More generically, for k = −1 and a fluid that dilutes over time for which γ > 2 3 , then, (1 + z) .
For γ > 2 3 , the dominant term as z → −1 is given by (9). Namely, the asymptotic evolution is towards a Dirac-Milne type evolution. On the contrary, for γ < 2 3 the universe becomes matter-dominated. Following our research program, in reference [47] case (II) was studied. Late-time attractors of the original and time-averaged systems for Bianchi I and flat FLRW are the following.
1. For 0 ≤ γ < 1 the late-time attractor is a matterdominated FLRW universe mimicking de Sitter, quintessence, or zero acceleration solutions with asymptotic metric 2. For 1 < γ ≤ 2 the late-time attractor is an equilibrium solution with asymptotic metric where c 1 and c 2 are constants. This solution can be associated with Einstein-de Sitter solution ( [55], Sec 9.1.1 (1)) with γ = 1).
This paper, which is the third of the series, is devoted to the case (III) KS and positively curved FLRW metrics. We will prove that the quantity where 3 R is the 3-Ricci curvature of spatial surfaces (if the congruence u is irrotational), plays the role of a time-dependent perturbation function which controls the magnitude of the error between the solutions of full and time-averaged problems. The analysis of the system is therefore reduced to study the corresponding time-averaged equations as the time-dependent perturbation function D goes monotonically to zero for a finite time interval. The region where the perturbation parameter D changes its monotony from monotonic decreasing to monotonic increasing is analyzed by a discrete symmetry and by defining the variable T = D/(1 + D) that maps [0, ∞) to a finite interval [0, 1). Consequently, the limit D → +∞ corresponds to T = 1 and the limit D → 0 corresponds to T = 0.
The paper is organized as follows. In section 2 we motivate our choice of potential and the topic of averaging in the context of differential equations. In section 3 we introduce the model under study. In section 4 we apply averaging methods to analyze periodic solutions of a scalar field with self-interacting potentials within the class of generalized harmonic potentials [54]. In section 4.1 KS model is studied by using D-normalization, rather than Hubble-normalization, because Hubble factor is not monotonic for closed universes. FLRW models with k = +1 (positive curvature) are studied in section 4.2. In section 5 we study the resulting timeaveraged systems for KS and positively curved FLRW models. In particular, in section 5.1 KS model is studied. The FLRW model with k = +1 is studied in section 5.2. In section 6 a regular dynamical system defined on a compact phase space is constructed. This allows to find global results for KS and closed FLRW models. Finally, in section 7 our main results are discussed. In Appendix A the proof of the main theorem is given. In Appendix B numerical evidence supporting the results of section 4 is presented.
Potentials (13) or (14) are related but not equal to the monodromy potential of [63] used in the context of loop-quantum gravity, which is a particular case of the general monodromy potential [64]. The potential studied in references [63,64] for p = 2, i.e., V (φ ) = is not good to describe the late-time FLRW universe driven by a scalar field as shown by references [51,52,54] because it has two symmetric local negative minima which are related to Anti-de Sitter solutions. Setting µ = √ 2 2 and bµ = 2 in eq. (13) we recover the potential that was studied by [51,54].
Potentials (15) and (16) provide non-negative local minima which can be related to a late-time accelerated universe. Generalized harmonic potentials (14), (15) and (16) belong to the class of potentials studied in [65].
µ is a parameter, is relevant for axion models [66]. In [67] axionic dark matter with modified periodic poten- where Φ * is a parameter describing the basic state of the axion field, has been studied in the framework of the axionic extension of Einstein-aether theory. This periodic potential has minima at φ = nΦ * , where n ∈ Z, whereas maxima are found when n → m + 1 2 . Near the minimum where m A is the axion rest mass. In reference [68] an axion model given by two canonical scalar fields φ 1 and φ 2 coupled via the potential was investigated by combining standard dynamical systems tools and averaging techniques to investigate oscillations in Einstein-KG equations. As in references [46,47] methods from the theory of averaging nonlinear dynamical systems allow to prove that timedependent systems and their corresponding time-averaged versions have the same future asymptotic behavior. Thus, oscillations entering a nonlinear system through KG equation can be controlled and smoothed out as the Hubble factor H tends monotonically to zero.

Simple example of averaging problem
One approximation scheme which can be used to approximately solve the ordinary differential equationẋ = f(x,t, ε) with ε ≥ 0 and f periodic in t is to solve the unperturbed problemẋ = f(x,t, 0) by setting ε = 0 at first and then, with the use of the approximated unperturbed solution, to formulate variational equations in standard form which can be averaged. For example, consider the usual equation of a damped harmonic oscillator with given φ (0) andφ (0), where ω 2 is the undamped angular frequency of the oscillator, and the parameter ε = ζ ω, with ζ the damping ratio, is considered as a small parameter. The unperturbed problem admits the solutioṅ where r 0 and Φ 0 are constants depending on initial conditions. Using the variation of constants we propose the solution for the perturbed problem (18) aṡ (21) such that This procedure is called the amplitude-phase transformation in chapter 11 of [45]. Then, eq. (18) becomeṡ From eq. (23) it follows that r and Φ are slowly varying functions of time, and the system takes the forṁ y = ε f (y). The idea is to consider only nonzero average of the right-hand-side keeping r and Φ fixed and leaving out the terms with average zero and ignoring the slow-varying dependence of r and Φ on t through the averaging process Replacing r and Φ by their averaged approximationsr andΦ we obtain the systeṁr Solving eq. (25) withr(0) = r 0 andΦ(0) = Φ 0 , we obtainφ = r 0 e −εt sin(ωt − Φ 0 ) which is an accurate approximation of the exact solution as ε → 0 + .

General class of systems with a time-dependent perturbation function
Let us consider for example the Einstein-KG system The similarity between (18) and (26) suggests to treat the latter as a perturbed harmonic oscillator as well and to apply averaging in an analogous way. However, in contrast to ε, H is time-dependent and it is governed by evolution equation (27). Then, a surprising feature of such approach is the possibility of exploiting the fact that H is strictly decreasing and goes to zero by promoting Hubble parameter H to a time-dependent perturbation function in (26) controlling the magnitude of the error between solutions of the full and timeaveraged problems. Hence, with strictly decreasing H the error should decrease as well. Therefore, it is possible to obtain information about the large-time behavior of the more complicated full system via an analysis of the simpler averaged system of equations by means of dynamical systems techniques. This result is based on the monotony of H and its sign invariance. With this in mind, in [69] the long-term behavior of solutions of a general class of spatially homogeneous cosmologies, when H is positive strictly decreasing in t and lim t→∞ H(t) = 0, was studied. However, this analysis is not valid when the Hubble parameter is not a monotonic function as in the case of this study.

Spatially homogeneous and anisotropic spacetimes
The spatially homogeneous but anisotropic spacetimes are known as either Bianchi or KS cosmologies. In Bianchi models the spacetime manifold is foliated along the time axis with three dimensional homogeneous hypersurfaces. On the other hand, the isometry group of KS spacetime is R × SO(3) and it does not act simply transitively on spacetime, nor does it possess a subgroup with simple transitive action. Hence, this model is spatially homogeneous but it does not belong to the Bianchi classification. KS model approaches a closed FLRW model [70][71][72][73] when it isotropizes. In GR the Hubble parameter H is always monotonic for Bianchi I and Bianchi III. For Bianchi I the anisotropy decays on time for H > 0. Therefore, isotropization occurs [74]. Generically, in KS as well as for closed FLRW, the Hubble parameter is non monotonic and anisotropies would increase rather than vanish as H changes the sign. We refer the reader to [55,58, and references therein for applications of KS models, spatially homogeneous, and LRS metrics. The typical behavior of KS metric for perfect fluids, Vlasov matter, etc., is that the generic solutions are past and future asymptotic to the non-flat LRS Kasner vacuum solution, which have a big-bang (or big-crunch). Moreover, there exist non-generic solutions which are past (future) asymptotic to the anisotropic Bianchi I matter solution and others to the flat Friedman matter solution. The qualitative properties of positive-curvature models and the KS models with a barotropic fluid and a non-interacting scalar field with exponential potential V (φ ) = V 0 e λ φ , being λ a constant, were examined, e.g., in [33]. The main results are the following. For positively curved FLRW models and for λ 2 > 2 all the solutions start from and recollapse to a singularity, and they are not generically inflationary. For λ 2 < 2 the universe can either recollapse or expand forever. The KS model exhibits similar global properties of the closed FLRW models. In particular, for λ 2 > 2, all initially expanding solutions reach a maximum expansion and after that recollapse. These solutions are not inflationary nor does they isotropize. For λ 2 < 2 the models generically recollapse or expand forever to power-law inflationary flat FLRW solution. Intermediate behavior of KS as compared with closed FLRW is rather different.

General relativistic 1 + 3 orthonormal frame formalism
In this section we follow references [55,56] where the 1 + 3 orthonormal frame formalism was presented. A cosmological model (M , g, u) (representing the universe at a particular averaged scale) is defined by specifying the spacetime geometry through a Lorentzian metric g defined on the manifold M and a family of fundamental observers whose congruence of worldlines is represented by the four velocity field u, which will usually be the matter four velocity.
, and δ a b is the Kronecker's delta which is equal to 1 if a = b or equal to zero if a = b. A system of units in which 8πG = c = = 1 is used. The following symbols are used. R is the scalar curvature of the spacetime, g µν are the metric components, g is the determinant of the metric, φ is the scalar field, and V (φ ) is the scalar field potential defined by (13). A semicolon ";" as well as ∇ a indicates covariant derivatives. It is common to describe cosmological models in terms of a basis of vector fields {e a } and a dual basis of 1forms {ω a }, a = 0, 1, 2, 3. Any vector X can be written as X = X a e a in terms of this basis. The components of the metric tensor g relative to this basis are given by The line element can be symbolically written as In any coordinate chart there is a natural basis, namely {e a } can be chosen to be a coordinate basis {∂ /∂ x i } with the dual basis being the coordinate 1-forms The general basis vector fields e a and 1-forms ω a can be written in terms of a coordinate basis as follows Thus, any vector field X can be interpreted as a differential operator which acts on scalars as In particular, In terms of a coordinate basis the components of the metric tensor g relative to this basis are given by The line element can be symbolically written as Another special type of basis is the orthonormal frame, in which the four vector fields e a are mutually orthogonal and of unit length with e 0 timelike. The vector fields, thus, satisfy g(e a , e b ) = η ab , where η is the Minkowski metric η ab = diag(−1, 1, 1, 1). Given any basis of vector fields e a , the commutators [e a , e b ] are vector fields and hence they can be written as a linear combination of the basis vectors [e a , e b ] = γ c ab e c . The coefficients γ c ab (x i ) are called the commutation functions (which are 24 functions). For a coordinate basis {∂ /∂ x i } the commutation functions are all zero. If we use an orthonormal frame, the 24 commutation functions are the basic variables and the gauge freedom is an arbitrary Lorentz transformation. Applying the Jacobi identity for vector fields to γ c ab we obtain a set of 16 identities (Eqs. 1.18 in [55]) where we denote by e a ( f ) the action of the vector field e a on a scalar f as the differential operator e a f = e i a ∂ f ∂ x i using coordinate basis {∂ /∂ x i }. The identities (35), in conjunction with Einstein field equations rewritten as 1 and conservation equations give first-order evolution equations for some of the commutation functions and for some stress-energy tensor components and also provide a set of constraints involving only spatial derivatives, which is referred as the orthonormal frame formalism (Sect. It is well-known that a unit timelike vector field u(x i ), u a u a = −1, determines a projection tensor which at each point projects into the 3-space orthogonal to u. It follows that Therefore, we can define two derivatives, one along the vector u a defined byṪ a..b c..d = u e ∇ e T a..b c..d , and a projected derivative defined as D e T a..b reduces to the usual time derivative of a function. The covariant derivative of a unit timelike vector field u(x i ), u a u a = −1 can be decomposed in its irreducible parts as follows (see [55] page 18, [57] page 70) where σ ab is symmetric and trace-free, ω ab is antisymmetric and σ ab u b = 0 = ω ab u b . Physically, a timelike 1 Where Λ is the cosmological constant. See eq. (23) of reference [56]. unit time vector is generally chosen as the four velocity of the fluid and the quantitiesu a , θ , ω ab , ω a , σ ab are called the acceleration vector, rate of expansion scalar, vorticity tensor, vorticity vector, and the rate of shear tensor, respectively. The magnitude of the shear tensor is It follows that θ := u a ;a , u a := u a;b u b ,u a u a = 0, where η abcd is the totally antisymmetric permutation tensor such that η 0123 = √ −g, η 0123 = −1/ √ −g, where g denotes the determinant of the spacetime metric tensor g with Lorentzian signature. In Cosmology it is useful to define a representative length along the worldlines of u describing the volume expansion (contraction) behavior of the congruence completely by the equatioṅ where H is the Hubble parameter defined by In the orthonormal frame formalism, the components of the connection are simplified to The spatial frame vectors are denoted by {e α }, where the indices chosen from the first half of the Greek alphabet run from 1 to 3 and ε αβ γ denotes the alternating symbol (ε 123 = +1). Notice that greek indices are raised and lowered with the spatial metric tensor g αβ = δ αβ . The commutation functions are decomposed into algebraically simple quantities, some of which have a direct physical or geometrical interpretation. Firstly, the commutation function γ c ab with one zero index can be expressed as functions of geometrical quantities (42) of the timelike congruence u = e 0 and the quantity Ω α = 1 2 ε α µν e i µ e νi; j u j , which is the local angular velocity of the spatial frame {e α } with respect to a Fermi-propagated spatial frame.
Eq. (1.61) of [55] gives where σ αβ , ω α , andu α are spatial components of σ ab , ω a andu a . Secondly, spatial components γ µ αβ are decomposed into a 2-index symmetric object n αβ and a 1-index object a α as follows In order to incorporate the variety of matter sources in Einstein field equations the standard decomposition of the stress-energy tensor T ab with respect to the timelike vector u is used, where µ is the total energy density, q a is the energy current density, p is the isotropic pressure and π ab is the anisotropic pressure tensor. It follows that q a u a = 0, π ab u b = 0, π a a = 0, π ab = π ba .
If the congruence u is irrotational (ω ab = 0), the curvature of the 3-spaces orthogonal to the congruence can be expressed as where Θ ac is the rate of expansion tensor given by and R abcd is the Riemann tensor defined by The trace-free spatial Ricci tensor is defined by Also, the Weyl conformal curvature tensor can be expressed as It is useful to define the electric part E ab and magnetic part H ab of the Weyl conformal curvature tensor relative to u according to These tensors are symmetric, trace-free, and satisfy E ab u b = 0 = H ab u b .
When performing the 1 + 3 decomposition for an irrotational congruence u 0 the components of Weyl tensor are reduced to where the magnitude of the shear tensor σ 2 is defined by (41).  [55]. This formalism was revisited in [56] where the authors discuss their applications in detail.
Bianchi identities for Weyl curvature tensor are presented in [56] in a fully expanded form, as they are given a central role in the extended formalism. By specializing the general 1 + 3 dynamical equations it was illustrated how a number of interesting problems can be obtained. In particular, the simplest choices of spatial frames for spatially homogeneous cosmological models, locally rotationally symmetric spacetime geometries, cosmological models with an Abelian isometry group G 2 , and "silent" dust cosmological models were discussed.

Specialization for spherically symmetric models
In the "resource" paper [58] the 1 + 3 orthonormal frame formalism (as developed in [55,56]) was specialized to write down the evolution equations for spherically symmetric models as a well-posed system of first order partial differential equations (PDEs) in two variables. This "resource" paper reviews a number of wellknown results properly cited in [58] and they are simply gathered together because of their functionality and in this context it serves to define all of the quantities. Therefore, we refer researchers interested in the formalism to reference [58] and references therein, and present essential equations now. The metric for the spherically symmetric models is given by where N, e 1 1 , and e 2 2 are functions of t and r and N is the lapse function. The Killing vector fields (KVF) in spherically symmetric spacetime are given by [57] ∂ ζ , cos ζ ∂ ϑ −sin ζ cot ϑ ∂ ζ , sin ζ ∂ ϑ +cos ζ cot ϑ ∂ ζ . Frame vectors in coordinate form are where e 3 3 = e 2 2 / sin ϑ . We see that frame vectors e 2 and e 3 tangent to the spheres are not group-invariant because commutators [e 2 , ∂ ζ ] and [e 3 , ∂ ζ ] are zero but not with the other two Killing vectors. Frame vectors e 0 and e 1 orthogonal to the spheres are group-invariant and the correspondingly commutator reads We use the symbol e a (in reference [55] ∂ ∂ ∂ a is used ) to denote the action of the vector field e a on a scalar f as a differential operator. Geometric objects of the 1 + 3 formalism (kinematic variables, spatial commutations functions) as well as the matter components are deduced as follows. First, we define the four velocity vector field by u = e 0 representing the congruence of worldlines. The representative length (t, r) along worldlines of u = e 0 describing the volume expansion (contraction) behavior of the congruence is reduced to [56] e 0 (t, r) where the Hubble parameter H(t, r) is brought to H(t, r) := − 1 3 e 0 ln e 1 1 (t, r)(e 2 2 (t, r)) 2 .
The anisotropic parameter σ + (t, r) in σ αβ is defined by The acceleration vector is calculated asu a = u a;b u b obtaining only one non-zero component in terms of the spatial derivatives of the lapse function given bẏ u 1 = e 1 ln N.
We have restrictions on spatial commutation functions (1-index objects and 2-index symmetric objects in eq. The dependence of a 2 and n 13 on ϑ is due to the fact that the chosen orthonormal frame is not groupinvariant. However, this is not a concern since dependence on ϑ will be hidden. On the matter components we have restrictions as follows, q α = (q 1 , 0, 0), π αβ = diag(−2π + , π + , π + ).
The frame rotation Ω αβ is zero. Furthermore, n 13 only appears in equations together with e 2 n 13 in the form of the Gaussian curvature of the spheres 2 K := 2(e 2 − 2n 13 )n 13 , which simplifies to Thus, the dependence on ϑ is hidden in the equations. In the following we will use 2 K as a dynamical variable instead of e 2 2 . Spatial curvatures also simplify to 3 S αβ = diag(−2 3 S, 3 S, 3 S), with 3 R and 3 S given by Weyl curvature components simplify to with E + given by To simplify notation we will write 2 K,u 1 , a 1 as K,u, a.
So far, there are no evolution equations for N, p and π + . They need to be specified by a temporal gauge (for N) and by a fluid model (for p and π + ). Recall that the total energy density and total isotropic pressure of the matter fields are µ and p, respectively.

Special cases with extra Killing vectors
Spherically symmetric models with more than three KVF are either spatially homogeneous or static. Let us discuss the spatially homogeneous cosmological models. Spatially homogeneous spherically symmetric models consist of two disjoint sets of models, KS models and FLRW models. Static and self-similar spherically symmetric models have been studied in [58,[98][99][100][101][102].

The Kantowski-Sachs models
Spatially homogeneous spherically symmetric models (that have four Killing vectors with the fourth being ∂ r ) are the so-called KS models [57]. The metric (60) simplifies to i.e., N, e 1 1 , and e 2 2 are now independent of r. Spatial derivative terms of type e 1 (·) in eqs. (69)-(70a) vanish and, as result, a = 0 =u. Sinceu = 0, the temporal gauge is synchronous and we can set N to any positive function of t. Spatial curvatures are given by 3 R = 2K, 3 S + = 1 3 K. Constraint (70b) restricts the source by q 1 = 0. Meanwhile, functions p, and π + are still unspecified. Evolution equations (69) for KS models with unspecified source reduce to The remaining constraint equation (70a) reduces to where we have substituted in eq. (70a) the relation 1 2 3.2.1.1 Kantowski-Sachs models for perfect fluid and homogeneous scalar field. In equations (72) and in the restriction (73) we can replace the expressions π + = 0,Λ = 0, µ = 1 Assuming that the energy-momentum of the scalar field and matter are separately conserved and setting N ≡ 1, we obtain the following equations for KS metric for perfect fluid and homogeneous scalar field Again, in eq. (79) we have substituted the relation 1 2 3 R = K valid for KS metric. As commented before, when u = ∂ /∂ t the dot derivative of a scalar f , given bẏ ∂t , denotes the usual time derivative.

The FLRW models
Spatially homogeneous spherically symmetric models, that are not KS, are FLRW models (with or without Λ ). The source must be a comoving perfect fluid (or vacuum). The metric has the form with f (r) = sin r, r, sinh r, for closed, flat, and open FLRW models, respectively.
(t) is the scale factor of the universe. The frame coefficients are given by e 1 , the temporal gauge is synchronous and we can set N to any positive function of t. The Hubble scalar H = e 0 ln (t) is also a function of t. For the spatial curvatures, 3 S vanishes because eq. (81) implies e 1 a = K, which is consistent with the fact that the frame vector e 1 is not group-invariant. The evolution equation (69d) for σ + and the constraint (70b) imply that π + = 0 = q 1 , i.e., the source is a comoving perfect fluid with unspecified pressure p. Also, note that µ and p only depend on t and p is not specified yet. From eq. (66) the Gaussian curvature of the two spheres is K = −2 f −2 . Meanwhile, from eq. (67) we obtain On the other hand, for closed (k = 1), flat (k = 0), and open FLRW (k = −1), respectively. We used k to indicate the choice of f in eq. (81). Then, The evolution equations simplify to The constraint becomes where we have substituted Minkowski spacetime (Λ = 0, k = 0) which is also static.
1. The vacuum model with Λ > 0, k = 1 is past asymptotic to a model with negative H, and it is future asymptotic to de Sitter model with positive H. Indeed, for Λ > 0 we have 3H 2 = −1/a 2 + Λ . So, the scale factor has to satisfy a ≥ a crit = 3/Λ . Additionally, the scale factor a is not monotonic because H changes its sign. The orbits generated at a negative value of H with a > a crit are such that the scale factor a decreases monotonically until reaching the value a crit = 3/Λ where H becomes zero. Then, H becomes positive because it is continuous and again the scale factor a starts growing from the smallest critical value a crit to infinite at an exponential rate (in a de Sitter phase).  (83) and (84) with k = 1, 0, −1 for closed, flat, and open FLRW, respectively.

Averaging scalar field cosmologies
For KS and positively curved FLRW metrics the Hubble scalar is not monotonic. This means that H cannot be used as a time-dependent perturbation function as in references [46,47]. However, the function (12) is monotonic in a finite time-interval before changing monotony. The region where the perturbation parameter changes its monotony can be analyzed by a discrete symmetry and by introducing the variable T = D/(1 + D) that brings infinite to a finite interval. For KS the normalization factor (12) becomes where K denotes Gaussian spatial curvature of the 2spheres. Using the calculation (82) we obtain for closed FLRW the normalization factor

Kantowski-Sachs
We define normalized variables where D, defined by (86), is the dominant quantity in eq. (79). The function Q is the Hubble factor normalized with D. Such a solution is classified as contracting if Q < 0, since H and Q have the same sign due to D > 0. The function Σ is a measure of the anisotropies of the model normalized with D. When Σ → 0, the solution isotropizes. Φ is the "phase" in the amplitude-phase transformation (21) and r is defined in (22). The function Ω 2 is the total energy of the harmonic oscillator (at the minimum) normalized with 3D 2 . Indeed, ω 2 r 2 /2 represents the total energy density of the pure harmonic oscillator with potential ω 2 2 φ 2 . This is the energy of the oscillator when oscillations are smoothed out and the scalar field approaches its minimum value.
Recall that the potential (14) Then, we obtaiṅ System (89) is invariant for the simultaneous change Assuming ω 2 > 2µ 2 and setting f = bµ 3 where with the time averaging (24) we obtain the time-averaged systeṁΩ Proceeding in an analogous way as in references [103,104] we implement a local nonlinear transformation where D is the normalization factor and its evolution equation is given by (99). Taking time derivative in both sides of (100) with respect to t we obtaiṅ where D x 0 g(D, x 0 ,t) is the Jacobian matrix of g(D, x 0 ,t) with respect to the vector x 0 . The function g(D, x 0 ,t) is conveniently chosen. By substituting eq. (92), which can be written aṡ along with eqs. (99) and (100) in eq. (102) we obtain where I 4 is the 4 × 4 identity matrix. Then, we obtaiṅ Using eq. (99), we haveḊ = O(D 2 ). Hence, The strategy is to use eq. (106) for choosing conveniently ∂ ∂t g(0, x 0 ,t) to prove thaṫ is unknown at this stage. By construction we neglect the dependence of ∂ g i /∂t and g i on D, i.e., we assume g = g(x 0 ,t) because the dependence of D is dropped out along with higher order terms in eq. (106). Next, we solve a partial differential equation for g(x 0 ,t) given by: where we have considered x 0 and t as independent variables. The right hand side of eq. (108) is almost periodic with period L = 2π ω for large times. Then, implementing the average process (24) on right hand side of eq. (108), where slow-varying dependence of quantities Defining the average (109) is zero so that g(x 0 ,t) is bounded. Finally, eq. (107) transforms tȯ and eq. (108) is simplified to Theorem 1 Let be defined the functionsΩ ,Σ ,Q,Φ, and D satisfying time-averaged equations (95), (96), (97), (98), and (99). Then, there exist continuously differentiable functions g 1 , g 2 , g 3 , and g 4 such that Ω , Σ , Q, and Φ are locally given by eq. According to eq. (93) (eq. (99), respectively), we have that D is monotonic decreasing when 0 respectively). Unfortunately, since these regions of full system or averaged system phase space are not invariant for the flow, the monotonicity of D is not guaranteed for all t. (89) and for the averaged equations (95), (96), (97), (98), and (99). Hence, although D(t) remains close to zero for t < t * , where t * satisfiesḊ(t * ) = 0, once the orbit crosses the initial region, D changes its monotony and it will strictly increase without bound for t > t * . Hence, Theorem 1 is valid on a time scale tD = O(1).

FLRW metric with positive curvature
In this section we will study the model with FLRW metric with positive curvature: For FLRW metric with positive curvature the field equations are obtained by setting k = +1 in eqs. (85). Using similar variables as in eqs. (88) with Σ = 0 and replacing the normalization factor D = H 2 + 3 R 6 , where 3 R denotes the 3-Ricci curvature of spatial surfaces calculated in (82) by eq. (87), we obtain the systeṁ The fractional energy density of matter Ω m := ρ m 3H 2 = ρ m 3Q 2 D 2 is parametrized by the equation Setting f = bµ 3 ω 2 −2µ 2 > 0, we obtain the series expansion near D = 0 where the vector function is defined as Systems (114) and (117) are invariant for the simul- the time averaging (24), we obtain for γ = 1 the following time-averaged systeṁΩ

Qualitative analysis of time-averaged systems
According to Theorem 1, in KS metrics and positively curved FLRW models the function (86) plays the role of a time-dependent perturbation function controlling the magnitude of error between solutions of the full and time-averaged problems with the same initial conditions as D → ∞. Thus, oscillations are viewed as perturbations as far as D is bounded. In the time-averaged system Raychaudhuri equation decouples. Therefore, the analysis of the system is reduced to study the corresponding time-averaged equations.

Kantowski-Sachs
With time variable η defined by d f dη = 1 D d f dt the timeaveraged system (95), (96), (97), (98) transforms to where we have definedΩ m as and it was interpreted as the time-averaged values of Ω m := ρ m 3H 2 . Then, the phase space is Furthermore, we have the auxiliary equationṡ Evaluating the averaged values Q =Q, Σ =Σ at eqs.
(126) and integrating the resulting equations, approximated solutions for e 1 1 and K as functions of t are obtained.
Recall that a set of non-isolated singular points is said to be normally hyperbolic if the only eigenvalues with zero real parts are those whose corresponding eigenvectors are tangent to the set. Since by definition any point on a set of non-isolated singular points will have at least one eigenvalue which is zero, all points in the set are non-hyperbolic. However, a set which is normally hyperbolic can be completely classified according to its stability by considering the signs of eigenvalues in the remaining directions (i.e., for a curve in the remaining n − 1 directions) (see [31], pp. 36).
In the special case γ = 1 there exist two lines of equilibrium points which are normally hyperbolic Recall, the subindex ± indicates whether they correspond to contracting ("−") or expanding ("+") solutions. A solution is classified as expanding if Q > 0 since H and Q have the same sign due to D > 0.
The equilibrium points of the guiding system (123a), (123b), (123c) are ii) It is nonhyperbolic for γ = 2 (contained in L − ). For P 1 we obtaiṅ Imposing initial conditions where t = 0, τ(0) = 0 is the current time and D 2 0 = H 2 0 + 1 3c 1 , we obtain by integration The asymptotic metric at P 1 is given by This point represents a non-flat LRS Kasner ( Sect. 6.2.2 and Sect. 9.1.1 (2)). This solution is singular at finite time t 0 = 1 3D 0 and is valid for t > t 0 .
Evaluating at P 9 for 0 ≤ γ < 2 3 we obtain the following: Imposing initial conditions (127), we obtain by integration The line element (71) becomes This solution is singular at finite time t 2 = 1 D 0 4 3γ − 1 and is valid for t > t 2 . For γ = 2 we havė Imposing initial conditions (127), we obtain by integration This solution is singular at finite time t 3 = 1 3D 0 and is valid for t > t 3 .
ii) It is nonhyperbolic for γ = 2 3 or γ = 2. Evaluating at P 10 for 0 ≤ γ < 2 3 we obtain the following: Imposing initial conditions (127), we obtain by integration , The line element (71) becomes For γ = 2 we havė Imposing initial conditions (127), we obtain by integration , To study the dynamics at the invariant boundaryΣ 2 + Ω 2 = 1 which corresponds to vacuum solutions, we introduce cylindrical coordinates The dynamics on the invariant surface is given by In Table 1, the equilibrium points of system (140) are presented. The eigenvalues are obtained evaluating the linearization matrix of (140) at each point.
In Figure 1 orbits in the invariant set (θ ,Q) with dynamics given by (140) are presented. There it is illustrated that P 7 and P 8 are saddle points. In Figure 2 some orbits in the phase space of the guiding system (123a), (123b), (123c) for γ = 0 corresponding to the CC are presented. In Figure 2(a) orbits in the phase space (Σ ,Q,Ω ) are displayed. In Figure 2(b) orbits in the invariant setΩ = 0 are shown. Points P 2 , P 4 , and P 5 are early-time attractors. Points P 1 , P 3 , and P 6 are late-time attractors. Points P 7 , P 8 , P 9 , and P 10 are saddle points. In Figure 3 some orbits of the phase space of the guiding system (123c) for γ = 2 3 are presented. In Figure  3(a) some orbits in the phase space (Σ ,Q,Ω ) are displayed. In Figure 3 with P 9 and P 6 coincides with P 10 . They are nonhyperbolic. Points P 2 and P 4 are early-time attractors. Points P 1 and P 3 are late-time attractors. Points P 5 , P 6 , P 7 , and P 8 are saddle points.
In Figure 4 some orbits in the phase space of the guiding system (123a), (123b), (123c) for γ = 1 which corresponds to dust are displayed. In Figure 4(a) some orbits in the phase space (Σ ,Q,Ω ) are presented. In Figure 4(b) some orbits in the invariant setΩ = 0 are presented. Points P 2 and P 4 are early-time attractors. Points P 1 and P 3 are late-time attractors. Points P 5 , P 6 , P 7 , and P 8 are saddle points.
In Figure 5 some orbits in the phase space of guiding system (123a), (123b), (123c) for γ = 4 3 corresponding to radiation are presented. In Figure 5(a) some orbits in the phase space (Σ ,Q,Ω ) are displayed. In Figure  5(b) some orbits in the invariant setΩ = 0 are shown. In both diagrams 5 and 6 points P 2 and P 4 are earlytime attractors. Points P 1 and P 3 are late-time attractors. Points P 5 , P 6 , P 7 , and P 8 are saddle points. Points P 9 and P 10 do not exist.
In Figure 6 some orbits in the phase space of guiding system (123a), (123b), (123c) for γ = 2 which corresponds to stiff matter are presented. In Figure 6(a) some orbits in the phase space (Σ ,Q,Ω ) are displayed. In Figure 6(b) some orbits in the invariant setΩ = 0 are shown. The line connecting points P 2 , P 4 , and P 6 (denoted by a dashed blue line) is invariant, and it is the early-time attractor. The line connecting points P 1 , P 3 , and P 5 (denoted by a dashed blue line) is invariant and it is the late-time attractor. Points P 7 and P 8 are saddle points. Points P 9 and P 10 do not exist.

late-time behavior in the reduce phase space
Now, we study the dynamics in the reduced phase space x = (Ω , Σ , Q) where the effect of D in the dynamics was neglected. The results from the linear stability analysis combined with Theorem 1 lead to the following local results. For global results when D ∈ [0, ∞) see section 6.
ii) It is a saddle for 2 3 < γ < 1. iii) It is nonhyperbolic for γ = 2 3 and γ = 1. iv) It is a sink for 1 < γ ≤ 2. This equilibrium point is related to the isotropic point P 5 of KS. The asymptotic metric at P 5 is given by The corresponding solution is a flat matter-dominated FLRW contracting solution withΩ m = 1.      ii) It is a saddle for 2 3 < γ < 1. iii) It is nonhyperbolic for γ = 2 3 and γ = 1. iv) It is a source for 1 < γ ≤ 2. This equilibrium point is related to the isotropic point P 6 of KS. The asymptotic metric at the equilibrium point P 6 is i) It is a sink for 0 ≤ γ < 1.
iii) It is a saddle for 1 < γ ≤ 2.   This equilibrium point is related to the isotropic point P 7 of KS. The line element (71) becomes i) It is a source for 0 ≤ γ < 1.
iii) It is a saddle for 1 < γ ≤ 2. This equilibrium point is related to the isotropic point P 8 of KS. The line element (71) becomes Hence, the equilibrium point can be associated with Einstein-de Sitter solution ([55], Sec 9.1.1 (1)) with γ = 1. It is an expanding solution.
In the special case γ = 2 3 there is one line of equilibrium points which are normally hyperbolic, M 0 : (Ω ,Q) = 0,Q c with eigenvalues −Q c 2 , 0 , which is a sink for Q c > 0 and a source forQ c < 0.
The first equation of guiding system (119), (120), (121) becomes trivial up to second order in the Dexpansion when γ = 1. Using Taylor expansion up to the fourth order in D the following time-averaged system is obtained From eq. (147b) we havē Substituting back eq. (149) in eq. (147) we have Using the method of the integrating factor we define v = e Q (η)dη D(η), to obtain The integration constant can be absorbed in the indefinite integral. Then, we have Substituting eq. (153) in eq. (149) we obtain Then, substituting eqs. (153) and (154) in eq. (147c) we obtain the quadraturē Substituting eqs. (153), (154) in eq. (147a) we obtain under assumptions (148) the differential equation Summarizing, system (147) admits the first integral (153), (154), (155), whereQ satisfies the differential equation (156). To analyze the asymptotic behavior of the solutions of eq. (156) we introduce the variables x =Q(η), y = Q (η) to obtain The origin is an equilibrium point with eigenvalues i √ 2, −i √ 2 . The dynamics of system (157) in the coordinates (x, y/ 1 + y 2 ) is presented in Figure 8 where the origin is a stable center. According to the analysis at first order of the guiding system (119), (120), (121), Einstein's static universe has coordinates E : (Ω ,Q) = 3γ−2 3γ−3 , 0 . When γ = 1, this point satisfiesΩ → ∞. However, by taking γ = 1 the equation forΩ is decoupled andΩ only takes arbitrary constant values. On the other hand, since γ = 1, we required higher order terms in Taylor expansion for obtaining system (147). Since extraΩ -coordinate is excluded from Figure 8  For global results when D ∈ [0, ∞) see section 6.

Regular dynamical system on a compact phase space for Kantowski-Sachs and closed FLRW models
According to Remark 1, Theorem 1 is valid on a finite time scale where D, given by eq. (86), remains close to zero, but at a critical time t * we have D (t * ) = 0 and D becomes strictly increasing when t > t * such that lim t→∞ D(t) = ∞. A lower bound for t * is estimated as

Kantowski-Sachs metric
In this section we analyze qualitatively the time-averaged system (123) as D → ∞ by introducing the variable that maps [0, ∞) to a finite interval [0, 1). Therefore, the limit D → +∞ corresponds to T = 1 and the limit D → 0 corresponds to T = 0. Then, we have the guiding system (123a), (123b), (123c) extended with equation We are interested in late-time or early-time attractors, and in discussing relevant saddle equilibrium points of the extended dynamical system (123a), (123b), (123c) (159). In this regard, we have the following results.
In the special case γ = 0, there exist four lines of equilibrium points which are normally hyperbolic:  3 2 , 0, 4 − 2Σ c is a source forΣ c < 2. As before, the subindex ± indicates if they correspond to contracting ("−") or to expanding ("+") solutions. Superindex ∞ refers to solutions with D → ∞. For KS model the extended phase is four dimensional. To illustrate the stability of the aforementioned normally hyperbolic lines we examine numerically some invariant sets. In Figure 9 the dynamics at the invariant setΩ = 0 corresponding to vacuum solutions is represented in the compact space (Σ , D/(1+D),Q) for γ = 0 and γ = 2. For γ = 1 andΣ = 0, the model is reduced to the closed FLRW metric with dust. The stability of the lines K ± and K ∞ ± on this invariant set corresponding to isotropic solutions is illustrated in Figure 10    Also, we refer to Figures 13,14,15,16,17,18,19, and 20, where projections of some orbits of the full system and the averaged system are presented. In these figures it is numerically confirmed the result of Theorem 1 for KS metric. That is to say, solutions of the full system (blue lines) follow the track of the solutions of the averaged system (orange lines) for the whole D-range.

FLRW metric with positive curvature
Using the variable (158) the time-averaged system for γ = 1 (141) becomes the guiding system (141a) and  (141b) extended with equation We are interested in late-time or early-time attractors and in discussing the relevant saddle equilibrium points of the extended system (141a), (141b) and (160). In this regard, we have the following results.
In the special case γ = 0 there exist two lines of equilibrium points which are normally hyperbolic, say,  Figure 10 the dynamics of the averaged system (141) for γ = 0, 2 3 and of the averaged system (147) for γ = 1 is represented in the compact space (Q, D/(1+D),Ω ). In Figures 13, 14, 15, 16, 17, 18, 19, and 20 projections of some solutions of the full system (89)  Now, we describe the regime Ω > 1. We define Substituting the form of potential (14) in constraint (85e) we obtain for closed FLRW models Then, Using definition (87) we obtain   11 Some solutions of the full system (114) (blue) and timeaveraged system (141) (orange) for the FLRW metric with positive curvature (k = +1) when γ = 4 3 . We have used for both systems the initial data sets presented in Table 2. Table 2 Here we present three initial data sets for the simulation of full system (114) and time-averaged system (141) for FLRW metric with positive curvature (k = +1) and for γ = 4 3 (radiation) and γ = 2 (stiff fluid). All the conditions are chosen in order to fulfill the inequalities Ω > 1 and 0 ≤ Q ≤ 1.
The function U(φ ) defined by eq. (161) satisfies U (0) = U (0) = U (0) = 0, and U (iv) (0) = − 2 f 2 , which implies that φ = 0 is a global degenerated maximum of order 2 for U(φ ). Therefore, U(φ ) ≤ 0. Then, from ω 2 > 2µ 2 and (164) it follows that Ω , Σ 2 and ρ m 3D 2 can be greater than 1. This implies that we can have solutions with Ω > 1 preserving the non-negativity of the   12 Some solutions of the full system (114) (blue) and timeaveraged system (141) (orange) for the FLRW metric with positive curvature (k = +1) when γ = 2. We have used for both systems the initial data sets presented in Table 2. energy densities. That is, even if U(φ ) ≤ 0, we have 2 φ 2 is dominant over the last term in eq. (164). This interesting behavior when Ω > 1 can be seen in Figures 7(d), 8, 11 and 12, where we show the solutions of the full system (114) and the time-averaged system (141) considering the values γ = 4 3 and γ = 2 and using as initial conditions the three data sets presented in Table 2. This dynamical behavior related to spiral tubes has been presented before in the literature in [113] and it is related to the fact that the line of equilibrium points E c (representing Einstein's static universes) has purely imaginary eigenvalues. In [113] a comprehensive dynamical description for closed cosmologies when the matter source admits Einstein's static model was pre-sented. Moreover, theorems about the global asymptotic behavior of solutions were established. Results in [113] and [114] disprove claims of non-predictability and chaos for models close to Einstein's model given in [115][116][117][118].
To illustrate the existence of spiral tubes we integrate the full system (114) and the time-averaged system (141) using the fixed constants µ = and γ = 2 for the barotropic index (cases where E c exists), and we use as initial conditions the three data sets presented in Table 2. In Figures 11 and 12

Conclusions
This is the last paper of the "Averaging generalized scalar field cosmologies" research program. We have used asymptotic methods and averaging theory to explore the solutions space of scalar field cosmologies with generalized harmonic potential (13) in vacuum or minimally coupled to matter. Different from references [46,47], here we have studied systems where Hubble scalar is not monotonic, but the systems admit a function D given by eq. (12) playing the role of a timedepending perturbation parameter, which is decreasing for a finite time scale t < t * where D (t * ) = 0. For t > t * monotony of the quantity D changes and this parameter increases without bound. We have proved Theorem 1 which states that late-time attractors of full and time-averaged systems are the same when the quantity D tends to zero. More specifically, according to Theorem 1 for KS metrics and positively curved FLRW models, the quantity D controls the magnitude of error between full time-dependent and time-averaged solutions as D → 0. Therefore, the analysis of the system is reduced to study the corresponding time-averaged equation as D → 0. However, for KS metric the initial region 0 < Q < 1, Σ 2 + Ω 2 < 1, Σ (1 − Q 2 + 3QΣ ) > 0 (and for closed FLRW the initial region Q > 0, respectively) is not invariant for the full system (89) and for the time-averaged equations (95), (96), (97), (98), and (99). According to Remark 1, Theorem 1 is valid on a time scale t < t * where D remains close to zero, but D(t) changes its monotony at a critical time t * becoming monotonic increasing without bound.
We have formulated Theorems 2 and 3 concerning the late-time behavior of our model valid when the evolution equation for D is decoupled, whose proofs are based on Theorem 1 and linear stability analysis. Hence, we can establish the stability of a periodic solution as it exactly matches the stability of the stationary solution of the time-averaged equation. In particular, for KS metric the local late-time attractor of full system (89) and time-averaged system (123) (where the evolution equation for D is decoupled) are the following.
We have commented that, although for t < t * , D(t) remains close to zero, once the orbit crosses the initial region, D changes its monotony, and it becomes strictly increasing without bound. Hence, Theorem 1 is valid on a time scale tD = O(1). To investigate the region D → ∞ we have used the transformation of coordinates (158) that maps [0, ∞) to a finite interval [0, 1). Therefore, the limit D → +∞ corresponds to T = 1 and the limit D → 0 corresponds to T = 0. This defines a regular dynamical system over a compact phase space that allows to obtain global results. We have studied the stability of the fixed points in a compactified phase space. These numerical results support the claim that late-time attractors in the extended phase space (x, T ), where T = D 1+D and x = (Ω , Σ , Q) for both the original system and the time-averaged are the same for KS. When the stability of the equilibrium point of the time-averaged is analyzed in extended phase space, we find for KS metric that the extra variable T introduces equilibrium points "at infinity", P ∞ 1 which is a non-flat LRS Kasner solution and P ∞ 3 which is Taub (flat LRS Kasner). They are contracting solutions and sink for 0 ≤ γ < 2 in the extended (global) phase space. Their analogous points P 1 and P 3 (with D = H = T = 0) become saddle along the T -axis in the extended phase space. The only equilibrium point that remains a sink for KS for 0 ≤ γ < 2 3 in the extended phase space is P 6 . Figures 13, 14, 15, 16, 17, 18, 19, and 20 are a numerical confirmation that main Theorem 1 presented in Section 4 is fulfilled for the KS metric. That is to say, the solutions of the full system (blue lines) follow the track of the solutions of the averaged system (orange lines) for the whole D-range. On the other hand, local late-time attractors of full system (114) and time-averaged system (141) (where the evolution equation for D is decoupled) for closed FLRW metric with positive curvature are the following.
(iii) The equilibrium point P 7 if 0 ≤ γ < 1. The equilibrium point can be associated with Einstein-de Sitter solution.
When the stability of the equilibrium point of the timeaveraged extended phase space for closed FLRW metric is analyzed in the extended phase space (x, T ), where x = (Ω , Q), for closed FLRW we find that the extra variable T introduces equilibrium points "at infinity", P ∞ 5 which is a sink for 1 < γ ≤ 2 and P ∞ 7 which is a sink for 0 ≤ γ < 1. As for KS, the only equilibrium point that remains a sink for KS for 0 ≤ γ < 2 3 in the extended phase space is P 6 . In Figures 21, 22, 23, and 24 we have presented projections of some solutions of the full system (114) and time-averaged system (141) for γ = 1. Also, system (147) for γ = 1 in the (Q, D/(1+D), Ω 2 ) space are presented with its respective projection when D = 0. Figures 21 and 25 show how solutions of the full system (blue lines) follow the track of solutions of the averaged system (orange lines) for the whole D-range. Figures  22, 23, 24, and 26 are evidence that main theorem presented in Section 4 is fulfilled for FLRW metric with positive curvature (k = +1) only when D is bounded. Namely, the solutions of the full system (blue lines) follow the track of solutions of the averaged system (orange lines) for the time interval tD = O(1). However, when D becomes infinite (T → 1) and for γ ≥ 1, the solutions of full system (blue lines) depart from solutions of averaged system (orange lines) as D becomes large. Then, different from KS, for the full system and given γ ≥ 1 the orbits (blue lines) do not follow the track of solutions of the averaged system, while P ∞ 8 and P ∞ 7 are early and late-time attractor, respectively, as D → ∞. This is a rather different behavior from time-averaged system, where they are saddle. This can be anticipated because when D becomes large, the approximation obtained under the assumption of small D fails. Additionally, for closed FLRW we have found by numerical tools the existence of spiral tubes confined in a finite region of the phase space when the line of equilib-rium points E c (representing Einstein's static universes) exist. This kind of dynamical structures have been presented before in reference [113] and they exist for any matter source that admits Einstein's static model. Results in the line of [113] are of interest since they disprove the claims of non-predictability and chaos for models close to Einstein's model related to the existence of infinitely many homoclinic orbits whose α and ω-limits are the same periodic orbit producing chaotic sets in the whole state space. Thus, a set of models in a neighborhood of Einstein's model were claimed to be unpredictable and characterized by "homoclinic chaos" [115][116][117][118]. However, the asserted "homoclinic phenomena", if they occur at all, must be confined to narrow regions of the phase space [113] (see also [114]). Now, we summarize the results of the "Averaging generalized scalar field cosmologies" research program. In reference [46] it was proved that for LRS Bianchi III the late-time attractors are a matter-dominated flat FLRW universe if 0 ≤ γ ≤ 2 3 (mimicking de Sitter, quintessence or zero acceleration solutions), a mattercurvature scaling solution if 2 3 < γ < 1, and Bianchi III flat spacetime for 1 ≤ γ ≤ 2. For FLRW metric with k = −1 the late-time attractors are flat matterdominated FLRW universe if 0 ≤ γ ≤ 2 3 and Milne solution if 2 3 < γ < 2. In all metrics, matter-dominated flat FLRW universe represents quintessence fluid if 0 < γ < 2 3 . In reference [47] for flat FLRW and LRS Bianchi I metrics it was obtained that late-time attractors of full and time-averaged systems are given by flat matterdominated FLRW solution and Einstein-de Sitter solution. It is interesting to note that for FLRW with negative or zero curvature and for Bianchi I metric when the matter fluid corresponds to a CC, H asymptotically tends to constant values depending on initial conditions. That is consistent with de Sitter expansion. In addition, for FLRW models with negative curvature for any γ < 2 3 and Ω k > 0, Ω k → 0, or when γ > 2 3 and Ω k > 0 the universe becomes curvature dominated (Ω k → 1). For flat FLRW and dust background from the qualitative analysis performed in paper [47] we have that lim τ→+∞Ω (τ) = const., and lim τ→+∞ H(τ) = 0. Also, it was numerically proved that as H → 0, the values ofΩ give an upper bound for the values Ω of the original system. Therefore, by controlling the error of the time-averaged higher order system the error of the original system can also be controlled. Finally, in the present research we have proved that in KS metric the global late-time attractors of full and time-averaged systems are two anisotropic contracting solutions, P ∞ 1 which is a non-flat LRS Kasner and P ∞ 3 which is a Taub (flat LRS Kasner) for 0 ≤ γ < 2, and P 6 which is a matter-dominated flat FLRW uni-verse if 0 ≤ γ ≤ 2 3 (mimicking de Sitter, quintessence or zero acceleration solutions). For FLRW metric with k = +1 global late-time attractors of time-averaged system are P ∞ 5 which is a flat matter-dominated contracting solution that is a sink for 1 < γ ≤ 2, P 6 which is a matter-dominated flat FLRW universe mimicking de Sitter, quintessence or zero acceleration solutions if 0 ≤ γ ≤ 2 3 , and the point P ∞ 7 which is an Einstein-de Sitter solution for 0 ≤ γ < 1 and large t. However, when D becomes infinite (T → 1) and for γ ≥ 1 solutions of the full system (blue lines) depart from solutions of averaged system (orange lines) as D is large. Then, different from KS, for the full system and given γ > 1 the orbits (blue lines) do not follow the track of the averaged system and P ∞ 8 and P ∞ 7 are the early and late-time attractor, respectively, as D → ∞. This is a rather different behavior from the time-averaged system, where they are saddle. Therefore, this analysis completes the characterization of the full class of homogeneous but anisotropic solutions and their isotropic limits with exception of LRS Bianchi V. Our analytical results were strongly supported by numerical results in Appendix B. We have shown that asymptotic methods and averaging theory are powerful tools to investigate scalar field cosmologies with generalized harmonic potential.

Lemma 5 (Mean value Theorem)
Let U ⊂ R n be open, f : U → R m continuously differentiable, and x ∈ U, h ∈ R m vectors such that the line segment x + z h, 0 ≤ z ≤ 1 remains in U. Then, we have where Df denotes the Jacobian matrix of f and the integral of a matrix is understood as a componentwise.
From eq. (93) D is a monotonic decreasing function of t whenever 0 < Q < 1, The last restriction holds, in particular, by choosing Q > 0, Σ > 0. In the last case we define the bootstrapping sequences .
This sequence, however, is finite; that is, t n < t * with t * such that lim t→t * D (t) = 0. If Q(t n ) < 0 or Σ (t n ) < 0 for some n we stop the integration because D changes its monotony to become a monotonic increasing function. However, (t, η, Σ , Q, Φ) → (−t, −η, −Σ , −Q, −Φ) leaves invariant the system. Therefore, the solution is completed by using the above symmetry. Given expansions (100), system (106) becomeṡ Explicit expressions for g i obtained by integration of (A.7) are given by where we have set four integration functions C i (D, Ω 0 , Σ 0 , Q 0 , Φ 0 ), i = 1, 2, 3, 4 to zero. The functions g i , i = 1, 2, 3, 4 are continuously differentiable such that their partial derivatives are bounded on t ∈ [t n ,t n+1 ].
The second order expansion around D = 0 of system (111) is written aṡ , (A.12) , (A.14) The last row corresponding to∆ Φ 0 was omitted. This vector function with polynomial components in variables (y 1 , y 2 , y 3 ) is continuously differentiable in all its components.
Using the same initial conditions for x 0 andx, we obtain by integration , for all t ∈ [t n ,t n+1 ]. Using Lemma 5 we havē Df denotes the Jacobian matrix off and the integral of a matrix is understood as a componentwise. We denote the matrix elements of A as where a i j are polynomial functions ofΩ , Ω 0 ,Σ , Σ 0 ,Q, Q 0 , and they are explicitly given by where the sup norm of a matrix | a i j | is defined by max{|a i j |, i = 1, 2, 3, j = 1, 2, 3} with a i j given by eqs. (A.21). Define L 1 = 3 max s∈[t n ,t n+1 ] |A(s)|, which is constant by continuity of Ω , Ω 0 ,Σ , Σ 0 ,Q, Q 0 in [t n ,t n+1 ]. Therefore, Using Gronwall's Lemma 4 we have for t ∈ [t n ,t n+1 ], due to t − t n ≤ t n+1 − t n = 1 D n . Then, Furthermore, defining which is finite by continuity ofΩ , Ω 0 ,Σ , Σ 0 , Φ 0 in the closed interval [t n ,t n+1 ]. We obtain from eq. (A.15) that due to t −t n ≤ t n+1 −t n = 1 D n . Finally, it follows that the functions Ω 0 , Σ 0 , Q 0 , Φ 0 andΩ ,Σ ,Q,Φ have the same asymptotics as D n → 0.

Appendix B: Numerical simulation
In this section we present the numerical evidence that supports the main theorem presented in section 4. Precisely, we numerically solve full and time-averaged systems obtained for KS and FLRW with positive curvature metrics. To this purpose an algorithm in the programming language Python was elaborated. This program solves the systems of differential equations using the solve ivp code provided by the SciPy open-source Python-based ecosystem. The used integration method was Radau, which is an implicit Runge-Kutta method of the Radau IIa family of order 5 with relative and absolute tolerances of 10 −4 and 10 −7 , respectively. All systems of differential equations were integrated with respect to η, instead of t, with an integration range of −40 ≤ η ≤ 10 for the original system and −40 ≤ η ≤ 100 for the time-averaged system. All of them were partitioned in 10000 data points. This dynamical behavior related to spiral tubes has been presented before in the literature in [113] and it is related to the fact that the line of equilibrium points E c (representing Einstein's static universes) has purely imaginary eigenvalues. We present numerical results in the same line of [113] in Figures 11 and 12.
As initial conditions we use the seven data set presented in Table 3 for a better comparison of both systems. In Figures 13,14,15,16,17,18,19, and 20 projections of some solutions of the full system (89) and timeaveraged system (123) in the (Σ , D/(1 + D), Ω 2 ) and (Q, D/(1 + D), Ω 2 ) space are presented with their respective projections when D = 0. Figures 13 and 14 show solutions for a fluid corresponding to CC (γ = 0). Figures 15 and 16 show solutions for a fluid corresponding to dust (γ = 1). Figures 17 and 18 show solu- Table 3 Here we present seven initial data sets for the simulation of the full system (89) and time-averaged system (123) for the KS metric. All the conditions are chosen in order to fulfill the inequalities Σ 2 (0) + Ω 2 (0) ≤ 1 and 0 ≤ Q ≤ 1.
Sol.  Table 4 Here we show ten initial data sets for the simulation of full system (114) and time-averaged system (141) for γ = 1 and of system (147) for γ = 1, for the FLRW metric with positive curvature (k = +1). All the conditions are chosen in order to fulfill the inequalities 0 ≤ Ω ≤ 1 and −1 ≤ Q ≤ 1.
Sol. For case γ = 1, for which the differential equation for Ω in time-averaged system (141) becomes trivial, we integrate time-averaged system (147). Independently of the value of γ, we use as initial conditions ten data set presented in Table 4, where data sets I, II, and V II are the symmetrical counterpart with respect to Q of data sets i, ii, and vii. In Figures 21,  22, 23, and 24 we present projections of some solutions of full system (114) and time-averaged system (141) for γ = 1 and of system (147) for γ = 1 in the (Q, D/(1 + D), Ω 2 /(1 + Ω 2 )) space with their respective projection when D = 0. For both systems the same initial data sets from Table 4 13 Some solutions of the full system (89) (blue) and time-averaged system (123) (orange) for the KS metric when γ = 0 in the projection Q = 0. We have used for both systems the initial data sets presented in Table 3.   (89) (blue) and time-averaged system (123) (orange) for the KS metric when γ = 0 in the projection Σ = 0. We have used for both systems the initial data sets presented in Table 3.   15 Some solutions of the full system (89) (blue) and time-averaged system (123) (orange) for the KS metric when γ = 1 in the projection Q = 0. We have used for both systems the initial data sets presented in Table 3.   16 Some solutions of the full system (89) (blue) and time-averaged system (123) (orange) for the KS metric when γ = 1 in the projection Σ = 0. We have used for both systems the initial data sets presented in Table 3.   17 Some solutions of the full system (89) (blue) and time-averaged system (123) (orange) for the KS metric when γ = 4 3 in the projection Q = 0. We have used for both systems the initial data sets presented in Table 3.   18 Some solutions of the full system (89) (blue) and time-averaged system (123) (orange) for the KS metric when γ = 4 3 in the projection Σ = 0. We have used for both systems the initial data sets presented in Table 3.   19 Some solutions of the full system (89) (blue) and time-averaged system (123) (orange) for the KS metric when γ = 2 in the projection Q = 0. We have used for both systems the initial data sets presented in Table 3.   20 Some solutions of the full system (89) (blue) and time-averaged system (123) (orange) for the KS metric when γ = 2 in the projection Σ = 0. We have used for both systems the initial data sets presented in Table 3.   Table 4.   22 Some solutions of the full system (114) (blue) and time-averaged system (147) (orange) for the FLRW metric with positive curvature (k = +1) when γ = 1. We have used for both systems the initial data sets presented in Table 4.  We have used for both systems the initial data sets presented in Table 4.
iv v (b) Projection in the space (Q, Ω 2 ). The black line represent the constraint Q = 0. Fig. 24 Some solutions of the full system (114) (blue) and time-averaged system (141) (orange) for the FLRW metric with positive curvature (k = +1) when γ = 2. We have used for both systems the initial data sets presented in Table 4.   D), Ω 2 /(1 + Ω 2 )). We have used for both systems the initial data sets presented in Table 4.   D), Ω 2 /(1 + Ω 2 )). We have used for both systems the initial data sets presented in Table 4.