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. [Reference Okubo, Chiang and Ebbesmeyer36]) and in swarms of flying insects (cf. [Reference Parrish and Turchin37]) are directed towards the centroid of such dynamically stable formations, the moving flea-beetles modify their acceleration in response to food patch quality (cf. [Reference Kareiva19]), individual fish in schools adjust their variation of velocity according to the difference between ambient and preferred temperatures (cf. [Reference Flierl, Grünbaum, Levins and Olson11]), and the average velocity is directed by the increasing individual reproduction rate in species clustering (cf. [Reference Grindrod12]). 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 [Reference Kareiva19] and Winder et al. [Reference Winder, Alexander, Holland, Woolley and Perry51], the following reaction–diffusion–preytaxis model was proposed in [Reference Arditi, Tyutyunov, Morgulis, Govorukhin and Senina6, Reference Sapoukhina, Tyutyunov and Arditi40]:
\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 [Reference Murdoch, Briggs and Nisbet33]. 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 [Reference Flierl, Grünbaum, Levins and Olson11] (see more modelling details in [Reference Arditi, Tyutyunov, Morgulis, Govorukhin and Senina6, Reference Sapoukhina, Tyutyunov and Arditi40]). The preytaxis system (1.1) with prey-induced acceleration can generate spatiotemporal patterns (cf. [Reference Arditi, Tyutyunov, Morgulis, Govorukhin and Senina6, Reference Sapoukhina, Tyutyunov and Arditi40]) qualitatively consistent with the observed spatiotemporal heterogeneity in experiments in [Reference Kareiva and Odell20, Reference Okubo, Chiang and Ebbesmeyer36, Reference Parrish and Turchin37]). 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. [Reference Czárán9, Reference Grünbaum13, Reference Kareiva and Odell20, Reference Turchin43] for instance), fail to achieve this explanatory power (cf. [Reference Jin and Wang18] 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 [Reference Lee, Hillen and Lewis25], pattern formation [Reference Chakraborty, Singh, Lucy and Ridland8, Reference He and Zheng15, Reference Lee, Hillen and Lewis26, Reference Li, Wang and Shao28, Reference Wang, Song and Shao48, Reference Wang, Wang and Zhang49], global solvability and stability [Reference Ahn and Yoon1, Reference Ainseba, Bendahmane and Noussair3, Reference Cai, Cao and Wang7, Reference Luo30, Reference Tao41, Reference Wang and Wang44, Reference Wang and Wang45, Reference Wu, Shi and Wu53, Reference Xiang54], 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. [Reference Arditi, Tyutyunov, Morgulis, Govorukhin and Senina6] 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 [Reference Sapoukhina, Tyutyunov and Arditi40] 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. [Reference Chakraborty, Singh, Lucy and Ridland8] 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. [Reference Thakur, Gupta and Upadhyay42] 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 [Reference Rai, Upadhyay and Thakur39] 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:
\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
\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:
(H1)
$F(v)\in C^{2}([0, \infty ))$
,
$F(0)=0$
and
$F(v)\gt 0$
in
$(0, \infty ).$
(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:
\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 )$
:
\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:
\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:
(H3)
$F^{\prime}(v)\gt 0$
in
$[0, \infty )$
.(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
\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
\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$
:
\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
\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. [Reference Leoni27, 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. [Reference Kikuchi and Liu22]).
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:
\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. [Reference Ahn and Yoon1, Reference Ahn and Yoon2, Reference Mishra and Wrzosek31, Reference Wang and Wang46, Reference Xing, Zheng and Pan55, Reference Zuo and Song56] 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 [Reference Amann4, Reference Amann5] (cf. also [Reference Wang and Hillen50, 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
\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:
\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:
\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 [Reference Amann4, Theorem 7.3 and Corollary 9.3] and the blow-up criteria (2.2) follows from [Reference Amann5, 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
\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 [Reference Lieberman29, 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 [Reference Quittner and Souplet38, Proposition 48.4*, 48.5 and 48.7*].
Lemma 2.3. ([Reference Quittner and Souplet38, 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.
(i) If
$1 \leq q\lt p \leq \infty$
, then
(2.5)holds for all
\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}
$z \in L^{q}(\Omega )$
.
(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}
(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 [Reference Winkler52, 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.
(i) If
$1 \leq q\leq p \leq \infty$
, then for any
$z \in L^{q}(\Omega )$
, it holds that
(2.8)and
\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}
(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}
(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}
(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 [Reference Ladyženskaja, Solonnikov and Uralćeva23] (see also [Reference Mora32, Theorem 2.2]), we can find constants
$C_1,C_2,C_3\gt 0$
such that
\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 [Reference DiBenedetto10, 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. [Reference DiBenedetto10, 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
\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
\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. [Reference Cai, Cao and Wang7, formula (1.12)] and [Reference Mora32, 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 [Reference Winkler52, Lemma 1.3 (iii)].
(iii) With (2.9) at hand, (2.11) can be proved by an argument similar to the proof of [Reference Winkler52, Lemma 1.3 (iv)]. Although [Reference Winkler52, 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. [Reference Lankeit24, 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
\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
\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
\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
\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. [Reference Winkler52, 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
\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
\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. [Reference Winkler52, Lemma 1.3]) shows that
\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*}
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 [Reference Jin and Wang17, 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
\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 [Reference Wang and Wang46, Theorem 2.1] for instance). For the reader’s convenience, we shall sketch the proof below. We rewrite (1.2) as:
\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*}
\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 [Reference Lieberman29, 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 [Reference Ladyženskaja, Solonnikov and Uralćeva23] 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. ([Reference Wang47, 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
\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
\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
\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
\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
\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
\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
\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 [Reference Wang and Wang46, 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
\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. [Reference Wang and Wang46, (3.12)], see also [Reference Hu16, Remark 6.1]) shows that
\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
\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
\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
\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
\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
\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
\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
\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:
\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:
(i) If
$\alpha \leq \mu$
, then
where the convergence is exponential if
\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*}
$\alpha \lt \mu$
.
(ii) If
$\alpha \gt \mu$
and
then
\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*}
where
\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*}
$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:
\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:
(i) If
$\alpha \leq 2\mu$
, then
where the convergence is exponential if
\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*}
$\alpha \lt 2\mu$
.
(i) If
$\alpha \gt 2\mu$
and
then
\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*}
where
\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*}
$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
\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:
\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:
\begin{equation} \rho ^2-(J_1+J_4)\rho +J_1J_4-J_2J_3=0. \end{equation}
By the Routh–Hurwitz criterion (cf. [Reference Murray35, 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
\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
\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. [Reference Sapoukhina, Tyutyunov and Arditi40, Appendix]):
\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
\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
\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
\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
\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
\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
\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. [Reference Murray35, 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
\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
\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 )$
:
Case 1:
$Re(\lambda _i)\lt 0$
,
$i=1,2,3$
. This case is equivalent to (4.15).Case 2:
$Re(\lambda _i)\leq 0$
,
$i=1,2,3,$
and at least one zero of
$P(\lambda )$
has zero real part.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 [Reference Kerlin and Upadhyaya21] 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
\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
${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
\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
\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
\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
\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
\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)Then, we have from (4.4) and (4.6) that
\begin{align} \alpha \gt \mu. \end{align}
which along with (4.14) and (4.24) implies 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*}
where
\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*}
$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)Then, we have from (4.4) and (4.6) that
\begin{align} \alpha \gt 2\mu. \end{align}
which alongside (4.14) and (4.25) implies 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*}
where
\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*}
$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. [Reference Hassell and Anderson14]). 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:
\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
\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
\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
\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 :
\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.
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.
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 [Reference Sapoukhina, Tyutyunov and Arditi40]. 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. [Reference Murdoch, Chesson and Chesson34]). 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).
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. [Reference Jin and Wang18]). 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 [Reference Kareiva19, Reference Kareiva and Odell20, Reference Winder, Alexander, Holland, Woolley and Perry51].
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 [Reference Sapoukhina, Tyutyunov and Arditi40], 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.
Acknowledgements
The authors are very thankful to the anonymous referees for their stimulating questions and valuable comments, which greatly help us improve the precision and exposition of this paper.
Funding
The research of C. Mu is partially supported by the NSF of China (grant nos. 11971082, 12271064), the Fundamental Research Funds for the Central Universities (grant nos. 2019CDJCYJ001. 2020CDJQY-Z001, 2021CDJZYJH-004, 2022CDJJCLK002), the NSF of Chongqing (grant nos. cstc2021jcyj-msxmX1051, cstc2022ycjh-bgzxm0169), Key Laboratory of Nonlinear Analysis and its Applications (Chongqing University), Ministry of Education, and Chongqing Key Laboratory of Analytic Mathematics and Applications. The research of W. Tao is partially supported by PolyU Postdoc Matching Fund Scheme Project ID P0030816/B-Q75G, 1-W15F and 1-YXBT, and the NSF of China (no. 12201082). The research of Z.-A. Wang is partially supported by Hong Kong RGC GRF grant no. PolyU 15307222 and Postdoc Matching Fund Scheme Project ID P0034904 (Primary Work Programme W15F).
Competing interests
The authors declare that they have no competing interests.
Open access
