Time-averaging axion-like interacting scalar fields models

In this paper, we study a cosmological model inspired in the axionic matter with two canonical scalar fields $\phi_1$ and $\phi_2$ interacting through a term added to its potential. Introducing novel dynamical variables, and a dimensionless time variable, the resulting dynamical system is studied. The main difficulties arising in the standard dynamical systems approach, where expansion normalized dynamical variables are usually adopted, are due to the oscillations entering the nonlinear system through the Klein-Gordon (KG) equations. This motivates the analysis of the oscillations using methods from the theory of averaging nonlinear dynamical systems. We prove that time-dependent systems, and their corresponding time-averaged versions, have the same late-time dynamics. Then, we study the time-averaged system using standard techniques of dynamical systems. We present numerical simulations as evidence of such behavior.


Introduction
Between 1998 and 1999 the independent projects Highz Supernova Search Team and Supernova Cosmology Project showed results that the universe went through a late-time accelerated expansion phase [1,2]. This behavior has been confirmed by many independent observations [3][4][5], turning it into the most fascinating puzzle in modern cosmology. The most simple model to a snilch@yzu.edu.cn b genly.leon@ucn.cl c esteban.gonzalez@uac.cl d wang b@sjtu.edu.cn describe this behavior of the universe is the Λ CDM model, in which the universe is currently dominated by two non-interacting fluids called dark matter (DM) and dark energy (DE). The DE component is given by the cosmological constant (CC) Λ driving the present epoch of the accelerated expansion of the universe, and the DM component is assumed to have negligible pressure (cold DM or CDM) and it is responsible for the large scale structure formation in the universe [6][7][8].
That is, the prevailing opinion assumes DE to be a CC whilst DM is modeled as a nonrelativistic fluid.
Even though the Λ CDM model is the most favored model by the observations [3,5,10,11], it has the following drawbacks from the theoretical point of view: 1) the value of the CC is between 60 and 120 orders of magnitude smaller than what it was estimated in Particle Physics. This is known as the CC problem [12], and 2) it suffers a "Cosmological Coincidence Problem" (CCP). If the energy density of DM evolves in terms of redshift as ρ m ∝ (1 + z) 3 , and the energy density of CC, ρ Λ = Λ is constant, why do their energy densities have the same order of magnitude today? and why in the near past (z eq 0.55), the matter density had dropped to the same value as the DE density? According to the Planck 2018 results [5] we have the current values Ω m0 ≈ 0.315, Ω Λ 0 ≈ 0.685, such that z eq ≈ 0.30.
The remarkable fact that the energy densities of DE and DM are of the same order of magnitude around the present time seems to indicate that we are living in a very special moment of cosmic history. As we mentioned before, within the standard model where the DE density is constant and the DM density scales with the inverse third power of the cosmic scale factor, this appears to be a coincidence since it requires extremely fine-tuned initial conditions in the early universe. Both in the very early universe and in the far future universe arXiv:2107.04651v3 [gr-qc] 27 Nov 2021 these energy densities differ by many orders of magnitude. CCP problem was first formulated by Steinhardt [9], and explored in [13,14].
An alternative to address these problems is to consider that the DM and DE can interact with each other (for example, see [15][16][17][18][19][20][21] and a review [22]), which is a reasonable assumption considering that both fluids currently dominate the energy content of the universe [23].
Keeping in mind that the nature of the DM and DE is still an open question, it is difficult to describe these components from physical first principles. In this sense, in different descriptions of the DE that are not the CC can be treated as a fluid (for example Chaplygin gas [24]), vectors field (for example Multi-Proca vector DE [25]), scalar field (for example k-essence [26]), etc. In the scalar field description of the DM, there are several candidates defined in terms of extension of the Standard Model like the axions, which were introduced to explain the CP violation [27].
Both axion and axion-like are pseudoscalars fields with the same properties, like the periodicity φ → φ + 2πC and the shift symmetry φ → φ +C, but in the axion model, their mass is related to their decay constant, while for axion-like scalars their mass and their decay constant are not linked [34], and therefore there is a family of axion-like pseudoscalar particles. Following this line, axionic models naturally appear in stringy motivated cosmology where two axionic scalar fields are employed. In particular, a generalization of this model is an assisted inflationary scenario, namely the so-called 'N-flation' [35] where N scalar fields are employed. The essential idea is that the 'inflaton', the scalar degree of freedom that drives the inflation of the early universe, is not a single scalar field but a collection of many axionic scalar fields. More recently, such a model has been employed to construct a late time interacting DM-DE scenario [36], where the DE and part of the DM that interacts with it were represented by two axionic scalar fields respectively.
It is important to mention that cosmology with several scalar fields has been studied for quite some time now. The theory of slow-roll inflation driven by an arbitrary number of scalar fields interacting through gravity only was developed in [37]. The detailed analysis of inflation driven by two massive scalar fields was presented in [38]. The case of exponential potential and interaction is well studied in [39,40]. In the pioneering paper [41], it was shown Einstein equations, with quan-tum one-loop contributions of conformally covariant matter fields, do admit a class of nonsingular isotropic homogeneous solutions, which corresponds to an early time de Sitter (inflationary) state. Several models are realizing the inflationary paradigm, most of which include scalar field(s). The scalar field models can be studied using local and global variables, providing a qualitative description of the solution space. In addition, it is possible to provide precise schemes to find analytical approximations of the solutions, as well as exact solutions or solutions in quadrature by choosing various approaches (see for example [21,). In particular, relevant information about the properties of the flow associated with an autonomous system of ordinary differential equations can be obtained by using qualitative techniques of dynamical systems. Textbooks related to qualitative theory of differential equations can be found in [126][127][128][129][130][131][132][133][134][135], with some applications in cosmology in [136][137][138][139][140][141].
The main difficulties that arise using standard dynamical system approaches in the study of scalar fields are due to the oscillations entering nonlinear system through the KG equations, motivating the study of theses system using perturbation and averaging methods.
Perturbations methods and averaging methods were used, for example, in [142][143][144][145][146][147][148][149][150]. One idea is to construct a time-averaged version of the original system, solving it; the oscillations of the original system are smoothed out [150]. This can be achieved for homogeneous metrics where the Hubble parameter H plays the role of a time dependent perturbation parameter which controls the magnitude of the error between the solutions of the full and the time-averaged problems whenever H is monotonic and sign invariant, H is positive strictly decreasing in t and lim t→∞ H(t) = 0 [147,148,151]. Therefore, it is possible to obtain information about the large-time behavior of more complicated systems via an analysis of the simpler averaged system equations using dynamical systems techniques.
In the references [147][148][149] a research program "Averaging generalized scalar field cosmologies" was started. It consisted of using asymptotic methods and averaging theory [158] to explore the solution's space of scalar field cosmologies with generalized harmonic potential in vacuum or minimally coupled to matter. As a difference with [147,148], where the Hubble parameter was used as a time-dependent perturbation parameter, in [149] systems where Hubble scalar is not monotonic were studied.
In this paper we generalize the program initiated in [147][148][149] to two-fluid cosmological models. To this goal, we analyze an axion-like coupled scalar fields model following Ref. [36], which have recently drawn significant interest among particle cosmologists (see Sec.2). In this model there are two canonical scalar fields φ 1 , φ 2 interacting through a term added to its po- . The decoupled part V decl. , satisfies periodicity property φ → φ + 2πC, shift symmetry φ → φ +C, and time-reversal symmetry t → −t. Time-reversal symmetry is not imposed to the coupling term V int . In the uncoupled case we identify and characterize the equilibrium points but in the coupled case we need to solve transcendental equations and a little progress is made.
Methods from the theory of averaging nonlinear dynamical systems allow us to prove that time-dependent systems and their corresponding time-averaged versions have the same late-time dynamics. The main result in this paper is that in the first approximations near H = 0 of matter density, normalized scalar field densities, and the phase variables Φ 1 and Φ 2 , defined as , where Ψ 1 and Ψ 2 are functions of φ 1 and φ 2 , through (21)), the full systems and their averaged values (with an averaged function to be properly defined) have the same limit as H → 0 (as τ → ∞). The averaged values of the phases, denoted bȳ Φ i , are zero, such that, on average,Ψ i (t) = r i sin (tω i ) for some constants r i . This is summarized in our Theorem 2. Therefore, simple time-averaged systems determine the future asymptotic behavior.
The paper is organized as follow: in section 2 we discuss a coupled effective axion-like model in flat FLRW cosmology. Then, in section 3 we introduce the model under study. We discuss a theorem based on energy density estimates in section 3.1. Numerical solutions confirming the analytical results are discussed in section 3.2. In section 4 we make a dynamical system analysis using suitably defined dimensionless dynamical variables and a dimensionless time variable. In section 5 we use the averaging techniques. In particular in section 5.1 we apply the variation of constants method to the model, and in section 5.2 we study the time-averaged system using standard techniques of dynamical systems. In section 6 we make a discussion of our results. Finally, in section 7 we present our conclusions and discuss further lines of research. Theorem 2 is proved in Appendix A. In section Appendix B numerical simulations as evidence of such behavior are presented.

Coupled effective axion-like model in flat FLRW cosmology
In this section we introduce the coupled effective axionlike model presented in [36], by considering the following Lagrangian density for two pseudoscalar fields, namely θ and σ , where the fields interact through the non-perturbative potential [152][153][154] V (θ , σ ) =µ 4 1 1 − cos where g 1,2 and h 1,2 are different axion-like decay constants, µ 1,2 are some scales non-perturbatively generated, and the unity is added in order to eradicate by hand any CC. We use units in which M −2 p = 8πG = 1. The potential (2) is a direct generalization of the single axion-like potential V (θ ) = µ 4 [1 − cos (θ /g)], which is of interest in the context of natural inflation [155], where the inflation is identified with an axion-like particle. This is because the shift symmetry θ → θ +C, C is a constant, protects the flatness of the potential from perturbative corrections. This flatness condition for the inflaton potential is necessary to the inflation take place and this matches with the primordial perturbations indicated by the CMB [152]. It is important to note that in [157] axionic DM model with a modified periodic potential for pseudoscalar field 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 small, where m A the axion rests mass. Variation of (1) with respect to θ and σ gives the Klein-Gordon equations By neglecting higher than linear order terms we havë near the minimum of the potential (θ , σ ) = (0, 0). Potential (2) has a mass matrix in the (θ , σ ) basis at (0, 0) given by that is not diagonal, and the flat direction exist when the condition on the axion decay constants g 1 h 2 = g 2 h 1 is fulfilled. The conditions g 1 h 1 µ 4 2 + g 2 h 2 µ 4 1 = 0 and Following [152] we assume g 1 = g 2 = g. Then, the condition on the axion decay constants g 1 h 2 − g 2 h 1 to be nearly zero, can be expressed as h 2 = h 1 + ε, ε 1. Now, we impose the conditions for some constants f 1 , f 2 , f 3 . First, notice that when ε = 0 and g 2 = g 1 = g, that above equations admits no solutions for σ and θ . Indeed, the linear matrix Assuming 0 < h 2 − h 1 := ε 1, and defining Then, the potential (2) takes the form where f 1 , f 2 , f 3 were chosen to have two new fields, the heavy field ψ and the light field ξ . That is, the field equations near the minimum becomes where higher order terms in ξ and ψ were dropped out. Note that we can adjust h 1 − h 2 to make the effective axion decay constant f 2 arbitrary large, and to make the real scalar field ψ a heavy field, whose evolution is dominated only by the first term in the potential (7) if µ 1 µ 2 and f 2 1, while the real scalar field ξ to be a light field, whose evolution is dominated only by the second term in the potential (7) with ψ ≈ 0. For µ 1 µ 2 we havë To deduce eq. (9), we assume Neglecting ψ in the second term of (7) and renaming ψ = φ 1 and ξ = φ 2 , we obtain the decoupled effective potential where In this paper, we study a modification of the potential (12) generated by the scalar fields φ 1 and φ 2 , following the references [36,156], in which the potential is written as where the interaction term between φ 1 and φ 2 is The heavy field is φ 1 and the light field is φ 2 , which interact through the third term presented in the potential (15). The interaction is turned on when n = 0. When n = 0 the first and third terms are merged by replacing µ 4 1 + µ 4 3 → µ 4 1 and potential (12) is recovered. To obtain the minima/maxima of V (φ 1 , φ 2 ) we have to solve the transcendental equations The set {(H, ρ, χ 1 , χ 2 , φ 1 , φ 2 ) ∈ X : H = 0} is invariant for the flow of (29) and H does not change sign. On the contrary, if there is an orbit with H(0) > 0 and H(t 1 ) < 0 for some t 1 > 0, this solution passes through the origin violating the existence and uniqueness of the solutions of a C 1 flow.

Local energy estimates
Using local energy estimates, we can prove the following theorem be the positive orbit that passes through the regular point at the time t = t 0 . Then, we have lim t→∞ (ρ(t), χ 1 (t), χ 2 (t)) = (0, 0, 0), along the positive orbit O + (x 0 ).

Numerical solutions
We integrate the equations (25), (26), (27) using the redshift z instead of the cosmic time t as the independent variable. Quantities t and z are related through the expression Moreover, where the deceleration parameter Hence, eqs. (27a)-(27b) becomes Raychaudhuri equation becomes where q can be expressed as Then, we obtain a system of differential equations for φ 1 (z), φ 2 (z), and E(z) given by (42a), (42b) and (42c) and integrate in terms of redshift z, from z = 100 to z = −1.
The parameter values f 1 = 0.1, f 2 = 0.1, µ 4 1 = 1.1H 2 0 , µ 4 2 = 10.75H 2 0 , µ 4 3 = 1.07 and n = 9 are chosen. As initial conditions for the fields we use φ 1 | z=100 = 0.155, φ 2 | z=100 = 0.7835, dφ 1 dz z=100 = 0 and dφ 2 dz z=100 = 0, such that V (φ 1 (z), φ 2 (z))| z=100 = 11.6351H 2 0 . As an initial condition for the Hubble parameter when z = 100, we take as an "educated guess" the value obtained for the Λ CDM model at z = 100. That is, we use the expression where Ω r0 = 2.469×10 −5 h −2 (1+0.2271N eff ) and The parametric curve (denoted by a blue line) V (φ 1 (z),φ 2 (z)) For the analysis we have defined Ω φ = ρ φ /3H 2 and with V (φ 1 , φ 2 ) given by equation (15), and from equation (25) we obtain the constraint Ω = 1 − Ω φ , where Ω = ρ/3H 2 . To compare with Λ CDM, we represent with dashed lines the functions Using (43), the deceleration parameter for Λ CDM model is . (47) We define the kinetic and potential terms of Ω φ as and Using theorem 1, it can be obtained Based on this, we obtain the effective CC constant Λ eff : where lim z→−1 φ 1 (z) = φ * 1 and lim z→−1 φ 2 (z) = φ * 2 should satisfy the conditions of extreme point and conditions for local minima, which are In Fig. 1 we present some numerical results obtained from the integration of the equations (42a), (42b) and (42c). In figure 1(a), we present the axion-like fields density parameter Ω φ and the matter density parameter Ω of our model and the density parameters Ω m , Ω Λ and Ω r of Λ CDM as a function of the redshift z. In figure 1(b), we present the of the effective barotropic index ω φ associated to the axion-like fields and ω Λ = −1 of Λ CDM as a function of the redshift z. In figure  1(c) the evolution of the kinetic and potential terms of Ω φ given by (48) are shown as functions of redshift z. In figure 1(d) a zoom of the evolution of kinetic terms (48a) and (48b) are presented as functions of redshift z. In figure 1(e), we present the evolution of the ratios Ω /Ω φ and Ω /Ω Λ as a function of the redshift z. In figure 1(f), we present the evolution of the dimensionless Hubble parameter E = H/H 0 for our model and for Λ CDM as a function of the redshift z. Finally, in figure 1(g), we present the evolution of the deceleration parameter q := −1 −Ḣ/H 2 for our model and for Λ CDM as a function of the redshift z. In figure 1(h) we present the evolution of the axion-like fields φ 1 and φ 2 as a function of the redshift z. According to figure 1(h) the values φ 1 = 0.061 and φ 2 = 0.572 can be used as initial points for finding numerically the roots of (52) which are found to be φ * 1 = 0.0614165, and φ * 2 = 0.572375. Evaluating (53) we have 176.099H 2 0 > 0, and 788143H 4 0 > 0. Then, as shown in figure 2, a local minimum with minimum

Dynamical systems analysis using H 0 -normalized variables
Although expansion normalized dynamical variables are commonly used for dynamical analysis of late time cosmologies, it turns out that for this model they do not lead to an autonomous dynamical system. To make an autonomous dynamical system for this model, we introduce the following H 0 -normalized variables together with the scalar field redefinition the dimensionless time variable and the following parameters In the above H 0 is the present value of the Hubble parameter. We will denote the derivative with respect to τ with a prime: (·) ≡ d/dτ. In terms of the above variables and parameters we can write down the dynamical system as We note that the dynamical system (58) is periodic in both y 1 , y 2 with a period of π as long as n is an integer. This means that for the decoupled limit (n = 0) of the model, as well as the interacting cases when n is an integer, we can confine our attention within the range −π/2 ≤ y 1 , y 2 ≤ π/2. The fixed points/lines of the above system are given by the conditions The isolated fixed points of the decoupled version of the model and their stability are listed in table 1. In the coupled case n = 0, the transcendental equations (59b) have to be solved numerically. For a given (y 1 , y 2 , E) = (y * 1 , y * 2 , E * ) which satisfies (59), it can be found univocally the values b 1 , b 2 and Table 1 Some isolated fixed points of the system (58) for n ∈ Z over (y 1 , Fixed point Coordinates (E, x 1 , x 2 , y 1 , y 2 ) Existence Stability Solution Saddle for n = 0, de Sitter depends on n otherwise Saddle for n = 0, de Sitter depends on n otherwise However, to proceed further is a hard task and little progress can be made.
In the first quadrant (y 1 ≥ 0, y 2 ≥ 0) there are two representations of P (0,±1) module 2π: Corresponding to the point P (0,0) we get a Minkowski solution because we have chosen the potential (15) to have the minimum value zero. Adding a cosmological constant term and provided it is positive, the stable fixed point P (0,0) will corresponds to a de Sitter solution.

Coupled case n = 0
For n = 0, the equilibrium points discussed in table 1 are equilibrium points of (58), satisfying conditions (59). The other possible equilibrium points are found by reducing the conditions (59), by solving the transcendental equations (59b) numerically using the parameter values: The conditions (59b) and (59c) becomes 107 sin(2(y 1 − ny 2 )) + 110 sin(2y 1 ) = 0, Therefore, in the coupled case we need to solve transcendental equations and a little progress is made. In Fig. 3 are presented the contour domains of the restrictions given by the equations (63) and (64) for the coupled case with n = 1, 2, 3 and 4. At the points where the solid red line and the dashed blue lines coincide are found the equilibrium points of the system.

Reduced dynamical system
All the fixed points of the system (58) must necessarily reside on the 3-dimensional hypersurface specified by (x 1 , x 2 ) = (0, 0), where we can express x 1 , x 2 in terms of the other variables as for x 1r = 0, x 2r = 0, E r = 0, where we use the subscript r to denote the value of a variable on the 3surface. Defining η = dτ 6E r , one can reduce the original 5-dimensional system (58) to a 3-dimensional system on this 3-surface: In the figure 4 is is represented a projection of the phase portrait for n = 0 (decoupled limit) and for n = 1, n = 2 and n = 3, on the 2-dimensional field space The red dot at the center being the origin corresponds to Minkowski solution. This extended domain is taken to show the periodic nature of the phase flow with a period π. The parameter values used for the plots are a 1 = 1 300 , a 2 = Using the fact that t (τ) = H −1 0 and defining τ 0 as the current value associated to t = t a ("age of the universe"), we have Finally, as t → ∞, we have → 0.

Oscillatory behavior: averaging
Let us assume n = 0. Defining the new fields (21) the field equations becomë This domain is taken to show the periodic nature of the phase flow with a period π 2 for

The Friedmann equation becomes
and the deceleration parameter is Imposing the conditions where ∆ = a 2 (b 1 + b 3 ) − a 1 b 2 + b 3 n 2 2 + 4a 1 a 2 b 2 3 n 2 , we obtain the decoupled oscillators such that Ψ 1 is a massive scalar field and Ψ 2 a lighter one. Defining the new variables θ 1 , θ 2 and T , through and the time derivative we obtain in the limit Ψ i → 0 the reduced dynamical system where we denoteω 1 = ω 1 H 0 ,ω 2 = ω 2 H 0 and we have used the fact that q → 1 2 as Ψ i → 0.
That means that the decoupled oscillators mimic dust matter. The dynamical system (70) is regular.

Variation of constants
According to the Raychaudhuri equation (26), H is a monotonic decreasing function. Additionally, as the Therefore, as H → 0, we obtain the decoupled oscillators (71). Motivated by the solution (72), we use the variation of constants to propose the solution of the full KG system as with inverse functions where c and ω 1 , ω 2 are undetermined constants.
The full system is given bẏ It is worth noticing that by conveniently choosing c, ω 1 , ω 2 as in (69), the equations up to linear order in r 1 , r 2 can be reduced to the usual equations for two decoupled harmonic oscillators. Let us define We obtain the full system: where the deceleration parameter is expressed as The Friedman equation is transformed to All the terms in the above equations are non-negative. Then, it follows that 0 ≤ Ω , Ω 1 , Ω 2 ≤ 1.
Imposing the conditions (69), or alternatively, imposing the conditions which admit as special solution (69), the undesired terms of zero order in H are eliminated. That is, by defining x = (Ω 1 , Ω 2 , Ω , Φ 1 , Φ 2 ) T , by expanding in Taylor's series around H = 0 we have the 6-dimensional systeṁ as H → 0, where Following arguments stressed in [151] but for the twoscalar field setting, since we obtain the equations of motion of two decoupled oscillators (71) as H → 0, we can study the full system as perturbed harmonic oscillators and apply averaging techniques for vector functions f n (t, x) which have two independent periods T n , n = 1, 2, where H is considered as a time-dependent and itself is governed by the evolution equation (80b). A surprising feature of such an 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; 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 the large-time behavior of a more complicated full system via an analysis of simpler averaged system equations. Indeed, from G [2] (t, x) ≥ 0, it follows that H is a decreasing function of time, that allows defining a decreasing sequence of parameters to construct asymptotic expansions. Therefore, although H is a function of time and it is not properly a "constant parameter", as t → 0, we can choose a sequence t n and for n large enough such that for H(t n ) = H n , which is actually a parameter. This procedure is supported by numerical simulations (see related work [150]). Similar arguments can be provided using the tools of [143] or, as is the case in this research, by implementing an analogous program as in the references [147][148][149] but for averaging multi-scalar fields cosmologies.
Using the first condition, the secular terms proportional to t are eliminated, and using the second condition we have c 3 + c 4 + c 5 = 0. These conditions, however, limit the applicability of the model.
In general, the regular asymptotic expansion fails in presence of resonant terms. One alternative can be using the method of multiples scales, e.g., a Poincarè-Lindstedt's -like method, where we set t 1 = (1 + ω 1 h + ω 2 h 2 +. . .)t, t 2 = ht, where h = H/H 0 1. This method would determine solutions of perturbed oscillators by suppressing resonant forcing terms that would yield spurious secular terms in the asymptotic expansions. The t 1 and t 2 time variables are introduced to keep a well ordered expansion, where t 1 is the regular (or "fast") time variable and t 2 is a new variable describing the "slow-time" dependence of the solution. The idea is to use any freedom that is in the t-dependence of t 1 and t 2 to minimize the approximation's error, and whenever is possible to remove unbounded or secular terms. To our knowledge, this method has not been implemented yet in the cosmological setup. However, basic examples of oscillators show that by implementing a time-averaged version of the model instead of multiple scales, the issue of secular terms is overcome; getting the same accuracy as in the two-timing method. We elaborate more on averaging techniques in subsection 5.1.2.

Time-averaging
If we have vector functions f n (t, x) which have N independent periods T n , n = 1, . . . , N, we take the averaginġ where where y is considered as a parameter that is kept constant during integration.
Assuming ω 1 = ω 2 , with H playing the role of ε, N = 2, T 1 = 2π ω 1 and T 2 = 2π ω 2 we can use the following averaging procedure: where the vector field f (t, y) in the right-hand-side of the equation must be the sum of two vector functions f 1 (t, y) and f 2 (t, y) where each of them is periodic with one period.
The averaged system obtained using such approach is  Proceeding in an analogous way as in references [143,147,148,159] we implement a local nonlinear transformation: Taking time derivative in both sides of (92) 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 (91a) and (92) in (94) we obtain where I 5 is the 5 × 5 identity matrix. Then we obtaiṅ Using eq. (91b), we haveḢ = O(H 2 ). Hence, The strategy is to use eq. (98) for choosing conveniently ∂ ∂t g(0, x 0 ,t) to prove thaṫ is unknown at this stage and it is independent of periods T 1 and T 2 due to it is independent of t. 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. (98). 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 and we have assumed that the function in the left hand side denoted by f(x 0 ,t) −f(x) can be separated in the sum of two vector functions with independent periods T 1 and T 2 . Hence, the right hand side of (100) is the sum of two almost periodic functions f 1 (x 0 ,t) −f 1 (x) and f 2 (x 0 ,t) −f 2 (x) of independent periods L 1 = 2π ω 1 and L 2 = 2π ω 2 for large times, respectively. Then, implementing the average process (90) on right hand side of (100), where slow-varying dependence of quantities x 0 andx on t are ignored through averaging process, we obtain Defining the average (102) is zero so that g(x 0 ,t) is bounded. Finally, eq. (99) transforms tȯ and eq. (100) is simplified to Theorem 2 establishes the existence of the vector (93).

Phase space analysis of the time-averaged systems
Introducing the time variable τ = τ, we obtain the guiding system: The equilibrium points of the guiding system (106) are the origin, with eigenvalues −3, −3, −3 and the points in the planeΩ 1 +Ω 2 +Ω = 1 with eigenvalues {3, 0, 0}. Therefore, that are nonhyperbolic, with a 1D unstable manifold and a 2D center manifold. Taking an arbitrary point Ω 1 ,Ω 2 ,Ω = Ω 1c ,Ω 2c , 1 −Ω 1c −Ω 2c , Ω > 0 lying on this plane, and defining the new variables that translates the equilibrium point to the origin, we obtain the system  Fig. 5 Orbits of the phase space of system (117) Hence, the unstable manifold is given locally by the graph Using the invariance of the unstable manifold we obtain for U = 0 that The solutions are the trivial h 1 (U) = h 2 (U) = 0 or The last solutions do not satisfy the condition h 1 (0) = 0, h 2 (0) = 0, h 1 (0) = 0, h 2 (0) = 0. Hence, the unstable manifold is the U-axis and the dynamics on the unstable manifold is given by U = 3U. The center manifold is given locally by the graph Using the invariance of the center manifold we obtain that h 3 (V 1 ,V 2 ) satisfies the quasilinear partial differential equation This equation have the trivial solution h 3 (V 1 ,V 2 ) = 0 and the general solution which do not satisfy the conditions h 3 (0, 0) = 0, and Dh 3 (0, 0) = 0. Then, the center manifold is the plane U = 0. The guiding equation can be studied by defininḡ x =Ω andȳ =Ω 1 +Ω 2 +Ω , we obtain the reduced 2D system This system (117) have the equilibrium point (x,ȳ) with eigenvalues −3, −3. Then it is a sink. Additionally, we have the line of equilibrium pointsȳ = 1 with eigenvalues 0, 3 which is normally hyperbolic and a source. In figure 5 some orbits of the phase space of system (117) are presented.

Discussions
In section 2 we presented a coupled effective axionlike model in flat FLRW cosmology. In particular, we discussed the mechanism to make the effective axion decay constant f 2 arbitrary large, and to make the real scalar field ψ a heavy field, whose evolution is dominated only by the first term in the potential (7) if µ 1 µ 2 and f 2 1, while the real scalar field ξ to be a light field, whose evolution is dominated only by the second term in the potential (7) with ψ ≈ 0.
For µ 1 µ 2 we havë Neglecting ψ in the second term of (7) and renaming ψ = φ 1 and ξ = φ 2 , we obtain the effective potential (12). This potential was generalized to (15), where the heavy field is φ 1 and the light field is φ 2 , which interact through the third term presented in the above potential. The interaction is turned on when n = 0. When n = 0 the first and third terms are merged by replacing µ 4 1 + µ 4 3 → µ 4 1 and potential (12) is recovered. Then, in section 3.1 we have discussed the theorem 1 based on energy density estimates. Using results (32) and (33) of theorem 1, we generically have That is, the kinetic terms in Ω φ and the matter-energy density tend to zero. Then, the cosmological solution is dominated by the potential term. Numerical simulations of the system (26)- (27) were discussed in section 3.2. Initial conditions φ 1 | z=100 = 0.155, φ 2 | z=100 = 0.7835, dφ 1 dz z=100 = 0 and dφ 2 dz z=100 = 0 were considered. The initial value H| z=100 is estimated from expression (43). We see that the Λ CDM is recovered and in this case a local (not zero) minimum of the potential is approached.
The parameter values and initial conditions were chosen for obtaining a ratio of the DM and DE densities at z = 0 equal to r| z=0 = 0.315 0.685 = 63 137 ≈ 0.459854. A late-times r → 0, Ω φ dominates, in particular That is, for this choice of parameters and initial conditions, our model it suffers of the CCP, since it is not indistinguishable from Λ CDM. A subtle difference, is in that lim z→−1 W int = 0.101566, which means that ≈ 1% of the total dimensionless energy density corresponds to interaction between the two fields φ 1 and φ 2 .
It is plausible to think of the light scalar field φ 2 as DE, and the heavy field φ 1 , as an axion-like DM which interacts through the interaction term V int (φ 1 , φ 2 ) defined by (16).
Increasing the contribution of V int at late-times, or introducing an explicit interaction term between the axion-like part and CDM, the CCP can be alleviated or solved (e.g., as in quintessence and phantom field scenarios in [15][16][17][18][19][20][21]. If V int (φ 1 , φ 2 ) ≡ 0, and introducing an effective interaction Q = 2 ln χ(φ 1 (τ), φ 2 (τ)), we obtain, from the balance of energy, the modified KG equations, and modified continuity equation given by: Then, defining τ = ln(a/a 0 ), we have As V 3H 2 → 1 (according to theorem 1 the potential energy dominates at late-times), we have Then, . Therefore, at latetimes r → r 0 = δ 2 1−δ 2 . Then, the CCP is solved due to Using Planck 2018 results [5] we have the current values Ω m0 ≈ 0.315, Ω Λ 0 ≈ 0.685, and, at late times, Setting, δ = 0.2 we obtain the equality DM-DE epoch approximately at z eq ≈ 0.547797. The approximations are valid for z < z crit = 1.2067. If we integrate the full equations (118) and (119), the exact expressions for r, Ω φ and Ω m ares less sharp than displayed in figure  6 for z ≈ 1, but the approximated solutions and the numerical exact ones matches very well as z → −1 where the potential energy start dominating. Then, (120) becomes an accurate approximation of (119). Figure 6 shows the evolution of (123) in terms of redshift, which give the asymptotic values of the ratio r and Ω φ , Ω m under the interacting scheme (118) with Q = 3(δ 2 + 1) ln(a/a 0 ) = −3(δ 2 + 1) ln(1 + z). This kind of interaction indicates a phenomenological solution to the CCP. Actually, once the universe reaches the stable Ω φ -dominated state with constant ratio r = r 0 , it will live in this state for a very long time. Then, it is not a coincidence to live in this long-living state where r 0 ∼ 1 implies ρ ∼ ρ φ . The full analysis is not covered in the present research. In a forthcoming paper we will develop arguments like in [15][16][17][18][19][20][21] for the two-field scenario.
Moreover, we introduced in section 4 novel dynamical variables and dimensionless time variables, which have not been used in analyzing these cosmological dynamics where expansion normalized dynamical variables are usually adopted. In the uncoupled case the equilibrium points can be fully identified and characterized. However, in the coupled case we need to solve transcendental equations and a little progress is made. The main difficulties that arise using standard dynamical systems approaches are due to the oscillations entering the nonlinear system through the Klein-Gordon (KG) equations. This motivated the analysis of the oscillations using averaging techniques in section 5. In particular, integrating (91) we use the definition Ω t :=Ω 1 +Ω 2 +Ω .
Then, we acquire the equations Then, assumingΩ t (0) < 1, we obtain . (126e) Hence, on average, H is decreasing from That is, as far asΩ t (0) < 1 the universe ends on a de Sitter inflationary phase.
Asymptotically we have dH dτ = − 3 2 H. Defining τ 0 as the current value associated to t = t a ("age of the universe"), which as before can be set to τ 0 = 0, we have H(τ) = H 0 e − 3 2 τ . Substituting this expression for H in (127) and integrating, we obtain t(τ) = 3H 0 t a +2e 3τ/2 −2 Finally, combining all we have H = 2H 0 2+3H 0 (t−t a ) which tends to zero as t → ∞. That is, the universe passes through a matter-dominated phase before reaching a Minkowski stage.
We have shown the application of methods from the theory of averaging nonlinear dynamical systems allows us 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. Then, we can study the time-averaged system using standard techniques of dynamical systems.

Conclusions
In this paper, we have analyzed the coupled axion-like model following Ref. [36] consisting of two canonical scalar fields φ 1 , φ 2 interacting via the potential (16). We have introduced dimensionless dynamical variables and dimensionless time variables. In the uncoupled case the equilibrium points were fully identified and characterized. However, in the coupled case we need to solve transcendental equations and a little progress is made. The main difficulties are due to the oscillations entering the nonlinear system through KG equations.
Using local energy estimates, we have proved theorem 1. This result shows that if O + (x 0 ) is the positive orbit that passes through the regular point x 0 defined as in (31), then, since H is positive and decreases along O + (x 0 ), the limit lim t→∞ H(t) exists, and it is a nonnegative number η. Then, we have lim t→∞ (ρ, χ 1 , χ 2 ) = (0, 0, 0) and 3η 2 = lim t→∞ V (φ 1 (t), φ 2 (t)). Depending on initial conditions, the solution tends to a constant , related with a local minimum φ * = (φ * 1 , φ * 2 ) with non zero minimum value of V (φ 1 , φ 2 ) satisfying eqs. (17). They correspond to de Sitter solutions a ∝ exp − t . However, choosing initial conditions with a small enough value of 3H 2 (t 0 ) the orbit is trapped by the basin of attraction of the global minimum φ * = (0, 0), H * = 0. Then, the Minkowski solution is the late-time attractor.
The main result in this paper is our Theorem 2. It states that in the first-order approximations of normalized scalar field densities, and the values of two scalar fields Φ 1 and Φ 2 , defined as Φ i = tω i − tan −1 ω i Ψ i Ψ i (called phase variables, where Ψ 1 and Ψ 2 are functions of φ 1 and φ 2 through (21)) near H = 0, the full systems and their averaged values (with an averaged function to be properly defined) have the same limit as H → 0 (as τ → ∞). The averaged values of the phases, denoted bȳ Φ i , are zero, such that, on average,Ψ i (t) = r i sin (tω i ) for some constants r i . Therefore, with this approach, oscillations entering the nonlinear system through the KG equation can be controlled and smoothed out as the Hubble factor H, acting as a time-dependent perturbation parameter, tends monotonically to zero. We have studied the time-averaged system using standard techniques of dynamical systems and numerical simulations are presented as evidence of such behavior.
This approach has potential applications in physical situations where it is necessary to consider the time variation of the fields in the vicinity of the potential minimum. One such situation is the reheating after inflation. During reheating the inflaton scalar field oscillates around the potential minimum. For nonzero H this gives rise to time-dependent oscillatory dynamics, which is responsible for particle production via quantum field theory. The relevant calculations depend on the form of the potential, and in particular, are quite complicated for harmonic potentials. The result presented here shows that one can "average out" the oscillations arising due to the harmonic functions, thus simplifying the problem. Indeed, by using some inverse transformations, one can find from theΩ i toφ i , i.e., the averaged version of the original field variables. We hope that this approach may help calculations of reheating in the context of N-inflation model [35].
This approach can also be useful if one wants to consider linear cosmological perturbations for coupled axion cosmological models near the potential minimum, which is an attractor. In the cosmological perturbation theory, cosmological perturbations at the linear level are governed by equations whose coefficients are composed of background quantities. Therefore, proper knowledge of the background dynamics is necessary for further perturbative analyses. The result presented in this paper simplifies the dynamical analysis of the background near this attractor, which in turn facilitates a subsequent analysis of cosmological perturbations using procedures similar to those used in [160] (chapt. 14) and in [84,[161][162][163][164][165]. ciación grant no. 11180126. Additionally, this research is funded by Vicerrectoría de Investigación y Desarrollo Tecnológico at Universidad Católica del Norte. The work of Bin Wang was partially supported by the key project of NNSFC under contract 11835009. Samuel Lepe is acknowledged for discussions. We thank an anonymous referee for valuable comments which helped improve our work.
Appendix A: Proof of Theorem 2 Lemma 3 (Gronwall's Lemma (Integral form)) Let be ξ (t) a nonnegative function, summable over [0, T ] which satisfies almost everywhere the integral inequality Then, ξ (t) ≤ C 2 e C 1 t , almost everywhere for t in 0 ≤ t ≤ T . In particular, if Lemma 4 (Mean value theorem) Let U ⊂ R n be open, f : U → R m continuously differentiable and z ∈ U, h ∈ R m vectors such that the line segment z + η h, 0 ≤ η ≤ 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.
for all t ∈ [t n ,t n+1 ] such that t − t n ≤ t n+1 − t n = 1 H n . Then Finally, taking the limit as n → ∞, we obtain H n → 0 and Result 1 and Result 2 follow.
(B.22) and the averaged system: The system (B.22) was integrated in the interval −40 ≤ τ ≤ 10, and the system (B.23) was integrated in the interval −40 ≤ τ ≤ 40, partitioned in 10000 data points. The values of ω 1 = √ 2 and ω 2 = √ 2/2, and the initial conditions presented in Table 2, were considered. Table 2 Five initial data sets for the simulation of the truncated system (B.22) and time-averaged system (B.23). All the conditions are chosen in order to fulfill the inequalityΩ 1 +Ω 2 +Ω ≤ 1.  and ω 2 = √ 2/2. We use as initial condition for both system the five data sets presented in Table 2. These plots are numerical evidence that the main Theorem 2 is fulfilled. That is, the solution of the truncated system follow the track of the solutions of the time-averaged system and the oscillations experimented by the truncated system are smoothed out as H → 0. In figure 7 are presented the numerical results of the integration of the truncated system (B.22) (blue lines) and time-averaged system (B.23) (orange lines), for the initial conditions presented in the Table 2. In figure 7(a) are presented the solutions in the projection (Ω , H, Ω t ) while in figure 7(b) are presented the solution in the projection (Ω , Ω t ), where Ω t = Ω 1 + Ω 2 + Ω . These figures represent the numerical evidence that the main Theorem 2, presented in section 5.1.2, is fulfilled. That is, the solution of the truncated system follows the track of the solutions of the time-averaged system, and the oscillations experimented by the truncated system are smoothed out as H → 0, having both systems of differential equations the same late-time dynamics. On average, H is decreasing from H = ∞ as τ → −∞ to H = H 0 1 −Ω t (0) as τ → +∞. That is, as far as Ω t (0) < 1 the universe ends on a de Sitter inflationary phase. Then, depending on initial conditions, both the non-oscillating curve and the oscillating one tends to the same constant value of H, and it is related to a de Sitter solution. Recall, each local minimum φ * = (φ * 1 , φ * 2 ) with non zero minimum value of V (φ 1 , φ 2 ), satisfying eqs. (17), corresponds to de Sitter solutions t . On the other hand, given that Ω t (0) = 1 is an invariant set, thenΩ =Ω (0),Ω = Ω (0),Ω =Ω (0) remains constant withΩ t (0) = 1. Finally, if we chose initial conditions for a small enough value of 3H 2 (t 0 ) such that the orbit is trapped by the basin of attraction of the global zero minimum, we obtain a Minkowski solution.