Jean-Marc Ginoux

Jaume Llibre

Canards Existence in Memristor's Circuits

Keywords: Geometric singular perturbation theory, singularly perturbed dynamical systems, canard solutions

INTRODUCTION

As recalled by Fruchard and Schäfke [20, p. 435]: "In the late 1970s, under the leadership of George Reeb, a group of young researchers, Jean-Louis Callot, Francine and Marc Diener, Albert Troesch, Emile Urlacher and then, Eric Benoît and Imme van den Berg, based in Strasbourg some in Oran and others in Tlemcen were given as research program to develop methods for "non-standard analysis 1 " for the study of singular perturbation problems. Georges Reeb had proposed to introduce a particular control parameter a in the original van der Pol's equation [START_REF] Van Der Pol | On relaxation-oscillations[END_REF].

εẍ + (1 -x 2 ) ẋ + x = a
The study of this equation has led this group to discover surprising solutions, which they named "ducks 2 ". Van der Pol relaxation oscillator is considered as the paradigm of slow -fast systems, i.e. two-dimensional singularly perturbed system with one slow variable and one fast. It is wellknown that for the control parameter value a = 1, a Hopf bifurcation takes place in this system 3 . So, as expected by this group of researchers, by setting ε constant, one would make the amplitude of the periodic solution (limit cycle) change very fast for values of a near (just below) 1. But, the results exceeded their expectations: at one critical value a = 0.9987404512 a very small change of this parameter's value produced an amplitude drop of about 80 %.

According to Marc Diener [13, p. 38]: "It was as if the existence of medium size solutions would be a "canard 4 "! Canard is now the name of a type of solution of a slow-fast differential system, to which the above "missing" medium-size solutions belong, that had previously been ignored.". Another interpretation of the denomination "canard" also given by Diener [13, p. 45] in the the same article, is that the periodic solution resembles, for this critical parameter value, to a duck (See Fig. 1.). In the beginning of the eighties, Benoît and Lobry [START_REF] Benoît | Les canards de R 3[END_REF], Benoît [START_REF] Benoît | Systèmes lents-rapides dans R 3 et leurs canards[END_REF] and then Benoît [START_REF] Benoît | Canards de R 3[END_REF] in his PhD-thesis studied canard solutions in R 3 . In the article entitled "Systèmes lents-rapides dans R 3 et leurs canards," Benoît [6, p. 170] proved the existence of canards solution for three-dimensional singularly perturbed systems with two slow variables and one fast variable while using "Non-Standard Analysis"according to a theorem which stated that canard solutions exist in such systems provided that the pseudo singu-lar point 5 of the slow dynamics, i.e., of the reduced vector field is of saddle type.

Nearly twenty years later, Szmolyan and Wechselberger [START_REF] Szmolyan | Canards in R 3[END_REF] extended "Geometric Singular Perturbation Theory6 " to canards problems in R 3 and provided a "standard version" of Benoît's theorem [START_REF] Benoît | Systèmes lents-rapides dans R 3 et leurs canards[END_REF]. Very recently, Wechselberger [START_REF] Wechselberger | À propos de canards[END_REF] generalized this theorem for n-dimensional singularly perturbed systems with k slow variables and m fast (Eq. ( 1)). The method used by Szmolyan and Wechselberger [START_REF] Szmolyan | Canards in R 3[END_REF] and Wechselberger [START_REF] Wechselberger | À propos de canards[END_REF] require to implement a "desingularization procedure" which can be summarized as follows: first, they compute the normal form of such singularly perturbed systems (see Eq. [START_REF] Itoh | Duality of Memristors[END_REF] for dimension three and Eq. (48) for dimension four) which is expressed according to some coefficients (a and b for dimension three and ã, b and cj for dimension four) depending on the functions defining the original vector field [START_REF] Argémi | Approche qualitative d'un problème de perturbations singulières dans R 4[END_REF] and their partial derivatives with respect to the variables. Secondly, they project the "desingularized vector field" (originally called "normalized slow dynamics" by Eric Benoît [6, p. 166]) of such a normal form on the tangent bundle of the critical manifold. Finally, they evaluate the Jacobian of the projection of this "desingularized vector field" at the folded singularity (originally called pseudo singular points by José Argémi [1, p. 336]. This leads Szmolyan and Wechselberger [41, p. 427] and Wechselberger [47, p. 3298] to a "classification of folded singularities (pseudo singular points)". Thus, they show that for three-dimensional singularly perturbed systems such folded singularities is of saddle type if the following condition is satisfied: a < 0 while for four-dimensional singularly perturbed systems such folded singularities is of saddle type if ã < 0. Then, Szmolyan and Wechselberger [41, p. 439] and Wechselberger [47, p. 3304] establish their Theorem 4.1. which state that "In the folded saddle and in the folded node case singular canards perturb to maximal canard for sufficiently small ε". However, in their works neither Szmolyan and Wechselberger [START_REF] Szmolyan | Canards in R 3[END_REF] nor Wechselberger [START_REF] Wechselberger | À propos de canards[END_REF] do not provide (to our knowledge) the expression of these constants (a and ã) which are necessary to state the existence of canard solutions in such systems.

So, the aim of this work is first to provide the expression of these constants and then to show that they can be directly determined starting from the normalized slow dynamics and not from the projection of the "desingularized vector field" of the normal form. This method enables to state a unique "generic" condition for the existence of "canard solutions" for such three and four-dimensional singularly perturbed systems which is based on the stability of folded singularities of the normalized slow dynamics deduced from a well-known property of linear algebra. This unique condi-tion which is completely identical to that provided by Benoît [START_REF] Benoît | Systèmes lents-rapides dans R 3 et leurs canards[END_REF] and then by Szmolyan and Wechselberger [START_REF] Szmolyan | Canards in R 3[END_REF] and finally by Wechselberger [START_REF] Wechselberger | À propos de canards[END_REF] is "generic" since it is exactly the same for singularly perturbed systems of dimension three and four with only one fast variable. So, it provides a path to many applications.

In the very beginning of the seventies, Leon Chua [START_REF] Chua | Memristor -The Missing Circuit Element[END_REF] considered the three basic building blocks of an electric circuit: the capacitor, the resistor and the inductor as well as the three laws linking the four fundamental circuit variables, namely, the electric current i, the voltage v, the charge q and the magnetic flux φ. He thus concluded from the logical as well as axiomatic points of view, that it is necessary, for the sake of completeness, to postulate the existence of a fourth circuit element to which he gave the name memristor since it behaves like a nonlinear resistor with memory. On April 30 th 2008, Stan Williams and co-workers [START_REF] Strukhov | The missing memristor found[END_REF] announced in the journal Nature that the missing circuit element, postulated thirty-seven years before by Leon Chua has been found [START_REF] Ginoux | The Singing Arc: The Oldest Memristor? in Chaos, CNN, Memristors and Beyond: A Festschrift for Leon Chua[END_REF]. Since, the memristor has been subject to many studies and applications [START_REF] Ventra | Circuit elements with memory: memristors, memcapacitors and meminductors[END_REF][START_REF] Pershin | Experimental demonstration of associative memory with memristive neural networks[END_REF]. More particularly, memristor-based circuits have been used by Itoh and Chua [START_REF] Itoh | Memristors oscillators[END_REF][START_REF] Itoh | Duality of Memristors[END_REF], Muthuswamy and Kokate [START_REF] Muthuswamy | Memristorbased chaotic circuits[END_REF], Muthuswamy [START_REF] Muthuswamy | Implementing memristor based chaotic circuits[END_REF], Muthuswamy and Chua [START_REF] Muthuswamy | Simplest chaotic circuit[END_REF] and Fitch et al. [START_REF] Fitch | Hyperchaos In A Memristor-Based Modified Canonical Chua's Circuit[END_REF][START_REF] Fitch | Development of Memristor Based Circuits[END_REF] to construct dynamical systems whose solutions exhibit chaotic and hyperchaotic behavior [START_REF] Fitch | Hyperchaos In A Memristor-Based Modified Canonical Chua's Circuit[END_REF].

In a paper entitled "Duality of Memristors Circuits", Itoh and Chua [28, p. 1330001-15] gave the memristor canonical Chua's circuit equation (69) in the three-dimensional flux-linkage and charge phase space. Differentiating this Eq. ( 69) with respect to time they obtained memristor-based canonical Chua's circuit equation (73) in the four-dimensional current-voltage phase space 7 . In both cases, the φ -q characteristic curve of these circuits has been represented by a piecewise-linear function (Eq. [START_REF] Strukhov | The missing memristor found[END_REF] in Itoh and Chua [27, p. 3189] and Eq. (70) in Itoh and Chua [28, p. 1330001-15]). In their works, Itoh and Chua [START_REF] Itoh | Memristors oscillators[END_REF][START_REF] Itoh | Duality of Memristors[END_REF] have shown that the dynamical systems modeling such circuits possess at least one eigenvalue with a large negative real part. This specific feature is of great interest since it enables to consider memristor-based canonical Chua's circuits as slow-fast dynamical systems. So, there exists in the phase-space a slow manifold on which trajectories, solution of the dynamical system modeling the memristor circuit, evolve slowly and toward which nearby orbits contract exponentially in time in the normal directions.

For these memristor-based canonical circuits Itoh and Chua [START_REF] Itoh | Memristors oscillators[END_REF][START_REF] Itoh | Duality of Memristors[END_REF] have used a classical piecewise-linear function for the φ -q characteristic curve. However, Muthuswamy [START_REF] Muthuswamy | Implementing memristor based chaotic circuits[END_REF] and Fitch et al. [START_REF] Fitch | Hyperchaos In A Memristor-Based Modified Canonical Chua's Circuit[END_REF] have proposed to replace this piecewise linear characteristic curve by a smooth cubic nonlinear function. This enables to exhibit the existence of generic "canard solutions" in such Memristor Based Chaotic Circuits.

The outline of this paper is as follows. In Sec. 2, definitions of singularly perturbed system, critical manifold, reduced system, "constrained system", canard cycles, folded singularities and pseudo singular points are recalled. The method proposed in this article is presented in Sec. 3 & 4 for the case of three and four-dimensional singularly perturbed systems with only one fast variable. Existence of canard solution for the third and fourth-order Chua's memristor is established according to this method in Sec. 5 & 6.

DEFINITIONS 2.1. Singularly perturbed systems

According to Tikhonov [START_REF] Tikhonov | On the dependence of solutions of differential equations on a small parameter[END_REF], Pontryagin [START_REF] Pontryagin | The asymptotic behaviour of systems of differential equations with a small parameter multiplying the highest derivatives[END_REF], Jones [START_REF] Jones | Geometric Singular Perturbation Theory in Dynamical Systems[END_REF] and Kaper [START_REF] Kaper | An Introduction to Geometric Methods and Dynamical Systems Theory for Singular Perturbation Problems[END_REF] singularly perturbed systems are defined as:

⃗ x ′ = ε ⃗ f (⃗ x, ⃗ y, ε) , ⃗ y ′ = ⃗ g (⃗ x, ⃗ y, ε) . ( 1
)
where

⃗ x ∈ R k , ⃗ y ∈ R m , ε ∈ R + ,
and the prime denotes differentiation with respect to the independent variable t ′ . The functions ⃗ f and ⃗ g are assumed to be

C ∞ functions 8 of ⃗ x, ⃗ y and ε in U × I, where U is an open subset of R k × R m and I is an open interval containing ε = 0.
In the case when 0 < ε ≪ 1, i.e. ε is a small positive number, the variable ⃗ x is called slow variable, and ⃗ y is called fast variable. Using Landau's notation: O (ε p ) represents a function f of u and ε such that f (u, ε)/ε p is bounded for positive ε going to zero, uniformly for u in the given domain.

In general we consider that ⃗ x evolves at an O (ε) rate; while ⃗ y evolves at an O (1) slow rate. Reformulating system (1) in terms of the rescaled variable t = εt ′ , we obtain

⃗ x = ⃗ f (⃗ x, ⃗ y, ε) , ε ⃗ y = ⃗ g (⃗ x, ⃗ y, ε) .
(

) 2
The dot represents the derivative with respect to the new independent variable t.

The independent variables t ′ and t 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 ε ̸ = 0, and they are labeled singular perturbation problems when 0 < ε ≪ 1. The label "singular" stems in part from the discontinuous limiting behavior in system (1) as ε → 0.

Reduced slow system

In such case system (2) leads to a differential-algebraic system (D.A.E.) called reduced slow system whose dimension decreases from k + m = n to m. Then, the slow variable ⃗ x ∈ R k partially evolves in the submanifold M 0 called the critical manifold 9 . The reduced slow system is

⃗ x = ⃗ f (⃗ x, ⃗ y, ε) , ⃗ 0 = ⃗ g (⃗ x, ⃗ y, ε) .
(3)

Slow Invariant Manifold

The critical manifold is defined by

M 0 := { (⃗ x, ⃗ y) : ⃗ g (⃗ x, ⃗ y, 0) = ⃗ 0 } . ( 4
)
Such a normally hyperbolic invariant manifold (4) of the reduced slow system (3) persists as a locally invariant slow manifold of the full problem (1) for ε sufficiently small. This locally slow invariant manifold is O(ε) close to the critical manifold.

When D ⃗ x ⃗ f is invertible, thanks to the Implicit Function Theorem, M 0 is given by the graph of a C ∞ function ⃗ x = ⃗ G 0 (⃗ y) for ⃗ y ∈ D, where D ⊆ R k is a compact, simply connected domain and the boundary of D is a (k -1)-dimensional C ∞ submanifold 10 .
According to Fenichel [START_REF] Fenichel | Persistence and smoothness of invariant manifolds for flows[END_REF][START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF] theory if 0 < ε ≪ 1 is sufficiently small, then there exists a function ⃗ G (⃗ y, ε) defined on D such that the manifold

M ε := { (⃗ x, ⃗ y) : ⃗ x = ⃗ G (⃗ y, ε) } , ( 5
)
is locally invariant under the flow of system (1). Moreover, there exist perturbed local stable (or attracting) M a and unstable (or repelling) M r branches of the slow invariant manifold M ε . Thus, normal hyperbolicity of M ε is lost via a saddle-node bifurcation of the reduced slow system (3). Then, it gives rise to solutions of "canard" type.

Canards, singular canards and maximal canards

A canard is a solution of a singularly perturbed dynamical system (1) following the attracting branch M a of the slow invariant manifold, passing near a bifurcation point located on the fold of this slow invariant manifold, and then following the repelling branch M r of the slow invariant manifold.

A singular canard is a solution of a reduced slow system (3) following the attracting branch M a,0 of the critical manifold, passing near a bifurcation point located on the fold of this critical manifold, and then following the repelling branch M r,0 of the critical manifold.

A maximal canard corresponds to the intersection of the attracting and repelling branches M a,ε ∩ M r,ε of the slow manifold in the vicinity of a non-hyperbolic point.

According to Wechselberger [47, p. 3302]: "Such a maximal canard defines a family of canards nearby which are exponentially close to the maximal canard, i.e. a family of solutions of (1) that follow an attracting branch M a,ε of the slow manifold and then follow, rather surprisingly, a repelling/saddle branch M r,ε of the slow manifold for a considerable amount of slow time. The existence of this family of canards is a consequence of the non-uniqueness of M a,ε and M r,ε . However, in the singular limit ε → 0, such a family of canards is represented by a unique singular canard."

Canards are a special class of solution of singularly perturbed dynamical systems for which normal hyperbolicity is lost. Canards in singularly perturbed systems with two or more slow variables (⃗ x ∈ R k , k

2) and one fast variable (⃗ y ∈ R m , m = 1) are robust, since maximal canards generically persist under small parameter changes 11 .

Constrained system

In order to characterize the "slow dynamics", i.e. the slow trajectory of the reduced slow system (3) (obtained by setting ε = 0 in (2)), Floris Takens [START_REF] Takens | Constrained equations, a study of implicit differential equations and their discontinuous solutions, in Structural stability, the theory of catastrophes and applications in the sciences[END_REF] introduced the "constrained system" defined as follows:

⃗ x = ⃗ f (⃗ x, ⃗ y, 0) , D ⃗ y ⃗ g. ⃗ y = -(D ⃗ x ⃗ g. ⃗ f ) (⃗ x, ⃗ y, 0) , ⃗ 0 = ⃗ g (⃗ x, ⃗ y, 0) . ( 6
)
11 See Benoît [START_REF] Benoît | Systèmes lents-rapides dans R 3 et leurs canards[END_REF][START_REF] Benoît | Perturbation singulière en dimension trois : Canards en un point pseudo singulier noeud[END_REF], Szmolyan and Wechselberger [START_REF] Szmolyan | Canards in R 3[END_REF] and Wechselberger [START_REF] Wechselberger | Existence and Bifurcation of Canards in R 3 in the case of a Folded Node[END_REF][START_REF] Wechselberger | À propos de canards[END_REF].

Since, according to Fenichel [START_REF] Fenichel | Persistence and smoothness of invariant manifolds for flows[END_REF][START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF], the critical manifold ⃗ g (⃗ x, ⃗ y, 0) may be considered as locally invariant under the flow of system (1), we have:

d⃗ g dt (⃗ x, ⃗ y, 0) = 0 ⇐⇒ D ⃗ x ⃗ g. ⃗ x + D ⃗ y ⃗ g. ⃗ y = ⃗ 0.

By replacing ⃗

x by ⃗ f (⃗ x, ⃗ y, 0) leads to:

D ⃗ x ⃗ g. ⃗ f (⃗ x, ⃗ y, 0) + D ⃗ y ⃗ g. ⃗ y = ⃗ 0.
This justifies the introduction of the constrained system. Now, let adj(D ⃗ y ⃗ g) denote the adjoint of the matrix D ⃗ y ⃗ g which is the transpose of the co-factor matrix D ⃗ y ⃗ g, then while multiplying the left hand side of (6) by the inverse matrix (D ⃗ y ⃗ g) -1 obtained by the adjoint method we have:

⃗ x = ⃗ f (⃗ x, ⃗ y, 0) , det(D ⃗ y ⃗ g) ⃗ y = -(adj(D ⃗ y ⃗ g).D ⃗ x ⃗ g. ⃗ f ) (⃗ x, ⃗ y, 0) , ⃗ 0 = ⃗ g (⃗ x, ⃗ y, 0) . (7)

Normalized slow dynamics

Then, by rescaling the time by setting t = -det(D ⃗ y ⃗ g)τ we obtain the following system which has been called by Eric Benoît [6, p. 166] "normalized slow dynamics":

⃗ x = -det(D ⃗ y ⃗ g) ⃗ f (⃗ x, ⃗ y, 0) , ⃗ y = (adj(D ⃗ y ⃗ g).D ⃗ x ⃗ g. ⃗ f ) (⃗ x, ⃗ y, 0) , ⃗ 0 = ⃗ g (⃗ x, ⃗ y, 0) . ( 8
)
where the overdot now denotes the time derivation with respect to τ . Let's notice that José Argémi [START_REF] Argémi | Approche qualitative d'un problème de perturbations singulières dans R 4[END_REF] proposed to rescale time by setting t = -det(D ⃗ y ⃗ g)sgn(det(D ⃗ y ⃗ g))τ in order to keep the same flow direction in [START_REF] Benoît | Canards et enlacements[END_REF] as in [START_REF] Benoît | Canards de R 3[END_REF].

Desingularized vector field

By application of the Implicit Function Theorem, let suppose that we can explicitly express from Eq. ( 4), say without loss of generality, x 1 as a function ϕ 1 of the other variables. This implies that M 0 is locally the graph of a function ϕ 1 : R k → R m over the base U = (⃗ χ, ⃗ y) where ⃗ χ = (x 2 , x 3 , ..., x k ). Thus, we can span the "normalized slow dynamics" on the tangent bundle at the critical manifold M 0 at the pseudo singular point. This leads to the so-called desingularized vector field :

⃗ χ = -det(D ⃗ y ⃗ g) ⃗ f (⃗ χ, ⃗ y, 0) , ⃗ y = (adj(D ⃗ y ⃗ g).D ⃗ x ⃗ g. ⃗ f ) (⃗ χ, ⃗ y, 0) . (9)

Pseudo singular points and pseudo singular manifolds

As recalled by Guckenheimer and Haiduc [24, p. 91], pseudo-singular points have been introduced by the late José Argémi [START_REF] Argémi | Approche qualitative d'un problème de perturbations singulières dans R 4[END_REF] for low-dimensional singularly perturbed systems and are defined as singular points of the "normalized slow dynamics" [START_REF] Benoît | Canards et enlacements[END_REF]. Twenty-three years later, Szmolyan and Wechselberger [41, p. 428] called such pseudo singular points, folded singularities. In a recent publication entitled "A propos de canards" Wechselberger [47, p. 3295] proposed to define such singularities for n-dimensional singularly perturbed systems with k slow variables and m fast as the solutions of the following system:

det(D ⃗ y ⃗ g) = 0, (adj(D ⃗ y ⃗ g).D ⃗ x ⃗ g. ⃗ f ) (⃗ x, ⃗ y, 0) = ⃗ 0, ⃗ g (⃗ x, ⃗ y, 0) = ⃗ 0. (10)
Thus, for dimensions higher than three, his concept encompasses that of Argémi. Moreover, Wechselberger [47, p. 3296] proved that folded singularities form a (k -2)-dimensional manifold. Thus, for k = 2 the folded singularities are nothing else than the pseudo singular points defined by Argémi [START_REF] Argémi | Approche qualitative d'un problème de perturbations singulières dans R 4[END_REF]. While for k 3 the folded singularities are no more points but a (k -2)-dimensional manifold. Moreover, let's notice on the one hand that the original system (1) includes n = k + m variables and on the other hand, that the system (10) comprises p = 2m + 1 equations. So, we are faced to a system of p equations with n unknowns. If p < n the system is triangular and will necessarily have an infinite number of solutions that will be able to express in terms of the last unknowns. Since in this work we are only interested in three and four-dimensional singularly perturbed systems with m = 1 fast variable and with k slow variables we have p = 3 and n = k + 1 and so, n -p = k -2. Thus, we will examine the case k = 2 and k 3.

Pseudo singular points -Case k = 2

If k = 2 the number of variables of system (1) is equal to n = 3 and the number of equations is also equal to p = 3. So, all the variables (unknowns) of system (10) can be determined. The solutions of such system are called pseudo singular points. An example of such situation is given by the thirdorder Memristor-Based canonical oscillator analyzed in Sec. 5 and for which (m, k) = (1, 2). We will see in the next Sec. 3 that the stability analysis of these pseudo singular points will give rise to a condition for the existence of canard solutions in such systems.

Pseudo singular manifolds -Case k 3

If k 3 the number of variables of system ( 1) is equal to n = k + 1 and the number of equations is still equal to p = 3. So, only three variables (unknowns) of system ( 10) can be determined while the remaining k -2 unknowns are undetermined. The solution of such system takes the form of a (k -2)-dimensional manifold that we call pseudo singular manifold. An example of such situation is given by the fourth-order Memristor-Based canonical oscillator analyzed in Sec. 6 and for which (m, k) = (1, 3). We will see in Sec. 4 that for k 3 the stability analysis of this pseudo singular manifold will give rise to a condition (represented by a domain) for the existence of canard solutions in such systems.

THREE-DIMENSIONAL SINGULARLY PERTURBED SYSTEMS

A three-dimensional singularly perturbed dynamical system (2) with k = 2 slow variables and m = 1 fast may be written as:

ẋ1 = f 1 (x 1 , x 2 , y 1 ) , ( 11a
) ẋ2 = f 2 (x 1 , x 2 , y 1 ) , ( 11b
)
ε ẏ1 = g 1 (x 1 , x 2 , y 1 ) , ( 11c
)
where

⃗ x = (x 1 , x 2 ) t ∈ R 2 , ⃗ y = (y 1 ) ∈ R, 0 < ε ≪ 1
and the functions f i and g i are assumed to be C 2 functions of (x 1 , x 2 , y 1 ).

Critical Manifold

The critical manifold equation of system [START_REF] Chua | Memristor -The Missing Circuit Element[END_REF] is defined by setting ε = 0 in Eq. (11c). Thus, we obtain:

g 1 (x 1 , x 2 , y 1 ) = 0. ( 12
)
By application of the Implicit Function Theorem, let suppose that we can explicitly express from Eq. ( 11), say without loss of generality, x 1 as functions of the others variables:

x 1 = ϕ (x 2 , y 1 )

Constrained system

The constrained system of system ( 11) is obtained by equating to zero the time derivative of g 1 (x 1 , x 2 , y 1 ):

dg 1 dt = ∂g 1 ∂x 1 ẋ1 + ∂g 1 ∂x 2 ẋ2 + ∂g 1 ∂y 1 ẏ1 = 0 (13)
By replacing ẋi by f i (x 1 , x 2 , y 1 ) with i = 1, 2, Eq. ( 13) reads:

ẏ1 = - ∂g 1 ∂x 1 f 1 + ∂g 1 ∂x 2 f 2 ∂g 1 ∂y 1 . ( 14
)
So, we have the following constrained system:

ẋ1 = f 1 (x 1 , x 2 , y 1 ) , ẋ2 = f 2 (x 1 , x 2 , y 1 ) , ẏ1 = - ∂g 1 ∂x 1 f 1 + ∂g 1 ∂x 2 f 2 ∂g 1 ∂y 1 , 0 = g 1 (x 1 , x 2 , y 1 ) . ( 15
)

Normalized slow dynamics

By rescaling the time by setting t = -∂g 1 ∂y 1 τ we obtain the "normalized slow dynamics":

ẋ1 = -f 1 (x 1 , x 2 , y 1 ) ∂g 1 ∂y 1 = F 1 (x 1 , x 2 , y 1 ) , ẋ2 = -f 2 (x 1 , x 2 , y 1 ) ∂g 1 ∂y 1 = F 2 (x 1 , x 2 , y 1 ) , ẏ1 = ∂g 1 ∂x 1 f 1 + ∂g 1 ∂x 2 f 2 = G 1 (x 1 , x 2 , y 1 ) , 0 = g 1 (x 1 , x 2 , y 1 ) . ( 16
)
where the overdot now denotes the time derivation with respect to τ12 .

Desingularized vector field

Then, since we have supposed that x 1 may be explicitly expressed as a function ϕ (x 2 , y 1 ) on the others variables (Eq. 12), it can be used to project the "normalized slow dynamics" ( 16) on the tangent bundle of the critical manifold. Thus we obtain the so-called "desingularized vector field"

ẋ2 = -f 2 (x 1 , x 2 , y 1 ) ∂g 1 ∂y 1 (x 1 , x 2 , y 1 ) , ẏ1 = ∂g 1 ∂x 1 f 1 (x 1 , x 2 , y 1 ) + ∂g 1 ∂x 2 f 2 (x 1 , x 2 , y 1 ) . ( 17
)
in which x 1 must be replaced by ϕ (x 2 , y 1 ).

Pseudo-Singular Points

Pseudo-singular points are defined as singular points of the "normalized slow dynamics" [START_REF] Fenichel | Asymptotic stability with rate conditions II[END_REF], i.e. as the set of points for which we have:

g 1 (x 1 , x 2 , y 1 ) = 0, (18a) ∂g 1 ∂y 1 = 0, ( 18b
)
∂g 1 ∂x 1 f 1 + ∂g 1 ∂x 2 f 2 = 0. ( 18c
)
Remark 1. According to Argémi [START_REF] Argémi | Approche qualitative d'un problème de perturbations singulières dans R 4[END_REF], pseudo singular points are singular points of (18) but not necessarily singular points of [START_REF] Chua | Memristor -The Missing Circuit Element[END_REF]. In the following, we do not consider the case for which

f 1 (x 1 , x 2 , y 1 ) = f 2 (x 1 , x 2 , y 1 ) = g 1 (x 1 , x 2 , y 1 ) = 0.
Let's notice that contrary to the previous works we don't use the "desingularized vector field" (17) but the "normalized slow dynamics" [START_REF] Fenichel | Asymptotic stability with rate conditions II[END_REF].

Thus, the Jacobian matrix of system ( 16) reads:

J (F1,F2,G1) =           ∂F 1 ∂x 1 ∂F 1 ∂x 2 ∂F 2 ∂y 1 ∂F 2 ∂x 1 ∂F 2 ∂x 2 ∂F 2 ∂y 1 ∂G 1 ∂x 1 ∂G 1 ∂x 2 ∂G 1 ∂y 1           (19)

Benoît's generic hypothesis

In his famous papers, Eric Benoît [START_REF] Benoît | Les canards de R 3[END_REF][START_REF] Benoît | Systèmes lents-rapides dans R 3 et leurs canards[END_REF][START_REF] Benoît | Canards et enlacements[END_REF] made the following assumptions without loss of generality. First, he supposed that by a "standard translation" the pseudo-singular point can be shifted at the origin and that by a "standard rotation" of y 1 -axis that the critical manifold (12) is tangent to (x 2 , y 1 )-plane, so he had:

f 1 (0, 0, 0) = g 1 (0, 0, 0) = ∂g 1 ∂x 2 (0,0,0) = ∂g 1 ∂y 1 (0,0,0) = 0. (20)
Then, he made the following assumptions for the non-degeneracy of the pseudo-singular point:

f 2 (0, 0, 0) ̸ = 0 ; ∂g 1 ∂x 1 (0,0,0) ̸ = 0 ; ∂ 2 g 1 ∂y 2 1 (0,0,0) ̸ = 0 (21)
According to Benoît's generic hypotheses Eqs. [START_REF] Fruchard | Sur le retard à la bifurcation[END_REF][START_REF] Ginoux | Flow curvature method applied to canard explosion[END_REF], the Jacobian matrix [START_REF] Fitch | Development of Memristor Based Circuits[END_REF] reads:

J (F1,F2,G1) =        0 0 0 -f 2 ∂ 2 g 1 ∂x 1 ∂y 1 -f 2 ∂ 2 g 1 ∂x 2 ∂y 1 -f 2 ∂ 2 g 1 ∂y 2 1 a 31 a 32 a 33        (22)
where

a 3i = ∂g 1 ∂x 1 ∂f 1 ∂x i + f 2 ∂ 2 g 1 ∂x 2 ∂x i for i = 1, 2, a 33 = ∂g 1 ∂x 1 ∂f 1 ∂y 1 + f 2 ∂ 2 g 1 ∂x 2 ∂y 1 .
Thus, we have the following Cayley-Hamilton eigenpolynomial associated with such a Jacobian matrix [START_REF] Ginoux | Canards from Chua's circuit[END_REF] evaluated at the pseudo singular point, i.e., at the origin:

λ 3 -σ 1 λ 2 + σ 2 λ -σ 3 = 0 ( 23
)
It appears that σ 3 = |J (F1,F2,G1) | = 0 since one row of the Jacobian matrix ( 22) is null. So, the Cayley-Hamilton eigenpolynomial reduces to:

λ ( λ 2 -σ 1 λ + σ 2 ) = 0 (24)
Let λ i be the eigenvalues of the eigenpolynomial ( 24) and let's denote by λ 3 = 0 the obvious root of this polynomial. We have:

σ 1 =T r(J (F1,F2,G1) ) = λ 1 + λ 2 = ∂g 1 ∂x 1 ∂f 1 ∂y 1 , σ 2 = 3 ∑ i=1 J ii (F1,F2,G1) = λ 1 λ 2 =f 2 2 ( ∂ 2 g 1 ∂x 2 2 ∂ 2 g 1 ∂y 2 1 - ( ∂ 2 g 1 ∂x 2 ∂y 1 ) 2 ) + f 2 ∂g 1 ∂x 1 ( ∂ 2 g 1 ∂y 2 1 ∂f 1 ∂x 2 - ∂ 2 g 1 ∂x 2 ∂y 1 ∂f 1 ∂y 1
) .

(

) 25
where σ 1 = T r(J (F1,F2,G1) ) = p is the sum of all first-order diagonal minors of J (F1,F2,G1) , i.e. the trace of J (F1,F2,G1) and σ 2 = ∑ 3 i=1 J ii (F1,F2,G1) = q represents the sum of all second-order diagonal minors of J (F1,F2,G1) . Thus, the pseudo singular point is of saddle-type iff the following conditions C 1 and C 2 are verified:

C 1 : ∆ = p 2 -4q > 0, C 2 : q < 0. ( 26
)
Condition C 1 is systematically satisfied provided that condition C 2 is verified. Thus, the pseudo singular point is of saddle-type iff q < 0.

Canard existence in R 3

In an article entitled "Systèmes lents-rapides dans R 3 et leurs canards", Benoît [6, p. 171] has stated in the framework of "non-standard analysis" a theorem that can be written as follows:

Benoît's theorem [1983]

If the desingularized vector field [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF] has a pseudo singular point of saddle type, then system (11) exhibits a canard solution which evolves from the attractive part of the slow manifold towards its repelling part.

Proof. See [START_REF] Benoît | Systèmes lents-rapides dans R 3 et leurs canards[END_REF].

In his work, Benoît [6, p. 168] computed the trace T and determinant D of the Jacobian matrix J (F2,G1) associated with the two-dimensional desingularized vector field [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF]. Taking into account his generic hypotheses Eqs. [START_REF] Fruchard | Sur le retard à la bifurcation[END_REF][START_REF] Ginoux | Flow curvature method applied to canard explosion[END_REF] he found that:

T =T r(J (F2,G1) ) = λ 1 + λ 2 = ∂g 1 ∂x 1 ∂f 1 ∂y 1 , D = J (F2,G1) = λ 1 λ 2 =f 2 2 ( ∂ 2 g 1 ∂x 2 2 ∂ 2 g 1 ∂y 2 1 - ( ∂ 2 g 1 ∂x 2 ∂y 1 ) 2 ) + f 2 ∂g 1 ∂x 1 ( ∂ 2 g 1 ∂y 2 1 ∂f 1 ∂x 2 - ∂ 2 g 1 ∂x 2 ∂y 1 ∂f 1 ∂y 1
) .

(

) 27
from which he established that the pseudo singular point is of saddle type provided that D < 0. Then, Benoît [6, p. 171] stated his theorem.

In a paper entitled "Canards et enlacements", Benoît [START_REF] Benoît | Canards et enlacements[END_REF] stated, while using a standard polynomial change of variables (see Appendix A), that the original system [START_REF] Chua | Memristor -The Missing Circuit Element[END_REF] can be transformed into the following "normal version":

ẋ1 = ax 2 + by 1 + O ( x 1 , ε, x 2 2 , x 2 y 1 , y 2 1
) ,

ẋ2 = 1 + O (x 1 , x 2 , y 1 , ε) , ε ẏ1 = - ( x 1 + y 2 1 ) + O ( εx 1 , εx 2 , εy 1 , ε 2 , x 2 1 y 1 , y 3 1 , x 1 x 2 y 1
) , [START_REF] Itoh | Duality of Memristors[END_REF] where he established that

a = 1 2 f 2 2 ( ∂ 2 g 1 ∂x 2 2 ∂ 2 g 1 ∂y 2 1 - ( ∂ 2 g 1 ∂x 2 ∂y 1 ) 2 ) + 1 2 f 2 ∂g 1 ∂x 1 ( ∂ 2 g 1 ∂y 2 1 ∂f 1 ∂x 2 - ∂ 2 g 1 ∂x 2 ∂y 1 ∂f 1 ∂y 1 ) , b = - ∂g 1 ∂x 1 ∂f 1 ∂y 1 ,
A few years later, Szmolyan and Wechselberger [START_REF] Szmolyan | Canards in R 3[END_REF] gave a "standard version" of Benoît's theorem [START_REF] Benoît | Systèmes lents-rapides dans R 3 et leurs canards[END_REF] (see Benoît's theorem above) for threedimensional singularly perturbed systems with k = 2 slow variables and m = 1 fast. While using "standard analysis" and blow-up technique, Szmolyan and Wechselberger [41, p. 427] stated in their Lemma 2.1, while using "a smooth change of coordinates" (see Appendix A), that the original system (11) can be transformed into the "normal form" (28) from which they deduced that the condition for the pseudo singular point to be of saddle type is a < 0. Then, they proved the existence of canard solutions for the original system (11) according to their Theorem 4.1(a) presented below.

Theorem 1.

Assume system [START_REF] Itoh | Duality of Memristors[END_REF]. In the folded saddle and in the folded node case singular canards perturb to maximal canards solutions for sufficiently small ε.

Proof. See Szmolyan and Wechselberger [START_REF] Szmolyan | Canards in R 3[END_REF].

As previously recalled, the method presented in this paper doesn't use the "desingularized vector field" (17) but the "normalized slow dynamics" [START_REF] Fenichel | Asymptotic stability with rate conditions II[END_REF]. So, we have the following proposition:

Proposition 1.
If the normalized slow dynamics [START_REF] Fenichel | Asymptotic stability with rate conditions II[END_REF] has a pseudo singular point of saddle type, i.e. if the sum σ 2 of all second-order diagonal minors of the Jacobian matrix of the normalized slow dynamics (16) evaluated at the pseudo singular point is negative, i.e. if σ 2 < 0 then, according to Theorem 1, system (11) exhibits a canard solution which evolves from the attractive part of the slow manifold towards its repelling part.

Proof.

According to Eqs. [START_REF] Hurwitz | Über die Bedingungen, unter welchen eine Gleichung nur Wurzeln mit negativen reellen Theilen besitzt[END_REF][START_REF] Itoh | Memristors oscillators[END_REF] it is easy to verify that:

σ 1 = T r(J (F1,F2,G1) ) = T r(J (F2,G1) ) = T = λ 1 + λ 2 = -b, σ 2 = 3 ∑ i=1 J ii (F1,F2,G1) = J (F2,G1) = D = λ 1 λ 2 = 2a. (29)
So, the condition for which the pseudo singular point is of saddle type, i.e. σ 2 < 0 is identical to that proposed by Benoît [1983, p. 171] in his theorem, i.e. D < 0 and also to that provided by Szmolyan and Wechselberger [START_REF] Szmolyan | Canards in R 3[END_REF], i.e. a < 0. So, Prop. 1 can be used to state the existence of canard solution for such systems.

Of course, in the three-dimensional case the proof is obvious. We will see in the next Sect. 4 that for four-dimensional singularly perturbed systems this is not the case. Application of Proposition 1 to the three-dimensional memristor canonical Chua's circuits, presented in Sec. 5, will enable to prove the existence of generic "canard solutions" in such Memristor Based Chaotic Circuits.

FOUR-DIMENSIONAL SINGULARLY PERTURBED SYSTEMS

A four-dimensional singularly perturbed dynamical system (2) with k = 3 slow variables and m = 1 fast may be written as:

ẋ1 = f 1 (x 1 , x 2 , x 3 , y 1 ) , ( 30a
) ẋ2 = f 2 (x 1 , x 2 , x 3 , y 1 ) , ( 30b
) ẋ3 = f 3 (x 1 , x 2 , x 3 , y 1 ) , ( 30c
)
ε ẏ1 = g 1 (x 1 , x 2 , x 3 , y 1 ) , ( 30d
)
where ⃗ x = (x 1 , x 2 , x 3 ) t ∈ R 3 , ⃗ y = (y 1 ) ∈ R, 0 < ε ≪ 1 and the functions f i and g i are assumed to be C 2 functions of (x 1 , x 2 , x 3 , y 1 ).

Critical Manifold

The critical manifold equation of system [START_REF] Kaper | An Introduction to Geometric Methods and Dynamical Systems Theory for Singular Perturbation Problems[END_REF] is defined by setting ε = 0 in Eq. (30d). Thus, we obtain:

g 1 (x 1 , x 2 , x 3 , y 1 ) = 0. ( 31
)
By application of the Implicit Function Theorem, let suppose that we can explicitly express from Eq. ( 31), say without loss of generality, x 1 as functions of the others variables:

x 1 = ϕ 1 (x 2 , x 3 , y 1 ) . ( 32
)

Constrained system

The constrained system is obtained by equating to zero the time derivative of g 1 (x 1 , x 2 , x 3 , y 1 ):

dg 1 dt = ∂g 1 ∂x 1 ẋ1 + ∂g 1 ∂x 2 ẋ2 + ∂g 1 ∂x 3 ẋ3 + ∂g 1 ∂y 1 ẏ1 = 0 ( 33
)
By replacing ẋi by f i (x 1 , x 2 , x 3 , y 1 ) with i = 1, 2, 3, Eqs. ( 33) may be written as:

ẏ1 = - ∂g 1 ∂x 1 f 1 + ∂g 1 ∂x 2 f 2 + ∂g 1 ∂x 3 f 3 ∂g 1 ∂y 1 . ( 34
)
So, we have the following constrained system:

ẋ1 = f 1 (x 1 , x 2 , x 3 , y 1 ) , ẋ2 = f 2 (x 1 , x 2 , x 3 , y 1 ) , ẋ3 = f 3 (x 1 , x 2 , x 3 , y 1 ) , ẏ1 = - ∂g 1 ∂x 1 f 1 + ∂g 1 ∂x 2 f 2 + ∂g 1 ∂x 3 f 3 ∂g 1 ∂y 1 , 0 = g 1 (x 1 , x 2 , x 3 , y 1 ) . ( 35
)

Normalized slow dynamics

By rescaling the time by setting t = -∂g 1 ∂y 1 τ we obtain the "normalized slow dynamics":

ẋ1 = -f 1 (x 1 , x 2 , x 3 , y 1 ) ∂g 1 ∂y 1 = F 1 (x 1 , x 2 , x 3 , y 1 ) , ẋ2 = -f 2 (x 1 , x 2 , x 3 , y 1 ) ∂g 1 ∂y 1 = F 2 (x 1 , x 2 , x 3 , y 1 ) , ẋ3 = -f 3 (x 1 , x 2 , x 3 , y 1 ) ∂g 1 ∂y 1 = F 3 (x 1 , x 2 , x 3 , y 1 ) , ẏ1 = ∂g 1 ∂x 1 f 1 + ∂g 1 ∂x 2 f 2 + ∂g 1 ∂x 3 f 3 = G 1 (x 1 , x 2 , x 3 , y 1 ) , 0 = g 1 (x 1 , x 2 , x 3 , y 1 ) . ( 36
)
where the overdot now denotes the time derivation with respect to τ .

Desingularized vector field

Then, since we have supposed that x 1 may be explicitly expressed as a function ϕ 1 of the others variables [START_REF] Muthuswamy | Implementing memristor based chaotic circuits[END_REF], it can be used to project the "normalized slow dynamics" (36) on the tangent bundle of the critical manifold. So, we have:

ẋ2 = -f 2 (x 1 , x 2 , x 3 , y 1 ) ∂g 1 ∂y 1 , ẋ3 = -f 3 (x 1 , x 2 , x 3 , y 1 ) ∂g 1 ∂y 1 , ẏ1 = ∂g 1 ∂x 1 f 1 + ∂g 1 ∂x 2 f 2 + ∂g 1 ∂x 3 f 3 . ( 37
)
in which x 1 must be replaced by ϕ 1 (x 2 , x 3 , y 1 ).

Pseudo singular manifold

Pseudo singular manifold is defined as singular solution of the "normalized slow dynamics" (36), so we have:

g 1 (x 1 , x 2 , x 3 , y 1 ) = 0, (38a) ∂g 1 ∂y 1 = 0, ( 38b
)
∂g 1 ∂x 1 f 1 + ∂g 1 ∂x 2 f 2 + ∂g 1 ∂x 3 f 3 = 0. ( 38c
)
Remark 2. In the case of a four-dimensional singularly perturbed system with k = 3 slow variables and m = 1 fast, pseudo singular manifold forms a (k -2)-dimensional manifold, i.e. a 1-dimensional manifold since the system [START_REF] Robinson | Nonstandard Analysis[END_REF] comprises p = 3 equations and n = 4 variables (unknowns). So, in spite of having a pseudo singular point (x 1 , x2 , x3 , ỹ1 ) we have pseudo singular manifold represented by, say without loss of generality, (x 1 , x 2 , x3 , ỹ1 ), where x 2 is undetermined.

Let's notice again that contrary to the previous works we don't use the "desingularized vector field" (37) but the "normalized slow dynamics" [START_REF] Pershin | Experimental demonstration of associative memory with memristive neural networks[END_REF].

The Jacobian matrix of system (36) reads:

J (F1,F2,F3,G1) =                ∂F 1 ∂x 1 ∂F 1 ∂x 2 ∂F 1 ∂x 3 ∂F 1 ∂y 1 ∂F 2 ∂x 1 ∂F 2 ∂x 2 ∂F 2 ∂x 3 ∂F 2 ∂y 1 ∂F 3 ∂x 1 ∂F 3 ∂x 2 ∂F 3 ∂x 3 ∂F 3 ∂y 1 ∂G 1 ∂x 1 ∂G 1 ∂x 2 ∂G 1 ∂x 3 ∂G 1 ∂y 1                (39)

Extension of Benoît's generic hypothesis

Without loss of generality, it seems reasonable to extend Benoît's generic hypotheses introduced for the three-dimensional case to the four-dimensional case. So, first let's suppose that by a "standard translation" the pseudo singular manifold can be transformed into (0, x 2 , 0, 0) and that by a "standard rotation" of y 1 -axis that the critical manifold (31) is tangent to (x 2 , x 3 , y 1 )hyperplane, so we have

f 1 (0, x 2 , 0, 0) = g 1 (0, x 2 , 0, 0) = 0 ∂g 1 ∂x 2 (0,x2,0,0) = ∂g 1 ∂x 3 (0,x2,0,0) = ∂g 1 ∂y 1 (0,x2,0,0) = 0 (40)
Then, let's make the following assumptions for the non-degeneracy of the folded singularity:

f 2 (0, x 2 , 0, 0) ̸ = 0 ; ∂g 1 ∂x 1 (0,x2,0,0) ̸ = 0 ; ∂ 2 g 1 ∂y 2 1 (0,x2,0,0) ̸ = 0. ( 41
)
According to these generic hypotheses Eqs. [START_REF] Strukhov | The missing memristor found[END_REF][START_REF] Szmolyan | Canards in R 3[END_REF], the Jacobian matrix (39) reads:

J (F1,F2,F3,G1) =         0 0 0 0 a 21 a
        (42)
where

a 2i = -f 2 ∂ 2 g 1 ∂x i ∂y 1 for i = 1, 2, 3,
a 24 = -f 2 ∂ 2 g 1 ∂y 2 1 a 3i = -f 3 ∂ 2 g 1 ∂x i ∂y 1 for i = 1, 2, 3,
a 34 = -f 3 ∂ 2 g 1 ∂y 2 1 a 4i =f 1 ∂ 2 g 1 ∂x 1 ∂x i + ∂g 1 ∂x 1 ∂f 1 ∂x i + f 2 ∂ 2 g 1 ∂x 2 ∂x i + ∂g 1 ∂x 2 ∂f 2 ∂x i + f 3 ∂ 2 g 1 ∂x 3 ∂x i + ∂g 1 ∂x 3 ∂f 3 ∂x i for i = 1, 2, 3, a 44 =f 1 ∂ 2 g 1 ∂x 1 ∂y 1 + ∂g 1 ∂x 1 ∂f 1 ∂y 1 + f 2 ∂ 2 g 1 ∂x 2 ∂y 1 + ∂g 1 ∂x 2 ∂f 2 ∂y 1 + f 3 ∂ 2 g 1 ∂x 3 ∂y 1 + ∂g 1 ∂x 3 ∂f 3 ∂y 1 .
In his paper Wechselberger [START_REF] Wechselberger | À propos de canards[END_REF] stated that the determinant of the Jacobian matrix associated to the "desingularized vector field" and evaluated at a folded singularity, i.e. on the pseudo singular manifold is always zero 13 . Thus, we have the following Cayley-Hamilton eigenpolynomial associated with such a Jacobian matrix (41) evaluated on the pseudo singular manifold, i.e., at (0, x 2 , 0, 0):

λ 4 -σ 1 λ 3 + σ 2 λ 2 -σ 3 λ + σ 4 = 0, (43)
where σ 1 = T r(J (F1,F2,F3,G1) ) is the sum of all first-order diagonal minors of J (F1,F2,F3,G1) , i.e., the trace of J (F1,F2,F3,G1) , σ 2 represents the sum of all second-order diagonal minors of J (F1,F2,F3,G1) and σ 3 represents the sum of all third-order diagonal minors of J (F1,F2,F3,G1) . It appears that σ 4 = |J (F1,F2,F3,G1) | = 0 since one row of the Jacobian matrix (42) is null. So, the Cayley-Hamilton eigenpolynomial reduces to:

λ ( λ 3 -σ 1 λ 2 + σ 2 λ -σ 3 ) = 0. ( 44
)
But, according to Wechselberger [START_REF] Wechselberger | À propos de canards[END_REF], σ 3 vanishes on the pseudo singular manifold. Let's prove it:

Proof.

The sum of all third-order diagonal minors of J reads:

σ 3 = -f 2 ∂ 2 g 1 ∂x 2 ∂y 1 -f 2 ∂ 2 g 1 ∂x 3 ∂y 1 -f 2 ∂ 2 g 1 ∂y 2 1 -f 3 ∂ 2 g 1 ∂x 2 ∂y 1 -f 3 ∂ 2 g 1 ∂x 3 ∂y 1 -f 3 ∂ 2 g 1 ∂y 2 1 a 42 a 43 a 44
Then, while using a Laplace's expansion to compute this determinant, it's easy to show that it vanishes. So, the Cayley-Hamilton eigenpolynomial [START_REF] Tsuneda | A Gallery Of Attractors From Smooth Chua's Equation[END_REF] is thus reduced to

λ 2 ( λ 2 -σ 1 λ + σ 2 ) = 0 (45)
Let λ i be the eigenvalues of the eigenpolynomial (45) and let's denote by λ 3,4 = 0 the obvious double root of this polynomial. We have:

σ 1 =T r(J (F1,F2,F3,G1) ) = λ 1 + λ 2 = - ∂g 1 ∂x 1 ∂f 1 ∂y 1 , σ 2 = 3 ∑ i=1 J ii (F1,F2,F3,G1) = λ 1 λ 2 =2f 2 f 3 ( ∂ 2 g 1 ∂x 2 ∂x 3 ∂ 2 g 1 ∂y 2 1 - ∂ 2 g 1 ∂x 2 ∂y 1 ∂ 2 g 1 ∂x 3 ∂y 1 ) + f 2 2 ( ∂ 2 g 1 ∂x 2 2 ∂ 2 g 1 ∂y 2 1 -( ∂ 2 g 1 ∂x 2 ∂y 1 ) 2
)

+ f 2 ∂g 1 ∂x 1 ( ∂ 2 g 1 ∂y 2 1 ∂f 1 ∂x 2 - ∂ 2 g 1 ∂x 2 ∂y 1 ∂f 1 ∂y 1 ) + f 2 3 ( ∂ 2 g 1 ∂x 2 3 ∂ 2 g 1 ∂y 2 1 -( ∂ 2 g 1 ∂x 3 ∂y 1 ) 2
)

+ f 3 ∂g 1 ∂x 1 ( ∂ 2 g 1 ∂y 2 1 ∂f 1 ∂x 3 - ∂ 2 g 1 ∂x 3 ∂y 1 ∂f 1 ∂y 1 ) , ( 46
)
where σ 1 = T r(J (F1,F2,F3,G1) ) = p is is the sum of all first-order diagonal minors of J (F1,F2,F3,G1) , i.e. the trace of the Jacobian matrix J (F1,F2,F3,G1)

and σ 2 = ∑ 3 i=1 J ii (F1,F2,F3,G1) = q represents the sum of all second-order diagonal minors of J (F1,F2,F3,G1) .

Thus, the pseudo singular manifold is of saddle-type iff the following conditions C 1 and C 2 are verified:

C 1 : ∆ = p 2 -4q > 0, C 2 : q < 0. ( 47
)
Condition C 1 is systematically satisfied provided that condition C 2 is verified. Thus, the pseudo singular manifold is of saddle-type iff q < 0. But, as recalled previously, one coordinate is undetermined, say x 2 without loss of generality. So, the eigenvalues (46) of the characteristic polynomial are also functions of the variable x 2 and of the parameters of system [START_REF] Kaper | An Introduction to Geometric Methods and Dynamical Systems Theory for Singular Perturbation Problems[END_REF]. Now, let suppose that one parameter, say without loss of generality α 2 (see Sec. 6), modifies the nature of the pseudo singular manifold. Condition C 2 , i.e. q < 0 is then represented in the space (x 2 , α 2 ) by a straight line defining a region within which the pseudo singular points are of saddle type. In other words, it means that by choosing a value of the coordinate x 2 inside this region ensures that the pseudo singular point would be of saddle type.

Canard existence in R 4

In a paper entitled "A propos de canards" Wechselberger [START_REF] Wechselberger | À propos de canards[END_REF] stated, while using a standard polynomial change of variables, that any n-dimensional singularly perturbed systems with k slow variables (k 2) and m fast (m 1) (1) can be transformed into the following "normal form" (see Appendix B):

ẋ1 =ãx 2 + by 1 + O ( x 1 , x 2 2 , x 2 y 1 , y 2 1 ) + εO (x 1 , x 2 , x k , y 1 ) , ẋ2 =1 + O (x 1 , x 2 , y 1 , ε) , ẋj =c j + O (x 1 , x 2 , y 1 , ε) , j = 3, . . . , k ε ẏ1 =x 1 + y 2 1 + x 1 y 1 O (x 2 , . . . , x k ) + y 2 1 O (x 1 , y 1 ) + εO (x 1 , x 2 , y 1 , ε) (48)
which is a generalization of system [START_REF] Itoh | Duality of Memristors[END_REF]. We will establish in Appendix B for any four-dimensional singularly perturbed systems [START_REF] Kaper | An Introduction to Geometric Methods and Dynamical Systems Theory for Singular Perturbation Problems[END_REF] with k = 3 slow variables and m = 1 fast variable that

ã = 1 2 f 2 2 ( ∂ 2 g 1 ∂x 2 2 ∂ 2 g 1 ∂y 2 1 -( ∂ 2 g 1 ∂x 2 ∂y 1 ) 2 ) + 1 2 f 2 ∂g 1 ∂x 1 ( ∂ 2 g 1 ∂y 2 1 ∂f 1 ∂x 2 - ∂ 2 g 1 ∂x 2 ∂y 1 ∂f 1 ∂y 1 ) + 1 2 f 2 3 ( ∂ 2 g 1 ∂x 2 3 ∂ 2 g 1 ∂y 2 1 -( ∂ 2 g 1 ∂x 3 ∂y 1 ) 2 ) + 1 2 f 3 ∂g 1 ∂x 1 ( ∂ 2 g 1 ∂y 2 1 ∂f 1 ∂x 3 - ∂ 2 g 1 ∂x 3 ∂y 1 ∂f 1 ∂y 1 ) + f 2 f 3 ( ∂ 2 g 1 ∂x 2 ∂x 3 ∂ 2 g 1 ∂y 2 1 - ∂ 2 g 1 ∂x 2 ∂y 1 ∂ 2 g 1 ∂x 3 ∂y 1 ) , b = - ∂g 1 ∂x 1 ∂f 1 ∂y 1 ,
Remark 3. Let's notice that by posing f 3 = 0 in ã we find again a given in Sec. 3.7.

Thus, in his article entitled "A propos de canards" Wechselberger [47, p. 3304] has provided in the framework of "standard analysis" a generalization of Benoît's theorem [START_REF] Benoît | Systèmes lents-rapides dans R 3 et leurs canards[END_REF] (see Benoît's theorem above) for any n-dimensional singularly perturbed systems with k slow variables (k 2) and m fast (m 1). According to his Theorem 4.1(b) presented below he proved the existence of canard solutions for the original system (1).

Theorem 2.

In the folded saddle case of system (48) singular canards perturb to maximal canards solutions for sufficiently small ε ≪ 1.

Proof. See Wechselberger [START_REF] Wechselberger | À propos de canards[END_REF].

As previously recalled, the method presented in this paper doesn't use the "desingularized vector field" (37) but the "normalized slow dynamics" [START_REF] Pershin | Experimental demonstration of associative memory with memristive neural networks[END_REF]. So, we have the following proposition:

Proposition 2.
If the normalized slow dynamics (36) has a pseudo singular point of saddle type, i.e. if the sum σ 2 of all second-order diagonal minors of the Jacobian matrix of the normalized slow dynamics (36) evaluated at a pseudo singular point is negative, i.e. if σ 2 < 0 then, according to Theorem 2, system (30) exhibits a canard solution which evolves from the attractive part of the slow manifold towards its repelling part.

Proof.

By making some smooth changes of time and smooth changes of coordinates (see Appendix B) we brought the system (30) to the following "normal form":

ẋ1 = ãx 2 + by 1 + O ( x 1 , ε, x 2 2 , x 2 y 1 , y 2 1
) ,

ẋ2 = 1 + O (x 1 , x 2 , y 1 , ε) , ẋ3 = 1 + O (x 1 , x 2 , y 1 , ε) , ε ẏ1 = x 1 + y 2 1 + O ( εx 1 , εx 2 , εy 1 , ε 2 , x 2 1 y 1 , y 3 1 , x 1 x 2 y 1
) ,

Then, we deduce that the condition for the pseudo singular point to be of saddle type is ã < 0. According to Eqs. [START_REF] Wechselberger | Existence and Bifurcation of Canards in R 3 in the case of a Folded Node[END_REF] it is easy to verify that

σ 1 = T r(J (F1,F2,F3,G1) ) = λ 1 + λ 2 = -b, σ 2 = 3 ∑ i=1 J ii (F1,F2,F3,G1) = λ 1 λ 2 = 2ã.
So, the condition for which the pseudo singular point is of saddle type, i.e. σ 2 < 0 is identical to that proposed by Wechselberger [47, p. 3298] in his theorem, i.e. ã < 0. So, Prop. 2 can be used to state the existence of canard solution for such systems. Application of Proposition 2 to the four-dimensional memristor canonical Chua's circuits, presented in Sec. 6, will enable to prove the existence of generic "canards solutions" in such Memristor Based Chaotic Circuits.

THIRD-ORDER MEMRISTOR-BASED CANONICAL OSCILLATOR

Let's consider the Memristor-Based canonical Chua's circuit [START_REF] Itoh | Memristors oscillators[END_REF][START_REF] Itoh | Duality of Memristors[END_REF] containing five circuits elements: two passive capacitors, one passive inductor, one negative resistor, and one active Chua's flux controlled memristor (see Fig. 2).

FIG. 2.

Memristor-Based canonical Chua's circuit [START_REF] Itoh | Duality of Memristors[END_REF].

The parameter values used by Itoh and Chua [28, p. 1330001-14] i.e. are:

C 1 = 1 10 , C 2 = 1 0.47 , G = L = 1, a = -2.0, b = 4.0.

Flux-linkage and charge phase space

Applying Kirchhoff's circuit laws to the nodes A, B and the loop C of the circuit Fig. 2, Itoh and Chua [START_REF] Itoh | Memristors oscillators[END_REF][START_REF] Itoh | Duality of Memristors[END_REF] obtained the following set of differential equations, i.e., the following memristor based chaotic circuit:

C 1 dφ 1 dt = q 3 -k (φ 1 ) , C 2 dφ 2 dt = -q 3 + Gφ 2 , L dq 3 dt = φ 2 -φ 1 . ( 49
)
where the φ -q characteristic curve of the Chua's memristor is given by the following piecewise-linear function:

q = k (φ) = bφ + a -b 2 (|φ + 1| -|φ -1|) (50) By setting x = φ 1 , y = q 3 , z = φ 2 , ε = C 1 , β = 1 C 2 , γ = G C 2 and L = 1
the memristor based chaotic circuit (50) can be written:

dx dt = 1 ε [y -k (x)] , dy dt = z -x, dz dt = -βy + γz. ( 51
)
Following the works of Tsuneda [START_REF] Tsuneda | A Gallery Of Attractors From Smooth Chua's Equation[END_REF], let's replace the φ -q characteristic curve of the Chua's memristor q(φ) which is given by the piecewise-linear function (51) by a smooth cubic nonlinear function k(φ) = c 1 φ 3 + c 2 φ for which the parameters c 1 and c 2 are determined while using the least squares method. The square error between k(φ) and k(φ) is defined by:

S = ∫ d -d [ k(φ) -k(φ) ] 2 dφ ( 52
)
where [-d, d] is an interval for approximation. Let's note that in our case d is considered as a parameter such that |d| > 1. Solving ∂S/∂c 1 = 0 and ∂S/∂c 2 = 0, we find

c 1 = - 35(a -b) ( -1 + d 2 ) 2 16d 7 , c 2 = (a -b) ( 21 -50d 2 + 45d 4 ) 16d 5 + b.
(53)

Piecewise-linear and cubic nonlinearity

While still using the same parameter values as Itoh and Chua [28, p. 1330001-14] i.e.

C 1 = 1 10 , C 2 = 1 0.47 , G = L = 1, a = -2.0, b = 4.0,
the coefficients c 1 and c 2 have been chosen such that the extrema of both piecewise-linear and cubic nonlinearity characteristic curves substantially coincides as exemplified on Fig. 3. This condition is realized for d = 3 and Now, let's make the following variable changes in Eqs. (54) in order to apply the method presented in Sec. 3:

c 1 = 280 729 ; c 2 = - 26 27 .
x → z, y → -x, z → y.
Thus, we have:

dx dt = z -y, dy dt = βx + γy, dz dt = 1 ε [-x -k (z)] . ( 55
)
This last transformation will enable to compare the condition (given below) for the existence of canard solutions in system (55) with those given in our previous works entitled "Canards from Chua's circuits" [START_REF] Ginoux | Canards from Chua's circuit[END_REF].

Finally, let's replace the variables (x, y, z) by (x 1 , x 2 , y 1 ) and let's apply the method presented in Sec. 3 to the following system

dx 1 dt = y 1 -x 2 , dx 2 dt = βx 1 + γx 2 , dy 1 dt = 1 ε [-x 1 -k (y 1 )] .
(56)

Critical manifold and constrained system

The critical manifold of this system (56) is given by -x 1 -k(y 1 ) = 0. According to Eq. ( 15) the constrained system on the critical manifold reads:

ẋ1 = y 1 -x 2 , ẋ2 = βx 1 + γx 2 , ẏ1 = y 1 -x 2 -(c 1 y 3 1 + c 2 y 1 ) , 0 = -x 1 - ( c 1 y 3 1 + c 2 y 1
) .

(57)

Normalized slow dynamics

Then, by rescaling the time by setting t = -∂g 1 ∂y 1 τ = (3c 1 y 2 1 + c 2 ) we obtain the "normalized slow dynamics":

ẋ1 = (y 1 -x 2 ) ( 3c 1 y 2 1 + c 2 ) = F 1 (x 1 , x 2 , y 1 ) , ẋ2 = (βx 1 + γx 2 ) ( 3c 1 y 2 1 + c 2 ) = F 2 (x 1 , x 2 , y 1 ) , ẏ1 = x 2 -y 1 = G 1 (x 1 , x 2 , y 1 ) , 0 = -x 1 - ( c 1 y 3 1 + c 2 y 1
) .

(58)

Pseudo singular points

According to Eq. ( 18), the pseudo singular points of system (56) are:

x1 = ± 2c 2 3 √ -c 2 3c 1 , x2 = ∓ √ -c 2 3c 1 , ỹ1 = ∓ √ -c 2 3c 1 . ( 59
)
Let's notice that these pseudo singular points are independent of the parameter γ. The Jacobian matrix of system (58) evaluated at the pseudo singular points reads:

J (F1,F2,G1) =      0 0 0 0 0 -2γc 2 + 4βc 2 2 3 0 1 -1      (60)
Remark 4. Although, the pseudo singular points have not been shifted at the origin Benoît's generic hypotheses [START_REF] Fruchard | Sur le retard à la bifurcation[END_REF][START_REF] Ginoux | Flow curvature method applied to canard explosion[END_REF] are satisfied.

Canard existence in third-order memristor Chua's circuit

According to Eqs. [START_REF] Hurwitz | Über die Bedingungen, unter welchen eine Gleichung nur Wurzeln mit negativen reellen Theilen besitzt[END_REF] we find that:

p = σ 1 = T r [J] = -1, q = σ 2 = 2 3 c 2 (3γ -2βc 2 ) .
Thus, according to Prop. 1, the pseudo singular points are of saddle-type if and only if:

2 3 c 2 (3γ -2βc 2 ) < 0. ∆ = p 2 -4q > 0 and q < 0.
So, we have the following conditions C 1 and C 2 :

C 1 : ∆ = 1 + 4(-2c 2 )(γ - 2βc 2
3 ) > 0,

C 2 : q = 2c 2 (γ - 2βc 2 3 ) < 0. ( 61
)
Since the pseudo singular points are independent of the parameter γ let's choose γ as the "canard parameter" or "duck parameter". Obviously, it appears that if the condition C 2 is verified then the condition C 1 is de facto satisfied 14 . Finally, the pseudo singular points are of saddle-type if and only if we have:

γ saddle-node = 2βc 2 3 < γ. ( 62
)
where γ saddle-node represents the critical value of the parameter γ for which one of the two remaining eigenvalues λ 1 or λ 2 of the eigenpoynomial associated with the Jacobian matrix (60) vanishes. With this set of parameters ε = 1/10,

β = 0.47, a = -2, b = 4, d = 3, c 1 = 280/729, c 2 = -26/27 γ saddle-node = 2βc 2 3 ≈ -0.3.

Fixed points stability and Routh-Hurwitz' theorem

However, as pointed out in our previous works entitled "Canards from Chua's circuits" [START_REF] Ginoux | Canards from Chua's circuit[END_REF] the system (56) admits, except the origin, two fixed points, the stability of which could preclude the existence of "canards solutions". So, let's compute the fixed points of system (56) and analyze their stability. System (56) admits except the origin the following fixed points:

x * 1 = ± γ β √ γ -c 2 β c 1 β , x * 2 = y * 1 = ∓ √ γ -c 2 β c 1 β . ( 63
)
The eigenpolynomial equation of the Jacobian matrix of system (56) evaluated at these fixed points (63) reads:

ελ 3 + λ 2 ( 3γ β -γε -2c 2 ) + λ(1 - 3γ 2 β + βε + 2γc 2 ) + 2(γ -βc 2 ) = 0 (64)
Let suppose that all the parameters are fixed except γ, i.e. the "duck parameter". There are two methods to analyze the stability of fixed points as functions of the "duck parameter" value. The first is to solve the above third degree eigenpolynomial equation ( 64) with the Cardano's method and the second consists in using the Routh-Hurwitz' theorem [START_REF] Routh | A Treatise on the Stability of a Given State of Motion: Particularly Steady Motion[END_REF][START_REF] Hurwitz | Über die Bedingungen, unter welchen eine Gleichung nur Wurzeln mit negativen reellen Theilen besitzt[END_REF]. This latter method enables to easier analyze the stability of the fixed points as functions of a parameter. According to Routh-Hurwitz' theorem, the eigenpolynomial equation can be written as:

a 3 λ 3 + a 2 λ 2 + a 1 λ + a 0 = 0.
It states that if D 1 = a 1 and D 2 = a 1 a 2 -a 0 a 3 are both positive then eigenpolynomial equation would have eigenvalues with negative real parts. In other words, if D 1 and D 2 are positive the fixed points will be stable. In the case of the eigenpolynomial equation (64) we have:

D 1 =1 - 3γ 2 β + βε + 2γc 2 , D 2 = -2ε (γ -βc 2 ) + ( 3γ β -γε -2c 2 ) ( 1 - 3γ 2 β + βε + 2γc 2
) .

(

By setting: ε = 1/10, β = 0.47, a = -2, b = 4, d = 3, c 1 = 280/729, c 2 = -26/27 and while considering that the "duck parameter" γ can vary, D 1 and D 2 are respectively polynomial equations of degree two and three in γ. These quadratic and cubic functions D 1 and D 2 have been plotted on Fig. 5. One can see that between the lower limit called γ saddle-node and, the upper limit called γ Hopf corresponding to the value of the parameter γ for which the real parts of both complex eigenvalues vanishes (see proof in Appendix C), D 1 and D 2 are strictly positive. So, for γ ∈ [γ saddle-node , γ Hopf ] (purple rectangle on Fig. 5) the fixed points are stable while for γ > γ Hopf they are unstable. With this set of parameters, γ Hopf ≈ 0.274. Thus, it appears from what precedes and from Prop. 1 that "canards solutions" may be observed in system (56) for γ Duck values such that:

γ saddle-node = 2βc 2 3 < γ Hopf < γ Duck ( 66
)
On Fig. 6, numerical "canards solutions" and slow manifold of system (56) have been plotted for the "duck parameter" γ Duck = 0.3275 (all other parameters are the same as indicated above). Due to the symmetry of the system (56), any of the two pseudo singular points plotted in green on Fig. 6 was chosen as initial condition.

Particular case

In a previous work entitled "Canards from Chua's circuits", Ginoux et al. [START_REF] Ginoux | Canards from Chua's circuit[END_REF] have studied the system (56) with the following particular parameters:

γ = β = α c 1 = 1 3 c 2 = -1
First, let's replace these parameters in the above conditions C 1 and C 2 (61). We have:

C 1 : ∆ = 1 + 40α 3 > 0, C 2 : q = - 10α 3 < 0.
Obviously, if α > 0, then both conditions C 1 and C 2 are verified. This is exactly the result provided by Itoh and Chua [START_REF] Itoh | Canards and chaos in nonlinear systems, Circuits and Systems[END_REF] as it has been noticed in Ginoux et al. [22, p. 1330010-4]. However, it has been also remarked in our same previous paper [START_REF] Ginoux | Canards from Chua's circuit[END_REF] that the system (56) admits, except the origin, two fixed points, the stability of which could preclude the existence of "canards solutions". By setting: γ = β = α, c 1 = 1 3 and c 2 = -1 in Eq. ( 63) we find again the fixed points obtained by Ginoux et al. [22, p. 1330010-4]:

x * 1 = ± √ 6, x * 2 = ∓ √ 6, y * 1 = ∓ √ 6.
Moreover, still using the Routh-Hurwitz' theorem and by setting:

γ = β = α, c 1 = 1
3 and c 2 = -1 in Eq. ( 65) we find that:

D 1 = 1 -5α + εα, D 2 = α 2 ( 5ε -ε 2 ) -25α + 5.
Functions D 1 and D 2 have been plotted on Fig. 7 on which one can see that between the lower limit called α saddle-node and, the upper limit called α Hopf corresponding to the value of the parameter α for which the real parts of both complex eigenvalues vanishes, D 1 and D 2 are strictly positive. So, for α ∈ [α saddle-node , α Hopf ] (purple rectangle on Fig. 7) the fixed points are stable while for α > α Hopf , i.e. α > 1/5 they are unstable.

Thus, it appears from what precedes and from Prop. 3 that "canards solutions" may be observed in system (56) provided that: This is exactly the result obtained by Ginoux et al. [22, p. 1330010-6]. The phase portrait of system (56) with this set of parameter values has already been published by Ginoux et al. [START_REF] Ginoux | Canards from Chua's circuit[END_REF].

α saddle-node = 0 < α Hopf < α Duck

FOURTH-ORDER MEMRISTOR-BASED CANONICAL OSCILLATOR

Let's consider again the Memristor-Based canonical Chua's circuit [START_REF] Itoh | Memristors oscillators[END_REF][START_REF] Itoh | Duality of Memristors[END_REF]. By adding an inductor in parallel with conductance -G, Fitch et al. [START_REF] Fitch | Hyperchaos In A Memristor-Based Modified Canonical Chua's Circuit[END_REF] have modified this circuit in order to obtain a fourth-order Memristor-Based canonical oscillator (see Fig. 8).

L 2 i L2 C 2 v 2 v 1 C 1 R i L1 -G L 1 ( ) W φ node1 node2 Loop1 FIG. 8.
Memristor canonical Chua's circuit [START_REF] Fitch | Hyperchaos In A Memristor-Based Modified Canonical Chua's Circuit[END_REF].

Flux-linkage and charge phase space

Applying Kirchhoff's circuit laws to the nodes 1, 2 and the loop 1 of the circuit Fig. 8, Fitch et al. [START_REF] Fitch | Hyperchaos In A Memristor-Based Modified Canonical Chua's Circuit[END_REF] obtained the following set of differential equations, i.e., the following memristor based chaotic circuit:

C 1 dφ 1 dt = Rq 1 -k (φ 1 ) , C 2 dφ 2 dt = -q 2 + Gφ 2 -q 1 , L 1 dq 1 dt = φ 2 -φ 1 -Rq 1 , L 2 dq 3 dt = φ 2 , (67)
where the classical piecewise-linear function k(φ) of the Chua's memristor (50) has been replaced by the cubic k(φ) = c 1 φ 3 + c 2 φ.

By setting x = φ 1 , y = q 1 , z = φ 2 , u = q 2 , C 1 = ε, C 2 = 1, β 1 = 1 L 1 , β 2 = 1 L 2 , G = -α 2 and R = 1 the memristor based chaotic circuit (67)
can be written:

ε dx dt = y -k (x) , dy dt = -u -α 2 z -y, dz dt = β 1 (z -x -y) , du dt = β 2 z. ( 68
)
Now, let's make the following variable changes in Eqs. (68) in order to apply the method presented in Sec. 4:

x → u, y → x, u → y.
Thus, we have:

dx dt = β 1 (z -x -u) , dy dt = β 2 z, dz dt = -y -α 2 z -x, ε du dt = x -k (u) . ( 69
)
Let's notice that system (69) is exactly identical to that studied by Ginoux et al. [START_REF] Ginoux | Canards from Chua's circuit[END_REF]. Thus, condition (we will provide below) for the existence of canard solutions in system (69) will be compared to that given in our previous works entitled "Canards from Chua's circuits" [START_REF] Ginoux | Canards from Chua's circuit[END_REF].

Finally, let's replace the variables (x, y, z, u) by (x 1 , x 2 , x 3 , y 1 ) and let's apply the method presented in Sec. 4 to the following system (56) where k(y 1 ) = c 1 y 3 1 + c 2 y 1 .

dx 1 dt = β 1 (x 3 -x 1 -y 1 ) , dx 2 dt = β 2 x 3 , dx 3 dt = -x 2 -α 2 x 3 -x 1 , ε dy 1 dt = x 1 -k (y 1 ) .
(70)

Critical manifold and contrained system

The critical manifold of this system (70) is given by x 1 -k(y 1 ) = 0. According to Eq. ( 35) the constrained system on the critical manifold reads:

dx 1 dt = β 1 (x 3 -x 1 -y 1 ) , dx 2 dt = β 2 x 3 , dx 3 dt = -x 2 -α 2 x 3 -x 1 , dy 1 dt = - β 1 (x 3 -x 1 -y 1 ) -(3c 1 y 2 1 + c 2 ) , 0 = x 1 -k (y 1 ) .
(71)

Normalized slow dynamics

Then, by rescaling the time by setting t = -∂g 1 ∂y 1 τ = (3c 1 y 2 1 + c 2 ) we obtain the "normalized slow dynamics":

dx 1 dt = β 1 (x 3 -x 1 -y 1 ) ( 3c 1 y 2 1 + c 2 ) = F 1 (x 1 , x 2 , x 3 , y 1 ) , dx 2 dt = β 2 x 3 ( 3c 1 y 2 1 + c 2 ) = F 2 (x 1 , x 2 , x 3 , y 1 ) , dx 3 dt = (-x 2 -α 2 x 3 -x 1 ) ( 3c 1 y 2 1 + c 2 ) = F 3 (x 1 , x 2 , x 3 , y 1 ) , dy 1 dt = β 1 (x 3 -x 1 -y 1 ) = G 1 (x 1 , x 2 , x 3 , y 1 ) , 0 = x 1 -k (y 1 ) .
(72)

Pseudo singular manifold

According to Eqs. [START_REF] Robinson | Nonstandard Analysis[END_REF], the pseudo singular manifold of system (70) is defined by:

(x 1 , x 2 , x3 , ỹ1 ) = ( ± 2c 2 3 √ -c 2 3c 1 , x 2 , ±( 2c 2 3 + 1) √ -c 2 3c 1 , ± √ -c 2 3c 1 ) (73)
Let's notice that x2 is undetermined. In "Canards from Chua's circuit", Ginoux et al. [START_REF] Ginoux | Canards from Chua's circuit[END_REF] have arbitrarily chosen x2 = 0. We will see in the following that this choice does not affect their results.

The Jacobian matrix of system (72) evaluated at (x ± 1 , x 2 , x± 3 , ỹ± 1 ) reads:

J (F1,F2,F3,G1) =         0 0 0 0 0 0 0 6β 2 c 1 x3± ỹ1± 0 0 0 -6c 1 ( α 2 x3± + x 2 + c 1 ỹ3 1± + c 2 ỹ1±
) ỹ1±

-β 1 0 β 1 -β 1         (74)
Remark 5. Although, the pseudo singular manifold has not been transformed into (0, x 2 , 0, 0) extension of Benoît's generic hypotheses [START_REF] Strukhov | The missing memristor found[END_REF][START_REF] Szmolyan | Canards in R 3[END_REF] are satisfied.

Canard existence in fourth-order memristor Chua's circuit

According to Eqs. [START_REF] Wechselberger | Existence and Bifurcation of Canards in R 3 in the case of a Folded Node[END_REF] we find that:

p = T r(J) = -β 1 , q = σ 2 = +6β 1 c 1 (α 2 x3± + x 2 + x1± ) ỹ1± (75)
Thus, the conditions C 1 and C 2 for (x ± 1 , x 2 , x± 3 , ỹ± 1 ) to be of saddle type reads:

C 1 : ∆ = β 1 [β 1 -24c 1 (α 2 x3± + x 2 + x1± ) ỹ1± ] > 0, C 2 : q = +6β 1 c 1 (α 2 x3± + x 2 + x1± ) ỹ1± < 0. ( 76
)
Then, due to the nature (±) of the pseudo singular manifold (73) we have two cases corresponding to the positive and negative values.

Positive case

Let's consider the positive case for which the pseudo singular manifold (73) can be written as:

( x+ 1 , x 2 , x+ 3 , ỹ+ 1 ) = ( + 2c 2 3 √ -c 2 3c 1 , x 2 , +( 2c 2 3 + 1) √ -c 2 3c 1 , + √ -c 2 3c 1
) .

Conditions C 1 and C 2 reads then:

C 1 : α 2 x3+ + x 2 + x1+ < β 1 24c 1 ỹ1+ , ( 77a
)
C 2 : α 2 x3+ + x 2 + x1+ < 0. ( 77b
)
Obviously, since the right hand side of the first inequality (77a) is positive (β 1 > 0, c 1 > 0 and ỹ1+ > 0), both conditions are satisfied provided that the condition C 2 is verified. So, to have pseudo singular manifold of saddle type, the straight line α 2 x3+ + x 2 + x1+ must verify:

α 2 x3+ + x 2 + x1+ < 0. ( 78
)
By taking into account the above preliminary result and while fixing all the parameters excepted α 2 , this straight line (D + ) can be plotted in the plane (x 2 , α 2 ) and reads:

(D + ) : α 2 ( 2c 2 3 + 1 ) √ -c 2 3c 1 + x 2 + 2c 2 3 √ -c 2 3c 1 < 0. ( 79
)
Let's notice that for:

α 2 = 0, x 2 = - 2c 2 3 √ -c 2 3c 1 ,
x 2 = 0, α 2 = -2c 2 2c 2 + 3 .

Negative case

Let's consider the negative case for which the pseudo singular manifold (73) can be written as:

( x- 1 , x 2 , x- 3 , ỹ- 1 ) = ( - 2c 2 3 √ -c 2 3c 1 , -( 2c 2 3 + 1) √ -c 2 3c 1 , - √ -c 2 3c 1
) .

Conditions C 1 and C 2 reads then:

C 1 : β 1 24c 1 ỹ1- < α 2 x3-+ x 2 + x1-, ( 80a
)
C 2 : 0 < α 2 x3-+ x 2 + x1-. ( 80b
)
Obviously, since the left hand side of the first inequality (80a) is negative (β 1 > 0, c 1 > 0 and ỹ1-< 0), both conditions are satisfied provided that the condition C 2 is satisfied. So, to have pseudo singular manifold of saddle type, the straight line α 2 x3-+ x 2 + x1-must verify:

α 2 x3-+ x 2 + c 1 ỹ3 1-+ c 2 ỹ1-> 0. (81)
By taking into account the above preliminary result and while fixing all the parameters excepted α 2 , this straight line (D -) can be plotted in the plane (x 2 , α 2 ) and reads:

(D -) : -α 2 ( 2c 2 3 + 1 ) √ -c 2 3c 1 + x 2 - 2c 2 3 √ -c 2 3c 1 < 0. ( 82
)
Let's notice that for:

α 2 = 0, x 2 = 2c 2 3 √ -c 2 3c 1 , x 2 = 0, α 2 = - 2c 2 2c 2 + 3 .
On Fig. 9a & 9b, the straight lines (D + ) and (D -) have been plotted in blue and in red (respectively). Thus, the region within which the pseudo singular points are of saddle type corresponds to the cyan triangle. A zoom of Fig. 9a is presented on Fig. 9b. Let's notice on the one hand that the point (α 2 = 0.9, x 2 = 0) arbitrarily chosen by Ginoux et al. , x 2 = 0) which corresponds exactly with the condition stated in Ginoux et al. [START_REF] Ginoux | Canards from Chua's circuit[END_REF] and for which canard solutions have been observed in Chua's system 4D (70) according to Prop. 2. In other words, to have a pseudo singular point of saddle type at x 2 = 0, α 2 < -2c2 2c2+3 . To confirm this fact, the two nonzero eigenvalues of the characteristic polynomial associated with the Jacobian matrix (74) evaluated at (x 1 , x 2 , x3 , ỹ1 ) (73) have been computed for α 2 = 0.9 and for the corresponding values of x 2 which have been taken equal to zero by Ginoux et al. [START_REF] Ginoux | Canards from Chua's circuit[END_REF] but which is in fact very small x 2 = ∓0.01. We have found that the two nonzero real eigenvalues are of opposite sign what corresponds to the case of pseudo singular points of saddle type.

So, the value of the "duck parameter" α 2 for which the pseudo singular points are of saddle-type is defined by:

α 2 < α 2saddle-node = - 2c 2 3 + 2c 2 . ( 83
)
where α 2saddle-node represents the critical value of the parameter α 2 for which one of the two remaining eigenvalues λ 1 or λ 2 of the eigenpoynomial associated with the Jacobian matrix (74) vanishes. With this set of parameters ε = 1/10.1428, β 1 = 0.121, β 2 = 0.0047, c 1 = 0.393781 and c 2 = -0.72357, α 2saddle-node = -2c 2 3 + 2c 2 ≈ 0.932.

Fixed points stability and Routh-Hurwitz' theorem

However, as pointed out in the previous Sect. 5.7 the system (70) admits the origin O(0, 0, 0, 0) as fixed point, the stability of which could preclude the existence of "canards solutions". The eigenpolynomial equation of the Jacobian matrix of system (70) evaluated at this fixed point reads:

a 4 λ 4 + a 3 λ 3 + a 2 λ 2 + a 1 λ + a 0 = 0 ( 84
)
where

a 0 = (1 + c 2 ) β 1 β 2 , a 1 = c 2 ((1 + α 2 ) β 1 + β 2 ) + β 1 (α 2 + εβ 2 ) , a 2 = (1 + ε + εα 2 ) β 1 + c 2 (α 2 + β 1 ) + εβ 2 , a 3 = c 2 + ε (α 2 + β 1 ) , a 4 = ε.
Let suppose that all the parameters are fixed except α 2 , i.e. the "duck parameter" and, let's make use again of the Routh-Hurwitz' theorem [START_REF] Routh | A Treatise on the Stability of a Given State of Motion: Particularly Steady Motion[END_REF][START_REF] Hurwitz | Über die Bedingungen, unter welchen eine Gleichung nur Wurzeln mit negativen reellen Theilen besitzt[END_REF]. Thus, it states that if

D 1 = a 1 , D 2 = a 1 a 2 -a 0 a 3 and D 3 = a 1 a 2 a 3 -a 0 a 2 3 -a 2
1 a 4 are all positive then eigenpolynomial equation would have eigenvalues with real negative parts. In other words, if D 1 , D 2 and D 3 are positive the fixed point will be stable.

By setting ε = 1/10.1428, β 1 = 0.121, β 2 = 0.0047, c 1 = 0.393781 and c 2 = -0.72357 and while considering that the "duck parameter" α 2 can vary, the functions D 1 , D 2 and D 3 have been plotted on Fig. 10. One can see that between the lower limit called α 2Hopf corresponding to the value of the parameter α 2 for which the real parts of both complex eigenvalues vanishes (see Proof in the Appendix D.) and the upper limit called α 2saddle-node , D 1 and D 3 are negative while D 2 is positive. So, in this interval, the fixed point is unstable. With this set of parameters, α 2Hopf ≈ 0.0451 and α 2saddle-node = -2c 2 3 + 2c 2 ≈ 0.932.

Thus, we deduce from what precedes and from Prop. 2 that "canards solutions" may be observed in system (70) provided that:

α 2Hopf < α 2Duck < α 2saddle-node = -2c 2 3 + 2c 2 (85)
On Figs. 11 & 12, numerical "canards solutions" and critical manifold of system (70) have been plotted for the "duck parameter" α 2Duck = 0.1 (all other parameters are the same as indicated above). Due to the symmetry of the system (70), any of the two pseudo singular points plotted in green on Figs. 11 & 12 was chosen as initial condition.

DISCUSSION

In this work we have proposed an alternative method for determining the condition of existence of "canard solutions" for three and four-dimensional singularly perturbed systems with only one fast variable in the folded saddle case. This method enables to highlight a unique generic condition (σ 2 < 0) for the existence of "canard solutions" for such three and fourdimensional singularly perturbed systems which is based on the stability of folded singularities of the normalized slow dynamics deduced from a well-known property of linear algebra. It has been stated that this unique generic condition was perfectly identical to that provided by Benoît [START_REF] Benoît | Systèmes lents-rapides dans R 3 et leurs canards[END_REF] and then by Szmolyan and Wechselberger [START_REF] Szmolyan | Canards in R 3[END_REF] and finally by Wechselberger [START_REF] Wechselberger | À propos de canards[END_REF]. Finally, it has been established that this condition is "generic" since it is exactly the same for singularly perturbed systems of dimension three and four with only one fast variable. Application of this method to the famous three and four-dimensional memristor canonical Chua's circuits for which the classical piecewise-linear characteristic curve has been replaced by a smooth cubic nonlinear function according to the least squares method has enabled to show the existence of "canards solutions" in such Memristor Based Chaotic Circuits. However, in this paper, only the case of pseudo singular points of saddletype has been analyzed. Of course, the case of pseudo singular points of node-type could be also studied with the same method. Moreover, this method could be successfully used for proving the existence of "canard solutions" in four-dimensional singularly perturbed systems with two fast variables such as the famous Hodgkin-Huxley model or in the so-called coupled FitzHugh-Nagumo system. In a future work we will state that the existence of canard solutions in such systems can be established according to the same unique generic condition (σ 2 < 0).

ACKNOWLEDGEMENTS

We would like to thank to Ernesto Pérez Chavela for previous discussions related with this work. The authors are partially supported by a MINECO/FEDER grant number MTM2008-03437. The second author is partially supported by a MICINN/FEDER grants numbers MTM2009-03437 and MTM2013-40998-P, by an AGAUR grant number 2014SGR-

A = 1 2 ∂g 1 ∂x ∂ 2 g 1 ∂z 2 , B = 1 4 [ ∂ 2 g 1 ∂y 2 ∂ 2 g 1 ∂z 2 - ( ∂ 2 g 1 ∂y∂z ) 2 ] , C = - 1 2 ∂ 2 g 1 ∂y∂z , D = - 1 2 ∂ 2 g 1 ∂z 2 .
(A-9)

Finally, we deduce: This is the result established by Benoît [START_REF] Benoît | Canards et enlacements[END_REF] and presented in Sec. 3.7.

a = 1 2 f 2 2 ( ∂ 2

B. Normal form of 4D singularly perturbed systems with one fast variable

Let's consider the four-dimensional singularly perturbed dynamical system [START_REF] Kaper | An Introduction to Geometric Methods and Dynamical Systems Theory for Singular Perturbation Problems[END_REF] with k = 3 slow variables and m = 1 fast and let's make the following change of variables:

x 1 = α 2 x, x 2 = αy, x 3 = αz, y 1 = αu where α ≪ 1.

(A-11)

By taking into account extension of Benoît's generic hypothesis Eqs. [START_REF] Strukhov | The missing memristor found[END_REF][START_REF] Szmolyan | Canards in R 3[END_REF] and while using Taylor series expansion the system [START_REF] Kaper | An Introduction to Geometric Methods and Dynamical Systems Theory for Singular Perturbation Problems[END_REF] and u by the right-hand-side of (A-14) and by identifying with the following system in which we have posed: (ε/α 2 ) = ϵ: This is the result we established in Sec. 4.7. Moreover, let's notice that by posing f 3 = 0 in ã we find again a given in Sec. 3.7.

Ẋ = ãY + bU + O ( X,
Routh-Hurwitz' theorem and their application to the determination of the Hopf bifurcation parameter-value in the case of three and four-dimensional singularly perturbed system are presented in this appendix.

C. Routh-Hurwitz's theorem for 3D systems

According to [START_REF] Ginoux | The Singing Arc: The Oldest Memristor? in Chaos, CNN, Memristors and Beyond: A Festschrift for Leon Chua[END_REF] the Cayley-Hamilton eigenpolynomial associated with the Jacobian of a three-dimensional singularly perturbed system (11) reads:

λ 3 -σ 1 λ 2 + σ 2 λ -σ 3 = 0 (A-21)
where

σ 1 = λ 1 + λ 2 + λ 3 , σ 2 = λ 1 λ 2 + λ 2 λ 3 + λ 1 λ 3 , σ 3 = λ 1 λ 2 λ 3 .
(A-22)

Let's rewrite the eigenpolynomial (A-21) as: a 3 λ 3 + a 2 λ 2 + a 1 λ + a 0 = 0 (a 0 > 0). Routh-Hurwitz' theorem [39, ?] states that the real parts of the eigenvalues of this eigenpolynomial are negative if and only if all the following determinants:

D 1 = a 1 ; D 2 = a 1 a 0 a 3 a 2 = a 1 a 2 -a 0 a 3 (A-23)
are positive. Now, let suppose that the eigenpolynomial (A-21) has one real eigenvalue λ 1 ̸ = 0 and two complex conjugated λ 2,3 = a + ıb (with a ̸ = 0 an b ̸ = 0). So, we have:

σ 1 = λ 1 + 2a, σ 2 = 2aλ 1 + a 2 + b 2 , σ 3 = λ 1 ( a 2 + b 2
) .

(A-24)

The determinant D 2 reads:

D 2 = -2a ( a 2 + b 2 + 2aλ 1 + λ 2 1 ) (A-25)
Moreover, if we consider that the real part of the complex conjugated eigenvalues λ 2,3 depends on a parameter, say µ, we have a = a (µ). Then, determinant D 2 vanishes at the location of the points where the real part a = a (µ). So, it can be used to determine the Hopf-parameter value.

D. Routh-Hurwitz's theorem for 4D systems

According to [START_REF] Tikhonov | On the dependence of solutions of differential equations on a small parameter[END_REF] the Cayley-Hamilton eigenpolynomial associated with the Jacobian of a four-dimensional singularly perturbed system (30) reads:

λ 4 -σ 1 λ 3 + σ 2 λ 2 -σ 3 λ + σ 4 = 0 (A-26)
where (A-27)

σ 1 = λ 1 + λ 2 + λ 3 + λ 4 ,
Let's rewrite the eigenpolynomial (A-26) as: a 4 λ 4 + a 3 λ 3 + a 2 λ 2 + a 1 λ + a 0 = 0 (a 0 > 0). Routh-Hurwitz' theorem [1877[START_REF] Hurwitz | Über die Bedingungen, unter welchen eine Gleichung nur Wurzeln mit negativen reellen Theilen besitzt[END_REF] states that the real parts of the eigenvalues of this eigenpolynomial are negative if and only if all the following determinants:

D 1 = a 1 ; D 2 =
a 1 a 0 a 3 a 2 = a 1 a 2 -a 0 a 3 ; D 3 = a 1 a 0 0 a 3 a 2 a 1 0 a 4 a 3 (A-28) are positive. Now, let suppose that the eigenpolynomial (A-26) has two real eigenvalues λ 1 , λ 2 with λ 1 ̸ = -λ 2 ̸ = 0 and two complex conjugated λ 3,4 = a + ıb (with a ̸ = 0 an b ̸ = 0). So, we have:

σ 1 = 2a + λ 1 + λ 2 , σ 2 = a 2 + b 2 + 2a (λ 1 + λ 2 ) + λ 1 λ 2 , σ 3 = 2aλ 1 λ 2 + ( a 2 + b 2 ) (λ 1 + λ 2 ) , σ 4 = ( a 2 + b 2 ) λ 1 λ 2 .
(A-29)

The determinant D 3 reads:

D 3 = 2a ( a 2 + b 2 + 2aλ 1 + λ 2 1 ) (λ 1 + λ 2 ) ( a 2 + b 2 + 2aλ 2 + λ 2 2 ) (A-30)
Moreover, if we consider that the real part of the complex conjugated eigenvalues λ 2,3 depends on a parameter, say µ, we have a = a (µ). Then, determinant D 3 vanishes at the location of the points where the real part a = a (µ). So, it can be used to determine the Hopf-parameter value.

FIG. 1 .

 1 FIG. 1. "Canard cycle" of the Van der Pol equation, Diener [13, p. 45].

FIG. 3 .

 3 FIG. 3. Piecewise-linear and cubic φ -q characteristic curves for parameter values: a = -2, b = 4 and d = 3.

D 2 FIG. 5 .

 25 FIG. 5. Routh-Hurwitz determinants of system (56). D 1 in blue, D 2 in red and the saddle-node axis γ = 2βc 2 /3 in green for parameter values: ε = 1/10, β = 0.47, a = -2, b = 4, d = 3, c 1 = 280/729 and c 2 = -26/27.

FIG. 6 .

 6 FIG. 6. Numerical "canards solutions" and slow manifold of system (56) for parameter values: ε = 1/10, β = 0.47, γ Duck = 0.3275, a = -2, b = 4, d = 3, c 1 = 280/729 and c 2 = -26/27.

D 2 FIG. 7 .

 27 FIG. 7.Routh-Hurwitz determinants of system (56). D 1 in blue, D 2 in red and the the saddle-node axis γ = 0 in green for parameter values: ε = 1/10, β = γ = α, c 1 = 1/3 and c 2 = -1.

FIG. 9 .

 9 FIG. 9.Region within which the pseudo singular points of fourth-order memristor Chua's circuit (70) are of saddle type.

  and plotted in yellow on Fig.9bbelongs to the cyan region within which the pseudo singular points are of saddle type. On the other hand, this cyan triangular region is limited on the right, at the top of the triangle, by the point of coordinate (α 2 = -2c 2 2c 2 + 3

2 D 1 , D 2 and D 3 FIG. 10 . 1 FIG. 12 .

 21310112 FIG. 10. Routh-Hurwitz determinants of system (70). D 1 in blue, D 2 in red, D 3 in cyan and the saddle-node axis α 2saddle-node = -2c 2 /(3 + 2c 2 ) in green for parameter values: ε = 1/10.1428, β 1 = 0.121, β 2 = 0.0047, c 1 = 0.393781 and c 2 = -0.72357.

σ 2 = λ 1 λ 2 + λ 1 λ 3 + λ 2 λ 3 + λ 1 λ 4 + λ 2 λ 4 + λ 3 λ 4 ,σ 3 = λ 1 λ 2 λ 3 + λ 1 λ 2 λ 4 + λ 1 λ 3 λ 4 + λ 2 λ 3 λ 4 ,σ 4 =

 2233444334444 λ 1 λ 2 λ 3 λ 4 .

  22 a 23 a 24 a 31 a 32 a 33 a 34 a 41 a 42 a 43 a 44

  ϵ, Y 2 , Y U, U 2 ) , Ẏ = 1 + O (X, Y, U, ϵ) ,

	Ż = 1 + O (X, Y, U, ϵ) , ϵ Ż = -( X + U 2 ) + O ( f 2 2 ( ∂ 2 g 1 ∂x 2 2 ∂ 2 g 1 ∂y 2 1 -( ∂ 2 g 1 ϵX, ϵY, ϵU, ϵ 2 , X 2 U, U 3 , XY U ∂x 2 ∂y 1 ) 2 ) + 1 2 f 2 ∂g 1 ∂x 1 ( ∂ 2 g 1 ∂y 2 1 ∂f 1 ∂x 2 -) f 2 3 ( ∂ 2 g 1 ∂x 2 3 ∂ 2 g 1 ∂y 2 1 -( ∂ 2 g 1 ∂x 3 ∂y 1 ) 2 ) + 1 2 f 3 ∂g 1 ∂x 1 ( ∂ 2 g 1 ∂y 2 1 ∂f 1 ∂x 3 -+ f 2 f 3 ã = 1 2 + 1 2 ( ∂ 2 g 1 ∂x 2 ∂x 3 ∂ 2 g 1 ∂y 2 -∂ 2 g 1 ∂x 2 ∂y 1 ) ∂ 2 g 1 , ∂x 3 ∂y 1 1 b = -∂g 1 ∂x 1 ∂f 1 ∂y 1 ,	(A-17) ) ∂f 1 ∂x 2 ∂y 1 ∂ 2 g 1 ∂y 1 ∂ 2 g 1 ) ∂f 1 ∂x 3 ∂y 1 ∂y 1 ,
																(A-20)
	we find:
	ã =A b = A G	( ∂y 1 ∂f 1 ∂x 2 , ∂f 1	-	E G	∂f 1 ∂y 1	)	f 2 + A	(	∂f 1 ∂x 3	-	F G	∂f 1 ∂y 1	)	+ 2Bf 2 2 + 2Cf 2 3 ,	(A-18)
	where
	A = B = C =	1 2 f 3 ∂g 1 ∂x [ 2f 2 1 4 [ ∂ 2 g 1 ∂ 2 g 1 ∂u 2 , ∂ 2 g 1 ∂u 2 ∂ 2 g 1 ∂y∂z ∂z 2 ∂ 2 g 1 ∂u 2 -	+ ( ∂z∂u ∂ 2 g 1 ∂y∂u ) 2 ] ∂ 2 g 1 ∂z∂u ∂ 2 g 1 ,	]	+	1 4	[	∂ 2 g 1 ∂u 2	∂ 2 g 1 ∂y 2 -	(	∂ 2 g 1 ∂y∂u	) 2 ]	,
	D = -	1 f 2	,
	E = -	1 2	∂ 2 g 1 ∂y∂u	,
	F = -	1 2	∂ 2 g 1 ∂z∂u	,
	G = -	1 2	∂ 2 g 1 ∂u 2 .									(A-19)
	Finally, we deduce:

For more details see Robinson[START_REF] Robinson | Nonstandard Analysis[END_REF] and Nelson[START_REF] Nelson | Internal Set Theory: a new approach to nonstandard analysis[END_REF].

Canards in French.1

See Callot et al.[START_REF] Callot | Le problème de la "chasse au canard[END_REF], Benoît et al.[START_REF] Benoît | Chasse au canard[END_REF], Benoît et al.[START_REF] Benoît | Tunnels et entonnoirs[END_REF], Benoît[4] and Ginoux et al.[START_REF] Ginoux | Flow curvature method applied to canard explosion[END_REF].

Canard = false report, from the old-French "vendre un canard moitié" (Sell the half of duck).

This concept has been originally introduced by José Argémi[START_REF] Argémi | Approche qualitative d'un problème de perturbations singulières dans R 4[END_REF]. See Sec. 2.7.

See Fenichel[START_REF] Fenichel | Persistence and smoothness of invariant manifolds for flows[END_REF][START_REF] Fenichel | Asymptotic stability with rate conditions[END_REF][START_REF] Fenichel | Asymptotic stability with rate conditions II[END_REF][START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF], O'Malley[START_REF] O'malley | Introduction to Singular Perturbations[END_REF], Jones[START_REF] Jones | Geometric Singular Perturbation Theory in Dynamical Systems[END_REF] and Kaper[START_REF] Kaper | An Introduction to Geometric Methods and Dynamical Systems Theory for Singular Perturbation Problems[END_REF].

Let's notice that Eq. (73) corresponds exactly to what Itoh and Chua[27, p. 3188] have called in their previous paper on "Memristor Oscillators" the fourth-order memristor-based canonical Chua's circuit equation[START_REF] O'malley | Introduction to Singular Perturbations[END_REF].

In certain applications these functions will be supposed to be C r , r 1.

It represents the approximation of the slow invariant manifold, with an error of O(ε).

The set D is overflowing invariant with respect to (2) when ε = 0. See Kaper[START_REF] Kaper | An Introduction to Geometric Methods and Dynamical Systems Theory for Singular Perturbation Problems[END_REF] and Jones[START_REF] Jones | Geometric Singular Perturbation Theory in Dynamical Systems[END_REF].

In the three-dimensional case det(D ⃗ y ⃗ g) = ∂g 1 /∂y 1 .

This result will be proved below.

Keep in mind that c 2 is generally negative so that the characteristic curve admits a negative slope.

FIG. 11.

Numerical "canards solutions" and critical manifold of system (70) in the (x 1 , x 3 , y 1 ) phase space for parameter values: ε = 1/10.1428, α 2 = 0.1, β 1 = 0.121, β 2 = 0.0047, c 1 = 0.393781 and c 2 = -0.72357.

APPENDIX

Change of coordinates leading to the normal forms of three and fourdimensional singularly perturbed systems with one fast variable are given in the following section.

A. Normal form of 3D singularly perturbed systems with one fast variable

Let's consider the three-dimensional singularly perturbed dynamical system [START_REF] Chua | Memristor -The Missing Circuit Element[END_REF] with k = 2 slow variables and m = 1 fast and let's make the following change of variables:

By taking into account Benoît's generic hypothesis Eqs. [START_REF] Fruchard | Sur le retard à la bifurcation[END_REF][START_REF] Ginoux | Flow curvature method applied to canard explosion[END_REF] and while using Taylor series expansion the system (11) becomes:

Then, let's make the standard polynomial change of variables:

From (A-3) we deduce that:

(A-4)

The time derivative of system (A-3) gives:

Then, multiplying the third equation of (A-5) by (ε/α 2 ) and while replacing in (A-5) ẋ, ẏ and ż by the right-hand-side of system (A-2) leads to:

Since ε/α 2 ≪ 1, the first term of the right-hand-side of the third equation of (A-6) can be neglected. Then, replacing in (A-6) x, y and z by the righthand-side of (A-4) and identifying with the following system in which we have posed: (ε/α 2 ) = ϵ:

we find:

where Then, let's make the standard polynomial change of variables:

From (A-13) we deduce that:

(A-14)

The time derivative of system (A-13) gives:

Then, multiplying the fourth equation of (A-15) by (ε/α 2 ) and while replacing in (A-15) ẋ, ẏ, ż and u by the right-hand-side of system (A-12) leads to:

) , (A-16)

Since ε/α 2 << 1, the two first terms of the right-hand-side of the fourth equation of (A-16) can be neglected. Then, by replacing in (A-16) x, y, z