# Study of a Discretized Fractional-order eco-epidemiological model with prey
infection

Shuvojit Mondal11footnotemark: 1 Xianbing Cao22footnotemark: 2 Nandadulal
Bairagi<EMAIL_ADDRESS>Department of Mathematics,
Rabindra Mahavidyalaya
Hooghly-712401, India.  Beijing Technology and Business University
Department of Mathematics, School of Science
Beijing, China.  Centre for Mathematical Biology and Ecology
Department of Mathematics, Jadavpur University
Kolkata-700032, India.

###### Abstract

In this paper, an attempt is made to understand the dynamics of a three-
dimensional discrete fractional-order eco-epidemiological model with Holling
type II functional response. We first discretize a fractional-order predator-
prey-parasite system with piecewise constant arguments and then explore the
system dynamics. Analytical conditions for the local stability of different
fixed points have been determined using the Jury criterion. Several examples
are given to substantiate the analytical results. Our analysis shows that
stability of the discrete fractional order system strongly depends on the
step-size and the fractional order. More specifically, the critical value of
the step-size, where the switching of stability occurs, decreases as the order
of the fractional derivative decreases. Simulation results explore that the
discrete fractional-order system may also exhibit complex dynamics, like
chaos, for higher step-size.

###### keywords:

Fractional order differential equation, Ecological model, Discrete fractional-
order system, Local stability, Bifurcation, Chaos

††journal:  cor1cor1footnotetext: A part of this work was done at BTBU,
Beijing, China when NB visited this university in February-March, 2019.

##  1 Introduction

The idea of fractional calculus has been known since the notion of fractional
derivative was firstly introduced by Leibniz in a letter dated 30th September,
1695, where half-order derivative was mentioned. Since then it was
subsequently developed by Euler, Lagrange, Laplace, Fourier, Abel, Liouville,
Riemann, Heaviside, Caputo along with many others [1]. Generally, fractional
calculus deals with the study of fractional order integral and derivative
operators over real or complex domains and their applications [2]. Specially
it is a generalization of classical differential and integral calculus of
integer order to arbitrary order. In recent times, fractional order
derivatives and fractional order differential equations have been applied in
several fields of science and engineering [3, 4, 5, 6, 7, 8, 9, 10], however,
initially it was treated as a topic of interest of pure mathematicians only
[11]. The reason is two-fold. First, fractional derivatives have an additional
degree of freedom over its integer order counterpart due to the additional
parameter that represents its order and more suitable for those systems having
higher order dynamics and complex nonlinear phenomena [12, 13]. Secondly and
more importantly, fractional order derivatives not only depend on the local
conditions but also on the history of the function [14] and, therefore,
fractional derivatives have become an efficient tool for those systems, where
the consideration of memory or hereditary properties of the function is
essential to represent the system, e.g., in the case of biological systems. In
the last two decades, fractional order calculus has found its many
applications in biological sciences [15, 16, 17, 18, 19, 20, 21, 22, 23].

Discrete system of a continuous time model has significant applications in
simulations and computations. Particularly, the population model can be better
described in discrete time and it is especially applicable in populations
having non-overlapping generations, e.g., insect populations. Such discrete
models are known to show richer dynamics compare to its continuous counterpart
[24, 25, 26, 27, 28, 29]. Thus, both from biological and dynamic point of
views, study of discrete system is important. Though the number of discrete
population models is numerous, the study of discrete fractional-order
population models is rarely observed in the literature. Recently, a
transformed method was employed to obtain solutions of a family of finite
fractional difference equations in [30]. The global and local existence of
solutions for the fractional difference equation are given in [31]. Numerical
solution was adopted to synchronize chaos in a discrete fractional logistic
map [32]. For an ecological system, Elsadany and Matouk [33] recently studied
a two-dimensional fractional-order Lotka-Volterra predator prey model with its
discretization and showed complex dynamics. The existence and uniqueness of
solutions in two-dimensional discrete fractional Lotka-Volterra model was
addressed in [34]. Selvam and Janagaraj [35] studied the dynamics of a two-
dimensional discrete fractional order predator-prey system with type II
functional response. A two-dimensional fractional-order predator-prey model
with saturated harvesting was studied in [36]. Local stability analysis and
bifurcation analysis have been done in a simple fractional-order SIR model by
Salman [37]. Similar analysis was done for a fractional-order discrete SI
epidemic model in [38]. Discretization has also been done in other fractional-
order biological and physical systems [39, 40, 41, 42, 43, 44, 45]. These
relatively simple models show that the fractional-order plays a critical role
in breaking the stability. By considering a simple two-dimensional epidemic
model, Raheem and Salman [46] have shown that the dynamic complexity depends
also on the step-size of the discrete system. Analytical complexity, however,
multiplies if the dimension of the system increases and/or strong nonliearity
is there. In this paper, we construct a discrete version of a highly nonlinear
three-dimensional fractional-order predator-prey-parasite (PPP) model and
reveal its dynamics. So far our knowledge goes, no study has been done on a
discrete fractional-order PPP system. We unveil the dynamic properties of this
highly nonlinear model both from analytical and numerical points of view. We
show analytically that the stability of different fixed points depends on both
the fractional order and step-size. It also shows complex dynamics like chaos
as shown by the continuous time integer order system for some parameter
values.

For this study, we assume the following continuous-time PPP model studied by
Chattopadhyay and Bairagi [47]:

$\displaystyle\frac{dX}{dt}$ $\displaystyle=$ $\displaystyle
rX\bigg{(}1-\frac{X+Y}{K}\bigg{)}-\lambda YX,$ $\displaystyle\frac{dY}{dt}$
$\displaystyle=$ $\displaystyle\lambda YX-\frac{mYZ}{a+Y}-\mu Y,$ (1)
$\displaystyle\frac{dZ}{dt}$ $\displaystyle=$ $\displaystyle\frac{\theta
YZ}{a+Y}-dZ.$

It says that a prey species is divided into two sub classes: susceptible prey
$(X)$ and infected prey $(Y)$ due to some micro-parasite infection.
Susceptible prey can give birth and grow logistically to its maximum value $K$
supported by the environment with intrinsic growth rate $r$. The infection
spreads by contact following mass action law with $\lambda$ as its force of
infection. Infected prey is unable to reproduce but share resources with
susceptible prey. Infected prey dies out at a rate $\mu$, which considers both
the natural and infection-related deaths. Predator population $Z$ consumes
infected preys only as they are weakened due to infection and can not escape
predation. Here $m$ is the maximum prey attack rate of predator and $a$ is the
half-saturation constant. The parameter $\theta$ ($0<\theta<1$) measures the
reproductive gain of predator and $d$ measures its natural mortality. All
parameters are assumed to be positive from biological point of view. Readers
are referred to [47] for more discussion about the model.

Table 1.1. Parameter descriptions and their units

Parameters | Description | Units
---|---|---
$X$ | Density of susceptible prey | Number per unit designated area
$Y$ | Density of infected prey | Number per unit designated area
$Z$ | Density of predator | Number per unit designated area
$r$ | Intrinsic growth rate of prey | Per day
$K$ | Carrying capacity | Number per unit designated area
$a$ | Half-saturation constant | Number per unit designated area
$\lambda$ | Force of infection | Per day
$m$ | Search rate for infected prey | Per day
$\mu$ | Total (natural + disease-induced) | Per day
| death rate of infected prey |
$\theta$ | Reproductive gain | Per day
$d$ | Death rate of predator | Per day
$\alpha$ | Fractional order | $0<\alpha<1$

Considering the fractional derivative in caputo sense, Mondal et al. [23]
extended the work of Chattopadhyay and Bairagi [47] and analyzed the following
fractional order PPP model (1):

${}^{c}_{0}D^{\alpha}_{t}X$ $\displaystyle=$ $\displaystyle
rX\bigg{(}1-\frac{X+Y}{K}\bigg{)}-\lambda YX,~{}X(0)>0,$
${}^{c}_{0}D^{\alpha}_{t}Y$ $\displaystyle=$ $\displaystyle\lambda
YX-\frac{mYZ}{a+Y}-\mu Y,~{}Y(0)>0,$ (2) ${}^{c}_{0}D^{\alpha}_{t}Z$
$\displaystyle=$ $\displaystyle\frac{\theta YZ}{a+Y}-dZ,~{}Z(0)>0,$

where ${}^{c}_{0}D^{\alpha}_{t}$ is the Caputo fractional derivative with
fractional order $\alpha$ $(0<\alpha<1)$. It is worth mentioning that the
systems (1) and (1) become identical for $\alpha=1$. They have shown that the
solutions of system (1) are positively invariant and uniformly bounded in
$R^{3}_{+}$ under some restrictions. Existence and uniqueness of solution of
system (1) have also been discussed. Local and global stability of all
biologically feasible equilibrium points of system (1) were also proved. Here
we discritize the model system (1) and explore its dynamics.

##  2 Discretized fractional-order model and its analysis

Following Elsadany and Matouk [33], discretization of model system (1) with
piecewise constant arguments can be done in the following manner:

${}^{c}_{0}D^{\alpha}_{t}X$ $\displaystyle=$ $\displaystyle
rX([t/s]s)\bigg{(}1-\frac{X([t/s]s)+Y([t/s]s)}{K}\bigg{)}-\lambda
X([t/s]s)Y([t/s]s),$ ${}^{c}_{0}D^{\alpha}_{t}Y$ $\displaystyle=$
$\displaystyle\lambda
X([t/s]s)Y([t/s]s)-\frac{mY([t/s]s)Z([t/s]s)}{(a+Y([t/s]s))}-\mu Y([t/s]s),$
${}^{c}_{0}D^{\alpha}_{t}Z$ $\displaystyle=$ $\displaystyle\frac{\theta
Y([t/s]s)Z([t/s]s)}{(a+Y([t/s]s))}-dZ([t/s]s),$

with initial condition $X(0)=X_{0}>0$, $Y(0)=Y_{0}>0,Z(0)=Z_{0}>0$ and $s(>0)$
is the step-size. Let $t\in[0,s)$, so that $t/s\in[0,1)$. In this case, we
have

${}^{c}_{0}D^{\alpha}_{t}X$ $\displaystyle=$ $\displaystyle
X_{0}\bigg{(}r(1-\frac{X_{0}+Y_{0}}{K})-\lambda Y_{0}\bigg{)},$
${}^{c}_{0}D^{\alpha}_{t}Y$ $\displaystyle=$ $\displaystyle
Y_{0}\bigg{(}\lambda X_{0}-\frac{mZ_{0}}{(a+Y_{0})}-\mu\bigg{)},$
${}^{c}_{0}D^{\alpha}_{t}Z$ $\displaystyle=$ $\displaystyle
Z_{0}\bigg{(}\frac{\theta Y_{0}}{(a+Y_{0})}-d\bigg{)},$

and the solution of this fractional differential equation can be written as

$\displaystyle X_{1}(t)$ $\displaystyle=$ $\displaystyle
X_{0}+J^{\alpha}_{0}\bigg{(}X_{0}\bigg{(}r(1-\frac{X_{0}+Y_{0}}{K})-\lambda
Y_{0}\bigg{)}\bigg{)}$ $\displaystyle=$ $\displaystyle
X_{0}+\frac{t^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}X_{0}\bigg{(}r(1-\frac{X_{0}+Y_{0}}{K})-\lambda
Y_{0}\bigg{)}\bigg{)},$ $\displaystyle Y_{1}(t)$ $\displaystyle=$
$\displaystyle Y_{0}+J^{\alpha}_{0}\bigg{(}Y_{0}\bigg{(}\lambda
X_{0}-\frac{mZ_{0}}{(a+Y_{0})}-\mu\bigg{)}\bigg{)}$ $\displaystyle=$
$\displaystyle
Y_{0}+\frac{t^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}Y_{0}\bigg{(}\lambda
X_{0}-\frac{mZ_{0}}{(a+Y_{0})}-\mu\bigg{)}\bigg{)},$ $\displaystyle Z_{1}(t)$
$\displaystyle=$ $\displaystyle
Z_{0}+J^{\alpha}_{0}\bigg{(}Z_{0}\bigg{(}\frac{\theta
Y_{0}}{(a+Y_{0})}-d\bigg{)}\bigg{)}$ $\displaystyle=$ $\displaystyle
Z_{0}+\frac{t^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}Z_{0}\bigg{(}\frac{\theta
Y_{0}}{(a+Y_{0})}-d\bigg{)}\bigg{)},$

where
$J^{\alpha}_{0}=\frac{1}{\Gamma(\alpha)}\int^{t}_{0}(t-\tau)^{\alpha-1}d\tau,\alpha>0$.

In the second step, we assume $t\in[s,2s)$ so that $t/s\in[1,2)$ and similarly
obtain

${}^{c}_{0}D^{\alpha}_{t}X$ $\displaystyle=$ $\displaystyle
X_{1}(s)\bigg{(}r(1-\frac{X_{1}(s)+Y_{1}(s)}{K})-\lambda Y_{1}(s)\bigg{)},$
${}^{c}_{0}D^{\alpha}_{t}Y$ $\displaystyle=$ $\displaystyle
Y_{1}(s)\bigg{(}\lambda X_{1}(s)-\frac{mZ_{1}(s)}{(a+Y_{1}(s))}-\mu\bigg{)},$
${}^{c}_{0}D^{\alpha}_{t}Z$ $\displaystyle=$ $\displaystyle
Z_{1}(s)\bigg{(}\frac{\theta Y_{1}(s)}{(a+Y_{1}(s))}-d\bigg{)}.$

The solution of this equation reads

$\displaystyle X_{2}(t)$ $\displaystyle=$ $\displaystyle
X_{1}(s)+J^{\alpha}_{s}\bigg{(}X_{1}(s)\bigg{(}r(1-\frac{X_{1}(s)+Y_{1}(s)}{K})-\lambda
Y_{1}(s)\bigg{)}\bigg{)}$ $\displaystyle=$ $\displaystyle
X_{1}(s)+\frac{(t-s)^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}X_{1}(s)\bigg{(}r(1-\frac{X_{1}(s)+Y_{1}(s)}{K})-\lambda
Y_{1}(s)\bigg{)}\bigg{)},$ $\displaystyle Y_{2}(t)$ $\displaystyle=$
$\displaystyle Y_{1}(s)+J^{\alpha}_{s}\bigg{(}Y_{1}(s)\bigg{(}\lambda
X_{1}(s)-\frac{mZ_{1}(s)}{(a+Y_{1}(s))}-\mu\bigg{)}\bigg{)}$ $\displaystyle=$
$\displaystyle
Y_{1}(s)+\frac{(t-s)^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}Y_{1}(s)\bigg{(}\lambda
X_{1}(s)-\frac{mZ_{1}(s)}{(a+Y_{1}(s))}-\mu\bigg{)}\bigg{)},$ $\displaystyle
Z_{2}(t)$ $\displaystyle=$ $\displaystyle
Z_{1}(s)+J^{\alpha}_{s}\bigg{(}Z_{1}(s)\bigg{(}\frac{\theta
Y_{1}(s)}{(a+Y_{1}(s))}-d\bigg{)}\bigg{)}$ $\displaystyle=$ $\displaystyle
Z_{1}(s)+\frac{(t-s)^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}Z_{1}(s)\bigg{(}\frac{\theta
Y_{1}(s)}{(a+Y_{1}(s))}-d\bigg{)}\bigg{)},$

where
$J^{\alpha}_{s}=\frac{1}{\Gamma(\alpha)}\int^{t}_{s}(t-\tau)^{\alpha-1}d\tau,\alpha>0$.
Repeating the discretization process $n$ times, we have

$\displaystyle X_{n+1}(t)$ $\displaystyle=$ $\displaystyle
X_{n}(ns)+\frac{(t-ns)^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}X_{n}(ns)\bigg{(}r(1-\frac{X_{n}(ns)+Y_{n}(ns)}{K})-\lambda
Y_{n}(ns)\bigg{)}\bigg{)},$ $\displaystyle Y_{n+1}(t)$ $\displaystyle=$
$\displaystyle
Y_{n}(ns)+\frac{(t-ns)^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}Y_{n}(ns)\bigg{(}\lambda
X_{n}(ns)-\frac{mZ_{n}(ns)}{(a+Y_{n}(ns))}-\mu\bigg{)}\bigg{)},$
$\displaystyle Z_{n+1}(t)$ $\displaystyle=$ $\displaystyle
Z_{n}(ns)+\frac{(t-ns)^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}Z_{n}(ns)\bigg{(}\frac{\theta
Y_{n}(ns)}{(a+Y_{n}(ns))}-d\bigg{)}\bigg{)},$

where $t\in[ns,(n+1)s)$.
Making $t\rightarrow(n+1)s$, the corresponding fractional discrete model of
the continuous fractional model (1) is obtained as

$\displaystyle X_{n+1}$ $\displaystyle=$ $\displaystyle
X_{n}+\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}X_{n}\bigg{(}r(1-\frac{X_{n}+Y_{n}}{K})-\lambda
Y_{n}\bigg{)}\bigg{)},$ $\displaystyle Y_{n+1}$ $\displaystyle=$
$\displaystyle
Y_{n}+\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}Y_{n}\bigg{(}\lambda
X_{n}-\frac{mZ_{n}}{(a+X_{n})}-\mu\bigg{)}\bigg{)},$ $\displaystyle Z_{n+1}$
$\displaystyle=$ $\displaystyle
Z_{n}+\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}Z_{n}\bigg{(}\frac{\theta
Y_{n}}{(a+Y_{n})}-d\bigg{)}\bigg{)}.$ (3)

For $\alpha=1$, one obtains the so called Euler forward discrete model as a
special case of this generalized discrete model.

###  2.1 Existence and stability of fixed points

In the following, we investigate the dynamics of the discretized fractional
order model (2). At the fixed point, we have $X_{n+1}=X_{n}$, $Y_{n+1}=Y_{n}$
and $Z_{n+1}=Z_{n}$. One can then easily compute that (2) has the same fixed
points as in the integer and fractional order systems (1) and (1). The system
(2) has four equilibrium points [23, 47], viz. $E_{0}=(0,0,0)$ as the trivial
equilibrium, $E_{1}=(K,0,0)$ as the axial equilibrium, $E_{2}=(X_{1},Y_{1},0)$
as the planar equilibrium, where $X_{1}=\frac{\mu}{\lambda}$ and
$Y_{1}=\frac{r(\lambda K-\mu)}{\lambda(r+\lambda K)}$ and
$E^{*}=(X^{*},Y^{*},Z^{*})$ as the interior equilibrium, where

$\displaystyle X^{*}=K-\bigg{(}1+\frac{\lambda
K}{r}\bigg{)}Y^{*},~{}~{}Y^{*}=\frac{ad}{\theta-d},~{}~{}Z^{*}=\frac{(a+Y^{*})(\lambda
X^{*}-\mu)}{m}.$ (4)

Note that the equilibria $E_{0}$ and $E_{1}$ always exist. The planar
equilibrium point $E_{2}$ exists if $R_{0}>1$, where $R_{0}=\frac{\lambda
K}{\mu}$. The interior equilibrium $E^{*}$ exists if $(i)$ $R_{0}>1$ and
$(ii)$ $\theta>\theta_{1}$, where $\theta_{1}=d+\frac{\lambda ad(r+\lambda
K)}{r(\lambda K-\mu)}$.

The Jacobian matrix of system (2) at any arbitrary fixed point $(X,Y,Z)$ reads

${J(X,Y,Z)}=\begin{pmatrix}a_{11}&a_{12}&0\\\ a_{21}&a_{22}&a_{23}\\\
0&a_{32}&a_{33}\end{pmatrix},$ (5)

where

$\displaystyle a_{11}$ $\displaystyle=$ $\displaystyle
1+\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}r(1-\frac{2X+Y}{K})-\lambda
Y\bigg{)},$ $\displaystyle a_{12}$ $\displaystyle=$
$\displaystyle-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}X(\lambda+\frac{r}{K})\bigg{)},$
$\displaystyle a_{21}$ $\displaystyle=$
$\displaystyle\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\lambda Y,$
$\displaystyle a_{22}$ $\displaystyle=$ $\displaystyle
1+\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}\lambda
X-\frac{mZ}{a+Y}-\mu\bigg{)}+\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{mYZ}{(a+Y)^{2}},$
$\displaystyle a_{23}$ $\displaystyle=$
$\displaystyle-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{mY}{a+Y},$
$\displaystyle a_{32}$ $\displaystyle=$
$\displaystyle\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{a\theta
Z}{(a+Y)^{2}},$ $\displaystyle a_{33}$ $\displaystyle=$ $\displaystyle
1+\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}\frac{\theta
Y}{a+Y}-d\bigg{)}.$

Let $\xi_{1}$, $\xi_{2}$ and $\xi_{3}$ be the eigenvalues of the Jacobian
matrix (5). Then we have the following definition and lemma.

Definition 1 [37, 45] A fixed point $(X,Y,Z)$ of system (2) is called a sink
which is locally asymptotically stable if $\mid\xi_{i}\mid<1$ $(i=1,2,3)$ for
all $i$. It is called a source which is unstable if $\mid\xi_{i}\mid>1$
$(i=1,2,3)$ for all $i$. It is called a saddle if any two eigenvalues follow
opposite inequality, and a non-hyperbolic fixed point if $\mid\xi_{i}\mid=1$
$(i=1,2,3)$ for at least one $i$.

Lemma 1 [48] Suppose the characteristic polynomial $F(\xi)$ of a second order
Jacobian matrix is given by $F(\xi)=\xi^{2}-a_{1}\xi+a_{2}$. Then the
solutions $\xi_{i},~{}i=1,2,$ of $F(\xi)=0$ satisfy $\mid\xi_{i}\mid<1$ for
$i=1,2,$ if the following conditions hold:

$(i)~{}a_{2}<1,~{}(ii)~{}1+a_{2}>\mid a_{1}\mid.$

Lemma 2 [48] Suppose the characteristic polynomial $p(\xi)$ of Jacobian matrix
(5) is given by $p(\xi)=\xi^{3}+a_{1}\xi^{2}+a_{2}\xi+a_{3}$. Then the
solutions $\xi_{i},~{}i=1,2,3,$ of $p(\xi)=0$ satisfy $\mid\xi_{i}\mid<1$ for
$i=1,2,3$ if the following conditions hold:
$(i)~{}p(1)=1+a_{1}+a_{2}+a_{3}>0,~{}(ii)~{}(-1)^{3}p(-1)=1-a_{1}+a_{2}-a_{3}>0,~{}\mbox{and}~{}(iii)~{}1-(a_{3})^{2}>\mid
a_{2}-a_{3}a_{1}\mid.$

###  2.2 Main result

We now prove our main theorem in relation to the stability of the discrete
fractional order system (2).

Theorem 1 (a) The fixed point $E_{0}$ is always unstable for any feasible
values of $\alpha$ and $s$.

(b) The fixed point $E_{1}$ is locally asymptotically stable for $0<\alpha\leq
1$ if $R_{0}=\frac{\lambda K}{\mu}<1$, $s<min\\{s_{2},s_{3},s_{4}\\}$, where
$s_{2}=\sqrt[\alpha]{\frac{2\alpha\Gamma(\alpha)}{d}}$,
$s_{3}=\sqrt[\alpha]{\frac{2\alpha\Gamma(\alpha)}{r}}$ and
$s_{4}=\sqrt[\alpha]{\frac{2\alpha\Gamma(\alpha)}{\mu-\lambda K}}$. When
$R_{0}>1$, the fixed point $E_{1}$ is always unstable for $0<\alpha\leq 1$ for
any step size $s$.

(c) The fixed point $E_{2}$ is locally asymptotically stable for $0<\alpha\leq
1$ if $R_{0}=\frac{\lambda K}{\mu}>1$, $s_{5}<s<min\\{s_{6},s_{7}\\}$ with
$d>d_{1}$, where $s_{5}=\sqrt[\alpha]{\frac{\alpha\Gamma(\alpha)}{\lambda
K-\mu}}$, $s_{6}=\sqrt[\alpha]{\frac{2\alpha\Gamma(\alpha)}{d-d_{1}}}$,
$s_{7}=\sqrt[2\alpha]{\frac{\lambda K(\alpha\Gamma(\alpha))^{2}}{\mu r(\lambda
K-\mu)}}$ and $d_{1}=\frac{\theta r(\lambda K-\mu)}{a\lambda(\lambda
K+r)+r(\lambda K-\mu)}$. When $d<d_{1}$, the fixed point $E_{2}$ is always
unstable for $0<\alpha\leq 1$.

(d) The fixed point $E^{*}$ is locally asymptotically stable for $0<\alpha\leq
1$ if $R_{0}>1$, $0<s<min\\{s_{8},s_{9}\\}$ with $\theta>\theta_{1}$, where
$s_{8}=\sqrt[\alpha]{\frac{2K\alpha\Gamma(\alpha)}{rX^{*}}}$,
$\theta_{1}=d+\frac{\lambda ad(r+\lambda K)}{r(\lambda K-\mu)}$ and $s_{9}$ is
such that (8) holds.

Proof At the fixed point $E_{0}$, the Jacobian matrix $J(E_{0})$ can be
obtained from (5) and it is given by

$J(E_{0})=\begin{pmatrix}1+\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}r&0&0\\\
0&1-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\mu&0\\\
0&0&1-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}d\end{pmatrix}.$

The eigenvalues are $\xi_{1}=1+\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}r$,
$\xi_{2}=1-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\mu$ and
$\xi_{3}=1-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}d$. Since $|\xi_{1}|>1$,
$E_{0}$ is always unstable for $0<\alpha\leq 1$ and any $s>0$.

The Jacobian matrix $J(E_{1})$ at the fixed point $E_{1}$ is evaluated as

$J(E_{1})=\begin{pmatrix}1-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}r&\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}(\lambda
K+r)&0\\\ 0&1+\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}(\lambda K-\mu)&0\\\
0&0&1-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}d\end{pmatrix}.$

The eigenvalues of $J(E_{1})$ are

$\xi_{1}=1-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}r,~{}\xi_{2}=1+\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}(\lambda
K-\mu),~{}\xi_{3}=1-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}d.$

As $0<\alpha\leq 1$, $|\xi_{1,2,3}|<1$ hold if

$s<min\bigg{\\{}\sqrt[\alpha]{\frac{2\alpha\Gamma(\alpha)}{r}},\sqrt[\alpha]{\frac{2\alpha\Gamma(\alpha)}{d}},\sqrt[\alpha]{\frac{2\alpha\Gamma(\alpha)}{\mu-\lambda
K}}\bigg{\\}},\mu>\lambda k.$

Therefore, $E_{1}$ is locally asymptotically stable for $0<\alpha\leq 1$ if
$R_{0}=\frac{\lambda K}{\mu}<1$, $s<min\\{s_{2},s_{3},s_{4}\\}$, where
$s_{2}=\sqrt[\alpha]{\frac{2\alpha\Gamma(\alpha)}{d}}$,
$s_{3}=\sqrt[\alpha]{\frac{2\alpha\Gamma(\alpha)}{r}}$ and
$s_{4}=\sqrt[\alpha]{\frac{2\alpha\Gamma(\alpha)}{\mu-\lambda K}}$. However,
if $R_{0}>1$ then $|\xi_{2}|>1$. So $E_{1}$ is always unstable for all
$0<\alpha\leq 1$ and for any step size $s$.

At the fixed point $E_{2}$, the Jacobian matrix $J(E_{2})$ takes the form

$J(E_{2})=\begin{pmatrix}1-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{r\mu}{\lambda
K}&-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{\mu(r+\lambda K)}{\lambda
K}&0\\\ \frac{r(\lambda K-\mu)}{(\lambda
K+r)}\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}&0&-\frac{mr(\lambda
K-\mu)}{a\lambda(\lambda K+r)+r(\lambda
K-\mu)}\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\\\
0&0&1+\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}\frac{\theta r(\lambda
K-\mu)}{a\lambda(\lambda K+r)+r(\lambda K-\mu)}-d\bigg{)}\end{pmatrix}.$

After some mathematical manipulations, one can have the following
characteristic equation

$\bigg{(}1+\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}\frac{\theta
r(\lambda K-\mu)}{a\lambda(\lambda K+r)+r(\lambda
K-\mu)}-d\bigg{)}-\xi\bigg{)}(\xi^{2}+A\xi+B)=0,$ (6)

where $A=\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{r\mu}{\lambda K}-1$ and
$B=\bigg{(}\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{)}^{2}\frac{r\mu(\lambda
K-\mu)}{\lambda K}$. Therefore, one eigenvalue is
$\xi_{1}=1+\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}(d_{1}-d)$, where
$d_{1}=\frac{\theta r(\lambda K-\mu)}{a\lambda(\lambda K+r)+r(\lambda K-\mu)}$
and the other two eigenvalues $\xi_{2,3}$ are the roots of $\xi^{2}+A\xi+B=0$.
Following Lemma $1$, we observe $\mid\xi_{2,3}\mid<1$ if
$\sqrt[\alpha]{\frac{\alpha\Gamma(\alpha)}{\lambda
K-\mu}}<s<\sqrt[2\alpha]{\frac{\lambda K(\alpha\Gamma(\alpha))^{2}}{\mu
r(\lambda K-\mu)}}$ with $R_{0}>1$. For $\xi_{1}$, we consider the following
two cases:

Case I: If $d>d_{1}$ then
$\xi_{1}=1-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}(d-d_{1})$ and
$\mid\xi_{1}\mid<1$ if
$0<s<\sqrt[\alpha]{\frac{2\alpha\Gamma(\alpha)}{d-d_{1}}}$. Therefore, the
equilibrium $E_{2}$ is locally asymptotically stable for $0<\alpha\leq 1$ if
$R_{0}=\frac{\lambda K}{\mu}>1$, $s_{5}<s<min\\{s_{6},s_{7}\\}$ with
$d>d_{1}$, where $s_{5}=\sqrt[\alpha]{\frac{\alpha\Gamma(\alpha)}{\lambda
K-\mu}}$, $s_{6}=\sqrt[\alpha]{\frac{2\alpha\Gamma(\alpha)}{d-d_{1}}}$,
$s_{7}=\sqrt[2\alpha]{\frac{\lambda K(\alpha\Gamma(\alpha))^{2}}{\mu r(\lambda
K-\mu)}}$ and $d_{1}=\frac{\theta r(\lambda K-\mu)}{a\lambda(\lambda
K+r)+r(\lambda K-\mu)}$.

Case II: If $d<d_{1}$ then
$\xi_{1}=1+\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}(d_{1}-d)>1$. So $E_{2}$ is
always unstable for all $0<\alpha\leq 1$ and for any step size $s$.

At the interior fixed point $E^{*}$, the Jacobian matrix $J(E^{*})$ is
evaluated as

${J(E^{*})}=\begin{pmatrix}1-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{rX^{*}}{K}&-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}X^{*}(\frac{r}{K}+\lambda)&0\\\
\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\lambda
Y^{*}&1+\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{mY^{*}Z^{*}}{(a+Y^{*})^{2}}&-\frac{mY^{*}}{a+Y^{*}}\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\\\
0&\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{\theta
aZ^{*}}{(a+Y^{*})^{2}}&1\end{pmatrix}.$

After some algebraic manipulation, the characteristic equation of $J(E^{*})$
can be expressed as

$p(\xi)=\xi^{3}+A_{1}\xi^{2}+A_{2}\xi+A_{3}=0,$ (7)

where

$\displaystyle A_{1}$ $\displaystyle=$
$\displaystyle\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{rX^{*}}{K}-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{mY^{*}Z^{*}}{(a+Y^{*})^{2}}-3,$
$\displaystyle A_{2}$ $\displaystyle=$ $\displaystyle
3+\frac{2s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}\frac{mY^{*}Z^{*}}{(a+Y^{*})^{2}}-\frac{rX^{*}}{K}\bigg{)}+\bigg{(}\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{)}^{2}\bigg{(}\frac{amdZ^{*}}{(a+Y^{*})^{2}}-\frac{rmX^{*}Y^{*}Z^{*}}{K(a+Y)^{2}}\bigg{)}$
$\displaystyle+$
$\displaystyle\bigg{(}\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{)}^{2}X^{*}Y^{*}\bigg{(}\frac{r\lambda}{K}+\lambda^{2}\bigg{)},$
$\displaystyle A_{3}$ $\displaystyle=$ $\displaystyle
1+\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}\frac{mY^{*}Z^{*}}{(a+Y^{*})^{2}}-\frac{rX^{*}}{K}\bigg{)}+\bigg{(}\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{)}^{2}\frac{amdZ^{*}}{(a+Y^{*})^{2}}\bigg{(}1-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{rX^{*}}{K}\bigg{)}$
$\displaystyle+$
$\displaystyle\bigg{(}\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{)}^{2}\bigg{(}\frac{rmX^{*}Y^{*}Z^{*}}{K(a+Y^{*})^{2}}+X^{*}Y^{*}\bigg{(}\frac{r\lambda}{K}+\lambda^{2}\bigg{)}\bigg{)}.$

We already know that the interior equilibrium $E^{*}$ exists if $(i)$
$R_{0}>1$ and $(ii)$ $\theta>\theta_{1}$, where $\theta_{1}=d+\frac{\lambda
ad(r+\lambda K)}{r(\lambda K-\mu)}$. Using these conditions, we prove the
conditions of Lemma $2$:

  * (i)

$p(1)=1+A_{1}+A_{2}+A_{3}\\\
=\bigg{(}\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{)}^{3}\bigg{(}\frac{adrmX^{*}Z^{*}}{K(a+Y^{*})^{2}}\bigg{)}$
$>0$ for any $s$.

  * (ii)

$(-1)^{3}p(-1)=1-A_{1}+A_{2}-A_{3}\\\
=4\bigg{(}2-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{rX^{*}}{K}\bigg{)}+\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{mZ^{*}}{(a+I^{*})^{2}}\bigg{(}\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}ad+2Y^{*}\bigg{)}\bigg{(}2-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{rX^{*}}{K}\bigg{)}\\\
+\bigg{(}\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{)}^{2}2X^{*}Y^{*}\bigg{(}\frac{r\lambda}{K}+\lambda^{2}\bigg{)}$
$>0$, if $s<\sqrt[\alpha]{\frac{2K\alpha\Gamma(\alpha)}{rX^{*}}}$.

  * (iii)

$(1-A_{3}^{2})-\mid
A_{2}-A_{3}A_{1}\mid=\bigg{[}1-\bigg{\\{}\bigg{(}1+\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{mZ^{*}}{(a+Y^{*})^{2}}\bigg{(}\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}ad+Y^{*}\bigg{)}\bigg{)}\bigg{(}1-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{rX^{*}}{K}\bigg{)}+\bigg{(}\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{)}^{2}X^{*}Y^{*}\bigg{(}\frac{r\lambda}{K}+\lambda^{2}\bigg{)}\bigg{\\}}^{2}\bigg{]}-\\\
\left|\bigg{\\{}3+\frac{2s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{(}\frac{mY^{*}Z^{*}}{(a+Y^{*})^{2}}-\frac{rX^{*}}{K}\bigg{)}+\bigg{(}\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{)}^{2}\bigg{(}\frac{amdZ^{*}}{(a+Y^{*})^{2}}-\frac{rmX^{*}Y^{*}Z^{*}}{K(a+Y^{*})^{2}}\bigg{)}\right.$
$\left.+\bigg{(}\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{)}^{2}X^{*}Y^{*}\bigg{(}\frac{r\lambda}{K}+\lambda^{2}\bigg{)}\bigg{\\}}-\bigg{(}\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{rX^{*}}{K}-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{mY^{*}Z^{*}}{(a+Y^{*})^{2}}-3\bigg{)}\right.$
$\left.\bigg{\\{}\bigg{(}1+\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{mZ^{*}}{(a+Y^{*})^{2}}\bigg{(}\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}ad+Y^{*}\bigg{)}\bigg{)}\bigg{(}1-\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\frac{rX^{*}}{K}\bigg{)}+\right.$
$\left.\bigg{(}\frac{s^{\alpha}}{\alpha\Gamma(\alpha)}\bigg{)}^{2}X^{*}Y^{*}\bigg{(}\frac{r\lambda}{K}+\lambda^{2}\bigg{)}\bigg{\\}}\right|$.

Due to the complexity of mathematical expression, it is difficult to find a
restriction on the step-size $s$ in terms of other system parameters so that
$(1-A_{3}^{2})-\mid A_{2}-A_{3}A_{1}\mid>0$. However, in the numerical
section, we show that there exists a value $s_{9}$ (say) of $s$ such that

$(1-A_{3}^{2})-\mid A_{2}-A_{3}A_{1}\mid>0,~{}\forall s<s_{9}.$ (8)

Thus, following Jury criteria, the fixed point $E^{*}$ is locally
asymptotically stable for $0<\alpha\leq 1$ if $R_{0}>1$, $\theta>\theta_{1}$
and $0<s<min\\{s_{8},s_{9}\\}$, where
$s_{8}=\sqrt[\alpha]{\frac{2K\alpha\Gamma(\alpha)}{rX^{*}}}$,
$\theta_{1}=d+\frac{\lambda ad(r+\lambda K)}{r(\lambda K-\mu)}$ and $s_{9}$ is
such that (8) holds. Hence the theorem.

##  3 Numerical Simulations

In this section, an extensive numerical simulations are performed for the
fractional-order discrete system (2) to validate our theoretical results for
different fractional values of $\alpha$ $(0<\alpha<1)$. Stability of fixed
points depends on the step-size $s$ and fractional-order $\alpha$ (see Theorem
1). To demonstrate the effect of fractional order $\alpha$ and step-size $s$
on the system dynamics, we cite different examples.

Table 3.1. Restrictions on the step-size, following Theorem 1(b,c), for the
stability of fixed points $E_{1}$ and $E_{2}$ for different fractional-orders
$\alpha$. For $E_{1}$, we consider $K=40$, $\lambda=0.005$, $\theta=0.189$ and
for $E_{2}$ we consider $K=200$, $\lambda=0.015$, $\theta=0.08$. Following
parameters are common for both the equilibrium points: $r=2.0$, $m=0.52$,
$\mu=0.28$, $a=15.0$, $d=0.09$.

| $E_{1}$ | $E_{2}$
---|---|---
Fractional order $\alpha$ |  Step-size $s<min(s_{2},s_{3},s_{4})$ |  Step-size $s_{5}<s<min(s_{6},s_{7})$
$\alpha=0.3$ | $s_{2}=21513$ | $s_{5}=0.0248$
| $s_{3}=0.6973$ | $s_{6}=1835600$
| $s_{4}=31857$ | $s_{7}=2.1578$
$\alpha=0.4$ | $s_{2}=1726.2$ | $s_{5}=0.0608$
| $s_{3}=0.7415$ | $s_{6}=48464$
| $s_{4}=2317.3$ | $s_{7}=1.7302$
$\alpha=0.6$ | $s_{2}=145.5959$ | $s_{5}=0.1564$
| $s_{3}=0.8289$ | $s_{6}=1344.9$
| $s_{4}=177.1754$ | $s_{7}=1.4582$
$\alpha=0.8$ | $s_{2}=44.1464$ | $s_{5}=0.2619$
| $s_{3}=0.9150$ | $s_{6}=233.9135$
| $s_{4}=51.1488$ | $s_{7}=1.3976$
$\alpha=0.95$ | $s_{2}=25.6082$ | $s_{5}=0.3414$
| $s_{3}=0.9788$ | $s_{6}=104.2794$
| $s_{4}=28.9884$ | $s_{7}=1.3984$

Example 1: We consider the following parameter values from [49]: $r=2.0$,
$K=40.0$, $\lambda=0.005$, $m=0.52$, $\mu=0.28$, $a=15.0$, $\theta=0.189$,
$d=0.09$ and initial point $S(0)=30,I(0)=5,P(0)=10$. With these parameter
values, we calculate different limits $s_{i}$ on the step-size for the
stability of fixed points $E_{1}$ and $E_{2}$ for different values of
fractional-order $\alpha$. For example, we compute $s_{2}=44.1464$,
$s_{3}=0.9150$ and $s_{4}=51.1488$ for $\alpha=0.8$. The equilibrium $E_{1}$
will then be locally asymptotically stable, following Theorem 1(b), if
$s<min(44.1464,0.9150,51.1488)=0.9150$. It is unstable if $s>0.9150$.
Similarly, $E_{2}$ will be locally asymptotically stable, following Theorem
1(c), if $0.2619<s<min(233.9135,1.3976)=1.3976$. Similar step-size
restrictions for the stability of the above two equilibrium points for
different fractional order $\alpha$ are tabulated in Table 3.1. Considering
$s=0.65(<0.9150)$ and $s=0.95(>0.9150)$ corresponding to $\alpha=0.8$, the
behavior of the solution trajectories of the system (1) are plotted in Fig. 1.
It shows that the predator-free equilibrium $E_{1}$ is stable for the step-
size $s=0.65$ and unstable for $s=0.95$.

Figure 1: Stable (Fig. a) and unstable (Fig. b) behavior of the fixed point
$E_{1}$ for the step-size $s=0.65$ and $s=0.95$ corresponding to the fraction
order $\alpha=0.8$. Parameter values are $r=2.0$, $K=40.0$, $\lambda=0.005$,
$m=0.52$, $\mu=0.28$, $a=15.0$, $\theta=0.189$, $d=0.09$.  Figure 2: Stable
(Fig. a) and unstable (Fig. b) behavior of the fixed point $E_{2}$ for the
step-size $s=0.85$ and $s=1.45$ corresponding to the fraction order
$\alpha=0.8$. Parameter values are $\lambda=0.005$, $K=200$, $\theta=0.08$.
Other parameters and initial values are as in Fig. 1.

Similarly, we observe (Figure 2) that the predator-free fixed point $E_{2}$ of
the system (1) is stable for $s=0.85$ and unstable for $s=1.45$ corresponding
to order $\alpha=0.8$ (see Theorem 1(c) and Table 3.1.).

Figure 3: Existence of the bound $s_{9}$ such that the condition
$1-A_{3}^{2}>$ $\mid A_{2}-A_{3}A_{1}\mid$ is hold. The bound is obtained as
$s_{9}=0.07798$ for the parameter values $r=15.0$, $K=40.0$, $\lambda=0.006$,
$m=14.5$, $\mu=0.0019$, $a=16.0$, $\theta=11.1$, $d=6.0$ and $\alpha=0.8$.

Table 3.2. Restriction on the step-size, following Theorem 1(d), for the
stability of fixed point $E^{*}$ for different fractional-order $\alpha$.
Parameter values are as in Fig. 3.

| $E^{*}$
---|---
Fractional order $\alpha$ |  Step-size $s<min(s_{8},s_{9})$
$\alpha=0.3$ | $s_{8}=0.0074$
| $s_{9}=0.00098$
$\alpha=0.4$ | $s_{8}=0.0245$
| $s_{9}=0.00538$
$\alpha=0.6$ | $s_{8}=0.0853$
| $s_{9}=0.03108$
$\alpha=0.8$ | $s_{8}=0.1662$
| $s_{9}=0.07798$
$\alpha=0.95$ | $s_{8}=0.2327$
| $s_{9}=0.1073$

Example 2: We consider another parameter set as in [47]: $r=15.0$, $K=40.0$,
$\lambda=0.006$, $m=14.5$, $\mu=0.0019$, $a=16.0$, $\theta=11.1$, $d=6.0$ and
initial point $S(0)=30,I(0)=5,P(0)=10$. First we numerically determine the
bound $s_{9}=0.07798$ (Fig. 3) for the stability of the interior fixed point
$E^{*}$ corresponding to a fractional order $\alpha=0.8$ so that the condition
$1-A_{3}^{2}>\mid A_{2}-A_{3}A_{1}\mid$ is satisfied (see Theorem 1(d)).

As before, we construct Table 3.2 to show the bounds on step-size for the
stability of the interior fixed point $E^{*}$ for different values of
fractional-order $\alpha$. One can also compute that $R_{0}=126.3158>1$ and
$\theta=11.1>\theta_{1}=8.4579$ for $\alpha=0.8$. Following Theorem 1(d), the
fixed point $E^{*}(20.8573,18.8235,$ $0.2962)$ is then stable (Fig. 4a) for
$s=0.05<min(s_{8},s_{9})=0.07798$ (see Table 3.2) and unstable for
$s=0.08(>0.07798)$ (Fig. 4b) for the fractional-order $\alpha=0.8$.

Figure 4: Stable and unstable behavior of the fixed point $E^{*}$,
respectively, for the step-sizes $s=0.05$ (Fig. a) and $s=0.08$ (Fig. b) when
$\alpha=0.8$. Here all parameters are as in Figure 3. Figure 5: Bifurcation
diagram of $X$ population of system (2) in the range $[0.055,0.085]$ with
step-size $s$ as the bifurcation parameter. System remains stable for
$s<s^{*}$ and becomes unstable for $s>s^{*}$, where $s^{*}=0.077$. Here
$\alpha=0.8$ and othere parameters are as in Figure 3.

Bifurcation diagram of $X$ population (Fig. 5) shows that population is stable
for $s<s^{*}$ and unstable for $s>s^{*}$, where $s^{*}=0.077$. These results
show that stability of the discrete fractional order system (1) strongly
depends on the step-size and the fractional order.

Example 3: In [50], the authors showed that the continuous system may show
chaotic dynamics for some parameter values. Here we consider the same
parameter set and initial value for which the system (1) shows chaotic
behavior: $r=22.0$, $K=300.0$, $\lambda=0.06$, $m=15.5$, $\mu=2.3$, $a=15.0$,
$\theta=10$, $d=8.3$ and $(S(0),I(0),P(0))=(30,5,10)$. As before, we determine
the bound $s_{9}=0.0150$ (see Fig. 6) for the stability of fixed point $E^{*}$
when $\alpha=0.85$ and tabulate (see Table 3.3) the bounds on step-size,
following Theorem 1(d), for different values of fractional-order $\alpha$.
Table 3.3 shows that the critical value of the step-size, where the switching
of stability occurs, decreases as the order of the fractional derivative
decreases. In Figure 7, we present the stable (Fig. a) and unstable (Fig. b)
behavior of the system around the interior fixed point $E^{*}=(166.8449,$
$73.2353,$ $43.8939)$ for the step-sizes $s=0.01$ and $s=0.04$, respectively,
when $\alpha=0.85$.

Figure 6: The bound $s_{9}=0.0150$ is determined such that the condition
$1-A_{3}^{2}>$ $\mid A_{2}-A_{3}A_{1}\mid$ holds. The parameter values are
$r=22.0$, $K=300.0$, $\lambda=0.06$, $m=15.5$, $\mu=2.3$, $a=15.0$,
$\theta=10$, $d=8.3$ and $\alpha=0.85$.

Table 3.3. Restriction on the step-size for the stability of fixed point
$E^{*}$ for different fractional-orders $\alpha$ with the parameters of Fig. 6
.

| $E^{*}$
---|---
Fractional order $\alpha$ |  Step-size $s<min(s_{8},s_{9})$
$\alpha=0.95$ |  $s_{8}=0.1455$, $s_{9}=0.0242$
$\alpha=0.85$ |  $s_{8}=0.1112$, $s_{9}=0.0150$
$\alpha=0.6$ |  $s_{8}=0.0405$, $s_{9}=0.0023$
$\alpha=0.55$ |  $s_{8}=0.0300$ , $s_{9}=0.0013$
$\alpha=0.45$ |  $s_{8}=0.0136$, $s_{9}=0.0003$

Figure 7: Stable and unstable behavior of the fixed point $E^{*}=(166.8449,$
$73.2353,$ $43.8939)$ for the step-sizes $s=0.01$ (Fig. a) and $s=0.04$ (Fig.
b) when $\alpha=0.85$. Here all parameters are as in Fig. 6. Figure 8:
Bifurcation diagrams of $X$ population of system (2) in the range
$[0.01,0.054]$ with respect to the step-size $s$. System is stable for
$s<s^{*}$ and unstable for $s>s^{*}$. Exchange of stability occurs at
$s^{*}=0.015$ (Fig. 8a). A magnified view of the bifurcation diagram for
$0.04\leq s\leq 0.054$ is presented in Fig. 8b. All parameters are as in Fig.
6 with $\alpha=0.85$.

Bifurcation diagram of system populations (Fig. 8a) shows that population is
stable for $s<s^{*}$ and losts its stability for $s>s^{*}$ corresponding to
fractional order $\alpha=0.85$, where $s^{*}=0.015$. As the step-size is
further increased, the system shows complex dynamics like chaos. An enlarged
view (Fig. 8b) of the bifurcation diagram for small range of higher step-size
is presented for visualisation of the chaotic dynamics. Existence of such
chaotic dynamics (see Fig. 9) may be verified by computing the largest
Lyapunov exponent [38].

Figure 9: Largest Lyapunov exponent (LLE) drawn for two values of $\alpha$.
The left figure shows that LLE becomes positive for $s>s^{*}$, where
$s^{*}=0.015$, $\alpha=0.85$ and other parameters are as in Fig. 8. The right
figure shows that LLE becomes positive for $s>s^{*}$, where $s^{*}=0.077$,
$\alpha=0.8$ and other parameters are as in Fig. 5.

##  4 Summary

The application of fractional calculus in the study of biological systems is
important because it considers past history of the state variables. On the
other hand, discretization of a continuous system is inevitable for numerical
computations and for population models, particularly where species generation
does not overlap. It is earlier shown that the stability of the discrete
fractional-order system depends on its fractional-order [33, 34, 35, 36, 37,
38]. It was shown for a simpler two-dimensional epidemic model that the
stability of a fractional-order discrete system also depends on the step-size
[46]. It is worth mentioning that the determination of explicit bounds would
be difficult if the nonlinearity and dimension of the system are high. The
novelty of this work is that we determined the analytical bounds on the step-
size for a highly nonlinear three-dimensional system so that the stability of
an equilibrium point is maintained. We considered a highly nonlinear three-
dimensional system of fractional-order differential equations that represent
the predator-prey interaction in the presence of infection. We revealed the
dynamics of the corresponding discrete fractional system. Discretization of
the fractional-order system was done with piecewise constant arguments and the
corresponding dynamics were explored around different feasible fixed points.
We used the Jury criterion for local stability of the discrete fractional-
order systems. It is observed that stability of the system depends on both the
step-size and fractional order. Different examples were cited to illustrate
the stability of predator-free, infection-free, and coexistence equilibrium
points. These simulation results perfectly agree with the analytical results.
More specifically, the critical value of the step-size, where the switching of
stability occurs, decreases as the order of the fractional derivative
decreases. Bifurcation diagrams with respect to step size show that the system
is stable and all populations coexist in the stable state if the step-size
does not exceed some threshold value, which is also dependent on the
fractional order, and unstable if it exceeds. Our simulation result also shows
that the system may show chaotic dynamics, as shown in the continuous-time
case, if the step-size is large. This study thus reveals that the dynamics of
a discrete fractional order system strongly depends on the step-size and the
fractional order.

## References

  * [1] S. Samko, A. Kilbas, O. Marichev, _Fractional Integrals and Derivatives: Theory and Applications_. Gordon and Breach Science Publishers, Yverdon, (1993).
  * [2] F.A. Rihan, S. Lakshmanan, A. H. Hashish, R. Rakkiyappan, E. Ahmed, Fractional-order delayed predator–prey systems with Holling type-II functional response. _Nonlinear Dynamics_. 80, 1-2 (2015), 777–789.
  * [3] J. T. Machado, V. Kiryakova, F. Mainardi, Recent history of fractional calculus. _Commun. Nonlinear Sci. Numer. Simulat._ 16, (2011), 1140–1153.
  * [4] I. Matychyn, V. Onyshchenko, Time-optimal control of fractional-order linear systems. _Fract. Calc. Appl. Anal._ 18, 3 (2015), 687–696.
  * [5] R. E. Gutierrez, J. M. Rosario, J. T. Machado, Fractional Order Calculus: Basic Concepts and Engineering Applications. _Math. Prob. Eng._ (2010); DOI: 10.1155/2010/375858.
  * [6] S. Abbas, M. Benchohra, G. M. N’Guerekata, _Topics in fractional differential equations_. Springer Science & Business Media, (2012).
  * [7] S. Das, _Introduction to fractional calculus for scientists and engineers_. Springer, (2011).
  * [8] M. L. Richard, _Fractional calculus in bioengineering_. Redding: Begell House, (2006).
  * [9] Y. Li and Q. Zhang, Blow-up and global existence of solutions for a time fractional diffusion equation. _Fract. Calc. Appl. Anal._ 21, 6 (2018), 1619–1640.
  * [10] J. Munkhammar, Chaos in a fractional order logistic map, _Fract. Calc. Appl. Anal._ 16, 3 (2013), 511–519.
  * [11] M. C. Tripathy, D. Mondal, K. Biswas, S. Sen, Experimental studies on realization of fractional inductors and fractional order bandpass filters. _Int. J. Circuit Theory and Applications._ 43, 9 (2015), 1183–1196.
  * [12] P. J. Torvik, R. L. Bagley, On the appearance of the fractional derivative in the behaviour of real materials. _J. Appl. Mechanics._ 51, 2 (1984), 294–298.
  * [13] J. A. Sabatier, O. P. Agrawal, J. T. Machado, _Advances in fractional calculus_. Dordrecht: Springer, (2007).
  * [14] A. Boukhouima, K. Hattaf, N. Yousfi, Dynamics of a fractional order HIV infection model with specific functional response and cure rate. _Int. J. Diff. Equations._ 2017, (2017).
  * [15] E. Ahmed, A. M. A. El-Sayed, H. A. A. El-Saka, Equilibrium points, stability and numerical solutions of fractional-order predator-prey and rabies models. _J. Math. Anal. Appl._ 325, (2007), 542–553.
  * [16] S. Ranaa, S. Bhattacharyaa, J. Pal, G. M. N’Guerekata, J. Chattopadhyay, Paradox of enrichment: A fractional differential approach with memory. _Physica A._ 392, (2013), 3610–3621.
  * [17] Z. Cui, Z. Yang, Homotopy perturbation method applied to the solution of fractional lotka-volterra equations with variable coefficients. _J. Mod. Meth. Numer. Math._ 5, (2014), 1–9.
  * [18] S. Mondal, N. Bairagi, G. M. N’Guerekata, Global stability of a Leslie-Gower-type fractional order tritrophic food chain model. _Fractional Differential Calculus._ 9, 1 (2019), 149–161.
  * [19] S. Mondal, N. Bairagi, A. Lahiri, A fractional calculus approach to Rosenzweig-MacArthur predator-prey model and its solution. _J. Mod. Meth. Numer. Math._ 8, 1-2 (2017), 66–76.
  * [20] J. Huo, H. Zhao, L. Zhu, The effect of vaccines on backward bifurcation in a fractional order HIV model. _Nonlinear Anal. RWA._ 26, (2015), 289–305.
  * [21] H. L. Li, L. Zhang, C. Hu, Y. L. Jiang, Z. Teng, Dynamical analysis of a fractional-order predator-prey model incorporating a prey refuge. _J. Appl. Math. Comput._ (2016); DOI: 10.1007/s12190-016-1017-8.
  * [22] C. Vargas-De-Leon, Volterra-type Lyapunov functions for fractional-order epidemic systems. _Commun. Nonlinear Sci. Numer. Simul._ 24, (2015), 75–85.
  * [23] S. Mondal, A. Lahiri, N. Bairagi, Analysis of a fractional order eco-epidemiological model with prey infection and type 2 functional response. _Math. Meth. Appl. Sci._ 40, 18 (2017), 6776–6789.
  * [24] K. Hattaf, A. A. Lashari, B. El Boukari,N. Yousfi, Effect of discretization on dynamical behavior in an epidemiological model. _Differential Equations and Dynamical Systems._ 23, 4 (2015), 403–413.
  * [25] M. Sekiguchi, E. Ishiwata, Global dynamics of a discretized SIRS epidemic model with time delay. _Journal of Mathematical Analysis and Applications._ 371, 1 (2010), 195–202.
  * [26] Q. Chen, Z. Teng, L. Wang, H. Jiang, The existence of codimension-two bifurcation in a discrete SIS epidemic model with standard incidence. _Nonlinear Dyn._ 71, (2013), 55–73.
  * [27] Z. Hu, Z. Teng, L. Wang, H. Jiang, Stability analysis in a class of discrete SIRS epidemic models. _Nonlinear Anal.: Real World Appl._ 13, (2012), 2017–2033.
  * [28] M. Biswas, N. Bairagi, Discretization of an eco-epidemiological model and its dynamic consistency. _Journal of Difference Equations and Applications._ 23, 5 (2017), 860–877.
  * [29] S. Ghorai, P. Chakraborty, S. Poria and N. Bairagi, Dispersal-induced pattern forming instabilities in host-parasitoid metapopulations. _Nonlinearity._ (2020), DOI.org/10.1007/s11071-020-05505-w.
  * [30] F. M. Atici, P. W. Eloe, A transform method in discrete fractional calculus. _Int. J. Difference Equations._ 2, 2 (2007), 165–176.
  * [31] C. Fulai, L. Xiannan, Z. Yong, Existence results for nonlinear fractional difference equation. _Advances in Difference Equations._ 2011, (2011); DOI: 10.1155/2011/713201.
  * [32] G. C. Wu, D. Baleanu, Chaos synchronization of the discrete fractional logistic map. _Signal Processing._ 102, (2014), 96–99.
  * [33] A. A. Elsadany, A. E. Matouk, Dynamical behaviors of fractional-order Lotka Volterra predator-prey model and its discretization. _J. Appl. Math. Comput._ 49, (2015), 269–283.
  * [34] J. Alzabut , T. Abdeljawad, D. Baleanu, Nonlinear delay fractional difference equations with applications on discrete fractional Lotka-Volterra competition model. _J. Comput. Anal. Appl._ 25, 5 (2018), 889–898.
  * [35] A.G. M. Selvam, R. Janagaraj, Numerical Analysis of a Fractional Order Discrete Prey – Predator System with Functional Response. _In. J. Engineering & Technology._ 7, 4.10 (2018), 681–684.
  * [36] Y Shi, Q. Ma and X. Ding, 2018. Dynamical behaviors in a discrete fractional-order predator-prey system. _Filomat._ 32, 17 (2018), 5857–5874.
  * [37] S. M. Salman, On a Discretized Fractional-Order SIR Model for Influenza. _Progress in Fractional Differentiation and Application._ 3, 2 (2017), 163–173.
  * [38] A. M. Abdelaziz, A. I. Ismail, F. A. Abdullah and M. H. Mohd, 2018. Bifurcations and chaos in a discrete SI epidemic model with fractional order. _Advances in Difference Equations._ 1 (2018), 1–19.
  * [39] A. Khan, J. F. Gomez-Aguilar, T. S. Khan, H. Khan, Stability analysis and numerical solutions of fractional order HIV/AIDS model. _Chaos, Solitons and Fractals._ 122, (2019), 119–128.

  * [40] A. Khan, T. Abdeljawad, J. F. Gomez-Aguilar, H. Khan, Dynamical study of fractional order mutualism parasitism food web module. _Chaos, Solitons and Fractals._ 134, (2020), 109–685.

  * [41] S. Mondal, M. Biswas, N. Bairagi, Local and global dynamics of a fractional-order predator-prey system with habitat complexity and the corresponding discretized fractional-order system. _J. Applied Mathematics and Computing._ (2020); DOI: https://doi.org/10.1007/s12190-020-01319-6.

  * [42] A. Khan, J. F. Gomez-Aguilar, T. Abdeljawad, H. Khan, Stability and numerical simulation of a fractional order plant-nectar-pollinator model. _Alexandria Engineering Journal._ 59, 1 (2020), 49–59.

  * [43] R. P. Agarwal, A. M. A. El-Sayed, S. M. Salman, Fractional-order Chua’s system: discretization, bifurcation and chaos. _Advances in Difference Equations._ 1, (2013), 320.

  * [44] A. M. A. El-Sayed, S. M. Salman, On a discretization process of fractional-order Riccati differential equation. _J. Fract. Calc. Appl._ 4, 2 (2013), 251–259.

  * [45] A. G. M. Selvam, D. A. Vianny, Behavior of a Discrete Fractional Order SIR Epidemic model. _International Journal of Engineering and Technology._ 7, 4.10 (2018), 675–680.
  * [46] Z. F. El Raheem, S. M. Salman, On a discretization process of fractional-order logistic differential equation. _Journal of the Egyptian Mathematical Society._ 22, 3 (2014), 407–412.
  * [47] J. Chattopadhyay, N. Bairagi, Pelicans at risk in Salton sea–an eco-epidemiological model. _Ecological Modelling._ 136, (2001), 103–112.
  * [48] S. N. Elaydi, _An Introduction to Difference Equations_. Springer, New York, NY, USA, (1999).
  * [49] N. Bairagi, R. R. Sarkar, J. Chattopadhyay, Impacts of incubation delay on the dynamics of an eco-epidemiological system–a theoretical study. _Bulletin of mathematical biology._ 70, 7 (2008), 2017–2038.
  * [50] R. K. Upadhyay, N. Bairagi, K. Kundu, J. Chattopadhyay, Chaos in eco-epidemiological problem of the Salton Sea and its possible control. _Applied Mathematics and Computation._ 196, 1 (2008), 392–401.


[b]A.Yu. Kotov

# Thermal QCD phase transition
with dynamical chiral fermions

Z. Fodor  K. Szabo

###### Abstract

We discuss properties of Quantum Chromodynamics at finite temperature obtained
by means of lattice simulations with overlap fermions. This fermion
discretization preserves chiral symmetry even at finite lattice spacing. We
present details of the lattice formulation, first results for the chiral
observables and discuss the behaviour of the system near the chiral thermal
phase transition.

##  1 Introduction

One of the most interesting results obtained by means of Lattice Quantum
Chromodynamics is the nature of the chiral phase transition at finite
temperature, which turned out to be a crossover, rather than a genuine phase
transition [1]. The major shortcoming of the [1] and subsequent studies of the
thermal QCD is that they are mainly performed with staggered fermions, which
possess only part of the chiral symmetry at finite lattice spacing, while the
full symmetry group is expected to be recovered only after a not always simple
continuum extrapolation. In this proceeding we discuss our current study of
the chiral phase transition in QCD with another fermion discretization, namely
overlap fermions [2], which have the exact chiral symmetry already at finite
lattice spacing. We present and discuss the first results for the chiral
observables and the spectrum of the Dirac operator.

##  2 Lattice details

Our study is based on the same lattice setup, which was used in the previous
studies with heavier pion masses [3, 4]. Now we extent the simulation details
to $N_{f}=2+1$ light quark flavours and physical pion mass. The overlap
implementation of odd number of flavours, which is used for the strange quark
sector, was introduced in [5]. We summarize the details of the action below:

  * •

tree-level Symanzik improved gauge action

  * •

two light quark flavours with the action given by overlap Dirac operator:

$aD_{\mathrm{ov}}=\frac{1}{2}\left(1+\gamma_{5}\mathrm{sign}(\gamma_{5}D_{W}(-m_{W}))\right),$
(1)

for the kernel we used Wilson Dirac operator $D_{W}(-m_{W})$ with negative
mass $-m_{W}=-1.3$

  * •

one strange quark flavour of the overlap Dirac operator with the same action
(1)

  * •

two fermion flavours with Wilson action, given by $D_{W}(-m_{W})$. They are
irrelevant in the continuum limit

  * •

two boson fields (also irrelevant in the continuum limit) with the mass
$m_{B}=0.54$ and the action:

$\phi^{{\dagger}}\left[D_{W}(-m_{W})+im_{B}\gamma_{5}\tau_{3}\right]\phi$ (2)

$N_{t}$ | $N_{s}$ | Aspect ratio $N_{s}/N_{t}$
---|---|---
8 | 16, 24, 32 | 2, 3, 4
10 | 20, 30, 40 | 2, 3, 4
12 | 24, 36, 48 | 2, 3, 4
Table 1: Summary of lattice sizes used in the current study

The main goal for the introduction of the auxiliary bosonic and fermion fields
is to suppress the tunnellings between different topological sectors [6]. In
this study we restrict simulations to the zero topological sector $Q=0$.
Simulations in the fixed topological sector lead to finite volume effects [6],
which disappear in the limit $V\to\infty$.

Generation of the configurations was done with the help of the standard Hybrid
Monte-Carlo algorithm. For the application of the overlap Dirac operator one
needs to calculate the sign function of a hermitian Wilson-Dirac operator
$X=\gamma_{5}D_{W}(-m_{W})$: $\mathrm{sign}X=\frac{X}{\sqrt{XX^{{\dagger}}}}$.
In order to calculate it we first computed 32-128 low-lying modes of operator
$X$, on which the sign function can be calculated explicitly, and for the rest
of the spectrum we used Chebyshev polynomial approximation of the sign
function. For the calculation of eigenvalues we used implicitly restarted
Lanczos algorithm and the inversion of the overlap Dirac operator was done
using flexible GMRES method with Wilson-Dirac operator as preconditioner.

We performed numerical calculations along the line of constant physics (LCP),
which was determined from calculations with larger pion masses in [5]. In Tab.
1 we summarize the ensembles which we use in our study. For each value of
$N_{t}$ and $N_{s}$ we performed simulations for eight temperatures:
$T=140,145,150,155,160,165,170,175$ MeV.

Figure 1: Chiral condensate as a function of temperature for all lattice sizes
used in the current study. Solid lines correspond to the continuum
extrapolation at fixed aspect ratio $N_{s}/N_{t}$.

##  3 Chiral phase transition

For the study of the chiral phase transition we measured the chiral condensate
and the chiral susceptibility. Lattice data at finite lattice spacing have to
be renormalized. In this study we employed the renormalization procedure,
based on the strange quark condensate. The renormalized observables are
defined then:

$\begin{split}\langle\bar{\psi}\psi\rangle_{r}&=2(m_{s}\langle\bar{\psi}\psi\rangle_{l}-m_{l}\langle\bar{\psi}\psi\rangle_{s})\\\
\chi_{r}&=m_{s}\partial_{m_{l}}M\end{split}$ (3)

This renormalization procedure does not eliminate all additive divergences,
logarithmic terms $\sim m^{3}\log ma$ remain in Eq. (3), however their
numerical values is expected to be small due to suppression by the quark mass
[7].

Figure 2: Chiral susceptibility as a function of temperature for all lattice
sizes used in the current study..

In Fig. 1 we present the chiral condensate for all ensembles used in the
current study, shown by dashed lines. One can clearly see that results
significantly depend both on the lattice spacing and on the spatial volume. In
addition, in Fig. 1 we also show the continuum extrapolation of the data for
each volume, obtained with a simple $\sim a^{2}$ behaviour. The results also
suggest large volume dependence, for aspect ratio $N_{s}/N_{t}=2$ there is no
clear inflection point in the studied temperature range $T\in[140,175]$ MeV,
while larger aspect ratios $N_{s}/N_{t}=3,4$ exhibit an expected threshold
behaviour with inflection point $T\sim 160$ MeV.

In Fig. 2 we show the chiral susceptibility for all ensembles. The behaviour
of the chiral susceptibility is consistent with the behaviour of the chiral
condensate: for small aspect ratio $N_{s}/N_{t}=2$ there is no clear peak and
consequently no indication of the phase transition $T\in[140,175]$ MeV, while
large aspect ratios $N_{s}/N_{t}=3,4$ exhibit a peak-like behaviour with
pseudo-critical temperature $T\sim 160$ MeV. We would like to note that recent
studies with the staggered fermions at finite volume [8] suggest that for
aspect ratio $N_{s}/N_{t}=2$ the pseudo-critical temperature is smaller $T\sim
140$ MeV than the pseudo-critical temperatue in the infinite volume $T\approx
155$ MeV, which is in agreement with our results obtained with the overlap
fermions. It is worth noting that our studies are done in the fixed
topological sector $Q=0$, which is another source of finite volume effects and
can also affect the results in the fixed volume. We plan to study how
simulations in the fixed topological sector $Q=0$ affect our data in the
future.

Figure 3: Continuum extrapolated spectrum of the Dirac operator
$\rho(\lambda/m_{q})$ versus $\lambda/m_{q}$ for aspect ratios
$N_{s}/N_{t}=2$, $3$, $4$.

##  4 Spectrum of the Dirac operator

Apart from the chiral condensate, we also studied the spectrum of the overlap
Dirac operator, described by the spectral function $\rho(\lambda)$. It is
related to the chiral condensate via Banks-Casher relation [9]:

$\langle\bar{\psi}\psi\rangle_{r}=2\int_{0}^{\infty}d\lambda\rho(\lambda)\frac{m}{m^{2}+\lambda^{2}}\xrightarrow[m\to
0]{}\pi\rho(\lambda).$ (4)

The behaviour of $\rho(\lambda)$ has been subject of several studies [10, 11,
12], which suggest rather nontrivial behaviour near origin with a formation of
the peak, possibly indicating some new phase in QCD. In Fig. 3 we present the
results for the continuum extrapolation of the spectrum for all values of
$N_{t}$ and aspect ratios used in our spectrum. In order to avoid issues with
the renormalization we used the ratio $\lambda/m_{q}$, in which
renormalization constants are cancelled [13].

We observe that the results for the spectral density $\rho$ almost do not
depend on the aspect ratio, except for the behaviour at small $\lambda$, where
for the smallest studied aspect ratio the spectral density
$\rho(\lambda/m_{q})$ is suppressed as compared to larger aspect ratios. It
can be attributed to the fact that low-lying large distance modes simply do
not fit into small lattice volumes. For this reason we also do not observe
standard signatures of the transition-like behaviour at small $N_{s}/N_{t}=2$.
However, data at larger $N_{s}/N_{t}=3,4$ show a nonzero intercept as
$\lambda\to 0$ which disappears if one increases temperature at $T\approx
160-165$ MeV. This behaviour is in a full agreement with the results extracted
from the chiral condensate and the chiral susceptibility. Our data also do not
suggest the existence of the large peak near origin for the temperatures and
volumes under study.

##  5 Conclusions

We presented the first results of our study of the chiral QCD phase transition
with dynamical overlap fermions. This fermion discretization possesses chiral
symmetry at finite lattice spacing. We restricted ourselves to the simulations
in the fixed toplogical sector $Q=0$. We showed the dependence of the chiral
condensate and the chiral susceptibility on the temperature in the vicinity of
the chiral phase transition $T\in[140,175]$ MeV for three values of
$N_{t}=8,10$ and $12$ and three spatial volumes corresponding to aspect ratios
$N_{s}/N_{t}=2,3,4$. Our results for the smallest aspect ratio $N_{s}/N_{t}=2$
most likely are subject to large finite volume effects and do not suggest the
existence of the phase transition or crossover in the temperature range
$T\in[140,175]$ MeV, while for two larger aspect ratios $N_{s}/N_{t}=3,4$ the
pseudo-critical temperature $T\approx 160$ MeV is in the agreement with
previous staggered studies. We also presented the results for the spectrum of
the Dirac operator, which is consistent with the same picture and the same
estimation of the pseudo-critical temperature.

##  6 Acknowledgements

The authors gratefully acknowledge computing time on the supercomputer
JURECA[14] at Forschungszentrum Jülich under grant no. qcdoverlap.

## References

  * [1] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo.  The Order of the quantum chromodynamics transition predicted by the standard model of particle physics.  Nature, 443:675–678, 2006.
  * [2] Herbert Neuberger.  Exactly massless quarks on the lattice.  Phys. Lett. B, 417:141–144, 1998.
  * [3] Szabolcs Borsanyi, Ydalia Delgado, Stephan Durr, Zoltan Fodor, Sandor D. Katz, Stefan Krieg, Thomas Lippert, Daniel Nogradi, and Kalman K. Szabo.  QCD thermodynamics with dynamical overlap fermions.  Phys. Lett. B, 713:342–346, 2012.
  * [4] Sz. Borsanyi, Z. Fodor, S. D. Katz, Stefan F. Krieg, T. Lippert, D. Nogradi, F. Pittler, K. K. Szabo, and B. C. Toth.  QCD thermodynamics with continuum extrapolated dynamical overlap fermions.  10 2015.
  * [5] Sz. Borsanyi et al.  Calculation of the axion mass based on high-temperature lattice quantum chromodynamics.  Nature, 539(7627):69–71, 2016.
  * [6] Hidenori Fukaya, Shoji Hashimoto, Ken-Ichi Ishikawa, Takashi Kaneko, Hideo Matsufuru, Tetsuya Onogi, and Norikazu Yamada.  Lattice gauge action suppressing near-zero modes of H(W).  Phys. Rev. D, 74:094505, 2006.
  * [7] H. T. Ding et al.  Chiral Phase Transition Temperature in ( 2+1 )-Flavor QCD.  Phys. Rev. Lett., 123(6):062002, 2019.
  * [8] Szabolcs Borsányi, Zoltán Fodor, Jana N. Guenther, Ruben Kara, Paolo Parotto, Attila Pásztor, and Chik Him Wong.  Finite volume effects near the chiral crossover.  1 2024.
  * [9] Tom Banks and A. Casher.  Chiral Symmetry Breaking in Confining Theories.  Nucl. Phys. B, 169:103–125, 1980.
  * [10] Andrei Alexandru and Ivan Horváth.  Possible New Phase of Thermal QCD.  Phys. Rev. D, 100(9):094507, 2019.
  * [11] Robin Kehr, Dominik Smith, and Lorenz von Smekal.  QCD Anderson transition with overlap valence quarks on a twisted-mass sea.  4 2023.
  * [12] Tamas G. Kovacs.  The fate of chiral symmetries in the quark-gluon plasma.  11 2023.
  * [13] Leonardo Giusti and Martin Luscher.  Chiral symmetry breaking and the Banks-Casher relation in lattice QCD with Wilson quarks.  JHEP, 03:013, 2009.
  * [14] Jülich Supercomputing Centre.  JURECA: Data Centric and Booster Modules implementing the Modular Supercomputing Architecture at Jülich Supercomputing Centre.  Journal of large-scale research facilities, 7(A182), 2021.