**4. Normal Form of Hopf Bifurcation**

In Section 3, we have shown that the equilibrium *E*<sup>2</sup> = (*x*(2), *y*(2)) is unstable when *τ* = 0, and the equilibrium *E*<sup>3</sup> = (*x*(3), *y*(3)) is locally asymptotically stable when *τ* = 0. To reflect the actual situation, we focus on the delay from technological innovation to practical production. Therefore, we consider the time-delay *τ* as a bifurcation parameter and denote the critical value *<sup>τ</sup>* <sup>=</sup> *<sup>τ</sup><sup>c</sup>* <sup>=</sup> *<sup>τ</sup>*(*j*) *<sup>n</sup>* , where *<sup>τ</sup>*(*j*) *<sup>n</sup>* is given in (23). When *<sup>τ</sup>* <sup>=</sup> *<sup>τ</sup>*(*j*) *<sup>n</sup>* , characteristic Equation (21) has a pair of pure imaginary roots *λ* = ±*iω*. Therefore, system (3) undergoes a Hopf bifurcation at equilibrium *E*3. In this section, we derive the normal form of Hopf bifurcation for the system (3) by using the multiple time scales method given in [23,24].

In order to normalize the delay, we first re-scale the time *t* by using *t* → *t*/*τ*, then translate the equilibrium *E*<sup>3</sup> = (*x*(3), *y*(3)) to the origin, that is,

$$\begin{cases} \overline{\mathfrak{X}} = \mathfrak{x} - \mathfrak{x}^{(3)},\\ \overline{\mathfrak{y}} = \mathfrak{y} - \mathfrak{y}^{(3)}, \end{cases}$$

,

For convenience, we still use *<sup>x</sup>* and *<sup>y</sup>* to represent *<sup>x</sup>*( and *<sup>y</sup>*( respectively, so Equation (3) is transformed into:

$$\begin{cases} \frac{dx}{dt} = \tau (a\_1 x - c\_1 x^{(3)} y - c\_1 x y - c\_1 y^{(3)} x - k x (t - 1)),\\ \frac{dy}{dt} = -\tau (y + y^{(3)}) (b\_2 x + c\_2 y). \end{cases} \tag{25}$$

Equation (25) can also be written as:

$$Z(t) = \tau N\_1 Z(t) + \tau N\_2 Z(t-1) + \tau F(Z(t), Z(t-1)),\tag{26}$$

where

$$Z(t) = (\mathbf{x}(t), \mathbf{y}(t))^T,\\ Z(t-1) = (\mathbf{x}(t-1), \mathbf{y}(t-1))^T,\\ F(Z(t), Z(t-1)) = (\mathbf{x}(t-1), \mathbf{y}(t-1))^T,\\ F(Z(t), Z(t-1)) = (\mathbf{x}(t), \mathbf{y}(t-1))^T$$

and

$$N\_1 = \begin{pmatrix} a\_1 - c\_1 y^{(3)} & -c\_1 x^{(3)} \\ -b\_2 y^{(3)} & -c\_2 y^{(3)} \end{pmatrix}, \\ N\_2 = \begin{pmatrix} -k & 0 \\ 0 & 0 \end{pmatrix}.$$

We let *h* be eigenvector corresponding to eigenvalue *λ* = *iωτ* of Equation (26), and *h*∗ be the eigenvector corresponding to eigenvalue *λ* = −*iωτ* of adjoint matrix of Equation (26), satisfying

$$ = \overline{h^\*}^T h = 1.\tag{27}$$

By calculating, we have

$$h = (1, -\frac{b\_2 y^{(3)}}{i\omega + c\_2 y^{(3)}})^T,\\ h^\* = d(\frac{i\omega - c\_2 y^{(3)}}{c\_1 x^{(3)}}, 1)^T,\\ d = \frac{c\_1 x^{(3)}}{c\_1 x^{(3)} + b\_2 y^{(3)}}.\tag{28}$$

We treat the delay *<sup>τ</sup>* as the bifurcation parameter, let *<sup>τ</sup>* <sup>=</sup> *<sup>τ</sup><sup>c</sup>* <sup>+</sup> *εμ*, where *<sup>τ</sup><sup>c</sup>* <sup>=</sup> *<sup>τ</sup>*(*j*) *n* (*j* = 0, 1, 2, ···) is the Hopf bifurcation critical value, *μ* is perturbation parameter, *ε* is dimensionless scale parameter. Suppose system (26) undergoes a Hopf bifurcation from the trivial equilibrium at the critical point *τ* = *τc*, and then, by the MTS method, the solution of (26) is assumed as follows:

$$Z(\mathbf{t}) = Z(T\_{0\prime}, T\_{1\prime}, T\_{2\prime}, \dots) = \sum\_{k=1}^{+\infty} \varepsilon^k Z\_k(T\_{0\prime}, T\_{1\prime}, T\_{2\prime}, \dots),\tag{29}$$

where

$$\begin{aligned} Z(T\_{0\prime}, T\_{1\prime}, T\_{2\prime}, \cdots) &= (\mathfrak{x}(T\_{0\prime}, T\_{1\prime}, T\_{2\prime}, \cdots), \mathfrak{y}(T\_{0\prime}, T\_{1\prime}, T\_{2\prime}, \cdots))^T, \\ Z\_k(T\_{0\prime}, T\_{1\prime}, T\_{2\prime}, \cdots) &= (\mathfrak{x}\_k(T\_{0\prime}, T\_{1\prime}, T\_{2\prime}, \cdots), \mathfrak{y}\_k(T\_{0\prime}, T\_{1\prime}, T\_{2\prime}, \cdots))^T, \end{aligned}$$

and the derivative with regard to *t* is transformed into

$$\frac{d}{dt} = \frac{\partial}{\partial T\_0} + \varepsilon \frac{\partial}{\partial T\_1} + \varepsilon^2 \frac{\partial}{\partial T\_2} + \cdots = D\_0 + \varepsilon D\_1 + \varepsilon^2 D\_2 + \cdots \cdot \tag{30}$$

where *Di* is differential operator, and

$$D\_i = \frac{\partial}{\partial T\_i} \left( i = 0, 1, 2, 3, \cdots \right).$$

From (26), we have

$$Z(t) = \varepsilon D\_0 Z\_1 + \varepsilon^2 D\_1 Z\_1 + \varepsilon^3 D\_2 Z\_1 + \varepsilon^2 D\_0 Z\_2 + \varepsilon^3 D\_1 Z\_2 + \varepsilon^3 D\_0 Z\_3 + \dotsb \dots \tag{31}$$

We expand *<sup>x</sup>*(*T*<sup>0</sup> − 1,*ε*(*T*<sup>0</sup> − <sup>1</sup>),*ε*2(*T*<sup>0</sup> − <sup>1</sup>), ···) at *<sup>x</sup>*(*T*<sup>0</sup> − 1, *<sup>T</sup>*1, *<sup>T</sup>*2, ···) by the Taylor expansion, that is,

$$\mathbf{x}(t-1) = \varepsilon \mathbf{x}\_{1,\tau\_{\varepsilon}} + \varepsilon^2 \mathbf{x}\_{2,\tau\_{\varepsilon}} + \varepsilon^3 \mathbf{x}\_{3,\tau\_{\varepsilon}} - \varepsilon^2 D\_1 \mathbf{x}\_{1,\tau\_{\varepsilon}} - \varepsilon^3 D\_2 \mathbf{x}\_{1,\tau\_{\varepsilon}} - \varepsilon^3 D\_1 \mathbf{x}\_{2,\tau\_{\varepsilon}} + \dotsb,\tag{32}$$

where *xj*,*τ<sup>c</sup>* = *xj*(*T*<sup>0</sup> − 1, *T*1, *T*2, ···), *j* = 1, 2, 3 ··· .

We substitute Formulas (29)–(32) into Equation (26), then comparing the coefficients of *ε*,*ε*<sup>2</sup> and *ε*<sup>3</sup> on both sides of the equation, respectively. Then, we obtain the following expressions, respectively,

$$\begin{aligned} D\_0 \mathbf{x}\_1 - \mathfrak{r}\_\varepsilon (a\_1 \mathbf{x}\_1 - c\_1 \mathbf{x}^{(3)} y\_1 - c\_1 \mathbf{x}\_1 y^{(3)} - k \mathbf{x}\_{1, \tau\_\varepsilon}) &= \mathbf{0}, \\ D\_0 y\_1 + \mathfrak{r}\_\varepsilon y^{(3)} (b\_2 \mathbf{x}\_1 + c\_2 y\_1) &= \mathbf{0}. \end{aligned} \tag{33}$$

$$\begin{split} &D\_0 \mathbf{x}\_2 - \mathfrak{x}\_\varepsilon (a\_1 \mathbf{x}\_2 - c\_1 \mathbf{x}^{(3)} y\_2 - c\_1 \mathbf{x}\_2 y^{(3)} - k \mathbf{x}\_{2, \tau\_\varepsilon}) \\ &= \mu (a\_1 \mathbf{x}\_1 - c\_1 \mathbf{x}^{(3)} y\_1 - c\_1 \mathbf{x}\_1 y^{(3)} - k \mathbf{x}\_{1, \tau\_\varepsilon}) - \mathfrak{x}\_\varepsilon (c\_1 \mathbf{x}\_1 y\_1 - k D\_1 \mathbf{x}\_{1, \tau\_\varepsilon}) - D\_1 \mathbf{x}\_{1, \tau\_\varepsilon} \\ &D\_0 y\_2 + \mathfrak{x}\_\varepsilon y^{(3)} (b\_2 \mathbf{x}\_2 + c\_2 y\_2) \\ &= - \left(\mu y^{(3)} + \mathfrak{x}\_\varepsilon y\_1\right) (b\_2 \mathbf{x}\_1 + c\_2 y\_1) - D\_1 y\_1. \end{split} \tag{34}$$

$$\begin{split} &D\_0 \mathbf{x}\_3 - \mathbf{c}\_c (a\_1 \mathbf{x}\_3 - c\_1 \mathbf{x}^{(3)} y\_3 - c\_1 \mathbf{x}\_3 y^{(3)} - k \mathbf{x}\_{3, \tau\_c}) \\ &= \mu (a\_1 \mathbf{x}\_2 - c\_1 \mathbf{x}\_2 y\_1 - c\_1 \mathbf{x}^{(3)} y\_2 - c\_1 \mathbf{x}\_2 y^{(3)} + k \mathbf{x}\_{2, \tau\_c} + k D\_1 \mathbf{x}\_{1, \tau\_c}) - \mathbf{x}\_c (c\_1 \mathbf{x}\_2 y\_1 + c\_1 \mathbf{x}\_1 y\_2 - k \mathbf{x}\_{3, \tau\_c}) \\ &+ k D\_2 \mathbf{x}\_{1, \tau\_c} - k D\_1 \mathbf{x}\_{2, \tau\_c}) - D\_1 \mathbf{x}\_2 - D\_2 \mathbf{x}\_{1, \tau\_c} \\ &- D\_0 y\_3 + \tau\_c y^{(3)} (b\_2 \mathbf{x}\_3 + c\_2 y\_3) \\ &= -\tau\_c (y\_2 (b\_2 \mathbf{x} - 1 + c\_2 y\_1) + y\_1 (b\_2 \mathbf{x}\_2 + c\_2 y\_2)) - \mu (y\_1 (b\_2 \mathbf{x}\_1 + c\_2 y\_1) + y^{(3)} (b\_2 \mathbf{x}\_2 + c\_2 y\_2)) \\ &- D\_1 y\_2 - D\_2 y\_1. \end{split} \tag{35}$$

Equation (33) has the solution with following form,

$$Z\_1 = G h e^{i\omega \tau\_c T\_0} + \overline{G} h e^{-i\omega \tau\_c T\_0} \, \tag{36}$$

where *h* is given by (28). Equation (34) is a linear non-homogeneous equation, and the non-homogeneous equation has a solution if and only if a solvability condition is satisfied. That is, the right-hand side of (34) should be orthogonal to every solution of the adjoint homogeneous problem. Thus, we solve (36) into the right part of equation (34), and the coefficient vector of *eiωτcT*<sup>0</sup> is noted as *m*1, by < *h*∗, *m*<sup>1</sup> >= 0, so we can solve *<sup>∂</sup><sup>G</sup> ∂T*<sup>1</sup> , namely,

$$\frac{\partial G}{\partial T\_1} = M\mu G\_\prime \tag{37}$$

where

$$\begin{aligned} M &= \frac{c\_1 \mathbf{x}^{(3)} (y^{(3)} b\_2 + y^{(3)} c\_2 b\_2) - (c\_2 y^{(3)} + i\omega)(-a\_1 + c\_1 \mathbf{x}^{(3)} h\_2 + c\_1 y^{(3)} + k e^{-i\omega \tau\_\xi})}{(c\_2 y^{(3)} + i\omega)(-\tau\_\varepsilon k e^{-i\omega \tau\_\xi} + 1) - c\_1 \mathbf{x}^{(3)} h\_2}, \\ h\_2 &= -\frac{b\_2 y^{(3)}}{i\omega + c\_2 y^{(3)}} \end{aligned}$$

with

$$\begin{aligned} x^{(3)} &= \frac{c\_1 a\_2 + (k - a\_1) c\_2 - \sqrt{(c\_1 a\_2 + (k - a\_1) c\_2)^2 - 4c\_1 c\_2 b\_2 km}}{2c\_1 b\_2}, \\ y^{(3)} &= \frac{c\_1 a\_2 - (k - a\_1) c\_2 + \sqrt{(c\_1 a\_2 + (k - a\_1) c\_2)^2 - 4c\_1 c\_2 b\_2 km}}{2c\_1 c\_2}. \end{aligned}$$

We solve Equation (34), as *μ* is a small disturbance coefficient, and we only consider the influence of *μ* on low-order terms, thus, we obtain its solutions with following form:

$$\begin{aligned} x\_2 &= g\_1 e^{2i\omega \tau\_\epsilon T\_0} + c.c. + l\_1, \\ y\_2 &= g\_2 e^{2i\omega \tau\_\epsilon T\_0} + c.c. + l\_2. \end{aligned} \tag{38}$$

where *c*.*c*. stands for the complex conjugate of the preceding terms, then, we substitute solutions (38) into (34), and we obtain

$$\begin{aligned} g\_1 &= K\_1 S G^2, \; l\_1 = P\_1 R G \overline{G}, \\ g\_2 &= K\_2 S G^2, \; l\_2 = P\_2 R G \overline{G}, \end{aligned} \tag{39}$$

with

$$\begin{aligned} K\_1 &= (2i\omega + c\_2 y^{(3)})(-c\_1 h\_2) + c\_1 x^{(3)} (c\_2 h\_2^2 + b\_2 h\_2), \\ K\_2 &= (2i\omega + c\_1 y^{(3)} - a\_1 + k e^{-2i\omega \tau\_4})(-c\_2 h\_2^2 - b\_2 h\_2) + c\_1 b\_2 y^{(3)} h\_2, \\ S &= ((2i\omega + c\_2 y^{(3)})(2i\omega + c\_1 y^{(3)} - a\_1 + k e^{-2i\omega \tau\_4}) - c\_1 b\_2 y^{(3)} x^{(3)})^{-1}, \\ P\_1 &= -c\_1 \mathbf{x}^{(3)} (b\_2 (\overline{h\_2} + h\_2) + 2 c\_2 h\_2 \overline{h\_2}) + c\_2 c\_1 y^{(3)} (\overline{h\_2} + h\_2), \\ P\_2 &= -c\_1 b\_2 y^{(3)} (\overline{h\_2} + h\_2) + (c\_1 y^{(3)} + k - a\_1) (b\_2 (\overline{h\_2} + h\_2) + 2 c\_2 h\_2 \overline{h\_2}), \\ R &= (-c\_2 y^{(3)} (k - a\_1 + c\_1 y^{(3)}) + c\_1 b\_2 \mathbf{x}^{(3)} y^{(3)})^{-1}, \ h\_2 = -\frac{b\_2 y^{(3)}}{i\omega + c\_2 y^{(3)}}. \end{aligned} \tag{40}$$

Next, substituting solution (36) and (38) into (35), and with the coefficient vector of *eiωτcT*<sup>0</sup> noted as *m*2, by solvability condition, we have < *h*∗, *m*<sup>2</sup> >= 0. Note that *μ* is a disturbance parameter, and *μ*<sup>2</sup> has little influence for small unfolding parameter, thus, we can ignore the *μ*2*G* term, then *<sup>∂</sup><sup>G</sup> ∂T*<sup>2</sup> , can be solved to yield

$$\frac{\partial G}{\partial T\_2} = HG^2 \overline{G}\_\prime \tag{41}$$

where

$$\begin{split} H &= \frac{(\dot{\imath}\omega + c\_{2}y^{(3)}) \left(-\tau\_{c}c\_{1}\right) (K\_{1}S\overline{h\_{2}} + K\_{2}S + P\_{1}Rh\_{2} + P\_{2}R)}{(\dot{\imath}\omega + c\_{2}y^{(3)})(k\tau\_{c}e^{-i\omega\tau\_{c}} - 1) + c\_{1}\mathbf{x}^{(3)}h\_{2}} \\ &+ \frac{\tau\_{c}b\_{2}c\_{1}\mathbf{x}^{(3)}(K\_{1}S\overline{h\_{2}} + K\_{2}S + P\_{1}Rh\_{2} + P\_{2}R) + 2\tau\_{c}c\_{2}c\_{1}\mathbf{x}^{(3)}(K\_{2}S\overline{h\_{2}} + P\_{2}Rh\_{2})}{(\dot{\imath}\omega + c\_{2}y^{(3)})(k\tau\_{c}e^{-i\omega\tau\_{c}} - 1) + c\_{1}\mathbf{x}^{(3)}h\_{2}}, \\ h\_{2} &= -\frac{b\_{2}y^{(3)}}{i\omega + c\_{2}y^{(3)}}, \end{split}$$

with

$$\begin{aligned} x^{(3)} &= \frac{c\_1 a\_2 + (k - a\_1) c\_2 - \sqrt{(c\_1 a\_2 + (k - a\_1) c\_2)^2 - 4c\_1 c\_2 b\_2 km}}{2c\_1 b\_2} \\ y^{(3)} &= \frac{c\_1 a\_2 - (k - a\_1) c\_2 + \sqrt{(c\_1 a\_2 + (k - a\_1) c\_2)^2 - 4c\_1 c\_2 b\_2 km}}{2c\_1 c\_2} \end{aligned}$$

Let *G* → *G*/*ε*, we obtain the normal form of Hopf bifurcation of system (26) truncated at the cubic order terms:

$$
\dot{G} = M\mu G + HG^2 \overline{G},
\tag{42}
$$

,

.

where *M* is given in (37), and *H* given in (41).

With the polar coordinate *G* = *re*i*θ*, substituting that expression into Equation (42), we obtain the amplitude equation of Equation (42) on the center manifold as

$$\begin{cases} \dot{r} = \text{Re}(M)\mu r + \text{Re}(H)r^3, \\ \dot{\theta} = \text{Im}(M)\mu + \text{Im}(H)r^2. \end{cases} \tag{43}$$

**Theorem 4.** *When* Re(*M*)*<sup>μ</sup>* Re(*H*) < 0*, system* (26) *has periodic solutions.*

*(1) If* Re(*M*)*μ* < 0*, then the periodic solution reduced on the center manifold is unstable;*

*(2) If* Re(*M*)*μ* > 0*, then the periodic solution reduced on the center manifold is stable.*
