Physica D 371 (2018) 28–47 Contents lists available at ScienceDirect Physica D journal homepage: www.elsevier.com/locate/physd Limit cycles via higher order perturbations for some piecewise differential systems Claudio A. Buzzi a, Maurício Firmino Silva Lima b, Joan Torregrosa c,* a Departamento de Matemática, Universidade Estadual Paulista, São José do Rio Preto, Brazil b Centro de Matemática Computação e Cognição. Universidade Federal do ABC, 09210-170. Santo Andre. S.P., Brazil c Departament de Matemàtiques, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Spain h i g h l i g h t s • The number of limit cycles bifurcating from linear centers. • Piecewise perturbations in two zones with a straight line of separation. • Upper bound for that number for arbitrary perturbation order. • The upper bounds are reached for perturbations of order one and two. • Some applications to concrete families. a r t i c l e i n f o Article history: Received 21 August 2017 Received in revised form 20 December 2017 Accepted 15 January 2018 Available online 31 January 2018 Communicated by J. Dawes Keywords: Non-smooth differential system Limit cycle in Melnikov higher order perturbation Liénard piecewise differential system a b s t r a c t A classical perturbation problem is the polynomial perturbation of the harmonic oscillator, (x′, y′) = (−y + εf (x, y, ε), x + εg(x, y, ε)). In this paper we study the limit cycles that bifurcate from the period annulus via piecewise polynomial perturbations in two zones separated by a straight line. We prove that, for polynomial perturbations of degree n, nomore thanNn−1 limit cycles appear up to a study of orderN . We also show that this upper bound is reached for orders one and two. Moreover, we study this problem in some classes of piecewise Liénard differential systems providing better upper bounds for higher order perturbation in ε, showing also when they are reached. The Poincaré–Pontryagin–Melnikov theory is the main technique used to prove all the results. © 2018 Elsevier B.V. All rights reserved. 1. Introduction In last decades, piecewise differential systems have been useful for modeling real processes and different modern devices. For simplicity also the linear piecewise differential systems provide adequate models with very accurate results close to the observed data. See for more details [1,2]. Although in recent years these systems have attracted a good deal of attention, the first stages of modeling with piecewise systems started with Andronov and coworkers, see [3]. In this paper, we study the number of isolated periodic orbits, the so called limit cycles, of a polynomial piecewise perturbation of degree n of a linear center, when the separation curve is a straight line Σ passing through the center. That is, we consider systems * Corresponding author. E-mail addresses: buzzi@ibilce.unesp.br (C.A. Buzzi), mauricio.lima@ufabc.edu.br (M.F.S. Lima), torre@mat.uab.cat (J. Torregrosa). written as (x′, y′) = ( − ∂H ∂y + N∑ i=1 εif ± i (x, y), ∂H ∂x + N∑ i=1 εig± i (x, y) ) , (1) such that the unperturbed system has a center at (xc, yc) with H a quadratic polynomial and the perturbations f ± i and g± i are polynomials of degree n defined in each side of Σ . Via an affine change of coordinates, if necessary, it is not restrictive to assume that we are studying the piecewise polynomial perturbation of the harmonic oscillator. Consequently, we consider the above system in the form (x′, y′) = ( −y + N∑ i=1 εif ± i (x, y), x + N∑ i=1 εig± i (x, y) ) , (2) with f ± i and g± i polynomials of degree n, defined inΣ± 0 = {(x, y) ∈ R2 : ±y ≥ 0}. This problem, for non piecewise perturbations, was studied by Iliev in [4], proving that the number of limit cycles is https://doi.org/10.1016/j.physd.2018.01.007 0167-2789/© 2018 Elsevier B.V. All rights reserved. https://doi.org/10.1016/j.physd.2018.01.007 http://www.elsevier.com/locate/physd http://www.elsevier.com/locate/physd http://crossmark.crossref.org/dialog/?doi=10.1016/j.physd.2018.01.007&domain=pdf mailto:buzzi@ibilce.unesp.br mailto:mauricio.lima@ufabc.edu.br mailto:torre@mat.uab.cat https://doi.org/10.1016/j.physd.2018.01.007 C.A. Buzzi et al. / Physica D 371 (2018) 28–47 29 bounded by [N(n − 1)/2]. In this paper, we prove by induction on the degree n that, in the piecewise case, the upper bound is approximately doubled. In the full paper [·] denotes the integer part function. The main tool used in this paper will be the Poincaré– Pontryagin–Melnikov technique. In this theory the limit cycles bifurcate from the level sets of a Hamiltonian H.More concretely, from the closed level sets that do not have equilibrium points. For system (2), the limit cycles γε(r) bifurcate from the level sets ϕr = {x2 + y2 = r2}. That is, γε(r) tends to ϕr when ε goes to 0. Here, this type of periodic orbits, with r > 0 and when they are isolated, are called limit cycles bifurcating from the period annulus of the center. It is clear that, for ε small enough, this type of solutions of (2) are well defined in the complement set of any neighborhood of the origin. We can say that these periodic orbits are away from the origin. Then, the Filippov convention is not necessary here, see [2]. In fact, there are other mechanism, like averaging method, to study the limit cycles bifurcating from the closed level sets of a Hamiltonian. But, in the plane and for an autonomous perturbation, they are equivalent, see [5] or [6]. In Section 2we recall the necessary results for proving ourmain goals. Next result provides upper bounds for the number of limit cy- cles bifurcating from the period annulus of the harmonic oscillator with piecewise polynomials of degree n. Theorem 1.1. The maximum number of limit cycles of system (1) bifurcating from the unperturbed period annulus is at most n up to order 1 and at most Nn − 1 up to order N ≥ 2, n ≥ 1. Moreover, there exist perturbation parameters such that for orders 1 and 2 the upper bounds are reached. The Poincaré–Pontryagin–Melnikov technique provides a func- tion, MN , for each perturbation order and the isolated periodic orbits are given from the simple zeros of it. The study of the number of zeros cannot be done in general because of the difficulties to find explicitly the functionsMN . There are few papers dealing with these objectives, see for example [4,7,8]. The upper bounds for the number of limit cycles up to order N are given by the number of simple zeros of MN and, from the upper bound presented in the above result, it seems that this number always increases with N. This is not the case for degree n = 1 perturbations in piecewise continuous and piecewise sewing, see [9] and [10], respectively. The study up to order seven and for general piecewise perturbation of degree 1 is done in [11], where it is proved that the increasing sequence of zeros is 1, 1, 2, 3, 3, 3, 3. Quadratic and cubic pertur- bations up to order 5 are studied in [12], where the authors prove that there exist polynomial perturbations such thatMN has 2, 3, 5, 6, and 8 simple zeros for n = 2 and 3, 5, 8, 11, and 13 for n = 3. Another example exhibiting 3 limit cycles for the quadratic family is given in [13]. It is clear that the study about how the number of limit cycles increaseswith the perturbation order is far to be solved, even in the case of small degrees. Only for some special families this problem can be solved. For example, if we consider the non-piecewise perturbation of the harmonic oscillator in the Liénard family, (x′, y′) = ( −y + N∑ i=1 εifi(x), x ) , with fi polynomials of degree n. In [8] it is proved that all the functions MN has the same number, [(n − 1)/2], of simple zeros asM1. For the classical Liénard family (x′, y′) = (−y+ f (x), x), Lins, de Melo and Pugh conjectured in [14] that the maximal number of limit cycles is exactly [(n − 1)/2]. Zuppa in [15] proved this conjecture for the limit cycles bifurcating from the origin via a degenerate Hopf bifurcation. The conjecture is also true for degree n = 4, see [16], but it fails for higher degrees, see [17] and [18]. As in the smooth case, the upper bounds will not be reached in general. But this is a question not solved yet even for the smooth case. The general study of the zeros of the Melnikov functions up to every order is also very difficult to be done for a polynomial family of degree n. So, it is quite natural to restrict our attention to study only some special families under piecewise perturbation. In what follows, we present some results about the number of limit cycles in some piecewise Liénard families. This kind of problems were also treated in [19] using Lyapunov quantities and in [20] for rational Liénard systems. Recently, Sheng in [21] usesMelnikov functions to study also the number of limit cycles for other piece- wise Liénard families. In this work line we present our main results. First, in Theo- rem 1.2, we deal with a Liénard family where the upper bounds given in the above theorem are decreased and achieved for pertur- bations of first and second orders. Next, adding some symmetry hypothesis, we can prove, in Theorem 1.3, that no limit cycles appear or the number of limit cycles up to any order of perturbation coincides with the ones bifurcating from the first order. Finally, changing the discontinuity line, from the x-axis to the y-axis, in Theorem 1.4, we prove that the Melnikov functions of any order have the same aspect and the upper bound is reached always. In this case we have proved that the number of limit cycles does not increase with the order of perturbation. Consider now, (x′, y′) = ( −y + N∑ i=1 εif ± i (x), x ) , (3) defined inΣ± 0 = {(x, y) ∈ R2 : ±y ≥ 0},where f ± i are polynomials of degree n. Theorem 1.2. The maximum number of limit cycles of system (3) bifurcating from the unperturbed period annulus is atmost [(n−1)/2] and n + [(n − 1)/2] for orders 1 and 2, respectively. Moreover, there exist polynomials f ± i such that for orders 1 and 2 these upper bounds are reached. Theorem 1.3. (a) When f ± i are even polynomials, then system (3) has a center at the origin for every ε. (b) When f ± i are odd polynomials, the maximum number of limit cycles of system (3) bifurcating from the period annulus is at most (n−1)/2 up to any order of perturbation. Moreover, there exist polynomials f ± i such that this upper bound is reached. Last result deals with system (3) but changing the separation line to the y-axis. In this case, we consider the system (x′, y′) = ( −y + N∑ i=1 εif ± i (x), x ) , (4) defined inΣ± 1 = {(x, y) ∈ R2 : ±x ≥ 0},where f ± i are polynomials of degree n. Theorem 1.4. The maximum number of limit cycles of family (4) bifurcating from the period annulus is at most n up to any order of perturbation. Moreover, there exist polynomials f ± i such that the upper bound is reached. The paper is organized as follows. In Section 2 we introduce the main tools for proving all the results. Section 3 is devoted to the proof of Theorem 1.1, first proving the upper bounds and second providing explicit examples exhibiting that number of limit cycles. The upper bounds are proved by induction on the degree of the perturbation. Sections 4 and 5 deal with the proofs of 30 C.A. Buzzi et al. / Physica D 371 (2018) 28–47 Theorems 1.2, 1.3, and 1.4. Finally, in Section 6 we study higher order perturbations for some fixed families, but for small values of n. First, we show that the number of limit cycles provided by Theorem 1.2 does not increase up to order 4. This is done in Proposition 6.1. Second, for a generalized Liénard family and up to order 5, we show that the number of limit cycles increases with the perturbation order but not as Theorem 1.1 says. This is done in Proposition 6.2. The conclusions of these results suggest that, when we fix the degree of the polynomial perturbation, the number of limit cycles will increase with the order of the perturbation up to a maximal number. We can say that there is a saturation in the growth of the number of limit cycles. This phenomenon can be observed in piecewise linear systems with a straight line of separation (see [11]),when the separation is non-regular (see [22]), and also in analytical families (see [8]). 2. The return map and the Poincaré–Pontryagin–Melnikov functions This section is devoted to present the main tools that we need to state and prove the results of this paper. We closely follow the presentation in [23], which decomposes an arbitrary one-form in polar coordinates. It is based on the decompositions introduced by Françoise in [24–26]. This procedure is also used in [11], but only for piecewise linear differential systems. The vector field X given in (1) with the separation line Σ0 = {(x, 0) : x ∈ R}, in polar coordinates (x, y) = (r cos θ, r sin θ ), can be written as (ṙ, θ̇ ) = ⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ ( N∑ i=1 εiR+ i (r, θ ), 1 + N∑ i=1 εiΘ+ i (r, θ ) ) if θ ∈ [0, π ),( N∑ i=1 εiR− i (r, θ ), 1 + N∑ i=1 εiΘ− i (r, θ ) ) if θ ∈ [π, 2π ), (5) where R± i ,Θ ± i are analytic functions in r , sin θ and cos θ . The above system can also be expressed as⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩ dH + N∑ i=1 εiω+ i = 0 if θ ∈ [0, π ), dH + N∑ i=1 εiω− i = 0 if θ ∈ [π, 2π ), (6) where H(r) = (x2 + y2)/2 = r2/2, and ω± i = ω± i (r, θ ) are analytic one-forms, 2π-periodic in θ and polynomial in r. Wedenote by X± each vector field (5) defined inΣ± 0 = {(x, y) ∈ R2 : ±y ≥ 0}. Clearly, we need to define X in the separation line Σ0.We say that a solution of the above vector field, X±, is of sewing type if when it crosses the separation line,Σ0, satisfies that X+ · (0, 1) and X− · (0, 1) have the same sign in the intersection points with Σ0.When this is not the case we can use the Filippov convention, see [2], to define the vector field over such special points in Σ0. These points define the so called sliding set. In this paper we are dealing with isolated periodic orbits bifurcating from a center. So, by continuity, we are interested only in the so called limit cycles of sewing type. The periodic orbits that cross the sliding set will be studied in future works. This is the reason why we do not detail more such special solutions. We define straightaway the main tool used in this paper: the Poincaré map. Let r+(θ, ρ) (resp. r−(θ, ρ)) be the solution of X such Fig. 1. (a) Return map and (b) half-return maps π+ X and (π− X )−1 = π+ R(X) of system (6). that r+(0, ρ) = ρ (resp. r−(π, ρ) = ρ). Then,we define the positive Poincaré half-return map as π+ X (ρ) = r+(π, ρ) and the negative Poincaré half-return map as π− X (ρ) = r−(2π, ρ). The complete Poincaré return map associated to X is given by the composition of these two maps πX (ρ) = π− X (π+ X (ρ)), see Fig. 1(a). We can also write them, for system (5), as a power series in ε as π± X (ρ, ε) = ρ + N∑ i=1 εip± i (ρ), and πX (ρ, ε) = ρ + N∑ i=1 εipi(ρ). It is important tomention that bothmaps, π+ X and π− X , are analytic in ρ and ε, so the complete Poincaré map associated to X is also analytic. The result below relates the half-return map on Σ− 0 of the vector field X with the half-return map onΣ+ 0 of the transformed vector field R(x) which is defined by the reversibility property. That is, for each solution γ (t) = (x(t), y(t)) of X then R(γ (t)) = (x(−t),−y(−t)) is a solution of R(X). The proof follows directly from the reversibility property. Moreover the study of the number of zeros of the displacement map, πX (ρ) − ρ, is easier. Lemma 2.1. The first non-zero term of the map πX (ρ) − ρ coincides with the first non-zero term of the map π+ X (ρ) − (π− X )−1(ρ) = π+ X (ρ) − π+ R(X)(ρ), where R(X) = (−P(x,−y),Q (x,−y)) for X = (P,Q ). See Fig. 1(b). From the above result we can study the general expression for r+ in the region Σ+ 0 and then, using the reversibility property, we can get r− from r+. In other words, we only need to consider θ ∈ [0, π ). Now, we recall the decomposition of a one-form given in [24] and the extension to piecewise differential forms, in polar coordi- nates, in [23] to the computation of the Poincaré map. The initial value problem (x′, y′) = (−y + εP(x, y, ε), x + εQ (x, y, ε)), (x(0, ε), y(0, ε)) = (ρ, 0), where P and Q are smooth functions, can be written, using usual polar coordinates (x, y) = (r cos θ, r sin θ ), as⎧⎪⎨⎪⎩dH + N∑ i=1 εiωi = 0, r(0, ε, ρ) = ρ. (7) Here H(r) = r2/2 and ωi = ωi(r, θ ) are smooth one-forms 2π-periodic in θ. The solution r(θ, ε, ρ) of (7)writes as r(θ, ε, ρ) =∑N i=0ri(θ, ρ)ε i and it can be easily checked that r(θ, 0, ρ) = r0(θ, ρ) ≡ ρ. See Fig. 2. The next result follows straightforward from the decomposi- tions of [24–26]. C.A. Buzzi et al. / Physica D 371 (2018) 28–47 31 Fig. 2. Solution of Eq. (7) for ε small enough. Lemma 2.2. Let Ω = α(r, θ )dr +β(r, θ )dθ be an arbitrary analytic one-form, 2π-periodic in θ and H(r) = r2/2. Then there exist functions h(r, θ ), S(r, θ ) and F (r) also 2π-periodic in θ and defined by F (r) = 1 2π ∫ 2π 0 β(r, ψ)dψ, S(r, θ ) = ∫ θ 0 β(r, ψ) dψ − F (r)θ and h(r, θ ) = ( α(r, θ ) − ∂S(r,θ ) ∂r ) /H ′(r), such that Ω = Ω0 + Ω1 where Ω0 = hdH + dS, Ω1 = F (r)dθ, and ∫ H=ρ2/2Ω 0 = 0,∫ H=ρ2/2Ω 1 = ∫ H=ρ2/2Ω. Next proposition allows us to obtain ri(θ, ρ) for every i from the decomposition of a one-form given in lemma behind, see [23]. Although it is a direct consequence of the results in [23], we have included it here by completeness and because this expression for the function ri is more explicit and clear than that of the aforemen- tioned reference. Proposition 2.3. Let r(θ, ε, ρ) = ∑N i=0ri(θ, ρ)ε i be the solution of (7). Assume that, from the decomposition given in Lemma 2.2, the functions Fi(r), hi(r, θ ), and Si(r, θ ) are defined inductively as h0 = 1 and −Ωi := − ∑i j=1ωjhi−j = hidH + dSi + Fidθ, for i = 1, . . . ,N. If the functions r0(θ, ρ) = ρ, and ri(θ, ρ), for i = 1, . . . ,N − 1 are known, then ρ rN (θ, ρ) = − 1 2 N−1∑ i=1 ri(θ, ρ)rN−i(θ, ρ) + G(S)(ρ, θ ) + ∫ θ 0 G(F )(ρ,ψ)dψ, (8) where S = (S1, . . . , SN ) and F = (F1, . . . , FN ). The operator G on f = (f1, . . . , fN ) is defined by G(f )(ρ, θ ) = fN (ρ, θ ) + N−1∑ i=1 N−i∑ j=1 Dj 1(fi)(ρ, θ ) × ∑ m1+2m2+···+ℓmℓ=N−i m1+m2+···+mℓ=j rm1 1 (θ, ρ)rm2 2 (θ, ρ) · · · rmℓℓ (θ, ρ) with Dj 1(g)(ρ, θ ) = ∂ jg(r,θ ) ∂r j ⏐⏐⏐ r=ρ . The above result applied to the computation of the first terms in the Taylor series in ε of r(θ, ρ, ε) gives ρr1(θ, ρ) = S1(ρ, θ ) + F1(ρ)θ, ρr2(θ, ρ) = − 1 2 r1(θ, ρ)2 + S2(ρ, θ ) + D1(S1)(ρ, θ )r1(θ, ρ) + F2(ρ)θ + F ′ 1(ρ) ∫ θ 0 r1(ψ, ρ)dψ, with D1(S1)(ρ, θ ) = ∂S1(r,θ ) ∂r ⏐⏐⏐ r=ρ . Proof of Proposition 2.3. From Theorem 2.4 in [23] and all the notation of the statement we have that for any N ∈ N, r(θ, ε, ρ) satisfies the following implicit equation r2(θ, ε, ρ) − ρ2 2 = N∑ i=1 εi (∫ θ 0 Fi(r(ψ, ε, ρ))dψ + Si(r(ψ, ε, ρ), ψ)|ψ=θ ψ=0 ) + OεN+1 . (9) First we substitute r(θ, ε, ρ) = ∑N i=0ri(θ, ρ)ε i in the above equa- tion. The proof follows, by induction, equating the coefficient of εn in both sides of Eq. (9). □ Remark 2.4. The functions h and S given in Lemma 2.2 are not unique. In fact, for a fixed pair h, S such thatΩ = hdH+dS+F (r)dθ we have Ω = h̃dH + d̃S + F (r)dθ with h̃ = h − f ′(r)/r and S̃ = S+ f (r) for any function f (r). The new decomposition satisfies also that h̃, S̃ are 2π-periodic in θ . In particular, the conclusions of the above proposition does not depend on which decomposition is used. Using the notions and results introduceduntil now,we can state the main result used in the next sections, the so called Poincaré– Pontryagin–Melnikov functions. The solution of system (6) in power series in ε writes as r±(θ, ε, ρ) = ∑N i=0r ± i (θ, ρ)εi, and it can be easily checked that r±(θ, 0, ρ) = r± 0 (θ, ρ) ≡ ρ. Lemma 2.1 and Proposition 2.3 allow us to compute r± i (θ, ρ) for every value of θ . Consequently the first non-vanishing term of the complete return map far from the sliding segment is obtained from the next result. Proposition 2.5. Denoting by M̂N (ρ) = r+ N (π, ρ) − r− N (−π, ρ), the Poincaré–Pontryagin–Melnikov functions, or Melnikov functions for shortness, of order N, for system (5), are given by M1(ρ) = M̂1(ρ) and MN (ρ) = M̂N (ρ) ⏐⏐ {Mk(ρ)≡0, k=1,...,N−1}, for N ≥ 2. Moreover, for ε small enough, each simple zero, ρ̂, of MN (ρ) gets a limit cycle of system (5). The limit cycle bifurcates from the level curve x2 + y2 = ρ̂2. We remark that the above proposition provides the functions M̂k(ρ), in terms of the coefficients of the perturbations, for every order k, but only coincides with the Melnikov function of a fixed order N when the previous ones vanish identically. This section ends with two technical results that they will be useful for the explicit computations in the next sections. Lemma 2.6. The next properties hold. (a) sk := 1 2π ∫ 2π 0 sinkψdψ = (k−1)!! k!! δ̂k where δ̂k = { 1, if k is even, 0, if k is odd. (b) ck := 1 2π ∫ 2π 0 coskψ dψ = sk. (c) Sk(θ ) := ∫ θ 0 sinkψ dψ − skθ has the form pk(cos θ ) if k odd or pk−1(cos θ ) sin θ if k even, where pj is a polynomial of degree j in cos θ. In particular, Sk(θ ) is 2π-periodic. (d) Ck(θ ) := ∫ θ 0 coskψ dψ − ckθ is a 2π-periodic function and it writes as pk−1(cos θ ) sin θ with pk−1 an even (odd) polynomial of degree k − 1 when k odd (even). (e) sk,ℓ := 1 2π ∫ 2π 0 sinkψ Sℓ(ψ) dψ = (k−1)!! k!! (ℓ−1)!! ℓ!! δ̂kδℓ where δℓ = { 0, if ℓ is even, 1, if ℓ is odd. (f) ck,ℓ := 1 2π ∫ 2π 0 coskψ Cℓ(ψ) dψ = 0. (g) Sk,ℓ(θ ) := ∫ θ 0 sinkψ Sℓ(ψ) dψ−sk,ℓθ is a 2π-periodic function. (h) Ck,ℓ(θ ) := ∫ θ 0 coskψ Cℓ(ψ) dψ is a 2π-periodic function and it writes as pk+ℓ(cos θ ) with pk+ℓ a polynomial of degree k + ℓ. (i) sck,ℓ := 1 2π ∫ 2π 0 sinkθ cos θ Sℓ(θ ) dθ = 1 k+1 (sk+1sℓ − sk+ℓ+1). 32 C.A. Buzzi et al. / Physica D 371 (2018) 28–47 (j) ∫ θ 0 sinkψ cosψ Sℓ(ψ) dψ − sk,ℓθ = 1 k+1 (Sℓ(θ )sin k+1θ − Sk+ℓ+1(θ ) + sℓSk+1(θ )) is a 2π-periodic function. Proof. The properties (a), (b), (e), and (f) follow directly using iteratively the integration by parts method. Alternatively, the first two can be done also integrating by quarters, see for example [27]. The periodicity properties (c), (d), (g), and (h) follow easily be- cause the correction by a multiple of θ converts the primitive of a periodic function into a periodic function. The second part of (d) follows also using iteratively the integration by parts method. Properties (i)–(j) follow by the integration by parts method and properties (a)–(b). □ Lemma 2.7. The next properties hold. (a) ∫ 2π 0 sinkθ Sk(θ ) dθ = 0. (b) ∫ ±π 0 Sk(ψ) dψ = ±π (k−1)!! k!! δk. (c) ŝk := Sk(π ) = Sk(−π ) = 2(k−1)!! k!! δk. (d) ŝk,k := Sk,k(π ) = Sk,k(−π ) = 1 2 (Sk(π ))2 = 2 ( (k−1)!! k!! )2 δk. (e) γ± k := ∫ ±π 0 sinkθdθ = ŝk ± πsk = (±1)k−1 (k−1)!! k!! σk, where σk = { π, if k is even, 2, if k is odd. (f) Ck,ℓ(π ) = ckC0,ℓ(π ) + cℓC0,k(π ) − Cℓ,k(π ), for k + ℓ odd. (g) Ck,ℓ(π ) = Ck,ℓ(−π ) and Ck,ℓ(π ) = 1−(−1)k+ℓ (k+ℓ)ℓ + ℓ−1 ℓ Ck,ℓ−2(π ). Moreover, Ck,ℓ(π ) = 0 for k+ ℓ even and Ck,ℓ(π ) ̸= 0 for k+ ℓ odd. Proof. Property (a) follows directly using Lemma 2.6(e) because δkδ̂k = 0 for every k. The other properties follow, as Lemma 2.6, using iteratively, if necessary, the integration by parts rule. □ 3. Polynomial perturbation This section is devoted to prove Theorem 1.1. The proof follows directly from the next three results. In the firstwe get the existence of the upper bounds for perturbations of every orderN. The second and the third give why they are reached for N = 1 and N = 2, respectively. Proposition 3.1. The maximum number of zeros of the Melnikov function MN associated to family (2) is n and Nn − 1 for N = 1 and N ≥ 2, respectively. Proposition 3.2. Consider the system X+ : (x′, y′) = ( − y + ε [ n−1 2 ]∑ j=0 a2j+1x2j+1, x + ε [ n 2 ]∑ j=0 a2jx2j ) , if y > 0, X− : (x′, y′) = (−y, x), if y < 0. (10) The function M1(ρ) is a complete polynomial of degree n with arbi- trary coefficients. In particular, there exist values of ak, k = 1, . . . , n, such that it has n simple zeros. Proposition 3.3. Consider the system X+ :(x′, y′) = ( −y − εxn + ε2 n∑ j=0 ajxj, x + ε n∑ j=2 bjxj ) , if x > 0, X− :(x′, y′) = ( −y − (−1)nεxn, x ) , if x < 0. The function M1(ρ) vanishes identically and the function M2(ρ) is a complete polynomial of degree 2n − 1. In particular, there exist coefficients aj, j = 0, . . . , n and bj, j = 2, . . . , n, such that the function M2(ρ) has exactly 2n − 1 simple zeros. Proof of Proposition 3.1. Following the procedure described in Section 2, we will prove first how are the expressions, as using polynomials in r , of the functions Fi, Si and hi, and second, how are the half return maps corresponding to the vector fields of both sides of the separation line. In polar coordinates, x = r cos θ and y = r sin θ , system (2) becomes⎧⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎩ rr ′ = N∑ i=1 εi(r cos θ f ± i (r cos θ, r sin θ ) + r sin θg± i (r cos θ, r sin θ )), r2θ ′ = 1 + N∑ i=1 εi(r cos θg± i (r cos θ, r sin θ ) − r sin θ f ± i (r cos θ, r sin θ )). (11) This system of equations writes, using differential one-forms, as ω± = rθ ′dr − rr ′dθ = rdr + N∑ i=1 εiω± i = 0. To simplify the notation, we use pn(r) to mean a polynomial of degree n in the variable r. We do not take into account the de- pendence in the other variables or parameters. Then, with this notation, we can write f ± i (r cos θ, r sin θ ) = pn(r), g± i (r cos θ, r sin θ ) = pn(r), and ω± i = pn(r)dr + rpn(r)dθ. From Lemma 2.2, and considering −ω± 1 = α± 1 (r, θ )dr + β± 1 (r, θ )dθ,we have F± 1 (r) = 1 2π ∫ 2π 0 β± 1 (r, ψ)dψ = rpn(r), S± 1 (r, θ ) = ∫ θ 0 β± 1 (r, ψ)dψ − F± 1 (r)θ = rpn(r), h± 1 (r, θ ) = 1 r ( α± 1 (r, θ ) − ∂S± 1 ∂r (r, θ ) ) = pn−1(r). Here we observe that, according to Remark 2.4, the functions F± k , S± k , and h± k given in Lemma 2.2 are not unique and can be chosen in such a way that they are polynomials. Now we prove by induction inN that the functions F± k , S± k and h± k associated to Eq. (11), of order N , can be written as F± k (r) = rpk(n−1)+1(r), S± k (r, θ ) = rpk(n−1)+1(r), h± k (r, θ ) = pk(n−1)(r), (12) for k = 1, . . . ,N. If we consider a perturbation of order N + 1 in (11), then the functions F± k , S± k , and h± k for k = 1, . . . ,N are, clearly, the same functions as the ones obtained for a perturbation of order N in (11). Consequently, we only need to show how they are for k = N + 1. In this case, −Ω± N+1 = α± N+1dr + β± N+1dθ = −(ω± N+1 + h± 1 ωN + · · · + h± Nω1) = (pn(r)dr + rpn(r)dθ )(1 + pn−1(r) + p2(n−1)(r) + · · · + pN(n−1)(r)) = pN(n−1)+n(r)dr + rpN(n−1)+n(r)dθ = p(N+1)(n−1)+1(r)dr + rp(N+1)(n−1)+1(r)dθ. C.A. Buzzi et al. / Physica D 371 (2018) 28–47 33 So, F± N+1(r) = 1 2π ∫ 2π 0 β± N+1(r, ψ)dψ = rp(N+1)(n−1)+1(r), S± N+1(r, θ ) = ∫ θ 0 β± N+1(r, ψ)dψ − F± N+1(r)θ = rp(N+1)(n−1)+1(r), h± N+1(r, θ ) = 1 r ( α± N+1(r, θ ) − ∂S± N+1 ∂r (r, θ ) ) = p(N+1)(n−1)(r). Now we can continue proving which are the expressions of the coefficients in ε of the half returnmap, that is, r± k , for k = 1, . . . ,N . For N = 1, the proof follows using Proposition 2.3, as we have r± 1 (θ, ρ) = 1 ρ (S± 1 (ρ, θ ) − F± 1 (ρ)θ ) = 1 ρ ρ pn(ρ) = pn(ρ). For N = 2 we write ρr± 2 (θ, ρ) = − 1 2 (r ± 1 (θ, ρ))2 + S± 2 (ρ, θ ) + D1S± 1 (ρ, θ )r± 1 (θ, ρ) + F± 2 (ρ)θ + (F± 1 )′(ρ) ∫ θ 0 r± 1 (ψ, ρ)dψ, with D1S± 1 (ρ, θ ) = ∂ ∂r S ± 1 (r, θ ) ⏐⏐ r=ρ . It can be seen easily, from (12), that r± 2 (θ, ρ) = p2n−1(ρ). Next we will prove by induction in N that r± N (θ, ρ) = pNn−1(ρ) ρN−2 . (13) We assume that (13) holds for k = 2, . . . ,N . According to Proposi- tion 2.3 we have that ρr± N+1(θ, ρ) is given as a sum of three terms. The first one is − 1 2 N∑ i=1 r± i (θ, ρ)r± N+1−i(θ, ρ) = N∑ i=1 pin−1(ρ) ρ i−2 p(N+1−i)n−1(ρ) ρ(N+1−i)−2 = p(N+1)n−2(ρ) ρN−3 = p(N+1)n−1(ρ) ρN−2 . The second term, using thenotationDj 1S ± i (ρ, θ )= ∂ jS± i ∂r j (r, θ ) ⏐⏐⏐ r=ρ , is given by S± N+1(ρ, θ ) + N∑ i=1 N+1−i∑ j=1 Dj 1S ± i (ρ, θ ) × ∑ m1+2m2+···+ℓmℓ=N+1−i m1+m2+···+mℓ=j (r± 1 (θ, ρ))m1 (r± 2 (θ, ρ))m2 · · · (r± ℓ (θ, ρ))mℓ = ρp(N+1)(n−1)+1(ρ) + pi(n−1)+2−j(ρ) × pm1n+m2(2n−1)+···+mℓ(ℓn−1)(ρ) rm3+2m4+3m5+···+(ℓ−2)mℓ = ρρN−2p(N+1)(n−1)+1(ρ) rN−2 + pi(n−1)+2−j+n(m1+2m2+···+ℓmℓ)−(m2+···+mℓ)(ρ) r (3m3+4m4+···+ℓmℓ)−2(m3+m4+···+mℓ) = p(N+1)n−1(ρ) ρN−2 + pi(n−1)+2−j+n(N+1−i)−(j−m1)(ρ) ρ(N+1−i−m1−2m2)−2(j−m1−m2) = p(N+1)n−1(ρ) ρN−2 + p(N+1)n+2+m1−i−2j(ρ) ρN+1+m1−i−2j = p(N+1)n−1(ρ) ρN−2 + ρ i+2j−m1−3 ρ i+2j−m1−3 p(N+1)n+2+m1−i−2j(ρ) ρN+1+m1−i−2j = p(N+1)n−1(ρ) ρN−2 . The third term has the same expression as above changing S by F and integrating in the variable θ, but there are no changes in the degrees as a function of r. So, as the three terms write in the same form, after division by ρ,we have r± N+1(θ, ρ) = p(N+1)n−1(ρ) ρ(N+1)−2 . Finally, from (13), the difference r+ N (π, ρ) − r− N (−π, ρ) also writes as pNn−1(ρ)/ρN−2. Then, fromProposition 2.5, themaximum number of positive zeros ofMN (ρ) is Nn − 1. □ Proof of Proposition 3.2. In polar coordinates, x = r cos θ , y = r sin θ, the half return map of the vector field X− is the identity. In particular r− 1 (−π, ρ) = 0 and M1(ρ) = r+ 1 (π, ρ). For simplify the reading we will omit in this proof the dependence of + because we only need to compute the half return map for θ ∈ (0, π ). Then, in polar coordinates, the vector field X+ of (10) writes⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ r ′ = ε ( no∑ j=0 a2j+1r2j+1cos2j+2θ + ne∑ j=0 a2jr2jcos2jθ sin θ ) , θ ′ = 1 + ε ( − no∑ j=0 a2j+1r2jcos2j+1θ sin θ + ne∑ j=0 a2jr2j−1cos2j+1θ ) , where no = [(n − 1)/2] and ne = [n/2]. In terms of differential one-forms we write ω = rθ ′dr − rr ′dθ = rdr + εω1, where − ω1 = α1(r, θ )dr + β1(r, θ )dθ, with α1(r, θ ) = no∑ j=0 a2j+1r2j+1cos2j+1θ sin θ − ne∑ j=0 a2jr2jcos2j+1θ, β1(r, θ ) = no∑ j=0 a2j+1r2j+2cos2j+2θ + ne∑ j=0 a2jr2j+1cos2jθ sin θ. According to Lemmas 2.2 and 2.6, and using ∫ 2π 0 cos2j+1ψ sinψ dψ = 0,we have that F1(r) = 1 2π ∫ 2π 0 β1(r, ψ)dψ = no∑ j=0 a2j+1c2j+2r2j+2, S1(r, θ ) = ∫ θ 0 β1(r, ψ)dψ − F1(r)θ = no∑ j=0 a2j+1r2j+2C2j+2(θ ) − ne∑ j=0 a2j 2j + 1 r2j+1(cos2j+1θ − 1). Clearly S1 is 2π-periodic in θ, because the functions C2j+2(θ ) also are. From Proposition 2.3 we have that r+ 1 (θ, ρ) = (S1(ρ, θ ) + F1(ρ)θ )/ρ and, from Lemma 2.7, r+ 1 (π, ρ) = ne∑ j=0 2a2j 2j + 1 ρ2j + no∑ j=0 a2j+1c2j+2ρ 2j+1, with c2j+2 ̸= 0. Consequently, there exist perturbation values such that M1(ρ) = r+ 1 (π, ρ) has n simple zeros. Because it is a polynomial of degree n = max{2no + 1, 2ne} with arbitrary coefficients. □ 34 C.A. Buzzi et al. / Physica D 371 (2018) 28–47 Proof of Proposition 3.3. First of all we consider the unified system X± : (x′, y′) = ( −y + δ±εxn + ε2 n∑ j=0 a± j x j, x + ε n∑ j=2 b± j x j, ) where δ+ = −1, a+ j = aj, b+ j = bj, δ− = −(−1)n, a− j = 0, and b− j = 0. In polar coordinates, x = r cosϕ, y = r sinϕ,we get⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ r ′ = ε ( δ±rncosn+1ϕ + n∑ j=2 b± j r jcosjϕ sinϕ ) + ε2 n∑ j=0 a± j r jcosj+1ϕ, ϕ′ = 1 − ε ( δ±rn−1cosnϕ sinϕ − n∑ j=2 b± j r j−1cosj+1ϕ ) − ε2 n∑ j=0 a± j r j−1cosjϕ sinϕ. (14) First we make a rotation of angle π/2 in order to have the discon- tinuity line as the horizontal axis. With the change θ = ϕ + π/2, cosϕ = sin θ and sinϕ = − cos θ, system (14) becomes⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ r ′ = ε ( δ±rnsinn+1θ − n∑ j=2 b± j r jsinjθ cos θ ) + ε2 n∑ j=0 a± j r jsinj+1θ, θ ′ = 1 + ε ( δ±rn−1sinnθ cos θ + n∑ j=2 b± j r j−1sinj+1θ ) + ε2 n∑ j=0 a± j r j−1sinjθ cos θ. or, using differential one-forms, we have ω± = rθ ′dr − rr ′dθ = rdr + εω± 1 + ε2ω± 2 = 0, where − ω± 1 = α± 1 (r, θ )dr + β± 1 (r, θ )dθ, with α± 1 (r, θ ) = −δ±rnsinnθ cos θ − n∑ j=2 b± j r jsinj+1θ, β± 1 (r, θ ) = δ±rn+1sinn+1θ − n∑ j=2 b± j r j+1sinjθ cos θ. The explicit expression of ω± 2 will appear later when we study the perturbation using the second order terms. For the first order study, according to Lemmas 2.2 and 2.6, we have F± 1 (r) = 1 2π ∫ 2π 0 β± 1 (r, ψ)dψ = δ±rn+1 2π ∫ 2π 0 sinn+1ψdψ = δ±rn+1sn+1, (15) S± 1 (r, θ ) = ∫ θ 0 β± 1 (r, ψ)dψ − F± 1 (r)θ = δ±rn+1Sn+1(θ ) − n∑ j=2 b± j 1 j + 1 r j+1sinj+1θ. (16) Clearly S± 1 is 2π-periodic because Sn+1(θ ) also is. From Proposi- tion 2.3 we have that r± 1 (θ, ρ) = (S± 1 (ρ, θ )+ F± 1 (ρ)θ )/ρ. It is clear, from Lemma 2.7, that S+ 1 (ρ, π )− S− 1 (ρ,−π ) = (δ+ − δ−)ŝn+1ρ n+1 and F+ 1 (ρ)π − F− 1 (ρ)(−π ) = πρn+1(δ+ + δ−)sn+1. We observe that, when n is even, δ+ = δ− and sn+1 = 0, and when n is odd, δ+ = −δ− and ŝn+1 = 0. So, the first Melnikov function is M1(ρ) = r+ 1 (π, ρ) − r− 1 (−π, ρ) = 0. The next step is the computation of the second Melnikov func- tionM2(ρ). In order to get it, first we compute the function h± 1 (r, θ ) using Lemma 2.2, h± 1 (r, θ ) = 1 r ( α± 1 (r, θ ) − ∂S± 1 ∂r (r, θ ) ) = −δ±rn−1(sinnθ cos θ + (n + 1)Sn+1(θ ) ) . The study of the second order terms, also using the decomposition of Lemma 2.2, uses − ω± 1 h ± 1 − ω± 2 = α± 2 (r, θ )dr + β± 2 (r, θ )dθ, with β± 2 (r, θ ) = n∑ j=0 a± j r j+1sinj+1θ + ( δ±rn+1sinn+1θ − n∑ j=2 b± j r j+1sinjθ cos θ ) × ( −δ±rn−1 ( sinnθ cos θ + (n + 1)Sn+1(θ ) )) = n∑ j=0 a± j r j+1sinj+1θ − (δ±)2r2n ( sin2n+1θ cos θ + (n + 1)sinn+1θSn+1(θ ) ) + n∑ j=2 δ±b± j r n+jsinn+jθcos2θ + (n + 1) × n∑ j=2 δ±b± j r n+jsinjθ cos θSn+1(θ ). We have not written the expression of α± because they are not needed for the computation of r± 2 (θ, ρ). Using again Lemma 2.2 and the properties of Lemma 2.6, and∫ 2π 0 sin2n+1θ cos θdθ = 0 we have that F± 2 (r) = 1 2π ∫ 2π 0 β± 2 (r, ψ)dψ = n∑ j=0 a± j r j+1sj+1 + n∑ j=2 δ±b± j r n+jsn+j − n∑ j=2 δ±b± j r n+jsn+j+2 + (n + 1) n∑ j=2 δ±b± j r n+j × 1 j + 1 (sn+1sj+1 − sn+j+2). So, F± 2 (r) = n∑ j=0 a± j r j+1sj+1 + n∑ j=2 δ±b± j r n+j × ( sn+j − sn+j+2 + n + 1 j + 1 (sn+1sj+1 − sn+j+2) ) . C.A. Buzzi et al. / Physica D 371 (2018) 28–47 35 Again from Lemmas 2.2 and 2.6, we have S± 2 (r, θ ) = ∫ θ 0 β± 2 (r, ψ)dψ − F± 2 (r)θ = n∑ j=0 a± j r j+1Sj+1(θ ) − (δ±)2r2n × ( 1 2n + 2 sin2n+2θ + (n + 1)Sn+1,n+1(θ ) ) + n∑ j=2 δ±b± j r n+j ( Sn+j(θ ) − Sn+j+2(θ ) ) + (n + 1) n∑ j=2 δ±b± j r n+j 1 j + 1 × ( Sn+1(θ )sinj+1θ − Sn+j+2(θ ) + sn+1Sj+1(θ ) ) . Finally, from Propositions 2.3 and 2.5, the second Melnikov function is M2(ρ) = r+ 2 (π, ρ) − r− 2 (−π, ρ). We compute this difference term by term. The first and third terms vanish because r+ 1 (π, ρ) − r− 1 (−π, ρ) = 0 and, from (15) and (16), ∂S+ 1 (r, π ) ∂r ⏐⏐⏐⏐ r=ρ − ∂S− 1 (r,−π ) ∂r ⏐⏐⏐⏐ r=ρ = (n + 1)(δ+ − δ−)ŝn+1 = 0. The above expression vanishes using the same argument explained in the paragraph before (16). From Lemma 2.7 we have 1 ρ (F+ 1 )′(ρ) ∫ π 0 r+ 1 (ψ, ρ)dψ = (δ+)2(n + 1)ρ2n−1sn+1 × (∫ π 0 Sn+1(ψ)dψ − sn+1 π2 2 ) + n∑ j=2 n + 1 j + 1 bjρn+j−1sn+1 × (ŝj+1 + πsj+1) = − (n + 1)π2 2 s2n+1ρ 2n−1 + n∑ j=2 n + 1 j + 1 bjsn+1γ + j+1ρ n+j−1. (17) In the last equality, the first term vanishes because, from the expressions of Lemma2.6, sn+1 is amultiple of δn+1 and the integral is a multiple of δ̂n+1. Now, doing similar simplifications, we get 1 ρ (F− 1 )′(ρ) ∫ −π 0 r− 1 (ψ, ρ)dψ = (δ−)2(n + 1)ρ2n−1sn+1 × (∫ −π 0 Sn+1(ψ)dψ − sn+1 π2 2 ) = − (n + 1)π2 2 s2n+1ρ 2n−1. (18) The contribution of the last term, in the expressions of r± 2 , toM2 is obtained from the difference of (17) and (18). Hence, the first part of the second Melnikov function writes as n∑ j=2 n + 1 j + 1 bjsn+1γ + j+1ρ n+j−1. (19) Now, we continue computing the second and fourth terms. Adding the following two expressions 1 ρ S+ 2 (ρ, π ) = n∑ j=2 bjρn+j−1 ( n + j + 2 j + 1 ŝn+j+2 − ŝn+j − n + 1 j + 1 sn+1ŝj+1 ) + n∑ j=0 ajρ jŝj+1 − (n + 1)ρ2n−1ŝn+1,n+1, π ρ F+ 2 (ρ) = n∑ j=2 πbjρn+j−1 ( n + j + 2 j + 1 sn+j+2 − sn+j − n + 1 j + 1 sn+1sj+1 ) + n∑ j=0 πajρ jsj+1, we obtain that the contribution of that terms, corresponding to the vector field X+, is n∑ j=0 ajρ jγj+1 + n∑ j=2 bjρn+j−1 ( n + j + 2 j + 1 γ+ n+j+2 − γ+ n+j − n + 1 j + 1 sn+1γ + j+1 ) − (n + 1)ρ2n−1ŝn+1,n+1. (20) The contribution corresponding to the vector field X− is the sum of the following two terms: 1 ρ S− 2 (ρ,−π ) = −(n + 1)ρ2n−1ŝn+1,n+1, − π ρ F− 2 (ρ) = 0, that is, after adding them, only − (n + 1)ρ2n−1ŝn+1,n+1. (21) Consequently the difference of Eqs. (20) and (21) provides the second part of the second Melnikov function: n∑ j=0 ajρ jγ+ j+1 + n∑ j=2 bjρn+j−1 ( n + j + 2 j + 1 γ+ n+j+2 − γ+ n+j − n + 1 j + 1 sn+1γ + j+1 ) . (22) Finally, adding (19) and (22) and using that σn+j+2 = σn+j and n + j + 1 j + 1 γ+ n+j − γ+ n+j = n j + 1 γ+ n+j, we obtain the expression for the second Melnikov function M2(ρ) = n∑ j=0 ajγ+ j+1ρ j + n∑ j=2 n j + 1 bjγ+ n+jρ n+j−1, with γ+ k ̸= 0, for every k, see also Lemma 2.7. The proof finishes because M2 is a complete polynomial of degree 2n − 1 with arbitrary coefficients. So, we can choose them such thatM2(ρ) has exactly 2n − 1 simple roots. □ 4. Perturbations of order one and two for someLiénard families This section is devoted to prove Theorem 1.2. It is a direct consequence of the next two results. Proposition 4.1. Consider the piecewise Liénard system (x′, y′) = ( −y + εf ± n (x), x ) , (23) 36 C.A. Buzzi et al. / Physica D 371 (2018) 28–47 defined inΣ± 0 = {(x, y) ∈ R2 ; ±y > 0}, where f ± n (x) = ∑n j=0a ± j x j. The first Melnikov function writes as M1,n(ρ) = ρ pm(ρ2) for a given polynomial pm of degree m = [(n − 1)/2]. Then, it has at most m positive zeros. Moreover, there exist values of the parameters a± j such that M1,n has exactly m positive simple zeros. Additionally, M1,n ≡ 0 if and only if a+ 2k−1 = −a− 2k−1 for k = 1, . . . ,m + 1. Proposition 4.2. Consider the piecewise Liénard system (x′, y′) = ( −y + εf ± n (x) + ε2g± n (x), x ) , (24) defined in Σ± 0 = {(x, y) ∈ R2 ; ±y > 0}, where f ± n (x) = ∑n j=0a ± j x j and g± n (x) = ∑n j=0b ± j x j, satisfying a+ 2k−1 = −a− 2k−1 for k = 1, . . . ,m+1, wherem = [(n−1)/2]. Then, the firstMelnikov function vanishes identically, M1,n ≡ 0, and the second Melnikov function writes as M2,n(ρ) = ρ pm(ρ2) + qn−1(ρ2) for given polynomials pm and qn−1 of degrees m = [(n − 1)/2] and n − 1, respectively. Then, it has at most m + n positive zeros. Moreover, there exist values of the parameters a± j and b± j such that M2,n has exactly m + n positive simple zeros. Proof of Proposition 4.1. We start writing system (23), with polar coordinates, as the piecewise one-form{ dH + εω+ 1,n = 0 if θ ∈ [0, π ), dH + εω− 1,n = 0 if θ ∈ [π, 2π ), where H = r2/2 and −ω± 1,n = α± 1,n(r, θ )dr + β± 1,n(r, θ )dθ,with α± 1,n(r, θ ) =f ± n (r cos θ ) sin θ = n∑ j=0 a± j r jcosjθ sin θ, β± 1,n(r, θ ) =f ± n (r cos θ )r cos θ = n∑ j=0 a± j r j+1cosj+1θ. (25) Using Lemma 2.2 and Proposition 2.3, from the above expressions, we get r± 1,n(θ ) := ρr± 1,n(θ ) = F± 1,n(ρ)θ + S± 1,n(ρ, θ ), (26) where F± 1,n(r) = 1 2π ∫ 2π 0 β± 1,n(r, θ )dθ and S± 1,n(r, θ ) = ∫ θ 0 β± 1,n(r, ψ)dψ − F± 1,n(r)θ. Now, we prove that the functions F± 1,n(r) and S± 1,n(r, θ ) are poly- nomials in r with constant coefficients and trigonometric func- tions, respectively. For n = 0, we have α± 1,0 = a± 0 sin θ and β± 1,0 = b± 0 r cos θ. Hence, using Proposition 2.5, M1,0(ρ) = 0, because F± 1,0(r) = 0 and S1,0(r, θ ) = a± 0 r sin θ . For any n, we can write F± 1,n+1(r) = [ n+1 2 ]∑ k=1 a± 2k−1c2kr 2k + a± n+1cn+2rn+2 = F± 1,n(r) + a± n+1cn+2rn+2, (27) because, from Lemma 2.6, we have cj = 0 for every j odd. We remark that the monomial in rn+2 only appears when n is even. Additionally, we can write S± 1,n+1(r, θ ) = n+1∑ j=0 a± j r j+1 ∫ θ 0 cosj+1ψdψ − F1,n+1(r)θ = ⎛⎝ n∑ j=0 a± j r j+1 ∫ θ 0 cosj+1ψdψ − F1,n(r)θ ⎞⎠ + ( a± n+1r n+2 ∫ θ 0 cosn+2ψdψ − a± n+1cn+2rn+2θ ) = S± 1,n(r, θ ) + a± n+1 (∫ θ 0 cosn+2ψdψ − cn+2θ ) rn+2 and, consequently, S± 1,n+1(r, θ ) = S± 1,n(r, θ ) + a± n+1 Cn+2(θ ) rn+2. (28) The above expression can also be written as S± 1,n(r, θ ) = r sin θ Gn(r, cos θ ), (29) where Gn(r, cos θ ) is a polynomial of degree n in r which coeffi- cients are polynomials in cos θ. This property follows by induction on n using also Lemma 2.6(d). Straightforward computations show, from (26), (27), and (28), that ρr± 1,n+1(θ, ρ) = ρr± 1,n(θ, ρ) + a± n+1(Cn+2(θ ) + cn+2θ )ρn+2, (30) with r± 1,0(θ ) = a± 0 sin θ, and, consequently, r± 1,n(θ, ρ) = n∑ j=0 a± j (Cj+1(θ ) + cj+1θ )ρ j. Therefore, using also Proposition 2.5, we can write the first Mel- nikov function as M1,n+1(ρ) = r+ 1,n+1(π, ρ) − r− 1,n+1(−π, ρ) = r+ 1,n(π, ρ) − r− 1,n(−π, ρ) + π (a+ n+1 + a− n+1)cn+2ρ n+1 = M1,n(ρ) + π (a+ n+1 + a− n+1)cn+2ρ n+1. Then, using Lemma 2.6, the function M1,n is a polynomial in ρ of degree 2m + 1, with m = [(n − 1)/2] , because M1,0(ρ) ≡ 0, cj = 0 for j odd, and cj ̸= 0 for j even. More concretely, it writes as M1,n(ρ) = ρ pm(ρ2) with pm(ρ) = π ( c2(a+ 1 + a− 1 ) + c4(a+ 3 + a− 3 )ρ + · · · + c2m+2(a+ 2m+1 + a− 2m+1)ρ m) , (31) where c2k ̸= 0 for k = 1, . . . ,m + 1. Clearly, from the above expression, all the statements follow. The function M1,n has at most m zeros and we can choose all the perturbation parameters in order to get these zeros as simple ones. Moreover,M1,n ≡ 0 if and only if a+ 2k−1 = −a− 2k−1 for k = 1, . . . ,m + 1. □ Proof of Proposition 4.2. First, we provide a short description of the proof, with the main properties. Then, we proceed to detail all the steps carefully. From Proposition 4.1 we know that the first Melnikov function vanishes identically when a+ 2k−1 + a− 2k−1 = 0. In the proof we will refer this property as (M1). The secondMelnikov functionwrites as a sum of two functions, M2,n(ρ) = M (1) 2,n(ρ) + M (2) 2,n(ρ), that they correspond to the perturbation parameters of first and second order, aj and bj, respectively in (24). We can also use C.A. Buzzi et al. / Physica D 371 (2018) 28–47 37 Proposition 4.1, changing ε to ε2, to see that M (2) 2,n(ρ) = ρ pm(ρ2) for some polynomial pm, of degreem = [(n−1)/2],with arbitrary coefficients depending linearly on bj. Hence, we fix our attention to studyM (1) 2,n(ρ).Wewill prove by induction thatM (1) 2,n(ρ) = qn−1(ρ2) for some polynomial qn−1, of degree n − 1, with coefficients de- pending quadratically on aj. This proves the first statements. That is, the polynomial M2,n(ρ) = ρpm(ρ2) + qn−1(ρ2) has m + n + 1 monomials and, consequently from the Descartes signs’ rule, the maximum number of positive zeros is m + n. The existence of such number of simple positive zeros is provedwithout the explicit computation of all their coefficients. Only themaximal degree term is necessary to be given. Because, at each step of the induction procedure, we add a complete polynomial in ρ2 that has a new monomial of higher degree and it is independent with respect to the smaller degree monomials. Therefore, the last statement follows. In fact, the coefficient of maximal degree, e2n−2, takes different expressions depending on the parity of the degree of the perturbation. More concretely, e2n−2 = { a+ n ( a+ n+1 + a− n+1 ) Kn, when n is even, a+ n+1 ( a+ n + a− n ) Kn, when n is odd, with positive constants Kn satisfying the recursive relation Kn = n − 1 n + 1 Kn−2 − 4 (4(n − 1)2 − 1)(n + 1) , (32) and K1 = 2, K2 = 2/9. The positiveness follows by induction proving that Kn > 1/(n + 1)2 > 0.We remark that the recurrence that defines Kn does not change with the parity of n, but the initial values does. We start writing system (24), with polar coordinates, in its equivalent piecewise one-form{ dH + εω+ 1,n + ε2ω+ 2,n = 0 if θ ∈ [0, π ), dH + εω− 1,n + ε2ω− 2,n = 0 if θ ∈ [π, 2π ), where H = r2/2, −ω± 1,n is defined in (25), and −ω± 2,n = α± 2,n(r, θ )dr + β± 2,n(r, θ )dθ,with α± 2,n(r, θ ) = f ± n (r cos θ ) sin θ = n∑ j=0 b± j r jcosjθ sin θ, β± 2,n(r, θ ) = f ± n (r cos θ )r cos θ = n∑ j=0 b± j r j+1cosj+1θ. The secondMelnikov function, from Proposition 2.5, is given by M2,n(ρ) = r+ 2,n(π, ρ) − r− 2,n(−π, ρ). (33) Furthermore, from Proposition 2.3, we get ρr± 2,n(θ, ρ) = − 1 2 r± 1,n(θ, ρ) 2 + S± 2,n(ρ, θ ) + D1(S± 1,n)(ρ, θ )r ± 1,n(θ, ρ) + F± 2,n(ρ)θ + (F± 1,n) ′(ρ) ∫ θ 0 r± 1,n(ψ, ρ)dψ, and D1(S± 1,n)(ρ, θ ) = ∂S± 1,n(r,θ ) ∂r ⏐⏐⏐⏐ r=ρ . Here the functions F± 1,n, S ± 1,n, and r± 1,n are defined in (27), (28), and (30), respectively, and the functions F± 2,n and S± 2,n are obtained applying Lemma 2.2 to the one-form −Ω± 2,n = −(ω± 2,n + h± 1,nω ± 1,n). (34) From the hypothesis (M1) and the proof of Proposition 4.1 we have r+ 1,n(π, ρ) = r− 1,n(−π, ρ). Moreover, from (29), S± 1,n(r,±π ) ≡ 0. Hence, Eq. (33) writes as M2,n(ρ) = 1 ρ (( (F+ 1,n) ′(ρ) ∫ π 0 r+ 1,n(ψ, ρ)dψ − (F− 1,n) ′(ρ) × ∫ −π 0 r− 1,n(ψ, ρ)dψ ) + (F+ 2,n(ρ) + F− 2,n(ρ))π + (S+ 2,n(ρ, π ) − S− 2,n(ρ,−π )) ) . (35) From (34) it is clear that each F± 2,n(ρ) and S± 2,n(ρ) can be written as the sum of two terms, the one that depends linearly on the coefficients of the second order (computed from ω± 2,n) and the ones that depend quadratically on the coefficients of the first order (computed from h± 1,nω ± 1,n). We denote them by F±,(2) 2,n+1(ρ), S±,(2) 2,n+1(ρ), F ±,(1) 2,n+1(ρ), and S±,(1) 2,n+1(ρ), respectively. With this notation (35) writes as M2,n(ρ) = M (2) 2,n(ρ) + M (1) 2,n(ρ), with M (2) 2,n(ρ) = (F+,(2) 2,n (ρ) + F−,(2) 2,n (ρ))π + (S+,(2) 2,n (ρ, π ) − S−,(2) 2,n (ρ,−π )) (36) and M (1) 2,n(ρ) = 1 ρ (( (F+ 1,n) ′(ρ) ∫ π 0 r+ 1,n(ψ, ρ)dψ − (F− 1,n) ′(ρ) ∫ −π 0 r− 1,n(ψ, ρ)dψ ) (F+,(1) 2,n (ρ) + F−,(1) 2,n (ρ))π + (S+,(1) 2,n (ρ, π ) − S−,(1) 2,n (ρ,−π )) ) . (37) Also fromProposition 4.1, (36)writes asM (2) 2,n(ρ) = ρ pm(ρ2) where pm is the polynomial (31) changing ai by bi.Werecall that its degree is m = [(n − 1)/2]. So, we continue the proof showing that (37) writes as M (1) 2,n(ρ) = qn−1(ρ2), for some polynomial qn−1 of degree n−1. As in the proof of Proposition 4.1, wewill prove, by induction in n, thatM (1) 2,n(ρ) function writes as M (1) 2,n+1(ρ) = M (1) 2,n(ρ) + Pn(ρ2), (38) where the polynomial Pn(ρ) has degree n. The proof will be done considering the three summands of Eq. (37) separately in the following parts (I), (II), and (III). (I) The contribution of the first term toM (1) 2,n+1 follows from the evaluation of the function (F± 1,n+1) ′(ρ) ∫ θ 0 r± 1,n+1(ψ, ρ)dψ at θ = ±π . Moreover, we will see that these expressions only provide even terms in ρ for the polynomial Pn(ρ). From Eqs. (27) and (30) we get (F± 1,n+1) ′(ρ) ∫ θ 0 r± 1,n+1(ψ, ρ)dψ = (F± 1,n+1) ′(ρ) ∫ θ 0 ( r± 1,n(ψ, ρ) + a± n+1 × (Cn+2(ψ) + cn+2ψ)ρn+1) dψ = (F± 1,n) ′(ρ) ∫ θ 0 r± 1,n(ψ, ρ)dψ + U± 1 (θ ) + U± 2 (θ ) + U± 3 (θ ), 38 C.A. Buzzi et al. / Physica D 371 (2018) 28–47 where F± 1,n(r) = ∑[ n+1 2 ] j=1 a± 2j−1c2jr 2j and U± 1 (θ ) = (F± 1,n) ′(ρ) ∫ θ 0 a± n+1(Cn+2(ψ) + cn+2ψ)ρn+1dψ, U± 2 (θ ) = a± n+1cn+2(n + 2) ∫ θ 0 r± 1,n(ψ, ρ)dψρ n+1, U± 3 (θ ) = a± n+1cn+2(n + 2) ∫ θ 0 a± n+1(Cn+2(ψ) + cn+2ψ)ρ2n+2dψ. Now, we will compute the evaluation of U± 1 (θ ), U± 2 (θ ), and U± 3 (θ ) at θ = ±π separately. First, from Lemma 2.6, we have U± 1 (θ ) = a± n+1 ( C0,n+2(θ ) + cn+2 2 θ2 ) [ n+1 2 ]∑ j=1 2ja± 2j−1c2jρ 2j+n and, due to hypothesis (M1) and C0,n+2(π ) = C0,n+2(−π ), U+ 1 (π ) − U− 1 (−π ) ρ = ( C0,n+2(π ) + cn+2 2 π2 ) ( a+ n+1 + a− n+1 ) × [ n+1 2 ]∑ j=1 2ja+ 2j−1c2jρ 2j+n−1. (39) This last expression is non identically zero when n+ 1 is even also due to hypothesis (M1) and it has only monomials of even degree. This is because, in this case, cn+2 = 0. So, we write the right hand side of (39) as ( a+ n+1 + a− n+1 ) [ n−1 2 ]∑ j=0 η (1) j a+ 2j+1ρ 2j+n+1, (40) for some real constants η(1)j , and moreover the term of higher degree in ρ is (n + 1)cn+1a+ n ( a+ n+1 + a− n+1 ) C0,n+2(π )ρ2n. (41) We remark that, in the above expression, cn+1C0,n+2(π ) ̸= 0. Second, from Lemma 2.6(h) and (30), we have U± 2 (θ ) = a± n+1cn+2(n + 2) n∑ j=0 a± j ( C0,j+1(θ ) + cj+1 2 θ2 ) ρn+j+1. So, as cn+2(a+ n+1 + a− n+1) = 0 and c2j+1 = 0, for all j and n, we get U+ 2 (π ) − U− 2 (−π ) ρ = cn+2(n + 2) [n/2]∑ j=0 (a+ n+1a + 2j − a− n+1a − 2j) × ( C0,2j+1(π ) + c2j+1 2 π2 ) ρ2j+n + cn+2(n + 2)(a+ n+1 + a− n+1) × [n/2]∑ j=0 a+ 2j+1 ( C0,2j+2(π ) + c2j+2 2 π2 ) ρ2j+1+n = cn+2(n + 2) [n/2]∑ j=0 (a+ n+1a + 2j − a− n+1a − 2j)C0,2j+1(π )ρ2j+n. (42) We have used Lemmas 2.7(g), 2.6(b) and (M1). The expression (42) is non identically zero only when n is even. Consequently, the last expression in (42) writes as a+ n+1 [n/2]∑ j=0 η (2) j (a+ 2j + a− 2j)ρ 2j+n, (43) for some real constants η(2)j . Clearly it has only monomials of even degree and the term of higher degree in ρ is (n + 2)cn+2a+ n+1(a + n + a− n )C0,n+1(π )ρ2n. (44) We also remark that, in the above expression, cn+2C0,n+1(π ) ̸= 0. Third, using also Lemma 2.6(h), we have U± 3 (θ ) = (a± n+1) 2cn+2(n + 2) ( C0,n+2(θ ) + cn+2 θ2 2 ) ρ2n+2 and, consequently, U+ 3 (π ) − U− 3 (−π ) ρ = (n + 2)cn+2 ( C0,n+2(π ) + cn+2 π2 2 ) × ( (a+ n+1) 2 − (a− n+1) 2 ) ρ2n+2 = 0. (45) Because, when n is even, by hypothesis (M1), a+ n+1 + a− n+1 = 0 and when n is odd cn+2 = 0. (II) The contribution of the second term toM (1) 2,n+1 in (35) follows from the computation of F±,(1) 2,n+1(ρ). In fact we will prove that both functions vanish identically. From the second summand of (34) we have F±,(1) 2,n+1(ρ) = 1 2π ∫ 2π 0 h± 1,n+1(ρ, θ )β ± 1,n+1(ρ, θ )dθ. (46) Straightforward computations show, from the definition of h1, see Lemma 2.2, and Eqs. (25) and (28), that h± 1,0(r, θ ) = 0 and h± 1,n+1(r, θ ) = h± 1,n(r, θ ) − a± n+1(n + 1)Cn(θ )rn. In particular, by induction, h± 1,n+1(r, θ ) = − n∑ j=0 a± j+1(j + 1)Cj(θ )r j. Finally, also from Lemma 2.6(d) and Eq. (25), it follows that (46) is zero. (III) The contribution of the third term toM (1) 2,n+1 in (35) follows from the evaluation of the functions S±,(1) 2,n+1(ρ, θ ) at θ = ±π . Moreover, we will see that these expressions only provide even terms in ρ for the polynomial Pn(ρ). From the second summand of (34) and the fact that F±,(1) 2,n+1(ρ) ≡ 0,we have S±,(1) 2,n+1(ρ, θ ) = ∫ θ 0 h± 1,n+1(ρ,ψ)β± 1,n+1(ρ,ψ)dψ = ∫ θ 0 (h± 1,n(ρ,ψ) − a± n+1(n + 1)Cn(ψ)ρn) × (β± 1,n(ρ,ψ) + a± n+1ρ n+2cosn+2ψ)dψ = S±,(1) 2,n (ρ,ψ) + V± 1 (θ ) + V± 2 (θ ) + V± 3 (θ ), where V± 1 (θ ) = −a± n+1(n + 1)ρn ∫ θ 0 β± 1,n(ρ,ψ)Cn(ψ)dψ, V± 2 (θ ) = a± n+1ρ n+2 ∫ θ 0 cosn+2ψh± 1,n(ρ,ψ)dψ, V± 3 (θ ) = −(a± n+1) 2(n + 1)ρ2n+2 ∫ θ 0 cosn+2ψCn(ψ)dψ. Now we will compute the evaluation of V± 1 (θ ), V± 2 (θ ), and V± 3 (θ ) at θ = ±π separately. C.A. Buzzi et al. / Physica D 371 (2018) 28–47 39 First, from Lemma 2.6(h) and the definition of β± 1,n in (25), we have that V± 1 (θ ) = −a± n+1(n + 1) n∑ j=0 a± j Cj+1,n(θ )ρn+j+1. So, from Lemma 2.7, V+ 1 (π ) − V− 1 (−π ) ρ = (n + 1) n∑ j=0 (a− n+1a − j − a+ n+1a + j ) × Cj+1,n(π )ρn+j. (47) The above polynomial has onlymonomials of even degree because, when n + j + 1 is even, the corresponding Cj+1,n(π ) = 0, see also Lemma 2.7(g). Moreover, the term of higher degree writes as (n + 1)(a− n+1a − n − a+ n+1a + n )Cn+1,n(π )ρ2n (48) and the right hand side of (47) writes as (40) or (43), changing the coefficients η(1)j and η(2)j , when n is odd or even, respectively. The maximal degree term is different from zero because Cn+1,n(π ) ̸= 0. Second, from Lemma 2.6(h), we have that V± 2 (θ ) = a± n+1ρ n+2 ∫ θ 0 cosn+2ψ ⎛⎝− n−1∑ j=0 a± j+1(j + 1)Cj(ψ)ρ j ⎞⎠ dψ = −a± n+1 n−1∑ j=0 a± j+1(j + 1)Cn+2,j(θ )ρn+j+2. So, again from Lemma 2.7, we have that V+ 2 (π ) − V− 2 (−π ) ρ = n−1∑ j=0 (a− n+1a − j+1 − a+ n+1a + j+1)(j + 1) × Cn+2,j(π )ρn+j+1. (49) The above polynomial has only monomials of even degree because when n + j + 2 is even the corresponding Cn+2,j(π ) = 0, see also Lemma 2.7(g). Moreover, the term of higher degree writes as n(a− n+1a − n − a+ n+1a + n )Cn+2,n−1(π )ρ2n (50) and, as before, the right hand side of (49) writes as (40) or (43), changing the coefficients η(1)j and η(2)j , when n is odd or even, respectively. Here, also the maximal degree term is different from zero because Cn+2,n−1(π ) ̸= 0. Third, from Lemma 2.6(h), we have that V± 3 (θ ) = −(a± n+1) 2(n + 1)ρ2n+2 ∫ θ 0 cosn+2ψCn(ψ)dψ = −(a± n+1) 2(n + 1)ρ2n+2Cn+2,n(θ ). So, using that Cn+2,n(π ) = Cn+2,n(−π ) = 0, V+ 3 (π ) − V− 3 (−π ) ρ = 0. (51) From all the computations done above we have proved, by induction, that (37) writes asM (1) 2,n(ρ) = qn−1(ρ2) for some polyno- mial qn−1 of degree n−1. This is due to the fact that the polynomial Pn(ρ2), in (38), follows adding the polynomials (39), (42), (45), (47), (49), and (51), all of them having only monomials of even degree. More concretely, it writes in the form (40) or (43) when n is odd or even, respectively. Moreover, the coefficient of maximal degree term, ρ2n, is obtained adding (41), (44), (48), and (50). In what follows we compute it explicitly because it is necessary to conclude the proof, showing also that it does not vanish. We remark that in both expressions, Pn(ρ2) and ρ2n, the n refers to the perturbation of degree n + 1, because of the induction procedure. Next, wewill write, to simplify the computations, the coefficient of maximal degree corresponding toM2,n+1(ρ) orM (1) 2,n+1(ρ), because they coincide. When n is odd, we have that the coefficient of maximal degree is a+ n ( a+ n+1 + a− n+1 ) K̃n with K̃n = (n + 1)cn+1C0,n+2(π ) − (n + 1)Cn+1,n(π ) − nCn+2,n−1(π ). This last expression can be simplified using properties (f) and (g) of Lemma 2.7. First, in the second and third summands, writing it as K̃n = 4(n + 1) (2n + 1)(n + 2) − n n + 2 Cn−1,n(π ), and then, in the second summand of the above expression, to write a recurrent expression for the coefficient of maximal degree, K̃n = n n + 2 K̃n−2 − 4 (4n2 − 1)(n + 2) . (52) The first value follows by direct computation and we obtain K̃1 = 2/9. All the coefficients are positive because we can prove, by induction, that K̃n > 1/(n+2)2 > 0.When n is even, the coefficient of maximal degree writes now as a+ n+1 ( a+ n + a− n ) K̃n with K̃n = (n + 2)cn+2C0,n+1(π ) − (n + 1)Cn+1,n(π ) − nCn+2,n−1(π ). Using similarly Lemma 2.7 as previously we have the same re- currence relation (52). But starting at K̃0 = 2. Clearly, K̃n is also positive because the same lower bound applies. We remark that the expression for (32) follows from (52) be- cause K̃n−1 = Kn. Consequently, Kn > 1/(n + 1)2 > 0. The last statement of the theorem follows writing the second Melnikov functionM2,n(ρ) asρpm(ρ2)+qn−1(ρ2), for a− 2k = 0when k = 0, . . . , [n/2] and b− 2k+1 = 0 when k = 0, . . .m. As we have mentioned before, the polynomial pm(ρ) writes as (31) changing a+ 2k+1 to b+ 2k+1 and a− 2k+1 to 0. The polynomial qn−1(ρ) is obtained by induction, using alternatively (40) and (43). So, we have M2,n(ρ) = m∑ k=0 ζkb+ 2k+1ρ 2k+1 + n−1∑ ℓ=0 ( ∑ i+j=2ℓ+1 i,j 0 small enough, because at this step a+ 3 = 0. The third bifurcates also from infinity adding b+ 3 < 0 small enough. At the fourth step, adding a+ 3 > 0 also small enough, the previous zeros remain, by continuity, close to the ones appeared in the previous step, but the coefficient of ρ2 changes slightly. The proof ends when all the coefficients of (53) have been added. □ From the last proof, it is clear that M2,n, see (53), has exactly m+n positive simple zeros and nomore. Moreover, only Kℓ and ζk should be different from zero. The next example shows the explicit 40 C.A. Buzzi et al. / Physica D 371 (2018) 28–47 expression forM2 when n = 6, M2,6(ρ) = 1774 24255 a+ 5 a + 6 ρ 10 + ( − 702 1225 a+ 3 a + 6 + 578 945 a+ 4 a + 5 ) ρ8 + ( − 1426 735 a+ 1 a + 6 + 26 21 a+ 2 a + 5 + 58 525 a+ 3 a + 4 ) ρ6 + 5π 16 b+ 5 ρ 5 + ( 2a+ 0 a + 5 − 74 75 a+ 1 a + 4 + 14 15 a+ 2 a + 3 ) ρ4 + 3π 8 b+ 3 ρ 3 + ( 2a+ 0 a + 3 + 2 9 a+ 1 a + 2 ) ρ2 + π 2 b+ 1 ρ + 2a+ 1 a + 0 . Then, it seems numerically that the eight positive simple zeros that the above polynomial can have, can be located arbitrarily in the positive real line. For example, if we fix the values a+ 0 = 1, a+ 1 = 1, a+ 2 = 0.02154486, a+ 3 = 2.6966906, a+ 4 = −0.1355754 · 10−3, a+ 5 = 0.3935304, a+ 6 = 0.4103265 · 10−5, b+ 1 = −3.351348, b+ 3 = −2.4176393, b+ 5 = −0.1400563,then M6 has only 8 real zeros at ρ = 1, 2, 3, . . . , 7, 8. Alternatively, choosing a+ n−1 = a+ n = 1, the zeros can be obtained bifurcating from the origin reordering in an adequateway the perturbation parameters. 5. Perturbations of any order for some classic Liénard families In the previous section we have seen that the study of the second Melnikov function for classical Liénard equations with a horizontal discontinuity line is more difficult than the first order study. A higher order analysis for perturbations of degree n can only be done adding some extra hypotheses or changing the slope of the discontinuity line. This is what we do in Theorems 1.3 and 1.4. In the first result we assume some symmetries on the perturbations and in the second we change the discontinuity line to the vertical axis. Their proofs follow directly from the next two propositions. Proposition 5.1. (a) When f ± i are even polynomials, system (3) has a center at the origin for every ε.Moreover, all the Melnikov functions vanish. (b) When f ± i are odd polynomials, the Melnikov function of order N, MN , associated to system (3) has at most (n − 1)/2 zeros. Moreover, there exist polynomials f ± i such that the zeros are realized as simple ones. Now we consider system (4). Under the change of coordinates {x = y, y = x} it becomes (x′, y′) = ( y,−x + N∑ i=1 εif ± i (y) ) , (54) where f ± i are polynomials of degree n defined now in Σ± 0 = {(x, y) ∈ R2 : ±y ≥ 0}. Observe that, for this new system, the discontinuity set is the x-axis. Proposition 5.2. All theMelnikov functions associated to system (54) have at most n zeros. Moreover, there exist polynomials f ± i such that all the zeros are realized as simple ones. Remark 5.3. In the above results every Melnikov function has exactly the same form. In fact all the higher order studies coincide with the first one.We can say that theMelnikov function stabilizes in the first step. Like in the Liénard smooth case, see [8]. Proof of Proposition 5.1. (a) Each systemdefined by (3), assuming the even property in f ±, is symmetricwith respect the y axis. In fact it is invariant via the change of variables (x, y, t) → (−x, y,−t). Hence, the system is globally invariant with respect to this change of variables and, consequently, the respective half-returnmaps are the identity ones for every ε: Π+(ρ) = Π−(ρ) = −ρ. Then the Melnikov functions vanish identically for each order and the proof finishes. (b) First, system (3), assuming the odd symmetry in f ± and writing f ± i (x) = xf̂ ± i (x2) = x ∑n j=0a ± 2j+1,i(x 2)j where a2j+1,i ∈ R, can be expressed as⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩ dH + N∑ i=1 εiω+ i = 0, if θ ∈ [0, π ), dH + N∑ i=1 εiω− i = 0, if θ ∈ [π, 2π ), (55) where −ω± i = α± i (r, θ )dr + β± i (r, θ ) with α± i (r, θ ) = n∑ j=0 a± 2j+1,ir 2j+1cos2j+1θ sin θ and β± i (r, θ ) = n∑ j=0 a± 2j+1,ir 2j+2cos2j+2θ. (56) In order to find the Melnikov functions of order i, Mi, we will write r±(θ, ε, ρ), the solution of (55), as r±(θ, ε, ρ) = N∑ i=0 εir± i (θ, ρ). Consequently, Mi(ρ) = ri(π, ρ) − r− i (−π, ρ). The first Melnikov function is given by M1 = r+ 1 (π ) − r− 1 (−π ) where, from Proposition 2.3, the functions r± 1 (θ ) := r± 1 (θ, ρ) satisfy ρr± 1 (θ ) = F± 1 (ρ)θ + S± 1 (ρ, θ ), (57) and F± 1 , S± 1 are obtained from the one-forms −Ω± 1 where h± 0 = 1 and −Ω± 1 = −ω± 1 h ± 0 = α± 1 (r, θ )dr + β± 1 (r, θ )dθ. With this notation, the functions F± 1 and S± 1 are given by the expressions F± 1 (r) = 1 2π ∫ 2π 0 β± 1 (r, θ )dθ and S± 1 (r, θ ) = ∫ θ 0 β± 1 (r, ψ)dψ − F± 1 (r)θ. So we get F± 1 (r) = n∑ j=0 a± 2j+1,1c2j+2r2j+2 and S± 1 (r, θ ) = n∑ j=0 a± 2j+1,1C2j+2(θ )r2j+2. (58) where, c2j+2 and C2j+2 are given in Lemma2.6 and the function C2j+2 has the form sin θ p2j+1(cos θ ) with p2j+1 a polynomial of degree 2j + 1 in cos θ.Moreover, h± 1 (r, θ ) = 1 r ( α± 1 (r, θ ) − ∂S± 1 (r, θ ) ∂r ) . (59) Observe that, from Eqs. (57) and (58), we have S± 1 (r,−θ ) = −S± 1 (r, θ ) and r± 1 (−θ ) = −r± 1 (θ ). C.A. Buzzi et al. / Physica D 371 (2018) 28–47 41 And, from the fact that S+ 1 (ρ, π ) = S− 1 (ρ,−π ) = 0, it follows from Eq. (57) that M1(ρ) = r+ 1 (π ) − r− 1 (−π ) = F+ 1 (ρ) + F− 1 (ρ) ρ π + S+ 1 (ρ, π ) − S− 1 (ρ,−π ) ρ = F+ 1 (ρ) + F− 1 (ρ) ρ π. Suppose that the first Melnikov function vanishes identically, that is,M1(ρ) ≡ 0. Then F+ 1 (ρ) = −F− 1 (ρ) and from (58) we obtain a− 2j+1,1 = −a+ 2j+1,1, j = 1, . . . , n. In particular, from Eqs. (56), (57), and (58) we get α− 1 (r, θ ) = − α+ 1 (r, θ ), β− 1 (r, θ ) = − β+ 1 (r, θ ), S+ 1 (r, θ ) = − S− 1 (r, θ ), r− 1 (θ ) = − r+ 1 (θ ). Also as α± 1 and S± 1 are odd functions with respect to the variable θ, it follows that h± 1 (r,−θ ) = −h± 1 (r, θ ). From the equalities above and the expressions (56) and (59) we can write F+ 1 (ρ) + F− 1 (ρ) = 0, S+ 1 (ρ, θ ) − S− 1 (ρ,−θ ) = 0, h+ 1 (ρ, θ ) − h− 1 (ρ,−θ ) = 0, α+ 1 (ρ, θ ) − α− 1 (ρ,−θ ) = 0, β+ 1 (ρ, θ ) + β− 1 (ρ,−θ ) = 0, r+ 1 (θ ) − r− 1 (−θ ) = 0. (60) Now taking the ε2-terms in the series expansion of Eq. (8) of Proposition 2.3 we get ρ r± 2 (θ ) = − 1 2 ( r± 1 (θ ) )2 + Oε (∫ θ 0 F± 1 (ρ + εr± 1 (ψ))dψ + S± 1 (ρ + εr± 1 (θ ), θ ) ) + F± 2 (ρ)θ + S± 2 (ρ, θ ), (61) where Oε(·) denotes the order ε terms of (·). The secondMelnikov function is given byM2 = r+ 2 (π )−r− 2 (−π ) and, from conditions (60), we have∫ −θ 0 F− 1 (ρ + εr− 1 (ψ))dψ = − ∫ θ 0 F− 1 (ρ + εr− 1 (−ψ))dψ = − ∫ θ 0 F− 1 (ρ + εr+ 1 (ψ))dψ = ∫ θ 0 F+ 1 (ρ + εr+ 1 (ψ))dψ. Also, S− 1 (ρ + εr− 1 (−θ ),−θ ) = −S− 1 (ρ + εr− 1 (−θ ), θ ) = −S− 1 (ρ + εr+ 1 (θ ), θ ) that implies S+ 1 (ρ + εr+ 1 (θ ), θ ) − S− 1 (ρ + εr− 1 (−θ ),−θ ) = S+ 1 (ρ + εr+ 1 (θ ), θ ) + S− 1 (ρ + εr+ 1 (θ ), θ ) = 0. Next, from (61), we write r+ 2 (θ ) − r− 2 (−θ ) = F+ 2 (ρ) + F− 2 (ρ) ρ θ + S+ 2 (ρ, θ ) − S− 2 (ρ,−θ ) ρ (62) where the functions F± 2 , S± 2 and h± 2 are obtained in a similar way as the case i = 1 but now using the one-forms −Ω± 2 = −ω± 2 − h± 1 ω ± 1 = α̂± 2 (r, θ )dr + β̂± 2 (r, θ )dθ, with α̂± 2 (r, θ ) = α± 2 (r, θ ) + h± 1 (r, θ )α ± 1 (r, θ ) = α± 2 (r, θ ) + α ±,(1) 2 (r, θ ), β̂± 2 (r, θ ) = β± 2 (r, θ ) + h± 1 (r, θ )β ± 1 (r, θ ) = β± 2 (r, θ ) + β ±,(1) 2 (r, θ ). (63) So, we have F± 2 (ρ) = F±,(2) 2 (ρ) + F±,(1) 2 (ρ), and S± 2 (ρ, θ ) = S±,(2) 2 (ρ, θ ) + S±,(1) 2 (ρ, θ ), where F±,(2) 2 , S±,(2) 2 depend on the coefficients of order ε2 of the one-forms (55) and are given by F±,(2) 2 (r) = 1 2π ∫ 2π 0 β± 2 (r, θ )dθ = n∑ j=0 a± 2j+1,2c2j+2r2j+1, S±,(2) 2 (r, θ ) = ∫ θ 0 β± 2 (r, ψ)dψ − F±,(2) 2 (r)θ = n∑ j=0 a± 2j+1,2C2j+2(θ )r2j+2, (64) where C2j+2 is given in Eq. (58). The functions F±,(1) 2 , S±,(1) 2 depend on the coefficients of order ε of the one-forms (55) and are given by F±,(1) 2 (r) = 1 2π ∫ 2π 0 β ±,(1) 2 (r, θ )dθ, S±,(1) 2 (r, θ ) = ∫ θ 0 β ±,(1) 2 (r, ψ)dψ − F±,(1) 2 (r)θ. Moreover, from the fact that β+,(1) 2 (r, θ ) = −β −,(1) 2 (r, θ ), these functions satisfy F+,(1) 2 (r) = −F−,(1) 2 (r), and S+,(1) 2 (r, θ ) = −S−,(1) 2 (r, θ ). (65) Now it follows from Eqs. (60) and (65) that S+,(1) 2 (r, θ ) − S−,(1) 2 (r,−θ ) = ∫ θ 0 β +,(1) 2 (r, ψ)dψ − ∫ −θ 0 β −,(1) 2 (r, ψ)dψ = ∫ θ 0 ( β +,(1) 2 (r, ψ) + β −,(1) 2 (r,−ψ) ) dψ = ∫ θ 0 ( β +,(1) 2 (r, ψ) + h− 1 (r,−ψ)β− 1 (r,−ψ) ) dψ = ∫ θ 0 ( β +,(1) 2 (r, ψ) − h+ 1 (r, ψ)β+ 1 (r, ψ) ) dψ = ∫ θ 0 ( β +,(1) 2 (r, ψ) − β +,(1) 2 (r, ψ) ) = 0. So we obtain from Eq. (62) the difference r+ 2 (θ ) − r− 2 (−θ ) = F+ 2 (ρ) + F− 2 (ρ) ρ θ + S+ 2 (ρ, θ ) − S− 2 (ρ, θ ) ρ = F+,(2) 2 (ρ) + F−,(2) 2 (ρ) ρ θ + S+,(2) 2 (ρ, θ ) − S−,(2) 2 (ρ,−θ ) ρ 42 C.A. Buzzi et al. / Physica D 371 (2018) 28–47 and, from the fact that S+,(2) 2 (ρ, π ) = S−,(2) 2 (ρ,−π ) = 0, we get the expression for the second Melnikov function M2(ρ) = r+ 2 (π ) − r− 2 (−π ) = F+,(2) 2 (ρ) + F−,(2) 2 (ρ) ρ π + S+,(2) 2 (ρ, π ) − S−,(2) 2 (ρ,−π ) ρ = F+,(2) 2 (ρ) + F−,(2) 2 (ρ) ρ π. (66) Now suppose that it vanishes identically, that is, M2(ρ) ≡ 0. Then from Eq. (66) we have F+,(2) 2 (ρ) = −F−,(2) 2 (ρ) and from Eqs. (64) the next relation holds: a− 2j+1,2 = −a+ 2j+1,2, j = 1, . . . , n. From Eq. (64) it follows that S+,(2) 2 (ρ, θ ) = −S−,(2) 2 (ρ, θ ). Also, using Eq. (65), we obtain F+ 2 (ρ) = F+,(2) 2 (ρ) + F+,(1) 2 (ρ) = −(F−,(2) 2 (ρ) + F−,(2) 2 (ρ)) = −F− 2 (ρ), and similarly S+ 2 (ρ, θ ) = −S− 2 (ρ, θ ). This implies that S+ 2 (ρ, θ ) − S− 2 (ρ,−θ ) = S+ 2 (ρ, θ ) + S− 2 (ρ, θ ) = 0. Now, from the form of α̂± 2 and β̂± 2 given in (63), we have α̂+ 2 (θ ) − α̂− 2 (−θ ) = α+ 2 (r, θ ) + h+ 1 (r, θ )α + 1 (r, θ ) − ( α− 2 (r,−θ ) + h− 1 (r,−θ )α − 1 (r,−θ ) ) = α+ 2 (r, θ ) + h+ 1 (r, θ )α + 1 (r, θ ) − ( α+ 2 (r, θ ) + h+ 1 (r, θ )α + 1 (r, θ ) ) = 0. Similarly we obtain β̂+ 2 (θ ) + β̂− 2 (−θ ) = 0 and, as h± 2 (r, θ ) = 1 r ( α̂± 2 (r, θ ) − ∂S± 2 ∂r ) , we get h+ 2 (r, θ ) − h− 2 (r,−θ ) = 0. Using all the above relations together we can write F+ 2 (ρ) + F− 2 (ρ) = 0, S+ 2 (ρ, θ ) − S− 2 (ρ,−θ ) = 0, h+ 2 (ρ, θ ) − h− 2 (ρ,−θ ) = 0, α+ 2 (ρ, θ ) − α− 2 (ρ,−θ ) = 0, β+ 2 (ρ, θ ) + β− 2 (ρ,−θ ) = 0, r+ 2 (θ ) − r− 2 (−θ ) = 0. Now suppose that for j = 1, . . . ,m − 1 we have F+ j (ρ) + F− j (ρ) = 0, S+ j (ρ, θ ) − S− j (ρ,−θ ) = 0, h+ j (ρ, θ ) − h− j (ρ,−θ ) = 0, α+ j (ρ, θ ) − α− j (ρ,−θ ) = 0, β+ j (ρ, θ ) + β− j (ρ,−θ ) = 0, r+ j (θ ) − r− j (−θ ) = 0. (67) Then the terms of order εm of Eq. (8) in Proposition 2.3 can be written as ρr± m (θ ) = F± m (ρ)θ + S± m (ρ, θ ) + Oεm × ( − 1 2 ( ρ + εr± 1 (θ ) + · · · + εm−1r± m−1(θ ) )2) + Oεm ( m−1∑ i=1 εi ∫ θ 0 F± i ( ρ + εr± 1 (ψ) + · · · + εm−1r± m−1(ψ) ) dψ ) + Oεm ( m−1∑ i=1 S± i ( ρ + εr± 1 (θ ) + · · · + εm−1r± m−1(θ ) )) . Now observe that as in the case i = 2 under the hypothesis of induction (67) we have r+ m (θ ) − r− m (−θ ) = F+ m (ρ) + F− m (ρ) ρ θ + S+ m (ρ, θ ) − S− m (ρ, θ ) ρ , (68) where the functions F± m , S ± m and h± m are obtained from the one- forms −Ω± m = − m∑ j=1 ω± j h ± m−j = −ω± m − m−1∑ j=1 ω± j h ± m−j = α̂± m (r, θ )dr + β̂± m (r, θ )dθ, with α̂± m (r, θ ) = α± m (r, θ ) + m∑ j=1 h± m−jα ± j (r, θ ) = α± m (r, θ ) + α±,(m−1) m (r, θ ), β̂± m (r, θ ) = β± m (r, θ ) + m∑ j=1 h± m−jβ ± j (r, θ ) = β± m (r, θ ) + β±,(m−1) m (r, θ ). (69) Furthermore, it follows that α+,(m−1) m (r, θ ) − α−,(m−1) m (r,−θ ) = m−1∑ j=1 ( h+ m−j(r, θ )α + j (r, θ ) − h− m−j(r,−θ )α − j (r,−θ ) ) = m−1∑ j=1 ( h+ m−j(r, θ )α + j (r, θ ) − h+ m−j(r, θ )α + j (r, θ ) ) = 0. Similarly, we obtain the vanishing relation β +,(m−1) m (r, θ ) + β −,(m−1) m (r,−θ ) = 0.Moreover, we can write F± m (r) = F±,(m) m (r) + F±,(m−1) m (r), and S± m (r, θ ) = S±,(m) m (r, θ ) + S±,(m−1) m (r, θ ), with F±,(m) m , S±,(m) m depending on the coefficients of order εm of the one-forms (55) and they are given by F±,(m) m (r) = − 1 2π ∫ 2π 0 β± m (r, θ )dθ = n∑ j=0 a± 2j+1,mc2j+2r2j+1, S±,(m) m (r, θ ) = ∫ θ 0 β± m (r, ψ)dψ − F±,(m) m (r)θ = n∑ j=0 a± 2j+1,mC2j+2(θ )r2j+2. (70) Here the functions F±,(m−1) m , S±,(m−1) m depend on the coefficients of order εj of the one-form (55), for every j = 1, . . . ,m − 1, and from the fact that β+,(m−1) m (r, θ ) = −β −,(m−1) m (r, θ ) these functions satisfy F+,(m−1) m (r) = −F−,(m−1) m (r), and S+,(m−1) m (r, θ ) = −S−,(m−1) m (r, θ ). (71) Now as S±,(m−1) m (r, θ ) = ∫ θ 0 β±,(m−1) m (r, ψ)dψ − F±,(m−1) m (r)θ, C.A. Buzzi et al. / Physica D 371 (2018) 28–47 43 we get, also using the symmetry of the functions β , S+,(m−1) m (r, θ ) − S−,(m−1) m (r,−θ ) = ∫ θ 0 β+,(m−1) m (r, ψ)dψ − ∫ −θ 0 β−,(m−1) m (r, ψ)dψ = ∫ θ 0 ( β+,(m−1) m (r, ψ) −β+,(m−1) m (r, ψ) ) dψ = 0. Then, from Eq. (68), we have r+ m (θ ) − r− m (−θ ) = F+ m (ρ) + F− m (ρ) ρ θ + S+ m (ρ, θ ) − S− m (ρ,−θ ) ρ = F+,(m) m (ρ) + F−,(m) m (ρ) ρ θ + S+,(m) m (ρ, θ ) − S−,(m) m (ρ,−θ ) ρ and, from the fact that S+,(m) m (ρ, π ) = 0 = S−,(m) m (ρ,−π ), we obtain Mm = r+ m (π ) − r− m (−π ) = F+,(m) m (ρ) + F−,(m) m (ρ) ρ π + S+,(m) m (ρ, π ) − S−,(m) m (ρ,−π ) ρ = F+,(m) m (ρ) + F−,(m) m (ρ) ρ π. This implies that, whenMm(ρ) ≡ 0,we have a− 2j+1,m = −a+ 2j+1,m, for j = 1, . . . ,m. From Eq. (70) it follows that S+,(m) m (ρ, θ ) = −S−,(m) m (ρ, θ ). Also, from Eq. (71), we obtain F+ m (ρ) = F+,(m) m (ρ) + F+,(m−1) m (ρ) = − ( F−,(m) m (ρ) + F−,(m−1) m (ρ) ) = −F− m (ρ), and similarly S+ m (ρ, θ ) = −S− m (ρ, θ ). This implies the next vanish- ing relation: S+ m (ρ, θ ) − S− m (ρ,−θ ) = S+ m (ρ, θ ) + S− m (ρ, θ ) = 0. Consequently, from Eq. (68), the next difference also vanishes: r+ m (θ ) − r− m (−θ ) = 0. Moreover, from the form of α̂± m and β̂± m given by (69), we have α̂+ m (θ ) − α̂− m (−θ ) = α+ m (r, θ ) + m−1∑ j=1 h+ j (r, θ )α + j (r, θ ) − ( α− m (r,−θ ) + m−1∑ j=1 h− j (r,−θ )α − j (r,−θ ) ) = α+ m (r, θ ) + m−1∑ j=1 h+ j (r, θ )α + j (r, θ ) − ( α+ m (r, θ ) + m−1∑ j=1 h+ j (r, θ )α + j (r, θ ) ) = 0. Similarly, we can obtain also that β̂+ m (θ ) + β̂− m (−θ ) = 0 and, as h± m(r, θ ) = 1 r ( α̂± m − ∂S± m ∂r ) ,we get h+ m(r, θ ) − h− m(r,−θ ) = 0. In this way we obtain Eqs. (67) for j = m and the result follows by induction. □ Proof of Proposition 5.2. Using the tools described in Section 2 and writing f ± i (y) = ∑n j=0a ± j,iy j, system (54) moves to its equiva- lent piecewise differential one-form⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩ dH + N∑ i=0 εiω+ i = 0, if θ ∈ [0, π ), dH + N∑ i=0 εiω− i = 0, if θ ∈ [π, 2π ), where H = r2/2 and −ω± i = α± i (r, θ )dr + β± i (r, θ )dθ,with α± i (r, θ ) = n∑ j=0 a± j,ir jsinjθ cos θ, and β± i (r, θ ) = − n∑ j=0 a± j,ir j+1sinj+1θ. Using these expressions we obtain F± 1 (r) = − [(n−1)/2]∑ j=0 a± 2j+1,1s2j+2r2j+2, S± 1 (r, θ ) = − [(n−1)/2]∑ j=0 a± 2j,1S2j+1(θ )r2j+1 − [(n−1)/2]∑ j=0 a± 2j+1,1S2j+2(θ )r2j+2. (72) Moreover, by Lemma 2.6, we can write also S2j+2(θ ) = sin θ p2j+1(cos θ ) and S2j+1(θ ) = p̃2j+1(cos θ ) for some p2j+1 and p̃2j+1 polynomials of degree 2j + 1 in cos θ. Now, as ρr± 1 (θ ) = F± 1 (ρ)θ + S± 1 (ρ, θ ), it follows from this equation and the previous functions that the first Melnikov function,M1(ρ) = r+ 1 (π ) − r− 1 (−π ),writes as M1(ρ) = F+ 1 (ρ) + F− 1 (ρ) ρ π + S+ 1 (ρ, π ) − S− 1 (ρ,−π ) ρ = [(n−1)/2]∑ j=0 (a− 2j,1 − a+ 2j,1)S2j+1(π )ρ2j − [(n−1)/2]∑ j=0 (a+ 2j+1,1 + a− 2j+1,1)s2j+2ρ 2j+1. In this case, M1(ρ) is a polynomial of degree n with monomials depending on the coefficients of the ε-perturbation of system (54). So (54) has at most n limit cycles bifurcating from the origin up to a first order study. Moreover, as all the coefficients of M1(ρ) are linearly independent, we can chose a system such that the corresponding functionM1 has exactly n simple zeros. In order to study the second Melnikov function we suppose M1(ρ) ≡ 0. In this case we get a+ 2j+1,1 = −a− 2j+1,1, and a+ 2j,1 = a− 2j,1. (73) Using (72), (73), and also the definition h± 1 (r, θ ) = 1 r ( α± 1 − ∂S± 1 ∂r (r, θ ) ) ,we can easily see that F+ 1 (ρ) + F− 1 (ρ) = 0, S+ 1 (ρ, θ ) − S− 1 (ρ,−θ ) = 0, α+ 1 (ρ, θ ) − α− 1 (ρ,−θ ) = 0, β+ 1 (ρ, θ ) + β− 1 (ρ,−θ ) = 0, r+ 1 (θ ) − r− 1 (−θ ) = 0, h+ 1 (ρ, θ ) − h− 1 (ρ,−θ ) = 0. 44 C.A. Buzzi et al. / Physica D 371 (2018) 28–47 Taking into account the previous equalities and using the same induction procedure used for proving Theorem1.3(b),we can show that Mk(ρ) has the same form of M1(ρ) for every k ∈ N. The proof of this result is quite similar to the proof of Theorem 1.3 and it is omitted. □ 6. First order perturbations of Liénard families of low degree In this section we present two families for which the number of limit cycles increase with the order of perturbation but not as fast as the upper bound given in Theorem 1.1. As an example of the difficulties found during the implementation of this procedure we will present here only different families of lower degree. In the first one, the number of limit cycles with a second order study is bigger than the first order, but we have no more limit cycles with a third order or higher. For the second one, the number of limit cycles increases with the order, but also in this case this growth stops at a given order, fourth for this family. Typically, for a fixed family, a kind of saturation growth occurs. This phenomenon was also showed for the piecewise linear family in [11]. It is clear that which is this maximum number and for which order it is given are very difficult questions that involve too many computations in our approach. The first family has been studied up to order 6 and the second one only up to order 5. We have not gone further because of the computational difficulties. Proposition 6.1. For n = 1, 2, 3, 4, consider the polynomial system, inΣ± 0 , X± : (x′, y′) = ( −y + N∑ i=1 εif ± i,n(x), x ) , (74) where f ± i,n are real polynomials of degree n. Then M1 has at most [(n − 1)/2] simple zeros and MN has at most n + [(n − 1)/2] simple zeros for N = 2, . . . , 6. Moreover there exist f ± i,n such that system (74), for ε small enough, has exactly these number of limit cycles. Proof. The proof is done in a case by case study. For each degree n, using the procedure described in Section 2, we obtain that MN (ρ) = Pn,N (ρ)/ρdN where Pn,N are polynomials of degree mn,N for N = 1, 2, . . . , 6. For simplicity we have not indicated, in Pn,N , the dependence on the coefficients of the polynomials f ± i,n.Here the values of dN are 0, 0, 1, 2, 2, 4 for N = 1, 2, . . . , 6 and the values ofmn,N are given in the next table: n N 1 2 3 4 5 6 1 1 1 2 3 3 5 2 1 2 4 6 7 10 3 3 5 8 11 13 17 4 3 6 10 14 17 22 For the case corresponding to degree n = 2, the first four polynomials P2,N , taking f ± i,2(x) = a± 0,i + a± 1,ix + a± 2,ix 2, are P2,1 = π 2 ( a+ 11 + a− 11 ) ρ, P2,2 = 2 9 ( a+ 11a + 21 − a− 11a − 21 ) ρ2 − π 8 ( π ((a− 11) 2 − (a+ 11) 2) − 4(a− 12 + a+ 12) ) ρ + 2 ( a+ 01a + 11 − a− 01a − 11 ) , P2,3 = 1 24 π ( a− 11(a − 21) 2 + a+ 11(a + 21) 2)ρ4 + 1 18 ( 3π ((a+ 11) 2a+ 21 + (a− 11) 2a− 21) + 4(a+ 11a + 22 − a− 11a − 22) + 4(a+ 12a + 21 − a− 12a − 21) ) ρ3 + 1 48 π ( π2((a− 11) 3 + (a+ 11) 3) − 12π (a− 11a − 12 − a+ 11a + 12) + 3((a− 11) 3 + (a+ 11) 3) + 24(a− 13 + a+ 13) ) ρ2 + 1 2 ( π (a+ 01(a + 11) 2 + a− 01(a − 11) 2) + 4(a+ 01a + 12 + a+ 02a + 11 − a− 01a − 12 − a− 02a − 11) ) ρ + 1 2 π ( (a− 01) 2a− 11 + (a+ 01) 2a+ 11 ) , P2,4 = − 206 2025 ( a− 11(a − 21) 3 − a+ 11(a + 21) 3)ρ6 − 1 648 ( 27(a− 11a − 21 − a+ 11a + 21)(a − 11a − 21 + a+ 11a + 21)π 2 − 27(2(a− 11a − 21a − 22 + a+ 11a + 21a + 22) + a− 12(a − 21) 2 + a+ 12(a + 21) 2)π + 32(a− 11a − 21 − a+ 11a + 21)× (a− 11a − 21 + a+ 11a + 21) ) ρ5 − 1 648 ( 45((a− 11) 3a− 21 − (a+ 11) 3a+ 21)π 2 − 108((a− 11) 2a− 22 + (a+ 11) 2a+ 22 + 2(a− 11a − 12a − 21 + a+ 11a + 12a + 21))π + 16(9(a− 01a − 11(a − 21) 2 − a+ 01a + 11(a + 21) 2) + 9(a− 11a − 23 − a+ 11a + 23) + 9(a− 12a − 22 − a+ 12a + 22) + 9(a− 13a − 21 − a+ 13a + 21) + 2((a− 11) 3a− 21 − (a+ 11) 3a+ 21)) ) ρ4 − 1 1152 ( 3(a− 11 − a+ 11)(a − 11 + a+ 11)((a − 11) 2 + (a+ 11) 2)π4 − 72((a− 11) 2a− 12 + (a+ 11) 2a+ 12)π 3 + 36((a− 11) 4 − (a+ 11) 4 + 8(a− 11a − 13 − a+ 11a + 13) + 4((a− 12) 2 − (a+ 12) 2))π2 − 72(3((a− 11) 2a− 12 + (a+ 11) 2a+ 12) + 8(a− 14 + a+ 14))π + 512(a− 01(a − 11) 2a− 21 − a+ 01(a + 11) 2a+ 21) ) ρ3 − 1 24 ( − 3(a− 01(a − 11) 3 − a+ 01(a + 11) 3)π2 + 12(2(a− 01a − 11a − 12 + a+ 01a + 11a + 12) + a− 02(a − 11) 2 + a+ 02(a + 11) 2)π + 16((a− 01) 2a− 11a − 21 − (a+ 01) 2a+ 11a + 21 + 3(a+ 01a + 13 − a− 01a − 13 + a+ 02a + 12 − a− 02a − 12 + a+ 03a + 11 − a− 03a − 11)) ) ρ2 − π 2 ( (a− 01) 2a− 12 + (a+ 01) 2a+ 12 + 2(a− 01a − 02a − 11 + a+ 01a + 02a + 11) ) ρ − 2 3 ( (a− 01) 3a− 11 − (a+ 01) 3a+ 11 ) . We have written the parameters ai,j as aij so write in short. More- over, the polynomials P2,5 and P2,6 are not written here because of the size of them. Instead of P2,4 has 62 monomials, P2,5 and P2,6 have 152 and 350, respectively. Clearly, P2,1 has no positive zeros, consequently there are no limit cycles bifurcating from the period annulus up to order 1. Let us continue assuming M1(ρ) ≡ 0. Under this assumption, P2,2 is a polynomial of degree 2 with arbitrary coefficients. Then there exists a choice of perturbed coefficients such that P2,2, and also M2, has exactly 2 simple zeros and from the period annulus we can get 2 limit cycles for system (74), for ε small enough. When M1(ρ) ≡ M2(ρ) ≡ 0, then all the coefficients of M3 of degree greater than 2 vanish. The same simplification occurs for M4, M5, andM6 when the previousMN vanish identically. Consequently, for ε small enough, no more than 2 limit cycles can bifurcate from the period annulus up to orders 3,4,5, and 6. It can be easily checked that for system X+ :(x′, y′) = ( −y + (9 2 x2 + x + e0 2 ) ε + 2e1xε2, x ) , X− :(x′, y′) = (−y − xε, x), C.A. Buzzi et al. / Physica D 371 (2018) 28–47 45 defined in Σ± 0 , we have M1(ρ) ≡ 0 and M2(ρ) = ρ2 + πe1ρ + e0. Using the same simplification procedure, the degrees of the polynomials Pn,N decrease. Next table shows which they are after the simplifications: n N 1 2 3 4 5 6 1 1 1 1 1 1 1 2 1 2 2 2 2 2 3 3 4 4 4 4 4 4 3 6 6 6 6 6 We remark that the above table does not give the number of limit cycles, only the degrees of Pn,N after the simplification procedure. The number of limit cycles given in the statement follows studying the independence of the parameters of Pn,N for each degree n and orderN . We remark that although P4,N has degree 6, the coefficient of degree 5 vanishes, consequently M2 can have only 5 positive simple zeros. Finally, we show explicit systems of degrees 1, 3 and 4 such that the corresponding functions M1 are identically zero and the functionsM2 have 1, 4 and 5 positive simple zeros, respectively. Using the described procedure we can prove that the linear system X+ :(x′, y′) = (−y + ε(e0x + 1) + ε2x, x), X− :(x′, y′) = (−y − εe0x, x), defined inΣ± 0 , hasM1(ρ) ≡ 0 andM2(ρ) = πρ/2+2e0. Similarly, the cubic system X+ :(x′, y′) = ( − y + ( 1 63e0 x3 + 135 2 e0x2 + x − 945 2 e20 + 63 2 e2e0 ) ε + 2e1xε2, x ) , X− :(x′, y′) = ( −y − ε ( 1 63e0 x3 + x ) + 8 3 e3x3ε2, x ) , defined inΣ± 0 , has M1(ρ) ≡ 0 and M2(ρ) = ρ4 + πe3ρ3 + e2ρ2 + πe1ρ − 63e0(15e0 − 2e2). We remark that there are values of ei, i = 0, . . . , 3 such thatM2 has 4 positive simple zeros. For example for e0 = 1, e1 = −417/(31π ), e2 = 467/31, e3 = −207/(31π ) we haveM4(ρ) = (ρ − 1)(ρ − 2)(ρ − 3)(ρ − 21/31). Finally, the quartic system X+ :(x′, y′) = ( − y + (525 58 x4 + x3 + e0x ) ε + (8 3 e3x3 + 2e1x ) ε2, x ) , X− :(x′, y′) = ( − y + ( − x3 + (555 58 e0 + 15 14 e4 ) x2 − e0x − 185 174 e20 − 5 42 e0e4 + 1 2 e2 ) ε, x ) , defined inΣ± 0 , has M1(ρ) ≡ 0 and M2(ρ) = ρ6 + e4ρ4 + πe3ρ3 + e2ρ2 + πe1ρ − e0(1295e20 + 145e0e4 − 609e2)/609. In general the above polynomial can have at most 5 positive zeros. Moreover, there exist values of e0, . . . , e4 such that it has 5 positive simple zeros. For example, when e0 = 3/2, e1 = −830235/(841π ), e2 = 644273/1682, e3 = 45885/(841π ), e4 = −68047/1682, we have that M2 has five positive simple zeros at 1, 2, 3, 4,−5 + √ 1534/58 and one negative at −5 − √ 1534/58. Consequently, there exist values of the perturbation parameters such that system (74), for ε small enough and n = 1, 2, 3, 4, has 1, 2, 4, 5 limit cycles. □ Proposition 6.2. For n = 1, 2, 3, 4, consider the polynomial system X± : (x′, y′) = ( −y + f ±(x, ε), x + g±(x, ε) ) = ( −y + N∑ i=1 εif ± i,n(x), x + N∑ i=1 εig± i,n(x) ) , (75) defined in Σ± 0 , where f ± i,n and g± i,n are polynomials of degree n. Then the functions MN have at most n, 2n − 1, 2n, 3n, 3n simple zeros for N = 1, 2, . . . , 5. Moreover there exist f ± i,n and g± i,n such that system (75), for ε small enough, has these number of limit cycles. Proof. The proof follows the same procedure aswe have described in the proof of Proposition 6.1. In fact the Melnikov functions take the same form, MN (ρ) = Pn,N (ρ)/ρdN , but with different degrees and coefficients. In this case the values of dN are 0, 0, 1, 2, 3 for N = 1, 2, . . . , 5 and the values of the degrees of the polynomials Pn,N are mn,1 = n andmn,N = Nn − 1 for N = 2, . . . , 5. We note that the maximum number of limit cycles appears up to order N = 4. Now we will only provide explicit examples exhibiting this maximum number of limit cycles, for each degree n in the statement. The studies of lower order can be done similarly. The linear Liénard system (75) with N = 4 given by X+ :(x′, y′) = ( −y + (e0 + x)ε + 2xε4, x − e1 2e0 ε ) , X− :(x′, y′) = ( − y + (e0 − x)ε + 1 2 (−4e0 + e2)ε3, x − e1 2e0 ε + 2e0ε2 ) , has M1(ρ) ≡ M2(ρ) ≡ M3(ρ) ≡ 0 and M4(ρ) = 1 ρ2 ( πρ3 + e2ρ2 + πe1ρ + 4 3 e30 ) . The quadratic system (75) with N = 4 and f +(x, ε) = (x2 + x + e0)ε + 81e4 + 180e0 + 4 18 x2ε3, g+(x, ε) = x2ε2 + 6e5x2ε3 + 9e2 − 4e20 − 36e0 18 ε4, f −(x, ε) = 1 3 (−5x2 − 3x + 3e0)ε − e1 e0 ε2 − 20e0 3 xε3 + 2(3e0e3 + 5e1) 3e0 xε4, g−(x, ε) = 1 9 (7x2 + 18e0)ε2 + 1 e0 (6e0e5x2 − e1)ε3, has M1(ρ) ≡ M2(ρ) ≡ M3(ρ) ≡ 0 and M4(ρ) = 1 ρ2 (19292 54675 ρ6 + πe5ρ5 + e4ρ4 + πe3ρ3 + e2ρ2 + πe1ρ + 4 3 e30 ) . The cubic system (75) with N = 4 and X+ : ⎧⎨⎩x′ = −y + (a8x2 + x + a0)ε + ( 53 120 x3 + a7x ) ε2, y′ = x + x2ε3, X− : ⎧⎪⎪⎪⎨⎪⎪⎪⎩ x′ = −y + ( − 571 1525 a8x2 − x + a0 ) ε + A1(x)ε2 + A2(x)ε3 + a5x3ε4, y′ = x + x3ε + ( 318 1525 a8x2 + 2a0 ) ε2 + A3(x)ε3 + a4x2ε4, 46 C.A. Buzzi et al. / Physica D 371 (2018) 28–47 where A1(x) = − 163 120 x3 + a6x2 − a7x + a1, A2(x) = ( − 343228 2325625 a28 − 11 12 a7 ) x3 − 2284 1525 a0a8x + a2, A3(x) = ( 318 1525 a7a8 + 42 5 a0 + 1 3 a6 + 1 ) x2 + a3x + 2a0a7 + a1, has M1(ρ) ≡ M2(ρ) ≡ M3(ρ) ≡ 0. And choosing a0 = e0, a1 = − 4 3 e1 a0 , a2 = 2 3 e2 + 318 1525 a20a8 − a1a7 − 2a0, a3 = 16 3 e3 − 4568 1525 a0a7a8 + 8a6a0 − 2284 1525 a1a8, a4 = −2e4 + 2 3 a0a28 + 23 4 a0a7 + 1 3 a7a6 + 163 40 a1 + 106 1525 a8, a5 = 32 9 e5 − 343228 2325625 a7a28 − 64316 22875 a0a8 − 376 13725 a8a6 + 636 1525 a8 + 1 72 , a6 = 48 61 e6 − 24 61 + 353934 2325625 a7a8 − 201099590064 5408531640625 a38 − 159 305 a0, a7 = − 2048 159 e7 + 3439367744 5546615625 a28, a8 = 448350000 150098827 e8, we obtain M4(ρ) = 1 ρ2 ( 13863 655360 πρ9 + e8ρ8 + πe7ρ7 + e6ρ6 + πe5ρ5 + e4ρ4 + πe3ρ3 + e2ρ2 + πe1ρ + e30 ) . The quartic system (75) with N = 4 and X+ : ⎧⎪⎪⎪⎨⎪⎪⎪⎩ x′ = −y + (a10x2 + a0 + x)ε + ( a9x4 + 4021 348 a11x3 + a1 ) ε2 + a6x4ε3 + a3xε4, y′ = x + A1(x)ε − 617 58 a11x3ε2 + A2(x)ε3 + ε4a5x3, X− : ⎧⎨⎩x′ = −y + A3(x)ε + − 617 58 a11x3ε2 + A4(x)ε3 + a5x3ε4, y′ = x + A5(x)ε2 + A6(x)ε3, with A1(x) = x4 + (11177145 373984 a211 + 2551 403 a10 ) x2 − x + a0, A2(x) = (41642856783675 139864032256 a411 + 1250598335 9419722 a10a211 + 21416500 1461681 a210 + 16 3 a0 − 11 12 a8 ) x3 + (11177145 93496 a0a211 + 10204 403 a0a10 ) x, A3(x) = x4 + (11177145 373984 a211 + 2551 403 a10 ) x2 − x + a0, A4(x) = (41642856783675 139864032256 a411 + 1250598335 9419722 a10a211 + 21416500 1461681 a210 + 16 3 a0 − 11 12 a8 ) x3 + 1 93496 a0 ( 11177145a211 + 2367328a10 ) x, A5(x) = 33 70 x4 + a8x3 + (3725715 373984 a211 + 2954 1209 a10 ) x2 + 2a0, A6(x) = (16091363085 21691072 a311 + 129441493 701220 a10a11 − 37 15 a9 ) x4 + (7375 116 a0a11 + a7 ) x2 + a1, has M1(ρ) ≡ M2(ρ) ≡ M3(ρ) ≡ 0. Choosing the values a0 = e0, a1 = e1 e0 , a2 = 1 2 e2 − 2a0 + 2954 1209 a20a10 + 3725715 373984 a20a 2 11, a3 = 2e3 + 11177145 186992 a1a211 + 5102 403 a1a10, a4 = 3 2 e4 − 2954 3627 a10 + 1 5 a20 − 1241905 373984 a211 − 2 3 a0a210 − 3673 116 a1a11 − 23 4 a0a8, a5 = 8 3 e5 + 41269745055 21691072 a311a0 + 916501 2418 a10a0a11 + 3725715 186992 a211a7 + 8 3 a0a9 + 5908 1209 a10a7 + 4a1 − 1 72 a11, a6 = 75 74 e6 + 277 111 − 957035384280675 28182602499584 a10a411 − 246061222155 60738367456 a210a 2 11 + 58101015 1495936 a211a8 + 40871 4836 a10a8 + 17041717574 65385376173 a310 − 25 37 a11a7 − 1397742487944051375 26153455119613952 a611 − 56001305 13837408 a0a211 + 686648 313131 a0a10, a7 = − 12 5 e7 − 416237010695545095 32448455483392 a511 − 993055562624097 174830040320 a311a10 + 82710873 3739840 a211a9 − 1751870985349 2825916600 a210a11 − 7375 116 a0a11 + 242861 60450 a10a9 + 5553 928 a11a8, a8 = 22050 43123 e8 + 408588973833 812413343612 a10a211 − 64879 43123 a9a11 − 1372 43123 a0 + 14949785585339325 1507839165743872 a411 − 12575899204 21010689921 a210, a9 = − 320 33 e9 − 1724620185 834272 a311 − 2467922575 4984056 a10a11, a10 = 1799445375 1023106108 e10 − 1075494673575 237360617056 a211, a11 = − 556800 211631 e11, we have that M4(ρ) = 1 ρ2 ( 328474 6670125 ρ12 + πe11ρ11 + e10ρ10 + πe9ρ9 + e8ρ8 + πe7ρ7 C.A. Buzzi et al. / Physica D 371 (2018) 28–47 47 + e6ρ6 + πe5ρ5 + e4ρ4 + πe3ρ3 + e2ρ2 + πe1ρ + 4 3 e30 ) . For all the studied cases, n = 1, 2, 3, 4, the first three Melnikov functions vanish identically and the fourth,M4, can bewrittenwith arbitrary coefficients, ej. Consequently up to a study of fourth order and for ε small enough, system (75) has 3n limit cycles, bifurcating from some chosen 3n circumferences, for n = 1, 2, 3, 4. □ Acknowledgments This work has been realized thanks to the MINECO grants MTM2013-40998-P and MTM2016-77278-P (FEDER), the AGAUR grant 2014 SGR568, the European Community grants FP7-PEOPLE- 2012-IRSES 316338 and 318999, and the Brazilian FAPESP grants 2012/18780-0, 2013/24541-0 and 2017/03352-6. References [1] M. di Bernardo, C.J. Budd, A.R. Champneys, P. Kowalczyk, Piecewise-Smooth Dynamical Systems, in: Applied Mathematical Sciences, vol. 163, Springer- Verlag London, Ltd., London, 2008. Theory and applications. [2] A.F. Filippov, Differential Equations with Discontinuous Righthand Sides, in: Mathematics and its Applications (Soviet Series), vol. 18, Kluwer Academic Publishers Group, Dordrecht, 1988. [3] A.A. Andronov, A.A. Vitt, 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] I.D. Iliev, The number of limit cycles due to polynomial perturbations of the harmonic oscillator, Math. Proc. Cambridge Philos. Soc. 127 (2) (1999) 317– 322. [5] A. Buică, On the equivalence of the Melnikov functions method and the averaging method, Qual. Theory Dyn. Syst. 16 (3) (2017) 547–560. [6] M. Han, V.G. Romanovski, X. Zhang, Equivalence of the Melnikov function method and the averaging method, Qual. Theory Dyn. Syst. 15 (2) (2016) 471– 479. [7] A. Buică, J. Giné, J. Llibre, Bifurcation of limit cycles from a polynomial degen- erate center, Adv. Nonlinear Stud. 10 (3) (2010) 597–609. [8] A. Gasull, J. Torregrosa, A relation between small amplitude and big limit cycles, Rocky Mountain J. Math. 31 (4) (2001) 1277–1303. [9] E. Freire, E. Ponce, F. Rodrigo, F. Torres, Bifurcation sets of continuous piecewise linear systemswith two zones, Internat. J. Bifur. Chaos Appl. Sci. Engrg. 28 (11) (1998) 2073–2097. [10] J.C. Medrado, J. Torregrosa, Uniqueness of limit cycles for sewing planar piece- wise linear systems, J. Math. Anal. Appl. 431 (1) (2015) 529–544. [11] C. Buzzi, C. Pessoa, J. Torregrosa, Piecewise linear perturbations of a linear center, Discrete Contin. Dyn. Syst. 33 (9) (2013) 3915–3936. [12] Y. Tang, J. Llibre, Limit cycles of discontinuous piecewise quadratic and cubic polynomial perturbations of a linear center, Preprint, 2016. [13] J. Llibre, A.C. Mereu, Limit cycles for discontinuous generalized Lienard poly- nomial differential equations, Electron. J. Differential Equations 195 (8) (2013). [14] A. Lins, W. de Melo, C.C. Pugh, On Liénard’s equation, in: Geometry and Topol- ogy (Proc. III Latin Amer. School of Math., Inst. Mat. Pura Aplicada CNPq, Rio de Janeiro, 1976), in: Lecture Notes in Math., vol. 597, Springer, Berlin, 1977, pp. 335–357. [15] C. Zuppa, Order of cyclicity of the singular point of Liénard’s polynomial vector fields, Bol. Soc. Brasil. Mat. 12 (2) (1981) 105–111. [16] C. Li, J. Llibre, Uniqueness of limit cycles for Liénard differential equations of degree four, J. Differential Equations 252 (4) (2012) 3142–3162. [17] P. De Maesschalck, F. Dumortier, Classical Liénard equations of degree n ≥ 6 can have [ n−1 2 ] + 2 limit cycles, J. Differential Equations 250 (4) (2011) 2162– 2176. [18] F. Dumortier, D. Panazzolo, R. Roussarie, More limit cycles than expected in Liénard equations, Proc. Amer. Math. Soc. 135 (6) (2007) 1895–1904. [19] B. Coll, A. Gasull, R. Prohens, First Lyapunov constants for non-smooth Liénard differential equations, in: Proceedings of the 2nd Catalan Days on Applied Mathematics, (Odeillo, 1995), in: Collect. Études, Presses Univ. Perpignan, Perpignan, 1995, pp. 77–83. [20] G.A. Leonov, Limit cycles of the Lienard equation with discontinuous coeffi- cients, Dokl. Akad. Nauk 426 (1) (2009) 47–50. [21] L. Sheng, Limit cycles of a class of piecewise smooth Liénard systems, Internat. J. Bifur. Chaos Appl. Sci. Engrg. 26 (1) (2016) 1650009 10. [22] P.T. Cardin, J. Torregrosa, Limit cycles in planar piecewise linear differential systems with nonregular separation line, Physica D 337 (2016) 67–82. [23] A. Gasull, J. Torregrosa, Center-focus problem for discontinuous planar differ- ential equations, Internat. J. Bifur. Chaos Appl. Sci. Engrg. 13 (7) (2003) 1755– 1765. Dynamical systems and functional equations (Murcia, 2000). [24] J.-P. Françoise, Successive derivatives of a first return map, application to the study of quadratic vector fields, Ergodic Theory Dynam. Systems 16 (1) (1996) 87–96. [25] J.-P. Françoise, The first derivative of the period function of a plane vector field, in: Proceedings of the Symposiumon Planar Vector Fields, Lleida, 1996, vol. 41, 1997, pp. 127–134. [26] J.-P. Françoise, The successive derivatives of the period function of a plane vector field, J. Differential Equations 146 (2) (1998) 320–335. [27] M. Abramowitz, I.A. Stegun, HandBook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, in: National Bureau of Stan- dards Applied Mathematics Series, vol. 55, U.S. Government Printing Office, Washington D.C., 1964. For sale by the Superintendent of Documents. http://refhub.elsevier.com/S0167-2789(17)30465-7/sb1 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb1 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb1 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb1 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb1 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb2 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb2 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb2 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb2 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb2 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb3 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb3 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb3 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb3 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb3 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb4 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb4 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb4 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb4 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb4 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb5 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb5 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb5 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb6 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb6 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb6 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb6 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb6 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb7 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb7 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb7 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb8 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb8 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb8 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb9 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb9 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb9 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb9 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb9 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb10 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb10 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb10 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb11 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb11 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb11 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb13 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb13 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb13 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb14 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb14 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb14 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb14 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb14 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb14 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb14 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb15 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb15 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb15 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb16 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb16 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb16 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb17 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb17 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb17 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb17 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb17 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb18 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb18 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb18 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb19 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb19 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb19 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb19 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb19 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb19 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb19 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb20 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb20 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb20 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb21 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb21 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb21 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb22 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb22 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb22 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb23 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb23 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb23 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb23 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb23 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb24 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb24 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb24 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb24 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb24 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb26 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb26 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb26 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb27 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb27 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb27 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb27 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb27 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb27 http://refhub.elsevier.com/S0167-2789(17)30465-7/sb27 Limit cycles via higher order perturbations for some piecewise differential systems Introduction The return map and the Poincare–Pontryagin–Melnikov functions Polynomial perturbation Perturbations of order one and two for some Lienard families Perturbations of any order for some classic Lienard families First order perturbations of Lienard families of low degree Acknowledgments References