1 Introduction

1.1 Development of the population model

In the ecosystem, the functional response function can reflect the optimal feeding behavior of the predator with the maximum energy intake per unit time in order to achieve the maximum growth capacity of the population. Holling proposed three kinds of different functional response functions, later, a Holling IV functional response function was proposed to describe the interaction between zooplankton and phytoplankton. In addition, scholars have also intensively studied on predator–prey models with Beddington–DeAngelis functional response function [1], a ratio-dependent functional response function [2], the Ivelev functional response function [3] and the Crowley–Martin functional response function [4].

At the same time, more and more biological effects are applied to the predator–prey systems when studying the stability of the equilibrium, such as Allee effect [5], prey refuge effect [611], habitat complexity effect [1214], and harvesting effect [28]. Habitat usually refers to the area where organisms live, in which the organism can find food, shelter and breed. Homogeneous habitat means that the habitat has the same resource level. However, in fact, the living environment of natural habitats is spatial heterogeneous, and existing studies have shown that most of habitats have complexity due to heterogeneity [1517]. A large number of experimental studies have shown that the complexity of habitats reduces the meeting rate between predator and prey, thus reduces the predation rate of predators [1823]. Therefore, the effect of habitat complexity on the interaction between predator and prey cannot be ignored. However, the habitat complexity effect reduces the predation probability of prey, and the prey is not absolutely safe, while prey with refuge effect are absolutely safe.

In the natural environment, because of the limited resources, the spatial distribution of the population is heterogeneous, the biology will search for food everywhere in order to survive, then migration and diffusion will occur. Therefore, considering the heterogeneity of spatial distribution of the population, the corresponding reaction–diffusion system can be obtained. Different boundary conditions represent different biological significance in the study of predator–prey systems with diffusion term. For example, the homogeneous Neumann boundary condition means that neither prey nor predator can cross the boundary, the homogeneous Dirichlet boundary means that the number of both prey and predator is zero at the boundary, the homogeneous Robin boundary condition means that prey or predator can cross the boundary.

1.2 Establishment of the model

In [24], the author introduced habitat complexity into the ordinary differential equation system with Holling type functional response function and delay. The model is as follows:

$$ \textstyle\begin{cases} \dot{x}(t)=rx(1-\frac{x}{K})- \frac{c(1-\beta ){{x}^{\alpha }}y}{1+ch(1-\beta ){{x}^{\alpha }}}, \\ \dot{y}(t)= \frac{ec(1-\beta ){{x}^{\alpha }}(t-\tau )y(t-\tau )}{1+ch(1-\beta ){x}^{\alpha } (t-\tau )}-dy, \\ x ( \xi )=\phi ( \xi )>0,\qquad y ( \xi )=\psi ( \xi )>0,\quad \xi \in (-\tau ,0] , \end{cases} $$
(1.1)

where \(x ( t )\), \(y ( t )\), respectively, represents the density of prey population and predator population at time t, and the other parameters are all positive. The biological significance is expressed as follows: r is the intrinsic growth rate of prey population; K is the maximum environmental capacity of prey population; \(\frac{c(1-\beta ){{x}^{\alpha }}}{1+ch(1-\beta ){{x}^{\alpha }}}\) represents Holling type functional response function, in which, \(\alpha \geq 1\) and represents a kind of aggregation efficiency, when \(\alpha =1\), it becomes Holling II type functional response, when \(\alpha =2\), it becomes Holling III type functional response; c represents the attack rate of the predator on the prey; h is the handing time; e (\(0< e<1\)) is the conversion efficiency; d represents per capita death rate of predators; β (\(0<\beta <1\)) represents the intensity of the habitat complexity effect.

Considering that the habitat is heterogeneous, we introduce the diffusion term in the model (1.1), and obtain the reaction–diffusion system with homogeneous Neumann boundary conditions, the model is as follows:

$$ \textstyle\begin{cases} \frac{\partial u}{\partial t}={{d}_{1}}\Delta u+ru(1-\frac{u}{K})- \frac{c(1-\beta ){{u}^{\alpha }}v}{1+ch(1-\beta ){{u}^{\alpha }}} , \\ \frac{\partial v}{\partial t}={{d}_{2}}\Delta v+ \frac{ec(1-\beta ){{u}^{\alpha }}(x,t-\tau )v(x,t-\tau )}{1+ch(1-\beta ){{u}^{\alpha }}(x,t-\tau )}-dv(x,t) , \\ {{u}_{x}}(0,t)={{v}_{x}}(0,t)=0,\qquad {{u}_{x}}(l\pi ,t)={{v}_{x}}(l\pi ,t)=0,\quad t>0 , \\ u(x,0)={{u}_{0}}(x)\ge 0,\qquad v(x,0)={{v}_{0}}(x)\ge 0,\quad x\in \Omega =(0,l \pi ). \end{cases} $$
(1.2)

Next, we will study the stability of the equilibrium point of the system, give the existence conditions of Hopf bifurcation, and derive the properties of Hopf bifurcation using the central manifold theory and normal form method proposed by Hassar [25], Wu [26], Faria [27], including the direction of Hopf bifurcation and the stability of bifurcating periodic solutions.

1.3 Existence of the constant equilibria

In order to ensure the biological significance of the model (1.2), we make the following hypothesis:

(\(\mathbf{H_{0}}\)):

\(h < e/d\) and \(\alpha \geq 1 \).

Three equilibria of the system (1.2) can be obtained:

$$ {{P}_{0}}= ( 0,0 ),\qquad {{P}_{1}}= ( K,0 ),\qquad {{P}^{*}}=\bigl({{u}^{*}},{{v}^{*}}\bigr), $$

where

$$ {{u}^{*}}=\sqrt[\alpha ]{\frac{1}{1-\beta }} \sqrt[\alpha ]{\frac{d}{c(e-dh)}}, \qquad {{v}^{*}}= \frac{er {{u}^{*}}}{d}\biggl(1-\frac{{{u}^{*}}}{k}\biggr), $$

\({{u}^{*}}\) can be regarded as a function of habitat complexity effect β, let \({{u}^{*}}=u(\beta )\). For convenience, denote \(\beta ^{*}=1-\frac{d}{ck^{\alpha }(e-dh)}\) and make the following assumption:

(\(\mathbf{H_{1}}\)):

\(\beta <\beta ^{*} \).

Theorem 1.1

If (\(\mathbf{H_{0}}\)) and (\(\mathbf{H_{1}}\)) hold, then the system (1.2) has a unique positive equilibrium point \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\).

2 Stability analysis of diffusion system without time delay

When \(\tau =0\), the system (1.2) becomes

$$ \textstyle\begin{cases} \frac{\partial u}{\partial t}={{d}_{1}}\Delta u+ru(1-\frac{u}{K})- \frac{c(1-\beta ){{u}^{\alpha }}v}{1+ch(1-\beta ){{u}^{\alpha }}} , \\ \frac{\partial v}{\partial t}={{d}_{2}}\Delta v+ \frac{ec(1-\beta ){{u}^{\alpha }}v}{1+ch(1-\beta ){{u}^{\alpha }}}-dv , \\ {{u}_{x}}(0,t)={{v}_{x}}(0,t)=0,\qquad {{u}_{x}}(l\pi ,t)={{v}_{x}}(l\pi ,t)=0,\quad t>0 , \\ u(x,0)={{u}_{0}}(x)\ge 0,\qquad v(x,0)={{v}_{0}}(x)\ge 0,\quad x\in \Omega =(0,l \pi ). \end{cases} $$
(2.1)

Define the real-valued Sobolev space

$$ X:= \bigl\{ {{ ( u,v )}^{T}}| u,v\in {{H}^{2}} ( 0,l\pi ), ( {{u}_{x}},{{v}_{x}} )| _{x=0,l \pi }= ( 0,0 ) \bigr\} , $$

and let the complexification of X be

$$ {{X}_{c}}:=X\oplus iX= \{ {{x}_{1}}+i{{x}_{2}}| {{x}_{1}},{{x}_{2}} \in X \} . $$

Let

$$ U= ( u,v )\in {{H}^{2}} ( 0,l\pi ),\qquad D=\operatorname{diag} ( {{d}_{1}},{{d}_{2}} ),\qquad F ( \beta ,U )= ( {{F}_{1}},{{F}_{2}} ), $$

then the system (2.1) can be written in the form of an abstract equation

$$ \dot{U} ( t )=D\Delta U ( t )+F ( \beta ,U ). $$

We use \(J ( F )\) to represent the Jacobian matrix of F, then the linearization of the steady state system corresponding to the system (2.1) at \(( \beta ,0,0 )\) is

$$ L ( \beta )=D\frac{{{\partial }^{2}}}{\partial {{x}^{2}}}+J ( F )| _{U\equiv 0} . $$

2.1 Stability of positive equilibrium

Let the Jacobian matrix corresponding to system (2.1) at the positive equilibrium point \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\) be

( a 1 a 2 a 3 a 4 ) ,

where

$$ \begin{aligned} &{{a}_{1}}= \biggl( r- \frac{2r{{u}^{*}}}{K} \biggr)+\alpha \biggl( 1- \frac{dh}{e} \biggr) \biggl( r-\frac{r{{u}^{*}}}{K} \biggr), \\ &{{a}_{2}}=- \frac{c ( 1-\beta ){{{u}^{*}}^{\alpha }}}{1+ch ( 1-\beta ){{{u}^{*}}^{\alpha }}}=- \frac{d}{e}, \\ &{{a}_{3}}=\alpha r \biggl( 1-\frac{{u}^{*}}{K} \biggr) ( e-dh ),\qquad {{a}_{4}}=0. \end{aligned} $$

We assume

L n = ( a 1 d 1 n 2 l 2 a 2 a 3 d 2 n 2 l 2 ) , E n = tr ( L n ) = a 1 + ( d 1 + d 2 ) n 2 l 2 , F n = | L n | = d 1 d 2 n 4 l 4 a 1 d 2 n 2 l 2 a 2 a 3 ,

then the characteristic equation of \({{L}_{n}}\) is

$$ {{\lambda }^{2}}+{{E}_{n}}\lambda +{{F}_{n}}=0, $$
(2.2)

the characteristic roots of Eq. (2.2) are

$$ \lambda _{1,2}^{ ( n )}= \frac{-{{E}_{n}}\pm \sqrt{{{E}_{n}}^{2}-4{{F}_{n}}}}{2},\quad n\in { \mathbb{N}_{0}}\triangleq \{0\}\cup \mathbb{N}. $$

Theorem 2.1

Suppose (\(\mathbf{H_{0}}\)) and (\(\mathbf{H_{1}}\)) are true, then we can obtain

  1. (i)

    If \(1-(1-{{\beta }^{*}}){{ ( 1+\frac{e}{e+\alpha (e-dh)} )}^{ \alpha }}<\beta <{{\beta }^{*}}\), then the system (2.1) is locally asymptotically stable at \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\).

  2. (ii)

    If \(\beta <1-(1-{{\beta }^{*}}){{ ( 1+\frac{e}{e+\alpha (e-dh)} )}^{\alpha }}\), then the system (2.1) is unstable at \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\).

2.2 Stability of boundary equilibria

At the equilibrium point \({{P}_{0}}= ( 0,0 )\), the corresponding characteristic roots are \(\lambda _{01}^{n}=r-{{d}_{1}}\frac{{{n}^{2}}}{{{l}^{2}}}\), \(\lambda _{02}^{n}=-d-{{d}_{1}} \frac{{{n}^{2}}}{{{l}^{2}}}<0\). At the equilibrium point \({{P}_{1}}= ( K,0 )\), the corresponding characteristic roots are \(\lambda _{11}^{n}=-r-{{d}_{1}}\frac{{{n}^{2}}}{{{l}^{2}}}<0\), \(\lambda _{12}^{n}= \frac{ec ( 1-\beta ){{K}^{\alpha }}}{1+ch ( 1-\beta ){{K}^{\alpha }}}-d-{{d}_{2}} \frac{{{n}^{2}}}{{{l}^{2}}}\).

Theorem 2.2

For the system (2.1), the following statements are true.

  1. (i)

    \({{P}_{0}}= ( 0,0 )\) is unstable.

  2. (ii)

    If \(\beta >{{\beta }^{*}}\), then \({{P}_{1}}= ( K,0 )\) is locally asymptotically stable; if \(\beta <{{\beta }^{*}}\), \({{P}_{1}}= ( K,0 )\) is unstable.

Theorem 2.3

Suppose (\(\mathbf{H_{0}}\)) and (\(\mathbf{H_{1}}\)) are true, for the system (2.1), if \(\beta >{{\beta }^{*}}\), then the equilibrium \({{P}_{1}}= ( K,0 )\) is globally asymptotically stable.

Proof

Under hypotheses (\(\mathbf{H_{0}}\)) and (\(\mathbf{H_{1}}\)), \(u\geq 0\), \(v \geq 0\). By the first equation of the system (2.1),

$$ \frac{\partial u}{\partial t}-{{d}_{1}}\Delta u=ru\biggl(1- \frac{u}{K}\biggr)- \frac{c(1-\beta ){{u}^{\alpha }}v}{1+ch(1-\beta ){{u}^{\alpha }}}\le ru\biggl(1- \frac{u}{K}\biggr). $$

Using the comparison principle, we know \(\lim_{t\to \infty } \max_{x\in [ 0,l\pi ]} u ( x,t )\le K \). From the second equation,

$$ \begin{aligned} \frac{\partial v}{\partial t}-{{d}_{2}} \Delta v&= \frac{ec(1-\beta ){{u}^{\alpha }}v}{1+ch(1-\beta ){{u}^{\alpha }}}-dv \\ &=v\biggl(\frac{ec(1-\beta )}{1/{{u}^{\alpha }}+ch(1-\beta )}-d\biggr) \\ &\le v\biggl(\frac{ec(1-\beta )}{1/{{K}^{\alpha }}+ch(1-\beta )}-d\biggr)< 0, \end{aligned} $$

we can obtain \(\frac{ec(1-\beta ){{u}^{\alpha }}v}{1+ch(1-\beta ){{u}^{\alpha }}}-dv<0\). Therefore, for any \(\varepsilon >0\), there exist \(T>0\), \(v ( x,t )\le \varepsilon \) so that

$$ \begin{aligned} \frac{\partial u}{\partial t}-{{d}_{1}}\Delta u&=ru\biggl(1-\frac{u}{K}\biggr)- \frac{c(1-\beta ){{u}^{\alpha }}v}{1+ch(1-\beta ){{u}^{\alpha }}} \\ &\ge u\biggl[r\biggl(1-\frac{u}{K}\biggr)-c(1-\beta ){{u}^{\alpha -1}}\varepsilon \biggr] \\ &\ge u\biggl[r\biggl(1-\frac{u}{K}\biggr)-c(1-\beta ){{K}^{\alpha -1}}\varepsilon \biggr]. \end{aligned} $$

Applying the comparison principle again, we have \(u ( x,t )\ge K ( 1- \frac{c ( 1-\beta ){{K}^{\alpha }}\varepsilon }{r} )\), \(t>T\), \(x\in [ 0,l\pi ]\), so \(\lim_{t\to \infty } \max_{x\in [ 0,l\pi ]} u ( x,t )=K\). That is, the equilibrium point \({{P}_{1}}= ( K,0 )\) is globally asymptotically stable. □

3 Hopf bifurcation property analysis of time-delay system

When \(\tau \ne 0\), we will study the time-delay effect on dynamic properties of the diffusion system (1.2).

3.1 Existence of Hopf bifurcation induced by time delay

Assuming that (\(\mathbf{H_{0}}\)) and (\(\mathbf{H_{1}}\)) are true, the system (1.2) has a unique positive equilibrium \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\). For convenience, we make transformations \(\hat{u}=u-{{u}^{*}}\), \(\hat{v}=v-{{v}^{*}}\) to move \({{P}^{*}}\) to the origin. We still use u, v to represent û, , and then the system (1.2) becomes

$$ \textstyle\begin{cases} \frac{\partial u}{\partial t}={{d}_{1}}\Delta u+r(u+{u}^{*})(1- \frac{u+{u}^{*}}{K})- \frac{c(1-\beta ){{ ( u+{u}^{*} )}^{\alpha }}(v+{v}^{*})}{1+ch(1-\beta ){{ ( u+{u}^{*} )}^{\alpha }}}, \\ \frac{\partial v}{\partial t}={{d}_{2}}\Delta v+ \frac{ec(1-\beta ){{[u(x,t-\tau )+{u}^{*}]}^{\alpha }}[v(x,t-\tau )+{v}^{*}]}{1+ch(1-\beta ){{[u(x,t-\tau )+{u}^{*}]}^{\alpha }}}-d(v+{v}^{*}), \\ {{u}_{x}}(0,t)={{v}_{x}}(0,t)=0,\qquad {{u}_{x}}(l\pi ,t)={{v}_{x}}(l\pi ,t)=0,\quad t>0, \\ u(x,0)={{u}_{0}}(x)\ge 0,\qquad v(x,0)={{v}_{0}}(x)\ge 0,\quad x\in \Omega =(0,l \pi ). \end{cases} $$
(3.1)

Let

$$ {{u}_{1}} ( t )=u ( \cdot ,t ),\qquad {{u}_{2}} ( t )=v ( \cdot ,t ),\qquad U={{ ( {{u}_{1}},{{u}_{2}} )}^{T}},\qquad X=C \bigl( [ 0,l\pi ],{{\mathbb{R}}^{2}} \bigr), $$

in phase space \({{\mathbb{C}}_{\tau }}=C ( [ -\tau ,0 ],X )\), (3.1) can be abstracted as

$$ \dot{U} ( t )=D\Delta U ( t )+L ( {{U}_{t}} )+F ( {{U}_{t}} ), $$
(3.2)

where \(\phi ={{ ( {{\phi }_{1}},{{\phi }_{2}} )}^{T}}\), \(D=\operatorname{diag}({{d}_{1}},{{d}_{2}}). L:{{\mathbb{C}}_{\tau }}\to X\), \(F:{{\mathbb{C}}_{\tau }}\to X\) are defined as follows:

L(ϕ)= ( a 1 a 2 0 d ) ( ϕ 1 ( 0 ) ϕ 2 ( 0 ) ) + ( 0 0 c 1 d ) ( ϕ 1 ( τ ) ϕ 2 ( τ ) ) ,F(ϕ)= ( F 1 ( ϕ ) F 2 ( ϕ ) ) ,

with

$$\begin{aligned}& \begin{aligned} {{F}_{1}}(\phi )={}&r \bigl( {{\phi }_{1}}(0)+{{u}^{*}} \bigr) \biggl( 1- \frac{{{\phi }_{1}}(0)+{{u}^{*}}}{k} \biggr)- \frac{c(1-\beta ){{ ( {{\phi }_{1}}(0)+{{u}^{*}} )}^{\alpha }}({{\phi }_{2}}(0)+{{v}^{*}})}{1+ch(1-\beta ){{ ( {{\phi }_{1}}(0)+{{u}^{*}} )}^{\alpha }}} \\ &{} -{{a}_{1}} {{\phi }_{1}}(0)-{{a}_{2}} {{ \phi }_{2}}(0), \end{aligned} \\& \begin{aligned} {{F}_{2}}(\phi )={}& \frac{ec(1-\beta ){{ ( {{\phi }_{1}}(-\tau )+{{u}^{*}} )}^{\alpha }}({{\phi }_{2}}(-\tau )+{{v}^{*}})}{1+ch(1-\beta ){{ ( {{\phi }_{1}}(-\tau )+{{u}^{*}} )}^{\alpha }}}-d\bigl({{ \phi }_{2}}(0)+{{v}^{*}}\bigr)+d{{\phi }_{2}}(0) \\ &{} -{{c}_{1}} {{\phi }_{1}}(-\tau )-d{{\phi }_{2}}(-\tau ), \end{aligned} \\& {{a}_{1}}= \biggl( r-\frac{2r{{u}^{*}}}{K} \biggr)+\alpha \biggl( 1- \frac{dh}{e} \biggr) \biggl( r-\frac{r{{u}^{*}}}{K} \biggr),\qquad {{a}_{2}}=- \frac{d}{e}, \\& {{c}_{1}}= \alpha r \biggl( 1-\frac{{u}^{*}}{K} \biggr) ( e-dh ). \end{aligned}$$

The linearized equation of (3.1) at the origin is

$$ \dot{U} ( t )=D\Delta U ( t )+{L} ( {{U}_{t}} ), $$
(3.3)

where

L( U t )= L 1 U+ L 2 U t , L 1 = ( a 1 a 2 0 d ) , L 2 = ( 0 0 c 1 d ) .

It is well known that the eigenvalue of \(-{\varphi }''=\mu \varphi \), \(x\in ( 0,l\pi ) \), \({\varphi }' ( 0 )={\varphi }' ( l\pi )=0 \) is \({{\mu }_{n}}={{{n}^{2}}}/{{{l}^{2}},n\in {{\mathbb{N}}_{0}}}\), the characteristic function is \({{\varphi }_{n}}=\cos \frac{n\pi }{l}\), \(n\in {{ \mathbb{N}}_{0}}\), λ is the eigenvalue of (3.3). Substituting y= n = 0 ( y 1 n y 2 n ) cos n π l into \(\lambda y-d\Delta y-L({{e}^{\lambda }}y)=0\), we can obtain

( a 1 d 1 μ n a 2 c 1 e λ τ d + d e λ τ d 2 μ n ) ( y 1 n y 2 n ) =λ ( y 1 n y 2 n ) ,n N 0 ,

the corresponding characteristic equation is

$$ \det \bigl( \lambda I+{{\mu }_{n}}D-{{L}_{1}}-{{L}_{2}} {{e}^{- \lambda \tau }} \bigr)=0, \quad n\in {{\mathbb{N}}_{0}}. $$

It is equivalent to

$$ {{f}_{n}} ( \lambda ,\tau )={{\lambda }^{2}}+{{A}_{n}} \lambda +{{B}_{n}}+{{C}_{n}} {{e}^{-\lambda \tau }}=0, $$
(3.4)

where

$$ \begin{aligned} &{{A}_{n}}= ( {{d}_{1}}+{{d}_{2}} ){{\mu }_{n}}-{{a}_{1}}+d, \\ &{{B}_{n}}={{d}_{1}} {{d}_{2}} {{\mu }_{n}}^{2}- ( {{a}_{1}} {{d}_{2}}-d{{d}_{1}} ){{\mu }_{n}}-{{a}_{1}}d, \\ &{{C}_{n}}={{a}_{1}}d-{{a}_{2}} {{c}_{1}}-d{{d}_{1}} {{\mu }_{n}}-d \lambda. \end{aligned} $$

Make the assumptions:

(\(\mathbf{H_{2}}\)):

\({{a}_{1}}<0\),

(\(\mathbf{H_{3}}\)):

\({{a}_{2}}{{c}_{1}}-2{{a}_{1}}d<0\),

(\(\mathbf{H_{4}}\)):

\({{a}_{2}}{{c}_{1}}-2{{a}_{1}}d\geq 0\).

We have the following lemmas.

Lemma 3.1

If (\(\mathbf{H_{0}}\))(\(\mathbf{H_{2}}\)) are true, the following conclusions can be drawn.

  1. (i)

    When \(\tau =0\), all the characteristic roots of Eq. (3.4) have negative real parts, and the system (3.1) is locally asymptotically stable at the positive equilibrium point \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\);

  2. (ii)

    \(\lambda =0\) is not the root of Eq. (3.4).

Lemma 3.2

Assume (\(\mathbf{H_{0}}\))(\(\mathbf{H_{2}}\)) are true, we have the following conclusions.

  1. (i)

    If (\(\mathbf{H_{3}}\)) holds, when \(n>{{N}_{1}}\), Eq. (3.4) has no pure imaginary roots. When \(0\leq n \leq {{N}_{1}}\), Eq. (3.4) has a pair of pure imaginary roots \(\pm i\omega _{n}^{+}\) at \(\tau =\tau _{n}^{j,+}\).

  2. (ii)

    If (\(\mathbf{H_{4}}\)) holds, Eq. (3.4) has no pure imaginary roots.

We have

$$ {{N}_{1}}= \textstyle\begin{cases} [\hat{N} ],&\hat{N}\notin \mathbb{N}, \\ [\hat{N} ]-1,&\hat{N}\in \mathbb{N}, \end{cases}\displaystyle \quad [\hat{N}] \text{ represents rounding } \hat{N}, $$

where

$$ \begin{aligned} &\hat{N}=l\sqrt{\frac{1}{2{{d}_{1}}{{d}_{2}}} \bigl[({{a}_{1}} {{d}_{2}}-2d{{d}_{1}})+ \sqrt{{{({{a}_{1}} {{d}_{2}}-2d{{d}_{1}})}^{2}}-4{{d}_{1}} {{d}_{2}}({{a}_{2}} {{c}_{1}}-2{{a}_{1}}d)} \bigr]}, \\ &\tau _{n}^{j,+}=\frac{1}{\omega _{n}^{+}}\arccos \frac{ ( {{D}_{n}}+d{{A}_{n}} ){{(\omega _{n}^{+})}^{2}}-{{D}_{n}}{{B}_{n}}}{{{D}_{n}}^{2}+{{d}^{2}}{{(\omega _{n}^{+})}^{2}}}+ \frac{2j\pi }{\omega _{n}^{+}},\quad j \in {{\mathbb{N}}_{0}}. \end{aligned} $$

Proof

We seek the critical value of τ which makes for the characteristic equation (3.4) exist a pair of pure imaginary roots. Let \(\lambda =i\omega\) (\(\omega >0 \)) be the root of (3.4), and for some \(n\in {{\mathbb{N}}_{0}}\), ω satisfies \(-{{\omega }^{2}}+i\omega {{A}_{n}}+{{B}_{n}}+[({{a}_{1}}-{{d}_{1}}{{ \mu }_{n}}-i\omega )d-{{a}_{2}}{{c}_{1}}] ( \cos \omega \tau -i \sin \omega \tau )=0\). Separate the real part and the imaginary part,

$$ \textstyle\begin{cases} ({{a}_{1}}d-d{{d}_{1}}{{\mu }_{n}}-{{a}_{2}}{{c}_{1}})\cos \omega \tau -d\omega \sin \omega \tau ={{\omega }^{2}}-{{B}_{n}}, \\ ({{a}_{1}}d-d{{d}_{1}}{{\mu }_{n}}-{{a}_{2}}{{c}_{1}})\sin \omega \tau +d\omega \cos \omega \tau ={{A}_{n}}\omega . \end{cases} $$
(3.5)

Let \({{D}_{n}}={{a}_{1}}d-d{{d}_{1}}{{\mu }_{n}}-{{a}_{2}}{{c}_{1}}\), then we have

$$ {{\omega }^{4}}+\bigl({{A}_{n}}^{2}-2{{B}_{n}}-{{d}^{2}} \bigr){{\omega }^{2}}+{{B}_{n}}^{2}-{{D}_{n}}^{2}=0. $$
(3.6)

Let \(z={{\omega }^{2}}\), then (3.6) can be rewritten as

$$ {{z}^{2}}+\bigl({{A}_{n}}^{2}-2{{B}_{n}}-{{d}^{2}} \bigr)z+{{B}_{n}}^{2}-{{D}_{n}}^{2}=0, $$
(3.7)

where

$$ \begin{aligned} &{{A}_{n}}^{2}-2{{B}_{n}}-{{d}^{2}}= \bigl({{d}_{1}}^{2}+{{d}_{2}}^{2} \bigr){{ \mu }_{n}}^{2}-2{{a}_{1}} {{d}_{1}} {{\mu }_{n}}^{2}+{{a}_{1}}^{2}+2d{{d}_{2}} {{ \mu }_{n}}>0, \\ &{{B}_{n}}+{{D}_{n}}={{d}_{1}} {{d}_{2}} {{\mu }_{n}}^{2}-{{a}_{1}} {{d}_{2}} {{ \mu }_{n}}-{{a}_{2}} {{c}_{1}}>0, \\ &{{B}_{n}}-{{D}_{n}}={{d}_{1}} {{d}_{2}} {{\mu }_{n}}^{2}-({{a}_{1}} {{d}_{2}}-2d{{d}_{1}}){{ \mu }_{n}}-2{{a}_{1}}d+{{a}_{2}} {{c}_{1}}, \end{aligned} $$

under hypothesis (\(\mathbf{H_{3}}\)), when \(0\le n\le {{N}_{1}}\), \({{B}_{n}}-{{D}_{n}}<0\), \({{B}_{n}}^{2}-{{D}_{n}}^{2}< 0\); when \(n>{{N}_{1}}\), \({{B}_{n}}-{{D}_{n}}\ge 0\), \({{B}_{n}}^{2}-{{D}_{n}}^{2} \geq 0\). If hypothesis (\(\mathbf{H_{4}}\)) holds, \({{B}_{n}}-{{D}_{n}}>0\), \({{B}_{n}}^{2}-{{D}_{n}}^{2}> 0\) for \(n\in {{\mathbb{N}}_{0}}\).

In summary, the conclusions are true. The roots of Eq. (3.7) are

$$ {{z}^{\pm }}= \frac{-({{A}_{n}}^{2}-2{{B}_{n}}-{{d}^{2}})\pm \sqrt{{{({{A}_{n}}^{2}-2{{B}_{n}}-{{d}^{2}})}^{2}}-4({{B}_{n}}^{2}-{{D}_{n}}^{2})}}{2}. $$

Thus Eq. (3.7) has one positive root \(z_{n}^{+}\). If all parameter values of the system (3.1) are given, the roots of Eq. (3.6) can be easily calculated, and \(\omega _{n}^{+}=\sqrt{z_{n}^{+}}\). □

Lemma 3.3

If hypothesis (\(\mathbf{H_{2}}\)) is true, then \({\alpha }' ( \tau _{n}^{j,+} )=\frac{d\lambda }{d\tau } | _{\tau =\tau _{n}^{j,+}}> 0\).

Proof

Taking the derivative of Eq. (3.4) with respect to τ, we can obtain

$$ {{ \biggl( \frac{d\lambda }{d\tau } \biggr)}^{-1}}= \frac{2\lambda +{{A}_{n}}-d{{e}^{-\lambda \tau }}}{\lambda {{C}_{n}}{{e}^{-\lambda \tau }}} - \frac{\tau }{\lambda }, $$

then

$$ \begin{aligned} &\operatorname{sign} \biggl\{ \operatorname{Re} {{ \biggl( \frac{d\lambda }{d\tau } \bigg| _{\tau =\tau _{n}^{j,+}} \biggr)}^{-1}} \biggr\} \\ &\quad =\operatorname{sign} {{ \biggl\{ \operatorname{Re} \biggl( \frac{2\lambda +{{A}_{n}}-d{{e}^{-\lambda \tau }}}{\lambda {{C}_{n}}{{e}^{-\lambda \tau }}}- \frac{\tau }{\lambda } \biggr) \biggr\} }_{{\tau =\tau _{n}^{j,+}}}} \\ &\quad =\operatorname{sign} \biggl\{ \frac{2{{\omega }^{2}}({{\omega }^{2}}-{{B}_{n}})+{{A}_{n}}^{2}{{\omega }^{2}}-{{d}^{2}}{{\omega }^{2}}}{{{d}^{2}}{{\omega }^{4}}+{{D}_{n}}^{2}{{\omega }^{4}}} \biggr\} \\ &\quad =\operatorname{sign} \biggl\{ \frac{2{{\omega }^{2}}-2{{B}_{n}}+{{A}_{n}}^{2}-{{d}^{2}}}{{{d}^{2}}{{\omega }^{2}}+{{D}_{n}}^{2}} \biggr\} \\ &\quad =\operatorname{sign} \biggl\{ \frac{\sqrt{{{({{A}_{n}}^{2}-2{{B}_{n}}-d)}^{2}}-4({{B}_{n}}^{2}-{{D}_{n}}^{2})}}{{{d}^{2}}{{\omega }^{2}}+{{D}_{n}}^{2}} \biggr\} >0. \end{aligned} $$

 □

Let \(\lambda ( \tau )=\alpha ( \tau )+i\beta ( \tau )\) be the root of Eq. (3.4) when τ is sufficiently close to \(\tau _{n}^{j,+}\) \((0\leq n\leq {{N}_{1}})\), which satisfies \(\alpha ( \tau _{n}^{j,+} )=0\), \(\beta ( \tau _{n}^{j,+} )=\omega _{n}^{+}\), \(j\in {{\mathbb{N}}_{0}}\). According to the Rouche theorem, as τ changes from a value less than \(\tau _{n}^{j,+}\) to a value greater than \(\tau _{n}^{j,+}\), the characteristic root of (3.4) transverses the imaginary axis. Therefore, when \(\tau =\tau _{n}^{j,+}\), the system (3.1) satisfies the condition for a Hopf bifurcation to occur.

Obviously, \(\tau _{n}^{0}=\min_{j\in \mathbb{N}_{0}} \{ \tau _{n}^{j,+} \} \), denote \(\tau _{*}^{0}=\min_{0\leq n\leq {{N}_{1}}} \{ \tau _{n}^{0} \} \), we have the following theorem.

Theorem 3.1

Assume that (\(\mathbf{H_{0}}\))(\(\mathbf{H_{2}}\)) are true, if (\(\mathbf{H_{3}}\)) (\(0\leq n \leq {{N}_{1}}\)) is also satisfied, we have the following conclusions for the system (3.1).

  1. (i)

    When \(\tau \in [0,\tau _{*}^{0})\), the equilibrium \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\) is locally asymptotically stable.

  2. (ii)

    When \(\tau >\tau _{*}^{0}\), the equilibrium \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\) is unstable.

  3. (iii)

    When \(\tau =\tau _{0}^{j,+}\), \(j\in {{\mathbb{N}}_{0}}\), Hopf bifurcation occurs at the equilibrium \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\) and the bifurcating periodic solutions are homogeneous; When \(\tau \in \{\tau _{n}^{j,+}:\tau _{n}^{j,+}\ne \tau _{m}^{i,+},m\ne n,0 \le m,n\leq {{N}_{1}},j,i\in {\mathbb{N}_{0}}\}/\{\tau _{0}^{k,+} | k\in {\mathbb{N}_{0}} \}\), the system undergoes Hopf bifurcation at \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\) and the bifurcating periodic solutions are inhomogeneous.

3.2 Direction and periodic solutions of Hopf bifurcation

In this section, we shall study the direction of Hopf bifurcation and the stability of bifurcating periodic solutions. For fixed \(j\in \mathbb{N}_{0}\) and \(0\le n \le N_{1}\), we denote \(\tilde{\tau }=\tau _{n}^{j,+}\), \(\omega_{n}=\omega_{n}^{+}\), and let \(\tilde{u}(x,t)=u(x,\tau t)-{{u}^{*}}\) and \(\tilde{v}(x,t)=v(x,\tau t)-{{v}^{*}}\). For convenience, we drop the tilde, then the system (1.2) can be transformed into

$$ \textstyle\begin{cases} \frac{\partial {u}}{\partial {t}}=\tau [d_{1} \Delta u +r(u+{{u}^{*}})(1- \frac{u+{{u}^{*}}}{k})- \frac{c(1-\beta )(u+{{u}^{*}})^{\alpha }(v+{{v}^{*}})}{1+ ch(1-\beta )(u+{{u}^{*}})^{\alpha }}], \\ \frac{\partial v}{\partial t}=\tau [d_{2} \Delta v+ \frac{ec(1-\beta )[u(t-1)+{{u}^{*}}]^{\alpha }[v(t-1)+{{v}^{*}}]}{1+ ch(1-\beta )[u(t-1)+{{u}^{*}}]^{\alpha }}-d(v+{{v}^{*}})], \end{cases} $$
(3.8)

for \(x \in (0, l\pi )\), and \(t>0\). Let

$$ \tau =\tilde{\tau }+\mu ,\qquad u_{1}(t)=u(\cdot ,t),\qquad u_{2}(t)=v(\cdot ,t)\quad \mbox{and} \quad U=(u_{1},u_{2})^{T}. $$

Then (3.8) can be rewritten in an abstract form in the phase space \({{\mathbb{C}}_{\tau }}:=C([-1,0],X)\)

$$ \frac{dU(t)}{dt}=\tilde{\tau }D\Delta U(t)+L_{\tilde{\tau }}(U_{t})+F(U_{t}, \mu ), $$
(3.9)

in which

$$\begin{aligned}& \begin{aligned} L_{\mu }(\phi )=\mu \begin{pmatrix} a_{1} \phi _{1}(0)+a_{2} \phi _{2}(0) \\ -d\phi _{2}(0)+c_{1}\phi _{1}(-1)+d\phi _{2}(-1) \end{pmatrix} \end{aligned} , \end{aligned}$$
(3.10)
$$\begin{aligned}& F(\phi ,\mu )=\mu D\Delta \phi +L_{\mu }(\phi )+f( \phi ,\mu ), \end{aligned}$$
(3.11)

with

$$\begin{aligned}& f(\phi ,\mu )=(\tilde{\tau }+\mu )\bigl(F_{1}(\phi ,\mu ),F_{2}(\phi ,\mu )\bigr)^{T}, \\& \begin{aligned} F_{1}(\phi ,\mu )&=r\bigl(\phi _{1}(0)+{{u}^{*}}\bigr)\biggl(1- \frac{\phi _{1}(0)+{{u}^{*}}}{k}\biggr)- \frac{c(1-\beta )(\phi _{1}(0)+{{u}^{*}})^{\alpha }(\phi _{2}(0)+{{v}^{*}})}{1+ ch(1-\beta )(\phi _{1}(0)+{{u}^{*}})^{\alpha }} \\ &\quad {}-a_{1}\phi _{1}(0)-a_{2}\phi _{2}(0), \end{aligned} \\& \begin{aligned} F_{2}(\phi ,\mu )&= \frac{ec(1-\beta )(\phi _{1}(-1)+{{u}^{*}})^{\alpha }(\phi _{2}(-1)+{{v}^{*}})}{1+ ch(1-\beta )(\phi _{1}(-1)+{{u}^{*}})^{\alpha }}-d\bigl( \phi _{2}(0)+{{v}^{*}}\bigr) +d\phi _{2}(0) \\ &\quad {}-c_{1}\phi _{1}(-1)-d\phi _{2}(-1), \end{aligned} \end{aligned}$$

respectively, for \(\phi =(\phi _{1}, \phi _{2})^{T} \in {{\mathbb{C}}_{\tau }}\). Consider the linearized equation

$$ \frac{dU(t)}{dt}=\tilde{\tau } D \Delta U(t)+L_{\tilde{\tau }}(U_{t}). $$
(3.12)

According to the results in Sect. 3.1, \(\pm i \omega _{n} \) are eigenvalues of the system (3.12) and the linearized functional differential equation is

$$ \frac{dz(t)}{dt}=-\tilde{\tau } D\frac{n^{2}}{l^{2}} z(t)+L_{ \tilde{\tau }}(z_{t}). $$
(3.13)

By the Riesz representation theorem, there exists a \(2\times 2\) matrix function \(\eta ^{n}(\sigma , \tilde{\tau })\), \(-1\le \sigma \le 0\), whose elements are bounded variation functions such that

$$ -\tilde{\tau } D\frac{n^{2}}{l^{2}} \phi (0)+L_{\tilde{\tau }}(\phi )= \int _{-1}^{0}d \eta ^{n}(\sigma , \tau ) \phi (\sigma ) $$

for \(\phi \in C([-1,0],{{\mathbb{R}}^{2}})\).

We select

$$ \eta ^{n}(\sigma , \tau )= \textstyle\begin{cases} \tau E, & \sigma =0, \\ 0, & \sigma \in (-1,0), \\ - \tau F, & \sigma =-1, \end{cases} $$
(3.14)

in which

$$ E= \begin{pmatrix} a_{1}-d_{1}\frac{n^{2}}{l^{2}} & a_{2} \\ 0 & -d-d_{2}\frac{n^{2}}{l^{2}} \end{pmatrix}, \qquad F= \begin{pmatrix} 0 & 0 \\ c_{1} & d \end{pmatrix} . $$
(3.15)

Denote by \(A(\tilde{\tau })\) the infinitesimal generators of semigroup included by the solutions of Eq. (3.13) and let \(A^{*}\) be the formal adjoint of \(A(\tilde{\tau })\) under the bilinear paring

$$ \begin{aligned} (\psi ,\phi )&=\psi (0) \phi (0)- \int _{-1}^{0} \int _{\xi =0}^{ \sigma } \psi (\xi -\sigma ) \,d\eta ^{n}(\sigma , \tilde{\tau }) \phi ( \xi ) \,d\xi \\ &=\psi (0) \phi (0)+ \tilde{\tau } \int _{-1}^{0} \psi (\xi +1) F \phi (\xi ) \,d\xi , \end{aligned} $$
(3.16)

for \(\phi , \psi \in C([-1,0],{{\mathbb{R}}^{2}})\). \(A(\tilde{\tau })\) and \(A^{*}\) both have a pair of simple purely imaginary eigenvalues \(\pm i \omega _{n} \tilde{\tau }\). Let P and \(P^{*}\) be the center subspace, i.e., the generalized eigenspace of \(A(\tilde{\tau })\) and \(A^{*}\) connected with \(\Lambda _{n}\), respectively. Then \(P^{*}\) is the adjoint space of P, and \(\operatorname{dim}P=\operatorname{dim}P^{*}=2\).

We can verify that \(p_{1}(\theta )=(1,\xi )^{T}e^{i\omega _{n} \tilde{\tau } \sigma } \) (\(\sigma \in [-1,0]\)), \(p_{2}(\sigma )=\overline{p_{1}(\sigma )} \) is a basis of \(A(\tilde{\tau })\) with \(\Lambda _{n}\) and \(q_{1}(r)=(1,\eta )e^{-i\omega _{n} \tilde{\tau } r} \) (\(r \in [0,1]\)), \(q_{2}(r)= \overline{q_{1}(r)}\) is a basis of \(A^{*}\) with \(\Lambda _{n}\), where

$$ \xi =\frac{i\omega _{n}-a_{1}+d_{1}\frac{n^{2}}{l^{2}}}{a_{2}},\qquad \eta = \frac{-i \omega _{n} -a_{1}+d_{1}\frac{n^{2}}{l^{2}}}{c_{1} e^{i \omega _{n} \tilde{\tau } }}. $$

Let \(\Phi =(\Phi _{1},\Phi _{2})\) and \(\Psi ^{*}=(\Psi ^{*}_{1},\Psi ^{*}_{2})^{T}\) with

$$ \begin{aligned} &\Phi _{1}(\sigma )=\frac{p_{1}(\sigma )+p_{2}(\sigma )}{2} = \begin{pmatrix} \cos (\omega _{n} \tilde{\tau } \sigma ) \\ \frac{-a_{1}+d_{1}\frac{n^{2}}{l^{2}}}{a_{2}} \cos \omega _{n} \tilde{\tau } \sigma - \frac{\omega _{n}}{a_{2}} \sin \omega _{n} \tilde{\tau } \sigma \end{pmatrix}, \\ &\Phi _{2}(\sigma )=\frac{p_{1}(\sigma )-p_{2}(\sigma )}{2i}= \begin{pmatrix} \sin (\omega _{n} \tilde{\tau } \theta ) \\ \frac{\omega _{n}}{a_{2}} \cos \omega _{n} \tilde{\tau } \sigma + \frac{-a_{1}+d_{1}\frac{n^{2}}{l^{2}}}{a_{2}} \sin \omega _{n} \tilde{\tau } \sigma \end{pmatrix}, \\ &\Psi ^{*} _{1}(r)=\frac{q_{1}(r)+q_{2}(r)}{2}= \begin{pmatrix} \cos (\omega _{n} \tilde{\tau } r ) \\ \frac{-a_{1}+d_{1}\frac{n^{2}}{l^{2}}}{c_{1}} \cos \omega _{n} \tilde{\tau }(r+1)-\frac{\omega _{n}}{c_{1}} \sin \omega _{n} \tilde{\tau } (r+1) \end{pmatrix}, \\ &\Psi ^{*} _{2}(r)=\frac{q_{1}(r)-q_{2}(r)}{2i}= \begin{pmatrix} -\sin (\omega _{n} \tilde{\tau } r ) \\ -\frac{\omega _{n}}{c_{1}} \cos \omega _{n} \tilde{\tau } (r+1)- \frac{-a_{1}+d_{1}\frac{n^{2}}{l^{2}}}{c_{1}} \sin \omega _{n} \tilde{\tau }(r+1) \end{pmatrix}, \end{aligned} $$

for \(\theta \in [-1,0]\), and \(r \in [0,1]\). According to (3.16), we can calculate

$$ D^{*}_{1}:=\bigl(\Psi ^{*} _{1},\Phi _{1}\bigr), \qquad D^{*}_{2}:= \bigl(\Psi ^{*} _{1}, \Phi _{2}\bigr),\qquad D^{*}_{3}:=\bigl(\Psi ^{*} _{2}, \Phi _{1}\bigr),\qquad D^{*}_{4}:=\bigl(\Psi ^{*} _{2},\Phi _{2}\bigr). $$

Define ( Ψ ,Φ)=( Ψ j , Φ k )= ( D 1 D 2 D 3 D 4 ) and construct a basis Ψ of \(P^{*}\) by

$$ \Psi =(\Psi _{1},\Psi _{2})^{T}=\bigl(\Psi ^{*},\Phi \bigr)^{-1}\Psi ^{*}, $$

then \((\Psi ,\Phi )=I_{2}\).

Furthermore, we define \(f_{n}:=(\beta ^{1}_{n},\beta ^{2}_{n})\) with

$$ \beta ^{1}_{n}= \begin{pmatrix} \cos \frac{n}{l}x \\ 0 \end{pmatrix},\qquad \beta ^{2}_{n}= \begin{pmatrix} 0 \\ \cos \frac{n}{l}x \end{pmatrix}, $$

and

$$ c\cdot f_{n}=c_{1} \beta ^{1}_{n}+c_{2} \beta ^{2}_{n}, \quad \mbox{for } c=(c_{1},c_{2})^{T} \in {{\mathbb{C}}_{\tau }}. $$

Thus the center subspace of the linear equation (3.12) is given by \(P_{CN} {{\mathbb{C}}_{\tau }} \oplus P_{S} {{\mathbb{C}}_{\tau }}\) and \(P_{S} {{\mathbb{C}}_{\tau }}\) is the complement subspace of \(P_{CN} {{\mathbb{C}}_{\tau }}\) in \({{\mathbb{C}}_{\tau }}\),

$$ < u,v>:=\frac{1}{l\pi } \int _{0}^{l\pi }u_{1} \overline{v_{1}}\,dx+ \frac{1}{l\pi } \int _{0}^{l\pi }u_{2} \overline{v_{2}}\,dx $$

for \(u=(u_{1},u_{2})\), \(v=(v_{1},v_{2})\), \(u,v\in X\) and \(\langle \phi ,f_{0}\rangle=(\langle \phi ,f^{1}_{0}\rangle, \langle\phi ,f^{2}_{0}\rangle)^{T}\).

Let \(A_{\tilde{\tau }}\) be the infinitesimal generator of an analytic semigroup induced by the linear system (3.12), and Eq. (3.8) can be rewritten as

$$ \frac{dU(t)}{dt}=A_{\tilde{\tau }} U_{t}+R(U_{t},\mu ), $$
(3.17)

where

$$ R(U_{t},\mu )= \textstyle\begin{cases} 0, & \theta \in [-1,0), \\ F(U_{t},\mu ), & \theta =0. \end{cases} $$
(3.18)

By the decomposition of \({{\mathbb{C}}_{\tau }}\), the solution above can be written as

$$ U_{t}=\Phi \begin{pmatrix} x_{1} \\ x_{2} \end{pmatrix}f_{n}+h(x_{1},x_{2}, \mu ), $$

where

$$ \begin{pmatrix} x_{1} \\ x_{2} \end{pmatrix}=\bigl(\Psi ,\langle U_{t},f_{n}\rangle \bigr) $$

and

$$ h(x_{1},x_{2},\mu ) \in P_{S} {{ \mathbb{C}}_{\tau }}, \qquad h(0,0,0)=0, \qquad Dh(0,0,0)=0. $$

Specially, the solution of (3.9) on the center manifold is given by

$$ U_{t}=\Phi \begin{pmatrix} x_{1}(t) \\ x_{2}(t) \end{pmatrix}f_{n}+h(x_{1},x_{2},0). $$
(3.19)

Let \(z=x_{1}-i x_{2}\), and because \(p_{1}=\Phi _{1}+i\Phi _{2}\), we have

$$ \Phi \begin{pmatrix} x_{1} \\ x_{2} \end{pmatrix}f_{n}= (\Phi _{1},\Phi _{2}) \begin{pmatrix} \frac{z+\overline{z}}{2} \\ \frac{i(z-\overline{z})}{2} \end{pmatrix}f_{n} = \frac{1}{2}(p_{1}z+ \overline{p_{1}} \overline{z})f_{n} $$

and

$$ h(x_{1},x_{2},0)=h\biggl(\frac{z+\overline{z}}{2}, \frac{i(z-\overline{z})}{2},0\biggr). $$

Therefore, Eq. (3.19) can be transformed into

$$ \begin{aligned} U_{t}&= \frac{1}{2}(p_{1}z+ \overline{p_{1}} \overline{z})f_{n}+ h\biggl( \frac{z+\overline{z}}{2}, \frac{i(z-\overline{z})}{2},0\biggr) \\ &=\frac{1}{2}(p_{1}z+ \overline{p_{1}} \overline{z})f_{n}+ W(z, \overline{z}), \end{aligned} $$
(3.20)

in which

$$ W(z,\overline{z})=h\biggl(\frac{z+\overline{z}}{2}, \frac{i(z-\overline{z})}{2},0\biggr), $$

z satisfies

$$ \dot{z}=i\omega _{n} \tilde{\tau } z +g(z,\overline{z}), $$
(3.21)

where

$$ g(z,\overline{z})=\bigl(\Psi _{1}(0)-i\Psi _{2}(0)\bigr)\bigl\langle F(U_{t},0),f_{n}\bigr\rangle . $$
(3.22)

Let

$$\begin{aligned}& W(z,\overline{z})=W_{20} \frac{z^{2}}{2}+W_{11} z \overline{z} +W_{02} \frac{\overline{z}^{2}}{2}+\cdots , \end{aligned}$$
(3.23)
$$\begin{aligned}& g(z,\overline{z})=g_{20} \frac{z^{2}}{2}+g_{11} z \overline{z} +g_{02} \frac{\overline{z}^{2}}{2}+\cdots , \end{aligned}$$
(3.24)

from Eqs. (3.20) and (3.23), we have

$$\begin{aligned}& u_{t}(0)=\frac{1}{2} (z+\overline{z}) \cos \biggl(\frac{n x}{l} \biggr)+W^{(1)}_{20}(0) \frac{z^{2}}{2}+W^{(1)}_{11}(0) z \overline{z}+W^{(1)}_{02}(0) \frac{\overline{z}^{2}}{2}+\cdots , \\& v_{t}(0)=\frac{1}{2} (\xi z+\overline{\xi }\overline{z}) \cos \biggl( \frac{n x}{l} \biggr)+W^{(2)}_{20}(0) \frac{z^{2}}{2}+W^{(2)}_{11}(0) z \overline{z}+W^{(2)}_{02}(0) \frac{\overline{z}^{2}}{2}+\cdots , \\& \begin{aligned} u_{t}(-1)={}&\frac{1}{2}\bigl(z e^{-i\omega _{n} \tilde{\tau }}+ \overline{z} e^{i \omega _{n} \tilde{\tau }}\bigr)\cos \biggl(\frac{n x}{l} \biggr)+W^{(1)}_{20}(-1) \frac{z^{2}}{2}+W^{(1)}_{11}(-1)z \overline{z} \\ &{} +W^{(1)}_{02}(-1)\frac{\overline{z}^{2}}{2}+\cdots , \end{aligned} \\& \begin{aligned} v_{t}(-1)={}&\frac{1}{2}\bigl(\xi z e^{-i\omega _{n} \tilde{\tau }}+ \overline{\xi }\overline{z} e^{i\omega _{n} \tilde{\tau }}\bigr)\cos \biggl( \frac{n x}{l}\biggr)+W^{(2)}_{20}(-1) \frac{z^{2}}{2}+W^{(2)}_{11}(-1)z \overline{z} \\ &{} +W^{(2)}_{02}(-1)\frac{\overline{z}^{2}}{2}+\cdots , \end{aligned} \end{aligned}$$

and

$$\begin{aligned}& \overline{F}_{1}(U_{t},0)= \frac{1}{\tilde{\tau }}F_{1}=\frac{1}{2} \alpha _{1} u_{t}^{2}(0)+\alpha _{2} u_{t}(0) v_{t}(0)+\frac{1}{2} \alpha _{3} v_{t}^{2}(0)+\cdots, \end{aligned}$$
(3.25)
$$\begin{aligned}& \begin{aligned} \overline{F}_{2}(U_{t},0)={}& \frac{1}{\tilde{\tau }}F_{2}=\frac{1}{2} \beta _{1} u_{t}^{2}(-1)+ \beta _{2} u_{t}(-1)v_{t}(-1)+\frac{1}{2} \beta _{3} v_{t}^{2}(-1) \\ &{}+\beta _{4} u_{t}(-1)v_{t}(0)+ \frac{1}{2} \beta _{5} v_{t}^{2}(0)+\cdots, \end{aligned} \end{aligned}$$
(3.26)

with

$$ \begin{aligned} &\alpha _{1}=-\frac{2 r}{k}- \frac{\alpha c(1-\beta ){{u}^{*}}^{\alpha -2}[\alpha -1-(\alpha +1)ch(1-\beta ){{u}^{*}}^{\alpha }]{{v}^{*}}}{[1+ch(1-\beta ){{u}^{*}}^{\alpha }]^{3}}, \\ &\alpha _{2}=- \frac{\alpha c(1-\beta ){{u}^{*}}^{\alpha -1}}{[1+ch(1-\beta ){{u}^{*}}^{\alpha }]^{2}}, \qquad \alpha _{3}=0, \\ &\beta _{1}=- \frac{\alpha ec (1-\beta ){{u}^{*}}^{\alpha -2} [\alpha -1-(\alpha +1)ch(1-\beta ){{u}^{*}}^{\alpha }]{{v}^{*}}}{[1+ch(1-\beta ){{u}^{*}}^{\alpha }]^{3}}, \\ &\beta _{2}= \frac{\alpha dr(1-\frac{{u}^{*}}{k})}{c(1-\beta ){{u}^{*}}^{\alpha }{{v}^{*}}},\qquad \beta _{3}=\beta _{4}=\beta _{5}=0. \end{aligned} $$

Therefore,

$$\begin{aligned}& \begin{aligned} \overline{F}_{1}(U_{t},0)={}& \cos ^{2}\biggl(\frac{n x}{l}\biggr) \biggl( \frac{z^{2}}{2} \chi _{20} +z\overline{z}\chi _{11} +\frac{\overline{z}^{2}}{2} \overline{\chi }_{20}\biggr) \\ &{}+\frac{z^{2}\overline{z}}{2} \cos \frac{nx}{l}\biggl[W_{11}^{(1)}(0) (\alpha _{1}+\xi \alpha _{2} )+W_{11}^{(2)}(0) \alpha _{2}+ W_{20}^{(1)}(0) \frac{\alpha _{1}+\overline{\xi } \alpha _{2}}{2} \\ &{}+W_{20}^{(2)}(0) \frac{\alpha _{2}}{2}\biggr]+\cdots , \end{aligned} \end{aligned}$$
(3.27)
$$\begin{aligned}& \begin{aligned} \overline{F}_{2}(U_{t},0)={}& \cos ^{2}\biggl(\frac{n x}{l}\biggr) \biggl( \frac{z^{2}}{2} \varsigma _{20} +z\overline{z}\varsigma _{11} + \frac{\overline{z}^{2}}{2}\overline{\varsigma }_{20} \biggr) \\ &{}+\frac{z^{2}\overline{z}}{2}\cos \frac{nx}{l}\biggl[W_{11}^{(1)}(-1) \bigl( e^{-i \tilde{\tau } \omega _{n}} \beta _{1}+\xi \beta _{2} \bigr) \\ &{}+W_{20}^{(1)}(-1) \biggl(\frac{1}{2} e^{i \tilde{\tau } \omega _{n}} \beta _{1}+\frac{1}{2}\overline{\xi } \beta _{2} \biggr)+ e^{-i \tilde{\tau } \omega _{n}} W_{11}^{(2)}(-1) \beta _{2} \\ &{}+\frac{1}{2} e^{i \tilde{\tau } \omega _{n}} W_{20}^{(2)}(-1) \beta _{2}\biggr] +\cdots , \end{aligned} \end{aligned}$$
(3.28)
$$\begin{aligned}& \begin{aligned} \bigl\langle F(U_{t},0),f_{n}\bigr\rangle ={}& \frac{z^{2}}{2}\tilde{\tau } \begin{pmatrix} \chi _{20} \\ \varsigma _{20} \end{pmatrix}\Gamma +z\overline{z} \tilde{\tau } \begin{pmatrix} \chi _{11} \\ \varsigma _{11} \end{pmatrix}\Gamma +\frac{\overline{z}^{2}}{2}\tilde{\tau } \begin{pmatrix} \overline{\chi }_{20} \\ \overline{\varsigma }_{20} \end{pmatrix}\Gamma \\ &{}+\frac{z^{2}\overline{z}}{2}\tilde{\tau } \begin{pmatrix} \kappa _{1} \\ \kappa _{2} \end{pmatrix} +\cdots , \end{aligned} \end{aligned}$$
(3.29)

with

$$\begin{aligned}& \Gamma =\frac{1}{l\pi } \int _{0}^{l\pi }\cos ^{3}\biggl( \frac{n x}{l}\biggr)\,dx, \\& \begin{aligned} \kappa _{1}={}&\biggl[(\alpha _{1}+\xi \alpha _{2})W^{(1)}_{11} (0) + \alpha _{2} W^{(2)}_{11} (0) +\frac{1}{2}( \alpha _{1}+\overline{\xi } \alpha _{2})W^{(1)}_{20}(0) \\ &{}+\frac{1}{2}\alpha _{1} W^{(2)}_{20} (0)\biggr]\frac{1}{l\pi } \int _{0}^{l \pi }\cos ^{2}\biggl( \frac{n x}{l}\biggr)\,dx, \end{aligned} \\& \begin{aligned} \kappa _{2}={}&\biggl[\bigl(e^{-i \tilde{\tau } \omega _{n}} \beta _{1}+\xi \beta _{2}\bigr)W^{1}_{11}(-1) + \biggl(\frac{1}{2} e^{i \tilde{\tau } \omega _{n}} \beta _{1}+ \frac{1}{2}\overline{\xi } \beta _{2} \biggr)W^{1}_{20}(-1) \\ &{}+ e^{-i \tilde{\tau } \omega _{n}} \beta _{2} W^{2}_{11}(0) + \frac{1}{2} e^{i \tilde{\tau } \omega _{n}} \beta _{2} W^{2}_{20}(0)\biggr] \frac{1}{l\pi } \int _{0}^{l\pi }\cos ^{2}\biggl( \frac{n x}{l}\biggr)\,dx, \end{aligned} \end{aligned}$$

and

$$ \begin{aligned} &\chi _{20}= \frac{1}{4} (\alpha _{1}+2 \xi \alpha _{2} ),\qquad \chi _{11}= \frac{1}{4} \bigl(\alpha _{1}+( \overline{\xi } +\xi ) \alpha _{2} \bigr), \\ &\overline{\chi }_{20}=\frac{1}{4} (\alpha _{1}+2 \overline{\xi } \alpha _{2} ), \qquad \varsigma _{20}= \frac{1}{4} e^{-2 i \tilde{\tau } \omega _{n}} (\beta _{1}+ 2\xi \beta _{2} ), \\ &\varsigma _{11}=\frac{1}{4} \beta _{1}+ \frac{1}{4} \overline{\xi } \beta _{2}+\frac{1}{4} \xi \beta _{2},\qquad \overline{\varsigma }_{20}= \frac{1}{4} e^{2 i \tilde{\tau } \omega _{n}} (\beta _{1}+ 2 \overline{\xi }\beta _{2} ). \end{aligned} $$
(3.30)

Denote

$$ \Psi _{1}(0)-i\Psi _{2}(0):=(\gamma _{1}, \gamma _{2}). $$

and because

$$ \frac{1}{l\pi } \int _{0}^{l\pi }\cos ^{3}\biggl( \frac{n x}{l}\biggr)\,dx=0, \quad n\in \mathbb{N}, $$

we can obtain

$$ \begin{aligned} &\bigl(\Psi _{1}(0)-i \Psi _{2}(0)\bigr)\bigl\langle F(U_{t},0),f_{n}\bigr\rangle \\ &\quad =\frac{z^{2}}{2} (\gamma _{1}\chi _{20} + \gamma _{2}\varsigma _{20}) \Gamma \tilde{\tau } +z\overline{z}( \gamma _{1}\chi _{11} + \gamma _{2} \varsigma _{11})\Gamma \tilde{\tau } +\frac{\overline{z}^{2}}{2}( \gamma _{1}\overline{\chi }_{20} + \gamma _{2} \overline{\varsigma }_{20}) \Gamma \tilde{\tau } \\ &\qquad {}+\frac{z^{2}\overline{z}}{2}\tilde{\tau } [\gamma _{1} \kappa _{1} + \gamma _{2} \kappa _{2} ]+\cdots , \end{aligned} $$
(3.31)

then, by (3.23), (3.24) and (3.31), we get \(g_{20}=g_{11}=g_{02}=0\), for \(n\in {{\mathbb{N}}}\). When \(n=0\), \(g_{20}=\gamma _{1}\tilde{\tau } \chi _{20} + \gamma _{2}\tilde{\tau } \varsigma _{20}\), \(g_{11}=\gamma _{1}\tilde{\tau }\chi _{11}+ \gamma _{2} \tilde{\tau }\varsigma _{11}\), \(g_{02}=\gamma _{1}\tilde{\tau } \overline{\chi }_{20} + \gamma _{2}\tilde{\tau }\overline{\varsigma }_{20} \). When \(n\in \mathbb{N}_{0}\), \(g_{21}=\tilde{\tau }( \gamma _{1} \kappa _{1} +\gamma _{2} \kappa _{2})\).

Next, we shall compute \(W_{20}(\theta )\) and \(W_{11}(\theta )\) (\(\theta \in [-1,0]\)) to get \(g_{21}\). We have

$$\begin{aligned}& \dot{W}(z,\overline{z})=W_{20}z\dot{z}+W_{11}\dot{z} \overline{z} +W_{11} z \dot{\overline{z}} +W_{02} \overline{z} \dot{\overline{z}}+\cdots , \\& A_{\tilde{\tau }} W(z,\overline{z})=A_{\tilde{\tau }} W_{20} \frac{z^{2}}{2}+A_{\tilde{\tau }} W_{11} z \overline{z} + A_{ \tilde{\tau }} W_{02} \frac{\overline{z}^{2}}{2}+\cdots , \end{aligned}$$

and \(W(z,\overline{z})\) satisfies

$$ \dot{W}(z,\overline{z})= A_{\tilde{\tau }} W+H(z,\overline{z}), $$

where

$$ \begin{aligned} H(z,\overline{z})&=H_{20} \frac{z^{2}}{2}+W_{11} z \overline{z} +H_{02} \frac{\overline{z}^{2}}{2}+\cdots \\ &=X_{0}F(U_{t},0)-\Phi \bigl(\Psi ,\bigl\langle X_{0} F(U_{t},0),f_{n}\bigr\rangle \cdot f_{n}\bigr). \end{aligned} $$
(3.32)

Hence, we can obtain

$$ (2i\omega _{n} \tilde{\tau }- A_{\tilde{\tau }})W_{20}=H_{20},\qquad -A_{ \tilde{\tau }}W_{11}=H_{11}, \qquad (-2i\omega _{n} \tilde{\tau }- A_{ \tilde{\tau }})W_{02}=H_{02}, $$
(3.33)

that is,

$$ W_{20}=(2i\omega _{n} \tilde{\tau }- A_{\tilde{\tau }})^{-1}H_{20},\qquad W_{11}=-A_{ \tilde{\tau }}^{-1} H_{11},\qquad W_{02}=(-2i\omega _{n} \tilde{\tau }- A_{ \tilde{\tau }})^{-1}H_{02}. $$
(3.34)

By (3.31), we obtain, for \(\theta \in [-1,0)\),

$$ \begin{aligned} H(z,\overline{z})&=-\Phi (0) \Psi (0)\bigl\langle F(U_{t},0),f_{n}\bigr\rangle \cdot f_{n} \\ &=-\biggl(\frac{p_{1}(\theta )+p_{2}(\theta )}{2}, \frac{p_{1}(\theta )-p_{2}(\theta )}{2i}\biggr) \begin{pmatrix} \Phi _{1}(0) \\ \Phi _{2}(0) \end{pmatrix} \bigl\langle F(U_{t},0),f_{n}\bigr\rangle \cdot f_{n} \\ &=-\frac{1}{2}\bigl[p_{1}(\theta ) \bigl(\Phi _{1}(0)-i\Phi _{2}(0)\bigr)+p_{2}( \theta ) \bigl(\Phi _{1}(0)+i\Phi _{2}(0)\bigr)\bigr] \bigl\langle F(U_{t},0),f_{n}\bigr\rangle \cdot f_{n} \\ &=-\frac{1}{2}\biggl[\bigl(p_{1}(\theta )g_{20}+p_{2}(\theta )\overline{g}_{02}\bigr) \frac{z^{2}}{2} + \bigl(p_{1}(\theta )g_{11}+p_{2}( \theta )\overline{g}_{11}\bigr)z \overline{z} \\ &\quad {}+ \bigl(p_{1}( \theta )g_{02}+p_{2}(\theta )\overline{g}_{20} \bigr) \frac{\overline{z}^{2}}{2}\biggr] \\ &\quad {}+\cdots . \end{aligned} $$

According to (3.32), for \(\theta \in [-1,0)\),

$$\begin{aligned}& H_{20}(\theta )= \textstyle\begin{cases} 0, & n\in \mathbb{N}, \\ -\frac{1}{2}(p_{1}(\theta )g_{20}+p_{2}(\theta )\overline{g}_{02}) \cdot f_{0}, & n=0, \end{cases}\displaystyle \\& H_{11}(\theta )= \textstyle\begin{cases} 0, & n\in \mathbb{N}, \\ -\frac{1}{2}(p_{1}(\theta )g_{11}+p_{2}(\theta )\overline{g}_{11}) \cdot f_{0}, & n=0, \end{cases}\displaystyle \\& H_{02}(\theta )= \textstyle\begin{cases} 0, & n\in \mathbb{N}, \\ -\frac{1}{2}(p_{1}(\theta )g_{02}+p_{2}(\theta )\overline{g}_{20}) \cdot f_{0}, & n=0, \end{cases}\displaystyle \end{aligned}$$

and

$$ H(z,\overline{z}) (0)=F(U_{t},0)-\Phi \bigl(\Psi ,\bigl\langle F(U_{t},0),f_{n}\bigr\rangle \bigr)\cdot f_{n}, $$

where

H 20 ( 0 ) = { τ ˜ ( χ 20 ς 20 ) cos 2 ( n x l ) , n N , τ ˜ ( χ 20 ς 20 ) 1 2 ( p 1 ( 0 ) g 20 + p 2 ( 0 ) g 02 ) f 0 , n = 0 , H 11 ( 0 ) = { τ ˜ ( χ 11 ς 11 ) cos 2 ( n x l ) , n N , τ ˜ ( χ 11 ς 11 ) 1 2 ( p 1 ( 0 ) g 11 + p 2 ( 0 ) g 11 ) f 0 , n = 0 .
(3.35)

By the definition of \(A_{\tilde{\tau }}\) and (3.33), for \(-1\le \theta <0\), we have

$$ \dot{W}_{20}=A_{\tilde{\tau }}W_{20}=2i\omega _{n} \tilde{\tau } W_{20}+ \frac{1}{2} \bigl(p_{1}(\theta )g_{20}+p_{2}(\theta ) \overline{g}_{02}\bigr) \cdot f_{n}. $$

That is,

$$ W_{20}(\theta )=\frac{i}{2\omega _{n} \tilde{\tau }}\biggl(g_{20} p_{1}( \theta )+\frac{\overline{g}_{02}}{3}p_{2}(\theta ) \biggr)\cdot f_{n} +E_{1} e^{2i \omega _{n} \tilde{\tau } \theta }, $$

in which

$$ E_{1}= \textstyle\begin{cases} W_{20}(0), & n\in \mathbb{N}, \\ W_{20}(0)-\frac{i}{2\omega _{n} \tilde{\tau }}(g_{20} p_{1}(\theta )+ \frac{\overline{g}_{02}}{3}p_{2}(\theta ))\cdot f_{0}, & n=0. \end{cases} $$

Also using the definition of \(A_{\tilde{\tau }}\) and (3.33), for \(-1\le \theta <0\), we have

$$ \begin{aligned} &{-}\biggl(g_{20}p_{1}(0)+ \frac{\overline{g}_{02}}{3}p_{2}(0)\biggr)\cdot f_{0} +2i \omega _{n} \tilde{\tau } E_{1}-A_{\tilde{\tau }}\biggl( \frac{i}{2\omega _{n} \tilde{\tau }} \biggl(g_{20}p_{1}(0) + \frac{\overline{g}_{02}}{3}p_{2}(0)\biggr) \cdot f_{0}\biggr) \\ &\qquad {}-A_{\tilde{\tau }} E_{1}-L_{\tilde{\tau }}\biggl( \frac{i}{2\omega _{n} \tilde{\tau }} \biggl(g_{20} p_{1}(0)+ \frac{\overline{g}_{02}}{3} p_{2}(0)\biggr) \cdot f_{n}+ E_{1} e^{2i\omega _{n} \tilde{\tau } \theta }\biggr) \\ &\quad =\tilde{\tau } \begin{pmatrix} \chi _{20} \\ \varsigma _{20} \end{pmatrix}-\frac{1}{2} \bigl(p_{1}(0)g_{20}+p_{2}(0) \overline{g}_{02}\bigr)\cdot f_{0}. \end{aligned} $$

According to

$$ A_{\tilde{\tau }} p_{1}(0)+L_{\tilde{\tau }}(p_{1} \cdot f_{0})=i \omega _{0} p_{1}(0) \cdot f_{0} $$

and

$$ A_{\tilde{\tau }} p_{2}(0)+L_{\tilde{\tau }}(p_{2} \cdot f_{0})=-i \omega _{0} p_{2}(0) \cdot f_{0}, $$

we can obtain

$$ 2i\omega _{n} E_{1}-A_{\tilde{\tau }} E_{1}-L_{\tilde{\tau }} E_{1} e^{2i \omega _{n}} = \tilde{\tau } \begin{pmatrix} \chi _{20} \\ \varsigma _{20} \end{pmatrix}\cos ^{2}\biggl( \frac{n x}{l}\biggr),\quad n\in {{\mathbb{N}}_{0}}. $$

That is,

$$ E_{1}=\tilde{\tau }E \begin{pmatrix} \chi _{20} \\ \varsigma _{20} \end{pmatrix}\cos ^{2} \biggl(\frac{n x}{l}\biggr), $$

in which

$$ E= \begin{pmatrix} 2i \omega _{n} \tilde{\tau }+ d_{1} \frac{n^{2}}{l^{2}} -a_{1} & -a_{2} \\ -c_{1} e^{-2i \omega _{n} \tilde{\tau }} & 2i \omega _{n} \tilde{\tau }+ d_{2} \frac{n^{2}}{l^{2}}+d-de^{-2i \omega _{n} \tilde{\tau }} \end{pmatrix}^{-1}. $$

Similarly, from (3.34), we can obtain

$$ -\dot{W}_{11}=\frac{i}{2\omega _{n} \tilde{\tau }}\bigl(p_{1}(\theta )g_{11}+p_{2}( \theta )\overline{g}_{11} \bigr)\cdot f_{n},\quad -1\le \theta < 0. $$

That is,

$$ W_{11}(\theta )=\frac{i}{2\omega _{n} \tilde{\tau }}\bigl(p_{1}(\theta ) \overline{g}_{11}-p_{1}(\theta )g_{11} \bigr)+E_{2}. $$

Similar to the procedure of computing \(W_{20}\), we have

$$ E_{2}=\tilde{\tau }E^{*} \begin{pmatrix} \chi _{11} \\ \varsigma _{11} \end{pmatrix} \cos ^{2}\biggl(\frac{n x}{l}\biggr), $$

where

$$ E^{*}= \begin{pmatrix} d_{1} \frac{n^{2}}{l^{2}} -a_{1} & -a_{2} \\ -c_{1} & d_{2} \frac{n^{2}}{l^{2}} \end{pmatrix}^{-1}. $$

Thus, the following quantities which determine the direction and the stability of bifurcating periodic solutions can be obtained:

$$ \textstyle\begin{cases} c_{1}(0)=\frac{i}{2\omega _{n} \tilde{\tau }}(g_{20}g_{11}-2 \vert g_{11} \vert ^{2}- \frac{ \vert g_{02} \vert ^{2}}{3}) +\frac{1}{2}g_{21}, \qquad \mu _{2}=- \frac{\operatorname{Re}(c_{1}(0))}{\operatorname{Re}(\lambda ' (\tau ^{j,+}_{n}))}, \\ T_{2}=-\frac{1}{\omega _{n} \tilde{\tau }}[\operatorname{Im}(c_{1}(0))+\mu _{2} \operatorname{Im}( \lambda '(\tau ^{j,+}_{n}))], \qquad \beta _{2}= 2 \operatorname{Re}(c_{1}(0)). \end{cases} $$
(3.36)

Theorem 3.2

For any critical value \(\tau _{n}^{j,+}\), we have

  1. (i)

    If \(\mu _{2}>0\) (\(\mu _{2}<0\)), then Hopf bifurcation is forward (backward), that is, the bifurcating periodic solutions exists for \(\tau >\tau _{n}^{j,+}\) (\(\tau <\tau _{n}^{j,+}\)).

  2. (ii)

    If \(\beta _{2}<0\) (\(\beta _{2}>0\)), then the bifurcating periodic solutions are orbitally asymptotically stable (unstable).

  3. (iii)

    If \(T_{2}>0\) (\(T_{2}<0\)), then the period increases (decreases).

4 Numerical simulation

In model (2.1), if \(\beta =0\), i.e., the habitat complexity effect is 0, it means that all prey may be caught by the predator. When \(\beta \to 1\), it means that only little of preys can be caught, i.e., most prey are protected.

1. Stability of the system without delay.

In the system (2.1), select the parameters as

$$\begin{aligned}& r=0.9,\qquad K=300,\qquad c=0.46,\qquad e=0.58, \\& h=0.053,\qquad d=0.6,\qquad \alpha =1. \end{aligned}$$

At \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\), when \(0<\beta <0.467\), the system (2.1) is unstable, when \(0.467<\beta <0.834\), the system is locally asymptotically stable. Fix \(\beta =0.5\), by calculation, \({{P}^{*}}=({{u}^{*}},{{v}^{*}})=(99.568,5.787)\) (see Figs. 14).

Figure 1
figure 1

\(\beta =0.5\), \({{d}_{1}}=1\), \({{d}_{2}}=0.3\) and the initial condition is \((99.568,5.787)\). \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\) is locally asymptotically stable

Figure 2
figure 2

\(\beta =0.9\), \({{d}_{1}}=1\), \({{d}_{2}}=0.3\) and the initial condition is \((10,40)\). \({{P}_{1}}=(300,0)\) shows global asymptotic stability

Figure 3
figure 3

\(\beta =0.3\), \({{d}_{1}}=1\), \({{d}_{2}}=0.3\) and the initial condition is \((99.568,5.787)\). The system produces periodic solutions

Figure 4
figure 4

\({{d}_{1}}=0.6888\), \({{d}_{2}}=0.16\), \(\beta =0.213\) and the initial condition is \((62.23,4.219)\). The system produces periodic solutions

2. Stability of the system with delay.

(1) For the system (2.1), let \(\alpha =1\), we choose other parameters:

$$\begin{aligned}& {{d}_{1}}=1, \qquad {{d}_{2}}=0.5, \qquad r=2.65, \qquad K=300, \qquad c=0.23, \\& e=0.115,\qquad h=0.054,\qquad d=1.16. \end{aligned}$$

By Theorem 2.1, we know that when \(0<\beta <0.0895\), the equilibrium \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\) is unstable; when \(0.0895<\beta <{{\beta }^{*}}=0.6789\), the equilibrium \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\) is locally asymptotically stable. We choose \(\beta =0.18\), by direct computation, \({{P}^{*}}=({{u}^{*}},{{v}^{*}})=(117.47,18.78)\), \(\tau _{0} ^{0} \approx 0.2236\). By the theorem, when \(\tau \in (0,\tau _{0} ^{0}]\), the equilibrium \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\) is locally asymptotically stable. When τ crosses \(\tau _{0} ^{0}\), the equilibrium \({{P} ^{*}} = ({{u} ^{*}}, {v ^{*}})\) loses stability and a Hopf bifurcation occurs (see Figs. 57).

Figure 5
figure 5

\(\tau =0.15<\tau _{0} ^{0}\) and the initial condition is \((117.47,18.78)\). \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\) is locally asymptotically stable

Figure 6
figure 6

\(\tau =0.3>\tau _{0} ^{0}\) and the initial condition is \((117.47,18.78)\). The system produces periodic solutions

Figure 7
figure 7

\(\tau = 0.15\), \(\beta = 0.06\) and the initial condition is \((117.47,18.78)\). The system produces periodic solutions

(2) For the system (2.1), let \(\alpha =2\), we choose the other parameters:

$$\begin{aligned}& {{d}_{1}}=1,\qquad {{d}_{2}}=0.5,\qquad r=3.3,\qquad K=500,\qquad c=0.045, \\& e=0.1,\qquad h=0.05,\qquad d=1.06. \end{aligned}$$

By Theorem 2.1, we know that when \(0<\beta <0.3743\), the equilibrium \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\) is unstable; when \(0.3743<\beta <{{\beta }^{*}}=0.998\), the equilibrium \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\) is locally asymptotically stable. We choose \(\beta =0.7\), by direct computation, \({{P}^{*}}=({{u}^{*}},{{v}^{*}})=(40.87,11.68)\), \(\tau _{0} ^{0} \approx 2.1563\). By the theorem, when \(\tau \in (0,\tau _{0} ^{0}]\), the equilibrium \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\) is locally asymptotically stable. When τ crosses \(\tau _{0} ^{0}\), the equilibrium \({{P} ^{*}} = ({{u} ^{*}}, {v ^{*}})\) loses stability and a Hopf bifurcation occurs (see Figs. 810).

Figure 8
figure 8

\(\tau =2<\tau _{0} ^{0}\) and the initial condition is \((40.87,11.68)\). \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\) is locally asymptotically stable

Figure 9
figure 9

\(\tau =2.2>\tau _{0} ^{0}\) and the initial condition is \((40.87,11.68)\). The system produces periodic solutions

Figure 10
figure 10

\(\tau = 2\), \(\beta = 0.16\) and the initial condition is \((40.87,11.68)\). The system produces periodic solutions

(3) In the system (2.1), let \(\alpha =3\), we choose other parameters:

$$\begin{aligned}& {{d}_{1}}=1,\qquad {{d}_{2}}=0.5,\qquad r=2.65,\qquad K=30,\qquad c=0.24, \\& e=0.07,\qquad h=0.054,\qquad d=1.16. \end{aligned}$$

By Theorem 2.1, we know that when \(0<\beta <0.6376\), the equilibrium \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\) is unstable; when \(0.6376<\beta <{{\beta }^{*}}=0.9757\), the equilibrium \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\) is locally asymptotically stable. We choose \(\beta =0.8\), by direct computation, \({{P}^{*}}=({{u}^{*}},{{v}^{*}})=(12.71,1.17)\), \(\tau _{0} ^{0} \approx 3.8624\). By the theorem, when \(\tau \in (0,\tau _{0} ^{0}]\), the equilibrium \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\) is locally asymptotically stable. When τ crosses \(\tau _{0} ^{0}\), the equilibrium \({{P} ^{*}} = ({{u} ^{*}}, {v ^{*}})\) loses stability and a Hopf bifurcation occurs (see Figs. 1113).

Figure 11
figure 11

\(\tau =3.8<\tau _{0} ^{0}\) and the initial condition is \((12.71,1.17)\). \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\) is locally asymptotically stable

Figure 12
figure 12

\(\tau =3.9>\tau _{0} ^{0}\) and the initial condition is \((12.71,1.17)\). The system produces periodic solutions

Figure 13
figure 13

\(\tau =3.9\), \(\beta =0.68\) and the initial condition is \((12.71,1.17)\). \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\) is locally asymptotically stable

5 Biological significance

From the biological standpoint, when the intensity of habitat complexity effect is higher, the magnitude of the predator population will decrease with the increase of habitat complexity effect. This is due to a lower predation rate causing predators to starve to death because of lacking sufficient food. However, when the intensity of habitat complexity effect is relatively low, the population equilibrium density of the predator will increase with the increase of habitat complexity effect, because when the habitat complexity effect is lower, the predator still has not enough food to survive continuously.

The existence of equilibrium \({{P}_{0}}= ( 0,0 )\) means the extinction of predator and prey populations. This is so because, when the intensity of habitat complexity effect is lower, the prey are quickly eaten by the predator, leading to a sharp reduction in the number of prey to extinction, and finally the predators are extinct without food.

The existence of equilibrium \({{P}_{1}}= ( K,0 )\) means the extinction of predators, which means that when the intensity of habitat complexity effect is higher, the predators cannot get food, the mortality rate of predators is higher than the growth rate, and the predators eventually die. The prey are absolutely safe and the number of prey eventually stabilizes at the maximum carrying capacity of the environment. Compared with the refuge effect, it is not difficult to find that they have the same effect on the equilibrium density of predator and prey populations. However, the difference is that the habitat complexity effect reduces the predation rate by reducing the meeting rate between predator and prey, prey are not absolutely safe. Under the refuge effect, prey are perfectly safe.

The stability of coexistence equilibrium \({{P}^{*}}=({{u}^{*}},{{v}^{*}})\) means that the system may have spatially homogeneous or inhomogeneous periodic solutions due to the existence of a diffusion term and time delay. That is to say, if the intensity of habitat complexity effect is higher, the predator’s ability is higher and there is a short digestion delay, then the predator and prey can coexist in time and space, and the population quantity will remain near the stable value. When the digest delay is close to the Hopf bifurcation value, the system may have stable periodic solutions. In this case, predators and preys can coexist, but the population quantity may show stable periodic solutions.