Global dynamics and spatiotemporal heterogeneity of a preytaxis model with prey-induced acceleration

Abstract

Conventional preytaxis systems assume that prey-tactic velocity is proportional to the prey density gradient. However, many experiments exploring the predator–prey interactions show that it is the predator’s acceleration instead of velocity that is proportional to the prey density gradient in the prey-tactic movement, which we call preytaxis with prey-induced acceleration. Mathematical models of preytaxis with prey-induced acceleration were proposed by Arditi et al. ((2001) Theor. Popul. Biol. 59(3), 207–221) and Sapoukhina et al. ((2003) Am. Nat. 162(1), 61–76) to interpret the spatial heterogeneity of predators and prey observed in experiments. This paper is devoted to exploring the qualitative behaviour of such preytaxis systems with prey-induced acceleration and establishing the global existence of classical solutions with uniform-in-time bounds in all spatial dimensions. Moreover, we prove the global stability of spatially homogeneous prey-only and coexistence steady states with decay rates under certain conditions on system parameters. For the parameters outside the stability regime, we perform linear stability analysis to find the possible patterning regimes and use numerical simulations to demonstrate that spatially inhomogeneous time-periodic patterns will typically arise from the preytaxis system with prey-induced acceleration. Noticing that conventional preytaxis systems are unable to produce spatial patterns, our results imply that the preytaxis with prey-induced acceleration is indeed more appropriate than conventional preytaxis to interpret the spatial heterogeneity resulting from predator–prey interactions.

1. Introduction

To precisely characterise the population dynamics for actively dispersing species, both random and directed movement (i.e., advection) should be considered. In the conventional framework of reaction–diffusion–advection models, the advective velocity of migrants is usually assumed to be proportional to the gradients of various biotic or abiotic stimuli, which is termed ‘taxis’, such as chemotaxis if the stimuli are chemical substances or preytaxis if the stimuli are food sources. However, there are many observations of the dependence of individual acceleration on the stimulus gradient. For instance, acceleration vectors of individuals in fish schools (cf. [36]) and in swarms of flying insects (cf. [37]) are directed towards the centroid of such dynamically stable formations, the moving flea-beetles modify their acceleration in response to food patch quality (cf. [19]), individual fish in schools adjust their variation of velocity according to the difference between ambient and preferred temperatures (cf. [11]), and the average velocity is directed by the increasing individual reproduction rate in species clustering (cf. [12]). In these observations, the directed movement of individuals is not determined by the velocity itself but by the velocity variation (i.e., acceleration) which is proportional to the gradients of stimuli. We shall call such directed movement preytaxis with prey-induced acceleration in the sequel. To understand the phenomenon of accelerated predator movement along the prey density gradient observed in Kareiva [19] and Winder et al. [51], the following reaction–diffusion–preytaxis model was proposed in [6, 40]:

(1.1) \begin{equation} \left \{\begin{array}{l@{\quad}l} u_{t}=d_{u} \Delta u-\nabla \cdot (u \mathbf{w} )+\alpha u g(u,v)-\mu u h(u), & x \in \Omega, t\gt 0, \\[5pt] v_{t}=d_{v} \Delta v+f(v)-\alpha u g(u,v), & x \in \Omega, t\gt 0, \\[5pt] \mathbf{w}_{t}=d_{w} \Delta \mathbf{w} +\gamma \nabla v, & x \in \Omega, t\gt 0, \end{array}\right. \end{equation}

where $\Omega \subset \mathbb{R}^n$ ( $n\geq 1$ ) is a bounded domain with smooth boundary, $u=u(x,t)$ and $v=v(x,t)$ represent predator and prey densities at location $x$ and at time $t$ , respectively, the vector-valued function $\mathbf{w}=(w_1,w_2,\cdots,w_n)$ denotes the velocity of predators, and $\Delta \mathbf{w}=(\Delta w_1,\Delta w_2,\cdots,\Delta w_n)$ . $d_u$ , $d_v$ , $d_w$ , $\alpha$ , $\mu$ and $\gamma$ are positive constants. $g(u,v)$ is called the trophic function (or functional response) representing the consumption of prey assuming a given number of predators, $h(u)$ is the mortality rate of predators and $f(v)$ denotes the growth function of the prey. There are many possible forms for the trophic functions $g(u,v)$ in different ecological applications, and the most commonly used include Holling type I (also called Lotka–Volterra) and Holling type II functional responses, and a summary of possible trophic functions can be found in a monograph [33]. The mortality rate function of predators is typically given by $h(u)=1+\theta u$ with $\theta \geq 0$ representing death due to intraspecific competitions. The prey growth function $f(v)$ is usually described by a logistic or bistable function. The first equation of (1.1) asserts that alongside the random diffusion the predator has an advection with the advective velocity $\mathbf{w}$ , where the temporal variation $\mathbf{w}_t$ (i.e., acceleration) is proportional to the prey density gradient (with a proportional constant $\gamma$ ) perturbed by a diffusion term $d_{w} \Delta \mathbf{w}$ (see the third equation of (1.1)) accounting for some social behaviours of species such as intraspecific competition for space or schooling effects equalising the speed and direction of neighbouring predators [11] (see more modelling details in [6, 40]). The preytaxis system (1.1) with prey-induced acceleration can generate spatiotemporal patterns (cf. [6, 40]) qualitatively consistent with the observed spatiotemporal heterogeneity in experiments in [20, 36, 37]). However, the conventional preytaxis systems featuring Lotka–Volterra interactions, where the advective velocity of the predator is directly proportional to the prey density gradient (i.e., $\mathbf{w} \sim \nabla v$ , cf. [9, 13, 20, 43] for instance), fail to achieve this explanatory power (cf. [18] or see discussions in Section 4.3.1).

There are plenty of interesting mathematical works carried out for the conventional preytaxis systems, for example, travelling wave solutions [25], pattern formation [8, 15, 26, 28, 48, 49], global solvability and stability [1, 3, 7, 30, 41, 44, 45, 53, 54], where different $g(u,v), h(u)$ and $f(v)$ may be used in different works. In contrast, the preytaxis system (1.1) with prey-induced acceleration determined by the prey density gradient was rarely studied in the literature except in a few preliminary studies as recalled below. First, by assuming that predators’ reproduction and mortality are negligible in comparison with the timescale of migration (i.e., $\alpha =\mu =0$ ), Arditi et al. [6] conduct the linear stability analysis of (1.1) with $g(u,v)=v$ (Holling type I trophic function) around the homogeneous equilibrium $(\bar{u}, 0, 0)$ with $\bar{u}= \frac{1}{|\Omega |}udx$ in a two-dimensional parallelepipedic box $\Omega$ with the zero-flux boundary condition and find that the model can produce spatial heterogeneity, in contrast to the conventional preytaxis system from which no spatial heterogeneity can arise. Later on, the work [40] performs the linear stability analysis for the model (1.1) with Holling type II trophic function $g(u,v)=\frac{v}{1+\delta v}$ , where $\delta \gt 0$ is a constant, and numerically finds the limit cycle (periodic patterns) in an interval $\Omega =[0, L]$ with the zero-flux boundary condition $u_x=v_x=w=0$ . When $h(u)=1$ , Chakraborty et al. [8] conduct the linear stability analysis for (1.1) in an interval with a variety of trophic functions $g(u,v)$ alongside numerical simulations showing the chaotic or cyclic patterns. For the Beddington–DeAngelis-type functional response, Thakur et al. [42] perform extensive simulations for $\Omega =[0, L]$ show that increasing the value of preytaxis coefficient $\gamma$ (from the bifurcation value) drives the system to exhibit chaotic behaviour, while increasing the value of random diffusion of the predator brings the system to recover from a disordered state to an ordered state. Similar observations in one dimension are observed in [39] for Holling type IV functional response (Michaelis–Menten inhibitory kinetics).

From the aforementioned results for the reaction–diffusion–advection system (1.1), we see that the existing analytical works are confined to the linear analysis and no results on the global or nonlinear dynamics seem to be available as far as we know. The goal of this paper will be to explore the global dynamics of (1.1) with prey-dependent trophic functions. Specifically, we shall consider system (1.1) in the following form:

(1.2) \begin{equation} \left \{\begin{array}{l@{\quad}l} u_{t}=d_{u} \Delta u-\nabla \cdot (u \mathbf{w} )+\alpha u F(v)-\mu u, & x \in \Omega, t\gt 0, \\[5pt] v_{t}=d_{v} \Delta v+f(v)-\alpha u F(v), & x \in \Omega, t\gt 0, \\[5pt] \mathbf{w} _{t}=d_{w} \Delta \mathbf{w} +\gamma \nabla v, & x \in \Omega, t\gt 0, \\[5pt] \nabla u \cdot \mathbf n =\nabla v \cdot \mathbf n =0,\ \mathbf{w}=\textbf{0},& x \in \partial \Omega, t\gt 0,\\[5pt] (u,v,\mathbf{w})(x,0)=(u_0,v_0,\mathbf{w}_0)(x),& x \in \Omega, \end{array}\right. \end{equation}

where $\Omega \subset \mathbb{R}^n$ ( $n\geq 1$ ) is a bounded domain with smooth boundary and $\mathbf n$ is the unit outward normal vector of $\partial \Omega$ . Writing $d_{u} \Delta u-\nabla \cdot (u \mathbf{w})=\nabla \cdot (d_u \nabla u-u \mathbf{w})$ , we find that the above boundary conditions give rise to the zero-flux boundary conditions meaning that both predator and prey live in a closed habitat and cannot cross the boundary.

In the present paper, the initial data $(u_0, v_0, \mathbf{w}_0)$ are supposed to satisfy

(1.3) \begin{eqnarray} u_{0}(x), v_{0}(x) \gvertneqq 0,\ u_{0}, v_{0} \in W^{1, \infty }(\Omega ),\ \mathbf{w}_0(x)\in [W^{1, \infty }(\Omega )]^n. \end{eqnarray}

Moreover, we suppose that the trophic function $F(v)$ and the growth function $f(v)$ satisfy the following hypotheses:

1. (H1) $F(v)\in C^{2}([0, \infty ))$ , $F(0)=0$ and $F(v)\gt 0$ in $(0, \infty ).$

2. (H2) $f \in C^{2}([0, \infty ))$ , $f(0)=0$ , and there exist two positive constants $\eta, K$ such that $f(v) \leq \eta v$ for $v \geq 0$ , $f(K)=0$ and $f(v)\lt 0$ for $v\gt K$ .

We remark that the above assumptions for $F(v)$ and $f(v)$ have covered a large class of typical examples such as:

\begin{equation*} \left \{ \begin {array}{l} F(v)=v \text { (Lotka-Volterra type or Holling type I)}, \\[10pt] F(v)=\dfrac {v}{\lambda +v}\ (\text {Holling type II}), \\[10pt] F(v)=\dfrac {v^{\kappa }}{\lambda ^{\kappa }+v^{\kappa }}\ (\text {Holling type III}), \end {array}\right. \end{equation*}

and

\begin{equation*} \left \{ \begin {array}{l} f(v)=\eta v\!\left (1-\dfrac {v}{K}\right ) \text { (Logistic type)}, \\[10pt] f(v)=\eta v\!\left (1-\dfrac {v}{K}\right )\left (\dfrac {v}{L}-1\right )\ (\text {Allee effect type)},\\ \end {array}\right. \end{equation*}

where $\lambda$ , $\kappa$ and $L$ are positive constants with $\kappa \gt 1$ , $0\lt L\lt K$ .

In this paper, we shall establish the global existence and stability of classical solutions of (1.2) with a generic prey-dependent trophic function, refine the linear analysis to identify the parameter regimes for pattern formations and use numerical simulations to demonstrate the spatiotemporal patterns generated by (1.2) implying that the preytaxis with prey-induced acceleration is more appropriate than the conventional one to interpret the field observation of spatially heterogeneous coexistence in predator–prey systems. The first main theorem stated below asserts that the problem (1.2) has a global-in-time classical solution which is uniformly bounded with respect to time in any spatial dimensions.

Theorem 1.1. Let $n\geq 1$ and the hypotheses (H1)–(H2) and (1.3) hold. Then the problem (1.2) has a unique global classical solution $(u, v, \mathbf{w})$ satisfying

\begin{equation*} {u(x,t)\geq 0\text { and }0\leq v(x,t)\leq m\quad \text {for all }(x,t)\in \bar \Omega \times (0,+\infty ),} \end{equation*}

and

\begin{align*}{\limsup \limits _{t\rightarrow +\infty }v(x,t)\leq K\quad \textit{for all }x\in \bar \Omega,} \end{align*}

where the positive constant $m$ is defined by:

(1.4) \begin{equation} m=\max \!\left \{\|v_{0}\|_{L^\infty (\Omega )}, K\right \}. \end{equation}

Moreover, we have

\begin{equation*} \|u\|_{L^{\infty }(\Omega )}+\|v(\cdot,t)\|_{W^{1, \infty }(\Omega )}+\|{\mathbf{w} }(\cdot,t)\|_{W^{1, \infty }(\Omega )}\leq C\quad \textit{for all }t\gt 0. \end{equation*}

where $C\gt 0$ is a constant independent of time $t$ .

In the following, we shall present the global stability of non-negative constant prey-only and coexistence steady states. The system (1.2) has three possible constant steady states $\left (u_s, v_s, \mathbf{w}_s\right )$ :

(1.5) \begin{equation} \left (u_s, v_s, \mathbf{w}_s\right )=\left \{\begin{array}{l@{\quad}l} (0,0,\textbf{0}),\left (0, 1, \textbf{0}\right ), & \textrm{ if } \alpha F(K)\leq \mu,\\[5pt] (0,0,\textbf{0}),\left (0, 1, \textbf{0}\right ),\left (u_*,v_*,\textbf{0}\right ), & \textrm{ if } \alpha F(K)\gt \mu, \end{array}\right. \end{equation}

where $(0,0,\textbf{0})$ is the trivial steady state (i.e., extinction steady state), $\left (0, 1, \textbf{0}\right )$ is the semi-trivial steady state (called the prey-only steady state) and $(u_*,v_*,\mathbf{0})$ is the coexistence steady state with $u_*$ , $v_*\gt 0$ determined by the following algebraic equations:

(1.6) \begin{equation} \alpha F(v_*)=\mu \quad \textrm{and}\quad \mu u_*=f(v_*). \end{equation}

Note that the non-negative constant steady states of $\mathbf{w}$ is $\textbf{0}$ due to $\mathbf{w}=\textbf{0}$ on $\partial \Omega$ .

For the global stability of steady states given in (1.5), except the hypotheses (H1) and (H2), we need additional assumptions for $F(v)$ and the compound function:

\begin{equation*}\Phi (v)=\frac {f(v)}{F(v)}\end{equation*}

as follows:

1. (H3) $F^{\prime}(v)\gt 0$ in $[0, \infty )$ .

2. (H4) $\Phi (v)\in C^1((0,\infty )),\ \Phi (0)=\lim _{v\rightarrow 0^+}\Phi (v)\gt 0$ and $\Phi ^{\prime}(v)\lt 0$ in $(0, \infty )$ .

It follows from Theorem 1.1 that there exists some $T_0\gt 0$ such that

(1.7) \begin{align} 0\lt v(x,t)\lt K+1\quad \textrm{for all }(x,t)\in \bar \Omega \times (T_0,+\infty ). \end{align}

Then, $F^{\prime}(v)$ reaches its positive minimum in the interval $[0,{K+1}]$ and $F^2(v)$ has an upper bound in the interval $[0,{K+1}]$ for $t\gt T_0$ . Let

(1.8) \begin{eqnarray} c_0=\inf _{v\in [0,{K+1}]}\Psi (v), \quad \textrm{where }\Psi (v)=\frac{F^{\prime}(v)}{F^2(v)}, \end{eqnarray}

then $c_0\gt 0$ is a constant independent of the initial data. We shall also use the following Poincaré constant for the domain $\Omega$ :

(1.9) \begin{equation} C_P(\Omega )=\inf \!\left \{C\gt 0\!\left |\|\boldsymbol{\varphi }\|_{L^{2}(\Omega )}\leq C \|\nabla \boldsymbol{\varphi }\|_{L^{2}(\Omega )}\ \textrm{for all }\boldsymbol{\varphi }\in \left(W^{1,2}_0(\Omega )\right)^n\right .\right \}. \end{equation}

In the case of $\alpha F(K)\gt \mu$ , we have global stability of the coexistence steady state.

Theorem 1.2. Assume that $\alpha F(K)\gt \mu$ . Let $n\geq 1$ and the hypotheses (H1)–(H4) and (1.3) hold, and let $(u, v, \mathbf{w} )$ be the solution of (1.2) obtained in Theorem 1.1. If

(1.10) \begin{equation} d_w\gt \left (\frac{u_*}{4d_u}+\frac{\alpha \gamma ^2}{d_v c_0\mu }\right )\frac{C_P^2(\Omega )}2, \end{equation}

then the coexistence steady state $(u_*,v_*,\textbf{0})$ is globally asymptotically stable.

Remark 1.1. The positive number $\left (\frac{u_*}{4d_u}+\frac{\alpha \gamma ^2}{d_v c_0\mu }\right )\frac{C_P^2(\Omega )}2$ on the right side of (1.10) is determined and can be estimated once the domain $\Omega$ , parameter values and the tropic function $F(v)$ are specified. For instance, if $F(v)=v$ (Holling type I) or $F(v)=\frac{v}{1+v}$ (Holling type II), then the constant $c_0\gt 0$ is determined by:

\begin{equation*} c_0=\inf _{v\in [0,{K+1}]}\Psi (v)=\inf _{v\in [0,{K+1}]}\frac {F^{\prime}(v)}{F^2(v)}=\inf _{v\in [0,{K+1}]}\frac 1{v^2}=\frac 1{{(K+1)}^2}. \end{equation*}

For an open bounded set $\Omega \subset \mathbb{R}^n$ , $C_P(\Omega )\leq \frac{2}{n} \Big(\inf \limits _{x \in \mathbb{R}^n} \sup \limits _{y \in \Omega }|x-y|\Big)^{\frac 12}$ (cf. [27, Exercise 13.22]). Moreover, $C_P$ can be specified in some special cases, such as $C_P=\frac 1{\pi }$ if $\Omega$ is a unit isosceles right triangle (cf. [22]).

Remark 1.2. Consider a special case: $\mathbf{w}$ is a conservative vector field with $\mathbf{w}=\nabla \phi$ for some scalar potential function $\phi \in C^1(\overline{\Omega })$ . Then the system (1.2) can be written as an indirect preytaxis system as follows:

(1.11) \begin{align} \left \{ \begin{array}{l@{\quad}l@{\quad}l} u_{t}=d_{u} \Delta u-\nabla \cdot (u \nabla \phi )+\alpha u F(v)-\mu u, & x \in \Omega, t\gt 0, \\[5pt] v_{t}=d_{v} \Delta v+f(v)-\alpha u F(v), & x \in \Omega, t\gt 0, \\[5pt] \phi _t=d_w\Delta \phi +\gamma v+\lambda _0, & x \in \Omega, t\gt 0, \\[5pt] \nabla u \cdot \mathbf n =\nabla v \cdot \mathbf n =\nabla \phi \cdot \mathbf n =0,& x \in \partial \Omega, t\gt 0,\\[5pt] (u,v,\phi )(x,0)=(u_0,v_0,\phi _0)(x),& x \in \Omega \end{array} \right. \end{align}

for some constant $\lambda _0 \in \mathbb{R}$ . Recently, there are some works for the above indirect preytaxis model with $\lambda _0=0$ for different functional response functions (cf. [1, 2, 31, 46, 55, 56] and references therein). In this sense, the indirect preytaxis system can be regarded as a special case of (1.2) when $\mathbf{w}$ is a gradient field.

In the case of $\alpha F(K)\leq \mu$ , we have the global stability of the prey-only steady state.

Theorem 1.3. Assume that $\alpha F(K)\leq \mu$ . Let $n\geq 1$ and the hypotheses (H1)–(H4) and (1.3) hold, and let $(u, v, \mathbf{w} )$ be the solution of (1.2) obtained in Theorem 1.1. Then for any positive parameters $d_u$ , $d_v$ , $d_w$ and $\gamma$ , the prey-only steady state $\left (0,K,\textbf{0}\right )$ is globally asymptotically stable. Moreover, if $\alpha F(K)\lt \mu$ , then $(0,K,\textbf{0})$ is exponentially stable, that is, there exist positive constants $C$ , $\lambda$ and $T_1$ such that

\begin{equation*} \|u\|_{L^\infty (\Omega )}+\|v-K\|_{L^\infty (\Omega )}+\|\mathbf{w} \|_{L^\infty (\Omega )} \leq C e^{-\lambda t}\quad \textit{for all }t\gt T_1. \end{equation*}

The rest of this paper is organised as follows. In Section 2, we establish the existence of globally bounded classical solutions of (1.2) by extending local solutions with the a priori estimates of solutions derived. In Section 3, we show the global stability of coexistence and prey-only steady states by constructing Lyapunov functionals along with some compactness arguments. In Section 4, we conduct linear stability analysis to identify the parameter regime for the pattern formation and perform numerical simulations to show that the preytaxis system (1.2) with prey-induced acceleration will typically generate spatially inhomogeneous time-periodic patterns which are well consistent with the experimental observations.

2. Global existence and uniform boundedness

In this section, we shall establish the global existence and boundedness of solutions to (1.2), which consists of local existence and some a priori estimates of solutions. Before proceeding, we introduce some notations frequently used in the paper.

Notations.

• • Without confusion, we shall abbreviate $\int _{\Omega } f d x$ as $\int _{\Omega } f$ for simplicity.

• • We denote by $C$ , $C_i$ ( $i=1,2,3,\cdots$ ) and $C_{\Omega }$ generic positive constants that may vary in the context, where $C$ and $C_i$ are independent of $\Omega$ , and $C_{\Omega }$ depends only on $\Omega$ .

2.1. Preliminary results

We first use Amann’s theorem of parabolic systems in [4, 5] (cf. also [50, Lemma 2.6]) to establish the existence of local-in-time classical solutions of the system (1.2).

Lemma 2.1. Let $n\geq 1$ and the hypotheses (H1), (H2) and (1.3) hold. Then there exists ${T_{\max}} \in (0, \infty ]$ such that (1.2) admits a unique classical solution $(u, v, \mathbf{w})$ on $[0,{T_{\max}})$ satisfying

\begin{equation*} \left \{\begin {array}{l@{\quad}c@{\quad}l} u, v \in C\!\left(\bar {\Omega } \times [0, {T_{max}})\right) \cap C^{2,1}\!\left(\bar {\Omega } \times \left(0, {T_{max}}\right)\right),\\[5pt] \mathbf {w} \in \left [C(\bar {\Omega } \times [0, {T_{max}})) \cap C^{2,1}\!\left(\bar {\Omega } \times \left(0, {T_{max}}\right)\right)\right ]^n \end {array}\right. \end{equation*}

and

(2.1) \begin{equation} u(x,t)\gt 0, \ 0\lt v(x,t) \leq m\textit{ in } \Omega \times (0,{T_{max}}), \end{equation}

where $m$ is given by (1.4). Moreover, there is a dichotomous criterion:

(2.2) \begin{equation} \textit{either ${T_{max}}=\infty,$}\ \textit{or}\ \lim _{t \rightarrow{T_{max}}} \sup \!\left (\|u(\cdot, t)\|_{L^\infty (\Omega )}+\|v(\cdot, t)\|_{L^\infty (\Omega )}+\|\mathbf{w} (\cdot, t)\|_{L^\infty (\Omega )}\right )=\infty. \end{equation}

Proof. Let $\boldsymbol \psi =\left(\psi _1,\psi _2,\cdots,\psi _{n+2}\right)^{T}=(u, v, \mathbf{w})^{T}=(u, v, w_1,w_2,\cdots,w_n)^{T}$ be a $(n+2)$ -dimensional vector-valued function, where $\boldsymbol{K}^{T}$ denotes the transpose of a matrix $\boldsymbol{K}$ . Denote $\mathbf{0}_{p\times q}$ by a $p$ -by- $q$ zero matrix with two positive integers $p$ and $q$ . Let

\begin{equation*}\mathbf {\xi }_i=\left(-\psi _{i+2},\textbf {0}_{1\times i},-\psi _1,\textbf {0}_{1\times (n-i)}\right),\quad i=1,2,\cdots,n,\end{equation*}

be a $(n+2)$ -dimensional vector-valued function, and

\begin{equation*} \mathbf {D}_i=\left ( \begin {array}{c} \mathbf {\xi }_i\\[5pt] \mathbf P_i \end {array} \right ),\quad i=1,2,\cdots,n,\end{equation*}

be a square matrix of order $(n+2)$ , where all elements of the $(n+1)$ -by- $(n+2)$ matrix $\mathbf P_i$ are 0 except the $(i+1)$ -by- $2$ element is $\gamma$ . Then the system (1.2) can be rewritten as:

(2.3) \begin{equation} \left \{\begin{array}{l@{\quad}l} \boldsymbol \psi _{t}=\boldsymbol{A}\cdot \Delta \boldsymbol \psi +\sum _{i=1}^n\boldsymbol{D}_i\cdot \partial _i\boldsymbol \psi +\boldsymbol{F}, & x \in \Omega, t\gt 0, \\[5pt] \mathcal{B}\boldsymbol \psi =0, & x \in \partial \Omega, t\gt 0, \\[5pt] \boldsymbol \psi (\cdot, 0)=\left (u_{0}, v_{0}, \mathbf{w}_0\right ), & x \in \Omega, \end{array}\right. \end{equation}

where

\begin{equation*} \boldsymbol{A}= \left ( \begin {array}{c@{\quad}c@{\quad}c} d_u & 0 & 0\\[5pt] 0 & d_v & 0\\[5pt] 0 & 0 & d_w\boldsymbol{E}_n\\[5pt] \end {array} \right ) \end{equation*}

is a constant square matrix of order $(n+2)$ with $\boldsymbol{E}_n$ being the identity matrix of order $n$ , and $\boldsymbol{F}$ is a $(n+2)$ -dimensional vector-valued function given by:

\begin{equation*}\boldsymbol{F}=\left ( \begin {array}{c} \alpha \psi _1 F(\psi _2)\boldsymbol -\mu \psi _1\\[5pt] f(\psi _2)-\alpha \psi _1 F(\psi _2)\\[5pt] \textbf {0}_{n\times 1} \end {array} \right ). \end{equation*}

Moreover, the boundary operator $\mathcal{B}$ is given by:

\begin{equation*}\mathcal {B}=\left ( \begin {array}{c@{\quad}c@{\quad}c} \partial _{\mathbf n} & & \\[5pt] & \partial _{\mathbf n} & \\[5pt] & & \boldsymbol{E}_n\\[5pt] \end {array} \right ),\end{equation*}

where $\partial _{\mathbf n}$ is the partial derivative with respect to $\mathbf n$ . Obviously, all eigenvalues of $\boldsymbol{A}$ are positive, and hence system (2.3) is uniformly parabolic. The local existence and uniqueness of classical solutions follow from Amann’s theorem [4, Theorem 7.3 and Corollary 9.3] and the blow-up criteria (2.2) follows from [5, Theorem 15.5].

The positivity of $u$ and $v$ for $t\in (0,T_{\max})$ follows from the strong maximum principle. To be precise, we rewrite the first equation of the system (1.2) as follows:

\begin{equation*} \left \{\begin{array}{l@{\quad}l} u_{t}-d_{u} \Delta u+\mathbf{w}\cdot \nabla u+q(x, t) u=0, & x \in \Omega, t \in \left(0, T_{\max}\right), \\[5pt] \nabla u \cdot{\mathbf n}=0, & x \in \partial \Omega, t \in \left(0, T_{\max}\right), \\[5pt] u(x, 0)=u_{0} \geq 0(\!\not \equiv 0), & x \in \Omega, \end{array}\right. \end{equation*}

where $q(x,t)=\nabla \cdot \mathbf{w}-\alpha F(v)+\mu$ . By the maximum principle, we know that $u\geq 0$ in $\Omega \times \left(0, T_{\max}\right)$ , we shall show that actually $u\gt 0$ in $\Omega \times \left(0, T_{\max}\right)$ . For any $(x_*,t_*)\in \Omega \times \left(0, T_{\max}\right)$ , we can find an open subset $\Omega _*\subset \Omega$ and $T_*\in \left(0, T_{\max}\right)$ such that

(2.4) \begin{equation} \left \{ \begin{array}{l@{\quad}c@{\quad}l} (x_*,t_*)\in \Omega _*\times (0,T_*)\,=\!:\,Q_{(x_*,t_*)}\subset \Omega \times \left(0, T_{\max}\right),\\[5pt] \exists \ x_0\in \Omega _*\textrm{ such that }u(x_0,0)\gt 0, \end{array}\right. \end{equation}

where the second condition can be satisfied according to (1.3). By the regularity of $q(x,t)$ , we can find some constant $R$ such that $R=\inf _{{Q_{(x_*,t_*)}}}q(x,t)$ , and hence $U(x, t)\,:\!=\,\textrm{e}^{R t} u(x, t)\geq 0$ for $(x,t)\in \Omega \times [0, T_{\max})$ satisfies

\begin{equation*} U_{t}-d_{u} \Delta U+\mathbf{w} \cdot \nabla U+(q(x, t)-R) U=0,\quad (x, t) \in Q_{(x_*,t_*)}. \end{equation*}

If $u(x_*,t_*)=0$ , then by $q(x, t)-R\geq 0$ , $U(x_*,t_*)=0$ and $U(x,t)\geq 0$ for $(x,t)\in Q_{(x_*,t_*)}$ , one can apply the strong maximum principle [29, Lemma 2.7] to obtain $U(x,t)\equiv 0$ for $(x,t)\in \Omega _*\times (0,t_*)$ . This together with the continuity of $U$ yields $U(x,0)\equiv 0$ for $x\in \Omega _*$ , which contradicts the second condition of (2.4). Hence, we have $u(x_*,t_*)\neq 0$ , that is, $u(x_*,t_*)\gt 0$ due to $u(x,t)\geq 0$ for $(x,t)\in \Omega \times (0,T_{\max})$ . Since $(x_*,t_*)\in \Omega \times \left(0, T_{\max}\right)$ is arbitrary, we have $u\gt 0$ in $\Omega \times \left(0, T_{\max}\right)$ . Similarly, using the strong maximum principle one can show that $0\lt v\leq m$ in $\Omega \times \left(0, T_{\max}\right)$ . Therefore, (2.1) is proved, and the proof is completed.

Lemma 2.2. For all $t\in (0,{T_{\max}})$ , there exists a constant $C\gt 0$ such that

\begin{eqnarray*} \|u (\cdot,t)\|_{L^1(\Omega )}\leq C. \end{eqnarray*}

Proof. It follows from the first two equations in (1.2) that

\begin{equation*} \frac d{dt}\int _\Omega (u+v)+\mu \int _\Omega u\leq \int _\Omega f(v)\quad \textrm{for all }t\in (0,{T_{\max}}). \end{equation*}

By the assumption (H2) and (2.1), we have

\begin{equation*} \frac d{dt}\int _\Omega (u+v)+\mu \int _\Omega (u+v)\leq \left(\max \limits _{v\in [0, m]}f(v)+\mu m\right)|\Omega |\leq C\quad \textrm{for all }t\in (0,{T_{\max}}), \end{equation*}

where $m$ is given by (1.4). Therefore, an application of Grönwall’s inequality completes the proof.

Since the boundary condition of $\mathbf{w}$ in (1.2) is Dirichlet type, we need the following $L^{p}$ - $L^{q}$ -estimates for the Dirichlet heat semigroup established in [38, Proposition 48.4^*, 48.5 and 48.7^*].

Lemma 2.3. ([38, Proposition 48.4^*, 48.5 and 48.7^*]) Let $n\geq 1$ and $\left (e^{t\Delta }\right )_{t \geq 0}$ be the Dirichlet heat semigroup in $\Omega$ and let $\lambda _1 \gt 0$ denote the first non-zero eigenvalue of $-\Delta$ in $\Omega$ under Dirichlet boundary conditions. We have the following properties.

1. (i) If $1 \leq q\lt p \leq \infty$ , then

   (2.5) \begin{equation} \left\|e^{t\Delta } z \right\|_{L^p(\Omega )} \leq (4 \pi t)^{-\frac n2\left(\frac 1q-\frac 1p\right)}\|z \|_{L^q(\Omega )}\quad \textit{for all }t\gt 0, \end{equation}
   holds for all $z \in L^{q}(\Omega )$ .

2. (ii) For all $1 \leq p \leq \infty$ and all $z \in L^{p}(\Omega ),$ it holds that

   (2.6) \begin{equation} \|e^{t\Delta } z\|_{L^p(\Omega )} \leq C_{\Omega } e^{-\lambda _1 t}\|z\|_{L^p(\Omega )}\quad \textit{for all }t\gt 0. \end{equation}

3. (iii) For all $z \in L^{\infty }(\Omega )$ , it holds that

   (2.7) \begin{equation} \|\nabla e^{t\Delta } z \|_{L^\infty (\Omega )} \leq C_{\Omega }\big(1+t^{-\frac 12}\big)\|z \|_{L^\infty (\Omega )}\quad \textit{for all }t\gt 0. \end{equation}

Next we give some further $L^{p}$ - $L^{q}$ -estimates for the Dirichlet heat semigroup, which can be deduced based on Lemma 2.3 and a similar argument as the proof of [52, Lemma 1.3]. Although some of the following results seem to be well known, we cannot find precise references in the literature that provide all estimates that we need for our purpose, and therefore we provide some proof as a complement.

Lemma 2.4. Let $e^{t{\Delta }}$ be the Dirichlet heat semigroup in $\Omega \subset \mathbb{R}^n$ $(n\geq 1)$ , $\lambda _1 \gt 0$ denote the first non-zero eigenvalue of $-\Delta$ in $\Omega$ under the Dirichlet boundary condition. Then the following properties hold.

1. (i) If $1 \leq q\leq p \leq \infty$ , then for any $z \in L^{q}(\Omega )$ , it holds that

   (2.8) \begin{eqnarray} \|e^{{t{\Delta }}} z \|_{L^p(\Omega )} &\leq & C_{\Omega }\!\left(1+t^{-\frac n2\left(\frac 1q-\frac 1p\right)}\right)e^{-\lambda _1 t}\|z \|_{L^q(\Omega )}\quad \textit{for all }t\gt 0, \end{eqnarray}
   and
   (2.9) \begin{eqnarray} \|\nabla e^{t\Delta }z\|_{L^p(\Omega )}&\leq & C_{\Omega }\!\left(1+t^{-\frac 12-\frac n2\left(\frac 1q-\frac 1p\right)}\right )e^{-\lambda _1 t}\|z\|_{L^q(\Omega )}\quad \textit{for all }t\gt 0. \end{eqnarray}

2. (ii) If $2\leq p\lt \infty$ , then for any $z\in W^{1,p}(\Omega )$ , it holds

   (2.10) \begin{eqnarray} \|\nabla e^{t\Delta }z\|_{L^p(\Omega )}\leq C_{\Omega } e^{-\lambda _1 t}\|\nabla z\|_{L^p(\Omega )}\quad \textit{for all }t\gt 0. \end{eqnarray}

3. (iii) If $1\lt q\leq p\leq \infty$ , then for $\textbf{z}\in \left (L^{q}(\Omega )\right )^{n}$ , one has

   (2.11) \begin{equation} \|e^{t\Delta }\nabla \cdot \textbf{z}\|_{L^p(\Omega )}\leq C_{\Omega }\!\left(1+t^{-\frac 12-\frac n2\left(\frac 1q-\frac 1p\right)}\right)e^{-\lambda _1 t}\|\textbf{z}\|_{L^q(\Omega )}\quad \textit{for all }t\gt 0. \end{equation}

Proof. (i) We first prove (2.8). For $1\leq q=p\leq \infty$ , (2.8) is a direct consequence of (2.6). For $1\leq q\lt p\leq \infty$ and $t\lt 2$ , (2.8) is a direct consequence of (2.5) since $e^{-\lambda _1 t}\gt e^{-2\lambda _1}$ for $t\lt 2$ . For $1 \leq q\lt p \leq \infty$ and $t\geq 2$ , using (2.5) and (2.6), we have

\begin{eqnarray*} \|e^{t \Delta } z\|_{L^p(\Omega )} &=& \|e^{(t-1) \Delta }e^{\Delta } z\|_{L^p(\Omega )}\\[5pt] &\leq & C_{\Omega } e^{-\lambda _1 (t-1)}\|e^{\Delta } z\|_{L^p(\Omega )} \\[5pt] &\leq & C_{\Omega } e^{-\lambda _1 t}\|z\|_{L^q(\Omega )}\quad \textrm{for all }t\geq 2. \end{eqnarray*}

This completes the proof of (2.8).

It remains to prove (2.9). Using the pointwise estimates for the spatial gradient of Green’s function of $e^{t\Delta }$ in [23] (see also [32, Theorem 2.2]), we can find constants $C_1,C_2,C_3\gt 0$ such that

(2.12) \begin{equation} \|\nabla e^{t\Delta } z\|_{L^{1}(\Omega )}\leq C_1 t^{-\frac{n+1}2}\int _\Omega e^{-C_2\frac{|x|^2}t}dx\|z\|_{L^1(\Omega )}\leq C_3t^{-\frac 12}\|z\|_{L^1 (\Omega )}\quad \textrm{for all }t\gt 0, \end{equation}

holds for all $z\in L^{1}(\Omega )$ . For any $t\in (0,+\infty )$ , define the map $\mathcal T_t$ on $L^{p}(\Omega )$ for $p\in [1,\infty ]$ by: $\mathcal T_t(z)=\nabla e^{t\Delta } z$ for $z\in L^{p}(\Omega )$ . It follows from (2.7) and (2.12) that for all $t\gt 0$ , $\mathcal T_t$ is of weak type $(1,1)$ and weak type $(\infty,\infty )$ (see [10, Section 9 of Chapter 9] for the definitions of weak type and $\|\cdot \|_{p,w}$ for $p\in [1,\infty ]$ ) since

\begin{equation*} \left \{ \begin{array}{l@{\quad}c@{\quad}l} \|\mathcal T_t z\|_{1,w}\leq \|\nabla e^{t\Delta } z\|_{L^{1}(\Omega )}\leq C_3t^{-\frac 12}\|z\|_{L^{1}(\Omega )}\leq C_3\big(1+t^{-\frac 12}\big)\|z\|_{L^{1}(\Omega )},\\[5pt] \|\mathcal T_t z\|_{\infty,w}=\|\nabla e^{t\Delta } z \|_{L^\infty (\Omega )} \leq C_{\Omega }\big(1+t^{-\frac 12}\big)\|z \|_{L^\infty (\Omega )}. \end{array}\right. \end{equation*}

Using the Marcinkiewicz interpolation theorem (cf. [10, Theorem 9.1]), we find that $\mathcal T_t$ is of strong type $(r,r)$ for $1\lt r\lt \infty$ : $\|\mathcal T_t z\|_r\leq C(r,t)\|z\|_{L^{r}(\Omega )}$ , where

\begin{equation*} C(r,t)=C\!\left (\frac r{r-1}\right )^{\frac 1r}\left (C_3\big(1+t^{-\frac 12}\big)\right )^\frac 1r\left (C_\Omega \big(1+t^{-\frac 12}\big)\right )^{1-\frac 1r}\,=\!:\, C(\Omega,r)\big(1+t^{-\frac 12}\big). \end{equation*}

That is, $\|\nabla e^{t\Delta } z \|_{L^r(\Omega )} \leq C(\Omega,r)\big(1+t^{-\frac 12}\big)\|z \|_{L^r(\Omega )}$ for $1\lt r\lt \infty$ , where the constant $C(\Omega,r)$ depends only on $\Omega$ and $r$ . Therefore, this along with (2.7) and (2.12) yields for $1\leq p\leq \infty$ that

(2.13) \begin{equation} \|\nabla e^{t\Delta } z \|_{L^p(\Omega )} \leq C(\Omega,p)\big(1+t^{-\frac 12}\big)\|z \|_{L^p(\Omega )}\quad \textrm{ for all } t\gt 0, \end{equation}

where the constant $C(\Omega,p)$ depends only on $\Omega$ and $p$ .

For $0\lt t\lt 2$ , we can use (2.8) and (2.13) to obtain

(2.14) \begin{eqnarray} \|\nabla e^{t \Delta } z \|_{L^p(\Omega )} &=&\|\nabla e^{\frac{t}{2} \Delta } e^{\frac{t}{2} \Delta }z \|_{L^p(\Omega )} \nonumber \\[5pt] &\leq & C\!\left (1+(t/2)^{-\frac{1}{2}}\right )\|e^{\frac{t}{2} \Delta }z \|_{L^p(\Omega )} \nonumber \\[5pt] &\leq & C_{\Omega }\!\left (1+(t/2)^{-\frac{1}{2}}\right )\left (1+\left ((t/2)\right )^{-\frac n2\left(\frac 1q-\frac 1p\right)}\right )e^{-\lambda _1 \frac t2}\|z \|_{L^q(\Omega )}\nonumber \\[5pt] &\leq & C_{\Omega }(t/2)^{-\frac{1}{2}-\frac n2\left(\frac 1q-\frac 1p\right)}e^{-\lambda _1 \frac t2}\|z \|_{L^q(\Omega )}\nonumber \\[5pt] &\leq & C_{\Omega } t^{-\frac 12-\frac n2\left(\frac 1q-\frac 1p\right)} e^{-\lambda _1 t}\|z \|_{L^q(\Omega )}, \end{eqnarray}

where we have used $1\leq \min \!\left \{\left (\frac{t}{2}\right )^{-\frac{1}{2}},\left (\frac t2\right )^{-\frac n2\left(\frac 1q-\frac 1p\right)}\right \}$ . This implies (2.9) for $t\lt 2$ . For $t\geq 2$ , we can use (2.8) and (2.13) to obtain

\begin{eqnarray*} \|\nabla e^{t \Delta } z \|_{L^p(\Omega )} &=&\|\nabla e^{\Delta } e^{(t-1) \Delta }z \|_{L^p(\Omega )} \leq C_{\Omega }\|e^{(t-1) \Delta }z \|_{L^p(\Omega )} \\[5pt] &\leq & C_{\Omega }\!\left (1+\left (t-1\right )^{-\frac n2\left(\frac 1q-\frac 1p\right)}\right )e^{-\lambda _1 (t-1)} \|z \|_{L^q(\Omega )} \\[5pt] &\leq & C_{\Omega }e^{-\lambda _1 (t-1)}\|z \|_{L^q(\Omega )}. \end{eqnarray*}

This together with (2.14) proves (2.9).

(ii) Recall that the following two inequalities hold also for the Dirichlet heat semigroup (cf. [7, formula (1.12)] and [32, formula (2.39)]),

\begin{equation*} \|\nabla e^{t \Delta } z\|_{L^2(\Omega )} \leq \|\nabla z\|_{L^2(\Omega )}\quad \text { for all } t \geq 0, \end{equation*}

and

\begin{equation*} \|\nabla e^{t \Delta } z\|_{L^\infty (\Omega )}\leq C\|\nabla z\|_{L^\infty (\Omega )}\quad \text { for all } t \in (0,1). \end{equation*}

Then (2.10) can be readily derived by a process similar to the proof of [52, Lemma 1.3 (iii)].

(iii) With (2.9) at hand, (2.11) can be proved by an argument similar to the proof of [52, Lemma 1.3 (iv)]. Although [52, Lemma 1.3 (iv)] is stated only for the case of $1\lt q\leq p\lt \infty$ , but the proof actually already covers the case $p=q=\infty$ (cf. [24, Lemma 3.1]).

2.2. A priori estimates (Proof of Theorem 1.1)

We first derive the $L^\infty$ -estimate for $\mathbf{w}$ , which is a direct consequence of the above $L^{p}$ - $L^{q}$ -estimates for the Dirichlet heat semigroup.

Lemma 2.5. For all $t\in (0,{T_{\max}})$ , there exists a constant $C\gt 0$ such that

(2.15) \begin{equation} \|\mathbf{w} (\cdot,t)\|_{L^\infty (\Omega )}\leq C. \end{equation}

Proof. By Duhamel’s principle, one has

\begin{eqnarray*} \mathbf{w}(t)&=& e^{d_w t\Delta }\mathbf{w}_0+\gamma \int ^t_0e^{d_w (t-s)\Delta }\nabla v(\cdot,s)ds,\\[5pt] &=& e^{d_w t\Delta }\mathbf{w}_0+\gamma \int ^t_0\nabla e^{d_w (t-s)\Delta } v(\cdot,s)ds\quad \textrm{for all }t\in (0,{T_{\max}}). \end{eqnarray*}

By Lemma 2.4 and (2.1), we have

\begin{eqnarray*} \|\mathbf{w}(\cdot,t)\|_{L^{\infty }(\Omega )}&\leq & C\|\mathbf{w}_0\|_{L^{\infty }(\Omega )}+C\int ^t_0\left(1+(t-s)^{-\frac 12}\right)e^{-d_w\lambda _1 (t-s)}\|v(\cdot,s)\|_{L^{\infty }(\Omega )} ds\\[5pt] &\leq & C\|\mathbf{w}_0\|_{L^{\infty }(\Omega )}+C\int ^t_0\left(1+(t-s)^{-\frac 12}\right)e^{-d_w\lambda _1 (t-s)}ds\\[5pt] &\leq & C\quad \textrm{for all }t\in (0,{T_{\max}}). \end{eqnarray*}

This completes the proof.

With (2.1) and (2.15) at hand, we can proceed to derive a priori $L^\infty$ -estimate of $u$ .

Lemma 2.6. For all $t\in (0,{T_{\max}})$ , there exists a constant $M\gt 0$ such that

(2.16) \begin{equation} \|u(\cdot,t)\|_{L^{\infty }(\Omega )}\leq M. \end{equation}

Proof. Multiplying the first equation of (1.2) by $pu^{p-1}$ ( $p\gt 1$ ), integrating the result with respect to $x$ over $\Omega$ and using (2.1), one has

(2.17) \begin{eqnarray} &&\frac d{dt}\int _\Omega u^p+p\mu \int _\Omega u^p+d_u (p-1)\frac 4{p}\int _\Omega \big|\nabla u^\frac p2\big|^2\nonumber \\[5pt] &=& p(p-1)\int _\Omega u^{p-1}\mathbf{w}\cdot \nabla u+\alpha p\int _\Omega F(v)u^p\nonumber \\[5pt] &\leq & 2(p-1)\int _\Omega u^\frac p2\mathbf{w}\cdot \nabla u^\frac p2+\alpha p \delta _1\int _\Omega u^p\quad \textrm{for all }t\in (0,{T_{\max}}), \end{eqnarray}

where $\delta _1\,:\!=\,\max \limits _{v\in [0, m]}F(v)$ is a positive constant due to the assumption (H1) (recall $m\gt 0$ is a constant given by (1.4)). By Lemma 2.5 and Young’s inequality, one can see that there exists a constant $C_1\gt 0$ satisfying

\begin{eqnarray*} 2(p-1)\int _\Omega u^\frac p2\mathbf{w}\cdot \nabla u^\frac p2&\leq & d_u (p-1)\frac 2{p}\int _\Omega \big|\nabla u^\frac p2\big|^2+2C_1\frac{p(p-1)}{d_u}\int _\Omega u^p \end{eqnarray*}

for all $t\in (0,{T_{\max}})$ . Inserting this inequality into (2.17), we obtain

(2.18) \begin{eqnarray} &&\frac d{dt}\int _\Omega u^p+p\mu \int _\Omega u^p+\frac{2d_u(p-1)}{p}\int _\Omega \big|\nabla u^\frac p2\big|^2\nonumber \\[5pt] &\leq & \left (2C_1\frac{p(p-1)}{d_u}+\alpha p \delta _1\right )\int _\Omega u^p\quad \textrm{for all }t\in (0,{T_{\max}}). \end{eqnarray}

By Lemma 2.2 and the Gagliardo–Nirenberg inequality, we know that

\begin{eqnarray*} \int _\Omega u^p&=&\big\|u^\frac p2\big\|_{L^{2}(\Omega )}^2\\[5pt] &\leq & C\!\left (\|\nabla u^\frac p2\|_{L^{2}(\Omega )}^{2\theta }\big\|u^\frac p2\big\|_{L^{\frac 2p}(\Omega )}^{2(1-\theta )}+\big\|u^\frac p2\big\|_{L^{\frac 2p}(\Omega )}^2\right )\\[5pt] &\leq & C(p)\|\nabla u^\frac p2\|_{L^{2}(\Omega )}^{2\theta }+C(p)\quad \textrm{for all }t\in (0,{T_{\max}}), \end{eqnarray*}

where $\theta =\frac{\frac p2-\frac 12}{\frac 1n+\frac p2-\frac 12}\in (0,1).$

Since $\theta \in (0,1)$ , it follows from Young’s inequality that

\begin{eqnarray*} \left (2C_1\frac{p(p-1)}{d_u}+\alpha p \delta _1\right )\int _\Omega u^p&\leq &\frac{d_u(p-1)}{p}\int _\Omega \big|\nabla u^\frac p2\big|^2+C(p) \end{eqnarray*}

$\textrm{for all }t\in (0,{T_{\max}})$ . Substituting this inequality into (2.18), we obtain

\begin{eqnarray*} \frac d{dt}\int _\Omega u^p+p\mu \int _\Omega u^p+\frac{d_u(p-1)}{p}\int _\Omega \big|\nabla u^\frac p2\big|^2\leq C(p)\quad \textrm{for all }t\in (0,{T_{\max}}). \end{eqnarray*}

Solving the above ordinary differential inequality, we get

\begin{equation*} \|u(\cdot,t)\|_{L^{p}(\Omega )}\leq C(p)\quad \text {for all }t\in \left(0, {T_{\max}}\right). \end{equation*}

Particularly, there exists a constant $C\gt 0$ independent of $p$ satisfying

\begin{equation*} \|u(\cdot,t)\|_{L^{2n}(\Omega )}\leq C\quad \text {for all }t\in \left(0, {T_{\max}}\right). \end{equation*}

By Duhamel’s principle, one has

\begin{equation*}u(\cdot,t)=e^{d_u t\Delta }u_0-\int ^t_0e^{d_u (t-s)\Delta }\nabla \cdot (u \mathbf {w} )ds+\int ^t_0e^{d_u (t-s)\Delta }\varphi (u,v)ds\quad \text {for all }t\in \left(0, {T_{\max}}\right),\end{equation*}

where $\varphi (u,v)=f(v)-\alpha uF(v)$ . By the $L^{p}$ - $L^{q}$ -estimates for the Neumann heat semigroup (cf. [52, Lemma 1.3]), we have

\begin{eqnarray*} \|u(\cdot,t)\|_{L^{\infty }(\Omega )}&\leq & C\|u_0\|_{L^\infty (\Omega )}+C\int ^t_0\left(1+(t-s)^{-\frac 12-\frac n2\cdot \frac 1{2n}}\right)e^{-d_u\rho _1 (t-s)}\|u\|_{L^{2n}(\Omega )}\|\mathbf{w}\|_{L^{\infty }(\Omega )} ds\\[5pt] &&+C\int ^t_0\left(1+(t-s)^{-\frac n2\cdot \frac 1{2n}}\right)e^{-d_u\rho _1 (t-s)}\|\varphi (u,v)\|_{L^{2n}(\Omega )}ds\\[5pt] &\leq & C\|u_0\|_{L^\infty (\Omega )}+C\int ^t_0\left(1+(t-s)^{-\frac 34}\right)e^{-d_u\rho _1 (t-s)}ds\\[5pt] &&+C\int ^t_0(1+(t-s)^{-\frac 14})e^{-d_u\rho _1 (t-s)}ds\\[5pt] &\leq & C\quad \textrm{for all }t\in (0,{T_{\max}}), \end{eqnarray*}

where $\rho _1 \gt 0$ denotes the first non-zero eigenvalue of $-\Delta$ in $\Omega$ under the Neumann boundary condition. This completes the proof.

In view of Lemmas 2.1, 2.5 and 2.6, we find $T_{\max}=\infty$ . We now derive a priori $L^\infty$ -estimate involving the gradient of $v$ and $\mathbf{w}$ .

Lemma 2.7. There exists a constant $C\gt 0$ such that

(2.19) \begin{equation} \|v(\cdot,t)\|_{W^{1, \infty }(\Omega )}+\|{\mathbf{w} }(\cdot,t)\|_{W^{1, \infty }(\Omega )}\leq C\quad \textit{for all }t\gt 0. \end{equation}

Proof. In view of Lemmas 2.5 and 2.6, we only need to prove

(2.20) \begin{equation} \|\nabla v(\cdot,t)\|_{L^{\infty }(\Omega )}+\|\nabla \mathbf{w}(\cdot,t)\|_{L^{\infty }(\Omega )}\leq C\quad \textrm{for all }t\gt 0. \end{equation}

Using the variation of constants representation of $v$ , we have

\begin{equation*} v(\cdot,t)=e^{d_v t\Delta }v_0+\int ^t_0e^{d_v (t-s)\Delta }\varphi (u(\cdot,s),v(\cdot,s))ds\quad \text {for all }t\gt 0, \end{equation*}

where $\varphi (u,v)=f(v)-\alpha uF(v)$ . It follows from (H1), (H2), Lemmas 2.1 and 2.6 that

\begin{equation*} \|\varphi (u(\cdot,t),v(\cdot,t))\|_{L^{\infty}(\Omega )} \leq C\quad \text {for all }t\gt 0, \end{equation*}

which together with the $L^{p}$ - $L^{q}$ -estimates for the Neumann heat semigroup (cf. [52, Lemma 1.3]) shows that

(2.21) \begin{eqnarray} \|\nabla v(\cdot,t)\|_{L^{\infty }(\Omega )}&\leq & C\|v_0\|_{W^{1,\infty }(\Omega )}+C\int ^t_0\left(1+(t-s)^{-\frac 12}\right)e^{-d_v\rho _1 (t-s)}\|\varphi \|_{L^{\infty}(\Omega )}ds\nonumber \\[5pt] &\leq & C\|v_0\|_{W^{1,\infty }(\Omega )}+C\int ^t_0\left(1+(t-s)^{-\frac 12}\right)e^{-d_v\rho _1 (t-s)}ds\nonumber \\[5pt] &\leq & C\quad \textrm{for all }t\gt 0, \end{eqnarray}

where $\rho _1 \gt 0$ denotes the first non-zero eigenvalue of $-\Delta$ in $\Omega$ under the Neumann boundary condition. Similarly, via the variation of constants formula of ${w_i}\ (i=1,2,\cdots,n)$ , we have

\begin{eqnarray*} \nabla{w_{i}}(t)=\nabla e^{d_w t\Delta }{w_{i}(\cdot,0)}+\gamma \int ^t_0 \nabla e^{d_w (t-s)\Delta }\partial _{x_i} v(\cdot,s)ds\quad \textrm{for all }t\gt 0. \end{eqnarray*}

By Lemma 2.4, for $i=1,2,\cdots,n$ , we have

\begin{eqnarray*} \|\nabla{w_{i}}(\cdot,t)\|_{L^{\infty }(\Omega )}&\leq & C\|{w_{i}(\cdot,0)}\|_{W^{1,\infty }(\Omega )}+C\int ^t_0\left(1+(t-s)^{-\frac 12}\right)e^{-d_w\lambda _1 (t-s)}\|\nabla v(\cdot,s)\|_{L^{\infty }(\Omega )} ds\\[5pt] &\leq & C\|{w_{i}(\cdot,0)}\|_{W^{1,\infty }(\Omega )}+C\int ^t_0\left(1+(t-s)^{-\frac 12}\right)e^{-d_w\lambda _1 (t-s)} ds\\[5pt] &\leq & C\quad \textrm{for all }t\gt 0. \end{eqnarray*}

This along with (2.21) proves (2.20).

Proof of Theorem 1.1. In veiw of Lemmas 2.1, 2.6 and 2.7, it remains to prove $\limsup \limits _{t\rightarrow +\infty }v(x,t)\leq K$ for all $x\in \bar \Omega$ . Based on the assumptions (H1)–(H4), this can be proved by using the strong maximum principle and the comparison principle, we omit the standard argument for brevity and refer readers to [17, Lemma 2.2].

3. Global stability

In this section, we are devoted to proving the global stability results stated in Theorems 1.2 and1.3 by constructing Lyapunov functionals. To this end, we first need some regularity results as follows.

Lemma 3.1. Let $(u, v, \mathbf{w} )$ be the unique global classical solution of (1.2), which is given by Theorem 1.1. Then for any $0\lt \theta \lt 1$ , there exists $C(\theta )\gt 0$ such that

(3.1) \begin{equation} \|u, v, \mathbf{w} \|_{C^{2+\theta, 1+\frac{\theta }{2}}(\bar{\Omega } \times [1, \infty ))} \leq C(\theta ). \end{equation}

Proof. This proof is based on a standard parabolic regularity for parabolic equations (see [46, Theorem 2.1] for instance). For the reader’s convenience, we shall sketch the proof below. We rewrite (1.2) as:

(3.2) \begin{equation} \left \{\begin{array}{l@{\quad}l} u_{t}-d_{u} \Delta u+\mathbf{w}\cdot \nabla u=H_1(u,v,\mathbf{w}), & x \in \Omega, t\gt 0, \\[5pt] v_{t}=d_{v} \Delta v+H_2(u,v), & x \in \Omega, t\gt 0, \\[5pt] \mathbf{w} _{t}=d_{w} \Delta \mathbf{w} +\gamma \nabla v, & x \in \Omega, t\gt 0, \\[5pt] \nabla u \cdot \mathbf n =\nabla v \cdot \mathbf n=0,\ \mathbf{w}=\mathbf 0,& x \in \partial \Omega, t\gt 0,\\[5pt] u(x,0)=u_0(x),\ v(x,0)=v_0(x)\ \mathbf{w}(x,0)=\mathbf{w}_0(x),& x \in \Omega, \end{array}\right. \end{equation}

where

\begin{equation*}H_1(u,v,\mathbf {w})=-u\nabla \cdot \mathbf {w} +\alpha u F(v)-\mu u\quad \text {and}\quad H_2(u,v)=f(v)-\alpha u F(v).\end{equation*}

By (2.16) and (2.19), one has

\begin{equation*}\|\nabla v, H_1(u,v,\mathbf {w}), H_2(u,v)\|_{L^\infty \left (\Omega \times \left (0,\infty \right )\right )}\leq C.\end{equation*}

For any $p\geq 1$ , applying the interior $L^{p}$ estimate [29, Theorems 7.30 and 7.35] to (3.2), we have

\begin{equation*}\|u,v,\mathbf {w}\|_{W^{2,1}_p\left (\Omega \times [1,\infty )\right )}\leq C.\end{equation*}

Taking $p$ appropriately large and using the Sobolev embedding theorem, we can find a positive constant $\theta \in (0,1)$ such that

\begin{equation*} \|u,v,\mathbf {w}\|_{C^{1+\theta, \frac {1+\theta }{2}}{\left (\bar {\Omega } \times [1,\infty )\right )}} \leq C. \end{equation*}

Then, it follows that

\begin{equation*}\|\nabla v, H_1(u,v,\mathbf {w}), H_2(u,v)\|_{C^{\theta, \frac {\theta }{2}}{\left (\bar {\Omega } \times [1,\infty )\right )}}\leq C.\end{equation*}

This along with an application of the interior Schauder estimate [23] to (3.2) shows that

\begin{equation*} \|u,v,\mathbf {w}\|_{C^{2+\theta, 1+\frac {\theta }{2}}(\bar {\Omega } \times [1, \infty ))} \leq C. \end{equation*}

This completes the proof.

To proceed, we recall two basic results.

Lemma 3.2. ([47, Lemma 1.1]) Let $\tau \geq 0, c\gt 0$ be constants, $\psi (t) \geq 0, \int _{\tau }^{\infty } \rho (t) \textrm{d} t\lt \infty$ . Assume that $\varphi \in C^{1}([\tau, \infty )), \varphi$ is bounded from below and satisfies

\begin{equation*} \varphi ^{\prime }(t) \leq -c \psi (t)+\rho (t) \quad\ \textit{in}\ [\tau, \infty ). \end{equation*}

If $\psi \in C^{1}([\tau, \infty ))$ and $\psi ^{\prime }(t) \leq k$ in $[\tau, \infty )$ for some $k\gt 0,$ or $\psi \in C^{\alpha }([\tau, \infty ))$ and $\|\psi \|_{C^{\alpha }([\tau, \infty ))} \leq k$ for some constants $0\lt \alpha \lt 1$ and $k\gt 0$ , then

\begin{equation*}\lim _{t \rightarrow \infty } \psi (t)=0.\end{equation*}

Lemma 3.3. Let $F$ satisfy the conditions in (H1) and (H3) and $(u, v,\mathbf{w} )$ be a solution of (1.2). Define a function for some $\xi \gt 0:$

\begin{equation*} \zeta (v)=\int _{\xi }^{v} \frac {F(s)-F(\xi )}{F(s)} d s. \end{equation*}

Then $\zeta (v)$ is a convex function such that $\zeta (v) \geq 0 .$ If we further assume that $v \rightarrow \xi$ as $t \rightarrow \infty$ , then there is a constant $T_2\gt 0$ such that for all $t \geq T_2$ it holds that

\begin{equation*} \frac {F^{\prime }(\xi )}{4 F(\xi )}(v-\xi )^{2} \leq \zeta (v)=\int _{\xi }^{v} \frac {F(s)-F(\xi )}{F(s)} d s \leq \frac {F^{\prime }(\xi )}{F(\xi )}(v-\xi )^{2}. \end{equation*}

Proof. The result immediately follows from the Taylor expansion of $\zeta (v)$ .

3.1. Global stability of the coexistence steady state: Proof of Theorem 1.2

Lemma 3.4. Let $\alpha F(K)\gt \mu$ and the conditions in Theorem 1.2 hold, if

(3.3) \begin{equation} d_w\gt \left (\frac{u_*}{4d_u}+\frac{\alpha \gamma ^2}{d_v c_0\mu }\right )\frac{C_P^2(\Omega )}2, \end{equation}

then the two functions $\mathcal E_1(t)$ and $\mathcal F_1(t)$ defined by:

\begin{eqnarray*} \mathcal E_1(t)&=&\int _\Omega \!\left (u-{u_*}-{u_*}\ln \frac u{u_*}\right )+\int _\Omega \int ^v_{v_*}\frac{F(s)-F({v_*})}{F(s)}ds+\int _\Omega |{\mathbf{w} }|^2,\quad t\gt 0,\\[5pt] \mathcal F_1(t)&=&\int _\Omega \!\left (\left |\nabla u\right |^2+|\nabla v|^2+| \mathbf{w} |^2\right )+\int _\Omega (v-{v_*})^2,\quad t\gt 0, \end{eqnarray*}

satisfy

(3.4) \begin{equation} \mathcal E_1^\prime (t) \leq -\varepsilon _1 \mathcal F_1 (t) \quad \textit{for all }t\gt{T_0}, \end{equation}

for some constant $\varepsilon _1\gt 0$ , where $T_0$ , $c_0$ and $C_P(\Omega )$ are given by (1.7), (1.8) and (1.9), respectively.

Proof. By (1.6), (2.16) and Theorem 1.1, with some $\beta \in (0,d_uu_*)$ to be specified below, for all $t\gt 0$ we have

(3.5) \begin{eqnarray} \frac d{dt}\int _\Omega \!\left (u-{u_*}-{u_*}\ln \frac u{u_*}\right ) &=&\int _\Omega \!\left (1-\frac{u_*} u\right )\left (d_{u} \Delta u-\nabla \cdot (u \mathbf{w} )+\alpha uF(v)-\mu u\right )\nonumber \\[5pt] &=&-d_u{u_*}\int _\Omega \frac{|\nabla u|^2}{u^2}+{u_*}\int _\Omega{\mathbf{w}\cdot \frac{\nabla u}u}+\int _\Omega \!\left (\alpha F(v)-\mu \right )(u-{u_*})\nonumber \\[5pt] &\leq &-\beta \int _\Omega \frac{|\nabla u|^2}{u^2}+\frac{{u_*}^2}{4(d_u{u_*}-\beta )}\int _\Omega \!\left |\mathbf{w}\right |^2+\alpha \int _\Omega \!\left (F(v)-F({v_*})\right )(u-{u_*})\nonumber \\[5pt] &\leq &-\frac \beta{M^2}\int _\Omega |\nabla u|^2+\frac{{u_*}^2}{4(d_u{u_*}-\beta )}\int _\Omega \!\left |\mathbf{w}\right |^2+\alpha \int _\Omega \!\left (F(v)-F({v_*})\right )(u-{u_*}),\nonumber\\ \end{eqnarray}

where $M$ is given by (2.16). Moreover, for all $t\gt 0$ we have

(3.6) \begin{eqnarray} \frac d{dt}\int _\Omega \int ^v_{v_*}\frac{F(s)-F({v_*})}{F(s)}ds &=&\int _\Omega \frac{F(v)-F({v_*})}{F(v)}\left (d_{v} \Delta v+f(v)-\alpha uF(v)\right )\nonumber \\[3pt] &=&-d_v F({v_*})\int _\Omega \Psi (v)\left |\nabla v\right |^2+\int _\Omega (F(v)-F({v_*}))\left (\Phi (v)-\alpha u\right )\nonumber \\[3pt] &=&-d_v F({v_*})\int _\Omega \Psi (v)\left |\nabla v\right |^2-\alpha \int _\Omega (F(v)-F({v_*}))(u-{u_*})\nonumber \\[3pt] &&+\int _\Omega (F(v)-F({v_*}))\left (\Phi (v)-\Phi (v_*)\right ), \end{eqnarray}

where we have used (1.6) in the last inequality. For any $\delta \gt 0$ , by Young’s inequality and (1.9) we have

(3.7) \begin{eqnarray} \frac d{dt}\int _\Omega |\mathbf{w}|^2&=&2d_w\int _\Omega \Delta \mathbf{w}\cdot \mathbf{w}+2\gamma \int _\Omega \mathbf{w}\cdot \nabla v\nonumber \\[3pt] &=&-2d_w\sum ^{n}_{i=1}\int _\Omega |\nabla w_i|^2+2\gamma \int _\Omega \mathbf{w}\cdot \nabla v\nonumber \\[3pt] &\leq &-2d_w\int _\Omega |\nabla \mathbf{w}|^2+\delta \int _\Omega |\nabla v|^2+\frac{\gamma ^2}{\delta }\int _\Omega |\mathbf{w}|^2\nonumber \\[3pt] &\leq &-\left (\frac{2d_w}{C_P^2(\Omega )}-\frac{\gamma ^2}{\delta }\right )\int _\Omega |\mathbf{w}|^2+\delta \int _\Omega |\nabla v|^2\quad \textrm{for all }t\gt 0. \end{eqnarray}

By (H1), (H3), (H4), and the mean value theorem, we can find $\eta _1, \eta _2\in (0,m)$ such that

(3.8) \begin{eqnarray} \int _\Omega (F(v)-F({v_*}))\left (\Phi (v)-\Phi (v_*)\right )=\int _\Omega F^{\prime}(\eta _1) \Phi ^{\prime}(\eta _2)(v-{v_*})^2\leq 0. \end{eqnarray}

Using (1.8) and (3.5)–(3.8), we obtain

\begin{eqnarray*} \mathcal E_1^\prime (t)&\leq &-\frac \beta{M^2}\int _\Omega |\nabla u|^2-\left (\frac{2d_w}{C_P^2(\Omega )}-\left (\frac{{u_*}^2}{4(d_u{u_*}-\beta )}+\frac{\gamma ^2}{\delta }\right )\right )\int _\Omega | \mathbf{w}|^2\\[5pt] &&-\left (d_v c_0F({v_*})-\delta \right )\int _\Omega \!\left |\nabla v\right |^2+\int _\Omega F^{\prime}(\eta _1) \Phi ^{\prime}(\eta _2)(v-{v_*})^2{\quad \textrm{for all }t\gt T_0}. \end{eqnarray*}

By continuity and (1.6), we know that if (3.3) holds, then we can pick $\beta \in (0,d_uu_*)$ small enough and $\delta \in (0,d_v c_0F({v_*}))$ closing to $d_v c_0F({v_*})$ such that

\begin{equation*} d_v c_0F({v_*})-\delta \gt 0\quad \text {and}\quad {\frac {2d_w}{C_P^2(\Omega )}-\left (\frac {{u_*}^2}{4(d_u{u_*}-\beta )}+\frac {\gamma ^2}{\delta }\right )}\gt 0. \end{equation*}

Therefore, we can arrive at (3.4) by taking

\begin{equation*}\varepsilon _1\,:\!=\,\min \!\left \{\frac {\beta }{M^2},\frac {2d_w}{C_P^2(\Omega )}-\left (\frac {{u_*}^2}{4(d_u{u_*}-\beta )}+\frac {\gamma ^2}{\delta }\right ), d_v c_0F({v_*})-\delta, -F^{\prime}(\eta _1) \Phi ^{\prime}(\eta _2)\right \}.\end{equation*}

Thus, the proof is completed.

Lemma 3.5. Let $\alpha F(K)\gt \mu$ and the conditions in Theorem 1.2 hold, for any $0\lt \theta \lt 1$ we have

(3.9) \begin{equation} \|u-{u_*}\|_{C^{2+\theta }(\bar{\Omega })}+\|v-{v_*}\|_{C^{2+\theta }(\bar{\Omega })} +\|\mathbf{w} \|_{C^{2+\theta }(\bar{\Omega })}\rightarrow 0 \textit{ as } t \rightarrow \infty. \end{equation}

Proof. With Lemma 3.1 at hand, the conclusion can be proved by a similar argument as in [46, Lemma 3.4]. For the reader’s convenience, we shall sketch the proof here. Let $\mathcal E_1 (t), \mathcal F_1 (t)$ be given by Lemma 3.4. For $\theta \in (0,1)$ , by (3.1) we have $\mathcal F_1 (t)\geq 0$ , $\mathcal F_1 (t) \in C^{\theta/ 2}([1, \infty ))$ and $\|\mathcal F_1 \|_{C^{\theta/ 2}([1, \infty ))} \leq k$ for some $k\gt 0$ . Clearly, $\mathcal E_1 (t) \in C^1([1, \infty ))$ . Using Taylor’s expansion, we can find $\tilde{v}$ between $v$ and $v_*$ such that

\begin{equation*}\int ^v_{v_*}\frac {F(s)-F({v_*})}{F(s)}ds=\frac {F({v_*}) F^{\prime }(\tilde {v})}{2 F^{2}(\tilde {v})}(v-{v_*})^{2} \geq 0,\end{equation*}

and $\tilde{u}$ between $u$ and $u_*$ such that

\begin{equation*}u-{u_*}-{u_*}\ln \frac u {u_*}=\frac {u_*}{2\tilde {u}^2}(u-{u_*})^2\geq 0.\end{equation*}

Therefore, $\mathcal E_1(t)$ is bounded from blow in $[1,\infty )$ with $\mathcal E_1(t)\geq 0$ for $t\in [1,\infty )$ . We now apply Lemma 3.2 to (3.4) and conclude that $\lim _{t \rightarrow \infty }\mathcal F_1(t)=0$ , which gives

(3.10) \begin{equation} \lim _{t \rightarrow \infty }\left (\|\nabla u\|_{L^2(\Omega )}+\|v-{v_*}\|_{L^2(\Omega )}+\|\mathbf{w}\|_{L^2(\Omega )}\right )=0. \end{equation}

Taking $0\lt \theta \lt{\theta ^{\prime}}\lt 1$ , by Lemma 3.1 we have

\begin{equation*} \|u, v, \mathbf {w}\|_{C^{2+\theta ^{\prime}, 1+\frac {\theta ^{\prime}}{2}}(\bar {\Omega } \times [1, \infty ))} \leq C(\theta ^{\prime}). \end{equation*}

This alongside the compact arguments and uniqueness of limits (cf. [46, (3.12)], see also [16, Remark 6.1]) shows that

(3.11) \begin{equation} \lim _{t \rightarrow \infty }\left (\|v-{v_*}\|_{C^{2+{\theta }}(\bar{\Omega })}+\|\mathbf{w}\|_{C^{2+{\theta }}(\bar{\Omega })}\right )=0. \end{equation}

It remains to prove

(3.12) \begin{equation} \|u-{u_*}\|_{C^{2+{\theta }}(\bar{\Omega })}\rightarrow 0 \textrm{ as } t \rightarrow \infty. \end{equation}

By the second equation of (1.2) and (1.6), we obtain

(3.13) \begin{eqnarray} \frac d{dt}\left (\frac 1{|\Omega |}\int _\Omega vdx\right )=\bar{v}^{\prime }(t) &=& \frac{1}{|\Omega |} \int _{\Omega }(f(v)-\alpha u F(v)) \textrm{d} x \nonumber \\[5pt] &=& \frac{1}{|\Omega |} \int _{\Omega }\left (f(v)-f\!\left ({v_*}\right )\right ) \textrm{d} x-\frac{\alpha }{|\Omega |} \int _{\Omega } u(F(v)-F\!\left ({v_*}\right )) \textrm{d} x \nonumber \\[5pt] &&-\frac{\alpha F\!\left ({v_*}\right )}{|\Omega |} \int _{\Omega }\!\left (u-{u_*}\right ) \textrm{d} x \nonumber \\[5pt] &\,=\!:\,& J_1 (t)+J_2 (t)+J_{3}(t)\quad \textrm{for all }t\gt 0, \end{eqnarray}

where $\bar{z}\,:\!=\,\frac{1}{|\Omega |} \int _{\Omega } z \textrm{d} x$ for $z \in L^{1}(\Omega )$ . It follows from (H1), (H3) and (3.11) that $J_1 (t)\rightarrow 0$ and $J_2 (t)\rightarrow 0$ as $t \rightarrow \infty$ . By (3.1) we know that $\|\bar{v}^{\prime }\|_{C^{{\theta }/ 2}([1, \infty ))} \leq C_2({\theta })$ , which together with (3.11) shows that $\bar{v}^{\prime }(t) \rightarrow 0$ as $t \rightarrow \infty$ . Therefore, we can infer from (3.13) that $J_{3}(t)\rightarrow 0$ as $t \rightarrow \infty$ , that is,

\begin{equation*} \bar {u}(t) \rightarrow {u_*} \text { as } t \rightarrow \infty. \end{equation*}

This yields, thanks to the Poincaré inequality and (3.10), that

\begin{align*} \|u-{u_*}\|_{L^2(\Omega )} \leq \|u-\bar{u}\|_{L^2(\Omega )}+\|\bar{u}-{u_*}\|_{L^2(\Omega )}\leq C\|\nabla u\|_2 +\|\bar{u}-{u_*}\|_{L^2(\Omega )}\rightarrow 0 \textrm{ as } t \rightarrow \infty. \end{align*}

Similar to the proof of (3.11), again by the compact arguments and uniqueness of limits, we have (3.12). This completes the proof.

Proof of Theorem 1.2. In view of Lemma 3.5, we obtain Theorem 1.2 immediately.

3.2. Global stability of the prey-only steady state: Proof of Theorem 1.3

In this subsection, we shall consider the stability of the prey-only steady state $\left (0, 1, \textbf{0}\right )$ which is given by (1.5) in the case of weak predation $\alpha F(K)\leq \mu$ . We shall show that $\left (0, 1, \textbf{0}\right )$ is globally asymptotically stable. Moreover, we shall establish the exponential stability of $\left (0, 1, \textbf{0}\right )$ under the condition of $\alpha F(K)\lt \mu$ . To this end, we construct an appropriate Lyapunov functional as below.

Lemma 3.6. Assume that $\alpha F(K)\leq \mu$ and the conditions in Theorem 1.3 hold. Let $D=\frac{2\gamma ^2 C_P^2(\Omega )}{d_vd_wF(K)c_0}$ . Then functions $\mathcal E_2(t)$ and $\mathcal F_2(t)$ defined by:

\begin{eqnarray*} \mathcal E_2(t)&=& D\int _\Omega u+D\int _\Omega \int ^v_K\frac{F(s)-F(K)}{F(s)}ds+\int _\Omega |\mathbf{w} |^2,\quad t\gt 0,\\[5pt] \mathcal F_2(t)&=& \int _\Omega \!\left (\left |\nabla v\right |^2+ (v-K)^2+|\mathbf{w} |^2\right ),\quad t\gt 0, \end{eqnarray*}

satisfy

(3.14) \begin{equation} \mathcal E_2^\prime (t) \leq -\varepsilon _2 \mathcal F_2 (t)-D(\mu -\alpha F(K))\int _\Omega u \quad \textit{for all }t\gt{T_0}, \end{equation}

for some constant $\varepsilon _2\gt 0$ , where $T_0$ is given by (1.7).

Proof. By the same argument as in the proof of Lemma 3.4, with some calculations, we can find $\eta _1, \eta _2$ between 1 and $v$ such that

(3.15) \begin{eqnarray} &&\frac d{dt}\int _\Omega \int ^v_{K}\frac{F(s)-F(K)}{F(s)}ds\nonumber \\[5pt] &=&-d_v F(K)\int _\Omega \Psi (v)\left |\nabla v\right |^2+\int _\Omega (F(v)-F(K))\left (\Phi (v)-\alpha u\right )\nonumber \\[5pt] &=&-d_v F(K)\int _\Omega \Psi (v)\left |\nabla v\right |^2-\alpha \int _\Omega u(F(v)-F(K))+\int _\Omega (F(v)-F(K))(\Phi (v)-\Phi (K))\nonumber \\[5pt] &=&-d_v F(K)\int _\Omega \Psi (v)\left |\nabla v\right |^2-\alpha \int _\Omega u(F(v)-F(K))+F^{\prime}(\eta _1) \Phi ^{\prime}(\eta _2)\int _\Omega (v-K)^2 \end{eqnarray}

for all $t\gt 0$ . Moreover, the integration of the first equation of (1.2) along with boundary conditions yields

\begin{eqnarray*} \int _\Omega u_t&=& \alpha \int _\Omega u(F(v)-F(K))+\int _\Omega u(\alpha F(K)-\mu )\quad \textrm{for all }t\gt 0. \end{eqnarray*}

Taking $\delta =\frac{\gamma ^2 C_P^2(\Omega )}{d_w}$ in (3.7) and using (3.15), we have

(3.16) \begin{eqnarray} \mathcal E_2^\prime (t)&\leq &-\int _\Omega \!\left ({d_v DF(K)}\Psi (v)-\delta \right )\left |\nabla v\right |^2+DF^{\prime}(\eta _1) \Phi ^{\prime}(\eta _2)\int _\Omega (v-K)^2\nonumber \\[5pt] &&-{\frac{d_w}{C_P^2(\Omega )}}\int _\Omega |\mathbf{w}|^2-D\int _\Omega (\mu -\alpha F(K)) u \quad \textrm{for all }t\gt 0. \end{eqnarray}

It follows from (H1), (H3), (H4) and (1.8) that $DF^{\prime}(\eta _1) \Phi ^{\prime}(\eta _2)\lt 0$ and

\begin{equation*}{{d_v DF(K)}\Psi (v)-\delta }=\frac {2\gamma ^2 C_P^2(\Omega )}{d_w}\cdot \frac {\Psi (v)}{c_0}-\frac {\gamma ^2 C_P^2(\Omega )}{d_w}\geq \frac {\gamma ^2 C_P^2(\Omega )}{d_w} {\quad \text {for all }t\gt T_0}.\end{equation*}

Therefore, (3.14) is a direct consequence of (3.16) with

\begin{equation*} \varepsilon _2\,:\!=\,\min \!\left \{\frac{\gamma ^2 C_P^2(\Omega )}{d_w}, -DF^{\prime}(\eta _1) \Phi ^{\prime}(\eta _2),\frac{d_w}{C_P^2(\Omega )}\right \}. \end{equation*}

This completes the proof.

Lemma 3.7. Assume that $\alpha F(K)\leq \mu$ and the conditions in Theorem 1.3 hold. For any $0\lt \theta \lt 1$ , we have

\begin{equation*} \|u\|_{C^{2+\theta }(\bar{\Omega })}+\|\mathbf{w} \|_{C^{2+\theta }(\bar{\Omega })}+\|v-K\|_{C^{2+\theta }(\bar{\Omega })} \rightarrow 0 \textit{ as } t \rightarrow \infty. \end{equation*}

Proof. The proof is similar to that of Lemma 3.5, hence we omit it here for brevity.

Lemma 3.8. Assume that $\alpha F(K)\lt \mu$ and the conditions in Theorem 1.3 hold. Then there exist positive constants $C$ , $\lambda$ and $T_1\geq 1$ such that

\begin{equation*} \|u\|_{L^\infty (\Omega )}+\|v-K\|_{L^\infty (\Omega )}+\|\mathbf{w} \|_{L^\infty (\Omega )} \leq C e^{-\lambda t}\quad \textit{for all } t\gt T_1. \end{equation*}

Proof. According to (3.9) and Lemma 3.3, we can find $T_1\geq 1$ such that

(3.17) \begin{equation} \frac{F^{\prime }(K)}{4 F(K)}\int _\Omega (v-K)^{2} \leq \int _\Omega \int _{K}^{v} \frac{F(s)-F(K)}{F(s)} d s \leq \frac{F^{\prime }(K)}{F(K)}\int _\Omega (v-K)^{2}\quad \textrm{for all }t\gt T_1. \end{equation}

Using the right inequality in (3.17), the definition of $\mathcal E_2(t)$ and $\mathcal F_2(t)$ , (3.14) and $\alpha F(K)\lt \mu$ , we can find two positive constants $C_1$ and $C_2$ such that

\begin{equation*} \mathcal E_2^\prime (t) \leq -C_1\!\left ( \mathcal F_2 (t)+\int _\Omega u\right )\leq -C_2 \mathcal E_2 (t)\quad \textrm{for all }t\gt T_1. \end{equation*}

Thus, there is a constant $C_3\gt 0$ such that

\begin{equation*}\mathcal E_2 (t) \leq C_{3} e^{-C_2 t}\quad \text {for all }t\gt T_1.\end{equation*}

This together with the definition of $\mathcal E_2(t)$ and the left inequality in (3.17) shows that

(3.18) \begin{equation} \|u\|_{L^1(\Omega )} +\|\mathbf{w}\|^2_{L^2(\Omega )}+\|v-K\|^2_{L^2(\Omega )}\leq C_{3} e^{-C_2 t}\quad \textrm{ for all } t\gt T_1. \end{equation}

We shall extend this result to the estimates of $L^{\infty }$ -norm. Indeed (3.9) implies that

\begin{equation*}\|u,v,\mathbf {w}\|_{W^{1,\infty }(\Omega )}\leq C\quad \text {for all }t\geq 1.\end{equation*}

This together with (3.18) and the Gagliardo–Nirenberg inequality for any $\psi \in W^{1, \infty }(\Omega )$ :

\begin{eqnarray*} \|\psi \|_{L^\infty (\Omega )} \leq C\|\psi \|_{W^{1, \infty }(\Omega )}^{\frac{n}{n+1}}\|\psi \|_{L^1(\Omega )}^{\frac{1}{n+1}}, \ \ \|\psi \|_{L^\infty (\Omega )} \leq C\|\psi \|_{W^{1, \infty }(\Omega )}^{\frac{n}{n+2}}\|\psi \|_{L^2(\Omega )}^{\frac{2}{n+2}}, \end{eqnarray*}

yields the following decay estimate for any $t\gt T_1$ (recalling $T_1\geq 1$ ):

\begin{eqnarray*} \left \{ \begin{array}{l@{\quad}l} \|u\|_{L^\infty (\Omega )} \leq C\|u\|_{W^{1, \infty }(\Omega )}^{\frac{n}{n+1}}\|u\|_{L^1(\Omega )}^{\frac{1}{n+1}} \leq C e^{-\frac{C_2 t}{n+1}},\\[5pt] \|v-K\|_{L^\infty (\Omega )} \leq C\|v-K\|_{W^{1, \infty }(\Omega )}^{\frac{n}{n+2}}\|v-K\|_{L^2(\Omega )}^{\frac{2}{n+2}} \leq C e^{-\frac{C_2 t}{n+2}},\\[5pt] \|\mathbf{w} \|_{L^\infty (\Omega )}\leq C\|\mathbf{w} \|_{W^{1, \infty }(\Omega )}^{\frac{n}{n+2}}\|\mathbf{w} \|_{L^2(\Omega )}^{\frac{2}{n+2}} \leq C e^{-\frac{C_2 t}{n+2}}. \end{array}\right. \end{eqnarray*}

This completes the proof by defining $\lambda =-\frac{C_2}{n+2}$ .

Proof of Theorem 1.3. Using Lemmas 3.7 and 3.8, we get Theorem 1.3 immediately.

4. Applications and spatiotemporal patterns

There are two main purposes in this section. The first is to apply Theorems 1.1, 1.2 and 1.3 to two specific trophic functions: Holling type I (i.e., Lotka–Volterra): $F(v)=v$ and Holling type II: $F(v)=\frac{v}{1+v}$ , and restate the results more explicitly. The second is to investigate whether the model (1.2) can generate spatial heterogeneous patterns comparable with experimental observations. To this end, we shall conduct the linear instability analysis to identify the possible parameter regimes of pattern formation and then use numerical simulations to illustrate the patterns. Throughout this section, for the sake of brevity, we shall take the carrying capacity $K=1$ and consider the growth function $f(v)$ in the case of the logistic type: $f(v)=v(1-v)$ .

4.1. Examples

The first example is the system (1.2) with the Lotka–Volterra type predator–prey interaction:

(4.1) \begin{equation} \left \{\begin{array}{l@{\quad}l} u_{t}=d_{u} \Delta u-\nabla \cdot (u \mathbf{w} )+\alpha u v-\mu u, & x \in \Omega, t\gt 0, \\[5pt] v_{t}=d_{v} \Delta v+v(1-v)-\alpha u v, & x \in \Omega, t\gt 0, \\[5pt] \mathbf{w} _{t}=d_{w} \Delta \mathbf{w} +\gamma \nabla v, & x \in \Omega, t\gt 0, \\[5pt] \nabla v \cdot \mathbf n =\nabla u \cdot \mathbf n =0, \mathbf{w}=\textbf{0},& x \in \partial \Omega, t\gt 0,\\[5pt] u(x,0)=u_0(x),\ v(x,0)=v_0(x),\ \mathbf{w}(x,0)=\mathbf{w}_0(x),& x \in \Omega. \end{array}\right. \end{equation}

Then, an application of Theorems 1.2 and 1.3 yields the following results on system (4.1).

Corollary 4.1. Let $\Omega \subset \mathbb{R}^n (n\geq 1)$ be a bounded domain with smooth boundary and assume the initial data $(u_0, v_0, \mathbf{w}_0)$ satisfy (1.3). Then the problem (4.1) has a unique global classical solution $(u,v,\mathbf{w} )$ satisfying

\begin{eqnarray*} \left \{\begin{array}{l@{\quad}c@{\quad}l} u, v \geq 0,\ u, v \in C(\bar{\Omega } \times [0, \infty )) \cap C^{2,1}(\bar{\Omega } \times (0, \infty )),\\[5pt] \mathbf{w} \in \left [C(\bar{\Omega } \times [0, \infty )) \cap C^{2,1}(\bar{\Omega } \times (0, \infty ))\right ]^n \end{array}\right. \end{eqnarray*}

with the following asymptotics:

1. (i) If $\alpha \leq \mu$ , then

   \begin{equation*} \|u\|_{L^\infty (\Omega )}+\|v-1\|_{L^\infty (\Omega )}+\|\mathbf {w} \|_{L^\infty (\Omega )} \rightarrow 0\textit { as }t\rightarrow \infty, \end{equation*}
   where the convergence is exponential if $\alpha \lt \mu$ .

2. (ii) If $\alpha \gt \mu$ and

   \begin{eqnarray*} d_w\gt \left (\frac{\alpha -\mu }{4\alpha ^2 d_u}+\frac{\alpha m^2\gamma ^2}{d_v \mu }\right )\frac{C_P^2(\Omega )}2, \end{eqnarray*}
   then
   \begin{equation*} \|u-u_*\|_{L^\infty (\Omega )}+\|v-v_*\|_{L^\infty (\Omega )}+\|\mathbf {w} \|_{L^\infty (\Omega )} \rightarrow 0\textit { as }t\rightarrow \infty, \end{equation*}
   where $m$ and $C_P(\Omega )$ are given by (1.4) and (1.9), respectively, and $(u_*,v_*)=\left (\frac{\alpha -\mu }{\alpha ^2},\frac \mu \alpha \right ).$

The second example is Holling type II predator–prey interaction, which specifies (1.2) as:

(4.2) \begin{equation} \left \{\begin{array}{l@{\quad}l} u_{t}=d_{u} \Delta u-\nabla \cdot (u \mathbf{w} )+\alpha \frac{uv}{1+v}-\mu u, & x \in \Omega, t\gt 0, \\[5pt] v_{t}=d_{v} \Delta v+v(1-v)-\alpha \frac{uv}{1+v}, & x \in \Omega, t\gt 0, \\[5pt] \mathbf{w} _{t}=d_{w} \Delta \mathbf{w} +\gamma \nabla v, & x \in \Omega, t\gt 0, \\[5pt] \nabla v \cdot \mathbf n =\nabla u \cdot \mathbf n =0, \mathbf{w}=\textbf{0},& x \in \partial \Omega, t\gt 0,\\[5pt] u(x,0)=u_0(x),\ v(x,0)=v_0(x),\ \mathbf{w}(x,0)=\mathbf{w}_0(x),& x \in \Omega. \end{array}\right. \end{equation}

Then the results of Theorems 1.2 and 1.3 imply the following results.

Corollary 4.2. Let $\Omega \subset \mathbb{R}^n (n\geq 1)$ be a bounded domain with smooth boundary and assume the initial data $(u_0, v_0, \mathbf{w}_0)$ satisfy (1.3). Then the problem (4.2) has a unique global classical solution $(u,v,\mathbf{w} )$ satisfying

\begin{eqnarray*} \left \{\begin{array}{l@{\quad}c@{\quad}l} u, v \geq 0,\ u, v \in C(\bar{\Omega } \times [0, \infty )) \cap C^{2,1}(\bar{\Omega } \times (0, \infty )),\\[5pt] \mathbf{w} \in \left [C(\bar{\Omega } \times [0, \infty )) \cap C^{2,1}(\bar{\Omega } \times (0, \infty ))\right ]^n, \end{array}\right. \end{eqnarray*}

and the following asymptotic behaviours hold:

1. (i) If $\alpha \leq 2\mu$ , then

   \begin{equation*} \|u\|_{L^\infty (\Omega )}+\|v-1\|_{L^\infty (\Omega )}+\|\mathbf {w} \|_{L^\infty (\Omega )} \rightarrow 0\textit { as }t\rightarrow \infty, \end{equation*}
   where the convergence is exponential if $\alpha \lt 2\mu$ .

2. (i) If $\alpha \gt 2\mu$ and

   \begin{eqnarray*} d_w\gt \left (\frac{\alpha -2\mu }{4d_u(\alpha -\mu )^2}+\frac{\alpha m^2\gamma ^2}{d_v \mu }\right )\frac{C_P^2(\Omega )}2, \end{eqnarray*}
   then
   \begin{equation*} \|u-u_*\|_{L^\infty (\Omega )}+\|v-v_*\|_{L^\infty (\Omega )}+\|\mathbf {w} \|_{L^\infty (\Omega )} \rightarrow 0\textit { as }t\rightarrow \infty, \end{equation*}
   where $m$ and $C_P(\Omega )$ are given by (1.4) and (1.9), respectively, and $(u_*,v_*)=\left (\frac{\alpha -2\mu }{(\alpha -\mu )^2},\frac \mu{\alpha -\mu }\right ).$

4.2. Linear instability analysis

In this subsection, we shall study the pattern formation possibly generated by the system (1.2). For this purpose, let us begin with the corresponding ODE system of (1.2). In this case, the third equation of (1.2) becomes $\mathbf{w}_t=\textbf{0}$ which together with $\left .\mathbf{w}\right |_{\partial \Omega }=\textbf{0}$ indicates that $\textbf{0}$ is the only possible steady state for $\mathbf{w}$ , and the first two equations of (1.2) in the absence of spatial effects become

(4.3) \begin{equation} \left \{\begin{array}{l@{\quad}l} u_{t}=\alpha u F(v)-\mu u, \\[5pt] v_{t}=v(1-v)-\alpha u F(v). \end{array}\right. \end{equation}

Clearly, (4.3) has three possible equilibria $(u_s, v_s)$ : $(0,0),\left (0, 1\right ),\left (u_*,v_*\right )$ , where $(u_*,v_*)$ is given by (1.6). Let $\boldsymbol{J}$ and $J_i$ $(i=1,2,3,4)$ be defined by:

(4.4) \begin{equation} \boldsymbol{J}=\left ( \begin{array}{c@{\quad}c} J_1 &{J_2} \\[5pt]{J_3} &{J_4} \\[5pt] \end{array} \right )=\left ( \begin{array}{c@{\quad}c} \alpha F(v_s)-\mu &\alpha u_s F^{\prime}(v_s)\\[5pt] -\alpha F(v_s) & 1-2v_s-\alpha u_s F^{\prime}(v_s) \\[5pt] \end{array} \right ). \end{equation}

The eigenvalue of $\boldsymbol{J}$ , denoted by $\rho$ , satisfies the following equation:

(4.5) \begin{equation} \rho ^2-(J_1+J_4)\rho +J_1J_4-J_2J_3=0. \end{equation}

By the Routh–Hurwitz criterion (cf. [35, Appendix B]) for second-order polynomials, we know that $(u_s, v_s)$ is linearly stable if and only if

\begin{equation*}-(J_1+J_4)\gt 0 \quad \text {and}\quad J_1J_4-J_2J_3\gt 0.\end{equation*}

Therefore, one can check that $(0,0)$ is linearly unstable for both Lotka–Volterra type and Holling type II predator–prey models, since $J_1J_4-J_2J_3=-\mu \lt 0$ . The steady state $(0,1)$ is linearly stable for $F(v)=v$ when $\alpha \lt \mu$ since $-(J_1+J_4)=1+\mu -\alpha \gt 0$ and $J_1J_4-J_2J_3=\mu -\alpha \gt 0$ and also linearly stable for $F(v)=\frac v{1+v}$ when $\alpha \lt 2\mu$ since $-(J_1+J_4)=1+\mu -\frac \alpha 2\gt 0$ and $J_1J_4-J_2J_3=\mu -\frac \alpha 2\gt 0$ . The homogeneous coexistence steady state

(4.6) \begin{align} (u_*,v_*)= \begin{cases} \left (\dfrac{\alpha -\mu }{\alpha ^{2}}, \dfrac{\mu }{\alpha }\right ),\quad &\textrm{if } F(v)=v,\\[12pt] \left (\dfrac{\alpha -2\mu }{(\alpha -\mu )^2},\dfrac \mu{\alpha -\mu }\right ),\quad &\textrm{if } F(v)=\dfrac{v}{1+v}, \end{cases} \end{align}

is linearly stable for $F(v)=v$ when $\alpha \gt \mu$ since $-(J_1+J_4)=\frac \mu \alpha \gt 0$ and $J_1J_4-J_2J_3=\frac{\mu (\alpha -\mu )}{\alpha }\gt 0$ and also linearly stable for $F(v)=\frac v{1+v}$ when $\alpha \gt 2\mu$ since $-(J_1+J_4)=\frac{2\mu ^2}{\alpha (\alpha -\mu )}\gt 0$ and $J_1J_4-J_2J_3=\frac{\mu (\alpha -2 \mu )}{\alpha }\gt 0$ .

Next, we proceed to consider the stability of $(0,1,\textbf{0})$ and $(u_*,v_*,\textbf{0})$ in the presence of spatial structures. For this purpose, we restrict our analysis to the one-dimensional domain $\Omega =(0,l)$ with $l\gt 0$ for simplicity and linearise the system (1.2) at the equilibrium $(u_s,v_s,0)$ to get

(4.7) \begin{equation} \left \{\begin{array}{l@{\quad}l} u_{t}=d_{u} u_{x x}-u_s w_x+J_1 u+J_2 v, & x \in (0,l),\ t\gt 0, \\[5pt] v_{t}=d_{v} v_{x x}+J_3 u+J_4 v, & x \in (0,l),\ t\gt 0, \\[5pt] w_{t}=d_{w} w_{x x}+\gamma v_x, & x \in (0,l),\ t\gt 0, \\[5pt] u_{x}=v_{x}=0,\ w=0, & x=0, l,\ t\gt 0. \end{array}\right. \end{equation}

The linearised system (4.7) has solutions in the form of (cf. [40, Appendix]):

(4.8) \begin{equation} \left \{\begin{array}{l} u(x, t)=\sum _{k\geq 0} U_{k} e^{\lambda t} \cos k x, \\[5pt] v(x, t)=\sum _{k\geq 0} V_{k} e^{\lambda t} \cos k x,\\[5pt] w(x, t)=\sum _{k\geq 0} W_{k} e^{\lambda t} \sin k x, \end{array}\right. \end{equation}

where the constants $U_k$ , $V_k$ and $W_k$ are determined by Fourier expansions of the initial conditions, $\lambda$ (depending on $k$ ) is the temporal growth rate and $k=\frac{N\pi }l$ is the wavenumber with the mode $N=0,1,2,\cdots$ . Substituting (4.8) into (4.7), we obtain

(4.9) \begin{equation} \left \{\begin{array}{l@{\quad}c@{\quad}l} \sum _{k\geq 0}\left (\lambda U_{ k}+d_uk^2U_{k}+k u_s W_k-J_1 U_k-J_2 V_k\right )e^{\lambda t} \cos k x=0,\\[8pt] \sum _{k\geq 0}\left (\lambda V_{ k}+d_vk^2V_{k}-J_3U_k-J_4V_k\right )e^{\lambda t} \cos k x=0,\\[8pt] \sum _{k\geq 0}\left (\lambda W_k+d_wk^2W_k+\gamma k V_k\right )e^{\lambda t} \sin k x=0. \end{array} \right. \end{equation}

Recall a basic fact

\begin{equation*} \int ^l_0\cos \!\left (\frac {N\pi }lx\right )\cos \!\left (\frac {M\pi }lx\right )dx= \left \{\begin {array}{l@{\quad}c@{\quad}l} l,\quad M=N=0,\\[5pt] \dfrac l2,\quad M=N\gt 0,\\[5pt] 0,\quad M\neq N. \end {array} \right. \end{equation*}

Multiplying the first two equations of (4.9) by $\cos \frac{N\pi }lx$ , integrating the results with respect to $x$ over $(0,l)$ and using $e^{\lambda t}\neq 0$ , one has

(4.10) \begin{equation} \left \{\begin{array}{l} \lambda U_{ k}+d_uk^2U_{k}+k u_s W_k-J_1 U_k-J_2 V_k=0,\\[5pt] \lambda V_{ k}+d_vk^2V_{k}-J_3U_k-J_4V_k=0. \end{array}\right. \end{equation}

Similarly, using the third equation of (4.9) and

\begin{equation*} \int ^l_0\sin \!\left (\dfrac {N\pi }lx\right )\sin \!\left (\frac {M\pi }lx\right )dx= \left \{\begin {array}{l@{\quad}c@{\quad}l} \frac l2,\quad M=N\gt 0,\\[5pt] 0,\quad M=N=0\text { or }M\neq N, \end {array} \right. \end{equation*}

one has

(4.11) \begin{equation} \lambda W_{ k}+d_wk^2W_{k}+\gamma{k} V_{ k}=0,\quad k=\frac{N\pi }l\neq 0,\ N=1,2,3,\cdots. \end{equation}

When $k=0$ , the third equation of (4.9) holds true for any $\lambda \in \mathbb C$ , and it follows from (4.10) that

\begin{equation*} \lambda \!\left ( \begin{array}{c}{U_k} \\[5pt]{V_k} \end{array} \right ) =\left ( \begin{array}{c@{\quad}c@{\quad}c}{J_1} & J_2 \\[5pt]{J_3} &J_4 \end{array} \right ) \left ( \begin{array}{c}{U_k} \\[5pt]{V_k} \end{array} \right ), \end{equation*}

which implies that $\lambda$ is the eigenvalue of matrix $\boldsymbol{J}$ satisfying (4.5). Therefore, when $k=0$ , by the discussion of the roots of (4.5) we know that for the prey-only steady state $(0,1,0)$ and the coexistence steady state $(u_*,v_*,0)$ , where $(u_*,v_*)$ is given by (4.6), the two roots of (4.5) have negative real parts for $F(v)=v$ with $\alpha \lt \mu$ and $F(v)=\frac v{1+v}$ with $\alpha \lt 2\mu$ , and hence it follows that

(4.12) \begin{equation} \left \{\begin{array}{l} \left .\lim _{t\rightarrow +\infty } U_{k} e^{\lambda t} \cos k x\right |_{k=0}=\lim _{t\rightarrow +\infty } U_{0} e^{\lambda t}=0, \\[5pt] \left .\lim _{t\rightarrow +\infty } V_{k} e^{\lambda t} \cos k x\right |_{k=0}=\lim _{t\rightarrow +\infty } V_{0} e^{\lambda t}=0,\\[5pt] \left .W_{k} e^{\lambda t} \sin k x\right |_{k=0}=0\quad \textrm{for all }t\gt 0. \end{array}\right. \end{equation}

This means $N=0$ is a stable mode.

When $k=\frac{N\pi }l\neq 0$ , $N=1,2,\cdots$ , it follows from (4.10) and (4.11) that

\begin{equation*} \lambda \left ( \begin{array}{c}{U_k} \\[5pt]{V_k} \\[5pt]{W_k} \\[5pt] \end{array} \right ) =\mathscr{A} \!\left ( \begin{array}{c}{U_k} \\[5pt]{V_k} \\[5pt]{W_k} \\[5pt] \end{array} \right ) \quad \textrm{with}\quad \mathscr{A}= \left ( \begin{array}{c@{\quad}c@{\quad}c}{-d_uk^2+J_1} & J_2& -k u_s \\[5pt]{J_3} &-d_vk^2+J_4 & 0\\[5pt]{0}&-k\gamma &-d_wk^2 \\[5pt] \end{array} \right ), \end{equation*}

which implies that  $\lambda$ is the eigenvalue of matrix $\mathscr{A}$ . Calculating the eigenvalue of matrix $\mathscr{A}$ , we obtain the eigenvalue $\lambda \big(k^2\big)$ depending on the wavenumber $k$ as the root of

(4.13) \begin{equation} P(\lambda )\,:\!=\,\lambda ^3+a_2\big(k^2\big)\lambda ^2+a_1\big(k^2\big) \lambda +a_0\big(k^2\big)=0, \end{equation}

where

(4.14) \begin{equation} \left \{\begin{array}{l} a_2=\left (d_u+d_v+d_w\right )k^2 -(J_1+J_4),\\[5pt] a_1=\left (d_ud_v+d_ud_w+d_wd_v\right )k^4-\left [(d_v+d_w)J_1+(d_u+d_w)J_4\right ]k^2+J_1J_4-J_2J_3,\\[5pt] a_0= d_ud_vd_w k^6-d_w(d_vJ_1+d_uJ_4)k^4+\left [d_w\!\left (J_1J_4-J_2J_3\right )-\gamma u_s J_3\right ] k^2. \end{array}\right. \end{equation}

It follows from the Routh–Hurwitz criterion (cf. [35, Appendix B]) for third-order polynomials that $(u_s, v_s, 0)$ is linearly stable (i.e., all roots of (4.13) have negative real parts) if and only if

(4.15) \begin{equation} a_0\gt 0,\ a_2\gt 0\quad \textrm{and}\quad \Gamma \,:\!=\,a_1a_2-a_0\gt 0. \end{equation}

Using the explicit expressions of $a_0$ , $a_1$ and $a_2$ in (4.14), we have

(4.16) \begin{eqnarray} \Gamma =\left [(d_u+d_v) k^{2}-\left (J_1+J_4\right )\right ]\left \{\left [(d_u+d_w) k^{2}-J_1\right ]\left [(d_v+d_w) k^{2}-J_4\right ]-J_2 J_3\right \}+J_3 u_s \gamma. \end{eqnarray}

We know that (4.13) has either one real root and a pair of complex conjugate roots or three real roots. Denote the zeros of $P(\lambda )$ by $\lambda _1$ , $\lambda _2$ and $\lambda _3$ . Then for each $k\neq 0$ , there are three cases for the zeros of $P(\lambda )$ :

1. Case 1: $Re(\lambda _i)\lt 0$ , $i=1,2,3$ . This case is equivalent to (4.15).

2. Case 2: $Re(\lambda _i)\leq 0$ , $i=1,2,3,$ and at least one zero of $P(\lambda )$ has zero real part.

3. Case 3: At least one zero of $P(\lambda )$ has a positive real part.

Remark 4.1. The equilibrium $(u_s,v_s,0)$ is linearly stable (or called locally asymptotically stable) in Case 1, marginally stable (or called neutrally stable, which is neither linearly stable nor linearly unstable, see [21] for instance) in Case 2, and linearly unstable in Case 3. When $\Gamma =0$ , the zeros of $P(\lambda )$ are $\lambda _1=-a_2$ and $\lambda _{2,3}=\pm i\sqrt{a_1}$ . When $\Gamma \lt 0$ and $a_i\gt 0$ for $i=0,1,2$ , then at least one zero of $P(\lambda )$ has a positive real part, and consequently the equilibrium $(u_s,v_s,0)$ is linearly unstable.

Remark 4.2. One can see from (4.14) and (4.15) that

\begin{equation*} \lim _{k\rightarrow \infty }\frac \Gamma{k^6}=(d_u+d_v)(d_v+d_w)(d_u+d_w)\gt 0. \end{equation*}

Therefore, (4.15) always holds for sufficient large $k$ when $a_0\gt 0$ and $a_2\gt 0$ . In the case of $a_i\gt 0$ for $i=0,1,2$ , one can see that if $J_3\lt 0$ and $u_s\gt 0$ , then the equilibrium $(u_s,v_s,0)$ is linearly unstable for large $\gamma \gt 0$ . Indeed, for any fixed $k\neq 0$ , let

(4.17) \begin{align} \gamma _*\big(k^2\big)=&-\frac 1{J_3 u_s k^2}\left [(d_u+d_v) k^{2}-\left (J_1+J_4\right )\right ]\times \nonumber \\[5pt] &\left \{\left [(d_u+d_w) k^{2}-J_1\right ]\left [(d_v+d_w) k^{2}-J_4\right ]-J_2 J_3\right \}. \end{align}

If $\gamma _*\big(k^2\big)\gt 0$ , it follows from (4.16) that $\Gamma \gt 0$ for $\gamma \lt \gamma _*\big(k^2\big)$ , $\Gamma =0$ for $\gamma =\gamma _*\big(k^2\big)$ and $\Gamma \lt 0$ for $\gamma \gt \gamma _*\big(k^2\big)$ . If $\gamma _*\big(k^2\big)\leq 0$ , then $\Gamma \lt 0$ for all $\gamma \gt 0$ . In view of Remark 4.1, $(u_s,v_s,0)$ is linearly unstable for $\gamma \gt \max\! \big\{\gamma _*\big(k^2\big),0\big\}$ , which gives the possibility of patterns bifurcating from $(u_s,v_s,0)$ .

We start by considering the equilibrium $(0,1,0)$ under the condition $\alpha F(1)\leq \mu$ . In this situation, we have

\begin{equation*} \mathscr{A}=\left ( \begin{array}{c@{\quad}c@{\quad}c} -d_u k^2+\alpha F(1)-\mu & 0 & 0 \\[5pt] -\alpha F(1) & -d_v k^2-1 & 0 \\[5pt] 0 & -\gamma k & -d_w k^2 \\[5pt] \end{array} \right ). \end{equation*}

One can see that all eigenvalues of $\mathscr{A}$ are negative for $k=\frac{N\pi }l$ ( $N=1,2,\cdots$ ), which together with (4.12) shows that the equilibrium $(0,1,0)$ is linearly stable. Thus, the pattern can only arise possibly from the homogeneous coexistence steady state $(u_*,v_*,0)$ .

Under the condition $\alpha F(1)\gt \mu$ , we next consider the stability of the equilibrium $(u_*,v_*,0)$ where $(u_*,v_*)$ is given by (4.6). For both Lotka–Volterra ( $F(v)=v$ ) and Holling type II $\big(F(v)=\frac v{1+v}\big)$ trophic functions, we have

(4.18) ${J_3} = - \mu \lt 0,\,{a_1} \gt 0,\,{a_2} \gt 0\,{\rm{and}}\,{a_0} \ge 0( \unicode{x201C}=\unicode{x201D} \,{\rm{holds}}\,{\rm{if}}\,{\rm{and}}\,{\rm{only}}\,{\rm{if}}\,k = 0).$

We shall give details of deriving (4.18) later in Remark 4.3. For the Lotka–Volterra trophic function $F(v)=v$ , we can deduce from (4.17) and Remark 4.3 that

(4.19) \begin{eqnarray} \gamma _*\big(k^2\big)= b_3k^4+b_2 k^2+b_1+\frac{b_0}{k^2}, \end{eqnarray}

where $b_i$ ( $i=0,1,2,3$ ) given by:

\begin{align*} \begin{cases} b_3= \displaystyle \frac{\alpha ^2 ({d_u}+{d_v}) ({d_u}+{d_w}) ({d_v}+{d_w})}{\mu (\alpha -\mu )},\\[16pt] b_2= \displaystyle \frac{\alpha ({d_u}+{d_w}) ({d_u}+2{d_v}+{d_w})}{\alpha -\mu },\\[10pt] b_1= \displaystyle \frac{{d_u} \!\left (\alpha ^2-\alpha \mu +\mu \right )+\alpha{d_v} (\alpha -\mu )+{d_w} \mu }{\alpha -\mu },\\[10pt] b_0=\mu, \end{cases} \end{align*}

are all positive constants since $\alpha \gt \mu$ . Denoting

(4.20) \begin{equation} \gamma _1=\inf _{k=\frac{N\pi }{l}\neq 0}\gamma _*\big(k^2\big)=\inf _{k=\frac{N\pi }{l}\neq 0}\left (b_3k^4+b_2 k^2+b_1+\frac{b_0}{k^2}\right ), \end{equation}

in view of (4.12), Remarks 4.1 and 4.2, we see that the equilibrium $\left ({\frac{\alpha -\mu }{\alpha ^2}},{\frac \mu \alpha },0\right )$ of the system (1.2) is

(4.21) \begin{equation} \left \{\begin{array}{l@{\quad}c@{\quad}l} \textrm{linearly stable},&&\gamma \lt \gamma _1,\\[5pt] \textrm{marginally stable},&&\gamma =\gamma _1,\\[5pt] \textrm{linearly unstable},&&\gamma \gt \gamma _1. \end{array}\right. \end{equation}

Similarly, for the Holling type II functional response function $F(v)=\frac{v}{1+v}$ , when $k\neq 0$ , we see that $a_0$ , $a_1$ and $a_2$ given by (4.14) are all positive thanks to $\alpha \gt 2\mu$ and it follows from (4.17) that

\begin{equation*} \gamma _*\big(k^2\big)= d_3k^4+d_2 k^2+d_1+\frac{d_0}{k^2}, \end{equation*}

where $d_i$ ( $i=0,1,2,3$ ) given by:

\begin{align*} \begin{cases} d_3= \displaystyle \frac{(\alpha -\mu )^2 ({d_u}+{d_v}) ({d_u}+{d_w}) ({d_v}+{d_w})}{\mu (\alpha -2 \mu )},\\[12pt] d_2= \displaystyle \frac{2 \mu (\alpha -\mu ) ({d_u}+{d_w}) ({d_u}+2{d_v}+{d_w})}{\alpha (\alpha -2 \mu )},\\[12pt] d_1= \displaystyle \frac{{d_u} \!\left (\alpha ^4-4 \alpha ^3 \mu +5 \alpha ^2 \mu ^2-2 \alpha \mu ^3+4 \mu ^3\right )+\alpha{d_v} (\alpha -\mu )^2 (\alpha -2 \mu )+4{d_w} \mu ^3}{\alpha ^2 (\alpha -2 \mu )},\\[12pt] d_0\displaystyle = \frac{2 \mu ^2 (\alpha -\mu )}{\alpha ^2}, \end{cases} \end{align*}

are all positives constants since $\alpha \gt 2\mu$ . Denoting

(4.22) \begin{equation} \gamma _2=\inf _{k=\frac{N\pi }{l}\neq 0}\gamma _*\big(k^2\big)=\inf _{k=\frac{N\pi }{l}\neq 0}\left (d_3k^4+d_2 k^2+d_1+\frac{d_0}{k^2}\right ), \end{equation}

it follows from (4.12), Remarks 4.1 and 4.2 that the equilibrium $\left (\frac{\alpha -2\mu }{(\alpha -\mu )^2},\frac \mu{\alpha -\mu },0\right )$ of the system (1.2) is

(4.23) \begin{equation} \left \{\begin{array}{l@{\quad}c@{\quad}l} \textrm{linearly stable},&&\gamma \lt \gamma _2,\\[5pt] \textrm{marginally stable},&&\gamma =\gamma _2,\\[5pt] \textrm{linearly unstable},&&\gamma \gt \gamma _2. \end{array}\right. \end{equation}

Remark 4.3. Below we present some details of obtaining (4.18). Recall $f(v)=v(1-v)$ and $K=1$ , and note that the homogeneous coexistence steady state $(u_*,v_*)$ given by (4.6) exists if and only if $\alpha F(1)\gt \mu$ . We next discuss into two cases: $F(v)=v$ (Lotka–Volterra) and $F(v)=\frac v{1+v}$ (Holling type II).

• • $F(v)=v$ . In this case, the condition $\alpha F(1)\gt \mu$ is equivalent to

  (4.24) \begin{align} \alpha \gt \mu. \end{align}
  Then, we have from (4.4) and (4.6) that
  \begin{equation*} \boldsymbol{J}\big |_{(u_s,v_s)=(u_*,v_*)}=\left ( \begin{array}{c@{\quad}c} J_1 &{J_2} \\[5pt]{J_3} &{J_4} \\[5pt] \end{array} \right )\Bigg |_{(u_s,v_s)=\left (\frac{\alpha -\mu }{\alpha ^{2}}, \frac{\mu }{\alpha }\right )}=\left ( \begin{array}{c@{\quad}c} 0 & 1-\frac{\mu }{\alpha } \\[5pt] -\mu & -\frac{\mu }{\alpha } \\[5pt] \end{array} \right ), \end{equation*}
  which along with (4.14) and (4.24) implies that
  \begin{align*} \begin{cases} a_2=\left (d_u+d_v+d_w\right )k^2 +\dfrac{\mu }{\alpha }\gt 0,\\[8pt] a_1=\left (d_ud_v+d_ud_w+d_wd_v\right )k^4+\dfrac{\mu (d_u+d_w)}\alpha k^2+\mu \left (1-\dfrac \mu \alpha \right )\gt 0,\\[10pt] a_0= d_ud_vd_w k^6+\dfrac{\mu d_ud_w}\alpha k^4+\dfrac{\mu (\alpha -\mu ) (\gamma +\alpha d_w)}{\alpha ^2} k^2\geq 0, \end{cases} \end{align*}
  where $a_0=0$ if and only if $k=0$ .

• • $F(v)=\frac v{1+v}$ . In this case, the condition $\alpha F(1)\gt \mu$ is equivalent to

  (4.25) \begin{align} \alpha \gt 2\mu. \end{align}
  Then, we have from (4.4) and (4.6) that
  \begin{equation*} \boldsymbol{J}\big |_{(u_s,v_s)=(u_*,v_*)}=\left ( \begin{array}{c@{\quad}c} J_1 &{J_2} \\[5pt]{J_3} &{J_4} \\[5pt] \end{array} \right )\Bigg |_{(u_s,v_s)=\left (\frac{\alpha -2\mu }{(\alpha -\mu )^2},\frac \mu{\alpha -\mu }\right )}=\left ( \begin{array}{c@{\quad}c} 0 & 1-\frac{2 \mu }{\alpha } \\[5pt] -\mu & -\frac{2 \mu ^2}{\alpha (\alpha - \mu ) } \\[5pt] \end{array} \right ), \end{equation*}
  which alongside (4.14) and (4.25) implies that
  \begin{align*} \begin{cases} a_2=\left (d_u+d_v+d_w\right )k^2 +\dfrac{2 \mu ^2}{\alpha (\alpha - \mu ) }\gt 0,\\[10pt] a_1=\left (d_ud_v+d_ud_w+d_wd_v\right )k^4+\dfrac{2 \mu ^2(d_u+d_w)}{\alpha (\alpha - \mu ) } k^2+\mu \!\left (1-\dfrac{2\mu }\alpha \right )\gt 0,\\[10pt] a_0= d_ud_vd_w k^6+\dfrac{2 \mu ^2 d_u d_w}{\alpha (\alpha - \mu ) } k^4+\dfrac{\mu (\alpha -2 \mu ) \left (\alpha \gamma +d_w (\alpha -\mu )^2\right )}{\alpha (\alpha -\mu )^2} k^2\geq 0, \end{cases} \end{align*}
  where $a_0=0$ if and only if $k=0$ .

4.3. Spatiotemporal patterns

The model (1.1) describes that predators respond to the heterogeneously distributed prey density by accelerating towards the localities where prey are abundant, which results in predator aggregation. When reaching a local maximum prey concentration, predators decelerate because the prey gradient reverses. Predator aggregations lead to local prey extinctions, while patches with low predator densities may provide partial refuges where prey densities grow. Then predators move actively to the newly formed prey clusters to aggregate to start a new cycle. As such one may expect time-periodic patterns with spatial heterogeneity promoting the persistence of predator–prey interactions (cf. [14]). In this subsection, we shall use numerical simulations to illustrate that the model (1.1) can generate the phenomena described above. In the previous subsection, by taking the preytaxis coefficient $\gamma$ as a bifurcation parameter, we show that spatial patterns may arise from the vicinity of the homogeneous coexistence steady state $\left (u_*,v_*,0\right )$ as $\gamma$ is greater than a critical value $\gamma _1$ (for Holling type 1 trophic function) or $\gamma _2$ (for Holling type II trophic function). Below we shall numerically illustrate the typical patterns generated by systems (4.1) and (4.2). Unless otherwise specified, in this section we take the value of parameters in all simulations as follows:

(4.26) \begin{equation} \alpha =0.3,\quad \mu =0.14\quad \textrm{and}\quad d_u=d_v=d_w=0.1. \end{equation}

We remark the specific values of the above parameters will not qualitatively affect the numerical pattern formations to be shown. Solving (1.6), we get

(4.27) \begin{equation} \left (u_*, v_*\right )=\left \{\begin{array}{l@{\quad}l} \left (\dfrac{16}9,\dfrac 7{15}\right ), & \textrm{ if } F(v)=v,\\[10pt] \left (\dfrac{25}{32},\dfrac 78\right ), & \textrm{ if } F(v)=\dfrac v{1+v}. \end{array}\right. \end{equation}

When $F(v)=v$ , it follows from (4.19) that

\begin{eqnarray*} \left .\gamma _*\big(k^2\big)\right |_{k=\frac{N\pi }{10}}&=&\left .\left (\frac{9 k^4}{280}+\frac{3 k^2}{20}+\frac{47}{200}+\frac{7}{50 k^2}\right )\right |_{k=\frac{N\pi }{10}}\\[5pt] &=& \frac{9 \pi ^4 N^4}{2,800,000}+\frac{3 \pi ^2 N^2}{2000}+\frac{47}{200}+\frac{14}{\pi ^2 N^2},\quad N=1,2,3,\cdots. \end{eqnarray*}

In order to specify $\gamma _1$ according to (4.20), we define a continuous function:

\begin{eqnarray*} H(s)\,:\!=\,\frac{9 \pi ^4 s^4}{2,800,000}+\frac{3 \pi ^2 s^2}{2000}+\frac{47}{200}+\frac{14}{\pi ^2 s^2},\quad s\gt 0. \end{eqnarray*}

One can check that $H^{\prime}(s)=\frac{9 \pi ^4 x^3}{700,000}+\frac{3 \pi ^2 x}{1000}-\frac{28}{\pi ^2 x^3}$ has only one zero $s_0$ in $(0,\infty )$ , and $H^{\prime}(s)$ is negative in $(0,s_0)$ and positive in $(s_0,\infty )$ . Hence, we know that $H(s)$ achieves its minimum at $s_0$ . Direct calculations show that $s_0\in (2,3)$ and

\begin{align*}H(2) & =\frac {47}{200}+\frac {7}{2 \pi ^2}+\frac {3 \pi ^2}{500}+\frac {9 \pi ^4}{175,000}\approx 0.653851, \\H(3)& =\frac {47}{200}+\frac {14}{9 \pi ^2}+\frac {27 \pi ^2}{2000}+\frac {729 \pi ^4}{2,800,000}\approx 0.5512.\end{align*}

Since $\gamma _*\big(k^2\big)$ is defined for discrete value $k\in (0,\infty )$ , $\gamma _*\big(k^2\big)$ attains its minimum at mode $N=3$ , namely at $k=\frac{3\pi }{10}$ . Hence, by (4.20) we have

(4.28) \begin{eqnarray} \gamma _1=\left .\gamma _*\big(k^2\big)\right |_{k=\frac{3\pi }{10}}=\frac{729 \pi ^4}{2,800,000}+\frac{27 \pi ^2}{2000}+\frac{47}{200}+\frac{14}{9 \pi ^2}\approx 0.5512. \end{eqnarray}

Similarly, when $F(v)=\frac v{1+v}$ , it follows from (4.19) that $\gamma _*\big(k^2\big)=\frac{64 k^4}{875}+\frac{224 k^2}{375}+\frac{6956}{5625}+\frac{392}{5625 k^2}.$ One can check that $\gamma _*\big(k^2\big)$ attains its minimum at mode $N=2$ . Hence, by (4.22) we have

(4.29) \begin{eqnarray} \gamma _2=\left .\gamma _*\big(k^2\big)\right |_{k=\frac{2\pi }{10}}=\frac{64 \pi ^4}{546,875}+\frac{224 \pi ^2}{9375}+\frac{6956}{5625}+\frac{392}{225 \pi ^2}\approx 1.6604. \end{eqnarray}

The initial data $(u_0,v_0,w_0)$ are taken as a small random perturbation of the coexistence steady state $\left (u_*,v_*,0\right )$ with 1% deviation :

(4.30) \begin{equation} (u_0,v_0,w_0)=(u_*+0.01\cdot R,v_*+0.01\cdot R,0.01\cdot R), \end{equation}

where $R$ is a random variable taking values in $(-1,1)$ generated by the Matlab and $(u_*, v_*)$ is given in (4.27).

We show numerical simulations in the following by differentiating the Holling type I and Holling type II trophic functions.

Figure 1. Numerical simulation of spatiotemporal patterns generated by (4.1) with $\gamma =1$ in the interval $[0,10]$ , where the initial value $(u_{0}, v_{0}, w_0)$ is given by (4.30) with $(u_*,v_*)=\left(\frac{16}9,\frac 7{15}\right)$ and other parameter values are chosen as in (4.26).

4.3.1. Spatiotemporal patterns for $F(v)=v$

With Holling type I (i.e., Lotka–Volterra) trophic function $F(v)=v$ , the system (1.2) becomes (4.1). Recall that $\gamma _1\approx 0.5512$ in view of (4.28). It follows from (4.21) that the equilibrium $\left (\frac{16}9,\frac 7{15},0\right )$ is linearly stable if $\gamma \lt \gamma _1$ , marginally stable if $\gamma =\gamma _1$ and linearly unstable if $\gamma \gt \gamma _1$ . Therefore, we expect patterns will arise in the supercritical case (i.e., $\gamma \gt \gamma _1$ ). In the critical case $\gamma =\gamma _1$ , in principle the equilibrium $(u_*,v_*, 0)$ may become unstable, stable or remain marginally stable upon a random small perturbation, and our analysis can not confirm which one will occur. Hence, it is also of interest to numerically explore the critical case $\gamma =\gamma _1$ to see whether patterns can develop from the marginally stable steady states upon a small perturbation.

Figure 2. Numerical simulation of spatiotemporal patterns generated by (4.1) with $\gamma =4$ in the interval $[0,10]$ , where the initial value $(u_{0}, v_{0}, w_0)$ is given by (4.30) with $(u_*,v_*)=\left(\frac{16}9,\frac 7{15}\right)$ and other parameter values are chosen as in (4.26).

The numerical spatial-temporal patterns generated by the model (4.1) for the supercritical case are plotted in Figure 1(a) where we observe the spatially inhomogeneous temporal periodic patterns arising from the vicinity of equilibrium $(u_*,v_*,0)$ , and an inhomogeneous limit cycle is eventually attained as shown in Figure 1(c). The spatial distributions of predators and prey at a fixed time plotted in Figure 1(b) show that predators and prey are nearly segregated in space to achieve a heterogeneous coexistence, where predators’ aggregation depresses local prey density while patches with low predators densities provide refuges for prey to promote the persistence at a desirable low level of prey densities. If prey are pests and predators are natural enemies of pests, the model (4.1) is relevant to a successful biological control without local outbreaks as discussed in [40]. If the value of prey-tactic coefficient $\gamma$ is increased to $\gamma =4$ , we still observe the spatially inhomogeneous temporal periodic patterns in the long run as shown in Figure 2(a), although the amplitude and periodicity of periodic patterns are different from the smaller prey-tactic coefficient. Moreover, predators with larger prey-tactic coefficients will be more concentrated in space (compare Figure 2(b) with Figure 1(b)). However, if the prey-tactic coefficient $\gamma$ is sufficiently large, then periodic patterns will be destroyed and chaotic dynamics will present (see Figure 3). However, the prey is not locally eradicated in space (see Figure 3(b)), which is because if predators are too active, they may easily move away from the locations with low prey density and hence save the prey from local eradication (cf. [34]). Next, we numerically explore the critical case $\gamma =\gamma _1$ and corresponding numerical patterns are shown in Figure 4 where spatially inhomogeneous temporal periodic patterns are observed as illustrated in Figure 4(a), a stable limit cycle is achieved as plotted in Figure 4(c) and predator and prey coexist with heterogeneous distributions nearly segregated in space as shown in Figure 4(b).

Figure 3. Numerical simulation of spatiotemporal patterns generated by (4.1) with $\gamma =200$ in the interval $[0,10]$ , where the initial value $(u_{0}, v_{0}, w_0)$ is given by (4.30) with $(u_*,v_*)=\left(\frac{16}9,\frac 7{15}\right)$ and other parameter values are chosen as in (4.26).

Figure 4. Numerical simulation of spatiotemporal patterns generated by (4.1) with $\gamma =\gamma _1=0.5512$ in the interval $[0,10]$ , where the initial value $(u_{0}, v_{0}, w_0)$ is given by (4.30) with $(u_*,v_*)=\left(\frac{16}9,\frac 7{15}\right)$ and other parameter values are chosen as in (4.26).

From the above numerical simulations, we find that spatially inhomogeneous but temporal periodic patterns will typically arise from the preytaxis system (1.2) with Lotka–Volterra trophic function, where predators and prey are nearly segregated in space and persist in time at a low level of prey densities when the preytaxis strength is moderate. If preytaxis is strong, then chaotic dynamics will develop with more pronounced local aggregation of predators but prey can be persistent. This implies preytaxis is a factor driving aggregation of predators but is unable to eradicate the prey. It was well known that if preytaxis in the predator-prey model is conventional (i.e., $\mathbf{w} \sim \nabla v$ ), no spatial patterns will arise (cf. [18]). However, our linear analysis in the previous subsection shows that if preytaxis is modelled with prey-induced acceleration by assuming that predator acceleration instead of velocity is proportional to prey density gradient, spatial patterns may arise as numerically shown in Figures 1, 2 and3. This is a significant difference between conventional and preytaxis with prey-induced acceleration. Hence, our results indicate that preytaxis with prey-induced acceleration is capable of promoting the spatiotemporal heterogeneity in the predator–prey systems and therefore is more appropriate to interpret the experimental observations as in [19, 20, 51].

Figure 5. Numerical simulation of spatiotemporal patterns generated by (4.2) with $\gamma =4$ in the interval $[0,10]$ , where the initial value $(u_{0}, v_{0}, w_0)$ is given by (4.30) with $(u_*,v_*)=\left(\frac{25}{32},\frac 7{8}\right)$ and other parameter values are chosen as in (4.26).

Figure 6. Numerical simulation of spatiotemporal patterns generated by (4.2) with $\gamma =\gamma _2\approx 1.6604$ in the interval $[0,10]$ , where the initial value $(u_{0}, v_{0}, w_0)$ is given by (4.30) with $(u_*,v_*)=\left(\frac{25}{32},\frac 7{8}\right)$ and other parameter values are chosen as in (4.26).

4.3.2. Spatiotemporal patterns $F(v)=\frac v{1+v}$

For the system (1.2) with Holling type II trophic function $F(v)=\frac v{1+v}$ , namely the model (4.2), recalling from (4.29) we have $\gamma _2\approx 1.6604$ . Then (4.23) asserts that the equilibrium $\left (\frac{25}{32},\frac 78,0\right )$ is linearly stable if $\gamma \lt \gamma _2$ , marginally stable if $\gamma =\gamma _2$ and linearly unstable if $\gamma \gt \gamma _2$ . Hence, pattern formations are expected when the prey-tactic coefficient $\gamma$ exceeds or is equal to the critical value $\gamma _2$ . The pattern formations of (1.2) with Holling type II trophic function are qualitatively similar to Holling type I (i.e., Lotka–Volterra) trophic function, and hence we only plot the spatiotemporal patterns formed for the supercritical case in Figure 5 and for the critical case in Figure 6. In [40], it was numerically illustrated that the chaotic dynamics will develop when the prey-tactic coefficient $\gamma$ is large. Our numerical simulations shown in Figure 7 not only demonstrate that the predator density $u$ has chaotic dynamics but also show that the prey density $v$ becomes periodic asymptotically in time. This is different from what is observed for the case of Holling type I trophic function shown in Figure 3 where the dynamics of both predator and prey density are chaotic. This is also the major difference between Holling type I and Holling type II trophic functions observed in numerical simulations.

Figure 7. Numerical simulation of spatiotemporal patterns generated by (4.2) with $\gamma =300$ in the interval $[0,10]$ , where the initial value $(u_{0}, v_{0}, w_0)$ is given by (4.30) with $(u_*,v_*)=\left(\frac{25}{32},\frac 7{8}\right)$ and other parameter values are chosen as in (4.26).