This is a preprint of: “Lower bounds for the local cyclicity of centers using high order developments and parallelization”, Luiz Fernando da Silva Gouveia, Joan Torregrosa, J. Differential Equations, vol. 271, 447–479, 2021. DOI: [10.1016/j.jde.2020.08.027] LOWER BOUNDS FOR THE LOCAL CYCLICITY OF CENTERS USING HIGH ORDER DEVELOPMENTS AND PARALLELIZATION LUIZ F. S. GOUVEIA AND JOAN TORREGROSA Abstract. We are interested in small-amplitude isolated periodic orbits, so-called limit cycles, surrounding only one equilibrium point, that we locate at the origin. We develop a parallelization technique to study higher order developments, with respect to the pa- rameters, of the return map near the origin. This technique is useful to study lower bounds for the local cyclicity of centers. We denote by M (n) the maximum number of limit cycles bifurcating from the origin via a degenerate Hopf bifurcation for a polyno- mial vector field of degree n. We get lower bounds for the local cyclicity of some known cubic centers and we prove that M (4) ≥ 20, M (5) ≥ 33, M (7) ≥ 61, M (8) ≥ 76, and M (9) ≥ 88. 1. Introduction Hilbert early last century presented a list of problems that almost all of them are solved. One problem that still is open is the second part of the 16th Hilbert’s problem: It consists in determine the maximal number H(n) of limit cycles, and their relative positions, of planar polynomial vector fields of degree n. In last years other related problems have been proposed. In 1977, Arnol’d ([2]) suggested a weakened version, focused on the study of the number of limit cycles bifurcating from the period annulus of Hamiltonian systems perturbing with polynomials of fixed degree n. We are interested here in another local version, that consists in providing the maximum number M (n) of small-amplitude limit cycles bifurcating from an elementary equilibrium point of center or focus type, clearly M (n) ≤ H(n). In other words, M (n) is an upper bound of the cyclicity of such equilibrium points. More details on related problems about the cyclicity of homoclinic and heteroclinic connections or of period annuli can be found in, for example, Roussarie’s book ([28]). There are very few general families for which this local number M (n) is completely determined. In 1954 Bautin ([5]) proved that M (2) = 3 and Sibirskiı̆ ([29]) and Blows and Lloyd ([6]) studied the special cubic family without quadratic terms. But the first complete proof that Mh (3) = 5 was obtained by Żola̧dek ([33]) in 1994. The first evidence that M (3) ≥ 11 was presented also by Żola̧dek ([34]) in 1995, providing a cubic center with very high local cyclicity. This problem was revisited by himself in 2016 ([36]). The first and simplest proof that this lower bound is reached for the cubic family was given by Christopher in 2005 ([12]). Although we will describe more precisely this mechanism later, Christopher’s idea is based in choosing a point on the center variety and, at this point, study the independence property of the first-order Taylor series with respect to (1) the perturbation parameters, denoted by Lk , of the Lyapunov constants. In fact, if the point is (generically) chosen on a component of the center variety of codimension r, then there exist perturbations which can produce r − 1 limit cycles. Apart from the fact that the center problem for polynomial vector fields of degree n is unsolved, the main problem is how to compute the codimension of each component of the center variety 2010 Mathematics Subject Classification. Primary 34C23, 34C25, 34C07. Key words and phrases. small-amplitude limit cycle, polynomial vector field, center cyclicity, Lyapunov constants, higher-order developments and parallelization. 1 2 LUIZ F. S. GOUVEIA AND JOAN TORREGROSA or, alternatively, how to find good center candidates to be perturbed for obtaining the highest (local) cyclicity value. This last procedure is the one mostly used to provide the best (highest) lower bounds for M (n). It is somewhat curious that both centers appear listed in the same work of Żola̧dek ([35]). We notice that the idea of studying only linear developments, with respect to the parameters, near centers appears in some previous works of Chicone and Jacobs ([10]) and of Han ([23]). In 2012, Giné ([19, 20]) conjectured that M (n) = n2 + 3n − 7. This suggests a very high value for M (n) for polynomial vector fields of low degree. In particular, for n = 4, it asserts that M (4) ≥ 21. The local cyclicity problem in the quartic polynomial class was studied in [19] using only second-order Taylor developments of the Lyapunov constants. But, as we will describe in Section 5, orders 2 and 3 are not enough to prove this result. Then, this lower bound remains to be proved. For degrees n = 5, 7, 8, 9 the best lower bounds for M (n) were obtained in [25] providing examples exhibiting 28, 54, 70, and 88 limit cycles, respectively. Next result improves all these values increasing the known lower bounds for M (n) for the respective degrees. Up to our knowledge, the highest value for n = 6 is M (6) ≥ 48 and it is proved in [4]. This suggests that the value for M (n) should be increased (at least) by one. The new conjectured value appears recently in [21]. A final comment about the Giné’s conjecture for n = 5 is that the value M (5) ≥ 33 suggested is reached in our main result. The proof will follow perturbing a quintic center with homogeneous nonlinearities, appearing in [18], in the complete polynomial class of degree 5. The best result for general polynomial vector fields of degree n is M (n) ≥ n2 − 2. A result of 2004 due to Movasati ([27]). Theorem 1.1. The number of limit cycles bifurcating from a singular monodromic point for vector fields of degree four, five, seven, eight and nine is at least M (4) ≥ 20, M (5) ≥ 33, M (7) ≥ 61, M (8) ≥ 76, and M (9) ≥ 88. A better description of the number of limit cycles in the cubic polynomial family can be obtained studying lower bounds for the local cyclicity of some cubic Darboux centers. We analyze here all those with codimension 12 and listed by Żola̧dek in [35] but having a real center. More concretely, we show that in almost all cases always 11 limit cycles bifurcate from each studied center. The proofs show why sometimes a very high-order study is necessary. For proving our results first we extend, in Section 3, the parallelization technique intro- duced by Liang and Torregrosa in [25] to higher-order developments. Second, following the ideas in [12] for studying higher-order developments, we use a blow-up procedure to get a complete unfolding of the return map near a polynomial center perturbing with polynomials of the same degree. We remark that the parallelization technique drastically reduces the computation time. In particular, in the proof of Theorem 1.1 when only developments of order two are necessary to be used, instead of more than one month of computation time, we need less than one hour. We have used a cluster of computers with ninety processors simultaneously. All the computations have been made with the Computer Algebra System Maple [26]. This paper is structured as follows. In Section 2, we recall the necessary definitions and algorithms to get the coefficients of the return map, the so-called Lyapunov constants among other preliminary results as the Poincaré–Miranda Theorem (see [24]) and the Gershgorin Theorem (see [17]) about localization of eigenvalues of a matrix. In Section 3, we present and prove the parallelization setting. The study of the local cyclicity of cubic Darboux centers is developed in Section 4. Finally, in Sections 5 and 6 we prove our main result Theorem 1.1. All the computations are analytic but, in some cases, we have used LOWER BOUNDS FOR THE LOCAL CYCLICITY OF CENTERS 3 numerical approximations that, with Computer Assisted Proofs techniques, have allowed us to obtain the analytical proof. 2. Degenerated Hopf bifurcation In this section, we recall how to obtain the Lyapunov constants or focal values, that is, the coefficients of the return map near an equilibrium point with differential matrix of an elementary center type. As usual, it is not restrictive to assume that the equilibrium point is located at the origin. So, following [1], we can consider the polynomial system of degree n  Xn     ẋ = λ x − y + Pk (x, y),  0 k=2 Xn (1)     ẏ = x + λ0 y +  Qk (x, y), k=2 being Pk and Qk homogeneous polynomials of degree k in the variables x, y. Writing the above system in polar coordinates, (r, θ) = (r cos θ, r sin θ), we get X∞ dr = λ0 r + Sk (θ)rk , (2) dθ k=2 being Sk (θ) trigonometric polynomials. Let r(θ, ρ) be the solution of system (2) such that r(0, ρ) = ρ. The stability of the origin is clearly determined, using Hartman–Grobman Theorem, when λ0 6= 0. When λ0 = 0 the stability problem is known as the center-focus problem and there are some classical tools to distinguish when the origin is stable or unstable. So, for λ0 = 0 and close to ρ = 0, we can develop this solution as a power series in ρ X∞ r(θ, ρ) = ρ + rk (θ)ρk , k=2 where rk (0) = 0 for all k ≥ 2. Then, the Poincaré return map, Π(ρ), can be obtained evaluating the above expression at θ = 2π, i.e. ∞ X Π(ρ) = r(2π, ρ) = ρ + rk (2π)ρk . (3) k=2 When VK̃ = rK̃ (2π) 6= 0 for some value K̃ we say that the origin of system (1) is a weak-focus, otherwise we say that the origin is a center. In this context, it is well known that the first nonzero value, when it exists, corresponds to an odd subscript K̃ = 2K + 1, see [1, 9, 16], and consequently the K-Lyapunov constant is defined as LK = V2K+1 when L1 = · · · = LK−1 = 0. Then, we say that the origin is a weak-focus of order K when LK 6= 0 and L1 = · · · = LK−1 = 0. These constants are polynomials in the coefficients of Pk and Qk defined in (1). For more details in the algebraic properties that they satisfy, we refer the reader to [13]. When the first Lyapunov constant is negative (positive), L1 6= 0, a small-amplitude stable (unstable) limit cycle bifurcates from the equilibrium when the trace parameter λ0 moves from zero to positive (negative). Because the stability of the equilibrium point changes and the limit cycle appears by using the Poincaré–Bendixson Theorem. This phenomenon is known as the classical Hopf bifurcation. The results in this work are dealing with the K-degenerated Hopf bifurcation. That is, when K small-amplitude limit cycles bifurcate from a weak-focus of order K, for more details see for example [11]. According Roussarie in [28], at most K limit cycles can bifurcate from a weak-focus of 4 LUIZ F. S. GOUVEIA AND JOAN TORREGROSA order K under analytic perturbations. In the context of (1), the main difficulty is how can we ensure the existence of polynomial perturbations such that the K limit cycles bifurcate from a weak-focus of order K. To unify notation, we denote the trace parameter in (1) by L0 := λ0 . The fact that the coefficients of the return map (3) corresponding to monomials of even degree in ρ do not play any role in the bifurcation phenomenon is due to the property that the Bautin ideal is generated only by the coefficients of odd degree, i.e. B = hL1 , . . . , Ln , . . .i = hV3 , V4 , . . . , V2n , V2n+1 , . . .i. See more details also in [28]. This property has been revisited recently in [14]. Another equivalent procedure to study the center-focus problem is to propose, for λ0 = 0, a function H(x, y) = x2 + y 2 + O(||(x, y)||3 ) such that, using (1), ∂H ∂H Ḣ = ẋ + ẏ = h4 r4 + h6 r6 + · · · + h2k r2k + · · · , ∂x ∂y with r2 = x2 + y 2 . Then, the first nonvanishing coefficient of the above derivative, which always has an even subscript, determines the stability of the origin of (1), being H a Lyapunov function. In fact, both coefficients, h2K+2 and L2K+1 , differ on a multiplicative nonzero constant when the previous vanish. Here the center property reads as h2k = 0 for all k. In complex variables, via the change of variables z = x + i y and for λ0 = 0, system (1) writes as n X ż = R(z, z̄) = i z + Rk (z, z̄), k=2 where Rk (z, z̄) are homogeneous polynomials of degree k in (z, z̄). Consequently, the above function H = z z̄ + O3 (z, z̄) satisfies X∞ ∂H ∂H Ḣ = ż + ˙ z̄ = gk (z z̄)k , ∂z ∂ z̄ k=2 being the coefficient gk+1 the k-Lyapunov constant. Following the Christopher’s idea ([12]), to find good (or high) lower bounds for M (n), instead of studying the local cyclicity of the origin of system (1), we consider perturbations of a fixed center. That is, taking Pc , Qc , P, and Q as polynomials of degree n in (x, y) and studying the perturbed system  ẋ = Pc (x, y) + P (x, y, λ), (4) ẏ = Qc (x, y) + Q(x, y, λ), with P, Q having only monomials of degree higher or equal than two. When λ = 0, P (x, y, 0) = Q(x, y, 0) = 0, and the unperturbed system must be written as (1), i.e. for λ0 = 0, the origin should be an elementary center. Thus, we are interested in finding limit cycles of small-amplitude bifurcating from the origin of system (4). In most cases, the explicit computation of the Lyapunov constants for system (1) is very hard and impossible to do by hand. So a Computer Algebra System1 is necessary to be used. Moreover, sometimes more specific algorithms can be developed to decrease the computation time. Here, with the aim to reduce not only the computation time but also the total memory requirements, we extend, in the next section, the parallelization algorithm proposed by Liang and Torregrosa in [25] to higher-order developments but for systems of type (4). Next two results, which are proved in [12], provide conditions to get the bifurcation of small-amplitude limit cycles near a polynomial center. The first uses only the linear developments of the Lyapunov constants. We notice that as we are perturbing centers, 1We have used Maple for all our computations ([26]) LOWER BOUNDS FOR THE LOCAL CYCLICITY OF CENTERS 5 the Taylor approximations vanish at the origin (λ = 0) in the parameter space. So, the first-order Taylor approximation is linear. As we have commented before, similar versions of this result can be found in [10] and [23]. The second result uses higher-order Taylor approximations. Theorem 2.1 ([12]). Suppose that s is a point on the center variety and that the first k Lyapunov constants, L1 , . . . , Lk , have independent linear parts (with respect to the ex- pansion of Li about s), then s lies on a component of the center variety of codimension at least k and there are bifurcations which produce k limit cycles locally from the center corresponding to the parameter value s. If, furthermore, we know that s lies on a compo- nent of the center variety of codimension k, then s is a smooth point of the variety, and the cyclicity of the center for the parameter value s is exactly k. In the latter case, k is also the cyclicity of a generic point on this component of the center variety. The scheme of the proof is as follows. Under the hypotheses of the above result, there exists a change of variables such that the first Lyapunov constants write as Li = ui + O2 (u1 , . . . , uk , . . . , um ), i = 1, . . . , k, (5) assuming that we have m ≥ k bifurcation parameters. Using the Implicit Function Theorem it is clear that we can write Li = vi , for i = 1, . . . , k. Then the first coefficients of the return map (3) are independent. It is clear that, when the trace parameter L0 = λ0 = 0 in (1) we get only k − 1 limit cycles. But adding the parameter L0 = λ0 we have an extra limit cycle by the classical Hopf bifurcation, obtaining in total k, as the above result ensures. In fact, this proves the existence of a variety, in the parameter space, of weak- foci of order k ({λ0 = u1 = u2 = · · · = uk−1 = 0} and uk 6= 0) that unfolds k hyperbolic limit cycles. Using the Weierstrass Preparation Theorem this is the maximal number near such curve. When the linear parts of the next Lyapunov constants are linear combination of the first k we can use the higher developments to obtain more limit cycles. This is the aim of the next result also proved by Christopher in [12, Theorem 3.1]. Theorem 2.2. Suppose that we are in a point s where Theorem 2.1 applies. After a change of variables if necessary, we can assume that L0 = L1 = · · · = Lk = 0 and the next Lyapunov constants Li = hi (u) + Om+1 (u), for i = k + 1, . . . , k + l, where hi are homogeneous polynomials of degree m ≥ 2 and u = (uk+1 , . . . , uk+l ). If there exists a line `, in the parameter space, such that hi (`) = 0, i = k + 1, . . . , k + l − 1, the hypersurfaces hi = 0 intersect transversally along ` for i = k + 1, . . . , k + l − 1, and hk+l (`) 6= 0, then there are perturbations of the center which produce k + l limit cycles. The above result is not written exactly as in the original Christopher paper because we have adapted to include also the conclusion of Theorem 2.1. We have also included here an alternative proof. Proof of Theorem 2.2. The first step uses Theorem 2.1 to assume that the first k Lyapunov constants together with the trace are the first k + 1 coordinates (u0 , u1 , . . . , uk ), that will be taken as zero. The proof continues taking the blow-up change of variable uj = vj uk+l , for j = k+1, . . . , k+l−1. Then v = (vk+1 , . . . , vk+l−1 ) and we can write hi (u) = um b k+l hi (v), for i = k + 1, . . . , k + l. Consequently, ∞ ! X Li (u) = um e m b j k+l Li (v) = uk+l hi (v) + gij (v)uk+l . j=1 The existence of a line ` as in the statement gets v ∗ such that b hi (v ∗ ) = 0, for i = k + 1, . . . , k + l − 1, b hk+l (v ∗ ) 6= 0, and the determinant of the Jacobian matrix of 6 LUIZ F. S. GOUVEIA AND JOAN TORREGROSA (b hk+1 , . . . , b hk+l−1 ) with respect to v does not vanish at v ∗ , then the Implicit Function Theorem applies, in a small neighborhood of v ∗ , and the change of variables wi = L ei (v) is well defined. The proof follows the same scheme explained in the comments before to state this result but changing uj = wj um m k+l . Because, now uk+l is a common factor of the complete return and the polynomial provided by the Weierstrass Preparation Theorem has b hk+l (v ∗ ) 6= 0 as the coefficient of maximal degree monomial and the other are the independent coefficients wj , for j = 1, . . . , k + l − 1. Finally, as above we can use the trace parameter to have a complete unfolding of k + l limit cycles.  We remark that, as in the previous result, this maximal number of limit cycles is obtained only near the weak-foci curve. Consequently, the previous results provide only lower bounds for the local cyclicity problem. Theorem 2.1 uses first order developments. After removing them (using the Implicit Function Theorem) Theorem 2.2 imposes that all the homogeneous principal parts at the origin of the remaining hi have the same degree. Our proof of Theorem 2.2 suggests that we can study, restricting the parameter space if necessary, how the intersection of the algebraic varieties SL = {L1 (u) = L2 (u) = · · · = Lm (u) = 0} is near u = 0. This is equivalent to know if the function f : (Rn , 0) → (Rn , 0) (u1 , . . . , un ) 7→ (L1 , L2 , . . . , Ln ) is locally surjective at the origin (see the interesting comment on this fact in [28, Page 71]). From the singularities classification theory, see [3], some properties used to compute the local multiplicity of a function in a point, µ0 [f ], are useful to study the local intersection of system SL . Proposition 2.3 ([3]). Let f : (Rn , 0) → (Rn , 0) be a finite map germ. Defining f = (f1 , . . . , fn ) and Sf = {f1 = f2 = · · · = fn = 0} the next properties hold, Q (i) if fi = hiQ+ H.O.T., where hi is homogeneous of degree ki , then µ0 [f ] ≥ ni=1 ki and µ0 [f ] = ni=1 ki if and only if the system hi = 0, i = 1, . . . , n, has only the trivial solution in C P. n (ii) if gi = fi + j 0 for each 1 ≤ i ≤ n, where Si± = {(x1 , . . . , xn ) ∈ S : xi = ±c}. So, there exists d ∈ S such that f (d) = 0. Theorem 2.5 ([17], Circles of Gershgorin). Let A = (ai,j ) ∈ Cn×n and αk its eigenvalues. Consider for each i = 1, . . . , n Di = {z ∈ C : |z − ai,i | ≤ ri }, X where ri = |ai,j |. So, for all k, each αk ∈ Di for some i. i6=j The Poincaré–Miranda’s Theorem was conjectured by Poincaré in the 19th century and proved by Miranda in last century. Note that this result is a generalization of the Bolzano’s Theorem for higher dimensions. The reader can get more details on Gershgorin Circle Theorem in [22]. 3. Parallelization In [25], Liang and Torregrosa present a parallelization mechanism to compute the Taylor linear approximation of the Lyapunov constants near centers. In this section, we extend this result to compute also the terms of higher degree. We start recalling the linearization result for completeness. Theorem 3.1 ([25]). Let p(z, z̄) and Qj (z, z̄), j = 1, . . . , s be polynomials with monomials (1) of degree higher or equal than two such that the origin of ż = i z+p(z, z̄) is a center. If Lk,j denotes the first-order Taylor approximation, with respect to λj ∈ R, of the k-Lyapunov constant of equation ż = i z + p(z, z̄) + λj Qj (z, z̄), j = 1, . . . , s, LOWER BOUNDS FOR THE LOCAL CYCLICITY OF CENTERS 9 then the first-order Taylor approximation of the k-Lyapunov constant, with respect to Λ = (λ1 , . . . , λN ) ∈ RN , of equation N X ż = i z + p(z, z̄) + λj Qj (z, z̄), j=1 (1) P N (1) is Lk = Lk,j . j=1 In the following Theorem 3.3, we show how to get the terms up to degree ` > 1 of the Lyapunov constants. The main idea is to decompose the global computation problem into a collection of simpler problems. The advantage of the previous result is that each perturbation parameter λj appears in only one problem. We remark that the authors were interested only in the linear approximation. Now, because we study developments up to degree `, we need to decompose the global problem in simpler problems having exactly ` parameters or monomials. But, as many parameters appear in more than one simple perturbation problem we need to correct (multiplying by a factor adequately chosen) the obtained coefficients of the developments of degree `. We will see that the idea of parallelization only makes sense if the number of parameters is greater than the degree up to which we want to compute. To further clarify the main idea of parallelization, we present the next proposition where we perturb a quadratic center with four parameters and we calculate the first three Lyapunov constants up to degree two. Proposition 3.2. Consider the quadratic perturbed equation ż = i z + z 2 + (a20 + i b20 )z 2 + (a11 + i b11 )z z̄. (9) Then the Taylor developments up to degree two of the Lyapunov constants of the above system are (2) L1 = −2a20 b11 − 2b20 a11 − 2b11 , (2) L2 = 36a20 b11 + 12b20 a11 + 32a11 b11 + 8b211 + 12b11 , (2) L3 = −540a20 b11 − 108b20 a11 − 582a11 b11 − 192b211 − 108b11 . Proof. We rename the four real parameters in the ordered list λ = (λ1 , λ2 , λ3 , λ4 ) = (a20 , b20 , a11 , b11 ). As system (9) is a holomorphic quadratic center when λ = 0, all the Lyapunov constants vanish. Then, as we are interested in degree two developments and we have 4 parameters, we will consider the six combinatorial (without repetition) pairs S = {(λ1 , λ2 ), (λ1 , λ3 ), (λ1 , λ4 ), (λ2 , λ3 ), (λ2 , λ4 ), (λ3 , λ4 )}, and 6 respective differential equations, one for each element in S, such that in (9) the parameters that are not in the chosen pair are taken to be zero. For example, the corre- sponding differential equation to S3 = (λ1 , λ4 ) is equation E3 : ż = i z + z 2 + λ1 z 2 + i λ4 z z̄. For each equation Ej corresponding to each pair Sj we compute, with the mechanism described in Section 2, the first three Lyapunov constants. We denote by Lk,j the k- (2) Lyapunov constant of equation Ej and by Lk,j the corresponding Taylor approximations (2) up to degree 2. Straightforward computations show Lk,j = 0, for k = 1, 2, 3 and j = 1, 2 10 LUIZ F. S. GOUVEIA AND JOAN TORREGROSA and (2) (2) L1,3 = −2λ1 λ4 − 2λ4 , L1,5 = −2λ4 , (2) (2) L2,3 = 36λ1 λ4 + 8λ24 + 12λ4 , L2,5 = 8λ24 + 12λ4 , (2) (2) L3,3 = −540λ1 λ4 − 192λ24 − 108λ4 , L3,5 = −192λ24 − 108λ4 , (2) (2) L1,4 = −2λ3 λ2 , L1,6 = −2λ4 , (2) (2) L2,4 = 12λ3 λ2 , L2,6 = 32λ3 λ4 + 8λ24 + 12λ4 , (2) (2) L3,4 = −108λ3 λ2 , L3,6 = −582λ3 λ4 − 192λ24 − 108λ4 . 6 X (2) (2) We notice that we can not compute Lk by Lk,j as in the linear case, because of j=1 the repeated terms. The monomials having only one parameter, λl or λ2l , appear more than once in different Ej , while the monomials having two, λl λm , only once. So, for (2) every fixed k, before considering the sum of all Lk,j we need to correct the corresponding coefficient adequately. For example, the monomials λ4 and λ24 appear in S3 , S5 , and S6 , that is, exactly 3 = 31 times. Hence, these monomials will be divided by 3. But, as the monomial λ1 λ4 only appears in S3 the corresponding coefficient remains unchanged. So, we need to multiply by a correction factor which depends on the number of repetitions. As we have four perturbation parameters, we should divide each repeated term by 3 because is the number of times that they appear in S. (2) Denoting by L̂k,j the corrected k-Lyapunov constant corresponding to the pair Sj , we X 6 (2) (2) (2) can obtain Lk = L̂k,j . The statement follows because, in our case, we have L̂k,j = 0, j=1 for k = 1, 2, 3 and j = 1, 2 and (2) 2 (2) 2 L̂1,3 = −2λ1 λ4 − λ4 , L̂1,5 = − λ4 , 3 3 (2) 8 (2) 8 L̂2,3 = 36λ1 λ4 + λ24 + 4λ4 , L̂2,5 = λ24 + 4λ4 , 3 3 (2) (2) L̂3,3 = −540λ1 λ4 − 64λ4 − 36λ4 , L̂3,5 = −64λ24 − 36λ4 , 2 (2) (2) 2 L̂1,4 = −2λ2 λ3 , L̂1,6 = − λ4 , 3 (2) (2) 8 L̂2,4 = 12λ2 λ3 , L̂2,6 = 32λ3 λ4 + λ24 + 4λ4 , 3 (2) (2) L̂3,4 = −108λ2 λ3 , L̂3,6 = −582λ3 λ4 − 64λ24 − 36λ4 .  Now we can state the main result of this section, the computation in a parallelized form of the Taylor developments of the Lyapunov constants up to degree `, assuming that the number of total parameters is ` ≤ N. Theorem 3.3. Let p(z, z̄) and Qj (z, z̄), j = 1, . . . , N be polynomials with monomials of degree higher or equal than two such that the origin of ż = i z + p(z, z̄) is a center. For (`) ` ≤ N, we denote by Lk the Taylor approximation of k-Lyapunov constant up to degree ` of equation N X ż = i z + p(z, z̄) + λj Qj (z, z̄), (10) j=1 LOWER BOUNDS FOR THE LOCAL CYCLICITY OF CENTERS 11 with Λ = (λ1 , . . . , λN ) ∈ RN . Let S be the set of all combinations of the components of Λ taken ` at a time. That is, S = {(λ1 , . . . , λ` ), (λ2 , . . . , λ`+1 ), . . . , (λ`−N , . . . , λN )} and having N` elements. For each element Sj ∈ S, we denote by σ(j, ζ) the subscript of (`) the parameters in Sj at position ζ, i.e. Sj = (λσ(j,1) , . . . , λσ(j,`) ), and we denote by Lk,j the Taylor approximation up to degree ` with respect to Λ of the k-Lyapunov constant of equation X̀ ż = i z + p(z, z̄) + λσ(j,l) Qσ(j,l) (z, z̄). (11) l=1 Then N X (`) (`) Lk = L̂k,j , l=1 (`) P µk,j,p where L̂k,j = p   Λpj , for Λpj = λpσ(j,1) 1 λpσ(j,2) 2 · · · λpσ(j,`) ` and p = (p1 , . . . , p` ) writing N −s(p) `−s(p) ( (`) P p P̀ 1, if x > 0, Lk,j = p µk,j,p Λj with s(p) = sgn(pl ) where sgn(x) = l=1 0, if x = 0. Proof. It is well known that the k-Lyapunov constant, Lk , of a differential equation (10) is a polynomial in the parameters Λ = (λ1 , . . . , λN ), see [13]. Moreover, in our case each (`) Lk vanishes when Λ = 0. Consequently, Lk is the Taylor polynomial of Lk up to degree (`) ` at Λ = 0 and Lk (0) = 0. We write it as (`) X Lk = µk,p Λp , (12) p P N where Λp = λp11 λp22 · · · λpNN , p = (p1 , . . . , pN ), and pl ≤ `. This last condition allows us l=1 to decompose the total sum (12) in partial sums of ` parameters in Λ. Each partial sum is in fact the Taylor polynomial of the k-Lyapunov constant up to degree ` of equation (11). As each monomial can appear more than once in each partial sum, we need to correct the corresponding coefficient with a multiplicative factor that controls how many times this monomial appears. This factor depends on the number of different λl appearing in each monomial besides the total number N and the degree `, in fact, the number  of times that N −s(p) it appears in the partial sums, that is the combinatorial number `−s(p) .  4. Applications to cubic centers In this section, we study lower bounds for the local cyclicity of some Darboux cubic centers, in particular those called of codimension 12 in [35], but having a real center. We use the results and procedures described in Sections 2 and 3. Here we fix n = 3, the center point is located at the origin and rescaling variables and time if necessary, we consider the perturbed system (1), that writes in complex coordinates as 3 X ż = (i +λ0 )z + p(z, z̄) + (akl + i bkl )z k z̄ l . (13) k+l=2 Moreover, the corresponding unperturbed system ż = i z + p(z, z̄) has a cubic Darboux center and we have, in general and among the trace parameter λ0 , 14 real parameters (a20 , a11 , a02 , . . . , b20 , b11 , b02 , . . .) ∈ R14 . 12 LUIZ F. S. GOUVEIA AND JOAN TORREGROSA In the proofs, we have denoted by λ = (u1 , u2 , . . . , u11 ) ∈ R11 the relevant parameters and there are three that have been chosen as zero to simplify the computations. This is because in each center it can be seen that there are three that do not play any role. Using Theorems 2.1 and 2.2, we first study the number of limit cycles appearing with these parameters, and then we add an extra one using the trace λ0 . In each result, we have detailed why, the unperturbed system, is a center and which is the labeled name following the notation in [35]. In some cases, we have described also what was known up to now. Most of the considered centers are 1-parameter families depending on a special parameter a. Following the ideas of Theorem 2.1, see also [12], most probably, the cyclicity results are generic. That is, are valid for almost every value a, but for simplicity, we have chosen only one value for it. Proposition 4.1. There exist cubic polynomial perturbations such that from the origin of system   343 3 7 2 49  ẋ = − x − x − xy − 2 y, 576 72 36 (14)  ẏ = 343 x3 − 343 x2 y + 49 x2 + 259 xy − 49 y 2 + x,  72 96 24 36 9 bifurcate at least 11 limit cycles of small-amplitude. 12 This family is labeled as CD11 in [35] and as CD45 in [36]. System (14) has 5764801 5 588245 4 132055 3 4 81 23887872 x + 331776 x + 20736 x + 20335 2592 x2 + 1715 648 xy + 245 72 x + 35 9 y +1 H=−  5 (15) 245 117649 4 x + 2401 3 x + 1715 x2 + 49 x + 28 y+1 331776 1152 288 18 9 as a first integral. In both previous works, Melnikov theory or order 2 and 3, respectively, is used to prove the same statement. The first integral proposed in [35] depends on one parameter, named a, that it is fixed to zero in [36]. For this value the corresponding system has a real saddle equilibrium point at the origin but as the computations done by Żola̧dek are all in complex he comments that, generically, the result for real cyclicity is also valid. Recently, in [30] this special parameter value is fixed to a = −3 where the system has a real center at the point (3/2, −11/4). The proof uses developments of seventh-order of the Lyapunov constants. The first integral (15) is obtained from the one proposed in [30] but moving the equilibrium point and doing an affine change of coordinates that writes the linear part simpler. This new expression reduces computation time. Our proof, which is different, shows also that a Taylor approximation up to degree 7 is the minimum necessary to unfold 11 limit cycles but we have used developments of the Lyapunov constants up to degree 10 for a better understanding of this local bifurcation phenomenon. Proof of Proposition 4.1. Clearly system (14) has a center at the origin because (15) is a first integral well defined at the origin and the corresponding level curves in a neighbor- hood of the origin are ovals. In fact H = −81/245 + x2 + 2y 2 + · · · . The first step is the computation of the Taylor approximations of the Lyapunov con- stants, Lk , corresponding to a cubic perturbation of system (14) having only quadratic and cubic terms. We will study which are the principal parts, near the origin in the parameters space, of each Lyapunov constant when the previous vanish. The second step is the study of linear parts of Lk . Unfortunately, we get only 9 linearly independent. More concretely, the 10th and the 11th are linearly dependent with respect to the first nine. Consequently, Theorem 2.1, adding the trace parameter, provides an unfolding with 9 limit cycles. To obtain the remaining two as it is stated, we need to look at higher-order developments. After a linear change of variables and doing the necessary LOWER BOUNDS FOR THE LOCAL CYCLICITY OF CENTERS 13 transformations following the scheme detailed in Section 2, we get Lk = uk + O2 (λ), for k = 1, . . . , 9, L10 = O2 (λ), L11 = O2 (λ). Here λ = (u1 , u2 , . . . , u11 ) ∈ R11 denotes the relevant parameters, the other are zero, and Ok (λ) contains all monomials of degree at least k in λ. In the third step, we study the higher-order developments, in particular, the ones corresponding to L10 and L11 . The Implicit Function Theorem ensures that there exists an analytical local change of variables in the parameter space, well defined in a neighborhood of the origin, such that Lk = vk , for k = 1, . . . , 9. With this change, doing an affine change of variables, and following the scheme detailed in Section 2 we can write the last two Lyapunov constants depending only on two parameters (v10 , v11 ), except positive multiplicative constants, as 3 L10 = v10 + O4 (v10 , v11 ), 3 L11 = −v10 + O4 (v10 , v11 ). Here, taking v11 = 0, it is clear that we have a curve, in the parameters space, of weak-foci of order 10, because L10 6= 0, that unfolds (using the linearity of v1 , . . . , v9 and the trace parameter) 10 limit cycles. This is also true taking v10 = 0, but with a higher development, 6 because L10 6= 0. In fact, except a positive multiplicative constant, L10 = v11 + O7 (v11 ). It is evident that, only with a third-order development, we can not unfold the eleventh limit cycle. The next step is the study, if they exist, of the different real branches near the origin that has the algebraic curve L10 = 0 and if there exists one such that L11 6= 0. This is done computing more terms in the development of L10 and L11 . We can follow [8] to use the Newton–Puiseux algorithm with different weights for the variables (v10 , v11 ). With weights (2, 1) the principal parts are 3 2 2 4 6 (2,1) L10 = v10 + α22 v10 v11 + α14 v10 v11 + α06 v11 + O7 (v10 , v11 ), 3 2 2 4 6 (2,1) L11 = −(v10 + α22 v10 v11 + α14 v10 v11 + α06 v11 ) + O7 (v10 , v11 ), (2,1) where Om denotes the monomials with degree higher or equal than m with respect to the weight (2, 1). With this weight, it is easy to see the different simple branches of the intersection at the origin. In fact, the principal part of L10 , with this specific choice of weights, decomposes as a product of two factors, one simple and one double, 2 2 2  (2,1) L10 (v10 , v11 ) = (v10 − a2 v11 ) v10 − b2 v11 + O7 (v10 , v11 ), (16) with 660160595890746 487045680336990 a2 = − and b2 = − . 37506906889 37506906889 [j] The study of the different branches, v10 (v11 ), of the curve L10 (v10 , v11 ) = 0 near the origin 2 is done using the weighted blow-up v10 = v11 w10 . Then the function (16), collecting in v11 , can be written as polynomials in v11 of degree m with coefficients polynomials in w10 of degree [m/2], for m ≥ 6, where [·] denotes the integer part function. Consequently, 6 dividing by v11 , (16) writes as ∞ X e10 (w10 , v11 ) = (w10 − a2 ) (w10 − b2 )2 + L m Wm (w10 )v11 , m≥1 14 LUIZ F. S. GOUVEIA AND JOAN TORREGROSA with Wm polynomials of degree [(m + 6)/2]. An equivalent expression, L e11 (w10 , v11 ) can be obtained for L11 but with different functions Wm (w10 ). Now, we can write Le10 = w e10 in a neighborhood of (w10 , v11 ) = (a2 , 0) using the Implicit Function Theorem, because L e10 (a2 , 0) = 0 and the partial derivative ∂L10 = (a2 − b2 )2 6= 0. ∂w10 (a2 ,0) Clearly, when we10 = 0, there exists an analytic branch, w10 = ω [1] (v11 ) = a2 + O1 (v11 ), e10 and also L10 vanish on it. Consequently, v10 = v 2 ω [1] (v11 ) = a2 v 2 +O3 (v11 ). such that L 11 11 6 6 With this change of variables, after multiplying by v11 , we write (16) as L10 = v11 we10 . Additionally, except a positive multiplicative constant, we get 7 L11 |L10 =0 = L11 |we10 =0 = −v11 + O8 (v11 ) 6= 0, 7 for v11 6= 0 small enough, and L11 = −v11 e10 , v11 )). We notice that the neighbor- (1 + O1 (w hood of (a2 , 0) for the variables (w10 , v10 ) is transformed to a neighborhood of the origin for the variables (we10 , v11 ). This proves that there exists a curve, in the parameters space, for v11 small enough, of weak-foci of order 11 that is born at the origin. The cubic perturbation mechanism described in Section 2 following the scheme of Roussarie, see [28], proves that only 11 limit cycles can bifurcate from the origin of system (14). This is because, from the Weierstrass Preparation Theorem, see [31], the Poincaré return map is a polynomial with coefficients the Lyapunov constants, that write, after all the changes of variables, as Lk = uk , for 6 7 k = 1, . . . , 9, L10 = v11 e10 and L11 = −v11 w (1 + O1 (w e10 , v11 )). Finally, using the weighted blow-up {u1 = z 7 z1 , u2 = z 7 z2 , . . . , u9 = z 7 z9 , w e10 = zz10 , u11 = z}, after dividing by z 7 , it is clear that we have constructed a versal unfolding, obtaining the maximal number of limit cycles. As v11 can vanish this upper bound is in fact a lower bound for the cyclicity of the center as we wanted to prove.  As we have explained in Section 2, the above proof is a kind of generalization of the result provided by Christopher in [12], see Theorem 2.2. After the simplification due to the linear developments of the first 9 Lyapunov constants, the next two has not the same order at the origin. It is also clear that the transversal straight line ` is now an analytic curve. We notice that, as we have proved that the curve L10 = 0 has a real branch associated to the simple factor, using the Weierstrass Preparation and Division Theorems, [31], we can write (16), except for a nonvanishing multiplicative function, as 2 (v10 − v11 ω [1] (v11 ))(v10 2 2 + v11 4 φ1 (v11 )v10 + v11 φ0 (v11 )), (17) where φ0 and φ1 are analytic functions that vanish at zero. As we have computed the Taylor approximations of the Lyapunov constants up to degree 10, using the Puiseux series (see [8]), we can compute four extra terms of this analytic simple branch v10 = 2 P ∞ v11 ω [1] (v11 ) = m am v11 . In this case, the Puiseux series has only natural exponents then m≥2 LOWER BOUNDS FOR THE LOCAL CYCLICITY OF CENTERS 15 it is, in fact, a Taylor series. Straightforward computations show that 16104570945819692121638226351 a3 = , 25209283713691672597112320 386258251571578220793485476718239267056083732367 a4 = , 123628442038275561958744770186426115741450240 7344527305232752838312438300220617202784745335855878366073914837 a5 = − , 10802169913777097568319537676224431667868634439810816147456000 1040667719410212727048282608984 · · · 195082403446045341049879843107 a6 = − . 2913611107447792035945464752572 · · · 844814820535324990043586560000 Using the above coefficients and the developments up to degree 10 of L10 we can compute the first terms of the Taylor series of the functions φ0 and φ1 in (17). It can be checked that they provide the double factor appearing in (16) and that the discriminant, with respect to v10 , writes as φ21 (v11 ) − 4φ0 (v11 ) = Av11 8 + · · · with A > 0. This proves that there 2 exist another two real branches, tangent to the double factor in (16), v10 = v11 ω [j] (v11 ) = P [j] m ∞ bm v11 , for j = 2, 3. In the above proof we have got the first coefficient, that coincides m≥2 [2] [3] for both branches, b2 = b2 = b2 . In fact, the second also coincides, [2] [3] 30834092507446246450289832 b3 = b 3 = b3 = − , 9847376450660809608247 [2] [3] because the discriminant starts with terms of degree 8. The next coefficients b4 and b4 , which are different, are the zeros of the quadratic equation β 2 + β1 β + β0 = 0, where 374029705710452551852715772722429520 β1 = , 369344631635866288497876795513583 8566777417085455703175077527180594435196203104777368650266555098833382208 β0 = . 33421786944967271392447900695765062189771656497217124612800686982805 We observe that we have needed approximations up to degree 10 to distinguish the other two real branches. Over them, we have checked that L11 vanishes up to this level of approximation. Proposition 4.2. There exist cubic polynomial perturbations such that from the origin of system   5 3 15 2 5  ẋ = x − x y− xy 2 − y, 32 64 32 (18)  ẏ = 15 x3 + 35 x2 y − 15 xy 2 − 15 y 3 + x,  64 32 16 32 bifurcate at least 11 limit cycles of small-amplitude. 12 This is the system CD12 in [35]. Christopher studies in [12] the local cyclicity for this family fixing the free parameter to a = 2. He also needs developments up to degree two, but with a high computational effort using Grobner Basis. The above system is obtained taking a = 3/5, but doing an adequate affine change of coordinates. The computations are simpler because the linear part at the equilibrium, which has been moved to the origin, is in the normal form of an elementary center type point. Proof of Proposition 4.2. System (18) has a center at the origin because it has the rational first integral 36 (625x4 + 1920x2 + 2560xy + 4096)3 H=− , 5 (78125x6 +360000x4 +480000x3 y+1044480x2 +737280xy+491520y 2 +786432)2 which is well defined at the origin and with Taylor series −4/5 + x2 + y 2 + · · · . 16 LUIZ F. S. GOUVEIA AND JOAN TORREGROSA Computing the Lyapunov constants, following the scheme explained in Section 3 up to degree 1, corresponding to the perturbed system (13), we have that the first five linear terms have rank five with respect to the parameters. Doing a linear change of (1) coordinates we have that they write Lk = uk , for k = 1, . . . , 5. Using the triangularization procedure describe before Proposition 2.3 we can simplify the next Lyapunov constants (1) to get Lk = 0, for k = 6, . . . , 11. Consequently, up to first-order only five limit cycles bifurcate from the origin. The second step is the computation of the second-order terms. Using again Proposi- tion 2.3 to eliminate the parameters uk , k = 1, . . . , 5, we obtain 6 homogeneous polyno- mials of degree 2 for the second-order terms of Lk for k = 6, . . . , 11. (2) 65234375 2 5695234375 187109375 1842265625 L6 = a11 − a02 a11 + a11 a20 + a11 b02 2491416576 226718908416 32388415488 75572969472 250390625 23046875 390625 2 1466328125 − a11 b11 + a11 b20 + b11 − a02 b11 5813305344 1937768448 899678208 75572969472 250390625 950546875 266796875 141015625 2 − a20 b11 − b02 b11 + b11 b20 − a02 75572969472 75572969472 25190989824 18893242368 21171875 20703125 708828125 390625 + a02 a20 − a02 b02 + a02 b20 + a220 4359979008 25190989824 75572969472 1349517312 204921875 23046875 241796875 2 47421875 11328125 2 − a20 b02 + a20 b20 + b02 + b02 b20 − b20 , 75572969472 25190989824 75572969472 37786484736 3598712832 (2) 804296875 2 62870546875 229296875 488059296875 L7 = − a11 + a02 a11 − a11 a20 − a11 b02 53150220288 4836670046208 76772540416 43530030415872 221214453125 240866796875 94140625 406132109375 + a11 b11 − b20 a11 − b211 + a02 b11 14510010138624 43530030415872 345476431872 43530030415872 8782421875 3529140625 26304296875 9333203125 2 + a20 b11 + b02 b11 − b11 b20 + a02 3348463878144 537407782912 4836670046208 3627502534656 1203359375 11367578125 163425546875 12109375 2 − a02 a20 + a02 b02 − a02 b20 − a20 329772957696 43530030415872 43530030415872 86369107968 39773359375 5510546875 126171875 2 10987578125 8826953125 + a20 b02 − a20 b20 + b20 − b02 b20 − b202 , 14510010138624 4836670046208 76772540416 7255005069312 14510010138624 (2) 205818359375 2 27010282421875 566685546875 L8 = a11 − a02 a11 + a11 a20 68341519613952 15786891030822912 2255270147260416 6561098046875 49375697265625 2022599609375 + a11 b02 + a11 b11 + b20 a11 15786891030822912 15786891030822912 15786891030822912 471892578125 6456421484375 4613353515625 + b211 − a02 b11 − a20 b11 5920084136558592 15786891030822912 5262297010274304 47627121484375 33662533203125 733501953125 − b02 b11 + b11 b20 + a202 47360673092468736 47360673092468736 1315574252568576 1488928515625 3438419921875 5852162890625 + a02 a20 + a02 b02 − a02 b20 1315574252568576 15786891030822912 15786891030822912 5248046875 121922265625 30685546875 − a220 − a20 b02 + a20 b20 140954384203776 103182294319104 53154515255296 3643357421875 509573828125 3731908203125 2 − b220 + b02 b20 − b02 , 15786891030822912 607188116570112 5262297010274304 (2) 4251095908203125 2 1670155486328125 120403642578125 L9 = a11 − a02 a11 + a11 a20 2742408499068665856 914136166356221952 210954499928358912 1678778216796875 43118063095703125 9782838056640625 + a11 b02 − a11 b11 + b20 a11 747929590655090688 8227225497205997568 8227225497205997568 144507529296875 9930148478515625 240089833984375 + b211 − a02 b11 − a20 b11 9598429746740330496 5235507134585634816 2742408499068665856 757186115234375 14849983759765625 377087236328125 2 − b02 b11 + b11 b20 − a02 914136166356221952 19196859493480660992 436292261215469568 143162263671875 1169718330078125 62189059427734375 + a02 a20 − a02 b02 + a02 b20 685602124767166464 4430044498495537152 57590578480441982976 2179873046875 9000880859375 210343720703125 + a220 + a20 b02 − a20 b20 28566755198631936 210954499928358912 2742408499068665856 1433194228515625 2 1701270740234375 1006989326171875 2 − b20 − b02 b20 + b02 , 6398953164493553664 9598429746740330496 1745169044861878272 LOWER BOUNDS FOR THE LOCAL CYCLICITY OF CENTERS 17 (2) 16403503182763671875 2 6019489875693359375 4110241175146484375 L10 = − a11 + a02 a11 − a11 a20 11057391068244860731392 4553043381042001477632 11057391068244860731392 226410773740234375 1537815485986328125 1853450635595703125 − a11 b02 + a11 b11 − b20 a11 180423630484181876736 781835732098121465856 2866731017693112041472 543782850048828125 27828407758955078125 719020527197265625 − b211 + a02 b11 + a20 b11 19350434369428506279936 25800579159238008373248 2866731017693112041472 50925810357353515625 43579580815966796875 1929508851220703125 2 + b02 b11 − b11 b20 + a02 77401737477714025119744 77401737477714025119744 6450144789809502093312 7638986422255859375 54847838720703125 103268798154296875 − a02 a20 + a02 b02 − a02 b20 19350434369428506279936 505893709004666830848 216811589573428641792 17619605615234375 6869360360732421875 279467075634765625 − a220 + a20 b02 − a20 b20 460724627843535863808 25800579159238008373248 2866731017693112041472 4340212115087890625 15956755576171875 3536245476904296875 + b220 − b02 b20 − b202 , 25800579159238008373248 159262834316284002304 25800579159238008373248 (2) 42718448752197265625 2 171788232009716796875 6210508852783203125 L11 = a11 − a02 a11 + a11 a20 95184091534832953196544 602832579720608703578112 95184091534832953196544 5259382308953564453125 2996301505219970703125 2217310253376220703125 + a11 b02 + a11 b11 + b20 a11 37978452522398348325421056 12659484174132782775140352 37978452522398348325421056 79495114489990234375 4526167344325244140625 1501926097403076171875 + b211 − a02 b11 − a20 b11 6329742087066391387570176 37978452522398348325421056 12659484174132782775140352 2113627250105517578125 8402376020751953125 170019216537353515625 − b02 b11 + b11 b20 + a202 12659484174132782775140352 67697776332260870455296 3164871043533195693785088 77972993629345703125 378160348073486328125 344733712266748046875 + a02 a20 + a02 b02 − a02 b20 452124434790456527683584 37978452522398348325421056 37978452522398348325421056 58214882568359375 2056445968159619140625 56233442866943359375 + a220 − a20 b02 + a20 b20 32294602485032609120256 12659484174132782775140352 744675539654869575008256 371336987060546875 666554572872607421875 930355010247314453125 − b220 + b02 b20 − b202 . 9568771106676328628224 6329742087066391387570176 12659484174132782775140352 (2) (2) We consider now the system {L6 = · · · = L10 = 0}. Doing, for example, the blow-up a02 = zv1 , a11 = zv2 , a20 = zv3 , b02 = z, b11 = zv4 , b20 = zv5 , we can solve a system of five equations of degree 2 with respect to 5 variables. Using a computer algebra system we get that vk = pk (α)/q(α) with pk and q polynomials with rational coefficients of degree 27 and α being a solution of a given polynomial, Q(α), also with rational coefficients, of degree 28. This polynomial has 20 simple real solutions, { − 1.460571830, −0.6718444255, −0.6670390163, −0.5158935998, −0.4874999611, − 0.3970369469, −0.3874233159, −0.2401990480, −0.02992848475, 0.02384205186, 0.03979267840, 0.08015376288, 0.2087598950, 0.2131172755, 0.2232320471, 0.2463926997, 0.2995004189, 0.3312032992, 0.3788127882, 1.397031032}. The next step is to check that the Jacobian of the five equations with respect to {v1 , . . . , v5 } is different from zero. The size of the polynomials does not make possible to compute directly the determinant of the Jacobian matrix in terms of α. Then, we do a Gauss elimination to get a triangular matrix. The diagonal elements are now of the form J11 = q1 (α)/q(α) and Jkk = qk (α)/(qk−1 (α)q(α)), for k = 2, . . . , 4, where all polynomials have also rational coefficients and degree 27 in α, in particular q1 = p1 . Consequently, (2) the determinant is q5 (α)/q 5 (α). Moreover, it can be checked also that L11 = p6 (α)/q 2 (α). The last step is the computation of the resultants of all polynomials pk , qk with Q with respect to α, checking that all are different from zero. So all the variables vk , among of (2) the determinant and the value of L11 are nonzero real numbers. It is clear that with the proposed blow-up, after dividing by z 2 we can apply the Implicit (2) (2) Function Theorem to find an analytic curve vk = Vk (z) where the varieties {L6 , . . . , L10 } (2) intersect transversally and L11 is nonzero. Then 6 extra limit cycles appear and the statement follows. The scheme of the proof is the same as the proof of Proposition 4.1. But, in fact, we could apply Theorem 2.2. But the work to find the line ` is the same. 18 LUIZ F. S. GOUVEIA AND JOAN TORREGROSA We notice that all the computations in this proof have been made with Maple on a personal computer in a few minutes. We have not shown the polynomials because of the size of them. The coefficients are rational numbers with numerators and denominators having more than 200 digits.  Remark 4.3. An alternative proof of the above result can be done computing numerically an approximation of a solution v1 = 0.1924453833548429, v2 = 0.2384205185830677, v3 = 0.1490077870024313, v4 = 0.7626068770346651, v5 = 1.5437258992144801, (2) with enough digits to ensure that L6 /z 2 ≈ 6.80504529555082 · 10−9 and the Jacobian, −4.45823335837756 · 10−10 , are nonzero real numbers. This can be obtained using a computer assisted proof with the Poincaré–Miranda Theorem, as we will do in some of the following results. The next proposition shows the difficulties to get more than 10 limit cycles. Our computations do not provide a better result for the cyclicity of the next cubic polynomial 12 system. The unperturbed system is labeled as CD29 in [35] but we have not considered it directly. To simplify the computations, we have made an affine change of coordinates. Proposition 4.4. There exist cubic polynomial perturbations such that the origin of sys- tem ( ẋ = 8x3 − 40x2 y + 2x2 − 30xy − 5y, 24 3 (19) ẏ = x + 24x2 y − 80xy 2 + 4x2 + 10xy − 10y 2 + x, 5 has cyclicity at least 10. Proof. System (19) has a center at the origin because it has the rational first integral (64x3 + 72x2 + 120xy + 30x + 30y + 5)4 , (128x4 + 192x3 + 320x2 y + 128x2 + 240xy + 40x + 40y + 5)3 which is well defined at the origin and the level curves are topologically circumferences. Doing a first-order analysis of the Lyapunov constants, only the first 6 are linearly (1) independent. Then, as the previous studies, we can write Lk = uk , for k = 1, . . . , 6. Applying the simplification algorithm described previously, using Proposition 2.3, we get (2) L7 = A7 u7 u8 , (3) L8 = A8 u7 u9 u10 , (3) L9 = u7 u9 (A9 u7 + B9 u10 ), (3) L10 = u7 u9 (A10 u7 + B10 u10 ), (3) L11 = u7 u9 (A10 u7 + B10 u10 ), (3) (3) where Ak , Bk are nonvanishing rational numbers. We notice that L11 = L10 . To study how is the local intersection of the varieties Lk , for k = 7, . . . , 11 we need to do an adequate weighted blow-up using a privileged parameter. Here we have chosen u7 = z, u8 = z 2 z8 , u9 = zz9 , u10 = zz10 , u11 = zz11 , the other three parameters are taken as zero. The Taylor LOWER BOUNDS FOR THE LOCAL CYCLICITY OF CENTERS 19 development with respect to z, dividing by a nonzero rational number, is X L7 = z 3 (z8 + p2 (z9 , z10 )) + W7j (z8 , z9 , z10 , z11 )z j , j≥4 X L8 = z 3 z9 z10 + W8j (z8 , z9 , z10 , z11 )z j , j≥4 X L9 = z 3 z9 p1 (z9 ) + W9j (z8 , z9 , z10 , z11 )z j , j≥4 X L10 = z 3 z92 + W10j (z8 , z9 , z10 , z11 )z j , j≥4 where p1 and p2 are polynomials of degree 1 and 2, respectively, p1 (0) 6= 0, and Wk,j are also polynomials. The statement follows as in the previous studies because there exists a transversal intersection of the varieties {z8 + p2 (z9 , z10 ) = 0, z9 z10 = 0, z9 p1 (z9 ) = 0} with z9 6= 0. Then, clearly, L10 is nonvanishing for z small enough.  Remark 4.5. We remark that we have tried to improve the above result unsuccessfully. We have used different weighted blow-ups, higher-order developments, and the study of the Newton polyhedron. Proposition 4.6. There exist cubic polynomial perturbations such that from the origin of system   338 3 4394 2 338 2 338 3 26 2 793 65  ẋ = − x + x y− xy + y − x − xy + y 2 + y, 441 9261 147 147 21 441 21   ẏ = 338 338 338 338 65 247 26 2 x3 − x2 y − xy 2 − y3 − x2 − xy − y − x, 3087 343 1323 1029 147 147 147 (20) bifurcate at least 11 limit cycles of small-amplitude. Proof. System (20) has a center at the origin because it has the rational first integral   3 2704 x2 5408 xy 2704 y 2 208 x 416 y 104 x 104 y 324135 3087 − 3087 + 1323 − 147 + 147 + 1 441 + 147 + 1 H=   2 86528 x2 y2 − 2704 15435 + 2704 1715 − 52 147 x + 5221y + 1 which is well defined at the origin, H = 324135/86528 + x2 + y 2 + · · · . We consider the perturbed system (13) being (20) the unperturbed center. Follow- ing the scheme explained in Section 3, the developments up to degree 1 of the first 8 Lyapunov constants are linearly independent. Then, after a linear change of the pertur- (1) bation parameters, we have that they write as Lk = uk , for k = 1, . . . , 8. Using the properties detailed in Proposition 2.3 we can simplify the next Lyapunov constants to (1) get Lk = 0, for k = 9, . . . , 11. At this point we assume, to simplify computations, that b03 = b12 = b30 = 0. Then, we can write (2) L9 = −Au9 u10 , (2) L10 = Bu9 u10 , (2) L11 = −Cu9 u10 , where A, B, and C are rational numbers having between 48 and 71 digits in the numerators and between 68 and 86 digits in the denominators. So, we see clearly here that we have at 20 LUIZ F. S. GOUVEIA AND JOAN TORREGROSA least 9 limit cycles using the trace parameter together with uk , for k = 1, . . . , 8. In fact, at most 9 with only order 2 developments. Hence, if we want more limit cycles, we should compute up to degree 4 developments. (3) (3) Because, after using again the simplification with L9 and Proposition 2.3, L10 = L11 = 0. Doing the blow-up u9 = z, u10 = zz10 , and u11 = zz11 , we can divide L9 by z 2 and we can use the Implicit Function Theorem to write z10 as a function of z and z11 . Then, (4) (4) L10 = z 4 p4 (z11 ) and L11 = z 4 q4 (z11 ) where p4 and q4 are polynomial with rational coefficients of degree 4 having both 4 different real roots. Moreover, z 4 is a common factor of the complete L10 and L11 and the resultant of p4 and q4 with respect to z11 is different from zero. Therefore, applying also the Implicit Function Theorem near the simple zeros of p4 , there exists values of the parameters, for small enough z, such that L10 = 0 and L11 6= 0. The statement follows because we have proved the existence of an analytic curve of weak-foci of order 11 such that, along it and as in the previous proofs, 11 limit cycles of small-amplitude bifurcate from the origin.  Proposition 4.7. There exist cubic polynomial perturbations such that from the origin of system ( 2809 ẋ = 32946 x3 − 22472 49419 5618 yx2 + 16473 69695 2 xy 2 + 98838 x − 66992 16473 xy + 636 289 y 2 + x − 2112 289 y, 22472 2809 86549 20988 151447 ẏ = 247095 xy 2 − 16473 y3 + 197676 xy − 27455 y2 + 247095 x − y, (21) bifurcate at least 11 limit cycles of small-amplitude. 12 This is the case labeled CD31 in [35]. Christopher in [12] provides with it the first analytic proof that 11 limit cycles of small-amplitude exist for a cubic polynomial vector field. Here we add it for completeness. Proof of Proposition 4.7. The corresponding first integral of (21), which is well defined at the origin, is  5 2 280 xy 342 y 2 22409 x 95760 y 7812755 xy + + + + + 53 53 2809 2809 148877 H=  3 , 342 2 x+ F6 (x, y) 53 where 700 4 342 5 406045 3 239400 4 30389450 2 F6 (x, y) = xy 5 +xy + y + xy + y + xy 53 53 5618 2809 148877 139611775 3 18788141215 10549512750 2 + y + xy + y 297754 63123848 7890481 150246782525 826646189040 26977377387858 + x+ y+ . 836390986 418195493 22164361129 The proof of the above proposition follows computing the linear terms of the Lyapunov constants and then using Theorem 2.1 to provide the complete unfolding of 11 limit cycles.  5. Order one studies to get lower bounds for M (8) and M (9) This section is devoted to proving the statement of Theorem 1.1 corresponding to local cyclicity of polynomial vector fields of degrees 8 and 9, using only linear developments. The proofs follow from Theorem 2.1 just computing the Taylor approximations up to degree LOWER BOUNDS FOR THE LOCAL CYCLICITY OF CENTERS 21 1 and computing how many are linearly independent. In both results, the unperturbed systems are centers of degrees 7 and 8 having a straight line of equilibrium points, that, for simplicity, we have fixed to {1 − x − y = 0}. We notice that the parallelization procedure described in Section 3 is indispensable to get the results. The total computation time, in both cases, is less than one hour. Proposition 5.1. Consider the perturbed system (4) of degree n = 8 with the center (ẋ, ẏ) = (Pc (x, y), Qc (x, y)) given by    ẋ = (1 − x − y)(− 2527 3 x6 y − 2968 3 x5 y 2 − 4186 3 x4 y 3 − 2800 3 x3 y 4 − 553x2 y 5     +56xy 6 + 184 x3 y + 88 x2 y 2 + 48xy 3 − y), 3 3   ẏ = (1 − x − y)(672x7 + 1484x6 y + 2219 x5 y 2 + 5684 x4 y 3 − 742 x3 y 4 + 1148 2 5 xy   3 3 3 3   −315y 6 − 28y 7 − 58x4 − 44x3 y − 104 3 x2 y 2 − 44 3 xy 3 + 10y 4 + x). There are perturbation parameters λ such that at least 76 limit cycles of small-amplitude bifurcate from the origin. The above system, without the straight line of equilibrium points, is a center because it has the first integral H(x(x2 + y 2 ), y(x2 + y 2 )) where (42x − 7y − 1)3 f3 (x, y) H(x, y) = (22) (448x2 + 336xy + 63y 2 − 44x − 12y + 1)3 (1183x2 − 68x + 1) and f3 (x, y) = 10752x3 + 29568x2 y + 17640xy 2 + 3024y 3 − 1600x2 − 2760xy − 576y 2 + 74x + 57y − 1. The rational first integral (22) corresponds to the cubic polynomial center provided by Bondar and Sadovskiı̆ in [7]. They prove that the cubic perturbations provide also 11 limit cycles using only up to degree 1 Taylor approximations as for system (21). Although we have also verified it, we have not included this result here. Proposition 5.2. Consider the perturbed system (4) of degree n = 9 with the center (ẋ, ẏ) = (Pc (x, y), Qc (x, y)) given by  54 8 18 7 54 6 2   ẋ = (1 − x − y)( 175 x + 35 x y − 175 x y + 894 175 66 3 5 x5 y 3 − 2x4 y 4 + 25 xy     26 2 6 − 35 x y − 342 xy 7 + 16 y 8 − y), 175 25   198 7 ẏ = (1 − x − y)(− 175 x y − 1254 x6 y 2 − 586 5 3 xy − 258 4 4 xy − 22 3 5 xy   175 175 35 5   18 2 6 + 25 x y − 382 175 162 8 xy 7 + 175 y + x). There are perturbation parameters λ such that at least 88 limit cycles of small-amplitude bifurcate from the origin. The proof that the above system, without the straight line of equilibrium points, is a center follows from an idea of Giné in [19]. We consider the center with homogeneous quartic nonlinearities given in [19, System (6), Pag. 8857] taking c = 4/5 and s = 3/5. The change of variables (x, y) = r3/7 (cos θ, sin θ) in such a quartic system gets a system of degree 8 having also a center at the origin. We notice that for other degrees, n = 3, . . . , 7, adding a straight line of equilibria to a center of degree n − 1, we have not obtained higher lower bounds for the local cyclicity than the ones obtained previously, nor better than the ones given in the results of the next section using higher-order Taylor series. For example, the best cubic system of Section 4 adding such curve provides a quartic system with only 19 limit cycles up to first and second-order studies. For degree 6, we have not found any system to improve the highest value found in [4], M (6) ≥ 48, which improves in 8 the best known value for M (6) in [25]. We remark that the number of total parameters for n = 6 is 50. 22 LUIZ F. S. GOUVEIA AND JOAN TORREGROSA 6. Higher order studies to get lower bounds for M (4), M (5), and M (7) This section is devoted to proving the statement of Theorem 1.1 for the local cyclicity of polynomial vector fields of degrees 4, 5, and 7. In the proofs, we will use higher-order Taylor developments. Proposition 6.1. Consider the perturbed system (4) of degree n = 4 with the center (ẋ, ẏ) = (Pc (x, y), Qc (x, y)),   1  ẋ = (32x4 − 168x2 y 2 + 32xy 3 + 24y 4 ) − y, 25 (23)  ẏ = 1 (−32x3 y − 192x2 y 2 − 24xy 3 + 64y 4 ) + x.  25 There are perturbation parameters λ such that at least 20 limit cycles of small-amplitude bifurcate from the origin. Proof. The system in the statement is presented in [19, 20]. In [19] it is proved that the origin is a center because it has the polynomial inverse integrating factor V (x, y) =(8xy − 4y 2 + 10y − 5)(64x3 − 96x2 y + 48xy 2 − 8y 3 − 25)× (64x2 y 2 − 64xy 3 + 16y 4 − 80xy 2 + 40y 3 + 40xy + 80y 2 + 50y + 25) We need to compute the Taylor developments up to degree 4 of the Lyapunov constants. Because we will see that up to degree 3 is not enough to prove the statement. As we have detailed in the proofs of Propositions 4.1 and 4.2, simplifying as we have described in Section 2 and up to multiplicative nonzero constants, we have the next triangularized expressions up to degree two: Lk = uk + O2 (λ), for k = 1, . . . , 16, L19 = O3 (λ), L17 = −u17 u19 + u218 + O3 (λ), L20 = O3 (λ), (24) L18 = u17 u18 + O3 (λ), L21 = u220 + O3 (λ). Here λ denotes the perturbation parameters and Ok (λ) are the monomials of degree at least k in λ. Clearly, using only up to degree two of Taylor approximations there are two Lyapunov constants that vanish. So, higher degrees must be computed. Using the Implicit Function Theorem and after using adequately the expressions of L1 , . . . , L16 , we can consider only Lk , for k = 17, . . . , 21 depending on (u17 , . . . , u24 ). With the Taylor approximation up to degree three, an adequate blow-up is u17 = z, u18 = z 2 w1 , u19 = z 2 w2 , u20 = zw3 , u21 = zw4 , u22 = zw5 , u23 = zw6 , u24 = zw7 . Then, after dividing by z 3 , we have e(3) 1 79 1 L 17 = w2 + w32 + w62 + w3 w6 , 3402000 4838400 241920 e(3) 387 9 1 1 L 18 = w1 − w6 w7 + w4 w6 + w5 w6 − w3 w7 39040 2440 1220 1600 22797130436674460681587504742109808816907484398453930617769 − w6 22599395017588000741825603692446552868781675590536155251513600 9999644540025045050531707316918826074133709626587611529 − w3 , 217079229805212609584646551861359050353307590841471778158750 e(3) 54197 · · · 00000 362945 · · · 0000 L 19 = w3 + w7 , 866211 · · · 28447 28873 · · · 76149 e(3) L 20 = 0. LOWER BOUNDS FOR THE LOCAL CYCLICITY OF CENTERS 23 From the above expressions is clear that only up to degree three is not enough to prove the statement. But, as the above nonvanishing three terms have rank 3 with respect to e(3) = zk , with w1 , w2 , and w3 , we can use also the Implicit Function Theorem to solve L k respect to w1 , w2 , w3 and compute the Taylor approximation up to degree four we get (3) (3) (3) (4) L17 = z 3 z17 , L18 = z 3 z18 , L19 = z 3 z19 , L20 = z 4 w62 . From which it follows that a curve of weak-foci of order 20 exists. Moreover, it unfolds 20 limit cycles of small-amplitude at it is indicated in the statement. The proof finishes as the proofs of Propositions 4.1 and 4.2.  Remark 6.2. In [19], it is studied the homogeneous nonlinearities perturbation of degree 4 of (23). It is proved, with approximations up to degree two, that 7 limit cycles exist in this special subclass. This proves that there exists a curve of weak-foci of order 21. This can be seen intuited from (24) but it is necessary to prove, using only homogeneous degree four perturbation terms, that L21 is nonvanishing. But it is not possible to find a complete unfolding of the 21 limit cycles with degree two because there are no free parameters, in (2) (2) fact L19 = L20 = 0. In the proof, we have shown that with degree three we get only 19 limit cycles and with degree four we get 20. The problem about the existence of a complete unfolding of 21 remains open. The computations to go further in the order are very hard. Proposition 6.3. Consider the perturbed system (4) of degree n = 5 with the center (ẋ, ẏ) = (Pc (x, y), Qc (x, y)),   1   ẋ = (42x5 − 12x4 y − 476x3 y 2 − 68x2 y 3 + 266xy 4 + 56y 5 ) − y, 25 1 (25)   ẏ = (−8x 5 − 26x 4 y − 28x 3 2 y − 4x 2 3 y − 132xy 4 + 6y 5 ) + x.  25 There are perturbation parameters λ such that at least 33 limit cycles of small-amplitude bifurcate from the origin. Proof. As the previous result, system (25) appears also in [19, 20]. The system has a center because V (x, y) =(64x8 + 1600x7 y + 13456x6 y 2 + 32000x5 y 3 − 99616x4 y 4 − 380800x3 y 5 + 320656x2 y 6 + 548800xy 7 + 153664y 8 + 20000x4 + 85000x3 y − 10000x2 y 2 − 155000xy 3 − 100000y 4 + 15625)× 1 (−64x4 + 192x3 y − 16x2 y 2 − 192xy 3 − 64y 4 + 25) 4 is an inverse integrating factor. Computing the Lyapunov constants up to degree 2, we can check that the first 17 linear parts are linearly independent. Then, up to a linear change of coordinates in the parameter (1) (1) space we can write Lk = uk , for k = 1, . . . , 17, and Lk = 0, for k = 18, . . . , 33. Conse- quently, using the same scheme as in the proof of Proposition 4.2, simplifying also with triangularization scheme detailed after Proposition 2.3, and using the Implicit Function Theorem we can restrict our study to see the intersection of 16 homogeneous polynomials of degree 2, (2) Lk = Lk (λ̂), for k = 18, . . . , 33, (26) with λ̂ = (u18 , u19 , . . . , u33 ) choosing three perturbation parameters as zero. We recall that for degree 5 perturbations we have 36 parameters, but here only 33 will be relevant. There are three that they do not play any role. 24 LUIZ F. S. GOUVEIA AND JOAN TORREGROSA The next step is to consider the blow-up uk = zzk , for k = 18, . . . , 32 and u33 = z in (26), writing Lk (λ̂) = z 2 Lk (ẑ) with ẑ = (z18 , z19 , . . . , z32 ). Then, we need to show that this system of 15 equations of degree 2 with respect to 15 variables has at least a transversal intersection real point, ẑ ∗ . Moreover, we should check that L33 (ẑ ∗ ) is nonvanishing. The proof finishes applying Theorem 2.2 to provide the complete unfolding of 33 limit cycles of small-amplitude. The main difference with respect to the proof of Proposition 4.2 is that here we can not obtain the explicit solution in terms of polynomials in one privileged variable. Because of the high number of variables and the size of the coefficients of the polynomials Lk . Numerically, we can get an approximate solution ∗ ∗ ∗ z18 ≈ 0.414467055443, z19 ≈ 0.977703106281, z20 ≈ 0.831273897080, ∗ ∗ ∗ z21 ≈ 10.87232453671, z22 ≈ 0.089114602089, z23 ≈ 5.803007782422, ∗ ∗ ∗ z24 ≈ −13.46886905316, z25 ≈ −2.653100632593, z26 ≈ 0.071920750628, (27) ∗ ∗ ∗ z27 ≈ 1.279070836650, z28 ≈ −0.963042490919, z29 ≈ −7.708796674748, ∗ ∗ ∗ z30 ≈ −0.27853535522, z31 ≈ −7.245147157590, z32 ≈ 2.513953010283. Then, L33 (ẑ ∗ ) = −5.28936073528 and the Jacobian matrix of (L18 , . . . , L32 ) with respect to ẑ at ẑ ∗ is 1.2572040284 · 1014 . We have solved numerically with different number of digits (up to 1000 digits) to ensure the convergence of the numerical solution. The above numerical approximation is shown with only 12 digits. Then we can use the Poincaré–Miranda Theorem, to prove analytically the existence of the point ẑ ∗ . This will be doing with an interval analysis for applying Theorem 2.4. Finally, we need to check that the Jacobian and L33 are nonvanishing at ẑ ∗ , also with another accurate interval analysis. This can be done with a Computer Assisted Proof mechanism using the technical lemmas of Section 8. The first step is to convert the approximate solution (27) to the adequate rational expression: ∗ 8289341108806467679 ∗ 97770310628309363091 z18 ≈ , z19 ≈ , 20000000000000000000 100000000000000000000 ∗ 83127389707517527887 ∗ 10872324536703550641 z20 ≈ , z21 ≈ , 100000000000000000000 1000000000000000000 ∗ 22278650522188043761 ∗ 14507519456029022899 z22 ≈ , z23 ≈ , 250000000000000000000 2500000000000000000 ∗ 6734434526579186771 ∗ 26531006325905618217 z24 ≈− , z25 ≈− , 500000000000000000 10000000000000000000 ∗ 35960375313628249487 ∗ 12790708366491310147 z26 ≈ , z27 ≈ , 500000000000000000000 10000000000000000000 ∗ 240760622728448541 ∗ 38543983373719519357 z28 ≈− , z29 ≈− , 250000000000000000 5000000000000000000 ∗ 1740845970112931499 ∗ 14490294315170358141 z30 ≈− , z31 ≈− , 6250000000000000000 2000000000000000000 ∗ 25139530102793502293 z32 ≈ . 10000000000000000000 Then, we consider an affine change of parameters such that the linear part of f = (L18 , . . . , L32 ) will be the new variables. Then, the Jacobian matrix at this point will be near the identity. Next, we apply Theorem 2.4 with n = 15, c = 10−9 to f = (L18 , . . . , L32 ). The conditions about the sign of the components of f on the faces Si± are obtained from Lemmas 8.1 and 8.2. The existence of ẑ ∗ is guaranteed because LOWER BOUNDS FOR THE LOCAL CYCLICITY OF CENTERS 25 fi (Si− ) ⊂ [−2.01 · 10−9 , −1.99 · 10−9 ] and fi (Si− ) ⊂ [1.99 · 10−9 , 2.01 · 10−9 ]. We notice that the numerator and denominators of Li , for i = 18, . . . , 32 are integer numbers hav- ing between 138 to 152 digits each. Moreover, L33 ∈ [−5.482536, −5.0966185] and its numerator and denominator have more than 120 digits each, so we have L33 < 0. Fi- nally, we must show that the determinant of the Jacobian matrix, Jf (ẑ), of f with respect to ẑ does not vanish at ẑ ∗ . This determinant, as it is a 15 × 15 matrix, needs a very high computational cost. Alternatively, we can use the Gershgorin result, Theo- rem 2.5, to show that its eigenvalues are in a ball centered at 1 with radius 10−3 . That is |λ − 1| < 10−3 . Calling Ji,j the i, j-element of Jf (ẑ ∗ ), and using again Lemmas 8.1 and 8.2, we can check that Ji,i ∈ [0.999919, 1.000009] and, for i 6= j, Ji,j ∈ [0.000018, 0.00212] or Ji,j ∈ [−0.0002926, −5.7414588 × 10−7 ] when Ji,j are positive or negative, respectively. Therefore, all eigenvalues are nonzero and, consequently, also the determinant of Jf (ẑ ∗ ). Then, the result follows.  Proposition 6.4. Consider the perturbed system (4) of degree n = 7 with the center (ẋ, ẏ) = (Pc (x, y), Qc (x, y)),     8   ẋ = (1 − x) − y + (2x − y)(24x 5 − 12x4 y − 32x3 y 2 + 12x2 y 3 − 42xy 4 − 5y 5 ) , 45   16    ẏ = (1 − x) 2x − y(2x − y)(28x4 + 66x3 y + 6x2 y 2 + 19xy 3 + 6y 4 ) . 45 (28) There are perturbation parameters λ such that at least 61 limit cycles of small-amplitude bifurcate from the origin. Proof of Proposition 6.4. A similar system as the unperturbed one, without the straight line of equilibria 1 − x = 0, appears in [19] in the study of local cyclicity for homogeneous nonlinearities perturbation. Giné proposes to start with a quartic system as (23) having a center at the origin because it has an integrating factor. Then, we should change to polar coordinates and transform variable r to a new radial variable R3/5 . With these changes the new system has degree 6 and we get system (28) after adding a straight line of equilibria. So, it has a center at the origin. Moreover, it has the next inverse integrating factor V (x, y) =(2x2 + y 2 )−2/3 (128x5 − 192x4 y + 160x3 y 2 − 112x2 y 3 + 48xy 4 − 8y 5 − 9)× (2048x7 y 3 − 3072x6 y 4 + 3584x5 y 5 − 3328x4 y 6 + 2048x3 y 7 − 1024x2 y 8 + 384xy 9 − 64y 10 + 864x3 y 2 + 432xy 4 − 27). Using only linear developments of the Lyapunov constants, see Theorem 2.1, we get only 58 limit cycles of small-amplitude because there exists a linear change of variables in (1) (1) the parameter space such that Lk = uk , for k = 1, . . . , 58 and Lk = 0, for k = 59, 60, 61. Computing up to degree two of Taylor approximations of the Lyapunov constants and doing the simplifications as in the previous proofs, we can remove uk for k = 1, . . . , 58 (2) from Lk , for k = 59, 60, 61. Vanishing the nonrelevant parameters and doing the adequate simple blow-up u59 = z z1 , u60 = z z2 , and u61 = z, we get (2) (2) (2) L59 = A59 z 2 L1 (z1 , z2 ), L60 = A60 z 2 L2 (z1 , z2 ), L61 = A61 z 2 L3 (z1 , z2 ), where Lk are polynomials of degree 2 and Ak rational nonvanishing numbers. These polynomials have rational coefficients with numerators and denominators of around 1900 digits each. Approximately, they write as L1 ≈ z12 + 77.576637z1 z2 + 22.493284z22 + 107.76288z1 + 1038.9032z2 + 1265.0912, L2 ≈ z12 − 2.6270001z1 z2 + 0.27770877z22 + 25.446941z1 − 35.950489z2 + 160.06265, L3 ≈ z12 + 3.4484543z1 z2 + 2.9923181z22 + 32.128183z1 + 44.721607z2 + 248.24137. 26 LUIZ F. S. GOUVEIA AND JOAN TORREGROSA The last step is to show that there exists at least a transversal real solution, z ∗ = (z1∗ , z2∗ ), of {L1 = 0, L2 = 0} such that L3 (z ∗ ) is nonvanishing. Then, Theorem 2.2 applies and the proof follows. With an algebraic manipulator we can find the solution of SL = {L1 = 0, L2 = 0}. It writes as (z1∗ , z2∗ ) = (p3 (α), α) where p3 (α) is a polynomial of degree 3 with rational coefficients and α is a real root of a given polynomial of degree 4, p4 (α). The polynomials p3 and p4 have rational coefficients with numerators and denominators of around 6000 and 4000 digits each, respectively. Approximately they write as p4 (α) ≈ α4 + 0.87600811 α3 − 1.97816765 α2 − 3.22558688 α − 1.29759793, p3 (α) ≈ 3.43467644 α3 − 0.51633002 α2 − 6.541427002 α − 17.7666993. As the polynomial p4 has only 2 real roots, the system SL has only two real solutions ∗ ∗ (z1a , z2a ) ≈ (−14.693346428044632240, −0.85248929481003427092), ∗ ∗ (z1b , z2b ) ≈ (−13.701549420630826548, 1.6907716352120896856). As a function of α, we can find explicitly the values det J(L1 , L2 )(z1∗ ,z2∗ ) and L3 (z1∗ , z2∗ ). They are also polynomials of degree 3 in α with rational coefficients that approximately write as det J(L1 , L2 )(z1∗ ,z2∗ ) ≈ 611.007749 α3 + 530.704211 α2 − 846.328099 α − 1100.76580, L3 (z1∗ , z2∗ ) ≈ 1.53111109 α3 + 0.80160181 α2 − 3.66450101 α − 3.44764247. It can be seen that all the polynomials of degree 3, p3 , det J(L1 , L2 )(z1∗ ,z2∗ ) , and L3 (z1∗ , z2∗ ) have no common zeros with p4 because their respective resultants, with respect to α, are nonzero rational numbers. This proves the transversality and that the last Lyapunov constant, L61 , is nonvanishing. In Figure 1, we have drawn the zero level curves of the polynomials Lk , for k = 1, 2, 3, in a neighborhood of the intersection points. Graphically, the transversality is also clear.  Figure 1. The level curves L1 = 0, L2 = 0, and L3 = 0 in red, green, and blue, respectively. The middle and right pictures are the corresponding zooms near the intersection points 7. Final comments Taking a look at all analyzed systems, it is clear that we need new good examples to get higher lower bounds for the local cyclicity. The main difficulty is to know how to get them to ensure that only with developments of first-order it is enough to get the value originally conjectured by Giné ([19, 20]) and recently updated in [21]. In the language of LOWER BOUNDS FOR THE LOCAL CYCLICITY OF CENTERS 27 Żola̧dek, see [32], this is equivalent to find systems with maximal codimension. It should also be noted that the solution of the center problem is still open, even for polynomial vector fields of degree 3. It is also clear that with this mechanism we will never provide upper bounds. We notice that the importance of Christopher work in [12] is that he pointed out that the computation of high-degree Taylor developments of the Lyapunov constants near a fixed center can be done without knowing their explicit and complete expressions. In no case have we calculated, because of the difficulties for their sizes, the constants first and then we have computed the Taylor developments. This fact has been crucial to perform all the computations made in this work and allow us to go further in determining the best lower bounds for M (n) for lower degrees n. In particular, to design our parallelization algorithm. For studying this local problem, the parallelization mechanism has been really a good tool. It has two computational advantages, the first is the decreasing of the total com- putation time, the second is the decreasing of memory necessities. Because partial com- putations require less time and less memory. Among these advantages, the difficulties now do not depend on the computation mechanisms. They are the size of the objects of higher developments, the knowledge of the local intersection of the varieties, and the high number of variables. Finally, numerical computations are also not easy. Because as the degrees are quite high, to have small approximation errors we need to work with very high precision. 8. Accurate interval analysis Next two technical results will help us to find upper and lower bounds for a polynomial of n variables in a n-dimensional cube. Their proofs can be found in [15]. Lemma 8.1 ([15]). Consider h > 0, p > 0, q real numbers such that p ∈ [p, p], with pp > 0, and q ∈ [q, q], with qq > 0. (i) Then, σ l (q, p) ≤ qp ≤ σ r (q, p), where  l qp, if q > 0, σ (q, p) =  qp, if q < 0, r qp, if q > 0, σ (q, p) = qp, if q < 0. (ii) If uj ∈ [−h, h], for j = 1, . . . , n, and denoting ui = ui11 · · · uinn , for i = (i1 , . . . , in ) 6= 0, we have X l (q, ui ) ≤ qui ≤ X r (q, ui ), where   0, if q > 0 and ik even for all k = 1, . . . , n, X l (q, ui ) = −qhi1 +···+in , if q > 0 and ik odd for some k = 1, . . . , n,  qhi1 +···+in , if q < 0, and   −qhi1 +···+in , if q > 0 and ik even for all k = 1, . . . , n, r i X (q, u ) = 0, if q < 0 and ik odd for some k = 1, . . . , n,  qhi1 +···+in , if q < 0. Furthermore, X l (q, 1) = q and X r (q, 1) = q. 28 LUIZ F. S. GOUVEIA AND JOAN TORREGROSA Lemma 8.2 ([15]). Let h > 0 and pj be positive nonrational numbers such that pj ∈ [pj , pj ] with pj , pj rational numbers satisfying pj , pj > 0, for j = 1, . . . , m. Consider the polynomial ! M X Xm U(u1 , . . . , un ) = Uj,i pj ui i1 +···+in =0 j=1 i with u = ui11 · · · uinn , for i = (i1 , . . . , in ), and Uj,i rational numbers. Then m X Uil ≤ Uj,i pj ≤ Uir , j=1 Pm Pm with Uil = j=1 Uj,i σ l (Uj,i , pj ) and Uir = j=1 Uj,i σ r (Uj,i , pj ). Moreover, if uj ∈ [−h, h], for j = 1, . . . , n, and Uil > Uir then M X M X U= X l (Uil , ui ) ≤ U(u1 , . . . , un ) ≤ X r (Uir , ui ) = U. i1 +···+in =0 i1 +···+in =0 9. Acknowledgments We want to thank Armengol Gasull and Francesc Mañosas for the interesting discus- sions that have allowed us to understand and develop the results of this work. We also thank Fernando Sanz and Maria Aparecida Soares Ruas for guiding us in the interesting geometry of the local intersection of varieties. Finally, we thank the referees for their indications to improve the clarity of the exposition of our results. This work has been realized thanks to the Catalan AGAUR 2017SGR1617 grant, the Spanish Ministerio de Ciencia, Innovación y Universidades - Agencia estatal de inves- tigación MTM2016-77278-P (FEDER) grant, the European Community H2020-MSCA- RISE-2017-777911 grant and the Brazilian CNPq 200484/2015-0 and FAPESP 2020/04717- 0 grants. References [1] A. A. Andronov, E. A. Leontovich, I. I. Gordon, and A. G. Maı̆er. Theory of bifurcations of dynamic systems on a plane. Halsted Press [A division of John Wiley & Sons], New York-Toronto, Ont.; Israel Program for Scientific Translations, Jerusalem-London, 1973. Translated from the Russian. [2] V. I. Arnol’d. Loss of stability of self-oscillations close to resonance and versal deformations of equivariant vector fields. Functional Analysis and Its Applications, 11(2):85–92, Apr 1977. [3] V. I. Arnol’d, S. M. Guseı̆n-Zade, and A. N. Varchenko. Singularities of differentiable maps. Vol. I, volume 82 of Monographs in Mathematics. Birkhäuser Boston, Inc., Boston, MA, 1985. The classifi- cation of critical points, caustics and wave fronts, Translated from the Russian by Ian Porteous and Mark Reynolds. [4] J. L. R. Bastos, C. A. Buzzi, and J. Torregrosa. Factor reversible systems with applications to planar centers. Preprint, May 2020. [5] N. N. Bautin. On the number of limit cycles which appear with the variation of coefficients from an equilibrium position of focus or center type. American Math. Soc. Translation, 1954(100):19, 1954. [6] T. Blows and N. Lloyd. The number of limit cycles of certain polynomial differential equations. Proc. Roy. Soc. Edinburgh Sect. A, 98(3-4):215–239, 1984. [7] Y. L. Bondar and A. P. Sadovskiı̆. On a theorem of Zoladek. Differ. Uravn., 44(2):263–265, 287, 2008. [8] E. Casas-Alvero. Singularities of plane curves, volume 276 of London Mathematical Society Lecture Note Series. Cambridge University Press, Cambridge, 2000. [9] C. Chicone. Ordinary differential equations with applications, volume 34 of Texts in Applied Mathe- matics. Springer, New York, second edition, 2006. [10] C. Chicone and M. Jacobs. Bifurcation of critical periods for plane vector fields. Trans. Amer. Math. Soc., 312(2):433–486, 1989. LOWER BOUNDS FOR THE LOCAL CYCLICITY OF CENTERS 29 [11] S.-N. Chow, C. Li, and D. Wang. Normal forms and bifurcation of planar vector fields. Cambridge University Press, Cambridge, 2008. Reprint of the 1994 original. [12] C. Christopher. Estimating limit cycle bifurcations from centers. In Differential equations with sym- bolic computation, Trends Math., pages 23–35. Birkhäuser, Basel, 2005. [13] A. Cima, A. Gasull, V. Mañosa, and F. Mañosas. Algebraic properties of the Liapunov and period constants. Rocky Mountain J. Math., 27(2):471–501, 1997. [14] A. Cima, A. Gasull, and F. Mañosas. A note on the Lyapunov and period constants. Qual. Theory Dyn. Syst., 19(33):1–13, 2020. [15] L. P. da Cruz, D. D. Novaes, and J. Torregrosa. New lower bound for the Hilbert number in piecewise quadratic differential systems. J. Differential Equations, 266(7):4170–4203, 2019. [16] F. Dumortier, J. Llibre, and J. C. Artés. Qualitative theory of planar differential systems. Universi- text. Springer-Verlag, Berlin, 2006. [17] S. Gerschgorin. Über die abgrenzung der eigenwerte einer matrix. Izv. Akad. Nauk SSSR Otd. Fiz.- Mat., page 749–754, 1931. [18] J. Giné. The center problem for a linear center perturbed by homogeneous polynomials. Acta Math. Sin. (Engl. Ser.), 22(6):1613–1620, 2006. [19] J. Giné. Higher order limit cycle bifurcations from non-degenerate centers. Appl. Math. Comput., 218(17):8853–8860, 2012. [20] J. Giné. Limit cycle bifurcations from a non-degenerate center. Appl. Math. Comput., 218(9):4703– 4709, 2012. [21] J. Giné, L. F. Gouveia, and J. Torregrosa. Lower bounds for the local cyclicity for families of centers. Preprint, March 2020. [22] G. H. Golub and C. F. Van Loan. Matrix Computations (3rd Ed.). Johns Hopkins University Press, Baltimore, MD, USA, 1996. [23] M. Han. Liapunov constants and Hopf cyclicity of Liénard systems. Ann. Differential Equations, 15(2):113–126, 1999. [24] W. Kulpa. The Poincaré–Miranda Theorem. The American Mathematical Monthly, 104(6):545–550, 1997. [25] H. Liang and J. Torregrosa. Parallelization of the Lyapunov constants and cyclicity for centers of planar polynomial vector fields. J. Differential Equations, 259(11):6494–6509, 2015. [26] Maplesoft, a division of Waterloo Maple Inc. Maple 2018. http://www.maplesoft.com/, Waterloo, Ontario. [27] H. Movasati. Center conditions: rigidity of logarithmic differential equations. J. Differential Equa- tions, 197(1):197–217, 2004. [28] R. Roussarie. Bifurcation of planar vector fields and Hilbert’s sixteenth problem, volume 164 of Progress in Mathematics. Birkhäuser Verlag, Basel, 1998. [29] K. S. Sibirskiı̆. On the number of limit cycles arising from a singular point of focus or center type. Dokl. Akad. Nauk SSSR, 161:304–307, 1965. [30] Y. Tian and P. Yu. Bifurcation of small limit cycles in cubic integrable systems using higher-order analysis. J. Differential Equations, 264(9):5950–5976, 2018. [31] C. T. C. Wall. Singular points of plane curves, volume 63 of London Mathematical Society Student Texts. Cambridge University Press, Cambridge, 2004. [32] H. Żola̧dek. The classification of reversible cubic systems with center. Topol. Methods Nonlinear Anal., 4(1):79–136, 1994. [33] H. Żola̧dek. On certain generalization of the Bautin’s Theorem. Nonlinearity, 7(1):233–279, 1994. [34] H. Żola̧dek. Eleven small limit cycles in a cubic vector field. Nonlinearity, 8(5):843–860, 1995. [35] H. Żola̧dek. Remarks on: “The classification of reversible cubic systems with center” [Topol. Methods Nonlinear Anal. 4 (1994), no. 1, 79–136]. Topol. Methods Nonlinear Anal., 8(2):335–342 (1997), 1996. [36] H. Żola̧dek. The CD45 case revisited. In Mathematical Sciences with Multidisciplinary Applications, pages 595–625, Cham, 2016. Springer International Publishing. Departament de Matemàtiques, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Catalonia, Spain 30 LUIZ F. S. GOUVEIA AND JOAN TORREGROSA Departamento de Matemática, Universidade Estadual Paulista, 15054-000 São José do Rio Preto, Brazil Email address: fernandosg@mat.uab.cat, fernando.gouveia@unesp.br Departament de Matemàtiques, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Catalonia, Spain Email address: torre@mat.uab.cat