Global well-posedness and Turing–Hopf bifurcation of prey-taxis systems with hunting cooperation

Abstract

This paper is concerned with a predator–prey system with hunting cooperation and prey-taxis under homogeneous Neumann boundary conditions. We establish the existence of globally bounded solutions in two dimensions. In three or higher dimensions, the global boundedness of solutions is obtained for the small prey-tactic coefficient. By using hunting cooperation and prey species diffusion as bifurcation parameters, we conduct linear stability analysis and find that both hunting cooperation and prey species diffusion can drive the instability to induce Hopf, Turing and Turing–Hopf bifurcations in appropriate parameter regimes. It is also found that prey-taxis is a factor stabilizing the positive constant steady state. We use numerical simulations to illustrate various spatiotemporal patterns arising from the abovementioned bifurcations including spatially homogeneous and inhomogeneous time-periodic patterns, stationary spatial patterns and chaotic fluctuations.

1. Introduction

Prey-taxis, the movement of predators toward regions of higher prey density, plays important roles in biological control such as regulating prey (pest) population to avoid incipient outbreaks of prey or forming large-scale aggregation for survival (cf. [13, 32, 44]). This mechanism was first applied to the predator–prey systems by Kareiva and Odell [26] to interpret the heterogeneous aggregative patterns observed in a field experiment for one predator and one prey involving area-restricted search strategies. The Kareiva–Odell prey-taxis model generally reads as

(1.1) \begin{align} \begin{cases} u_t=d_1 \Delta u-\chi \nabla \cdot \left (u\nabla v\right )+\beta G(v,u) u-\theta u,&x\in \Omega ,\ t\gt 0,\\ v_t=d_2 \Delta v+F(v) v-G(v,u) u,&x\in \Omega ,\ t\gt 0,\\ \frac {\partial u}{\partial \nu }=\frac {\partial v}{\partial \nu }=0,\quad & x\in \partial \Omega , \end{cases} \end{align}

where $u(x,t)$ and $v(x,t)$ denote predator and prey density at position $x\in \Omega$ and time $t$ , respectively, and $\Omega \subset {\mathbb {R}}^n(n\geq 1)$ is a bounded domain with smooth boundary. The function $G(v, u)$ is the functional response that describes the predator’s consumption rate of prey, and $F(v)$ represents the per capita prey growth rate in the absence of predators. $d_1$ and $d_2$ are diffusion rates of the predator and the prey, respectively. The parameters $\beta$ and $\theta$ represent the conversion rate of captured prey into predator and the predator mortality rate, respectively. All parameters are positive unless otherwise stated.

The functional response function $G(u,v)$ is a crucial factor shaping various dynamical behaviors in predator–prey systems [53]. Over the past century, different functional responses have been identified based on various biological applications, including Holling type [16–18], Hassell-Varley type [15], Beddington–DeAngelis type [4, 10], ratio-dependent type [3] and Crowley–Martin type [9], among others [53]. The prey growth function $F(v)$ typically takes the form of a logistic type

(1.2) \begin{align} F(v)=\sigma (1-{v}/{K}) \end{align}

with the prey intrinsic growth rate $\sigma$ and the environmental carrying capacity $K$ for the prey. Apart from the logistic type, $F(v)$ could be bistable (or Allee effect) type [55] or fear effect type [56], and so on. In the sequel, we shall assume $F(v)$ is of logistic type (1.2) unless otherwise stated.

This paper will be concerned with the functional response function employing the hunting cooperation strategy, which has attracted extensive attention recently. Cooperative hunting is a prominent behaviour among large social carnivores, which enhances their ability to capture prey in their natural habitats [34, 58]. This behaviour enables predators to subdue large prey and improve their hunting success. For instance, hyenas and wolves usually hunt alone when pursuing small-size prey such as gazelles and sheep. In contrast, they prefer to hunt in packs when hunting large-size prey such as zebras and deer [27, 33, 46]. This communal hunting phenomenon, leading to increased foraging efficiency, has been observed in higher predator densities of terrestrial carnivores, including lions [45], wolves [46] and African wild dogs [8]. To model this phenomenon of hunting cooperation, several functional responses that increase with respect to the predator density, also known as hunting cooperation functional responses, have been proposed, including

(1.3) \begin{align} G(v,u)= \left \{ \begin{aligned} &(\lambda +\alpha u) v,& \quad &\text {(Type I [51])},\\ &\frac {c e (\lambda +\alpha u) v}{1+c h e (\lambda +\alpha u) v},& \quad &\text {(Type II [5,7])},\\ &\frac {c e (\lambda +\alpha u) v^2}{1+c h e (\lambda +\alpha u) v^2},& \quad &\text {(Type III [53])}, \end{aligned} \right . \end{align}

where $\lambda \gt 0$ represents the predation rate on prey and $\alpha \geq 0$ characterizes the level of predator cooperation during hunting, $h$ represents the handling time per prey item, $e_0$ is the encounter rate per predator per prey unit time, and $c$ is a fraction of a prey item killed per predator per encounter. When $\alpha =0$ , the Type I, II and III hunting cooperation response functions reduce to the well-known Holling I, II and III functional response functions, respectively.

Among other things, this paper is focused on the Type I hunting cooperation:

(1.4) \begin{align} G(v,u)=(\lambda +\alpha u) v \end{align}

which was first proposed in [52] and the corresponding temporal system (namely system (1.1) with (1.4) and $d_1=d_2=\chi =0$ ) was numerically studied in [52]. It was found that hunting cooperation can benefit predator populations by increasing attack rates, and large values of $\alpha$ can induce Hopf bifurcations.

The system (1.1) with $\chi =0$ is referred to as the diffusive predator–prey system. It was shown that (1.1) with (1.4) will not have Turing instability without hunting cooperation (i.e. $\alpha =0$ ) (cf. [6, 60]), while it will do with hunting cooperation (i.e. $\alpha \gt 0$ ) only if the predator spreads slower than the prey (i.e. $d_1\gt d_2$ ) (cf. [6, 47, 60]). Later, the Hopf bifurcation was investigated in [29] by the centre manifold and the normal form theory. If the Neumann boundary condition was replaced by the homogeneous Dirichlet boundary condition $u|_{\partial \Omega }=v|_{\partial \Omega }=0$ and $n\lt 6$ , the stationary solution of the system (1.1) with (1.4) was studied in [40] showing that sufficiently large $\alpha$ leads to the extinction of the predator species. For more related works on the predator–prey models with hunting cooperation, we refer readers to [14, 20, 28, 50, 54] for Allee effects in prey, [35] for fear effects in prey, [12] for group defense in prey, [48] for three-species food chain, [11, 31, 41, 51, 54] for hunting cooperation functional response functions other than Type I, and references therein.

Though a considerable body of literature has studied the dynamics of predator–prey systems with hunting cooperation in the absence of prey-taxis (i.e. $\chi =0)$ as mentioned above, only a few results are available to the predator–prey systems with prey-taxis and hunting cooperation (cf. [38, 42, 43, 62]). When the hunting cooperation functional response function $G(u,v)$ is Type II as given in (1.3), the existence of non-constant positive steady states of prey-taxis system (1.1) was obtained by the bifurcation theory in one dimension in [37] where it was also shown that small prey-tactic sensitivity $\chi \gt 0$ can induce the Turing instability. Recently the work [62] established the existence of globally bounded classical solutions of (1.1) with Type II hunting cooperation in any dimension, and showed that negative prey-tactic sensitivity $\chi \lt 0$ (i.e. repulsive prey-taxis) can induce the Turing instability. When the hunting cooperation functional response function $G(u,v)$ is uniformly bounded for all $u,v\geq 0$ , the global boundedness of solutions was established [42, 43] where the global stability of constant positive steady states and existence/nonexistence of non-constant positive solutions were further obtained if $G(u,v)$ is Type II. We remak the assumption in [42, 43] rules out the Type I hunting cooperation functional response function (1.4) since it is obviously not uniformly bounded for $u,v\geq 0$ .

It is seen from the abovementioned existing results that the prey-taxis system (1.1) with Type I hunting cooperation functional response function $G(u,v)$ has not been studied, which is the most difficult case to study from a mathematical point of view among three types given in (1.3) since $G(u,v)$ is not uniformly bounded for $u,v\geq 0$ if it is of Type I. When $G(u,v)$ is uniformly bounded for all $u,v\geq 0$ like Type II or Type III, the regularity of the Neumann heat semigroup (cf. [36]) can be directly applied to establish global boundedness of solutions, while Type I forfeits this advantage. Apart from this, it is also unclear whether the Type I hunting cooperation may induce the spatiotemporal patterns. This paper will address these questions and we therefore consider the following prey-taxis system (1.1) with Type I hunting cooperation

(1.5) \begin{align} \begin{cases} u_t=d_1 \Delta u-\chi \nabla \cdot \left (u\nabla v\right )+u\left [\beta \left (\lambda +\alpha u\right )v-\theta \right ],&x\in \Omega ,\ t\gt 0,\\ v_t=d_2 \Delta v+\sigma v\left (1-\frac vK\right )-\left (\lambda +\alpha u\right )u v,&x\in \Omega ,\ t\gt 0,\\ \frac {\partial u}{\partial \nu }=\frac {\partial v}{\partial \nu }=0,&x\in \partial \Omega ,\\ u(x, 0)=u_0(x),\ v(x, 0)=v_0(x),&x\in \Omega , \end{cases} \end{align}

where the parameters $d_1,d_2,\beta ,\lambda ,\alpha ,\theta ,\sigma ,K\gt 0$ and $\chi \geq 0$ have the same biological interpretations as mentioned above. To reduce the number of parameters, we introduce the following rescalings

\begin{align*} \tilde d_1=\frac {d_1}\theta , \quad \tilde d_2=\frac {d_2}\theta , \quad \tilde \chi =\frac {\chi }{\beta \lambda }, \quad \tilde \sigma =\frac \sigma \theta , \quad \tilde K=\frac {\beta \lambda K}\theta , \quad \tilde \alpha =\frac {\alpha \theta }{\lambda ^2}, \end{align*}

and

\begin{align*} \tilde u=\frac {\lambda u}\theta , \quad \tilde v=\frac {\beta \lambda v}\theta , \quad \tilde t=\theta t. \end{align*}

Substituting the above rescalings into (1.5) and dropping the tildes for brevity, we obtain the following rescaled system

(1.6) \begin{align} \begin{cases} u_t=d_1 \Delta u-\chi \nabla \cdot \left (u\nabla v\right )+\left (1+\alpha u\right )uv-u,&x\in \Omega ,\ t\gt 0,\\ v_t=d_2 \Delta v+\sigma v\left (1-\frac vK\right )-\left (1+\alpha u\right )uv,&x\in \Omega ,\ t\gt 0,\\ \frac {\partial u}{\partial \nu }=\frac {\partial v}{\partial \nu }=0,&x\in \partial \Omega ,\\ u(x, 0)=u_0(x),\ v(x, 0)=v_0(x),&x\in \Omega , \end{cases} \end{align}

where the parameters $d_1,d_2,\alpha ,K,\sigma$ are positive and $\chi \gt 0$ . Throughout the paper, we shall often use the following notations

(1.7) \begin{align} f(u,v)\,:\!=\, \left (1+\alpha u\right )uv-u \quad \text {and}\quad g(u,v)\,:\!=\, \sigma v\left (1-\frac vK\right )-\left (1+\alpha u\right )uv \end{align}

for brevity.

The main analytical result of this paper is the following.

Theorem 1.1. Let $\Omega \subset \mathbb {R}^{n}$ ( $n\geq 2$ ) be a bounded domain with smooth boundary, and let $(u_{0}, v_{0}) \in [W^{1, q}(\Omega )]^2$ with some $q\gt n$ and $u_{0}(x), v_{0}(x) \geq 0(\!\not \equiv 0)$ . Then there exists $\chi _*\in (0,+\infty ]$ such that for all $\chi \in {(}0,\chi _*)$ , the system (1.6) admits a unique classical solution $(u,v)\in [C^0(\overline {\Omega } \times [0, \infty )) \cap C^{2,1}(\overline {\Omega } \times (0, \infty ))]^2$ satisfying $u,v\gt 0$ for all $t\gt 0$ . Moreover, there exists a constant $C\gt 0$ independent of $t$ such that

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

In particular, $\chi _*=+\infty$ if $n=2$ .

Before ending this section, we outline the main difficulties and proof strategies.

1.1 Sketch of proof strategies

Without prey-taxis (i.e. $\chi =0$ ), problem (1.6) can be regarded as a special one considered in [38]. Then the global boundedness of solutions to (1.6) with $\chi =0$ in any dimension ( $n\geq 1$ ) is a consequence of [38, Theorem 3.1] (see also [19]) based on the $L^p$ -duality method. Without hunting cooperation (i.e. $\alpha =0$ ), the global boundedness of solutions to (1.6) has been established in [23, 61]. Now with $\alpha \gt 0$ , the source term contains a cubic polynomial $\alpha u^2 v$ and many arguments and estimates in [23, 61] for the quadratic polynomial are no longer applicable. Hence, new technical ingredients are needed to deal with the higher-order nonlinearity. Below we briefly describe our strategies used to obtain the global boundedness of solutions for $n=2$ and $n\geq 3$ .

1. (i) $n=2$ . In this case, by fully exploiting the model structural feature, we use cancellation ideas to handle the troubling prey-taxis term while simultaneously controlling the cubic nonlinear source term $\alpha u^2v$ (see the proof of Lemma 3.1), and finally obtain the crucial uniform-in-time bound of $\|u\|_{L^{2}(\Omega )}$ with sophisticated coupling estimates by establishing a Grownwall type inequality for the linear combination of following four terms (see (3.19) and (3.27))

   \begin{align*} \frac d{dt}\int _\Omega u^2,\quad \frac d{dt}\int _\Omega \frac {|\nabla v|^2}{v},\quad \frac d{dt}\int _\Omega u\ln u,\quad \frac d{dt}\int _\Omega uv. \end{align*}
   Then we further perform $L^p$ estimates to obtain the uniform-in-time bound of $\|u\|_{L^{p}(\Omega )}$ for $p\gt 2$ , which yields the global boundedness of solutions by the boundedness criterion in Lemma 2.6 and the extension criterion in Lemma 2.1. We stress that here we obtain the $L^2$ -estimates directly with the sophisticated coupling estimates by fully exploiting the advantage of cubic decay term $-\alpha u^2v$ . This is different from the estimates in the existing literature for the prey-taxis model without hunting cooperation such as [23, 61] where the $L^2$ -estimate of $u$ is established based on the estimate of $\|u\ln u\|_{L^1}$ and $\|\nabla v\|_{L^2}$ which are first obtained separately.

2. (ii) $n\geq 3$ . In this case, we first establish a uniform-in-time bound of the weighted $L^p$ -norm $\int _\Omega u^p\varphi (v)$ ( $p\gt 1$ ) by constructing an appropriate weight function $\varphi (v)$ (see Lemma 4.1 and Lemma 4.2), and then obtain a bound of $\|u\|_{L^{p}(\Omega )}$ under a smallness condition on $\chi \gt 0$ .

The rest of the paper is organized as follows. In Sect. 2, we establish the local existence of solutions to the system (1.6) and derive some preliminary results. We prove Theorem 1.1 for $n=2$ in Sect. 3 and for $n\geq 3$ in Sect. 4. Finally, in Sect. 5, we conduct linear stability analysis to show that both Turing bifurcations and Turing–Hopf bifurcations may arise from the system (1.6), and use numerical simulations to demonstrate that (1.6) may generate various complex spatiotemporal patterns such as spatially homogeneous time-periodic patterns, stationary patterns, spatiotemporal periodic patterns and chaos.

2. Local existence and preliminaries

In this section, we establish the local existence of solutions to the system (1.6) and provide some preliminary results. Hereafter, we use $C$ and $C_i$ ( $i=1,2,3,\cdots$ ) to denote generic positive constants that may vary from line to line in the context. We begin with the local existence of solutions.

Lemma 2.1. Suppose that the assumptions in Theorem 1.1 hold. Then there exists $T_{\max } \in (0, +\infty ]$ such that the system (1.6) has a unique classical solution

\begin{align*} (u, v) \in \left [C^0 (\overline {\Omega } \times (0, T_{\max }))\cap C^{2,1}(\overline {\Omega } \times (0, T_{\max }))\right ]^2 \end{align*}

satisfying

(2.1) \begin{align} u\gt 0\quad \text {and}\quad 0\lt v\leq M_0\,:\!=\, \max \left \{\|v_0\|_{L^{\infty }(\Omega )},K\right \} \quad \text {for all }(x,t)\in \Omega \times (0,T_{\max }). \end{align}

Moreover,

(2.2) \begin{align} \text {either}\quad T_{\max }=+\infty \quad \text {or}\quad \limsup _{t \nearrow T_{\max }}\left (\|u(\cdot ,t)\|_{L^{\infty }(\Omega )}+\|v(\cdot ,t)\|_{L^\infty (\Omega )}\right )=+\infty . \end{align}

Proof. Denote $\omega =(u, v)$ . Then the system (1.6) can be written as

\begin{align*} \begin{cases} \omega _t=\nabla \cdot (\mathbb {A}(\omega ) \nabla \omega )+\Phi (\omega ), & x \in \Omega , t\gt 0, \\ \frac {\partial \omega }{\partial \nu }=0, & x \in \partial \Omega , \\ \omega (x, 0)=\left (u_0(x), v_0(x)\right ), & x \in \Omega , \end{cases} \end{align*}

where $\mathbb {A}(\omega )=\left [\begin{array}{cc} d_1 & -\chi u \\ 0 & d_2 \end{array}\right ]$ and $ \Phi (\omega )= \left ( \begin{array}{c} f(u,v) \\ g(u,v) \end{array} \right )$ with $f(u,v)$ and $g(u,v)$ given by (1.7). The upper triangular matrix $\mathbb {A}(\omega )$ is positively definite. Hence, the system (1.6) is normally parabolic. Then the local existence and uniqueness of classical solutions follow from Amann’s theorem [1, Theorem 7.3 and Corollary 9.3] and the blow-up criteria (2.2) follows from [2, Theorem 15.5]. Moreover, $u,v\gt 0$ for $(x, t) \in \Omega \times \left (0, T_{\max }\right )$ can be established by the strong maximum principle along with $u_0,v_0\geq 0(\not \equiv 0)$ . These arguments are standard and we refer readers to [57, Lemma2.6] and [24, Lemma2.1] for details. Finally, an application of [23, Lemma2.2] proves $v\leq M_0$ for $(x, t) \in \Omega \times \left (0, T_{\max }\right )$ .

The following result can be easily obtained.

Lemma 2.2. Suppose that the assumptions in Theorem 1.1 hold and $(u, v)$ is the solution of the system (1.6) . Then

(2.3) \begin{align} \|u(\cdot ,t)\|_{L^{1}(\Omega )}\leq M_1\,:\!=\, \max \left \{\|u_0+v_0\|_{L^{1}(\Omega )},\frac {K(1+\sigma )^2}{4\sigma }|\Omega |\right \}\quad \text {for all }t\in (0,T_{\max }). \end{align}

Proof. Adding the first two equations of (1.6), integrating the resulting equation by parts and using Young’s inequality $(1+\sigma )v\leq \frac \sigma K v^2+\frac {K(1+\sigma )^2}{4\sigma }$ , we have

\begin{align*} \frac {d}{dt}\int _\Omega (u+v) = \sigma {\int _\Omega }v\left (1-\frac vK\right )-{\int _\Omega }u\leq -{\int _\Omega }(u+v)+\frac {K(1+\sigma )^2}{4\sigma }{|\Omega |} \quad \text {for all }t\in (0,T_{\max }), \end{align*}

which alongside $u,v\geq 0$ implies

\begin{align*} \int _\Omega (u+v)\leq \max \left \{\|u_0+v_0\|_{L^{1}(\Omega )},\frac {K(1+\sigma )^2}{4\sigma }|\Omega |\right \}. \end{align*}

The proof is completed.

We shall recall some inequalities that will be used later.

Lemma 2.3 ([22, Lemma 2.3]). Let $\alpha , \beta , T\gt 0$ , $\tau \in (0, T)$ , and suppose that ${\phi }:[0, T) \rightarrow$ $[0, \infty )$ is absolutely continuous satisfying

\begin{align*} {\phi }^{\prime }(t)+{\phi }^{1+\sigma }(t) \leq {\varphi }(t) {\phi }(t)+h(t), \quad t\gt 0, \end{align*}

where $\sigma \gt 0$ is a constant, ${\varphi }(t), h(t) \geq 0$ with ${\varphi }(t), h(t) \in L_{l o c}^1([0, T))$ and

\begin{align*} \sup _{t \in [\tau , T)} \int _{t-\tau }^t {\varphi }(s) d s \leq \alpha , \quad \sup _{t \in [\tau , T)} \int _{t-\tau }^t h(s) d s \leq \beta . \end{align*}

Then for any $t\gt t_0$ , we have

(2.4) \begin{align} {\phi }(t) \leq {\phi }\left (t_0\right ) e^{\int _{t_0}^t {\varphi }(s) d s}+\int _{t_0}^t h(\tau ) e^{\int _\tau ^t {\varphi }(s) d s} d \tau \end{align}

and

(2.5) \begin{align} \sup _{t \in (0, T)} {\phi }(t) \leq \sigma \left (\frac {2 A}{1+\sigma }\right )^{\frac {1+\sigma }{\sigma }}+2 B, \sup _{t \in [\tau , T)} \int _{t-\tau }^t {\phi }^{1+\sigma }(s) d s \leq (1+\alpha ) \sup _{t \in (0, T)}\{{\phi }(t)\}+\beta , \end{align}

where

(2.6) \begin{align} A=\tau ^{-\frac {1}{1+\sigma }}(1+\alpha )^{\frac {1}{1+\sigma }} e^{2 \alpha }, \quad B=\tau ^{-\frac {1}{1+\sigma }} \beta ^{\frac {1}{1+\sigma }} e^{2 \alpha }+2 \beta e^{2 \alpha }+{\phi }(0) e^\alpha . \end{align}

Lemma 2.4 ([59, Lemma 3.3]). Suppose that $H \in C^1((0, \infty ),{\mathbb {R}}_+)$ and $\Theta (s)\,:\!=\, \int _1^s \frac {d \sigma }{H(\sigma )}$ for $s\gt 0$ . Then for all $\varphi \in C^2(\overline {\Omega })$ fulfilling $\frac {\partial \varphi }{\partial v}=0$ on $\partial \Omega$ , it holds that

\begin{align*} \int _{\Omega } \frac {H^{\prime }(\varphi )}{H^3(\varphi )}|\nabla \varphi |^4 \leq (2+\sqrt {n})^2 \int _{\Omega } \frac {H(\varphi )}{H^{\prime }(\varphi )}\left |D^2 \Theta (\varphi )\right |^2. \end{align*}

Lemma 2.5 ([30, Lemma 4.2]). Assume that $\Omega$ is a bounded domain, and let $w \in C^2(\overline {\Omega })$ satisfy $\frac {\partial w}{\partial \nu }=0$ on $\partial \Omega$ . Then we have

\begin{align*} \frac {\partial |\nabla w|^2}{\partial \nu } \leq 2 \kappa |\nabla w|^2, \end{align*}

where $\kappa =\kappa (\Omega )$ is an upper bound of the curvatures of $\partial \Omega$ .

Now we establish a boundedness criterion.

Lemma 2.6. Suppose that the assumptions in Theorem 1.1 hold and $(u, v)$ is the solution of the system (1.6) . If there exists $p \gt \frac {3n}{2}$ ( $n\geq 2$ ) such that

(2.7) \begin{align} \sup _{t \in \left (0, T_{\max }\right )}\|u(\cdot , t)\|_{L^{p }(\Omega )} \leq K_1 \end{align}

for a positive constant $K_1$ , there exists a positive constant $K_2$ independent of $t$ such that

\begin{align*} \|u(\cdot , t)\|_{L^{\infty }(\Omega )} \leq K_2\quad \text {for all }t\in (0,T_{\max }). \end{align*}

Proof. Assume that (2.7) holds for some $p \gt \frac {3n}{2}$ ( $n\geq 2$ ). Then there exists a constant $C_1\gt 0$ such that

\begin{align*} \sup _{t \in \left (0, T_{\max }\right )}\|u^2(\cdot , t)\|_{L^{\frac {p }2}(\Omega )} \leq C_1. \end{align*}

Since $\frac {p }2\gt \frac {3n}4$ , one can use a standard argument based on the smooth property of the Neumann heat semigroup (cf. [36]) to find two constants $r_1\gt 3n$ and $C_2\gt 0$ such that

(2.8) \begin{align} \sup _{t \in \left (0, T_{\max }\right )}\|\nabla v(\cdot , t)\|_{L^{r_1}(\Omega )} \leq C_2. \end{align}

Given $t\in (0,T_{\max })$ , we let $t_0\,:\!=\, (t-1)_+$ . By Duhamel’s principle, $u$ can be represented as

(2.9) \begin{align} u(\cdot ,t)&= e^{d_1(t-{t_0})\Delta }u(\cdot ,t_0) -\chi \int ^t_{t_0} e^{d_1(t-s)\Delta }\nabla \cdot \left [u(\cdot ,s)\nabla v(\cdot ,s)\right ]ds+\int ^t_{t_0} e^{d_1(t-s)\Delta } f(u(\cdot ,s),v(\cdot ,s))ds\nonumber \\ &=\!:\,I_1+I_2+I_3\quad \text {for all }t\in (0,T_{\max }). \end{align}

By (2.3), the maximum principle and the smooth property of the Neumann heat semigroup again, for all $t\in (0, T_{\max })$ , we can find two positive constants $C_3$ and $C_4$ such that

(2.10) \begin{align} \|I_1\|_{L^{\infty }(\Omega )}= \begin{cases} \|e^{d_1t\Delta }u(\cdot ,0)\|_{L^{\infty }(\Omega )}\leq \|u_0\|_{L^{\infty }(\Omega )},\quad &\text {if }t\leq 1,\\ \|e^{d_1\Delta }u(\cdot ,t-1)\|_{L^{\infty }(\Omega )}\leq C_3 \|u(\cdot ,t-1)\|_{L^{1}(\Omega )}\leq C_4,&\text {if }t\gt 1. \end{cases} \end{align}

Let $r_2\,:\!=\, \frac { {p }r_1}{ {p }+r_1}$ . Then $r_2\gt n$ due to $\frac 1{r_2}=\frac 1{p }+\frac 1{r_1}\lt \frac 1n$ . By (2.8) and Hölder’s inequality, for all $t\in (0,T_{\max })$ , we have

\begin{align*} \sup _{t \in \left (0, T_{\max }\right )}\|u(\cdot , t)\nabla v(\cdot , t)\|_{L^{r_2}(\Omega )} \leq \sup _{t \in \left (0, T_{\max }\right )}\|u(\cdot , t)\|_{L^{p }(\Omega )}\|\nabla v(\cdot , t)\|_{L^{r_1}(\Omega )} \leq K_1C_2, \end{align*}

which alongside the smooth property of the Neumann heat semigroup gives two positive constants $C_5$ and $C_6$ such that

(2.11) \begin{align} \|I_2\|_{L^{\infty }(\Omega )} \leq C_5\int ^t_{t_0} \left (1+(t-s)^{-\frac 12-\frac {n}{2r_2}}\right ) e^{-\lambda _1d_1(t-s)}\|u(\cdot ,s)\nabla v(\cdot ,s)\|_{L^{r_2}(\Omega )}ds \leq C_6 \end{align}

for all $t\in (0,T_{\max })$ , where we have used the fact $\frac 12+\frac {n}{2r_2}\lt 1$ , and $\lambda _1\gt 0$ denotes the first nonzero eigenvalue of $-\Delta$ in $\Omega$ under the homogeneous Neumann boundary condition.

It follows from (1.7), (2.1), (2.7), Young’s inequality and Hölder’s inequality that there exist two constants $C_7,C_8\gt 0$ such that

(2.12) \begin{align} \int _\Omega |f(u(\cdot ,t),v(\cdot ,t))|^{\frac {3n}4}\leq C_7\int _\Omega \left (u^{\frac {3n}2}+1\right )\leq C_8\quad \text {for all }t\in (0,T_{\max }). \end{align}

Let $\overline f(t)\,:\!=\, \frac 1{|\Omega |}\int _\Omega f(u(\cdot ,t),v(\cdot ,t))$ for all $t\in (0,T_{\max })$ . Notice that $\frac {3n}4\gt 1$ for $n\geq 2$ . Then (2.12) and Hölder’s inequality imply

(2.13) \begin{align} |\overline f(t)|\leq \frac {1}{|\Omega |}\int _\Omega (|f(u(\cdot ,t),v(\cdot ,t))|^{\frac {3n}4}+1)\leq \frac {C_8}{|\Omega |}+1 \quad \text {for all }t\in (0,T_{\max }). \end{align}

Using (2.12), (2.13), $t-t_0\leq 1$ , and the $L^p$ - $L^q$ estimates of the Neumann heat semigroup, we can find two positive constants $C_9$ and $C_{10}$ satisfying

(2.14) \begin{align} \|I_3\|_{L^{\infty }(\Omega )}\leq &\, \int ^t_{t_0} \|e^{d_1(t-s)\Delta }f(\cdot ,s)\|_{L^{\infty }(\Omega )} ds\nonumber \\ &\leq \int _{t_0}^{t}\| e^{d_1(t-s)\Delta }(f (\cdot ,s)-\overline f (s)) \|_{L^{\infty }(\Omega )} d s+ \int _{t_0}^{t}\| e^{d_1(t-s)\Delta }\overline f (s) \|_{L^{\infty }(\Omega )} d s\nonumber \\ &\leq C_9 \int _{t_0}^{t} (1+(t-s)^{-\frac 23})e^{-\lambda _1d_1(t-s)}\|f (\cdot ,s)-\overline f (s)\|_{L^{\frac {3n}4}(\Omega )} d s+\int _{t_0}^{t} \left (\frac {C_8}{|\Omega |}+1\right ) ds\nonumber \\ &\leq C_{10}\quad \text {for all }t\in (0,T_{\max }). \end{align}

The combination of (2.9), (2.10), (2.11) and (2.14) yields that

\begin{align*} \sup _{t \in \left (0, T_{\max }\right )}\|u(\cdot , t)\|_{L^{\infty }(\Omega )} \leq C_{11}\quad \text {for all }t\in (0,T_{\max }) \end{align*}

with a constant $C_{11}\gt 0$ . The proof is completed.

3. Proof of Theorem1.1 for $n=2$

This section is devoted to proving Theorem1.1 for $n=2$ . To this end, we construct an appropriate functional in the form of a linear combination of the four terms

\begin{align*} \frac d{dt}\int _\Omega \frac {|\nabla v|^2}{v},\quad \frac d{dt}\int _\Omega u\ln u,\quad \frac d{dt}\int _\Omega u^2,\quad \frac d{dt}\int _\Omega uv \end{align*}

to derive some a priori estimates including the time-independent boundedness of $\|u\|_{L^2(\Omega )}$ . This idea was first used in [59] for a chemotaxis-fluid model, and then developed for predator–prey models with prey-taxis [23, 25].

3.1 $L^2$ -estimate via sophisticated coupling estimates

Lemma 3.1 Suppose that the assumptions in Theorem 1.1 hold with $n=2$ and $(u, v)$ is the solution of the system (1.6) . Then for $\tau \in \left (0,T_{\max }\right )$ , we have

(3.1) \begin{align} \int _\Omega \frac {|\nabla v|^2}{v} +\int _\Omega u^2\leq C \quad \text {for all }t\in (0,T_{\max }) \end{align}

and

(3.2) \begin{align} \int _{t}^{t+\tau }\int _\Omega \left ( \frac {\left |D^2 v\right |^2}{v}+\frac {|\nabla v|^4}{v^3} + \frac {|\nabla u|^2}{u}+ |\nabla u|^2+u^2|\nabla v|^2 \right )(\cdot ,s)ds \leq C(\tau +1) \end{align}

for all $t\in (0,T_{\max }-\tau )$ , where $D^2 v$ denotes the Hessian matrix of $v$ .

Proof. We split the proof into several steps.

Step 1. An inequality for $\frac d{dt}\int _\Omega \frac {|\nabla v|^2}{v}$ . For all $t\in (0,T_{\max })$ , it holds that

(3.3) \begin{align} \frac 12\frac d{dt}\int _\Omega \frac {|\nabla v|^2}{v} =&\, -\frac 12 \int _\Omega \frac {|\nabla v|^2}{v^2}v_t + \int _\Omega \frac {\nabla v\cdot \nabla v_t}{v}\nonumber \\[3pt] =&\, \frac 12 \int _\Omega \frac {|\nabla v|^2}{v^2}v_t - \int _\Omega \frac {\Delta v }{v}v_t\nonumber \\[3pt] =&\, \underbrace {d_2\left (\frac {1}2 \int _\Omega \frac {|\nabla v|^2}{v^2}\Delta v -\int _\Omega \frac {|\Delta v|^2}{v}\right )}_{=\!:\,I_4} +\frac 12 \int _\Omega \frac {|\nabla v|^2}{v^2} g(u,v) - \int _\Omega \frac {\Delta v }{v}g(u,v). \end{align}

We have from (1.7) and (2.1) that

(3.4) \begin{align} \frac 12 \int _\Omega \frac {|\nabla v|^2}{v^2} g(u,v) \leq \frac 12 \int _\Omega \frac {|\nabla v|^2}{v^2} \left (\sigma v - \alpha u^2v\right ) \leq \frac \sigma 2 \int _\Omega \frac {|\nabla v|^2}{v} -\frac \alpha {2 M_0} \int _\Omega u^2|\nabla v|^2 \end{align}

for all $t\in (0,T_{\max })$ . For the term $I_4$ , we claim that there exists a positive constant $C_1$ such that

(3.5) \begin{align} I_4\leq &\, -\frac {3A_1}{4} \int _\Omega \left (\frac {|D^2 v|^2}{v}+\frac {|\nabla v|^4}{v^3}\right )+C_1\int _\Omega \frac {|\nabla v|^2}{v} \quad \text {for all }t\in (0,T_{\max }), \end{align}

where $A_1=\frac {d_2}{3(2+\sqrt {2})^2+2}$ . Then the combination of (3.3)-(3.5) indicates that

(3.6) \begin{align} &\frac 12\frac d{dt}\int _\Omega \frac {|\nabla v|^2}{v} +\frac 12 \int _\Omega \frac {|\nabla v|^2}{v} +\frac {3A_1}{4}\int _\Omega \left (\frac {|D^2 v|^2}{v}+\frac {|\nabla v|^4}{v^3}\right )\nonumber \\[3pt] \leq &\, -\frac \alpha {2 M_0} \int _\Omega u^2|\nabla v|^2 +\left (\frac {1+\sigma }2+C_1\right ) \int _\Omega \frac {|\nabla v|^2}{v} - \int _\Omega \frac {\Delta v }{v}g(u,v) \quad \text {for all }t\in (0,T_{\max }). \end{align}

Using (2.1) and Young’s inequality, we obtain

(3.7) \begin{align} \left (\frac {1+\sigma }2+C_1\right ) \int _\Omega \frac {|\nabla v|^2}{v} \leq \frac {A_1}{4}\int _\Omega \frac {|\nabla v|^4}{v^3}+\frac {\left (\frac {1+\sigma }2+C_1\right )^2}{A_1}\int _\Omega v \leq \frac {A_1}{4}\int _\Omega \frac {|\nabla v|^4}{v^3}+C_2 \end{align}

for all $t\in (0,T_{\max })$ , where $C_2\,:\!=\, \frac {M_0|\Omega |}{A_1}\left (\frac {1+\sigma }2+C_1\right )^2$ . By (1.7), we have

(3.8) \begin{align} -\int _\Omega \frac {\Delta v }{v}g(u,v) = -\int _\Omega \Delta v\left (\sigma \left (1-\frac vK\right )-\left (1+\alpha u\right )u\right ) = -\frac \sigma K\int _\Omega |\nabla v|^2 -I_5 \leq -I_5, \end{align}

where $I_5\,:\!=\, \int _\Omega \left (1+2\alpha u\right )\nabla u\cdot \nabla v$ . Substituting (3.7) and (3.8) into (3.6), we have

(3.9) \begin{align} \frac 12\frac d{dt}\int _\Omega \frac {|\nabla v|^2}{v} +\frac 12 \int _\Omega \frac {|\nabla v|^2}{v} +\frac {A_1}{2}\int _\Omega \left (\frac {|D^2 v|^2}{v}+\frac {|\nabla v|^4}{v^3}\right ) \leq C_2-\frac \alpha {2 M_0} \int _\Omega u^2|\nabla v|^2 -I_5 \end{align}

for all $t\in (0,T_{\max })$ .

We next prove (3.5) by similar arguments to that of [23, 3.2 and Lemma 3.3] (see also the proof of [25, 3.2 and Lemma 3.3]). For completeness, we sketch the proof here. Applying [23, Equs. (3.19) and (3.27)] with $F(v)=v$ therein, we have

(3.10) \begin{align} I_4 =&\, \frac {d_2}2\int _{\partial \Omega }\frac 1v\cdot \frac {\partial |\nabla v|^2}{\partial \nu }-d_2\int _\Omega v|D^2 \ln v|^2\quad \text {for all }t\in (0,T_{\max }) \end{align}

and

(3.11) \begin{align} \int _\Omega \frac {|D^2 v|^2}{v} \leq 2\left ((2+\sqrt {2})^2+1\right )\int _\Omega v|D^2 \ln v|^2\quad \text {for all }t\in (0,T_{\max }). \end{align}

Using Lemma2.4 with $H(v)=v$ yields

(3.12) \begin{align} \int _\Omega \frac {|\nabla v|^4}{v^3} \leq (2+\sqrt {2})^2\int _\Omega v|D^2 \ln v|^2\quad \text {for all }t\in (0,T_{\max }). \end{align}

To proceed, we recall the trace inequality [39, Remark 52.9]: for any $\varepsilon \gt 0$ , there exists a constant $C(\varepsilon )\gt 0$ such that

(3.13) \begin{align} \|w\|_{L^2(\partial \Omega )} \leq \varepsilon \|\nabla w\|_{L^2(\Omega )}+C(\varepsilon )\|w\|_{L^2(\Omega )}\quad \text {for }w\in W^{1,2}(\Omega ). \end{align}

Denote $\varphi _0(v)=\frac {|\nabla v|}{\sqrt v}$ . Then for all $t\in (0,T_{\max })$ , it follows from Young’s inequality that

\begin{align*} |\nabla \varphi _0(v)|^2 =\left (\frac {D^2v\cdot \nabla v}{\sqrt v|\nabla v|}-\frac {|\nabla v|\nabla v}{2v^\frac 32}\right )^2 \leq \frac {|D^2 v|^2}{v} +\frac {|D^2v|\cdot |\nabla v|^2}{v^2} +\frac {|\nabla v|^4}{4v^3} \leq \frac 54\left (\frac {|D^2 v|^2}{v}+\frac {|\nabla v|^4}{v^3}\right ). \end{align*}

This alongside Lemma 2.5, (3.13) with $\varepsilon =\sqrt {\frac {A_1}{10\kappa d_2}}$ , and the fact that $(a+b)^2\leq 2(a^2+b^2)$ for $a,b\in {\mathbb {R}}$ gives

(3.14) \begin{align} \frac {d_2}{2} \int _{\partial \Omega } \frac {1}{v} \cdot \frac {\partial |\nabla v|^2}{\partial \nu } \leq &\,\frac {d_2}{2} \int _{\partial \Omega } \frac {2\kappa |\nabla v|^2}{v} \nonumber \\[3pt] = &\, \kappa d_2 \|\varphi _0(v)\|^2_{L^2(\partial \Omega )} \nonumber \\[3pt] \leq &\, \kappa d_2 \left [\varepsilon \|\nabla \varphi _0(v)\|_{L^2(\Omega )}+C(\varepsilon )\|\varphi _0(v)\|_{L^2(\Omega )}\right ]^2\nonumber \\[3pt]\leq &\, \kappa d_2\left [2\varepsilon ^2\|\nabla \varphi _0(v)\|^2_{L^2(\Omega )}+2C(\varepsilon )^2\|\varphi _0(v)\|^2_{L^2(\Omega )}\right ]\nonumber \\[3pt]\leq &\, \kappa d_2\left [2\varepsilon ^2\frac 54\int _\Omega \left (\frac {|D^2 v|^2}{v}+\frac {|\nabla v|^4}{v^3}\right )+2C(\varepsilon )^2\int _{\Omega } \frac {|\nabla v|^2}{v}\right ]\nonumber \\[3pt]\leq &\, \frac {A_1}{4} \int _{\Omega }\left (\frac {\left |D^2 v\right |^2}{v}+\frac {|\nabla v|^4}{v^3}\right )+C_1 \int _{\Omega } \frac {|\nabla v|^2}{v}\quad \text {for all }t\in (0,T_{\max }), \end{align}

where $C_1\gt 0$ is a constant. Now the combination of (3.10)-(3.12) and (3.14) indicates

\begin{align*} I_4 \leq &\, \left [\frac {A_1}{4}-\frac {d_2}{3(2+\sqrt {2})^2+2}\right ] \int _\Omega \left (\frac {|D^2 v|^2}{v}+\frac {|\nabla v|^4}{v^3}\right )+C_1\int _\Omega \frac {|\nabla v|^2}{v}\quad \text {for all }t\in (0,T_{\max }), \end{align*}

which gives (3.5) by recalling $A_1=\frac {d_2}{3(2+\sqrt {2})^2+2}$ . Therefore, the claim (3.5) is proved, and hence (3.9) is obtained.

Step 2. Cancellation of the nonlinear term $\int _\Omega (1+2\alpha u)\nabla u\cdot \nabla v$ . We now deal with the nonlinear term $I_5=\int _\Omega \left (1+2\alpha u\right )\nabla u\cdot \nabla v$ appearing in (3.9). For all $t\in (0,T_{\max })$ , using (1.6) and integration by parts, we have

(3.15) \begin{align} \frac d{dt} \int _\Omega u\ln u =&\, -d_1\int _\Omega \frac {|\nabla u|^2}{u}+\chi \int _\Omega \nabla u\cdot \nabla v+\int _\Omega (1+\ln u)f(u,v), \end{align}

(3.16) \begin{align} \frac 12\frac d{dt} \int _\Omega u^2 =&\, -d_1\int _\Omega |\nabla u|^2+\chi \int _\Omega u\nabla u\cdot \nabla v+\int _\Omega u f(u,v), \end{align}

(3.17) \begin{align} \frac d{dt} \int _\Omega uv =&-\left (d_1+d_2\right )\int _\Omega \nabla u\cdot \nabla v +\chi \int _\Omega u|\nabla v|^2+\int _\Omega vf(u,v)+\int _\Omega ug(u,v). \end{align}

Let $A_2\,:\!=\, \frac {2(d_1+d_2)}\chi +\frac 1{2\alpha }$ . Then $\chi A_2+\chi u-2(d_1+d_2)=\frac \chi {2\alpha }(1+2\alpha u)$ , which alongside (3.15)-(3.17) gives

(3.18) \begin{align} & \frac d{dt}\int _\Omega u\left (A_2\ln u+\frac u2+2v\right ) +d_1A_2\int _\Omega \frac {|\nabla u|^2}{u} +d_1\int _\Omega |\nabla u|^2 \nonumber \\ =&\, \frac \chi {2\alpha }I_5 +2\chi \int _\Omega u|\nabla v|^2 +\underbrace {\int _\Omega u\left (f(u,v)+2g(u,v)\right )+A_2\int _\Omega (1+\ln u)f(u,v) +2\int _\Omega vf(u,v)}_{=\!:\,I_6} \end{align}

for all $t\in (0,T_{\max })$ . If we multiply (3.18) by $\frac {2\alpha }\chi$ and add the result to (3.9), then the nonlinear term $I_5$ can be canceled.

Step 3. A Grönwall-type inequality. Define the function

(3.19) \begin{equation} y(t)\,:\!=\, \frac 12 \int _\Omega \frac {|\nabla v|^2}{v} +\frac {2\alpha }\chi \int _\Omega u\left (A_2\ln u+\frac u2+2v\right )\quad \text {for all }t\in (0,T_{\max }). \end{equation}

Then for all $t\in (0,T_{\max })$ , the combination of (3.9) and (3.18) multiplied by $\frac {2\alpha }\chi$ indicates that

(3.20) \begin{align} & y'(t)+y(t) +\frac {A_1}{2}\int _\Omega \left (\frac {|D^2 v|^2}{v}+\frac {|\nabla v|^4}{v^3}\right ) +\frac {2\alpha }\chi \left (d_1A_2\int _\Omega \frac {|\nabla u|^2}{u} +d_1\int _\Omega |\nabla u|^2\right ) \nonumber \\ \leq &\, C_2 +\underbrace {4\alpha \int _\Omega u|\nabla v|^2-\frac \alpha {2 M_0} \int _\Omega u^2|\nabla v|^2}_{=\!:\,I_7} +\underbrace {\frac {2\alpha }\chi \left (\int _\Omega u\left (A_2\ln u+\frac u2+2v\right )+I_6\right )}_{=\!:\,I_8}. \end{align}

We next estimate the term $I_7$ . Indeed, Young’s inequality and (2.1) yield

\begin{align*} 4\alpha \int _\Omega u|\nabla v|^2 \leq &\, \int _\Omega \left (\frac \alpha {4 M_0} u^2+16M_0\alpha \right )|\nabla v|^2 \\ \leq &\, \frac \alpha {4 M_0} \int _\Omega u^2|\nabla v|^2 +16M_0\alpha \int _\Omega \left (\frac {A_1}{64M_0\alpha }\frac {|\nabla v|^4}{v^3}+\frac {16M_0\alpha }{A_1}v^3\right ) \\ \leq &\, \frac \alpha {4 M_0} \int _\Omega u^2|\nabla v|^2 +\frac {A_1}{4}\int _\Omega \frac {|\nabla v|^4}{v^3}+\frac {256M_0^5\alpha ^2}{A_1}|\Omega |\quad \text {for all }t\in (0,T_{\max }), \end{align*}

which implies

(3.21) \begin{align} I_7\leq -\frac \alpha {4 M_0} \int _\Omega u^2|\nabla v|^2+\frac {A_1}{4}\int _\Omega \frac {|\nabla v|^4}{v^3}+\frac {256M_0^5\alpha ^2}{A_1}|\Omega |\quad \text {for all }t\in (0,T_{\max }). \end{align}

It remains to estimate the term $I_8$ . Using (1.7), (2.1), (2.3), Young’s inequality and the fact that $-s\ln s\leq \frac 1e$ and $\ln s\lt s$ for $s\gt 0$ , we have

\begin{align*} \int _\Omega u(f(u,v)+2g(u,v)) =&\,\int _\Omega u\left [-(1+\alpha u)uv-u+2\sigma v\left (1-\frac vK\right )\right ]\\ \leq &\, 2\sigma M_0M_1-\int _\Omega (1+\alpha u)u^2v\quad \text {for all }t\in (0,T_{\max }). \end{align*}

For the term $I_6$ included in $I_8$ , it holds that

(3.22) \begin{align} I_6 \leq &\, 2\sigma M_0M_1-\int _\Omega (1+\alpha u)u^2v \nonumber \\ &+A_2\int _\Omega (1+\ln u)\left (1+\alpha u\right )uv- A_2\int _\Omega (1+\ln u)u+2M_0\int _\Omega \left (1+\alpha u\right )uv\nonumber \\ \leq &\, 2\sigma M_0M_1 -\int _\Omega \left [u-A_2(1+\ln u)-2M_0\right ] (1+\alpha u)uv+\frac {A_2}e|\Omega |\quad \text {for all }t\in (0,T_{\max }). \end{align}

Moreover, for all $t\in (0,T_{\max })$ , using (2.1), (2.3) and the fact that $s^2\geq s\ln s$ for $s\geq 0$ , we obtain

(3.23) \begin{align} \int _\Omega u\left (A_2\ln u+\frac u2+2v\right ) \leq &\, \int _\Omega u\left (A_2u+\frac u2+2M_0\right ) \leq \left (A_2+1\right )\int _\Omega u^2 +2 M_0M_1. \end{align}

Using (2.3), the Gagliardo–Nirenberg inequality and Young’s inequality, for all $t\in (0,T_{\max })$ , we can find two positive constants $C_3$ and $C_4$ such that

(3.24) \begin{align} \frac {2\alpha }\chi \left (A_2+1\right )\int _\Omega u^2 \leq C_3\left (\|\nabla u\|_{L^{2}(\Omega )}^\frac 12\|u\|_{L^{1}(\Omega )}^\frac 12+\|u\|_{L^{1}(\Omega )}\right )^2 \leq \frac {d_1\alpha }{2\chi }\int _\Omega |\nabla u|^2+C_4. \end{align}

The combination of (3.22)-(3.24) shows that

(3.25) \begin{align} I_8 \leq &\, C_5+ \frac {d_1\alpha }{2\chi }\int _\Omega |\nabla u|^2-\frac {2\alpha }\chi \int _\Omega \left [u-A_2(1+\ln u)-2M_0\right ] (1+\alpha u)uv\quad \text {for all }t\in (0,T_{\max }), \end{align}

where $C_5\,:\!=\, C_4+\frac {2\alpha }\chi \left [2 (1+\sigma )M_0M_1+\frac {A_2}e|\Omega |\right ]$ . Define the function

\begin{align*} \varphi (s)\,:\!=\, -(s-A_2(1+\ln s)-2M_0)\qquad \text {for all }\ s\gt 0. \end{align*}

Then $\varphi '(s)=\frac {A_2}s-1$ and hence

\begin{align*} \varphi (s)\leq \varphi (A_2)=2M_0+A_2\ln A_2\leq 2M_0+A_2|\ln A_2|=\!:C_6. \end{align*}

Therefore, by (2.1), (2.3), (3.25) and the same argument as in deriving (3.24), we know that there exists a constant $C_7\gt 0$ such that

(3.26) \begin{align} I_8\leq &\, C_5+\frac {d_1\alpha }{2\chi }\int _\Omega |\nabla u|^2+\frac {2\alpha }\chi C_6\int _\Omega (1+\alpha u)uv\nonumber \\ \leq &\, C_5+\frac {d_1\alpha }{2\chi }\int _\Omega |\nabla u|^2+\frac {2\alpha }\chi C_6 \left (M_0M_1+\alpha M_0\int _\Omega u^2\right )\nonumber \\ \leq &\, C_7+\frac {d_1\alpha }\chi \int _\Omega |\nabla u|^2 \quad \text {for all }t\in (0,T_{\max }). \end{align}

Substituting (3.21) and (3.26) into (3.20) yields

(3.27) \begin{align} &y'(t)+y(t) +\frac {A_1}{4}\int _\Omega \left (\frac {|D^2 v|^2}{v}+\frac {|\nabla v|^4}{v^3}\right ) +\frac {2d_1A_2\alpha }\chi \int _\Omega \frac {|\nabla u|^2}{u} +\frac {d_1\alpha }\chi \int _\Omega |\nabla u|^2 +\frac \alpha {4 M_0} \int _\Omega u^2|\nabla v|^2 \nonumber \\ \leq &\, C_2 +\frac {64M_0^5\alpha ^2}{A_1}|\Omega | +C_7 \quad \text {for all }t\in (0,T_{\max }). \end{align}

Finally, an application of Grönwall’s inequality along with the facts $u,v\geq 0$ and $u\ln u\geq -\frac 1e$ yields (3.1). Furthermore, the integration of (3.27) with respect to $t$ over $[t,t+\tau ]$ gives (3.2).

3.2 The uniform-in-time estimate of $\|u(\cdot ,t)\|_{L^{p}(\Omega )}$ for $p\geq 2$

Lemma 3.2. Suppose that the assumptions in Theorem 1.1 hold with $n=2$ , and $(u, v)$ is the solution of the system (1.6) . Then for any $p\geq 2$ , there exists a constant $C(p)\gt 0$ independent of $t$ such that

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

Proof. Using the first equation of (1.6) and integration by parts, for all $t\in (0,T_{\max })$ , we obtain

(3.29) \begin{align} \frac 1p\frac d{dt}\int _\Omega u^p =&\, \int _\Omega u^{p-1}\left (d_1 \Delta u-\chi \nabla \cdot \left (u\nabla v\right )+f(u,v)\right ) \nonumber \\ =&\, -d_1(p-1)\int _\Omega u^{p-2} |\nabla u|^2 +\chi (p-1)\int _\Omega u^{p-1} \nabla u\cdot \nabla v +\int _\Omega u^{p-1}f(u,v). \end{align}

For all $t\in (0,T_{\max })$ , it follows from Young’s inequality that

(3.30) \begin{align} \chi (p-1)\int _\Omega u^{p-1} \nabla u\cdot \nabla v \leq &\, \frac {d_1(p-1)}4\int _\Omega u^{p-2} |\nabla u|^2 +\frac {\chi ^2(p-1)}{d_1}\int _\Omega u^p|\nabla v|^2\nonumber \\ =&\, \frac {d_1(p-1)}{p^2}\|\nabla u^{\frac p2}\|_{L^{2}(\Omega )}^2 +\frac {\chi ^2(p-1)}{d_1}\int _\Omega u^p|\nabla v|^2. \end{align}

For the last term in the right-hand side of the above inequality, we use Hölder’s inequality $\|u^p|\nabla v|^2\|_{L^{1}(\Omega )}\leq \|u^p\|_{L^{2}(\Omega )}\||\nabla v|^2\|_{L^{2}(\Omega )}$ , the Gagliardo–Nirenberg inequality and Young’s inequality to find a constant $C_1\gt 0$ such that

(3.31) \begin{align} \frac {\chi ^2(p-1)}{d_1}\int _\Omega u^p|\nabla v|^2 \leq &\, \frac {\chi ^2(p-1)}{d_1}\|u^{\frac p2}\|_{L^{4}(\Omega )}^2\|\nabla v\|_{L^{4}(\Omega )}^2\nonumber \\ \leq &\, \frac {\chi ^2(p-1)}{d_1} C_1\left (\|\nabla u^{\frac p2}\|_{L^{2}(\Omega )} \|u^{\frac p2}\|_{L^{2}(\Omega )} +\|u^{\frac p2}\|_{L^{2}(\Omega )}^2\right )\|\nabla v\|_{L^{4}(\Omega )}^2\nonumber \\ \leq &\, \frac {d_1(p-1)}{p^2}\|\nabla u^{\frac p2}\|_{L^{2}(\Omega )}^2 +\frac {p^2}{4 d_1(p-1)}\left (\frac {\chi ^2(p-1)}{d_1} C_1\|u^{\frac p2}\|_{L^{2}(\Omega )} \|\nabla v\|_{L^{4}(\Omega )}^2\right )^2\nonumber \\ &+\frac {\chi ^2(p-1)}{d_1}C_1\|u\|_{L^{p}(\Omega )}^p \|\nabla v\|_{L^{4}(\Omega )}^2\nonumber \\ \leq &\, \frac {d_1(p-1)}{p^2}\|\nabla u^{\frac p2}\|_{L^{2}(\Omega )}^2 +\frac {p^2(p-1)\chi ^4}{4d_1^3}C_1^2\|u\|_{L^{p}(\Omega )}^p \|\nabla v\|_{L^{4}(\Omega )}^4\nonumber \\ &+\frac {\chi ^2(p-1)}{d_1}C_1\|u\|_{L^{p}(\Omega )}^p \left (\|\nabla v\|_{L^{4}(\Omega )}^4+1\right )\nonumber \\ =&\, \frac {d_1(p-1)}{p^2}\|\nabla u^{\frac p2}\|_{L^{2}(\Omega )}^2 +C_2 \|u\|_{L^{p}(\Omega )}^p \|\nabla v\|_{L^{4}(\Omega )}^4+C_3\int _\Omega \left (u^{p+1}+1\right ) \end{align}

for all $t\in (0,T_{\max })$ , where

\begin{align*} C_2\,:\!=\, \frac {p^2(p-1)\chi ^4}{4d_1^3}C_1^2+\frac {\chi ^2(p-1)}{d_1}C_1 \quad \text {and}\quad C_3\,:\!=\, \frac {\chi ^2(p-1)}{d_1}C_1. \end{align*}

For the last term in the right-hand side of (3.29), we have from (1.7), (3.1) and Young’s inequality that

\begin{align*} \int _\Omega u^{p-1}f(u,v) \leq \int _\Omega u^{p-1}\left (1+\alpha u\right )uv \leq M_0\int _\Omega \left (u^p+\alpha u^{p+1}\right ) \leq M_0\left (|\Omega |+(\alpha +1)\int _\Omega u^{p+1}\right ) \end{align*}

for all $t\in (0,T_{\max })$ . Let ${\phi }(t)\,:\!=\, \frac 1p\|u(\cdot ,t)\|_{L^{p}(\Omega )}^p$ for all $t\in (0,T_{\max })$ . Then Hölder’s inequality implies

(3.32) \begin{align} {\phi }(t)^{\frac {p+1}p}=\left (\frac 1p\right )^{\frac {p+1}p}\|u\|_{L^{p}(\Omega )}^{p+1}\leq \left (\frac 1p\right )^{\frac {p+1}p}|\Omega |^{\frac 1p}\|u\|_{L^{p+1}(\Omega )}^{p+1}\quad \text {for all }t\in (0,T_{\max }). \end{align}

For all $t\in (0,T_{\max })$ , the combination of (3.29)-(3.32) implies that

(3.33) \begin{align} {\phi }'(t)+{\phi }(t)^{\frac {p+1}p} +\frac {2d_1(p-1)}{p^2}\|\nabla u^{\frac p2}\|_{L^{2}(\Omega )}^2 \leq C_2 \|u\|_{L^{p}(\Omega )}^p \|\nabla v\|_{L^{4}(\Omega )}^4 +C_4\|u\|_{L^{p+1}(\Omega )}^{p+1}+C_3+M_0|\Omega |, \end{align}

where $C_4\,:\!=\, {C_3+(\alpha +1)M_0}+\left (\frac 1p\right )^{\frac {p+1}p}|\Omega |^{\frac 1p}.$ Again applying the Gagliardo–Nirenberg inequality and Young’s inequality, and using (3.1), one can find two positive constants $C_5$ and $C_6$ such that

(3.34) \begin{align} C_4\|u\|_{L^{p+1}(\Omega )}^{p+1} =C_4\|u^{\frac p2}\|_{L^{\frac {2(p+1)}p}(\Omega )}^{\frac {2(p+1)}p} \leq &\, C_5\left (\|\nabla u^{\frac p2}\|_{L^{2}(\Omega )}^{\frac {2(p+1)}p\cdot \frac {p-1}{p+1}} \|u^{\frac p2}\|_{L^{\frac 4p}(\Omega )}^{\frac {2(p+1)}p\cdot \frac 2{p+1}}+\|u^{\frac p2}\|_{L^{\frac 4p}(\Omega )}^{\frac {2(p+1)}p}\right )\nonumber \\ \leq &\, C_6\left (\|\nabla u^{\frac p2}\|_{L^{2}(\Omega )}^{\frac {2(p-1)}p}+1\right )\nonumber \\ \leq &\, \frac {d_1(p-1)}{p^2}\|\nabla u^{\frac p2}\|_{L^{2}(\Omega )}^2+C_7 \quad \text {for all }t\in (0,T_{\max }), \end{align}

where

\begin{align*} C_7\,:\!=\, \frac 1p\left (\frac p{p-1}\frac {d_1(p-1)}{p^2}\right )^{-(p-1)}C_6^p+C_6=\frac 1p\left (\frac p{d_1}\right )^{p-1}C_6^p+C_6. \end{align*}

Substituting (3.34) into (3.33), for all $t\in (0,T_{\max })$ , we arrive at

(3.35) \begin{align} {\phi }'(t)+{\phi }(t) ^{\frac {p+1}p} +\frac {d_1(p-1)}{p^2}\|\nabla u^{\frac p2}\|_{L^{2}(\Omega )}^2 \leq pC_2 {\phi }(t)\|\nabla v\|_{L^{4}(\Omega )}^4 +C_3+M_0|\Omega |+C_7. \end{align}

Let $\tau =\min \left \{1,\frac {T_{\max }}2\right \}\leq 1$ . Then it follows from (2.1) and (3.2) that there exists a constant $C_8\gt 0$ such that

(3.36) \begin{align} pC_2 \int _{t}^{t+\tau }\|\nabla v\|_{L^{4}(\Omega )}^4 =&\, pC_2 \int _{t}^{t+\tau } {\int _\Omega }\frac {|\nabla v|^4}{v^3}v^3 \leq pC_2 M_0^{3}\int _{t}^{t+\tau } {\int _\Omega }\frac {|\nabla v|^4}{v^3} \leq C_8(\tau +1) \leq 2C_8 \end{align}

for all $t\in (0,T_{\max }-\tau )$ .

We are now in a position to prove (3.28). We shall discuss two cases: $T_{\max }\lt 2$ and $T_{\max }\geq 2$ . If $T_{\max }\lt 2$ , then $\tau \lt 1$ , and (2.4) implies

(3.37) \begin{align} \phi (t) \leq \phi (0) e^{\int _{0}^t pC_2 \|\nabla v(\cdot ,s)\|_{L^{4}(\Omega )}^4 d s}+(C_3+M_0|\Omega |+C_7)\int _{0}^t e^{\int _\tau ^t pC_2 \|\nabla v(\cdot ,s)\|_{L^{4}(\Omega )}^4 d s} d \tau \end{align}

for all $t\in (0,T_{\max })$ . By an argument similar to that used to derive (3.36), we have

\begin{align*} e^{\int _{\tau }^t pC_2 \|\nabla v(\cdot ,s)\|_{L^{4}(\Omega )}^4 d s}\leq e^{\int _{0}^t pC_2 \|\nabla v(\cdot ,s)\|_{L^{4}(\Omega )}^4 d s}\leq e^{\int _{0}^{T_{\max }} pC_2 \|\nabla v(\cdot ,s)\|_{L^{4}(\Omega )}^4 d s}\leq C_9(T_{\max }+1)\leq 3C_9. \end{align*}

This alongside (3.37) gives

\begin{align*} \phi (t) \leq 3h(0) C_9+3(C_3+M_0|\Omega |+C_7)\int _{0}^2 C_9 d \tau \leq 3h(0)C_9 +6(C_3+M_0|\Omega |+C_7)C_9, \end{align*}

which proves (3.28) in the case of $T_{\max }\lt 2$ . If $T_{\max }\geq 2$ , then $\tau =1$ . By (2.5) and (2.6), we get

(3.38) \begin{align} \sup _{t \in (0, T_{\max })} \phi (t) =\sup _{t \in (0, T_{\max })}\frac 1p\|u(\cdot ,t)\|_{L^{p}(\Omega )}^p\leq \frac 1p\left (\frac {2A}{1+\frac 1p}\right )^{p+1}+2B, \end{align}

where $A=(1+2C_8)^{\frac p{1+p}}e^{4C_8}$ and

\begin{align*} B=(C_3+M_0|\Omega |+C_7)^{\frac p{1+p}}e^{4C_8}+2(C_3+M_0|\Omega |+C_7) e^{4C_8}+\frac 1p\|u_0(\cdot )\|_{L^{p}(\Omega )}^pe^{2C_8} \end{align*}

are independent of $T_{\max }$ . Therefore, (3.38) yields (3.28) in the case of $T_{\max }\geq 2$ . The proof is completed.

3.2.1 Proof of Theorem1.1 for $n=2$

Lemma 3.2 gives $\|u\|_{L^{4}(\Omega )}\leq C$ . Then Theorem1.1 with $n=2$ follows from (2.2) and Lemma 2.6 immediately. $\square$

4. Proof of Theorem1.1 for $n\geq 3$

This section is devoted to proving Theorem1.1 for $n\geq 3$ by the weighted $L^p$ -estimates ( $p\gt 1$ ). For $M_0$ given by (2.1) and $p\gt 1$ , we define the following positive constants

(4.1) \begin{align} \begin{cases} \Gamma _1\,:\!=\,\dfrac { d_2 p (p-1) }{p (d_1 ^2+ d_2 ^2)+2 d_1 d_2 },\quad \Gamma _2\,:\!=\,\dfrac {1}{2d_2} \sqrt {p\left (p+\frac {2d_2}{d_1}\right )},\\[10pt] \Gamma _3\,:\!=\,\dfrac {2 d_1 d_2 (p-1)}{p (d_1 ^2+ d_2 ^2)+2 d_1 d_2 },\quad \chi _p\,:\!=\,\dfrac \pi {2\Gamma _2 M_0}. \end{cases} \end{align}

We then construct the following function to be used later.

Lemma 4.1. Let $p\gt 1$ , $\Gamma _1, \Gamma _2, \Gamma _3$ , $\chi _p$ be given by (4.1) and $\chi \in (0,\chi _p)$ . Then for all $s\in [0,M_0]$ , the function

(4.2) \begin{align} \varphi (s)\,:\!=\, e^{\Gamma _1 \chi s}\left [\cos (\chi \Gamma _2 s)\right ]^{-\Gamma _3} \end{align}

satisfies

(4.3) \begin{align} &1\leq \varphi (s)\leq \varphi (M_0)\lt +\infty , \end{align}

(4.4) \begin{align} &0\lt \chi \Gamma _1\leq \frac {\varphi '(s)}{\varphi (s)} =\chi \left [\Gamma _1+\Gamma _2\Gamma _3\tan (\chi \Gamma _2 s)\right ] \leq \chi \left [\Gamma _1+\Gamma _2\Gamma _3\tan (\chi \Gamma _2 M_0)\right ]=\!: k_{\chi ,p}\lt +\infty , \end{align}

(4.5) \begin{align} &\frac {d_2}p\varphi ''(s)-\chi \varphi '(s)-\frac {\left [\left (\chi (p-1)\varphi (s)-(d_1+d_2)\varphi '(s)\right )\right ]^2}{2d_1(p-1)\varphi (s)}=0. \end{align}

Proof. Clearly, (4.1) and $\chi \in (0,\chi _p)$ imply

\begin{align*} 0\lt \chi \Gamma _2 s\leq \chi \Gamma _2 M_0\lt \frac \pi 2\quad \text {for all }s\in [0,M_0]. \end{align*}

Hence (4.1) and (4.2) indicate that $\varphi (s)$ increases in $s\in [0,M_0]$ and (4.3) holds. By a simple calculation, we have

\begin{align*} \frac {\varphi '(s)}{\varphi (s)} =\chi \left [\Gamma _1+\Gamma _2\Gamma _3\tan (\chi \Gamma _2 s)\right ]\quad \text {for all }s\in [0,M_0], \end{align*}

and then (4.4) is obvious. We next prove (4.5). By tedious calculations, we arrive at

(4.6) \begin{align} &2d_1(p-1)\varphi (s)\left (\frac {d_2}p\varphi ''(s)-\chi \varphi '(s)\right )-\left [\left (\chi (p-1)\varphi (s)-(d_1+d_2)\varphi '(s)\right )\right ]^2\nonumber \\ =&-\frac 1p\left (\frac {\varphi (s) }{\cos (\chi \Gamma _2 s)}\right )^2 \left \{ B_1 \cos ^2 (\chi \Gamma _2 s) +\Gamma _2 \Gamma _3 \chi \sin (\chi \Gamma _2 s ) \left [ 2 B_2 \cos (\chi \Gamma _2 s ) +B_3 \Gamma _2 \chi \sin (\chi \Gamma _2 s ) \right ]\right \}, \end{align}

where

\begin{align*} \begin{cases} B_1\,:\!=\,\Gamma _1^2 \chi ^2 \left [p (d_1 ^2+ d_2 ^2)+2 d_1 d_2 \right ] -2 d_1 d_2 \Gamma _2^2 \Gamma _3 (p-1) \chi ^2+(p-1) p \chi \left [(p-1) \chi -2 d_2 \Gamma _1 \chi \right ],\\ B_2\,:\!=\, \Gamma _1 \chi \left [p (d_1 ^2+ d_2 ^2)+2 d_1 d_2 \right ]- d_2 p (p-1) \chi ,\\ B_3\,:\!=\, \Gamma _3\left [p (d_1 ^2+ d_2 ^2)+2 d_1 d_2 \right ]-2 d_1 d_2 (p-1). \end{cases} \end{align*}

Clearly, it follows from (4.1) that $B_2=B_3=0$ and

\begin{align*} B_1= (p-1) \chi ^2 \left [p(p-1) -p d_2\Gamma _1 - 2 d_1d_2 \Gamma _2^2 \Gamma _3 \right ] =\frac { d_1 (p-1)^2 \chi ^2 \left [ d_1 \left (p^2-4 d_2^2 \Gamma _2^2\right )+2 d_2 p\right ]}{ p (d_1 ^2+ d_2 ^2)+2 d_1 d_2 }=0. \end{align*}

This along with (4.3), (4.6) and $p\gt 1$ proves (4.5). The proof is completed.

Now we are in a position to use the above auxiliary function as a weight function to derive the uniform bound of $\|u(\cdot ,t)\|_{L^{p}(\Omega )}$ for all $t\in (0, T_{\max })$ .

Lemma 4.2. Suppose that the assumptions in Theorem 1.1 hold and $(u, v)$ is the solution of the system (1.6) . Let $p\gt 1$ , $\Gamma _1, \Gamma _2, \Gamma _3$ , $\chi _p$ be given by (4.1) , $\chi \in (0,\chi _p)$ and $k_{\chi ,p}$ be given by (4.4) . Then there exists a constant $C(p)\gt 0$ independent of $t$ such that

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

Proof. Let $\varphi (s)$ be defined by (4.2) for $s\in [0,M_0]$ . Then (2.1) and (4.3) show that

(4.7) \begin{align} 1\leq \varphi (v(x,t))\leq \varphi (M_0)\quad \text {for all } (x,t)\in \Omega \times (0,T_{\max }). \end{align}

With integration by parts, one has

(4.8) \begin{align} \frac 1p\frac d{dt}\int _\Omega u^p\varphi (v) =&\, \int _\Omega u^{p-1}\varphi (v)\left (d_1 \Delta u-\chi \nabla \cdot \left (u\nabla v\right )+f(u,v)\right ) +\frac 1p\int _\Omega u^p\varphi '(v)\left (d_2 \Delta v+g(u,v)\right )\nonumber \\ =&\, -d_1(p-1)\int _\Omega u^{p-2}\varphi (v)|\nabla u|^2 -\int _\Omega \left (\frac {d_2}p\varphi ''(v)-\chi \varphi '(v)\right )u^p|\nabla v|^2\nonumber \\ &\quad \ +\underbrace {\int _\Omega \left [\chi (p-1)\varphi (v)-(d_1+d_2)\varphi '(v)\right ] u^{p-1}\nabla u\cdot \nabla v}_{=\!:\,I_9}\nonumber \\ &\quad \ +\underbrace {\int _\Omega u^{p-1}\varphi (v)\left (f(u,v)+\frac {\varphi '(v)}{p\varphi (v)} ug(u,v)\right )}_{=\!:\,I_{10}}\quad \text {for all }t\in (0,T_{\max }). \end{align}

Using Young’s inequality, for all $t\in (0,T_{\max })$ , we obtain

(4.9) \begin{align} I_9\leq \frac {d_1(p-1)}2\int _\Omega u^{p-2}\varphi (v)|\nabla u|^2+\int _\Omega \frac {\left [\left (\chi (p-1)\varphi (v)-(d_1+d_2)\varphi '(v)\right )\right ]^2}{2d_1(p-1)\varphi (v)}u^p|\nabla v|^2. \end{align}

For the term $I_{10}$ , it follows from (1.7) and (4.4) that

(4.10) \begin{align} I_{10}=&\, \int _\Omega u^{p-1}\varphi (v) \left \{\left (1+\alpha u\right )uv-u+\frac {\varphi '(v)}{p\varphi (v)} u\left [\sigma v\left (1-\frac vK\right )-\left (1+\alpha u\right )uv\right ]\right \}\nonumber \\ \leq &\, \int _\Omega u^{p-1}\varphi (v) \left [\left (1+\alpha u\right )uv+\frac {\sigma k_{\chi ,p}}p uv-\frac {\chi \Gamma _1}p\left (1+\alpha u\right )u^2v\right ]\nonumber \\ =&\, \int _\Omega \underbrace {u^{p}\varphi (v)v \left [-\left (1+\alpha u\right )\left (\frac {\chi \Gamma _1}p u-1\right )+\frac {\sigma k_{\chi ,p}}p\right ]} _{=\!:\,I_{11}}\quad \text {for all }t\in (0,T_{\max }). \end{align}

For the term $I_{11}$ , (2.1) and (4.7) imply that

\begin{align*} \int _{\left \{u\lt \frac p{\chi \Gamma _1}\right \}}I_{11} \leq |\Omega |\left (\frac p{\chi \Gamma _1}\right )^{p}\varphi (M_0) M_0 \left [ \left (1+\alpha \frac p{\chi \Gamma _1}\right ) +\frac {\sigma k_{\chi ,p}}p \right ]=\!:C_1 \end{align*}

and

\begin{align*} \int _{\left \{u\geq \frac p{\chi \Gamma _1}\right \}}I_{11} \leq \frac {\sigma k_{\chi ,p} M_0}p \int _{\left \{u\geq \frac p{\chi \Gamma _1}\right \}}u^{p}\varphi (v) \leq \frac {\sigma k_{\chi ,p} M_0}p \int _\Omega u^{p}\varphi (v), \end{align*}

which along with (4.10) indicates that

(4.11) \begin{align} I_{10} \leq C_1+\frac {\sigma k_{\chi ,p} M_0}p \int _\Omega u^{p}\varphi (v)\quad \text {for all }t\in (0,T_{\max }). \end{align}

The combination of (4.5), (4.8), (4.9) and (4.11) yields

(4.12) \begin{align} \frac 1p\frac d{dt}\int _\Omega u^p\varphi (v)+\frac 1p\int _\Omega u^p\varphi (v) +\frac {d_1(p-1)}2\int _\Omega u^{p-2}|\nabla u|^2\varphi (v) \leq &\, C_1+C_2 \int _\Omega u^{p}\varphi (v) \end{align}

for all $t\in (0,T_{\max })$ , where $C_2\,:\!=\, \frac {1+\sigma k_{\chi ,p} M_0}p$ . For the last term in the right-hand side of the above inequality, it follows from (4.7) and the Gagliardo–Nirenberg inequality that there exists a constant ${C_{GN}}\gt 0$ such that

(4.13) \begin{align} C_2 \int _{\Omega } u^p \varphi \left (v\right ) \leq C_2 \varphi (M_0) \|u^{\frac {p}{2}}\|_{L^2(\Omega )}^2 \leq C_2 \varphi (M_0) {C_{GN}} \left (\|\nabla u^{\frac {p}{2}}\|_{L^2(\Omega )}^{2 \theta }\|u^{\frac {p}{2}}\|_{L^{\frac {2}{p}}(\Omega )}^{2(1-\theta )}+\|u^{\frac {p}{2}}\|_{L^{\frac {2}{p}}(\Omega )}^2\right ), \end{align}

where

\begin{align*} \theta =\frac {\frac {p}{2}-\frac {1}{2}}{\frac {p}{2}-\frac {1}{2}+\frac {1}{n}} \in (0,1). \end{align*}

Since $0\lt 2\theta \lt 2$ , by Young’s inequality, (2.3) and (4.7), we obtain from (4.13) that

(4.14) \begin{align} C_2 \int _{\Omega } u^p \varphi \left (v\right ) \leq &\, C_2 \varphi (M_0) {C_{GN}}\left (M_1+1\right )^p\left (\|\nabla u^{\frac {p}{2}}\|_{L^2(\Omega )}^{2 \theta }+1\right )\nonumber \\ \leq &\, C_2 \varphi (M_0) {C_{GN}}\left (M_1+1\right )^p\left (\varepsilon _1\|\nabla u^{\frac {p}{2}}\|_{L^2(\Omega )}^2+\varepsilon _1^{-\frac \theta {1-\theta }}+1\right )\nonumber \\ \leq &\, \frac {d_1(p-1)}{p^2}\int _\Omega |\nabla u^{\frac {p}{2}}|^2\varphi (v)+C_3\nonumber \\ =&\, \frac {d_1(p-1)}4\int _\Omega u^{p-2}|\nabla u|^2\varphi (v)+C_3, \end{align}

where

\begin{align*} \varepsilon _1\,:\!=\, \frac {d_1(p-1)}{p^2C_2 \varphi (M_0) {C_{GN}}\left (M_1+1\right )^p} \quad \text {and}\quad C_3\,:\!=\, C_2 \varphi (M_0) {C_{GN}}\left (M_1+1\right )^p\left (\varepsilon _1^{-\frac \theta {1-\theta }}+1\right ). \end{align*}

Substituting (4.14) into (4.12), one has

\begin{align*} \frac 1p\frac d{dt}\int _\Omega u^p\varphi (v)+\frac 1p\int _\Omega u^p\varphi (v) +\frac {d_1(p-1)}4\int _\Omega u^{p-2}|\nabla u|^2\varphi (v) \leq &\, C_1+C_3. \end{align*}

Solving the above inequality, using (4.7), $u_0(x)\in W^{1,q}(\Omega )$ (recall $q\gt n$ ) and the Sobolev embedding $W^{1,q}(\Omega )\hookrightarrow L^p(\Omega )$ , we can find two positive constants $C_4$ and $C_5$ satisfying

\begin{align*} \|u(\cdot ,t)\|_{L^{p}(\Omega )}^p\leq \int _\Omega u^p\varphi (v)\leq \max \left \{\varphi (M_0)\|u_0\|_{L^{p}(\Omega )}^p,p(C_1+C_3)\right \} \leq C_4\left (\|u_0\|^p_{W^{1,q}(\Omega )}+1\right ) \leq C_5 \end{align*}

for all $t\in (0,T_{\max })$ . This completes the proof.

4.2.1 Proof of Theorem1.1 for $n\geq 3$

By taking $p =2n$ in (2.7) with Lemma 4.2, the result of Theorem1.1 for $n\geq 3$ follows from Lemma 2.6 and (2.2). $\square$

5. Linear stability analysis and spatiotemporal patterns

Spatiotemporal patterns are important to understand the population distribution of biological systems. Predator–prey systems with prey-taxis can produce spatial patterns, as observed in experiments [26]. For the corresponding temporal predator–prey system of (1.1) with (1.2), hunting cooperation (1.4) can induce the Hopf bifurcations [49, 52, 60]. By incorporating diffusions for both predator and prey species, the system has Turing instability only if there is hunting cooperation and the prey species spreads faster than the predator species ( $d_2\gt d_1$ ). Moreover, the prey-taxis system (1.5) has no Turing patterns without hunting cooperation (i.e., $\alpha =0$ ) [23]. We expect that prey-taxis will stabilize the system (1.5), similar to that observed in [37] and [62]. This section is devoted to investigating the effects of the interaction of prey-taxis and hunting cooperation (1.4) on the spatiotemporal distribution of population dynamics described by (1.6). We conduct linear stability analysis and perform numerical simulations to illustrate possible spatiotemporal patterns.

Without spatial structures, (1.6) becomes the following ordinary differential system

(5.1) \begin{align} \begin{cases} u_t=\left (1+\alpha u\right )uv-u,\\ v_t=\sigma v\left (1-\frac vK\right )-\left (1+\alpha u\right )uv, \end{cases} \end{align}

which has a trivial equilibrium $E_0=(0,0)$ , a prey-only equilibrium $E_1=(0,K)$ , and possible positive equilibria $E_*=(u_*, v_*)$ solving

(5.2) \begin{align} v_*=\frac 1{1+\alpha u_* }\lt 1\quad \text {and}\quad (1+\alpha u_* )u_*=\sigma \left (1-\frac { v_*}K\right ){,} \end{align}

where $u_*$ is the positive root of the equation

\begin{align*} \alpha ^2 u^3+2 K \alpha u^2+K(1-\alpha \sigma ) u+\sigma (1-K)=0. \end{align*}

The number of positive equilibria and the linear stability of equilibria of the system (5.1) are well studied in [49, 60]. We reorganize these results below. Define the positive constant

(5.3) \begin{align} K_*=\frac {27 \alpha \sigma }{2+9 \alpha \sigma +2(1+3 \alpha \sigma ) \sqrt {1+3 \alpha \sigma }}\in (0,1], \end{align}

and the sets

\begin{align*} \begin{cases} S_0=(0,1]\times (0,1]\cup (1,\infty )\times (0,K_*),\\ S_1=(0,\infty )\times (1,\infty )\cup (1,\infty )\times \{1\},\\ S_2=(1,\infty )\times (K_*,1). \end{cases} \end{align*}

Then the number of positive equilibria of the system (5.1), denoted by $N_0$ , is

\begin{align*} N_0=i\quad \text {if}\quad (\alpha \sigma ,K)\in S_i,\quad i=0,1,2. \end{align*}

To be precise, if $(\alpha \sigma ,K)\in S_2$ , then we denote the two positive equilibria by

\begin{align*} \left (u_{1 *}, v_{1 *}\right ) \text { and }\left (u_{2 *}, v_{2 *}\right )\text { with }0\lt v_{1*}\lt v_{2*}\lt 1. \end{align*}

It is easy to show that $E_0=(0,0)$ is a saddle, $E_1=({0, K})$ is linearly stable for $0\lt K\lt 1$ and linearly unstable for $K\gt 1$ (see also [49, Theorem 2.1]). The following results on the linear stability of positive equilibria hold.

Lemma 5.1 ([49, Theorem 2.2]). Let $K_*$ be given by (5.3) and

\begin{align*} &\widetilde K_*=\frac {1-3 \sigma +\sqrt {\sigma ^2+6 \sigma +1}}{2}\lt 1, \qquad \widetilde {\alpha }_*=\frac {(\widetilde {K}_*+\sigma )(\widetilde {K}_*+3\sigma )}{\sigma \widetilde {K}_*^2},\\ &\widetilde \alpha _1(K)=\frac {\left (K+\sigma \right )^2}{K^2(K+\sigma -1)}\ \text { if }K\gt 1-\sigma . \end{align*}

1. (a) If $(\alpha \sigma ,K)\in S_1$ , then the unique positive equilibrium $(u_*, v_*)$ is linearly stable for $0\lt \alpha \lt \widetilde {\alpha }_1(K)$ and linearly unstable for $\alpha \gt \widetilde {\alpha }_1(K)$ , and the Hopf bifurcation arise from $(u_*, v_*)$ at $\alpha =\widetilde {\alpha }_1(K)$ .

2. (b) If $(\alpha \sigma ,K)\in S_2$ , then $\left (u_{2 *}, v_{2 *}\right )$ is linearly unstable, and for $\left (u_{1 *}, v_{1 *}\right )$ , we have the following results:

   1. (i) when $\sigma \geq \frac {3}{2},\left (v_{1 *}, u_{1 *}\right )$ is linearly stable for $\alpha \lt \widetilde {\alpha }_1(K)$ and linearly unstable for $\alpha \gt \widetilde {\alpha }_1(K)$ , and the Hopf bifurcation arise from $(u_*, v_*)$ at $\alpha =\widetilde {\alpha }_1(K)$ ;

   2. (ii) when $0\lt \sigma \lt \frac {3}{2}$ and $\widetilde {K}^*\lt K\lt 1$ , $\left (v_{1 *}, u_{1 *}\right )$ is linearly stable for $\alpha \lt \widetilde {\alpha }_1(K)$ and linearly unstable for $\alpha \gt \widetilde {\alpha }_1(K)$ , and the Hopf bifurcation arise from $(u_*, v_*)$ at $\alpha =\widetilde {\alpha }_1(K)$ ;

   3. (iii) when $0\lt \sigma \lt \frac {3}{2}$ and $K_*\lt K\leq \widetilde {K}^*\lt 1$ , $\left (v_{1 *}, u_{1 *}\right )$ is linearly unstable.

5.1 Linear stability of (1.6)

Let $\Omega \subset \mathbb {R}^{n}$ ( $n\geq 1$ ) be a bounded domain with smooth boundary. We now consider the stability of the constant steady states $E_1=(0, K)$ and $E_*=(u_*,v_*)$ in the presence of spatial structure. The linearized system of (1.6) at a constant steady state $E_s=(u_s,v_s)$ is

(5.4) \begin{align} \begin{cases} \Phi _t=\mathcal {A} \Delta \Phi +\mathcal {J} \Phi , & x \in \Omega , t\gt 0, \\ (\nu \cdot \nabla ) \Phi =0, & x \in \partial \Omega , t\gt 0, \\ \Phi (x, 0)=(u_0-u_s, v_0-v_s)^{\mathcal {T}}, & x \in \Omega , \end{cases} \end{align}

where $\mathcal {T}$ represents the transpose,

\begin{align*} \Phi =\left (\begin{array}{c} u-u_s \\ v-v_s \end{array}\right ), \quad \mathcal {A}=\left ( \begin{array}{cc} d_1 & -\chi u_s \\ 0 & d_2 \end{array} \right ), \end{align*}

and $\mathcal {J}$ is the corresponding Jacobian matrix of (5.1)

(5.5) \begin{align} \mathcal {J}= \left ( \begin{array}{cc} 2 \alpha u_s v_s+v_s-1 & u_s (1+\alpha u_s) \\ -v_s (1+2 \alpha u_s) & \sigma \left (1-\frac {2 v_s}{K}\right ) -u_s (1+\alpha u_s) \\ \end{array} \right ) =\!: \left ( \begin{array}{cc} J_1 & J_2 \\ J_3 & J_4 \\ \end{array} \right ). \end{align}

By the method of separation of variables, the linear system (5.4) has solutions in the form of

\begin{align*} \Phi (x, t)=\sum _{k\geq 0} (U_k,V_k)^{\mathcal {T}} e^{\rho t} \psi _k(x), \end{align*}

where $\psi _k(x)$ is the eigenfunction of the Neumann eigenvalue problems

\begin{align*} \begin{cases} -\Delta \psi _k(x)=k^2 \psi _k(x), & x \in \Omega ,\\ \frac {\partial \psi _k(x)}{\partial \nu }=0, & x \in \partial \Omega , \end{cases} \end{align*}

with the wave number $k$ , the coefficients $U_k$ and $V_k$ are given by $(U_k,V_k)^{\mathcal {T}}=\int _\Omega \Phi (x, 0) \psi _k(x)$ , and $\rho$ is the temporal eigenvalue satisfying

\begin{align*} \rho \psi _k(x)=-k^2 \mathcal A \psi _k(x)+\mathcal J \psi _k(x). \end{align*}

Using the above equation and the fact that the sequence $\left \{\psi _k\right \}_{k\geq 0}$ forms an orthonormal basis of $L^2(\Omega )$ , we know that the two eigenvalues of the matrix $\mathcal M_k\,:\!=\, -k^2\mathcal A+\mathcal J$ are the roots of the equation

(5.6) \begin{align} \rho ^2+P_k\rho +Q_k=0, \end{align}

where

(5.7) \begin{align} P_k\,:\!=\, (d_1+d_2)k^2-\beta _1,\quad Q_k\,:\!=\, d_1d_2k^4-\beta _2 k^2+\beta _3, \end{align}

and

(5.8) \begin{align} \beta _1\,:\!=\, J_1+J_4,\quad \beta _2\,:\!=\, d_1J_4+d_2J_1+\chi u_s J_3,\quad \beta _3\,:\!=\, J_1J_4-J_2J_3. \end{align}

The two roots of (5.6), denoted by $\rho _\pm$ , are given by

(5.9) \begin{align} \rho _\pm =\frac {-P_k\pm \sqrt {\Delta _{\rho _k}}}{2},\quad \Delta _{\rho _k}\,:\!=\, P_k^2-4Q_k. \end{align}

Lemma 5.2. Let $E_s=(u_s,v_s)$ be a constant steady state of the system (1.6) . Then the following stability results hold for $E_s=(u_s,v_s)$ .

1. (i) $E_s$ is linearly stable if and only if $\min \limits _{k^2\geq 0}\left \{P_k,Q_k\right \}\gt 0$ , and $E_s$ is linearly unstable if and only if $\min \limits _{k^2\geq 0}\left \{P_k,Q_k\right \}\lt 0$ .

2. (ii) Turing instability arises if and only if

   \begin{align*} \beta _1\lt 0,\ \beta _2\gt 0,\ \beta _3\gt 0 \quad \text {and}\quad \min \limits _{k^2\gt 0} Q_k\lt 0. \end{align*}

Proof. The first conclusion is obvious. Moreover, Turing instability arises if $E_s=(u_s,v_s)$ is linearly stable in the ODE system (5.1), while it is unstable in the PDE system (1.6). Therefore, Turing instability occurs if and only if

\begin{align*} P_0=-\beta _1\gt 0,\ Q_0=\beta _3\gt 0 \quad \text {and}\quad \min \limits _{k^2\gt 0} \left \{P_k, Q_k\right \}\lt 0. \end{align*}

Given (5.7), $d_1,d_2\gt 0$ and $k^2\geq 0$ , the second conclusion of Lemma 5.2 follows immediately and the proof is completed.

At $E_1=(0, K)$ , $J_1=K-1$ , $J_2=0$ , $J_3=-K\lt 0$ , $J_4=-\sigma \lt 0$ , which imply

\begin{align*} \beta _1=K-1-\sigma ,\ \beta _2=-d_1\sigma +d_2(K-1),\ \beta _3=-\sigma (K-1). \end{align*}

This indicates $\beta _2\lt 0$ if $\beta _3\gt 0$ . Therefore, Lemma 5.2(ii) implies that Turing instability can never arise from $E_1=(0, K)$ . We next investigate whether Turing instability can arise from the positive constant steady state $E_*=(u_*, v_*)$ . Let $E_s=E_*$ . Then (5.5) implies

\begin{align*} \mathcal {J} = \left ( \begin{array}{cc} \alpha u_* v_* & \frac { u_*}{ v_*}\\ v_* -2 & -\frac { \sigma v_*}{K} \\ \end{array} \right ) = \left ( \begin{array}{cc} 1- v_* & \sigma \left (1-\frac { v_*}K\right )\\ v_* -2 & -\frac { \sigma v_*}{K} \\ \end{array} \right ). \end{align*}

Clearly, (5.2) implies $J_3=v_*-2\lt -1$ . Hence (5.8) implies

\begin{align*} \lim _{\chi \rightarrow +\infty }\beta _2=\lim _{\chi \rightarrow +\infty }(d_1J_4+d_2J_1+\chi u_s J_3)=-\infty . \end{align*}

This alongside Lemma 5.2(ii) indicates that Turing instability cannot arise from $E_*=(u_*,v_*)$ for sufficiently large $\chi \gt 0$ .

5.2 Spatiotemporal patterns

In this subsection, we present a specific example to illustrate the above analysis. For definiteness, we let

(5.10) \begin{align} \sigma =K=d_1=1,\ d_2\gt 0,\ \alpha \gt 0,\ \chi \geq 0. \end{align}

Then the system (1.6) with (5.10) admits a unique positive constant steady state $E_*=(u_*,v_*)=\left (\frac {\sqrt {\alpha }-1}{\alpha },\frac {1}{\sqrt {\alpha }}\right )$ if and only if $\alpha \gt 1$ . Hence, we assume $\alpha \gt 1$ below. We have

(5.11) \begin{align} P_k=(d_2+1) k^2-\beta _1,\quad Q_k=d_2 k^4-\beta _2 k^2+\beta _3, \end{align}

where

(5.12) \begin{align} \beta _1=1-\frac {2}{\sqrt {\alpha }},\ -\beta _2=\frac {1-(\sqrt \alpha -1)d_2}{\sqrt \alpha }+\frac {(2 \alpha -3 \sqrt {\alpha }+1)\chi }{\alpha ^{3/2}},\quad \beta _3=\frac {2 \left (\sqrt {\alpha }-1\right )^2}{\alpha }\gt 0. \end{align}

Therefore, steady-state bifurcations occur if

(5.13) \begin{align} \beta _2\gt 0,\quad \varphi _Q\left (\frac {\beta _2}{2d_2}\right )\lt 0 \quad \text {and}\quad Q_k\lt 0\text { for some wave number }k. \end{align}

where

(5.14) \begin{align} \varphi _Q(s)\,:\!=\, d_2 s^2-\beta _2 s+\beta _3,\quad s\geq 0. \end{align}

Lemma 5.3. Let $d_2\gt 0$ , $\alpha \gt 1$ , $\chi \geq 0$ , (5.12) hold, and the function $\varphi _Q(s)$ be given by (5.14) . Then $\beta _2\gt 0$ and $\varphi _Q\left (\frac {\beta _2}{2d_2}\right )\lt 0$ if and only if

(5.15) \begin{align} d_2\gt d^*_+(\alpha ,\chi ), \end{align}

where

(5.16) \begin{align} d^*_+(\alpha ,\chi )\,:\!=\, -\frac {D_1}{2D_2}-\frac {\sqrt {D_1^2-4D_2D_0}}{2D_2}\gt 0 \end{align}

with

(5.17) \begin{align} \begin{cases} D_2= -\left (\sqrt {\alpha }-1\right )^2 \alpha ^2\lt 0,\\ D_1= 2 \alpha \left (\sqrt {\alpha }-1\right ) \left [\left (2 \alpha -3 \sqrt {\alpha }+1\right ) \chi +\alpha \left (4 \sqrt {\alpha }-3\right )\right ]\gt 0,\\ D_0= -\left [\alpha +(2 \alpha -3 \sqrt {\alpha }+1)\chi \right ]^2\lt 0. \end{cases} \end{align}

Moreover, it holds that

(5.18) \begin{align} \frac {\partial }{\partial \alpha }d^*_+(\alpha ,\chi )\lt 0,\quad \frac {\partial }{\partial \chi }d^*_+(\alpha ,\chi )\gt 0, \end{align}

and for any $\chi \geq 0$ , we have

(5.19) \begin{align} \lim _{\alpha \rightarrow 1}d^*_+(\alpha ,\chi )=+\infty ,\quad \lim _{\alpha \rightarrow +\infty }d^*_+(\alpha ,\chi )=8. \end{align}

Proof. Let $d_2\gt 0$ , $\alpha \gt 1$ and $\chi \geq 0$ . First, (5.12) implies that $\beta _2\gt 0$ is equivalent to

(5.20) \begin{align} d_2\gt \frac {\alpha +(2 \alpha -3 \sqrt {\alpha }+1)\chi }{\alpha (\sqrt \alpha -1)}=\!:d^*_{\alpha ,\chi }\gt 0. \end{align}

Moreover, we have from (5.12), (5.14) and (5.17) that

(5.21) \begin{align} \varphi _Q\left (\frac {\beta _2}{2d_2}\right ) = \frac {\varphi _*(d_2)}{4 \alpha ^3 d_2}, \qquad \varphi _*(d_2)\,:\!=\, D_2d_2^2+D_1d_2+D_0. \end{align}

It holds that $D_1^2-4D_2D_0=32 \left (\sqrt {\alpha }-1\right )^3 \left (2 \sqrt {\alpha }-1\right ) \alpha ^3 \left ((\sqrt {\alpha }-1) \chi +\alpha \right )\gt 0$ , which along with $D_2,D_0\lt 0$ and $D_1\gt 0$ implies that $\varphi _*(d_2)$ has two positive roots

(5.22) \begin{align} d^*_-(\alpha ,\chi )\,:\!=\, -\frac {D_1}{2D_2}+\frac {\sqrt {D_1^2-4D_2D_0}}{2D_2},\quad d^*_+(\alpha ,\chi )\,:\!=\, -\frac {D_1}{2D_2}-\frac {\sqrt {D_1^2-4D_2D_0}}{2D_2}. \end{align}

Clearly, $0\lt d^*_-(\alpha ,\chi )\lt d^*_+(\alpha ,\chi )$ . We next show that

(5.23) \begin{align} 0\lt d^*_-(\alpha ,\chi )\lt d^*_{\alpha ,\chi }\lt d^*_+(\alpha ,\chi ). \end{align}

Indeed, (5.17) and (5.22) imply

\begin{align*} d^*_{\alpha ,\chi }+\frac {D_1}{2D_2} =&\,\frac {\alpha +(2 \alpha -3 \sqrt {\alpha }+1)\chi }{\alpha (\sqrt \alpha -1)}+ \frac {2 \alpha \left (\sqrt {\alpha }-1\right ) \left [\left (2 \alpha -3 \sqrt {\alpha }+1\right ) \chi +\alpha \left (4 \sqrt {\alpha }-3\right )\right ]}{-2\left (\sqrt {\alpha }-1\right )^2 \alpha ^2}\\ =&\,\frac {\alpha +(2 \alpha -3 \sqrt {\alpha }+1)\chi -\left [\left (2 \alpha -3 \sqrt {\alpha }+1\right ) \chi +\alpha \left (4 \sqrt {\alpha }-3\right )\right ]}{\alpha (\sqrt \alpha -1)} =-4, \end{align*}

hence

\begin{align*} -2 D_2\left(d_{\alpha, \chi}^*+\frac{D_1}{2 D_2}\right)=8D_2\lt 0\lt \sqrt {D_1^2-4D_2D_0}. \end{align*}

This along with (5.17) indicates $d^*_{\alpha ,\chi }\lt d^*_+(\alpha ,\chi )$ . Moreover, it follows from (5.17) and (5.22) that

\begin{align*} \left(D_1^2-4 D_2 D_0\right)-\left[2 D_2\left(d_{\alpha, \chi}^*+\frac{D_1}{2 D_2}\right)\right]^2 =&\, 32 \left (\sqrt {\alpha }-1\right )^3 \left (2 \sqrt {\alpha }-1\right ) \alpha ^3 \left ((\sqrt {\alpha }-1) \chi +\alpha \right )-(8D_2)^2\\ =&\,32 \left (\sqrt {\alpha }-1\right )^3 \alpha ^3 \left [\left (2 \alpha -3 \sqrt {\alpha }+1\right ) \chi +\alpha \right ]\\ \gt &\,0, \end{align*}

which proves $d^*_-(\alpha ,\chi )\lt d^*_{\alpha ,\chi }$ . Therefore, (5.23) holds. Then the combination of (5.20)-(5.23) proves that $\beta _2\gt 0$ and $\varphi _Q\left (\frac {\beta _2}{2d_2}\right )\lt 0$ if and only if (5.15) holds. We next prove (5.18). Since $D_2\lt 0$ , we have $d^*_+(\alpha ,\chi )=-\frac {D_1}{2D_2}+\frac 12\sqrt {\frac {D_1^2-4D_2D_0}{D_2^2}}$ . Elementary calculations (omitted for brevity) show that

\begin{align*} \frac {\partial }{\partial \alpha }\left (\frac {D_1}{D_2}\right )\gt 0,\quad \frac {\partial }{\partial \alpha }\left (\frac {D_1^2-4D_2D_0}{D_2^2}\right )\lt 0 \quad \text {and}\quad \frac {\partial }{\partial \chi }\left (\frac {D_1}{D_2}\right )\lt 0,\quad \frac {\partial }{\partial \chi }\left (\frac {D_1^2-4D_2D_0}{D_2^2}\right )\gt 0, \end{align*}

which along with $D_1^2-4 D_2 D_0\gt 0$ proves (5.18). Finally, we have

\begin{align*} \lim \limits _{\alpha \rightarrow +\infty } \frac {D_1}{D_2}=-8,\quad \lim \limits _{\alpha \rightarrow +\infty }\left (\frac {D_1^2-4D_2D_0}{D_2^2}\right )=64, \end{align*}

and it is elementary to show that $\lim \limits _{\alpha \rightarrow 1}d^*_+(\alpha ,\chi )=+\infty$ . Therefore (5.19) is obtained and the proof is completed.

Clearly, (5.11) implies that Hopf bifurcations may occur if $\alpha \geq 4$ and will never occur if $\alpha \in (1,4)$ . Lemma 5.3 alongside (5.13) indicates that steady-state bifurcations may occur if $d_2\gt d^*_+(\alpha ,\chi )$ . Moreover, steady-state bifurcations never occur if $0\lt d_2\lt d^*_+(\alpha ,\chi )$ . Indeed, the proof of Lemma 5.3 shows that $\beta _2\leq 0$ if $0\lt d_2\leq d^*_{\alpha ,\chi }$ and $\varphi _Q\left (\frac {\beta _2}{2d_2}\right )\gt 0$ if $d^*_{\alpha ,\chi }\lt d_2\lt d^*_+(\alpha ,\chi )$ , this alongside $d_2,\beta _3\gt 0$ and (5.11) indicates that $Q_k\gt 0$ for all $k$ if $0\lt d_2\lt d^*_+(\alpha ,\chi )$ , which means that a steady-state bifurcation is impossible. In particular, a Turing–Hopf bifurcation occurs (cf. [21, Definition 2.1 and Remark 2.2]) if $(\alpha ,d_2)=(4,d_\chi ^*)$ , where

(5.24) \begin{align} d_\chi ^*\,:\!=\, d^*_+(4,\chi )=5+\frac {3 \chi }{4}+ \sqrt {6(\chi +4)} \geq d_0^*=5+2\sqrt 6. \end{align}

The stability of $E_*$ can be fully classified in the $\alpha$ - $d_2$ plane, as the following.

Figure 1. (a) Plots of $\textrm{Re}(\rho _+)$ and $\textrm { Im}(\rho _+)$ in terms of $k^2\geq 0$ at $(\alpha ,d_2)=(4,d_1^*)$ , where $(d_1^*,k_1^*)\approx (11.2272,0.2110)$ is given by (5.24) and (5.25). (b) Turing-Hopf bifurcation diagram in the $\alpha$ - $d_2$ plane, where the state-steady bifurcation curve $d_2=d^*_+(\alpha ,1)$ given by (5.16) and the Hopf bifurcation curve $\alpha =4$ divide the region $\{(\alpha ,d_2):\ \alpha \gt 1, d_2\gt 0\}$ into four parts: I (stable), II (Hopf), III (Turing–Hopf) and IV (Turing). Other parameters are given by (5.10) with $\chi =1$ .

Figure 2. Numerical simulation of spatially homogeneous time-periodic patterns generated by the system (1.6) in the interval $\Omega =(0,20\pi )$ with (5.10), $\chi =1$ and $\alpha =d_2=5$ . The initial value $(u_0, v_0)$ is set as a small random perturbation of $(u_*,v_*)=(0.0.2472,0.4472)$ .

Lemma 5.4. Let (5.10) hold with $\alpha \gt 1$ and $d_\chi ^*$ be given by (5.24) . Then the system (1.6) has a unique positive constant steady state $E_*=(\frac {\sqrt {\alpha }-1}{\alpha },\frac {1}{\sqrt {\alpha }})$ . Moreover, the following conclusions hold.

1. (i) The ODE system (5.1) undergoes a Hopf bifurcation at $E_*=(\frac 14,\frac 12)$ when $\alpha =4$ .

2. (ii) The system (1.6) undergoes a Turing–Hopf bifurcation at $E_*=(\frac 14,\frac 12)$ when $(\alpha ,d_2)=(4,d_\chi ^*)$ . At $(\alpha ,d_2)=(4,d_\chi ^*)$ , $\mathcal {M}_k$ has the eigenvalues

   (5.25) \begin{align} \begin{cases} \rho _\pm =\pm \frac {i}{\sqrt {2}},\quad & \text {if }k^2=0,\\ \rho _-=-\frac {(8+3\chi ) \sqrt {6(\chi +4)}+12 \chi }{4 (3 \chi +4)}\lt 0,\ \rho _+=0,\quad & \text {if }k^2=k_\chi ^*\,:\!=\, \frac {2}{ \sqrt {6(\chi +4)}+4},\\ \textrm{Re}(\rho _\pm )\lt 0,\quad & \text {if }k^2\neq 0,k_\chi ^*. \end{cases} \end{align}

Proof. Let $k^2=0$ in (5.6). Then $\rho _\pm =\pm \sqrt {Q_0}=\pm \sqrt {\beta _3}=\pm \frac {i}{\sqrt {2}}$ . The first conclusion is obtained. We next prove (ii). When $(\alpha ,d_2)=(4,d_\chi ^*)$ , we have from the proof of Lemma 5.3 that $Q_k\geq 0$ , where “=” holds if and only if $k^2=\frac {\beta _2}{2d_2}=\frac {2}{ \sqrt {6(\chi +4)}+4}=k_\chi ^*$ . If $k^2=k_\chi ^*$ , then $P_k=(d_2+1) k^2=(d_\chi ^*+1) k_\chi ^*\gt 0$ and $Q_k=0$ , which implies $\rho _-=-P_k=-\frac {(8+3\chi ) \sqrt {6(\chi +4)}+12 \chi }{4 (3 \chi +4)}\lt 0$ and $\rho _+=0$ . If $k^2\neq 0,k_\chi ^*$ , then $P_k,Q_k\gt 0$ , and hence $\textrm{Re}(\rho _\pm )\lt 0$ .

Figure 3. Plots of $\textrm{Re}(\rho _+)$ and $\textrm { Im}(\rho _+)$ for $k=\frac m{20}$ with $m=0,1,2,\cdots$ under the parameter setting (5.10) with $\chi =1$ , $d_2=20$ and different values of $\alpha$ : (a) $\alpha =3$ ; (b) $\alpha =6$ ; (c) $\alpha =30$ .

Figure 4. Numerical simulations generated by the system (1.6) in the interval $\Omega =(0,20\pi )$ with (5.10), $\chi =1$ , $d_2=20$ , and different values of $\alpha$ shown in the columns: (a) $\alpha =3$ ; (b) $\alpha =6$ ; (c) $\alpha =30$ . The initial data $(u_0, v_0)$ are set as a small random perturbation of the positive constant steady state $(u_*,v_*)=(\frac {\sqrt {\alpha }-1}{\alpha },\frac {1}{\sqrt {\alpha }})$ .

For clarity, we shall first discuss a special case $\chi =1$ and then turn to the general case $\chi \geq 0$ .

Remark 5.1. Fig. 1 gives an illustration for Lemma 5.4 with $\chi =1$ . Since (5.9) implies $\textrm{Re}(\rho _-)\leq \textrm{Re}(\rho _+)$ , we only show the real and imaginary parts of $\rho _+$ in Fig. 1 (a). Moreover, in the $\alpha$ - $d_2$ plane, the Hopf bifurcation curve $\alpha =4$ and the steady state bifurcation curve $d_2=d^*_+(\alpha ,1)$ given by Lemma 5.3 divide the region $\left \{(\alpha ,d_2):\ \alpha \gt 1, d_2\gt 0\right \}$ into four parts: I (stable), II (Hopf), III (Turing–Hopf) and IV (Turing).

Figure 5. Turing–Hopf bifurcation diagram within the region $\{(\alpha ,d_2):\ \alpha \gt 1, d_2\gt 0\}$ in the $\alpha$ - $d_2$ plane for the system (1.6) with (5.10) and $\chi \in \{0,1,5,10\}$ . When $(\alpha , d_2)$ is above the steady state bifurcation curve $d_2=d^*_+(\alpha ,\chi )$ , which given by (5.16) is higher as $\chi \gt 0$ is larger, the Turing instability (resp. Hopf–Turing instability) will occur if $(\alpha , d_2)$ is left (resp. right) to the vertical Hopf bifurcation line $\alpha =4$ .

For numerical simulations, we let

\begin{align*} \Omega =(0,L)\quad \text {with}\quad L=20\pi . \end{align*}

We begin with the subcritical case $d_2\lt 8$ , for which the steady-state bifurcation will never occur (see Figure 1(b) where the horizontal asymptote $d_2=8$ is below the steady state bifurcation curve $d_2=d^*_+(\alpha ,1)$ ). Without loss of generality, we take $d_2=5$ . Then for $\alpha \in (0,4)$ , $(\alpha ,5)$ belongs to the stable region I. In this case, $(u_*,v_*)$ is linearly stable and no patterns will arise. For the Hopf region II, we take $\alpha =d_2=5$ , and then the system (1.6) can generate spatially homogeneous time-periodic patterns, as shown in Figure 2. We next consider the supercritical case $d_2\gt 8$ . Without loss of generality, we take $d_2=20$ and $\alpha =3,6,30$ . This gives

\begin{align*} (\alpha ,d_2)=(3,20)\in \textrm{IV},\quad {(\alpha ,d_2)=(6,20)\in \textrm{III},\quad (\alpha ,d_2)=(30,20)\in \textrm{III}}. \end{align*}

For these three values of $(\alpha ,d_2)$ , the eigenvalues $\rho _+$ of $\mathcal {M}_k$ with $k=\frac {m\pi }L=\frac m{20}$ , $m=0,1,2,\cdots$ , are shown in Figure 3, and the corresponding numerical simulations are shown in Figure 4. Clearly, for $(\alpha ,d_2)=(3,20)$ in the Turing instability region IV, Figure 3(a) shows that the eigenvalues of $\mathcal M_k$ with $m=0,1,2,3$ are a pair of complex conjugate numbers with negative real parts, and all other eigenvalues of $\mathcal M_k$ with $m=4,5,6,\cdots$ are real, and $\rho _+$ changes from negative to positive at $m=5$ . This indicates that there is a steady-state bifurcation but no Hopf bifurcation. The numerical simulations in Figure 4(a) show that non-constant stationary patterns arise, which is aligned with our analysis. For $(\alpha ,d_2)$ in the Turing–Hopf instability region III, the system (1.6) is unstable under temporal perturbation due to Hopf bifurcation and complex eigenvalues of $\mathcal M_0$ have positive real parts (see Figure 3(b) and Figure 3(c) at $m=0$ ), and the system (1.6) also undergoes Turing instability with real parts of an eigenvalue changing from negative to positive for some $k\neq 0$ (see Figure 3(b) with $m=6$ and Figure 3(c) with $m=7$ for instance). The interaction of Hopf and state-steady bifurcations will result in various complex patterns in the Turing–Hopf instability region III. For $(\alpha ,d_2)=(6,20)\in \textrm{III}$ , spatiotemporal periodic patterns are observed, as shown in Figure 4(b). However, if we keep the diffusion rate $d_2=20$ and increase the strength of hunting cooperation to $\alpha =30$ , then chaotic oscillatory patterns appear as shown in Figure 4(c).

The above analysis is performed for $\chi =1$ . We next consider $\chi \geq 0$ and explain how prey-taxis stabilizes the positive constant steady state. Clearly, (5.18) and (5.19) indicate that $d^*_+(\alpha ,\chi )$ strictly increases with respect to $\chi \in [0,+\infty )$ , with two asymptotes: $\alpha =1$ and $d_2=8$ in the $\alpha$ - $d_2$ plane. Figure 5 gives an illustration for $\chi =0,1,5,10$ , where the solid blue curve $d_2={d^*_+(\alpha ,1)}$ and the red vertical line $\alpha =4$ have been shown in Figure 1(b). The vertical coordinate of the intersection of $d_2={d^*_+(\alpha ,\chi )}$ and $\alpha =4$ is $d_\chi ^*$ given by (5.24). The steady state bifurcation curve $d^*_+(\alpha ,\chi )$ strictly increases in $\chi \in [0,+\infty )$ for $\alpha \gt 1$ (see (5.18)), while the area of stable region increases with respect to $\chi \in [0,+\infty )$ , see Figure 5 for $\chi \in \{0,1,5,10\}$ for instance. This indicates that prey-taxis plays a role in stabilizing the positive constant steady states.

As demonstrated above, the system (1.6) with $\chi \geq 0$ and $\alpha \gt 0$ may generate various spatiotemporal patterns, including spatially homogeneous time-periodic patterns, non-constant stationary patterns, spatially inhomogeneous time-periodic patterns and chaos. The former two types of patterns are also found in [6, 47, 60], where the latter two types of patterns result from the interaction of Hopf and state-steady bifurcations in the Turing–Hopf region. Moreover, we find that prey-taxis plays a stabilizing role in terms of pattern formation, which is similar to the case that the hunting cooperation functional response function adopts the Type II form [37, 62].