<figure><img src="image_1.png" /></figure>

Canards, canard cascades and black swans

Vladimir Sobolev* — Elena Shchepakina**

Department of Technical Cybernetics
Samara State Aerospace University (SSAU)
Samara 443086
Russia

*hsablem@gmail.com **shchepakina@yahoo.com

**ABSTRACT.** The paper is devoted to the investigation of the slow integral manifolds of variable stability. The existence of non periodic canards, canard cascades and black swans is stated. The theoretical developments are illustrated by several examples.

**RÉSUMÉ.** L'article est consacré à l'étude de variétés lentes intégrales de stabilité variable. L'existence de canards non périodiques, de cascades de canard et des cygnes noirs est établie. Les developments théoriques sont illustrés par plusieurs exemples.

**KEYWORDS** : Invariant manifold, stability, canard, canard cascade, black swan.

**MOTS-CLÉS** : Variété invariante, stabilité, canard, cascade de canards, cygne noir.
# 1. Introduction

In this paper we outline an approach for constructing invariant surfaces of variable stability. As a one-dimensional example of such objects it is possible to consider a canard. The term “canard” (or “duck-trajectory”) originally has been introduced in the context of the van der Pol equation [1]–[4], but later have been investigated for other types of singularly perturbed systems including the ones of higher dimensions, i.e. when the phase variables are vectors (see, for example, [5]–[10] and the references therein). In the first papers devoted to canards non-standard analysis was the main tool of investigations, the Gevrey version of matched asymptotic expansions was used in [11]–[13]

From the viewpoint of the geometrical theory of singular perturbations a canard may be considered as a result of gluing stable and unstable slow invariant manifolds at one point of the breakdown surface due to the availability of an additional scalar parameter in the differential system. This approach was first proposed in [14, 15] and was then applied in [7], [16]–[20].

Other example of invariant surfaces of variable stability may be a canard cascade [21, 22] which is a result of gluing stable (attractive) and unstable (repulsive) slow invariant manifolds at several points of the breakdown surface. For the possibility of this gluing several additional parameters are required. If we take an additional function of a vector variable parameterizing the breakdown surface, we can glue the stable and unstable slow invariant manifolds at all points of the breakdown surface simultaneously. As a result we obtain the continuous stable/unstable invariant surface or black swan [8, 18, 22, 23, 24]. Such surfaces are considered as a multidimensional analogue of the notion of a canard.

From the viewpoint of analytic representation, an invariant surface of variable stability is closely connected to a solution that is bounded on the whole axis of the corresponding differential system. Furthermore, the continuity condition for the invariant surface yields the analytic representation of the gluing function.

# 2. Theoretical Background

We briefly recall the basic definitions and notions concerning the invariant surfaces of variable stability of singularly perturbed systems. The main object of our consideration is the autonomous singularly perturbed system of the form

$$ \frac{dx}{dt} = f(x, y, \alpha, \varepsilon), \qquad [1] $$

$$ \varepsilon \frac{dy}{dt} = g(x, y, \alpha, \varepsilon), \qquad [2] $$

where x and y are vectors in Euclidean spaces $R^n$ and $R^m$, respectively, ε is a small positive parameter, vector-functions f and g are sufficiently smooth and their values are
comparable to unity. The slow and fast subsystems are represented by (1) and (2), respectively.

The degenerate equation $g(x, y, \alpha, 0) = 0$ describes a slow surface $S$ (or slow curve in the case of the scalar variable $x$). The subset $\Gamma$ of $S$ given by

$$ \det \left| \frac{\partial g}{\partial y}(x, y, \alpha, 0) \right| = 0 \quad [3] $$

determines a *breakdown surface* (breakdown, turning or jump points in the scalar case) [22]. The surface $\Gamma$ divides $S$ into foliations on which

$$ \det \left| \frac{\partial g}{\partial y}(x, y, \alpha, 0) \right| \neq 0. $$

By the implicit function theorem, a leaf of the slow surface is determined by a well-defined vector-function: $y = \varphi(x, \alpha)$. A slow surface can consist of several foliations determined by different functions $y = \varphi_i(x, \alpha)$, domains of which can intersect, depending on the structure of the slow surface. It should be noted that the conditions of the implicit function theorem do not hold on the leaf boundary given by (3).

The leaf $S^s$ of the slow surface is *stable* (or *attractive*) if the spectrum of the Jacobian matrix

$$ \frac{\partial g}{\partial y}(x, \varphi(x, \alpha), \alpha, 0) $$

is located in the left open complex half-plane. If there is at least one eigenvalue of the Jacobian matrix with a positive real part then the leaf of $S$, say $S^u$, is *unstable* (or *repulsive*).

In an $\varepsilon$-neighborhood of $S^s$ ($S^u$) there exists a stable (unstable) slow invariant manifold. The slow invariant manifold is defined as an invariant surface of slow motions.

**Definition 1** A smooth surface $S_\varepsilon$ in $R^m \times R^n$ is called an invariant manifold of the system (1), (2) if any trajectory of the system that has at least one point in common with $S_\varepsilon$ lies entirely on $S_\varepsilon$.

The only invariant manifolds of system (1), (2) discussed here are those of dimension $m$ (the dimension of the slow variable $x$) that can be represented as graphs of vector-valued functions $y = h(x, \alpha, \varepsilon) = \varphi(x, \alpha) + O(\varepsilon)$, $\varepsilon \to 0$. Here, it is assumed that $h(x, \alpha, \varepsilon)$ is a sufficiently smooth function of $\varepsilon$. Such integral manifolds are called manifolds of slow motions — the origin of this term lies in nonlinear mechanics.

The notions of (un)stability for slow invariant manifold become clear if we rescale the system (1), (2) via the change of variables $\tau = t/\varepsilon$. Hence we obtain the equivalent system

$$ \frac{dx}{d\tau} = \varepsilon f(x, y, \alpha, \varepsilon), \quad [4] $$
$$ \frac{dy}{d\tau} = g(x, y, \alpha, \varepsilon). \qquad [5] $$

For simplicity we consider the case when the variables $x$ and $y$ are scalar. The stable $S^s$ and unstable $S^u$ subsets of the slow curve $S$ for which $g_y < 0$ and $g_y > 0$, respectively, are separated by the jump point(s) at which $g_y = 0$. When $\varepsilon = 0$, $S$ consists of equilibrium points of the system (4), (5) which are stable if $g_y < 0$ and unstable if $g_y > 0$. When $\varepsilon > 0$ and sufficiently small, the manifold of equilibria is replaced by a boundary layer close to $S$ with slow dynamics. There is an attractive trajectory in the boundary layer close to the stable branch $S^s$, the *stable invariant manifold* $S_\varepsilon^s$, and there is a repulsive trajectory close to an unstable branch $S^u$, the *unstable invariant manifold* $S_\varepsilon^u$ [25].

The presence of the additional scalar parameter $\alpha$ provides the possibility of gluing the stable and unstable invariant manifolds at the jump point to form a single trajectory, the canard. The term “canard” (or duck-trajectory) had been originally given by French mathematicians to the intermediate periodic trajectories of the van der Pol equation between the small and the large orbits due to their special shapes. However, in our work a canard is a one-dimensional slow invariant manifold of variable stability.

**Definition 2** *Trajectories which at first move along the stable slow invariant manifold and then continue for a while along the unstable slow invariant manifold are called canards or duck-trajectories.*

**Definition 3** *Trajectories which at first move along the unstable slow invariant manifold and then continue for a while along the stable slow invariant manifold are called false canard trajectories.*

Therefore it is possible to consider a straight line as a basic model of a canard.

**Example 1** As the simplest system with a canard we propose

$$ \dot{x} = 1, \quad \varepsilon\dot{y} = xy. \qquad [6] $$

Here $g(x, y) = xy$. It is clear that the trajectory $y = 0$ is a canard. The left part ($x < 0$) is attractive since $\partial g/\partial y < 0$ and the right part ($x > 0$) is repulsive. These two parts are divided by a turning point, which separates stable and unstable parts of the slow curve, at $x = 0$.

Note that $y = 0$ is global (the longest) canard and the unique solution that is bounded on the whole axis of system (6). Others solutions of system (6) which are local canards, are infinitely increasing (decreasing) as $t \to \infty$ or $t \to -\infty$.

**Example 2** Some generalization of (6) is the system

$$ \dot{x} = 1, \quad \varepsilon\dot{y} = xy + \phi(x, \varepsilon) + \alpha, $$
the unique solution that bounded on the whole axis of which has the representation

$$y = \begin{cases} \frac{1}{\varepsilon} \int_{-\infty}^{x} \exp\left(\frac{x^2 - s^2}{2\varepsilon}\right) [\phi(s, \varepsilon) + \alpha] ds, & x \le 0, \\ -\frac{1}{\varepsilon} \int_{x}^{\infty} \exp\left(\frac{x^2 - s^2}{2\varepsilon}\right) [\phi(s, \varepsilon) + \alpha] ds, & x > 0. \end{cases} \quad [7]$$

The “canard value” of the parameter $\alpha$ can be obtained from the continuity condition for (7) at $x = 0$:

$$\varepsilon^{-1} \int_{-\infty}^{\infty} e^{-s^2/2\varepsilon} [\phi(s, \varepsilon) + \alpha] ds = 0$$

or

$$\alpha_{\phi} = -\frac{1}{\sqrt{2\pi\varepsilon}} \int_{-\infty}^{\infty} e^{-s^2/2\varepsilon} \phi(s, \varepsilon) ds.$$

This technique of canard construction is also applicable to further generalizations.

### 3. Canard as a unique bounded on the whole axis solution

Will apply the technique of canard construction described in the previous section for the planar system

$$\dot{x} = 1, \quad \varepsilon\dot{y} = 2xy + F(x, y, \alpha, \varepsilon), \quad [8]$$

where

$$F(x, y, \alpha, \varepsilon) = f(x, y) + \varepsilon g(x, y, \varepsilon) + \varepsilon \alpha (1 + h(x, y, \varepsilon)).$$

We assume that

1) The functions $f$, $g$ and $h$ are sufficiently smooth in the domain

$$\Omega = \{x \in \mathbb{R}, y \in \mathbb{R}, \alpha \in \mathbb{R}, \varepsilon \in [0, \varepsilon_0]\}.$$

2) The functions $f$, $g$ and $h$ satisfy the following inequalities in the domain $\Omega$

$$|f(x, y)| \le C|y|^2, \quad [9]$$

$$|g(x, y, \varepsilon)| \le C\mu, \quad [10]$$

$$|h(x, y, \varepsilon)| \le C\mu, \quad [11]$$

$$|f(x, y) - f(x, \bar{y})| \le \frac{L}{2} (|y| + |\bar{y}|) |y - \bar{y}|, \quad [12]$$
$$|g(x, y, \varepsilon) - g(x, \bar{y}, \varepsilon)| \leq L|y - \bar{y}|, \qquad [13]$$

$$|h(x, y, \varepsilon) - h(x, \bar{y}, \varepsilon)| \leq L|y - \bar{y}|, \qquad [14]$$

where C and L are the positive constants, μ is a sufficiently small positive constant.

Let us reduce the system (8) to the form

$$\varepsilon \frac{dy}{dx} = 2xy + F(x, y, \alpha, \varepsilon). \qquad [15]$$

Due to (9) the straight line $y = 0$ is the slow curve of (15), the part of which with $x < 0$ ($x > 0$) is stable (unstable). Our goal is to prove that for some value of the parameter $\alpha$ there exists a one-dimensional invariant manifold of variable stability (a canard) located in a $\varepsilon$-neighborhood of the slow curve $y = 0$.

We introduce the metric space $\Phi$. Its elements are the scalar functions $\phi(x, \varepsilon)$, bounded by the value $\sqrt{\varepsilon q\mu}$, i.e.

$$|\phi(x, \varepsilon)| \leq \sqrt{\varepsilon q\mu}, \qquad [16]$$

and continuous in the domain

$$\Omega_1 = \{x \in \mathbb{R}, \mu \in (0, \mu_0), \varepsilon \in (0, \varepsilon_0)\}.$$

$\Phi$ is a complete space with the metric

$$d(\phi, \bar{\phi}) = \sup_{\Omega_1} |\phi(x, \varepsilon) - \bar{\phi}(x, \varepsilon)|, \forall \phi \in \Phi, \forall \bar{\phi} \in \Phi.$$

On the space $\Phi$ we define an integral operator $T$ by the formula

$$T(\phi)(x) = \begin{cases} \frac{1}{\varepsilon} \int_{-\infty}^{x} \exp\left(\frac{x^2 - s^2}{\varepsilon}\right) F(s, \phi(s, \varepsilon), \alpha_{\phi}, \varepsilon) ds, & x \leq 0, \\ -\frac{1}{\varepsilon} \int_{x}^{\infty} \exp\left(\frac{x^2 - s^2}{\varepsilon}\right) F(s, \phi(s, \varepsilon), \alpha_{\phi}, \varepsilon) ds, & x > 0. \end{cases} \qquad [17]$$

It should be noted that we use a modification of the usual technique of the integral manifold theory [26], and the operator in the upper (lower) line in (17) is usually used for proving the existence of stable (unstable) invariant surfaces.

The continuity condition for $T(\phi)(x)$ at $x = 0$ defines the parameter value $\alpha_{\phi}$ corresponding to the function $\phi(x, \varepsilon)$ and implies the equation

$$\varepsilon^{-1} \int_{-\infty}^{\infty} e^{-s^2/\varepsilon} F(s, \phi(s, \varepsilon), \alpha_{\phi}, \varepsilon) ds = 0$$
or

$$\alpha_{\phi} = -\frac{I_1}{I_0},$$

where

$$I_1 = \frac{1}{\varepsilon\sqrt{\varepsilon}} \int_{-\infty}^{\infty} e^{-s^2/\varepsilon} [f(s, \phi(s, \varepsilon)) + \varepsilon g(s, \phi(s, \varepsilon), \varepsilon)] ds,$$

$$I_0 = \frac{1}{\sqrt{\varepsilon}} \int_{-\infty}^{\infty} e^{-s^2/\varepsilon} [1 + h(s, \phi(s, \varepsilon), \varepsilon)] ds.$$

From (9), (10), (16) we obtain

$$|I_1| \le \frac{1}{\varepsilon\sqrt{\varepsilon}} \int_{-\infty}^{\infty} e^{-s^2/\varepsilon} |C\varepsilon\mu^2q^2 + \varepsilon C\mu| ds = \sqrt{\pi}C(\mu^2q^2 + \mu),$$

and, taking into account (11), we get

$$|I_0| \ge \frac{1}{\sqrt{\varepsilon}} \int_{-\infty}^{\infty} e^{-s^2/\varepsilon} (1 - C\mu) ds = \sqrt{\pi}(1 - C\mu).$$

Last two inequalities yield the estimate

$$|\alpha_{\phi}| \le \frac{C\mu(\mu^2q^2 + 1)}{1 - C\mu}. \qquad [18]$$

It is clear that the operator $T(\phi)(x)$ satisfies the equation

$$\varepsilon \frac{dT(\phi)}{dx} = 2xT(\phi) + F(x, \phi(x, \varepsilon), \alpha_{\phi}, \varepsilon).$$

Hence, the existence problem of the canard reduces to the existence problem of the solution of the operator equation $T(\phi) = \phi$.

In accordance with (16) to check that the operator $T(\phi)$ transforms the complete metric space $\Phi$ into itself it is sufficiently to verify that

$$|T(\phi)(x)| \le \sqrt{\varepsilon q\mu} \qquad [19]$$

for all $\phi \in \Phi$. Applying (9)-(11) and (18) and taking into account the estimate

$$\int_{x}^{\infty} \exp\left(\frac{x^2 - s^2}{\varepsilon}\right) ds < \frac{\sqrt{\pi}}{2} \sqrt{\varepsilon}, \quad x > 0,$$
we get

$$|T(\phi)(x)| \le \sqrt{\varepsilon} \frac{\sqrt{\pi}}{2} C \mu \left( \mu q^2 + 1 + \frac{1 + \mu q^2}{1 - C \mu} \right).$$

The similar inequality is obtained for $x < 0$.

The numbers $\mu$ and $q$ will be chosen to obtain the following inequality

$$\frac{\sqrt{\pi}}{2} C \left( \mu q^2 + 1 + \frac{1 + \mu q^2}{1 - C \mu} \right) \le q, \quad [20]$$

which means that (19) is fulfilled.

For any $\phi, \bar{\phi} \in \Phi$ from (9)–(14) and (18) we obtain

$$d(T(\phi), T(\bar{\phi})) \le \frac{\sqrt{\pi}}{2} \left[ L\mu q + L\sqrt{\varepsilon} + \frac{C\mu L\sqrt{\varepsilon}(1 + \mu q^2)}{1 - C\mu} \right] d(\phi, \bar{\phi}).$$

The numbers $\varepsilon_0, \mu$ and $q$ will be chosen to obtain both (20) and the inequality

$$\frac{\sqrt{\pi}}{2} \left[ L\mu q + L\sqrt{\varepsilon} + \frac{C\mu L\sqrt{\varepsilon}(1 + \mu q^2)}{1 - C\mu} \right] < 1$$

for sufficiently small $\varepsilon$ ($0 < \varepsilon \le \varepsilon_0$). In this case the contracting operator $T(p)$ transforms the complete metric space $\Phi$ into itself. Hence, due the Banach Contraction Principle the operator $T$ has a unique fixed point $\phi(x, \varepsilon) \in \Phi$ corresponding to the canard of (15).

The above argument permits us to formulate the following statement.

**Theorem 1** Let the conditions (i),(ii) hold. Then there are numbers $\varepsilon_0 > 0$, $q$ and $\mu$ such that, for all $\varepsilon \in (0, \varepsilon_0)$, there exist $\alpha = \alpha_\phi$ and $\phi(x, \varepsilon) \in \Phi$ such that $y = \phi(x, \varepsilon)$ is a canard of (8).

**Remark.** Usually the conditions (9)–(14) are fulfilled for $|x| \le r_1$, $\|y\| \le r_2$ only. In this case we deal with a local canard.

**Example 3** Originally the existence of canards for the van der Pol system

$$\frac{dx}{dt} = y - \alpha, \quad [21]$$

$$\varepsilon \frac{dy}{dt} = y - y^3/3 - x \quad [22]$$

has been proven by the non-standard analysis methods [1, 2]. We show that this statement also follows from the Theorem 1.

The change of variables

$$x = x_1 - \frac{2}{3}, \quad y = y_1 - 1, \quad \alpha = -1 + \varepsilon a,$$
transforms (21), (22) to the form

$$ 
\begin{aligned}
\dot{x}_1 &= y_1 - \varepsilon a, \\
\varepsilon \dot{y}_1 &= -x_1 + y_1^2 - \frac{1}{3} y_1^3.
\end{aligned}
$$

Applying a rectifying transformation

$$ \varepsilon z = -x_1 + y_1^2 - \frac{1}{3}y_1^3 + \frac{\varepsilon}{y_1 - 2} $$

and eliminating the variable $t$ we obtain

$$ \varepsilon \frac{dz}{dy_1} = \frac{y_1(y_1 - 2)^2 z}{1 - z(y_1 - 2)} - \frac{\varepsilon}{(y_1 - 2)^2} - \frac{\varepsilon a(y_1 - 2)}{1 - z(y_1 - 2)}. $$

By the next change of the parameter $a = 1/8 + a_1$ we get

$$ \begin{aligned} \varepsilon \frac{dz}{dy_1} = & y_1(y_1-2)^2 z \left[ 1 + z(y_1-2) + \frac{z^2(y_1-2)^2}{1-z(y_1-2)} \right] \\ & - \varepsilon \left[ \frac{1}{(y_1-2)^2} + \frac{y_1-2}{8-8z(y_1-2)} \right] - \varepsilon a_1 \frac{y_1-2}{1-z(y_1-2)}, \end{aligned} $$

or, finally,

$$ \begin{aligned} \varepsilon \frac{dz}{dy_1} = & 4y_1z + 4y_1^2 \left( \frac{y_1}{4} - 1 \right) z + y_1(y_1-2)^3 \frac{z^2}{1-z(y_1-2)} \\ & - \varepsilon \left[ \frac{1}{(y_1-2)^2} + \frac{y_1-2}{8-8z(y_1-2)} \right] - \varepsilon a_1 \frac{y_1-2}{1-z(y_1-2)}. \end{aligned} \quad [23] $$

It easy to see that (23) has the form of equation (8) and the conditions of the Theorem 1 are fulfilled in some neighborhood of $y_1 = 0$, outside of which the functions $f, g, h$ should be extended with preservation of the properties (9)–(14).

The existence theorem for the canards can be extended to the case of vector variables $x$ and $y$ [7] and can be obtained as a corollary of the existence theorem for the black swans, see Subsection 6.

## 4. Further generalizations

The described method for constructing of the invariant surfaces of variable stability in the form of bounded on the whole axis solutions of corresponding differential systems assumes various generalizations, some of which are discussed in this section. A distinctive feature in these cases lies in replacing a single gluing parameter with a gluing function to ensure continuity of the invariant surface.
## 4.1. Canards of Three-Dimensional Systems

In this subsection we discuss the existence of canards for some special types of three-dimensional systems.

We consider the autonomous system of three ordinary differential equations:

$$ \dot{x} = f(x, y, z, \varepsilon), \qquad [24] $$

$$ \dot{y} = g(x, y, z, \mu, \varepsilon), \qquad [25] $$

$$ \varepsilon \dot{z} = p(x, y, z, \mu, \varepsilon), \qquad [26] $$

where the dot denotes the first derivative with respect to time, f, g, p are scalar functions, and μ and ε are scalars. The set of points

$$ S = \{(x, y, z) : p(x, y, z, \mu, 0) = 0\} $$

is the slow surface of the system (24)–(26). The stable ($p_z < 0$) and the unstable ($p_z > 0$) parts of the slow surface are separated by a breakdown curve (i.e. the subset of the slow surface where $p_z = 0$).

The question we address is whether the system (24)–(26) has a canard. To answer this question, we investigate the two-dimensional system obtained from (24)–(26) by eliminating the variable t. It is assumed that this two-dimensional system can be represented as:

$$ y' = g/f = Y(x, y, z, \mu, \varepsilon), \qquad [27] $$

$$ \varepsilon z' = p/f = 2xz + Z(x, y, z, \mu, \varepsilon), \qquad [28] $$

where μ is a scalar and the function $Z(x, y, z, \mu, \varepsilon)$ has the following form:

$$ Z(x, y, z, \mu, \varepsilon) = Z_{1}(x, y, z) + \varepsilon(C + \mu C_{0}) + \varepsilon Z_{2}(x, y, z, \mu, \varepsilon), \qquad [29] $$

and prime represents a derivative with respect to x. Here $C, C_0$ are constants, functions $Y(x, y, z, \mu, \varepsilon)$, $Z_1(x, y, z)$ and $Z_2(x, y, z, \mu, \varepsilon)$ are defined, bounded and continuous in

$$ \Omega = \{x \in \mathbb{R}, y \in \mathbb{R}, |\mu + CC_{0}^{-1}| \le v, \varepsilon \in [0, \varepsilon_{0}], v > 0, \varepsilon_{0} > 0, $$

and satisfy the following conditions in $\Omega$:

$$ |Y(x, y, z, \mu, \varepsilon) - Y(x, \bar{y}, \bar{z}, \bar{\mu}, \varepsilon)| \le M (|y - \bar{y}| + |z - \bar{z}|) + v|\mu - \bar{\mu}|, \qquad [30] $$

$$ |Z_{1}(x, y, z)| \le M|z|^{2}, \qquad [31] $$

$$ |Z_{1}(x, y, z) - Z_{1}(x, \bar{y}, \bar{z})| \le M (|z| + |\bar{z}|)^{2} |y - \bar{y}| + \frac{M}{2} (|z| + |\bar{z}|) |z - \bar{z}|, \qquad [32] $$

$$ |Z_{2}(x, y, z, \mu, \varepsilon)| \le M\nu, \qquad [33] $$
$$|Z_2(x, y, z, \mu, \varepsilon) - Z_2(x, \bar{y}, \bar{z}, \bar{\mu}, \varepsilon)| \le M (|y - \bar{y}| + |z - \bar{z}|) + \nu |\mu - \bar{\mu}|, \quad [34]$$

where $M$ is a positive constant and $\nu$ is a sufficiently small positive constant. These conditions are required for the existence of a unique canard which is bounded for all $x$ and $y$.

The slow surface of system (27), (28) is defined by the equation $z = 0$, due to (29), (31) and the identity

$$\{2xz + Z(x, y, z, \mu, 0)\}_{z=0} \equiv 0.$$

We know that in an $\varepsilon$-neighborhood of stable and unstable foliations of the slow surface there are stable and unstable slow integral manifolds

$$z = h(x, y, \mu, \varepsilon).$$

The parameter $\mu$ ensures the existence of a gluing point of these integral manifolds, see for instance Example 4 at the end of this section. By fixing the jump point $(0, y^*)$, we can single out the trajectory

$$y = \phi(x, \mu) \quad (\phi(0, \mu) = y^*)$$

on the integral manifold $h(x, y, \mu, \varepsilon)$ which passes along the stable leaf to the jump point and then continues for a while along the unstable leaf. For convenience we use the same term 'trajectory' for both systems (24)–(26) and (27), (28). The following theorem holds.

If some natural conditions for $M$, $\nu$ and $\upsilon$ hold, then there is $\varepsilon_0$ such that for every $\varepsilon \in (0, \varepsilon_0)$ there exist $\mu = \mu^*(\varepsilon)$ and a canard corresponding to this parameter value which passes through the point $(0, y^*)$. The reader is referred to [7] for an exact statement and its proof.

### 4.1.1. Asymptotic Expansions for Canards

In this subsection the asymptotic expansions for the canards of the system (27), (28) are obtained.

It is assumed that functions $Y$ and $Z$ in (27), (28) have sufficient continuous and bounded partial derivatives with respect to all variables. For simplicity we exclude the $\varepsilon$-dependence of the functions $Y$ and $Z_2$. Then the canard and the parameter value $\mu^*$ (corresponding to this trajectory) allow asymptotic expansions in powers of the small parameter $\varepsilon$:

$$\mu^* = \sum_{i \ge 0} \varepsilon^i \mu_i,$$

$$y = \phi(x, \mu^*) = \sum_{i \ge 0} \varepsilon^i \phi_i(x), \quad [35]$$

$$z = \psi(x, \mu^*, \varepsilon) = h \left( x, \sum_{i \ge 0} \varepsilon^i \phi_i(x), \sum_{i \ge 0} \varepsilon^i \mu_i, \varepsilon \right) = \sum_{i \ge 0} \varepsilon^i \psi_i(x).$$

We can calculate these asymptotic expansions from (27), (28).
Note that this statements can be generalized to the cases $y \in R^n, z \in R$ and $y \in R^n, z \in R^m$.

**Example 4** As a very simple example of (24)–(26) consider the system

$$ \dot{x} = 1, \quad \dot{y} = 0, \quad \varepsilon \dot{z} = 2xz + \mu - y. $$

Since $p_z = 2x$, the slow surface $2xz + \mu - y = 0$ is divided by the breakdown curve $x = 0$ into the stable part ($x < 0$) and the unstable one ($x > 0$). If $\mu$ is a parameter then the different canards are determined by

$$ \dot{x} = 1, \quad y = y_0, \quad z = 0, $$

that is, they pass through the unique gluing point $x = 0, y = y_0, z = 0$ on the breakdown curve for $\mu = y_0$.

## 4.2. Case of a multiple turning point

In this Subsection we deal with a multiple turning point.

**Example 5** Consider a modification of the Example 1 given by

$$ \dot{x} = 1, \quad \varepsilon \dot{y} = x^k y, $$

or

$$ \varepsilon \frac{dy}{dx} = x^k y. \qquad [36] $$

For $k = 3, 5, \dots$ the straight line $y = 0$ is a canard.

**Example 6** Some generalization of (36) is the equation

$$ \varepsilon \frac{dy}{dx} = x^k y + f(x, y, \alpha). $$

In this case the gluing function can be found in the form

$$ \alpha = \alpha(x, \varepsilon) = \alpha_1(\varepsilon)x^{k-1} + \dots + \alpha_{k-1}(\varepsilon)x + \alpha_k(\varepsilon). $$

The other approach for investigation of canards in the case of a multiple turning point has been proposed in [27].
### 4.3. Model with distributed parameters

Consider nonlinear singularly perturbed parabolic system

$$ \epsilon \frac{\partial \theta}{\partial \tau} = \eta(1 - \eta) \exp \left( \frac{\theta}{1 + \beta\theta} \right) + \frac{1}{\delta} D_{\xi}\theta, \quad [37] $$

$$ \epsilon \frac{\partial \eta}{\partial \tau} = \epsilon \eta(1 - \eta) \exp \left( \frac{\theta}{1 + \beta\theta} \right) + \frac{1}{\rho} D_{\xi}\eta, \quad [38] $$

where

$$ D_{\xi}(\cdot) = \frac{\partial^2(\cdot)}{\partial\xi^2} + \frac{n}{\xi} \frac{\partial(\cdot)}{\partial\xi}, \quad n = 0, 1, 2, $$

boundary conditions

$$ \left. \frac{\partial \theta}{\partial \xi} \right|_{\xi=0} = 0, \quad \left. \theta \right|_{\xi=1} = 0, \quad \left. \frac{\partial \eta}{\partial \xi} \right|_{\xi=0} = 0, \quad \left. \frac{\partial \eta}{\partial \xi} \right|_{\xi=1} = 0 \quad [39] $$

and initial conditions

$$ \theta \Big|_{\tau=0} = 0, \quad \eta \Big|_{\tau=0} = \eta_0 / (1 + \eta_0). \quad [40] $$

This is a mathematical model of the problem of thermal explosion given heat transfer and diffusion [20]. Here $\theta$ is the dimensionless temperature; $\eta$ is the dimensionless depth of conversion; $\tau$ is the dimensionless time; $\epsilon$ and $\beta$ are small positive parameters; $\delta$ is a Frank-Kamenetsky criterion, that is the scalar parameter, characterizing initial state of the system. Depending on its value, reaction either is explosive or proceeds slowly. The value of parameter $\delta$ separating slow and explosive regimes is called critical.

Using the method of integral manifolds, the critical value of $\delta$ is calculated as an asymptotic series with respect to degrees of the small parameter $\epsilon$

$$ \delta = \delta_0(1 + \delta_1\epsilon) + O(\epsilon^2), \quad [41] $$

where critical regimes are modelled by the canards. For $n = 0$ (plane-parallel reactor), $n = 1$ (cylindrical reactor), $n = 2$ (spherical reactor) corresponding values of $\delta_0$ and $\delta_1$ are calculated: $\delta_0 = 3.513828$, 8 and 13.32; $\delta_1 = 15.58$, 1.9232 and 1.74.

## 5. Canard cascades

If it is necessary to glue stable and unstable slow invariant manifolds at several break-down points, we need several additional parameters and as a result we obtain a canard cascade.

**Definition 4** The continuous one-dimensional slow invariant manifold of (1), (2) which contains at least two canards or false canards is called a canard cascade.
In the case of a planar system, if we take an additional function whose arguments are a vector parameter and a slow variable, we can glue the stable (attractive) and unstable (repulsive) slow invariant manifolds at all breakdown points at the same time. As a result we obtain a canard cascade.

**Example 7** The differential system

$$ \dot{x} = 1, \quad \varepsilon \dot{y} = x(x-1)y $$

gives a simplest canard cascade $y = 0$ which consists of two repulsive parts ($x < 0$ and $x > 1$ since $\partial g/\partial y = x(x-1) > 0$) and one attractive part ($0 < x < 1$ since $\partial g/\partial y = x(x-1) < 0$) with two breakdown points $x = 0$ and $x = 1$.

**Example 8** An example of a system with an infinite number of breakdown points is

$$ \dot{x} = 1, \quad \varepsilon \dot{y} = y \sin x. $$

**Example 9** The slow curve of the van der Pol system (21), (22) has two jump points. To construct a canard cascade we need to have two independent parameters. We can consider $\alpha$ as a gluing function $\alpha = \alpha(x, \lambda)$ depending on the slow variable $x$ and the vector $\lambda = (\lambda_1, \lambda_2)$. The vector $\lambda$ is a function of $\varepsilon$: $\lambda = \lambda(\varepsilon)$. The consideration of a variant where the function $\alpha$ is a polynomial in $x$, i.e., $\alpha = \lambda_1 x + \lambda_2$ seems quite natural. Of course, we can use a trigonometric polynomial or the linear combination of any linearly independent functions in $x$ with the coefficients $\lambda_1, \lambda_2$. We will search the canard cascade in the form of a polynomial

$$ x = y - y^3/3 + \varepsilon q_m(y, \varepsilon), \quad [42] $$

where $q_m(y, \varepsilon)$ is an $m$th-degree polynomial in $y$. From (21), (22) and (42) we have

$$ \frac{\partial}{\partial y} \left( y - \frac{y^3}{3} + \varepsilon q_m(y, \varepsilon) \right) \left( -q_m(y, \varepsilon) \right) = y - \alpha. $$

Balancing the degrees of polynomials on both sides of this equality implies that $m = 1$ since the highest power of $y$ in $x$ is $y^3$, and therefore, $q_m(y, \varepsilon) = ay + b$. After some calculations we obtain [21, 22] the representation for gluing function

$$ \alpha = \alpha(x, \varepsilon) = 3\mu(\varepsilon)x, $$

where

$$ \mu(\varepsilon) = \frac{\sqrt{1 + 2\varepsilon} - 1}{2\varepsilon}, $$

and a canard cascade

$$ x = y - \frac{y^3}{3} + \varepsilon \mu(\varepsilon)y. $$
The canard cascade (the one-dimensional slow invariant manifold) passes near the slow curve of the system. We glue the stable and the unstable slow invariant manifolds at two jump points $A_1$ and $A_2$. In this case we have used only one additional parameter to obtain a canard cascade.

Note that the change of variable

$$x = -\varepsilon z + y - y^3/3 + \varepsilon\mu(\varepsilon)y$$

transforms (21), (22) to the system

$$\dot{y} = z - \mu y, \quad \varepsilon\dot{z} = (y - 1)(y + 1)z - 4\varepsilon\mu z$$

with the canard cascade $z = 0$ which plays the role of a bounded on the whole axes solution.

In the general case, when the slow curve of a system has $k$ jump points, it is necessary to use $k$ additional parameters to construct a canard cascade.

As was mentioned above, we consider control functions depending on the slow variable only. But the special case of the control function depending on the fast variable is of interest also. Éric Benoît is the author of the following statement.

**Remark.** The system

$$\dot{x} = \mu(y), \quad \varepsilon\dot{y} = p(y) - x$$

has a canard cascade with an invariant manifold $x = p(y) + \varepsilon q(y)$ if we choose the polynomial control $\mu(y) = -q(y)(p'(y) + \varepsilon q'(y))$. The direction of slow dynamics can be controlled by the sign of $q(y)$ to choose the sequence of true and false canards.

# 6. Black swans

We return now to the more general case and consider the multidimensional invariant surfaces of variable stability or black swans [23, 24]. These surfaces are considered as natural generalizations of the notion of a canard.

We consider a non-autonomous system obtained from an original autonomous system of type (1), (2) with a scalar fast variable by eliminating the variable $t$. It is assumed that this system can be represented as:

$$\frac{dy}{dx} = Y(x, y, z, \varepsilon), \quad y \in \mathbb{R}^n, \quad x \in \mathbb{R}; \qquad [43]$$

$$\varepsilon\frac{dz}{dx} = 2xz + a + Z(x, y, z, a, \varepsilon), \quad |z| \le r, \qquad [44]$$

where $r$ and $a_0$ are positive constants. It is supposed that the functions $Y, Z$ are continuous and satisfy the following inequalities for $x \in \mathbb{R}, y \in \mathbb{R}^n, |z| \le r, |a| \le a_0, \varepsilon \in [0, \varepsilon_0]$:

$$\|Y(x, y, z, \varepsilon)\| \le k, \quad |Z(x, y, z, a, \varepsilon)| \le M (\varepsilon^2 + \varepsilon|z| + |z|^2), \qquad [45]$$
$$
\begin{align}
& \|Y(x, y, z, \varepsilon) - Y(x, \bar{y}, \bar{z}, \varepsilon)\| \le M(\|y - \bar{y}\| + |z - \bar{z}|), \tag{46} \\
& |Z(x, y, z, a, \varepsilon) - Z(x, \bar{y}, \bar{z}, \bar{a}, \varepsilon)| \le M \{(\varepsilon + |\tilde{z}|)|z - \bar{z}| \nonumber \\
& \qquad + (\varepsilon^2 + \varepsilon|\tilde{z}| + |\tilde{z}|^2)\|y - \bar{y}\| + \varepsilon|a - \bar{a}|\}, \quad |\tilde{z}| = \max\{|z|, |\bar{z}|\}, \tag{47}
\end{align}
$$

where $\| \cdot \|$ denotes the usual norm in $R^n$ and $| \cdot |$ denotes the absolute value of a scalar, $k$ and $M$ are positive constants.

We take $a$ as a function: $a = a(y, \varepsilon)$.

Let $F$ be the complete metric space of functions $a(y, \varepsilon)$ continuous with respect to $y$
and satisfying

$$
|a(y, \varepsilon)| \leq \varepsilon^2 K, \quad |a(y, \varepsilon) - a(\bar{y}, \varepsilon)| \leq \varepsilon^2 L \|y - \bar{y}\|,
$$

for $\varepsilon \in (0, \varepsilon_0]$, where $K$ and $L$ are positive constants, with the metric defined by

$$
\rho(a, \bar{a}) = \sup_{y \in \mathbb{R}^n} |a(y, \varepsilon) - \bar{a}(y, \varepsilon)|.
$$

Let $H$ be the complete metric space of functions $h(x, y, \varepsilon)$ mapping $R \times R^n$ to $R$
continuous with respect to $x, y$ and satisfying

$$
\begin{align*}
|h(x, y, \varepsilon)| &\le \varepsilon^{\frac{3}{2}} q, \\
|h(x, y, \varepsilon) - h(x, \bar{y}, \varepsilon)| &\le \varepsilon^{\frac{3}{2}} \delta \|y - \bar{y}\|,
\end{align*}
$$

for $\varepsilon \in (0, \varepsilon_0]$, where $q$ and $\delta$ are positive constants, with the metric

$$
\rho(h, \bar{h}) = \sup_{x \in R, y \in R^n} |h(x, y, \varepsilon) - \bar{h}(x, y, \varepsilon)|.
$$

On the space *H* we define an operator *T* by the formula

$$
Th(x, y, \varepsilon) = \begin{cases} \frac{1}{\varepsilon} \int_{-\infty}^{x} \exp\left(\frac{x^2 - s^2}{\varepsilon}\right) [Z(\cdot) + a(\varphi(s, x), \varepsilon)] ds, & x \le 0, \\ -\frac{1}{\varepsilon} \int_{x}^{-\infty} \exp\left(\frac{x^2 - s^2}{\varepsilon}\right) [Z(\cdot) + a(\varphi(s, x), \varepsilon)] ds, & x > 0, \end{cases}
$$

where $Z(\cdot) = Z(s, \varphi(s, x), h(s, \varphi(s, x), \varepsilon), a(\varphi(s, x), \varepsilon), \varepsilon)$, and $\varphi(s, x)$ is defined as
follows. For any element $h \in H$ the initial value problem

$$
\begin{align*}
\frac{d\varphi}{ds} &= Y(s, \varphi, h(s, \varphi, \varepsilon), \varepsilon), \\
\varphi(x) &= y,
\end{align*}
$$
is considered. This problem is obtained from (43) if we replace $z$ by an element $h \in H$. The solution of this problem is denoted by $\Phi(s, x, y, \varepsilon|h) = \varphi(s, x)$. When the operator $T$ possesses a fixed point $h(x, y, \varepsilon)$ in $H$ then the surface $z = h(x, y, \varepsilon)$ is a slow stable/unstable invariant manifold (a black swan). It should be noted that the operator $T$ defines a manifold of the solutions which are bounded with respect to $z$.

For the detailed proof, which is multidimensional modification of the proof of the Theorem 1, the readers are referred to [23, 24]. The following theorem gives sufficient conditions for the system (43), (44) to have a black swan.

**Theorem 2** Let the conditions (45)–(47) be satisfied. Then there are numbers $\varepsilon_0 > 0$ and $K, L, q, \delta$ such that, for all $\varepsilon \in (0, \varepsilon_0)$, there exist functions $a(y, \varepsilon) \in F$ and $h(x, y, \varepsilon) \in H$ such that $z = h(x, y, \varepsilon)$ is a slow integral manifold (the black swan).

**Remark.** An original singularly perturbed differential system can be reduced to the form (43), (44) by a multidimensional analogue of the rectifying transformation described in the Example 3 (see, for instance, [7]).

**Remark.** Usually the conditions (45)–(47) are fulfilled for $|x| \le r_1$, $\|y\| \le r_2$ only. Integral manifolds in this case are local.

Note that Theorem 2 can be generalized to the case when $z$ is a vector variable.

Let us consider now the system

$$ \frac{dy}{dx} = Y(x, y, z_1, z_2, \varepsilon), \quad y \in \mathbb{R}^n, \quad x \in \mathbb{R}, \qquad [48] $$

$$ \varepsilon \frac{dz_1}{dx} = 2xz_1 + a(y, \varepsilon) + Z_1(x, y, z_1, z_2, \varepsilon), \quad z_1 \in \mathbb{R}, \qquad [49] $$

$$ \varepsilon \frac{dz_2}{dx} = A(x)z_2 + a(y, \varepsilon)B + Z_2(x, y, z_1, z_2, \varepsilon), \quad z_2 \in \mathbb{R}^m, \qquad [50] $$

where $A(x)$ is a bounded matrix, satisfying a Lipschitz condition, with eigenvalues $\lambda_i(x)$:

$$ \mathrm{Re}\lambda_i(x) \le -2\beta < 0 \quad (i = 1, 2, \dots, m), $$

$B$ is a constant vector, the continuous functions $Y, Z_1, Z_2, a$ satisfy the inequalities

$$ \|Y(x, y, z_1, z_2, \varepsilon)\| \le k, \qquad [51] $$

$$ |Z_1(x, y, z_1, z_2, \varepsilon)| \le M (\varepsilon^2 + \varepsilon\|z\| + \|z\|^2), \qquad [52] $$

$$ \|Z_2(x, y, z_1, z_2, \varepsilon)\| \le M (\varepsilon^2 + \varepsilon\|z\| + \|z\|^2), \qquad [53] $$

$$ \|Y(x, y, z_1, z_2, \varepsilon) - Y(x, \bar{y}, \bar{z}_1, \bar{z}_2, \varepsilon)\| \le M (\|y - \bar{y}\| + \|z - \bar{z}\|), \qquad [54] $$

$$ |Z_1(x, y, z_1, z_2, \varepsilon) - Z_1(x, \bar{y}, \bar{z}_1, \bar{z}_2, \varepsilon)| \le M \{(\varepsilon + \|\tilde{z}\|)\|z - \bar{z}\| + (\varepsilon^2 + \varepsilon\|\tilde{z}\| + \|\tilde{z}\|^2)\|y - \bar{y}\|\}, \qquad [55] $$
$$ \|Z_2(x, y, z_1, z_2, \varepsilon) - Z_2(x, \bar{y}, \bar{z}_1, \bar{z}_2, \varepsilon)\| \le M \{(\varepsilon + \|\tilde{z}\|)\|z - \bar{z}\| \\ + (\varepsilon^2 + \varepsilon\|\tilde{z}\| + \|\tilde{z}\|^2)\|y - \bar{y}\|\}, \quad [56] $$

where

$$ z = \begin{pmatrix} z_1 \\ z_2 \end{pmatrix}, \quad \bar{z} = \begin{pmatrix} \bar{z}_1 \\ \bar{z}_2 \end{pmatrix}, \quad \|\tilde{z}\| = \max\{\|z\|, \|\bar{z}\|\}, k > 0, M > 0. $$

We consider $a(y, \varepsilon)$ to be a continuous function satisfying the inequalities

$$ |a(y, \varepsilon)| \le \varepsilon^2 K, \quad |a(y, \varepsilon) - a(\bar{y}, \varepsilon)| \le \varepsilon^2 L \|y - \bar{y}\|, \quad [57] $$

for $\varepsilon \in (0, \varepsilon_0]$ with some positive constants $K$ and $L$.

**Theorem 3** Let the conditions (51)–(56) be satisfied. Then there are numbers $\varepsilon_0 > 0$ and $K, L, q, \delta$ such that, for all $\varepsilon \in (0, \varepsilon_0)$, there exist functions $a(y, \varepsilon)$ satisfying (57) and a continuous function $h(x, y, \varepsilon)$ satisfying

$$ |h_1(x, y, \varepsilon)| \le \varepsilon^{3/2} q, \quad |h_1(x, y, \varepsilon) - h_1(x, \bar{y}, \varepsilon)| \le \varepsilon^{3/2} \delta \|y - \bar{y}\|, \quad [58] $$

$$ \|h_2(x, y, \varepsilon)\| \le \varepsilon^2 q, \quad \|h_2(x, y, \varepsilon) - h_2(x, \bar{y}, \varepsilon)\| \le \varepsilon^2 \delta \|y - \bar{y}\|, \quad [59] $$

$$ h(x, y, \varepsilon) = \begin{pmatrix} h_1(x, y, \varepsilon) \\ h_2(x, y, \varepsilon) \end{pmatrix}, \quad h_1 \in \mathbb{R}, \quad h_2 \in \mathbb{R}^m, $$

such that $z = h(x, y, \varepsilon)$ is a black swan of the system (48)–(50).

The proof of Theorem 3 follows the lines of the proof of Theorem 2, but now $H$ is the complete metric space of functions $h(x, y, \varepsilon)$ mapping $R \times R^n$ to $R^{m+1}$, continuous with respect to $x, y$ and satisfying (58), (59) for $\varepsilon \in (0, \varepsilon_0]$. The operator $T$ is defined by the formula $T = \begin{pmatrix} T_1 \\ T_2 \end{pmatrix}$:

$$ T_1 h(x, y) = \begin{cases} -\varepsilon^{-1} \int_{-\infty}^{\infty} e^{(x^2-s^2)/\varepsilon} [Z_1(\cdot) + a(\varphi(s, x), \varepsilon)] ds, & x \ge 0, \\ \varepsilon^{-1} \int_{-\infty}^{\infty} e^{(x^2-s^2)/\varepsilon} [Z_1(\cdot) + a(\varphi(s, x), \varepsilon)] ds, & x < 0, \end{cases} $$

$$ T_2 h(x, y) = \varepsilon^{-1} \int_{-\infty}^{x} W(x, s, \varepsilon) [Z_2(\cdot) + Ba(\varphi(s, x), \varepsilon)] ds, $$

where

$$ Z_{1,2}(\cdot) = Z_{1,2}(s, \varphi(s, x), h_1(s, \varphi(s, x), \varepsilon), h_2(s, \varphi(s, x), \varepsilon), \varepsilon) $$
and $W(x, s, \varepsilon)$, $W(s, s, \varepsilon) = I$, is the fundamental matrix of the equation

$$ \varepsilon \frac{dz_2}{dx} = A(x)z_2. $$

**Remark.** The proof of the Theorem 3 can be extended to the case when the system (48)–(50) has the form

$$ \frac{dy}{dx} = Y(x, y, z_1, z_2, \varepsilon), \quad y \in \mathbb{R}^n, \quad x \in \mathbb{R}, \qquad [60] $$

$$ \varepsilon \frac{dz_1}{dx} = 2xb(x, y)z_1 + a(y, \varepsilon) + Z_1(x, y, z_1, z_2, \varepsilon), \quad z_1 \in \mathbb{R}, \qquad [61] $$

$$ \varepsilon \frac{dz_2}{dx} = A(x, y)z_2 + a(y, \varepsilon)B(x, y, \varepsilon) + Z_2(x, y, z_1, z_2, \varepsilon), \quad z_2 \in \mathbb{R}^m. \qquad [62] $$

Here the matrix $A(x, y)$ and the functions $Y, Z_1, Z_2$, and $a$ satisfy the conditions of the Theorem 3. Moreover, the matrix $A(x, y)$, the scalar function $b(x, y)$ ($0 < c_1 \le b(x, y) \le c_2$) and the vector function $B(x, y, \varepsilon)$ are continuous and bounded, and satisfy a Lipschitz condition.

**Corollary 1** The existence of canards for the systems (48)–(50) and (60)–(62) follows from the Theorem 3. Indeed, the black swan $z = h(x, y, \varepsilon)$ can be considered as a surface consisting entirely of the canards, each of which passes through a certain point of the breakdown surface. If we take the gluing function has to be a constant, $a(y, \varepsilon) = a(y_0, \varepsilon)$, the stable and unstable invariant manifolds can be glued at one point $y = y_0$ of the breakdown surface only. The canard passes only through this point.

## 7. Summary

In this paper we have looked from unified positions at the invariant surfaces of variable stability both from a geometric viewpoint, when dealing with gluing the stable and the unstable slow invariant manifolds, and in terms of their analytic representations as the unique bounded on the whole axis solutions. The invariant surfaces of variable stability are interesting not only from a mathematical viewpoint, they are successfully used for modeling of critical phenomena of diverse nature. It should be noted that the considered cases do not exhaust the possible generalizations.
## Acknowledgments

The authors would like to gratefully acknowledge Professor Eric Benoit for helpful discussions on this subject, comments and recommendations.

This work is supported in part by the Russian Foundation for Basic Research (grants 13-01-97002-p, 14-01-97018-p) and the Ministry of education and science of the Russian Federation in the framework of the implementation of Program of increasing the competitiveness of SSAU for 2013–2020 years.

## 8. References

[1] É. BENOÎT, J. L. CALLOT, F. DIENER, M. DIENER, "Chasse au canard", Collect. Math., vol. 31-32, p. 37-119, 1981-1982.

[2] É. BENOÎT, "Enlacements de canards", Comptes-Rendus de l'Académie des Sciences de Paris, serie I, 300:225-230, 1985.

[3] É. BENOÎT, editor, "Dynamic Bifurcations" (Lect. Notes in Math. vol 1493), Berlin-Heidelberg-London: Springer, 1991.

[4] É. BENOÎT, C.LOBRY, "Les canards de R³", Comptes-Rendus de l'Académie des Sciences de Paris, serie I, 294:483-488, 1982.

[5] É. BENOÎT, "Systèmes lents-rapides dans R³ et leurs canards", Société Mathématique de France, Astérisque, vol. 109-110, p. 159-191, 1983.

[6] K. TCHIZAWA, "On relative stability in 4-dimensional duck solutions", Journal of Mathematics and System Science, vol. 2 (9), p. 558-563, 2012.

[7] V. A. SOBOLEV, E. A. SHCHEPAKINA, "Duck trajectories in a problem of combustion theory", Differential Equations, vol. 32, p. 1177-1186, 1996.

[8] E. A. SHCHEPAKINA, V. A. SOBOLEV, "Black Swans and Canards in Laser and Combustion Models", in Singular Perturbation and Hysteresis (eds. M. P. Mortell et al.), Philadelphia: SIAM, 2005, pp. 207-255.

[9] F. MARINO, F. MARIN, S. BALLE, O. PIRO, "Chaotically spiking canards in an excitable system with 2D inertial fast manifolds", Phys. Rev. Lett., vol. 98, 074104, 2007.

[10] M. DESROCHES, B. KRAUSKOPF, H. M. OSINGA, "Numerical continuation of canard orbits in slow-fast dynamical systems", Nonlinearity, vol. 23, p. 739-765, 2010.

[11] É. BENOÎT, A.EL HAMIDI, A. FRUCHARD, "On combined asymptotic expansions in singular perturbations", Electronic Journal of Differential Equations, vol. 2002 (51), p. 1-27, 2002.

[12] É. BENOÎT, A.EL HAMIDI, A. FRUCHARD, "Combined asymptotic expansions", in Singular Perturbation and Hysteresis (eds. M. P. Mortell et al.), Philadelphia: SIAM, 2005, pp. 101-1109.

[13] A. FRUCHARD, R. SCHÄFKE, "Composite Asymptotic Expansions" (Lect. Notes in Math. vol 2066), Berlin-Heidelberg-London: Springer, 2013.
[14] G. N. GORELOV, V. A. SOBOLEV, “Mathematical modelling of critical phenomena in thermal explosion theory”, *Combust. Flame*, vol. 87, p. 203–210, 1991.

[15] G. N. GORELOV, V. A. SOBOLEV, “Duck-trajectories in a thermal explosion problem”, *Appl. Math. Lett.*, vol. 5, p. 3–6, 1992.

[16] V. A. SOBOLEV, E. A. SHCHEPAKINA, “Self-ignition of laden medium”, *J. Combustion, Explosion and Shock Waves*, vol. 29, p. 378–381, 1993.

[17] V. GOL'DSHTEIN, A. ZINOVIEV, V. SOBOLEV, E. SHCHEPAKINA, “Criterion for thermal explosion with reactant consumption in a dusty gas”, *Proc. London Roy. Soc. Ser. A.*, vol. 452, p. 2103–2119, 1996.

[18] E. A. SHCHEPAKINA, “Black swans and canards in self-ignition problem”, *Nonlinear Analysis: Real Word Applications*, vol. 4, p. 45–50, 2003.

[19] K. SCHNEIDER, E. SHCHEPAKINA, V. SOBOLEV, “A new type of travelling wave”, *Mathematical Methods in the Applied Sciences*, vol. 26, p. 1349–1361, 2003.

[20] G. N. GORELOV, V. A. SOBOLEV, E. SHCHEPAKINA, “Canards and critical behaviour in autocatalytic combustion models”, *Journal of Engineering Mathematics*, vol. 56, p. 143–160, 2006.

[21] V. A. SOBOLEV, “Canard Cascades”, *Discr. and Cont. Dynam. Syst. B*, vol. 18, p. 513–521, 2013.

[22] E. SHCHEPAKINA, V. SOBOLEV, M. P. MORTELL, “Singular Perturbations: Introduction to system order reduction methods with applications” (Lect. Notes in Math. vol 2114), *Berlin-Heidelberg-London: Springer*, 2014.

[23] E. SHCHEPAKINA, V. SOBOLEV, “Attracting/repelling invariant manifolds”, *Stab. Control Theory Appl.*, vol. 3, p. 263–274, 2000.

[24] E. SHCHEPAKINA, V. SOBOLEV, “Integral manifolds, canards and black swans”, *Nonlinear Analysis A*, vol. 44, p. 897–908, 2001.

[25] V. V. STRYGIN, V. A. SOBOLEV, “Effect of geometric and kinetic parameters and energy dissipation on orientation stability of satellites with double spin”, *Cosmic Research*, vol. 14, p. 331–335, 1976.

[26] YU. L. DALETSKII, M. G. KREIN, “Stability of Solutions of Differential Equations in Banach Space”, *New York: AMS*, 1974.

[27] T. FORGET, “Asymptotic study of planar canard solutions”, *Bull. Belg. Math. Soc. Simon Stevin*, vol. 15, p. 809–824, 2008.