Multi-stable and spatiotemporal staggered patterns in a predator-prey model with predator-taxis and delay

1.   Introduction

It is common for population dynamics to be presented by using random interactions of species in space ^[1,2,3]. Often, however, changes in the interactions and spatial distributions of populations can cause the population dynamics to change fundamentally ^[4,5,6,7]. A characteristic of living systems is their ability to respond to changes in the environment and approach or move away from environmental stimuli. The importance of chemotaxis in the modeling of a variety of biological and ecological processes has been recognized, including the spread of epidemics, aggregations of cellular slime mold, dynamics of planktonic communities and dynamic of insect populations ^[8,9,10,11]. For predator and prey populations in space, the predator usually spends more time in places with more abundant food resources, and such behavior has a clear selective advantage. It directly leads to the relative aggregation of the predator in areas with a high prey density ^[12]. Although experiments on spatial models are limited, this aggregative behavior has been inscribed in some predator-prey models since as early as the 1970s ^{[13,14,15,16]}. The results suggest that the aggregation of the predator in regions of high prey density is important to the dynamics of predator-prey interactions. In 1974, Hassell and May ^[17] used laboratory experiments and field observations of birds, insect predators and insect parasites to show the effect of predator aggregation on the dynamic properties of the two populations. It turns out that stability is increased when the prey distribution becomes more clumped. It was not until 1987 that Kareiva and Odell ^[18] explicitly characterized the aggregative behavior of the predator in terms of the non-random distribution of predators in space, yielding the first ecological model to describe the directional movement of species using prey-taxis. Since then, some mathematical models have been proposed and the effects of prey-taxis on the spatiotemporal distributions of populations have been considered; see, for example ^[19,20,21] and the references therein. It has been shown that a strong prey-taxis can prevent population overcrowding. That is to say, prey-taxis has a stabilizing effect on population survival, which also supports the experimental results in ^[17]. Corresponding to the chemotactic movement of the predator is the chemotactic movement of the prey (namely predator-taxis), which refers to the fact that the prey can sense the risk of being preyed upon and will thus move in the opposite direction of the predator distribution gradient. The existing research about the directed movements of prey have also provided new insights into the emergence of spatial non-uniform distribution of species. For example, Wang and Zou ^[22] have demonstrated that predator-taxis is a positive anti-predator behavior exhibited by the prey, because it can inhibit the formations of spatial heterogeneity. And in ^[23], Wu et al. analogized the avoidance of infected populations by susceptible populations to the phenomenon of predator-taxis, demonstrating that chemotaxis reduces the probability of spatial heterogeneity and thus reduces the risk of influenza transmission. In addition, the time delay caused by the conversion of the predator's capture behavior into the predator's growth has also often been taken into account ^[24,25], and it is often considered as one of the mechanisms leading to the emergence of time-periodic oscillations in the system.

Hence, incorporating chemotaxis and time delay into the classical predator-prey model is necessary and reasonable, and the corresponding model is written as follows

where $u_{\tau} = u(x, t-\tau), v_{\tau} = v(x, t-\tau)$, $\Omega\subset\mathbb{R}^{n}$ and $\vec{n}$ is the outer unit normal vector on $\partial\Omega$. $u(x, t)$ and $v(x, t)$ represent the densities of the prey and predator at location $x$ and time $t$ respectively, $d_{u}$ and $d_{v}$ denote the random dispersal rates of the prey and predator respectively, and thus are positive, and $c > 0$ reflects the conversion rate. The function $f_{1}(u)$ describes the growth rate of the prey, with the most common cases being the constant growth rate $f_{1}(u) = r_{0} > 0$, logistic growth rate $f_{1}(u) = r_{0}(1-\frac{u}{N})$ with $N > 0$, etc. $f_{2}(v)$ describes the mortality rate of the predator, with the most common cases being the constant mortality rate $f_{2}(v) = m_{1} > 0$, linear mortality rate $f_{2}(v) = m_1+m_{2}v$ with $m_{2} > 0$, etc. Moreover, $p(u, v)$ is the functional response function, and some common cases are as follows:

where $b_{1} > 0$ and $b_{2} > 0$. In particular, the time delay $\tau\geqslant0$ is regarded as the contribution of the predation that occurred in the past to the current growth of the predator ^[26,27]. $\nabla\cdot(\alpha\zeta_{1}(u)\nabla v)$ describes the movement of prey toward a lower density of predator, while $-\nabla\cdot(\eta\zeta_{2}(v)\nabla u)$ describes the predator movement toward a higher density of prey, where $\alpha, \; \eta\geqslant0$ are chemotaxis coefficients, and the sensitivity function $\zeta_{1}(u)$ can be chosen as in ^[28], e.g.,

where $\epsilon > 0, \; m\geqslant1$ and $\zeta_{2}(v)$ can be selected analogously.

For model (1.1) with $\zeta_{1}(u) = u, \; \zeta_{2}(v) = v$ and $\tau = 0$, spatial patterns induced by prey-taxis and predator-taxis have been discussed by Wang et al. ^[29], and they showed that spatial patterns can be eliminated by strong predator-taxis and prey-taxis. Results on Turing instability can be also found in a paper by Cao and Wu ^[30], where sufficient conditions for the emergence of spatial patterns were provided and the authors further revealed that the appearances of spatial patterns does not necessarily require that the predator self-diffusion rate be greater than that of the prey in the presence of chemotaxis.

As for model (1.1) with only prey-taxis ($\alpha = 0, \; \eta > 0$) and $\tau = 0$, the global existence, asymptotic behavior or blow-up of solutions in general parabolic-parabolic systems with prey-taxis have been widely studied; for example, see ^[31,32] and the references therein. Pattern formations induced by prey-taxis with $\zeta_{2}(v) = v$ have also been discussed in detail in ^[33] for different $f_1, \; f_{2}$ and $p$. Specifically speaking, Lee et al. ^[33] showed that strong prey-taxis tends to stabilize the coexistence equilibrium when $f_1, \; f_{2}$ and $p$ correspond to the logistic growth rate, constant mortality rate and ratio-dependent forms, respectively. Moreover, Gao and Guo ^[34] demonstrated that the local stability of the constant steady state is enhanced by the presence of prey-taxis. Subsequently, Qiu et al. ^[35] showed that prey-taxis can suppress the global asymptotic stability of the coexistence steady state, and they pointed out that due to the effect of prey-taxis, periodic solutions bifurcating from the coexistence steady state via Hopf bifurcation can be spatially inhomogeneous.

Also, there has been some research on model (1.1) with only predator-taxis ($\alpha > 0, \; \eta = 0$) and $\tau = 0$. For example, Wu et al. ^[23] discussed the general model with predator-taxis and proved that strong predator-taxis can make the spatial patterns caused by self-diffusion disappear.

Particularly, when considering the logistic growth rate $f_{1}(u)$, constant mortality rate $f_{2}(v)$, ratio-dependent functional response $p(u, v)$ and linear sensitivity function $\zeta_{1}(u)$, model (1.1) can be written as follows:

where $\Omega = (0, \pi)$. $r_{0}$ and $\frac{r_{0}}{a}$ denote the intrinsic growth rate of prey and carrying capacity, respectively, $m_{1}$ is the mortality rate of the predator, which is independent of its density, and $b_{1}$ and $b_{2}$ denote the capturing rate and half saturation constant, respectively. The initial functions $u_{0}$ and $v_{0}$ are non-negative, continuous and satisfy

For model (1.2) with $\tau = 0$, by replacing $\nabla\cdot(\alpha u\nabla v)$ with $\nabla\cdot(\alpha(1-\frac{u}{M})u\nabla v)$, Wang and Zou ^[22] investigated the formations of spatial patterns and found that weak predator-taxis can lead to spatial patterns, where $M$ measures the maximum number of prey that a unit volume can accommodate. As for (1.2) with $\tau = 0$ and a general form of predator-taxis, Gao ^[36] proved the global existence and uniform boundedness of the classical solutions. For the case with no chemotaxis, that is, for $\alpha = 0$, Song et al. ^[37] considered the delayed model (1.2) and discussed the existence and stability of the delay-induced spatially homogeneous periodic orbit.

In this paper, we would like to reveal the diversity of populations' survival patterns as caused by predator-taxis and conversion time delay from the perspective of studying Turing bifurcation, Hopf bifurcation, Turing-Turing bifurcation and Turing-Hopf bifurcation for system (1.2).

First, we prove the well-posedness of (1.2). Second, for system (1.2) without conversion time delay, we establish conditions for the existence of Turing bifurcation and Turing-Turing bifurcation. Particularly, we describe the first Turing bifurcation curve precisely in a $(d_u, \alpha)$-plane, which is piecewise smooth with the segment points being Turing-Turing bifurcation points. Hence, the latent steady states with different wave frequencies that system (1.2) may exhibit can be forecasted by virtue of the first Turing bifurcation curve. Compared to ^[22], we give a larger range of parameters when the stability of the positive constant steady state is broken by Turing bifurcation. Moreover, for system (1.2), we theoretically reveal that strong predator-taxis will suppress the appearance of spatial patterns caused by random diffusion.

Afterward, the conditions for the existences of Hopf bifurcation and Hopf-Hopf bifurcation are given under the conditions of considering the effects of conversion time delay. We determine the finite range of wave frequencies for the case of Hopf bifurcation occurring in system (1.2). Especially, by describing the finite Hopf bifurcation curves, we also attain the first Hopf bifurcation curve and determine the stable periodic solutions with different wave frequencies that system (1.2) may exhibit intuitively. It is shown that the positive constant steady state will be destabilized via Hopf-Hopf bifurcation at the segment points. Unlike most previous studies on reaction-diffusion systems with only self-diffusion, our results show that when the positive constant steady state is destabilized, system (1.2) will exhibit stable periodic patterns with non-zero wave frequencies via Hopf bifurcation, as well as some abundant periodic patterns via Hopf-Hopf bifurcation, such as transient quasi-periodic patterns and bistable periodic patterns, i.e., two stable periodic patterns with different spatial wave frequencies coexisting.

Apart from those, the instability of a positive constant steady state may also be brought about by some other bifurcations, such as Turing-Hopf bifurcation, Turing-Turing-Hopf bifurcation, Turing-Hopf-Hopf bifurcation and so on. We also establish the conditions for the existence of the corresponding bifurcations. And it is shown from a numerical perspective that a pair of stable spatially inhomogeneous synchronous time-periodic patterns appear via Turing-Hopf bifurcation.

This paper is organized as follows. Global existence of the classical solutions is given in Section 2. In Sections 3 and 4, we state the main results about Turing instability and Hopf bifurcation for system (1.2). Some numerical simulations are also given to illustrate the theoretical results. At last, we finish our study with conclusions in Section 5.

Throughout the paper, $\mathbb{N}$ is the set of all positive integers, and $\mathbb{N}_{0} = \mathbb{N}\cup\left\{0\right\}$. As for the definitions of the mode-$k_{1}$ Turing bifurcation and mode-$(k_{1}, k_{1}+1)$ Turing-Turing bifurcation, as well as mode-$k_{2}$ Hopf bifurcation, mode-$(k_{1}, k_{2})$ Turing-Hopf bifurcation, mode-$(k_{2}, \tilde{k}_{2})$ Hopf-Hopf bifurcation and so on, which will be mentioned later, the reader may refer to ^[38,39].

2.   Global existence of the classical solutions

Before starting to formally investigate the long-term dynamics, we first prove the well-posedness of system (1.2), and the result is inspired by ^[40].

Theorem 1. System (1.2) admits a unique and positive classical solution $(u, v)\in(C((0, \infty), W^{1, p}(\Omega\times[-\tau, 0], \mathbb{R}^{+}))\cap C^{2, 1}((0, \infty)\times\overline{\Omega}, \mathbb{R}^{+}))^{2}$ with $p > 1$ if the initial values satisfy (1.3).

Proof. When $\tau = 0$, it follows from ^[36] that the conclusion is valid.

When $\tau > 0$, for $0\leqslant t\leqslant\tau$, from the first (second) equation of (1.2) and the comparison principle for parabolic equations, we have that $u\geqslant0$ $(v\geqslant0)$. Further by ^[41], we have that $u > 0\; (v > 0)$ for $0\leqslant t\leqslant\tau$. Moreover, it follows from ^[42] that (1.2) has a unique classical solution $(C((0, T), W^{1, p}(\Omega\times[-\tau, 0], \mathbb{R}^{+}))\cap C^{2, 1}((0, T)\times\overline{\Omega}, \mathbb{R}^{+}))^{2}$. Noting that

the classical solution can be extended to $t\in[0, \tau]$ and $(u, v)\in(C((0, \tau), W^{1, p}(\Omega\times[-\tau, 0], \mathbb{R}^{+}))\cap C^{2, 1}((0, \tau)\times\overline{\Omega}, \mathbb{R}^{+}))^{2}$. Now by repeating the above proof for $t\in[\tau, 2\tau]$ and $t\in[n\tau, (n+1)\tau]\; (n\geqslant2)$, we have the analogous results. Thus, the global existence of the classical positive solution of (1.2) can be obtained, which is obviously positive.

3.   Turing instability and multi-stable spatial patterns

First, we analyze the existence and stability of a positive constant steady state for the homogeneous system corresponding to system (1.2).

In fact, there is a unique positive equilibrium $\bar{E}: = (\bar{u}, \bar{v})$ for system (1.2) if

where

Then, the linearized equations of system (1.2) at $\bar{E}$ are given by

Let $\{\mu_{k} = k^2: k\in\mathbb{N}_{0}\}$ be the eigenvalues of the operator $-\frac{\partial^{2}}{\partial x^{2}}$ on $(0, \pi)$, as subject to Neumann boundary conditions. And for the sake of convenience, denote

then the characteristic equations of (3.2) yield

where

For $\tau = 0$, (3.3) turns into

Denote $DET_{k}: = \sigma_{k}+q_{k}(\alpha)$, $TR_{k}: = -(p_{k}+s_{k})$ for $k\in\mathbb{N}_{0}$; then,

If

then $TR_{0} < 0$; hence, $\bar{E}$ is locally asymptotically stable for a local ODE system since $DET_{0} > 0$. So in the following discussions, we always assume that (3.6) is satisfied.

Next, we devote ourselves to discussing the occurrences of Turing bifurcation for system (1.2).

For convenience, denote

and it is obvious that $\overline{r}_{0} > \underline{r_{0}}$.

When $r_{0}\geqslant\overline{r}_{0}$, then $a\bar{u}-\varpi\geqslant0$; hence, $TR_{k} < 0$ and $DET_{k} > 0$ for $k\in\mathbb{N}_{0}$, which indicates that $\bar{E}$ is locally asymptotically stable.

When $r_{0} < \overline{r}_{0}$, then $a\bar{u}-\varpi < 0$. According to (3.6), $TR_{k} < 0$ for $k\in\mathbb{N}_{0}$ also holds. However, there exists some $k\in\mathbb{N}$ such that $DET_{k} < 0$, meaning that Turing bifurcation can occur for system (1.2). Then, for any fixed $d_{v} > 0$, define

where

Also, denote

Then, to describe the critical conditions for the instability of the positive equilibrium, we first discuss the monotonicity of $\bar{\alpha}(k)$ with respect to the wave number $k$.

Lemma 2. For any $a, b_{1}, b_{2}, m_{1}, d_{v} > 0$, provided that $c > \frac{m_{1}}{b_{1}}$ and $\underline{r_{0}} < r_{0} < \overline{r}_{0}$, it holds that $\bar{\alpha}(k) > 0$ is monotonically increasing in $k$ for $k\in\Lambda$.

A direct calculation from $\alpha(k, d_{u}) = 0$ yields

which is a critical value of $d_{u}$ when $\bar{E}$ is destabilized.

For any $d_{v} > 0$, let

then $\bar{k}$ is the critical wave number of the nonconstant steady states arising from $\bar{E}$ through Turing instability, where

and $\lfloor\cdot\rfloor$ is the floor function.

The next lemma concerns the monotonicity of $d_{u}^{k}$ in $k$ and $\alpha(k, d_{u})$ in $d_{u}$, respectively.

Lemma 3. For any $a, b_{1}, b_{2}, m_{1}, d_{v} > 0$, provided that $c > \frac{m_{1}}{b_{1}}$ and $\underline{r_{0}} < r_{0} < \overline{r}_{0}$, it holds that

(ⅰ) $d_{u}^{k}$ is monotonically increasing with respect to $k\in\mathbb{N}$ for $\hat{k} < k < \bar{k}$, and monotonically decreasing with respect to $k\in\mathbb{N}$ for $k\geqslant\bar{k}$.

(ⅱ) for any $k\in\Lambda$, $\alpha(k, d_{u})$ is linear and monotonically decreasing with respect to $d_{u}$ for $d_{u} > 0$. Particularly, $\alpha(k, d_{u}) > 0$ for $0 < d_{u} < d_{u}^{k}$, $k\in\Lambda$.

Proof. We only prove the case of $d_{u}^{\lfloor k_{*} \rfloor}\leqslant d_{u}^{\lfloor k_{*} \rfloor+1}$, where $\bar{k} = \lfloor k_{*} \rfloor+1$, and the other cases can be proved by following similar arguments.

For part (ⅰ), define $\Theta_{1}(x)$ by

So by a direct calculation, there exists $x^{*} = k_{*}^{2}$ satisfying that $\Theta_{1}(x)$ is increasing in $x$ on $(0, x^{*})$, and decreasing in $x$ on $(x^{*}, \infty)$. And by (3.8), it follows that

that is, $k_{*} > \hat{k}$; thus, $\lfloor k_{*}\rfloor+1\geqslant k_{*} > \hat{k}$, which implies that $d_{u}^{k}$ is decreasing for $\lfloor k_{*}\rfloor+1\leqslant k\in\mathbb{N}$. On the other hand, for $k\in\{k\in\mathbb{N}: \hat{k} < k < \lfloor k_{*}\rfloor+1\}$, $d_{u}^{k}$ is monotonically increasing with respect to $k$.

For part (ⅱ), for any $k\in\Lambda$, we rewrite the expression of $\alpha(k, d_{u})$ in (3.7) as follows

then $\alpha(k, d_{u})$ is a linear function of $d_{u}$; thus, for $d_{u} > 0$, $\alpha(k, d_{u})$ is monotonically decreasing with respect to $d_{u}$, implying that $\alpha(k, d_{u}) > 0$ for $0 < d_{u} < d_{u}^{k}$, $k\in\Lambda$.

As required by the model background, we focus on the case in which the predator-taxis coefficient is non-negative. To this end, for any $d_{v} > 0$, define

It is well known that the first critical value of Turing bifurcation determines the stability of the positive constant steady state ^[43]. So to describe the critical region in which the stability of $\bar{E}$ changes in the $(d_{u}, \alpha)$-plane, we next we discuss the intersections of $\alpha = \alpha_{*}(k, d_u)$ for $k\in\Lambda, 0 < d_{u} < d_{u}^{\bar{k}}$.

Lemma 4. For any $a, b_{1}, b_{2}, m_{1}, d_{v} > 0$, provided that $c > \frac{m_{1}}{b_{1}}$ and $\underline{r_{0}} < r_{0} < \overline{r}_{0}$, it follows that for any $k\in\Lambda$,

(ⅰ) the equation

has a unique root $d_{u}^{k, k+1}\in(0, d_{u}^{k+1})$, which is denoted by

(ⅱ) for $k\geqslant\bar{k}$,

for $d_{u}^{k, k+1} < d_{u} < d_{u}^{k-1, k}$, where $d_{u}^{\bar{k}-1, \bar{k}}: = d_{u}^{\bar{k}}$;

(ⅲ) for $\hat{k} < k < \bar{k}$

Proof. The part (ⅰ) is obvious. As for parts (ⅱ) and (ⅲ), below we only give the proof when $d_{u}^{\lfloor k_{*} \rfloor}\leqslant d_{u}^{\lfloor k_{*} \rfloor+1}$, with $\bar{k} = \lfloor k_{*} \rfloor+1$, and the other cases can be proved similarly.

For part (ⅱ), if $k\geqslant\bar{k}$, then for the first sequence, if $d_{u} > d_{u}^{k, k+1}$, it follows that

that is, $\alpha_{*}(k, d_{u}) > \alpha_{*}(k+1, d_{u})$. In addition, if $d_{u} > d_{u}^{k+1} > d_{u}^{k, k+1}$, together with the definition of $\alpha_{*}(k, d_{u})$ in (3.13), then $\alpha_{*}(k+1, d_{u}) = \alpha_{*}(k+2, d_{u}) = \cdots = 0$, implying that for $d_{u} > d_{u}^{k, k+1}$, $\alpha_{*}(k, d_{u}) > \alpha_{*}(k+1, d_{u})\geqslant\alpha_{*}(k+2, d_{u})\geqslant\cdots$. Similarly, if $d_{u} > d_{u}^{k+1, k+2}$, $\alpha_{*}(k+1, d_{u}) > \alpha_{*}(k+2, d_{u})\geqslant\alpha_{*}(k+3, d_{u})\geqslant\cdots$ holds, and the rest can be done in the same manner. And by (3.14), $d_{u}^{\lfloor k_{*} \rfloor+1} > d_{u}^{k, k+1} > d_{u}^{k+1, k+2} > \cdots > 0$ also holds. Thus, when $d_{u} > d_{u}^{k, k+1}$, $\alpha_{*}(k, d_u) > \alpha_{*}(k+1, d_u)\geqslant\alpha_{*}(k+2, d_u)\geqslant\cdots$.

If $d_{u} < d_{u}^{k-1, k}$, it analogously follows that

that is, $\alpha_{*}(k, d_{u}) > \alpha_{*}(k-1, d_{u})$. So if $d_{u} < d_{u}^{k-2, k-1}$, then $\alpha_{*}(k-1, d_{u}) > \alpha_{*}(k-2, d_{u})$ and so on. By (3.14) again, $0 < d_{u}^{k-1, k} < d_{u}^{k-2, k-1} < \cdots < d_{u}^{\lfloor k_{*} \rfloor+1}$ holds. Hence, when $d_{u} < d_{u}^{k-1, k}$, $\alpha_{*}(k, d_u) > \alpha_{*}(k-1, d_u) > \alpha_{*}(k-2, d_u) > \cdots > \alpha_{*}(\lfloor k_{*} \rfloor+1, d_u)$.

For part (ⅲ), according to Lemma 2 and Lemma 3, if $\hat{k} < k < \bar{k}$, then $\alpha_{*}(\lfloor k_{*} \rfloor+1, d_{u})\geqslant\alpha_{*}(\lfloor k_{*} \rfloor, d_{u})\geqslant\cdots\geqslant\alpha_{*}(\lfloor \hat{k} \rfloor+1, d_{u})$, $0 < d_{u} < d_{u}^{\bar{k}}$.

Synthesizing the above discussions, for $k\in\Lambda$ and $0 < d_{u} < d_{u}^{\bar{k}}$, define

with $\bar{k}$ as defined in (3.11). Then

is called the first Turing bifurcation curve, which is the critical curve for the Turing instability of $\bar{E}$.

Lemma 5. For any $a, b_{1}, b_{2}, m_{1}, d_{v} > 0$, provided that $c > \frac{m_{1}}{b_{1}}$ and $\underline{r_{0}} < r_{0} < \overline{r}_{0}$, it holds that, for any $d_{u}\in(0, d_{u}^{\bar{k}})$, there must exist some integer $k_{1}\geqslant\bar{k}$ such that $d_{u}\in[d_{u}^{k_{1}, k_{1}+1}, d_{u}^{k_{1}-1, k_{1}})$. Specifically,

(ⅰ) if $d_{u}\in(d_{u}^{k_{1}, k_{1}+1}, d_{u}^{k_{1}-1, k_{1}})$ and $\alpha = \alpha_{*}(d_{u})$, $0$ is a simple root of characteristic equation (3.5) with $k = k_{1}$, and the other roots of (3.5) have strictly negative real parts; also, let $\lambda_{1}(k_{1}, 0, \alpha)$ be the root of (3.5) satisfying $\lambda_{1}(k_{1}, 0, \alpha_{*}(d_{u})) = 0$; then,

(ⅱ) if $d_{u} = d_{u}^{k_{1}, k_{1}+1}$ and $\alpha = \alpha_{*}(d_{u}^{k_{1}, k_{1}+1})$, then $0$ is a simple root of characteristic equation (3.5) for both $k_{1}$ and $k_{1}+1$, and the other roots of (3.5) have strictly negative real parts.

Proof. Inspired by the proof of ^[43], we only give the proof of part (ⅰ), as the second assertion can be analogously proved.

Note that in (3.5), $DET_{k} = 0$ if and only if $\alpha = \alpha_{*}(k, d_{u})$ for $k\in\Lambda$ and $0 < d_u < d_{u}^{\bar{k}}$. So, for any $k_{1}\in\Lambda$, $\lambda_{1}(k_1, 0, \alpha) = 0$ is always a root of (3.5) with such a $k_{1}$ when $\alpha = \alpha_{*}(k_{1}, d_{u})$. By the definitions of $\alpha_{*}(d_{u})$ and $\alpha_{*}(k, d_{u})$, if $d_{u}\in(d_{u}^{k_{1}, k_{1}+1}, d_{u}^{k_{1}-1, k_{1}})$ and $\alpha = \alpha_{*}(d_{u})$, then $0$ is a root of (3.5) with $k = k_{1}$. Moreover, it follows from

that $\lambda = 0$ is simple. The condition (3.6) ensures that $TR_{k} < 0$ for all $k\in\mathbb{N}_{0}$ and $DET_{k} > 0$ for $k\geqslant\bar{k}, k\in\mathbb{N}$ and $k\neq k_{1}$. Thus, all other roots of (3.5) have strictly negative real parts. Differentiating (3.5) with respect to $\alpha$, the transversality conditions in (3.19) can be obtained directly.

In Figure 1, we show the schematic diagram of Turing bifurcation curves for $\alpha = \alpha_{*}(k, d_{u})$ when $0 < d_{u} < d_{u}^{\bar{k}}$ for different values of $k\in\Lambda$ in the $(d_{u}, \alpha)$-plane to illustrate the properties presented in Lemma 4. The corresponding first Turing bifurcation curve $\alpha = \alpha_{*}(d_{u}), 0 < d_{u} < d_{u}^{\bar{k}}$ is plotted in Figure 1, and the non-smooth points of $\alpha = \alpha_{*}(d_{u}), 0 < d_{u} < d_{u}^{\bar{k}}$, $T_{\bar{k}, \bar{k}+1}$, $T_{\bar{k}+1, \bar{k}+2}$, $\cdots$ are Turing-Turing bifurcation points.

So in summary, we have the following theorem.

Theorem 6. For model (1.2) with $\tau = 0$ and any $a, b_{1}, b_{2}, m_{1}, d_{v} > 0$, provided that $c > \frac{m_{1}}{b_{1}}$, it follows that

(1) when $r_{0}\geqslant\overline{r}_{0}$, for $d_{u} > 0$ and $\alpha > 0$, the positive constant steady state $\bar{E}$ is always locally asymptotically stable;

(2) when $\underline{r_{0}} < r_{0} < \overline{r}_{0}$,

$(a)$ if $d_{u}\geqslant d_{u}^{\bar{k}}$, $\bar{E}$ is locally asymptotically stable for $\alpha > 0$;

$(b)$ when $0 < d_{u} < d_{u}^{\bar{k}}$,

(ⅰ) $\bar{E}$ is locally asymptotically stable for $\alpha > \alpha_{*}(d_{u})$, and unstable for $0 < \alpha < \alpha_{*}(d_{u})$;

(ⅱ) if $d_{u}\in(d_{u}^{k_{1}, k_{1}+1}, d_{u}^{k_{1}-1, k_{1}})$ for some $k_{1}\geqslant\bar{k}$ and $k_{1}\in\mathbb{N}$, (1.2) undergoes mode-$k_1$ Turing bifurcation when $\alpha = \alpha_{*}(d_{u})$;

(ⅲ) if $d_{u} = d_{u}^{k_{1}, k_{1}+1}$ for some $k_{1}\geqslant \bar{k}$ and $k_{1}\in\mathbb{N}$, (1.2) undergoes mode-$(k_{1}, k_{1}+1)$ Turing-Turing bifurcation when $\alpha = \alpha_{*}(d_{u}^{k_{1}, k_{1}+1})$.

Remark 1. (Turing patterns and Turing-Turing patterns) Note that sufficient conditions for the formations of spatial patterns have been provided in ^[22], while the condition that $0 < \alpha < \alpha_{*}(d_{u}), \; 0 < d_{u} < d_{u}^{\bar{k}}$ given in Theorem 6 is not only sufficient but also necessary for pattern formations.

As shown in Figure 1$(b)$, the first Turing bifurcation curve $\alpha = \alpha_{*}(d_{u}), \; 0 < d_{u} < d_{u}^{\bar{k}}$, is formed by connecting Turing bifurcation curves of mode-$\bar{k}$, mode-$\bar{k}+1$, mode-$\bar{k}+2$, $\cdots$ with Turing-Turing bifurcation points $T_{\bar{k}, \bar{k}+1}$, $T_{\bar{k}+1, \bar{k}+2}$, $\cdots$, which indicates that the positive constant steady state can be destabilized by mode-$k$ Turing bifurcation; thus, system (1.2) will harvest the spatially inhomogeneous steady states shaped like $\cos kx$, and it also suggests that the positive constant steady state can be destabilized by a mode-$(k, k+1)$ Turing-Turing bifurcation; thus, system (1.2) will harvest more complex spatial patterns, such as multiple superposition solutions of $\cos kx$ and $\cos (k+1)x$, or the coexistence of multi-stable spatial patterns shaped like $\cos kx$ and $\cos (k+1)x$, where the possible values of $k$ are $\bar{k}, \bar{k}+1, \cdots$.

Moreover, it can be concluded that strong predator-taxis has the effects of eliminating spatial patterns which arise due to the random dispersal of predator and prey. Besides, it follows from Lemma 3 and (3.17) that $\alpha = \alpha_{*}(d_{u})$ is decreasing with respect to $d_{u}$ for $0 < d_{u} < d_{u}^{\bar{k}}$, which means that when the self-diffusion $d_u$ of prey is small, the heterogeneity of the spatial distributions of the two populations can be eliminated by increasing the predator-taxis coefficient $\alpha$, that is to say, strong predator-taxis can supplement weak self-diffusion of prey.

Next, we provide some examples to intuitively explain the above theoretical analysis.

Here we refer to the parameter selections in ^[44] and let

Then (3.6) is satisfied; thus, $\bar{E} = (0.1877, 0.1276)$ is locally asymptotically stable for the local ODE system.

Further, when $d_v = 0.2, \; 1.8$, respectively, the Turing bifurcation sets are as given in Figure 2, where $\alpha = \alpha_{*}(k, d_{u})$ for $k = 1, 2, 3, 4, \cdots$ are Turing bifurcation curves and the regions in which multiple spatial patterns exist are marked. For ease of citation, stable regions are denoted as $D$, the mode-$2$ Turing bifurcation region is denoted as $D_{2}$ and the mode-$(1, 2)$ and mode-$(2, 3)$ Turing-Turing bifurcation regions are denoted as $D_{1, 2}$ and $D_{2, 3}$, respectively.

Example 1. (Turing bifurcation and spatial patterns)

Let $d_{v} = 0.2$; we want to show the effects of predator-taxis on pattern formations.

It follows from (3.11) that $\bar{k} = 2$. Due to (3.10), $d_{u}^{\bar{k}} = d_{u}^{2} = 0.0241$. And by (3.14), $d_{u}^{2, 3} = 0.0031$.

When $d_{u}\in(d_{u}^{2, 3}, d_{u}^{2})$, Theorem 6 indicates that system (1.2) undergoes a mode-$2$ Turing bifurcation at $\alpha = \alpha_{*}(d_{u})$. And the eigenfunction corresponding to the eigenvalue $\mu_{2}$ is

Bifurcation theory indicates that when $(d_u, \alpha)\in D_{2}$, there will be spatially inhomogeneous patterns and the waveforms are consistent with those of eigenfunctions given by (3.20).

(1) Let $d_{u} = 0.02$; then, $\alpha_{*}(0.02) = 0.1384$. Let $\alpha = 0 < \alpha_{*}(0.02)$ and $\alpha = 0.087 < \alpha_{*}(0.02)$ respectively; then, $(d_u, \alpha)\in D_{2}$ (see Figure 2 $(a)$). Thus, a mode-$2$ Turing bifurcation occurs and $\bar{E}$ becomes unstable. Meanwhile, spatially inhomogeneous patterns shaped like $\cos2x$ are theoretically expected to appear. Specifically, when choosing $(0.1877+0.1\cos2x, 0.1276+0.1\cos2x)$ as initial values, the numerical results are as presented in Figure 3, where the third column of Figure 3 illustrates that the waveforms of spatially inhomogeneous patterns are consistent with the waveforms of the corresponding eigenfunctions given by (3.20).

(2) If we further choose a larger predator-taxis coefficient $\alpha = 0.5 > \alpha_{*}(0.02) = 0.1384$, then $(d_u, \alpha)\in D$ (see Figure 2 $(a)$). Thus $\bar{E}$ is locally asymptotically stable. Let $(0.1877+0.1\cos2x, 0.1276+0.1\cos2x)$ be initial values, and the numerical results are as presented in Figure 4.

Consequently, it is numerically confirmed by comparing Figure 3 with Figure 4 that strong predator-taxis contributes to eliminating spatial patterns resulting from self-diffusion.

Remark 2. Actually, when the system (1.2) undergoes a mode-$k$ Turing bifurcation, the appearance of the spatially inhomogeneous patterns shaped like $\cos kx$ is independent of the choice of initial values. By taking the mode-$2$ Turing bifurcation of the system as an example, we additionally perform several sets of numerical examples for different initial values, with the other parameters chosen to be the same as in Figure 3 $(b)$. As shown in Figure 5, we specially plot the shapes of the initial values $(u_{0}, v_{0})$, indicating that the solutions of the system (1.2) will tend to the same spatial patterns regardless of the selection of initial values.

Example 2. (Turing-Turing bifurcation and multi-stable spatial patterns)

(1) Let $d_{v} = 0.2$. Compared with Example 1(2), we further choose a smaller self-diffusion coefficient of prey $d_{u} = 0.0055$. It follows from (3.10), (3.13) and (3.14) that

thus, $d_{u}\in(d_{u}^{2, 3}, d_{u}^{3})$ and $(d_u, \alpha)\in D_{2, 3}$, that is, $(d_u, \alpha)$ is selected in the mode-$(2, 3)$ Turing-Turing bifurcation region (see Figure 2 $(a)$), where superposition patterns of $\cos2x$ and $\cos3x$ or multi-stable spatial patterns shaped like $\cos2x$ and $\cos3x$ are theoretically expected to emerge. When we choose $(0.1877+0.1\cos2x, 0.1276+0.1\cos2x)$, $(0.1877-0.1\cos2x, 0.1276-0.1\cos2x)$, $(0.1877+0.1\cos3x, 0.1276+0.1\cos3x)$ and $(0.1877-0.1\cos3x,$ $0.1276-0.1\cos3x)$ as initial values respectively, as shown in Figure 6, a pair of spatially inhomogeneous patterns with wave frequency 2 and a pair of spatially inhomogeneous patterns with wave frequency 3 coexist, which are multi-stable spatial patterns.

(2) Let $d_{u} = 0.02$. Compared with Example 1(2), we further choose a larger self-diffusion coefficient for the predator, $d_{v} = 1.8$. It follows from (3.10), (3.13) and (3.14) that

thus $d_{u}\in(d_{u}^{1, 2}, d_{u}^{2})$ and $(d_u, \alpha)\in D_{1, 2}$, that is, $(d_u, \alpha)$ is selected in the mode-$(1, 2)$ Turing-Turing bifurcation region (see Figure 2 $(b)$), where superposition patterns of $\cos x$ and $\cos2x$ or multi-stable spatial patterns shaped like $\cos x$ and $\cos2x$ are theoretically expected to emerge. When we choose $(0.1877+0.1\cos x, 0.1276+0.1\cos x)$, $(0.1877-0.1\cos x,$ $0.1276-0.1\cos x)$, $(0.1877+0.2\cos2x,$ $0.1276+0.1\cos2x)$ and $(0.1877-0.2\cos2x, 0.1276-0.1\cos2x)$ as initial values respectively, as shown in Figure 7, a pair of spatially inhomogeneous patterns with wave frequency $1$ and a pair of spatially inhomogeneous patterns with wave frequency $2$ coexist, which are multi-stable spatial patterns.

Consequently, the results of Example 2 show that, on the one hand, the positive constant steady state can be destabilized via Turing-Turing bifurcation, resulting in multi-stable spatial patterns. On the other hand, by comparing the results of Figure 6 and Figure 7 with those of Figure 4, respectively, it is numerically illustrated that either minimal self-diffusion of the prey or a large amount of self-diffusion of the predator will encourage spatial patterns to appear, which suggests that predator-taxis can balance the spatial heterogeneity caused by self-diffusion.

For convenience, we summarize the results of Example 1 and Example 2 in Table 1 and make the following comparisons.

The comparison of cases Ⅰ and Ⅱ shows that strong predator-taxis can eliminate spatial patterns caused by self-diffusion, that is to say, the more sensitive the prey to predation (i.e., strong predator-taxis), the more evenly distributed are the two populations.

Beyond that, a reasonable phenomenon is explained, that is, when only self-diffusion is considered, fast predator dispersal leads to the appearance of spatial patterns, which will eventually disappear when further considering the chemotaxis behavior of prey avoiding predator, suggesting that predator-taxis can partially cancel out the non-uniform advantage caused by self-diffusion of the predator.

The comparison of cases Ⅱ and Ⅲ shows that spatial patterns will exist as long as the value of self-diffusion of the prey is small enough. From another perspective, it also shows that predator-taxis and self-diffusion of the prey are complementary, that is, predator-taxis is equivalent to accelerating the self-diffusion of the prey, and the spatial distributions of the two populations will finally reach homogeneity if predator-taxis is strong enough.

The comparison of cases Ⅱ and Ⅳ shows that the spatial patterns will exist as long as the value of self-diffusion of the predator is large enough. From another perspective, the accelerated self-diffusion of predator is equivalent to the weakened ability of the prey to avoid the risk of being targeted, thus increasing the probability of emergence of spatial inhomogeneous patterns.

Case Ⅲ and Case Ⅳ show that spatially inhomogeneous steady states with different wave frequencies can coexist, which in turn indicates that the constant stationary can be destabilized by Turing-Turing bifurcation, thus giving rise to multi-stable spatial patterns.

4.   Hopf bifurcation and spatiotemporal staggered periodic patterns

In this part, we discuss the effects of time delay on the stability of the positive constant steady state of system (1.2).

First, the following result concerning the nonoccurence of Hopf bifurcation for system (1.2) with $\tau = 0$ holds.

Theorem 7. For model (1.2) with $\tau = 0$ and $a, b_{1}, b_{2}, m_{1}, d_{v} > 0$, provided that $c > \frac{m_{1}}{b_{1}}$ and $r_{0} > \underline{r_{0}}$, it follows that, for $d_{u}, \alpha > 0$, there is no Hopf bifurcation.

Next, we establish the conditions for the occurrence of Hopf bifurcation for system (1.2) with $\tau > 0$.

Let $\lambda = \pm \mathrm{i}\omega_{k}(\alpha)$ with $\omega_{k}(\alpha) > 0$ be a pair of purely imaginary roots of (3.3). For the sake of convenience, denote $\omega_{k}(\alpha)\triangleq\omega_{k}$. Then, for $k\in\mathbb{N}_{0}$,

Separating the real and imaginary parts yields

which implies that

Denote

So, we first discuss the sign of $\sigma_{k}^{2}-q_{k}(\alpha)^{2} = (\sigma_{k}+q_{k}(\alpha))(\sigma_{k}-q_{k}(\alpha))$. Since $\sigma_{k}+q_{k}(\alpha) > 0$ is guaranteed by Theorem 7, the sign of $\sigma_{k}^{2}-q_{k}(\alpha)^{2}$ coincides with that of $\sigma_{k}-q_{k}(\alpha)$, and

Because $\sigma_{k}-q_{k}(\alpha) = 0$ with respect to $k$ has only one positive root, denoted as $K^{0}$ for simplicity, where

with $\Gamma: = d_{v}(a\bar{u}-\varpi)-d_{u}cb_{2}\varpi-\alpha cb_{2}\bar{v}\varpi$, we can conclude that

And, the following theorem deals with the case concerning the roots of (4.3).

Theorem 8. For $a, b_{1}, b_{2}, m_{1}, d_{v} > 0$, provided that $c > \frac{m_{1}}{b_{1}}$ and $r_{0} > \underline{r_{0}}$, for $d_{u}, \alpha > 0$ and $k\in\mathbb{N}_{0}$, it follows that

(1) when $0 \leqslant k < K^{0}$, $\omega_{k}^{+}$ is the unique positive root of (4.3);

(2) when $k \geqslant K^{0}$, (4.3) has no positive root.

Proof. For $0 \leqslant k < K^{0}$, since $\sigma_{k}^{2}-q_{k}^{2}(\alpha) < 0$, $\omega_{k}^{+}$ is the only positive root of (4.3), regardless of the positivity or negativity of $s_{k}^{2}-p_{k}^{2}+2\sigma_{k}$.

For $k\geqslant K^{0}$, $\sigma_{k}^{2}-q_{k}^{2}(\alpha)\geqslant0$ holds, and it follows from (4.5) that for $k = K^{0}$,

Note that

thus,

$(a)$ if $\underline{r_{0}} < r_{0} < \overline{r}_{0}$, that is to say, $a\bar{u}-\varpi < 0$, one can show that $s_{k}^{2}-p_{k}^{2}+2\sigma_{k} = 0$ with respect to $k$ has a unique positive root, denoted as $K^{+}$. And by (4.6), we have

so it can be verified that

which yields that $s_{K^{0}}^{2}-p_{K^{0}}^{2}+2\sigma_{K^{0}} < 0$, that is to say, $s_{k}^{2}-p_{k}^{2}+2\sigma_{k} < 0$ for $k\geqslant K^{0} > K^{+}$. Hence, from (4.4), we know that (4.3) has no positive root for $k\geqslant K^{0}$;

$(b)$ if $r_{0}\geqslant\overline{r}_{0}$, that is to say, $a\bar{u}-\varpi\geqslant0$, then $d_{u}(K^{0})^2+(a\bar{u}-\varpi) > 0$. It follows from (4.7) that $d_{v}(K^{0})^2-cb_{2}\varpi > 0$. By a similar argument to $(a)$, (4.3) has no positive root for $k\geqslant K^{0}$ either.

Thus, due to Theorem 8, if the positive root $\omega_{k} = \omega_{k}^{+}$ of (4.3) exists, denote the root of (4.2) in $(0, 2\pi]$ as $\tau_{k}(\alpha)$. For $a, b_{1}, b_{2}, m_{1}, d_{v} > 0$, provided that $c > \frac{m_{1}}{b_{1}}$ and $r_{0} > \underline{r_{0}}$, given that $d_{u}, \alpha > 0$, denote the critical values for $\tau$ by

Suppose that $\lambda_{2}(k, \tau, \alpha) = \xi(k, \tau, \alpha)\pm \mathrm{i}\omega(k, \tau, \alpha)$ are a pair of conjugated complex roots of the characteristic equation $D_{k}(\lambda, \tau, \alpha) = 0$ near $\tau = \tau_{k}^{(j)}(\alpha)$ with $\xi(k, \tau_{k}^{(j)}(\alpha), \alpha) = 0$ and $\omega(k, \tau_{k}^{(j)}(\alpha), \alpha) = \omega_{k}^{+} > 0$. According to ^[45,46], we obtain the conclusion regarding the transversality conditions.

Theorem 9. For $a, b_{1}, b_{2}, m_{1}, d_{v} > 0$, provided that $c > \frac{m_{1}}{b_{1}}$ and $r_{0} > \underline{r_{0}}$, it follows that for $d_{u}, \alpha > 0$,

Proof. For $0\leqslant k < K^{0}$ and $k\in\mathbb{N}_{0}$, denote

Obviously, $P(\mathrm{i}\omega_{k}^{+})+Q(\mathrm{i}\omega_{k}^{+})\neq0$. For $\forall\omega_{k} > 0$, define

Hence, $F(\omega_{k}) = 0$ implies that (4.3) holds, and its roots are given by (4.4). If we denote

then for $\omega_{k} = \omega_{k}^{+}$, it follows that

By (4.12) and (4.13), we have

By (4.2), define

with

It follows from (4.11) that

which is independent of $\tau$. Therefore,

Then, by ^[45], we have

Thus, we have completed the proof.

Finally, we have the conclusion on the Hopf bifurcation of system (1.2).

Theorem 10. For $a, b_{1}, b_{2}, m_{1}, d_{v} > 0$, provided that $c > \frac{m_{1}}{b_{1}}$ and $r_{0} > \underline{r_{0}}$, it follows that, for $d_{u}, \alpha > 0$, there exists $k_{2}\in[0, K^{0})$ such that

So,

(1) when $\tau = \tau_{k_{2}}(\alpha)$, $D_{k_{2}}(\lambda, \tau, \alpha) = 0$ has a pair of purely imaginary roots for $\lambda$ and all other roots of (3.3) have strictly negative real parts for $k\in\mathbb{N}_{0}$. When $\tau\in\left[0, \tau_{k_{2}}(\alpha)\right)$, the steady state $\bar{E}$ is locally asymptotically stable;

(2) if $k_{2}$ is unique, system (1.2) undergoes a mode-$k_{2}$ Hopf bifurcation when $\tau = \tau_{k_{2}}(\alpha)$; and if there are exactly $k_{2}$ and $\tilde{k}_{2}$ satisfying (4.16), with $k_{2}, \tilde{k}_{2}\in[0, K^{0})$ and $k_{2}\neq\tilde{k}_{2}$, system (1.2) will undergo a mode-$(k_{2}, \tilde{k}_{2})$ Hopf-Hopf bifurcation when $(\tau, \alpha) = (\tau_{k_{2}}(\alpha_{H}), \alpha_{H})$, where $\alpha_{H}$ denotes the root of $\tau_{k_{2}}(\alpha) = \tau_{\tilde{k}_{2}}(\alpha)$.

Remark 3. (Hopf patterns and Hopf-Hopf patterns) In Theorem 8, the range of wave numbers in which system (1.2) undergoes a Hopf bifurcation is determined, which is upper bounded by $K^{0}$.

Hopf bifurcation sets are given in Figure 8 $(a)$, where some Hopf bifurcation curves for different modes are plotted and the regions in which multiple periodic patterns exist are marked, including Hopf bifurcation regions $D_{0}^{H}, \; D_{1}^{H}, \; D_{2}^{H}$ where stable spatially homogeneous or inhomogeneous periodic patterns are expected to appear, and Hopf-Hopf bifurcation regions $D_{0, 1}^{H}, \; D_{1, 2}^{H}, \; D_{2, 3}^{H}$ where periodic patterns with different spatial wave frequencies coexist or quasi-periodic patterns are anticipated to emerge.

And the corresponding first Hopf bifurcation curve is plotted in Figure 8 $(b)$, which was formed by connecting Hopf bifurcation curves of mode-$0$, mode-$1$, mode-$2$, $\cdots$, mode-$\lfloor K^{0}\rfloor$ with Hopf-Hopf bifurcation points $H_{0, 1}$, $H_{1, 2}$, $\dots$, $H_{\lfloor K^{0}\rfloor-1, \lfloor K^{0}\rfloor}$, indicating that the positive constant steady state can be destabilized by mode-$k$ Hopf bifurcation; thus, system (1.2) will harvest the spatially homogeneous or inhomogeneous staggered periodic solutions shaped like $\cos kx\cos \omega t$, also suggesting that the positive constant steady state can be destabilized by mode-$(k, k+1)$ Hopf-Hopf bifurcation. Thus, system (1.2) will harvest more complex spatiotemporal patterns, such as the coexistence of periodic solutions with two different spatial wave frequencies, or some quasi-periodic solutions, where $\omega$ is a positive constant and $0\leqslant k < K^{0}, k\in\mathbb{N}_{0}$. This is completely different from systems with only self-diffusion ^[37,47], where the positive constant steady state can only be destabilized by a mode-$0$ Hopf bifurcation, producing a stable spatially homogeneous periodic solution.

Combining Theorems 6 and 10, it is found that predator-taxis $\alpha$ and delay $\tau$ can induce other high codimensional bifurcations under certain conditions, such as Turing-Hopf bifurcation.

Theorem 11. For $a, b_{1}, b_{2}, m_{1}, d_{v} > 0$, provided that $c > \frac{m_{1}}{b_{1}}$ and $\underline{r_{0}} < r_{0} < \overline{r}_{0}$, it follows that,

(1) if $d_{u}\in(d_{u}^{k_{1}, k_{1}+1}, d_{u}^{k_{1}-1, k_{1}})$ for some $k_{1}\geqslant\bar{k}$ and $k_{1}\in\mathbb{N}$, system (1.2) undergoes a mode-$(k_{1}, k_{2})$ Turing-Hopf bifurcation near $\bar{E}$ when $(\tau, \alpha) = (\tau_{k_{2}}(\alpha_{*}(d_{u})), \alpha_{*}(d_{u}))$, where $\bar{k}$ is defined as in (3.11) and $k_{2}$ is uniquely determined by (4.16);

(2) if $d_{u}\in(d_{u}^{k_{1}, k_{1}+1}, d_{u}^{k_{1}-1, k_{1}})$ for some $k_{1}\geqslant\bar{k}$ and $k_{1}\in\mathbb{N}$, system (1.2) will undergo a mode-$(k_{1}, k_{2}, \tilde{k}_{2})$ Turing-Hopf-Hopf bifurcation near $\bar{E}$ when $(d_{u}, \tau, \alpha) = (d_{u}^{H}, \tau_{k_{2}}(\alpha_{H}), \alpha_{H})$, where $d_{u}^{H}$ denotes the root of $\alpha_{*}(d_{u}) = \alpha_{H}$, $\alpha_{H}$ is defined in Theorem 10, $\bar{k}$ is defined as in (3.11) and $k_{2}$, $\tilde{k}_{2}$ are determined by Theorem 10;

(3) if $d_{u} = d_{u}^{k_{1}, k_{1}+1}$ for some $k_{1}\geqslant\bar{k}$ and $k_{1}\in\mathbb{N}$, system (1.2) will undergo a mode-$(k_{1}, k_{1}+1, k_{2})$ Turing-Turing-Hopf bifurcation near $\bar{E}$ when $(d_{u}, \tau, \alpha) = (d_{u}^{k_{1}, k_{1}+1},$ $\tau_{k_{2}}(\alpha_{*}(d_{u}^{k_{1}, k_{1}+1})), \alpha_{*}(d_{u}^{k_{1}, k_{1}+1}))$, where $\bar{k}$ is defined as in (3.11) and $k_{2}$ is determined by (4.16).

Next, we illustrate some examples to support and extend our analytical results.

With reference to the parameter selections in ^[22], let

Then, (3.6) is satisfied and $\bar{E} = (0.2412, 0.2533)$ is locally asymptotically stable for the local ODE system.

In the following discussions, we observe that the spatiotemporal distribution types of the two populations are always analogous. Therefore, we only take the prey distribution as examples to show distribution patterns. Besides, the distribution patterns in $(u, x, t)$-plane and $(x, t)$-plane are presented separately in each of the following figures.

Example 3. (Hopf bifurcation and spatially inhomogeneous staggered periodic patterns, Hopf-Hopf bifurcation and bistable periodic patterns)

Let

It follows from (3.11) that $\bar{k} = 3$ and $d_{u} > d_{u}^{\bar{k}} = 0.0284$; then, due to Theorem 6, for $\alpha > 0$, the positive steady state $\bar{E}$ is locally asymptotically stable. According to Theorem 10, Hopf bifurcation curves in the $(\alpha, \tau)$-plane can be characterized for different modes as in Figure 8. And,

When

respectively, either $\bar{E}$ is locally asymptotically stable or system (1.2) may accordingly generate stable spatially homogeneous periodic patterns or spatially inhomogeneous periodic patterns with spatial wave frequency $1$ or $2$, which are bifurcated from $\bar{E}$ through mode-$0$, mode-$1$ and mode-$2$ Hopf bifurcation respectively.

Further, when

respectively, in addition to the above mentioned periodic patterns, system (1.2) may also generate quasi-periodic patterns or the coexisting periodic patterns with different spatial wave frequencies, which are bifurcated from $\bar{E}$ through mode-$(0, 1)$ and mode-$(1, 2)$ Hopf-Hopf bifurcation, respectively. Specifically speaking, suppose we have the following.

(1) Let $(\alpha, \tau) = (1.2, 0.72)$. By (4.9), we have that $K^{0} = 5.1978$ and

Thus, $(\alpha, \tau)\in D_{2}^{H}$. It follows from Theorem 10 that system (1.2) undergoes a mode-$2$ Hopf bifurcation when $\tau = \tau_{2}^{(0)}(1.2) = 0.7004$ and spatially inhomogeneous staggered periodic patterns shaped like $\cos2x\cos \omega t$ are theoretically expected to appear, where $\omega$ is a positive constant. The numerical results are shown in Figure 9.

Let $(\alpha, \tau) = (0.7, 0.87)$; then, by a similar process, $(\alpha, \tau)\in D_{1}^{H}$; thus, system (1.2) undergoes a mode-$1$ Hopf bifurcation when $\tau = \tau_{1}^{(0)}(0.7) = 0.8607$ and spatially inhomogeneous staggered periodic patterns shaped like $\cos x\cos \omega t$ are theoretically expected to appear, where $\omega$ is a positive constant; see Figure 10.

(2) Let $(\alpha, \tau) = (0.61, 0.91364)$. By (4.9), we have that $K^{0} = 4.4040$, and

Thus $(\alpha, \tau)\in D_{0, 1}^{H}$. Then bistable spatiotemporal periodic patterns shaped like $\cos \omega t$ and $\cos x\cos \omega t$ or their superposition are theoretically expected to appear, where $\omega$ is a positive constant. When choosing $(0.2412, 0.2533)$ and $(0.2412+0.01\cos x, 0.2533+0.01\cos x)$ as initial values respectively, as shown in the first two columns of Figure 11, spatially homogeneous and inhomogeneous staggered periodic patterns coexist. When choosing $(0.2412, 0.2533+0.01\cos x)$ as an initial value, transient quasi-periodic patterns are as presented in the third column of Figure 11. Hence, these phenomena are bistable periodic patterns.

Let $(\alpha, \tau) = (0.87, 0.8415)$; then, by a similar process, we know that $(\alpha, \tau)\in D_{1, 2}^{H}$. Then spatiotemporal periodic patterns shaped like $\cos x\cos \omega t$ and $\cos 2x\cos \omega t$ or their superposition are theoretically expected to appear, where $\omega$ is a positive constant. When choosing $(0.2412+0.01\cos x, 0.2533+0.01\cos x)$ and $(0.2412+0.01\cos 2x, 0.2533+0.01\cos 2x)$ as initial values respectively, as shown in the first two columns of Figure 12, spatially inhomogeneous staggered periodic patterns with wave frequencies $1$ and $2$ coexist. When further choosing $(0.2412+0.01\cos x, 0.2533+0.01\cos 2x)$ as an initial value, transient quasi-periodic patterns are as presented in the third column of Figure 12. That is, bistable periodic patterns are also found, but the spatial wave frequencies are different from those shown in Figure 11.

Example 4. (Turing-Hopf bifurcation and bistable spatially inhomogeneous synchronous periodic patterns)

Let

then $d_{u}\in(d_{u}^{3, 4}, d_{u}^{1, 2}) = (0.0084, 0.0334)$, and by Theorems 6 and 10, we conclude that $k_{1} = 3$ and $k_{2} = 0$. It follows from Theorem 11 that system (1.2) undergoes a mode-$(3, 0)$ Turing-Hopf bifurcation when $(\tau, \alpha) = (\tau_{0}(\alpha_{*}(0.113)), \alpha_{*}(0.113)) = (0.887,$ $0.1712)$.

We show the corresponding bifurcation sets for system (1.2) in Figure 13 $(a)$, and the intersection of the mode-$3$ Turing bifurcation curve $\alpha = \alpha_{*}(0.0113) = 0.1712$ and mode-$0$ Hopf bifurcation curve $\tau = \tau_{*} = 0.8870$ is denoted by $TH_{0}$. We have marked the Turing-Hopf bifurcation region and Turing bifurcation region as $D_{TH_{0}}$ and $D_{T}$, respectively.

In particular, we have plotted the corresponding first bifurcation curve in Figure 13 $(b)$, which was formed by connecting the first Turing bifurcation curve $\alpha = \alpha_{*}(d_{u}), \; 0 < d_{u} < d_{u}^{\bar{k}}$, Hopf bifurcation curves of mode-$0$, mode-$1$, mode-$2$, $\dots$, mode-$\lfloor K^{0}\rfloor$ with Turing-Hopf bifurcation point $TH_{0}$ and Hopf-Hopf bifurcation points $H_{0, 1}$, $H_{1, 2}$, $\dots$, $H_{\lfloor K^{0}\rfloor-1, \lfloor K^{0}\rfloor}$, indicating that the positive constant steady state can not only be destabilized by a mode-$k_{2}$ Hopf bifurcation or mode-$(k_{2}, k_{2}+1)$ Hopf-Hopf bifurcation, but also by a mode-$(k_{1}, 0)$ Turing-Hopf bifurcation, where system (1.2) will harvest the spatially inhomogeneous synchronous periodic solutions shaped like $h_{1}\cos k_{1}x\pm h_{2}\cos \omega t$, where $h_{1}$ and $h_{2}$ are constants, $\omega$ is a positive constant, $k_{1} = \bar{k}, \bar{k}+1, \cdots$ and $0\leqslant k_{2} < K^{0}, k_{2}\in\mathbb{N}_{0}$.

Further, let $(\alpha, \tau) = (0.06, 0.9)\in D_{TH_{0}}$, and we find a pair of stably coexisting spatially inhomogeneous synchronous periodic patterns, as shown in Figure 14. Unlike in Figure 9 and Figure 10, the spatial non-homogeneity is caused by the occurrence of Turing bifurcation in system (1.2). Actually, Turing-Hopf bifurcation can also reveal other spatiotemporal patterns ^[39,43], which we will explore later.

Remark 4. (Turing-Hopf patterns and other patterns) Although it is almost impossible to theoretically prove that for system (1.2) with chemotaxis, a Turing-Hopf bifurcation of mode-$(k, 0)$ usually occurs first when the positive constant steady state is destabilized via Turing-Hopf bifurcation, some numerical simulations indicate that stable spatially inhomogeneous synchronous periodic solutions can only be generated via mode-$(k, 0)$ Turing-Hopf bifurcation.

In Theorem 11, we assert that Turing-Turing-Hopf bifurcation can occur theoretically, but it is challenging to reveal spatiotemporal patterns resulting from this type of bifurcation. In addition, although we cannot rule out the existence of Turing-Hopf-Hopf bifurcation theoretically, we never numerically found Turing-Hopf-Hopf bifurcation points after exhaustive numerical practices. But obviously, if Turing-Hopf bifurcation points such as $TH_1, \; TH_2, \cdots$ do not exist, Turing-Hopf-Hopf bifurcation must not occur.

Summarizing the results about Turing bifurcation in Section 3 and Hopf bifurcation in Section 4 in Table 2, we show the possible bifurcations for system (1.2) when the parameters are chosen in different ranges. For convenience, we use notations T-T, H-H, T-H, T-H-H and T-T-H instead of Turing-Turing, Hopf-Hopf, Turing-Hopf, Turing-Hopf-Hopf and Turing-Turing-Hopf bifurcations, respectively.

5.   Conclusions

The formation mechanisms for populations' survival patterns in a predator-taxis model with conversion delay, such as multi-stable spatial patterns and spatially staggered periodic patterns, have been discussed by studying the Turing bifurcation, Turing-Turing bifurcation, Hopf bifurcation, Hopf-Hopf bifurcation and Turing-Hopf bifurcation of system (1.2), among others.

In the absence of a delayed effect of conversion of capture behavior into predator growth, the critical condition $0 < \alpha < \alpha_{*}(d_{u}), \; 0 < d_{u} < d_{u}^{\bar{k}}$ that the positive constant steady state loses stability and system (1.2) exhibits spatial patterns is provided, which is not only sufficient, but is also necessary, and thus can be regarded as a supplement to the sufficient conditions in ^[22]. With the aid of the condition, one can predict potential spatial patterns with arbitrary wave frequencies that system (1.2) could exhibit. In addition to this, it is also theoretically expected that the positive constant steady state will be destabilized via Turing-Turing bifurcation, resulting in superposition spatial patterns or the coexistence of spatially inhomogeneous patterns with different wave frequencies; see Theorem 6 and Figures 6 and 7. Moreover, we suggest that the amplitudes of changes in the spatial distributions of the two populations are consistent with the corresponding eigenfunctions when Turing bifurcation occurs.

The above results theoretically show that with the increase of the sensitivity of the prey to predation, the spatial distributions of the two populations will gradually transition from spatially heterogeneous patterns to spatially homogeneous patterns; see Figures 3 and 4. Conversely, when the random movement speed of the prey is sufficiently low or the random movement speed of the predator is sufficiently high, the spatial distributions' homogeneity of the two populations will disappear and heterogeneity will appear; see Figures 3 and 6 or Figures 3 and 7. In other words, relatively strong predator-taxis can be regarded as accelerating the self-diffusion of the prey by complementing the self-diffusion of the prey, making two populations achieve their steady states eventually. And, a higher self-diffusion rate for the predator is equivalent to lowering the prey's ability to avoid the risk of predation, causing the spatial distributions of the two populations to be heterogeneous; however, the relatively strong predator-taxis counteracts the spatially heterogeneous distributions caused by the self-diffusion of the predator, bringing the spatial distributions of the two populations back to homogeneity. In a word, reinforcing the sensitivity of the prey to predation is beneficial for maintaining the stable survival states of the populations.

When taking a delayed effect of the conversion of capture behavior into predator growth into account, by choosing time delay $\tau$ as a parameter, we can establish the critical conditions when $\bar{E}$ is destabilized, which demonstrate that Hopf bifurcation, Hopf-Hopf bifurcation, Turing-Hopf bifurcation, etc., can also destabilize the positive constant steady state, as described by Theorems 10 and 11, leading to various spatiotemporal periodic patterns, such as stable spatially inhomogeneous staggered periodic patterns (Figures 9, 10), bistable spatiotemporal periodic patterns in which two stable periodic patterns with different spatial wave frequencies coexist (Figures 11, 12), stable spatially inhomogeneous synchronous periodic patterns (see Figure 14) and so on. We assert that the existence of stable spatially inhomogeneous staggered periodic patterns and coexistence of bistable spatiotemporal periodic patterns are due to the addition of a predator-taxis term.

These phenomena are different from the results in ^[37] for the ratio-dependent predator-prey system only with time delay, as well as the corresponding system with only predator-taxis in ^[22]. When only the delay effect of capture behavior translated into predator growth is considered, the survival of the populations is only described in a spatially homogeneous form with a single type of time-periodic oscillations. Influenced by the avoidance behavior of prey, the survival states of the populations can also be represented as time-periodic with spatially staggered oscillations. Therefore, the chemotactic behavior of the prey makes the populations' survival patterns more diverse.

Further, there are still some issues to be resolved in theory. One is how to provide a definite order for these finite Hopf bifurcation curves when the system generates stable spatially inhomogeneous periodic solutions via mode-$k$ $(0\leqslant k < K^{0})$ Hopf bifurcation, as has been shown for the Turing bifurcation curves associated with Lemma 4. Moreover, since, numerically, we only found that it is always a mode-$(k, 0)$ Turing-Hopf bifurcation that destabilizes the positive constant steady state, how to prove this finding theoretically is still a problem. Finally, it is worth mentioning that some other new spatiotemporal patterns can also arise via Turing-Turing bifurcation, Hopf-Hopf bifurcation, Turing-Hopf bifurcation and Turing-Turing-Hopf bifurcation as expected by theory, so it will be of interest to establish the existence of other coexistence patterns analytically. We leave them as possible future work.

Use of AI tools declaration

We declare that we have not used artificial intelligence tools in the creation of this article.

Acknowledgments

The research was supported by the National Natural Science Foundation of China (No. 12371165). We are grateful to the anonymous referees for their valuable comments.

Conflict of interest

We declare that there is no conflict of interest.