application’ similar to the Van der Pol oscillator [69, 70] for continuous-
time fast-slow ODE systems. Sections 5.2-5.3 are more theoretical, and
demonstrate the utility of the theory for analysing discretizations and
Poincaré maps induced by fast-slow ODE systems respectively.

###  5.1 A map-based neural model

We consider a map-based model for neuronal bursting originally due to Chialvo
[5] and considered further in e.g. [33, 68, 71]. The model is given by

$\begin{split}w\mapsto&\ \bar{w}=aw-bv+c,\\\ v\mapsto&\
\bar{v}=v^{2}\exp(w-v)+k,\end{split}$ (75)

where $v$ denotes membrane potential voltage, $w$ is a recovery variable,
$k\geq 0$ models a time-dependent perturbation of the voltage, and
$a\in(0,1)$, $b\in(0,1)$, $c>0$ are parameters relating to the recovery
process; see [5] for a physical interpretation.

Recently in [68], the authors considered a reduced $1-$dimensional model
obtained by considering the recovery variable as a parameter. Although not
mathematically justified, this assumption is motivated by the observation in
[5] that (75) exhibits fast-slow dynamics in particular parameter regimes.
This assumption can be justified in suitable parameter regimes using the
theory developed in Sections 2-3, i.e. using DGSPT. Specifically, the map in
which $w$ is treated as a parameter can be viewed as a layer map obtained in
the singular limit over regions of $(a,b,c)-$parameter space such that

$aw-bv+c=w+\varepsilon g(w,v,\varepsilon),$ (76)

where $0<\varepsilon\ll 1$ and
$g:\mathbb{R}\times\mathbb{R}\times[0,\varepsilon_{0})\to\mathbb{R}$ is
$C^{r}-$smooth. In order for the condition (76) to be satisfied for
$\varepsilon\to 0$ uniformly in $(w,v)-$phase space, we impose the following
model assumption.

######  Assumption 4.

(Restriction to fast-slow parameter regime) The parameters $a\in(0,1)$,
$b\in(0,1)$ and $c>0$ depend smoothly on $\varepsilon$ and satisfy

$\lim_{\varepsilon\to 0}a(\varepsilon)=1^{-},\qquad\lim_{\varepsilon\to
0}b(\varepsilon)=0^{+},\qquad\lim_{\varepsilon\to 0}c(\varepsilon)=0^{+},$

such that in particular we have

$a(\varepsilon)=1-\tilde{a}\varepsilon+O(\varepsilon^{2}),\qquad
b(\varepsilon)=\tilde{b}\varepsilon+O(\varepsilon^{2}),\qquad
c(\varepsilon)=\tilde{c}\varepsilon+O(\varepsilon^{2}),$

for some $\tilde{a}>0$, $\tilde{b}>0$ and $\tilde{c}>0$.

Thus we consider

$\begin{split}w\mapsto&\hskip 2.0pt\bar{w}=w+\varepsilon
g(w,v,\varepsilon),\\\ v\mapsto&\hskip 2.5pt\bar{v}\hskip 1.0pt=v+f(w,v),\\\
\end{split}$ (77)

where

$f(w,v)=-v+v^{2}\exp(w-v)+k,\qquad
g(w,v,\varepsilon)=\tilde{c}-\tilde{b}v-\tilde{a}w+O(\varepsilon).$

The map (77) is in the standard form for fast-slow maps, and can be written in
the general form (10) after setting $N(v,v)=(0,1)^{\textnormal{T}}$ and
$G(w,v,\varepsilon)=(g(w,v,\varepsilon),0)^{\textnormal{T}}$ so that (77)
becomes

$\begin{pmatrix}w\\\ v\end{pmatrix}\mapsto\begin{pmatrix}\bar{w}\\\
\bar{v}\end{pmatrix}=\begin{pmatrix}w\\\ v\end{pmatrix}+\begin{pmatrix}0\\\
1\end{pmatrix}f(w,v)+\varepsilon\begin{pmatrix}g(w,v,\varepsilon)\\\
0\end{pmatrix}.$ (78)

In what follows, we refer to (77) or equivalently (78) as the fast-slow
Chialvo map.

######  Remark 5.1.

The authors in [68] point out that the map (75) is fast-slow for
$0<1-a=b=c:=\varepsilon\ll 1$. This corresponds to the particular case with
$\tilde{a}=\tilde{b}=\tilde{c}=1$ in (77).

######  Remark 5.2.

Further examples of fast-slow maps for which the theory developed herein
applies include the so-called Izhikevich model [28] obtained by Euler
discretization of a continuous-time model in [27], different variants of the
Rulkov model [58, 59, 62], and the Courbage-Nekorkin-Vdovin model appearing in
[50] and later modified in [6]. We refer to the review article [26] and the
many references therein.

####  5.1.1 Layer map and slow manifolds

Due to the explicit separation of ‘time-scales’ in (77), the layer map

$\begin{pmatrix}w\\\ v\end{pmatrix}\mapsto\begin{pmatrix}\bar{w}\\\
\bar{v}\end{pmatrix}=\begin{pmatrix}w\\\ v\end{pmatrix}+\begin{pmatrix}0\\\
1\end{pmatrix}f(w,v)$ (79)

reduces to the $1-$dimensional map $v\mapsto v+f(w,v)$ with parameters
$w\in\mathbb{R}$ and $k\geq 0$ considered in detail in [68]. The critical
manifold $S=\\{(w,v)\in\mathbb{R}^{2}:f(w,v)=0\\}$ can be written as a graph

$S=\left\\{(\varphi_{0}(v),v):v>k\right\\},\qquad\varphi_{0}(v):=v+\ln\left(\frac{v-k}{v^{2}}\right),$

and by Proposition 2.8, the (unique) non-trivial multiplier along $S$ is given
by

$\mu(v)=1+DfN|_{S}=1+D_{v}f|_{S}=\frac{(v-k)(2-v)}{v},\qquad v>k.$ (80)

Non-normally hyperbolic points on $S$ are points $(\varphi_{0}(v),v)$ such
that $|\mu(v)|=1$. Since $\mu(v)\in\mathbb{R}$, this can only occur for fold-
type singularities with $\mu(v)=1$, or flip-type singularities with
$\mu(v)=-1$.

We obtain the following result, which can also (for the most part) be found in
[68]; we simply restate it here in a manner more suited to our formalism.

######  Lemma 5.3.

Fix $k\in(0,3-2\sqrt{2})$. Then the critical manifold decomposes like

$S=S_{-}^{a}\cup F_{-}\cup S_{-}^{r}\cup F_{+}\cup S_{+}^{a}\cup Q\cup
S_{+}^{r},$

where $S_{-}^{a}$, $S_{+}^{a}$ are normally hyperbolic and attracting,
$S_{-}^{r}$, $S_{+}^{r}$ are normally hyperbolic and repelling, $F_{-}$,
$F_{+}$ and $Q$ are respectively fold, fold and supercritical flip points for
the layer map (79).

The fold points are located at $F_{\pm}:(\varphi_{0}(v_{\pm}),v_{\pm})$, where
$v_{\pm}=(1+k\mp\sqrt{k^{2}-6k+1})/2$, the flip point is located at
$Q:(\varphi_{0}(v_{flip}),v_{flip})$, where
$v_{flip}=(3+k+\sqrt{k^{2}-2k+9})/2$, and

$\displaystyle S_{-}^{a}$ $\displaystyle=S\cap\\{v\in(k,v_{-})\\},\qquad$
$\displaystyle S_{-}^{r}=S\cap\\{v\in(v_{-},v_{+})\\},$ $\displaystyle
S_{+}^{a}$ $\displaystyle=S\cap\\{v\in(v_{+},v_{flip})\\},\qquad$
$\displaystyle S_{+}^{r}=S\cap\\{v>v_{flip}\\},$

see Figure 9.

If $k=0$, there is only one fold point $F_{+}:(1,1)$, and the critical
manifold decomposes like $S=S_{-}^{r}\cup F_{+}\cup S_{+}^{a}\cup Q\cup
S_{+}^{r}$.

###### Proof.

The existence of fold and flip bifurcations can be verified directly by
solving $\mu(v)=1$ and $\mu(v)=-1$ respectively for $v$, where $\mu(v)$ is
given by (80), and by checking the relevant non-degeneracy conditions in e.g.
[45, Theorem 4.1] and [45, Theorem 4.3] respectively. Details and explicit
calculations are given in [68, Theorems 2.2-2.3]. ∎

In the following we restrict to $k\in(0,3-2\sqrt{2})$, such that the critical
manifold is ‘$S-$shaped’ as in Figure 9. For further details on the dynamics
in the layer problem for $k=0$ and $k\geq 3-2\sqrt{2}$, we refer again to
[68].

Persistence of (compact submanifolds of) normally hyperbolic branches of $S$
as slow manifolds for $0<\varepsilon\ll 1$ follows by the results in Section
3.

######  Proposition 5.4.

Fix $k\in(0,3-2\sqrt{2})$. For $0<\varepsilon\ll 1$ sufficiently small,
compact submanifolds of $S_{-}^{a}$, $S_{+}^{a}$, $S_{-}^{r}$ and $S_{+}^{r}$
persist as nearby slow manifolds $S_{-,\varepsilon}^{a}$,
$S_{+,\varepsilon}^{a}$, $S_{-,\varepsilon}^{r}$ and $S_{+,\varepsilon}^{r}$
respectively. Each slow manifold is given by the graph of

$\varphi_{\varepsilon}(v)=v+\ln\left(\frac{v-k}{v^{2}}\right)-\varepsilon\left(\frac{\tilde{c}-(\tilde{b}+\tilde{a})v-\tilde{a}\ln\left(\frac{v-k}{v^{2}}\right)}{\mu(v)-1}\right)+O(\varepsilon^{2}),$
(81)

over a suitable compact $v-$interval, e.g.
$S_{+,\varepsilon}^{a}=\\{(\varphi_{\varepsilon}(v),v):v\in[v_{+}+\delta,v_{flip}-\delta]\\}$
for an arbitrarily small but fixed $\delta>0$.

###### Proof.

This follows directly from Theorem 3.2 with notation $(x,y)=(v,w)$, and $N$,
$f$ and $G$ defined as in (78). Substituting the relevant expressions into
(38) gives

$\begin{split}\varphi_{\varepsilon}(v)&=\varphi_{0}(v)-\varepsilon\left[(D_{w}f)^{-1}(DfN)^{-1}DfG\right](\varphi_{0}(v),v)+O(\varepsilon^{2})\\\
&=\varphi_{0}(v)-\varepsilon\left(\frac{g(\varphi_{0}(v),v,0)}{\mu(v)-1}\right)+O(\varepsilon^{2}).\end{split}$

The expression (81) follows after substituting the expressions for $g$ and
$\varphi_{0}(v)$. ∎

Figure 9: Geometry and dynamics of the fast-slow Chialvo map (77) in the
singular limit $\varepsilon=0$, with parameter values $\tilde{a}=1$,
$\tilde{b}=5$, $\tilde{c}=3.5$ and $k=0.035$. The critical manifold $S$ has a
cubic-like profile, with normally hyperbolic and attracting (repelling)
critical manifolds $S^{a}_{\pm}$ ($S^{r}_{\pm}$) shown in blue (red) separated
by non-hyperbolic fold points $F_{\pm}$ (orange) and a supercritical flip
point $Q$ (green); see Lemma 5.3. Iterates of the layer map (79) along
vertical fast fibers (shaded grey) are sketched along with the direction of
iteration, in order to illustrate the fast dynamics. The reduced map, i.e.
(82) truncated at $O(\varepsilon^{2})$, has a unique unstable fixed point
$p_{\ast}\in S_{-}^{r}$, see Proposition 5.8. The direction of ‘slow
iteration’ under the reduced map is indicated with single arrows. Also shown
in cyan is the (numerically computed) curve of period-2 points $S^{2}$
bifurcating from the flip point $Q$, as well as numerically identified flip
bifurcations $Q_{\pm}^{2}$ in the second iterate map $v\mapsto\bar{v}^{2}$. A
period-doubling route to chaos occurs with increasing $w$; see [68] and Remark
5.5.

The slow manifolds $S_{-,\varepsilon}^{a}$, $S_{+,\varepsilon}^{a}$,
$S_{-,\varepsilon}^{r}$ and $S_{+,\varepsilon}^{r}$ identified in Proposition
5.4 have all the properties described in Theorems 3.1, 3.2, 3.3 and 3.4.

######  Remark 5.5.

Classical bifurcation theory implies the existence of a locally quadratic,
normally hyperbolic and attracting critical manifold $S^{2}$ for the second
iterate (layer) map $v\mapsto\bar{v}^{2}$. The manifold $S^{2}$ intersects $S$
transversally at the supercritical flip point $Q$, see Figure 9. Dynamically,
$S^{2}$ forms a manifold of stable period-2 cycles for the layer map (79). In
fact there is a period-doubling cascade for increasing $w$ which leads to
chaotic dynamics in the layer map. We do not consider chaotic dynamics in
detail in this work, but refer to [68] for details on chaotic dynamics in the
layer map (79) and [2] for details on the slow passage through a period-
doubling cascade.

####  5.1.2 Reduced map and dynamics on slow manifolds

To leading order in $\varepsilon$, the dynamics on slow manifolds
$S_{-,\varepsilon}^{a}$, $S_{+,\varepsilon}^{a}$, $S_{-,\varepsilon}^{r}$ and
$S_{+,\varepsilon}^{r}$ is determined by the reduced map (24). This is
obtained directly using Proposition 2.14.

######  Lemma 5.6.

Slow iteration along $S_{-,\varepsilon}^{a}$, $S_{+,\varepsilon}^{a}$,
$S_{-,\varepsilon}^{r}$ and $S_{+,\varepsilon}^{r}$ is governed by the map

$\begin{pmatrix}w\\\ v\end{pmatrix}\mapsto\begin{pmatrix}\bar{w}\\\
\bar{v}\end{pmatrix}=\begin{pmatrix}\varphi_{0}(v)\\\
v\end{pmatrix}+\varepsilon\begin{pmatrix}1\\\
-\tfrac{v-k}{\mu(v)-1}\end{pmatrix}\left(\tilde{c}-(\tilde{b}+\tilde{a})v-\tilde{a}\ln\left(\frac{v-k}{v^{2}}\right)\right)+O(\varepsilon^{2}),$
(82)

for all $(w,v)=(\varphi_{0}(v),v)\in S\setminus(F_{-}\cup F_{+})$ restricted
to the corresponding (compact) $v-$interval. By invariance, this can be
represented by the $1-$dimensional map

$v\mapsto\bar{v}=v-\varepsilon\left(\frac{v-k}{\mu(v)-1}\right)\left(\tilde{c}-(\tilde{b}+\tilde{a})v-\tilde{a}\ln\left(\frac{v-k}{v^{2}}\right)\right)+O(\varepsilon^{2}).$
(83)

###### Proof.

Let $S_{n}:=S\setminus(F_{-}\cup F_{+})$. The form of the map (82) is obtained
via equation (21) in Proposition 2.14, noting that the projection operator
$\Pi_{\mathcal{N}}^{S_{n}}$ is given by

$\Pi_{\mathcal{N}}^{S_{n}}=I_{2}-N(DfN)^{-1}Df|_{S_{n}}=\begin{pmatrix}0&1\\\
0&-(D_{v}f)^{-1}D_{w}f\end{pmatrix}\bigg{|}_{S_{n}}=\begin{pmatrix}0&1\\\
0&-\tfrac{v-k}{\mu(v)-1}\end{pmatrix}\bigg{|}_{S_{n}},$

where we used the fact that $\mu(v)=1+D_{v}f|_{S_{n}}$. ∎

Fixed points on $S_{-,\varepsilon}^{a}$, $S_{+,\varepsilon}^{a}$,
$S_{-,\varepsilon}^{r}$ and $S_{+,\varepsilon}^{r}$ occur for values
$v_{\ast,\varepsilon}=v_{\ast}+O(\varepsilon)$, where $v_{\ast}$ satisfies

$g(\varphi_{0}(v_{\ast}),v_{\ast},0)=\tilde{c}-(\tilde{b}+\tilde{a})v_{\ast}-\tilde{a}\ln\left(\frac{v_{\ast}-k}{v_{\ast}^{2}}\right)=0.$
(84)

The intermediate value theorem guarantees that equation (84) has at least one
solution $v_{\ast}\in(k,\infty)$, since
$g(\varphi_{0}(v),v,0)\sim-\tilde{a}\ln((v-k)/v^{2})>0$ as $v\to k^{+}$ and
$g(\varphi_{0}(v),v,0)\sim-(\tilde{b}+\tilde{a})v<0$ as $v\to\infty$. For
simplicity we restrict to the case in which $v_{\ast}$ is the only solution
and, correspondingly, the reduced map (82) has precisely one fixed point.

######  Assumption 5.

(Unique equilibrium on $S$) The parameters $\tilde{a}>0$, $\tilde{b}>0$,
$\tilde{c}>0$ and $k\in(0,3-2\sqrt{2})$ are chosen such that the function
$g(\varphi_{0}(v),v,0)$ is strictly decreasing, i.e.
$D_{v}g(\varphi_{0}(v),v,0)<0$ for all $v>k$. This is true for parameter
values satisfying

$-2k\tilde{a}+(\tilde{a}+\tilde{a}k+\tilde{b}k)v-(\tilde{a}+\tilde{b})v^{2}<0$
(85)

for all $v>k$.

Assumption 5 is satisfied for the parameter values in Figure 9, for example,
on the condition that $-0.07+1.21v-6v^{2}<0$ for all $v>k=0.035$. This is
indeed the case, since this expression is negative at e.g. $v=1$, and the
discriminant $\Delta\approx-0.22<0$.

######  Remark 5.7.

For other parameter choices the graph of $g(\varphi_{0}(v),v,0)$ has a ‘cubic’
profile, with a unique local minimum and a unique local maximum. In this case
equation (84) can have either $1$, $2$ or $3$ solutions depending on whether
the local minimum is positive, zero or negative respectively. A complete
description of the partitioning of the phase portrait in
$(\tilde{a},\tilde{b},\tilde{c},k)-$space is beyond the scope of this article.

The following result classifies the stability of the perturbed fixed point
$p_{\ast,\varepsilon}$ depending on whether
$p_{\ast}:(\varphi_{0}(v_{\ast}),v_{\ast})$ lies on $S_{-}^{a}$, $S_{+}^{a}$,
$S_{-}^{r}$ or $S_{+}^{r}$.

######  Proposition 5.8.

Consider the fast-slow Chialvo map (77) with parameter values in the set
defined by Assumption 5. The following assertions hold for all
$0<\varepsilon\ll 1$ sufficiently small:

  1. (i)

For parameter values $\tilde{a},\tilde{b},\tilde{c}$ and $k$ such that
$p_{\ast}\in S_{-}^{a}$, $S_{+}^{a}$ or $S_{+}^{r}$, there exists a unique
fixed point $p_{\ast,\varepsilon}\in S_{-,\varepsilon}^{a}$,
$S_{+,\varepsilon}^{a}$ or $S_{+,\varepsilon}^{r}$ respectively, which is
$O(\varepsilon)-$close to $p_{\ast}$ and asymptotically stable within
$S_{-,\varepsilon}^{a}$, $S_{+,\varepsilon}^{a}$ or $S_{+,\varepsilon}^{r}$
respectively.

  2. (ii)

For parameter values $\tilde{a},\tilde{b},\tilde{c}$ and $k$ such that
$p_{\ast}\in S_{-}^{r}$, there exists a unique fixed point
$p_{\ast,\varepsilon}\in S_{-,\varepsilon}^{r}$ which is
$O(\varepsilon)-$close to $p_{\ast}$ and asymptotically unstable within
$S_{-,\varepsilon}^{r}$.

###### Proof.

First we show that if $p_{\ast}$ is a normally hyperbolic point on $S$ then it
persists as an $O(\varepsilon)-$close fixed point of the fast-slow Chialvo map
(77) for $0<\varepsilon\ll 1$. By Theorem 3.1, $p_{\ast}$ perturbs to a point
$p_{\ast,\varepsilon}$ on one of the slow manifolds $S_{-,\varepsilon}^{a}$,
$S_{+,\varepsilon}^{a}$, $S_{-,\varepsilon}^{r}$ or $S_{+,\varepsilon}^{r}$.
It therefore suffices to prove the result for the map (83). Since
$g(\varphi_{0}(v_{\ast}),v_{\ast},0)=0$ and by Assumption 5
$D_{v}g(\varphi_{0}(v_{\ast}),v_{\ast},0)<0$, the implicit function theorem
implies persistence of $p_{\ast}$ as an $O(\varepsilon)-$close fixed point
$p_{\ast,\varepsilon}$ as required.

It remains to prove the stability claims in statements (i)-(ii). The Jacobian
associated to the right-hand-side of (83) has a unique multiplier

$\mu(v_{\ast,\varepsilon})=1-\varepsilon\left(\frac{v_{\ast}-k}{\mu(v_{\ast})-1}\right)D_{v}g(\varphi_{0}(v_{\ast}),v_{\ast},0)+O(\varepsilon^{2}),$

where $v_{\ast,\varepsilon}=v_{\ast}+O(\varepsilon)$ is the $v-$coordinate of
$p_{\ast,\varepsilon}$. Since $v_{\ast}-k>0$ and by Assumption 5
$D_{v}g(\varphi_{0}(v_{\ast}),v_{\ast},0)<0$, the stability claims in
statements (i)-(ii) follow from the fact that $\mu(v)-1<0$ on $S_{-}^{a}\cup
S_{+}^{a}\cup S_{+}^{r}$ and $\mu(v)-1>0$ on $S_{-}^{r}$. ∎

We note that each possibility in Proposition 5.8 is realised on open parameter
sets. For the particular choice of parameters in Figure 9, a unique
asymptotically unstable fixed point $p_{\ast}\in S_{-}^{r}$ perturbs to a
unique asymptotically unstable fixed point on $S_{-,\varepsilon}^{r}$ for
$0<\varepsilon\ll 1$.

####  5.1.3 Global dynamics for $0<\varepsilon\ll 1$

Using the preceding analysis, we can provide a partial geometric description
of the global dynamics of the fast-slow Chialvo map in suitable regions of
phase space, specifically, away from non-normally hyperbolic points on $S$. We
consider four possibilities, distinguished on the ‘singular level’ according
to (i) which branch the fixed point $p_{\ast}$ is on, and (ii) the relative
positioning of the lower fold point $F_{-}$ and the flip point $Q$. We
restrict to the case of a single fixed point on $S$ in accordance with
Assumption 5 throughout. Numerical computations are carried out using the
software package MatContM [46].

######  Remark 5.9.

A full classification of bursting mechanisms in the fast-slow Chialvo map (77)
is beyond the scope of this work. We refer to [28] for more on the
classification of bursting mechanisms in two-dimensional fast-slow maps.

(a) Case I: Excitability.

(b) Case II: Relaxation.

(c) Case III: Non-chaotic bursting.

(d) Case IV: Chaotic bursting.

Figure 10: Global dynamics of the fast-slow Chialvo map (77) for
$0<\varepsilon\ll 1$. Singular structure for $\varepsilon=0$ including
critical manifolds $S_{\pm}^{a}$, $S_{\pm}^{r}$, the fold, flip and
equilibrium points $F_{\pm}$, $Q$ and $p_{\ast}$ respectively are shown
overlaid with numerically computed iterates of (77) with initial condition
$(w,v)=(1/4,2)$, $\varepsilon=10^{-3}$, $\tilde{a}=1$, $\tilde{b}=5$ and
varying values of $\tilde{c},k$. The ‘slow nullcline’ $\\{g(w,v,0)=0\\}$ is
also shown in shaded grey. Dynamics away from the non-normally hyperbolic
points $F_{\pm}$ and $Q$ is described by the theory in Sections 2 and 3. We
identify four distinct behaviours depending on $(\tilde{c},k)$, which are
described in the text. (a) Case I: Excitability, with
$(\tilde{c},k)=(7,0.07)$. (b) Case II: Relaxation, with
$(\tilde{c},k)=(3.5,0.07)$. (c) Case III: Non-chaotic bursting, with
$(\tilde{c},k)=(3.5,0.035)$. (d) Case IV: Chaotic bursting, with
$(\tilde{c},k)=(3.5,0.02)$.

#### Case I: Excitability

Excitable dynamics occurs for parameters such that $p_{\ast}\in S_{+}^{a}$
with $w_{\ast}=\varphi_{0}(v_{\ast})<w_{-}$ so that it lies to the left of the
lower fold point $F_{-}$, and $w_{flip}>w_{-}$ so that the flip point $Q$ lies
to the right of $F_{-}$. Figure 10(a) shows the singular dynamics for a
particular choice of parameter values satisfying these conditions, overlaid
with iterates of the map (77) for $\varepsilon=10^{-3}$.

Starting from the point $(w,v)=(1/4,2)$, subsequent iterates are exponentially
attracted to the attracting slow manifold $S^{a}_{-,\varepsilon}$ along stable
fibers according to Theorem 3.4. After reaching a neighbourhood of
$S_{-,\varepsilon}^{a}$, iterates are well approximated by the map (82) (up to
an error of $O(\varepsilon^{2})$), and move to the right for an
$O(\varepsilon^{-1})$ but finite number of iterates before reaching a
neighbourhood of $F_{-}$. Unfortunately the dynamics near the fold point near
$F_{-}$ are not covered by the normally hyperbolic theory in Sections 2-3. We
therefore omit this part in our description; see however Remark 5.10 below.
After leaving a neighbourhood of $F_{-}$, iterates are exponentially attracted
along stable fibers to $S^{a}_{+,\varepsilon}$, before moving to the left as
governed by the map (82). After $O(\varepsilon^{-1})$ but finitely many
iterations, iterates reach (and remain thereafter in) a neighbourhood of the
asymptotically stable fixed point $p_{\ast,\varepsilon}\in
S_{+,\varepsilon}^{a}$ which perturbs from $p_{\ast}\in S_{+}^{a}$ as
described by Proposition 5.8.

######  Remark 5.10.

In particular cases, fold bifurcations in maps arising after discretization of
ODE systems with a (continuous-time) fold bifurcation have been treated in
detail using geometric techniques in e.g. [54, 55]. In general, however, a
detailed geometric description of the dynamics near the discrete fold
bifurcation remains a topic of future work.

#### Case II: Relaxation

Relaxation-type dynamics occurs for parameter values such that $p_{\ast}\in
S^{r}_{-}$ and $w_{flip}>w_{-}$, for which flip point $Q$ lies to the right of
the lower fold point $F_{-}$. Figure 10(b) shows the singular dynamics for
particular choice of parameter values, overlaid with iterates of the map (77)
for $\varepsilon=10^{-3}$.

As before, we take an initial condition $(w,v)=(1/4,2)$. The geometric
description of the dynamics is analogous to the excitable case I, up to the
point where iterates reach a neighbourhood of $S^{a}_{+,\varepsilon}$. In the
relaxation case, iterates continues to track iterates of (82) until they reach
a neighbourhood of the upper fold $F^{+}$. As with the lower fold $F_{-}$, the
normally hyperbolic theory of Sections 2-3 breaks down near $F_{+}$, and we
omit a proper treatment of the dynamics here. Normally hyperbolic theory
applies once more after iterates leave a neighbourhood of $F_{+}$, after which
iterates are exponentially attracted along stable fibers to
$S_{-,\varepsilon}$. This provides the basic mechanism for relaxation-type
dynamics observed in Figure 10(b).

#### Case III: Non-chaotic bursting

Here we consider parameter values such that $p_{\ast}\in S^{r}_{-}$,
$w_{flip}\in(w_{+},w_{-})$ and $|w_{flip}-w_{-}|<\delta$ is sufficiently small
that the flip points $Q^{2}_{\pm}$ in the second iterate layer map
$v\mapsto\bar{v}^{2}$ lie to the right of $F_{-}$ as in Figure 9. Figure 10(c)
shows the singular dynamics for the same choice of parameter values as Figure
9, overlaid with iterates of the map (77) for $\varepsilon=10^{-3}$.

Starting from $(w,v)=(1/4,2)$, iterates approach $S_{-,\varepsilon}^{a}$ and
subsequently follow iterates of the map (82) up to a neighbourhood of $F_{-}$,
similarly to cases I and II. After leaving a neighbourhood of $F_{-}$,
iterates are exponentially attracted to an attracting slow manifold
$S_{+,\varepsilon}^{a,2}$ which perturbs from (a suitable submanifold of) the
normally hyperbolic and attracting manifold of period-2 points $S^{2}$ in
Figure 9; the existence and properties of the slow manifold
$S_{+,\varepsilon}^{a,2}$ and its stable manifold follow after applying the
results in Section 3 to the second iterate map
$(w,v)^{\textnormal{T}}\mapsto(\bar{w}^{2},\bar{v}^{2})^{\textnormal{T}}$.
Iterates move leftward along $S_{+,\varepsilon}^{a,2}$ in a ‘flip-like’
orientation-reversing manner about the repelling slow manifold
$S_{+,\varepsilon}^{r}$. They subsequently traverse a neighbourhood of the
flip point $Q$, before being attracted to $S_{+,\varepsilon}^{a}$ after a
delay effect which is described in detail in [2]. Finally, iterates track
iterates of the map (82) along $S_{+,\varepsilon}^{a}$ towards to fold
$F_{+}$, and return to a neighbourhood of $S_{-,\varepsilon}^{a}$ after
leaving a neighbourhood of $F_{+}$. Assuming a simple jump-type dynamics at
$F_{\pm}$ (as is indeed observed numerically), the geometric sequence repeats,
providing a mechanism for non-chaotic bursting.

######  Remark 5.11.

Similar (though more complicated) non-chaotic dynamics is expected as long as
$w_{flip}$ and $w_{-}$ are sufficiently close, i.e. $|w_{flip}-w_{-}|<\delta$
for $\delta>0$ sufficiently small.

#### Case IV: Chaotic bursting

Lastly, we consider parameter values such that $p_{\ast}\in S^{r}_{-}$,
$w_{flip}\in(w_{+},w_{-})$ and $|w_{flip}-w_{-}|>\delta$, where the latter
condition implies that the layer map is chaotic for $w-$values close to
$w_{-}$. Figure 10(d) shows the singular dynamics for particular choice of
parameter values satisfying these conditions, overlaid with iterates of the
map (77) for $\varepsilon=10^{-3}$.

Starting again from $(w,v)=(1/4,2)$, iterates are initially attracted to
$S_{-,\varepsilon}^{a}$ and subsequently track iterates of the map (82)
towards the fold $F_{-}$, similarly to cases I, II and III. After leaving a
neighbourhood of the fold, complicated dynamics corresponding to a slow
(leftward) drift through part of a period-doubling cascade occurs. The details
of this transition are beyond the scope of this work, however the reader is
referred to [68] for details on the period-doubling cascade in the
corresponding layer problem, and to [2] for details on the slow drift through
a period-doubling cascade in a particular ‘normal form’. Once iterates are
close enough to $w_{flip}$, the dynamics becomes non-chaotic and orientation-
reversing as in case III, and eventually approach a neighbourhood of
$S_{+,\varepsilon}^{a}$ after a delay effect (again described by [2]). The
subsequent tracking of iterates along $S_{+,\varepsilon}^{a}$ and the observed
return to $S_{-,\varepsilon}^{a}$ is similar to that described for cases II
and III. Assuming a simple-jump like behaviour near the folds $F_{\pm}$, this
geometric sequence repeats itself and provides a mechanism for chaotic
bursting behaviour.

###  5.2 Euler-discretized fast-slow ODEs in non-standard form

We consider the existence and properties of invariant manifolds in Euler-
discretizations of singularly perturbed ODEs in the general non-standard form

$z^{\prime}=N(z)f(z)+\varepsilon G(z,\varepsilon),\qquad
z\in\mathbb{R}^{n},\qquad 0<\varepsilon\ll 1,$ (86)

where the right-hand-side is $C^{r}-$smooth in $(z,\varepsilon)$ and assumed
to have a normally hyperbolic critical manifold for $\varepsilon=0$. More
precisely, we impose the following assumptions:

######  Assumption 6.

(Restriction to fast-slow ODEs) The ODE (86) satisfies the geometric
definition of a $C^{r}-$smooth singularly perturbed ODE in e.g. [17, 72], i.e.
the level set

$C_{ode}:=\left\\{z\in\mathbb{R}^{n}:N(z)f(z)=O_{n}\right\\}$

contains a $(k<n)-$dimensional regularly embedded critical manifold

$S_{ode}:=\left\\{z\in\mathbb{R}^{n}:f(z)=O_{n-k}\right\\}.$

We assume for simplicity that $S\subseteq C\subset\mathbb{R}^{n}$ is
connected.

######  Assumption 7.

(Regularity of the matrix $N(z)$) We assume that the $n\times(n-k)$ matrix
$N(z)$ has full column rank for all $z\in S_{ode}$, and that equilibria of
$N(z)f(z)$ for all $z\notin S_{ode}$, if they exist, are isolated.

Assumptions 6-7 are the defining assumptions for fast-slow ODEs in [72], which
are directly analogous to the defining Assumptions 1-2 for fast-slow maps in
Section 2.

In the following we are concerned with the dynamics of maps induced by Euler
discretization of general fast-slow ODEs (86) under Assumptions 6-7.
Specifically, Euler discretization in time yields the map

$z\mapsto\bar{z}=z+hN(z)f(z)+\varepsilon hG(z,\varepsilon),$ (87)

where $h>0$ is the associated step-size. It is straightforward to verify that
the map (87) has critical manifold $S=\\{z\in\mathbb{R}^{n}:f(z)=O_{n-k}\\}$
and satisfies Assumptions 1-2; this follows directly from Assumptions 6-7 on
the ODE (86). Hence, the formalism of Section 2 applies.

######  Remark 5.12.

All results obtained for the map (87) extend and are similar in nature to
their pre-existing counterparts in the context of Euler-discretized fast-slow
ODEs in standard form

$\begin{split}x^{\prime}&=\varepsilon\tilde{g}(x,y,\varepsilon),\\\
y^{\prime}&=\tilde{f}(x,y,\varepsilon),\end{split}\qquad(x,y)\in\mathbb{R}^{k}\times\mathbb{R}^{n-k},\qquad
0<\varepsilon\ll 1,$ (88)

since this system can be written in the general form (86) by setting
$z=(x,y)^{\textnormal{T}}$ and

$N(z)=\begin{pmatrix}O_{k,n-k}\\\ I_{n-k}\end{pmatrix},\qquad
f(z)=\tilde{f}(x,y,0),\qquad
G(z,\varepsilon)=\begin{pmatrix}\tilde{g}(x,y,\varepsilon)\\\
\tilde{f}_{rem}(x,y,\varepsilon)\end{pmatrix},$

where
$\tilde{f}_{rem}(x,y,\varepsilon)=\tilde{f}(x,y,\varepsilon)-\tilde{f}(x,y,0)$.
See [52, 53] and the references therein for more on existing results for
discretizations of fast-slow ODEs in standard form (88).

#### Layer map

The layer map is given by

$z\mapsto\bar{z}=z+hN(z)f(z).$ (89)

By Proposition 2.8 and Corollary 2.11, the $n-k$ non-trivial multipliers
$\mu_{j}(z)$ along the critical manifold $S$ are given by

$\mu_{j}(z)=1+\lambda_{j}(z),\qquad i=1,\ldots,n-k,$

where $\lambda_{j}(z)$ are eigenvalues of the matrix

$hDf(z)N(z),\qquad z\in S.$

For the ODE, non-trivial eigenvalues $\lambda_{j,ode}(z)$ along $S_{ode}$ are
given by eigenvalues to the matrix $Df(z)N(z)$. It follows that

$\mu_{j}(z)=1+\lambda_{j}(z)=1+h\lambda_{j,ode}(z).$ (90)

We obtain the following stability conditions for $S$, depending on the real
part of the ODE eigenvalues along $S_{ode}$ and the step-size $h$ associated
to the discretization.

######  Lemma 5.13.

Consider the layer map (89) with step-size $h>0$. The critical manifold $S$ is
normally hyperbolic at $z\in S$ if

$2\textup{Re}\,\lambda_{j,ode}(z)\neq-h|\lambda_{j,ode}(z)|^{2},$

for all $j=1,\ldots,n-k$. A normally point $z\in S$ is attracting if

$2\textup{Re}\,\lambda_{j,ode}(z)<-h|\lambda_{j,ode}(z)|^{2},$

for all $j=1,\ldots,n-k$, repelling if

$2\textup{Re}\,\lambda_{j,ode}(z)>-h|\lambda_{j,ode}(z)|^{2},$

for all $j=1,\ldots,n-k$, and saddle-type otherwise.

###### Proof.

This follows directly from equation (90) and Corollary 2.11. ∎

It follows that normal hyperbolicity breaks down at fold-type singularities
$z\in S$ where $\lambda_{j,ode}(z)=0$ for some $j=1,\ldots,n-k$, independently
of $h$. On the other hand, Hopf-type non-normally hyperbolic singularities in
the ODE have eigenvalues such that $\textup{Im}\,\lambda_{j,ode}(z)\neq 0$,
and their location and may be shifted in by $O(h)$ distances in the
discretization. The discretization also introduces additional singularities
that are not present in the ODE, with multipliers that cross through the left
hand side of the unit circle in $\\{z\in\mathbb{C}:\textup{Im}\,z\leq 0\\}$,
however the location of the corresponding singularities depends on $h$, and
tend to infinity in the limit $h\to 0^{+}$.

#### Reduced map

We now consider the reduced map for the Euler discretization (87). If we
assume that $S_{ode}$ is normally hyperbolic, it follows from Lemma 5.13 that
for sufficiently small step-sizes $h>0$, compact submanifolds of $S$ are also
normally hyperbolic since any non-hyperbolic points on $S$ induced by the
discretization are pushed to infinity as $h\to 0^{+}$. Therefore, the oblique
projection operator in (20) and (22) is well-defined on compact submanifolds
of $S$ for sufficiently small $h>0$. In fact, the projection is invariant with
respect to the discretization in the sense that

$\Pi_{\mathcal{N}}^{S}=I_{n}-hNh^{-1}(DfN)^{-1}Df|_{S}=I_{n}-N(DfN)^{-1}Df|_{S}=\Pi_{\mathcal{N},ode}^{S},$
(91)

where $\Pi_{\mathcal{N},ode}^{S}$ denotes the projection operator used to
define the reduced problem for the ODE (86), see [72, Def. 3.8]. It follows by
Proposition 2.14 that the reduced map is

$z\mapsto\bar{z}=z+\varepsilon h\Pi_{\mathcal{N}}^{S}G(z,0)|_{S}=z+\varepsilon
h\Pi_{\mathcal{N},ode}^{S}G(z,0)|_{S},$

with $\Pi_{\mathcal{N}}^{S}=\Pi_{\mathcal{N},ode}^{S}$ given by (91). As noted
already in Remark 2.18, in this case it is possible to rescale the step-size
via $h=\varepsilon^{-1}\tilde{h}$ in order to obtain a well-defined reduced
map on $S$ in the limit $\varepsilon\to 0$. Explicitly, we obtain

$z\mapsto\bar{z}=z+\varepsilon
h\Pi_{\mathcal{N}}^{S}G(z,0)|_{S}=z+\tilde{h}\Pi_{\mathcal{N},ode}^{S}G(z,0)|_{S},$

which is $\varepsilon-$independent when formulated in terms of the rescaled
step-size parameter $\tilde{h}=\varepsilon h$.

#### Dynamics for $0<\varepsilon\ll 1$

The results of Section 3 apply to the Euler discretization (87) under
conditions of normally hyperbolicity. Moreover, Lemma 5.13 shows that the
normally hyperbolicity of the critical manifold $S$ after discretization is
entirely determined by the non-trivial eigenvalues along the critical manifold
of the ODE (86) and the step-size $h>0$. This allows for the statements on the
existence and properties of slow manifolds and corresponding stable/unstable
manifolds in the discretized system, purely in terms of ODE properties and the
step-size parameter $h$.

It will be helpful for the statement of our results to introduce additional
notation for spectral bounds in the ODE. Given a $k-$dimensional normally
hyperbolic critical manifold $S_{ode}$, denote the $n_{a}\leq n-k$ non-trivial
eigenvalues such that $\textrm{Re}\,\lambda_{j,ode}(z)<0$ by
$\lambda_{a,j,ode}(z)$, and the $n_{r}=n-k-n_{a}$ non-trivial eigenvalues such
that $\textrm{Re}\,\lambda_{j,ode}(z)>0$ by $\lambda_{r,j,ode}(z)$. We define
the corresponding spectral bounds by

$\nu_{A,ode}:=\sup_{z\in
S_{ode},j=1,\ldots,n_{a}}|\textrm{Re}\,\lambda_{a,j,ode}(z)|,\qquad\nu_{R,ode}:=\inf_{p\in
S_{ode},j=1,\ldots,n_{r}}|\textrm{Re}\,\lambda_{r,j,ode}(z)|.$ (92)

Our main results for the Euler discretization (87) are as follows.

######  Theorem 5.14.

(Existence and characterisation of slow and stable/unstable manifolds in Euler
discretized ODEs) Assume that the fast-slow ODE system (86) satisfies
Assumptions 6-7, with compact, connected and normally hyperbolic critical
manifold $S_{ode}$. Then for the Euler discretization (87) with
$\varepsilon=0$ we have $S=S_{ode}$ and coincident transverse linear fiber
bundles $\mathcal{N}=\mathcal{N}_{ode}$. Moreover, there exists an
$\varepsilon_{0}>0$ and a $h_{0}>0$ such that for all
$\varepsilon\in(0,\varepsilon_{0})$ and $h\in(0,h_{0})$ the following
assertions are true:

  1. (i)

$S$ and $S_{ode}$ perturb to nearby $C^{r}-$smooth locally invariant slow
manifolds $S_{\varepsilon,h}$ and $S_{ode,\varepsilon}$ respectively, such
that

$\textup{dist}(S_{\varepsilon,h},S_{ode,\varepsilon})=O(\varepsilon^{2}h).$
(93)

For each fixed $h$ the manifold $S_{\varepsilon,h}$ has all the properties
described in Theorem 3.1, and the manifold $S_{ode,\varepsilon}$ is a Fenichel
slow manifold.

  2. (ii)

$S_{\varepsilon,h}$ and $S_{ode,\varepsilon}$ have local graph representations

$S_{\varepsilon,h}=\left\\{(x,\varphi_{\varepsilon,h}(x)):x\in\mathcal{K}\right\\},\qquad
S_{ode,\varepsilon}=\left\\{(x,\varphi_{ode,\varepsilon}(x)):x\in\mathcal{K}\right\\},$

where

$\begin{split}\varphi_{\varepsilon,h}(x)&=\varphi_{0}(x)-\varepsilon(D_{y}f)^{-1}(DfN)^{-1}DfG|_{S}+O(\varepsilon^{2}h),\\\
\varphi_{ode,\varepsilon}(x)&=\varphi_{0}(x)-\varepsilon(D_{y}f)^{-1}(DfN)^{-1}DfG|_{S_{ode}}+O(\varepsilon^{2}).\end{split}$
(94)

  3. (iii)

$W^{s/u}_{loc}(S)$ and $W^{s/u}_{loc}(S_{ode})$ perturb to $C^{r}-$smooth
stable/unstable manifolds $W^{s/u}_{loc}(S_{\varepsilon,h})$ and
$W^{s/u}_{loc}(S_{\varepsilon,ode})$ respectively, such that

$\textup{dist}(W^{s/u}_{loc}(S_{\varepsilon,h}),W^{s/u}_{loc}(S_{ode,\varepsilon}))=O(\varepsilon^{2}h).$
(95)

For each fixed $h$, the manifolds $W^{s/u}_{loc}(S_{\varepsilon,h})$ have all
the properties described in Theorem 3.3, while the manifolds
$W^{s/u}_{loc}(S_{ode,\varepsilon})$ are described by Fenichel theory [17, 34,
42, 73].

  4. (iv)

The manifolds $W^{s/u}_{loc}(S_{\varepsilon,h})$ and
$W^{s/u}_{loc}(S_{\varepsilon,ode})$ admit foliations by stable/unstable
fibers $w^{s/u}_{loc}(S_{\varepsilon,h})$ and
$w^{s/u}_{loc}(S_{\varepsilon,ode})$ respectively. For each fixed $h$, the
foliations of $W^{s/u}_{loc}(S_{\varepsilon,h})$ are described by Theorem 3.4,
while the foliations of $W^{s/u}_{loc}(S_{ode,\varepsilon})$ are described by
Fenichel theory. Contraction along stable fibers
$w^{s}_{loc}(S_{\varepsilon,h})$ occurs with rate $\chi_{A}(h)$ fixed and
bounded in the interval

$1-\nu_{A,ode}h+c_{A}h^{2}<\chi_{A}(h)<1,$ (96)

for a constant $c_{A}$, and repulsion along unstable fibers
$w^{u}_{loc}(S_{\varepsilon,h})$ occurs with rate $\chi_{R}(h)$ fixed and
satisfying

$\chi_{R}(h)>1+\nu_{R,ode}h+c_{R}h^{2}>1,$ (97)

for a constant $c_{R}$.

###### Proof.

The equality $S=S_{ode}$ is immediate from the form of the layer map (89).
Coincidence of the transverse linear fiber bundles $\mathcal{N}=\cup_{z\in
S}\mathcal{N}_{z}$ and $\mathcal{N}_{ode}=\cup_{z\in
S_{ode}}\mathcal{N}_{ode,z}$ follows from the pointwise equality

$\mathcal{N}_{z,ode}=\text{span}\\{N^{i}(z):1,\ldots,n-k\\}=\text{span}\\{hN^{i}(z):1,\ldots,n-k\\}=\mathcal{N}_{z},$

where $N^{i}(z)$ are the columns of $N(z)$, which holds for all $h>0$ and
$z\in S_{ode}=S$.

Except for the statement pertaining to the distance between slow manifolds in
(93), statement (i) follows directly from Fenichel theory [17, 34, 42, 73] and
Theorem 3.1. The distance in (93) follows from statement (ii) after extending
the local estimate

$|\varphi_{\varepsilon,h}(x)-\varphi_{ode,\varepsilon}(x)|=O(\varepsilon^{2}h)$

obtained from the equations in (94) via compactness and a partition of unity.
Thus it suffices to prove the local graph representations in (94). The
expansion for $\varphi_{ode,\varepsilon}(x)$ follows directly from [72, eqns
(3.33), (3.36)], and the expansion for $\varphi_{\varepsilon,h}(x)$ up to and
including $O(\varepsilon)$ follows directly from Theorem 3.2. Carrying out the
matching argument in the proof of Proposition 2.14 for the map (87) and
expanding in powers of $\varepsilon h$ instead of $\varepsilon$ leads to the
higher order correction $O(\varepsilon^{2}h)$ (a factor of $h$ cancels out
since $h$ is a common factor of $\bar{z}-z$).

Except for the statement pertaining to the distance between slow manifolds in
(95), statement (iii) follows directly from Fenichel theory [17, 34, 42, 73]
and Theorem 3.3. The estimate (95) follows from the following facts:

  * •

$W_{loc}^{s/u}(S)$ and $W_{loc}^{s/u}(S_{ode})$ are linearly tangent along
$S=S_{ode}$ since $\mathcal{N}=\mathcal{N}_{ode}$;

  * •

By assertion (i),
$\text{dist}(S_{\varepsilon,h},S_{ode,\varepsilon})=O(\varepsilon^{2}h)$;

  * •

$S_{ode,\varepsilon}\subseteq W_{loc}^{s/u}(S_{ode,\varepsilon})$ and
$S_{\varepsilon,h}\subseteq W_{loc}^{s/u}(S_{\varepsilon,h})$, since both slow
manifolds are contained within the intersection of their corresponding
perturbed stable/unstable manifolds by Theorem 3.3.

It remains to prove assertion (iv). Statements pertaining to foliations of
$W^{s/u}_{loc}(S_{ode,\varepsilon})$ follows directly from Fenichel theory
[17, 34, 42, 73]. Those pertaining to foliations of
$W^{s/u}_{loc}(S_{\varepsilon,h})$ for fixed $h$ follow from Theorem 3.4. To
see that the bounds on the contraction rate $\chi_{A}(h)$ are given by (96),
notice that by (90) we have

$|\mu_{a,j}(z)|=|1+h\lambda_{a,j,ode}(z)|=1+h\textrm{Re}\,\lambda_{a,j,ode}(z)+O(h^{2}),$

as $h\to 0^{+}$, for all $z\in S=S_{ode}$. Combining this with (92) and
$\nu_{A}=\sup_{z\in S,j=1,\ldots,n_{a}}|\mu_{a,j}(z)|$ yields the bounds on
the contraction rate in (96) via Theorem 3.4 (iii). Similar arguments yield
the bounds for the expansion rate $\chi_{R}(h)$ in (97) via Theorem 3.4 (iv).
∎

######  Remark 5.15.

The latter part of assertion (iii) as well as assertion (iv) in Theorem 5.14
hold only for fixed $h>0$ and asymptotically small $\varepsilon\sim 0$. They
do not account for asymptotic dependence in the case that both $\varepsilon\to
0$ and $h\to 0$, in which case one has an instance of the general class of
double singular limit problems described in [43].

###  5.3 Fast-slow Poincaré maps

Poincaré maps associated to parameter-dependent systems of ODEs with
hyperbolic limit cycles naturally give rise to fast-slow Poincaré maps if the
parameter is allowed to evolve slowly in time. In this section we apply the
formalism and results of Sections 2-3 in order to study the dynamics of such
systems. We begin by characterising the dynamics of the Poincaré map, before
considering the implications in the corresponding (higher dimensional) ODE.

Consider the ODE system

$\begin{split}z^{\prime}&=\tilde{f}(z,\alpha,\varepsilon),\\\
\alpha^{\prime}&=\varepsilon\tilde{g}(z,\alpha,\varepsilon),\end{split}$ (98)

where $(z,\alpha)\in\mathbb{R}^{n}\times\mathbb{R}$, $0<\varepsilon\ll 1$ and
the right-hand-side is $C^{r\geq 1}-$smooth in $(z,\varepsilon)$. System (98)
is in the standard form for fast-slow systems, with $n$ fast variables $z$ and
a single slow variable $\alpha$. For the systems considered here, however, we
do not require the existence of a critical manifold
$\\{(z,\alpha)\in\mathbb{R}^{n+1}:f(z,\alpha)=0\\}$. Thus (98) is not
necessarily slow-fast in the geometric sense of [17], see also [72, Definition
3.2]. Nevertheless, we shall refer to the limiting system

$\begin{split}z^{\prime}&=\tilde{f}(z,\alpha,0),\\\
\alpha^{\prime}&=0,\end{split}$ (99)

as $\varepsilon\to 0$ as the layer problem for system (98). Our key assumption
is the following.

######  Assumption 8.

(Existence of a hyperbolic limit cycle) The $\alpha-$family of ODEs defined by
the layer problem (99) has a hyperbolic limit cycle $\Gamma_{\alpha_{\ast}}$
of period $T_{\alpha_{\ast}}$ passing through the point $z=z_{\ast}$ for
$\alpha=\alpha_{\ast}$.

Since $\Gamma_{\alpha_{\ast}}$ is hyperbolic, it persists for nearby $\alpha$
in a sufficiently small open interval $V_{\alpha}\ni\alpha_{\ast}$ in
$\mathbb{R}$, i.e. the layer problem (99) has a hyperbolic limit cycle
$\Gamma_{\alpha}$ of period $T_{\alpha}$ for each fixed $\alpha\in
V_{\alpha}$. It follows that

$\mathcal{M}:=\left\\{\Gamma_{\alpha}:\alpha\in
V_{\alpha}\right\\}\subseteq\mathbb{R}^{n}\times V_{\alpha}$

defines a $2-$dimensional manifold of limit cycles topologically equivalent to
a cylinder in $\mathbb{R}^{n}\times V_{\alpha}$.

In the following we describe the Poincaré map
$\mathcal{P}_{\Delta}:\Delta\to\Delta$ induced by the flow of the perturbed
system (98) on a fixed codimension$-1$ cross-section $\Delta$ which is assumed
to satisfy the usual non-degeneracy conditions with respect to $\mathcal{M}$,
i.e.

  * •

$S_{\Delta}:=\Delta\cap\mathcal{M}$ defines a $1-$dimensional submanifold of
$\mathbb{R}^{n}\times V_{\alpha}$;

  * •

$T_{p}\mathcal{M}\oplus T_{p}\Delta=T_{p}\mathbb{R}^{n+1}$ at each $p\in
S_{\Delta}$;

see Figure 11.

Figure 11: Setup for $\varepsilon=0$. The layer problem (99) has a manifold of
hyperbolic limit cycles $\mathcal{M}$, shown here in shaded purple, which
contains the limit cycle $\Gamma_{\alpha_{\ast}}$ and extends over the open
$\alpha-$interval $V_{\alpha}\ni\alpha_{\ast}$. A Poincaré map
$P_{\Delta}:\Delta\to\Delta$ on the cross-section $\Delta$ shown in shaded
green is induced by the flow of the layer problem (99). Hyperbolicity of the
limit cycles $\Gamma_{\alpha}$ implies that the 1-dimensional submanifold
$S_{\Delta}=\mathcal{M}\cap\Delta$ (bold green) is a normally hyperbolic
critical manifold for $P_{\Delta}$, see Lemma 5.17.

Without loss of generality, coordinates
$z=(x,y)^{\textnormal{T}}\in\mathbb{R}^{n-1}\times\mathbb{R}$ and
$\tilde{f}=(\tilde{f}_{1},\tilde{f}_{2})^{\textnormal{T}}$ can be chosen such
that system (98) can be rewritten as

$\begin{split}x^{\prime}&=\tilde{f}_{1}(x,y,\alpha,\varepsilon),\\\
y^{\prime}&=\tilde{f}_{2}(x,y,\alpha,\varepsilon),\\\
\alpha^{\prime}&=\varepsilon\tilde{g}(x,y,\alpha,\varepsilon),\end{split}$
(100)

and the section $\Delta$ can be written as a graph

$\Delta=\left\\{(x,Y(x,\alpha),\alpha):x\in V_{x},\alpha\in
V_{\alpha}\right\\},$ (101)

where $V_{x}$ is a sufficiently small neighbourhood of the point $x_{\ast}$
defined by $z_{\ast}=(x_{\ast},y_{\ast})$, and
$Y:\mathbb{R}^{n-1}\times\mathbb{R}\to\mathbb{R}$ is $C^{r}-$smooth function.
It follows that $S_{\Delta}$ has graph representation

$S_{\Delta}=\left\\{(\varphi_{0}(\alpha),Y(\varphi_{0}(\alpha)),\alpha):\alpha\in
V_{\alpha}\right\\},$ (102)

for a $C^{r}-$smooth function $\varphi_{0}:\mathbb{R}\to\mathbb{R}^{n-1}$.

######  Remark 5.16.

In the following we shall often consider $S_{\Delta}$ only in
$(x,\alpha)-$coordinates on $\Delta$, writing simply
$S_{\Delta}=\\{(\varphi_{0}(\alpha),\alpha):\alpha\in V_{\alpha}\\}$, while
still identifying it with the intersection $\Delta\cap\mathcal{M}$. The
existence of the parameterization in (102) prevents this slight abuse of
notation from causing confusion.

We have the following result.

######  Lemma 5.17.

The Poincaré map $\mathcal{P}_{\Delta}$ is a fast-slow map satisfying
Assumptions 1-2 with normally hyperbolic critical manifold $S_{\Delta}$.

###### Proof.

The structure of the singular set and its spectrum is invariant with respect
to the choice of cross-section $\Delta$, as well as the choice of local
coordinates on $\Delta$; see e.g. [45, Lemma 1.3]. Thus it suffices to show
the result for a particular choice of $\Delta$. The derivation of the Poincaré
map is similar to derivation in the context 3-dimensional systems in [64].

We work with system (100) and assume without loss of generality that

$z_{\ast}=(x_{\ast},0)^{\textnormal{T}},\qquad\tilde{f}_{2}(x_{\ast},0,\alpha_{\ast},0)\neq
0,$

i.e. that the limit cycle $\Gamma_{\alpha_{\ast}}$ in the layer problem (99)
intersects the $n-$dimensional hyperplane

$\Sigma:=\left\\{(x,0,\alpha):x\in V_{x},\alpha\in V_{\alpha}\right\\},$

transversally at $(x_{\ast},0,\alpha_{\ast})$, for a suitable (sufficiently
small) neighbourhood $x_{\ast}\ni V_{x}\subset\mathbb{R}^{n-1}$. The
neighbourhoods $V_{x}$ and $V_{\alpha}$ can be chosen such that $\Sigma$ is a
valid choice of cross-section, i.e. such that the choice $\Delta=\Sigma$
satisfies the non-degeneracy conditions for a cross-section given in bullet
points prior to the statement of Lemma 5.17. Note that with this choice,
$Y(x,\alpha)\equiv 0$ on $V_{x}\cap V_{\alpha}$.

System (100) induces a flow

$\Phi=(\Phi_{1},\Phi_{2},\Phi_{3}):\mathbb{R}^{n-1}\times\mathbb{R}\times\mathbb{R}\times[0,\varepsilon_{0})\times\mathbb{R}\to\mathbb{R}^{n-1}\times\mathbb{R}\times\mathbb{R}$

which by Assumption 8 satisfies

$\Phi_{2}(x_{\ast},0,\alpha_{\ast},0,T_{\alpha_{\ast}})=0\qquad\text{and}\qquad\Phi_{2}^{\prime}(x_{\ast},0,\alpha_{\ast},0,T_{\alpha_{\ast}})\neq
0.$

Therefore, the implicit function theorem implies the existence of a smooth and
(locally) invertible function $\tilde{t}:V_{x}\times V_{y}\times
V_{\alpha}\times[0,\varepsilon_{0})\to\mathbb{R}$ such that
$\tilde{t}(x_{\ast},0,\alpha_{\ast},0)=T_{\alpha_{\ast}}$ and

$\Phi_{2}(x,y,\alpha,\varepsilon,\tilde{t}(x,y,\alpha,\varepsilon))\equiv 0$

on $V_{x}\times V_{y}\times V_{\alpha}\times[0,\varepsilon_{0})$, where
$V_{y}$ is a neighbourhood about $y_{\ast}=0$ in $\mathbb{R}$. In terms of
earlier notation, a limit cycle $\Gamma_{\alpha}\subset\mathcal{M}$ such that
$\Gamma_{\alpha}\cap S_{\Sigma}=\\{(\varphi_{0}(\alpha),0,\alpha)\\}$ has
period $T_{\alpha}=\tilde{t}(\varphi_{0}(\alpha),0,\alpha,0)$.

The preceding argument guarantees the existence of a well-defined Poincaré map
$\mathcal{P}_{\Sigma}:\Sigma\to\Sigma$ such that
$\mathcal{P}_{\Sigma}(S_{\Sigma})=S_{\Sigma}$. Explicitly,

$\mathcal{P}_{\Sigma}:\begin{pmatrix}x\\\
\alpha\end{pmatrix}\mapsto\begin{pmatrix}\bar{x}\\\
\bar{\alpha}\end{pmatrix}=\begin{pmatrix}x\\\
\alpha\end{pmatrix}+\begin{pmatrix}\widehat{f}(x,\alpha,\varepsilon)\\\
\varepsilon g(x,\alpha,\varepsilon)\end{pmatrix},$ (103)

with

$\begin{split}\widehat{f}(x,\alpha,\varepsilon)&:=\Phi_{1}(x,0,\alpha,\varepsilon,\tilde{t}(x,0,\alpha,\varepsilon))-x,\\\
g(x,\alpha,\varepsilon)&:=D_{\varepsilon}\Phi_{3}(x,0,\alpha,0,\tilde{t}(x,0,\alpha,0))+O(\varepsilon),\end{split}$
(104)

where the latter expression follows after expansion in $\varepsilon$. In
particular, the function $g(x,\alpha,\varepsilon)$ can be expressed in terms
of the original function
$\tilde{g}(z,\alpha,\varepsilon)=\tilde{g}(x,y,\alpha,\varepsilon)$ via the
integral formula

$g(x,\alpha,\varepsilon)=\int_{0}^{\tilde{t}(x,0,\alpha,0)}\tilde{g}\left(\Phi_{1}(x,0,\alpha,0,s),0,\alpha,0\right)ds+O(\varepsilon).$

We need to verify Assumptions 1-2 for the map (103), as well as the normal
hyperbolicity of $S_{\Sigma}$. Notice first that (103) can be rewritten in the
general form (10) after defining
$\widehat{f}(x,\alpha,\varepsilon)=:\widehat{f}(x,\alpha,0)+f_{rem}(x,\alpha,\varepsilon)$
and

$N(x,\alpha):=\begin{pmatrix}I_{n-1}\\\ O_{1,n-1}\end{pmatrix},\qquad
f(x,\alpha):=\widehat{f}(x,\alpha,0),\qquad\tilde{G}(x,\alpha,\varepsilon):=\begin{pmatrix}f_{rem}(x,\alpha,\varepsilon)\\\
g(x,\alpha,\varepsilon)\end{pmatrix}.$

We have that

$f(x_{\ast},\alpha_{\ast})=\Phi_{1}(x_{\ast},0,\alpha_{\ast},0,T)-x_{\ast}=\Phi_{1}(x_{\ast},0,\alpha_{\ast},0,0)-x_{\ast}=x_{\ast}-x_{\ast}=0,$

and, by Assumption 8, that the $n-1$ characteristic multipliers $\mu_{j}$ of
the matrix $I_{n-1}+D_{x}f(x_{\ast},\alpha_{\ast})$ satisfy $|\mu_{j}|\neq 1$.
It follows that the $n-1$ eigenvalues $\lambda_{j}$ of the matrix
$D_{x}f(x_{\ast},\alpha_{\ast})$ satisfy $\text{Re}\ \lambda_{j}\neq 0$, i.e.
it is regular, so that the implicit function theorem implies the existence of
a $C^{r}-$smooth function $\tilde{\varphi}_{0}:V_{\alpha}\to V_{x}$ such that
$\tilde{\varphi}_{0}(\alpha_{\ast})=x_{\ast}$ and
$f(\varphi_{0}(\alpha),\alpha)\equiv O_{n-1}$ on $V_{\alpha}$. The fact that
$\tilde{\varphi}_{0}=\varphi_{0}$, i.e. that
$S_{\Sigma}=\\{(\varphi_{0}(\alpha),\alpha):\alpha\in
V_{\alpha}\\}=\\{(x,\alpha):f(x,\alpha)=O_{n-1}\\}$ defines the critical
manifold of (103), follows from the local uniqueness of $\varphi$ and
$\tilde{\varphi}_{0}$ and their coincidence at $(x_{\ast},\alpha_{\ast})$.
Thus, (103) is a fast-slow map satisfying Assumption 1. Assumption 2 is
immediate since we consider a neighbourhood with $(x,\alpha)\in V_{x}\times
V_{\alpha}$ within which the only zeroes of $N(x,\alpha)f(x,\alpha)$ are those
identified above for $f(x,\alpha)=O_{n-1}$ (uniqueness follows from the
implicit function theorem).

It remains to verify that $S_{\Sigma}$ is normally hyperbolic. By Proposition
2.8, non-trivial multipliers along $S_{\Sigma}$ are in $1-1$ correspondence
with eigenvalues of the $(n-1)\times(n-1)$ square matrix

$I_{n-1}+N(\varphi_{0}(\alpha),\alpha)Df(\varphi_{0}(\alpha),\alpha)=I_{n-1}+D_{x}f(\varphi_{0}(\alpha),\alpha).$
(105)

For each fixed $\alpha\in V_{\alpha}$ the $n-1$ multipliers $\mu_{j}(\alpha)$
of the matrix (105) are $1-1$ correspondence with the characteristic
multipliers associated a hyperbolic limit cycle
$\Gamma_{\alpha}\subset\mathcal{M}$, and therefore satisfy
$|\mu_{j}(\alpha)|\neq 1$ for all $\alpha\in V_{\alpha}$. Thus $S_{\Sigma}$ is
normally hyperbolic. ∎

Lemma 5.17 implies that the results of Section 3 can be applied directly to
the Poincaré map $\mathcal{P}_{\Delta}$. The derivation of the Poincaré map on
$\Delta=\Sigma$ generalises straightforwardly to the case where $\Delta$ is
given by the graph (101); one must simply restrict to $y=Y(x,\alpha)$ instead
in the defining expressions (104). Explicitly,

$\mathcal{P}_{\Delta}:\begin{pmatrix}x\\\
\alpha\end{pmatrix}\mapsto\begin{pmatrix}\bar{x}\\\
\bar{\alpha}\end{pmatrix}=\begin{pmatrix}x\\\
\alpha\end{pmatrix}+\begin{pmatrix}\widehat{f}(x,\alpha,\varepsilon)\\\
\varepsilon g(x,\alpha,\varepsilon)\end{pmatrix},$ (106)

with

$\begin{split}\widehat{f}(x,\alpha,\varepsilon)&:=\Phi_{1}(x,Y(x,\alpha),\alpha,\varepsilon,\tilde{t}(x,Y(x,\alpha),\alpha,\varepsilon))-x,\\\
g(x,\alpha,\varepsilon)&:=D_{\varepsilon}\Phi_{3}(x,Y(x,\alpha),\alpha,0,\tilde{t}(x,Y(x,\alpha),\alpha,0))+O(\varepsilon),\end{split}$
(107)

where in particular we have

$g(x,\alpha,\varepsilon)=\int_{0}^{\tilde{t}(x,Y(x,\alpha),\alpha,0)}\tilde{g}\left(\Phi_{1}(x,Y(x,\alpha),\alpha,0,s),Y(x,\alpha),\alpha,0\right)ds+O(\varepsilon).$
(108)

Our main results on the persistence of the critical manifold $S_{\Delta}$, its
local stable/unstable manifolds $W^{s/u}_{loc}(S_{\Delta})$ and their
foliations for $0<\varepsilon\ll 1$, are stated for the map
$\mathcal{P}_{\Delta}$ defined in equations (106)-(107).

######  Theorem 5.18.

Consider the Poincaré map $\mathcal{P}_{\Delta}$ defined by equations
(106)-(107). There exists and $\varepsilon_{0}>0$ such that following
properties hold for all $\varepsilon\in(0,\varepsilon_{0})$:

  1. (i)

The critical manifold $S_{\Delta}$ persists as a locally invariant slow
manifold $S_{\Delta,\varepsilon}$ described by Theorems 3.1-3.2. In
particular,

$S_{\Delta,\varepsilon}=\\{(\varphi_{\varepsilon}(\alpha),\alpha):\alpha\in
V_{\alpha}\\},$

where

$\varphi_{\varepsilon}(\alpha)=\varphi_{0}(\alpha)-\varepsilon\left((D_{x}f)f_{rem}+(D_{\alpha}g)\right)\big{|}_{S_{\Delta}}+O(\varepsilon^{2}),$

and slow iteration along $S_{\varepsilon}$ is governed by

$\mathcal{P}_{\Delta}|_{S_{\Delta,\varepsilon}}:\begin{pmatrix}x\\\
\alpha\end{pmatrix}\mapsto\begin{pmatrix}\bar{x}\\\
\bar{\alpha}\end{pmatrix}=\begin{pmatrix}\varphi_{0}(\alpha)\\\
\alpha\end{pmatrix}+\varepsilon\begin{pmatrix}-(D_{x}f)^{-1}(D_{\alpha}f)g\\\
g\end{pmatrix}\bigg{|}_{S_{\Delta}}+O(\varepsilon^{2}),$ (109)

as long as iterates remain in $V_{x}\times V_{\alpha}$.

  2. (ii)

The stable/unstable manifolds $W_{loc}^{s/u}(S_{\Delta})$ persist as locally
invariant manifolds $W^{s/u}_{loc}(S_{\Delta,\varepsilon})$ described by
Theorem 3.3.

  3. (iii)

Foliations of the stable/unstable manifolds $W_{loc}^{s/u}(S_{\Delta})$
persist as locally invariant foliations of
$W_{loc}^{s/u}(S_{\Delta,\varepsilon})$ described by Theorem 3.4. The spectral
bounds $\nu_{A}$ and $\nu_{R}$ determining the contraction and repulsion rates
in Theorem 3.4 (iii) and (iv) respectively are given by

$\nu_{A}=\sup_{\alpha\in
V_{\alpha},j=1,\ldots,n-1}\\{|\mu_{j}(\alpha)|:|\mu_{j}(\alpha)|<1\\},\ \
\nu_{R}=\sup_{\alpha\in
V_{\alpha},j=1,\ldots,n-1}\\{|\mu_{j}(\alpha)|:|\mu_{j}(\alpha)|>1\\},$

where $\mu_{j}(\alpha)$, $j=1,\ldots,n-1$ denote the characteristic
multipliers associated to the limit cycle $\Gamma_{\alpha}\subset\mathcal{M}$
in the layer problem (99), which are in $1-1$ correspondence with multipliers
of the matrix $I_{n-1}+D_{x}f{(\varphi_{0}(\alpha),\alpha)}$.

###### Proof.

It follows by Lemma 5.17 that the results of Sections 2-3 can be applied
directly to the map $P_{\Delta}$.

The existence and form of the parameterization of the slow manifold
$S_{\Delta,\varepsilon}$ is assertions (i) follows directly from Theorems
3.1-3.2, and the form of the map describing slow iteration along
$S_{\Delta,\varepsilon}$ follows from Proposition 2.14.

Assertions (ii) and (iii) follow directly from Theorem 3.3 and Theorem 3.4
respectively. ∎

Figure 12: Perturbation of the singular geometry and dynamics for
$\varepsilon=0$ (left) for $0\leq\varepsilon\ll 1$ (right). We sketch the
3-dimensional case of a normally hyperbolic and attracting limit cycle
manifold $\mathcal{M}$. Theorem 5.18 describes the persistence properties for
the Poincaré map $P_{\Delta}$ on an arbitrary cross-section $\Delta$. Three
representative stable fibers on $\Delta$ for $\varepsilon=0$ are shown in blue
on the left, which perturb to invariant stable fibers for $0<\varepsilon\ll 1$
(also in blue) on the right, e.g. $w_{loc}^{s}(z)$ perturbs to
$w_{loc}^{s}(z_{\varepsilon})$. A number of example iterates are shown on
$\Delta$ for both $\varepsilon=0$ and $0<\varepsilon\ll 1$. Note that the
arrows connecting points on e.g. $w_{loc}^{s}(z_{\varepsilon})$ and
$w_{loc}^{s}({\bar{z}_{\varepsilon}}^{5})$ only show the direction of
iteration (they do not indicate a single iteration). The persistence of the
manifold of limit cycles $\mathcal{M}$ as a nearby locally invariant ‘slow
manifold’ $\mathcal{M}_{\varepsilon}$ as described by Theorem 5.18 is also
shown ($\mathcal{M}_{\varepsilon}$ is shown in a darker shade or purple).

The situation is sketched in Figure 12.

Our final result for this section characterises the persistence of the limit
cycle manifold $\mathcal{M}$ in the ODE system (98) for $0<\varepsilon\ll 1$.
We also provide a criterion for the existence of hyperbolic limit cycles in
the perturbed system (98).

######  Theorem 5.19.

Consider the fast-slow ODE system (98) under Assumption 8. Then there exists
and $\varepsilon_{0}>0$ such that for all $\varepsilon\in(0,\varepsilon_{0})$,
the limit cycle manifold $\mathcal{M}$ perturbs to a locally invariant
manifold $\mathcal{M}_{\varepsilon}$ which is $O(\varepsilon)-$close,
$C^{r}-$smooth and diffeomorphic to $\mathcal{M}$.

The manifold $\mathcal{M}_{\varepsilon}$ contains a hyperbolic limit cycle of
the ODE (98) intersecting the point with local coordinates
$(x_{lc},y_{lc},\alpha_{lc})\in\mathcal{M}_{\varepsilon}$ if and only if

$g(\varphi_{0}(\alpha),\alpha,0)=0,\qquad\text{and}\qquad
D_{\alpha}g(\varphi_{0}(\alpha),\alpha,0)\neq 0,$ (110)

where $g$ is defined in the $(x,y,\alpha)-$coordinates of the Poincaré map
(107) via the integral formula (108).

###### Proof.

The limit cycle manifold $\mathcal{M}$ can be parameterised by a one-parameter
family of submanifolds $S_{\sigma}$ defined via the flow of the layer problem
(99). Specifically, define

$S_{\sigma}:=\Phi_{\sigma
T_{\alpha}}(S)=\\{\Phi(\varphi_{0}(\alpha),Y(\varphi_{0}(\alpha),\alpha),\alpha,0,\sigma
T_{\alpha}):\alpha\in V_{\alpha}\\},\quad\sigma\in[0,1),$

where we set $S_{0}:=S=\mathcal{M}\cap\Delta$, such that
$\mathcal{M}=\cup_{\sigma\in[0,1)}S_{\sigma}$. Since it is always possible to
associate a non-degenerate cross-section $\Delta_{\sigma}$ such that
$S_{\sigma}=\Delta_{\sigma}\cap\mathcal{M}$, a family of Poincaré maps
$\\{\mathcal{P}_{\Delta_{\sigma}}:\sigma\in[0,1)\\}$ can be obtained by
constructing each
$\mathcal{P}_{\Delta_{\sigma}}:\Delta_{\sigma}\to\Delta_{\sigma}$ analogously
to the map (103) described in Theorem 5.18. In particular, $S_{\sigma}$ is a
normally hyperbolic critical manifold for $\mathcal{P}_{\sigma}$, for each
$\sigma\in[0,1)$.

By Theorem 5.18 (i), each $S_{\sigma}$ perturbs to a nearby slow manifold
$S_{\sigma,\varepsilon}\subset\Delta_{\sigma}$. Thus for any fixed choice of
$\sigma$, we may define the manifold $\mathcal{M}_{\varepsilon}$ by

$\mathcal{M}_{\varepsilon}:=\Phi\left(S_{\sigma,\varepsilon},\varepsilon,[-\tilde{T},\tilde{T}]\right)\cap\\{(x,y,\alpha):\alpha\in
V_{\alpha}\\}$

where $\tilde{T}>0$ is fixed larger than the maximum of
$T_{\alpha,\varepsilon}$ on $V_{\alpha}\times[0,\varepsilon_{0})$. Since the
flow map $\Phi$ is $C^{r}-$smooth, $\mathcal{M}_{\varepsilon}$ is exactly as
smooth as $S_{\sigma,\varepsilon}$, which by Theorem 5.18 is exactly as smooth
as $S_{\sigma}$, which is exactly as smooth as $\mathcal{M}$. The fact that
$\mathcal{M}_{\varepsilon}$ is $O(\varepsilon)-$close and diffeomorphic to
$\mathcal{M}$ follows from the corresponding results for
$S_{\sigma,\varepsilon}$ and $S_{\sigma}$ implied by Theorem 5.17 (i) and
regular perturbation theory.

It remains to prove the statement pertaining to limit cycles in
$\mathcal{M}_{\varepsilon}$. We work in local coordinates $(x,y,\alpha)\in
V_{x}\times V_{y}\times V_{\alpha}$ for which we have
$x_{lc}=\varphi_{\varepsilon}(\alpha_{lc})$ and
$y_{lc}=Y(\varphi_{0}(\alpha_{lc}),\alpha_{lc})$ since we may assume that
$(x_{lc},y_{lc},\alpha_{lc})\in
S_{\Delta,\varepsilon}=\Delta\cap\mathcal{M}_{\varepsilon}$. Limit cycles of
(98) in $\mathcal{M}_{\varepsilon}$ are in 1$-$1 correspondence with fixed
points of the restricted Poincaré map on
$P_{\Delta}|_{S_{\Delta,\varepsilon}}$ given by (109). Since the matrix
$D_{x}f$ is locally regular under normally hyperbolic conditions, fixed points
and their hyperbolicity are completely determined by the restricted,
1-dimensional map

$\alpha\mapsto\bar{\alpha}=\alpha+\varepsilon
g(\varphi_{\varepsilon}(\alpha),\alpha,\varepsilon)=\alpha+\varepsilon
g(\varphi_{0}(\alpha),\alpha,0)+O(\varepsilon^{2}).$ (111)

It follows from the implicit function theorem that the map (111) has a
hyperbolic fixed point for all $0<\varepsilon\ll 1$ sufficiently small if and
only if the conditions in (110) are satisfied. This concludes the proof. ∎

Figure 12 shows both $\mathcal{M}$ and it’s perturbation
$\mathcal{M}_{\varepsilon}$, which intersects $\Delta$ along the slow manifold
$S_{\Delta,\varepsilon}$ of the Poincaré map $P_{\Delta}$.

######  Remark 5.20.

Equation (108) implies an integral formulation of the conditions in (110) in
terms of the expression

$g(\varphi_{0}(\alpha),\alpha,0)=\int_{0}^{T_{\alpha}}\tilde{g}\left(\Phi_{1}(\varphi_{0}(\alpha),Y(\varphi_{0}(\alpha),\alpha),\alpha,0,s),Y(\varphi_{0}(\alpha),\alpha),\alpha,0\right)ds,$
(112)

where
$T_{\alpha}=\tilde{t}(\varphi_{0}(\alpha),Y(\varphi_{0}(\alpha),\alpha),\alpha,0)$
is the period associated to the limit cycle
$\Gamma_{\alpha}\subset\mathcal{M}$ of the layer problem (99) which passes
through the point
$(\varphi_{0}(\alpha),Y(\varphi_{0}(\alpha),\alpha),\alpha)\in S_{\Delta}$. In
classical averaging theory, the integral (112) defines the so-called averaged
equation [22, 60].

######  Remark 5.21.

Although the contraction and expansion along invariant foliations under the
Poincaré map $P_{\Delta,\varepsilon}$ are described in detail by Theorem 5.18,
the persistence of stable/unstable manifolds $W^{u/s}_{loc}(\mathcal{M})$ and
their corresponding foliations in the ODE (98) is not described in Theorem
5.19. A detailed proof of persistence is left for future work.

##  6 Outlook

Discrete multi-scale dynamical systems induced by maps arise in a wide variety
applied and theoretical settings. However, to the best of our knowledge, a
complete geometric theory for their analysis analogous to the established GSPT
for multi-scale continuous-time systems does not yet exist. The aim of the
present manuscript is to provide such a theory, referred to herein as DGSPT,
for the class of fast-slow maps defined by (3) or, more precisely, general
maps (10) under Assumptions 1 and 2.

The singular theory is presented in Section 2, where the layer map is
introduced and used to define a discrete notion of normal hyperbolicity in
terms of a spectral condition on the non-trivial multipliers along the
critical manifold $S$, recall Definition 2.6. By Proposition 2.8, the problem
of checking for normal hyperbolicity of $S$ reduces to a direct evaluation of
the spectrum associated to the matrix $I_{n-k}+DfN|_{S}$, which is given
purely in terms of $N$ and $f$. Under normally hyperbolic assumptions, we
introduced the reduced map (24) in Definition 2.16. Due to the fact that there
is no discrete analogue of the time rescaling $\tau=\varepsilon t$ in the
continuous-time setting, the reduced map depends necessarily on $\varepsilon$
in the leading order. By Proposition 2.14, however, it nevertheless
approximates slow iteration along perturbed slow manifolds $S_{\varepsilon}$
up to an accuracy of $O(\varepsilon^{2})$, which is typically all that is
required in applications. An expression for the $m$’th iterate map on
$S_{\varepsilon}$ was also derived using the asymptotic self-similarity of the
map near $S$, see again Proposition 2.19. Interestingly, this map can to
leading order in $\varepsilon$ be related to the Euler discretization of the
corresponding continuous-time reduced problem. It is also possible that the
reduced $m$’th iterate map is well-defined on $S$ under a suitable formulation
of the dual limit $(\varepsilon,m)\to(0,\infty)$, recall Remark 2.20, however
the details remain for future work.

The persistence of key dynamic and geometric features for $0<\varepsilon\ll 1$
under normally hyperbolic conditions is described in detail by the Fenichel-
like theorems provided in Section 3. The persistence of a compact, normally
hyperbolic critical manifold as a nearby locally invariant slow manifold
$S_{\varepsilon}$ is described by Theorems 3.1 and 3.2, with the latter
providing an explicit parameterization for $S_{\varepsilon}$ up to
$O(\varepsilon^{2})$ in local coordinates. The persistence of local
stable/unstable manifolds $W_{loc}^{s/u}(S)$ and their corresponding
asymptotic rate foliations are described by Theorems 3.3 and 3.4 respectively.
Theorem 3.4 in particular gives quantitative estimates for the contraction and
expansion rates close to $S_{\varepsilon}$ in terms of the scalar quantities
$\nu_{A}$ and $\nu_{R}$ defined in (41), which define the (annular) spectral
gap about the unit circle associated to the matrix $I_{n-k}+DfN|_{S}$. While
the invariant manifold theorems in Section 3 are in some regards less general
than their counterparts in [54] (e.g. we require $C^{1}-$smoothness as opposed
to Lipschitz continuity), the ‘Fenichel-like’ formulation of the results in
Section 3 are expected to be advantageous in many applications due the fact
that they (i) apply directly to fast-slow maps (3) without the need for prior
transformations into special coordinates, and (ii) apply directly near compact
submanifolds of $S$ in the locally invariant setting.

The utility of DGSPT was demonstrated in Section 5 in the context of three
non-trivial applications. In Section 5.1 we used DGSPT in order to provide a
geometric analysis of the Chialvo map-based model for neuronal bursting in a
fast-slow parameter regime, thereby validating and extending previous work in
[5, 33, 68]. The DGSPT formalism provided a systematic approach to the
identification of distinct singular structure corresponding to four types of
neuronal dynamics: excitability, relaxation, regular (non-chaotic) bursting
and chaotic bursting, recall Figure 10. In Section 5.2 we used DGSPT to prove
results on the relationship between the geometry and dynamics of fast-slow
ODEs in the general (non-standard) form (86) and its corresponding Euler
discretization, see Theorem 5.14. These results extend pre-existing results on
the discretization of fast-slow ODEs in standard form (2) in e.g. [52, 53].
Finally in Section 5.3, we showed that fast-slow ODE systems with a single
slow variable (which often arise in standard form (2) after allowing a
parameter to evolve slowly in time) give rise to fast-slow Poincaré maps with
a normally hyperbolic critical manifold if the layer problem has a hyperbolic
limit cycle, see Lemma 5.17. This critical manifold, which lies within the
intersection of a limit cycle manifold $\mathcal{M}$ of the (continuous-time)
layer problem and a transverse section $\Delta$, perturbs to a nearby slow
manifold for the Poincaré map for $0<\varepsilon\ll 1$. The geometry and
dynamics of the Poincaré map are described in detail in Theorem 5.18, which
follows after a direct application of the invariant manifold theorems in
Section 3. Properties of the Poincaré map were then applied in order to prove
Theorem 5.19, which describes the persistence of the limit cycle manifold
$\mathcal{M}$ as a nearby locally invariant manifold
$\mathcal{M}_{\varepsilon}$ in the ODE and, additionally, provides a criteria
for the existence and hyperbolicity of limit cycles within
$\mathcal{M}_{\varepsilon}$.

The present work aims to provides a mathematical framework for the analysis of
a large class of multi-scale discrete dynamical systems, thereby paving the
way for a wealth of future work. In particular, the question of how to
supplement DGSPT with geometric methods for studying the non-normally
hyperbolic regime frequently arose in this work. The singular theory of
Section 2 provides the means for classifying non-normally hyperbolic
singularities on $S$ in terms of ‘$\varepsilon=0$ conditions’ on the layer map
(8) and the function $\mathcal{G}(z):=\Pi_{\mathcal{N}}^{S}G(z,0)|_{S}$ which
determines the dynamics of the reduced map (24). A detailed treatment of the
codimension-1 singularities, i.e. the fast-slow fold, flip/period-doubling and
Neimark-Sacker/torus singularities, in the low dimensional setting using DGSPT
in combination with adaptations of the geometric blow-up technique constitutes
very recent, ongoing and future work.

##  7 Acknowledgements

SJ and CK acknowledge funding from the SFB/TRR 109 Discretization and Geometry
in Dynamics grant. SJ would like to thank Prof. Martin Wechselberger for many
helpful discussions pertaining to geometric singular perturbation theory
“beyond the standard form” in general. CK thanks the VolkswagenStiftung for
support via a Lichtenberg Professorship.

## References

  * [1] L. Arcidiacono, M. Engel, and C. Kuehn, Discretized fast–slow systems near pitchfork singularities, Journal of Difference Equations and Applications, 25 (2019), pp. 1024–1051.
  * [2] C. Baesens, Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation, Physica D: Nonlinear Phenomena, 53 (1991), pp. 319–375.
  * [3] E. Benoit, Dynamic bifurcations: Proceedings of a conference held in Luminy, France, March 5-10, 1990, Springer, 2006.
  * [4] P. Carter and B. Sandstede, Unpeeling a homoclinic banana in the Fitzhugh–Nagumo system, SIAM Journal on Applied Dynamical Systems, 17 (2018), pp. 236–349.
  * [5] D. R. Chialvo, Generic excitable dynamics on a two-dimensional map, Chaos, Solitons & Fractals, 5 (1995), pp. 461–479.
  * [6] M. Courbage, V. Nekorkin, and L. Vdovin, Chaotic oscillations in a map-based model of neural activity, Chaos: An Interdisciplinary Journal of Nonlinear Science, 17 (2007), p. 043109.
  * [7] P. De Maesschalck, T. S. Doan, and J. Wynen, Intrinsic determination of the criticality of a slow–fast hopf bifurcation, Journal of Dynamics and Differential Equations, (2020), pp. 1–17.
  * [8] P. De Maesschalck, F. Dumortier, and R. Roussarie, Canard cycles, vol. 73, Springer, 2021.
  * [9] P. De Maesschalck and M. Wechselberger, Neural excitability and singular bifurcations, The Journal of Mathematical Neuroscience (JMN), 5 (2015), pp. 1–32.
  * [10] Y. Do and J. M. Lopez, Slow passage through multiple bifurcation points, Discrete & Continuous Dynamical Systems-B, 18 (2013), p. 95.
  * [11] F. Dumortier and R. Roussarie, Canard cycles and center manifolds, no. 577 in Memoirs of the American Mathematical Society, American Mathematical Society, 1996.
  * [12] M. Engel and H. Jardón-Kojakhmetov, Extended and symmetric loss of stability for canards in planar fast-slow maps, SIAM Journal on Applied Dynamical Systems, 19 (2020), pp. 2530–2566.
  * [13] M. Engel and C. Kuehn, Discretized fast-slow systems near transcritical singularities, Nonlinearity, 32 (2019), p. 2365.
  * [14] M. Engel, C. Kuehn, M. Petrera, and Y. Suris, Discretized fast–slow systems with canards in two dimensions, Journal of Nonlinear Science, 32 (2022), pp. 1–41.
  * [15] N. Fenichel, Asymptotic stability with rate conditions, Indiana University Mathematics Journal, 23 (1974), pp. 1109–1137.
  * [16] , Asymptotic stability with rate conditions, ii, Indiana University Mathematics Journal, 26 (1977), pp. 81–93.
  * [17] , Geometric singular perturbation theory for ordinary differential equations, Journal of Differential Equations, 31 (1979), pp. 53–98.
  * [18] N. Fenichel and J. Moser, Persistence and smoothness of invariant manifolds for flows, Indiana University Mathematics Journal, 21 (1971), pp. 193–226.
  * [19] A. Fruchard and R. Schäfke, A survey of some results on overstability and bifurcation delay, Discrete & Continuous Dynamical Systems-S, 2 (2009), p. 931.
  * [20] A. Goeke and S. Walcher, A constructive approach to quasi-steady state reductions, Journal of Mathematical Chemistry, 52 (2014), pp. 2596–2626.
  * [21] J. Guckenheimer, Towards a global theory of singularly perturbed dynamical systems, in Nonlinear Dynamical Systems and Chaos, Springer, 1996, pp. 213–225.
  * [22] J. Guckenheimer and P. Holmes, Nonlinear oscillations, dynamical systems, and bifurcation of vector fields, Springer-Verlag, Berlin Heidelberg New York, 1983.
  * [23] M. G. Hayes, T. J. Kaper, P. Szmolyan, and M. Wechselberger, Geometric desingularization of degenerate singularities in the presence of fast rotation: A new proof of known results for slow passage through Hopf bifurcations, Indagationes Mathematicae, 27 (2016), pp. 1184–1203.
  * [24] M. W. Hirsch, C. C. Pugh, and M. Shub, Invariant manifolds, Bulletin of the American Mathematical Society, 76 (1970), pp. 1015–1019.
  * [25] A. Huber and P. Szmolyan, Geometric singular perturbation analysis of the Yamada model, SIAM Journal on Applied Dynamical Systems, 4 (2005), pp. 607–648.
  * [26] B. Ibarz, J. M. Casado, and M. A. Sanjuán, Map-based models in neuronal dynamics, Physics reports, 501 (2011), pp. 1–74.
  * [27] E. M. Izhikevich, Simple model of spiking neurons, IEEE Transactions on neural networks, 14 (2003), pp. 1569–1572.
  * [28] E. M. Izhikevich and F. Hoppensteadt, Classification of bursting mappings, International Journal of Bifurcation and Chaos, 14 (2004), pp. 3847–3854.
  * [29] H. Jardón-Kojakhmetov and C. Kuehn, A survey on the blow-up method for fast-slow systems, arXiv preprint arXiv:1901.01402, (2019).
  * [30] H. Jardón-Kojakhmetov, J. M. Scherpen, and D. del Puerto-Flores, Stabilization of a class of slow–fast control systems at non-hyperbolic points, Automatica, 99 (2019), pp. 13–21.
  * [31] S. Jelbart, N. Pages, V. Kirk, J. Sneyd, and M. Wechselberger, Process-oriented geometric singular perturbation theory and calcium dynamics, arXiv preprint arXiv:2104.07304, to appear in SIAM Journal on Applied Dynamical Systems, (2021).
  * [32] S. Jelbart and M. Wechselberger, Two-stroke relaxation oscillators, Nonlinearity, 33 (2020), p. 2364.
  * [33] Z. Jing, J. Yang, and W. Feng, Bifurcation and chaos in neural excitable system, Chaos, Solitons & Fractals, 27 (2006), pp. 197–215.
  * [34] C. K. Jones, Geometric singular perturbation theory, in Dynamical systems, vol. 1609 of Lecture Notes in Mathematics, Springer, 1995, pp. 44–118.
  * [35] P. I. Kaleda, Singular systems on the plane and in space, Journal of Mathematical Sciences, 179 (2011), pp. 475–490.
  * [36] T. J. Kaper, Systems theory for singular perturbation problems, Analyzing multiscale phenomena using singular perturbation methods: American Mathematical Society short course, January 5–6, 1998, Baltimore, Maryland, 56 (1999), p. 85.
  * [37] I. Kosiuk and P. Szmolyan, Geometric singular perturbation analysis of an autocatalator model, Discrete and Continuous Dynamical Systems, 2 (2009), pp. 783–806.
  * [38] , Scaling in singular perturbation problems: blowing up a relaxation oscillator, SIAM Journal on Applied Dynamical Systems, 10 (2011), pp. 1307–1343.
  * [39] , Geometric analysis of the Goldbeter minimal model for the embryonic cell cycle, Journal of mathematical biology, 72 (2016), pp. 1337–1368.
  * [40] N. Kruff and S. Walcher, Coordinate-independent singular perturbation reduction for systems with three time scales, Mathematical Biosciences and Engineering, 16 (2019), pp. 5062–5091.
  * [41] M. Krupa and P. Szmolyan, Extending geometric singular perturbation theory to nonhyperbolic points—fold and canard points in two dimensions, SIAM Journal on Mathematical Analysis, 33 (2001), pp. 286–314.
  * [42] C. Kuehn, Multiple time scale dynamics, vol. 191 of Applied Mathematical Sciences, Springer, 2015.
  * [43] C. Kuehn, N. Berglund, C. Bick, M. Engel, T. Hurth, A. Iuorio, and C. Soresina, A general view on double limits in differential equations, arXiv preprint arXiv:2106.01160, (2021).
  * [44] C. Kuehn and P. Szmolyan, Multiscale geometry of the Olsen model and non-classical relaxation oscillations, Journal of Nonlinear Science, 25 (2015), pp. 583–629.
  * [45] Y. A. Kuznetsov, Elements of applied bifurcation theory, vol. 112 of Applied Mathematical Sciences, Springer Science & Business Media, 2013.
  * [46] Y. A. Kuznetsov and H. G. Meijer, Numerical bifurcation analysis of maps: From theory to software, vol. 34, Cambridge University Press, 2019.
  * [47] I. Lizarraga, R. Marangell, and M. Wechselberger, Slow unfoldings of contact singularities in singularly perturbed systems beyond the standard form, Journal of Nonlinear Science, 30 (2020), pp. 3161–3198.
  * [48] I. Lizarraga, B. Rink, and M. Wechselberger, Multiple timescales and the parametrisation method in geometric singular perturbation theory, Nonlinearity, 34 (2021), p. 4163.
  * [49] I. Lizarraga and M. Wechselberger, Computational singular perturbation method for nonstandard slow-fast systems, SIAM Journal on Applied Dynamical Systems, 19 (2020), pp. 994–1028.
  * [50] V. Nekorkin and L. Vdovin, Diskretnaya model neyronnoy aktivnosti, Izvest vys ucheb zaved Prikladnaya nelinejnaya dinamika, 15 (2007), pp. 36–60.
  * [51] K. Nipp and D. Stoffer, Attractive invariant manifolds for maps: existence, smoothness and continuous dependence on the map, no. 92,11 in Research report / Eidgenössische Technische Hochschule Zürich, Seminar für Angewandte Mathematik, ETH, Zürich, 1992.
  * [52] , Invariant manifolds and global error estimates of numerical integration schemes applied to stiff systems of singular perturbation type–part i: Rk-methods, Numerische Mathematik, 70 (1995), pp. 245–257.
  * [53] , Invariant manifolds and global error estimates of numerical integration schemes applied to stiff systems of singular perturbation type–part ii: Linear multistep methods, Numerische Mathematik, 74 (1996), pp. 305–323.
  * [54] , Invariant Manifolds in Discrete and Continuous Dynamical Systems, vol. 21 of EMS Tracts in Mathematics, European Mathematical Society, 2013.
  * [55] K. Nipp, D. Stoffer, and P. Szmolyan, Graph transform and blow-up in singular perturbations, in AIP Conference Proceedings, American Institute of Physics, 2009, pp. 861–868.
  * [56] W. F. Osgood, Lehrbuch der Funktionentheorie, no. 20,1 in B. G. Teubners Sammlung von Lehrbüchern auf dem Gebiete der Mathematischen Wissenschaften mit Einschluss ihrer Anwendungen, Teubner, Leipzig, 2. aufl. ed., 1912.  Literaturangaben.
  * [57] C. Pötzsche, Slow and fast variables in non-autonomous difference equations, Journal of Difference Equations and Applications, 9 (2003), pp. 473–487.
  * [58] N. F. Rulkov, Regularization of synchronized chaotic bursts, Physical Review Letters, 86 (2001), p. 183.
  * [59] , Modeling of spiking-bursting neural behavior using two-dimensional map, Physical Review E, 65 (2002), p. 041922.
  * [60] J. A. Sanders, F. Verhulst, and J. Murdock, Averaging methods in nonlinear dynamical systems, vol. 59, Springer, 2007.
  * [61] S. Schecter and P. Szmolyan, Persistence of rarefactions under Dafermos regularization: blow-up and an exchange lemma for gain-of-stability turning points, SIAM Journal on Applied Dynamical Systems, 8 (2009), pp. 822–853.
  * [62] A. L. Shilnikov and N. F. Rulkov, Subthreshold oscillations in a map-based neuron model, Physics Letters A, 328 (2004), pp. 177–184.
  * [63] M. Shub, Global stability of dynamical systems, Springer Science & Business Media, 2013.
  * [64] M. Stiefenhofer, Singular perturbation with limit points in the fast dynamics, Zeitschrift für angewandte Mathematik und Physik ZAMP, 49 (1998), pp. 730–758.
  * [65] P. Szmolyan, Transversal heteroclinic and homoclinic orbits in singular perturbation problems, Journal of differential equations, 92 (1991), pp. 252–281.
  * [66] P. Szmolyan and M. Wechselberger, Relaxation oscillations in $\mathbb{R}^{3}$, Journal of Differential Equations, 200 (2004), pp. 69–104.
  * [67] A. N. Tikhonov, Systems of differential equations containing small parameters in the derivatives, Matematicheskii sbornik, 73 (1952), pp. 575–586.
  * [68] F. L. Trujillo, J. Signerska-Rynkowska, and P. Bartłomiejczyk, Periodic and chaotic dynamics in a map-based neuron model, arXiv preprint arXiv:2111.14499, (2021).
  * [69] B. van der Pol, A theory of the amplitude of free and forced triode vibrations, Radio Review, 1 (1920), pp. 701–710.
  * [70] , On “relaxation-oscillations”, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2 (1926), pp. 978–992.
  * [71] F. Wang and H. Cao, Mode locking and quasiperiodicity in a discrete-time chialvo neuron model, Communications in Nonlinear Science and Numerical Simulation, 56 (2018), pp. 481–489.
  * [72] M. Wechselberger, Geometric singular perturbation theory beyond the standard form, in Frontiers in Applied Dynamical Systems: Reviews and Tutorials, Springer International Publishing, 2020.
  * [73] S. Wiggins, Normally hyperbolic invariant manifolds in dynamical systems, vol. 105, Springer Science & Business Media, 1994.

The corresponding author S. Jelbart can be contacted via<EMAIL_ADDRESS>
C. Kuehn can be contacted via<EMAIL_ADDRESS>