The dynamical analysis of a nonlocal predator–prey model with cannibalism

Abstract

Cannibalism is often an extreme interaction in the animal species to quell competition for limited resources. To model this critical factor, we improve the predator–prey model with nonlocal competition effect by incorporating the cannibalism term, and different kernels for competition are considered in this model numerically. We give the critical conditions leading to the double Hopf bifurcation, in which the gestation time delay and the diffusion coefficient were selected as the bifurcation parameters. The innovation of the work lies near the double Hopf bifurcation point, and the stable homogeneous and inhomogeneous periodic solutions can coexist. The theoretical results of the extended centre manifold reduction and normal form method are in good agreement with the numerical simulation.

1. Introduction

Cannibalism is the biological phenomenon of consuming the same species and helping to provide a food source. It is commonly seen that some multiparous animals and some primates chew the young animals in the first few days after delivery [5, 20, 27]. For example, a variety of spiders, mantis, chironomids and sea slugs have the special habit of females eating males. As a special natural phenomenon, cannibalism often occurs in social insects [13], spiders [18] and plankton [23], which has attracted the attention of many scholars. The mutual killing between brothers and sisters for survival may happen in the birds’ or small animals’ nests. For the phenomenon of sexual cannibalism [11, 26, 29], some people believe that it is due to various complex evolutionary reasons, but after research, it is found that the motivation for this creepy cannibalism is straightforward, often due to the size of male and female. Large female spiders eat their weak mates for two reasons: hunger and the ability to eat smaller male spiders [28].

In our work [7], we investigated the following predator–prey model with cannibalism effect and time delay:

(1.1) \begin{eqnarray} \begin{cases} \displaystyle\frac{\partial }{\partial t}u(x,t)=d_{1} \Delta u(x,t)+ru(x,t)\Big [1-\frac{u(x,t)}{K}\Big ]-\frac{au(x,t)v(x,t)}{1+ahu(x,t)}, \\[8pt] \displaystyle \frac{\partial }{\partial t}v(x,t)=d_{2} \Delta v(x,t)+v(x,t)\Big [\frac{bau(x,t-\tau )}{1+ahu(x,t-\tau )}-d-\frac{a_{0}v(x,t)}{1+ahu(x,t)}\Big ], & x\in \Omega, t\gt 0, \end{cases} \end{eqnarray}

where $u(x,t)$ and $v(x,t)$ denote the population density of the prey and the predator at location $x$ and time $t$ , respectively. $\tau$ is the gestation time of predator. $d_{1}$ and $d_{2}$ are diffusion coefficients. $r$ is the growth rate and $K$ is the carrying capacity of prey. $a$ and $a_0$ are the scalings of the predator–prey and predator–predator encounter rate, respectively. $h$ is the handling time and $d$ is death rate of the predator. $b$ is the food-to-newborn conversion factor. For the diffusive system (1.1), Sun et al. [25] considered spatial pattern formation induced by predator cannibalism in two-dimensional spatial domain when $\tau =0$ . Li et al. [14] investigated instability and Hopf bifurcation of system (1.1) induced by time delay with normal form method and centre manifold theory. Based on [14], we have studied the double Hopf bifurcation induced by two different parameters and obtained the coexistence of periodic solutions caused by cannibalism.

In nature, individual species compete for common resources not only in the neighbouring regions but also in the entire spatial domain. Therefore, nonlocal competition effect can be introduced into the biological model to make the system more reasonable [2, 9, 17, 24]. In [2], a single population model with nonlocal competition effect is first studied, where Britton introduced the nonlocal competition term $\int _{\Omega }W(x-y)u(y,t)dy$ with $\Omega =({-}\infty,\infty )$ and $W(x-y)$ is general kernel function. In the simple case, Furter and Grinfeld [10] replaced $W(x-y)$ with $W(x,y)$ and took $W(x,y)=1/\left |\Omega \right |$ when $\Omega$ is bounded and investigated the steady-state bifurcation of a single population model with Neumann boundary conditions. Since then, the nonlocal competition has been broadly studied in competitive population [17, 21] and predator–prey models [1, 3, 15, 16]. Many results show that nonlocal prey competition can make the system more likely to generate spatiotemporal dynamics and affect the stability of wave trains.

What kind of interesting dynamic behaviours will be generated, when we consider the influence of nonlocal effects on the system (1.1)? In nature, the prey is often small with abundant amounts and competition for intraspecies resources, so the nonlocal competition is more intense and common for the prey, while the number of predators is less than that of the prey, and the size is larger, hence the intraspecific competition is not considered normally. Keep this in mind, we propose a diffusive model with nonlocal intraspecific competition in prey:

(1.2) \begin{eqnarray} \begin{cases} \displaystyle\frac{\partial u(x,t)}{\partial t}=d_{1} \Delta u(x,t)+ru(x,t)\Big [1-\frac{1}{K}\int _{\Omega }W(x,y)u(y,t)dy\Big ]-\frac{au(x,t)v(x,t)}{1+ahu(x,t)}, \\[12pt] \displaystyle\frac{\partial v(x,t)}{\partial t}=d_{2} \Delta v(x,t)+v(x,t)\Big [\frac{bau(x,t-\tau )}{1+ahu(x,t-\tau )}-d-\frac{a_{0}v(x,t)}{1+ahu(x,t)}\Big ], \,\,x\in \Omega, t\gt 0, \\[12pt] \displaystyle\frac{\partial u(x,t)}{\partial \overrightarrow{n}} =\frac{\partial v(x,t)}{\partial \overrightarrow{n}}=0, \,\,x\in \partial \Omega, t\gt 0,\\[12pt] \displaystyle u(x,t)=u_{0}(x,t)\geqslant 0, v(x,t)=v_{0}(x,t)\geqslant 0, \,\,x\in \overline{\Omega },t\in [{-}\tau,0]. \end{cases} \end{eqnarray}

For simplicity of notation, we denote $\hat{u}(x,t)=\int _{\Omega }W(x,y)u(y,t)dy$ . To facilitate the study of the distribution of eigenvalues, we choose $\Omega =(0,l\pi )$ , $l\in \mathbb{R}^{+}$ . If the kernel is a classical Dirac function, that is, $W(x,y)=\delta (x-y)$ , the system (1.2) is reduced to the system (1.1). Some other common kernel functions can be used in the model, such as symmetric kernel function, average kernel function, etc, which are commonly used to describe the probability distribution of random variables [9]. Therefore, we choose the following particular kernels:

\begin{eqnarray*} &&\textrm{(I)}\,W(x,y)=\frac{1}{l\pi }+\mu \cos\!(x-y),\\ &&\textrm{(II)}\, W(x,y)=g_1(x)g_2(y). \end{eqnarray*}

Compared with the research on the reaction–diffusion system without nonlocal competition effect [7], we first show the effects of different kernel functions on the system dynamics. Through intuitive simulation, it is found that the system will exhibit constant/nonconstant steady-state solutions and spatial homogeneous/nonhomogeneous periodic solutions. This means that due to the different forms of kernel functions, there are significant differences in the final dynamics of the population. Second, we find the existence of double Hopf bifurcations for the average kernel function. In addition, we can provide the theoretical results near the bifurcation points of double Hopf by using the generalised normal form method and demonstrate the complex dynamic behaviour that the system may present numerically.

After the introduction, in Section 2, we show, by selecting different nonuniform kernel functions, that the system can generate a stable constant/nonconstant steady-state solution, and a stable nonconstant periodic solution, etc. In Section 3, we mainly discuss the existence of double Hopf bifurcation near the positive steady state. In Section 4, near the double Hopf bifurcation point, we provide the extended calculation methods of normal form for the nonlocal predator–prey system with time delay and predator diffusion. Furthermore, we study the coexistence behaviour of spatially nonhomogeneous and homogeneous periodic oscillations near the double Hopf bifurcation point based on the influence of diffusion and cannibalism on the system.

In summary, we establish a framework to calculate the normal form of the centre manifold of system (1.2) at double Hopf bifurcation points. Some recursive transformations related to the variables are used to derive the codimension two normal form, which is a tedious but effective derivation process. We provide the range of related parameters through a normal form bifurcation set, and the selection of parameters is determined by several different bifurcation curves. As for the initial function, due to the local dynamical analysis, we choose it near the steady-state solution. The numerical simulation results provide good support for our theoretical analysis.

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 [7]:

(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.

Lemma 1. Under the assumption (S ₁ ), the positive equilibrium $E_{\ast }$ is asymptotically stable in the system (1.2) when $\tau =0$ .

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 [19], 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.

Lemma 2. Suppose that (S ₁ ) holds.

$\textrm{(a)}$ If (S ₃ ) and (S ₅ ) hold, then $\operatorname{Re} \lambda '(\tau ^{j+})\gt 0$ , $\operatorname{Re} \lambda '(\tau _{n}^{j+})\gt 0$ for $j\in \mathbb{N}_{0}$ , $n\in \mathbb{N}_{+}$ .

$\textrm{(b)}$ If (S ₂ ) and (S ₄ ) hold, then $\operatorname{Re} \lambda '(\tau ^{j+})\gt 0$ , $\operatorname{Re} \lambda '(\tau ^{j-})\lt 0$ , $\operatorname{Re} \lambda '(\tau _{n}^{j+})\gt 0$ , $\operatorname{Re} \lambda '(\tau _{n}^{j-})\lt 0$ for $j\in \mathbb{N}_{0}$ , $n\in \mathbb{N}_{+}$ .

Remark 3.1. The possible cases are $({\textbf{S}}_{\textbf{2}})$ and $({\textbf{S}}_{\textbf{5}})$ , $({\textbf{S}}_{\textbf{3}})$ and $({\textbf{S}}_{\textbf{4}})$ , a conclusion similar to Lemma 2 can be obtained as well.

Remark 3.2. Under assumptions $({\textbf{S}}_{\textbf{1}})$ , $({\textbf{S}}_{\textbf{3}})$ and $({\textbf{S}}_{\textbf{5}})$ , system (1.2) may undergo a double Hopf bifurcation at the positive constant stationary solution $E_{\ast }$ when $(\tau,d_{2})=(\tau ^{\ast },d_{2}^{\ast })$ . In order to better analyse the dynamical behaviour near the bifurcation point, we will give the normal form analysis of $E_{\ast }$ about system (1.2).

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 [22], 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 [4, 8], 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 [22]):

(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 [6, 22], 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 [12]).

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$ .

6. Conclusions

Compared with the model without nonlocal prey competition, we find that the predator–prey model with the nonlocal term can generate new dynamic behaviours, such as nonhomogeneous stable periodic patterns. We find that the dynamic behaviour of the predator and the prey is directly affected by the shape of the kernel function. Different types of kernel functions can affect the population’s dynamic behaviour density over time and space. In order to further investigate the dynamics with different kernel functions, we take uniform kernel functions as an example and provide theoretical analysis. The main derivation is the form of the double Hopf bifurcation normal form near the positive equilibrium point of the system, and the simulation results page displays the complex dynamic behaviour of the system. The highlights are mainly divided into two parts. First, the coexistence of spatially homogeneous and nonhomogeneous stable periodic oscillations is one of the highlights of this paper. Another highlight is the derivation of the normal form at the double Hopf bifurcation point. This normal form algorithm can be extended to the general reaction–diffusion model, which can be our next work.