This is a preprint of: “Melnikov analysis in nonsmooth differential systems with nonlinear switch- ing manifold”, Jéfferson L. R. Bastos, Claudio A. Buzzi, Jaume Llibre, Douglas D. Novaes, J. Differential Equations, vol. 267, 3748–3767, 2019. DOI: [10.1016/j.jde.2019.04.019] MELNIKOV ANALYSIS IN NONSMOOTH DIFFERENTIAL SYSTEMS WITH NONLINEAR SWITCHING MANIFOLD Jéfferson L. R. Bastos, Claudio A. Buzzi Universidade Estadual Paulista, IBILCE-UNESP Av. Cristovão Colombo, 2265, 15.054-000, S. J. Rio Preto, SP, Brasil jeferson@ibilce.unesp.br and buzzi@ibilce.unesp.br Jaume Llibre Universitat Autònoma de Barcelona, UAB Edifici C Facultat de Ciències, 08193 Bellaterra, Barcelona, Spain jllibre@mat.uab.cat Douglas D. Novaes Universidade Estadual de Campinas, IMECC-UNICAMP R. Sérgio Buarque de Holanda, 651, 13.083-859, Campinas, SP, Brasil ddnovaes@ime.unicamp.br Abstract. We study the family of piecewise linear differential systems in the plane with two pieces separated by a cubic curve. Our main result is that 7 is a lower bound for the Hilbert num- ber of this family. In order to get our main result, we develop the Melnikov functions for a class of nonsmooth differential sys- tems, which generalizes, up to order 2, some previous results in the literature. Whereas the first order Melnikov function for the nonsmooth case remains the same as for the smooth one (i.e. the first order averaged function) the second order Melnikov function for the nonsmooth case is different from the smooth one (i.e. the second order averaged function). We show that, in this case, a new term depending on the jump of discontinuity and on the geometry of the switching manifold is added to the second order averaged function. 2010 Mathematics Subject Classification. Primary 34C07, 37G15, 34A36. Key words and phrases. Melnikov theory, averaging theory, nonsmooth differ- ential systems, piecewise linear differential systems, nonlinear switching manifold, limit cycles, Hilbert number. 1 2 J. BASTOS, C. BUZZI, J. LLIBRE, D. D. NOVAES 1. Introduction and statement of the main results In recent years, there has been a growing interest in studying non- smooth differential systems, motivated mainly by their engineering ap- plications. In particular, piecewise linear differential systems have been used to model many real processes and different modern devices, for more details see, for instance, [7] and the references therein. In the smooth differential systems many results have been derived under con- venient smoothness assumptions. Thus, the following natural questions arise: What results from the smooth differential systems theory extend to the nonsmooth ones? How we must modify the results from the smooth differential systems theory in order that they work for the non- smooth ones? These questions are not merely academic because the nonsmooth differential systems appear in a natural way in the context of many applications, see for instance [5, 27, 30, 31]. Here, we are mainly interested in studying the existence of limit cy- cles for piecewise linear differential systems with two pieces separated by a nonlinear switching curve. The Averaging Theory is a classical method to attack this problem. It has been developed also for non- smooth differential systems (see [15, 18, 19, 23]). However, in the pre- vious works it is assumed some strong conditions on the switching set, which do not hold when it is a cubic, for instance. So here, we also de- velop the bifurcation functions of first and second order for computing periodic solutions of a wider class of nonsmooth differential systems. Usually, these bifurcation functions are called Melnikov Functions. Ad- ditionally, we compare them with the Melnikov functions of first and second order for smooth systems. Usually, piecewise linear differential systems have been considered when a straight line separates the plane in two half–planes. In each of the half–planes we have a linear differential system. If both linear systems coincide in the switching line we call it a continuous piecewise linear differential system. Otherwise we call it a discontinuous piece- wise linear differential system. In recent years, many authors have studied intensively the number of limit cycles of discontinuous piece- wise linear differential systems with two zones separated by a straight line, see for instance [1, 4, 6, 8, 9, 10, 12, 13, 14, 17, 21, 22, 24, 25, 26] and the references quoted in these papers. Up to now all results of these papers provide examples of at most 3 crossing limit cycles for this class of discontinuous systems. It remains the open question: Is 3 the maximum number of crossing limit cycles that such discontinuous differential systems can exhibit? MELNIKOV ANALYSIS IN NONSMOOTH DIFFERENTIAL SYSTEMS 3 In this paper, we consider the family Lh of piecewise linear differ- ential systems in the plane with two pieces separated by the curve Σ = h−1 (0), h(x, y) = y − x3 . More precisely, we shall study the class of discontinuous piecewise linear differential systems obtained by per- turbing up to order N ∈ N in the small parameter ε the linear center ẋ = y, ẏ = −x, i.e.   PN   y+ εi Pi+ (x, y), y ≥ x3 i=1 ẋ =   PN  y+ εi Pi− (x, y), y ≤ x3 i=1 (1)   P N   −x + εi Q+ i (x, y), y ≥ x3 i=1 ẏ =   P N  −x + εi Q− i (x, y), y ≤ x3 i=1 where the functions Pi± (x, y) and Q± i (x, y) are the following polynomi- als of degree one Pi+ (x, y) = a0i + a1i x + a2i y, Pi− (x, y) = b0i + b1i x + b2i y, Q+i (x, y) = α0i + α1i x + α2i y, Q−i (x, y) = β0i + β1i x + β2i y. We shall assume that N = 1, 2. As usual, we denote by HLh the maximum number of limit cycles that piecewise linear differential systems in Lh can have. This number is called Hilbert number. Previous results show that HLA ≥ 3, when A : R2 → R is a linear function. Our main result provides a lower bound for HLh , when h(x, y) = y − x3 . Theorem A. For |ε| = 6 0 sufficiently small, there exist piecewise linear differential systems of kind (1) admitting at least 7 limit cycles, which bifurcate from the periodic orbits of the linear center. Consequently, for h(x, y) = y − x3 we have HLh ≥ 7. Theorem A is proved in Section 3. Its proof is mainly based on the Chebyshev Theory, recently extended for Chebyshev systems with positive accuracy (see [29]), and Melnikov Theory. In order to prove Theorem A we develop the Melnikov functions for a class of nonsmooth differential systems (see Theorem B in Section 2), which generalizes, up to order 2, the results of [15, 19]. 4 J. BASTOS, C. BUZZI, J. LLIBRE, D. D. NOVAES It is important to mention that at this moment in the literature the maximum number of limit cycles exhibited by piecewise linear differen- tial systems in two zones separated by a straight line is three. Theorem A shows that the nonlinearity of the switching curve Σ is responsible for the increasing of the number of limit cycles. In fact, this phenom- ena has already been observed in [2, 28]. Indeed, in [28] it was proved that the number of limit cycles can increase arbitrarily with the num- ber of oscillations of the switching curve Σ. By oscillations of Σ, one can think as the transversal intersections of Σ with the straight line x = 0. For instance, in [28] the authors considered the switching curve Σ = {(0, y) if y ≤ 0} ∪ {(h(y), y) if y > 0}, where   k sin(πy), 0 ≤ y ≤ 2n+1 2 , h(y) =  (−1)n k, y > 2n+1 2 . It was shown that the number of limit cycles of a particular piecewise linear differential system corresponds to the number of zeros of the function h, which corresponds to the intersections between Σ and the straight line x = 0. 2. Melnikov functions for a class of nonsmooth systems In order to state our results for the nonsmooth differential systems, we need some notations and definitions. Let D be an open set of Rd , S1 = R/T for some period T > 0, and N a positive integer. We consider a finite set of C r (r ≥ N + 1) functions θi : D → S1 for i = 1, 2, . . . , M, satisfying 0 < θ1 (x) < θ2 (x) < · · · < θM (x) < T for all x ∈ D. For completeness, we take θ0 (x) ≡ 0 and θM +1 (x) ≡ T. We consider the following nonsmooth differential system N X (2) ẋ = εi Fi (t, x) + εN +1 R(t, x, ε), i=1 with    F 0 (t, x), 0 < t < θ1 (x),  i1 Fi (t, x), θ1 (x) < t < θ2 (x), Fi (t, x) =   ...  F M (t, x), θ (x) < t < T, i M and    R0 (t, x, ε), 0 < t < θ1 (x),  1 R (t, x, ε), θ1 (x) < t < θ2 (x), R(t, x, ε) =   ...  RM (t, x, ε), θ (x) < t < T, M MELNIKOV ANALYSIS IN NONSMOOTH DIFFERENTIAL SYSTEMS 5 where Fij : S1 × D → Rd , Rj : S1 × D × (−ε0 , ε0 ) → Rd , for i = 1, 2, . . . , N and j = 0, 1, . . . , M, are C r functions, r ≥ N + 1, and T -periodic in the variable t. Notice that, the switching manifold is given by Σ = {(θi (x), x) : x ∈ D, i = 0, 1, . . . , M, M + 1}. Through- PN i j out the paper we shall also denote F j (t, x, ε) = i=1 ε Fi (t, x) + N +1 j ε R (t, x, ε). The Melnikov Theory and the Averaging Theory are classical tools to investigate the existence of periodic solutions of perturbative systems of kind (2). In short, both theories provide a sequence of functions ∆i , i = 1, 2, . . . , k, which “control” the existence of isolated periodic solu- tions of (2). These functions are called Melnikov Functions or Bifurca- tion Functions. In [19], the averaging theory at any order was developed for systems of kind (2) assuming that the functions θi , i = 1, 2, . . . , M, were constant. It was shown that in this case the bifurcation func- tions coincide with the averaged functions, denoted by fi ’s, for smooth systems (see [20]). In particular, for i = 1, 2, we have Z T f1 (x) = F1 (s, x)ds and 0 (3) Z  Z  T s f2 (x) = Dx F1 (s, x) F1 (t, x)dt + F2 (s, x) ds. 0 0 Nevertheless, in [18] it was observed that higher order averaged func- tions, fi , i ≥ 2, do not always control the bifurcation of isolated periodic solutions for nonsmooth differential systems, and a strong degenerate condition on the switching manifold was stablished in order that f2 stands as the bifurcation function of order 2 (see hypotheses (Hb2) of [18, Theorem B]). As the main result of this section, Theorem B shows that the Mel- nikov functions of first and second order, ∆1 and ∆2 , of system (2) write (4) ∆1 (x) = f1 (x) and ∆2 (x) = f2 (x) + f2∗ (x), where f1 and f2 are given by (3), and M  X  Z θj (x) j−1 j f2∗ (x) = F1 (θj (x), x) − F1 (θj (x), x) Dx θj (x) F1 (s, x)ds. j=1 0 This generalizes up to order 2 the results of [15, 19] for a wider class of nonsmooth differential equations. In what follows, J∆i (x) denotes the corresponding Jacobian matrix of ∆i evaluated at x. 6 J. BASTOS, C. BUZZI, J. LLIBRE, D. D. NOVAES Theorem B. Consider the nonsmooth differential system (2) and the functions ∆1 and ∆2 defined in (4). So, the following statements hold. i. (First Order) Assume that a∗ ∈ D satisfies ∆1 (a∗ ) = 0 and det (J∆1 (a∗ )) 6= 0. Then, for |ε| = 6 0 sufficiently small, there exists a unique T -periodic solution x(t, ε) of system (2) such that x(0, ε) → a∗ as ε → 0. ii. (Second Order) Assume that ∆1 (x) ≡ 0 and that a∗ ∈ D satifies ∆2 (a∗ ) = 0 and det (J∆2 (a∗ )) 6= 0. Then, for |ε| 6= 0 sufficiently small, there exists a unique T -periodic solution x(t, ε) of system (2) such that x(0, ε) → a∗ as ε → 0. From statement (i) of Theorem B, we see that the first order Mel- nikov function ∆1 (x) for nonsmooth differential systems coincides with the corresponding one f1 (x) for smooth differential systems. This had already been observed in [18]. The second order Melnikov function for smooth differential systems is equal to f2 (x). Indeed, for smooth differ- ential systems the Melnikov functions are, in some sense, equivalent to the averaged functions at any order (for more details, see [11]). Nev- ertheless, from statement (ii) of Theorem B, the corresponding second order Melnikov function for nonsmooth differential systems of kind (2) has the extra term f2∗ (x), which depends on the jump of discontinuity and on the geometry of the switching manifold measured, respectively, by the expressions Z θj (x) j−1 j F1 (θj (x), x) − F1 (θj (x), x) and Dx θj (x) F1 (s, x)ds. 0 We notice that, in particular, if either the vector field F1 is con- tinuous, or if the the switching manifold Σ is “vertical”, i.e. Dx θj (x) vanishes identically for j = 1, 2, . . . , M, then the extra term f2∗ van- ishes. Under one of these two last assumptions Theorem B provides the known results given in [3, 20] for smooth systems, and in [15, 18, 19, 23] for nonsmooth differential systems. In what follows, before the proof of Theorem B, we provide a general geometric interpretation for the increment f2∗ . In short, the next result states that f2∗ 6= 0 stands as a kind of transversal condition between the switching manifold Σ and the average of the perturbation F1 , which is highly related to hypotheses (Hb2) of [18, Theorem B]. Accordingly, given p = (tp , xp ) ∈ Σ, we define the vector  Z  1 tp v(p) = 0, F1 (t, xp )dt . tp 0 MELNIKOV ANALYSIS IN NONSMOOTH DIFFERENTIAL SYSTEMS 7 Notice that the second component of v(p) is the average of the first order perturbation F1 of (2) computed on the backward orbit {(t, xp ) : t ∈ [0, tp ]} of the unperturbed differential equation (2) ε=0 with initial condition (t(0), x(0)) = p. Proposition 1. If f2∗ 6= 0, then there exists p ∈ Σ such that v(p) is transversal to Σ at p. Proof. We know that, for p = (tp , xp ) ∈ Σ, tp = θj (xp ) for some j ∈ {0, 1, . . . , M, M + 1}. Notice that Z θj (xp ) Dx θj (xp ) F1 (s, x)ds = 0 0 is equivalent to v(p) ∈ Tp Σ. Consequently, v(p) ∈ Tp Σ for every p ∈ Σ implies that f2∗ = 0.  Proof of Theorem B. We denote the solution of (2) by    ϕ (t, x, ε), 0 ≤ t ≤ α1 (x, ε);  0 ϕ1 (t, x, ε), α1 (x, ε) ≤ t ≤ α2 (x, ε); (5) ϕ(t, x, ε) =   ···  ϕ (t, x, ε), α (x, ε) ≤ t ≤ T. M M In the above expression, αj (x, ε) is the flying time that the trajectory ϕj−1(·, x, ε), starting at ϕj−1(αj−1 (x, ε), x, ε) ∈ D for t = αj−1 (x, ε), reaches the manifold {(θj (x), x) : x ∈ D} ⊂ Σ (see Figure 1), that is, αj (x, ε) = θj (ϕj−1(αj (x, ε), x, ε)). for j = 1, 2, . . . , M, α0 (x, ε) = 0, and αM +1 (x, ε) = T. Moreover, ∂ϕj (t, x, ε) = F j (t, ϕj (t, x, ε), ε), for j = 0, 1, . . . , M, where ∂t (6)  ϕ0 (0, x, ε) = x, ϕj (αj (x, ε), x, ε) = ϕj−1(αj (x, ε), x, ε), for j = 1, 2, . . . , M. From (6) and from the differential dependence of the solutions on the initial conditions and on the parameter, we can see inductively that αj (x, ε) and ϕj (t, x, ε) are C r functions, for j = 0, 1, . . . , M. We notice that the recurrence (6) describes initial value problems. Therefore, it is equivalent to the following recurrence: Z t ϕ0 (t, x, ε) = x + F 0 (s, ϕ0 (s, x, ε), ε)ds, 0 Z t ϕj (t, x, ε) = ϕj−1(αj (x, ε), x, ε) + F j (s, ϕj (s, x, ε), ε)ds, αj (x,ε) 8 J. BASTOS, C. BUZZI, J. LLIBRE, D. D. NOVAES x θ1 (x) · · · θj−1 (x) θj (x) · · · θM (x) ϕ(T, x0 , ε) ϕj−1 (t, x0 , ε) ϕ0 (t, x0 , ε) x0 ϕM (t, x0 , ε) t 0 αj−1(x0 , ε) αj (x0 , ε) T Figure 1. Illustration of a solution (5) of system (2) starting at x0 and crossing the switching manifold for j = 1, 2, . . . , M. So, we may compute XM Z αj (x,ε) ϕM (t, x, ε) = x + F j−1(s, ϕj−1(s, x, ε), ε)ds j=1 αj−1 (x,ε) Z t + F M (s, ϕM (s, x, ε), ε)ds. αM (x,ε) The displacement function is given by (7) ∆(x, ε) = ϕ(T, x, ε) − x = ϕM (T, x0 , ε) − x. Notice that, from the above comments, ∆(x, ε) is a C r function. Then, expanding (7) around ε = 0 we have ∆(x, ε) = ∆0 (x) + ε∆1 (x) + ε2 ∆2 (x) + O(ε3 ). We observe that ∆0 (x) ≡ 0, because ∆0 (x) = ∆(x, 0) = ϕ(T, x, 0) − x. The functions ∆1 (x) and ∆2 (x) are called Melnikov functions of first and second order, respectively. The displacement function reads M X +1 Z αj (x,ε) (8) ∆(x, ε) = ϕM (T, x, ε) − x = F j−1(s, ϕj−1(s, x, ε), ε)ds. j=1 αj−1 (x,ε) MELNIKOV ANALYSIS IN NONSMOOTH DIFFERENTIAL SYSTEMS 9 In what follows, we compute the expansion, around ε = 0, of the M + 1 summands of (8). For the first one we obtain Z α1 (x,ε) F 0 (s, ϕ0 (s, x, ε), ε)ds 0 Z α1 (x,ε) = [εF10(s, ϕ0 (s, x, ε)) + ε2 F20 (s, ϕ0 (s, x, ε)) + O(ε3 )]ds 0 Z θ1 (x) ! Z θ1 (x)  (9) 0 2 ∂ϕ0 =ε F1 (s, x)ds + ε DF10 (s, x) (s, x, 0) 0 0 ∂ε  ! ∂α 1 + F20 (s, x) ds + F10 (θ1 (x), x) (x, 0) + O(ε3 ). ∂ε In order to finish the computation of the second summand we have to compute the terms ∂ϕ ∂ε 0 (s, x, 0) and ∂α ∂ε 1 (x, 0). Expanding the equation Z t ϕ0 (t, x, ε) = x + F 0 (s, ϕ0 (s, x, ε), ε)ds 0 around ε = 0 we obtain Z t ∂ϕ0 x+ε (t, x, 0) + O(ε2 ) = x + ε F10 (s, x)ds + O(ε2 ). ∂ε 0 So, we get Z t ∂ϕ0 (10) (t, x, 0) = F10 (s, x)ds. ∂ε 0 In the same way, expanding the equation α1 (x, ε) = θ1 (ϕ0 (α1 (x, ε), x, ε)) around ε = 0 we obtain ∂α1 α1 (x, 0) + ε (x, 0) + O(ε2 ) = ∂ε   ∂ϕ0 ∂α1 ∂ϕ0 θ1 (x) + εDx θ1 (x) (θ1 (x), x, 0) (x, 0) + (θ1 (x), x, 0) + O(ε2 ). ∂t ∂ε ∂ε Since ϕ0 (t, x, 0) = x for all t we get Z θ1 (x) ∂α1 (11) (x, 0) = Dx θ1 (x) F10 (s, x)ds. ∂ε 0 10 J. BASTOS, C. BUZZI, J. LLIBRE, D. D. NOVAES Substituting (10) and (11) in (9) we have the expression for the first summand of (8) Z Z ! α1 (x,ε) θ1 (x) F 0 (s, ϕ0 (s, x, ε), ε)ds = ε F10 (s, x)ds 0 0 Z θ1 (x)  Z s  2 (12) +ε Dx F10 (s, x) F10 (t, x)dt + F20 (s, x) ds 0 0 Z ! θ1 (x) + F10 (θ1 (x), x)Dx θ1 (x) F10 (s, x)ds + O(ε3 ). 0 The other summands of (8) can be computed in a similar way and they are given by Z αj (x,ε) Z θj (x) ! F j−1(s, ϕj−1(s, x, ε), ε)ds = ε F1j−1 (s, x)ds αj−1 (x,ε) θj−1 (x) Z θj (x)  Z s  +ε 2 Dx F1j−1(s, x) F1 (t, x)dt + F2j−1 (s, x) ds θj−1 (x) 0 (13) Z θj (x) + F1j−1(θj (x), x)Dx θj (x) F1 (s, x)ds 0 Z ! θj−1 (x) − F1j−1(θj−1 (x), x)Dx θj−1 (x) F1 (s, x)ds + O(ε3 ). 0 From (12) and (13), we have that the Melnikov functions of first and second order are given by Z T ∆1 (x) = F1 (s, x)ds 0 and ∆2 (x) = f2 (x) + f2∗ (x), where Z T  Z s  f2 (x) = Dx F1 (s, x) F1 (t, x)dt + F2 (s, x) ds 0 0 and M  X  Z θj (x) j−1 j f2∗ (x) = F1 (θj (x), x) − F1 (θj (x), x) Dx θj (x) F1 (s, x)ds. j=1 0 MELNIKOV ANALYSIS IN NONSMOOTH DIFFERENTIAL SYSTEMS 11 From the definition of the displacement function in (7), it is clear that the T -periodic solutions x(t, ε) of system (2), satisfying x(0, ε) = x, are in one-to-one correspondence to the zeros of the equation ∆(x, ε) = 0. b Now, define ∆(x, ε) = ∆(x, ε)/ε = ∆1 (x) + O(ε). From hypotheses, we have that ! ∂ b ∆ b ∗ , 0) = ∆1 (a∗ ) = 0, and det ∆(a (a∗ , 0) = det (J∆1 (a∗)) 6= 0. ∂x Therefore, from the Implicit Function Theorem, we get the existence of a unique C r function a(ε) ∈ D such that a(0) = a∗ and ∆(a(ε), ε) = b ∆(a(ε), ε) = 0, for every |ε| = 6 0 sufficiently small. This completes the proof of statement (i). e Finally, assuming that ∆1 (x) ≡ 0, define ∆(x, ε) = ∆(x, ε)/ε2 = ∆2 (x) + O(ε). So, the proof of statement (ii) follows analogously to the proof of statement (i).  3. Bifurcation of limit cycles In this section we shall apply the Melnikov functions developed in the previous section to study the bifurcation of limit cycles of system (1) for |ε| = 6 0 sufficiently small. Then, Theorem A will be a direct consequence of Proposition 5. A useful tool to study the number of isolated zeros of a function is the Chebyshev Theory. We recall that the Wronskian of the ordered k + 1 functions u0 , . . . , uk is  (14) Wk (x) = Wk (u0 , . . . , uk )(x) = det M(u0 , . . . , uk )(x) , where   u0 (x) ··· uk (x)  u′0 (x) ··· u′k (x)    M(u0 , . . . , uk )(x) =  .. .. .. .  . . .  (k) (k) u0 (x) · · · uk (x) We say that F = [u0 , u1, . . . , un ] is an Extended Complete Chebyshev system or an ECT-system on a closed interval [a, b] if and only if W (u0, u1 , . . . , uk )(x) 6= 0 on [a, b] for 0 ≤ k ≤ n, see [16]. In what follows, Span(F ) denotes the set of all functions given by linear combi- nations of the functions of F . A classical result concerning ECT-system is the following: Theorem 2 ([16]). Let F = [u0 , u1 , . . . , un ] be an ECT-system on a closed interval [a, b]. Then, the number of isolated zeros for every 12 J. BASTOS, C. BUZZI, J. LLIBRE, D. D. NOVAES element of Span(F ) does not exceed n. Moreover, for each configuration of m ≤ n zeros, taking into account their multiplicity, there exists F ∈ Span(F ) with this configuration of zeros. In [29], the above result was extended for Chebyshev systems with positive accuracy. In what follows, we say that a zero of a function is simple if its Jacobian matrix evaluated at the zero is a nonsingular matrix. Theorem 3 ([29]). Let F = [u0 , u1 , . . . , un ] be an ordered set of func- tions on [a, b]. Assume that all the Wronskians Wk (x), k = 0, . . . , n−1, are nonvanishing except Wn (x), which has exactly one zero on (a, b) and this zero is simple. Then, the number of isolated zeros for every ele- ment of Span(F ) does not exceed n+1. Moreover, for any configuration of m ≤ n + 1 zeros there exists F ∈ Span(F ) realizing it. As a first application of Theorem B we consider piecewise linear differential systems (1) with N = 1. The next result shows that in this case the first order Melnikov function ∆1 has at most 3 isolated zeros. From statement (i) of Theorem B, this means that the differential system (1) has at most 3 limit cycles bifurcating from the periodic orbits of the linear center. Moreover, we shall see that this maximum number is reached by particular examples. Proposition 4. Consider the piecewise linear differential system (1) with N = 1, and let ∆1 be its first order Melnikov function (4). Then, the number of isolated zeros of ∆1 does not exceed 3. Moreover, there exist linear differential systems (1) such that ∆1 has exactly 3 simple zeros. Proof. In order to apply the Melnikov method, we must write the dif- ferential system (1) in polar coordinates. Performing the change of variables (x, y) = (r cos θ, r sin θ) and taking θ as the new independent variable, the differential system (1) is transformed into the differential equation dr = εF1± (θ, r) + ε2 F2± (θ, r) + O(ε3 ), dθ where  F1+ (θ, r) = −(a01 cos θ + b01 sin θ) − 2r a11 + b21 −  r 2 (a11 − b21 ) cos(2θ) + (a21 + b11 ) sin(2θ) and  F1− (θ, r) = −(α01 cos θ + β01 sin θ) − 2r α11 + β21 −  r 2 (α11 − β21 ) cos(2θ) + (α21 + β11 ) sin(2θ) . MELNIKOV ANALYSIS IN NONSMOOTH DIFFERENTIAL SYSTEMS 13 The switching curve Σ is given in polar coordinates by sin θ = r 2 cos3 θ. It means that for each r > 0 there exist two angles θ1 (r) ∈ [0, π2 ] and θ2 (r) = θ1 (r)+π ∈ [π, 3π 2 ] such that the circle x2 +y 2 = r 2 in- tersects Σ at the points (r cos θ1 (r), r sin θ1 (r)) and (r cos θ2 (r), r sin θ2 (r)) (see Figure 2). Σ r x3 θ1 x Figure 2. Identities involving x, r and θ1 . According to the theory developed in Section 2, the Melnikov func- tion of first order is given by Z θ1 (r) Z θ1 (r)+π Z 2π ! 1 − + − ∆1 (r) = F1 (θ, r)dθ + F1 (θ, r)dθ + F1 (θ, r)dθ 2π 0 θ1 (r) θ1 (r)+π = γ0 cos θ1 (r) + γ1 r + γ2 sin θ1 (r), where γ0 = −2(b01 − β01 ), (15) γ1 = −π 2 (a11 + b21 + α11 + β21 ), γ2 = a01 − α01 . From Theorem B, each isolated zero of ∆1 (r) = 0 corresponds to a limit cycle of (1). We prove now that the maximum number of positive zeros ∆1 (r) = 0 is three. As we can see in Figure 2, the variables x, r and θ1 are related by the identities x x3 (16) r 2 = x2 + x6 , cos θ1 = , and sin θ1 = . r r 14 J. BASTOS, C. BUZZI, J. LLIBRE, D. D. NOVAES For r > 0 we have that ∆1 (r) = 0 if and only if r∆1 (r) = 0. So, the maximum number of limit cycles of (1) obtained by the first or- der Melnikov method is the maximum number of positive zeros of the polynomial (17) p1 (x) = γ0 + γ1 (x + x5 ) + γ2 x2 . Consider the ordered set of functions F = [u0 , u1 , u2 ], where u0 (x) = 1, u1 (x) = x + x5 and u2 (x) = x2 . We have that the Wronskians, defined in (14), are W0 (x) = 1, W1 (x) = 1 + 5x4 and W2 (x) = 2 − 30x4 . The Wronskians W0 (x) and W1 (x) do not vanish in R and W2 (x) has exactly one positive zero. So, according to Theorem 3, three is the maximum number of positive roots of polynomial (17) and there is a choice for γ0 , γ1 and γ2 such that p1 (x) has exactly three zeros. From (15), it is clear that for any choice of γ’s it is possible to obtain the corresponding perturbation. This concludes the proof of Proposition 4.  Example 1. In what follows, we exhibit a concrete perturbation pre- senting three limit cycles for |ε| = 6 0 sufficiently small. For γ0 = −0.05, γ1 = 1 and γ2 = −2, the polynomial (17) has three positive roots (see Figure 3). It is clear that the linear system 0.2 0.1 -0.2 0.2 0.4 0.6 0.8 1.0 -0.1 -0.2 -0.3 Figure 3. Graph of p1 (x) = x5 − 2x2 + x − 0.05. −2(b01 − β01 ) = −0.05, −π 2 (a11 + b21 + α11 + β21 ) = 1, a01 − α01 = −2, has infinity many solutions. For example, we can take a01 = −2, a11 = −2/π, b01 = 1/40 and all other coefficients equal to zero. So, for |ε| = 6 0 MELNIKOV ANALYSIS IN NONSMOOTH DIFFERENTIAL SYSTEMS 15 sufficiently small, system  y − 2ε(1 + π1 x), if y ≥ x3 , ẋ = 1 y + 40 ε, if y ≤ x3 , ẏ = −x has three limit cycles. As a second application of Theorem B, we consider piecewise linear differential systems (1) with N = 2. We impose conditions on the coefficients of P1± and Q± 1 in order that the first order Melnikov function ∆1 vanishes identically. The next result shows that in this case the second order Melnikov function ∆2 has at most 7 isolated zeros. From statement (ii) of Theorem B, this means that the differential system (1) has at most 7 limit cycles bifurcating from the periodic orbits of the linear center. Moreover, we shall see that this maximum number is reached by particular examples. This proves Theorem A. Proposition 5. Consider the piecewise linear differential system (1) with N = 2, and let ∆1 and ∆2 be its first and second order Melnikov functions (4), respectively. In addition, assume that ∆1 = 0. Then, the number of isolated zeros of ∆2 does not exceed 7. Moreover, there exist linear differential systems (1) such that ∆2 has exactly 7 simple zeros. Proof. In order to apply the Melnikov method of second order it is necessary that the Melnikov function of first order ∆1 be identically zero. So, we impose that the γ’s defined in (15) are identically zero. It is enough to take a11 = −(b21 + α11 + β21 ), (18) b01 = β01 , a01 = α01 . Again, we use polar coordinates (x, y) = (r cos θ, r sin θ) and take θ as the new independent variable to transform system (1) into the differ- ential equation dr = εF1± (θ, r) + ε2 F2± (θ, r) + ε3 R2± (θ, r, ε), dθ where F1± (θ, r) = f11 ± ± (θ) + rf12 (θ), 1 ± F2± (θ, r) = f21 ± (θ) + f22 ± (θ) + rf23 (θ), r 16 J. BASTOS, C. BUZZI, J. LLIBRE, D. D. NOVAES with ± f11 (θ) = c± ± 101 cos θ + s101 sin θ, ± f12 (θ) = c± ± ± 110 + c112 cos(2θ) + s112 sin(2θ), ± f21 (θ) = c± ± 202 cos(2θ) + s202 sin(2θ), ± f22 (θ) = c± ± ± ± 211 cos θ + s211 sin θ + c213 cos(3θ) + s213 sin(3θ), ± f23 (θ) = c± ± ± ± ± 220 + c222 cos(2θ) + c224 cos(4θ) + s222 sin(2θ) + s224 sin(4θ), and 1 c+ 101 = −α01 , s+ 101 = −β01 , c+ 110 = (α11 + β21 ), 2 1 1 c+ 112 = (2b21 + α11 + β21 ), s+112 = (−a21 − b11 ), 2 2 1  c+ 202 = −α01 β01 , s+ 202 = α2 − β012 , 2 01 1 c+ 211 = (−2a02 + a21 α01 − b11 α01 + α11 β01 + β01 β21 ), 2 1 s+ 211 = (a21 β01 − 2b02 − b11 β01 − α01 α11 − α01 β21 ), 2 1 c+ 213 = (−a21 α01 − b11 α01 + 2b21 β01 + α11 β01 + β01 β21 ), 2 1 s+ 213 = (−a21 β01 − b11 β01 − 2b21 α01 − α01 α11 − α01 β21 ), 2 1 c+ 220 = (−2a12 − a21 α11 − a21 β21 + b11 α11 + b11 β21 − 2b22 ), 4 1 c+ 222 = (−a12 − a21 b21 + b11 b21 + b11 α11 + b11 β21 + b22 ), 2 1 2  s+ 222 = a21 − 2a22 − b211 − 2b12 + 2b21 α11 + 2b21 β21 + (α11 + β21 )2 , 4 1 c+ 224 = (a21 + b11 )(2b21 + α11 + β21 ), 4 1  s+ 224 = −a221 − 2a21 b11 − b211 + (2b21 + α11 + β21 )2 , 8 1 c− 101 = −α01 , s− 101 = −β01 , c− 110 = (−α11 − β21 ), 2 1 1 c− 112 = (β21 − α11 ), s− 112 = (2b21 + α11 + β21 ), 2 2 1  c− 202 = −α01 β01 , s−202 = 2 α01 2 − β01 , 2 1 c− 211 = (α01 α21 − α01 β11 − 2α02 − α11 β01 − β01 β21 ), 2 MELNIKOV ANALYSIS IN NONSMOOTH DIFFERENTIAL SYSTEMS 17 1 s− 211 = (α01 (α11 + β21 ) + α21 β01 − β01 β11 − 2β02 ), 2 1 c− 213 = (β01 (β21 − α11 ) − α01 (α21 + β11 )), 2 1 s− 213 = (α01 (α11 − β21 ) − β01 (α21 + β11 )), 2 1 c− 220 = (α11 α21 − α11 β11 − 2α12 + α21 β21 − β11 β21 − 2β22 ), 4 1 c− 222 = (−α11 β11 − α12 − α21 β21 + β22 ), 2 1 2  s− 222 = 2 α11 + α21 − 2α22 − β11 2 − 2β12 − β21 2 , 4 1 c− 224 = − (α11 − β21 )(α21 + β11 ), 4 1  s− 224 = 2 α11 − 2α11 β21 − α21 2 − 2α21 β11 − β11 2 2 + β21 . 8 According to Section 2 the Melnikov function of second order ∆2 is given by ∆2 (r) = f2 (r) + f2∗ (r), where Z θ1 (r) Z θ  1 − − f2 (r) = Dr F1 (θ, r) F1 (t, r)dt + F2 (θ, r) dθ 2π 0 0 Z θ2 (r) Z θ  1 + + + Dr F1 (θ, r) F1 (t, r)dt + F2 (θ, r) dθ 2π θ1 (r) 0 Z 2π  Z θ  1 − − + Dr F1 (θ, r) F1 (t, r)dt + F2 (θ, r) dθ. 2π θ2 (r) 0 and   Z θ1 (r) f2∗ (r) = F1− (θ1 (r), r) − F1+ (θ1 (r), r) ′ θ1 (r) F1 (s, r)ds 0   Z θ2 (r) + F1+ (θ2 (r), r) − F1− (θ2 (r), r) θ2′ (r) F1 (s, r)ds. 0 We recall that in our context θ2 (r) = θ1 (r) + π. So, we have θ1′ (r) = θ2′ (r). For simplicity in the following formulas we denote θ1 (r) by θ1 . After some computations we obtain f2 (r) = δ1 r + δ2 cos θ1 + δ3 sin θ1 + δ4 cos3 θ1 + δ5 sin3 θ1 , and r 2 cos2 θ1   f2∗ (r) = µ 1 r + µ 2 cos θ1 + µ 3 sin θ1 h(θ1 ), 2 + 6r 2 cos θ1 sin θ1 18 J. BASTOS, C. BUZZI, J. LLIBRE, D. D. NOVAES where h(θ1 ) = η1 + η2 cos2 θ1 + η3 cos θ1 sin θ1 , π δ1 = − 2a12 + 2b22 + a21 α11 − b11 α11 + 2α12 − α11 α21 4  + α11 β11 + a21 β21 − b11 β21 − α21 β21 + β11 β21 + 2β22 ,  δ2 = 2 − b02 + 2b21 α01 + a21 β01 − α21 β01 + β02 − 2α01 β21 ,  δ3 = 2 a02 + b11 α01 − α02 − 2b21 β01 − 4α11 β01 − α01 β11 − 2β01 β21 ,   δ4 = α01 − 4b21 − 4α11 + β01 − 2a21 − 2b11 + 2α21 + 2β11 ,   δ5 = 4b21 + 4α11 β01 + α01 − 2a21 − 2b11 + 2α21 + 2β11 ,  µ1 = π(α11 + β21 ), η1 = 2 β21 − b21 ,  µ2 = −4β01 , η2 = 4 α11 + b21 ,  µ3 = 4α01 , η3 = 2 α21 + β11 − a21 − b11 . Again, using the identities (16) equation ∆2 (r) = 0 is transformed into the polynomial equation x3 (1 + x4 )p2 (x) = 0, where 7 X (19) p2 (x) = λk uk (x), k=0 with u0 (x) = x5 , u1 (x) = x4 , u2(x) = x6 , u3 (x) = x7 , u4 (x) = 1, u5 (x) = x2 , u6(x) = x3 − x7 , u7 (x) = x + 3x9 , and coefficients  λ0 = − 8π a12 + b22 + α12 + (a21 − α21 )(α11 + β21 ) + β22 ,  λ1 =8 3β02 − 3b02 + 4b21 α01 + 3a21 β01 − 3α21 β01 − 4α01 β21 ,  λ2 =24 a02 − a21 α01 − α02 + α01 α21 − 2β01 (α11 + β21 ) , λ3 =8π(α11 + β21 )2 ,  (20) λ4 = − 8 b02 + b11 β01 − β02 − β01 β11 + 2α01 (α11 + β21 ) ,  λ5 =8 a02 + b11 α01 − α02 − α01 β11 − 4β01 (b21 + 2α11 + β21 ) ,   λ6 =4π α11 + β21 b21 + 2α11 + β21 , λ7 = − π 2a12 + 2b22 + (a21 − b11 − α21 + β11 )(α11 + β21 )  + 2(α12 + β22 ) . MELNIKOV ANALYSIS IN NONSMOOTH DIFFERENTIAL SYSTEMS 19 Considering the ordered set of functions F = [u0 , u1, . . . , u7 ] and computing the Wronskians defined in (14) we obtain W0 (x) = x5 , W1 (x) = −x8 , W2 (x) = −2x12 , W3 (x) = −12x16 , W4 (x) = −10080x12 , W5 (x) = −2419200x9 , W6 (x) = −174182400x6 , W7 (x) = −125411328000(1 + 189x8 ). From Theorem 2, the maximum number of positive roots of poly- nomial p2 is seven and the upper bound is reached. From (20), it is clear that for any choice of λ’s it is possible to obtain the corresponding perturbation. This concludes the proof of Proposition 5.  Example 2. In what follows, we exhibit a concrete perturbation pre- senting seven limit cycles for |ε| = 6 0 sufficiently small. Choosing 434699860 λ0 = 1, λ1 = − , 124987243 18912094 5527397195 (21) λ2 = − , λ3 = , 124987243 874910701 1929240 613409686 λ4 = , λ5 = − , 1582117 124987243 11037105503 6534 λ6 = , λ7 = − , 1749821402 874910701 the polynomial (19) becomes 7 (39204x2 + 1097712x + 423361) Y p2 (x) = − (x − k). 1749821402 k=1 So, p2 (x) has seven positive roots at 1, 2, . . . , 7 (see Figure 4). 0.6 0.4 0.2 1 2 3 4 5 6 7 -0.2 - -0.4 - -0.6 - Figure 4. Graph of p2 (x). 20 J. BASTOS, C. BUZZI, J. LLIBRE, D. D. NOVAES The system composed by equations (20) with the parameters given by (21) has infinitely many solutions. One of these solutions is given by a01 = a01 , a11 = a11 , a21 = a21 , b01 = b01 , b21 = b21 , α01 = α01 , α11 = α11 , β01 = β01 , a02 = a02 , a12 = a12 , b02 = b02 , where q 5527397195π 822816127288564767 1749821402 a01 = , r 9860485382218537051 1 5527397195 17688887 a11 = − √ , 2 1749821402π 2 9671957909165767390π 874936837 a21 =√ , 9671957909165767390π q 1014941851585915217 11054794390π 874910701 b01 =− , 29581456146655611153 17688887 b21 = √ , 2 9671957909165767390π q 822816127288564767 5527397195π 1749821402 α01 = , r 9860485382218537051 1 5527397195 α11 =− , 2 1749821402π q 1014941851585915217 11054794390π 874910701 β01 =− , 29581456146655611153 26106731320711553594487715355 a02 = , 103524530135484878237818593012 874962973 a12 = , 6999285608π 1918068234277679948089129635 b02 = , 17254088355914146372969765502 and all other coefficients equal to zero. Notice that conditions (18) have been taken into account. So, for |ε| = 6 0 sufficiently small, system  y + (a01 + a11 x + a21 y)ε + (a02 + a12 x)ε2 if y ≥ x3 , ẋ = y + (b01 + b21 y)ε + b02 ε2 if y ≤ x3 ,  −x + (α01 + α11 x)ε if y ≥ x3 , ẏ = −x + β01 ε if y ≤ x3 has seven limit cycles. MELNIKOV ANALYSIS IN NONSMOOTH DIFFERENTIAL SYSTEMS 21 Acknowledgements We thank to the referee for the helpful comments and suggestions. JLRB has been partially supported by the Brazilian FAPESP grant 2013/24541-0. CAB has been partially supported by the Brazilian FAPESP grant 2013/24541-0 and the Brazilian Capes grant 88881.068462/2014- 01. JL has been partially supported by a FEDER-MINECO grant MTM2016-77278-P, a MINECO grant MTM2013-40998-P, and an AGAUR grant number 2014SGR-568. DDN has been partially supported by the Brazilian FAPESP grant 2016/11471-2. References [1] J. Artes, J. Llibre, J. C. Medrado, and M. A. Teixeira, Piecewise linear differential systems with two real saddles., Math. Comput. Simulation, to appear (2017). [2] D. D. C. Braga and L. F. Mello, More than three limit cycles in dis- continuous piecewise linear differential systems with two zones in the plane, Internat. J. Bifur. Chaos Appl. Sci. Engrg., 24 (2014), pp. 1450056, 10. [3] A. Buica and J. Llibre, Averaging methods for finding periodic orbits via Brouwer degree, Bull. Sci. Math., 128 (2004), pp. 7–22. [4] C. Buzzi, C. Pessoa, and J. Torregrosa, Piecewise linear perturbations of a linear center, Discrete Contin. Dyn. Syst., 33 (2013), pp. 3915–3936. [5] A. Colombo, M. Jeffrey, J. T. Lázaro, and J. M. Olm, Extended Abstracts Spring 2016 - Nonsmooth Dynamics, vol. 8 of Trends in Mathemat- ics Research Perspectives CRM Barcelona, Birkhäuser, Springer International Publishing, 2017. [6] D. de Carvalho Braga and L. F. Mello, Limit cycles in a family of discontinuous piecewise linear differential systems with two zones in the plane, Nonlinear Dynam., 73 (2013), pp. 1283–1288. [7] M. di Bernardo, C. J. Budd, A. R. Champneys, and P. Kowalczyk, Piecewise-smooth dynamical systems, vol. 163 of Applied Mathematical Sci- ences, Springer-Verlag London, Ltd., London, 2008. [8] E. Freire, E. Ponce, and F. Torres, Canonical discontinuous planar piecewise linear systems, SIAM J. Appl. Dyn. Syst., 11 (2012), pp. 181–211. [9] , The discontinuous matching of two planar linear foci can have three nested crossing limit cycles, Publ. Mat., 58 (2014), pp. 221–253. [10] F. Giannakopoulos and K. Pliete, Planar systems of piecewise linear differential equations with a line of discontinuity, Nonlinearity, 14 (2001), pp. 1611–1632. [11] M. Han, V. G. Romanovski, and X. Zhang, Equivalence of the Melnikov function method and the averaging method, Qual. Theory Dyn. Syst., 15 (2016), pp. 471–479. [12] M. Han and W. Zhang, On Hopf bifurcation in non-smooth planar systems, J. Differential Equations, 248 (2010), pp. 2399–2416. 22 J. BASTOS, C. BUZZI, J. LLIBRE, D. D. NOVAES [13] S.-M. Huan and X.-S. Yang, Existence of limit cycles in general planar piecewise linear systems of saddle-saddle dynamics, Nonlinear Anal., 92 (2013), pp. 82–95. [14] , On the number of limit cycles in general planar piecewise linear systems of node-node types, J. Math. Anal. Appl., 411 (2014), pp. 340–353. [15] J. Itikawa, J. Llibre, and D. D. Novaes, A new result on averaging theory for a class of discontinuous planar differential systems with applications, To appear in R. Matemática Iberoamericana, (2017). [16] S. Karlin and W. J. Studden, Tchebycheff systems: With applications in analysis and statistics, Pure and Applied Mathematics, Vol. XV, Interscience Publishers John Wiley & Sons, New York-London-Sydney, 1966. [17] L. Li, Three crossing limit cycles in planar piecewise linear systems with saddle-focus type, Electron. J. Qual. Theory Differ. Equ., (2014), pp. No. 70, 14. [18] J. Llibre, A. C. Mereu, and D. D. Novaes, Averaging theory for discon- tinuous piecewise differential systems, J. Differential Equations, 258 (2015), pp. 4007–4032. [19] J. Llibre, D. D. Novaes, and C. A. B. Rodrigues, Averaging theory at any order for computing limit cycles of discontinuous piecewise differential systems with many zones, Phys. D, 353/354 (2017), pp. 1–10. [20] J. Llibre, D. D. Novaes, and M. A. Teixeira, Higher order averaging the- ory for finding periodic solutions via Brouwer degree, Nonlinearity, 27 (2014), pp. 563–583. [21] , Limit cycles bifurcating from the periodic orbits of a discontinuous piece- wise linear differentiable center with two zones, Internat. J. Bifur. Chaos Appl. Sci. Engrg., 25 (2015), pp. 1550144, 11. [22] , Maximum number of limit cycles for certain piecewise linear dynamical systems, Nonlinear Dynam., 82 (2015), pp. 1159–1175. [23] , On the birth of limit cycles for non-smooth dynamical systems, Bull. Sci. Math., 139 (2015), pp. 229–244. [24] J. Llibre and E. Ponce, Three nested limit cycles in discontinuous piecewise linear differential systems with two zones, Dyn. Contin. Discrete Impuls. Syst. Ser. B Appl. Algorithms, 19 (2012), pp. 325–335. [25] J. Llibre, M. A. Teixeira, and J. Torregrosa, Lower bounds for the maximum number of limit cycles of discontinuous piecewise linear differential systems with a straight line of separation, Internat. J. Bifur. Chaos Appl. Sci. Engrg., 23 (2013), pp. 1350066, 10. [26] J. Llibre and X. Zhang, Limit cycles for discontinuous planar piecewise linear differential systems. Preprint, 2016. [27] O. Makarenkov and J. S. W. Lamb, Dynamics and bifurcations of non- smooth systems: a survey, Phys. D, 241 (2012), pp. 1826–1844. [28] D. D. Novaes and E. Ponce, A simple solution to the Braga-Mello con- jecture, Internat. J. Bifur. Chaos Appl. Sci. Engrg., 25 (2015), pp. 1550009, 7. [29] D. D. Novaes and J. Torregrosa, On extended Chebyshev systems with positive accuracy, J. Math. Anal. Appl., 448 (2017), pp. 171–186. [30] D. J. W. Simpson, Bifurcations in piecewise-smooth continuous systems, vol. 70 of World Scientific Series on Nonlinear Science. Series A: Monographs and Treatises, World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2010. MELNIKOV ANALYSIS IN NONSMOOTH DIFFERENTIAL SYSTEMS 23 [31] M. A. Teixeira, Perturbation theory for non-smooth systems, in Mathematics of complexity and dynamical systems. Vols. 1–3, Springer, New York, 2012, pp. 1325–1336.