# Stability and Bifurcation Analysis of Two-Term Fractional Differential
Equation with Delay

Sachin Bhalekar<EMAIL_ADDRESS>Deepa Gupta<EMAIL_ADDRESS>School of Mathematics and Statistics, University of Hyderabad, Hyderabad,
Telangana, India- 500046

###### Abstract

This manuscript deals with the stability and bifurcation analysis of the
equation $D^{2\alpha}x(t)+cD^{\alpha}x(t)=ax(t)+bx(t-\tau)$, where
$0<\alpha<1$ and $\tau>0$. We sketch the boundaries of various stability
regions in the parameter plane under different conditions on $\alpha$ and $b$.
First, we provide the stability analysis of this equation with $\tau=0$.
Change in the stability of the delayed counterpart is possible only when the
characteristic roots cross the imaginary axis. This leads to various delay-
independent as well as delay-dependent stability results. The stability
regions are bifurcated on the basis of the following behaviors with respect to
the delay $\tau$ viz. stable region for all $\tau>0$, unstable region, single
stable region, stability switch, and instability switch.

> Fractional derivatives and the delay are the popular tools used by applied
> scientists to model the accurate behaviors of natural systems. As the
> fractional derivative is a nonlocal operator, one has to provide the
> complete history of the state during its computation. On the other hand, the
> delay uses the short history of the state. These equations have proved
> useful in modeling the systems in applied sciences and engineering. Though
> these equations look very similar to the classical ordinary differential
> equations, the analysis is, however, not simple. The characteristic
> equations of these equations are of a transcendental nature and admit
> infinitely many roots, in general. Therefore, the simple conditions in terms
> of the parameters, as provided in this manuscript, are very useful to the
> researchers. The complex dynamics of these systems include not only stable
> and unstable behaviors but the delay-dependent features such as a single
> stable region provided by the delay interval and stability/instability
> switches.

##  I Introduction

Differential equations (DE) can be treated as the heart of mathematical
analysis. These equations emerge as models in natural systems. As the
derivative represents the rate of change of a quantity, the DE models are
widely used by Scientists, Engineers, Economists, and so on
simmons2016differential ; khalil2002control ; strogatz2018nonlinear . The
second-order ordinary differential equations are found in classic examples
such as Newton’s law of gravitation, simple harmonic motion (damped undamped
and forced oscillations), Kepler’s laws of planetary motion, and LRC circuit
in electronics simmons2016differential . B van der Pol van1926lxxxviii  used
these equations to model the reaction oscillators.

Though researchers widely use the ODE models, they are unsuitable for modeling
the memory properties in natural systems. Volterra volterra  proposed the
integral equations containing the nonlocal integral operator to model the
biological systems. Furthermore, the delay in the model also has the ability
to include the history of the state in the model lakshmanan2011dynamics ;
hale2012functional . If the rate of change of the present state depends on the
state at past times then the resulting equation is called the delay
differential equation (DDE). Smith smith2011introduction  discussed the
applications of the DDEs in life science. The stability and oscillation theory
of the DDEs arising in the population dynamics is developed by Gopalsamy
gopalsamy2013stability . Various applications of these equations in fluid
dynamics, economics, mechanical engineering, life science, chemistry, and
physics are presented by Erneux in the book erneux2009applied .

The constant delay in the equation is adequate for modeling the “short-time
memory". However, we need nonlocal operators such as “fractional order
derivative" to include the “long-term memory" in the model. The order of the
derivative in the “fractional order differential equation (FDE)" is “non-
integer". As one can expect, the generalization to fractional derivative (FD)
can be done in many ways podlubny1998fractional ; kilbas2006theory . However,
we stick to the definition provided by Caputo because it is more suitable to
real-life problems.

Thus, there are improvements in the model involving both “fractional order
derivative" and the “delay" so as to describe the “real behavior" of the
natural systems. Bhalekar, Daftardar-Gejji, and coworkers propose numerical
methods to solve fractional-order delay differential equations (FDDE) in
bhalekar2011predictor ; daftardar2015solving . Stability analysis of the
scalar FDDE is proposed by Bhalekar in bhalekar2016stability  where the
regions of the stability are provided in the parameter plane. Stability,
bifurcations and chaos in FDDEs are investigated in bhalekar2011fractional ;
bhalekar2010fractional ; baleanu2015chaos ; daftardar2012dynamics ;
bhalekar2012dynamical ; bhalekar2013stability ; bhalekar2022stability ;
gu2017artificial ; shi2022chaos . The detailed bifurcation analysis of scalar
FDDE is presented by Bhalekar and Gupta in bhalekar2023can .

The DDEs involving the second-order derivative are useful in modeling a wide
range of systems. Second-order DDE with distributed delay appears as a model
of hereditary dynamics (see volterra , pp. 191, eq. (3)). Milton and Longtin
milton1990evaluation  used these equations to study the human pupil cycle.
Campbell et al. campbell1995complex  proved the existence of limit cycles,
two-tori, and multi-stability in the damped harmonic oscillator with delayed
feedback. Hopf bifurcations in the delayed Duffing oscillator are studied in
erneux2023short . Hybrid bistable device vallee1987second  in the optics can
also be modeled by using these equations.

This discussion motivated us to analyze the stability of the “two-term
fractional order delay differential equation," which is the generalization of
second-order DDE. We provide the complete bifurcation analysis of this
equation under various conditions on parameters and fractional order. We
observed various stability behaviors of this system.

The paper is organized as follows: Section (II) deals with some basic
definitions and Theorems. Section (III) provides stability analysis of non-
delayed case of equation (1). In Section (IV) we give the sufficient condition
under which our equation (8) have the non-existence of critical values of
delay. Section (V) deals with the conditions on the parameter for the
existence of positive root. Some conditions under which the characteristic
equation has roots with positive real part are described in Section (VI). The
open problems are given in Section (VII). Section (VIII) will provide the
region in the parametric plane where the stability depends on delay parameter.
In section (IX), we have some examples which validate our results. Section (X)
summarizes the results given in the paper.

##  II Preliminaries

In this section, we provide some basic definitions described in the literature
podlubny1998fractional ; lakshmanan2011dynamics ; kilbas2006theory ;
bhalekar2013stability ; diethelm2002analysis .

######  Definition II.1.

$\mathcal{L}_{p}[a,b]:=\\{f:[a,b]\rightarrow\mathbb{R};\textit{ f is
measurable on $[a,b]$}\textit{ and }\int_{a}^{b}|f(x)|^{p}dx<\infty\\}$.

######  Definition II.2 (Fractional Integral).

For any $f\in\mathcal{L}_{1}(0,b)$ the Riemann-Liouville fractional integral
of order $\upmu>0$, is given by

$\textit{I}^{\upmu}f(t)=\dfrac{1}{\Gamma(\upmu)}\int_{0}^{t}(t-\tau)^{\upmu-1}f(\tau)d\tau,\quad
0<t<b.$

######  Definition II.3 (Caputo Fractional Derivative).

For $f^{m}\in\mathcal{L}_{1}(0,b)$, $0<t<b$ and $m-1<\upmu\leq m$,
$m\in\mathbb{N}$, the Caputo fractional derivative of function $f$ of order
$\upmu$ is defined by,

$\textit{D}^{\upmu}f(t)=\begin{cases}\frac{d^{m}}{dt^{m}}f(t),\textit{ if
}\quad\upmu=m\\\ \textit{I}^{m-\upmu}\dfrac{d^{m}f(t)}{dt^{m}},\textit{ if
}\quad m-1<\upmu<m.\end{cases}$

Note that for $m-1<\upmu\leq m$, $m\in\mathbb{N},$

$\textit{I}^{\upmu}\textit{D}^{\upmu}f(t)=f(t)-\sum_{k=0}^{m-1}\dfrac{d^{k}f(0)}{dt^{k}}\dfrac{t^{k}}{k!}.$

######  Definition II.4.

vcermak2015stability  The two term FDE

$D^{\alpha}x(t)+cD^{2\alpha}x(t)=a_{1}x(t),\quad 0<\alpha<1$ (1)

where $a_{1}\in\mathbb{R}$ is said to be
(a) stable if all its solutions are bounded as $t\rightarrow\infty$.
(b) asymptotically stable if all its solutions tend to zero as
$t\rightarrow\infty$;
(c) $t^{-\gamma}$ asymptotically stable if there is a real scalar $\gamma>0$
such that any solution of equation (1) tends to zero like $O(t^{-\gamma})$ as
$t\rightarrow\infty$.

######  Theorem II.1.

vcermak2015stability  The FDE (1) is asymptotically stable if and only if all
the zeros $s$ of the polynomial

$p(s)=s^{2}+\frac{1}{c}s-\frac{a_{1}}{c}$ (2)

satify

$|arg(s)|>\alpha\pi/2.$ (3)

More precisely, the condition (3) is necessary and sufficient for the
$t^{-\gamma}$ asymptotic stability of equation (1), where

$\displaystyle\gamma=\begin{cases}&\alpha,\textit{ if }0<\alpha\leq 1/2\\\
&2\alpha-1,\textit{ if }1/2<\alpha<1.\\\ \end{cases}$

##  III Stability region for equation (1) in the $a_{1}c$-plane

Roots of equation (2) are

$s_{1}=\frac{-1+\sqrt{1+4a_{1}c}}{2c},\quad
s_{2}=\frac{-1-\sqrt{1+4a_{1}c}}{2c}.$

Case 1 If $a_{1}c>-1/4$ then $s_{1},s_{2}\in\mathbb{R}$.
The condition for stability of (3) becomes $s_{1}<0$ and $s_{2}<0$. If $c<0$
then $s_{2}>0$ and the system is unstable. Therefore we assume that $c>0$. In
this case, $s_{2}<0$. Further, $s_{1}<0$ if $\frac{-1}{4c}<a_{1}<0.$ Thus, the
stability condition (3) becomes

$c>0,\quad\frac{-1}{4c}<a_{1}<0.$ (4)

Case 2 If $a_{1}c<\frac{-1}{4}$ then the roots of equation (2) are

$\displaystyle s_{1}=\frac{-1+i\sqrt{-1-4a_{1}c}}{2c},\quad
s_{2}=\frac{-1-i\sqrt{-1-4a_{1}c}}{2c}.$

Therefore, $|arg(s_{1})|=|arg(s_{2})|$.
Subcase 2.1 $c>0$.
So, $arg(s_{1})=\pi-\arctan(\sqrt{-1-4a_{1}c}).$

Since,

$\displaystyle\frac{-\pi}{2}\leq-\arctan(\sqrt{-1-4a_{1}c})\leq 0,$
$\displaystyle\pi-\frac{\pi}{2}\leq\pi-\arctan(\sqrt{-1-4a_{1}c})\leq\pi.$

Therefore, we get $arg(s_{1})\geq\frac{\pi}{2}>\frac{\alpha\pi}{2}$ for any
$\alpha\in(0,1)$. Thus, the stability condition (3) reduces to

$c>0\textit{ and }a_{1}c<\frac{-1}{4}.$ (5)

Subcase 2.2 $c<0.$

We have, $|arg(s_{1})|=|\arctan(\sqrt{-1-4a_{1}c})|$.

If $a_{1}c<\dfrac{-\tan^{2}(\alpha\pi/2)-1}{4}$ then, we get

$4a_{1}c+1<-\tan^{2}(\alpha\pi/2)$.

$\Rightarrow\sqrt{-(4a_{1}c+1)}>\tan(\alpha\pi/2)$.

So, $\arctan(\sqrt{-(4a_{1}c+1)})>\alpha\pi/2$.

$\Rightarrow|arg(s_{1})|>\alpha\pi/2.$

Hence, by Theorem (II.1), equation (1) is stable if

$c<0\textit{ and }a_{1}c<\frac{-\tan^{2}(\alpha\pi/2)-1}{4}.$ (6)

Using the conditions (4), (5) and (6), we sketch the stability region of
equation (1) in $a_{1}c$ plane as shown in Figure(1).

Note that, in the first and third quadrant both the roots $s_{1}$ and $s_{2}$
are real and one of the root is positive. So, the first and third quadrants
are unstable. The second quadrant is divided in two parts by the curve
$\Gamma_{1}:a_{1}c=-1/4$. on the right side of $\Gamma_{1}$ the roots $s_{j}$
of (5) are real and negative (Case 1) whereas on the left side they are
complex and satisfying (3) (subcase 2.1). Therefore, we get the stable
solutions of equation (1) for all $a_{1}<0$ and $c>0$.

In the fourth quadrant on the left side of the curve $\Gamma_{1}$ both the
roots are positive and hence the system is unstable. There exists one more
curve $\Gamma_{2}:a_{1}c=\frac{-\tan^{2}(\alpha\pi/2)-1}{4}$ such that on the
right side of it, the $s_{j}$ are complex and their argument greater than
$\alpha\pi/2$ (refer (6)). Between $\Gamma_{1}$ and $\Gamma_{2}$, the roots
$s_{j}$ are complex but violating the condition (6). Thus, this bounded region
is unstable.
Note that, for any fixed $a_{1}>0$,
$\Gamma_{2}(a_{1},c)<\Gamma_{1}(a_{1},c),\forall c<0.$

Furthermore, $\Gamma_{1}$ and $\Gamma_{2}$ will not intersect each other.

Figure 1: Stability region for the non-delayed equation (1), where S=Stable
Region and U=Unstable Region.

##  IV Two-term fractional delay differential equation (FDDE)

The stability analysis of an autonomous delay differential equation

$D^{\alpha}x(t)+cD^{\beta}x(t)=ax(t)+bx(t-\tau)$ (7)

with $\alpha>\beta$ involving two fractional-order derivatives is given in
bhalekar2019analysis . In this work, we provide detailed stability and
bifurcation analysis of the following FDDE:

$D^{\alpha}x(t)+cD^{2\alpha}x(t)=ax(t)+bx(t-\tau)$ (8)

and sketch the stable/unstable regions in the parameter planes. This system
(8) reduces to the system (1) when $\tau=0$ with $a_{1}=a+b.$
The characteristics equation bhalekar2019analysis  of the given FDDE (8) is

$\lambda^{\alpha}+c\lambda^{2\alpha}=a+b\exp(-\lambda\tau).$ (9)

There is a change in stability when the characteristics root $\lambda$ crosses
the imaginary axis $\lambda=iv$, $(v>0)$ bhalekar2019analysis .
So, by substituting $\lambda=iv$ in the equation (9), as in
bhalekar2019analysis , we get

$(iv)^{\alpha}+c(iv)^{2\alpha}=a+b\exp(-iv\tau_{*}),$ (10)

where $\tau_{*}$ is the critical value of delay where change in stability can
occur.
Separating the real and imaginary parts, we get

$v^{\alpha}\cos\left(\frac{\alpha\pi}{2}\right)+cv^{2\alpha}\cos(\alpha\pi)-a=b\cos(v\tau_{*}),$
(11)
$v^{\alpha}\sin\left(\frac{\alpha\pi}{2}\right)+cv^{2\alpha}\sin(\alpha\pi)=-b\sin(v\tau_{*}).$
(12)

Squaring and adding (11) and (12), we get

$v^{2\alpha}+c^{2}v^{4\alpha}+2cv^{3\alpha}\cos(\alpha\pi/2)-2av^{\alpha}\cos(\alpha\pi/2)-2acv^{2\alpha}\cos(\alpha\pi)+a^{2}-b^{2}=0.$
(13)

Furthermore, equation (11) gives

$\tau_{+}(n)(v)=\frac{2n\pi+\arccos(\frac{v^{\alpha}\cos(\frac{\alpha\pi}{2})+cv^{2\alpha}\cos(\alpha\pi)-a}{b})}{v^{1/\alpha}},\quad
n=0,1,2,3,\ldots$ (14)

and

$\tau_{-}(n)(v)=\frac{2n\pi-\arccos(\frac{v^{\alpha}\cos(\frac{\alpha\pi}{2})+cv^{2\alpha}\cos(\alpha\pi)-a}{b})}{v^{1/\alpha}},\quad
n=1,2,\ldots.$ (15)

which are the critical values of delay for each root $v$ of the equation (13).
Note that, for given root $v$ of equation (13), either $\tau_{+}(n)$ or
$\tau_{-}(n)$ will be the critical value that will be decided by the equation
(12). So, we can take $\tau_{+}(n)(v)$ if
$\frac{v^{\alpha}\sin(\frac{\alpha\pi}{2})+cv^{2\alpha}\sin(\alpha\pi)}{-b}$
lies between $2k\pi$ to $(2k+1)\pi$ and $\tau_{-}(n)(v)$ are the critical
values if LHS of (12) lies between $(2k+1)\pi$ to $(2k+2)\pi$,
$k=0,1,2,\ldots$.

###  IV.1  The conditions for the non existence of $\tau_{*}$

Critical value $\tau_{*}$ exists if and only if equations (11) and (12) are
satisfied. Let us define the curves
$L(v)=v^{\alpha}\cos\left(\frac{\alpha\pi}{2}\right)+cv^{2\alpha}\cos(\alpha\pi)-a$
and
$R(v)=b\cos(v\tau_{*})$.
The curves $L(v)$ and $R(v)$ intersect each other if and only if the equation
(11) is satisfied. Note that if one of the equations (11) and (12) is not
satisfied then the stability will remain same as the non-delayed case i.e. if
the FDE system (1) is stable (respectively, unstable), then the FDDE system
(8) will also be stable (respectively, unstable) for all $\tau\geq 0$.
So, the non-existence of $\tau_{*}$ will give the existence of delay
independent stability regions.
In this section, we provide some (sufficient) conditions for the nonexistence
of the critical value $\tau_{*}$.

######  Theorem IV.1.

If $c>0$, $0<\alpha\leq 1/2$ and $a<-|b|$, then there does not exists critical
value $\tau_{*}$ and the stability of system (8) is independent of delay
$\tau.$

###### Proof.

For $c>0$ and $0<\alpha\leq 1/2$ we have
$v^{\alpha}\cos\left(\frac{\alpha\pi}{2}\right)+cv^{2\alpha}\cos(\alpha\pi)>0.$
So,
$-a<-a+(v^{\alpha}\cos\left(\frac{\alpha\pi}{2}\right)+cv^{2\alpha}\cos(\alpha\pi)).$

$\Rightarrow-a<L(v).$ (16)

Also by the assumption we have,

$|b|<-a.$ (17)

Therefore, from (16) and (17) we have,
$|b|<L(v).$
Further, the range of $R(v)$ is $(-|b|,|b|)$. Hence, there is no intersection
point between the two curves $L(v)$ and $R(v)$ for any $v$, $\alpha$, $c$, $a$
and $b$.
Therefore equation (11) is not satisfied. Hence $\tau_{*}$ does not exists. ∎

Note 1: Recall, $a_{1}=a+b$. In Theorem (IV.1), the condition $a<-|b|$ is
equivalent to $a_{1}<0$, if $b>0$ and $a_{1}<2b,$ if $b<0$.
This will be used to sketch the delay-independent stability region in Figure
(2).

######  Theorem IV.2.

For $c<0$ and $0<\alpha<1/2$, there will be no delay dependent stability
region if $|b|<a+\frac{1}{4c}\cos^{2}({\alpha\pi/2})\sec(\alpha\pi).$

###### Proof.

Suppose $0<\alpha<1/2$ and $c<0$. In this case, the curve $L(v)$ has local
maxima at
$v=(\frac{-\cos(\frac{\alpha\pi}{2})\sec(\alpha\pi)}{2c})^{\frac{1}{\alpha}}$
and the maximum value is
$max=-a-\frac{\cos^{2}(\frac{\alpha\pi}{2})\sec(\alpha\pi)}{4c}.$
By assumption, $max<0.$
Therefore, $L(v)$ is negative for all $v>0.$
Further, by assumption, $-|b|>max$, where $-|b|$ is the minimum value of
$R(v)$.
This shows that there is no intersection between the curves $L(v)$ and $R(v).$
Therefore, there does not exists any $\tau_{*}$ satisfying equation (11).
Hence there will be no delay dependent stability region. ∎

######  Note IV.1.

Now, we utilize Theorem (IV.2) and sketch the delay independent stability
region in $a_{1}c$-plane as below:
The condition $|b|<a+\frac{1}{4c}\cos^{2}({\alpha\pi/2})\sec(\alpha\pi)$ in
Theorem (IV.2) is equivalent to

$\displaystyle\begin{cases}a_{1}>\frac{-\cos^{2}(\frac{\alpha\pi}{2})\sec(\alpha\pi)}{4c},\textit{
if }b<0&\\\ \textsl{ and
}a_{1}>\frac{-\cos^{2}(\frac{\alpha\pi}{2})\sec(\alpha\pi)}{4c}+2b,\textit{ if
}b>0.\end{cases}$ (18)

We define the boundary curves of the regions in (18) as

$\Gamma_{3}:a_{1}=\frac{-\cos^{2}(\frac{\alpha\pi}{2})\sec(\alpha\pi)}{4c}$
and
$\Gamma_{4}:a_{1}=\frac{-\cos^{2}(\frac{\alpha\pi}{2})\sec(\alpha\pi)}{4c}+2b$,
respectively.
These regions are in the fourth quadrant of $a_{1}c$-plane because $c<0$ and
$0<\alpha<1/2$.

Recall that, the boundary of stable region of non-delayed equation (1) in this
region is described by the curve $\Gamma_{2}$ (cf. Figure (1)).

We need to find the positions of $\Gamma_{3}$ and $\Gamma_{4}$ relative to
$\Gamma_{2}$ to discuss the delay independent stability of equation (8).

$\bullet$ $\Gamma_{2}$ and  $\Gamma_{3}$:
Since $-(\tan^{2}(\alpha\pi/2)+1)>-\cos^{2}(\alpha\pi/2)\sec(\alpha\pi)$,
the curve $\Gamma_{3}$ cannot intersect $\Gamma_{2}$ and is on the right side
of $\Gamma_{2}$ in the $a_{1}c$-plane.
Thus, the intersection between the stable region bounded by $\Gamma_{2}$ and
the delay independent stable region bounded by $\Gamma_{3}$ is the region
defined by the first inequality in (18). This is sketched in Figure (2)(Fourth
quadrant). In this region, the system (8) is stable for all $\tau\geq 0$.

$\bullet$ $\Gamma_{2}$ and  $\Gamma_{4}$:
Using the similar arguments, the intersection between the regions bounded by
$\Gamma_{2}$ and $\Gamma_{4}$ is given by the second inequality in (18). This
delay independent stable region is sketched in the fourth quadrant of Figure
(2).

Now, we provide some results for the case $1/2<\alpha<1.$

######  Theorem IV.3.

The stability of system (8) is independent of delay $\tau$ if $c<0$,
$1/2<\alpha<1$ and $a<-|b|$.

###### Proof.

Note that
$v^{\alpha}\cos\left(\frac{\alpha\pi}{2}\right)+cv^{2\alpha}\cos(\alpha\pi)$
is always positive for $c<0$ and $1/2<\alpha<1.$

So, $-a<L(v).$
Therefore, $|b|<L(v)$ by assumption.
Since, $|b|$ is maximum value of $R(v),$ the two curves $L(v)$ and $R(v)$ will
never intersect each other under the assumptions of this Theorem.
Hence, the stability region of system (8) is independent of delay $\tau$. ∎

######  Theorem IV.4.

If $c>0$, $1/2<\alpha<1$ then, the stability region is delay independent if
$|b|<a+\frac{\cos^{2}(\alpha\pi/2)\sec(\alpha\pi)}{4c}$.

Note Proof of Theorem (IV.4) is analogous with that of Theorem (IV.2).

Using the similar arguments as in Note (1) and Note (2), we sketch the delay
independent stability regions of system (8) in Figure (2) using Theorem (IV.3)
and Theorem (IV.4) respectively. Note that, the system (8) is unstable for all
$\tau\geq 0$ if the conditions of Theorem (IV.3) or Theorem (IV.4) are
satisfied.

(a)

(b)

(c)

(d)

Figure 2:

##  V Geometrical method to find the stability

If

$P(\lambda)=\dfrac{\lambda^{\alpha}+c\lambda^{2\alpha}-a}{b}$ (19)

and

$Q(\lambda,\tau)=\exp{(-\lambda\tau)}$ (20)

then the characteristic equation (9) becomes

$P(\lambda)=Q(\lambda,\tau).$ (21)

Note that $\exists$ a characteristic root $\lambda_{0}$ if the graph of
$P(\lambda)$ intersects the graph of $Q(\lambda,\tau)$ at $\lambda_{0}$, for
some $\tau>0$. Moreover, the image set
$\\{Q(\lambda,\tau)\quad|\lambda\in\mathbb{C},\tau>0\\}$ will be a punctured
unit disc in $\mathbb{C}$.
If $\lambda\in\mathbb{R}$ (i.e. $Im(\lambda)=0$) then the graphs of
$\mathcal{P}$ and $Q$ are subsets of $\mathcal{R}^{2}$. In this case, the
intersection of these graphs at $\lambda>0$ will be sufficient condition for
the instability of system (8). This will be discussed in Section (V.1).
If $Im(\lambda)\neq 0$ then the graphs of $\mathcal{P}$ and $Q$ are in
$\mathcal{C}^{2}$ and we cannot observe these intersections. However, we can
compare the image sets of $\mathcal{P}$ and $Q$ with $Re(\lambda)>0$ and get
some information on the stability of equation (8). This will be provided in
Section (VI).

###  V.1 Existence of positive real root $\lambda$

Now, for any real $\lambda>0$ and $\tau>0$, $0<exp(-\lambda\tau)\leq 1$, i.e.
the range of $Q(\lambda,\tau)=(0,1]$.
Furthermore, if $\exists$ $\lambda>0$ such that $P(\lambda)\in(0,1]$ then we
can find some $\tau$ such that $P(\lambda)=Q(\lambda,\tau)$. Therefore, this
$\lambda$ is positive real root of characteristic equation.

Note that the existence of a positive real characteristic root is sufficient
for the instability of system (8). If $\exists$ $\lambda>0$ such that
$P(\lambda)\in(0,1]=Range(Q)$ then we can find some $\tau$ such that
$P(\lambda)=Q(\lambda,\tau)$. Therefore, this $\lambda$ is positive real root
of characteristic equation.

Figure 3:  Graph of $P(\lambda)$ and $Q(\lambda,\tau)$ with respect to
$\lambda$

  * •

Various conditions under which the characteristic equation has a positive real
root are described below.

    1. (1)

If Range($P$)$\supset(0,1]$ then system (8) is unstable $\forall\tau>0$.

    2. (2)

If $P(0)=-a/b<1$ and $P(\lambda)=1$ for some $\lambda>0$ then system (8) is
unstable $\forall$ $\tau\geq 0$ even if $(0,1]$ is not the subset of Range(P)
as shown in Figure (4).

(a) $\frac{-a}{b}<0<1$

(b) $0<\frac{-a}{b}<1$ and $P(\lambda)\shortrightarrow-\infty$

(c) $0<\frac{-a}{b}<1$ and $P(\lambda)\shortrightarrow\infty$

Figure 4:  $-a/b<1$ and $P(\lambda)=1$ has one positive root $\lambda$ for
some $\lambda$

  1. (3)

If $P(0)=-a/b>1$ and $P(\lambda)=0$ for some $\lambda>0$ then also
Range($P$)$\supset(0,1]$ (cf. Figure (5)) and system (8) is unstable $\forall$
$\tau\geq 0$.

(a) $\frac{-a}{b}>1$ and $P(\lambda)\shortrightarrow-\infty$

(b) $\frac{-a}{b}>1$ and $P(\lambda)$ has a local maxima

(c) $\frac{-a}{b}>1$ and $P(\lambda)$ has a local minima

Figure 5:  $-a/b>1$ and $P(\lambda)=0$ for some $\lambda$

  1. (4)

If $P(0)=-a/b>1$ and the curve $P(\lambda)$ has local minima between $0$ to
$1$ then there exists critical value $\tau_{*}$ such that
$Q(\lambda,\tau_{*})$ touches $\mathcal{P}(\lambda)$ and equation (8) is
unstable for $0<\tau<\tau_{*}$(Figure (6)).

Figure 6:  $-a/b>1$ and $P(\lambda)$ for some $\lambda>0$ has a local minima
between 0 to 1

  1. (5)

If $P(0)=-a/b<1$ and the curve $P(\lambda)$ for some $\lambda>0$ has local
maxima lies between 0 to 1 then also there exists critical value $\tau_{*}$
such that for all $\tau>\tau_{*}$ the curve $P(\lambda)$ intersects with
$Q(\lambda,\tau)$ (Figure(7)). Hence, equation (8) is unstable for
$\tau>\tau_{*}$.

Figure 7:  $-a/b<1$ and $P(\lambda)$ for some $\lambda>0$ has a local maxima
between 0 to 1

######  Theorem V.1.

The equation (8) is unstable $\forall\tau\geq 0$ if any one of the following
conditions hold:

  * (i)

$c>0$, $b<0$ and $\frac{-a}{b}>1$.

  * (ii)

$c>0$, $b>0$ and $\frac{-a}{b}<1$.

  * (iii)

$c<0$, $b<0$ and $\frac{-a}{b}<1$.

  * (iv)

$c<0$, $-b<a<\frac{-1}{4c}$ and $b<0$.

###### Proof.

(i) Given that $P(0)>1$. For $b<0$ and $\frac{-a}{b}>1$ we have $a>0$.
Since $c>0$ and $a>0$, $1+4ac>1$. Hence, one root of $P(\lambda)=0$ i.e.
$\big{(}\frac{-1+\sqrt{1+4ac}}{2c}\big{)}^{\frac{1}{\alpha}}$ is positive.
Therefore, Range($P)\supset(0,1]$. So, $P(\lambda)$ intersects with
$Q(\lambda,\tau)$ for every $\tau\geq 0$ as shown in Figure (5). Therefore,
system (8) is unstable $\forall$ $\tau\geq 0$.

(ii) We have $-a/b<1$ and $b>0$ so $a+b>0$. Hence for $c>0$, $1+4(a+b)c>1$.
Therefore, $P(\lambda)=1$ has one positive root viz.
$\Big{(}\frac{-1+\sqrt{1+4c(a+b))}}{2c}\Big{)}^{\frac{1}{\alpha}}$. So, the
curve $P(\lambda)$ intersects $Q(\lambda,\tau)$ $\forall\tau\geq 0$ by the
Figure (4) and system (8) is unstable.

(iii) If $P(\lambda)=1$ has any positive root and $P(0)=\frac{-a}{b}<1$ then,
from Figure(4) the curve $P(\lambda)$ will intersect with $Q(\lambda,\tau)$
for every $\tau\geq 0$ and $\lambda>0$. Also, by solving $P(\lambda)=1$ we get

$\lambda=\Big{(}\frac{-1\pm\sqrt{1+4c(a+b)}}{2c}\Big{)}^{\frac{1}{\alpha}}.$

Since $-a/b<1$, we get $a+b<0$. Furthermore, by assumption $c<0$. So we have
$1+4c(a+b)>0$. Hence, $\frac{-1-\sqrt{1+4c(a+b)}}{2c}>0$.
This implies that $P(\lambda)=1$ has a positive root. Hence, system (8) is
unstable $\forall\tau>0$.

(iv) We have $-b<a<\frac{-1}{4c}$
$\Rightarrow 1+4ac>0$.
Hence, for $c<0$ one root of $P(\lambda)=0$ i.e. $\frac{-1-\sqrt{1+4ac}}{2c}$
is always positive. Hence, by the Figure(5) the curve $P(\lambda)$ intersects
with $Q(\lambda,\tau)$ for every $\tau\geq 0$. So, the system (8) is unstable
$\forall\tau\geq 0$.

∎

######  Theorem V.2.

When $c<0$, $b>0$ and $-\infty<a<\frac{-1}{4c}-b$ then system (8) is unstable
$\forall\tau\geq 0$

###### Proof.

Step-1 We have $P(0)=\frac{-a}{b}>1$. If we could prove that $P(\lambda)=0$
has a positive real root $\lambda$ then from Figure(5) the curve $P(\lambda)$
will intersects with $Q(\lambda,\tau)$ for every $\tau\geq 0$. If $c<0$ then
$\frac{-1-\sqrt{1+4ac}}{2c}$ is positive because $a<0$ and $c<0$ so $1+4ac>1$.
Therefore, the root
$\lambda=\Big{(}\frac{-1-\sqrt{1+4ac}}{2c}\big{)}^{\frac{1}{\alpha}}$ of $P$
is positive. So, Range($P\supset(0,1]$.
Hence, in this case $P(\lambda)$ intersects with $Q(\lambda,\tau)$ for every
$\tau\geq 0$.
Step-2 Consider $a_{1}=0$ or $a=-b$ in this case $P(0)=1$. Also $a=-b$ implies
that $a$ is a negative number. So, $\frac{-1-\sqrt{1+4ac}}{2c}$ is positive
for $c<0$. Hence, $P(\lambda)=0$ has one positive root. So, the system (8) is
unstable for all $\tau\geq 0$ as the graph of $P(\lambda)$ intersects with the
curves $Q(\lambda,\tau)$ for every $\tau\geq 0$.
Step-3 We have $-b<a$ so $P(0)=\frac{-a}{b}<1$. Also, $a<\frac{-1}{4c}-b$
therefore, $1+4ca_{1}>0$. Hence, $\frac{-1-\sqrt{1+4c(a+b)}}{2c}>0$ for $c<0$.
So, one root of $P(\lambda)=1$ i.e.
$\Big{(}\frac{-1-\sqrt{1+4c(a+b)}}{2c}\Big{)}^{\frac{1}{\alpha}}$ is positive.
Therefore, by the Figure (4) the curve $P(\lambda)$ intersects with the curves
$Q(\lambda,\tau)$ for every $\tau\geq 0$. Hence, system (8) is unstable for
all $\tau\geq 0$. ∎

######  Theorem V.3.

If $0<a+b<\frac{-1}{4c}<a$ and $b<0$ then $\exists$
$\tau_{*}=\frac{-\log(\frac{-1-4ac}{4bc})}{\Big{(}\frac{-1}{2c}\Big{)}^{\frac{1}{\alpha}}}$
such that the system is unstable if $0\leq\tau<\tau_{*}.$

###### Proof.

Since $b<0$, the function $P$ has local minima at
$\lambda=\Big{(}\frac{-1}{2c}\Big{)}^{\frac{1}{\alpha}}$ and the minimum value
of $P$ is $\frac{-1-4ac}{4bc}$.
Also, by the assumption, $P(0)=-a/b>1$, we can expect such $\tau_{*}$ if the
minimum value will lies from 0 to 1.
We have $\frac{-1}{4c}<a$ and $c<0$ so $-1-4ac>0$. Therefore, for $b<0$,
$\frac{-1-4ac}{4bc}>0$ i.e. the minimum value is positive.
Also, $0<a+b<\frac{-1}{4c}$ we have $4c(a+b)>-1$. Therefore,
$\frac{-1-4ac}{4bc}<1$ which shows that the minimum of $P$ lies between $0$ to
$1$. Hence, there exists $\tau_{*}$ such that

$\frac{-1-4ac}{4bc}=e^{-\Big{(}\frac{-1}{2c}\Big{)}^{1/\alpha}\tau_{*}}.$

Taking log both side we get,
$\tau_{*}=\frac{-\log(\frac{-1-4ac}{4bc})}{\Big{(}\frac{-1}{2c}\Big{)}^{\frac{1}{\alpha}}}$.
So, for every $0\leq\tau<\tau_{*}$ $P(\lambda)$ intersects $Q(\lambda,\tau)$
for some $\lambda>0$. Hence, the system will be unstable $\forall
0\leq\tau<\tau_{*}$ by Figure(6). ∎

######  Theorem V.4.

If $c<0$, $b>0$, $a>-b$ and $\frac{-1}{4c}-b<a<\frac{-1}{4c}$ then

$\tau_{*}=\frac{-\log(\frac{-1-4ac}{4bc})}{\Big{(}\frac{-1}{2c}\Big{)}^{\frac{1}{\alpha}}}$
(22)

such that $\tau>\tau_{*}$ the system (8) is unstable.

###### Proof.

If $b>0$ the curve $P(\lambda)$ gives local maxima at
$\lambda=(\frac{-1}{2c})^{-1/\alpha}$ and the maximum value is
$\frac{-1-4ac}{4bc}$.
Since, $a<\frac{-1}{4c}$ so for $c<0$ we get $-1-4ac<0$. Therefore,
$\frac{-1-4ac}{4bc}>0$ for $b<0$ and $c>0$. Hence, the maximum value is
greater than 0.
Also, by assumption, $\frac{-1}{4c}<a+b$ so for $c<0$,
$\Rightarrow-1>4c(a+b)$
$\Rightarrow\frac{-1-4ac}{4bc}<1$.
So, the maximum value lies in between 0 to 1. Hence, here also there exists a
$\tau_{*}$ given by equation (22) such that for all $\tau>\tau_{*}$,
$P(\lambda)$ intersects $Q(\lambda,\tau)$ for some $\lambda>0$ from Figure(7).
So, equation (8) is unstable $\forall$ $\tau>\tau_{*}$. ∎

Using these Theorems we get the delay independent stability results as:

(a)

(b)

(c)

(d)

Figure 8:

##  VI Complex root $\lambda$ with positive real part

Now, we consider $\lambda=u+iv\in\mathbb{C}$, with $u>0$. Note that, the
boundary of the set
$\widehat{Q}=\\{Q(\lambda,\tau)\quad|\quad\lambda=u+iv,u>0,\tau\geq 0\\}$ is
the unit circle $x^{2}+y^{2}=1$ in $\mathbb{C}$.
If $|P(\lambda)|>1$, $\forall\lambda=u+iv$, $u>0$ then the image set
$\mathcal{P}=\\{P(\lambda)|\lambda=u+iv\in\mathbb{C},u>0\\}$ will not
intersect the set $\widehat{Q}$. Therefore, the characteristic equation will
not have any root $\lambda$ with positive real part.
Therefore, the bifurcation curve in the parameter plane is obtained when the
$\mathcal{P}$ touches the $\widehat{Q}$, (where $\partial\mathcal{P}$ and
$\partial\widehat{Q}$ are boundaries of the sets $\mathcal{P}$ and
$\widehat{Q}$ respectively.)
In this case, the sets $\mathcal{P}$ and $\widehat{Q}$ have a common tangent
in $\mathbb{C}$.

###  VI.1 Method to find the common tangent

Consider $P(\lambda)$ given by equation (19). For $\lambda=u+iv\in\mathbb{C}$,
${\lim_{u\to\infty}}Re(P(u+iv))$ goes to $\infty$ or $-\infty$ depending on
the sign of $b$ and $c$. If $b$ and $c$ are of same sign then
${\lim_{u\to\infty}}Re(P(u+iv))\to\infty$ and the region of $\mathcal{P}$ is
unbounded at the right end and bounded at the left as shown in Figure (12)
whereas, if $b$ and $c$ are of opposite sign then
${\lim_{u\to\infty}}Re(P(u+iv))\to-\infty$ and $\mathcal{P}$ is unbounded at
the left and bounded at the right(cf. Figure(11)).

The boundary of the set $\mathcal{P}$ is at $u=0$ i.e.
$\partial\mathcal{P}=\\{P(iv)|v\in\mathbb{R}\\}$.
Similarly, $\partial\widehat{Q}=\\{Q(iv\tau)|v\in\mathbb{R},\tau>0\\}$. The
unit circle $x^{2}+y^{2}=1$ is the boundary of $\widehat{Q}$. Since,
$\partial\mathcal{P}$ is symmetric about x-axis, we may take $v\geq 0$.
Moreover, note that the initial point of $\partial\mathcal{P}$ i.e.
$P(0)=\frac{-a}{b}$. depends on $a$ and $b$ only. Therefore, the region
$\mathcal{P}$ moves in horizontal direction if we change the parameter values
$a$ and $b$ as shown in Figures (9) and (10).

(a) $\partial\mathcal{P}$ touches $\widehat{Q}$ from the left side

(b) $\partial\mathcal{P}$ goes away from $\widehat{Q}$ as we increase $a$

Figure 9: Region $\mathcal{P}$ is on the left side of $\partial\mathcal{P}$

(a) If we increase $a$, $\partial\mathcal{P}$ goes away $\widehat{Q}$

(b) $\partial\mathcal{P}$ touches $\widehat{Q}$ from the right side

Figure 10: Region $\mathcal{P}$ is on the right side of $\partial\mathcal{P}$

If $P(iv)=x+iy$ then by separating the real and imaginary part we get

$\displaystyle
x=\dfrac{v^{\alpha}\cos(\frac{\alpha\pi}{2})+cv^{2\alpha}\cos(\alpha\pi)-a}{b}\quad\text{and}$
(23) $\displaystyle
y=\dfrac{v^{\alpha}\sin(\frac{\alpha\pi}{2})+cv^{2\alpha}\sin(\alpha\pi)}{b}.$
(24)

Since $\partial\mathcal{P}$ and $\partial\widehat{Q}$ (unit circle) have a
point say $(x,y)$ in common, it should satisfy the equations (23), (24) and
the equation of the unit circle

$x^{2}+y^{2}=1.$ (25)

In this case, at the point $(x,y)$, the slopes of the tangents to
$\partial\mathcal{P}$ and $\partial\widehat{Q}$ should match.
We find the expressions for the tangent curve in different cases, by
eliminating $v$ between equations (23) and (24) and equating the slopes.

###  VI.2 Bifurcation analysis for various values of $b$ and $\alpha$

In this section, we find the bifurcation curves separating the delay
independent stable/ unstable region with the delay dependent stable/ unstable
region. We consider four cases viz. $b>0$ and $0<\alpha<1/2$, $b<0$,
$0<\alpha<1/2$, $b>0$ and $1/2<\alpha<1$ and $b<0$ and $1/2<\alpha<1$. In the
section (VIII), we provide bifurcation curves separating different behaviors
in the delay dependent stable/ unstable regions.

####  VI.2.1  When $b>0$ and $\alpha$ between $0$ to $1/2$

Solving (24), we get
$v^{\alpha}=\dfrac{\csc(\alpha\pi)\Big{(}-\sin(\frac{\alpha\pi}{2})\pm\sqrt{\sin^{2}(\frac{\alpha\pi}{2})+4bcy\sin(\alpha\pi)}\Big{)}}{2c}.$
Note that, the stability analysis of system (8) in the I, II and III quadrants
of $a_{1}c$-plane for $b>0$ and $0<\alpha<1/2$ is provided in Theorems (IV.1),
(V.1)(ii) and (V.2). The only remaining part is fourth quadrant i.e. $a_{1}>0$
and $c<0$.

So, for $c$ negative

$v^{\alpha}=\dfrac{\csc(\alpha\pi)\Big{(}-\sin(\frac{\alpha\pi}{2})-\sqrt{\sin^{2}(\frac{\alpha\pi}{2})+4bcy\sin(\alpha\pi)}\Big{)}}{2c}.$
(26)

Putting this value of $v^{\alpha}$ in equation (23), we get the equation of
$\partial\mathcal{P}$ as

$\begin{split}\dfrac{1+4ac+4bcx+\csc(\frac{\alpha\pi}{2})\sqrt{\sin^{2}(\frac{\alpha\pi}{2})+4bcy\sin(\alpha\pi)}}{4c}-\\\
\dfrac{\cot(\alpha\pi)\csc(\alpha\pi)\big{(}\sin(\frac{\alpha\pi}{2})+\sqrt{\sin^{2}(\frac{\alpha\pi}{2})+4bcy\sin(\alpha\pi)}\big{)}^{2}}{4c}=0.\end{split}$
(27)

Equating the slopes of $\partial\widehat{Q}$ and $\partial\mathcal{P}$ using
(25) and (27), we get the equation of their common tangent as

$\frac{x}{y}+\frac{\csc(\frac{\alpha\pi}{2})\sqrt{\sin^{2}(\frac{\alpha\pi}{2})+4bcy\sin(\alpha\pi)}}{-\cot(\frac{\alpha\pi}{2})+\cot(\alpha\pi)\Big{(}1+\csc(\frac{\alpha\pi}{2})\sqrt{\sin^{2}(\frac{\alpha\pi}{2})+4bcy\sin(\alpha\pi)}\Big{)}}=0.$
(28)

Furthermore, we want the tangency conditions in terms of parameters $a_{1}$
and $c$. So, we eliminate $x$ and $y$ from the equations (25), (27) and (28)
to find the bifurcation curve $\Gamma_{5}$ in the $a_{1}c$-plane. On the right
side of $\Gamma_{5}$, there will not be any intersection between $\mathcal{P}$
and $\widehat{Q}$ (cf. Figure (11(a))) and the system (8) is stable,
$\forall\tau\geq 0$. The region bounded by the curves $\Gamma_{2}$ and
$\Gamma_{5}$ (see Figure (17) ), is stable for $\tau=0$. We observed that,
there exists a positive root of $v$ by the equation (13) in this region.
Corresponding to that positive root of $v$ we have critical values of $\tau$
given by equations (14) and (15).Note that at the smallest critical value of
$\tau_{*}$ $Re(\frac{d\lambda}{d\tau}|_{u}=0)>0$. If
$Re(\frac{d\lambda}{d\tau}|_{u}=0)<0$ at the smallest critical value
$\tau_{*}$ then we must have some characteristic root which is moving from
right to left half plane at $\tau_{*}$. Since $\nexists$ any characteristic
root in the right half plane it cannot possible. This shows that, the system
(8) is unstable for some $\tau\geq 0$. There are two behaviors viz SSR and SS.
The details will be given in Section (VIII).
First, we solve (27) for $y$. For $c<0$, the expression of y is given in the
data set-1 accompanying this article. All the data sets in this paper are also
available at
https://drive.google.com/drive/folders/147KhyNARmYlQhIt5caLDqGK9GHqvz2Yb?usp=sharing.
Further, we put this y in equation (25) and solve it for $x$. At the end, put
these values of $x$ and $y$ in the tangent equation (28) to obtain the
required bifurcation curve $\Gamma_{5}$ in terms of $a_{1}$, $b$, $c$ and
$\alpha$ for $b>0$ and $0<\alpha<1/2$ is given in the data set 1.
In the Note (IV.1) after Theorem (IV.2), we obtained a curve $\Gamma_{4}$ in
the fourth quadrant of $a_{1}c$-plane such that the system (8) is stable on
the right side of $\Gamma_{4}$. Now, we show that, the curve $\Gamma_{5}$
improves that stable region i.e. the exact value of bifurcation occurs at
$\Gamma_{5}$ which is on the left of $\Gamma_{4}$. We observed that the curve
$\Gamma_{5}$ is always on the left side of $\Gamma_{4}$. Therefore, we can
ignore $\Gamma_{4}$ and consider $\Gamma_{5}$ as the bifurcation curve in the
fourth quadrant, in this case.

(a) System (8) is stable for all $\tau\geq 0$

(b) System is stable for all $\tau$ upto the tangent equation

(c) Stability depends on delay parameter

Figure 11:  Different behavior of $P(\lambda)$ when $\lambda$ is a complex
root and $b$ and $c$ are of opposite sign

####  VI.2.2 We consider $b<0$ and $0<\alpha<1/2$.

  * Case 1:

$c>0$

The Theorem (IV.1) in Section (IV) and Theorem (V.1)(a) in Section (V.1)
provide the stability analysis except for $2b<a_{1}<0$ or $b<a<-b$. We discuss
this region in following subcases.

  * Subcase 1.1:

$b<a<0$

Subcase: When $b<a<0$

$Re\Big{[}\dfrac{1}{\alpha}\dfrac{d\lambda}{d\tau}\Big{|}_{u=0}\Big{]}=v^{2\alpha}+(2)c^{2}v^{4\alpha}-av^{\alpha}\cos(\frac{\alpha\pi}{2})-2acv^{2\alpha}\cos(\alpha\pi)+3cv^{3\alpha}\cos(\frac{\alpha\pi}{2})$
(29)

The sign of $\dfrac{1}{\alpha}\dfrac{d\lambda}{d\tau}$ bhalekar2019analysis
is positive for $a<0$, $c>0$ and $0<\alpha<1/2$.

  * Subcase 1.2:

$0<a<-b$

We have from equation (16) given in bhalekar2019analysis

$-2acv^{2\alpha}\cos(\alpha\pi)-av^{\alpha}\cos(\frac{\alpha\pi}{2})=-v^{2\alpha}-c^{2}v^{4\alpha}-2cv^{3\alpha}\cos(\frac{\alpha\pi}{2})+av^{\alpha}\cos(\frac{\alpha\pi}{2})+b^{2}-a^{2}$
(30)

Putting (30) in (29) we get

$Re\Big{[}\dfrac{1}{\alpha}\dfrac{d\lambda}{d\tau}\Big{|}_{u=0}\Big{]}=c^{2}v^{4\alpha}+cv^{3\alpha}\cos(\frac{\alpha\pi}{2})+av^{\alpha}\cos(\frac{\alpha\pi}{2})+b^{2}-a^{2}>0$

Therefore,
$Re\Big{[}\dfrac{1}{\alpha}\dfrac{d\lambda}{d\tau}\Big{|}_{u=0}\Big{]}>0$ in
both cases. Hence, by Theorem 3.1 given in bhalekar2019analysis  system (8)
undergoes Hopf bifurcation at

$\tau_{*}=\frac{\arccos(\frac{v^{\alpha}\cos(\frac{\alpha\pi}{2})+cv^{2\alpha}\cos(\alpha\pi)-a}{b})}{v}$
(31)

such that the system (8) is stable for all $0<\tau<\tau_{*}$ and unstable for
$\tau>\tau_{*}$ in the region $2b<a_{1}<0$ when $b<0$, $c>0$ and
$0<\alpha<1/2$ (cf. Figure (18)).

  * Case 2:

$c<0$

Theorem (V.1)(c) shows that the 3rd quadrant is unstable $\forall\tau\geq 0$.
Let us consider the fourth quadrant in the $a_{1}c$-plane. Since, $b<0$ and
$c<0$, we have ${\lim_{\lambda\to\infty}}Re(P(\lambda))\to\infty$ i.e. the
region bounded by $P(iv)$ is on the right side of $\partial\mathcal{P}$ as
shown in Figure (12). As in Section (VI.1), we have a bifurcation curve
$\Gamma_{6}$ which is obtained by using the condition for the curves
$\partial\mathcal{P}$ and $\partial\widehat{Q}$ to have a common tangent.
The expression for $\Gamma_{6}$ is given in the data set 2 accompanying this
paper.
If we take the parameter values $(a_{1},c)$ on the right of $\Gamma_{6}$ then
$\mathcal{P}$ does not intersect $\widehat{Q}$ and the system (8) is stable
for all $\tau\geq 0$. We observe that, the bifurcation curve $\Gamma_{3}$
obtained in Section (IV) is always on the right of $\Gamma_{6}$ (as shown in
Figure (17)). We already have shown that the system (8) is stable
$\forall\tau\geq 0$ on the right of $\Gamma_{3}$. Therefore, $\Gamma_{6}$
provides a better estimate for the bifurcation curve and we can ignore
$\Gamma_{3}$.

(a) System (8) is stable for all $\tau\geq 0$

(b) System (8) is stable for all $\tau\geq 0$ upto the tangent equation

Figure 12:  Different behavior of $P(\lambda)$ when $\lambda$ is a complex
root and $b$ and $c$ are of same sign

####  VI.2.3 We consider $b>0$ and $1/2<\alpha<1$.

  * Case 1:

$c>0$

For $b>0$ and $c>0$ the region of $\mathcal{P}$ will be on the right side of
$\partial\mathcal{P}$ as ${\lim_{u\to\infty}}Re(P(u+iv))\to\infty$ and if we
fix $b$, $\alpha$, $c$ and change the parameter $a$, the region of
$\mathcal{P}$ gets translated from right side to the left side in the
$a_{1}c$-plane (cf. Figure (12)). So, here also we have to find the tangent
equation as given in section (VI.1). So, note that for $c>0$, we get

$v^{\alpha}=\dfrac{\csc(\alpha\pi)\Big{(}-\sin(\frac{\alpha\pi}{2})+\sqrt{\sin^{2}(\frac{\alpha\pi}{2})+4bcy\sin(\alpha\pi)}\Big{)}}{2c}.$
(32)

By putting $v^{\alpha}$ in equation (23), we get the expression for
$\partial\mathcal{P}$

$\begin{split}\dfrac{1+4ac+4bcx-\csc(\frac{\alpha\pi}{2})\sqrt{\sin^{2}(\frac{\alpha\pi}{2})+4bcy\sin(\alpha\pi)}}{4c}-\\\
\dfrac{\cot(\alpha\pi)\csc(\alpha\pi)\big{(}\sin(\frac{\alpha\pi}{2})+\sqrt{\sin^{2}(\frac{\alpha\pi}{2})+4bcy\sin(\alpha\pi)}\big{)}^{2}}{4c}=0.\end{split}$
(33)

If we equate slope of circle (25) and slope of $\partial\mathcal{P}$ (33), we
get the expression for common tangent as

$\frac{x}{y}+\frac{\csc(\frac{\alpha\pi}{2})\sqrt{\sin^{2}(\frac{\alpha\pi}{2})+4bcy\sin(\alpha\pi)}}{\cot(\frac{\alpha\pi}{2})+\cot(\alpha\pi)\Big{(}-1+\csc(\frac{\alpha\pi}{2})\sqrt{\sin^{2}(\frac{\alpha\pi}{2})+4bcy\sin(\alpha\pi)}\Big{)}}=0.$
(34)

Using equations (33), (34) and (25), we can eliminate $x$ and $y$ as given in
the section (VI.1). We get the bifurcation curve $\Gamma_{7}$ given in data
set 3 in the $a_{1}c$-plane. On the left side of $\Gamma_{7}$, we don’t have
any intersection between $\mathcal{P}$ and $\widehat{Q}$ (as shown in Figure
(12(a))) and system (8) is stable $\forall\tau\geq 0$. Note that, the
bifurcation curve $\Gamma_{7}$ intersects the $c-$axis in the $a_{1}c$ plane
at some value $c=c_{2}$. The expression for $c_{2}$ in terms of $b>0$ and
$1/2<\alpha<1$ is given in the data set 4 which is given in this article.
The curve $\Gamma_{7}$ provides a bifurcation in second quadrant of
$a_{1}c$-plane in this case.
Note that, Theorem (V.1)(b) shows that the first quadrant of $a_{1}c$-plane in
this case is unstable region $\forall\tau\geq 0$.

  * Case 2:

$c<0$

Note that for $c<0$, the region of $\mathcal{P}$ is on the left side of the
boundary $\partial\mathcal{P}$ as in Figure (11(a)). For $c<0$, the value of
$v^{\alpha}$ is given by the equation (26) and the common tangent equation is
given by (28). As in Section (VI.1), we eliminate $x$ and $y$ from (26), (28)
and (25) to get the bifurcation curve. Let us call the equation of the
bifurcation curve for $b>0$, $c<0$ and $1/2<\alpha<1$ as $\Gamma_{8}$ (see
data set 5). On the right side of $\Gamma_{8}$, we don’t have any intersection
between $\mathcal{P}$ and $\widehat{Q}$ and the system (8) is stable for all
$\tau\geq 0$.
By Theorem (V.1)(c) in Section (V.1) the $3rd$ quadrant of $a_{1}c$-plane is
unstable $\forall$ $\tau\geq 0$, for this case.

####  VI.2.4 Consider $b<0$ and $1/2<\alpha<1$.

  * Case 1:

$c>0$.

If $b<0$, $c>0$ and $1/2<\alpha<1$, then $\mathcal{P}$ is bounded from right
and unbounded from left as ${\lim_{u\to\infty}}Re(P(u+iv))\to-\infty$ (cf.
Figure (11)). If we fixed $b$, $\alpha$, $c$ and increase $a$ the region of
$\mathcal{P}$ goes away from $\widehat{Q}$ (cf. Figure (11(a))). If we
decrease $a$ then we observed that the boundary $\partial\mathcal{P}$ and
boundary $\widehat{Q}$ intersect as shown in Figure (13(a)) for a fix $b$, $c$
and $\alpha$. It means both $\partial\mathcal{P}$ and $\partial\widehat{Q}$
have a common tangent. For this, we have to find the expression for
$\partial\mathcal{P}$. So, for $c>0$ we have $v^{\alpha}$ is equation (32)
putting this $v^{\alpha}$ in equation (23) we get $\partial\mathcal{P}$ as
equation (33). Therefore, the expression for common tangent is obtained by
equation (34). So, as in Section (VI.1) we find a bifurcation curve
$\Gamma_{9}$ (see data set 6) such that on the left side of $\Gamma_{9}$ we
don’t have any intersection between $\mathcal{P}$ and $\widehat{Q}$ and system
(8) is stable $\forall$ $\tau\geq 0$. Note that this bifurcation curve
$\Gamma_{9}$ intersect with line $a_{1}=2b$ at $c_{7}$. The expression for
$c_{7}$ is given in data set 7 attached with this manuscript.

(a) Bifurcation curve when $a_{1}<2b<0$

(b) Bifurcation curve when $2b<a_{1}<0$

  * Case 2:

$c<0$.
For $c<0$ and $b<0$ the region of $\mathcal{P}$ is unbounded from right and
bounded from left (given as Figure (12)). So, here also we have to find the
condition such that the boundary $\partial\mathcal{P}$ touches to boundary of
$\widehat{Q}$ i.e. unit circle. For that we need to equate the tangent of
$\partial\mathcal{P}$ and $\partial\widehat{Q}$. Therefore, for $c<0$ the
value of $v^{\alpha}$ is given by the equation (26). The equation of
$\partial\mathcal{P}$ and equation of common tangent are (27) and (28)
respectively. So, we obtain the bifurcation curve $\Gamma_{10}$ (see data set
8) as given in Section (VI.1) such that on the right side of $\Gamma_{10}$ we
don’t have any common points between $\mathcal{P}$ and $\widehat{Q}$ and
system (8) is stable for all $\tau\geq 0$.

Note that the first and third quadrant in the $a_{1}c$-plane is unstable for
all $\tau\geq 0$ by Theorems (V.1)(a) and (V.1)(c) given in Section (V.1).

##  VII Few more unstable regions

Based on our observations for sufficient number of parameter values, we
conjecture that the following regions are unstable:

  * •

The regions bounded by the curve $\Gamma_{1}$ and $\Gamma_{2}$ in the Figure
(17) and (18).

It is open problem to prove these conjectures or to provide counter examples.
If any of these regions is not unstable then we can expect the instability
switches.

##  VIII Delay dependent stability/ instability:

###  VIII.1 The regions SS and SSR

We observed following types of bifurcations which depend on delay.
Single stable region (SSR): In this case, $\exists$ $\tau_{*}>0$ such that
$0<\tau<\tau_{*}$ $\Longrightarrow$ the system is stable and $\tau>\tau_{*}$
$\Longrightarrow$ the unstable behavior (cf. Figure (15)).
Stability Switch (SS): The stability switches are observed if $\exists$
$\tau_{0*}=0$ and positive constants $\tau_{1*}$, $\tau_{2*}$,
$\tau_{3*},\ldots\tau_{k*}$ such that
$\tau_{2jk*}<\tau<\tau_{(2j+1)*}$ $\Longrightarrow$ stable,
$j=0,1,\ldots,\frac{k-1}{2}$,
$\tau_{(2j+1)*}<\tau<\tau_{(2j+2)*}$ $\Longrightarrow$ unstable,
$j=0,1,\ldots,\frac{k-3}{2}$ and
$\tau>\tau_{k}$ $\Longrightarrow$ unstable.
Note that, $k$ is an odd number $\geq 3$ e.g. if $k=5$, then we get the
stability properties as shown in Figure (14) for the delay value $\tau$.

Figure 14: Stability switches for $k=5$

If we allow $k=1$ then this reduces to the case SSR and the properties will be
as in the Figure (15):

Figure 15: Single stable region for $k=1$

We provide more details and the conditions for such switches in further
sections.
In the previous Sections, we provided some delay independent stable/unstable
regions. These are some regions in Figure((17)-(20)) which were not discussed.
Now, we show that the stability properties in these regions depend on the
delay. The regions

  * (i)

bounded by $\Gamma_{2}$ and $\Gamma_{5}$ for $b>0$, $c<0$, $0<\alpha<1/2$ in
Figure (17)

  * (ii)

bounded by $a_{1}=2b$ and $2b<a_{1}<0$, for $c>0$, $b<0$, $0<\alpha<1/2$ in
Figure (18)

  * (iii)

bounded by $\Gamma_{6}$ and $\Gamma_{2}$ for $c<0$, $b<0$, $0<\alpha<1/2$ in
Figure (19)

  * (iv)

bounded by $\Gamma_{7}$ and the line $a_{1}=0$ for $b>0$, $c>0$,
$1/2<\alpha<1$ in Figure (19).

  * (v)

bounded by $\Gamma_{2}$ and $\Gamma_{8}$ for $b>0$, $c<0$ and $1/2<\alpha<1$
in Figure (19)

  * (vi)

bounded by $\Gamma_{9}$ and the line $a_{1}=0$ for $c>0$, $b<0$ and
$1/2<\alpha<1$ in Figure (20)

  * (vii)

bounded by $\Gamma_{2}$ and $\Gamma_{10}$ for $c<0$, $b<0$ and $1/2<\alpha<1$
in Figure (20) are subsets of the stable region at $\tau=0$ given in Figure
(1).
The system (8) is stable region at $\tau=0$ in all the above regions.

We observed that $\exists$ a real root of $v^{\alpha}$ of equation (13) in all
these regions. Since, the system is stable at $\tau=0$, $\nexists$ any
characteristic root in the right half complex plane. As $exists$ v, we can
have critical values of $\tau$ defined by (14) and (15). If
$Re\Big{(}\frac{d\lambda}{d\tau}\Big{|}_{u=0}\Big{)}<0$ at the smallest
critical values $\tau_{*}$ then we must have some characteristic root which is
moving from right half plane to left half plane at $\tau_{*}$. Since,
$\nexists$ any characteristic root in the right half this cannot be true.
Therefore, we must have
$Re\Big{(}\frac{d\lambda}{d\tau}\Big{|}_{u=0}\Big{)}>0$ at $\tau_{*}$.
This shows that $exists$ some $\tau_{0}>\tau_{*}$ such that system (8) becomes
unstable at $\tau_{0}$.

Instability switches (IS): In this case, we will have $\tau_{0*}=0$ and
positive constants $\tau_{1*}$, $\tau_{2*},\ldots,\tau_{k*}$ such that
$\tau_{2j*}<\tau<\tau_{(2j+1)*}\Longrightarrow$ Unstable,
$j=0,1,\ldots,(k-2)/2$,
$\tau_{(2j+1)*}<\tau<\tau_{(2j+2)*}\Longrightarrow$ stable,
$j=0,1,\ldots,(k-2)/2$
and $\tau>\tau_{k}\Longrightarrow$ Unstable.
E.g. if $k=6$ then the stability properties are as in Figure (16).

Figure 16: Instability switches for $k=6$

In both the cases SS and IS, the last unbounded region for $\tau$ gives
instability. As we go arbitrary close to some bifurcation values, the number
of switches increase.
We consider four cases viz. $b>0$ and $0<\alpha<1/2$, $b<0$ and
$0<\alpha<1/2$, $b>0$ and $1/2<\alpha<1$ and $b<0$ and $1/2<\alpha<1$. In each
cases, we give the expressions for the bifurcation curves that separate SSR,
SS, IS and the unstable region.

###  VIII.2 We consider $b>0$ and $0<\alpha<1/2$.

In the Sections (VI) and (V), we provided the stability analysis of all the
points in the $a_{1}c$-plane except those in the region bounded by the curves
$\Gamma_{2}$ and $\Gamma_{5}$ in the fourth quadrant. Solving equation (13)
under these conditions on parameters $b$ and $\alpha$, we get two roots viz.
$v_{1}$ and $v_{2}$ with $0<v_{1}<v_{2}$(see data set 9). These values will
produce the following critical values of delay viz. $\tau_{-}(n)(v_{1})$ and
$\tau_{+}(n)(v_{2})$ (see (14) and (15) for the expressions).

We observed that, $Re\big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\big{)}$ is negative
at $v_{1}$ and positive at $v_{2}$ for these values of parameters. Therefore,
the characteristic roots move from left to right at $\tau=\tau_{+}(n)$ and
from right to left at $\tau=\tau_{-}(n)$. Since, the system (8) is stable at
$\tau=0$ in the fourth quadrant of $a_{1}c$-plane and $\tau_{+}(0)$ is
smallest critical value, the system becomes unstable at $\tau=\tau_{+}(0)$.
Therefore, $0<\tau<\tau_{+}(0)\Longrightarrow$ stability.
This can be the only stable region (SSR case) or we can have stability
switches (SS). We provide the bifurcation analysis below.
Let $\delta_{1}=\tau_{-}(n+1)-\tau_{-}(n)=\frac{2\pi}{v_{1}}$,
$\delta_{2}=\tau_{+}(n+1)-\tau_{+}(n)=\frac{2\pi}{v_{2}}$
and $\upmu=\tau_{-}(1)-\tau_{+}(0)$.
On the right side of the curve $\Gamma_{2}$ system (8) is stable at $\tau=0$
by the Section (III). Note that $\tau_{+}(0)$ is the closest critical value to
$0$ where $Re\left(\dfrac{d\lambda}{d\tau}|_{\tau_{+}(0)}\right)>0$.
If $\upmu<\delta_{2}$ then,
$0<\tau_{+}(0)<\tau_{-}(1)<\tau_{+}(1)<\ldots<\tau_{+}(k)$ for some $k>0$ and
we have stability switches (SS) upto $\tau_{+}(k)$. So, the switches will get
disappear if $\exists$ smallest $k$ such that
$\upmu_{k}=\tau_{-}(k+1)-\tau_{+}(k)>\delta_{2}$. Moreover, system become
unstable for all $\tau>\tau_{+}(k)$.
If $\upmu>\delta_{2}$ then, $0<\tau_{+}(0)<\tau_{+}(1)<\tau_{-}(1)\ldots$ and
we have only single stable region (SSR) $[0,\tau_{+}(0))$. Therefore, the
equation $\upmu=\delta_{2}$ is the bifurcation curve $\Gamma_{11}$ in the
$a_{1}c$-plane. On the right side of $\Gamma_{11}$ we have SS and on the left
we have SSR.
For $b>0$ and $0<\alpha<1/2$ the expression of curve $\Gamma_{11}$ is attached
with this manuscript in the data set 10.
Using all this analysis, we are able to provide the stability region of system
(8) in the $a_{1}c$-plane for $b>0$ and $0<\alpha<1/2$ (cf. Figure (17)).

###  VIII.3 We consider $b<0$ and $0<\alpha<1/2$.

For $b<0$ and $0<\alpha<1/2$, we provided the stability results in Section
(VI.2.2) and Theorem (V.1)(d) (Section (V.1)) except the region bounded by the
curve $\Gamma_{0}$, vertical axis with $c<0$ and $\Gamma_{7}$ in $a_{1}c$
plane. We now propose the stability properties of system (8) in this remaining
part. If we solve equation (13) with $b<0$ and $0<\alpha<1/2$, we get two
values of $v^{\alpha}$ viz. $v_{1}$ and $v_{2}$ with $0<v_{1}<v_{2}$ (see data
set 9). The critical values of $\tau$ corresponding to $v_{1}$ and $v_{2}$ are
$\tau_{-}(n)(v_{1})$ and $\tau_{-}(n)(v_{2})$.
We proceed as in Section (VIII.2). If
$\delta_{1}=\tau_{-}(n+1)(v_{1})-\tau_{-}(n)(v_{1})=\frac{2\pi}{v_{1}}$,
$\delta_{2}=\tau_{-}(n+1)(v_{2})-\tau_{-}(n)(v_{2})=\frac{2\pi}{v_{2}}$
and $\upmu=\tau_{-}(1)(v_{2})-\tau_{-}(1)(v_{1})$
then $\upmu=\delta_{2}$ is the bifurcation curve in the $a_{1}c$-plane if
system (8) is stable at $\tau=0$ i.e. the right side of the curve
$\Gamma_{2}$. Let us call the bifurcation curve $\upmu=\delta_{2}$ as
$\Gamma_{12}$ whose expression in the $a_{1}c$-plane for $b<0$ and
$0<\alpha<1/2$ is given in data-set 11. Moreover, the curve $\upmu=\delta_{2}$
intersects the curve $\Gamma_{2}$ at $c_{0}$ (Figure (18)). The expression for
$c_{0}$ is given in data set 12. On the right side of the curve $\Gamma_{12}$,
$\upmu<\delta_{2}$. So, we have
$0<\tau_{-}(1)(v_{2})<\tau_{-}(1)(v_{1})<\ldots\tau_{-}(k)(v_{2})$ and we have
SS upto $\tau_{-}(k)(v_{2})$. Note that we expect SS only upto
$\tau_{-}(k)(v_{2})$ where $k$ is the smallest positive number such that
$\upmu_{k}=\tau_{-}(k)(v_{1})-\tau_{-}(k)(v_{2})>\delta_{2}$. System is
unstable for all $\tau>\tau_{-}(k)(v_{2})$. On the left side of the curve
$\upmu=\delta_{2}$ we have $0<\tau_{-}(1)(v_{2})<\tau_{-}(2)(v_{2})<\ldots$
and it will give SSR from $[0,\tau_{-}(1)(v_{2}))$.

####  VIII.3.1 Condition for the Instability switch (IS):

Note that region bounded by $\Gamma_{0}$ and $\Gamma_{2}$ is still remaining
for $b<0$ and $0<\alpha<1/2$. By the Section (III) and Figure (1), the region
on the left side of the curve $\Gamma_{2}$ is unstable at $\tau=0$. If we
solve equation (13) with $b<0$ and $0<\alpha<1/2$ we get $v_{1}$ and $v_{2}$
as two roots of $v^{\alpha}$ with $0<v_{1}<v_{2}$. In this region, the
critical values of delay are $\tau_{+}(n)(v_{1})$ and $\tau_{-}(n)(v_{2})$.

We notice that $Re\big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\big{)}$ is negative at
$v_{1}$ and positive at $v_{2}$. It means that the roots of characteristic
equation shift from right half plane to left half plane at $\tau=\tau_{+}$ and
from left to right at $\tau=\tau_{-}$. Since, the system (8) is unstable at
$\tau=0$ and if $\tau_{-}(0)$ is the smallest critical value then multiplicity
of positive root in the right half plane is increased by two. On the other
hand if $\tau_{+}(0)$ is the smallest critical value then the multiplicity of
positive root in the right half plane is decreased by two. The former will
provide IS and the later will give unstable solutions $\forall\tau>0$.
Therefore, the curve bifurcating the unstable region with the IS region is
provided by $\tau_{-}(0)=\tau_{+}(0)$ i.e.

$\frac{\arccos(\frac{v_{1}\cos(\frac{\alpha\pi}{2})+cv_{1}\cos(\alpha\pi)-a}{b})}{v_{1}^{1/\alpha}}-\frac{2\pi-\arccos(\frac{v_{2}\cos(\frac{\alpha\pi}{2})+cv_{2}\cos(\alpha\pi)-a}{b})}{v_{2}^{1/\alpha}}=0.$

Let us call this curve as $\Gamma_{13}$ (see data set 13). On the right side
of $\Gamma_{13}$ we have $\tau_{+}(0)<\tau_{-}(1)$. So, the smallest critical
delay is $\tau_{+}(0)$. By the Section (III), we have only two roots of the
characteristic equation at $\tau=0$ in the right half plane and those two
roots shift from right half plane to left half plane as we increase $\tau$
from $\tau_{+}(0)$. Since, we don’t have any root on the right half plane, the
system becomes stable. As we further increase $\tau>\tau_{-}(1)$ the
characteristic roots are shifted from left half plane to the right half plane
and the system will become unstable. So, we have instability switches if
$\tau_{+}(0)$ is the smallest critical value.
Let $\delta_{1}=\tau_{+}(n+1)-\tau_{+}(n)=\frac{2\pi}{v_{1}}$,
$\delta_{2}=\tau_{-}(n+1)-\tau_{-}(n)=\frac{2\pi}{v_{2}}$
and $\upmu=\tau_{-}(1)-\tau_{+}(0)$.

  * (a)

If $\tau_{+}(0)$ is the smallest critical value and $\upmu<\delta_{1}$ then we
have instability switches.

  * (b)

If, in addition, $\upmu>\delta_{1}-\delta_{2}$ i.e.
$\upmu<\delta_{1}<\upmu+\delta_{2}$ then we have switch of the form $U-S-U-
S-U$. Otherwise (i.e. if $\upmu<\delta_{1}-\delta_{2}$) we have switch of the
form $U-S-U$.

So, we have instability switch upto smallest $k$ such that
$0<\tau_{+}(0)<\tau_{-}(1)<\tau_{+}(1)\ldots\tau_{-}(k)$ where
$\upmu_{k}=\tau_{-}(k+1)-\tau_{+}(k)>\delta_{1}$.
On the left side of $\Gamma_{13}$, $\tau_{-}(1)>\tau_{+}(0)$ so the
multiplicity of positive root on the right side increase by two and system
becomes unstable $\forall\tau\geq 0$.
Using this analysis, we are able to sketch the stability region in the
complete $a_{1}c$-plane for $b<0$ and $0<\alpha<1/2$ (cf. Figure (18)).

###  VIII.4  Consider $b>0$ and $1/2<\alpha<1$.

  * Case 1:

$c>0$.
If we solve equation (13), we get two positive values of $v^{\alpha}$ namely
$v_{1}$ and $v_{2}$ with $0<v_{1}<v_{2}$ as in previous sections.
Corresponding to those $v_{1}$ and $v_{2}$ the critical values of delay $\tau$
are $\tau_{-}(n)(v_{1})$ and $\tau_{-}(n)(v_{2})$.
We observed that, $Re\big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\big{)}$ is negative
at $v_{1}$ and positive at $v_{2}$ for these values of parameters. Now, we
define
$\delta_{1}=\tau_{-}(n+1)(v_{1})-\tau_{-}(n)(v_{1})=\frac{2\pi}{v_{1}}$,
$\delta_{2}=\tau_{-}(n+1)(v_{2})-\tau_{-}(n)(v_{2})=\frac{2\pi}{v_{2}}$
and $\upmu=\tau_{-}(1)(v_{2})-\tau_{-}(1)(v_{1})$.
Note that $\upmu=\delta_{2}$ is the bifurcation curve which we call as
$\Gamma_{14}$ (see data set 14). On the left side of $\Gamma_{14}$, we have
$\upmu<\delta_{2}$ and
$0<\tau_{-}(1)(v_{2})<\tau_{-}(1)(v_{1})<\ldots\tau_{-}(k)(v_{2})$ which gives
SS upto $\tau_{-}(k)(v_{2})$. Moreover, $k$ is the smallest positive integer
where $\upmu_{k}=\tau_{-}{(k)}(v_{2})-\tau_{-}{(k)}(v_{1})>\delta_{2}$ and if
$\tau>\tau_{-}(k)$ then system becomes unstable. On the right of
$\Gamma_{14}$, we have $\upmu>\delta_{2}$ so we get only SSR from
$[0,\tau_{-}(1)(v_{2}))$ as the critical values occur as
$0<\tau_{-}(1)(v_{2})<\tau_{-}(2)(v_{2})<\ldots$. Note that this bifurcation
curve $\Gamma_{14}$ intersects the vertical axis at $c_{1}$. The expression
for $c_{1}$ is given in data set 15 of this manuscript.

  * Case 2:

$c<0$.
The critical values of delay $\tau$ with respect to $v_{1}$ and $v_{2}$ are
$\tau_{-}(n)$ and $\tau_{+}(n)$ respectively. Note that
$Re\big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\big{)}$ is negative at $v_{1}$ and
positive at $v_{2}$ for these values of parameters. We define,
$\delta_{1}=\tau_{-}(n+1)-\tau_{-}(n)=\frac{2\pi}{v_{1}}$,
$\delta_{2}=\tau_{+}(n+1)-\tau_{+}(n)=\frac{2\pi}{v_{2}}$
and $\upmu=\tau_{+}(0)-\tau_{-}(1)$.
The bifurcation curve $\Gamma_{15}$ is obtained by solving $\upmu=\delta_{2}$
(see dataset-16). Note that the curve $\Gamma_{15}$ meets the curve
$\Gamma_{2}$ at $c_{3}$. The expression for $c_{3}$ is also given in data set
17. On the right side of $\Gamma_{15}$ we have SS and on the left side we have
SSR.

###  VIII.5  Consider $b<0$ and $1/2<\alpha<1$.

  * Case 1:

$c>0$.
If we solve equation (13) we get two values of $v^{\alpha}$ as $v_{1}$ and
$v_{2}$ with $0<v_{1}<v_{2}$. The critical values of $\tau$ with respect to
$v_{1}$ are $\tau_{+}(n)(v_{1})$ and $\tau_{+}(n)(v_{2})$. Here we observed
that $Re\Big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\Big{)}$ is negative at $v_{1}$
and positive at $v_{2}$. We set
$\delta_{1}=\tau_{+}(n+1)(v_{1})-\tau_{+}(n)(v_{1})=\frac{2\pi}{v_{1}}$,
$\delta_{2}=\tau_{+}(n+1)(v_{2})-\tau_{+}(n)(v_{2})=\frac{2\pi}{v_{2}}$
and $\upmu=\tau_{+}(0)(v_{2})-\tau_{+}(0)(v_{2})$.
When $\upmu<\delta_{2}$, we have
$0<\tau_{+}(0)(v_{2})<\tau_{+}(0)(v_{1})<\tau_{+}(1)(v_{2})\ldots\tau_{+}(k)(v_{2})$
and we have SS upto $\tau_{+}(k)(v_{2})$ where $k$ is the smallest positive
integer such that
$\upmu_{k}=\tau_{+}(k)(v_{2})-\tau_{+}(k)(v_{1})>\delta_{2}$. Moreover, when
$\upmu>\delta_{2}$, we get
$0<\tau_{+}(0)(v_{2})<\tau_{+}(1)(v_{2})<\tau_{+}(0)(v_{1})\ldots$ and we get
SSR $\tau\in[0,\tau_{+}(0)(v_{2}))$. So, $\upmu=\delta_{2}$ is the bifurcation
curve $\Gamma_{16}$ for $b<0$ and $1/2<\alpha<1$ (see Figure (20)). See the
dataset-18 for the expression of $\Gamma_{16}$. Note that the curve
$\Gamma_{16}$ and the tangent curve $\Gamma_{9}$ will meet at $(a_{5},c_{5})$.
For $b=-1$ and $\alpha=0.8$ we get $(a_{5},c_{5})=(-1.97444,2.00571)$.

  * Case 2:

$c<0$.
Let us consider the region bounded by the curve $\Gamma_{2}$ and the tangent
curve $\Gamma_{10}$. In this region, we have two values $v_{1}$ and $v_{2}$
with $0<v_{1}<v_{2}$ if we solve equation (13). Corresponding to $v_{1}$, we
have critical values of delay $\tau$ as $\tau_{-}(n)(v_{1})$ and with respect
to $v_{2}$ the critical value of $\tau$ as $\tau_{-}(n)(v_{2})$. We define,
$\delta_{1}=\tau_{-}(n+1)(v_{1})-\tau_{-}(n)(v_{1})=\frac{2\pi}{v_{1}}$,
$\delta_{2}=\tau_{-}(n+1)(v_{2})-\tau_{-}(n)(v_{2})=\frac{2\pi}{v_{2}}$
and $\upmu=\tau_{-}(1)(v_{2})-\tau_{-}(1)(v_{1})$.
We have $\upmu=\delta_{2}$ as the bifurcation curve $\Gamma_{17}$ given in
data-set 19 (see Figure (20)). On the right side of $\Gamma_{17}$ the system
is stable at $\tau=0$ and the smallest critical value of $\tau$ is
$\tau_{-}(1)(v_{2})$. So, we have
$0<\tau_{-}(1)(v_{2})<\tau_{-}(1)(v_{1})<\tau_{-}(2)(v_{2})\ldots\tau_{-}(k)(v_{2})$
and we get stability switches upto $\tau_{-}(k)(v_{2})$. On the left side of
$\Gamma_{17}$, we have SSR $[0,\tau_{-}(1)(v_{2}))$. Note that the curve
$\Gamma_{17}$ meets with the curve $\Gamma_{2}$ at $c_{6}$. Note that some
region which is on the left side of the curve $\Gamma_{2}$ is yet to be done
(Figure (20)). On the left side of the curve $\Gamma_{2}$, system (8) is
unstable at $\tau=0$. If we solve equation (13) in the region bounded by
$\Gamma_{2}$, vertical axis with $c<0$ and the curve $\Gamma_{0}$ we get two
values of $v$ say $v_{1}$ and $v_{2}$ with $0<v_{1}<v_{2}$ given in the data
set 9. These $v_{1}$ and $v_{2}$ give the critical values of delay $\tau$ as
$\tau_{+}(n)(v_{1})$ and $\tau_{+}(n)(v_{2})$.
We notice that $Re\Big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\Big{)}$ is negative at
$v_{1}$ and positive at $v_{2}$. So, if
$\tau_{+}(0)(v_{1})<\tau_{+}(0)(v_{2})$ then two roots are going from right
half plane to the left half plane at $\tau+(0)(v_{1})$ and no roots are on the
right side. On the other hand, the roots are moving from the left to the right
half plane at $\tau_{+}(0)(v_{2})$ and the system becomes unstable. So, system
(8) is stable for $\tau_{+}(0)(v_{1})<\tau<\tau_{+}(0)(v_{2})$. If
$\tau_{+}(0)(v_{1})>\tau_{+}(0)(v_{2})$ it means that the system (8) unstable
at $\tau=0$ and some more roots of characteristic equation (9) goes from left
half plane to the right half plane. So, the system becomes unstable for
$\forall\tau\geq 0$. So, the bifurcation curve $\Gamma_{18}$ given in data set
20 that discriminate IS from unstable $\forall\tau$ is

$\tau_{+}(0)(v_{1})=\tau_{+}(0)(v_{2}).$

In the region bounded by $\Gamma_{2}$ and $\Gamma_{18}$ we have
$\tau_{+}(0)(v_{1})<\tau_{+}(0)(v_{2})$. So,
$0<\tau_{+}(0)(v_{1})<\tau_{+}(0)(v_{2})<\tau_{+}(1)(v_{1})\ldots\tau_{+}(k)(v_{2})$
and we have IS upto smallest positive $k$ such that
$\upmu_{k}=\tau_{+}(k)(v_{2})-\tau_{+}(k)(v_{1})>\delta_{1}$ where
$\delta_{1}=\tau_{+}(n+1)(v_{1})-\tau_{+}(n)(v_{1})=\frac{2\pi}{v_{1}}$. On
the left side of $\Gamma_{18}$, we have
$\frac{\arccos(\frac{v_{1}\cos(\frac{\alpha\pi}{2})+cv_{1}\cos(\alpha\pi)-a}{b})}{v_{1}^{1/\alpha}}>\frac{\arccos(\frac{v_{2}\cos(\frac{\alpha\pi}{2})+cv_{2}\cos(\alpha\pi)-a}{b})}{v_{2}^{1/\alpha}}$
and system (8) is unstable for all $\tau\geq 0$ (See Figure (20)).

Using all this analysis we provide the stability region for all the cases viz.
$b>0$, $0<\alpha<1/2$, $b<0$, $0<\alpha<1/2$, $b>0$, $1/2<\alpha<1$ and $b<0$,
$1/2<\alpha<1$.

Figure 17: Graph for $b$ positive and $0<\alpha<1/2$ in complete
$a_{1}c-$plane  Figure 18: Graph for $b$ negative and $0<\alpha<1/2$ in
complete $a_{1}c-$plane  Figure 19: Graph for $b$ positive and $1/2<\alpha<1$
in complete $a_{1}c-$plane  Figure 20: Graph for $b$ negative and
$1/2<\alpha<1$ in complete $a_{1}c-$plane

##  IX Examples for the FDDE (8)

In this subsection, we verify the results derived in the manuscript and
presented in Figures (17), (18), (19) and (20). We utilize the numerical
methods proposed in bhalekar2019analysis  to solve these examples.
In the regions S and U of Figures (17)-(20), we verified the results by taking
a wide range of parameter values. As a representative example, we include a
set of parameter values in these regions and provide the solution curve. The
details are given in the Following tables (1) to (4).

Parameter values |  Nature of solution |  Figure
---|---|---
$a=0.45,b=1,c=4,\alpha=0.3,\tau=1.1$ (First quadrant in Figure (17)) |  Unstable |  Figure (21)(a)
$a=-9,b=4,c=4,\alpha=0.4,\tau=0.17$(Second quadrant in Figure (18)) |  Stable |  Figure (21)(b)
$a=-8.6,b=1.1,c=-4,\alpha=0.49,\tau=0.17$ (Third quadrant in Figure (17)) |  Unstable |  Figure (21)(c)
$a=-1,b=1.3,c=-0.8,\alpha=0.35,\tau=1.7$ (In the fourth quadrant region bounded by $c<0$ and the curve $\Gamma_{2}$ in Figure (17) ) |  Unstable |  Figure (21)(d)
$a=3.8,b=1.2,c=-4,\alpha=0.3,\tau=1$ (Region on the left side of $\Gamma_{6}$) |  Stable |  Figure(21)(e)
Table 1: System (8) $b>0$ and $0<\alpha<1/2$

(a)

(b)

(c)

(d)

(e)

Figure 21:  Figures representing Table (1) Parameter values |  Nature of solution |  Figure
---|---|---
$a=4.5,b=-3,c=9,\alpha=0.45,\tau=0.3$ (First quadrant in Figure (18)) |  Unstable |  Figure (22)(a)
$a=-8,b=-3.5,c=10,\alpha=0.29,\tau=0.8$ ((When $a_{1}<2b<0$ in Figure (18)) |  Stable |  Figure (22)(b)
$a=-1.1,b=-3.5,c=-1,\alpha=0.38,\tau=1$ (Third quadrant in Figure (19)) |  Unstable |  Figure (22)(c)
$a=0.3,b=-0.1,c=-0.2,\alpha=0.38,\tau=0.3$ (Region bounded by $c<0$ and $\Gamma_{13}$ in Figure (18)) |  Unstable |  Figure (22)(d)
$a=7,b=-0.1,c=-2,\alpha=0.38,\tau=0.3$ (Region on the right side of $\Gamma_{6}$ ) |  Stable |  Figure (22)(e)
Table 2: When $b<0$ and $0<\alpha<1/2$

(a)

(b)

(c)

(d)

(e)

Figure 22: Table (2) for $b<0$ and $0<\alpha<1/2$ Parameter values |  Nature of solution |  Figure
---|---|---
$a=1.06,b=8,c=4,\alpha=0.8,\tau=0.09$ (First quadrant in Figure (19)) |  Unstable |  Figure (23(a))
$a=-8,b=3.5,c=4,\alpha=0.75,\tau=0.6$ (Region on the left side of $\Gamma_{7}$ in Figure (19)) |  Stable |  Figure (23 (b))
$a=-8,b=5,c=-9,\alpha=0.93,\tau=0.9$ (Unstable region in the third quadrant in Figure (19)) |  Unstable |  Figure (23 (c))
$a=-4.9,b=5,c=-2,\alpha=0.96,\tau=0.9$ (Region on the left side of $\Gamma_{2}$ in Figure (19)) |  Unstable |  Figure (23 (d))
$a=200,b=5,c=-2,\alpha=0.96,\tau=2.3$ (Region on the right side of $\Gamma_{8}$ in Figure (19)) |  Stable |  Figure (23 (e))
Table 3: System (8) when $b>0$ and $1/2<\alpha<1$

(a)

(b)

(c)

(d)

(e)

Figure 23: Table for $b>0$ and $1/2<\alpha<1$ Parameter values |  Nature of solution |  Figure
---|---|---
$a=4,b=-3.3,c=19,\alpha=0.85,\tau=3.5$ (First quadrant in Figure (20)) |  Unstable |  Figure (24 (a))
$a=-7,b=-3.3,c=19,\alpha=0.85,\tau=1.8$ (Second quadrant in Figure (20)) |  Stable |  Figure (24 (b))
$a=-7,b=-3.3,c=-5,\alpha=0.79,\tau=1.3$ (Third quadrant in Figure (20)) |  Unstable |  Figure (24 (c))
$a=3.002,b=-3,c=-1,\alpha=0.7,\tau=1.3$ (Region on the left side of $\Gamma_{18}$ in Figure (20)) |  Unstable |  Figure (24 (d))
$a=10,b=-3,c=-2,\alpha=0.7,\tau=1.7$ (Region on the right side of $\Gamma_{10}$ in Figure (20)) |  Stable |  Figure (24 (e))
Table 4: System (8) when $b<0$ and $1/2<\alpha<1$

(a)

(b)

(c)

(d)

(e)

Figure 24: Table for $b<0$ and $1/2<\alpha<1$

Now, Let us consider the examples in each cases where the stability depend on
delay parameter.

######  Example IX.1.

Consider $b>0$ and $0<\alpha<1/2$ (Figure (17)). Only region that depends on
$\tau$ is the region bounded by $\Gamma_{2}$ and $\Gamma_{5}$. Let us fix
$c=-0.4$, $\alpha=0.3$ and $b=1$. At $\Gamma_{2}$, $\Gamma_{11}$ and
$\Gamma_{5}$, the bifurcation values of $a_{1}$ are $0.78726$, $2.78762$ and
$2.8219$ or $a=-0.21274$, $1.78762$ and $1.8219$ respectively. So, the system
(8) will be unstable for $a<-0.21274$. We will have SSR for
$-0.21274<a<1.78762$ and SS for $1.78762<a<1.8219$ from Figure (17). System
becomes stable $\forall\tau\geq 0$ if $a>1.8219$. The unstable solution for
$a=-0.7$ and $\tau=0.1$ is given in Figure (25)(a). For $a=1.7$ we get two
positive values of $v^{\alpha}$ as $1.12581$ and $2.12454$. Corresponding to
$v^{\alpha}=1.12581$ the critical values of delay are
$2.18454,6.41737,10.6502,\ldots$. Furthermore, $v^{\alpha}=2.12454$ gives
critical $\tau$ as $0.212729,0.722411,1.23209,\ldots$. The SSR is given by
$\tau\in[0,0.212729)$. Note that, at the critical values $0.722411$ and
$1.23209$, $Re\Big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\Big{)}>0$ as described in
Section (VIII). Therefore, the system (8) remains unstable for
$\tau>0.212729$. The stable solution for $\tau=0.1$ is shown in Figure (25)(b)
and the unstable solution for $\tau=0.3$ is given in Figure (25)(c). Now let
us consider $a=1.8$ in the SS region. For this parameter the two positive
values of $v^{\alpha}$ are $1.46349$ and $1.88464$. Corresponding to
$v^{\alpha}=1.46349$, we get the critical values of $\tau$ as
$0.874706,2.64025,4.40578,6.17132,\ldots$ and corresponding to
$v^{\alpha}=1.88464$ we get $0.343884,1.10378,1.86368,\ldots$. We have
$Re\Big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\Big{)}>0$ at $v^{\alpha}=1.88464$ and
$Re\Big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\Big{)}<0$ at $v^{\alpha}=1.46349$.
So, the characteristics roots will shift from left to right half plane at
$\tau=0.343884,1.10378,1.86368,\ldots$ and those roots will again shift back
at $\tau=0.874706,2.64025,4.40578,6.17132,\ldots$. We get stable solution for
$\tau=0.3$(cf. Figure (25)(d)), unstable solution for $\tau=0.4$(cf. Figure
(25)(e)) and again stable solution for $\tau=1$(cf. Figure (25)(f)).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 25:

######  Example IX.2.

Consider the case $b<0$ and $0<\alpha<1/2$ (Figure (18)). Let us verify the
region $c>0$ and $2b<a_{1}<0$. We have $a=-0.3$, $b=-10$, $\alpha=0.3$ and
$c=0.5$. We get only one positive value of $v^{\alpha}$ as $3.59088$ and
corresponding to that $v^{\alpha}$ we have
$Re\Big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\Big{)}>0$. So, the smallest critical
value of delay from equation (31) is $\tau=0.0336727$. The stable and unstable
solution for $\tau=0.02$ and $\tau=0.04$ are given in Figure (26)(a) and
(26)(b), respectively.

(a)

(b)

Figure 26: Figure of Example (IX.2)

######  Example IX.3.

Now let us consider the fourth quadrant of Figure (18) that depend on delay
$\tau$. When $b=-1$, $\alpha=0.45$ we have $c_{0}=-0.195086$.

  * ❖

So, when $c=-0.3<c_{0}$ then at bifurcation curve $\Gamma_{13}$,
$a_{1}=0.674322$ ($a=1.674322$), at $\Gamma_{2}$, $a_{1}=1.44121$
($a=2.44121$), at $\Gamma_{12}$, $a_{1}=1.50111$ ($a=2.50111$) and at
$\Gamma_{6}$, $a_{1}=1.62678$ ($a=2.62678$).

    * $\bullet$

If $a=2.35$ the two positive values of $v^{\alpha}$ are $2.03726$ and
$3.1629$. We have $Re\Big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\Big{)}>0$ at
$v^{\alpha}=3.1629$ and $Re\Big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\Big{)}<0$ at
$v^{\alpha}=2.03726$. The critical values of $\tau$ corresponding to $2.03726$
are $0.0192191,1.31166,2.60411,\ldots$ and for $3.1629$ we have
$0.397771,0.884042,1.37031,\ldots$. So, we get unstable solution for
$\tau\in[0,0.0192191)$, stable solution for $\tau\in(0.0192191,0.397771)$ and
it will remain unstable $\forall\tau>0.397771$. The unstable solution for
$\tau=0.01$ (Figure (27)(a)), stable solution for $\tau=0.02$ (Figure (27)(b))
and again unstable solution for $\tau=0.5$ (Figure (27)(c)) (Instability
switch).

    * $\bullet$

When $a=2.48$, we are in the SSR region from Figure (18). We get positive
values of $v^{\alpha}$ as $2.26432$ and $3.08213$. Note that we have
$Re\Big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\Big{)}<0$ at $v^{\alpha}=2.26432$ and
$Re\Big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\Big{)}>0$ at $v^{\alpha}=3.08213$.
Critical values of $\tau$ are $1.01403,2.03598,3.05793,\ldots$ corresponding
to $2.264320$ and $0.450092,0.965138,1.48018,\ldots$ corresponding to
$3.08213$. So, we have stable solution for $\tau\in[0,0.450092)$ and unstable
for all $\tau>0.450092$. The stable solution for $\tau=0.3$ is given in Figure
(27)(d) and unstable solution for $\tau=0.5$ is given in Figure (27)(e).

    * $\bullet$

Now, let us take $a=2.6$ so that we are in the SS region from Figure (18). So,
by solving equation (13) for $v^{\alpha}$ we get $2.56269$ and $2.91154$. We
have $Re\Big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\Big{)}<0$ at
$v^{\alpha}=2.56269$ and $Re\Big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\Big{)}>0$ at
$v^{\alpha}=2.91154$. The critical values of $\tau$ are
$0.740921,1.51711,2.2933,\ldots$ corresponding to $2.56269$ and $0.522197,\\\
1.10671,1.69123,\ldots$ corresponding to $2.91154$. We get a stability switch.
The stable solution for $\tau=0.4$ is given in Figure (27)(f), unstable
solution for $\tau=0.6$ is given in Figure (27)(g) and stable solution for
$\tau=1$ is given in Figure (27)(h).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 27:

######  Example IX.4.

Consider $b>0$ and $1/2<\alpha<1$ in Figure (19). Note that when $c>0$ the
region bounded by $\Gamma_{7}$ and the vertical axis with $c>0$ is the only
region that depends on delay $\tau$. Let us fix $b=1$ and $\alpha=0.8$ then
$c_{1}=2.75575$ and $c_{2}=2.52097$. If $c=4>c_{1}$ then at the bifurcation
curve $\Gamma_{7}$ we have $a_{1}=-0.0726$ or $a=-1.0726$ and at
$\Gamma_{16}$, $a_{1}=-0.0418347$ or $a=-1.0418347$. Now, if we take $a=-1.06$
we get two positive values of $v^{\alpha}$ as $0.28028999$ and $0.385268362$.
Note that $Re\Big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\Big{)}<0$ at $0.28028999$
and $Re\Big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\Big{)}>0$ at $0.385268362$.
Corresponding to $0.28028999$ the critical values of $\tau$ are
$28.5127,59.3212,90.1297,120.938,\ldots$ and corresponding to $0.385268362$
critical values of $\tau$ are $18.0739,38.7741,59.4744,80.1746,\ldots$. So, we
get stability switches here i.e. for $\tau\in[0,18.0739)$ we get stable
solution (cf. Figure (28)(a) for $\tau=17$), $\tau\in(18.0739,28.5127)$ gives
unstable solution (cf. Figure (28)(b) for $\tau=20$),
$\tau\in(28.5127,38.7741)$ gives stable solution (cf. Figure (28)(c)) for
$\tau=34$) and so on. Now, if we take $a=-1.02$, we get two positive value of
$v^{\alpha}$ as $0.18462087812$ and $0.420968$ and
$Re\Big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\Big{)}<0$ at $0.18462087812$ and
$Re\Big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\Big{)}>0$ at $0.420968$.
Corresponding to $0.18462087812$ we get critical values of $\tau$ are
$49.7825,101.702,153.621,205.54,257.46,\ldots$ and corresponding to $0.420968$
we get critical $\tau$ are $15.7097,34.2394,52.7691,71.2988,\ldots$. This
gives SSR region. So, we get stable solution for $\tau\in[0,15.7097)$ and it
will remain unstable $\forall\tau>15.7097$.

(a)

(b)

(c)

Figure 28:

######  Example IX.5.

Consider $c<0$ in Figure (19). Note that the value of $c_{3}$ is $-0.751566$
for $\alpha=0.8$ and $b=1$. Let us fix $c=-0.4>c_{3}$ then at $\Gamma_{2}$,
$a_{1}=6.54508$ ($a=5.54508$) and at $\Gamma_{8}$, $a_{1}=10.5099$ or
$a=9.5099$. If we choose $a=6$ then we are in the region bounded by the curves
$\Gamma_{2}$ and $\Gamma_{8}$ from Figure (19). With these parameters, we get
two positive values of $v^{\alpha}$ viz. $3.52068$ and $4.19432$.
Corresponding to $v^{\alpha}=3.52068$, we have
$Re\Big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\Big{)}<0$ and corresponding to
$v^{\alpha}=4.19432$ we have
$Re\Big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\Big{)}>0$. The critical values of
$\tau$ for $3.52068$ are $0.744534,2.04739,3.35024,\ldots$ and for $4.19432$
these are $0.0246074,1.07138,2.11816,\ldots$. We get stability switch
$[0,0.0246074)(S)$; $(0.0246074,0.744534)(U)$; $(0.744534,1.07138)(S)$,
$\ldots$. The stable solution for $\tau=0.01$ is given in Figure (29)(a). The
unstable solution for $\tau=0.03$ is given in Figure (29)(b). Figure (29)(c)
shows stable solution for $\tau=0.9$. So, we are in the SS region as shown in
Figure (19).

(a)

(b)

(c)

Figure 29:

######  Example IX.6.

Let us take $c>0$ in Figure (20). We have two bifurcation values of $c$ viz.
$c_{7}$ and $c_{5}$. Let us deal with each cases. When $b=-1$ and $\alpha=0.8$
we have $c_{7}=2.52097$ and $c_{5}=2.00571$.

  * ❖

When $c=9>c_{7}$ then at the bifurcation curve $\Gamma_{9}$ we have
$a_{1}=-2.2168$ or $a=-1.2168$ and at $\Gamma_{16}$, $a_{1}=-2.08023$ or
$a=-1.08023$.

    * $\bullet$

When $a<-1.2168$ then no positive value of $v$ exists and since system is
stable at $\tau=0$ it will remain stable $\tau\geq 0$.

    * $\bullet$

When $-1.2168<a<-1.08023$ then we have SS from Figure (20). So, let us take
$a=-1.2166$ and the two positive values of $v^{\alpha}$ are $v_{1}=0.28027$
and $v_{2}=0.286791$. Note that
$Re\big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\big{)}>0$ at $v_{2}$ and
$Re\big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\big{)}<0$ at $v_{1}$. Critical values
of delay corresponding to these $v^{\alpha}$ are given in the table (5). So,
we get stable solution for $\tau\in[0,11.2217)$ (cf. Figure (30)(a) for
$\tau=10$), unstable solution for $\tau\in(11.2217,11.7248)$ (cf. Figure
(30)(b) for $\tau=11.4$) and stable solution for $\tau\in(11.7248,41.1597)$
(cf. Figure (30)(c) for $\tau=20$) so on.

    * $\bullet$

Now, if we take $a>-1.08023$ then we are in SSR region. So, let us take
$a=-1.07$ we get two values of $v^{\alpha}$ as $v_{1}=0.13912$ and
$v_{2}=0.342681$. $Re\big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\big{)}<0$ at
$v_{1}$ and $Re\big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\big{)}>0$ at $v_{2}$. The
critical values corresponding to $v_{1}$ and $v_{2}$ are given in Table (6).
So, we get stable region for $\tau\in[0,7.23693)$ and unstable solution for
$\tau>7.23693$. Figure (30)(d) shows stable solution for $\tau=7$ and Figure
(30)(e) shows unstable solution for $\tau=9$.

  * ❖

Now, let us take $c=2.2$ which is in between $c_{7}$ and $c_{5}$. We have
$a_{1}=-2$ at $a_{1}=2b$ or $a=-1$, at $\Gamma_{9}$, $a_{1}=-1.9788$ or
$a=-0.9788$ and at $\Gamma_{16}$, $a_{1}=-1.98325$ or $a=-0.98325$.

    * $\bullet$

Let $a=-0.99$ we are in the SSR region from Figure (20) and there exist only
one positive value of $v^{\alpha}=0.0384923$ and
$Re\big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\big{)}>0$ at $v^{\alpha}$. The
critical values of $\tau$ corresponding to $v^{\alpha}$ are $182.001$,
$550.522$, $919.044$, $1287.56$ etc. So, we have the stable solution of system
(8) for $\tau\in[0,182.001)$. The stable solution for $\tau=9$ is shown in
Figure (30)(e) and one of the complex root with positive real part for
$\tau=184$ is given by $3.79477*10^{-7}+0.0168663I$. We know that one positive
root is sufficient for the instability.

    * $\bullet$

Now, let us take $a=-0.98$ then from Figure (20) we are in SS region. We get
three positive value of $v^{\alpha}$ namely $v_{1}=0.114518$, $v_{2}=0.208169$
and $v_{3}=0.357206$. Note that
$Re\big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\big{)}>0$ at $v_{1}$,
$Re\big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\big{)}<0$ at $v_{2}$ and
$Re\big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\big{)}>0$ at $v_{3}$. So, the
critical values of $\tau$ corresponding to each $v_{1}$, $v_{2}$ and $v_{3}$
are given in Table (7). We get stability switch $[0,9.46047)(S)$;
$(9.46047,20.5158)(U)$; $(20.5158,32.2131)(S)$, $\ldots$.

  * ❖

If $c=1<c_{5}$ then we have only one bifurcation value of $a_{1}=-2$, $(a=-1)$
from figure (20). So,

    * $\bullet$

if we take $a=-1.1$ then no positive $v^{\alpha}$ exist and system become
stable $\forall\tau\geq 0$

    * $\bullet$

if we take $a=-0.9$ then only one positive $v^{\alpha}=0.364932$ exist and
critical values of $\tau$ are $9.52706,31.6791,53.8312,\ldots$. We get stable
solution for $\tau\in[0,9.52706)$ (cf. Figure (31)(a) for $\tau=9$) and
unstable for $\tau>9.52706$(cf. Figure (31)(b) for $\tau=10$).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 30:  $v_{1}=0.28027$ |  $11.7248$ |  $42.536$ |  $73.3473$ |  $104.159$
---|---|---|---|---
$v_{2}=0.286791$ |  $11.2217$ |  $41.1597$ |  $71.0978$ |  $101.036$
Table 5: value of $v^{\alpha}$ and critical values of delay $\tau$ $v^{\alpha}=0.13912$ |  $34.1872$ |  $108.138$ |  $182.089$ |  $256.04$
---|---|---|---|---
$v^{\alpha}=0.342681$ |  $7.23693$ |  $31.2013$ |  $55.1657$ |  $79.1301$
Table 6: Value of $v^{\alpha}$ and critical values of delay $\tau$ $v_{1}=0.114518$ |  $45.264$ |  $139.581$ |  $233.898$ |  $328.215$
---|---|---|---|---
$v_{2}=0.208169$ |  $20.5158$ |  $65.2005$ |  $109.885$ |  $154.57$
$v_{3}=0.357206$ |  $9.46047$ |  $32.2131$ |  $54.9656$ |  $77.7182$
Table 7: Value of $v^{\alpha}$ and critical values of delay $\tau$

(a)

(b)

Figure 31:

######  Example IX.7.

Let us take $c<0$. The region bounded between $\Gamma_{18}$ and $\Gamma_{10}$
depends on delay $\tau$ (cf. Figure (20)). If we fix $b=-1$ and $\alpha=0.8$
then the intersection of $\Gamma_{17}$ and $\Gamma_{2}$ is $c_{6}=-1.5598$.

  * ❖

Let us first take $c=-0.9>c_{6}$ and then at $\Gamma_{18}$, $a_{1}=0.0220747$
or $a=1.0220747$ at $\Gamma_{2}$, $a_{1}=2.90893$ or $a=3.90893$ and at
$\Gamma_{10}$, $a_{1}=4.695$ or $a=5.695$.

    * $\bullet$

If we fix $a=1.03$ we get two positive value of $v^{\alpha}$ as $0.0889235$
and $1.45654$. The critical values of delay corresponding to $0.088923516$ are
$1.65726,131.05,260.442,\ldots$ and corresponding to $1.4565436$ we get
$1.79705,5.72373,9.65041,\ldots$. So, we have IS from Figure (20) i.e
$\tau\in[0,1.65726)(U)$; $(1.65726,1.79705)(S)$ and again unstable
$\forall\tau>1.79705$. The unstable solution for $\tau=1.5$ is given in Figure
(32)(a). The stable solution for $\tau=1.7$ is given in Figure (32)(b).
Unstable solution for $\tau=1.9$ shown in Figure (32)(c).

    * $\bullet$

If we take $a=4.1$ we are in the SS region from Figure (20). We get two
positive values of $v^{\alpha}$ as $1.86275$ and $2.36524$. The critical
values of delay corresponding to $1.86275$ and $2.36524$ are
$2.85785,5.74513,8.6324,\ldots$ and $1.34019,3.48227,5.62435,\ldots$
respectively. So, we get stable solution for $\tau\in[0,1.34019)$ (Figure
(32)(d) for $\tau=1.2$), unstable solution for $\tau\in(1.34019,2.85785)$
(Figure (32)(e) for $\tau=2.4$), again stable solution for
$\tau\in(2.85785,3.48227)$ (Figure (32)(f) for $\tau=3.3$) and so on.

    * –

Now for $a>5.695$ we get stable solution for all $\tau\geq 0$ because there
does not exist any positive value of $v^{\alpha}$ and system is stable at
$\tau=0$ so it will remain stable $\forall\tau\geq 0$.

  * ❖

Now let us take $c=-3<c_{6}$. At the bifurcation curve $\Gamma_{18}$,
$a_{1}=0.00178628$ or $a=1.00178628$, at $\Gamma_{2}$, $a_{1}=0.872678$ or
$a=1.872678$, at $\Gamma_{17}$, $a_{1}=1.24228$ or $a=0.24228$ and at
$\Gamma_{10}$ we have $a_{1}=2.34685$ or $a=3.34685$.

    * $\bullet$

If $a=1.8$ then we are in IS region from Figure (20). The two positive values
of $v^{\alpha}$ viz. $v_{1}=0.514075$ and $v_{2}=0.955756$. We have
$Re\big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\big{)}>0$ at $v_{2}$ and
$Re\big{(}\dfrac{d\lambda}{d\tau}|_{u=0}\big{)}<0$ at $v_{1}$. The critical
values of delay are given in Table (8). So, we have $[0,0.0526283)(U)$ (cf.
Figure (33)(a) for $\tau=0.04$); $(0.0526283,4.14761)(S)$ (cf. Figure (33)(b)
for $\tau=0.07$) and $\tau>4.14761$ we have unstable (cf. Figure (33)(c) for
$\tau=5$).

    * $\bullet$

If we take $a=2$ then from figure (20) we are in the SSR region. We get two
positive values of $v^{\alpha}$ as $0.581679$ and $0.981791$. The critical
values of $\tau$ are given in Table (9). So, we get SSR from
$\tau\in[0,4.10744)$. The stable solution for $\tau=3$ is shown in Figure
(34)(a) and unstable solution for $\tau=5$ is shown in Figure (34)(b).

    * $\bullet$

If we take $a=2.8$ then we are in SS region as shown in Figure (20). We get
two positive values of $v^{\alpha}$ as $0.820205$ and $1.05854$. The
corresponding critical values of $\tau$ are given in Table (10). The stable
solution for $\tau=3$, unstable solution for $\tau=5$ and stable solution for
$\tau=9$ are shown in Figures (33)(d), (e) and (f) respectively.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 32:  $v_{1}=0.514075$ |  $0.0526283$ |  $14.487$ |  $28.9213$ |  $43.3556$
---|---|---|---|---
$v_{2}=0.955756$ |  $4.14761$ |  $10.7965$ |  $17.4453$ |  $24.0941$
Table 8: Value of $v^{\alpha}$ and critical values of delay $\tau$ $v_{1}=0.581679$ |  $12.2832$ |  $24.652$ |  $37.0207$ |  $49.3894$
---|---|---|---|---
$v_{2}=0.981791$ |  $4.10744$ |  $10.5366$ |  $16.9658$ |  $23.395$
Table 9: Value of $v^{\alpha}$ and critical values of delay $\tau$ $v_{1}=0.820204$ |  $7.51375$ |  $15.5634$ |  $23.613$ |  $31.6627$
---|---|---|---|---
$v_{2}=1.058537$ |  $4.15684$ |  $10.0087$ |  $15.8607$ |  $21.7126$
Table 10: Value of $v^{\alpha}$ and critical values of delay $\tau$

(a)

(b)

(c)

(d)

(e)

(f)

Figure 33:

(a)

(b)

Figure 34:

##  X Conclusions

The transcendental nature of the characteristic equation
$\lambda^{2\alpha}x(t)+c\lambda^{\alpha}x(t)-a-be^{-\lambda\tau}=0$ makes the
two-term fractional-order delay differential equation a more complex than the
usual ODE or a one-term FDDE. We provided various conditions under which the
stability of the proposed equation does not depend on the delay. Further, we
used the hypothesis that the “change in stability can occur only when the
characteristic root crosses the imaginary axis in the complex plane". This
leads us to provide the boundary of the stable region in the parameter plane.
Furthermore, we worked on the stability switch and an instability switch too.
This made the work complete in all aspects. We hope this work will be very
useful to the scientists dealing with the systems involving memory and
hereditary properties. The conditions provided by us are very simple and
depend only on the parameter values.

## Acknowledgments

S. Bhalekar acknowledges the University of Hyderabad for Institute of
Eminence-Professional Development Fund (IoE-PDF) by MHRD (F11/9/2019-U3(A)).
D. Gupta thanks University Grants Commission for financial support
(No.F.82-44/2020(SA-III)).

## References

  * (1) G. F. Simmons, Differential equations with applications and historical notes.  CRC Press, 2016.
  * (2) H. K. Khalil, Control of nonlinear systems.  Prentice Hall, New York, NY, 2002.
  * (3) S. H. Strogatz, Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering.  CRC press, 2018.
  * (4) B. Van der Pol, “Lxxxviii. on “relaxation-oscillations”,” The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, vol. 2, no. 11, pp. 978–992, 1926.
  * (5) V. Volterra, Theory of Functionals and of Integral and Integro-differential Equations.  Dover Publications, INC, 1959.
  * (6) M. Lakshmanan and D. V. Senthilkumar, Dynamics of nonlinear time-delay systems.  Springer Science & Business Media, 2011.
  * (7) J. Hale, Functional differential equations, vol. 3.  Springer Science & Business Media, 2012.
  * (8) H. L. Smith, An introduction to delay differential equations with applications to the life sciences, vol. 57.  springer New York, 2011.
  * (9) K. Gopalsamy, Stability and oscillations in delay differential equations of population dynamics, vol. 74.  Springer Science & Business Media, 2013.
  * (10) T. Erneux, Applied delay differential equations.  Springer, 2009.
  * (11) I. Podlubny, Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications.  elsevier, 1998.
  * (12) A. Kilbas, Theory and applications of fractional differential equations.
  * (13) S. Bhalekar and V. Daftardar-Gejji, “A predictor-corrector scheme for solving nonlinear delay differential equations of fractional order,” J. Fract. Calc. Appl, vol. 1, no. 5, pp. 1–9, 2011.
  * (14) V. Daftardar-Gejji, Y. Sukale, and S. Bhalekar, “Solving fractional delay differential equations: a new approach,” Fractional Calculus and Applied Analysis, vol. 18, pp. 400–418, 2015.
  * (15) S. Bhalekar, “Stability and bifurcation analysis of a generalized scalar delay differential equation,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 26, no. 8, 2016.
  * (16) S. Bhalekar, V. Daftardar-Gejji, D. Baleanu, and R. Magin, “Fractional bloch equation with delay,” Computers & Mathematics with Applications, vol. 61, no. 5, pp. 1355–1365, 2011.
  * (17) S. Bhalekar and V. Daftardar-Gejji, “Fractional ordered liu system with time-delay,” Communications in Nonlinear Science and Numerical Simulation, vol. 15, no. 8, pp. 2178–2191, 2010.
  * (18) D. Baleanu, R. L. Magin, S. Bhalekar, and V. Daftardar-Gejji, “Chaos in the fractional order nonlinear bloch equation with delay,” Communications in Nonlinear Science and Numerical Simulation, vol. 25, no. 1-3, pp. 41–49, 2015\.
  * (19) V. Daftardar-Gejji, S. Bhalekar, and P. Gade, “Dynamics of fractional-ordered chen system with delay,” Pramana, vol. 79, pp. 61–69, 2012.
  * (20) S. Bhalekar, “Dynamical analysis of fractional order uçar prototype delayed system,” Signal, Image and Video Processing, vol. 6, pp. 513–519, 2012.
  * (21) S. B. Bhalekar, “Stability analysis of a class of fractional delay differential equations,” Pramana, vol. 81, no. 2, pp. 215–224, 2013.
  * (22) S. Bhalekar and D. Gupta, “Stability and bifurcation analysis of a fractional order delay differential equation involving cubic nonlinearity,” Chaos, Solitons & Fractals, vol. 162, p. 112483, 2022.
  * (23) W. Gu, Y. Yu, and W. Hu, “Artificial bee colony algorithm based parameter estimation of fractional-order chaotic system with time delay,” IEEE/CAA Journal of Automatica Sinica, vol. 4, no. 1, pp. 107–113, 2017.
  * (24) J. Shi, K. He, and H. Fang, “Chaos, hopf bifurcation and control of a fractional-order delay financial system,” Mathematics and Computers in Simulation, vol. 194, pp. 348–364, 2022.
  * (25) S. Bhalekar and D. Gupta, “Can a fractional order delay differential equation be chaotic whose integer-order counterpart is stable?,” in 2023 International Conference on Fractional Differentiation and Its Applications (ICFDA), pp. 1–6, IEEE, 2023.
  * (26) J. G. Milton and A. Longtin, “Evaluation of pupil constriction and dilation from cycling measurements,” Vision research, vol. 30, no. 4, pp. 515–525, 1990.
  * (27) S. A. Campbell, J. Belair, T. Ohira, and J. Milton, “Complex dynamics and multistability in a damped harmonic oscillator with delayed negative feedback,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 5, no. 4, pp. 640–645, 1995.
  * (28) T. Erneux, A. V. Kovalev, and E. A. Viktorov, “Short delay limit of the delayed duffing oscillator,” Physical Review E, vol. 108, no. 6, p. 064201, 2023.
  * (29) R. Vallée, P. Dubois, M. Coté, and C. Delisle, “Second-order differential-delay equation to describe a hybrid bistable device,” Physical Review A, vol. 36, no. 3, p. 1327, 1987.
  * (30) K. Diethelm, “The analysis of fractional differential equations,” 2004.
  * (31) J. Čermák and T. Kisela, “Stability properties of two-term fractional differential equations,” Nonlinear Dynamics, vol. 80, no. 4, pp. 1673–1684, 2015.
  * (32) S. Bhalekar, “Analysis of 2-term fractional-order delay differential equations,” in Fractional Calculus and Fractional Differential Equations, pp. 59–75, Springer, 2019.
  * (33) https://drive.google.com/drive/folders/147KhyNARmYlQhIt5caLDqGK9GHqvz2Yb?usp=sharing.