Cyclicity of rigid centres on centre manifolds of three-dimensional systems

Abstract

We work with polynomial three-dimensional rigid differential systems. Using the Lyapunov constants, we obtain lower bounds for the cyclicity of the known rigid centres on their centre manifolds. Moreover, we obtain an example of a quadratic rigid centre from which is possible to bifurcate 13 limit cycles, which is a new lower bound for three-dimensional quadratic systems.

1. Introduction

In 1900, the International Mathematical Congress took place in Paris. In this event, the acclaimed mathematician David Hilbert proposed a list of open problems among which is the problem of determining and locating the limit cycles for a given planar polynomial differential system. This is known as Hilbert's 16th problem and it is considered to be the hardest challenge within the qualitative theory of differential equations, since it remains open even for the simpler cases, i.e. quadratic polynomial vector fields (see [4, 13]). There are several weaker versions of this problem, for instance, determining the quantity $M(n)$ of limit cycles that bifurcate from an elementary centre or focus for polynomial vector fields of degree $n$ [6, 18, 23]. A singular point of a planar polynomial vector field is an elementary centre or focus if the eigenvalues of its Jacobian matrix evaluated at the singular point are purely imaginary. Any differential system associated to a vector field having an elementary centre or focus can be written, after the proper change of variables and time rescaling, in the form:

\begin{align*} \dot{x}& ={-}y+P(x,y),\\ \dot{y}& =x+Q(x,y), \end{align*}

where $P$ and $Q$ are polynomials with no linear or constant terms.

The standard method to estimate the quantity of limit cycles bifurcating from an elementary centre or focus is to compute the Lyapunov coefficients (or equivalently the focal coefficients) when small perturbations of the system are considered (see [2, 5, 18] and the references therein). Christopher [5] developed a simple approach to obtain lower bounds for the quantity of limit cycles bifurcating from an elementary centre or focus, i.e. the small amplitude limit cycles. The idea consists in considering the linear part of the power series expansion of the first $k$ Lyapunov coefficients. If those $k$ linear parts are linearly independent, then it is possible by small perturbations to obtain $k$ limit cycles, that is to say that the cyclicity of the considered singular point is at least $k$.

We consider three-dimensional differential systems with a Hopf singular point, that is, systems having a singular point for which the Jacobian matrix has a pair of purely imaginary eigenvalues and a non-zero real eigenvalue. These systems can be put in the following canonical form, by means of a linear change of variables and time rescaling:

(1.1)\begin{equation} \begin{aligned} \dot{x}={-}y+P(x,y,z),\\ \dot{y}=x+Q(x,y,z),\\ \dot{z}={-}\lambda z+R(x,y,z), \end{aligned} \end{equation}

where $P,Q,R$ are polynomials with no linear or constant terms and $\lambda \neq 0$. The centre problem and Hilbert's 16th can be naturally stated for vector fields (1.1), since there exists an invariant bidimensional $C^r$-manifold $W^c$ tangent to the $xy$-plane at the origin for every $r\geqslant 3$. This result is known as the centre manifold theorem (see [12, 21]).

Recently, these problems have attracted the attention of researchers in the qualitative theory of ordinary differential equations. In [7], the authors prove that the set of systems of form (1.1) having a centre on the local centre manifold at the origin, to a fixed value of $\lambda$, corresponds to a variety in the space of admissible coefficients. In [17] the authors solve the centre problem for the Lü system and Mahdi in [14] considers the centre problem for quadratic systems of the form (1.1) obtained from a third order ODE. The authors of [15] complete this study using a new hybrid symbolic-numerical approach.

Regarding the cyclicity problem, the results of Christopher are originally stated for planar systems, but they can be extended for system (1.1) (see [8]), since, for these systems, the Lyapunov coefficients are well-defined and play the same role as their bidimensional counterpart [7, 8]. In [9, 10, 19, 20, 22], the authors use these ideas to study the cyclicity of centres on centre manifolds of several families of three-dimensional systems.

In [16], the authors extended the concept of rigid systems and rigid centres for three-dimensional systems (1.1). In essence, system (1.1) has a rigid centre at the origin when its restriction to the centre manifold has a rigid centre at the origin. Moreover, in [16], the authors solved the centre problem for several families of rigid systems.

Our goal in this work is to obtain lower bounds for the quantity of small amplitude limit cycles which bifurcate from the rigid centres identified in [16]. The structure of this paper is as follows: In § 2, we present the preliminary results on the investigation of the cyclicity for three-dimensional systems as well as an overview of the concept of rigid systems studied in [16]. In § 3 we study the cyclicity of the rigid centres considered in [16] obtaining lower bounds for the cyclicity of each system under quadratic perturbations. Also in this section, we prove the main result of this work in which we present an example of a quadratic three-dimensional system from which is possible to obtain $13$ small limit cycles from quadratic perturbations. This is the highest known lower bound for quadratic three-dimensional systems.

2. Fundamental concepts

2.1 Lyapunov coefficients and limit cycles

We work with the following perturbation of system (1.1):

(2.1)\begin{equation} \begin{aligned} \dot{x} & ={-}y+P(x,y,z)+G_1(x,y,z),\\ \dot{y} & =x+Q(x,y,z)+G_2(x,y,z),\\ \dot{z} & ={-}\lambda z+R(x,y,z)+G_3(x,y,z), \end{aligned} \end{equation}

where $P,Q,R,G_1,G_2,G_3$ are polynomials with no constant or linear terms. Let $\Lambda =(\lambda _1,\ldots,\lambda _m)$ denote the perturbation parameters, i.e. the coefficients of $G_i$, $i=1,2,3$. We assume that for $\Lambda =0$, the origin is a centre on the centre manifold of system (2.1).

First we need to define a displacement map to search for limit cycles. Via the change of variables $x=\rho \cos \theta$, $y=\rho \sin \theta$ and $z=\rho \omega$, the solution curves of (2.1) are determined by the following equations:

(2.2)\begin{equation} \begin{aligned} \dfrac{{\rm d}\rho}{{\rm d}\theta} & =\dfrac{\psi(\rho,\theta,\Lambda)}{1+\phi(\rho,\theta,\Lambda)},\\ \dfrac{{\rm d}\omega}{{\rm d}\theta} & =\dfrac{-\mu\omega+\Omega(\rho,\theta,\Lambda)}{1+\phi(\rho,\theta,\Lambda)}. \end{aligned} \end{equation}

For each $(\rho _0,\omega _0)$ with sufficiently small $\|(\rho _0,\omega _0)\|$, let

\[ \varphi(\theta,\rho_0,\omega_0)=(\rho(\theta,\rho_0,\omega_0),\omega(\theta,\rho_0,\omega_0)), \]

be the solution of (2.2) with initial conditions $(\rho (0,\rho _0,\omega _0),\omega (0,\rho _0,\omega _0))=(\rho _0,\omega _0)$.

Definition 2.1 For system (1.1), the map

\[ {\rm d}(\rho_0,\omega_0)=\varphi(2\pi,\rho_0,\omega_0)-(\rho_0,\omega_0) \]

is called displacement map.

The displacement map ${\rm d}(\rho _0,\omega _0)=({\rm d}_1(\rho _0,\omega _0),{\rm d}_2(\rho _0,\omega _0))$ is a bidimensional map, unlike the planar counterpart. However, there exists a unique analytical function $\tilde {\omega }(\rho _0)$ defined in a neighbourhood $V$ of $\rho _0=0$ such that ${\rm d}_2(\rho _0,\tilde {\omega }(\rho _0))\equiv 0$ (see [3]). The analytical function $\mathbf {d}(\rho _0)={\rm d}_1(\rho _0,\tilde {\omega }(\rho _0))$, is called reduced displacement map. We write its power series as:

\[ \mathbf{d}(\rho_0)=v_1(\Lambda)\rho_0+v_2(\Lambda)\rho_0^2+v_3(\Lambda)\rho_0^3+v_4(\Lambda)\rho_0^4+\cdots. \]

It is known that the first nonzero coefficient $v_k(\Lambda )$ is the coefficient of an odd power of $\rho _0$ (see [3]). Thus, we define the coefficients $l_k=v_{2k+1}(\Lambda )$ as the Lyapunov coefficients. The isolated zeros of the reduced displacement map correspond to limit cycles of system (2.1). Furthermore, the origin is a centre on the centre manifold $W^c$ if and only if $l_k=0$ for all $k\in \mathbb {N}$.

For computational purposes, working with the Lyapunov coefficients is not very efficient. Thus, another method is applied to study the displacement function.

Following the steps of [1, 7, 8], we construct a formal series

\[ H(x,y,z)=x^2+y^2+\sum_{j+k+l\geqslant 3}p_{jkl}x^jy^kz^l, \]

with real coefficients $p_{jkl}$ yet to determine. Let $X$ denote the vector field associated to system (2.1). It is always possible to choose $p_{jkl}$ such that

(2.3)\begin{equation} XH=\sum_{k\geqslant2}L_{k-1}(x^2+y^2)^k. \end{equation}

Moreover, the quantities $L_{k}$ are rational functions on the parameters of system (2.1) including $\Lambda$ whose denominators depend only on $\lambda$ (see [8]). The coefficients $L_{k}$ are called focal coefficients for system (1.1). The focal coefficients and the Lyapunov coefficients have the following relationship: given a positive integer $k>2$, we have

\[ l_1=\pi L_{1}\;\text{ and }\;l_{k}=\pi L_{k}\quad \bmod \langle L_{1},L_{2},\ldots,L_{k-1}\rangle. \]

This result is proved in [8]. Since the computations involving the expressions given in (2.3) are algebraic, one can use symbolic mathematics software (e.g. Maple and Mathematica) to compute a large amount of focal coefficients in less computational time.

Once a sufficient amount of focal coefficients is computed, we can use the next result to estimate the cyclicity of the centre of (2.1) $\mid _{\Lambda =0}$.

Theorem 2.2 For system (2.1), if

\[ {\rm rank}\left[\dfrac{\partial (L_{1},L_{2},\ldots,L_{k})}{\partial (\lambda_1,\ldots,\lambda_m)}\right]_{\Lambda=0}=k, \]

then it is possible to choose parameters $\Lambda$ such that at least $k-1$ small amplitude limit cycles bifurcate from the origin.

The proof of this theorem can be found in [8]. It is a simple computational task to compute the matrix in theorem 2.2. In fact, since the perturbation of the centre in system (2.1) is given through the functions $G_i$, $i=1,2,3$, it is sufficient to compute the linear term of the Taylor expansion at $\Lambda =0$ of each focal coefficient. The idea of using the linear parts of the Lyapunov coefficients (or focal coefficients) to estimate the cyclicity of centres was first introduced by Christopher in [5] for the planar case.

It is known in the literature (see, for instance, [8]) that an additional limit-cycle can be obtained by making a perturbation of the trace of the linear part of system (1.1), i.e. the following perturbation:

(2.4)\begin{equation} \begin{aligned} \dot{x} & =\lambda_0x-y+P(x,y,z)+ G_1(x,y,z),\\ \dot{y} & =x+\lambda_0y+Q(x,y,z)+ G_2(x,y,z),\\ \dot{z} & ={-}\lambda z+R(x,y,z)+ G_3(x,y,z). \end{aligned} \end{equation}

In [11], the authors consider an extension of Christopher's ideas. Their approach is based in the next theorem, proved in [11], which consists in the study of the higher order developments of the Lyapunov coefficients.

Theorem 2.3 Suppose that system (2.1) satisfies the hypothesis of theorem 2.2 and the linear parts of the Lyapunov constants $l_{i}$, $i=k+1,\ldots,k+l$, are linear combination of the first k. Then, by a change of variables if necessary, we can assume that $l_i=u_i+O(|u|^{2})$, for $i=0,\ldots,k$, and the next Lyapunov coefficients are written as

\[ l_i=h_i(u)+O(|u|^{m+1}),\ \text{for}\ i=k+1,\ldots,k+l, \]

where $h_i$ are homogeneous polynomials of degree $m\geqslant 2$ and $u=(u_{k+1},\ldots, u_{k+l})$. If there exists a line $\eta$, in the parameter space, such that $h_i(\eta )=0,$ for $i=k+1,\ldots,k+l-1$, the hypersurfaces $h_i=0$ intersect transversally along $\eta$ for $i=k+1,\ldots, k+l-1$, and $h_{k+l}(\eta )\neq 0$, then there are perturbations of the centre which produce $k+l$ limit cycles.

Theorem 2.3 has been successfully applied to obtain new lower bounds for the cyclicity of several polynomial systems, for instance in papers [10, 11].

2.2 Rigid systems in $\mathbb {R}^3$

We exhibit a quick overview of rigid systems in $\mathbb {R}^3$, which are studied with more detail in [16].

Definition 2.4 (Rigid systems in $\mathbb {R}^3$)

We say that the three-dimensional system (1.1) is rigid when its restriction to a centre manifold is rigid, i.e. the restricted system has constant angular speed. Moreover, if the origin is a centre on the centre manifold, we say that the origin is a rigid centre.

In order to avoid restricting system (1.1) to a centre manifold, which is not effective computationally, the authors in [16] introduced a subclass of rigid systems.

Definition 2.5 We say that system (1.1) is rigid by cylindrical coordinates when, in cylindrical coordinates $(x,y,z)\to (\rho \cos \theta,\rho \sin \theta,z)$, the transformed system satisfies $\dot {\theta }=1$.

The orbits of rigid systems by cylindrical coordinates rotate around the $z$-axis with constant angular speed. The following result classify rigid systems (1.1).

Proposition 2.6 (Corollary 2.5 in [16])

Consider system (1.1). Then,

1. (a) it is a rigid system if and only if

   \[ \left(xQ(x,y,z)-yP(x,y,z)\right)\vert_{z=h(x,y)}\equiv 0, \]
   where $h(x,y)$ is the local expression of a local centre manifold$;$

2. (b) it is a rigid system by cylindrical coordinates if and only if $xQ(x,y,z)-yP(x,y,z)\equiv 0$.

3. Cyclicity of rigid centres in $\mathbb {R}^3$

As mentioned in § 2.1 the focal coefficients are rational functions on the parameters of system (1.1). For a fixed value of $\lambda$, these functions become polynomials which make the computations more manageable.

Consider the following rigid system:

(3.1)\begin{equation} \begin{aligned} \dot{x} & ={-}y+xF_n(x,y,z),\\ \dot{y} & =x+yF_n(x,y,z),\\ \dot{z} & ={-}z+R_m(x,y), \end{aligned} \end{equation}

where $F_n$ is a homogeneous polynomial of degree $n$ and $R_m$ is a homogeneous polynomial of degree $m$. The parameter values for which the above system has a rigid centre at the origin are called centre conditions and were obtained in [16, theorem 4.1]. For each centre condition we perform the perturbation (2.1), with $G_i,i=1,2,3$ being quadratic polynomials to study the cyclicity. More precisely, we obtained the following result.

Theorem 3.1 Consider system (3.1) with $F_n(x,y,z)=\sum _{j+k+l=n}a_{jkl}x^jy^kz^l$ and $R_m(x,y)=\sum _{j+k=m}b_{jk}x^jy^k$.

1. (a) For $n=1$ and $m=2$, the two centre conditions are $\{a_{001}=0\}$ and $R_2\equiv 0$. The cyclicity of the centre is at least $6$ for the first centre condition and at least $4$ for the last;

2. (b) For $n=1$ and $m=3$, the two centre conditions are $\{a_{001}=0\}$ and $R_3\equiv 0$. The cyclicity of the centre is at least $7$ for the first centre condition and at least $4$ for the last;

3. (c) For $n=1$ and $m=4$, the two centre conditions are $\{a_{001}=0\}$ and $R_4\equiv 0$. The cyclicity of the centre is at least $8$ for the first centre condition and at least $4$ for the last;

4. (d) For $n=2$ and $m=2$, with the centre condition $\{a_{200}=-a_{020}, a_{101}=a_{011}=a_{002}=0\}$, the cyclicity of the centre is at least $4$.

Proof. Consider the following perturbation of system (3.1):

\begin{align*} \dot{x}& ={-}y+xF_n(x,y,z)+G_1(x,y,z),\\ \dot{y}& =x+yF_n(x,y,z)+G_2(x,y,z),\\ \dot{z}& ={-}z+R_m(x,y)+G_3(x,y,z), \end{align*}

with $G_i$, $i=1,2,3$ being homogeneous quadratic polynomials. For each centre condition described above, we compute the first 12 focal coefficients of the above perturbed system and the rank of their linear parts with respect to the coefficients of $G_i$, the perturbation parameters. Hence, applying theorem 2.2 and performing a perturbation of the trace, the result holds. For instance, for the first centre condition in statement (a), the rank of the linear part of the focal coefficients is 6. For the other statements, we have:

1. (a) Case $n=1,m=2$:

   • • $a_{001}=0$: rank $6$;

   • • $R_2\equiv 0$: rank $4$.

2. (b) Case $n=1,m=3$:

   • • $a_{001}=0$: rank $7$;

   • • $R_3\equiv 0$: rank $4$.

3. (c) Case $n=1,m=4$:

   • • $a_{001}=0$: rank $8$;

   • • $R_4\equiv 0$: rank $4$.

4. (d) Case $n=2,m=2$:

   • • $a_{200}=-a_{020}, a_{101}=a_{011}=a_{002}=0$: rank $4$.

By theorem 2.2, the result holds.

Remark 3.2 For case $n=2$ and $m=2$ in statement (d) of theorem 3.1, there is another centre condition given by $\{a_{200}=-a_{020}, b_{20}=b_{11}=b_{02}=0\}$. For this centre condition, we were not able to obtain lower bounds for the cyclicity using quadratic perturbations since the linear parts of the focal coefficients are all null. The same situation occurs for the case $n=2$ and $m=3$ with its respective centre conditions described in [16, theorem 4.1]. We believe that it is possible to obtain better lower bounds than those obtained in theorem 3.1 by using higher order perturbations. For instance, for the case $n=2$ and $m=3$ with centre condition $\{a_{200}=-a_{020}, a_{101}=a_{011}=a_{002}=0\}$, using cubic perturbations, the cyclicity is at least $15$. Here, we chose quadratic perturbations to simplify the computations since, even with cubic perturbations, the computation of the linear parts of the focal coefficients is computationally a hard task.

The rigid system

(3.2)\begin{equation} \begin{aligned} \dot{x} & ={-}y+xF_n(x,y,z),\\ \dot{y} & =x+yF_n(x,y,z),\\ \dot{z} & ={-}z+R_m(x,y,z), \end{aligned} \end{equation}

where $F_n$ is a homogeneous polynomial of degree $n$ and $R_m$ is a homogeneous polynomial of degree $m$ has also been considered in [16, theorems 4.3 and 4.4] where they obtained centre conditions for particular values of $m$ and $n$.

Theorem 3.3 Consider system (3.2) with $F_n(x,y,z)=\sum _{j+k+l=n}a_{jkl}x^jy^kz^l$ and $R_m(x,y,z)=\sum _{j+k+l=m}b_{jkl}x^jy^kz^l$. For $n=1$ and $m=2$, the known centre conditions are$:$

1. (a) $\{a_{001}=0\};$

2. (b) $R_2(x,y,0)\equiv 0;$

3. (c) $\{a_{100}=a_{010}=0$, $b_{101}=b_{011}=0$, $b_{002}=2a_{001}$, $b_{020}=-b_{200}\};$

4. (d) $\{b_{101}=2a_{100}+a_{010}$, $b_{020}=-b_{200}$, $b_{011}=-a_{100}+2a_{010}, b_{002}=2a_{001}\};$

5. (e) $\{b_{101}=2a_{100}$, $b_{020}=-b_{200}$, $b_{011}=2a_{010}$, $\eqalign{& b_{002}=2a_{001}, \cr & (a_{100}^2-a_{010}^2)(b_{200}-b_{110})+a_{100}a_{010}(b_{110}+4b_{200})=0\}}$;

6. (f) $\eqalign{& \left \{b_{200}=\dfrac {a_{100}a_{010}(b_{002}-a_{001})}{a_{001}^2}, b_{110}=\dfrac {(a_{100}^2-a_{010}^2)(a_{001}-b_{002})}{a_{001}^2},\right. \cr & \left. b_{101}=\dfrac {b_{002}(a_{100} + a_{010})-a_{010}a_{001}}{a_{001}}, b_{020}=\dfrac {a_{100}a_{010}(a_{001}-b_{002})}{a_{001}^2}, \right. \cr & \left. b_{011}=\dfrac {b_{002}(a_{010}-a_{100}) + a_{100}a_{001}}{a_{001}}\right \}}$.

The cyclicity of the centre at the origin of system (3.2) is at least $8$ for condition $($a$)$, at least $5$ for condition $($b$)$, at least $3$ for condition $($c$)$, at least $9$ for condition $($d$)$, at least $10$ for condition $($e$)$ and at least $11$ for condition $($f$)$.

Theorem 3.4 Consider system (3.2) with $F_n(x,y,z)=\sum _{j+k+l=n}a_{jkl}x^jy^kz^l$ and $R_m(x,y,z)=\sum _{j+k+l=m}b_{jkl}x^jy^kz^l$. For $n=m=2$ and $F_2=R_2$, the two centre conditions are $\{a_{200}=a_{110}=a_{020}=0\}$ and $\{a_{020}=-a_{200}, a_{101}=a_{011}=a_{002}=0\}$. The cyclicity of the centre is at least 3 for the first centre condition and at least 2 for the last.

Proofs of theorems 3.3 and 3.4. Proofs of theorems 3.3 and 3.4

The proofs are analogous to that of theorem 3.1. More precisely, for each centre condition, we compute the first 12 focal coefficients and using theorem 2.2 and performing a perturbation of the trace, we estimate the cyclicity by computing the rank of their linear parts.

The following rigid system

(3.3)\begin{equation} \begin{aligned} & \dot{x}={-}y+xF(z),\\ & \dot{y}=x+yF(z),\\ & \dot{z}={-}z+R_2(x,y,z), \end{aligned} \end{equation}

where $F(z)=\sum _{i=1}^{9}a_iz^i$ and $R_2$ is a homogeneous polynomial of degree $2$ was studied in [16, theorem 5.1].

Theorem 3.5 Consider system (3.3) with $F(z)=\sum _{i=1}^{9}a_iz^i$ and $R_2(x,y,z)=\sum _{j+k+l=2}b_{jkl}x^jy^kz^l$. For the three centre conditions: $F(z)\equiv 0$, $R_2(x,y,0)\equiv 0$ and $\{a_2=a_4=a_5=\cdots =a_9=0, b_{101}=b_{011}=a_3=0$,$b_{002}=2a_1, b_{020}=-b_{200}\}$, the cyclicity of the centre is at least $5$, $3$ and $3$, respectively for each centre condition.

Proof. Analogous to the proof of the previous theorems. For the first and third centre conditions, we compute the first 12 focal coefficients as well as the rank of their linear parts. For the second, we compute the first 10 focal coefficients. The conclusion follows from theorem 2.2 and by performing a perturbation of the trace.

Finally, we consider the following rigid system:

(3.4)\begin{equation} \begin{aligned} \dot{x} & ={-}y+xF(z),\\ \dot{y} & =x+yF(z),\\ \dot{z} & ={-}z+R_3(x,z), \end{aligned} \end{equation}

where $F(z)=a_1z+a_2z^2+a_3z^3$ and $R_3$ is a homogeneous polynomial of degree $3$ which was studied in [16, theorem 5.2].

Theorem 3.6 Consider system (3.4) with $R_m(x,z)=\sum _{j+l=3}b_{jl}x^jz^l$. For the two centre conditions: $\{b_{30}=0\}$ and $\{a_2=3a_1^2$, $a_3=0$, $b_{21}=b_{12}=0, b_{03}=9a_1^2\}$, the cyclicity of the centre is at least $3$ for the first centre condition and at least 14 for the last.

Proof. We proceed analogously to the previous theorems. For the first centre condition, we compute the first 12 focal coefficients as well as the rank of their linear part. For the second centre condition, we compute 17 focal coefficients. Once again, the conclusion follows from theorem 2.2 and by performing a perturbation of the trace.

Remark 3.7 In [16, theorem 5.2], the authors also present the centre condition for system (3.4) given by $F(z)\equiv 0$. Since under quadratic perturbations, the linear parts of the focal coefficients are all null, we were not able to obtain a lower bound for the cyclicity.

In [19], the authors obtained a lower bound for the number of small limit cycles bifurcating from the centre for polynomial three-dimensional systems with a Hopf singular point for degrees 2 through 5. For quadratic systems, their lower bound was 11 limit cycles. In [10], the authors applied the study of the higher order terms of the focal coefficients to find an example of a quadratic system that bifurcates, under quadratic perturbations, 12 limit cycles from the centre. Using the same technique for system (3.2) with $n=1$ and $m=2$, which consists in applying theorem 2.3, we were able to surpass this lower bound, as we prove the next result.

Theorem 3.8 Consider the rigid system

(3.5)\begin{equation} \begin{aligned} & \dot{x}={-}y+xF_1(x,y,z),\\ & \dot{y}=x+yF_1(x,y,z),\\ & \dot{z}={-}z+R_2(x,y,z), \end{aligned} \end{equation}

with $F_1(x,y,z)=x+2y+3z$ and $R_2(x,y,z)=\tfrac {2}{3}x^2+xy+4xz-\tfrac {2}{3}y^2+3yz+6z^2$. The origin is a rigid centre on the centre manifold and at least 13 limit cycles bifurcate from it under quadratic perturbations.

Proof. Note that system (3.5) satisfies the centre condition described in item (d) from theorem 3.3 and therefore the origin is a centre on the centre manifold. Considering a quadratic perturbation (2.1) with

\[ G_1=\sum_{j+k+l=2}a_{jkl}x^jy^kz^l,\quad G_2=\sum_{j+k+l=2}b_{jkl}x^jy^kz^l,\quad G_3=\sum_{j+k+l=2}c_{jkl}x^jy^kz^l \]

for system (3.5), we compute the first 13 focal coefficients. The rank of their linear part is 9 and by theorem 2.2, it is possible to obtain 9 limit cycles bifurcating from the origin using the perturbation parameters $(a_{jkl},b_{jkl},c_{jkl})$ as well as a perturbation of the trace (2.4). To simplify the computations, we assign $c_{011}=c_{020}=c_{101}=c_{110}=c_{200}=0$ since these parameters do not play a role when we apply theorems 2.2 and 2.3 to obtain a higher estimate for the cyclicity. Due to the size of the expressions, we only exhibit the linear part of $L_1,\ldots, L_3$:

\begin{align*} L_1^1& =\frac{22}{45}a_{{011}}+4a_{{020}}-{\frac{4}{45}} a_{{101}}+\frac{2}{3}a_{{110}}+\frac{4}{3}a_{{200 }}+\frac{4}{45}b_{{011}}-\frac{2}{3}b_{{020}}\\ & \quad+\frac{22}{45}b_{{101}}-\frac{4}{3}b_{{110}}-2b_{{200}},\\ L_2^1& =\frac{3713}{25}b_{{200}}+\frac{4}{3}c_{{002}}-\frac{2}{9}a_{{002}}-\frac{39\,937}{1125}a_{{011}}-\frac{7297}{25}a_{{020}}+\frac{6409}{1125}a_{{101}}\\ & \quad-\frac{3389}{75}a_{{110}}-\frac{1463}{15}a_{{200}}-\frac {9409}{1125}b_{{011}}+\frac{3839}{75}b_{{020}}-\frac{37\,687}{1125}b_{{101}}+\frac{1433}{15}b_{{110}},\\ L_3^1& =\frac{44\,084}{1785}a_{{002}}+\frac{28\,460\,414}{7875}a_{{011}}+\frac{442\,099\,341}{14\,875}a_{{020}}-\frac{76\,527\,541}{133\,875}a_{{101}}- \frac{5168}{35}c_{{002}}\\ & \quad+\frac{203\,415\,874}{44\,625}a_{{110}}+\frac{444\,002\,837}{44\,625}a_{{200}}+\frac{5707}{1785}b_{{002}}+\frac{116\,062\,741}{133\,875}b_{{011}}\\ & \quad-\frac{233\,857\,774}{44\,625}b_{{020}}+ \frac{26\,669\,714}{7875}b_{{101}}-\frac{433\,328\,537}{44\,625}b_{{110}}-\frac{225\,465\,802}{14\,875}b_{{200}}, \end{align*}

where $L_k^j$ denotes the homogeneous part of degree $j$ of $L_k$. We now turn to theorem 2.3. After a suitable change of variables in the perturbation parameters $(a_{002}, a_{011}, a_{020}, a_{101}, a_{110}, a_{200},$ $b_{002}, b_{011}, b_{020})$, we can write

\[ L_i^1=u_i,\quad i=1,\ldots,9. \]

Considering the new variables $\Lambda =(u_1,\ldots,u_9,b_{101}, b_{110}, b_{200}, c_{002})$, we find a solution $\Lambda =s$ for the equations $L_{10}^2(\Lambda )=L_{11}^2(\Lambda )=L_{12}^2(\Lambda )=0$ such that $L_{13}^2(s)\neq 0$. Since

\begin{align*} L_{10}^2& =\!{-}{\frac{42463\ldots00071}{62328\ldots85264}}b_{{101}}b_{{200}}\,{-}\,{\frac{40320\ldots92721}{31164\ldots42632}}b_{{110}}b_{{101}}\,{-}\,{\frac{65360\ldots49715}{93493\dots27896}}b_{{101}}c_{{002}}\\ & \quad-{\frac{49690\ldots18261}{12119\ldots05468}}b_{{110}}b_{{200}}-{\frac{91271\ldots35547}{81806\ldots36909}}b_{{101}}^{2}-{\frac{46335\ldots13063}{36358\ldots16404}}b_{{110}}^{2}\\ & \quad+{\frac{58144\ldots10977}{14543\ldots65616}}b_{{200}}^{2}-{\frac{27934\ldots30905}{34627\ldots15848}}c_{{002}}b_{{200}}-{\frac{20443\ldots06675}{13754\ldots68284}}c_{{002}}b_{{110}},\\ L_{11}^2& ={\frac{20939\ldots26951}{30054\ldots46160}}b_{{101}}b_{{200}}+{\frac{19880\ldots86461}{15027\ldots23080}}b_{{110}}b_{{101}}+{\frac{64454\ldots32783}{90164\ldots33848}}b_{{101}}c_{{002}}\\ & \quad+{\frac{40002\ldots58387}{95411\ldots04640}}b_{{110}}b_{{200}}+{\frac{25716\ldots92939}{22541\ldots84620}}{b_{{101}}}^{2}+{\frac{65259\ldots07273}{50091\ldots74360}}{b_{{110}}}^{2}\\ & \quad-{\frac{16384\ldots09829}{40072\ldots94880}}{b_{{200}}}^{2}+{\frac{33058\ldots21677}{40072\ldots59488}}c_{{002}}b_{{200}}+{\frac{57083\ldots74365}{37899\ldots60712}}c_{{002}}b_{{110}},\\ L_{12}^2& =\!{-}{\frac{81825\ldots04719}{11696\ldots02500}}b_{{101}}b_{{200}}\,{-}\,{\frac{14372\ldots72644}{10818\ldots73125}}b_{{110}}b_{{101}}\,{-}\,{\frac{93191\ldots38049}{12982\dots67750}}b_{{101}}c_{{002}}\\ & \quad-{\frac{16194\ldots76093}{38467\ldots60000}}b_{{110}}b_{{200}}-{\frac{40196\ldots79299}{35088\ldots07500}}{b_{{101}}}^{2}-{\frac{37739\ldots55031}{28850\ldots95000}}{b_{{110}}}^{2}\\ & \quad+{\frac{94762\ldots70553}{23080\ldots60000}}{b_{{200}}}^{2}-{\frac{19118\ldots35489}{23080\ldots76000}}c_{{002}}b_{{200}}-{\frac{65905\ldots19663}{43657\ldots69800}}c_{{002}}b_{{110}}, \end{align*}

a solution $\Lambda =s$ is given by:

\begin{align*} b_{101}& =1,\\ b_{200}& =\frac{\alpha}{3},\\ b_{110}& ={\frac{99478\ldots56893{\alpha}^{2}- 15595\ldots75578\alpha-61109\ldots90912}{74402\ldots67872{\alpha}^{2}+13857\ldots14244\alpha-17563\ldots30264}},\\ c_{002}& ={-}\frac{\begin{array}{l}14326\ldots21024{\alpha}^{4}+41672\ldots01633{\alpha}^{3}-54197\ldots40400{\alpha}^{2}\\ \quad- 12687\ldots07756\alpha+19077\ldots08656\end{array}} {\begin{array}{l}65412\ldots11040{\alpha}^{3}+29236\ldots40360{\alpha}^{2}\\\quad +16155\ldots29280\alpha-40172\ldots14560\end{array}}, \end{align*}

and

\[ L_{13}^2(s)=\frac{60274\ldots32675{\alpha}^{2}-33282\ldots47954\alpha-87785\ldots71976}{19996\ldots30400}, \]

where $\alpha$ is a real root of the cubic polynomial

\[ 37499\ldots25384x^3 + 97562\ldots80503x^2 - 13970\ldots27210x - 27339\ldots94856. \]

Note that $L_{13}^2(s)\neq 0$. Moreover, since the Jacobian determinant

\begin{align*} & \det\dfrac{\partial (L_{10}^2(s),L_{11}^2(s),L_{12}^2(s))}{\partial (b_{110}, b_{200}, c_{002})}\\& \quad ={\frac{15316\ldots73503{\alpha}^{2}- 11623\ldots71962\alpha-78417\ldots92584}{20165\ldots60960}}, \end{align*}

is non-zero, we have that the hypersurfaces $L_{i}^2=0, i=10,11,12$ intersect transversally along $s$. Therefore, by theorem 2.3, there are perturbations of the centre at the origin which yield $13$ limit cycles.