This is a preprint of: “Lower bounds for the number of limit cycles in a generalized Rayleigh- Liénard oscillator”, Rodrigo D. Euzébio, Jaume Llibre, Durval J. Tonon, Nonlinearity, vol. 35, 3883–3906, 2022. DOI: [10.1088/1361-6544/ac7691] LOWER BOUNDS FOR THE NUMBER OF LIMIT CYCLES IN A GENERALIZED RAYLEIGH-LIÉNARD OSCILLATOR RODRIGO D. EUZÉBIO1 , JAUME LLIBRE2 AND DURVAL J. TONON1 Abstract. In this paper a generalized Rayleigh-Liénard oscillator is consider and lower bounds for the number of limit cycles bifurcating from weak focus equilibria and saddle connections are pro- vided. By assuming some open conditions on the parameters of the considered system the existence of up to twelve limit cycles is provided. More precisely, the approach consists in perform suitable changes in the sign of some specific parameters and apply Poincaré-Bendixson Theorem for assure the existence of limit cycles. In particular, the algorithm for obtaining the limit cycles through the referred approach is explicitly exhibited. The main techniques applied in this study are the Lya- punov constants and the Melnikov method. The obtained results contemplate the simultaneity of limit cycles of small amplitude and medium amplitude, the former emerging from a weak focus equilibrium and the latter from homoclinic or heteroclinic saddle connections. 1. Introduction 1.1. Historical facts and equations of Rayleigh and Liénard. Ordinary differential equa- tions (ODEs) have been largely studied in mathematics since the invention of Calculus back in 17th century. Since that the theory have proved to be very accurate to model real problems from mechanics movements and chemical reactions to social and financial sciences. The interest by ODEs gained even more attraction after the remarkable work of Poincaré entitled Mémoire sur les courbes définies par une équation différentielle, see [17]. This paper, dated 1882, is consider one of the start- ing points of the so called qualitative theory of ODEs. In particular, Poincaré formally introduced the concept of limit cycle, an isolated periodic orbit inside the set of all periodic orbits of an ODE, and exhibited an ad hoc example of ODE presenting a limit cycle without any connection to some concrete problem. However, the first reported case of a limit cycle surging from a real model ODE was probably provided by Rayleigh in 1877 in his study on the oscillations of a violin string, see [18]. Posteriorly in 1908 another example of limit cycle emerged from a series of works of Poincaré addressing wireless telegraphy, although the most recognized example of a limit cycle is due to Van der Pol on the electrical circuits in 1927, see [21]. The goal of this paper is to study the existence and the number of limit cycles in a system which is a generalization of the Rayleigh and Van der Pol systems. The equation proposed by Rayleigh which is nowadays known as Rayleigh equation is (1) ẍ + ax + ε(c3 + c4 ẋ2 )ẋ = 0 where ε is a small parameter. We also point to [13] for some historical facts about Rayleigh work. 2010 Mathematics Subject Classification. 34C23 and 34C25 and 37G15. Key words and phrases. Rayleigh Liénard oscillator and limit cycles and Lyapunov constants and Melnikov function. 1 2 R. D. EUZÉBIO, J. LLIBRE, D. J. TONON For the sake of applications of non-linear systems it is usually interesting to assume that the unperturbed part of (1) has a potential of the form V (x) = ax2 +bx4 so the total energy is ẋ2 +V (x). This is achieved by adding the term 2bx3 to the last equation obtaining then the equation (2) ẍ + ax + 2bx3 + ε(c3 + c4 ẋ2 )ẋ = 0, which is the generalized Rayleigh system. By replacing the term c4 ẋ2 in the Rayleigh equations by c2 x2 + c1 x4 we obtain the famous generalized Van der Pol or Liénard equation. Therefore, by combining generalized Rayleigh and Liénard equations we obtain (3) ẍ + ax + 2bx3 + ε(c3 + c2 x2 + c1 x4 + c4 ẋ2 )ẋ = 0. Last equation is common referred in the literature as mixed generalized Rayleigh-Liénard equation. In the general case, one can study the problem ẍ + ax + 2bx3 = εf (x, ẋ), see for instance the work of Guckenheimer and Holmes in [10]. In this paper we study the limit cycles for the case f (x, ẋ) = (c3 +c2 x2 +c1 x4 +c4 ẋ2 +c5 x6 +c6 ẋ4 )ẋ so that the mixed generalized Rayleigh-Liénard equation becomes a particular case of the equation we deal with. The equation we study is then (4) ẍ + ax + 2bx3 − ε(c3 + c2 x2 + c1 x4 + c4 ẋ2 + c5 x6 + c6 ẋ4 )ẋ = 0, that is equivalent to the system ẋ = y, (5) ẏ = −ax − 2bx3 + εQ(x, y). where Q(x, y) = (c3 + c2 x2 + c1 x4 + c4 y 2 + c5 x6 + c6 y 4 )y. Our purpose is to obtain a lower bound for the number of limit cycles of (5). The problem of finding limit cycles involves several methods and approaches, some of them are briefly summarized in what follows. In [2] the authors study equation (3) by using the harmonic balance and Krylov-Bogoliubov methods to obtain up to two limit cycles. System (3) was also studied by Lynch in [14] by using Lyapunov constants obtaining three limit cycles. Then, the same authors of [2] using harmonic balance method and elliptic functions presented an example with seven limit cycles for the case a < 0 and b > 0, see [3]. In [22] the authors prove that system (3) can have eight limit cycles, improving the lower bound obtained in [3]. The authors use both Lyapunov constants as well as Melnikov method to get the results for the case ab < 0. It also worth to mention the recent work [5] on which the authors study the global dynamics from a particular case of system (2), namely, assuming ε = c4 = 2b = 1. The main result of that paper provides conditions on a and c3 so that system (2) presents pitchfork, Hopf, homoclinic and double limit cycle bifurcations, which would be a very difficult - if possible - task considering arbitrary parameters or a more general system (3). The study of lower and upper bounds for the number of limit cycles is widely considered in the literature. An usual approach is to consider suitable perturbations on structural unstable systems and then to look at the bifurcations taking place. In this case, the bounds for the numbers of LIMIT CYCLES IN A GENERALIZED RAYLEIGH LIÉNARD OSCILLATOR 3 limit cycles that can bifurcate provide an idea about the co-dimension of the problem. In the local scenario, these bifurcations usually occur at degenerated Hopf bifurcations equilibria and they can be studied through the Lyapunov constants, we mention the work of Torregrosa and coauthors, see for instance [8] and [12]. The non-local case can be studied through some techniques, an effective one being the Melnikov function which can be applied in different ways. One of them is to consider global bifurcations from non hyperbolic limit cycles or saddle connections, see for instance [22] and [19]. Another way to apply Melnikov method is through the analysis of critical and regular level curves associated to some Hamiltonian function. This is a more general approach because it contemplates both local and global scenarios. In particular, the study of Melnikov functions emerging from perturbations of Hamiltonian sys- tems have been well established in the work of Petrov in a series of papers published some decades ago, see for instance [15], [16] and references therein. In those papers Petrov provides some up- per bounds for the number of limit cycles bifurcating from distinct configurations of Hamiltonian functions in terms of the degree of the perturbations. In particular, Petrov use Tchebyshev spaces to associate the number of zeros of the Melnikov function to the number of zeros of some linear combinations involving elliptic integrals, and they represent a fine advance in the study of limit cycles. Nevertheless, a sharp bound is obtained in [7] where the authors use an approach similar to Petrov’s to provide a new upper bound for the number of limit cycles. 1.2. The main result of the paper. In this paper we study system (5) which generalizes equation (3). We consider the case where a and b have opposite sign and apply both Lyapunov and Melnikov methods as the authors did in [22]. The approach consists in to consider suitable perturbations on the coefficients of the system to produce changes in the sign of the Lyapunov constants (derivatives of the Poincaré map) and Melnikov function. The limit cycles essentially bifurcate from weak focus equilibria and from heteroclinic and homoclinic loops containing those weak focus. We highlight that the lower bounds provided in this paper agree with Petrov’s statements mentioned before. We distinguish between the limit cycles that bifurcate from the weak focus equilibria from the ones that bifurcate from saddle connections, denoting them as small amplitude limit cycles, or medium amplitude limit cycles, respectively. Therefore we say that system (5) presents a configuration (i, j) of limit cycles if there exists i limit cycles of small amplitude and j limit cycles of medium amplitude. The main contributions of this paper can be summarized as follows: first, we consider a quite general system without assuming any hypotheses on the parameters except the condition ab < 0. So we are able to obtain from one to twelve limit cycles in a region of the phase portrait containing the equilibria. On the other hand, we explicitly exhibit the algorithm to obtain those limit cycles and we provide the conditions for the realization of those number of limit cycles. This important part of the process for obtaining the limit cycles is in general omitted in the literature. This is the case in [22], where the authors only pointed out the approach which a priori does not guarantee the realization of the limit cycles as claimed. We stress, however, that although paper [22] contain some minor miscalculations, assuming c5 = c6 = 0 our results points to an upper bound of eight limit cycles, which is the quantity obtained in that paper. Besides then, we emphasize that the procedure considered in the present paper, based in the Lyapunov constants and Melnikov method, allow us to consider the simultaneous existence of limit cycles bifurcating from 4 R. D. EUZÉBIO, J. LLIBRE, D. J. TONON the equilibrium of center type (small amplitude) and from the homoclinic/heteroclinic connection (medium amplitude) that represents a critical point of the Melnikov function. The main result of this paper is the following. Theorem 1. Consider system (5). Then there exists suitable values of parameters realizing the following configurations of limit cycles: (a) (s, m) if a > 0 and b < 0, where m ∈ {0, 1, 2, 3}, s ∈ {0, 1, 2, 3, 4, 5} and s + m ≤ 5. (b) (2s, 3m + k) if a < 0 and b > 0, where m ∈ {0, 1, 2}, k ∈ {0, 1, 2}, s ∈ {0, 1, 2, 3, 4} and 2s + 3m + k ≤ 12. The rest paper is organized as follows: In Section 2 the canonical forms that we consider in this paper are presented. In Section 3 we present the general facts about Lyapunov constants and Melnikov method. The Lyapunov constants and the Melnikov functions related with the generalized Rayleigh-Liénard system are given in Sections 4 and 5, respectively. Finally, in Section 6 we state and prove some auxiliary results and Theorem 1. 2. Canonical forms for the generalized mixed Rayleigh-Liénard oscillator In the following we apply a linear change of coordinates to equation (4), providing a simpler expression for the system that models a generalized mixed Rayleigh-Liénard oscillator. Lemma 1. Consider system (5). The following statements holds. (a) If a > 0 and b < 0 then system (5) is topologically equivalent to system ẋ = y, (6) 2b 3 ẏ = −x − x + εQ1 (x, y), a2 √ √ where Q1 (x, y) = (d3 + d2 x2 + d1 x4 + d4 y 2 + d5 x6 + d6 y 4 )y, d1 = c1 / 5 a, d2 = c2 / 3 a, /, d3 = √ √ √ √ c3 / a, d4 = c4 / a, d5 = c5 / 7 a and d6 = c6 / a. (b) If a < 0 and b > 0 then system (5) is topologically equivalent to system ẋ = y, (7) √ 3 b 2 b ẏ = −x + x − 2 x3 + εQ2 (x, y), 2a 2a where Q2 (x, y) = (e3 + e7 x + e2 x2 + e8 x3 + e1 x4 + e9 x5 + e5 x6 + e4 y 2 + e6 y 4 )y, e1 = 2c1 b − 15ac5 15a2 c5 − 12ac1 b + 4c2 b2 a3 c5 − 2a2 c1 b + 4ac2 b2 − 8c3 b3 √ , e2 = √ , e3 = − √ √ , e4 = 8 2(−a)5/2 b 8 2(−a)3/2 b2 8 2 −ab3 c c c6 √ 4 , e5 = √ 5 and e6 = √ . Besides then the parameters e7 , e8 and e9 are −2a 8 2(−a) 7/2 −2a linearly dependent with respect  to the parameters e1 , . . . , e6 and satisfies 2a 48a4 e5 − 4a2 be1 + b2 e2 40a3 e5 − 4abe1 6ae5 e7 = − 5/2 , e8 = 3/2 and e9 = − √ . b b b Proof. Consider system (5). We initially consider the case a > 0 and b < 0. Applying the change of √ √ e = ax, ye = y and e coordinates x t = at we obtain system (6), where we have removed the tilde in the expression of the system. LIMIT CYCLES IN A GENERALIZED RAYLEIGH LIÉNARD OSCILLATOR 5 r  −a For the case a < 0 and b > 0, we initially translate the point p1 = , 0 to the origin by r 2b −a e=x− the change of coordinates x and ye = y. So system (5) becomes 2b ẋ = y, (8) √  ẏ = x 3 −2abx + 2a − 2bx2 + εQ(x, y), y  √  √  where Q(x, y) = 3 8 −2ababx c1 + 5c5 x2 −6 −2aba2 c5 x − a3 c5 + 2a2 b c1 + 15c5 x2 √ 8b  −8 −2abb2 x 2c x2 + c + 3c x4 − 4ab2 (6c x2 + c +15c x4 ) + 8b3 (c x4 + c x2 + c + c y 2 + c x6  1 2 5 1 2 5 1 2 3 4 5 +c6 y 4 ) and as in the previous case, we remove the tildes in the expression of the system. e √ √ After we consider the rescaling x e = −2ae x, e ye = ye and e t = −2at. Applying this change of coordinates to system (8), we get system (7), where we remove the tildes in the expression of the system. This completes the proof of the lemma.  3. Main techniques to analyze small and medium amplitude limit cycles To analyze the types and the number of limit cycles that can bifurcate near to the perturbed center and from the homoclinic/heteroclinic connections we make use of two main techniques, namely, Lyapunov constants and Melnikov functions. We apply these two powerful techniques separately and posteriorly we combine them to obtain simultaneity of small and medium amplitude limit cycles. In this way, exhibiting the Lyapunov constants with the local basis and improving the Bautin’s algorithm we are able to analyze the number of small amplitude limit cycles, that bifurcates from the perturbed center when system (5) presents a homoclinic connection (a < 0, b > 0) or a heteroclinic one (a > 0, b < 0). Next, we present the Melnikov functions that characterize the existence of the homoclinic/hete- roclinic connections, as well as the stability of these connections. According to this approach of changing the stability of the connections and breaking them ultimately, we are able to provide con- ditions for the existence of limit cycles that bifurcates from the homoclinic/heteroclinic connections, which we call medium amplitude limit cycles. All these results are given in terms of the param- eters of the systems. Finally, we analyze the simultaneous existence of the bifurcations of small and medium limit cycles, near to the perturbed center and the homoclinic/heteroclinic connection, respectively. In the following we present these techniques. 3.1. Lyapunov constants. In the following we briefly present the approach to deal with the limit cycles emerging from changing signs in the Lyapunov constants. It can be founded for instance in [1, 4, 6]. Consider differential systems (9) ẋ = λx − y + P (x, y), ẏ = x + λy + Q(x, y), where P and Q are polynomials without constant and linear terms. Then the origin of system (9) is a weak focus if λ = 0. The limit cycles that bifurcate from a weak focus are called small amplitude limit cycles. When λ = 0 we denote a local Lyapunov function by V , defined in a neighborhood of the origin. Then the origin is a weak focus it is stable or unstable if V̇ < 0 or V̇ > 0, respectively, where 6 R. D. EUZÉBIO, J. LLIBRE, D. J. TONON V̇ denotes the rate of changes of V along the trajectories of (9). The expression of V can be constructed, see [4] and [9], and it is of the form V̇ = η2 (x2 + y 2 ) + η4 (x2 + y 2 )2 + · · · + η2k (x2 + y 2 )2k + . . . , where η2k is a polynomial in the coefficients of the polynomials P and Q. We define η2k the k-th focal value. In the following we get that the weak focus is stable, unstable if the first non-zero focal value is negative, positive, respectively, see [4]. When the origin is a weak focus it is a center if and only if η2k = 0 for all k. Moreover the stability of the origin is determined by the sign of the first non–zero focal value. As η2k is relevant only when η2l = 0 for l < k, we put η2 = η4 = . . . = η2k−2 = 0 in the expression for η2k . Instead of working with the η2k following to many authors we prefer to work with V2k+1 = 2πη2k , and call it the n–th Lyapunov constant. For ε > 0 and small consider the interval J = {(x, 0) : 0 ≤ x < ε} and the Poincaré return map x 7→ h(x) defined from J → {(x, 0) : 0 ≤ x}. It assigns to x the abscissa h(x) of the point where the orbit of the differential system (9) starting at the point (x, 0) ∈ J first returns to the positive half-axis {(x, 0) : 0 ≤ x}. Then the displacement function is defined as x 7→ d(x) from J → R by d(x) = h(x) − x. Therefore the orbit of system (9) through the point (x, 0) 6= (0, 0) is periodic if and only if x is a zero of the displacement function. Clearly the Lyapunov constants are related with the coefficients of the displacement function, because the origin of system (9) is a center if and only if the displacement function is identically zero, if and only if the Lyapunov constants V2k+1 = 0 for k ≥ 1. In fact it is known that     X∞ X X∞ d(x) = λx a0 + a1j xj  + V2k+1 x2k+1 1 + a2k+1 j xk  , j=1 k≥1 j=1 for |x| < ε, where e2πλ − 1 a0 = = 2π + O(λ), λ and the aj2k+1 for k = 1, 2, . . . are analytic functions in λ and in the coefficients of the polynomials P and Q. For more details in the displacement function see [20]. Since P and Q are polynomials, by the Hilbert basis theorem there is a constant m such that V2k+1 = 0 for all k ≥ 1 if and only if V2k+1 = 0 if k = 1, . . . , m. Therefore it is necessary to compute only a finite number of the Lyapunov constants, though with few exceptions for any given case it is unknown a priori how many are required. Definition 1. We say that the origin is a weak focus of order k of system (9) if η2 = · · · = η2k = 0 and η2k+2 6= 0. Remark 1. By the previous definition, the origin is a weak focus of order k of system (9) if V3 = · · · = V2k−1 = 0 and V2k+1 6= 0. In [4] the authors stated that if a system presents a weak focus at the origin of order k at most k small amplitude limit cycles can bifurcate from the origin under perturbation of the system. In LIMIT CYCLES IN A GENERALIZED RAYLEIGH LIÉNARD OSCILLATOR 7 the following we describe briefly how to provide a convenient perturbation of the original system presenting a weak focus of order k at the origin to produce k small amplitude limit cycles. The expressions η2k provides the Lyapunov constants V2k−1 . We assume that V1 = 1, V3 = · · · = V2k−1 = 0 and V2k+1 6= 0. Without loss of generality we assume that V2k+1 < 0. Therefore the origin is stable. Consider Γ1 be a level curve of V which is sufficiently close to the origin. So the flow of system (9) is inward across it. Now consider a suitable perturbation S1 of system (9) such that the Lyapunov constants satisfy V3 = · · · = V2k−3 = 0 and V2k−1 > 0. Therefore, the origin now is unstable for S1 . As S1 is sufficiently close to system (9) then the flow remains inward across Γ1 . Consider W1 the Lyapunov function of S1 , take Γ2 a level curve of W1 , inside of the region limited by Γ1 and sufficiently close to the origin in such a way that the flow of S1 is outward across of Γ2 . Therefore by the Poincaré-Bendixson Theorem we conclude that there exists a limit cycle of S1 between Γ1 and Γ2 . In the following working in a similar way we consider a convenient perturbation S2 of S1 with analogous properties. In this way at most k limit cycles can be produced by convenient perturbations of system (9). Note that the Lyapunov constants must satisfy: |V2i−1 |  |V2i+1 | and V2i−1 · V2i+1 < 0, for i = 3, . . . , k. If all the constants V2k+1 are zero then the origin is a center. 3.2. Melnikov method. In order to present the Melnikov method, note that system (5) with ε = 0 is a Hamiltonian system with the Hamiltonian function 1 H(x, y) = (y 2 + ax2 + bx4 ). 2 If ab > 0 then the origin isr the unique  equilibrium point r and if ab < 0 there exists three equilibrium −a −a points: O = (0, 0), p1 = , 0 and p2 = − , 0 . In this way if a > 0 and b < 0, then 2b 2b O = (0, 0) is a center point and p1 , p2 are saddle points, and if a < 0 and b > 0 then O = (0, 0) is saddle point and p1 , p2 are center points. As system (5) with ε = 0 is Hamiltonian, then for the case a > 0 and b < 0 we have a heteroclinic loop between the two saddle points p1 and p2 . On the other hand in the case a < 0 and b > 0 the system presents a homoclinic loop. It is important to note that system (5) is invariant under the transformation (x, y) 7→ (−x, −y), so its phase portrait is symmetric with respect to the origin, see Figure 1. p2 O p1 p2 p1 O (a) (b) Figure 1. The topological structure of system (5) with ε = 0. In (a) we get a > 0, b < 0, and in (b) we have a < 0, b > 0. 8 R. D. EUZÉBIO, J. LLIBRE, D. J. TONON We note that for ε 6= 0 system (5) is no longer Hamiltonian so if ε is sufficiently small generically the saddle connections are broken and the center structure is destroyed. Therefore in both scenarios we can eventually have the birth of limit cycles emerging from those equilibria or loops. n1 Aε1 Aε2 p Aε1 p n1 p2 O p1 O p1 Aε2 (a) (b) Figure 2. In figures (a) and (b) are presented the geometric approach of the con- struction of Melnikov functions of systems (6) and (7), respectively. More precisely, consider Aε1 the unstable manifold of p2 and Aε2 the stable manifold of p1 , see Figure 2. Let Aε be the heteroclinic, homoclinic loop of systems (6), (7), resp., and p ∈ Aε . Consider (Hx (p), Hy (p)) n1 = , ||(Hx (p), Hy (p))|| then Melnikov function is −−−→ d(ε, M1 ) = −hn1 , Aε1 Aε2 i, see Figure 2, and for more details about the homoclinic and heteroclinic loops see [19]. We notice that system (6) has the same equilibrium points than system (5) but in this case the Hamiltonian function associated to its unperturbed part is given by b 1 H1 (x, y) = ax2 + x4 + y 2 . 2 2 Since the orbits of system (6) with ε = 0 lie on the h-levels of function H1 , we see that for h = −a2 /8b it has a heteroclinic loop A formed by the saddle equilibrium points and two orbits A1 and A2 connecting them, that is, A = {p1 } ∪ A1 ∪ {p2 } ∪ A2 . From the expression of H1 we obtain the analytical expression of these arcs: r r −b 2 a A1,2 : y = y(x) = ∓ x ± , a −4b r a where |x| ≤ − . In this case we note that A surrounds the center equilibrium O = (0, 0). 2b     a 2a On the other hand, the points O = (0, 0), q1 = √ , 0 and q2 = √ , 0 are equilibrium b b points for system (7) being √ x2 b 3 b 1 H2 (x, y) = − x + 2 x4 + y 2 2 2a 8a 2 LIMIT CYCLES IN A GENERALIZED RAYLEIGH LIÉNARD OSCILLATOR 9 the Hamiltonian function of the system for ε = 0. Doing a translation of the original system for obtaining system (7), q1 is a saddle when ε = 0, the equilibria O and q2 are centers. The level h = a2 /8b of Hamiltonian function H2 contains two homoclinic loops Ll and Lr and the saddle point centered at q1 . We call L = Ll ∪ {q1 } ∪ Lr . In this case Ll and Lr are determined by the following relations: s √ √ a2 + 2a bx − bx2 Ll : y = sign (y)(−a + bx) , a2 b √ a(1 + 2) a if √ ≤ x ≤ √ and b b s √ √ a2 + 2a bx − bx2 Lr : y = sign (y)(a − bx) , a2 b √ a a(1 − 2) if √ ≤ x ≤ √ . Observe that the homoclinic loop Li surrounds the equilibria O and q2 . b b Now we denoted by Q2 (x, y) the perturbed part of system (7), that is, Q2 (x, y) = (e3 + e7 x + e2 x2 + e8 x3 + e1 x4 + e9 x5 + e4 y 2 + e5 x6 + e6 y 4 )y. 4. Preliminary to obtain small limit cycles using Lyapunov constants In this section we present the Lyapunov constants for each system. Lemma 2. Consider system (6) with d3 = 0. The first Lyapunov constant is V1 = 1, then 1 (10) V3 = − πε(d2 + 3d4 ). 4 If V3 = 0, then πε(a2 d1 + 5a2 d6 + 6bd4 ) (11) V5 = − . 8a2 If V3 = V5 = 0, then  πε −6a2 d34 ε2 + 5a2 d5 + 80bd6 (12) V7 = − . 64a2 If V3 = V5 = V7 = 0, then  3πε 15a4 d24 d6 ε2 + 12a2 bd34 ε2 − 35b2 d6 (13) V9 = 160a4 If V3 = V5 = V7 = V9 = 0, then 3πd34 ε3 12 6 6 8 2 4 4 4 4 2 2 6  (14) V11 = − 2 315a d4 ε + 2644a b d4 ε − 9065a b d4 ε − 7350b . 3200a4 3a4 d24 ε2 − 7b2 Proof. Consider the polynomial system ẋ = y + p(x, y), ẏ = −x + q(x, y). 10 R. D. EUZÉBIO, J. LLIBRE, D. J. TONON Applying the change of coordinates x = r cos θ and y = sin θ, we can write dr X (15) = vi (θ)ri , dθ i≥1 where vi (θ) are trigonometric polynomials in the variables cos θ and sin θ. Denoting by r(θ, r0 ) the solution of (15) satisfying r(0, r0 ) = r0 , then in a neighborhood of r = 0 we obtain X r(θ, r0 ) = u1 (θ)r0 + ui (θ)r0i , i≥2 with ui (0) = 0 for all i ≥ 2, then the Poincaré return map is X Π(r0 ) = r(2π, r0 ) = u1 (2π)r0 + ui (2π)r0i . i≥2 As the Poincaré map is analytic the condition Π(r0 ) ≡ r0 is equivalent to the fact that system (6) presents a center at the origin. Note that Π(r0 ) ≡ r0 if and only if ui (2π) = 0 for all i ≥ 2. As stated in [1] and by direct computations we obtain that u2k (2π) = 0 for all k a positive integer. At this moment we are able to obtain the Lyapunov constants for system (6). Consider system (6) and assuming that d3 = 0, then following the steps described previously, we get that the first Poincaré constant u1 (2π) = 1. and then we obtain 1 u3 (2π) = − πε(d2 + 3d4 ). 4 Assuming that d2 = −3d4 we obtain πε(a2 d1 + 5a2 d6 + 6bd4 ) u5 (2π) = − . 8a2 6bd4 If d2 = −3d4 and d1 = − − 5d6 , then a2  πε −6a2 d34 ε2 + 5a2 d5 + 80bd6 u7 (2π) = − . 64a2 6bd4 6d34 ε2 16bd6 If d2 = −3d4 , d1 = − − 5d6 and d5 = − , then a2 5 a2  3πε 15a4 d24 d6 ε2 + 12a2 bd34 ε2 − 35b2 d6 u9 (2π) = . 160a4 6bd4 6d34 ε2 16bd6 12a2 bd34 ε2 If d2 = −3d4 , d1 = − − 5d6 , d5 = − and d6 =  then a2 5 a2 5 7b2 − 3a4 d24 ε2 3πd34 ε3 12 6 6 8 2 4 4 4 4 2 2 6  u11 (2π) = − 2 315a d4 ε + 2644a b d4 ε − 9065a b d4 ε − 7350b . 3200a4 3a4 d24 ε2 − 7b2 As ε is a small parameter then u11 (2π) = 0 if and only if d4 = 0. So if besides d3 = 0 we assume 6bd4 6d34 ε2 16bd6 12a2 bd34 ε2 d2 = −3d4 , d1 = − 2 − 5d6 , d5 = − , d 6 =  and d4 = 0 then we a 5 a2 5 7b2 − 3a4 d24 ε2 obtain that all the others coefficients di are null and therefore system (6) is a Hamiltonian system, LIMIT CYCLES IN A GENERALIZED RAYLEIGH LIÉNARD OSCILLATOR 11 and therefore the origin is a center. In what follows we denote the Lyapunov constant ui (2π) by Vi .  4a2 be1 − 48a4 e5 Lemma 3. Consider system (7) with e3 = 0 and e2 = . The first Lyapunov constant b2 is V1 = 1, then  πε −48a4 e5 + 4a2 be1 + 3b2 e4 (16) V3 = − . 4b2 If V3 = 0, then   1 3be4 8a2 e5 (17) V5 = πε + − 5e6 . 8 a2 b If V3 = V5 = 0, then   1 24b2 e4 (18) V7 = πε + 3e34 ε2 + 44e5 32 a4 If V3 = V5 = V7 = 0, then  9πe4 ε 2a12 e64 ε6 + 7a8 b2 e44 ε4 − 25a4 b4 e24 ε2 + 35b6 (19) V9 = . 1760a6 b3 4a2 be1 − 48a4 e5 Proof. Considering now system (7) with e3 = 0 and e2 = . Proceeding in a similar b2 way as in the proof of Lemma 2, the first Lyapunov constant is u1 (2π) = 1 and  we conclude that the 3 b2 e4 − 16a4 e5 expression of u3 (2π) is given in (16). Assuming that e1 = − we conclude that the  4a2 b 3 b2 e4 − 16a4 e5 8a4 e5 + 3b2 e4 expression of u5 (2π) is given in (17). If e1 = − and e6 = we conclude 4a2 b  5a2 b 3 b2 e4 − 16a4 e5 8a4 e5 + 3b2 e4 that the expression of u7 (2π) is given in (18). If e1 = − , e 6 = and  4a2 b 5a2 b 4 3 2 3 a e4 ε + 8b e42 e5 = − we conclude that the expression of u9 (2π) is given in (19). 44a4 Similarly we done in the previous case, u9 (2π) = 0 if and only if e4 = 0. In this way, assuming e4 = 0 and replacing this condition in the others parameters we conclude that all ei ’s are null and therefore system (7) is a Hamiltonian system. So, the origin is a center. As we done in the previous case, we denote the Lyapunov constant ui (2π) by Vi .  5. Preliminary to obtain medium limit cycles using Melnikov function We now study the persistence and stability of the saddle connections of systems (6) and (7). We start analyzing on which (non generic) conditions the homoclinic and heteroclinic loops persist for these systems with ε 6= 0. This is provided by the first order Melnikov functions associated to each system. Once we have established conditions to persistence, we will deal with the stability of such loops. 5.1. Persistence and stability of heteroclinic loop in system (6). Denote by Q1 (x, y) the perturbed part of system (6), that is, Q1 (x, y) = (d3 + d2 x2 + d1 x4 + d4 y 2 + d5 x6 + d6 y 4 )y. For system (6) we have the following result: 12 R. D. EUZÉBIO, J. LLIBRE, D. J. TONON Lemma 4. System (6) has a heteroclinic loop Aε close to A if and only if d2 = ϕ1 (d1 , d3 , d4 , d5 , d6 , ε) 12d4 5a4 d5 a2 (99d1 + 160d6 ) 10bd3 =− − + + + O(ε). 7 84b2 462b a2 We write Aε = {p1 } ∪ Aε1 ∪ {p2 } ∪ Aε2 , where Aε1,2 denotes the orbits coming from A1,2 for ε 6= 0. In what follows we prove Lemma 4 only for the arc A1 , because by symmetry A1 is broken if and only if A2 is broken. Proof. In order to obtain the conditions for the persistence of the heteroclinic loop of system (6) with ε = 0 it is sufficient to study the zeros of the first order Melnikov function M1 associated to the heteroclinic loop. As done in [19] this function depends on the parameters of the system and it is given by Z M1 = Q1 (x, y) dx, A1 where we highlight that M1 = M1 (a, b, d1 , . . . , d6 ). Using the expression of A1 obtained in Subsection 3.2 we get Z M1 = (d3 + d2 x2 + d1 x4 + d4 y 2 + d5 x6 + d6 y 4 )y dx A1 Z √a −2b = 2 (d3 + d2 x2 + d1 x4 + d4 y(x)2 + d5 x6 + d6 y(x)4 )y(x) dx 0 r r −b 2 a2 where y(x) = − x + . Replacing M1 into the expression of the integral we get a2 −4b Z √a   −2b a2 + 2bx2 2 4 6 d4 2 2 2 d6 2 2 4 − √ (d3 + d2 x + d1 x + d5 x − 2 (a + 2bx ) + (a + 2bx ) ) dx. 0 2 −a2 b 4a b 16a4 b2 Integrating we obtain 1 M1 = √ (198a6 bd1 − 924a4 b2 d2 + 9240a2 b3 d3 − −1584a4 b2 d4 − 55a8 d5 + 320a6 bd6 ) 13860 2b4 and finally solving this last equation with respect to the parameter d2 we are done.  Since system (5) is analytic the heteroclinic loop is isolated, so one can study its stability. Lemma 5. Assume that d2 = ϕ1 . The heteroclinic loop Aε is stable (respect. unstable) if the value 1 div(p2 ) =(−22a6 d5 + a4 b(33d1 − 40d6 ) + 198a2 b2 d4 − 924d3 b3 )ε. 231b3 in negative (respect. positive). Proof. The proof is straightforward from the computation of the divergence at the saddle point p2 .  If the condition of Lemma 4 occurs we have a simple loop. However if that divergence is zero in order to obtain the stability of the heteroclinic loop we must analyze higher orders of the return map, see [11]. In this direction we have the following result. LIMIT CYCLES IN A GENERALIZED RAYLEIGH LIÉNARD OSCILLATOR 13 Lemma 6. Assume that d2 = ϕ1 and div(p2 ) = 0, that is, 40d6 2a2 d5 6d4 b 28d3 b2 d1 = ϕ2 (d3 , d4 , d5 , d6 ) = + − 2 + . 33 3b a a4 Then the heteroclinic loop Aε is stable (resp. unstable) if the value d5 is 6b(−100a4 d6 + 231a2 bd4 + 308d3 b2 ) ϕ3 (d3 , d4 , d6 , ε) = ε 11a6 is negative (resp. positive). Proof. In [19] we see that the stability of Aε is given by the sign of the expression Z Z  ∂Q1 ∂Q1 ε dt = ε dt + O(ε) Aε ∂y A ∂y Thus, since ε 6= 0 is sufficiently small it is sufficient to study the sign of the right hand side of the last equality. Using the particular values of d1 and d2 we get ∂Q1 d3 (a2 + 2bx2 )(a2 + 14bx2 ) x2 (a2 + 2bx2 )(80bd6 + 11d5 (a2 + 6bx2 )) = + ∂y a4 132b2   2bx4 +5d6 y 4 + 3d4 −x2 − 2 + y 2 . a Due to the symmetry Z and using the parametrization of arc A1 obtained in Subsection 3.2 we get ∂Q1 the expression for ε dt is A1 ∂y ε −11a7 d5 − 600a5 bd6 + 1386a3 b2 d4 + 1848ab3 d3 ε √ 3564 −2b7 This last equation vanishes for d5 = ϕ3 (d3 , d4 , d6 , ε), so we are done.  In summary the function ϕ1 controls the existence or not of the heteroclinic loop. If the loop remains for ε 6= 0, the stability of it is determined by the sign of ϕ2 if the loop is simple, or by the sign of ϕ2 if not. Now we perform the same computations for the case a < 0 and b > 0. 5.2. Persistence and stability of homoclinic loop in system (7). Analogously to Subsection 5.1 we study the Melnikov function and the stability of the homoclinic loop when it persists for ε 6= 0. Lemma 7. System (7) has a homoclinic loop Lε close to L if and only if 5be3 6e4 1417a4 e5 a2 (1287e1 + 40e6 ) e2 = φ1 (e1 , e3 , e4 , e5 , e6 , ε) = + − + + O(ε). a2 7 21b2 231b Now we call Lε = Lεl ∪ {q1 } ∪ Lεr , where in this case Lεl,r denotes the orbit arcs coming from Ll,r for ε 6= 0 sufficiently small. Proof. The conditions for the persistence of the homoclinic loop is provided by the zeros of the first order Melnikov function M2 associated to the homoclinic loop, that is Z M2 = Q2 (x, y) dx. Lr 14 R. D. EUZÉBIO, J. LLIBRE, D. J. TONON From the expression of the arc Lr we have M2 into the form √ s √ √ Z a−a√ 2 h b 1 (a − bx)2 (a2 + 2a bx − bx2 )  6ae5 2 e3 + e2 x2 + e1 x4 − √ x5 + e5 x6 √a 2 a b b b 2 be √ √ +b2 e2 +48a4 e5 ) 3 e4 (a− bx)2 (a2 +2a bx−bx2 ) − 2a(−4a 1√ x+ −4abe√ 1 +40a e5 3 x + 2 4a b b5 b3 √ √ i bx)2 (a2 +2a bx−bx2 )2 + e6 (a− 2 16a b dx. By integrating M2 we get √ 2 2a M2 = (1155b3 e3 + 33a2 b2 (−7e2 + 6e4 ) − 15587a6 e5 + a4 b(1287e1 + 40e6 )). 3465b4 Then solving this last equation with respect to the parameter e2 the result follows.  The proof of the next result is straightforward. Lemma 8. Assume that e2 = φ1 . The homoclinic loop Lε is stable (resp. unstable) if the value 2 div(q1 ) = − (462b3 e3 + 99a2 b2 e4 − 748a6 e5 + a4 b(66e1 + 20e6 ))ε 231b3 is negative (resp. positive). Lemma 9. Assume that e2 = φ1 and div(q1 ) = 0, that is, −99b2 e4 + 748a4 e5 − 20a2 be6 e1 = φ2 (e3 , e4 , e5 , e6 ) = . 66a2 b Then the homoclinic loop Lε is stable (resp. unstable) if the value 3b(308b2 e3 − 231a2 be4 − 100a4 e6 ) e5 = φ3 (e3 , e4 , e6 , ε) = − 44a6 is negative (resp. positive). Z ∂Q2 Proof. As proceeded in the proof of Lemma 6 we must study the sign of ε dt. For fixed L1 ∂y values of e1 and e2 according with Lemmas 8 and 9, we get √ √ √ ∂Q2 x(a − bx)2 (2a − bx)(99b2 e4 + 110a4 e5 + 132a3 be5 x + 2a2 b(10e6 − 33e5 x2 )) = √ ∂y 66a2 b5 √ √ bx)2 (a2 + 14a bx − 7bx2 ) e3 (a − + + 3e4 y 2 + 5e6 y 4 . a4 As before, using symmetry and the parametrization of L1 obtained in Subsection 3.2 we obtain Z √ ∂Q2 2(924ab3 e3 − 693a3 b2 e4 + 44a7 e5 − 300a5 be6 ) dt = √ ε. L1 ∂y 3465 b7 The last equation vanishes for e5 = φ3 (e3 , e4 , e6 , ε), so we are done.  Function φ1 has the same role than the function ϕ, but now it controls the existence and stability of the double homoclinic loop. LIMIT CYCLES IN A GENERALIZED RAYLEIGH LIÉNARD OSCILLATOR 15 6. Proof of the main result In this section we prove the main result of the paper. Since Theorem 1 contemplates several cases for different values of s, m and k, we split the proof in some lemmas, see next subsection. Once the lemmas are stated and proved, the proof of the main theorem follows quite straightforwardly. However we emphasize that lemmas are also important by themselves because they explicitly exhibit the approach of choosing suitable perturbations of the parameters to obtain limit cycles, as described in Section 3. This is important because such calculations sometimes are neglected or exhibited for weak focus of low multiplicity as occurs in Hopf Theorem. More precisely, we start considering the values in the parameters space for which the systems present a center. In the following we perform suitable perturbations of these parameters to obtain in each step a change in the sign of the Lyapunov constants. In other words, if for the values of the parameters (α1∗ , . . . , αm ∗ ) the system presents a center, then we start choosing a parameter, that we denote by αm , that determinate the sign of V2k+1 (the last Lyapunov constant), i.e., if αm > α m ∗ , (α < α∗ ) then V m m 2k+1 > 0 (V2k+1 < 0), resp.. Fixing this value of αm , in the following we consider a convenient perturbation of the next parameter, we denote by αm−1 ∗ , that vanishes ∗ the previous Lyapunov constant V2k−1 . We consider a convenient perturbation of αm−1 , denoted by ∗ αm−1 (ε) such that |αm−1 − αm−1 (ε)| << ε and satisfies V2k−1 · V2k+1 < 0 and |V2k−1 | << |V2k+1 |. We repeat this procedure until reach the first Lyapunov constant V1 , where we consider a convenient perturbation such that the equilibrium point be a unstable or stable focus at the system. 6.1. Technical lemmas for the heteroclinic connection occurring in system (6). In this subsection we state and prove some lemmas where we explicitly apply Bautin’s algorithm, see [1, 6], combined with Melnikov techniques. That is actually the very essence of the proof of Theorem 1. We start exhibiting Bautin’s algorithm for obtaining five limit cycles which bifurcate from the center at the origin. Lemma 10. There exists at least five limit cycles bifurcating from the origin for system (6). Proof. In the proof of Lemma 2 under the hypothesis d3 = 0 we obtained that: if the parameters 6bd4 6d3 ε2 16bd6 12a2 bd34 ε2 di ’s satisfy d2 = −3d4 , d1 = − 2 − 5d6 , d5 = 4 − , d 6 =  and d4 = 0, a 5 a2 5 7b2 − 3a4 d24 ε2 then system (6) is Hamiltonian and therefore the origin is a center. Our objective here is to consider small perturbations in these values of parameters to obtain a variation of the signs of the Lyapunov constants and the stability of the origin. First we consider the expression of V11 given in (14). Note that if d4 > 0, d4 < 0 then the origin is a weak repeller focus, weak attractor focus, resp.. Without loss of generality, we assume that d4 > 0. Consider now the expression of V9 given in (13) and the perturbation of the parameter d6 given 12a2 d34 36a6 d5  by d6 (ε) = d62 ε2 + d63 ε3 , with d62 = . Note that |d6 − d6 (ε)| = 245b34 ε4 + O ε6 << ε. 35b Then we obtain the new expression of 21πb2 d63 ε4 27πa2 d54 ε5 9 V9 = − 4 + + πd24 d63 ε6 . 32a 280b 32 16 R. D. EUZÉBIO, J. LLIBRE, D. J. TONON So the sign of V9 is given by the sign of d63 . Then |V9 |  |V11 |, and if d4 > 0 and d63 > 0 then V9 · V11 < 0. In this way we obtain at least one limit cycle bifurcating from the origin under these assumptions. Consider the expression of V7 given in (12), replacing d6 by d6 (ε) = d62 ε2 + d63 ε3 and d5 by 12a2 d34 1  16bd63 d5 (ε) = d52 ε2 + d53 ε3 + d54 ε4 with d62 = , d52 = − 30d34 and d53 = − . We have 35b 7 a2 |d5 − d5 (ε)| = |d54 |ε4 << ε. Therefore we obtain −5 V7 = πd54 ε5 . 64 As in the previous case we obtain |V7 |  |V9 |  |V11 |, and if d4 > 0, d63 > 0 and d54 < 0 then V9 · V11 < 0 and V7 · V9 < 0. So we obtain that at least two limit cycles bifurcating from the origin. Working similarly we consider the expression V5 given in (11), replacing d6 by d6 (ε) = d62 ε2 + d63 ε3 , d5 by d5 (ε) = d52 ε2 + d53 ε3 + d54 ε4 and d1 by d1 (ε) = d10 + d11 ε + d12 ε2 + d13 ε3 + d14 ε4 + d15 ε5 12a2 d34 1  16bd63 6bd4 12a2 d34 with d62 = , d52 = − 30d34 , d53 = − , d10 = − , d11 = 0, d12 = − , d13 = 35b 7 a2 a2 7b 5 −5d63 and d14 = 0, we get |d1 − d1 (ε)| = |d15 |ε << ε and 1 V5 = − πd15 ε6 . 8 So we obtain that |V5 |  |V7 |  |V9 |  |V11 |, and if d4 > 0, d63 > 0, d54 < 0 and d15 > 0 then V9 · V11 < 0, V7 · V9 < 0 and V5 · V7 < 0. So we obtain that at least three limit cycles bifurcating from the origin. Consider now the expression of V3 given in (10), replacing d6 by d6 (ε) = d62 ε2 + d63 ε3 , d5 by d5 (ε) = d52 ε2 + d53 ε3 + d54 ε4 , d1 by d1 (ε) = d10 + d11 ε + d12 ε2 + d13 ε3 + d14 ε4 + d15 ε5 and d2 by 12a2 d34 1  d2 (ε) = d20 + d21 ε + d22 ε2 + d23 ε3 + d24 ε4 + d25 ε5 + d26 ε6 with d62 = , d52 = − 30d34 , d53 = 35b 7 16bd63 6bd4 12a2 d34 − , d10 = − 2 , d11 = 0, d12 = − , d13 = −5d63 , d14 = 0, d20 = −3d4 , d21 = 0, d22 = a2 a 7b 0, d23 = 0, d24 = 0 and d25 = 0, we obtain |d2 − d2 (ε)| = |d26 |ε6 << ε and 1 V3 = − πd26 ε7 . 4 Therefore we get that |V3 |  |V5 |  |V7 |  |V9 |  |V11 |, and if d4 > 0, d63 > 0, d54 < 0, d15 > 0 and d26 < 0 then V9 · V11 < 0, V7 · V9 < 0, V5 · V7 < 0 and V3 · V5 < 0. So we obtain that at least four limit cycles bifurcating from the origin. Finally, replacing d6 by d6 (ε) = d62 ε2 + d63 ε3 , d5 by d5 (ε) = d52 ε2 + d53 ε3 + d54 ε4 , d1 by d1 (ε) = d10 +d11 ε+d12 ε2 +d13 ε3 +d14 ε4 +d15 ε5 , d2 by d2 (ε) = d20 +d21 ε+d22 ε2 +d23 ε3 +d24 ε4 +d25 ε5 +d26 ε6 12a2 d34 1  16bd63 6bd4 and d3 by d3 (ε) = d37 ε7 , with d62 = , d52 = − 30d34 , d53 = − 2 , d10 = − 2 , d11 = 35b 7 a a 12a2 d34 0, d12 = − , d13 = −5d63 , d14 = 0, d20 = −3d4 , d21 = 0, d22 = 0, d23 = 0, d24 = 0 and d25 = 0, 7b we get that |d3 − d3 (ε)| = |d37 |ε7 << ε and the origin is a stable focus if d37 < 0 and unstable focus if d37 > 0. Hence we obtain that |V3 |  |V5 |  |V7 |  |V9 |  |V11 |, and if d4 > 0, d63 > 0, d54 < 0, d15 > 0, d26 < 0 and d37 > 0 then V9 · V11 < 0, V7 · V9 < 0, V5 · V7 < 0 and V3 · V5 < 0. In this way we LIMIT CYCLES IN A GENERALIZED RAYLEIGH LIÉNARD OSCILLATOR 17 obtain one more limit cycle bifurcating from the origin. Therefore we get at least five limit cycles bifurcating from the origin, concluding the proof.  Remark 2. One could argue about the existence of more than five limit cycles inside the heteroclinic loop. Although that situation could be realizable, usually the multiplicity of the weak focus determines an upper bound for the number of limit cycles inside the loop. For instance, in the previous lemma it is easy to see that the first limit cycle bifurcation from the origin, call γ, is unstable. Moreover we are able to compare this information with the persistence or not of the heteroclinic loop. Indeed by fixing the values of parameters for which we obtain five limit cycles, the expression of M1 from Subsection 5.1 we get a8 d34 M1 = √ ε2 + O(ε3 ). 2310 2b4 Since we are assuming d4 > 0 we get that M1 > 0. From Subsection 5.1 the orbit leaving p1 goes away from γ, so by using the Poincaré-Bendixson Theorem a convenient annular region it follows that system (6) has either no limit cycles between γ and the saddle point p1 or it appears in pairs. So we conjecture that five is an upper bound for the number of limit cycles on the region located between the weak focus and the saddle point. Lemma 11. There exist a suitable choice of parameters such that system (6) has simultaneously 3 limit cycles bifurcating from the heteroclinic loop and 2 limit cycles bifurcating from the origin. 12d4 5a4 d5 Proof. As we see in Subsection 5.1 system (6) has a heteroclinic loop Aε1 if d2 = − − + 7 84b2 a2 (99d1 + 160d6 ) 10bd3 40d6 2a2 d5 6d4 b 28d3 b2 + 2 +O(ε). Moreover, if d1 = + − 2 + then the stability 462b a 33 3b a a4 of that loop is established by the sign of d5 = 6b(−100a d6 + 231a bd4 + 308d3 b )/11a6 . The first 4 2 2 part of the proof consists in changing the stability of the loop to get 2 limit cycles, and finally to destroy the loop in order to obtain one more limit cycle. We start assuming d5 < 0 so Aε1 is stable. Now we write d5 (ε) = d50 + εd51 with d50 = (1848b3 d3 + 1386a2 b2 d4 − 600a4 bd6 )/11a6 so Z ∂Q1 a7 d51 (20) M1 = dt = − √ ε2 . A1 ∂y 315 −2b 7 Consequently since a > 0 the loop Aε1 now is defined by the sign of d51 . Indeed, if d51 < 0 then the loop changes its stability from stable to unstable and by applying Poincaré-Bendixson Theorem a stable limit cycle emerges. In order to obtain a new limit cycle we change again the stability of Aε1 . For doing this, consider the ∂Q1 expression of div(p2 ) = considering d1 (ε) a convenient perturbation of the parameter d1 given ∂y by d1 (ε) = d10 +εd11 +ε2 d12 with d11 = 2a2 d51 /3b and d10 = (66b(70bd3 +39a2 d4 )−1160a4 d6 )/33a4 , we get a2 d12 3 (21) div(p2 ) = ε . 7b2 Therefore by choosing d12 < 0 the heteroclinic loop changes stability from unstable to stable and again, by the Poincaré-Bendixson theorem, a limit cycle bifurcates from Aε1 . More precisely, that limit cycles is unstable. Finally to obtain the third limit cycle we choose the parameters so that the loop Aε1 is broken. We do that considering d2 (ε) = d20 + εd21 + ε2 d22 + ε3 d23 where now we choose 18 R. D. EUZÉBIO, J. LLIBRE, D. J. TONON d22 = 3a2 d12 /14b, d21 = a4 d51 /12b2 and d20 = (1980a2 b2 d3 + 495a2 bd4 − 260a4 bd6 )/66a2 b. Then we obtain a4 d23 M 1 = − √ ε3 . 15 2b2 We see that assuming d23 6= 0 the loop is broken. Moreover, from Section 5 and since the second limit cycle is unstable, by choosing d23 > 0 we get M1 < 0. Consequently, by the Poincaré-Bendixson Theorem we obtain a third limit cycle which is stable. Therefore we obtain three limit cycles by using Melnikov method. Moreover, for ε positive and sufficiently small we have that the expression of equations (20) and (21) have opposite sign and |div(p2 )|  |M1 |. It means that the change of stability of Aε1 are only local so we can apply the Poincaré-Bendixson Theorem to obtain the limit cycles. The proof of the lemma follows by applying the same approach used in Lemma 10 having now fixed the values d5 (ε) = d50 + εd51 , d1 (ε) = d10 + εd11 + ε2 d12 and d2 (ε) = d20 + εd21 + ε2 d22 + ε3 d23 . That is, it can be obtained two more limit cycles but now bifurcating from the origin of system (6) so we get 5 limit cycles for such a system. The simultaneity occurs because at each step the obtained limit cycles are hyperbolic so they remain by assuming perturbations of the parameters involving higher orders of ε.  Lemma 12. There exist a suitable choice of parameters such that system (6) has s limit cycles bifur- cating from the origin and m limit cycles bifurcating from the heteroclinic loop with s ∈ {0, 1, 3, 4, 5}, m ∈ {0, 1, 2} and s + m ≤ 5. Proof. The proof of this Lemma 12 is straightforward by using the same construction of Lemmas 10 and 11.  6.2. Technical lemmas for the homoclinic connection occurring in system (7). We now state similar results but concerning system (7). Again we obtain limit cycles bifurcating from the center and from the homoclinic loops considering convenient values of the parameters of system (7). In addition we also consider some symmetry in system (7) so that limit cycles may emerge in pairs. Lemma 13. There exists at least four limit cycles bifurcating from the origin of system (7). 4a2 be1 − 48a4 e5 Proof. In the proof of Lemma 3 under the hypothesis e3 = 0 and e2 = we obtain  b2 2 4 3 b e4 − 16a e5 8a e5 + 3b2 e4 4 that: if the parameters ei satisfies the conditions e1 = − , e6 = , e5 =  4a2 b 5a2 b 4 3 2 2 3 a e4 ε + 8b e4 − and e4 = 0 then the system (7) is Hamiltonian and therefore the origin is a 44a4 center. Similarly we did in the proof of Lemma 10, we will consider small perturbations in these values of parameters to obtain a variation of the sign of the Lyapunov constants and the stability of the origin. The signal of V9 , given in (19), is given by the sign of e4 . Therefore, if e4 > 0, e4 < 0, resp., then the origin is a weak repeller focus, weak attractor focus, resp.. Without loss of generality we assume that e4 > 0. Consider the expression of V7 given in (18) and the perturbation of the parameter e5 LIMIT CYCLES IN A GENERALIZED RAYLEIGH LIÉNARD OSCILLATOR 19 given by e5 (ε) = e50 + e52 ε2 , with 6b2 e4 e50 = − ,   11a4 3e3 we obtain |e5 − e5 (ε)| = − 4 − e52 ε2 << ε and 44 1 3  V7 = πε 3e34 + 44e52 . 32 3 3 So the sign of V7 is given by the sign of 3e34 + 44e52 . Then |V7 |  |V9 | and if e4 > 0 and e52 < − 44 e4 then V7 · V9 < 0. Therefore we conclude that at least one limit cycle bifurcating from the origin with these assumptions. Consider the expression of V5 given in (17) and replacing e6 by e6 (ε) = e60 + e62 ε2 + e63 ε3 and e5 6b2 e4 3be4 8a2 e52 by e50 +e52 ε2 , with e50 = − , e60 = − and e 62 = . We get |e6 −e6 (ε)| = |e63 |ε3 << ε 11a4 11a2 5b and 5 V5 = − πe63 ε4 . 8 3 3 Then |V5 |  |V7 |  |V9 |, and if e4 > 0, e52 < − 44 e4 and e63 < 0 then V7 · V9 < 0 and V5 · V7 < 0. So we obtain that at least two limit cycles bifurcating from the origin. Consider the expression of V3 given in (16), replacing e6 by e6 (ε) = e60 + e62 ε2 + e63 ε3 , e5 by 6b2 e4 3be4 8a2 e52 e5 (ε) = e50 + e52 ε2 and e1 by e1 (ε) = e10 + e14 ε4 , with e50 = − , e60 = − and e 62 = 11a4 11a2 5b 321be4 4 and e10 = − therefore |e 1 − e 1 (ε)| = |e14 |ε << ε and 44a2 πa2 e14 ε5 V3 = − . b 3 3 So we obtain that |V3 |  |V5 |  |V7 |  |V9 | and if e4 > 0, e52 < − 44 e4 , e63 < 0 and e14 > 0 then V7 · V9 < 0, V5 · V7 < 0 and V3 · V5 < 0, then we obtain that at least three limit cycles bifurcating from the origin. Finally, replacing e6 by e6 (ε) = e60 +e62 ε2 +e63 ε3 , e5 by e5 (ε) = e50 +e52 ε2 , e1 by e1 (ε) = e10 +e14 ε4 6b2 e4 3be4 8a2 e52 321be4 and e3 by e3 (ε) = e36 ε6 , with e50 = − , e60 = − and e62 = and e10 = − , we 11a4 11a2 5b 44a2 get that |e3 − e3 (ε)| = |e36 |ε6 << ε and the origin is a stable focus if e36 < 0 and unstable focus if e36 > 0. 3 3 Hence we obtain that |V3 |  |V5 |  |V7 |  |V9 |  |V11 |, and if e4 > 0, e52 < − 44 e4 , e63 < 0, e14 > 0 and e36 > 0 then V7 · V9 < 0, V5 · V7 < 0, V3 · V5 < 0 with V3 < 0 and the origin an unstable focus. In this way we obtain one more limit cycle bifurcating from the origin. Therefore we get at least four limit cycles bifurcating from the origin, concluding the proof.  Lemma 14. There exist a suitable choice of parameters such that system (7) has 4 limit cycles bifurcating from the homoclinic loop and 2 limit cycles bifurcating from the origin. Proof. The proof of the lemma is similar to the proof of Lemma 11, so we only highlight some minor differences between them. Indeed, we perform replacements of the parameters e5 , e1 and e2 in terms of order 1, 2 and 3 in ε, analogously to what we have done for the parameters d5 , d1 and d2 in Lemma 11. Through the replacement of e5 we change the stability of the homoclinic loop Lε1 . 20 R. D. EUZÉBIO, J. LLIBRE, D. J. TONON Applying the Poincaré-Bendixson Theorem in the convenient annular regions, one of them internal to Lε1 and other one external to it, we obtain a limit cycle in each annulus. Proceeding in a complete analogous way we obtain two more limit cycles from the replacement of the parameter e1 . Therefore, we get four limit cycles. In order to obtain the fifth limit cycle, we destroy the loop of Lε1 doing the referred replacement in the parameter e2 . In this case whatever is the sign of M2 (see Section 5.1) we obtain a limit cycle, being internal if M2 < 0 and external to it otherwise. In any case we obtain five limit cycles bifurcating from the homoclinic loop. To obtain the two limit cycle bifurcating from the weak focus we proceed as in Lemma 13.  6.3. Proof of Theorem 1. As commented before, the proof of Theorem 1 consists in apply the lemmas from Subsection 6.1. We remark that in the proof of the lemmas we obtain limit cycles of opposite stability for the cases a > 0, b < 0 and a > 0, b < 0 assuming that d3 = 0 and e3 = 0, respectively. Then one more limit cycle bifurcates by suitable perturbations of d3 and e3 , respec- tively, preserving the already obtained limit cycles. • Proof of statement (a): Assume that a > 0 and b < 0. The configurations (s, 0) with s ∈ {0, 1, . . . , 5} follows from Lemma 10. The case s = 2 and m = 3 follows from Lemma 11. The remaining configurations (s, m), m ∈ {1, 2, 3} and s + m ≤ 5 follows from Lemma 12. The proof of the case a > 0 and b < 0 is then completed since the configurations of limit cycles of systems (5) and (6) are the same. Realization of the maximal ciclicity and simultaneity of limit cycles: We notice that the sharper lower bound of 5 limit cycles when a > 0 and b < 0 can be obtained by the configurations (s, m) given by (5, 0), (4, 1), (3, 2) and (2, 3). The simultaneity is achieved in four out of the five situations. • Proof of statement (b): Now assume that a < 0 and b > 0. We proceed in a similar way to the previous case but now the double homoclinic loop plays a role. In what concerns the configurations of limit cycle bifurcating from the weak focus, after a translation of p1 to the origin, configurations (2s, k) with s ∈ {0, 1, . . . , 4}, k ∈ {0, 1, 2} follows from both Lemma 13 and by symmetry with respect to the y−axis, see Remark 3. Now we prove the other configurations of limit cycles. First, the configuration (2s, 3m), m ∈ {1, 2} is obtained by Lemma 14 because every internal limit cycle bifurcating from the homoclinic loop appears pairwise due to the symmetry. Therefore, each one of the first two steps of the proof of Lemma 14 generates three limit cycles, being one external and two internal, so we get 3m medium limit cycles with m = 1, 2. The 2s small limit cycles of the configuration (2s, 3m) follows also from Lemma 14 and symmetry. Finally the configuration (2s, 3m + k) is obtained from the previous cases and observing that the value of k is determined from the break of the homoclinic loop. That break could generate no limit cycles (k = 0), one external medium limit cycle (k = 1), or two internal medium limit cycles by using symmetry (k = 2), see Lemma 14. Realization of the maximal ciclicity and simultaneity of limit cycles: In this case the maximal ciclicity of 12 limit cycles in obtained in only one situation which corresponds to a configuration of LIMIT CYCLES IN A GENERALIZED RAYLEIGH LIÉNARD OSCILLATOR 21 simultaneity between small and medium limit cycles. Indeed, the configuration 2s + 3m + k = 12 can be obtained with s = m = k = 2, represented by the pair (2, 10). This is precisely the situation of Lemma 14 after applying the symmetry. The cases 2s + 3m + k < 12 are obtained similarly from the other lemmas from the current section. Finally we notice that the configurations of limit cycles after and before the translation of system (7) are clearly preserved. The same results can be obtained for system (5) with a < 0 and b > 0. So we are done. Remark 3. We finish remarking some aspects of Theorem 1 that must be clarified. (I) It does not provide any information about configurations having more than three medium limit cycles. Indeed to obtain such configurations one should take into account higher orders of Melnikov function, which is not considered in this paper. (II) The value k is determined by breaking the loop, then an extra limit cycle can emerge. When a limit cycle takes place, if such a loop is broken in such way that a limit cycle appears internally, then by the symmetry we have a second limit cycle and in this case k = 2. Otherwise, we get an externally limit cycle so in this case k = 1. Acknowledgments. The authors are grateful to the reviewers for their thoughtful comments. The first and third authors are supported by CNPq-Brazil grant 443302/2014-6, PROCAD/CAPES grant 88881.0 68462/2014-01 and PRONEX/FAPEG grant 2017/10267000-508. The second author is partially supported by the Ministerio de Ciencia, Innovación y Universidades, Agencia Estadual de Investigación grant MJM 2016-77278-8 (FEDER), the Agència de Gestió d’Ajuts Universitaris i de Recerca grant 2017SGR1617, and the H2020 European Research Council grant MSCA-RISE- 2017-777911. References [1] Bautin, N. N. On the number of limit cycles which appear with the variation of coefficients from an equilibrium position of focus or centre type, Mat. Sb. vol. 30, 181-196, 1952 (in Russian); Amer. Math. Soc. Transl. vol. 100 (1954). [2] Bejarano, J. D., Garcia-Margallo, L., The limit cycles of the generalized Rayleigh-Liénard oscillator, J. Sound Vib, vol. 156, pp 283–301, 1992. [3] Bejarano, J. D., Garcia-Margallo, L., The greatest number of limit cycles of the generalized Rayleigh-Liénard oscillator, J. Sound Vib, vol. 221, pp 133–142, 1999. [4] Blows, T. R. and Lloyd, N. G. The number of limit cycles of certain polynomial differential equations, Proc. Roy. Soc. Edinburgh Sect. A, (Vol. 98), 1984. [5] Chen, H., Zou, L., Global study of Rayleigh-Duffing oscillators, J. Physics A: Mathematical and Theoretical, vol. 49, pp 165–202, 2016. [6] Dumortier F., Llibre J. and Artés J. C., Qualitative theory of planar differential systems, Springer-Verlag, Univer- sitext, Berlin, 2006. [7] , Gasull, A., Li, W., Llibre, J and Zhang, Z., Chebyshev property of complete elliptic integrals and its application to Abelian integrals, Pacific Journal of Mathematics, vol. 2, 341–361, 2002. [8] Giné, J., Gouveia, L. F. S. and Torregrosa, J., Lower bounds for the local cyclicity for families of centers, J. Differential Equations, vol. 275, 309–331, 2021. 22 R. D. EUZÉBIO, J. LLIBRE, D. J. TONON [9] Gobber F. and Willamowski K.D., Liapunov approach to multiple Hopf bifurcation, J. Math. Anal. Appl., vol. 71, pp. 333–350, 1979. [10] Guckenheimer J. and Holmes P., Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. Revised and corrected reprint of the 1983 original, Applied Mathematical Sciences 42, Springer-Verlag, New York, 1990. [11] Iliev, I. D., On second order bifurcations of limit cycles, J. London Math. Soc, vol. 58, pp. 353–366, 1998. [12] , Liang, H. and Torregrosa, J., Parallelization of the Lyapunov constants and cyclicity for centers of planar polynomial vector fields, J. Differential Equations, vol. 259, 6494–6509, 2015. [13] Leissa, A. W., The historical bases of the Rayleigh and Ritz methods, J. Sound Vib, vol. 287, pp. 961–978, 2005. [14] Lynch, S., Small amplitude limit cycles of the generalized mixed Rayleigh-Liénard oscillator, J. Sound Vib, vol. 178, pp. 615–620, 1994. [15] Petrov, G. S., Elliptic integrals and their nonoscillation, Functional Analysis and Its Applications, 20, 37–40, (1986). [16] Petrov, G. S. Complex zeros of an elliptic integral, Functional Analysis and Its Applications, 23, 160–161, (1989). [17] Poincaré, H. Mémoire sur les courbes définies par une équation différentielle, Journal de Mathématiques Pures et Appliquées, vol. 8, pp 251-296, 1882. [18] Lord Rayleigh, The Theory of Sound, vol. 1, The Macmillan Company, 1877. [19] Roussarie, R., On the number of limit cycles which appear by perturbation of separatrix loop of planar vector fields, Bol. Soc. Brasil Mat., vol. 17, pp. 67–101, 1986. [20] Roussarie, R., Bifurcations of planar vector fields and Hilbert’s sixteenth problem, reprint of the 1998 edition. Modern Birkhäuser Classics. Birkhäuser/Springer, Basel, 2013. [21] Van der Pol, B., Forced oscillations in a circuit with nonlinear resistance, London, Edinburgh and Dublin Phil. Mag., 3, 1927. [22] Wu Y., Han M. and Chen X., On the study of limit cycles of the generalized Rayleigh-Liénard oscillator, Intern. J. Bifurc. and Chaos in Appl. Sci. and Eng., vol. 14, 2906–2914, 2004. 1 Institute of Mathematics and Statistics of Federal University of Goiás, Avenida Esperança s/n, Campus Samambaia, 74690-900, Goiânia, Goiás, Brazil 2 Departament de Matemàtiques, Universitat Autònoma de Barcelona, 08193 Bellaterra, Catalunya- Barcelona, Spain. Email address: euzebio@ufg.br Email address: jllibre@mat.uab.cat Email address: djtonon@ufg.br