#  The spatially inhomogeneous Hopf bifurcation induced by memory delay in a
memory-based diffusion system111Partially supported by the Natural Science
Foundation of China (Nos.11971143, 12071105), and Natural Science Foundation
of Zhejiang Province of China (No.LY19A010010).

Yongli Song1, Yahong Peng2222Corresponding author, Email:
<EMAIL_ADDRESS>, Tonghua Zhang3
1 Department of Mathematics, Hangzhou Normal University, Hangzhou 311121,
China
2 Department of Mathematics, Donghua University, Shanghai, 201620, China
3 Department of Mathematics, Swinburne University of Technology, Hawthorn,
VIC3122, Australia

###### Abstract

The memory-based diffusion systems have wide applications in practice. Hopf
bifurcations are observed from such systems. To meet the demand for computing
the normal forms of the Hopf bifurcations of such systems, we develop an
effective new algorithm where the memory delay is treated as the perturbation
parameter. To illustrate the effectiveness of the algorithm, we consider a
diffusive predator-prey system with memory-based diffusion and Holling type-II
functional response. By employing this newly developed procedure, we
investigate the direction and stability of the delay-induced mode-1 and mode-2
Hopf bifurcations. Numerical simulations confirm our theoretical findings,
that is the existence of stable spatially inhomogeneous periodic solutions
with mode-1 and mode-2 spatial patterns, and the transition from the unstable
mode-$2$ spatially inhomogeneous periodic solution to the stable mode-$1$
spatially inhomogeneous periodic solution.

Keywords: Spatial memory; delay; Hopf bifurcation; normal form; periodic
solution
MSC2010: 35B10; 37L10; 37G05

##  1 Introduction

The study of delay-induced Hopf bifurcations is a very active research area in
recent years, and has attracted many researchers’ attention due to its
importance. In the study of nonlinear dynamical systems that may exist rich
dynamics, such as Hopf bifurcation, the normal form plays an import role as it
can be used to determine the direction and stability of the Hopf bifurcation.
Thus, many researchers have dedicated themselves to the development of ways
calculating normal forms of the nonlinear dynamical systems. For example,
Hassard, Kazarinoff and Wan [9] and Faria and Magalhaes [7] proposed an
algorithm for computing the normal form of the Hopf bifrucations for delay
differential systems (without diffusion), respectively; Faria [6] developed an
algorithm for reaction-diffusion systems with time delays, but which only
appear in the reaction terms; most recently Song and his colleagues [28, 23]
developed some algorithms for reaction-diffusion systems with both time delay
and nonlocal spatial average. It is now known that in the classical delay
reaction-diffusion systems, the first delay-induced Hopf bifurcation is often
homogeneous and the bifurcating periodic oscillation is spatially homogeneous
[32, 33, 31, 27, 29, 26, 30, 21, 10, 17, 14, 3, 2, 5, 15, 12, 13, 11]. And the
aforementioned algorithms were developed for such scenarios.

However, in applications, the delay reaction-diffusion systems may not be
homogeneous in space. For example, very recently, to better understand how
memory influences the animal movement, Shi et al. [19] proposed a modified
Fick’s law where in addition to the negative gradient of the density
distribution function at the present time, there is a directed movement toward
the negative gradient of the density distribution function at past time, and
then based on the modified Fick’s law, they formulated a single species
diffusive model with spatial memory. Since then, there have recently been an
increasing activities and interests on the study of the dynamics of the
reaction-diffusion equations with spatial memory [18, 19, 1, 16, 20, 22]. For
this kind of equations, the phenomenon that the first delay-induced Hopf
bifurcation is inhomogeneous has been observed by [18, 24]. Based on the
assumption that the prey has no memory or cognition, we introduced the spatial
memory to the diffusive predator-prey model and proposed the following model
with random and memory-based diffusions subject to Neumann boundary conditions
on one-dimensional spatial domain $(0,\ell\pi)$ with $\ell\in\mathbb{R}^{+}$
is

$\begin{cases}\frac{\partial u(x,t)}{\partial
t}=d_{11}u_{xx}(x,t)+f(u(x,t),v(x,t)),&0<x<\ell\pi,t>0,\vspace{0.1cm}\\\
\frac{\partial v(x,t)}{\partial
t}=d_{22}v_{xx}(x,t)-d_{21}(v(x,t)u_{x}(x,t-\tau))_{x}&\\\ \hskip
43.36243pt+g(u(x,t),v(x,t)),&0<x<\ell\pi,t>0,\vspace{0.1cm}\\\
u_{x}(0,t)=u_{x}(\ell\pi,t)=v_{x}(0,t)=v_{x}(\ell\pi,t)=0,&t\geq
0,\end{cases}$ (1.1)

where $d_{11}$ and $d_{22}$ are the Fickian diffusion coefficients, $d_{12}$
and $d_{21}$ are the memory-based diffusion coefficients, the time delay,
denoted by $\tau>0$ represents the averaged memory period of predator. It has
been shown that delay-induced first Hopf bifurcations are often spatially
inhomogeneous [22].

Compared with the classical reaction-diffusion system with delay, system (1.1)
has two main characteristics: $(i)$ the delay appears in the diffusion term;
$(ii)$ the diffusion term is not linear. Therefore, algorithms, in the
existing literatures such as [9, 6, 28, 23], that were developed for analysing
the Hopf bifurcation of the reaction-diffusion systems, can not be applied to
the system (1.1). To bridge the gap, in this paper, we aim to develop an
algorithm that can calculate the normal form for analysing the Hopf
bifurcation which may occur in a spatially inhomogeneous scenario and be
induced by the memory delay in (1.1). And then we apply this algorithm to a
diffusive predator-prey system with memory-based diffusion and Holling type-II
functional response. More precisely, we organize the rest of our paper as
follows. In Section 2, we derive the algorithm of calculating normal form of
spatially Hopf bifurcation induced by memory delay for system (1.1). In
Section 3, we study the direction and stability of delay-induced mode-1 and
mode-2 Hopf bifurcations in the diffusive predator-prey system with memory-
based diffusion and Holling type-II functional response and the numerical
simulations are used to illustrate the theoretical results. Finally, we give a
brief discussion and conclude our study in Section 4. For the convenience of
discussion, throughout the paper, we let $\mathbb{N}$ represent the set of
natural numbers, and $\mathbb{N}_{0}=\mathbb{N}\cup\\{0\\}$ represent the set
of all nonnegative integers.

##  2 Calculation of the normal form of the Hopf bifurcation

###  2.1 Characteristic equation at the equilibrium and some basic assumptions

Without loss of generality, assume that $E_{*}=(u_{*},v_{*})$ is a constant
coexistence (positive) equilibrium of system (1.1). Then straightforward
calculation gives the linearization of (1.1) at $(u_{*},v_{*})$:

$\left(\begin{array}[]{c}u_{t}(x,t)\\\
v_{t}(x,t)\end{array}\right)=D_{1}\left(\begin{array}[]{c}u_{xx}(x,t)\\\
v_{xx}(x,t)\end{array}\right)+D_{2}\left(\begin{array}[]{c}u_{xx}(x,t-\tau)\\\
v_{xx}(x,t-\tau)\end{array}\right)+A\left(\begin{array}[]{c}u(x,t)\\\
v(x,t)\end{array}\right),$ (2.1)

where

$D_{1}=\left(\begin{array}[]{cc}d_{11}&0\\\
0&d_{22}\end{array}\right),~{}D_{2}=\left(\begin{array}[]{cc}0&0\\\
-d_{21}v_{*}&0\end{array}\right),~{}A=\left(\begin{array}[]{cc}a_{11}&a_{12}\\\
a_{21}&a_{22}\\\ \end{array}\right),$ (2.2)

and

$a_{11}=\frac{\partial f(u_{*},v_{*})}{\partial u},~{}a_{12}=\frac{\partial
f(u_{*},v_{*})}{\partial v},~{}a_{21}=\frac{\partial g(u_{*},v_{*})}{\partial
u},~{}a_{22}=\frac{\partial g(u_{*},v_{*})}{\partial v}.$

Thus, the characteristic equation of (2.1) is

$\prod\limits_{n\in\mathbb{N}_{0}}\Gamma_{n}(\lambda)=0,$ (2.3)

where $\Gamma_{n}(\lambda)=\mathrm{det}\left(\mathcal{M}_{n}(\lambda)\right)$.
Notice

$\mathcal{M}_{n}(\lambda)=\lambda
I_{2}+(n/\ell)^{2}D_{1}+(n/\ell)^{2}e^{-\lambda\tau}D_{2}-A,$ (2.4)

with $I_{2}$ is the identity matrix of size 2. Then we have

$\Gamma_{n}(\lambda)=\mathrm{det}\left(\mathcal{M}_{n}(\lambda)\right)=\lambda^{2}-T_{n}\lambda+\widetilde{J}_{n}(\tau)=0,$
(2.5)

where

$\begin{split}T_{n}&=Tr(A)-Tr(D_{1})(n/\ell)^{2},\\\
\widetilde{J}_{n}(\tau)&=d_{11}d_{22}(n/\ell)^{4}-\left(d_{11}a_{22}+d_{22}a_{11}+d_{21}v_{*}a_{12}e^{-\lambda\tau}\right)(n/\ell)^{2}+Det(A),\end{split}$
(2.6)

with $Tr(A)=a_{11}+a_{22}$, $Tr(D_{1})=d_{11}+d_{22}$ and
$Det(A)=a_{11}a_{22}-a_{12}a_{21}$.

For the Hopf bifurcation, we assume that at $\tau=\tau_{c}$, Eq.(2.5) has a
pair of purely imaginary roots $\pm i\omega_{n_{c}},\omega_{n_{c}}>0$ for
$n=n_{c}\in\mathbb{N}$ and all other eigenvalues have negative real part. Let
$\lambda(\tau)=\alpha(\tau)+i\beta(\tau)$ be a pair of roots of Eq.(2.5) near
$\tau=\tau_{c}$ satisfying $\alpha\left(\tau_{c}\right)=0$ and
$\beta\left(\tau_{c}\right)=\omega_{n_{c}}$. Then, it is easy to verify the
transversality condition

$\left.\frac{d\mathrm{Re}(\lambda(\tau))}{d\tau}\right|_{\tau=\tau_{c}}\neq
0.$

###  2.2 Algorithm for computing the normal form

We first introduce some conventional notations: define the two-dimensional
real-valued Sobolev space

$\mathscr{X}=\left\\{U=\left(U_{1},U_{2}\right)^{T}\in\left(W^{2,2}(0,\ell\pi)\right)^{2},\frac{\partial
U_{1}}{\partial x}=0,\frac{\partial U_{2}}{\partial
x}=0,~{}x=0,\ell\pi\right\\},$

with the inner product defined by

$\left[U,V\right]=\int_{0}^{\ell\pi}U^{T}Vdx,~{}\mbox{for}~{}U,V\in\mathscr{X},$

and let $\mathscr{C}:=C\left([-1,0];\mathscr{X}\right)$ be the Banach space of
continuous mappings from $[-1,0]$ to $\mathscr{X}$ .

Then, we take a small perturbation of $\tau_{c}$ by setting
$\tau=\tau_{c}+\mu,|\mu|\ll 1$ such that $\mu=0$ corresponds to the Hopf
bifurcation value for Eq.(1.1). Here, $D_{j},j=1,2$, are given by (2.2).

Now shift $E_{*}$ to the origin by setting

$U(x,t)=\left(U_{1}(x,t),U_{2}(x,t)\right)^{T}=(u(x,t),v(x,t))^{T}-(u_{*},v_{*})^{T},$

normalize the delay by rescaling the time variable $t\rightarrow t/\tau$, and
rewrite $U(t)$ for $U(x,t)$ and $U_{t}\in\mathscr{C}$ for
$U_{t}(\theta)=U(x,t+\theta),-1\leq\theta\leq 0$. Then, the system (1.1)
becomes the compact form

$\frac{dU(t)}{dt}=d(\mu)\Delta(U_{t})+L(\mu)(U_{t})+F(U_{t},\mu),$ (2.7)

where for
$\varphi=\left(\varphi^{(1)},\varphi^{(2)}\right)^{T}\in\mathscr{C}$,
$d(\mu)\Delta,L(\mu):\mathscr{C}\rightarrow\mathscr{X}$,
$F:\mathscr{C}\times\mathbb{R}^{2}\rightarrow\mathscr{X}$ are given,
respectively, by

$d(\mu)\Delta(\varphi)=d_{0}\Delta(\varphi)+F^{d}(\varphi,\mu),~{}L(\mu)(\varphi)=(\tau_{c}+\mu)A\varphi(0),$

and

$F(\varphi,\mu)=\left(\tau_{c}+\mu\right)\left(\begin{array}[]{c}f\left(\varphi^{(1)}(0)+u_{*},\varphi^{(2)}0)+v_{*}\right)\vspace{0.2cm}\\\
g\left(\varphi^{(1)}(0)+u_{*},\varphi^{(2)}(0)+v_{*}\right)\end{array}\right)-L(\mu)(\varphi),$
(2.8)

where

$d_{0}\Delta(\varphi)=\tau_{c}D_{1}\varphi_{xx}(0)+\tau_{c}D_{2}^{c}\varphi_{xx}(-1),$
$\begin{array}[]{lll}F^{d}(\varphi,\mu)&=&-d_{21}\left(\tau_{c}+\mu\right)\left(\begin{array}[]{c}0\vspace{0.2cm}\\\
\varphi^{(1)}_{x}(-1)\varphi^{(2)}_{x}(0)+\varphi^{(1)}_{xx}(-1)\varphi^{(2)}(0)\end{array}\right)\vspace{0.3cm}\\\
&&+\mu\left(\begin{array}[]{c}d_{11}\varphi^{(1)}_{xx}(0)\vspace{0.2cm}\\\
-d_{21}v_{*}\varphi^{(1)}_{xx}(-1)+d_{22}\varphi^{(2)}_{xx}(0)\end{array}\right).\end{array}$
(2.9)

In what follows, we assume that $F(\varphi,\mu)$ is $C^{k},k\geq 3,$ smooth
with respect to $\varphi$ and $\mu$. Noticing that $\mu$ is the perturbation
parameter and treated as a variable in the calculation of normal forms, we
denote $L_{0}(\varphi)=\tau_{c}A\varphi(0)$ and rewrite (2.7) as the following

$\frac{dU(t)}{dt}=d_{0}\Delta(U_{t})+L_{0}(U_{t})+\widetilde{F}(U_{t},\mu),$
(2.10)

by separating the linear terms from the nonlinear terms,

$\widetilde{F}(\varphi,\mu)=\mu
A\varphi(0)+F(\varphi,\mu)+F^{d}(\varphi,\mu).$ (2.11)

Then the characteristic equation for the linearized system of (2.10)

$\frac{dU(t)}{dt}=d_{0}\Delta(U_{t})+L_{0}(U_{t})$ (2.12)

is

$\prod\limits_{n\in\mathbb{N}_{0}}\widetilde{\Gamma}_{n}(\lambda)=0,$ (2.13)

where
$\widetilde{\Gamma}_{n}(\lambda)=\mathrm{det}\left(\widetilde{\mathcal{M}}_{n}(\lambda)\right)$
with

$\widetilde{\mathcal{M}}_{n}(\lambda)=\lambda
I_{2}+\tau_{c}(n/\ell)^{2}D_{1}+\tau_{c}(n/\ell)^{2}e^{-\lambda}D_{2}-\tau_{c}A.$
(2.14)

Comparing (2.14) with (2.4), we know that Eq.(2.13) has a pair of purely
imaginary roots $\pm i\omega_{c}$ for $n=n_{c}\in\mathbb{N}$, and all other
eigenvalues have negative real parts, where
$\omega_{c}=\tau_{c}\omega_{n_{c}}$.

It is well known that the eigenvalue problem

$-\gamma^{\prime\prime}=\mu\gamma,~{}~{}x\in(0,\ell\pi);~{}~{}\gamma^{\prime}(0)=\gamma^{\prime}(\ell\pi)=0$

has eigenvalues $\mu_{n}=(n/\ell)^{2},n\in\mathbb{N}_{0}$, with the
corresponding normalized eigenfunctions

$\gamma_{n}(x)=\frac{\cos\left(\frac{nx}{\ell}\right)}{\|\cos\left(\frac{nx}{\ell}\right)\|_{2,2}}=\left\\{\begin{array}[]{ll}\frac{1}{\sqrt{\ell\pi}},&~{}\mathrm{for}~{}n=0,\vspace{0.2cm}\\\
\frac{\sqrt{2}}{\sqrt{\ell\pi}}\cos\left(\frac{nx}{\ell}\right),&~{}\mathrm{for}~{}n\neq
0.\end{array}\right.$

Let $\beta_{n}^{(j)}=\gamma_{n}(x)e_{j},j=1,2$, where $e_{j}$ are the unit
coordinate vectors of $\mathbb{R}^{2}$. Then the eigenfunctions
$\left\\{\beta_{n}^{(j)}\right\\}_{n=1}^{\infty}$ form an orthonormal basis
for $\mathscr{X}$.

Let
$\mathscr{B}_{n}=\mbox{span}\left\\{\left[v(\cdot),\beta_{n}^{(j)}\right]\beta_{n}^{(j)}|~{}v\in\mathscr{C},j=1,2\right\\}$.
Then it is easy to verify that

$L_{0}(\mathscr{B}_{n})\subset\mbox{span}\left\\{\beta_{n}^{(1)},\beta_{n}^{(2)}\right\\},n\in\mathbb{N}_{0}.$

Assume that $z_{t}(\theta)\in C=C\left([-1,0],\mathbb{R}^{2}\right)$ and

$z_{t}^{T}(\theta)\left(\begin{array}[]{c}\beta_{n}^{(1)}\vspace{0.1cm}\\\
\beta_{n}^{(2)}\end{array}\right)\in\mathscr{B}_{n}.$

Then, on $\mathscr{B}_{n}$, the linearized equation (2.12) is equivalent to
the following functional differential equation (FDE) in $C$:

$\dot{z}(t)=L_{0}^{d}\left(z_{t}(\theta)\right)+L_{0}(z_{t}(\theta)),$ (2.15)

where

$L_{0}^{d}\left(z_{t}(\theta)\right)=\tau_{c}\left(\begin{array}[]{cc}-d_{11}(n/\ell)^{2}&0\\\
0&-d_{22}(n/\ell)^{2}\end{array}\right)z_{t}(0)+\tau_{c}\left(\begin{array}[]{cc}0&0\\\
d_{21}v_{*}(n/\ell)^{2}&0\end{array}\right)z_{t}(-1).$

The characteristic equation of linear system (2.15) is the same as the one
given in (2.13).

Define $\eta_{n}(\theta)\in BV([-1,0],\mathbb{R}^{2})$ such that

$\int^{0}_{-1}d\eta_{n}(\theta)\varphi(\theta)=L_{0}^{d}(\varphi(\theta))+L_{0}(\varphi(\theta)),~{}~{}\varphi\in
C,$

and use the adjoint bilinear form on $C^{*}\times C$,
$C^{*}=C([0,1],\mathbb{R}^{2*})$, where $\mathbb{R}^{2*}$ is a 2-dimensional
space of row vectors, as follows

$\langle\psi(s),\varphi(\theta)\rangle_{n}=\psi(0)\varphi(0)-\int^{0}_{-1}\int^{\theta}_{0}\psi(\xi-\theta)d\eta_{n}(\theta)\varphi(\xi)d\xi,\
\ \ \textmd{for}\ \psi\in C^{*},\varphi\in C.$

Let $\Lambda=\\{i\omega_{c},-i\omega_{c}\\}$. Denote the generalized
eigenspace of (2.15) associated with $\Lambda$ by $P$ and the corresponding
adjoint space by $P^{*}$. Then, by the adjoint theory of functional
differential equation [8], $C$ can be decomposed as $C=P\oplus Q$, where
$Q=\\{\varphi\in C:\langle\psi,\varphi\rangle=0,\forall\psi\in P^{*}\\}$.
Choose the bases $\Phi(\theta)$ and $\Psi(s)$ of $P$ and $P^{*}$,
respectively, as follows

$\Phi(\theta)=\left(\phi(\theta),\overline{\phi}(\theta)\right),~{}\Psi(s)=\mbox{col}\left(\psi^{T}(s),\overline{\psi}^{T}(s)\right),$

such that $\langle\Psi,\Phi\rangle_{n_{c}}=I_{2},$ where

$\phi(\theta)=\left(\begin{array}[]{c}\phi_{1}(\theta)\vspace{0.2cm}\\\
\phi_{2}(\theta)\end{array}\right)=\phi
e^{i\omega_{c}\theta},~{}\psi(s)=\left(\begin{array}[]{c}\psi_{1}(\theta)\vspace{0.2cm}\\\
\psi_{2}(\theta)\end{array}\right)=\psi e^{-i\omega_{c}s},$

and

$\phi=\left(\begin{array}[]{c}1\vspace{0.2cm}\\\
\frac{i\omega_{n_{c}}+(n_{c}/\ell)^{2}d_{11}-a_{11}}{a_{12}}\end{array}\right),~{}~{}\psi=\eta\left(\begin{array}[]{c}1\vspace{0.2cm}\\\
\frac{a_{12}}{i\omega_{n_{c}}+(n_{c}/\ell)^{2}d_{22}-a_{22}}\end{array}\right),$

with

$\eta=\frac{i\omega_{n_{c}}+(n_{c}/\ell)^{2}d_{22}-a_{22}}{i\omega_{n_{c}}+(n_{c}/\ell)^{2}d_{11}-a_{11}+i\omega_{n_{c}}+(n_{c}/\ell)^{2}d_{22}-a_{22}+\tau_{c}a_{12}d_{21}v_{*}(n_{c}/\ell)^{2}e^{-i\omega_{c}}}.$

Using the decomposition $C=P\oplus Q$, the phase space $\mathscr{C}$ for (2.7)
can be decomposed as

$\mathscr{C}=\mathcal{P}\oplus\mathcal{Q},~{}\mathcal{P}=\mbox{Im}\pi,~{}\mathcal{Q}=\mbox{Ker}\pi,$

where $\pi:\mathscr{C}\rightarrow\mathcal{P}$ is the projection operator
defined by

$\begin{array}[]{lll}\pi(\phi)&=&\Phi(\theta)\left\langle\Psi(0),\left(\begin{array}[]{c}\left[\phi(\cdot),\beta_{n_{c}}^{(1)}\right]\vspace{0.2cm}\\\
\left[\phi(\cdot),\beta_{n_{c}}^{(2)}\right]\end{array}\right)\right\rangle\gamma_{n_{c}}(x).\end{array}$

Following [6] and [25], we define
$\mathscr{C}_{0}^{1}=\left\\{\phi\in\mathscr{C}:\dot{\phi}\in\mathscr{C},\phi(0)\in\mbox{dom}(d\Delta)\right\\}$
and let

$z_{x}=\left(z_{1}(t)\gamma_{n_{c}}(x),z_{2}(t)\gamma_{n_{c}}(x)\right)^{T},\Phi(\theta)=\left(\phi_{n_{c}}(\theta),\overline{\phi}_{n_{c}}(\theta)\right).$

For $\varphi(\theta)\in\mathscr{C}_{0}^{1}$, we have the following
decomposition

$\varphi(\theta)=\Phi(\theta)z_{x}+w,~{}~{}w=(w^{(1)},w^{(2)})^{T}\in\mathscr{C}_{0}^{1}\cap\mbox{Ker}\pi:=\mathscr{Q}^{1}.$

Following the notations in [6], we define

$X_{0}(\theta)=\left\\{\begin{array}[]{ll}0,&-1\leq\theta<0,\vspace{0.2cm}\\\
1,&\theta=0.\end{array}\right.$

and then

$\pi\left(X_{0}(\theta)\widetilde{F}_{2}\left(\Phi(\theta)z_{x},0\right)\right)\vspace{0.2cm}\\\
=\Phi(\theta)\Psi(0)\left(\begin{array}[]{c}\left[\widetilde{F}_{2}\left(\Phi(\theta)z_{x},0\right),\beta_{n_{c}}^{(1)}\right]\vspace{0.2cm}\\\
\left[\widetilde{F}_{2}\left(\Phi(\theta)z_{x},0\right),\beta_{n_{c}}^{(2)}\right]\end{array}\right)\gamma_{n_{c}}(x).$
(2.16)

Let $z=(z_{1}(t),z_{2}(t))^{T}.$ Then system (2.10) can be decomposed as a
system of abstract ODEs on $\mathbb{R}^{4}\times\mbox{Ker}\pi$:

$\begin{cases}\dot{z}=Bz+\Psi(0)\left(\begin{array}[]{c}\left[\widetilde{F}\left(\Phi(\theta)z_{x}+w,\mu\right),\beta_{n_{c}}^{(1)}\right]\vspace{0.1cm}\\\
\left[\widetilde{F}\left(\Phi(\theta)z_{x}+w,\mu\right),\beta_{n_{c}}^{(2)}\right]\end{array}\right),\vspace{0.1cm}\\\
\dot{w}=A_{\mathcal{Q}^{1}}w+(I-\pi)X_{0}(\theta)\widetilde{F}\left(\Phi(\theta)z_{x}+w,\mu\right),\end{cases}$
(2.17)

where $B=\mbox{diag}\left\\{i\omega_{c},-i\omega_{c}\right\\},$
$A_{\mathcal{Q}^{1}}:\mathcal{Q}^{1}\rightarrow\mbox{Ker}\pi$ is defined by

$A_{\mathscr{Q}^{1}}w=\dot{w}+X_{0}(\theta)\left(L_{0}(w)+L_{0}^{d}(w)-\dot{w}(0)\right).$

Consider the formal Taylor expansion

$\widetilde{F}(\varphi,\mu)=\sum\limits_{j\geq
2}\frac{1}{j!}\widetilde{F}_{j}(\varphi,\mu),~{}~{}F(\varphi,\mu)=\sum\limits_{j\geq
2}\frac{1}{j!}F_{j}(\varphi,\mu),~{}~{}F^{d}(\varphi,\mu)=\sum\limits_{j\geq
2}\frac{1}{j!}F_{j}^{d}(\varphi,\mu).$

From (2.11), we have

$\widetilde{F}_{2}(\varphi,\mu)=2\mu
A\varphi(0)+F_{2}(\varphi,\mu)+F_{2}^{d}(\varphi,\mu)$ (2.18)

and

$\widetilde{F}_{j}(\varphi,\mu)=F_{j}(\varphi,\mu)+F_{j}^{d}(\varphi,\mu),j=3,4,\cdots.$
(2.19)

Then (2.17) is written as

$\left\\{\begin{array}[]{l}\dot{z}=Bz+\sum\limits_{j\geq
2}\frac{1}{j!}f^{1}_{j}(z,w,\mu),\vspace{0.2cm}\\\
\dot{w}=A_{\mathscr{Q}^{1}}w+\sum\limits_{j\geq
2}\frac{1}{j!}f^{2}_{j}(z,w,\mu),\end{array}\right.$

where

$f^{1}_{j}(z,w,\mu)=\Psi(0)\left(\begin{array}[]{c}\left[\widetilde{F}_{j}\left(\Phi(\theta)z_{x}+w,\mu\right),\beta_{n_{c}}^{(1)}\right]\vspace{0.1cm}\\\
\left[\widetilde{F}_{j}\left(\Phi(\theta)z_{x}+w,\mu\right),\beta_{n_{c}}^{(2)}\right]\end{array}\right),$
(2.20)
$f^{2}_{j}(z,w,\mu)=(I-\pi)X_{0}(\theta)\widetilde{F}_{j}\left(\Phi(\theta)z_{x}+w,\mu\right).$
(2.21)

In terms of the normal form theory of partial functional differential
equations [6], after a recursive transformation of variables of the form

$(z,w)=(\widetilde{z},\widetilde{w})+\frac{1}{j!}\left(U_{j}^{1}(\widetilde{z},\mu),U_{j}^{2}(\widetilde{z},\mu)(\theta)\right),j\geq
2,$ (2.22)

where $z,\widetilde{z}\in\mathbb{R}^{2},w,\widetilde{w}\in\mathscr{Q}^{1}$ and
$U_{j}^{1}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{2},U_{j}^{2}:\mathbb{R}^{3}\rightarrow\mathscr{Q}^{1}$
are homogeneous polynomials of degree $j$ in $\widetilde{z}$ and $\mu,$ the
flow on the local center manifold for (2.10) can be written as

$\dot{z}=Bz+\sum\limits_{j\geq 2}\frac{1}{j!}g_{j}^{1}(z,0,\mu),$ (2.23)

which is the normal form as in the usual sense for ODEs.

Following [6] and [7], we have

$g_{2}^{1}(z,0,\mu)=\mbox{Proj}_{\mbox{Ker}(M_{2}^{1})}f_{2}^{1}(z,0,\mu),$

and

$g_{3}^{1}(z,0,\mu)=\mbox{Proj}_{\mbox{Ker}(M_{3}^{1})}\widetilde{f}_{3}^{1}(z,0,\mu)=\mbox{Proj}_{S}\widetilde{f}_{3}^{1}(z,0,0)+O(\mu^{2}|z|),$
(2.24)

where $\widetilde{f}_{3}^{1}(z,0,\mu)$ is vector, where elements are
polynomials of degree 3 in $(z,\mu)$ obtained from (2.22) after performing the
change of variables, and is determined by (2.34),

$\mbox{Ker}\left(M_{2}^{1}\right)=\mbox{Span}\left\\{\left(\begin{array}[]{l}\mu
z_{1}\\\ 0\end{array}\right),\left(\begin{array}[]{l}0\\\ \mu
z_{2}\end{array}\right)\right\\},$
$\mbox{Ker}\left(M_{3}^{1}\right)=\mbox{Span}\left\\{\left(\begin{array}[]{l}z_{1}^{2}z_{2}\\\
0\end{array}\right),\left(\begin{array}[]{l}\mu^{2}z_{1}\\\
0\end{array}\right),\left(\begin{array}[]{l}0\\\
z_{1}z_{2}^{2}\end{array}\right),\left(\begin{array}[]{l}0\\\
\mu^{2}z_{2}\end{array}\right)\right\\}$

and

$S=\mbox{Span}\left\\{\left(\begin{array}[]{l}z_{1}^{2}z_{2}\\\
0\end{array}\right),\left(\begin{array}[]{l}0\\\
z_{1}z_{2}^{2}\end{array}\right)\right\\}.$

For notational convenience, in what follows we let

$\mathcal{H}\left(\alpha
z_{1}^{q_{1}}z_{2}^{q_{2}}\mu\right)=\left(\begin{array}[]{c}\alpha
z_{1}^{q_{1}}z_{2}^{q_{2}}\mu\vspace{0.2cm}\\\
\overline{\alpha}z_{1}^{q_{2}}z_{2}^{q_{1}}\mu\end{array}\right),~{}\alpha\in\mathbb{C}.$

We then calculate $\bf{g_{j}^{1}(z,0,\mu)}$.

####  2.2.1 Calculation of $\bf{g_{2}^{1}(z,0,\mu)}$

From (2.9), we have

$F^{d}_{2}(\varphi,\mu)=F^{d}_{20}(\varphi)+\mu F^{d}_{21}(\varphi),$ (2.25)

and

$F^{d}_{3}(\varphi,\mu)=\mu
F^{d}_{31}(\varphi),~{}F^{d}_{j}(\varphi,\mu)=(0,0)^{T},~{}j=4,5,\cdots.$
(2.26)

where

$\left\\{\begin{array}[]{l}F^{d}_{20}(\varphi)=-2d_{21}\tau_{c}\left(\begin{array}[]{c}0\vspace{0.2cm}\\\
\varphi^{(1)}_{x}(-1)\varphi^{(2)}_{x}(0)+\varphi^{(1)}_{xx}(-1)\varphi^{(2)}(0)\end{array}\right),\vspace{0.2cm}\\\
F^{d}_{21}(\varphi)=2D_{1}\varphi_{xx}(0)+2D_{2}\varphi_{xx}(-1),\vspace{0.2cm}\\\
F^{d}_{31}(\varphi)=-6d_{21}\left(\begin{array}[]{c}0\vspace{0.2cm}\\\
\varphi^{(1)}_{x}(-1)\varphi^{(2)}_{x}(0)+\varphi^{(1)}_{xx}(-1)\varphi^{(2)}(0)\end{array}\right),\vspace{0.2cm}\\\
\end{array}\right.$ (2.27)

It is easy to verify that

$\left(\begin{array}[]{c}\left[2\mu
A\left(\Phi(0)z_{x}\right),\beta_{n_{c}}^{(1)}\right]\vspace{0.3cm}\\\
\left[2\mu
A\left(\Phi(0)z_{x}\right),\beta_{n_{c}}^{(2)}\right]\end{array}\right)=2\mu
A\left(\Phi(0)\left(\begin{array}[]{c}z_{1}\\\
z_{2}\end{array}\right)\right),$ (2.28)
$\begin{array}[]{lll}&&\left(\begin{array}[]{c}\left[\mu
F^{d}_{21}\left(\Phi(\theta)z_{x}\right),\beta_{n_{c}}^{(1)}\right]\vspace{0.3cm}\\\
\left[\mu
F^{d}_{21}\left(\Phi(\theta)z_{x}\right),\beta_{n_{c}}^{(2)}\right]\end{array}\right)\vspace{0.3cm}\\\
&=&-2(n_{c}/\ell)^{2}\mu\left(D_{1}\left(\Phi(0)\left(\begin{array}[]{c}z_{1}\\\
z_{2}\end{array}\right)\right)+D_{2}\left(\Phi(-1)\left(\begin{array}[]{c}z_{1}\\\
z_{2}\end{array}\right)\right)\right).\end{array}$ (2.29)

In addition, from (2.8), we have for all $\mu\in\mathbb{R}$,

$F_{2}(\Phi(\theta)z_{x},\mu)=F_{2}(\Phi(\theta)z_{x},0),$

It follows from (2.20) that

$f^{1}_{2}(z,0,\mu)=\Psi(0)\left(\begin{array}[]{c}\left[\widetilde{F}_{2}\left(\Phi(\theta)z_{x},\mu\right),\beta_{n_{c}}^{(1)}\right]\vspace{0.1cm}\\\
\left[\widetilde{F}_{2}\left(\Phi(\theta)z_{x},\mu\right),\beta_{n_{c}}^{(2)}\right]\end{array}\right).$
(2.30)

This, together with (2.18), (2.25), (2.28) and (2.29), yields to

$g^{1}_{2}(z,0,\mu)=\mbox{Proj}_{\tiny{\mbox{Ker}(M_{2}^{1})}}f_{2}^{1}(z,0,\mu)=\mathcal{H}\left(B_{1}\mu
z_{1}\right),$ (2.31)

where

$B_{1}=2\psi^{T}\left(A\phi(0)-(n_{c}/\ell)^{2}\left(D_{1}\phi(0)+D_{2}\phi(-1)\right)\right)=2i\omega_{n_{c}}\psi^{T}\phi.$

####  2.2.2 Calculation of $\bf{g_{3}^{1}(z,0,\mu)}$

In this subsection, we calculate the third term $g^{1}_{3}(z,0,0)$ in terms of
(2.24). Notice that $\frac{1}{3!}\widetilde{f}_{3}^{1}$ in (2.24) is the term
of order 3 obtained after the changes of variables in previous step. Denote

$f^{(1,1)}_{2}(z,w,0)=\Psi(0)\left(\begin{array}[]{c}\left[F_{2}\left(\Phi(\theta)z_{x}+w,0\right),\beta_{n_{c}}^{(1)}\right]\vspace{0.2cm}\\\
\left[F_{2}\left(\Phi(\theta)z_{x}+w,0\right),\beta_{n_{c}}^{(2)}\right]\end{array}\right),$
(2.32)
$f^{(1,2)}_{2}(z,w,0)=\Psi(0)\left(\begin{array}[]{c}\left[F^{d}_{2}\left(\Phi(\theta)z_{x}+w,0\right),\beta_{n_{c}}^{(1)}\right]\vspace{0.2cm}\\\
\left[F^{d}_{2}\left(\Phi(\theta)z_{x}+w,0\right),\beta_{n_{c}}^{(2)}\right]\end{array}\right).$
(2.33)

In addition, it follows from (2.31) that $g_{2}^{1}(z,0,0)=(0,0)^{T}$. Then
$\widetilde{f}_{3}^{1}(z,0,0)$ is determined by

$\begin{array}[]{lll}&&\widetilde{f}_{3}^{1}(z,0,0)\vspace{0.3cm}\\\
&=&f_{3}^{1}(z,0,0)+\frac{3}{2}\left[\left(D_{z}f_{2}^{1}(z,0,0)\right)U_{2}^{1}(z,0)+\left(D_{w}f_{2}^{(1,1)}(z,0,0)\right)U_{2}^{2}(z,0)(\theta)\right.\vspace{0.3cm}\\\
&&\left.+\left(D_{w,w_{x},w_{xx}}f^{(1,2)}_{2}(z,0,0)\right)U_{2}^{(2,d)}(z,0)(\theta)\right],\end{array}$
(2.34)

where $f_{2}^{1}(z,0,0)=f_{2}^{(1,1)}(z,0,0)+f_{2}^{(1,2)}(z,0,0)$,

$D_{w,w_{x},w_{xx}}f^{(1,2)}_{2}(z,0,0)=\left(D_{w}f^{(1,2)}_{2}(z,0,0),D_{w_{x}}f^{(1,2)}_{2}(z,0,0),D_{w_{xx}}f^{(1,2)}_{2}(z,0,0)\right),$
$U_{2}^{1}(z,0)=\left(M_{2}^{1}\right)^{-1}\mbox{Proj}_{Im\left(M_{2}^{1}\right)}f_{2}^{1}(z,0,0),~{}~{}U_{2}^{2}(z,0)(\theta)=\left(M_{2}^{2}\right)^{-1}f_{2}^{2}(z,0,0),$
(2.35)

and

$U_{2}^{(2,d)}(z,0)(\theta)=\mbox{col}\left(U_{2}^{2}(z,0)(\theta),U_{2x}^{2}(z,0)(\theta),U_{2xx}^{2}(z,0)(\theta)\right).$
(2.36)

Next, we compute $\mbox{Proj}_{S}\widetilde{f}_{3}^{1}(z,0,0)$ step by step
according to (2.34). The calculation is divided into the following four steps.

Step 1: The calculation of $\bf{\mbox{Proj}_{S}f_{3}^{1}(z,0,0)}$

Let

$F_{3}(\Phi(\theta)z_{x},0)=\gamma_{n_{c}}^{3}(x)\left(\sum_{q_{1}+q_{2}=3}A_{q_{1}q_{2}}z_{1}^{q_{1}}z_{2}^{q_{2}}\right),~{}q_{1},q_{2}\in\mathbb{N}_{0}.$
(2.37)

From (2.19) and (2.26), we have
$\widetilde{F}_{3}(\Phi(\theta)z_{x},0)=F_{3}(\Phi(\theta)z_{x},0)$. Then it
follows from (2.20) and (2.37) that

$f_{3}^{1}(z,0,0)=\Psi(0)\left(\sum\limits_{q_{1}+q_{2}=3}A_{q_{1}q_{2}}\int_{0}^{\ell\pi}\gamma_{n_{c}}^{q_{1}+q_{2}+1}(x)dxz_{1}^{q_{1}}z_{2}^{q_{2}}\right),$

which, together with the fact that
$\int_{0}^{\ell\pi}\gamma_{n_{c}}^{4}(x)dx=\frac{3}{2\ell\pi},$ implies that

$\mbox{Proj}_{S}f_{3}^{1}(z,0,0)=\mathcal{H}\left(B_{21}z_{1}^{2}z_{2}\right)$
(2.38)

where

$B_{21}=\frac{3}{2\ell\pi}\psi^{T}A_{21}.$ (2.39)

Step 2: The calculation of
$\bf{\mbox{Proj}_{S}\left(\left(D_{z}f_{2}^{1}\right)(z,0,0)U_{2}^{1}(z,0)\right)}$

Form (2.18) and (2.25), we have

$\widetilde{F}_{2}(\Phi(\theta)z_{x},0)=F_{2}(\Phi(\theta)z_{x},0)+F^{d}_{20}\left(\Phi(\theta)z_{x}\right)$
(2.40)

By (2.8), we write

$\begin{array}[]{lll}&&F_{2}(\Phi(\theta)z_{x}+w,\mu)=F_{2}(\Phi(\theta)z_{x}+w,0)\vspace{0.2cm}\\\
&=&\gamma_{n_{c}}^{2}(x)\left(\sum\limits_{q_{1}+q_{2}=2}A_{q_{1}q_{2}}z_{1}^{q_{1}}z_{2}^{q_{2}}\right)+\mathcal{S}_{2}(\Phi(\theta)z_{x},w)+O\left(|w|^{2}\right),~{}q_{1},q_{2}\in\mathbb{N}_{0},\end{array}$
(2.41)

where $\mathcal{S}_{2}(\Phi(\theta)z_{x},w)$ is the second cross terms of
$\Phi(\theta)z_{x}$ and $w$. In addition, by (2.25) and (2.27), we write

$\begin{array}[]{lll}&&F^{d}_{2}\left(\Phi(\theta)z_{x},0\right)=F^{d}_{20}\left(\Phi(\theta)z_{x}\right)\vspace{0.2cm}\\\
&=&(n_{c}/\ell)^{2}\left(\xi^{2}_{n_{c}}(x)-\gamma_{n_{c}}^{2}(x)\right)\left(\sum\limits_{q_{1}+q_{2}=2}A^{d}_{q_{1}q_{2}}z_{1}^{q_{1}}z_{2}^{q_{2}}\right),\end{array}$
(2.42)

where

$\xi_{n_{c}}(x)=\frac{\sqrt{2}}{\sqrt{\ell\pi}}\sin\left(\frac{n_{c}x}{\ell}\right),$
$\left\\{\begin{array}[]{l}A_{20}^{d}=-2d_{21}\tau_{c}\left(\begin{array}[]{c}0\vspace{0.2cm}\\\
\phi_{1}(-1)\phi_{2}(0)\end{array}\right)=\overline{A_{02}^{d}},\vspace{0.3cm}\\\
A_{11}^{d}=-2d_{21}\tau_{c}\left(\begin{array}[]{c}0\vspace{0.2cm}\\\
2\mbox{Re}\left\\{\phi_{1}(-1)\overline{\phi_{2}}(0)\right\\}\end{array}\right).\end{array}\right.$
(2.43)

It is easy to verify that

$\int_{0}^{\ell\pi}\gamma_{n_{c}}^{3}(x)dx=\int_{0}^{\ell\pi}\xi_{n_{2}}^{2}(x)\gamma_{n_{c}}(x)dx=0,$

Then, by (2.40)-(2.42), we have

$f^{1}_{2}(z,0,0)=\Psi(0)\left(\begin{array}[]{c}\left[\widetilde{F}_{2}\left(\Phi(\theta)z_{x},0\right),\beta_{n_{c}}^{(1)}\right]\vspace{0.2cm}\\\
\left[\widetilde{F}_{2}\left(\Phi(\theta)z_{x},0\right),\beta_{n_{c}}^{(2)}\right]\end{array}\right)=\left(\begin{array}[]{c}\
0\\\ 0\end{array}\right)$

Hence,

$\mbox{Proj}_{S}\left[\left(D_{z}f_{2}^{1}\right)(z,0,0)U_{2}^{1}(z,0)\right]=\left(\begin{array}[]{c}\
0\\\ 0\end{array}\right).$ (2.44)

Step 3: The calculation of
$\bf{\mbox{Proj}_{S}\left(\left(D_{w}f_{2}^{(1,1)}(z,0,0)\right)U_{2}^{2}(z,0)(\theta)\right)}$

Let

$U_{2}^{2}(z,0)(\theta)\triangleq
h(\theta,z)=\sum\limits_{n\in\mathbb{N}_{0}}h_{n}(\theta,z)\gamma_{n}(x),$
(2.45)

where

$h_{n}(\theta,z)=\sum\limits_{q_{1}+q_{2}=2}h_{n,q_{1}q_{2}}(\theta)z_{1}^{q_{1}}z_{2}^{q_{2}}.$

Then, from (2.32) and (2.45), we obtain

$\begin{array}[]{lll}&&\left(D_{w}f^{(1,1)}_{2}(z,0,0)\right)U_{2}^{2}(z,0)(\theta)\vspace{0.2cm}\\\
&=&\Psi(0)\left(\begin{array}[]{c}\left[\left.D_{w}F_{2}\left(\Phi(\theta)z_{x}+w,0\right)\right|_{w=0}\left(\sum\limits_{n\in\mathbb{N}_{0}}h_{n}(\theta,z)\gamma_{n}(x)\right),\beta_{n_{c}}^{(1)}\right]\vspace{0.2cm}\\\
\left[\left.D_{w}F_{2}\left(\Phi(\theta)z_{x}+w,0\right)\right|_{w=0}\left(\sum\limits_{n\in\mathbb{N}_{0}}h_{n}(\theta,z)\gamma_{n}(x)\right),\beta_{n_{c}}^{(2)}\right]\end{array}\right).\end{array}$

By (2.41), we obtain

$\left.D_{w}F_{2}\left(\Phi(\theta)z_{x}+w,0\right)\right|_{w=0}\left(\sum\limits_{n\in\mathbb{N}_{0}}h_{n}(\theta,z)\gamma_{n}(x)\right)=\mathcal{S}_{2}\left(\Phi(\theta)z_{x},\sum\limits_{n\in\mathbb{N}_{0}}h_{n}(\theta,z)\gamma_{n}(x)\right)$

and

$\begin{array}[]{lll}&&\left(\begin{array}[]{c}\left[\mathcal{S}_{2}\left(\Phi(\theta)z_{x},\sum\limits_{n\in\mathbb{N}_{0}}h_{n}(\theta,z)\gamma_{n}(x)\right),\beta_{n_{c}}^{(1)}\right]\vspace{0.2cm}\\\
\left[\mathcal{S}_{2}\left(\Phi(\theta)z_{x},\sum\limits_{n\in\mathbb{N}_{0}}h_{n}(\theta,z)\gamma_{n}(x)\right),\beta_{n_{c}}^{(2)}\right]\end{array}\right)\vspace{0.2cm}\\\
&=&\sum\limits_{n\in\mathbb{N}_{0}}b_{n}\left(\mathcal{S}_{2}\left(\phi(\theta)z_{1},h_{n}(\theta,z)\right)+\mathcal{S}_{2}\left(\overline{\phi}(\theta)z_{2},h_{n}(\theta,z)\right)\right),\end{array}$

where

$b_{n}=\int_{0}^{\ell\pi}\gamma^{2}_{n_{c}}(x)\gamma_{n}(x)dx=\left\\{\begin{array}[]{ll}\frac{1}{\sqrt{\ell\pi}},&n=0,\vspace{0.2cm}\\\
\frac{1}{\sqrt{2\ell\pi}},&n=2n_{c},\vspace{0.2cm}\\\
0,&\mbox{otherwise}.\end{array}\right.$ (2.46)

Hence,

$\begin{array}[]{lll}&&\left(D_{w}f_{2}^{(1,1)}(z,0,0)\right)U^{2}_{2}(z,0)(\theta)\vspace{0.2cm}\\\
&=&\Psi(0)\left(\begin{array}[]{l}\sum\limits_{n=0,2n_{c}}b_{n}\left(\mathcal{S}_{2}\left(\phi(\theta)z_{1},h_{n}(\theta,z)\right)+\mathcal{S}_{2}\left(\overline{\phi}(\theta)z_{2},h_{n}(\theta,z)\right)\right)\end{array}\right).\end{array}$

Then, we have

$\textmd{Proj}_{S}\left(D_{w}f_{2}^{(1,1)}(z,0,0)U^{2}_{2}(z,0)(\theta)\right)=\mathcal{H}\left(B_{22}z_{1}^{2}z_{2}\right),$
(2.47)

where

$\begin{array}[]{lll}B_{22}&=&\frac{1}{\sqrt{\ell\pi}}\psi^{T}\left(\mathcal{S}_{2}(\phi(\theta),h_{0,11}(\theta))+\mathcal{S}_{2}(\overline{\phi}(\theta),h_{0,20}(\theta))\right)\vspace{0.2cm}\\\
&&+\frac{1}{\sqrt{2\ell\pi}}\psi^{T}\left(\mathcal{S}_{2}(\phi(\theta),h_{2n_{c},11}(\theta))+\mathcal{S}_{2}(\overline{\phi}(\theta),h_{2n_{c},20}(\theta))\right).\end{array}$
(2.48)

Step 4: The calculation of
$\bf{\mbox{Proj}_{S}\left(\left(D_{w,w_{x},w_{xx}}f^{(1,2)}_{2}(z,0,0)\right)U_{2}^{(2,d)}(z,0)(\theta)\right)}$

Denote
$\varphi(\theta)=\left(\varphi^{(1)},\varphi^{(2)}\right)^{T}=\Phi(\theta)z_{x}$
and

$\begin{array}[]{lll}&&F_{2}^{d}(\varphi(\theta),w,w_{x},w_{xx})=F_{2}^{d}(\varphi(\theta)+w,0)=F_{20}^{d}(\varphi(\theta)+w)\vspace{0.3cm}\\\
&=&-2d_{21}\tau_{c}\left(\begin{array}[]{c}0\vspace{0.2cm}\\\
\left(\varphi^{(1)}_{xx}(-1)+w_{xx}^{(1)}(-1)\right)\left(\varphi^{(2)}(0)+w^{(2)}(0)\right)\end{array}\right)\vspace{0.3cm}\\\
&&-2d_{21}\tau_{c}\left(\begin{array}[]{c}0\vspace{0.2cm}\\\
\left(\varphi^{(1)}_{x}(-1)+w_{x}^{(1)}(-1)\right)\left(\varphi^{(2)}_{x}(0)+w_{x}^{(2)}(0)\right)\end{array}\right),\end{array}$

and for
$\phi(\theta)=\left(\phi_{1}(\theta),\phi_{2}(\theta)\right)^{T},y(\theta)=\left(y_{1}(\theta),y_{2}(\theta)\right)^{T}\in
C\left([-1,0],\mathbb{R}^{2}\right)$,

$\left\\{\begin{array}[]{l}\mathcal{S}_{2}^{(d,1)}(\phi(\theta),y(\theta))=-2d_{21}\tau_{c}\left(\begin{array}[]{c}0\vspace{0.2cm}\\\
\phi_{1}(-1)y_{2}(0)\end{array}\right),\vspace{0.3cm}\\\
\mathcal{S}_{2}^{(d,2)}(\phi(\theta),y(\theta))=-2d_{21}\tau_{c}\left(\begin{array}[]{c}0\vspace{0.2cm}\\\
\phi_{1}(-1)y_{2}(0)\end{array}\right)-2d_{21}\tau_{c}\left(\begin{array}[]{c}0\vspace{0.2cm}\\\
\phi_{2}(0)y_{1}(-1)\end{array}\right),\vspace{0.3cm}\par\\\
\mathcal{S}_{2}^{(d,3)}(\phi(\theta),y(\theta))=-2d_{21}\tau_{c}\left(\begin{array}[]{c}0\vspace{0.2cm}\\\
\phi_{2}(0)y_{1}(-1)\end{array}\right).\end{array}\right.$ (2.49)

Then, from (2.33), (2.36) and (2.45), we have

$\begin{array}[]{lll}&&\left(D_{w,w_{x},w_{xx}}f^{(1,2)}_{2}(z,0,0)\right)U_{2}^{(2,d)}(z,0)(\theta)\vspace{0.3cm}\\\
&=&\Psi(0)\left(\begin{array}[]{c}\left[D_{w,w_{x},w_{xx}}F_{2}^{d}(\varphi(\theta),w,w_{x},w_{xx})U_{2}^{(2,d)}(z,0)(\theta),\beta_{n_{c}}^{(1)}\right]\vspace{0.2cm}\\\
\left[D_{w,w_{x},w_{xx}}F_{2}^{d}(\varphi(\theta),w,w_{x},w_{xx})U_{2}^{(2,d)}(z,0)(\theta),\beta_{n_{c}}^{(2)}\right]\end{array}\right)\end{array}$

and then we obtain

$\textmd{Proj}_{S}\left(\left(D_{w,w_{x},w_{xx}}f^{(1,2)}_{2}(z,0,0)\right)U_{2}^{(2,d)}(z,0)(\theta)\right)=\mathcal{H}\left(B_{23}z_{1}^{2}z_{2}\right),$
(2.50)

where

$\begin{array}[]{lll}B_{23}&=&-\frac{1}{\sqrt{\ell\pi}}\left(n_{c}/\ell\right)^{2}\psi^{T}\left(\mathcal{S}_{2}^{(d,1)}\left(\phi(\theta),h_{0,11}(\theta)\right)+\mathcal{S}_{2}^{(d,1)}\left(\overline{\phi}(\theta),h_{0,20}(\theta)\right)\right)\vspace{0.2cm}\\\
&&+\frac{1}{\sqrt{2\ell\pi}}\psi^{T}\sum\limits_{j=1,2,3}b_{2n_{c}}^{(j)}\left(\mathcal{S}_{2}^{(d,j)}\left(\phi(\theta),h_{2n_{c},11}(\theta)\right)+\mathcal{S}_{2}^{(d,j)}\left(\overline{\phi}(\theta),h_{2n_{c},20}(\theta)\right)\right)\end{array}$
(2.51)

with

$b_{2n_{c}}^{(1)}=-(n_{c}/\ell)^{2},~{}b_{2n_{c}}^{(2)}=2(n_{c}/\ell)^{2},~{}b_{2n_{c}}^{(3)}=-(2n_{c}/\ell)^{2}.$

####  2.2.3 Normal form of the Hopf bifurcation and the corresponding
coefficients

According to the above calculations, we have obtained the normal form of the
Hopf bifurcation in the following form

$\dot{z}=Bz+\frac{1}{2}\left(\begin{array}[]{c}B_{1}z_{1}\mu\\\
\bar{B}_{1}z_{2}\mu\\\
\end{array}\right)+\frac{1}{3!}\left(\begin{array}[]{c}B_{2}z_{1}^{2}z_{2}\\\
\bar{B}_{2}z_{1}z^{2}_{2}\\\ \end{array}\right)+O(|z|\mu^{2}+|z^{4}|),$ (2.52)

where

$B_{1}=2i\omega_{n_{c}}\psi^{T}\phi,~{}B_{2}=B_{21}+\frac{3}{2}\left(B_{22}+B_{23}\right),$

and $B_{2j},j=1,2,3,$ are determined by (2.39), (2.48) and (2.51). Through the
change of variables $z_{1}=w_{1}-iw_{2},z_{2}=w_{1}+iw_{2}$ and
$w_{1}=\rho\cos\xi$, $w_{2}=\rho\sin\xi$, the normal form (2.52) can be
written as the following form in polar coordinates

$\begin{array}[]{rcl}\dot{\rho}&=&K_{1}\mu\rho+K_{2}\rho^{3}+O(\mu^{2}\rho+|(\rho,\mu)|^{4}),\\\
\end{array}$

with $K_{1}=\frac{1}{2}\textrm{Re}\left(B_{1}\right)$,
$K_{2}=\frac{1}{3!}\textrm{Re}\left(B_{2}\right)$. According to [4], the sign
of $K_{1}K_{2}$ determines the direction of the bifurcation (supercritical for
$K_{1}K_{2}<0$ and subcritical for $K_{1}K_{2}>0$), and the sign of $K_{2}$
determines the stability of the Hopf bifurcating periodic orbits (stable for
$K_{2}<0$ and unstable for $K_{2}>0$).

From (2.39), (2.48) and (2.51), it is obvious that in order to obtain the
value of $K_{2}$, we still need to compute
$h_{0,20}(\theta),h_{0,11}(\theta),h_{2n_{c},20}(\theta),h_{2n_{c},11}(\theta)$,
and $A_{ij}$.

From [6] , we have

$\begin{array}[]{lll}&&M_{2}^{2}\left(h_{n}(\theta,z)\gamma_{n}(x)\right)\vspace{0.2cm}\\\
&=&D_{z}\left(h_{n}(\theta,z)\gamma_{n}(x)\right)Bz-
A_{\mathcal{Q}^{1}}\left(h_{n}(\theta,z)\gamma_{n}(x)\right),\end{array}$

which leads to

$\begin{array}[]{lll}&&\left(\begin{array}[]{c}\left[M_{2}^{2}\left(h_{n}(\theta,z)\gamma_{n}(x)\right),\beta_{n}^{(1)}\right]\vspace{0.2cm}\\\
\left[M_{2}^{2}\left(h_{n}(\theta,z)\gamma_{n}(x)\right),\beta_{n}^{(2)}\right]\end{array}\right)\vspace{0.2cm}\\\
&=&2i\omega_{c}\left(h_{n,20}(\theta)z_{1}^{2}-h_{n,02}(\theta)z_{2}^{2}\right)\vspace{0.2cm}\\\
&&-\left(\dot{h}_{n}(\theta,z)+X_{0}(\theta)\left(\mathscr{L}_{0}\left(h_{n}(\theta,z)\right)-\dot{h}_{n}(0,z)\right)\right),\end{array}$
(2.53)

where

$\mathscr{L}_{0}\left(h_{n}(\theta,z)\right)=-\tau_{c}(n/\ell)^{2}\left(D_{1}h_{n}(0,z)+D_{2}h_{n}(-1,z)\right)+\tau_{c}Ah_{n}(0,z).$

By (2.21), we get

$\begin{array}[]{lll}&&f_{2}^{2}(z,0,0)\vspace{0.2cm}\\\
&=&X_{0}(\theta)\widetilde{F}_{2}\left(\Phi(\theta)z_{x},0\right)-\pi\left(X_{0}(\theta)\widetilde{F}_{2}\left(\Phi(\theta)z_{x},0\right)\right)\vspace{0.2cm}\\\
&=&X_{0}(\theta)\widetilde{F}_{2}\left(\Phi(\theta)z_{x},0\right)-\Phi(\theta)\Psi(0)\left(\begin{array}[]{c}\left[\widetilde{F}_{2}\left(\Phi(\theta)z_{x},0\right),\beta_{n_{c}}^{(1)}\right]\vspace{0.2cm}\\\
\left[\widetilde{F}_{2}\left(\Phi(\theta)z_{x},0\right),\beta_{n_{c}}^{(2)}\right]\end{array}\right)\gamma_{n_{c}}(x).\end{array}$
(2.54)

By (2.16), (2.40), (2.41) and (2.42), we have

$\begin{array}[]{lll}&&\left(\begin{array}[]{c}\left[f_{2}^{2}(z,0,0),\beta_{n}^{(1)}\right]\vspace{0.1cm}\\\
\left[f_{2}^{2}(z,0,0),\beta_{n}^{(2)}\right]\end{array}\right)\vspace{0.2cm}\\\
&=&\begin{cases}\begin{array}[]{l}\frac{1}{\sqrt{\ell\pi}}X_{0}(\theta)\left(A_{20}z_{1}^{2}+A_{02}z_{2}^{2}+A_{11}z_{1}z_{2}\right),\end{array}&n=0,\vspace{0.3cm}\\\
\frac{1}{\sqrt{2\ell\pi}}X_{0}(\theta)\left(\widetilde{A}_{20}z_{1}^{2}+\widetilde{A}_{02}z_{2}^{2}+\widetilde{A}_{11}z_{1}z_{2}\right),&n=2n_{c},\vspace{0.3cm}\\\
\end{cases}\end{array}$ (2.55)

where $\widetilde{A}_{j_{1}j_{2}}$ is defined as follows

$\left\\{\begin{array}[]{l}\widetilde{A}_{j_{1}j_{2}}=A_{j_{1}j_{2}}-2(n_{c}/\ell)^{2}A_{j_{1}j_{2}}^{d},\vspace{0.2cm}\\\
j_{1},j_{2}=0,1,2,\quad j_{1}+j_{2}=2,\end{array}\right.$

where $A_{j_{1}j_{2}}^{d}$ is determined by (2.43) and $A_{j_{1}j_{2}}$ will
be determined in the following.

Hence, from (2.35), (2.53)-(2.55) and matching the coefficients of
$z_{1}^{2},z_{1}z_{2}$, we have

$n=0,~{}\left\\{\begin{array}[]{ll}z_{1}^{2}:&\begin{cases}\dot{h}_{0,20}(\theta)-2i\omega_{c}h_{0,20}(\theta)=(0,0)^{T},\vspace{0.1cm}\\\
\dot{h}_{0,20}(0)-L_{0}(h_{0,20}(\theta))=\frac{1}{\sqrt{\ell\pi}}A_{20},\end{cases}\vspace{0.2cm}\\\
z_{1}z_{2}:&\begin{cases}\dot{h}_{0,11}(\theta)=(0,0)^{T},\vspace{0.1cm}\\\
\dot{h}_{0,11}(0)-L_{0}(h_{0,11}(\theta))=\frac{1}{\sqrt{\ell\pi}}A_{11},\end{cases}\end{array}\right.$
(2.56)
$n=2n_{c},~{}\left\\{\begin{array}[]{ll}z_{1}^{2}:&\begin{cases}\dot{h}_{2n_{c},20}(\theta)-2i\omega_{c}h_{2n_{c},20}(\theta)=(0,0)^{T},\vspace{0.1cm}\\\
\dot{h}_{2n_{c},20}(0)-\mathscr{L}_{0}(h_{2n_{c},20}(\theta))=\frac{1}{\sqrt{2\ell\pi}}\widetilde{A}_{20},\end{cases}\vspace{0.2cm}\\\
z_{1}z_{2}:&\begin{cases}\dot{h}_{2n_{c},11}(\theta)=(0,0)^{T},\vspace{0.1cm}\\\
\dot{h}_{2n_{c},11}(0)-\mathscr{L}_{0}(h_{2n_{c},11}(\theta))=\frac{1}{\sqrt{2\ell\pi}}\widetilde{A}_{11},\end{cases}\end{array}\right.$
(2.57)

Solving (2.56), (2.57), we obtain

$\left\\{\begin{array}[]{l}h_{0,20}(\theta)=\frac{1}{\sqrt{\ell\pi}}\left(\widetilde{\mathcal{M}}_{0}(2i\omega_{c})\right)^{-1}A_{20}e^{2i\omega_{c}\theta},\vspace{0.2cm}\\\
h_{0,11}(\theta)=\frac{1}{\sqrt{\ell\pi}}\left(\widetilde{\mathcal{M}}_{0}(0)\right)^{-1}A_{11},\end{array}\right.$

and

$\left\\{\begin{array}[]{lll}h_{2n_{c},20}(\theta)&=&\frac{1}{\sqrt{2\ell\pi}}\left(\widetilde{\mathcal{M}}_{2n_{c}}(2i\omega_{c})\right)^{-1}\widetilde{A}_{20}e^{2i\omega_{c}\theta}\vspace{0.2cm}\\\
h_{2n_{c},11}(\theta)&=&\frac{1}{\sqrt{2\ell\pi}}\left(\widetilde{\mathcal{M}}_{2n_{c}}(0)\right)^{-1}\widetilde{A}_{11},\end{array}\right.$

where the matrix $\widetilde{\mathcal{M}}_{n}(\lambda)$ is defined by (2.14).

######  Remark 2.1.

Although the normal form of the delay-induced Hopf bifurcation has extensively
studied, the existing procedures can not be applied to the memory-based
diffusion system because of the nonlinearity of the diffusion terms and the
existence of the delay in the diffusion terms. Compared with the existing
procedures, the procedure in our Step 4 is the main characteristic difference.

##  3 Application to a predator-prey model with Holling type-II functional
response

We now apply our newly developed algorithm in Section 2 to a a predator-prey
model, where in the model (1.1)

$f(u,v)=u\left(1-\frac{u}{a}\right)-\frac{buv}{1+u},~{}~{}g(u,v)=\frac{buv}{1+u}-cv.$

Then, (1.1) becomes the following predator-prey model with Holling type II
functional response:

$\left\\{\begin{array}[]{ll}\frac{\partial u(x,t)}{\partial
t}=d_{11}u_{xx}(x,t)+u\left(1-\frac{u}{a}\right)-\frac{buv}{1+u},&0<x<\ell\pi,t>0,\vspace{0.2cm}\\\
\frac{\partial v(x,t)}{\partial
t}=-d_{21}(v(x,t)u_{x}(x,t-\tau))_{x}+d_{22}v_{xx}(x,t)-cv+\frac{buv}{1+u},&0<x<\ell\pi,t>0,\vspace{0.2cm}\\\
u_{x}(0,t)=u_{x}(\ell\pi,t)=v_{x}(0,t)=v_{x}(\ell\pi,t)=0,&t\geq
0.\end{array}\right.$ (3.1)

###  3.1 Stability and bifurcation analysis

System (3.1) has the positive constant steady state $E_{*}(u_{*},v_{*})$,
where

$u_{*}=\frac{c}{b-c},~{}~{}v_{*}=\frac{(a-u_{*})(1+u_{*})}{ab},$

provided that $b>\frac{c(1+a)}{a}$ holds. For $E_{*}(u_{*},v_{*})$, we have

$\begin{array}[]{c}a_{11}=\frac{\gamma(a-1-2\gamma)}{a(1+\gamma)}\left\\{\begin{array}[]{ll}\leq
0,&\frac{a-1}{2}\leq\gamma<a,\vspace{0.2cm}\\\
>0,&0<\gamma<\frac{a-1}{2},\end{array}\right.\vspace{0.2cm}\\\
a_{12}=-c<0,~{}~{}a_{21}=\frac{a-\gamma}{a(1+\gamma)}>0,~{}~{}a_{22}=0,\end{array}$
(3.2)

where $\gamma=\frac{c}{b-c}$. Let

$J_{n}=d_{11}d_{22}(n/\ell)^{4}-\left(d_{11}a_{22}+d_{22}a_{11}\right)(n/\ell)^{2}+Det(A),$
(3.3)

From (3.2) and (3.3), it is easy to verify that $T_{n}<0$ and $J_{n}>0$
provided that

$(C_{0})\quad\quad\frac{a-1}{2}<\gamma<a.$

This implies that when $d_{21}=0$ and the condition $(C_{0})$ holds, the
positive steady state $E_{*}$ is asymptotically stable for $d_{11}\geq 0$ and
$d_{22}\geq 0$ In what follows, we always assume that the condition $(C_{0})$
holds.

Since $J_{n}>0$ under the condition $(C_{0})$, we have
$\Gamma(0)=J_{n}-d_{21}v_{*}a_{12}>0$. This implies that $\lambda=0$ is not a
root of Eq.(2.5).

Let $\lambda=i\omega~{}(\omega>0)$ be a root of (2.5). Substituting it along
with expressions in (3.2) into (2.5) and separating the real from the
imaginary parts, we have

$\begin{cases}\omega^{2}-J_{n}=(n/\ell)^{2}cd_{21}v_{*}\cos\left(\omega\tau\right),\vspace{0.1cm}\\\
-T_{n}\omega=(n/\ell)^{2}cd_{21}v_{*}\sin\left(\omega\tau\right),\end{cases}$
(3.4)

which yields

$\omega^{4}+P_{n}\omega^{2}+Q_{n}=0,$ (3.5)

where

$\begin{array}[]{lll}P_{n}&=&T_{n}^{2}-2J_{n}\vspace{0.2cm}\\\
&=&\left(d_{11}^{2}+d_{22}^{2}\right)(n/\ell)^{4}-2\left(d_{11}a_{11}+d_{22}a_{22}\right)(n/\ell)^{2}+a_{11}^{2}+a_{22}^{2}+2a_{12}a_{21}.\end{array}$
(3.6)

and

$Q_{n}=\left(J_{n}+cd_{21}v_{*}(n/\ell)^{2}\right)\left(J_{n}-cd_{21}v_{*}(n/\ell)^{2}\right).$
(3.7)

and

$a_{11}^{2}+a_{22}^{2}+2a_{12}a_{21}\left\\{\begin{array}[]{ll}\leq 0,&c\geq
c_{*},\vspace{0.2cm}\\\ >0,&c<c_{*},\end{array}\right.$

where

$c_{*}=\frac{\gamma^{2}(a-1-2\gamma)^{2}}{2a(1+\gamma)(a-\gamma)}.$ (3.8)

For simplification, we assume that $c<c_{*}$. Then, $P_{n}>0$ for any
$n\in\mathbb{N}_{0}$. Define

$d_{21}^{(n)}=\frac{J_{n}}{cv_{*}(n/\ell)^{2}}=\frac{1}{cv_{*}}\left(d_{11}d_{22}(n/\ell)^{2}+\frac{Det(A)}{(n/\ell)^{2}}-\left(d_{11}a_{22}+d_{22}a_{11}\right)\right)>0,$
(3.9)

Then, for fixed $n$, by (3.7) we have

$Q_{n}\left\\{\begin{array}[]{ll}>0,&0<d_{21}<d_{21}^{(n)},\vspace{0.2cm}\\\
=0,&d_{21}=d_{21}^{(n)},\vspace{0.2cm}\\\
<0,&d_{21}>d_{21}^{(n)}.\end{array}\right.$ (3.10)

Thus, when $d_{21}>d_{21}^{(n)}$, Eq.(3.5) has one positive root $\omega_{n}$,
where

$\omega_{n}=\sqrt{\frac{-P_{n}\pm\sqrt{P_{n}^{2}-4Q_{n}}}{2}}.$ (3.11)

Notice that $T_{n}<0$. From (3.4), we have

$\sin\left(\omega_{n}\tau\right)=\frac{-T_{n}\omega_{n}}{(n/\ell)^{2}cd_{21}v_{*}}>0$

Thus, set

$\tau_{n,j}=\frac{1}{\omega_{n}}\left\\{\arccos\left\\{\frac{\omega_{n}^{2}-J_{n}}{d_{21}v_{*}c(n/\ell)^{2}}\right\\}+2j\pi\right\\},~{}~{}j\in\mathbb{N}_{0},~{}n\in\mathbb{N},$
(3.12)

then (2.5) has a pair of purely imaginary roots $\pm\omega_{n}i$ at
$\tau=\tau_{n,j}^{+}$. And it is easy to verify the transversality condition
satisfies
$\left.\frac{d\mathrm{Re}(\lambda(\tau))}{d\tau}\right|_{\tau=\tau_{n,j}}>0$.

Let

$d_{21}^{*}=\min\limits_{n\in\mathbb{N}}\left\\{d_{21}^{(n)}\right\\}>0.$
(3.13)

From (3.9), it is easy to verify that $d_{21}^{~{}\left(n\right)}$ is
decreasing for $n<\ell\sqrt[4]{\frac{Det(A)}{d_{11}d_{22}}}$, is increasing
for $n>\ell\sqrt[4]{\frac{Det(A)}{d_{11}d_{22}}}$ and
$d_{21}^{~{}\left(n\right)}\to\infty$ as $n\to\infty$. This implies that
$d_{21}^{*}$ exists.

For fixed $d_{21}>d_{21}^{*}$, define an index set

$U(d_{21})=\left\\{n\in\mathbb{N}:~{}d_{21}^{(n)}<d_{21}\right\\}.$

By the above discussion, we obtain the following results on the stability and
Hopf bifurcation of system (3.1).

######  Theorem 3.1.

Assume that the condition $(C_{0})$ holds and $c<c_{*}$. Then we have the
following:

(a)

when $0<d_{21}\leq d_{21}^{*}$, the positive constant steady state
$(u_{*},v_{*})$ of system (3.1) is locally asymptotically stable for any
$\tau\geq 0$;

(b)

when $d_{21}>d_{21}^{*}$, there exists a critical value $\tau_{*}(d_{21})$ of
the delay such that

  * (b1)

the positive constant steady state $(u_{*},v_{*})$ of system (3.1) is locally
asymptotically stable for $0\leq\tau<\tau_{*}(d_{21})$ and unstable for
$\tau>\tau_{*}(d_{21})$, where

$\tau_{*}(d_{21})=\min\limits_{n\in U(d_{21})}\left\\{\tau_{n,0}\right\\}.$

  * (b2)

system (3.1) undergoes mode$-n$ Hopf bifurcation at $\tau=\tau_{n,j},n\in
U(d_{21})$, and multiple Hopf bifurcations occur when these Hopf bifurcation
curves $\tau=\tau_{n,j},n\in U(d_{21})$ have interaction point on the
$d_{21}-\tau$ plane.

###  3.2 Direction and stability of the Hopf bifurcation

We now numerically investigate the bifurcation direction and stability of the
bifurcation. To this end, we set the parameters as follows

$a=1,~{}b=\frac{3}{10},~{}c=\frac{1}{10},~{}d_{11}=\frac{3}{5},~{}d_{22}=\frac{4}{5},~{}\ell=2.$

Then, we have $(u_{*},v_{*})=(1/5,5/2)$, $\gamma=1/2$,

$a_{11}=-\frac{1}{3},a_{12}=-\frac{1}{10},a_{21}=\frac{1}{3},a_{22}=0$

and $Tr(A)=-1/2,Det(A)=1/30$. It follows from (3.6) and (3.9) that

$P_{n}=\frac{1}{16}n^{4}+\frac{1}{10}n^{2}+\frac{2}{45}>0,$

and

$d_{21}^{(n)}=\frac{12n^{2}}{25}+\frac{8}{15n^{2}}+\frac{16}{15}.$ (3.14)

Notice that $P_{n}>0$ for any $n\in\mathbb{N}$, which together with (3.10),
implies that for a fixed $n$, Eq.(3.5) has no positive root for
$d_{21}<d_{21}^{(n)}$ and has only one positive root for $d_{21}\geq
d_{21}^{(n)}$. From (3.14), it is easy to verify that
$d_{21}^{(n)}<d_{21}^{(n+1)}$ for any $n\in\mathbb{N}$, and

$d_{21}^{(1)}=2.08<d_{21}^{(2)}=3.12<d_{21}^{(3)}\doteq 5.4459.$

Therefore, by (3.13), we have $d_{21}^{*}=d_{21}^{(1)}=2.08$. It follows from
(3.8) that $c_{*}=1/6$. By Theorem (3.1), we have the following stability
result.

######  Proposition 3.2.

For system (3.1) with the parameters
$a=1,~{}b=\frac{3}{10},~{}c=\frac{1}{10},~{}d_{11}=\frac{6}{10},~{}d_{22}=\frac{8}{10},~{}\ell=2$,
when $0\leq d_{21}<d_{21}^{(1)}=2.08$, the positive constant steady state
$(u_{*},v_{*})=(1/5,5/2)$ is locally asymptotically stable for any $\tau\geq
0$;

For fixed $d_{21}>d_{21}^{(1)}$, the positive constant steady state
$(u_{*},v_{*})$ is asymptotically stable for $\tau<\tau^{*}(d_{21})$ and
unstable for $\tau>\tau^{*}(d_{21})$. Fig.1 illustrates the stability region
and the Hopf bifurcation curves in the $d_{21}-\tau$ plane for $3\leq
d_{21}\leq 8$ and $0\leq\tau\leq 6$. The Hopf bifurcation curves
$\tau=\tau_{1,0}$ and $\tau=\tau_{2,0}$ intersect at the point
$P_{1}(4.1354,4.0292)$, which is the double Hopf bifurcation point. For the
point $P_{1}(4,2)$ located in the stability region, Fig.2 illustrates the
evolution of the solution of system (3.1) starting from the initial values
$u_{0}(x)=0.2+0.1\cos(x/2),~{}v_{0}(x)=2.5+0.1\cos(x/2)$, finally converging
to the constant steady state $(u_{*},v_{*})$.

From Fig.1, it is obvious to see that

$\tau^{*}(d_{21})=\left\\{\begin{array}[]{ll}\tau_{1,0},&d_{21}^{(1)}<d_{21}<4.1354,\vspace{0.2cm}\\\
\tau_{2,0},&4.1354<d_{21}<8.\end{array}\right.$

Figure 1: Stability region and Hopf bifurcation curves in $d_{21}-\tau$ plane.
The dotted region is the stability region and $\tau=\tau_{k,0},k=1,2,3$, are
Hopf bifurcation curves. Hopf bifurcation curves $\tau=\tau_{1,0}$ and
$\tau=\tau_{2,0}$ intersect at the point $P(4.1354,4.0292)$. The points
$P_{1}(4,2)$, $P_{2}(3.6,5.3)$, $P_{3}(6,2)$ and $P_{4}(4.3,5.2)$ are chosen
for the numerical simulations.

For $d_{21}=3.6$, it follows from (3.11) and (3.12) that

$\tau_{1,0}\doteq 5.1033<\tau_{2,0}\doteq 6.6493.$

For the critical mode-$1$ Hopf bifurcation at $\tau_{1,0}\doteq 5.1033$, the
direction and stability of this Hopf bifurcation can be determined by the
procedure in the previous section with $\tau_{c}=\tau_{1,0}\doteq 5.1033$ and
$d_{21}^{c}=3.6$. A directional calculation shows that

$K_{1}\doteq 0.0597>0,~{}~{}K_{2}\doteq-1.5624<0,$

which implies that the mode$-1$ spatially imhomogeneous Hopf bifurcation at
$\tau_{1,0}$ is supercritical and stable. For $\tau=5.2>\tau_{1,0}$,
Fig.3$(a)-(b)$ illustrate the existence of the spatially inhomogeneous
periodic solution with mode-$1$ spatial pattern.

Figure 2: The spatial-temporal dynamics of system (3.1) with the parameters
$a=1,~{}b=\frac{3}{10},~{}c=\frac{1}{10},~{}d_{11}=\frac{6}{10},~{}d_{22}=\frac{8}{10},~{}\ell=2$
and $\left(d_{21},\tau\right)$ being chosen as the point $P_{1}(4,2)$ of
Fig.1. The constant steady state is stable is asymptotically stable. The
initial values are $u_{0}(x)=0.2+0.1\cos(x/2),~{}v_{0}(x)=2.5+0.1\cos(x/2)$.

For $d_{21}=6$, it follows from (3.11) and (3.12) that

$\tau_{2,0}\doteq 1.8398<\tau_{1,0}\doteq 2.3542.$

For the critical Hopf bifurcation value $\tau_{2,0}\doteq 1.8398$, we have

$K_{1}\doteq 0.1733>0,~{}~{}K_{2}\doteq-2.2283<0,$

which implies that the mode$-2$ spatially imhomogeneous Hopf bifurcation at
$\tau_{2,0}$ is also supercritical and stable. For $\tau=2>\tau_{2,0}$,
Fig.3$(c)-(d)$ illustrate the existence of the spatially inhomogeneous
periodic solution with mode-$2$ spatial pattern.

(a)  (b)  (c)  (d)

Figure 3: The spatial-temporal dynamics of system (3.1) with the parameters
$a=1,~{}b=\frac{3}{10},~{}c=\frac{1}{10},~{}d_{11}=\frac{6}{10},~{}d_{22}=\frac{8}{10},~{}\ell=2$.
$(a)-(b)$ For the point $P_{2}(3.6,5.3)$ of Fig.1, there exists a spatially
inhomogeneous periodic solution with mode-1 spatial pattern. The initial
values are $u_{0}(x)=0.2+0.1\cos(x/2),~{}v_{0}(x)=2.5+0.1\cos(x/2)$. $(c)-(d)$
For the point $P_{3}(6,2)$ of Fig.1, there exists a spatially inhomogeneous
periodic solution with mode-2 spatial pattern. The initial values are
$u_{0}(x)=0.2+0.1\cos(x),~{}v_{0}(x)=2.5+0.1\cos(x)$.

We would like to mention that the interaction of mode-$1$ and mode-$2$ Hopf
bifurcations leads to more complex dynamics. To investigate the dynamical
classification near the double Hopf point $P$ in detail, the normal form for
double Hopf bifurcation should be further developed. For the point
$P_{4}(4.3,5.2)$ of Fig.1 far from the Hopf bifurcation curves
$\tau=\tau_{k,0},k=1,2$, Fig.4 shows the pattern transition from the mode-$2$
spatially inhomogeneous periodic solution (Figs.4$(a)$ and $(d)$ ) to mode-$1$
spatially inhomogeneous periodic solution (Figs.4$(c)$ and $(f)$). Figs.4$(b)$
and $(e)$ also illustrate the transiently unstable quasi-periodic patterns.

Figure 4: The spatial-temporal dynamics of system (3.1) with the parameters
$a=1,~{}b=\frac{3}{10},~{}c=\frac{1}{10},~{}d_{11}=\frac{6}{10},~{}d_{22}=\frac{8}{10},~{}\ell=2$
and $\left(d_{21},\tau\right)$ being the point $P_{4}(4.3,5.2)$ of Fig.1 far
from the Hopf bifurcation curves. Pattern transition from the mode-$2$
spatially inhomogeneous periodic solution to mode-$1$ spatially inhomogeneous
periodic solution. The initial values are
$u_{0}(x)=0.2+0.1\cos(x),~{}v_{0}(x)=2.5+0.2\cos(x)$.

##  4 Conclusion and discussion

In this paper, we have developed an algorithm for computing the normal form of
the Hopf bifurcation for the reaction-diffusion systems with memory diffusion.
Because of the nonlinearity of the diffusion term and the the presence of the
delay in the diffusion term, the traditional algorithms for computing the
normal form of the Hopf bifurcation for the reaction-diffusion system is not
applied to this system. To fill this gap, we generalized the existing
algorithms for the reaction-diffusion systems where the diffusion terms are
linear and the delay only appears in the reaction terms, to the case where the
diffusion terms are not linear and the delay appears not only in the reaction
terms but also in the diffusion terms.

As an illustration of this newly developed algorithm, we considered a
diffusive predator-prey system with memory-based diffusion and Holling type-II
functional response. The memory delay-induced spatially inhomogeneous Hopf
bifurcations and the double Hopf bifurcations due to their interactions are
observed.

We would also like to mention that the delay-induced double Hopf bifurcation
often occurs in reaction-diffusion systems with memory diffusion and may lead
to more complex dynamics like two/three invariant torus. To determine the
dynamical classifications near the double Hopf bifurcation points, new
algorithms for computing the normal form of the double Hopf bifurcation of the
reaction-diffusion systems with memory diffusion are desired.

## References

  * [1] Qi An, Chuncheng Wang, and Hao Wang.  Analysis of a spatial memory model with nonlocal maturation delay and hostile boundary condition.  Discret. Contin. Dyn. Syst., 40(10):5845–5868, OCT 2020\.
  * [2] Jianxin Chen, Tonghua Zhang, and Yongwu Zhou.  Dynamics of a risk-averse newsvendor model with continuous-time delay in supply chain financing.  Math. Comput. Simul., 169:133–148, MAR 2020.
  * [3] Mengxin Chen, Ranchao Wu, Biao Liu, and Liping Chen.  Spatiotemporal dynamics in a ratio-dependent predator-prey model with time delay near the Turing-Hopf bifurcation point.  Commun. Nonlinear Sci. Numer. Simul., 77:141–167, OCT 2019\.
  * [4] S-N Chow and J K Hale.  Methods of Bifurcation Theory.  Springer-Verlag, New York, 1982.
  * [5] Huan Dai, Yuying Liu, and Junjie Wei.  Stability analysis and Hopf bifurcation in a diffusive epidemic model with two delays.  Math. Biosci. Eng., 17(4):4127–4146, 2020.
  * [6] T Faria.  Normal forms and Hopf bifurcation for partial differential equations with delays.  Trans. Am. Math. Soc., 352(5):2217–2238, 2000.
  * [7] T Faria and L T Magalhães.  Normal forms for retarded functional differential equations with parameters and applications to Hopf singularity.  J. Differ. Equ., 122(2):181–200, NOV 1 1995.
  * [8] Jack K Hale.  Theory of Functional Differential Equations.  Springer-Verlag, Berlin, 1977.
  * [9] B.D. Hassard, N.D. Kazarinoff, and Y.H. Wan.  Theory and Applications of Hopf Bifurcation.  Cambridge University Press, Cambridge, 1981.
  * [10] Weihua Jiang, Hongbin Wang, and Xun Cao.  Turing instability and Turing-Hopf bifurcation in diffusive Schnakenberg systems with gene expression time delay.  J. Dyn. Differ. Equ., 31(4):2223–2247, DEC 2019.
  * [11] Xin Jiang, Ran Zhang, and Zhikun She.  Dynamics of a diffusive predator-prey system with ratio-dependent functional response and time delay.  Int. J. Biomath., 13(6):2050036, AUG 2020.
  * [12] Longyue Li, Yingying Mei, and Jianzhi Cao.  Hopf Bifurcation Analysis and Stability for a Ratio-Dependent Predator-Prey Diffusive System with Time Delay.  Int. J. Bifurcation Chaos, 30(3):2050037, MAR 15 2020\.
  * [13] Jiazhe Lin, Rui Xu, and Liangchen Li.  Turing-Hopf bifurcation of reaction -diffusion neural networks with leakage delay.  Commun. Nonlinear Sci. Numer. Simul., 85:105241, JUN 2020\.
  * [14] Kalyan Manna and Malay Banerjee.  Stability of Hopf-bifurcating limit cycles in a diffusion-driven prey-predator system with Allee effect and time delay.  Math. Biosci. Eng., 16(4):2411–2446, 2019.
  * [15] Xin-You Meng and Jiao-Guo Wang.  Dynamical analysis of a delayed diffusive predator-prey model with schooling behaviour and Allee effect.  J. Biol. Dyn., 14(1):826–848, JAN 1 2020.
  * [16] Karen A. Oliveira and Juliana M. Berbert.  Crossover in spreading behavior due to memory in population dynamics.  Math. Biosci., 324:108346, JUN 2020.
  * [17] Zuolin Shen and Junjie Wei.  Bifurcation Analysis in a Diffusive Mussel-Algae Model with Delay.  Int. J. Bifurcation Chaos, 29(11): 10.1142/S021812741950144X, OCT 2019.
  * [18] Junping Shi, Chuncheng Wang, and Hao Wang.  Diffusive spatial movement with memory and maturation delays.  Nonlinearity, 32(9):3188–3208, SEP 2019.
  * [19] Junping Shi, Chuncheng Wang, Hao Wang, and Xiangping Yan.  Diffusive Spatial Movement with Memory.  J. Dyn. Differ. Equ., 32(2):979–1002, JUN 2020.
  * [20] Qingyan Shi, Junping Shi, and Hao Wang.  Spatial movement with distributed delay.  J. Math. Biol., 28(4):https://doi.org/10.1007/s00285–021–01588–0, 2021.
  * [21] Y. Song, H. Jiang, and Y. Yuan.  Turing-Hopf bifurcation in the reaction-diffusion system with delay and application to a diffusive predator-prey model.  J. Appl. Anal. Comput., 9(3):1132–1164, JUN 2019.
  * [22] Yongli Song, Junping Shi, and Hao Wang.  Spatiotemporal dynamics of a diffusive consumer-resource model with explicit spatial memory.  , page submitted for review, 2021.
  * [23] Yongli Song and Qingyan Shi.  Stability and spatially non-homogeneous Hopf bifurcation in a diffusive delayed predator-prey model with spatial average.  , page submitted for review, 2020.
  * [24] Yongli Song, Shuhao Wu, and Hao Wang.  Spatiotemporal dynamics in the single population model with memory-based diffusion and nonlocal effect.  J. Differ. Equ., 267(11):6316–6351, NOV 15 2019.
  * [25] Yongli Song, Tonghua Zhang, and Yahong Peng.  Turing-Hopf bifurcation in the reaction-diffusion equations and its applications.  Commun. Nonlinear Sci. Numer. Simul., 33:229–258, APR 2016\.
  * [26] Xiaosong Tang, Heping Jiang, Zhiyun Deng, and Tao Yu.  Delay induced subcritical Hopf bifurcation in a diffusive predator-prey model with herb behavior and hyperbolic mortality .  J. Appl. Anal. Comput., 7(4):1385–1401, 2017.
  * [27] Xiaosong Tang and Yongli Song.  Stability, Hopf bifurcations and spatial patterns in a delayed diffusive predator-prey model with herd behavior.  Appl. Math. Comput., 254:375–391, MAR 1 2015.
  * [28] Shuhao Wu and Yongli Song.  Spatiotemporal dynamics of a diffusive predator-prey model with nonlocal effect and delay.  Commun. Nonlinear Sci. Numer. Simul., 89:105310, OCT 2020\.
  * [29] Ruizhi Yang, Ming Liu, and Chunrui Zhang.  A delayed-diffusive predator-prey model with a ratio-dependent functional response.  Commun. Nonlinear Sci. Numer. Simul., 53:94–110, DEC 2017\.
  * [30] Ruizhi Yang, Chunrui Zhang, and Yazhuo Zhang.  A delayed diffusive predator-prey system with Michaelis-Menten type predator harvesting.  Int. J. Bifurcation Chaos, 28(8):1850099, JUL 2018.
  * [31] Jia-Fang Zhang and Xiang-Ping Yan.  Effects of Delay and Diffusion on the Dynamics of a Leslie-Gower Type Predator-Prey Model.  Int. J. Bifurcation Chaos, 24(4):1450043, APR 2014.
  * [32] Wenjie Zuo and Junjie Wei.  Stability and Hopf bifurcation in a diffusive predator-prey system with delay effect.  Nonlinear Anal.-Real World Appl., 12(4):1998–2011, AUG 2011.
  * [33] Wenjie Zuo and Junjie Wei.  Stability and bifurcation in a ratio-dependent Holling-III system with diffusion and delay.  Nonlinear Anal.-Model Control, 19(1):132–153, 2014.