# Mixed-mode oscillations in slow-fast delayed optoelectronic systems

Jimmi Hervé Talla Mbe, Alain Francis Talla, Géraud Russel Goune Chengui, Aurélien Coillet, Laurent Larger, Paul Woafo, Kouomou Yanne Chembo

► To cite this version:

Jimmi Hervé Talla Mbe, Alain Francis Talla, Géraud Russel Goune Chengui, Aurélien Coillet, Laurent Larger, et al.. Mixed-mode oscillations in slow-fast delayed optoelectronic systems. Physical Review E, 2015, 91 (1), pp.012902 (6). 10.1103/PhysRevE.91.012902. hal-02991659

HAL Id: hal-02991659

https://hal.science/hal-02991659

Submitted on 20 May 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Mixed-mode oscillations in slow-fast delayed optoelectronic systems

Jimmi H. Talla Mbé,¹,²,* Alain F. Talla,¹,² Geraud R. Goune Chengui,¹,² Aurélien Coillet,³,† Laurent Larger,³ Paul Woafo,¹ and Yanne K. Chembo³

¹Laboratory of Modelling and Simulation in Engineering, Biomimetics and Prototypes, University of Yaoundé I, Department of Physics, P.O. Box 812, Yaoundé, Cameroon

²African Center of Excellence for Information and Communication Technologies (CETIC), Polytechnic School of Yaoundé, P.O. Box 8390 Yaoundé, Cameroon

³FEMTO-ST Institute [CNRS UMR6174], Optics Department, 15B Avenue des Montboucons, 25030 Besançon Cedex, France
(Received 9 May 2014; published 5 January 2015)

In this article, we investigate the dynamical behavior of breathers in optoelectronic oscillators from the standpoint of mixed-mode oscillations. In the phase space, these breathers are composite oscillations that are damped to the attractive branches of an invariant manifold. Our study shows that the emergence of breather dynamics is linked to the apparition of inflection points in the phase space, and we develop an analytical framework based on the Liénard reduction form in order to provide an analytical insight into this phenomenology. Our theoretical results are in excellent agreement with experimental measurements.

DOI: 10.1103/PhysRevE.91.012902

PACS number(s): 05.45.Jn, 42.65.Sf

# I. INTRODUCTION

The optoelectronic oscillator (OEO) is one of the most practical representations of an Ikeda system [1], described by a delay-differential equation (DDE). It consists of a closed-loop oscillator characterized by an optical path with local nonlinearity and a linear frequency-filtered electrical path. This paradigmatic system is an ideal benchmark for the investigation of delay-based infinite-dimensional systems, which provide a higher complexity than low-dimensional nonlinear systems ruled by ordinary differential equations (see Ref. [2] and references therein). The OEO mainly finds its numerous applications in two distinct and complementary configurations. In the first one, a narrowband filter is inserted inside the electrical path. In this case, the OEO can be modeled using a complex-valued slowly varying envelope equation: the output signal generally has a single frequency and the main technological application is ultrastable microwave generation [3–6]. In the second configuration, the electrical path is subject to a broadband filter and yields high-dimensional hyperchaos when the feedback gain is large enough. Here the main technological applications are optical communications and information processing [7–14].

The latter case is the one with the most diverse set of complex behaviors. It has for example been shown that, before the state of fully developed hyperchaos, a hybrid state displaying the features of slow-fast dynamics can emerge as the gain is increased beyond the first Hopf bifurcation [7]. These time-domain structures, which are referred to as *breathers*, can be either periodic or chaotic [7]. These breathers are the result of an interaction between a slow variable and a fast other variable. Generally, they appear in one or two successive regular packages that are alternated (up and down).

However, the question of the mechanism leading to the emergence of breathers has still not been fully addressed,

despite some noticeable advances [7,8]. Generally, slow-fast dynamical systems can exhibit bursting oscillations [15–20] or deterministic mixed-mode oscillations (MMOs) if they are at least three-dimensional [21–23], or two-dimensional with noise [24]. MMOs can therefore be considered as trajectories of a dynamical system in which there are alternated oscillations of distinct large and small amplitudes (see Ref. [25] and references therein).

It is remarkable that, by definition, breathers have a structure that is very similar to that of MMOs. For instance, the particular model of chemical reactor reported by Koper [26] displays a time evolution that is very similar to that of breathers. Noticeably, in the time domain, the fast variable has almost two plateaus, one on the top and the second on the bottom. Both plateaus present small-amplitude oscillations (SAOs) that are followed by a large-amplitude oscillation (LAO), thereby generating an alternated sequence [25].

From a dynamical point of view, MMOs arise in systems where there are slow and fast variables, and their phase portrait is characterized by canard cycles [27]. Desroches and Jeffrey [27] investigated the boundary of the singular perturbation parameter that favors the presence of canard explosion in slow-fast dynamical systems. They used the geometric singular perturbation theory (also known as *Fenichel's theory* [28]) to predict the appearance of inflection points that emerge in the repulsive domain of the invariant critical manifold. The reader can for example refer to Refs. [25,27] to gain more insight on theory of MMOs, which are observed in many systems [29], such as chemical reactions [30], biological systems [31,32], and electrical circuits [33,34].

In this article, we use the theory of MMOs to achieve a deeper understanding of the dynamics of breathers in an OEO. The paper is organized as follows. In Sec. II, we present the slow-fast dynamical model of an OEO. Section IV is devoted to the study of the symmetry properties of the breathers. The last section concludes the article.

<footnote>*Corresponding author: jhtallam@yahoo.fr</footnote>

<footnote>†Now at NIST, Boulder, CO 80305, USA.</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. (Color online) Experimental setup of an OEO. The laser used is a DFB telecom laser diode (~1550 nm) with a threshold injection current $I_{th}$ = 15.2 mA. The thermalized 4-km optical fiber induces a delay time of $T_D$ = 19.6 µs. The cutoff frequencies are $f_L$ = 3.1 Hz ($\theta$ = 51.3 ms) and $f_H$ = 480 kHz ($\tau$ = 332 ns). The LiNbO$_3$ Mach-Zehnder modulator is characterized by half-wave voltages $V_{dc}$ = 6.2 V and $V_{rf}$ = 3.9 V.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 2. Temporal evolution of the breathers. The x variable exhibits breathers while the y variable performs a slow periodic oscillation. For this plot, the parameters are $\epsilon$ = 4 × 10⁻⁴, $\nu$ = 2 × 10⁻², $\phi$ = -0.57, and $\beta$ = 2.0.</figcaption></figure>

## II. THE SYSTEM

The system under study is presented in Fig. 1. It generally consists of a Mach-Zehnder (MZ) modulator characterized by radio-frequency (rf) and direct-current (dc) half-wave voltages, $V_{rf}$ and $V_{dc}$, respectively. This modulator is seeded optically by a continuous-wave (cw) semiconductor laser of power P and biased electrically with a constant voltage $V_B$ on the dc electrode. The light at the exit of the modulator passes through an optical delay line of delay $T_D$, and the light intensity is later on converted into electrical current by a photodiode with conversion factor S. The electric signal is then filtered by an electronic filter in which the low- and high-cutoff frequencies are $f_L$ and $f_H$, respectively. The loop is closed by returning the filtered voltage $V(t)$ back into the modulator via the rf electrode, after amplification with overall gain G.

The dynamics of this OEO can be modeled by the following delay integro-differential equation:

$$x + \tau \frac{dx}{dt} + \frac{1}{\theta} \int_{0}^{t} x ds = \beta \{\cos^2[x(t - T_D) + \phi] - \cos^2 \phi\}, \quad (1)$$

where $x(t) = \pi V(t)/2V_{rf}$ is the dimensionless variable standing for the voltage at the rf input of the MZ modulator, while the parameters of the equation are the high cutoff time $\tau = 1/2\pi f_H$, the low cutoff time $\theta = 1/2\pi f_L$, the offset phase $\phi = \pi V_B/2V_{dc}$, and the normalized loop-gain $\beta = \pi SGP/2V_{rf}$.

It is mathematically convenient to rewrite the above equation under the form of a flow of first-order coupled delay differential equations. If we introduce the variable $y = -\frac{1}{\theta} \int_{0}^{t} x(s) ds$ and the dimensionless time $\zeta = t/\theta$, then Eq. (1) can be rewritten as

$$\epsilon \frac{dx}{d\zeta} = y - x + \beta \{\cos^2[x_v + \phi] - \cos^2 \phi\}, \quad (2)$$

$$\frac{dy}{d\zeta} = -x, \quad (3)$$

where $\epsilon = \tau/\theta$ is the frequency cutoff ratio, $\nu = T_D/\theta$ is the normalized delay, and $x_v = x(\zeta - \nu)$ is the delayed variable. Written in the form of Eqs. (2) and (3), the OEO has the structure of a slow-fast dynamical system with the singular perturbation parameter $\epsilon < 1$, one fast variable x involving the fast time scale $\zeta/\epsilon$, and one slow variable y involving the slow time scale $\zeta$.

In the case of a large bandwidth filter, $\tau$ and $\theta$ are several orders of magnitude apart, so that we actually have $\epsilon \ll 1$. It is already known that the OEO displays breathers in this configuration [7]. Breathers generally arise when $\sin 2\phi$ is negative, and they are of two kinds: nonalternated and alternated. Nonalternated breathers are a repetition of one LAO and $s \ge 1$ SAOs that arise in a repetitive way. Alternated breathers are characterized by $k \ge 1$ and $s \ge 1$ SAOs appearing separately on the top and the bottom plateaus, respectively; each group of SAOs is preceded by one large transition. Alternated breathers are plotted in Fig. 2 after simulating Eqs. (2) and (3). The variable x performs k damped oscillations on the top and s damped oscillations on the bottom. In general, the number of oscillations on the top is not always equal to the number of oscillations on the bottom. For example, in the case of Fig. 2, $k \simeq 11$ while $s \simeq 7$. During these oscillations, the variable y displays a slow relaxation evolution. More explanations about this phenomenon are given in the next section.

We also depict in Fig. 3 both experimental and numerical results obtained with $\epsilon = 6.47 \times 10^{-6}$ and $\nu = 3.82 \times 10^{-4}$. The experimental parameter $\beta = (I - I_{th})/(I_H - I_{th})$ is the relative value of the laser injection current (I) calculated over the solitary measured threshold ($I_{th}$) and the closed-loop Hopf threshold at $\phi = -\pi/4$ ($I_H$). The breathers start to appear in the system at $\beta = 1.46$. From that value, they progressively grow as the pump increases. Their inner structure reveals a square-wave oscillation at a period equal to $2T_D = 39.2$ µs. There are approximately 263 oscillations ($k \simeq 263$) on the top against 365 oscillations on the bottom ($s \simeq 365$) in Figs. 3(c) and 3(d) giving a total of $k + s \simeq 628$. When the pump becomes very large, breathers disappear, the square-wave oscillations are no longer present, and, consequently, the
<figure><img src="image_4.png" /><figcaption>FIG. 3. (Color online) Experimental and numerical time traces demonstrating the evolution of breathers as a function of the gain. These results are obtained with $\varepsilon = 6.47 \times 10^{-6}$, $\nu = 3.82 \times 10^{-4}$, and $\phi = -0.97$. For the numerical results, from top to bottom, the values of $\beta$ are 1.5, 2.04, 3.0, and 3.42. For experimental figures, from top to bottom, the values of $\beta$ are 1.48, 1.86, 2.78, and 3.42.</figcaption></figure>

system becomes totally chaotic [see Figs. 3(g) and 3(h)]. The chaotic behavior is confirmed through the evaluation of the maximal Lyapunov exponent which is positive over the time interval of Fig. 3(h) ($\lambda = 0.0166$). For the sake of the clarity of the figures and the phenomenology that we want to highlight in this article, the parameters $\varepsilon$ and $\nu$ hereafter keep the values used in Fig. 2.

The set of Eqs. (2) and (3) can be rewritten under the following general form:

$$\varepsilon \frac{dx}{d\zeta} = y - f(x, x_{\nu}, \beta), \quad (4)$$

$$\frac{dy}{d\zeta} = g(x), \quad (5)$$

and alternatively, it can also be written in the form of

$$\frac{dx}{d\xi} = y - f(x, x_{\mu}, \beta), \quad (6)$$

$$\frac{dy}{d\xi} = \varepsilon g(x), \quad (7)$$

if the dimensionless time is rather defined as $\xi = t/\tau = \zeta/\varepsilon$, while the dimensionless delay is redefined as $\mu = T_D/\tau$. The functions $f$ and $g$ are, respectively, defined as $f = x - \beta\{\cos^2[x_{\nu} + \phi] - \cos^2 \phi\}$ and $g = -x$ in our present case.

The first form [set of Eqs. (4) and (5)] is generally referred to as the slow subsystem (or the reduced system), while the second form [set of Eqs. (6) and (7)] is the fast subsystem (also known as the layer problem). These two subsystems are able to give insight into the full dynamics of the system in the limit of small $\varepsilon$ as we will see later.

If we cancel the delay, i.e., $x_{\nu} \to x$ in the above sets of equations, one obtains for each pair the so-called Liénard system [35]. Recent studies on the Liénard equation have revealed the existence of critical values of $\varepsilon$ for the appearance of canard cycles [27]. Despite the fact that Liénard forms are of course more complex in DDEs than in ordinary differential equations (ODEs), we will see that the theories developed for the case of low-dimensional systems are still useful to understand the phenomenology observed in high-dimensional delayed systems.
### III. CRITICAL MANIFOLD AND PHASE SPACE OF THE OEO EQUATION

According to singular perturbation theory, the critical manifold is the curve in the $(x, y)$ plane defined in the limit of $\epsilon \to 0$. In other terms, following Eqs. (4) and (5), this manifold is explicitly defined as

$$ 
\begin{aligned}
S: y &= f(x,x_v,\beta) \\
&= x - \beta\{\cos^2[x_v + \phi] - \cos^2 \phi\}.
\end{aligned} 
\quad (8) $$

The variables $x$ and $x_v$ are two variables that bring out the same information but at two different times, $\zeta$ and $\zeta - v$. Investigating the dynamics of the system in an infinite-dimensional phase space is a very complicated task: Here, we show that the key manifolds obtained in the limit $v \to 0$ (null delay) still provide an insightful perspective into the mechanisms leading to the emergence of breathers when $v \neq 0$.

In the limit $v \to 0$, the critical (or Liénard) manifold is explicitly defined as

$$ 
\begin{aligned}
S_L: y &= f(x,x,\beta) \\
&= x - \beta\{\cos^2[x + \phi] - \cos^2 \phi\},
\end{aligned} 
\quad (9) $$

where the subscript L stands for Liénard. The fold points of this critical manifold as defined in Ref. [27] are solutions of the equation

$$ \frac{df}{dx} = 1 + \beta \sin(2x + 2\phi) = 0. \quad (10) $$

Restricting the study in the range of $x \in [-\pi/2, \pi/2]$, there exist two fold points, namely:

$$ x_1 = -\frac{1}{2} \arcsin \left[ \frac{1}{\beta} \right] - \phi, \quad (11) $$

$$ x_2 = -\frac{\pi}{2} + \frac{1}{2} \arcsin \left[ \frac{1}{\beta} \right] - \phi. \quad (12) $$

The middle point between $x_1$ and $x_2$ is given by $x_0 = -\pi/4 - \phi$. It should also be noted that $x_1$ and $x_2$ only exist for $\beta > 1$: Hence, we conclude that $df/dx < 0$ if $x \in [x_1, x_2]$, and $df/dx > 0$ outside of this interval. Therefore, based on the definition of attractive and repulsive branches given in Refs. [25,27], the Liénard critical manifold $S_L$ has one repulsive branch located within $x_1$ and $x_2$, that is,

$$ S_L^r := x \in [x_2, x_1], \quad (13) $$

while two attractive branches exist beyond $x_1$ and $x_2$ following

$$ S_L^a := x \in [-\pi/2, x_2] \cup [x_1, \pi/2]. \quad (14) $$

The points $x_1$ and $x_2$ are marked by large dots in the figure. These repulsive and attractive branches are represented in Fig. 4. In the general case, it should be noted that the critical manifold is a function of the bifurcation parameter $\beta$, as well as of the delay. This situation makes the OEO more complex than other usual systems that exhibit MMOs such as the Van der Pol oscillator, as reported in Ref. [25]. In addition, the S-shaped critical manifold shows the possibility of a return mechanism in the nullcline which is favorable for the creation of canard cycles in OEOs.

<figure><img src="image_5.png" /><figcaption>FIG. 4. (Color online) Phase portrait of Liénard form of the OEO with and without delay (black and light gray curves, respectively). The critical manifold from the Liénard without delay is displayed in blue (dark gray). The fold points $x_1$ (on the right) and $x_2$ (on the left) are marked by large dots. The value of the loop gain is $\beta = 1.4$.</figcaption></figure>

Figure 4 presents the nullcline plotted after the numerical simulation of Eqs. (4) and (5), when considering the delay (black curve) and in the absence of the delay ($x_v \to x$, light gray). In fact for low values of $\beta$, starting from an initial condition away from the critical manifold, the OEO is quickly (almost horizontally) pushed towards the attractive branch of the critical manifold $S_L^a$. The dynamics is then dominated by the fast subsystem at this stage. By the time the trajectory passes close to an $\epsilon$ neighborhood of $S_L^a$, the slow dynamics becomes dominant. The orbit now evolves slowly along the attractive branch until it reaches a first fold point; hence it reenters into the fast transition and jumps to the other attractive branch. Once on that branch, the slow motion along the branch occurs until the trajectory reaches the second fold point. Later on, it is reinjected to the first $S_L^a$ and the cycle starts again. This corresponds to a canard cycle, as illustrated by the Liénard form.

<figure><img src="image_6.png" /><figcaption>FIG. 5. (Color online) Phase portrait of breathers (black) and of the Liénard form without delay (light gray). The MMOs are induced by the delay. The value of the loop gain is $\beta = 2.0$. The points in the inset mark some inflection points of the oscillations.</figcaption></figure>
When $\beta$ becomes large (enough to generate breathers), a difference occurs at the moment the orbit reaches an $\varepsilon$ neighborhood of the attractive branch (see Fig. 5). The amplitudes of oscillation become larger, and the slow motion along the attractive branches is first dominated by the damping behavior towards these branches. In this way, the orbit accomplishes $k$ damping SAO before decaying down to $S_L^a$. Once these oscillations are damped, the system moves slowly along $S_L^a$ until it meets the fold point of the branch (domination of the slow subsystem) where it is subjected to the effect of fast dynamics. The orbit is suddenly projected towards another attractive branch where it behaves like on the previous attractive branch; however, it rather accomplishes $s$ damping SAO in this case. It is also noteworthy that $S_L$ is still a good approximation of $S$, and no damping oscillation is observed with the phase space from the Liénard form (see Fig. 5).

#### IV. SYMMETRY OF BREATHERS

The symmetry of breathers depends on the offset phase operating point $\phi$. According to the value of $\phi$, oscillations can be more important on one attractive branch than another. Symmetric oscillations are observed only when the middle

<figure><img src="image_7.png" /><figcaption>FIG. 6. Frontiers between breathers and stability corresponding to $x_1 = 0$ and $x_2 = 0$ are marked by light gray. The vertical dark gray line represents the boundary $\phi = -\pi/4$ for which the breathers are symmetric. (a) Symmetry breathers with $k = s \approx 9$. (b) and (d) Asymmetric breathers. (c) Decay towards the trivial fixed point. All the time traces are plotted with the same value of $\beta = 2.0$. (a) $\phi = -\pi/4$, (b) $\phi = -0.275$, (c) $\phi = -0.12$, (d) $\phi = -1.30$. The point of Fig. 2 is indicated.</figcaption></figure>

point $x_0 = 0$, which corresponds to $\phi = -\pi/4$. In that case, $s = k$ as can be seen in Fig. 6(a). This assertion is identical to the condition of symmetry observed for the square-wave oscillation in the same system [14].

When $\phi$ moves away from $-\pi/4$, the shift has a significant influence for the number of oscillations in the attractive branches. In fact as $x_0$ moves towards the left (this corresponds to $\phi$ moving from $-\pi/4$ to the upper boundary 0), the fold points $x_1$ and $x_2$ decrease in absolute value. The fixed point $x_1$ progressively comes closer to zero, whereas $x_2$ moves away from zero. When $x_1$ coincides with the origin, $\beta$ and $\phi$ satisfy the condition $\beta \sin(-2\phi) = 1$, which is a good approximation of the critical value around which the oscillations arise through a Hopf bifurcation [7,36].

Hence, when $x_1$ moves to zero, it dynamically means that one fold point ($x_1$ in this case) is approaching the trivial fixed point of the system. The result is that more damping oscillations will appear in the attractive branch connected to $x_1$ [see Fig. 6(b)]. From that value of $\phi$, any increase immediately leads the system to the fixed point [see Fig. 6(c)]. On the contrary, if $\phi$ decreases from $-\pi/4$ to the lower boundary $-\pi/2$, it is now $x_2$ that approaches the origin [note that at the origin, $x_2 = 0$ also corresponds to $\beta \sin(-2\phi) = 1$]: In this case, SAO are now mainly localized on the bottom as shown in Fig. 6(d).

The frontiers corresponding to $x_1 = 0$ and $x_2 = 0$ are plotted in Fig. 6, and therefore three domains are created: The domains 1 and 2 correspond to oscillations mainly located on the bottom and the top, respectively, while the domain 3 represents the trivial fixed point. The domains 1 and 2 are separated by the point $\phi = -\pi/4$ (the vertical dark gray line).

#### V. CONCLUSION

In this work, we have proposed an analytical framework for the study of the dynamics of breathers in slow-fast optoelectronic systems. We have demonstrated that the dynamics of these breathers consists of oscillations on attractive branches of a critical manifold. These oscillations are connected one to another by a fast transition that characterizes the repulsive branch of the critical manifold. Using an appropriate Liénard form of the model, we have analyzed the mechanisms leading to the emergence of these breathers in the phase space. We have shown that the bifurcation diagram of an OEO is subdivided into several areas, including the stability basin of the trivial fixed points, canard limit-cycles, breathers which build up with an increase of the number of inner damping oscillations, and finally hyperchaos. Another major conclusion of this work is that delayed systems can also be investigated within the framework of mixed-mode oscillations, and future work will focus on the investigation of other architectures, such as multiple-delay systems.

#### ACKNOWLEDGMENTS

J.H.T.M. acknowledges financial support from AGNES (African-German Network of Excellence in Science). The authors would like to acknowledge financial support from the
ERC project NextPhase, from the project ORA of the *Agence Nationale de la Recherche* (ANR), from the project SHYRO

of the *Centre National d'Etudes Spatiales* (CNES), and from
the Labex ACTION.

[1] K. Ikeda, Opt. Commun. **30**, 257 (1979).

[2] L. Larger, Philos. Trans. R. Soc., A **371**, 20120464 (2013).

[3] X. S. Yao and L. Maleki, J. Opt. Soc. Am. B **13**, 1725 (1996).

[4] Y. K. Chembo, L. Larger, H. Tavernier, R. Bendoula, E. Rubiola, and P. Colet, Opt. Lett. **32**, 2571 (2007).

[5] Y. K. Chembo, L. Larger, and P. Colet, IEEE J. Quantum Electron. **44**, 858 (2008).

[6] Y. K. Chembo, A. Hmima, P. A. Lacourt, L. Larger, and J. M. Dudley, J. Lightwave Technol. **27**, 5160 (2009).

[7] Y. C. Kouomou, P. Colet, L. Larger, and N. Gastaud, Phys. Rev. Lett. **95**, 203903 (2005).

[8] M. Peil, M. Jacquot, Y. K. Chembo, L. Larger, and T. Erneux, Phys. Rev. E **79**, 026208 (2009).

[9] A. B. Cohen, B. Ravoori, T. E. Murphy, and R. Roy, Phys. Rev. Lett. **101**, 154102 (2008).

[10] K. E. Callan, L. Illing, Z. Gao, D. J. Gauthier, and E. Scholl, Phys. Rev. Lett. **104**, 113901 (2010).

[11] B. Ravoori, A. B. Cohen, J. Sun, A. E. Motter, T. E. Murphy, and R. Roy, Phys. Rev. Lett. **107**, 034102 (2011).

[12] L. Weicker, T. Erneux, O. D'Huys, J. Danckaert, M. Jacquot, Y. Chembo, and L. Larger, Phys. Rev. E **86**, 055201(R) (2012).

[13] R. Martinenghi, S. Rybalko, M. Jacquot, Y. K. Chembo, and L. Larger, Phys. Rev. Lett. **108**, 244101 (2012).

[14] L. Weicker, T. Erneux, O. D'Huys, J. Danckaert, M. Jacquot, Y. K. Chembo, and L. Larger, Philos. Trans. R. Soc., A **371**, 20120459 (2013).

[15] S. M. Baer, T. Erneux, and J. Rinzel, SIAM J. Appl. Math. **49**, 55 (1989).

[16] E. M. Izhikevich, Int. J. Bifurcation Chaos Appl. Sci. Eng. **10**, 1171 (2000).

[17] Y. Ji and Q. Bi, Phys. Lett. A **374**, 1434 (2010).

[18] H. Simo and P. Woafo, Mech. Res. Commun. **38**, 537 (2011).

[19] L. T. Abobda and P. Woafo, Commun. Nonlinear Sci. Numer. Simulat. **17**, 3082 (2012).

[20] S. T. Kingni, G. S. M. Ngueuteu, and P. Woafo, Nonlinear Dyn. **76**, 1169 (2014).

[21] M. Wechselberger, SIAM J. Appl. Dyn. Syst. **4**, 101 (2005).

[22] M. Krupa, N. Popovic, N. Kopell, and H. G. Rotstein, Chaos **18**, 015106 (2008).

[23] A. Milik and P. Szmolyan, in *Multiple Time-Scale Dynamical Systems*, The IMA Volumes in Mathematics and its Applications Vol. 122, edited by C. K. R. T. Jones and A. I. Khibnik (Springer, New York, 2001), pp. 117–140.

[24] C. B. Muratov and E. Vanden-Eijnden, Chaos **18**, 015111 (2008).

[25] M. Desroches, J. Guckenheimer, B. Krauskopf, C. Kuehn, H. M. Osinga, and M. Wechselberger, SIAM Rev. **54**, 211 (2012).

[26] M. T. M. Koper, Physica D **80**, 72 (1995).

[27] M. Desroches and M. R. Jeffrey, Proc. R. Soc London, Ser. A **467**, 2404 (2011).

[28] N. Fenichel, J. Differ. Equations **31**, 53 (1979).

[29] M. Brons, T. J. Kaper, and H. G. Rotstein, J. Chem. Phys. **97**, 6191 (1992).

[30] V. Petrov, S. K. Scott, and K. Schowalter, Chaos **18**, 015101 (2008).

[31] T. Hauck and F. W. Schneider, J. Phys. Chem. **97**, 391 (1993).

[32] T. Vo, R. Bertram, J. Tabak, and M. Wechselberger, J. Comput. Neurosci. **28**, 443 (2010).

[33] M. Sekikawa, N. Inaba, T. Yoshinaga, and T. Hikihara, Phys. Lett. A **374**, 3745 (2010).

[34] J. Duhram and J. Moehlis, Chaos **18**, 015110 (2008).

[35] S. H. Strogatz, in *Nonlinear Dynamics and Chaos* (Persus Books, Reading, MA, 1994), p. 210.

[36] G. R. G. Chengui, A. F. Talla, J. H. Talla Mbé, A. Coillet, K. Saleh, L. Larger, P. Woafo, and Y. K. Chembo, J. Opt. Soc. Am. B **31**, 2310 (2014).