# Two-Stroke Relaxation Oscillators

Samuel Jelbart & Martin Wechselberger
School of Mathematics & Statistics, University of Sydney,
Camperdown, NSW 2006, Australia

###### Abstract

Two-stroke relaxation oscillations consist of two distinct phases per cycle –
one slow and one fast – which distinguishes them from the well-known van der
Pol-type ‘four-stroke’ relaxation oscillations. This type of oscillation can
be found in singular perturbation problems in non-standard form, where the
slow-fast timescale splitting is not necessarily reflected in a slow-fast
variable splitting. The existing literature on such non-standard problems has
developed primarily through applications – we compliment this by providing a
general framework for the application of geometric singular perturbation
theory in this non-standard setting and illustrate its applicability by
proving existence and uniqueness results on a general class of two-stroke
relaxation oscillators. We apply this non-standard geometric singular
perturbation toolbox to a collection of examples arising in the dynamics of
nonlinear transistors and models for mechanical oscillators with friction.

##  1 Introduction

The term relaxation oscillation was coined by Balthasar van der Pol in the
1920s [1, 2] to distinguish nonlinear from harmonic oscillations observed in
electronic circuits. He also derived a prototypical mathematical model, the
now well-known van der Pol (vdP) oscillator model,

$\ddot{x}+r_{vdp}(x)\dot{x}+x=0$ (1.1)

with nonlinear differential ‘resistance’

$r_{vdp}(x)=\mu(x^{2}-1)\,,$ (1.2)

where $\mu\geq 0$ denotes the main system parameter that measures the ratio of
the two characteristic timescales of the electronic circuit model under study.
The overdot denotes differentiation with respect to time $\tau$. This
oscillator model can be recast as a dynamical system (in Liénard form),

$\displaystyle\begin{array}[]{lcl}\dot{\bar{y}}=-x,\\\
\dot{x}=\bar{y}-R_{vdp}(x),\end{array}$ (1.5)

with

$R_{vdp}(x)=\mu\bigg{(}\frac{x^{3}}{3}-x\bigg{)}\,,$ (1.6)

i.e. $R^{\prime}_{vdp}(x)=r_{vdp}(x)$. The function $R_{vdp}(x)$ denotes the
(non-linear) characteristic of the oscillator model which takes the form of a
cubic.

######  Remark 1.1.

The (dimensionless) current-voltage (I-V) characteristic (1.6) derives from
the presence of a tunnel diode (an active, non-linear element). The
requirement that the characteristic has negative slope
$R^{\prime}_{vdp}(x)=r_{vdp}(x)<0$ in parts of the phase space, here for
$|x|<1$, gives rise to an effective negative resistance $r_{vdp}(x)$ that
allows energy to be pumped back (relaxed) into the system which is necessary
for oscillatory behaviour.

We are interested in the relaxation case $\mu\gg 1$ of the vdP oscillator
model (1.5) and define a new variable $y=\bar{y}/\mu$, change to a new (fast)
timescale $dt=\mu d\tau$, and set $\epsilon:=1/\mu^{2}\ll 1$ to obtain

$\displaystyle\begin{array}[]{lcl}x^{\prime}=y+x-\frac{x^{3}}{3},\\\
y^{\prime}=-\epsilon x,\end{array}$ (1.9)

where the dash notation denotes differentiation with respect to the new (fast)
time $t$.

(a) Time traces.

(b) Phase space and limit cycle.

Figure 1.1: Four-stroke relaxation oscillation in the vdP oscillator model
(1.9) with $\epsilon=10^{-2}$. In 1(b) the characteristic (1.6) which
represents the $x$-nullcline of system (1.5) is shown in green, and the unique
equilibrium of system (1.5) is shown in red. The same colour convention
applies without reference in subsequent figures.

Figure 1.1 shows the typical vdP relaxation oscillator time traces of system
(1.9) as well as the associated limit cycle and cubic-shaped characteristic in
phase space. In particular, we notice that the time trace $x(t)$ is
characterised by alternations between two slow and two fast motions over the
course of a single relaxation cycle (‘slow-fast-slow-fast’) while the time
trace $y(t)$ is uniformly slow. This is the hallmark of vdP-type relaxation
oscillations, and the phase space representation shown in Figure 1(b) reveals
the four distinct phases of the limit cycle.

######  Remark 1.2.

All oscillators considered in this work will be recast as dynamical systems in
such a way that the corresponding ‘characteristic’ can be identified with a
nonlinear nullcline. The specific form of the characteristic depends on the
manner in which the problem is recast as a dynamical system. Typically, this
will be clear from the context.

From a dynamical systems point of view, system (1.9) represents a singular
perturbation problem in standard form

$\displaystyle\begin{array}[]{lcl}x^{\prime}=f(x,y,\epsilon),\\\
y^{\prime}=\epsilon g(x,y,\epsilon),\end{array}$ (1.12)

where the variable $y$ is considered slow relative to the fast variable $x$,
which is a consequence of the order of magnitude difference in the right hand
sides of (1.12) caused by the singular perturbation parameter $\epsilon\ll 1$.
The mathematical theory for such singular perturbation problems in standard
form is well established and a wealth of results for such standard (vdP-type)
relaxation oscillators have been derived with a variety of mathematical tools;
we refer the reader to the book by Kuehn [3] and the many references therein.

######  Remark 1.3.

We want to emphasise that the mathematical theory developed for the study of
singular perturbation problems does not depend on the special (standard) form
(1.12). Significant advances have been made in the study of non-standard
relaxation oscillators models; see e.g. [4, 5, 6, 7, 8].

In the 1960s Le Corbeiller [9], motivated by the study of electronic
oscillators, termed the vdP-type oscillations ‘four-stroke’ by reference to
the four distinct phases undergone in the relaxation cycle, in order to
distinguish these oscillations from ‘two-stroke’ oscillations, which consist
of only two distinct phases per cycle. Such two-stroke oscillations arise not
only in the context of electronic oscillators, but also in mechanical
oscillators with friction [10, 11, 12, 13, 14, 15], models of the trade cycle
in economics [16], aircraft-ground dynamics [17, 18], climate models [19],
discontinuous plastic deformation in metals [20] and cell-signalling models
[21]. We refer to Section 2 were we review some of these two-stroke oscillator
models.

Firstly, let us introduce a representative two-stroke oscillator model

$x^{\prime\prime}+\bigg{(}\frac{\epsilon}{1-x^{\prime}}-x^{\prime}\bigg{)}+x=0,\qquad
0<\epsilon\ll 1,$ (1.13)

which we recast as a dynamical system

$\displaystyle\begin{array}[]{lcl}x^{\prime}=1-y,\\\
y^{\prime}=x-1+y+\frac{\epsilon}{y}\,.\end{array}$ (1.16)

(a) Time trace.

(b) Phase space.

Figure 1.2: Two-stroke oscillation in (1.16) with $\epsilon=10^{-2}$.

Figure 2(a) shows the observed two-stroke oscillations in the time trace of
$y(t)$, i.e. the limit cycle consist of two distinct phases: a ‘static’ and a
‘dynamic’ phase. The distinct phases can be recognised as segments in phase
space relative to the characteristic

$R_{ts}(y):=1-y-\frac{\epsilon}{y}\,,$ (1.17)

with the static phase comprising the part of the limit cycle which follows
closely the part of the characteristic which asymptotes along $y=0$, and the
dynamic phase comprising the part of the limit cycle which is ‘off’ the
characteristic (or crosses it); see Figure 1.2. It is the $\epsilon$-dependent
rational term $\epsilon/y$ that causes this distinct two-stroke behaviour.
This term does not represent a uniformly small perturbation throughout phase
space, which distinguishes this model from standard perturbation problems such
as the vdP-type four-stroke oscillator models.

While Figure 1.2 mimics a slow-fast problem in phase space by following the
characteristic for a distinct part of the two-stroke oscillations, the
corresponding time trace shown in Figure 2(a) does not resemble an appreciable
separation of timescales, i.e. the static and dynamic phase evolve on
comparable timescales. The underlying ‘relaxation’ structure in this two-
stroke model (1.16) is revealed via a state-dependent time transformation, a
desingularisation,

$dt=yd\bar{t}$ (1.18)

which gives

$\displaystyle\begin{array}[]{lcl}x^{\prime}=(1-y)y,\\\
y^{\prime}=(x-1+y)y+\epsilon,\end{array}$ (1.21)

where with a slight abuse of notation the dash refers now to differentiation
with respect to the new time $\bar{t}$. System (1.21) is equivalent to system
(1.16) for $y>0$ (and up to a change of orientation for $y<0$).111Equivalence
is a basic topological concept in dynamical systems (see, e.g., [22, 23]) that
is very useful to resolve dynamics near certain types of singularities such as
poles of rational functions. Importantly, we have obtained a polynomial vector
field in (1.21) including a uniformly small perturbation term, and this system
produces now relaxation-type two-stroke oscillations as shown in Figure 1.3,
i.e. the static and dynamic phase of the two-stroke oscillator can now be
clearly identified as slow and fast segments in the corresponding time traces.
Hence, the desingularisation (1.18) has allowed us to extract a singular
perturbation problem in the form of a slow-fast system which preserves the
dynamical properties of the original two-stroke oscillator (1.13).

Figure 1.3: Two-stroke relaxation oscillation in (1.21); compare with Figure
2(a).

We emphasise that the singular perturbation problem (1.21) is not given in the
standard form (1.12). While we can clearly distinguish slow and fast motions,
there is no distinction between slow and fast variables in this model, i.e.
both time traces shown in Figure 1.3 consist of slow and fast segments. System
(1.21) is part of a more general class of perturbation problems,

$\begin{pmatrix}{x^{\prime}}\\\
{y^{\prime}}\end{pmatrix}=\begin{pmatrix}N_{1}(x,y)\\\
N_{2}(x,y)\end{pmatrix}f(x,y)+\epsilon\begin{pmatrix}G_{1}(x,y,\epsilon)\\\
G_{2}(x,y,\epsilon)\end{pmatrix}\,,$ (1.22)

with the specific choice of

$\begin{pmatrix}N_{1}(x,y)\\\ N_{2}(x,y)\end{pmatrix}=\begin{pmatrix}1-y\\\
x-1+y\end{pmatrix}\,,\quad
f(x,y)=y\,,\quad\begin{pmatrix}G_{1}(x,y,\epsilon)\\\
G_{2}(x,y,\epsilon)\end{pmatrix}=\begin{pmatrix}0\\\ 1\end{pmatrix},$ (1.23)

for system (1.21). Systems of the form (1.22) are identified as singular
perturbation problems whenever the set

$S_{0}=\\{z\in\mathbb{R}^{2}:N(z)f(z)=0\\}$ (1.24)

contains a smooth one-dimensional submanifold of $\mathbb{R}^{2}$, which
reflects the geometric definition of a singular perturbation problem.222A
general treatment of singular perturbation theory beyond the standard form in
arbitrary dimensions can be found in [8, 24]. The main aim of this article is
to consider this more general class of singular perturbation problems (1.22)
and provide a framework for the application of geometric singular perturbation
theory (GSPT) for which a (global) separation of slow and fast variables is
not required. This will allow us to prove existence, uniqueness and stability
results for a class of two-stroke relaxation oscillations, in a manner
conceptually analogous to the derivation of similar results in vdP-type
oscillators; see, e.g., [25].

In general, a standard form (1.12) of a singular perturbation problem where
the variables reflect the scale separation can only be achieved locally but
not globally. Thus from an application point-of-view it is desirable to
provide tools to analyse singularly perturbation problems of the more general
form (1.22) in a coordinate-independent manner. As pointed out by Fenichel in
his seminal work on GSPT [26], a global standard ‘form is not natural,
however, because it depends on the choice of special coordinates’, and he
clearly explains how to deal with a more general form. The general GSPT
framework we present here goes beyond Fenichel’s work and deals also with loss
of normal hyperbolicity, a necessary ingredient for two-stroke (or any
relaxation-type) oscillations.

######  Remark 1.4.

In the context of quasi steady-state reduction (QSSR) techniques commonly used
in chemical reaction systems, Goeke & Walcher [24] provide a general framework
that applies to systems (1.22) in the normally hyperbolic case as outlined by
Fenichel. Kaleda [27] shows the first results on two-stroke oscillations with
a focus on the bifurcation of a slow-fast separatrix loop. We would like to
point out that our results on relaxation oscillations (Section 4) can be
derived from a more general result found in [28], Theorem 2. We present an
alternative argument in this work in order to illustrate the scope, simplicity
and applicability of the methods developed herein.

The manuscript is structured as follows. In Section 2 we motivate our study by
introducing examples of two-stroke oscillations arising in models for
transistor oscillations and mechanical oscillations with friction, and show
how they can be appropriately phrased as singular perturbation problems in the
general form (1.22). In Section 3 we develop a framework for the application
of GSPT to this more general class of singular perturbation problems (1.22),
with an emphasis on the presentation of global, relaxation oscillatory results
from classic GSPT in a coordinate-independent setting. In Section 4 we define
minimal assumptions on singular perturbation problems in the general form
(1.22) that enable two-stroke relaxation oscillations, and prove existence,
uniqueness and stability results. Section 5 collects and contrasts dynamic
features associated with two- and four-stroke relaxation oscillations, and we
discuss possible transitions between two and four-stroke relaxation
oscillations. We also discuss the onset of two-stroke oscillations in a
mechanical oscillator. Finally, in Section 6 we conclude and outline future
work.

##  2 Two-stroke oscillators in applications

We motivate our work by providing examples of two-stroke oscillators which we
draw from the study of nonlinear transistor oscillators and mechanical
oscillators with friction.

###  2.1 An electronic two-stroke oscillator model

In [29], the author uses the Ebers-Moll large-signal approximation to show
that a class of nonlinear transistor oscillators including tuned-collector,
tuned-base, and Hartley transistor oscillators, can be described by a Lord-
Rayleigh type equation

$\ddot{x}+r_{em}(\dot{x})+x=0,$ (2.1)

where $x$ denotes (dimensionless) ‘current’,

$r_{em}(\dot{x})=-\mu\left(e^{a\dot{x}}-\kappa e^{(a+b)\dot{x}}\right),$

$\mu,\kappa,a,b>0$ are positive constants, and we assume

$\kappa<\frac{a}{a+b}<1$ (2.2)

is sufficiently small. Equation (2.1) can be recast as a dynamical system

$\displaystyle\begin{array}[]{lcl}\dot{x}=-\bar{y},\\\
\dot{\bar{y}}=x-\bar{R}_{em}(\bar{y}),\end{array}$ (2.5)

with nonlinear characteristic

$\bar{R}_{em}(\bar{y})=-r_{em}(-\bar{y})=\mu e^{-a\bar{y}}(1-\kappa
e^{-b\bar{y}}),$ (2.6)

This characteristic $\bar{R}_{em}(\bar{y})$ has a (unique) turning point at

$(x_{\ast},{y}_{\ast})=\left(\mu\bigg{(}\frac{a}{(a+b)\kappa}\bigg{)}^{a/b}\bigg{(}\frac{b}{a+b}\bigg{)},-\frac{1}{b}\ln\bigg{(}\frac{a}{(a+b)\kappa}\bigg{)}\right)\,,$
(2.7)

where ${y}_{\ast}<0$ due to (2.2). The characteristic approaches zero as
$\bar{y}\to\infty$ and grows exponentially towards $-\infty$ for
$\bar{y}<{y}_{\ast}$. For later convenience, we make a coordinate change
$\bar{y}=y+y_{\ast}$ which shifts the position of the turning point of the
characteristic to the $x$-axis, i.e. we obtain

$\displaystyle\begin{array}[]{lcl}\dot{x}=-y-{y}_{\ast},\\\
\dot{y}=x-R_{em}(y).\end{array}$ (2.10)

with characteristic

$R_{em}(y)=x_{*}\bigg{(}\frac{a+b}{b}\,e^{-ay}-\frac{a}{b}\,e^{-(a+b)y}\bigg{)}\,.$
(2.11)

Figure 1(a) shows the observed two-stroke oscillations in system (2.10), i.e.
the time trace $y(t)$ consists of a static and a dynamic phase. The
corresponding phase portrait, Figure 1(b), shows that the limit cycle follows
closely the characteristic in the static phase while it is off the
characteristic in the dynamic phase, i.e. it shows the same qualitative
features of the two-stroke oscillator shown in Figure 1.2. The main
(mathematical) difference between the part of the two characteristics (2.11)
and (1.17) that determine the static phase of the limit cycle, is exponential
growth versus unlimited growth due to a pole of a rational function.

(a) Time trace.

(b) Phase space.

Figure 2.1: Two-stroke oscillation in an Ebers-Moll model for a nonlinear
transistor (2.10), with parameter values $(\mu,\kappa,a,b)=(1,10^{-2},4,6)$.

Figure 2.2: Comparison of the approximation (2.12) (green), and the actual
characteristic (2.11) (blue). For the characteristic (2.11), the same
parameter values as Figure 2.1 were used; note the turning point is located at
$(x_{*},0)$. For the characteristic (2.12) we used $\epsilon=10^{-3}$.

This observation motivates us to model and replace the exponential growth of
the characteristic (2.11) for $y<0$ by a rational term, and we define the
following approximation

$R_{em,\epsilon}(y)=x_{*}e^{-ay}-\frac{\epsilon}{y},$ (2.12)

where $0<\epsilon\ll 1$ is a sufficiently small parameter. The rational term
$-\epsilon/y$ approximates the lower branch near $y=0$ while the exponential
profile of the vertical asymptote is preserved by the $x_{*}e^{-ay}$ term;
Figure 2.2 compares the two characteristics (2.11) and (2.12).

######  Remark 2.1.

The characteristic (2.12) has two turning points for $\epsilon>0$: one at
$(x,y)\sim(x_{*},0)$, and the other at $(x,y)\sim(0,\hat{y})$, with
$\hat{y}>0$. Both characteristics (2.11) and (2.12) have a vertical asymptote
along $x=0$, though the asymptote is approached from the right in (2.11) and
from the left in (2.12). These minor differences in the characteristics have
no significant effect on the dynamics. For similar reasons, we also refrain
from shifting the characteristic (2.12) slightly to the right to remove
(unphysical, yet very small) negative values for $y>0$.

(a) Time trace.

(b) Phase space.

Figure 2.3: Two-stroke oscillation in system (2.15), for the same parameter
values as Figure 2.1 and $\epsilon=10^{-3}$; compare with Figure 2.1.

The corresponding approximating system is given by

$\displaystyle\begin{array}[]{lcl}x^{\prime}=-{y}_{\ast}-y,\\\
y^{\prime}=x-R_{em,\epsilon}(y),\end{array}$ (2.15)

and the observed two-stroke oscillation and corresponding limit cycle are
shown in Figure 2.3; compare with Figure 2.1.

Figure 2.4: Two-stroke relaxation oscillation in (2.18) with the same
parameters as Figure 2.1 and $\epsilon=10^{-3}$. Compare with Figure 3(a), and
note the time-scale difference.

In order to reveal the underlying relaxation structure, we make the same
state-dependent time desingularisation $dt=y\,d\bar{t}$ as in (1.18), which
leads to the system

$\displaystyle\begin{array}[]{lcl}x^{\prime}=-({y}_{\ast}+y)y,\\\
y^{\prime}=(x-x_{\ast}e^{-ay})y+\epsilon,\end{array}$ (2.18)

where the dash notation now denotes differentiation with respect to $\bar{t}$.
System (2.18) is equivalent to (2.15) on $\\{y>0\\}$ (and up to a change of
orientation on $\\{y<0\\}$). The effect of the desingularisation (1.18) is
that we observe relaxation-type two-stroke oscillations in system (2.18); see
Figure 2.4. This system is in the general form of a singular perturbation
problem (1.22) with

$\begin{pmatrix}N_{1}(x,y)\\\
N_{2}(x,y)\end{pmatrix}=\begin{pmatrix}-y_{\ast}-y\\\
x-x_{\ast}e^{-ay}\end{pmatrix}\,,\quad
f(x,y)=y\,,\quad\begin{pmatrix}G_{1}(x,y,\epsilon)\\\
G_{2}(x,y,\epsilon)\end{pmatrix}=\begin{pmatrix}0\\\ 1\end{pmatrix}\,.$ (2.19)

Section 3 will provide the GSPT tools for the analysis of the observed two-
stroke relaxation oscillations.

###  2.2  A stick-slip oscillator model

Mechanical systems commonly exhibit two-stroke oscillations as a consequence
of the so-called stick-slip phenomenon due to friction. The corresponding
observations range from earthquake faulting or the sound of a violin, to the
unwanted screeching of chalk on the chalk-board or sliding of machine parts
(see [10] and references therein).

Many of the key dynamical features occurring in such complex mechanical
systems with friction are captured by the simple spring-mass system shown in
Figure 2.5, in which one considers a mass $m$ on a conveyor moving with
constant velocity $v_{0}$, attached to a wall by a spring of stiffness $k$.
For small displacements $x$, the mass moves with the conveyor. In this case,
the relative velocity between mass and belt is $v_{r}=\dot{x}-v_{0}=0$. This
is referred to as the ‘stick mode’, or ‘static phase’. As the mass moves with
the belt, the restoring force of the spring increases linearly in accordance
with Hookes law, and the mass starts to slip once this restoring force
balances the maximum static friction: this is the ‘stick-slip transition’.
Once slipping begins, we are in the ‘slip mode’ or ‘dynamic phase’. Finally,
the spring counteracts the sliding motion until static friction takes hold
again, and the process starts over.

Figure 2.5: Schematic of a spring-mass oscillator.

This mechanical system is described by the following (dimensionless) equation
of motion

$\ddot{x}+F_{f}(v_{r})+x=0,$ (2.20)

where $x$, $t$ and $v_{r}=\dot{x}-v_{0}$ denote displacement, time and
relative velocity, respectively. The stick-slip problem is usually modelled as
a discontinuous system, with static and dynamic phases treated independently
in accordance with

$F_{f}(v_{r})=\begin{cases}-x,&v_{r}=0,\\\ sgn(v_{r})\mu(v_{r}),&v_{r}\neq
0,\end{cases}$ (2.21)

where the function $\mu(v_{r})$ denotes the coefficient of friction. The
transition from stick to slip is determined by the stiction law, which asserts
that the stick phase ($v_{r}=0$) is maintained as long as

$|F_{f}(0)|=|x|\leq\mu_{s},$ (2.22)

where $\mu_{s}$ denotes the maximal value of static friction capable of
preventing the onset of the slipping motion. The function $\mu(v_{r})$ defines
the dynamic friction law, while the $sgn(v_{r})$ term ensures that the force
due to friction opposes the direction of relative motion. As is typical, we
assume $\mu(v_{r})$ is an even and strictly positive function, noting that in
general the specific form depends on the application. It is crucial for the
existence of oscillations that $F_{f}^{\prime}(v_{r})<0$ for small $|v_{r}|$.
This is known as the negative slope criterion in the stick-slip literature,
and is observed in many applications exhibiting the so-called ‘Stribeck
effect’ [10]. Following [15], we consider the two forms for $\mu(v_{r})$ shown
in Figure 2.6 as important examples:

  * •

$\mu(v_{r})$ decays exponentially toward a minimum value $\mu_{m}$, as in
Figure 6(a). Such a dependence is typical for the case of dry friction between
solid surfaces, and an appropriate form is suggested in [10] (see also [15])
as follows,

$\mu(v_{r})=\mu_{m}+(\mu_{s}-\mu_{m})e^{-a|v_{r}|},$ (2.23)

where $a>0$ is a fitting parameter which controls the slope of the
characteristic.

  * •

$\mu(v_{r})$ is modelled as a polynomial function which decays initially
towards a minimum value $\mu_{m}$ at relative velocity $v_{m}$, and increases
for $|v_{r}|>v_{m}$, as in Figure 6(b). This is suitable for systems
exhibiting dry friction for small $|v_{r}|$, and liquid or ‘viscous’ friction
for larger $|v_{r}|$. An appropriate form appears in, e.g. [14, 30, 31] (see
also [15, 32]):

$\mu(v_{r})=\mu_{s}-\frac{3(\mu_{s}-\mu_{m})}{2v_{m}}|v_{r}|+\frac{(\mu_{s}-\mu_{m})}{2v_{m}^{3}}|v_{r}|^{3}\,.$
(2.24)

(a) Exponential-type.

(b) Polynomial-type.

Figure 2.6: Friction-velocity curves for exponential and polynomial-type laws
in (a) and (b) respectively. The thick line for $v_{r}=0$, $\mu<\mu_{s}$
indicates the regime where the static friction dominates.

######  Remark 2.2.

Note that positivity of $\mu(v_{r})$ and the presence of the $sgn(v_{r})$ term
in (2.21) leads to a jump discontinuity at $v_{r}=0$. This (discontinuous)
modelling choice serves as an approximation of the real (smooth) mechanical
system. Two-stroke cycles in such discontinuous models are obtained by phase
space constructions in which segments from static and dynamic phases are
concatenated, recalling that the transition from stick to slip is enforced
when the threshold in (2.22) is reached; see e.g. [10, 13, 15, 32] and the
references therein.

In the absence of external forcing, the velocity of the mass never exceeds the
belt velocity, and so stick-slip oscillations occur only for $v_{r}\leq 0$
[14]. We intend to ‘smooth’ or ‘regularise’ the discontinuous dynamical system
(2.20) in the domain relevant for our analysis, i.e for $v_{r}\leq 0$. This is
achieved by replacing (2.21) with the $\epsilon$-dependent characteristic

$F_{f,\epsilon}(v_{r})=sgn(v_{r})\mu_{\epsilon}(v_{r})=sgn(v_{r})\bigg{(}\mu(v_{r})-\frac{\epsilon}{|v_{r}|}\bigg{)},\qquad
0<\epsilon\ll 1.$ (2.25)

This approximation effectively smooths out the corner at
$(v_{r},\mu(v_{r}))=(0,\mu_{s})$ and incorporates an asymptote along the line
$v_{r}=0$. Figure 2.7 shows that a sufficiently accurate approximation is
obtained in the case of both polynomial and exponential-type characteristics.
We emphasise that the regularisation in (2.25) is valid only for $v_{r}<0$ (or
$v_{r}>0$). In fact, this smooth approximation leads to a ‘stick mode’
$v_{r}\approx 0$ since the asymptotic speed $v_{r}=0$ cannot be reached for
$0<\epsilon\ll 1$.

######  Remark 2.3.

For background on the theory of regularisation in non-smooth problems, we
refer the reader to [33, 34] and the many references therein. Our choice of
regularisation (2.25) valid only over the half-plane $v_{r}<0$ simplifies the
analysis of this autonomous ‘stick-slip’ problem sufficiently. On the other
hand, we would like to emphasise that non-smooth problems involving external
forcing or bifurcations occurring near $v_{r}=0$ require regularisations which
are valid in an entire neighbourhood of $v_{r}=0$, and the analysis becomes
(necessarily) more involved. We refer the reader to, e.g., [11, 35, 36, 37,
38, 39] for rigorous treatments of such problems using techniques from GSPT
and blow-up.

(a)

(b)

Figure 2.7: In (a): the characteristic (2.25) with (2.23) plotted against the
exponential-type (piecewise) characteristic. Parameters:
$(\epsilon,v_{0},\mu_{m},\mu_{s},a)=(10^{-3},0.5,1,2,3)$. In (b): the
characteristic (2.25) with (2.24) plotted against the polynomial-type
(piecewise) characteristic. Parameters:
$(\epsilon,v_{0},v_{m},\mu_{m},\mu_{s})=(10^{-3},0.25,1,0.5,1)$. In each case
the original characteristic is plotted in blue, and the approximating
characteristic in green. Note that $y=-v_{r}$.

We introduce a new variable $y=-v_{r}=v_{0}-\dot{x}$ in the friction model
(2.20) to obtain the following dynamical system,

$\displaystyle\begin{array}[]{lcl}\dot{x}=v_{0}-y,\\\
\dot{y}=x-\mu(y)+\frac{\epsilon}{y},\end{array}$ (2.28)

where we used the smooth approximation (2.25) of the friction characteristic
(2.21) restricted to $v_{r}<0$ (i.e. $y>0$) together with the symmetry
$\mu(y)=\mu(-y)$. Again, after making the same time desingularisation
$dt=y\,d\bar{t}$ as in (1.18) we obtain the system

$\displaystyle\begin{array}[]{lcl}x^{\prime}=(v_{0}-y)y,\\\
y^{\prime}=(x-\mu(y))y+\epsilon,\end{array}$ (2.31)

where dash denotes differentiation with respect to the new time $\bar{t}$.
System (2.31) is in the general form of a singularly perturbed system (1.22)
with

$N(x,y)=\begin{pmatrix}{v_{0}-y}\\\ {x-\mu(y)}\end{pmatrix},\qquad
f(x,y)=y,\qquad G(x,y,\epsilon)=\begin{pmatrix}{0}\\\ {1}\end{pmatrix},$
(2.32)

and it is equivalent to system (2.28) in the relevant domain $y>0$. With the
specific choice of $\mu(v_{r})$ given by (2.23) or (2.24) we observe two-
stroke oscillations in (2.28), and these show up as two-stroke relaxation
oscillations in the corresponding system (2.31). Figure 2.8 shows the time
trace and corresponding relaxation cycle for the polynomial-type
characteristic (2.24).

######  Remark 2.4.

The exponential-type characteristic (2.23) has exactly the same features as
the characteristic (2.12) describing two-stroke relaxation oscillation in the
transistor oscillator model. Consequently, the time trace and corresponding
relaxation cycle are similar to those in Figures 2.4 and 3(b), respectively.

(a) Time trace.

(b) Phase space.

Figure 2.8: Two-stroke relaxation oscillation in (2.31) with polynomial-type
dynamic friction (2.24), with parameters identical to those in Figure 7(b).

######  Remark 2.5.

The dynamic friction coefficients (2.23) and (2.24) have a common
linearisation near $v_{r}=0$,

$\mu(v_{r})=a-b|v_{r}|,$

with positive coefficients $a,b$. Since the representative model (1.21) is in
the form (2.31) with

$\mu(y)=1-y,\qquad v_{0}=1,$ (2.33)

it can be viewed as a local minimal ‘stick-slip’ model. More importantly, it
serves as a mathematical ‘canonical’ model for two-stroke relaxation
oscillations.

##  3 A general GSPT framework

We consider planar perturbation problems of the general form

$z^{\prime}=H(z;\epsilon)=N(z)f(z)+\epsilon G(z;\epsilon),\qquad
z\in\mathbb{R}^{2},\qquad 0<\epsilon\ll 1,$ (3.1)

where $N(z)=(N_{1}(z),N_{2}(z))^{T}$ and
$G(z;\epsilon)=(G_{1}(z;\epsilon),G_{2}(z;\epsilon))^{T}$ are sufficiently
smooth vector fields, $f:\mathbb{R}^{2}\to\mathbb{R}$ is a sufficiently smooth
function, and the ′ notation denotes differentiation with respect to (fast)
time $t$. We will frequently denote variables componentwise by $z=(x,y)^{T}$.
We also consider the equivalent problem on a slow timescale $\tau=\epsilon t$:

$\dot{z}=\frac{1}{\epsilon}H(z;\epsilon)=\frac{1}{\epsilon}N(z)f(z)+G(z;\epsilon),$
(3.2)

where the dot notation denotes differentiation with respect to (slow) time
$\tau=\epsilon t$. Notice that (3.1) and (3.2) are equivalent for
$\epsilon>0$, but not in the limit $\epsilon\to 0$. In order to study systems
of form (3.1) respectively (3.2), we require coordinate independent analogues
for all the notions of standard GSPT [3, 40].

######  Remark 3.1.

We emphasise that standard form problems (1.12) can always be written in the
general form (3.1):

$\begin{pmatrix}{x^{\prime}}\\\
{y^{\prime}}\end{pmatrix}=\begin{pmatrix}{f(x,y,\epsilon)}\\\ {\epsilon
g(x,y,\epsilon)}\end{pmatrix}=\begin{pmatrix}{1}\\\
{0}\end{pmatrix}f_{0}(x,y)+\epsilon\begin{pmatrix}{f_{R}(x,y,\epsilon)}\\\
{g(x,y,\epsilon)}\end{pmatrix},$ (3.3)

where $f_{0}(x,y)=f(x,y,0)$. In the remainder of this work, we refer to (3.3)
when referring to problems in the ‘standard form’.

The converse, however, is not true. In particular, if the vector field $N$
contains isolated singularities (assumed to be bounded away from the set
$\\{f=0\\}$) then there is no global flowbox transformation rectifying orbits
of $N$ and hence, no global transformation of (3.1) into standard form (3.3).
All our example systems, (1.21), (2.18) and (2.31), are of this kind.

###  3.1 Layer problem

######  Definition 3.1.

(Layer problem). The system

$z^{\prime}=H(z,0)=h(z)=N(z)f(z)$ (3.4)

obtained from (3.1) in the limit $\epsilon\to 0$ is called the layer problem.

We impose the following on the system (3.4) for the remainder of this work.

######  Assumption 3.1.

The set of equilibria for the layer problem (3.4) takes the form of a disjoint
union

$S_{0}=S\cup V_{0},$

where

$S=\big{\\{}z\in\mathbb{R}^{2}|f(z)=0\big{\\}}$ (3.5)

is a one-dimensional regularly embedded submanifold of $\mathbb{R}^{2}$ called
the critical manifold, and

$V_{0}=\big{\\{}z\in\mathbb{R}^{2}|N(z)=(0,0)^{T}\big{\\}}$ (3.6)

is the (possibly empty) set of singularities of the vector field $N$. Without
loss of generality, we assume that $V_{0}$ contains only isolated
singularities.

######  Remark 3.2.

The existence of a smooth critical manifold $S$ (3.5) defines the system (3.1)
as a singular perturbation problem in the GSPT sense; see [26]. Accordingly,
we will often refer to $\epsilon\to 0$ as the singular limit. The assumption
that $V_{0}$ is discrete simplifies the analysis without being too restrictive
for our purposes.

The GSPT literature also deals with the case that the ‘critical manifold’ $S$
is an immersed submanifold, allowing for self-intersections. By assuming that
$S$ is a regularly embedded submanifold in Assumption 3.1, we are ruling out
the possibility of self-intersections. Thus, $Df$ is assumed to be non-
vanishing $\forall z\in S$ in this work.

######  Lemma 3.1.

Orbits of the layer problem (3.4) coincide with orbits of the auxiliary system
$z^{\prime}=N(z)$ on $\mathbb{R}^{2}\setminus S$.

###### Proof.

These systems are equivalent via the time desingularisation $dt=f(z)d\bar{t}$
modulo a reversal of orientation on $\\{z\in\mathbb{R}^{2}|f(z)<0\\}$. ∎

Evaluating the Jacobian for the system (3.4) along $S$ gives

$Dh\big{|}_{S}=NDf\big{|}_{S}=\begin{pmatrix}{N_{1}D_{x}f}&{N_{1}D_{y}f}\\\
{N_{2}D_{x}f}&{N_{2}D_{y}f}\end{pmatrix}\Bigg{|}_{S},$

which has a single trivial eigenvalue $\lambda_{0}=0$ since $\det
Dh\big{|}_{S}=0$, and a single non-trivial eigenvalue at $z\in S$ given by

$\lambda(z)=\text{Tr}\big{(}NDf\big{|}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}z}}\big{)}=\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla
f},N\rangle\big{|}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>(3.7)

The corresponding eigenspace of the trivial eigenvalue spans the tangent space
$T_{z}S$ at $z\in S$, i.e. it is orthogonal to the gradient of $f$,

$(T_{z}S_{n})^{\perp}=\text{span}\,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla
f}\big{|}_{z}\,.$

The corresponding eigenspace of the nontrivial eigenvalue is spanned by
$N(z)$, since

$(Dh)N=(NDf)N=N\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f,N\rangle=N\lambda\,,\qquad\forall
z\in S\,.$

######  Definition 3.2.

(Normal Hyperbolicity). We say that $z\in S$ is normally hyperbolic if the
non-trivial eigenvalue $\lambda(z)\neq 0$, and likewise call any submanifold
$S_{n}\subseteq S$ normally hyperbolic if $\lambda(z)\neq 0$, $\forall z\in
S_{n}$. A normally hyperbolic submanifold $S_{n}$ is called attracting if
$\lambda(z)<0$ $\forall z\in S_{n}$, and repelling if $\lambda(z)>0$, $\forall
z\in S_{n}$.

Let $S_{n}\subseteq S$ be a normally hyperbolic submanifold of $S$ and $z\in
S_{n}$. The inner product in (3.7), which is non-zero $\forall z\in S_{n}$,
induces the pointwise splitting

$T_{z}\mathbb{R}^{2}\big{|}_{S_{n}}=T_{z}S_{n}\oplus\mathcal{N}_{z}.$

Here $\mathcal{N}_{z}$ denotes the linear transverse fiber which has base at
$z\in S_{n}$ and is spanned by $N(z)$. The collection of all such fibers forms
a linear transverse fiber bundle $\mathcal{N}$, leading to the splitting

$T\mathbb{R}^{2}\big{|}_{S_{n}}=TS_{n}\oplus\mathcal{N},$ (3.8)

where $TS_{n}$ is the tangent bundle similarly obtained by collecting all the
tangent spaces $T_{z}S_{n}$ with base at $z\in S_{n}$. We note that more
generally, one can consider the existence of such a splitting as the defining
property of normal hyperbolicity.

######  Definition 3.3.

(Contact point). A point $F\in S$ such that

$\lambda(F)=\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla
f},N\rangle\big{|}_{F}=0$ (3.9)

is called a contact point.

A contact point $F\in S$ indicates a loss of normal hyperbolicity of the
critical manifold $S$.

######  Definition 3.4.

(Contact order). Let $U\subset\mathbb{R}^{2}$ be a neighbourhood of a contact
point $F$ and denote by $\mathcal{F}$ a segment of the corresponding layer
orbit through $F$ in $U$. Given local parameterisations $c(s)$ and $p(s)$,
$s\in[-s_{0},s_{0}]$, of $\mathcal{F}$ and $S$, respectively, such that
$c(0)=p(0)=F$. The layer flow has contact order $\sigma_{F}\in\mathbb{N}^{+}$
at $F$ with $S$ if

$c^{(m)}(0)=p^{(m)}(0)\text{ for }m=1,\ldots\sigma_{F},\qquad
c^{(\sigma_{F}+1)}(0)\neq p^{(\sigma_{F}+1)}(0).$

A contact point $F\in S$ is a point of tangency between the layer flow and
$S$, and the contact order $\sigma_{F}$ describes the degree of the tangency.

######  Proposition 3.2.

Let $F\in S$ denote a contact point of system (3.1), and assume without loss
of generality that $S$ is given as a graph $y=\varphi(x)$ locally near $F$.
Then the contact order at $F$ is given by

$\sigma_{F}=\text{min}\big{\\{}n\in\mathbb{N}^{+}\big{|}D^{(n)}_{x}\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f(x,\varphi(x)),N(x,\varphi(x))\rangle\big{|}_{F}\neq
0\big{\\}}.$ (3.10)

###### Proof.

We begin by introducing coordinates for which $S$ is locally rectified along
the $x-$axis.

######  Lemma 3.3.

Given $z=(x,y)\in S$ such that $D_{y}f|_{z}\neq 0$. There exists a smooth
change of coordinates such that the system (3.1) can be written locally as

$\begin{split}\begin{pmatrix}{x^{\prime}}\\\
{u^{\prime}}\end{pmatrix}&:=\tilde{N}(x,u)u+\epsilon\tilde{G}(x,u,\epsilon)\\\
&=\begin{pmatrix}{N_{1}(x,M(x,u))}\\\
{\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f(x,M(x,u),N(x,M(x,u)){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0})}\rangle}\end{pmatrix}u+\epsilon\begin{pmatrix}{G_{1}(x,M(x,u),\epsilon)}\\\
{\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f(x,M(x,u)),G(x,M(x,u),\epsilon)\rangle}\end{pmatrix}\,.\end{split}$
(3.11)

The critical manifold $S=\\{(x,u)\in\mathbb{R}^{2}\,|\,u=0\\}$ is straightened
locally along the $x$-axis.

###### Proof.

Define $u=f(x,y)$, which has a locally well-defined inverse $y=M(x,u)$, since
$D_{y}f|_{z}\neq 0$. Transforming (3.1) into the new $(x,u)$-coordinate system
gives the system (3.11). ∎

For system (3.11) with a flat manifold $S$, the condition for tangency between
the layer flow and $S$ at $F$ is

$\frac{du}{dx}\bigg{|}_{F}=\frac{\tilde{N}_{2}}{\tilde{N}_{1}}\bigg{|}_{F}=0\,,$

where $\tilde{N}_{2}|_{F}=\lambda(F)=0$ and $\tilde{N}_{1}|_{F}\neq 0$ (by
Assumption 3.1). The contact order is

$\sigma_{F}=\text{min}\bigg{\\{}n\in\mathbb{N}^{+}\bigg{|}D^{(n)}_{x}\bigg{(}\frac{du}{dx}\bigg{)}\bigg{|}_{F}\neq
0\bigg{\\}},$

which can be simplified by noting that

$D^{(n)}_{x}\bigg{(}\frac{du}{dx}\bigg{)}\bigg{|}_{F}=0\qquad\iff\qquad
D^{(n)}_{x}\tilde{N}_{2}(x,\varphi(x))\big{|}_{F}=0,$

and hence

$\sigma_{F}=\text{min}\big{\\{}n\in\mathbb{N}^{+}\big{|}D^{(n)}_{x}\tilde{N}_{2}(x,\varphi(x))\big{|}_{F}\neq
0\big{\\}}.$ (3.12)

In terms of the original coordinates, we have

$\sigma_{F}=\text{min}\big{\\{}n\in\mathbb{N}^{+}\big{|}D^{(n)}_{x}\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f(x,\varphi(x)),N(x,\varphi(x))\rangle\big{|}_{F}\neq
0\big{\\}},$

as required. ∎

######  Remark 3.3.

Lemma 3.3 holds locally about any $z\in S$, i.e. one can locally rectify $S$
near contact points or normally hyperbolic points.

######  Remark 3.4.

An order one contact point $F$ implies that the non-trivial eigenvalue
$\lambda$ (3.7) switches sign along $S$ as one ‘crosses’ $F$. The converse is
also true, i.e. a switch in stability of $S$ can occur only via a contact with
the layer flow.

######  Remark 3.5.

In standard form problems (3.3), order one contact points are fold points. A
fold point is generic if it satisfies the non-degeneracy and transversality
conditions,

$D^{2}_{x}f_{0}\big{|}_{F}\neq 0\qquad\text{and}\qquad
D_{y}f_{0}\big{|}_{F}\neq 0.$ (3.13)

Since

$D_{x}^{2}f_{0}\big{|}_{F}=D_{x}\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f,N\rangle\big{|}_{F}\neq
0,$

generic fold points are order one contact points, i.e. $\sigma_{F}=1$.
Conversely, if $F$ is an order one contact point in (3.3), then
$\lambda(F)=D_{x}f_{0}|_{F}=0$ which implies $D_{y}f_{0}|_{F}\neq 0$ (c.f.
Remark 3.2), and
$D_{x}\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f,N\rangle|_{F}=D_{x}^{2}f_{0}|_{F}\neq
0$, since $\sigma_{F}=1$. Hence order one contact points are generic folds in
standard form problems (3.3).

######  Application.

(Two-stroke relaxation oscillator model (2.31)). Recall that all our two-
stroke oscillator models introduced are of the general form (3.1) with $N(z)$,
$f(z)$ and $G(z)$ defined in (2.32). The distinguishing feature is the choice
of $v_{0}$ and $\mu(y)$ as highlighted in Table 3.1.

model | $v_{0}$ | $\mu(y)$
---|---|---
minimal model (1.23) | $1$ | $1-y$
electronic model (2.19) | $-y_{\ast}$ | $x_{\ast}e^{-ay}$
stick-slip model (2.32), exponential -type (2.23) | $v_{0}$ | $\mu_{m}+(\mu_{s}-\mu_{m})e^{-ay}$
stick-slip model (2.32), polynomial-type (2.24) | $v_{0}$ | $\mu_{s}-\frac{3(\mu_{s}-\mu_{m})}{2v_{m}}y+\frac{(\mu_{s}-\mu_{m})}{2v_{m}^{3}}y^{3}$
Table 3.1: The two-stroke oscillator (2.31) for the different models.

For all these models, the critical manifold is given by

$S=\big{\\{}(x,y)\in\mathbb{R}^{2}\;|\;y=0\big{\\}},$

and the set $V_{0}$ contains the single point

$p_{0}=\big{(}\mu(v_{0}),v_{0}\big{)}\,.$

Assumption 3.1 is satisfied for $v_{0}\neq 0$, since all models assume
$v_{0}>0$. The Jacobian of the layer problem at $p_{0}\notin S$ evaluates to

$Dh(p_{0})=DN(p_{0})f(p_{0})=\begin{pmatrix}{0}&{-v_{0}}\\\
{v_{0}}&{-v_{0}\mu^{\prime}(v_{0})}\end{pmatrix}.$

We have $\det Dh(p_{0})=v_{0}^{2}>0$, and the trace is given by

$\text{tr}Dh(p_{0})=-v_{0}\mu^{\prime}(v_{0}).$ (3.14)

The expression (3.14) is always positive for all models except the stick-slip
oscillator (2.31) with polynomial-type characteristic (2.24), for which (3.14)
is positive only for $v_{0}\in(0,v_{m})$. We restrict to values in this regime
in this work (see Section 5.2). Thus the equilibrium $p_{0}$ is an unstable
node or focus (in all cases). The non-trivial eigenvalue along $S$ is given by

$\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f,N\rangle\big{|}_{y=0}=x-\mu(0),$

and so the critical manifold decomposes into two normally hyperbolic branches

$S^{a}=\big{\\{}(x,0)\,\big{|}\,x<\mu(0)\big{\\}},\qquad
S^{r}=\big{\\{}(x,0)\,\big{|}\,x>\mu(0)\big{\\}},$

which are attracting and repelling, respectively. For all models, we have
$\mu(0)>0$. We also identify a single contact point

$F=(x_{F},y_{F})=(\mu(0),0),$

which we can classify as order one by noting that

$D_{x}\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f,N\rangle\big{|}_{F}=1\neq
0\qquad\implies\qquad\sigma_{F}=1.$ (3.15)

The layer problem dynamics near the contact point $F$ are sketched in Figure
3.1 for the case of the electronic two-stroke relaxation oscillator system
(2.18), which is obtained by substituting $\mu(y)=x_{*}e^{-ay}$ into (2.31);
see Table 3.1.

Figure 3.1: Singular limit dynamics near the jump-off point $F$ of system
(2.18). Parameter values are the same as Figure 2.2 with $\epsilon=0$.

###  3.2 Reduced problem

Consider system (3.2) which evolves on the slow timescale $\tau$. Taking the
singular limit $\epsilon\to 0$ becomes a non-trivial task:

$\dot{z}=\lim_{\epsilon\to
0}\bigg{(}\frac{1}{\epsilon}N(z)f(z)+G(z;\epsilon)\bigg{)}.$ (3.16)

Observe that in order for (3.16) to be well defined, the phase space must be
restricted to $S$, i.e. we need $f(z)$ to vanish, which implies that the
vector field of (3.16) must lie in the tangent bundle $TS$ of the critical
manifold $S$. Note that the vector $G(z;0)\in T\mathbb{R}^{2}$ does, in
general, not lie in $TS$. Hence one must determine the component of the vector
$G(z;0)$ in $TS$ to define the corresponding vector field.

Given a normally hyperbolic submanifold $S_{n}\subseteq S$, the existence of
the splitting (3.8) provides the means to define an appropriate vector field
in $TS_{n}$, because it guarantees the existence of a unique projection
operator

$\Pi^{S_{n}}_{N}:T\mathbb{R}^{2}\big{|}_{S_{n}}=TS_{n}\oplus\mathcal{N}\to
TS_{n}\,,$

i.e. $\Pi^{S_{n}}_{N}$ projects a vector with base point $z\in S_{n}$ along
$\mathcal{N}_{z}$ onto $T_{z}S_{n}$; see Figure 3.2. This allows for the
following definition, which is originally due to Fenichel [26] (see also
[24]).

######  Definition 3.5.

(Reduced Problem). For $z\in S_{n}$, the singular limit problem of system
(3.2) is defined by

$\dot{z}=\Pi_{N}^{S_{n}}\frac{\partial}{\partial\epsilon}H(z;\epsilon)\bigg{|}_{S_{n}\times\\{0\\}}=\Pi_{N}^{S_{n}}G(z;0)\bigg{|}_{S_{n}}\,,$
(3.17)

and it is called the reduced problem for (3.2).

######  Proposition 3.4.

The projection operator in (3.17) is given by

$\Pi_{N}^{S_{n}}=I_{2}-\frac{NDf}{\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f,N\rangle}\bigg{|}_{S_{n}}.$
(3.18)

An equivalent formulation of the reduced problem (3.17) is given by

$\dot{z}=\left[\frac{\det(N|G)}{\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f,N\rangle}\begin{pmatrix}{-D_{y}f}\\\
{D_{x}f}\end{pmatrix}\right]\Bigg{|}_{S_{n}},$ (3.19)

where $\det(N|G)$ denotes the determinant of the matrix with columns $N$ and
$G$.

Figure 3.2: Oblique projection of the vector $G(z;0)$ along $\mathcal{N}_{z}$
onto
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}T_{z}S_{n}}$
at a normally hyperbolic point $z\in S_{n}$.

###### Proof.

The form of the projection operator $\Pi_{N}^{S_{n}}$ reflects the definition
of an oblique projection as shown in Figure 3.2, i.e. based on the splitting
(3.8), the vector $G(z;0)$ has a unique representition

$G(z;0)=\Pi_{N}^{S_{n}}G(z;0)+\alpha N(z)$

where $\alpha\in\mathbb{R}$ is obtained by noting that

$\frac{||\text{proj}_{\nabla f(z)}G(z;0)||}{||\text{proj}_{\nabla
f(z)}N(z)||}=\alpha=\frac{\langle\nabla f(z),G(z;0)\rangle}{\langle\nabla
f(z),N(z)\rangle}.$

Note that $\text{proj}_{\nabla f(z)}N(z)\neq 0$ given normal hyperbolicity and
Assumption 3.1. Rearranging gives

$\begin{split}\Pi_{N}^{S_{n}}G(z;0)&=G(z;0)-\left(\frac{\langle\nabla
f(z),G(z;0)\rangle}{\langle\nabla f(z),N(z)\rangle}\right)N(z)\\\
&=\left[I_{2}-\frac{N(z)Df(z)}{\langle\nabla
f(z),N(z)\rangle}\right]G(z;0),\end{split}$

and the expression in (3.18) follows. Equivalence of (3.17) and (3.19) then
follows from

$\begin{split}\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f,N\rangle\big{|}_{S_{n}}\dot{z}&=\left[\big{(}\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f,N\rangle
I_{2}-NDf\big{)}G\right]\big{|}_{S_{n}}\\\
&=\left[\begin{pmatrix}{\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f,N\rangle-
N_{1}D_{x}f}&{-N_{1}D_{y}f}\\\
{-N_{2}D_{x}f}&{\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f,N\rangle-
N_{2}D_{y}f}\end{pmatrix}\begin{pmatrix}{G_{1}}\\\
{G_{2}}\end{pmatrix}\right]\Bigg{|}_{S_{n}}\\\
&=\left[\det(N|G)\begin{pmatrix}{-D_{y}f}\\\
{D_{x}f}\end{pmatrix}\right]\Bigg{|}_{S_{n}},\end{split}$

as required. ∎

######  Remark 3.6.

For a slow-fast system in standard form (3.3),

$\Pi_{N}^{S_{n}}=\left[\begin{pmatrix}{1}&{0}\\\
{0}&{1}\end{pmatrix}-\frac{1}{D_{x}f_{0}}\begin{pmatrix}{D_{x}f_{0}}&{D_{y}f_{0}}\\\
{0}&{0}\end{pmatrix}\right]\Bigg{|}_{S_{n}}=\begin{pmatrix}{0}&{-(D_{x}f_{0})^{-1}D_{y}f_{0}}\\\
{0}&{1}\end{pmatrix}\Bigg{|}_{S_{n}},$

and so the reduced problem is

$\begin{array}[]{lcl}\dot{x}=-(D_{x}f_{0})^{-1}(D_{y}f_{0})g(x,y,0),\\\
\dot{y}=g(x,y,0),\end{array}$ (3.20)

with $(x,y)\in S_{n}$. System (3.20) is the expression for the reduced vector
field for standard form slow-fast problems; see e.g. [3, 40].

######  Remark 3.7.

Since $D\\!f|_{z}\neq(0,0)$ $\forall z\in S$, it is clear from the form of
system (3.19) that equilibria in the reduced problem occur if and only if
$\det(N|G)=0$. Thus, the reduced dynamics on $S$ can be entirely characterised
via the scalar functions
$\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f,N\rangle$
and $\det(N|G)$.

Note that the projection operator (3.18) is not defined where the splitting
(3.8) breaks down, i.e. at contact points $F\in S$ where $S$ loses normal
hyperbolicity. In order to study the reduced problem (3.17) near contact
points $F$, we make a time desingularisation
$d\tau=-\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f,N\rangle|_{S}\,d\bar{\tau}$,
obtaining the desingularised problem,

$\dot{z}=-\left[\big{(}\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f,N\rangle
I_{2}-NDf\big{)}G\right]\big{|}_{S}=\left[\det(N|G)\begin{pmatrix}{D_{y}f}\\\
{-D_{x}f}\end{pmatrix}\right]\Bigg{|}_{S},$ (3.21)

where the overdot notation now denotes differentiation with respect to
$\bar{\tau}$. The desingularised problem (3.21) is equivalent to the reduced
problem (3.17) modulo a reversal of orientation when
$\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f,N\rangle|_{S}>0$,
i.e. on repelling submanifolds $S^{r}$ of $S$. Importantly, the desingularised
problem (3.21) is well defined in a neighbourhood of a contact point $F$,
which makes it a valuable tool for analysing the reduced problem in the case
of loss of normal hyperbolicity.

######  Definition 3.6.

(Regular Contact Point). Let $F\in S$ be a contact point with $\sigma_{F}=1$.
$F$ is called a regular contact point if it satisfies the following equivalent
conditions:

$\det(N|G)\big{|}_{(F,0)}\neq
0\qquad\iff\qquad\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f,G\rangle\big{|}_{(F,0)}\neq
0.$ (3.22)

(a) Jump-off

(b) Jump-on

Figure 3.3: The reduced flow near a regular contact point $F$ of order one.

A regular contact point implies that solutions of the desingularised problem
(3.21) reach a contact point $F$ in finite (forward or backward) time which
implies a finite (forward or backward) time blow-up of solutions near $F$ in
the reduced problem (3.17), i.e. solutions of the reduced problem cease to
exist. Since a regular contact point is also an order one contact point, the
stability property of the critical manifold changes near $F$; see Remark 3.4.
Hence, the reduced flow is either towards or away from a regular contact point
$F$ as shown in Figure 3.3.

######  Definition 3.7.

(Jump-Off/On Point). A regular contact point $F$ is called a jump-off point if
the reduced flow is towards $F$ (see, e.g., Figure 3(a)) or a jump-on point if
the reduced flow is away from $F$ (see, e.g., Figure 3(b)).

######  Remark 3.8.

Jump-off points play an important role in relaxation oscillations insofar as
they mark a transition from slow to fast motion, which is a defining feature
for relaxation oscillations.

######  Remark 3.9.

A point $F\in S$ in standard form problems (3.3) is called a regular fold
point if it satisfies $g|_{(F,0)}\neq 0$ in addition to the fold conditions
$f_{0}|_{F}=0$, $D_{x}f_{0}|_{F}=0$ and the nondegeneracy conditions in (3.13)
(see, e.g. [3]). Hence at a regular fold point $F$ we have

$g|_{(F,0)}\neq
0\qquad\implies\qquad\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f,G\rangle\big{|}_{(F,0)}=(D_{y}f_{0})g\big{|}_{(F,0)}\neq
0,$

i.e. regular fold points are regular contact points with $\sigma_{F}=1$.

Conversely, let $F$ be a regular contact point with $\sigma_{F}=1$ in (3.3).
Then by analogous reasoning as Remark 3.5 we have $D_{x}f_{0}|_{F}=0$,
$D_{y}f_{0}|_{F}\neq 0$ and $D_{x}^{2}f_{0}|_{F}\neq 0$. Moreover,
$\det(N|G)|_{(F,0)}=g|_{(F,0)}\neq 0$. Hence regular contact points with
$\sigma_{F}=1$ in system (3.3) are regular fold points.

######  Application.

(Two-stroke relaxation oscillator model (2.31) continued). For this model, we
obtain a projection operator (3.18) of the form

$\Pi_{N}^{S}=\begin{pmatrix}{1}&{0}\\\
{0}&{1}\end{pmatrix}-\frac{1}{x-\mu(0)}\begin{pmatrix}{0}&{v_{0}}\\\
{0}&{x-\mu(0)}\end{pmatrix}=\begin{pmatrix}{1}&{-v_{0}/(x-\mu(0))}\\\
{0}&{0}\end{pmatrix}\,,$ (3.23)

and a reduced problem (3.17),

$\begin{pmatrix}\dot{x}\\\
\dot{y}\end{pmatrix}=\Pi_{N}^{S}G(z;0)|_{S}=\begin{pmatrix}{-v_{0}/(x-\mu(0))}\\\
{0}\end{pmatrix}\,.$ (3.24)

Since $S$ loses normal hyperbolicity at the regular contact point $F$ for
$x=\mu(0)$, we study the corresponding desingularised problem (3.21),

$\begin{pmatrix}\dot{x}\\\ \dot{y}\end{pmatrix}=v_{0}\begin{pmatrix}1\\\
0\end{pmatrix}\,,$ (3.25)

with $v_{0}>0$. Thus $F$ is a regular jump-off point by Definition 3.3, i.e.
the reduced flow is towards $F$ and a finite (forward) time blow-up of
solutions occurs. The combined layer and reduced problem dynamics near the
jump-off point $F$ are shown in Figure 3.1 for the case of system (2.18).

###  3.3 Local GSPT results in non-standard form

While any general slow-fast problem (3.1) can be locally transformed to a
standard form problem (3.3) as shown in Appendix A, it is only a theoretical
result in nature, i.e. in most cases this transformation cannot be calculated
explicitly. In fact, it is not desirable to make any coordinate
transformations in applications when one can apply GSPT directly to a non-
standard problem (3.1). Thus in the following, we present local results which
are well-known in standard GSPT literature in their non-standard form.

######  Theorem 3.5.

(Fenichel Theorems [26, 41, 42]. See also [3, 43]). Consider (3.1) and suppose
$S_{n}\subseteq S$ is normally hyperbolic and compact. Then
$\exists\epsilon_{0}>0$ such that $\forall\epsilon\in(0,\epsilon_{0})$ the
following assertions are true:

  * (F1)

There exists a $C^{r}$-smooth locally invariant manifold $S_{n,\epsilon}$,
called the slow manifold, which is $\mathcal{O}(\epsilon)$ from $S_{n}$ in the
Hausdorff distance.

  * (F2)

The flow on $S_{n,\epsilon}$ converges to the reduced flow on $S_{n}$ as
$\epsilon\to 0$.

  * (F3)

The manifold $S_{n,\epsilon}$ is normally hyperbolic, and exhibits the same
stability properties with respect to the fast dynamics as $S_{n}$ (i.e.
$S_{n,\epsilon}$ is attracting or repelling if $S_{n}$ is attracting or
repelling respectively).

  * (F4)

The manifold $S_{n,\epsilon}$ is usually not unique, but all manifolds
satisfying (F1)-(F3) lie within a Hausdorff distance which is
$\mathcal{O}(e^{-K/\epsilon})$ for some constant $K>0$.

  * (F5)

Statements (F1)-(F4) also hold locally for the stable and unstable manifolds
given by the foliations

$W^{s}_{loc}(S_{n})=\bigcup_{z\in S_{n}}W^{s}_{loc}(z),\qquad
W^{u}_{loc}(S_{n})=\bigcup_{z\in S_{n}}W^{u}_{loc}(z),$

where $W^{s}_{loc}(z)$ and $W^{u}_{loc}(z)$ denote the local stable and
unstable one-dimensional manifolds associated with $z\in S_{n}$. The manifolds
$W^{s}_{loc}(S_{n})$ and $W^{s}_{loc}(S_{n})$ persist as manifolds
$W^{s}_{loc}(S_{n,\epsilon})$ and $W^{u}_{loc}(S_{n,\epsilon})$ with
foliations consisting of $C^{r}$-smooth leaves $W_{loc}^{s}(z)$ and
$W_{loc}^{s}(z)$ with base points $z\in S_{n,\epsilon}$. In particular, we
have the following:

    * (i)

$W^{s}_{loc}(S_{n,\epsilon})=\bigcup_{z\in
S_{n,\epsilon}}W^{s}_{loc}(z),\qquad W^{u}_{loc}(S_{n,\epsilon})=\bigcup_{z\in
S_{n,\epsilon}}W^{u}_{loc}(z).$ (3.26)

    * (ii)

The foliations in (3.26) are positively and negatively invariant respectively,
i.e. $W_{loc}^{s}(z)\cdot t\subset W_{loc}^{s}(z\cdot t)$ $\forall t\geq 0$
such that $z\cdot t\in S_{n,\epsilon}$, where $\cdot t$ denotes the forward
evolution of $z$ in time, and similarly $W^{u}_{loc}(z)\cdot t\subset
W^{u}_{loc}(z\cdot t)$ $\forall t\leq 0$ such that $z\cdot t\in
S_{n,\epsilon}$.

    * (iii)

If $S_{n}$ is attracting and $\lambda(z)<\alpha_{s}<0$ $\forall z\in S_{n}$,
there exists a constant $\kappa_{s}>0$ such that if $z\in S_{n,\epsilon}$ and
$v\in W^{s}_{loc}(z)$, then

$\|v\cdot t-z\cdot t\|\leq\kappa_{s}e^{\alpha_{s}t}$

$\forall t\geq 0$ such that $z\cdot t\in S_{n,\epsilon}$. Similarly, if
$S_{n}$ is repelling and $\lambda(z)>\alpha_{u}>0$ $\forall z\in S_{n}$, there
exists a constant $\kappa_{u}>0$ such that if $z\in S_{n,\epsilon}$ and $v\in
W^{u}_{loc}(z)$, then

$\|v\cdot t-z\cdot t\|\leq\kappa_{u}e^{\alpha_{u}t}$

$\forall t\leq 0$ such that $z\cdot t\in S_{n,\epsilon}$.

######  Remark 3.10.

Although Theorem 3.5 actually implies existence of an entire family of slow
manifolds, by (F4) all such slow manifolds are exponentially close in
$\epsilon$, so fixing a choice of slow manifold is rarely problematic in
calculations.

######  Remark 3.11.

Fenichel’s original work in [26] does not depend on the system being
expressible in the standard form (3.3), and can be applied directly for
general systems (3.1).

In the case that normal hyperbolicity breaks down, Theorem 3.5 no longer
applies. Thus, we still require a description of the perturbed dynamics near a
contact point $F\in S$. We consider here only the least degenerate case, i.e.
the dynamics near a regular contact point $F$. Let
$U_{F}\subset\mathbb{R}^{2}$ denote a neighbourhood of a regular jump-off
point $F=(x_{F},y_{F})$, and let $\mathcal{F}$ denote the layer problem orbit
segment in $U_{F}$ that has contact of order one with $S$ at $F$. We define
$\Sigma_{in/out}\in U_{F}$,

$\Sigma_{in}=\big{\\{}(x_{F}-\rho,y)|y\in
J\big{\\}}\,,\qquad\Sigma_{out}=\big{\\{}(x_{F}+\rho,y)|y\in J\big{\\}}\,,$
(3.27)

and assume without loss of generality that these vertical sections are
transverse to both $S$ and the layer flow for sufficiently small $\rho>0$ and
a suitably defined real interval $J$; see Figure 3.4.

Figure 3.4: Behaviour near a regular contact point $F$, as described by
Theorem 3.6.

######  Theorem 3.6.

Let $F\in S$ be a jump-off point of system (3.1), and assume without loss of
generality that $D_{y}f|_{F}\neq 0$. Then $\exists\epsilon_{0}>0$ such that
$\forall\epsilon\in(0,\epsilon_{0}]$ the following assertions hold:

  * (1)

The attracting slow manifold $S^{a}_{\epsilon}$ leaves the neighbourhood of
the contact point via $\Sigma_{out}$, and

$|y_{s}-y_{l}|=\mathcal{O}(\epsilon^{2/3})\,,$

where $y_{s}$ respectively $y_{l}$ denote the $y$-coordinate of
$S_{\epsilon}^{a}\cap\Sigma_{out}$ respectively $\mathcal{F}\cap\Sigma_{out}$.

  * (2)

The transition map $\pi:\Sigma_{in}\to\Sigma_{out}$ is a contraction with
contraction rate $\mathcal{O}(e^{-c/\epsilon})$, for some constant $c>0$.

###### Proof.

See Appendix B. ∎

######  Remark 3.12.

The distance $\mathcal{O}(\epsilon^{2/3})$ in Theorem 3.6 is consistent with
the result in [44], where it is shown that the slow manifold
$S_{\epsilon}^{a}$ leaves a neighbourhood of a regular fold point at a
Hausdorff distance which is $\mathcal{O}(\epsilon^{2/3})$ from the critical
fiber. The only difference between Theorem 3.6 and the result in [44] is that
the distance $\mathcal{O}(\epsilon^{2/3})$ in Theorem 3.6 is stated for
general systems (3.1) in terms of simple transversals in the original
coordinates.

##  4 Existence of two-stroke relaxation oscillations

In this section we present existence, uniqueness and stability results for
two-stroke relaxation oscillations in non-standard singular perturbation
problems (3.1) which apply to our model systems (1.21), (2.18), and (2.31).

First, we note that the presence of a regular contact point and the associated
finite time blow-up in the reduced problem allows one to concatenate segments
of layer and reduced problems.

######  Definition 4.1.

(Reciprocal point, cf. [27]). A point $L_{z}\in S$ is reciprocal to $z\in S$
if they belong to the endpoints of a heteroclinic orbit of the layer problem.
We also say that the pair $(z,L_{z})$ is reciprocal.

######  Remark 4.1.

For standard form problems (3.3), if $z=(x_{z},y_{z})\in S$ has a reciprocal
point $L_{z}=(x_{L_{z}},y_{L_{z}})$, then $y_{z}=y_{L_{z}}$ since the fast
fibers are parallel to the $x$-axis. For general systems (3.1) however, the
relationship between $z\in S$ and a corresponding reciprocal point $L_{z}$ is
non-trivial.

######  Definition 4.2.

(Singular relaxation cycle). A closed singular orbit consisting of at least
one segment from the layer and the reduced problem is called a singular
relaxation cycle; see Figure 4.1.

Figure 4.1: Two-stroke singular relaxation cycle
$\Gamma=\Gamma^{R}\cup\Gamma^{L}$ for a system (3.1) satisfying Assumption
4.1. Also shown are the segments $\Sigma_{i}$ defined for the proof of Theorem
4.1, and the relaxation oscillation $\Gamma_{\epsilon}$.

We are now able to provide minimal conditions on system (3.1) for the
existence of a singular two-stroke singular relaxation cycle and state the
main result about the persistence of a two-stroke relaxation cycle under
sufficiently small perturbations $\epsilon\ll 1$.

######  Assumption 4.1.

System (3.1) has the following properties:

  * ($A1$)

The critical manifold $S$ has precisely one jump-off point $F$, and therefore
decomposes $S=S^{a}\cup\\{F\\}\cup S^{r}$, where $S^{a}$ ($S^{r}$) is
attracting (repelling) and normally hyperbolic. Without loss of generality,
$D_{y}f|_{F}\neq 0$ is satisfied.

  * ($A2$)

The jump-off point $F$ has a reciprocal point $L_{F}\in S^{a}$, implying the
existence of a singular relaxation cycle $\Gamma=\Gamma^{L}\cup\Gamma^{R}$,
where the segment $\Gamma^{R}$ is a trajectory segment of the reduced problem
from $L_{F}\in S^{a}$ to $F$ (see Figure 4.1).

######  Theorem 4.1.

Given system (3.1) under Assumptions 3.1 and 4.1. Let $U$ denote a fixed
tubular neighbourhood of the singular relaxation cycle $\Gamma$ such that
$U\cap V_{0}=\O$. Then $\exists\epsilon_{0}>0$ such that
$\forall\epsilon\in(0,\epsilon_{0})$, system (3.1) has a unique two-stroke
relaxation cycle $\Gamma_{\epsilon}\subset U$. The relaxation cycle
$\Gamma_{\epsilon}$ is attracting with Floquet exponent bounded above by
$-K/\epsilon$ for some constant $K>0$, and converges to $\Gamma$ in Hausdorff
distance as $\epsilon\to 0$.

###### Proof.

Under Assumption 4.1 and assuming a sufficiently small tubular neighbourhood
$U$ of the singular relaxation cycle $\Gamma$, we define cross sections
$\Sigma_{i}\in U$, $i\in\\{1,2,3\\}$, as shown in Figure 4.1. Let
$\Pi:\Sigma_{3}\to\Sigma_{3}$ denote a global return map defined by the
composition

$\Pi=\Pi_{23}\circ\Pi_{12}\circ\Pi_{31}$

where $\Pi_{31}:\Sigma_{3}\to\Sigma_{1}$, $\Pi_{12}:\Sigma_{1}\to\Sigma_{2}$
and $\Pi_{23}:\Sigma_{2}\to\Sigma_{3}$ are the corresponding transition maps.
We show that $\Pi$ is a contraction, from which the existence of a unique
stable limit cycle follows by the contraction mapping principle.

First, consider the map $\Pi_{31}$ which describes the flow near the normally
hyperbolic attracting branch $S^{a}$ of the critical manifold $S$. By Fenichel
Theorem 3.5 (F1)-(F3), $S^{a}$ perturbs to an attracting slow manifold
$S_{\epsilon}^{a}$. Moreover, Theorem 3.5 (F5) implies that initial conditions
in $\Sigma_{3}$ are exponentially attracted to the slow manifold
$S_{\epsilon}^{a}$ with rate $\mathcal{O}(e^{-c_{1}/\epsilon})$ for some
$c_{1}>0$, and they follow their base points on $S^{a}_{\epsilon}$ until
exiting through $\Sigma_{1}$. Hence an (open) interval
$I_{3}\subset\Sigma_{3}$ of size $O(1)$ about $\Gamma\cap\Sigma_{3}$ is mapped
to an interval $\Pi_{31}(I_{3})\subset\Sigma_{2}$ of width
$\mathcal{O}(e^{-c_{1}/\epsilon})$ about $S_{\epsilon}^{a}\cap\Sigma_{1}$,
i.e. the map $\Pi_{31}$ is exponentially contracting.

Second, consider the map $\Pi_{12}$ which describes the flow past the contact
point $F$ where normal hyperbolicity is lost. This transition map is covered
by Theorem 3.6 which states that an (open) interval $I_{1}\subset\Sigma_{1}$
of size $O(1)$ about $\Gamma\cap\Sigma_{1}$ is mapped to an interval
$\Pi_{12}(I_{1})\subset\Sigma_{2}$ of width $\mathcal{O}(e^{-c_{2}/\epsilon})$
about $S_{\epsilon}^{a}\cap\Sigma_{2}$ for some $c_{2}>0$, i.e. the map
$\Pi_{12}$ is also exponentially contracting.

Third, consider the map $\Pi_{23}$ which describes regular flow from
$\Sigma_{2}$ to $\Sigma_{3}$ in $U$. The flow-box theorem (see e.g. [45, 46,
47]) and regular perturbation theory imply that $\Pi_{23}$ is a diffeomorphism
with at most algebraic growth, i.e. there exists an (open) interval
$I_{2}\subset\Sigma_{2}$ of size $O(1)$ about $\Gamma\cap\Sigma_{2}$ such that
$\Pi_{23}(I_{2})\subset\Sigma_{3}$.

Finally, we take the composition of the three transition maps. Since
$\Pi_{31}$ and $\Pi_{12}$ are exponentially contracting while $\Pi_{23}$ has
(at most) algebraic growth, the return map $\Pi(I)\subset I$ is an interval of
width $\mathcal{O}(e^{-c/\epsilon})$ for some $c>0$, i.e. $\Pi$ is a
contraction with rate $\mathcal{O}(e^{-c/\epsilon})$. This guarantees the
existence of a unique fixed point corresponding to a unique stable limit cycle
$\Gamma_{\epsilon}\in U$.

By Theorem 3.5 (F1)-(F2), $S^{a}_{\epsilon}$ is $\mathcal{O}(\epsilon)$ from
$S^{a}$ and converges to $S^{a}$ in the Hausdorff distance as $\epsilon\to 0$.
By Theorem 3.6 (i), $S_{\epsilon}^{a}\cap\Sigma_{2}$ is
$\mathcal{O}(\epsilon^{2/3})$ from the intersection of $\Sigma_{3}\cap\Gamma$
and converging to $\Gamma\cap\Sigma_{3}$ as $\epsilon\to 0$. Hence
$\Gamma_{\epsilon}$ approaches $\Gamma$ in the Hausdorff distance as
$\epsilon\to 0$. ∎

######  Application.

(Application of Theorem 4.1 to two-stroke oscillator model (1.16)).

(a)

(b)

(c)

Figure 4.2: Singular limit dynamics and singular relaxation cycle
$\Gamma=\Gamma^{L}\cup\Gamma^{R}$ for (a) system (1.21); (b) system (2.31)
with characteristic (2.23); (c) system (2.31) with polynomial-type
characteristic (2.24). Parameter values (except $\epsilon$) in (b) and (c) are
the same as in Figure 2.7.

System (1.16) is equivalent to system (1.21). The existence of a jump-off
point in system (1.21) was shown in Section 3. What is left to show is that
system (1.21) has a singular relaxation cycle $\Gamma$, i.e. it remains to
find a reciprocal point $L_{F}\in S^{a}$ of $F$. Note that the auxiliary layer
problem of system (1.21),
$(x^{\prime},y^{\prime})^{T}=N(x,y)=(1-y,x-1+y)^{T}$, is linear with an
unstable focus at $p_{0}=(0,1)$. Let $(x(t),y(t))$ denote the unique solution
for this auxiliary system with $(x(0),y(0))=(1,0)=F$. Since there are no other
singularities of $N(x,y)$, the expansion and rotation due to the unstable
focus at $p_{0}$ guarantees the existence of a (unique) reciprocal point
$L_{F}=(x(T),0)$ for some time $T>0$, where $x(T)<1$. Numerically we obtain an
estimate $L_{F}\approx(-11.2,0)$ (Figure 2(a)). Hence we can construct a
singular relaxation cycle $\Gamma=\Gamma^{R}\cup\Gamma^{L}$, where
$\Gamma^{R}=\big{\\{}(x,0)\in S^{a}|x\in[x(T),1]\big{\\}}$ and
$\Gamma^{L}=\\{(x(t),y(t))|t\in[0,T]\\}$. Hence Assumption 4.1 holds and by
Theorem 4.1 there exists a strongly attracting two-stroke relaxation cycle
$\Gamma_{\epsilon}$ converging to $\Gamma$ in the Hausdorff distance as
$\epsilon\to 0$. Thus we have proven existence and stability of the observed
two-stroke oscillations in system (1.16).

######  Remark 4.2.

We omit the details involved with proving the existence of reciprocal points
$L_{F}$ for system (2.31) with either polynomial or exponential-type
characteristics, which is clear from numerical simulations. See Figures 2(b)
and 2(c), where for our specific choice of parameters we find
$L_{F}\approx(-6.86,0)$ and $L_{F}\approx(-0.09,0)$, respectively.

For completeness, we also include results on singular relaxation cycles for
which $F$ is a jump-on point; see Figure 3(b).

######  Assumption 4.2.

System (3.1) has the following properties:

  * ($\bar{A1}$)

The critical manifold $S$ has precisely one jump-on point $F$, and therefore
decomposes $S=S^{a}\cup\\{F\\}\cup S^{r}$, where $S^{a}$ ($S^{r}$) is
attracting (repelling) and normally hyperbolic. Without loss of generality,
$D_{y}f|_{F}\neq 0$ is satisfied.

  * ($\bar{A2}$)

The jump-on point $F$ has a reciprocal point $L_{F}\in S^{r}$, implying the
existence of a singular relaxation cycle $\Gamma=\Gamma^{L}\cup\Gamma^{R}$,
where the segment $\Gamma^{R}$ is a trajectory segment of the reduced problem
from $L_{F}\in S^{r}$ to $F$.

######  Theorem 4.2.

Given system (3.1) under Assumptions 3.1 and 4.2. Let $U$ denote a fixed
tubular neighbourhood of the singular relaxation cycle $\Gamma$ such that
$U\cap V_{0}=\O$. Then $\exists\epsilon_{0}>0$ such that
$\forall\epsilon\in(0,\epsilon_{0})$, system (3.1) has a unique two-stroke
relaxation cycle $\Gamma_{\epsilon}\subset U$. The relaxation cycle
$\Gamma_{\epsilon}$ is repelling with Floquet exponent bounded below by
$K/\epsilon$ for some constant $K>0$, and converges to $\Gamma$ in Hausdorff
distance as $\epsilon\to 0$.

###### Proof.

Reversing time and applying the same arguments given in the proof of Theorem
4.1 yields the desired result. ∎

##  5 Comparison of two and four-stroke relaxation oscillators

We reiterate that the models presented in this work show that in general, two-
stroke oscillations can occur given an $S$-shaped characteristic with two
turning points (e.g. system (2.31) with polynomial-type characteristic), a
characteristic with one turning point and a vertical asymptote (e.g. system
(2.5)), or, minimally, a $C$-shaped characteristic with a single turning point
and no vertical asymptote (e.g. system (1.21)). Thus, in general, only one
turning point is necessary for two-stroke oscillation, in contrast to the
four-stroke oscillation typified by the vdP oscillator, which requires two.
These observations are summarised in Table 5.1.

| Two-stroke | Four-stroke
---|---|---
Characteristic | (S), (SC), or (C) | (S) only
# contact points | $1$ | $2$
# equilibria on $S$ | $0$ | $1$
# equilibria of $N$ | $1$ | $0$
Table 5.1: Comparison of minimal features for two and four-stroke relaxation
oscillation. Here the notation ($C$), ($SC$), ($S$) refers to the shape of the
characteristic: ($C$) means one turning point and no vertical asymptote (e.g.
system (1.21)); ($SC$) means one turning point and a vertical asymptote (e.g.
system (2.18)); ($S$) means two turning points (e.g. (1.5) or (2.31) with
polynomial-type characteristic).

We emphasise that the singularly perturbed two-stroke oscillator models
presented here cannot be globally put into the standard form (1.12), i.e. they
are genuine singular perturbation problems of the more general form (3.1).

We would like to point out, though, that two-stroke relaxation oscillations
have been studied in singularly perturbed problems in the standard form (1.12)
in the context of a model for aircraft-ground dynamics [18] (see also [17]),
as well as in a model of discontinuous plastic deformation of metals [20]. The
characteristics in these models have one turning point and a vertical
asymptote, similar to the characteristic of system (2.10) shown in Figure
1(b). What distinguishes these two-stroke relaxation oscillations from those
presented in this work is that the amplitude of the relaxation oscillation
approaches infinity in the singular limit $\epsilon\to 0$. The unbounded
growth in cycle amplitude in these systems is a consequence of the global
separation of slow and fast variables which implies a layer flow along
straight fast fibers and, hence, no return mechanism is possible on a compact
domain. Instead, the return mechanism in the case of the oscillators presented
in [17, 18] occurs due to a ‘loss of normal hyperbolicity at infinity’, i.e.
an alignment of the fast fibers with the vertical asymptote of the
characteristic. A similar feature occurs in an autocatolator model, which
exhibits three-timescale relaxation oscillations, and has been studied using
non-standard form techniques in detail in [5].

It is worth noting that system (2.10) could be analysed as a standard singular
perturbation problem by artificially introducing a perturbation parameter
$\epsilon\ll 1$ so that (2.10) becomes

$\displaystyle\begin{array}[]{lcl}\dot{x}=\epsilon(-y-{y}_{\ast}),\\\
\dot{y}=x-R_{em}(y),\end{array}$ (5.3)

where $R_{em}(y)$ is given by (2.11). System (5.3) also ‘loses normal
hyperbolicity at infinity’ due to an alignment of the critical curve with the
fast fibers (see [17] for details on how to deal with this). This approach has
the advantage that the characteristic does not need to be approximated, but
the disadvantage of extreme sensitivity of the cycle amplitude as a function
of $\epsilon$ (as mentioned above, the amplitude tends to infinity as
$\epsilon\to 0$). Furthermore, one must be able to understand perturbations
all the way up to $\epsilon=1$ (homotopy argument) in order to completely
capture the original two-stroke oscillation observed in (2.10).

Table 5.2 categorises the different approaches, showing which characteristic
types allow for a standard form analysis, whether it can be undertaken on a
compact domain, and those cases for which one must move beyond the standard
form. Note that two-stroke relaxation oscillation in systems with a $C$-shaped
characteristic cannot occur in standard form problems.

| Two-stroke (C) | Two-stroke (SC) | Two-stroke (S) | Four-stroke (S)
---|---|---|---|---
Stnd compact | N | N | N | Y
Stnd noncompact | N | Y | N | Y
Non-stnd compact | Y | Y | Y | Y
Table 5.2: Different relaxation oscillation types, and whether or not they can
be analysed as standard form problems on a compact or non-compact domain. As
in Table 5.1, the notation ($C$), ($SC$), ($S$) refers to the shape of the
characteristic

###  5.1 Transition from two to four-stroke

Consider a stick-slip oscillator model with a pole and polynomial-type
approximation for the characteristic:

$\ddot{X}+F_{f,\epsilon}(v_{r})+X=0,\qquad
F_{f,\epsilon}(v_{r})=\tilde{N}\bigg{(}-\mu_{s}+a_{1}|v_{r}|-a_{3}|v_{r}|^{3}+\frac{\epsilon}{|v_{r}|}\bigg{)},$

where $\tilde{N}$ is the (dimensionless) normal force associated with the
mass. As before, $F_{f,\epsilon}(v_{r})$ coincides with a common choice of
dynamic friction law for $v_{r}\neq 0$ (see, e.g. [14, 15, 32]), and the limit
$\epsilon,|v_{r}|\to 0$ approximates the stick phase. Restricting to
$v_{r}<0$, setting $X=\tilde{N}x$, $y=-v_{r}$, and applying the time
desingularisation

$d\tau=\frac{y}{\tilde{N}}dt,$

we obtain the dynamical system

$\displaystyle\begin{array}[]{lcl}x^{\prime}=\delta(v_{0}-y)y,\\\
y^{\prime}=(x-\mu_{s}+a_{1}y-a_{3}y^{3})y+\epsilon,\end{array}$ (5.6)

where $\delta:=\tilde{N}^{-2}$. System (5.6) exhibits two and four-stroke
relaxation oscillation in different (limiting) regions of
$(\epsilon,\delta)$-parameter space; see Figures 5.1 and 5.2:

(a)

(b)

Figure 5.1: Two-stroke relaxation oscillation in (5.6) for
$(\epsilon,\delta)=(10^{-2},5)$, and
$(v_{0},\mu_{s},a_{1},a_{3})=(2,9,4,10^{-1})$.

(a)

(b)

Figure 5.2: Four-stroke relaxation oscillation in (5.6) for
$(\epsilon,\delta)=(5,10^{-2})$, with all other parameters the same as in
Figure 5.1.

  * R1.

$\epsilon\ll 1$, $\delta=\mathcal{O}(1)$. Two-stroke relaxation oscillations;

  * R2.

$\epsilon=\mathcal{O}(1)$, $\delta\ll 1$. Four-stroke relaxation oscillations.

Hence, there is a transition from two- to four-stroke relaxation oscillations
as one traverses a path in $(\epsilon,\delta)$-parameter space from region R1
to R2. This example shows that a single oscillator can exhibit both two- and
four-stroke relaxation oscillations if there are multiple small perturbation
parameters in the model which allows for different singular limits.

######  Remark 5.1.

In the context of stick-slip oscillations, the case $0<\delta\ll\epsilon\ll 1$
is non-physical: the limit $\epsilon\to 0$ approximates a discontinuity, so
one should have $\epsilon<\delta$ asymptotically.

###  5.2 Timescales and the two/four-stroke distinction

Consider system (5.6) with $0<\epsilon\ll 1$ and $\delta=1$. If we define the
parameters $a_{1}$, $a_{3}$ as in (2.24), this system coincides with the two-
stroke oscillator model (2.31) with polynomial characteristic. Our analysis in
Section 3.1 showed that the equilibrium $p_{0}$ is stable for $v_{0}>v_{m}$,
in which case oscillations are not possible. Physically, large belt speeds
$v_{0}>v_{m}$ mean that the mass cannot ‘stick’: its position stays fixed
while the belt slides underneath it in an equilibrium state known as ‘steady
sliding’ [14, 31, 32]. As one decreases the belt speed, small amplitude
oscillations known as ‘pure-slip’ oscillations appear for $v_{0}<v_{m}$. These
oscillations are not of stick-slip type, and exist only in a narrow parameter
regime $(v_{ss},v_{m})$; see e.g. [48]. Stick-slip oscillations occur only
once the belt speed is decreased below $v_{ss}$; see Figures 5.3 and 5.4. The
value of $v_{m}$ can be identified as a supercritical Andronov-Hopf
bifurcation, but the value $v_{ss}$ is typically harder to identify and known
only for a few specific cases. For example, the authors in [14] show that
under the assumption $\mu_{s}-\mu_{m}\ll 1$,

$v_{ss}=\sqrt{\frac{4}{5}}v_{m},$ (5.7)

for the stick-slip oscillator with characteristic (2.24).

Figure 5.3: Behaviour for different belt speeds $v_{0}$. Region $I$: stick-
slip oscillation; region $II$: pure-slip oscillation; region $III$: steady
sliding (no oscillation).

(a)

(b)

(c)

Figure 5.4: System (5.6) with $\delta=1$ and
$(\mu_{s},a_{1},a_{3},\epsilon)=(1,3/4,1/4,10^{-3})$. Shown are representative
dynamics for different belt speeds: (a) steady sliding for $v_{0}=1.1$; (b)
pure-slip oscillation for $v_{0}=0.96$; (c) stick-slip oscillation for
$v_{0}=0.86$.

(a)

(b)

(c)

Figure 5.5: System (5.6) with
$\delta={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}10^{-2}}$
and $(\mu_{s},a_{1},a_{3},\epsilon)=(1,3/4,1/4,10^{-3})$. Shown are
representative dynamics for different belt speeds: (a) steady sliding for
$v_{0}=1.1$; (b) pure-slip oscillation for $v_{0}=0.997166$; (c) stick-slip
oscillation for $v_{0}=0.997165$.

Now consider the case $0<\epsilon\ll\delta\ll 1$, which corresponds physically
to the case of a large normal force $(\delta=\tilde{N}^{-2})$. Figure 5.5
shows that by decreasing $\delta$ solutions begin to track the upper branch of
the characteristic in phase space. The presence of two singular perturbation
parameters indicates the presence of three timescales, however, which
distinguishes the relaxation oscillations observed for $0<v_{0}<v_{ss}$ from
vdP-type relaxation oscillations.

The additional timescale allows for an interpretation of the transition from
steady sliding to stick-slip using standard GSPT. The layer problem obtained
in the limit $\epsilon\to 0$ has Jacobian at $p_{0}$ given by

$J=\begin{pmatrix}{0}&{-\delta v_{0}}\\\
{v_{0}}&{(a_{1}-3a_{3}v_{0}^{2})v_{0}}\end{pmatrix},$

which has eigenvalues $\pm i\delta v_{0}^{2}$ at the Andronov-Hopf value
$v_{0}=\sqrt{a_{1}/3a_{3}}$, i.e. there is a singular Andronov-Hopf
bifurcation [25, 49, 50] as $\delta\to 0$. Moreover, one observes a rapid
decrease in the width of the interval $(v_{ss},v_{m})$ as $\delta\to 0$. This
can be seen by comparing Figure 5.5, which shows the dynamics for different
belt speeds with $\delta=10^{-2}$, with Figure 5.4. Pure-slip and stick-slip
cases in Figure 5.5 are separated by a change in belt speed $v_{0}$ which is
$\mathcal{O}(10^{-6})$. Both observations point to the existence of a _canard
explosion_ in system (5.6) with $0<\epsilon\ll\delta\ll 1$; a dynamic
phenomena characterised by a smooth transition from small-amplitude
oscillations born in a singular Andronov-Hopf bifurcation to large-amplitude
relaxation oscillation over an exponentially small interval in parameter
space. In particular, determination of the parameter value $v_{ss}$ separating
pure-slip and stick-slip solutions reduces to the problem of locating the
solution which tracks the inner branch of the characteristic for the longest
time.333In the language of standard GSPT, the value $v_{ss}$ corresponds to
the location of the maximal canard, which separates canard cycles with and
without head. A complete study of the dynamics associated with system (5.6)
constitutes future/ongoing work.

##  6 Conclusion and future work

Two-stroke oscillation is an important phenomenon occurring across nature, as
well as in a host of engineering problems. We have illustrated that two-stroke
oscillators can be well described via an underlying singularly perturbed
problem featuring two-stroke relaxation oscillation. Geometric singular
perturbation theory can be adapted for the study of two-stroke relaxation
oscillations, as we showed in the context of a number of applications deriving
from the study of electronic oscillators, and mechanical oscillators with
friction. Our study showed that relaxation oscillations can occur under
minimal conditions which differ from those associated with the (four-stroke)
vdP-type relaxation oscillations. This observation is made also in [17, 18,
20], where two-stroke relaxation oscillations are observed in singularly
perturbed problems in standard form (1.12). The models proposed in these works
describe two-stroke relaxation oscillations as perturbations of singular
relaxation cycles containing segments at infinity. Our approach shows that by
relaxing the requirement that the system is globally expressible in the
standard form (1.12) and hence allowing for non-zero curvature of the layer
flow, two-stroke relaxation cycles can be described as perturbations of
singular relaxation cycles residing in a compact region of phase space.
Moreover, our approach is applicable to the analysis of systems like (1.21),
which cannot be understood as a standard form problem with a return mechanism
at infinity.

The manuscript also raises a number of interesting questions in both theory
and application. The remaining discussion focuses on a small number of these
in turn. Each constitutes either current or future work.

Onset of stick-slip by canard explosion in a three-timescale problem. As
outlined in Section 5.2, the transition from ‘steady sliding’ to ‘stick-slip’
oscillation in the model (5.6) with $0<\epsilon\ll\delta\ll 1$ appears to
occur via a canard explosion, with the transition from small-amplitude sliding
oscillations to large-amplitude stick-slip oscillations occurring at the
transition from canard cycles without head to canard cycles with head. The
observed dynamics is, more generally, a feature of all stick-slip oscillators
with a characteristic $\mu(v_{r})$ that has a turning point bounded away from
$v_{r}=0$. We are interested in a complete study of the canard explosions
exhibited by these systems, and in particular, with the transition from (two-
timescale) canard cycles without head to (three-timescale) canard cycles with
head; see Figure 5.5.

More general friction forces and exponential nonlinearities. A natural
progression from this work is to study two-stroke oscillators for which the
friction has ‘rate-and-state’ dependence. Such a generalised friction force
$F_{f}(x,v_{r},t)$ models rate-and-state dependence in a wide range of
physical phenomena and engineering applications [10, 48]. From a modelling
point of view, incorporating rate dependence in particular leads one into the
realm of non-autonomous dynamical systems (‘open’ models). If the friction
force changes only slowly in time then singular perturbation techniques can
still provide answers to understand the underlying non-autonomous dynamics;
see e.g. [51].

In addition to the difficulties associated with non-autonomous dynamical
systems, many of the existing models involving exponential or logarithmic type
(friction) terms lead, in the context of singularly perturbed problems, to
more general problems surrounding the difficulties associated with the
breakdown of normal hyperbolicity at greater than algebraic rates. Significant
progress has been made in this area by in [17, 52, 53] by adapting the blow-up
method in order to deal with such difficulties. These techniques are also
applicable in the study of certain two-stroke oscillations with exponential
nonlinearities occurring in the aformentioned electrical oscillators
originally introduced by Hester [29] which is ongoing work.

Canards and bifurcations two-stroke relaxation oscillators. In the case of
vdP-type (four-stroke) relaxation oscillation, the mechanism responsible for
the onset of relaxation oscillation is a canard explosion, i.e. a rapid
transition from small to large amplitude relaxation oscillations under an
exponentially small variation in parameter space. Canard explosion in the vdP
oscillator with constant forcing is well understood [25], and known to
correlate with the passage of an equilibrium over a fold under additional
parameter variation. Although the canard explosion itself is a global
phenomenon, it is driven by the local dynamics, in particular the occurrence
of a (singular) Andronov-Hopf bifurcation, and the presence of small (and
large) canard cycles.

Given the local equivalence between general singularly perturbed problems
(3.1) and standard form problems (1.12) discussed in Section 3.3, then, a
natural question is the following: what is the global effect of introducing
analogous local dynamics in the case of a two-stroke relaxation oscillator?
Canard explosions observed in general systems (3.1) capable of generating two-
stroke relaxation oscillation may differ qualitatively from the vdP-type
canard explosion in a number of ways (see, e.g., [8]): among other
distinguishing features, canard cycles can have arbitrarily large amplitudes
which exceed that of the relaxation cycle significantly, and the onset of
stable two-stroke relaxation oscillation is possible only if the associated
(singular) Andronov-Hopf bifurcation is subcritical. In addition to the study
of canard explosion under ‘minimal’ conditions, we currently consider the
possibilities for global dynamics associated with two-stroke relaxation
oscillators with up to two additional equilibria on the slow manifold. The
occurrence of a local (singular) Bogdanov-Takens bifurcation (see [50, 54]) in
this scenario leads to a wide variety of possibilities for global bifurcations
and dynamics.

Regularised boundary equilibrium bifurcation in stick-slip. The bifurcations
discussed in Section 5.2 are known to occur in general stick-slip problems for
belt speeds $v_{0}$ bounded uniformly away from zero. Bifurcations occurring
in a neighbourhood of the regularised region $v_{r}=0$ are necessarily more
complex, and are associated with boundary equlibrium bifurcations in the
piecewise-smooth (PWS) literature [33, 34, 55], which occur when an
equilibrium collides with a discontinuity surface under variation of a
parameter. In an ongoing work, we consider the case of a (regularised)
‘boundary focus’ bifurcation which occurs in stick-slip problems as $v_{0}\to
0^{+}$. In the PWS case, it is known that the amplitude of stick-slip cycles
shrinks to zero and terminates at the point of collision. In the regularised
analogue, we identify an Andronov-Hopf bifurcation as the bifurcation
responsible for the termination of oscillations for $0<\epsilon\ll 1$, and
describe the growth of these cycles to order one stick-slip type cycles
observed for belt speeds $v_{0}$ bounded uniformly away from zero.

The general framework for $\mathbb{R}^{n}$. The framework for GSPT applied to
general systems (3.1) developed in this work extends beyond the study of
planar two-stroke oscillations. In particular, the framework presented in this
work generalises to higher dimensional singular perturbation problems in the
general form

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

where $N(z)$ denotes a matrix of dimension $n\times(n-k)$ and $f(z)$ is a
$(n-k)$-dimensional (vector-valued) function, and the existence of an
$k$-dimensional  regularly embedded submanifold
$S=\\{z\in\mathbb{R}^{n}|f(z)=0\\}\subset\mathbb{R}^{n}$ is assumed, $1\leq
k<n$. This work is developed in [8].

## Acknowledgments

This work was supported by the Australian Research Council DP180103022 grant.
The authors would also like to thank Peter Szmolyan, Ilona Kosiuk, Kristian
Uldall Kristiansen and the referees for helpful comments on an early version
of the manuscript. In particular, we would like to acknowledge the work of
Peter Szmolyan and Ilona Kosiuk as motivation for this work.

## References

  * [1] van der Pol B.  A theory of the amplitude of free and forced triode vibrations.  Radio Review. 1920;1(1920):701–710.
  * [2] van der Pol B.  On “relaxation-oscillations”.  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science. 1926;2(11):978–992.
  * [3] Kuehn C.  Multiple time scale dynamics. vol. 191.  Springer; 2015.
  * [4] Huber A, Szmolyan P.  Geometric singular perturbation analysis of the Yamada model.  SIAM Journal on Applied Dynamical Systems. 2005;4(3):607–648.
  * [5] Gucwa I, Szmolyan P.  Geometric singular perturbation analysis of an autocatalator model.  Discrete and Continuous Dynamical Systems. 2009;2(4):783–806.
  * [6] Kuehn C, Szmolyan P.  Multiscale geometry of the Olsen model and non-classical relaxation oscillations.  Journal of Nonlinear Science. 2015;25(3):583–629.
  * [7] Kosiuk I, Szmolyan P.  Geometric analysis of the Goldbeter minimal model for the embryonic cell cycle.  Journal of mathematical biology. 2016;72(5):1337–1368.
  * [8] Wechselberger M.  Geometric singular perturbation theory beyond the standard form.  To appear in Frontiers in Applied Dynamics: Reviews and Tutorials. 2019;.
  * [9] Le Corbeiller P.  Two-stroke oscillators.  IRE Transactions on Circuit Theory. 1960;7(4):387–398.
  * [10] Berger E.  Friction modeling for dynamic system simulation.  Applied Mechanics Reviews. 2002;55(6):535–577.
  * [11] Bossolini E, Brøns M, Kristiansen KU.  Canards in stiction: On solutions of a friction oscillator by regularization.  SIAM Journal on Applied Dynamical Systems. 2017;16(4):2233–2258.
  * [12] Pomeau Y, Berre ML.  Critical speed-up vs critical slow-down: a new kind of relaxation oscillation with application to stick-slip phenomena.  arXiv preprint arXiv:11073331. 2011;.
  * [13] Popp K, Stelter P.  Stick-slip vibrations and chaos.  Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. 1990;332(1624):89–105.
  * [14] Thomsen JJ, Fidlin A.  Analytical approximations for stick–slip vibration amplitudes.  International Journal of Non-Linear Mechanics. 2003;38(3):389–403.
  * [15] Won HI, Chung J.  Stick–slip vibration of an oscillator with damping.  Nonlinear Dynamics. 2016;86(1):257–267.
  * [16] Puu T, Sushko I.  Business cycle dynamics: Models and tools.  Springer-Verlag; 2006.
  * [17] Kristiansen KU.  Blowup for flat slow manifolds.  Nonlinearity. 2017;30(5):2138.
  * [18] Rankin J, Desroches M, Krauskopf B, Lowenberg M.  Canard cycles in aircraft ground dynamics.  Nonlinear Dynamics. 2011;66(4):681–688.
  * [19] Omta AW, Kooi BW, van Voorn GA, Rickaby RE, Follows MJ.  Inherent characteristics of sawtooth cycles can explain different glacial periodicities.  Climate dynamics. 2016;46(1-2):557–569.
  * [20] Brøns M.  Relaxation oscillations and canards in a nonlinear model of discontinuous plastic deformation in metals at very low temperatures.  In: Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. vol. 461. The Royal Society; 2005. p. 2289–2302.
  * [21] Goldbeter A.  Biochemical oscillations and cellular rhythms: the molecular bases of periodic and chaotic behaviour.  Cambridge university press; 1997.
  * [22] Guckenheimer J, Holmes P.  Nonlinear oscillations, dynamical systems, and bifurcation of vector fields.  Springer-Verlag, Berlin Heidelberg New York; 1983.
  * [23] Kuznetsov YA.  Elements of applied bifurcation theory. vol. 112.  Springer Science & Business Media; 2013.
  * [24] Goeke A, Walcher S.  A constructive approach to quasi-steady state reductions.  Journal of mathematical chemistry. 2014;52(10):2596–2626.
  * [25] Krupa M, Szmolyan P.  Relaxation oscillation and canard explosion.  Journal of Differential Equations. 2001;174(2):312–368.
  * [26] Fenichel N.  Geometric singular perturbation theory for ordinary differential equations.  Journal of Differential Equations. 1979;31(1):53–98.
  * [27] Kaleda PI.  Singular systems on the plane and in space.  Journal of Mathematical Sciences. 2011;179(4):475–490.
  * [28] De Maesschalck P, Dumortier F, Roussarie R.  Cyclicity of common slow–fast cycles.  Indagationes Mathematicae. 2011;22(3-4):165–206.
  * [29] Hester D.  The nonlinear theory of a class of transistor oscillators.  IEEE Transactions on Circuit Theory. 1968;15(2):111–117.
  * [30] Ibrahim R.  Friction-induced vibration, chatter, squeal, and chaos—part II: dynamics and modeling.  Applied Mechanics Reviews. 1994;47(7):227–253.
  * [31] Panovko YG, Gubanova II.  Stability and oscillations of elastic systems, paradoxes, fallacies and new concepts.  Consultants Bureau; 1965.
  * [32] Chen GS.  Handbook of friction-vibration interactions.  Elsevier; 2014.
  * [33] Bernardo M, Budd C, Champneys AR, Kowalczyk P.  Piecewise-smooth dynamical systems: theory and applications. vol. 163\.  Springer Science & Business Media; 2008.
  * [34] Jeffrey MR, Jeffrey MR, Chernyk.  Hidden Dynamics.  Springer; 2018.
  * [35] Bonet Revés C, Martínez-Seara Alonso MT.  Regularization of sliding global bifurcations derived from the local fold singularity of Filippov systems.  Discrete and continuous dynamical systems Series A. 2016;36(7):3545–3601.
  * [36] Kristiansen KU, Hogan SJ.  Regularizations of two-fold bifurcations in planar piecewise smooth systems using blowup.  SIAM Journal on Applied Dynamical Systems. 2015;14(4):1731–1786.
  * [37] Kristiansen KU, Hogan SJ.  On the use of blowup to study regularizations of singularities of piecewise smooth dynamical systems in $\mathbb{R}^{3}$.  SIAM Journal on Applied Dynamical Systems. 2015;14(1):382–422.
  * [38] Kristiansen KU, Hogan SJ.  Resolution of the Piecewise Smooth Visible–Invisible Two-Fold Singularity in $\mathbb{R}^{3}$ Using Regularization and Blowup.  Journal of Nonlinear Science. 2019;29(2):723–787.
  * [39] Kristiansen KU.  The regularized visible fold revisited.  arXiv preprint arXiv:190806781. 2019;.
  * [40] Jones CK.  Geometric singular perturbation theory.  In: Dynamical systems. Springer; 1995. p. 44–118.
  * [41] Fenichel N.  Asymptotic stability with rate conditions.  Indiana University Mathematics Journal. 1974;23(12):1109–1137.
  * [42] Fenichel N.  Asymptotic stability with rate conditions, II.  Indiana University Mathematics Journal. 1977;26(1):81–93.
  * [43] Tikhonov AN.  Systems of differential equations containing small parameters in the derivatives.  Matematicheskii sbornik. 1952;73(3):575–586.
  * [44] Krupa M, Szmolyan P.  Extending geometric singular perturbation theory to nonhyperbolic points—fold and canard points in two dimensions.  SIAM journal on mathematical analysis. 2001;33(2):286–314.
  * [45] Arnold V.  Ordinary Differential Equations.  Translated, edited by Richard Silverman, editors. MIT Press, Cambridge; 1973.
  * [46] Chicone C.  Ordinary Differential Equations with Applications; 1999.
  * [47] Hirsch MW, Smale S, Devaney RL.  Differential equations, dynamical systems, and an introduction to chaos.  Academic press; 2012.
  * [48] Dankowicz H, Nordmark AB.  On the origin and bifurcations of stick-slip oscillations.  Physica D: Nonlinear Phenomena. 2000;136(3-4):280–302.
  * [49] Dumortier F, Roussarie R.  Canard cycles and center manifolds. vol. 577.  American Mathematical Soc.; 1996.
  * [50] De Maesschalck P, Wechselberger M.  Neural excitability and singular bifurcations.  The Journal of Mathematical Neuroscience (JMN). 2015;5(1):16.
  * [51] Wechselberger M, Mitry J, Rinzel J.  Canard theory and excitability.  In: Nonautonomous dynamical systems in the life sciences. Springer; 2013\. p. 89–132.
  * [52] Bossolini E, Brøns M, Kristiansen KU.  Singular limit analysis of a model for earthquake faulting.  Nonlinearity. 2017;30(7):2805.
  * [53] Kristiansen KU.  A new type of relaxation oscillation in a model with rate-and-state friction.  arXiv preprint arXiv:190312232. 2019;.
  * [54] De Maesschalck P, Dumortier F.  Slow-fast Bogdanov-Takens bifurcations.  Journal of differential equations. 2011;250(2):1000–1025.
  * [55] Kuznetsov YA, Rinaldi S, Gragnani A.  One-parameter bifurcations in planar Filippov systems.  International Journal of Bifurcation and chaos. 2003;13(08):2157–2188.
  * [56] Mishchenko E, Rozov N Kh.  Differential equations with small parameters and relaxation oscillations.  Nauka, Moscow. (transl.: Plenum Press, New York, 1980); 1975.

##  Appendix A Local transformation of system (3.1) to standard form

Due to the existence of a local fibration in a sufficiently narrow tubular
neighbourhood $\mathcal{B}$ of $S$, the study of general slow-fast systems
(3.1) is related to the study of standard form problems (3.3) by local
(topological) equivalence. This is emphasised in the following theoretical
result, which provides the means for locally rectifying the fibers near $S$.

######  Lemma A.1.

Given system (3.1) with critical manifold $S$. Then there exists a local
coordinate transformation such that system (3.1) takes the standard form

$\begin{pmatrix}{x^{\prime}}\\\
{v^{\prime}}\end{pmatrix}=\begin{pmatrix}{1}\\\
{0}\end{pmatrix}\tilde{f}(x,v)+\epsilon\begin{pmatrix}{g_{1}(x,v,\epsilon)}\\\
{g_{2}(x,v,\epsilon)}\end{pmatrix}$ (A.1)

in a local tubular neighbourhood $\mathcal{B}$ of $S$.

###### Proof.

Consider system (3.1) with coordinates $z=(x,y)^{T}$. Assume without loss of
generality that fast fibers in a local tubular neighbourhood $\mathcal{B}$ of
$S$ can be described as level sets $L(x,y)=c\in I\subset\mathbb{R}$ for some
smooth real-valued function $L(x,y)$ with $D_{y}L(x,y)|_{\mathcal{B}}\neq 0$,
i.e. each fast fiber can be locally written as a graph $y=Y^{c}(x)$. We define
a new local coordinate $v=L(x,y)$ which has a locally well-defined inverse
$y=\tilde{M}(x,v)$, and obtain the standard form system

$\begin{pmatrix}{x^{\prime}}\\\
{v^{\prime}}\end{pmatrix}=\begin{pmatrix}{1}\\\
{0}\end{pmatrix}N_{1}(x,\tilde{M}(x,v))f(x,\tilde{M}(x,v))+\epsilon\begin{pmatrix}{G_{1}(x,\tilde{M}(x,v),\epsilon)}\\\
{\langle
DL(x,\tilde{M}(x,v)),G(x,\tilde{M}(x,v),\epsilon)\rangle}\end{pmatrix},$

where we have used the fact that $\langle DL,N\rangle|_{\mathcal{B}}=0$ is the
defining condition for a local fast foliation. ∎

##  Appendix B Proof of Theorem 3.6

Without loss of generality, assume that $F=(x_{F},y_{F})=(0,0)$, and
$D_{y}f|_{F}>0$ (one can always ensure this by adjusting $N$ if necessary).
Making the preliminary transformation $u=f(x,y)$ and denoting the inverse of
the transformation by $y=M(x,u)$ as in proof of Lemma 3.3, we obtain the
system

$\begin{split}\begin{pmatrix}{x^{\prime}}\\\
{u^{\prime}}\end{pmatrix}&:=\tilde{N}(x,u)u+\epsilon\tilde{G}(x,u,\epsilon)\\\
&=\begin{pmatrix}{N_{1}(x,M(x,u))}\\\
{\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f(x,M(x,u),N(x,M(x,u)){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0})}\rangle}\end{pmatrix}u+\epsilon\begin{pmatrix}{G_{1}(x,M(x,u),\epsilon)}\\\
{\langle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla}f(x,M(x,u)),G(x,M(x,u),\epsilon)\rangle}\end{pmatrix},\end{split}$
(B.1)

and the sections $\Sigma_{in/out}$ defined in (3.27) with
$(x_{F},y_{F})=(0,0)$ are mapped to

$\tilde{\Sigma}_{in}=\big{\\{}(-\rho,M(-\rho,u))\big{|}u\in\tilde{J}\big{\\}},\qquad\tilde{\Sigma}_{in}=\big{\\{}(\rho,M(\rho,u))\big{|}u\in\tilde{J}\big{\\}},$

which can in turn be rewritten as

$\tilde{\Sigma}_{in}=\big{\\{}(-\rho,u)\big{|}u\in\hat{J}_{in}\big{\\}},\qquad\tilde{\Sigma}_{in}=\big{\\{}(\rho,u)\big{|}u\in\hat{J}_{out}\big{\\}},$
(B.2)

since $M(\rho,u)=\mathcal{O}(u)$. One can verify that the presence of a jump-
off point at $(x,y)=(0,0)$ implies the presence of a jump-off point at
$(x,u)=(0,0)$ in the system (B.1).444Since $D_{y}f|_{F}>0$, the coordinate
transformation preserves orientation. We consider the extended system obtained
from (B.1) by adding the trivial equation $\epsilon^{\prime}=0$ and expanding
about $(x,u,\epsilon)=(0,0,0)$,

$\begin{pmatrix}{x^{\prime}}\\\ {u^{\prime}}\\\
{\epsilon^{\prime}}\end{pmatrix}=\begin{pmatrix}{a_{0}+a_{1}x+a_{2}u+\cdots}\\\
{b_{1}x+b_{2}u+\cdots}\\\
{0}\end{pmatrix}u+\epsilon\begin{pmatrix}{c_{0}+c_{1}x+c_{2}u+c_{3}\epsilon+\cdots}\\\
{d_{0}+d_{1}x+d_{2}u+d_{3}\epsilon+\cdots}\\\ {0}\end{pmatrix},$ (B.3)

where $a_{0}\neq 0$, $b_{1}\neq 0$, $d_{0}\neq 0$ by the definition of regular
contact (Definition 3.6). We assume in the following that $a_{0}$, $b_{0}$,
$d_{0}$ are positive; the proof is analogous for different choices provided
the relative orientations are consistent with $F$ being a jump-off point. We
now define the blow-up $\Phi:S^{2}\times[0,\rho]\to\mathbb{R}^{3}$ by the
mapping

$(x,u,\epsilon)=\big{(}r\bar{x},r^{2}\bar{u},r^{3}\bar{\epsilon}\big{)},\qquad\big{(}\bar{x},\bar{u},\bar{\epsilon},r\big{)}\in
S^{2}\times[0,\rho],$ (B.4)

where in particular, since ${\epsilon}\geq 0$, we need only consider the
dynamics on and near the hemisphere

$(\bar{x},\bar{u},\bar{\epsilon})\in
S^{2}_{+}=\\{(\bar{x},\bar{u},\bar{\epsilon})|\bar{x}^{2}+\bar{u}^{2}+\bar{\epsilon}^{2}=1,\bar{\epsilon}\geq
0\\}.$

The dynamics are studied in different coordinate charts. We define entry and
exit charts $K_{1}:\bar{x}=-1$ respectively $K_{3}:\bar{x}=1$, and a family
rescaling chart $K_{2}:\bar{\epsilon}=1$. In line with conventions, we denote
the image of an object $\gamma$ under the blow-up map (B.4) by $\bar{\gamma}$,
and its image in a specific chart $K_{j}$ by $\gamma_{j}$, for
$j\in\\{1,2,3\\}$. Charts $K_{1}$ and $K_{3}$ allow one to describe the
extension of the (extended) slow manifolds $\bar{S^{a}_{\epsilon}}$ and
$\bar{S^{r}_{\epsilon}}$ respectively into the neighbourhood of $S_{+}^{2}$,
while the flow on (and near) the interior of $S_{+}^{2}$ is understood in
chart $K_{2}$.

In order to put the technicalities into context, we briefly summarise the key
findings in charts $K_{1}$ and $K_{3}$ before considering a more detailed
analysis. The main dynamical features are sketched in Figure B.1, which shows
a birds-eye perspective of the dynamics in the blow-up, restricted to the
hemisphere $S_{+}^{2}$. We refer to Figure B.1 throughout the proof; for now
it suffices to note the existence of four singularities on the equator
$S^{1}$. Those denoted $p_{a}$ and $p_{r}$ are partially hyperbolic and
attracting/repelling along $S^{1}$ respectively, and $q_{in}$ and $q_{out}$
are hyperbolic singularities corresponding to the intersection with the
(extended) critical fiber $\mathcal{\bar{F}}$.

Figure B.1: Dynamics on $S^{2,+}$.

We now present an analysis of the dynamics in charts $K_{i}$,
$i\in\\{1,2,3\\}$ in turn. Most of the analysis is conceptually analogous to
the work in [44], in which the authors study the dynamics of a regular fold
point in standard form problems (3.3). We cite corresponding and relevant
results contained in this work as the analysis proceeds.

Chart $K_{2}$. After an additional desingularisation
$dt=r_{2}^{-1}d\bar{t}_{2}$, we obtain the following in the family rescaling
chart $K_{2}$:

$\displaystyle\begin{array}[]{lcl}x_{2}^{\prime}=a_{0}u_{2}+\mathcal{O}(r_{2}),\\\
u_{2}^{\prime}=d_{0}+b_{1}x_{2}u_{2}+\mathcal{O}(r_{2}),\\\
r_{2}^{\prime}=0,\end{array}$ (B.8)

where the subscript notation indicates the use of chart-specific coordinates

$(x,u,\epsilon)=(r_{2}x_{2},r_{2}^{2}u_{2},r_{2}^{3}).$

System (B.8) is a regular perturbation problem with perturbation parameter
$r_{2}$, and so the dynamics are governed by the limiting system

$\displaystyle\begin{array}[]{lcl}x_{2}^{\prime}=a_{0}u_{2},\\\
u_{2}^{\prime}=d_{0}+b_{1}x_{2}u_{2}.\end{array}$ (B.11)

In the following we define the segments

$\Sigma_{2}^{in}=\big{\\{}(x_{2},u_{2},r_{2})\big{|}x=-\delta^{-1/3},r_{2}\in[0,\delta^{1/3}\rho]\big{\\}},\qquad\Sigma_{2}^{out}=\big{\\{}(x_{2},u_{2},r_{2})\big{|}x=\delta^{-1/3},r_{2}\in[0,\delta^{1/3}\rho]\big{\\}},$

(see Figure B.1) and denote the transition map defined by the flow of (B.8) by
$\Pi_{2}:\Sigma_{2}^{in}\to\Sigma_{2}^{out}$.

######  Proposition B.1.

(cf. Proposition 2.3 and 2.4 in [44]). There exists a solution $\gamma_{c,2}$
for (B.11) such that for the initial condition
$q\in\Sigma_{2}^{in}\cap\gamma_{c,2}$, we have the following:

  * (i)

$\Pi_{2}(q)=\bigg{(}\delta^{-1/3},\bigg{(}\frac{b_{1}}{2a_{0}}\bigg{)}\delta^{-2/3}+\bigg{(}\frac{2d_{0}^{2}}{a_{0}b_{1}}\bigg{)}^{1/3}\Omega_{0}-\bigg{(}\frac{2d_{0}}{b_{1}}\bigg{)}\delta^{1/3}+\mathcal{O}(\delta),0\bigg{)}^{T},$

  * (ii)

A neighbourhood of $q$ is mapped diffeomorphically onto a neighbourhood of
$\Pi_{2}(q)$.

###### Proof.

The key observation is that (B.11) is a Ricatti equation with special solution
$\gamma_{c,2}$ (see Figure B.1). In particular, $\gamma_{c,2}$ can be
parameterised as $(x_{2},\zeta(x_{2}))$, where $\zeta(x_{2})$ has known
asymptotic behaviour described by

$\zeta(x_{2})=\bigg{(}\frac{2d_{0}^{2}}{a_{0}b_{1}}\bigg{)}^{1/3}\Omega_{0}+\bigg{(}\frac{b_{1}}{2a_{0}}\bigg{)}x_{2}^{2}-\bigg{(}\frac{2d_{0}}{b_{1}}\bigg{)}\frac{1}{x_{2}}+\mathcal{O}\bigg{(}\frac{1}{x_{2}^{3}}\bigg{)},\qquad
x_{2}\to\infty,$

where $\Omega_{0}$ is a known positive constant,555$\Omega_{0}$ is the
smallest positive solution to $J_{-1/3}(2z^{3/2}/3)+J_{1/3}(2z^{3/2}/3)$,
where $J_{-1/3}$ and $J_{1/3}$ are Bessel functions of the first kind. See
Remark 2.4. in [44] and also [56]. and

$\zeta(x_{2})=-\bigg{(}\frac{d_{0}}{b_{1}}\bigg{)}\frac{1}{x_{2}}+\mathcal{O}\bigg{(}\frac{1}{x_{2}^{4}}\bigg{)},\qquad
x_{2}\to-\infty.$ (B.12)

The result (i) follows, and (ii) follows by standard results from regular
perturbation theory. ∎

Chart $K_{1}$. After an additional desingularisation
$dt=r_{1}^{-1}d\bar{t}_{1}$, we obtain the following in the entry chart
$K_{1}$:

$\displaystyle\begin{array}[]{lcl}u_{1}^{\prime}=-b_{1}u_{1}+d_{0}\epsilon_{1}+2u_{1}\theta_{1}(u_{1},\epsilon_{1},r_{1})+r_{1}((b_{3}+b_{2}u_{1})u_{1}-d_{1}\epsilon_{1})+\mathcal{O}(r_{1}^{2}),\\\
\epsilon_{1}^{\prime}=3\epsilon_{1}\theta_{1}(u_{1},\epsilon_{1},r_{1}),\\\
r_{1}^{\prime}=-r_{1}\theta(u_{1},\epsilon_{1},r_{1}),\end{array}$ (B.16)

where the subscript notation indicates the use of chart-specific coordinates

$(x,u,\epsilon)=(-r_{1},r_{1}^{2}u_{1},r_{1}^{3}\epsilon_{1}),$

and

$\theta_{1}(u_{1},\epsilon_{1},r_{1})=a_{0}u_{1}+r_{1}(-a_{1}u_{1}+c_{0}\epsilon_{1})+\mathcal{O}(r_{1}^{2}).$

The portion of the equator $S^{1}$ visible in chart $K_{1}$ can be identified
with the invariant line $\\{(u_{1},0,0)|u_{1}\in\mathbb{R}\\}$, along which we
identify the partially hyperbolic singularity $p_{a}=(0,0,0)$, and hyperbolic
saddle $q_{in}=(b_{1}/2a_{0},0,0)$ shown in Figure B.1. Our main task here is
to understand the dynamics near $p_{a}$, so we restrict our analysis to the
set

$D_{1}:=\big{\\{}(u_{1},\epsilon_{1},r_{1})\big{|}(\epsilon_{1},r_{1})\in[0,\delta]\times[0,\rho]\big{\\}}.$

The situation is sketched in Figure B.2, which shows all the relevant objects
in the analysis near $p_{a}$.

Figure B.2: Setup for the local analysis near $p_{a}$.

######  Proposition B.2.

(cf. Proposition 2.6 in [44]). Given $\delta,\rho$ sufficiently small, we have
the following:

  * (i)

There exists a locally invariant, two-dimensional attracting centre manifold
$W_{a,1}^{loc}$, which is tangent to

$E^{c}=span\\{(0,0,1)^{T},(d_{0},b_{1},0)^{T}\\}$

at $p_{a}$ and given by a graph $u_{1}=g(\epsilon_{1},r_{1})$.

  * (ii)

$W_{a,1}^{loc}$ contains two locally invariant, one-dimensional centre
manifolds as restrictions

$N^{a}_{1}=W_{a,1}^{loc}\big{|}_{r_{1}=0}\qquad\text{and}\qquad
S_{1}^{a}=W_{a,1}^{loc}\big{|}_{\epsilon_{1}=0}.$

The manifold $N^{a}_{1}$ is unique in $D_{1}$ and coincides with the image of
the special trajectory $\gamma_{c}$ in chart $K_{1}$, denoted $\gamma_{c,1}$.
The manifold $S^{a}_{1}$ coincides with image of the attracting (extended)
critical manifold $\bar{S^{a}}$ in chart $K_{1}$.

  * (iii)

There exists a stable invariant foliatation $\mathcal{F}^{s}$ with base
$W_{a,1}^{loc}$ and one-dimensional fibers. Moreover, we have that $\forall
c>-b_{1}$, $\exists\delta,\rho>0$ such that the contraction rate along
$\mathcal{F}^{s}$ during the time interval $[0,T]$ is greater than $e^{cT}$.

###### Proof.

We prove statements (i) and (ii). Statement (iii) follows from centre manifold
theory (see, e.g, [22]). The Jacobian at $p_{a}$ is given by

$J_{p_{a}}=\begin{pmatrix}{-b_{1}}&{d_{0}}&{0}\\\ {0}&{0}&{0}\\\
{0}&{0}&{0}\end{pmatrix},$

which has eigenvalues $-b_{1}<0$, $0$ and $0$. Centre manifold theory implies
the existence of a locally invariant two-dimensional attracting centre
manifold $W^{loc}_{a,1}$ given as a graph $u_{1}=g(\epsilon_{1},r_{1})$, which
is tangent to the centre eigenspace
$E^{c}=span\\{(d_{0},b_{1},0)^{T},(0,0,1)^{T}\\}$ at $p_{a}$. Making a power
series ansatz for $g(\epsilon_{1},r_{1})$ and matching terms gives a local
expression for the manifold:

$W^{loc}_{a,1}=\bigg{\\{}(u_{1},\epsilon_{1},r_{1})\bigg{|}u_{1}=g(\epsilon_{1},r_{1})=\frac{d_{0}}{b_{1}}\epsilon_{1}+\bigg{(}\frac{b_{3}d_{0}-b_{1}d_{1}}{b_{1}^{2}}\bigg{)}r_{1}\epsilon_{1}-\frac{a_{0}d_{0}^{2}}{b_{1}^{3}}\epsilon_{1}^{2}+\mathcal{O}(3)\bigg{\\}},$

where
$\mathcal{O}(3):=\mathcal{O}(\epsilon_{1}^{3},\epsilon_{1}^{2}r_{1},\epsilon_{1}r_{1}^{2},r_{1}^{3})$.
Dynamics on $W_{a,1}^{loc}$ are determined by restricting (B.16):

$\displaystyle\begin{array}[]{lcl}u_{1}^{\prime}=\frac{3a_{0}d_{0}^{2}}{b_{1}^{2}}\epsilon_{1}^{2}+\mathcal{O}(3),\\\
\epsilon_{1}^{\prime}=\frac{3a_{0}d_{0}}{b_{1}}\epsilon_{1}^{2}+\mathcal{O}(3),\\\
r_{1}^{\prime}=-\frac{a_{0}d_{0}}{b_{1}}r_{1}\epsilon_{1}+\mathcal{O}(3).\\\
\end{array}$ (B.20)

Restricting to the invariant plane $\\{r_{1}=0\\}$ gives

$N^{a}_{1}=\bigg{\\{}(u_{1},\epsilon_{1},0)\bigg{|}u_{1}=g(\epsilon_{1},0)=\frac{d_{0}}{b_{1}}\epsilon_{1}-\frac{a_{0}d_{0}^{2}}{b_{1}^{3}}\epsilon_{1}^{2}+\mathcal{O}(\epsilon_{1}^{3})\bigg{\\}},$

and restricting to the invariant plane $\\{\epsilon_{1}=0\\}$ gives
$S^{a}_{1}=\\{(u_{1},0,r_{1})|u_{1}=\mathcal{O}(r_{1}^{3})\\}$. One can verify
that $N^{a}_{1}$ and $\gamma_{c,1}$ agree near $p_{a}$ by using equation
(B.12) and the form of the blow-up transformation (B.4) to obtain

$\gamma_{c,1}=\big{(}\epsilon_{1}^{2/3}\zeta(-\epsilon_{1}^{-1/3}),\epsilon_{1},0\big{)}=\bigg{(}\frac{d_{0}}{b_{1}}\epsilon_{1}+\mathcal{O}(\epsilon_{1}^{2}),\epsilon_{1},0\bigg{)},$

in the limit $\epsilon_{1}\to 0$. Uniqueness of $N^{a}_{1}$ follows from the
equation $\epsilon_{1}^{\prime}|_{N_{a,1}}>0$. ∎

Entry and exit segments $\Sigma_{1}^{in/out}$ shown in Figure B.2 are defined
by

$\Sigma_{1}^{in}=\\{(u_{1},\epsilon_{1},r_{1})|r_{1}=\rho\\},\qquad\Sigma_{1}^{out}=\\{(u_{1},\epsilon_{1},r_{1})|\epsilon_{1}=\delta\\}.$

We also define the rectangle
$R_{1}=\\{|u_{1}|\leq\omega\\}\cap\Sigma_{1}^{in}$ for some sufficiently small
$\omega>0$, and let $I(\epsilon_{0})=R_{1}|_{\epsilon_{1}=\epsilon_{0}}$ for
each fixed $\epsilon_{0}\in[0,\delta]$ (see Figure B.2). The following result
summarises the dynamics near $p_{a}$.

######  Proposition B.3.

(cf. Proposition 2.8 in [44]). Given $\delta,\rho,\omega$ sufficiently small,
the transition map $\Pi_{1}:\Sigma_{1}^{in}\to\Sigma_{1}^{out}$ has the
following properties:

  * (i)

$\Pi_{1}(R_{1})$ is wedge-shaped in $\Sigma_{1}^{out}$.

  * (ii)

$\forall\epsilon_{1}\in(0,\delta]$ and $c<b_{1}$ fixed, $\exists K>0$ such
that $\Pi_{1}|_{I(\epsilon_{1})}$ is a contraction with contraction rate
bounded below by

$K\exp\bigg{[}-\frac{cb_{1}}{3a_{0}d_{0}}\bigg{(}\frac{1}{\epsilon_{1}}-\frac{1}{\delta}\bigg{)}\bigg{]}.$

###### Proof.

We need an expression for the transition time $T$ taken for solutions to pass
from $\Sigma_{1}^{in}$ to $\Sigma_{1}^{out}$ (cf. Lemma 2.7 in [44]). Direct
integration of the equation for $\epsilon_{1}^{\prime}$ in (B.20) gives the
leading order estimate

$T\sim\frac{b_{1}}{3a_{0}d_{0}}\bigg{(}\frac{1}{\epsilon_{1}}-\frac{1}{\delta}\bigg{)}.$
(B.21)

Statements (i) and (ii) follow from the expression (B.21) and Proposition B.2.
∎

Chart $K_{3}$. After an additional desingularisation
$dt=r_{3}^{-1}d\bar{t}_{3}$, we obtain the following in exit chart $K_{3}$:

$\displaystyle\begin{array}[]{lcl}u_{3}^{\prime}=b_{1}u_{3}+d_{0}\epsilon_{3}-2u_{3}\theta_{3}(u_{3},\epsilon_{3},r_{3})+r_{3}((b_{3}+b_{2}u_{3})u_{3}+d_{1}\epsilon_{3})+\mathcal{O}(r_{3}^{2}),\\\
\epsilon_{3}^{\prime}=-3\epsilon_{3}\theta_{3}(u_{3},\epsilon_{3},r_{3}),\\\
r_{3}^{\prime}=r_{3}\theta_{3}(u_{3},\epsilon_{3},r_{3}),\end{array}$ (B.25)

where the subscript notation indicates the use of chart-specific coordinates

$(x,u,\epsilon)=(r_{3},r_{3}^{2}u_{3},r_{3}^{3}\epsilon_{3}),$

and

$\theta_{3}(u_{3},\epsilon_{3},r_{3})=a_{0}u_{3}+r_{3}(a_{1}u_{3}+c_{0}\epsilon_{3})+\mathcal{O}(r_{3}^{2}).$

The portion of the equator $S^{1}$ visible in chart $K_{3}$ can be identified
with the invariant line $\\{(u_{3},0,0)|u_{3}\in\mathbb{R}\\}$, along which we
identify the partially hyperbolic singularity $p_{r}=(0,0,0)$, and hyperbolic
saddle $q_{out}=(b_{1}/2a_{0},0,0)$ shown in Figure B.1. We are interested in
the manner by which solutions leave the neighbourhood of $S^{2}_{+}$ near
$q_{out}$. Accordingly, we restrict to the set

$D_{3}:=\big{\\{}(u_{3},\epsilon_{3},r_{3})\big{|}(\epsilon_{3},r_{3})\in[0,\delta]\times[0,\rho]\big{\\}},$

and define the sections

$\Sigma_{3}^{in}=\big{\\{}(u_{3},r_{3},\epsilon_{3})\big{|}r_{3}\in[0,\rho],\epsilon_{3}=\delta\big{\\}},\qquad\Sigma_{3}^{out}=\big{\\{}(u_{3},r_{3},\epsilon_{3})\big{|}r_{3}=\rho,\epsilon_{3}\in[0,\delta]\big{\\}}.$

######  Proposition B.4.

(cf. Propostion 2.11 in [44]). The transition map
$\Pi_{3}:\Sigma_{3}^{in}\to\Sigma_{3}^{out}$ has form

$\Pi_{3}(u_{3},\delta,r_{3})=\bigg{(}\Pi_{31}(u_{3},\delta,r_{3}),\delta\bigg{(}\frac{r_{3}}{\rho}\bigg{)}^{3},\rho\bigg{)}^{T},$

where

$\Pi_{31}(u_{3},\delta,r_{3})\sim\frac{b_{1}}{2a_{0}}+\bigg{(}u_{3}-\frac{b_{1}}{2a_{0}}\bigg{)}\bigg{(}\frac{r_{3}}{\rho}\bigg{)}^{2}.$
(B.26)

In particular, if $q=\gamma_{c,3}\cap\Sigma_{3}^{in}$, then

$\Pi_{31}(q)\sim\frac{b_{1}}{2a_{0}}+\Omega_{0}\bigg{(}\frac{2d_{0}^{2}}{a_{0}b_{1}}\bigg{)}^{1/3}\bigg{(}\frac{r_{3}}{\rho}\bigg{)}^{2}\delta^{2/3}.$
(B.27)

###### Proof.

Although $q_{out}$ is hyperbolic, it has a resonance preventing a local
transformation into the leading order linear system. An explicit solution for
the leading order dynamics near $q_{out}$ can be obtained directly from
(B.25), though. Since $\theta_{3}(u_{3},\epsilon_{3},r_{3})\sim b_{1}/2>0$
near $q_{out}$ we can rescale
$d\bar{t}_{3}=(\theta_{3}(u_{3},\epsilon_{3},r_{3}))^{-1}d\tilde{t}$,
obtaining

$\displaystyle\begin{array}[]{lcl}u_{3}^{\prime}=\frac{b_{1}u_{3}+d_{0}\epsilon_{3}}{\theta_{3}(u_{3},\epsilon_{3},r_{3})}-2u_{3},\\\
\epsilon_{3}^{\prime}=-3\epsilon_{3},\\\ r_{3}^{\prime}=r_{3},\end{array}$
(B.31)

where the dash notation now denotes differentiation with respect to
$\tilde{t}_{3}$. The time $T$ taken for solutions of (B.31) to travel from
$\Sigma_{3}^{in}$ to $\Sigma_{3}^{out}$ can be determined from the equation
for $\epsilon_{3}^{\prime}$. Direct integration gives

$T=\ln\bigg{(}\frac{\rho}{r_{3}}\bigg{)},$

and the equation for $\Pi_{31}(u_{3},\delta,r_{3})$ is determined by expanding
the right hand side of the expression for $u_{3}^{\prime}$ near $q_{out}$,
which gives

$u_{3}^{\prime}\sim\frac{b_{1}}{a_{0}}-2u_{3}.$

Integrating and evaluating at $T$ yields equation (B.26). Finally, expression
(B.27) follows from the form for the image of $\gamma_{c}$ in chart $K_{3}$,
given by

$\gamma_{c,3}=\big{(}\epsilon_{3}^{2/3}\zeta(\epsilon_{3}^{-1/3}),\epsilon_{3},0\big{)}=\bigg{(}\frac{b_{1}}{2a_{0}}+\Omega_{0}\bigg{(}\frac{2d_{0}^{2}}{a_{0}b_{1}}\bigg{)}^{1/3}\epsilon_{3}^{2/3}+\mathcal{O}(\epsilon_{3}),\epsilon_{3},0\bigg{)}.$

∎

######  Remark B.1.

One can apply center manifold theory to prove the existence of a locally
invariant, two-dimensional repelling center manifold $W_{r,3}^{loc}$ tangent
to the center subspace at $p_{r}$, and appeal to similar arguments to those
provided in the preceding proof to understand the leading order dynamics near
$q_{in}$. We omit the details, since the dynamics near $p_{r}$ and $q_{in}$
are not of primary importance here.

Proof of Theorem 3.6. The proof follows conceptually the same arguments as
those given in Section 2.8 of [44]. We denote the transition maps between
charts $K_{i}$ and $K_{j}$ by $\kappa_{ij}$, and note that smoothness of the
maps $\kappa_{ij}$ follows from the fact that $S^{2}$ is a manifold.

###### Proof.

We need to track $W_{a,1}^{loc}$ and $R_{1}$ under the flow. Define a map
$\Pi:\Sigma_{1}^{in}\to\Sigma_{3}^{out}$ by the composition

$\Pi=\Pi_{3}\circ\kappa_{23}\circ\Pi_{2}\circ\kappa_{12}\circ\Pi_{1}.$

The map $\pi:\Sigma_{in}\to\Sigma_{out}$ from the statement of Theorem 3.6 is
given by blow-down $\pi=\Phi\circ\Pi\circ\Phi^{-1}$, for $\epsilon>0$.

By Proposition B.3 and smoothness of the transition maps $\kappa_{ij}$, the
image $\kappa_{12}\circ\Pi_{1}(R_{1}\cap W_{a,1}^{loc})$ is a smooth curve in
$\Sigma_{2}^{in}$, transverse to $\\{r_{2}=0\\}$. Since
$\kappa_{12}\circ\Pi_{1}(W_{a,1}^{loc}\cap\\{r_{1}=0\\})=\gamma_{c,2}\cap\Sigma_{2}^{in}$,
it follows from Proposition B.1 that
$\Pi_{2}\circ\kappa_{12}\circ\Pi_{1}(W_{a,1}^{loc}\cap R_{1})$ is a smooth
curve in $\Sigma_{2}^{out}$ of form

$\\{(\delta^{-1/3},s_{2}^{out}(r_{2}),r_{2})|r_{2}\in[0,\delta^{1/3}\rho]\\},$

where $s_{2}^{out}:[0,\delta^{1/3}\rho]\to\mathbb{R}$ is a smooth function
such that $\\{s_{2}^{out}(0)\\}=\gamma_{c,2}\cap\Sigma_{2}^{out}$. In
particular, this curve is transverse to $\\{r_{2}=0\\}$, and so the curve
$\kappa_{23}\circ\Pi_{2}\circ\kappa_{12}\circ\Pi_{1}(W_{a,1}^{loc}\cap R_{1})$
is transverse to $\\{r_{3}=0\\}$ in $\Sigma_{3}^{in}$. Proposition B.4 implies
that the image of the curve under $\Pi_{3}$ is of form

$\\{(s_{3}^{out}(\epsilon_{3}),\epsilon_{3},\rho)|\epsilon_{3}\in[0,\delta]\\},$

with

$s_{3}^{out}(\epsilon_{3})=\frac{b_{1}}{2a_{0}}+\mathcal{O}(\epsilon_{3}^{2/3}).$

Hence after applying the blow-down transformation we obtain

$|u_{s}-u_{l}|=\mathcal{O}(\epsilon^{2/3})\qquad\implies\qquad|M(\rho,u_{s})-M(\rho,u_{l})|=|y_{s}-y_{l}|=\mathcal{O}(\epsilon^{2/3}),$

where $u_{s}$, $u_{l}$ denote the $u$-coordinate of the intersections
$S^{a}_{\epsilon}\cap\Sigma_{out}$, $\mathcal{F}\cap\Sigma_{out}$
respectively, and we have used the fact that $M(\rho,u)=\mathcal{O}(u)$. This
proves assertion (1) in Theorem 3.6. Assertion (2) follows exactly as in [44],
only with Propositions B.3 and B.4 in place of Propositions 2.8 and 2.11
respectively in the cited work. ∎

######  Remark B.2.

Strictly speaking, one must also verify that the flow is regular across
$\\{(\bar{x},\bar{u},\bar{\epsilon})|\bar{x}=0,\bar{\epsilon}\geq
0\\}\times[0,\rho]$. This can be done in additional charts $\bar{u}=\pm 1$.

Figure B.3: Regular fold point dynamics on $S^{2}_{+}$, as studied in [44].
Compare with Figure B.1.

######  Remark B.3.

In [44], the authors identify a Ricatti equation in chart $K_{2}$ of the form

$\displaystyle\begin{array}[]{lcl}\tilde{x}_{2}^{\prime}=\tilde{x}_{2}^{2}-\tilde{y}_{2},\\\
\tilde{y}_{2}^{\prime}=-1,\end{array}$ (B.34)

in the limit $\tilde{r}_{2}\to 0$, where we have used the tilde notation to
distinguish their coordinates from ours. This system is related to the Ricatti
equation in (B.11) by the coordinate transformation

$(\tilde{x}_{2},\tilde{u}_{2})=(\tilde{x}_{2},\tilde{x}_{2}^{2}-\tilde{y}_{2})$

followed by

$\big{(}x_{2},u_{2},t_{2}\big{)}=(\alpha^{-1}\tilde{x}_{2},\beta^{-1}\tilde{u}_{2},\gamma^{-1}\tilde{t}_{2}),$

where

$\alpha=\bigg{(}\frac{b_{1}^{2}}{4a_{0}d_{0}}\bigg{)}^{1/3},\qquad\beta=\bigg{(}\frac{a_{0}b_{1}}{2d_{0}^{2}}\bigg{)}^{1/3},\qquad\gamma=\bigg{(}\frac{a_{0}b_{1}d_{0}}{2}\bigg{)}^{1/3}.$

The dynamics in the blow-up observed in [44] are sketched in B.3, which should
be compared with Figure B.1. Note that the segments chosen for the statement
and proof of Theorem 3.6 are not the same as those in [44].