# SLOW INVARIANT MANIFOLDS OF
SLOW-FAST DYNAMICAL SYSTEMS

Jean-Marc GINOUX  Aix Marseille Univ, Université de Toulon, CNRS, CPT,
Marseille, France<EMAIL_ADDRESS>
###### Abstract

Slow-fast dynamical systems, i.e., singularly or non-singularly perturbed
dynamical systems possess slow invariant manifolds on which trajectories
evolve slowly. Since the last century various methods have been developed for
approximating their equations. This paper aims, on the one hand, to propose a
classification of the most important of them into two great categories:
singular perturbation-based methods and curvature-based methods, and on the
other hand, to prove the equivalence between any methods belonging to the same
category and between the two categories. Then, a deep analysis and comparison
between each of these methods enable to state the efficiency of the Flow
Curvature Method which is exemplified with paradigmatic Van der Pol singularly
perturbed dynamical system and Lorenz slow-fast dynamical system.

##  I Introduction

A great number of phenomena can generally be modeled with dynamical systems,
i.e., with sets of nonlinear ordinary differential equations. Among their huge
variety, the evolution of some of them are characterized by the existence of
at least two time scales: a slow time and a fast time. The metaphorical
example of the Tantalus cup is often used to emphasize this kind of evolution
since it fills slowly and quickly empties Le Corbeiller (1931). Such slow-fast
evolution is transcribed by the presence of at least a small multiplicative
parameter $\varepsilon$ in the velocity vector field of these dynamical
systems which have been called since singularly perturbed dynamical systems.

The classical geometric theory of differential equations was originally
developed by Poincaré in his famous memoirs entitled “On the curves defined by
differential equations” Poincaré (1881, 1882, 1885, 1886) in which he
formalized the concept of dynamical systems consisting of set of ordinary
differential equations (ODE) as well as that of phase space (For a History of
Nonlinear Oscillations Theory, see Ginoux Ginoux (2017)). Few years later,
perturbation method, which had been introduced by Siméon Denis Poisson in
1830, were greatly improved by Poincaré in his New Methods of Celestial
Mechanics Poincaré (1892). According to Minorsky (Minorsky, 1962, p. 213):

> “Poincaré opens the new period in the perturbation theory in that the
> periodicity begins to play the primary role, and this question became of
> fundamental importance in the theory of oscillations.”

Indeed, this idea was emphasized in the very beginning of the chapter III of
the first volume of his “New Methods” by Poincaré (Poincaré, 1892, p. 82)
through this famous sentence:

> “In addition, these periodic solutions are so valuable for us because they
> are, so to say, the only breach by which we may attempt to enter an area
> heretofore deemed inaccessible.”

In the 1930-1950s Andronov & Chaikin Andronov & Chaikin (1937), Levinson
Levinson (1949) and Tikhonov Tikhonov (1952) generalized Poincaré’s ideas and
stated that singularly perturbed systems possess invariant manifolds on which
trajectories evolve slowly, and toward which nearby orbits contract
exponentially in time (either forward or backward) in the normal directions.
These manifolds have been called asymptotically stable (or unstable) slow
invariant manifolds. Then, Fenichel theory Fenichel (1971, 1974, 1977, 1979)
for the persistence of normally hyperbolic invariant manifolds enabled to
establish the local invariance of slow invariant manifolds that possess both
expanding and contracting directions and which were labeled slow invariant
manifolds. The theory of invariant manifolds for an ordinary differential
equation is based on the work of Hirsch, et al. Hirsch et al. (1977).

During the last century, various methods have been developed to compute the
slow invariant manifold or, at least an asymptotic expansion in power of
$\varepsilon$.

The seminal works of Wasow Wasow (1965), Cole Cole (1968), O’Malley O’Malley
(1974, 1991) and Fenichel Fenichel (1971, 1974, 1977, 1979) to name but a few,
gave rise to the so-called Geometric Singular Perturbation Theory. According
to this theory, existence as well as local invariance of the slow invariant
manifold of singularly perturbed dynamical systems has been stated. Then, the
determination of the slow invariant manifold equation turned into a regular
perturbation problem in which one generally expected the asymptotic validity
of such expansion to breakdown O’Malley (1991). In the framework of the
Geometric Singular Perturbation Theory, the zero-order approximation in
$\varepsilon$ of the slow invariant manifold is called singular approximation
Rossetto (1987) or critical manifold Guckenheimer et al. (2000). At the end of
the 1980s, Rossetto Rossetto (1986, 1987) developed the Successive
Approximations Method to approximate the slow invariant manifold of singularly
perturbed dynamical systems. Then, Gear et al. Gear et al. (2005) and Zagaris
et al. Zagaris et al. (2009) used the Zero-Derivative Principle for the same
purpose.

Beside, singularly perturbed dynamical systems there are dynamical systems
without an explicit timescale splitting. In this case the Geometric Singular
Perturbation Theory can not be applied anymore. Some of these systems, which
have been called slow-fast dynamical systems, have the following property:
their Jacobian matrix has at least a fast eigenvalue, i.e. with the largest
absolute value of the real part. Let’s notice that a singularly perturbed
dynamical system is slow-fast while a slow-fast dynamical system is not
necessary singularly perturbed (For a proof of this statement see Rossetto et
al. Rossetto et al. (1998)). Thus, in the beginning of the 1990s various
approaches have been proposed in order to approximate slow manifolds of such
slow-fast dynamical systems. As recalled by Lebiedz and Unger Lebiedz & Unger
(2016) “In 1992, Maas and Pope introduced the intrinsic low-dimensional
manifold (ILDM) method, which has become very popular and widely used in the
reactive flow community, in particular in combustion applications.” Two years
later, Brøns and Bar-Eli Brøns & Bar-Eli (1994) presented an approach, only
applicable to two-dimensional systems, based on the inflection line, i.e. the
locus of points where the curvature of trajectories vanishes, to provide an
approximation of the slow manifold. Then, Rossetto et al. Rossetto et al.
(1998) proposed to use the Tangent Linear System Approximation for the same
purpose. Since 2005, a new approach of $n$-dimensional singularly perturbed
dynamical systems or slow-fast dynamical systems based on the location of the
points where the curvature of trajectory curves vanishes, called Flow
Curvature Method has been developed by Ginoux et al. Ginoux & Rossetto
(2006a); Ginoux et al. (2008) and Ginoux Ginoux (2009). So, the aim of this
paper is to propose a classification of the most important method of
approximation of slow manifolds into two great categories:

  * •

singular perturbation-based methods and

  * •

curvature-based methods.

Thus, it is stated that Geometric Singular Perturbation Method encompasses
Successive Approximations Method and Zero-Derivative Principle which are
perfectly identical and belong to the first category. Then, we show on the one
hand that both Intrinsic Low-Dimensional Manifold Method (IDLM) and Tangent
Linear System Approximation (TLSA) are also identical. On the other hand we
prove that the Flow Curvature Method is a generalization of the Inflection
Line Method (ILM) which cannot in anyway be extended to higher dimensions than
two. The reasons of such impossibility which mainly consists in the arc length
parametrization are fully detailed. Finally, we recall the identity between
Flow Curvature Method and Tangent Linear System Approximation established in
our previous works more than ten years ago. It follows that Flow Curvature
Method encompasses these three other methods (IDLM, TLSA, ILM) which belong to
the second category. After having divided these methods of approximation of
slow manifolds into two categories, we also recall the identity between Flow
Curvature Method and Geometric Singular Perturbation Method for
$n$-dimensional singularly perturbed dynamical systems up to suitable order in
$\varepsilon$. At last, a deep analysis and comparison between each of these
various methods enable to state the efficiency of the Flow Curvature Method
what is exemplified with paradigmatic singularly perturbed dynamical systems
and slow-fast dynamical systems of dimension two and three.

##  II Singular Perturbation-Based Methods

###  II.1 Singularly perturbed dynamical systems

Following the works of Jones Jones (1994) and Kaper Kaper (1999) some
fundamental concepts and definitions for systems of ordinary differential
equations with two time scales, i.e., for singularly perturbed dynamical
systems are briefly recalled. Thus, we consider systems of differential
equations of the form:

$\left\\{\begin{aligned}
{\vec{x}}^{\prime}&=\vec{f}\left({\vec{x},\vec{y},\varepsilon}\right),\hfill\\\
{\vec{y}}^{\prime}&=\varepsilon\vec{g}\left({\vec{x},\vec{y},\varepsilon}\right),\hfill\\\
\end{aligned}\right.$ (1)

where $\vec{x}\in\mathbb{R}^{m}$, $\vec{y}\in\mathbb{R}^{p}$,
$\varepsilon\in\mathbb{R}^{+}$, and the prime denotes differentiation with
respect to the independent variable $t$. The functions $\vec{f}$ and $\vec{g}$
are assumed to be $C^{\infty}$ functions (In certain applications these
functions will be supposed to be $C^{r}$, $r\geqslant 1$) of $\vec{x}$,
$\vec{y}$ and $\varepsilon$ in $U\times I$, where $U$ is an open subset of
$\mathbb{R}^{m}\times\mathbb{R}^{p}$ and $I$ is an open interval containing
$\varepsilon=0$.

In the case when $\varepsilon\ll 1$, i.e., is a small positive number, the
variable $\vec{x}$ is called fast variable, and $\vec{y}$ is called slow
variable. Using Landau’s notation: $O(\varepsilon^{l})$ represents a real
polynomial in $\varepsilon$ of $l$ degree, with $l\in\mathbb{Z}$, it is used
to consider that generally $\vec{x}$ evolves at an $O\left(1\right)$ rate;
while $\vec{y}$ evolves at an $O\left(\varepsilon\right)$ slow rate.

Reformulating the system (1) in terms of the rescaled variable
$\tau=\varepsilon t$, we obtain:

$\left\\{\begin{aligned}
\varepsilon\dot{\vec{x}}&=\vec{f}\left({\vec{x},\vec{y},\varepsilon}\right),\hfill\\\
\dot{\vec{y}}&=\vec{g}\left({\vec{x},\vec{y},\varepsilon}\right).\hfill\\\
\end{aligned}\right.$ (2)

The dot $\left(\cdot\right)$ represents the derivative with respect to the new
independent variable $\tau$. The independent variables $t$ and $\tau$ are
referred to the fast and slow times, respectively, and (1) and (2) are called
the fast and slow systems, respectively. These systems are equivalent whenever
$\varepsilon\neq 0$, and they are labeled singular perturbation problems when
$\varepsilon\ll 1$, i.e., is a small positive parameter. The label “singular”
stems in part from the discontinuous limiting behavior in the system (1) as
$\varepsilon\to 0$. In such case, the system (1) reduces to an $m$-dimensional
system called reduced fast system, with the variable $\vec{y}$ as a constant
parameter:

$\left\\{\begin{aligned}
{\vec{x}}^{\prime}&=\vec{f}\left({\vec{x},\vec{y},0}\right),\hfill\\\
{\vec{y}}^{\prime}&=\vec{0}.\hfill\\\ \end{aligned}\right.$ (3)

System (2) leads to the following differential-algebraic system called reduced
slow system which dimension decreases from $m+p$ to $p$:

$\left\\{\begin{aligned}
\vec{0}&=\vec{f}\left(\vec{x},\vec{y},0\right),\hfill\\\
\dot{\vec{y}}&=\vec{g}\left(\vec{x},\vec{y},0\right).\hfill\\\
\end{aligned}\right.$ (4)

By exploiting the decomposition into fast and slow reduced systems (3) and
(4), the geometric approach reduced the full singularly perturbed system to
separate lower-dimensional regular perturbation problems in the fast and slow
regimes, respectively.

###  II.2 Fenichel geometric theory

Fenichel geometric theory for general systems (1), i.e., a theorem providing
conditions under which normally hyperbolic invariant manifolds in system (1)
persist when the perturbation is turned on, i.e., when $0<\varepsilon\ll 1$ is
briefly recalled in this subsection. This theorem concerns only compact
manifolds with boundary.

####  II.2.1 Normally hyperbolic manifolds

Let’s make the following assumptions about system (1):

(H1) The functions $\vec{f}$ and $\vec{g}$ are $C^{\infty}$ in $U\times I$,
where U is an open subset of $\mathbb{R}^{m}\times\mathbb{R}^{p}$ and I is an
open interval containing $\varepsilon=0.$

(H2) There exists a set $M_{0}$ that is contained in
$\left\\{{\left({\vec{x},\vec{y}}\right):\vec{f}\left({\vec{x},\vec{y},0}\right)={\vec{0}}}\right\\}$
such that $M_{0}$ is a compact manifold with boundary and $M_{0}$ is given by
the graph of a $C^{\infty}$ function $\vec{y}=\vec{Y}_{0}\left(\vec{x}\right)$
for $\vec{x}\in D$, where $D\subseteq\mathbb{R}^{p}$ is a compact, simply
connected domain and the boundary of D is an $\left({p-1}\right)$ dimensional
$C^{\infty}$ submanifold. Finally, the set D is overflowing invariant with
respect to (2) when $\varepsilon=0.$

(H3) $M_{0}$ is normally hyperbolic relative to (3) and in particular it is
required for all points $\vec{p}\in M_{0}$, that there are k (resp. l)
eigenvalues of $D_{\vec{y}}\vec{f}\left({\vec{p},0}\right)$ with positive
(resp. negative) real parts bounded away from zero, where $k+l=m.$

####  II.2.2 Fenichel persistence theory for singularly perturbed systems

For compact manifolds with boundary, Fenichel’s persistence theory states
that, provided the hypotheses $(H_{1})-(H_{3})$ are satisfied, the system (1)
has a slow (or center) manifold, and this slow manifold has fast stable and
unstable manifolds.

Theorem for compact manifolds with boundary:

Let system (1) satisfy the conditions (H${}_{1})-(H_{3}$). If $\varepsilon>0$
is sufficiently small, then there exists a function
$\vec{Y}\left({\vec{x},\varepsilon}\right)$ defined on D such that the
manifold
$M_{\varepsilon}=\\{{\left({\vec{x},\vec{y}}\right):\vec{y}=\vec{Y}\left({\vec{x},\varepsilon}\right)}\\}$
is locally invariant under (1). Moreover,
$\vec{Y}\left({\vec{x},\varepsilon}\right)$ is $C^{r}$ for any $r<+\infty$,
and $M_{\varepsilon}$ is $C^{r}O\left(\varepsilon\right)$ close to $M_{0}$. In
addition, there exist perturbed local stable and unstable manifolds of
$M_{\varepsilon}$. They are unions of invariant families of stable and
unstable fibers of dimensions l and k, respectively, and they are
$C^{r}O\left(\varepsilon\right)$ close for all $r<+\infty$ , to their
counterparts.

###### Proof.

For proof of this theorem see Fenichel Fenichel (1971, 1974, 1977, 1979). ∎

The label slow manifold is attached to $M_{\varepsilon}$ because the magnitude
of the vector field restricted to $M_{\varepsilon}$ is
$O\left(\varepsilon\right)$, in terms of the fast independent variable $t$. So
persistent manifolds are labeled slow manifolds, and the proof of their
persistence is carried out by demonstrating that the local stable and unstable
manifolds of $M_{0}$ also persist as locally invariant manifolds in the
perturbed system, i.e., that the local hyperbolic structure persists, and then
the slow manifold is immediately at hand as a locally invariant manifold in
the transverse intersection of these persistent local stable and unstable
manifolds.

###  II.3 Geometric Singular Perturbation Method

Earliest geometric approaches to singularly perturbed systems have been
developed by Cole Cole (1968), O’Malley O’Malley (1974, 1991), Fenichel
Fenichel (1971, 1974, 1977, 1979) for the determination of the slow manifold
equation. Generally, Fenichel theory enables to turn the problem for
explicitly finding functions
$\vec{y}=\vec{Y}\left({\vec{x},\varepsilon}\right)$ whose graphs are locally
invariant slow manifolds $M_{\varepsilon}$ of system (1) into regular
perturbation problem (O’Malley, 1974, p. 112). Invariance of the manifold
$M_{\varepsilon}$ implies that $\vec{Y}\left({\vec{x},\varepsilon}\right)$
satisfies:

$D_{\vec{x}}\vec{Y}\left({\vec{x},\varepsilon}\right)\vec{f}(\vec{x},{\vec{Y}\left({\vec{x},\varepsilon}\right),\varepsilon})=\varepsilon\vec{g}(\vec{x},{\vec{Y}\left({\vec{x},\varepsilon}\right),\varepsilon}).$
(5)

According to Guckenheimer et al. (Guckhenheimer & Holmes, 1983, p.131), this
(partial) differential equation for
$\vec{Y}\left({\vec{x},\varepsilon}\right)$ cannot be solved exactly. So, its
solution can be approximated arbitrarily closely as a Taylor series at
$\left({\vec{x},\varepsilon}\right)=(\vec{0},0)$. Then, the following
perturbation expansion is plugged:

$\vec{Y}\left({\vec{x},\varepsilon}\right)=\vec{Y}_{0}\left(\vec{x}\right)+\varepsilon\vec{Y}_{1}\left(\vec{x}\right)+O\left({\varepsilon^{2}}\right)$
(6)

into (5) to solve order by order for
$\vec{Y}\left({\vec{x},\varepsilon}\right)$. The Taylor series expansion
O’Malley (1991) for
$\vec{f}(\vec{x},{\vec{Y}\left({\vec{x},\varepsilon}\right),\varepsilon})$ and
$\vec{g}(\vec{x},{\vec{Y}\left({\vec{x},\varepsilon}\right),\varepsilon})$ up
to terms of order two in $\varepsilon$ reads:

$\displaystyle\vec{f}(\vec{x},{\vec{Y}\left({\vec{x},\varepsilon}\right),\varepsilon})=\vec{f}(\vec{x},{\vec{Y}_{0}\left(\vec{x}\right),0})+\varepsilon\left[D_{\vec{y}}\vec{f}(\vec{x},{\vec{Y}_{0}\left(\vec{x}\right),0})\vec{Y}_{1}\left(\vec{x}\right)+\frac{\partial\vec{f}}{\partial\varepsilon}(\vec{x},{\vec{Y}_{0}\left(\vec{x}\right),0})\right],$
$\displaystyle\vec{g}(\vec{x},{\vec{Y}\left({\vec{x},\varepsilon}\right),\varepsilon})=\vec{g}(\vec{x},{\vec{Y}_{0}\left(\vec{x}\right),0})+\varepsilon\left[D_{\vec{y}}\vec{g}(\vec{x},{\vec{Y}_{0}\left(\vec{x}\right),0})\vec{Y}_{1}\left(\vec{x}\right)+\frac{\partial\vec{g}}{\partial\varepsilon}(\vec{x},{\vec{Y}_{0}\left(\vec{x}\right),0})\right].$

• At order $\varepsilon^{0}$, Eq. (5) gives:

$D_{\vec{x}}\vec{Y_{0}}\left(\vec{x}\right)\vec{f}(\vec{x},{\vec{Y}_{0}\left(\vec{x}\right),0})=\vec{0},$
(7)

which defines $\vec{Y}_{0}\left(\vec{x}\right)$ due to the invertibility of
$D_{\vec{y}}\vec{f}$ and the Implicit Function Theorem.

• The next order $\varepsilon^{1}$ provides:

$D_{\vec{x}}\vec{Y}_{0}\left(\vec{x}\right)\left[D_{\vec{y}}\vec{f}(\vec{x},{\vec{Y}_{0}\left(\vec{x}\right),0})\vec{Y_{1}}\left(\vec{x}\right)+\frac{\partial\vec{f}}{\partial\varepsilon}\right]=\vec{g}(\vec{x},{\vec{Y}_{0}\left(\vec{x}\right),0}),$
(8)

which yields $\vec{Y}_{1}\left(\vec{x}\right)$ and so forth.

So, regular perturbation theory makes it possible to build an approximation of
locally invariant slow manifolds $M_{\varepsilon}$. Thus, in the framework of
the Geometric Singular Perturbation Method, three conditions are needed to
characterize the slow manifold associated with singularly perturbed system:
existence, local invariance and determination. Existence and local invariance
of the slow manifold are stated according to Fenichel theorem for compact
manifolds with boundary while asymptotic expansions provide its equation up to
the order of the expansion.

####  II.3.1 Van der Pol system

The Van der Pol system Van der Pol (1926) was introduced in the middle of the
twenties for modelling relaxation oscillations. This model presented below is
nowadays considered as a paradigm of singularly perturbed dynamical systems.

$\left\\{\begin{aligned}
x^{\prime}&=f\left(x,y,\varepsilon\right)=x+y-\dfrac{x^{3}}{3},\hfill\\\
y^{\prime}&=g\left(x,y,\varepsilon\right)=-\varepsilon x.\hfill\\\
\end{aligned}\right.$ (9)

According to Eq. (5) invariance of the manifold $M_{\varepsilon}$ reads:

$\left(\dfrac{\partial Y}{\partial
x}(x,\varepsilon)\right)f\left(x,Y(x,\varepsilon),\varepsilon\right)=\varepsilon
g\left(x,Y(x,\varepsilon),\varepsilon\right).$ (10)

By application of the Geometric Singular Perturbation Method and by solving
equation (10) order by order provides at:

Order $\varepsilon^{0}$

$\dfrac{\partial Y_{0}}{\partial
x}(x)f(x,Y_{0}(x),0)=0\quad\Leftrightarrow\quad
f(x,Y_{0}(x),0)=0\quad\Leftrightarrow\quad Y_{0}(x)=\dfrac{x^{3}}{3}-x.$ (11)

Order $\varepsilon^{1}$

$\dfrac{\partial Y_{0}}{\partial
x}(x)Y_{1}(x)=g(x,Y_{0}(x),0)\quad\Leftrightarrow\quad
Y_{1}(x)=\dfrac{x}{1-x^{2}}.$ (12)

Thus, according to Eq. (6), and following this method an approximation up to
order $\varepsilon^{3}$ of the slow invariant manifold of the Van der Pol
singularly perturbed dynamical system is given by:

$y=Y(x,\varepsilon)=\dfrac{x^{3}}{3}-x+\varepsilon\dfrac{x}{1-x^{2}}+\varepsilon^{2}\dfrac{x(1+x^{2})}{(1-x^{2})^{4}}+O\left({\varepsilon^{3}}\right).$
(13)

###  II.4 Successive Approximations Method

At the end of the 1980s, Rossetto Rossetto (1986, 1987) proposed the
Successive Approximations Method to approximate the slow invariant manifold of
singularly perturbed dynamical systems. He considered that the singular
approximation, i.e., the zero-order approximation in $\varepsilon$ of the slow
invariant manifold, i.e., the critical manifold $M_{0}$ defined by
$\vec{f}\left({\vec{x},\vec{y},\varepsilon}\right)=\vec{0}$ is overflowing
invariant with respect to (1) when $\varepsilon=0$. Invariance of the manifold
$M_{0}$ implies that $\vec{f}\left({\vec{x},\vec{y},\varepsilon}\right)$
satisfies:

$\dfrac{d\vec{f}}{dt}\left({\vec{x},\vec{y},\varepsilon}\right)=\vec{0}\quad\Leftrightarrow\quad(D_{\vec{x}}\vec{f})\vec{x}^{\prime}+(D_{\vec{y}}\vec{f})\vec{y}^{\prime}=\vec{0}.$
(14)

According to (1), this leads to:

$(D_{\vec{x}}\vec{f})\vec{f}({\vec{x},\vec{y},\varepsilon})+\varepsilon(D_{\vec{y}}\vec{f})\vec{g}({\vec{x},\vec{y},\varepsilon})=\vec{0}.$
(15)

By application of the Implicit Function Theorem, we have:

$\vec{f}(\vec{x},\vec{y},\varepsilon)=\vec{0}\quad\Leftrightarrow\quad\vec{y}=\vec{Y}\left({\vec{x},\varepsilon}\right),$
(16)

implies that:

$D_{\vec{x}}\vec{Y}=-(D_{\vec{y}}\vec{f})^{-1}(D_{\vec{x}}\vec{f}).$ (17)

Left multiplication of (15) by $(D_{\vec{y}}\vec{f})^{-1}$ gives:

$(D_{\vec{y}}\vec{f})^{-1}(D_{\vec{x}}\vec{f})\vec{f}({\vec{x},\vec{y},\varepsilon})+\varepsilon\vec{g}({\vec{x},\vec{y},\varepsilon})=\vec{0}.$
(18)

Taking into account (16) & (17), we find the following equation:

$(D_{\vec{x}}\vec{Y})\vec{f}({\vec{x},\vec{Y}(\vec{x},\varepsilon),\varepsilon})=\varepsilon\vec{g}({\vec{x},\vec{Y}(\vec{x},\varepsilon),\varepsilon}),$
(19)

which is absolutely identical to (5). Thus, identity between Geometric
Singular Perturbation Method and Successive Approximations Method is stated.

###  II.5 Zero-Derivative Principle

In 2005, an iterative model reduction method called Zero-Derivative Principle
(ZDP) was presented by Gear et al. Gear et al. (2005). Four years later,
Zagaris et al. Zagaris et al. (2009) applied the ZDP to singularly perturbed
dynamical systems (1) to approximate their slow invariant manifold. They
considered that if a function is locally invariant under (1) the zero set of
its time derivative provides an even more accurate approximation of the slow
invariant manifold (SIM) of (1). In other words, it means that the $k^{th}$
approximation of the slow invariant manifold of (1) is given by the following
equation:

$\dfrac{d^{k+1}\vec{x}}{dt^{k+1}}=\dfrac{d^{k}\vec{f}}{dt^{k}}\left({\vec{x},\vec{y},\varepsilon}\right)=\vec{0}.$
(20)

According to Benoît et al. Benoît et al. (2015): “estimates [of the SIM] can
be obtained in terms of approximations of the slow manifold to order
$\varepsilon^{k}$ when looking at the zeros of the $k^{th}$ derivative of the
fast components of the original singularly perturbed system.” First of all,
let’s notice that such result has been already stated and published by
Rossetto Rossetto (1986, 1987) and has not been quoted in any reference. In
his works, Rossetto Rossetto (1986, 1987) presented his Successive
Approximations Method and stated that:

> “The $m^{th}$ order approximation of slow motion is solution of the
> dynamical system $\vec{f}^{(m)}=\vec{0}$.”

Thus, according to the ZDP, if
$\vec{x}^{\prime}=\vec{f}\left({\vec{x},\vec{y},\varepsilon}\right)=\vec{0}$
provides the zero-order approximation in $\varepsilon$ of the slow invariant
manifold of (1), $\vec{x}^{\prime\prime}=\vec{0}$ will define its first-order
approximation in $\varepsilon$. So, according to (1), we have the following
equation:

$\vec{x}^{\prime\prime}=(D_{\vec{x}}\vec{f})\vec{x}^{\prime}+(D_{\vec{y}}\vec{f})\vec{y}^{\prime}=\vec{0},$
(21)

which reads:

$(D_{\vec{x}}\vec{f})\vec{f}\left({\vec{x},\vec{y},\varepsilon}\right)+\varepsilon(D_{\vec{y}}\vec{f})\vec{g}\left({\vec{x},\vec{y},\varepsilon}\right)=\vec{0}.$
(22)

This equation (22) is absolutely identical to Eq. (15). This proves that both
Successive Approximations Method and Zero-Derivative Principle are identical.
Moreover, since we have already stated above that both Geometric Singular
Perturbation Method and Successive Approximations Method are identical, it
follows by transitivity that these three methods are absolutely identical and
so belong to the first category called: Singular Perturbation-Based Methods.

##  III Curvature-Based Methods

As previously recalled, beside, singularly perturbed dynamical systems (1)
there are dynamical systems without an explicit timescale splitting. In this
case the Singular Perturbation-Based Methods, i.e., GSPM, SAM and ZDP can not
be applied anymore. Some of these systems, called slow-fast dynamical systems,
may be defined by the following set of nonlinear ordinary differential
equations:

$\frac{d\vec{X}}{dt}=\overrightarrow{\Im}(\vec{X}),$ (23)

with $\vec{X}=\left[{x_{1},x_{2},...,x_{n}}\right]^{t}\in
E\subset{\mathbb{R}}^{n}$ and
$\overrightarrow{\Im}(\vec{X})=\left[{f_{1}(\vec{X}),f_{2}(\vec{X}),...,f_{n}(\vec{X})}\right]^{t}\in
E\subset{\mathbb{R}}^{n}$. The vector $\overrightarrow{\Im}(\vec{X})$ defines
a velocity vector field $\dot{\vec{X}}=\vec{V}$ in $E$ whose components
$f_{i}$ which are supposed to be continuous and infinitely differentiable with
respect to all $x_{i}$ and $t$, i.e. are $C^{\infty}$ functions in $E$ and
with values included in $\mathbb{R}$, satisfy the assumptions of the Cauchy-
Lipschitz theorem. For more details, see for example Coddington & Levinson
Coddington & Levinson (1955). A solution of this system is a trajectory curve
$\vec{X}\left(t\right)$ tangent (Except at the fixed points) to
$\overrightarrow{\Im}$ whose values define the states of the dynamical system
described by the Eq. (23). According to Ginoux Ginoux (2009), non-singularly
perturbed dynamical systems (23) can be considered as slow-fast dynamical
systems if their Jacobian matrix has at least one fast eigenvalue, i.e. with
the largest absolute value of the real part over a huge domain of the phase
space.

###  III.1 Intrinsic Low-Dimensional Manifold Method

In 1992, Maas and Pope introduced the Intrinsic Low-Dimensional Manifold
(ILDM) method Maas & Pope (1992). Within the framework of application of the
Tikhonov’s theorem Tikhonov (1952), this method uses the fact that in the
vicinity of the slow manifold the eigenmode associated with the fast
eigenvalue is evanescent. According to Mass and Pope Maas & Pope (1992):

> “If a subspace in the composition space can be found where the system is “in
> equilibrium” with respect to its smallest $m$ eigenvectors, this subspace
> defines a low-dimensional manifold that is characterized by the fact that
> movements along it are associated with slow time scales and can be used to
> simplify chemical kinetics.”

Thus, IDLM method provides an approximation of the slow invariant manifold of
the slow-fast dynamical system (23) by considering that its velocity vector
field is coplanar to the slow eigenvectors of the Jacobian of
$\overrightarrow{\Im}(\vec{X})$.

###  III.2 Tangent Linear System Approximation Method

At the end of the nineties, Rossetto et al. Rossetto et al. (1998) presented
their Tangent Linear System Approximation Method. The application of this
method requires that the slow-fast dynamical system (23) satisfies the
following assumptions:

  * (H1)

The components $f_{i}$, of the velocity vector field
$\overrightarrow{\Im}(\vec{X})$ defined in $E$ are continuous, $C^{\infty}$
functions in $E$ and with values included in $\mathbb{R}$.

  * (H2)

The slow-fast dynamical system (23) satisfies the nonlinear part condition,
i.e., that the influence of the nonlinear part of the Taylor series of the
velocity vector field $\overrightarrow{\Im}(\vec{X})$ of this system is
overshadowed by the fast dynamics of the linear part.

$\overrightarrow{\Im}(\vec{X})=\overrightarrow{\Im}(\vec{X}_{0})+(\vec{X}-\vec{X}_{0})\left.D_{\vec{X}}\overrightarrow{\Im}(\vec{X})\right|_{\vec{X}_{0}}+O\left((\vec{X}-\vec{X}_{0})^{2}\right).$
(24)

To the slow-fast dynamical system (23) is associated a tangent linear system
defined as follows:

$\frac{d\delta\vec{X}}{dt}=J(\vec{X}_{0})\delta\vec{X},$ (25)

where $\delta\vec{X}=\vec{X}-\vec{X}_{0}$, $\vec{X}_{0}=\vec{X}(t_{0})$ and
$J(\vec{X}_{0})=\left.D_{\vec{X}}\overrightarrow{\Im}(\vec{X})\right|_{\vec{X}_{0}}$
is the Jacobian matrix. Thus, the solution of the tangent linear system (25)
reads:

$\delta\vec{X}=e^{J\left({\vec{X}_{0}}\right)\left({t-t_{0}}\right)}\delta\vec{X}\left({t_{0}}\right).$
(26)

So,

$\delta\vec{X}=\sum\limits_{i=1}^{n}{a_{i}}\vec{Y}_{\lambda_{i}}.$ (27)

where $n$ is the dimension of the eigenspace, $a_{i}$ represents coefficients
depending explicitly on the co-ordinates of space and implicitly on time and
$\vec{Y}_{\lambda_{i}}$ the eigenvectors associated with the Jacobian matrix
of the tangent linear system. Thus, according to the TLSA method we have the
following proposition:

######  Proposition 1.

In the vicinity of the slow manifold the velocity of the slow-fast dynamical
system (23) and that of the tangent linear system (25) merge.

$\frac{d\delta\vec{X}}{dt}=\vec{V}_{T}\approx\vec{V},$ (28)

where $\vec{V}_{T}$ represents the velocity vector associated with the tangent
linear system.

Hence, the Tangent Linear System Approximation method consists in the
projection of the velocity vector field $\vec{V}$ on the eigenbasis associated
with the eigenvalues of the functional Jacobian of the tangent linear system.
Indeed, by taking account of (25) and (27) we have according to (28):

$\frac{d\delta\vec{X}}{dt}=J(\vec{X}_{0})\delta\vec{X}=J(\vec{X}_{0})\sum\limits_{i=1}^{n}{a_{i}}\vec{Y}_{\lambda_{i}}=\sum\limits_{i=1}^{n}{a_{i}}J(\vec{X}_{0})\vec{Y}_{\lambda_{i}}=\sum\limits_{i=1}^{n}{a_{i}}\lambda_{i}\vec{Y}_{\lambda_{i}}.$
(29)

Thus, Proposition 1 leads to:

$\frac{d\delta\vec{X}}{dt}=\vec{V}_{T}\approx\vec{V}=\sum\limits_{i=1}^{n}{a_{i}}\lambda_{i}\vec{Y}_{\lambda_{i}}.$
(30)

The equation (30) constitutes in dimension two (resp. dimension three) a
condition called collinearity (resp. coplanarity) condition which provides the
analytical equation of the slow manifold of the slow-fast dynamical system
(23). An alternative proposed by Rossetto et al. Rossetto et al. (1998) uses
the “fast” eigenvector on the left associated with the “fast” eigenvalue of
the transposed functional Jacobian of the tangent linear system. In this case
the velocity vector field $\vec{V}$ is then orthogonal with the “fast”
eigenvector on the left. This constitutes a condition called orthogonality
condition which provides the analytical equation of the slow manifold of the
slow-fast dynamical system (23) (See Rossetto et al. Rossetto et al. (1998)
and Ginoux et al. Ginoux & Rossetto (2006a); Ginoux et al. (2008)).

For a two-dimensional slow-fast dynamical system (23) the projection of the
velocity vector field $\vec{V}$ on the eigenbasis reads thus:

$\frac{d\delta\vec{X}}{dt}=\vec{V}_{T}\approx\vec{V}=\sum\limits_{i=1}^{2}{a_{i}}\lambda_{i}\vec{Y}_{\lambda_{i}}=\alpha\vec{Y}_{\lambda_{1}}+\beta\vec{Y}_{\lambda_{2}},$

where $\alpha$ and $\beta$ represent coefficients depending explicitly on the
co-ordinates of space and implicitly on time and where $\vec{Y}_{\lambda_{1}}$
represents the “fast” eigenvector and $\vec{Y}_{\lambda_{2}}$ the “slow”
eigenvector. The existence of an evanescent mode in the vicinity of the slow
manifold implies according to Tikhonov’s theorem Tikhonov (1952): $\alpha\ll
1$. We deduce:

$\vec{V}\approx\beta\vec{Y}_{\lambda_{2}}.$

So, the following collinearity condition between the velocity vector field
$\vec{V}$ of a two-dimensional slow-fast dynamical system (23) and the “slow”
eigenvector $\vec{Y}_{\lambda_{2}}$ provides an approximation of its slow
manifold:

$\vec{V}\wedge\vec{Y}_{\lambda_{2}}=\vec{0}.$ (31)

For a three-dimensional slow-fast dynamical system (23) the projection of the
velocity vector field $\vec{V}$ on the eigenbasis reads thus:

$\frac{d\delta\vec{X}}{dt}=\vec{V}_{T}\approx\vec{V}=\sum\limits_{i=1}^{3}{a_{i}}\lambda_{i}\vec{Y}_{\lambda_{i}}=\alpha\vec{Y}_{\lambda_{1}}+\beta\vec{Y}_{\lambda_{2}}+\delta\vec{Y}_{\lambda_{3}},$

where $\alpha$, $\beta$ and $\delta$ represent coefficients depending
explicitly on the co-ordinates of space and implicitly on time and where
$\vec{Y}_{\lambda_{1}}$ represents the “fast” eigenvector and
$\vec{Y}_{\lambda_{2}},\vec{Y}_{\lambda_{3}}$ the “slow” eigenvectors. The
existence of an evanescent mode in the vicinity of the slow manifold implies
according to Tikhonov’s theorem Tikhonov (1952): $\alpha\ll 1$. We deduce:

$\vec{V}\approx\beta\vec{Y}_{\lambda_{2}}+\delta\vec{Y}_{\lambda_{3}}.$

So, the following coplanarity condition between the velocity vector field
$\vec{V}$ of a three-dimensional slow-fast dynamical system (23) and the
“slow” eigenvectors $\vec{Y}_{\lambda_{2}}$ and $\vec{Y}_{\lambda_{3}}$
provides an approximation of its slow manifold:

$\vec{V}.(\vec{Y}_{\lambda_{2}}\wedge\vec{Y}_{\lambda_{3}})=0.$ (32)

This brief presentation enables to state that both Intrinsic Low-Dimensional
Manifold method and Tangent Linear System Approximation method are absolutely
identical. Let’s notice that these methods can be also applied to singularly
perturbed dynamical systems (1). However, they presented a major drawback
since they required the computation of eigenvalues and eigenvectors of the
Jacobian which could only be done numerically for dynamical systems of
dimension greater than two. Moreover, according to the nature of the “slow”
eigenvalues (real or complex conjugated) the plot of their slow manifold
analytical equation was difficult even impossible. So, to solve this problem
it was necessary to make the slow manifold analytical equation independent of
the “slow” eigenvalues. That’s was the aim of the Flow Curvature Method.

###  III.3 Flow Curvature Method

Fifteen years ago, Ginoux et al. Ginoux & Rossetto (2006a) published the first
article in which the Flow Curvature Method was presented. They explained that
by multiplying the slow manifold analytical equation of a two dimensional
dynamical system by a “conjugated” equation, that of a three dimensional
dynamical system by two “conjugated” equations make the slow manifold
analytical equation independent of the “slow” eigenvalues of the tangent
linear system.

For a two-dimensional slow-fast dynamical system (23), let’s multiply equation
(31) by its “conjugated” equation, i.e., by an equation in which the
eigenvalue $\lambda_{2}$ is replaced by the eigenvalue $\lambda_{1}$. Let’s
notice that the “conjugated” equation of the equation (31) corresponds to the
collinearity condition between the velocity vector field $\vec{V}$ and the
fast eigenvector $\vec{Y}_{\lambda_{1}}$. The product of the equation (31) by
its “conjugated” equation reads:

$(\vec{V}\wedge\vec{Y}_{\lambda_{1}})\cdot(\vec{V}\wedge\vec{Y}_{\lambda_{2}})=0.$
(33)

It has been proved in Ginoux et al. Ginoux & Rossetto (2006a) that:

$(\vec{V}\wedge\vec{Y}_{\lambda_{1}})\cdot(\vec{V}\wedge\vec{Y}_{\lambda_{2}})=\left\|\vec{\gamma}\wedge\vec{V}\right\|^{2}=0,$
(34)

where $\vec{\gamma}=\dot{\vec{V}}=\ddot{\vec{X}}$ represents the acceleration
vector field, i.e., the time derivative of the velocity vector field (23).
But, in the framework of Differential Geometry, curvature of trajectory curve
integral of slow-fast dynamical system (23) is defined by:

$\dfrac{1}{\mathcal{R}}=\dfrac{\left\|\vec{\gamma}\wedge\vec{V}\right\|}{\left\|\vec{V}\right\|^{3}}$
(35)

where $\mathcal{R}$ represents the radius of curvature. Thus, it followed from
this result that the slow manifold of a two-dimensional slow-fast dynamical
system (23) can be directly approximated from the collinearity condition
between the velocity vector field $\vec{V}=\dot{\vec{X}}$ (23) and its
acceleration vector field $\vec{\gamma}=\dot{\vec{V}}=\ddot{\vec{X}}$. More
than ten years ago, Ginoux et al. Ginoux et al. (2008) proved that this
collinearity condition can be written as:

$\det(\dot{\vec{X}},\ddot{\vec{X}})=0.$ (36)

This led Ginoux et al. Ginoux et al. (2008) to the following proposition:

######  Proposition 2.

The location of the points where the local curvature of the trajectory curves
integral of a two-dimensional dynamical system defined by (23) vanishes,
directly provides the slow manifold analytical equation associated to this
system.

For a three-dimensional slow-fast dynamical system (23), equation (32) must be
multiplied by two “conjugated” equations obtained by circular shifts of the
eigenvalues. The product of equation (32) by its “conjugated” equation reads:

$\left[\vec{V}\cdot(\vec{Y}_{\lambda_{1}}\wedge\vec{Y}_{\lambda_{2}})\right]\cdot\left[\vec{V}\cdot(\vec{Y}_{\lambda_{2}}\wedge\vec{Y}_{\lambda_{3}})\right]\cdot\left[\vec{V}\cdot(\vec{Y}_{\lambda_{1}}\wedge\vec{Y}_{\lambda_{3}})\right]=0.$
(37)

It has been proved in Ginoux et al. Ginoux & Rossetto (2006a) that:

$\left[\vec{V}\cdot(\vec{Y}_{\lambda_{1}}\wedge\vec{Y}_{\lambda_{2}})\right]\cdot\left[\vec{V}\cdot(\vec{Y}_{\lambda_{2}}\wedge\vec{Y}_{\lambda_{3}})\right]\cdot\left[\vec{V}\cdot(\vec{Y}_{\lambda_{1}}\wedge\vec{Y}_{\lambda_{3}})\right]=\dot{\vec{\gamma}}\cdot(\vec{\gamma}\wedge\vec{V})=0,$
(38)

where $\dot{\vec{\gamma}}$ represents the over-acceleration vector field. But,
in the framework of Differential Geometry, torsion of trajectory curve
integral of slow-fast dynamical system (23) is defined by:

$\frac{1}{\mathcal{T}}=-\frac{\dot{\vec{\gamma}}\cdot\left({\vec{\gamma}\wedge\vec{V}}\right)}{\left\|{\vec{\gamma}\wedge\vec{V}}\right\|^{2}},$
(39)

where $\mathcal{T}$ represents the radius of torsion. Thus, it followed from
this result that the slow manifold of a three-dimensional slow-fast dynamical
system (23) can be directly approximated from the coplanarity condition
between the velocity vector field $\vec{V}=\dot{\vec{X}}$ (23), its
acceleration vector field $\vec{\gamma}$ and its over-acceleration
$\dot{\vec{\gamma}}$. More than ten years ago, Ginoux et al. Ginoux et al.
(2008) proved that this coplanarity condition can be written as:

$\det(\dot{\vec{X}},\ddot{\vec{X}},\dddot{\vec{X}})=0.$ (40)

This led Ginoux et al. Ginoux et al. (2008) to the following proposition:

######  Proposition 3.

The location of the points where the local torsion of the trajectory curves
integral of a three-dimensional dynamical system defined by (23) vanishes,
directly provides the slow manifold analytical equation associated to this
system.

Then, Ginoux et al. Ginoux et al. (2008) and Ginoux Ginoux (2009) generalized
these results to any $n$-dimensional slow-fast dynamical system (23) as well
as any $n$-dimensional singularly perturbed dynamical system (1). This led to
the following general proposition which encompasses proposition 2 & 3:

######  Proposition 4.

The location of the points where the curvature of the flow, i.e., the
curvature of the trajectory curves of any $n$-dimensional dynamical system (1)
or (23) vanishes, directly provides its $(n-1)$-dimensional slow invariant
manifold analytical equation which reads:

$\phi(\vec{X})=\dot{\vec{X}}\cdot(\ddot{\vec{X}}\wedge\dddot{\vec{X}}\wedge\ldots\wedge\mathop{\vec{X}}\limits^{\left(n\right)})=det(\dot{\vec{X}},\ddot{\vec{X}},\dddot{\vec{X}},\ldots,\mathop{\vec{X}}\limits^{\left(n\right)})=0.$
(41)

###### Proof.

For proof of this proposition 4, see Ginoux et al. Ginoux et al. (2008) and
Ginoux Ginoux (2009). ∎

Let’s notice that within the framework of Differential Geometry,
$n$-dimensional smooth curves, i.e., smooth curves in Euclidean $n-$space are
defined by a regular parametric representation in terms of arc length also
called natural representation or unit speed parametrization. According to
Herman Gluck Gluck (1966) local metrics properties of curvatures may be
directly deduced from curves parametrized in terms of time and so natural
representation is not necessary. This fundamental result allows applying the
Grahm-Schmidt orthogonalization process, with a collection of $i=1,\ldots,n$
vectors $\vec{u}_{i}(t)$ depending explicitly on time (not on arc length) and
forming an orthogonal basis, to deduce the expression of curvatures of
trajectory curves integral of dynamical systems (1) or (23) in Euclidean
$n-$space.

Then, to prove that the slow manifold (41) is locally invariant by the flow of
(23), we use Darboux invariance theorem. According to Schlomiuk Schlomiuk
(1993) and Llibre et al. Llibre & Medrado (2007) it seems that Gaston Darboux
Darboux (1878) has been the first to define the concept of invariant manifold
which can be presented as follows:

######  Proposition 5.

The manifold defined by $\phi(\vec{X})=0$ where $\phi$ is a $C^{1}$ in an open
set U is invariant with respect to the flow of (1) or (23) if there exists a
$C^{1}$ function denoted $k(\vec{X})$ and called cofactor which satisfies:

$L_{\vec{V}}\phi(\vec{X})=k(\vec{X})\phi(\vec{X}),$ (42)

for all $\vec{X}\in U$ and with the Lie derivative operator defined as:

$L_{\vec{V}}\phi=\vec{V}\cdot\vec{\nabla}\phi=\sum\limits_{i=1}^{n}{\frac{\partial\phi}{\partial
x_{i}}\dot{x}_{i}}=\frac{d\phi}{dt}.$

###### Proof.

To state the invariance of slow manifold (42), it is first necessary to recall
the following identities (43)-(44) the proof of which is obvious:

$\ddot{\vec{X}}=J\dot{\vec{X}},$ (43)

where $J$ is the Jacobian matrix associated with any $n$-dimensional slow-fast
dynamical system (23).

$J\vec{a}_{1}.\left({\vec{a}_{2}\wedge\ldots\wedge\vec{a}_{n}}\right)+\vec{a}_{1}.\left({J\vec{a}_{2}\wedge\ldots\wedge\vec{a}_{n}}\right)+\ldots+\vec{a}_{1}.\left({\vec{a}_{2}\wedge\ldots\wedge
J\vec{a}_{n}}\right)=Tr\left(J\right)\vec{a}_{1}.\left({\vec{a}_{2}\wedge\ldots\wedge\vec{a}_{n}}\right).$
(44)

Thus, the Lie derivative of the slow manifold (42) reads:

$L_{\vec{V}}\phi({\vec{X}})=\dot{\vec{X}}\cdot({\ddot{\vec{X}}\wedge\dddot{\vec{X}}\wedge\ldots\wedge\mathop{\vec{X}}\limits^{\left({n+1}\right)}}),$
(45)

From the identity (43) we find that:

$\mathop{\vec{X}}\limits^{\left({n+1}\right)}=J^{n}\dot{\vec{X}}\mbox{\quad if
\quad}\frac{dJ}{dt}=0,$ (46)

where $J^{n}$ represents the $n^{th}$ power of $J$. Then, it follows that:

$\mathop{\vec{X}}\limits^{\left({n+1}\right)}=J\mathop{\vec{X}}\limits^{\left(n\right)}.$
(47)

Replacing $\mathop{\vec{X}}\limits^{\left({n+1}\right)}$ in Eq. (45) with Eq.
(47) we have:

$L_{\vec{V}}\phi({\vec{X}})=\dot{\vec{X}}\cdot({\ddot{\vec{X}}\wedge\dddot{\vec{X}}\wedge\ldots\wedge
J\mathop{\vec{X}}\limits^{\left(n\right)}}).$ (48)

The right hand side of this Eq. (48) can be written:

$J\dot{\vec{X}}\cdot({\ddot{\vec{X}}\wedge\dddot{\vec{X}}\wedge\ldots\wedge\mathop{\vec{X}}\limits^{\left(n\right)}})+\dot{\vec{X}}\cdot(J{\ddot{\vec{X}}\wedge\dddot{\vec{X}}\wedge\ldots\wedge\mathop{\vec{X}}\limits^{\left(n\right)}})+\ldots+\dot{\vec{X}}\cdot({\ddot{\vec{X}}\wedge\dddot{\vec{X}}\wedge\ldots\wedge
J\mathop{\vec{X}}\limits^{\left(n\right)}}).$

According to Eq. (47) all terms are null except the last one. So, by taking
into account identity (44), we find:

$L_{\vec{V}}\phi({\vec{X}})=Tr\left(J\right)\dot{\vec{X}}\cdot({\ddot{\vec{X}}\wedge\dddot{\vec{X}}\wedge\ldots\wedge\mathop{\vec{X}}\limits^{\left(n\right)}})=Tr\left(J\right)\phi({\vec{X}})=k({\vec{X}})\phi({\vec{X}}),$

where $k({\vec{X}})=Tr\left(J\right)$ represents the trace of the Jacobian
matrix.

So, the invariance of the slow manifold analytical equation of any
$n$-dimensional slow-fast dynamical system (23) is established provided that
the Jacobian matrix is locally stationary. ∎

More than ten years ago, Ginoux Ginoux (2009) established for any
$n$-dimensional singularly perturbed dynamical systems (1) the identity
between Fenichel’s invariance (5) and Darboux invariance (42) (see also Ginoux
Ginoux (2014)).

Thus, it follows from what precedes that both Tangent Linear System
Approximation Method and Flow Curvature Method are absolutely identical.

###  III.4 Inflection Line Method

In 1994, Brøns and Bar-Eli Brøns & Bar-Eli (1994) presented their Inflection
Line Method for planar or two-dimensional slow-fast dynamical system (23)
which can be written as follows:

$\left\\{\begin{aligned}
\varepsilon\dot{x}&=f\left(x,y,\varepsilon\right),\hfill\\\
\dot{y}&=g\left(x,y,\varepsilon\right).\hfill\\\ \end{aligned}\right.$ (49)

According to Brøns and Bar-Eli Brøns & Bar-Eli (1994):

> “The curvature $\kappa$ of a smooth curve measures the rate of turn of the
> tangent vector with respect to arc length $s$. The inflection line $l_{i}$
> associated with (49) is the locus of points where $\kappa_{t}=0$. The
> inflection line typically consists of one or more smooth curves. For our
> purposes, it suffices to consider systems (49) in a region where
> $f(x,y,\varepsilon)\neq 0$. Then time can be eliminated to give the
> equivalent equation for the trajectories,
>
>
> $\dfrac{dy}{dx}=\dfrac{g\left(x,y,\varepsilon\right)}{f\left(x,y,\varepsilon\right)}.$
> (50)
>
> For curves in the ($x,y$)-plane of the form $y=y(x)$, the curvature can be
> computed as
>
> $\kappa=\dfrac{y^{\prime\prime}(x)}{\left(1+y^{\prime}(x)\right)^{3/2}},$
> (51)
>
> where the positive root in the denominator is used. Thus, the inflection
> line for (50) is determined by $y^{\prime\prime}(x)=0$, and an equation for
> it can be found by differentiating (50).”

For such two-dimensional dynamical systems (49), the issue of the
parametrization highlighted by Gluck Gluck (1966) just above can be easily
solved by dividing the two components of the velocity vector field (49) as
done by Brøns and Bar-Eli Brøns & Bar-Eli (1994). Nevertheless, such
simplification is no more possible in dimension higher than two. Let’s notice
that for three-dimensional dynamical systems the use of curvature instead of
torsion has already been carried out by Ginoux et al. Ginoux & Rossetto
(2006b) and by Gilmore et al. Gilmore et al. (2010). However, the use of
curvature or more precisely, the first curvature require to solve a system of
two nonlinear equations since in dimension three a curve is defined as the
intersection of two surfaces. Nowadays, this kind of problem cannot be solved
analytically in the general case. Thus, such considerations precludes any
generalization of the Inflection Line Method proposed by Brøns and Bar-Eli
Brøns & Bar-Eli (1994) at least from an analytical point of view.

Now let’s prove that the Inflection Line Method is just a particular case of
the more general Flow Curvature Method. According to Brøns and Bar-Eli Brøns &
Bar-Eli (1994) “the inflection line for (50) is determined by
$y^{\prime\prime}(x)=0$, and an equation for it can be found by
differentiating (50).” So, let’s differentiate equation (50):

$\dfrac{d}{dx}\left(\dfrac{dy}{dx}\right)=\dfrac{d}{dx}\left(\dfrac{dy}{dt}\dfrac{dt}{dx}\right)=\dfrac{d}{dt}\left(\dfrac{\dot{y}}{\dot{x}}\right)\dfrac{dt}{dx}=\dfrac{\ddot{y}\dot{x}-\ddot{x}\dot{y}}{\dot{x}^{3}}=\dfrac{1}{\dot{x}^{3}}\begin{vmatrix}\dot{x}&\ddot{x}\\\
\dot{y}&\ddot{y}\end{vmatrix}=0.$ (52)

So, the condition $y^{\prime\prime}(x)=0$ given by Eq. (52) is absolutely
identical to the condition $\det(\dot{\vec{X}},\ddot{\vec{X}})=0$ where
$\vec{X}=(x,y)^{t}$ provided by (36). As a consequence, it is thus proved that
the Inflection Line Method is just a particular case of the more general Flow
Curvature Method.

To resume, we have proved in this section that the Flow Curvature Method
encompasses the Intrinsic Low-Dimensional Manifold Method, the Tangent Linear
System Approximation Method and the Inflection Line Method which all belong to
the second category called Curvature-Based Methods.

Now, a last issue arises as follows: for $n$-dimensional singularly perturbed
dynamical systems (1), does the Flow Curvature Method provide the same
approximation of their slow invariant manifold as the Geometric Singular
Perturbation Method?

##  IV Singular Perturbation and Curvature Based Methods

While the slow invariant manifold analytical equation (6) given by the
Geometric Singular Perturbation Theory is an explicit equation, the slow
invariant manifold analytical equation (41) obtained according to the Flow
Curvature Method is an implicit equation. So, in order to compare the latter
with the former it is necessary to plug the following perturbation expansion:
$\vec{Y}\left({\vec{x},\varepsilon}\right)=\vec{Y}_{0}\left(\vec{x}\right)+\varepsilon\vec{Y}_{1}\left(\vec{x}\right)+O\left({\varepsilon^{2}}\right)$
into (41). Thus, solving order by order for
$\vec{Y}\left({\vec{x},\varepsilon}\right)$ will transform (41) into an
explicit analytical equation enabling the comparison with (6). The Taylor
series expansion for
$\phi(\vec{X},\varepsilon)=\phi(\vec{x},\vec{Y}\left({\vec{x},\varepsilon}\right),\varepsilon)$
up to terms of order one in $\varepsilon$ reads:

$\phi(\vec{X},\varepsilon)=\phi(\vec{x},\vec{Y}_{0}\left(\vec{x}\right),0)+\varepsilon
D_{\vec{y}}\phi(\vec{x},\vec{Y}_{0}\left(\vec{x}\right),0)\vec{Y}_{1}\left(\vec{x}\right)+\varepsilon\frac{\partial\phi}{\partial\varepsilon}(\vec{x},\vec{Y}_{0}\left(\vec{x}\right),0).$
(53)

• At order $\varepsilon^{0}$, Eq. (53) gives:

$\phi\left(\vec{x},\vec{Y}_{0}\left(\vec{x}\right),0\right)=0,$ (54)

which defines $\vec{Y}_{0}\left(\vec{x}\right)$ due to the invertibility of
$D_{\vec{y}}\phi$ and application of the Implicit Function Theorem.

• The next order $\varepsilon^{1}$, provides:

$D_{\vec{y}}\phi(\vec{x},\vec{Y}_{0}\left(\vec{x}\right),0)\vec{Y}_{1}\left(\vec{x}\right)+\frac{\partial\phi}{\partial\varepsilon}(\vec{x},\vec{Y}_{0}\left(\vec{x}\right),0)=\vec{0},$
(55)

which yields $\vec{Y}_{1}\left(\vec{x}\right)$ and so forth.

In order to prove that this equation is completely identical to Eq. (8), let’s
rewrite it as follows:

$\vec{Y}_{1}\left(\vec{x}\right)=-\left[D_{\vec{y}}\phi(\vec{x},\vec{Y}_{0}\left(\vec{x}\right),0)\right]^{-1}\frac{\partial\phi}{\partial\varepsilon}(\vec{x},\vec{Y}_{0}\left(\vec{x}\right),0).$

By application of the chain rule, i.e., the derivative of
$\phi(\vec{x},\vec{Y}_{0}\left(\vec{x}\right),0)$ with respect to the variable
$\vec{y}$ and then with respect to $\varepsilon$, it can be stated that:

$\vec{Y}_{1}\left(\vec{x}\right)=-\left[(D_{\vec{x}}\vec{f})(D_{\vec{y}}\vec{f})\right]^{-1}(D_{\vec{y}}\vec{f})\vec{g}(\vec{x},\vec{Y}_{0}(\vec{x}),0)-\left[D_{\vec{y}}\vec{f}\right]^{-1}D_{\varepsilon}\vec{f}(\vec{x},\vec{Y}_{0}\left(\vec{x}\right),0).$

But, according to the Implicit Function Theorem we have:

$(D_{\vec{x}}\vec{f})=-(D_{\vec{y}}\vec{f})(D_{\vec{x}}\vec{y})=-(D_{\vec{y}}\vec{f})(D_{\vec{x}}\vec{Y}_{0}(\vec{x})).$

Then, by replacing into the previous equation we find:

$\vec{Y}_{1}\left(\vec{x}\right)=\left[(D_{\vec{y}}\vec{f})(D_{\vec{x}}\vec{Y}_{0}(\vec{x}))(D_{\vec{y}}\vec{f})\right]^{-1}(D_{\vec{y}}\vec{f})\vec{g}(\vec{x},\vec{Y}_{0}(\vec{x}),0)-\left[D_{\vec{y}}\vec{f}\right]^{-1}D_{\varepsilon}\vec{f}(\vec{x},\vec{Y}_{0}\left(\vec{x}\right),0).$

After simplifications, we have:

$\vec{Y_{1}}\left(\vec{x}\right)=\left[D_{\vec{x}}\vec{Y}_{0}\left(\vec{x}\right)D_{\vec{y}}\vec{f}(\vec{x},{\vec{Y}_{0}\left(\vec{x}\right),0})\right]^{-1}\vec{g}(\vec{x},{\vec{Y}_{0}\left(\vec{x}\right),0})-\left[D_{\vec{y}}\vec{f}(\vec{x},{\vec{Y}_{0}\left(\vec{x}\right),0})\right]^{-1}\frac{\partial\vec{f}}{\partial\varepsilon}.$

Finally, Eq. (55) may be written as:

$D_{\vec{x}}\vec{Y}_{0}\left(\vec{x}\right)\left[D_{\vec{y}}\vec{f}(\vec{x},{\vec{Y}_{0}\left(\vec{x}\right),0})\vec{Y_{1}}\left(\vec{x}\right)+\frac{\partial\vec{f}}{\partial\varepsilon}\right]=\vec{g}(\vec{x},{\vec{Y}_{0}\left(\vec{x}\right),0}).$

Thus, identity between the slow invariant manifold equation given by the
Geometric Singular Perturbation Theory and by the Flow Curvature Method is
proved up to first order term in $\varepsilon$.

Let’s notice that the slow invariant manifold equation (41) associated with
$n$–dimensional singularly perturbed systems defined by the Flow Curvature
Method is a tensor of order $n$. As a consequence, it can only provide an
approximation of $n$-order in $\varepsilon$ of the slow invariant manifold
equation (6). In order to emphasize this result let’s come back to the example
of the Van der Pol system Van der Pol (1926).

###  IV.1 Van der Pol system

Let’s consider the Van der Pol system (9) which is a well-known two-
dimensional singularly perturbed dynamical system. By application of the
Geometric Singular Perturbation Method an approximation up to order
$\varepsilon^{3}$ of the slow invariant manifold of the Van der Pol singularly
perturbed dynamical system (9) is given by Eq. (13):

$y=Y(x,\varepsilon)=\dfrac{x^{3}}{3}-x+\varepsilon\dfrac{x}{1-x^{2}}+\varepsilon^{2}\dfrac{x(1+x^{2})}{(1-x^{2})^{4}}+O\left({\varepsilon^{3}}\right).$
(56)

By using the Flow Curvature Method the slow manifold equation associated with
this system reads:

$\phi({\vec{X}})=\phi\left({x,y,\varepsilon}\right)=9y^{2}+\left({9x+3x^{3}}\right)y+6x^{4}-2x^{6}+9x^{2}\varepsilon=0.$
(57)

So, by plugging the perturbation expansion (56) into (57) and then, solving
order by order leads to:

$y=\frac{x^{3}}{3}-x+\varepsilon\frac{x}{1-x^{2}}+\varepsilon^{2}\frac{x}{\left({1-x^{2}}\right)^{3}}+O\left({\varepsilon^{3}}\right).$
(58)

|  |
---|---|---
(a) |  | (b)
|  |
(c) |  | (d)
Figure 1: Slow manifold of Van der Pol system (9).

So, both slow invariant manifolds approximated by Eqs. (56) & (58) are
completely identical up to order one in $\varepsilon$. At order
$\varepsilon^{2}$ a difference appears which is due to the fact that the slow
invariant manifold obtained with the Flow Curvature Method is defined, for a
two-dimensional dynamical system, by the second order tensor of curvature,
i.e. by a determinant (36) involving the first and second time derivatives of
$\vec{X}$. If one makes the same computation as previously but with the Lie
derivative of this determinant we obtain a determinant containing the first
and third time derivatives of $\vec{X}$, i.e. the third order tensor
$\det(\dot{\vec{X}},\dddot{\vec{X}})=0$. Then there is no more difference
between order two in $\varepsilon$ and slow invariant manifold given by both
methods are exactly the same. These results, highlighted in figures 1 and more
particularly in Fig. 1d, had been already found by Bruno Rossetto Rossetto
(1986, 1987) by using his Successive Approximations Method and then, by Gear
et al. Gear et al. (2005) and Zagaris et al. Zagaris et al. (2009) with their
Zero-Derivative Principle. In figures 1, we have plotted in blue the slow
invariant manifold approximated with the Flow Curvature Method (see Fig. 1a)
and in cyan dotted line the slow invariant manifold approximated with the
Geometric Singular Perturbation Method (see Fig. 1a & 1c). Both slow invariant
manifolds are compared in figure 1b in which we emphasize that Flow Curvature
Method only provides an approximation up to order one in $\varepsilon$. In
figures 1c & 1d, we have plotted the time derivative of the slow invariant
manifold obtained with the Flow Curvature Method, i.e., given by the third
order tensor $\det(\dot{\vec{X}},\dddot{\vec{X}})=0$ in magenta dotted line.
Then, as observed on Fig. 1d, there is no difference between both slow
invariant manifolds. According to Benoît et al. Benoît et al. (2015):

> “However, the global solution sets are very complicated and have branches
> which do not approximate any slow manifold. These branches are the ghosts.
> Close to the fold points of the critical manifolds the branches which
> approximate a Fenichel manifold merge with a ghost in a fold. Furthermore,
> we see that the solution to $G=0$ [$\phi({\vec{X}})=0$, according to the
> Flow Curvature Method] has a double point at the equilibrium where a ghost
> intersects the approximation to the Fenichel manifold transversally.”

As highlighted in Fig. 1a & 1c, the slow invariant manifold approximated with
the Flow Curvature Method obviously contains some “ghost parts” near the fold.
Nevertheless, let’s notice that such “ghost parts” do also exist for the
Geometric Singular Perturbation Theory as emphasized in Fig. 1a & 1c by the
dark points representing various initial conditions. Moreover, the “double
point at the equilibrium where a ghost intersects the approximation to the
Fenichel manifold transversally”, i.e. the curve with the shape of a
lemniscate of Bernoulli, provides in fact the eigendirections of both slow and
fast eigenvectors of Van der Pol system at the equilibrium point, i.e. the
origin.

###  IV.2 Lorenz system

In the beginning of the sixties a young meteorologist of the M.I.T., Edward N.
Lorenz, was working on weather prediction. He elaborated a model derived from
the Navier-Stokes equations with the Boussinesq approximation, which described
the atmospheric convection. Although his model lost the correspondence to the
actual atmosphere in the process of approximation, chaos appeared from the
equation describing the dynamics of the nature. Let’s consider the Lorenz
model Lorenz (1963):

$\left\\{\begin{aligned} \dot{x}&=\sigma\left(y-x\right),\hfill\\\
\dot{y}&=-xz+rx-y,\hfill\\\ \dot{z}&=xy-bz,\hfill\\\ \end{aligned}\right.$
(59)

with the following parameters set $\sigma=10$, $b=8/3$ and by setting $r=28$
this dynamical system exhibits the famous Lorenz butterfly. Although the
Lorenz’s system (59) has no singular approximation or no critical manifold, it
has been numerically stated by Rossetto et al. Rossetto et al. (1998) that its
Jacobian matrix possesses at least a large and negative real eigenvalue in a
large domain of the phase space. So, it is considered as a slow fast dynamical
system but not as a singularly perturbed dynamical system. Thus, none of the
singular perturbation-based methods can provide the slow invariant manifold
associated with Lorenz system. However, it can be obtained with any curvature-
based methods (see Rossetto et al. Rossetto et al. (1998)) and more
particularly with the Flow Curvature Method (see Ginoux et al. Ginoux &
Rossetto (2006a) and Ginoux Ginoux (2009)) which gives the following equation:

$\displaystyle\phi(x,y,z)=-12954816x^{4}+21168x^{6}+13331304x^{3}y-2772x^{5}y\hfill$
$\displaystyle-15012x^{4}y^{2}+27x^{6}y^{2}-6554520xy^{3}+13410x^{3}y^{3}+340200y^{4}\hfill$
$\displaystyle+15120x^{2}y^{4}-8100xy^{5}+15906240x^{2}z+1311744x^{4}z-1512x^{6}z\hfill$
$\displaystyle+5112720xyz-1454430x^{3}yz-45x^{5}yz-5680800y^{2}z-720x^{2}z^{4}\hfill$
$\displaystyle-790440x^{2}y^{2}z+540x^{4}y^{2}z+456180xy^{3}z-810x^{3}y^{3}z+1800y^{4}z\hfill$
$\displaystyle-1686544x^{2}z^{2}-45252x^{4}z^{2}+27x^{6}z^{2}-317920xyz^{2}+58320x^{3}yz^{2}\hfill$
$\displaystyle+372800y^{2}z^{2}+15750x^{2}y^{2}z^{2}-8100xy^{3}z^{2}+59040x^{2}z^{3}+540x^{4}z^{3}\hfill$
$\displaystyle+9297666x^{2}y^{2}+6480xyz^{3}-810x^{3}yz^{3}-7200y^{2}z^{3}=0.\hfill$
(60)

In fact, it can been numerically stated that in the vicinity of this slow
manifold, the Jacobian matrix is quasi-stationary and so $dJ/dt=0$. So, the
invariance of this slow manifold (IV.2) is stated according to Proposition 5.
Let’s notice that such slow invariant manifold exhibits the symmetry of
Lorenz’s system (59), i.e. $\phi(-x,-y,z)=\phi(x,y,z)$.

Figure 2: Slow invariant manifold of the Lorenz system (59).

On figure 2, we observe some kind of “horns” near the two fixed points. Such
“horns” are the three-dimensional analog to the “ghost parts” in dimension
two. As previously, they provide the eigendirection of both slow and fast
eigenvectors of Lorenz system at the fixed points. For the origin, their
shapes correspond to the nature of this equilibrium point, i.e., a saddle-
node. For the two others, these “horns” represent the cone of attraction of
the eigendirection of the negative real eigenvalue.

Let’s notice that Lorenz system is paradigmatic since it also implies chaos in
lasers due to the fact that a class C laser is equivalent to the Lorenz model.
So, the determination of the slow invariant manifold equation of such system
is of great importance. Beyond this, class B lasers, ruled only by two
equations, are slow-fast dynamical systems and the addition of a third
equation or a time modulation on a system parameter leads to chaos Meucci et
al. (2021).

##  V Discussion

During the last century, the analysis of a great number of phenomena modeled
with dynamical systems, i.e., with sets of nonlinear ordinary differential
equations, highlighted the existence of at least two time scales for their
evolution: a slow time and a fast time. This was transcribed by the presence
of at least a small multiplicative parameter $\varepsilon$ in the velocity
vector field of these dynamical systems. They were thus called singularly
perturbed dynamical systems and it was proved that they possess slow invariant
manifolds. Then, various methods were developed to compute such slow invariant
manifolds or, at least an asymptotic expansion in power of $\varepsilon$.
Thus, determination of the slow invariant manifold equation turned into a
regular perturbation problem and the seminal works of Wasow Wasow (1965), Cole
Cole (1968), O’Malley O’Malley (1974, 1991) and Fenichel Fenichel (1971, 1974,
1977, 1979) gave rise in the 1960s-1970s to the so-called Geometric Singular
Perturbation Theory. At the end of the 1980s, Rossetto Rossetto (1986, 1987)
developed the Successive Approximations Method to approximate the slow
invariant manifold of singularly perturbed dynamical systems. In 2005, Gear et
al. Gear et al. (2005) and then, Zagaris et al. Zagaris et al. (2009) used the
Zero-Derivative Principle for the same purpose. We have established in this
work that Geometric Singular Perturbation Theory, Successive Approximations
Method and Zero-Derivative Principle are absolutely identical, i.e., provide
exactly the same approximation of slow invariant manifold equation and so
belong to the first category we have called: Singular Perturbation-Based
Methods. However, according to O’Malley O’Malley (1974), Rossetto Rossetto
(1987) and Benoît et al. Benoît et al. (2015), the main drawback of these
methods is that the validity of the asymptotic expansion in power of
$\varepsilon$, approximating the slow invariant manifold equation, is expected
to breakdown near the fold or, near non-hyperbolic regions.

Beside, singularly perturbed dynamical systems there are dynamical systems
without an explicit timescale splitting. Some of these systems, which have
been called slow-fast dynamical systems, have the following property: their
Jacobian matrix has at least a fast eigenvalue, i.e. with the largest absolute
value of the real part. For such systems, Singular Perturbation-Based Methods
can not be applied anymore. Thus, in the beginning of the 1990s various
approaches have been proposed in order to approximate slow manifolds of such
slow-fast dynamical systems. In 1992, Maas and Pope introduced the Intrinsic
Low-Dimensional Manifold Method Maas & Pope (1992) and two year later, Brøns
and Bar-Eli Brøns & Bar-Eli (1994), their Inflection Line Method only
applicable to two-dimensional slow fast dynamical systems. In 1998, Rossetto
et al. Rossetto et al. (1998) proposed the Tangent Linear System Approximation
and then, the Flow Curvature Method was developed by Ginoux et al. Ginoux &
Rossetto (2006a); Ginoux et al. (2008) and Ginoux Ginoux (2009). We have
proved in this work that the Flow Curvature Method encompasses the Intrinsic
Low-Dimensional Manifold Method, the Tangent Linear System Approximation
Method and the Inflection Line Method which all belong to the second category
we have called Curvature-Based Methods. At last, we have also established the
identity between the slow invariant manifold equation given by the Geometric
Singular Perturbation Method and by the Flow Curvature Method up to suitable
order in $\varepsilon$. Then, one of the main criticism made against the Flow
curvature Method, i.e., the existence and significance of “ghost parts” in the
slow invariant manifold has been completely clarified and moreover, it has
been proved that such “ghost parts” do also exist for the Geometric Singular
Perturbation Theory. The existence of such “ghost parts” is due to the loss of
normal hyperbolicity near the fold or, near non-hyperbolic regions.

According to Heiter et al. Heiter & Lebiedz (2018), the Flow Curvature Method
“uses extrinsic curvature of curves in hyperplanes (codimension-1 manifolds)
in order to derive a determinant criterion for computing slow manifold
points.” Last year, Poppe et al. Poppe & Lebiedz (2019) generalized the Flow
Curvature Method based upon intrinsic curvature and reformulated it in the
framework of Riemannian Geometry. With these authors, we “share the opinion
that the field of differential geometry is an appropriate frame to gain
further insight in order to adequately define SIMs as slow attracting phase
space structures exploited for model reduction purposes.”

## Acknowledgments

Author would like to thank Pr. Alain Goriely who convinced him of the
necessity to write this paper and his wife Dr. Roomila Naeck who contributed
to state the proof presented in Sec. 4.

## References

  * Andronov & Chaikin (1937) Andronov, A. A. & Chaikin, S. E., Theory of Oscillators, I. Moscow, 1937; English transl., Princeton Univ. Press, Princeton, N.J., 1949.
  * Benoît et al. (2015) Benoît, E., Brøns, M., Desroches, M. & Krupa, M., “Extending the zero-derivative principle for slow-fast dynamical systems,” Z. Angew. Math. Phys. 66, (2015) pp. 2255–2270.
  * Brøns & Bar-Eli (1994) Brøns, M. & Bar-Eli, K., “Asymptotic analysis of canards in the EOE equations and the role of the inflection line,” Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 445, (1994) pp. 305–322.
  * Coddington & Levinson (1955) Coddington, E. A. & Levinson, N., Theory of Ordinary Differential Equations, Mac Graw Hill, New York, 1955.
  * Cole (1968) Cole, J. D., Perturbation Methods in Applied Mathematics, Blaisdell, Waltham, MA, 1968.
  * Darboux (1878) Darboux, G., “Sur les éuations différentielles algébriques du premier ordre et du premier degré” Bull. Sci. Math. 2(2), (1878) pp. 60–96, 123–143 & 151–200.
  * Fenichel (1971) Fenichel, N., “Persistence and smoothness of invariant manifolds for flows,” Ind. Univ. Math. J. 21, (1971) pp. 193–225.
  * Fenichel (1974) Fenichel, N., “Asymptotic stability with rate conditions,” Ind. Univ. Math. J. 23, (1974) pp. 1109–1137.
  * Fenichel (1977) Fenichel, N., “Asymptotic stability with rate conditions II,” Ind. Univ. Math. J. 26, (1977) pp. 81–93.
  * Fenichel (1979) Fenichel, N., “Geometric singular perturbation theory for ordinary differential equations,” J. Diff. Eq. 31, (1979) pp. 53–98.
  * Gear et al. (2005) Gear, C. W., Kaper, T. J., Kevrekidis, I. G. & Zagaris, A., “Projecting to a slow manifold: singularly perturbed systems and legacy codes,” SIAM J. Applied Dynamical Systems, Mathematics 4(3), (2005) pp. 711–732.
  * Gilmore et al. (2010) Gilmore, R., Ginoux, J. M., Jones, T., Letellier, C. & Freitas, U. S., “Connecting curves for dynamical systems,” J. Phys. A: Math. Theor. 43 (2010) 255101 (13pp).
  * Ginoux & Rossetto (2006a) Ginoux, J. M. & Rossetto, B., “Differential geometry and mechanics applications to chaotic dynamical systems,” International Journal of Bifurcation & Chaos 4(16), (2006) pp. 887–910.
  * Ginoux & Rossetto (2006b) Ginoux, J. M. & Rossetto, B., “Slow manifold of a neuronal bursting model,” In Emergent Properties in Natural and Articial Dynamical Systems (eds M.A. Aziz-Alaoui & C. Bertelle), Springer-Verlag, (2006) pp. 119–128.
  * Ginoux et al. (2008) Ginoux, J. M., Rossetto, B. & Chua, L. O., “Slow invariant manifolds as curvature of the flow of dynamical systems,” International Journal of Bifurcation & Chaos 11(18), (2008) pp. 3409–3430.
  * Ginoux (2009) Ginoux, J. M., Differential Geometry applied to Dynamical Systems, World Scientific Series on Nonlinear Science, Series A 66 (World Scientific, Singapore) 2009.
  * Ginoux (2014) Ginoux, J. M., “The Slow invariant manifold of the Lorenz-Krishnamurthy model,” Qualitative Theory of Dynamical Systems 13(1), (2014) pp. 19–37.
  * Ginoux (2017) Ginoux, J. M., History of Nonlinear Oscillations Theory, Archimede, New Studies in the History and Philosophy of Science and Technology, vol. 49, Springer, New York, 2017.
  * Gluck (1966) Gluck, H., “Higher curvatures of curves in Euclidean space,” The American Mathematical Monthly, 73(7), (1966) pp. 699–704.
  * Guckhenheimer & Holmes (1983) Guckhenheimer, J. & Holmes, Ph., Nonlinear Oscillations, Dynamical Systems and Bifurcation of Vector Fields, Springer-Verlag, New York, 1983.
  * Guckenheimer et al. (2000) Guckenheimer, J., Hoffman, K. & Weckesser, W., “Numerical Computation of Canards,” International Journal of Bifurcation & Chaos 10(12), (2000) pp. 2669–2687.
  * Heiter & Lebiedz (2018) Heiter, P. & Lebiedz, D., “Towards differential geometric characterization of slow invariant manifolds in extended phase space: Sectional curvature and flow invariance,” SIAM Journal on Applied Dynamical Systems 17, (2018) pp. 732–753.
  * Hirsch et al. (1977) Hirsch, M. W., Pugh, C. C. & Shub, M., Invariant Manifolds, Springer-Verlag, New York, 1977.
  * Jones (1994) Jones, C. K. R. T., Geometric Singular Pertubation Theory, In Dynamical Systems, Montecatini Terme, L. Arnold, Lecture Notes in Mathematics, vol. 1609, Springer-Verlag, New-York, (1994) pp. 44–118.
  * Kaper (1999) Kaper, T., “An introduction to geometric methods and dynamical systems theory for singular perturbation problems,” In Analyzing multiscale phenomena using singular perturbation methods, (Baltimore, MD, 1998), pp. 85–131. Amer. Math. Soc., Providence, RI, 1999.
  * Lebiedz & Unger (2016) Lebiedz, D. & Unger, J., “On unifying concepts for trajectory-based slow invariant attracting manifold computation in kinetic multiscale models,” Mathematical and Computer Modelling of Dynamical Systems 22(2), (2016) pp. 87–112.
  * Le Corbeiller (1931) Le Corbeiller, Ph., “Les Systèmes Auto-entretenus et les Oscillations de Relaxation,” Actualités scientifiques et industrielles, vol. 144 (1931) Hermann, Paris.
  * Levinson (1949) Levinson, N., “A second order differential equation with singular solutions,” Ann. Math. 50, (1949) pp. 127–153.
  * Llibre & Medrado (2007) Llibre, J. & Medrado, J. C., “On the invariant hyperplanes for d-dimensional polynomial vector fields” J. Phys. A: Math. Theor. 40, (2007) pp. 8385–8391.
  * Lorenz (1963) Lorenz, E. N., “Deterministic nonperiodic flow,” J. Atmos. Sci. 20, (1963) pp. 130-141.
  * Maas & Pope (1992) Maas, U. & Pope, S. B., “Simplifying chemical kinetics: Intrinsic low-dimensional manifolds in composition space”, Combust. Flame, 88, (1992) pp. 239-264.
  * O’Malley (1974) O’Malley, R. E., Introduction to Singular Perturbations, Academic Press, New York, 1974.
  * O’Malley (1991) O’Malley, R. E., Singular Perturbation Methods for Ordinary Differential Equations, Springer-Verlag, New York, 1991.
  * Meucci et al. (2021) Meucci R., Euzzor S., Arecchi F.T. & Ginoux J.M., “Minimal universal model for chaos in laser with feedback,” International Journal of Bifurcation & Chaos, (2021) in press, https://arxiv.org/abs/2010.06526.
  * Minorsky (1962) Minorsky, N., Nonlinear Oscillations, Princeton/Toronto/New York: D. Van Nostrand, 1962.
  * Poincaré (1881) Poincaré, H., “Sur les courbes définies par une équation différentielle,”
Journal de mathématiques pures et appliquées, Série III 7, (1881) pp. 375–422.

  * Poincaré (1882) Poincaré, H., “Sur les courbes définies par une équation différentielle,”
Journal de mathématiques pures et appliquées, Série III 8, (1882) pp. 251–296.

  * Poincaré (1885) Poincaré, H., “Sur les courbes définies par une équation différentielle,”
Journal de mathématiques pures et appliquées, Série IV 1, (1885) pp. 167–244.

  * Poincaré (1886) Poincaré, H., “Sur les courbes définies par une équation différentielle,”
Journal de mathématiques pures et appliquées, Série IV 2, (1886) pp. 151–217.

  * Poincaré (1892) Poincaré, H., Les méthodes Nouvelles de la Mécanique Céleste, Vol. I, II & III, Gauthier-Villars, Paris, 1892-93-99.
  * Poppe & Lebiedz (2019) Poppe, J. & Lebiedz, D., “Towards differential geometric characterization of slow invariant manifolds in extended phase space II: geodesic stretching,” arXiv: Dynamical Systems:1912.00676, 2019.
  * Rossetto (1986) Rossetto, B., “Trajectoires lentes des syst‘emes dynamiques lents-rapides,” In: Analysis and Optimization of System, Springer Berlin Heidelberg, (1986) pp. 680–695.
  * Rossetto (1987) Rossetto, B., “Singular approximation of chaotic slow-fast dynamical systems,” In: The Physics of Phase Space Nonlinear Dynamics and Chaos Geometric Quantization, and Wigner Function, Springer Berlin Heidelberg, (1987) pp. 12–14.
  * Rossetto et al. (1998) Rossetto, B., Lenzini, T., Ramdani, S. & Suchey, G., “Slow fast autonomous dynamical systems,” International Journal of Bifurcation & Chaos 8(11), (1998) pp. 2135–2145.
  * Schlomiuk (1993) Schlomiuk, D., “Elementary first integrals of differential equations and invariant algebraic curves,” Expositiones Mathematicae 11, (1993) pp. 433–454.
  * Tikhonov (1952) Tikhonov, A. N., “Systems of differential equations containing small parameters in the derivatives,” Mat. Sbornik N.S. 31(73)3, (1952) pp. 575–586.
  * Van der Pol (1926) Van der Pol, B., “On “relaxation-oscillations”,” The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, (VII), 2, (1926) pp. 978–992.
  * Wasow (1965) Wasow, W. R., Asymptotic Expansions for Ordinary Differential Equations, Wiley-Interscience, New York, 1965.
  * Zagaris et al. (2009) Zagaris, A., Gear, C. W., Kaper, T. J. & Kevrekidis, Y. G., “Analysis of the accuracy and convergence of equation-free projection to a slow manifold,” ESAIM: Math. Model. Num. 43, (2009) pp. 757–784.