\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2016 (2016), No. 233, pp. 1--32.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2016 Texas State University.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2016/233\hfil 
Dynamical bifurcation in coupled oscillators]
{Dynamical bifurcation in a system of coupled oscillators with slowly varying
parameters}

\author[I. Parasyuk, B. Repeta \hfil EJDE-2016/233\hfilneg]
{Igor Parasyuk, Bogdan Repeta}

\address{Igor Parasyuk \newline
Faculty of Mechanics and Mathematics,
Taras Shevchenko National University of Kyiv,
64/13, Volodymyrska Street,
City of Kyiv, 01601, Ukraine}
\email{pio@univ.kiev.ua}

\address{Bogdan Repeta \newline
Faculty of Mechanics and Mathematics,
Taras Shevchenko National University of Kyiv,
64/13, Volodymyrska Street,
City of Kyiv, 01601, Ukraine}
\email{bogdan.repeta@gmail.com}

\thanks{Submitted May 2, 2016. Published August 25, 2016.}
\subjclass[2010]{34C15, 34C23, 34C46, 37D10}
\keywords{Dynamical bifurcations; transient processes; invariant tori;
\hfill\break\indent multi-frequency oscillations; coupled oscillators;
fast-slow systems}

\begin{abstract}
 This paper deals with a fast-slow system representing $n$ nonlinearly coupled
 oscillators with slowly varying parameters. We find conditions which
 guarantee that all $\omega$-limit sets near the slow surface of the system
 are equilibria and invariant tori of all dimensions not exceeding $n$,
 the tori of dimensions less then $n$ being hyperbolic. We show that a
 typical trajectory demonstrates the following transient process:
 while its slow component is far from the stationary points of the slow
 vector field, the fast component exhibits damping oscillations;
 afterwards, the former component enters and stays in a small neighborhood
 of some stationary point, and the oscillation amplitude of the latter begins
 to increase; eventually the trajectory is attracted by an $n$-dimesional
 invariant torus and a multi-frequency oscillatory regime is established.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{corollary}[theorem]{Corollary}
\allowdisplaybreaks


\section{Introduction}

The coupled oscillators theory plays a significant role in understanding
various patterns of collective behavior occurring in physical,
chemical, biological and social systems (see, e.~g., \cite{HopIzh97,ACN2016}
and references therein). The variety of behaviors exhibited by systems of
coupled oscillators (SCO) ranges from synchronization to complex chaotic motions.
In many cases, transient processes in SCO eventually turn into self-excited
 multi-frequency (quasiperiodic) oscillations on toric attractors.
Such a type of behavior in non-conservative systems was observed as early
as in the 20s--30s of the XX century and since that time was intensively
studied (see, e.~g.
 \cite{KryBog34, BMS69, EKST2013, Gilsinn,  Hale61, NayMoo79, SamMit76, 
SamPet04, Sam91,    Shch04}). In the middle of the XX century,
there was discovered a phenomenon of a 2-dimensional torus bifurcation
accompanying the stability loss of a limit cycle \cite{Kuz98, Nei59, Sack64}.
Later, studies on bifurcations of invariant tori and quasiperiodic
oscillations were conducted by many authors (see, e. g.,
\cite{Bib81,Bib90,Bib91,KKE2012, Han91,IooJos90,SchMar84} and references therein)
and the actual toolkit for qualitative investigation of such bifurcations
was developed in~\cite{ BMS69, Fen71, Hale61, SamMit76, Sack65,Sam91}.

The aforementioned results concern static bifurcation theory which
deals with systems dependent on time-constant parameters. Within the
framework of this theory, the birth of a stable $k$-dimensional invariant
torus from an equilibrium of the system
\begin{gather}
\dot{x}=f(x,u),\quad x\in\mathbb{R}^{d},\quad
(\dot{x}:=\frac{\mathrm{d}x}{\mathrm{d}t}),\label{eq:par-sys}
\end{gather}
dependent on the $m$-dimensional parameter $u$ can be ensured by the
following conditions: there exists a sufficiently smooth curve (equilibrium
curve) $x=\xi(s),u=\upsilon(s)$, $s\in(-1,1)$, such that $f(\xi(s),\upsilon(s))=0$
for all $s\in(-1,1)$; the eigenvalues of the Jacobi matrix
$f_{x}'(\xi(s),\upsilon(s))$
have negative real parts for all $s\in(-1,0)$ and positive real parts
for all $s\in(0,1)$; there exists a sufficiently smooth mapping
$\Xi(\cdot,\cdot)\colon\mathbb{T}^{k}\times[0,1)\to\mathbb{R}^n$
(here $\mathbb{T}^{k}:=\mathbb{R}^{k}/2\pi\mathbb{Z}^{k}$ denotes
the standard torus), such that $\Xi(\varphi,\upsilon(0))=\xi(0)$ for
all $\varphi\in\mathbb{T}^{k}$ and
$\operatorname{rank}\Xi_{\varphi}'(\varphi,\upsilon(s))=k$
for all $\varphi\in\mathbb{T}^{k}$, $s\in(0,1)$; finally, for any
$s\in(0,1)$ the toroidal surface $x=\Xi(\varphi,\upsilon(s))$ is
a local attractor of the flow generated by the system
$\dot{x}=f(x,\upsilon(s))$.
Under such conditions, when the parameter $u$, restricted to the
curve $u=\upsilon(s)$, passes through the point $u=\upsilon(0)$,
we observe the stability loss of equilibrium and the birth of a stable
invariant torus. It should be stressed that here the verb ``passes''
does not have any relation to a parameter motion over time.

On the contrary to the theory of static bifurcations, dynamical
bifurcation theory deals with systems which depend on slowly varying
in time parameters (fast-slow systems). Dynamical bifurcation
theory focuses on qualitative behavioral transformations which happen
in fast-slow systems due to the slow passage of parameters through
certain critical points in the parameter space. The origin of this theory
can be found in papers on relaxation oscillations (see the review
\cite{AAISh94}), although the term ``dynamical bifurcation''
appeared later, in the 80s of the XX century. The papers~\cite{Shi73,
Nei85, Nei87,Nei88, Ben91} gave start to studies of actually dynamic bifurcations
in fast-slow systems. During the last several decades many important
results concerning the considered area were obtained
\cite{Car84,IzhiHoppen,BuNeSch04,ShShT05,RachSch05,Ano05,KMR08,GKR08,Lie11}.
Some of the most peculiar features of fast-slow systems, such as the
delayed loss of stability, the synchronization, the existence of the
canard solutions and the blue-sky catastrophe can be of great importance
in the real-world applications. Nevertheless, some phenomena have
not yet been fully understood. In particular, as it was noted in \cite{AnoDV05},
this can be said about the emergence of multi-frequency oscillations
as a result of parameters evolution in fast-slow systems.

The present paper grounds on our previous results \cite{SamParRep,ParRep16}
and aims to fill the gap above. Here we consider the SCO governed by the
$n$-dimensional second order system
\begin{equation}\label{eq:osc-sys}
\ddot{w} + \Omega_0^2(u) w
= 2 \varepsilon \Lambda(u) \dot{w}+ F(w, \dot{w}, u, \mu),
\end{equation}
dependent on external (environmental) parameters
 $u = (u_1, \ldots, u_{m})$ and small positive parameters
$\varepsilon, \mu \ll 1$. Here
\begin{gather*}
\Lambda(u):= \operatorname{diag}(\lambda_1(u), \ldots, \lambda_{n}(u)), \quad
\lambda_{j}(\cdot) \in \mathrm{C}^{\infty} (\mathbb{R}^{m}; \mathbb{R}), \\
\Omega_0^2(u)
:= \operatorname{diag}(\omega_{01}^2(u), \ldots, \omega_{0n}^2(u)), \quad
 \omega_{0j}(\cdot)
\in \mathrm{C}^{\infty} (\mathbb{R}^{m}; \mathbb{R}_{++}),
 \\
 F(\cdot)
 \in \mathrm{C}^{\infty} (
 \mathbb{R}^{2n + m + 1}; \mathbb{R}^n)
\end{gather*}
and
\begin{equation}\label{eq:F=00003DO(w2p2)}
 F(w, p, u, 0)
  = O ( \| w \|^2 + \| p \|^2  ),
 \quad
 \| w \|^2 + \| p \|^2 \to 0,
\end{equation}
where $\| \cdot \| := \sqrt{\langle \cdot, \cdot \rangle}$ stands
for the Euclidean norm associated with the standard dot product in the 
coordinates vector space. Hereafter, we will also use a norm 
$| \cdot |$ defined as the sum of the absolute values of vector components.
We study the case where the slow evolution of parameters $u$ is governed 
by the system
\begin{equation}\label{eq:sys-evol-u}
 \dot{u}  = \mu G(w, \dot{w}, u, \mu)
\end{equation}
where $G(\cdot) \in \mathrm{C}^{ \infty} ( \mathbb{R}^{2n + m + 1}; \mathbb{R}^{m} )$. Thus, the parameter $\mu$ plays the role of the so-called slowness parameter.
When $\mu=0$, the system
\begin{equation}\label{eq:sys-osc-m=0}
 \ddot{w} + \Omega_0^2(u) w
  = 2 \varepsilon \Lambda(u) \dot{w}
  + F(w, \dot{w}, u, 0),\quad \dot{u}=0
\end{equation}
has an invariant surface given by the equations $w=\dot{w}=0$ and 
the parameter $\varepsilon$ is responsible for the oscillations damping
 rate (oscillations growth rate) of variables $w,\dot{w}$ near this surface 
for those $u$ that belong to the stability zone where all $\lambda_j(u)<0$ 
(complete instability zone where all $\lambda_j(u)>0$). 
The presence of two small parameters in fast-slow systems is a rather usual case.
 Initially these small parameters are completely independent, however, 
later we will impose a restriction that $\mu \propto \varepsilon$. 
All the functions involved in systems \eqref{eq:osc-sys} and \eqref{eq:sys-evol-u} 
may also continuously depend on $\varepsilon$, but we will not show this explicitly.

By the terminology of \cite{AAISh94}, the equations $w=\dot{w}=0$ define the 
so-called slow surface in the phase space $\mathbb{R}^{2n+m}$, the vector
 $g(u):=G(0,0,u,0)$ is called the slow velocity vector, and, in this way, 
we obtain the slow system on the slow surface
\begin{equation*}
 \dot{u}=g(u).
\end{equation*}
In \cite{IzhiHoppen,SamParRep,ParRep16}, there was considered the case when 
the slow system has a unique equilibrium attracting all its other trajectories. 
Here we study a more general situation allowing multiple equilibria, among 
which are stable, hyperbolic and completely unstable ones, but require the 
slow vector field to be gradient-like. This means that there exists a 
Morse function $V(\cdot) \in \mathrm{C}^2 (\mathbb{R}^{m} ; \mathbb{R}_{+} )$
 such that $\dot{V}_{g}(u) := \langle \nabla V(u), g(u) \rangle < 0$ 
for any non-stationary point $u$ of $V(\cdot)$. We find additional conditions 
under which a neighborhood of the slow surface is forward invariant under 
the semi-flow of system \eqref{eq:osc-sys}--\eqref{eq:sys-evol-u} and the 
set of all $\omega$-limit points contained in this neighborhood consists 
of equilibria and invariant tori of all dimensions less than or equal to $n$. 
We show that a typical forward semi-trajectory starting at $(w_0,\dot{w}_0,u_0)$,
 where $u_0$ belongs to the instability zone of the fast 
system \eqref{eq:sys-osc-m=0}, demonstrates the following transient process: 
while the slow component $u(t)$ is far from the stationary points of the Morse 
function $V$, the fast component $(w(t),\dot{w}(t))$ exhibits damping 
oscillations; afterwards, this component eventually enters and stays in a 
small neighborhood of some stationary point, and the oscillation amplitude 
of the fast component begins to increase. Since the trajectory is attracted 
by an invariant torus, eventually a multi-frequency oscillatory regime is 
established. Such behavior can be naturally interpreted as the dynamical 
bifurcation of multi-frequency oscillations.

In fact, we will also be able to categorize the solutions by their ultimate 
behavior near the slow surface. It will be shown that in a small neighborhood 
of the slow surface most of the system's trajectories, in terms of the
 Lebesgue measure, are attracted to trajectories on the stable $n$-dimensional 
invariant torus, while the rest ones lie on the stable manifolds of 
hyperbolic tori of dimension less than $n$.

The current article is organized as follows. 
Section \ref{sec:conditions} provides the key hypotheses regarding 
system \eqref{eq:osc-sys}--\eqref{eq:sys-evol-u} and the statement of 
the main theorem. Then, in section~\ref{sec:aux-lemmas} we introduce
auxiliary lemmas, which enable us to state in
section~\ref{sec:preliminary-results} certain preliminary results on 
the system's dynamics near its slow surface, and, consequently, 
describe the solutions behavior and classify them in 
sections~\ref{sec:ultimate1..n}, \ref{sec:ultimatei1..is}. 
After that, we provide an example depicting oscillation excitation 
in a circuit of two coupled oscillators which have components 
with temperature dependent properties. Finally, the paper ends with 
an addendum containing information on the normal forms method for 
systems with slowly varying parameters.


\section{Main theorem}
\label{sec:conditions}

Let us describe the conditions imposed on system 
\eqref{eq:osc-sys}--\eqref{eq:sys-evol-u}

We will assume that the slow vector field 
$u \mapsto g(u) := G(0, 0, u, 0)$ satisfy the following conditions
\begin{itemize}

\item[(H1)] 
 there exists a Morse function $V(\cdot) \in \mathrm{C}^2 (
 \mathbb{R}^{m}; \mathbb{R}_{+} )$ with the properties:
\begin{enumerate}
  \item   $V(\cdot)$ has a non-empty finite set of stationary points;

  \item   $V(u) \to +\infty$ when $\| u \| \to \infty$;

  \item   $\dot{V}_{g}(u) := \langle \nabla V(u), g(u) \rangle < 0$ for any
  non-stationary point $u$ of $V(\cdot)$;

  \item   the Hesse matrix $\frac{\partial^2}{\partial u^2} \dot{V}_{g}(u)$
  is negative definite at any stationary point of $V(\cdot)$.
 \end{enumerate}
\end{itemize}
Then, according to this hypothesis, any level set of $V(\cdot)$ is compact 
and if $V^{\ast} > 0$ is sufficiently large, then the sub-level set
\begin{equation*}
 \mathcal{V}  := V^{-1} ( [0, V^{\ast})  )
  = \{
 u \in \mathbb{R}^n \colon V(u) < V^{\ast}
  \}
\end{equation*}
contains the set
\begin{equation*}
 \mathcal{W} := \{
  u \in \mathbb{R}^{m}  \colon
  \nabla V(u) = 0
 \}
\end{equation*}
of all stationary points of $V(\cdot)$. Moreover, there exist such 
$\nu^{\ast} \geq \nu_{\ast} > 0$ and $\delta > 0$, that for any stationary 
point $u_{\ast}$ the following inequalities hold
\begin{gather}\label{eq:est_dotV_g}
  - \nu^{\ast} \| u - u_{\ast} \|^2 < \dot{V}_{g}(u) \leq - \nu_{\ast}
 \| u - u_{\ast} \|^2,
  \quad
  \| \nabla V(u) \|
   \leq \nu^{\ast} \| u - u_{\ast} \|
\end{gather}
 for all $u \colon \| u - u_{\ast} \| \leq \delta$,
and $\delta$-neighborhoods of any two points of $\mathcal{W}$ do not intersect. 
Obviously, $\mathcal{W}$ thereby coincides with the set of all singular points 
of the vector field $g$.

Now, for such a number $V^{\ast}$, let us adopt certain non-resonant
 hypotheses which are necessary for construction of the system's normal 
form in the whole domain $\mathcal{V}$, as well as in a vicinity of the set 
$\mathcal{W}$.
\begin{itemize}
 \item[(H2)]  if $u \in \operatorname{cl}(\mathcal{V})$, then
 \begin{equation*}
  \omega_{0i}(u) \neq \omega_{0j}(u),
  \quad
  \omega_{0i}(u) \neq \omega_{0j}(u) + \omega_{0k}(u),
  \quad
  \omega_{0i}(u) \neq \omega_{0j}(u) + \omega_{0k}(u) \pm \omega_{0l}(u)
 \end{equation*}
 for all $i, j, k, l \in \{ 1, 2, \ldots, n \}$.

 \item[(H3)]
 there exists such a number $N \geq 5$, that for any $u_{\ast} \in \mathcal{W}$ the
 equality
 \begin{equation*}
  \sum_{j = 1}^n (q_{j} - q_{j + n}) \omega_{0j}(u_{\ast})
 = \sigma \omega_{0i}(u_{\ast}),
 \end{equation*}
 where $\sigma \in \{ 0, 1 \}$, $i \in \{ 1,
 \ldots, n \}$, $\mathbf{q} = (q_1, \ldots, q_{2n}) \in
 \mathbb{Z}_{+}^{2n}$, $4 \leq \sum_{j = 1}^{2n} q_{j} \leq N$, is valid iff
 $q_i = q_{i + n} + \sigma$ and $q_{j} = q_{j + n}$ for all $j \in
 \{ 1, \ldots, n \} \setminus \{ i
 \}$.
\end{itemize}

Furthermore, we may assume that for all $\varepsilon \in [0, \varepsilon_0]$,
 with $\varepsilon_0 > 0$ being small enough, and for all
$u \in \operatorname{cl} (\mathcal{V})$ the  frequencies
\begin{equation*}
 \omega_{j}(u, \varepsilon)
  := \sqrt{\omega_{0j}^2(u)
 - (  \varepsilon\lambda_{j}(u) ) ^2}
  = \omega_{0j}(u) + O(\varepsilon^2),
 \quad
 j \in \{ 1, \ldots, n \}.
\end{equation*}
are correctly defined and satisfy hypotheses (H2), (H3) where each 
$\omega_{0j}(u)$ is replaced with $\omega_{j}(u, \varepsilon)$. Hereafter, 
to simplify our notations, we will omit explicit dependencies of functions 
on $\varepsilon$, as long as it does not lead to confusion. Thus, we will use 
the notation $\omega_{j}(u)$ instead of $\omega_{j}(u, \varepsilon)$ and so on.

Proceeding to the new variables $w, \dot{w} \mapsto x = (x_1, \ldots, x_{2n})$
by means of a substitution
\begin{equation}\label{eq:new-var-x}
 w_{j} = x_{2j - 1},
 \quad
 \dot{w}_{j}
  = \varepsilon\lambda_{j}(u) x_{2j - 1}
  + \omega_{j}(u) x_{2j},
 \quad
 j \in \{ 1, \ldots, n\},
\end{equation}
we come to an equivalent system
\begin{equation}\label{eq:main1}
 \dot{x} = J(u) x + \hat{F}(x, u, \mu),
 \quad
 \dot{u} = \mu \hat{G}(x, u, \mu),
\end{equation}
where
\begin{gather*}
 J(u)
  := \operatorname{diag}\Big[
 \begin{pmatrix}
  \varepsilon\lambda_1(u) & \omega_1(u)
  \\
  - \omega_1(u) & \varepsilon\lambda_1(u)
    \end{pmatrix},
 \ldots,
 \begin{pmatrix}
  \varepsilon\lambda_{n}(u) & \omega_{n}(u)
  \\
  - \omega_{n}(u) & \varepsilon\lambda_{n}(u)
 \end{pmatrix}
  \Big],
 \\
 \hat{G}(x, u, \mu)
  := 
 G(w, \dot{w}, u, \mu)
  \big|_{w, \dot{w} \mapsto x},
 \\
 \hat{F}(x, u, \mu)
  = (
 \hat{F}_1(x, u, \mu), \ldots, \hat{F}_{2n}(x, u, \mu)
  ),
 \\
 \hat{F}_{2j - 1}(x, u, \mu) \equiv 0,
 \\
\begin{aligned}
 \hat{F}_{2j}(x, u, \mu)
  :=& 
 \frac{1}{\omega_{j}(u)}
 F(w, \dot{w}, u, \mu)
  \big|_{w, \dot{w} \mapsto x}
  -\mu\varepsilon \sum_{i = 1}^{m}
  \frac{
 \partial\lambda_{j}(u)
  }{
 \partial u_i
  } \hat{G}_i(x, u, \mu) x_{2j - 1}\\
 &- \mu \sum_{i = 1}^{m}
  \frac{
 \partial\omega_{j}(u)
  }{
 \partial u_i
  } \hat{G}_i(x, u, \mu) x_{2j}.
\end{aligned}
\end{gather*}

In view of \eqref{eq:F=00003DO(w2p2)}, when $\mu = 0$, 
system~\eqref{eq:main1} has a slow invariant manifold of equilibria
 $\mathcal{M}_0$ defined by the equation $x = 0$.
Alike static bifurcation theory \cite{Bib90, Bib91}, we will study 
the behavior of system~\eqref{eq:main1} in a neighborhood of this manifold. 
And to do so, our first step will be finding conditions that guarantee 
the forward semi-invariance of such a neighborhood. 
This can be achieved by transforming the $N$-jet of system~\eqref{eq:main1} 
to the normal form with respect to variables $x$.

Let $\mathcal{D}_{N} \subset \mathbb{R}^{m}$ denote a
domain (or a collection of domains), such that for a fixed natural $N$ 
and for any $k \in \{ 1, \ldots, n \}$ and
$\sigma \in \{ 0, 1\}$ the equality
\begin{equation*}
 \min_{u \in \operatorname{cl} \mathcal{D}_{N}}
  \big|
 \sum_{l = 1}^n
  (q_{l} - q_{l + n} - \sigma\delta_{kl}) \omega_{0l}(u)
  \big| = 0,
\end{equation*}
where $\mathbf{q} \in \mathbb{Z}_{+}^{2n}, 1 \leq | \mathbf{q} | :=
\sum_{k = 1}^{2n} q_{k} \leq N$, fulfills iff
\begin{equation*}
 q_{l} - q_{l + n} - \sigma\delta_{kl} = 0
 \quad
 \forall l \in \{ 1, \ldots, n \}.
\end{equation*}

\begin{remark}\label{remark:proposition_homol_equ_1} \rm
Hypothesis (H2) guarantees that $\mathcal{D}_{3} = \mathcal{V}$, 
and under hypothesis (H3) the domain $\mathcal{D}_{N}$ is non-empty and 
contains some neighborhood of the set $\mathcal{W}$.
\end{remark}

If we introduce a vector of complex coordinates 
$z = (z_1, \ldots, z_{n}) \in \mathbb{C}^n$ and notations
\begin{equation*}
 \overrightarrow{| z |}
  := (
 | z_1 |,
 \ldots,
 | z | _{n}
  ),
 \quad
 \overrightarrow{| z | }^2
  := (
 | z_1 |^2,
 \ldots,
 | z | _{n}^2
 ),
 \quad
 \overrightarrow{| z |}^{\mathbf{p}}
  := \prod_{k = 1}^n | z_{k} |^{p_{k}},
\end{equation*}
where $\mathbf{p} := (p_1, \ldots, p_{n}) \in \mathbb{Z}_{+}^n$,
then, as it is shown in the Addendum, for all sufficiently small $\mu$ 
and $\varepsilon$ there exists a smoothly diffeomorphic change of 
variables $(x, u) \mapsto (\operatorname{Re} z, \operatorname{Im} z, v)$ which for $v \in \mathcal{D}_{N}$
transforms system~\eqref{eq:main1} into
\begin{equation}\label{eq:normform_2}
\begin{gathered}
\begin{aligned}
 \dot{z}_{k}
 &= \Big[
 \varepsilon\lambda_{k}(v)
 + \mathrm{i}\omega_{k}(v)
 + \sum_{3 \leq 2 | \mathbf{p} | + 1 \leq N}
  h_{0, k, \mathbf{p}}(v)
  \overrightarrow{| z |}^{2\mathbf{p}}
  \Big] z_{k}
  + \sum_{j = 1}^{P}
 \mu^{j}\eta_{j, k}(v)z_{k}
 \\
 &\quad +\sum_{j = 1}^{P}
   \mu^{j} \Big[
 \sum_{3 \leq 2 | \mathbf{p} | + 1 \leq N}
 h_{j, k, \mathbf{p}}(v)
 \overrightarrow{| z |}^{2\mathbf{p}}
   \Big] z_{k}
   + O (
 \| z \|^{N + 1} + \mu^{P + 1}
   ),
\end{aligned}  \\
  k \in \{ 1, \ldots, n \},
  \\
\begin{aligned}
 \dot{v}
 &= \mu \Big[
 g(v) + \sum_{2 \leq 2 | \mathbf{p} | \leq N}
  g_{0, \mathbf{p}}(v)
  \overrightarrow{| z |}^{2\mathbf{p}}
  + \sum_{j = 1}^{P}
   \mu^{j} \sum_{0 \leq 2 | \mathbf{p} | \leq N}
 g_{j, \mathbf{p}}(v)
 \overrightarrow{| z |}^{2\mathbf{p}}
 \Big]
 \\
&\quad +\mu O (
  \| z \|^{N + 1} + \mu^{P + 1}
 ).
\end{aligned}
\end{gathered}
\end{equation}
Here $P \geq N/2$ is an arbitrary fixed natural number, 
$\eta_{j, k}(\cdot), h_{j, k, \mathbf{p}}(\cdot)$ are smooth complex-valued 
functions in $\mathcal{D}_{N}$ and $g_{j, \mathbf{p}}(\cdot)$ are 
smooth $\mathbb{C}^n$-valued functions in $\mathcal{D}_{N}$. 
Besides that, all these functions smoothly depend on the parameter 
$\varepsilon$. The remainder terms are smooth in the sense of real calculus 
on the set
\begin{equation*}
  \| z \| < \delta,
  \quad
  v \in \mathcal{D}_{N},
  \quad
  \mu \in [0, \mu_0],
  \quad
  \varepsilon \in [0, \varepsilon_0]
\end{equation*}
with sufficiently small positive numbers $\delta$, $\mu_0$ and $\varepsilon_0$,
and are uniform with respect to $v \in \mathcal{D}_{N}$ and 
$\varepsilon \in [0, \varepsilon_0]$.

Further, we will also denote
\begin{gather*}
 \lambda(v) = (
  \lambda_1(v), \ldots, \lambda_{n}(v)
 ),
 \quad
 A(v) = \{ a_{kl}(v)\}_{k, l = 1}^n,
 \\
 a_{kl}(v) := - \mathrm{Re}\, h_{0, k, \mathbf{e}_{l}}(v),
 \quad
 b_{kl}(v) := - \mathrm{Im}\, h_{0, k, \mathbf{e}_{l}}(v),
\end{gather*}
where $\mathbf{e}_{l} \in \mathbb{Z}_{+}^n$ is a vector having its $l$-th
component equal $1$ and all other equal $0$. As we will see later, 
the functions $\lambda(v)$ and $A(v)$ play the key role in emergence of 
the bifurcation phenomenon, which is why we require them to satisfy 
additional constraints.
\begin{itemize}

\item[(H4)]
 the symmetrical part of the matrix $A(v)$ is positive definite on $\operatorname{cl}
 (\mathcal{V})$, and all non-diagonal elements $a_{ij}(v)$, $i \neq j$, are
 non-positive at any stationary point $v \in \mathcal{W}$.

\item[(H5)]
 the set $\mathcal{V}$ admits representation as a union
of three nonempty sets
\begin{gather*}
 \mathcal{V}_{+} := \{
  v \in \mathcal{V}
  \colon
  \lambda_{j}(v) > 0
  \quad
  \forall j \in \{1, \ldots, n \}
 \},
 \\
 \mathcal{V}_{-} := \{
  v \in \mathcal{V}
  \colon
  \lambda_{j}(v) < 0
  \quad
  \forall j \in \{1, \ldots, n \}
 \},
 \\
 \mathcal{V_{\ast} = \mathcal{V}\setminus} \left[
  \mathcal{V}_{+} \cup \mathcal{V}_{-}
 \right],
\end{gather*}
and each function $\lambda_{j}(\cdot)$, $j\in\{ 1,\ldots,n\} $,
is positive at any stationary point of $V(\cdot)$.
\end{itemize}

Note, that $\mathcal{W} \subset \mathcal{V}_{+}$ and for system~\eqref{eq:main1} 
with $\mu = 0$ the submanifold $\{(x,u)\in\mathcal{M}_0:u \in \mathcal{V}_{-}\}$ 
is a local attractor, while $\{(x,u)\in\mathcal{M}_0:u \in \mathcal{V}_{-}\}$ 
is a local repeller.

Let us fix sufficiently small $\kappa>0$ in such a way that
\begin{equation*}
\mathcal{V}_{-}^{\kappa}:=\{ v\in\mathcal{V}_{-}:
\lambda_{j}(v)\le-\kappa\quad\forall j\in\{ 1,\ldots,n\} \} \ne\emptyset.
\end{equation*}

Now we are in a position to state our main result.

\begin{theorem}\label{thm:main}
There exist $\rho_0>0$, $\varsigma_0>0$, and for any
$\varsigma_{\ast}\in(0,\varsigma_0)$
there is $\varepsilon_0>0$, such that once $\varepsilon\in(0,\varepsilon_0)$
and $\mu\in(\varsigma_{\ast}\varepsilon,\varsigma_0\varepsilon)$,
the following statements are true:
\begin{enumerate}
\item system~\eqref{eq:normform_2} generates a forward semi-flow on the
set
\begin{equation*}
 \mathfrak{A}
  := \{
 (z,v) \in \mathbb{C}^n \times \mathbb{R}^{m}
 \colon
 \| z\| \le \rho_0,
 v\in\mathcal{V}
  \};
\end{equation*}

\item if $\mathbb{R}_{+}\ni t \mapsto (z(t),u(t))$ is a solution of system~\eqref{eq:normform_2},
such that $(z(0),v(0))\in\mathfrak{A}$ and $(t_0,t_1)\in\mathbb{R}_{+}$
is an interval, such that $\| z(t)\| >\mu^{P}$ and
$v(t)\in\mathcal{V}_{-}^{\kappa}$ for all $t\in(t_0,t_1)$, then
\begin{gather*}
\| z(t)\| \le\mathrm{e}^{-\varepsilon\kappa(t-t_0)/2}\| z(t_0)\| \quad\forall t\in(t_0,t_1);
\end{gather*}

\item to any point of the set $\mathcal{W}$, one can put into correspondence
a finite collection of invariant tori belonging to $\mathfrak{A}$,
and each such collection contains tori of all dimensions from $1$
to $n$; in addition, any torus of a dimension less then $n$ is truly
hyperbolic, while any $n$-dimensional torus corresponding to a local
minimum of the Morse function $V(\cdot)$ is a local attractor of
system~\eqref{eq:normform_2};

\item any non-equilibrium forward semi-trajectory of this system lying in
$\mathfrak{A}$ is attracted by one of the invariant tori, and those
trajectories that are attracted by $n$-dimensional tori form the
set of the full Lebesgue measure in $\mathfrak{A}$;

\item each forward semi-trajectory approaching the $n$-dimensional invariant
torus is attracted by a forward semi-trajectory lying on this torus.

\end{enumerate}
\end{theorem}

The rest of this article will be devoted to proving and illustrating this theorem.


\section{Auxiliary lemmas}
\label{sec:aux-lemmas}

\begin{lemma}\label{lemma:W_f < 0}
For $\delta > 0$, set 
$\bar{B}_{\delta}^{d} := \{ x \in \mathbb{R}^{d} \colon \| x \| \leq \delta\}$. 
Suppose that there exists a Morse function
 $W(\cdot) \in \mathrm{C}^2 ( \bar{B}_{\delta}^{d}; \mathbb{R} )$ 
having a unique stationary point $x_{\ast} = 0$ and a vector field 
$f(\cdot) \in \mathrm{C}^{1} ( \bar{B}_{\delta}^{d} ; \mathbb{R}^{d} )$, such that
\begin{equation*}
  \langle \nabla W(x), f(x) \rangle
   \leq -\theta \| x \| ^2,
  \quad
  \| f'(x) \|
   \leq \theta L
  \quad
  \forall x \in \bar{B}_{\delta}^{d}(0),
\end{equation*}
where $L$ and $\theta$ are some positive constants. Define
\begin{equation*}
 K := \max_{\| x \| \leq \delta}
  \| W''(x) \|
\end{equation*}
and let $M$ and $\epsilon$ be arbitrary positive numbers satisfying
\begin{equation}\label{eq:M^0}
 M \geq M_{\ast}(K, L)
  := \frac{
 1 + KL + \sqrt{1 + 2K^{3}L^{3}}
  }{
 L
  },
 \quad
 M \epsilon < \delta.
\end{equation}
Then for any $f_1(\cdot) \in \mathrm{C}^{1} ( [\tau_0, \infty)
\times \bar{B}_{\delta}^{d} ; \bar{B}_{\theta}^{d} )$ and any 
$x_0 \in B_{\delta}^{d}$ a solution $x(t)$ of the initial problem
\begin{equation}\label{eq:f + f_1}
 \dot{x} = f(x) + \varepsilon f_1(t, x),
 \quad
 x(\tau_0) = x_0
\end{equation}
meets the following alternative: either there exists such
 $\tau^{\ast} > \tau_0$, that $\| x ( \tau^{\ast} ) \| = \delta$, or
there exists such $\tau_{\ast} \geq \tau_0$, that $\| x(t) \| < M\epsilon$
for all $t > \tau_{\ast}$.

Additionally, if $0 < \epsilon < \delta / N(K, L, M)$, where
\begin{equation}\label{eq:d^0}
 N(K, L, M) := \frac{
  1 + KL + \sqrt{\left[ 1 + (K - M)L \right]^2 + 2KL^{3}M^2}
 }{
  L
 } > M,
\end{equation}
and if the first scenario takes place, but at some instant of time the 
solution belongs to $\bar{B}_{M \epsilon}^{d}$, then
\begin{equation*}
 W ( x(\tau^{\ast}) )
  < W(0) - \frac{
 K(M\epsilon)^2
  }{
 2
  } \leq \min
 \{ W(x) \colon \| x \| = M\epsilon\}.
\end{equation*}

Furthermore, in the case when the stationary point $x_{\ast} = 0$ is elliptic, numbers $\beta \in (0, \delta)$ and $\epsilon > 0$ are such, that
\begin{equation}\label{eq:maxW < minW}
 \max
  \{ W(x) \colon \| x \| = \beta\}
 < \min
  \{ W(x) \colon \| x \| = \delta\},
 \quad
 0 < \epsilon < \beta / K,
\end{equation}
and the solution of~\eqref{eq:f + f_1} at some moment of time belongs to $\bar{B}_{\beta}^{d}$, then the second scenario fulfills.
\end{lemma}

\begin{proof}
Since $\nabla W(0) = 0$, it is true that $\| \nabla W(x) \| \leq K \| x \|$
and
\begin{equation*}
 \langle \nabla W(x), f(x) + \varepsilon f_1(t, x) \rangle
  \leq \theta \left[
 - \| x \| + \epsilon K
  \right] \| x \|
 \quad
 \forall x \in \bar{B}_{\delta}^{d}.
\end{equation*}
Besides,
\begin{equation*}
 \max
  \{ W(x) \colon \| x \| = \varrho\}
  \leq W(0) + \frac{K\varrho^2}{2},
 \quad
 \min
  \{ W(x) \colon \| x \| = \varrho\}
  \geq W(0) - \frac{K\varrho^2}{2}
\end{equation*}
for any $\varrho \in (0, \delta]$. As the Hessian of $W(\cdot)$ at $x = 0$ 
is non-degenerate, we have $f(0) = 0$ and $\| f(x) \| \leq \theta L \| x \|$. Hence,
\begin{equation*}
 \| f(x) + \varepsilon f_1(t, x) \|
  \leq \theta (
 L \| x \| + \epsilon
  )
 \quad
 \forall x \in \bar{B}_{\delta}^{d}.
\end{equation*}

Now let us demonstrate how to choose $M$. At first, we will require only that 
$M \geq K$ and $M\epsilon < \delta$. Take an arbitrary 
$\varrho \in ( K\epsilon, M\epsilon )$. If the moment $\tau^{\ast}$ does not exists, 
i.~e. $\| x(t) \| < \delta$ for all $t \geq \tau_0$, then the function
$W(x(t))$ strictly decreases
until $x(t)$ reaches the sphere $\| x \| = \varrho$ at an instant of time 
$\tau_{\ast} \geq \tau_0$. The moment $\tau_{\ast}$ necessarily exists,
since otherwise $W(x(t))$ would decrease unboundedly, which is impossible.

Suppose that $x(t)$ reaches the sphere $\| x \| = M\epsilon$ after the moment 
$\tau_{\ast}$. Then there exist $\tau_{2} > \tau_{\ast}$ and 
$\tau_1 \in [\tau_{\ast}, \tau_{2})$, such that
\begin{equation*}
 \| x(\tau_1) \| = \varrho,
 \quad
 \varrho< \| x(t) \| < M \epsilon
 \quad
 \forall t \in (\tau_1, \tau_{2}),
 \quad
 \| x(\tau_{2}) \| = M\epsilon.
\end{equation*}
For $t \in [\tau_1, \tau_{2}]$, we have
\begin{equation*}
 \frac{
  \mathrm{d} \| x(t) \|
 }{
  \mathrm{d}t
 } \leq 
  \frac{
 \langle x, f(x) + \varepsilon f_1(t, x) \rangle
  }{
 \| x \|
  } \big|_{x = x(t)}
 \leq \theta (
  L \| x(t) \| + \epsilon
 ),
\end{equation*}
which implies
\begin{equation*}
 \frac{
  \mathrm{d} \| x(t) \| / \mathrm{d}t
 }{
  \theta (
 L \| x(t) \| + \epsilon
  )
 } \leq 1
\end{equation*}
and
\begin{align*}
W(x(\tau_{2})) - W(x(\tau_1))
& = \int_{\tau_1}^{\tau_{2}}
 \langle
   \nabla W(x), f(x) + \varepsilon f_1(t, x)
  \rangle
 \big|_{x = x(t)}
 \mathrm{d} t
 \\
&\leq\int_{\tau_1}^{\tau_{2}}
   \theta \left[
 - \| x(t) \| + K\epsilon
   \right] \| x(t) \| \frac{
 \mathrm{d} \| x(t) \| / \mathrm{d}t
   }{
 \theta (
  L \| x(t) \| + \epsilon
 )
   } \mathrm{d} t \\
&= \int_{\varrho}^{M\epsilon}
   \frac{
 (
  - s + K\epsilon
 ) s
   }{
 Ls + \epsilon
   } \mathrm{d} s.
\end{align*}
Taking into account that $W(x(\tau_1)) \leq W(0) + K\varrho^2 / 2$
and making $\varrho$ tend to $K\epsilon$, we obtain the estimate
\begin{align*}
 W(x(\tau_{2}))
&\leq W(0) + \frac{
 K^{3} \epsilon^2
  }{
 2
  } + \big[
 - \frac{s^2}{2L}
 + \frac{\epsilon(1 + KL)s}{L^2}
 - \frac{
  \epsilon^2(1 + KL)
 }{
  L^{3}
 } \ln \frac{
  Ls + \epsilon
 }{
  L
 }
  \big]_{K\epsilon}^{M\epsilon}
 \\
&< W(0)
   + \frac{K^{3}\epsilon^2}{2}
   - \frac{\epsilon^2}{2L} \big[
 M^2 - \frac{2(1 + KL)M}{L}
  + \frac{2(1 + KL)K}{L} - K^2
 \big].
\end{align*}

If we introduce a quadratic polynomial of variable $\xi$,
\begin{equation*}
 \mathcal{P}(\xi; \epsilon, \eta)
  := \xi^2
  - \frac{2(1 + KL)\epsilon}{L} \xi
  + \frac{2(1 + KL)\eta\epsilon^2}{L}
  - ( 1 + 2KL ) \eta^2\epsilon^2,
\end{equation*}
where $\epsilon$ and $\eta$ are positive parameters, 
one can verify that $M_{\ast}(K, L)$
is the greatest root of $\mathcal{P}(\xi; 1, K)$. 
Thus, for any $M \geq M_{\ast}(K, L)$
and any $\epsilon < \delta / M$, we obtain
$\mathcal{P}(M \varepsilon;\epsilon, K) = \varepsilon^2 \mathcal{P}(M; 1, K) \geq 0$. 
It means that
\begin{gather*}
 M^2 - \frac{2(1 + KL)M}{L} + \frac{2(1 + KL)K}{L} - K^2 \geq 2 K^{3} L, \\
 W(x(\tau_{2}))
  < W(0) - \frac{K^{3}\epsilon^2}{2}
  \leq \min
 \{
  W(x)
  \colon
  \| x \| = K\epsilon
 \} .
\end{gather*}
Hence, the function $W(x(t))$ keeps on decreasing and satisfies for 
$t > \tau_{2}$ the inequality $W(x(t)) < W(x(\tau_{2}))$. 
This means that $x(t)$ never reaches the sphere $\| x \| = K\epsilon$, and moreover,
\begin{equation*}
 \inf_{t \geq \tau_{2}}
  \left[
   \| x(t) \| - K\epsilon
  \right] > 0.
\end{equation*}
As a consequence, $W(x(t)) \to -\infty$ as $t \to \infty$, and we come to a 
contradiction. Therefore, such a choice of $M$ and $\epsilon$ guarantees 
the validity of the inequality $\| x(t) \| < M\epsilon$ for all $t > \tau_{\ast}$.

In a similar way, one can show that if $\| x(\tau') \| \leq M\epsilon$,
but there exists $\tau^{\ast} > \tau'$, such that $\| x(\tau^{\ast}) \| = \delta$, 
then
\begin{align*}
& W ( x(\tau^{\ast}) ) \\
&\leq W(0) + \frac{KM^2\epsilon^2}{2}
  - \big[
 \frac{s^2}{2L}
 - \frac{\epsilon(1 + KL)s}{L^2}
 \frac{\epsilon^2(1 + KL)}{L^{3}} \ln \frac{Ls + \epsilon}{L}
  \big]_{M\epsilon}^{\delta}
 \\
&\leq W(0)
  + \frac{KM^2\epsilon^2}{2}
  - \frac{1}{2L} \big[
 \delta^2
 - \frac{2(1 + KL)\epsilon}{L}\delta
 + \frac{2(1 + KL)M\epsilon^2}{L}
 - (M\epsilon)^2
  \big].
\end{align*}
Since $N(K, L, M)$ is the greatest root of $\mathcal{P}(\xi; 1, M)$, 
once $0 < N(K, L, M)\epsilon < \delta$, we have 
$\mathcal{P}(\delta; \epsilon, M) \geq 0$ and
\begin{equation*}
 W ( x(\tau^{\ast}) )
  < W(0) - \frac{K ( M\epsilon )^2}{2}
  \leq \min \{
 W(x)
 \colon
 \| x \| = M\epsilon
  \}.
\end{equation*}

Finally, if the point $0$ is elliptic and inequalities~\eqref{eq:maxW < minW} are
fulfill, then
\begin{equation}\label{eq:W' < 0}
 \langle \nabla W(x), f(x) + \varepsilon f_1(t, x) \rangle
  \leq \left[
 -\theta \| x \| + \epsilon K
  \right] \| x \|
  < 0,
\end{equation}
as soon as $\beta \leq \| x \| \leq \delta$. Let the solution belong to 
$\bar{B}_{\beta}^{d}$ at some moment of time. If by reasoning ad absurdum
 we supposed that there existed such $\tau^{\ast} > \tau_0$,
that $\| x(\tau^{\ast}) \| = \delta$, then there would exist $\tau'' < \tau^{\ast}$, 
such that
\begin{equation*}
 \| x(\tau'') \| = \beta,
 \quad
 \beta < \| x(t) \|
  < \delta
 \quad
 \forall t \in (\tau'', \tau^{\ast}).
\end{equation*}
Thereby, $W ( x(\tau'') ) \leq \max_{ \| x \| = \beta} W(x) < \min_{ \| x \| 
= \delta} W ( x(\tau^{\ast}) )$, which is impossible, 
since~\eqref{eq:W' < 0} yields that the function $W(x(t))$ is decreasing 
on $(\tau', \tau^{\ast})$.
\end{proof}

\begin{lemma}\label{lemma:W'f < 0 - 2}
Let $D \subset \mathbb{R}^{d}$ be a bounded domain with a $\mathrm{C}^2$-boundary, 
$\bar{D} := \operatorname{cl} (D)$ and let $W(\cdot) \in \mathrm{C}^2 ( \bar{D} ; \mathbb{R} )$
be a Morse function with a finite set of stationary points $\mathfrak{W} \subset D$. 
Define
\begin{equation*}
 K := \max_{x \in \bar{D}}
  \{
 \| \nabla W(x) \|,
 \| W''(x) \|
  \}
\end{equation*}
and choose sufficiently small $\delta > 0$ and $\beta \in (0, \delta)$ 
that meet the following requirements:
\begin{enumerate}
 \item
 the $\delta$-neighborhood of \textup{$\mathfrak{W}$} belongs to $D$;

 \item
 for any $x_{\ast}', x_{\ast}'' \in \mathfrak{W}$,
 such that $W(x_{\ast}') > W(x_{\ast}'')$ the following
 inequality holds:
 \begin{equation*}
  \min \{
 W(x)
 \colon
 \| x - x_{\ast}' \| = \delta
  \}
  > \max \{
 W(x)
 \colon
 \| x - x_{\ast}'' \| = \delta
  \};
 \end{equation*}

 \item
 for any elliptic point $x_{\ast} \in \mathfrak{W}$ there holds the inequality
 \begin{equation*}
  \max \{
 W(x)
 \colon
 \| x - x_{\ast} \| = \beta
  \}
  < \min \{
 W(x)
 \colon
 \| x - x_{\ast} \| = \delta
  \}.
\end{equation*}
\end{enumerate}

Also, suppose that there is such a vector field 
$f(\cdot) \in \mathrm{C}^{1} ( \bar{D}; \mathbb{R}^{d} )$, 
that for some positive constants $\theta, L$ such inequalities fulfill:
\begin{gather*}
 \langle \nabla W(x), f(x) \rangle
  < -\theta\delta^2
 \quad
 \forall x \in D
  \colon
  \min_{y \in \mathfrak{W}}
 \| x - y \| > \delta,
 \\
 \langle \nabla W(x), f(x) \rangle
  \leq -\theta \| x - y \| ^2,
 \quad
 \| f'(x) \|
  \leq \theta L
 \quad
 \forall y \in \mathfrak{W},
  \forall x \colon \| x - y \| \leq \delta.
\end{gather*}
Then, with the corresponding functions defined by formulae~\eqref{eq:M^0} 
and~\eqref{eq:d^0}, for any $M \geq M_{\ast}(K, L)$ and all 
$\epsilon \in ( 0, \epsilon_0(K, L, M) )$, where
\begin{equation*}
 \epsilon_0(K, L, M)
  := \min \{
 \frac{\delta^2}{K}, \frac{\beta}{N(K, L, M)}
  \},
\end{equation*}
the following assertion is correct. 
If $f_1(\cdot) \in \mathrm{C}^{1} ( [\tau_0, \infty) \times \bar{D} ;
\bar{B}_{\theta}^{d} )$, then for any such $x_0 \in D$, that the corresponding
solution $x(\cdot)$ of initial problem~\eqref{eq:f + f_1} is defined on 
$[\tau_0, \infty)$ and takes values in $D$, there exist
$x_{\ast} \in \mathfrak{W}$ and
$t_{\ast} > \tau_0$, such that $\| x(t) - x_{\ast} \| < M\epsilon$
for all $t > t_{\ast}$.
\end{lemma}

\begin{proof}
Under the conditions of this lemma, we have
\begin{equation}\label{eq:nablaW, f + f1 < 0}
 \langle \nabla W(x), f(x) + \varepsilon f_1(t, x) \rangle < 0
 \quad
 \forall x \colon \min_{y \in \mathfrak{W}}
  \| x - y \| \geq M\epsilon.
\end{equation}
Therefore, if the solution $x(\cdot)$ is defined on $[\tau_0, \infty)$
and takes values in $D$, then there exist $x_{\ast}^{1} \in \mathfrak{W}$ 
and $t_1 \geq \tau_0$, such that $\| x(t_1) - x_{\ast}^{1} \| < M\epsilon$.
Indeed, otherwise, the function $W(x(t))$ would decrease unboundedly when 
$t \to \infty$, which is impossible.

By Lemma~\ref{lemma:W_f < 0}, the solution $x(\cdot)$ meets the following
alternative: either for all $t \geq t_1$ we have 
$\| x(t) - x_{\ast} \| \leq M\epsilon$, or there exists $t_{2} > t_1$, 
such that $\| x(t_{2}) - x_{\ast}^{1} \| = \delta$ and
\begin{equation*}
 W ( x(t_{2}) )
  < W(x_{\ast}^{1}) - \frac{K(M\epsilon)^2}{2}
  \leq \min \{
 W(x)
 \colon
 \| x - x_{\ast}^{1} \| = M\epsilon
  \}.
\end{equation*}
The first case always takes place if $x_{\ast}^{1}$ is elliptic. In
the second one, on account of choice of $\delta$ and~\eqref{eq:nablaW, f + f1 < 0}, 
there exist such $x_{\ast}^2 \in \mathfrak{W}$ and $t_{3} > t_{2}$, that 
$\| x(t_{3}) - x_{\ast}^2 \| < M\epsilon$, and $W(x_{\ast}^2) < W(x_{\ast}^{1})$. 
Now, it is clear that eventually the solution enters and then never leaves an 
$M\epsilon$-neighborhood of some point $x_{\ast} \in \mathfrak{W}$.
\end{proof}

For the sake of completeness of our exposition, let us represent the following
 result from the theory of non-negative invertible matrices.

\begin{lemma}\label{lemma:A inv - pos}
Let a real matrix $P = \{ p_{ij} \}_{i, j = 1}^{d}$ have such properties:
\begin{enumerate}
 \item \label{lst:PisPosDef}
 the matrix $P + P^{\top}$ is positive definite;

 \item \label{lst:Pij>0}
 $p_{ij} \leq 0$ for any $i, j \in \{ 1, \ldots, d\} $, $i \neq j$.
\end{enumerate}
Then for any vector $y = (y_1, \ldots, y_{d})$ with the positive elements 
each component of the vector $x := P^{-1}y$ satisfies the inequalities
\begin{gather}
 \label{eq:x - k > y - k / p - k}
 x_i \geq \frac{y_i}{p_{ii}},
 \quad
 i \in \{ 1, \ldots, d\},
 \\
 \label{eq:|x| < }
 | x |
  \leq \frac{
 \max_{1 \leq i \leq d}y_i
  }{
 p_{+}
  },
\end{gather}
where $p_{+} := \min \{ \langle P\xi, \xi \rangle :\xi \in \mathbb{R}_{+}^{d},
 | \xi | = 1, \}$
\end{lemma}

\begin{proof}
The system $\dot{x} = - Px$, in which $P$ has the property~\ref{lst:PisPosDef} 
is asymptotically stable. Hence, all eigenvalues of $-P$ have negative real parts. 
By \cite[Theorem 5, XIII]{Gantmacher}, all of the main minors of $P$ are positive. 
Thus, according to \cite[Theorem 6.3]{Nikaido} the matrix $P^{-1}$ has
 non-negative elements, and any row of this matrix contains at least one 
positive element. Obviously, then $x := P^{-1}y \in \mathbb{R}_{++}^{d}$, 
once $y \in \mathbb{R}_{++}^{d}$. Moreover, since
\begin{equation*}
 \sum_{j = 1, j \neq i}^{d}
  p_{ij}x_{j} \geq 0,
\end{equation*}
the components $x_i$ satisfy~\eqref{eq:x - k > y - k / p - k},
and the inequalities
\begin{equation*}
 p_{+} | x |^2
  \leq \langle Px, x \rangle
  = \langle y, x \rangle
  \leq \max_{1 \leq i \leq d}y_i | x |
\end{equation*}
yield~\eqref{eq:|x| < }.
\end{proof}


\section{Preliminary results on the behavior of the normalized system}
\label{sec:preliminary-results}

Having introduced the aforementioned lemmas, we may proceed to the
investigation of system~\eqref{eq:normform_2} dynamics. This section 
provides the general description of the solutions behavior and suggests 
a way to classify them. Later we will refine this information.
Define
\begin{equation*}
\begin{gathered}
 a_{+} := \min \{
  \langle A(v) r, r \rangle
  \colon
  \forall r \in \mathbb{R}_{+}^n,
  | r | = 1,
  v \in \operatorname{cl} (\mathcal{V})
 \},
 \\
 \lambda^{+} := \max \{
  \langle \lambda(v), r \rangle
  \colon
  \forall r \in \mathbb{R}_{+}^n,
  | r | = 1,
  v \in \operatorname{cl} (\mathcal{V})
 \}.
\end{gathered}
\end{equation*}

\begin{proposition}\label{proposition:sqrt - mu - neighbor}
Assume hypotheses {\rm (H1), (H2), (H4)}  to be true and $N = 3$. 
Then there exist a sufficiently small $\rho_0 > 0$ and sufficiently large
$\rho_{\ast} > 0$, $R_{\ast} > 0$, such that for any $\rho > \rho_{\ast}$,
 $R > R_{\ast}$ one can choose $\varepsilon_0 > 0$ in such a way, that
once $\varepsilon \in (0, \varepsilon_0]$
and $\mu \in [0, \varepsilon]$, system~\eqref{eq:normform_2} generates a 
forward semi-flow on the set 
$\mathfrak{A} := \{(z, v) \in \mathbb{C}^n \times \mathbb{R}^{m} 
\colon \| z \| \le\rho_0, v \in \mathcal{V}\} $. Furthermore, for any
solution $ ( z(t), v(t) ) $ of~\eqref{eq:normform_2} with the initial 
values $ ( z(0), v(0) ) \in \mathfrak{A}$, there exist such a stationary 
point $v_{\ast}$ of $V(\cdot)$ and an instant of time $t_{\ast}$, that
\begin{equation*}
 \| z(t) \| < \rho \sqrt{\varepsilon} < \rho_0,
 \quad
 \| v(t) - v_{\ast} \| < R \varepsilon
 \quad
 \forall t \geq t_{\ast}.
\end{equation*}
\end{proposition}

\begin{proof}
There is a constant $c_1 > 0$, such that if $ \| z \| \leq \rho_0$, then a
quadratic form $\| z \|^2 = \sum_{k = 1}^n | z_{k} |^2$ and the function 
$V(\cdot)$ admit the following estimates for their directional derivatives 
along the vector field of system~\eqref{eq:normform_2}
\begin{gather}
\begin{aligned}
 \| z \|^2 \big| _{\eqref{eq:normform_2}}'
&\leq 2 \varepsilon \sum_{k = 1}^n  \lambda_{k}(v) | z_{k} | ^2
  - 2 \sum_{k, l = 1}^n  a_{kl}(v) | z_{k} | ^2 | z_{l} | ^2
 \\
&\quad +  c_1 (
 \mu \| z \|^2 + \| z \|^{5} + \| z \| \mu^{P + 1}
  ) \\
&\leq  \big[
  (
 2 \varepsilon\lambda^{+} + c_1\mu
  ) \| z \| - (
 2 a_{+} - c_1 \rho_0
  ) \| z \|^{3} + c_1 \mu^{P + 1}
 \big] \| z \|,
\end{aligned}  \nonumber \\
  V(v) \big| _{\eqref{eq:normform_2}}'
  \leq \mu \big[
 \dot{V}_{g}(v) + c_1 \| \nabla V(v) \| (
  \| z \|^2 + \mu
 )
  \big]. \label{eq:est_derV}
\end{gather}
It is easily seen that one can choose the positive numbers 
$\rho_0 < 2a_{+} / c_1$, and $\rho_{\ast} > 0$ in such a way, that for any
$\rho > \rho_{\ast}$, $\mu \in [0, \varepsilon]$, $v \in \operatorname{cl} (\mathcal{V})$,
$\varepsilon \in (0, \varepsilon_1]$, where 
$\varepsilon_1 := \min\{1, \rho_0^2 / \rho^2, \mu_0\}$, the inequality
\begin{equation*}
  \| z \|^2 \big| _{\eqref{eq:normform_2}}' < 0,
\end{equation*}
holds as long as $z$ satisfies $\rho \sqrt{\varepsilon} \leq \| z \| \leq \rho_0$.
It means that $\mathfrak{A}$ is forward semi-invariant and, moreover, there 
exists $t_0 > 0$, such that $\| z(t) \| < \rho\sqrt{\varepsilon}$ for all
 $t > t_0$. Now we can regard $v(t)$ as a solution of a system
$\dot{v} = \mu [ g(v) + \tilde{g}(t, v, \varepsilon, \mu) ] $ defined on 
$[t_0, \infty) \times \operatorname{cl} (\mathcal{V})$ and obtained from the
 corresponding sub-system of~\eqref{eq:normform_2} via the substitution $z = z(t)$.
 Obviously, there exists such a constant $c_{2} > 0$, that
\begin{equation*}
 \| \tilde{g}(t, v, \varepsilon, \mu) \|
  \leq c_{2} \varepsilon
 \quad
 \forall (t, v, \varepsilon, \mu)
  \in [t_0, \infty)
 \times \operatorname{cl} (\mathcal{V})
 \times(0, \varepsilon_1]
 \times[0, \varepsilon].
\end{equation*}

On account of (H1), \eqref{eq:est_dotV_g} and \eqref{eq:est_derV}, after 
an appropriate additional correction of $\delta$ in \eqref{eq:est_dotV_g}, 
the final part of the proof follows from Lemma~\ref{lemma:W'f < 0 - 2}
in the case when $f = \mu g$, $f_1 = \mu\varepsilon\tilde{g}$, $W = V$, 
$\epsilon = \varepsilon c_{2}$, $\theta\propto\mu$. In particular, 
if we find $M_{\ast}$ and $\epsilon_0$ from the lemma, we can set
$R_{\ast} = M_{\ast}$, $R = M$ and 
$\varepsilon_0 = \min\{ \epsilon_0 / c_{2}, \varepsilon_1 \}$.
\end{proof}

\begin{corollary}
Let $0<\mu< \varepsilon\kappa / (2c_1)$ and let $(t_0,t_1)\in\mathbb{R}_{+}$
be such an interval, that
\begin{gather*}
\mu^{P}<\| z(t)\| <\rho_0,\quad v(t)\in\mathcal{V}_{-}^{\kappa}\quad
\forall t\in(t_0,t_1)
\end{gather*}
Then $\| z(t)\| \le\mathrm{e}^{-\varepsilon\kappa(t-t_0)/2}\| z(t_0)\| $
for all $t\in(t_0,t_1)$.
\end{corollary}

In fact, for $\mu^{P}\le\| z\| \le\rho_0$, $v\in\mathcal{V}_{-}^{\kappa}$
and $0<\mu<\varepsilon\kappa/(2c_1)$, in the same way as in the
proof of Proposition~\ref{proposition:sqrt - mu - neighbor}, we obtain the
inequality
\begin{align*}
\| z\| ^2\big|_{\eqref{eq:normform_2}}'
&\leq\big[(-2\varepsilon\kappa+c_1\mu)\| z\| 
 -(2a_{+}-c_1\rho_0)\| z\| ^{3}+c_1\mu^{P+1}\big]\| z\| \\
&\le\big[-\frac{3\varepsilon\kappa}{2}\| z\| 
+\frac{\varepsilon\kappa}{2}\mu^{P}\big]\| z\| 
\le-\frac{\varepsilon\kappa}{2}\| z\| ^2.
\end{align*}
The above corollary proves the statement (2) of the main theorem.

Hypothesis (H3) and Proposition~\ref{proposition:sqrt - mu - neighbor}
allow us to focus on system~\eqref{eq:normform_2} defined on the set
\begin{equation*}
 \{
  (r, v)
  \colon
  \| z \| < \rho\sqrt{\varepsilon},
  \| v - v_{\ast} \| < R \varepsilon
 \}.
\end{equation*}
Hereafter, we will require the numbers $\rho$ and $R$ to be large enough.

Without loss of generality we may suppose that $v_{\ast} = 0$ in 
Proposition~\ref{proposition:sqrt - mu - neighbor}.
Then, having applied the scaling $z \mapsto \sqrt{\varepsilon}z$, 
$v \mapsto \varepsilon v$ to~\eqref{eq:normform_2}, we obtain the system
\begin{equation}\label{eq:scalingNF}
\begin{gathered}
\begin{aligned}
 \dot{z}_{k}
& = \Big[
 \mathrm{i} \omega_{k} (\varepsilon v)
 + \varepsilon\lambda_{k}(\varepsilon v)
 - \varepsilon \sum_{l = 1}^n (
  a_{kl}(\varepsilon v)
  + \mathrm{i}b_{kl}(\varepsilon v)
 ) | z_{l} | ^2
  \Big] z_{k}
 \\
&\quad +\sum_{5 \leq 2 | \mathbf{p} | + 1 \leq N}
  \varepsilon^{| \mathbf{p} |}
  h_{0, k, \mathbf{p}}(\varepsilon v)
  \overrightarrow{| z |}^{2\mathbf{p}} z_{k}
 + \sum_{j = 1}^{P}
  \mu^{j} \eta_{j, k}(\varepsilon v) z_{k}
 \\
&\quad +\sum_{j = 1}^{P}
 \mu^{j} \Big[
  \sum_{3 \leq 2 | \mathbf{p} | + 1 \leq N}
   \varepsilon^{| \mathbf{p} |}
   h_{j, k, \mathbf{p}}(\varepsilon v)
   \overrightarrow{| z |}^{2\mathbf{p}}
  \Big] z_{k}
  + O \big(
 \varepsilon^{N / 2} \| z \| ^{N + 1}
 + \mu^{P}
  \big),
\end{aligned} \\
 k \in \{1, \ldots, n\},
 \\
 \begin{aligned}
\dot{v}
&= \frac{\mu}{\varepsilon} \Big[
 g(\varepsilon v)
 + \sum_{2 \leq 2 | \mathbf{p} | \leq N}
  \varepsilon^{ | \mathbf{p} | }
  g_{0, \mathbf{p}}(\varepsilon v)
  \overrightarrow{| z |}^{2\mathbf{p}}
  \Big]
 \\
&\quad +\frac{\mu}{\varepsilon} \sum_{j = 1}^{P}
  \mu^{j} \sum_{0 \leq 2 | \mathbf{p} | \leq N}
 \varepsilon^{| \mathbf{p} |}
 g_{j, \mathbf{p}}(\varepsilon v)
  \overrightarrow{| z | }^{2\mathbf{p}}
 + O \big(
  \mu \varepsilon^{(N - 1) / 2} \| z \| ^{N + 1}
  + \mu^{P + 2} / \varepsilon
 \big)
\end{aligned}
\end{gathered}
\end{equation}
defined on the set
\begin{equation*}
 \{
  ( z, v, \varepsilon, \mu )
  \colon
  \| z \| \leq \rho,
  \| v \| \leq R,
  \varepsilon \in (0, \varepsilon_0],
  \mu \in [0, \varepsilon]
 \}.
\end{equation*}
In this system we constrain the parameter $\mu$ to be 
$\mu = \varepsilon\varsigma$ with $\varsigma$ being an arbitrary 
fixed number satisfying
\begin{equation*}
 0 < \varsigma \leq \varsigma_0
  := \frac{1}{2} \min_{1 \leq j \leq n}
 \frac{
  \lambda_{j}(0)
 }{
  | \eta_{1, j}(0) |
 }.
\end{equation*}
Note that such a condition ensures the validity of the inequality
\begin{equation*}
 \alpha_{k} := \lambda_{k}(0) + \varsigma \operatorname{Re} \eta_{1, k}(0)
  \geq \lambda_{k}(0) / 2 > 0
 \quad
 \forall k \in \{1, \ldots, n \}.
\end{equation*}

Now, if we recall the earlier imposed condition $P\ge N/2$ and 
if for a vector $x = (x_1, \ldots, x_{d})$ we define 
$D[x] := \operatorname{diag} ( x_1, \ldots, x_{d} )$, then 
system~\eqref{eq:scalingNF} can be presented in the form
\begin{gather}
 \label{eq:vect - norm - form - z}
\begin{aligned}
 \dot{z} &= D [
  \mathrm{i} (
 \omega_0 + \varepsilon (
  \beta + \Omega' v - B \overrightarrow{| z |}^2
 )
  )] z
 \\
&\quad +[ \varepsilon (
 \alpha - A \overrightarrow{| z |}^2
  ) + \varepsilon^2 \hat{h} (
 \overrightarrow{| z |}^2, v
  )
 ] z + O ( \varepsilon^{N / 2} ),
\end{aligned} \\
 \label{eq:vect - norm - form - v}
 \dot{v} = \varepsilon\varsigma [
  \Gamma v
  + \Upsilon \overrightarrow{| z |}^2
  + \varepsilon \hat{g}(\overrightarrow{| z |}^2, v)
  + O(\varepsilon^{(N + 1) / 2})
 ]
\end{gather}
where, for the sake of notations simplicity, we assign
\begin{gather*}
 \omega_0 := (
  \omega_{01}(0), \ldots\omega_{0n}(0)
 ),
 \quad
 \eta_1 := (
  \eta_{1, 1}(0), \ldots, \eta_{1, n}(0)
 ),
 \quad
 \alpha := \lambda(0) + \varsigma \operatorname{Re} \eta_1(0),
 \\
 \beta := \zeta \operatorname{Im} \eta_1(0),
 \quad
 A := \big\{ a_{kl}(0) \big\}_{k, l = 1}^n,
 \quad
 B := \big\{ b_{kl}(0) \big\}_{k, l = 1}^n,
 \\
 \Omega' := \Big\{
  \frac{
 \partial \omega_{0k}(0)
  }{
 \partial v_{j}
  }
\Big\}_{k = 1, }^n
 \quad
 \Gamma := g'(0),
 \quad
 \Upsilon x := \sum_{| \mathbf{p} | = 1}
  g_{0, \mathbf{p}}(0)x^{\mathbf{p}}.
\end{gather*}
The definitions of the remainder terms $\hat{h} ( \cdot ) $ and 
$\hat{g}(\cdot)$ inside the square brackets in the right-hand sides 
of \eqref{eq:vect - norm - form - z}, \eqref{eq:vect - norm - form - v} 
is obvious. Recall that we agreed not to mention $\varepsilon$ directly 
as functions arguments.

\begin{proposition}\label{proposition:hyp - equil}
For all $\varepsilon \in (0, \varepsilon_0]$, with $\varepsilon_0 > 0$
being sufficiently small, 
system~\eqref{eq:vect - norm - form - z}--\eqref{eq:vect - norm - form - v} 
has an equilibrium $(z^0, v^0) = ( z^0(\varepsilon), v^0(\varepsilon) )$,
such that
\begin{equation}\label{eq:equilibz0v0}
 z^0(\varepsilon) = O ( \varepsilon^{N / 2} ),
 \quad
 v^0(\varepsilon) = O(\varepsilon).
\end{equation}
If $( z(t), v(t) )$ is a solution 
of~\eqref{eq:vect - norm - form - z}--\eqref{eq:vect - norm - form - v}, 
such that $\| z(t) \| < \sqrt{\varepsilon}$ and $\| v(t) \| < R$ for all $t > 0$, 
then
\begin{equation}\label{eq:z - to - z0}
 \lim_{t \to \infty} \left[
  \| z(t) - z^0 \| + \| v(t) - v^0 \|
 \right] = 0,
\end{equation}
and the set of all such solutions forms a manifold whose real dimension 
equals the number of eigenvalues of $\Gamma$ with negative real parts. 
In the case when all eigenvalues of $\Gamma$ have positive real parts the 
only solution with the stated property is the equilibrium $( z^0, v^0 )$.
\end{proposition}

\begin{proof}
The existence of the equilibrium $( z^0(\varepsilon), v^0(\varepsilon) )$ 
satisfying~\eqref{eq:equilibz0v0} directly follows from the implicit 
function theorem. If $( z(t), v(t) )$ is a solution
of \eqref{eq:vect - norm - form - z}--\eqref{eq:vect - norm - form - v}, 
such that $\| z(t) \| < \sqrt{\varepsilon}$ and $\| v(t) \| < R$ for all $t > 0$, 
then the functions
\begin{equation*}
 w(t) := \exp \Big[
  -\mathrm{i} \int_0^{t} D \Big[
 \Big(
  \omega_0 + \varepsilon (
   \Omega' v(s)
   - B \overrightarrow{| z(s) |}^2
  )
 \Big)
  \Big] \mathrm{d} s
 \Big] z(t)
\end{equation*}
and $v(t)$ satisfy a system of the form
\begin{gather}
 \label{eq:eq - w}
 \dot{w} = D \left[
  \varepsilon (
 \alpha - A \overrightarrow{| w |}^2
  ) + \varepsilon^2 \hat{h} (
 \overrightarrow{ | w |}^2, v
  )
 \right] w + O ( \varepsilon^{N / 2} ),
 \\
 \label{eq:eq - v}
 \dot{v} = \varepsilon\varsigma \left[
  \Gamma v
  + \Upsilon \overrightarrow{| w |}^2
  + \varepsilon\hat{g}(\overrightarrow{| w | }^2, v)
  + O(\varepsilon^{(N + 1) / 2})
 \right],
\end{gather}
which also has a solution
\begin{gather*}
 w = w^0(t) := \exp \Big[
  -\mathrm{i} \int_0^{t} D \big[
 (
  \omega_0 + \varepsilon (
   \Omega' v(s)
   - B \overrightarrow{| z(s) |}^2
  )
 )
  \big] \mathrm{d} s
 \Big] z^0,
 \\
 v = v^0,
 \quad
 t \geq 0.
\end{gather*}
Hence, the realification of system~\eqref{eq:eq - w}--\eqref{eq:eq - v} 
has a pair of solutions
\begin{equation*}
 \xi(t) := (
  \operatorname{Re} w(t), \operatorname{Im} w(t), v(t)
 ),
 \quad
 \xi^0(t) := (
  \operatorname{Re} w^0(t), \operatorname{Im} w^0(t), v^0
 ).
\end{equation*}

One can consider the difference $\eta(t) := \xi(t) - \xi^0(t)$ to be a bounded 
solution for $t \geq 0$ of the linear system 
$\dot{\eta} = \varepsilon [ \mathcal{A} + \mathcal{A}_1(t;\varepsilon)] \eta$ where
\begin{equation*}
 \mathcal{A} = \operatorname{diag}[
  D[\alpha], D[\alpha], \varsigma\Gamma],
\end{equation*}
and $\sup_{t \geq 0} \| \mathcal{A}_1(t; \varepsilon) \| = O(\sqrt{\varepsilon})$. 
Such a system is hyperbolic on $[0, \infty)$ if $\varepsilon$ is small enough, 
and each of its bounded solutions approaches zero as $t \to \infty$. 
This yields \eqref{eq:z - to - z0}.

Next, it is not hard to see that at $( z^0, v^0 )$ the Jacobi matrix of the 
right-hand side of 
system~\eqref{eq:vect - norm - form - z}--\eqref{eq:vect - norm - form - v}
 realification has the form
\begin{equation*}
 \pm \mathrm{i} \omega_{0k}
 + \varepsilon (
  \alpha_{k} \pm \mathrm{i} \beta_{k}
 )
 + o ( \varepsilon ),
 \quad
 k \in \{ 1, \ldots, n \},
 \quad
 \varepsilon\varsigma\gamma_{j} + o(\varepsilon),
 \quad
 j \in \{ 1, \ldots, m \}
\end{equation*}
where $\gamma_1, \ldots, \gamma_{m}$ are the eigenvalues of $\Gamma$ 
counted according to multiplicities. All of these numbers have non-zero 
real parts. It is well known that in this case the set of all solutions 
of the realification 
of~\eqref{eq:vect - norm - form - z}--\eqref{eq:vect - norm - form - v} 
which approach the equilibrium $\xi^0(\varepsilon)$ forms the so-called 
stable manifold, whose dimension equals the number of eigenvalues of 
$\mathcal{A}(\varepsilon)$ with negative real parts (see, e.~g. \cite{Har70}). 
In the case when all of the eigenvalues of $\Gamma$ have positive real parts, 
the equilibrium $( \operatorname{Re} z^0, \operatorname{Im} z^0, v^0 )$ is completely unstable.
\end{proof}

\begin{proposition}\label{proposition:important}
There exists a constant $c_{3}>0$, such that for sufficiently small 
$\varepsilon_0 > 0$ and for any $\varepsilon \in (0, \varepsilon_0]$
the following assertion is valid. If $(z(t), v(t))$ is a solution 
of~\eqref{eq:vect - norm - form - z}--\eqref{eq:vect - norm - form - v}, 
such that $\| z(0) \| \leq \rho$, $| z_{k}(0) | \geq \sqrt{c_{3}\varepsilon^{N - 2}}$ 
for some $k \in \{ 1, \ldots, n \}$ and $\| v(t) \| < R$ for all $t>0$, then
\begin{equation*}
 | z_{k}(t) |
  > \sqrt{c_{3}\varepsilon^{N - 2}},
 \quad
 \| z(t) \| < \rho
 \quad
 \forall t > 0,
\end{equation*}
and there exists such $t_{\ast}>0$, that
\begin{equation*}
 | z_{k}(t) |
  > \sqrt{\frac{\alpha_{k}}{2a_{kk}}},
 \quad
 \| z(t) \| < \sqrt{\frac{
  2\max_{1 \leq k \leq n}
 \alpha_{k}
  }{
 a_{+}^0
  }}
 \quad
 \forall t>t_{\ast},
\end{equation*}
where
\begin{equation*}
 a_{+}^0 := \min \{
  \langle Ar, r \rangle
  \colon
  r \in \mathbb{R}_{+}^n, | r | = 1
 \}.
\end{equation*}
\end{proposition}

\begin{proof}
It is sufficient to note that there is such $c_{4}>0$, that for sufficiently 
small $\varepsilon_0 > 0$ and for any $\varepsilon \in (0, \varepsilon_0]$
the following inequalities are satisfied.
\begin{align*}
 \| z \|^2 \big|_{\eqref{eq:vect - norm - form - z}}'
&\leq 2 \varepsilon \langle
 \overrightarrow{| z |}^2,
 \alpha - A \overrightarrow{| z |}^2
  \rangle + c_{4} (
 \varepsilon^2 \| z \|^2
 + \varepsilon^{N / 2} \| z \|
  )
 \\
&\leq2 \varepsilon \big[
  \max_{1 \leq k \leq n} \alpha_{k}
  - a_{+}^0 \| z \|^2
  + c_{4}\varepsilon
 \big] \| z \|^2
 + c_{4}\varepsilon^{N / 2} \| z \| < 0
\end{align*}
if $2 \max_{1 \leq k \leq n}\alpha_{k} / a_{+}^0 \leq \| z \|^2 \leq \rho^2$,
and, in view of (H4),
\begin{equation*}
  | z_{k} |^2 \big|_{\eqref{eq:vect - norm - form - z}}'
  \geq 2 \varepsilon | z_{k} |^2 \left[
 \alpha_{k} - a_{kk} | z_{k} |^2
 - c_{4}\varepsilon
 - c_{4}\varepsilon^{(N - 2) / 2} / | z_{k} |
  \right] > 0
\end{equation*}
if $c_{3}\varepsilon^{N - 2} \leq | z_{k} |^2 \leq \alpha_{k} / (2 a_{kk})$,
 where $c_{3}$ is a sufficiently large positive constant.
\end{proof}

\begin{definition}\label{def:sol - type} \rm
Let $\{ i_1, \ldots, i_{s} \}$ be an ordered collection of distinct natural 
numbers not exceeding $n$. We will say that a solution $( z(t), v(t) )$ 
of~\eqref{eq:vect - norm - form - z}--\eqref{eq:vect - norm - form - v} 
is of type $( i_1, \ldots, i_{s} )$ if
\begin{equation*}
 \| v(t) \| < R,
 \quad
 | z_{k}(t) | < \varepsilon
 \quad
 \forall k \not \in \{ i_1, \ldots, i_{s}\},
 \quad
 \forall t > 0
\end{equation*}
and there exists such a moment of time $t^{\ast} > 0$, that
\begin{equation*}
 | z_{k}(t) |
  > \sqrt{\frac{\alpha_{k}}{2a_{kk}}}
 \quad
 \forall k \in \{ i_1, \ldots, i_{s} \},
 \quad
 \| z(t) \| < \sqrt{\frac{
  2 \max_{1 \leq k \leq n}
 \alpha_{k}
 }{
 a_{+}^0
 }}
 \quad
 \forall t > t^{\ast}.
\end{equation*}
\end{definition}

As a consequence of Propositions~\ref{proposition:hyp - equil}
and~\ref{proposition:important}, we obtain the following result.

\begin{proposition}
Let $(z(t), v(t))$ be a solution of
\eqref{eq:vect - norm - form - z}--\eqref{eq:vect - norm - form - v}, 
such that $\| z(0) \| \leq \rho$ and $\| v(t) \| \leq R$ for all 
$t \geq 0$. If $\varepsilon_0$ is sufficiently small and
$\varepsilon \in (0, \varepsilon_0]$, then either this solution tends to
 the equilibrium $( z^0, v^0 )$ as $t \to \infty$, or it is of type 
$( i_1, \ldots, i_{s} )$ for some $s \in \{ 1, \ldots, n \}$.
\end{proposition}

\begin{proof}
If $\| z(t) \| < \sqrt{\varepsilon}$ for all $t > 0$,
then $( z(t), v(t) )$ is a solution of
\eqref{eq:vect - norm - form - z}--\eqref{eq:vect - norm - form - v}
 from Proposition~\ref{proposition:hyp - equil}, and therefore tends to
 $( z^0, v^0 )$ as $t\to\infty$. If $(z(t), v(t))$ does not possess the 
aforementioned property, then there exist $k \in \{ 1, \ldots, n\} $ and 
$t_1 \geq 0$, such that 
$| z_{k}(t_1) | ^2 > \varepsilon / n \geq c_{3}\varepsilon^{N - 2}$. Hence, 
the solution $( z(t - t_1), v(t - t_1) )$ satisfies the conditions of 
Proposition~\ref{proposition:important}. Now it becomes apparent that in such case,
basing on Proposition~\ref{proposition:important}, one can decompose a set
$\{ 1, \ldots, n\}$ into two ordered subsets, $\{ i_1, \ldots, i_{s}\} \ni k$ 
and $\{ j_1, \ldots, j_{n - s}\} \subset \{ 1, \ldots, n\} 
\setminus \{ i_1, \ldots, i_{s}\}$, and choose $t^{\ast} > 0$ in such a way, 
that $| z_i(t) | ^2 > \alpha_i / (2 a_{ii})$ for all $t > t^{\ast}$ if $i$
belongs to the first subset, whereas 
$| z_{j}(t) | ^2 \leq c_{3}\varepsilon^{N - 2} < \varepsilon^2$ for all 
$t > 0$ if $j$ belongs to the second subset, with the latter being empty when
 $s = n$. Besides that, $\| z(t) \|$ meets the imposed requirements for all 
$t > t^{\ast}$.
\end{proof}


\section{Ultimate behavior of solutions of type $(1, \ldots, n)$}
\label{sec:ultimate1..n}

To study the final behavior of solutions of type $(1, \ldots, n )$, 
introduce the polar-like coordinates 
$z_{k} = \sqrt{r_{k}} \mathrm{e}^{\mathrm{i} \varphi_{k}}$, $k = 1, \ldots, n$, 
and set $r = ( r_1, \ldots, r_{n} )$. This transforms 
system~\eqref{eq:vect - norm - form - z}--\eqref{eq:vect - norm - form - v} into
\begin{gather}
 \label{eq:polar - r}
 \dot{r} = 2\varepsilon D[r] \left[
  \alpha - Ar + \varepsilon\hat{a}(r, v)
 \right]
 + \varepsilon^{N / 2} D^{1 / 2}[r] \tilde{a}(r, v, \varphi),
 \\
 \label{eq:polar - v}
 \dot{v} = \varepsilon\varsigma \left[
  \Upsilon r + \Gamma v + \varepsilon\hat{g}(r, v)
  + \varepsilon^{(N + 1) / 2} \tilde{g}(r, v, \varphi)
 \right],
 \\
 \label{eq:polar - phi}
 \dot{\varphi} = \omega_0 + \varepsilon (
  \beta + \Omega' v - B r
 )
 + \varepsilon^2 \hat{b}(r, v)
 + \varepsilon^{N / 2} D^{-1 / 2}[r] \tilde{b}(r, v, \varphi)
\end{gather}
where $\hat{a}(r, v) = \operatorname{Re} \hat{h} ( r, v )$, $\hat{b}(r, v) = \operatorname{Im} \hat{h} ( r, v )$,
and $\tilde{a}(r, v, \varphi)$, $\tilde{b}(r, v, \varphi)$, 
$\tilde{g}(r, v, \varphi)$
are defined by the remainder terms of
\eqref{eq:vect - norm - form - z}--\eqref{eq:vect - norm - form - v}. 
On ground of Lemma~\ref{lemma:A inv - pos} and
Proposition~\ref{proposition:sqrt - mu - neighbor}, we can consider
$\rho$ and $R$ to be so large, that
\begin{equation}\label{eq:A^ - 1alpha < }
 | A^{-1}\alpha |
  \leq \frac{
 \max_{1 \leq k \leq n}\alpha_{k}
  }{
 a_{+}^0
  } < \frac{\rho}{2},
 \quad
 \| \Gamma^{-1}\Upsilon A^{-1}\alpha \| < \frac{R}{2}.
\end{equation}

\begin{proposition}\label{proposition:r - * v - *}
For sufficiently small $\varepsilon_0 > 0$ and for any
$\varepsilon \in (0, \varepsilon_0]$, the system
\begin{equation*}
 \alpha - Ar + \varepsilon\hat{a}(r, v) = 0,
 \quad
 \Upsilon r + \Gamma v + \varepsilon\hat{g}(r, v) = 0
\end{equation*}
has a solution
\begin{equation*}
 r = r_{\ast}(\varepsilon) := A^{-1}\alpha + O(\varepsilon),
 \quad
 v = v_{\ast}(\varepsilon) := - \Gamma^{-1}\Upsilon A^{-1}\alpha + O(\varepsilon),
\end{equation*}
such that
\begin{equation*}
 r_{\ast k}(\varepsilon) > \frac{2\alpha_{k}}{3a_{kk}},
 \quad
 k \in \{ 1, \ldots, n\},
 \quad
 | r_{\ast}(\varepsilon) |
  < \frac{
 3\max_{1 \leq k \leq n}\alpha_{k}
  }{
 2a_{+}^0
  },
 \quad
 \| v_{\ast}(\varepsilon) \| < \frac{2R}{3}.
\end{equation*}
\end{proposition}

\begin{proof}
Taking into account hypothesis~(H4), Lemma~\ref{lemma:A inv - pos} and
inequalities~\eqref{eq:A^ - 1alpha < }, the desired result follows from the 
implicit function theorem.
\end{proof}

\begin{proposition}\label{proposition:|r - r_8 * |}
There exist such positive numbers $c_{5}$, $\varsigma_0$ and $\varepsilon_0$,
that for any $\varsigma \in ( 0, \varsigma_0 )$,
$\varepsilon \in (0, \varepsilon_0)$ the following assertion is true.
If $( z(t), v(t) )$ is a solution of type $( 1, \ldots, n )$ of
 system~\eqref{eq:vect - norm - form - z}--\eqref{eq:vect - norm - form - v} 
and $r(t) := {\overrightarrow{ | z(t) | }}^2$ then there is such $t^{\ast} > 0$, 
that
\begin{equation*}
 \sqrt{\| r(t) - r_{\ast}(\varepsilon) \|^2
  + \| v(t) - v_{\ast}(\varepsilon) \| ^2}
 < \frac{c_{5} \varepsilon^{(N - 2) / 2}}{\varsigma}
 \quad
 \forall t > t^{\ast}.
\end{equation*}
\end{proposition}

\begin{proof}
Let $( r(t), v(t), \varphi(t) )$ represent the solution $( z(t), v(t) )$ of type 
$(1, \ldots, n)$ in the polar-like coordinates. On account of 
\eqref{eq:polar - r}--\eqref{eq:polar - phi} and 
Proposition~\ref{proposition:r - * v - *}, the pair $( r(t), v(t) )$
satisfies the system
\begin{gather*}
\begin{aligned}
 \dot{r} &= 2 \varepsilon D[r] \big[
  (
 - A + \varepsilon\hat{A}_{r}(r, v)
  ) (r - r_{\ast})
  + \varepsilon\hat{A}_{v}(r, v)(v - v_{\ast})
 \big]
 \\
&\quad +\varepsilon^{N / 2} D^{1 / 2}[r] \tilde{a}(r, v, \varphi(t)),
\end{aligned} \\
\begin{aligned}
 \dot{v} &= \varepsilon\varsigma \big[
  (
 \Upsilon + \varepsilon\hat{G}_{r}(r, v)
  ) (r - r_{\ast})
  + (
 \Gamma + \varepsilon\hat{G}_{v}(r, v)
  ) (v - v_{\ast})
 \\
&\quad
  + \varepsilon^{(N + 1) / 2} \tilde{g}(r, v, \varphi(t))
 \big],
\end{aligned}
\end{gather*}
where
\begin{gather*}
 \hat{A}_{r}(r, v)
   := \int_0^{1}  \frac{
 \partial\hat{a}(r, v)
  }{
 \partial r
  } \Big|_{\begin{smallmatrix}
 r \mapsto sr + (1 - s) r_{\ast}
 \\
 v \mapsto sv + (1 - s)v_{\ast}
  \end{smallmatrix}} \mathrm{d}s,
 \\
 \hat{A}_{v}(r, v)
   := \int_0^{1}  \frac{
 \partial\hat{a}(r, v)
  }{
 \partial v
  } \Big|_{\begin{smallmatrix}
 r \mapsto sr + (1 - s) r_{\ast}
 \\
 v \mapsto sv + (1 - s)v_{\ast}
  \end{smallmatrix}} \mathrm{d}s,
 \\
 \hat{G}_{r}(r, v)
   := \int_0^{1}  \frac{
 \partial\hat{g}(r, v)
  }{
 \partial r
  } \Big|_{\begin{smallmatrix}
 r \mapsto sr + (1 - s) r_{\ast}
 \\
 v \mapsto sv + (1 - s)v_{\ast}
  \end{smallmatrix}} \mathrm{d}s,
 \\
 \hat{G}_{v}(r, v)
   := \int_0^{1}  \frac{
 \partial\hat{g}(r, v)
  }{
 \partial v
  } \Big|_{\begin{smallmatrix}
 r \mapsto sr + (1 - s) r_{\ast}
 \\
 v \mapsto sv + (1 - s)v_{\ast}
  \end{smallmatrix}} \mathrm{d}s.
\end{gather*}

By Definition \ref{def:sol - type}, there exists $t^{\ast} > 0$, such that for 
all $t > t^{\ast}$ a point $(r(t), v(t))$ belongs to the domain
\begin{align*}
 \mathcal{D} := \Big\{&
  ( r, v ) \in \mathbb{R}_{+}^n
 \times \mathbb{R}^{m}
  \colon
  r_{k} > \frac{\alpha_{k}}{2a_{kk}}
  \;  \forall k \in \{ 1, \ldots, n\}, \\
 & | r | < \frac{
 2\max_{1 \leq k \leq n} \alpha_{k}
  }{
 a_{+}^0
  },
  \quad
  \| v \| < R
 \Big\},
\end{align*}
which contains the unique equilibrium $(r_{\ast}, v_{\ast})$ of the system
\begin{equation}\label{eq:trunc - sys}
\begin{gathered}
 \dot{r} = 2\varepsilon D[r] \left[
  \big(
 - A + \varepsilon\hat{A}_{r}(r, v)
  \big) (r - r_{\ast})
  + \varepsilon \hat{A}_{v}(r, v) (v - v_{\ast})
 \right],
 \\
 \dot{v} = \varepsilon\varsigma \left[
  (
 \Upsilon + \varepsilon\hat{G}_{r}(r, v)
  ) (r - r_{\ast}) + \big(
 \Gamma + \varepsilon\hat{G}_{v}(r, v)
  \big) (v - v_{\ast})
 \right] ,
\end{gathered}
\end{equation}
This equilibrium is the unique stationary point of the Morse function
\begin{equation*}
 W(r, v) := \sum_{i = 1}^n
  \big(
 r_i + r_{\ast i} \ln (
  \frac{r_{\ast i}}{r_i}
 ) - r_{\ast i}
  \big)
  + \frac{1}{2} \langle
 V''(0)(v - v_{\ast}), v - v_{\ast}
  \rangle
\end{equation*}
in $\operatorname{cl} (\mathcal{D}$).

The first inequality~\eqref{eq:est_dotV_g} yields 
$\langle V''(0)\Gamma v, v \rangle \leq - \nu_{\ast} \| v \| ^2$
for any $v \in \mathbb{R}^{m}$. Therefore, there is such a constant $c_{6} > 0$, 
that for sufficiently small $\varepsilon_0$ and for all
 $\varepsilon \in (0, \varepsilon_0)$, $(r, v) \in \operatorname{cl} (\mathcal{D})$
the following inequality holds.
\begin{align*}
& W(r, v) |_{\eqref{eq:trunc - sys}}' \\
 & \leq 2\varepsilon \langle
 r - r_{\ast}, \big(
  - A + \varepsilon\hat{A}_{r}(r, v)
 \big) (r - r_{\ast})
 + \varepsilon\hat{A}_{v}(r, v)(v - v_{\ast})
  \rangle
 \\
&\quad  +\varepsilon\varsigma \langle
  V''(0)(v - v_{\ast}), \big(
 \Upsilon + \varepsilon\hat{G}_{r}(r, v)
  \big) (r - r_{\ast}) + (
 \Gamma + \varepsilon\hat{G}_{v}(r, v)
  ) (v - v_{\ast})
 \rangle
 \\
&\leq- \varepsilon a^0 \| r - r_{\ast} \|^2
 + 2 c_{6} (
 \varepsilon^2 + \varepsilon\varsigma
  )
  \| r - r_{\ast} \|
  \| v - v_{\ast} \|
  - \frac{\varepsilon\varsigma\nu_{\ast}}{2}
  \| v - v_{\ast} \|^2,
\end{align*}
where
\begin{equation*}
 a^0 := \min \{
  \langle A \zeta, \zeta \rangle
  \colon
  \zeta \in \mathbb{R}^n, \| r \| = 1
 \}.
\end{equation*}

Now, observe that for sufficiently small positive $\varsigma_0$ and
$\varepsilon_0$, and for any $\varsigma \in (0, \varsigma_0)$ and
$\varepsilon \in (0, \varepsilon_0)$ the smallest eigenvalue of the matrix
\begin{equation*}
 \begin{pmatrix}
  a^0 & - c_{6}(\varepsilon + \varsigma)
  \\
  - c_{6}(\varepsilon + \varsigma) & \varsigma\nu_{\ast} / 2
 \end{pmatrix}
\end{equation*}
exceeds $\varsigma\kappa$, where
\begin{equation*}
 \kappa_0 := \frac{1}{2} \cdot \frac{
  a^0\nu_{\ast}
 }{
  2a^0 + \varsigma\nu_{\ast}
 }.
\end{equation*}
Hence,
\begin{equation*}
W(r, v) \big|_{\eqref{eq:trunc - sys}}'
  \leq - \varepsilon\varsigma\kappa_0 \big[
 \| r - r_{\ast} \|^2
 + \| v - v_{\ast} \|^2
  \big],
\end{equation*}
as soon as $(r, v) \in \operatorname{cl} (\mathcal{D})$, $\varsigma \in (0, \varsigma_0)$
and $\varepsilon \in (0, \varepsilon_0)$. It only remains to apply
Lemma~\ref{lemma:W'f < 0 - 2} in the case of unique stationary point of
Morse function with $\theta\propto\varepsilon\varsigma$ and 
$\epsilon\propto\varepsilon^{N / 2}$ to finish this proof.
\end{proof}

Next, we are going to utilize results on the existence of invariant tori 
obtained in~\cite{Hale61}. To do so, we introduce a new vector variable 
$\xi \in \mathbb{R}^{n + m}$ via the formula
\begin{equation*}
 (r, v) = ( r_{\ast}, v_{\ast} ) + \varepsilon^{\sigma}\xi
\end{equation*}
with $0 < \sigma < 1 / 2$, and define the block matrices
\begin{gather*}
 \mathcal{B} := \begin{pmatrix}
  - 2D[r_{\ast}]A & 0
  \\
  \varsigma\Upsilon & \varsigma\Gamma
 \end{pmatrix},
 \quad
 \mathcal{B}_1(\varepsilon^{\sigma}\xi) := 
  \begin{pmatrix}
 - 2D[r - r_{\ast}]A & 0
 \\
 0 & 0
  \end{pmatrix}
 \Big|_{(r, v) = (r_{\ast}, v_{\ast}) + \varepsilon^{\sigma}\xi}\,,
 \\
 \hat{\mathcal{B}}(\varepsilon^{\sigma}\xi) := 
  \begin{pmatrix}
 2D[r]\hat{A}_{r}(r, v) & 2D[r]\hat{A}_{v}(r, v)
 \\
 \varsigma\hat{G}_{r}(r, v) & \varsigma\hat{G}_{v}(r, v)
  \end{pmatrix}
 \Big|_{(r, v) = (r_{\ast}, v_{\ast}) + \varepsilon^{\gamma}\xi},
\end{gather*}
the vector functions
\begin{gather*}
 \tilde{f}(\varepsilon^{\sigma}\xi, \varphi) := 
  \big(
 D^{1 / 2}[r]\tilde{a}(r, v, \varphi),
 \varepsilon^{3 / 2}\varsigma\tilde{g}(r, v, \varphi)
  \big)
 \big|_{(r, v) = (r_{\ast}, v_{\ast})
 + \varepsilon^{\sigma}\xi}\,,
 \\
 \hat{\Phi}(\varepsilon^{\sigma}\xi, \varphi)
  :=  \hat{b}(r, v) \big|_{
 (r, v) = (r_{\ast}, v_{\ast}) + \varepsilon^{\sigma}\xi
  },
 \\
 \tilde{\Phi}(\varepsilon^{\sigma}\xi, \varphi)
  := 
 \varepsilon^{(N - 4) / 2} D^{-1 / 2}[r] \tilde{b}(r, v, \varphi)
  \big|_{(r, v) = (r_{\ast}, v_{\ast}) + \varepsilon^{\sigma}\xi},
\end{gather*}
the vector $\omega_{\ast} := \omega_0 + \varepsilon(\beta + \Omega'v_{\ast}
 - B r_{\ast})$ and, finally, the 
$(n \times (n + m))$-matrix $\Xi := [ - B;\Omega']$.
After that, system~\eqref{eq:polar - r}--\eqref{eq:polar - phi} takes the 
following  form in the variables $(\xi, \varphi)$.
\begin{gather}
 \label{eq:dotxi}
 \dot{\xi} = \varepsilon [
  \mathcal{B} + \mathcal{B}_1(\varepsilon^{\sigma}\xi)
  + \varepsilon\hat{\mathcal{B}}(\varepsilon^{\sigma}\xi)
 ] \xi
 + \varepsilon^{N / 2 - \gamma} \tilde{f}(\varepsilon^{\sigma}\xi, \varphi),
 \\
 \label{eq:dotphi}
 \dot{\varphi} = \omega_{\ast} + \varepsilon^{1 + \gamma}\Xi\xi
  + \varepsilon^2\hat{\Phi}(\varepsilon^{\sigma}\xi)
  + \varepsilon^{N / 2}\tilde{\Phi}(\varepsilon^{\sigma}\xi, \varphi).
\end{gather}
Assuming that $\varepsilon$ is small enough for $\omega_{\ast}$ to be positive, 
we hereby arrive at the perturbation problem for the trivial invariant torus 
$\xi = 0$ of the system
\begin{equation*}
 \dot{\xi} = \varepsilon\mathcal{B}\xi,
 \quad
 \dot{\varphi} = \omega_{\ast}.
\end{equation*}
Note, that the quadratic forms $\langle D[r_{\ast}(0)]A\cdot, \cdot \rangle$ 
and $\langle V''(0)\Gamma\cdot, \cdot \rangle$ are positive and negative definite
 respectively, which means that the eigenvalues of $\mathcal{B}$ have non-zero 
real parts. The Lipschitz constants for the right-hand sides of 
\eqref{eq:dotxi}--\eqref{eq:dotphi} are $o(\varepsilon)$. Under such 
circumstances, one can easily verify that 
system~\eqref{eq:dotxi}--\eqref{eq:dotphi} satisfies all conditions of 
\cite[Lemma 2.1]{Hale61}
(There is a misprint in condition~(ii) on page~507: in 
$\mathrm{Lip}(\theta, y, z;\Sigma_{\sigma, \mu};\eta(\varepsilon, \sigma, \mu))$ 
the symbol $\eta$ should be replaced by $\lambda$.)
according to which, there exists sufficiently small $\varepsilon_0 > 0$,
such that for any $\varepsilon \in (0, \varepsilon_0)$
system~\eqref{eq:dotxi}--\eqref{eq:dotphi} has an invariant torus given
 by the equation $\xi = \tilde{\xi}(\varphi; \varepsilon)$, where
 $\tilde{\xi}(\cdot;\varepsilon) \colon \mathbb{T}^n \to \mathbb{R}^{n + m}$ 
is a Lipschitzian mapping, such that 
$\max_{\varphi \in \mathbb{T}^n} \| \xi(\varphi; \varepsilon) \| \to 0$ 
as $\varepsilon \to 0$. Therefore, we have shown that 
system~\eqref{eq:polar - r}--\eqref{eq:polar - phi} has an invariant torus
 $\mathcal{T}_{\varepsilon}^n$ located in $\mathcal{D}$.

Also, if $\xi(t)$ is a solution of 
system~\eqref{eq:polar - r}--\eqref{eq:polar - phi} corresponding to 
$(r(t), v(t))$, then
\begin{equation*}
 \| \xi(t) \|
  < \frac{
 c_{5}\varepsilon^{(N - 2) / 2 - \gamma}
  }{
 \varsigma
  }
 \quad
 \forall t > t^{\ast}.
\end{equation*}
Lemma~2.3 in~\cite{Hale61} claims that the trajectory corresponding to the 
solution $\xi(t)$ belongs to the stable invariant manifold of the invariant torus. 
This yields the inequality
\begin{equation*}
 \| \tilde{\xi}(\varphi; \varepsilon) \|
  \leq \frac{
 c_{5}\varepsilon^{(N - 2) / 2 - \gamma}
  }{
 \varsigma
  }
 \quad
 \forall\psi \in \mathbb{T}^n.
\end{equation*}

Furthermore, if $N \geq 5$ and $V''(0)$ is positive definite, then it was 
shown in \cite{SamParRep} that not only does the trajectory corresponding 
to the solution from Proposition~\ref{proposition:|r - r_8 * |} approach the
torus $\mathcal{T}_{\varepsilon}^n$ as $t\to\infty$, but it is also attracted 
by some trajectory on $\mathcal{T}_{\varepsilon}^n$ when $t\to\infty$.

Hence, we have proved the part of statements (3)--(5) from the main theorem
concerning $n$-dimensional tori.

\section{Ultimate behavior of solutions of type $(i_1, \ldots, i_{s})$}
\label{sec:ultimatei1..is}

Now, we will cover the case when solutions are of type $(1, \ldots, s)$. 
Obviously, results that we are going to obtain will be applicable to other 
possible types $(i_1, \ldots, i_{s})$, too.

We introduce new variables $p := (p_1, \ldots, p_{s}) \in \mathbb{R}_{+}^{s}$, 
$\vartheta := (\vartheta_1, \ldots, \vartheta_{s}) \in \mathbb{T}^{s}$ and
 $\zeta = ( \zeta_1, \ldots, \zeta_{n - s} ) \in \mathbb{C}^{n - s}$ 
via the formulae
\begin{equation*}
 z_{k} = \sqrt{p_{k}} \mathrm{e}^{\mathrm{i}\vartheta_{k}},
 \quad
 k \in \{ 1, \ldots, s\},
 \quad
 z_{k + s} = \varepsilon\zeta_{k},
 \quad
 k \in \{ 1, \ldots, n - s\}.
\end{equation*}
If we denote
\begin{gather*}
 \alpha'
  = ( \alpha_1, \ldots, \alpha_{s} ),
 \quad
 \alpha''
  = ( \alpha_{s + 1}, \ldots, \alpha_{n} ),
 \\
 \beta'
  = ( \beta_1, \ldots, \beta_{s} ),
\quad
 \beta''
  = ( \beta_{s + 1}, \ldots, \beta_{n} ),
 \\
 \omega_0'
  = ( \omega_{01}, \ldots, \omega_{0s} ),
 \quad
 \omega_0''
  = ( \omega_{0, s + 1}, \ldots, \omega_{0n} )
\end{gather*}
and decompose into blocks, the matrices
\begin{equation*}
 A := \begin{pmatrix}
  A_{11} & A_{12}
  \\
  A_{21} & A_{22}
 \end{pmatrix},
 \quad
 B := \begin{pmatrix}
  B_{11} & B_{12}
  \\
  B_{21} & B_{22}
 \end{pmatrix},
 \quad
 \Omega' = \begin{pmatrix}
  \Omega_1'
  \\
  \Omega_{2}'
 \end{pmatrix},
 \quad
 \Upsilon = \begin{pmatrix}
  \Upsilon_1 & \Upsilon_{2}
 \end{pmatrix},
\end{equation*}
where $\dim A_{11} = \dim B_{11} = s \times s$, $\dim \Omega_1' = s \times m$,
$\dim \Upsilon_1 = m \times s$,  we will get a smooth system of the form
\begin{gather}
 \label{eq:eq - polar - p}
\begin{aligned}
 \dot{p} &= 2\varepsilon D[p] \left[
  \alpha'
  - A_{11}p
  + \varepsilon\bar{a} \big(
 p, \overrightarrow{| \varepsilon\zeta |}^2, v
 \big)
  - \varepsilon^2 A_{12} \overrightarrow{| \zeta |}^2
 \right]
 \\
&\quad +  \varepsilon^{N / 2} D^{1 / 2}[p]
  \breve{a}(p, \operatorname{Re} \zeta, \operatorname{Im} \zeta, v, \vartheta),
\end{aligned} \\
 \label{eq:eq - v - 1}
\begin{aligned}
 \dot{v} &= \varepsilon\varsigma \big[
  \Upsilon_1p
  + \Gamma v
  + \varepsilon\bar{g} (
 p, \overrightarrow{| \varepsilon\zeta |}^2, v
  )
  + \varepsilon^2 \Upsilon_{2} \overrightarrow{| \zeta |}^2
 \\
&\quad + \varepsilon^{(N + 1) / 2} \breve{g}(p, \operatorname{Re}\zeta, \operatorname{Im} \zeta, v, \vartheta)
 \big],
\end{aligned} \\
 \nonumber
\begin{aligned}
 \dot{\zeta} &= D \big[
  \mathrm{i} (
 \omega_0''
 + \varepsilon (
  \beta''
  + \Omega_{2}'v
  - B_{21}p
 )
  )
  + \varepsilon (
 \alpha''
 - A_{21}p
  )
 \big] \zeta
 \\
&\quad   + \varepsilon^2D \big[
  - (
 \mathrm{i}B_{22} + A_{22}
  )
  \overrightarrow{| \zeta |}^2
  + \bar{h} (
 p, \overrightarrow{| \varepsilon\zeta |}^2, v
  )
 \big] \zeta
 + \varepsilon^{(N - 2) / 2}
  \breve{h}(p, \operatorname{Re}\zeta, \operatorname{Im} \zeta, v, \vartheta),
\end{aligned} \\
 \nonumber
\begin{aligned}
 \dot{\vartheta} &= \omega_0'
  + \varepsilon (
 \beta'
 + \Omega_1'v
 - B_{11}p
  )
  +\varepsilon^2 \big[
  - B_{12} \overrightarrow{| \zeta |}^2
  + \bar{b} (
 p, \overrightarrow{| \varepsilon\zeta |}^2, v
  )
 \big] \\
&\quad + \varepsilon^{N / 2} D^{-1 / 2}[p]
  \breve{b}(p, \operatorname{Re}\zeta, \operatorname{Im} \zeta, v, \vartheta),
\end{aligned}
\end{gather}
whose domain contains the set 
$\operatorname{cl} (\mathcal{D}_1) \times \{ \zeta \in \mathbb{C}^{n - s}
 \colon \| \zeta \| \leq 1\} \times \mathbb{T}^{s}$, 
with $\mathcal{D}_1 \subset \mathbb{R}_{+}^{s} \times \mathbb{R}^{m}$ 
being defined by the inequalities
\[
 p_{k} > \frac{\alpha_{k}}{2a_{kk}}, \quad
  k \in \{ 1, \ldots, s\},
 \quad
 | p | < \frac{
  2\max_{1 \leq k \leq s}\alpha_{k}
 }{
  a_{+}^0
 },
 \quad
 \| v \| < R.
\]

Let $( p(t), v(t), \zeta(t), \vartheta(t) )$ be a solution of this system 
which corresponds to a solution of type $(1, \ldots, s)$, so that
\begin{equation*}
 ( p(t), v(t), \zeta(t) )
  \in \mathcal{D}_1 \times \{
 \zeta \in \mathbb{C}^{n - s}
 \colon
 \| \zeta \| \leq 1
  \}
 \quad
 \forall t > t_{\ast}
\end{equation*}
if $t_{\ast} > 0$ is large enough.

The following is an analogue of Proposition~\ref{proposition:r - * v - *}.

\begin{proposition}
For sufficiently small $\varepsilon_0 > 0$ and for any
$\varepsilon \in (0, \varepsilon_0]$, the system
\begin{equation*}
 \alpha' - A_{11}p + \varepsilon\bar{a} ( p, 0, v ) = 0,
 \quad
 \Upsilon_1p + \Gamma v + \varepsilon\bar{g} ( p, 0, v )
\end{equation*}
has a solution
\begin{equation*}
 p = \bar{p}_{\ast}(\varepsilon)
  := A_{11}^{-1} \alpha' + O(\varepsilon),
 \quad
 v = \bar{v}_{\ast}(\varepsilon)
  := - \Gamma^{-1} \Upsilon_1 A_{11}^{-1} \alpha'
  + O(\varepsilon),
\end{equation*}
such that
\begin{equation*}
 \bar{p}_{\ast k}(\varepsilon)
  > \frac{2\alpha_{k}}{3a_{kk}},
 \quad
 k \in \{ 1, \ldots, s\},
 \quad
 | \bar{p}_{\ast}(\varepsilon) |
  < \frac{
 3 \max_{1 \leq k \leq n} \alpha_{s}
  }{
 2a_{+}^0
  },
 \quad
 \| \bar{v}_{\ast}(\varepsilon) \| < \frac{2R}{3}.
\end{equation*}
\end{proposition}

Note that the equilibrium $(\bar{p}_{\ast}, \bar{v}_{\ast})$ of the system
\begin{gather*}
 \dot{p} = 2\varepsilon D[p] [
  \alpha'
  - A_{11}p
  + \varepsilon\bar{a} ( p, 0, v )
 ],
 \\
 \dot{v} = \varepsilon\varsigma [
  \Upsilon_1p
  + \Gamma v
  + \varepsilon\bar{g} ( p, 0, v )
]
\end{gather*}
is the unique stationary point in $\mathcal{D}_1$ of the Morse function
\begin{equation*}
 \bar{W}(p, v) := \sum_{i = 1}^{s} \Big(
  p_i + \bar{p}_{\ast i} \ln (
 \frac{\bar{p}_{\ast i}}{p_i}
  ) - \bar{p}_{\ast i}
 \Big) + \frac{1}{2} \langle
  V''(0) (v - \bar{v}_{\ast}), v - \bar{v}_{\ast}
 \rangle.
\end{equation*}

\begin{proposition}
There are positive numbers $c_{7}$, $\varsigma_0$ and $\varepsilon_0$,
such that for any $\varsigma \in ( 0, \varsigma_0 )$, 
$\varepsilon \in (0, \varepsilon_0)$ the following assertion is valid. 
If $( z(t), v(t) )$ is a solution of type $( 1, \ldots, s )$ of 
system~\eqref{eq:vect - norm - form - z}--\eqref{eq:vect - norm - form - v} and
\begin{equation*}
 p(t) := \big(  | z_1(t) |^2, \ldots, | z_{s}(t) |^2
 \big),
 \quad
 \zeta(t) := \frac{1}{\varepsilon} ( z_{s + 1}(t), \ldots, z_{n}(t) ),
\end{equation*}
then there exists  $t^{\ast} > 0$, such that
\begin{equation*}
 \sqrt{
  \| p(t) - \bar{p}_{\ast}(\varepsilon) \|^2
  + \| v(t) - \bar{v}_{\ast}(\varepsilon) \|^2
 } < \frac{
  c_{7}\varepsilon^{(N - 2) / 2}
 }{
  \varsigma
 },
 \quad
 \| \zeta(t) \|
  \leq c_{7} \varepsilon^{(N - 4) / 2}
 \quad
 \forall t > t^{\ast}.
\end{equation*}
\end{proposition}

\begin{proof}
Since all elements of the matrix $A_{21}$ are non-positive, then the minimal 
element of $\alpha'' - A_{21}p$ is not less than the minimal element of $\alpha''$. 
Thus, there exists a constant $c_{8} > 0$, such that for all sufficiently small 
$\varepsilon > 0$ there hold the inequalities
\begin{equation*}
 \frac{\mathrm{d}}{\mathrm{d}t} \| \zeta(t) \| ^2
  \geq \| \zeta(t) \| \big[
 \varepsilon \min_{s + 1 \leq k \leq n}
  \alpha_{k} \| \zeta(t) \|
 - c_{8} \varepsilon^{(N - 2) / 2}
  \big],
 \quad
 \| \zeta(t) \| \leq 1  \quad
 \forall t \geq t^{\ast},
\end{equation*}
and consequently,
\begin{equation*}
 \| \zeta(t) \|
  \leq \frac{
 c_{8} \varepsilon^{(N - 4) / 2}
  }{
 \min_{s + 1 \leq k \leq n} \alpha_{k}
  }
\end{equation*}
for all $t \geq t^{\ast}$.

Now observe that for
\begin{equation*}
 \| \zeta \|
  \leq \frac{
 c_{8} \varepsilon^{(N - 3) / 2}
  }{
 \min_{s + 1 \leq k \leq n} \alpha_{k}
  }
\end{equation*}
sub-system~\eqref{eq:eq - polar - p}--\eqref{eq:eq - v - 1} takes the form
\begin{gather*}
 \dot{p} = 2\varepsilon D[p] \left[
  \alpha'
  - A_{11}p
  + \varepsilon\bar{a} ( p, 0, v )
 \right] + O ( \varepsilon^{N / 2} ),
 \\
 \dot{v} = \varepsilon\varsigma \left[
  \Upsilon_1p
  + \Gamma v
  + \varepsilon\bar{g} ( p, 0, v )
 \right] + O ( \varepsilon^{(N + 3) / 2} ).
\end{gather*}
Using Lemma~\ref{lemma:W'f < 0 - 2} for the function $\bar{W}(p, v)$,
one can obtain the inequality for $p(t)$ and $v(t)$ precisely in the 
same way as we did in Proposition~\ref{proposition:|r - r_8 * |}.
\end{proof}

If we introduce new variables $\eta \in \mathbb{R}^{s + m}$ and 
$\xi \in \mathbb{R}^{s + m} \times \mathbb{C}^{n - s}$ via the formulae
\begin{equation*}
 (p, v) = (\bar{p}_{\ast}, \bar{v}_{\ast}) + \varepsilon^{\sigma}\eta,
\end{equation*}
then again we come to the perturbation problem for the trivial hyperbolic
 invariant torus of the system
\begin{gather*}
 \dot{\eta} = \varepsilon \begin{pmatrix}
  - 2D[\bar{p}_{\ast}(0)] A_{11} & 0
  \\
  \varsigma\Upsilon_1 & \varsigma\Gamma
 \end{pmatrix} \eta,
 \\
 \dot{\zeta} = D \left[
  \mathrm{i} (
 \omega_0''
 + \varepsilon (
  \beta''
  + \Omega_{2}'\bar{v}_{\ast}(0)
  - B_{21} \bar{p}_{\ast}(0)
 )
  )
  + \varepsilon(\alpha''
  - A_{12} \bar{p}_{\ast}(0))
 \right] \zeta,
 \\
 \dot{\vartheta} = \omega_0' + \varepsilon \left[
  \beta'
  + \Omega_1'\bar{v}_{\ast}(0)
  - B_{11} \bar{p}_{\ast}(0)
 \right].
\end{gather*}
The real parts of the eigenvalues of this system's matrix are non-zero, 
since they are defined by the eigenvalues of the matrix
\begin{equation*}
 \varepsilon\operatorname{diag}(
  - 2D[\bar{p}_{\ast}(0)] A_{11},
  \,
  \varsigma\Gamma,
  \,
  \alpha'' - A_{12} \bar{p}_{\ast}(0)
 ).
\end{equation*}
It is not hard to see that the perturbation terms satisfy the conditions 
of \cite[Lemma 2.1]{Hale61}. Using the same arguments as above for the
 realification of the perturbed system in the variables 
$\eta, \operatorname{Re}\zeta, \operatorname{Im} \zeta, \vartheta$, we can prove the existence of an
$s$-dimensional truly hyperbolic invariant torus 
$\mathcal{T}_{\varepsilon}^{s}(1, \ldots, s)$, which attracts the 
trajectories corresponding to the solutions of type $(1, \ldots, s)$. 
In the coordinates $(p, v, \zeta, \vartheta)$ this torus is given by the equations
\begin{equation*}
 (p, v) = (\bar{p}_{\ast}(\varepsilon), \bar{v}_{\ast}(\varepsilon))
  + \varepsilon^{\sigma} \eta_{\ast}(\vartheta; \varepsilon),
 \quad
 \zeta = \zeta_{\ast}(\vartheta; \varepsilon),
\end{equation*}
where $\eta_{\ast}(\cdot; \varepsilon) \colon \mathbb{T}^{s} \to \mathbb{R}^{s + m}$ 
and $\zeta_{\ast}(\cdot; \varepsilon) \colon \mathbb{T}^{s} \to \mathbb{C}^{n - s}$
 are Lipschitzian mappings satisfying the conditions
\begin{equation*}
 \max_{\vartheta \in \mathbb{T}^{s}}
  \| \eta_{\ast}(\vartheta; \varepsilon) \|
  \leq \frac{c_{8}\varepsilon^{(N - 2) / 2 - \gamma}}{\varsigma},
 \quad
 \max_{\vartheta \in \mathbb{T}^{s}}
  \| \zeta_{\ast}(\vartheta; \varepsilon) \|
  \leq c_{8}\varepsilon^{(N - 4) / 2}.
\end{equation*}

Hence, we have proved the part of statements (3)--(5) from the main 
theorem concerning the tori of dimensions less then $n$. Note, that those forward
semi-trajectories in $\mathfrak{A}$ that are not attracted by stable $n$-dimensional 
tori lie on stable manifolds of truly hyperbolic tori, and thus form the set 
of zero Lebesgue measure.

\section{Excitation of two-frequency oscillations in a system of two coupled 
generators due to slow cooling}

Let us now provide an example of a hypothetical device where the described 
phenomenon can be observed. Consider a system of two coupled generators 
(Figure~\ref{fig:1}).

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.8\textwidth]{fig1} % circuit
\end{center}
 \caption{Coupled generators}
 \label{fig:1}
\end{figure}
Here, the $i$-th generator consists of the following elements: a thermistor
$R_i$, a magneto resistor $\rho_i$, magnetically connected inductors
$L_i$, $L_i'$ having a negative mutual induction $M_i=\gamma_i\sqrt{L_iL_i'}$, 
a capacitor $C_i$, and an
active feedback element $A_i$. We suppose that the I-V characteristics $\chi_i$
of $A_i$ is a smooth function of voltage difference $V_i$ on $L_i'$ and admits
the expansion
\begin{equation*}\label{eq:chi_i=00003D}
 \chi_i = \chi_i ( V_i)
  := \chi_{i1} V_i - \chi_{i3} V_i^{3}
  + O ( V_i^{4} ),
 \quad
 \chi_{i1} > 0,
 \quad
 \chi_{i3} > 0.
\end{equation*}
We assume that the $i$-th thermistor has a positive temperature coefficient
and that its resistance depends on the thermistor's temperature $T_i$ 
by the relation
\begin{equation*}\label{eq:R_i=00003D}
 R_i = R_i(T_i)
  := \varepsilon R_{i0} + \mu R_{i1}T_i,
 \quad
 i \in \{ 1, 2 \},
\end{equation*}
where $R_{i0}$, $R_{i1}$ are positive constants and $\varepsilon$
and $\mu$ are small parameters. We also suppose the inductances $L_i$, 
$L_i'$ and the capacities $C_i$ to be constant.

The resistance of the magneto resistor depends on an external magnetic field.
If the latter is orthogonal to the direction of the current in the resistor,
then the change of its resistance is approximately proportional to the
square of the magnetic field magnitude~\cite{Kir78}.

The generators interact with each other in the following
way. The current $I_i$ in the inductor $L_i$ produces a magnetic field of
a magnitude proportional to $I_i$, and this magnetic field influences
the resistance of the magneto resistor $\rho_{j}$. Thus, it is natural
to consider the case where
\begin{equation}\label{eq:rho_i}
 \rho_{j} = \rho_{j}(I_i) = \rho_{j0} + \rho_{j2} ( L_i I_i )^2,
 \quad
 \rho_{j0} > 0,
 \quad \rho_{j2} > 0,
 \quad  i \in \{ 1, 2 \}.
\end{equation}
Taking into account the Newton law of cooling and the Joule --
Lenz law of ohmic heating and assuming the environment temperature
to be zero, we adopt the next equation for the resistor $R_i$
temperature change
\begin{equation}\label{eq:temper}
 \dot{T_i} = -\mu k_iT_i + K_i I_i^2 R_i,
\end{equation}
where $k_i$ and $K_i$ are some positive coefficients.

Our goal is to show that under an appropriate choice of the generators
parameters one can observe such a phenomenon. If the described
device of the coupled generators with the sufficiently high initial temperatures
$T_i(0)$ of the resistances $R_i$ is placed into the environment, then
first for a long period of time this device remains in the sleep mode.
However, when the resistances temperatures drop almost to zero,
the device wakes up and after a transient process, in general, it starts producing
two-frequency oscillations. Such a kind of behavior can be treated as a 
phenomenon of dynamical bifurcation.

Denote by $q_i$ the charge of the capacitor $C_i$, and let 
$I_i, I_{C_i}, I_{\rho_i}$ and $I_{A_i}$ stand for the currents through the 
resistor $R_i$, the capacitor $C_i$, the magneto resistor $\rho_i$ and the active
element $A_i$ respectively. The Kirchhoff laws yield that
\begin{gather}
 \label{eq:Kirch-1}
 I_i + I_{C_i} + I_{\rho_i} = I_{A_i} = \chi_i,
 \\
 \label{eq:Kirch-2}
 L_i \dot{I}_i + R_i I_i = \frac{q_i}{C_i} = \rho_iI_{\rho_i}.
\end{gather}
Since $\dot{q}_i=I_{C_i}$, by differentiation of~(\ref{eq:Kirch-2}) we obtain
\begin{equation*}
 L_i \ddot{I}_{R_i} + R_i \dot{I}_i + \dot{R}_i I_i
 = \frac{I_{C_i}}{C_i}.
\end{equation*}
But from~(\ref{eq:Kirch-2}) and~(\ref{eq:Kirch-1}) we can find
\begin{gather*}
 I_{\rho_i} = [
  L_i\dot{I}_i + R_iI_i ]/\rho_i,
 \\
 I_{C_i} = \chi_i - I_i - \frac{1}{\rho_i} [
  L_i\dot{I}_i + R_iI_i ].
\end{gather*}
Hence,
\begin{equation*}
 \ddot{I}_i + \frac{R_i}{L_i} \dot{I}_i + \frac{I_i}{L_iC_i} 
+ \frac{1}{\rho_iL_iC_i} [ L_i\dot{I}_i + R_iI_i ] 
= \frac{\chi_i}{L_iC_i} - \frac{\dot{R}_i}{L_i} I_i,
\end{equation*}
and on account of~\eqref{eq:temper}
\begin{equation*}
 \ddot{I}_i + [
  \frac{R_i}{L_i} + \frac{1}{\rho_iC_i}
 ] \dot{I}_i + \frac{1}{L_iC_i} [
  1 + \frac{R_i}{\rho_i} ] I_i 
= \frac{\chi_i}{L_iC_i} - \frac{\mu R_{i1}}{L_i} [
  -\mu k_iT_i + K_iI_i^2 R_i
] I_i.
\end{equation*}
Introducing new variables $I_i=\sqrt{\rho_{j0}}w_i$, 
$T_i=R_{i0}u_i/(\varsigma R_{i1})$, imposing the constraint 
$\varsigma\varepsilon=\mu$ with the constant
$\varsigma$ playing the same role as in section~\ref{sec:preliminary-results} 
and taking into account~\eqref{eq:rho_i} and~\eqref{eq:temper}, we finally 
reach the system
\begin{gather*}
\begin{aligned}
& \ddot{w}_i + \Big[
  \frac{\varepsilon R_{i0}(1+u_i)}{L_i}
  +\frac{1}{C_i\rho_{i0}(1+\rho_{i2}L_{j}^2w_{j}^2)}
 \Big] \dot{w}_i
 +\frac{1}{L_iC_i} \Big[
  1 + \frac{\varepsilon R_{i0}(1+u_i)}{\rho_{i0} (
 1+\rho_{i2}L_{j}^2w_{j}^2  )} \Big] w_i  \\
&=  \frac{1}{L_iC_i\sqrt{\rho_{j0}}} \chi_i(
  M_i\sqrt{\rho_{j0}}\dot{w}_i )
 -\frac{\mu\varepsilon R_{i0}}{L_i}[
  - k_i u_i + K_i R_{i1} (
 1 + u_i
  ) \rho_{j0} w_i^2
] w_i,
\end{aligned} \\
 \dot{u_i} = \mu [
  - k_i u_i + R_{i1}K_i \rho_{j0} (
 1 + u_i
  ) w_i^2],
 \quad
 i \in \{ 1, 2 \},
 j \in \{ 1, 2 \},
 \quad
 i \ne j.
\end{gather*}
This system can be represented in the form \eqref{eq:main1} if the
 parameters of the generators are chosen in such a way, that
\begin{equation*}
 \frac{\chi_{i1}M_i}{L_iC_i}
 - \frac{1}{C_i\rho_{i0}}
 = \varepsilon b_i,
 \quad 
 i\in\{ 1,2\},
\end{equation*}
where $b_1,$ $b_{2}$ are positive constants satisfying the inequalities
\begin{equation}\label{eq:b_i}
 b_i > R_{i0} / L_i,
 \quad
 i\in\{ 1,2\}.
\end{equation}
In fact, we have
\begin{equation*}
 \omega_{0i}(u) = \sqrt{\frac{1}{L_iC_i}}+O(\varepsilon),
 \quad
 \lambda_i(u) = b_i - \frac{R_{i0}(1+u_i)}{L_i},
 \quad
 i\in\{ 1,2\},
\end{equation*}
and
\begin{align*}
& F_i(w_1,w_{2},\dot{w}_1,\dot{w}_{2},u,0)
 \\
&  = \sum_{k+l+m+n=3}
 f_{i,jklm}w_1^{j} \dot{w}_1^{k} w_{2}^{l} \dot{w}_{2}^{m}
  + O(\| w \|^{4} + \| \dot{w} \|^{4}),
 \quad
 i \in \{ 1,2\},
\end{align*}
where the only non-zero coefficients are
\begin{gather*}
f_{1,0300} = -\frac{\chi_{13}\rho_{20}M_1^{3}}{L_iC_i},
 \quad
 f_{1,0120} = \frac{\rho_{12}L_1^2}{C_1\rho_{10}},
 \\
 f_{2,2001} = \frac{\rho_{22}L_{2}^2}{C_{2}\rho_{20}},
 \quad
 f_{2,0003} = -\frac{\chi_{23}\rho_{10}M_{2}^{3}}{L_iC_i}.
\end{gather*}
Therefore, performing the change of variables~\eqref{eq:new-var-x}, we obtain
\begin{equation*}
 \hat{F}(x,u,\varsigma\varepsilon) = \begin{pmatrix}
  0
  \\
  \frac{1}{\omega_{01}(u)} F_1 (
 x_1,x_{3},\omega_{01}(u)x_{2},\omega_{02}(u)x_{4},u,0
  )
  \\
  0
  \\
  \frac{1}{\omega_{02}(u)} F_{2} (
 x_1,x_{3},\omega_{01}(u)x_{2},\omega_{02}(u)x_{4},u,0
  )
 \end{pmatrix} + O(\varepsilon).
\end{equation*}

It is easily seen the eigenvectors of the matrix $J(u)$ are
\begin{equation*}
 s_1^{+} = \begin{pmatrix}
  -\mathrm{i}\\
  1\\
  0\\
  0
 \end{pmatrix},
 \quad
 s_{2}^{+} = \begin{pmatrix}
  0\\
  0\\
  -\mathrm{i}\\
  1
 \end{pmatrix},
 \quad
 s_1^{-} = \begin{pmatrix}
  \mathrm{i}\\
  1\\
  0\\
  0
 \end{pmatrix},
 \quad
 s_{2}^{-} = \begin{pmatrix}
  0\\
  0\\
  \mathrm{i}\\
  1
 \end{pmatrix}.
\end{equation*}
If we now introduce new complex variables $z_1,z_{2}\in\mathbb{C}$ via
\begin{equation*}
 x_1=-\mathrm{i}z_1+\mathrm{i}\bar{z}_1,\quad x_{2}=z_1+\bar{z}_1,
 \quad
 x_{3}=-\mathrm{i}z_{2}+\mathrm{i}\bar{z}_{2},\quad x_{4}=z_{2}+\bar{z}_{2},
\end{equation*}
we will be able to find the elements of the matrix $A(u)$ by extraction
 of the resonant terms from cubic nonlinearities:
\begin{align*}
& \frac{1}{2 \omega_{01}(u)}
 F_1 (
  -\mathrm{i}z_1+\mathrm{i}\bar{z}_1,
  -\mathrm{i}z_{2}+\mathrm{i}\bar{z}_{2},
  \omega_{01}(u)(z_1+\bar{z}_1),
  \omega_{02}(u)(z_{2}+\bar{z}_{2}),u,0
 )
 \\
&=- a_{11}(u) |z_1|^2 z_1
 - a_{12}(u) |z_{2}|^2 z_1
 + [\text{nonresonant terms}]
\\
&\frac{1}{2 \omega_{02}(u)}
 F_{2} (
  -\mathrm{i}z_1+\mathrm{i}\bar{z}_1,
  -\mathrm{i}z_{2}+\mathrm{i}\bar{z}_{2},
  \omega_{01}(u)(z_1+\bar{z}_1),
  \omega_{02}(u)(z_{2}+\bar{z}_{2}),u,0
  )
 \\
 &=- a_{21}(u)|z_1|^2z_{2}
 - a_{22}(u) |z_{2}|^2z_{2}
 + [\text{nonresonant terms}].
\end{align*}
It turns out that
\begin{gather*}
 a_{11}(u)
  = - \frac{1}{2} f_{1,2100} - \frac{3}{2} \omega_{01}^2(u) f_{1,0300},
 \\
 \nonumber
 a_{12}(u)
  = - f_{1,0120} - \omega_{02}^2(u)f_{1,0102},
 \\
 a_{21}(u)
  = - f_{2,2001} - \omega_{01}^2(u)f_{2,0201},
 \\
 a_{22}(u)
  = - \frac{1}{2} f_{2,0021} - \frac{3}{2} \omega_{02}^2(u) f_{2,0003}.
\end{gather*}
In our case, when $\varepsilon=0$, these elements does not depend on $u$:
\begin{gather*}
 a_{11} = \frac{3\chi_{13}\rho_{20}M_1^{3}}{2(C_1L_1)^2},
 \quad
 a_{12} = -\frac{\rho_{12}L_{2}^2}{C_1\rho_{10}},
 \\
 a_{21} = -\frac{\rho_{22}L_1^2}{C_{2}\rho_{20}},
 \quad
 a_{22} = \frac{3\chi_{23}\rho_{10}M_{2}^{3}}{2(C_{2}L_{2})^2}.
\end{gather*}

Hence, when $\varepsilon=0$, the positive definiteness condition of the 
symmetric part of the matrix $A$ takes the form
\begin{equation}\label{eq:posdepA}
 9\chi_{13}\chi_{23} (
  \rho_{10}\rho_{20} M_1M_{2}
 )^{3} > \left[
  L_1L_{2} (
 \rho_{12}\rho_{20} C_{2}L_{2}^2
 + \rho_{10}\rho_{22} C_1L_1^2
  )
 \right]^2.
\end{equation}

Since $g(u) = ( -k_1 u_1,- k_{2} u_{2} )$, the Morse function can be chosen 
as $V(u) = u_1^2 + u_{2}^2$. It has a unique stationary point $u_{\ast}=(0,0)$. 
The instability and stability sets are defined as follows
\begin{gather*}
\mathcal{V}_{+}= \{ (u_1,u_{2}) \in\mathbb{R}^2
  \colon
  u_1^2 + u_{2}^2 < V^{\ast}, u_i < b_iL_{i0} / R_{i0} - 1,
  i \in \{ 1,2 \} \},
 \\
 \mathcal{V}_{-}= \{ (u_1,u_{2}) \in \mathbb{R}^2
  \colon
  u_1^2 + u_{2}^2 < V^{\ast}, u_i > b_iL_{i0} / R_{i0} - 1,
  i \in \{ 1,2 \} \},
\end{gather*}
where $V^{\ast} > 0$ is large enough.

Thus, if the numbers $\sqrt{C_1L_1}$, $\sqrt{C_{2}L_{2}}$ are rationally 
independent and conditions~\eqref{eq:b_i} and~\eqref{eq:posdepA} fulfill, 
then hypotheses~(H1)--(H5) hold, and Theorem~\ref{thm:main}
implies that under the appropriate choice of $\varsigma_0>0$, 
$0<\varsigma_{\ast}<\varsigma_0$, $\varepsilon_0>0$ the aforementioned 
changes in behavior of the coupled generators can actually be observed, 
once $\varepsilon\in (0,\varepsilon_0)$, 
$\mu\in (\varsigma_{\ast}\varepsilon,\varsigma_0\varepsilon)$.


\section{Addendum}
\label{sec:addendum}

Consider the formal system
\begin{equation}\label{eq:formal_eps}
\begin{gathered}
 \dot{x} = J(u)x + \sum_{i \geq 0} \mu^{i} F_i(x, u),
 \\
 \dot{u} = \mu \Big[
  g(u) + \sum_{i \geq 0} \mu^{i} G_i(x, u)
 \Big]
\end{gathered}
\end{equation}
obtained from~\eqref{eq:main1} by expanding its right-hand sides into the 
Taylor series expansions in powers of $\mu$. Our goal is to simplify this 
system with aid of the formal change of variables
\begin{equation}\label{eq:transform}
 x = y + \sum_{i \geq 0} \mu^{i}X_i(y, v),
 \quad
 u = v + \mu\sum_{i \geq 0} \mu^{i}U_i(y, v).
\end{equation}

Construction of such a transformation consists in solving homological equations 
of the form
\begin{gather} \label{eq:homol_equ_Y}
 \mathfrak{L}_{J(v)y} [  \mathcal{Y}_{j}(v) y^{j}] 
= \mathcal{P}_{j}(v)y^{j} - \mathcal{R}_{j}(v)y^{j},
 \\
 \label{eq:homol_equ_Z}
 \partial_{J(v)y} [
  \mathcal{Z}_{j}(v) y^{j}] 
= \mathcal{Q}_{j}(v)y^{j} - \mathcal{S}_{j}(v)y^{j}.
\end{gather}
Here $\mathcal{P}_{j}(v)y^{j}$ and $\mathcal{Q}_{j}(v)y^{j}$ are known $j$-th 
order homogeneous forms which take values in $\mathbb{R}^{2n}$ and 
$\mathbb{R}^{m}$ respectively and smoothly depend on the parameter $v$, 
whereas $\partial_{J(v)y}$ and $\mathfrak{L}_{J(v)y}$ are, respectively, 
the directional and the Lie derivatives along the vector field $J(v)y$. Namely,
\begin{gather*}
 \partial_{J(v)y}Z(y) := \frac{\partial Z(y)}{\partial y} J(v)y
 \quad
 \forall Z(\cdot) \in \mathrm{C}^{1}
  ( \mathbb{R}^{2n}; \mathbb{R}^{m} ),
 \\
 \mathcal{L}_{J(v)y}Y(y)
  := \frac{\partial Y(y)}{\partial y} J(v)y - J(v)Y (y)
 \quad
 \forall Y(\cdot) \in \mathrm{C}^{1}
  ( \mathbb{R}^{2n}; \mathbb{R}^{2n} ).
\end{gather*}
The forms $\mathcal{Y}_{j}(v), \mathcal{R}_{j}(v)y^{j}$, $\mathcal{Z}_{j}(v)$ 
and $\mathcal{S}_{j}(v)y^{j}$ are determined in such a way, that they smoothly 
depend on $v$ and satisfy the corresponding equations.

If $j = 0$, equation~\eqref{eq:homol_equ_Y} has obvious solutions
\begin{equation}\label{eq:R0 = 00003D0}
 \mathcal{R}_0(v) \equiv 0,
 \quad
 \mathcal{Y}_0(v) = - J^{-1}(v)\mathcal{P}_0(v),
 \quad
 \mathcal{S}_0(v) = \mathcal{Q}_0(v),
 \quad
 \mathcal{Z}_0(v) = 0.
\end{equation}
In the case when $j \geq 1$, in the same manner as in \cite{SamParRep}, 
we can introduce a suitable basis in the space of vector-valued polynomial forms. 
Note, that the matrix $J(v)$ has constant linearly independent eigenvectors 
$s_{j}^{\pm} \in \mathbb{C}^{2n}$, $j = 1, \ldots, n$, such that
\begin{equation*}
 J(v) s_{j}^{\pm} = [
  \varepsilon\lambda_{j}(v) \pm \mathrm{i}\omega_{j}(v)
 ] s_{j}^{\pm},
\end{equation*}
with vectors $s_{j}^{-}$ and $s_{j}^{+}$ being complex conjugate for all 
$j = 1, \dots, n$. Denote by $S$ the matrix with the columns 
$s_1^{+}, \ldots, s_{n}^{+}, s_1^{-}, \ldots, s_{n}^{-}$, and define 
the homogeneous forms
\begin{equation*}
 \mathfrak{s}_{\mathbf{q}}(y)
  := [S^{-1}y]^{\mathbf{q}},
 \quad
 e_{j, \mathbf{q}}^{\pm}(y)
  = \mathfrak{s}_{\mathbf{q}}(y)s_{j}^{\pm}, \label{eq:basform}
\end{equation*}
where $\mathbf{q} := (q_1, \ldots, q_{2n}) \in \mathbb{Z}_{+}^{2n}$ and
$x^{\mathbf{q}} := x_1^{q_1}\cdots x_{2n}^{q_{2n}}$ for $x = (x_1, \ldots, x_{2n})$.
This gives us the following expansions.
\begin{gather*}
 \mathcal{P}_{j}(v) y^{j}
  = \sum_{k = 1}^n \sum_{| \mathbf{q} | = j}
 [
  P_{k, \mathbf{q}}^{+}(v) e_{k, \mathbf{q}}^{+}(y)
  + P_{k, \mathbf{q}}^{-}(v) e_{k, \mathbf{q}}^{-}(y)
 ],
 \quad
 \mathcal{Q}_{j}(v)y^{j}
  = \sum_{| \mathbf{q} | = j}
 \mathfrak{s}_{\mathbf{q}}(y) Q_{\mathbf{q}}(v),
 \\
 \mathcal{R}_{j}(v) y^{j}
  = \sum_{k = 1}^n \sum_{| \mathbf{q} | = j}
 [
  R_{k, \mathbf{q}}^{+}(v) e_{k, \mathbf{q}}^{+}(y)
  + R_{k, \mathbf{q}}^{-}(v) e_{k, \mathbf{q}}^{-}(y)
 ],
 \quad
 \mathcal{S}_{j}(v) y^{j}
  = \sum_{| \mathbf{q} | = j}
 \mathfrak{s}_{\mathbf{q}}(y)S_{\mathbf{q}}(v).
\end{gather*}

\begin{proposition}\label{proposition:homol_equ}
Let $\mathcal{D}_{N} \subset \mathbb{R}^{m}$ be such a domain, that for any 
$k \in \{ 1, \ldots, n\} $ and $\sigma \in \{ 0, 1\}$ the equality
\begin{equation*}
 \min_{v \in \operatorname{cl} (\mathcal{D}_{N})}
  \big|
 \sum_{l = 1}^n (q_{l} - q_{l + n}
 - \sigma\delta_{kl}) \omega_{0l}(v)
  \big| = 0,
\end{equation*}
where $\mathbf{q} \in \mathbb{Z}_{+}^{2n}$, 
$1 \leq | \mathbf{q} | := \sum_{k = 1}^{2n}q_{k} \leq N$ and 
$\delta_{kl}$ is Kronecker's delta, holds if and only if
\begin{equation*}
 q_{l} - q_{l + n} - \sigma\delta_{kl} = 0
 \quad  \forall l \in \{1, \ldots, n \}.
\end{equation*}
Suppose also that the forms $\mathcal{P}_{j}(v)y^{j}$ and $\mathcal{Q}_{j}(v)y^{j}$ 
smoothly depend on $v \in \operatorname{cl} (\mathcal{D}_{N})$ and assign
\begin{gather*}
 R_{k, \mathbf{q}}^{\pm}(v)
  = \begin{cases}
 P_{k, \mathbf{q}}^{\pm}(v)
 & \text{if }
 \sum_{l = 1}^n
  | q_{l} - q_{l + n} \mp \delta_{kl} | = 0,
 \\
 0
 & \text{if }
 \sum_{l = 1}^n
  | q_{l} - q_{l + n}\mp\delta_{kl} | \neq 0,
  \end{cases}
\\
 S_{\mathbf{q}}(v)
  = \begin{cases}
 Q_{\mathbf{q}}(v)
 & \text{if }
 \sum_{l = 1}^n
  | q_{l} - q_{l + n} | = 0,
 \\
 0
 & \text{if }
 \sum_{l = 1}^n
  | q_{l} - q_{l + n} | \neq 0.
  \end{cases}
\end{gather*}
Then, for sufficiently small $\varepsilon_0 > 0$ and for all 
$\varepsilon \in [0, \varepsilon_0]$, there exist an $\mathbb{R}^{2n}$-valued form 
$\mathcal{Y}_{j}(v)y^{j}$
and an $\mathbb{R}^{m}$-valued form $\mathcal{Z}_{j}(v)y^{j}$, which satisfy 
equations~ \eqref{eq:homol_equ_Y} and \eqref{eq:homol_equ_Z} respectively. 
The coefficients of these forms are smooth functions in
 $\operatorname{cl} (\mathcal{D}_{N})$.
\end{proposition}

\begin{proof}
Since
\begin{align*}
 S^{-1}J(v)S
& = \operatorname{diag}\big[
  \varepsilon \lambda_1(v)  + \mathrm{i}\omega_1(v),
  \ldots,
  \varepsilon \lambda_{n}(v)  + \mathrm{i}\omega_{n}(v),
  \varepsilon \lambda_1(v) \\
&\quad - \mathrm{i}\omega_1(v),
  \ldots,
  \varepsilon \lambda_{n}(v) - \mathrm{i}\omega_{n}(v)
 \big],
\end{align*}
we have
\begin{align*}
 \partial_{J(v)y} \mathfrak{s}_{\mathbf{q}}(y)
&=  \frac{\mathrm{d}}{\mathrm{d}t} \big|_{t = 0} [
 (
  S^{-1} \mathrm{e}^{J(v)t}y
 )^{\mathbf{q}}
 ]
 \\
& =\Big[
   \varepsilon \sum_{l = 1}^n (q_{l} + q_{l + n}) \lambda_{l}(v)
   + \mathrm{i} \sum_{l = 1}^n (q_{l} - q_{l + n}) \omega_{l}(v)
 \Big] \mathfrak{s}_{\mathbf{q}}(y)
\end{align*}
and
\begin{align*}
 \mathfrak{L}_{J(v)y} e_{k, \mathbf{q}}^{\pm}(y)
&=  \frac{\mathrm{d}}{\mathrm{d}t} \big|_{t = 0}
  \mathrm{e}^{-J(v)t} e_{k, \mathbf{q}}^{\pm} \big(
 \mathrm{e}^{J(v)t}y
  \big)
 \\
&=\Big[
  \varepsilon \sum_{l = 1}^n
 (q_{l} + q_{l + n} - \delta_{kl}) \lambda_{l}(v)
  + \mathrm{i} \sum_{l = 1}^n
 (q_{l} - q_{l + n}\mp\delta_{kl}) \omega_{l}(v)
 \Big] e_{k, \mathbf{q}}^{\pm}(y).
\end{align*}
After expanding the forms
\begin{gather*}
\mathcal{Y}_{j}(v) y^{j}
  = \sum_{k = 1}^n \sum_{| \mathbf{q} | = j} \big[
 Y_{k, \mathbf{q}}^{+}(v) e_{k, \mathbf{q}}^{+}(y)
 + Y_{k, \mathbf{q}}^{-}(v) e_{k, \mathbf{q}}^{-}(y)
  \big],
 \\
 \mathcal{Z}_{j}(v) y^{j}
  = \sum_{| \mathbf{q} | = j}
 \mathfrak{s}_{\mathbf{q}}(y) Z_{\mathbf{q}}(v),
\end{gather*}
the homological equations are reduced to
\begin{gather*}
\begin{aligned}
&\Big[
  \varepsilon \sum_{l = 1}^n
 (q_{l} + q_{l + n} - \delta_{kl}) \lambda_{l}(v)
  + \mathrm{i} \sum_{l = 1}^n
 (q_{l} - q_{l + n} \mp \delta_{kl}) \omega_{l}(v)
 \Big] Y_{k, \mathbf{q}}^{\pm}(v)
 \\
&=P_{k, \mathbf{q}}^{\pm}(v)
 - R_{k, \mathbf{q}}^{\pm}(v),
\end{aligned} \\
 \Big[
  \varepsilon \sum_{l = 1}^n
 (q_{l} + q_{l + n}) \lambda_{l}(v)
  + \mathrm{i} \sum_{l = 1}^n
 (q_{l} - q_{l + n}) \omega_{l}(v)
 \Big] Z_{\mathbf{q}}(v)
 = Q_{\mathbf{q}}(v) - S_{\mathbf{q}}(v).
\end{gather*}
Taking into account the definitions of $R_{k, \mathbf{q}}^{\pm}(v)$ and 
$S_{\mathbf{q}}(v)$, these equations are soluble for any 
$\varepsilon \in [0, \varepsilon_0]$ with sufficiently small $\varepsilon_0 > 0$, 
and, as a consequence, the same is true for equations~\eqref{eq:homol_equ_Y}.
 Although the coefficients $Y_{k, \mathbf{q}}^{\pm}(v)$ and $Z_{\mathbf{q}}(v)$ 
are complex-valued, one can make the forms $\mathcal{Y}_{j}(v)y^{j}$ and 
$\mathcal{Z}_{j}(v)y^{j}$ take values, respectively, in $\mathbb{R}^{2n}$ and 
$\mathbb{R}^{m}$ by setting $Y_{k, \mathbf{q}}^{\pm}(v)\equiv0$ for any $\mathbf{q}$,
such that $\sum_{l = 1}^n | q_{l} - q_{l + n}\mp\delta_{kl} | = 0$, and
 $Z_{\mathbf{q}}(v)\equiv0$ for any $\mathbf{q}$, such that
 $\sum_{l = 1}^n | q_{l} - q_{l + n} | = 0$. The smoothness properties of 
$\mathcal{Y}_{j}(v)y^{j}, \mathcal{Z}_{j}(v)y^{j}$ are obvious.
\end{proof}

\begin{remark}\label{remark:e_k^pm} \rm
For $\mathbf{q} \in \mathbb{Z}_{+}^{2n}$ with $| \mathbf{q} | = 1$,
the equality $\sum_{l = 1}^n | q_{l} - q_{l + n}\mp\delta_{kl} | = 0$ is 
satisfied if and only if $\mathbf{q} = \mathbf{e}_{k}^{\pm}$, where 
$\mathbf{e}_{k}^{+} \in \mathbb{Z}_{+}^{2n}$ 
($\mathbf{e}_{k}^{-} \in \mathbb{Z}_{+}^{2n}$) is a vector whose
 $k$-th ($(k + n)$-th) coordinate equals $1$ while the other are $0$.
 Hence, in such case,
\begin{equation*}
 R_{k, \mathbf{q}}^{\pm}(v) \neq 0
 \quad
 \text{if and only if}
 \quad
 \mathbf{q} = \mathbf{e}_{k}^{\pm}.
\end{equation*}
\end{remark}

Let
\begin{gather*}
\dot{y} =  J(v)y + \sum_{i \geq 0}\mu^{i}H_i(y, v), \\
\dot{v} =  \mu \Big[ g(v) + \sum_{i \geq 0}\mu^{i}C_i(y, v) \Big] ,
\end{gather*}
be the system obtained from \eqref{eq:formal_eps} by means of the formal 
change of variables~\eqref{eq:transform}. In view of \eqref{eq:F=00003DO(w2p2)} 
and the definition of $g(v)$, we have
\begin{equation*}
 F_0(y, v) = O(y^2),  \quad
 G_0(y, v) = O(y).
\end{equation*}
Therefore, we require that
\begin{equation*}
 X_0(y, v) = O(y^2), 
 \quad
 H_0(y, v) = O(y^2),
 \quad
 C_0(y, v) = O(y),
 \quad
 U_0(y, v) = O(y).
\end{equation*}
Substituting~\eqref{eq:transform} in~\eqref{eq:formal_eps} and equating 
coefficients near like powers of $\mu$, we obtain the following chain of
 homological equations for the unknown coefficients
\begin{gather*}
 \mathfrak{L}_{J(v)y}X _0(y, v)
  = F_0 ( y + X_0(y, v), v )
  - \frac{\partial X_0(y, v)}{\partial y} H_0(y, v)
  - H_0(y, v),
 \\
 \partial_{J(v)y} U_0(y, v)
  = G_0 ( y + X_0(y, v), v )
  - \frac{\partial U_0(y, v)}{\partial y} H_0(y, v)
  - C_0(y, v),
\end{gather*}
and, for $i > 0$,
\begin{align*}
& \mathfrak{L}_{J(v)y} X_i(y, v) \\
&=   \frac{\partial X_{i - k}(y, v)}{\partial y} H_{k}(y, v)
  - \frac{\partial X_{i - 1}(y, v)}{\partial v} g(v)
 \\
&\quad -\sum_{k = 0}^{i - 1}
  \frac{\partial X_{i - k - 1}(y, v)}{\partial v} C_{k}(y, v)
 + \big[
  \frac{\partial F_0(x, v)}{\partial x} X_i(y, v)
 \big]_{x = y + X_0(y, v)}
 \\
&\quad  +    \frac{1}{i!} \frac{\partial^{i}}{\partial\mu^{i}}
  \big|_{\mu = 0} \Big[
   J \Big(
 v + \mu\sum_{j = 0}^{i - 1} \mu^{j}U_{j}(y, v)
   \Big) \Big(
 y + \sum_{k = 0}^{i - 1} \mu^{k} X_{k}(y, v)
   \Big)
  \Big]
  \\
&\quad +
  \frac{1}{i!} \frac{\partial^{i}}{\partial\mu^{i}}
 \big|_{\mu = 0} \sum_{l = 0}^{i} \mu^{l}F_{l} \Big(
  y + \sum_{k = 0}^{i - 1} \mu^{k} X_{k}(y, v),
  v + \mu \sum_{j = 0}^{i - 1} \mu^{j} U_{j}(y, v)
 \Big) - H_i(y, v),
\end{align*}
\begin{align*}
& \partial_{J(v)y} U_i(y, v)\\
&=   - \sum_{k = 0}^{i} 
 \frac{\partial U_{i - k}(y, v)}{\partial y} H_{k}(y, v)
  - \frac{\partial U_{i - 1}(y, v)}{\partial v} g(v)
 \\
&\quad -\sum_{k = 0}^{i - 1}
  \frac{\partial U_{i - k - 1}(y, v)}{\partial v} C_{k}(y, v)
 + \big[
  \frac{\partial G_0(v, x)}{\partial x} X_i(y, v)
 \big]_{x = y + X_0(y, v)}
 \\
&\quad +  \frac{1}{i!} \frac{\partial^{i}}{\partial\mu^{i}}
 \big|_{\mu = 0} g \Big(
  v + \mu \sum_{j = 0}^{i - 1} \mu^{j} U_{j}(y, v)
 \Big)
 \\
&\quad +  \frac{1}{i!} \frac{\partial^{i}}{\partial\mu^{i}}
 \big|_{\mu = 0}
 \sum_{l = 0}^{i}
  \mu^{l}G_{l} \Big(
 y + \sum_{k = 0}^{i - 1} \mu^{k}X_{k}(y, v),
 v + \mu \sum_{j = 0}^{i - 1} \mu^{j} U_{j}(y, v)
  \Big)
 - C_i(y, v).
\end{align*}
On account of~\eqref{eq:F=00003DO(w2p2)} and the definition of $g(v)$, 
we have the following expansions.
\begin{gather*}
 F_0(y, v) = \sum_{j \geq 2} F_{0j}(v)y^{j},
 \quad
 F_i(y, v) = \sum_{j \geq 0}F_{ij}(v)y^{j},
 \quad i \geq 1,
 \\
 G_0(y, v) = \sum_{j \geq 1} G_{0j}(v)y^{j},
 \quad
 G_i(y, v) = \sum_{j \geq 0} G_{ij}(v)y^{j},
 \quad
 i \geq 1.
\end{gather*}
Thus, the unknown functions can be sought in the form
\begin{gather*}
 X_0(y, v) = \sum_{j \geq 2} X_{0j}(v)y^{j},
 \quad
 X_i(y, v) = \sum_{j \geq 0} X_{ij}(v)y^{j},
 \quad
 i \geq 1,
 \\
 H_0(y, v) = \sum_{j \geq 2} H_{0j}(v)y^{j},
 \quad
 H_i(y, v) = \sum_{j \geq 0} H_{ij}(v)y^{j},
 \quad i \geq 1,
\\
 U_0(y, v) = \sum_{j \geq 1} U_{0j}(v)y^{j},
 \quad
 U_i(y, v) = \sum_{j \geq 0} U_{ij}(v)y^{j},
 \quad
 i \geq 1,
 \\
 C_0(y, v) = \sum_{j \geq 1} C_{0j}(v)y^{j},
 \quad
 C_i(y, v) = \sum_{j \geq 0} C_{ij}(v)y^{j},
 \quad
 i \geq 1,
\end{gather*}
It is easily seen that the coefficients of these expansions satisfy 
the chain of equations
\begin{equation*}
 \mathfrak{L}_{J(v)y} [ X_{ij}(v)y^{j} ]
  = P_{ij}(v)y^{j} - H_{ij}(v)y^{j},
 \quad
 \partial_{J(v)y} [ U_{ij}(v)y^{j} ]
  = Q_{ij}(v)y^{j} - C_{ij}(v)y^{j},
\end{equation*}
where $P_{ij}(v)$ and $Q_{ij}(v)$ can be determined subsequently. 
In fact, $P_{02}(v) := F_{02}(v)$, $Q_{01}(v) := G_{01}(v)$, and now,
 if $v \in \mathcal{D}_{N}$, one can use Proposition~\ref{proposition:homol_equ}
to get $H_{02} = 0$, $X_{02}$, $C_{01} = 0$, $U_{01}$, and then 
subsequently find $P_{0j}$, $H_{0j}$, $X_{0j}$ for $j = 3, \ldots, N$, 
and $Q_{0j}$, $C_{0j}$, $U_{0j}$ for $j = 2, \ldots, N$. 
If $0 \leq k < i$, $l \leq N$ and the coefficients $X_{kl}(v)$, 
$H_{kl}(v)$, $U_{kl}(v)$, $C_{kl}(v)$ are already known, then one can determine 
$P_{i0}(v)$ and subsequently find $H_{ij}(v)$, $X_{ij}(v)$, 
$P_{i, j + 1}(v)$ for $j = 0, \ldots, N$. (Note that 
\eqref{eq:R0 = 00003D0} yields $H_{i0}(v) = 0$.) After that, $Q_{i0}(v)$ 
can be determined, and subsequently $C_{ij}(v)$, $U_{ij}(v)$, $Q_{i, j + 1}(v)$ 
may be found for $j = 0, \ldots, N$. At last, for any $j > N$, we assign 
$H_{ij}(v) \equiv P_{ij}(v)$, $X_{ij}(v) \equiv 0$, $C_{ij}(v) \equiv Q_{ij}(v)$ 
and $U_{ij}(v) \equiv 0$.

This result can be summed up as the following proposition.

\begin{proposition}\label{proposition:exis_N_norm_form}
Suppose that $P \geq 2$, $N \geq 3$ and $\mathcal{D}_{N}$ satisfy the 
conditions of Proposition~\ref{proposition:homol_equ}.
Then there exist $\delta_0 > 0$ and $\mu_0 > 0$, such that for any 
$\varepsilon \in [0, \varepsilon_0]$ the smooth diffeomorphic 
change of variables
\begin{gather*}
 x = y + \sum_{j = 2}^{N} X_{0j}(v)y^{j}
  + \sum_{i = 1}^{P} \mu^{i} \sum_{j = 0}^{N} X_{ij}(v) y^{j},
 \\
 u = v + \sum_{j = 1}^{N} U_{0j}(v) y^{j}
  + \sum_{i = 1}^{P} \mu^{i} \sum_{j = 0}^{N} U_{ij}(v) y^{j}
\end{gather*}
defined on the set 
$\{ (y, v, \mu) \colon \| y \| < \delta_0, v \in \mathcal{D}_{N}, 
\mu \in [0, \mu_0]\}$ transforms system~\eqref{eq:main1} into
\begin{equation}\label{eq:norm_form_N}
\begin{gathered}
 \dot{y} = J(v)y
  + \sum_{j = 3}^{N} H_{0j}(v) y^{j}
  + \sum_{i = 1}^{P} \mu^{i} \sum_{j = 1}^{N} H_{ij}(v) y^{j}
  + O (
 \| y \| ^{N + 1} + \mu^{P + 1}
  ),
 \\
 \dot{v} = \mu \Big[
  g(v) + \sum_{j = 2}^{N} C_{0j}(v) y^{j}
  + \sum_{i = 1}^{P} \mu^{i} \sum_{j = 0}^{N} C_{ij}(v) y^{j}
  + O (
 \| y \| ^{N + 1} + \mu^{P + 1}
  )
 \Big].
\end{gathered}
\end{equation}
Here the homogeneous forms in the right-hand sides admit the expansions
\begin{gather*}
 H_{ij}(v) y^{j} = \sum_{k = 1}^n \Big[
  {\sum_{| \mathbf{q} | = j}}^{+}
 h_{i, k, \mathbf{q}}^{+}(v) e_{k, \mathbf{q}}^{+}(y)
  + {\sum_{ | \mathbf{q} | = j}}^{-}
 h_{i, k, \mathbf{q}}^{-}(v)e_{k, \mathbf{q}}^{-}(y)
 \Big],
 \\
 C_{ij}(v)y^{j} = \sum_{k = 1}^n
  {\sum_{| \mathbf{q} | = j}}'
 \mathfrak{s}_{\mathbf{q}}(y) c_{i, \mathbf{q}}(v),
\end{gather*}
where $h_{i, k, \mathbf{q}}^{\pm}(\cdot) \in \mathrm{C}^{\infty} 
( \mathcal{D}_{N} ; \mathbb{C} )$, $c_{i, \mathbf{q}}(\cdot) 
\in \mathrm{C}^{\infty} ( \mathcal{D}_{N} ; \mathbb{C}^{m} )$ 
and the summations ${\sum_{ | \mathbf{q} | = j}^{+}}$, 
${\sum_{ | \mathbf{q} | = j}^{-}}$, ${\sum_{ | \mathbf{q} | = j}'}$ 
are performed on all vectors $\mathbf{q} \in \mathbb{Z}_{+}^{2n}$ with 
$| \mathbf{q} | = j$ whose components satisfy, respectively, the equalities 
$q_{l} = q_{l + n} + 1$, $q_{l + n} = q_{l} + 1$ and $q_{l} = q_{l + n}$ 
for all $l \in \{ 1, \ldots, n\}$.
\end{proposition}

\begin{remark}\label{remark:proposition:exis_N_norm_form} \rm
In view of Remark~\ref{remark:proposition_homol_equ_1}, under hypothesis (H2)
the assertion of Proposition~\ref{proposition:exis_N_norm_form} is true for
$N = 3$ and $\mathcal{D}_{3} = \mathcal{V}$. If hypothesis (H3) is valid, 
then this proposition is correct in a sufficiently small neighborhood
of any stationary point $v_{\ast} \in \mathcal{W}$.
\end{remark}

On account of Remark \ref{remark:e_k^pm},
\begin{gather*}
 \mathfrak{s}_{\mathbf{q}} (
  S ( z, \bar{z} )
 ) = ( z, \bar{z} )^{\mathbf{q}},
 \quad
 S^{-1} e_{k, \mathbf{q}}^{+} (Sz)
  = ( z, \bar{z} )^{\mathbf{q}} \mathbf{e}_{k}^{+},
 \quad
 S^{-1} e_{k, \mathbf{q}}^{+} (Sz)
  = ( z, \bar{z} )^{\mathbf{q}} \mathbf{e}_{k}^{-},
 \\
 S^{-1} H_{i1}(v) (
  S ( z, \bar{z} )
 )^{1} = \sum_{k = 1}^n \left[
  z_{k} h_{i, k, \mathbf{e}_{k}^{+}}^{+}(v) \mathbf{e}_{k}^{+}
  + \bar{z}_{k} h_{i, k, \mathbf{e}_{k}^{-}}^{-}(v) \mathbf{e}_{k}^{-}
 \right],
\end{gather*}
which means that the matrices of each of the linear form
$S^{-1}H_{i1}(v) ( S ( z, \bar{z} ) ) ^{1}$ are diagonal and the change 
of variables
\begin{equation*}
 y = S ( z, \bar{z} )
  = \sum_{k = 1}^n (
 z_{k}s_{k} + \bar{z}_{k}\bar{s}_{k}
  ),
 \quad
 \| z \|
  < \| S \| \delta_0 =: \delta_1,
\end{equation*}
reduces system~\eqref{eq:norm_form_N} to the form~\eqref{eq:normform_2}, where
\begin{equation*}
 \eta_{j, k}(v) := h_{i, k, \mathbf{e}_{k}^{+}}^{+}(v),
 \quad
 h_{j, k, \mathbf{p}}(v)
  := h_{j, k, (\mathbf{p, \mathbf{p}})
  + \mathbf{e}_{k}^{+}}^{+}(v),
 \quad
 g_{j, \mathbf{p}}(v)
  := c_{j, (\mathbf{p, \mathbf{p}})}(v),
\end{equation*}
and $(\mathbf{p, \mathbf{p}}) := (p_1, \ldots, p_{n}, p_1, \ldots, p_{n})$.


\section{Summary}

In this article we have examined a system of oscillators with weak and 
slow coupling which demonstrates a dynamic bifurcation of multi-frequency 
oscillations. Having adopted results of the static bifurcation theory, 
we have shown that when the system's parameters slowly evolve and the 
static parameters are sufficiently small, then certain general conditions 
guarantee occurrence of the following transient process for typical forward 
trajectories within a small neighborhood of the slow surface. 
While the slow component $u(t)$ is far from the stationary points of the 
Morse function $V$ and lies inside the stability zone of the fast subsystem, 
the fast component $(w(t),\dot{w}(t))$ exhibits damping oscillations with the
 amplitude approaching zero. However, as soon as $u(t)$ leaves the stability
 zone and later enters the zone of instability of the fast subsystem, 
the amplitude starts to grow and eventually the forward trajectory is 
attracted by an invariant torus, which means establishment of some multi-frequency 
oscillatory regime.

It was shown that almost all forward trajectories, in terms of the Lebesgue measure, 
starting from the neighborhood of the slow surface demonstrate such a behavior. 
More than that, they are attracted by trajectories on the stable $n$-dimensional 
invariant tori, whereas all other forward trajectories of the system lie 
on the stable manifolds of hyperbolic tori of dimensions less than $n$. 
This enables us to easily categorize the trajectories by the type of their 
ultimate behavior.

At last, we have also considered a practical example which depicts occurrence 
of the multi-frequency bifurcation in a circuit of two coupled oscillators.


\begin{thebibliography}{10}

\bibitem{AAISh94}
V.~S.~Afrajmovich, V.~I.~Arnold, Yu.~S.~Il'yashenko, L.~P.~Shil'nikov; 
\emph{Dynamical Systems V. Encyclopaedia of Mathematical Sciences.} 
Springer-Verlag, 1994.

\bibitem{AnoDV05}
D.~V.~Anosov; 
\emph{On the development of the theory of dynamical systems during the past 
quarter century.} Surveys in modern mathematics, London Math. Soc. 
Lecture Note Ser., 321, Cambridge Univ. Press, Cambridge, 2005, 70--185.

\bibitem{Ano05}
O.~D.~Anosova; \emph{Invariant manifolds and dynamic bifurcations.} 
Russ. Math. Surv. 60 (1) (2005) 151--153.

\bibitem{ACN2016} P.~Ashwin, S.~Coombes, R.~Nicks; 
\emph{Mathematical frameworks for oscillatory network dynamics in neuroscience.} 
J. Math. Neurosci. 6:2 (2016).

\bibitem{Ben91}
\'{E}~Beno\^{\i}t (ed.); 
\emph{Dynamic bifurcations.} Lecture Notes in Mathematics. 1493. Springer-Verlag, 
Berlin, 1991.

\bibitem{Bib81}
Yu.~N.~Bibikov; 
\emph{Hopf bifurcation for quasiperiodic motion.} Differ. Equ., 16 (1981) 977--982.

\bibitem{Bib90}
Yu.~N.~Bibikov; 
\emph{Bifurcation of a stable invariant torus from an equilibrium.} 
Math. Notes 48 (1) (1990), 632--635.

\bibitem{Bib91}
Yu.~N.~Bibikov; 
\emph{Multi-frequency non-linear oscillations and their bifurcations (in Russian).}
 Leningrad Gos. Univ., Leningrad, 1991.

\bibitem{KryBog34}
N.~N.~Bogolyubov, N.~M.~Krylov; 
\emph{Application of the methods of non-linear mechanics to the theory of stationary
 oscillations (in Russian).} Kyiv: Publ. Acad. of Sci. Ukr. SSR, 1934.

\bibitem{BMS69}
N.~N.~Bogolyubov, Yu.~A.~Mitropol'skij, A.~M.~Samojlenko; 
\emph{Methods of accelerated convergence in nonlinear mechanics.} 
Delhi, India: Hindustan Publishing Corporation; Berlin-Heidelberg-New York: 
Springer-Verlag, 1976.

\bibitem{BuNeSch04}
V.~F.~Butuzov, N.~N.~Nefedov, K.~R.~Schneider; 
\emph{Singularly perturbed problems in case of exchange of stabilities.} 
J. Math. Sci., 121 (1) (2004), 1973--2079.

\bibitem{Car84}
P.~Cartier; 
\emph{Singular perturbations of ordinary differential equations and 
non-standard analysis (in Russian).} Uspekhi Mat. Nauk., 39 (2) (1984) 57--76.

\bibitem{EKST2013}
Yu.~P.~Emelianova, A.~P.~Kuznetsov, I.~R.~Sataev, L.~V.~Turukina; 
\emph{Synchronization and multi-frequency oscillations in the low-dimensional 
chain of the self-oscillators.}
Physica D: Nonlinear Phenomena, 244 (1) (2013), 36--49.

\bibitem{KKE2012}
T.~Endo, K.~Kamiyama, M.~Komuro; 
\emph{Bifurcation of quasi-periodic oscillations in mutually coupled hard-type 
oscillators: demonstration of unstable quasi-periodic orbits.} 
Internat. J. Bifur. Chaos Appl. Sci. Engrg., 22 (6) (2012).

\bibitem{Fen71}
N.~Fenichel; 
\emph{Persistence and smoothness of invariant manifolds and flows.}
 Indiana Univ. Math., 21 (3) (1971), 193--226.

\bibitem{Gantmacher}
F.~R.~Gantmacher; 
\emph{The Theory of Matrices.} American Mathematical Society, 1998.

\bibitem{Gilsinn}
D.~E.~Gilsinn; 
\emph{Constructing invariant tori for two weakly coupled van der Pol oscillators}
 Nonlinear Dynamics, 4 (3) (1993), 289--308.

\bibitem{GKR08}
S.~D.~Glyzin, A.~Yu.~Kolesov, N.~Kh.~Rozov; 
\emph{Blue sky catastrophe in relaxation systems with one fast and two slow 
variables.} Differ. Equ., 44 (2) (2008), 161--175.

\bibitem{Hale61}
J.~K.~Hale; 
\emph{Integral manifolds of perturbed differential systems.} Ann. 
of Math. Second Series, 73 (3) (1961), 496--531.

\bibitem{Han91}
H.~Hanssmann; 
\emph{A Survey on Bifurcations of Invariant Tori.} New Advances in Celestial 
Mechanics and Hamiltonian Systems. HAMSYS-2001 (2004), 109--121.

\bibitem{Har70}
P.~Hartman; 
\emph{Ordinary differential equations.} John Wiley and Sons, 
New York-London-Sydney, 1964.

\bibitem{HopIzh97}
F.~C.~Hoppensteadt, E.~M.~Izhikevich; 
\emph{Weakly connected neural networks.}
Applied Mathematical Sciences, 126. New York, NY: Springer. xvi, 1997.

\bibitem{IzhiHoppen}
F.~C.~Hoppensteadt, E.~M.~Izhikevich; 
\emph{Slowly coupled oscillators: phase dynamics and synchronization.} 
SIAM J. Appl. Math., 63 (6) (2003), 1935--1953.

\bibitem{IooJos90}
G.~Iooss, D.~D.~Joseph; 
\emph{Elementary Stability and Bifurcation Theory.} Springer-Verlag New York, 1990.

\bibitem{Kir78}
P.~S.~Kireev; 
\emph{Semiconductor physics, 2nd ed.} Mir Publishers, Moscow, 1978.

\bibitem{KMR08}
A.~Yu.~Kolesov, E.~F.~Mishchenko, N.~Kh.~Rozov; 
\emph{On one bifurcation in relaxation systems.} 
Ukranian Math. J., 60 (1) (2008), 66--77.

\bibitem{Kuz98}
Yu.~A.~Kuznetsov; 
\emph{Elements of applied bifurcation theory.} Springer, New York, 1998.

\bibitem{Lie11}
S.~Liebscher; 
\emph{Dynamics near manifolds of equilibria of codimension one and bifurcation 
without parameters.} Electron. J. Differential Equations, 2011 (63) (2011), 1--12.

\bibitem{SchMar84}
J.~Marsden, J.~Scheurle; 
\emph{Bifurcation to quasi-periodic tori in the interaction of steady 
state and hopf bifurcations.} SIAM J. Math. Anal., 15 (6) (1984), 1055--1074.

\bibitem{NayMoo79}
D.~T.~Mook, A.~H.~Nayfeh; 
\emph{Nonlinear oscillations.} John Wiley \& Sons, New York, 1979.

\bibitem{SamMit76}
Yu.~A.~Mitropol'skij, A.~M.~Samojlenko; 
\emph{Multifrequency oscillations of weakly nonlinear second-order systems.}
 Ukrainian Math. J., 28 (1977), 572--586.

\bibitem{Nei59}
Yu.~I.~Neimark; 
\emph{On some cases of dependency of periodic motions on parameters (in Russian).}
 Dokl. AN SSSR., 129 (4) (1959), 736--739.

\bibitem{Nei85}
A.~I.~Nejshtadt;
 \emph{Asymptotic investigation of the loss of stability by an equilibrium as
 a pair of eigenvalues slowly cross the imaginary axis.}
 Uspkhi Mat. Nauk, 40 (5) (1985), 190--191

\bibitem{Nei87}
A.~I.~Nejshtadt; 
\emph{Persistence of stability loss for dynamical bifurcations. I.} Differ. 
Equ., 23 (12) (1987), 1385--1391.

\bibitem{Nei88}
A.~I.~Nejshtadt; 
\emph{On deferring the loss of stability for dynamic bifurcations. II.} 
Differ. Equ., 24 (2) (1988), 171--176.

\bibitem{Nikaido}
H.~Nikaido; \emph{Convex Structures and Economic Theory.} Academic Press, 
New York and London, 1968.

\bibitem{SamParRep}
I.~O.~Parasyuk, B.~V.~Repeta, A.~M.~Samoilenko; 
\emph{Dynamical bifurcation of multifrequency oscillations in a fast-slow system.}
 Ukrainian Math. J., 67 (7) (2015), 1008--1037.

\bibitem{ParRep16}
I.~O.~Parasyuk, B.~V.~Repeta; 
\emph{Hyperbolic invariant tori of a fast-slow system exhibiting the dynamical
 bifurcation of multi-frequency oscillations (in Ukrainian).}
 Nelinijni Kolyvannya, 19 (1) (2016), 101--121.

\bibitem{SamPet04}
R.~Petryshyn, A.~Samoilenko; 
\emph{Multifrequency oscillations of nonlinear systems.
Mathematics and its Applications.} Kluwer Academic Publishers, Dordrecht, 2004.

\bibitem{RachSch05}
D.~Rachinskii, K.~Schneider; 
\emph{Dynamic Hopf bifurcations generated by nonlinear terms.} 
J. Differential Equations, 210 (1) (2005), 65--86.

\bibitem{Sack64}
R.~Sacker; 
\emph{On invariant surfaces and bifurcation of periodic solutions of ordinary 
differential equations.} Report IMM-NYU 333, New York University, 1964.

\bibitem{Sack65}
R.~Sacker; 
\emph{A new approach to perturbation theory of invariant surfaces.} 
Comm. Pure Appl. Math., 18 (1965) 717--732.

\bibitem{Sam91}
A.~M.~Samoilenko; 
\emph{Elements of mathematical theory of multi-frequency oscillations.} 
Kluwer Academic Publishers, Dordrecht, 1991.

\bibitem{Sam97}
A.~M.~Samoilenko; 
\emph{Perturbation theory of smooth invariant tori of dynamical systems.}
 Nonlinear Anal., 30 (5) (1997), 3121--3133.

\bibitem{Shch04}
S.~A.~Shchegolev; 
\emph{On some classes of solutions of two-dimensional multifrequency quasilinear 
differential systems.} Differ. Equ., 40 (10) (2004), 1413--1418.

\bibitem{ShShT05}
A.~Shilnikov, L.~Shilnikov, D.~Turaev; 
\emph{Blue-sky catastrophe in singularly perturbed systems.} 
Mosc. Math. J., 5 (1) (2005), 269--282.

\bibitem{Shi73}
M.~A.~Shishkova; 
\emph{Examination of a system of differential equations with a small parameter 
in the highest derivatives.} Sov. Math., Dokl., 14 (1973), 483--487.

\end{thebibliography}

\end{document}