Averaging with a time-dependent perturbation parameter

The theory of averaging in nonlinear dynamical systems has its origins in the 18th century, when perturbative methods for differential equations became important to test the predictions of Newtons theory of gravitation for celestial mechanics against observational data. The motion of a planet is described by a Keplerian ellipse only in good approximation, since it represents the solution to the two-body problem planet-sun. The actual trajectory of a planet is however also influenced by all other celestial bodies such as their own satellites. The influence of these bodies results in (generally periodic) perturbations from Keplerian motion. The necessity to take into account these effects due to the advances in accuracy of astronomical observations lead to the development of perturbative methods for differential equations, and eventually to the birth of averaging methods in nonlinear dynamical systems; cf [22, appendix A].

Nowadays these methods represent widely used tools in science and engineering, backed by an increasingly well developed mathematical theory; cf e.g. [22]. The core idea is to approximate a perturbed full system by an unperturbed averaged system. The latter is thereby obtained from the former by an integration over a scale associated with the perturbation; e.g. a period in the case of a periodic perturbation. The basic theorems then quantify the accuracy of the averaged approximation in terms of error estimates; cf e.g. [22, p 31, theorem 2.8.1 or p 101, theorem 5.5.1]. Typically the perturbation is controlled by a constant parameter, and the error estimates together with their range of validity are expressed in terms of orders of it.

With respect to modern gravitational physics, averaging often refers to spatial averaging in inhomogenous cosmology, a subject studied in particular with regards to the backreaction source terms appearing in the averaged equations, which potentially provide a natural explanation for dark energy; cf e.g. [9, 10]. Here we however consider spatially homoegenous cosmology and temporal averaging over oscillatory perturbations. The latter occur for instance when the spacetime is populated by matter which satisfies a wave equation, such as scalar fields (Klein–Gordon equation) or Yang–Mills fields; cf e.g. [3, 5, 6, 12, 13]. Yet such oscillatory behaviour can even be present in perfect-fluid cosmologies, as for instance demonstrated by the cases which exhibit a phenomenon called Weyl curvature dominance or asymptotic self-similarity breaking; cf [15, 23, 24].

A focus point of spatially homogenous cosmology is the investigation of the past and future asymptotics of these systems by analytical means, or, expressed in the language of the analysis of PDEs, their stability; cf [1, 2, 4, 8, 11, 16, 18, 21, 25] for general reference and also [7, 14] for the stability of some prominent spacetimes. In other words, one is interested in how the cosmological models evolved in a distant past (e.g. towards a big-bang singularity) or distant future (e.g. towards a big-crunch singularity or in a forever expanding state). The present work is motivated by [12] which offers an analysis of the future asymptotics of LRS Bianchi III cosmologies with a massive scalar field. Type III thereby refers to the Bianchi type of the spatial homogeneity while LRS (locally rotationally symmetric) refers to an additional spatial rotation symmetry; cf [21, 25]. The scalar field ϕ = ϕ(t) satisfies the Klein–Gordon equation

$$\begin{equation}\ddot {\phi }+\phi =H\left(-3\dot {\phi }\right),\end{equation} \tag{ 1.1 }$$

in which the mass is set to 1 and where a dot denotes differentiation with respect to coordinate time t. H = H(t) is the Hubble scalar—a measure of the overall isotropic rate of spatial expansion; cf [16, 25]. Evidently, this cosmology is an example for one of the systems mentioned above, in which oscillations enter via a matter wave equation.

From the point of view of (1.1), H represents a time-dependent damping of the harmonically oscillating matter field ϕ. Its time dependence is governed by the coupling to the Einstein equations, which can be expressed as a set of constraints and first-order evolution equations for H and other geometric quantities ³ . From the explicit form of the Raychaudhuri equation (i.e. the evolution equation for H), it is easy to show that (for an initially expanding Universe) H is strictly decreasing and furthermore

$$\begin{equation}\underset{t\to \infty }{\mathrm{lim}}\enspace H\left(t\right)=0;\end{equation} \tag{ 1.2 }$$

cf [12, lemma 3] ⁴ . This guaranteed smallness of H for large times was motivation in [12] to associate (1.1) with the simpler problem

$$\begin{equation}\ddot {\phi }+\phi ={\epsilon}\enspace \left(-3\dot {\phi }\right)\end{equation} \tag{ 1.3 }$$

in which H(t) is replaced by a small perturbation constant > 0. For problems of this kind the theory of (periodic) averaging in non-linear dynamical systems offers powerful theorems; cf [22]. For (1.3) these provide error estimates between ϕ and the solution of an associated averaged equation in terms of orders of the perturbation constant . In the light of (1.2) and the association of (1.1) with (1.3) it then seems natural to conjecture that the error between the solution of (1.1) and its associated time averaged equation decays and goes to zero. If this holds, then one can obtain the future asymptotics of the full oscillatory system by analysing the much simpler averaged system, which is free of oscillations.

The above paragraph summarises the core idea entertained in [12], in which a respective conjecture is formulated and supported. A closely related approach has also been taken independently in [3, 5, 6], where [3, 5] even provide theorems concerning the future asymptotics. However, the applicability of these theorems is limited to the concrete cosmological models and matter models discussed in these sources, as is also the conjecture of [12]. Yet oscillatory problems of this kind are a reoccurring theme in mathematical cosmology, at least (but not only) w. r. t. matter fields which satisfy a wave equation. Thus, instead of having to prove a separate theorem for each symmetry class and matter model, one would ideally like to have only a short number of theorems available which are however applicable to a wider range of models. With the present work we aim to make a contribution in this spirit. With theorem 3.1 we provide a rigorous result concerning the large-time behaviour of a wide class of cosmological models abstracted from the system of [12], i.e. from the LRS Bianchi III Einstein-Klein-Gordon system. The statement of the theorem is stronger than that of the theorems of [3, 5] and the conjecture of [12] w. r. t. its generality and range of applicability. It is weaker in that it does not provide the solutions in the limit for t to infinity, but rather gives error estimates for finite time intervals. However these time intervals increase, and the estimated errors decrease, with the time of truncation from the full to the averaged system. Since the latter can be chosen freely, our theorem 3.1 should be a versatile and useful tool to investigate the future asymptotics of cosmologies with oscillatory behaviour by analytical means.

Lastly, we note that while we do not explicitly exclude a cosmological constant for the class of models which satisfy the premises of our theorem, it is likely that the requirement that H → 0 does not hold in these cases. The importance of this decay of H for controlling the error between the full and the averaged solutions is also why we restrict to the future asymptotics here, since towards the past, H and thus the error is increasing in the motivating model of [12].

The paper is structured as follows. In section 2 we give some background on periodic averaging in non-linear dynamical systems, and we relate our case to the standard problems of this theory. We then formulate our result in section 3: a theorem on the large-time behaviour of a general class of systems, tailored to problems in spatially homogenous cosmology with oscillatory behaviour. We provide the proof in section 4 and conclude with a summary in section 5.

Throughout the paper bold case quantities denote vectors whereas quantities in italics denote scalars. We denote the 0th to ith-order coefficients and the coefficient of the ith-order remainder term in the Taylor expansion of $f\left(\mathbf{x},t,{\epsilon}\right)$ near = 0 as ${f}^{0}\left(\mathbf{x},t\right),\dots ,{f}^{i}\left(\mathbf{x},t\right)$ and f^[i+1](x, t, ), respectively, such that

$$\begin{equation*}f\left(\mathbf{x},t,{\epsilon}\right)=\sum\limits _{j=0}^{i}{f}^{j}\left(\mathbf{x},t\right){{\epsilon}}^{j}+{f}^{\left[i+1\right]}\left(\mathbf{x},t,{\epsilon}\right){{\epsilon}}^{i+1};\end{equation*}$$

cf also [22, p 13, notation 1.5.2]. The norm || ⋅|| denotes the standard discrete ℓ¹ norm

$$\begin{equation*}{\Vert}\mathbf{u}{\Vert}\enspace \text{def}{=}\enspace \sum\limits _{i=1}^{n}\vert {u}_{i}\vert \end{equation*}$$

for $\mathbf{u}\in {\mathbb{R}}^{n}$. We also denote as ${L}_{\mathbf{x},t}^{\infty }$ (or simply L^∞) the standard L^∞ space in both x and t variables with the corresponding norm defined as

$$\begin{equation*}{\Vert}f{{\Vert}}_{{L}_{\mathbf{x},t}^{\infty }}\enspace \text{def}{=}\enspace {\mathrm{sup}}_{\mathbf{x},t}\vert f\left(\mathbf{x},t\right)\vert .\end{equation*}$$

For our purpose it is instructive to compare (1.1) to a general class of Van der Pol equations

$$\begin{equation}\ddot {\phi }+\phi ={\epsilon}\enspace g\left(\phi ,\dot {\phi }\right),\end{equation} \tag{ 2.1 }$$

with a sufficiently smooth function g and where is a constant, in contrast to H = H(t) in (1.1); cf [22, section 2.2]. For small , (2.1) represents a perturbed harmonic oscillator which suggests a parametrisation of the solution by variation of arbitrary constants in terms of an amplitude r = r(t) and phase φ = φ(t):

$$\begin{equation}\phi =r\enspace \mathrm{sin}\left(t-\varphi \right)\quad \text{and}\quad \dot {\phi }=r\enspace \mathrm{cos}\left(t-\varphi \right).\end{equation} \tag{ 2.2 }$$

This transforms (2.1) into two first-order equations of the form ${\left[\dot {r},\dot {\varphi }\right]}^{T}={\epsilon}\enspace \mathbf{f}\left(r,\varphi ,t\right)$, with f2π-periodic in t. From [22, chapter 2] this is a periodic averaging problem of the more general standard form

$$\begin{equation}\dot {\mathbf{x}}={\epsilon}\enspace {\mathbf{f}}^{1}\left(\mathbf{x},t\right)+{{\epsilon}}^{2}\enspace {\mathbf{f}\enspace }^{\left[2\right]}\left(\mathbf{x},t,{\epsilon}\right),\quad \mathbf{x}\left(0\right)=\mathbf{a},\end{equation} \tag{ 2.3 }$$

with f T-periodic in t.

For this kind of problem the theory of averaging provides theorems which estimate the error between the solution x(t) of (2.3) and a z(t) which solves the corresponding averaged equation

$$\begin{equation*}\dot {\mathbf{z}}={\epsilon}\enspace {\bar{\mathbf{f}}}^{1}\left(\mathbf{z}\right),\quad \mathbf{z}\left(0\right)=\mathbf{a},\quad \text{where}\quad {\bar{\mathbf{f}}}^{1}\left(\mathbf{z}\right){:=}\frac{1}{T}{\int }_{0}^{T}{\mathbf{f}}^{1}\left(\mathbf{z},s\right)\enspace \mathbf{d}s.\end{equation*}$$

For instance [22, p 31, theorem 2.8.1] states that the error between x(t) and z(t) is O() on timescales of O(⁻¹); cf also lemma 4.2 below. Note that both the error estimate and its domain of validity are given in terms of orders of powers of the perturbation constant .

Remark 1. The averaging can be interpreted as a time dependent coordinate transform, a nonlinear transform for nonlinear problems, and many of the results apply to deterministic time-dependent perturbations; cf [20, chapter 3] and [19, section 5].

The similarity between (1.1) and (2.1) suggests to treat the former as a perturbed harmonic oscillator as well, and to apply averaging in an analogous way. Care has to be taken however since in contrast to , H is time-dependent and itself subject to an evolution equation of the coupled system. If valid, then a striking feature of such an approach is the possibility to exploit the fact (1.2) that H is strictly decreasing and going to zero. In (1.1) H takes the role of the perturbation parameter in (2.1), which controls the error and the domain of validity of the error estimate. Hence, with strictly decreasing H the error should decrease as well. One can obtain the information about the large-time behaviour of the more complicated full system via an analysis of the simpler averaged system.

This ideas just summarised were already entertained in [12] at the example of LRS Bianchi III. Likewise [3, 5, 6] took a closely related approach for the concrete models of these sources. [3, 5] thereby give theorems concerning the respective future asymptotics.

The purpose of the present work is to formulate a theorem concerning the large-time behaviour which is applicable to a wide range of problems in spatially homogeneous cosmology. We therefore define our problem as general as possible while also as restrictive as necessary to exploit certain key features.

The systems that we consider in the following have the general form

$$\begin{equation}\left[\begin{matrix}\hfill \dot {H}\hfill \\ \hfill \dot {\mathbf{x}}\hfill \end{matrix}\right]=H{\mathbf{F}}^{1}\left(\mathbf{x},t\right)+{H}^{2}{\mathbf{F}}^{\left[2\right]}\left(\mathbf{x},t\right)=H\left[\begin{matrix}\hfill 0\hfill \\ \hfill {\mathbf{f}}^{1}\left(\mathbf{x},t\right)\hfill \end{matrix}\right]+{H}^{2}\left[\begin{matrix}\hfill {f}^{\left[2\right]}\left(\mathbf{x},t\right)\hfill \\ \hfill 0\hfill \end{matrix}\right],\end{equation} \tag{ 2.4 }$$

where H = H(t) is positive, is strictly decreasing in t, and

$$\begin{equation}\underset{t\to \infty }{\mathrm{lim}}\enspace H\left(t\right)=0.\end{equation} \tag{ 2.5 }$$

Example 1. Our prime example for a system which fits the general form (2.4) is the LRS Bianchi type III Einstein–Klein–Gordon system introduced in [12], which is given by

$$\begin{equation}\dot {H}={H}^{2}\left[-\left(1+q\right)\right],\end{equation} \tag{ 2.6 }$$

$$\begin{equation}{\dot {{\Sigma}}}_{+}=H\left[-\left(2-q\right){{\Sigma}}_{+}+1-{{\Sigma}}_{+}^{2}-{\Omega}\right],\end{equation} \tag{ 2.7 }$$

$$\begin{equation}\ddot {\phi }+\phi =H\left[-3\dot {\phi }\right],\end{equation} \tag{ 2.8 }$$

with the deceleration parameter

$$\begin{equation*}q=2{{\Sigma}}_{+}^{2}+\frac{1}{6{H}^{2}}\left(2{\dot {\phi }}^{2}-{\phi }^{2}\right).\end{equation*}$$

Here H is the Hubble scalar, and we see that in this example the first component of (2.4) is the the Raychaudhuri equation (2.6). ϕ is the scalar field satisfying the Klein–Gordon equation (2.8). Σ₊ is the only independent component of the Hubble normalized shear tensor, and Ω is the Hubble normalized energy density. The Hamiltonian constraint $1-{{\Sigma}}_{+}^{2}-{\Omega}{ >}0$ and the non-vacuum condition Ω > 0 give the uniform bounds of Σ₊ ∈ (−1, 1) and Ω ∈ (0, 1). The Klein–Gordon equation (2.8) has the form of a Van der Pol equation. Via an amplitude-phase transformation (2.2) it can be formulated as the following two first-order equations [12, section 2.2]:

$$\begin{equation}\begin{aligned}\hfill \dot {{\Omega}}& =H\left[2{\Omega}\left(\right.1+q-3\enspace \mathrm{cos}{\left(t-\varphi \right)}^{2}\right],\hfill \\ \hfill \dot {\varphi }& =H\left[-3\enspace \mathrm{sin}\left(t-\varphi \right)\mathrm{cos}\left(t-\varphi \right)\right],\hfill \end{aligned}\end{equation} \tag{ 2.9 }$$

which also implies

$$\begin{equation*}q=2{{\Sigma}}_{+}^{2}+{\Omega}\left(3\enspace \mathrm{cos}{\left(t-\varphi \right)}^{2}-1\right).\end{equation*}$$

Then H and $\mathbf{x}\enspace \text{def}{=}\enspace {\left({{\Sigma}}_{+},{\Omega},\varphi \right)}^{\top }$ satisfy the system (2.4); i.e., the components of $\dot {\mathbf{x}}$ stand for both the remaining Einstein evolution equation (2.7) and the two first-order Klein–Gordon equation (2.9). Note that f¹ and f^[2] are smooth with respect to x. Also, both f¹ and f^[2] are in ${L}_{\mathbf{x},t}^{\infty }$.

Example 2. In addition to LRS Bianchi III also the Einstein Klein–Gordon systems of the spatially homogenous, LRS cosmologies of Bianchi types I, II, VIII and IX as well as Kantowski–Sachs can be brought into the form (2.4). The reader can verify that by starting from the formulation of [17], however using a Klein–Gordon field as matter model, and then performing the usual amplitude-phase transformation (2.2). For type IX it is however expected that the condition (2.5) fails, since it is a closed cosmology, and as such it tends to re-collapse; cf [25].

Example 3. The systems treated in [3, 5, 6] are also examples of the form (2.4). In [3] the matter source is a Yang–Mills field, which demonstrates that averaging techniques in mathematical cosmology are by no means limited to Einstein–Klein–Gordon cosmology.

Consider the system (2.4) with positive and strictly decreasing H = H(t), and with lim_t→∞  H(t) = 0. We then obtain the following large-time behaviour on the system (2.4):

Theorem 3.1 (Large-time behaviour). Suppose that H = H(t) > 0 is strictly decreasing in t, and

$$\begin{equation*}\underset{t\to \infty }{\mathrm{lim}}\enspace H\left(t\right)=0.\end{equation*}$$

Fix any ɛ > 0 with ɛ < H(0) and define t_* > 0 such that ɛ = H(t_*). Suppose that

$$\begin{equation*}{\Vert}{\mathbf{f}}^{1}{{\Vert}}_{{L}_{\mathbf{x},t}^{\infty }},\enspace {\Vert}{f}^{\left[2\right]}{{\Vert}}_{{L}_{\mathbf{x},t}^{\infty }}{< }\infty ,\end{equation*}$$

and that $\mathbf{f}$ ¹($\mathbf{x}$, t) is Lipschitz continuous and f^[2] is continuous with respect to $\mathbf{x}$ for all t ⩾ t_*. Also, assume that $\mathbf{f}$ ¹ and f^[2] are T-periodic for some T > 0. Then for all t > t_* with $t={t}_{{\ast}}+O\left(H{\left({t}_{{\ast}}\right)}^{-\gamma }\right)$ for any given γ ∈ (0, 1), we have

$$\begin{equation}\mathbf{x}\left(t\right)-\mathbf{z}\left(t\right)=O\left(H{\left({t}_{{\ast}}\right)}^{\mathrm{min}\left\{1,2-2\gamma \right\}}\right),\end{equation} \tag{ 3.1 }$$

where $\mathbf{x}$ is the solution of the system (2.4) with the initial condition $\mathbf{x}$(0) = $\mathbf{x}$ ₀ and $\mathbf{z}$(t) is the solution of the averaged equation

$$\begin{equation*}\dot {\mathbf{z}}=H\left({t}_{{\ast}}\right){\bar{\mathbf{f}}}^{1}\left(\mathbf{z}\right),\quad \mathrm{f}\mathrm{o}\mathrm{r}\enspace t{ >}{t}_{{\ast}}\end{equation*}$$

with the initial condition $\mathbf{z}$(t_*) = $\mathbf{x}$(t_*) where the average ${\bar{\mathbf{f}}}^{1}$ is defined as

$$\begin{equation*}{\bar{\mathbf{f}}}^{1}\left(\mathbf{z}\right)=\frac{1}{T}{\int }_{{t}_{{\ast}}}^{{t}_{{\ast}}+T}{\mathbf{f}}^{1}\left(\mathbf{z},s\right)\mathrm{d}s.\end{equation*}$$

Proof. Cf sections 4.1 and 4.2. □

Remark 2. Example 1 of section 2.4 satisfies the assumptions of Theorem 3.1. We expect the same to hold for examples 2 and 3, with the exception of LRS Bianchi IX which as a closed cosmology is expected to re-collapse and thus to fail the assumption (2.5); cf [25].

Remark 3. We emphasise the following key properties of theorem 3.1, which are connected to the choice of truncation time t_* to the first order system:

• The longer one waits with truncation, the better the estimate gets.
• Furthermore, the longer one waits with truncation, for the longer is the estimate valid.

These are attractive properties for the analysis of the large-time behaviour of such systems. To the best of our knowledge, theorem 3.1 is the first of its kind, which gives the large-time behaviour for a general class of systems, with a time dependent perturbation parameter, and which naturally appear in spatially homogenous cosmology with oscillating matter models.

In order to approximate the solution x(t) of the system (2.4), we consider the following truncated first-order system; for t > t_* with some given time t_*, we consider

$$\begin{equation}\left[\begin{matrix}\hfill \dot {\mathcal{H}}\hfill \\ \hfill \dot {\mathbf{y}}\hfill \end{matrix}\right]=\mathcal{H}\left(t\right)\left[\begin{matrix}\hfill 0\hfill \\ \hfill {\mathbf{f}}^{1}\left(\mathbf{y},t\right)\hfill \end{matrix}\right],\end{equation} \tag{ 4.1 }$$

where $\mathcal{H}\left({t}_{{\ast}}\right)=H\left({t}_{{\ast}}\right)$ and y(t_*) = x(t_*); i.e., the initial data coincide. Then since $\dot {\mathcal{H}}=0$, we have $\mathcal{H}\left(t\right)\equiv H\left({t}_{{\ast}}\right)$. We hereby denote ${H}_{{\ast}}\enspace \text{def}{=}\enspace H\left({t}_{{\ast}}\right)$. Then we introduce the following proposition on the error estimates for the first-order approximation:

Proposition 4.1. Fix any t_* > 0. Suppose that x and y are the solutions of (2.4) and (4.1), respectively, with the initial conditions x(0) = x₀ and y(t_*) = x(t_*). Suppose further that H = H(t) > 0 is strictly decreasing in t, and

$$\begin{equation*}\underset{t\to \infty }{\mathrm{lim}}\enspace H\left(t\right)=0.\end{equation*}$$

In addition, suppose that

$$\begin{equation*}{\Vert}{\mathbf{f}}^{1}{{\Vert}}_{{L}_{\mathbf{x},t}^{\infty }},\enspace {\Vert}{f}^{\left[2\right]}{{\Vert}}_{{L}_{\mathbf{x},t}^{\infty }}{< }\infty ,\end{equation*}$$

and that $\mathbf{f}$ ¹($\mathbf{x}$, t) is Lipschitz continuous in $\mathbf{x}$ for all t ⩾ t_*. Then for all t > t_* with $t={t}_{{\ast}}+O\left(H{\left({t}_{{\ast}}\right)}^{-\gamma }\right)$ for some γ ∈ (0, 1), we have

$$\begin{equation*}\mathbf{x}\left(t\right)-\mathbf{y}\left(t\right)=O\left(H{\left({t}_{{\ast}}\right)}^{2-2\gamma }\right).\end{equation*}$$

Proof. Since $\dot {\mathcal{H}}=0$ in (4.1), we have $\mathcal{H}\left(t\right)={H}_{{\ast}}$, for all t ⩾ t_*. Now, we subtract the equation for y in (4.1) from that for x in (2.4) to obtain that

$$\begin{equation*}\dot {\mathbf{x}}-\dot {\mathbf{y}}=H\left(t\right){\mathbf{f}}^{1}\left(\mathbf{x},t\right)-{H}_{{\ast}}{\mathbf{f}}^{1}\left(\mathbf{y},t\right).\end{equation*}$$

By integrating this identity, we have

$$\begin{equation}\mathbf{x}\left(t\right)-\mathbf{y}\left(t\right)={\int }_{{t}_{{\ast}}}^{t}\left(H\left(\tau \right){\mathbf{f}}^{1}\left(\mathbf{x},\tau \right)-{H}_{{\ast}}{\mathbf{f}}^{1}\left(\mathbf{y},\tau \right)\right)\mathrm{d}\tau ,\end{equation} \tag{ 4.2 }$$

by the initial condition y(t_*) = x(t_*). Now we consider the mean-value theorem for H(τ) and observe that, for some t_c = t_c (τ) ∈ (t_*, τ), we have

$$\begin{equation*}H\left(\tau \right)-{H}_{{\ast}}=\dot {H}\left({t}_{c}\right)\left(\tau -{t}_{{\ast}}\right),\enspace \quad \mathrm{f}\mathrm{o}\mathrm{r}\enspace \tau \in \left({t}_{{\ast}},t\right).\end{equation*}$$

Then (4.2) implies that

$$\begin{align*}\hfill \mathbf{x}\left(t\right)-\mathbf{y}\left(t\right)& ={\int }_{{t}_{{\ast}}}^{t}\left(\left({H}_{{\ast}}+\dot {H}\left({t}_{c}\right)\left(\tau -{t}_{{\ast}}\right)\right){\mathbf{f}}^{1}\left(\mathbf{x},\tau \right)-{H}_{{\ast}}{\mathbf{f}}^{1}\left(\mathbf{y},\tau \right)\right)\mathrm{d}\tau \hfill \\ \hfill & ={H}_{{\ast}}{\int }_{{t}_{{\ast}}}^{t}\left({\mathbf{f}}^{1}\left(\mathbf{x},\tau \right)-{\mathbf{f}}^{1}\left(\mathbf{y},\tau \right)\right)\mathrm{d}\tau +{\int }_{{t}_{{\ast}}}^{t}\dot {H}\left({t}_{c}\right)\left(\tau -{t}_{{\ast}}\right){\mathbf{f}}^{1}\left(\mathbf{x},\tau \right)\mathrm{d}\tau \hfill \\ \hfill & ={H}_{{\ast}}{\int }_{{t}_{{\ast}}}^{t}\left({\mathbf{f}}^{1}\left(\mathbf{x},\tau \right)-{\mathbf{f}}^{1}\left(\mathbf{y},\tau \right)\right)\mathrm{d}\tau \hfill \\ \hfill & \quad +{\int }_{{t}_{{\ast}}}^{t}{H}^{2}\left({t}_{c}\right){f}^{\left[2\right]}\left(\mathbf{x},{t}_{c}\right)\left(\tau -{t}_{{\ast}}\right){\mathbf{f}}^{1}\left(\mathbf{x},\tau \right)\mathrm{d}\tau ,\hfill \end{align*}$$

where the last line is by the first equation on $\dot {H}$ in the system (2.4). Therefore, we have

$$\begin{align*}\hfill & {\Vert}\mathbf{x}\left(t\right)-\mathbf{y}\left(t\right){\Vert}\hfill \\ \hfill & \qquad {\leqslant}{H}_{{\ast}}{\int }_{{t}_{{\ast}}}^{t}{\Vert}{\mathbf{f}}^{1}\left(\mathbf{x},\tau \right)-{\mathbf{f}}^{1}\left(\mathbf{y},\tau \right){\Vert}\mathrm{d}\tau \hfill \\ \hfill & \qquad +{\int }_{{t}_{{\ast}}}^{t}{H}^{2}\left({t}_{c}\right)\vert {f}^{\left[2\right]}\left(\mathbf{x},{t}_{c}\right)\vert \left(\tau -{t}_{{\ast}}\right){\Vert}{\mathbf{f}}^{1}\left(\mathbf{x},\tau \right){\Vert}\mathrm{d}\tau \hfill \\ \hfill & \qquad {\leqslant}{H}_{{\ast}}{\int }_{{t}_{{\ast}}}^{t}{\Vert}{\mathbf{f}}^{1}\left(\mathbf{x},\tau \right)-{\mathbf{f}}^{1}\left(\mathbf{y},\tau \right){\Vert}\mathrm{d}\tau +\frac{{\left(t-{t}_{{\ast}}\right)}^{2}}{2}{H}_{{\ast}}^{2}{\Vert}{f}^{\left[2\right]}{{\Vert}}_{{L}^{\infty }}{\Vert}{\mathbf{f}}^{1}{{\Vert}}_{{L}^{\infty }}\hfill \\ \hfill & \qquad {\leqslant}{H}_{{\ast}}{c}_{L}{\int }_{{t}_{{\ast}}}^{t}{\Vert}\mathbf{x}\left(\tau \right)-\mathbf{y}\left(\tau \right){\Vert}\mathrm{d}\tau +\frac{{\left(t-{t}_{{\ast}}\right)}^{2}}{2}{H}_{{\ast}}^{2}{\Vert}{f}^{\left[2\right]}{{\Vert}}_{{L}^{\infty }}{\Vert}{\mathbf{f}}^{1}{{\Vert}}_{{L}^{\infty }},\hfill \end{align*}$$

because H(⋅) is decreasing, where c_L is the Lipschitz constant and the third inequality is by the Lipschitz continuity of f¹. Now we apply the Grönwall inequality and obtain

$$\begin{equation*}{\Vert}\mathbf{x}\left(t\right)-\mathbf{y}\left(t\right){\Vert}{\leqslant}\frac{{\left(t-{t}_{{\ast}}\right)}^{2}}{2}{H}_{{\ast}}^{2}{\Vert}{f}^{\left[2\right]}{{\Vert}}_{{L}^{\infty }}{\Vert}{\mathbf{f}}^{1}{{\Vert}}_{{L}^{\infty }}\enspace \mathrm{exp}\left({H}_{{\ast}}{c}_{L}\left(t-{t}_{{\ast}}\right)\right).\end{equation*}$$

Therefore, for the time-scale $t={t}_{{\ast}}+O\left({H}_{{\ast}}^{-\gamma }\right)$ for any γ ∈ (0, 1) we have

$$\begin{equation*}{\Vert}\mathbf{x}\left(t\right)-\mathbf{y}\left(t\right){\Vert}{\leqslant}O\left({H}_{{\ast}}^{2-2\gamma }\right)\mathrm{exp}\left(O\left({H}_{{\ast}}^{1-\gamma }\right)\right)=O\left({H}_{{\ast}}^{2-2\gamma }\right),\end{equation*}$$

as long as H_* is finite. This completes the proof. □

Equipped with the first-order approximation and the error estimates, the next ingredient is a result from periodic averaging theory.

First we remind ourselves of the periodic averaging problem in standard form (2.3) which we discussed in section 2.2. For convenience we rewrite it here for the special case f^[2] = 0 as

$$\begin{equation}\dot {\mathbf{y}}=\varepsilon {\mathbf{f}}^{1}\left(\mathbf{y},t\right),\enspace \mathbf{y}\left(0\right)={\mathbf{y}}_{0},\quad \mathrm{f}\mathrm{o}\mathrm{r}\enspace t{ >}0\end{equation} \tag{ 4.3 }$$

with $\mathbf{y},\enspace {\mathbf{f}}^{1}\left(\mathbf{y},t\right)\in {\mathbb{R}}^{n}$, where f¹ is T-periodic and ɛ is a small parameter. Now we note that the equation for y of our first order system (4.1) is of this form with ɛ = H_*. Then we fix a sufficiently large time t_* > 0, and consider the first-order associated averaged system

$$\begin{equation}\dot {\mathbf{z}}=\varepsilon {\bar{\mathbf{f}}}^{1}\left(\mathbf{z}\right),\quad \mathrm{f}\mathrm{o}\mathrm{r}\enspace t{ >}{t}_{{\ast}}\end{equation} \tag{ 4.4 }$$

with the initial condition z(t_*) = y(t_*) where ɛ = H(t_*) and the average ${\bar{\mathbf{f}}}^{1}$ is defined as

$$\begin{equation*}{\bar{\mathbf{f}}}^{1}\left(\mathbf{z}\right)=\frac{1}{T}{\int }_{{t}_{{\ast}}}^{{t}_{{\ast}}+T}{\mathbf{f}}^{1}\left(\mathbf{z},s\right)\mathrm{d}s.\end{equation*}$$

Let $D\subset {\mathbb{R}}^{n}$ be a connected, bounded open set (with compact closure) containing the initial value y(t_*) = z(t_*). Fix L > 0 and ɛ₀ > 0 such that the solutions y(t, ɛ) and z(t, ɛ) with 0 ⩽ ɛ ⩽ ɛ₀ remain in D for t_* ⩽ t ⩽ t_* + L/ɛ. Then there is a basic result on the asymptotic error estimates:

Lemma 4.2 (theorem 2.8.1 of [22]). Let f¹ be Lipschitz continuous. Also, assume that f¹ is T-periodic for some T > 0. Let ɛ₀, D, L be as above. Then there exists a constant c > 0 such that the solutions y and z of (4.3) and (4.4) satisfy

$$\begin{equation*}{\Vert}\mathbf{y}\left(t,\varepsilon \right)-\mathbf{z}\left(t,\varepsilon \right){\Vert}{< }c\varepsilon ,\end{equation*}$$

for 0 ⩽ ɛ ⩽ ɛ₀ and ${t}_{{\ast}}{\leqslant}t{\leqslant}{t}_{{\ast}}+\frac{L}{\varepsilon }$, where || ⋅|| denotes the norm

$$\begin{equation*}{\Vert}\mathbf{u}{\Vert}\enspace \mathrm{def}{=}\enspace \sum\limits _{i=1}^{n}\vert {u}_{i}\vert \end{equation*}$$

for $\mathbf{u}\in {\mathbb{R}}^{n}$.

Then proposition 4.1 and lemma 4.2 together imply theorem 3.1 by letting ɛ = H(t_*), using the triangle inequality and noting that the time scale $O\left(H{\left({t}_{{\ast}}\right)}^{-\gamma }\right)$ is smaller than $O\left(H{\left({t}_{{\ast}}\right)}^{-1}\right)$ for a sufficiently large time t_* such that H(t_*) = ɛ ≪ 1. This completes the proof of theorem 3.1.