# Pattern alternations induced by nonlocal interactions

Swadesh Pal<EMAIL_ADDRESS>Roderick Melnik<EMAIL_ADDRESS>Malay Banerjee
<EMAIL_ADDRESS>MS2Discovery Interdisciplinary Research Institute, Wilfrid
Laurier University, Waterloo, Canada  BCAM - Basque Center for Applied
Mathematics, E-48009, Bilbao, Spain  Department of Mathematics and Statistics,
IIT Kanpur, Kanpur, India

###### Abstract

Pattern formation is a visual understanding of the dynamics of complex
systems. Patterns arise in many ways, such as the segmentation of animals,
bacterial colonies during growth, vegetation, chemical reactions, etc. In most
cases, the long-range diffusion occurs, and the usual reaction-diffusion (RD)
model can not capture such phenomena. The nonlocal RD model, on the other
hand, can fill the gap. Analytical derivation of the amplitude equations (AE)
for an RD system is a valuable tool to predict the pattern selections, in
particular, the stationary Turing patterns when they occur. In this paper, we
analyze the conditions for the Turing bifurcation for the nonlocal model and
also derive the AE for the nonlocal RD model near the Turing bifurcation
threshold to describe the reason behind the pattern selections. This
derivation of the AE is not only limited to the nonlocal prey-predator model,
as shown in our representative example but also can be applied to other
nonlocal models near the Turing bifurcation threshold. The analytical
prediction agrees with numerical simulation near the Turing bifurcation
threshold. Moreover, the analytical and numerical results fit each other well
even more remote from the Turing bifurcation threshold for the small values of
the nonlocal parameter but not for the higher values.

###### keywords:

Nonlocal model , Kernel function , Turing bifurcation , Spatial-Hopf
bifurcation , Amplitude equations , Weakly nonlinear analysis

##  1 Introduction

Ordinary differential equation models of interacting populations assume their
homogeneous distribution over their habitat [1]. However, in reality, the
distribution of individuals of different species is heterogeneous and guided
by the uneven distribution of favourable resources. A reaction-diffusion model
of an interacting population can capture the heterogeneous distribution of the
individuals of constituent species and their random movement within their
habitat. The main focus of the study with reaction-diffusion models is to
understand the species’ self-organized distribution (patterns) due to the
intra- and inter-specific interaction for favourable resources and survival.
Non-homogeneous stationary distributions are classified as spatial self-
organization patterns, and time-dependent heterogeneous distributions are
known as dynamic patterns.

In the ecological systems, populations interact with the other individuals at
their spatial location and in nearby locations. The usual reaction-diffusion
system can not capture such type of phenomena; however, a reaction-diffusion
system with nonlocal interaction can. Different type of nonlocal models have
been studied for prey-predator interactions in different circumstances, e.g.,
nonlocal dispersal [2], nonlocal consumption of resources [3], nonlocal
intraspecific competition [4, 5, 6, 7], etc. The use of reaction-diffusion
equations with nonlocal interaction terms is a newly emerged area of research,
whereas this modelling approach is not limited by the models of population
interactions. Rather this approach is well accepted for the models in
biomedical applications, the study of various diseases, nanotechnology and
neuroscience [8, 9, 10, 11, 12, 13, 14].

One of our motivations in this paper is derived from the fact that resource-
consumer interactions are one of the main building blocks of several food webs
and food chains. The dynamics exhibited by the resource-consumer models
capture the complex interaction between several trophic levels in natural
ecosystems. These models include the dynamical relationships between
autotroph-herbivore, prey-predator, host-parasites, etc. According to Abrams
[15], ‘eating is a necessity for all heterotrophic organisms when the foods
are themselves living organisms, the interaction between consumer and food is
predation’. A wide variety of mathematical models are available in the
literature for prey-predator interactions, but essentially they can be
classified based on the Lotka-Volterra type model, Gause type model, and
Kolmogorov model [16]. The predator population exerts negative feedback on the
prey growth through predation, and at the same time, the prey population
contributes to the development of predators by supplying energy through food,
which helps in the production of new off-springs [17]. These two mechanisms
are modelled through the functional and numerical responses, respectively. In
literature, we find mainly two types of functional responses, namely prey-
dependent and predator-dependent [18, 19]. For prey-dependent functional
response, the rate and amount of consumption of the prey biomass by their
predators depend solely upon the prey population density. On the other hand,
when the functional response depends upon both the prey and predator
population densities, it is called predator-dependent functional response
without any ambiguity. However, a numerical response is measured through
scaling of the functional response. Various types of biological mechanisms are
responsible for the influence of predator population density to shape the
functional response.

Mutual interference among the predators, difficulty in getting food at a low
prey-to-predator ratio, forming groups to enhance the success in catching and
handling large prey – are some well-known mechanisms responsible for the
inclusion of predator population density in the functional response. Hunting
cooperation describes the cooperative mechanism among the predator individuals
to have success in searching, catching, and handling prey [20]. The
incorporation of hunting cooperation in the modelling approach results in the
inclusion of predator population density in functional and numerical response
terms, and hence the functional response is predator-dependent. As a result,
the dynamics described by the model become a little bit more complex compared
to its counterpart with the prey-dependent functional response. For prey-
predator models with specialist predators, the hunting cooperation among the
predators can enhance the success in catching prey; however, it exploits the
resource excessively, resulting in predator extinction due to scarcity of food
at high population density. This negative feedback in the growth of specialist
predators is known as the component Allee effect in predators.

Two types of spatial patterns occur in a reaction-diffusion model: stationary
and non-stationary. There are several mechanisms for generating various non-
stationary or dynamic patterns. However, the existence of a stationary pattern
is related to Turing instability. This condition gives a sufficient analytical
prerequisite for the existence of the stationary pattern formation [21, 22].
The challenging issue is to obtain the analytical prediction about the types
of the stationary patterns for the parameter values away from the Turing
bifurcation thresholds. On the other hand, near the Turing bifurcation
threshold, a finite number of Fourier modes can capture such type of pattern
forming scenarios for periodic boundary conditions over a spatial domain [23,
24]. In this case, the solution of the partial differential equation (PDE)
model can be approximated by a solution of a system of ordinary differential
equations, called amplitude equations. In 1965, Eckhaus used it for the first
time to reduce the Navier-Stokes equations of fluid mechanics to a system of
ordinary differential equations for the amplitudes [25]. This reduction
technique has been used in the theory of nonlinear wave interactions to
predict the patterns of reaction-diffusion equations. Later on, Segal and
Levin derived the amplitude equations for prey-predator systems by the
technique known as the weakly nonlinear analysis [26]. This result is a normal
form of the spatio-temporal model near the Turing bifurcation threshold.

In this work, we use a basic form of hunting cooperation in prey-predator
interaction for the temporal dynamics [20]. We study the existence of the non-
trivial equilibrium points of the model and also their stability behaviours
through different temporal bifurcations. The most popular weakly nonlinear
analysis helps predict the underlying stationary Turing patterns analytically
for the local and nonlocal reaction-diffusion models. In this work, we analyze
the nonlocal interaction in intraspecific prey competition. For the kernel
function, we restrict our attention to “thin-tailed” (Gaussian) function [2].
To the best of our knowledge, the amplitude equation for the spatio-temporal
models with nonlocal interactions remains unexplored and it is the primary
focus of this work. In general, the weakly nonlinear analysis is a multiscale
analysis with respect to a bifurcation parameter and the time derivative. In
this method, the bifurcation parameter and the time derivative are expanded in
terms of a small parameter [27, 28, 29, 30].

##  2 Mathematical Model

For the time $t>0$, suppose $u(t)$ and $v(t)$ are the prey and predator
populations, respectively. Following [20], we consider the simplest form of a
prey-predator interaction with hunting cooperation among the predators as:

$\displaystyle\frac{du}{dt}$ $\displaystyle=\eta
u\left(1-\frac{u}{\kappa}\right)-(1+\alpha v)uv\equiv N_{1}(u,v),$ (1a)
$\displaystyle\frac{dv}{dt}$ $\displaystyle=(1+\alpha v)uv-v\equiv
N_{2}(u,v),$ (1b)

with non-negative initial conditions. Here, $\eta$ is the per capita intrinsic
growth rate and $\kappa$ is the environmental carrying capacity. The parameter
$\alpha$ represents the predator hunting cooperation, and $a=0$ corresponds to
the prey-predator model without hunting cooperation [20]. All the parameters
involved in the model are assumed to be positive.

We extend the temporal model (1) into the spatio-temporal model to account for
the random movements of the species in a two-dimensional habitat. In the
mathematical model, the diffusion term captures such types of movements
phenomena, and it has been widely accepted in the spatio-temporal prey-
predator model. Along with this, researchers have been studied different type
of nonlocal interactions in the prey-predator model [2, 3, 4, 5, 6]. Following
[4], we consider the nonlocal interaction in the intraspecific competition of
the prey population, and in this case, the corresponding integro-differential
reaction-diffusion model is given by

$\displaystyle\frac{\partial u}{\partial t}$ $\displaystyle=\Delta u+\eta
u\left(1-\frac{\psi_{\sigma}\ast u}{\kappa}\right)-(1+\alpha v)uv,$ (2a)
$\displaystyle\frac{\partial v}{\partial t}$ $\displaystyle=d\Delta
v+(1+\alpha v)uv-v,$ (2b)

with non-negative initial conditions and periodic boundary conditions. We have
chosen the non-dimensional diffusion coefficients for the prey and predator
species to be $1$ and $d$, respectively. The convolution term
$\psi_{\sigma}\ast u$ is defined as

$(\psi_{\sigma}\ast
u)(x,y,t)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\psi_{\sigma}(x-w,y-z)u(w,z,t)dwdz,$

where $\psi_{\sigma}$ is the kernel function. The parameter $\sigma$ has a
prominent role in the nonlocal interactions, it captures the effective area of
nonlocal interactions for the nonlocal model (2). We assume that the kernel
function is non-negative, even, normalized and exponentially bounded in
$\mathbb{R}^{2}$. Different type of kernel functions are available in the
literature [2, 4, 5], but, in this work, we restrict our attention to the
Gaussian kernel because it is the most widely used kernel functions in
ecological models [2, 31]. It is given by

${}\psi_{\sigma}(x,y)=\frac{1}{2\pi\sigma^{2}}e^{-\frac{x^{2}+y^{2}}{2\sigma^{2}}}.$
(3)

This specific choice of the kernel function satisfies all the above-mentioned
assumptions. Now, we first analyze the temporal dynamics of the model, and
then we move toward the local and nonlocal models.

##  3 Analysis of the temporal model

In this section, we study the non-spatial model (1). For analyzing the
temporal model, we generally study the equilibrium points and their
stabilities, which forwards to the model’s bifurcation analysis. Now, the
model (1) admits a trivial equilibrium point $E_{0}=(0,0)$ and an axial
equilibrium point $E_{1}=(\kappa,0)$. Furthermore, the coexisting equilibrium
point (or points) of the system (1) is (are) the point (or points) of the
intersection of the non-trivial nullclines

$\displaystyle u=\frac{\kappa}{\eta}(\eta-(1+\alpha v)v)\equiv
n_{1}(v)\leavevmode\nobreak\ \mbox{and}\leavevmode\nobreak\
u=\frac{1}{1+\alpha v}\equiv n_{2}(v),$ (4)

inside the first quadrant of the $uv$-plane. Suppose, $E_{*}=(u_{*},v_{*})$
denotes the components of the coexisting equilibrium point, then $u_{*}(>0)$
satisfies the polynomial equation

$\displaystyle\phi(u)\equiv\alpha\eta u^{3}-\alpha\kappa\eta u^{2}-\kappa
u+\kappa=0$ (5)

and $v_{*}$ satisfies

${}v_{*}=\frac{1}{\alpha}\left(\frac{1}{u_{*}}-1\right).$ (6)

Note that $u_{*}$ has to be less than $1$ to satisfy the positivity of $v_{*}$
(see eq. (6)). The non-trivial prey nullcline $u=n_{1}(v)$, lying in the first
quadrant of the $uv$-plane, is a monotone decreasing function in $v$, concave
with respect to $u$-axis, and passes through the point $(\kappa,0)$ in the
$uv$-plane. On the other hand, the non-trivial predator nullcline
$u=n_{2}(v)$, lying in the first quadrant, is a segment of the hyperbola,
passes through the point $(1,0)$ in the $uv$-plane. The geometry of these two
nullclines ensure the existence of the unique coexisting equilibrium point
$E_{*}$ for $\kappa>1$, e.g., see Fig. 1(1(a)). Furthermore, two coexisting
equilibrium points exist for $\kappa<1$, and in this case, two nullclines
intersect each other at two different points in the first quadrant of the
$uv$-plane, e.g., see Fig. 1(1(b)). These points of intersections can not be
determined explicitly because of the involvement of the cubic polynomial
equation (5), and they also depend on the magnitudes of other parameters.

As we have mentioned earlier, all the parameters involved in the cubic
equation (5) are positive. From the Descartes rule of sign, the equation (5)
possesses at most two positive roots, and the feasibility of coexisting
equilibrium point demands $u_{*}<1$. Depending on the parameter values, the
coexisting equilibrium point can bifurcate from the axial equilibrium point or
can be generated through saddle-node bifurcation. We first present the
analytical conditions for the existence and the stability behaviours of these
two bifurcations and verify them numerically later on. In most cases, the
analytical conditions are implicit in nature.

(a)

(b)

Figure 1: (Color online) Two different scenarios of the non-trivial nullclines
for the system (1). The dotted magenta color curve represents the prey
nullcline, and the solid blue color curve represents the predator nullcline.

Before moving towards the other equilibrium points, we first summarize the
stability behaviours of the trivial equilibrium point. With the help of linear
stability analysis, we can obtain the local asymptotic stability of all the
equilibrium points of the system (1). As we can see, the trivial equilibrium
point $E_{0}$ is always unstable, and in particular, it is a saddle point
having unstable and stable sub-spaces along the $u$-axis and $v$-axis,
respectively. Next, we summarize the stability behaviour of the predator-free
equilibrium point $E_{1}$ in the following theorem.

######  Theorem 1.

The system (1) undergoes a transcritical bifurcation at $\kappa_{TC}=1$. The
equilibrium point $E_{1}$ is stable for $\kappa<1$ and is unstable for
$\kappa>1$.

###### Proof.

The Jacobian matrix for the system (1), evaluated at $E_{1}=(\kappa,0)$, is
given by

$\displaystyle\textbf{J}_{1}=\left[\begin{array}[]{cc}-\eta&-\kappa\\\
0&-1+\kappa\\\ \end{array}\right].$

Therefore, the equilibrium point $E_{1}$ is stable for $\kappa<1$ and unstable
for $\kappa>1$. Furthermore, the matrix $\textbf{J}_{1}$ has a simple zero-
eigenvalue for $\kappa=1$. Now, at $\kappa=1$, the eigenvectors of
$\textbf{J}_{1}$ and $\textbf{J}_{1}^{T}$ associated with the simple zero-
eigenvalues are $z=(1,-\eta)^{T}$ and $w=(0,1)^{T}$, respectively. Finally, we
obtain the transversality conditions for the transcritical bifurcation [32] at
$\kappa=\kappa_{TC}$, and they are given by

$w^{T}\mathbf{F}_{\kappa}(E_{1},\kappa_{TC})=0,$
$w^{T}\left[D\mathbf{F}_{\kappa}(E_{1},\kappa_{TC})z\right]=0,$
$w^{T}\left[D^{2}\mathbf{F}(E_{1},\kappa_{TC})(z,z)\right]=2\eta(\alpha\eta-1)\neq
0,$

where $\mathbf{F}=(N_{1},N_{2})^{T}$. In particular, the transcritical
bifurcation is degenerate as
$w^{T}\left[D\mathbf{F}_{\kappa}(E_{1},\kappa_{TC})z\right]=0$. ∎

######  Theorem 2.

The system (1) undergoes a saddle-node bifurcation at $\eta_{SN}$ for
$u_{*}<1$, and the bifurcation threshold is defined implicitly as

$\eta_{SN}=\frac{\kappa(2-u_{*})}{\alpha u_{*}^{3}}.$

###### Proof.

As we have discussed earlier, two interior equilibrium points can be generated
through a saddle-node bifurcation when two non-trivial nullclines
$f(u,v)=\eta(1-u/\kappa)-(1+\alpha v)v$ and $g(u,v)=(1+\alpha v)u-1$ touch
(excluding the intersecting case) each other at $E_{*}=(u_{*},v_{*})$ in the
first quadrant. In this case, at $E_{*}$, both the nullclines share a common
tangent. Therefore, at $E_{*}$, we must have

${}\left.-\frac{f_{u}}{f_{v}}\right|_{E_{*}}=\left.-\frac{g_{u}}{g_{v}}\right|_{E_{*}}.$
(8)

Now, the Jacobian matrix for the system (1) at $E_{*}$ is given by

$\displaystyle\textbf{J}_{*}=\left[\begin{array}[]{cc}uf_{u}&uf_{v}\\\
vg_{u}&vg_{v}\\\ \end{array}\right]_{E_{*}}.$

By using (8), we obtain $\det(\textbf{J}_{*})=0$ and
$\textrm{tr}(J_{*})\,=\,u_{*}(\alpha v_{*}-\eta/\kappa)\neq 0$. Therefore, the
matrix $\textbf{J}_{*}$ has a zero eigenvalue with multiplicity one. If we
solve $\det(\textbf{J}_{*})=0$ implicitly for $\eta$, we can find the
threshold for $\eta_{SN}$ by which the saddle-node bifurcation exists. To
verify the transversality conditions for the saddle-node bifurcation, we
consider $\eta$ as the bifurcation parameter denoted by $\eta_{SN}$. The
saddle-node bifurcation threshold $\eta_{SN}$ is a positive root of the
implicit equation $\alpha\eta u_{*}=\kappa(1+\alpha v_{*})(1+2\alpha v_{*})$.
Hence, at $\eta=\eta_{SN}$, the matrix $\textbf{J}_{*}$ has a simple zero
eigenvalue, while its eigenvectors, associated with zero eigenvalue of the
matrices $\textbf{J}_{*}$ and $\textbf{J}_{*}^{T}$, are

$\displaystyle z=\left[\begin{array}[]{c}1\\\ -\frac{1+\alpha v_{*}}{\alpha
u_{*}}\\\ \end{array}\right]\,\,\textrm{and}\,\,w=\left[\begin{array}[]{c}1\\\
\frac{1+2\alpha v_{*}}{\alpha v_{*}}\\\ \end{array}\right],$

respectively. The transversality conditions for the saddle-node bifurcation
are as follows:

$w^{t}\mathbf{F}_{\eta}(E_{*},\eta_{SN})\,=\,u_{*}\left(1-\frac{u_{*}}{\kappa}\right)\,\neq\,0,$
$w^{t}\left[D^{2}\mathbf{F}(E_{*},\eta_{SN})(z,z)\right]\,=\,-\frac{2\eta}{\kappa}-\frac{2}{\alpha
u_{*}^{2}}(1+\alpha v_{*})\,\neq\,0.$

∎

The considered model admits either one or two interior equilibrium points, but
none of the components can be determined explicitly in either of these cases.
This prevents us from finding the exact condition for local asymptotic
stability of the coexisting equilibrium point. However, we summarize the local
asymptotic stability condition for a typical coexisting equilibrium point
$E_{*}=(u_{*},v_{*})$ in an implicit way in the following theorem.

######  Theorem 3.

The coexisting equilibrium point $E_{*}=(u_{*},v_{*})$ of the system (1) is
locally asymptotically stable if $\alpha>\kappa v_{*}/\eta$ and
$\kappa(1+\alpha v_{*})(1+2\alpha v_{*})>\eta\alpha u_{*}$ with $u_{*}<1$.

###### Proof.

For the feasibility of $v_{*}$, we assume $u_{*}<1$. The Jacobian matrix for
(1) at $E_{*}$ is given by

$\textbf{J}_{*}=\begin{bmatrix}-\frac{\eta u_{*}}{\kappa}&-(1+2\alpha
v_{*})u_{*}\\\ (1+\alpha v_{*})v_{*}&\alpha u_{*}v_{*}\end{bmatrix}.$ (11)

Therefore, the equilibrium point $E_{*}$ is locally asymptotically stable if
$\mbox{tr}(\textbf{J}_{*})<0$ and $\det(\textbf{J}_{*})>0$ hold, i.e.,
$\alpha>\kappa v_{*}/\eta$ and $\kappa(1+\alpha v_{*})(1+2\alpha
v_{*})>\eta\alpha u_{*}$. ∎

Depending on the parameter values, sometimes the implicit condition
$\kappa(1+\alpha v_{*})(1+2\alpha v_{*})>\eta\alpha u_{*}$ holds, but the
other condition $\alpha>\kappa v_{*}/\eta$ is violated, and in this case, the
coexisting equilibrium point loses its stability through the Hopf bifurcation.
We summarize all the conditions for the Hopf bifurcation in the following
theorem.

######  Theorem 4.

The coexisting equilibrium point $E_{*}=(u_{*},v_{*})$ with $u_{*}<1$
undergoes a Hopf bifurcation when $\alpha$ crosses the implicit threshold
$\alpha_{H}=\kappa v_{*}/\eta$, while maintaining the inequality
$\kappa(1+\alpha v_{*})(1+2\alpha v_{*})>\eta\alpha u_{*}$.

###### Proof.

The Jacobian matrix for (1) at $E_{*}$ is given in (11) and the trace of the
matrix $\textbf{J}_{*}$ is equal to $(\alpha v_{*}-\eta/\kappa)u_{*}$. The
trace of the Jacobian matrix $\textbf{J}_{*}$ equated to zero gives the Hopf
bifurcation threshold in terms of $\alpha$ implicitly as $\alpha_{H}=\kappa
v_{*}/\eta$. The other condition $\kappa(1+\alpha v_{*})(1+2\alpha
v_{*})>\eta\alpha u_{*}$ ensures that the determinant of the Jacobian matrix
$\textbf{J}_{*}$ is positive. The transversality condition for the Hopf
bifurcation is given by

$\left.\frac{d}{d\alpha}(\textrm{Re}(\lambda))\right|_{E_{*};\alpha_{H}}\neq
0,$ (12)

where $\lambda$ is a root of the Jacobian matrix $\textbf{J}_{*}$. ∎

Until now, we have explained the generation of coexistence equilibrium points
through saddle-node bifurcation and their destabilization through Hopf
bifurcation. Due to the implicit involvement of the non-trivial equilibrium
point with the parameters, we have chosen $\eta$ and $\alpha$ as bifurcation
parameters for finding the bifurcation thresholds. We can also find saddle-
node and Hopf bifurcation curves in the two-dimensional parameter plane by
considering $\eta$ and $\alpha$ as bifurcation parameters. Therefore, it is
expected that the two local bifurcation curves will intersect at a Bogdanov-
Takens bifurcation point of co-dimension two.

######  Theorem 5.

The coexisting equilibrium point $E_{*}$ undergoes a Bogdanov-Takens
bifurcation when the implicit parametric conditions $\alpha_{BT}=\kappa
v_{*}/\eta$ and $\eta_{BT}=\kappa(2-u_{*})/\alpha u_{*}^{3}$ are satisfied.

###### Proof.

At $E_{*}$, the trace and determinant of $\textbf{J}_{*}$ are given by

$\textrm{tr}(\textbf{J}_{*})=-\frac{\eta u_{*}}{\kappa}+\alpha
u_{*}v_{*},\,\,\textrm{det}(\textbf{J}_{*})\,=\,u_{*}v_{*}\left[(1+\alpha
v_{*})(1+2\alpha v_{*})-\frac{\eta\alpha u_{*}}{\kappa}\right].$

The Bogdanov-Takens bifurcation can be obtained by solving
$\textrm{tr}(\textbf{J}_{*})=0$ and $\textrm{det}(\textbf{J}_{*})=0$
simultaneously for $\eta$ and $\alpha$ with the facts $u_{*},v_{*}>0$ and
$1+\alpha v_{*}=1/u_{*}$. Therefore, from these conditions, we find

$\alpha_{BT}=\frac{\kappa
v_{*}}{\eta},\,\,\eta_{BT}\,=\,\frac{\kappa(2-u_{*})}{\alpha u_{*}^{3}}.$

This bifurcation ensures that the matrix $\textbf{J}_{*}$ has zero as an
eigenvalue with multiplicity two. Moreover, satisfying all the parametric
restrictions mentioned above, the Jordan canonical form of $\textbf{J}_{*}$ at
$(\alpha,\,\eta)\,=\,(\alpha_{BT},\,\eta_{BT})$ can be obtained as
$\left[\begin{array}[]{cc}0&1\\\ 0&0\\\ \end{array}\right]$. ∎

##  4 Analysis of the nonlocal spatio-temporal model

In this section, we first find the Turing bifurcation conditions for the local
and nonlocal models and then move towards the weakly nonlinear analysis for
both models. Note that the equilibrium points of the temporal model are the
homogeneous solutions for both local and nonlocal models. This invariance of
homogeneous steady-state is due to the choice of the kernel function and the
boundary conditions.

###  4.1 Linear stability analysis

We first assume that the homogeneous steady-state $E_{*}=(u_{*},v_{*})$ of the
system (1) is locally asymptotically stable, i.e.,
$\mbox{tr}(\mathbf{J}_{*})<0$ and $\mbox{det}(\mathbf{J}_{*})>0$. Now, we
perturb the homogeneous steady-state by
$u=u_{*}+\epsilon\tilde{u}_{1}\exp{(\lambda t+i(k_{x}x+k_{y}y))}$ and
$v=v_{*}+\epsilon\tilde{v}_{1}\exp{(\lambda t+i(k_{x}x+k_{y}y))}$, where
$|\epsilon|\ll 1$. Substitution of $u$ and $v$ into the system (2) and
linearization leads to

${}\begin{bmatrix}a_{11}-(\eta/\kappa)u_{*}\widehat{\psi}_{\sigma}(k_{x},k_{y})-(k_{x}^{2}+k_{y}^{2})-\lambda&a_{12}\\\
a_{21}&a_{22}-d(k_{x}^{2}+k_{y}^{2})-\lambda\end{bmatrix}\begin{bmatrix}\tilde{u}_{1}\\\
\tilde{v}_{1}\end{bmatrix}\equiv\mathbf{M}\begin{bmatrix}\tilde{u}_{1}\\\
\tilde{v}_{1}\end{bmatrix}=\begin{bmatrix}0\\\ 0\end{bmatrix},$ (13)

where $a_{11}=0$, $a_{12}=-(1+2\alpha v_{*})u_{*}$, $a_{21}=(1+\alpha
v_{*})v_{*}$, $a_{22}=\alpha u_{*}v_{*}$, and
$\widehat{\psi}_{\sigma}(k_{x},k_{y})=\mbox{exp}(-\sigma^{2}(k_{x}^{2}+k_{y}^{2})/2)$
is the Fourier transform of the kernel function $\psi_{\sigma}(x,y)$ in two
variables. For the non-trivial solution of the matrix equation (13), the
determinant of the matrix $\mathbf{M}$ has to be equal to $0$, and it leads to
the characteristic equation

${}\lambda^{2}-\mathcal{T}(k)\lambda+\mathcal{D}(k)=0,$ (14)

where
$\mathcal{T}(k)=a_{11}-(\eta/\kappa)u_{*}\exp{(-\sigma^{2}k^{2}/2)}+a_{22}-(1+d)k^{2}$
and
$\mathcal{D}(k)=(a_{11}-(\eta/\kappa)u_{*}\exp{(-\sigma^{2}k^{2}/2)}-k^{2})(a_{22}-dk^{2})-a_{12}a_{21}$
with $k^{2}=k_{x}^{2}+k_{y}^{2}$ with $k$ being the wave number.

From (14), we obtain

${}\lambda_{\pm}(k)=\frac{\mathcal{T}(k)\pm\sqrt{(\mathcal{T}(k))^{2}-4\mathcal{D}(k)}}{2}.$
(15)

For a fixed $\sigma$, we see that $\mathcal{T}(k)<0$ holds for all $k>0$ and
$d$, as $\mbox{tr}(\mathbf{J}_{*})<0$, but we can not conclude anything about
the sign of $\mathcal{D}(k)$. If $\mathcal{D}(k)>0$ holds for all $k>0$ at
given values of $d$ and $\sigma$, then the homogeneous steady-state
$(u_{*},v_{*})$ is stable under the heterogeneous perturbations. Violating the
condition $\mathcal{D}(k)>0$ for some $k>0$ causes an instability in the
homogeneous steady-state $(u_{*},v_{*})$, called Turing instability. In this
case, a Turing pattern can be observed in the system (2) for such parameter
values.

So, for a fixed value of $\sigma$, we focus on finding the critical value
$d=d_{T}$ and the corresponding critical value $k=k_{T}>0$ for which the
Turing bifurcation occurs. At the Turing bifurcation threshold,
$\mathcal{D}(k)=0$ holds for a unique $k$, and in this case $\mathcal{D}(k)$
satisfies

${}\mathcal{D}(k)=0\leavevmode\nobreak\ \mbox{and}\leavevmode\nobreak\
\frac{\partial\mathcal{D}(k)}{\partial k}=0.$ (16)

From the first equation of (16), we get

${}d=\frac{a_{22}(a_{11}-u_{*}e^{-k^{2}\sigma^{2}/2}-k^{2})-a_{12}a_{21}}{k^{2}(a_{11}-u_{*}e^{-k^{2}\sigma^{2}/2}-k^{2})}.$
(17)

After eliminating $d$ between the first and second equations of (16), we
obtain

$\displaystyle{}2a_{22}(a_{11}-u_{*}e^{-\sigma^{2}k^{2}/2}-k^{2})^{2}-a_{12}a_{21}\bigg{(}2a_{11}-4k^{2}-u_{*}(2-k^{2}\sigma^{2})e^{-k^{2}\sigma^{2}/2}\bigg{)}=0.$
(18)

For a fixed $\sigma$, we find the solution $k=k_{T}$ of the equation (18)
numerically, By substituting it into the equation (17), we find the critical
diffusion coefficient $d_{T}$ for the Turing bifurcation. The Turing
bifurcation threshold for the local model can be obtained by calculating
$d=d_{T}$ with $\sigma=0$.

###  4.2 Weakly Nonlinear Analysis

Near the Turing bifurcation threshold, the dynamics of the system (2) change
slowly due to the small variation in the parameter values, and we can study
the pattern formation for such parameters with the help of the amplitude
equations. In this case, we consider three active dominant resonant pairs of
eigenmodes $(\mathbf{k}_{j},-\mathbf{k}_{j})\leavevmode\nobreak\ (j=1,2,3)$
making angles of $2\pi/3$ with $|\mathbf{k}_{j}|=k_{T}$. The solution of the
nonlocal model (2) near the Turing bifurcation threshold $d=d_{T}$ can be
written as the sum of the Fourier plane waves:

${}\begin{pmatrix}u\\\ v\end{pmatrix}=\begin{pmatrix}u_{*}\\\
v_{*}\end{pmatrix}+\sum_{j=1}^{3}[\textbf{A}_{j}(t)\exp(i\mathbf{k}_{j}\cdot\mathbf{r})+\overline{\textbf{A}}_{j}(t)\exp(-i\mathbf{k}_{j}\cdot\mathbf{r})],$
(19)

where $\textbf{A}_{j}=(A_{j,u},A_{j,v})^{T}$ and
$\overline{\textbf{A}}_{j}=(\overline{A}_{j,u},\overline{A}_{j,v})^{T}$ are
the amplitudes associated with the eigenmodes $\mathbf{k}_{j}$ and
$-\mathbf{k}_{j}$, respectively.

Now, we derive the amplitude equations for $\textbf{A}_{j}(j=1,2,3)$. We
expand the time in terms of different time-scales with respect to a small
parameter $\epsilon$ and also expand the bifurcation parameter $d$ and the
other variables $u$ and $v$ as follows:

$\displaystyle t$ $\displaystyle=t_{0}+\epsilon
t_{1}+\epsilon^{2}t_{2}+\cdots,$ (20a) $\displaystyle d$
$\displaystyle=d_{T}+\epsilon d^{(1)}+\epsilon^{2}d^{(2)}+\cdots,$ (20b)
$\displaystyle u$ $\displaystyle=u_{*}+\epsilon
u_{1}+\epsilon^{2}u_{2}+\epsilon^{3}u_{3}+\cdots,$ (20c) $\displaystyle v$
$\displaystyle=v_{*}+\epsilon
v_{1}+\epsilon^{2}v_{2}+\epsilon^{3}v_{3}+\cdots.$ (20d)

The amplitude $\textbf{A}_{j}$ of the spatial pattern evolves on a slow
temporal scale at the initial stage. The derivative $\frac{\partial}{\partial
t_{0}}$ does not have an effect on $\textbf{A}_{j}$ because it corresponds to
the fast time. Therefore, we separate the fast and slow time scales as

${}\frac{\partial}{\partial t}=\epsilon\frac{\partial}{\partial
t_{1}}+\epsilon^{2}\frac{\partial}{\partial t_{2}}+o(\epsilon^{3}).$ (21)

We substitute (20) and (21) into (2) and equate the coefficients of
$\epsilon$, $\epsilon^{2}$ and $\epsilon^{3}$. Comparing the first order
coefficients of $\epsilon$, we obtain

${}\mathbf{L}_{T}\begin{pmatrix}u_{1}\\\
v_{1}\end{pmatrix}\equiv\begin{pmatrix}f_{10}+\Delta-(\eta/\kappa)u_{*}\psi_{\sigma}\ast&f_{01}\\\
g_{10}&g_{01}+d_{T}\Delta\end{pmatrix}\begin{pmatrix}u_{1}\\\
v_{1}\end{pmatrix}=\begin{pmatrix}0\\\ 0\end{pmatrix},$ (22)

where $f_{10}=a_{11}$, $f_{01}=a_{12}$, $g_{10}=a_{21}$, and $g_{01}=a_{22}$.

The solution of the system (22) can be written in the form

${}\begin{pmatrix}u_{1}\\\ v_{1}\end{pmatrix}=\begin{pmatrix}f_{1}\\\
g_{1}\end{pmatrix}\bigg{(}\sum_{j=1}^{3}W_{j}\mbox{exp}(i\mathbf{k}_{j}\cdot\mathbf{r})\bigg{)}+\mbox{c.c.},$
(23)

where

$f_{1}=-f_{01}/f_{c}\leavevmode\nobreak\ \mbox{and}\leavevmode\nobreak\
g_{1}=(f_{10}-(\eta/\kappa)u_{*}e^{-k_{T}^{2}\sigma^{2}/2}-k_{T}^{2})/f_{c}$

with
$f_{c}=\sqrt{(f_{10}-(\eta/\kappa)u_{*}\exp{(-k_{T}^{2}\sigma^{2}/2)}-k_{T}^{2})^{2}+f_{01}^{2}}$.
Here, for given $j$, $W_{j}$ is the modulus of the first order disturbance
term, and c.c. denotes the complex conjugate.

Now, comparing the second order coefficients of $\epsilon$, we get

$\displaystyle{}\mathbf{L}_{T}\begin{pmatrix}u_{2}\\\
v_{2}\end{pmatrix}=\frac{\partial}{\partial t_{1}}$
$\displaystyle\begin{pmatrix}u_{1}\\\ v_{1}\end{pmatrix}-\begin{pmatrix}0&0\\\
0&d^{(1)}\Delta\end{pmatrix}\begin{pmatrix}u_{1}\\\ v_{1}\end{pmatrix}$
$\displaystyle-\begin{pmatrix}f_{20}u_{1}^{2}+f_{11}u_{1}v_{1}+f_{02}v_{1}^{2}-(\eta/\kappa)u_{1}\psi_{\sigma}\ast
u_{1}\\\
g_{20}u_{1}^{2}+g_{11}u_{1}v_{1}+g_{02}v_{1}^{2}\end{pmatrix}=\begin{pmatrix}F_{u}\\\
F_{v}\end{pmatrix},$ (24)

where $f_{20}=0$, $f_{11}=-(1+2\alpha v_{*})$, $f_{02}=-\alpha u_{*}$,
$g_{20}=0$, $g_{11}=(1+2\alpha v_{*})$, and $g_{02}=\alpha u_{*}$.

Our next target is to solve the system (4.2) for $(u_{2},v_{2})^{T}$. The
Fredholm solvability condition ensures the existence of a nontrivial solution
of the non-homogeneous problem (4.2). According to that condition, the right-
hand side of the equation (4.2) must be orthogonal to the zero eigenvectors of
the operator $\mathbf{L}_{T}^{+}$ (the adjoint operator of the operator
$\mathbf{L}_{T}$). Now, the zero eigenvector of the operator $\mathbf{L}_{T}$
is

${}\begin{pmatrix}f_{2}\\\
g_{2}\end{pmatrix}\mbox{exp}(-i\mathbf{k}_{j}\cdot\mathbf{r})+\mbox{c.c.},$
(25)

where

$f_{2}=-g_{10}/g_{c}\leavevmode\nobreak\ \mbox{and}\leavevmode\nobreak\
g_{2}=(f_{10}-(\eta/\kappa)u_{*}e^{-k_{T}^{2}\sigma^{2}/2}-k_{T}^{2})/g_{c},$

with
$g_{c}=\sqrt{(f_{10}-(\eta/\kappa)u_{*}\exp{(-k_{T}^{2}\sigma^{2}/2)}-k_{T}^{2})^{2}+g_{10}^{2}}$.

From the orthogonality condition, we obtain

$(f_{2},g_{2})\begin{pmatrix}F_{u}^{j}\\\
F_{v}^{j}\end{pmatrix}=0,\leavevmode\nobreak\ \leavevmode\nobreak\ (j=1,2,3),$
(26)

where $F_{u}^{j}$ and $F_{v}^{j}$ are the coefficients of
$\mbox{exp}(i\mathbf{k}_{j}\cdot\mathbf{r})$ in $F_{u}$ and $F_{v}$
respectively.

For $j=1$, we get

$\begin{pmatrix}F_{u}^{1}\\\ F_{v}^{1}\end{pmatrix}=\begin{pmatrix}f_{1}\\\
g_{1}\end{pmatrix}\frac{\partial W_{1}}{\partial t_{1}}+\begin{pmatrix}0\\\
g_{1}k_{T}^{2}d^{(1)}\end{pmatrix}W_{1}-2\begin{pmatrix}F_{1}\\\
G_{1}\end{pmatrix}\overline{W}_{2}\overline{W}_{3},$ (27)

where
$F_{1}=f_{20}f_{1}^{2}+f_{11}f_{1}g_{1}+f_{02}g_{1}^{2}-(\eta/\kappa)f_{1}^{2}\exp{(-k_{T}^{2}\sigma^{2}/2)}$
and $G_{1}=g_{20}f_{1}^{2}+g_{11}f_{1}g_{1}+g_{02}g_{1}^{2}$.

After applying the solvability condition, we get

$(f_{1}f_{2}+g_{1}g_{2})\frac{\partial W_{1}}{\partial
t_{1}}=-g_{1}g_{2}k_{T}^{2}d^{(1)}W_{1}+2(f_{2}F_{1}+g_{2}G_{1})\overline{W}_{2}\overline{W}_{3}.$
(28)

Similarly, for $j=2$ and $3$, we obtain the following equations:

$(f_{1}f_{2}+g_{1}g_{2})\frac{\partial W_{2}}{\partial
t_{1}}=-g_{1}g_{2}k_{T}^{2}d^{(1)}W_{2}+2(f_{2}F_{1}+g_{2}G_{1})\overline{W}_{1}\overline{W}_{3},$
(29) $(f_{1}f_{2}+g_{1}g_{2})\frac{\partial W_{3}}{\partial
t_{1}}=-g_{1}g_{2}k_{T}^{2}d^{(1)}W_{3}+2(f_{2}F_{1}+g_{2}G_{1})\overline{W}_{1}\overline{W}_{2}.$
(30)

The solution of the system (4.2) can be written as:

$\displaystyle{}\begin{pmatrix}u_{2}\\\
v_{2}\end{pmatrix}=\begin{pmatrix}X_{0}\\\ Y_{0}\end{pmatrix}$
$\displaystyle+\sum_{j=1}^{3}\begin{pmatrix}X_{j}\\\
Y_{j}\end{pmatrix}\mbox{exp}(i\mathbf{k}_{j}\cdot\mathbf{r})+\sum_{j=1}^{3}\begin{pmatrix}X_{jj}\\\
Y_{jj}\end{pmatrix}\mbox{exp}(i2\mathbf{k}_{j}\cdot\mathbf{r})$
$\displaystyle+\begin{pmatrix}X_{12}\\\
Y_{12}\end{pmatrix}\mbox{exp}(i(\mathbf{k}_{1}-\mathbf{k}_{2})\cdot\mathbf{r})+\begin{pmatrix}X_{23}\\\
Y_{23}\end{pmatrix}\mbox{exp}(i(\mathbf{k}_{2}-\mathbf{k}_{3})\cdot\mathbf{r})$
$\displaystyle+\begin{pmatrix}X_{13}\\\
Y_{13}\end{pmatrix}\mbox{exp}(i(\mathbf{k}_{1}-\mathbf{k}_{3})\cdot\mathbf{r})+\mbox{c.c.}$
(31)

Substituting (4.2) and (23) into (4.2), and collecting the constant terms and
the coefficients of $\mbox{exp}(i\mathbf{k}_{j}\cdot\mathbf{r})$,
$\mbox{exp}(i2\mathbf{k}_{j}\cdot\mathbf{r})$ and
$\mbox{exp}(i(\mathbf{k}_{j}-\mathbf{k}_{m})\cdot\mathbf{r})$ ($1\leq j,m\leq
3$ and $j\neq m$), we find

$\displaystyle\begin{pmatrix}X_{0}\\\ Y_{0}\end{pmatrix}$
$\displaystyle=-2\begin{pmatrix}f_{10}-(\eta/\kappa)u_{*}&f_{01}\\\
g_{10}&g_{01}\end{pmatrix}^{-1}\begin{pmatrix}F_{1}\\\
G_{1}\end{pmatrix}(|W_{1}|^{2}+|W_{2}|^{2}+|W_{3}|^{2})$
$\displaystyle=\begin{pmatrix}\xi_{u0}\\\
\xi_{v0}\end{pmatrix}(|W_{1}|^{2}+|W_{2}|^{2}+|W_{3}|^{2}),$ (32a)
$\displaystyle\begin{pmatrix}X_{j}\\\ Y_{j}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}f_{1}\\\ g_{1}\end{pmatrix}Z_{j},$ (32b)
$\displaystyle\begin{pmatrix}X_{jj}\\\ Y_{jj}\end{pmatrix}$
$\displaystyle=-\begin{pmatrix}f_{10}-4k_{T}^{2}-(\eta/\kappa)u_{*}e^{-4k_{T}^{2}\sigma^{2}/2}&f_{01}\\\
g_{10}&g_{01}-4d_{T}k_{T}^{2}\end{pmatrix}^{-1}\begin{pmatrix}F_{1}\\\
G_{1}\end{pmatrix}W_{j}^{2}$ $\displaystyle=\begin{pmatrix}\xi_{u1}\\\
\xi_{v1}\end{pmatrix}W_{j}^{2},$ (32c) $\displaystyle\begin{pmatrix}X_{jm}\\\
Y_{jm}\end{pmatrix}$
$\displaystyle=-2\begin{pmatrix}f_{10}-3k_{T}^{2}-(\eta/\kappa)u_{*}e^{-3k_{T}^{2}\sigma^{2}/2}&f_{01}\\\
g_{10}&g_{01}-3d_{T}k_{T}^{2}\end{pmatrix}^{-1}\begin{pmatrix}F_{1}\\\
G_{1}\end{pmatrix}W_{j}\overline{W}_{m}$
$\displaystyle=\begin{pmatrix}\xi_{u2}\\\
\xi_{v2}\end{pmatrix}W_{j}\overline{W}_{m}.$ (32d)

Now, comparing the third order coefficients of $\epsilon$, we obtain

$\displaystyle{}\mathbf{L}_{T}\begin{pmatrix}u_{3}\\\ v_{3}\end{pmatrix}=$
$\displaystyle\begin{pmatrix}\frac{\partial u_{2}}{\partial
t_{1}}+\frac{\partial u_{1}}{\partial t_{2}}\\\ \frac{\partial v_{2}}{\partial
t_{1}}+\frac{\partial v_{1}}{\partial
t_{2}}\end{pmatrix}-\begin{pmatrix}0&0\\\
0&d^{(2)}\Delta\end{pmatrix}\begin{pmatrix}u_{1}\\\
v_{1}\end{pmatrix}-\begin{pmatrix}0&0\\\
0&d^{(1)}\Delta\end{pmatrix}\begin{pmatrix}u_{2}\\\ v_{2}\end{pmatrix}$
$\displaystyle-\begin{pmatrix}2f_{20}u_{1}u_{2}+f_{11}(u_{1}v_{2}+u_{2}v_{1})+2f_{02}v_{1}v_{2}-(\eta/\kappa)u_{1}\psi_{\sigma}\ast
u_{2}-(\eta/\kappa)u_{2}\psi_{\sigma}\ast u_{1}\\\
2g_{20}u_{1}u_{2}+g_{11}(u_{1}v_{2}+u_{2}v_{1})+2g_{02}v_{1}v_{2}\end{pmatrix}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
-\begin{pmatrix}f_{30}u_{1}^{3}+f_{21}u_{1}^{2}v_{1}+f_{12}u_{1}v_{1}^{2}+f_{03}v_{1}^{3}\\\
g_{30}u_{1}^{3}+g_{21}u_{1}^{2}v_{1}+g_{12}u_{1}v_{1}^{2}+g_{03}v_{1}^{3}\end{pmatrix},$
(33)

where $f_{30}=0$, $f_{21}=0$, $f_{12}=-\alpha$, $f_{03}=0$, $g_{30}=0$,
$g_{21}=0$, $g_{12}=\alpha$, and $g_{03}=0$.

After applying the Fredholm solvability condition to the system (4.2) we
obtain:

$\displaystyle(f_{1}f_{2}+g_{1}g_{2})\bigg{(}\frac{\partial W_{1}}{\partial
t_{2}}+\frac{\partial Z_{1}}{\partial t_{1}}\bigg{)}=$ $\displaystyle-
g_{1}g_{2}k_{T}^{2}(d^{(1)}Z_{1}+d^{(2)}W_{1})+2(f_{2}F_{1}+g_{2}G_{1})(\overline{W}_{2}\overline{Z}_{3}+\overline{W}_{3}\overline{Z}_{2})$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+((f_{2}F_{2}+g_{2}G_{2})|W_{1}|^{2}+(f_{2}F_{3}+g_{2}G_{3})(|W_{2}|^{2}+|W_{3}|^{2}))W_{1},$
(34a) $\displaystyle(f_{1}f_{2}+g_{1}g_{2})\bigg{(}\frac{\partial
W_{2}}{\partial t_{2}}+\frac{\partial Z_{2}}{\partial t_{1}}\bigg{)}=$
$\displaystyle-
g_{1}g_{2}k_{T}^{2}(d^{(1)}Z_{2}+d^{(2)}W_{2})+2(f_{2}F_{1}+g_{2}G_{1})(\overline{W}_{1}\overline{Z}_{3}+\overline{W}_{3}\overline{Z}_{1})$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+((f_{2}F_{2}+g_{2}G_{2})|W_{2}|^{2}+(f_{2}F_{3}+g_{2}G_{3})(|W_{1}|^{2}+|W_{3}|^{2}))W_{2},$
(34b) $\displaystyle(f_{1}f_{2}+g_{1}g_{2})\bigg{(}\frac{\partial
W_{3}}{\partial t_{2}}+\frac{\partial Z_{3}}{\partial t_{1}}\bigg{)}=$
$\displaystyle-
g_{1}g_{2}k_{T}^{2}(d^{(1)}Z_{3}+d^{(2)}W_{3})+2(f_{2}F_{1}+g_{2}G_{1})(\overline{W}_{1}\overline{Z}_{2}+\overline{W}_{2}\overline{Z}_{1})$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+((f_{2}F_{2}+g_{2}G_{2})|W_{3}|^{2}+(f_{2}F_{3}+g_{2}G_{3})(|W_{1}|^{2}+|W_{2}|^{2}))W_{3},$
(34c)

where

$\displaystyle F_{2}$
$\displaystyle=2f_{20}f_{1}(\xi_{u0}+\xi_{u1})+f_{11}(f_{1}(\xi_{v0}+\xi_{v1})+g_{1}(\xi_{u0}+\xi_{u1}))+2f_{02}g_{1}(\xi_{v0}+\xi_{v1})$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
-(\eta/\kappa)f_{1}(\xi_{u0}+\xi_{u1}e^{-4k_{T}^{2}\sigma^{2}/2})-(\eta/\kappa)f_{1}(\xi_{u0}+\xi_{u1})e^{-k_{T}^{2}\sigma^{2}/2}+3F_{4},$
$\displaystyle G_{2}$
$\displaystyle=2g_{20}f_{1}(\xi_{u0}+\xi_{u1})+g_{11}(f_{1}(\xi_{v0}+\xi_{v1})+g_{1}(\xi_{u0}+\xi_{u1}))+2g_{02}g_{1}(\xi_{v0}+\xi_{v1})+3G_{4},$
$\displaystyle F_{3}$
$\displaystyle=2f_{20}f_{1}(\xi_{u0}+\xi_{u2})+f_{11}(f_{1}(\xi_{v0}+\xi_{v2})+g_{1}(\xi_{u0}+\xi_{u2}))+2f_{02}g_{1}(\xi_{v0}+\xi_{v2})$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
-(\eta/\kappa)f_{1}(\xi_{u0}+\xi_{u2}e^{-3k_{T}^{2}\sigma^{2}/2})-(\eta/\kappa)f_{1}(\xi_{u0}+\xi_{u2})e^{-k_{T}^{2}\sigma^{2}/2}+6F_{4},$
$\displaystyle G_{3}$
$\displaystyle=2g_{20}f_{1}(\xi_{u0}+\xi_{u2})+g_{11}(f_{1}(\xi_{v0}+\xi_{v2})+g_{1}(\xi_{u0}+\xi_{u2}))+2g_{02}g_{1}(\xi_{v0}+\xi_{v2})+6G_{4},$
$\displaystyle F_{4}$
$\displaystyle=f_{1}^{3}f_{30}+f_{1}^{2}g_{1}f_{21}+f_{1}g_{1}^{2}f_{12}+g_{1}^{3}f_{03},$
$\displaystyle G_{4}$
$\displaystyle=f_{1}^{3}g_{30}+f_{1}^{2}g_{1}g_{21}+f_{1}g_{1}^{2}g_{12}+g_{1}^{3}g_{03}.$

From the equations (19), (20), (23), and (4.2), we find the relation between
the amplitudes as

${}\begin{pmatrix}A_{j,u}\\\
A_{j,v}\end{pmatrix}=\epsilon\begin{pmatrix}f_{1}\\\
g_{1}\end{pmatrix}W_{j}+\epsilon^{2}\begin{pmatrix}f_{1}\\\
g_{1}\end{pmatrix}Z_{j}+o(\epsilon^{3}).$ (35)

Now, we focus on the amplitude equations corresponding to the $u$-component.
For the notational simplicity, we denote $A_{j,u}$ as $A_{j}$. Therefore, from
(35), we have

$A_{j}=\epsilon
f_{1}W_{j}+\epsilon^{2}f_{1}Z_{j}+o(\epsilon^{3}),\leavevmode\nobreak\ 1\leq
j\leq 3.$ (36)

Then the amplitude equation with three nodes $A_{j}(j=1,2,3)$ is given by

$\displaystyle\tau_{0}\frac{\partial A_{1}}{\partial t}$ $\displaystyle=\mu
A_{1}+h_{0}\overline{A}_{2}\overline{A}_{3}-(m_{1}|A_{1}|^{2}+m_{2}(|A_{2}|^{2}+|A_{3}|^{2}))A_{1},$
(37a) $\displaystyle\tau_{0}\frac{\partial A_{2}}{\partial t}$
$\displaystyle=\mu
A_{2}+h_{0}\overline{A}_{1}\overline{A}_{3}-(m_{1}|A_{2}|^{2}+m_{2}(|A_{1}|^{2}+|A_{3}|^{2}))A_{2},$
(37b) $\displaystyle\tau_{0}\frac{\partial A_{3}}{\partial t}$
$\displaystyle=\mu
A_{3}+h_{0}\overline{A}_{1}\overline{A}_{2}-(m_{1}|A_{3}|^{2}+m_{2}(|A_{1}|^{2}+|A_{2}|^{2}))A_{3},$
(37c)

where $\mu=(d_{T}-d)/d_{T}$ is a normalized distance to onset of the Turing
bifurcation threshold,
$\tau_{0}=(f_{1}f_{2}+g_{1}g_{2})/d_{T}k_{T}^{2}g_{1}g_{2}$ is a typical
relaxation time,
$h_{0}=2(f_{2}F_{1}+g_{2}G_{1})/d_{T}k_{T}^{2}f_{1}g_{1}g_{2}$,
$m_{1}=-(f_{2}F_{2}+g_{2}G_{2})/d_{T}k_{T}^{2}f_{1}^{2}g_{1}g_{2}$ and
$m_{2}=-(f_{2}F_{3}+g_{2}G_{3})/d_{T}k_{T}^{2}f_{1}^{2}g_{1}g_{2}$. Note that
the amplitude equations corresponding to the $v$-component are also like the
equations for the $u$-component. For the case of $v$-component, the expression
of the coefficients $\mu$ and $\tau_{0}$ remain unaltered and the other
coefficients are give by
$h_{0}=2(f_{2}F_{1}+g_{2}G_{1})/d_{T}k_{T}^{2}g_{1}^{2}g_{2}$,
$m_{1}=-(f_{2}F_{2}+g_{2}G_{2})/d_{T}k_{T}^{2}g_{1}^{3}g_{2}$ and
$m_{2}=-(f_{2}F_{3}+g_{2}G_{3})/d_{T}k_{T}^{2}g_{1}^{3}g_{2}$.

Now, each of the amplitudes in equation (37) can be decomposed into the mode
$\rho_{j}=|A_{j}|(j=1,2,3)$ and a corresponding phase $\phi_{j}$. Substituting
$A_{j}=\rho_{j}\mbox{exp}(i\phi_{j})$ into equations of (37) and separating
the real and imaginary parts, we obtain the following differential equations
in real variables:

$\displaystyle\tau_{0}\frac{\partial\Phi}{\partial t}$
$\displaystyle=-h_{0}\frac{\rho_{1}^{2}\rho_{2}^{2}+\rho_{2}^{2}\rho_{3}^{2}+\rho_{3}^{2}\rho_{1}^{2}}{\rho_{1}\rho_{2}\rho_{3}}\sin\Phi,$
(38a) $\displaystyle\tau_{0}\frac{\partial\rho_{1}}{\partial t}$
$\displaystyle=\mu\rho_{1}+h_{0}\rho_{2}\rho_{3}\cos\Phi-
m_{1}\rho_{1}^{3}-m_{2}(\rho_{2}^{2}+\rho_{3}^{2})\rho_{1},$ (38b)
$\displaystyle\tau_{0}\frac{\partial\rho_{2}}{\partial t}$
$\displaystyle=\mu\rho_{2}+h_{0}\rho_{1}\rho_{3}\cos\Phi-
m_{1}\rho_{2}^{3}-m_{2}(\rho_{1}^{2}+\rho_{3}^{2})\rho_{2},$ (38c)
$\displaystyle\tau_{0}\frac{\partial\rho_{3}}{\partial t}$
$\displaystyle=\mu\rho_{3}+h_{0}\rho_{1}\rho_{2}\cos\Phi-
m_{1}\rho_{3}^{3}-m_{2}(\rho_{1}^{2}+\rho_{2}^{2})\rho_{3},$ (38d)

where $\Phi=\phi_{1}+\phi_{2}+\phi_{3}$. The system (38) has always an
equilibrium point with the components $\rho_{1}=\rho_{2}=\rho_{3}=0$. If this
trivial equilibrium point is stable, then following (19), we can not find any
non-homogeneous stationary pattern (Turing pattern) for the nonlocal model
(2). Therefore, for the generation of a non-homogeneous stationary pattern,
the trivial equilibrium point of the system (38) has to be unstable.

Now, we find all the equilibrium points of the system (38) and determine their
stabilities. For non-zero amplitudes $\rho_{1}$, $\rho_{2}$, $\rho_{3}$ and
$\tau_{0}>0$, the solution corresponding to $\Phi=0$ ($H_{0}$ pattern) is
stable if $h_{0}>0$, and the solution corresponding to $\Phi=\pi$ ($H_{\pi}$
pattern) is stable if $h_{0}<0$. The mode equations for the stable solutions
of the equation (38a) is given by

$\displaystyle\tau_{0}\frac{\partial\rho_{1}}{\partial t}$
$\displaystyle=\mu\rho_{1}+|h_{0}|\rho_{2}\rho_{3}-m_{1}\rho_{1}^{3}-m_{2}(\rho_{2}^{2}+\rho_{3}^{2})\rho_{1},$
(39a) $\displaystyle\tau_{0}\frac{\partial\rho_{2}}{\partial t}$
$\displaystyle=\mu\rho_{2}+|h_{0}|\rho_{1}\rho_{3}-m_{1}\rho_{2}^{3}-m_{2}(\rho_{1}^{2}+\rho_{3}^{2})\rho_{2},$
(39b) $\displaystyle\tau_{0}\frac{\partial\rho_{3}}{\partial t}$
$\displaystyle=\mu\rho_{3}+|h_{0}|\rho_{1}\rho_{2}-m_{1}\rho_{3}^{3}-m_{2}(\rho_{1}^{2}+\rho_{2}^{2})\rho_{3}.$
(39c)

######  Theorem 6.

Suppose $\tau_{0}>0$. Then system (2) admits four kinds of solutions, and they
are

  * (I)

Homogeneous solution: $\rho_{1}=\rho_{2}=\rho_{3}=0$. It is stable for
$\mu<\mu_{2}=0$ and unstable for $\mu>\mu_{2}$.

  * (II)

Stripe pattern: $\rho_{1}=\sqrt{\mu/m_{1}}\neq 0,\rho_{2}=\rho_{3}=0$. It is
stable for $\mu>\mu_{3}=h_{0}^{2}m_{1}/(m_{2}-m_{1})^{2}$ and unstable for
$\mu<\mu_{3}$.

  * (III)

Hexagonal pattern ($H_{0}$ or $H_{\pi}$):

$\rho_{1}=\rho_{2}=\rho_{3}=\frac{|h_{0}|\pm\sqrt{h_{0}^{2}+4(m_{1}+2m_{2})\mu}}{2(m_{1}+2m_{2})}\equiv\rho^{\pm}.$

These amplitudes exist if $\mu>\mu_{1}=-h_{0}^{2}/(4(m_{1}+2m_{2}))$. The
solution $\rho^{+}$ is stable for
$\mu<\mu_{4}=h_{0}^{2}(2m_{1}+m_{2})/(m_{2}-m_{1})^{2}$, and $\rho^{-}$ is
always unstable.

  * (IV)

Mixed pattern:

$\rho_{1}=\frac{|h_{0}|}{m_{2}-m_{1}},\leavevmode\nobreak\
\leavevmode\nobreak\ \rho_{2}=\rho_{3}=\sqrt{\frac{\mu-
m_{1}\rho_{1}^{2}}{m_{1}+m_{2}}},$

with $m_{2}>m_{1}$, $\mu>m_{1}\rho_{1}^{2}$, and they are always unstable.

###### Proof.

Let $(\delta\rho_{1},\delta\rho_{2},\delta\rho_{3})$ be a perturbation of
$(\rho_{1},\rho_{2},\rho_{3})$. Putting these perturbations in the system (39)
and ignoring the second and higher order terms, we obtain the following matrix
equation

$\tau_{0}\frac{d}{dt}\begin{bmatrix}\delta\rho_{1}\\\ \delta\rho_{2}\\\
\delta\rho_{3}\end{bmatrix}=\begin{bmatrix}\omega_{11}&\omega_{12}&\omega_{13}\\\
\omega_{21}&\omega_{22}&\omega_{23}\\\
\omega_{31}&\omega_{32}&\omega_{33}\end{bmatrix}\begin{bmatrix}\delta\rho_{1}\\\
\delta\rho_{2}\\\
\delta\rho_{3}\end{bmatrix}\equiv\mathbf{W}\begin{bmatrix}\delta\rho_{1}\\\
\delta\rho_{2}\\\ \delta\rho_{3}\end{bmatrix},$ (40)

where $\omega_{11}=\mu-3m_{1}\rho_{1}^{2}-m_{2}(\rho_{2}^{2}+\rho_{3}^{2})$,
$\omega_{12}=|h_{0}|\rho_{3}-2m_{2}\rho_{1}\rho_{2}$,
$\omega_{13}=|h_{0}|\rho_{2}-2m_{2}\rho_{1}\rho_{3}$,
$\omega_{21}=|h_{0}|\rho_{3}-2m_{2}\rho_{1}\rho_{2}$,
$\omega_{22}=\mu-3m_{1}\rho_{2}^{2}-m_{2}(\rho_{1}^{2}+\rho_{3}^{2})$,
$\omega_{23}=|h_{0}|\rho_{1}-2m_{2}\rho_{2}\rho_{3}$,
$\omega_{31}=|h_{0}|\rho_{2}-2m_{2}\rho_{1}\rho_{3}$,
$\omega_{32}=|h_{0}|\rho_{1}-2m_{2}\rho_{2}\rho_{3}$ and
$\omega_{33}=\mu-3m_{1}\rho_{3}^{2}-m_{2}(\rho_{1}^{2}+\rho_{2}^{2})$.

Case (I): If $\rho_{1}=\rho_{2}=\rho_{3}=0$, then $\mathbf{W}$ becomes

$\mathbf{W}=\begin{bmatrix}\mu&0&0\\\ 0&\mu&0\\\ 0&0&\mu\end{bmatrix}$

Therefore, the homogeneous solution is stable if $\mu<\mu_{2}=0$ and unstable
for $\mu>\mu_{2}$.
Case (II): If $\rho_{1}=\sqrt{\mu/m_{1}}\neq 0,\rho_{2}=\rho_{3}=0$, then the
eigenvalues of the matrix $\mathbf{W}$ are

$-2\mu,\mu\bigg{(}1-\frac{m_{2}}{m_{1}}\bigg{)}\pm|h_{0}|\sqrt{\frac{\mu}{m_{1}}}.$

Therefore, the stripe pattern is stable if $m_{2}>m_{1}>0$ and
$\mu>\mu_{3}=h_{0}^{2}m_{1}/(m_{2}-m_{1})^{2}$. Similarly, we can prove the
rest of the cases. ∎

##  5 Numerical Results

This section validates all the obtained theoretical results through numerical
computations. First, we present some bifurcations for the temporal model, then
the spatio-temporal pattern for the local model, followed by the patterns for
the nonlocal model. For the numerical simulations of the local and nonlocal
models, we have used the SHARCNET (www.sharcnet.ca) high performance
computational facilities to minimize the time in computations.

Figure 2: (Color online) Three parametric bifurcation diagram for the temporal
model (1). Bifurcation surfaces are transcritical (TC), saddle-node (SN), and
Hopf (HB). The intersection of two surfaces SN and HB is the Bogdanov-Takens
(BT) bifurcation.

###  5.1 Temporal bifurcations

Figure 2 represents a bifurcation diagram for the temporal model in three-
dimensional parameter space. We find three local bifurcation surfaces, namely
saddle-node (SN), transcritical (TC) and Hopf (HB) bifurcation surfaces. Among
them, the transcritical bifurcation surface is a plane $\kappa=1$. The saddle-
node bifurcation surface lies below the transcritical bifurcation plane. The
Hopf bifurcation surface lies both above and below the plane $\kappa=1$. The
Hopf bifurcation and saddle-node bifurcation surfaces touch each other along
the Bogdanov-Takens (BT) bifurcation curve (marked in red). There is no
coexistence equilibrium point for parameter values lying below the SN surface.
We find two coexistence equilibrium points for parameter values lying between
SN and TC surfaces. One of the coexistence equilibrium points disappears
through a transcritical bifurcation as $\kappa$ crosses the value 1 from
below.

SN, TC, and HB are the global bifurcations for the temporal model; however,
our primary goal is to match the analytical pattern selections with the
numerical non-homogeneous stationary patterns in the Turing domain. Here, HB
plays an essential role as for the existence of Turing bifurcation, the non-
trivial equilibrium point of the temporal model has to be Hopf-stable. So, we
choose the parameter values so that the non-trivial equilibrium point of the
temporal model satisfies the necessary condition for the Turing bifurcation.
We discuss all the parameter values in the next section.

###  5.2 Patterns for the local model

In this subsection, we report extensive numerical simulations for the local
and nonlocal models and present some qualitative results. We use periodic
boundary conditions for both the models (local and nonlocal) in a square
domain $[-25,25]\times[-25,25]$ with space steps $dx=dy=0.25$ and a time step
$dt=0.01$. We have verified the numerical artifact of the patterns by choosing
smaller space and time discretization steps for the simulations, and the
presented patterns are unaltered. For finding the numerical solution of the
local model, we have used Euler’s method for the time derivative and a second-
order finite-difference scheme to approximate the diffusion terms.
Heterogeneous perturbations around the homogeneous steady-state have been used
for the initial conditions for the species $u$ and $v$. We run the simulations
in each of the upcoming patterns until the solution characteristics stop
changing over time. We choose $\kappa=0.65$ and $\alpha=10$ to be fixed
parameters (unless stated otherwise) and $\eta$ as a bifurcation parameter.

We first analyze the pattern formation for the local model. Taking the fixed
parameter values, as mentioned above, we choose the bifurcation parameter
value $\eta=0.92$. For this specific choice of $\eta$, the temporal system has
two coexisting equilibrium points: one is a saddle point and the other is
locally asymptotically stable. As mentioned earlier, we choose the locally
asymptotically stable equilibrium point for finding the Turing bifurcation.
For the local model, i.e., $\sigma=0$, we solve the equation (18) for $k$, and
we find the critical wave number as $k_{T}=0.871$ and the corresponding Turing
bifurcation threshold as $d_{T}=0.2715$. For the local model, the region
$d<d_{T}$ is the Turing domain, where the stationary Turing pattern exists,
and in the other region $d>d_{T}$, the homogeneous solution remains stable
under the heterogeneous perturbation.

Figure 3: (Color online) Bifurcation diagram of the patterns obtained through
weakly nonlinear analysis for the local model with the parameter values
$\kappa=0.65$, $\alpha=10$ and $\eta=0.92$. $HS$: homogeneous solution;
$H_{0}$: hexagonal patterns with $\Phi=0$; $H_{\pi}$: hexagonal patterns with
$\Phi=\pi$; solid curve: stable state; dashed line: unstable state.

Now, for $\eta=0.92$, we find all the parameter values associated with the
weakly nonlinear analysis and the values rounded up to three decimal places
are $m_{1}=6.035$, $m_{2}=121.566$, $h_{0}=3.050$, $\tau_{0}=3.536$,
$\mu_{1}=-0.009$, $\mu_{2}=0$ $\mu_{3}=0.004$, $\mu_{4}=0.093$. Figure 3
depicts the stability of the stationary patterns for the local model obtained
through the weakly nonlinear analysis. We find the diffusion coefficients
$d_{j}$ ($j=1,2,3$ and $4$) corresponding to each $\mu_{j}$ by using the
relation $d_{j}=(1-\mu_{j})d_{T}$ and the values are $d_{1}=0.274$,
$d_{2}=d_{T}$, $d_{3}=0.2704$, and $d_{4}=0.2462$. Next, we discuss all the
results based on $d_{j}$. From Fig. 3, we can see that the homogeneous
solution ($HS$) is stable for $d>d_{2}$. On the other hand, for
$d_{3}<d<d_{2}$, two types of hexagonal patterns ($H_{0}$ and $H_{\pi}$) may
occur in the system. Since, $\tau_{0}$ and $h_{0}$ both are positive, the
hexagonal pattern $H_{0}$ (hot-spot) is stable, and the other one $H_{\pi}$
(cold-spot) is unstable (by Theorem 6).

We find a hexagonal Turing pattern ($H_{0}$) for $d=0.271$ near the Turing
bifurcation threshold in the Turing domain [see Fig. 4(4(a))]. Here, the
diffusion parameter $d$ lies in the bistable region $(d_{3},d_{2})$, hence the
hexagonal pattern ($H_{0}$) is stable, and the other patterns are unstable, as
predicted through weakly nonlinear analysis. Figure 4(4(a)) depicts the
stationary Turing pattern for the prey population. In this case, we observe a
cold-spot pattern ($H_{\pi}$) corresponding to the predator populations
(excluded in this paper). The relation of the patterns depends on the sign of
$f_{1}$ and $g_{1}$, obtained through the weakly nonlinear analysis. If both
of them are of the same sign, then the sign of $h_{0}$ remains the same for
both components, and hence by Theorem 6 both populations follow the same
hexagonal patterns. Otherwise, the populations follow the opposite patterns
because of the opposite signs in $h_{0}$. In our considered parameter values,
$f_{1}$ is positive, and $g_{1}$ is negative [see Table 1], and hence the
weakly nonlinear analysis is consistent with the pattern obtained through
numerical simulation. The same type of negative correlation has been observed
for all the upcoming stationary patterns. Moreover, we can find the amplitude
equation corresponding to the predator population and it also predicts the
same trend as discussed here.

(a)

(b)

(c)

Figure 4: (Color online) Turing patterns for the local model with $\eta=0.92$,
$\kappa=0.65$ and $\alpha=10$: (a) $d=0.271$, (b) $d=0.25$, and (c) $d=0.2$.
In all the plots, $x$-axis and $y$-axis are horizontal and vertical,
respectively.

According to Theorem 6, the hexagonal pattern $H_{0}$ is stable for
$d_{4}<d<d_{2}$ and the stripe pattern is stable for $d<d_{3}$ [see Fig. 3].
Therefore, the system has a bistable region $d_{4}<d<d_{3}$, where both
patterns (hot-spot and stripe) are stable. We choose the diffusion parameter
$d=0.26$ in the bistable region $(d_{4},d_{3})$, and the corresponding stripe
pattern is shown in Fig. 4(4(b)). We will find the hot-spot pattern for
different initial conditions.

Theorem 6 predicts a stable stripe pattern for $d<d_{3}$. In the meantime,
Fig. 4(4(c)) depicts a cold-spot ($H_{\pi}$) pattern for $d=0.2$ in the prey
population. This is called a reentry of a hexagonal pattern, and it occurs for
the parameter values far from a Turing bifurcation threshold. In this case,
some other primary slave modes become active, along with the critical
wavenumber obtained through the linear stability analysis. We cannot neglect
them in the derivation of the amplitude equations; rather, we should consider
them in the amplitude equations [33], and it would be an interesting extension
of this work in the future. Combining all the results, with a decrease in the
diffusion parameter value $d$, the number of stationary patches in a spatial
domain of fixed size increases.

Table 1: Computed parameter values of the local and nonlocal models for $\eta=0.92$ corresponding to the prey population. $\sigma$ | $d_{T}$ | $k_{T}$ | $f_{1}$ | $g_{1}$ | $m_{1}$ | $m_{2}$ | $h_{0}$ | $\tau_{0}$ | $d_{3}$ | $d_{4}$
---|---|---|---|---|---|---|---|---|---|---
0.00 | 0.2715 | 0.871 | 0.757 | -0.654 | 6.035 | 121.566 | 3.050 | 3.536 | 0.2704 | 0.2462
0.25 | 0.2665 | 0.879 | 0.757 | -0.654 | 7.892 | 123.659 | 2.981 | 3.536 | 0.2651 | 0.2419
0.50 | 0.2521 | 0.906 | 0.756 | -0.654 | 14.398 | 130.878 | 2.712 | 3.529 | 0.2502 | 0.2303
0.75 | 0.2310 | 0.953 | 0.753 | -0.657 | 27.485 | 144.896 | 2.099 | 3.498 | 0.2289 | 0.2162
1.00 | 0.2076 | 1.028 | 0.745 | -0.667 | 49.105 | 166.414 | 0.927 | 3.410 | 0.2070 | 0.2042
1.25 | 0.1881 | 1.122 | 0.728 | -0.686 | 74.819 | 189.901 | -0.670 | 3.260 | 0.1876 | 0.1860
1.50 | 0.1758 | 1.208 | 0.707 | -0.707 | 95.247 | 208.345 | -2.071 | 3.110 | 0.1702 | 0.1523

###  5.3 Patterns for the nonlocal model

In this subsection, we focus on the numerical results of the analytical
findings and the stationary Turing patterns for the nonlocal model. We have
used the trapezoidal rule to find the numerical value of the convolution term.
Figure 5 depicts the Turing bifurcation curves for the nonlocal model for
different values of $\sigma$. With an increase in the parameter value of
$\sigma$, the Turing bifurcation curve shifts downwards, however, it converges
to a curve inside the considered domain in Fig. 5. Theoretically, one can find
the asymptotic curve of this Turing bifurcation curve (cf. [7]), but here we
restrict the pattern formation for $\sigma\leq 1.5$. A detailed study of the
pattern formation of the considered nonlocal model (2) for higher values of
$\sigma$ is one of the important issues for future work.

Figure 5: (Color online) Turing and temporal-Hopf bifurcation curves for the
local and nonlocal models.

(a)

(b)

(c)

Figure 6: (Color online) Turing patterns for the nonlocal model with
$\eta=0.92$, $\kappa=0.65$, $\alpha=10$ and $\sigma=0.5$: (a) $d=0.251$, (b)
$d=0.22$ and (c) $d=0.19$. In all the plots, $x$-axis and $y$-axis are
horizontal and vertical, respectively.

We fix $\eta=0.92$ and $\sigma=0.5$. For this choice of the parameter values,
the Turing bifurcation threshold becomes $d_{T}=0.2521$ and the critical
wavenumber $k_{T}=0.906$ [see Table 1]. Some of the computational parameter
values for the weakly nonlinear analysis are listed in Table 1. Following
Theorem 6, the cold spot pattern is stable for $d_{4}<d<d_{2}(=d_{T})$, the
stripe pattern is stable for $d<d_{3}$. Numerical solutions of the nonlocal
model for three different diffusion parameter values are plotted in Fig. 2,
and they show a similar trend to the one obtained with the local model.
However, the amplitudes of each of the patterns are different from the local
model. Due to a similar reason discussed earlier, we can not predict the
reentry of the hexagonal pattern (hot-spot) through the weakly nonlinear
analysis described in this article.

(a)

(b)

(c)

Figure 7: (Color online) Turing patterns for the nonlocal model with
$\eta=0.92$, $\kappa=0.65$, $\alpha=10$ and $\sigma=1.5$: (a) $d=0.172$, (b)
$d=0.16$, and (c) $d=0.1$. In all the plots, $x$-axis and $y$-axis are
horizontal and vertical, respectively.

Further, if we increase the parameter values of $\sigma$ in the mentioned
range with fixing $\eta=0.92$, the parameter $h_{0}$ changes its sign from
positive to negative [see Table 1], but the other parameters remain with the
same sign. This sign change in the parameter $h_{0}$ is reflected in the
resulting pattern, e.g., we consider $\sigma=1.5$, and the computed parameter
values involved in the weakly nonlinear analysis is mentioned in Table 1. The
numerical simulation result shows a cold-spot pattern for the prey population
for the diffusion parameter $d=0.172$ [see Fig. 7(7(a))], near the Turing
bifurcation threshold in the Turing domain. Figures 7(7(b)) and (7(c)) show
the cold-spot pattern for the lower diffusion parameters far from the Turing
bifurcation threshold. Numerical simulation results show only the cold-spot
pattern for the nonlocal model with $\sigma=1.5$.

From Table 1, we see that the interval $(d_{4},d_{3})$ for the existence of
the labyrinthine pattern becomes more narrow with an increase in the magnitude
of $\sigma$. As a result, we observe that the thresholds corresponding to the
hexagonal patterns ($d_{3}$ and the threshold for reentry of the hexagonal
pattern) come closer, and then the labyrinthine pattern disappears. As a
matter of fact, the amplitude equations fail to predict the selection of
stripe pattern theoretically; however, it predicts the stationary hexagonal
patterns near the Turing bifurcation threshold.

##  6 Conclusions

In this paper, we have studied the pattern formation for the local and
nonlocal models, exemplifying our results on prey-predator models with hunting
cooperation. Implicit analytical conditions have been derived for the
existence of different temporal bifurcations for the diffusionless model, in
particular, transcritical, saddle-node, and Hopf bifurcations. Moreover, a
three-dimensional bifurcation diagram has been plotted to show these
bifurcations numerically. The Turing bifurcation condition for the nonlocal
model has been discussed analytically. Fixing the temporal parameters and the
kernel parameter $\sigma$, we have derived the analytical equation for finding
the critical wave number corresponding to the Turing bifurcation and this
helps in finding the critical diffusion threshold. In particular, the Turing
bifurcations for the local model are obtainable by taking the kernel parameter
$\sigma=0$.

The non-homogeneous stationary patterns exist for the local and nonlocal
models in the Turing domain. The results of pattern selections can be
identified through the weakly nonlinear analysis described in this manuscript.
We have derived the amplitude equations for the nonlocal model near the Turing
bifurcation threshold. Likewise, in the Turing bifurcation, the amplitude
equations corresponding to the local model can be derived by taking
$\sigma=0$. We have compared the numerical simulation results with the
theoretical results for the local model. They agreed near the Turing
bifurcation threshold but differed far from the threshold. If we move from
being close to the Turing bifurcation threshold to being further away the
stationary pattern corresponding to the prey population changes from cold-spot
$\rightarrow$ labyrinthine $\rightarrow$ hot-spot. The reentry of the
hexagonal pattern (hot-spot) can not be predicted through the weakly nonlinear
analysis. We have derived the amplitude equations for both prey and predator
populations. With the help of weakly nonlinear analysis, we have shown that
the prey and predator populations are negatively correlated with each other,
and hence, the predator population shifts their patterns from hot-spot
$\rightarrow$ labyrinthine $\rightarrow$ cold-spot when the parameter shifts
being close to the Turing bifurcation threshold to being further away.

The pattern selection of the nonlocal model has also been discussed in this
manuscript. With an increase in the nonlocal parameter $\sigma$, the Turing
bifurcation curve shifts downwards but does not disappear; rather, it
saturates to a curve. The weakly nonlinear analysis also predicts a tapered
behaviour in the region where a labyrinthine pattern exists for the nonlocal
model. Furthermore, numerical simulation results suggest the thresholds
corresponding to reentry of the hexagonal pattern move towards the Turing
bifurcation curve and eventually coincide with it. As a result, the stripe
labyrinthine pattern disappears, and only the cold-spot hexagonal patterns
exist for the nonlocal model. The reentry of the hexagonal pattern can not be
identified through our described weakly nonlinear analysis due to neglecting
some other primary slave modes, and it could be a possible extension of this
work in the future. Finally, we note that the number of non-homogeneous
stationary patches in a spatial domain of fixed size corresponding to a small
value of the nonlocal parameter $\sigma$ is smaller compared to the higher
values of $\sigma$.

## Acknowledgements

SP and RM are grateful to the NSERC and the CRC Program for their support. RM
is also acknowledging support of the BERC 2022-2025 program and Spanish
Ministry of Science, Innovation and Universities through the Agencia Estatal
de Investigacion (AEI) BCAM Severo Ochoa excellence accreditation
SEV-2017-0718 and the Basque Government fund AI in BCAM EXP. 2019/00432. This
research was enabled in part by support provided by SHARCNET (www.sharcnet.ca)
and Compute Canada (www.computecanada.ca).

## References

  * [1] M. Kot, Elements of mathematical ecology, Cambridge University Press, 2001.
  * [2] J. A. Sherratt, Periodic traveling waves in integrodifferential equations for nonlocal dispersal, SIAM Journal on Applied Dynamical Systems 13 (4) (2014) 1517–1541.
  * [3] M. Banerjee, V. Volpert, Prey-predator model with a nonlocal consumption of prey, Chaos: An Interdisciplinary Journal of Nonlinear Science 26 (8) (2016) 083120\.
  * [4] S. Pal, S. Ghorai, M. Banerjee, Analysis of a prey-predator model with non-local interaction in the prey population, Bulletin of Mathematical Biology 80 (2018) 906–925.
  * [5] A. Bayliss, V. Volpert, Complex predator invasion waves in a Holling-Tanner model with nonlocal prey interaction, Physica D: Nonlinear Phenomena 346 (2017) 37–58.
  * [6] S. Pal, S. Petrovskii, S. Ghorai, M. Banerjee, Spatiotemporal pattern formation in 2d prey-predator system with nonlocal intraspecific competition, Communications in Nonlinear Science and Numerical Simulation 93 (2021) 105478\.
  * [7] S. Pal, S. Ghorai, M. Banerjee, Effect of kernels on spatio-temporal patterns of a non-local prey-predator model, Mathematical Biosciences 310 (2019) 96–107.
  * [8] M. Banerjee, M. Kuznetsov, O. Udovenko, V. Volpert, Nonlocal reaction-diffusion equations in biomedical applications, Acta Biotheoretica 70 (2) (2022) 1–28.
  * [9] S. Pal, R. Melnik, Nonlocal models in the analysis of brain neurodegenerative protein dynamics with application to Alzheimer’s disease, Preprint: Scientific Reports (2022).
  * [10] R. Eftimie, G. de Vries, M. Lewis, Weakly nonlinear analysis of a hyperbolic model for animal group formation, Journal of mathematical biology 59 (1) (2009) 37–74.
  * [11] D. Sytnyk, R. Melnik, Mathematical models with nonlocal initial conditions: An exemplification from quantum mechanics, Mathematical and Computational Applications 26 (4) (2021) 73.
  * [12] V. Volpert, Elliptic partial differential equations, Vol. 2, Springer, 2014.
  * [13] F. Paquin-Lefebvre, On the weakly nonlinear analysis of coupled bulk-surface reaction-diffusion systems: theory, numerics, applications, Ph.D. thesis, University of British Columbia (2020).
  * [14] S. Pal, R. Melnik, The role of astrocytes in Alzheimer’s disease progression, IWBBIO-2022 (2022).
  * [15] W. W. Murdoch, C. J. Briggs, R. M. Nisbet, Consumer-resource dynamics (mpb-36), in: Consumer-Resource Dynamics (MPB-36), Princeton University Press, 2013.
  * [16] H. I. Freedman, Deterministic mathematical models in population ecology, Vol. 57, Marcel Dekker Incorporated, 1980.
  * [17] J. M. Fryxell, A. R. Sinclair, G. Caughley, Wildlife ecology, conservation, and management, John Wiley & Sons, 2014.
  * [18] P. A. Abrams, L. R. Ginzburg, The nature of predation: prey dependent, ratio dependent or neither?, Trends in Ecology & Evolution 15 (2000) 337–341.
  * [19] M. Lu, C. Xiang, J. Huang, H. Wang, Bifurcations in the diffusive bazykin model, Journal of Differential Equations 323 (2022) 280–311.
  * [20] M. T. Alves, F. M. Hilker, Hunting cooperation and Allee effects in predators, Journal of theoretical biology 419 (2017) 13–22.
  * [21] A. M. Turing, The chemical basis of morphogenesis, Bulletin of Mathematical Biology 52 (1) (1952) 153–197.
  * [22] J. D. Murray, Mathematical Biology I. An Introduction, Springer, 2002.
  * [23] C. B. Price, Computing with active media: pattern formation and analysis by reaction-diffusion systems, neural fields and coupled map lattices, K.U. Leuven (1993).
  * [24] C. B. Price, Turing hexagons in pure cubic nonlinear fields, Physics Letters A 194 (1994) 385–390.
  * [25] W. Eckhaus, Studies in non-linear stability theory, Springer, Berlin, Heidelberg (1965).
  * [26] L. A. Segel, S. A. Levin, Application of nonlinear stability theory to the study of the effects of diffusion on predator-prey interactions, in: AIP conference proceedings, Vol. 27, American Institute of Physics, 1976, pp. 123–152.
  * [27] X.-C. Zhang, G.-Q. Sun, Z. Jin, Spatial dynamics in a predator-prey model with Beddington-DeAngelis functional response, Physical Review E 85 (2012) 021924.
  * [28] M. Ipsen, F. Hynne, P. Sørensen, Amplitude equations for reaction-diffusion systems with a Hopf bifurcation and slow real modes, Physica D: Nonlinear Phenomena 136 (2000) 66–92.
  * [29] G. H. Gunaratne, Q. Ouyang, H. L. Swinney, Pattern formation in the presence of symmetries, Physical Review E 50 (1994) 2802.
  * [30] X. Li, G. Hu, Z. Feng, Pattern dynamics in a spatial predator-prey model with nonmonotonic response function, International Journal of Bifurcation and Chaos 28 (2018) 1850077.
  * [31] F. Lutscher, E. Pachepsky, M. A. Lewis, The effect of dispersal patterns on stream populations, Siam Review 47 (4) (2005) 749–772.
  * [32] L. Perko, Diferential equations and dynamical systems, Springer, USA (2000).
  * [33] Q. Ouyang, Pattern formation in reaction-diffusion systems, Shanghai: Shanghai Sci-Tech Education Publishing House)(in Chinese) (2000).