Impact of prey-taxis on a harvested intraguild predation predator–prey model

Abstract

In this paper, we report the spatiotemporal dynamics of an intraguild predation (IGP)-type predator–prey model incorporating harvesting and prey-taxis. We first discuss the local and global existence of the classical solutions in N-dimensional space. It is found that the model has a global classical solution when controlling the prey-taxis coefficient in a certain range. Thereafter, we focus on the existence of the steady-state bifurcation. Moreover, we theoretically investigate the properties of the bifurcating solution near the steady-state bifurcation critical threshold. As a consequence, the spatial pattern formation of this model can be theoretically confirmed. Importantly, by means of rigorous theoretical derivation, we provide discriminant criteria on the stability of the bifurcating solution. Finally, the complicated patterns are numerically displayed. It is demonstrated that the harvesting and prey-taxis significantly affect the pattern formation of this IGP-type predator–prey model. Our main results of this paper reveal that (i) The repulsive prey-taxis could destabilize the spatial homogeneity, while the attractive prey-taxis effect and self-diffusion will stabilize the spatial homogeneity of this model. (ii) Numerical results suggest that over-harvesting for prey or predators is not advisable, it can lead to an ecological imbalance due to a significant reduction in population numbers. However, harvesting within a certain range is a feasible approach.

1. Introduction

Intraguild predation (IGP) is ubiquitous in the natural environment, and it describes an interaction in which two or more species compete for shared resources and consume each other. Typically, prey promotes the growth of predator density due to the consumption of prey by predators. However, the impact of prey on predators resulting from competition for the same resource is rarely considered in some existing literature, see Refs. [7, 11, 16, 28]. As a consequence, there is interest in studying the dynamics of the predator–prey model with IG prey and IG predators. In fact, some scholars have devoted great attention to the dynamics of IGP-type predator–prey models. Ji et al. [19] investigated the well-posedness, properties of the solution semiflow, and spatiotemporal dynamics of a three-dimensional IGP-type predator–prey model with homogeneous Neumann boundary conditions. By employing a delayed IGP model, Shu et al. [34] demonstrated that delays could induce the stability switch, multitype bistability, and chaos phenomena. Bl $\acute {\textrm {e}}$ et al. [4] reported on the Hopf and Bautin bifurcations of an intraguild predation model with general functional responses for the predators and a significantly growing rate functions for the prey. The longtime behavior of solutions, the existence of biologically meaningful equilibria, and the linear and nonlinear stability of equilibria in an intraguild predator–prey model with a Holling type II functional response were investigated by Capone et al. in [5]. Please refer to Refs. [18, 29, 30, 31] for more experimental and theoretical results regarding IGP-type predator–prey models.

In this paper, we investigate the following IGP-type predator–prey model incorporating prey-taxis and linear harvesting:

(1) \begin{align} \left \{ {\begin{array}{l@{\quad}l} \displaystyle \frac {\partial P}{\partial t}=d_1\Delta P-\nabla \cdot (\xi \phi (P)\nabla Q)+P\left (\frac {bc}{cP+eQ}+dQ-\alpha \right ),\ &x\in \Omega, \ t\gt 0,\\[6pt] \displaystyle \frac {\partial Q}{\partial t}=d_2\Delta Q+Q\left (\frac {be}{cP+eQ}-dP-\beta \right )-hQ,\ &x\in \Omega, \ t\gt 0,\\[6pt] \displaystyle \frac {\partial P}{\partial \nu }=\frac {\partial Q}{\partial \nu }=0,\ &x\in \partial \Omega, \ t\geq 0,\\[6pt] P(x,0)=P_0(x)\geq 0,\ Q(x,0)=Q_{0}(x)\geq 0,\ &x\in \Omega, \end{array}} \right . \end{align}

where $P=P(x,t)$ and $Q=Q(x,t)$ are the densities of the IG predator and IG prey at spatial location $x$ and time $t$ , respectively. The domain $\Omega \subset \mathbb {R}^N$ is a bounded region with $N\geq 1$ , $\nu$ is the outward unit normal vector along the smooth $\partial \Omega$ , and $\Delta$ is the Laplacian operator. The parameters $d_1$ and $d_2$ describe the movement speeds of the predator $P$ and prey $Q$ , respectively. The terms $\frac {c}{cP+eQ}$ and $\frac {e}{cP+eQ}$ represent the per capita share of resources accruing to the predator $P$ and prey $Q$ , respectively. The parameter $b$ measures the consumption of the resources by predator and prey, while $\alpha$ and $\beta$ are the natural death rates of the predator $P$ and prey $Q$ , respectively. The term $-hQ$ explains the linear harvesting of the $Q$ species with the harvesting constant $h$ . Furthermore, the term $-\nabla \cdot (\xi \phi (P)\nabla Q)$ represents the prey-taxis with the sensitivity coefficient $\xi$ . This means that the predator species $P$ moves toward higher gradient directions of prey species $Q$ . The prey-taxis can be attractive or repulsive when $\xi \gt 0$ or $\xi \lt 0$ , respectively. $\phi (P)$ is a density functions related to population $P$ . This density function can take different forms. For instance, linear form: $\phi (P)=P$ , saturated form: $\phi (P)=\frac {P}{1+\epsilon P^m}$ with $\epsilon \gt 0$ and $m\geq 1$ , Ricker form: $\phi (P)=Pe^{-\epsilon P}$ with $\epsilon \gt 0$ , monotonic non-increasing form: $\phi (P)=\frac {1}{1+P}$ (or $\phi (P)=\frac {1}{(1+P)^2}$ ), among others. The parameters $b,e,c,d,h,\alpha, \beta, d_1,d_2$ are positive constants and prey-taxis sensitivity parameter $\xi \gt 0$ or $\xi \lt 0$ for its different biological meanings. We would like to mention that the prey-taxis term in the model (1) is similar to the chemotaxis term in some population models, see the references [21, 22, 32], for instance. When the prey-taxis coefficient $\xi =0$ and the harvesting constant $h=0$ , the model (1) degenerates into the classical IGP model, which was proposed by Holt and Polis in [14]. There are recent works focused on the dynamics of the IGP-type predator–prey model (1) with $\xi =0$ or $h=0$ . Ma et al. [25] reported spatiotemporal patterns in the model with delay and cross-fractional diffusion, showing that cross-fractional diffusion can induce Turing pattern formation. If choosing the density function $\phi (P)=P$ and the harvesting constant $h=0$ , Wang and Wang [35] showed the boundedness of classical solutions and the global stability of the positive equilibrium. The existence of global-in-time solutions and the Hopf bifurcation of the model with Schoener’s kinetic and indirect taxis have been reported by Mishra and Wrzosek in [26].

Let us state our tasks in this paper about the IGP-type predator–prey model (1). The first aim of this paper is to explore the solution profiles of the model (1). To be more specific, we want to study the local and global existence of the classical solution $(P(x,t),Q(x,t))$ in an $N$ -dimensional space. We can show that the IGP-type predator–prey model (1) admits a unique non-negative local-in-time classical solution $(P(x,t),Q(x,t))\in [C([0,T_{max});\, W^{1,p}(\Omega ))\cap C^{2,1}(\overline {\Omega }\times (0,T_{max}))]^2$ , with its maximal existence time $T_{max}$ by virtue of the Amann’s theorem [2]. The global existence of the classical solution $(P(x,t),Q(x,t))$ for the IGP-type predator–prey model (1) can be obtained by using estimates and the Neumann heat semigroup theory [13, 37]. Here, we can explain that the prey-taxis sensitivity coefficient $\xi$ can govern the global existence of the classical solution $(P(x,t),Q(x,t))$ . Our theoretical results show that if $ 0\lt \xi \leq \frac {d_1d_2}{3c_0(2+N)C_1(d_1+d_2)},$ where $C_1=\mathrm {max}\left \{\|Q_0(x)\|_{L^\infty (\Omega )},\frac {b}{h+\beta }\right \}$ is valid, then the IGP-type predator–prey model (1) possesses a unique non-negative global classical solution $(P(x,t),Q(x,t))\in [C([0,\infty );\, W^{1,p}(\Omega ))\cap C^{2,1} (\overline {\Omega }\times (0,\infty ))]^2$ and $\|P(\cdot, t)\|_{L^{\infty }(\Omega )} +\|Q(\cdot, t)\|_{L^{\infty }(\Omega )} \leq M$ , where $M$ is a positive constant dependents on $P_0(x)$ and $Q_{0}(x)$ for $P_{0}(x),\ Q_0(x)\geq 0(\not \equiv 0)$ .

Using bifurcation theory, the exploration of spatiotemporal pattern formation in ecological models is still a hot research area. Consequently, our next task is to explore the existence of steady-state bifurcation and the stability of the bifurcating solutions for the spatial local system of the system (1) when $\Omega =(0,L\pi )$ . This is

(2) \begin{align} \left \{ {\begin{array}{l@{\quad}l} \displaystyle d_1\Delta P-\nabla \cdot (\xi \phi (P)\nabla Q)+P\left (\frac {bc}{cP+eQ}+dQ-\alpha \right )=0,\ &x\in \Omega, \\[6pt] \displaystyle d_2\Delta Q+Q\left (\frac {be}{cP+eQ}-dP-\beta \right )-hQ=0,\ &x\in \Omega, \\[6pt] \displaystyle \frac {\partial P}{\partial \nu }=\frac {\partial Q}{\partial \nu }=0,\ &x\in \partial \Omega . \end{array}} \right . \end{align}

One difficulty is how to determine the stability of these bifurcating solutions of the system (2). Typically, scholars have adopted some existing techniques to investigate the stability of the bifurcating solutions. For instance, they use weakly nonlinear analysis method (or multiple time scale) [27, 3] and normal form theory [20, 12]. In these approaches, the authors derived the amplitude equations and normal forms so that the stability of the bifurcating solution can be established. In contrast to the previously mentioned technique, we will apply the Crandall–Rabinowitz local bifurcation theory [6, 8, 9, 33, 36] to demonstrate the existence and stability of the bifurcating solution (i.e., the nonconstant steady state) around the threshold of the steady-state bifurcation. By choosing the prey-taxis sensitivity coefficient $\xi$ as the steadystate bifurcation parameter, we can theoretically demonstrate that the repulsive prey-taxis ( $\xi \lt 0$ ) could destabilize the spatial homogeneity of this IGP-type predator–prey model, while the attractive prey-taxis ( $\xi \gt 0$ ) effect will stabilize the spatial homogeneity. Naturally, we conduct extensive numerical simulations to confirm our theoretical results by choosing different density functions $\phi (P)$ . For example, considering linear form $\phi (P)=P$ , saturated form $\phi (P)=\frac {P}{1+P}$ , and the Ricker form $\phi (P)=Pe^{-P}$ , we can observe the pattern formations in 1D and 2D domains, and on spherical and torus surfaces. We also investigate the influence of the harvesting effects on pattern formation. It is shown that extensive harvesting of IG prey will lead to the disappearance of spatial patterns. This phenomenon reminds us that over-harvesting for prey or predators is not advisable because of the drastic reduction in population numbers from the point of view of ecology.

In this paper, we require $(P_{0}(x),Q_0(x))$ and the density function $\phi (P)$ to fulfill the following.

1. (H1) $(P_{0}(x),Q_0(x))\in [W^{1,p}(\Omega )]^2$ with $p\gt N$ and $P_{0}(x),\ Q_{0}(x)\geq 0\ (\not \equiv 0)$ .

2. (H2) There is a $c_0$ such that $\phi (P)\leq c_0P$ for $\forall P\geq 0$ and $x\in \overline {\Omega }$ , where $\phi :[0,\infty )\rightarrow [0,\infty )$ is continuously differentiable and $\phi (0)=0$ . Moreover, we suppose that

3. (H3) $\frac {\beta +h}{ed}\lt \frac {b}{e\alpha -c(\beta +h)}\lt \frac {\alpha }{cd}$ .

Now we can release our main results of this article. The first result is concerned with the global existence of the classical solution $(P(x,t),Q(x,t))$ of the system (1) with the assumptions (H1) and (H2).

Theorem 1.1. (Global existence of the classical solution) Suppose $\Omega \subset \mathbb {R}^N$ with the smooth boundary $\partial \Omega$ and the initial conditions $(P_0(x),Q_{0}(x))\in [W^{1,p}(\Omega )]^2$ with $p\gt N$ and $P_0(x)\geq 0,Q_{0}(x)\geq 0$ for $x\in \overline {\Omega }$ . If

\begin{align*} 0\lt \xi \leq \frac {d_1d_2}{3c_0(2+N)C_1(d_1+d_2)}, \end{align*}

where $C_1=\mathrm {max}\left \{\|Q_0(x)\|_{L^\infty (\Omega )},\frac {b}{h+\beta }\right \}$ , then system (1) enjoys a unique global solution $(P(x,t),Q(x,t))\in [C([0,\infty );\,W^{1,p}(\Omega ))\cap C^{2,1} (\overline {\Omega }\times (0,\infty ))]^2$ and

\begin{equation*}\|P(\cdot, t)\|_{L^{\infty }(\Omega )} +\|Q(\cdot, t)\|_{L^{\infty }(\Omega )} \leq M,\end{equation*}

where $M$ is a positive constant depending on $P_0(x)$ and $Q_{0}(x)$ for $P_{0}(x),Q_0(x)\geq 0(\not \equiv 0)$ .

Remark 1.1. Theorem 1.1 shows the global existence of the classical solutions when the prey-taxis is attractive-type. We shall point out that a similar result could also be obtained as the prey-taxis is repulsive-type in (1).

Our following goal is to explore the existence and stability of the bifurcating solution induced by the steady state bifurcation. We need to mention that if the spatial dimension is high, namely, $N\geq 2$ , the analysis of bifurcation is very difficult, especially to discuss the stability criterion of the bifurcating solution. Therefore, to finish our goal, we restrict $N=1$ and choose $\Omega =(0,L\pi )$ with $L\gt 0$ .

If the assumption (H3) is true and fix $\Omega =(0,L\pi )$ with $L\gt 0$ , then system (1) has a unique positive equilibrium $E_{*}=(P_{*},Q_{*})=\left (\frac {be}{e\alpha -c(\beta +h)}-\frac {\beta +h}{d},\frac {\alpha }{d}-\frac {bc}{e\alpha -c(\beta +h)}\right )$ . Define

(3) \begin{align} \xi _k^S=\frac {d_1d_2\delta _k^4-(f_Pd_2 +g_Qd_1)\delta _k^2+d^2P_{*}Q_{*}}{\phi (P^{*})g_P\delta _k^2}\lt 0,\ k\in \mathbb {N}_0=\{0,1,2,\cdot \cdot \cdot \}, \end{align}

where $\delta _k=\frac {k}{L}\gt 0$ and $f_P=-\frac {bc^2P_{*}}{(cP_{*}+eQ_{*})^2}, g_P=-dQ_{*}-\frac {bceQ_{*}}{(cP_{*}+eQ_{*})^2}, g_Q= -\frac {be^2Q_{*}}{(cP_{*}+eQ_{*})^2}$ . Also, if there is a $k_0\in \mathbb {N}_0\backslash \{0\}$ satisfying

\begin{align*} k_0= \left \{ {\begin{array}{l@{\quad}l} \left [\hat {k}_0\right ]+1, \ &\textrm {if}\ \xi _{[\hat {k}_0]}^S\leq \xi _{[\hat {k}_0]+1}^S,\\[5pt] \left [\hat {k}_0\right ], \ &\textrm {if}\ \xi _{[\hat {k}_0]}^S\gt \xi _{[\hat {k}_0]+1}^S \end{array}} \right . \end{align*}

with $\hat {k}_0=L\sqrt {d \sqrt {\frac {P_{*}Q_{*}}{d_1d_2}}}$ , then $\xi _k^S$ has its maximum $\xi _{k_0}^S$ at $k=k_0$ , where $[\cdot ]$ is the integer function.

In this fashion, we can establish the stability result of the constant steady state $E_{*}$ .

Theorem 1.2. (Local stability of the constant steady state $E_{*}$ ) Suppose that (H2)–(H3) are satisfied and take $\Omega =(0,L\pi )$ with $L\gt 0$ .

1. (i) If $\xi \geq 0$ , $E_{*}$ is locally asymptotically stable;

2. (ii) If $\xi =\xi _k^S$ , then system (1) suffers from the steady-state bifurcation. Moreover, $E_{*}$ is locally asymptotically stable as $\xi _k^S\lt \xi \lt 0$ and it becomes unstable when $\xi \lt \xi _k^S\lt 0$ ;

3. (iii) If $0\lt \xi ^2\lt \frac {4d_1d_2Q_{*}}{c_0^2P_{*}C_1^2}$ , then $E_{*}$ is globally asymptotically stable.

Remark 1.2. Clearly, if $k=k_0$ , then system (1) will undergo the steady-state bifurcation at the threshold $\xi =\xi _{k_0}^S$ . We will later discuss the existence and stability of the nonconstant steady state (bifurcating solution) at this onset.

Remark 1.3. From (i)–(ii) of Theorem1.2, we infer that the repulsive prey-taxis (i.e., $\xi \lt 0$ ) could destabilize the spatial homogeneity of the IGP predator–prey model (1). On the contrary, the attractive prey-taxis effect (i.e., $\xi \gt 0$ ) and self-diffusion (i.e., $\xi =0$ ) will stabilize the spatial homogeneity.

Our third result implies that system (2) exhibits nonconstant steady state around $(P_{*},Q_{*},\xi _k^S)$ for $k\in \mathbb {N}_0\backslash \{0\}$ in $\mathbf {X}=\{u\in H^2(0,L\pi )|u^{\prime }(0)=u^{\prime }(L\pi )=0\}$ . To do so, define

(4) \begin{align} \mathcal {F}(P,Q,\xi )=\left ( \begin{array}{c@{\quad}c@{\quad}c} \displaystyle (d_1P^{\prime }-\xi \phi (P) Q^{\prime })^{\prime }+P\left (\frac {bc}{cP+eQ}+dQ-\alpha \right )\\[3pt] \displaystyle d_2 Q^{\prime \prime }+Q\left (\frac {be}{cP+eQ}-dP-\beta \right )-hQ \end{array} \right ) \end{align}

and the Fr $\acute {\textrm {e}}$ chet derivative $D_{(P,Q)}\mathcal {F}(\breve {P},\breve {Q},\xi )(P,Q)$ of the operator $\mathcal {F}(P,Q,\xi )$ . Then, for any $(\breve {P},\breve {Q},\xi )\in \mathbf {X}\times \mathbf {X}\times \mathbb {R}$ , we deduce

(5) \begin{align} &D_{(P,Q)}\mathcal {F}(\breve {P},\breve {Q},\xi )(P,Q)\\[3pt] \nonumber \displaystyle = & \left ( \begin{array}{c@{\quad}c@{\quad}c} d_1P^{\prime \prime }-\xi (\phi ^{\prime }(\breve {P})P \breve {Q}^{\prime }+\phi (\breve {P}) Q^{\prime })^{\prime } +\left [\frac {bce\breve {Q}}{(c\breve {P}+e\breve {Q})^2}+d\breve {Q}-\alpha \right ]P +\left [d\breve {P}-\frac {bce\breve {P}}{(c\breve {P}+e\breve {Q})^2}\right ]Q\\[6pt] \displaystyle d_2 Q^{\prime \prime }-\left [d\breve {Q}+\frac {bce\breve {Q}}{(c\breve {P}+e\breve {Q})^2}\right ]P +\left [\frac {bce\breve {P}}{(c\breve {P}+e\breve {Q})^2}-d\breve {P}-\beta -h\right ]Q \end{array} \right ). \end{align}

Let $(\breve {P},\breve {Q},\xi )=(P_{*},Q_{*},\xi )$ , we obtain

(6) \begin{align} &D_{(P,Q)}\mathcal {F}(P_{*},Q_{*},\xi )(P,Q)\\[3pt] \nonumber \displaystyle =& \left ( \begin{array}{c@{\quad}c@{\quad}c} d_1P^{\prime \prime }-\xi (\phi (P_{*}) Q^{\prime })^{\prime } +\left [\frac {bceQ_{*}}{(cP_{*}+eQ_{*})^2}+dQ_{*}-\alpha \right ]P +\left [dP_{*}-\frac {bceP_{*}}{(cP_{*}+eQ_{*})^2}\right ]Q\\[6pt] \displaystyle d_2 Q^{\prime \prime }-\left [dQ_{*}+\frac {bceQ_{*}}{(cP_{*}+eQ_{*})^2}\right ]P +\left [\frac {bceP_{*}}{(cP_{*}+eQ_{*})^2}-dP_{*}-\beta -h\right ]Q\\[6pt] \displaystyle \end{array} \right )\\[6pt] \displaystyle \nonumber =&\left ( \begin{array}{c@{\quad}c@{\quad}c} d_1P^{\prime \prime }-\xi (\phi (P_{*}) Q^{\prime })^{\prime } -\frac {bc^2P_{*}}{(cP_{*}+eQ_{*})^2}P +\left [dP_{*}-\frac {bceP_{*}}{(cP_{*}+eQ_{*})^2}\right ]Q\\[6pt] \displaystyle d_2 Q^{\prime \prime }-\left [dQ_{*}+\frac {bceQ_{*}}{(cP_{*}+eQ_{*})^2}\right ]P -\frac {be^2Q_{*}}{(cP_{*}+eQ_{*})^2}Q \end{array} \right ). \end{align}

We can establish the following.

Theorem 1.3. (Existence of the nonconstant steady state) Suppose that (H1)–(H3) are satisfied and take $\Omega =(0,L\pi )$ with $L\gt 0$ , $\xi _j^S\neq \xi _k^S$ for $j\neq k$ and $k\in \mathbb {N}_0\setminus \{0\}$ , where $\xi _k^S$ is given by (3). Then system (2) admits a spatially inhomogeneous solution which resulted from $(P_{*},Q_{*})$ when $\xi =\xi _k^S$ for $k\in \mathbb {N}_0\setminus \{0\}$ . Moreover, in the vicinity of the onset $(P_{*},Q_{*},\xi _k^S)$ , there exists a bifurcation branch $\mathcal {S}_k(\varepsilon )=(P_{k}(\varepsilon, x),Q_{k}(\varepsilon, x))$ that satisfies

(7) \begin{align} \left \{ {\begin{array}{*{20}{l}} \xi _k^S(\varepsilon )=\xi _k^S+\mathcal {O}(\varepsilon ),\\[3pt] (P_{k}(\varepsilon, x),Q_{k}(\varepsilon, x))=(P_{*},Q_{*}) +\varepsilon (\widehat {P}_{k},\widehat {Q}_{k})+\mathcal {O}(\varepsilon ) \end{array}} \right . \end{align}

for any $\varepsilon \in ({-}\varrho, \varrho )$ and $\varrho$ is a small positive constant. Also, $(P_{k}(\varepsilon, x),Q_{k}(\varepsilon, x))-(P_{*},Q_{*}) -\varepsilon (\widehat {P}_{k},\widehat {Q}_{k})=\mathcal {O}(\varepsilon )\in \mathcal {K}$ with $\mathcal {K}$ is a closed complement of $\mathcal {N}(D_{(P, Q)}\mathcal {F}(P_{*},Q_{*},\xi ))$ and it admits

(8) \begin{align} \mathcal {K}=\left \{(P,Q)\in \mathbf {X}\times \mathbf {X}\Big | \int _0^{L\pi }(P\widehat {P}_{k}+Q\widehat {Q}_{k})dx=0\right \}, \end{align}

where $\mathcal {N}$ is null space and

(9) \begin{align} \widehat {P}_{k}=\mathrm {cos}\frac {k x}{L},\ \widehat {Q}_{k}=\alpha _k\mathrm {cos}\frac {k x}{L} \end{align}

with

\begin{align*} \alpha _k=-\frac {dQ_{*}(cP_{*}+eQ_{*})^2+bceQ_{*}}{be^2Q_{*}+d_2\delta _k^2(cP_{*}+eQ_{*})^2}\lt 0,\ \ k\in \mathbb {N}_0\backslash \{0\}. \end{align*}

Benefiting form (7) of Theorem1.3, we can set $\xi _k^S(\varepsilon )=\xi _k^S+\varepsilon \xi _1+\varepsilon ^2\xi _2+\cdot \cdot \cdot, $ where $\xi _1$ and $\xi _2$ are undetermined constants. Let $\xi _{k_0}^S=\mathrm {max}_{k\in \mathbb {N}_0\backslash \{0\}}\xi _k^S$ . Accordingly, our fourth result shows that $\xi _1=0$ and the sign of $\xi _2$ uniquely determines the stability of the bifurcating solution $(P_{k_0}(\varepsilon, x),Q_{k_0}(\varepsilon, x))$ for $\varepsilon \in ({-}\varrho, \varrho )$ .

Theorem 1.4. (Local stability of the nonconstant steady state) Suppose that the conditions (H1)–(H3) hold and fix $\Omega =(0,L\pi )$ with $L\gt 0$ . Then we can compute the first perturbation term $\xi _1=0$ in $\xi _{k}^S(\varepsilon )$ . In addition, when $k=k_0$ , near $(P_{*},Q_{*},\xi _{k_0}^S)$ , the bifurcating solution $\mathcal {S}_{k_0}(\varepsilon )=(P_{k_0}(\varepsilon, x),Q_{k_0}(\varepsilon, x))$ is asymptotically stable when $\xi _2\lt 0$ and it is unstable as $\xi _2\gt 0$ for $\varepsilon \in ({-}\varrho, \varrho )$ .

Remark 1.4. The results presented in Theorem1.4 show that the stability of the bifurcating solution (namely, nonconstant steady state) completely depends on the symbol of the second perturbation term $\xi _2$ in $\xi _{k}^S(\varepsilon )$ for $\varepsilon \in ({-}\varrho, \varrho )$ .

2. Existence and boundedness of classical solution

Lemma 2.1. Suppose that $\Omega \subset \mathbb {R}^N$ with the smooth boundary $\partial \Omega$ and $(P_0(x),Q_0(x))\in [W^{1,p}(\Omega )]^2$ with $p\gt N$ fulfilling $P_0(x)\geq 0,Q_{0}(x)\geq 0$ for $x\in \overline {\Omega }$ . Then, we can yield the following.

1. (i) System ( 1 ) enjoys a unique nonnegative classical solution $(P(x,t),Q(x,t))$ satisfying $(P(x,t),Q(x,t))\in [C([0,T_{max});\, W^{1,p}(\Omega ))\cap C^{2,1}(\overline {\Omega }\times (0,T_{max}))]^2$ . Also, we have

   (10) \begin{align} P(x,t)\gt 0,\ Q(x,t)\leq C_1,\ x\in \overline {\Omega },\ t\in [0,T_{max}) \end{align}
   where $C_1=\mathrm {max}\left \{\|Q_0(x)\|_{L^\infty (\Omega )},\frac {b}{h+\beta }\right \}$ and $T_{max}\gt 0$ implies that the maximal existence time.

2. (ii) There are $C_2\gt 0$ and $C_3\gt 0$ such that

   \begin{align*} \|Q(x,t)\|_{L^1(\Omega )}\leq C_2,\ \|P(x,t)\|_{L^1(\Omega )}\leq C_3,\ t\in (0,T_{max}), \end{align*}
   where
   \begin{align*} C_2=\mathrm {max}\left \{\int _\Omega Q_0(x)dx,\frac {b|\Omega |}{h+\beta }\right \},\ C_3=\mathrm {max}\left \{\int _\Omega \left (P_0(x)+Q_{0}(x)\right )dx,\frac {b|\Omega |}{\min \{\alpha, h+\beta \}}\right \}. \end{align*}

3. (iii) If for any $T\gt 0$ , there exists $C(T)$ such that

   \begin{align*} \sup _{0\leq t\leq T}\|P(x,t),Q(x,t)\|_{L^\infty (\Omega )}\leq C(T),\ 0\lt T\lt \min \{1,T_{max}\}, \end{align*}
   then there holds $T_{max}=+\infty$ , where $C(T)$ depends on $T$ and $\|P_{0}(x),Q_0(x)\|_{W^{1,p}(\Omega )}$ .

Proof. The local-in-time existence of the nonnegative classical solution $(P(x,t),Q(x,t))$ in (i) can be confirmed by employing Amann’s theorem [2]. Next, using the $P-$ equation of (1), we obtain

\begin{align*} \left \{ {\begin{array}{l@{\quad}l} \frac {\partial P}{\partial t}=d_1\Delta P-\xi \phi ^{\prime }(P)\nabla P\nabla Q-\xi \phi (P)\Delta Q+P\Gamma (P,Q),\ \ &x\in \Omega, \ t\in (0,T_{max}),\\[3pt] \frac {\partial P}{\partial \nu }=0,\ \ &x\in \partial \Omega, \ t\in (0,T_{max}),\\[3pt] P(x,0)=P_{0}(x)\geq 0,\ \ &x\in \Omega, \end{array}} \right . \end{align*}

where $\Gamma (P,Q)=\frac {bc}{cP+eQ}+dQ-\alpha$ . It follows from the maximum principle that $0$ is a lower solution for the above equation. Thus, it follows that $P(x,t)\geq 0$ for all $(x,t)\in \Omega \times (0,T_{max})$ . Using the strong maximum principle and the initial data $P_{0}(x)\geq 0(\not \equiv 0)$ , one can claim that $P(x,t)\gt 0$ is true. Next, from the $Q-$ equation of (1), one can derive

\begin{align*} \left \{ {\begin{array}{l@{\quad}l} \displaystyle \frac {\partial Q}{\partial t}-d_2\Delta Q\leq Q\left (\frac {b}{Q}-\beta \right )-hQ,\ \ &x\in \Omega, \ t\gt 0,\\[6pt] \displaystyle \frac {\partial Q}{\partial \nu }=0,\ \ &x\in \partial \Omega, \ t\geq 0,\\[6pt] Q(x,0)=Q_{0}(x)\geq 0,\ \ &x\in \Omega . \end{array}} \right . \end{align*}

Therefore, the maximum principle gives that $Q(x,t)\leq \frac {b}{h+\beta }$ for any $(x,t)\in \Omega \times (0,T_{max})$ . For (ii), integrating the $Q-$ equation of (1) over $\Omega$ , we get

\begin{align*} \frac {d }{d t}\int _\Omega Qdx\leq b|\Omega |-(h+\beta )\int _\Omega Qdx. \end{align*}

Accordingly, one has

\begin{align*} \int _\Omega Qdx\leq \mathrm {max}\left \{\int _\Omega Q_0(x)dx,\frac {b|\Omega |}{h+\beta }\right \}\,:\!=\, C_2. \end{align*}

On the other hand, we can obtain

\begin{align*} \frac {d }{d t}\int _\Omega \left (P+Q\right )dx= b|\Omega |-\alpha \int _\Omega Pdx-(h+\beta )\int _\Omega Qdx \leq b|\Omega |-\min \{\alpha, (h+\beta )\}\int _\Omega (P+ Q)dx. \end{align*}

This gives

\begin{align*} \int _\Omega Pdx\leq \mathrm {max}\left \{\int _\Omega \left (P_0(x)+Q_{0}(x)\right )dx,\frac {b|\Omega |}{\min \{\alpha, h+\beta \}}\right \}\,:\!=\, C_3. \end{align*}

Finally, conclusion (iii) can be directly obtained by using Theorem 15.5 in [1]. This ends the proof.

Lemma2.1 shows the local-in-time existence of the classical solution $(P(x,t),Q(x,t))$ , our following goal is exploring its global existence. To obtain the global existence of the classical solution $(P(x,t),Q(x,t))$ , we introduce some existing results.

Lemma 2.2. (Lemma2.6 of [22]) Suppose that $z(t)$ satisfies

\begin{align*} \left \{ {\begin{array}{*{20}{l}} z^{\prime }(t)\leq -a_1z^\ell (t)+a_2z(t)+a_3,\\[3pt] z(0)=z_0\gt 0, \end{array}} \right . \end{align*}

where $a_1,a_2,a_3\gt 0$ and $\ell \gt 1$ . Then,

\begin{align*} z(t)\leq \mathrm {max}\{C_4(z_0),C_5(a_1,a_2,a_3,\ell )\}. \end{align*}

Lemma 2.3. For $n\gt 1$ , $p(x)\geq 0$ and $q(x)\geq 0$ , the following inequality holds

\begin{align*} \int _\Omega p^{n-1}\varphi (q)dx\leq C_6\int _\Omega p^n\varphi (q)dx+C_7, \end{align*}

where $C_6$ and $C_7$ are positive constants and $\varphi (q)$ is bounded with respect to $q$ .

Proof. By employing $\varepsilon -$ Young inequality, we get

\begin{align*} \int _\Omega p^{n-1}\varphi (q)dx \leq &\varepsilon \int _\Omega (p^{n-1}\varphi (q))^{\frac {n}{n-1}}dx+C_\varepsilon |\Omega | \\[3pt]=&\varepsilon \int _\Omega (\varphi (q))^{\frac {1}{n-1}}(p^{n}\varphi (q))dx+C_\varepsilon |\Omega | \leq C_6\int _\Omega p^n\varphi (q)dx+C_7. \end{align*}

This ends the proof.

Lemma 2.4. (Lemma2.3 of [22]) Suppose $m\in \{0,1\},p\in [1,\infty )$ , and $q\in (1,\infty )$ . Then, there is a $C_{8}\gt 0$ such that

(11) \begin{align} \|u\|_{W^{m,p}(\Omega )}\leq C_{8}\|({-}\Delta +1)^k u\|_{L^{q}(\Omega )}, \end{align}

for $u\in D(({-}\Delta +1)^k)$ with ${D(({-}\Delta +1)^k)}=\left \{\zeta \in W^{2,p}(\Omega ):\zeta _\nu =0\ \textrm {over}\ \partial \Omega \right \}$ and $k\in (0,1)$ satisfies

\begin{align*} m-\frac {N}{p}\lt 2k-\frac {N}{q}. \end{align*}

In addition, if $q\geq p$ is satisfied, then $C_9\gt 0$ and $\gamma \gt 0$ exist such that

(12) \begin{align} \|({-}\Delta +1)^k e^{-t({-}\Delta +1)}u\|_{L^{q}(\Omega )}\leq C_9t^{-k-\frac {N}{2}\left (\frac {1}{p}-\frac {1}{q}\right )}e^{-\gamma t}\|u\|_{L^{q}(\Omega )}, \end{align}

for $u\in L^{p}(\Omega )$ , where the diffusion semigroup $\{e^{-t({-}\Delta +1)}\}_{t\geq 0}$ maps $L^{p}(\Omega )$ into $D(({-}\Delta +1)^k)$ . Moreover, for any $p\in (1,\infty )$ and $\varepsilon \gt 0$ , there are $C_{10}\gt 0$ and $\mu \gt 0$ satisfying

(13) \begin{align} \|({-}\Delta +1)^k e^{t\Delta }\nabla \cdot u\|_{L^{q}(\Omega )}\leq C_{10}t^{-k-\frac {1}{2}-\varepsilon } e^{-\mu t}\|u\|_{L^{p}(\Omega )} \end{align}

for $u\in L^{p}(\Omega )$ .

Now we can prove the following results.

Lemma 2.5. Assume that $\Omega \subset \mathbb {R}^N$ with the smooth boundary $\partial \Omega$ , the initial condition $(P_0(x),Q_{0}(x))\in [W^{1,p}(\Omega )]^2$ with $p\gt N$ and $P_0(x)\geq 0,\ Q_{0}(x)\geq 0$ for $x\in \overline {\Omega }$ . If

\begin{align*} 0\lt \xi \leq \frac {d_1d_2}{3c_0(2+N)C_1(d_1+d_2)}, \end{align*}

then there is a positive constant $C_{11}$ such that

(14) \begin{align} \|P(\cdot, t)\|_{L^{N+2}(\Omega )}\leq C_{11},\ t\in (0,T_{max}). \end{align}

Proof. Let $n=N+2$ and define an auxiliary function $\varphi (Q)=e^{(\sigma Q)^2}$ for $0\leq Q(x,t)\leq C_1$ , where $\sigma$ satisfies

\begin{equation*}\sigma =\frac {1}{C_1(d_1+d_2)}\sqrt {\frac {d_1d_2(n-1)}{6n}}\gt 0.\end{equation*}

Accordingly, multiplying $P-$ equation by $P^{n-1}\varphi (Q)$ and integrating it over $\Omega$ , one yields

\begin{align*} &\frac {1}{n}\frac {d}{dt}\int _\Omega P^n\varphi (Q)dx\\[3pt]=&\int _\Omega P^{n-1}\varphi (Q)P_tdx+\frac {1}{n}\int _\Omega P^n\varphi ^{\prime }(Q)Q_tdx \\[3pt] =&d_1\int _\Omega P^{n-1}\varphi (Q)\Delta Pdx-\int _\Omega P^{n-1}\varphi (Q)\nabla (\xi \phi (P)\nabla Q)dx\\[3pt]&+\int _\Omega P^{n-1}\varphi (Q)\left (\frac {bc P}{cP+eQ}+dPQ-\alpha P\right )dx\\[3pt]&+\frac { d_2}{n}\int _\Omega P^n\varphi ^{\prime }(Q)\Delta Qdx+\frac {1}{n}\int _\Omega P^n\varphi ^{\prime }(Q)\left [\frac {beQ}{cP+eQ}-dPQ-(\beta +h)Q\right ]dx\\[3pt]\leq &-d_1(n-1)\int _\Omega P^{n-2}\varphi (Q)|\nabla P|^2dx-d_1\int _\Omega P^{n-1}\varphi ^{\prime }(Q)\nabla P\nabla Qdx\\[3pt]&+\xi \int _\Omega P^{n-1}\varphi ^{\prime }(Q)\phi (P)|\nabla Q|^2dx+\xi (n-1)\int _\Omega P^{n-2}\varphi (Q)\phi (P)\nabla P\nabla Qdx \\[3pt]&+b\int _\Omega P^{n-1}\varphi (Q)dx+dC_1\int _\Omega P^n\varphi (Q)dx-\frac { d_2}{n}\int _\Omega P^n\varphi ^{\prime \prime }(Q)|\nabla Q|^2dx \\[3pt]&-d_2\int _\Omega P^{n-1}\varphi ^{\prime }(Q)\nabla P\nabla Qdx+\frac {2b\sigma ^2 C_1}{n}\int _\Omega P^{n-1}\varphi (Q) dx. \end{align*}

Recalling the assumption $\phi (P)\leq c_0P$ in (H2), we have

(15) \begin{align} &\frac {1}{n}\frac {d}{dt}\int _\Omega P^n\varphi (Q)dx+d_1(n-1)\int _\Omega P^{n-2}\varphi (Q)|\nabla P|^2dx+\frac { d_2}{n}\int _\Omega P^n\varphi ^{\prime \prime }(Q)|\nabla Q|^2dx \\[3pt] \nonumber \leq &-(d_1+d_2)\int _\Omega P^{n-1}\varphi ^{\prime }(Q)\nabla P\nabla Qdx+c_0\xi \int _\Omega P^{n}\varphi ^{\prime }(Q)|\nabla Q|^2dx\\[3pt]&\nonumber +c_0\xi (n-1)\int _\Omega P^{n-1}\varphi (Q)\nabla P\nabla Qdx +dC_1\int _\Omega P^n\varphi (Q)dx \\[3pt]&\nonumber +\left (b+\frac {2b\sigma ^2 C_1}{n}\right )\int _\Omega P^{n-1}\varphi (Q) dx. \end{align}

Now, by employing Young’s inequality and note that $P^{n-1}=P^{\frac {n-2}{2}}P^{\frac {n}{2}}$ , one obtains

\begin{align*} &-(d_1+d_2)\int _\Omega P^{n-1}\varphi ^{\prime }(Q)\nabla P\nabla Qdx \\[3pt] \leq &(d_1+d_2)\int _\Omega P^{n-1}\varphi ^{\prime }(Q)|\nabla P||\nabla Q|dx \\[3pt] =&{\int _\Omega \left (\sqrt {\frac {(n-1)d_1\varphi (Q)}{2}}P^{\frac {n-2}{2}}|\nabla P|\right )\left (\frac {\sqrt {2}(d_1+d_2)}{\sqrt {(n-1)d_1\varphi (Q)}}P^{\frac {n}{2}}\varphi ^{\prime }(Q)|\nabla Q|\right )dx} \\[3pt] \leq &\frac {(n-1)d_1}{4}\int _\Omega P^{n-2}\varphi (Q)|\nabla P|^2dx +\frac {(d_1+d_2)^2}{(n-1)d_1}\int _\Omega \frac {P^n\varphi ^{\prime 2}(Q)}{\varphi (Q)}|\nabla Q|^2dx \end{align*}

and, similarly, one yields

\begin{align*} &c_0\xi (n-1)\int _\Omega P^{n-1}\varphi (Q)\nabla P\nabla Qdx \\[3pt] \leq &\frac {(n-1)d_1}{4}\int _\Omega P^{n-2}\varphi (Q)|\nabla P|^2dx +\frac {c_0^2\xi ^2(n-1)}{d_1}\int _\Omega P^n\varphi (Q)|\nabla Q|^2dx. \end{align*}

Consequently, putting these into (15), we get

(16) \begin{align} &\frac {1}{n}\frac {d}{dt}\int _\Omega P^n\varphi (Q)dx+\frac {d_1(n-1)}{2}\int _\Omega P^{n-2}\varphi (Q)|\nabla P|^2dx+\frac { d_2}{n}\int _\Omega P^n\varphi ^{\prime \prime }(Q)|\nabla Q|^2dx\\[3pt] \nonumber \leq &\frac {(d_1+d_2)^2}{(n-1)d_1}\int _\Omega \frac {P^n\varphi ^{\prime 2}(Q)}{\varphi (Q)}|\nabla Q|^2dx+c_0\xi \int _\Omega P^n\varphi ^{\prime }(Q)|\nabla Q|^2dx\\[3pt]&\nonumber +\frac {c_0^2\xi ^2(n-1)}{d_1}\int _\Omega P^n\varphi (Q)|\nabla Q|^2dx +dC_1\int _\Omega P^n\varphi (Q)dx \\[3pt]&\nonumber +\left (b+\frac {2b\sigma ^2 C_1}{n}\right )\int _\Omega P^{n-1}\varphi (Q) dx. \end{align}

Let

\begin{align*} &\omega _1(Q)=\frac {(d_1+d_2)^2}{(n-1)d_1}\frac {\varphi ^{\prime 2}(Q)}{\varphi (Q)},\ \omega _2(Q)=c_0\xi \varphi ^{\prime }(Q),\\[3pt] &\omega _3(Q)=\frac {c_0^2\xi ^2(n-1)}{d_1}\varphi (Q),\ \omega _4(Q)=\frac { d_2}{n}\varphi ^{\prime \prime }(Q). \end{align*}

As a consequence, (16) becomes

(17) \begin{align} &\frac {1}{n}\frac {d}{dt}\int _\Omega P^n\varphi (Q)dx+\frac {d_1(n-1)}{2}\int _\Omega P^{n-2}\varphi (Q)|\nabla P|^2dx+\int _\Omega P^n\omega _4(Q)|\nabla Q|^2dx\\[3pt] \nonumber \leq &\int _\Omega P^n\omega _1(Q)|\nabla Q|^2dx+\int _\Omega P^n\omega _2(Q)|\nabla Q|^2dx+\int _\Omega P^n\omega _3(Q)|\nabla Q|^2dx \\[3pt]&\nonumber +dC_1\int _\Omega P^n\varphi (Q)dx +\left (b+\frac {2b\sigma ^2 C_1}{n}\right )\int _\Omega P^{n-1}\varphi (Q) dx. \end{align}

Recalling that $\varphi (Q)=e^{(\sigma Q)^2}$ for $0\leq Q(x,t)\leq C_1$ , one obtains

\begin{align*} &\omega _1(Q)=\frac {4\sigma ^4(d_1+d_2)^2Q^2}{(n-1)d_1}\varphi (Q),\ \omega _2(Q)=2\sigma ^2Qc_0\xi \varphi (Q),\\[3pt] &\omega _3(Q)=\frac {c_0^2\xi ^2(n-1)}{d_1}\varphi (Q),\ \omega _4(Q)=\frac { d_2}{n}(2\sigma ^2\varphi (Q)+4\sigma ^4Q^2\varphi (Q)). \end{align*}

For $0\leq Q(x,t)\leq C_1$ , we take

\begin{align*} 0\lt \xi \leq \frac {d_1d_2}{3c_0nC_1(d_1+d_2)}. \end{align*}

Using this approach, we have

\begin{align*} &\frac {3\omega _1(Q)}{\omega _4(Q)}\leq \frac {6n\sigma ^2(d_1+d_2)^2Q^2}{(n-1)d_1d_2}\leq \frac {6n\sigma ^2(d_1+d_2)^2C_1^2}{(n-1)d_1d_2}=1,\\[3pt] &\frac {3\omega _2(Q)}{\omega _4(Q)}\leq \frac {d_1}{d_1+d_2}\lt 1, \end{align*}

and

\begin{align*} \frac {3\omega _3(Q)}{\omega _4(Q)}\leq \frac {3c_0^2\xi ^2n(n-1)}{2d_1d_2\sigma ^2}\leq 1. \end{align*}

Therefore, one has

\begin{align*} \int _\Omega (\omega _1(Q)+\omega _2(Q)+\omega _3(Q))P^n|\nabla Q|^2dx \leq \int _\Omega \omega _4(Q) P^n|\nabla Q|^2dx. \end{align*}

As such, (17) takes the form

(18) \begin{align} &\frac {1}{n}\frac {d}{dt}\int _\Omega P^n\varphi (Q)dx+\frac {d_1(n-1)}{2}\int _\Omega P^{n-2}\varphi (Q)|\nabla P|^2dx \\[3pt] &\nonumber \leq dC_1\int _\Omega P^n\varphi (Q)dx +\left (b+\frac {2b\sigma ^2 C_1}{n}\right )\int _\Omega P^{n-1}\varphi (Q) dx. \end{align}

In light of Lemma2.3, we get

(19) \begin{align} &\frac {1}{n}\frac {d}{dt}\int _\Omega P^n\varphi (Q)dx+\frac {d_1(n-1)}{2}\int _\Omega P^{n-2}\varphi (Q)|\nabla P|^2dx \\[3pt] &\nonumber \leq (dC_1+C_{12})\int _\Omega P^n\varphi (Q)dx+C_{13}. \end{align}

On the other hand, owing to $\varphi (Q)=e^{(\sigma Q)^2}$ , so $\varphi (Q)\leq e^{\sigma ^2C_1^2}$ for $0\leq Q(x,t)\leq C_1$ . Thereby, by utilizing the first two inequalities on page 55 of [15] and (ii) of Lemma2.1, we get

\begin{align*} \int _\Omega P^n\varphi (Q)dx\leq & e^{\sigma ^2C_1^2}\int _\Omega P^ndx=e^{\sigma ^2C_1^2}\left \|P^{\frac {n}{2}}\right \|_{L^2(\Omega )}^2 \\[3pt] \leq &e^{\sigma ^2C_1^2}C\left \|P^{\frac {n}{2}}\right \|_{W^{1,2}(\Omega )}^{2\nu } \left \|P^{\frac {n}{2}}\right \|_{L^{\frac {2}{n}}(\Omega )}^{2(1-\nu )} \\[3pt] \leq &e^{\sigma ^2C_1^2}C\left (\left \|\nabla P^{\frac {n}{2}}\right \|_{L^{2}(\Omega )}+\left \|P^{\frac {n}{2}}\right \|_{L^{\frac {2}{n}}(\Omega )}\right )^{2\nu } \left \|P^{\frac {n}{2}}\right \|_{L^{\frac {2}{n}}(\Omega )}^{2(1-\nu )} \\[3pt] =&e^{\sigma ^2C_1^2}C\left (\left \|\nabla P^{\frac {n}{2}}\right \|_{L^{2}(\Omega )}+\left \|P\right \|_{L^{1}(\Omega )}^{\frac {n}{2}}\right )^{2\nu } \left \|P\right \|_{L^{1}(\Omega )}^{n(1-\nu )} \\[3pt] \leq &C_{14}\left (\left \|\nabla P^{\frac {n}{2}}\right \|_{L^{2}(\Omega )}^2+1\right )^{\nu } \end{align*}

where $C$ is positive constant and $ \nu =\frac {\frac {nN}{2}-\frac {N}{2}}{\frac {nN}{2}+1-\frac {N}{2}}\in (0,1).$ Accordingly, for $n\gt 2$ and $0\lt \nu \lt 1$ , we get

(20) \begin{align} \int _\Omega P^{n-2}\varphi (Q)|\nabla P|^2dx\geq \int _\Omega P^{n-2}|\nabla P|^2dx{=\frac {4}{n^2}\int _\Omega |\nabla P^{\frac {n}{2}}|^2dx}\geq \frac {4C_{14}^{-\frac {1}{\nu }}}{n^2} \left (\int _\Omega P^n\varphi (Q)dx\right )^{\frac {1}{\nu }}-\frac {4}{n^2}. \end{align}

Thereby, putting (20) into (19), we have

\begin{align*} \frac {1}{n}\frac {d}{dt}\int _\Omega P^n\varphi (Q)dx \leq (dC_1+C_{12})\int _\Omega P^n\varphi (Q)dx -\frac {2d_1(n-1)C_{14}^{-\frac {1}{\nu }}}{n^2} \left (\int _\Omega P^n\varphi (Q)dx\right )^{\frac {1}{\nu }}+\frac {2d_1(n-1)}{n^2}. \end{align*}

By using Lemma2.2, there is a $C_{15}\gt 0$ such that

\begin{align*} \int _\Omega P^n\varphi (Q)dx\leq C_{15}. \end{align*}

This implies

\begin{align*} \|P(\cdot, t)\|_{L^n(\Omega )}\leq C_{11} \end{align*}

is valid. We finish the proof.

The following result means that the solution $P(x,t)$ admits the $L^\infty$ -bound.

Lemma 2.6. Suppose that $\Omega \subset \mathbb {R}^N$ with the smooth boundary $\partial \Omega$ , the initial conditions $(P_0(x),Q_{0}(x))\in [W^{1,p}(\Omega )]^2$ with $p\gt N$ and $P_0(x)\geq 0,\ Q_{0}(x)\geq 0$ for $x\in \overline {\Omega }$ . If

\begin{align*} 0\lt \xi \leq \frac {d_1d_2}{3c_0(2+N)C_1(d_1+d_2)}, \end{align*}

then there is a positive constant $C_{16}$ such that

(21) \begin{align} \|P(\cdot, t)\|_{L^{\infty }(\Omega )}\leq C_{16},\ t\in (0,T_{max}). \end{align}

Proof. Rewrite $Q-$ equation of (1) as follows.

\begin{align*} \frac {\partial Q}{\partial t}=d_2\Delta Q-Q+Q\left (\frac {be}{cP+eQ}-dP-\beta -h\right )+Q. \end{align*}

Then, we can compute

\begin{align*} Q(\cdot, t)=e^{-t({-}d_2\Delta +1)}Q_0+\int _0^te^{-(t-s)({-}d_2\Delta +1)} \left [Q\left (\frac {be}{cP+eQ}-dP-\beta -h\right )+Q\right ]ds. \end{align*}

Let $\tau \in (0,T_{max}),0\lt \tau \lt 1,q\gt N$ , $\frac {1}{2}(1+\frac {N}{q})\lt k\lt 1$ . Then, using (11), (12) in Lemma2.4 and (14) in Lemma2.5, one gets

\begin{align*} &\|Q(\cdot, t)\|_{W^{1,\infty }(\Omega )}\\[3pt] &\leq C_8{\left \|({-}d_2\Delta +1)^k {\left [e^{-t({-}d_2\Delta +1)}Q_0+\int _0^te^{-(t-s)({-}d_2\Delta +1)}(b+(dP+\beta +h+1)C_1)ds\right ]}\right \|}_{L^{q}(\Omega )}\\[3pt]&\leq C_8C_9t^{-k}e^{-\gamma t}\|Q_0\|_{L^{q}(\Omega )}+C_8C_9\int _0^t(t-s)^{-k}e^{-\gamma (t-s)}(b+(d\|P\|_{L^{q}(\Omega )}+\beta +h+1)C_1)ds \\[3pt]&\leq C_{16}t^{-k}+C_{17}\int _0^t(t-s)^{-k}e^{-\gamma (t-s)}ds \\[3pt]&\leq C_{16}t^{-k}+C_{17}\int _0^\infty {\varrho ^{-k}}e^{-\gamma \varrho }d\varrho \\[3pt]&\leq C_{16}\tau ^{-k}+C_{17}\Gamma (1-k)\,:\!=\, K(\tau ),\ t\in (\tau, T_{max}), \end{align*}

where $\Gamma (\cdot )$ is a Gamma function and $\Gamma (1-k)\gt 0$ due to $\frac {1}{2}(1+\frac {N}{q})\lt k\lt 1$ . Therefore, we get

(22) \begin{align} \|\nabla Q(\cdot, t)\|_{L^{\infty }(\Omega )}\leq K(\tau ),\ t\in (\tau, T_{max}). \end{align}

Now, the variation of the constant formula to the $P-$ equation of (1) shows

\begin{align*} P(\cdot, t)=&e^{-t({-}d_1\Delta +1)}P_{0}-\xi \int _0^te^{-(t-s)({-}d_1\Delta +1)}\nabla (\phi (P)\nabla Q)ds \\[3pt]&+\int _0^te^{-(t-s)({-}d_1\Delta +1)}\left [P\left (\frac {bc}{cP+eQ}+dQ-\alpha \right )+P\right ]ds. \end{align*}

Let

\begin{align*} P_1(\cdot, t)=e^{-t({-}d_1\Delta +1)}P_{0},P_2(\cdot, t) =-\xi \int _0^te^{-(t-s)({-}d_1\Delta +1)}\nabla (\phi (P)\nabla Q)ds, \end{align*}

and

\begin{align*}P_3(\cdot, t)=\int _0^te^{-(t-s)({-}d_1\Delta +1)}\left [P\left (\frac {bc}{cP+eQ}+dQ-\alpha \right )+P\right ]ds. \end{align*}

Therefore, $P(\cdot, t)=P_1(\cdot, t)+P_2(\cdot, t)+P_3(\cdot, t)$ . That is to say, we must give $L^\infty$ -bounds of $P_1(\cdot, t),P_2(\cdot, t)$ and $P_3(\cdot, t)$ to obtain $\|P(\cdot, t)\|_{L^{\infty }(\Omega )}$ .

Now, for $P_1(\cdot, t)$ , by using (11) and (12), we have

\begin{align*} \|P_1(\cdot, t)\|_{L^{\infty }(\Omega )}=&\|e^{-t({-}d_1\Delta +1)}P_{0}\|_{L^{\infty }(\Omega )} \\[3pt] \leq & C_8\|({-}d_1\Delta +1)^k e^{-t({-}d_1\Delta +1)}P_{0}\|_{L^{q}(\Omega )} \\[3pt] \leq &C_8C_9t^{-k}e^{-\gamma t}\|P_{0}\|_{L^{q}(\Omega )} \\[3pt] \leq &C_8C_9{\tau ^{-k}}e^{-\gamma t}\|P_{0}\|_{L^{q}(\Omega )} \\[3pt] \leq &C_{18}\|P_{0}\|_{L^{\infty }(\Omega )},\ \ t\in (\tau, T_{max}) \end{align*}

for $m=0,p=\infty, 0\lt \tau \lt 1,{\frac {N}{2q}\lt k\lt 1},q\gt N$ and $\gamma \gt 0$ .

For $P_2(\cdot, t)$ , one takes $\frac {N}{2q}\lt k\lt \frac {1}{2}$ , so $0\lt \varepsilon \lt 1/2-k$ . Then, by employing (11)–(13) and (22), we obtain

\begin{align*} \|P_2(\cdot, t)\|_{L^{\infty }(\Omega )}\leq & C_8{\left \|({-}d_1\Delta +1)^k \xi \int _0^te^{-(t-s)({-}d_1\Delta +1)}|\nabla (\phi (P)\nabla Q)|ds\right \|}_{L^{q}(\Omega )} \\[3pt] \leq &\xi C_{8} \int _0^t{\left \|({-}d_1\Delta +1)^k e^{-(t-s)({-}d_1\Delta +1)}|\nabla (\phi (P)\nabla Q)|\right \|}_{L^{q}(\Omega )}ds \\[3pt] \leq &C_{19} \int _0^t(t-s)^{-k-\varepsilon -1/2} e^{-(\mu +1)(t-s)}ds \\[3pt] \leq &C_{19} \int _0^\infty \varrho ^{-k-\varepsilon -1/2} e^{-(\mu +1)\varrho }d\varrho \\[3pt] \leq &C_{19} \Gamma (1/2-k-\varepsilon ),\ \ t\in (\tau, T_{max}), \end{align*}

where $\Gamma (1/2-k-\varepsilon )\gt 0$ due to $0\lt \varepsilon \lt 1/2-k$ .

For $P_3(\cdot, t)$ , in a similar approach, we have

\begin{align*} \|P_3(\cdot, t)\|_{L^{\infty }(\Omega )}\leq & C_{8} {\left \|({-}d_1\Delta +1)^k \int _0^te^{-(t-s)({-}d_1\Delta +1)}(b+(dC_1+1+\alpha )P)ds\right \|}_{L^{q}(\Omega )} \\[3pt] \leq & C_{8} \int _0^t{\left \|({-}d_1\Delta +1)^k e^{-(t-s)({-}d_1\Delta +1)}(b+(dC_1+1+\alpha )P)\right \|}_{L^{q}(\Omega )}ds \\[3pt] \leq &C_{8}C_{9} \int _0^t(t-s)^{-k} e^{-(t-s)\gamma }(b+(dC_1+1+\alpha )\left \|P\right \|_{L^{q}(\Omega )})ds \\[3pt] \leq & C_{20} \int _0^t(t-s)^{-k} e^{-(t-s)\gamma }ds \\[3pt] \leq & C_{20} \int _0^\infty \varrho ^{-k} e^{-\varrho \gamma }d\varrho \\[3pt] \leq & C_{20} \Gamma (1-k),\ \ t\in (\tau, T_{max}), \end{align*}

where $\Gamma (1-k)\gt 0$ since $0\lt k\lt 1$ . Therefore, the result performed in (21) is valid. The proof readily follows.

Proof of Theorem 1.1. By employing $Q(x,t)\leq C_1$ in Lemma2.1 and $\|P(\cdot, t)\|_{L^{\infty }(\Omega )}\leq C_{16}$ in Lemma2.6, where $C_{16}$ is a positive constant and $C_1=\mathrm {max}\left \{\|Q_0(x)\|_{L^\infty (\Omega )},\frac {b}{h+\beta }\right \}$ , we can infer that there exists a positive constant $M$ depending on $P_0(x)$ and $Q_{0}(x)$ for $P_{0}(x),Q_0(x)\geq 0(\not \equiv 0)$ such that $\|P(\cdot, t)\|_{L^{\infty }(\Omega )} +\|Q(\cdot, t)\|_{L^{\infty }(\Omega )} \leq M$ is fulfilled. The proof is finished.

3. Steady-state bifurcation

In this section, we shall establish the existence and stability of the nonconstant steady state resulting from the steady-state bifurcation near the positive equilibrium of the system (1). To achieve this, let

\begin{align*} \left \{ {\begin{array}{*{20}{l}} \displaystyle f(P,Q)=P\left (\frac {bc}{cP+eQ}+dQ-\alpha \right ),\\[3pt] \displaystyle g(P,Q)=Q\left (\frac {be}{cP+eQ}-dP-\beta \right )-hQ. \end{array}} \right . \end{align*}

3.1 Stability analysis

Taking $\Omega =(0,L\pi )$ with $L\gt 0$ . Then at $E_{*}$ , the linearization form of the system (1) is given by

(23) \begin{align} \left \{ {\begin{array}{*{20}{l}} \displaystyle \frac {\partial P}{\partial t}=d_1\Delta P-\xi \phi (P^{*})\Delta Q-\frac {bc^2P_{*}}{(cP_{*}+eQ_{*})^2}P +\left [dP_{*}-\frac {becP^{*}}{(cP_{*}+eQ_{*})^2}\right ]Q,\\[11pt] \displaystyle \frac {\partial Q}{\partial t}=d_2\Delta Q-\left [dQ_{*}+\frac {bceQ_{*}}{(cP_{*}+eQ_{*})^2}\right ]P -\frac {be^2Q_{*}}{(cP_{*}+eQ_{*})^2}Q. \end{array}} \right . \end{align}

Considering the eigenvalue problem

(24) \begin{align} \left \{ {\begin{array}{*{20}{l}} \displaystyle d_1\zeta _{xx}-\xi \phi (P^{*})\eta _{xx}-\frac {bc^2P_{*}}{(cP_{*}+eQ_{*})^2}\zeta +\left [dP_{*}-\frac {becP^{*}}{(cP_{*}+eQ_{*})^2}\right ]\eta =\lambda _k\zeta, \\[7pt] \displaystyle d_2\eta _{xx}-\left [dQ_{*}+\frac {bceQ_{*}}{(cP_{*}+eQ_{*})^2}\right ]\zeta -\frac {be^2Q_{*}}{(cP_{*}+eQ_{*})^2}\eta =\lambda _k\eta, \\[7pt] \displaystyle \frac {\partial \zeta }{\partial \nu }=\frac {\partial \eta }{\partial \nu }=0, \end{array}} \right . \end{align}

where $\lambda _k$ denotes the eigenvalue of the problem (24). For the no-flux boundary conditions, one takes the form of $(\zeta (x),\eta (x))$ as follows

\begin{align*} \zeta (x)=\sum _{k=0}^\infty a_k\mathrm {cos}\frac {k x}{L},\ \eta (x)=\sum _{k=0}^\infty b_k\mathrm {cos}\frac {k x}{L}, \end{align*}

where $a_k$ and $b_k$ are constants. By using (24), one has

\begin{align*} \sum _{k=0}^\infty (J_k-\lambda _k I)\left ( \begin{array}{c} a_k\\[3pt] b_k \end{array} \right )\mathrm {cos}\frac {k x}{L}=0, \end{align*}

where

\begin{align*} J_k=\left ( \begin{array}{c@{\quad}c} \displaystyle -\frac {bc^2P_{*}}{(cP_{*}+eQ_{*})^2}-d_1\delta _k^2 & \displaystyle dP_{*}-\frac {becP^{*}}{(cP_{*}+eQ_{*})^2}+\xi \phi (P^{*})\delta _k^2 \\[9pt] \displaystyle -dQ_{*}-\frac {bceQ_{*}}{(cP_{*}+eQ_{*})^2} & \displaystyle -\frac {be^2Q_{*}}{(cP_{*}+eQ_{*})^2}-d_2\delta _k^2 \end{array} \right ), \end{align*}

with $\delta _k=\frac {k}{L}\gt 0$ . Consequently, we have the following characteristic equation at $E_{*}$

(25) \begin{align} \lambda _k^2-T_k(\xi )\lambda _k+D_k(\xi )=0, \ \textrm {for}\ k\in \mathbb {N}_0=\{0,1,2,\cdot \cdot \cdot \}, \end{align}

where

\begin{align*} \left \{ {\begin{array}{*{20}{l}} T_k(\xi )=-(d_1+d_2)\delta _k^2+f_P+g_Q,\\[5pt] D_k(\xi )=d_1d_2\delta _k^4-\left [f_Pd_2+g_Qd_1+\phi (P^{*})g_P\xi \right ]\delta _k^2+d^2P_{*}Q_{*}, \end{array}} \right . \end{align*}

with

\begin{align*} f_P=-\frac {bc^2P_{*}}{(cP_{*}+eQ_{*})^2}, \ f_Q=dP_{*}-\frac {becP^{*}}{(cP_{*}+eQ_{*})^2}, \ g_P=-dQ_{*}-\frac {bceQ_{*}}{(cP_{*}+eQ_{*})^2}, \ g_Q= -\frac {be^2Q_{*}}{(cP_{*}+eQ_{*})^2}. \end{align*}

It is noticed that $f_P\lt 0$ and $g_Q\lt 0$ . As a consequence, we know that $T_k(\xi )\lt 0$ for any $k\in \mathbb {N}_0$ . This implies that the stability of the unique positive equilibrium $E_{*}$ completely depends on the sign of $D_k(\xi )$ . By direct calculation, we can show that $D_k(\xi )=D_k(\xi _k^S)=0$ when $\xi =\xi _k^S$ , where

\begin{align*} \xi _k^S=\frac {d_1d_2\delta _k^4-(f_Pd_2 +g_Qd_1)\delta _k^2+d^2P_{*}Q_{*}}{\phi (P^{*})g_P\delta _k^2}\lt 0. \end{align*}

We establish the following.

Lemma 3.1. If there is a $k_0\in \mathbb {N}_0\backslash \{0\}$ satisfying

\begin{align*} k_0= \left \{ {\begin{array}{l@{\quad}l} \left [\hat {k}_0\right ]+1, \ &\textrm {if}\ \xi _{[\hat {k}_0]}^S\leq \xi _{[\hat {k}_0]+1}^S,\\[6pt] \left [\hat {k}_0\right ], \ &\textrm {if}\ \xi _{[\hat {k}_0]}^S\gt \xi _{[\hat {k}_0]+1}^S, \end{array}} \right . \end{align*}

where $\hat {k}_0=L\sqrt {d \sqrt {\frac {P_{*}Q_{*}}{d_1d_2}}}$ . Then $\xi _k^S$ has its maximum $\xi _{k_0}^S$ at $k=k_0$ , where $[\cdot ]$ is the integer function.

Proof. Since

\begin{align*} \xi _k^S=\frac {d_1d_2\delta _k^2}{\phi (P^{*})g_P}+\frac {d^2P_{*}Q_{*}}{\phi (P^{*})g_P\delta _k^2}-\frac {f_Pd_2 +g_Qd_1}{\phi (P^{*})g_P}, \end{align*}

we define

\begin{align*} F(z)=\frac {d_1d_2z}{\phi (P^{*})g_P}+\frac {d^2P_{*}Q_{*}}{\phi (P^{*})g_Pz}-\frac {f_Pd_2 +g_Qd_1}{\phi (P^{*})g_P}. \end{align*}

Taking the derivative of $F(z)$ with respect to $z$ , one has

\begin{align*} F^\prime (z)=\frac {d_1d_2}{\phi (P^{*})g_P}-\frac {d^2P_{*}Q_{*}}{\phi (P^{*})g_Pz^2}. \end{align*}

Let $F^\prime (z)=0$ , then we have $z=z_0=d \sqrt {\frac {P_{*}Q_{*}}{d_1d_2}}$ . As a consequence, $F^\prime (z)\lt 0$ as $z\gt z_0$ and $F^\prime (z)\gt 0$ as $0\lt z\lt z_0$ . Moreover, $\lim _{z \rightarrow +\infty }F(z)=-\infty$ since $g_P\lt 0$ . Therefore, $F(z)$ could achieve its maximum at $z=z_0$ . This is

\begin{align*} \mathrm {max}_{z\gt 0}F(z)=F(z_0)=\frac {2d\sqrt {d_1d_2P_{*}Q_{*}}}{\phi (P^{*})g_P}-\frac {f_Pd_2 +g_Qd_1}{\phi (P^{*})g_P}\lt 0. \end{align*}

Recalling that $z=\delta _k^2=(k/L)^2\gt 0$ and the definition of $F(z)$ , we infer that there is a $k_0$ satisfying

\begin{align*} k_0= \left \{ {\begin{array}{l@{\quad}l} \left [\hat {k}_0\right ]+1, \ &\textrm {if}\ \xi _{[\hat {k}_0]}^S\leq \xi _{[\hat {k}_0]+1}^S,\\[6pt] \left [\hat {k}_0\right ], \ &\textrm {if}\ \xi _{[\hat {k}_0]}^S\gt \xi _{[\hat {k}_0]+1}^S \end{array}} \right . \end{align*}

with $\hat {k}_0=L\sqrt {d \sqrt {\frac {P_{*}Q_{*}}{d_1d_2}}}$ such that $\xi _k^S$ has its maximum at $k=k_0$ . The proof is completed.

1.2. Proof of Theorem

Clearly, $f_P\lt 0,\ g_Q\lt 0,$ and $g_P\lt 0$ are valid conditions. As a consequence, we know that $T_k(\xi )\lt 0$ for any $k\in \mathbb {N}_0$ . This implies that the stability of the unique positive equilibrium $E_{*}$ completely depends on the sign of $D_k(\xi )$ . If $\xi \geq 0$ , it immediately follows that $D_k(\xi )\gt 0$ for all $k\in \mathbb {N}_0$ . This shows that all eigenvalues of the characteristic equation (25) with negative real parts. Therefore, (i) is true. If $\xi =\xi _k^S$ is valid, we can check that $D_k(\xi )=0$ , namely, $0$ is an eigenvalue of the characteristic equation (25). Hence, system (1) admits the steady-state bifurcation as $\xi =\xi _k^S$ . Now for $D_k(\xi )$ , taking its derivation with respect to $\xi$ , one yields $D^{\prime }_k(\xi )=-\phi (P^{*})g_P\delta _k^2\gt 0$ . Therefore, $D_k(\xi )$ is strictly increasing about $\xi \in ({-}\infty, 0)$ . Keeping this in mind, if $\xi _k^S\lt \xi \lt 0$ , we have $0=D_k(\xi _k^S)\lt D_k(\xi )$ . Clearly, $E_{*}$ is locally asymptotically stable as $\xi \gt \xi _k^S$ for any $k\in \mathbb {N}_0$ . However, if $\xi \lt \xi _k^S\lt 0$ is valid, we infer that $D_k(\xi )\lt D_k(\xi _k^S)=0$ . This implies that there is at least one eigenvalue of the characteristic equation (25) with a positive real part. In this case, $E_{*}$ is unstable, and (ii) is valid. For (iii), define the following time-evolution Lyapunov function

\begin{align*} V(t)=&\int _\Omega \left (P(\cdot, t)-P_{*}-P_{*} \ln \frac {P(\cdot, t)}{P_{*}}\right )dx+\int _\Omega \left (Q(\cdot, t)-Q_{*}-Q_{*} \ln \frac {Q(\cdot, t)}{Q_{*}}\right )dx. \end{align*}

Then, one yields

\begin{align*} \dot {V}(t)=&\int _\Omega \left (1-\frac {P_{*}}{P}\right )P_{t}dx+\int _\Omega \left (1-\frac {Q_{*}}{Q}\right )Q_{t}dx \\[3pt]=&\int _\Omega (P-P_{*})\left (\frac {bc}{cP+eQ}+dQ-\alpha \right )dx +\int _\Omega (Q-Q_{*})\left (\frac {be}{cP+eQ}-dP-\beta -h\right )dx \\[3pt]&-d_1P_{*}\int _\Omega \frac {\vert \nabla P\vert ^2}{P^2}dx-d_2Q_{*}\int _\Omega \frac {\vert \nabla Q\vert ^2}{Q^2}dx +\int _\Omega \frac {\xi P_{*}\phi (P)\nabla P\nabla Q}{P^2}dx \\[3pt]=&V_1(t)+V_2(t), \end{align*}

where

\begin{align*} V_1(t)=&\int _\Omega (P-P_{*})\left (\frac {bc}{cP+eQ}+dQ-\alpha \right )dx +\int _\Omega (Q-Q_{*})\left (\frac {be}{cP+eQ}-dP-\beta -h\right )dx \\[3pt] =&-b\int _\Omega \frac {[c(P-P_{*})+e(Q-Q_{*})]^2}{(cP+eQ)(cP_{*}+eQ_{*})}dx \\[3pt] \lt &0, \end{align*}

and

\begin{align*} V_2(t)=&-d_1P_{*}\int _\Omega \frac {\vert \nabla P\vert ^2}{P^2}dx-d_2Q_{*}\int _\Omega \frac {\vert \nabla Q\vert ^2}{Q^2}dx +\int _\Omega \frac {\xi P_{*}\phi (P)\nabla P\nabla Q}{P^2}dx \\[3pt] \leq &-d_1P_{*}\int _\Omega \frac {\vert \nabla P\vert ^2}{P^2}dx-d_2Q_{*}\int _\Omega \frac {\vert \nabla Q\vert ^2}{Q^2}dx +\int _\Omega \frac {c_0\xi P_{*}\nabla P\nabla Q}{P}dx \\[3pt] =&-\int _\Omega XAX^Tdx, \end{align*}

where we define $ X(x,t)=(|\nabla P(x,t)|,|\nabla Q(x,t)|)$ in $\Omega \times (0,\infty )$ and

\begin{align*} A=\left ( \begin{array}{c@{\quad}c} \displaystyle \frac {d_1P_{*}}{P^2} & \displaystyle -\frac {c_0\xi P_{*}}{2P}\\[3pt] \displaystyle -\frac {c_0\xi P_{*}}{2P} & \displaystyle \frac {d_2Q_{*}}{Q^2} \end{array} \right ). \end{align*}

Therefore, if $0\lt \xi ^2\lt \frac {4d_1d_2Q_{*}}{c_0^2P_{*}C_1^2}$ is valid, one obtains

\begin{align*} Trace (A)=\frac {d_1P_{*}}{P^2}+\frac {d_2Q_{*}}{Q^2}\gt 0,\ \ Det (A)=\frac {P_{*}}{P^2}\left (\frac {d_1d_2Q_{*}}{Q^2}-\frac {c_0^2\xi ^2 P_{*}}{4}\right )\gt 0. \end{align*}

Hence, $A$ is a positive definite matrix which implies that $V_2(t)=-\int _\Omega XAX^Tdx\lt 0$ is true. Thereby, $\dot {V}(t)=V_1(t)+V_2(t)\lt 0$ is valid, namely, the unique positive equilibrium $E_{*}$ is globally asymptotically stable. This finishes the proof.

3.2 Bifurcating solution: nonconstant steady state

3.2.1 Existence

In this subsection, we explore the existence and stability of the nonconstant steady states around the steady state bifurcation onset $\xi =\xi _k^S$ for $k\in \mathbb {N}_0\backslash \{0\}$ . Define two Hilbert spaces: $\mathbf {X}=\{u\in H^2(0,L\pi )|u^{\prime }(0)=u^{\prime }(L\pi )=0\}\ \textrm {and}\ \mathbf {Y}=L^2(0,L\pi )$ . Rewrite system (2) as follows.

(26) \begin{align} \left \{ {\begin{array}{l@{\quad}l} \displaystyle 0=(d_1P^{\prime }-\xi \phi (P) Q^{\prime })^{\prime }+P\left (\frac {bc}{cP+eQ}+dQ-\alpha \right ),\ \ &x\in \Omega .\\[3pt] \displaystyle 0=d_2 Q^{\prime \prime }+Q\left (\frac {be}{cP+eQ}-dP-\beta \right )-hQ,\ \ &x\in \Omega .\\[3pt] P^{\prime }= Q^{\prime }=0,\ \ &x\in \partial \Omega . \end{array}} \right . \end{align}

Recalling the operator $\mathcal {F}(P,Q,\xi )$ in (4), then system (26) is equivalent to $\mathcal {F}(P,Q,\xi )=0$ and $\mathcal {F}(P,Q,\xi ): \mathbf {X}\times \mathbf {X}\times \mathbb {R}\longrightarrow \mathbf {Y}\times \mathbf {Y}$ is analytic for $(P,Q,\xi )\in \mathbf {X}\times \mathbf {X}\times \mathbb {R}$ . Now, at the onset $\xi =\xi _k^S$ , we can confirm that $\mathcal {N}(D_{(P,Q)}\mathcal {F}(P_{*},Q_{*},\xi _k^S))\neq \{0\}$ , where $\mathcal {N}$ is the null space and $D_{(P,Q)}\mathcal {F}(P_{*},Q_{*},\xi )$ has been appeared in (6). Benefiting from (6), we infer that the null space $\mathcal {N}$ consists of solutions to the problem

(27) \begin{align} \left \{ {\begin{array}{*{20}{l}} \displaystyle 0=d_1P^{\prime \prime }-\xi (\phi (P_{*}) Q^{\prime })^{\prime } -\frac {bc^2P_{*}}{(cP_{*}+eQ_{*})^2}P +\left [dP_{*}-\frac {bceP_{*}}{(cP_{*}+eQ_{*})^2}\right ]Q,\\[6pt] \displaystyle 0=d_2 Q^{\prime \prime }-\left [dQ_{*}+\frac {bceQ_{*}}{(cP_{*}+eQ_{*})^2}\right ]P -\frac {be^2Q_{*}}{(cP_{*}+eQ_{*})^2}Q,\\[6pt] \displaystyle P^{\prime }= Q^{\prime }=0. \end{array}} \right . \end{align}

Let

\begin{align*} P(x)=\sum _{k=0}^\infty a^{\prime }_k\mathrm {cos}\frac {k x}{L},\ Q(x)=\sum _{k=0}^\infty b^{\prime }_k\mathrm {cos}\frac {k x}{L}. \end{align*}

Then, putting them into (27), we get

\begin{align*} \left ( \begin{array}{c@{\quad}c} \displaystyle -\frac {bc^2P_{*}}{(cP_{*}+eQ_{*})^2}-d_1\delta _k^2 & \displaystyle dP_{*}-\frac {becP^{*}}{(cP_{*}+eQ_{*})^2}+\xi \phi (P^{*})\delta _k^2 \\[3pt] \displaystyle -dQ_{*}-\frac {bceQ_{*}}{(cP_{*}+eQ_{*})^2} & \displaystyle -\frac {be^2Q_{*}}{(cP_{*}+eQ_{*})^2}-d_2\delta _k^2 \end{array} \right )\left ( \begin{array}{c} a^{\prime }_k \\[3pt] b^{\prime }_k \end{array} \right )=\left ( \begin{array}{c} 0 \\[3pt] 0 \end{array} \right ). \end{align*}

Consequently, if $\xi =\xi _k^S$ , we have $ \mathcal {N}(D_{(P,Q)}\mathcal {F}(P_{*},Q_{*},\xi _k^S))=\textrm {span}\{\widehat {P}_{k},\widehat {Q}_{k}\},$ where $\widehat {P}_{k}$ and $\widehat {Q}_{k}$ can be found in (9). Moreover, utilizing Theorem 3.3 of [33] or Lemma2.3 of [36], we know that $D_{(P,Q)}\mathcal {F}(P_{*},Q_{*},\xi )$ is a Fredholm operator with index $0$ and $codim\mathcal {R}(D_{(P,Q)}\mathcal {F}(P_{*},Q_{*},\xi ))=1$ is true, where $\mathcal {R}$ is the range of the operator.

Now we can show the validity of Theorem1.3.

1.3. Proof of Theorem

Owing to the Crandall–Rabinowitz bifurcation theory [9], we only need to prove the following transversality condition

(28) \begin{align} \frac {d}{d\xi }(D_{(P,Q)}\mathcal {F}(P_{*},Q_{*},\xi ))(\widehat {P}_{k},\widehat {Q}_{k})|_{\xi =\xi _k^S} \notin \mathcal {R}(D_{(P,Q)}\mathcal {F}(P_{*},Q_{*},\xi )) \end{align}

is true, where $\mathcal {R}$ denotes the range of the operator. Now, we assume that (28) fails, then from (6), we can set

(29) \begin{align} \left \{ {\begin{array}{*{20}{l}} \displaystyle d_1P^{\prime \prime }-\xi _k^S(\phi (P_{*}) Q^{\prime })^{\prime } -\frac {bc^2P_{*}}{(cP_{*}+eQ_{*})^2}P +\left [dP_{*}-\frac {bceP_{*}}{(cP_{*}+eQ_{*})^2}\right ]Q=\phi (P_{*})\delta _k^2\mathrm {cos}\frac {k x}{L},\\[6pt] \displaystyle d_2 Q^{\prime \prime }-\left [dQ_{*}+\frac {bceQ_{*}}{(cP_{*}+eQ_{*})^2}\right ]P -\frac {be^2Q_{*}}{(cP_{*}+eQ_{*})^2}Q=0,\\[3pt] P^{\prime }= Q^{\prime }=0. \end{array}} \right . \end{align}

Then, multiplying (29) by $\mathrm {cos}\frac {k x}{L}$ and integrating it over $(0,L\pi )$ , one obtains

(30) \begin{align} \left ( \begin{array}{c@{\quad}c} \displaystyle -\frac {bc^2P_{*}}{(cP_{*}+eQ_{*})^2}-d_1\delta _k^2 & \displaystyle dP_{*}-\frac {becP^{*}}{(cP_{*}+eQ_{*})^2}+\xi _k^S\phi (P^{*})\delta _k^2 \\[3pt] \displaystyle -dQ_{*}-\frac {bceQ_{*}}{(cP_{*}+eQ_{*})^2} & \displaystyle -\frac {be^2Q_{*}}{(cP_{*}+eQ_{*})^2}-d_2\delta _k^2 \end{array} \right )\left ( \begin{array}{c} \displaystyle {\int _0^{L\pi } P\mathrm {cos}\frac {k x}{L}dx}\\[3pt] \displaystyle {\int _0^{L\pi } Q\mathrm {cos}\frac {k x}{L}dx} \end{array} \right )=\left ( \begin{array}{c} \displaystyle \frac {\pi \delta _k^2\phi (P_{*})L}{2} \\[3pt] 0 \end{array} \right ). \end{align}

Because there is the steady-state bifurcation when $\xi =\xi _k^S$ , we obtain

\begin{align*} \left | \begin{array}{c@{\quad}c} \displaystyle -\frac {bc^2P_{*}}{(cP_{*}+eQ_{*})^2}-d_1\delta _k^2 & \displaystyle dP_{*}-\frac {becP^{*}}{(cP_{*}+eQ_{*})^2}+\xi _k^S\phi (P^{*})\delta _k^2 \\[3pt] \displaystyle -dQ_{*}-\frac {bceQ_{*}}{(cP_{*}+eQ_{*})^2} & \displaystyle -\frac {be^2Q_{*}}{(cP_{*}+eQ_{*})^2}-d_2\delta _k^2 \end{array} \right |=0. \end{align*}

This leads to a contradiction from (30), and thereby (28) is valid. We end the proof.

3.2.2 Stability

In this subsection, we want to ensure the stability of the bifurcating solution $(P_{k}(\varepsilon, x),Q_{k}(\varepsilon, x))$ in Theorem1.3. To this end, from (26), we know that the bifurcating solution $(P_{k}(\varepsilon, x),Q_{k}(\varepsilon, x))$ admits

(31) \begin{align} \left \{ {\begin{array}{*{20}{l}} \displaystyle 0=(d_1P^{\prime }_k(\varepsilon, x)-\xi _k^S(\varepsilon )\phi (P_k(\varepsilon, x)) Q^{\prime }_k(\varepsilon, x))^{\prime }+P_k(\varepsilon, x)\left (\frac {bc}{cP_k(\varepsilon, x)+eQ_k(\varepsilon, x)}+dQ_k(\varepsilon, x)-\alpha \right ),\\[6pt] \displaystyle 0= d_2 Q_k^{\prime \prime }(\varepsilon, x)+Q_k(\varepsilon, x)\left (\frac {be}{cP_k(\varepsilon, x)+eQ_k(\varepsilon, x)}-dP_k(\varepsilon, x)-\beta -h\right ),\\[6pt] P_k^{\prime }= Q_k^{\prime }=0. \end{array}} \right . \end{align}

In the sequel, let us expand the critical threshold $\xi _k^S(\varepsilon )$ and bifurcating solution $(P_{k}(\varepsilon, x),Q_{k}(\varepsilon, x))$ as below:

(32) \begin{align} \left \{ {\begin{array}{*{20}{l}} \xi _k^S(\varepsilon )=\xi _k^S+\varepsilon \xi _1+\varepsilon ^2\xi _2+\cdot \cdot \cdot, \\[6pt] P_{k}(\varepsilon, x)=P_{*}+\varepsilon \mathrm {cos}\frac {k x}{L}+\varepsilon ^2P_1(x)+\varepsilon ^3P_2(x)+\cdot \cdot \cdot, \\[6pt] Q_{k}(\varepsilon, x)=Q_{*}+\varepsilon \alpha _k\mathrm {cos}\frac {k x}{L}+\varepsilon ^2Q_1(x)+\varepsilon ^3Q_2(x)+\cdot \cdot \cdot, \end{array}} \right . \end{align}

where $\xi _1,\xi _2$ will be computed later and $(P_j,Q_j)\in \mathcal {K}$ for $j=1,2$ , where $\mathcal {K}$ has been defined in (8). For the density function $\phi (P_{k}(\varepsilon, x))$ , we set

\begin{align*} \phi (P_{k}(\varepsilon, x))=\phi (P_{*})+\phi _{P_{k}(\varepsilon, x)}(P_{*})P_{k}(\varepsilon, x) +\frac {1}{2}\phi _{P_{k}(\varepsilon, x)P_{k}(\varepsilon, x)}(P_{*})P_{k}^2(\varepsilon, x)+\cdot \cdot \cdot, \end{align*}

Using the second perturbation of (32), we get

(33) \begin{align} \phi (P_{k}(\varepsilon, x))=\phi (P_{*})+\varepsilon \phi _{P_{k}}(P_{*})\mathrm {cos}\frac {k x}{L}+\varepsilon ^2\left (\phi _{P_{k}}(P_{*})P_1(x)+ \frac {1}{2}\phi _{P_{k}P_{k}}(P_{*})\mathrm {cos}^2\frac {k x}{L}\right )+\cdot \cdot \cdot . \end{align}

Then, submitting (32)–(33) into (31), we obtain

(34) \begin{align} \left \{ {\begin{array}{*{20}{l}} 0=(d_1 P^{\prime }_{k}-\xi _k^S(\varepsilon )\phi (P_{k}) Q_k^{\prime })^{\prime }+f_PP_{*}+f_QQ_{*}+\mathcal {R}_0+\varepsilon \mathcal {R}_1(x)+\varepsilon ^2 \mathcal {R}_2(x)+\varepsilon ^3 \mathcal {R}_3(x)+\cdot \cdot \cdot, \\[6pt] 0=d_2 Q_k^{\prime \prime }+g_PP_{*}+g_QQ_{*}+\mathcal {V}_0+\varepsilon \mathcal {V}_1(x)+\varepsilon ^2 \mathcal {V}_2(x)+\varepsilon ^3 \mathcal {V}_3(x)+\cdot \cdot \cdot, \\[6pt] P^{\prime }_k= Q^{\prime }_k=0, \end{array}} \right . \end{align}

where

\begin{align*} &(d_1 P^{\prime }_{k}(\varepsilon, x)-\xi _k^S(\varepsilon )\phi (P_{k}(\varepsilon, x)) Q^{\prime }_k(\varepsilon, x))^{\prime }\\[3pt] =&d_1 P^{\prime \prime }_{k}(\varepsilon, x)-\xi _k^S(\varepsilon )(\phi ^{\prime }(P_{k}(\varepsilon, x))Q^{\prime }_k(\varepsilon, x)+\phi (P_{k}(\varepsilon, x)) Q^{\prime \prime }_k(\varepsilon, x)) \\[3pt]=&\varepsilon \left [\delta _k^2\xi _k^S\phi (P_{*})\alpha _k\mathrm {cos}\frac {kx}{L}-\delta _k^2d_1\mathrm {cos}\frac {kx}{L}\right ] +\varepsilon ^2\left [d_1P^{\prime \prime }_1(x)+\delta _k^2\xi _1\phi (P_{*})\alpha _k\mathrm {cos}\frac {kx}{L} \right .\\[3pt] &\left .+\delta _k^2\xi _k^S\phi _{P_{1k}}(P_{*})\alpha _k\mathrm {cos}\frac {2kx}{L}-\xi _k^S\phi (P_{*})Q^{\prime \prime }_1(x)\right ] +\varepsilon ^3\left [d_1P^{\prime \prime }_2(x)+\delta _k^2\xi _1\phi _{P_{k}}(P_{*})\alpha _k\mathrm {cos}\frac {2kx}{L}\right .\\[3pt] &\left .-\xi _1\phi (P_{*})Q^{\prime \prime }_1(x)-\xi _k^S\phi _{P_{k}}(P_{*})\mathrm {cos}\frac {kx}{L}Q^{\prime \prime }_1(x)-\xi _k^S\phi (P_{*})Q^{\prime \prime }_2(x) +\xi _k^S\phi _{P_{k}}(P_{*})\delta _k\sin \frac {kx}{L}Q^{\prime }_1(x)\right .\\[3pt] &\left .+\delta _k^2\xi _k^S\alpha _k\mathrm {cos}\frac {kx}{L}\left (\phi _{P_{k}}(P_{*})P_1(x)+\frac {1}{2}\phi _{P_{k}P_{k}}(P_{*})\mathrm {cos}^2\frac {k x}{L}\right )+\delta _k\xi _k^S\alpha _k\sin \frac {kx}{L}\right .\\[3pt] &\left .\times \left (\phi _{P_{k}}(P_{*})P^{\prime }_1(x)-\frac {1}{2}\phi _{P_{k}P_{k}}(P_{*})\delta _k\sin \frac {2k x}{L}\right ) +\delta _k^2\xi _2\phi (P_{*})\alpha _k\mathrm {cos}\frac {kx}{L}\right ], \end{align*}

and $\mathcal {R}_0,\mathcal {R}_j(x),\mathcal {V}_0,\mathcal {V}_j(x)$ for $j=1,2,3$ can be found in Appendix A and B, respectively.

1.4. Proof of Theorem

To obtain the desired results, we should first determine the values of $\xi _1$ and $\xi _2$ , respectively. From the perturbation equation (34), we can get $\mathcal {O}(\varepsilon ^2)$ term as below.

(35) \begin{align} \left \{ {\begin{array}{*{20}{l}} 0=d_1P^{\prime \prime }_1(x)+\Theta (x)+ \mathcal {R}_2(x),\\[3pt] 0=d_2 Q_1^{\prime \prime }(x)+ \mathcal {V}_2(x),\\[3pt] P^{\prime }_1(x)=Q_1^{\prime }(x)=0, \end{array}} \right . \end{align}

where $ \Theta (x)=\delta _k^2\xi _1\phi (P_{*})\alpha _k\mathrm {cos}\frac {kx}{L} +\delta _k^2\xi _k^S\phi _{P_{k}}(P_{*})\alpha _k\mathrm {cos}\frac {2kx}{L}-\xi _k^S\phi (P_{*})Q^{\prime \prime }_1(x).$ Multiplying (35) by $\mathrm {cos}\frac {kx}{L}$ and integrating it over $(0,L\pi )$ , we have

(36) \begin{align} \frac {\delta _k^2L\pi }{2}\xi _1\phi (P_{*})\alpha _k=\left (d_1\delta _k^2-\mathcal {R}_{21}\right )\int _0^{L\pi } P_1(x)\mathrm {cos}\frac {k x}{L}dx -\left (\xi _k^S\phi (P_{*})\delta _k^2+\mathcal {R}_{22}\right )\int _0^{L\pi } Q_1(x)\mathrm {cos}\frac {k x}{L}dx, \end{align}

and

(37) \begin{align} 0=\mathcal {V}_{21}\int _0^{L\pi } P_1(x)\mathrm {cos}\frac {k x}{L}dx +\left (\mathcal {V}_{22}-d_2\delta _k^2\right )\int _0^{L\pi }Q_1(x)\mathrm {cos}\frac {k x}{L}dx, \end{align}

where

\begin{align*} \mathcal {R}_{21}=&f_P+f_{PP}P_{*}+f_{PQ}Q_{*}+f_{PPQ}Q_{*}P_{*}+\frac {f_{PQQ}}{2}Q_{*}^2+\frac {f_{PPP}}{2}P_{*}^2,\\[3pt] \mathcal {R}_{22}=&f_Q+f_{PQ}P_{*}+f_{QQ}Q_{*}+\frac {f_{QQQ}}{2}Q_{*}^2+\frac {f_{PPQ}}{2}P_{*}^2+f_{PQQ}P_{*}Q_{*},\\[3pt] \mathcal {V}_{21}=&g_P+g_{PP}P_{*}+g_{PQ}Q_{*}+g_{PPQ}Q_{*}P_{*}+\frac {g_{PQQ}}{2}Q_{*}^2+\frac {g_{PPP}}{2}P_{*}^2,\\[3pt] \mathcal {V}_{22}=&g_Q+g_{PQ}P_{*}+g_{QQ}Q_{*}+\frac {g_{QQQ}}{2}Q_{*}^2+\frac {g_{PPQ}}{2}P_{*}^2+g_{PQQ}Q_{*}P_{*}. \end{align*}

Moreover, in light of (8) and $(P_1,Q_1)\in \mathcal {K}$ , we have

(38) \begin{align} 0=\int _0^{L\pi } P_1(x)\mathrm {cos}\frac {k x}{L}dx+\alpha _k\int _0^{L\pi } Q_1(x)\mathrm {cos}\frac {k x}{L}dx, \end{align}

where $\alpha _k$ can be found in (9). Consequently, by using (37) and (38), one obtains

\begin{align*} \left ( \begin{array}{c@{\quad}c} \mathcal {V}_{21}&\mathcal {V}_{22}-d_2\delta _k^2\\[3pt] 1&\alpha _k \end{array} \right )\left ( \begin{array}{c} {\int _0^{L\pi } P_1(x)\mathrm {cos}\frac {k x}{L}dx}\\[3pt] {\int _0^{L\pi } Q_1(x)\mathrm {cos}\frac {k x}{L}dx} \end{array} \right )=\left ( \begin{array}{c} 0\\[3pt] 0 \end{array} \right ). \end{align*}

It is noticed that

\begin{align*} \left | \begin{array}{c@{\quad}c} \mathcal {V}_{21}&\mathcal {V}_{22}-d_2\delta _k^2\\[3pt] 1&\alpha _k \end{array} \right |=\alpha _k\mathcal {V}_{21}-\mathcal {V}_{22}+d_2\delta _k^2\neq 0. \end{align*}

This gives

(39) \begin{align} {\int _0^{L\pi } P_1(x)\mathrm {cos}\frac {k x}{L}dx}= {\int _0^{L\pi } Q_1(x)\mathrm {cos}\frac {k x}{L}dx}=0. \end{align}

Putting (39) into (36), we infer that $\xi _1=0$ .

Our next task is to determine $\xi _2$ in the first perturbation equation of (32). To this end, we investigate the $\mathcal {O}(\varepsilon ^3)$ term of (34). This is

(40) \begin{align} \left \{ {\begin{array}{*{20}{l}} 0=d_1P^{\prime \prime }_2(x)+\Phi (x)+ \mathcal {R}_3(x),\\[3pt] 0=d_2 Q_2^{\prime \prime }(x)+ \mathcal {V}_3(x),\\[3pt] P^{\prime }_2(x)=Q_2^{\prime }(x)=0, \end{array}} \right . \end{align}

where

\begin{align*} \Phi (x)=&\delta _k^2\xi _1\phi _{P_{k}}(P_{*})\alpha _k\mathrm {cos}\frac {2kx}{L} -\xi _1\phi (P_{*})Q^{\prime \prime }_1(x)-\xi _k^S\phi _{P_{k}}(P_{*})\mathrm {cos}\frac {kx}{L}Q^{\prime \prime }_1(x)-\xi _k^S\phi (P_{*})Q^{\prime \prime }_2(x) \\[3pt]&+\xi _k^S\phi _{P_{k}}(P_{*})\delta _k\sin \frac {kx}{L}Q^{\prime }_1(x) +\delta _k^2\xi _k^S\alpha _k\mathrm {cos}\frac {kx}{L}\left (\phi _{P_{k}}(P_{*})P_1(x)+\frac {1}{2}\phi _{P_{k}P_{k}}(P_{*})\mathrm {cos}^2\frac {k x}{L}\right )\\[3pt]&+\delta _k\xi _k^S\alpha _k\sin \frac {kx}{L} \left (\phi _{P_{k}}(P_{*})P^{\prime }_1(x)-\frac {1}{2}\phi _{P_{k}P_{k}}(P_{*})\delta _k\sin \frac {2k x}{L}\right ) +\delta _k^2\xi _2\phi (P_{*})\alpha _k\mathrm {cos}\frac {kx}{L}. \end{align*}

Let us multiply (40) by $\mathrm {cos}\frac {kx}{L}$ and integrating over $(0,L\pi )$ , one yields

(41) \begin{align} -\frac {\delta _k^2\pi L\phi (P*)\alpha _k}{2}\xi _2=&\frac {\delta _k^2\xi _k^S}{2}\phi _{P_{k}}(P^{*})\alpha _k\int _0^{L\pi } P_1(x)\left (1-\mathrm {cos}\frac {2k x}{L}\right )dx+\frac {\mathcal {R}_{31}}{2}\int _0^{L\pi } P_1(x)\left (1+\mathrm {cos}\frac {2k x}{L}\right )dx \\[3pt] \nonumber &+\delta _k^2\xi _k^S\phi _{P_{k}}(P^{*})\int _0^{L\pi } Q_1(x)\mathrm {cos}\frac {2k x}{L}dx +\frac {\mathcal {R}_{32}}{2}\int _0^{L\pi } Q_1(x)\left (1+\mathrm {cos}\frac {2k x}{L}\right )dx \\[3pt]\nonumber &+\left (\mathcal {R}_{33}-d_1\delta _k^2\right )\int _0^{L\pi }P_2(x)\mathrm {cos}\frac {k x}{L}dx +\left (\delta _k^2\xi _k^S\phi _{P_{k}}(P^{*})+\mathcal {R}_{34}\right )\int _0^{L\pi } Q_2(x)\mathrm {cos}\frac {k x}{L}dx \\[3pt]\nonumber &+\frac {\delta _k^2L\pi }{16}\phi _{P_{k}P_{k}}(P*)\xi _k^S\alpha _k+\frac {3L\pi }{8}\mathcal {R}_{35}, \end{align}

and

(42) \begin{align} 0=&\left (\mathcal {V}_{34}-d_2\delta _k^2\right )\int _0^{L\pi }Q_2(x)\mathrm {cos}\frac {k x}{L}dx +\frac {\mathcal {V}_{31}}{2}\int _0^{L\pi }P_1(x)\left (1+\mathrm {cos}\frac {2k x}{L}\right )dx\\[3pt] \nonumber & +\frac {V_{32}}{2}\int _0^{L\pi }Q_1(x)\left (1+\mathrm {cos}\frac {2k x}{L}\right )dx +\mathcal {V}_{33}\int _0^{L\pi }P_2(x)\mathrm {cos}\frac {k x}{L}dx +\frac {3L\pi }{8}\mathcal {V}_{35}, \end{align}

where

\begin{align*} \mathcal {R}_{31}=&f_{PP}+f_{PQ}\alpha _k+f_{PPQ}(Q_{*}+\alpha _kP_{*})+f_{PQQ}Q_{*}\alpha _k +\frac {5 f_{PPP}}{6}P_{*},\\[3pt]\mathcal {R}_{32}=&f_{PQ}+f_{QQ}\alpha _k+\frac {5 \alpha _k f_{QQQ}}{6}Q_{*}+f_{PPQ}P_{*}+f_{PQQ}(\alpha _kP_{*}+Q_{*}), \end{align*}

\begin{align*} \mathcal {R}_{33}=&f_{PP}P_{*}+f_{PQ}Q_{*}+f_{PPQ}P_{*}Q_{*}+\frac {f_{PQQ}}{2}Q_{*}^2+\frac {f_{PPP}}{2}P_{*}^2+f_P,\\[3pt]\mathcal {R}_{34}=&f_{PQ}P_{*}+f_{QQ}Q_{*}+\frac {f_{QQQ}}{2}Q_{*}^2+\frac {f_{PPQ}}{2}P_{*}^2+f_{PQQ}P_{*}Q_{*}+f_Q,\\[3pt]\mathcal {R}_{35}=&\frac {\alpha _k^3 f_{QQQ}}{3!}+\frac {\alpha _kf_{PPQ}}{2}+\frac {\alpha _k^2f_{PQQ}}{2}+\frac {f_{PPP}}{3!}, \end{align*}

\begin{align*} \mathcal {V}_{31}=&g_{PP}+g_{PQ}\alpha _k+g_{PPQ}(Q_{*}+P_{*}\alpha _k)+g_{PQQ}P_{*}\alpha _k +\frac {5 g_{PPP}}{6}P^{*},\\[3pt] \mathcal {V}_{32}=&g_{PQ}+g_{QQ}\alpha _k+\frac {5\alpha _k g_{QQQ}}{6}Q_{*}+g_{PPQ}P_{*}+g_{PQQ}(\alpha _kP_{*}+Q_{*}),\\[3pt] \mathcal {V}_{33}=&g_{PP}P_{*}+g_{PQ}Q_{*}+g_{PPQ}P_{*}Q_{*}+\frac {g_{PQQ}}{2}Q_{*}^2+\frac {g_{PPP}}{2}P_{*}^2+g_P,\\[3pt] \mathcal {V}_{34}=&g_{PQ}P_{*}+g_{QQ}Q_{*}+\frac {g_{QQQ}}{2}Q_{*}^2+\frac {g_{PPQ}}{2}P_{*}^2 +g_{PQQ}P_{*}Q_{*}+g_Q,\\[3pt] \mathcal {V}_{35}=&\frac {g_{QQQ}\alpha _k^3}{3!}+\frac {\alpha _k g_{PPQ}}{2}+\frac {\alpha _k^2 g_{PQQ}}{2} +\frac {g_{PPP}}{3!}. \end{align*}

On the other hand, by using (8) and $(P_2,Q_2)\in \mathcal {K}$ , we get

(43) \begin{align} 0=\int _0^{L\pi } P_2(x)\mathrm {cos}\frac {k x}{L}dx+\alpha _k\int _0^{L\pi } Q_2(x)\mathrm {cos}\frac {k x}{L}dx. \end{align}

Obviously, to get the expression of $\xi _2$ in (41), we have to compute

\begin{align*} {\int _0^{L\pi } P_2(x)\mathrm {cos}\frac {k x}{L}dx},{\int _0^{L\pi } Q_2(x)\mathrm {cos}\frac {k x}{L}dx},{\int _0^{L\pi } P_1(x)dx},\\[3pt] {\int _0^{L\pi } Q_1(x)dx},{\int _0^{L\pi } P_1(x)\mathrm {cos}\frac {2k x}{L}dx},{\int _0^{L\pi } Q_1(x)\mathrm {cos}\frac {2k x}{L}dx}. \end{align*}

To this end, we utilize three steps to finish this task.

Step 1: Compute $\int _0^{L\pi } P_2(x)\mathrm {cos}\frac {k x}{L}dx$ and $\int _0^{L\pi } Q_2(x)\mathrm {cos}\frac {k x}{L}dx$ .

In light of (42)–(43), one gets

(44) \begin{align} \left ( \begin{array}{c@{\quad}c} \mathcal {V}_{33}&\mathcal {V}_{34}-d_2\delta _k^2\\[3pt] 1&\alpha _k \end{array} \right )\left ( \begin{array}{c} {\int _0^{L\pi } P_2(x)\mathrm {cos}\frac {k x}{L}dx}\\[3pt] {\int _0^{L\pi } Q_2(x)\mathrm {cos}\frac {k x}{L}dx} \end{array} \right )=\left ( \begin{array}{c} w\\[3pt] 0 \end{array} \right ), \end{align}

where

\begin{align*} w=-\frac {\mathcal {V}_{31}}{2}\int _0^{L\pi }P_1(x)\left (1+\mathrm {cos}\frac {2k x}{L}\right )dx -\frac {\mathcal {V}_{32}}{2}\int _0^{L\pi }Q_1(x)\left (1+\mathrm {cos}\frac {2k x}{L}\right )dx-\frac {3L\pi }{8}\mathcal {V}_{35}. \end{align*}

By solving (44), we have

(45) \begin{align} {\int _0^{L\pi } P_2(x)\mathrm {cos}\frac {k x}{L}dx}=\frac {\left | \begin{array}{c@{\quad}c} w&\mathcal {V}_{34}-d_2\delta _k^2\\[3pt] 0&\alpha _k \end{array} \right |}{\left | \begin{array}{c@{\quad}c} \mathcal {V}_{33}&\mathcal {V}_{34}-d_2\delta _k^2\\[3pt] 1&\alpha _k \end{array} \right |},\ {\int _0^{L\pi }Q_2(x)\mathrm {cos}\frac {k x}{L}dx}=\frac {\left | \begin{array}{c@{\quad}c} \mathcal {V}_{33}&w\\[3pt] 1&0 \end{array} \right |}{\left | \begin{array}{c@{\quad}c} \mathcal {V}_{33}&\mathcal {V}_{34}-d_2\delta _k^2\\[3pt] 1&\alpha _k \end{array} \right |}. \end{align}

Step 2: Compute $\int _0^{L\pi } P_1(x)dx$ and $\int _0^{L\pi } Q_1(x)dx$ .

Integrating (35) over $(0,L\pi )$ , we have

(46) \begin{align} \left \{ {\begin{array}{*{20}{l}} 0=d_1{\int _0^{L\pi }\Theta (x)dx+ \int _0^{L\pi }\mathcal {R}_2(x)dx},\\[3pt] 0=d_2 { \int _0^{L\pi }\mathcal {V}_2(x)dx},\\[3pt] P^{\prime }_1(x)=Q_1^{\prime }(x)=0, \end{array}} \right . \end{align}

where we employ

\begin{align*} \int _0^{L\pi }P^{\prime \prime }_1(x)dx= \int _0^{L\pi }Q_1^{\prime \prime }(x)dx=0. \end{align*}

In addition, we can calculate

(47) \begin{align} \int _0^{L\pi }\Theta (x)dx= \delta _k^2\xi _k^S\phi _{P_{k}}(P^{*})\alpha _k\int _0^{L\pi }\mathrm {cos}\frac {2kx}{L}dx-\xi _k^S\phi (P^{*})\int _0^{L\pi }Q^{\prime \prime }_1(x)dx=0, \end{align}

(48) \begin{align} \int _0^{L\pi }\mathcal {R}_2(x)dx=\mathcal {R}_{21}\int _0^{L\pi }P_1(x)dx+\mathcal {R}_{22}\int _0^{L\pi }Q_1(x)dx+\frac {L\pi }{2}\widetilde {\mathcal {R}}_2, \\[6pt] \nonumber \end{align}

and

(49) \begin{align} \int _0^{L\pi }\mathcal {V}_2(x)dx=\mathcal {V}_{21}\int _0^{L\pi }P_1(x)dx +\mathcal {V}_{22}\int _0^{L\pi }Q_1(x)dx+\frac {L\pi }{2}\widetilde {\mathcal {V}}_2, \end{align}

where

\begin{align*} \widetilde {\mathcal {R}}_2=\frac {f_{PP}}{2}+f_{PQ}\alpha _k+\frac {f_{QQ}\alpha _k^2}{2}+\frac {f_{QQQ}\alpha _k^2}{2}Q_{*} +\frac {f_{PPQ}}{2}(2\alpha _kP_{*}+1) +\frac {f_{PQQ}}{2}(2\alpha _kQ_{*}+P_{*}\alpha _k^2) +\frac {f_{PPP}}{2}P_{*},\\[3pt] \widetilde {\mathcal {V}}_2=\frac {g_{PP}}{2}+g_{PQ}\alpha _k+\frac {g_{QQ}\alpha _k^2}{2}+\frac {g_{QQQ}\alpha _k^2}{2}Q_{*} +\frac {g_{PPQ}}{2}(2\alpha _kP_{*}+1) +\frac {g_{PQQ}}{2}(2\alpha _kQ_{*}+P_{*}\alpha _k^2)+\frac {g_{PPP}}{2}P_{*}. \end{align*}

Consequently, putting (47)–(49) into (46), we can get

(50) \begin{align} \left \{ {\begin{array}{*{20}{l}} 0={\mathcal {R}_{21}\int _0^{L\pi }P_1(x)dx+\mathcal {R}_{22}\int _0^{L\pi }Q_1(x)dx+\frac {L\pi }{2}\widetilde {\mathcal {R}}_2},\\[3pt] 0={\mathcal {V}_{21}\int _0^{L\pi }P_1(x)dx +\mathcal {V}_{22}\int _0^{L\pi }Q_1(x)dx+\frac {L\pi }{2}\widetilde {\mathcal {V}}_2},\\[3pt] P^{\prime }_1(x)=Q_1^{\prime }(x)=0. \end{array}} \right . \end{align}

This is

(51) \begin{align} \left ( \begin{array}{c@{\quad}c} \mathcal {R}_{21}&\mathcal {R}_{22}\\[3pt] \mathcal {V}_{21}&\mathcal {V}_{22} \end{array} \right )\left ( \begin{array}{c} {\int _0^{L\pi } P_1(x)dx}\\[3pt] {\int _0^{L\pi } Q_1(x)dx} \end{array} \right )=\left ( \begin{array}{c} -\frac {L\pi }{2}\widetilde {\mathcal {R}}_2\\[3pt] -\frac {L\pi }{2}\widetilde {\mathcal {V}}_2 \end{array} \right ). \end{align}

By solving (51), we obtain

(52) \begin{align} {\int _0^{L\pi } P_1(x)dx}=\frac {\left | \begin{array}{c@{\quad}c} -\frac {L\pi }{2}\widetilde {\mathcal {R}}_2&\mathcal {R}_{22}\\[3pt] -\frac {L\pi }{2}\widetilde {\mathcal {V}}_2&\mathcal {V}_{22} \end{array} \right |}{\left | \begin{array}{c@{\quad}c} \mathcal {R}_{21}&\mathcal {R}_{22}\\[3pt] \mathcal {V}_{21}&\mathcal {V}_{22} \end{array} \right |},\ {\int _0^{L\pi } Q_1(x)dx}=\frac {\left | \begin{array}{c@{\quad}c} \mathcal {R}_{21}&-\frac {L\pi }{2}\widetilde {\mathcal {R}}_2\\[3pt] \mathcal {V}_{21}&-\frac {L\pi }{2}\widetilde {\mathcal {V}}_2 \end{array} \right |}{\left | \begin{array}{c@{\quad}c} \mathcal {R}_{21}&\mathcal {R}_{22}\\[3pt] \mathcal {V}_{21}&\mathcal {V}_{22} \end{array} \right |}. \end{align}

Step 3: Compute $\int _0^{L\pi } P_1(x)\mathrm {cos}\frac {2k x}{L}dx$ and $\int _0^{L\pi } Q_1(x)\mathrm {cos}\frac {2k x}{L}dx$ .

Multiplying (35) by $\mathrm {cos}\frac {2k x}{L}$ , one yields

(53) \begin{align} \left \{ {\begin{array}{*{20}{l}} 0=d_1{\int _0^{L\pi }P^{\prime \prime }_1(x)\mathrm {cos}\frac {2k x}{L}dx +\int _0^{L\pi }\Theta (x)\mathrm {cos}\frac {2k x}{L}dx+ \int _0^{L\pi }\mathcal {R}_2(x)\mathrm {cos}\frac {2k x}{L}dx},\\[3pt]\ 0=d_2 {\int _0^{L\pi }Q_1^{\prime \prime }(x)+ \int _0^{L\pi }\mathcal {V}_2(x)\mathrm {cos}\frac {2k x}{L}dx},\\[3pt] P^{\prime }_1(x)=Q_1^{\prime }(x)=0. \end{array}} \right . \end{align}

We can obtain

(54) \begin{align} d_1\int _0^{L\pi }P^{\prime \prime }_1(x)\mathrm {cos}\frac {2k x}{L}dx=-4d_1\delta _k^2\int _0^{L\pi }P_1(x)\mathrm {cos}\frac {2k x}{L}dx, \end{align}

(55) \begin{align} d_2\int _0^{L\pi }Q^{\prime \prime }_1(x)\mathrm {cos}\frac {2k x}{L}dx=-4d_2\delta _k^2\int _0^{L\pi }Q_1(x)\mathrm {cos}\frac {2k x}{L}dx, \end{align}

(56) \begin{align} \int _0^{L\pi }\Theta (x)\mathrm {cos}\frac {2k x}{L}dx=4\xi _k^S\phi (P_{*})\delta _k^2\int _0^{L\pi }Q_1(x)\mathrm {cos}\frac {2k x}{L}dx +\frac {\delta _k^2L\pi }{2}\xi _k^S\phi _{P_{k}}(P_{*})\alpha _k, \end{align}

(57) \begin{align} \int _0^{L\pi }\mathcal {R}_2(x)\mathrm {cos}\frac {2k x}{L}dx=\mathcal {R}_{21}\int _0^{L\pi }P_1(x)\mathrm {cos}\frac {2k x}{L}dx+\mathcal {R}_{22}\int _0^{L\pi }Q_1(x)\mathrm {cos}\frac {2k x}{L}dx+\frac {L\pi }{4}\widetilde {\mathcal {R}}_2, \end{align}

and

(58) \begin{align} \int _0^{L\pi }\mathcal {V}_2(x)\mathrm {cos}\frac {2k x}{L}dx=\mathcal {V}_{21}\int _0^{L\pi }P_1(x)\mathrm {cos}\frac {2k x}{L}dx +\mathcal {V}_{22}\int _0^{L\pi }Q_1(x)\mathrm {cos}\frac {2k x}{L}dx+\frac {L\pi }{4}\widetilde {\mathcal {V}}_2. \end{align}

Then, submitting (54)–(58) into (53), one gets

\begin{align*} 0 = & \left (\mathcal {R}_{21}-4d_1\delta _k^2\right ) \int _0^{L\pi }P_1(x)\mathrm {cos}\frac {2k x}{L}dx +\left (\mathcal {R}_{22}+4\xi _k^S\phi (P_{*})\delta _k^2\right ) \\& \quad \int _0^{L\pi }Q_1(x)\mathrm {cos}\frac {2k x}{L}dx +\frac {\delta _k^2L\pi }{2}\xi _k^S\phi _{P_{k}}(P_{*})\alpha _k+\frac {L\pi }{4}\widetilde {\mathcal {R}}_2, \end{align*}

and

\begin{align*} 0=\mathcal {V}_{21}{\int _0^{L\pi }P_1(x)\mathrm {cos}\frac {2k x}{L}dx +\left (\mathcal {V}_{22}-\frac {4d_2k^2}{L^2}\right )\int _0^{L\pi }Q_1(x)\mathrm {cos}\frac {2k x}{L}dx+\frac {L\pi }{4}\widetilde {\mathcal {V}}_2}. \end{align*}

We thereby obtain

\begin{align*} \left ( \begin{array}{c@{\quad}c@{\quad}c} \mathcal {R}_{21}-4d_1\delta _k^2&\mathcal {R}_{22}+4\xi _k^S\phi (P_{*})\delta _k^2\\[3pt] \mathcal {V}_{21}&\mathcal {V}_{22}-4d_2\delta _k^2 \end{array} \right )\left ( \begin{array}{c@{\quad}c@{\quad}c} {\int _0^{L\pi } P_1(x)\mathrm {cos}\frac {2k x}{L}dx}\\[3pt] {\int _0^{L\pi } Q_1(x)\mathrm {cos}\frac {2k x}{L}dx} \end{array} \right )=\left ( \begin{array}{c} \displaystyle -\frac {\delta _k^2L\pi }{2}\xi _k^S\phi _{P_{k}}(P_{*})\alpha _k-\frac {L\pi }{4}\widetilde {\mathcal {R}}_2\\[3pt] \displaystyle -\frac {L\pi \widetilde {\mathcal {V}}_2}{4} \end{array} \right ). \end{align*}

Therefore, one achieves

(59) \begin{align} \int _0^{L\pi } P_1(x)\mathrm {cos}\frac {2k x}{L}dx=\frac {\left | \begin{array}{c@{\quad}c} \displaystyle -\frac {\delta _k^2L\pi }{2}\xi _k^S\phi _{P_{k}}(P_{*})\alpha _k-\frac {L\pi }{4}\widetilde {\mathcal {R}}_2 & \displaystyle R_{22}+4\xi _k^S\phi (P_{*})\delta _k^2\\[3pt] \displaystyle -\frac {L\pi \widetilde {\mathcal {V}}_2}{4}&\mathcal {V}_{22}-4d_2\delta _k^2 \end{array} \right |}{\left | \begin{array}{c@{\quad}c@{\quad}c} \mathcal {R}_{21}-4d_1\delta _k^2&\mathcal {R}_{22}+4\xi _k^S\phi (P_{*})\delta _k^2\\[3pt] \mathcal {V}_{21}&\mathcal {V}_{22}-4d_2\delta _k^2 \end{array} \right |}, \end{align}

and

(60) \begin{align} \int _0^{L\pi } Q_1(x)\mathrm {cos}\frac {2k x}{L}dx=\frac {\left | \begin{array}{c@{\quad}c} \mathcal {R}_{21}-4d_1\delta _k^2 & \displaystyle -\frac {\delta _k^2L\pi }{2}\xi _k^S\phi _{P_{k}}(P_{*})\alpha _k-\frac {L\pi }{4}\widetilde {\mathcal {R}}_2\\[3pt] \mathcal {V}_{21} & \displaystyle -\frac {L\pi \widetilde {\mathcal {V}}_2}{4} \end{array} \right |}{\left | \begin{array}{c@{\quad}c@{\quad}c} \mathcal {R}_{21}-4d_1\delta _k^2&\mathcal {R}_{22}+4\xi _k^S\phi (P_{*})\delta _k^2\\[3pt] \mathcal {V}_{21}&\mathcal {V}_{22}-4d_2\delta _k^2 \end{array} \right |}. \end{align}

Clearly, $\xi _2$ could be obtained by submitting (45), (52), (59), and (60) into (41).

Let $\xi _{k_0}^S=\mathrm {max}_{k\in \mathbb {N}_0\backslash \{0\}}\xi _k^S$ (see also Lemma3.1), then the validity of the second part of Theorem1.4 can be confirmed. Now, owing to Corollary 1.13 of [9], there exists an interval $I$ with $\xi _{k_0}^S\in I$ and $C^1-$ smooth function $(\xi, \varepsilon ):I\times ({-}\varrho, \varrho )\longrightarrow (\lambda _1(\xi ),\lambda _2(\varepsilon ))$ with $\lambda _1(\xi _{k_0}^S)=\lambda _2(0)=0$ and

(61) \begin{align} D_{(P,Q)}\mathcal {F}(P_{*},Q_{*},\xi )(P,Q)=\lambda _1(\xi )(P,Q), \end{align}

and

(62) \begin{align} D_{(P,Q)}\mathcal {F}(P_{k}(\varepsilon, x),Q_{k}(\varepsilon, x),\xi _k^S(\varepsilon ))(P,Q) =\lambda _2(\varepsilon )(P,Q), \end{align}

for $(P,Q)\in \mathbf {X}\times \mathbf {X}$ . It is not difficult to find that $\lambda _1(\xi )$ and $\lambda _2(\varepsilon )$ are eigenvalues of (61) and (62), respectively. The eigenfunction of the problem (61) could be represented by $(P(x,\xi ),Q(x,\xi ))$ and is uniquely described by ${(P(x,\xi _{k_0}^S),Q(x,\xi _{k_0}^S))}=\left (\mathrm {cos}\frac {k_0 x}{L},\alpha _{k_0}\mathrm {cos}\frac {k_0 x}{L}\right )$ . Also, $(P(x,\xi ),Q(x,\xi ))-\left (\mathrm {cos}\frac {k_0 x}{L},\alpha _{k_0}\mathrm {cos}\frac {k_0 x}{L}\right )\in \mathcal {K}$ is valid. Now from (6), we know that (61) takes the form

(63) \begin{align} \left \{ {\begin{array}{*{20}{l}} \displaystyle d_1P^{\prime \prime }-\xi (\phi (P_{*}) Q^{\prime })^{\prime } -\frac {bc^2P_{*}}{(cP_{*}+eQ_{*})^2}P +\left [dP_{*}-\frac {bceP_{*}}{(cP_{*}+eQ_{*})^2}\right ]Q=\lambda _1(\xi )P,\\[10pt] \displaystyle d_2 Q^{\prime \prime }-\left [dQ_{*}+\frac {bceQ_{*}}{(cP_{*}+eQ_{*})^2}\right ]P -\frac {be^2Q_{*}}{(cP_{*}+eQ_{*})^2}Q=\lambda _1(\xi )Q,\\[10pt] P^{\prime }=Q^{\prime }=0. \end{array}} \right . \end{align}

Differentiating (63) with respect to $\xi$ and then setting $\xi =\xi _{k_0}^S$ , one has

(64) \begin{align} \left \{ {\begin{array}{*{20}{l}} \displaystyle d_1\dot {P}^{\prime \prime }-\alpha _{k_0}\phi (P_{*}) \left (\mathrm {cos}\frac {k_0 x}{L}\right )^{\prime \prime }-\xi _{k_0}^S\phi (P_{*}) \dot {Q}^{\prime \prime } +f_P\dot {P}+f_Q\dot {Q}=\dot {\lambda }_1(\xi _{k_0}^S)\mathrm {cos}\frac {k_0 x}{L},\\[3pt] \displaystyle d_2\dot {Q}_2^{\prime \prime }+g_P\dot {P}+g_Q\dot {Q}=\dot {\lambda }_1(\xi _{k_0}^S)\alpha _{k_0}\mathrm {cos}\frac {k_0 x}{L},\\[3pt] \dot {P}^{\prime }=\dot {Q}^{\prime }=0, \end{array}} \right . \end{align}

where

\begin{align*} f_P=-\frac {bc^2P_{*}}{(cP_{*}+eQ_{*})^2}, \ f_Q=dP_{*}-\frac {becP^{*}}{(cP_{*}+eQ_{*})^2},\ g_P=-dQ_{*}-\frac {bceQ_{*}}{(cP_{*}+eQ_{*})^2}, \ g_Q= -\frac {be^2Q_{*}}{(cP_{*}+eQ_{*})^2}. \end{align*}

As a result, multiplying (64) by $\mathrm {cos}\frac {k_0 x}{L}$ , we have

(65) \begin{align} \left ( \begin{array}{c@{\quad}c} \displaystyle f_P-d_1\delta _{k_0}^2 & \displaystyle f_Q+\delta _{k_0}^2\xi _{k_0}^S\phi (P_{*})\\[6pt] \displaystyle g_P & \displaystyle g_Q-d_2\delta _{k_0}^2 \end{array} \right )\left ( \begin{array}{c} \displaystyle {\int _0^{L\pi } \dot {P}\mathrm {cos}\frac {k_0 x}{L}dx}\\[6pt] \displaystyle {\int _0^{L\pi } \dot {Q}\mathrm {cos}\frac {k_0 x}{L}dx} \end{array} \right )=\left ( \begin{array}{c} \displaystyle \dot {\lambda }_1(\xi _{k_0}^S)\frac {L\pi }{2}-\frac {\delta _{k_0}^2\alpha _{k_0}\phi (P_{*})L\pi }{2}\\[6pt] \displaystyle \dot {\lambda }_1(\xi _{k_0}^S)\alpha _{k_0}\frac {L\pi }{2} \end{array} \right ). \end{align}

Obviously, the coefficient matrix in (65) is singular since $\xi =\xi _{k_0}^S$ . This implies that

\begin{align*} \frac {f_P-d_1\delta _{k_0}^2}{g_P}= \frac {\dot {\lambda }_1(\xi _{k_0}^S)-\delta _{k_0}^2\alpha _{k_0}\phi (P_{*})} {\dot {\lambda }_1(\xi _{k_0}^S)\alpha _{k_0}}. \end{align*}

Consequently, one obtains

\begin{align*} \dot {\lambda }_1(\xi _{k_0}^S)=-\frac {\delta _{k_0}^2\alpha _{k_0}\phi (P_{*})g_P}{(f_P-d_1\delta _{k_0}^2)\alpha _{k_0}-g_P}\lt 0. \end{align*}

Using Theorem 1.16 of [9], $\lambda _2(\varepsilon )$ and $-\varepsilon \xi ^{\prime S}_{k_0}(\varepsilon )\dot {\lambda }_1(\xi _{k_0}^S)$ have the same sign near $\varepsilon =0$ . As a result, we can compute $\textrm {Sign}(\lambda _2(\varepsilon ))=\textrm {Sign}({-}\varepsilon \xi ^{\prime S}_{k_0}(\varepsilon )\dot {\lambda }_1(\xi _{k_0}^S)) =\textrm {Sign}({-}2\varepsilon ^2\xi _2\dot {\lambda }_1(\xi _{k_0}^S))=\textrm {Sign}(\xi _2)$ . This means that the bifurcating solution $\mathcal {S}_{k_0}(\varepsilon )=(P_{k_0}(\varepsilon, x),Q_{k_0}(\varepsilon, x))$ is asymptotically stable when $\xi _2\lt 0$ and it is unstable when $\xi _2\gt 0$ for $\varepsilon \in ({-}\varrho, \varrho )$ . The proof readily follows.

4. Numerical simulations

In this section, we will describe the numerical solution algorithms to solve the IGP-type predator–prey model (1). We shall conduct various computational simulations to confirm the validity of Theorem1.4. Our main aim is to perform simulations for the stable nonconstant steady states around the steady state bifurcation onset $\xi =\xi _{k_0}^S$ . More precisely, we want to find not only the nonconstant steady states in traditional 1D and 2D domains but also on spherical and torus surfaces. Now, some specific parameters in (1) are fixed

(66) \begin{align} e=1,\alpha =1.5,c=1, \beta =0.2,h=0.05,d=0.85,b=0.65, d_1=0.85, d_2=0.5. \end{align}

4.1 Nonconstant steady states exist in 1D space

Let us consider a one-dimensional space $\Omega = (0,L_x)$ . By using a cell-centered grid, the uniform discrete computational domain is defined as $\Omega _d = \big \{ x_i\vert (i-0.5)\Delta x,\,1\leq i \leq N_x\big \}$ , where $\Delta x = L_x/N_x$ and $N_x$ is the number of discrete points. Numerical approximations $P(x_i,n\Delta t)$ and $Q(x_i,n\Delta t)$ are denoted by $P_i^n$ and $Q_i^n$ , respectively. Here, $\Delta t$ is time step size and $n=0,1,\cdots$ . On the discrete computational domain $\Omega _d$ , the IGP-type predator–prey model (1) can be discretized using the explicit Euler method as follows:

\begin{eqnarray*} \left \{ {\begin{array}{*{20}{l}} \displaystyle \frac {P_i^{n+1}-P_i^n}{\Delta t} = d_1 \Delta _d P_i^n - \nabla _d \cdot \big (\xi \phi (P_i^n)\nabla _d Q_i^n\big ) + P_i^n \left (\frac {bc}{cP_i^n+eQ_i^n}+dQ_i^n-\alpha \right ),\\[10pt] \displaystyle \frac {Q_i^{n+1}-Q_i^n}{\Delta t} = d_2 \Delta _d Q_i^n + P_i^n \left (\frac {be}{cP_i^n+eQ_i^n}-dP_i^n-\beta \right )-hQ_i^n, \end{array}} \right . \end{eqnarray*}

where $ \Delta _d P_i^n = (P_{i+1}^n-2P_i^n+P_{i-1}^n)/\Delta x^2$ and $\Delta _d Q_i^n=(Q_{i+1}^n-2Q_i^n+Q_{i-1}^n)/\Delta x^2$ are the discrete Laplacian operators [23]. We use the conservative discretization for the term $\nabla _d \cdot \big (\xi \phi (P_i^n)\nabla _d Q_i^n\big )$ as follows:

\begin{eqnarray*} \nabla _d \cdot \big (\xi \phi (P_i^n)\nabla _d Q_i^n\big ) = \frac {\xi }{\Delta x^2}\bigg [\phi (P_{i+\frac {1}{2}}^n)(Q_{i+1}^n-Q_i^n)-\phi (P_{i-\frac {1}{2}}^n)(Q_i^n-Q_{i-1,j}^n)\bigg ], \end{eqnarray*}

where $P_{i+\frac {1}{2}} = (P_{i+1,j}+P_i)/2$ and $P_{i-\frac {1}{2}} = (P_{i}+P_{i-1})/2.$

We numerically solve the model (1) to validate the nonconstant steady-state solution in a one-dimensional space $\Omega = (0,8\pi )$ with $N_x=256$ points, a uniform spatial grid size of $\Delta x=8\pi /N_x$ , and a time step of $\Delta t=0.2\Delta x^2$ . First, let us take the density function $\phi (P)=P$ . Clearly, the assumption (H2) is satisfied. Next, we choose the parameters in (66) and the spatial domain $\Omega =(0,8\pi )$ . Then, we know that (H3) holds, namely, there is a unique positive equilibrium $E_{*}=(0.2259,1.2447)$ and

\begin{align*} &\xi _1^S\approx -41.2976,\ \xi _2^S\approx -11.2173,\ \xi _3^S\approx -5.7382,\ \xi _4^S\approx -3.9282,\ \xi _5^S\approx -3.2087,\\[3pt] &\xi _6^S\approx -2.9433,\ \xi _7^S\approx -2.9140,\ \xi _8^S\approx -3.0297,\ \xi _9^S\approx -3.2469,\ \xi _{10}^S\approx -3.5427,\\[3pt] &\xi _{11}^S\approx -3.9040,\ \xi _{12}^S\approx -4.3231,\ \xi _{13}^S\approx -4.7950,\ \xi _{14}^S\approx -5.3165, \ \cdots . \end{align*}

Consequently, we have $\xi _{k_0}^S=\mathrm {max}_{k\in \mathbb {N}_0\backslash \{0\}} \xi _k^S=-2.9140$ . To display the stable nonconstant steady states, we choose $-3.5=\xi \lt \xi _{k_0}^S$ and the initial data $(P_0(x),Q_0(x))=(0.2259-0.02 \mathrm {cos}\left (\frac {7 x}{8}\right ),1.2447-0.02\mathrm {cos}\left (\frac {7 x}{8}\right ))$ . Our numerical simulations indicate that there are stable nonconstant steady states, see Figure 1, where pictures (a)–(b) for predator $P(x,t)$ and prey $Q(x,t)$ and picture (c) for their space series diagrams.

In the sequel, we assume that the density function $\phi (P)$ satisfies the saturated effect, i.e., one takes $\phi (P)=\frac {P}{1+P}$ . We can show that $\phi (P)=\frac {P}{1+P}\lt P$ and the condition (H2) is satisfied. For the parameters performed in (66) and the spatial length $\Omega =(0,8\pi )$ . Then, we know that (H3) holds, namely, there exists a unique positive equilibrium $E_{*}=(0.2259,1.2447)$ and

\begin{align*} &\xi _1^S\approx -50.6260,\ \xi _2^S\approx -13.7511,\ \xi _3^S\approx -7.0343,\ \xi _4^S\approx -4.8155,\ \xi _5^S\approx -3.9335,\\[3pt] &\xi _6^S\approx -3.6082,\ \xi _7^S\approx -3.5722,\ \xi _8^S\approx -3.7141,\ \xi _9^S\approx -3.9803,\ \xi _{10}^S\approx -4.3429,\\[3pt] &\xi _{11}^S\approx -4.7858,\ \xi _{12}^S\approx -5.2996,\ \xi _{13}^S\approx -5.8781,\ \xi _{14}^S\approx -6.5174,\ \cdots . \end{align*}

Accordingly, $\xi _{k_0}^S=\mathrm {max}_{k\in \mathbb {N}_0\backslash \{0\}}\xi _k^S=-3.5722$ . To display the stable nonconstant steady states, we choose $-3.75=\xi \lt \xi _{k_0}^S$ and the initial data $(P_0(x),Q_0(x))=(0.2259-0.02\mathrm {cos}\left (\frac {7 x}{8}\right ),1.2447-0.02\mathrm {cos}\left (\frac {7 x}{8}\right ))$ . Then, the numerical simulations show that there are stable nonconstant steady states because the density function $\phi (P)$ takes the saturated form, see Figure 2, where pictures (a)–(b) for predator $P(x,t)$ and prey $Q(x,t)$ and picture (c) for their space series diagrams.

Figure 1. Taking the density function $\phi (P)=P,\xi =-3.5$ and the other parameters are fixed in (66), system (1) admits the stable nonconstant steady states, where the initial data $(P_0(x),Q_0(x))=(0.2259-0.02\mathrm {cos}(\frac {7 x}{8}),1.2447-0.02\mathrm {cos}(\frac {7 x}{8}))$ .

Figure 2. Taking the density function $\phi (P)=\frac {P}{1+P},\xi =-3.75$ and the other parameters are fixed in (66), system (1) admits the stable nonconstant steady states, where the initial data $(P_0(x),Q_0(x))=(0.2259-0.02\mathrm {cos}(\frac {7 x}{8}),1.2447-0.02\mathrm {cos}(\frac {7 x}{8}))$ .

Next, let us suppose that the density function $\phi (P)$ with the Ricker effect, specifically, $\phi (P)=Pe^{-P}$ . Clearly, the condition (H2) is satisfied. Now, we maintain the same parameters and the spatial length as in Figures 1 and 2. Then, the unique positive equilibrium $E_{*}=(0.2259,1.2447)$ and

\begin{align*} &\xi _1^S\approx -51.7636,\ \xi _2^S\approx -14.0601,\ \xi _3^S\approx -7.1924,\ \xi _4^S\approx -4.9237,\ \xi _5^S\approx -4.0219,\\[3pt] &\xi _6^S\approx -3.6892,\ \xi _7^S\approx -3.6525,\ \xi _8^S\approx -3.7975,\ \xi _9^S\approx -4.0697,\ \xi _{10}^S\approx -4.4405,\\[3pt] &\xi _{11}^S\approx -4.8934,\ \xi _{12}^S\approx -5.4187,\ \xi _{13}^S\approx -6.0102,\ \xi _{14}^S\approx -6.6639,\ \cdots . \end{align*}

Figure 3. Taking the density function $\phi (P)=Pe^{-P},\xi =-3.85$ and the other parameters are fixed in (66), system (1) admits the stable nonconstant steady states, where the initial data $(P_0(x),Q_0(x))=(0.2259-0.02\mathrm {cos}(\frac {7 x}{8}),1.2447-0.02\mathrm {cos}(\frac {7 x}{8}))$ .

It is found that $\xi _{k_0}^S=\mathrm {max}_{k\in \mathbb {N}_0 \backslash \{0\}} \xi _k^S=-3.6525$ . We choose $-3.85=\xi \lt \xi _{k_0}^S$ and the initial data $(P_0(x),Q_0(x))=(0.2259-0.02\mathrm {cos}\left (\frac {7 x}{8}\right ),1.2447-0.02 \mathrm {cos}\left (\frac {7 x}{8}\right ))$ . Then, numerical simulations show that there are stable nonconstant steady states, as shown in Figure 3, where pictures (a)–(b) for predator $P(x,t)$ and prey $Q(x,t)$ and picture (c) for their space series diagrams.

4.2 Nonconstant steady states exist in 2D space

In this subsection, we investigate nonconstant steady states in the two-dimensional computational domain $\Omega = (0,L_x)\times (0,L_y)$ . We define the discrete computational domain $\Omega _d = \big \{(x_i,y_j)\vert \big ( (i-0.5)\Delta x,\,(j-0.5)\Delta y),\,1\leq i\leq N_x,\,1\leq j \leq N_y\big \}$ , where $\Delta x =L_x/N_x$ and $\Delta y = L_y/N_y$ ; and $N_x$ and $N_y$ are the number of grid points in the $x$ - and $y$ -directions, respectively. Numerical approximations $P(x_i,y_j,n\Delta t)$ and $Q(x_i,y_j,n\Delta t)$ are denoted by $P_{ij}^n$ and $Q_{ij}^n$ , respectively. Here, $\Delta t$ is time step size and $n=0,1,\cdots$ . By using the explicit Euler method, the IGP-type predator– prey model (1) in the two-dimensional domain can be discretized as follows:

(67) \begin{eqnarray} \left \{ {\begin{array}{*{20}{l}} \displaystyle \frac {P_{ij}^{n+1} - P_{ij}^n}{\Delta t} = d_1 \Delta _d P_{ij}^n - \nabla _d \cdot (\xi \phi (P_{ij}^n)\nabla _d Q_{ij}^n) + P_{ij}^n\left (\frac {bc}{cP_{ij}^n+eQ_{ij}^n}+dQ_{ij}^n-\alpha \right ),\\[3pt] \displaystyle \frac {Q_{ij}^{n+1} - Q_{ij}^n}{\Delta t} = d_2 \Delta _d Q_{ij}^n + Q_{ij}^n\left (\frac {be}{cP_{ij}^n+eQ_{ij}^n}-dP_{ij}^n-\beta \right )-hQ_{ij}^n, \end{array}} \right . \end{eqnarray}

where the two-dimensional discrete Laplacian operators are defined as follows [24]:

\begin{eqnarray*} \Delta _d P_{ij}^n = \frac {P_{i+1,j}^n-2P_{ij}^n+P_{i-1,j}^n}{\Delta x^2}+\frac {P_{i,j+1}^n-2P_{ij}^n+P_{i,j-1}^n}{\Delta y^2} \end{eqnarray*}

and

\begin{eqnarray*} \Delta _d Q_{ij}^n = \frac {Q_{i+1,j}^n-2Q_{ij}^n+Q_{i-1,j}^n}{\Delta x^2}+\frac {Q_{i,j+1}^n-2Q_{ij}^n+Q_{i,j-1}^n}{\Delta y^2}. \end{eqnarray*}

We use the conservative form to define the term $\nabla _d \cdot (\xi \phi (P_{ij}^n)\nabla _d Q_{ij}^n)$ as follows:

\begin{eqnarray*} \nabla _d \cdot (\xi \phi (P_{ij}^n)\nabla _d Q_{ij}^n) &=&\xi \bigg [\frac {1}{\Delta x^2}\Big \{\phi (P_{i+\frac {1}{2},j}^n)(Q_{i+1,j}^n-Q_{ij}^n)-\phi (P_{i-\frac {1}{2},j}^n)(Q_{ij}^n-Q_{i-1,j}^n)\Big \}\\[3pt] &&+\frac {1}{\Delta y^2}\Big \{\phi (P_{i,j+\frac {1}{2}}^n)(Q_{i,j+1}^n-Q_{ij}^n)-\phi (P_{i,j-\frac {1}{2}}^n)(Q_{ij}^n-Q_{i,j-1}^n)\Big \}\bigg ], \end{eqnarray*}

where $ P_{i+\frac {1}{2},j} = \left (P_{i+1,j}+P_{ij}\right )/2, \ P_{i,j+\frac {1}{2}} = \left (P_{i,j+1}+P_{ij}\right )/2, \ P_{i-\frac {1}{2},j} = \left (P_{ij}+P_{i-1,j}\right )/2, \ P_{i,j-\frac {1}{2}} = \left (P_{ij}+P_{i,j-1}\right )/2.$ From Eqs. (67), we can obtain numerical solutions as

\begin{eqnarray*} \left \{ {\begin{array}{*{20}{l}} \displaystyle P_{ij}^{n+1}= P_{ij}^n +{\Delta t}\left [ d_1 \Delta _d P_{ij}^n - \nabla _d \cdot (\xi \phi (P_{ij}^n)\nabla _d Q_{ij}^n) + P_{ij}^n\left (\frac {bc}{cP_{ij}^n+eQ_{ij}^n}+dQ_{ij}^n-\alpha \right )\right ],\\[3pt] \displaystyle Q_{ij}^{n+1} = Q_{ij}^n + {\Delta t}\left [ d_2 \Delta _d Q_{ij}^n + Q_{ij}^n\left (\frac {be}{cP_{ij}^n+eQ_{ij}^n}-dP_{ij}^n-\beta \right )-hQ_{ij}^n\right ]. \end{array}} \right . \end{eqnarray*}

For the two-dimensional numerical simulations, discrete $l_2$ -errors for $P$ and $Q$ are defined as

\begin{eqnarray*} E_P^{n} = \sqrt {\frac {1}{N_x N_y} \sum _{i=1}^{N_x} \sum _{j=1}^{N_y} (P_{ij}^{n}-P_{ij}^{n-1})^2}\ \textrm {and}\ E_Q^{n} = \sqrt {\frac {1}{N_x N_y} \sum _{i=1}^{N_x} \sum _{j=1}^{N_y} (Q_{ij}^{n}-Q_{ij}^{n-1})^2}, \end{eqnarray*}

respectively. We define the numerical steady states $P^s$ and $Q^s$ when the average of errors is less than a tolerance; $0.5(E_P^s+E_Q^s)\lt tol$ . In the following numerical experiments, we use a tolerance of $tol=1.953e$ - $9$ . Now, we consider the following random perturbed initial condition:

(68) \begin{eqnarray} \left \{ {\begin{array}{*{20}{l}} P(x,y,0) = P_{*} + 0.02\textrm {rand}(x,y),\\[3pt] Q(x,y,0) = Q_{*} + 0.02\textrm {rand}(x,y), \end{array}} \right . \end{eqnarray}

where $\textrm {rand}(x,y)$ is a random variable between $-1$ and 1. We use a uniform mesh with $N_x=N_y=128$ , $\Delta x = \Delta y = 25/128$ and time step $\Delta t=0.2\Delta x^2$ on the computational domain $\Omega = (0,25) \times (0,25)$ .

In Figure 4, we choose the density function $\phi (P)=P$ with the parameters in (66). As a consequence, we have the critical value of the steady-state bifurcation as $\xi _{k_0}^S=-2.9140$ . Next, we set the prey-taxis sensitivity constant $\xi =-3.5$ around the critical value $\xi _{k_0}^S$ . It is found that the spotted pattern occupies the bounded domain $\Omega = (0,25)\times (0,25)$ as time progresses.

Next, we suppose that the density function with the saturated form, $\phi (P)=\frac {P}{1+P}$ in the IGP-type predator–prey model (1). To observe the nonconstant steady state of the predator–prey model (1) under this density function, we set the parameters as in (66). In this case, the critical value of the onset of steady state bifurcation is $\xi _{k_0}^S=-3.5722$ . Now, we set the prey-taxis sensitivity coefficient as $\xi =-3.75$ . Using these parameters, a combination of stripe and spot patterns (mixed patterns) can be found in the bounded domain $\Omega = (0,25)\times (0,25)$ , as shown in Figure 5.

Figure 4. From top to bottom rows, temporal evolutions of patterns for $P$ and $Q$ . Taking the density function $\phi (P)=P,\xi =-3.5$ and the other parameters are fixed in (66), system (1) admits the nonconstant steady states with the initial data (68).

Figure 5. From top to bottom rows, temporal evolutions of patterns for $P$ and $Q$ . Taking the density function $\phi (P)=\frac {P}{1+P},\xi =-3.75$ and the other parameters are fixed in (66), system (1) admits the nonconstant steady states with the initial data (68).

Figure 6. From top to bottom rows, temporal evolutions of patterns for $P$ and $Q$ . Taking the density function $\phi (P)=Pe^{-P},\xi =-3.85$ and the other parameters are fixed in (66), system (1) admits the nonconstant steady states with the initial data (68).

We would like to mention that a similar pattern formation can be shown in Figure 6, where we choose the density function $\phi (P)=Pe^{-P}$ when choosing the parameter (66) and $\xi =-3.85$ .

4.3 Nonconstant steady states on spherical and torus surfaces

In this subsection, we will perform the nonconstant steady states of the IGP-type predator–prey model (1) on both spherical and torus surfaces. To this end, let us first illustrate the discrete computational domains for the spherical and torus surfaces. On a closed smooth surface $\mathcal {S}$ , to numerically investigate pattern formations of the governing system, a triangular surface mesh $\mathcal {S}_d$ is used, see Figure 7 (a). We discretize the Laplace–Beltrami operator using an approach in [10, 17]. We define a surface point set $ \left \{\textbf {x}_{i}\right \}_{i=1}^{m} = \left \{\textbf {x}_1, \textbf {x}_2, \textbf {x}_3, \ldots, \textbf {x}_{m} \right \}$ on a triangular surface mesh $\mathcal {S}_{d}$ . Then, each vertex point $\textbf {x}_{i}$ has one-ring triangular surface points with an index set $I(i)=\{i_{1}, i_{2}, \cdots, i_p \}$ with $i_{1}=i_p$ , see Figure 7(b).

Figure 7. Schematic visualizations: (a) triangulated mesh of discretized spherical surface $\mathcal {S}_d$ , (b) surrounding one-ring surface points set for $\textbf {x}_i$ , (c) triangles $T_{j}$ and $T_{j+}$ featuring the angles $\alpha _{ij}$ and $\beta _{ij_{+}}$ and (d) vertex ${\textbf x}_{i}$ and its corresponding area $ {\mathcal A}({\textbf x}_{i})$ .

The discrete numerical approximation is denoted as $P_i^n = P(\textbf {x}_i,n\Delta t)$ , where $\Delta t$ is the time step size. We discretize the IGP-type predator–prey model (1) as follows:

\begin{eqnarray*} \left \{ {\begin{array}{*{20}{l}} \displaystyle \frac {P_i^{n+1}-P_i^n}{\Delta t} = d_1 \Delta _{\mathcal {S}} P_i^n - \nabla _{\mathcal {S}} \cdot \big (\xi \phi (P_i^n)\nabla _{\mathcal {S}} Q_i^n\big ) + P_i^n \left (\frac {bc}{cP_i^n+eQ_i^n}+dQ_i^n-\alpha \right ),\\[3pt] \displaystyle \frac {Q_i^{n+1}-Q_i^n}{\Delta t} = d_2 \Delta _{\mathcal {S}} Q_i^n + Q_i^n \left (\frac {be}{cP_i^n+eQ_i^n}-dP_i^n-\beta \right )-hQ_i^n. \end{array}} \right . \end{eqnarray*}

Here, we consider the discrete Laplace–Beltrami operator defined as

\begin{align*} \Delta _{\mathcal S} P_i=\frac {3}{\mathcal {A}(\textbf {x}_i)}\sum _{j\in I(i)} \frac {\cot \alpha _{ij}+\cot \beta _{ij}}{2} (P_j-P_i) \textrm { and } \Delta _{\mathcal S} Q_i=\frac {3}{\mathcal {A}(\textbf {x}_i)}\sum _{j\in I(i)} \frac {\cot \alpha _{ij}+\cot \beta _{ij}}{2} (Q_j-Q_i), \end{align*}

where ${\mathcal A}(\textbf {x}_{i})$ is the cumulative area for the individual triangles $T_{j}$ centered around surface point $\textbf {x}_{i}$ (Figure 7(d)):

\begin{eqnarray*} {\mathcal A}(\textbf {x}_{i})= \sum _{j \in I(i)} \frac {\sqrt {||\textbf {x}_{j}-\textbf {x}_{i}||^{2}||\textbf {x}_{j_{+}}-\textbf {x}_{i}||^{2}-\left (\textbf {x}_{j}-\textbf { x}_{i},\textbf {x}_{j_{+}}-\textbf {x}_{i}\right )^{2}}}{2}. \end{eqnarray*}

The discrete divergence term $\nabla _{\mathcal {S}} \cdot \left ( \xi \phi (P_i)\nabla _{\mathcal {S}} Q_i\right )$ is approximated using a conservative form as follows:

\begin{eqnarray*} \nabla _{\mathcal {S}} \cdot \left ( \xi \phi (P_i)\nabla _{\mathcal {S}} Q_i\right ) = \frac {3\xi }{\mathcal {A}(\textbf {x}_i)}\sum _{j\in I(i)} \frac {\cot \alpha _{ij}+\cot \beta _{ij}}{2}\phi \bigg (\frac {P_j+P_i}{2}\bigg ) (Q_j-Q_i). \end{eqnarray*}

In the following numerical experiments, for the numerical simulations on triangular surfaces, we use the following randomly perturbed initial condition:

(69) \begin{eqnarray} \left \{ {\begin{array}{*{20}{l}} P(\textbf {x}_i,0) = P_{*} + 0.02\textrm {rand}(\textbf {x}_i),\\[3pt] Q(\textbf {x}_i,0) = Q_{*} + 0.02\textrm {rand}(\textbf {x}_i), \end{array}} \right . \end{eqnarray}

where $\textrm {rand}(\textbf {x}_i)$ is the uniformly distributed random perturbation between $-1$ and $1$ .

4.3.1 Nonconstant steady states on the spherical surface

We consider a triangulated spherical surface mesh $\mathcal {S}_d$ with a radius value $r=15$ and the number of triangulated spherical surface points is $16590$ .

Figure 8. From top to bottom rows, temporal evolutions for patterns of $P$ and $Q$ on the spherical surface. Taking the density function $\phi (P)=P,\xi =-3.65$ and the other parameters are fixed in (66), system (1) admits the nonconstant steady states with the initial data (69).

In Figure 8, we choose the density function $\phi (P)=P$ and parameter (66) in the IGP-type predator–prey model (1). Based on the theoretical analysis, we have the critical value of the steady state bifurcation to be $\xi _{k_0}^S=-2.9140$ . Accordingly, we take the prey-taxis sensitivity parameter $\xi =-3.65$ . Considering these known parameters, we can observe that spot patterns can be formed on the spherical surface as time progresses.

Figure 9. From top to bottom rows, temporal evolutions for patterns of $P$ and $Q$ on the spherical surface. Taking the density function $\phi (P)=\frac {P}{1+P},\xi =-3.75$ and the other parameters are fixed in (66), system (1) admits the nonconstant steady states with the initial data(69).

Next, let us consider the density function $\phi (P)=\frac {P}{1+P}$ in the IGP predator–prey model (1) and fix the parameters in (66) and $d_2=0.5$ . Through direct calculation, we have the steady-state bifurcation threshold as $\xi _{k_0}^S=-3.5722$ . Thus, we take $\xi =-3.75$ in proximity to the onset $\xi _{k_0}^S=-3.5722$ . Our numerical simulations demonstrate that the IGP-type predator–prey model (1) exhibits mixed stripe and spot patterns on the spherical surface, as shown in Figure 9.

Finally, one takes the Ricker form density function $\phi (P)=Pe^{-P}$ in the IGP predator–prey model (1). Moreover, the parameters are set in (66). As a consequence, we obtain the threshold for the steady- state bifurcation to be $\xi _{k_0}^S=-3.6525$ . Taking the prey-taxis coefficient to $\xi =-3.85$ , our numerical simulation results suggest that the IGP predator–prey model (1) exhibits mixed patterns on the spherical surface, as shown in Figure 10.

Figure 10. From top to bottom rows, temporal evolutions for pattern formation of $P$ and $Q$ on the spherical surface. Taking the density function $\phi (P)=Pe^{-P},\xi =-3.85$ and the other parameters are fixed in (66), system (1) admits the nonconstant steady states with the initial data (69).

4.3.2 Nonconstant steady states on the torus surface

Now, we will explore the nonconstant steady states of the IGP-type predator–prey model (1) on the tour’s surface. To this end, let us consider a triangulated torus surface mesh $\mathcal {S}_d$ , which has a major radius (the distance from the center of the tube to the center of the torus) value of $R=15$ , a minor radius (radius of the tube) value of $r=20$ , and the number of triangulated spherical surface points is $16544$ .

Figure 11. From top to bottom rows, temporal evolutions for patterns of $P$ and $Q$ on the torus surface. Taking the density function $\phi (P)=P,\xi =-3.65$ and the other parameters are fixed in (66), system (1)admits the nonconstant steady states with the initial data (69).

In Figure 11, we take the density function $\phi (P)=P$ and the parameters fixed in (66) in the IGP predator–prey model (1). We can get the critical value of the steady-state bifurcation as $\xi _{k_0}^S=-2.9140$ . Now, let us keep the prey-taxis sensitivity coefficient $\xi =-3.65$ , then there exist the spot patterns of the IGP predator–prey model (1) on the torus surface.

Figure 12. From top to bottom rows, temporal evolutions for patterns of $P$ and $Q$ on the torus surface. Taking the density function $\phi (P)=\frac {P}{1+P},\xi =-3.75$ and the other parameters are fixed in (66), system (1) admits the nonconstant steady states with the initial data (69).

Figure 12 suggests that the IGP predator–prey model (1) possesses the mixed patterns on the torus surface when choosing the saturated form density function $\phi (P)=\frac {P}{1+P}$ and the parameter values in (66) and $\xi =-3.75$ .

Figure 13. From top to bottom rows, temporal evolutions for patterns of $P$ and $Q$ on the torus surface. Taking the density function $\phi (P)=Pe^{-P},\xi =-3.85$ and the other parameters are fixed in (66), system (1) admits the nonconstant steady states with the initial data (69).

Similar pattern formations on the torus surface can be found in Figure 13, where one adopts the Ricker-type density function $\phi (P)=Pe^{-P}$ and the parameter values in (66) and $\xi =-3.85$ .

4.4 Influence of the harvesting coefficient $h$

Now, we keep the same density function $\phi (P)$ and parameters in Figure 3 (see also Figures 6,10,13) but change the harvesting coefficient $h$ to display how the harvesting coefficient $h$ will affect the pattern formation dynamics of the IGP-type predator–prey model (1). When there is no harvesting effect, namely, $h=0$ , system (1) exhibits the nonconstant steady states (stripe patterns), see picture (a) of Figures 14–16, respectively. Similar pictures can be found in pictures (b) and (c) of Figures 14– 16. A clear fact is that the stripe patterns (nonconstant steady states) gradually diminish as the harvesting coefficient $h$ progressively increases.

As the harvesting coefficient $h$ increases, the stripe patterns gradually disappear, as shown in Figures 14(d)–(f), 15(d)–(f), and 16(d)–(f). In fact, with the continuous increase of the harvesting coefficient $h$ , $\xi =-3.85$ gradually moves further away from the steady-state bifurcation threshold (see Figures 3, 10, 13), this leads to the change of the patterns in Figure 14, Figure 15, and Figure 16. Consequently, prey harvesting plays an important role in inducing spatial patterns. Ecologically, over-harvesting for prey or predators is not advisable, it can lead to an ecological imbalance due to a significant reduction in population numbers. However, harvesting within a certain range is a feasible approach. This harvesting strategy is consistent with reality.

Figure 14. Influence of the harvesting coefficient $h$ when $\phi (P)=Pe^{-P}$ , where we choose  $e=1.0,\alpha =1.5,c=1,\beta =0.2,d=0.85,b=0.65,d_1=0.85,d_2=0.5,\xi =-3.85$ and the initial data are $(P_0(x),Q_0(x))=(P_{*}-0.02\mathrm {cos}(\frac {7 x}{8}),Q_{*}-0.02\mathrm {cos}(\frac {7 x}{8}))$ .

Figure 15. Influence of the harvesting coefficient $h$ on the spherical surface when $\phi (P)=Pe^{-P}$ , where $e=1.0,\alpha =1.5,c=1,\beta =0.2,d=0.85,b=0.65,d_1=0.85,d_2=0.5,\xi =-3.85$ and the initial data is (69).

Figure 16. Influence of the harvesting coefficient $h$ on the torus surface when $\phi (P)=Pe^{-P}$ , where $e=1.0,\alpha =1.5,c=1,\beta =0.2,d=0.85,b=0.65,d_1=0.85,d_2=0.5,\xi =-3.85$ and the initial data is (69).

In summary, we have displayed the emergence of nonconstant steady states in 1D and 2D spaces, as well as on the spherical and torus surfaces. These numerical results are in good agreement with our previous theoretical analysis. Moreover, we can conclude that prey-taxis and harvesting effects will induce wealthy pattern dynamics for the IGP-type predator–prey model.

5. Discussions

This paper reports the existence of the classical solution and spatiotemporal dynamics of an IGP-type predator–prey model incorporating both harvesting and prey-taxis. First, we discuss the local-in time and global existence of the classical solution $(P(x,t),Q(x,t))$ of the IGP-type predator–prey model (1) in $N$ -dimensional space by virtue of some estimates, Amann’s theorem, and Neumann heat semigroup theory. Especially, it is found that the classical solution $(P(x,t),Q(x,t))$ of IGP-type predator–prey model (1) exists for small prey-taxis sensitivity coefficient $\xi$ as the dimension of space $N$ is large, see Theorem1.1. Thereafter, we explore the steady-state bifurcation of the model (1). To this end, the stability analysis of the positive equilibrium $E_{*}$ is first discussed, see Theorem1.2. Our theoretical result demonstrates that the unique positive equilibrium $E_{*}$ is locally asymptotically stable for any $\xi \geq 0$ and the predator–prey model (1) suffers from the steady-state bifurcation when $\xi =\xi _k^S$ , where $\xi _k^S\lt 0$ for all $k\in \mathbb {N}_0$ . Accordingly, the repulsive prey-taxis could destabilize the spatial homogeneity of the IGP-type predator–prey model (1), while the attractive prey-taxis effect will stabilize the spatial homogeneity. It is not difficult to find that there only is self-diffusion involved in the system (1) when $\xi =0$ . Self-diffusion means that the movements of both predators and prey are random. For this case, the unique positive equilibrium $E_{*}$ always keeps its local asymptotic stability. This implies that predators and prey will coexist and random movement can not change the stable structure of the system (1). However, such a homogeneous stable status can be destroyed by integrating the prey-taxis into the system (1).

In what follows, with the help of the Crandall–Rabinowitz local bifurcation theory, we respectively establish the existence and stability of the bifurcating solution, which resulted from the steady state bifurcation, see Theorem1.3 and Theorem1.4, respectively. Finally, numerical experiments are conducted to verify our theoretical analysis by choosing the different density functions $\phi (P)$ . In light of the theoretical results, we find the stripe patterns of the IGP predator–prey model (1) in the 1D domain (see Figures 1–3) and the spot patterns and stripe-spot mixed patterns in the 2D domain (see Figures 4–6). Interestingly, these complicated pattern formations can also be observed on the spherical surface (see Figures 8–10) and torus surface (see Figures 11–13). These numerical results are in good agreement with the theoretical analysis. Of course, one plots Figures 14–16 to explore the influence of the harvesting effect of the IGP-type predator–prey model (1). It is found that the spatial patterns will gradually disappear with the continuous increase of the harvesting coefficient $h$ . This phenomenon may enlighten us that over-harvesting for prey or predators is not advisable, which will lead to ecological imbalance due to the drastic reduction in population numbers. Overall, we can conclude that this IGP-type predator–prey model with the prey-taxis and harvesting effects will perform the wealthy and interesting dynamic profiles. These results may be useful for exploring and understanding the complex dynamical evolution among different populations in a harvesting and prey-taxis environment.