# Normal form for singular Bautin bifurcation in a slow-fast system with
Holling type III functional response

Tapan Saha  Pranali Roy Chowdhury  Pallav Jyoti Pal  Malay Banerjee

###### Abstract

Over the last few decades, complex oscillations of slow-fast systems have been
a key area of research. In the theory of slow-fast systems, the location of
singular Hopf bifurcation and maximal canard is determined by computing the
first Lyapunov coefficient. In particular, the analysis of canards is based on
the genericity condition that the first Lyapunov coefficient must be non-zero.
This manuscript aims to further extend the results to the case where the first
Lyapunov coefficient vanishes. For that, the analytic expression of the second
Lyapunov coefficient and the investigation of the normal form for
codimension-2 singular Bautin bifurcation in a predator-prey system is done by
explicitly identifying the locally invertible parameter-dependent
transformations. A planar slow-fast predator-prey model with Holling type III
functional response is considered here, where the prey population growth is
affected by the weak Allee effect and the prey reproduces much faster than the
predator. Using geometric singular perturbation theory, normal form theory of
slow-fast systems, and blow-up technique, we provide a detailed mathematical
investigation of the system to show a variety of rich and complex nonlinear
dynamics including but not limited to the existence of canards, relaxation
oscillations, canard phenomena, singular Hopf bifurcation, and singular Bautin
bifurcation. Additionally, numerical simulations are conducted to support the
theoretical findings.

##  1 Introduction

Slow-fast systems or singularly perturbed systems of ordinary differential
equations evolving on two or more timescales are represented in the form

$\displaystyle\dot{x}\,=\,f\left(x,y,\epsilon\right),\,\,\dot{y}\,=\,\epsilon
g\left(x,y,\epsilon\right),$ (1)

where $x\in\mathbb{R}^{m}$ and $y\in\mathbb{R}^{n}$ ($m,n\geq 1$) are the fast
and the slow variables respectively, $f$ and $g$ are sufficiently smooth
functions, $0<\epsilon\ll 1$ is the singular perturbation parameter, and the
over dot $(~{}\dot{}~{})$ stands for derivative w.r.t. $t\in\mathbb{R}$. They
are known to exhibit a wide range of complicated oscillations of many physical
and applied systems including, but not limited to, canard cycles, relaxation
oscillations, and mixed mode oscillations [21]. The investigation of these
systems has received considerable attention in recent years [7, 10, 14, 17,
18, 20, 24, 31, 32].

In this article, we concentrate on slow-fast predator-prey systems of the type
(1) with $m=n=1$. Some examples of recent research in this field are as
follows: The author in [11] discussed and analysed the geometric singular
perturbation theory, three major Fenichel theorems, and their relevance to
biological modelling applications. In [16], the traditional Rosenzweig-
MacArthur (RM) model and the Mass Balance (MB) chemostat model are
investigated in the slow-fast framework. According to their analysis, the RM
model captures the dynamically rich canard explosion phenomenon, but the MB
model does not. Consequently, in the slow-fast setting, the predictions of the
MB Model are quite different from that of the RM model. The canard phenomenon
and number of canard cycles for a predator-prey system with Holling type III
and IV functional response was investigated by Li and Zhu [23]. In [2] the
author studies the dynamics of a multiple timescale predator-prey model where
the predator is a generalist who feeds on both the focal prey and the
functional response is Holling type III. With the help of the normal form
theory of slow-fast systems, the geometric singular perturbation theory, and
blow-up technique, the author [19] has investigated the existence of
homoclinic orbits, heteroclinic orbits, canard limit cycles, and relaxation
oscillations bifurcating from the singular homoclinic cycles.

In order to systematically analyze a slow-fast system, it is important to
recognize the different subprocesses occurring at the various time scales,
comprehend them, and then attempt to characterize the complete dynamics of the
full system based on the dynamical behaviour of the subsystems. The geometric
analysis pioneered by Fenichel [6] is based on this approach. Normal
hyperbolicity is a crucial property that critical manifolds may have. The
mathematical theory, known as geometric singular perturbation theory (GSPT) is
used to analyze the systems of the form (1) [5, 6, 20, 21]. To analyse slow-
fast systems with hyperbolic points, one may make use of the geometrical tools
and methods outlined by Fenichel theory [6]. However, other geometrical
approaches, such as the blow-up method, introduced by Dumortier and Roussarie
[5] and developed by Krupa and Szmolyan [19, 20], and slow-fast normal form
theory [1], are commonly used to explore the dynamics of a slow-fast system
with non-hyperbolic singularities. The blow-up technique includes de-
singularization of non-hyperbolic singularities. The key idea of such blow-up
transform is to reduce a slow-fast system into the form

$\displaystyle x^{\prime}=f(x,y,\mu),~{}~{}y^{\prime}=\epsilon
g(x,y,\mu),~{}~{}\epsilon^{\prime}=0,$ (2)

such that a non-hyperbolic equilibrium, say, $(x,y,\epsilon)=(0,0,0)$ (where
$D_{x}f(0,0,0)$ has all its eigenvalues equal to zero), is blown-up, for
example, to a sphere, to desingularize the equilibrium $(0,0,0)$. Generally,
different (overlapping) directional charts are used to investigate the
dynamics for the blown-up system and finally, a qualitative description of the
original system is obtained by blow-down. Using the above-mentioned technique,
one can identify the existence of canards, relaxation oscillations, canard
phenomena, singular Hopf bifurcation, singular Bautin bifurcation.

In this article, we consider a slow-fast predator-prey system with a weak
Allee effect in the prey ($x$) population growth and sigmoid functional
response of the form $f(x)=\frac{qx^{2}}{x^{2}+c}$, called Holling type-III,
where $c>0$ is the half-saturation constant for the predator $y$ that gauges
how abruptly $f(x)$ changes [9, 15]. Allee effect comprises both weak and
strong Allee effects [3, 29]. We note that the weak Allee effect growth
function, in contrast to the strong Allee effect, is always positive and has
no threshold value [4, 26, 30]. In this article, the weak Allee effect is
incorporated into the prey population growth, allowing the prey to grow even
at low population densities. Based upon the above considerations, the model
can be described as follows:

$\displaystyle\frac{dx}{dT}$ $\displaystyle=$ $\displaystyle
rx(1-\frac{x}{K})(x+m)-\frac{qx^{2}y}{x^{2}+c},$ (3a)
$\displaystyle\frac{dy}{dT}$ $\displaystyle=$
$\displaystyle\frac{px^{2}y}{x^{2}+c}-dy,$ (3b)

with $(c,d,K,m,p,q,r)\in\mathbb{R}^{7}_{+}$, where $x=x(t)$ and $y=y(t)$
denote the prey and predator population densities at time $t>0$, $c>0$ is the
half-saturation constant for the predator, $d$ is the predator’s natural death
rate, $K$ is the environmental carrying capacity, $r$ is the intrinsic per
capita growth rate of prey, $q$ denotes the maximum per capita consumption
rate, and $p$ denotes the conversion efficiency of consumed prey into new
predators [9, 12, 13]. Considering that the prey reproduces substantially much
faster than the predator and using the following scaling of the variables
given by,

$x=Ku,~{}~{}y=\frac{rK^{2}}{q}v,~{}~{}T=\frac{t}{rK},$

the model system (3) can then be rewritten in the following slow-fast
dimensionless form:

$\displaystyle\frac{du}{dt}$ $\displaystyle=$
$\displaystyle\left((1-u)(u+\theta)-\frac{uv}{u^{2}+\eta}\right)u,$ (4a)
$\displaystyle\frac{dv}{dt}$ $\displaystyle=$
$\displaystyle\epsilon\left(\frac{u^{2}}{u^{2}+\eta}-\delta\right)v,$ (4b)

where $\mu=(\delta,\theta,\eta)\in\mathbb{R}_{+}\times(0,1)^{2}$ and
$0<\epsilon\ll 1$ and subjected to the initial conditions $u(0)>0$ and
$v(0)>0$ such that $\delta=\frac{d}{p}$, $\theta=\frac{m}{K}$,
$\eta=\frac{c}{K^{2}}$ and $\epsilon=\frac{d}{p}$ are the dimensionless
parameters.

The main aim of this article is to analyze the phase-space of predator–prey
systems (4) in the context of a slow-fast framework and offer a detailed
description of the vast range of rich and complex local and global dynamical
behaviour, including, but not just restricted to, singular Hopf bifurcations,
singular Bautin bifurcations, canard cycles, canard phenomenon and occurrence
of relaxation oscillations via canard cycles known as “boom and bust cycle”
[27, 28]. It suggests that both populations may coexist with a predictable
pattern of population explosions and contractions. The main motivation of this
work is the derivation of the analytic form of the second Lyapunov coefficient
and a thorough study of the singular Bautin bifurcation for a system under
slow-fast framework [22, 25]. The normal form for singular Bautin bifurcation
has been derived here by explicitly finding the locally invertible parameter-
dependent transformations. So far as our knowledge goes, the explicit
derivation of locally invertible parameter-dependent transformations is not
available in the literature.

The remaining portion of the paper is organized as follows: Section 2 reports
the invariance, boundedness of the solution, existence, and local stability
analysis of equilibrium points of the system (4). The slow-fast system is
analyzed by separating system (4) into the fast and slow limiting subsystems
in Section 3. The singular bifurcations (singular Hopf and singular Bautin)
results are described in 4. We investigate the canard explosion and relaxation
oscillation in Section 5. In Section 6, numerical simulations are performed to
verify the major theoretical predictions. In the last Section 7, we discuss
our findings and provide our conclusions.

##  2 Positivity and boundedness of the solution

We state the following lemmas to ensure the positivity and boundedness of the
solution for the system (4). The proofs of the lemmas are omitted as it
directly follows from the approach discussed in [8].

######  Lemma 2.1.

The first quadrant $\mathbb{R}_{+}=\\{(u,v)\in\mathbb{R}^{2}|u\geq 0,v\geq
0\\}$ is invariant under the flow generated by the vector field
$\mathcal{F}=\hat{f}\frac{\partial}{\partial
u}+\epsilon\hat{g}\frac{\partial}{\partial v}$ where
$\hat{f}=u(1-u)(u+\theta)-\frac{u^{2}v}{u^{2}+\eta}$ and
$\hat{g}=\frac{u^{2}v}{u^{2}+\eta}-\delta v.$

######  Lemma 2.2.

All the solutions of the model system (4) initiated from an interior point of
$\mathbb{R}^{2}_{+}$ are bounded.

###  2.1 Equilibrium points and their local stability

The equilibria of the system (4) are given by the interaction of the prey and
predator nullcline as follows,

$(1-u)(u+\theta)(u^{2}+\eta)\,=\,uv,\,\,u=0,$

and

$u^{2}/(u^{2}+\eta)=\delta,\,\,v=0.$

The system has $E_{0}=(0,0)$ and $E_{1}=(1,0)$ as the boundary equilibria. The
interior equilibrium points are obtained when the non-trivial prey and
predator nullclines intersect in the interior of $\mathbb{R}^{2}_{+}.$ This
gives us at most one feasible interior equilibrium point $E_{*}=(u_{*},v_{*})$
where

$\displaystyle
u_{*}=\sqrt{\delta\eta/(1-\delta)},\,\,\,\,\,\,v_{*}=(1-u_{*})(u_{*}+\theta)(u_{*}^{2}+\eta)/u_{*}.$

The equilibrium $E_{*}$ is feasible if $\delta<\frac{1}{1+\eta}$. Linearizing
(4) in the neighbourhood of $E_{0}$ we obtain that the eigenvalues of the
Jacobian matrix are $0$ and $-\epsilon\delta.$ Because of the zero eigenvalue,
we cannot directly conclude about the dynamics of the system (4) near the
origin. The eigenvector corresponding to $0$ is the $u$ axis and the
eigenvector corresponding to $-\epsilon\delta$ is the $v$ axis. Therefore, any
trajectory starting very close to the $v$ axis gets attracted toward the
origin and repels away from the vicinity of the origin. The origin is thus an
attracting saddle node. We further linearize (4) in the neighbourhood of
$E_{1}$ and the eigenvalues of the Jacobian matrix are $-1-\theta<0$ and
$\epsilon\Big{(}\frac{1}{1+\eta}-\delta\Big{)}.$ For $\frac{1}{1+\eta}<\delta$
the equilibrium $E_{1}$ is a stable node, for $\frac{1}{1+\eta}>\delta$ the
equilibrium is a saddle point, and for $\delta=\frac{1}{1+\eta}$ the
equilibrium $E_{1}$ is an attracting saddle node. The eigenvector
corresponding to the eigenvalue $-1-\theta$ is the $u$ axis. Thus, combining
we get, the trajectory starting near the vertical axis gets attracted towards
the origin, and then gets attracted towards $E_{1}.$ The stability of the
equilibrium point $E_{*}$ depends on the geometry of the non-trivial prey
nullcline and the position of $E_{*}$ on it. If there does not exist any
feasible coexistence equilibrium point $E_{*}$ then the boundary equilibrium
$E_{1}$ is the global attractor, otherwise, it is a saddle point. The system
(4) thus undergoes a transcritical bifurcation at $\delta=\frac{1}{1+\eta}$.
The stability of the interior equilibrium point will be discussed in the
upcoming sections.

##  3 Analysis of slow-fast system

We consider a topologically equivalent system of the slow-fast system (4) by
re-scaling the time $dt\rightarrow(u^{2}+\eta)dt,$ such that $(u^{2}+\eta)>0$
as follows

$\displaystyle\frac{du}{dt}$
$\displaystyle=u(1-u)(u+\theta)(u^{2}+\eta)-u^{2}v=f(u,v),$ (5)
$\displaystyle\frac{dv}{dt}$ $\displaystyle=\epsilon(u^{2}v-\delta
v(u^{2}+\eta))=\epsilon g(u,v).$

The parameters are positive with $0<\epsilon\ll 1$. With the transformation
$\tau:=\epsilon t$, the system (5) is equivalent to

$\displaystyle\epsilon\frac{du}{d\tau}$
$\displaystyle=\left((1-u)(u+\theta)(u^{2}+\eta)-uv\right)u,$ (6)
$\displaystyle\frac{dv}{d\tau}$
$\displaystyle=\left(u^{2}-\delta(u^{2}+\eta)\right)v.$

The systems (5) and (6) are referred to as fast and slow systems with
$`t^{\prime}$ as the fast timescale and $`\tau^{\prime}$ as the slow
timescale. The variables $u$ and $v$ are known as fast and slow variables. In
the limiting sense of the fast system, that is, as $\epsilon=0,$ we obtain the
layer systems as

$\displaystyle\frac{du}{dt}=\left((1-u)(u+\theta)(u^{2}+\eta)-uv\right)u,\,\,\frac{dv}{dt}=0,$
(7)

and the slow system is reduced to

$\displaystyle
0=\left((1-u)(u+\theta)(u^{2}+\eta)-uv\right)u,\,\,\frac{dv}{d\tau}=\left(u^{2}-\delta(u^{2}+\eta)\right)v.$
(8)

The fast flow is along the horizontal line $v=$ constant, whereas the slow
flow is confined to the set

$\displaystyle M_{10}$ $\displaystyle=\\{(u,v)\in\mathbb{R}^{2}_{+}:u=0\\},$
$\displaystyle M_{20}$
$\displaystyle=\\{(u,v)\in\mathbb{R}^{2}_{+}:v=\phi(u)=\frac{1}{u}(1-u)(u+\theta)(u^{2}+\eta)\\},$

known as critical manifold. The critical manifold $M_{20}$ is a one-
dimensional curve that is normally hyperbolic whenever $\phi^{\prime}(u)>0$ or
$<0.$ The points at which $\phi^{\prime}(u)=0$ are the fold points of the
critical manifold. From $\phi^{\prime}(u)=0,$ we obtain,

$3u^{4}-2(1-\theta)u^{3}-(\theta-\eta)u^{2}+\eta\theta=0.$ (9)

Therefore, the fold points of the critical manifold are the roots of the above
quartic equation (if any). We claim that under any of the following parametric
restrictions P1: $\eta<\theta\leq 1,$ P2: $\eta\leq\theta<1,$ P3:
$\theta>\max\\{1,\eta\\},$ and P4: $\theta<\min\\{1,\eta\\}$ the critical
manifold $M_{20}$ can have at most two-feasible fold points where the
$u$-components are given by the positive solution of the quartic equation (9)
and $v=\frac{1}{u}(1-u)(u+\theta)(u^{2}+\eta).$ Let
$F(u)=3u^{4}-2(1-\theta)u^{3}-(\theta-\eta)u^{2}+\eta\theta,$ then
$F(0)=\eta\theta>0$ and $F(1)=(\theta+1)(\eta+1)>0.$ We argue that if the
function $F$ has a minimum for some $u_{\mathrm{min}}\in(0,1)$ and
$F(u_{\mathrm{min}})<0,$ then by continuity of the function $F$ we can say
that there exist exactly two roots of the equation $F=0$ which corresponds to
two fold points. For this we differentiate $F$ with respect to $u$ and from
$\frac{dF}{du}=0$ we obtain,

$u=0,\,\frac{1-\theta}{4}\pm\frac{\sqrt{9\theta^{2}+6\theta+9-24\eta}}{12}.$

We first consider the case (P1) and assume
$\Gamma:=9\theta^{2}+6\theta+9-24\eta>0$. This implies

$0<\frac{1-\theta}{4}+\frac{\sqrt{\Gamma}}{12}<1\,\,\text{and}\,\,\frac{1-\theta}{4}-\frac{\sqrt{\Gamma}}{12}<0.$

We then claim that the function $F$ attains its minimum at
$u_{\mathrm{min}}=\frac{1-\theta}{4}+\frac{\sqrt{\Gamma}}{12}.$ This gives

$\frac{d^{2}F}{du^{2}}\Big{|}_{u_{\mathrm{min}}}=\frac{\Gamma}{6}+\frac{(1-\theta)}{2}\sqrt{\Gamma}>0,\,\,\text{and}\,\,F(u_{\mathrm{min}})=\Lambda_{1}-\Lambda_{2}\sqrt{\Gamma},$
$\text{where}\,\,\Lambda_{1}=\frac{\eta}{8}(\theta^{2}+\frac{22}{3}\theta+1)-\frac{1}{96}(3(1+\theta^{4})+2(\theta^{2}+4\eta^{2})),\,\,\Lambda_{2}=\frac{1-\theta}{288}(3\theta^{2}+2\theta+3-8\eta).$

Therefore there exist exactly two feasible roots of the polynomial $F(u)$ in
the parametric region $\mathcal{R}_{1},$ where

$\mathcal{R}_{1}=\\{(\eta,\theta):0<\eta<\theta\leq
1,\,\,\frac{\sqrt{\Gamma}}{3}-3<\theta<1+\frac{\sqrt{\Gamma}}{3},\,\,\Gamma>\left(\frac{\Lambda_{1}}{\Lambda_{2}}\right)^{2}\\},$

and the two extrema of the curve $\phi(u)$ correspond to two fold points of
the critical manifold $M_{20}$. The case (P2) is similar to case (P1). Under
the restriction (P3) we consider sub-cases, either $\eta<1<\theta$ or
$1<\eta<\theta.$ We consider the first subcase. Since for
$\eta<1<\theta,\,\,\frac{1-\theta}{4}-\frac{\sqrt{\Gamma}}{12}<0.$ Thus two
fold points exist for $0<\frac{1-\theta}{4}+\frac{\sqrt{\Gamma}}{12}<1$ which
implies $\frac{\sqrt{\Gamma}}{3}-3<\theta<1+\frac{\sqrt{\Gamma}}{3},$ and the
minimum value must be negative, that is,
$F(u_{\mathrm{min}})=\Lambda_{1}+\Lambda_{2}\sqrt{\Gamma}<0,$ that is
$\Gamma<\Big{(}\frac{\Lambda_{1}}{\Lambda_{1}}\Big{)}^{2}.$ Therefore, in the
parametric region

$\mathcal{R}_{2}=\Big{\\{}(\eta,\theta):0<\eta<1<\theta,\,\,\frac{\sqrt{\Gamma}}{3}-3<\theta<1+\frac{\sqrt{\Gamma}}{3},\,\,\Gamma<\left(\frac{\Lambda_{1}}{\Lambda_{2}}\right)^{2}\Big{\\}},$

we obtain two fold points in the interior of $\mathbb{R}_{+}^{2}.$ Similarly
under the restriction (P4) there can be two sub-cases, either $\theta<1<\eta$
or $\theta<\eta<1.$ We consider the first sub-case and assume that for
$\theta<\eta,$ $\Gamma>0.$ This gives us

$\theta<\eta<\frac{1}{8}(3\theta^{2}+2\theta+3).$

Since $u_{\mathrm{min}}\in(0,1)$ such that $F({u_{\mathrm{min}}})<0,$ we must
have

$\frac{\sqrt{\Gamma}}{3}-3<\theta<1,\,\,\text{and}\,\,\Gamma>\
\left(\frac{\Lambda_{1}}{\Lambda_{2}}\right)^{2}.$

Therefore under the restriction (P4) we obtain exactly two fold points for

$\mathcal{R}_{3}=\Big{\\{}(\eta,\theta):0<\theta<1<\eta,\,\,\frac{\sqrt{\Gamma}}{3}-3<\theta<1,\,\,\Gamma>\left(\frac{\Lambda_{1}}{\Lambda_{2}}\right)^{2}\Big{\\}}.$

We now consider region $\mathcal{R}$ such as
$\mathcal{R}=\mathcal{R}_{1}\cup\mathcal{R}_{2}\cup\mathcal{R}_{3}.$ Thus for
$(\eta,\theta)\in\mathcal{R}$ the critical manifold is a $S$-shaped curve with
exactly two fold points. Let the two fold points be $P(u_{m},v_{m})$ (local
minimum of the critical manifold $M_{20}$) and $Q(u_{M},v_{M})$ (local maximum
of the critical manifold $M_{20}$). Depending on the structure of the critical
manifold $M_{20},$ that as $\lim_{u\rightarrow 0}\phi(u)\rightarrow\infty$ and
$\lim_{u\rightarrow 1}\phi(u)\rightarrow 0$ we can infer that $u_{m}<u_{M}$.
It consists of three branches, namely, $S_{0}^{l}$, $S_{0}^{m}$ and
$S_{0}^{r}$, where

$\displaystyle S_{0}^{l}$
$\displaystyle=M_{20}\cap\left\\{(u,v)\in\mathbb{R}^{2}_{+}\,\middle|\,0<u<u_{m}\right\\},$
$\displaystyle S_{0}^{m}$
$\displaystyle=M_{20}\cap\left\\{(u,v)\in\mathbb{R}^{2}_{+}\,\middle|\,u_{m}<u<u_{M}\right\\},$
$\displaystyle S_{0}^{r}$
$\displaystyle=M_{20}\cap\left\\{(u,v)\in\mathbb{R}^{2}_{+}\,\middle|\,u_{M}<u<1\right\\}.$

From the geometry of the critical manifold $M_{20}$, we conclude that the sub-
manifolds $S_{0}^{l}$ and $S_{0}^{r}$ are normally hyperbolic attracting
whereas $S_{0}^{m}$ is normally hyperbolic repelling. When the non-trivial
predator nullcline $u=\sqrt{\frac{\delta\eta}{1-\delta}}$ intersects the non-
trivial prey nullcline $M_{20}$ on the left and right branch of $M_{20},$ that
is on $S^{l}_{0}$ or $S^{r}_{0}$ then the equilibrium is stable. Whereas if
$E_{*}$ lies on $S^{m}_{0}$, the equilibrium is unstable. We show the shape of
the critical manifold along with the fold points in Fig. 1.

Figure 1: A representation of the critical manifold $M_{20}$ with the left
($S^{l}_{0}$), middle ($S^{m}_{0}$), and right ($S^{r}_{0}$) branch. The blue
dots $P$ and $Q$ are the fold points of the critical manifold.

##  4 Singular Bifurcations

In this section, we focus our study on bifurcation structures in the
singularly perturbed system (5). Throughout the analysis, we assume that
$(\eta,\theta)\in\mathcal{R}$ so that the critical manifold $M_{20}$ is an
$S$-shaped curve having two non-degenerate fold points $P(u_{m},v_{m})$ and
$Q(u_{M},v_{M})$. We will show that when the stable equilibrium $E_{*}$ on the
left branch $S_{0}^{l}$ crosses the fold point $P$ for some threshold
parameter value, the system (5) undergoes a Hopf bifurcation which is singular
in the sense that at the Hopf bifurcation threshold, the eigenvalues of the
Jacobian matrix become singular in the singular limit $\epsilon\to 0$.
Subsequently, we will show that there will be a change in the criticality of
the singular Hopf bifurcation and the system undergoes a codimension-2 Bautin
bifurcation which we termed here as singular Bautin bifurcation as the
corresponding Hopf bifurcation is singular.

###  4.1 Singular Hopf Bifurcation

We choose $\delta$ as the bifurcation parameter. We assume that for
$\delta=\delta_{*}$, the equilibrium $E_{*}$ passes through the fold point
$P(u_{m},v_{m})$ of the critical manifold $M_{20}$, where $\delta_{*}$ is
implicitly given by the equation $F(u_{*})=0.$ Consequently, we have the
following

$\displaystyle
f(u_{m},v_{m},\delta_{*})=0,\,\,g(u_{m},v_{m},\delta_{*})=0,\,\,\text{and},\,\frac{\partial
f}{\partial u}(u_{m},v_{m},\delta_{*})=0.$ (10)

We further assume that

$\displaystyle\frac{\partial f}{\partial v}(u_{m},v_{m},\delta_{*})\neq
0,\,\,\frac{\partial^{2}f}{\partial u^{2}}(u_{m},v_{m},\delta_{*})\neq
0,\,\,\frac{\partial g}{\partial u}(u_{m},v_{m},\delta_{*})\neq
0\,\,\,\frac{\partial g}{\partial\delta}(u_{m},v_{m},\delta_{*})\neq 0.$ (11)

With the above assumption, the fold point $P$ is now the non-degenerate canard
point or the singular contact point of the system. Using the transformation,
$\bar{u}=u-u_{m},\,\,\bar{v}=v-v_{m},\,\,\lambda=\delta-\delta_{*},$ we write
the system (5) into the following form

$\displaystyle\frac{d\bar{u}}{dt}$
$\displaystyle=a_{01}\bar{v}+\frac{a_{20}}{2}\bar{u}^{2}+a_{11}\bar{u}\bar{v}+\frac{a_{30}}{6}\bar{u}^{3}+\frac{a_{21}}{2}\bar{u}^{2}\bar{v}+\frac{a_{40}}{24}\bar{u}^{4}+\frac{a_{50}}{120}\bar{u}^{5},$
(12) $\displaystyle\frac{d\bar{v}}{dt}$
$\displaystyle=\epsilon(b_{10}\bar{u}+\frac{b_{20}}{2}\bar{u}^{2}+b_{11}\bar{u}\bar{v}-v_{m}(u_{m}^{2}+\eta)\lambda$
$\displaystyle-2u_{m}v_{m}\lambda\bar{u}-(u_{m}^{2}+\eta)\lambda\bar{v}+\frac{b_{21}}{2}\bar{u}^{2}\bar{v}-v_{m}\lambda\bar{u}^{2}-\lambda\bar{u}^{2}\bar{v}),$

where $\displaystyle a_{ij}=\left.\frac{\partial^{i+j}f}{\partial
u^{i}\partial v^{j}}\right|_{(u_{m},v_{m},\delta_{*})}$ and $\displaystyle
b_{ij}=\left.\frac{\partial^{i+j}g}{\partial u^{i}\partial
v^{j}}\right|_{(u_{m},v_{m},\delta_{*})}$, i.e.,

$\displaystyle a_{10}=0,\,a_{01}=-u_{m}^{2},\,a_{50}=-120,\,b_{01}=0,$
$\displaystyle
a_{20}=\frac{4\delta_{*}\eta}{1-\delta_{*}}(3-3\theta-5u_{m})+2\eta(1-3u_{m}-\theta)+6\theta
u_{m}-2v_{m},\,a_{11}=-2u_{m},$ $\displaystyle
a_{30}=24(1-\theta)u_{m}+6(\theta-\eta)-\frac{60\delta_{*}\eta}{1-\delta_{*}},\,a_{21}=-2,\,a_{40}=24(1-\theta)-120u_{m},$
$\displaystyle
b_{10}=2(1-\delta_{*})u_{m}v_{m},\,\,b_{20}=2(1-\delta_{*})v_{m},\,\,b_{11}=2(1-\delta_{*})u_{m},\,\,b_{21}=2(1-\delta_{*}).$

One of the fundamental tools to study the flow of (12) near the canard point
is to use the blow-up transformation

$\Phi:\mathbb{S}^{3}\times[0,\rho]\to\mathbb{R}^{4}$

given by $\bar{u}=r^{\prime}u^{\prime},~{}~{}\bar{v}=r^{\prime
2}v^{\prime},~{}~{}\lambda=r^{\prime}\lambda^{\prime}~{}~{}\epsilon=r^{\prime
2}\epsilon^{\prime},~{}\text{where}~{}(u^{\prime},v^{\prime},\lambda^{\prime},\epsilon^{\prime})\in\mathbb{S}^{3}.$
The blow-up transformation, desingularizes the vector field (12). One can then
use the tools of dynamical systems to investigate the dynamics of the blown-up
vector field in different charts of the manifold
$\mathbb{S}^{3}\times[0,\rho]$ and connecting the results in different charts,
the final results are derived by “blow down” [20, 21]. The dynamics in the
chart $K_{2}$ which describes a neighbourhood of the upper half-sphere are
only necessary to study singular Hopf bifurcation. The blow-up transformation
in the chart $K_{2}$ is given by

$\displaystyle\bar{u}=r_{2}u_{2},~{}~{}~{}\bar{v}=r_{2}^{2}v_{2},~{}~{}~{}\lambda=r_{2}\lambda_{2}~{}~{}~{}\epsilon=r_{2}^{2}.$
(13)

Using the above transformation and dividing by $r_{2}$, the system (12)
transforms into the following

$\displaystyle\frac{du_{2}}{dt}$
$\displaystyle=a_{01}v_{2}+\frac{a_{20}}{2}u_{2}^{2}+\left(a_{11}u_{2}v_{2}+\frac{a_{30}}{6}u_{2}^{3}\right)r_{2}+\left(\frac{a_{21}}{2}u_{2}^{2}v_{2}+\frac{a_{40}}{24}u_{2}^{4}\right)r_{2}^{2}+\frac{a_{50}}{120}r_{2}^{3}u_{2}^{5},$
(14) $\displaystyle\frac{dv_{2}}{dt}$
$\displaystyle=b_{10}u_{2}-v_{m}(u_{m}^{2}+\eta)\lambda_{2}-2u_{m}v_{m}r_{2}\lambda_{2}u_{2}+\frac{b_{20}}{2}r_{2}u_{2}^{2}+b_{11}r_{2}^{2}u_{2}v_{2}$
$\displaystyle-(u_{m}^{2}+\eta)r_{2}^{2}\lambda_{2}v_{2}-v_{m}r_{2}^{2}\lambda_{2}u_{2}^{2}+\frac{b_{21}}{2}r_{2}^{3}u_{2}^{2}v_{2}-r_{2}^{4}\lambda_{2}u_{2}^{2}v_{2}.$

The system (14) has an equilibrium $P_{2}(u_{2e},v_{2e})$ with

$\displaystyle u_{2e}$
$\displaystyle=\frac{v_{m}(u_{m}^{2}+\eta)}{b_{10}}\lambda_{2}+\frac{1}{b_{10}^{2}}\left(2u_{m}v_{m}^{2}(u_{m}^{2}+\eta)-\frac{b_{20}}{2b_{10}}v_{m}^{2}(u_{m}^{2}+\eta)^{2}\right)r_{2}\lambda_{2}^{2}+\mathcal{O}(4),$
$\displaystyle v_{2e}$
$\displaystyle=-\frac{a_{20}}{2a_{01}b_{10}^{2}}v_{m}^{2}(u_{m}^{2}+\eta)^{2}\lambda_{2}^{2}+\mathcal{O}(4),$

where $\mathcal{O}(4)=\mathcal{O}(|\lambda_{2}+r_{2}|^{4})$. The Jacobian
matrix at $P_{2}$ is given by

$\displaystyle J=\left[\begin{array}[]{cc}\alpha_{11}&\alpha_{12}\\\
\alpha_{21}&\alpha_{22}\end{array}\right]$

where

$\displaystyle\alpha_{11}$
$\displaystyle=\frac{a_{20}}{b_{10}}v_{m}(u_{m}^{2}+\eta)\lambda_{2}+\frac{a_{20}}{b_{10}^{2}}\left(2u_{m}v_{m}^{2}(u_{m}^{2}+\eta)-\frac{b_{20}}{2b_{10}}v_{m}^{2}(u_{m}^{2}+\eta)^{2}\right)r_{2}\lambda_{2}^{2}$
$\displaystyle\hskip
113.81102pt+\left(a_{30}-\frac{a_{20}a_{11}}{a_{01}}\right)\frac{v_{m}^{2}(u_{m}^{2}+\eta)^{2}}{2b_{10}^{2}}r_{2}\lambda_{2}^{2}+\mathcal{O}(4),$
$\displaystyle\alpha_{12}$
$\displaystyle=a_{01}+\frac{a_{11}}{b_{10}}v_{m}(u_{m}^{2}+\eta)r_{2}\lambda_{2}+\mathcal{O}(4),$
$\displaystyle\alpha_{21}=b_{10}-2u_{m}v_{m}r_{2}\lambda_{2}+\frac{b_{20}}{b_{10}}v_{m}(u_{m}^{2}+\eta)r_{2}\lambda_{2}+\mathcal{O}(4),$
$\displaystyle\alpha_{22}$
$\displaystyle=\frac{b_{11}}{b_{10}}v_{m}(u_{m}^{2}+\eta)r_{2}^{2}\lambda_{2}-(u_{m}^{2}+\eta)r_{2}^{2}\lambda_{2}+\mathcal{O}(4).$

Now, using the transformation $x=u_{2}-u_{2e}$, $y=v_{2}-v_{2e}$ and expanding
into Taylor’s series, we write the system (14) as

$\displaystyle\frac{dx}{dt}$
$\displaystyle=a_{10}^{\prime}x+a_{01}^{\prime}y+\frac{a_{20}^{\prime}}{2}x^{2}+a_{11}^{\prime}xy+\frac{a_{30}^{\prime}}{6}x^{3}+\frac{a_{21}^{\prime}}{2}x^{2}y+\frac{a_{40}^{\prime}}{24}x^{4}+\frac{a_{50}^{\prime}}{120}x^{5},$
(15a) $\displaystyle\frac{dy}{dt}$
$\displaystyle=b_{10}^{\prime}x+b_{01}^{\prime}y+\frac{b_{20}^{\prime}}{2}x^{2}+b_{11}^{\prime}xy+\frac{b_{21}^{\prime}}{2}x^{2}y,$
(15b)

where,

$\displaystyle a_{10}^{\prime}$
$\displaystyle=\alpha_{11},\,a_{01}^{\prime}=\alpha_{12},\,b_{10}^{\prime}=\alpha_{21},\,b_{01}^{\prime}=\alpha_{22},\,a_{20}^{\prime}=a_{20}+\frac{a_{30}}{b_{10}}v_{m}(u_{m}^{2}+\eta)r_{2}\lambda_{2}+\mathcal{O}(4),$
$\displaystyle a_{11}^{\prime}$
$\displaystyle=a_{11}r_{2}+\frac{a_{21}}{b_{10}}v_{m}(u_{m}^{2}+\eta)r_{2}\lambda_{2}+\mathcal{O}(4),$
$\displaystyle a_{30}^{\prime}$
$\displaystyle=a_{30}r_{2}+\frac{a_{40}}{b_{10}}v_{m}(u_{m}^{2}+\eta)r_{2}\lambda_{2}+\mathcal{O}(4),$
$\displaystyle a_{21}^{\prime}$
$\displaystyle=a_{21}r_{2}^{2},\,a_{40}^{\prime}=a_{40}r_{2}^{2}+\mathcal{O}(4),\,a_{50}^{\prime}=a_{50}r_{2}^{3},\,b_{20}^{\prime}=b_{20}r_{2}-2v_{m}r_{2}\lambda_{2}+\mathcal{O}(4),$
$\displaystyle b_{11}^{\prime}$
$\displaystyle=b_{11}r_{2}^{2}+\mathcal{O}(4),\,\text{and}\,b_{21}^{\prime}=b_{21}r_{2}^{3}+\mathcal{O}(4).$

For the Hopf bifurcation, Trace $J=0$ gives rise to $\lambda_{2}=0$ and in
terms of original parameters this transforms into

$\delta=\delta_{H}=\delta_{*}.$

Thus, we observe that the Hopf bifurcation threshold does not depend on
$\epsilon$. For $|\delta-\delta_{H}|$ small, the complex conjugate eigenvalues
of the Jacobian matrix $J$ are given by $\mu$ and $\bar{\mu}$, where

$\displaystyle\mu=\alpha+i\beta,\,\,\,\alpha=\frac{\alpha_{11}+\alpha_{22}}{2},\,\,\,\beta=\frac{1}{2}\sqrt{4(\alpha_{11}\alpha_{22}-\alpha_{12}\alpha_{21})-(\alpha_{11}+\alpha_{22})^{2}}.$

Consequently, at $\delta=\delta_{H}$, we have
$\alpha(\delta_{H})=0,\,\,\beta(\delta_{H})=\beta_{0}+\mathcal{O}(4)$ where
$\beta_{0}=\sqrt{-a_{01}b_{10}}$. We also have the following

$\displaystyle\frac{d}{d\delta}\left(Re(\mu(\delta))\right)\bigg{|}_{\delta=\delta_{H}}$
$\displaystyle\neq 0.$

Let $q=\left(q_{1},q_{2}\right)^{T}\in\mathbb{C}$ be an eigenvector of $J$
corresponding to the eigenvalue $\mu$ and $p=\left(p_{1},p_{2}\right)^{T}$ be
an eigenvector of $J^{T}$ corresponding to $\bar{\mu}$ such that $\langle
p,q\rangle=\bar{p}_{1}q_{1}+\bar{p}_{2}q_{2}=1$. We then have

$\displaystyle\left[\begin{array}[]{c}q_{1}\\\
q_{2}\end{array}\right]=\left[\begin{array}[]{c}\alpha_{12}\\\
\mu-\alpha_{11}\end{array}\right],~{}~{}~{}\left[\begin{array}[]{c}p_{1}\\\
p_{2}\end{array}\right]=\frac{1}{2\alpha_{12}\bar{\mu}-\alpha_{12}(\alpha_{11}+\alpha_{22})}\left[\begin{array}[]{c}\bar{\mu}-\alpha_{22}\\\
\alpha_{12}\end{array}\right]$

With the help of the transformation
$z=\bar{p}_{1}x+\bar{p}_{2}y,\,\,\,z\in\mathbb{C},$ the system (15) can be
written as

$\displaystyle\dot{z}=\mu
z+\sum_{k+l=2}^{5}\frac{g_{kl}}{k!l!}z^{k}\bar{z}^{l},$ (16)

where the expressions of $g_{kl}$’s are given by

$\displaystyle g_{11}$
$\displaystyle=\bar{p}_{1}(a_{20}^{\prime}q_{1}\bar{q}_{1}+a_{11}^{\prime}(q_{1}\bar{q}_{2}+q_{2}\bar{q}_{1})+\bar{p}_{2}(b_{20}^{\prime}q_{1}\bar{q}_{1}+b_{11}^{\prime}(q_{1}\bar{q}_{2}+q_{2}\bar{q}_{1}),\,g_{05}=a_{50}^{\prime}\bar{p}_{1}\bar{q}_{1}^{5},$
$\displaystyle g_{20}$
$\displaystyle=\bar{p}_{1}(a_{20}^{\prime}q_{1}^{2}+2a_{11}^{\prime}q_{1}q_{2})+\bar{p}_{2}(b_{20}^{\prime}q_{1}^{2}+2b_{11}^{\prime}q_{1}q_{2}),$
$\displaystyle g_{02}$
$\displaystyle=\bar{p}_{1}(a_{20}^{\prime}\bar{q}_{1}^{2}+2a_{11}^{\prime}\bar{q}_{1}\bar{q}_{2})+\bar{p}_{2}(b_{20}^{\prime}\bar{q}_{1}^{2}+2b_{11}^{\prime}\bar{q}_{1}\bar{q}_{2}),$
$\displaystyle g_{21}$
$\displaystyle=\bar{p}_{1}(a_{30}^{\prime}q_{1}^{2}\bar{q}_{1}+a_{21}^{\prime}(q_{1}^{2}\bar{q}_{2}+2q_{1}\bar{q}_{1}q_{2}))+b_{21}^{\prime}\bar{p}_{2}(q_{1}^{2}\bar{q}_{2}+2q_{1}\bar{q}_{1}q_{2}),$
$\displaystyle g_{22}$
$\displaystyle=a_{40}^{\prime}\bar{p}_{1}q_{1}^{2}\bar{q}_{1}^{2},\,\,g_{23}=a_{50}^{\prime}\bar{p}_{1}q_{1}^{2}\bar{q}_{1}^{3},$
$\displaystyle g_{12}$
$\displaystyle=\bar{p}_{1}(a_{30}^{\prime}q_{1}\bar{q}_{1}^{2}+a_{21}^{\prime}(\bar{q}_{1}^{2}q_{2}+2q_{1}\bar{q}_{1}\bar{q}_{2}))+b_{21}^{\prime}\bar{p_{2}}(\bar{q}_{1}^{2}q_{2}+2q_{1}\bar{q}_{1}\bar{q}_{2}),$
$\displaystyle g_{13}$
$\displaystyle=a_{40}^{\prime}\bar{p}_{1}q_{1}\bar{q}_{1}^{3},\,\,g_{14}=a_{50}^{\prime}\bar{p}_{1}q_{1}\bar{q}_{1}^{4},\,\,g_{30}=\bar{p}_{1}(a_{30}^{\prime}q_{1}^{3}+3a_{21}^{\prime}q_{1}^{2}q_{2})+3b_{21}^{\prime}\bar{p}_{2}q_{1}^{2}q_{2},$
$\displaystyle g_{03}$
$\displaystyle=\bar{p}_{1}(a_{30}^{\prime}\bar{q_{1}}^{3}+3a_{21}^{\prime}\bar{q}_{1}^{2}\bar{q}_{2})+3b_{21}^{\prime}\bar{p}_{2}\bar{q}_{1}^{2}\bar{q}_{2},\,\,g_{40}=a_{40}^{\prime}\bar{p}_{1}q_{1}^{4},\,\,g_{04}=a_{40}^{\prime}\bar{p}_{1}\bar{q}_{1}^{4},$
$\displaystyle g_{50}$
$\displaystyle=a_{50}^{\prime}\bar{p}_{1}q_{1}^{5},\,\,g_{41}=a_{50}^{\prime}\bar{p}_{1}q_{1}^{4}\bar{q}_{1},\,\,g_{32}=a_{50}^{\prime}\bar{p}_{1}q_{1}^{3}\bar{q}_{1}^{2},\,g_{31}=a_{40}^{\prime}\bar{p}_{1}q_{1}^{3}\bar{q}_{1}.$

Using the following locally invertible parameter-dependent change of complex
coordinate,

$\displaystyle
z=w+\frac{h_{20}}{2}w^{2}+h_{11}w\bar{w}+\frac{h_{02}}{2}\bar{w}^{2}+\frac{h_{30}}{6}w^{3}+\frac{h_{12}}{2}w\bar{w}^{2}+\frac{h_{03}}{6}\bar{w}^{3},$
(18)

the system (16) can be written in a neighbourhood of $\delta=\delta_{H}$ as

$\displaystyle\dot{w}=\mu
w+c_{1}w^{2}\bar{w}+\mathcal{O}\left(|w|^{4}\right),$ (19)

where

$\displaystyle c_{1}$
$\displaystyle=\frac{g_{20}g_{11}\left(2\mu+\bar{\mu}\right)}{2|\mu|^{2}}+\frac{|g_{11}|^{2}}{\mu}+\frac{|g_{02}|^{2}}{2\left(2\mu-\bar{\mu}\right)}+\frac{g_{21}}{2},$
(20) $\displaystyle
h_{20}=\frac{g_{20}}{\mu},\,\,h_{11}=\frac{g_{11}}{\bar{\mu}},\,\,\,h_{02}=\frac{g_{02}}{2\bar{\mu}-\mu},$
$\displaystyle h_{30}$
$\displaystyle=\frac{3}{\mu}\left(\frac{g_{20}h_{20}}{2}+\frac{g_{11}\overline{h}_{02}}{2}+\frac{g_{30}}{6}\right),\,\,h_{12}=\frac{1}{2\bar{\mu}}\left(g_{20}h_{02}+2g_{11}h_{11}+g_{11}\overline{h}_{20}+g_{12}\right),$
$\displaystyle h_{03}$
$\displaystyle=\frac{1}{3\bar{\mu}-\mu}\left(g_{03}+3g_{11}h_{02}+3g_{02}\overline{h}_{20}\right),$

and the first Lyapunov coefficient using [22] is given by

$\displaystyle L_{1}(\delta_{H})$
$\displaystyle=\frac{\Re\left(c_{1}\right)}{\beta}\bigg{|}_{\delta=\delta_{H}}=\frac{1}{2\beta^{2}}\Re\left(ig_{20}g_{11}+\beta_{0}g_{21}\right)\bigg{|}_{\delta=\delta_{H}},$
(21)

where $\Re(\cdot)$ stands for the real part of $(\cdot)$.

At $\delta=\delta_{H}$,

$\displaystyle\left[\begin{array}[]{c}q_{1}(\delta_{H})\\\
q_{2}(\delta_{H})\end{array}\right]=\left[\begin{array}[]{c}a_{01}\\\
i\beta_{0}\end{array}\right],~{}~{}~{}\left[\begin{array}[]{c}p_{1}(\delta_{H})\\\
p_{2}(\delta_{H})\end{array}\right]=\left[\begin{array}[]{c}\frac{1}{2a_{01}}\\\
\frac{i}{2\beta_{0}}\end{array}\right],$

and thus, the first Lyapunov coefficient is

$\displaystyle L_{1}(\sqrt{\epsilon})$
$\displaystyle=-\frac{a_{01}}{4\beta_{0}b_{10}}A\sqrt{\epsilon},$ (22)

where

$\displaystyle A$
$\displaystyle=a_{01}a_{20}b_{20}-a_{01}a_{30}b_{10}+a_{11}a_{20}b_{10}{\color[rgb]{1,0,0}.}$
(23)

Thus, we see that all the conditions of the Hopf bifurcation theorem are
satisfied and consequently, the system (5) undergoes a Hopf bifurcation at
$\delta=\delta_{H}(\sqrt{\epsilon})=\delta_{*}.$ It also follows that for
$0<\epsilon\ll 1$, the leading order coefficient of $L_{1}$ i.e., $A$
determines the criticality of the singular Hopf bifurcation. The criticality
changes if $A$ changes its sign from negative to positive and consequently,
the singular Hopf bifurcation is degenerate if $A=0$. We now summarize the
above-mentioned results by the following theorem.

######  Theorem 4.1.

Let, $(\theta,\eta)\in\mathcal{R}$ and the conditions (10) and (11) hold. Then
$\exists$ $\epsilon_{0}>0$ and $\delta_{0}>0$ such that for
$0<\epsilon<\epsilon_{0}$ and $|\delta-\delta_{*}|<\delta_{0}$, the system (5)
has an equilibrium point $P_{2}$ in a neighbourhood of the fold point $P$
which converges to $P$ as $(\epsilon,\delta)\to(0,\delta_{*})$. The system (5)
undergoes a singular Hopf bifurcation at

$\displaystyle\delta_{H}(\sqrt{\epsilon})=\delta_{*}+\mathcal{O}(\epsilon^{\frac{5}{2}})$
(24)

The Hopf bifurcation is non-degenerate when $A\neq 0$. It is supercritical if
$A<0$ and sub-critical if $A>0$ where $A$ is given by (23).

###  4.2 Singular Bautin Bifurcation

We observe that for $\theta=\theta_{B}$, $A=0$ where $\theta_{B}$ is given by

$\displaystyle\theta=\theta_{B}=-\frac{{u_{m}}\left(\left(-1+\delta^{*}\right){u_{m}}^{3}+\left(-5\delta^{*}+5\right){u_{m}}^{2}-\eta\left(-1+\delta^{*}\right){u_{m}}-6\delta^{*}\eta\right)}{\left(5\delta^{*}-5\right){u_{m}}^{3}+\left(-1+\delta^{*}\right){u_{m}}^{2}+6\delta^{*}\eta{u_{m}}-\eta\left(-1+\delta^{*}\right)}.$
(25)

Consequently, it follows from theorem 4.1 that the Hopf bifurcation is
degenerate i.e., there is a change in the criticality of the singular Hopf
bifurcation as $\theta$ passes through $\theta=\theta_{B}$. We now proceed to
show that in such a case, the system undergoes a codimension-2 Bautin
bifurcation by computing the second Lyapunov coefficient, which determines the
criticality of the Bautin bifurcation.

To derive the normal form for Bautin bifurcation, we repeat the same process
as carried out in the previous section, i.e., we use the following locally
invertible parameter-dependent change of complex coordinate

$\displaystyle z$
$\displaystyle=w+\frac{h_{20}}{2}w^{2}+h_{11}w\bar{w}+\frac{h_{02}}{2}\bar{w}^{2}+\frac{h_{30}}{6}w^{3}+\frac{h_{12}}{2}w\bar{w}^{2}+\frac{h_{03}}{6}\bar{w}^{3}+\frac{h_{40}}{24}w^{4}$
$\displaystyle\hskip
28.45274pt+\frac{h_{31}}{6}w^{3}\bar{w}+\frac{h_{22}}{4}w^{2}\bar{w}^{2}+\frac{h_{13}}{6}w\bar{w}^{3}+\frac{h_{04}}{24}\bar{w}^{4}+\frac{h_{50}}{120}w^{5}$
$\displaystyle\hskip
56.9055pt+\frac{h_{41}}{24}w^{4}\bar{w}+\frac{h_{23}}{12}w^{2}\bar{w}^{3}+\frac{h_{14}}{24}w\bar{w}^{4}+\frac{h_{05}}{120}\bar{w}^{5},$

to eliminate all the quadratic terms and bi-quadratic terms of the equation
(16) and left with only resonant cubic and fifth-degree terms. The equation
(16) then transforms into,

$\displaystyle\dot{w}=\mu
w+c_{1}w^{2}\bar{w}+c_{2}w^{3}\bar{w}^{2}+\mathcal{O}\left(|w|^{6}\right),$
(26)

where the expressions for $c_{1}$, $h_{ij},2\leq i+j\leq 3$ are given by (20)
and the expressions for $c_{2}$, $h_{ij},i+j=4$ are given below. We have not
mentioned the expressions for $h_{ij},i+j=5$ here, as these coefficients are
no longer required to compute the second Lyapunov coefficient.

$\displaystyle c_{2}$
$\displaystyle=\frac{1}{2}g_{20}\left(\frac{h_{20}h_{12}}{2}+\frac{h_{22}}{2}+\frac{h_{30}h_{02}}{6}\right)$
$\displaystyle\hskip
56.9055pt+g_{11}\left(\frac{\bar{h}_{22}}{4}+\frac{h_{11}\bar{h}_{12}}{2}+\frac{h_{02}\bar{h}_{03}}{12}+\frac{h_{30}\bar{h}_{20}}{12}+\frac{h_{12}\bar{h}_{02}}{4}+\frac{h_{31}}{3}\right)$
$\displaystyle\hskip
56.9055pt+\frac{g_{02}}{2}\left(\frac{\bar{h}_{20}\bar{h}_{03}}{6}+\bar{h}_{11}\bar{h}_{12}+\frac{\bar{h}_{13}}{3}\right)+\frac{1}{6}g_{30}\left(3{h_{11}}^{2}+\frac{3h_{12}}{2}+\frac{3h_{20}h_{02}}{2}\right)$
$\displaystyle\hskip
56.9055pt+\frac{1}{2}g_{21}\left(\frac{|h_{20}|^{2}}{2}+2|h_{11}|^{2}+h_{20}h_{11}+\frac{|h_{02}|^{2}}{2}\right)$
$\displaystyle\hskip
56.9055pt+\frac{1}{2}g_{12}\left(\frac{h_{30}}{6}+\bar{h}_{11}^{2}+\bar{h}_{11}h_{20}+h_{11}\bar{h}_{02}+\bar{h}_{12}\right)$
$\displaystyle\hskip 56.9055pt\hskip
5.69046pt+\frac{1}{6}g_{03}\left(\frac{\bar{h}_{03}}{2}+3\bar{h}_{11}\bar{h}_{02}\right)+\frac{g_{40}h_{02}}{12}+\frac{g_{31}}{6}\left(3h_{11}+\frac{\bar{h}_{20}}{2}\right)$
$\displaystyle\hskip
56.9055pt+\frac{g_{22}}{4}\left(2\bar{h}_{11}+h_{20}\right)+\frac{g_{13}\bar{h}_{02}}{4}+\frac{g_{32}}{12},$

$\displaystyle h_{40}$
$\displaystyle=\frac{8}{\mu}\left[\frac{g_{20}}{2}\left(\frac{h_{20}^{2}}{4}+\frac{h_{30}}{3}\right)+g_{11}\left(\frac{\bar{h}_{03}}{6}+\frac{h_{20}\bar{h}_{02}}{4}\right)+\frac{g_{02}\bar{h}_{02}^{2}}{8}+\frac{g_{30}h_{20}}{4}+\frac{g_{21}\bar{h}_{02}}{4}+\frac{g_{40}}{24}\right],$
$\displaystyle h_{31}$
$\displaystyle=\frac{6}{2\mu+\bar{\mu}}\left[\frac{g_{20}}{2}\left(h_{21}+h_{20}h_{11}\right)+g_{11}\left(\frac{\bar{h}_{12}}{2}+\frac{h_{20}\bar{h}_{11}}{2}+\frac{h_{11}\bar{h}_{02}}{2}+\frac{h_{30}}{6}\right)\right.$
$\displaystyle\hskip
5.69046pt+\frac{g_{02}}{2}\left(\frac{\bar{h}_{03}}{3}+\bar{h}_{11}\bar{h}_{02}\right)+\left.\frac{g_{30}h_{11}}{2}+\frac{g_{21}}{2}\left(h_{20}+\bar{h}_{11}\right)+\frac{g_{12}\bar{h}_{02}}{2}+\frac{g_{31}}{6}-c_{1}h_{20}\right],$
$\displaystyle h_{22}$
$\displaystyle=\frac{4}{\mu+2\bar{\mu}}\left[\frac{g_{20}}{2}\left(h_{11}^{2}+h_{12}+\frac{h_{20}h_{02}}{2}\right)+g_{11}\left(\frac{|h_{20}|^{2}}{4}+|h_{11}|^{2}+\frac{|h_{02}|^{2}}{4}\right)+\frac{g_{21}}{2}\left(2h_{11}+\frac{\bar{h}_{20}}{2}\right)\right.$
$\displaystyle\hskip
5.69046pt+\left.\frac{g_{02}}{2}\left(\frac{\bar{h}_{20}\bar{h}_{02}}{2}+\bar{h}_{11}^{2}+\bar{h}_{12}\right)+\frac{g_{30}h_{02}}{4}+\frac{g_{12}}{2}\left(2\bar{h}_{11}+\frac{h_{20}}{2}\right)+\frac{g_{03}\bar{h}_{02}}{4}+\frac{g_{22}}{4}-2h_{11}\Re(c_{1})\right],$
$\displaystyle h_{13}$
$\displaystyle=\frac{2}{\bar{\mu}}\left[\frac{g_{20}}{2}\left(h_{11}h_{02}+\frac{h_{03}}{3}\right)+g_{11}\left(\frac{\bar{h}_{30}}{6}+\frac{h_{12}}{2}+\frac{h_{11}\bar{h}_{20}}{2}+\frac{\bar{h}_{11}h_{02}}{2}\right)+\frac{g_{02}}{2}\left(\bar{h}_{21}+\bar{h}_{11}\bar{h}_{20}\right)\right.$
$\displaystyle\hskip
5.69046pt+\left.\frac{g_{21}h_{02}}{2}+\frac{g_{12}\left(h_{11}+\bar{h}_{20}\right)}{2}+\frac{g_{03}\bar{h}_{11}}{2}+\frac{g_{13}}{6}-\bar{c}_{1}h_{02}\right],$
$\displaystyle h_{04}$
$\displaystyle=\frac{24}{4\bar{\mu}-\mu}\left[\frac{g_{20}h_{02}^{2}}{8}+g_{11}\left(\frac{h_{03}}{6}+\frac{h_{02}\bar{h}_{20}}{4}\right)+\frac{g_{02}}{2}\left(\frac{\bar{h}_{20}^{2}}{4}+\frac{\bar{h}_{30}}{3}\right)+\frac{g_{12}h_{02}}{4}+\frac{g_{03}\bar{h}_{20}}{4}+\frac{g_{04}}{24}\right].$

Assuming,

$\displaystyle
L_{1}(\delta=\delta_{H},\theta=\theta_{B})=\alpha(\delta=\delta_{H},\theta=\theta_{B})=0,$
(27)

the second Lyapunov coefficient at $(\delta=\delta_{H},\theta=\theta_{B})$ is
given by

$\displaystyle
L_{2}(\delta=\delta_{H},\theta=\theta_{B})=\frac{\Re\left(c_{2}\right)}{\beta}\bigg{|}_{(\delta=\delta_{H},\theta=\theta_{B})}.$
(28)

An expression for $L_{2}(\delta=\delta_{H},\theta=\theta_{B})$ as mentioned in
[22] in compact form is given by the following,

$\displaystyle 12L_{2}(\delta=\delta_{H},\theta=\theta_{H})$
$\displaystyle\hskip
28.45274pt=\frac{\Re(g_{32})}{\beta}+\frac{1}{\beta^{2}}\Im\left[g_{20}\bar{g}_{31}-g_{11}\left(4g_{31}+3\bar{g}_{22}\right)-\frac{1}{3}g_{02}\left(g_{40}+\bar{g}_{13}\right)-g_{30}g_{12}\right]$
$\displaystyle\hskip
28.45274pt+\frac{1}{\beta^{3}}\left\\{\Re\left[g_{20}\left(\bar{g}_{11}(3g_{12}-\overline{g_{30}})+g_{02}(\bar{g}_{12}-\frac{g_{30}}{3})+\frac{1}{3}\bar{g}_{02}g_{03}\right)\right.\right.$
$\displaystyle\hskip
28.45274pt\left.\left.+~{}g_{11}\left(\bar{g}_{02}(\frac{5}{3}\bar{g}_{30}+3g_{12})+\frac{1}{3}g_{02}\bar{g}_{03}-4g_{11}g_{30}\right)\right]+3\Im\left(g_{20}g_{11}\right)\Im(g_{21})\right\\}$
$\displaystyle\hskip
28.45274pt+\frac{1}{\beta^{4}}\left\\{\Im\left[g_{11}\bar{g}_{02}(\bar{g}_{20}^{2}-3\bar{g}_{20}g_{11}-4g_{11}^{2})\right]+\Im(g_{20}g_{11})\left[3\Re(g_{20}g_{11})-2|g_{02}|^{2}\right]\right\\},$

where all the expressions of $g_{kl}$ and $\beta$ are evaluated at the point
$(\delta=\delta_{H},\theta=\theta_{B})$, and $\Im(\cdot)$ stands for the
imaginary part of $(\cdot)$. After some computation, we have

$\displaystyle L_{2}(\sqrt{\epsilon})$
$\displaystyle=B\epsilon^{\frac{3}{2}},$ (29)

with

$\displaystyle B=$
$\displaystyle\frac{a_{01}^{4}b_{10}}{144\beta_{0}^{7}}\Big{(}18a_{11}^{3}a_{20}b_{10}^{2}+\left(27a_{20}b_{20}-2a_{30}b_{10}\right)a_{01}a_{11}^{2}b_{10}+\left(7a_{30}b_{10}b_{20}-27b_{20}^{2}a_{20}-10b_{10}^{2}a_{40}\right)a_{01}^{2}a_{11}$
(30)
$\displaystyle+2\left(10a_{40}b_{20}-3a_{50}b_{10}\right)a_{01}^{3}b_{10}+\left(48b_{20}a_{20}^{2}b_{11}+36a_{21}a_{30}b_{10}^{2}-30a_{20}^{2}b_{21}b_{10}\right)a_{01}^{2}\Big{)}$

We have $A=0$ for $\theta=\theta_{B}$ and consequently, for fixed
$0<\epsilon\ll 1$ the criticality of the second Lyapunov coefficient $L_{2}$
at the point $(\delta=\delta_{H},\theta=\theta_{B})$ is determined by $B$.

Thus, assuming the following conditions,

$\displaystyle B(\delta=\delta_{H},\theta=\theta_{B})\neq
0,\,\,\text{and}\,\,{\rm
Det}\left[\begin{array}[]{cc}\frac{\partial\alpha}{\partial\delta}&\frac{\partial\alpha}{\partial\theta}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr\frac{\partial
L_{1}}{\partial\delta}&\frac{\partial
L_{1}}{\partial\theta}\end{array}\right]_{\delta=\delta_{H},\
\theta=\theta_{B}}\neq 0,$ (31)

we observe that all the conditions of the Bautin bifurcation hold [22] and the
system undergoes a codimension-2 Bautin bifurcation at
$(\delta,\theta)=(\delta_{H},\theta_{B})$. We now summarize the above results
by the following theorem.

######  Theorem 4.2.

Assume $(\theta,\eta)\in\mathcal{R}$ and the conditions (10), (11), (27) and
(31) holds. Then $\exists$ $\epsilon_{0}>0$, $\delta_{0}>0$ and $\theta_{0}>0$
such that for $0<\epsilon<\epsilon_{0}$ and $|\delta-\delta_{*}|<\delta_{0}$,
$|\theta-\theta_{*}|<\theta_{0}$ the system (5) has an equilibrium point
$P_{2}$ in a neighbourhood of the fold point $P$ which converges to $P$ as
$(\epsilon,\delta,\theta)\to(0,\delta_{*},\theta_{*})$. The system (5)
undergoes a codimension-2 Bautin bifurcation at
$(\delta,\theta)=(\delta_{H},\theta_{B})$.

The Bautin bifurcation shows that for fixed $\theta$, $\theta<\theta_{B}$, the
system undergoes a subcritical singular Hopf bifurcation at
$\delta=\delta_{H}$. Subsequently, under the variation of parameter $\theta$
in a neighbourhood of $\theta=\theta_{B}$, the singular Hopf bifurcation is
then accompanied by a saddle-node bifurcation of canard cycles showing the
appearance of unstable and stable canard cycles, coalescent and disappearance
of the cycles.

##  5 Canard Explosion and Relaxation Oscillation

We have shown in the previous section that the model system (5) undergoes a
singular Hopf bifurcation at $\delta=\delta_{H}$ assuming the conditions (10),
(11) and $A\neq 0$. Now, the small $\mathcal{O}(\sqrt{\epsilon})$ amplitude
family of limit cycles generated through the singular Hopf bifurcation at
$\delta=\delta_{H}$ changes abruptly within an exponentially small parameter
interval of the bifurcation parameter $\delta$. This phenomenon in literature
is known as “canard explosion” which connects the small amplitude limit cycles
generated through the singular Hopf bifurcation to large relaxation cycles in
a continuous fashion. Hence, based on the theorems 3.2 and 3.3 of [20] and
theorem 8.4.3 of [21], we have the following results for the existence of
maximal canard and canard cycles for the system (5).

######  Theorem 5.1.

Let $(\theta,\eta)\in\mathcal{R}$ and the conditions (10) and (11) hold. Then
$\exists$ $\epsilon_{0}>0$ and $\delta_{0}>0$ such that for
$0<\epsilon<\epsilon_{0}$ and $|\delta-\delta_{*}|<\delta_{0}$, the system (5)
has an equilibrium point $P_{2}$ in a neighbourhood of the fold point $P$
which converges to $P$ as $(\epsilon,\delta)\to(0,\delta_{*})$. The system (5)
admits a maximal canard at

$\displaystyle\delta_{C}=\delta_{*}-\left(\frac{b_{10}}{2a_{20}^{3}v_{m}(u_{m}^{2}+\eta)}A\right)\epsilon+\mathcal{O}(\epsilon^{2}).$
(32)

######  Theorem 5.2.

Let $(\theta,\eta)\in\mathcal{R}$, the conditions (10) and (11) hold and
$\delta_{C}$ is the maximal canard value. Then there exists a smooth
parameterized family of canard cycles growing from
$\mathcal{O}(\sqrt{\epsilon})$ amplitude to a relaxation oscillation within an
exponentially small parameter interval. However, if $A<0,$ then the singular
Hopf bifurcation is supercritical and the family of canard cycles uniformly
close to the canard point is stable. Whereas, if $A>0,$ then the family of
canard cycles is unstable and there exists a unique parameter value
$\delta_{\mathrm{SNL}}$ where the cycles undergo a saddle-node bifurcation of
limit cycles.

Canard explosion and relaxation oscillation are both considered to be global
phenomena that can be observed in systems with multiple timescales, as they
characterize transitions from a small amplitude to a large amplitude canard
cycle. We next prove the existence of relaxation oscillation in the system
(5).

######  Theorem 5.3.

Assume that $(\eta,\theta)\in\mathcal{R}$ and $(u_{m},v_{m})$ and
$(u_{M},v_{M})$ be the two fold points of the critical manifold $M_{20}$ where
$u_{m}<u_{M}.$ Let $u_{m}<u_{*}<u_{M},$ then, for $\epsilon>0$ sufficiently
small, the system has a unique relaxation oscillation $\gamma_{\epsilon}$.

###### Proof.

Let $(u_{l},v_{M})$ be the point of intersection of the horizontal line
$v=v_{M}$ with the critical manifold $S^{l}_{0}$ and $(u_{r},v_{m})$ be the
point of intersection of the line $v=v_{m}$ with the critical manifold
$S^{r}_{0}.$ We then define the singular trajectory $\gamma_{0}$ as a union of
two alternative slow and fast flows. Let us define the horizontal segments
joining $(u_{m},v_{m})$ to $(u_{r},v_{m})$ as $l_{1}$ and $(u_{M},v_{M})$ to
$(u_{l},v_{M})$ as $l_{2}.$ Then the fast flows are along the line segments
$l_{1}$ and $l_{2}.$ Let us define the slow flow along the sub-manifold
$S^{l}_{0}$ from $(u_{l},v_{M})$ to $(u_{m},v_{m})$ as $c_{l}$ and along
$S^{r}_{0}$ from $(u_{r},v_{m})$ to $(u_{M},v_{M})$ as $c_{r}.$ We then define
the singular trajectory as

$\gamma_{0}=l_{1}\cup c_{r}\cup l_{2}\cup c_{l}.$

Let the unique coexistence equilibrium $E_{*}=(u_{*},v_{*})$ lies on the
normally hyperbolic repelling branch $S^{m}_{0}.$ From the slow-fast dynamics
of the system (4), the equilibrium $E_{*}$ is unstable. The two fold points
$(u_{m},v_{m})$ and $(u_{M},v_{M})$ are the jump points of the system. We now
consider a small horizontal section $\Delta$ transversal to $S^{l}_{0},$ where

$\Delta=\\{(u,v_{0}):u\in[a,a+\rho],\,\,0<a<u_{m},\,\,v_{m}<v_{0}<v_{M},\,\,\rho>\epsilon>0\\},$

and track two trajectories $\xi_{\epsilon}^{1},\,\xi_{\epsilon}^{2}$ starting
on $\Delta.$ From Fenichel’s theorem [6], for $\epsilon>0,$ the normally
hyperbolic sub-manifolds $S^{l}_{0}$ and $S^{r}_{0}$ perturbs to
$S^{l}_{\epsilon}$ and $S^{r}_{\epsilon}$ respectively. Therefore, for
$\epsilon>0$ the trajectories starting at $\Delta$ get attracted towards
$S^{l}_{\epsilon}$ with an exponential rate. From the vicinity of the fold
point $(u_{m},v_{m}),$ it then jumps to the other attracting manifold
$S^{r}_{\epsilon},$ and follows until the vicinity of another jump point
$(u_{M},v_{M}).$ The trajectories then jump and get exponentially attracted
towards $S^{l}_{\epsilon}$ and return to $\Delta.$ We then define a return map
$\Pi:\Delta\rightarrow\Delta.$ From [19], it shows that the map $\Pi$ is a
contraction map with an exponential contracting rate. Therefore, from the
contraction mapping theorem, there exists a unique fixed point. This fixed
point is attracting, which gives rise to a unique limit cycle
$\gamma_{\epsilon}.$ By [6, 19] we can conclude that
$\gamma_{\epsilon}\rightarrow\gamma_{0}$ as $\epsilon\rightarrow 0.$ ∎

Figure 2: A schematic diagram of the singular trajectory $\gamma_{0}.$

######  Theorem 5.4.

Assume that $(\eta,\theta)\in\mathcal{R}$ and $P(u_{m},v_{m}),Q(u_{M},v_{M})$
be the two fold points of the critical manifold $M_{20}$. Let $u_{*}>u_{M},$
then, for $\epsilon>0$ sufficiently small, the coexistence equilibrium
$E_{*}(u_{*},v_{*})$ is globally stable.

###### Proof.

Let the predator nullcline of the system (5)
$u=\sqrt{\frac{\delta\eta}{1-\delta}}$ coincides with the sub-manifold
$S^{r}_{0}$ of the critical manifold $M_{20}$ at the coexistence equilibrium
$E_{*}.$ From the geometry of the $S^{r}_{0},$ we can infer that $E_{*}$ is
locally stable. Let us take a vertical line $u=u_{M}$ which divides the
interior of the first quadrant $\mathbb{R}^{2}_{+}$ into two regions $D_{1}$
and $D_{2}$, where

$\displaystyle D_{1}=\\{(u,v):0<u\leq
u_{M},v>0\\}\,\,\text{and}\,\,D_{2}=\\{(u,v):u>u_{M},v>0\\}.$

In the region $D_{2}$, we define the Dulac function $D:D_{2}\to\mathbb{R}$ by
$D(u,v)=\frac{1}{u^{2}v}$. It then follows from the nature of the critical
manifold $M_{20}$ that

$\displaystyle\frac{\partial(fD)}{\partial u}+\frac{\partial(\epsilon
gD)}{\partial v}=\phi^{\prime}(u)<0,\,\,\forall\,\,u,v\in D_{2}.$

Hence, by Dulac criterion, there exists no periodic orbit in $D_{2}$. Thus,
any trajectory starting in the region $D_{2}$ converges to $E_{*}.$ For global
stability we now track the trajectories starting in $D_{1}.$ We consider two
trajectories: $\xi_{1}$ starting above the critical manifold $M_{20}$ and
$\xi_{2}$ starting below $M_{20}.$ Since $S^{l}_{0}$ is a normally hyperbolic
attracting sub-manifold, $\xi_{1}$ gets attracted to the vicinity of
$S^{l}_{0}$ till it encounters the fold point $P.$ Note that both the fold
points $P,Q$ are jump points and the fast flow moves away from the critical
manifold $M_{20}$. It then gets attracted to the vicinity of $S^{r}_{0}$ which
lies in the region $D_{2}.$ Thus it converges to $E_{*}.$ However, the
trajectory $\xi_{2}$ starting below $M_{20}$ and in between two fold points
$P$ and $Q$ gets attracted in the region $D_{1}.$ Thus, it converges to
$E_{*}.$ Therefore any trajectory starting in $D_{1}$ or $D_{2}$ will approach
to the equilibrium $E_{*}$ as $t\to\infty$. ∎

######  Theorem 5.5.

Assume that $(\eta,\theta)\in\mathcal{R}$ and the quantity $A$ given in (23)
is positive. Let $P(u_{m},v_{m})$, $Q(u_{M},v_{M})$ be the two fold points of
the critical manifold $M_{20},$ and $\delta_{\mathrm{SNL}}<\delta<\delta_{H}$
then for $0<\epsilon\ll 1$ sufficiently small, the system has two relaxation
oscillations, where the outer cycle is a stable relaxation oscillation and the
inner cycle is an unstable relaxation oscillation.

###### Proof.

The proof of this theorem is a direct consequence of Theorem (5.2). The
stability of the canard cycles can be computed numerically by using the
constant $A$ from (23). ∎

##  6 Numerical Illustration

For $\epsilon$ sufficiently small i.e., $0<\epsilon\ll 1$, the asymptotic
expansion of the first Lyapunov coefficient $L_{1}$ in the blow-up coordinates
(13) is computed in eqn. (22). It has been shown that the sign of the leading
order term $A$ in the expansion of $L_{1}$ determines the criticality of the
singular Hopf bifurcation for $\epsilon\to 0$. The family of canard cycles
which emerge due to singular Hopf bifurcation grow from
$\mathcal{O}(\epsilon)$ amplitude to large amplitude relaxation oscillation
within an exponential small change of the bifurcation parameter. Now, the
criticality of the singular Hopf bifurcation will be changed for $\epsilon>0$
sufficiently small, if $A$ changes its sign and hence, the singular Hopf
bifurcation will be degenerate if $A=0$.

In such a case, the system undergoes a codimension-2 Bautin bifurcation and
the Bautin bifurcation threshold is given by
$(\delta,\theta)=(\delta_{H},\theta_{B})$. Now, keeping all other parameters
fixed, in a small neighborhood of the Bautin threshold point
$(\delta,\theta)=(\delta_{H},\theta_{B})$ in the $\delta-\theta$ plane, the
existence of a single canard cycle or two canard cycles (unstable and stable)
and the coalescent and disappearance of the cycles through a saddle-node
bifurcation of cycles can be observed. We call this singular Bautin
bifurcation in the sense that the associated Hopf bifurcation is singular; and
for $\epsilon>0$ sufficiently small, the inner unstable canard cycle is the
small canard cycle but, the outer stable canard cycle grows to a big
relaxation oscillation.

To numerically validate the results obtained We fix the parameters
$\eta,\,\theta$ throughout our simulations such that P4:
$\theta<\min\\{1,\eta\\}.$ We choose the hypothetical values from the above
domain as $\theta=0.05$ and $\eta=0.176.$ For this set of parameter values, we
obtain two real roots of the quartic equation (9) which are the $u$ components
of the fold points of the critical manifold. Thus, the fold points of $M_{20}$
are

$P(u_{m},v_{m})=(0.2375,0.2145),\,\text{and}\,\,Q(u_{M},v_{M})=(0.5359,0.235).$

The parameter $\delta$ is chosen as the bifurcation parameter throughout the
text to track the different singular bifurcations. Therefore, depending on the
value of $\delta$ the coexistence equilibrium $E_{*}$ either lies on
$S^{l}_{0},$ $S^{m}_{0}$ or on $S^{r}_{0}$ part of critical manifold,
$M_{20}.$ If $E_{*}$ coincides with either of the fold points, for some value
of $\delta,$ then those thresholds determine a critical singular threshold.
That is, for $\delta_{*}=0.2426879409,$ the predator nullcline intersects the
critical manifold $M_{20}$ at $P(u_{m},v_{m})=(0.2375,0.2145).$ The Jacobian
matrix evaluated at $Q$ has purely imaginary complex conjugate eigenvalues.
Thus, from the condition of Hopf bifurcation, the system (5) undergoes a
singular Hopf bifurcation, and the quantity $A$ given by the equation (23) is

$A\approx 2.796\times 10^{-7}.$

This implies that the Hopf bifurcation is degenerate and subcritical.
Therefore, the value of $\delta$ gives the threshold for singular Bautin
bifurcation. As a result, the canard cycle originating from singular Hopf
bifurcation is unstable, and the maximal canard curve is given by the equation
5.1

$\displaystyle\delta_{C}$
$\displaystyle\approx\delta_{*}-\mathcal{O}(10^{-9}).$ (33)

Since the transition from small amplitude unstable canard cycle to unstable
maximal canard occurs in an extremely small parameter interval, it is
difficult to catch this phenomenon through numerical simulation. However, we
check for the second Lyapunov coefficient $L_{2}$ to confirm that this is a
point of singular Bautin bifurcation. From the expression (29) and (30) we
obtain the leading order coefficient of $L_{2}$ as

$B=-0.004.$

(a)

(b)

Figure 3: (a) The bifurcation diagram of the system for
$\eta=0.176,\,\theta=0.05,\,\epsilon=0.005.$ The zoomed-in diagram of the
dotted rectangular box is shown in (b). The red broken and continuous line
corresponds to the unstable equilibrium and unstable canard cycle,
respectively. The blue continuous line represents the maxima and minima of the
stable component (equilibrium and limit cycle). The black dot is the Singular
Hopf point and the broken line is the threshold for saddle-node bifurcation of
limit cycles $(\delta_{\mathrm{SNL}}).$

The parameter $\delta$ is again varied such that the predator nullcline
coincides with the prey nullcline at the point $Q(u_{M},v_{M})$ when
$\delta=0.62.$ At this point the system undergoes a singular Hopf bifurcation
which is supercritical, since

$A=-0.1055<0.$

Thus, a stable canard cycle originates from this point and exists for
$0.24268<\delta<0.62.$ Since the direction of the Hopf bifurcation is reversed
in the maximum and minimum point of the critical manifold so is the direction
of the canard explosion. Therefore, in this case, we observe the transition
from small amplitude canard cycle to relaxation oscillation whenever
$\delta<\delta_{*}$ and it continues in this direction. From theorem (5.2), we
conclude that the family of canard cycles that emerges from this point are
stable. Therefore, theorem (5.4) states that there exists a unique relaxation
oscillation. In Fig. (3), we plot the maximum and minimum components of the
variable $u$ with respect to the parameter $\delta.$ The bifurcation diagram
shows two Hopf bifurcation thresholds (cf. Fig.(3)a). The zoomed-in diagram of
the extreme left dotted box is shown in Fig. (3)(b). It shows that in an
extremely small neighborhood of the Bautin bifurcation,
$\delta_{\mathrm{SNL}}<\delta<\delta_{*},$ the small unstable canard cycles
(red) and the stable canard cycles (blue) exists and they coalesce at the
saddle-node bifurcation of the limit cycle. The threshold for $\delta$ is
denoted as $\delta_{\mathrm{SNL}}$ and marked by a black broken line.

##  7 Discussion

In this paper, we have considered a singularly perturbed planar predator-prey
system with Holling type III functional response, where the prey population
growth is affected by the weak Allee effect, and the prey reproduces at a
significantly higher rate than the predator. The chosen model in the slow-fast
framework admits explicit components for coexistence equilibrium, which
destabilizes through supercritical and subcritical Hopf bifurcation for
specific parameters. The Bautin bifurcation point separates the subcritical,
and supercritical branches of the Hopf bifurcation curve. The stable and
unstable limit cycles bifurcated through supercritical and subcritical Hopf
bifurcation disappear through a saddle-node bifurcation of limit cycles.

We provide a thorough mathematical analysis of the system by using tools
namely geometric singular perturbation theory, normal form theory of slow-fast
systems, and blow-up technique to study a wide range of rich and complex
nonlinear dynamics, such as singular Hopf bifurcations, singular Bautin
bifurcations, canard cycles, canard explosion in the neighbourhood of
subcritical and supercritical Hopf bifurcation. We have identified the
occurrence of at most two relaxation oscillations via canard cycles known as
the “boom and bust cycle” [27, 28]. Important ecological ramifications
stemming from the existence of relaxation oscillation and the onset of canard
explosion are also an important part of the discussion. Using the entry-exit
function, we have investigated that relaxation oscillations indeed occur in
this slow-fast system, which suggests the possibility of the coexistence of
both populations with a predictable pattern of rapid population growth and
contractions over a significantly short time interval. On the other hand, the
canard explosion is a surprising and fascinating event that happens in a tiny
exponential region of the parameter $\delta$. From a biological perspective,
this canard explosion can be seen as an indicator of an impending regime shift
brought on by an exponentially modest change in parameter $\delta$.

The main contribution of this work is two–fold:

  1. 1.

We have derived the analytic form of the second Lyapunov coefficient and
presented a thorough study of the singular Bautin bifurcation for a system
under a slow-fast framework for the first time in the literature (see [19, 20,
21]). The normal form for singular Bautin bifurcation has been derived by
explicitly finding the locally invertible parameter-dependent transformations.
To the best of our knowledge, no literature provides an explicit derivation of
locally invertible parameter-dependent transformations.

  2. 2.

We have derived the expansions of the first and second Lyapunov coefficients.
It has been observed that the leading order coefficient of $L_{1}$, i.e., the
coefficient of $\sqrt{\epsilon}$ involved with the expression of $L_{1}$,
determines the criticality of the singular Hopf bifurcation, whereas the
leading order coefficient of $L_{2}$, i.e., the coefficient of
$\epsilon^{\frac{3}{2}}$ in $L_{2}$, determines the criticality of the
singular Bautin bifurcation.

It is important to mention here that this approach would be applicable to
investigate Bautin bifurcation for a wide range of slow-fast systems,
including but not limited to biological, physical, ecological, and medical
fields. Further, how this approach can be extended for higher dimensional
slow-fast systems is a non-trivial open problem.

## References

  * [1] V. I. Arnold, Dynamical systems V: bifurcation theory and catastrophe theory, vol. 5 (Encyclopaedia of Mathematical Sciences), Springer-Verlag, Berlin, 1994.
  * [2] A. Atabaigi, Canard explosion, homoclinic and heteroclinic orbits in singularly perturbed generalist predator–prey systems, International Journal of Biomathematics, 14 (2021), p. 2150003.
  * [3] D. Bai, Y. Kang, S. Ruan, and L. Wang, Dynamics of an intraguild predation food web model with strong Allee effect in the basal prey, Nonlinear Analysis: Real World Applications, 58 (2021), p. 103206.
  * [4] F. Courchamp, L. Berec, and J. Gascoigne, Allee effects in ecology and conservation, OUP Oxford, 2008.
  * [5] F. Dumortier and R. Roussarie, Canard cycles and center manifolds, vol. 577, American Mathematical Soc., 1996.
  * [6] N. Fenichel, Geometric singular perturbation theory for ordinary differential equations, Journal of Differential Equations, 31 (1979), pp. 53–98.
  * [7] R. FitzHugh, Mathematical models of threshold phenomena in the nerve membrane, The Bulletin of Mathematical Biophysics, 17 (1955), pp. 257–278.
  * [8] B. Garrett and R. Gian-Carlo, Ordinary Differential Equations, Wiley Boston, 1982.
  * [9] W. M. Getz, A hypothesis regarding the abruptness of density dependence and the growth rate of populations, Ecology, 77 (1996), pp. 2014–2026.
  * [10] X. Han and Q. Bi, Slow passage through canard explosion and mixed-mode oscillations in the forced van der pol’s equation, Nonlinear Dynamics, 68 (2012), pp. 275–283.
  * [11] G. Hek, Geometric singular perturbation theory in biological practice, Journal of Mathematical Biology, 60 (2010), pp. 347–386.
  * [12] S.-B. Hsu, T.-W. Hwang, and Y. Kuang, Global analysis of the Michaelis–Menten-type ratio-dependent predator-prey system, Journal of Mathematical Biology, 42 (2001), pp. 489–506.
  * [13] S.-B. Hsu, T.-W. Hwang, and Y. Kuang, Global dynamics of a predator-prey model with Hassell-Varley type functional response, Discrete Contin. Dyn. Syst. Ser. B, 10 (2008), pp. 857–871.
  * [14] T.-H. Hsu and S. Ruan, Relaxation oscillations and the entry-exit function in multidimensional slow-fast systems, SIAM Journal on Mathematical Analysis, 53 (2021), pp. 3717–3758.
  * [15] J. Huang, S. Ruan, and J. Song, Bifurcations in a predator–prey system of Leslie type with generalized Holling type III functional response, Journal of Differential Equations, 257 (2014), pp. 1721–1752.
  * [16] B. Kooi and J. Poggiale, Modelling, singular perturbation and bifurcation analyses of bitrophic food chains, Mathematical Biosciences, 301 (2018), pp. 93–110.
  * [17] K. U. Kristiansen, Geometric singular perturbation analysis of a dynamical target mediated drug disposition model, Journal of Mathematical Biology, 79 (2019), pp. 187–222.
  * [18] M. Krupa, N. Popović, and N. Kopell, Mixed-mode oscillations in three time-scale systems: a prototypical example, SIAM Journal on Applied Dynamical Systems, 7 (2008), pp. 361–420.
  * [19] M. Krupa and P. Szmolyan, Extending geometric singular perturbation theory to nonhyperbolic points—fold and canard points in two dimensions, SIAM Journal on Mathematical Analysis, 33 (2001), pp. 286–314.
  * [20] M. Krupa and P. Szmolyan, Relaxation oscillation and canard explosion, Journal of Differential Equations, 174 (2001), pp. 312–368.
  * [21] C. Kuehn, Multiple time scale dynamics, vol. 191 of Applied Mathematical Sciences, Springer Cham, 2015.
  * [22] Y. A. Kuznetsov, Elements of applied bifurcation theory, vol. 112 of Applied Mathematical Sciences, Springer, 1998.
  * [23] C. Li and H. Zhu, Canard cycles for predator–prey systems with Holling types of functional response, Journal of Differential Equations, 254 (2013), pp. 879–910.
  * [24] J. Li, T. Quan, and W. Zhang, Bifurcation and number of subharmonic solutions of a 4D non-autonomous slow–fast system and its application, Nonlinear Dynamics, 92 (2018), pp. 721–739.
  * [25] M. Lu and J. Huang, Global analysis in Bazykin’s model with Holling II functional response and predator competition, Journal of Differential Equations, 280 (2021), pp. 99–138.
  * [26] P. J. Pal, T. Saha, M. Sen, and M. Banerjee, A delayed predator–prey model with strong Allee effect in prey population growth, Nonlinear Dynamics, 68 (2012), pp. 23–42.
  * [27] S. Rinaldi and S. Muratori, Slow-fast limit cycles in predator-prey models, Ecological Modelling, 61 (1992), pp. 287–308.
  * [28] S. Sadhu, Analysis of the onset of a regime shift and detecting early warning signs of major population changes in a two-trophic three-species predator-prey model with long-term transients, Journal of Mathematical Biology, 85 (2022), pp. 1–33.
  * [29] G. A. Van Voorn, L. Hemerik, M. P. Boer, and B. W. Kooi, Heteroclinic orbits indicate overexploitation in predator–prey systems with a strong Allee effect, Mathematical Biosciences, 209 (2007), pp. 451–469.
  * [30] J. Wang, J. Shi, and J. Wei, Predator–prey system with strong Allee effect in prey, Journal of Mathematical Biology, 62 (2011), pp. 291–331.
  * [31] Y. Xia, Z. Zhang, and Q. Bi, Relaxation oscillations and the mechanism in a periodically excited vector field with pitchfork–Hopf bifurcation, Nonlinear Dynamics, (Online) (2020).
  * [32] L. Yaru and L. Shenquan, Canard-induced mixed-mode oscillations and bifurcation analysis in a reduced 3D pyramidal cell model, Nonlinear Dynamics, (2020), pp. 1–37.