New lower bounds for the Hilbert numbers using reversible centers

In this paper we provide the best lower bounds, that are known up to now, for the Hilbert numbers of polynomial vector fields of degree N, , for small values of N. These limit cycles appear bifurcating from symmetric Darboux reversible centers with very high simultaneous cyclicity. The considered systems have, at least, three centers, one on the reversibility straight line and two symmetric outside it. More concretely, the limit cycles are in a three nests configuration and the total number of limit cycles is at least 2n  +  m, for some values of n and m. The new lower bounds are obtained using simultaneous degenerate Hopf bifurcations. In particular, and


Introduction
We consider two-dimensional differential systemṡ x = P(x, y),ẏ = Q(x, y), (1) in which P and Q are polynomials of degree N. The maximum possible number, H(N), of limit cycles of system (1) is known as the Hilbert number. As usual we define limit cycle as every isolated periodic orbit.
As it follows from former definition, the Hilbert number refers to the total amount of limit cycles that system (1) exhibits and, in this sense, it is a global concept. For instance, Shi (see [22]) and Chen and Wang (see [4]) proved that H(2) 4. For cubics, Li, Liu, and Yang (see [17]) and Li and Liu (see [18]) proved that H(3) 13. In our work we prove, among other new best lower bounds for the Hilbert numbers, that H(4) 28.
The aim of this work is to obtain new lower bounds values for H(N). This goal is attained by simultaneously perturbing some reversible centers. We proceed studying simultaneous degenerate Hopf bifurcations of reversible centers and, to overcome heavy computations, we use an efficient way (parallelization method) for the Lyapunov quantities calculation. This work strongly relies on the results of Christopher, [6], and Han, [11] and [12]. The idea is to obtain good lower bounds for the number of small limit cycles bifurcating from an equilibrium point of center-focus type. This is done from the series expansions of the Lyapunov quantities at points on the center variety. In particular, if the first r linear terms of the Lyapunov quantities are independent, then the Hilbert number is bigger than or equal to r, if we consider also the trace as another independent parameter. The germ of this linearization idea can be glimpsed at the work of Chicone and Jacobs, see [5], where, for higher order limit cycles bifurcations, the Lyapunov quantities are obtained from the Taylor series expansion of the perturbation parameters. We will use this bifurcation phenomena simultaneously in a family with more than one center. We want to remark that, in the paper [6], Christopher studies this symmetric simultaneous bifurcation, only up to first order perturbation, of a quartic system obtaining 22 limit cycles in two symmetric nests of 6 cycles and one nest of 10. For short, we say that these limit cycles are in configuration 6,10,6 . In section 3 we give the precise definition of the configuration concept and we describe the complete bifurcation mechanism.

Theorem 1. H(4)
The linear computational approach to estimate the above lower bounds has, at least, two advantages: The first one is computational, and it is based on the fact that, by shortening the length of the expressions to be manipulated, more fast computations can be performed. The second one has to do with the functional independence argument of the Lyapunov quantities, since this argumentation is replaced by a study, in most of the cases, on the independence of their linear parts.
Nevertheless, sometimes, linear terms of the Lyapunov quantities are not independent and, for studying the simultaneous bifurcation, we must take into account higher order terms. This disadvantage was considered by Christopher and, to overcome the setback that the computations are no longer linear and become unmanageable, under some generic assumptions, something specific can be proved, see [6, theorem 3.1].
The improvement for degree 4, with respect to Christopher's work, is that we obtain 28 limit cycles in two nests of 8 cycles and one nest of 12, i.e. in a 8,12,8 configuration. The key point in obtaining this new best estimate of H (4) is the use of fifth order, instead of linear order, approximation of the Lyapunov quantities. See section 5. It is also worth mentioning that, despite having tried to use this technique in the cubic case, we have not been able to improve the current lower bound of H (3). This is done in section 4.
Concerning the parallelization method, roughly speaking, it is as follows. We, one-byone, compute the linear part of a fixed number of Lyapunov quantities of simpler differential equations having only one perturbation monomial. Then, by taking the sum of all the results, we obtain the corresponding linear part of the Lyapunov quantity of the complete differential equation. If we proceed in this way it can be checked that the computations in each simpler equation are shorter in size and time. See [19] for more details.
In our paper, we follow former approach and latter technique applied to families of polynomial reversible symmetric differential systems (1) having two centers outside the symmetry line, which will be the x-axis, plus one more extra center at the origin. See figure 1. Then, we force the existence of two symmetric weak foci of order n where the symmetric centers were located, outside of the x-axis. We aim to create n limit cycles, by using only reversible perturbations, around each one of them; in this way, 2n simultaneous limit cycles will be achieved. After that, by perturbing again the system, but without the reversibility property, we can generate m more limit cycles surrounding a weak focus of order m at the origin. Summing up, the total number of limit cycles we have generated is 2n + m, in configuration n, m, n . The proof of theorem 1 follows from propositions 6-8, while theorem 2 is proved in section 8. We note that the supplied lower bounds for the Hilbert numbers, H(N), for N = 4, . . . , 10, in theorem 1 significantly improve the previous best known Hilbert numbers. See for instance [13,14,19]. We also give a new general lower bound for these Hilbert numbers. See theorem 2.
Even taking into account that the parallelization method has more facilities than using other ways to proceed in finding new better lower bounds of H(N), some obstacles are found in to go ahead for values of N higher than those considered. Next, we comment some of them. The first obstruction concerns the computer memory. Despite of the usage of linearizations, the computations of Lyapunov quantities for higher values of N requires, increasing N, more and more memory. Another arising problem is the necessity to get reversible differential systems providing good seeds for bifurcating the required number of limit cycles. Finally, another problem, for our computational approach, is how to obtain good first integrals in such a way that, when writing it in a complex normal form, the coefficients in the differential equations are all rational numbers. This paper is structured as follows. In section 2 we recall, not only some of the known best estimates of lower bounds for the Hilbert numbers, but also the corresponding known best lower bounds for the cyclicity of non-degenerated monodromic equilibrium points. In section 3 we describe the methodology used to get the main results. In section 4 we provide a detailed description of the method applied to a cubic vector field, getting a characterization of all cubic systems having a reversible center with two extra symmetric centers outside the symmetry axis. In section 5, we study some reversible quartic systems and, by considering fifth order terms in the Lyapunov quantities Taylor series, and we get the new best lower bound H(4) 28. Section 6 is devoted to present, among other lower bounds for the Hilbert numbers, the new best lower bounds H(5) 37 and H(6) 53. In section 7 we obtain the new best lower bounds H(7) 74, H(8) 96, H(9) 120, and H(10) 142. Finally, in section 8 we prove theorem 2.

Previous lower bounds of H(N)
Before to start with the global concept of the Hilbert number, a local concept concerning the bifurcation of small limit cycles around a specific equilibrium point is also considered. This is the case for the cyclicity of an equilibrium point.
Roughly speaking, a monodromic equilibrium point for a family of vectors fields in (1) has finite cyclicity M(N) if, in Hausdorff distance, the maximum number of limit cycles (taking into account its multiplicity) is M(N), and they appear in a close neighbourhood of the equilibrium point, when the perturbation tends to zero whereas the neighbourhood shrinks to the equilibrium point. From previous definitions it is easy to conclude that H(N) M(N), for all N.
In this section we present the known lower bounds, up to now, of M(N) and H(N). As we prove in this work, some of them will be improved from our results.
A classical bifurcation mechanism to study the existence of periodic orbits around an equilibrium point of monodromic type is the non-degenerate Hopf bifurcation method. This mechanism studies the limit cycles that emerge from an equilibrium point of center-focus type when, by changing the sign of the trace in a suitable way, its stability changes from stable to unstable or vice-versa. The non-degenerate Hopf bifurcation method can be generalized, by computing the Lyapunov quantities, to more degenerate critical points. This generalization is known as the degenerated Hopf bifurcation method and it is used to give estimates from below of M(N), see [1] for instance. Besides former reference, see also [12,29] for more results on bifurcation theory of limit cycles. Among these specialized references, a general approach to qualitative theory of differential equation can be found in [8].
In the last decades some improvements of the lower bounds of H(N) have been achieved by using better estimates of M(N). In fact, this is the case for some special small values of N. At this point it is worth to mention that only for quadratic polynomial vector fields we know the exact value of it, that is M(2) = 3. This fact was proved by Bautin in [2]. Concerning M(N) for N 3, up to now, only lower bounds for these numbers has been obtained, as detailed below.
In finding lower bounds of H(N) it is worth mentioning that one of the most used techniques to create limit cycles is by perturbing Z q -equivariant polynomial Hamiltonian vector fields. Frequently, this technique is combined with other procedures such as: Hopf bifurcation method, homoclinic and heteroclinic bifurcation methods, Abelian integrals and bifurcation from infinity, ..., see [8,12,21].
The known best result on M(3) was obtained by Żołądek in [30], getting M(3) 11. This work has been amended in [32] although the same bound for M(3) is maintained. The proof was made, after a cubic perturbation of a cubic differential equation having a rational first integral, by bounding the zeros of the displacement map. This displacement map, in first approximation, is given by a linear Poincaré-Pontriaguin integral which is expressed in terms of twelve Abelian or Melnikov type integrals (first order Melnikov integrals) generating the aforementioned eleven limit cycles. Linked with the first paper of Żołądek, in 2005, Christopher in [6] confirmed the lower bound for cubic polynomial vector fields, that is M(3) 11. He uses a family of cubic systems, C (12) 31 in accordance with Żołądek's more recent classification [31], having also a rational first integral and showing that the linear part of the first Lyapunov quantities have maximal rank 11. In [3], Bondar' and Sadovski also proved, with the same technique but with a different Darboux center, that M(3) 11. In the case N = 4, that is on quartics systems, in [10], Giné proved that M(4) 21, by studying small limit cycles bifurcating from an elementary center-focus type equilibrium point. In the case N = 5, i.e. concerning lower bounds of M(5), the known best estimate is the one given by Giné in [10], where M(5) 26 was proved. These limit cycles also bifurcate from an elementary center-focus type equilibrium point at the origin. In this paper, Giné conjectured that the number of functionally independent focal values, i.e. the minimum number of ideal generators of the ideal generated by the focal values, of system (1) with an elementary center-focus type equilibrium point at the origin will be N 2 + 3N − 7. From this value he also conjectured that M(N) = N 2 + 3N − 7. We want to point out that, regarding Giné's conjecture, we would have M(4) = 21 and M(5) = 33.
In the cases N = 6, . . . , 13, Liang and Torregrosa in [19] proved that the cyclicity of some holomorphic Darboux centers for system (1) is, at least, N 2 + N − 2, for 4 N 13. Hence, we have the lower bound H(N) M(N) N 2 + N − 2, for 4 N 13. We point out that in that work, the authors think that the best center candidate to produce high cyclicity is of Darboux type. One consequence of their paper is that this result gives the highest known up to now lower bound of M(6), M(7), . . . , M (13), and another is that these numbers are smaller than those conjectured before. In particular, for N = 6, 8, 10 are also the best lower bound for Hilbert numbers, i.e. H(6) 40, H(8) 70, and H(10) 108.
The general lower bound M(N) N 2 − 2 was done by Movasati, see [20, corollary 4.1]. The author study the perturbation of holomorphic centers in the class of polynomial vector fields of degree N. Showing that there are system with not less than N 2 − 2 limit cycles. In regards to upper bounds of M(N), which it is a much more difficult question, up to our knowledge there are no results.
Next we summarize some of the recent progress in finding lower bounds of H(N). On H(2), Shi (see [22]) and Chen and Wang (see [4]) proved that H(2) 4. They obtained these four limit cycles in two nests, one with 3 and another with 1.
On H(3), in 2009, Li, Liu, and Yang (see [17]), by counting the number of zeros of some Abelian integrals, demonstrated that an specific planar cubic system has H(3) 13. Later, Li and Liu (see [18]) also proved that H(3) 13, using Z 2 -equivariant cubic perturbations, Hopf bifurcation and changing the stability of infinity. These 13 limit cycles appear in two nests of 6 and one coming from infinity.
About H(4), Christopher (see [6]) in 2005 provided a quartic system with H(4) 22. As we have mentioned before, he used the linear parts of the Lyapunov quantities for center bifurcations with symmetries, the limit cycles appear in two nests of 6 and one nest of 10. Recently, in 2011, Johnson constructed (see [14]) a quartic system that has 24 limit cycles, using Lyapunov quantities, in four nests of 6. Moreover, the existence of 2 more, having each one inside two of these nests, is shown numerically. Then, H(4) 26. In the study of H(5), in 2008, Wu, Gao, and Han in [27] found that H(5) 28, in four different configurations, by using the perturbation of a Z 2 -equivariant symmetric quintic Hamiltonian system. Later, in 2010, Wu, Wang, and Tian, in [28] proved also H(5) 28 with a Z 4 -equivariant symmetric quintic vector field, by using simultaneously Hopf and polycycle bifurcation methods. Also in 2010, Johnson and Tucker, in [16], studied the limit cycle bifurcation of a Z 2 -equivariant quintic planar Hamiltonian vector field under Z 2 -equivariant quintic perturbation and, by computer aided approach, proved that H(5) 27. In 2015, Sun and Han in [23] also proved that H(5) 28 using Z 4 -equivariant perturbations.
Concerning H(6), in 2005, Wang and Yu, showed that H(6) 35, by using global and local bifurcations, see [25]. In this paper the authors perturb a fifth-degree Z 2 -equivariant symmetric Hamiltonian system by adding a sixth-degree polynomial system. More recently, in 2015, Liang and Torregrosa, by computing the required independent linear parts of the Lyapunov quantities, proved that H(6) M(6) 40, see [19].
On H(N), for N 7, it is worth mentioning the relevant paper of Han and Li [13] in 2012. The results in this work improved almost all existing lower bounds of H(N) for such values of N, up to that moment. More concretely, they proved that H(7) 65, H(9) 98, H(11) 153, H(12) 157, . . . In 2010, Johnson and Tucker in [15] proved that H(7) 53, perturbing a Z 2 -equivariant planar Hamiltonian vector field of degree 7. The proof follows a computer aided approach. For N = 11, before the paper of Han and Li, the best one was the one given by Wang and Yu in 2006, in [26], proving that H(11) 121 using a Z 12 -equivariant vector field. With respect to H(8), until 2012 the best lower bound for this Hilbert number was the one given by Han and Li,as H(8) 67. Nevertheless, in 2015, Liang and Torregrosa, in [19], obtained a new better one, that is H(8) M(8) 70. In this paper it was also proved that H(10) M(10) 108.
Concerning the growth rate of H(N) as N increases, one of the first results in finding a general lower bound we want to point out, is the work of Christopher and Lloyd in 1995, see [7], where they introduced a recurrence mechanism to provide a lower bound for H(N) that grows like KN 2 log N, for some real constant K. In fact, they prove that Recently, in 2012, Han and Li improved this mechanism by providing a better lower bound of the same kind, see [13]. Concretely, for N big enough, they showed that H(N) grows at least as fast as (N + 2) 2 log(N + 2)/ log 4.

The bifurcation mechanism
The bifurcation mechanism is as follows. First, we consider a reversible system (symmetric with respect to the x-axis) having a center at the origin and two more at two symmetric points (x c , ±y c ), with x c y c = 0, located at the vertices of a triangle. Second, we consider a first perturbation obtaining simultaneously n hyperbolic limit cycles in each one of the two nests while the origin remains as a center. Finally, we consider a second perturbation of the system giving us m limit cycles, around this third nest, increasing the total number of limit cycles obtained. Then, we will say that the 2 × n + m limit cycles appear in configuration n, m, n .
The limit cycles appear using degenerated Hopf bifurcations of some order k and all the unperturbed polynomial systems have rational first integrals. Concerning the Lyapunov quanti ties computation we use linear or higher order developments of them, as in the works of Christopher ( [6]) and Han ([11]). In this paper, reversible differential systems in the (x, y)-plane means invariant under the change of variables (x, −y, −t) → (x, y, t). Moreover, all the centers are non-degenerate.
As all the computations are quite hard to do and the obtained expressions are too big to be written here, we will describe only the bifurcation mechanism more in detail. This bifurcation scheme is the same followed by Christopher (see also [6]) but we use higher instead of first order. In some simple cases and for lower degrees we will also show which are the Taylor series of the Lyapunov quantities.

Degenerated Hopf bifurcation
An analytic system with a non-degenerate equilibrium point of center-focus type can be transformed to the system where f and g are analytic functions with lower order terms of degree 2 in x, y, being λ a real parameter and with the equilibrium point located at the origin. When λ = 0 we have a weak focus or a center at the origin. The problem to distinguish its stability is the well-known center-focus problem. There are some classical approaches to solve it, see [1] for instance. Here we will use the Andronov-Poincaré method, that is the computation of the obstruction conditions for its integrability, but using complex notation. Hence, following the approach of [24], system (2) writes, in complex coordinates (z = x + iy), as where F j are homogeneous polynomials of degree j in z,z. Now, if λ = 0, we look for the existence of a formal first integral of the form with H j homogeneous polynomials of degree j in z,z. This can be done imposing that the level curves of H contain solutions of equation (3). Straightforward computations show that the coefficients of H 3 can be uniquely determined in terms of the coefficients of F 2 , by solving the diagonal linear system of equations obtained vanishing all the coefficients in z,z of the homogeneous polynomial of degree 3 We can not follow exactly in the same way, to fix the coefficients of H 4 , because the linear system of equations obtained vanishing the coefficients of the homogeneous polynomial of degree 4,z has not maximum rank five. This obstruction can be removed by adding an extra term, L 1 z 2z2 for example, in the above homogeneous identity. This constant L 1 is known as the first Lyapunov quantity. Clearly, H is a Lyapunov function when L 1 = 0 and gets the stability of the origin. When L 1 = 0, we need to compute higher degree terms of the function H for distinguishing when the equilibrium point is a center or a focus. This constant L 1 can be thought as the first obstruction for (3) to be integrable. In the recursive procedure to get H j , for j = 3, . . . , it can be proved that all the linear systems giving the coefficients of H 2k+1 , k = 1, 2, . . ., have maximum rank. Consequently H 2k+1 can be uniquely determined. But in each even step, i.e. to determine H 2k+2 , k = 1, 2, . . ., we need to add a correction term of the form L k z k+1zk+1 . These quantities L k are the well known Lyapunov quantities associated to equation (3). Clearly, the first non-vanishing Lyapunov quantity, L k , gives the stability of the origin. In this case, as usual, we say that the origin is a weak focus of order k. We would remark that the main computational difficulties to obtain the Lyapunov quantities, applying this procedure or in general other methods, are only because of their huge expressions. For polynomial systems, the Lyapunov quantities can be used, besides to determine the stability of the origin and the center conditions of system (2), to obtain periodic orbits from the origin. For example when L 1 < 0 and λ = 0 the origin is stable and for λ > 0, small enough, the stability of the origin changes to unstable and an small stable hyperbolic periodic orbit bifurcates from the origin. This is the well-known Hopf bifurcation. An analogous procedure can be used to create k limit cycles from a weak focus of order k by a degenerate Hopf bifurcation. These limit cycles are also known as small limit cycles. See below. Now we will describe how the small limit cycles appear in a degenerate Hopf bifurcation. Suppose that system (2) is a polynomial system, that we call system S, having, without loss of generality, a stable weak focus of order k at the origin. We will think the trace parameter λ as the zero-order Lyapunov quantity. Hence, we write L 0 = λ. Then, the Lyapunov quantities satisfy L 0 = L 1 = · · · = L k−1 = 0 and L k < 0. Let Γ be a level curve of H which is sufficiently near the origin that the flow is inward across it. Now we perturb system S in such a way that the perturbed system, we call it system S 1 , has L 0 = L 1 = · · · = L k−2 = 0 and L k−1 > 0. Below we comment on the conditions ensuring this is possible. The origin is thus unstable. If S 1 is sufficiently near S, the flow remains inward across Γ and, if H 1 is the Lyapunov function corresponding to S 1 , there exists a level curve Γ 1 of H 1 inside Γ and sufficiently near the origin that the flow is outward across Γ 1 . By the Poincaré-Bendixson theorem, there is a limit cycle between Γ and Γ 1 . The next step is to take a perturbation, S 2 , of system S 1 so that L 0 = L 1 = · · · = L k−3 = 0 and such that L k−2 < 0. In this case, if the perturbation is sufficiently small, the flow remains inward across Γ and outward across Γ 1 . Hence, S 2 has a limit cycle between Γ and Γ 1 and, since the origin is stable for S 2 , there is also a limit cycle inside Γ 1 . Proceeding in this way, k limit cycles can be generated. The necessary conditions are The above necessary conditions are not simple to be proved in a fixed family of polynomial systems (2) of a given degree. The existence of parameters ensuring the independence of the Lyapunov quantities to get exactly k periodic orbits near a weak focus of order k is only guarantied for general analytical perturbations, where always a weak focus of order k unfolds k hyperbolic limit cycles, see [21].

Cyclicity of centers
In this paper we select polynomial center families of degree N having cyclicity at least k. We recall that the cyclicity of an equilibrium point can be defined as the maximum number of isolated periodic orbits (taking into account its multiplicity) bifurcating from it. In our approach this is equivalent to say that, in the parameter space, near the origin and after perturbing the center with polynomials also of degree N, there will be a curve of weak foci of order k passing though the origin and unfolding k hyperbolic periodic orbits. We remark that, in the parameter space, the origin corresponds to the case in which the system, not perturbed, has a center. As we have explained before, the trace parameter, λ, can be seen as the zero-order Lyapunov quantity, L 0 , and we always obtain a limit cycle with it. Consequently, we will restrict our analysis to the class of zero trace perturbations. That is, to the perturbations whose lower degree terms are of degree two. As we have mentioned before, to get the existence of these curves in the parameter space and the corresponding unfolding from the Taylor developments, we follow the ideas of Christopher (see [6]) and Han (see [11]). The existence of a curve of weak focus is obtained studying the intersection of the Taylor developments of the varieties L j = 0, for j = 1, . . . , k − 1 ensuring that along it L k is non-vanishing. The unfolding of k limit cycles (adding the trace parameter) and their hyperbolicity is obtained from the transversality of the intersection of the varieties. Clearly, using the implicit function theorem, if the matrix defined by the linear terms of L 1 , . . . , L k have rank k we have (adding λ) at least k hyperbolic limit cycles bifurcating from the center point. Consider the equation in complex variables ( z = x + iy), with F c a fixed polynomial that is well-chosen starting with second order terms, having a center at the origin. The Lyapunov quantities, L j , associated with equation where F is a polynomial perturbation, also starting with second order terms in z,z, can be written, using its Taylor series in the parameter space, as where L ( ) j are homogeneous polynomials of degree in the coefficients of F. All the computations in this paper are done adapting the algorithm described above, see also [19], for obtaining directly the Taylor series of L j up to order , generally = 1, and in some cases = 2 or higher, at the selected center.

Simultaneous Hopf bifurcations
The use of previous described center bifurcation method, combined with symmetry techniques (reversible systems), is a useful tool in obtaining lower bounds for the number of limit cycles in polynomial systems, see [7]. We use this tool, with two perturbations, to obtain limit cycles appearing simultaneously in three nests of limit cycles, two symmetric outside the symmetry line and one on it. We are going to present the idea of combining these two perturbations to generate limit cycles.
Let f and g be polynomials in x, y. Consider a polynomial reversible system with an equilibrium point of center type at the origin having two more symmetric equilibrium points, (x c , ±y c ), of center type outside the symmetry line. This configuration is depicted in figure 1. We consider a general perturbation of degree N of system (7) in two independent steps. We start studying the limit cycles bifurcating from the symmetric centers at (x c , ±y c ) but preserving the center at the origin. In order to use the degenerated Hopf bifurcation explained in section 3.1 we should move such center equilibrium point to the origin. Then, we take a perturbation of system (7) keeping the reversibility symmetry at the origin, In particular the origin remains as a critical point of center-focus type. In the perturbation above, we have considered also linear perturbations but with zero trace in order to have an extra parameter in the next step. Then, we translate system (8) in order that the center (x c , y c ) moves to the origin. Additionally, we should restrict the perturbation parameters f k, such that the perturbation terms in the translated system start with terms of degree 2. We will see in the examples that these conditions reduce in 5 the number of free parameters. Then, the method of small limit cycles bifurcating from the origin, described in the first part of this section, applies. Let n be the maximum number of limit cycles obtained around the origin with this technique. Then, we have n small limit cycles by a simultaneous degenerated Hopf bifurcation surrounding each one of the equilibrium points (x c , ±y c ).
We remark that the origin of system (8) remains a center when the 2n limit cycles have been created. Clearly, with this method we can obtain, at most, as limit cycles as the number of free parameters, that is n (N + 5)(N − 2)/2. This number is obtained counting the parameters in the perturbation in (8), N(N + 3)/2, removing the 5 commented above. In the next sections it is shown that in our best results we never reach this number but, nevertheless, we are close to it.
The second step is to consider a perturbation that completely breaks the symmetry at the origin, Let m be the maximum number of small limit cycles obtained using, once again, the previous described technique to study the degenerate Hopf bifurcation, but now applied to system (9). As above, we can obtain, at most, as many limit cycles as the number of free parameters; that is, m (N + 4)(N − 1)/2. Finally, from the above arguments, and since the perturbation parameters in systems (8) and (9) are independent, we have that the 2n small limit cycles surrounding the equilibria (x c , ±y c ) are hyperbolic and remain after the other perturbation, i.e. they remain while m limit cycles emerge from the origin. This method is used in detail, in the next section, for a cubic system.
Considering all together, at least 2n + m small limit cycles has been obtained; m small limit cycles surrounding the origin and n small limit cycles surrounding each one of the two symmetric equilibrium points. We will say that the 2 × n + m limit cycles are in configuration n, m, n .

Reversible cubic centers with two extra symmetric centers outside the symmetry line in a triangular configuration
The aim of this section is to present the method introduced in section 3, of simultaneous occurrence of limit cycles via a degenerate Hopf bifurcation, through a particular case. We introduce this method by studying some cubic reversible family of centers. We want to remark that using other families of cubic centers and other bifurcation techniques, more limit cycles can be obtained, see [17,18]. First, see proposition 3, we classify all reversible cubic systems having a non-degenerate equilibrium point of center type at the origin plus two more non-degenerate centers outside the reversibility symmetry line when the three centers are not aligned. Clearly, up to an affine change of variables and a time rescaling if necessary, it is not restrictive to assume that the symmetry axis is the straight line y = 0. Hence, the reversible cubic centers family that we consider is x = −y + a 11 xy + a 21 x 2 y + a 03 y 3 , Second, see proposition 5, following the procedure described in the previous section, we study the simultaneous degenerate Hopf bifurcation of limit cycles in configuration 1, 7, 1 for a specific family in (10).
As we will see in the proof we can not provide better examples using this bifurcation technique for cubic vector fields. We recall that the family of cubic systems can exhibit up to 13 limit cycles, as it is proved in [17,18].
We finish this section showing some systems of family (10), with different phase portraits, exhibiting the same configuration 1, 7, 1 of limit cycles. All having the three centers in the vertices of a triangle. This is why we have added the condition x c y c = 0 in the next result.

Proposition 3.
After a rescaling, if necessary, the reversible cubic system (10) has three non-degenerate centers in a symmetric triangular configuration (one at the origin and two at (x c , ±y c ) with x c y c = 0) if and only if, the determinant of the Jacobian matrix of system (10) at (x c , y c ), α, is positive, a 03 = 1, and c 1 = 0 or c 2 = 0, with Moreover, system (10) either has a polynomial inverse integrating factor when c 1 = 0 or a polynomial first integral of degree four when c 2 = 0.
Proof. Straightforward computations show that when a 03 = 0 it is not possible to have a weak focus at (x c , y c ). Hence, after a re-scaling of the (x, y) variables if necessary, we have a 03 = 1. Then, three of the conditions (11) follow imposing that system (10) has a non-degenerate weak focus at (x c , y c ). That is vanishing the two components of the vector field at (x c , y c ) and imposing that the trace is also zero. The fourth appears because we have added an extra parameter, that we call α, to simplify the computations. It is the determinant of the Jacobian matrix at such equilibrium point. Additionally, it should be positive. Following the procedure described in section 3, the first Lyapunov quantity is where c 1 and c 2 are given by expressions (12) and (13), respectively. Hence, conditions c 1 = 0 or c 2 = 0 are necessary for the existence of a center at (x c , y c ). Straightforward computations show that the next two Lyapunov quantities are zero.
The proof finishes showing that the above necessary conditions are also sufficient. Direct calculations show that, when c 2 = 0, system (10) is a Hamiltonian system with a degree four polynomial first integral. The case c 1 = 0 is more intricate. Generically the associated system has no invariant straight lines but it has two invariant quadratic curves, that in general are complex. From them, it can be obtained the real polynomial inverse integrating factor where □ Remark 4. In the above proof it is explained that when the first Lyapunov quantity is zero the next two also vanish. This fact could indicate that only one limit cycle can be obtained around the equilibrium point lying outside the symmetry axis.
In the proof of the next result it can be seen the two invariant quadratic curves that provides the inverse integrating factor (14).

Proposition 5. The cubic system
has a center at the origin and two symmetric centers at (1, ±2). Moreover, there exists a cubic polynomial perturbation such that it unfolds 9 small limit cycles by a simultaneous degenerate Hopf bifurcation in configuration 1, 7, 1 . In detail, after perturbing, 7 limit cycles surround the origin while each one of the 2 other limit cycles surrounds one of the perturbed symmetric centers.
Proof. The first statement follows from proposition 3 checking that system (15) satisfies c 1 = 0 in (12) taking (x c , ±y c ) = (1, ±2) and α = 1176 2 . Moreover, the polynomial inverse integrating factor (14) factorizes as Hence, we can compute the rational first integral The second statement follows from considering the perturbation monomials in two steps, first the reversible perturbation monomials with respect to the origin, and then the non reversible ones. We can proceed in this way because both perturbations are independent. Let us expose these arguments in detail.
First, we perturb system (15) only with the cubic reversible perturbation monomials as it is given in expression (8), that is (17) We perturb in this way because this perturbation does not break the symmetry at the origin and, hence, we keep the center at that point. Adding the next conditions, 12 , the points (1, ±2) remain as equilibrium points of system (17) with zero trace. That is, they are either centers or weak foci. Now, we will study the degenerate Hopf bifurcation and the maximum order of the equilibria (1, ±2) as weak foci.
Following the notation introduced in (6), the linearization of the first Lyapunov quantity write as L (1) with y 1 = (244f 12 − 61f 30 + 488f 03 )/153 664. Clearly, using this independent quantity, and after changing the trace at (1, ±2), we can get one small limit cycle surrounding each equilibrium by using the double symmetric Hopf bifurcation. Moreover, the origin remains as a center. As the linearization of the next two Lyapunov quantities write as the use of only first order terms does not provide more limit cycles.
Computing higher order terms, in the Taylor series expansion of the Lyapunov quantities with respect to the perturbation coefficients, we can get no more limit cycles surrounding these equilibrium points. An alternative way to see this for our family is using the proof of proposition 3. It has been shown that, when the trace and L 1 vanish we have already a center.
From the above first order study, using the implicit function theorem, up to a change of variables (in the parameter space) if necessary, we can write Finally, as the above homogeneous parts of degree three are also multiple one from the other and since L (4) 7 = 0, we need to compute up to order five the first seven Lyapunov quantities. If we proceed in this way, doing the same simplification procedure as in the previous steps, we get Then, after some necessary arrangements because of the linear dependence of the first order terms, the Lyapunov quantities up to order five can be written as The polynomial p 2 has two simple real zeros, λ 1 and λ 2 , and the resultant of p 2 and p 4 is different from zero. Hence, for every x 6 small enough and different from zero, there exist two curves near x 1 = x 2 = · · · = x 5 = 0 and x 7 = λ x 6 , = 1, 2 such that L j = 0, j = 1, . . . , 6 and L 7 = 0. Moreover, over these curves, the intersection of the varieties L j is transversal. Consequently, as in [6], they define two families of weak foci of order seven that unfold, perturbing also the parameter that defines the trace at the origin, seven small limit cycles via a degenerate Hopf bifurcation. Hence, by computing the Lyapunov quantities up to order five in terms of the perturbation parameters, we have proved the existence of cubic polynomial perturbations of system (15) unfolding 9 small limit cycles by a simultaneous degenerate Hopf bifurcation in configuration 1, 7, 1 , as we wanted prove. Finally, we would like to remark that the homogeneous terms of degree 2 and 4 for all the computed Lyapunov quantities corresponding to the origin of system (18) vanish identically. This explains why, in the above computations, there only appear the terms of degrees 1, 3, and 5. Additionally, we have not written down all the digits of all the involved integer numbers in the previous expressions because of their size. □ Under the assumptions of proposition 3 there are other cubic systems with the same configuration of two symmetric centers and one more on the symmetry line and exhibiting the same bifurcation phenomenon. All of them with 9 limit cycles, but having different type of Darboux first integrals and global phase portraits. For example: The corresponding phase portraits are depicted in figure 2. The finite equilibrium points are all of them of center or saddle type, meanwhile at infinity the equilibrium points, when they exist, all are of node or saddle type. Finally, we remark that system correspoding to (22) belongs to the family studied in [17] where 13 limit cycles appear bifurcating simultaneously from different period annulus.

Reversible quartic systems
Christopher, in [6], study the simultaneous bifurcation of limit cycles in configuration n, m, n for a quartic Darboux center given by Żołądek in [31]. Using the method described in section 3, up to first order, i.e. calculating the linear parts of the Lyapunov quantities, he provides 22 limit cycles in configuration 6,10,6 .
In this section we extend this study to some reversible quartic systems having rational first integrals, using linear and higher order developments of the Lyapunov quantities. More concretely, let us consider rational first integrals of the type The given polynomials p N,d−1 and q N,d are of degree d − 1 and d, respectively. The corresponding differential system, that it has degree N, writes as for a suitable inverse integrating factor V N,d (x, y). All the reversible systems presented in this section are symmetric with respect to the x-axis having a center at the origin and two more at (x c , ±y c ). All of them are non-degenerate.
In the first subsection, using the linearization of the first Lyapunov quantities and simultaneous bifurcations, we obtain different configurations with 16, 19, and 22 limit cycles. Clearly, the total number increase with d. The last one, H 4,5 , is in the same class as Christopher's center and we recover his result. In the second subsection, using fifth order developments, we get the new lower bound H(4) 28, in a simultaneous bifurcation and in a 8,12,8 configuration. See proposition 6. As in the cubic family studied in the previous section, only the developments of odd order are useful for obtaining more limit cycles.

Configurations found in the class H 4,d
Next table shows the results on simultaneous bifurcation of limit cycles for system (24), of degree four, and d = 3, 4, 5, using only the linear parts of the Lyapunov quantities. We indicate the found configuration and the total number of small limit cycles. Because of the size of the expressions of the Lyapunov quantities and as we are following exactly the procedure described in section 3 we only list the polynomials that define the rational first integrals (23):

Higher order studies for a quartic system
In proposition 5 we have seen how higher order developments of the Lyapunov quantities can be used to increase the number of small limit cycles. In that case only from the center on the symmetry line we can get more limit cycles, passing from the configuration 1, 5, 1 , with a first order study, to the configuration 1, 7, 1 , with a fifth order development.
In this subsection, we get the new lower bound of H(4) 28 by proceeding in a manner similar to the previously exposed. Due to the difficulties on the computations we have only checked how the number of limit cycles increases for a system of degree four which is of type H 4,5 . As we have mentioned above, we have selected a Darboux center which is in the same class as the one studied by Christopher. Proposition 6. The quartic differential system associated to the rational first integral has one reversible center at the origin and two more centers at the points (1, ±2). Perturbing with polynomials of degree four, from these centers bifurcate 22, 27, and 28 limit cycles in configurations 6, 10, 6 , 8, 11, 8 , and 8,12,8 , using the series expansion of the Lyapunov quantities on the perturbation coefficients up to order 1, 3, and 5, respectively.
Proof. It is easy to check that the reversible differential equation associated to the rational first integral (25) has a center at (0, 0) on the symmetry axis and two more symmetric centers at (1, ±2). The limit cycles appear via a simultaneous degenerate Hopf bifurcation, following the procedure described in section 3. More concretely, let us denote by L ( ) j and L ( ) j the Taylor series expansion of order of the j Lyapunov quantities at the equilibrium points obtained perturbing the differential equation, i.e. near the centers at (0, 0) and (1, ±2), respectively. We notice that, by the symmetric perturbation procedure, the j Lyapunov quantities at (1,2) and (1, −2) coincide.
We will follow the same scheme described in the proof of proposition 5. Straightforward computations provide, after some linear changes in the parameter space, the first Lyapunov quantities up to order one. We can write them as L (1) j = x j + O 2 (x) for j = 1, . . . , 10 and L (1) j = y j + O 2 (y) for j = 1, . . . , 6. We have denoted by O 2 (x) and O 2 ( y ) all the terms of degree two in variables x j and y j , respectively. Moreover, we observe that L (1) 11 , L (1) 12 , L (1) 7 , and L (1) 8 are linear combinations of the first ten and six Lyapunov quantities up to order one, respectively. From these expressions, adding the perturbation parameters that control the traces at the critical points after the perturbation, we get the same configuration of limit cycles, 6, 10, 6 , given by Christopher in [6]. For obtaining more limit cycles, before using the 'trace parameters', we need to compute the Taylor series expansion of higher order of the next Lyapunov quantities L j for j = 11, 12 and L j for j = 7, 8.
The second order terms in the Taylor developments of the Lyapunov quantities are zero when the parameters that control first order developments vanish, i.e. when x 1 = · · · = x 10 = y 1 = · · · = y 6 = 0. Then, third order Taylor developments of the Lyapunov quantities are necessary to be computed so that more limit cycles could be obtained. If we compute them, we obtain the expressions L (3) for some constants C 1 , C 2 , C 3 , and C 4 and polynomials We want to remark that we have four functions but only three polynomials. From previous expressions, we only obtain, up to order three, a final extra limit cycle surrounding the origin. This is so because L (3) 11 is different from zero and the polynomial p 2 appears in the two expressions L (3) 11 and L (3) 12 . Surrounding each one of the critical points (1, ±2) we have two more limit cycles, because the polynomials p 3 and p 3 are of degree three and they have no common zeros, as it follows from the fact that the resultant of them is not zero. Consequently, up to order three, we have the configuration 8, 11, 8 of limit cycles.
We have not written completely the above polynomial of degree four, because each coefficient is an integer number with more than 200 digits. Moreover, the resultant of the polynomials p 3 and p 4 is non zero, then the existence of an extra hyperbolic limit cycle is guaranteed. Consequently, we get at least twenty-eight small limit cycles in configuration 8, 12, 8 . This finishes the proof. □ We have not computed higher order Lyapunov quantities to find more limit cycles because to get the maximum, 28 limit cycles, we have exhausted all the perturbation parameters. To get more limit cycles, if they exist, other mechanisms should be applied.

Simultaneous bifurcation of higher degree systems
This section is devoted to the perturbation of vector field (24) for degrees N = 5, . . . , 9. Here we provide the new best lower bounds H(5) 37 and H(6) 53, among others lower bounds. To reach them we use linearizations of Lyapunov quantities, following the way introduced in the previous sections. Next proposition proves the lower bounds given in theorem 1 for N = 5 and N = 6.
For the next values of N, similar computations, as the ones we have performed for degree 4 in the previous section, can be done. It can be seen that the sequence of the number of limit cycles in each configuration is increasing up to the maximal value obtained for d = 2N − 3. Higher values of d do not improve the results obtained for such specific degree. Moreover, when we increase the degree N we get better results in next sections using other type of systems. These, together with the computational difficulties when N increases, are the main reasons why we stop the computations at these values.
All the systems are chosen so that the involved rational numbers be as small as possible. We observe that the method described in section 3 needs a normal form transformation where the square root of the determinant of the Jacobian matrix at the equilibrium points appears in the computations. In most of the cases we achieve that the determinant of the Jacobian matrix at the equilibrium points be a perfect square. This fact simplifies significantly the involved calculations. Alternatively, when a square root of an integer number appears we look for a system that has this number as small as possible. Proof. The proof follows as the proof of proposition 6. But only computing the linearization of the Lyapunov quantities at the origin and at the other two symmetric centers. For the sake of brevity, we omit the computations because of their size. The polynomials defining the corresponding rational first integral (23) are: Proof. The proof follows using the same procedure described in proposition 6 but studying, as in proposition 7, only the linearization of the Lyapunov quantities. We only indicate the polynomial system having the three centers in a triangular symmetric configuration with respect to the x-axis such that the 2n + m simultaneous small limit cycles will bifurcate.
The key point is to consider all systems in the proof of proposition 7 adding a polynomial curve of equilibria avoiding the center points. This curve will increase the degree of the vector field and, in a natural way, also the number of limit cycles. Hence, choosing F M (x, y) as a polynomial of degree M, the polynomial system In all cases the centers remain at (0, 0) and (1, ±2) because the curves F 2 and F 3 do not vanish at such points. □ Straightforward computations show that changing the polynomials F M the configurations and the total number of limit cycles are the same.
We remark that the results of the above proposition are the best that we have obtained choosing different vector fields and also different polynomials F M . The best strategy has not been to select only one vector field and adding curves F M increasing the degree. For example, for H 5,7 and M = 1, 2, 3, 4 we obtain polynomial vector fields with degrees 6,7,8,9 having at least 51, 74, 93, 113 limit cycles, respectively. Their configurations are 15, 21, 15 , 23, 28, 23 , 29, 35, 29 , 36, 41, 36 . In this case only for N = 7 we have obtained the highest value. In fact, the values obtained for N = 6, 8, 9 with other Darboux first integrals and other polynomials F M are higher. In any case we think that, for every degree N, highest values of the lower bounds for the Hilbert numbers will be obtained perturbing good Darboux systems, instead of choosing other Darboux systems with lower degree multiplied by a line of equilibria. This is the case for degrees 5 and 6.

Lower bounds for the Hilbert numbers using double symmetry
This section is devoted to prove theorem 2. Its proof follows part of the approach introduced in [7] to get an asymptotic lower bound of the form H(N) > C N 2 log N, for some constant C. Even though there are lower bounds that asymptotically behave like this one, and that we do not improve, they do not give better results for small values of N, see for instance [13].
Let us consider a polynomial degree N differential system, (x , y ) = ( f (x, y), g(x, y)), with