This is a preprint of: “Limit cycles in 4-star-symmetric planar piecewise linear systems”, Claudio A. Buzzi, Joao Carlos Medrado, Joan Torregrosa, J. Differential Equations, vol. 268(5), 2414– 2434, 2020. DOI: [https://doi.org/10.1016/j.jde.2019.09.008] LIMIT CYCLES IN 4-STAR-SYMMETRIC PLANAR PIECEWISE LINEAR SYSTEMS CLAUDIO A. BUZZI, JOAO C. MEDRADO, AND JOAN TORREGROSA Abstract. Our interest is centered in the study of the number of limit cycles for nonsmooth piecewise linear vector fields on the plane when the switching curve is xy = 0. We consider the symmetric case. That is, one vector field defined in the odd quadrants and the other in the even ones. We deal with equilibrium points of center-focus type, with matrices in real Jordan form, in each vector field when the infinity is monodromic. In this case, we provide the center classification at infinity, we prove that the maximum order of a weak focus is five. Moreover, we show the existence of systems exhibiting five limit cycles bifurcating from infinity. 1. Introduction These last years a big interest appeared for studying nonsmooth piecewise differ- ential systems, motivated mainly by their biological and engineering applications. In particular, piecewise linear differential systems have been used to model many real processes and different modern devices, for more details see for instance [1, 11] and the references therein. In piecewise-smooth dynamical systems, the phase space is divided in regions with different smooth vector fields in each one. At the boundary between smooth regions, the vector field could be neither continuous nor differentiable. But it can be defined using the Filippov convention, see [14]. Such systems arise in several applications, and their analysis constitutes an exciting field of current research. In fact, this nonsmoothness may cause a whole new rich family of bifurcations that do not appear in smooth systems. Due to the nondifferentiability, the analysis of such bifurcations prevents the use of standard local methods (see for example [2]), and thus demands the development of new techniques. In general, many works consider the piecewise smooth vector fields when the switching manifold (or discontinuity manifold), Σ = f −1 (p), is a submanifold given by the preimage of a regular value, p, of a differentiable real function. Although there are few works considering the case when the switching manifold has a singularity (it is the preimage of a nonregular value of a function), we can found some real phenomena modeled by this type of systems. Different planar models with f (x, y) = xy can be found in the literature. For example, a model of a steam-engine ([3]), a periodically forced double-belt friction oscillator ([13]), or a two-mass model of the human vocal folds used to simulate their oscillatory motion at phonation ([22]). There are models in higher dimensions, (x1 , x2 , . . . , xn ), where the nonlinearities depend also on |xi |, for i = 1, . . . , n. These piecewise differential equations are defined with discontinuity varieties having also a nonregular intersection, that is 2010 Mathematics Subject Classification. Primary: 34C25, 34C07, 37C27. Key words and phrases. Non-smooth differential system, 4-cross-symmetric, cyclicity, limit cy- cles, center-focus problem. 1 2 C. BUZZI, J. C. MEDRADO, AND J. TORREGROSA x1 x2 · · · xn = 0. For example, among other bifurcations, the study of periodic orbits by a Hopf bifurcation in a model for ship maneuvering is studied in [4]. An interesting phenomenon that appears in the analysis of the models described by (smooth or not) differential equations, is the existence of periodic behaviors. In particular the study of the number, stability, location,. . . of isolated periodic orbits. When we restrict to planar polynomial vector fields the problem of finding the maximum number of limit cycles remains open. It is part of the famous 16th- Hilbert problem. This problem make no sense for linear vector fields because they have no limit cycles. But as there exist limit cycles crossing Σ for piecewise linear vector fields we can consider the generalized 16th-Hilbert problem. We call Hilbert number to this maximum when it exists. The Hilbert number for piecewise linear vector fields defined in two zones separated by a straight line is only known for some particular cases. For the continuous case, this number was conjectured to be one by Lum and Chua in [23] and proved by Freire et al. in [15]. For the discontinuous case without sliding motion, see [14], in [24] was proved that the Hilbert number is also one. When there is sliding motion, the best known results show that the Hilbert number is bigger or equal than three (see [6, 16, 17, 21]). In this paper we consider the class of 4-star-symmetric planar piecewise linear systems. That is, when Σ = {(x, y) : xy = 0} and there are only two piecewise linear systems defined in the odd and even quadrants, respectively. Our main result shows that the Hilbert number in this class is at least five. The same lower bound is also obtained in [7] studying the case having one linear system defined in the first quadrant and another in the other quadrants. We deal with piecewise linear systems of the form  0 x = α± x − β ± y + γ ± , (1) y 0 = β ± x + α± y + δ ± , defined in Σ± = {(x, y) ∈ R2 : ±xy > 0} when β ± 6= 0. On Σ = {(x, y) ∈ R2 : xy = 0}\{(0, 0)} we use the Filippov convention, see [14]. As we are interested in periodic motions, we will study only the case β + β − > 0. By means of a time-rescaling, if necessary, and defining a± = α± /β ± , b± = γ ± /β ± , and c± = δ ± /β ± , system (1) becomes  0 x = a± x − y + b ± , (2) y 0 = x + a± y + c± . It is not restrictive, doing the change (x, y, t) 7→ (−x, y, −t) if necessary, to assume also that a+ > 0. See Figure 1. When b± = c± = 0 the origin is the unique equilibrium point, which is of center- focus type, for each linear system in (2) defined in Σ± , respectively. Consequently, the origin is an equilibrium point of monodromic type for system (2). Its stability is given by the sign of a+ + a− , see [28]. That is, stable (unstable) when it is negative (positive). Moreover, the Poincaré return map is globally linear and the origin is a global focus (a+ + a− 6= 0) or a global center (a+ = −a− ). In this case, we have no limit cycles. Otherwise, some sliding motion can appear near the origin of the vector field (2), in particular the monodromy property could fail. But this property is kept near the infinity when we extend (2) to the Poincaré disk, see [12]. In Figure 1, the solutions are depicted far from the origin. We present results about the existence of limit cycles and its maximum number for system (2) near the infinity. In this neighborhood all the trajectories are of crossing type. This is done by an adequate LIMIT CYCLES IN 4-STAR-SYMMETRIC PLS 3 Figure 1. The phase portrait of system (2) near infinity. polar coordinates change, (x, y) = (R−1 cos θ, R−1 sin θ), that moves the infinity to the origin. After this change, this problem is related with the classical center-focus problem. More concretely, we study the series expansion of the return map near the origin and we compute the Lyapunov quantities. More details about this theory for smooth systems can be seen in [2]. Our computations deal near a focus-focus case. The first stability condition for a general system with an equilibrium point at the origin of focus type and defined in the regions depicted in Figure 1 was done in [19, 28]. Previously, but when the discontinuity curve is a straight line, the first terms of the series expansion were obtained also for the fold-fold and fold-focus cases in [10]. The approach that we will use is closer to the classical smooth one. For a general introduction to the center-focus problem for smooth systems see [26]. In this paper we extend this theory to be applied in our class of piecewise vector fields. Next, we present the main results of this work. In the first we provide a lower bound for the maximum number of isolated periodic orbits (limit cycles) for sys- tem (2) near the infinity. In the second we classify the centers at infinity. We observe that the listed families are not disjoint. Theorem 1.1. There exist values of the parameters such that system (2) has five limit cycles bifurcating from infinity. Theorem 1.2. System (2) has a center at infinity if and only if the parameters belong to one of the next families: F1 ={b− + b+ = c− − c+ = 0}, F2 ={(Aa − 1)b− − (A + a)c− = (Aa − 1)b+ + (A + a)c+ = 0}, F3 ={b− − b+ = c− + c+ = 0}, F4± ={(A ± a)b− − (A ∓ a)b+ = (A ± a)c− − (Aa ± 1)b+ = (A ± a)c+ + (Aa ∓ 1)b+ = 0}, C1 ={a = b− − b+ = c− − c+ = 0}, C2 ={a = b− − c− = b+ + c+ = 0}, C3± ={a = b− ± c− = b+ ± c+ = 0}, π where a = a+ = −a− and A = e−a 2 . 4 C. BUZZI, J. C. MEDRADO, AND J. TORREGROSA This paper is structured as follows. In Section 2 we introduce the technical results to define, via the coefficients of the series expansion of the solution near the origin, the center conditions (Lyapunov quantities). They are used, in Section 3, to classify the centers for proving Theorem 1.2. The study of the isolated periodic orbits is done in Section 4 with a generalization of the classical degenerate Hopf bifurcation. We study when the infinity is a weak focus and which is its maximal order. Besides the proof of Theorem 1.1 we prove that there are no weak foci of order greater than five. The five limit cycles appear proving the existence of a weak focus of order five, via Poincaré–Miranda’s Theorem ([20]), together with a generic unfolding, using the Implicit Function Theorem, in family (2). To prove the existence of a generic unfolding we closely follow the approach used in [8, 9, 18]. As this problem has only six parameters (a± , b± , c± ), it is quite natural that the number of limit cycles be five. We notice that the difficulties for proving the existence of such limit cycles are devoted to the fact that our perturbation should be restricted to be inside the considered family. When the perturbation is analytic the explicit computations are not necessary to be done. See, for example, [5]. The computations are done with the Computer Algebra System MAPLE and with a rigorous rational interval analysis. 2. Canonical form and the Poincaré return map As we are interested in the periodic orbits bifurcating from infinity, by the mon- odromy property, the natural coordinates to study the problem are (x, y) = (R−1 cos θ, R−1 sin θ). (3) Then, a neighborhood of infinity of (2) is moved to a neighborhood of the origin of equation dR −aR − (b cos θ + c sin θ)R2 R0 = = f (R, θ) = , (4) dθ 1 + (c cos θ − b sin θ)R where (a, b, c) = (a+ , b+ , c+ ) if θ ∈ (0, π/2) ∪ (π, 3π/2) and (a, b, c) = (a− , b− , c− ) if θ ∈ (π/2, π) ∪ (3π/2, 2π). Clearly, for R small enough the solution is well defined for θ ∈ [0, 2π]. We consider, for each i = 1, 2, 3, 4, and ρ small enough, the solution Ri (θ, ρ) of (4) satisfying the initial condition Ri ((i − 1)π/2, ρ) = ρ. Then we define the first return map on each quadrant i = 1, 2, 3, 4 as Πi (r) = Ri (iπ/2, ρ) and the total return map as the ordered composition of all of them, Π(ρ) = Π4 (Π3 (Π2 (Π1 (ρ)))). We observe that Π is analytic, in a neighborhood of the origin, because it is defined by composition of analytic maps. Consequently, the periodic orbits of (2) are the fixed points of Π or, alternatively, the zeros of a displacement function ∆(ρ) = (Π2 ◦ Π1 )(ρ) − (Π−1 −1 3 ◦ Π4 )(ρ). Next result relates the first return maps in each quadrant with the defined in the first one. Lemma 2.1. Let R(θ, ρ) be the solution of (4) satisfying the initial condition R(0, ρ) = ρ, then we write Π(a,b,c) (ρ) = R(π/2, ρ). Then, Π1 (ρ) = Π(a+ ,b+ ,c+ ) (ρ), Π2 (ρ) = Π(a− ,c− ,−b− ) (ρ), Π3 (ρ) = Π(a+ ,−b+ ,−c+ ) (ρ), and Π4 (ρ) = Π(a− ,−c− ,b− ) (ρ). The proof follows doing an adequate translation in the angle variable that moves the respective quadrant to the first one. From this lemma we only need the expres- sion of Π(a,b,c) (ρ). The series expansion in ρ can be obtained from the series expan- sion of the solution R(θ, ρ) of (4) with R(0, ρ) = ρ. That is, substituting R(θ, ρ) = LIMIT CYCLES IN 4-STAR-SYMMETRIC PLS 5 P ∞ Ri (θ)ρi , with R1 (0) = 1 and Ri (0) = 0, for i = 2, . . . , in (4). Straightforward i=1 computations show that R1 (θ) = e−aθ and R2 (θ) = e−2aθ (b sin θ − c cos θ) + e−aθ c. We do not show the other terms because of the size of them. From the previous definition we can write  2  −a π2 a c − 2ab − c −a π a2 b + 2ac − b −aπ 2 Π(a,b,c) (ρ) =R(π/2, ρ) = e ρ+ e 2 + e ρ a2 + 1 a2 + 1  4 2 a c − 6a3 bc + 7a2 b2 − 5a2 c2 + 10abc − b2 + 2c2 −a π + e 2 2(a2 + 1) 2(a2 c − 2ab − c)(a2 b + 2ac − b) −aπ + e a2 + 1  3a4 b2 + 10a3 bc − 3a2 b2 + 9a2 c2 − 6abc + 2b2 + c2 −3a π 3 + e 2 ρ + O(ρ4 ). 2(a2 + 1) Using Lemma 2.1 we can get the first terms of the series development of the complete Poincaré map Π(ρ), or equivalently, the displacement function ∆(ρ) = P∞ i i=1 R̂i ρ , with + +a− ) π + +a− ) π R̂1 = e−(a 2 − e(a 2 , (a− )2 c− − 2a− b− − c− −(a− +a+ )π − + π  R̂2 = − 2 e − e(a +a ) 2 (a ) + 1 (a ) c − 2a+ b+ − c+ −(a− +a+ ) π + 2 + (a− +a+ )π  + e 2 −e (a+ )2 + 1 1 − ((a− )2 (a+ )2 − 1)(b− − b+ ) ((a ) + 1)((a− )2 + 1) + 2 + ((a− )2 − (a+ )2 )(b− + b+ ) + 2a− a+ (a+ c− − a− c+ )  − + π − + π  + 2(a− c− − a+ c+ ) e−(a +2a ) 2 − e(2a +a ) 2 . The stability of the origin for system (4) is given by the sign of the first nonvanishing coefficient R̂i . The next proposition is equivalent to the classical result for the center- focus problem in the smooth scenario. Proposition 2.2. The first non-vanishing term R̂i corresponds to an odd index i. Proof. The proof is the same as the one provided in [2, Chap. IX, Sect. 24] for the smooth context. Because, it uses only the fact the function f, in (4), satisfies the symmetry relation f (−R, θ + π) = −f (R, θ).  We notice that, in the general case with four quadrants without symmetry, the above result is not valid. It is natural to define the first nonzero value when the previous ones are zero. Then, using the above result, we define the Lyapunov quantity of order i of the origin of system (4) as + +a− ) π + +a− ) π L0 = R̂1 = e−(a 2 − e(a 2 (5) or Li = R̂2i+1 |{L1 =···=Li−1 =0} , for i = 1, 2, . . . 6 C. BUZZI, J. C. MEDRADO, AND J. TORREGROSA So, as L0 = 0 if and only if a+ + a− = 0, straightforward computations show that (a+ c− − a+ c+ + b− + b+ )(a+ c− + a+ c+ + b− − b+ ) L1 = (a+ )2 + 1 (6) (a+ b− + a+ b+ − c− + c+ )(a+ b− − a+ b+ − c− − c+ ) −πa+ − e . (a+ )2 + 1 Although in the next section we will need the explicit expression of the first five Lyapunov quantities, we do not present all of them because of their size. This symmetric case reproduces the same properties as the smooth scenario. Then we say that the origin of (4) is a hyperbolic focus when L0 6= 0. When L0 = L1 = · · · = Li−1 = 0 and Li 6= 0 the origin of (4) is a weak focus of order i. In fact, the case L0 = 0 corresponds, for an analytic vector field, to a monodromic equilibrium point with a Jacobian matrix having trace equal to zero and positive determinant. We will say that the origin is an equilibrium point of (4) of center-focus type. The next result corresponds to the bifurcation of a limit cycle from the origin by a Hopf bifurcation when we perturb a system having a weak focus of order 1. Proposition 2.3. Let us consider equation (4) with a+ + a− = ε and when ε = 0 we assume that L1 < 0(> 0) in (6). Then, for ε = 0 the origin is a stable (unstable) weak focus of order 1 of (4). Moreover, for ε < 0(> 0) small enough, a stable (unstable) limit cycle bifurcates from the origin. The degenerate Hopf bifurcation, obtaining k limit cycles from a weak focus of order k, can be also studied similarly as in the classical theory of analytic systems. See [2, 27]. 3. The center problem The studies of the necessary and sufficient conditions in order that system (4) has a center at the origin or system (2) has a center at infinity are equivalent, due to the change of coordinates introduced in previous section. The proof of Theorem 1.2 fol- lows from the results of this section. Proposition 3.2 proves the necessary conditions and Proposition 3.3 proves the sufficiency of them. From the results of the above section we can assume L0 = 0, that is a+ = −a− = a, see (5). From the expression L1 in (6) we will introduce a change of variables in the parameter space, because the size of the expressions of the Lyapunov quantities decrease. This is done in the next lemma, which provides the structure of the first five Lyapunov quantities for system (4). π Lemma 3.1. Consider system (4) with a+ = a, a− = −a, A = e−a 2 and B1 =a(b− + b+ ) − (c− − c+ ), C1 =a(c− − c+ ) + b− + b+ , (7) B2 =a(b− − b+ ) − (c− + c+ ), C2 =a(c− + c+ ) + b− − b+ . Then, the origin is an equilibrium point of center-focus type of (4), L0 = 0, and the first five Lyapunov quantities write as Mk (a, A, B1 , B2 , C1 , C2 ) Lk = (8) (a2 + 1)2k−1 for some polynomials Mk with rational coefficients of degree 6k − 2. Moreover, they are also homogeneous polynomials of degree 2k in the variables (B1 , B2 , C1 , C2 ). LIMIT CYCLES IN 4-STAR-SYMMETRIC PLS 7 More concretely, M1 = − (A2 B1 B2 − C1 C2 ), M2 = aA3 (3A2 − 8aA + 9)B13 B2 − A3 (4aA − 3)B13 C2 + 3A2 (3a2 A − A − 4a)B12 B2 C1 + aA3 (3A2 + 1)B1 B23 (9) − 3A3 (4aA + 3a2 − 1)B1 B22 C2 + aA(A2 + 3)B1 B2 C12 + aA(9A2 + 8aA + 3)B1 B2 C22 − (3A + 4a)B1 C23  + a2 A3 B23 C1 − a2 AB2 C13 /12. We have not written all the polynomials Mk because of their size. In particular, M3 , M4 , and M5 have 303, 1175, and 3398 monomials, respectively. π Proposition 3.2. Let a+ = a, a− = −a, A = e−a 2 and B1 , B2 , C1 , and C2 defined in (7). The conditions F1 = {B1 = C1 = 0}, C1 = {a = B1 = C2 = 0}, F2 = {C1 = AB1 , C2 = AB2 }, C2 = {a = 0, B1 = −C1 , B2 = −C2 }, F3 = {B2 = C2 = 0}, C3± = {a = 0, B1 = ±C2 , B2 = ±C1 }, F4± = {B1 = C2 = 0, C1 = ±AB2 }, are necessary for system (4) has a center at the origin. We notice that the families defined in the above result are exactly the ones ap- pearing in Theorem 1.2. Proof of Proposition 3.2. The proof follows studying when all the Lyapunov quan- tities that appear in Lemma 3.1 vanish. First we consider the condition B2 = 0 and then B2 6= 0. We start with the condition B2 = 0. From the expression of M1 in (9) it is clear that if L1 = 0 then C2 = 0 or C1 = 0. The first case is in F3 . For the second case, that is C1 = 0 and C2 6= 0, computing the next Lyapunov quantities, see (8), we obtain Mk = B1 C2 Pk (a, A, B1 , C2 ), for k = 2, 3, 4, 5, with Pk a polynomial with rational coefficients of degree 6k − 5. Clearly, B1 = 0 vanish all the quantities and defines family F1 . Only remains the study of vanishing all the polynomials Pk . Computing the crossed resultants of Pk with P2 we obtain R2,k = a2 A8 (A2 + 1)2 C24k−4 Qk (a, A)2 , for k = 3, 4, 5, with Qk polynomials with rational coefficients of degrees 6, 13, and 20, respectively. When a = 0 all R2,k vanish identically and we obtain C3± . The proof of the case B2 = 0 finishes showing that, when a 6= 0, B1 6= 0, C2 6= 0, and C1 = 0, the system has a weak focus of order 4 or 3, depending if Q3 (a, A) is zero or not. The existence of a focus of order 4 implies that Q4 is not zero when Q3 is. This is because, if the curves Q3 (a, A) = 0 and Q4 (a, A) = 0 have an intersection point (a∗ , A∗ ), a∗ will be a root of the resultant of Q3 and Q4 with respect to A, which is a polynomial with rational coefficients of degree 26 in a2 . Similarly, A∗ will be a root of another polynomial, also with rational coefficients, of degree 27 in A2 . Consequently, both a∗ and A∗ will be algebraic numbers and, from the definition of ∗ A, we have that (A∗ )2 = e−πa . This last assertion contradicts Gelfond–Scheneider Theorem, see [25, Pag. 134], which says that “If α and β are algebraic number with α 6= 0, α 6= 1, and β 6∈ Q then αβ is transcendental.” Because, as eπ = i−2i , we will ∗ ∗ have (A∗ )2 = e−πa = i2ia is transcendental but it is algebraic. 8 C. BUZZI, J. C. MEDRADO, AND J. TORREGROSA We continue assuming now that B2 6= 0. From the expression of M1 in (9) we have B1 = C1 C2 /(A2 B2 ), then the next Lyapunov quantities, by (8), write Mk = C1 (C2 − AB2 )B21−2k A5−4k Pk (a, A, B2 , C1 , C2 ), (10) for k = 2, 3, 4, 5, with each Pk is a polynomial with rational coefficients of degree 12(k − 1). Consequently, there are two simple cases, C1 = 0 and C2 = AB2 , that define families F1 and F2 , respectively. Only remains the study of vanishing all the polynomials Pk when B2 6= 0, C1 6= 0, and C2 6= AB2 . The intersection points of the varieties Pk = 0 can be studied also from the two by two crossed resultants. Straightforward computations show that res(P2 , Pk , C1 ) = a2 A2k (A2 + 1)2 B24k−4 C22 (AB2 − C2 )4 R2,k (a, A, B2 , C2 ), (11) for k = 3, 4, 5, for some R2,k polynomials with rational coefficients of degree 17k−27. Here we have also two simple cases, a = 0 or C2 = 0. In the first case, Pk = C22 (B2 + C2 )(B22 − C12 )Qk (B2 , C1 , C2 ), for k = 2, 3, 4, 5, and Qk some polynomials with rational coefficients with degree 4(k − 2). Then, as Q2 is a non vanishing constant, the factors of P2 define all the center families, that is, C1 , C2 , and C3± . In the second case, Pk = a2 A4k−5 B22k−1 (A2 B22 − C12 )Qk (B2 , C1 , C2 ), for k = 2, 3, 4, 5, and Qk some polynomials with rational coefficients with degree 6(k − 2). Then, also from the factors of P2 we have the center families C1 and F4± . After all the above discussions we can restrict our parameters region to the set R = {a 6= 0, B2 6= 0, C1 6= 0, C2 6= 0, C2 6= AB2 }. The proof finishes showing that, computing some crossed resultants as in the previous studies, there are no centers in R. That is, the system of algebraic varieties S1 = {P1 = 0, · · · , P5 = 0}, defined in (10), has no solution. From (11) any point in S also is in S2 = {R2,3 = 0, R2,4 = 0, R2,5 = 0}. Straightforward computations show that the two by two resultants of the polynomials in S2 write as res(R2,3 , R2,k , C2 ) = a10k−8 A38k−42 (A2 + 1)3k−2 B240k−56 (a + A)4k (aA − 1)4k R3,k (a, A), for k = 4, 5, with R3,k polynomials with rational coefficients of degree 70k − 164. When a + A = 0 or aA − 1 = 0 we compute again the above resultant and we get similar expression but with R3,k polynomials with rational coefficients in a2 and degrees 21, 34 and 21, 32, respectively. When a + A 6= 0 and aA − 1 6= 0 we compute the last resultants res(R3,4 , R3,5 , a) =A2574 (16A2 + 1)(A2 + 16)(A2 + 9)4 (9A2 + 1)4 (A2 + 1)1362 × S382 (A2 )S4843 (A2 ), res(R3,4 , R3,5 , A) =a738 (25a2 + 1)2 (a2 + 25)2 (a2 + 9)16 (9a2 + 1)32 (a2 + 1)1283 × (1369a4 + 1102a2 − 171)2 (1038a4 − 171a2 − 17)2 × (321a4 + 48a2 + 31)8 (T6 (a2 ))2 (T̂6 (a2 ))6 T13 (a2 )(T191 (a2 ))2 T4843 (a2 ), for some polynomials with rational coefficients Sk , Tk , and T̂k of degree k. Conse- quently, all the possible values for a and A in S1 ∩ R are both algebraic. Using the same argument as in the first part of this proof, about the contradiction with the  ∗ ∗ fact that (A∗ )2 = e−πa = i2ia is transcendental, the set S1 ∩ R is empty. Proposition 3.3. System (2) has a center at infinity if the parameters a± , b± , c± belong to the families in Theorem 1.2. LIMIT CYCLES IN 4-STAR-SYMMETRIC PLS 9 Proof. All subsystems in families C1 , C2 , and C3± are linear centers. Family C1 has a global center because the subsystems defined in regions Σ± coincide. The singulari- ties of subsystems corresponding to family C2 are visible and located in the straight lines x − y = 0 and x + y = 0. Moreover, each one of them is time reversible with respect to the corresponding line. From both symmetries and using the definitions introduced in Section 2 we can check that the first return maps, in each quadrant, are all the identity. So, the composition is also the identity. The system associated to families C3± is time reversible with respect to the straight line y ± x = 0. For the parameters associated to families F1 and F3 , the full system (2) is time reversible with respect to y = 0 and x = 0, respectively. Then the infinity is a center because in a neighborhood of it the vector field is monodromic. The systems in family F4− can be moved to family F4+ changing θ to θ + π/2 as we have done in Lemma 2.1. Then, only remains the proofs of families F2 and F4+ . Both cases are proved similarly in three steps following the ideas of Section 2. First we compute the flying time maps, second we compute the return maps in each quadrant and finally we check that the full composition is the identity. In the family F2 , all the flying times are π/2 instead of the family F4+ that only two of them are. In this last family, we use also that the flying times in Σ+ are π/2 and for systems defined in Σ− the sum is π. As we have explained in the introduction we consider only a+ = a > 0. Family F2 can be expressed, as A + a 6= 0, also writing Aa − 1 ± c± = ∓ b . (12) A+a The return map in the first quadrant can be computed from the parametrized solu- tion of equation (2) in polar coordinates   b+ at Ab+ at (x(t), y(t)) = − + r e cos t, + r e sin t . (13) A+a A+a We notice that the equilibrium point belongs to the straight line Ax + y = 0 and, consequently, it is an invisible singularity for system (2) defined in Σ+ . Let τ0 , τ0 + τ1 , be the values of the time such that the solution starts at (ρ, 0) and ends at (0, ρ1 ). Substituting the values of cos τ0 and sin τ0 , from the equations x(τ0 ) = ρ and y(τ0 ) = 0, in the equation x(τ0 + τ1 ) = 0 we obtain that τ1 satisfies fρ (τ1 ) = 0 with b+ fρ (τ ) = ρ eaτ cos τ + (A eaτ sin τ + eaτ cos τ − 1). (14) A+a We observe that fρ (π/2) = 0 and, if this value is the smallest positive zero of fρ , we can obtain ρ b+ (A2 + 1) ρ1 = y(τ0 + π/2) = + . (15) A A(A + a) The solution in Σ+ defining Π1 can be seen in Figure 2 (right). The property fρ (τ ) 6= 0 for τ ∈ [0, π/2) follows from ρ (Aa − 1)b+ fρ (0) = ρ, fρ0 (0) = aρ + b+ , fρ0 (π/2) = − + , A A(A + a) 10 C. BUZZI, J. C. MEDRADO, AND J. TORREGROSA Σ− (0, ρ) (0, ρ1 ) Σ+ π 2 (ρ, 0) π 2 τ0 τ0 (ρ1 , 0) Figure 2. Representation of the return maps Π1 (right) and Π2 (left). and the fact that the derivative eaτ fρ0 (τ ) = eaτ (a cos τ − sin τ )ρ + b+ ((Aa − 1) sin τ + (A + a) cos τ ) A+a vanishes only at τ = 0 for a = 0 and has exactly one positive zero for a > 0. See Figure 3 (left). π 2 ρ ρ π 2 Figure 3. The graphs of fρ for ρ > 0 (left) and ρ < 0 (right), for different values of a > 0, a = 0 (dash). As ρ1 > 0, from the definition in Section 2 and (15), we get that the return map in the first quadrant is the affine map ρ b+ (A2 + 1) Π1 (ρ) = + . (16) A A(A + a) The study for the return map in the third quadrant coincides with above discus- sion, taking the solution passing through (ρ, 0) and (0, ρ1 ), with ρ < 0 and ρ1 < 0. Because the above procedure works also for negative values of ρ. See the plot of the function fρ for this case in Figure 3 (right). The definition of Π3 in Section 2 needs the positiveness of ρ and ρ1 . Changing adequately the signs in (15) we get ρ b+ (A2 + 1) Π3 (ρ) = − . (17) A A(A + a) The return maps in the second and fourth quadrant can be obtained similarly taking the solution of equation (2) passing through (0, ρ) and (ρ1 , 0). The equilibrium LIMIT CYCLES IN 4-STAR-SYMMETRIC PLS 11 point, that belongs to the straight line −Ax+y, is also invisible for system (2) defined in Σ− . Equations (13), (14), and (16) change to  −  b −at Ab− −at (x(t), y(t)) = + r e cos t, + r e sin t , A+a A+a b− gρ (τ ) = ρ e−aτ cos τ − (A e−aτ cos τ + e−aτ sin τ − A) A+a and b− (A2 + 1) ρ1 = y(τ0 + π/2) = −Aρ + . (18) A+a The solution in Σ− defining Π2 can be seen in Figure 2 (left). Using the definitions of the return maps in Section 2 but with distances, changing adequately the signs in (18), we get b− (A2 + 1) b− (A2 + 1) Π2 (ρ) = Aρ − , Π4 (ρ) = Aρ + . (19) A+a A+a The proof of this family finishes straightforward checking, using (16), (17), and (19), that the composition Π4 ◦ Π3 ◦ Π2 ◦ Π1 is the identity map. The study of the centers in family F4+ starts obtaining the equivalent relations A−a + Aa ∓ 1 + b− = b , c± = ∓ b . A+a A+a We notice that the definition of b+ coincides with the one in (12). Then, the return maps Π1 and Π3 , are the same as the previous family, defined in (16) and (17). Consequently, only the computation of Π2 and Π4 are necessary to be done. In a similar way, the solution of (2) in the second quadrant writes as   b+ −at Ab+ −at (x(t), y(t)) = − + r e cos t, + r e sin t . (20) A+a A+a The main difference with respect to family F2 is in the flying times from a point in x = 0 to y = 0, that they are not π/2. For this family the equilibrium point of system (2) defined in Σ− is visible and it coincides with the one, which is invisible, defined in Σ+ also for system (2). Let us consider the solution of (2), defined in the second quadrant, starting at (0, π1 ) and finishing at (π2 , 0), being τ1 the corresponding flying time. See Figure 4. Straightforward computations get e−aτ1 (A + a)ρ cos τ1 + b+ (A e−aτ1 sin τ1 + e−aτ1 cos τ1 + A2 ) = 0 (21) and 1 −aτ1 b+ π2 = − e ρ sin τ1 − (e−aτ1 sin τ1 − A e−aτ1 cos τ1 + A). A A(A + a) Using the function Π3 , defined in (17), we can compute the image, (0, π3 ), of the point (π2 , 0) with the third return map, 1 −aτ1 b+ π3 = − e ρ sin τ1 − (e−aτ1 sin τ1 − A e−aτ1 cos τ1 − A3 ). A2 A2 (A + a) 12 C. BUZZI, J. C. MEDRADO, AND J. TORREGROSA Σ− (0, π1 ) (0, π1 ) Σ+ (π2 , 0) (ρ, 0) π τ1 τ0 2 τ0 Figure 4. Representation of the return maps Π1 (right) and Π2 (left). . For the last step, we consider also the solution (20) from the point (0, π3 ) to (π4 , 0) with flying time τ2 . The simplified equivalent equation (14) for the flying τ2 , using that 2π is the total flying time and so, τ1 + τ2 = π, is (A + a) sin(τ1 − τ2 )ρ + b+ (sin(τ1 − τ2 ) − A cos(τ1 − τ2 ) − 2 e−aτ2 sin τ2 − A) = 0, (22) and 1 b+ π4 = (cos(τ1 − τ2 ) + 1)ρ + (A sin(τ1 − τ2 ) + cos(τ1 − τ2 ) + 2 e−aτ2 cos τ2 − 1). 2 2(A + a) Finally, we check, using the equations (21) and (22) and also that τ1 + τ2 = π, that π4 = ρ. As the composition of the four return maps is the identity, family F4+ defines a center at infinity as we wanted to prove.  4. Five limit cycles in a 4-star-symmetric system This section is devoted to the proof of Theorem 1.1. From Proposition 2.3 we can restrict our analysis to a+ + a− = 0. Under this condition, Proposition 4.2 will show that system (2) exhibits a weak focus of order 5 that unfolds 4 limit cycles. When this condition is broken an extra limit cycles appear and Theorem 1.1 follows. As in the previous section we denote a+ = −a− = a. We notice that with the change of variables (3), the infinity is moved to the origin. With the same procedure, we could obtain other results showing that there exist systems with weak foci of order k bifurcating k limit cycles for k = 2, 3, 4. From the proof of Proposition 3.2 it follows the next corollary. Corollary 4.1. There are no parameter values (a± , b± , c± ) such that system (4) has a weak focus of order higher than five at the origin. Proof. The proof follows checking all the cases appearing in the proof of Proposi- tion 3.2. More concretely, under the change of parameters (7), when B2 = 0 and C2 6= 0 the order of the weak focus is one. If B2 = C1 = 0 then the maximum weak focus order is four. Only when B2 6= 0 weak foci of order 5 appear. But when the first five Lyapunov quantities obtained in Lemma 3.1 vanish we get centers of (2).  Next result exhibits the maximum number of limit cycles bifurcating from the origin for system (4) under the condition a+ + a− = 0. LIMIT CYCLES IN 4-STAR-SYMMETRIC PLS 13 Proposition 4.2. Let us consider system (4) under the conditions a+ = −a− = a. Then, there exist values of the parameters a, b± , and c± such that system (4) has a weak focus at the origin of order 5 that unfolds 4 limit cycles. The proof follows fixing an specific value for one of the five free parameters and proving that the varieties {L1 = L2 = L3 = L4 = 0} intersects in, at least, one point in R4 . This is proved using the Poincaré–Miranda Theorem that we recall here by completeness. This technical result is a generalization of Bolzano’s Theorem for higher dimensions. Moreover, it is also necessary to prove that, at this point, the intersection is transversal and L5 is different from zero. Theorem 4.3 ([20]). Let f : Rn → Rn , f = (f1 , . . . , fn ), be a continuous map such that for each 1 ≤ i ≤ n, fi (Ii− ) ⊂ (−∞, 0] and fi (Ii+ ) ⊂ [0, ∞). Then there exists a point ζ ∈ I n such that f (ζ) = 0. Where I n = [−h, h]n , and Ii± = {(x1 , . . . , xn ) ∈ I n : xi = ±h}, for some positive value h. Next we will prove the main result of this section. Proof of Proposition 4.2. We start doing the change in the parameters space (7) as in Proposition 3.2 and fixing the value B2 = 10. Then, we have only four parameters (a, B1 , C1 , C2 ) and we define the function f = (L1 , L2 , L3 , L4 ). The proof is given in two steps. First, we provide a box I 4 in Theorem 4.3 such that, it exists a point, ζ ∗ = (a∗ , B1∗ , C1∗ , C2∗ ), satisfying f (ζ ∗ ) = 0. Second, we check that the next Lyapunov quantity and the determinant of the Jacobian matrix of f, at every point in the box I 4 , are different from zero. In particular, L5 (ζ ∗ ) 6= 0 and Df(a,B1 ,C1 ,C2 ) (ζ ∗ ) 6= 0. We will check such conditions first numerically and then analytically. More concretely, we do a computer assisted proof to prove the existence of a transversal intersection point by means of an accurate rational interval analysis. From Lemma 3.1 we compute the first Lyapunov quantities Lk , k = 1, . . . , 5. Then, as a+ + a− = 0, we have L0 = 0. As in the previous results, we denote a+ = −a− = a π and A = e− 2 a . Solving L1 = L2 = 0 in (9), with the condition C1 (10A − C2 ) 6= 0, we get C1 C2 AQ1 B1 = , C12 = , (23) 10A2 Q2 where Q1 = Q1 (a, A, C2 ) and Q2 = Q2 (a, A, C2 ) are the polynomials Q1 =102 (3A + 4a)C23 − 103 A(9A2 a + 8Aa2 − 3A − a)C22 + 104 A2 a(3A2 + Aa + 1)C2 + 105 A4 a2 , Q2 =(−4Aa + 3)C23 − 10(A2 a + 8Aa2 − 3A − 9a)C22 − 102 Aa(A2 − Aa + 3)C2 + 103 A3 a2 . Clearly necessary conditions are that Q2 in (23) should be different from zero and the quotient Q1 /Q2 be positive. Using the above relations we get L3 = C1 C2 aA−4 (A2 + 1)(10A − C2 )3 Q−2 2 P3 (a, A, C2 ), L4 = C1 C2 aA−7 (A2 + 1)(10A − C2 )3 Q−3 2 P4 (a, A, C2 ), (24) −10 L5 = C1 C2 aA 2 (A + 1)(10A − C2 ) 3 Q−4 2 P5 (a, A, C2 ), 14 C. BUZZI, J. C. MEDRADO, AND J. TORREGROSA where Pk are polynomials with rational coefficients of degrees 16, 28, 40 for k = 3, 4, 5, respectively. We do not write here the expressions of such polynomials be- cause of their size. In fact, they have 104, 489, and 1313 monomials, respectively. Consequently, we can use (L1 , L2 , P3 , P4 , P5 ) instead of (L1 , L2 , L3 , L4 , L5 ) and the equivalent system of equations Se = {L1 = L2 = P3 = P4 = 0}. The other factors in (24) are not considered because they vanish identically all the equations. As A = e−aπ/2 , then a weak focus point of order five is given, from (24), by a solution (a∗ , C2∗ ) of the system of two equations {P3 = P4 = 0} (25) ∗ with only two unknowns (a, C2 ) such that P5 (a , A , C2∗ ) 6= 0 ∗ and the necessary ∗ ∗ ∗ ∗ ∗ ∗ ∗ −a∗ π/2 condition Q1 (a , A , C2 )Q2 (a , A , C2 ) > 0, where A = e . From (23) we get the corresponding values for C1∗ and B1∗ . Straightforward computations show that the Jacobian matrix of the function P : (a, B1 , C1 , C2 ) 7→ (L1 , L2 , P3 , P4 ) is DP(a∗ , B1∗ , C1∗ , C2∗ ) = 2A2 C1∗ Q2 (a∗ , A∗ , C2∗ )D(P3 , P4 )(a,C2 ) (a∗ , A∗ , C2∗ ). Then, the transversality of the varieties defined in system Se is a direct consequence of the nonvanishing Jacobian J = D(P3 , P4 )(a,C2 ) (a∗ , A∗ , C2∗ ) 6= 0. We point out that the above function is a polynomial in (a, A, C2 ) of degree 43 with 3180 monomials. In Figure 5 we have drawn the zero level curves {P3 = 0} and {P4 = 0} in different zones of the reduced parameter space (a, C2 ). It is shown also some zooms to see better, at least, one transversal intersection satisfying the necessary conditions for the existence of a transversal intersection point in the complete parameter space (a, B2 , C1 , C2 ), together with the condition Q1 Q2 > 0, Figure 5. The drawing of the zero level curves {P3 = 0} and {P4 = 0} in red and blue, respectively, in the parameter space (a, C2 ) Numerically, working with enough accuracy, we can get (a∗ , A∗ , B1∗ , C1∗ , C2∗ ) ≈ (0.322766, 0.602300, 8.496686, 33.902075, 0.909179). (26) Additionally, Q1 (a∗ , C2∗ ) = 4.31954 · 103 , Q2 (a∗ , C2∗ ) = 2.26359, (27) P5 (a∗ , C2∗ ) = −2.03342 · 1010 , J(a∗ , C2∗ ) = 3.30868 · 1012 . We remark that the chosen point also satisfies the necessary conditions such that B1 and C1 , defined in (23), exist. This is not true for all intersection points. LIMIT CYCLES IN 4-STAR-SYMMETRIC PLS 15 We can recover the original parameters, with the inverse change of (7), a(B1 ∓ B2 ) + C1 ∓ C2 ∓a(C1 + C2 ) ± B1 − B2 b± = , c± = . 2(a2 + 1) 2(a2 + 1) Then, the origin of system (4) for (a∗ , b+ + − − ∗ , c∗ , b∗ , c∗ ) is a weak focus of order 5. Substituting (26) in the above expression we get (a∗ , b+ + − − ∗ , c∗ , b∗ , c∗ ) ≈ (0.3227657, 14.72030, −5.502868, 18.46683, −3.287880). (28) So, at this point, numerically the first four Lyapunov quantities vanish transversally and L5 (a∗ , b+ + − − 8 ∗ , c∗ , b∗ , c∗ ) ≈ −2.0197688 · 10 . Next, we prove analytically the existence of the point (26), equivalently (28), satisfying all the conditions described above. We start writing (eaπ/2 )−1 instead e−aπ/2 in all the expressions involved in the computations because they are simpler when the exponential function is increasing. In order to apply Theorem 4.3 to system (25) we perform the change of variables 363394237984 14774107617 53184397 a= − u− v, 1125876185023 61873882348312 13103427838970537 (29) 2250212532359 2109559243 20674157 v C2 = − u− v. 2474994238126 18668085007969 6487787928613874 Then, we get P3 (u, v) ≈ u + O2 (u, v) and P4 (u, v) ≈ v + O2 (u, v). With this change of variables, the involved functions, P3 , P4 , P5 , Q1 , Q2 , and J, are all poly- nomials with rational coefficients in the variables π, u, v, and eπ(x0 +x1 u+x2 v) , with x0 = −181697118992/1125876185023, x1 = 14774107617/123747764696624, and x2 = 53184397/26206855677941074. The existence of an intersection point, (u∗ , v ∗ ), follows finding an appropriate value for h satisfying the hypotheses of Theorem 4.3. For h = 10−7 , it can be checked numerically that max P3 (−h, v) < 0, min P3 (h, v) > 0 for v ∈ [−h, h], (30) max P4 (u, −h) < 0, min P4 (u, h) > 0 for u ∈ [−h, h]. Additionally, min Q1 (u, v) > 0, min Q2 (u, v) > 0, max P5 (u, v) < 0, min J(u, v) > 0, (31) for (u, v) ∈ [−h, h]2 . We can recover the existence of the point (26) from (29) and (23). The transversality and the existence of a weak focus of order 5, at least numerically, is also clear. The analytical proof of the above properties can be done using arithmetic interval analysis with h = 10−9 . We will use the technical results of the last section with the polynomial approximation for the exponential function given by m = 21, see Proposition 5.1. For π, we take 21053343141 1783366216531 π= ≤π≤ =π 6701487259 567663097408 from its continuous fraction expression with 22 and 23 partial quotients. First we prove (30) and second (31). With the change (29) we write − − + + P3 (−h, v) = g1 (v, eπ(y0 +y1 v) ) =: ĝ1 (v), P3 (h, v) = g2 (v, eπ(y0 +y1 v) ) =: ĝ2 (v), − − + + P4 (u, −h) = g3 (u, eπ(z0 +z1 u) ) =: ĝ3 (u), P4 (u, h) = g4 (u, eπ(z0 +z1 u) ) =: ĝ4 (u), 16 C. BUZZI, J. C. MEDRADO, AND J. TORREGROSA with y0− = x0 − x1 h, y1− = x2 , y0+ = x0 + x1 h, y1+ = x2 , z0− = x0 − x2 h, z1− = x1 , z0+ = x0 + x2 h, and z0+ = x1 all rational numbers. We do not detail the explicit expressions of the functions g1 (x, y), g2 (x, y), g3 (x, y), and g4 (x, y) because of their size. They are polynomials with rational coefficients of degrees 16, 16, 28, and 28 with 118, 118, 332, and 332 monomials in (x, y), respectively. Then applying Proposition 5.2 we have that ĝ1 ([−h, h]) ⊂ [η1 , ξ1 ] ≈ [−1.0084558 · 10−9 , −9.9154416 · 10−10 ], ĝ2 ([−h, h]) ⊂ [η2 , ξ2 ] ≈ [9.9154416 · 10−10 , 1.0084558 · 10−9 ], ĝ3 (−h) ∈ [η3 , ξ3 ] ≈ [−1.09345834416 · 10−9 , −9.06375115265 · 10−10 ], ĝ3 (−h)+2hĝ30 (−h) ∈ [η4 , ξ4 ] ≈ [−1.09345830474 · 10−9 , −9.06375075847 · 10−10 ], ĝ300 ([−h, h]) ⊂ [η5 , ξ5 ] ≈ [−19.7094899356371570, −19.7091168953367466], ĝ4 (−h) ∈ [η6 , ξ6 ] ≈ [9.065416558399 · 10−10 , 1.093624884733879 · 10−9 ], ĝ4 (h) ∈ [η7 , ξ7 ] ≈ [9.065416558414 · 10−10 , 1.093624884735517 · 10−9 ], ĝ400 ([−h, h]) ⊂ [η8 , ξ8 ] ≈ [−19.7094899356371588, −19.7091168953367484]. (32) The above expressions prove (30). We remark that for ĝ3 and ĝ4 we have used the extreme values of the interval of definition, the sign of the second derivative and the value of the straight line tangent at x = −h evaluated at x = h. This is because the lower and upper bounds obtained with the results in the last section for ĝ3 ([−h, h]) and ĝ4 ([−h, h]) do not prove directly (30). Next we prove (31). As above, with the change (29) we write Q1 (u, v) = g5 (u, v, eπ(x0 +x1 u+x2 v) ), Q2 (u, v) = g6 (u, v, eπ(x0 +x1 u+x2 v) ), P5 (u, v) = g7 (u, v, eπ(x0 +x1 u+x2 v) ), J(u, v) = g8 (u, v, eπ(x0 +x1 u+x2 v) ), where g5 (x, y, z), g6 (x, y, z), and g7 (x, y, z) are polynomials with rational coefficients of degrees 6, 5, and 40 and 56, 41, 8494 monomials, respectively. The last g8 (x, y, z) is also a polynomial, of degree 43 and 20695 monomials, but the coefficients are polynomials of degree 1 in π also with rational coefficients. Applying Proposition 5.2 we have Q1 ([−h, h]2 ) ⊂ [η9 , ξ9 ] ≈ [4.3195435234081 · 103 , 4.3195435234277 · 103 ], Q2 ([−h, h]2 ) ⊂ [η10 , ξ10 ] ≈ [2.263592961126, 2.263592961413], P5 ([−h, h]2 ) ⊂ [η11 , ξ11 ] ≈ [−2.03344613 · 1010 , −2.03340080 · 1010 ], J([−h, h]2 ) ⊂ [η12 , ξ12 ] ≈ [3.30867903 · 1012 , 3.30868520 · 1012 ]. We notice that the numerical values (27) agree with the above approximations. Finally, we remark that with higher values for h we have not been able to prove the properties (30) and (31). Moreover, with smaller values for h the lower and upper bounds for the images of ĝ3 ([−h, h]) and ĝ4 ([−h, h]) obtained with Proposition 5.2 do not prove the sign properties (30) for P4 . This is why we have used all properties (32)  5. Appendix. Bounding a exponential-polynomial functions Let us denote u = (u1 , . . . , un ), b = (b1 , . . . , bn ), and ` = (`0 , `1 , . . . , `n , `n+1 ). In this section we provide rational lower and upper bounds for polynomials in the LIMIT CYCLES IN 4-STAR-SYMMETRIC PLS 17 variables u, π, and eπ(a+b·u) with rational coefficients when a < 0 and bi ≥ 0. That is, each monomial can be written as A` M` (u) = A` eπ`0 (a+b·u) u`11 · · · u`nn π `n+1 , P n where b·u = bi ui . We notice that u ∈ [−h, h]n for a rational number h. Moreover, i=1 the value a, the coefficients A` and the components of b are all rational numbers and the components of ` are nonnegative integer numbers. The following straightforward results provide rational lower and upper bounds for the principal functions in the above monomial and for the complete polynomial when n = 0, 1, and 2. They can be easily implemented in any computer algebra system. Proposition 5.1. Let m be a natural number, h > 0 be a rational number, and π, π be two rational numbers such that π < π < π. Then the next properties hold: (i) For −1 ≤ aπ ≤ 0, we have E(a) := (Em (−πa, 3))−1 ≤ eπa ≤ (Em (−πa, 0))−1 =: E(a), m xj P xm+1 where Em (x, τ ) = +τ . j=0 j! (m + 1)! (ii) When M` = eπ`0 a π `1 and −1 ≤ aπ ≤ 0, then M` := (E(a))`0 π `1 ≤ M` ≤ (E(a))`0 π `1 =: M` . (iii) When M` (u) = eπ`0 (a+b1 u1 ) u1`1 π `2 , −1 ≤ (a + b1 h)π ≤ 0, −1 ≤ (a − b1 h)π ≤ 0, and `1 − h`0 b1 > 0 then M` (u) ≤ M` (u) ≤ (E(a + b1 h))`0 h`1 π `2 =: M` (u), where   ` ` (E(a − b1 h)) 0 π 2 `1 = 0, M` (u) = −(E(a − b1 h))`0 h`1 π `2 `1 odd,  0 `1 > 0, even. (iv) When M` (u) = eπ`0 (a+b1 u1 +b2 u2 ) u`11 u`22 π `3 , −1 ≤ (a + b1 h + b2 h)π ≤ 0, −1 ≤ (a − b1 h − b2 h)π ≤ 0, `1 − h`0 b1 > 0 and `2 − h`0 b2 > 0 then M` (u) ≤ M` (u) ≤ (E(a + b1 h + b2 h))`0 h`1 +`2 π `3 =: M` (u), where    (E(a − b1 h − b2 h))`0 π `3 `1 = `2 = 0,     (`1 , `2 both odd and b2 ≥ b1 )  −(E(a − (b1 − b2 )h))`0 h`1 +`2 π `3 or (`1 odd , `2 even), M` (u) =   (`1 , `2 both odd and b1 ≥ b2 )  −(E(a − (b2 − b1 )h))`0 h`1 +`2 π `3   or (`1 even , `2 odd),   0 `1 > 0, `2 > 0, both even. Proposition 5.2. Let A` be rational numbers and M` and M` (u) be monomials, satisfying the hypotheses of Proposition 5.1. If a < 0 and bi ≥ 0, for i = 1, . . . , n, then there exists small enough rational number h > 0 such that the next properties hold. 18 C. BUZZI, J. C. MEDRADO, AND J. TORREGROSA P d (i) If F = A` M` where M` = eπ`0 a π `1 , `=0 " d d # X X F ∈ [η, ξ] := A` L(A` , M` ), A` U(A` , M` ) . `=0 `=0 P d (ii) If F (u) = A` M` (u) where M` (u) = eπ`0 (a+b·u) u`11 · · · u`nn π `n+1 and n = 1 or `=0 n = 2, then " d d # X X F ([−h, h]n ) ⊂ [η, ξ] := A` L(A` , M` (u)), A` U(A` , M` (u)) . `=0 `=0 Here, for every 0 < x < x < x, being x, x rational numbers, we have    x if A > 0,  x if A > 0, L(A, x) = 0 if A = 0, and U(A, x) = 0 if A = 0,  x if A < 0,  x if A < 0. We notice that h is chosen in such way that the maximal and minimal values of each monomial are taken in the vertices of [−h, h]n . Acknowledgements This work has received funding from the Ministerio de Economı́a, Industria y Competitividad - Agencia Estatal de Investigación (FEDER Grant MTM2016-77278- P), the Agència de Gestió d’Ajuts Universitaris i de Recerca (Grant 2017 SGR 1617), the Brazilian agencies FAPESP (Grants 2013/24541-0 and 2017/03352-6), CAPES (Grant PROCAD 88881.068462/2014-01), CNPq (Grant 308006/2015-1), FAPEG (29199/2018), and FAPEG/CNPq (Grant PRONEX 2017-10267000508), and the European Union’s Horizon 2020 research and innovation programme (Grant Dynamics-H2020-MSCA-RISE-2017-777911). References [1] V. Acary, O. Bonnefon, and B. Brogliato. Nonsmooth modeling and simulation for switched circuits, 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., 1973. [3] A. A. Andronov, A. A. Vitt, and S. E. Khaı̆kin. Theory of oscillators. Dover Publications, Inc., New York, 1987. Translated from the Russian by F. Immirzi, Reprint of the 1966 translation. [4] M. Apri, N. Banagaaya, B. van den J. B, R. Brussee, D. Bourne, S. Fatima, F. Irzal, J. Rademacher, B. Rink, F. Veerman, S. Verpoort, and J. Hulshof. Analysis of a model for ship maneuvering. In Proceedings of the 79 th European Study Group Mathematics with Industry, pages 83–116. Vrije Univ. Amst., 2011. [5] C. A. Buzzi, T. de Carvalho, and M. A. Teixeira. Birth of limit cycles bifurcating from a nonsmooth center. J. Math. Pures Appl. (9), 102(1):36–47, 2014. [6] C. A. Buzzi, C. Pessoa, and J. Torregrosa. Piecewise linear perturbations of a linear center. Discrete Contin. Dyn. Syst., 33(9):3915–3936, 2013. [7] P. T. Cardin and J. Torregrosa. Limit cycles in planar piecewise linear differential systems with nonregular separation line. Phys. D, 337:67–82, 2016. [8] C. Chicone and M. Jacobs. Bifurcation of critical periods for plane vector fields. Trans. Amer. Math. Soc., 312(2):433–486, 1989. LIMIT CYCLES IN 4-STAR-SYMMETRIC PLS 19 [9] C. Christopher. Estimating limit cycle bifurcations from centers. In Differential equations with symbolic computation, Trends Math., pages 23–35. Birkhäuser, Basel, 2005. [10] B. Coll, R. Prohens, and A. Gasull. The center problem for discontinuous Liénard differ- ential equation. Internat. J. Bifur. Chaos Appl. Sci. Engrg., 9(9):1751–1761, 1999. Discrete dynamical systems. [11] M. di Bernardo, C. J. Budd, A. R. Champneys, and P. Kowalczyk. Piecewise-smooth dynamical systems. Theory and applications, volume 163 of Applied Mathematical Sciences. Springer- Verlag London, Ltd., London, 2008. [12] F. Dumortier, J. Llibre, and J. C. Artés. Qualitative theory of planar differential systems. Universitext. Springer-Verlag, Berlin, 2006. [13] J. Fan, S. Xue, and G. Chen. On discontinuous dynamics of a periodically forced double-belt friction oscillator. Chaos Solitons Fractals, 109:280–302, 2018. [14] A. F. Filippov. Differential equations with discontinuous righthand sides, volume 18 of Math- ematics and its Applications (Soviet Series). Kluwer Academic Publishers Group, Dordrecht, 1988. [15] E. Freire, E. Ponce, F. Rodrigo, and F. Torres. Bifurcation sets of continuous piecewise linear systems with two zones. Internat. J. Bifur. Chaos Appl. Sci. Engrg., 28(11):2073–2097, 1998. [16] 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. [17] E. Freire, E. Ponce, and F. Torres. A general mechanism to generate three limit cycles in planar Filippov systems with two zones. Nonlinear Dynam., 78(1):251–263, 2014. [18] M. Han. Liapunov constants and Hopf cyclicity of Liénard systems. Ann. Differential Equa- tions, 15(2):113–126, 1999. [19] S. Huan and X.-S. Yang. Generalized Hopf bifurcation emerged from a corner in general planar piecewise smooth systems. Nonlinear Anal., 75(17):6260–6274, 2012. [20] W. Kulpa. The Poincaré-Miranda theorem. Amer. Math. Monthly, 104(6):545–550, 1997. [21] J. Llibre and E. Ponce. Three nested limit cycles in discontinuous piecewise linear differen- tial systems with two zones. Dyn. Contin. Discrete Impuls. Syst. Ser. B Appl. Algorithms, 19(3):325–335, 2012. [22] J. C. Lucero and C. A. Gajo. Oscillation region of a piecewise-smooth model of the vocal folds. Commun. Math. Sci., 4(2):453–469, 2006. [23] R. Lum and L. O. Chua. Global properties of continuous piecewise-linear vector fields. Part I: Simplest case in R2 , 1990. Memorandum UCB/ERL M90/22, University of California at Berkeley. [24] J. C. Medrado and J. Torregrosa. Uniqueness of limit cycles for sewing planar piecewise linear systems. J. Math. Anal. Appl., 431(1):529–544, 2015. [25] I. Niven. Irrational numbers. The Carus Mathematical Monographs, No. 11. The Mathematical Association of America. Distributed by John Wiley and Sons, Inc., New York, N.Y., 1956. [26] V. G. Romanovski and D. S. Shafer. The center and cyclicity problems: a computational algebra approach. Birkhäuser Boston, Inc., Boston, MA, 2009. [27] R. Roussarie. Bifurcations of planar vector fields and Hilbert’s sixteenth problem. Modern Birkhäuser Classics. Birkhäuser/Springer, Basel, 1998. [28] Y. Zou and T. Küpper. Generalized Hopf bifurcation emanated from a corner for piecewise smooth planar systems. Nonlinear Anal., 62(1):1–17, 2005. Departamento de Matemática, Universidade Estadual Paulista, 15054-000 São José do Rio Preto, Brazil Email address: claudio.buzzi@unesp.br Instituto de Matemática e Estatı́stica, Universidade Federal de Goiás, 74001-970 Goiânia, Goiás, Brazil Email address: medrado@ufg.br Departament de Matemàtiques, Universitat Autònoma de Barcelona, 08193 Bel- laterra, Barcelona, Spain Email address: torre@mat.uab.cat