2.1. The formal series method of computing nondegenerate singular point quantities on center manifold

Considering the Jacobian matrix A at the origin of system (1) has a pair of purely imaginary eigenvalues and a negative one, then by certain nondegenerate transformation, the system (1) can be changed into the following system:

$$\begin{cases} \frac{\mathbf{dx}}{\mathbf{dt}} = -y + \sum\_{k+j+l=2}^{\infty} A\_{kjl} \mathbf{x}^k y^l u^l = \mathbf{X}(\mathbf{x}, y, u),\\ \frac{\mathbf{dy}}{\mathbf{dt}} = \mathbf{x} + \sum\_{k+j+l=2}^{\infty} B\_{kjl} \mathbf{x}^k y^l u^l = \mathbf{Y}(\mathbf{x}, y, u),\\ \frac{\mathbf{du}}{\mathbf{dt}} = -d\_0 u + \sum\_{k+j+l=2}^{\infty} \tilde{d}\_{kjl} \mathbf{x}^k y^j u^l = \tilde{\mathbf{U}}(\mathbf{x}, y, u) \end{cases} \tag{4}$$

where <sup>x</sup>; <sup>y</sup>; <sup>u</sup>; Akjl; Bkjl; <sup>~</sup>dkjl <sup>∈</sup> <sup>R</sup> <sup>ð</sup>k;j;l∈N<sup>Þ</sup> and <sup>d</sup><sup>0</sup> <sup>&</sup>gt; 0.

Here, we recall first the calculation method of the singular point quantities on center manifold for the above real three-dimensional nonlinear dynamical systems. By means of transformation

$$z = x + y\mathbf{i}, \quad w = x - y\mathbf{i}, \quad u = u, \quad T = \mathbf{i}t, \quad \mathbf{i} = \sqrt{-1} \tag{5}$$

system (4) is also transformed into the following complex system:

$$\begin{cases} \frac{\text{d}z}{\text{d}T} = z + \sum\_{k+j+l=2}^{\circ} a\_{kjl} z^k w^l u^l = Z(z, \varpi, u),\\ \frac{\text{d}w}{\text{d}T} = -\varpi - \sum\_{k+j+l=2}^{\circ} b\_{kjl} w^k z^j u^l = -\mathcal{W}(z, \varpi, u),\\ \frac{\text{d}u}{\text{d}T} = \text{id}u + \sum\_{k+j+l=2}^{\circ} d\_{kjl} z^k w^j u^l = \mathcal{U}(z, \varpi, u) \end{cases} \tag{6}$$

where z; w; T; akjl; bkjl; dkjl∈C ðk; j; l ∈ NÞ, the systems (4) and (6) are called concomitant.

Theorem 1 (see [16]). For system (6), using the program of term by term calculations, we can determine a formal power series:

$$F(z, \varpi, \boldsymbol{\mu}) = z\varpi + \sum\_{\alpha + \beta + \gamma = 3}^{\infty} c\_{\alpha\beta\gamma} z^{\alpha} \varpi^{\beta} \boldsymbol{\mu}^{\gamma} \tag{7}$$

such that

$$\frac{\partial F}{\partial T} = \frac{\partial F}{\partial z}Z - \frac{\partial F}{\partial y}\mathcal{W} + \frac{\partial F}{\partial u}\mathcal{U} = \sum\_{m=1}^{\circ} \mu\_m (zw)^{m+1} \tag{8}$$

where c<sup>110</sup> ¼ 1; c<sup>101</sup> ¼ c<sup>011</sup> ¼ c<sup>200</sup> ¼ c<sup>020</sup> ¼ 0; ckk<sup>0</sup> ¼ 0, k ¼ 2, 3, ⋯.

Definition 1. The μ<sup>m</sup> in the expression (8) is called the mth singular point quantity at the origin on center manifold of system (6) or (4), m ¼ 1; 2;⋯.

Theorem 2 (see [16, 34]). For the mth singular point quantity and the mth focal value at the origin on center manifold of system (4), i.e., μ<sup>m</sup> and v2mþ<sup>1</sup>; m ¼ 1; 2;⋯, we have the following relation:

$$\left(\upsilon\_{2m+1}(2\pi) = i\pi\mu\_m + i\pi\sum\_{k=1}^{m-1} \xi\_m^{(k)}\mu\_k\right. \tag{9}$$

where ξðk<sup>Þ</sup> <sup>m</sup> ðk ¼ 1; 2;⋯;m − 1Þ are polynomial functions of coefficients of system (6). Usually, it is called algebraic equivalence and written as v2mþ<sup>1</sup>~iπμm.

Based on the previous work in Ref. [16], we have developed the calculation method of the focal values on the center manifold for real four-dimensional nonlinear dynamical systems in Ref. [35]. In fact, here Theorem 1 can be generalized in the n-dimensional real systems as follows

$$\begin{cases} \frac{\text{dx}}{\text{dt}} = \text{-}y + \text{h.o.t.} = X(\text{x}, \text{y}, \text{u}),\\ \frac{\text{dy}}{\text{dt}} = \text{x} + \text{h.o.t.} = Y(\text{x}, \text{y}, \text{u}),\\ \frac{\text{du}\_{i}}{\text{dt}} = -d\_{i}u\_{i} + \text{h.o.t.} = \check{\mathcal{U}}\_{i}(\text{x}, \text{y}, \text{u}), \ i = 1, 2, \cdots, n-2 \end{cases} \tag{10}$$

where u ¼ ðu1;u2;⋯;un<sup>−</sup>2Þ, h.o.t denotes the terms in x;y;u1;u2;⋯;un<sup>−</sup><sup>2</sup> with orders greater than or equal to 2, and all di > 0.

By means of transformation of Eq. (5), system (10) can be transformed into the following complex system

$$\begin{cases} \frac{d\mathbf{z}}{dT} = \mathbf{z} + \sum\_{k+j+1=2}^{\circ\circ} a\_{kj} z^k w^j \mathbf{u}^1 = Z(z, w, \mathbf{u}),\\ \frac{d\mathbf{u}}{dT} = -\boldsymbol{\nu} - \sum\_{k+j+1=2}^{\circ\circ} b\_{kj} w^k z^j \mathbf{u}^1 = -\mathcal{W}(z, \boldsymbol{\nu}, \mathbf{u}),\\ \frac{d\mathbf{u}\_i}{dT} = \mathbf{i} d\_i \, u\_i + \sum\_{k+j+1=2}^{\circ\circ} d\_{kj} z^k w^j \mathbf{u}^1 = \mathcal{U}\_i(z, \boldsymbol{\nu}, \mathbf{u}), \ \mathbf{i} = 1, 2, \cdots, n-2 \end{cases} \tag{11}$$

where the subscript "kj1" denotes "kjl1⋯ln<sup>−</sup>2", <sup>u</sup><sup>1</sup> <sup>¼</sup> ul<sup>1</sup> <sup>1</sup> ul<sup>2</sup> <sup>2</sup> <sup>⋯</sup>uln<sup>−</sup><sup>2</sup> <sup>n</sup>−2, and l ¼ ∑ n−2 i¼1 li, all ui ∈ R, z; w; T; akj1; bkj1; dkj<sup>1</sup> ∈ C ðk; j; li ∈ NÞ, we call that system (10) and system (11) are concomitant.

Theorem 3. For system (11), using the program of term by term calculations, we can determine a formal power series:

$$F(z, w, \mathfrak{u}) = zw + \sum\_{\alpha + \beta + \ell = 3}^{\infty} c\_{\alpha\beta} \, \iota \, z^{\alpha} w^{\beta} \mathfrak{u}^{\ell} \tag{12}$$

such that

Theorem 1 (see [16]). For system (6), using the program of term by term calculations, we can

Definition 1. The μ<sup>m</sup> in the expression (8) is called the mth singular point quantity at the origin

Theorem 2 (see [16, 34]). For the mth singular point quantity and the mth focal value at the origin on

Based on the previous work in Ref. [16], we have developed the calculation method of the focal values on the center manifold for real four-dimensional nonlinear dynamical systems in Ref. [35]. In fact, here Theorem 1 can be generalized in the n-dimensional real systems as follows

<sup>d</sup><sup>t</sup> <sup>¼</sup> <sup>−</sup>diui <sup>þ</sup> <sup>h</sup>:o:t: <sup>¼</sup> <sup>U</sup><sup>~</sup> <sup>i</sup>ðx;y;uÞ, <sup>i</sup> <sup>¼</sup> <sup>1</sup>; <sup>2</sup>;⋯;n−<sup>2</sup>

where u ¼ ðu1;u2;⋯;un<sup>−</sup>2Þ, h.o.t denotes the terms in x;y;u1;u2;⋯;un<sup>−</sup><sup>2</sup> with orders greater than

By means of transformation of Eq. (5), system (10) can be transformed into the following

z; w; T; akj1; bkj1; dkj<sup>1</sup> ∈ C ðk; j; li ∈ NÞ, we call that system (10) and system (11) are concomitant.

<sup>u</sup><sup>1</sup> <sup>¼</sup> <sup>Z</sup>ðz;w;uÞ,

<sup>u</sup><sup>1</sup> <sup>¼</sup> <sup>−</sup>Wðz;w;uÞ,

<sup>u</sup><sup>1</sup> <sup>¼</sup> Uiðz;w;uÞ, <sup>i</sup> <sup>¼</sup> <sup>1</sup>; <sup>2</sup>;⋯;n−<sup>2</sup>

<sup>n</sup>−2, and l ¼ ∑

n−2 i¼1

<sup>1</sup> ul<sup>2</sup> <sup>2</sup> <sup>⋯</sup>uln<sup>−</sup><sup>2</sup>

<sup>d</sup><sup>t</sup> <sup>¼</sup> <sup>−</sup><sup>y</sup> <sup>þ</sup> <sup>h</sup>:o:t: <sup>¼</sup> <sup>X</sup>ðx;y;uÞ,

<sup>d</sup><sup>t</sup> <sup>¼</sup> <sup>x</sup> <sup>þ</sup> <sup>h</sup>:o:t: <sup>¼</sup> <sup>Y</sup>ðx;y;uÞ,

center manifold of system (4), i.e., μ<sup>m</sup> and v2mþ<sup>1</sup>; m ¼ 1; 2;⋯, we have the following relation:

v2mþ<sup>1</sup>ð2πÞ ¼ iπμ<sup>m</sup> þ iπ ∑

∞ αþβþγ¼3

> U ¼ ∑ ∞ m¼1

> > m−1 k¼1 ξðk<sup>Þ</sup>

<sup>m</sup> ðk ¼ 1; 2;⋯;m − 1Þ are polynomial functions of coefficients of system (6). Usually, it is called

μmðzwÞ

cαβγz<sup>α</sup>w<sup>β</sup>u<sup>γ</sup> (7)

<sup>m</sup>þ<sup>1</sup> (8)

<sup>m</sup> μ<sup>k</sup> (9)

(10)

(11)

li, all ui ∈ R,

Fðz; w; uÞ ¼ zw þ ∑

dF <sup>d</sup><sup>T</sup> <sup>¼</sup> <sup>∂</sup><sup>F</sup> ∂z Z− ∂F ∂y W þ ∂F ∂u

on center manifold of system (6) or (4), m ¼ 1; 2;⋯.

algebraic equivalence and written as v2mþ<sup>1</sup>~iπμm.

dx

8 >>>>><

>>>>>:

dz

8 >>>>><

>>>>>:

dw

dui

<sup>d</sup><sup>T</sup> <sup>¼</sup> <sup>z</sup> <sup>þ</sup> <sup>∑</sup><sup>∞</sup>

<sup>d</sup><sup>T</sup> <sup>¼</sup> <sup>−</sup>w−∑<sup>∞</sup>

<sup>d</sup><sup>T</sup> <sup>¼</sup> <sup>i</sup>di ui <sup>þ</sup> <sup>∑</sup><sup>∞</sup>

where the subscript "kj1" denotes "kjl1⋯ln<sup>−</sup>2", <sup>u</sup><sup>1</sup> <sup>¼</sup> ul<sup>1</sup>

kþjþ1¼2

kþjþl¼2

akj<sup>1</sup> zk wj

bkj1w<sup>k</sup> zj

> dkj<sup>1</sup> z<sup>k</sup> wj

kþjþ1¼2

or equal to 2, and all di > 0.

complex system

dy

dui

where c<sup>110</sup> ¼ 1; c<sup>101</sup> ¼ c<sup>011</sup> ¼ c<sup>200</sup> ¼ c<sup>020</sup> ¼ 0; ckk<sup>0</sup> ¼ 0, k ¼ 2, 3, ⋯.

determine a formal power series:

6 Manifolds - Current Research Areas

such that

where ξðk<sup>Þ</sup>

$$\frac{\partial F}{\partial T} = \frac{\partial F}{\partial z} Z - \frac{\partial F}{\partial y} W + \sum\_{i=1}^{n-2} \frac{\partial F}{\partial u\_i} lI\_i = \sum\_{m=1}^{\circ} \mu\_m (zw)^{m+1} \tag{13}$$

where the subscript "αβℓ" denotes "αβγ1⋯γn−2", <sup>u</sup><sup>ℓ</sup> <sup>¼</sup> <sup>u</sup><sup>γ</sup><sup>1</sup> <sup>1</sup> u γ2 <sup>2</sup> ⋯u γn−<sup>2</sup> <sup>n</sup>−<sup>2</sup> , and <sup>ℓ</sup> <sup>¼</sup> <sup>∑</sup> n−2 i¼1 γi , and more setting <sup>c</sup>αβ<sup>ℓ</sup> <sup>¼</sup> <sup>0</sup> with <sup>0</sup>≤<sup>α</sup> <sup>þ</sup> <sup>β</sup> <sup>þ</sup> <sup>ℓ</sup>≤<sup>2</sup> except for c<sup>110</sup> <sup>¼</sup> <sup>1</sup>, and ckk<sup>0</sup> <sup>¼</sup> <sup>0</sup> with k≥2.

Proof. It is very similar to the proving course of Theorem 1.3.1 in [16], by computing carefully and comparing the above power series with the two sides of (13), we can obtain the expression of μm.

Definition 2. The μ<sup>m</sup> in the expression (13) is called the mth singular point quantity at the origin on center manifold of system (11) or (10), m ¼ 1; 2;⋯.

Remark 1. Similar to Theorem 2, there exists a equivalence between μ<sup>m</sup> and v2mþ1, namely, if μ<sup>1</sup> ¼ μ<sup>2</sup> ¼ ⋯ ¼ μm−<sup>1</sup> ¼ 0; μm≠0, then v<sup>3</sup> ¼ v<sup>5</sup> ¼ ⋯ ¼ v2m−<sup>1</sup> ¼ 0; v2mþ<sup>1</sup> ¼ iπμm; m ¼ 1; 2;⋯, and vice versa.

Corollary 1. The origin of system (10) or (11) is a center restricted to the center manifold if and only if μ<sup>m</sup> ¼ 0 for all m.

Remark 2. From the relation given by Remark 1 and Corollary 1, the center-focus problem and Hopf bifurcation of equilibrium point restricted to the center manifold can be figured out by the calculation of singular point quantities for system (10).

#### 2.2. An example of four-dimensional system

Recently, the study of chaos has become a hot research topic, and the attention of many researchers is turning to 4D systems from 3D dynamical systems, for example, the authors of Ref. [36] investigated Hopf bifurcation of a 4D-hyoerchaotic system by applying the normal form theory in 2012, but its multiple Hopf bifurcation on the center manifold have not been considered. Here, we will investigate the system further by computing the singular point quantities of its equilibrium point, which takes the following form

$$\begin{cases}
\dot{\mathbf{x}}\_1 = a(\mathbf{x}\_2 - \mathbf{x}\_1) \\
\dot{\mathbf{x}}\_2 = c\mathbf{x}\_1 - \mathbf{x}\_2 + \mathbf{x}\_4 - \mathbf{x}\_1 \mathbf{x}\_3 \\
\dot{\mathbf{x}}\_3 = \mathbf{x}\_1 \mathbf{x}\_2 - b\mathbf{x}\_3 + e \mathbf{x}\_1^2 \\
\dot{\mathbf{x}}\_4 = -\mathbf{K} \mathbf{x}\_2
\end{cases} \tag{14}$$

where a; b; c;e; K∈R. Obviously, system (14) has only one isolated equilibrium: Oð0; 0; 0; 0Þ when K≠0. Therefore, we only need to consider O. The Jacobian matrix of system (14) at O is

$$A = \begin{pmatrix} -a & a & 0 & 0 \\ c & -1 & 0 & 1 \\ 0 & 0 & -b & 0 \\ 0 & -K & 0 & 0 \end{pmatrix}$$

with the characteristic equation:

$$(\lambda + b)(\lambda^3 + (a+1)\lambda^2 + (a \text{--} ac + K)\lambda + aK = 0. \tag{15}$$

To guarantee that A has a pair of purely imaginary eigenvalues �i ωðω > 0Þ and two negative real eigenvalues λ1;λ2, we let its characteristic equation take the form

$$(\lambda^2 + \alpha^2)(\lambda - \lambda\_1)(\lambda - \lambda\_2) = 0.1$$

Thus, we obtain the critical condition of Hopf bifurcation at O:

$$a^2(c-1) = a^2, \; K = a(a+1)(c-1), \; \lambda\_1 = -b, \; \lambda\_2 = -a-1 \tag{16}$$

where <sup>a</sup> <sup>&</sup>gt; <sup>−</sup>1; <sup>b</sup> <sup>&</sup>gt; <sup>0</sup>; <sup>c</sup> <sup>&</sup>gt; 1, namely, <sup>c</sup> <sup>¼</sup> <sup>a</sup>2þω<sup>2</sup> <sup>a</sup><sup>2</sup> ; <sup>K</sup> <sup>¼</sup> <sup>ð</sup>aþ1Þω<sup>2</sup> <sup>a</sup> . Under the conditions (16), one can find a nondegenerate matrix

$$P = \begin{pmatrix} -\frac{\text{i}a^2}{(a+1)(a+i\omega)\omega} & \frac{\text{i}a^2}{(a+1)(a-i\omega)\omega} & 0 & -\frac{a^2}{\omega^2} \\ -\frac{\text{i}a}{(a+1)\omega} & \frac{\text{i}a}{a\omega+\omega} & 0 & \frac{a}{\omega^2} \\ 0 & 0 & 1 & 0 \\ 1 & 1 & 0 & 1 \end{pmatrix}$$

such that

$$P^{-1}AP = \begin{pmatrix} \omega \mathbf{i} & 0 & 0 & 0 \\ 0 & -\omega \mathbf{i} & 0 & 0 \\ 0 & 0 & -b & 0 \\ 0 & 0 & 0 & -a-1 \end{pmatrix} \tag{17}$$

Namely, we can use the nondegenerate transformation and the time rescaling: T ¼ itω to make the system (14) become the following same form as the complex system (11) with n ¼ 4:

$$\begin{cases} \frac{d\mathbf{z}}{dT} = \mathbf{z} + \sum\_{k+l+l=n=2}^{2} a\_{ljk} z^k w^l u^l v^n = \mathbf{Z}(z, \omega v, \mathbf{u}, \mathbf{v}),\\ \frac{d\mathbf{u}}{dT} = -\mathbf{w} - \sum\_{k+j+l=n=2}^{2} b\_{ljk} w^k z^j u^l v^n = -\mathcal{W}(z, \omega v, \mathbf{u}, \mathbf{v}),\\ \frac{d\mathbf{u}}{dT} = \frac{\mathbf{b}\mathbf{i}}{\omega} + \sum\_{k+j+l=n=2}^{2} d\_{ljk} z^k w^l u^l v^n = \mathcal{U}(z, \omega v, \mathbf{u}, \mathbf{v}),\\ \frac{d\mathbf{v}}{dT} = \frac{(a+1)\mathbf{i}}{\omega} \mathbf{v} + \sum\_{k+j+l=n=2}^{2} e\_{ljk} z^k w^l u^l v^n = V(z, \omega v, \mathbf{u}, \mathbf{v}) \end{cases} \tag{18}$$

where u∈R, z; w; T∈C, and all akjln ¼ bkjln ¼ dkjln ¼ ekjln ¼ 0 except the following coefficients

$$\begin{split} a\_{0011} &= \frac{a^3 + a^2(1+\mathbf{i}\omega) + \mathbf{i}a\omega}{2\omega^2(a+\mathbf{i}\omega+1)}, \quad a\_{0110} = \frac{a(\omega-\mathbf{i}a)}{2\omega\left(a^2+a+\omega(\omega-\mathbf{i})\right)}, \\ b\_{ij\mathbf{i}a} &= \overline{a}\_{ij\mathbf{i}a} \quad (ikjl = 0011, 0110), \\ d\_{0002} &= \frac{\mathbf{i}a^3(1-a)e}{a^5}, \quad d\_{0101} = -\frac{a^4(2e+1) - a^3(1+\mathbf{i}\omega)}{(a+1)\omega^4(a-\mathbf{i}\omega)}, \\ d\_{0200} &= \frac{a^3\omega+\mathbf{i}a^4(e+1)}{(a+1)^2a^3(a-\mathbf{i}\omega)^2}, \quad d\_{1001} = \frac{a^4(2e+1) + a^3(\mathbf{i}\omega-1)}{(a+1)\omega^4(a+\mathbf{i}\omega)}, \\ d\_{1100} &= -\frac{2\mathbf{i}a^4(e+1)}{(a+1)^2a^3(a^2+\omega^2)}, \quad d\_{2000} = -\frac{a^3\omega-\mathbf{i}a^4(e+1)}{(a+1)^2a^3(a+\mathbf{i}\omega)^2}, \\ e\_{0011} &= -\frac{\mathbf{i}a(a+1)}{\omega(a^2+2a+\omega^2+1)}, \quad e\_{0110} = -\frac{a}{(a-\mathbf{i}\omega)(a^2+2a+\omega^2+1)}, \\ e\_{1010} &= \frac{a}{(a+\mathbf{i}\omega)(a^2+2a+\omega^2+1)} \end{split}$$

where akjln denotes the conjugate complex number of akjln. According to Theorem 3, we obtain the recursive formulas of cαβγ and μm.

Theorem 5. For system (18), setting cαβγλ ¼ 0 with 0≤α þ β þ γ þ λ≤2 except for c<sup>1100</sup> ¼ 1, and ckk<sup>00</sup> ¼ 0 with k≥2, we can derive successively and uniquely the terms of the following formal series (12) with n <sup>¼</sup> <sup>4</sup>, such that (13) with n <sup>¼</sup> <sup>4</sup> holds and if <sup>α</sup>≠<sup>β</sup> or <sup>α</sup> <sup>¼</sup> <sup>β</sup>; <sup>λ</sup><sup>2</sup> <sup>þ</sup> <sup>γ</sup><sup>2</sup>≠0, cαβγλ is determined by following recursive formula:

$$\begin{array}{l} \omega\\ c\_{a\beta\gamma} = \frac{\omega}{a(\alpha-\beta) + i(b\ \gamma+(a+1)\ \lambda)}\\ \{-d\_{2000}(1+\gamma)c[\alpha-2.\beta,\gamma+1,\lambda]-d\_{1100}(\gamma+1)c[\alpha-1,\beta-1,\gamma+1,\lambda]-e\_{010}\}\\ c\_{1010}(\lambda+1)c[\alpha-1,\beta,\gamma-1,\lambda+1]-d\_{1001}(\gamma+1)c[\alpha-1,\beta,\gamma+1,\lambda-1]+\\ b\_{0110}(\beta+1)c[\alpha-1,\beta+1,\gamma-1,\lambda]-d\_{0200}(\gamma+1)c[\alpha,\beta-2,\gamma+1,\lambda]-\\ e\_{0110}(\lambda+1)c[\alpha,\beta-1,\gamma-1,\lambda+1]-d\_{0101}(\gamma+1)c[\alpha,\beta-1,\gamma+1,\lambda-1]-e\_{010}\}\\ e\_{0011}\lambda c[\alpha,\beta,\gamma-1,\lambda]-d\_{0002}(\gamma+1)c[\alpha,\beta,\gamma+1,\lambda-2]+\\ b\_{0011}(\beta+1)c[\alpha,\beta+1,\gamma-1,\lambda-1]-a\_{0110}(\alpha+1)c[\alpha+1,\beta-1,\gamma-1,\lambda]-\\ a\_{0011}(\alpha+1)c[\alpha+1,\beta,\gamma-1,\lambda-1] \end{array} \tag{19}$$

and for any positive integer m; μ<sup>m</sup> is determined by following recursive formula:

$$\begin{aligned} \mu\_m &= d\_{2000}c[-2+m,m,1,0] \\ &+ d\_{1100}c[-1+m,-1+m,1,0] + d\_{0200}c[m,-2+m,1,0] \end{aligned} \tag{20}$$

and when α < 0 or β < 0 or γ < 0 or λ < 0 or α ¼ β; γ ¼ λ ¼ 0, we have let cαβγλ ¼ 0, and where each c½α;β;γ;λ� denotes cαβγλ .

By applying the above formulas in the Mathematica symbolic computation system, we figure out easily the first two singular point quantities of the origin of system (18):

$$\begin{aligned} \mu\_1 &= \text{i} \mathfrak{a} f\_1 \left[ \begin{smallmatrix} |a| \ b \ c \end{smallmatrix} \, \mathfrak{a} \, \begin{smallmatrix} a+1 \end{smallmatrix} \right. \mathfrak{d}\_0 \right]^{-1}, \\ \mu\_2 &= 108 \text{i} \mathfrak{a}^3 b^4 f\_2 \, f\_3^2 \, f\_4 \left[ |a| \ c^2 d\_0 \mathrm{d}\_1^{-2} \mathrm{d}\_2^{-4} \mathrm{d}\_3 \right]^{-1} \end{aligned} \tag{21}$$

where

A ¼

real eigenvalues λ1;λ2, we let its characteristic equation take the form

Thus, we obtain the critical condition of Hopf bifurcation at O:

<sup>−</sup> <sup>i</sup>a<sup>2</sup>

P<sup>−</sup><sup>1</sup> AP ¼

<sup>−</sup> <sup>i</sup><sup>a</sup> ða þ 1Þω

ða þ 1Þða þ iωÞω

0

BB@

kþjþlþn¼2

kþjþlþn¼2

kþjþlþn¼2

<sup>v</sup> <sup>þ</sup> <sup>∑</sup><sup>2</sup>

<sup>u</sup> <sup>þ</sup> <sup>∑</sup><sup>2</sup>

akjln zk wj ul

bkjlnw<sup>k</sup> zj ul

> dkjln zk wj ul

kþjþlþn¼2

where u∈R, z; w; T∈C, and all akjln ¼ bkjln ¼ dkjln ¼ ekjln ¼ 0 except the following coefficients

ekjln zk wj ul

where <sup>a</sup> <sup>&</sup>gt; <sup>−</sup>1; <sup>b</sup> <sup>&</sup>gt; <sup>0</sup>; <sup>c</sup> <sup>&</sup>gt; 1, namely, <sup>c</sup> <sup>¼</sup> <sup>a</sup>2þω<sup>2</sup>

0

BBBBBBB@

dz

8

>>>>>>>>><

>>>>>>>>>:

dw

du <sup>d</sup><sup>T</sup> <sup>¼</sup> <sup>b</sup><sup>i</sup> ω

dv

<sup>d</sup><sup>T</sup> <sup>¼</sup> <sup>z</sup> <sup>þ</sup> <sup>∑</sup><sup>2</sup>

<sup>d</sup><sup>T</sup> <sup>¼</sup> <sup>−</sup>w−∑<sup>2</sup>

<sup>d</sup><sup>T</sup> <sup>¼</sup> <sup>ð</sup><sup>a</sup> <sup>þ</sup> <sup>1</sup>Þ<sup>i</sup> ω

P ¼

find a nondegenerate matrix

such that

<sup>ð</sup>λ<sup>2</sup> <sup>þ</sup> <sup>ω</sup><sup>2</sup>

with the characteristic equation:

8 Manifolds - Current Research Areas

0

BBB@

−a a 0 0 c −1 01 0 0 −b 0 0 −K 0 0

To guarantee that A has a pair of purely imaginary eigenvalues �i ωðω > 0Þ and two negative

Þðλ−λ1Þðλ−λ2Þ ¼ 0:

<sup>a</sup><sup>2</sup> ; <sup>K</sup> <sup>¼</sup> <sup>ð</sup>aþ1Þω<sup>2</sup>

1

CCCA

<sup>ð</sup><sup>λ</sup> <sup>þ</sup> <sup>b</sup>Þðλ<sup>3</sup> þ ð<sup>a</sup> <sup>þ</sup> <sup>1</sup>Þλ<sup>2</sup> þ ða−ac <sup>þ</sup> <sup>K</sup>Þ<sup>λ</sup> <sup>þ</sup> aK <sup>¼</sup> <sup>0</sup>: (15)

<sup>a</sup><sup>2</sup>ðc−1Þ ¼ <sup>ω</sup><sup>2</sup>; <sup>K</sup> <sup>¼</sup> <sup>a</sup>ð<sup>a</sup> <sup>þ</sup> <sup>1</sup>Þðc−1), <sup>λ</sup><sup>1</sup> <sup>¼</sup> <sup>−</sup>b; <sup>λ</sup><sup>2</sup> <sup>¼</sup> <sup>−</sup>a−<sup>1</sup> (16)

ia2 ða þ 1Þða−iωÞω

> ia aω þ ω

0 0 10 1 1 01

> ωi 00 0 0 −ωi 0 0 0 0 −b 0 000 −a−1

Namely, we can use the nondegenerate transformation and the time rescaling: T ¼ itω to make the system (14) become the following same form as the complex system (11) with n ¼ 4:

<sup>a</sup> . Under the conditions (16), one can

1

CCCCCCCA

(17)

(18)

<sup>0</sup> <sup>−</sup> <sup>a</sup><sup>2</sup> ω2

<sup>0</sup> <sup>a</sup> ω2

1

CCA

vn <sup>¼</sup> <sup>Z</sup>ðz;w;u;vÞ,

vn <sup>¼</sup> <sup>−</sup>Wðz;w;u;vÞ,

vn <sup>¼</sup> <sup>U</sup>ðz;w;u;vÞ,

vn <sup>¼</sup> <sup>V</sup>ðz;w;u;v<sup>Þ</sup>

$$\begin{aligned} f\_1 &= 8a^3ce + 8a^3c - 8a^3e - 8a^2 - 2a^2bc + 2a^2be + 8a^2ce + 8a^2c + 8a^2c \\ &- 8a^2c - 8a^2 + ab^2c + 3ab^2e + 2ab^2 + 2abc - 2ab + 3b^2e + 3b^2, \\ f\_2 &= (2a+b+2)^3(2ae + 2a^2b)(e+1), \\ f\_3 &= 4a^2e + 4a^2 - 3abc - 2ab + 4ae + 4a + b, \\ f\_4 &= 8a^5c^2 - 16a^5c + 8a^5 - 2a^4bc^2 + 2a^4bc + 8a^4c^2 - 16a^4c + 8a^4 + 2a^3b^2c \\ &- 2a^3b^2 - 4a^3bc + 4a^3b - 5a^2b^3c + 4a^2b^3 + 2a^2b^2c \\ &- 2a^2b^2 - 2a^2bc + 2a^2b - 2ab^3 - b^3, \\ d\_0 &= (a^2c + 2a + 1)(4a^2c - 4a^2 + b^2)(c-1)^{3/2}, \\ d\_1 &= 8a^3c - 8a^2 - 2ab^2 + 8a^2c - 8a^2 + 3ab^2 + 3b^2, \\ d\_2 &= 8a^2c + 8a^2 - 2abc + 8a + b^2 + 2b, \\ d\_3 &= 9a^2c - 8a^2 + 2a + 1, \end{aligned}$$

and the above expression of μ<sup>2</sup> is obtained under the condition of μ<sup>1</sup> ¼ 0. From Remark 1 and the singular point quantities (21), we have

Theorem 6. For the flow on center manifold of the system (14), the first 2 focal values of the origin are as follow

$$
\upsilon\_3 = i\pi\mu\_1, \qquad \upsilon\_5 = i\pi\mu\_2 \tag{22}
$$

where the expression of v<sup>5</sup> is obtained under the condition of v<sup>3</sup> ¼ 0.

Remark 3. In contrast to the result and process in [36], one can easily see that our first quantity is basically consistent with its characteristic exponent of bifurcating periodic solutions, and our algorithm is easy to realize with computer algebra system due to the linear recursion formulas, and more convenient to investigate the multiple Hopf bifurcation on center manifold.

Considering its Hopf bifurcation form of Theorem 6, we have the following:

Theorem 7. At least two small limit cycles can be bifurcated from the origin of the 4D-hyoerchaotic system (14), which lie in the neighborhood of the origin restricted to the center manifold.

The rigorous proof of the above theorem is very similar to the previous ones in [14, 16], namely, by calculating the Jacobian determinant with respect to the functions v3; v<sup>5</sup> and its variables, which will not be given here.
