2. The influence of different kernel function

In this section, we study the influence of spatial nonlocality on the dynamic behaviour by selecting different kernel functions of the trigonometric function class. In the following, we discuss the impact of individual resource distribution on population size in several scenarios. In order to compare the dynamic behaviours between the systems (1.2) and (1.1), we fix the following same parameters as in [Reference Duan, Niu and Wei7]:

(2.1) \begin{align} d_{1}=0.3, K=21, a_0=0.28, r=0.85, h=0.9, a=0.28, b=0.998, d=0.1, l=1.2, \end{align}

and choose $d_{2}=0.8$ , the initial functions $u(x,t)=6.52+0.01\cos x$ and $v(x,t)=5.53+0.01\cos x$ .

Case I: When $W(x,y)=\frac{1}{l\pi }+\mu \cos\!(x-y)$ , $\mu \in \big [{-}\frac{1}{l\pi },\frac{1}{l\pi }\big ]$ . When $\mu =0$ , $W(x,y)=\frac{1}{l\pi }$ , implying that the intensity of competition is the same for the prey group. If $\mu \in \big [{-}\frac{1}{l\pi },0)\cup (0,\frac{1}{l\pi }\big ]$ , $W(x,y)$ is a symmetric kernel due to the form of $\cos\!(x-y)$ , the nonlocal effect is a function of the distance. Taking $\mu =0$ and $\mu =\frac{1}{l\pi }$ , respectively, we can observe the significant difference in the prey population (see Figure 1). With $\mu =0$ , there exits stable constant state for small delay ( $\tau =2$ ), then spatially nonhomogeneous periodic solutions ( $\tau =8$ ), and spatially homogeneous periodic solutions ( $\tau =12$ ) when the delay is increased gradually (Figure 1 (a)(c)(e)), while a stable state converges to a stable cosine form when $\mu =\frac{1}{l\pi }$ for small delay ( $\tau =2$ ) and the spatial ‘cos-type’ shape periodic solutions( $\tau =8$ ) and then steady-state solutions ( $\tau =12$ ) exist with the increasing of delay $\tau$ (Figure 1 (b)(d)(f)). Notice that in Figure 1, we only show the variation of prey and the predator evolves in a similar form. Figure 1(d) shows that the prey distribution is concentrated at both ends rather than in the middle. Figure 1(b)(f) illustrate that without temporal vibration, the spatial distribution of species is uneven, with one concentrated at both ends and one concentrated near $x=0$ .

Figure 1. The spatiotemporal patterns of prey with $\tau =2,8,12$ . (a,c,e): $\mu =0$ . (b,d,f): $\mu =\frac{1}{l\pi }$ .

Case II: It should be pointed out that in (1.2), substituting $W(x,y)=g_1(x)g_2(y)$ yields the partial form:

(2.2) \begin{align} \frac{1}{K/g_1(x)}\int _{\Omega }g_2(y)u(y,t)dy. \end{align}

Note that $K/g_1(x)$ , $g_2(y)$ correspond to the environmental carrying capacity (resource) at $x$ and the consumption capacity at $y$ , respectively. By changing the time delay, when $g_1(x)=1-\sin (\frac{1}{l}x)$ , $g_2(y)=\sin \frac{1}{l}y$ , we can observe spatially nonhomogeneous steady-state solutions and spatially nonhomogeneous periodic solutions (Figure 2), indicating that the variation of $g_1(x)$ in the kernel function $W(x,y)$ can directly reflect the density distribution of prey. When $g_1(x)=\sin \frac{1}{l}x$ , $g_2(y)=1-\sin (\frac{1}{l}y)$ , we can observe that the prey ultimately exhibits steady-state and periodic solutions related to its distribution position (Figure 3). From Figures 2 and 3, it can be seen that the density distribution of prey is closely related to the form of $g_1(x)$ . If we project Figures 2 and 3 onto the $x-u$ plane, the shape of the solution is mainly influenced by $[g_1(x)]^{-1}$ . More specially, as $x$ increases, the population density in Figure 2 shows a trend of first increasing and then decreasing, while in Figure 3, the opposite is true. The density of prey decreases at first and then increases.

Figure 2. The spatiotemporal diagram of prey. (a) $\tau =2$ . (b) $\tau =6$ .

Figure 3. The spatiotemporal diagram of prey. (a) $\tau =2$ . (b) $\tau =12$ .

In particular, if $g_1(x)=C$ remains constant, $g_2(y)=\sin \frac{1}{l}y$ , $W(x,y)=g_1(x)g_2(y)$ , (2.2) becomes

\begin{equation*}\frac {C}{K}\int _{\Omega }\sin \frac {y}{l}u(y,t)dy.\end{equation*}

The system (1.2) eventually exhibits a stable equilibrium solution, and the smaller the value of $C$ , the larger the value of the solution; see Figure 4. Figure 4(c) demonstrates that as $C$ increases, the steady-state solution $u_{\ast }$ gradually decreases in contrast. If we keep $g_2(y)$ at the constant, and change $g_1(x)$ as $g_1(x)=\sin \frac{1}{l}x$ , and $\tau =2$ . By changing the value of $C$ , it is found that the number of prey eventually converges to a stable nonhomogeneous steady-state solution in space, with the projection curve opening upwards, which is exactly opposite to the direction of $g_1(x)$ opening; see Figure 5.

Figure 4. The prey ultimately reaches a stable steady-state solution with $\tau =2$ . (a) $C=0.5$ . (b) $C=2.5$ . (c) $u_{\ast }$ monotonically decreases with $C$ .

Figure 5. The prey eventually converges to a stable nonhomogeneous steady-state solution. (a) $C=0.5$ . (b) $C=6.5$ .

Intuitively, through the numerical simulations, we can observe that the system exhibits spatially homogeneous or nonhomogeneous periodic solutions, constant or nonconstant steady-state solutions, etc., and the shapes of the solutions are closely related to the choice of kernel functions. Next, we analyse the dynamical properties mathematically. If the kernel function is uniformly distributed, the system may generate constant/nonconstant steady-state solutions or homogeneous/inhomogeneous periodic oscillations. The case of nonuniform distribution usually leads to the occurrence of extreme steady-state solutions in the system. This further indicates that different choices of kernel functions can lead to different dynamics in the system.

3. Existence of double Hopf bifurcation for the uniform kernel

It is straightforward to see that, same as that in the system (1.1), there is a unique constant solution $E_{\ast }=(u_{\ast },v_{\ast })$ for the system (1.2) when $d-a(b-d h)K\lt 0$ , with

\begin{equation*}v_{\ast }=\frac {r}{a}(1+ahu_{\ast })\big (1-\frac {u_{\ast }}{K}\big ).\end{equation*}

As $u_{\ast }$ gradually increases, the trend of $v_{\ast }$ is increases first and then decreases. When $u_{\ast }=\frac{ahK-1}{2ah}$ , $v_{\ast }$ reaches its maximum value. As seen from Figure 6, when the predator cannibalism factor reaches a certain point, the predator population reaches a maximum by selecting appropriate parameters, indicating that cannibalism can promote species reproduction. Cannibalism can regulate population size and benefit cannibalism individuals and their relatives, as additional shelter, territory, and food resources are released.

Figure 6. Relationship between population size and cannibalism factors.

Mathematically, let $\Omega =(0,l\pi )$ , linearising system (1.2) at $E_{\ast }$ , we obtain

\begin{eqnarray*} \frac{\partial }{\partial t}W(x,t) =D \Delta W(x,t)+A_0W(x,t)+B_0\hat{W}(x,t) +CW(x,t-\tau ), \end{eqnarray*}

where

\begin{equation*}W(x,t)=\begin {pmatrix} u(x,t) \\[6pt] v(x,t) \\[6pt] \end {pmatrix}, \hat {W}(x,t)= \begin {pmatrix} \hat {u}(x,t) \\[6pt] \hat {v}(x,t) \\[6pt] \end {pmatrix}, W(x,t-\tau )=\begin {pmatrix} u(x,t-\tau ) \\[6pt] v(x,t-\tau ) \\[6pt] \end {pmatrix},\end{equation*}

(3.1) \begin{equation} D=\begin{pmatrix} d_{1} &\quad 0 \\[6pt] 0 & \quad d_{2} \\[6pt] \end{pmatrix}, A_{0}=\begin{pmatrix} \alpha _{11} & \quad \alpha _{12} \\[6pt] \alpha _{21} &\quad \alpha _{22} \\[6pt] \end{pmatrix}, B_{0}=\begin{pmatrix} \beta _{11} &\quad 0 \\[6pt] 0 &\quad 0 \\[6pt] \end{pmatrix}, C=\begin{pmatrix} 0 &\quad 0 \\[6pt] c_{21} & 0 \\[6pt] \end{pmatrix}, \end{equation}

and

\begin{eqnarray*} &&\alpha _{11}=\frac{a^{2}h u_{\ast } v_{\ast }}{(1+a h u_{\ast })^{2}}, \alpha _{12}=-\frac{au_{\ast }}{1+a h u_{\ast }}, \alpha _{21}=\frac{a a_{0}h v_{\ast }^{2}}{(1+a h u_{\ast })^{2}},\\ &&\alpha _{22}=-\frac{a_{0} v_{\ast }}{1+a h u_{\ast }}, \beta _{11}=-\frac{ru_{\ast }}{K}, c_{21}=\frac{a b v_{\ast }}{(1+a h u_{\ast })^{2}}. \end{eqnarray*}

The characteristic problem

\begin{eqnarray*} -\Delta \varphi =\sigma \varphi,x\in (0,l\pi ), \varphi '(0)=\varphi '(l\pi )=0 \end{eqnarray*}

results the eigenvalues $\frac{n^2}{l^2}$ $(n\in \mathbb{N}_{0})$ , and the normalised eigenfunctions

(3.2) \begin{eqnarray} \begin{split} \xi _{n}(x) = \frac{\textrm{cos}\frac{n}{l}x}{\|\textrm{cos}\frac{n}{l}x\|}=\left \{ \begin{array}{ll} \sqrt{\frac{1}{l\pi }}, &n=0,\\[6pt] \sqrt{\frac{2}{l\pi }}\textrm{cos}\frac{n}{l}x, & n\geqslant 1. \end{array} \right. \end{split} \end{eqnarray}

Let $\beta _n^i(x)=\xi _{n}(x)e_i$ , $i=1,2$ , where $e_i$ is the unit coordinate vector of $\mathbb{R}^2$ . First, we define the real-valued Hilbert space

\begin{eqnarray*} X:\!=\big \{(u,v)\in H^{2}(0,l\pi )\times H^{2}(0,l\pi)\,{:}\,\partial _{x}u(0,t)=\partial _{x}u(l\pi,t)=0, \partial _{x}v(0,t)=\partial _{x}v(l\pi,t)=0\big \}, \end{eqnarray*}

and the complexification of $X$ is $X_{\mathbb{C}}:\!=\{x_{1}+ix_{2}, x_{1},x_{2}\in X \}$ . $\{\beta _n^i\}_{n\in \mathbb{N}_0}$ forms an orthonormal basis of $X_{\mathbb{C}}$ . For $U\in X_{\mathbb{C}}$ and $\beta _n=(\beta _n^1,\beta _n^2)$ , we define $\langle \beta _n,U\rangle =(\langle \beta _n^1,U\rangle,\langle \beta _n^2,U\rangle )^{T}$ , and

\begin{equation*}\mathfrak {B}_n=\textrm {span}\{\langle \beta _n^i,U\rangle \beta _n^i |U\in X_{\mathbb {C}},i=1,2\}.\end{equation*}

The sequence of characteristic equations is

(3.3) \begin{align} \lambda ^{2}+A\lambda +\tilde{B}-\alpha _{12}c_{21}e^{-\lambda \tau }=0, \end{align}

(3.4) \begin{align} \lambda ^{2}+D_{n}\lambda +E_{n}-\alpha _{12}c_{21}e^{-\lambda \tau }=0, \,n\in \mathbb{N}_{+}, \end{align}

with $A=-(\alpha _{11}+\beta _{11}+\alpha _{22})$ , $\tilde{B}=\alpha _{22}(\alpha _{11}+\beta _{11}) -\alpha _{12}\alpha _{21}$ , $D_{n}=(d_{1}+d_{2})n^{2}/l^{2}-\alpha _{11}-\alpha _{22}$ , $E_{n}= (d_{1}n^{2}/l^{2}-\alpha _{11})(d_{2}n^{2}/l^{2}-\alpha _{22})-\alpha _{12}\alpha _{21}$ . It is easy to see that

\begin{equation*}\tilde {B}-\alpha _{12}c_{21}=\frac {\big [ra_0(1+a h u_{\ast })^{2}+Ka^2b\big ]u_{\ast }v_{\ast }}{1+a h u_{\ast }}\gt 0.\end{equation*}

Thus, if

\begin{equation*}({\textbf{S}}_{\textbf{1}})\,A\gt 0, \,D_{n}\gt 0,\,E_{n}-\alpha _{12}c_{21}\gt 0,\end{equation*}

the following conclusion can be drawn.

To find the critical value of $\tau$ such that the equations (3.3) or (3.4) has a pair of simple purely imaginary roots $\pm i\omega$ $(\omega \gt 0)$ under $({\textbf{S}}_{\textbf{1}})$ , by using the similar method in [Reference Ruan and Wei19], and let $\lambda =i\omega$ be solutions of (3.3), we have

(3.5) \begin{eqnarray} \begin{cases} \displaystyle\cos \omega \tau =\frac{ -\omega ^{2}+\tilde{B}}{\alpha _{12}c_{21}}:\!=C(\omega ), \\[12pt] \displaystyle\sin \omega \tau = \frac{-A\omega }{\alpha _{12}c_{21}} :\!=S(\omega ) \gt 0. \end{cases} \end{eqnarray}

That is,

(3.6) \begin{eqnarray} \omega ^{4}+(A^{2}-2\tilde{B})\omega ^{2}+\tilde{B}^{2}-\alpha _{12}^{2}c_{21}^{2}=0. \end{eqnarray}

The two positive roots of equation (3.6) are represented as:

(3.7) \begin{eqnarray} \omega ^{\pm } =\bigg [\frac{1}{2}\Big ({-}A^{2}+2\tilde{B}\pm \sqrt{A^{4}-4A^{2}\tilde{B}+4\alpha _{12}^{2}c_{21}^{2}}\Big )\bigg ]^{1/2} \end{eqnarray}

if

(3.8) \begin{eqnarray} ({\textbf{S}}_{\textbf{2}})\,-A^{2}+2\tilde{B}\gt 0,\,\tilde{B}+\alpha _{12}c_{21}\gt 0,\,A^{4}-4A^{2}\tilde{B}+4\alpha _{12}^{2}c_{21}^{2}\gt 0 \end{eqnarray}

holds. In particular, equation (3.6) has a unique positive root $\omega ^{+}$ if

\begin{equation*}({\textbf{S}}_{\textbf{3}})\,\tilde {B}+\alpha _{12}c_{21}\lt 0.\end{equation*}

Similarly, when

(3.9) \begin{eqnarray} \omega _{n}^{\pm } =\bigg [\frac{1}{2}\Big ({-}D_{n}^{2}+2E_{n}\pm \sqrt{D_{n}^{4}-4D_{n}^{2}E_{n}+4\alpha _{12}^{2}c_{21}^{2}}\Big )\bigg ]^{1/2}, \end{eqnarray}

$i\omega _{n}^{\pm }$ are the solution in (3.4) if

(3.10) \begin{eqnarray} ({\textbf{S}}_{\textbf{4}})\,-D_{n}^{2}+2E_{n}\gt 0,\,E_{n}+\alpha _{12}c_{21}\gt 0,\,D_{n}^{4}-4D_{n}^{2}E_{n}+4\alpha _{12}^{2}c_{21}^{2}\gt 0. \end{eqnarray}

With

\begin{equation*}({\textbf{S}}_{\textbf{5}})\,E_{n}+\alpha _{12}c_{21}\lt 0,\end{equation*}

$\omega _{n}^{+}$ exists. The critical values of $\tau$ , created to $\omega ^{\pm }$ / $\omega _{n}^{\pm }$ , are

\begin{eqnarray*} \begin{cases} \displaystyle\tau ^{j\pm }=(\!\arccos C(\omega )+2j\pi )/\omega ^{\pm },\\[6pt] \displaystyle \tau _{n}^{j\pm }=(\!\arccos C_{n}(\omega )+2j\pi )/\omega _{n}^{\pm }, \,j\in \mathbb{N}_{0},\,n\in \mathbb{N}_{+}, \end{cases} \end{eqnarray*}

where $C_{n}(\omega )=({-}\omega ^{2}+E_{n})/(\alpha _{12}c_{21})$ , respectively.

Compared with the previous study, the nonlocal competition term in model (1.2) can improve the chances of coexistence between prey and predators to a certain extent. In this section, we have established the existence conditions of double Hopf bifurcation. In the next section, theoretically, we will calculate the normal form near the bifurcation point of the double Hopf, mainly by calculating the relevant coefficients, which are significantly different from the formula derived by previous researchers. Subsequently, based on the positivity and negativity of coefficient symbols, dynamic classification can be explored.

4. Normal form of double Hopf bifurcation

To study the dynamic behaviour in the system (1.2) near the double Hopf bifurcation point, we then calculate the normal form. Based on the normal form of the reaction–diffusion equation without time delay in [Reference Shen, Liu and Wei22], we give a new algorithm of the normal form for system (1.2) with time delay and nonlocal term. Let $u(t)\rightarrow u({\cdot},\tau t)-u_{\ast }$ , $v(t)\rightarrow v({\cdot},\tau t)-v_{\ast }$ , and $\hat{u}(t)\rightarrow \hat{u}({\cdot},\tau t)-u_{\ast }$ , then we obtain

(4.1) \begin{equation} \begin{cases} \displaystyle \dot{u}(t)&= \tau \big [d_{1}{\Delta } u(t)+\alpha _{11} u(t)+\alpha _{12} v(t)+ \beta _{11} \hat{u}(t)+\alpha _{1} u^{2}(t) + \tilde{\alpha }_{1}u(t)\hat{u}(t)\\[6pt] &\quad\displaystyle +\, \alpha _{2}u(t)v(t)+\alpha _{3}u^{3}(t)+\alpha _{4}u^{2}(t)v(t)+O(4) \big ], \\[6pt] \displaystyle\dot{v}(t)&=\tau \big [d_{2}{\Delta } v(t)+\alpha _{21}u(t)+\alpha _{22}v(t)+c_{21}u(t-1)+\beta _{1}u^{2}(t) + \beta _{2} u(t)v(t)\\[6pt] &\quad\displaystyle+\,\beta _{3}u^{2}(t-1)+\beta _{4}u(t-1)v(t)+\beta _{5}v^{2}(t)+\beta _{6}u^{3}(t)+\beta _{7}u^{2}(t)v(t)\\[6pt] &\quad\displaystyle+\,\beta _{8}u(t)v^{2}(t)+\beta _{9}u^{3}(t-1)+\beta _{10}u^{2}(t-1)v(t)+O(4)\big ], \end{cases} \end{equation}

where

\begin{eqnarray*} \begin{split} \displaystyle&\alpha _{1}=\frac{a^{2}h v_{\ast }}{(1+a h u_{\ast })^{3}}, \tilde{\alpha }_{1}=-\frac{r}{K}, \alpha _{2}=\frac{-a}{(1+a h u_{\ast })^{2}},\alpha _{3}=\frac{-a^{3}h^{2} v_{\ast }}{(1+a h u_{\ast })^{4}}, \alpha _{4} = \frac{a^{2}h }{(1+a h u_{\ast })^{3}},\\[6pt] &\displaystyle\beta _{1}=-\frac{a^{2}a_{0}h^{2} v_{\ast }^{2}}{(1+a h u_{\ast })^{3}}, \beta _{2} = \frac{2aa_{0}h v_{\ast }}{(1+a h u_{\ast })^{2}}, \beta _{3}=\frac{-a^{2}b h v_{\ast }}{(1+a h u_{\ast })^{3}}, \beta _{4} = \frac{ab}{(1+a h u_{\ast })^{2}}, \beta _{5}=\frac{-a_{0}}{1+a h u_{\ast }},\\[6pt] &\displaystyle\beta _{6} = \frac{a^{3}a_{0}h^{3}v_{\ast }^{2}}{(1+a h u_{\ast })^{4}}, \beta _{7} = -\frac{2a^{2}a_{0}h^{2}v_{\ast }}{(1+a h u_{\ast })^{3}}, \beta _{8} = \frac{a a_{0}h}{(1+a h u_{\ast })^{2}}, \beta _{9} = \frac{a^{3}bh^{2}v_{\ast }}{(1+a h u_{\ast })^{4}}, \beta _{10}= \frac{-a^{2}bh}{(1+a h u_{\ast })^{3}}. \end{split} \end{eqnarray*}

Under assumptions $({\textbf{S}}_{\textbf{1}})$ , $({\textbf{S}}_{\textbf{3}})$ and $({\textbf{S}}_{\textbf{5}})$ , there exists $d_{2}^{\ast }$ that causes two Hopf bifurcation curves to intersect, with the intersection point represented as $(\tau ^{\ast },d_{2}^{\ast })$ , known as a double Hopf bifurcation point. Setting $\tau =\tau ^{\ast }+\mu _{1}, d_{2}=d_{2}^{\ast }+\mu _{2}$ , $U(t)=(u(t),v(t))^{T}$ , $\hat{U}(t)=(\hat{u}(t),\hat{v}(t))^{T}$ , $U_{t}({\cdot} )=U(t+\!{\cdot} )$ , and $\hat{U}_{t}({\cdot} )=\hat{U}(t+\!{\cdot} )$ . Then (4.1) becomes

(4.2) \begin{eqnarray} \begin{split} \frac{d }{d t}U(t)=&(\tau ^{\ast }+\mu _{1})\big [D_{0}\Delta U(t)+A_{0} U_{t}(0)+B_{0} \hat{U}_{t}(0)+C U_{t}({-}1)\big ]+\mu _{2}\tau ^{\ast }D_{1}\Delta U(t)\\ &+[f_{1}(\mu _{1},\mu _{2},U_{t},\hat{U}_{t}),f_{2}(\mu _{1},\mu _{2},U_{t},\hat{U}_{t})]^{T}, \end{split} \end{eqnarray}

where $D_{0}=\operatorname{diag}\{d_{1},d_{2}^{\ast }\}$ , $D_{1}=\operatorname{diag}\{0, 1\}$ , $\hat{U}_{t}=\frac{1}{l\pi }\int _{0}^{l\pi }U_{t}dx$ , and

\begin{eqnarray*} f_{1}(\mu _{1},\mu _{2},U_{t},\hat{U}_{t}) &=& (\tau ^{\ast }+\mu _{1})\big [\alpha _{1}u_{t}^{2}(0)+\tilde{\alpha }_{1}u_{t}(0)\hat{u}_{t}(0)+\alpha _{2}u_{t}(0)v_{t}(0)\\ &+&\alpha _{3}u_{t}^{3}(0)+\alpha _{4}u_{t}^{2}(0)v_{t}(0)\big ], \\[6pt] f_{2}(\mu _{1},\mu _{2},U_{t},\hat{U}_{t}) &=& (\tau ^{\ast }+\mu _{1})\big [\beta _{1}u_{t}^{2}(0)+\beta _{2}u_{t}(0)v_{t}(0)+\beta _{3}u_{t}^{2}({-}1)+\beta _{4}u_{t}({-}1)v_{t}(0)+\beta _{5}v_{t}^{2}(0)\\ &+&\beta _{6}u_{t}^{3}(0)+\beta _{7}u_{t}^{2}(0)v_{t}(0)+\beta _{8}u_{t}(0)v_{t}^{2}(0)+\beta _{9}u_{t}^{3}({-}1) +\beta _{10}u_{t}^{2}({-}1)v_{t}(0)\big ]. \end{eqnarray*}

At the double Hopf bifurcation point, the linearisation of system (4.2) has purely imaginary eigenvalues $\Lambda =\{\pm i\omega _{1}\tau ^{\ast },\pm i\omega _{2}\tau ^{\ast }\}$ , and all eigenvalues except $\Lambda$ have negative real parts under assumptions $({\textbf{S}}_{\textbf{1}})$ , $({\textbf{S}}_{\textbf{3}})$ and $({\textbf{S}}_{\textbf{5}})$ . To discuss double Hopf bifurcation, we assume

$({\textbf{H}}_{\textbf{1}})$ $\omega _{1}\lt \omega _{2}$ and $\omega _{1}\,{:}\,\omega _{2}$ is irrational number.

We choose

\begin{eqnarray*} \begin{array}{l} \displaystyle\phi _{1}(\theta ) =(1,p_{1})^{T} e^{i\omega _{1}\tau ^{\ast }\theta }, \phi _{2}(\theta ) =(1,p_{2})^{T} e^{i\omega _{2}\tau ^{\ast }\theta }, (\theta \in [{-}1,0]), \end{array} \end{eqnarray*}

as the column eigenvectors, corresponding to the eigenvalue $i\omega _{1}\tau ^{\ast }$ and $i\omega _{2}\tau ^{\ast }$ , respectively:

\begin{eqnarray*} \begin{array}{l} \displaystyle\psi _{1}(s) =M_{1}(1,q_{1}) e^{-i\omega _{1}\tau ^{\ast }s}, \psi _{2}(s) =M_{2}(1,q_{2})e^{-i\omega _{2}\tau ^{\ast }s}, (s\in [0,1]), \end{array} \end{eqnarray*}

as the row eigenvectors, corresponding to the eigenvalue $i\omega _{1}\tau ^{\ast }$ and $i\omega _{2}\tau ^{\ast }$ , respectively. By direct calculations, we obtain

\begin{eqnarray*} \begin{split} & p_{1} =(i\omega _{1}-\alpha _{11}-\beta _{11})/\alpha _{12},\,\, p_{2}=(i\omega _{2}+d_{1}/l^{2}-\alpha _{11})/\alpha _{12},\\[6pt] &q_{1}=\alpha _{12}/(i\omega _{1}-\alpha _{22}),\,\, q_{2}=\alpha _{12}/(i\omega _{2}+d_{2}/l^{2}-\alpha _{22}), \\[6pt] &M_{1}=(1+p_{1}q_{1}+\tau ^{\ast }q_{1}c_{21}e^{-i\omega _{1}\tau ^{\ast }})^{-1},\,\, M_{2}=(1+p_{2}q_{2}+\tau ^{\ast }q_{2}c_{21}e^{-i\omega _{2}\tau ^{\ast }})^{-1}. \end{split} \end{eqnarray*}

Decomposing the phase space $X_{\mathbb{C}}$ as:

\begin{eqnarray*} X_{\mathbb{C}}=P\oplus \textrm{ Ker}\pi, \end{eqnarray*}

where $P=\textrm{Im}\pi$ , and $\pi\,{:}\,X_{\mathbb{C}}\rightarrow P$ is the projection:

\begin{eqnarray*} \pi (U)=\sum _{k=1}^{2}\Phi _{k}(\Psi _{k},\langle \beta _n,U\rangle )\xi _{n}, \end{eqnarray*}

$U\in X_{\mathbb{C}}$ have the form:

(4.3) \begin{eqnarray} \begin{split} U=\sum _{k=1}^{2}(\Phi _{k}\tilde{z}_k(t))\xi _{n_k}+w=\Phi z^{x}+w, \end{split} \end{eqnarray}

with $\tilde{z}_k(t)=(\Psi _{k},\langle \beta _n,U\rangle )\in \mathbb{C}^2$ , $\Psi =(\Psi _{1},\Psi _{2})$ , $z^{x}=(z_1\xi _{n_1},z_2\xi _{n_1},z_3\xi _{n_2},z_4\xi _{n_2})^{T}$ , and $w\in \textrm{Ker}\pi$ . Let $z(t)=(z_1(t),z_2(t),z_3(t),z_4(t))^{T}\in \mathbb{C}^4$ and

\begin{eqnarray*} \widetilde{F}(z,w,\hat{w},\mu )=\widetilde{F}\Big (\sum _{k=1}^{2}(\Phi _{k}\tilde{z}_k(t))\xi _{n_k}+w, \frac{1}{l\pi }\int _{0}^{l\pi }\Big (\sum _{k=1}^{2}(\Phi _{k}\tilde{z}_k(t))\xi _{n_k}+w\Big )dx,\mu \Big ). \end{eqnarray*}

We use $(z,w,\hat{w},\mu )$ in place of $(U,\hat{U},\mu )$ , (4.2) can be re-represented with the following equation in $\mathbb{C}^{4}\times \textrm{Ker}\pi$ :

(4.4) \begin{eqnarray} \begin{cases} \displaystyle\dot{z}=Bz+\bar{\Psi }\begin{pmatrix}{c} \langle \beta _{n_1},\widetilde{F}(z,w,\hat{w},\mu )\rangle \\[6pt] \displaystyle\langle \beta _{n_2},\widetilde{F}(z,w,\hat{w},\mu )\rangle \\[6pt] \end{pmatrix},\\[6pt] \displaystyle\frac{d}{dt}w=\mathcal{L}_1w+(I-\pi )\widetilde{F}(z,w,\hat{w},\mu ), \end{cases} \end{eqnarray}

where $B=\textrm{diag}(B_1,B_2)$ , $B_1=\textrm{diag}(i\omega _1,-i\omega _1)$ , $B_2=\textrm{diag}(i\omega _2,-i\omega _2)$ .

Using the Taylor expansions:

\begin{eqnarray*} \widetilde{F}(U,\hat{U},\mu )=\sum _{j\geqslant 2}\frac{1}{j!}\widetilde{F}_j(U,\hat{U},\mu ), \mu \in \mathbb{R}^{2}, U\in X_{\mathbb{C}}, \end{eqnarray*}

where $\widetilde{F}_j$ is the $j$ th Fr $\acute{e}$ chet derivation of $\widetilde{F}$ , (4.4) can be recorded as:

(4.5) \begin{eqnarray} \begin{cases} \displaystyle\dot{z}=Bz+\mathop{\sum}\nolimits _{j\geqslant 2}\frac{1}{j!}f_{j}^{1}(z,w,\hat{w},\mu ),\\[6pt] \displaystyle\frac{d}{dt}w=\mathcal{L}_1w+\mathop{\sum}\nolimits_{j\geqslant 2}\frac{2}{j!}f_{j}^{1}(z,w,\hat{w},\mu ), \end{cases} \end{eqnarray}

where $\hat{w}=\frac{1}{l\pi }\int _{0}^{l\pi }wdx\in \textrm{Ker}\pi$ , and

(4.6) \begin{eqnarray} &&f_{j}^{1}(z,w,\hat{w},\mu )=\bar{\Psi }\begin{pmatrix} \langle \beta _{n_1},\widetilde{F}_{j}(z,w,\hat{w},\mu )\rangle \\[6pt] \langle \beta _{n_2},\widetilde{F}_{j}(z,w,\hat{w},\mu )\rangle \\[6pt] \end{pmatrix}, \end{eqnarray}

(4.7) \begin{align} f_{j}^{2}(z,w,\hat{w},\mu )=(I-\pi )\widetilde{F}_{j}(z,w,\hat{w},\mu ), \end{align}

From [Reference Chow and Hale4, Reference Faria8], it can be seen that through some variable transformations, we can obtain the normal form of (4.5). Let

\begin{eqnarray*} (z,w,\mu )=(\tilde{z},\tilde{w},\mu )+\frac{1}{j!}(U_{j}^1(\tilde{z},\mu ),U_{j}^2(\tilde{z},\mu )), j\geqslant 2, \end{eqnarray*}

and $U_{j}=(U_{j}^1,U_{j}^2)\in V_{j}^{4+2}(\mathbb{C}^4)\times V_{j}^{4+2}(\textrm{Ker}\pi )$ . $V_{j}^{4+2}(Y)$ is denoted as:

(4.8) \begin{eqnarray} V_{j}^{4+2}=\Big \{\sum _{|(p,l)=j|}c_{(p,l)}z^p\mu ^l\,{:}\,(p,l)\in \mathbb{N}_0^{4+2},c_{(p,l)}\in Y\Big \}, \end{eqnarray}

and $p=(p_1,p_2,p_3,p_4)\in \mathbb{N}_0^{4}$ , $l=(l_1,l_2)\in \mathbb{N}_0^{2}$ , $\sum _{i=1}^{4}p_i+\sum _{i=1}^{2}l_i=j$ , $z^p=z_1^{p_1}z_2^{p_2}z_3^{p_3}z_4^{p_4}$ and $\mu ^{l}=\mu ^{l_1}\mu ^{l_2}$ . The operators $M_j=(M_j^1,M_j^2)$ , $j\geqslant 2$ is defined as:

(4.9) \begin{eqnarray} \begin{split} &M_j^1\,{:}\,V_{j}^{4+2}(\mathbb{C}^4)\rightarrow V_{j}^{4+2}(\mathbb{C}^4),\\[6pt] &(M_j^1p)(z,\mu )=D_{z}p(z,\mu )Bz-Bp(z,\mu ),\\[6pt] &M_j^2\,{:}\,V_{j}^{4+2}(\textrm{Ker}\pi )\rightarrow V_{j}^{4+2}(\textrm{Ker}\pi ),\\[6pt] &(M_j^2h)(z,\mu )=D_{z}h(z,\mu )Bz-\mathcal{L}_1p(z,\mu ). \end{split} \end{eqnarray}

Through simplification, (4.5) can be written as:

(4.10) \begin{eqnarray} \begin{cases} \displaystyle\dot{z}=Bz+\mathop{\sum}\nolimits _{j\geqslant 2}\frac{1}{j!}g_{j}^{1}(z,w,\hat{w},\mu ),\\[6pt] \displaystyle\frac{d}{dt}w=\mathcal{L}_1w+\mathop{\sum}\nolimits_{j\geqslant 2}\frac{2}{j!}g_{j}^{1}(z,w,\hat{w},\mu ), \end{cases} \end{eqnarray}

with $g_{j}=(g_{j}^{1},g_{j}^{2})$ , $j\geqslant 2$ ,

\begin{equation*}g_{j}(z,w,\hat {w},\mu )=\bar {f}_{j}(z,w,\hat {w},\mu )-M_jU_j(z,\mu ).\end{equation*}

$U_j\in V_{j}^{4+2}(\mathbb{C}^4)\times V_{j}^{4+2}(\textrm{Ker}\pi )$ can be calculated by:

(4.11) \begin{eqnarray} U_j(z,\mu )=(M_j)^{-1}{\textbf{P}}_{\textrm{Im},j}\bar{f}_{j}(z,0,0,\mu ), \end{eqnarray}

where $(M_j)^{-1}$ is the inverse of $M_j$ . ${\textbf{P}}_{\textrm{Im},j}=({\textbf{P}}_{\textrm{Im},j}^1,{\textbf{P}}_{\textrm{Im},j}^2)$ is the projection operator related to previous decomposition of $V_{j}^{4+2}(\mathbb{C}^4)\times V_{j}^{4+2}(\textrm{Ker}\pi )\rightarrow \textrm{Im}(M_j^1)\times \textrm{Im}(M_j^2)$ .

4.1. Formula derivation of normal forms

To obtain the explicit form in the normal form, by (4.9), we have

(4.12) \begin{eqnarray} \begin{split} &M_j^1(z^p\mu ^le_k)=D_z(z^p\mu ^le_k)Bz-Bz^p\mu ^le_k=\\[6pt] &\begin{cases} (({-}1)^ki\omega _1+i\omega _1p_1-i\omega _1p_2+i\omega _2p_3-i\omega _2p_4)z^p\mu ^le_k, k=1,2,\\[6pt] (({-}1)^ki\omega _2+i\omega _1p_1-i\omega _1p_2+i\omega _2p_3-i\omega _2p_4)z^p\mu ^le_k, k=3,4, \end{cases} \end{split} \end{eqnarray}

where $e_k\in \mathbb{R}^4$ represents the $k$ th standard basic vector, $z^p=z_1^{p_1}z_2^{p_2}z_3^{p_3}z_4^{p_4}$ , and $\mu ^{l}=\mu ^{l_1}\mu ^{l_2}$ . Hence,

\begin{eqnarray*} \textrm{Ker}(M_2^1)=\textrm{span}\{\mu _iz_1e_1,\mu _iz_2e_2,\mu _iz_3e_3,\mu _iz_4e_4\}, i=1,2. \end{eqnarray*}

By (4.6), we have

(4.13) \begin{eqnarray} \frac{1}{2!}f_{2}^1(z,0,0,\mu )=\frac{1}{2!}\bar{\Psi }\begin{pmatrix} \langle \beta _{n_1},\widetilde{F}_{2}(z,0,0,\mu )\rangle \\[6pt] \langle \beta _{n_2},\widetilde{F}_{2}(z,0,0,\mu )\rangle \\[6pt] \end{pmatrix}. \end{eqnarray}

Therefore, on the centre manifold of $(0,0)$ near $(\tau ^{\ast }, d_{2}^{\ast })$ , for (4.2) the normal forms of second order is

\begin{eqnarray*} \dot{z}=Bz+\frac{1}{2!}g_{2}^1(z,0,0,\mu )+\ldots, \end{eqnarray*}

with $g_{2}^1(z,0,0,\mu )=\textrm{Proj}_{\textrm{Ker}(M_2^1)}f_{2}^1(z,0,0,\mu )$ . By further computation,

(4.14) \begin{eqnarray} \frac{1}{2!}g_{2}^1(z,0,0,\mu )=\frac{1}{2!}\textrm{Proj}_{\textrm{Ker}(M_2^1)}f_{2}^1(z,0,0,\mu )= \begin{pmatrix} E_{11}\mu _1z_1+E_{21}\mu _2z_1 \\[6pt] \overline{E_{11}}\mu _1z_2+\overline{E_{21}}\mu _2z_2 \\[6pt] E_{13}\mu _1z_3+E_{23}\mu _2z_3 \\[6pt] \overline{E_{13}}\mu _1z_4+\overline{E_{23}}\mu _2z_4 \\[6pt] \end{pmatrix}, \end{eqnarray}

where $E_{11}$ , $E_{21}$ , $E_{13}$ and $E_{23}$ have the following representations:

\begin{eqnarray*} \begin{split} &E_{11}=\psi _{1}(0)\big (A_{0}\phi _{1}(0)+B_{0}\phi _{1}(0)+C\phi _{1}({-}1)\big ),\\[6pt] &E_{21}=0, \\[6pt] &E_{13}=\psi _{2}(0)\big ({-}D_{0}\phi _{2}(0)/l^{2}+A_{0}\phi _{2}(0)+C\phi _{2}({-}1)\big ),\\[6pt] &E_{23}=-\psi _{2}(0)\tau ^{\ast }D_{1}\phi _{2}(0)/l^{2}, \end{split} \end{eqnarray*}

where $A_{0}$ , $B_{0}$ and $C$ are given by (3.1), $D_{0}=\operatorname{diag}\{d_{1},d_{2}^{\ast }\}$ and $D_{1}=\operatorname{diag}\{0, 1\}$ .

Similarly, by (4.5) and (4.12), we can obtain the expression of the third-order normal forms:

\begin{eqnarray*} \dot{z}=Bz+\frac{1}{2!}g_{2}^1(z,0,0,\mu )+\frac{1}{3!}g_{3}^1(z,0,0,0)+\ldots, \end{eqnarray*}

by noting that

\begin{equation*}\textrm {Ker}(M_{3}^{1})=\textrm {span}\{e_1z_1^2z_2,e_1z_1z_3z_4,e_2z_1z_2^2,e_2z_2z_3z_4,e_3z_3^2z_4, e_3z_1z_2z_3,e_4z_3z_4^2,e_4z_1z_2z_4\}.\end{equation*}

Here, $g_{3}^1(z,0,0,0)=\textrm{Proj}_{\textrm{Ker}(M_3^1)}\bar{f}_{3}^1(z,0,0,0)$ , $\bar{f}_{3}^1(z,0,0,0)$ can be obtained by the formula (see [Reference Shen, Liu and Wei22]):

(4.15) \begin{eqnarray} \begin{split} \bar{f}_{3}^1(z,0,0,0)&=f_{3}^1(z,0,0,0)+\frac{3}{2}\Big (D_zf_2^1(z,0,0,0)U_2^1(z,0)+D_wf_2^1(z,0,0,0)U_2^2(z,0)\\ &+D_{\hat{w}}f_2^1(z,0,0,0)\hat{U}_2^2(z,0)-D_zU_2^1(z,0)g_{2}^1(z,0,0,0)\Big ), \end{split} \end{eqnarray}

where $U_2^1$ and $U_2^2$ are given in (4.11), and

\begin{equation*}U_2^2=\frac {1}{l\pi }\int _{0}^{l\pi }U_2^2dx.\end{equation*}

It is straightforward to have $g_{2}^1(z,0,0,0)=0$ by (4.14); next, we need to calculate the remaining four parts in formula (4.15):

\begin{eqnarray*} &&(a)\,\,\textrm{Proj}_{\textrm{Ker}(M_3^1)}f_{3}^1(z,0,0,0),\\ &&(b)\,\,\textrm{Proj}_{\textrm{Ker}(M_3^1)}\big (D_z\,f_2^1(z,0,0,0)U_2^1(z,0)\big ),\\ &&(c)\,\,\textrm{Proj}_{\textrm{Ker}(M_3^1)}\big (D_w\,f_2^1(z,0,0,0)U_2^2(z,0)\big ),\\ &&(d)\,\,\textrm{Proj}_{\textrm{Ker}(M_3^1)}\big (D_{\hat{w}}\,f_2^1(z,0,0,0)\hat{U}_2^2(z,0)\big ), \end{eqnarray*}

which are provided in the appendix.

Using the algorithms similar to those in [Reference Du, Niu, Guo and Wei6, Reference Shen, Liu and Wei22], we can get the following normal form of double Hopf bifurcation, up to the third order:

(4.16) \begin{eqnarray} \begin{cases} \begin{split} &\dot{z}_{1}=i\omega _{1}\tau ^{\ast }z_{1}+E_{11}\mu _{1}z_{1}+E_{21}\mu _{2}z_{1}+E_{2100}z_{1}^{2}z_{2}+E_{1011}z_{1}z_{3}z_{4},\\[6pt] &\dot{z}_{2}=-i\omega _{1}\tau ^{\ast }z_{2}+\overline{E_{11}}\mu _{1}z_{2}+\overline{E_{21}}\mu _{2}z_{2} +\overline{E_{2100}}z_{1}z_{2}^{2}+\overline{E_{1011}}z_{2}z_{3}z_{4},\\[6pt] &\dot{z}_{3}=i\omega _{2}\tau ^{\ast }z_{3}+E_{13}\mu _{1}z_{3}+E_{23}\mu _{2}z_{3}+E_{0021}z_{3}^{2}z_{4}+E_{1110}z_{1}z_{2}z_{3},\\[6pt] &\dot{z}_{4}=-i\omega _{2}\tau ^{\ast }z_{4}+\overline{E_{13}}\mu _{1}z_{4}+ \overline{E_{23}}\mu _{2}z_{4}+\overline{E_{0021}}z_{3}z_{4}^{2}+\overline{E_{1110}}z_{1}z_{2}z_{4}. \end{split} \end{cases} \end{eqnarray}

By the polar coordinate transformations, let $z_1=\tilde{r}_1e^{i\theta _1}$ , $z_2=\tilde{r}_2e^{i\theta _2}$ , and

\begin{eqnarray*} &&\epsilon _1=\textrm{Sign}(\textrm{Re}(E_{2100})),\epsilon _2=\textrm{Sign}(\textrm{Re}(E_{0021})),\\ &&r_1=\tilde{r}_1\sqrt{|E_{2100}|},r_2=\tilde{r}_2\sqrt{|E_{0021}|},\tilde{t}=t\epsilon _1, \end{eqnarray*}

the system (4.16) can be rewritten as follows:

(4.17) \begin{eqnarray} \begin{cases} \begin{array}{ll} \displaystyle\frac{d r_{1}}{d \tilde{t}}=r_{1}(\upsilon _{1}+r_{1}^{2}+b_{0}r_{2}^{2}), \\[6pt] \displaystyle\frac{d r_{2}}{d \tilde{t}}=r_{2}(\upsilon _{2}+c_{0}r_{1}^{2}+d_{0}r_{2}^{2}), \end{array} \end{cases} \end{eqnarray}

where

(4.18) \begin{eqnarray} &&\upsilon _{1}=\epsilon _1\big (\textrm{Re}(E_{11})\mu _1+\textrm{Re}(E_{21})\mu _2\big ), \nonumber \\ &&\upsilon _{2}=\epsilon _1\big (\textrm{Re}(E_{13})\mu _1+\textrm{Re}(E_{23})\mu _2\big ),\nonumber \\ &&b_0=\frac{\epsilon _1\epsilon _2\textrm{Re}(E_{1011})}{\textrm{Re}(E_{0021})}, c_0=\frac{\textrm{Re}(E_{1110})}{\textrm{Re}(E_{2100})},d_0=\epsilon _1\epsilon _2. \end{eqnarray}

Obviously, there is a zero equilibrium $E_1(0,0)$ . The three nonnegative equilibria:

(4.19) \begin{eqnarray} &&E_2=(\sqrt{-\upsilon _{1}},0), \textrm{when}\,\upsilon _{1}\lt 0, \nonumber \\ &&E_3=\Big (0,\sqrt{-\frac{\upsilon _{2}}{d_0}}\Big ), \textrm{when}\,d_0\upsilon _{2}\lt 0,\nonumber \\ &&E_4=\Big (\sqrt{\frac{b_0\upsilon _{2}-d_0\upsilon _{1}}{d_0-b_0c_0}}, \sqrt{\frac{c_0\upsilon _{1}-\upsilon _{2}}{d_0-b_0c_0}}\Big ), \textrm{when}\,\frac{b_0\upsilon _{2}-d_0\upsilon _{1}}{d_0-b_0c_0}\gt 0,\frac{c_0\upsilon _{1}-\upsilon _{2}}{d_0-b_0c_0}\gt 0. \end{eqnarray}

Corresponding to the original system (1.2), $E_1$ associates with the positive steady state; $E_2$ and $E_3$ represent the spatially periodic solutions, and $E_4$ represents the spatially quasi-periodic solution. Due to the possible different sign in the coefficients $b_0$ , $c_0$ , $d_0$ and $d_{0}-b_{0}c_{0}$ , there are 12 types of unfolding in (4.17) (see chapter 7.5 in [Reference Guckenheimer and Holmes12]).

5. Numerical simulations

Throughout this section, we always fix the parameters in (2.1) and vary $d_{2}$ , $\tau$ and $a_0$ to explore the dynamics with respect to diffusion, time delay and predator cannibalism on the spatial distribution of prey. Predators can also exhibit similar pattern structures.

Figure 7. The bifurcation diagrams on the $\tau -d_{2}$ plane.

Figure 8. The dynamical phenomenon of different regions.

5.1. The dynamical phenomenon in different regions

We first draw the bifurcation diagram with respect to $d_{2}$ and $\tau$ in Figure 7, where we can see with the increase of the diffusion coefficient $d_{2}$ , the double Hopf bifurcation points appear. We can calculate $(\tau ^{\ast },d_{2}^{\ast })=(8.1471,0.3499)$ (denoted by “ $\textrm{HH}_{\textrm{1}}$ ”), and

\begin{eqnarray*} &&E_{11}=0.0099+0.0961i, E_{21}=0, \\ &&E_{13}=0.0526+0.0980i, E_{23}=-0.6613-0.1267i,\\ &&E_{2100}=-0.0021-0.002i, E_{1011}= -0.0128+0.0103i,\\ &&E_{0021}=-0.0093-0.0022i, E_{1110}=-0.0083+0.0009i. \end{eqnarray*}

It follows from (4.16) and (4.17) that

(5.1) \begin{eqnarray} \epsilon _{1}=-1, b_{0}=1.3696, c_{0}=3.9506, d_{0}=1, d_{0}-b_{0}c_{0}=-4.4110. \end{eqnarray}

Near the double Hopf bifurcation point $(\tau ^{\ast },d_{2}^{\ast })$ , the dynamics of system (1.2) is topologically equivalent to (4.17) near $(\mu _{1},\mu _{2})=(0,0)$ , where $(\mu _{1},\mu _{2})=(\tau -\tau ^{\ast },d_{2}-d_{2}^{\ast })$ . We can divide $\mu _{1}-\mu _{2}$ plane into six parts by lines $l_{1}$ ( $\upsilon _{1}=0$ ), $l_{2}$ ( $\upsilon _{2}=0$ ), and half-lines $l_{3}$ ( $b_0\upsilon _{2}-d_0\upsilon _{1}=0$ ), $l_{4}$ ( $c_0\upsilon _{1}-\upsilon _{2}=0$ ), obtained from (4.18) and (4.19), as depicted in Figure 8. The four red dots marked in Figure 8 are the relative positions of the values. We then calculate that

\begin{eqnarray*} &&l_{1}\,{:}\, \mu _{1}=0,\,l_{2}\,{:}\,\mu _{1}=12.5722\mu _{2},\\ &&l_{3}\,{:}\,\mu _{1}=14.5638\mu _{2}\,(\mu _{2}\geqslant 0),\\ &&l_{4}\,{:}\,\mu _{1}=49.0891\mu _{2}\, (\mu _{2}\geqslant 0). \end{eqnarray*}

Near the double Hopf bifurcation point $(\tau ^{\ast },d_{2}^{\ast })$ , the values of time delay and diffusion coefficient taken in the numerical simulation are consistent with the theoretical results. More specially, the dynamical classifications in each region separating by $l_i$ $(i=1,2,3,4)$ are shown in Figure 9.

Figure 9. The dynamical classifications in regions 1, 2, 4 and 5.

Figure 10. When $\tau =6\lt \tau ^{\ast }$ , the positive constant stationary solution $E_{\ast }$ of system (1.2) is locally asymptotically stable.

5.2. The effect of predator diffusion on pattern formation

From Figure 7, we know that the positive equilibrium remains stable if $\tau \lt \operatorname{min}\{\tau ^{0+},\tau _{1}^{0+}\}$ , spatially homogeneous or nonhomogeneous periodic solutions appear if $\tau \gt \operatorname{min}\{\tau ^{0+},\tau _{1}^{0+}\}$ . We then choose the time delay $\tau$ and the diffusion coefficient $d_{2}$ as double Hopf bifurcation parameters to study the dynamic behaviour of the system (1.2). With the same initial functions $u(x,t)=6.5+0.01\cos x$ , $v(x,t)=5.5+0.01\cos x$ , when we fix $d_{2}=0.8$ , we can observe different patterns as $\tau$ varies. In Figure 10, the positive equilibrium is locally asymptotically stable. In Figures 11 and 12 ( $d_{2}=0.34$ ), there are stable spatially homogeneous and nonhomogeneous periodic solutions. In Figure 13, we find stable spatially homogeneous and nonhomogeneous periodic oscillations could coexist with different initial functions. One initial function is $u(x,t)=6.48+0.01\cos x$ , $v(x,t)=5.53$ , and the other is $u(x,t)=6.52$ , $v(x,t)=5.53+0.01\cos x$ . With the increase of time, the unstable quasi-periodic solution disappears gradually.

Figure 11. When $\tau =13.9\gt \tau ^{\ast }$ , spatially homogeneous periodic solutions are stable of system (1.2).

Figure 12. When $\tau =18.7\gt \tau ^{\ast }$ , spatially nonhomogeneous periodic solutions are stable of system (1.2).

Figure 13. Stable spatially homogeneous periodic solution and stable spatially nonhomogeneous periodic solution with different initial values coexist when $\tau =18.7$ .

Remark 5.1. Under the assumptions $({\textbf{S}}_{\textbf{1}})$ , $({\textbf{S}}_{\textbf{2}})$ and $({\textbf{S}}_{\textbf{5}})$ , there exists $d_{2}^{\ast }$ such that $\tau ^{0+}=\tau _{1}^{0+}$ , system (1.2) undergoes a double Hopf bifurcation at the positive equilibrium $E_{\ast }$ , when $d_{2}=d_{2}^{\ast }$ and $\tau ^{\ast }=\tau ^{0+}=\tau _{1}^{0+}$ . Notice that $d_{2}=0.5$ and $l=1.2$ , the relationship of $b$ and $\tau$ is shown in Figure 14. The double Hopf bifurcation point $(\tau ^{\ast },b^{\ast })=(9.2949,0.9925)$ is denoted by “HH”. In fact, the bifurcation set near HH has the form given in Figure 8, which can be proved by using the normal form method in the coming section.

Figure 14. The bifurcation diagrams on the $b-\tau$ plane.

5.3. The effect of cannibalism on pattern formation

In this section, we discuss the effect of cannibalism. Fix $a=0.28$ , $\tau =18.7$ . When the intraspecies self-killing rate is less than the interspecies capture rate ( $a_0\lt a$ ), the system exhibits a spatially homogeneous periodic solution. When the intraspecies self-killing rate and interspecies capture rate are equal ( $a_0=a$ ), the system exhibits spatially nonhomogeneous periodic solutions. The increase in $a_0$ results in a spatial distribution of species with certain patterns. When the intraspecies self-killing rate exceeds the interspecies capture rate, that is, $a_0\gt a$ , the system exhibits a stable equilibrium solution; see Figure 15(c). Predators eventually reach a stable number, indicating that cannibalism may lead to a decline in the population, but it may enable stronger individuals to survive. This is beneficial in an environment where food is relatively scarce, ensuring that a small number of excellent individuals can master sufficient resources to breed the next generation.

Figure 15. The spatiotemporal diagram of predator. Left, $a_0=0.2$ ; centre, $a_0=0.28$ ; right, $a_0=0.3$ .