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 [Reference Shi and Wang33] or Lemma2.3 of [Reference Wang and Xu36], 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 [Reference Crandall and Rabinowitz9], 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 [Reference Crandall and Rabinowitz9], 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 [Reference Crandall and Rabinowitz9], $\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.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 [Reference Kwak, Kang, Ham, Hwang, Lee and Kim23]. 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 [Reference Lee, Kim and Kwak24]:

\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 [Reference Desbrun, Meyer, Schroder and Barr10, Reference Hwang, Ham, Lee, Lee, Kang and Kim17]. 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.