Over the years, predator-prey models described by the ordinary differential equations (ODEs), which were proposed and studied widely due to the pioneering theoretical works by Lotka [1] and Volterra [2].

A most crucial element in these models is the functional response, the function that describes the number of prey consumed per predator per unit time for given quantities of prey $x$ and predator $y.$

Arditi and Ginzburg [3] suggested that the essential properties of predator-dependence could be rendered by a simpler form which was called ratio-dependence. The trophic function is assumed to depend on the single variable $\frac{x}{y}$ rather than on the two separate variables $x$ and $y.$ Generally, a ratio-dependent predator-prey model of Arditi and Ginzburg [3] is

$ \left\{ \begin{array}{ccl} \dot x(t)&=&xf(x)-yp(x/y), \\ \dot y(t)&=&syq(x/y)-\delta y. \\ \end{array} \right. $

(1.1)

In this paper, we will focus our attention on the ratio-dependent type predator-prey model with Michaelis-Menten type functional response, which takes the form of

$ \left\{ \begin{array}{ccl} \dot x(t)&=&\alpha x-\frac{\beta xy}{my+x}, \\ \dot y(t)&=&y(-d_{3}+\frac{\beta_{1} x}{my+x}), \\ \end{array} \right. $

(1.2)

where $\alpha, \beta, m, d_{3}$ and $\beta_{1}$ are positive constants. $d_{3}, \beta, m$ and $\beta_{1}$ stand for the predator death rate, capturing rate, half saturation constant, and conversion rate, respectively. The dynamics of predator-prey system was studied extensively [4-11].

In order to reflect that the dynamical behaviors of models that depend on the past history of the system, it is often necessary to incorporate time-delays into the models. Suppose that in a certain environment there are the prey and predator species with respective population densities $x(t)$ and $y(t)$ at time $t.$

Based on the above discussion, by incorporating age-structure of prey and time delay into system (1.2) and suppose that the predator species captures only the adult prey species, we establish the following model

$ \left\{ \begin{array}{ccl} \dot x_{1}(t)&=&\alpha x_{2}(t)-d_{1}x_{1}(t)-\gamma x_{2}(t-\tau), \\ \dot x_{2}(t)&= &\gamma x_{2}(t-\tau)-\frac{\beta x_{2}(t)y(t)}{my(t)+x_{2}(t)}-d_{2}x_{2}(t)-\mu x_{2}^{2}(t), \\ \dot y(t)&= &\frac{\beta_{1} x_{2}(t)y(t)}{my(t)+x_{2}(t)}-d_{3}y(t), \\ \end{array} \right. $

(1.3)

where $x_{1}(t), x_{2}(t), y(t)$ represents the densities of immature prey, mature prey and predator, respectively, $\alpha, d_{1}, d_{2}, d_{3}, \beta, \gamma, m, \mu$ and $\beta_{1}$ stand for the birth rate of prey, the immature prey death rate, the mature prey death rate, the predator death rate, capturing rate, the conversion rate of immature prey to mature prey, half saturation constant, the density dependence rate of the mature prey and conversion rate, respectively. $\tau$ is called the maturation time of the prey species.

In this section, we will discuss the local stability of a positive equilibrium and the existence of Hopf bifurcations in system (1.3). It is easy to deduce that system (1.3) has a unique positive equilibrium $E^{*}=(x_{1}^{*}, x_{2}^{*}, y^{*})$ if the following holds.

(H1) $\alpha>\gamma, \beta_{1}>d_{3}$ and $(\gamma-d_{2})(\beta_{1}-d_{3})>\frac{\beta(\beta_{1}-d_{3})^{2}}{m\beta_{1}}, $ where

$ x_{1}^{*}=\frac{\alpha-\gamma}{d_{1}}x_{2}^{*}, x_{2}^{*}=\frac{md_{3}}{\beta_{1}-d_{3}}y^{*}, y^{*}=\frac{(\gamma-d_{2})(\beta_{1}-d_{3})-\frac{\beta(\beta_{1}-d_{3})^{2}}{m\beta_{1}}}{\mu md_{3}}. $

Set

$ \left. \begin{array}{ccl} f^{(1)}&=&\alpha x_{2}(t)-d_{1}x_{1}(t)-\gamma x_{2}(t-\tau),\\ f^{(2)}&=&\gamma x_{2}(t-\tau)-\frac{\beta x_{2}(t)y(t)}{my(t)+x_{2}(t)}-d_{2}x_{2}(t)-\mu x_{2}^{2}(t),\\ f^{(3)}&=&\frac{\beta_{1} x_{2}(t)y(t)}{my(t)+x_{2}(t)}-d_{3}y(t),\\ p_{1}&=&-d_{1}, p_{2}=\alpha, p_{3}=\gamma, p_{4}=d_{2}-2\gamma +\frac{\beta}{m}(1-\frac{d_{3}^{2}}{\beta_{1}^{2}}), \\ p_{5}&=&-\frac{\beta d_{3}^{2}}{\beta_{1}^{2}}, p_{6}=\frac{(\beta_{1}-d_{3})^{2}}{m\beta_{1}}, p_{7}=d_{3}(\frac{d_{3}}{\beta_{1}}-1). \end{array} \right. $

(2.1)

Then we can get the linear part of system (1.3),

$ \left\{ \begin{array}{ccl} \dot x_{1}(t)&=&p_{1}x_{1}(t)+ p_{2}x_{2}(t)-p_{3}x_{2}(t-\tau), \\ \dot x_{2}(t)&= &p_{3}x_{2}(t-\tau)+ p_{4}x_{2}(t)+p_{5}y(t), \\ \dot y(t)&=&p_{6}x_{2}(t)+p_{7}y(t), \\ \end{array} \right. $

(2.2)

where $p_{i}~(i=1, 2, 3, 4, 5, 6, 7)$ are defined in (2.1).

Therefore the corresponding characteristic equation [18] of system (2.2) is

$ (\lambda-p_{1})[\lambda^{2}-(p_{4}+p_{7})\lambda-p_{3}(\lambda-p_{7})e^{-\lambda\tau}+p_{4}p_{7}-p_{5}p_{6}]=0. $

(2.3)

It is well-known that the zero steady state of system (2.2) is asymptotically stable if all roots of eq. (2.3) have negative real parts, and is unstable if eq.(2.3) has a root with positive real part. In the following, we will study the distribution of roots of eq.(2.3).

Obviously, $\lambda=p_{1}=-d_{1} < 0$ is a negative root of eq.(2.3).

If $i\omega(\omega>0)$ is a root of eq.(2.3), then

$ \begin{aligned} \lambda^{2}-(p_{4}+p_{7})\lambda-p_{3}(\lambda-p_{7})e^{-\lambda\tau}+p_{4}p_{7}-p_{5}p_{6}=0, \\-\omega^{2}-i(p_{4}+p_{7})\omega-p_{3}(i\omega-p_{7})e^{-i\omega\tau}+p_{4}p_{7}-p_{5}p_{6}=0.\end{aligned} $

(2.3i)

Separating the real and imaginary parts of above formula gives the following equations:

$ \left. \begin{array}{ccl} -\omega^{2}=p_{5}p_{6}-p_{4}p_{7}+p_{3}\omega\sin \omega\tau -p_{3}p_{7}\cos \omega\tau, \\ -(p_{4}+p_{7})\omega =p_{3}\omega\cos \omega\tau +p_{3}p_{7}\sin \omega\tau, \\ \end{array} \right. $

(2.4)

which implies that

$ \omega^{4}+(p_{4}^{2}+p_{7}^{2}-p_{3}^{2}+2p_{5}p_{6})\omega^{2}+(p_{4}p_{7}-p_{5}p_{6})^{2}-p_{3}^{2}p_{7}^{2}=0. $

(2.5)

Let

$ \triangle=p_{4}^{4}+p_{7}^{4}+p_{3}^{4}-2p_{4}^{2}p_{7}^{2}-2p_{3}^{2}p_{4}^{2}+4p_{4}^{2}p_{5}p_{6} +2p_{3}^{2}p_{7}^{2}+4p_{5}p_{6}p_{7}^{2}-4p_{3}^{2}p_{5}p_{6}+8p_{4}p_{5}p_{6}p_{7}, $

further note that if

$ |\frac{\beta d_{3}}{m}(1+\frac{d_{3}^{2}}{\beta_{1}^{2}})(2\gamma-d_{2})d_{3}(1-\frac{d_{3}}{\beta_{1}}) |<\gamma d_{3}|1-\frac{d_{3}}{\beta_{1}} |, \hspace{10pt} \rm {(H2)} $

then $(p_{4}p_{7}-p_{5}p_{6})^{2}-p_{3}^{2}p_{7}^{2}<0, $ we will derive $\triangle>0.$

Positive root of the equation (2.4) is given by

$ \omega_{0}^{2}=\frac{(p_{3}^{2}-p_{4}^{2}-p_{7}^{2}-2p_{5}p_{6})+\sqrt{(p_{3}^{2}-p_{4}^{2}-p_{7}^{2}-2p_{5}p_{6})^{2} +4[p_{3}^{2}p_{7}^{2}-(p_{4}p_{7}-p_{5}p_{6})^{2}]}}{2}. $

(2.6)

Again from (2.4), we get

$ \tan \omega_{0}\tau_{0}=\frac{\omega_{0}(\omega_{0}^{2}+p_{5}p_{6}+p_{7}^{2} )}{p_{4}\omega_{0}^{2}+ p_{4}p_{7}^{2}-p_{5}p_{6}p_{7}}. $

Solving for $\tau_{0}, $ we get

$ \tau_{0n}=\frac{1}{\omega_{0}}\arctan[\frac{\omega_{0}(\omega_{0}^{2}+p_{5}p_{6}+p_{7}^{2} )}{p_{4}\omega_{0}^{2}+ p_{4}p_{7}^{2}-p_{5}p_{6}p_{7}}]+\frac{2n\pi}{\omega_{0}}, $

(2.7)

where $n=0, 1, 2, 3, \cdots.$

The smallest $\tau_{0}$ is obtained by choosing $n=0, $ then from (2.5), we get

$ \tau_{0}=\frac{1}{\omega_{0}}\arctan[\frac{\omega_{0}(\omega_{0}^{2}+p_{5}p_{6}+p_{7}^{2} )}{p_{4}\omega_{0}^{2}+ p_{4}p_{7}^{2}-p_{5}p_{6}p_{7}}]. $

then $(\tau_{0n}, \omega_{0})$ solves eq. (2.4). This means that when $\tau=\tau_{0n}, $ eq. (2.3). has a pair of purely imaginary roots $\pm i\omega_{0}.$

Now let us consider the behavior of roots of eq. (2.3) near $\tau_{0n}.$ Denote $\lambda(\tau)=\alpha(\tau)+i\omega(\tau)$ as the root of eq. (2.3) such that

$ \alpha(\tau_{0n})=0, \hspace{10pt} \omega(\tau_{0n})=\omega_{0}. $

Substituting $\lambda(\tau)$ into eq. (2.3i) and differentiating both sides of it with respect to $\tau, $ we have

$ (\frac{d\lambda(\tau)}{d\tau})^{-1}=\frac{(2\lambda-p_{4}-p_{7})e^{\lambda\tau} }{-p_{3}\lambda(\lambda-p_{7})}+\frac{1}{\lambda(\lambda-p_{7})}-\frac{\tau}{\lambda}. $

Noting that $\lambda=\pm i\omega_{0}, \omega_{0}$ satisfy (2.5), therefore, we have

$ \left. \begin{array}{ccl} {\rm Re}(\frac{d\lambda(\tau)}{d\tau})^{-1}_{\tau=\tau_{0n}}&=&{\rm Re}\{\frac{(2\lambda-p_{4}-p_{7})e^{\lambda\tau} }{-p_{3}\lambda(\lambda-p_{7})}\}_{\tau=\tau_{0n}} +{\rm Re}\{\frac{1}{\lambda(\lambda-p_{7})}-\frac{\tau}{\lambda}\}_{\tau=\tau_{0n}} \\ &=&{\rm Re}\{\frac{1}{p_{3}\omega_{0}^{2}+ip_{3}p_{7}\omega_{0}}[-(p_{4}+p_{7})\cos \omega_{0}\tau_{0n}-2\omega_{0}\sin \omega_{0}\tau_{0n}\\&+&i(2\omega_{0}\cos \omega_{0}\tau_{0n}- (p_{4}+p_{7})\sin \omega_{0}\tau_{0n} )]\} + {\rm Re}\{\frac{-p_{3}}{p_{3}\omega_{0}^{2}+ip_{3}p_{7}\omega_{0}}\}\\ &=&\frac{1}{p_{3}^{2}\omega_{0}^{4}+p_{3}^{2}p_{7}^{2}\omega_{0}^{2}} -p_{3}^{2}\omega_{0}^{2}-p_{3} (p_{4}+p_{7})\omega_{0}^{2} \cos \omega_{0}\tau_{0n}-2\omega_{0}^{3}p_{3} \sin \omega_{0}\tau_{0n}\\ &+& 2\omega_{0}^{2}p_{3}p_{7} \cos \omega_{0}\tau_{0n} -p_{3} p_{7} \omega_{0}(p_{4}+p_{7})\sin \omega_{0}\tau_{0n}\\ &=& \frac{1}{\Gamma}\{-(p_{4}+p_{7})\omega_{0}[p_{3}\omega_{0}\cos \omega_{0}\tau_{0n}+p_{3} p_{7}\sin \omega_{0}\tau_{0n}]\\ &-&2\omega_{0}^{2}[p_{3}\omega_{0}\sin \omega_{0}\tau_{0n} -p_{3} p_{7}\cos \omega_{0}\tau_{0n}]- p_{3}^{2}\omega_{0}^{2} \}\\ &=& \frac{1}{\Gamma}{(p_{4}+p_{7})^{2}\omega_{0}^{2}-2\omega_{0}^{2}(-\omega_{0}^{2}-p_{5}p_{6}+p_{4}p_{7})-p_{3}^{2}\omega_{0}^{2}}\\ &=&\frac{\omega_{0}^{2}}{\Gamma}{(p_{4}+p_{7})^{2}+2\omega_{0}^{2}+2(p_{5}p_{6}-p_{4}p_{7})-p_{3}^{2} }\\ &=&\frac{\omega_{0}^{2}}{\Gamma} {\sqrt{(p_{3}^{2}-p_{4}^{2}- p_{7}^{2}-2p_{5}p_{6} )^{2}+4[p_{3}^{2}p_{7}^{2}-(p_{4}p_{7}-p_{5}p_{6})^{2}]}},\\ \end{array} \right. $

where we have used eqs. (2.4), (2.6) and $\Gamma=p_{3}^{2}\omega_{0}^{4}+p_{3}^{2}p_{7}^{2}\omega_{0}^{2}>0.$ Hence, it follows that

$ \left. \begin{array}{l}{\rm sign}{{\rm Re}(\frac{d\lambda(\tau)}{d\tau})_{\tau=\tau_{0n}}} ={\rm sign}{{\rm Re}(\frac{d\lambda(\tau)}{d\tau})^{-1}_{\tau=\tau_{0n}}}\\ ={\rm sign}\{\frac{\omega_{0}^{2}}{\Gamma} {\sqrt{(p_{3}^{2}-p_{4}^{2}- p_{7}^{2}-2p_{5}p_{6} )^{2}+4[p_{3}^{2}p_{7}^{2}-(p_{4}p_{7}-p_{5}p_{6})^{2}]}}\}>0. \end{array} \right. $

(2.8)

Therefore, when the delay $\tau$ near $\tau_{0n}$ is increased, the root of eq. (2.3) crosses the imaginary axis from left to right. In addition, note that when $\tau=0, $ eq. (2.3) has roots with negative real parts only if

(H3) $p_{3}+p_{4}+p_{7}<0$ and $p_{3} p_{7}+p_{4} p_{7}-p_{5} p_{6}>0.$

Thus, summarizing the above remarks and the well-known Rouche theorem, we have the following results on the distribution of roots of eq. (2.3).

Lemma 2.1. Let $\tau_{0n} (n=0, 1, 2, \cdots)$ be defined as in (2.7). Then all roots of eq.(2.3) have negative real parts for all $\tau\in [0, \tau_{0}).$ However, eq. (2.3) has at least one root with positive real part when $\tau> \tau_{0}, $ and eq.(2.3) has a pair of purely imaginary root $\pm i\omega_{0}, $ when $\tau= \tau_{0}.$ More detail, for $\tau\in (\tau_{0n}, \tau_{0n+1}](n=0, 1, 2, \cdots), $ eq. (2.3) has $2(n+1)$ roots with positive real parts. Moreover, all roots of eq. (2.3) with $\tau=\tau_{0n}(n=0, 1, 2, \cdots) $ have negative real parts except $\pm i\omega_{0}.$

Applying Lemma 2.1, Theorem 11.1 developed in [12], we have the following results.

Theorem 2.1. Let (H1), (H2) and (H3) hold. Let $\omega_{0}$ and $\tau_{0n}(n=0, 1, 2, \cdots) $ be defined as in (2.6) and (2.7), respectively.

(ⅰ) The positive equilibrium $E^{*}$ of system (1.2) is asymptotically stable for all $\tau\in [0, \tau_{0})$ and unstable for $\tau>\tau_{0}.$

(ⅱ) System (1.3) undergoes a Hopf Bifurcation at the positive equilibrium $E^{*}$ when $\tau=\tau_{0n}(n=0, 1, 2, \cdots) $.

In section 2, we have proven that system (1.3) has a series of periodic solutions bifurcating from the positive equilibrium $E^{*}$ at the critical values of $\tau.$ In this section, we derive explicit formulae to determine the properties of the Hopf bifurcation at critical values $\tau_{0n}$ by using the normal form theory and center manifold reduction[13].

Without loss of generality, denote the critical values $\tau_{0n}$ by $\bar{\tau}, $ and set $\tau=\bar{\tau}+\mu.$ Then $\mu=0$ is a Hopf bifurcation value of system (1.3). Thus, we can work in the phase space $C=C([-\bar{\tau}, 0], R^{3}).$

Let $u_{1}(t)=x_{1}(t)-x_{1}^{*}, u_{2}(t)=x_{2}(t)-x_{2}^{*}, u_{3}(t)=y(t)-y^{*}.$ Then system (1.3) is transformed into

$ \left\{ \begin{array}{ccl} \dot u_{1}(t)&=&p_{1}u_{1}(t)+ p_{2}u_{2}(t)-p_{3}u_{2}(t-\tau)+\sum\limits_{i+j+l\geq 2}\frac{1}{i!j!l!}f^{1}_{ijl}u_{1}^{i}(t)u_{2}^{j}(t)u_{2}^{l}(t-\bar{\tau}), \\ \dot u_{2}(t)&= &p_{3}u_{2}(t-\tau)+ p_{4}u_{2}(t)+p_{5}u_{3}(t)+\sum\limits_{i+j+l\geq 2}\frac{1}{i!j!l!}f^{2}_{ijl}u_{2}^{i}(t-\bar{\tau})u_{2}^{j}(t)u_{3}^{l}(t), \\ \dot u_{3}(t)&=&p_{6}u_{2}(t)+p_{7}u_{3}(t)+\sum\limits_{i+j\geq 2}\frac{1}{i!j!}f^{3}_{ij}u_{2}^{i}(t)u_{3}^{j}(t), \\ \end{array} \right. $

(3.1)

where

$ f^{1}_{ijl}=\frac{\partial^{i+j+l}f^{(1)}}{\partial x_{1}^{i} \partial x_{2}^{j} \partial x_{2}^{l}(t-\bar{\tau})}\mid_{(x_{1}^{*}, x_{2}^{*}, y^{*})}, f^{2}_{ijl}=\frac{\partial^{i+j+l}f^{(2)}}{\partial x_{2}^{i}(t-\bar{\tau}) \partial x_{2}^{j} \partial y^{l}}\mid_{(x_{1}^{*}, x_{2}^{*}, y^{*})}, \\ f^{3}_{ij}=\frac{\partial^{i+j}f^{(3)}}{\partial x_{2}^{i} \partial y^{j}}\mid_{(x_{1}^{*}, x_{2}^{*}, y^{*})}, i, j, l\geq 0, $

here $f^{(1)}$, $f^{(2)}$ and $f^{(3)}$ are defined in (2.1).

For the simplicity of notations, we rewrite (3.1) as

$ \dot{u}(t)=L_{\mu}(u_{t})+f(\mu, u_{t}), $

(3.2)

where $u(t)=(u_{1}(t), u_{2}(t), u_{3}(t))^{T}\in R^{3}, u_{t}(\theta)\in C $ is defined by $u_{t}(\theta)=u(t+\theta), $ and $L_{\mu}: C\rightarrow R, f: R\times C\in R$ are given by

$ L_{\mu}\phi=\left( \begin{array}{ccc} p_1&p_{2}& 0 \\ 0&p_{4}&p_{5} \\ 0 &p_{6}& p_{7} \end{array} \right)\phi(0)+ \left( \begin{array}{ccc} 0&-p_{3}& 0 \\ 0&p_{3}&0 \\ 0&0& 0 \end{array} \right)\phi(-\bar{\tau}) $

(3.3)

and

$ f(\mu, \phi)=\left( \begin{array}{c} \sum\limits_{i+j+l\geq 2}\frac{1}{i!j!l!}f^{1}_{ijl}\phi_{1}^{i}(0)\phi_{2}^{j}(0)\phi_{2}^{l}(-\bar{\tau})\\ \sum\limits_{i+j+l\geq 2}\frac{1}{i!j!l!}f^{2}_{ijl}\phi_{2}^{i}(-\bar{\tau})\phi_{2}^{j}(0)\phi_{3}^{l}(0)\\ \sum\limits_{i+j\geq 2}\frac{1}{i!j!}f^{3}_{ij}\phi_{2}^{i}(0)\phi_{3}^{j}(0) \end{array} \right) , $

(3.4)

respectively. By Riesz representation theorem, there exists a function $\eta(\theta, \mu)$ of bounded variation for $\theta\in [-\bar{\tau}, 0]$ such that

$ L_{\mu}\phi=\int_{-\bar{\tau}}^{0}d\eta(\theta, 0)\phi(\theta) $

(3.5)

for $\phi \in C.$

In fact, we can choose

$ \eta(\theta, \mu)=\left( \begin{array}{ccc} p_1&p_{2}& 0 \\ 0&p_{4}&p_{5} \\ 0 &p_{6}& p_{7} \end{array} \right)\delta(\theta)- \left( \begin{array}{ccc} 0&-p_{3}& 0 \\ 0&p_{3}&0 \\ 0&0& 0 \end{array} \right)\delta(\theta+\bar{\tau}), $

(3.6)

where $\delta$ is the vector, whose components are the Dirac delta functions. For $\phi\in C^{1}([-\bar{\tau}, 0], R^{3}), $ define

$ A(\mu)\phi=\left\{ \begin{array}{cc} \frac{d\phi(\theta)}{d\theta},&\theta\in [-\bar{\tau}, 0), \\ \int_{-\bar{\tau}}^{0}d\eta(\mu, s)\phi(s),&\theta=0 \\ \end{array} \right. $

and

$ R(\mu)\phi=\left\{ \begin{array}{cc} 0,&\theta\in [-\bar{\tau}, 0), \\ f(\mu, \phi),&\theta=0. \\ \end{array} \right. $

Then system (3.2) is equivalent to

$ \dot{u}_{t}=A(\mu) u_{t}+R(\mu)u_{t}, $

(3.7)

where $x_{t}(\theta)=x(t+\theta)$ for $ \theta\in [-\bar{\tau}, 0].$

For $\psi\in C^{1}([0, \bar{\tau}], (R^{3})^{*}), $ define

$ A^{*}\psi(s)=\left\{ \begin{array}{cc} -\frac{d\psi(s)}{ds},&s\in (0, \bar{\tau}], \\ \int_{-\bar{\tau}}^{0}d\eta^{T}(t, 0)\psi(-t),&s=0 \\ \end{array} \right. $

and a bilinear inner product

$ 〈\psi(s), \phi(\theta)〉=\bar{\psi}(0)\phi(0)-\int_{-\bar{\tau}}^{0}\int_{\varphi=0}^{\theta}\bar{\psi}(\varphi-\theta)d\eta(\theta)\phi(\xi)d\xi, $

(3.8)

where $\eta(\theta)=\eta(\theta, 0).$ Then $A(0)$ and $A^{*}$ are adjoint operators. By discussions in Section 2 and foregoing assumption, we know that $\pm i\omega_{0}$ are eigenvalues of $A(0).$ Thus, they are also eigenvalues of $A^{*}.$ We first need to compute the eigenvector of $A(0)$ and $A^{*}$ corresponding to $i\omega_{0}$ and $-i\omega_{0}, $ respectively.

Suppose that $q(\theta)=(\rho_{1}, 1, \rho_{2})^{T}e^{i\omega_{0}\theta}$ is the eigenvector of $A(0), $ (3.5) and (3.6) that

$ \left( \begin{array}{ccc} p_1-i\omega_{0}&p_{2}-p_{3}e^{-i\omega_{0}\bar{\tau}} &0 \\ 0&p_{4}+p_{3}e^{-i\omega_{0}\bar{\tau}}-i\omega_{0}&p_{5} \\ 0 &p_{6}& p_{7}-i\omega_{0} \end{array} \right)q(0)= \left( \begin{array}{c} 0 \\ 0 \\ 0 \end{array} \right). $

We therefore derive that

$ q(0)=(\rho_{1}, 1, \rho_{2})^{T}=( \frac{p_{2}+ p_{4}-i\omega_{0}+\frac{p_{5}p_{6}}{i\omega_{0}-p_{7}}}{i\omega_{0}-p_{1}}, 1, \frac{p_{6}}{i\omega_{0}-p_{7}})^{T}. $

On the other hand, suppose that $q^{*}(s)=D(\sigma_{1}, 1, \sigma_{2})e^{i\omega_{0}s}$ is the eigenvector of $A^{*}$ corresponding to $-i\omega_{0}.$ From the definition of $A^{*}, $ (3.5) and (3.6), we have

$ \left( \begin{array}{ccc} p_1+i\omega_{0}&0 &0 \\ p_{2}-p_{3}e^{i\omega_{0}\bar{\tau}}&p_{4}+p_{3}e^{i\omega_{0}\bar{\tau}}+i\omega_{0}&p_{6} \\ 0 &p_{5}& p_{7}+i\omega_{0} \end{array} \right)(q^{*}(0))^{T}= \left( \begin{array}{c} 0 \\ 0 \\ 0 \end{array} \right), $

which yields

$ q^{*}(0)=D(\sigma_{1}, 1, \sigma_{2})=D( 0, 1, -\frac{p_{5}}{i\omega_{0}+p_{7}})^{T}. $

In order to assure $〈q^{*}(s), q(\theta)〉=1, $ we need to determine the value of D. From (3.7), we have

$ \left. \begin{array}{ccl} &-&\omega^{2}=p_{5}p_{6}-p_{4}p_{7}+p_{3}\omega\sin \omega\tau -p_{3}p_{7}\cos \omega\tau, \\ &-&(p_{4}+p_{7})\omega =p_{3}\omega\cos \omega\tau +p_{3}p_{7}\sin \omega\tau, \\ &&〈 q^{*}(s), q(\theta)〉\\ &=&\bar{D}{(0, 1, \bar{\sigma}_{2})(\rho_{1}, 1, \rho_{2})^{T}-\int_{-\bar{\tau}}^{0}\int_{\xi=0}^{\theta}(0, 1, \bar{\sigma}_{2})e^{-i(\xi-\theta)\omega_{0}} d\eta(\theta)(\rho_{1}, 1, \rho_{2})^{T}e^{i\xi\omega_{0}}d\xi}\\ &=&\bar{D}{1+\rho_{2} \bar{\sigma}_{2}- \int_{-\bar{\tau}}^{0}(0, 1, \bar{\sigma}_{2})\theta e^{i\theta\omega_{0}}d\eta(\theta)(\rho_{1}, 1, \rho_{2})^{T}}\\ &=&\bar{D}{1+\rho_{2} \bar{\sigma}_{2}+p_{3}\bar{\tau}e^{-i\bar{\tau}\omega_{0}}}.\\ \end{array} \right. $

Thus, we can choose

$ D=\frac{1}{1+\rho_{2} \bar{\sigma}_{2}+p_{3}\bar{\tau}e^{i\bar{\tau}\omega_{0}}}, $

such that $〈q^{*}(s), q(\theta)〉 =1, $ $〈q^{*}(s), \bar{q}(\theta)〉 =0.$

In the remainder of this section, we will compute the coordinates to describe the center manifold $C_{0}$ at $\mu=0.$ Let $u_{t}$ be the solution of eq.(3.2) with $\mu=0.$

Define

$ z(t)=〈q^{*}, u_{t}〉, W(t, \theta)=u_{t}(\theta)-2{\rm Re}\{z(t)q(\theta)\}. $

(3.9)

On the center manifold $C_{0}$ we have $W(t, \theta)=W(z(t), \bar{z}(t), \theta), $ where

$ W(z(t), \bar{z}(t), \theta)=W_{20}(\theta)\frac{z^{2}}{2}+W_{11}(\theta)z\bar{z}+W_{02}(\theta)\frac{\bar{z}^{2}}{2}+\cdots, $

$z$ and $\bar{z}$ are local coordinates for center manifold $C_{0}$ in the direction of $q^{*}$ and $\bar{q}^{*}.$

Note that $W$ is real if $u_{t}$ is real. We only consider real solutions. For the solution $u_{t}\in C_{0}$ of (3.2), since $\mu=0, $ we have

$ \left. \begin{array}{l}\dot{z}=i\omega_{0}z+〈\bar{q}^{*}(\theta), f(0, W(z(t), \bar{z}(t), \theta)+2{\rm Re}\{zq(\theta)\})〉\\ =i\omega_{0}z+\bar{q}^{*}(\theta)f(0, W(z(t), \bar{z}(t), \theta)+2{\rm Re}\{zq(\theta)\})\\ =i\omega_{0}z+\bar{q}^{*}(0)f(0, W(z(t), \bar{z}(t), 0)+2{\rm Re}\{zq(0)\})\\ : =i\omega_{0}z+\bar{q}^{*}(0)f_{0}(z, \bar{z}). \end{array} \right. $

(3.10)

We rewrite (3.10) as $\dot{z}=i\omega_{0}z+g(z, \bar{z})$ with

$ g(z,\bar{z})=\bar{q}^{*}(0)f_{0}(z,\bar{z})=g_{20}\frac{z^{2}}{2}+g_{11}z\bar{z}+g_{02}\frac{\bar{z}^{2}}{2}+\cdots. $

(3.11)

Noting that

$ u_{t}(\theta)=(u_{1t}(\theta), u_{2t}(\theta), u_{3t}(\theta))=W(t, \theta)+zq(\theta)+\bar{z}\bar{q}(\theta) $

and $q(\theta)=(\rho_{1}, 1, \rho_{2})^{T}e^{i\theta\omega_{0}}, $ we have

$ u_{1t}(0)=\rho_{1}z+\bar{\rho}_{1}\bar{z}+W^{(1)}_{20}(0)\frac{z^{2}}{2}+W^{(1)}_{11}(0)z\bar{z}+W^{(1)}_{02}(0)\frac{\bar{z}^{2}}{2}+\cdots, \\u_{2t}(0)=z+\bar{z}+W^{(2)}_{20}(0)\frac{z^{2}}{2}+W^{(2)}_{11}(0)z\bar{z}+W^{(2)}_{02}(0)\frac{\bar{z}^{2}}{2}+\cdots, \\u_{2t}(-\bar{\tau})=e^{-i\omega_{0}\bar{\tau}}z+e^{i\omega_{0}\bar{\tau}}\bar{z}+W^{(2)}_{20}(-\bar{\tau})\frac{z^{2}}{2}+W^{(2)}_{11}(-\bar{\tau})z\bar{z} +W^{(2)}_{02}(-\bar{\tau})\frac{\bar{z}^{2}}{2}+\cdots, \\u_{3t}(0)=\rho_{2}z+\bar{\rho}_{2}\bar{z}+W^{(3)}_{20}(0)\frac{z^{2}}{2}+W^{(3)}_{11}(0)z\bar{z}+W^{(3)}_{02}(0)\frac{\bar{z}^{2}}{2}+\cdots, $

thus, it follows from (3.4) and (3.11) that

$ \left. \begin{array}{l} g(z, \bar{z})=\bar{q}^{*}(0)f_{0}(z, \bar{z})=\bar{D}(\bar{\sigma}_{1}, 1, \bar{\sigma}_{2})\left( \begin{array}{c} \sum\limits_{i+j+l\geq 2}\frac{1}{i!j!l!}f^{(1)}_{ijl}u_{1t}^{(i)}(0)u_{2t}^{(j)}(0)u_{2t}^{(l)}(-\bar{\tau})\\ \sum\limits_{i+j+l\geq 2}\frac{1}{i!j!l!}f^{(2)}_{ijl}u_{2t}^{(i)}(0)u_{2t}^{(j)}(-\bar{\tau})u_{3t}^{(l)}(0)\\ \sum\limits_{i+j\geq 2}\frac{1}{i!j!}f^{(3)}_{ij}u_{2t}^{(i)}(0)u_{3t}^{(j)}(0) \end{array} \right)\\ =\bar{D}(\bar{\sigma}_{1}(f^{(1)}_{110}\rho_{1}+f^{(1)}_{011}e^{-i\omega_{0}\bar{\tau}}+\frac{f^{(1)}_{020}}{2}+\frac{f^{(1)}_{200}}{2}\rho_{1}^{2} )z^{2}+(f^{(2)}_{110}e^{-i\omega_{0}\bar{\tau}}+f^{(2)}_{101}\rho_{2} +\frac{f^{(2)}_{200}}{2}+\frac{f^{(2)}_{002}}{2}\rho_{2}^{2})z^{2}\\ +\bar{\sigma}_{2}(f^{(3)}_{11}\rho_{2}+\frac{f^{(3)}_{20}}{2}+\frac{f^{(3)}_{02}}{2}\rho_{2}^{2} )z^{2} +\bar{\sigma}_{1}(2f^{(1)}_{110}Re{\rho_{1}}+2f^{(1)}_{011}Re{e^{-i\omega_{0}\bar{\tau}}}+f^{(1)}_{020}+f^{(1)}_{200}\rho_{1}\bar{\rho}_{1} )z\bar{z}\\ +(2f^{(2)}_{110}Re{e^{-i\omega_{0}\bar{\tau}}}+2f^{(2)}_{101}Re{\rho_{2}} +f^{(2)}_{200}+f^{(2)}_{002}\rho_{2}\bar{\rho}_{2})z\bar{z} +\bar{\sigma}_{2}(2f^{(3)}_{11}Re{\rho_{2}}+f^{(3)}_{20}+f^{(3)}_{02}\rho_{2}\bar{\rho}_{2} )z\bar{z}\\ +\bar{\sigma}_{1}(f^{(1)}_{110}\bar{\rho}_{1}+f^{1}_{011}e^{i\omega_{0}\bar{\tau}}+\frac{f^{(1)}_{020}}{2}+\frac{f^{(1)}_{200}}{2}\bar{\rho}_{1}^{2} )\bar{z}^{2}+(f^{(2)}_{110}e^{i\omega_{0}\bar{\tau}}+f^{(2)}_{101}\bar{\rho}_{2} +\frac{f^{(2)}_{200}}{2}+\frac{f^{(2)}_{002}}{2}\bar{\rho}_{2}^{2})\bar{z}^{2}\\ +\bar{\sigma}_{2}(f^{(3)}_{11}\bar{\rho}_{2}+\frac{f^{(3)}_{20}}{2}+\frac{f^{(3)}_{02}}{2}\bar{\rho}_{2}^{2} )\bar{z}^{2} +\bar{\sigma}_{1}[f^{(1)}_{110}(\frac{W_{20}^{(1)}(0)}{2}+\frac{W_{20}^{(2)}(0)}{2}\bar{\rho}_{1}+W_{11}^{(2)}(0)\rho_{1}+W_{11}^{(1)}(0))\\ +f^{(1)}_{011}(W_{11}^{(2)}(-\bar{\tau})+ \frac{W_{20}^{(2)}(-\bar{\tau})}{2}+\frac{W_{20}^{(2)}(0)}{2}e^{i\omega_{0}\bar{\tau}}+W_{11}^{(2)}(0)e^{i\omega_{0}\bar{\tau}}) +f^{(1)}_{020}(W_{11}^{(2)}(0)+ \frac{W_{20}^{(2)}(0)}{2})\\ +f^{(1)}_{200}(\rho_{1}W_{11}^{(1)}(0)+ \frac{\bar{\rho}_{1}W_{20}^{(1)}(0)}{2})]z^{2}\bar{z} +[f^{(2)}_{110}(W_{11}^{(2)}(-\bar{\tau})+ \frac{W_{20}^{(2)}(-\bar{\tau})}{2}+\frac{W_{20}^{(2)}(0)}{2}e^{i\omega_{0}\bar{\tau}}\\+W_{11}^{(2)}(0)e^{i\omega_{0}\bar{\tau}}) +f^{(2)}_{101}(W_{11}^{(3)}(0)+ \frac{W_{20}^{(3)}(0)}{2}+\frac{W_{20}^{(2)}(0)}{2}\bar{\rho}_{2}+W_{11}^{(2)}(0)\rho_{2}) +f^{(2)}_{200}(W_{11}^{(2)}(0)\\+ \frac{W_{20}^{(2)}(0)}{2}) +f^{(2)}_{002}(\rho_{2}W_{11}^{(3)}(0) + \frac{\bar{\rho}_{2}W_{20}^{(3)}(0)}{2}) ]z^{2}\bar{z} +\bar{\sigma}_{2}[ f^{(3)}_{11}(W_{11}^{(3)}(0)+\frac{W_{20}^{(3)}(0)}{2}+\frac{W_{20}^{(2)}(0)}{2}\bar{\rho}_{2}\\ +W_{11}^{(2)}(0)\rho_{2}) +f^{(3)}_{20}(W_{11}^{(2)}(0) + \frac{W_{20}^{(2)}(0)}{2}) +f^{(3)}_{02}(\rho_{2}W_{11}^{(3)}(0)+ \frac{\bar{\rho}_{2}W_{20}^{(3)}(0)}{2}) ]z^{2}\bar{z}+\cdots ).\\ \end{array} \right. $

Comparing the coefficients in (3.11), we get

$ \left. \begin{array}{l} g_{20}=(\bar{\sigma}_{1}(2f^{(1)}_{110}\rho_{1}+2f^{(1)}_{011}e^{-i\omega_{0}\bar{\tau}}+f^{(1)}_{020}+f^{(1)}_{200}\rho_{1}^{2})+( 2f^{(2)}_{110}e^{-i\omega_{0}\bar{\tau}}+2f^{(2)}_{101}\rho_{2}+f^{(2)}_{200}+f^{(2)}_{002}\rho_{2}^{2})\\ +\bar{\sigma}_{2}(2f^{(3)}_{11}\rho_{2}+f^{(3)}_{20}+f^{(3)}_{02}\rho_{2}^{2} ) )\bar{D}, \\ g_{11}=(\bar{\sigma}_{1}(2f^{(1)}_{110}Re{\rho_{1}}+2f^{(1)}_{011}Re{e^{i\omega_{0}\bar{\tau}}}+f^{(1)}_{020}+f^{(1)}_{200}\rho_{1}\bar{\rho}_{1})\\ +( 2f^{(2)}_{110}Re{e^{i\omega_{0}\bar{\tau}}}+2f^{(2)}_{101}Re{\rho_{2}}+f^{(2)}_{200}+f^{(2)}_{002}\rho_{2}\bar{\rho}_{2}) +\bar{\sigma}_{2}(2f^{(3)}_{11}Re{\rho_{2}}+f^{(3)}_{20}+f^{(3)}_{02}\rho_{2}\bar{\rho}_{2} ) )\bar{D}, \\ g_{02}=(\bar{\sigma}_{1}(2f^{(1)}_{110}\bar{\rho}_{1}+2f^{(1)}_{011}e^{i\omega_{0}\bar{\tau}}+f^{(1)}_{020}+f^{(1)}_{200}\bar{\rho}_{1}^{2})+( 2f^{(2)}_{110}e^{i\omega_{0}\bar{\tau}}+2f^{(2)}_{101}\bar{\rho}_{2}+f^{(2)}_{200}+f^{(2)}_{002}\bar{\rho}_{2}^{2})\\ +\bar{\sigma}_{2}(2f^{(3)}_{11}\bar{\rho}_{2}+f^{(3)}_{20}+f^{(3)}_{02}\bar{\rho}_{2}^{2} ) )\bar{D}, \\ g_{21}=(\bar{\sigma}_{1}[f^{(1)}_{110}(W_{20}^{(1)}(0)+W_{20}^{(2)}(0)\bar{\rho}_{1}+2W_{11}^{(2)}(0)\rho_{1}+2W_{11}^{(1)}(0)) +\\f^{(1)}_{011}(2W_{11}^{(2)}(-\bar{\tau})+ W_{20}^{(2)}(-\bar{\tau})\\+W_{20}^{(2)}(0)e^{i\omega_{0}\bar{\tau}} +2W_{11}^{(2)}(0)e^{i\omega_{0}\bar{\tau}}) +f^{(1)}_{020}(2W_{11}^{(2)}(0)+ W_{20}^{(2)}(0)) +\\f^{(1)}_{200}(2\rho_{1}W_{11}^{(1)}(0)+\bar{\rho}_{1}W_{20}^{(1)}(0))]\\ +[f^{(2)}_{110}(2W_{11}^{(2)}(-\bar{\tau})+ W_{20}^{(2)}(-\bar{\tau})+W_{20}^{(2)}(0)e^{i\omega_{0}\bar{\tau}}+2W_{11}^{(2)}(0)e^{i\omega_{0}\bar{\tau}})\\ +f^{(2)}_{101}(2W_{11}^{(3)}(0)+ W_{20}^{(3)}(0)+W_{20}^{(2)}(0)\bar{\rho}_{2}+2W_{11}^{(2) }(0)\rho_{2}) +f^{(2)}_{200}(2W_{11}^{(2)}(0)\\+ W_{20}^{(2)}(0)) +f^{(2)}_{002}(2\rho_{2}W_{11}^{(3)}(0)+ \bar{\rho}_{2}W_{20}^{(3)}(0))] +\bar{\sigma}_{2}[ f^{(3)}_{11}(2W_{11}^{(3)}(0)+W_{20}^{(3)}(0)\\ +W_{20}^{(2)}(0)\bar{\rho}_{2}+2W_{11}^{(2)}(0)\rho_{2}) +f^{(3)}_{20}(2W_{11}^{(2)}(0)+ W_{20}^{(2)}(0)) \\+f^{(3)}_{02}(2\rho_{2}W_{11}^{(3)}(0)+\bar{\rho}_{2}W_{20}^{(3)}(0))])\bar{D}, \end{array} \right. $

(3.12)

We now compute $W_{20}(\theta)$ and $W_{11}(\theta).$ It follows from (3.7) and (3.9) that

$ \dot{W}=\dot{u}_{t}-\dot{z}q-\dot{\bar{z}}\bar{q}=\left\{ \begin{array}{l} AW-2{\rm Re}\{\bar{q}^{*}(0)f(0)q(\theta)\}, \theta\in (0, \bar{\tau}], \\ AW-2{\rm Re}\{\bar{q}^{*}(0)f(0)q(0)\}+f_{0}, \theta=0 \\ \end{array} \right.=AW+H(z, \bar{z}, \theta), $

(3.13)

where

$ H(z, \bar{z}, \theta)=H_{20}(\theta)\frac{z^{2}}{2}+H_{11}(\theta)z\bar{z}+H_{02}(\theta)\frac{\bar{z}^{2}}{2}+\cdots. $

(3.14)

On the other hand, on $C_{0}$ near the origin

$ \dot{W}=W_{z}\dot{z}+ W_{\bar{z}}\dot{\bar{z}}, $

(3.15)

we derive from (3.13)-(3.15) that

$ (A-2i\omega_{0})W_{20}(\theta)=-H_{20}(\theta), AW_{11}(\theta)=-H_{11}(\theta), \cdots. $

(3.16)

It follows from (3.11) and (3.13) that for $\theta\in [-\bar{\tau}, 0), $

$ H(z, \bar{z}, \theta)=-\bar{q}^{*}(0)f(0)q(\theta)-q^{*}(0)\bar{f}(0)\bar{q}(\theta)=-gq(\theta)-\bar{g}\bar{q}(\theta). $

(3.17)

Comparing the coefficients in (3.14) and formula of $g(z, \bar{z})$ gives that for $\theta\in [-\bar{\tau}, 0), $

$ H_{20}(\theta)=-g_{20}q(\theta)-\bar{g}_{02}\bar{q}(\theta) $

(3.18)

and

$ H_{11}(\theta)=-g_{11}q(\theta)-\bar{g}_{11}\bar{q}(\theta). $

(3.19)

We derive from (3.16) and (3.18) and the definition of A that

$ \dot{W}_{20}(\theta)=2i\omega_{0}W_{20}(\theta)+g_{20}q(\theta)+\bar{g}_{02}\bar{q}(\theta). $

Note that $q(\theta)=q(0)e^{i\omega_{0}\theta}, $ it follows that

$ W_{20}(\theta)=\frac{ig_{20}}{\omega_{0}}q(0)e^{i\omega_{0}\theta} +\frac{i\bar{g}_{02}}{3\omega_{0}}\bar{q}(0)e^{-i\omega_{0}\theta}+C_{1}e^{2i\omega_{0}\theta}, $

(3.20)

where $C_{1}=(C_{1}^{(1)}, C_{1}^{(2)}, C_{1}^{(3)})\in R^{3}$ is a constant vector.

Similarly, from (3.16) and (3.19), we can obtain

$ W_{11}(\theta)=-\frac{ig_{11}}{\omega_{0}}q(0)e^{i\omega_{0}\theta} +\frac{i\bar{g}_{11}}{\omega_{0}}\bar{q}(0)e^{-i\omega_{0}\theta}+C_{2}, $

(3.21)

where $C_{2}=(C_{2}^{(1)}, C_{2}^{(2)}, C_{2}^{(3)})\in R^{3}$ is a constant vector.

In the following, we seek appropriate $C_{1}$ and $C_{2}.$ From the definition of A and (3.16), we obtain

$ \int_{-\bar{\tau}}^{0}d\eta(\theta)W_{20}(\theta)=2i\omega_{0}W_{20}(0)-H_{20}(0) $

(3.22)

and

$ \int_{-\bar{\tau}}^{0}d\eta(\theta)W_{11}(\theta)=-H_{11}(0), $

(3.23)

where $\eta(\theta)=\eta(0, \theta).$ From (3.13), it follows that

$ H_{20}(0)=-g_{20}q(0)-\bar{g}_{02}\bar{q}(0)+\left( \begin{array}{c} 2f^{(1)}_{110}\rho_{1}+2f^{(1)}_{011}e^{-i\omega_{0}\bar{\tau}}+f^{(1)}_{020}+f^{(1)}_{200}\rho_{1}^{2}\\ 2f^{(2)}_{110}e^{-i\omega_{0}\bar{\tau}}+2f^{(2)}_{101}\rho_{2}+f^{(2)}_{200}+f^{(2)}_{002}\rho_{2}^{2}\\ 2f^{(3)}_{11}\rho_{2}+f^{(3)}_{20}+f^{(3)}_{02}\rho_{2}^{2} \end{array} \right), $

(3.24)

$ H_{11}(0)=-g_{11}q(0)-\bar{g}_{11}\bar{q}(0)+\left( \begin{array}{c} 2f^{(1)}_{110}Re{\rho_{1}}+2f^{(1)}_{011}Re{e^{-i\omega_{0}\bar{\tau}}}+f^{(1)}_{020}+f^{(1)}_{200}\rho_{1}\bar{\rho}_{1}\\ 2f^{(2)}_{110}Re{e^{-i\omega_{0}\bar{\tau}}}+2f^{(2)}_{101}Re{\rho_{2}}+f^{(2)}_{200}+f^{(2)}_{002}\rho_{2}\bar{\rho}_{2}\\ 2f^{(3)}_{11}Re{\rho_{2}}+f^{(3)}_{20}+f^{(3)}_{02}\rho_{2}\bar{\rho}_{2} \end{array} \right). $

(3.25)

Substituting (3.20) and (3.24) into (3.22) and noticing that

$ (i\omega_{0}I-\int_{-\bar{\tau}}^{0}e^{i\omega_{0}\theta} d\eta(\theta))q(0)=0 $

and

$ (-i\omega_{0}I-\int_{-\bar{\tau}}^{0}e^{-i\omega_{0}\theta} d\eta(\theta))\bar{q}(0)=0, $

we obtain

$ (2i\omega_{0}I-\int_{-\bar{\tau}}^{0}e^{2i\omega_{0}\theta} d\eta(\theta))C_{1}=\left( \begin{array}{c} 2f^{(1)}_{110}\rho_{1}+2f^{(1)}_{011}e^{-i\omega_{0}\bar{\tau}}+f^{(1)}_{020}+f^{(1)}_{200}\rho_{1}^{2}\\ 2f^{(2)}_{110}e^{-i\omega_{0}\bar{\tau}}+2f^{(2)}_{101}\rho_{2}+f^{(2)}_{200}+f^{(2)}_{002}\rho_{2}^{2}\\ 2f^{(3)}_{11}\rho_{2}+f^{(3)}_{20}+f^{(3)}_{02}\rho_{2}^{2} \end{array} \right), $

which leads to

$ \left( \begin{array}{ccc} 2i\omega_{0}-p_1 &-p_{2}+p_{3}e^{-2i\omega_{0}\bar{\tau}}&0 \\ 0&2i\omega_{0}-p_{4}-p_{3}e^{-2i\omega_{0}\bar{\tau}} &-p_{5} \\ 0 &-p_{6}&2i\omega_{0}- p_{7} \end{array} \right) C_{1}\\=\left( \begin{array}{c} 2f^{(1)}_{110}\rho_{1}+2f^{(1)}_{011}e^{-i\omega_{0}\bar{\tau}}+f^{(1)}_{020}+f^{(1)}_{200}\rho_{1}^{2}\\ 2f^{(2)}_{110}e^{-i\omega_{0}\bar{\tau}}+2f^{(2)}_{101}\rho_{2}+f^{(2)}_{200}+f^{(2)}_{002}\rho_{2}^{2}\\ 2f^{(3)}_{11}\rho_{2}+f^{(3)}_{20}+f^{(3)}_{02}\rho_{2}^{2} \end{array} \right). $

It follows that

$ C_{1}=\left( \begin{array}{ccc} 2i\omega_{0}-p_1 &-p_{2}+p_{3}e^{-2i\omega_{0}\bar{\tau}}&0 \\ 0&2i\omega_{0}-p_{4}-p_{3}e^{-2i\omega_{0}\bar{\tau}}&-p_{5} \\ 0 &-p_{6}&2i\omega_{0}- p_{7} \end{array} \right)^{-1} \\\left( \begin{array}{c} 2f^{(1)}_{110}\rho_{1}+2f^{(1)}_{011}e^{-i\omega_{0}\bar{\tau}}+f^{(1)}_{020}+f^{(1)}_{200}\rho_{1}^{2}\\ 2f^{(2)}_{110}e^{-i\omega_{0}\bar{\tau}}+2f^{(2)}_{101}\rho_{2}+f^{(2)}_{200}+f^{(2)}_{002}\rho_{2}^{2}\\ 2f^{(3)}_{11}\rho_{2}+f^{(3)}_{20}+f^{(3)}_{02}\rho_{2}^{2} \end{array} \right). $

Similarly, substituting (3.21) and (3.25) into (3.23), we get

$ \left( \begin{array}{ccc} -p_1&-p_{2}+p_{3}&0 \\ 0&-p_{4}-p_{3}& -p_{5} \\ 0 &-p_{6}&- p_{7} \end{array} \right)C_{2}=\left( \begin{array}{c} 2f^{(1)}_{110}Re{\rho_{1}}+2f^{(1)}_{011}Re{e^{-i\omega_{0}\bar{\tau}}}+f^{(1)}_{020}+f^{(1)}_{200}\rho_{1}\bar{\rho}_{1}\\ 2f^{(2)}_{110}Re{e^{-i\omega_{0}\bar{\tau}}}+2f^{(2)}_{101}Re{\rho_{2}}+f^{(2)}_{200}+f^{(2)}_{002}\rho_{2}\bar{\rho}_{2}\\ 2f^{(3)}_{11}Re{\rho_{2}}+f^{(3)}_{20}+f^{(3)}_{02}\rho_{2}\bar{\rho}_{2} \end{array} \right). $

Thus, we can determine $W_{20}(\theta)$ and $W_{11}(\theta)$ from $(3.20)$ and $(3.21).$ Furthermore, we can determine $g_{21}.$ Therefore, each $g_{ij}$ in (3.12) is determined by the parameters and delay in system (3.1). Thus, we can compute the following values

$ \begin{array}{*{20}{l}} {{c_1}(0) = \frac{i}{{2{\omega _0}}}({g_{11}}{g_{20}} - 2|{g_{11}}{|^2} - \frac{{|{g_{02}}{|^2}}}{3}) + \frac{{{g_{21}}}}{2},}\\ {{\mu _2} = - \frac{{{\rm{Re}}\{ {c_1}(0)\} }}{{{\rm{Re}}\{ \lambda '(\bar \tau )\} }},{\beta _2} = 2{\rm{Re}}\{ {c_1}(0)\} ,}\\ {{T_2} = - \frac{{{\rm{Im}}{c_1}(0) + {\mu _2}Im\lambda '(\bar \tau )}}{{{\omega _0}}},} \end{array} $

(3.26)

which determine the quantities of bifurcating periodic solutions in the center manifold at the critical value $\bar{\tau}, $ i.e., $\mu_{2}$ determines the direction of the Hopf bifurcation: if $\mu_{2}>0(\mu_{2}<0), $ then the Hopf bifurcation is supercritical(subcritical) and the bifurcating periodic solutions exist for $\tau>\bar{\tau}(\tau<\bar{\tau});$ $\beta_{2}$ determines the stability of the bifurcating periodic solutions: the period increase (decrease) if $T_{2}>0(T_{2}<0).$

From what has been discussed above, we could determine the stability and direction of periodic solutions bifurcating from the positive equilibrium $E^{*}$ at the critical pint $\tau_{0n}.$