Averaging generalized scalar field cosmologies I: locally rotationally symmetric Bianchi III and open Friedmann–Lemaître–Robertson–Walker models

Scalar field cosmologies with a generalized harmonic potential and a matter fluid with a barotropic Equation of State (EoS) with barotropic index γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} for locally rotationally symmetric (LRS) Bianchi III metric and open Friedmann–Lemaître–Robertson–Walker (FLRW) metric are investigated. Methods from the theory of averaging of nonlinear dynamical systems are used to prove that time-dependent systems and their corresponding time-averaged versions have the same late-time dynamics. Therefore, simple time-averaged systems determine the future asymptotic behavior. Depending on values of barotropic index γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} late-time attractors of physical interests for LRS Bianchi III metric are Bianchi III flat spacetime, matter dominated FLRW universe (mimicking de Sitter, quintessence or zero acceleration solutions) and matter-curvature scaling solution. For open FLRW metric late-time attractors are a matter dominated FLRW universe and Milne solution. With this approach, oscillations entering nonlinear system through Klein–Gordon (KG) equation can be controlled and smoothed out as the Hubble factor H – acting as a time-dependent perturbation parameter – tends monotonically to zero. Numerical simulations are presented as evidence of such behaviour.

This research program -named "Averaging Generalized Scalar Field Cosmologies"-has three steps according to three cases of study: (I) Bianchi III and open FLRW model, (II) Bianchi I and flat FLRW model and (III) Kantowski-Sachs (KS) and closed FLRW. This paper is devoted to case I, and cases II and III will be studied in two companion papers [164,165]. The main aspect in the present work is the interaction of KG fields and field equations. The paper is organized as follows. In Sect. 2 we discuss the class of generalized harmonic potentials in which we are interested. In Sect. 3 we introduce the models. In Sect. 4 we apply averaging methods to analyze periodic solutions of a scalar field with self-interacting potentials within the class of generalized harmonic potentials [148]. In particular, in Sect. 4.3 LRS Bianchi III metric is studied. In Sect. 4.4 is investigated FLRW metrics with negative curvature. In Sect. 5 we study averaged systems. Bianchi III metric is studied in Sects. 5.1, and 5.2 is devoted to FLRW metric with negative curvature (open FLRW). Finally, in Sect. 6 our main results are discussed. In Appendix A is proved our main Theorem. In Appendix B center manifold calculations for nonhyperbolic equilibrium points are presented. In Appendix C numerical evidences supporting the results of Sect. 4 are presented.

Generalized harmonic potential
Chaotic inflation is a model of cosmic inflation which takes the potential term V of a hypothetical inflaton field φ as V (φ) = m 2 φ φ 2 2 , the so-called harmonic potential (φ 2interaction) [201][202][203][204]. Whereas other inflationary models assume a monotonic decreasing potential with φ; assuming in an ad hoc way that inflaton field has a large amplitude "at Big Bang", then slowly "roll down" the potential. The idea of [201] is that instead of inflaton rolls down and sits on its potential minimum at equilibrium, quantum fluctuations stochastically ("chaotically") drive it out of its minimum back and forward. Wherever this happens cosmic infla-tion sets in and blows up the region of ambient spacetime in which inflaton happened to fluctuate out of its equilibrium. Relevant experimental results disfavoring φ 2 -interaction are due to [205,206]. These results state that chaotic inflation generically predicts large values of tensor-to-scalar ratio r . In contrast to recent measurements which show low upper bounds on r . Notwithstanding, we investigate variations of φ 2 -potential and we do not refer to tensor-to-scalar ratio issue for the potential (4).
The action integral of interest is It is expressed in a system of units in which 8π G = c = = 1 where L m is the Lagrangian density of matter, R is the curvature scalar, φ is the scalar field, ∇ α is the covariant derivative and the scalar field potential V (φ) of interest in this research is given by It is related but not equal to monodromy potential of [162] used in the context of loop-quantum gravity, which is a particular case of general monodromy potential [163]. In references [148][149][150] were proved that potential of [162,163] for p = 2, say V (φ) = μ 3 φ 2 μ + b f cos φ f , b = 0 is not good to describe the late-time FLRW universe driven by a scalar field because it has two symmetric local negative minimums which are related to Anti-de-Sitter solutions.
Therefore, in [148,149] we have studied the potential obtained by setting μ = √ 2 2 and bμ = 2 in Eq. (2). The potential (3) provides non-negative local minimums which can be related to a late-time accelerated universe. In Sect. 2.4 of [149] a scalar field cosmology with potential (3) nonminimally coupled to matter with coupling function χ = χ 0 e λφ 4−3γ was studied, where λ is a constant and the barotropic index satisfies 0 ≤ γ ≤ 2, γ = 4 3 for FLRW metrics with k = −1, 0 and Bianchi I metric. The late time attractors are associated to equilibrium points with φ = φ * whenever φ * is a local non zero minimum of V (φ). For FLRW metrics, global minimum is unstable to curvature perturbations for γ > 2 3 . Therefore, the result in [107] is confirmed, that for γ > 2 3 the curvature has a dominant effect on late evolution of the universe and it will eventually dominate both perfect fluid and scalar field energy densities. For Bianchi I model, the global minimum with V (0) = 0 is unstable to shear perturbations.  [207]. In [208] axionic dark matter model with a modified periodic potential in the framework of axionic extension of Einstein-aether theory was studied. This periodic potential has minima at φ = nΦ * , n ∈ Z, whereas maxima when n → m + 1 2 are found. Near the minimum when φ = nΦ * + ψ and |ψ| is where m A the axion rests mass. The previous statements justify the study of potential (2), which can be expressed as by introducing an angular frequency ω ∈ R through conditions bμ 3 + 2 f μ 2 − f ω 2 = 0 and ω 2 − 2μ 2 > 0. Their applicability will be revealed in Sect. 4. The generalized harmonic potentials (3) and (4) belong to the class of potentials studied by [127]. Potential (4) has the following generic features: There is a finite number of values φ c = 0 satisfying 2μ 2 φ c + f ω 2 − 2μ 2 sin φ c f = 0, which are extreme points of V (φ). They are local maximums or local minimums depending on whether V (φ c ) The function V has no upper bound but it has a lower bound equal to zero.
The asymptotic features of potential (4) are the following. Near global minimum φ = 0, we have V (φ) ∼ That is, ω 2 can be related to the mass of the scalar field near its global minimum. As φ → ±∞ cosine-correction is bounded, then, This makes it suitable to describe oscillatory behavior in cosmology. Setting f , the potential of [148][149][150] (3) is recovered. Although (3) can be derived as a particular case of our study, cases of interest with f 1 leads to complex frequency ω = | f −1 f |i, in contradiction to ω ∈ R. Therefore, potential (4) will not contain potential studied in [148][149][150], unless we set f > 1 and ω = f −1 f . In

Spatially homogeneous scalar field cosmologies
In general relativity (GR) the SH but anisotropic spacetimes are known as either Bianchi or KS metrics. In Bianchi models, the spacetime manifold is foliated along the time axis with three dimensional homogeneous hypersurfaces. Bianchi spacetimes contain many important cosmological models that have been used to study anisotropies of primordial universe and its evolution towards current observed isotropy [166][167][168][169]. The list includes standard FLRW model in the limit of isotropization; Bianchi III isotropizing to open FLRW models and Bianchi I isotropizing to flat FLRW models. Hubble parameter H is always monotonic for Bianchi I and Bianchi III. For Bianchi I anisotropy decays on time for H > 0 and isotropization occurs [170].
There is an interesting hierarchy in Bianchi models [121,[171][172][173]. In particular, LRS Bianchi I model naturally appears as a boundary subset of LRS Bianchi III model. The last one is an invariant boundary of the LRS Bianchi type VIII model as well. Additionally, LRS Bianchi type VIII can be viewed as an invariant boundary of LRS Bianchi type IX model [174][175][176][177][178][179][180]. Bianchi spacetimes in presence of a scalar field were studied in [181]. It was proved that an initial anisotropic universe isotropizes into a FLRW universe for specific initial conditions if the scalar field potential has a large positive value. An exact solution of field equations for an exponential potential in some particular Bianchi spacetimes has been found in [182][183][184]. These exact solutions lead to isotropic homogeneous spacetimes as it was found in [185,186]. An anisotropic solution of special interest is Kasner spacetime. Kasner solution is essential for the description of BKL singularity when the contribution of Ricci scalar of three-dimensional spatial hypersurface in the field equations is negligible [187]. For other applications of Kasner universe and Bianchi I spacetimes in gravitational physics see [188][189][190][191][192][193][194][195] and references therein. In [196] the conformal algebra of Bianchi III and Bianchi V spacetimes which admit a proper conformal Killing vector were studied. In [197] method of Lie symmetries was applied for the Wheeler-De Witt equation in Bianchi class A cosmologies for minimally coupled scalar field gravity and hybrid gravity in GR. Using these symmetries several invariant solutions were determined and classified according to the form of scalar field potential.

LRS Bianchi III, Bianchi I and Kantowski-Sachs models
Due to the metric element for LRS Bianchi III, Bianchi I and KS models can be written as [198] where e 1 1 , e 2 2 and e 3 3 = √ ke 2 2 / sin( √ kϑ) are functions of t which are components of the frame vectors [200]: e 0 = ∂ t , e 1 = e 1 1 ∂ r , e 2 = e 2 2 ∂ ϑ , e 3 = e 3 3 ∂ ζ . That is, we obtain LRS Bianchi III for k = −1, Bianchi I for k = 0 and KS for k = +1 [62]. Comparing with reference [198] we have settled parameters a = f = 0 and e 1 1 (t) = D 2 (t) −1 , e 2 2 (t) = D 2 (t) −1 in their metric and we have used the identifications (ϑ, ζ ) = (y, z). The line elements for spatially selfsimilar LRS models have been given by Wu in [199]. We only focus on SH but anisotropic class with the exception of SH LRS Bianchi V, that is: LRS Bianchi III, Bianchi I and KS.
It is useful to define a representative length along worldlines of the 4-velocity vector u = ∂ t describing the volume expansion (contraction) of the congruence, denoted (t) and defined bẏ where the Hubble parameter H (t) and the anisotropic parameter σ + (t) are given by Taking variation of (1) for the 1-parameter family of metrics (9) leads to [62]: For modeling matter in our model we use a perfect fluid with a barotropic EoS p m = (γ − 1)ρ m with pressure p m , energy density ρ m and barotropic index γ ∈ [0, 2]. The Gauss curvature of spatial 2-space and 3-curvature scalar are Furthermore, evolution of K iṡ while evolution for e 1 1 is given by [200]: From Eqs. (13), (14) the shear equation is obtained: Equations (12), (13), (14), (18) give the Raychaudhuri equation: Finally, the matter and KG equations are: In this paper we will study LRS Bianchi III model. By convenience we write where g H 2 = dϑ 2 + sinh 2 (ϑ)dζ 2 denotes the 2-metric of negative constant curvature on hyperbolic 2-space. The functions A(t) and B(t), interpreted as scale factors, are defined as A(t) = e 1 1 (t) −1 and B(t) 2 = K (t) −1 .

FLRW models
The general line element for spherically symmetric (SS) models can be written as [200] SH-SS models that are not Kantowski-Sachs are FLRW models where the metric can be written as with f (r ) = sin r, r, sinh r (25) for closed, flat and open FLRW models, respectively. In comparison with metric (23), frame coefficients here are given by e 1 1 = a −1 (t) and e 2 2 = a −1 (t) f −1 (r ) where a(t) is the scale factor. Anisotropic parameter σ + = 1 3 ∂ ∂t ln(e 1 1 /e 2 2 ) vanishes and Hubble parameter (10) can be written as for closed, flat and open FLRW models, respectively. Therefore, evolution/constraint equations are reduced tö By settingφ = ρ m = 0 and V (φ) = Λ we obtain vacuum cases (with or without cosmological constant Λ), which are de Sitter model (Λ > 0, k = 0), the model with Λ > 0, k = 1, the model with Λ > 0, k = −1, Milne model (Λ = 0, k = −1) and Minkowski spacetime (Λ = 0, k = 0), which is also static. The model with Λ > 0, k = 1 is past asymptotic to de Sitter model with negative H, and, is future asympotic to de Sitter model with positive H . The model with Λ > 0, k = −1 (and positive H ) is past asymptotic to Milne model and it is future asympotic to de Sitter model with positive H . In this paper we will study open FLRW model. By convenience we write the metric (28) for k = −1 as where S −1 (r ) = sinh(r ), dΩ 2 = dϑ 2 + sin 2 ϑ dζ 2 .

Averaging scalar field cosmologies
Given the differential equationẋ = f(t, x, ε) with f periodic in t. One approximation scheme which can be used to solve the full problem is solving the unperturbed probleṁ x = f(t, x, 0) by setting ε = 0 and then use the approximated unperturbed solution to formulate variational equations in standard form which can be averaged. The term averaging is related to approximation of initial value problems in ordinary differential equations which involves perturbations [161, chapter 11].

Simple example
For example, consider this simple equation with φ(0) andφ(0) given. The unperturbed problemφ + ω 2 φ = 0 admits solutionφ(t) = r 0 ω cos(ωt − Φ 0 ), φ(t) = r 0 sin(ωt − Φ 0 ), where r 0 and Φ 0 are constants depending on the initial conditions. Let be defined the amplitude-phase transformation [161, chapter 11]: such that Then, Eq. (29) becomeṡ From (32) it follows that r and Φ are slowly varying with time, and the system takes the formẏ = ε f (y). The idea is consider only nonzero average of right-hand-sides keeping r and Φ fixed and leaving out terms with average zero and ignoring slow-varying dependence of r and Φ on t through averaging process: Replacing r and Φ by their averaged approximationsr and Φ we obtain the systeṁ Solving (34) with initial conditionsr (0) = r 0 andΦ(0) = Φ 0 , we obtainφ = r 0 e −εωt sin(ωt−Φ 0 ), which is an accurate approximation of the exact solution

General class of systems with a time-dependent perturbation parameter
Generalizing (29), now let us consider the KG system The similarity between (29) and (35) suggests treating the latter as a perturbed harmonic oscillator as well, and applying averaging in an analogous way; taking into consideration that, in contrast to ε, H is time-dependent and itself is governed by evolution Eq. (36). If it is valid, then a surprising feature of such approach is the possibility of exploiting the fact that H is strictly decreasing and goes to zero, therefore, it can be promoted to a time-dependent perturbation parameter in (35); controlling the magnitude of the error between solutions of full and time-averaged problems. Hence, with strictly decreasing H the error should decrease as well. Therefore, it is possible to obtain information about large-time behaviour of more complicated full system via an analysis of simpler averaged system equations.
With this in mind, in [153] the long-term behavior of solutions of a general class of systems in standard form was studied; where H > 0 is strictly decreasing in t and lim t→∞ H (t) = 0. The following Theorem by [153] gives local-in-time asymptotics for system (37). Let the norm · denotes the standard discrete 1 -norm u := n i |u i | for u ∈ R n . Let also L ∞ x,t denotes the standard L ∞ space in both t and x variables with norm defined as f L ∞ x,t := sup x,t |f(x, t)|. .
is Lipschitz continuous and f [2] is continuous with respect to x for all t ≥ t * . Also, assume that f 1 and f [2] are T -periodic for some Therefore, Hubble parameter H can be used as a timedependent perturbation parameter.
In this paper we study systems which are not in the standard form (37) but can be expressed as a series with center in H = 0 according to the equation depending on a parameter ω which is a free frequency that can be tuned to make f 0 (x, t) = 0. Therefore, systems can be expressed in the standard form (37).

LRS Bianchi III
Firstly, we consider LRS Bianchi III metrics for the generalized harmonic potential (2) minimally coupled to matter with field equations: Defining we obtain the full systeṁ and the Raychauhuri equation iṡ where the deceleration parameter is given by (46) can be symbolically written as a Taylor series of the form (38).
Notice that the term in expression (38) is eliminated imposing the condition bμ 3 + 2 f μ 2 − f ω 2 = 0, which defines an angular frequency ω ∈ R. Then, order zero terms in the series expansion around H = 0 are eliminated assuming ω 2 > 2μ 2 and setting f = bμ 3 ω 2 −2μ 2 , which is equivalent to tune ω. Hence, we obtain: where T andf(y) given by time-averaging (33) we obtain: Proceeding in analogous way as in references [128,137] we implement a local nonlinear transformation: Taking time derivative in both sides of (57) with respect to t we obtaiṅ where is the Jacobian matrix of g(H, x 0 , t) for the vector x 0 . The function g(H, x 0 , t) is conveniently chosen. By substituting (49) and (57) in (59) we obtain where Then we obtaiṅ Using Eq. (50), we haveḢ = O(H 2 ). Hence, The strategy is to use Eq. (63) for choosing conveniently ∂ ∂t g(0, x 0 , t) to prove thaṫ is unknown at this stage. By construction we neglect dependence of ∂g i /∂t and g i on H , i.e., assume g = g(x 0 , t) because dependence of H is dropped out along with higher order terms Eq. (63). 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 (65) is almost periodic of period L = 2π ω for large times. Then, implementing the average process (33) on right hand side of (65), where slow-varying dependence of quantities Ω 0 , Σ 0 , Ω k0 , Φ 0 andΩ,Σ,Ω k ,Φ on t are ignored through averaging process, we obtain Defining the average (66) is zero so that g(x 0 , t) is bounded. Finally, Eq. (64) transforms tȯ and Eq. (65) is simplified to Theorem 2 establish the existence of the vector (58). Secondly, we consider FLRW metric (24) with k = −1 for generalized harmonic potential (2) minimally coupled to matter with field equations (27) with the substitution k = −1. Defining the full system is deduced from (46) by setting Σ = 0. Then, with deceleration parameter Setting f = bμ 3 ω 2 −2μ 2 > 0, we obtain: where T and using the time averaging (33) we obtain the time-averaged system: Theorem 2 applies to Bianchi III, and the invariant set Σ = 0 corresponds to negatively curved FLRW models.

Qualitative analysis of averaged systems
Theorem 2 proved in Appendix A implies that Ω, Σ, Ω k and Φ evolve according to time-averaged system (52), (53), (54), (55) as H → 0. Hence, the full equations of time-dependent system (46) and their corresponding time-averaged versions have the same late-time dynamics as H → 0. Therefore, the simplest time-averaged system determines the future asymptotic of full system. In particular, depending on values of barotropic index γ , the generic late-time attractors of physical interests are found. With this approach, the oscillations entering full system through KG equation can be controlled and smoothed out as Hubble factor H -acting as a timedependent perturbation parameter -tends monotonically to zero. These results are supported by numerical simulations given in Appendix C.
Evaluating Raychaudhuri equation (46e) at Q and solving it we obtain Hence, Gauss equation (16) and evolution equation (17) becomė anḋ Then, by integration, That is, line element (22) becomes Therefore, the corresponding solution can be expressed in form of non-flat LRS Kasner ( Hence, Gauss equation (16) and evolution equation (17) becomė Hence, That is, line element (22) becomes Setting ψ(t) = tω − Φ(t) and evaluating Raychaudhuri equation (46e) at F we obtaiṅ Therefore, for large t. In average, Φ is constant, setting for simplicity Φ = 0 and integrating we obtain where H 0 is the current value of H (t). Finally, H (t) ∼ 2 3t for large t.
(iii) It is nonhyperbolic for γ = 2 3 or γ = 1 or γ = 2. The corresponding solution is a matter-curvature scaling solution withΩ m = 3(1−γ ). We obtain same expressions in (104) for , H and for ds 2 we obtain the expression: In Fig. 2a some orbits in the phase space of the guiding system (78a), (78b), (78c) for γ = 0 corresponding to cosmological constant are presented. The attractor is F 0 where scalar field mimics a cosmological constant. The equilibrium point D is a saddle.
In Fig. 2b some orbits of the phase space of the guiding system (78a), (78b), (78c) for γ = 2 3 are presented. The point F 0 coincides with MC; it is asymptotically stable as proved in Appendix B.1 by means of Center Manifold theory. The equilibrium point D is a saddle.
In Fig. 3a, b some orbits in the phase space of the guiding system (78a), (78b), (78c) for γ = 0.8 and γ = 0.9 are presented. In both cases MC is a stable node and D is a saddle.
According to the center manifold analysis in Appendix B.2 and supported by Fig. 10 for γ = 1, it is shown that D is unstable (saddle type) for 1 8  if we restrict the analysis to 1 8 (4Σ + 4Ω k − 5) < 0 which is the physical region of the phase space, D is asymptotically stable and behaves as a local attractor.
In Fig. 5a some orbits in the phase space of guiding system (78a), (78b), (78c) for γ = 4 3 corresponding to radiation are presented. In Fig. 5b some orbits in the phase space of guiding system (78a), (78b), (78c) for γ = 2 which corresponds to stiff matter are presented. In both figures, the attractor on the invariant setΩ k = 0 is F. For γ > 1, D is locally asymptotically stable according to the center mani- In Table 1 exact solutions associated with the equilibrium points of reduced averaged system (78a), (78b) and (78c) are summarized. A(t) and B(t) denote scale factors of the metric (22) where c 1 , c 2 , a 0 ∈ R + . Table 1 Exact solutions associated with equilibrium points of reduced averaged system (78a), (78b) and (78c). A(t) and B(t) denote scale factors of metric (22).
Bianchi III form of flat spacetime Matter-curvature scaling solution

Late-time behaviour
The results from the linear stability analysis combined with Theorem 2 (for Σ = 0) lead to:

Theorem 4
The late time attractors of the full system (71) and the averaged system (107) are: (i) The matter dominated FLRW universe F 0 with line element (112) for 0 < γ ≤ 2 3 . F 0 represents a quintessence fluid or a zero-acceleration model for γ = 2 3 . In the limit γ = 0 we have de Sitter solution.

Conclusions
In this paper we have used asymptotic methods and averaging theory to explore the solution's space of scalar field cosmologies with generalized harmonic potential (2) in vacuum or minimally coupled to matter. We have studied systems that can be expressed in the standard form (37), where H is the Hubble parameter and is positive strictly decreasing in t and lim t→∞ H (t) = 0. Defining Z = 1 − Σ 2 − Ω 2 − Ω k which is monotonic as H → 0, due to there is a continuous function α such thatŻ = −H Zα + O H 2 it was proved that sign of 1−Σ 2 −Ω 2 −Ω k is invariant as H → 0. Then, from Eq. (50) it was proved that H is a monotonic decreasing function of t if 0 < Ω 2 + Σ 2 + Ω k < 1. This implies lim t→∞ H (t) = 0 based on the invariance of initial surface for large t. Hence, from Theorem 2 was deduced that Ω, Σ, Ω k , and Φ evolve according to the time-averaged system (52), (53), (54), (55) as H → 0. Then, the stability of a periodic solution can be established as it matches exactly the stability of stationary solution of time-averaged system. We have given a rigorous demonstration of Theorem 2 in Appendix A based on the construction of a smooth local near-identity nonlinear transformation, well-defined as H tends to zero. We have used properties of the sup norm, and the theorem of the mean values for a vector functionf : R 3 −→ R 3 . We have explained preliminaries of the method in Sect. 4.3.
In particular, in LRS Bianchi III late time attractors of full system (46) and averaged system (78) for Bianchi III line element are: (i) The matter dominated FLRW universe F 0 with line element (105) if 0 < γ ≤ 2 3 . F 0 represents a quintessence fluid or a zero-acceleration model for γ = 2 3 . In the limit γ = 0 we have de Sitter solution.
For FLRW metric with k = −1, late time attractors of full system (71) and averaged system are: (i) The matter dominated FLRW universe F 0 with line element (112) for 0 < γ ≤ 2 3 . F 0 represents a quintessence fluid or a zero-acceleration model for γ = 2 3 . In the limit γ = 0 we have de Sitter solution.
Summarizing, in LRS Bianchi III late-time attractors are: a matter dominated flat FLRW universe if 0 ≤ γ ≤ 2 3 (mimicking de Sitter, quintessence or zero acceleration solutions), a matter-curvature scaling solution if 2 3 < γ < 1 and Bianchi III flat spacetime for 1 ≤ γ ≤ 2. For FLRW metric with k = −1 late time attractors are: the matter dominated FLRW universe if 0 ≤ γ ≤ 2 3 (mimicking de Sitter, quintessence or zero acceleration solutions) and Milne solution if 2 3 < γ < 2. In all metrics, matter dominated flat FLRW universe represents quintessence fluid if 0 < γ < 2 3 . Continuing the program "Averaging generalized scalar field cosmologies", the cases: (II) Bianchi I and flat FLRW model and (III) KS and closed FLRW are studied in two companion papers [164,165], respectively. In [164] using the same approach used here we found the oscillations entering the full system through KG equation can be controlled and smoothed out when the Hubble factor H is used as a timedependent parameter since it tends monotonically to zero and preserves its sign during the evolution. However, in case (III) this approach is not valid given that H is not necessarily monotonically decreasing to zero, and it can change its sign. Therefore, we have developed an alternative procedure in [165]. For LRS Bianchi I and flat FLRW metrics as well as for LRS Bianchi III and open FLRW, we use Taylor expansion with respect to H near H = 0 such that the resulting system can be expressed in standard form (37) after selecting a convenient angular frequency ω in the transformation (31). Next, we have taken the time-averaged of previous system, obtaining a system that can be easily studied using dynamical systems tools. Using the last approach, we have formulated Theorems 3 and 4 about late-time behavior of our model, whose proofs are based on Theorem 2 center manifold calculations and linear stability analysis.
As in paper [152], our analytical results were strongly supported by numerics in Appendix C as well. We showed that asymptotic methods and averaging theory are powerful tools to investigate scalar field cosmologies with generalized harmonic potential, which have evident advantages, e.g., it is not needed to analyze the full dynamics to determine the stability of the full oscillation, but only the late-time behavior of the time-averaged (simpler) system has to be analyzed. Interestingly, for LRS Bianchi III and open FLRW model, when matter fluid corresponds to a cosmological constant, H tends asymptotically to constant values depending on the initial conditions which is consistent to de Sitter expansion (see Figs. 13a, 21a). In addition, for open FLRW and any γ < 2 3 and Ω k > 0, Ω k → 0. On the other hand, when γ > 2 3 and Ω k > 0 the universe becomes curvature dominated asymptotically (Ω k → 1). In Appendix C, evidences that the main theorem of Sect. 4 is valid for LRS Bianchi III and for FLRW metrics with negative curvature are presented.

Acknowledgements The research of Genly Leon and Esteban
González was funded by Agencia Nacional de Investigación y Desarrollo-ANID through the program FONDECYT Iniciación grant no. 11180126. Alfredo D. Millano was supported by Agencia Nacional de Investigación y Desarrollo -ANID -Subdirección de Capital Humano/ Doctorado Nacional/año 2020-folio 21200837. Claudio Michea was supported by Agencia Nacional de Investigación y Desarrollo -ANID -Subdirección de Capital Humano/Doctorado Nacional/año 2021-folio 21211604. Vicerrectoría de Investigación y Desarrollo Tecnológico at Universidad Católica del Norte is acknowledged by financial support. Genly Leon thanks Bertha Cuadros-Melgar for her useful comments. Ellen de los Milagros Fernández Flores is acknowledged for proofreading this manuscript and improving the English. We thank anonymous referee for his/her comments which have helped us improve our work.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: This is a theoretical research project.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4. Then, almost everywhere for t in 0 ≤ t ≤ T . In particular, if where Df denotes the Jacobian matrix of f and the integral of a matrix is understood as componentwise.

Proof of Theorem 2 Defining
From Eq. (50) it follows that H is a monotonic decreasing function of t if 0 < Ω 2 + Σ 2 + Ω k < 1. These allow to define recursively bootstrapping sequences , such that lim n→∞ H n = 0 y lim n→∞ t n = ∞.
Keeping the terms of second order in H , system (68) becomeṡ (A.8d) can be written as a 3-dimensional system: plus Eq. (A.8d), where the vector functionf is explicitly given (last row corresponding toΔΦ 0 was omitted) by: It is a vector function with polynomial components in variables (y 1 , y 2 , y 3 ). Therefore, it is continuously differentiable in all its components. Let be Δx Using same initial conditions for x 0 andx we obtain by integration: Using Lemma 6 we havē where Df denotes the Jacobian matrix off and the integral of a matrix is understood as componentwise. Omitting the dependence on s we calculate the components of which are Taking sup norm |Δx 0 | = max |Ω 0 −Ω|, |Σ 0 −Σ|, |Ω k0 −Ω k | and the sup norm of a matrix |A| defined by max{|a i j | : i = 1, 2, 3, j = 1, 2, 3}, where a i j are given in (A.11) we have By continuity of polynomials a i j Ω 0 , Σ 0 , Ω k0 , Φ 0 ,Ω,Σ, Ω k ,Φ given in (A.11) and by continuity of functions ] the following finite constants are found: |A(t)|, , and such that for all t ∈ [t n , t n+1 ]: Using Gronwall's Lemma 5, we have for t ∈ [t n , t n+1 ]: Furthermore, from eq. (A.8d) we have Finally, taking limit as n → ∞, we obtain H n → 0. Then, as H n → 0, functions Ω 0 , Σ 0 , Ω k0 , Φ 0 andΩ,Σ,Ω k ,Φ have the same limit as τ → ∞.
for some δ > 0. Therefore, we can use Taylor series to define The following quasilinear differential equations (B.42) Fig. 12 Some solutions of the full system (46) (blue) and time-averaged system (78) (orange), corresponding to LRS Bianchi III metric, when for which the results of [152] are recovered. We have used initial data sets presented in the Table 3 (a) (b) must be solved for a i , b i up to order six to approximate the center manifold. We notice that the even terms are zero and . Hence, we obtain Letting γ = 4 3 and defining new variables we obtain new equations  (78) (orange) for the LRS Bianchi III metric when γ = 0, in the projection Ω k = 0. We have used for both systems initial data sets presented in the Table 3 (a) (b)  (46) (blue) and time-averaged system (78) (orange) for the LRS Bianchi III metric when γ = 0, in the projection Σ = 0. We have used for both systems initial data sets presented in the Table 3 (a) Letting γ = 3 2 and defining new variables we obtain new equations where 15 Some solutions of the full system (46) (blue) and time-averaged system (78) (orange) for the LRS Bianchi III metric when γ = 1, in the projection Ω k = 0. We have used for both systems initial data sets presented in the Table 3 (a) (b)  17 Some solutions of the full system (46) (blue) and time-averaged system (78) (orange) for the LRS Bianchi III metric when γ = 4/3, in the projection Ω k = 0. We have used for both systems initial data sets presented in the Table 3 (a)

(b)
Therefore, the origin is locally asymptotically stable as shown in Fig. 9.
To analyze the stability of D : (Ω,Σ,Ω k ) = (0, 1 2 , 3 4 ) we define new variables to obtain new equations  Table 3 (a) We propose the Taylor expansion h(y, z) = a 1 y 2 + a 2 yz + a 3 z 2 where O(4) denotes terms of fourth order in the vector norm. Therefore, Eq. (B.78) can be expressed as (B.80) Equating the terms of the same power in y, z we have a solution a 1 = − 1 6 , a 2 = 0, a 3 = 0, b 1 arbitrary, b 2 = Fig. 21 Some solutions of the full system (71) (blue) and time-averaged system (107) (orange) for the FLRW metric with negative curvature (k = −1) when γ = 0. We have used for both systems initial data sets presented in the Table 4 (a) The dynamics at the center manifold is given by where is shown that the origin is unstable (saddle type) for y = 0 is presented. However, if we restrict the analysis to y < 0, D is asymptotically stable and behaves as a local attractor.
To analyze the stability of an arbitrary point D * : (Ω,Σ,Ω k ) = (Ω * , 1 2 , 3 4 ) withΩ * = 0 we define new variables 22 Some solutions of the full system (71) (blue) and time-averaged system (107) (orange) for the FLRW metric with negative curvature (k = −1) when γ = 1. We have used for both systems initial data sets presented in the Table 4 (a) to obtain new equations

Appendix C: Numerical simulation
In this section we present numerical evidence that supports the main Theorem of Sect. 4 by solving numerically full and time-averaged systems obtained for each metric, namely LRS Bianchi III and open FLRW. For this purpose an algorithm in the programming language Python was implemented. The systems of differential equations were solved using the solve_ivp code provided by the SciPy open-source Python-based ecosystem. The integration method used was Radau that is an implicit Runge-Kutta method of the Radau IIa family of order 5 with a 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 the range of integration −40 ≤ τ ≤ 10 for original systems and −40 ≤ τ ≤ 100 for averaged systems. All of them partitioned in 10000 data points. Furthermore, each full and time-averaged systems were solved considering only one matter component, these are cosmological constant (γ = 0), non relativistic matter or dust (γ = 1), radiation (γ = 4/3) and stiff fluid (γ = 2). Thereby the vacuum solutions corresponds to those where Ω = Ω m ≡ 0 and the solutions without matter component corresponds to Ω m ≡ 0. Finally we have considered fixed constants μ = √ 2/2, b = √ 2/5 and ω = √ 2, that leads to a value of f = bμ 3 ω 2 −2μ 2 = 1/10, that fulfills condition f ≥ 0. With these values a generalized harmonic potential of the form (1 − cos(10φ)) (C.103) is obtained.
As initial conditions we use seven data set presented in Table 3 as initial conditions for a better comparison of both systems. For data set vii current values of Ω m (0) = 0.315 and Ω k (0) = 0.001 according to [213] were considered.
It is important to mention that the first six initial conditions correspond to initial conditions presented in Table 2 of [152] and the additional data set vii is obtained considering current values of Ω m (0) = 0.315 and Ω k (0) = 0.001 according to [213]. Even more, the model presented in [152] is contained in our model when γ = b = f = 0, ω 2 = 2, μ = 1, Ω m = 0 and Ω k = 1 − Σ 2 − Ω 2 (then γ does not appear in the model presented in [152]) with the identification Ω 2 → Ω. As can be seen in Fig. 12a, b where some solutions of the full system (46) and time-averaged system (78) are presented, showing that our results are in complete agreement with results presented in [152] for the limiting case.
Seven data set presented in Table 4 were used as initial conditions. For data set vii current values of Ω m (0) = 0.315 and Ω k (0) = 0.001 according to [213] were considered.
In Figs. 21, 22, 23 and 24 projections of some solutions of the full system (71) and time-averaged system (107) for the FLRW metric with negative curvature (k = −1) in the (Ω k , H, Ω 2 ) space with their respective projection when H = 0 and considering in both systems the same initial data sets from Table 4 are presented. Figure 21 show solutions for matter fluid corresponding to cosmological constant (γ = 0). Figure 22 show solutions for matter fluid corresponding to dust (γ = 1). Figure 23 show solutions for matter fluid corresponding to radiation (γ = 4/3). Figure 24 show solutions for matter fluid corresponding to stiff fluid (γ = 2). It is interesting to note that in the FLRW with negative curvature case, when the matter fluid corresponds to a cosmological constant, H tends asymptotically to constant values depending on the initial conditions which is consistent to de Sitter expansion. In addition, for any γ < 2 3 and Ω k > 0, Ω k → 0. On the other hand, when γ > 2 3 and Ω k > 0 the universe becomes curvature dominated asymptotically (Ω k → 1). These figures are evidence that the main Theorem presented in Sect. 4 is fulfilled for FLRW metric with negative curvature.