Averaging Generalized Scalar Field Cosmologies II: Locally Rotationally Symmetric Bianchi I and flat Friedmann-Lema\^itre-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 $\gamma$ for the Locally Rotationally Symmetric (LRS) Bianchi I and flat Friedmann-Lema\^itre-Robertson-Walker (FLRW) metrics 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, the simplest time-averaged system determines the future asymptotic behavior. Depending on the values of $\gamma$, the late-time attractors of physical interests are flat quintessence dominated FLRW universe and Einstein-de Sitter solution. With this approach, the oscillations entering the system through the Klein-Gordon (KG) equation can be controlled and smoothed out as the Hubble parameter $H$ - acting as time-dependent perturbation parameter - tends monotonically to zero. Numerical simulations are presented as evidence of such behavior.


Introduction
Mathematical methods have been widely used in cosmology. For example, in reference [1] the method of a genly.leon@ucn.cl b sebastian.cuellar01@ucn.cl c esteban.gonzalezb@usach.cl d samuel.lepe@pucv.cl e claudio.ramirez@ce.ucn.cl f alfredo.millano@alumnos.ucn.cl Lie symmetries was applied to 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 determined and classified. In reference [2] a modelindependent criterion based on first integrals of motion was used; and in reference [3] dynamical symmetries of the field equations were used to classify Dark Energy (DE) models in the context of scalar field (quintessence or phantom) FLRW cosmologies. Using Noether symmetries in [2] the system was simplified and its integrability was determined. For the exponential potential as well as some types of hyperbolic potentials, extra Noether symmetries apart of the conservation law were found; suggesting that these potentials should be preferred along the hierarchy of scalar field potentials. In [3] under the requirement that field equations admit dynamical symmetries resulted in two potentials, one of them is the well known Unified Dark Matter (UDM) potential and another hyperbolic model. In reference [4] a mathematical approach to reconstruct the EoS and the inflationary potential of the inflaton field from observed spectral indices for the density perturbations and the tensor-to-scalar ratio (based on the constraints system) was implemented. In reference [5] an algorithm to generate new solutions of the scalar field equations in FLRW universes was used. Solutions for pure scalar fields with various potentials in absence and in presence of spatial curvature and other perfect fluids were obtained. A series of generalizations of Chaplygin gas and bulk viscous cosmological solutions for inflationary universes were found. In reference [6] the f (T ) cosmological scenario was studied. In particular, analytical solutions for isotropic, homogeneous universe containing dust fluid and radiation, and for an empty anisotropic Bianchi I universe were found. The method of movable singularities of differential equations was used. For the isotropic universe, the solutions are expressed in terms of Laurent expansion, while for anisotropic universe a family of exact Kasner-like solutions in vacuum is found. In reference [7] the symmetry classification of the KG equation in Bianchi I spacetime was performed. A geometric method which relates the Lie symmetries of the KG equation with the conformal algebra of the underlying geometry was applied. Furthermore, by means of Lie symmetries that follow from the conformal algebra (which are also Noether symmetries for the KG equation) all the potentials in which the KG equation admits Lie and Noether symmetries were determined. The Lie admitted symmetries are useful to determine the corresponding invariant solution of the KG equation for specific potentials. Additionally, the classification problem of Lie/Noether point symmetries of the wave equation in Bianchi I spacetime was solved, and invariant solutions of the wave equation were determined. In reference [8] a new method to classify Bianchi I spacetimes which admits Conformal Killing Vectors (CKV) was developed. The method is useful to study the conformal algebra of Kasner spacetime and other Bianchi type I matter solutions of GR.
Inspired in [16][17][18][19] we have started the "Averaging Generalized Scalar Field Cosmologies" program which consists in using asymptotic methods and averaging theory to obtain relevant information about the solution's space of scalar field cosmologies with generalized harmonic potential in presence of matter (with a barotropic EoS with barotropic index γ) minimally coupled to a scalar field. This research program has three steps according to the three cases of study: (I) Bianchi III and open FLRW model [20], (II) Bianchi I and flat FLRW model (the present case) and (III) Kantowski-Sachs and closed FLRW [21]. In reference [19] relevant results for the aforementioned program were presented. In particular, interacting scalar field cosmologies with generalized harmonic potentials for flat and negatively curved FLRW, and for Bianchi I metrics were studied. Using asymptotic and averaging methods stability conditions for several solutions of interest as H → 0 were obtained. This analysis suggests that the asymptotic behavior of the time-averaged model is independent of the coupling function and the geometry. Following analogous procedures in references [20] and [21] the cases (I) and (III) of the program were studied.
This paper is devoted to case (II). It 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 [17] where in particular, Section 4.1 is devoted to LRS Bianchi I model and Section 4.2 is devoted to flat FLRW metric. In Section 5 we study the resulting timeaveraged systems where, in particular, Section 5.1 is devoted to LRS Bianchi I models and Section 5.2 is devoted the flat FLRW metric. Finally, in Section 6 our main results are discussed. In Appendix A the proof of our main theorem is given and in Appendix B numerical evidence supporting the results of Section 4 is presented.

The generalized harmonic potential
Scalar fields are relevant in the physical description of the universe, particularly, in inflationary scenario [49][50][51][52][53]. For example, chaotic inflation is a model of cosmic inflation in which the potential term takes the form of [50][51][52][53]. In this research we consider the generalized harmonic potential which incorporates cosine-like corrections with µ 3 b f 1. Introducing a new parameter ω through the equation (1) can be reexpressed as The applicability of this re-parametrization will be discussed at the end of section 2.3. Potential (2) 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 local maximums or local minimums depending on The function V has no upper bound but it has a lower bound equal to zero.
The asymptotic features of potential (2) are the following. Near the global minimum φ = 0, we have V (φ ) ∼ That is, ω 2 can be related to the mass of the scalar field near its global minimum. As φ → ±∞ the cosine-correction is bounded, then V (φ ) ∼ µ 2 φ 2 +O (1) as φ → ±∞. This makes it suitable to describe oscillatory behavior in cosmology. Potential (1) or (2) is related but not equal to the monodromy potential of [54] used in the context of loop-quantum gravity, which is a particular case of the general monodromy potential [55]. In references [17][18][19] it was proved that the potential of [54,55] for 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 [17,18] the following potential was studied that is obtained by setting µ = √ 2 2 and bµ = 2 in eq. (1). On the other hand, setting µ = √ 2 2 and ω = √ 2, we have The potentials (3) and (4) provide non-negative local minimums which can be related to a late-time accelerated universe. The generalized harmonic potentials (2), (3) and (4) belong to the class of potentials studied by [23]. Additionally, potentials like V (φ ) = are of interest in the context of axion models [56]. In [57] axionic dark matter with modified periodic potential for the pseudoscalar field has been studied in the framework of the axionic extension of the Einsteinaether theory. This periodic potential has minima at φ = nΦ * , n ∈ Z, whereas maxima are found when n → m + 1 2 . Near the minimum, i.e., φ = nΦ * + ψ with |ψ| a small value, V → where m A the axion rests mass.

Simple example
Given the ordinary differential equationẋ = f(x,t, ε) with ε ≥ 0 and f periodic in t. One approximation scheme which can be used to solve the full problem is the resolution of 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. The term averaging is related to the approximation of initial value problems involving perturbations (chapter 11, [15]). For example, consider the initial value problem: with φ (0) andφ (0) prescribed. The unperturbed problemφ + ω 2 φ = 0 admits the solutionφ (t) = r 0 ω cos(ωt − Φ 0 ), φ (t) = r 0 sin(ωt − Φ 0 ), where r 0 and Φ 0 are constants depending on initial conditions. Let be defined the amplitude-phase transformation ( [15], chapter 11): such that Then, eq. (5) becomes, From (8) 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 the righthand-sides 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 = −εωr,˙Φ = 0.

General class of systems with a time-dependent perturbation parameter
Let us consider for example the KG system The similarity between (5) and (11) suggests to treat the latter as a perturbed harmonic oscillator as well, and to apply averaging in an analogous way. However, care has to be taken because in contrast to ε, H is timedependent and itself is governed by the evolution equation (12). If it is valid, then a surprising feature of such approach is the possibility of exploiting the fact that it is strictly decreasing and goes to zero by promoting Hubble parameter H to a time-dependent perturbation parameter in (11) controlling the magnitude of the error between solutions of the full and time-averaged problems. Hence, with strictly decreasing H the error should decrease as well. Therefore, it is possible to obtain the information about the large-time behavior of the more complicated full system via an analysis of the simpler averaged system equations by means of dynamical systems techniques . With this in mind, in [43] the long-term behavior of solutions of a general class of systems in standard form was studied: where H is positive strictly decreasing in t and lim t→∞ H(t) = 0. In this paper we study systems which are not in the standard form (13) 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 (13). In particular, assuming ω 2 > 2µ 2 and setting f = bµ 3 (which is equivalent to tune the angular frequency ω) the undesired terms evolving as ∝ H 0 are eliminated in the series expansion around H = 0.

The model
It is well-known that there is an interesting hierarchy in Bianchi models [86][87][88][89]. 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 LRS Bianchi type VIII model as well. Additionally, LRS Bianchi type VIII can be viewed as an invariant boundary of LRS Bianchi type IX models [90][91][92][93][94][95]. Bianchi spacetimes contain many important cosmological models that have been used to study anisotropies of primordial universe and its evolution towards the observed isotropy of the present epoch [96][97][98][99][100]. The list includes FLRW model in the limit of the isotropization.
In GR the Hubble parameter is always monotonic for Bianchi I and anisotropies decay for H > 0. Therefore, isotropization occurs [101,102]. The exact solutions of field equations have been found in some particular Bianchi spacetimes for an exponential potential [103][104][105]. These exact solutions lead to isotropic homogeneous spacetimes as it was found in references [106,107]. An anisotropic solution of special interest is Kasner spacetime [108][109][110][111][112][113][114][115], essential for the description of BKL singularity [116]. The action integral of interest given by is expressed in a system of units in which 8πG = c = = 1. In eq. (15) R is the scalar curvature of the spacetime, L m is the Lagrangian density of matter, φ is the scalar field, ∇ α is the covariant derivative and V (φ ) is the scalar field potential defined by (1).

LRS Bianchi III, Bianchi I and Kantowski-Sachs models
Considering that the metric element for LRS Bianchi III, Bianchi I and Kantowski-Sachs models can be written as [117] 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 [118]: e 0 = ∂ t , e 1 = e 1 1 ∂ r , e 2 = e 2 2 ∂ ϑ , e 3 = e 3 3 ∂ ζ . Comparing with reference [117] we have settled the parameters a = f = 0 and e 1 1 (t) = D 2 (t) −1 , e 2 2 (t) = D 2 (t) −1 and we have used the identifications (ϑ , ζ ) = (y, z). The line elements for spatially homogeneous self-similar LRS models have been given by Wu in [119]. We concentrate only in the spatially homogeneous but anisotropic class with the exception of spatially homogeneous LRS Bianchi V, that is: LRS Bianchi III (k = −1), Bianchi I (k = 0) and Kantowski-Sachs (k = +1) [120]. It is useful to define a representative length (t) along worldlines of u = ∂ t for describing the volume expansion (contraction) behavior of the congruence completely by [121] (t) where dots denote derivatives with respect to time t, H(t) is the Hubble parameter in terms of (t) its time derivative. The anisotropic parameter σ + (t) is defined by The variation of (15) for the 1-parameter family of metrics (19) leads to [120]: where for the matter component we use barotropic EoS p m = (γ − 1)ρ m with p m the pressure of the fluid, ρ m is the energy density and the barotropic index is a constant γ which satisfies 0 ≤ γ ≤ 2.
The Gauss curvature of the spatial 2-space and 3curvature scalar are [118] Furthermore, the evolution equation of the Gauss curvature of the spatial 2-space iṡ while the evolution for e 1 1 is given by [118] From eqs. (23) and (24) the shear equatioṅ is obtained. Eqs. (22), (23), (24) and (28) give the Raychaudhuri equatioṅ Finally, the matter and KG equations arė In this paper we will focus our study in LRS Bianchi I model. Therefore, using eq. (17) the metric (19) reduces to where the functions A(t) and B(t) are interpreted as the scale factors: A(t) = e 1 1 (t) −1 and B 2 (t) = K(t) −1 .

FLRW models
The general line element for spherically symmetric models can be written as [118] Spatially homogeneous spherically symmetric models, that are not Kantowski-Sachs, are the FLRW models where the metric can be written as with f (r) = sin r, r, sinh r, for closed, flat and open FLRW models, respectively. In comparison with metric (33), the frame coefficients are given by e 1 is the scale factor. The anisotropic parameter σ + = for closed, flat and open FLRW, respectively. Therefore, evolution/constraint equations reduce tö

Averaging scalar field cosmologies
As in reference [42] we construct a time-averaged version of the original system and prove that it shares the same late-time dynamics of the original system.

Bianchi I metric
In this section averaging methods are applied for Bianchi I metrics for the generalized harmonic potential (1) minimally coupled to matter.
Setting k = 0 in Eqs. (22), (28) and (29) we obtaiṅ Using the characteristic length scale along worldlines of the 4-velocity field such that H =˙ , defining 0 the current value of such that (t) and denoting by convention t = 0 the current time, Using the definition (20) and integrating (39) we obtain σ + = σ 0 3 0 / 3 , where σ 0 is an integration constant, which is the value of σ + when = 0 . The term G 0 ( ) = σ 2 0 6 0 / 6 , which corresponds to anisotropies in Bianchi I metric, does not correspond to a fluid component in the model. However, it can be interpreted as a stiff-matter fluid for flat FLRW metric with scale factor a(t) = (t). The term σ 2 + dilutes very fast with expansion, isotropizing if H > 0. The evolution equation for matter and the KG equation do not depend on k. Therefore, the field equations are deduced: Now, we define Hubble normalized variables along with r and Φ which are defined in (7) and σ + = σ 0 3 0 / 3 is obtained by integrating (39). Then, we obtain the systeṁ where the deceleration parameter q is given by (44) can be symbolically written in form (14). Notice that using the condi- in the eq. (14) becomes trivial. Hence, we obtain: where (48) is Raychaudhuri equation and Replacingẋ = Hf(x,t) with f(x,t) as defined in (49) bẏ y = Hf(y) with y = Ω ,Σ ,Φ T andf as defined by (9), we obtain the averaged system:Ω Proceeding in analogous way as in references [24,33] but for 3 dimensional systems instead of a 1dimensional one, we implement a local nonlinear transformation which in vector form can be written as Taking derivative of (55) with respect to t, we obtaiṅ where is the Jacobian matrix of g(H, x 0 ,t) for the vector x 0 . The function is conveniently chosen. By substituting (47) and (55) in (57) we obtain where Then we obtaiṅ Using eq. (48), we haveḢ = O(H 2 ). Hence, The strategy is to use eq. (61) for choosing conveniently ∂ ∂t g(0, x 0 ,t) in order 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 in eq. (61). Next, we solve the differential equation for g(x 0 ,t): where we have considered x 0 and t as independent variables. The right hand side of (63) is almost periodic of period L = 2π ω for large times. Then, implementing the average process (9) on right hand side of (63), where the slowvarying dependence of quantities x 0 = (Ω 0 , Σ 0 , Φ 0 ) T andx = (Ω ,Σ ,Φ) T on t are ignored through the averaging process, we obtain Defining and eq. (63) is simplified to Theorem 1 establishes the existence of the vector (56).
Proof. The proof is given in Appendix A. Theorem 1 implies that Ω , Σ and Φ evolve according to the averaged equations (50), (51), (52) as H → 0 because (54) is a formal near-identity (this means that D x 0 ψ| H=0 = I) nonlinear change of coordinates, the first order solutions Ω 0 , Σ 0 , Φ 0 and averaged solutions Ω ,Σ ,Φ have the same limit when t → ∞ by monotony and non-negativity of H and the limit H → 0.

Flat FLRW metric.
In this case the field equations are obtained from (37) by setting k = 0. We obtainΩ ,Φ andḢ by substituting Σ = 0 in (44) and (45). Finally, we obtain the Taylor expansion: Replacingẋ = Hf(x,t) and f(x,t) as defined by (69) withẏ = Hf (y) where y = Ω ,Φ T with the time averaging (9), we obtain the following time-averaged system:Ω The time-averaged Raychaudhuri equation for flat FLRW metric is obtained by settingΣ = 0 in eq. (53). Theorem 1 applies for Bianchi I and the invariant set Σ = 0 corresponds to flat FLRW models. dt . Therefore, the analysis of the original system is reduced to study the corresponding averaged equations.

Bianchi I metric
The averaged system (50), (51), (52) and (53) is transformed to by defining the logarithmic time τ through dt dτ = 1/H. We investigate the 2D guiding system (72a)-(72b). The functionΩ m = 1 −Σ 2 −Ω 2 is interpreted as an averaged Hubble-normalized density parameter for the matter component. Therefore, imposing the energy conditionΩ m ≥ 0 the phase space is Before the discussion we introduce the following concept. 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 as per its stability by considering the signs of eigenvalues in the remaining directions (i.e., for a curve, in the remaining n − 1 directions) (see [67], pp. 36). The resulting 2D guiding system (72a)-(72b) has the following equilibrium points: 1. T : (Ω ,Σ ) = (0, −1) with eigenvalues 3 2 , 3(2 − γ) . i) It is a source for 0 ≤ γ < 2, ii) It is nonhyperbolic for γ = 2.

Evaluating Raychaudhuri equation (44d) at equilibrium point F we havė
Therefore, for large t. In average, Φ is a constant, setting Φ = 0 for simplicity and integrating we obtain where H 0 is the current value of H(t). Finally, H(t) ∼ 2 3t for large t. Eqs. (26) and (27) becomė with general solutioṅ Then, line element (32) becomes For large t, F can be associated with Einstein-de Sitter solution ( [89], Sec 9.1.1 (1)) with γ = 1).
In Figure 1 a phase plane for system (72a), (72b) for different choices of γ is presented. System (72a), (72b) when γ = 1 reduces to From the above system we obtain d dτ ln We have assumed that the orbit passes by (Ω ,Σ ) = (Ω 0 ,Σ 0 ) at time τ 0 = 0. These values are identified with the current epoch. Using equations (92) and the chain rule we have Then,   Due to (Ω ,Σ ) = (Ω 0 ,Σ 0 ) at the time τ 0 = 0 by solving (95) with the variables separation method it follows that Finally, the orbits of system (92) are given bȳ In table 1 exact solutions are associated with equilibrium points of the reduced averaged system (72a), (72b) are summarized, where A(t) and B(t) denote the scale factors of the metric (32) and c 1 , c 2 , 0 ∈ R + are integration constants.

Late-time behavior
Results from the linear stability analysis which are combined with Theorem 1 lead to:  Table 1 Exact solutions associated with equilibrium points of the reduced averaged system (72a), (72b). A(t) and B(t) denote the scale factors of metric (32) and c 1 , c 2 , 0 ∈ R + .

A(t) B(t) Solution
Flat matter dominated FLRW universe For γ = 1, F 0 and F are stable because they belong to stable normally hyperbolic line of equilibrium points. For γ = 2, F is asymptotically stable (see figures 1 (b) and 1 (d)).

Late-time behavior
The results from the linear stability analysis are combined with Theorem 1 (for Σ = 0) and lead to: Observe in figure 7(b) that as H → 0 the values of Ω lying on the orange line (solution of time-averaged system at fourth order) give an upper bound to the value Ω of the original system. Therefore, by controlling the error of averaged higher order system, one can also control the error in the original one.

Conclusions
This is the second paper of the "Averaging Generalized Scalar Field Cosmologies" program that was initiated in reference [20]. This program consists in using asymptotic methods and averaging theory to obtain relevant information about solution's space of scalar field cosmologies in presence of a matter fluid with EoS with barotropic index γ minimally coupled to a scalar field with generalized harmonic potential (1).
According to this research program, in paper I [20] was proved that late-time attractors in LRS Bianchi III model 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 if 1 ≤ γ ≤ 2. Late-time attractors in FLRW metric with k = −1 are: a 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, the matter dominated flat FLRW universe represents quintessence fluid if 0 < γ < 2 3 . For LRS Bianchi I and flat FLRW metrics as well as for LRS Bianchi III and open FLRW, we can use Taylor expansion with respect to H near H = 0. Hence, the resulting system can be expressed in standard form (13) after selecting a convenient angular frequency ω in the transformation (7). Next, we have taken the timeaveraged of previous system obtaining a system that can be easily studied using dynamical system's tools.
In particular, we have proved in Theorem 1 that latetime attractors of full and time-averaged systems are the same for some homogeneous metrics. Theorem 1 implies that Ω , Σ and Φ evolve according to timeaveraged equations (50), (51) and (52) as H → 0. Therefore, we can establish the stability of a periodic solution as it matches exactly the stability of a stationary solution of averaged equation.
We have given a rigorous demonstration of Theorem 1 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 2 −→ R 2 . We have explained preliminaries of the method of proof in Section 4.1. As in paper [42], our analytical results were strongly supported by numerics in Appendix B as well.
More specific, according to Theorem 1 for Bianchi I and flat FLRW metrics, Hubble parameter H plays the role of a time-dependent perturbation parameter controlling the magnitude of error between solutions of full and time-averaged systems. Therefore, analysis of system is reduced to study time-averaged equations.
In this regard, we have formulated theorems 2 and 3 concerning to the late-time behavior of our model. For LRS Bianchi I late-time attractors of full system (44) and averaged system (72) are: (i) The matter dominated FLRW Universe F 0 with line element (85) if 0 < γ < 1. F 0 represents a quintessence fluid or a zero-acceleration model for γ = 2 3 . In the limit γ = 0 we have a de Sitter solution.
(ii) The scalar field dominated solution F with line element (91) if 1 < γ ≤ 2. For large t the equilibrium point can be associated with Einstein-de Sitter solution.
For flat FLRW metric late-time attractors of full system (44) with Σ = 0 and averaged system are: (i) The matter dominated FLRW Universe F 0 with line element (100) if 0 ≤ γ < 1. F 0 represents a quintessence fluid or a zero-acceleration model for γ = 2 3 . In the limit γ = 0 we have a de Sitter solution.
(ii) The scalar field dominated solution F with line element (102) if 1 < γ ≤ 2. For large t the equilibrium point can be associated with Einsteinde Sitter solution.
It is interesting to note that for LRS Bianchi I and flat FLRW cases when matter fluid is a cosmological constant, H tends asymptotically to constant values depending on initial conditions which is consistent to de Sitter expansion (see figures 3(a) and 7(a)). For dust γ = 1 in flat FLRW metric, we have from qualitative analysis in Section 5.2 thatΩ (τ) tends to a constant and H(τ) tends to zero as τ → +∞. Observe in figure  7(b) that as H → 0 the values ofΩ lying on the orange line (solution of time-averaged system at fourth order) give an upper bound to values of Ω in the original system. Therefore, by controlling the error of averaged higher order system, one can also control the error of the original one. We have illustrated that asymptotic methods and averaging theory are powerful tools to investigate scalar field cosmologies with generalized harmonic potential. One evident advantage is that to determine stability of full oscillation it is not needed to analyze the full dynamics, but only the late-time behavior of timeaveraged (simpler) system has to be analyzed. In all metrics, the matter dominated flat FLRW universe represents quintessence fluid if 0 < γ < 2 3 .
where we set three integration functions C i (Ω 0 , Σ 0 , Φ 0 ), i = 1, 2, 3 to zero. The g i , i = 1, 2, 3 are continuously differentiable, such that their partial derivatives are bounded on t ∈ [t n ,t n+1 ]. The second order expansion around H = 0 of system (66) is written as: 7a) and (A.7b) are reduced to: where the vector functionf is explicitly given (the last row corresponding to eq. (A.7c) was omitted) by: It is a vector function with polynomial components in variables (y 1 , y 2 ). Therefore, it is continuously differentiable in all its components.
Using same initial conditions for x 0 andx we obtain by integration: The is finite by continuity ofΩ , Ω 0 ,Σ , Σ 0 , Φ 0 in the closed interval [t n ,t n+1 ]. Using Lemma 5 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 matrix elements of A as (A.14) Taking the sup norm |∆ : On the other hand where the sup norm of a matrix a b c d is defined by max{|a|, |c|, |b|, |d|}.
Step 2: For FLRW the second order expansion around H = 0 of system (66) is written as: Then, By continuity ofΩ , Ω 0 and Φ 0 in [t n ,t n+1 ] the following finite constants are found: such that the terms proportional to H 2 in eq. (A.16a) are bounded in absolute value by M 3 H 2 n and the terms proportional to H 2 in eq. (A.16b) are bounded in absolute value by M 4 H 2 n in the interval [t n ,t n+1 ]. Then, |∆ Ω 0 (s)|ds + M 3 H n due to t − t n ≤ t n+1 − t n = 1 H n .
Summarizing, according to Theorem 1 for Bianchi I and flat FLRW metrics, Hubble parameter H plays the role of a time-dependent perturbation parameter controlling the magnitude of the error between solutions of full and time-averaged systems. Therefore, the analysis is reduced to study time-averaged equations.

Appendix B: Numerical simulation
In this Section we present numerical evidence that supports the main results in Section 4 by solving numerically full and averaged systems for each metric, namely LRS Bianchi I and flat FLRW. For this purpose, an algorithm in the programming language Python was elaborated where systems of differential equations were solved using the solve ivp code that is provided by the SciPy open-source Pythonbased ecosystem. The integration method 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 τ in an integration range of −40 ≤ τ ≤ 10 for the original systems and in an integration range of −40 ≤ τ ≤ 100 for the averaged systems. All of them are partitioned in 10000 data points. Furthermore, each full and timeaveraged 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). Vacuum solutions correspond to Ω = Ω m ≡ 0 and solutions without matter fluid correspond to Ω m ≡ 0. Finally, we have considered these constants: µ = √ 2/2, b = √ 2/5 and ω = √ 2, that lead to the value of f = bµ 3 ω 2 −2µ 2 = 1/10 which fulfills the condition f ≥ 0. These values are translated into a generalized harmonic potential: 1. The full system (44). 2. The time-averaged system (72).
As initial conditions we use seven data sets which are presented in Table 3. Table 3 Seven initial data sets for the simulation of full system (44) and averaged system (72) for Bianchi I metric are displayed. All initial conditions satisfyΣ 2 (0) +Ω 2 (0) +Ω m (0)=1. In figures 3(a)-6(b) projections of some solutions in the (Σ , H, Ω 2 ) space of full system (44) and time-averaged system (72) along with their respective projection in the subspace H = 0 are presented. Both systems were integrated using the same initial data sets from Table  3.   3 Some solutions of the full system (44) (blue) and time-averaged system (72) (orange) for Bianchi I metric when γ = 0. We have used for both systems the initial data sets that are presented in Table 3.   4 Some solutions of the full system (44) (blue) and time-averaged system (72) (orange) for Bianchi I metric when γ = 1. The line denoted by red diamonds corresponds to the attracting line of equilibrium points. We have used for both systems the initial data sets that are presented in Table 3.   5 Some solutions of the full system (44) (blue) and time-averaged system (72) (orange) for Bianchi I metric when γ = 4/3. We have used for both systems the initial data sets that are presented in Table 3.   6 Some solutions of the full system (44) (blue) and time-averaged system (72) (orange) for Bianchi I metric when γ = 2. We have used for both systems the initial data sets that are presented in Table 3.   Table 4. which corresponds to dust (γ = 1). Figures 5(a)-5(b) show solutions for a matter fluid which corresponds to radiation (γ = 4/3). Figures 6(a)-6(b) show solutions for a matter fluid which corresponds to stiff fluid (γ = 2). These figures numerically support the main theorem that is presented in Section 4 for Bianchi I metric. As an interesting point, in this example, for γ = 0 when the matter fluid corresponds to a cosmological constant, H tends asymptotically to a constant, which is consistent to de Sitter expansion. 1. The full system given by (44) with Σ = 0. 2. The time-averaged system (98) if γ = 1 or the system (103) truncated at fourth order if γ = 1.
As initial conditions we use eight data sets that are presented in Table 4. Table 4 Eight initial data sets for the simulation of full system (44) (44) with Σ = 0 and time-averaged system (98) for flat FLRW metric (k = 0) are presented. Both systems were integrated using the same initial data sets from Table 4. Figure 7(a) shows solutions for a matter fluid which corresponds to cosmological constant (γ = 0). Figure 7(b) shows solutions for a matter fluid which corresponds to dust (γ = 1). Observe that as H → 0 the values ofΩ lying on the orange line (solution of time-averaged system at fourth order) give an upper bound to the value of Ω of the original system. Figure 7(c) shows solutions for a matter fluid which corresponds to radiation (γ = 4/3). Figure 7(d) shows solutions for a matter fluid which corresponds to stiff fluid (γ = 2). These figures support numerically the main theorem that is presented in Section 4 for flat FLRW metric. It is interesting to note that in the flat FLRW case when the matter fluid corresponds to a cosmological constant, H tends asymptotically to constant values depending on initial conditions which is consistent to de Sitter expansion (see figure 7(a)). With our approach, the oscillations which enter the system through the KG equation can be controlled and smoothed out as the Hubble parameter H tends monotonically to zero. This fact was analytically proved in Appendix A and we have presented numerical simulations as evidence of such behavior in Appendix B.