This is a preprint of: “Local cyclicity in low degree planar piecewise polynomial vector fields”, Luiz Fernando da Silva Gouveia, Joan Torregrosa, Nonlinear Anal. Real World Appl., vol. 60, 103278, 2021. DOI: [10.1016/j.nonrwa.2020.103278] LOCAL CYCLICITY IN LOW DEGREE PLANAR PIECEWISE POLYNOMIAL VECTOR FIELDS LUIZ F. S. GOUVEIA AND JOAN TORREGROSA Abstract. In this work, we are interested in isolated crossing periodic orbits in planar piecewise polynomial vector fields defined in two zones separated by a straight line. In particular, in the number of limit cycles of small amplitude. They are all nested and surrounding one equilibrium point or a sliding segment. We provide lower bounds for the local cyclicity for planar piecewise polynomial systems, Mpc (n), with degrees 2, 3, 4, and 5. More concretely, Mpc (2) ≥ 13, Mpc (3) ≥ 26, Mpc (4) ≥ 40, and Mpc (5) ≥ 58. The computations use parallelization algorithms. 1. Introduction The study of piecewise or nonsmooth differential systems was started by the school of Andronov, see for example [3]. Many problems of engineering can be modeled by this class of systems, see [1]. Recently, they also appear modeling different situations in physics and biology, see [5]. One of the most studied situations in the plane is given by two vector fields defined in two half-planes separated by a straight line. As in the case of the classical qualitative theory of polynomial systems, the study of the number and location of the isolated periodic orbits, also called limit cycles, have received special attention. See for example [7, 14, 20, 25, 28, 31]. In particular, it can be seen as an extension of the 16th-Hilbert problem for planar piecewise polynomial vector fields. Ilyashenko, in [27], presented an updated summary of the status of this problem, proposed by Hilbert more than one hundred years ago. In this work, our main interest is the study of the number of limit cycles bifurcating from the origin in the class of piecewise polynomial differential equations defined in two zones separated by a straight line. In particular, the ones that write as ( (ẋ, ẏ) = (P + (x, y, λ+ ), Q+ (x, y, λ+ )), when y ≥ 0, (1) (ẋ, ẏ) = (P − (x, y, λ− ), Q− (x, y, λ− )), when y < 0, where P ± (x, y, λ± ) and Q± (x, y, λ± ) are real polynomials of degree n in (x, y) and all the 2 coefficients of their monomials define the parameter vector λ = (λ+ , λ− ) ∈ R2n +6n+4 . The straight line Σ = {y = 0} divides the plane in two half-planes Σ± = {(x, y) : ±y > 0} and the trajectories on Σ are defined following the Filippov convention, see [16]. The so-called crossing limit cycles are the ones that, when they pass through the separation line Σ, both vector fields point out in the same direction. For polynomial vector fields of degree n, as usual, we denote by M (n) the maximum number of limit cycles bifurcating from a monodromic singular point and by H(n) the maximum number of limit cycles. For piecewise polynomial vector fields, also of degree n, we denote respectively both numbers by Mpc (n) and Hpc (n). The upper index c means crossing limit cycles. Clearly, M (n) ≤ H(n) and Mpc (n) ≤ Hpc (n). Clearly, linear 2010 Mathematics Subject Classification. Primary 34C07, 34C23, 37C27. Key words and phrases. piecewise vector field, piecewise center cyclicity, Lyapunov quantities. 1 2 LUIZ F. S. GOUVEIA AND J. TORREGROSA systems have no limit cycles, then H(1) = M (1) = 0. While they appear in piecewise linear differential systems. For the class of piecewise linear differential systems defined in two zones separated by a straight line, Huan and Yang in 2012 ([26]) were the first showing the numerical evidence that Hpc (1) ≥ 3 The first analytic proof appeared the same year in [31]. In 2013, using the averaging technique of high order, this lower bound was reobtained in [7]. The same number was confirmed by Freire et al in 2014 ([18]), where a detailed study of the way that the crossing limit cycles appear is done using the full return map. More concretely, two appear near the origin and the other one far from it. But the three limit cycles appear nested and surrounding one sliding segment. In fact, the two limit cycles of small amplitude appearing from an equilibrium point provide the lower bound Mpc (1) ≥ 2. This value can be proved with the results in [17] and we will also show it in the next section. The first proof of the existence of four limit cycles in quadratic vector fields was done in 1980 ([36]). Then H(2) ≥ 4. But only 3 can bifurcate from the origin, that is, M (2) = 3. This fact was proved by Bautin in 1954 ([4]). For piecewise quadratic systems it is not proved yet which will be that local maximum number. Moreover, there are few works providing good lower bounds. Using averaging theory of order five, and perturbing the linear center, Llibre and Tang in [32] proved that Hpc (2) ≥ 8. Recently, da Cruz and et al. in [15] provide a better lower bound, Hpc (2) ≥ 16. These limit cycles appear using also averaging method up to order 2 and perturbing some quadratic isochronous systems. The new lower bound is higher than what it can be expected a priori. Its value is more than the double (because we have two vector fields) of the maximal value aforementioned for quadratic vector fields. In cubic and quartic systems, up to our knowledge, the best known lower bound for the number of limit cycles is H(3) ≥ 13 and H(4) ≥ 28, see [29] and [33], respectively. But for the local cyclicity, the best results, M (3) ≥ 12 and M (4) ≥ 21, are given in [21]. In piecewise polynomial vector fields there are no so many results studying Hpc (3) nor the local cyclicity problem Mpc (3). In the very recent works [22, 24] the reader can find the best-known values for these numbers. In the first, it is proved that Hpc (3) ≥ Mpc (3) ≥ 24. Previously, in the second, studying the phenomenon of simultaneous bifurcation, it was shown that there are 18 limit cycles in a configuration of two nests of 9 each. The main result of this paper deals with new lower bounds for the local cyclicity of piecewise polynomial vector fields of degrees three, four, and five. Theorem 1.1. The local cyclicity for piecewise polynomial vector fields of degrees n = 3, 4, and 5 is Mpc (3) ≥ 26, Mpc (4) ≥ 40, and Mpc (5) ≥ 58, respectively. In particular, Hpc (3) ≥ Mpc (3) ≥ 26, Hp (4) ≥ Mpc (4) ≥ 40 and Hp (5) ≥ Mpc (5) ≥ 58. But, of course, these are not so relevant because using the known values for H(3) ≥ 13 and H(4) ≥ 28, we can provide directly better lower bounds: Hp (3) ≥ 26 and Hp (4) ≥ 56. Because, as our piecewise polynomial vector fields are constructed from 2 polynomial vector fields the most direct lower bound for Hp (n) is, at least, the double than the best known value for H(n). Moreover, we also provide a quadratic system exhibiting at least 13 limit cycles of small amplitude. Furthermore, the limit cycles are all of crossing type and in only one nest, surrounding the same equilibrium point. In this paper, we study the limit cycles of small amplitude bifurcating from an equilibrium point of center-focus type. This bifurcation mechanism is known as degenerate Hopf bifurcation and is based on the precise study of the return map defined in a cross section containing the origin. More concretely, computing the first coefficients of the Taylor series of the return map at the LIMIT CYCLES IN POLYNOMIAL PIECEWISE SYSTEMS 3 origin. These coefficients allow us to define the (generalized) Lyapunov constants. We provide the explicit definition of them in the next section. The key point is the study of polynomial perturbations of centers using high-order Taylor series of the Lyapunov constants with respect to the perturbation parameters. This is based on the Implicit Function Theorem as was firstly used by Chicone and Jacobs in [9] for a similar prob- lem but only with a first-order analysis. Our work is the piecewise extension of the one done by Christopher in [11] and it can be considered as the continuation of [22], where only a first-order study was presented for a piecewise cubic system, proving that Mpc (3) ≥ 24. However, as the computations are quite hard, we have implemented the same parallelization mechanism introduced in [23, 30]. The paper is structured as follows. In Section 2 we present how the Lyapunov con- stants can be computed and how to use them to study degenerate Hopf and pseudo-Hopf bifurcations. In particular, we prove the existence of generic unfolding of 2n + 1 limit cycles of small amplitude bifurcating from a weak-focus of such order. We detail also the differences in the order between a weak-focus of analytic and piecewise analytic vector fields. Section 3 is devoted to show how the proposed technique works to find good lower bounds for the local cyclicity in quadratic, cubic, quartic, and quintic piecewise poly- nomial vector fields. From them we get our main result. We perturb quadratic, cubic, quartic, and quintic centers with piecewise polynomial systems, of the corresponding degrees 2, 3, 4, and 5, defined in two zones separated by a straight line. We finish this work, in Section 4, detailing the computational difficulties to get further in the study of number of limit cycles of small amplitude in higher degree piecewise polynomial vector fields. 2. Degenerated Hopf bifurcations in piecewise vector fields Let us introduce the concepts of crossing and sliding segments which are necessary for the study of the dynamic of a piecewise differential system. The existence of an sliding segment is used in this work only to provide one extra limit cycle. Mainly, we will use perturbations that maintain the (isolated) equilibrium point on the separation line. Let be a planar differential system defined as ( Z + (x, y), when f (x, y) ≥ 0, Z(x, y) = (2) Z − (x, y), when f (x, y) < 0, with f : R2 → R a C 1 function such that 0 is a regular value. The discontinuity curve is given by Σ = f −1 (0), and Z ± = (X ± , Y ± ). Following the notation introduced by Filippov in [16], in Σ when both vector fields meet, we can have on Σ three different behaviors, that are of crossing, escaping and sliding type. We will denote the respective segments by ΣC , ΣE and ΣS . Given a point p ∈ Σ, we say that p ∈ ΣC if, and − + S fS(p) · Z f (p) > only if Z 0 where Z ± f (p) = h∇f (p), Z ± (p)i. Consequently, we have p∈Σ E Σ if, only if Z f (p) · Z − f (p) < 0. Figure 1 illustrates how is the vector field + near these three regions when f (x, y) = y. Figure 1. Crossing, escaping, and sliding segments 4 LUIZ F. S. GOUVEIA AND J. TORREGROSA As in this work we are interested in piecewise differential systems of type (1), we will consider only piecewise analytic vector fields (2) with f (x, y) = y. That is, the ones defined by ( Z + (x, y), when y ≥ 0, Z(x, y) = (3) Z − (x, y), when y < 0, being Z ± analytic vector fields. Now, we detail the algorithm that we have implemented to compute the coefficients of the return map, Π(ρ), near the origin, when it is of monodromic type, in a planar piece- wise vector field. As, from (3), we have two vector fields, we will have two half-return maps, Π± (ρ), and the global one can be defined by composition, Π(ρ) = Π− (Π+ (ρ)), and equivalently, the periodic orbits can be obtained from the zeros of the displacement map, Π(ρ) − ρ. But, for simplicity, instead of this map we will compute the equivalent one ∆(ρ) = (Π− )−1 (ρ) − Π+ (ρ). Here, the coefficients of the Taylor series of the differ- ence map ∆ at the origin are also called Lyapunov constants for piecewise polynomial vector fields. See Figure 2. Π+ (ρ) Π− (Π+ (ρ)) Π+ (ρ) ρ (Π− )−1 (ρ) ρ Figure 2. Crossing, escaping, and sliding segments As in this work we are dealing with polynomial vector fields having nondegenerate centers and the polynomial perturbations will be different on y > 0 and on y < 0, we write our system (1), or the vector field (3), as  ẋ = −y + P ± (x, y, λ± ), (4) ẏ = x + Q± (x, y, λ± ), being P ± , Q± polynomials without constant or linear terms. That is, we restrict our analysis to the case that, for every λ± , the trace of the Jacobian matrix of (4) at the origin is zero and there are no sliding or scaping segments. Writing (4) in the usual polar coordinates we have  X∞  dr =   Ri+ (θ, λ+ )ri , θ ∈ [0, π],  dθ i=2 X∞ (5)   dr − − i   = Ri (θ, λ )r , θ ∈ [π, 2π], dθ i=2 where Ri± (θ, λ± ) are trigonometrical polynomials in θ. In a neighborhood of the origin, the Taylor series of the solution r± (θ, ρ) of (5), such that r± (0, ρ) = ρ, can be written as X∞ ± ± r (θ, ρ, λ ) = ρ + ri± (θ, λ± )ρi , (6) i=2 LIMIT CYCLES IN POLYNOMIAL PIECEWISE SYSTEMS 5 with ri± (0) = 0 for i ≥ 2. As our piecewise systems are defined separated by the straight line {y = 0}, the half-return maps close to the origin are given by ∞ X + Π (ρ) = ρ + ri+ (π, λ+ )ρi , i=2 ∞ X (Π− )−1 (ρ) = ρ + ri− (−π, λ− )ρi . i=2 Therefore, as we have mentioned above, the Taylor series at ρ = 0 of the difference map can be written as X∞ ∞ − −1 + − − + +  i X ∆(ρ) = (Π ) (ρ) − Π (ρ) = ri (−π, λ − ri (π, λ )) ρ = Li ρ i . (7) i=2 i=2 As the systems (4) are analytic, the half-return maps Π± and consequently the return map Π and the difference ∆ as well. The coefficients Li in (7) are known as the generalized Lyapunov constants associated to system (4). Consequently, the first nonvanishing Lk provides the stability of the origin and we will say that the origin is a generalized weak-focus of order k and so ∆(ρ) = Lk ρk + O(ρk+1 ) in the piecewise context. For simplicity, when the piecewise context perturbation is clear, we will avoid the word “generalized”. The above procedure follows closely the classical Lyapunov algorithm described, for example, in [2] for analytic vector fields. For them, we only have one differential equation and so ri+ = ri− , λ+ = λ− = λ, and the first nonvanishing coefficient has always an odd subscript. Therefore, the difference map in the analytic context is ∆(ρ) = Lk ρ2k+1 + O(ρ2k+2 ) and we say that the origin is a weak-focus of order k. This property among others are proved in [2]. We recall that these constants appear solving the analytic center-focus problem for nondegenerate centers. In [12] is shown that the Lyapunov constants Lk , after dividing by π, are polynomials in the parameters λ with rational coefficients. Clearly, they are defined when the previous vanish. The former property is also valid in the piecewise context, but for the latter, it is a little different. The constants are polynomials in π with rational coefficients. We emphasize that the main difference in the Taylor series (7) between the analytic study versus the piecewise one is the fact that in the first, the Lyapunov constants with even indices are zero while in the second not. This difference will be shown in the examples of next sections. In the classical Hopf bifurcation for analytic vector fields only one limit cycle bi- furcates from a weak-focus of order 1 located at the origin, see [2]. In this case, the trace of the Jacobian matrix of the vector field at the equilibrium point is zero, the determinant is positive, and the Taylor development of the displacement map writes as ∆(ρ) = L1 ρ3 +O(ρ4 ). The limit cycle bifurcates from the origin adding a new parameter that changes the stability of the equilibrium point. This can be done with a parame- ter that controls the trace of such Jacobian matrix because in the Taylor series of the displacement map, the coefficient in ρ appears. The next result, stated and proved in [22], is the generalization of this property to piecewise analytic vector fields, where two crossing limit cycles of small amplitude bifurcate from the origin. We have also added here its proof by completeness. It follows from the study of the return map near the origin in a similar way as the procedure described in [17]. In this case, the displacement map writes as ∆(ρ) = L1 ρ2 + O(ρ3 ) and the two extra crossing limit cycles appear when two new essential parameters are added. One is again the trace parameter and another is the constant term. By convenience, the first is added in the upper differential 6 LUIZ F. S. GOUVEIA AND J. TORREGROSA system while the second in the lower one. As it can be seen in the proof, with these new essential parameters, we can control the coefficient of ρ and the constant term of the displacement map of the perturbed piecewise differential system. Proposition 2.1 ([22]). Consider the perturbed system    P∞  P ∞  2  ẋ = −(1 + c )y + + k ` ak` x y ,   ẋ = −y + a− k ` k` x y , k+`=2 k+`=2 P ∞ P∞ (8)   b+ k `   b− k `  ẏ = x + 2cy + k` x y ,  ẏ = d + x + k` x y , k+`=2 k+`=2 − − − for y ≥ 0 and y < 0, respectively. If a+ + + 11 + 2b02 + b20 6= a11 + 2b02 + b20 then there exist c and d small enough such that two crossing limit cycles of small amplitude bifurcate from the origin. Proof. When c = d = 0 the piecewise differential equation in the statement write as (4) with λ± = (a± ± 20 , . . . , b20 , . . .). Writing it in polar coordinates in the form (5), we can get the first term in the sum of (6) as 1 r2± (θ, λ± ) = −(a11 − b02 + b20 ) cos3 θ − (a02 − a20 + b11 ) sin θ cos2 θ 3  − b02 cos θ + (a02 + 2 a20 + b11 ) sin θ + a11 + 2b02 + b20 . Hence, evaluating at ±π we get r2± (±π, λ± ) = 2(a± ± ± 11 + 2b02 + b20 )/3 and, consequently, the first term of the difference map (7) is 2 − − −  L2 = − (a+ + + 11 + 2b02 + b20 ) − (a11 + 2b02 + b20 ) . 3 Therefore, from the condition in the statement, the origin is stable or unstable because the above coefficient is nonvanishing. So, for c, d small enough, computing the return map similarly as in the proof of Proposition 7.3 of [17] we can write, for ρ 6= 0, ∆(ρ) = ∆0 (c, d) + ∆1 (c, d)ρ + ∆2 (c, d)ρ2 + O(ρ3 ), where ∆0 (c, d) = 2d, ∆1 (c, 0) = 1 − eπc , and ∆2 (0, 0) = L2 . As c and d are arbitrary parameters, two crossing limit cycles of small amplitude bifurcate from the origin.  In the above result, when there is no sliding segment (d = 0 and ∆0 = 0), a unique crossing limit cycle of small amplitude bifurcates (for c small enough) from the origin when ∆1 and ∆2 have opposite sign. This bifurcation is described in [13] where the return map of a focus-focus point is studied without the existence of a sliding segment. This bifurcation phenomenon is exactly the same as the classical Hopf bifurcation (see for example [2]) in the nonpiecewise scenario because, after the perturbation due to the parameter c, inside the limit cycle of small amplitude there is only a unique equilibrium point. After this bifurcation and due to the parameter d, a second limit cycle of small amplitude can bifurcate from the equilibrium point. If this is the case, inside the limit cycle of small amplitude, there is only a sliding segment. This second bifurcation mechanism has recently called pseudo-Hopf bifurcation, see more details in [8, 15]. It was previously presented in [17, 18] but without using this name. Although the proof published in [18] was for piecewise linear differential equations, it can be used also for the general case because it does not depend on the higher degree terms of the piecewise perturbation. The described notation for the Lyapunov constants was introduced in [11] in the analytic scenario, where the first nonvanishing coefficient in the return map is labeled as L1 and the trace of the Jacobian matrix is named as L0 to have a unified notation. LIMIT CYCLES IN POLYNOMIAL PIECEWISE SYSTEMS 7 In the piecewise scenario we have also used this unified notation: The first nonvanishing Lyapunov constant is labeled with a subindex 2, that is L2 , and it is only defined when there are no constant terms (∆0 = 0) and the linear part has zero trace and positive determinant (∆1 = 0). As in [11], if necessary and making abuse of notation, we also denote the first two terms of the displacement map of the piecewise perturbed system as L0 and L1 . The main advantage of the above notations and definitions is that the subscript of the first nonvanishing Lyapunov constant determines not only the weak-focus order, it also provides the maximal number of limit cycles of small amplitude that bifurcate from an equilibrium point of monodromic type with finite order. Hence, we can easily present what is called the degenerate Hopf bifurcation, that describes the birth of k limit cycles of small amplitude from a weak-focus of order k. This result is well-known for analytic vector fields, see [2, 35], and we will extend it to piecewise perturbations in Theorem 2.4 and Corollary 2.6. We notice that although the results seem the same, because in both cases k limit cycles will bifurcate, the definitions of weak-focus order are not. From an analytic weak-focus of order k bifurcates k limit cycles with (generic) analytic perturbations while 2k + 1 when the (generic) perturbations are piecewise analytic of type (4). Because the weak-focus order of an analytic differential system depends on which scenario is considered. In qualitative theory of differential equations a very interesting and difficult problem is the study of the maximal number of limit cycles of small amplitude bifurcating from centers. The so-called cyclicity of a center. In this work we are mainly interested in this problem but for nondegenerate centers and using piecewise perturbations. Due to the difficulty to find upper bounds for this number, we will present here new results to provide lower bounds of it. More concretely, to find upper bounds of the number of limit cycles of small amplitude bifurcating from the center itself, the so-called local cyclicity. We will consider only piecewise families of type (4) with perturbations of the form  ẋ = −y + Pc (x, y) + P ± (x, y, λ± ), (9) ẏ = x + Qc (x, y) + Q± (x, y, λ± ), being Pc , Qc , P ± , Q± polynomials without constant or linear terms such that when the perturbation parameters λ = (λ+ , λ− ) vanish, the system has a center at the origin. So P ± (x, y, 0) = Q± (x, y, 0) ≡ 0 and the polynomial system (ẋ, ẏ) = (−y + Pc (x, y), x + Qc (x, y) is a nonpiecewise polynomial nondegenerate center. Clearly, the Lyapunov constants Lk of system (9) depend polynomially on λ and vanish at λ = 0. Therefore, (`) denoting by Lk (λ) the Taylor development of Lk up to order ` at λ = 0, we have that (`) also Lk (0) = 0 for every positive integer `. Before presenting the general results of this section, we will show a simple example of how acts the degenerate Hopf bifurcation perturbing a quadratic center with a specific piecewise quadratic vector field. Moreover, it will be an illustration of how are the usual computations appearing in the other results that we will not show because of the size of the expressions. In Proposition 2.2 we will explain the main mechanism to get a weak-focus of high-order at the origin and how the complete unfolding of limit cycles of small amplitude is obtained. We notice that when the perturbation is polynomial the unfolding may not be complete while it is for analytic perturbations. As the Lyapunov constants can be very difficult to manage for providing such unfoldings, (`) we will show in Proposition 2.3 how the Taylor developments Lk (λ) of the Lyapunov constants near fixed centers can be used to study lower bounds for the local cyclicity. 8 LUIZ F. S. GOUVEIA AND J. TORREGROSA We will see that the number of limit cycles bifurcating from the origin, in family (10), is less when the parameters are big (Proposition 2.2) than when they are small enough (Proposition 2.3). This is the mechanism recovered by Christopher in [11] but originally developed by Chicone and Jacobs in [10]. Proposition 2.2. Consider the piecewise quadratic perturbed system  ẋ = −y − c2 y + x2 − y 2 + a1 x2 + a2 xy, for y ≥ 0, ẏ = x + 2cy + 2xy,  (10) ẋ = −y + x2 − y 2 + a3 xy, for y ≤ 0. ẏ = d + x + 2xy + a4 x2 , Then, there exist parameters a1 , a2 , a3 , a4 , c, and d such that 6 crossing limit cycles of small amplitude bifurcate from the origin. Moreover, when d = 0 the origin is a center if and only if c = a2 = a3 = a4 = 0 or c = a1 = a2 − a3 = a4 = 0. Proof. We start restricting our analysis to the case c = d = 0. After computing and studying the displacement map we will characterize the center families. Then we will prove that there exists a∗ = (a∗1 , a∗2 , a∗3 , a∗4 ) such that the origin is a weak-focus of order 6 for which a generic unfolding inside (10) provides 4 limit cycles. The other 2 limit cycles are provided by Proposition 2.1 choosing adequately small enough values for c and d. System (10), fixing c = d = 0, can be written in polar coordinates as (5) with 1 R2+ (θ) = ((3a1 + 4) cos θ + a2 sin θ + a1 cos 3θ + a2 sin 3θ) , 4 1 R3+ (θ) = 4a1 a2 + 2a1 a2 cos 2θ + (5a21 + 3a22 − 16) sin 2θ 32  − 4a1 a2 cos 4θ + 4a21 sin 4θ − 2a1 a2 cos 6θ + (a21 − a22 ) sin 6θ , 1 R2− (θ) = (4 cos θ + (a3 + a4 ) sin θ + (a3 + a4 ) sin 3θ) , 4 1 R3− (θ) = −16a4 − 16a4 cos 2θ + (3a23 − 2a3 a4 − 5a24 − 16) sin 2θ 32  − 4a4 (a3 + a4 ) sin 4θ − (a23 + 2a3 a4 + a24 ) sin 6θ . The next functions Ri are not written here because of their size. Straightforward computations show that the first coefficients in ρ of the solution (6) are 1  r2+ (θ) = 4a2 − 3a2 cos θ + 3(3a1 + 4) sin θ − a2 cos 3θ + a1 sin 3θ , 12 1 r3+ (θ) = 72a1 a2 θ + 230a21 + 108a2 2 + 432a1 + 144 + −96a22 cos θ 576 + 96(3a1 + 4)a2 sin θ − 3(57a21 − a22 + 128a1 + 48) cos 2θ − 6(11a1 + 16)a2 sin 2θ − 32a22 cos 3θ + 32a1 a2 sin 3θ − 6(9a21 − 2a22 + 8a1 ) cos 4θ − 6(11a1 + 8)a2 sin 4θ  − 5(a21 − a22 ) cos 6θ − 10a1 a2 sin 6θ , 1  r2− (θ) = 4(a3 + a4 ) − 3(a3 + a4 ) cos θ + 12 sin θ − (a3 + a4 ) cos 3θ , 12 1 r3− (θ) = −288a4 θ + 108a23 + 126a3 a4 + 18a24 + 144 − 96(a3 + a4 )2 cos θ 576 + 384(a3 + a4 ) sin θ + 3(a23 + 26a3 a4 + 25a24 − 48) cos 2θ LIMIT CYCLES IN POLYNOMIAL PIECEWISE SYSTEMS 9 − 48(2a3 + 5a4 ) sin 2θ − 32(a3 + a4 )2 cos 3θ + 6(2a23 + 7a3 a4 + 5a24 ) cos 4θ  − 48(a3 + a4 ) sin 4θ + 5(a3 + a4 )2 cos 6θ . Again, because of the size of them, the next functions ri are not explicitly written here. The next step is its evaluation, respectively at ±π, and use the definition (7) to calculate the first Lyapunov constants. Which can be written as L2 = −2(a2 − a3 − a4 )/3, L3 = −π(a1 a2 − 4a2 + 4a3 )/8, L4 = −(14a21 a2 − 12a32 + 30a22 a3 − 18a2 a23 − 30a2 + 30a3 )/45, L5 = −π(a3 − a2 )(422a22 + 25a2 a3 − 840a23 + 2161)/9408, (11) L6 = −(4320a52 − 42180a42 a3 + 108168a32 a23 − 111780a22 a33 + 41472a2 a43 − 139510a21 a3 − 30704a32 + 103396a22 a3 − 119702a2 a23 + 47010a33 + 129480a2 − 129480a3 )/4862025. Here we have used that each one is defined when the previous vanish. Solving the system of algebraic equations {L2 = L3 = L4 = L5 = L6 = 0} we get only the two families of centers in the statement. The first is a center because up and down differential systems are invariant with respect to the change (x, y, t) → (−x, y, −t). More concretely, the solutions and the respective half-return maps are symmetric with respect to x = 0. For the second family, (10) is a nonpiecewise system because up and down ones coincide, we have a Lotka-Volterra center. See the quadratic classification of centers in [38]. Therefore, the center characterization is finished. The weak-foci of order 6 are obtained solving {L2 = L3 = L4 = L5 = 0} and removing the two families of centers previously studied. Straightforward computations show that, in the restricted 4-parameter space, for each simple zero α of ϕ1 (β) = 27β 6 +1344β 4 −34158β 2 +138304, we have that L2 (a∗ ) = L3 (a∗ ) = L4 (a∗ ) = L5 (a∗ ) = 0 and L6 (a∗ ) = −αϕ2 (α)/18234720000 6= 0. Being a∗ = (a∗1 , a∗2 , a∗3 , a∗4 ), a∗1 = (27α4 + 1680α2 − 9678)/3384216, a∗2 = α, a∗3 = −α(27α4 + 1680α2 − 20398)/10720, a∗4 = 3α(9α4 + 560α2 − 3226)/10720. and ϕ2 (β) = 4291351539β 4 +277259653344β 2 −1572537916894. The condition L6 (a∗ ) 6= 0 follows easily because the polynomials ϕ1 and ϕ2 have no common roots and α 6= 0. In fact, the resultant of ϕ1 and ϕ2 with respect to β is a nonzero rational number. The complete unfolding, and so the existence of four limit cycles of small ampli- tude, is guaranteed by the Implicit Function Theorem. In particular, the determi- nant of the Jacobian matrix of (L2 , L3 , L4 , L5 ) with respect to (a1 , a2 , a3 , a4 ) at a∗ is αϕ3 (α)π 2 /2500761600 6= 0 where ϕ3 (β) = 32113869β 4 + 2044570656β 2 − 11609103778. The nonvanishing condition follows because the polynomials ϕ1 and ϕ3 also have no common roots and α 6= 0. In fact the polynomial ϕ1 has exactly 4 simple real roots. Approximately they are ±2.2903752145 and ±3.7579140809.  Proposition 2.3. Consider the piecewise quadratic perturbed system (10). Then, the local cyclicity, considering λ = (a1 , a2 , a3 , a4 , c, d) small enough parameters, is at least 4 using Taylor series up to order 20 with respect to λ at λ = 0. 10 LUIZ F. S. GOUVEIA AND J. TORREGROSA Proof. As the previous proof, we will prove only that 2 limit cycles bifurcate from the origin under the restriction c = d = 0. Finally, applying Proposition 2.1, the statement follows. From the proof of Proposition 2.2 the Lyapunov constants of system (10) are written as (11) and their first-order developments are (1) (1) L2 = −2(a2 − a3 − a4 )/3, L5 = 2161π(a2 − a3 )/9408, (1) (1) L3 = π(a2 − a3 )/2, L6 = 8632(a2 − a3 )/15435. (1) L4 = 2(a2 − a3 )/3, We remark that only the first two are linearly independent. Hence, up to first-order (1) development and with an adequate change of parameters, we can write L2 = u2 and (1) L3 = u3 considering that the parameters are now λ = (a1 , u2 , a3 , u3 ). Clearly, when (1) (1) (1) (1) (1) L2 = L3 = 0 we have that also L4 = L5 = L6 = 0. Using the Implicit Function Theorem, there is a variety (in the parameter space) of weak-foci of order 3 when u2 = 0 and u3 6= 0. Moreover, a limit cycle of small amplitude bifurcates from the origin when u2 , u3 are small enough, 0 < |u2 |  |u3 |, and u2 u3 < 0. If we consider up to second-order developments we have (2) (2) L2 = −2(a2 − a3 − a4 )/3, L5 = 2161π(a2 − a3 )/9408, (2) (2) L3 = π(a2 − a3 )/2 − πa1 a2 /8, L6 = 8632(a2 − a3 )/15435. (2) L4 = 2(a2 − a3 )/3, (2) In this case, the term a1 a2 in L3 is the key point to get an extra limit cycle because higher-order developments provide weak-foci of order 4 but not higher, at least studying up to Taylor series of order 20 in the perturbation parameters. One way to see this property, for each 2 ≤ ` ≤ 20, is using the Implicit Function Theorem for writing (`) (`) (a2 , a3 , a4 ) as functions of new parameters (u2 , u3 , u4 ) such that L2 = u2 , L3 = u3 , (`) (`) and L4 = a1 u4 . The next step is to check that when u2 = u3 = 0 we have L4 = a1 u4 (`) and Lk = a1 u4 M (a1 , u4 ), for k = 5, . . . , 8, with M (0, 0) 6= 0. These computations can be done easily with the help of a Computer Algebra System. Consequently, in a neighborhood of λ = 0, taking u2 = u3 = 0 and a1 u4 6= 0, we have a variety of weak-foci of order 4 passing through the origin of the parameter space that unfolds four limit cycles of small amplitude because of the independence of u2 and u3 from a1 and u4 . In addition, there are no weak-foci of higher-order because when a1 or u4 is zero we have (`) also Lk = 0, for k = 5, . . . , 8.  An alternative way to do the higher-order study presented in the last proof is the simplification mechanism described in [23]. Moreover, we have no computed more Lya- punov constants because, with the center characterization presented in Proposition 2.2, they are not necessary to be used. In addition, although higher-order study could be developed, our interest here is not the study of the local cyclicity of such a fixed qua- dratic system because it has nothing special. It is just a simple example of how the two approaches work. In fact, the results of the following sections improve this local cyclicity value for the quadratic polynomial piecewise family. The following results are the natural extension for piecewise analytic vector fields of the corresponding to analytic vector fields. LIMIT CYCLES IN POLYNOMIAL PIECEWISE SYSTEMS 11 Theorem 2.4. Consider the class of piecewise analytic systems (3) without sliding segment and such that both Z ± have equilibria at the origin. If Z has a (generalized) weak-focus of order k at the origin then, the local cyclicity is at most k − 1. Moreover, there are piecewise analytic perturbations inside the same class (3) without constant terms, such that k − 1 hyperbolic crossing limit cycles of small amplitude bifurcate from the origin. Remark 2.5. As we have commented before, Theorem 2.4 also includes the case when the unperturbed system is analytic and so k is always an odd number (k = 2` + 1). We notice that the piecewise perturbation exhibits k − 1 = 2` limit cycles of small amplitude bifurcating from the origin instead of the analytic one that only bifurcate `. See [35]. Corollary 2.6. Consider the class of piecewise analytic systems (3) and such that both Z ± have equilibria at the origin. If Z has a (generalized) weak-focus of order k at the origin then, at least k hyperbolic crossing limit cycles of small amplitude bifurcate from the origin. The above corollary is a direct consequence of Theorem 2.4 together with the psedo- Hopf bifurcation provided in Theorem 2.1. Because we allow the existence of an sliding segment in the perturbed system having also constant terms. So, we only need to prove the theorem. In fact, the proof follows closely the birth of k limit cycles of small amplitude bifurcating from a weak-focus point of order k when the unperturbed and the perturbed vector field are both analytic. This is clearly described in [2, 35]. We only need to check that the proof can be adapted to our perturbed piecewise analytic vector fields. Proof of Theorem 2.4. We take a piecewise analytic vector field Z as in (3) having a (generalized) weak-focus of order k. Then, it can be written in the form (4) for some fixed values of the parameters λ = (λ+ , λ− ). We can assume that λ = 0, doing a translation in the parameters if necessary. For using the mechanism described at the beginning of this section we write it in polar coordinates as (5) with λ = 0. For simplicity, we will denote it by Z̃0 , just to indicate that the unperturbed system is written in polar coordinates and λ = 0. By hypothesis, the difference map (7) can be written as ∆(ρ) = Lk ρk + · · · , with Lk 6= 0. Considering a general (analytic) perturbation (3) written in polar coordinates (5) and denoted by Z̃λ we have that the difference map of the perturbed piecewise differential system can be written as ∆(ρ, λ) = f1 (λ)ρ + f2 (λ)ρ2 + · · · + fk−1 (λ)ρk−1 + Lk ρk + O(ρk+1 ). Clearly ∆(0, λ) = 0, because we have no sliding segment by hypothesis. We have also indicated the dependence on the parameters λ. Using the Weierstrass Preparation 6= 0, there exist analytic functions f˜k and F such that k Theorem, because ∂∂ k∆ρ (0,0) ∆(ρ, λ) = (f˜1 (λ)ρ + f˜2 (λ)ρ2 + · · · + f˜k−1 (λ)ρk−1 + f˜k (λ)ρk )F (ρ, λ), where F (0, 0) 6= 0 and f˜k (0) = Lk . Then, clearly, the function ∆ can have at most k − 1 nonvanishing zeros as the first statement ensures. Next, we will prove that there exists a precise perturbation such that k − 1 limit cycles bifurcate from the origin. Let us consider now λ = (λ1 , . . . , λk−1 ) and a specific perturbation Z̃λ = Z̃0 + W̃λ being W̃λ the piecewise vector field associated to the 12 LUIZ F. S. GOUVEIA AND J. TORREGROSA piecewise differential equation, in polar coordinates,  dr   = λ1 r + λ2 r2 + · · · + λk−1 rk−1 , θ ∈ [0, π], dθ (12)   dr = 0, θ ∈ [π, 2π]. dθ We observe that the perturbation only affects the system on y > 0. Therefore, for the unperturbed system, λ = 0, we have ∆(ρ, 0) = Lk ρk + O(ρk+1 ) and for the perturbed system, λ ≈ 0, and using the mechanism described at the beginning of the section only the half-return map Π+ defined on y > 0 is necessary to be computed. The half-return map Π− , defined on y < 0, remains unchanged. Finally, we will explain how the values λj should be chosen. We start assuming, rescaling time if necessary, that the origin is an attracting equi- librium point (Lk < 0). So there exists a neighborhood Uk of the origin such that all the orbits in Uk go, in forward time, to the origin. Next we take λj = 0, for j = 1, . . . , k − 2 and λk−1 < 0 small enough. Then, considering a small neighborhood Uk−1 ⊂ Uk , we can compute the positive half-return map following the previously described mecha- nism and, as the negative half-return map does not change, the displacement map is ∆(ρ) = −λk−1 πρk−1 + O(ρk ). So, the origin is a repelling equilibrium point and an stable crossing limit cycle of small amplitude inside Uk−1 bifurcates from the origin as λk−1 changes from 0 to a negative small value. This phenomenon is qualitatively equal to a Hopf bifurcation. Repeating this reasoning for every j up to j = 1 and alternating the signs of λj , we can obtain a new crossing limit cycle of “Hopf type” at each step j. So the sequence of parameters in (12) must be decreasing in absolute value and satisfying λj λj+1 < 0. We remark that in this mechanism, the birth of each limit cycle is controlled by each λj . Then the second part of the statement is proved.  As we have shown in the proof of Proposition 2.2, computing the first-order terms of the Lyapunov constants we can use the Implicit Function Theorem to study the bifur- cation of hyperbolic limit cycles of small amplitude from the origin. Next two results generalize Theorems 2.1 and 3.1 given by Christopher in [11] to piecewise polynomial vector fields. Although the proofs of Christopher follow in this piecewise context, an alternative proof of the second one can be found in [23]. That is, with a direct use of the Implicit Function Theorem or also using previously a specific blowup. An scheme of the proof of this second result can be follow in the proofs of some of the results of next section. Of course, in both cases, we will use at the end again Proposition 2.1. That is, first studying the hyperbolic limit cycles bifurcating from the origin with c = d = 0 in (8) and then adding 2 extra hyperbolic limit cycles varying these two new essential parameters. Theorem 2.7. Consider the perturbed piecewise polynomial differential system of degree n of the form (9), with Pc , Qc without linear or constant terms and the unperturbed vector field, (ẋ, ẏ) = (−y + Pc (x, y), x + Qc (x, y)), has a center at the origin. If the first- order Taylor developments with respect to the perturbation parameters at the origin (1) (1) of the first k − 1 Lyapunov constants, (L2 , . . . , Lk ), associated to (9) are linearly independent, then there exist perturbation parameters (a± ± k` , bk` ) such that system     Pn Pn   (ẋ, ẏ) = −y + Pc (x, y) + + k ` ak` x y , x + Qc (x, y) + + k ` bk` x y for y ≥ 0,  k+`=0 k+`=0  (13)  Pn Pn   (ẋ, ẏ) = −y + Pc (x, y) + − k ` ak` x y , x + Qc (x, y) + − k ` bk` x y for y < 0, k+`=0 k+`=0 LIMIT CYCLES IN POLYNOMIAL PIECEWISE SYSTEMS 13 has at least k crossing limit cycles of small amplitude bifurcating from the origin. Theorem 2.8. Consider the perturbed system of the form (9) with Pc , Qc without linear or constant terms and the unperturbed vector field, (ẋ, ẏ) = (−y +Pc (x, y), x+Qc (x, y)), has a center at the origin such that, after a change of variables in the parameter space if necessary, the first k − 1 Lyapunov constants vanish (L2 = · · · = Lk = 0) and the next l Lyapunov constants write as 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 ) are the new parameter values. 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 perturbation parameters (a± ± k` , bk` ) such that system (13) has at least k + l crossing limit cycles of small amplitude bifurcating from the origin. 3. Lower bounds for the local cyclicity In this section, we illustrate how the degenerated Hopf bifurcation explained in the previous section provides a good mechanism to obtain new lower bounds for the number of limit cycles of small amplitude bifurcating from a monodromic equilibrium point that, without loss of generality, we have located at the origin. That is, to get lower bounds for the local cyclicity Mpc (n). We present the results for low degrees n = 2, 3, and 4. Proposition 3.1. Consider the perturbed system of the form (13) with n = 2 and Pc (x, y) = 18x2 + 8xy − 8y 2 and Qc (x, y) = 4x2 + 14xy − 4y 2 . Then, there exist small enough values of the parameters a± ± k` , bk` such that at least 13 hyperbolic crossing limit cycles of small amplitude bifurcate from the origin. Proof. The origin of the unperturbed system (13) is a Darboux center with the rational first integral, well defined at the origin, (80x3 − 480x2 y + 960xy 2 − 640y 3 + 120xy − 240y 2 − 30y − 1)2 H(x, y) = . (20x2 − 80xy + 80y 2 + 20y + 1)3 First we restrict our analysis to the case that the constant and linear perturba- tion monomials are zero in (13), i.e. a perturbation of type (9). Under this as- sumption and changing to polar coordinates, the algorithm described at the begin- ning of Section 2 provides the coefficients Lk of the Taylor series of the displacement map, (7), with respect to ρ. As we have already mentioned before, the first two co- efficients in ρ, the corresponding to constant and linear ones, vanish. We will see in the following that only L2 , . . . , L13 are necessary to be computed. Moreover, as we are perturbing a center, the Taylor developments, with respect to the parameters − − − − − − λ = (a+ 20 , a11 , a02 , b20 , b11 , b02 , a20 , a11 , a02 , b20 , b11 , b02 ) ∈ R , start with linear terms and, + + + + + 12 after some straightforward computations, we can write the first seven as (1) 2 − − − L2 = − (a+ + 11 − a11 + 2(b02 − b02 ) + b20 − b20 ), + 3 (1) 5π L3 = − (−4a− + − 11 + 6a11 − 13b02 − 10b20 + 7b02 ), − + 4 (1) 8 − − − L4 = − (100(a+ + 02 − a02 ) + 54(a20 − a20 ) − 250(b02 − b02 ) + 15 − + 2(b+ + 11 − b11 ) − 261(b20 − b20 )), + (1) 125π − − − L5 = − (110a+ + 02 − 90a02 + 60a20 − 48a20 − 250(b02 − b02 ) + (14) 6 14 LUIZ F. S. GOUVEIA AND J. TORREGROSA − 4b− + − 11 − 261(b20 − b20 )), (1) 32 − − − L6 = − (78120(a+ + + 02 − a02 ) + 42540(a20 − a20 ) − 203043(b20 − b20 ) 21 − + 194470(b+ 02 − b02 )), (1) 625π − − L7 =− (3675a+ 02 − 152565a02 − 85080a20 − 173200b02 + 2127 + 215740b− + 02 − 139233b20 + 266853b20 ), − (1) 2560 − − L8 =− (1942830(a+ + 02 − a02 ) − 21484712(b02 − b02 ) 402003 − + 1142229(b+ 20 − b20 )). We do not write here the next six because of their size but they have similar expressions (1) (1) as the above ones. We point out that L2 , . . . , L8 are linearly independent with respect (1) (1) (1) (1) to λ and L9 , . . . , L13 can be written as linear combinations of L2 , . . . , L8 . From Theorem 2.7 we can only ensure the existence of 8 limit cycles of small ampli- tude using a perturbation of type (13). So, we need to use second-order analysis in the parameter space for proving the statement using Theorem 2.8. With this aim, we change some of the old parameters a± ± k` , bk` by some new ones, namely {u2 , u3 , . . . , u8 }, such that, (1) after a linear change of variables in the parameter space we can write Lk = uk , for k = 2, . . . , 8. To unify notation, we denote by {u9 , . . . , u13 } the remaining old parameters a± ± k` , bk` that were not changed. Naming the parameters by u = (u2 , u3 , . . . , u13 ) ∈ R , 12 we have that the first Lyapunov constants can be written as Lk = uk + O2 (u), for k = 2, . . . , 8. Where O2 (u) denotes all monomials in the components of u with degree at least 2. In the next step we take a new change of coordinates again in the param- eter space, which is analytic and introduces new parameters {v2 , v3 , . . . , v8 } instead of {u2 , u3 , . . . , u8 }, such that, using the Implicit Function Theorem, we get Lk = vk , for k = 2, . . . , 8. Naming vk = uk for k = 9, . . . , 13, the perturbation parameters are now v = (v2 , . . . , v13 ) ∈ R12 . In fact, after this analytic change in the parameter space, the first seven Lyapunov constants are the first seven components in v. (2) The next step uses the second-order Taylor developments, Lk , of Lk with respect to v, for k = 9, . . . , 13. Using vk = 0 for k = 2, . . . , 8, and simplifying as in the proof (2) of Proposition 2.3, we obtain that Lk for k = 9, . . . , 13, are homogeneous polynomials of degree 2 on the relevant parameters v̂ = (v9 , . . . , v13 ) ∈ R5 . Note that under the vanishing conditions on the first components of v, we have reduced the problem to a 5-dimensional parameter space. The statement follows using Theorem 2.8 proving (2) that the varieties {Lk = 0, k = 9, . . . , 13}, intersect transversally along a straight line Υ ∈ R5 passing through the origin. Parameterizing the straight line Υ by v̂ = (w9 z, w10 z, w11 z, w12 z, z), we have that (2) Lk = z 2 Lk (w9 , . . . , w12 ), for k = 9, . . . , 13. Then, straightforward computations show that there exist two real solutions of {L9 = L10 = L11 = L12 = 0} such that L13 6= 0. These solutions write as w9∗ = α, w10 ∗ ∗ = p4 (α)/q3 (α), w11 = p2 (α)/q1 (α), ∗ w12 = p̂1 (α)/q̂1 (α), where pj , qj , p̂j , q̂j are polynomials of degree j with coefficients poly- nomials in π with rational coefficients, and α is a real root of a given polynomial φ of degree 2 with coefficients polynomials of degree 2 in π also with rational coefficients. We have not write here the polynomials pj , qj , p̂j , q̂j , or the polynomial φ because of their size. Moreover, there are only two solutions because φ has only two real solu- tions, which are simple. Additionally, we need to compute L∗13 = L13 (w9∗ , w10 ∗ ∗ , w11 ∗ , w12 ) LIMIT CYCLES IN POLYNOMIAL PIECEWISE SYSTEMS 15 and the determinant, D∗ , of the Jacobian matrix of (L9 , L10 , L11 , L12 ) with respect to ∗ (w9 , w10 , w11 , w12 ) evaluated at (w9∗ , w10 ∗ , w11 ∗ , w12 ) which are both rational functions of α. We notice that the involved rational functions are well defined because all the de- nominators are polynomials not vanishing at α. Even the numerators of L∗13 and D∗ . These conditions are satisfied by computing all the resultants with φ, with respect to α, and checking that are nonzero. Hence, the transversality of each straight line Υ is proved and so all the conditions for using Theorem 2.8 are satisfied. We finally remark that the two straight lines in the restricted parameter space are trans- formed into two analytic curves of weak-foci of order 13 that, with a versal unfolding, 13 limit cycles of small amplitude bifurcate from the origin. So the proof is complete.  Remark 3.2. In the last proof we have seen that from the origin of the quadratic center in Proposition 3.1, under a perturbation of type (9) having 12 parameters, bifurcate 11 limit cycles of small amplitude from the origin. As with a perturbation of type (13) we have more parameters but only two are relevant (see Proposition 2.1) we obtain in total 13 limit cycles of small amplitude. So, we think that we have obtained the maximum number of limit cycles that bifurcate from the origin because we have used (after a rescaling) all the essential parameters. Remark 3.3. From the expressions of the first-order Taylor developments of the Lya- − punov constants (14) it is clear that when the perturbation is not piecewise (a+ k` = ak` + − and bk` = bk` ) the Lyapunov constants having an even subscript vanish. Consequently, (1) (1) (1) we have only 3 linearly independent (L3 , L5 , L7 ) and three limit cycles of small am- plitude bifurcate from the origin as it is well-known. This is one simple example showing (1) that the upper bound M (2) = 3 is reached. In this case, L8 vanishes. But this is not the case when the perturbation is of piecewise type as we have already seen in the above proof. The first-order analysis of the system studied in the next proposition was done in [22]. Here we improve the result considering higher-order analysis, increasing the local cyclicity of the center perturbing also with piecewise polynomials of degree three. Proposition 3.4. Consider the perturbed system of the form (13) with n = 3 and Pc (x, y) = −yx(−68 + 1183x) and Qc (x, y) = −58x2 − 44xy + 30y 2 + 672x3 + 1484x2 y − 945xy 2 −84y 3 . Then, there exist small enough values of the parameters a± ± k` , bk` such that at least 26 hyperbolic crossing limit cycles of small amplitude bifurcate from the origin. Proof. The rational first integral (42x − 7y − 1)3 h(x, y) H(x, y) = , (448x2 + 336xy + 63y 2 − 44x − 12y + 1)3 (1183x2 − 68x + 1) with h(x, y) = 10752x3 + 29568x2 y + 17640xy 2 + 3024y 3 − 1600x2 − 2760xy − 576y 2 + 74x + 57y − 1, well defined at the origin, associated to the unperturbed system (13) provides the existence of a Darboux center at the origin. The proof follows the same scheme as the proof of Proposition 3.1 but here we will need to compute more Lyapunov constants because more limit cycles bifurcate from the origin. First we start restricting the analysis to a perturbation of type (9) and we will need to work with L2 , . . . , L26 . But due the difficulty of the computations we will not compute the complete expressions. Only the Taylor developments up to second- order with respect to the perturbative parameters, as we have done in the proof of Proposition 3.1. Here, again as we are perturbing centers, the Lyapunov constants vanish when all the perturbation parameters vanish. 16 LUIZ F. S. GOUVEIA AND J. TORREGROSA In [22] we explain that, straightforward computations show that the first-order Tay- (1) (1) lor series of the first Lyapunov constants, L2 , . . . , L24 , are linearly independent with (1) (1) respect to the perturbation parameters. In fact it can be proved that L25 and L26 are (1) (1) linear combinations of L2 , . . . , L24 . Hence, Theorem 2.7 ensures only the existence of (1) 24 limit cycles. We do not show here the complete expressions of Lk because of their size. To complete the proof of the statement we need a higher-order analysis. Using the same strategy than in the proof of Proposition 3.1 and using the linearity property of the paragraph above we can use a linear change of coordinates, from pa- rameters a± ± k` , bk` to u = (u2 , . . . , u29 ) ∈ R 28 such that Lk = uk + O2 (u) for k = 2, . . . , 24. (1) (1) Because we have 28 perturbative parameters in (9) but only L2 , . . . , L24 are linearly independent. Then, there is a small neighborhood of the origin in the parameter space such that, using the Implicit Function Theorem for introducing the new parameters v, Lk = vk , for k = 2, . . . , 24, after an analytic change of coordinates. In fact, as we will use the second-order developments of L25 and L26 this last change of variables is necessary to be obtained also up to second-order. Up to this change of coordinates and using that vk = 0, for k = 2, . . . , 24, we can restrict our analysis to v̂ = (v25 , v26 ) taking two extra conditions v27 = v28 = 0 to (2) (2) obtain L25 and L26 as homogeneous polynomials of degree 2 depending only on the two relevant parameters named v25 and v26 . Then we parameterize the straight line Υ (2) by v̂ = (w25 z, z) and we have Lk = z 2 Lk (w25 ), for k = 25, 26 with Lk polynomials of degree 2 in w25 with coefficients polynomials of degree 2 in π with rational coefficients. Moreover, L25 has simple real zeros where L26 is nonvanishing. Hence, the conditions of Theorem 2.8 are satisfied and the statement follows.  We notice that the cubic center that we have used to perturb in Proposition 3.4 is of Darboux type and it provides 11 limit cycles in the non-piecewise scenario. In fact, the limit cycles also appear only computing the linear parts of the Lyapunov constants and also considering the trace as another perturbation parameter. See [6]. As in Proposition 3.1, the local cyclicity in the piecewise scenario is higher but second-order analysis is also required. In the following results we have only used first-order Taylor developments of the corresponding Lyapunov constants due to the difficulties in the computations. Proposition 3.5. Consider the perturbed system of the form (13) with n = 4 and Pc (x, y) = x4 −6x2 y 2 +y 4 +x3 −3xy 2 +x2 −y 2 and Qc (x, y) = 4x3 y−4xy 3 +3x2 y−y 3 +2xy. Then, there exist small enough values of the parameters a± ± k` , bk` such that (13) has at least 36 hyperbolic crossing limit cycles of small amplitude bifurcating from the origin. Proof. The unperturbed system (13) with Pc and Qc as in the statement writes as the holomorphic differential system ż = i z + z 2 + z 3 + z 4 in complex coordinates z = x + i y. By the results in [19] it has an integrating factor and, consequently, a center at the origin. The proof follows using Theorem 2.7 from the computation of the first-order Taylor developments of the corresponding Lyapunov constants Lk , for k = 2, . . . , 36. Because (1) (1) L2 , . . . , L36 are linearly independent. We show only the expressions of the first three Lyapunov constants, because of their size, (1) 2 L2 = − (a+ − a− + − 11 + b20 − b20 ), 3 11 LIMIT CYCLES IN POLYNOMIAL PIECEWISE SYSTEMS 17 (1) 1 − − − − L3 = − π(a+ + + + 02 + a02 + a12 + a12 + 3a30 + 3a30 + 3b03 + 3b03 8 − − − 4b+ + 20 − 4b20 + b21 + b21 ), (1) − − − − L4 = −2(4b+ + + + + 20 − 4b20 + 2a02 − 2a02 − 3b21 + 3b21 + 2b22 − 2b22 − a11 + a− + − + − + − + 11 + 2a12 − 2a12 + 2a13 − 2a13 − 8a20 + 8a20 + 3a30 − 3a− + − + − + − + − 30 + 3a31 − 3a31 − 2b03 + 2b03 + 8b04 − 8b04 − 7b11 + 7b11 ).  Proposition 3.6. Consider the perturbed system of the form (13) with n = 4 and Pc (x, y) = (1 − x − y)y(1183x2 − 68x + 1) − y, and Qc (x, y) = (1 − x − y)(−672x3 − 1484x2 y + 945xy 2 + 84y 3 + 58x2 + 44xy − 30y 2 − x) + x Then, there exist small enough values of the parameters a± ± k` , bk` such that (13) has at least 40 hyperbolic crossing limit cycles of small amplitude bifurcating from the origin. Proof. The unperturbed system (13) where the polynomials Pc and Qc are the ones defined in the statement is the same cubic system given in Proposition 3.4 but multiplied by a straight line of equilibrium points that does not passes through the origin. Then, we have also a center at the origin. Straightforward computations show that the first-order Taylor developments with respect to the parameters of the first Lyapunov constants, Lk , for k = 2, . . . , 40, are linearly independent. Hence, the statement follows as the previous results using Theo- rem 2.7. Here, we show only the expressions of the first three Lyapunov constants up to a first-order Taylor development, because of their size, (1) 2 L2 = − (a+ − a− + − + − 11 + 2b02 − 2b02 + b20 − b20 ), 3 11 (1) π L3 = (129a− − − − − + 0,2 + 46a1,1 + a1,2 + 187a2,0 + 3a3,0 + 129a0,2 8 − − − 44a+ + + + 1,1 + a1,2 + 187a2,0 + 3a3,0 + 137b0,2 + 3b0,3 + 29b− − − + + + + 1,1 + 90b2,0 + b2,1 − 43b0,2 + 3b0,3 + 29b1,1 + b2,1 ), (1) 2 L4 = (−9b− + + + + + + 4,0 − 4164a0,2 + 774a0,3 − 82a1,2 + 6a1,3 + 11318a2,0 + 381a2,1 45 + 141a+ + + + + + + 3,0 + 9a3,1 − 17988b0,2 + 24b0,3 + 24b0,4 + 4096b1,1 − 6b1,2 − − − − 1013b+ + + + 2,0 + 182b2,1 + 111b3,0 + 9b4,0 + 4164a0,2 − 774a0,3 + 82a1,2 − 6a− − − − − − − 1,3 − 11318a2,0 − 381a2,1 − 141a3,0 − 9a3,1 + 17988b0,2 − 24b0,3 − 24b− − − − − − 0,4 − 4096b1,1 + 6b1,2 + 1013b2,0 − 182b2,1 − 111b3,0 )  Proposition 3.7. Consider the perturbed system of the form (13) with n = 5 and Pc (x, y) = x5 − 10x3 y 2 + 5xy 4 + x4 − 6x2 y 2 + y 4 + x3 − 3xy 2 + x2 − y 2 , and Qc (x, y) = 5x4 y − 10x2 y 3 + y 5 + 4x3 y − 4xy 3 + 3x2 y − y 3 + 2xy. Then, there exist small enough values of the parameters a± ± k` , bk` such that (13) has at least 58 hyperbolic crossing limit cycles of small amplitude bifurcating from the origin. Proof. The unperturbed system (13) where the polynomials Pc and Qc are the ones defined in the statement, in complex coordinates z = x + i y, writes as ż = i z + z 2 + z 3 + z 4 + z 5 . As in Proposition 3.5, it is an holomorphic planar differential system with an equilibrium point of monodromic type. Consequently, it has a center at the origin. See also [19]. 18 LUIZ F. S. GOUVEIA AND J. TORREGROSA The statement follows as the previous results, using again Theorem 2.7, because first- order Taylor developments of the first Lyapunov constants, L2 , . . . , L58 , are linearly independent. The computations are straightforward. Because of their size, we only show the first-order developments of the first three Lyapunov constants. (1) 2 − − − L2 = (−(a+ + + 11 − a11 ) − 2(b02 − 2b02 ) − (b20 − b20 )), 3 (1) 1 − − − − L3 = − π(a+ + + + 12 + a12 + 3(a30 + a30 ) − (b02 + b02 ) + 3(b03 + b03 ) 8 − − − 4(b+ + 20 + b20 ) + b21 + b21 ), (1) 2 L4 = − (3b+ − 3b− − + + − − + + − 40 + 3b21 − 3b21 + 4b20 − 4b20 + 2b03 − 2b03 − 7b11 + 7b11 15 40 − 3a− + − + + − + − + 30 + 3a30 − 8b04 + 8b04 − 8a20 + 8a20 + 3a31 − 3a31 + 2a13 − 2a13 − − 6b− + + − + − + − + − 02 + 6b02 + 2a12 − 2a12 + 2b22 − 2b22 − a11 + a11 − 12a02 + 12a02 ).  4. How to deal with the computational difficulties The main computational difficulty found in this paper is directly related to the mech- anism itself to obtain the Lyapunov constants. That is, the computation of the series expansion of the solution in polar coordinates, with respect to the initial condition in the x-axis. This method needs a high computational effort because of the huge trigonometrical expressions appearing in the functions ri± (θ) defined in (6) solving the differential equation (5). We have shown some of them in the proof of Proposition 2.2. This obstacle can be skipped in analytic vector fields using other mechanisms like the Poincaré or the normal form methods, see more details of them in, for example, the book of Romanovsky and Shafer [34]. But, up to our knowledge, for piecewise differential type systems, the Lyapunov mechanism is the only valid. Another hurdle is the huge size that the Lyapunov constants have, even for ana- lytic vector fields. In fact, the Lyapunov constants for cubic polynomial systems have been obtained very recently in [37] and they are unknown even for piecewise quadratic polynomial systems. In our context, as it can be seen in Section 2, both the size and the number of constants needed to be calculated are doubled. As our objective is to study the local cyclicity, in fact lower bounds of it, we study center perturbations. So, as Christopher showed in [11], we do not need to compute the complete Lyapunov constants. Only the first or second-order Taylor developments with respect to the per- turbative parameters is necessary to be obtained. We have used this advantage in [21] for proving the best lower bounds for the local cyclicity problem in polynomial vector fields of low degree. The computation of the Taylor developments can be also simplified using the par- allelization procedure described in [23]. The main goal of this technique, originally developed only for first-order developments in [30], is to decompose the total compu- tation in a collection of simpler problems. We will explain here the case when we have a perturbation with m parameters and such that only first-order developments are needed. Denoting by λ = (λ1 , . . . , λm ) the perturbation parameters, the first-order Taylor development of each k-Lyapunov constant can be written as a linear combina- (1) Pm tion of the components of λ, that is Lk (λ) = αi λi . We notice that we have used i=1 that, as we are perturbing centers, Lk (0) = 0. Then each αi can be obtained considering LIMIT CYCLES IN POLYNOMIAL PIECEWISE SYSTEMS 19 the perturbation problem that has only the perturbation monomial corresponding to λi , vanishing all the others, and computing the first-order Taylor development with respect to this λi . The strategy decrease the total computation time because, having a cluster of computers, we can get all αi simultaneously. Even if we use this important improve- ment, we have only worked with piecewise polynomial vector fields up to degree n = 4 and only the Taylor series up of first and second-orders, because both the computation time and the memory requirements are very high. For example, in the proof of the last result, the total computation time has been around 4 days, working with 5 computers using the parallelization method. In particular, the memory requirements for the com- putation of the first-order Taylor series forces us to use only 4 perturbative monomials in each computer. The advantage of parallelization here was clear because without it the calculations would have been impossible to obtain. Therefore, to go further in a higher degree analysis, it is necessary to develop a new mechanism for calculating the coefficients of the return map. The last trick has been the computation of r(θ, ρ, λ) and then evaluating at ±π but changing λ by λ± . Considering r± (θ, ρ, λ± ) in (6) as two different expressions the computation time would be doubled. 5. Acknowledgments This work has been realized thanks to the AGAUR 2017 SGR 1617 grant, the Ministerio de Ciencia, Innovación y Universidades - Agencia Estatal de Investigación MTM2016-77278-P (FEDER) and PID2019-104658GB-I00 grants, the European Com- munity H2020-MSCA-RISE-2017-777911 grant, and the Brazilian CNPq 200484/2015-0 and FAPESP 2020/04717-0 grants. We would like to thank the referees for their careful reading. With their comments we have improved the presentation and the proofs of the results. References [1] V. Acary, O. Bonnefon, and B. Brogliato. Nonsmooth modeling and simulation for switched cir- cuits, volume 69 of Lecture Notes in Electrical Engineering. Springer, Dordrecht, 2011. [2] 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. [3] A. A. Andronov, A. A. Vitt, and S. E. Khaikin. Theory of oscillators. Translated from the Russian by F. Immirzi; translation edited and abridged by W. Fishwick. Pergamon Press, Oxford-New York-Toronto, Ont., 1966. [4] 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. [5] M. di Bernardo, C. J. Budd, A. R. Champneys, and P. Kowalczyk. Piecewise-smooth dynamical systems, volume 163 of Applied Mathematical Sciences. Springer-Verlag London, Ltd., London, 2008. Theory and applications. [6] Y. L. Bondar and A. P. Sadovskiı̆. On a theorem of Zoladek. Differ. Uravn., 44(2):263–265, 287, 2008. [7] C. Buzzi, C. Pessoa, and J. Torregrosa. Piecewise linear perturbations of a linear center. Discrete Contin. Dyn. Syst., 33(9):3915–3936, 2013. [8] J. Castillo, J. Llibre, and F. Verduzco. The pseudo-Hopf bifurcation for planar discontinuous piecewise linear differential systems. Nonlinear Dynam., 90(3):1829–1840, 2017. [9] C. Chicone and M. Jacobs. Bifurcation of critical periods for plane vector fields. Trans. Amer. Math. Soc., 312(2):433–486, 1989. [10] C. Chicone and M. Jacobs. Bifurcation of limit cycles from quadratic isochrones. J. Differential Equations, 91(2):268–326, 1991. 20 LUIZ F. S. GOUVEIA AND J. TORREGROSA [11] C. Christopher. Estimating limit cycle bifurcations from centers. In Differential equations with symbolic computation, Trends Math., pages 23–35. Birkhäuser, Basel, 2005. [12] 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. [13] B. Coll, A. Gasull, and R. Prohens. Degenerate Hopf bifurcations in discontinuous planar systems. J. Math. Anal. Appl., 253(2):671–690, 2001. [14] S. Coombes. Neuronal networks with gap junctions: a study of piecewise linear planar neuron models. SIAM J. Appl. Dyn. Syst., 7(3):1101–1129, 2008. [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] A. F. Filippov. Differential equations with discontinuous righthand sides, volume 18 of Mathemat- ics and its Applications (Soviet Series). Kluwer Academic Publishers Group, Dordrecht, 1988. Translated from the Russian. [17] E. Freire, E. Ponce, and F. Torres. Canonical discontinuous planar piecewise linear systems. SIAM J. Appl. Dyn. Syst., 11(1):181–211, 2012. [18] E. Freire, E. Ponce, and F. Torres. The discontinuous matching of two planar linear foci can have three nested crossing limit cycles. Publ. Mat., 58(suppl.):221–253, 2014. [19] A. Garijo, A. Gasull, and X. Jarque. Local and global phase portrait of equation ż = f (z). Discrete Contin. Dyn. Syst., 17(2):309–329, 2007. [20] F. Giannakopoulos and K. Pliete. Planar systems of piecewise linear differential equations with a line of discontinuity. Nonlinearity, 14(6):1611–1632, 2001. [21] J. Giné, L. F. S. Gouveia, and J. Torregrosa. Lower bounds for the local cyclicity for families of centers. To appear in J. Differential Equations. [22] L. F. S. Gouveia and J. Torregrosa. 24 crossing limit cycles in only one nest for piecewise cubic systems. Applied Mathematics Letters, 103:106189, 2020. [23] L. F. S. Gouveia and J. Torregrosa. Lower bounds for the local cyclicity of centers using high order developments and parallelization. J. Differential Equations, 271:447–479, 2021. [24] L. Guo, P. Yu, and Y. Chen. Bifurcation analysis on a class of Z2 -equivariant cubic switching systems showing eighteen limit cycles. J. Differential Equations, 266(2-3):1221–1244, 2019. [25] M. Han and W. Zhang. On Hopf bifurcation in non-smooth planar systems. J. of Differential Equations, 248:2399–2416, 2010. [26] S. M. Huan and X. S. Yang. On the number of limit cycles in general planar piecewise linear systems. Discrete Contin. Dyn. Syst., 32(6):2147–2164, 2012. [27] Y. Ilyashenko. Centennial history of Hilbert’s 16th problem. Bull. Amer. Math. Soc. (N.S.), 39(3):301–354, 2002. [28] R. Leine and D. van Campen. Discontinuous bifurcations of periodic solutions. Mathematical and Computer Modelling, 36(3):259 – 273, 2002. [29] C. Li, C. Liu, and J. Yang. A cubic system with thirteen limit cycles. J. Differential Equations, 246(9):3609–3619, 2009. [30] 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. [31] J. Llibre and E. Ponce. Three nested limit cycles in discontinuous piecewise linear differential systems with two zones. Dyn. Contin. Discrete Impuls. Syst. Ser. B Appl. Algorithms, 19(3):325– 335, 2012. [32] J. Llibre and Y. Tang. Limit cycles of discontinuous piecewise quadratic and cubic polynomial perturbations of a linear center. Discrete Contin. Dyn. Syst. Ser. B, 24(4):1769–1784, 2019. [33] R. Prohens and J. Torregrosa. New lower bounds for the Hilbert numbers using reversible centers. Nonlinearity, 32(1):331–355, 2019. [34] V. G. Romanovski and D. S. Shafer. The center and cyclicity problems: a computational algebra approach. Birkhäuser Boston, Ltd., Boston, MA, 2009. [35] R. Roussarie. Bifurcation of planar vector fields and Hilbert’s sixteenth problem, volume 164 of Progress in Mathematics. Birkhäuser Verlag, Basel, 1998. [36] S. L. Shi. A concrete example of the existence of four limit cycles for plane quadratic systems. Sci. Sinica, 23(2):153–158, 1980. [37] I. Sánchez-Sánchez and J. Torregrosa. New advances on the lyapunov constants of some families of planar differential systems. Research Perspectives CRM, pages 161–167. Springer, Barcelona, 2019. LIMIT CYCLES IN POLYNOMIAL PIECEWISE SYSTEMS 21 [38] H. Żoladek. Quadratic systems with center and their perturbations. J. Differential Equations, 109(2):223–273, 1994. Departament de Matemàtiques, Universitat Autònoma de Barcelona, 08193 Cer- danyola del Vallès, Barcelona (Spain) 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 Cer- danyola del Vallès, Barcelona (Spain) Centre de Recerca Matemàtica, Campus de Bellaterra, 08193 Cerdanyola del Vallès, Barcelona (Spain) Email address: torre@mat.uab.cat