Bifurcation of small limit cycles in cubic integrable systems using higher-order analysis
Yun Tian Pei Yu




































    theoremTheoremlemma[theorem]LemmarmkRemark=20cm
        = 0.2in
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        J. Differential Equation
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        ^aDepartment of Mathematics, Shanghai Normal University, 
        Shanghai, 200234, P. R. China 
        ^bDepartment of Applied Mathematics, Western University, 
        London, Ontario, Canada N6A 5B7
        In this paper, we present a method of higher-order analysis on bifurcation 
        of small limit cycles around an elementary center of integrable systems 
        under perturbations. This method is equivalent to higher-order Melinikov 
        function approach used for studying bifurcation of limit cycles 
        around a center but simpler.
        Attention is focused on planar cubic polynomial systems 
        and particularly it is shown that the system studied
        by H. Żoła̧dek in the article 
         Eleven small limit cycles in a cubic vector field 
        (Nonlinearity 8, 843–860, 1995) can indeed 
        have eleven limit cycles under 
        perturbations at least up to 7th order. Moreover, 
        the pattern of numbers of limit cycles produced near 
        the center is discussed up to 39th-order perturbations,
        and no more than eleven limit cycles are found.
        
        
        
        
        Bifurcation of limit cycles; Higher-order analysis; Darboux integral; 
        Focus value.
        
        
        
        
         34C07  34C23
        
        
        
        
        
        
        § INTRODUCTION
        
        
        
        Bifurcation theory of limit cycles is important for both theoretical 
        development of qualitative analysis
        and applications in solving real problems. It is closely related 
        to the well-known Hilbert's 16th problem <cit.>,
        whose second part asks for the upper bound, 
        called Hilbert number H(n), on the number of limit cycles 
        that the following system,
        
            dx/dt = P_n(x,y),     dy/dt = Q_n(x,y),
        
        can have, where P_n(x,y) and Q_n(x,y)
        represent n^ th-degree polynomials in x and y.
        This problem has motivated many mathematicians 
        and researchers in other disciplines to 
        develop mathematical theories and methodologies in the areas of 
        differential equations and dynamical systems. 
        However, this problem has not been completely solved even for 
        quadratic systems since Hilbert proposed the problem in 
        the Second Congress of World Mathematicians in 1900. 
        The maximal number of limit cycles obtained for some quadratic systems is 
        4 <cit.>. 
        However, whether H(2)=4 is still open.
        For cubic polynomial systems, many
        results have been obtained on the lower bound of the number of limit cycles.
        So far, the best result for cubic systems is
        H(3) ≥ 13 <cit.>. Note that the 13 limit cycles
        obtained in <cit.>
        are distributed around several singular points.
         
        When the problem is restricted to consider the maximum number of 
        small-amplitude limit cycles, denoted by M(n), bifurcating from a 
        focus or a center in system (<ref>),
        one of the best-known results is M(2)=3, which was obtained by
        Bautin in 1952 <cit.>.
        For n=3, a number of results in this research direction have been obtained. 
        So far the best result for the number of small limit cycles around a 
        focus is 9 <cit.>, 
        and that around a center is 
        12 <cit.>. 
        
        One of powerful tools used for analyzing local bifurcation of limit cycles 
        around a focus or a center
        is normal form theory (e.g., see <cit.>). 
        Suppose system (<ref>) has an elementary focus
        or an elementary center at the origin.
        With the computation methods using computer algebra systems
        (e.g., see <cit.>),
        we obtain the normal form expressed in polar coordinates as
        
            [    dr/dt = r  ( v_0 + v_1   r^2 + v_2   r^4 
            + ⋯ + v_k   r^2k + ⋯),; d θ/dt 
            = ω_c + τ_0 + τ_1   r^2 + τ_2   r^4 + ⋯ 
            + τ_k   r^2k + ⋯ , ]
        
        where r and θ represent the amplitude and phase of motion,
        respectively.  v_k   (k=0,1,2, ⋯) is called the kth-order
        focus value. v_0 and τ_0 are obtained from linear analysis.
        The first equation of (<ref>) can be used for studying
        bifurcation and stability of limit cycles, while the
        second equation can be used to determine the frequency of the bifurcating
        periodic motion. Moreover, the coefficients τ_j can be used to
        determine the order or critical periods of a center
        (when v_j=0,  j ≥ 0).
        
        
        A particular attention has been paid to 
        near-integrable polynomial systems, described in the form of 
        
            [   dx/dt = M^-1(x,y,μ)H_y(x,y,μ) 
            + ε  p(x,y,ε,δ),; dy/dt = -M^-1(x,y,μ) H_x(x,y,μ) 
            + ε  q(x,y,ε,δ), ]
        
        where 0 < ε≪ 1,
        μ and δ are vector parameters;
        H(x,y,μ) is an analytic function in x, y and μ;
        p (x,y,ε,δ) and q(x,y,ε,δ) are
        polynomials in x and y,
        and analytic in δ and ε.
        M(x,y,μ) is an integrating factor of the unperturbed 
        system (<ref>)|_ε=0.
        
        Suppose the unperturbed system (<ref>)|_ε=0 
        has an elementary center.
        Then, considering limit cycles bifurcation in system (<ref>) 
        around the center, we may use the normal form theory to 
        obtain the first equation of (<ref>) as follows: 
        
            [ dr/dt 
            = r [v_0(ε) + v_1(ε)r^2 
            + v_2(ε)r^4 + ⋯ +v_i(ε)r^2i+⋯], ]
         
        where 
        
            v_i(ε) = ∑_k=1^∞ε^k V_ik, 
               i=0,1,2, …,
         
        in which V_ik denotes the ith ε^k-order focus value, 
        and will be used throughout this paper. 
        Note that v_i(ε)=O(ε)
        since the unperturbed system (<ref>)|_ε=0 
        is an integrable system.
        Further, because system (<ref>) is analytic in ε,
        we can rearrange the terms in (<ref>), and obtain 
        
            dr/dt = V_1(r) ε + V_2(r) ε^2 + ⋯ 
            + V_k(r) ε^k + ⋯ ,
         
        where 
        
            V_k(r) = ∑_i=0^∞ V_ik  r^2i+1,    k=1,2, ….
        
        
        Similarly, for the normal form of system (<ref>) 
        we have the θ differential equation, given by 
        
            dθ/dt = T_0(r) + O(ε),
         
        with T_0(0)≠0, and thus 
        
            dr/dθ=V_1(r) ε + V_2(r) ε^2 
            + ⋯ + V_k(r) ε^k + ⋯/T_0(r) + O(ε).
        
        Assume the solution r(θ,ρ,ε) of (<ref>), satisfying
        the initial condition r(0,ρ,ε)=ρ, is given 
        in the form of 
        
            r(θ,ρ,ε)=r_0(θ,ρ)+r_1(θ,ρ)ε
            +r_2(θ,ρ)ε^2 + ⋯ 
            + r_k(θ,ρ)ε^k + ⋯,
        
        with 0<ρ≪1. Then, r_0(0,ρ)=ρ and r_k(0,ρ)=0, 
        for k≥1.
        
        If there exists a positive integer K such that V_k(r)≡0, 1≤ k< K,
        and V_K(r)≢0, then it follows from (<ref>) that 
        
            r_0(θ,ρ)=ρ,   r_k(θ,ρ)=0,
               1≤ k< K,     r_K(θ,ρ)=V_K(ρ)/T_0(ρ)θ.
        
        Thus, the displacement function d(ρ) of system (<ref>) 
        can be written as
        
            d(ρ)=r(2π,ρ,ε)-ρ=2πV_K(ρ)/T_0(ρ)ε^K+O(ε^K+1).
        
        Therefore, if we want to determine the number of small-amplitude limit cycles 
        bifurcating from the center in system (<ref>), we only need to study 
        the number of isolated zeros of V_K(ρ) for 0<ρ≪ 1, 
        and have to obtain the expression of the first non-zero coefficient
        V_K(r) in (<ref>) by computing V_iK, for i≥0.  
        
        The above discussions show that the basic idea of using focus values 
        is actually the same as that of the Melnikov function method.
        Using H(x,y)=h to parameterize the section (i.e. the 
        Poincaré map),
        we obtain the displacement function of (<ref>), given by
        
            d(h)= M_1(h)ε + M_2(h)ε^2 + ⋯ +
            M_k(h)ε^k + ⋯,
        
        where
        
            M_1(h) = ∮_H(x,y,μ)=h M(x,y,μ) 
            [ q(x,y,0,δ)   dx - p(x,y,0,δ)   dy ],
        
        evaluated along closed orbits H(x,y,μ)=h for h∈(h_1,h_2).
        Then, we can study the first non-zero Melnikov function
        M_k(h) in (<ref>) to determine the number of 
        limit cycles in system (<ref>). 
        In the following, we remark on the comparison of the Melnikov function 
        method and the method of normal forms (or focus values). 
        
        
           (1)
        Let H=h,  0<h-h_1≪1 define closed orbits around 
        the center of system (<ref>)|_ε=0.
        It is easy to see that for any integer K≥1,
        equation (<ref>) holds
        if and only if
        M_k(h)≡0, 1≤ k<K and M_K(h)≢0 in (<ref>).
        Moreover, V_K(ρ) for 0<ρ≪1 and M_K(h) for 0<h-h_1 ≪ 1 
        have the same maximum number of isolated zeros. 
        
        
           (2)
        As we can see, V_k(r) can be obtained by the computation of normal forms 
        or focus values.
        
        
           (3)
        In particular, when the original system is not a Hamiltonian system but
        an integrable system, then even computing the coefficients of the
        first-order Melnikov function is much more involved than the computation
        of using the method of normal forms. 
        
        
           (4)
        However, the method of normal forms (or focus values) is restricted to Hopf and 
        generalized Hopf bifurcations,
        while the Melnikov function method can also be applied to study bifurcation
        of limit cycles from homoclinic/heteroclinic loops or any closed orbits.
        
        
        When we apply the method of normal form computation,
        some unnecessary perturbation parameters are involved 
        in the computation of high-order focus values,
        which could be extremely computation demanding (in both time and memory), 
        and makes it much more difficult to solve the problem. 
        Meanwhile, before we use the first non-zero coefficient 
        V_K(r) in (<ref>) to find limit cycles,
        we need to prove V_k(r)≡0, 1≤ k < K.
        The unnecessary parameters involved could greatly 
        increase the difficulty of proving that.
        
        In this paper, without loss of limit cycles,
        we introduce a linear transformation to eliminate unnecessary 
        parameters from system (<ref>). 
        With less parameters in (<ref>), we can use the approximation of 
        first integrals to prove V_k(r)≡0.
        The idea will be illuminated in Section 2.
        
        We will apply our method to study the bifurcation of small-amplitude limit 
        cycles in the system 
        
            [                           dx/dt = a + 5/2  x + xy + x^3 
            + ∑_k=1^n ε^k p_k(x,y),; dy/dt = -2 a x + 2 y - 3 x^2 + 4 y^2 - a x^3 + 6 x^2 y 
            +  ∑_k=1^n ε^k q_k(x,y), ]
         
        where 
        
            p_k(x,y) = a_00k + ∑_i+j=1^3 a_ijk  x^i y^j,      
            q_k(x,y) = b_00k + ∑_i+j=1^3 b_ijk  x^i y^j,
         
        in which a_ijk and b_ijk are ε^kth-order coefficients 
        (parameters). 
        The unperturbed system (<ref>)|_ε=0 has a rational 
        Darboux integral <cit.>, 
        
            H_0=f_1^5/f_2^4=(x^4+4x^2+4y)^5/(x^5+5x^3+5xy+5x/2+a)^4,
        
        with the integrating factor M=20f_1^4f_2^-5.
        It can be shown that for a<-2^5/4, system (<ref>)|_ε=0
        has a center at E_0 =(- a/2,- a^2+2/4). 
        The system (<ref>)|_ε=0 was proposed 
        in <cit.>, and it was claimed that this system 
        could have 11 limit 
        cycles around the center by studying the second-order Melnikov 
        function. Later, Yu and Han applied the normal form computation method and 
        got only 9 limit cycles around E_0 <cit.>
        by analyzing the ε- and ε^2-order focus values. 
        Recently, it has been shown <cit.> 
        that errors are made in <cit.> for choosing 12 integrals as 
        the basis of the linear space of corresponding Melnikov functions of system 
        (<ref>)|_ε=0.
        In fact, among the 12 chosen integrals, 
        two of them can be expressed as linear combinations of the other ten integrals,
        and therefore only 9 limit cycles can exist, agreeing with that shown in <cit.>.   
        
        The rest of the paper is organized as follows. 
        In the next section, we consider system (<ref>), and 
        construct a transformation to reduce the number of perturbation 
        parameters, which greatly simplifies the analysis in the following section. 
        Section 3 is devoted to the computation of higher 
        ε^k-order focus values and the existence of 
        11 limit cycles in system (<ref>), which needs 
        computing at least ε^7-order focus values.
        Finally, conclusion is drawn in Section 4. 
        
         
        
        
        
        
        § PRELIMINARIES
        
        
        
        The method of focus values (or normal forms) is one of important and powerful 
        tools for the study of small-amplitude limit cycles generated from Hopf 
        bifurcation. In general, a sufficient number of focus values would be needed 
        if one wants to find more small-amplitude limit cycles.
        One main challenge is that the computation of focus values becomes
        more and more difficult as the order of focus values goes up.
        That is why computer algebra systems such as Maple and Mathematica have 
        been used for computing the focus values 
        to improve the computational efficiency (e.g. see <cit.>).
        Another approach is 
        to eliminate certain parameters from the system, 
        which is the method we shall develop here for near-integrable systems.
        
        
        In most studies of near-integrable systems,
        full perturbations like those polynomials p(x,y,ε,δ) 
        and q(x,y,ε,δ) given in system (<ref>) are considered. 
        The parameter vector δ usually represents the coefficients
        in p and q. When normal forms are used to study small limit cycles,
        it is easy to get and solve the focus values of ε order
        (coefficients in V_1(r)), because they are linear functions of the 
        system parameters, namely the coefficients in p(x,y,0,δ) 
        and q(x,y,0,δ).
        For the ε^k-order focus values (coefficients in V_k(r)), 
        more parameters would be involved in the computation.
        One can observe that some parameters are not necessary 
        for obtaining the maximum number of limit cycles, 
        and they only increase the difficulty in finding limit cycles.
        
        When  the first n functions V_k(r) in (<ref>), 
        1≤ k≤ n
        are applied to studying 
        bifurcation of limit cycles, in order to 
        remove unnecessary parameters 
        without reducing the number of limit cycles,  
        we may use the following transformation: 
        
            {   x→ x+e_1(ε)x+e_2(ε)y+e_3(ε),
               y→ y+e_4(ε)x+e_5(ε)y+e_6(ε),
               t→ t+e_7(ε)t,
               μ→μ+e_8(ε),
            .
        
        where   
        
            e_i(ε)=e_i1ε+e_i2ε^2+⋯ 
            +e_inε^n,   i=1,⋯,8.
        
        Note that (<ref>)|_ε=0 is an identity map.
        Thus, (<ref>) keeps the unperturbed system of (<ref>) unchanged.
        Furthermore, the new system obtained by using (<ref>) can
        be still written in the same form of (<ref>).
        So we only need to find proper e_i(ε)'s to get 
        simpler perturbation functions without loss of generality.
        
        
        To illustrate how to obtain e_i(ε), 
        we take system (<ref>) as an example.
        The coefficients a_ijk and b_ijk in (<ref>) are the parameters.
        Substituting the transformation (<ref>) 
        into system (<ref>) yields
        
            [                            dx/dt = a + 5/2  x + xy + x^3 
            + ∑_k=1^n ε^k p̃_k(x,y)+o(ε^n),; dy/dt = -2 a x + 2 y - 3 x^2 + 4 y^2 - a x^3 + 6 x^2 y 
            +  ∑_k=1^n ε^k q̃_k(x,y) +o(ε^n), ]
         
        where
        
            p̃_k(x,y) 
            = ã_00k + ∑_i+j=1^3 ã_ijk  x^i y^j,     q̃_k(x,y) 
            = b̃_00k + ∑_i+j=1^3 b̃_ijk  x^i y^j.
         
        
        
        Obviously, the coefficients ã_ijk and b̃_ijk 
        in (<ref>) are linear in e_mk, m=1,…,8.
        Let E_k=(e_1k,e_2k,⋯,e_8k)^T.
        For any 1≤ k≤ n, ã_ijk and 
        b̃_ijk can be written in the form of 
        
            ã_ijk = A_ijE_k +η_ijk,  b̃_ijk = B_ijE_k +ζ_ijk,
        
        where A_ij and B_ij are 1×8 matrices, 
        and η_ijk and ζ_ijk, given by 
        
            η_ijk = η_ijk(E_1,⋯,E_k-1,
            a_ml1,⋯,a_mlk, b_ml1,⋯,b_mlk),
               ζ_ijk = ζ_ijk(E_1,⋯,E_k-1,
            a_ml1,⋯,a_mlk, b_ml1,⋯,b_mlk) ,
        
        are polynomials in e_ml, 1≤ l  ≤ k-1, 
        and the coefficients in the perturbation functions (<ref>). 
        
        Note that A_ij and B_ij are not dependent on k.
        We hope that we can find some proper values for e_ik 
        to make some of the coefficients 
        ã_ijk and b̃_ijk vanish or satisfy some conditions,
        so that the computation of the focus values would become easier.
        For instance, we can choose for 1≤ k≤ n, 
        
            [                                                                         ã_10k= ã_01k= ã_20k= 
            ã_11k= ã_02k= ã_30k= 0,;                                                            ã_p_k≜p̃_k(-a/2,-a^2+4/4)=0,  ã_q_k≜q̃_k(-a/2,-a^2+4/4)=0. ]
        
        The last two equations in (<ref>) keep the equilibrium 
        of system (<ref>) in a neighborhood of 
        E_0 with radius o(ε^n).
        A direct computation yields 
        
            [                                                      ã_10k =      2ae_2k+e_6k+5/2e_7k+η_10k,  ã_01k = 1/2e_2k+e_3k+η_01k,;                                                      ã_20k =            3e_2k+3e_3k+e_4k+η_20k,  ã_11k = e_5k+e_7k+η_11k,;                                                      ã_02k =               -3e_2k+η_02k,  ã_30k = 2e_1k+ae_1k+e_7k+η_30k,;                                                      ã_p_k =   -1/4a(4+a^2)e_1k 
            - 1/8(4+a^2)(2+a^2)e_2k + 1/4(4+a^2)e_3k;                                                                  +1/4a^2e_4k + 1/8a(2+a^2)e_5k 
            - 1/2ae_6k + e_8k +η̃_k ,;                                                      ã_q_k =          -1/8a^2(16+3 a^2)e_1k 
            - 1/16a(16+3 a^2)(2+a^2)e_2k;                                                              + 1/4a(16+3a^2)e_3k
            +1/4a(4+a^2)e_4k + 1/8(4+a^2)(2+a^2)e_5k;                                                                                   -1/4(4+a^2)e_6k
            + 1/8a(a^2+8)e_8k+ζ̃_k,;                                                              ]
        
        where η̃_k and ζ̃_k are also functions 
        in η_ijl and ζ_ijl with 1≤ l  ≤ k-1, 
        respectively.
        
        Because
        
            [∂(ã_10k,ã_01k,ã_20k,
            ã_11k, ã_02k,ã_30k,ã_p_k,ã_q_k)/∂(e_1k,e_2k,e_3k,e_4k,e_5k,e_6k,e_7k,e_8k)]
            =3/4(32-a^4)<0
        
        for a<-2^-5/4,
        we can solve (<ref>) for e_mk to obtain 
        
            e_mk=e_mk(η_10k,η_01k,η_20k,η_11k,η_02k,
            η_30k,η̃_k,ζ̃_k),   
            1≤ m≤ 8,
        
        which can be rewritten by using (<ref>) as
        
            e_mk=ẽ_mk(E_1,⋯,E_k-1,
            a_ij1,⋯,a_ijk, b_ij1,⋯,b_ijk).
        
        Note that e_m1 only depends on a_ij1 and b_ij1.
        Therefore, for all 1≤ m≤ 8, 1≤ k≤ n, e_mk can 
        be expressed as a polynomial
        in a_ijl and b_ijl, 1≤ l≤ k.
        In other words, (<ref>) has solutions for all 1≤ k≤ n.
        
        Thus, without loss of generality, we assume that (<ref>)
        takes the following form, 
        
            [                                                                p_k(x,y) =                           a_00k + 
            a_21k x^2 y + a_12k x y^2 + a_03k y^3,;                                                                q_k(x,y) = b_00k 
            + b_10k x   + b_01k y
                      +b_20k x^2 + b_11k x y + b_02k y^2;                                                                                                 +b_30k x^3 + b_21k x^2 y + b_12k x y^2 + b_03k y^3, ]
         
        with
        
            [                                                                                          a_00k=                             1/64 (a^2+2)
            [ (a^2+2)^2  a_03k + 2 a (a^2+2)  a_12k +4 a^2 a_21k],;                                                                                         b_00k = 1/64{ 8 a^3 b_30k+ 16 a (2 b_10k- a b_20k)
                   + 4 (a^2 + 2) (4 b_01k- 2 a b_11k+ a^2 b_21k);                                                                                                                                             - (a^2+2)^2 [ 4 b_02k -2 a b_12k -(a^2+2) b_03k]
            }. ]
         
        
        
        
        As mentioned in Section 1,
        to find limit cycles around E_0 in system (<ref>),
        we apply normal form theory to compute the focus values and then solve 
        the multivariate polynomial equations based on the focus values. 
        Particularly, we have
        
            [                                                   b_01k = 1/16[ 4 a  (2 b_11k-a b_21k) - (a^2+2)^2 (a_12k +3 b_03k);                                                                         +  4 (a^2+2) (2 b_02k -a a_21k -a b_12k) ], ]
         
        solved from the zeroth-order focus value V_0k=0, where  
        
            [                                                                     V_0k = 1/32{ 
            16  b_01k - 4 a  (2 b_11k-a b_21k) + (a^2+2)^2 (a_12k +3 b_03k);                                                                  [1.0ex]                              -  4 (a^2+2) (2 b_02k -a a_21k -a b_12k) ]. ]
         
        Higher-order focus values are relatively complex, and we shall
        study them in Section 3. 
        
        When we want to use focus values V_iK in V_K(r), i=0,1,2,⋯,
        to study limit cycles, we first need to show V_k(r)≡0, 
        1≤ k<K, or 
        dr/dt  = O(ε^K) in (<ref>). In order to prove this, 
        we use the approximation of first integrals, and claim that 
        if there exists an analytic function H_K(x,y,ε) such that
        
            (M^-1 H_y+ε p) ∂ H_K/∂ x
            + (- M^-1 H_x+ε q)∂ H_K/∂ y
            = O(ε^K),
         
        then dr/dt=O(ε^K).
        This claim can be easily proved 
        by using the closed contour H_K=h as the parameter
        to express the displacement function.
        
        
        Usually, like that considered in <cit.> 
        the method of focus values is used only to prove
        how many limit cycles around the equilibrium point
        that system (<ref>) can have. 
        Combining it with the approximation of first integrals,
        we can obtain
        the maximal number of small limit cycles for parameters 
        in a neighborhood of critical conditions.
        Furthermore, if the focus values are linear functions in  
        parameters, we have a global result as follows. 
        
        
        Consider system (<ref>) and assume V_k(r)≡0, 1≤ k<K. 
        Suppose that for an integer m≥ 1,
        each V_iK, 0≤ i<m is linear in δ,
        and further the following two conditions hold:
        
        
         (i)[ ∂(V_0K,⋯,V_m-1,K)/∂(δ_1,⋯,δ_m)] =m,
        
        
         (ii)V_K(r)≡0,  if V_iK=0, i=0,1,⋯,m-1. 
        
        
         
        Then, for any given N>0, there exist ε_0>0 and a 
        neighborhood V of the origin such that system (<ref>) 
        has at most m-1 limit cycles in V for 
        0<|ε|<ε_0, |δ|≤ N. Moreover, m-1 limit 
        cycles can appear in an arbitrary neighborhood of the origin for some 
        values of (ε,δ). 
        
        
        The above theorem can be proved following the proof given 
        in <cit.> with a minor modification. So the proof is omitted here. 
        
        
        
        
        
        
        § HIGHER-ORDER ANALYSIS LEADING TO 11 LIMIT CYCLES 
        IN SYSTEM (<REF>)
        
        
        
        In this section, we focus on system (<ref>)
        and show that it can have 11 limit cycles by using 
        perturbations at least up to 7th order. In the following,  
        we will use the transformed system (<ref>) with the 
        simplified perturbations given in (<ref>) for the analysis. 
        
        In order to compute the focus values of this system, we first shift 
        the equilibrium of system (<ref>), 
        (x,y)=(-a/2+o(ε^n),
        -a^2+2/4+o(ε^n)) 
        to the origin and then use a computer algebra system
        and software package (e.g., the Maple program in <cit.>) to 
        obtain the focus values in terms of the parameters 
        a, a_ijk and b_ijk. 
        We shall give detailed analysis for the first few lower-order 
        focus values, 
        and then summarize the results obtained from 
        higher-order analysis. 
        
        For convenience, define the vectors:  
        
            [                                                                           W_k^8                                                        = (V_1k, V_2k, ⋯, V_8k),;                                                                           W_k^9                                                        = (V_1k, V_2k, ⋯, V_9k),;                                                                          W_k^10                                                       = (V_1k, V_2k, ⋯, V_10k),;                                                                           S_k^8                  =  
            (b_10k, b_20k, b_11k, b_02k, b_30k, b_21k, b_12k, 
            b_03k),;                                                                           S_k^9          =   
            (b_10k, b_20k, b_11k, b_02k, b_30k, b_21k, b_12k, 
            b_03k, a_03k),;                                                                          S_k^10 =  
            (b_10k, b_20k, b_11k, b_02k, b_30k, b_21k, b_12k, 
            b_03k, a_03k, a_12(3m)), ]
         
        where in S_k^10, k=7m for Case (A) and k=13m for Case (B) 
        (m ≥ 1, integer) to be considered in Sections 3.4 and 3.5;  
        and the determinants:  
        
            [                                                                            _k^8 = [ ∂ W_k^8/∂ S_k^8 ],
              _k^9 = [ ∂ W_k^9/∂ S_k^9],   _k^10 = [ ∂ W_k^10/∂ S_k^10]; ]
         
        and the functions:  
        
            [                                                                                       F_1 =                                             - 373423834799904305184768/5  a^36 (a^4-32)^8 ,;                                                                                       F_2 =                                         3013505105717894236809449177088/5  a^45 (a^4-32)^9,;                                                                                       F_3 =                                  - 57397219210893210316046010501071634432/a^55 (a^4-32)^10,;                                                                                       F_4 =                           - 279638476916415193342384256641414767487418/a^66 ( a^4 - 32)^11,;                                                                                       G_1 =                                                               - 258237837/32  a^9 (a^4-32),;                                                                                       G_2 =                            23476167/64  a^11 (a^4-32)^2 
                  (57697 a^4-35728 a^2-88704),;                                                                                       G_3 = - 23476167/1024  a^13 (a^4-32)^3  
                  (2304313595 a^8 - 1702233920 a^6 - 11829269248 a^4;                                                                                                                                                         -  39211065344 a^2+8642101248),;                                                                                       G_4 =                                                                  - 75246080/ a^10 (a^4-32),;                                                                                       G_5 =                               9405760/3  a^12(a^4-32)^2 
                  (75767 a^4-46944 a^2-96768),;                                                                                       G_6 =         - 180880/3  a^14 (a^4-32)^3 
                  (11681524055 a^8-8555309984 a^6 -56944147200 a^4;                                                                                     [0.5ex]                                                           -  204210659328 a^2+30640177152),;                                                                                       G_7 =                                                   2006968901247765/2883584  
            a^11 (a^4-32),;                                                                                       G_8 =                     - 154382223172905/2883584  
            a^13(a^4-32)^2 (48667 a^4-30160 a^2-52416),;                                                                                       G_9 =           66163809931245/46137344 a^15 
            (a^4-32)^3 
                     (6314158847 a^8-4591849024 a^6;                                                                                                                                      - 29599122432 a^4 - 112639700992 a^2+11915624448), ]
          
        Note that F_i  0,  i=1,2,3, and 
        G_i  0,   i=1,2, …, 9, since 
        a^4 - 32 > 0 for a < -2^-5/4. 
        
        
        
        
         §.§ ε- and ε^2-order analysis
         
        
        
        The ε-order focus values 
        V_11, V_21, ⋯, V_111 are obtained 
        by using the algorithm and Maple program developed in <cit.>.
        Their expressions are lengthy, and here we only present the first one 
        for brevity, 
        
            [                                                         V_11              = - 1/64 a (a^4-32){ 6912 b_101 -5760 a   b_201;                                                              +16 (a^4-36 a^2+40)   b_111 
            +48 a (a^4 +36 a^2 +160)  b_021;                                                                             +3456 a^2  b_301 -24 a (a^4-12 a^2+40)  b_211;                                                                                      -16 (3 a^6+68 a^4+300 a^2+20)  b_121;                                                                                   -24 a (3 a^6+65 a^4+300 a^2+224)  b_031;                                                                             +27 (a^2+2) (7 a^6+82 a^4+320 a^2+128)  a_031;                                                                                +8 a (24 a^6+223 a^4+1140 a^2+1180)  a_121;                                                                                   +8 (21 a^6-73 a^4+480 a^2+320)  a_211}. ]
         
        It is noted that all V_i1's are linear polynomials in
        a_ij1 and b_ij1. It can be shown that 
        
            _1^8 = F_1  0,   _1^9 = F_2  0,   _1^10 = 0.
        
        In fact, with the solution of S_1^8 solved from W_1^8=0, we obtain 
        
            [ V_91 = G_1   a_031,    
            V_101 = G_2   a_031,    
            V_111 = G_3   a_031, ]
         
        where G_i's are given in (<ref>). Noticing 
        G_1  0 for a < - 2^5/4, we have 
        V_91 0 if a_031 0. 
        Moreover, _1^8  0 and (<ref>)
        indicate that perturbing W_1^8 and V_01 
        around the solutions S_1^8 and b_011 (see (<ref>)) does yield 
        9 small limit cycles around the equilibrium E_0. 
        
        It is seen from (<ref>) that 
        V_91= V_101= V_111= 0 for a_031= 0. 
        For convenience, define the critical condition S_1c^8, 
        satisfying (<ref>), W_1^8 = 0 and a_031= 0, as 
        
            S_1c^8 : {[                                                                     b_011 =                                                   C_1   a_121 -9/8a^3a_211,;                                                                     a_031 = b_211 = 0,   b_121 = 7/2  a_211,  b_021 = - 6  a_121,   b_031 = 8/3  a_121,;                                                                     b_111 =               9 a  a_121 + 9/2  a_211,  
            b_101 = C_2   a_121 + C_3   a_211,;                                                                     b_201 =           C_4   a_121 + C_5  a_211,     
            b_301 = C_6   a_121 + C_7   a_211, ].
        
        where C_i's are given in Appendix A.  
        
        We have the following result.  
        
         
        
        The equilibrium E_0 of system (<ref>)
        is a center up to ε-order, 
        i.e. all ε-order focus values vanish  
        if and only if the condition S_1c^8 holds.
        Furthermore, there exist at most 9 small limit cycles
        around E_0 for all parameters a_ij1 and 
        b_ij1, and 9 small limit cycles can be obtained for some parameter
        values near S_1c^8. 
         
        The existence of 9 small limit cycles has been shown under 
        the solution S_1^8 with 
        a_031 0 and det_1^8  0. 
        It is obvious that the critical condition S_1c^8 
        is  necessary for all ε-order focus values 
        to vanish since V_i1 = 0,  0 ≤ i ≤ 11 under this 
        condition. To prove  sufficiency,  
        under the critical condition S_1c^8, 
        we use (<ref>) to obtain 
        the following ε-order approximation of the first integral, 
        
            H_1(x,y,ε)=f_1+ε f_11/f_2+ε f_21,
        
        where f_1 and f_2 are given in (<ref>), 
        f_11=a_121r_1+a_211r_2 and
        f_21=a_121r_3+a_211r_4 with
        
            [                                                                                 r_1 =                                   -1/48[a^2(3a^2+4)(5+2y+2x^2+x^4)
            +220 -192ax + 280y;                                                                                                                                   +120x^2 -64y^2 +128ax^3 +64x^2y +76x^4 ],;                                                                                 r_2 =                              -1/8a(5a^2-4) +5x -1/8a(a^2-4)(2y+2x^2+x^4)
             +2xy -2x^3,;                                                                                 r_3 =                        1/192[a^2(3a^2 + 4)(4a - 15x 
            + 10xy + 10x^3) +304a+16a^3-180x;                                                                                                            -  40 ( 16ay - 8ax^2 + 5xy - 23x^3
            - 8xy^2 + 16ax^4 + 8x^3y ) ],;                                                                                  r_4= a^2/8(a^2 - 1)
            - a(15/32a^2 +5/8)x 
            +5/2x^2(3/2+ y - x^2) 
            +5/16a(a^2 - 4)x(y + x^2). ]
         
        This implies that setting 
        the first 10 focus values V_i1=0, i=0,⋯,9
        yields V_i1=0 for all i ≥ 10.
        Moreover, due to that all V_i1 are linear in  
        all parameters a_ij1 and b_ij1,   
        by Theorem <ref> at most 9 
        small limit cycles can be obtained for this case. 
        The proof is complete. 
         
        
        
        
        Now suppose the condition S_1c^8 holds and so all 
        ε-order 
        focus values vanish, we then need to use the ε^2-order focus 
        values to study bifurcation of limit cycles. 
        With an almost exact same procedure as that used in 
        the ε-order analysis, 
        we can find a solution S_2^8 such that W_2^8=0, and then 
        
            [ V_92 = G_1   a_032,    
            V_102 = G_2   a_032,    
            V_112 = G_3   a_032,   _2^8 = F_1  0, ]
         
        where F_1 and G_i's are given in (<ref>). 
        Note that the above equations are exactly the same as those given in 
        (<ref>) and (<ref>), if we replace k = 1 by k = 2 
        in (<ref>) and (<ref>). 
        This clearly shows that there can exist 9 limit cycles around the 
        equilibrium E_0 when all ε-order focus values vanish. 
        It is also noted that all V_i2 are linear polynomials 
        in a_ij2 and b_ij2. 
        
        Similarly, we see that setting a_032= 0 in (<ref>) yields 
        V_92= V_102=  V_112= 0, 
        implying that the solution S_2^8 with a_032= 0
        and b_012 given in (<ref>) defines 
        a necessary condition for all ε^2-order focus values to 
        vanish. This critical condition is given below: 
        
            S_2c^8: {[                                                                                       b_012 =               9/64 a ^4 a_211^2
            - 9/8 a^3 a_212  + C_1 a_122 + C_8 a_121a_211
            + C_9 a_121^2 ,;                                                                                       a_032 =            0,   
            b_032 = 8/3  a_122 +5  a_121^2 ,   
            b_212 = a/2  a_121 (5 a  a_121-9 a_211),;                                                                                       b_122 =                                             7/2  a_212 
            - 1/4  a_121 (31 a  a_121-45 a_211) ,;                                                                                       b_102 =                C_2  a_122 + C_3   a_212 
            + C_10  a_121^2 + C_11  a_211^2 + C_12  a_121 a_211,;                                                                                       b_202 =               C_4  a_122 + C_5   a_212 
            +  C_13  a_121^2 + C_14  a_211^2
            + C_15  a_121 a_211,;                                                                                       b_112 = 9 a  a_122+ 9/2  a_212 
                   + 9a^3/32 a_211^2 
                   + C_16  a_121^2 + C_17  a_121 a_211,;                                                                                       b_022 =                                               - 6 a_122 + C_18  a_121^2 
            + C_19  a_121 a_211,;                                                                                       b_302 =               C_6   a_122 + C_7   a_212 + C_20  a_121^2 
            + C_21  a_211^2 + C_22  a_211 a_121, ].
         
        where C_i's are given in Appendix A.  
        
        We have the following theorem. 
        
        
        
        Assume S_1c^8 holds.
        The equilibrium E_0 of system (<ref>)
        is a center up to ε^2-order, 
        if and only if S_2c^8 holds.
        Furthermore, there exist at most 9 small limit cycles
        around E_0 for all parameters a_ij2 
        and b_ij2, and 9 small limit cycles exist for some parameter values 
        near S_2c^8.  
        
        Similarly, we only need to prove  sufficiency.
        With S_1c^8 and S_2c^8 holding,
        we can use (<ref>) to find 
        the following ε^2-order approximation of the first integral, 
        
            H_2(x,y,ε)=f_1+ε f_11+ε^2 f_12/f_2+ε f_21+ε^2f_22,
        
        where f_11 and f_21 are given in H_1(x,y,ε) (see 
        Eq. <ref>), and 
        
            [ f_21=a_122r_1+a_212r_2+a_121^2s_1+a_211^2s_2+a_121a_211s_3,; f_22=a_122r_3+a_212r_4+a_121^2s_4+a_211^2s_5+a_121a_211s_6, ]
         
        in which r_i,   i=1,2,3,4 are given in (<ref>), 
        and s_i,  i=1,2, …, 8, are listed in Appendix B. 
        
        The existence of 9 small limit cycles is easily seen from 
        V_92 0 and _2^8  0 when a_032 0 
        under the critical condition S_2c^8. 
        On the other hand, the above results show that setting 
        V_i2, 0≤ i≤ 9 results in 
        V_i2 =  0 for all i ≥ 10.
        Further, all V_i2's are linear in a_ij2 and b_ij2, 
        and S_2c^8 is the unique solution of V_i2=0, 0≤ i≤ 9.
        Then by Theorem <ref>, 
        at most 9 small limit cycles can be obtained around E_0
        for all parameters a_ij2 and b_ij2.
        
        
         §.§ ε^3-order analysis
         
        
        In this section, we assume 
        the critical condition { S_1c^8, S_2c^8 }, 
        which stands for that both 
        the critical conditions S_1c^8 and S_2c^8 hold, 
        under which all ε- and ε^2-order focus values vanish. 
        Thus, we use ε^3-order focus values V_i3 to 
        study bifurcation of limit cycles around the equilibrium 
        E_0. With a similar procedure, but for this order, 
        we solve 9 equations  W_3^9=0 to obtain the solution S_3^9 
        for which  
        
            [ V_103 = G_4   a_121^3,    
            V_113 = G_5   a_121^3,    
            V_123 = G_6   a_121^3,   _3^9 = F_2  0, ]
         
        where F_2 and G_i's are given in (<ref>). 
        Note that for this order, there is one more independent coefficient
        a_033 in S_3^9 for solving W_3^9 =0, 
        compared to the solutions S_1^8 and S_2^8 which have only 
        8 independent coefficients to be used for solving the first 
        8 focus value equations. 
        The equations in (<ref>) show that when all 
        ε- and ε^2-order focus values vanish, 
        the ε^3-order focus values can have solutions such that 
        V_i3=0,   i = 0, 1, ⋯, 9 but V_103 0, as well 
        as _3^9 0, implying that 10 small limit cycles 
        can bifurcate from the equilibrium E_0. 
        
        Setting a_121 = 0 in (<ref>), 
        we have V_103 = V_113 = V_123 =0, implying that 
        under the solution S_3^9 with a_121=0 and 
        b_013 given in (<ref>), the equilibrium E_0 
        might be a center up to ε^3 order. 
        This critical condition is given by 
        
            S_3c^9 : {[                                                                                   b_013 =      9/32 a^4 a_211 a_212 
            - 9/8 a^3 a_213 + C_1 a_123 + C_8 a_122 a_211 
            + C_23 a_211^3,;                                                                                   a_121 =                a_033 = 0,    
            b_033 = 8/3 a_123,    
            b_023 = - 6 a_123 + C_19 a_122 a_211;                                                                                   b_213 = - 9a/16 a_211 (8 a_122+ a_211^2), 
                b_123 = 7/2 a_213+45/32 a_211 (8 a_122+ a_211^2 ),;                                                                                   b_113 =    9 a  a_123 + 9/2 a_213 
                    + a_211[ 9/16 a^3 a_212 + C_17 a_122 
            + C_24 a_211^2 ];                                                                                   b_103 =             C_2 a_123 + C_3 a_213  
            + a_211[ 2 C_11  a_212 +C_12 a_122 
            + C_25 a_211^2 ],;                                                                                   b_203 =                 C_4 a_123 + C_5 a_213
            + a_211[ 2 C_14 a_212 +C_15 a_122 
            + C_26 a_211^2 ];                                                                                   b_303 =              C_6 a_123 + C_7 a_213 
            + a_211[ 2 C_21 a_212 + C_22 a_122 
            + C_27 a_211^2 ], ].
        
        under which the critical conditions S_1c^8 and S_2c^8 
        are simplified. Here, C_i's are given in Appendix A.  
        
        We have the following theorem. 
        
        
        Let { S_1c^8, S_2c^8} hold.
        The equilibrium E_0 of system (<ref>)
        is a center up to ε^3-order, 
        if and only if the condition S_3c^9 holds.
        Furthermore, there exist 10 small limit cycles
        around E_0 for some parameter values of a_ij3 and b_ij3 
        near the critical value defined by S_3c^9 when V_103≠ 0. 
        
        Similarly again, we only need to prove  sufficiency.
        Under the condition { S_1c^8, S_2c^8, S_3c^9}, 
        we obtain the following ε^3-order
        approximation of first integral, 
        
            H_3(x,y,ε)=f_1+ε a_211r_1
            +ε^2 (a_122r_1+a_212r_2+a_211^2s_2)
            +ε^3 f_31/f_2+ε a_211r_4
            +ε^2 (a_122r_3+a_212r_4+a_211^2s_5)
            +ε^3 f_32,
        
        where 
        
            f_31=a_123r_1+a_213r_2
            +a_211(a_122t_1+a_212t_2+a_211^2t_3),
               f_32=a_123r_3+a_213r_4
            +a_211(a_122t_4+a_212t_5+a_211^2t_6),
        
        in which r_i,  i=1,2,3,4 are given in (<ref>), 
        and s_2, s_5 and t_i,  i=1,2, …, 6 are listed in Appendix B. 
        This implies that setting 
        V_i3=0, 0≤ i≤ 10 yields V_i3 = 0 for all i ≥ 11. 
        Then, there exist at most 10 small limit cycles 
        for this case. On the other hand, 10 small limit cycles exist 
        since when a_121 0, V_101 0 and 
        _3^9  0. 
        
        
        
        
         §.§ ε^4–ε^6-order analysis
         
        
        The analyses for ε^4-,  ε^5- 
        and ε^6-order are  similar to that 
        of ε^1-, ε^2- and ε^3-order, 
        respectively. 
        
        Let { S_1c^8, S_2c^8, S_3c^9} hold.
        Following the same procedure used in the previous sections, we 
        can solve the equations W_4^8=0 to obtain a solution S_4^8 such that 
        
            V_94 = G_1   a_034,    
            V_104 = G_2   a_034,    
            V_114 = G_3   a_034,   _4^8 = F_1  0,
         
        which has the exactly same form of the equations as those given in 
        (<ref>) and (<ref>), implying that perturbing the 
        ε^4-order focus values from the solution S_4^8
        and b_014 (see (<ref>)) 
        can yield 9 limit cycles around the equilibrium E_0. 
        Similarly, the solution S_4^8 and b_014 with a_034 = 0  
        yields a critical condition S_4c^8, under which 
        the equilibrium E_0 is a center up to ε^4 order.  
        
        Then let { S_1c^8, S_2c^8, S_3c^9,S_4c^8} hold.
        In the same line, we can solve the equations 
        W_5^8=0 to obtain a solution S_5^8 such that 
        
            [                                                                       V_95 = G_1  A_035,   
            V_105 = G_2   A_035,   
            V_115 = G_3   A_035,   _5^8 = F_1  0, ]
        
        where 
        
            A_035 = a_035
            + 1/48  a_122  a_211  (140   a_122+35   a_211^2).
        
        This shows that perturbing the ε^5-order 
        focus values near the solution S_5^8 
        and b_015 given in (<ref>) 
        can also yield 9 limit cycles around the equilibrium E_0. 
        It is easy to see that the solution of A_035 = 0, 
        
            a_035 = - 35/48  a_122  a_211  (4   a_122+ a_211^2),
         
        yields V_95= V_105 = V_115 = 0.
        Now, we combine the solution S_5^8, b_015 and a_035 
        to obtain the critical condition S_5c^8, under which
        the equilibrium E_0 becomes 
        a center up to ε^5 order.  
        
        The lengthy critical conditions S_4c^8 and S_5c^8 are 
        omitted here for brevity. 
        Summarizing the above results leads to the following theorem. 
        
        
        
        System (<ref>) can have maximal 9 limit cycles 
        around the equilibrium E_0 under 
        the condition { S_1c^8, S_2c^8, S_3c^9} 
        by perturbing the ε^4-order focus values
        around the critical value S_4c^8; 
        and under the critical condition 
        { S_1c^8, S_2c^8, S_3c^9, S_4c^8} 
        by perturbing the ε^5-order 
        focus values near the critical point S_5c^8. 
        The equilibrium E_0 becomes a center up to ε^4 order  
        under the condition { S_1c^8, S_2c^8, S_3c^9, S_4c^8}, 
        and a center up to ε^5 order under the condition 
        { S_1c^8, S_2c^8, S_3c^9, S_4c^8, S_5c^8}.  
         
        The proof for the center conditions in Theorem <ref> 
        is similar to that in proving Theorems <ref>, <ref>
        and <ref> by finding 
        the ε^4-order and ε^5-order approximations 
        of the first integrals. This is the major and tedious part. 
        For higher-order analysis, the proofs are similar. 
        We omit the detailed proofs in the following analysis for brevity. 
         
        
        Next, suppose the condition 
        { S_1c^8, S_2c^8, S_3c^9, S_4c^8, S_5c^8} 
        is satisfied, then all 
        ε^k,  k=1,2, …, 5, order focus values vanish. 
        Following a similar analysis as that for ε^3 order, we 
        solve the equations W_6^9=0 
        to obtain a solution S_6^9 such that 
        
            [                                                    V_106 =                  G_4   
            a_122^2 ( a_122 + 9/8  a_211^2 ) ,;                                                    V_116 =                  G_5   
            a_122^2 ( a_122 + 9/8  a_211^2 ) ,;                                                    V_126 = G_6   
            a_122^2 ( a_122 + 9/8  a_211^2 ) ,   _6^9 = F_2  0, ]
         
        which indeed shows the existence of 10 limit cycles around the 
        equilibrium E_0, generated from perturbing the ε^6-order
        focus values near the solution S_6^9 
        under the condition 
        { S_1c^8, S_2c^8, S_3c^9, S_4c^8, S_5c^8}. 
        Moreover, when a_122= -9/8  a_211^2 or a_122=0, 
        we have V_106= V_116= V_126= 0, 
        indicating that the solution S_6^9 with either a_122= 
        -9/8  a_211^2 or a_122= 0, plus 
        b_016 given by (<ref>)|_k=6, yields a 
        critical condition S_6c^9a (corresponding to 
        the former) or S_6c^9b (corresponding to the latter) 
        under which all ε^6-order focus values vanish. 
        Thus, under the critical condition 
        { S_1c^8, S_2c^8, S_3c^9, S_4c^8, S_5c^8, S_6c^9}
        (S_6c^9 
        equals either S_6c^9a or S_6c^9b), 
        the equilibrium E_0 becomes a center 
        up to ε^6 order. 
        
        We have the following theorem for this order. 
        
        
        System (<ref>) can have maximal 10 limit cycles 
        bifurcating from the equilibrium E_0 under 
        the condition 
        { S_1c^8, S_2c^8, S_3c^9, S_4c^8, S_5c^8} 
        by perturbing the ε^6-order 
        focus values near the critical point S_6c^9a 
        or S_6c^9b. 
        Further, the equilibrium E_0 becomes a center 
        up to ε^6 order under the condition 
        { S_1c^8, S_2c^8, S_3c^9, S_4c^8, S_5c^8, S_6c^9 }, 
        for which all ε^k-order  (k=1,2, …, 6 ) 
        focus values vanish. 
         
        
        
        Suppose the condition 
        { S_1c^8, S_2c^8, S_3c^9, S_4c^8, S_5c^8, S_6c^9 } holds. 
        Then, all the ε^k-order 
        (k=1,2, …, 6) focus values vanish. We have two cases 
        for higher-order analysis, defined as 
        
            [ Case  (A)   { S_1c^8, S_2c^8, S_3c^9, S_4c^8, S_5c^8, S_6c^9a},;  Case  (B)   {S_1c^8, S_2c^8, S_3c^9, S_4c^8, S_5c^8, S_6c^9b}. ]
        
        
         §.§ Higher-order analysis for Case (A)
         
        
        First we consider Case (A), under which 
        we will show that 11 limit cycles can bifurcate from the equilibrium 
        E_0 based on the ε^7-order focus values. 
        
        
        
          §.§.§ ε^7-order analysis
         
        
        Under the condition (A) defined in (<ref>) 
        with a_122 = - 9/8 a_211^2, we obtain 
        
            _7^10 = F_3  a_211^4 ,   _7^11 = F_4  a_211^10,
          
        which shows that _7^10 0 and 
        _7^11 0 when a_211 0, implying that 
        we may have solutions such that the first ten focus values vanish
        but V_117 0 and so 11 small limit cycles may be obtained.   
        Indeed, we can solve the first ten focus values 
        equations: W_k^10 = 0 to obtain a solution S_7^10 such that  
        
            V_117 = G_7   a_211^7,    
            V_127 = G_8   a_211^7,    
            V_137 = G_9   a_211^7,
         
        which clearly shows that V_117 0 if a_211 0. 
        In addition, due to _7^10 0 when a_211 0,
        implying that 11 small limit cycles exist. 
        
        Letting a_211= 0, we have V_117= V_127= V_137= 0, leading to a critical condition S_7c^10,
        defined by 
        
        
            S_7c^10: {[                                                                 b_017 =            C_1 a_127-9/8 a^3 a_217+ C_8 C_28+9/32 a^4 C_29+ C_23 C_30 ,;                                                                 a_211 =                                a_123 = a_037 = 0,  b_037 
            = 8/3  a_127,;                                                                 b_027 =            -6 a_127  +  C_19 C_28,    
            b_217 = -9a/16 (8 C_28 + C_30 ),;                                                                 b_127 =                                   7/2 a_217 
            + 45/32 ( 8 C_28 + C_30 ),;                                                                 b_107 =           C_2 a_127 + C_3 a_217 
            + 2 C_11 C_29 + C_12 C_28 + C_25 C_30,;                                                                 b_207 = C_4 a_127 + C_5 a_217 
            + 2 C_14 C_29 + C_15 C_28 
                    + C_26 C_30 ,;                                                                 b_117 =       9 a  a_127+ 9/2  a_217 
            + 9 a^3/16 C_29 + C_17 C_28 + C_24 C_30 ,;                                                                 b_307 = C_6 a_127 + C_7 a_217  
            + 2 C_21 C_29 + C_22 C_28 
                    + C_27 C_30, ].
         
        where C_i's are given in Appendix A. 
        
        We have the following result. 
        
        
        Let { S_1c^8, 
        S_2c^8, S_3c^9, S_4c^8, S_5c^8, S_6c^9a} hold.
        The equilibrium E_0 of (<ref>) becomes a center 
        up to ε^7 order under 
        S_7c^10 for which all ε^7-order 
        focus values vanish. Furthermore, there exist 
        11 small limit cycles around E_0 
        for parameter values of a_ij7 and b_ij7 near the critical point 
        S_7c^10. 
        
        
          §.§.§ Higher-order analysis
         
        
        
        For higher-order analysis (k ≥ 8), 
        we first briefly list the results for a few orders to see the 
        patterns and then summarize the results in a table for higher orders.  
        
        The analysis on ε^k (k=8,9,10,11) orders 
        show the same pattern, giving 
        9 limit cycles for each order, as follows: 
        
            [                                                                                                                               [      Order k:; (k=8,9,10,11) ]                                                                                                                                                 {S_k^8,W_k^8}, {[                                V_9k = G_1   A_03k,    
            V_10k = G_2   A_03k,;                     [0.7ex]  
            V_11k =           G_3   A_03k,   _k^8 = F_1 , ]. ]
         
        where {S_k^m,W_k^m} denotes the solution S_k^m solved from 
        W_k^m=0, and 
        
            [                                                                               A_038 =                                                            a_038,      A_039 = a_039,;                                                                              A_0310 =                                  a_0310 + 35/48   a_124  a_212  (4   a_124+ a_212^2),;                                                                              A_0311 = a_0311 
            + 35/48[ a_125 a_212 (8 a_124+a_212^2) 
            + a_124 a_213 (4 a_124 +3 a_212^2) ]. ]
         
        This clearly shows that for each order of k=8,9,10,11, one can 
        solve A_03k=0 to get a unique solution for a_03k 
        under which (together with the solutions S_k^m 
        and b_01k obtained in the previous 
        orders and the current order)
        the equilibrium E_0 becomes a center up to that order.  
        
        When the equilibrium E_0 is a center up to 11th order, 
        as given in (<ref>)
        we obtain the following result for order 12:  
        
            [                                                                                                                                                                                                                                                                                                                 Order 12:                                                                                                                                                                                                                                                                                                          {S_12^9,W_12^9}, {[                                           V_1012 =                 G_4  a_124^2 (a_124+9/8 a_212^2 ),;                                 [1.0ex]  
            V_1112 =                G_5  a_124^2 (a_124+9/8 a_212^2 ) ,;                                  [1.0ex] 
            V_1212 = G_6  a_124^2 (a_124+9/8 a_212^2 ),   _12^9 = F_2 , ]. ]
         
        which has the exactly the same pattern as order 6, shown in 
        (<ref>), indicating that 10 limit cycles can be obtained from 
        this order, and there are two solutions from the 
        equations V_1012 = V_1112= V_1212 = 0: 
        a_124= -9/8 a_212^2 and a_124= 0, which are 
        again similar to that as 
        in order 6. When a_124 = 0, it will be shown 
        in Section 4.7 that 
        it yields the same pattern as that for Case (B) in higher orders.    
        So in this section, we choose a_124= -9/8 a_212^2, 
        like we chose a_122= -9/8 a_212^2 in order 
        6 to obtain the center condition. 
        
        Let a_124= -9/8 a_212^2, under which (together with 
        the solutions obtained from previous orders and this order) the equilibrium 
        E_0 becomes a center up to ε^12 order. 
        Then, we have the result for ε^13 order: 
        
            [                                                                                                                                                                                                                                                                                                                                                                                   Order 13:                                                                                                                                                                                                                                                                                                                                                                            {S_13^9,W_13^9}, {[                                                      V_1013 =                 81/64 G_4  a_212^4 
            (a_125+ 9/4 a_212 a_213),;                                                      V_1113 =                 81/64 G_5  a_212^4 
            (a_125+ 9/4 a_212 a_213),;                                                      V_1213 = 81/64 G_6  a_212^4 
            (a_125+ 9/4 a_212 a_213),   _13^9 = F_2 , ]. ]
          
        which shows that perturbing ε^13-order focus values 
        can also yield 10 small limit cycles around the equilibrium E_0. 
        It can be seen from (<ref>) that either a_212= 0 or 
        a_125= -9/4 a_212 a_213 leads to 
        the equilibrium E_0 being a center. However, it can be shown 
        that setting a_212= 0 at this order will not yield 11 
        small limit cycles at the next order though it will resume the same pattern 
        at higher orders. 
        
        So let a_125= -9/4 a_212 a_213. Then, we obtain 
        the following result for ε^14 order: 
        
            [                                                                                                                                                                        Order   14:                                                                                                                                                                 {S_14^10,W_14^10}, {[                                   V_1114 = G_7   a_212^7,   
            V_1214 = G_8   a_212^7 ,;                          [0.7ex] 
            V_1314 =          G_9   a_212^7,   _14^10 = F_3  0, ]. ]
         
        which shows that perturbing ε^14-order focus values
        can yield 11 limit cycles around the equilibrium E_0, 
        and setting a_212 = 0 leads to E_0 being a center 
        up to ε^14 order. 
        It has been noted that choosing a_212= 0 at order 13 or 14 makes
        differences. More precisely, as shown in Table <ref>, if taking 
        a_125= -9/4 a_212 a_213 at 
        order 13, we have small limit cycles 11,   9,   9,   9,  9 
        for the orders 14–18; while if taking a_212= 0 at order 13, 
        then the limit cycles obtained for the orders 14–18  
        are 9,  10,   9,   9,   10, and then the two different choices 
        merge into the same pattern from order 19. Note that the choice 
        a_212= 0 at order 13 does not yield 11 small limit cycles 
        at order 14, but 
        gives two more 10 small limit cycles at orders 15 and 18. 
        However, it returns to the general pattern at order 19.
        So we treat a_212= 0 as a special case of the 
        case a_125= -9/4 a_212 a_213. 
        
        Summarizing the above results we have the following pattern: 11 
        limit cycles are obtained from ε^7 order, then 9 limit cycles 
        from four consecutive ε^k orders (k=8,9,10,11), and then 
        10 limit cycles from two consecutive ε^k orders (k=12,13), 
        and finally return to 11 limit cycles at ε^14 order.    
        This pattern, starting from order 8, four 
        9 limit cycles, followed by two 10 limit cycles, and then 
        11 limit cycle, has been verified up to ε^35 order. 
        We call this as 9^4-10^2-11^1 generic pattern, and the corresponding 
        solution (or center condition) is called  generic solution 
        (or generic center condition). 
        By  generic we mean that 
        one should always choose a non-zero solution (if it exists) 
        when one solves the center conditions at each order.  
        Other types of solutions are called  non-generic. 
        For example, as discussed above, if choosing the non-generic solution 
        a_212= 0 at order 13, then 11 limit cycles will be missed  
        at order 14 but the solution procedure will return to the 
        generic 9^4-10^2-11^1 pattern at order 19.  
        However, it should be noted that a non-generic solution in Case (A) 
        does not always lead to the generic 9^4-10^2-11^1 pattern. 
        For instance, choosing the non-generic solution 
        a_124= 0 at order 12 will generate solutions 
        in the form of generic patter of Case (B) at a higher order, 
        as shown in the next section.  
        
        
        
        It has been observed from the above analysis, the values of the 
        parameter a in the Hamiltonian function does not affect the number of 
        limit cycles. In other words, a can not be used to increase the 
        number of bifurcating limit cycles. Thus, to simply the computations 
        in higher order analysis, we set a=- 3 in higher-order 
        (k ≥ 15) computations, which greatly simplify the computations. 
         
        
        We summarize the results of Case (A) in Table <ref>, 
        where k is the order of ε^k focus values, 
        (S_k^m,W_k^m) represents the solution S_k^m solved from 
        W_k^m = 0, and LC denotes limit cycles 
        around the equilibrium E_0 obtained by perturbing the 
        ε^k-order focus values. 
        The “Condition for Center” in each row only lists the condition for 
        the current row, which assumes that the conditions in the previous 
        rows hold. For example, when k= 4, S_4c^8 only gives the 
        center condition for k = 4, which should be combined with the 
        conditions given in the previous rows: S_1c^8, 
        S_2c^8 and S_3c^9 to get a complete 
        center condition { S_1c^8,S_2c^8,S_3c^9,S_4c^8}. 
        Note that the critical condition S_kc^8 contains the solutions 
        S_k^8, the  b_01k given in (<ref>) 
        and a particular coefficient. For example, S_2c^8={S_2^8, b_012, a_032}, S_3c^9={S_3^9,
         b_013 , a_121}, and 
        S_7c^10={S_7^10,
        b_017 , a_211}, etc. 
        The solutions of these key coefficients are given below. 
        
            [                                                                          9  LC:                                                                     k=1,2,4,8,9                                                                         a_03k=0;                                                                                                                                                             k=5                                     a_035=-35/48 a_122a_211 (4 a_122 + a_211^2);                                                                         [0.5ex]                                                                            k=10                                   a_0310=-35/48 a_124a_212 
            (4 a_124 + a_212^2);                                                                         [0.5ex]                                                                            k=11                                  a_0311 = -35/48[ a_212 a_125 (8 a_124+a_212^2);                                                                                                                                                                                                              +a_213 a_124 (4 a_124+3 a_212^2) ];                                                                                                                                                            k=15                                                        a_0315= -735/256 a_213^5;                                                                         [0.5ex]                                                                            k=16                           a_0316= 35/768 a_213^3  
            (128 a_127 - 27 a_214 a_213);                                                                                                                                                            k=17                       a_0317 = -35/768 a_213[ 
            32 a_127 (2 a_127-3 a_214 a_213);                                                                                                                                                                                            -a_213^2 (128 a_128+ 54 a_214^2
            - 27 a_215  a_213) ];                                                                         [0.5ex]                                                                            k=18 a_0318 = 35/6 a_213^3 a_129 
                    -35/24 a_213 a_128 (4 a_127-3 a_214 a_213);                                                                                                                                                                                - 35/48 a_127[ a_214 (4 a_127+ 3 a_214 a_213)
            - 6 a_213^2 a_215];                                                                                                                                                                                           + 3/4096 a_213^2 [ 1120 a_214 (a_214^2+6 a_215 a_213);                                                                                                                                                                                                           +3 a_213^2 (2269 a_213^2-560 a_216) ];                                                                         [0.5ex]                                                                            k=22              a_0322= 35/768 
            a_214^3 (128 a_1210- 486 a_215^2 - 27 a_216 a_214);                                                                                                                                                            k=23                   a_0323 = 35/768 a_214^2 [ 
            a_214 (128 a_1211 -27 a_217 a_214);                                                                         [0.5ex]                                                                                                             +a_215 (384 a_1210 -108 a_214 a_216 -198 a_215^2) ];                                                                         [0.5ex]                                                                         k=24,25                                                                       a_03k = ⋯;                                                                                                                                                         k=29–32                                                                       a_03k = ⋯;                                                               [1.0ex] 
             10  LC:                                                                             k=3                                                                       a_121 = 0;                                                                                                                                                             k=6                                                             a_122 =-9/8 a_211^2;                                                                                                                                                            k=12                                                             a_124 =-9/8 a_212^2;                                                                                                                                                            k=13                                                         a_125 =-9/4 a_213 a_212;                                                                                                                                                            k=19                                                         a_127 =-9/4  a_213a_214;                                                                                                                                                            k=20                                             a_128 =-9/8 (a_214^2+2 a_215 a_213);                                                                                                                                                            k=26                                            a_1210 =-9/8 (a_215^2+2 a_216 a_214);                                                                                                                                                            k=27                                          a_1211 =-9/4 (a_217 a_214+a_216 a_215);                                                                                                                                                            k=33                                          a_1213 =-9/4 (a_218 a_215+a_217 a_216);                                                                                                                                                            k=34                          a_1214 =-9/8 (a_217^2 +2 a_219 a_215 
            + 2 a_218 a_216);                                                               [1.0ex] 
             11  LC:                                                                            k=7m                                                                                ;                                                                                                                                                          m= 1–5                                                                       a_21m = 0 ]
        
        where `⋯' 
        represents the omitted lengthy expressions for brevity. 
        In addition, in Table <ref>, the blue and red colors denote 
        the solutions and center conditions corresponding to the 
        10 and 11 small limit cycle, respectively. 
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
         
        
        
        
         §.§ Higher-order analysis for Case (B)
         
        
        We now turn to Case (B) for which we choose a_122= 0 at 
        ε^6 order. 
        Thus, the results starting from ε^6 order 
        are different from those given in Table <ref>. 
        Now under the condition a_122= 0, together with the 
        conditions obtained in previous 
        orders, the equilibrium E_0 becomes a center up to 
         ε^6 order. Then for ε^7-order  
        focus values we solve W_7^8 =0 to obtain S_7^8 and then 
        
            V_97 = G_1   A_037,    
            V_107 = G_2   A_037,    
            V_117 = G_3   A_037,    _7^8 = F_1,
         
        where 
        
            A_037 = 
            a_037+35/768 a_211[ a (a^2 - 4) a_123 a_211^3
            + 16 a_124 a_211^2 + 16 a_123 (4 a_123+ 
            3 a_212 a_211) ],
         
        which shows that for Case (B) only 9 small limit cycles can be obtained from 
        ε^7-order. 
        Then, solving A_037= 0 gives a unique solution for a_037, 
        under which, together with the conditions obtained in the previous 
        orders as well as S_7^8 and b_017, 
        the equilibrium E_0 becomes a center up to 
        ε^7 order.  
        
        Next, the ε^8-order analysis shows that 10 limit cycles 
        can be obtained by solving W_8^9 =0 to have the solution S_8^9, 
        under which higher-order focus values become 
        
            V_108 = 9/8 G_4 a_211^2 a_123^2,    
            V_118 = 9/8 G_5 a_211^2 a_123^2,    
            V_128 = 9/8 G_6 a_211^2 a_123^2,   _8^9 = F_2.
         
        This clearly indicates that either a_211= 0 or a_123= 0, 
        together with b_018,  
        leads to the equilibrium E_0 being a center up to 
        ε^8 order. 
        If taking a_211= 0, then we again obtain 10 small limit cycles 
        from ε^9 order by solving W_9^9 =0 to obtain the solution 
        S_9^9 and 
        
            V_109 = G_4  a_123^3,      
            V_119 = G_5  a_123^3,      
            V_129 = G_6  a_123^3,   _9^9 = F_2.
         
        Thus, for the equilibrium E_0 being a center up to 
        ε^9 order, a_123 must be taken zero (with 
        b_019), yielding the same result as that generated from Case (A) at 
        order 9 (and so the result at order 8 also becomes same). 
        In other words, choosing the non-generic solution a_211= 0 
        at order 8 makes the higher-order solutions (k ≥ 9) follow the 
        generic pattern of Case (A). 
        
        Now we consider the choice a_123= 0 at ε^8 order. 
        It can be shown that under this condition only 9 limit cycles 
        exist for ε^9 order. Then for the ε^10 order, 
        we solve W_10^9=0 to obtain the solution S_10^9 and then get 
        
            [                                                                             V_1010 =              9/8  G_4  
            a_211^2 (a_124^2- 5/16 a_211^4 a_124+ 429/40960 
             a_211^8 ),;                                                                             V_1110 =                9/8  G_5 
            a_211^2 (a_124^2- 5/16 a_211^4 a_124+ 429/40960
             a_211^8 ),;                                                                             V_1210 = 9/8  G_6 
            a_211^2 (a_124^2- 5/16 a_211^4 a_124+ 429/40960
             a_211^8 ),   _10^9 = F_2, ]
         
        which gives two solutions leading to a center at E_0, one of them 
        is a_211= 0, which yields the same solution as that 
        obtained in Case (A) at order 10. Thus, choosing the non-generic solution 
        a_211= 0 at this order leads to the generic pattern of Case (A) 
        starting from  ε^11 order (i.e., for k ≥ 11). 
        The second solution, given by 
        
            a_124 =1/64(10 ±1/10√(5710)) a_211^4,
          
        is a generic solution for Case (B), different from the 
        generic pattern of Case (A).  
        Then, following a similar computation procedure as that used in 
        Case (A), we obtain the generic solutions up to ε^39 order. 
        The results are given in Table <ref>, showing a 
        9^6-10^6-11^1 generic pattern, starting from order 14. 
        The notations used in this table are the same as that used in 
        Table <ref>.
        For each k, the key coefficient used to obtain the center condition is 
        give below. 
        
            [                                                                 9  LC:                                                                    k=7             a_037= - 35/768 a_211[
            64 a_123^2
            -a_211 (15 a_123 a_211^2;                                                                                                                                                                                   -16 a_124 a_211 -48 a_123 a_212) ];                                                                                                                                           k=9                                                               a_039= ⋯;                                                                [0.5ex]                                                                   k=14                                          a_0314= - 35/48 a_212^3 a_128;                                                                [0.5ex]                                                                   k=15                  a_0315= - 35/48 a_212^2 (a_129 a_212 
            +3 a_128 a_213);                                                                [0.5ex]                                                                   k=16    a_0316=  - 35/768 a_212[48 a_213^2 a_128  
            +a_212 ( 16 a_1210 a_212;                                                                                                                                                              +48 a_129 a_213
            +48 a_128 a_214
            - 15 a_128 a_212^2 
            ) ];                                                                [0.5ex]                                                                k=17–19                                                               a_03k= ⋯;                                                                [0.5ex]                                                                k=27–32                                                               a_03k= ⋯;                                                      [1.5ex] 
             10  LC:                                                                  k = 8                                                              a_123 = 0;                                                                                                                                          k=10                                      a_124 = 100 ±√(5710)/640  a_211^4;                                                                [0.5ex]                                                                   k=11 a_125 = 100 ±√(5710)/10240  
            a_211^3 [ 64 a_212 + a ( a^2-4) a_211^2 ];                                                                [0.5ex]                                                                   k=12                                                              a_126 = ⋯;                                                                [0.5ex]                                                                   k=20                                      a_128 = 100 ±√(5710)/640  a_212^4;                                                                [0.5ex]                                                                   k=21                                a_129 = 100 ±√(5710)/160  a_212^3 a_213;                                                                [0.5ex]                                                                k=22–25                                                         a_12(k-12) = ⋯;                                                                [0.5ex] ]
        
            [                                                                                                                                  k=33 a_1215 = - 100 ±√(5710)/10240 a_213[a_213^2 (15 a_213^2 -64 a_216);                                                                                                                                                                   - 64 a_214 ( a_214^2 +3 a_213 a_215) ];                                                            [0.5ex]                                                            k=34–38                                                     a_12(k-18) = ⋯;                                                            11  LC:                                                              k=13m                                                                   ;                                                                                                                              m= 1,2,3                                                          a_21m = 0 ]
        
        
         §.§ Non-generic solutions
         
        
        Couple of non-generic solutions have been discussed in Case (B) 
        (see Section 3.5), showing that setting a_211= 0 at order 
        8 or 10 (see Eqns. (<ref>) and (<ref>))
        leads to the 9^4-10^2-11^1 generic pattern of Case (A) 
        for orders greater than 10 or 11.
        These two examples give a route from Case (B) to Case (A).  
        In this section, we present several more non-generic solutions to 
        show other possibilities that they eventually return to either 
        the 9^4-10^2-11^1 generic pattern of Case (A) or 
        9^6-10^6-11^1 generic pattern of Case (B). 
        Other cases can be similarly discussed. 
        Since the discussions for different cases are similar, 
        we will not give the details but list the cases below and 
        summarize the results in Table <ref>. 
        
        
          (A1)  
        In Case (A), at order 13: a_212=0, leading to Case (A). 
        
        
          (A2)  
        In Case (A), at order 12: a_124=0, leading to Case (B). 
        
        
          (B1)  
        In Case (B), at order 11: a_211=0, leading to Case (B).  
        
         
        
        
        For each k, the key coefficient used to obtain the center condition is 
        given below. 
        
            [                                                      Case  (A1)                                                            k=13                                                        a_212= 0;                                                                                                                            k=14 a_0314= -35 a_125/48[ 4 a_125 a_214+ a_213 (8 a_126+ a_213^2) ];                                                         [0.5ex]                                                            k=16                 a_0316= - 35/48[ a_127 a_213 
            (8 a_126+a_213^2);                                                                                                                                                              +a_126 a_214 (4 a_126 +3 a_213^2)];                                                         [0.5ex]                                                            k=17                    a_0317= -35/48 a_128 a_213 (a_213^2+8 a_126);                                                                                                                                     -35/48 a_127 (4 a_213 a_127+ 8 a_214 a_126+ 3 a_214a_213^2);                                                                                                                                 - 35/48 a_126 (3 a_214^2 a_213+ 4 a_126 a_215+ 3 a_215 a_213^2);                                                         [0.5ex]                                                            k=15                                                        a_125= 0;                                                         [0.5ex]                                                            k=18                                            a_126= -9/8  a_213^2 ]
        
            [                                                               Case  (B1)                                                                     k=11                                                                 a_211= 0;                                                                                                                                              k=12    a_0312= -35/48 a_212 
            (a_126 a_212^2 + 3 a_125 a_212 a_213+ 4a_125^2);                                                                                                                                        k=13,15,17                                                                a_03k = ⋯;                                                                                                                                        k=14,16,18                                                         a_12(k/2-2) =  0;                                                     [1.5ex] 
             Case  (A2)                                                                     k=13 a_0313= -35/48 a_212[ a_127 a_212^2 
            +  a_126 (3 a_213 a_212+ 8 a_125) ];                                                                                                                                                                              +35/768 a_125 
            ( 15 a_212^4 -48 a_212^2 a_214;                                                                                                                                                                                        -48 a_212 a_213^2 -64 a_125 a_213 );                                                                                                                                           k=15,17                                                                a_03k = ⋯;                                                                                                                                     k=12,14,16,18                                                          a_12(k/2-2) = 0 ]
        
        
        Therefore, there are four possible routes for the non-generic solutions: 
        from Case (A) to Case (A) or Case (B); and 
        from Case (B) to Case (A) or Case (B). 
        
        
        
        
         §.§ Summary of this section
         
        
        Summarizing the results obtained in sections 
        3.4, 3.5 and 3.6, we have the following theorem. 
        
          
        For system (<ref>), based on the higher-order focus values, 
        there exist two generic patterns: One is 9^4-10^2-11^1 pattern 
        starting from order 8 with four consecutive 9 limit cycles, 
        followed by two consecutive 10 limit cycles, and then 
        one 11 limit cycles up to ε^35 order; 
        and the other is 9^6-10^6-11^1 pattern, 
        starting from order 14 with six consecutive 9 limit cycles, 
        followed by six consecutive 10 limit cycles, and then 
        one 11 limit cycles up to ε^39 order. 
        Other non-generic solutions deviate from the current pattern for 
        certain orders and eventually return to either 
        the 9^4-10^2-11^1 pattern or the 9^6-10^6-11^1 pattern.   
        
        
        Finally, we propose a conjecture on the number of limit cycles 
        around E_0 for system (<ref>). 
        
        
         Conjecture. For the perturbed system (<ref>), 
        the maximal number of small limit cycles which can bifurcate from the
        equilibrium E_0 is 11,
        
        
        
        
        
        § CONCLUSION
        
        
        In this paper, we have applied high-order focus value computation to 
        prove that system (<ref>)
        can have 11 limit cycles around the equilibrium of 
        (<ref>), obtained by 
        perturbing at least ε^7-order focus values. 
        Moreover, no more than 11 limit cycles can be found up to 
        ε^39-order analysis. It is believed that 
        system (<ref>) can have maximal 11 small limit cycles 
        around the equilibrium. 
        
        
        
        
        § ACKNOWLEDGMENT
        
        This work was supported by the National Natural Science Foundation of 
        China (NSFC No. 11501370), and the  
        Engineering Research Council of Canada (NSERC No. R2686A02). 
        
        
        
        
        
        
        
        
        
        
        
        
        
        § APPENDIX A
         
        
        The coefficients C_i's
        in (<ref>), (<ref>) and (<ref>) are given 
        below. 
        
            [                                       C_1 =                  - 3/16 (3 a^4 + 4 a^2 +44)           C_2 = - a/48 (3 a^4+12 a^2 + 116);                                       C_3 =                       - 1/8 (a^4 +2 a^2 +5)         C_4 = - 1/16 (9 a^4 - 20 a^2 + 172);                                       C_5 =                         - 3 a/8 (3 a^2 -8 )                   C_6 = - a/12 ( a^2 -120 );                                       C_7 =                           - 1/8 (3 a^2 -80)              C_8 = 3/128 a (7a^4+68a^2-900);                                       C_9 =       1/64(3 a^6 + 40 a^4 - 860 a^2 - 1600) C_10 = - a/576 ( 303 a^4 + 1596 a^2 + 8096);                                      C_11 =                      -  a/64 (55 a^2 - 256)      C_12 = - 1/384 (569 a^4-1660 a^2+1420);                                      C_13 =    1/192 (21 a^6 + 80 a^4 - 2996 a^2 -9040)           C_14 = 3/64 (7 a^4 + 8 a^2 + 160);                                      C_15 =              a/128 (49 a^4 + 220 a^2 - 380)                C_16 = a/32 (9a^4+36a^2-172);                                      C_17 =                    3/64 (15 a^4+28 a^2-916)         C_18 = - 1/16 (3 a^4 + 4 a^2 + 340) ]
        
            [                                      C_19 =                           - 3a/8 (a^2 - 16)         C_20 = a/96(45 a^4 + 266 a^2 - 644);                                      C_21 =                            a/32 (41 a^2+72)        C_22 = 1/192 (303 a^4+1728 a^2-3620);                                      C_23 =              -9/512 a ( a^4 + 12 a^2 + 200)                    C_24 = - 9/256 (a^4+160);                                      C_25 = - 1/512 (70 a^6 - 471 a^4 + 128 a^2 
            - 300)      C_26 = - 9/512 a (a^4 - 108 a^2 + 696);                                      C_27 =              3/256 (21 a^4 + 58 a^2 - 1880)              C_28 = a_125 a_212+a_124 a_213;                                      C_29 =                     a_215 a_212+a_214 a_213                      C_30 = 3 a_213 a_212^2 ]
        
        
        § APPENDIX B
         
        
        The coefficients s_i's involved in H_2 (see Eq. (<ref>)) 
        and t_i's involved in H_3 (see Eq. (<ref>)) are given 
        as follows. 
        
            [                                                                                  s_1 =                         1/3072a^6(3a^2+16)(10+4y+4x^2+3x^4)
            -1/192(2850+1824a^2-85a^4);                                                                                                  +a/18(6a^2 - 319)x
            -1/288(5902 + 1568a^2 + 45a^4)y
            -a/9(53+10x^2)xy
            +13/6y^2;                                                                                                       -1/72a^2(3a^2+4)y(y-x^2)
            -1/96(1074+200a^2+45a^4)x^2
            -29/18x^2y
            -2/9y^3;                                                                                                         +a/36(12-4a^2+3a^4)x^3
            -1/1152(6746+3140a^2+91a^4+3a^6)x^4
            +2/3x^2y^2;                                                                                  s_2 =                  1/128a^2(24-10a^2+5a^4)
            -1/64(1120-24a^2+ 2a^4 -a^6)(y+x^2)
            
            -1/2x^2y;                                                                                        +a/16(4 + 5a^2)x
            -25/4x^2
            -a/256(a^2 - 4)x [32(y - x^2) - 3a^3x^3]
            -1/8(73 - 2a^2 )x^4;                                                                                  s_3 =                             1/384a(-4716+376a^2+25a^4+15a^6)
            -1/96(2260-164a^2-15a^4)x;                                                                                                                +1/192a^5(5+3a^2)(y+x^2)
            -1/12a(a^2-4)y(y-x^2)
            -2/3xy^2
            +2x^3y;                                                                                                                     -1/48[a(1151+10a^2)y
            -a(1511+40a^2)x^2
            +(100-4a^2-3a^4)xy;                                                                                                                   +(140 -76a^2 +11a^4)x^3]
            -1/768a(10216+100a^2-2a^4-9a^6)x^4;                                                                                  s_4 =                                             1/9216a(96656+17952a^2+3640a^4+24a^6+9a^8);                                                                                                                 -5/4608(2256+7272a^2-56a^4+6a^6+9a^8)x
            -5/72a(224+8a^2+3a^4)y;                                                                                                                   -5/144a(642-8a^2-3a^4)x^2
            -5/1152(368-1628a^2-92a^4+3a^6)xy;                                                                                                                     -5/2ay^2
            +5/1152(5272+86a^4+2964a^2-3a^6)x^3
            -275/36ax^2y;                                                                                                                       +5/288(372 +4a^2 +3a^4) xy^2
            -5/144a(140+12a^2+3a^4)x^4;                                                                                                                     -5/288(268+3a^4+4a^2)x^3y
            -5/18xy^3-25/18ax^4y
            +5/6x^3y^2;                                                                                  s_5 =                               1/256a(-16a^2-8a^4+a^6+1120)
            +5/256(560+9a^4-2a^6-4a^2)x;                                                                                                             -5/64a(20-11a^2)x^2
            -5/128a^2(3a^2+4)x(y+x^2)
            +175/8xy
            +265/16x^3;                                                                                                                                              +5/32a(a^2-4)x^2(y-x^2)
            -5/8x^3y;                                                                                        ]
        
            [                                                                            s_6 =    1/768a^2(540a^2 + 3904 - 8a^4 + 3a^6)
            -5/1536a(484 - 100a^2 - 23a^4 + 12a^6)x;                                                                                  -5/12a^2(a^2 - 1)y
            -35/384(244 - 4a^2 - 7a^4)x^2
            +5/768a(5132 + 84a^2 - 13a^4)xy;                                                                                                           -5/768a(-3756-76a^2+13a^4)x^3
            +5/192(380+3a^4+4a^2)x^2y;                                                                                              +5/48a(a^3-4)xy(y-x^2)
            -5/192(11a^4+20a^2+12)x^4
            -5/6x^2y^2+5/2x^4y. ]
        
            [                                                                          t_1 =                     1/384a(-4716+376a^2+25a^4+15a^6)
            -1/96(2260-164a^2-15a^4)x;                                                                                              +1/192a^5(5 +3a^2)(y+x^2)
            -1/48a( 1151 +10a^2 )y
            -2/3xy^2
            +2x^3y;                                                                                    -1/48a( 1511 + 40a^2 )x^2
            +1/48(100 - 4a^2 - 3a^4)xy
            -1/12a(a^2-4)y(y-x^2);                                                                                              +1/48( 140 -76a^2 +11a^4)x^3
            +1/768a(-100a^2-10216+9a^6+2a^4)x^4;                                                                          t_2 =                 1/64a^2( 24 -10a^2 +5a^4)
            -1/32(1120 -24a^2 +2a^4 -a^6)(y+x^2);                                                                                                            +1/8a(5a^2+4)x
            -25/2x^2
            -1/4a(a^2-4)x(y-x^2)
            -x^2y;                                                                                                                           -1/128(2336 -64a^2 +12a^4 -3a^6)x^4;                                                                          t_3 =                         -1/2048a^5(20+a^4)(5+2y+2x^2+2x^4)
            -1/128a(550-343a^2);                                                                                -1/256(6800 - 24a^2 - 10a^4 + 5a^6)x
            -3/64a(25a^2 + 34)y
            -1/32a(45a^2 + 61)x^2;                                                                                            -1/1024a^4(2-a^2)x(8y-8x^2+3ax^3)
            -3/16(50-a^2)xy
            +3/16(80-a^2)x^3;                                                                                                             +1/32a(a^2-4)x^2y
            -1/64a( 197 -2a^2 )x^4
            +1/4x^3y;                                                                          t_4 =  1/768a^2(3904 + 540a^2 - 8a^4 + 3a^6)
            -5/1536a(484 - 100a^2 - 23a^4 + 12a^6)x;                                                                                                                    -5/12a^2(a^2-1)y
            -35/384(244-7a^4-4a^2)x^2;                                                                                                    -5/768a(-5132-84a^2+13a^4)xy
            -5/768a(-3756-76a^2+13a^4)x^3;                                                                                ]
        
            [                                                                                                                +5/192(380+4a^2+3a^4)x^2y
            +5/48a(a^2-4)xy(y-x^2);                                                                                                                   -5/192(12+20a^2+11a^4)x^4
            -5/6x^2y^2
            +5/2x^4y;                                                                           t_5 =                         1/128a(1120-16a^2-8a^4+a^6)
            +5/128(560-4a^2+9a^4-2a^6)x;                                                                                                           +5/32a(11a^2-20)x^2
            -5/64a^2(4+3a^2)x(y+x^2)
            +175/4xy;                                                                                                                      +265/8x^3
            +5/16a(a^2-4)x^2(y-x^2)
            -5/4x^3y;                                                                           t_6 = -1/2048a^2(2544 - 2336a^2 + 8a^4 + 3a^6)
            -5/2048a(2632 + 276a^2 + 3a^4 - 8a^6)x;                                                                                                    +5/512(1240-4a^2-41a^4+4a^6)x^2
            +5/1024a(-1832+892a^2+a^4)xy;                                                                                                       +5/1024a(-1368+696a^2+a^4)x^3
            -5/256a^2(3a^2+4)x^2(y-x^2);                                                                                                              +375/32x^2y
            -65/16x^4
            -5/128a(a^2-4)x^3y
            +5/16x^4y ]
        99tocsectionBibliographyHan2013
        M. Han, Bifurcation theory of limit cycles, Science Press, Beijing, 2013.
        
        Hilbert1900
        D. Hilbert,
        Mathematical problems,
        (M. Newton, Transl.) Bull. Amer. Math. 8 (1902) 437–479.
        
        Shi1980
        S. Shi,
        A concrete example of the existence of four limit cycles for
        plane quadratic systems,
        Sci. Sinica, 23 (1980) 153–158.
        
        ChenWang1979
        L. Chen, M. Wang,
        The relative position, and the number, of limit cycles of a
        quadratic differential system,
        Acta. Math. Sinica, 22 (1979) 751–758.
        
        LiLiu2010
        J. Li, Y. Liu,
        New results on the study of Z_q-equivariant planar polynomial
        vector fields,
        Qualitative Theory of Dynamical Systems, 9(1-2) (2010) 167–219.
        
        LLY09
        C. Li, C. Liu, J. Yang,
        A cubic system with thirteen limit cycles,
        J. Diff. Eqns. 246 (2009) 3609–3619.
        
        Arnold83
        V. I. Arnold,
        Geometric Methods in the Theory of Ordinary Differential Equations,
        Springer-Verlag, New York, 1983.
        
        Han2006
        M. Han,
        Bifurcation of limit cycles of planar systems,
        Handbook of Differential Equations, Ordinary Differential Equations,
        Vol. 3 (Eds. A. Canada, P. Drabek and A. Fonda), Elsevier, 2006.
        
        HanYu2012
        Han M., Yu, P.,
        Normal Forms, Melnikov Functions and Bifurcations of Limit Cycles.
        Sringer-Verlag, New York, 2012.
        
        Bautin1952
        N. Bautin,
        On the number of limit cycles appearing from an equilibrium point of the
        focus or center type under varying coefficients,
        Matem. Sb. 30 (1952) 181–196.
        
        YC2009
        P. Yu, R. Corless,
        Symbolic computation of limit cycles associated with Hilbert's 16th problem,
        Communications in Nonlinear Science and Numerical Simulation, 14(12)
        (2009) 4041–4056.
        
        
        Lloyd2012
        N. Lloyd, J. Pearson,
        A cubic differential system with nine limit cycles,
        Journal of Applied Analysis and Computation 2 (2012) 293–304.
        
        
        CCMYZ2013
        C. Chen, R. Corless, M. Maza, P. Yu, Y. Zhang,
        A modular regular chains method and its application to dynamical systems,
        Int. J. Bifurcation and Chaos 23(9) (2013) 1350154 (21 pages).
        
        YuTian2014 
        P. Yu, Y. Tian, 
        Twelve limit cycles around a singular point in a planar cubic-degree
        polynomial system,
        Commun. Nonlinear Sci. Numer. Simulat. 19 (2014) 2690–2705.
        
        
        
        
        GuckenheimerHolmes1993
        J. Guckenhermer, P. Holmes,
        Nonlinear Oscillations, Dynamical Systems, and Bifurcations of
        Vector Fields (4th Ed.).
        Springer, New York, 1993.
        
        ChowLiWang1994
        S. N. Chow, C. C. Li, D. Wang,
        Normal Forms and Bifurcation of Planar Vector Fields.
        Cambridge University Press, Cambridge, 1994.
        
        Kuznetsov1998
        Yuri A. Kuznetsov,
        Elements of Applied Bifurcation Theory (2nd Ed.).
        Springer, New York, 1998.
        
        GazorYu2012
        M. Gazor, P. Yu,
        Spectral sequences and parametric normal forms.
        J. Diff. Eqns. (2012) 252: 1003-1031.
        
        Yu1998
        P. Yu,
        Computation of normal forms via a perturbation technique,
        J. Sound and Vib. 211 (1) (1998) 19–38.
        
        YuLeung2003
        P. Yu, A. Y. T. Leung,
        The simplest normal form of Hopf bifurcation.
        Nonlinearity (2003) 16: 277-300.
        
        TianYu2013
        Y. Tian, P. Yu,
        An explicit recursive formula for computing the normal form and
        center manifold of n-dimensional differential systems
        associated with Hopf bifurcation.
        Int. J. Bifur. Chaos (2013) 23: 1350104 (18 pages).
        
        TianYu2014
        Y. Tian, P. Yu,
        An explicit recursive formula for computing the normal forms associated
        with semisimple cases.
        Commun Nonlinear Sci. Numer. Simul. (2014) 19: 2294-2308.
        
        Zoladek1995
        H. Żoła̧dek,
        Eleven small limit cycles in a cubic vector field,
        Nonlinearity 8 (1995) 843–860.
        
        YuHan2011
        P. Yu, M. Han, 
        A study on Żoła̧dek's example, 
        J. Appl. Anal. Comput. 1 (2011) 143–153. 
        
        TianYu2016
        Y. Tian, P. Yu, 
        Bifurcation of ten small-amplitude limit cycles by perturbing a quadratic 
        Hamiltonian system with cubic polynomials, 
        J. Diff. Eqns. 260 (2016) 971–990.  
        
        
        
        
        
        § APPENDIX B
        
            [                                                                                       a_0311=                   -35/48[ a_125 a_212 (8 a_124+a_212^2)  
            + a_124 a_213 (4 a_124+3 a_212^2) ];                                                                                      a_0317 =                                             -35/768 a_213[ 
            32 a_127 (2 a_127 -3 a_214 a_213);                                                                                                                                         - a_213^2 (128 a_128 +54 a_214^2-27 a_215 a_213) ],;                                                                                      a_0318 =                               35/6 a_129 a_213^3 
            - 35/24 a_213 a_128 (4 a_127-3 a_214 a_213);                                                                                                                              - 35/48 a_127[ a_214 (4 a_127 +3 a_214 a_213)-6 a_213^2 a_215];                                                                                               + 3/4096 a_213^2 
            [ 1120 a_214 (a_214^2+6 a_215 a_213)
            +3 a_213^2 (2269 a_213^2-560 a_216) ], ]
        James1991
        E. M. James, N. G. Lloyd,
        A cubic system with eight small-amplitude limit cycles,
        IMA, J. of Appl. Math. 47 (1991) 163–171.
        
        Ilyashenko1991
        Yu. S. Ilyashenko, S. Yakovenko,
        Finitely smooth normal forms of local families of diffeomorphismes and
        vector fields,
        Russ. Math. Surv. 46 (1991) 3–19.
        
        Ecalle1992
        J. Écalle,
        Introduction aux Fonctions Analysables et Preuve Constructive
        de la Conjecture de Dulac (French),
        Hermann, Paris, 1992.
        
        YH2012
        P. Yu, M. Han,
        Four limit cycles from perturbing quadratic integrable systems by
        quadratic polynomials,
        Int. J. Bifurcation & Chaos, 22(10) (2012) 1250254 (28 pages).
        
        Li2003
        J. Li,
        Hilbert's 16th problem and bifurcations of planar polynomial vector fields,
        Int. J. Bifurcation and Chaos  13 (2003) 47-106.
        
        HorozovIliev1994
        E. Horozov, I. D. Iliev,
        On the number of limit cycles in perturbations of 
        quadratic Hamiltonian systems,
        Proc. Lond. Math. Soc. 69 (1994) 198–224.
        
        Iliev2000
        I. D. Iliev,
        On the limit cycles available from polynomial perturbations of the
        Bogdanov-Takens Hamiltonian,
        Israel J. of Math. 115 (2000) 269–284.
        
        Gavrilov2001
        L. Gavrilov,
        The infinitesimal 16th Hilbert problem in the quadratic case,
        Invent. Math. 143 (2001) 449–497.
        
        HorozovIliev1998
        E. Horozov, I. D. Iliev,
        Linear estimate for the number of zeros of Abelian integrals with cubic
        Hamiltonians,
        Nonlinearity 11 (1998) 1521–1537.
        
        Petrov1986
        G. S. Petrov,
        Elliptic integrals and their non-oscillatorness,
        Funct. Anal. Appl. 20 (1986) 46–49 (in Russian).
        
        Gavrilov1998
        L. Gavrilov,
        Nonoscillation of elliptic integrals related to cubic polynomials with
        symmetry of order three,
        Bull. Lond. Math. Soc. 30 (1998) 267–273.
        
        LiZhang1995
        B. Li, Z. Zhang,
        A note on a result of G. S. Petrov about the weakened 16th Hilbert problem,
        Journal of Mathematical Analysis and Applications 190 (1995) 489–516.
        
        Francoise1996
        J. P. Françoise,
        Successive derivatives of a first return map, application to the study of
        quadratic vector fields,
        Ergodic Theory and Dynam. Syst. 16 (1996) 87–96.
        
        Jebrane2007
        A. Jebrane, H. Żoła̧dek,
        A Note on higher order Melnikov functions,
        Qualitative Theory of Dynamical Systems 6 (2005) 273–287.
        
        
        Uribe2006
        M. Uribe,
        Principal Poincaré-Pontryagin function of polynomial perturbations of the Hamiltonian triangle,
        J. Dyn. Control Syst. 12 (1) (2006) 109–134
        
        Arnold1988
        V. I. Arnold, A. N. Varchenko, S.M. Gusein-Zade,
        Singularities of Differentiable Maps,
        v. 2, Monographs in Mathematics, Birkhäuser, Boston, 1988.
