Discretized Fast-Slow Systems near Pitchfork Singularities

TL;DR

This paper analyzes how discretized fast-slow systems behave near pitchfork singularities, showing the extension of the slow manifold and estimating transition contraction rates using a discrete blow-up method.

Contribution

It adapts the blow-up technique to discrete systems and provides precise estimates for the dynamics near singularities, advancing understanding of numerical discretizations of fast-slow systems.

Findings

Slow manifold extends beyond the singularity in discretized systems.

Transition mapping exhibits quantifiable contraction rates.

The discrete blow-up method effectively analyzes singular behavior.

Abstract

Motivated by the normal form of a fast-slow ordinary differential equation exhibiting a pitchfork singularity we consider the discrete-time dynamical system that is obtained by an application of the explicit Euler method. Tracking trajectories in the vicinity of the singularity we show, how the slow manifold extends beyond the singular point and give an estimate on the contraction rate of a transition mapping. The proof relies on the blow-up method suitably adapted to the discrete setting where a key technical contribution are precise estimates for a cubic map in the central rescaling chart.

Equations

$$P:\left(\begin{array}{c}x\\ y\end{array}\right)\mapsto \left(\begin{array}{c}\bar{x}\\ \bar{y}\end{array}\right)=\left(\begin{array}{c}x+h\left(x(y-x^2)+\lambda \epsilon \right)\\ y+h\epsilon \end{array}\right)$$

$$\displaystyle\begin{array}[]{rcrcl}\frac{\mathrm{d}x}{\mathrm{d}\tau}&=&\dot{x}&=&x(y-x^{2})+\lambda\varepsilon\,,\\ \frac{\mathrm{d}y}{\mathrm{d}\tau}&=&\dot{y}&=&\varepsilon\,,\end{array}$$

$$\displaystyle\begin{array}[]{rcrcl}\frac{\mathrm{d}x}{\mathrm{d}t}&=&x^{\prime}&=&f(x,y,\varepsilon)\,,\\ \frac{\mathrm{d}y}{\mathrm{d}t}&=&y^{\prime}&=&\varepsilon g(x,y,\varepsilon)\,,\quad\ x\in\mathbb{R}^{m},\quad y\in\mathbb{R}^{n},\quad 0<\varepsilon\ll 1\,,\end{array}$$

$$\displaystyle\begin{array}[]{rcrcl}\varepsilon\frac{\mathrm{d}x}{\mathrm{d}\tau}&=&\varepsilon\dot{x}&=&f(x,y,\varepsilon)\,,\\ \frac{\mathrm{d}y}{\mathrm{d}\tau}&=&\dot{y}&=&g(x,y,\varepsilon)\,.\end{array}$$

$$\begin{array}{rl}{x}^{\prime }& =f(x,y,0),\\ y^{\prime }& =0,\end{array}$$

$$\begin{array}{rl}0& =f(x,y,0),\\ y^{\prime }& =g(x,y,0),\end{array}$$

$$S_0:=\{(x,y)\in {\mathbb{R}}^m\times {\mathbb{R}}^n\text{leavevmode}|\text{leavevmode}\text{leavevmode}f(x,y,0)=0\},$$

$$\begin{array}[]{rr}f(0,0,0)=0,&\qquad\partial_{x}f(0,0,0)=0,\\ \partial_{y}f(0,0,0)=0,&\qquad\partial_{xx}f=0,\\ \partial_{xxx}f(0,0,0)\neq 0,&\qquad\partial_{xy}f(0,0,0)\neq 0.\end{array}$$

$${\partial }_{xxx}f(0,0,0)<0,{\partial }_{xy}f(0,0,0)>0\text{leavevmode},$$

$$\begin{array}{rl}\dot{x}& =x(y-x^2)+\lambda \epsilon +h_1(x,y,\epsilon ),\\ \dot{y}& =\epsilon (1+h_2(x,y,\epsilon )),\end{array}$$

$$\begin{array}{rl}\dot{x}& =x(y-x^2)+\lambda \epsilon ,\\ \dot{y}& =\epsilon .\end{array}$$

$$S_0=\{(x,y)\in {\mathbb{R}}^2:y=x^2\}\cup \{(x,y)\in {\mathbb{R}}^2:x=0\}\text{leavevmode}.$$

$$P:\left(\begin{array}{c}x\\ y\end{array}\right)\mapsto \left(\begin{array}{c}\bar{x}\\ \bar{y}\end{array}\right)=\left(\begin{array}{c}x+h\left(x(y-x^2)+\lambda \epsilon \right)\\ y+h\epsilon \end{array}\right).$$

$$S_0=\{(x,y)\in {\mathbb{R}}^2:y=x^2\}\cup \{(x,y)\in {\mathbb{R}}^2:x=0\},$$

$$S_a^0=\{x=0,y<0\},S_a^{+}=\{y=x^2,x>0\},$$

$$S_a^{-}=\{y=x^2,x<0\},S_r^0=\{x=0,y>0\}.$$

$$|1+hy|\ne 1.$$

$${\Delta }_{\text{in}}=\{(x,y)\in {\mathbb{R}}^2\text{leavevmode}|\text{leavevmode}y=-{\rho }^2,x\in J\}\text{leavevmode}.$$

$${\Delta }_{\text{out}}^0=\{(x,y)\in {\mathbb{R}}^2\text{leavevmode}|\text{leavevmode}y={\tilde{\rho }}^2,x\in J\}\text{leavevmode},\text{leavevmode}\text{leavevmode}{\Delta }_{\text{out}}^{\pm }=\{(x,y)\in {\mathbb{R}}^2\text{leavevmode}|\text{leavevmode}y={\tilde{\rho }}^2,x\mp \tilde{\rho }\in J\}$$

$$P:\text{leavevmode}\text{leavevmode}{\mathbb{R}}^4\to {\mathbb{R}}^4\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\left(\begin{array}{c}x\\ y\\ \epsilon \\ h\end{array}\right)\mapsto \left(\begin{array}{c}\bar{x}\\ \bar{y}\\ \bar{\epsilon }\\ \bar{h}\end{array}\right)=\left(\begin{array}{c}x+h\left(x(y-x^2)+\lambda \epsilon \right)\\ y+h\epsilon \\ \epsilon \\ h\end{array}\right),$$

$$\Phi :B:={\mathbb{S}}^2\times (0,\infty )\times (0,\infty )\to {\mathbb{R}}^4,\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}(\tilde{x},\tilde{y},\tilde{\epsilon },\tilde{h},r)\mapsto (x,y,\epsilon ,h),$$

$$x=r\tilde{x},y=r^2\tilde{y},\epsilon =r^4\tilde{\epsilon },h=r^{-2}\tilde{h}.$$

$$\displaystyle K_{1}:D_{1}:=\mathbb{R}\times\mathbb{R}^{+}\times\mathbb{R}^{+}_{0}\times\mathbb{R}^{+}\to\leavevmode\nobreak\ \mathbb{R}\times\mathbb{R}^{-}\times\mathbb{R}^{+}_{0}\times\mathbb{R}^{+}\leavevmode\nobreak\$$

$$\displaystyle K_{2}:D_{2}:=\mathbb{R}\times\mathbb{R}\times\mathbb{R}^{+}\times\mathbb{R}^{+}\leavevmode\nobreak\ \leavevmode\nobreak\ \to\leavevmode\nobreak\ \mathbb{R}\times\mathbb{R}\times\mathbb{R}^{+}\times\mathbb{R}^{+}\leavevmode\nobreak\$$

$$\displaystyle K_{3}:D_{3}:=\mathbb{R}\times\mathbb{R}^{+}\times\mathbb{R}^{+}_{0}\times\mathbb{R}^{+}\to\leavevmode\nobreak\ \mathbb{R}\times\mathbb{R}^{+}\times\mathbb{R}^{+}_{0}\times\mathbb{R}^{+}\leavevmode\nobreak\$$

$$K_1:x$$

$$K_2:x$$

$$K_3:x$$

$${\kappa }_{12}:x_2$$

$${\kappa }_{21}:x_1$$

$${\kappa }_{23}:x_3$$

$${\kappa }_{32}:x_2$$

$$\lambda >0.$$

$${\Sigma }_1^{\text{in}}=\left\{\left(\frac{x}{\rho },\rho ,\delta ,\nu \right):x\in J\right\},$$

Full text

Keywords: pitchfork bifurcation, slow manifolds, invariant manifolds, loss of normal hyperbolicity, blow-up method, discretization.

Mathematics Subject Classification (2010): 34E15, 37M99, 37G10, 34C45, 39A99.

Acknowledgments

This work was supported by the German Research Foundation (DFG) via the SFB-TRR 109. Luca Arcidiacono acknowledges support from the graduate program TopMath of the Elite Network of Bavaria and the TopMath Graduate Center of TUM Graduate School at Technische Universität München.

1 Introduction

We study the dynamical system generated by the two-dimensional cubic polynomial map

[TABLE]

for small $\varepsilon,h>0$. The parameter $\varepsilon$ implies a time scale separation between the fast variable $x$ and the slow variable $y$, while the parameter $h$ is viewed as the stepsize in an explicit Euler method for the ordinary differential equation (ODE)

[TABLE]

which is the normal form of a fast-slow system exhibiting a pitchfork singularity at the origin $(x,y)=(0,0)$. Indeed, for $\varepsilon=0$, we can view $y$ as a bifurcation parameter for the flow in the $x$-variable: when $y<0$, equation (1.4) has a hyperbolic sink at $(0,y)$ and, when $y>0$, there are three equilibria with $(0,y)$ unstable and the other two, $(-\sqrt{y},y)$ and $(\sqrt{y},y)$, locally asymptotically stable. The origin $(x,y)=(0,0)$ is called singular since hyperbolicity is lost at this point, and this is also the case for the map (1.1). We will analyze the dynamics close to the origin for small $\varepsilon,h>0$. Since we focus on the local behaviour around the singularity, we will neglect potential higher order nonlinearities in the majority of this work, but will show how to adapt the proof when including those.

For our analysis we will make use of the blow-up method [1, 2], which has turned out to be a successful tool for treating singular points of fast-slow systems. It was first applied to fast-slow systems by Dumortier and Roussarie [3] to gain insight in the dynamics around non-hyperbolic equilibria. The method uses a non-injective transformation that maps a higher dimensional object such as a sphere onto the non-hyperbolic equilibrium constituting the singularity. The dynamics on this larger, blown up version of the singularity may then be desingularized by an appropriate rescaling of time and exhibit (partially) hyperbolic behaviour. Then one can use dynamical systems techniques to analyze the dynamics in blown-up space. Finally, a typical result allows one to extend invariant manifolds past the singular point in the blown down system; see e.g. [13, Chapter 7] for an introduction and [16, 15, 17, 7, 10, 12, 14] for an, of course non-exhaustive, list of different applications to planar fast-slow systems.

By means of the blow-up method Krupa and Szmolyan [10, 11] analyze different kinds of singularities in fast-slow ODEs, i.e. fold, canard, transcritical and pitchfork singularities, and show how certain invariant manifolds, so-called slow manifolds, extend around the singular points for small $\varepsilon>0$. In the case of fold points, Nipp and Stoffer [18] transform the blow-up technique to the corresponding explicit Runge-Kutta, in particular Euler, discretization and prove the extension of slow manifolds for the discrete time system around the singularity. Whereas they apply an abstract existence theory for invariant manifolds developed in [18], Engel and Kuehn [4] use direct estimates in the blow-up charts to prove the extension of slow manifolds for transcritical singularities. In both cases, a crucial aspect of the discretized blow-up lies in finding the right rescaling of the step size $h$.

In a similar spirit, we investigate, how trajectories of (1.1) behave near the origin and show how the slow manifold may be continued beyond the pitchfork singularity in the discrete setting. We prove that, depending on the sign of $\lambda$, trajectories starting in the vicinity of the single slow branch near $\{(0,y)\in\mathbb{R}^{2}\,:\,y<0\}$ are attracted exponentially to one of the parabolic branches near $\{(\pm\sqrt{y},y)\in\mathbb{R}^{2}\,:\,y>0\}$. Furthermore, for $\lambda=0$, we show that canard-type orbits can track the unstable branch $\{(0,y)\in\mathbb{R}^{2}\,:\,y>0\}$. Our analysis uses three charts that cover different parts of the blown up space around the singularity. We track trajectories through several checkpoints along a curve of fixed points of the cubic map and give estimates on the contraction of the transition mappings. In this way, we also give an alternative way of proof to the result in [11] for the ODE case, by letting the step size $h\to 0$ in (1.1).

This paper is structured as follows. After giving a short introduction to continuous time fast-slow systems and pitchfork singularities in Section 2, we formulate the setup and main results of this paper at the start of Section 3. The major part of Section 3 is dedicated to the proof of the main theorem which is divided into several steps. We start with $\lambda\neq 0$. The relevant coordinate changes are discussed in Section 3.2. Sections 3.3 and 3.5 describe the dynamics in the vicinity of the branches of the critical manifold, which allows to define the slow manifolds and control contractivity of the transition map. In Section 3.4 we describe the continuation of a slow manifold through the blown-up singularity by direct estimates on the trajectories. Finally Section 3.6 combines the results in a blown down version, which finishes the proof for $\lambda\neq 0$. The required modifications to cover the canard case $\lambda=0$ are outlined in Section 3.7. Section 3.8 shows how the previous results can be adapted to a more general setting.

2 Continuous-time fast-slow systems with pitchfork singularity

2.1 Fast-slow ODEs

Fast slow systems occur in various fields of science such as neurobiology or chemistry and are usually found as a system of ODEs with two time scales, this means, they are of the form

[TABLE]

where $f,g,$ are $C^{k}$-functions with $k\geq 3$. The small parameter $\varepsilon$ consitutes the separation between two time scales. The variables $x$ and $y$ are often called the fast variable(s) and the slow variable(s) respectively. The time variable in (2.3), denoted by $t$, is termed the fast time scale. By a change of variables, one can also consider the slow time scale $\tau=\varepsilon t$ and rewrite (2.3) as

[TABLE]

The singular limit $\varepsilon=0$ can be seen from two different perspectives corresponding with the two time scales. Setting $\varepsilon=0$ in (2.3) yields

[TABLE]

which is called the layer problem (or fast subsystem), since we can view the equation layer-wise parametrized by the constant $y$. Setting $\varepsilon=0$ in (2.6) gives the differential algebraic equations

[TABLE]

called the reduced problem (or slow subsystem). The flow of (2.8), the so-called slow flow, is restricted to the set

[TABLE]

which consists of equilibria of the layer problem (2.7). We refer to this set as the critical set or often also critical manifold, in case it is a manifold.

A subset $S\subset S_{0}$ is called normally hyperbolic if the matrix $\mathrm{D}_{x}f(x,y,0)\in\mathbb{R}^{m\times m}$ has no eigenvalue with vanishing real part for all $(x,y)\in S$. In the vicinity of normally hyperbolic submanifolds of $S_{0}$, the dynamics can be very well described for sufficiently small $\varepsilon>0$: Fenichel Theory [6, 9, 13, 19] gives the existence of a locally invariant manifold, the slow manifold $S_{\varepsilon}$, which lies close to $S_{0}$ and maintains the stability properties of the layer problem (2.7). Furthermore, the restriction of (2.6) to $S_{\varepsilon}$ is a regular perturbation of the reduced problem (2.8).

However, points $p\in\mathbb{R}^{m}\times\mathbb{R}^{n}$, which do not satisfy normal hyperbolicity are called singularities in this context and are more delicate to handle. From the view point of the layer equation (2.7) singularities often correspond to bifurcations of the fast subsystem, and the breakdown of normal hyperbolicity is typically associated with the intersection of multiple parts of $S_{0}$ at the point $p$ where degeneracy of $\mathrm{D}_{x}f(p)$ follows from the Implicit Function Theorem. In the case of a pitchfork singularity, which we will consider in the following, we are precisely in such a situation.

2.2 Pitchfork singularity in continuous time

We consider a two-dimensional fast-slow system of the form (2.3) where the critical manifold resembles the shape of a pitchfork, and we call the associated non-hyperbolic singularity a pitchfork singularity. Such a situation occurs when the vector field $f$ satisfies

[TABLE]

These conditions guarantee that, for $\varepsilon=0$, there is a non-hyperbolic equilibrium at the origin where the critical manifold has a transversal self-intersection and one part of the branches crosses the other tangentially to the $x$-direction. In particular, we assume that

[TABLE]

such that the singularity is supercritical. Furthermore, we assume $g(0,0,0)>0$ such that the slow dynamics pass through the origin in positive $y$-direction. In other words, we consider the problem of how the slow dynamics behave in the vicinity of a splitting into three critical branches (see [11, Figure 4]). The case of a subcritical pitchfork singularity or the situation of $g(0,0,0)<0$ are less challenging, since the dynamics only heads into the direction of one critical branch, and will therefore not be treated in this paper.

There is a linear change of coordinates (see [11]) which brings the system into the normal form

[TABLE]

where $h_{1}$ and $h_{2}$ satisfy $h_{1}(x,y,\varepsilon)=\mathcal{O}\left(x^{2}y,xy^{2},\varepsilon x,\varepsilon y,\varepsilon^{2}\right),\leavevmode\nobreak\ h_{2}(x,y,\varepsilon)=\mathcal{O}\left(x,y,\varepsilon\right)$. Since we are mainly interested in the local dynamics around the origin, we may initially ignore the higher order terms and only consider the system

[TABLE]

The critical manifold is given as

[TABLE]

For negative $y$ there is only one stable equilibrium at $x=0$, while for positive $y$ we have an unstable equilibrium at $x=0$ and two locally asymptotically stable ones at $x=\pm\sqrt{y}$.

3 Pitchfork singularity in discrete time

A time-discretization of equation (2.10) by the explicit Euler method with time step size $h>0$ yields the map

[TABLE]

As in continuous time, the system induced by (3.1) clearly possesses the critical manifold

[TABLE]

consisting of fixed points of (3.1) for $\varepsilon=0$. We split the set $S_{0}$ into the four branches

[TABLE]

By linearization we see that these four branches are normally hyperbolic as long as for $(x,y)\in S_{0}\setminus\{(0,0)\}$ we have

[TABLE]

Since we want to restrict the analysis to the non-hyperbolic singularity at the origin $(x,y)=(0,0)$, we always assume that $h$ is chosen small enough so that (3.3) holds as well as the same stability properties as in the time-continuous case. For example, for a fixed initial condition with $y_{0}<0$, we have to ensure $1+hy_{0}>-1$ which yields the restriction $h<2/|y_{0}|$, which then implies that $S_{a}^{0}$ is normally hyperbolic and locally attracting. Note that we shall still use the notation in (3.2) in this context. In contrast to the continuous case the fixed points on $S_{a}^{0}$ are not globally stable, but only inside an interval around zero of size $\sqrt{\frac{2}{h}+y}$. Outside this interval solutions diverge in growing oscillations. Compare also with the reasoning of Lemma 3.5 in Section 3.3. Due to normal hyperbolicity and according to [8, Theorem 4.1], for $\varepsilon$ sufficiently small, there exist corresponding invariant slow manifolds $S_{a,\varepsilon,h}^{0}$ $S_{a,\varepsilon,h}^{+}$ $S_{a,\varepsilon,h}^{-}$ $S_{r,\varepsilon,h}^{0}$. However, the four branches of the critical manifold $S_{0}$ intersect at the origin $(x,y)=(0,0)$, where we have $\mathrm{D}_{x}P(0,0)=1$, i.e. we observe the loss of normal hyperbolicity as in the ODE case.

3.1 Main result

We want to investigate where points around $S_{a,\varepsilon,h}^{0}$ get mapped to by iterations of $P$ in order to find the continuation of $S_{a,\varepsilon,h}^{0}$ beyond the singularity. For that purpose, fix some $\rho>0$, let $J\subset\mathbb{R}$ be a small interval containing [math] and define the section $\Delta_{\textnormal{in}}$ around the point $(0,-\rho^{2})$ on the critical branch $S_{a}^{0}$ by

[TABLE]

In particular, we always assume that the initial condition is chosen on $\Delta_{\textnormal{in}}$ and $J$ is sufficiently small so that trajectories effectively start close to the attracting slow manifold $S_{a,\varepsilon,h}^{0}$.

We are going to follow trajectories of (1.1) starting in $\Delta_{\textnormal{in}}$ up to height $y=\rho^{2}$. Since the line $\{y=\rho^{2}\}$ can only be reached in case $\frac{2\rho^{2}}{\varepsilon h}\in\mathbb{N}$, we introduce $\tilde{\rho}$ as the closest reachable height, which then satisfies $|\rho^{2}-\tilde{\rho}^{2}|<\varepsilon h$. This allows us to define the sections

[TABLE]

around the points $(\tilde{\rho},\tilde{\rho}^{2}),(-\tilde{\rho},\tilde{\rho}^{2})$ or $(0,\tilde{\rho}^{2})$ on the branches $S_{a}^{+},S_{a}^{-}$ or $S_{r}^{0}$ respectively. Depending on the sign of $\lambda$ we will show that the transition maps $\Pi^{\pm}:\leavevmode\nobreak\ \Delta_{\textnormal{in}}\to\Delta_{\textnormal{out}}^{\pm}$ or $\Pi^{0}:\Delta_{\textnormal{in}}\to\Delta_{\textnormal{out}}^{0}$ if $\lambda=0$ are well defined. The transition maps are given by the $N$-fold of $P$ where $N=[\frac{2\rho^{2}}{\varepsilon h}]$ is the closest integer to $\frac{2\rho^{2}}{\varepsilon h}$. For the discrete setting, induced by the map (3.1), we have the following main result (see Figure 1 for an illustration of the case $\lambda>0$).

Theorem 3.1**.**

Fix $\rho>0$ and let $\lambda\in\mathbb{R}$ . Then there are $\varepsilon_{0},h_{0}>0$, depending on $\lambda$, such that for all $\varepsilon\in(0,\varepsilon_{0}]$ and all $h\in(0,h_{0}]$ the following holds.

1. (T1)

If $\lambda>0$, the set $\Delta_{\textnormal{in}}$ (including the point $\Delta_{\textnormal{in}}\cap S_{a,\varepsilon,h}^{0}$) is mapped by $\Pi^{+}$ into a subset of $\Delta_{\textnormal{out}}^{+}$ that contains the point $\Delta_{\textnormal{out}}^{+}\cap S_{a,\varepsilon,h}^{+}$ and has a width of order $\mathcal{O}\left(\left(1-c\cdot h\right)^{\frac{K}{\varepsilon h}}\right)$ for some constants $c,K>0$. 2. (T2)

If $\lambda<0$, an analogous statement holds with $\Delta_{\textnormal{out}}^{-}$ , $S_{a,\varepsilon,h}^{-}$ and $\Pi^{-}$ instead of $\Delta_{\textnormal{out}}^{+}$, $S_{a,\varepsilon,h}^{+}$ and $\Pi^{+}$. 3. (T3)

If $\lambda=0$ the slow manifolds $S_{a,\varepsilon,h}^{0}$ and $S_{r,\varepsilon,h}^{0}$ coincide with the critical branches $S_{a}^{0}$ and $S_{r}^{0}$ and are connected by a canard solution, i.e. $\Pi^{0}(\Delta_{\textnormal{in}}\cap S_{a,\varepsilon,h}^{0})\in S_{r,\varepsilon,h}^{0}$. The set $\Delta_{\textnormal{in}}$ gets mapped by $\Pi^{0}$ into $\Delta_{\textnormal{out}}^{0}$ and its image has a width less than \Big{(}1-h^{2}(\tfrac{\rho^{2}}{2})^{2}\Big{)}^{\tfrac{\rho^{2}}{2\varepsilon h}}\cdot|J|.

Remark 3.2**.**

An analogous result for the continuous-time system (2.10) has been shown in [11]. The sections are chosen to be $\Delta_{\textnormal{out}}^{\pm}=\{(x,y)\in\mathbb{R}^{2}\leavevmode\nobreak\ |\leavevmode\nobreak\ x=\rho^{2},y\mp\rho^{4}\in J\}$ such that the transition maps $\Pi^{\pm}:\leavevmode\nobreak\ \Delta_{\textnormal{in}}\to\Delta_{\textnormal{out}}^{\pm}$ are well defined if $\textnormal{sign}(\lambda)=\pm 1$. The image $\Pi^{\pm}(\Delta_{\textnormal{in}})\subset\Delta_{\textnormal{out}}^{\pm}$ is an interval about the corresponding point on the slow manifold $S_{a,\varepsilon}^{\pm}$ and has a size of $\mathcal{O}\left(\mathrm{e}^{-\frac{C}{\varepsilon}}\right)$ for some constant $C>0$.

Note that we will assume $h_{0}\varepsilon_{0}\ll\rho$ to obtain meaningful time lengths for the dynamical analysis. For further details on the choice of $h_{0}$ see Section 3.6, but notice that it immediately implies the stability restriction $h<\frac{2}{\rho^{2}}$ from the discussion below (3.3).

3.2 Transformation to the chart coordinates

The proof of Theorem 3.1 uses the blow-up method for the dynamical system induced by the map

[TABLE]

where the fast-slow separation parameter $\varepsilon$ and the stepsize $h$ are also seen as variables.

The quasi-homogeneous blow-up transformation around the pitchfork singularity is given as

[TABLE]

The transformation of the $(x,y,\varepsilon)$-coordinates is the same as in the continuous-time case (see [11]). The change of variables in $h$ is chosen such that the map is desingularized in the relevant charts. We exclude [math] from the domain of $\tilde{h}$ since at $\tilde{h}=0$ every point is a neutral fixed point. Due to the transformation $h=\tilde{h}/r$ we have to exclude [math] from the domain of $r$ as well.

The transformation $\Phi$ induces a map $\bar{P}:=\Phi^{-1}\circ P\circ\Phi$ on the manifold $B$. We analyse the dynamics of $\bar{P}$ by using the charts $K_{i}$, $i=1,2,3$,

[TABLE]

which are given by

[TABLE]

To switch between different chart coordinates we use the following coordinate changes

[TABLE]

and

[TABLE]

For the proof of Theorem 3.1 we will proceed as follows. Transforming (3.4) using the coordinate changes $K_{i}\leavevmode\nobreak\ (i=1,2,3)$ induces a dynamical system on $D_{i}\subset\mathbb{R}^{4}$ in chart coordinates. Trajectories of (3.4) are analyzed via their corresponding transformed versions in each of the charts. In every chart $K_{i}$ we will define sets $\Sigma_{i}^{\textnormal{in}}$ and $\Sigma_{i}^{\textnormal{out}}$, show that the transition maps $\Pi_{i}:\Sigma_{i}^{\textnormal{in}}\to\Sigma_{i}^{\textnormal{out}}$ are well defined and study their contractivity. The mappings $\Pi^{\pm}$ are then built by connecting the three chart-wise transition maps, which are elaborated in Section 3.6 in more detail.

We refer to $K_{1}$ as the entering chart, as we start our analysis in this chart and trajectories are brought closer to the origin. Charts of the type of $K_{2}$ are often called rescaling charts, since the transformation is basically a rescaling with suitable powers of the fast-slow-separation constant $\varepsilon$. In this chart, the dynamics arbitrarily close to the origin are analyzed. Finally the exiting chart $K_{3}$ is used to describe, how trajectories exit the vicinity of the origin and is crucial for the contractivity statement of Theorem 3.1.

3.3 Dynamics in the entering chart

Fix some $\rho>0$ and also consider the case

[TABLE]

from now on until Section 3.7. The case $\lambda<0$ can be treated analogously, see also Section 3.6 for more details. Further take $\varepsilon,h>0$ sufficiently small. During the next sections we will specify what sufficiently small means for $\varepsilon$ and $h$ such that $\varepsilon_{0}$ and $h_{0}$ are determined. In the coordinates $(x_{1},r_{1},\varepsilon_{1},h_{1})$ of the first chart $K_{1}$, the set $\Delta_{\textnormal{in}}$ is given as

[TABLE]

for which we define

[TABLE]

We investigate the dynamics within the domain

[TABLE]

In order to find an expression for the map (3.4) in terms of the entering chart $K_{1}$, we first rewrite $\bar{y}=y+\varepsilon h$ in $K_{1}$-coordinates as

[TABLE]

This yields

[TABLE]

The remaining three equations of (3.4) in $K_{1}$-coordinates read as

[TABLE]

Hence, by using (3.5), we obtain the maps

[TABLE]

The dynamics in $r_{1},\varepsilon_{1}$ and $h_{1}$ can be calculated explicitly for the first chart.

Lemma 3.3**.**

For $\xi_{0}:=\frac{1}{h_{1}(0)\varepsilon_{1}(0)}>0$, the trajectories of (3.10) in $r_{1},\varepsilon_{1},h_{1}$ are given by

[TABLE]

for $n\in\mathbb{N},n<\xi_{0}.$

Proof.

We define $\eta:=h_{1}\varepsilon_{1}$. Multiplying the last two equations of (3.10) gives the relation

[TABLE]

Given the initial condition $\eta(0)=\xi_{0}^{-1}>0$, we obtain the solution $\eta(n)=(\xi_{0}-n)^{-1}$ . Hence, we can compute

[TABLE]

Similarly, we obtain the formulas for $r_{1}$ and $\varepsilon_{1}$. ∎

Furthermore, we can observe from equations (3.10) that the set $\{\varepsilon_{1}=0\}$ is invariant for the dynamics and, for given $r_{1}^{*},h_{1}^{*}$, consists of the two-parameter family of invariant one-dimensional lines

[TABLE]

Each of these lines has a fixed point located at $(x_{1},r_{1},\varepsilon_{1},h_{1})=(0,r_{1}^{*},0,h_{1}^{*})$, which has a three-dimensional center eigenspace and a one-dimensional stable eigenspace in $x_{1}$-direction with eigenvalue $1-h_{1}^{*}$ (recall that $h_{1}^{*}\leq h\rho^{2}<\frac{1}{3}<2$). In other words, the two-dimensional plane

[TABLE]

is an invariant manifold in $D_{1}$ only consisting of fixed points, attracting in the $x_{1}$-direction and neutral in the other directions, and corresponding to the branch $S_{{\textnormal{a}}}^{0}$ of the critical manifold. In particular, for each $h_{1}\geq 0$ we have the fixed point

[TABLE]

We obtain the following statement:

Proposition 3.4**.**

The invariant manifold $S_{{\textnormal{a}},1}^{0}$ extends to a center-stable invariant manifold $M_{{\textnormal{a}},1}^{0}$ (at $p_{a}^{0}(0)$) which is given in $D_{1}$ by a graph $x_{1}=l_{1}(r_{1},\varepsilon_{1},h_{1})$ for some mapping $l_{1}$. Furthermore, $\varepsilon_{1}$ is increasing in $D_{1}$ (and thereby, in particular, on $M_{{\textnormal{a}},1}^{0}$), whereas $h_{1},r_{1}$ are decreasing in $D_{1}$ (and thereby also on $M_{{\textnormal{a}},1}^{0}$).

Proof.

This is an immediate consequence of the considerations above and classical center manifold theory (cf. e.g. [8, Chapter 5A]). From Lemma 3.3, we can see immediately that as long as $r_{1}(n),h_{1}(n)>0$, we have $n<\xi_{0}$. Hence, the claim follows from the formulas in Lemma 3.3. ∎

Note that, on $\{r_{1}>0,\varepsilon_{1}>0,h_{1}>0\}$, the manifold $M_{a,1}^{0}$ corresponds to the union of the slow manifolds $S_{a,\varepsilon,h}^{0}$. Assume we iterate system (3.10) until $r_{1}(n)$ reaches a value less or equal to $\frac{\rho}{2}$. (This means that the trajectory of (3.4) has reached a point above the line $y=-\rho^{2}/4$ in $x,y$-coordinates.) As already remarked before, when defining $\Delta_{\textnormal{out}}^{\pm}$, the level $\frac{\rho}{2}$ might not be hit precisely. However, let us assume for simplicity that we are in the situation of reaching $r_{1}=\frac{\rho}{2}$ after $n^{*}\in\mathbb{N}$ iterates. (Also in the next sections we will assume in a similar way that specific values are hit since the small errors do clearly not change our results.) This means we have $\frac{\xi_{0}-n^{*}}{\xi_{0}}=\frac{1}{4}$, and reversely, $n^{*}=\frac{3}{4}\xi_{0}=\frac{3}{4\nu\delta}$.

Denote by $\Pi_{1}:\Sigma_{1}^{\textnormal{in}}\to\mathbb{R}^{4}$ the transition map after $\frac{3}{4\nu\delta}$ iterations. We can deduce the following Lemma:

Lemma 3.5**.**

For $\delta>0$ sufficiently small, we have

[TABLE]

and $M_{{\textnormal{a}},1}^{0}\cap\Sigma_{1}^{\textnormal{in}}$ as well as $M_{{\textnormal{a}},1}^{0}\cap\Sigma_{1}^{\textnormal{out}}$ are non-empty sets.

Proof.

From the explicit solutions in Lemma 3.3, we obtain

[TABLE]

For $\varepsilon=0$ the $x_{1}$-equation of (3.10) reads

[TABLE]

Consequently we have

[TABLE]

Additionally a direct computation yields that for $|x_{1}|\leq\sqrt{\frac{2}{h}-2}$ we have

[TABLE]

This means that $\bar{x}_{1}$ lies in the cone bordered by $\pm(1-h_{1})x_{1}$ for all $x_{1}$ that satisfy $|x_{1}|\leq\sqrt{\frac{2}{h_{1}}-2}$ (see also Figure 2). Thus contraction towards the fixed point at 0 is guaranteed for all initial values of $x_{1}$ in that domain and the speed is at least the linear rate $1-h_{1}$.

Note that $h_{1}(n)$ stays inside the interval $\left[\frac{\nu}{4},\nu\right]$ for $n\leq n^{*}$ and the stable eigenvalue in $x_{1}$-direction is given by $1-h_{1}$ along the manifold $S_{{\textnormal{a}},1}^{0}$. We may choose $\nu$ small enough so that $\rho^{-1}J\subset\left[-\sqrt{\frac{2}{\nu}-2},\sqrt{\frac{2}{\nu}-2}\leavevmode\nobreak\ \right]$. Hence it follows from standard perturbation arguments that, for $\delta$ sufficiently small, the map $\Pi_{1}$ is a contraction with rate $(1-c)^{\frac{3}{4\nu\delta}}$ for some constant $c<\frac{\nu}{4}$. Therefore we have $\Pi_{1}(\Sigma_{1}^{\textnormal{in}})\subset\Sigma_{1}^{\textnormal{out}}$ for a sufficiently small choice of $\delta$.

Since for $\delta\to 0$ the point $M_{{\textnormal{a}},1}^{0}\cap\{\varepsilon_{1}=\delta\}$ approaches $(0,\rho,0,\nu)$, we can deduce that $M_{{\textnormal{a}},1}^{0}\cap\{\varepsilon_{1}=\delta\}$ lies in $\Sigma_{1}^{\textnormal{in}}$ for $\delta$ small enough. Additionally, since the manifold $M_{{\textnormal{a}},1}^{0}$ is invariant under the forward iterations of (3.10), we have that $M_{{\textnormal{a}},1}^{0}\cap\Sigma_{1}^{\textnormal{out}}$ is non-empty as well. ∎

3.4 Dynamics in the rescaling chart

We use $\kappa_{12}$ to transfer the set $\Sigma_{1}^{\textnormal{out}}$ to the second chart and define

[TABLE]

For the transformation of (3.4) via the chart $K_{2}$, first observe that since $\bar{\varepsilon}=\varepsilon$ we have

[TABLE]

Using the coordinates of $K_{2}$, the remaining equations from (3.4) become

[TABLE]

which can be simplified with (3.12) so that we get in total

[TABLE]

Since $r_{2}$ and $h_{2}$ stay constant in this chart, we can plug in the values from $\Sigma_{2}^{\textnormal{in}}$ to write the maps as

[TABLE]

We denote this two-dimensional map by

[TABLE]

Furthermore, we abbreviate

[TABLE]

such that we have $f(x_{2},y_{2})=x_{2}+\nu\delta^{\frac{1}{2}}\left(x_{2}(y_{2}-x_{2}^{2})+\tilde{\lambda}\right)$. We will see that the appropriate exiting set in chart $K_{2}$ is given by

[TABLE]

Recall from chart $K_{1}$ the attracting center manifold $M_{{\textnormal{a}},1}^{0}$ (see Proposition 3.4). Similarly to [4], this manifold corresponds with the global manifold $M_{{\textnormal{a}}}^{0}$ on the blow-up manifold $B$. In chart $K_{2}$ we therefore have the attracting center manifold $M_{{\textnormal{a}},2}^{0}=\kappa_{12}\left(M_{{\textnormal{a}},1}^{0}\right)$ (with, reversely, $M_{{\textnormal{a}},1}^{0}=\kappa_{21}\left(M_{{\textnormal{a}},2}^{0}\right)$), whose behaviour is described in the following main result of this section (recall that we are still in the setting $\lambda>0$).

Theorem 3.6**.**

The set

[TABLE]

gets mapped by iteration of $F$ into the set

[TABLE]

where $\mu=\sqrt{2\lambda\rho}\delta^{\frac{1}{4}}$. In particular, the transition mapping  $\Pi_{2}:\Sigma_{2}^{\textnormal{in}}\to\Sigma_{2}^{\textnormal{out}}$   is well defined and maps $M_{{\textnormal{a}},2}^{0}\cap\Sigma_{2}^{\textnormal{in}}$ to $M_{{\textnormal{a}},2}^{0}\cap\Sigma_{2}^{\textnormal{out}}$.

For the proof of this theorem we will analyze the evolution of our starting set at different heights (see Figure 3), treated in the propositions below. Before we do that, we collect some properties of the function $f$ in the following Lemmas 3.7–3.9.

Firstly, we characterize the positive fixed points of $f(\cdot,y_{2})$ on $y_{2}$-fibres.

Lemma 3.7**.**

For any $y_{2}\in\mathbb{R}$, the mapping $x_{2}\mapsto f(x_{2},y_{2})$ has precisely one positive fixed point $x_{2}^{*}(y_{2})$, satisying the equation

[TABLE]

As long as $y_{2}\in[-\frac{1}{4}\delta^{-\frac{1}{2}},\frac{1}{4}\delta^{-\frac{1}{2}}]$, the family of fixed points $x_{2}^{*}(y_{2})$ is monotonically increasing and satisfies $x_{2}^{*}(y_{2})<\frac{1}{2}\delta^{-\frac{1}{4}}+1$ .

Proof.

Fixed points of $f(\cdot,y_{2})$ are characterized by the equation

[TABLE]

As $x_{2}=0$ is not a solution, we can solve this for $y_{2}$ and obtain

[TABLE]

The graph of $x_{2}\mapsto x_{2}^{2}-\frac{\tilde{\lambda}}{x_{2}}$ consists of two branches. Notice that $x_{2}^{2}-\frac{\tilde{\lambda}}{x_{2}}>0$ for $x_{2}<0$. For $x_{2}>0$ the graph is monotonically increasing (since the summands $x_{2}^{2}$ and $-\frac{\tilde{\lambda}}{x_{2}}$ are). Thus we have a curve of positive fixed points, that we can also parametrize by $y_{2}$, call them $x_{2}^{*}(y_{2})$. It remains to verify the bound of $x_{2}^{*}(y_{2})$. Using $\tilde{\lambda}<1$, we have

[TABLE]

Consequently the point $\left(\frac{1}{2}\delta^{-\frac{1}{4}}+1,\frac{1}{4}\delta^{-\frac{1}{2}}\right)$ lies below the graph of $x_{2}\mapsto x_{2}^{2}-\frac{\tilde{\lambda}}{x_{2}}$, so right of the curve $\{\left(x_{2}^{*}(y_{2}),y_{2}\right)\}$ of fixed points. Since the curve is increasing, all other fixed points $x_{2}^{*}(y_{2})$ with $y_{2}\leq\frac{1}{4}\delta^{-\frac{1}{2}}$ also satisfy the estimate $x_{2}^{*}(y_{2})\leq\frac{1}{2}\delta^{{}^{-}\frac{1}{4}}+1$. ∎

Secondly, we find $f(\cdot,y_{2})$ to be monotonically increasing on a suitable interval and, by using Lemma 3.7, we find invariant sets under $f(\cdot,y_{2})$.

Lemma 3.8**.**

Let $y_{2}\in[-\frac{1}{4}\delta^{-\frac{1}{2}},\frac{1}{4}\delta^{-\frac{1}{2}}]$. Then the following holds:

The function $f(\cdot,y_{2})$ is monotonically increasing on $[-\frac{1}{2}\delta^{-\frac{1}{4}}-1,\frac{1}{2}\delta^{-\frac{1}{4}}+1]$ for all $\nu$ sufficiently small. 2. 2.

The set $[0,\frac{1}{2}\delta^{-\frac{1}{4}}+1]$ is (positively) invariant under $f(\cdot,y_{2})$. Furthermore, when $y_{2}\leq 0$, the set $[-\frac{1}{2}\delta^{-\frac{1}{4}},\frac{1}{2}\delta^{-\frac{1}{4}}]$ is (positively) invariant under $f(\cdot,y_{2})$.

Proof.

We compute the derivative

[TABLE]

which shows that the cubic function $f(\cdot,y_{2})$ may have two stationary points located at $x_{2}=\pm\sqrt{\frac{1}{3}(\nu^{-1}\delta^{-\frac{1}{2}}+y_{2})}$ and then is monotonically increasing in between these.

By choosing $\nu$ sufficiently small, we can achieve that the stationary points exist for all $y_{2}\in[-\frac{1}{4}\delta^{-\frac{1}{2}},\frac{1}{4}\delta^{-\frac{1}{2}}]$ and that $\sqrt{\frac{1}{3}(\nu^{-1}\delta^{-\frac{1}{2}}+y_{2})}\geq\frac{1}{2}\delta^{-\frac{1}{4}}+1$. The choice of $\nu$ can be made independently of $\delta$, for all arbitrarily small $\delta$. E. g. for $\delta\leq 1$ the inequality is fulfilled for all $\nu\leq\frac{2}{15}$. Hence, the first claim follows.

Since we have precisely one positive fixed point by Lemma 3.7, the graph of the continuous function $f(\cdot,y_{2})$ hits the diagonal precisely once. As $f(0,y_{2})=\nu\delta^{\frac{3}{4}}\lambda\rho>0$ and $\lim_{x_{2}\to\infty}f(x_{2},y_{2})=-\infty$ we can conclude that it crosses the diagonal at $x_{2}^{*}$ and $f(\cdot,y_{2})$ lies above the diagonal (i.e. $f(x_{2},y_{2})>x_{2}$) for $x_{2}\in[0,x_{2}^{*}(y_{2}))$ and below the diagonal (i.e. $f(x_{2},y_{2})<x_{2}$) for $x_{2}\in(x_{2}^{*}(y_{2}),\infty)$. If $y_{2}\leq 0$ we can widen the estimate $f(x_{2},y_{2})<x_{2}$ to $x_{2}\in\left(-\infty,x_{2}^{*}(y_{2})\right)$, because there is overall only one fixed point.

We have already shown that $f(\cdot,y_{2})$ is increasing on $[-\frac{1}{2}\delta^{-\frac{1}{4}}-1,\frac{1}{2}\delta^{-\frac{1}{4}}+1]$. To prove the claimed invariance of $[0,\frac{1}{2}\delta^{-\frac{1}{4}}+1]$ and $[-\frac{1}{2}\delta^{-\frac{1}{4}},\frac{1}{2}\delta^{-\frac{1}{4}}]$, it thus suffices to consider only the endpoints of the intervals.

Recall from Lemma 3.7 that $x_{2}^{*}(y_{2})\leq\frac{1}{2}\delta^{{}^{-}\frac{1}{4}}+1$ and, hence, we can conclude with the above that $f(\frac{1}{2}\delta^{-\frac{1}{4}}+1,y_{2})<\frac{1}{2}\delta^{-\frac{1}{4}}+1$. Since additionally $f(0,y_{2})>0$, the claimed positive invariance of $[0,\frac{1}{2}\delta^{-\frac{1}{4}}+1]$ follows.

Regarding the invariance of $[-\frac{1}{2}\delta^{-\frac{1}{4}},\frac{1}{2}\delta^{-\frac{1}{4}}]$ in case of $y_{2}\leq 0$, first note that the unique fixed point $x_{2}^{*}(0)$ is given by $x_{2}^{*}(0)=\tilde{\lambda}^{\frac{1}{3}}=\left(\delta^{\frac{1}{4}}\lambda\rho\right)^{\frac{1}{3}}$, which can be directly checked with the defining equation (3.14). As the curve of fixed points is increasing, we can estimate

[TABLE]

for sufficiently small $\delta$ and all $y_{2}\leq 0$. With the considerations above, we thus have

[TABLE]

yielding the invariance of the interval $[-\frac{1}{2}\delta^{-\frac{1}{4}},\frac{1}{2}\delta^{-\frac{1}{4}}]$. ∎

Finally, we can show the following contraction property of $f$ on relevant intervals.

Lemma 3.9**.**

Let $y_{2}\in[-\frac{1}{4}\delta^{-\frac{1}{2}},-\frac{1}{8}\delta^{-\frac{1}{2}}]$. The function $f(\cdot,y_{2})$ restricted to the interval $[-\frac{1}{2}\delta^{-\frac{1}{4}},\frac{1}{2}\delta^{-\frac{1}{4}}]$ is a contraction with constant $1-\tfrac{\nu}{8}$.

Proof.

Recall from the proof of Lemma 3.8 that the derivative reads

[TABLE]

As stated in Lemma 3.8, the map $f(\cdot,y_{2})$ is monotonically increasing on $[-\frac{1}{2}\delta^{-\frac{1}{4}},\frac{1}{2}\delta^{-\frac{1}{4}}]$ for $y_{2}\in[-\frac{1}{4}\delta^{-\frac{1}{2}},-\frac{1}{8}\delta^{-\frac{1}{2}}]$, and, hence, the derivative $\frac{\partial}{\partial x_{2}}f(x_{2},y_{2})$ is non-negative for all $(x_{2},y_{2})\in E:=[-\frac{1}{2}\delta^{-\frac{1}{4}},\frac{1}{2}\delta^{-\frac{1}{4}}]\times[-\frac{1}{4}\delta^{-\frac{1}{2}},-\frac{1}{8}\delta^{-\frac{1}{2}}]$. Thus we obtain

[TABLE]

Hence, the claim follows by a standard application of the mean-value theorem. ∎

We now turn to showing the transitions from $I_{i}$ to $I_{i+1}$ for $i=1,2,3$ (see Figure 3). Each transition is formulated in one of the following Propositions.

Proposition 3.10**.**

The set $I_{1}=\left[-\frac{1}{2}\delta^{-\frac{1}{4}},\frac{1}{2}\delta^{-\frac{1}{4}}\right]\times\{-\frac{1}{4}\delta^{-\frac{1}{2}}\}$ gets mapped by iterations of (3.13) into the set $I_{2}:=\left[2\delta^{\frac{3}{4}}\lambda\rho-\Delta(\delta),8\delta^{\frac{3}{4}}\lambda\rho+\Delta(\delta)\right]\times\{-\frac{1}{8}\delta^{-\frac{1}{2}}\}$, where $\Delta(\delta):=\frac{1}{2}\delta^{-\frac{1}{4}}\mathrm{e}^{-\frac{1}{64}\delta^{-1}}$.

Proof.

It takes $\frac{1}{8}\delta^{-\frac{1}{2}}\cdot(\nu\delta^{\frac{1}{2}})^{-1}=\frac{1}{8}\nu^{-1}\delta^{-1}$ iterations to get from $y_{2}=-\frac{1}{4}\delta^{-\frac{1}{2}}$ to $y_{2}=-\frac{1}{8}\delta^{-\frac{1}{2}}$ with steps of size $\nu\delta^{\frac{1}{2}}$. Since $[-\frac{1}{2}\delta^{-\frac{1}{4}},\frac{1}{2}\delta^{-\frac{1}{4}}]$ is invariant under $f(\cdot,y_{2})$ for all $y_{2}\in[-\frac{1}{4}\delta^{-\frac{1}{2}},-\frac{1}{8}\delta^{-\frac{1}{2}}]$ by Lemma 3.8 and $f(\cdot,y_{2})$ is a contraction with constant $1-\frac{\nu}{8}$ by Lemma 3.9, the image of $I_{1}$ under iterations of (3.13) has a width of at most $\delta^{-\frac{1}{4}}\left(1-\frac{\nu}{8}\right)^{\frac{1}{8}\nu^{-1}\delta^{-1}}$. Since the exponential limit

[TABLE]

is attained monotonically from below, we can bound this width via

[TABLE]

Hence, any trajectory starting in $I_{1}$ goes through the set

[TABLE]

where $x_{2,\min}^{*}$ and $x_{2,\max}^{*}$ denote the minimal and the maximal value, respectively, that the fixed point $x_{2}^{*}(y_{2})$ attains while $y_{2}$ varies over $[-\frac{1}{4}\delta^{-\frac{1}{2}},-\frac{1}{8}\delta^{-\frac{1}{2}}]$. Since, according to Lemma 3.7, the curve of fixed points is increasing, we have $x_{2,\min}^{*}=x_{2}^{*}(-\frac{1}{4}\delta^{-\frac{1}{2}})$ and $x_{2,\max}^{*}=x_{2}^{*}(-\frac{1}{8}\delta^{-\frac{1}{2}})$.

It remains to verify that $x_{2}^{*}(-\frac{1}{4}\delta^{-\frac{1}{2}})\geq 2\delta^{\frac{3}{4}}\lambda\rho$ and $x_{2}^{*}(-\frac{1}{8}\delta^{-\frac{1}{2}})\leq 8\delta^{\frac{3}{4}}\lambda\rho$. To that purpose, we plug these bounds into equation (3.14), satisfied by the fixed points, which gives

[TABLE]

where the last inequality holds for all sufficiently small $\delta$, and

[TABLE]

Since the curve of fixed points is increasing (see Figure 4), the claim follows. ∎

Proposition 3.11**.**

The set $I_{2}$ gets mapped by iterations of (3.13) into the set $I_{3}:=\left[0,\tilde{\lambda}^{\frac{1}{3}}\right]\times\{0\}$.

Proof.

Recall that the fixed point $x_{2}^{*}(0)$ is given by $x_{2}^{*}(0)=\tilde{\lambda}^{\frac{1}{3}}=\left(\delta^{\frac{1}{4}}\lambda\rho\right)^{\frac{1}{3}}$. Furthermore, we have

[TABLE]

for all sufficiently small $\delta$. Analogously to Lemma 3.9, we observe that, for all $y_{2}\in[-\frac{1}{8}\delta^{-\frac{1}{2}},0)$, the map $f(\cdot,y_{2})$ is a contraction on $\left[0,\tilde{\lambda}^{\frac{1}{3}}\right]$ to the fixed point $x_{2}^{*}(y_{2})$. Since $x_{2}^{*}(y_{2})\in\left[0,\tilde{\lambda}^{\frac{1}{3}}\right]$ for all $y_{2}<0$, we can deduce that $\left[0,\tilde{\lambda}^{\frac{1}{3}}\right]$ is invariant under $f(\cdot,y_{2})$ for all $y_{2}\in\left[-\frac{1}{8}\delta^{-\frac{1}{2}},0\right)$. This implies the claim. ∎

The final transition is described in the following proposition.

Proposition 3.12**.**

The set $I_{3}$ gets mapped by iterations of (3.13) into the set

[TABLE]

where $\mu=\sqrt{2\lambda\rho}\delta^{\frac{1}{4}}$.

Proof.

We first consider the set $S_{1}$ (see Figure 6), defined by

[TABLE]

Let $(x_{2,0},y_{2,0})$ be in $S_{1}$. Now, since $f(\cdot,y_{2,0})$ is increasing and also the curve of fixed points is increasing, we have

[TABLE]

Hence, the point $(\bar{x}_{2,0},\bar{y}_{2,0})$ lies in $S_{1}$, as long as $y_{2,0}\leq\frac{1}{4}\delta^{-\frac{1}{2}}-\nu\delta^{\frac{1}{2}}$. This means that trajectories leave $S_{1}$ through the upper boundary at $y_{2}=\frac{1}{4}\delta^{-\frac{1}{2}}$.

Moreover, we consider the sets

[TABLE]

as depicted in Figure 6. For points $(x_{2},y_{2})\in S_{3}$, we have $\bar{x}_{2}=f(x_{2},y_{2})>x_{2}+\nu\delta^{\frac{3}{4}}\lambda\rho=x_{2}+\nu\delta^{\frac{1}{2}}\tilde{\lambda}$ as well as $\bar{y}_{2}=y_{2}+\nu\delta^{\frac{1}{2}}$. Thus they will be mapped into $S_{2}\cup S_{3}$, as long as $y_{2}\leq\frac{1}{4}\delta^{-\frac{1}{2}}-\nu\delta^{\frac{1}{2}}$. Points on the left-hand boundary of $S_{2}$, given by $\{y_{2}=x_{2}^{2}\}$, clearly satisfy the same estimate such that they are also mapped into $S_{2}\cup S_{3}$. By monotonicity of $f(\cdot,y_{2})$, see Lemma 3.7, this holds for all points in $S_{2}$. Hence, we deduce that trajectories leave $S_{2}\cup S_{3}$ through the upper boundary at $y_{2}=\frac{1}{4}\delta^{-\frac{1}{2}}$.

In order to find a right-hand bound for the set $S_{1}$, we compute

[TABLE]

where we used that $\mu=\sqrt{2\tilde{\lambda}}\leavevmode\nobreak\ \delta^{\frac{1}{8}}=\sqrt{2\lambda\rho}\leavevmode\nobreak\ \delta^{\frac{1}{4}}$. Since the curve of fixed points is increasing (see Figure 4), this means that we can estimate $x_{2}^{*}(\frac{1}{4}\delta^{-\frac{1}{2}})\leq\frac{1}{2}\delta^{-\frac{1}{4}}+\mu$. We observe from the proof of Proposition 3.11 that $\left(0,\tilde{\lambda}^{\frac{1}{3}}\right)\times\{0\}$ is a subset of $S_{2}\cup S_{3}$. Furthermore, it is now easy to see that

[TABLE]

and

[TABLE]

This concludes the proof. ∎

Theorem 3.6 is now an immediate consequence of combining Proposition 3.10, Proposition 3.11 and Proposition 3.12.

3.5 Dynamics in the exiting chart

Transforming $\Sigma_{2}^{\textnormal{out}}$ to the coordinates of the third chart gives

[TABLE]

Furthermore, we define

[TABLE]

with $\theta=\max\{1-\tfrac{\lambda\rho}{2},\tfrac{1}{2}\}$, such that $\kappa_{23}(\Sigma_{2}^{\textnormal{out}})\subset\Sigma_{3}^{\textnormal{in}}$ for sufficiently small $\delta$. As already noted in section 3.3, we may assume, due to the controllably small error, that specific levels are hit by the trajectories and therefore, in particular, we assume $\rho=\tilde{\rho}$, where $\tilde{\rho}$ was used to define $\Delta_{\textnormal{out}}^{+}$ (3.1).

Similarly to the situation in $K_{1}$, we consider the domain

[TABLE]

for the chart $K_{3}$, and we obtain the map

[TABLE]

corresponding with (3.4). Similarly to the system obtained in $K_{1}$, the special structure of (3.20) again allows to explicitly determine solutions of the induced dynamical system in the components $r_{3},\varepsilon_{3}$ and $h_{3}$.

Lemma 3.13**.**

For $\zeta_{0}:=\frac{1}{h_{3}(0)\varepsilon_{3}(0)}>0$, the trajectories of (3.20) in $r_{3},\varepsilon_{3},h_{3}$ are given by

[TABLE]

for $n\in\mathbb{N}$.

Proof.

Let $\vartheta:=h_{3}\varepsilon_{3}$. Multiplying both sides of the equations for $\varepsilon_{3}$ and $h_{3}$ in (3.20) yields

[TABLE]

Solving this recursion for some initial condition $\vartheta(0)=\zeta_{0}^{-1}>0$ gives

[TABLE]

We use this observation to calculate

[TABLE]

and analogously

[TABLE]

This finishes the proof. ∎

We observe from (3.20) that the hyperplane $\{\varepsilon_{3}=0\}$ is an invariant set for the system (3.20), foliating into the invariant lines $\{\varepsilon_{3}=0,r_{3}=r_{3}^{*},h_{3}=h_{3}^{*}\}$ for all $r_{3}^{*},h_{3}^{*}\geq 0$. Each of these lines has three fixed points, located at $(x_{3},r_{3},\varepsilon_{3},h_{3})=(s,r_{3}^{*},0,h_{3}^{*})$, $s\in\{-1,0,1\}$. Linearizing around each of these, for $h_{3}^{*}>0$, we see that in $x_{3}$-direction the fixed point at $x_{3}=0$ is unstable with eigenvalue $1+h_{3}^{*}$ while those at $x_{3}=-1$ and $x_{3}=1$ are stable with eigenvalue $1-2h_{3}^{*}$ (recall that $h_{3}^{*}\leq h\rho^{2}<\frac{1}{3}<1$).

Since, in our considerations, we enter $K_{3}$ via $\Sigma_{3}^{\textnormal{in}}$, our main interest lies in the family of stable fixed points at $(x_{3},r_{3},\varepsilon_{3},h_{3})=(1,r_{3}^{*},0,h_{3}^{*})$, corresponding with the branch $S_{a}^{+}$ of the critical manifold for $r_{3}^{*}>0$. Each of these fixed points has a three-dimensional center eigenspace as well as a one-dimensional stable eigenspace in $x_{3}$-direction with eigenvalue $1-2h_{3}^{*}$. The union of these fixed points forms an invariant manifold, which we call

[TABLE]

In particular, for each $h_{3}\geq 0$ it contains the fixed point

[TABLE]

which has gained hyperbolicity due to the desingularization of the origin. In analogy to Proposition 3.4 we get the following:

Proposition 3.14**.**

The invariant manifold $S_{{\textnormal{a}},3}^{+}$ extends to a center-stable invariant manifold $M_{{\textnormal{a}},3}^{+}$ (at $p_{a}^{+}(0)$) which is given in $D_{3}$ by a graph $x_{3}=l_{3}(r_{3},\varepsilon_{3},h_{3})$ for a smooth mapping $l_{3}$. Furthermore, $\varepsilon_{3}$ is decreasing in $D_{3}$ (and thereby, in particular, on $M_{{\textnormal{a}},3}^{+}$), whereas $h_{3},r_{3}$ are increasing in $D_{3}$ (and thereby on $M_{{\textnormal{a}},3}^{+}$).

Proof.

This is an immediate consequence of the considerations above and classical center manifold theory. Similarly to the proof of Proposition 3.4, the second claim follows from Lemma 3.13. ∎

Note that, on $\{r_{3}>0,\varepsilon_{3}>0,h_{3}>0\}$, the manifold $M_{a,3}^{+}$ corresponds to the union of the slow manifolds $S_{a,\varepsilon,h}^{+}$.

We will follow the iterations of (3.20) until the values $r_{3}=\rho,\varepsilon_{3}=\delta,h_{3}=\nu$ are reached. Lemma 3.13 tells us that this is the case when $\frac{\zeta_{0}+n^{*}}{\zeta_{0}}=4$, which means that the number of iterations equals $n^{*}=3\zeta_{0}=\frac{3}{4\nu\delta}$. Let $\Pi_{3}:\Sigma_{3}^{\textnormal{in}}\to\mathbb{R}^{4}$ be the transition map induced by $n^{*}$ iterations of (3.20) and consider the set

[TABLE]

Note that with this choice we have $K_{3}(\Sigma_{3}^{\textnormal{out}})=\Delta_{\textnormal{out}}^{+}$. The following Lemma summarizes the properties of the transition map $\Pi_{3}$.

Lemma 3.15**.**

For $\delta>0$ sufficiently small, we have $\Pi_{3}(\Sigma_{3}^{\textnormal{in}})\subset\Sigma_{3}^{\textnormal{out}}$, where $\Pi_{3}(\Sigma_{3}^{\textnormal{in}})$ has a $x_{3}$-width of at most $\left(1-C^{*}\frac{\nu}{4}\right)^{\frac{3}{4\nu\delta}}$ for some constant $C^{*}>0$. Furthermore, the intersections $M_{{\textnormal{a}},3}^{+}\cap\Sigma_{3}^{\textnormal{in}}$ as well as $M_{{\textnormal{a}},3}^{+}\cap\Pi_{3}(\Sigma_{3}^{\textnormal{in}})$ are non-empty sets.

Proof.

From the explicit solutions in Lemma 3.13, we directly see that

[TABLE]

Next, we consider the $x_{3}$-equation of (3.20) for $\varepsilon_{3}=0$ and $h_{3}>0$, which reads

[TABLE]

The cubic function $g(\cdot,h_{3})$ is increasing between the stationary points at

[TABLE]

Taking $h_{3}\leq\frac{1}{3}$, we achieve that $\kappa\geq 2>1+\theta$. Since the fixed point at $x_{3}=0$ is unstable we can deduce that the set $(0,\kappa]$ gets attracted to the fixed point at $x_{3}=1$, and so does $\left[1-\theta,1+\theta\right]\subset(0,\kappa]$.

We turn to giving estimates on the contraction rate towards $x_{3}=1$ for different starting values in $[1-\theta,1+\theta].$ Note that the map $g(\cdot,h_{3})$ is increasing and concave on $(0,\kappa]$. Since it is increasing, the subintervals $[1-\theta,1)$ and $(1,1+\theta]$ get mapped by $g(\cdot,h_{3})$ into themselves. Due to its concavity, the function $g(\cdot,h_{3})$ lies below the tangent at $x_{3}=1$, so that the contraction rate for values in $\left(1,1+\theta\right]$ is at least as strong as the linear rate $1-2h_{3}$ coming from the linearization around $x_{3}=1$. Concavity also yields that $g(\cdot,h_{3})$ lies above the secant on the interval $[1-\theta,1]$. Thus the slope of the corresponding secant gives an estimate for the contraction rate of points in $[1-\theta,1)$ to the fixed point at $x_{3}=1$.

In more detail, linear interpolation between the points $\left(1-\theta,g(1-\theta,h_{3})\right)$ and $(1,1)$ yields the slope

[TABLE]

where

[TABLE]

Note that $0<C^{*}(\theta)<2$ and $C^{*}(\theta)\to 0$ as $\theta\to 1$.

We conclude that, on the interval $[1-\theta,1+\theta]$, the map $g(\cdot,h_{3})$ is contracting with constant $1-h_{3}C^{*}$. The transition map $\Pi_{3}$ defined on

[TABLE]

is in $x_{3}$-direction a perturbation of the $n^{*}$-fold of $g(\cdot,h_{3})$, as we now have $r_{3},\varepsilon_{3}>0$.

Also note that during the $n^{*}=\frac{3}{4\nu\delta}$ iterations that define $\Pi_{3}$, the variable $h_{3}$ lies in the interval $[\tfrac{\nu}{4},\nu]$. Thus, if we choose $\delta$ sufficiently small, we can achieve that $\Pi_{3}$ is a contraction (in $x_{3}$-direction) with a rate of $\frac{1}{2}\left(1-C^{*}\frac{\nu}{4}\right)^{\frac{3}{4\nu\delta}}$. Consequently the image $\Pi_{3}(\Sigma_{3}^{\textnormal{in}})$ has a $x_{3}$-width of at most

[TABLE]

Since for $\delta\to 0$ the point $M_{{\textnormal{a}},3}^{+}\cap\left\{r_{3}=\frac{\rho}{2},\varepsilon_{3}=16\delta,h_{3}=\frac{\nu}{4}\right\}$ approaches $(1,\frac{\rho}{2},0,\frac{\nu}{4})$, so especially its $x_{3}$ coordinate approaches $1$, we can deduce that $M_{{\textnormal{a}},3}^{+}\cap\left\{r_{3}=\frac{\rho}{2},\varepsilon_{3}=16\delta,h_{3}=\frac{\nu}{4}\right\}$ lies in $\Sigma_{3}^{\textnormal{in}}$ for $\delta$ small enough. Moreover, since the manifold $M_{{\textnormal{a}},3}^{+}$ is invariant under the forward iterations of (3.20), the point $M_{{\textnormal{a}},3}^{+}\cap\left\{r_{3}=\rho,\varepsilon_{3}=\delta,h_{3}=\nu\right\}$ lies within $\Pi_{3}(\Sigma_{3}^{\textnormal{in}})$.

Because $M_{{\textnormal{a}},3}^{+}\cap\left\{r_{3}=\rho,\varepsilon_{3}=\delta,h_{3}=\nu\right\}$ converges to $(1,\rho,0,\nu)$ as $\delta\to 0$, we infer that for small $\delta$ the point $M_{{\textnormal{a}},3}^{+}\cap\left\{r_{3}=\rho,\varepsilon_{3}=\delta,h_{3}=\nu\right\}$ and the exponentially small set $\Pi_{3}(\Sigma_{3}^{\textnormal{in}})$ around this point lie within $\Sigma_{3}^{\textnormal{out}}=\left\{x_{3}\in 1+\rho^{-1}J,\leavevmode\nobreak\ r_{3}=\rho,\leavevmode\nobreak\ \varepsilon_{3}=\delta,\leavevmode\nobreak\ h_{3}=\nu\right\}$. ∎

3.6 Blown-down dynamics

As a last step, we will connect the individual results, obtained in each of the three charts, to prove Theorem 3.1.

Proof of Theorem 3.1 (T1) and (T2).

Throughout the previous proofs we needed to make $\delta$ sufficiently small; now we choose $\delta_{0}>0$ such that all of the statements hold true for every $\delta\in(0,\delta_{0}]$. This gives the value $\varepsilon_{0}:=\delta_{0}\rho^{4}$.

As a consequence of condition (3.21) in the proof of Lemma 3.15, we assumed $\nu\leq\frac{1}{3}$. Even stronger, Lemma 3.8.1 gives $\nu\leq\frac{2}{15}$. For Lemma 3.5 we required $\nu\leq\left({1+\frac{|J|^{2}}{8\rho^{2}}}\right)^{-1}$. Thus we can take $h_{0}:=\min\left\{\frac{2}{15\rho^{2}},\frac{1}{\rho^{2}+\frac{1}{8}|J|^{2}}\right\}$.

In the previous sections we have seen that

[TABLE]

so that we can concatenate the three transition maps $\Pi_{1},\Pi_{2}$ and $\Pi_{3}$ , using the appropriate coordinate changes in between, and define the map $\tilde{\Pi}$ from $\Sigma_{1}^{\textnormal{in}}$ to $\Sigma_{3}^{\textnormal{out}}$

[TABLE]

In particular, we have seen from the analysis in the charts $K_{1}$, $K_{2}$ and $K_{3}$ (see Lemma 3.5, Theorem 3.6 and Lemma 3.15) that $M_{{\textnormal{a}}}^{0}$, equalling $M_{{\textnormal{a}},1}^{0}$ in $K_{1}$ and $M_{{\textnormal{a}},2}^{0}$ in $K_{2}$, is continued to $\Sigma_{3}^{\textnormal{out}}$ in the vicinity of $M_{{\textnormal{a}},3}^{+}$. Additionally, we have $\Delta_{\textnormal{in}}=K_{1}(\Sigma_{2}^{\textnormal{in}})$ and $\Delta_{\textnormal{out}}^{+}=K_{3}(\Sigma_{3}^{\textnormal{out}})$. Hence, by reverting to the original coordinates $(x,y,\varepsilon,h)$, we obtain the transition map $\Pi^{+}$ from $\Delta_{\textnormal{in}}$ to $\Delta_{\textnormal{out}}^{+}$ given by

[TABLE]

We have shown in Lemma 3.15 that the set $\Pi_{3}(\Sigma_{3}^{\textnormal{in}})\subset\Sigma_{3}^{\textnormal{out}}$ has a $x_{3}$-width of $\left(1-C^{*}\frac{\nu}{4}\right)^{\frac{3}{4\nu\delta}}$. Hence, by transforming back to the blown-down coordinates, it is easy to see that the width in $x$- direction of $K_{3}(\Pi_{3}(\Sigma_{3}^{\textnormal{in}}))\subset\Delta_{\textnormal{out}}^{+}$ is $\mathcal{O}\left(\left(1-C\cdot h\right)^{\frac{K}{\varepsilon h}}\right)$ when we take $K=\frac{3\rho^{2}}{4}$ and $C=C^{*}(\theta)\frac{\rho^{2}}{4}$ (with $\theta$ and $C^{*}$ from (3.15) and (3.22) respectively).

Furthermore, Lemma 3.5 implies that $K_{1}(\Sigma_{1}^{\textnormal{in}})$ contains a point of the slow manifold $S_{a,\varepsilon,h}^{0}$ . On the other hand, Lemma 3.15 shows that the exponentially small set $K_{3}(\Pi_{3}(\Sigma_{3}^{\textnormal{in}}))$ contains a point of the slow manifold $S_{a,\varepsilon,h}^{+}$ , and, due to (3.23) and the definition of $\Pi^{+}$, we also have

[TABLE]

This completes the proof of the statement in Theorem 3.1 for the case $\lambda>0$.

Finally, when $\lambda<0$, observe that under the change of variables $x\mapsto-x$ the $x$-equation in system (3.4) gets transformed into

[TABLE]

which is equivalent to

[TABLE]

Hence, the analysis is the same as for positive $\lambda$, with the same outcome under symmetric change of variables. ∎

3.7 Canard Case

The analysis in the case of $\lambda=0$ may be carried out without a blow-up transformation. Hence, we treat the proof of Theorem 3.1 (T3) separately here.

Proof of Theorem 3.1 (T3).

One observes, that the system

[TABLE]

keeps $\{x=0\}$ invariant for all $\varepsilon\geq 0$, since there is no dependency on $\varepsilon$ in the $x$-equation. This means that we have $S_{a,\varepsilon,h}^{0}=S_{a}^{0}$ and $S_{r,\varepsilon,h}^{0}=S_{r}^{0}$ on the domain of our analysis (see discussion around (3.3)) so that $S_{a,\varepsilon,h}^{0}$ and $S_{r,\varepsilon,h}^{0}$ are connected. The connecting trajectory $\gamma$ starting in $\Delta_{\textnormal{in}}$ is explicitly given by \gamma(n)=\big{(}x(n),y(n)\big{)}=\big{(}0,-\rho^{2}+nh\varepsilon\big{)}.

The linearization along the trajectory $\gamma$ is characterized by the variational equation

[TABLE]

While the fixed point $w=(0,1)^{\top}$ of (3.24) corresponds with the centre-direction along $\gamma$, the solution of (3.24) starting at $v(0)=(1,0)^{\top}$ corresponds with the transversal hyperbolic direction and can be explicitly solved to be

[TABLE]

Let us for simplicity assume that we have $\frac{2\rho^{2}}{\varepsilon h}=:N\in 4\mathbb{N}$. In particular, we set $\Delta_{\textnormal{out}}^{0}=\{(x,\rho^{2})|\leavevmode\nobreak\ x\in J\}$, where clearly $\gamma(N)\in\Delta_{\textnormal{out}}^{0}$. We have already seen in Section 3.3 that due to the cubic structure the contraction rates in $x$-direction towards the locally stable fixed point for $y<0$ are at least as strong as the linear rates achieved by linearization around $x=0$. Similarly one also observes that the linear rates give a bound for the expansion rate for positive values of $y$. Hence, the contraction and expansion of trajectories from $\Delta_{\textnormal{in}}$ to a neighborhood of $\Delta_{\textnormal{out}}^{0}$ can be estimated from above by the linear rate $\mu$ along the trajectory $\gamma$ which satisfies

[TABLE]

Hence, we can give the bound

[TABLE]

meaning that the transition map $\Pi^{0}$ is contractive for the canard case since the contraction rates along $S_{a}^{0}$ prevail over the expansion rates along $S_{r}^{0}$. Hence, the claim follows. ∎

However, note that

[TABLE]

Thus, as expected, in the limit $h\to 0$ one obtains the stability behaviour of the corresponding continuous-time system where contraction and expansion exactly compensate each other. It is still remarkable that the Euler method not only preserves the stability behaviour for trajectories close to the canard but even enhances stability as compared to the continuous-time case for sufficiently small $h>0$. We already observed this surprising effect in the case of transcritical canards (cf. [4]) but emphasize that in other similar situations, like the folded canard (cf. [5]), the Euler method has clearly unfavourable stability properties.

3.8 Higher Order Terms

We briefly discuss how our results can be generalized when higher order terms $h_{1}(x,y,\varepsilon)=\mathcal{O}\left(x^{2}y,xy^{2},\varepsilon x,\varepsilon y,\varepsilon^{2}\right)$ and $h_{2}(x,y,\varepsilon)=\mathcal{O}\left(x,y,\varepsilon\right)$ from (2.9) are included. The corresponding discretized dynamical system reads

[TABLE]

Note that due to the dependence of $\bar{y}$ on $x$, points in the image of $\Delta_{\textnormal{in}}$ under iteration of $P$ will not share the same $y$-coordinate. Thus we cannot define the transition mappings $\Pi^{\pm,0}$ by just a fixed number of iterations of $P$, but instead pointwise for each initial value $(x_{0},-\rho^{2})\in\Delta_{\textnormal{in}}$ by

[TABLE]

We see that for a sufficiently small choice of $\rho$ every trajectory will get $\varepsilon h$ close to $\{y=\rho^{2}\}$. In [4] this had already to be taken into account for the case of a normal form without higher order terms.

Firstly, let us discuss the problem for fixed $\lambda\neq 0$. It is an important benefit of the blow-up method that in entering and exiting charts higher order terms have no significant impact. In more detail, we transform system (3.25) by $K_{1}$, proceeding as in Section 3.3, to obtain

[TABLE]

This relation yields

[TABLE]

which simplifies to

[TABLE]

Consequently the transformed system in $K_{1}$ can be written as

[TABLE]

For $r_{1}=0$, this system is identical to (3.10). Hence, for sufficiently small $r_{1}$, we still obtain the existence of a center-stable manifold $M_{a,1}^{0}$ at the point $p_{a}^{0}(0)$ and the consequences thereof. A small choice of $r_{1}$ means that we have to restrict $\rho$ to sufficiently small values. For the exiting chart $K_{3}$, the situation is similar.

In the rescaling chart, however, the higher order terms may not be bypassed that easily, but the strategy from Section 3.4 can be adapted. As in Section 3.4, we still have $\bar{r}_{2}=r_{2}$ and $\bar{h}_{2}=h_{2}$. The remaining equations of (3.25) transform to

[TABLE]

and

[TABLE]

which can be simplified and desingularized into

[TABLE]

[TABLE]

The following arguments will not only require small $\varepsilon$ but also sufficiently small $\rho$, so that the impact of normal form higher order terms can be controlled and the dynamics are determined by the remaining terms. Since the small parameter $\delta$ incorporates $\rho$ and $\varepsilon$, it is more apparent in the original not blown up coordinates $(x,y,\varepsilon,h)$ how the choice of $\rho$ determines the considered neighbourhood of the origin. Note that in the following small $\rho$ and $\varepsilon$ mean that the statements hold for sufficiently small fixed $\rho$ and for all positive $\varepsilon$ below some sufficiently small threshold. Restricting $\rho$ and $\varepsilon$, we can assure upwarded movement in $y_{2}$-direction taking $\mathcal{O}\left(\nu^{-1}\delta^{-1}\right)$ steps to travel through the domain considered in the second chart.

Our approach in Section 3.4 relied heavily on the curve of fixed points $x_{2}^{*}(y_{2})$ introduced in Lemma 3.7. In the more general setting involving the higher order terms, a curve corresponding to $x_{2}^{*}(y_{2})$ persists. In other words, one can show that for fixed, negative values of $y_{2}$ there is exactly one positive fixed point for equation (3.31), given that $\rho$ and $\varepsilon$ (and thus $\delta$) are sufficiently small. Note that this can only be accomplished for $y_{2}$ outside an interval of size $\mathcal{O}\left(r_{2}^{2}\right)$, since otherwise the term $x_{2}y_{2}$ is of order $\mathcal{O}\left(x_{2}r_{2}^{2}\right)$ and therefore does not dominate additional terms of that order any longer. Using the curve of fixed points a result analogous to Proposition 3.10 can be shown, which implies that trajectories will enter the quadrant $\{x>0,y<0\}$. Equivalent statements to those of Lemma 3.8 can also be obtained for $\rho$ and $\varepsilon$ small enough. This means that for sufficiently small $\nu$ the mapping $x_{2}\mapsto\bar{x}_{2}$ is monotone. Using the monotonicity one then easily checks that $\{x_{2}\geq 0\}$ is invariant under (3.31), simply by plugging in $x_{2}=0$. In a similar manner one can ensure that trajectories leave the rectangular set $\{x_{2}\in[0,\delta^{-\frac{1}{4}}],y_{2}\in[-\frac{1}{4}\delta^{-\frac{1}{2}},\frac{1}{4}\delta^{-\frac{1}{2}}]\}$ only through its upper boundary.

Moreover, one can show that in the quadrant $\{x\geq 0,y\geq 0\}$ trajectories are bounded away from the $y_{2}$-axis, independently from $\varepsilon$. The corresponding result is found in Proposition 3.12. Hence, we may deduce that all trajectories will end up at a $y_{2}$-height close to $\frac{1}{4}\delta^{-\frac{1}{2}}$ with positive $x_{2}$-value, bounded away from 0 and smaller than $\delta^{-\frac{1}{4}}$. In other words, we obtain a result similar to Theorem 3.6. Transforming the exiting set into $K_{3}$-coordinates allows to proceed similarly as in Section 3.5. As we have already seen in chart $K_{1}$, the higher order terms do not change the behaviour for sufficiently small $\rho$. Summarizing, we deduce that Theorem 3.1, (T1) and (T2), can be transferred to the general setting of (3.25).

Furthermore, note that in the general case of the normal form (2.9) including higher order terms, the value of $\lambda$ close to [math] giving a canard changes with the value of $\varepsilon$. For continuous time, this phenomenon is studied in detail for canards in folds in [10] and discussed for the transcritical and pitchfork case in [11, Remark 2.2 and Remark 4.1]. Using a Melnikov computation, one may show the existence of a function $\lambda_{c}(\varepsilon^{1/4})$ with $\lambda_{c}(0)=0$ such that for $\lambda=\lambda_{c}(\varepsilon^{1/4})$ the slow manifold $S_{{\textnormal{a}},\varepsilon}^{-}$ extends to $S_{{\textnormal{r}},\varepsilon}^{+}$ for sufficiently small $\varepsilon$. In Theorem 3.1 (T3), we only treated the case $\lambda_{c}(\varepsilon^{1/4})\equiv 0$ since we did not take into account perturbations from higher order terms. In order to obtain an analogous result to the ODE case, a treatment of the more general problem (3.25) about $\lambda=0$ requires a discrete Melnikov computation, which is more complicated. Therefore, we are going to treat the general canard problem in the separate study [5].