can be analysed directly after a second positive transformation of time
$F(\rho_{3},\theta_{3},y_{3},\varepsilon_{3})\mathrm{d}t_{3}=\mathrm{d}\tilde{t}_{3}$,
which amounts to division of the right-hand side in system (54) by
$F(\rho_{3},\theta_{3},y_{3},\varepsilon_{3})$ (which is strictly positive in
$\mathcal{D}_{3}$). This leads to the orbitally equivalent system

$\displaystyle\rho_{3}^{\prime}$ $\displaystyle=\rho_{3},$ (55)
$\displaystyle\theta_{3}^{\prime}$
$\displaystyle=\frac{\rho_{3}^{\alpha-1}\varepsilon_{3}^{\alpha}}{b(\theta_{3})-a(\theta_{3})y_{3}+\mathcal{O}(\rho_{3})}=\frac{\rho_{3}^{\alpha-1}\varepsilon_{3}^{\alpha}}{b(\theta_{3})}+\mathcal{O}(\rho_{3}^{\alpha-1}\varepsilon_{3}^{\alpha}y_{3},\rho_{3}^{\alpha}\varepsilon_{3}^{\alpha}),$
$\displaystyle y_{3}^{\prime}$
$\displaystyle=-2y_{3}+\frac{\varepsilon_{3}^{3}(-c(\theta_{3})+\mathcal{O}(\rho_{3}))}{b(\theta_{3})-a(\theta_{3})y_{3}+\mathcal{O}(\rho_{3})}=-2y_{3}-\varepsilon_{3}^{3}\frac{c(\theta_{3})}{b(\theta_{3})}+\mathcal{O}(\varepsilon_{3}^{3}y_{3},\rho_{3}\varepsilon_{3}^{3}),$
$\displaystyle\varepsilon_{3}^{\prime}$ $\displaystyle=-\varepsilon_{3},$

where by a slight abuse of notation the prime notation now refers to
differentiation with respect to $\tilde{t}_{3}$.

System (55) is used to obtain the following result.

Figure 18: Geometry and dynamics within the invariant hyperplane
$\\{\rho_{3}=0\\}$ with $\alpha=2$, projected into
$(y_{3},\theta_{3},\varepsilon_{3})$-space; c.f. Figure 17.

.

Figure 19: Geometry and dynamics within the invariant hyperplane
$\\{\varepsilon_{3}=0\\}$ with $\alpha=2$, projected into
$(y_{3},\theta_{3},\rho_{3})$-space; c.f. Figures 17 and 19.

######  Proposition 4.20.

Fix $\nu>0$ and $E,R>0$ sufficiently small. Then the map
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\Pi_{3}^{(\alpha)}}}:\Sigma_{3}^{\textup{in}}\to\Sigma_{3}^{\textup{out}}$
is well-defined, at least $C^{1}$-smooth and given by

$\Pi_{3}^{(\alpha)}(\rho_{3},\theta_{3},y_{3},E)=\left(R,h_{\theta_{3}}^{(2)}(\rho_{3},\theta_{3},y_{3}),h_{y_{3}}^{(\alpha)}(\rho_{3},\theta_{3},y_{3}),\frac{E}{R}\rho_{3}\right).$

We have that
$h_{\theta_{3}}^{(\alpha)}(\rho_{3},\theta_{3},y_{3})=\tilde{h}_{\theta_{3}}^{(\alpha)}(\rho_{3},\theta_{3},y_{3})\mod{1}$,
where

$\tilde{h}_{\theta_{3}}^{(\alpha)}(\rho_{3},\theta_{3},y_{3})=\begin{cases}\theta_{3}+\mathcal{O}\left(\ln\left(\rho_{3}^{-1}\right)\right),&\alpha=1,\\\
\theta_{3}+\mathcal{O}\left({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\rho_{3}}}\ln\left(\rho_{3}^{-1}\right)\right),&\alpha=2.\end{cases}$

Moreover, there exists constants $\sigma_{-}<\sigma_{+}$ such that in case
$\alpha=1$ we have

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\left(-\nu+\sigma_{-}E^{3}\right)\left(1+\mathcal{O}(\rho_{3})\right)\left(\frac{\rho_{3}}{R}\right)^{2}\leq
h_{y_{3}}^{(1)}(\rho_{3},\theta_{3},y_{3})\leq\left(\nu+\sigma_{+}E^{3}\right)\left(1+\mathcal{O}(\rho_{3})\right)\left(\frac{\rho_{3}}{R}\right)^{2}.}}$
(56)

In case $\alpha=2$ we have

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{h_{y_{3}}^{(2)}}}(\rho_{3},\theta_{3},y_{3})=\left(y_{3}-E^{3}\right)\left(\frac{\rho_{3}}{R}\right)^{2}+\mathcal{O}\left(\rho_{3}^{3}\ln\left(\rho_{3}^{-1}\right)\right).$

###### Proof.

We start with the case $\alpha=2$, and consider solutions of system (55) which
satisfy
$(\rho_{3},\theta_{3},y_{3},\varepsilon_{3})(0)=(\rho_{3}^{\ast},\theta_{3}^{\ast},y_{3}^{\ast},E)\in\Sigma_{3}^{\textup{in}}$
and $(R,\theta_{3},y_{3},\varepsilon_{3})(T_{3})\in\Sigma_{3}^{\textup{out}}$
for some $T_{3}>0$. The equations for $\rho_{3}$ and $\varepsilon_{3}$ can be
solved directly, leading to
$\rho_{3}(\tilde{t}_{3})=\rho_{3}^{\ast}\mathrm{e}^{\tilde{t}_{3}}$ and
$\varepsilon_{3}(\tilde{t}_{3})=E\mathrm{e}^{-\tilde{t}_{3}}$, and the
boundary constraint $\rho_{3}(T_{3})=R$ leads to the following expression for
the transition time:

$T_{3}=\ln\left(\frac{R}{\rho_{3}^{\ast}}\right).$ (57)

It follows that $\varepsilon_{3}(T_{3})=E\rho_{3}^{\ast}/R$, as required. We
need to estimate $y_{3}(T_{3})$ and $\theta_{3}(T_{3})$ at
$\tilde{t}_{3}=T_{3}$.  Since $\alpha=2$, we have that

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\theta_{3}^{\prime}=\rho_{3}^{\ast}\left(E^{2}\mathrm{e}^{-\tilde{t}_{3}}(b(\theta_{3})-a(\theta_{3})y_{3}+\mathcal{O}(\rho_{3}))\right)^{-1}.}}$

One can show with direct estimates that the term in parentheses is uniformly
bounded by a positive constant, which implies that

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\rho_{3}^{\ast}\zeta_{-}}}\tilde{t}_{3}\leq\tilde{\theta}(\tilde{t}_{3})-\theta_{3}^{\ast}\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\rho_{3}^{\ast}\zeta_{+}}}\tilde{t}_{3},$

for all $\tilde{t}_{3}\in[0,T_{3}]$, where
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\zeta_{+}}}\geq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\zeta_{-}}}>0$
and $\tilde{\theta}_{3}$ is defined via
$\theta_{3}(\tilde{t}_{3})=\tilde{\theta}_{3}(\tilde{t}_{3})\mod 1$. Using
$\rho_{3}(\tilde{t}_{3})=\rho_{3}^{\ast}\mathrm{e}^{\tilde{t}_{3}}$, it
follows that $\tilde{\theta}_{3}$ can be written in terms of $\rho_{3}$,
specifically, we obtain

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\tilde{\theta}_{3}(\rho_{3})=\theta_{3}^{\ast}+\mathcal{O}\left(\rho_{3}^{\ast}\ln\left(\frac{\rho_{3}}{\rho_{3}^{\ast}}\right)\right),}}$
(58)

where we permit a slight abuse of notation in switching the argument from
$\tilde{t}_{3}$ to $\rho_{3}$. Substituting
$\rho_{3}(\tilde{t}_{3})=\rho_{3}^{\ast}\mathrm{e}^{\tilde{t}_{3}}$ and
evaluating the resulting expression at $\tilde{t}_{3}=T_{3}$ yields the
desired estimate for
$\theta_{3}(T_{3})=\tilde{h}_{\theta_{3}}^{(2)}(\rho_{3},\theta_{3},y_{3})\mod
1$.

It remains to estimate $y_{3}(T_{3})$. It will be helpful to use (58) in order
to write

$a(\theta_{3})=a_{3}+\mathcal{O}\left({\rho_{3}^{\ast}}\ln\left(\frac{\rho_{3}}{\rho_{3}^{\ast}}\right)\right),\
\
b(\theta_{3})=b_{3}+\mathcal{O}\left({\rho_{3}^{\ast}}\ln\left(\frac{\rho_{3}}{\rho_{3}^{\ast}}\right)\right),\
\
c(\theta_{3})=c_{3}+\mathcal{O}\left({\rho_{3}^{\ast}}\ln\left(\frac{\rho_{3}}{\rho_{3}^{\ast}}\right)\right),$
(59)

where $a_{3}=a(\theta_{3}^{\ast})$, $b_{3}=b(\theta_{3}^{\ast})$ and
$c_{3}=c(\theta_{3}^{\ast})$. From here, we use and adaptation of the ‘partial
linearisation’ method used in the proof of [39, Prop. 2.11]. Define a new
variable $\tilde{\varepsilon}_{3}=\varepsilon_{3}^{3}$, and consider the
resulting equations within the invariant hyperplane $\\{\rho_{3}=0\\}$:

$\begin{split}\theta_{3}^{\prime}&=0,\\\
y_{3}^{\prime}&=-2y_{3}-\frac{c_{3}}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{b_{3}}}}\tilde{\varepsilon}_{3}+\mathcal{O}(y_{3}\tilde{\varepsilon}_{3}),\\\
\tilde{\varepsilon}_{3}^{\prime}&=-3\tilde{\varepsilon}_{3},\end{split}$ (60)

which may be considered as a $\theta_{3}$-family of planar systems in the
$(y_{3},\tilde{\varepsilon}_{3})$ variables, with parameter-dependent
constants $a_{3}$, $b_{3}$ and $c_{3}$. Since the Jacobian at the equilibrium
$(0,0)$ is hyperbolic and non-resonant with the eigenvalues $-2$ and $-3$,
there exists a parameter-dependent, near-identity transformation of the form

$\tilde{y}_{3}=\psi(y_{3},\tilde{\varepsilon}_{3})=y_{3}+\mathcal{O}(y_{3}\tilde{\varepsilon}_{3}),\qquad
y_{3}=\tilde{\psi}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\tilde{y}_{3}}},\tilde{\varepsilon}_{3})=\tilde{y}_{3}+\mathcal{O}(\tilde{y}_{3}\tilde{\varepsilon}_{3}),$
(61)

such that the transformed system has been linearized, i.e.

$\begin{split}\tilde{y}_{3}^{\prime}&=-2\tilde{y}_{3}-\frac{c_{3}}{b_{3}}\tilde{\varepsilon}_{3},\\\
\tilde{\varepsilon}_{3}^{\prime}&=-3\tilde{\varepsilon}_{3}.\end{split}$ (62)

Applying the transformation (61) to system (55) and using (59) leads to

$\displaystyle\rho_{3}^{\prime}$ $\displaystyle=\rho_{3},$ (63)
$\displaystyle\theta_{3}^{\prime}$
$\displaystyle=\frac{\rho_{3}\varepsilon_{3}^{2}}{b(\theta_{3})}+\mathcal{O}(\rho_{3}\varepsilon_{3}^{2}\tilde{y}_{3},\rho_{3}^{2}\varepsilon_{3}^{2}),$
$\displaystyle\tilde{y}_{3}^{\prime}$
$\displaystyle=-2\tilde{y}_{3}-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\frac{c_{3}}{b_{3}}}}\varepsilon_{3}^{3}+\varepsilon_{3}^{3}\rho_{3}H(\rho_{3},\theta_{3},\tilde{y}_{3},\varepsilon_{3}),$
$\displaystyle\varepsilon_{3}^{\prime}$ $\displaystyle=-\varepsilon_{3},$

where the remainder function
$H(\rho_{3},\theta_{3},\tilde{y}_{3},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\varepsilon_{3}}})$
is smooth and uniformly bounded over $\tilde{t}_{3}\in[0,T_{3}]$. Substituting
the solutions for $\rho_{3}(\tilde{t}_{3})$ and
$\varepsilon_{3}(\tilde{t}_{3})$ into the equations for $\theta_{3}$ and
$\tilde{y}_{3}$ in system (63) leads to the planar non-autonomous system

$\displaystyle\theta_{3}^{\prime}$
$\displaystyle=\frac{{\rho_{3}^{\ast}}E^{2}}{b(\theta_{3})}\mathrm{e}^{-\tilde{t}_{3}}+{\rho_{3}^{\ast}}\mathcal{O}(\tilde{y}_{3}\mathrm{e}^{-\tilde{t}_{3}},\rho_{3}^{\ast}),$
$\displaystyle\tilde{y}_{3}^{\prime}$
$\displaystyle=-2\tilde{y}_{3}-\frac{c_{3}}{b_{3}}E^{3}\mathrm{e}^{-3\tilde{t}_{3}}+\rho_{3}^{\ast}E^{3}\mathrm{e}^{-2\tilde{t}_{3}}H(\rho_{3}^{\ast}\mathrm{e}^{\tilde{t}_{3}},\theta_{3},\tilde{y}_{3},E\mathrm{e}^{-\tilde{t}_{3}}).$

Now introduce a new time-dependent variable $z_{3}$ via

$\tilde{y}_{3}(\tilde{t}_{3})=\left(y_{3}^{\ast}-E^{3}+z_{3}(\tilde{t}_{3})\right)\mathrm{e}^{-2\tilde{t}_{3}}+\frac{c_{3}}{b_{3}}E^{3}\mathrm{e}^{-3\tilde{t}_{3}}.$

Differentiating $\tilde{y}_{3}(\tilde{t}_{3})$ and using
$\rho_{3}(\tilde{t}_{3})=\rho_{3}^{\ast}\mathrm{e}^{\tilde{t}_{3}}$ leads to
the following ODE in $z_{3}$:

$z_{3}^{\prime}=\rho_{3}^{\ast}E^{3}\tilde{H}\left(\rho_{3}^{\ast},\theta_{3},y_{3}^{\ast},z_{3},\tilde{t}_{3}\right),$

with
$\tilde{H}(\rho_{3}^{\ast},\theta_{3},y_{3}^{\ast},z_{3},\tilde{t}_{3})\coloneqq
H(\rho_{3}^{\ast}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathrm{e}^{\tilde{t}_{3}}}},\theta_{3},(y_{3}^{\ast}-E^{3}+z_{3}(\tilde{t}_{3}))\mathrm{e}^{-2\tilde{t}_{3}}+(c_{3}/b_{3})E^{3}\mathrm{e}^{-3\tilde{t}_{3}},E\mathrm{e}^{-\tilde{t}_{3}})$
and initial condition $z_{3}(0)=0$. One can show that
$\tilde{H}\left(\rho_{3}^{\ast},\theta_{3},y_{3}^{\ast},z_{3},\tilde{t}_{3}\right)$
is uniformly bounded over $\tilde{t}_{3}\in[0,T_{3}]$, so that integration
gives

$z_{3}(\tilde{t}_{3})=\mathcal{O}(\rho_{3}^{\ast}\tilde{t}_{3})\quad\implies\quad
z_{3}({T}_{3})=\mathcal{O}\left(\rho_{3}^{\ast}\ln\left(\frac{1}{\rho_{3}^{\ast}}\right)\right).$

Changing variables back to $\tilde{y}_{3}$ and evaluating at
$\tilde{t}_{3}=T_{3}$ yields

$\tilde{y}_{3}({T}_{3})=(y_{3}^{\ast}-E^{3}+z_{3}({T}_{3}))\mathrm{e}^{-2{T}_{3}}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\frac{c_{3}}{b_{3}}}}E^{3}\mathrm{e}^{-3{T}_{3}}\\\
=(y_{3}^{\ast}-E^{3})\left(\frac{\rho_{3}^{\ast}}{R}\right)^{2}+\mathcal{O}\left((\rho_{3}^{\ast})^{3}\ln\left(\frac{1}{\rho_{3}^{\ast}}\right)\right),$

and therefore

$y_{3}(T_{3})=(y_{3}^{\ast}-E^{3})\left(\frac{\rho_{3}^{\ast}}{R}\right)^{2}+\mathcal{O}\left((\rho_{3}^{\ast})^{3}\ln\left(\frac{1}{\rho_{3}^{\ast}}\right)\right)$

as required.

Now fix $\alpha=1$. The fact that
$\tilde{\theta}_{3}(T_{3})=\mathcal{O}(\ln{\rho_{3}^{\ast}}^{-1})$ follows
from direct estimates using the boundedness of the right-hand side of the
equation for $\theta_{3}^{\prime}$ and the expression for the transition time
$T_{3}$ in (57). It remains to bound $y_{3}(T_{3})$. Note that we cannot apply
the same partial linearisation approach as we did for $\alpha=2$, since
$\theta_{3}$ is dynamic (i.e. not constant) in $\\{\rho_{3}=0\\}$. We can,
however, obtain a coarser estimate via a more direct approach. Notice that for
sufficiently small $E$ we have $y_{3}^{\prime}|_{y_{3}=-\nu}>0$ and
$y_{3}^{\prime}|_{y_{3}=\nu}<0$. It follows that $\Pi_{3}^{(1)}$ is well-
defined, since solutions are inflowing along the faces of $\mathcal{D}_{3}$
defined by $\\{y=\pm\nu\\}$ for all $\tilde{t}_{3}\in[0,T_{3}]$. Moreover,
there exists constants $\sigma_{-}<\sigma_{+}$ such that

$-2y_{3}+\sigma_{-}\varepsilon_{3}^{3}\leq
y_{3}^{\prime}\leq-2y_{3}+\sigma_{+}\varepsilon_{3}^{3},$

for all $y\in[-\nu,\nu]$ and $\varepsilon_{3}\in[0,E]$. Substituting
$\varepsilon_{3}(\tilde{t}_{3})=E\mathrm{e}^{-\tilde{t}_{3}}$ and solving the
first order linear equations

${y_{3}^{\pm}}^{\prime}+2y_{3}^{\pm}=\sigma_{\pm}E^{3}\mathrm{e}^{-3\tilde{t}_{3}},$

we obtain lower and upper solutions on $\tilde{t}_{3}\in[0,T_{3}]$ with the
property that

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{y_{3}^{-}(\tilde{t}_{3})=\left(-\nu+\sigma_{-}E^{3}\right)\mathrm{e}^{-2\tilde{t}_{3}}-\sigma_{-}E^{3}\mathrm{e}^{-3\tilde{t}_{3}}\leq
y_{3}(\tilde{t}_{3})\leq
y_{3}^{+}(\tilde{t}_{3})=\left(\nu+\sigma_{+}E^{3}\right)\mathrm{e}^{-2\tilde{t}_{3}}-\sigma_{+}E^{3}\mathrm{e}^{-3\tilde{t}_{3}}.}}$

Evaluating these expressions at $\tilde{t}_{3}=T_{3}=\ln(R/\rho_{3}^{\ast})$
yields

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\left(-\nu+\sigma_{-}E^{3}\right)\left(1+\mathcal{O}(\rho_{3}^{\ast})\right)\left(\frac{\rho_{3}^{\ast}}{R}\right)^{2}\leq
y_{3}(T_{3})\leq\left(\nu+\sigma_{+}E^{3}\right)\left(1+\mathcal{O}(\rho_{3}^{\ast})\right)\left(\frac{\rho_{3}^{\ast}}{R}\right)^{2}.}}\qed$

Note that systems (54) and (55) are orbitally equivalent. Since the result in
Proposition 4.20 only depends on the initial condition in
$\Sigma_{3}^{\textup{in}}$, and not on the transition time, it holds for both
systems (54) and (55).

######  Remark 4.21.

Similarly to the arguments used to describe the map
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\Pi_{1}^{(\alpha)}}}$
in proof of Proposition 4.9, the arguments used to describe the map
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\Pi_{3}^{(\alpha)}}}$
in Proposition 4.20 are (strictly speaking) only valid for initial conditions
with $\rho_{3}\in(0,R]$. This is not problematic since we aim to derive
results for $\varepsilon>0$. Moreover,
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\Pi_{3}^{(\alpha)}}}$
can be smoothly and uniquely extended to $\rho_{3}=0$.

###  4.5 Proof of Theorem 3.2

We now combine the results obtained in the charts $K_{1}$, $K_{2}$ and $K_{3}$
in order to prove the Theorem 3.2. The idea is to describe the extended map
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\pi_{\varepsilon}^{(\alpha)}}}:\Delta_{\varepsilon}^{\textup{in}}\rightarrow\Delta_{\varepsilon}^{\textup{out}}$
defined in Section 4.1, the first three components of which coincide with
components of
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\pi^{(\alpha)}}}:\Delta^{\textup{in}}\rightarrow\Delta^{\textup{out}}$,
via its preimage in the blown-up space:

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\Pi^{(\alpha)}\coloneqq\Pi_{3}^{(\alpha)}\circ\kappa_{23}\circ\Pi_{2}^{(\alpha)}\circ\kappa_{12}\circ\Pi_{1}^{(\alpha)}:\Sigma_{1}^{\textup{in}}\to\Sigma_{3}^{\textup{out}}.}}$

This can be done explicitly using the change of coordinate maps $\kappa_{ij}$
in Lemma 4.1 and the characterisation of $\Pi_{1}^{(\alpha)}$,
$\Pi_{2}^{(\alpha)}$ and $\Pi_{3}^{(\alpha)}$ in Propositions 4.9, 4.15 and
4.20 respectively. Note that these results immediately imply that
$\Pi^{(\alpha)}$ is at least $C^{1}$-smooth, since each $\kappa_{ij}$ and
$\Pi_{i}^{(\alpha)}$ is at least $C^{1}$-smooth. The geometry is sketched in
Figure 20.

In order to prove Assertions (a)-(c), we need to derive the form of the map
$\Pi^{(\alpha)}$. Initial conditions
$(r_{1},\theta_{1},R,\varepsilon_{1})\in\Sigma_{1}^{\textup{in}}$ are mapped
to $\Sigma_{1}^{\textup{out}}$ under $\Pi_{1}^{(\alpha)}$ as described by
Proposition 4.9. Since
$\Sigma_{2}^{\textup{in}}=\kappa_{12}(\Sigma_{1}^{\textup{out}})$, we obtain
the following input for the map $\Pi_{2}^{(\alpha)}$:

$\displaystyle(r_{2},\theta_{2},y_{2},\rho_{2})$
$\displaystyle=\kappa_{12}\left(\Pi_{1}^{(\alpha)}(r_{1},\theta_{1},R,\varepsilon_{1})\right)$
$\displaystyle=\left(\frac{1}{E}\left(h_{r_{1}}^{(\alpha)}\left(h_{\theta_{1}}^{(\alpha)}(r_{1},\theta_{1},\varepsilon_{1}),\frac{R}{E}\varepsilon_{1},E\right)+\mathcal{O}\left(\mathrm{e}^{-\tilde{\varrho}/\varepsilon_{1}^{3}}\right)\right),h_{\theta_{1}}^{(\alpha)}(r_{1},\theta_{1},\varepsilon_{1}),\frac{1}{E^{2}},R\varepsilon_{1}\right),$
(64)

where we recall that
$h_{\theta_{1}}^{(\alpha)}(r_{1},\theta_{1},\varepsilon_{1})=\tilde{h}_{\theta_{1}}^{(\alpha)}(r_{1},\theta_{1},\varepsilon_{1})\mod{1}$,

$\tilde{h}_{\theta_{1}}^{(\alpha)}(r_{1},\theta_{1},\varepsilon_{1})=\theta_{1}+\frac{(R\varepsilon_{1})^{\alpha-1}}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{c_{0}}}}\left(\frac{1}{{\varepsilon_{1}}^{2}}+\mathcal{O}(1)\right),$
(65)

and $c_{0}\coloneqq\int_{0}^{1}c(\theta_{1}(t_{1}))\,\mathrm{d}t_{1}$. Using
(4.5), Lemma 4.1 and Proposition 4.15, we obtain an expression for the input
to the map $\Pi_{3}^{(\alpha)}$, namely

$\displaystyle(\rho_{3},\theta_{3},y_{3},\varepsilon_{3})$
$\displaystyle=\kappa_{23}\left(\Pi_{2}^{(\alpha)}(r_{2},\theta_{2},y_{2},\rho_{2})\right)$
$\displaystyle=\left(\frac{R}{E}\varepsilon_{1},h_{\theta_{2}}^{(\alpha)}\left(r_{2},\theta_{2},\rho_{2}\right),E^{2}h_{y_{2}}^{(\alpha)}(E^{-1},\theta_{2})+\mathcal{O}(\varepsilon_{1}),E\right),$
(66)

where we used the fact that
$r_{2}-r_{2,c}=\mathcal{O}(\mathrm{e}^{-\tilde{\varrho}/\varepsilon_{1}^{3}})$
in the derivation of the third component in the right-hand side, and
$h_{\theta_{2}}^{(\alpha)}(r_{2},\theta_{2},\rho_{2})=\tilde{h}_{\theta_{2}}^{(\alpha)}(r_{2},\theta_{2},\rho_{2})\mod
1$ with

$\tilde{h}_{\theta_{2}}^{(\alpha)}(r_{2},\theta_{2},\rho_{2})=\begin{cases}\tilde{h}_{\theta_{1}}^{(1)}(r_{1},\theta_{1},\varepsilon_{1})+\mathcal{O}(1),&\alpha=1,\\\
\tilde{h}_{\theta_{1}}^{(2)}(r_{1},\theta_{1},\varepsilon_{1})+\frac{R\varepsilon_{1}}{c(\theta_{2})}\left(E^{-2}+\left(\frac{c^{2}}{ab}\right)^{1/3}\Omega_{0}-E\right)+\mathcal{O}\left(\varepsilon_{1}^{2}\right),&\alpha=2,\end{cases}$
(67)

where
$c(\theta_{2})=c(h_{\theta_{1}}^{(2)}(r_{1},\theta_{1},\varepsilon_{1}))$.
With regards to the third component in (4.5), if $\alpha=1$ then we have

$-\nu\leq E^{2}h_{y_{2}}^{(1)}(E^{-1},\theta_{2})\leq\nu,$

where $\nu>0$ is the constant appearing in the definition of
$\Sigma_{2}^{\textup{out}}$ (recall Proposition 4.13). If $\alpha=2$, then by
Lemma 4.16 it holds

$h_{y_{2}}^{(2)}(E^{-1},\theta_{2})={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{-\left(\frac{c(\theta_{2})^{2}}{a(\theta_{2})b(\theta_{2})}\right)^{1/3}}}\Omega_{0}+E,$

where
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a(\theta_{2})}}=a\left(h_{\theta_{1}}^{(\alpha)}(r_{1},\theta_{1},\varepsilon_{1})\right)$,
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{b(\theta_{2})}}=b\left(h_{\theta_{1}}^{(\alpha)}(r_{1},\theta_{1},\varepsilon_{1})\right)$
and
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{c(\theta_{2})}}=c\left(h_{\theta_{1}}^{(\alpha)}(r_{1},\theta_{1},\varepsilon_{1})\right)$.

Substituting the right-most expression in (4.5) into the expression for
$\Pi_{3}^{(\alpha)}(\rho_{3},\theta_{3},y_{3},E)$ in Proposition 4.20 yields
an expression for $\Pi^{(\alpha)}(r_{1},\theta_{1},R,\varepsilon_{1})$. We
obtain

$\Pi^{(\alpha)}(r_{1},\theta_{1},R,\varepsilon_{1})=\left(R,h_{\theta_{3}}^{(\alpha)}(\rho_{3},\theta_{3},y_{3}),h_{y_{3}}^{(\alpha)}(\rho_{3},\theta_{3},y_{3}),\varepsilon_{1}\right),$
(68)

and specify the second and third components in the following. We have
$h_{\theta_{3}}^{(\alpha)}(\rho_{3},\theta_{3},y_{3})=\tilde{h}_{\theta_{3}}^{(\alpha)}(\rho_{3},\theta_{3},y_{3})\mod
1$, with

$\tilde{h}_{\theta_{3}}^{(\alpha)}(\rho_{3},\theta_{3},y_{3})=\begin{cases}\tilde{h}_{\theta_{2}}^{(1)}(r_{1},\theta_{1},\varepsilon_{1})+\mathcal{O}\left(\ln\rho_{3}^{-1}\right)=\theta_{1}+\frac{1}{c_{0}\varepsilon_{1}^{2}}+\mathcal{O}(\ln\varepsilon_{1}^{-1}),&\alpha=1,\\\
\tilde{h}_{\theta_{2}}^{(2)}(r_{1},\theta_{1},\varepsilon_{1})+\mathcal{O}\left(\rho_{3}\ln\rho_{3}^{-1}\right)=\theta_{1}+\frac{R}{c_{0}\varepsilon_{1}}+\mathcal{O}(\varepsilon_{1}\ln\varepsilon_{1}^{-1}),&\alpha=2,\end{cases}$

where we used (65) and (67) in order to obtain the right-hand side. It remains
to specify the third component in the right-hand side of (68). If $\alpha=1$,
we can use the bounds in (56) to estimate

$\left(-\nu+\sigma_{-}E^{3}\right)\left(\frac{\varepsilon_{1}}{E}\right)^{2}+\mathcal{O}\left(\varepsilon_{1}^{3}\right)\leq
h_{y_{3}}^{(1)}(\rho_{3},\theta_{3},y_{3})\leq\left(\nu+\sigma_{+}E^{3}\right)\left(\frac{\varepsilon_{1}}{E}\right)^{2}+\mathcal{O}\left(\varepsilon_{1}^{3}\right)$

for constants $\nu>0$ and $\sigma_{-}<\sigma_{+}$, such that
$h_{y_{3}}^{(1)}(\rho_{3},\theta_{3},y_{3})=\mathcal{O}(\varepsilon_{1}^{2})$.
If $\alpha=2$, then

$\displaystyle h_{y_{3}}^{(2)}(\rho_{3},\theta_{3},y_{3})$
$\displaystyle=\left(E^{2}h_{y_{2}}^{(2)}(E^{-1},\theta_{2})-E^{3}\right)\left(\frac{\varepsilon_{1}}{E}\right)^{2}+\mathcal{O}\left(\varepsilon_{1}^{3}\ln\varepsilon_{1}^{-1}\right)$
$\displaystyle=-\left(\frac{c(\theta_{3})^{2}}{a(\theta_{3})b(\theta_{3})}\right)^{1/3}\Omega_{0}\varepsilon_{1}^{2}+\mathcal{O}\left(\varepsilon_{1}^{3}\ln\varepsilon_{1}^{-1}\right),$

where

$\displaystyle a(\theta_{3})$
$\displaystyle=a\left(h_{\theta_{3}}^{(2)}(\rho_{3},\theta_{3},y_{3})\right)=a(\theta_{2})+\mathcal{O}(\varepsilon_{1}\ln\varepsilon_{1}^{-1}),$
$\displaystyle b(\theta_{3})$
$\displaystyle=b\left(h_{\theta_{3}}^{(2)}(\rho_{3},\theta_{3},y_{3})\right)=b(\theta_{2})+\mathcal{O}(\varepsilon_{1}\ln\varepsilon_{1}^{-1}),$
$\displaystyle c(\theta_{3})$
$\displaystyle=c\left(h_{\theta_{3}}^{(2)}(\rho_{3},\theta_{3},y_{3})\right)=c(\theta_{2})+\mathcal{O}(\varepsilon_{1}\ln\varepsilon_{1}^{-1}),$

which follows from
$a(\theta_{2})=a(h_{\theta_{1}}^{(2)}(r_{1},\theta_{1},\varepsilon_{1}))$,
$b(\theta_{2})=b(h_{\theta_{1}}^{(2)}(r_{1},\theta_{1},\varepsilon_{1}))$,
$c(\theta_{2})=c(h_{\theta_{1}}^{(2)}(r_{1},\theta_{1},\varepsilon_{1}))$
together with
$\tilde{h}_{\theta_{3}}^{(2)}(\rho_{3},\theta_{3},y_{3})=\tilde{h}_{\theta_{1}}^{(2)}(r_{1},\theta_{1},\varepsilon_{1})+\mathcal{O}(\varepsilon_{1}\ln\varepsilon_{1}^{-1})$.

Applying the blow-down transformations defined by the local coordinate
formulae in (31) to (68) yields the map

$\pi_{\varepsilon}^{(\alpha)}(r,\theta,R^{2},\varepsilon)=\left(R,h_{\theta}^{(\alpha)}(r,\theta,\varepsilon),h_{y}^{(\alpha)}(r,\theta,\varepsilon),\varepsilon\right),$
(69)

for all $\varepsilon\in(0,\varepsilon_{0}]$ (the limiting value
$\varepsilon=0$ is omitted because the blow-down transformation is only
defined for $\varepsilon>0$). The estimates above ‘blow-down’ to those given
in Assertion (b) in Theorem 3.2, since the map $\pi$ is identified with
$\pi_{\varepsilon}$ after omitting the trivial component
$\varepsilon\mapsto\varepsilon$. Assertion (a) is a direct consequence of
Assertion (b).

Figure 20: Global geometry and dynamics in the blown-up space projected into
$(\bar{r},\bar{y},\bar{\varepsilon})$-space. The three-dimensional manifold
$M^{\textup{a}}$ which extends from $\Delta_{\varepsilon}^{\textup{in}}$ up to
$\Delta_{\varepsilon}^{\textup{out}}$ (which are identified with
$\Sigma^{\textup{in}}_{1}$ and $\Sigma^{\textup{out}}_{3}$ respectively) is
shown again in shaded green. Theorem 3.2 is proven by combining results
obtained in charts $K_{1}$, $K_{2}$ and $K_{3}$: The flow from
$\Sigma_{1}^{\textup{in}}$ to $\Sigma_{1}^{\textup{out}}$ is described by
Proposition 4.9, the flow from $\Sigma_{2}^{\textup{in}}$ to
$\Sigma_{2}^{\textup{out}}$ is described by Proposition 4.15 and the flow from
$\Sigma_{3}^{\textup{in}}$ to $\Sigma_{3}^{\textup{out}}$ is described by
Proposition 4.20.

It remains to prove the strong contraction property in Assertion (c). This can
be done by differentiating the expression for $\Pi^{(\alpha)}$. Letting
$\mathcal{Q}^{(\alpha)}\coloneqq\Pi_{3}^{(\alpha)}\circ\kappa_{23}\circ\Pi_{2}^{(\alpha)}\circ\kappa_{12}$,
we obtain

$\frac{\partial\Pi^{(\alpha)}}{\partial
r_{1}}(r_{1},\theta_{1},R,\varepsilon_{1})=\frac{\partial\mathcal{Q}^{(\alpha)}}{\partial
r_{1}}\left(\Pi_{1}^{(\alpha)}(r_{1},\theta_{1},R,\varepsilon_{1})\right)\frac{\partial\Pi_{1}^{(\alpha)}}{\partial
r_{1}}(r_{1},\theta_{1},R,\varepsilon_{1}),$

after applying the chain rule. Direct calculations using Lemma 4.1,
Propositions 4.15 and 4.20 show that $(\partial\mathcal{Q}^{(\alpha)}/\partial
r_{1})(\Pi_{1}^{(\alpha)}(r_{1},\theta_{1},R,\varepsilon_{1}))$ is (at worst)
algebraically growing in $r_{1}$ and $\varepsilon_{1}$. Since the second term
$(\partial\Pi_{1}^{(\alpha)}/\partial
r_{1})(r_{1},\theta_{1},R,\varepsilon_{1})$ is exponentially small due to the
strong contraction in $K_{1}$, recall Proposition 4.9 Assertion (b), it
follows that

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\frac{\partial\Pi^{(\alpha)}}{\partial
r_{1}}(r_{1},\theta_{1},R,\varepsilon_{1})=\mathcal{O}\left(\mathrm{e}^{-\tilde{\varrho}/\varepsilon_{1}^{3}}\right),}}$

for a constant $\tilde{\varrho}>0$, assuming the constant $E$ bounding
$\varepsilon_{1}$ is sufficiently small. Assertion (c) in Theorem 3.2 follows
for all $\varepsilon\in(0,\varepsilon_{0}]$ with $\varepsilon_{0}$
sufficiently small after applying the blow-down transformation (in particular
$\varepsilon_{1}=\varepsilon/R$). ∎

##  5 Applications

In this section we consider semi-oscillatory dynamics near folded limit cycle
manifolds in two different applications. Specifically, we consider

  1. 1.

Periodically forced oscillators of Liénard type;

  2. 2.

A toy model ‘normal form’ for the study of tipping phenomena in climate
systems proposed in [63].

The first of these will allow us to illustrate the application of the main
results in Section 3. The second is an example of a model exhibiting a folded
limit cycle bifurcation in the partial singular limit $\varepsilon_{1}>0$,
$\varepsilon_{2}\to 0$ (recall system (5)), which violates both Assumptions
1-2.

###  5.1 Periodically forced oscillators

Periodically forced oscillators in mechanical and electrical engineering can
often be modelled using Liénard equations of the form

$x^{\prime\prime}+\mu f(x)x^{\prime}+g(x)=A(\omega t),$ (70)

where $A(\omega t)=A(\omega t+1)$ for all $t\in\mathbb{R}$, the
friction/resistence function $f(x)$ is smooth, and the parameters $\mu$ and
$\omega$ determine the magnitude of the friction/resistance and the frequency
of the (periodic) forcing term $A$, respectively.

######  Remark 5.1.

A very well-known and historically important example of a Liénard equation is
the van der Pol equation [58, 59]. We refer to [8, 21] for detailed geometric
studies of the (local and global) dynamics of the periodically forced vdP
equation in a number of frequency regimes.

We recast equation (70) as an autonomous dynamical system by (i) rescaling
time via $\tau=\mu^{-1}t$, and (ii) introducing new variables via the Liénard
transformation

$y\coloneqq-\mu^{-1}x^{\prime}+K(x),\qquad s=\omega\mu\tau,$

where $K(x)=-\int_{0}^{x}f(\xi)\,\mathrm{d}\xi$. This leads to

$\begin{split}\dot{x}&=\mu^{2}\left(-y+K(x)\right),\\\ \dot{s}&=\mu\omega,\\\
\dot{y}&=g(x)-A(s),\end{split}$

where the overdot denotes now differentiation with respect to $\tau$.
Exploiting the periodicity property $A(s+1)=A(s)$ in order to replace $s$ with
an angular variable $\theta\in\mathbb{R}/\mathbb{Z}$ and rewriting the system
on the (new) fast time-scale defined by $\tilde{t}=\mu^{2}\tau$ leads to

$\begin{split}x^{\prime}&=-y+K(x),\\\ \theta^{\prime}&=\mu^{-1}\omega,\\\
y^{\prime}&=\mu^{-2}\left(g(x)-A(\theta)\right),\end{split}$

where by a slight abuse of notation, we allow the dash to denote
differentiation with respect to the new time $\tilde{t}$. Our analysis applies
in particular scaling regimes with
$0<{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mu^{-1}\omega}},\mu^{-2}\ll
1$. Specifically, we derive results for scaling regimes defined by

$\varepsilon^{\alpha}=\mu^{-1}\omega\vartheta^{-1},\qquad\varepsilon^{3}=\mu^{-2},$
(71)

where $\alpha\in\mathbb{N}_{+}$ and $\vartheta>0$ is a parameter such that
$\vartheta=\mathcal{O}(1)$ as $\varepsilon\to 0$. Note that this is equivalent
to considering the scaling regimes defined by

$\omega=\vartheta\mu^{1-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\tfrac{2}{3}}}\alpha},\qquad\alpha=1,2,\ldots,\qquad
0<\mu^{-1}\ll 1.$ (72)

After one more change of time-scale $\hat{t}=\vartheta\tilde{t}$, our system
becomes

$\begin{split}x^{\prime}&=-\vartheta^{-1}y+\vartheta^{-1}K(x),\\\
\theta^{\prime}&=\varepsilon^{\alpha},\\\
y^{\prime}&=\varepsilon^{3}\left(\vartheta^{-1}g(x)-\vartheta^{-1}A(\theta)\right),\end{split}$

where by a slight abuse of notation, we now allow the dash to denote
differentiation with respect to the time $\hat{t}$. Finally, we assume the
following properties of the functions $K$, $g$ and $A$.

######  Assumption 3.

There exists a point $x_{F}$ such that

$\frac{\partial K}{\partial x}(x_{F})=0,\qquad\frac{\partial^{2}K}{\partial
x^{2}}(x_{F})>0,\qquad g(x_{F})-A(\theta)<0,$

for all $\theta\in\mathbb{R}/\mathbb{Z}$.

Assumption 3 implies that the critical manifold, given by
$S_{0}=\\{(x,\theta,K(x)):x\in\mathbb{R},\theta\in\mathbb{R}/\mathbb{Z}\\}$,
contains a regular folded cycle
$S_{0}^{\textup{c}}=\\{(x_{F},\theta,K(x_{F})):\theta\in\mathbb{R}/\mathbb{Z}\\}$.
The set $S_{0}^{\textup{c}}$ can be moved to the ‘origin’ via the coordinate
translation $(\tilde{x},\tilde{y})=(x-x_{F},y-K(x_{F}))$, which yields the
system

$\begin{split}\tilde{x}^{\prime}&=-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a(\theta)\tilde{y}}}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{b(\theta)}}\tilde{x}^{2}+\mathcal{O}(\tilde{x}^{3}),\\\
\theta^{\prime}&=\varepsilon^{\alpha},\\\
\tilde{y}^{\prime}&=\varepsilon^{3}\left(-c(\theta)+\mathcal{O}(\tilde{x})\right),\end{split}$
(73)

which is – up to relabelling $(\tilde{x},\tilde{y})\leftrightarrow(r,y)$ – in
the general form (22) with

$a(\theta)\equiv\vartheta^{-1},\qquad
b(\theta)\equiv\vartheta^{-1}\frac{1}{2}\frac{\partial^{2}K}{\partial
x^{2}}(x_{F})=:\vartheta^{-1}\frac{1}{2}K^{\prime\prime}(x_{F}),\qquad
c(\theta)=-\vartheta^{-1}\left(g(x_{F})-A(\theta)\right),$ (74)

all of which are strictly positive due to Assumption 3. Using the same
relabelling, we define entry and exit sections $\Delta^{\textup{in}}$ and
$\Delta^{\textup{out}}$ analogously to those in (23) and (24), respectively,
and let $\pi^{(\alpha)}:\Delta^{\textup{in}}\rightarrow\Delta^{\textup{out}}$
denote the transition map induced by the forward flow of system (73) for a
given $\alpha\in\mathbb{N}_{+}$.

######  Theorem 5.2.

Consider system (73) with fixed $\alpha\in\mathbb{N}_{+}$. There exists an
$\varepsilon_{0}>0$ such that for all $\varepsilon\in(0,\varepsilon_{0}]$, the
map $\pi^{(\alpha)}:\Delta^{\textup{in}}\rightarrow\Delta^{\textup{out}}$ is
described by Theorem 3.2 after relabelling
$(\tilde{x},\tilde{y})\leftrightarrow(r,y)$. In particular, we have

$\tilde{h}_{\theta}^{(\alpha)}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\tilde{x}}},\theta,\varepsilon)=\begin{cases}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\theta+\frac{R^{2}}{c_{0}}\varepsilon^{-2}+\mathcal{O}(\ln\varepsilon)}},&\alpha=1,\\\
\theta+\frac{R^{2}}{c_{0}}\varepsilon^{-1}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathcal{O}(\varepsilon\ln\varepsilon)}},&\alpha=2,\\\
\psi(\theta)+\mathcal{O}(\varepsilon^{3}\ln\varepsilon),&\alpha=3,\\\
\theta+\mathcal{O}(\varepsilon^{3}\ln\varepsilon),&\alpha\geq 4,\end{cases}$

where

$c_{0}=-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\vartheta}}^{-1}g(x_{F})+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\vartheta}}^{-1}\int_{0}^{1}A(\theta)\mathrm{d}\theta,\qquad\psi(\theta)=\theta+\frac{1}{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\vartheta}}\frac{K^{\prime\prime}(x_{F})}{A(\theta)-g(x_{F})}R+\mathcal{O}(R^{2}),$

and

$h_{\tilde{y}}^{(\alpha)}(\tilde{\theta},\varepsilon)=\begin{cases}\mathcal{O}(\varepsilon^{2}),&\alpha=1,\\\
-\left(\frac{2(g(x_{F})-A(\tilde{\theta}))^{2}}{K^{\prime\prime}(x_{F})}\right)^{1/3}\Omega_{0}\varepsilon^{2}+\mathcal{O}(\varepsilon^{3}\ln\varepsilon),&\alpha\geq
2.\end{cases}$

###### Proof.

This is a direct application of Theorem 3.2, using the notation of (73) and
the identities in (74). ∎

Note that via (71), the results in Theorem 5.2 can be reformulated and
restated for the scaling regimes defined by (72). As for our main results in
Section 3, the results for scalings corresponding to semi-oscillatory dynamics
($\alpha\in\\{1,2\\}$) cannot be obtained using established theory for two
time-scale systems. These results are new, to the best of our knowledge, even
for the forced van der Pol example (the intermediate frequency regime
considered in [8] covers the case $\alpha=3/2$).

###  5.2 Tipping phenomena in climate systems

The authors in [63] proposed the following simple normal form for studying
early warning signs for saddle-node induced tipping in climate systems:

$z^{\prime}(t)=\tilde{a}(t)-z^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{-}}\mathcal{A}\sin(2\pi\omega
t),\qquad\tilde{a}(t)=\tilde{a}(0)-\tilde{\delta}t,$

where $0<\tilde{\delta}\ll 1$ and the parameters $\mathcal{A},\omega>0$
determine the amplitude, frequency of the sinusoidal periodic forcing term
respectively. This can recast as an autonomous system of the form

$\begin{split}z^{\prime}&=\tilde{a}-z^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{-}}\mathcal{A}\sin(2\pi\theta),\\\
\theta^{\prime}&=\omega,\\\ \tilde{a}^{\prime}&=-\tilde{\delta},\end{split}$
(75)

where $\theta\in\mathbb{R}/\mathbb{Z}$. We shall be interested in the dynamics
when $0<\tilde{\delta},\omega\ll 1$, and in particular the semi-oscillatory
case $0<\tilde{\delta}\ll\omega\ll 1$. In this case, system (75) is
oscillatory with respect to the partial singular limit $\omega>0$,
$\tilde{\delta}\to 0$, but stationary with respect to the double singular
limit $(\omega,\tilde{\delta})\to(0,0)$.

We are particularly interested in system (75) as a model example of a system
which violates Assumptions 1-2.

######  Proposition 5.3.

Consider system (75) with $\omega>0$ fixed and $0<\tilde{\delta}\ll 1$. This
system has a regular fold of limit cycles at $z=\tilde{a}=0$. In particular,
the necessary and sufficient conditions on the Poincaré map in (9) and (10)
are satisfied, however Assumptions 1 and 2 are not.

###### Proof.

The fact that Assumptions 1-2 are violated follows from direct calculations on
the system obtained from (75) in the partial singular limit $\omega>0$,
$\tilde{\delta}\to 0$. The details are omitted for brevity.

In order to verify the existence of a regular fold of cycles, we derive an
expression for the Poincaré map $P$ induced on a section contained in the
plane $\\{\theta=0\\}$. Using the fact that

$\frac{\mathrm{d}z}{\mathrm{d}\theta}=\omega^{-1}\left(\tilde{a}-z^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{-}}\mathcal{A}\sin(2\pi\theta)\right),\qquad\frac{\mathrm{d}\tilde{a}}{\mathrm{d}\theta}=-\omega^{-1}\tilde{\delta},$

we obtain

$P(z,\tilde{a},\omega,\tilde{\delta})=\begin{pmatrix}P_{z}(z,\tilde{a},\omega,\tilde{\delta})\\\
P_{\tilde{a}}(z,\tilde{a},\omega,\tilde{\delta})\end{pmatrix}=\begin{pmatrix}z\\\
\tilde{a}\end{pmatrix}+\omega^{-1}\begin{pmatrix}\int_{0}^{1}\left(\tilde{a}(\theta)-z(\theta)^{2}\right)\mathrm{d}\theta\\\
-\tilde{\delta}\end{pmatrix}.$

This expression can be used to check the necessary and sufficient conditions
in (9) and (10) for a regular fold of cycles directly. ∎

Our aim in what follows is to demonstrate that problems of this kind can (to a
large extent) be analysed with existing theory for two time-scale systems. We
shall restrict attention to scaling regimes defined by

$\tilde{\delta}=\varepsilon^{3}\nu,\qquad\omega=\varepsilon^{\alpha},$

where $\nu=\mathcal{O}(1)$ as $\varepsilon\to 0$ and
$\alpha\in\mathbb{N}_{+}$. We shall also simplify the fast equation by
defining a new variable
$a\coloneqq\tilde{a}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{-}}\mathcal{A}\sin(2\pi\theta)$.
This leads to

$\begin{split}z^{\prime}&=a-z^{2},\\\
\theta^{\prime}&=\varepsilon^{\alpha},\\\
a^{\prime}&=-2\pi\varepsilon^{\alpha}\mathcal{A}\cos(2\pi\theta)-\varepsilon^{3}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>(76)

The main results for this section will describe the dynamics of system (76).
Letting $\varepsilon\to 0$ yields the layer problem

$\begin{split}z^{\prime}&=a-z^{2},\\\ \theta^{\prime}&=0,\\\
a^{\prime}&=0.\end{split}$ (77)

The critical manifold is given by
$S_{0}=\\{(z,\theta,z^{2}):z\in\mathbb{R},\theta\in\mathbb{R}/\mathbb{Z}\\}$,
and the non-trivial eigenvalue of the linearisation along $S_{0}$ is
$\lambda=-2z$. This implies a folded critical manifold structure
$S_{0}=S_{0}^{\textup{a}}\cup S_{0}^{\textup{c}}\cup S_{0}^{\textup{r}}$,
where $S_{0}^{\textup{a}}=S_{0}\cap\\{z>0\\}$
($S_{0}^{\textup{r}}=S_{0}\cap\\{z<0\\}$) is normally hyperbolic and
attracting (repelling), and
$S_{0}^{\textup{c}}=\\{(0,\theta,0):\theta\in\mathbb{R}/\mathbb{Z}\\}$ is of
regular fold type.

The reduced problem on $S_{0}$ will differ depending on the scaling, i.e.
depending on the value of $\alpha$. We are primarily interested in the cases
$\alpha\in\\{1,2\\}$, for which
$0<{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\tilde{\delta}\ll\omega}}\ll
1$; see however Remark 5.8 below for the cases $\alpha\geq 3$. Fixing
$\alpha\in\\{1,2\\}$, rewriting system (76) on the slow time-scale
$\tau=\varepsilon^{\alpha}t$, and taking the limit $\varepsilon\to 0$ yields
the reduced problem

$\displaystyle\begin{split}0&=a-z^{2},\\\ \dot{\theta}&=1,\\\
\dot{a}&=-2\pi\mathcal{A}\cos(2\pi\theta).\end{split}$ (78)

This leads to the following reduced vector field on $S_{0}$, expressed in the
$(z,\theta)$-coordinate chart:

$\begin{split}\dot{z}&=-\frac{\pi\mathcal{A}}{z}\cos(2\pi\theta),\\\
\dot{\theta}&=1.\end{split}$ (79)

It is typical for problems of this kind to consider the so-called
desingularised reduced problem [42, 55, 61]. This may be obtained after a
(singular) time transformation which amounts to multiplication of the vector
field by $z$ (see again Remark 2.2), leading to

$\begin{split}\dot{z}&=-\pi\mathcal{A}\cos(2\pi\theta),\\\
\dot{\theta}&=z.\end{split}$ (80)

System (80) is orbitally equivalent to system (79) on $S_{0}\setminus
S_{0}^{\textup{c}}$, but the orientation of orbits on $S_{0}^{\textup{r}}$
(where $z<0)$ is reversed. Direct calculations reveal the presence of two
equilibria along $S_{0}^{\textup{c}}$, namely

$p_{s}:(0,1/4),\qquad p_{c}:(0,3/4),$

which are of neutral folded saddle and folded center type respectively (see
[55] for definitions), see Figure 21. Moreover, system (80) is Hamiltonian,
with Hamiltonian function

$H(z,\theta)=\frac{\mathcal{A}}{2}\sin(2\pi\theta)+\frac{z^{2}}{2}.$

Except for the points $p_{s}$ and $p_{c}$, there are three different types of
orbits shown in Figure 21. These are identified in the following result.

Figure 21: Dynamics for the reduced problem (79) on $S_{0}$, projected into
the $(\theta,z)$-plane. We sketch the case with $\mathcal{A}=2$ and $R=1$,
where we recall that $R>0$ is the constant used to define
$\Delta^{\textup{in}}$. The normally hyperbolic and attracting (repelling)
branch $S_{0}^{\textup{a}}$ ($S_{0}^{\textup{r}}$) is indicated in shaded
turquoise (gray). The fold cycle $S_{0}^{\textup{c}}$ which separates these
two branches is sketched in blue. There are two folded singularities on
$S_{0}^{\textup{c}}$: a neutral folded saddle $p_{s}$, and a folded center
$p_{c}$. The singular true and faux canards through $p_{s}$ coincide in this
case (denoted $\gamma_{c}$ and shown in green), forming a limit cycle of
period $2$, and the region bounded outside of $\gamma_{c}$ is foliated by
periodic orbits of period $1$. Finally, we sketch the forward evolution of an
initial condition in
$\Delta^{\textup{in}}\cap{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{S_{0}^{\textup{a}}}}$
with initial angle $\theta^{\ast}\in\mathcal{I}$ (the interval $\mathcal{I}$
is shown in blue), i.e. inside the region bounded by $\gamma_{c}$. Such points
reach a regular jump point on $S_{0}^{\textup{c}}$ with angle
$\theta_{e}\in(3/4,1)\cup[0,1/4)$ in finite time. For $\mathcal{A}=2$ and
$R=1$ we have $\theta_{l}=5/12$, $\theta_{r}=1/12$. The initial condition in
red has angle $\theta^{\ast}=5/6$, for which our results imply a jump angle of
$\theta_{e}(\theta^{\ast})={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{11/12}}$.

######  Lemma 5.4.

Consider the reduced problem (79). Orbits are contained within constant level
sets $H(z,\theta)=K\geq-\mathcal{A}/2$, where

  1. 1.

$K=-\mathcal{A}/2$ corresponds to the folded center $p_{c}$.

  2. 2.

$K\in(-\mathcal{A}/2,\mathcal{A}/2)$ correspond to orbits that connect to the
folded cycle $S_{0}^{\textup{c}}$ in both forward and backward time.

  3. 3.

$K=\mathcal{A}/2$ corresponds to two saddle homoclinic orbits in the
desingularised reduced problem (80), and a limit cycle $\gamma_{c}$ of period
$\tau=2$ along the true and faux canards through the neutral folded saddle
$p_{s}$ in the reduced problem (79).

  4. 4.

$K>\mathcal{A}/2$ corresponds to limit cycles of period $\tau=1$ outside of
the region bounded by $\gamma_{c}$.

###### Proof.

This follows from direct calculations with the desingularised reduced problem
(80) and a reversal of time on $S_{0}^{\textup{r}}$. ∎

We can use the information obtained above to describe the local dynamics near
$S_{0}^{\textup{c}}$. Specifically, we want to describe the transition map
$\pi^{(\alpha)}:\Delta^{\textup{in}}\to\Delta^{\textup{out}}$ induced by the
forward flow of system (76), where

$\Delta^{\textup{in}}\coloneqq\left\\{(z,\theta,R^{2}):|{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{z}}-R|\leq\beta,\theta\in\mathbb{R}/\mathbb{Z}\right\\},\qquad\Delta^{\textup{out}}\coloneqq\left\\{(-R,\theta,a):\theta\in\mathbb{R}/\mathbb{Z},a\in[-a_{0},a_{0}]\right\\},$

for small but fixed $R,\beta,a_{0}>0$. A singular transition map
$\pi_{0}^{(\alpha)}:\Delta^{\textup{in}}\to\Delta^{\textup{out}}$ for
$\varepsilon=0$ can be constructed explicitly, by concatenating solutions to
the layer and reduced problems described above. If we choose $R>0$ so that

$R^{2}<2\mathcal{A},$ (81)

then $\Delta^{\textup{in}}\cap S_{0}^{\textup{a}}$ intersects the singular
canard solution $\gamma_{c}$ at the two points $(R,\theta_{l})$ and
$(R,\theta_{r})$, where $\theta_{l}\in(1/4,3/4]$ and
$\theta_{r}\in[3/4,0)\cup[0,1/4)$ solve

$\sin(2\pi\theta_{l,r})=1-\frac{R^{2}}{\mathcal{A}},$

see Figure 21. Let $\mathcal{I}\subset\mathbb{R}/\mathbb{Z}$ denote the
interior of the interval between $\theta_{l}$ and $\theta_{r}$ which is
contained within the region bounded by $\gamma_{c}$. Different types of
singular orbits can be constructed depending on whether or not one starts with
an initial angle $\theta\in\mathcal{I}$,
$\theta\in\\{\theta_{l},\theta_{r}\\}$ or
$\theta\notin\overline{\mathcal{I}}$.

Consider an initial condition
$(z^{\ast},\theta^{\ast},R^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{)\in\Delta^{\textup{in}}}}$
with $\theta^{\ast}\in\mathcal{I}$. This is mapped by the layer flow of system
(77) to the point $(R,\theta^{\ast},R^{2})\in\Delta^{\textup{in}}\cap
S_{0}^{\textup{a}}$. Since $\theta^{\ast}\in\mathcal{I}$, the point
$(R,\theta^{\ast},R^{2})$ lies on an orbit of the desingularized reduced
problem (80) which is contained within a level curve
$H(z,\theta)=K\in(-\mathcal{A}/2,\mathcal{A}/2)$. In fact, one can show
directly that

$K=K(\theta^{\ast})\coloneqq\frac{\mathcal{A}}{2}\sin(2\pi\theta^{\ast})+\frac{R^{2}}{2}.$
(82)

Lemma 5.4 implies that the forward evolution of $(R,\theta^{\ast},R^{2})$
under the reduced problem (79) reaches $S_{0}^{\textup{c}}$ at a regular jump
point $(0,\theta_{e}(\theta^{\ast}),0)$, where
$\theta_{e}(\theta^{\ast})\in(3/4,1)\cup[0,1/4)$. In fact, the value of
$\theta_{e}$ can be caluclated explicitly using

$H(0,\theta_{e})=\frac{\mathcal{A}}{2}\sin(2\pi\theta_{e})=K(\theta^{\ast}),\qquad\theta_{e}(\theta^{\ast})\in(3/4,1)\cup[0,1/4).$
(83)

Note that the containment condition on the right ensures that
$\theta_{e}(\theta^{\ast})$ is uniquely determined, i.e. that the function
$\theta_{e}:\mathcal{I}\to(3/4,1)\cup[0,1/4)$ is well-defined. Since
$(0,\theta_{e}(\theta^{\ast}),0)\in S_{0}^{\textup{c}}$ is a (classical)
regular jump point, we can connect to the layer flow of system (77) and map it
to the point $(-R,\theta_{e}(\theta^{\ast}),0)\in\Delta^{\textup{out}}$. Thus
we have shown that for all initial conditions with
$(z^{\ast},\theta^{\ast},R^{2})$ with $\theta^{\ast}\in\mathcal{I}$, the
singular transition map is

$\pi_{0}^{(\alpha)}(z^{\ast},\theta^{\ast},R^{2})=\left(-R,\theta_{e}(\theta^{\ast}),0\right),$

where $\theta_{e}(\theta^{\ast})$ is uniquely determined by the expressions
(82) and (83).

Now consider an initial condition
$(z^{\ast},\theta^{\ast},R^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{)\in\Delta^{\textup{in}}}}$
with $\theta^{\ast}\notin\mathcal{I}$. In this case the point
$(z^{\ast},\theta^{\ast},R^{2})$ is mapped via the layer flow to
$(R,\theta^{\ast},R^{2})\in S_{0}^{\textup{a}}$, but this point lies on an
orbit of the reduced problem with $H(z,\theta)\geq\mathcal{A}/2$.
Specifically, if $\theta^{\ast}\in\\{\theta_{l},\theta_{r}\\}$, then
$(R,\theta^{\ast},R^{2})\in\gamma_{c}$. The map $\pi_{0}^{(\alpha)}$ is multi-
valued in this case, since every point on $\gamma_{c}\cap S_{0}^{\textup{r}}$
can be concatenated with an orbit segment of the layer problem. In particular,

$\pi_{0}^{(\alpha)}(z^{\ast},\theta_{l,r},R^{2})=\left\\{\left(-R,\theta,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathcal{A}}}\left(1-\sin(2\pi\theta)\right)\right):\theta\in\mathbb{R}/\mathbb{Z}\right\\}.$

If $\theta^{\ast}\notin\overline{\mathcal{I}}$, then $(R,\theta^{\ast},R^{2})$
lies on a periodic orbit of the reduced problem which is bounded above
$\gamma_{c}$. In this case solutions never leave $S_{0}^{\textup{a}}$, and the
map $\pi_{0}^{(\alpha)}$ is not defined.

The following result describes the perturbation of singular orbits with
initial conditions in the sector defined by requiring that
$\theta\in\mathcal{I}$.

######  Proposition 5.5.

Consider system (76) with $\alpha\in\\{1,2\\}$, assume the constant $R$ which
defines $\Delta^{\textup{in}}$ satisfies (81), and let

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\tilde{\Delta}^{\textup{in}}}}\coloneqq\\{(z,\theta,R^{2})\in\Delta^{\textup{in}}:\theta\in\mathcal{J}_{1}\\},$

for a closed interval
$\mathcal{J}_{1}\subset{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathcal{I}}}$.
There exists an $\varepsilon_{0}>0$ such that for all
$\varepsilon\in(0,\varepsilon_{0})$, the restricted map
$\tilde{\pi}^{(\alpha)}\coloneqq\pi^{(\alpha)}|_{\tilde{\Delta}^{\textup{in}}}:\tilde{\Delta}^{\textup{in}}\to{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\Delta^{\textup{out}}}}$
is well-defined and given by

$\tilde{\pi}^{(\alpha)}(z,\theta,R^{2})=\left(-R,\tilde{\pi}^{(\alpha)}_{\theta}(z,\theta,\varepsilon),\tilde{\pi}^{(\alpha)}_{a}(z,\theta,\varepsilon)\right),$

where

$\tilde{\pi}_{\theta}^{(\alpha)}(z,\theta,\varepsilon)=\theta_{e}(\theta)+\mathcal{O}(\varepsilon^{\alpha}\ln\varepsilon),\qquad\tilde{\pi}_{a}^{(\alpha)}(z,\theta,\varepsilon)=\mathcal{O}(\varepsilon^{2\alpha/3}),$

$\tilde{\pi}^{(\alpha)}\to\pi_{0}^{(\alpha)}|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\tilde{\Delta}^{\textup{in}}}}}$
as $\varepsilon\to 0$, and $\theta_{e}(\theta)$ is uniquely determined by
equations (82) and (83). Explicitly,

$\theta_{e}(\theta)=\frac{1}{2\pi}\arcsin\left(\sin(2\pi\theta)+\frac{R^{2}}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathcal{A}}}}\right),\qquad\theta\in(3/4,1)\cup[0,1/4).$
(84)

###### Proof.

Consider the forward evolution of initial conditions
$(z^{\ast},\theta^{\ast},R^{2})\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\tilde{\Delta}^{\textup{in}}}}$
under system (76). By Fenichel theory, such solutions are initially attracted
at an exponential rate to the Fenichel slow manifold
$S_{\varepsilon}^{\textup{a}}$ which perturbs from $S_{0}^{\textup{a}}$, after
which they remain $\mathcal{O}(\varepsilon^{\alpha})$-close to solutions of
the reduced vector field induced by system (78) until reaching a neighbourhood
of the regular jump point $(-R,\theta_{e}(\theta^{\ast}),0)$. The local
dynamics near $(-R,\theta_{e}(\theta^{\ast}),0)$ are governed by the system

$\begin{split}z^{\prime}&=a-z^{2},\\\
\theta^{\prime}&=\varepsilon^{\alpha},\\\
a^{\prime}&=\varepsilon^{\alpha}\left(-2\pi\mathcal{A}\cos(2\pi\theta_{e}(\theta^{\ast}))+\mathcal{O}((\theta-\theta_{e}(\theta^{\ast}))^{2},\varepsilon^{3-\alpha})\right),\end{split}$

where $\cos(2\pi\theta_{e}(\theta^{\ast}))>0$ since
$\theta_{e}(\theta^{\ast})\in(3/4,1)\cup[0,1/4)$. Direct calculations verify
that this is a stationary two time-scale problem with a regular fold curve
passing through the point $(0,\theta^{\ast},0)$. This observation allows for a
derivation of the desired result using a slight adaptation of the well-known
result in [56]. ∎

Proposition 5.5 describes regular jump type dynamics for solutions with
initial conditions in $\tilde{\Delta}^{\textup{in}}$, which can be chosen
arbitrarily close to $\Delta^{\textup{in}}|_{\theta\in\mathcal{I}}$ by
decreasing $\varepsilon_{0}$ if necessary. It is worthy to note that the
leading order approximation for the ‘jump angle’
$\tilde{\pi}^{(\alpha)}_{\theta}(z,\theta,\varepsilon)\sim\theta_{e}(\theta)$
can be calculated explicitly as a function of the initial angle using (84). It
is also worthy to note that for $\alpha\in\\{1,2\\}$ we have
$\tilde{\pi}^{(\alpha)}_{a}(z,\theta,\varepsilon)=\mathcal{O}(\varepsilon^{2\alpha/3})$,
in contrast to the $\mathcal{O}(\varepsilon^{2})$ estimates for the parameter
drift in Theorem 3.2. This is a consequence of the fact that the ‘parameter
drift’ in system (76) in $a$ occurs on the intermediate, rotational time-
scale. This is not the case for systems satisfying Assumptions 1-2.

We may also consider the case in which initial conditions are chosen ‘higher
up’, in an annular section

$\Sigma^{\textup{in}}\coloneqq\left\\{(z,\theta,\rho^{2}):|{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{z}}-R|\leq\beta,\theta\in\mathbb{R}/\mathbb{Z}\right\\}$

where the constant $\rho>0$ satisfies

$\rho^{2}>R^{2},$ (85)

c.f. equation (81). A subset of solutions with initial conditions in
$\Sigma^{\textup{in}}$ are described by the following result. We formulate it
in terms of the transition map
$\Pi^{(\alpha)}:\Sigma^{\textup{in}}\to\Delta^{\textup{out}}$.

######  Proposition 5.6.

Consider system (76) with $\alpha\in\\{1,2\\}$, assume the constant $\rho$
which defines $\Sigma^{\textup{in}}$ satisfies (85), and let

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widetilde{\Sigma}^{\textup{in}}\coloneqq\\{(z,\theta,\rho^{2})\in\Sigma^{\textup{in}}:\theta\in\mathcal{J}_{2}\\}}},$

where
$\mathcal{J}_{2}:=\left\\{\theta\in\mathbb{R}/\mathbb{Z}:T(\theta)\in\mathcal{J}_{1}\right\\}$,
where $\mathcal{J}_{1}\subset\mathcal{I}$ is the interval from Proposition
5.5, and $T(\theta)$ is the minimal positive solution of the transcendental
equation

$\rho^{2}+\mathcal{A}\sin(2\pi\theta)-R^{2}=\varepsilon^{3}\nu
T+\mathcal{A}\sin(2\pi(\theta+\varepsilon^{\alpha}T)).$ (86)

Then there exists an $\varepsilon_{0}>0$ such that for all
$\varepsilon\in(0,\varepsilon_{0})$, the restricted transition map
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widetilde{\Pi}^{(\alpha)}\coloneqq\Pi^{(\alpha)}|_{\widetilde{\Sigma}^{\textup{in}}}:}}$
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widetilde{\Sigma}^{\textup{in}}\rightarrow\Delta^{\textup{out}}}}$
is given by

$\widetilde{\Pi}^{(\alpha)}(z,\theta,\rho^{2})=\left(-R,\theta_{e}(\theta+\varepsilon^{\alpha}T\mod
1)+\mathcal{O}(\varepsilon^{\alpha}\ln\varepsilon),\mathcal{O}(\varepsilon^{2\alpha/3})\right).$
(87)

###### Proof.

The idea is to consider $\widetilde{\Pi}^{(\alpha)}$ as a composition
$\widetilde{\Pi}^{(\alpha)}=\tilde{\pi}^{(\alpha)}\circ\hat{\pi}^{(\alpha)}$,
where $\tilde{\pi}^{(\alpha)}$ is the local transition map described in
Proposition 5.5 above and
$\hat{\pi}^{(\alpha)}:{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widetilde{\Sigma}^{\textup{in}}}}\to\tilde{\Delta}^{\textup{in}}$.
The $(\theta^{\prime},a^{\prime})$ system decouples from system (76), and can
be solved directly. We obtain

$\theta(t)=\theta^{\ast}+\varepsilon^{\alpha}t\mod
1,\qquad{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a(t)=\rho^{2}+\mathcal{A}\left(\sin(2\pi\theta^{\ast}){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{-}}\sin(2\pi(\theta^{\ast}+\varepsilon^{\alpha}t))\right)-\varepsilon^{3}\nu
t}},$

for all $t\in\mathbb{R}$, where $\theta^{\ast}=\theta(0)$. The transition time
$T=T(\theta^{\ast})>0$ is determined implicitly as the first positive solution
of $a(T)=R^{2}$, i.e. of the transcendental equation (86). The estimate for
$z(T)$ follows from Fenichel theory, which implies that the attracting slow
manifold $S_{\varepsilon}^{\textup{a}}$ which perturbs from
$S_{0}^{\textup{a}}$ can be written as a graph
$z=\varphi_{\varepsilon}(\theta,a)=\sqrt{a}+\mathcal{O}(\varepsilon^{\alpha})$
in this regime (normal hyperbolicity is guaranteed by the fact that $a(t)\geq
R^{2}>0$ for all $t\in[0,T]$).

The preceding arguments imply that

$\hat{\pi}^{(\alpha)}(z,\theta,\rho^{2})=\left({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\varphi_{\varepsilon}(\theta(T),R^{2})}}+\mathcal{O}(\mathrm{e}^{-\kappa/\varepsilon^{\alpha}}),\theta(T),R^{2}\right).$

Since the domain of $\hat{\pi}^{(\alpha)}$ is
$\widetilde{\Sigma}^{\textup{in}}$, we have that
$\theta(T)=\theta(T(\theta^{\ast}))\in\mathcal{J}_{1}$, where
$\mathcal{J}_{1}\subset\mathcal{I}$ is the interval from Proposition 5.5. Thus
the composition
$\widetilde{\Pi}^{(\alpha)}=\tilde{\pi}^{(\alpha)}\circ\hat{\pi}^{(\alpha)}$
is well-defined and, using Proposition 5.5, given by (87). ∎

######  Remark 5.7.

Due to the transcendental nature of equation (86), a closed form expression
for the function $T=T(\theta)$ can only be given in terms of special
functions. Nevertheless, for any given $\theta$, $T(\theta)$ be calculated to
arbitrary precision with standard numerical methods (e.g. Newton or bisection
methods).

Propositions 5.5 and 5.6 describe solutions of ‘regular jump type’, but they
do not explain what happens to solutions with initial conditions in
$\Delta^{\textup{in}}\setminus\mathcal{J}_{1}$. This corresponds to a set of
initial conditions that is ‘small’ in these sense that the diameter of the
interval $[0,1)\setminus\overline{\mathcal{I}}$ tends to zero as $R\to 0$, but
for fixed $R>0$ it is still $\mathcal{O}(1)$ with respect to $\varepsilon\to
0$. Moreover, this set contains solutions close to the singular canard
$\gamma_{c}$, which are expected to play an important role in determining the
qualitative dynamics [55]. A detailed description of the dynamics in this case
is left for future work.

######  Remark 5.8.

For $\alpha\geq 3$, the dynamics can be analysed and understood entirely with
established theory for two time-scale systems. The case $\alpha=3$ is non-
trivial, and features a degenerate bifurcation in the desingularized reduced
problem as the folded saddle $p_{s}$ and neutral folded center $p_{c}$ collide
and annihilate at a particular value of the forcing
$\mathcal{A}=\mathcal{A}_{c}>0$. For $\mathcal{A}<\mathcal{A}_{c}$, the folded
cycle $S_{0}^{\textup{c}}$ is regular, i.e. there are no folded singularities.
For $\alpha\geq 4$, $S_{0}^{\textup{c}}$ is always regular. A direct
application of the results in [56] in this case yields an explicit expression
for the transition map:

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\pi^{(\alpha)}}}\left(z,\theta,R^{2}\right)=\left(-R,\theta+\mathcal{O}(\varepsilon^{\alpha}\ln\varepsilon),\mathcal{O}(\varepsilon^{2\alpha/3})\right),\qquad\alpha\geq
4.$

##  6 Summary and outlook

GSPT is an established and powerful tool for the analysis of stationary fast-
slow systems, particularly when combined with the blow-up method. The
corresponding theory for oscillatory fast-slow systems which possess a limit
cycle manifold instead of (or in addition to) a critical manifold is less
developed, despite the fact that oscillatory systems are common in
applications. One reason for this appears to be that the theory for stationary
fast-slow systems, notably Fenichel-Tikhonov theory and the blow-up method,
depend on quasi-equilibrium properties possessed by stationary but not
oscillatory fast-slow systems. The main purpose of this article has been to
show that both Fenichel-Tikhonov theory and blow-up method can be applied to
study the dynamics of a class of multiple time-scale systems with three or
more time-scales, as long as the angular/rotational dynamics are sufficiently
slow relative to the radial dynamics. The class of systems considered, namely
those in the general form (4), contains a ‘semi-oscillatory’ class of systems
that are oscillatory with respect to the partial singular limit
$\varepsilon_{1}>0,\ \varepsilon_{2}\to 0$, but stationary with respect to the
double singular limit $(\varepsilon_{1},\varepsilon_{2})\to(0,0)$. Although
our approach is not applicable in the purely oscillatory case, our analysis
showed that Fenichel-Tikhonov theory and the blow-up method can be applied
directly in the semi-oscillatory case.

More concretely, we focused on the dynamics near a (three time-scale) global
singularity corresponding to a kind of folded cycle bifurcation in the layer
problem obtained in the partial singular limit. In order to do so we derived a
prototypical system for the three time-scale global singularity near the non-
hyperbolic cycle $S_{0}^{\textup{c}}$, recall Proposition 2.4, and studied the
transition map induced by the flow. Our main result is Theorem 3.2, which is
stated for system (22), which provides a rigorous characterisation of the
asymptotic and strong contraction properties of the flow near
$S_{0}^{\textup{c}}$. The asymptotics
$h_{y}^{(\alpha)}(\theta,\varepsilon)=\mathcal{O}(\varepsilon^{2})$ and the
strong contraction property are the same as for the stationary regular fold
point studied in [39, 48, 56], not only in the semi-oscillatory cases
$\alpha\in\\{1,2\\}$ considered herein in detail but for all
$\alpha\in\mathbb{N}_{+}$. The semi-oscillatory cases $\alpha\in\\{1,2\\}$ are
distinguished from the cases $\alpha\geq 3$ by the fact that the leading order
estimates for $h_{y}^{(\alpha)}(\theta,\varepsilon)$ are, in general, non-
constant functions of $\theta$. We derived an explicit expression for the
leading order estimate in terms of the functions $a(\theta)$, $b(\theta)$ and
$c(\theta)$ appearing in system (22) when $\alpha=2$ in particular. In
addition, we provided detailed estimates for the exit angle, recall (25), as
well as estimates for the total number of rotations about the $y$-axis as
solutions traverse the region between $\Delta^{\textup{in}}$ and
$\Delta^{\textup{out}}$, recall Corollary 3.3. Theorem 3.2 is proved using the
blow-up method in Section 4, thereby demonstrating the applicability of
geometric blow-up techniques to study global singularities of limit cycle type
in suitable classes of three time-scale systems. We focused on a detailed
presentation of the proof in the semi-oscillatory cases $\alpha\in\\{1,2\\}$.
The proof for $\alpha\geq 3$ can be obtained either via straightforward
adaptations of the proof for $\alpha=2$, or via straightforward adaptations of
the proof of the established results for two time-scale systems in [39, 48,
56]. The primary complication in cases $\alpha\in\\{1,2\\}$ stems from the
fact that the dynamics remain ‘global’ in $\theta$ after blow-up near
$S_{0}^{\textup{c}}$. Consequently, one has to consider a blow-up of the
entire cycle (as opposed to a local blow-up about a point on
$S_{0}^{\textup{c}}$). We also demonstrated the applicability of our results
by using them to derive detailed geometric and asymptotic information about
the passage through a folded cycle singularity in a class of periodically
forced Liénard equations; recall Section 5.1 and in particular Theorem 5.2.
Finally, in Section 5.2, we considered a toy model for the study of tipping
phenomena. Due to the periodic forcing in the fast equation, this model
violates both of our primary Assumptions 1-2. The partial geometric analysis
presented herein provides reason to argue that cases of this kind can for the
most part be analysed using classical two time-scale theory.

In summary, the current manuscript provides a precedent for the effectiveness
of stationary methods including Fenichel-Tikhonov theory and the blow-up
method in semi-oscillatory multiple time-scale systems with at least three
time-scales. There are many interesting directions which one might take from
here. Our analysis of the toy model for tipping considered in Section 5.2 is
incomplete, and motivates the analysis of a more general class of systems with
regular folded cycle manifolds which violate our main Assumptions 1-2. One
might also consider the applicability of the methods developed herein for
other global cycle singularities for which there is no direct and lower
dimensional analogue in stationary fast-slow systems, e.g. near a non-
hyperbolic flip/period-doubling type cycle in three time-scale systems of the
form (4). These and other related problems are left for future work.

### Acknowledgements

SJ, CK and SVK acknowledge funding from the SFB/TRR 109 Discretization and
Geometry in Dynamics, i.e. Deutsche Forschungsgemeinschaft (DFG - German
Research Foundation) - Project-ID 195170736 - TRR109. CK thanks the
VolkswagenStiftung for support via a Lichtenberg Professorship. SVK also
acknowledges funding from the DFG priority program SPP 2298 Theoretical
Foundations of Deep Learning as well as the Munich Data Science Institute. SJ,
CK and SVK would also like to thank an anonymous referee for valuable feedback
and suggestions, which led to significant improvements in the article.

### Data Availability

No datasets were generated or analyzed during the current study.

### Conflict of Interest

The authors have no competing interests.

## References

  * [1] O. D. Anosova, On invariant manifolds in singularly perturbed systems, Journal of Dynamical and Control Systems, 5 (1999), pp. 501–507.
  * [2] , Invariant manifolds in singularly perturbed systems, Trudy Matematicheskogo Instituta Imeni VA Steklova, 236 (2002), pp. 27–32.
  * [3] C. Baesens, Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation, Physica D: Nonlinear Phenomena, 53 (1991), pp. 319–375.
  * [4] , Gevrey series and dynamic bifurcations for analytic slow-fast mappings, Nonlinearity, 8 (1995), p. 179.
  * [5] E. Baspinar, D. Avitabile, and M. Desroches, Canonical models for torus canards in elliptic bursters, Chaos: An Interdisciplinary Journal of Nonlinear Science, 31 (2021), p. 063129.
  * [6] G. N. Benes, A. M. Barry, T. J. Kaper, M. A. Kramer, and J. Burke, An elementary model of torus canards, Chaos: An Interdisciplinary Journal of Nonlinear Science, 21 (2011), p. 023131.
  * [7] R. Bertram, M. J. Butte, T. Kiemel, and A. Sherman, Topological and phenomenological classification of bursting oscillations, Bulletin of Mathematical Biology, 57 (1995), pp. 413–439.
  * [8] J. Burke, M. Desroches, A. Granados, T. J. Kaper, M. Krupa, and T. Vo, From canards of folded singularities to torus canards in a forced van der Pol equation, Journal of Nonlinear Science, 26 (2016), pp. 405–451.
  * [9] P. Cardin, P. da Silva, and M. Teixeira, Three time scale singular perturbation problems and nonsmooth dynamical systems, Quarterly of Applied Mathematics, 72 (2014), pp. 673–687.
  * [10] P. T. Cardin and M. A. Teixeira, Fenichel theory for multiple time scale singular perturbation problems, SIAM Journal on Applied Dynamical Systems, 16 (2017), pp. 1425–1452.
  * [11] M. L. Cartwright and J. E. Littlewood, On non-linear differential equations of the second order II, Annals of Mathematics, (1947), pp. 472–492.
  * [12] P. De Maesschalck, E. Kutafina, and N. Popović, Three time-scales in an extended Bonhoeffer–van der Pol oscillator, Journal of Dynamics and Differential Equations, 26 (2014), pp. 955–987.
  * [13] M. Desroches, J. Burke, T. J. Kaper, and M. A. Kramer, Canards of mixed type in a neural burster, Physical Review E, 85 (2012), p. 021920.
  * [14] M. Desroches and V. Kirk, Spike-adding in a canonical three-time-scale model: Superslow explosion and folded-saddle canards, SIAM Journal on Applied Dynamical Systems, 17 (2018), pp. 1989–2017.
  * [15] F. Dumortier and R. Roussarie, Canard cycles and center manifolds, no. 577 in Memoirs of the American Mathematical Society, American Mathematical Society, 1996.
  * [16] B. Ermentrout and D. H. Terman, Mathematical foundations of neuroscience, vol. 35, Springer, 2010.
  * [17] N. Fenichel, Geometric singular perturbation theory for ordinary differential equations, Journal of Differential Equations, 31 (1979), pp. 53–98.
  * [18] 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.
  * [19] J. Grasman, Asymptotic methods for relaxation oscillations and applications, Springer-Verlag, 1987.
  * [20] J. Guckenheimer, Towards a global theory of singularly perturbed dynamical systems, in Nonlinear Dynamical Systems and Chaos, Springer, 1996, pp. 213–225.
  * [21] J. Guckenheimer, K. Hoffman, and W. Weckesser, The forced forced van der Pol equation i: the slow flow and its bifurcations, SIAM Journal on Applied Dynamical Systems, 2 (2003), pp. 1–35.
  * [22] J. Guckenheimer, J. H. Tien, and A. R. Willms, Bifurcations in the fast dynamics of neurons: implications for bursting, in Bursting: The Genesis of Rhythm in the Nervous System, World Scientific, 2005, pp. 89–122.
  * [23] Y. Ilyashenko, Embedding theorems for local maps, slow-fast systems and bifurcation from Morse-Smale to Smale-Williams, Translations of the American Mathematical Society-Series 2, 180 (1997), pp. 127–140.
  * [24] J. Jalics, M. Krupa, and H. G. Rotstein, Mixed-mode oscillations in a three time-scale system of ODEs motivated by a neuronal model, Dynamical Systems, 25 (2010), pp. 445–482.
  * [25] H. Jardón-Kojakhmetov and C. Kuehn, A survey on the blow-up method for fast-slow systems, arXiv preprint arXiv:1901.01402, (2019).
  * [26] S. Jelbart and C. Kuehn, Discrete geometric singular perturbation theory, arXiv preprint arXiv:2201.06996, (2022).
  * [27] C. K. Jones, Geometric singular perturbation theory, in Dynamical systems, vol. 1609 of Lecture Notes in Mathematics, Springer, 1995, pp. 44–118.
  * [28] P. Kaklamanos and N. Popović, Complex oscillatory dynamics in a three-timescale El Niño southern oscillation model, arXiv preprint arXiv:2207.03230, (2022).
  * [29] P. Kaklamanos, N. Popović, and K. U. Kristiansen, Bifurcations of mixed-mode oscillations in three-timescale systems: An extended prototypical example, Chaos: An Interdisciplinary Journal of Nonlinear Science, 32 (2022), p. 013108.
  * [30] , Geometric singular perturbation analysis of the multiple-timescale Hodgkin-Huxley equations, arXiv preprint arXiv:2203.02948, (2022).
  * [31] S. Y. Kirillov and V. Nekorkin, Dynamic saddle-node bifurcation of the limit cycles in the model of neuronal excitability, Radiophysics and Quantum Electronics, 57 (2015), pp. 837–847.
  * [32] I. Kosiuk and P. Szmolyan, Geometric singular perturbation analysis of an Autocatalator model, Discrete and Continuous Dynamical Systems, 2 (2009), pp. 783–806.
  * [33] , Scaling in singular perturbation problems: Blowing up a relaxation oscillator, SIAM Journal on Applied Dynamical Systems, 10 (2011), pp. 1307–1343.
  * [34] , Geometric analysis of the Goldbeter minimal model for the embryonic cell cycle, Journal of Mathematical Biology, 72 (2016), pp. 1337–1368.
  * [35] M. A. Kramer, R. D. Traub, and N. J. Kopell, New dynamics in cerebellar purkinje cells: Torus canards, Physical Review Letters, 101 (2008), p. 068103.
  * [36] K. U. Kristiansen and P. Szmolyan, Relaxation oscillations in substrate-depletion oscillators close to the nonsmooth limit, Nonlinearity, 34 (2021), p. 1030.
  * [37] 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.
  * [38] M. Krupa, N. Popović, and N. Kopell, Mixed-mode oscillations in three time-scale systems: A prototypical example, SIAM Journal on Applied Dynamical Systems, 7 (2008), pp. 361–420.
  * [39] 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.
  * [40] , Extending slow manifolds near transcritical and pitchfork singularities, Nonlinearity, 14 (2001), p. 1473.
  * [41] , Relaxation oscillation and canard explosion, Journal of Differential Equations, 174 (2001), pp. 312–368.
  * [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, Physica D: Nonlinear Phenomena, 431 (2022), p. 133105.
  * [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] B. Letson, J. E. Rubin, and T. Vo, Analysis of interacting local oscillation mechanisms in three-timescale systems, SIAM Journal on Applied Mathematics, 77 (2017), pp. 1020–1046.
  * [47] I. Lizarraga, B. Rink, and M. Wechselberger, Multiple timescales and the parametrisation method in geometric singular perturbation theory, Nonlinearity, 34 (2021), p. 4163.
  * [48] E. Mishchenko and N. Rozov, Kh, Differential equations with small parameters and relaxation oscillations, Nauka, Moscow. (transl.: Plenum Press, New York, 1980), 1975.
  * [49] P. Nan, Y. Wang, V. Kirk, and J. E. Rubin, Understanding and distinguishing three-time-scale oscillations: Case study in a coupled Morris–Lecar system, SIAM Journal on Applied Dynamical Systems, 14 (2015), pp. 1518–1557.
  * [50] K. Nipp and D. Stoffer, Invariant manifolds in discrete and continuous dynamical systems, vol. 21 of EMS Tracts in Mathematics, European Mathematical Society, 2013.
  * [51] 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.
  * [52] L. S. Pontryagin and L. V. Rodygin, Approximate solution of a system of ordinary differential equations involving a small parameter in the derivatives, in Doklady Akademii Nauk, vol. 131(2), Russian Academy of Sciences, 1960, pp. 255–258.
  * [53] J. Rinzel, A formal classification of bursting mechanisms in excitable systems, in Mathematical Topics in Population Biology, Morphogenesis and Neurosciences, Springer, 1987, pp. 267–281.
  * [54] K.-L. Roberts, J. E. Rubin, and M. Wechselberger, Averaging, folded singularities, and torus canards: Explaining transitions between bursting and spiking in a coupled neuron model, SIAM Journal on Applied Dynamical Systems, 14 (2015), pp. 1808–1844.
  * [55] P. Szmolyan and M. Wechselberger, Canards in $\mathbb{R}^{3}$, Journal of Differential Equations, 177 (2001), pp. 419–453.
  * [56] , Relaxation oscillations in $\mathbb{R}^{3}$, Journal of Differential Equations, 200 (2004), pp. 69–104.
  * [57] A. N. Tikhonov, Systems of differential equations containing small parameters in the derivatives, Matematicheskii Sbornik, 73 (1952), pp. 575–586.
  * [58] B. van der Pol, A theory of the amplitude of free and forced triode vibrations, Radio Review, 1 (1920), pp. 701–710.
  * [59] , On “relaxation-oscillations”, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2 (1926), pp. 978–992.
  * [60] T. Vo, Generic torus canards, Physica D: Nonlinear Phenomena, 356 (2017), pp. 37–64.
  * [61] M. Wechselberger, Geometric singular perturbation theory beyond the standard form, in Frontiers in Applied Dynamical Systems: Reviews and Tutorials, Springer International Publishing, 2020.
  * [62] S. Wiggins, Normally hyperbolic invariant manifolds in dynamical systems, vol. 105, Springer Science & Business Media, 1994.
  * [63] J. Zhu, R. Kuske, and T. Erneux, Tipping points near a delayed saddle node bifurcation with periodic forcing, SIAM journal on applied dynamical systems, 14 (2015), pp. 2030–2068.