Available online at www.sciencedirect.com ScienceDirect Journal of Differential Equations 309 (2022) 291–310 www.elsevier.com/locate/jde First-order perturbation for multi-parameter center families Jackson Itikawa a , Regilene Oliveira b , Joan Torregrosa c,d,∗ a Departament of Mathematics, Universidade Federal de Rondônia, Porto Velho, RO 76801–0590, Brazil b Departamento de Matemática, ICMC-Universidade de São Paulo, Avenida Trabalhador Sãocarlense, 400, São Carlos, SP, 13566–590, Brazil c Departament de Matemàtiques, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Spain d Centre de Recerca Matemàtica, Campus de Bellaterra, 08193 Bellaterra, Barcelona, Spain Received 9 April 2021; revised 9 November 2021; accepted 23 November 2021 Available online 30 November 2021 Abstract In the weak 16th Hilbert problem, the Poincaré-Pontryagin-Melnikov function, M1 (h), is used for ob- taining isolated periodic orbits bifurcating from centers up to a first-order analysis. This problem becomes more difficult when a family of centers is considered. In this work we provide a compact expression for the first-order Taylor series of the function M1 (h, a) with respect to a, being a the multi-parameter in the unperturbed center family. More concretely, when the center family has an explicit first integral or inverse integrating factor depending on a. We use this new bifurcation mechanism to increase the number of limit cycles appearing up to a first-order analysis without the difficulties that higher-order studies present. We show its effectiveness by applying it to some classical examples. © 2021 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). MSC: primary 34C07, 34C29; secondary 34C25, 37G15 Keywords: Limit cycle; Melnikov functions; Averaging theory; Multi-parameter perturbation * Corresponding author. E-mail addresses: itikawa@unir.br (J. Itikawa), regilene@icmc.usp.br (R. Oliveira), torre@mat.uab.cat (J. Torregrosa). https://doi.org/10.1016/j.jde.2021.11.035 0022-0396/© 2021 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). J. Itikawa, R. Oliveira and J. Torregrosa Journal of Differential Equations 309 (2022) 291–310 1. Introduction and statement of the main results A limit cycle is an isolated periodic orbit in the set of all periodic orbits of a differential system. The maximal number of limit cycles that a polynomial differential system in R2 of degree n might have is denoted by H(n) and it is called the 16th Hilbert’s problem. For further information see [8–10] and the references therein. The determination of H(n), even for the simplest case n = 2, remains open. A classic strategy to investigate lower bounds for H(n) is to apply the averaging theory of first and/or higher-order in the study of the number of limit cycles that bifurcate from centers. In fact, the problem is reduced to the study of the number of positive simple zeros of a specific non-linear function. For autonomous planar vector fields, it coincides with the well- known Poincaré-Pontryagin-Melnikov function, see more details in [2]. When we perturb an integrable system, this function  Q(x, y) dx − P (x, y) dy M1 (h) = (1) V (x, y) H (x,y)=h is obtained from the Taylor expansion of the displacement function (h, ε) = εM1 (h) + ε 2 M2 (h) + · · · , (2) associated to the differential equation ⎧ ⎪ P (x, y) ⎪  ⎨ x = −Hy + ε , V (x, y) (3) ⎪ ⎪ Q(x, y) ⎩ y  = Hx + ε . V (x, y) Here we have denoted by V the inverse integrating factor defined in an open set of R2 that contains the singular point which is a center for the unperturbed system and, as usual, the level sets of the first integral, H (x, y) = h, are closed curves for h ∈ [0, h0 ). The 0 and h0 level sets of H define the center equilibrium point and the boundary of the period annulus, respectively. We remark that the function (1) is called an Abelian integral when V is a constant and the functions H, P , Q are polynomials. In the literature, the investigation of the number of zeros of the function M1 is also called the weakened 16th Hilbert’s problem and it was proposed by Arnold, see [1]. The Implicit Function Theorem provides the relation of both problems. More concretely, for each positive simple zero h∗ of M1 there exists a limit cycle ε of (3) such that ε tends to the level curve {H = h∗ } when ε goes to zero. Usually, as the center defined by the first integral H and the inverse integrating factor V are fixed, the function M1 (h), defined in (1), depends only on the perturbation parameters defined in the perturbative functions P and Q. Clearly, when all the involved functions are polynomials, this dependence is linear in the coefficients of P and Q. In this paper we will perturb families of systems depending on a multi-parameter a = (a1 , a2 , . . . , a ) ∈ R . That is, the associated first integral and inverse integrating factor are respectively H = H (x, y, a) and V = V (x, y, a) in (1) and (3). So, the displacement function (2) is given by (h, a, ε) = εM1 (h, a) + ε 2 M2 (h, a) + · · · . (4) 292 J. Itikawa, R. Oliveira and J. Torregrosa Journal of Differential Equations 309 (2022) 291–310 The main obstruction on the study of the zeros of  Q(x, y) dx − P (x, y) dy M1 (h, a) = (5) V (x, y, a) H (x,y,a)=h is the difficulty to work with the above integral for any a. Most of the works study special values of a for which the integrals on the level curves of H (x, y, 0) = h are explicitly obtained although they were defined via elliptic or hyperelliptic integrals. When a family is analyzed, it is not restrictive to put these special cases at a = 0. Then, near this special point, the study of the zeros of M1 can be developed via the Taylor series of M1 with respect to a parameter. This idea has been used previously in some examples, see for example [7, Chapter 9]. The goal of our main result is to get a closed formula for its linear Taylor development. We will show how, with this mechanism, we obtain more limit cycles than working only at a = 0. Obviously, a second-order approach might provide better lower bounds for H(n) but with a higher computational effort. But it is well-known that a higher-order analysis not always can be developed. As we will see in the proofs, although this bifurcation technique is essentially a first-order bifurcation analysis for a family of centers, it can be partially interpreted as a second-order analysis of a fixed center. We could say that we are seeing a piece of the second-order function M2 in (2). Hence, roughly speaking we could say that we are presenting an averaging function of order one and a half. In the following, we present our main results, which will be proved in Section 2. Firstly, we detail the expression of the first-order Taylor series of M1 (h, a) with respect to a. Secondly, we define the new Poincaré–Pontryagin–Melnikov parametric function and how its positive simple zeros provide limit cycles for system (3). Theorem 1.1. For a ∈ R in a neighborhood of the origin, we can write (5) as   M1 (h, a) = M0 (h) − ai (Mi (h) + Li (h)) + O(a2 ), (6) i=1 where  Q(x, y) dx − P (x, y) dy M0 (h) = , V (x, y, 0) H (x,y,0)=h   ∂V (x, y, a)  Q(x, y) dx − P (x, y) dy Mi (h) =  , ∂ai a=0 (V (x, y, 0))2 H (x,y,0)=h   Q(x, y) dx − P (x, y) dy  Li (h) = LXHa  , i V (x, y, a)  a=0 H (x,y,0)=h 1 ∂H with XHai = ∇H . Here ||∇H || ∂ai 2 293 J. Itikawa, R. Oliveira and J. Torregrosa Journal of Differential Equations 309 (2022) 291–310 ∂α ∂α ∂f ∂g ∂β ∂β ∂f ∂g LX ω = f +g +α +β dx + f +g +α +β dy ∂x ∂y ∂x ∂x ∂x ∂y ∂y ∂y ∂ is the Lie derivative of the one-form ω = α dx + β dy with respect to the vector field X = f + ∂x ∂ g . ∂y The expression for M1 (h, a) defined above is a consequence of the proof of Kelvin’s circu- lation theorem, for more details see [5]. In the following we present a direct application of this first-order development study which proposes a new first-order Melnikov function, providing a better understanding of how the number of limit cycles bifurcates from a period annulus change depending on the parameters of a center family. As we have commented before, we are interested in periodic orbits in planar polynomial differential equations. In this case, the perturbative functions P (x, y, λ) and Q(x, y, λ) de- pend linearly on the monomial coefficients λ. Hence, also the function M1(h) defined in (1) depends linearly on λ. Under this assumption, the next result present a new Poincaré-Pontryagin- Melnikov function that can be used in the study of limit cycles bifurcating from center families. Theorem 1.2. Let M1 (h, a, λ) be the function defined in (6) associated to the perturbed differ- ential equations (3) with H = H (x, y, a), V = V (x, y, a) and being λ = (λ0 , λ1 ) ∈ Rm0 ⊕ Rm1 the vector of coefficients of the perturbative functions P (x, y, λ), Q(x, y, λ). We assume that [j ] [j ] [j ] m0 [j ] [j ] there exist functions f0 , f1,i , and g0 , such that M0 (h, λ) = λ0 f0 (h) and Mi (h, λ) + j =1  m0 [j ] [j ]  m1 [j ] [j ] Li (h, λ) = λ0 gi (h) + λ1 f1,i (h), for i = 1, . . . , . Then, taking λ0 = δμ0 , λ1 = μ1 , j =1 j =1 and a = δb, for each positive simple zero h∗ of  m0 [j ] [j ]    m1 [j ] [j ] N[1] (h, μ) = μ0 f0 (h) + bi μ1 f1,i (h) (7) j =1 i=1 j =1 and for a, ε small enough, there exists a limit cycle ε,h∗ of the perturbed system (3) such that ε,h∗ tends to the closed curve {H = h∗ } when (h, a) goes to (h∗ , 0). We remark that all the involved function M0 and Mi , Li , for i = 1, . . . ,  in (6) are defined via integrals over the closed level curves of a non-perturbed system (i.e. (3) with ε = 0). More- over, the above result can be used only in the bifurcation study of limit cycles of small (Hopf) and medium (Melnikov) amplitudes. We can not consider the bifurcation of limit cycles from the outer boundary of the period annulus, because it could change with a. This last problem requires a more delicate analysis, similar to the one developed recently in [14]. The key point in Theo- rem 1.2 is the blow-up type change of parameters. This mechanism is also used in [4] to analyze the limit cycles bifurcating from the outer boundary of center when it is of heteroclinic cycle type. In [4] the number of zeros of the Abelian integral is less than the ones obtained studying directly the return map. This bifurcation phenomenon is similar to the one described in Theo- rem 1.2, because the Abelian integral M0 obtained for a = 0 changes to (7) adding, taking a 294 J. Itikawa, R. Oliveira and J. Torregrosa Journal of Differential Equations 309 (2022) 291–310 small enough, some extra terms that increases the number of limit cycles in some families. This blow-up type change is also used recently in [6] for studying the local cyclicity problem when center polynomial families are perturbed also with polynomials but with a fixed small degree class. We notice that in (7) we have an extra dependence on b that not always is used in the study of the number of zeros of (7). Clearly, this number increases when some of the functions [j ] f1,i are independent with respect to the set {f0[1] , . . . , fm[1]0 }. This is the case in all the analyzed applications. Some applications of the above theorem are presented in Section 3. In such section, we also discuss how the presence of a parameter a in the system can increase the number of limit cycles on it, i.e., why the function (6) has more zeros with a = 0 than when a = 0. We notice that, as can be seen in the definition of the function N[1] , the dependency with respect to the parameters (a, λ) is actually of degree two. We could say that we are doing a partial second-order study. In general, we are not seeing completely the second-order Melnikov function but we are looking at part of it. As we are particularly interested in planar polynomial vector fields, we have presented the results considering that the first integral H is of Darboux type and, as it is usual in these cases, the inverse integrating factor as well as the perturbation functions P and Q are polynomials. Clearly we only need a C 1 dependence on the multi-parameter a in the unperturbed vector field. Hence, our last application of first-order analysis deals with piecewise polynomial perturbations of a classical and simple Hamiltonian. The first-order analysis obtained from Theorems 1.1 and 1.2, using the function N[1] , is mainly based on the fact that we have a closed formula for the first derivative with respect to the parameters in the center family. A generalization of Theorem 1.2 is presented in Section 4 when higher-order developments of (5) can be obtained. The blow-up procedure is used up to kth-order developments for special center families such that only V depends on a, while H not. The higher-order analysis is studied from a collection of functions N[i] , for i = 1, . . . , k, presented in Section 4. To be consistent in the notation we notice that M0 = N[0] . 2. First averaging function depending on parameters In this section we prove our main results Theorems 1.1 and 1.2. We remark that if system (3) has a parameter a then the Melnikov function of first-order M1 (h, a), defined in (5), depends on the first-order Taylor series of such function with respect to the multi-parameter a ∈ R at a = 0. In Theorem 1.1 we present a computable expression for the first-order Taylor series in a of the first-order Melnikov function M1 (h, a). Since the expression of the function M1 (h, a) involves an integration of a one-form over a closed curve parametrized by a, in order to compute the partial derivatives with respect to the components of a, we need a multivariate version of the Leibniz integral rule. The next technical lemma provides the derivative of the integral over a closed curve of a one-form when both depend on a parameter that, for simplicity, we will take only one dimensional. For the sake of completeness, we provide its proof. For further information on this topic, see for instance [5]. Lemma 2.1. Consider b ∈ R. Let ω = ω(x, y, b) be a one-form, γb = γ (x, y, b) a closed curve, and X the vector field associated to the flow h that describes the variation of γb with respect to b, that is γb+h = h (γb ) for all h ≈ 0. Then    d ∂ω ω= + LX ω, (8) db ∂b γb γb γb 295 J. Itikawa, R. Oliveira and J. Torregrosa Journal of Differential Equations 309 (2022) 291–310 where LX ω is the Lie derivative of the one-form ω with respect to the vector field X . Proof. From the definition of the first derivative with respect to b, the left hand side of (8) can be written as ⎛ ⎞    d 1⎝ ω = lim ω(x, y, b + h) − ω(x, y, b)⎠ db h→0 h γb γb+h γb ⎛ ⎞   1⎜ ⎟ = lim ⎝ ω(x, y, b + h) − ω(x, y, b)⎠ h→0 h h (γb ) γb ⎛ ⎞   1⎝ = lim ∗ h ω(x, y, b + h) − ω(x, y, b)⎠ , h→0 h γb γb where ∗h ω(x, y, b + h) denotes the pullback of the one-form ω along h , so the last equality above is the substitution rule for integrals in several variables. This expression can be computed as follows: first we swap the limit and the integral; second we  ∂ replace ω(x, y,b + h) by 2its Taylor series at h = 0 which is ω(x, y, b + h) = ω(x, y, b) + b ∂h ω(x, y, b + h) h=0 + O(h ); and third we apply the linearity of ∗ to obtain the following    1 ∗ ∗ ∂ ∗ h w(x, y, b) + h w(x, y, b + h) + h O(h ) − ω(x, y, b) 2 lim h . h→0 h ∂h h=0 γb ∗ h O(h ) = 0 2 Since lim and h→0   ∗ ∂ ∂ ∂ lim h ω(x, y, b + h) = lim ω(x, y, b + h) = ω(x, y, b + h) , h→0 ∂h h→0 ∂h ∂h h=0 the last equality becomes     ∗ ω(x, y, b + h) − ω(x, y, b) ∂ ω(x, y, b + h) + lim h . ∂h h=0 h→0 h γb γb The proof finishes because the limit in the integral is by definition the Lie derivative of the one- form ω with respect to the vector field X , that is LX ω.  Proof of Theorem 1.1. Clearly the evaluation at a = 0 of M1 (h, a) provides (6), that is M0 (h) = M1 (h, 0). For the computation of the first derivative with respect to each component of the multi- parameter a = (a1 , a2 , . . . , a ) we can apply (8) taking b = ai , for i = 1, . . . , . The curve γa is defined implicitly by H (x, y, a) = h and the one-form ω by 296 J. Itikawa, R. Oliveira and J. Torregrosa Journal of Differential Equations 309 (2022) 291–310 Q(x, y) dx − P (x, y) dy ω(x, y, a) = . V (x, y, a) The derivative of the level curve H (x, y, a) = h with respect to time provides the relation ∇H · ∂H XHai + = 0, where ∇ denotes the gradient operator. Since XHai ⊥ H , then we can write ∂ai XHai = F ∇H, ∂H for a function F . Hence, as ∇H · F ∇H + = 0, we can get the function F that writes as ∂ai −1 ∂H F= . This is the provided expression of the vector field −XHai in the statement. ||∇H ||2 ∂ai The proof finishes substituting all the above expressions in (8) and considering the linearity of the Lie derivative we obtain the expression for the first derivative with respect to a, that is   − ai (Mi (h) + Li (h)). The functions Mi (h) and Li (h) are the ones provided in the state- i=1 ment.  Proof of Theorem 1.2. From the hypotheses of the statement and equation (6), the displacement function (h, a, λ, ε) defined in (4) writes as ⎛ ⎛ ⎞ ⎞ m0 [j ] [j ]   m0 [j ] [j ]  m1 [j ] [j ] ε⎝ λ0 f0 (h) + ai ⎝ λ0 gi (h) + λ1 f1,i (h)⎠ + O(a2 )⎠ + O(ε 2 ). (9) j =1 i=1 j =1 j =1 With the blow-up type change in the parameter space, writing λ0 = δ λ̃0 , a = δ ã, and ε = δ 2 . Then (9) writes as ⎛ ⎛ ⎞ ⎞ m0 [j ] [j ]   m0 [j ] [j ]  m1 [j ] [j ] δ2 ⎝ δ λ̃ f (h) + 0 δ ãi ⎝ 0 δλ g (h) + λ f (h)⎠ + O(δ 2 )⎠ + O(δ 4 ) 0 i 1 1,i j =1 i=1 j =1 j =1 ⎛ ⎞ m0 [j ] [j ]    m1 [j ] [j ] =δ 3 ⎝ λ̃ f (h) + 0 ãi 0 λ f (h)⎠ + O(δ 4 ). 1 1,i j =1 i=1 j =1 The coefficient of δ 3 in the above function defines the function N[1] in the statement where, for simplicity, the ˜ in the parameters λ0 and a has been removed. The statement follows directly from the Implicit Function Theorem dividing the above function by δ 3 , as in the classical result the displacement function is divided by ε done by Pontryagin in [16].  3. More limit cycles from developments of first-order In this section, we present four different applications of Theorems 1.1 and 1.2. The aim is to study the existence of limit cycles emphasizing how the number of zeros of the Abelian inte- gral for a = 0 changes when a is a non-vanishing but small multi-parameter. More concretely, showing which is the influence of this distinguished parameter. The first and simpler example is 297 J. Itikawa, R. Oliveira and J. Torregrosa Journal of Differential Equations 309 (2022) 291–310 the perturbation of the quadratic potential with quadratic polynomials getting, as usual, only one limit cycle, see Proposition 3.1. We remark that, up to a first-order analysis with a = 0, no limit cycles bifurcate. The second center family has three parameters, it is the quartic potential, which we perturb inside the quartic polynomial family, see Proposition 3.2. These first two examples are Hamiltonian, so, as they have constant inverse integrating factors, the functions Mi vanish for i = 0. The relevance of our main results is shown in Proposition 3.3 where we perturb a center family having a rational first integral and, obviously, an inverse integrating factor. Both depend- ing on a, which is chosen in R2 to detail better how Theorems 1.1 and 1.2 apply. Although we have not detailed either in the statements or in the proofs, both main results can be applied when we do a piecewise perturbation of a smooth family of centers. The only difference is how we compute the integrals defined in (6), see for example [11]. This last application is described in Proposition 3.4. Proposition 3.1. Consider the perturbed differential system (x  , y  ) = (−y, x + ax 2 ) + ε(P (x, y, λ), Q(x, y, λ)), (10) where P and Q are quadratic polynomials in (x, y), a ∈ R, and ε = 0 sufficiently small. Then, for each small enough a = 0, there exist λ such that system (10) has a limit cycle bifurcating from the period annulus up to a first-order analysis in ε. Moreover, when a = 0 no limit cycles bifurcate from the period annulus up to a first-order analysis in ε.  2  2 Proof. We consider system (10) with P (x, y, λ) = pij x i y j , Q(x, y, λ) = qij x i y j . i+j =0 i+j =0 As it is not restrictive to assume that P (0, 0, λ) = Q(0, 0, λ) = 0, the parameter space defined by the coefficients of the perturbative polynomials P and Q, is denoted by λ ∈ R10 . We notice that for ε = 0 system (10) has the following Hamiltonian x2 + y2 x3 H (x, y, a) = +a , 2 3 taking V (x, y, a) = 1 as a trivial inverse integrating factor. In order to calculate the function M1 (h, a), defined in (6), we compute the following vector field x 4 (ax + 1) x3y X Ha = ∂x + ∂y 3(a 2 x 4 + 2ax + x + y ) 3 2 2 3(a x + 2ax 3 + x 2 + y 2 ) 2 4 and we obtain the expression of LXHa (P (x, y) dy − Q(x, y) dx), which we omit here because of its size. Applying a polar change of coordinates x = h cos θ , y = h sin θ , with h ∈ (0, +∞), we obtain the following expressions M0 (h) = −π(p10 + q01 )h2 , M1 (h) = 0, since V (x, y, a) = 1, 298 J. Itikawa, R. Oliveira and J. Torregrosa Journal of Differential Equations 309 (2022) 291–310 p q11  4 20 L1 (h) = − π + h . 2 4 Clearly, the function M0 has no positive simple zeros and the second part of the statement follows. Considering all the functions we have p q11  4 20 M1 (h, a) = −π(p10 + q01 )h2 + aπ + h + O(a 2 ). 2 4 The proof follows directly applying Theorem 1.2. The function (7) writes as N[1] (h, λ) = h2 (λ0 + λ1 h2 ), being m0 = m1 = 1, b = 1, λ0 = −π(p10 + q01 ), λ1 = π(2p20 + q11 )/4, f0 (h) = h2 , g0 (h) = 0, and f1 (h) = h4 . The superscript [1] has been removed to simplify the reading.  We observe that from the differential system (10) when a = 0 (linear center) no limit cycle bifurcates up to a first-order analysis in the quadratic polynomial class. Therefore, using the function presented in Theorem 1.1, which is written so that Theorem 1.2 applies, we have a limit cycle due to the special role that the small additional parameter a = 0 has. Proposition 3.2. Consider the perturbed differential system (x  , y  ) = (−y, x + a1 x 2 + a2 x 3 + a3 x 4 ) + ε(P (x, y, λ), Q(x, y, λ)), (11) where P and Q are quartic polynomials, a = (a1 , a2 , a3 ) ∈ R3 , and ε = 0 sufficiently small. Then, for small enough a = 0, there exist λ such that system (11) has at least three limit cycles bifurcating from the period annulus up to a first-order analysis in ε. Moreover, when a = 0 only one limit cycle bifurcates from the period annulus up to a first-order analysis in ε. Proof. Note that 1 a1 a2 a3 H (x, y, a1 , a2 , a3 ) = (x 2 + y 2 ) + x 3 + x 4 + x 5 2 3 4 5 is a first integral of system (11) and V (x, y, a1 , a2 ) = 1 is an inverse integrating factor, so such system has a non-degenerate center at the origin. Applying polar coordinates x = h cos θ, y = h sin θ and performing the calculations to get the integrals described in Theorem 1.1, we obtain 1 M0 (h) = − π (q21 + 3q03 + 3p30 + p12 ) h4 − π (q01 + p10 ) h2 , 4 1 1 L1 (h) = − π (5q31 + 3q13 + 20p40 + 2p22 ) h6 − π (q11 + 2p20 ) h4 , 24 4 1 3 L2 (h) = − π (5q21 + 3q03 + 15p30 + p12 ) h6 − π (q01 + p10 ) h4 , 32 16 1 1 L3 (h) = − π (7q31 + 3q13 + 28p40 + 2p22 ) h8 − π (q11 + 2p20 ) h6 . 64 8 299 J. Itikawa, R. Oliveira and J. Torregrosa Journal of Differential Equations 309 (2022) 291–310 Since V (x, y, a) = 1 we have that Mi (h) = 0, i = 1, 2, 3. It is easy to check that the coefficients of M0 are linearly independent and it has only one positive simple zero, therefore, the second part of the statement follows. After performing a linear change of variables in the parameter space, we obtain, by using the notation introduced in Theorem 1.2,   M1 (h, a) = λ[2] 0 h4 + λ [1] 2 0 h + a 1 λ [3] 6 1 h + λ [1] 4 1 h + λ[2] 3λ[1] 3λ[3] λ[1] a2 λ[2] 1 − 0 h6 − 0 h4 + a3 λ[4] 1 + 1 h8 + 1 h6 8 16 8 2 and, reordering monomials in h and taking b1 = b2 = b3 = 1, 3μ[3] μ[1]  [1]  4 N[1] (h, μ) = μ[4] 1 + 1 h + 8 1 + μ[2] [3] 1 + μ1 h6 + μ[2] [1] 2 0 + μ1 h + μ0 h . 8 2 The proof finishes taking for example μ[1] [3] 1 = μ1 = 0 and using Theorem 1.2. Because the remaining four parameters are independent and therefore the polynomial N[1] will have three positive simple zeros. To simplify the reading, we have not written explicitly the functions [j ] [j ] [j ] f0 , gi , and f0 .  In the last proof, as it can be seen analyzing the function M1 , the number of zeros increases from one, when a1 = a2 = a3 = 0, to three when a1 a2 a3 = 0. But, clearly, the same result can be obtained with a1 = a2 = 0 and a3 = 0. Proposition 3.3. Consider the perturbed differential system    x  = − (−1 + 2 x) 2 a2 x 2 y + 6 a2 y 3 + a1 x 2 + 5 a1 y 2 + 4 y + εP (x, y, λ),    (12) y  = 4 x 2 − y 2 − x a2 y 2 + a1 y + 1 + εQ(x, y, λ), where P and Q are quartic polynomials, a = (a1 , a2 ) ∈ R2 , and ε = 0 sufficiently small. Then, for small enough a = 0, there exist λ such that system (12) has at least four limit cycles bifur- cating from the period annulus up to a first-order analysis in ε. Moreover, when a = 0 only two limit cycles bifurcate from the period annulus up to a first-order analysis in ε.  4  4 Proof. We consider system (12) with P (x, y, λ) = pij x i y j and Q(x, y, λ) = qij x i y j . i+j =0 i+j =0 As it is not restrictive to assume that P (0, 0, λ) = Q(0, 0, λ) = 0, the vector defined by the coefficients pij and qij provides λ ∈ R28 . We notice that for ε = 0 system (12) has the following first integral  (x 2 + y 2 ) a2 y 2 + a1 y + 1 H (x, y, a) = 1 − 2x  and the inverse integrating factor V (x, y, a) = 2 (1 − 2 x)2 a2 y 2 + a1 y + 1. 300 J. Itikawa, R. Oliveira and J. Torregrosa Journal of Differential Equations 309 (2022) 291–310 Since the first integral H is analytic in a neighborhood of the origin, by the well-known Lyapunov–Poincaré Theorem [15] system (12) with ε = 0 has a non-degenerate center at the origin. The vector field XHai , for i = 1, 2, defined in the statement of Theorem 1.1 is given by XHa2 = yXHa1 , where  XHa1 = S(x, y) 4(x 2 − y 2 − x)(a2 y 2 + a1 y + 1)∂/∂x  − (1 − 2x)(2a2 x 2 y + 6a2 y 3 + a1 x 2 + 5a1 y 2 + 4y)∂/∂y , with  S(x, y) = y(1 − 2x) (4x 4 + 52x 2 y 2 + 16y 4 − 4x 3 − 68xy 2 + x 2 + 25y 2 )a12 + 4y(4x 4 + 36x 2 y 2 + 8y 4 − 4x 3 − 44xy 2 + x 2 + 15y 2 )a1 a2 + 8y(8x 2 + 4y 2 − 12x + 5)a1 + 4y 2 (4x 4 + 24x 2 y 2 + 4y 4 − 4x 3 − 28xy 2 + x 2 + 9y 2 )a22 −1 + 16y 2 (6x 2 + 2y 2 − 8x + 3)a2 + 16(x 2 + y 2 − 2x + 1) . We apply the following change of coordinates, defined for h ∈ [0, 1), h cos θ − h2 h sin θ x= , y= , 2h cos θ − h2 − 1 2h cos θ − h2 − 1 to the integrals (6). By direct computation we have that Mr (h), r = 0, 1, 2 and Ls (h), s = 1, 2 are written as sum of Abelian integrals multiplied by rational functions in h. More concretely, each of them writes as  Rk,m k,m j (h, λ)Ij (h), where the sum is defined for j = 1, . . . , 7, m = 0, 1, and k = 0, . . . , j − m + 1. Moreover, Rk,m j (h, λ) are rational functions on h depending linearly on the perturbation parameters λ and the Abelian integrals 2π cosk θ sinm θ dθ Ijk,m (h) = (2h cos θ − h2 − 1)j 0 can be explicitly computed. Therefore the expression of M1 (h, a) defined in (6) is given as fol- lows M1 (h, a1 , a2 ) = M0 (h) + a1 (M1 (h) + L1 (h)) + a2 (M2 (h) + L2 (h)) + O(a2 ), with 301 J. Itikawa, R. Oliveira and J. Torregrosa Journal of Differential Equations 309 (2022) 291–310 h2 P2,1 (h2 ) h4 P2,2 (h2 ) h4 P3,1 (h2 ) M0 (h) = , M1 (h) = , M2 (h) = , (h2 − 1)3 (h2 − 1)4 (h2 − 1)5 h4 P2,3 (h2 ) h4 P3,2 (h2 ) L1 (h) = , L2 (h) = , (h2 − 1)4 (h2 − 1)5 being Pi,j polynomials of degree i. Moreover, their coefficients are polynomials of degree 1 with rational coefficients in the parameters λ. Applying a suitable linear change of variables in the parameter space, we get h2 (λ[3] [2] 2 [1] 0 h + λ0 h + λ0 ) 4 h4 (λ[3] [2] 2 [1] 1 h + λ1 h + λ1 ) 4 M1 (h, a) =  3 + a1 h2 − 1 (h2 − 1)4 h4 (λ[4] [3] 4 [2] 2 [1] 2 h + λ2 h + λ2 h + λ2 ) 6 + a2 , (h2 − 1)5 that, with the adequate rescaling and taking b1 = b2 = 1, the function (7) is h2 N[1] (h, μ) = ((μ[4] [3] [3] 8 [3] [3] [2] [2] [3] 6 2 + μ1 + μ0 )h + (μ2 − μ1 + μ1 + μ0 − 2μ0 )h + (h2 − 1)5 (μ[2] [2] [2] [3] [2] [1] 4 [1] [1] [2] [1] 2 [1] 2 − μ1 + μ1 + μ0 − 2μ0 + μ0 )h + (μ2 − μ1 + μ0 − 2μ0 )h + μ0 ). The proof of the first part of the statement finishes as the previous proofs just taking zero all the parameters except μ[1] [2] [3] [3] [4] 0 , μ0 , μ0 , μ1 , and μ2 . We notice that the remaining five coefficients in N[1] are linearly independent, so the four limit cycles described in the statement can be obtained easily. Finally, the number of positive simple zeros of M1 taking a1 = a2 = 0 is only two and, consequently, the first part of the statement also follows.  In the above proof we could choose a rational first integral and a rational inverse integrating factor for system (12) with ε = 0, but the way that we have presented shows that our approach also works when the first integral is non-rational. Proposition 3.4. Consider the piecewise perturbed differential system    (P + (x, y, λ), Q+ (x, y, λ)), if y > 0, (x , y ) = (−y, ax + ax + x) + ε 3 2 (13) (P − (x, y, λ), Q− (x, y, λ)), if y < 0, where P + , P − , Q+ , and Q− are cubic polynomials, a ∈ R and ε = 0 sufficiently small. Then, for small enough a = 0, there exist λ such that system (13) has at least four limit cycles bifurcating from the period annulus up to a first-order analysis in ε. Moreover, when a = 0 only two limit cycles bifurcate from the period annulus up to a first-order analysis in ε. Proof. System (13), when ε = 0, has 1 a a H (x, y, a) = (x 2 + y 2 ) + x 3 + x 4 2 3 4 302 J. Itikawa, R. Oliveira and J. Torregrosa Journal of Differential Equations 309 (2022) 291–310 as a first integral, being V (x, y, a) = 1 the corresponding inverse integrating factor. We note that system (13) with ε = 0 has a non-degenerate center at the origin for any a. In this case λ ∈ R40 and it is not restrictive to assume that P ± (0, 0, λ) = Q± (0, 0, λ) = 0. Doing the polar change of coordinates x = h cos θ, y = h sin θ we perform the calculations similar to the ones in the previous examples, except to the integrals along the closed curves are done in two parts, the upper half-plane (θ ∈ [0, π]) and the lower half-plane (θ ∈ [π, 2π]) since system (13) has the straight line y = 0 of discontinuity. Then we get 1 + − + − + − + − M0 (h) = − (3(q03 + q03 ) + (p12 + p12 ) + (q21 + q21 ) + 3(p30 + p30 ))h4 8 2 − + − + − + 1 + − + − + ((p11 − p11 ) + (q20 − q20 ) + 2(q02 − q02 ))h3 − ((p1,0 + p1,0 ) + (q0,1 + q0,1 ))h2 , 3π 2 1 + − + − + − + + L1 (h) = (15(p30 + p30 ) + 5(q21 + q21 ) + (p12 + p12 ) + 3(q03 + q03 ))h6 64 1 + − + − + − + − + (8(q12 − q12 ) + 20(q30 − q30 ) + 8(p21 − p21 ) + 15(q20 − q20 ) 30π + − + − 1 + − + − + 3(p11 − p11 ) + 6(q02 − q02 ))h5 + ((4q11 + q11 ) + 3(p10 + p10 ) 32 + − + − 2 + − 3 + 8(p20 + p20 ) + 3(q01 + q01 ))h4 + (q − q10 )h . 3π 10 Since V (x, y, a) = 1 we have that M1 (h) = 0. Applying a suitable linear change of variables in the parameter space, we get  [1] 3 [4] 6  M1 (h, a) =λ[1] [2] 3 [3] 4 [2] 4 [3] 5 0 h + λ0 h + λ0 h + a λ 1 h + λ1 h + λ1 h + λ1 h , 2 that, with the adequate rescaling and taking b1 = 1, the function (7) is   [2]  [1]   [3] [2]  2 N[1] (h, μ) =h2 μ[1] 0 + μ0 + μ1 h + μ0 + μ1 h + μ1 h + μ1 h . [3] 3 [4] 4 Then the proof of the first part of the statement follows. We notice that the coefficients of N[1] are linearly independent, so the four limit cycles can be obtained easily. We finally remark that when a = 0, the linear independence of the coefficients of M1 provides the second part of the statement.  4. Higher-order developments This section is devoted to describing a partial higher-order analysis. It is known that the com- putation of the higher-order Melnikov function is usually a difficult task, sometimes impossible to be obtained explicitly. But, as we have previously explained, the study of the first Melnikov func- tion for center families gets some information about the complete higher-order analysis. More concretely, Theorem 4.1 generalizes Theorem 1.2 when the Taylor series of the first Melnikov function with respect to the parameters of the center family can be computed up to kth-order. We show an application to quadratic vector fields in Proposition 4.2 and two to cubic vector fields in Propositions 4.3 and 4.4. 303 J. Itikawa, R. Oliveira and J. Torregrosa Journal of Differential Equations 309 (2022) 291–310 Theorem 4.1. Let M1 (h, a, λ) be the function (5) associated to the perturbed differential equa- tions (3) with H = H (x, y, a), V = V (x, y, a), and P (x, y, λ), Q(x, y, λ), being a ∈ R and λ ∈ Rm . Let  k M[k] (h, a, λ) = Mα (h, λ)a α (14) |α|=0 be the Taylor series of M1 with respect to the multi-parameter a up to kth-order, where a α = a1α1 a2α2 · · · aα , α = (α1 , . . . , α ), and |α| = α1 + · · · + α . We assume that there exist functions [j ] [j ] fα and gα,β such that, for each α,  mα  m  α,β [j ] [j ] Mα (h, λ) = λ[j ] [j ] α fα (h) + λβ gα,β (h), j =1 |β|<|α| j =1 [j ] [j ] with β = (β1 , . . . , β ), |β| = β1 + · · · + β , and mα , mα,β ∈ N. Then taking λα = δ k−|α| μα , a = δb, for each positive simple zero h∗ of  k  mα N[k] (h, μ) = μ[j ] α [j ] α b fα (h), (15) |α|=0 j =1 and for a, ε small enough, there exists a limit cycle ε,h∗ of the perturbed system (3) such that ε,h∗ tends to the closed curve {H = h∗ } when (h, a) goes to (h∗ , 0). Proof. The displacement function (h, a, λ, ε) defined in (4) writes as (h, a, λ, ε) = ε(M[k] (h, a, λ) + O(ak+1 )) + O(ε 2 ), (16) with M[k] as in (14). Then, using the blow-up rescaling proposed in the statement together with ε = δ k+1 we have that (16), up to straightforward computations, writes as δ 2k+1 N[k] (h, μ) + O(δ 2k+2 ), with the function N[k] defined in (15). We observe that, from a = δb we have a α = δ |α| bα . The proof finishes as the proof of Theorem 1.2, applying the Implicit Function Theorem.  We remark that in the following applications we have taken concrete values of b, such as b = 1 and b = (1, 2). But examples where the number of zeros of (15) depends on the chosen b may exist. In fact we could use a weighted blow-up so the role of taken other b can help to study the new averaged function N[1] . The simplest cases where the above result applies are the ones for which H does not depend on a, only the inverse integrating factor V . Hence, the function (5) writes as   Q(x, y, λ) dx − P (x, y, λ) dy M1 (h, a) = = ω(x, y, λ, a). (17) V (x, y, a) H (x,y)=h H (x,y)=h 304 J. Itikawa, R. Oliveira and J. Torregrosa Journal of Differential Equations 309 (2022) 291–310 We can use Theorem 4.1 just integrating the Taylor series of the one-form ω(x, y, λ, a) over the ovals of H = h. The following applications are centers with an explicit parametrization of the level curves such that a curve of equilibrium points has been added to increase (artificially) the degree of the family. Obviously, the curve never passes through the center point. As is usual in these studies, when we have fixed the degree, the number of limit cycles under polynomial perturbations increases. The next example is the quadratic perturbation of the linear center with a straight line of equilibria. The first- and second-order analysis of this problem for degree n perturbations can be found in [12] and [13], respectively. The key point of both works is the explicit computation of the integrals. This fact is also observed in [17] where, as application, the quadratic perturbation is again analyzed, obtaining two limit cycles up to a second-order analysis. In Proposition 4.2 we show how with our approach also two limit cycles can be obtained without the difficulties of the second-order study. Because all our integrals are computed over circles. In particular we can see the existence of a second limit cycle for (simpler) third- and higher-order studies. We will see also that the optimal result is obtained analyzing third-order developments. Proposition 4.2. Consider the perturbed system (x  , y  ) = (ax + 1)(−y, x) + ε(P (x, y, λ), Q(x, y, λ)), (18) where P and Q are quadratic polynomials and ε = 0 is sufficiently small. The first new Melnikov functions (15) associated to system (18) can be written as N[0] (h, μ) = μ0 h2 , N[1] (h, μ) = μ0 h2 + μ1 h4 , (19) N[2] (h, μ) = μ0 h2 + (μ1 + μ2 )h4 , N[k] (h, μ) = μ0 h2 + (μ1 + μ2 )h4 + μ3 h6 , for k = 3, 4, 5, 6. Therefore, for ε and a small enough, there exist λ such that system (18) has no limit cycles up to order 0, one limit cycle up to first- and second-order and two limit cycles for kth-order studies with k = 3, 4, 5, 6. All limit cycles bifurcate from the period annulus up to a first-order analysis in ε.  2  2 Proof. We write P (x, y, λ) = pij x i y j , Q(x, y, λ) = qij x i y j . As usual, it is not re- i+j =0 i+j =0 strictive to assume P (0, 0, λ) = Q(0, 0, λ) = 0. Hence, λ ∈ R10 is the parameter space defined by the coefficients of the perturbative polynomials. System (18) admits H (x, y) = (x 2 + y 2 )/2 as a first integral and V (x, y, a, b) = ax + 1 as an inverse integrating factor. Clearly, system (18) when ε = 0 has a non-degenerate center at the origin. In fact, this system is topologically equivalent to the linear center adding the straight line ax + 1 = 0 of equilibria. Applying the usual polar change of coordinates (x, y) = (h cos θ, h sin θ ), by direct computa- tion we have that (14) can be written as 305 J. Itikawa, R. Oliveira and J. Torregrosa Journal of Differential Equations 309 (2022) 291–310  3 Si,j (λ)hi+j Ji,j (h), i+j =2 where Si,j (λ) are homogeneous polynomials of degree 1 in λ and 2π cosi θ sinj θ Ji,j (h) = dθ. (20) ah cos θ + 1 0 Instead of computing explicitly Ji,j (h) we get the Taylor series of it with respect to a and then we do the integral. So, the expression of M1 (h, a) defined in (17) is given as follows 1 M1 (h, a) = λ0 h2 + λ1 h4 a + λ2 h4 a 2 + λ3 h6 a 3 − h6 (λ0 − 8λ2 )a 4 8 (21) 5 5 − h8 (λ1 − 4λ3 )a 5 − h8 (λ0 − 6λ2 )a 6 + O(a 7 ). 16 32 To simplify reading, we have removed the superscript [1] in the above expression. Therefore, Theorem 4.1 provides the functions (19) and the conclusion of the number of limit cycles for each k. We remark that the obtained functions for k = 3, 4, 5, 6 are equal because the parameters of the coefficients of a 4 , a 5 , and a 6 in (21) are linear combinations of the independent ones which have appeared up to k = 3. Hence, all the new powers in h appear at the higher-order terms in δ 2k+2 and do not play any role when k increases.  The next example was analyzed firstly in [3] for degree n perturbations. Here we only recover partially some of the results because our perturbation is in the cubic family. We show how The- orem 4.1 can be applied to higher-order developments to simple families with two parameters. Proposition 4.3. Consider the perturbed system (x  , y  ) = (a1 x + 1)(a2 y + 1)(−y, x) + ε(P (x, y, λ), Q(x, y, λ)), (22) where P and Q are polynomials of degree 3 in (x, y) and ε = 0 sufficiently small. There exist linear changes of coordinates in the parameters λ such that the first functions (15) in Theorem 4.1 write as N[0] (h, μ) = h2 (η01 h2 + η00 ), N[1] (h, μ) = h2 (η11 h2 + η10 ), N[2] (h, μ) = h2 (η22 h4 + η21 h2 + η20 ), N[3] (h, μ) = h2 (η32 h4 + η31 h2 + η30 ), N[k] (h, μ) = h2 (η43 h6 + η42 h4 + η41 h2 + η40 ), for k = 4, 5, 6. Consequently, for ε and a = (a1 , a2 ) small enough, there exist values of the parameters λ such that system (22) has one limit cycle up to orders 0 and 1, two up to orders 2 and 3 and three up to 306 J. Itikawa, R. Oliveira and J. Torregrosa Journal of Differential Equations 309 (2022) 291–310 orders 4, 5, and 6. All limit cycles bifurcate from the period annulus up to a first-order analysis in ε.  3  3 Proof. As the previous proofs we take P (x, y, λ) = pij x i y j , Q(x, y, λ) = qij x i y j i+j =0 i+j =0 and P (0, 0, λ) = Q(0, 0, λ) = 0. Hence λ ∈ R18 belongs to the parameter space defined by the coefficients of the perturbative polynomials. System (22) has H (x, y) = (x 2 + y 2 )/2 and V (x, y, a1 , a2 ) = (a1 x + 1)(a2 y + 1) as first integral and inverse integrating factor, respectively. So, when ε = 0, we have a non-degenerate linear center at the origin and two straight lines of equilibria: a1x + 1 = 0 and a2 y + 1 = 0. Applying the usual polar change of coordinates (x, y) = (h cos θ, h sin θ ), by direct computa- tion we have that (14) is written as  3 Si,j (λ)hi+j Ji,j (h), i+j =2 where Si,j (λ) are homogeneous polynomials of degree 1 in λ and 2π cosi θ sinj θ Ji,j (h) = dθ. (23) (a1 h cos θ + 1)(a2 h sin θ + 1) 0 Integrating after computing the Taylor series of the integrand with respect to a and applying a linear change of coordinates, we obtain the expression of M1(h, a) defined in (17) as follows M[5] (h, a, λ) = h2 (λ[2] [1] [1] 4 [1] 4 00 h + λ00 ) + λ10 h a1 + λ01 h a2 2 1 + h4 (2(λ[2] [2] 2 [1] [1] 2 00 + 2λ20 )h + 3λ00 + 4λ20 )a1 4 1 4 + h4 (λ[2] [1] [2] [2] 2 [1] [1] 2 11 h + λ11 )a1 a2 + h (2(λ00 − 2λ20 )h + λ00 − 4λ20 )a2 2 4 h6 h6 h6 + (2λ[1] [1] 3 30 + λ10 )a1 + (2λ[1] [1] 2 21 + λ01 )a1 a2 − (2λ[1] [1] 30 − λ10 )a1 a2 2 2 2 2 h6 h6 − (2λ[1] 21 − λ [1] 3 )a 01 2 + ((5λ[2] [2] [1] 2 [1] [1] 4 00 + 20λ20 + 16λ40 )h + 10λ00 + 16λ20 )a1 2 16 h6 h6 + ((3λ[2] [1] 2 [1] 3 11 + 8λ31 )h + 4λ11 )a1 a2 + ((3λ[2] [2] [1] 2 [1] 2 2 00 − 4λ20 − 16λ40 )h + 2λ00 )a1 a2 8 16 h6 h6 + ((5λ[2] [1] 2 [1] 11 − 8λ31 )h + 4λ11 )a1 a2 + 3 ((5λ[2] [2] [1] 2 [1] [1] 4 00 − 12λ20 + 16λ40 )h + 2λ00 − 16λ20 )a2 8 16 5h8 h8 h8 + (4λ[1] 30 + λ [1] 5 10 )a 1 + (12λ [1] 21 + 5λ [1] 4 01 )a 1 a 2 + (−4λ[1] [1] 3 2 30 + 3λ10 )a1 a2 16 16 16 h8 h8 5h8 + (4λ[1] 21 + 3λ [1] 2 3 )a 01 1 2a − (12λ [1] 30 − 5λ [1] 10 )a 1 a 4 2 − (4λ[1] [1] 5 21 − λ01 )a2 . 16 16 16 307 J. Itikawa, R. Oliveira and J. Torregrosa Journal of Differential Equations 309 (2022) 291–310 The last step is to establish an adequate blow-up to use Theorem 4.1. We have chosen (a1 , a2 ) = δ(1, 2) and, straightforward computations, the new Melnikov functions (15) are N[0] (h, μ) = h2 (μ[2] [1] 00 h + μ00 ), 2 N[1] (h, μ) = h2 ((μ[2] [1] [1] 2 [1] 00 + 2μ01 + μ10 )h + μ00 ), N[2] (h, μ) = h2 ((2μ[2] [2] 4 [2] [1] [1] [1] [1] 2 [1] 11 − 3μ20 )h + (μ00 + 2μ01 + μ10 + 2μ11 − 3μ20 )h + μ00 ), N[3] (h, μ) = h2 ((2μ[2] [2] [1] [1] 4 11 − 3μ20 − 6μ21 − 3μ30 )h + (μ[2] [1] [1] [1] [1] 2 [1] 00 + 2μ01 + μ10 + 2μ11 − 3μ20 )h + μ00 ), N[k] (h, μ) = h2 ((−6μ[1] [1] 6 [2] [2] [1] [1] 4 31 + 13μ40 )h + (2μ11 − 3μ20 − 6μ21 − 3μ30 )h + (μ[2] [1] [1] [1] [1] 2 [1] 00 + 2μ01 + μ10 + 2μ11 − 3μ20 )h + μ00 ), for k = 4, 5. We have computed also the sixth-order function and it coincides with the fourth-order above. We have not written here the complete development because of its size and no better results can be obtained. Renaming the coefficients of the functions N[k] the result follows.  It can be checked also that if in the last proof we take (a1 , a2 ) = δ(1, 1) not all the independent [j ] parameters μkl appear, but although the functions N[k] in the proof are different, the ones given in the statement are the same. Hence, the same number of zeros are obtained. As in the previ- ous quadratic family, the increasing number of zeros stop when there are no new independent parameters in the higher-order developments. Our last application of Theorem 4.1 is to the same family of the previous result but fixing a2 = 1. In particular, all the integrals over the closed paths can be explicitly computed and the result is quite better than the obtained in Proposition 4.3. As we will see in the proof, the number of zeros increases up to stabilizes when no more new independent parameters appear in the higher-order developments. Proposition 4.4. Consider the perturbed system (x  , y  ) = (ax + 1)(y + 1)(−y, x) + ε(P (x, y, λ), Q(x, y, λ)), (24) where P and Q are polynomials of degree 3 in (x, y) and ε = 0 sufficiently small. √ There exists a linear change of coordinates in the parameters λ such that, taking r = 1 − h2 , the first functions (15) in Theorem 4.1 can be written as follows: N[0] (r, η) = (1 − r)(η0[4] r 3 + η0[3] r 2 + η0[2] r + η0[1] )/r, N[1] (r, η) = (1 − r)(η1[5] r 4 + η1[4] r 3 + η1[3] r 2 + η1[2] r + η1[1] )/r, N[2] (r, η) = (1 − r)(η2[5] r 4 + η2[4] r 3 + η2[3] r 2 + η2[2] r + η2[1] )/r, N[k] (r, η) = (1 − r)(η3[6] (1 + r)r 5 + η3[5] r 4 + η3[4] r 3 + η3[3] r 2 + η3[2] r + η3[1] )/r, for k = 3, 4, 5. Consequently, for ε and a small enough, there exist λ such that system (24) has three limit cycles up to order 0, four up to orders 1 and 2 and five up to orders 3, 4, and 5. All limit cycles bifurcate from the period annulus up to a first-order analysis in ε. 308 J. Itikawa, R. Oliveira and J. Torregrosa Journal of Differential Equations 309 (2022) 291–310 Proof. The proof follows exactly as the proof of Proposition 4.3 taking a1 = a and a2 = 1 up to getting the expression of the function M1 (h, a) defined in (17). In the Taylor development of integrals of type (23) with a2 = 1 with respect to a1 = a appear integrals of type (20) but with√h sin θ + 1 in the denominator. After straightforward computations and taking k = 5 and h = 1 − r 2 we get, for 0 < r < 1, (1 − r)  [4] 3 [1]   [1]  M[5] (r, a, λ) = λ0 r + λ[3] [2] 0 r + λ0 r + λ0 2 + (1 − r)2 λ[3] [2] 1 r + λ1 r + λ1 a 2 r (1 − r)2  [4] 3 + 4λ0 r + (4λ[1] [3] [4] 2 [1] 2 + 3λ0 + 2λ0 )r + 2(4λ2 + 2λ0 + λ0 )r + 2λ0 [2] [3] [1] 4 [3]  2 (1 − r)3  + 4λ[1] [2] 2 + 2λ0 + λ0 a + 2(2λ[1] [3] 3 [1] 3 + λ1 )r + 2(6λ3 + 2λ1 − λ1 )r [2] [3] 2 4 [3]  3 + 3(λ[1] [1] [2] [3] [1] 1 + 4λ3 + λ1 − λ1 )r + λ1 + 4λ3 + λ1 − λ1 a [1] [2] (1 − r)3  4 [4] + 8r λ0 + (8λ[1] [3] [4] 3 [1] [2] [3] [4] 2 2 + 5λ0 + 9λ0 )r + (24λ2 + 8λ0 + 7λ0 + 3λ0 )r 8 [3]  4 + (3λ[1] [1] [2] [3] [1] [1] [2] 0 + 24λ2 + 9λ0 + 6λ0 )r + λ0 + 8λ2 + 3λ0 + 2λ0 a (1 − r)4  + 5(4λ[1] [3] 4 [1] [2] [3] 3 3 + λ1 )r + 4(20λ3 + 4λ1 − 3λ1 )r 16 + 2(5λ[1] [1] [2] [3] 2 [1] [1] [2] [3] 1 + 60λ3 + 12λ1 − 14λ1 )r + 4(2λ1 + 20λ3 + 4λ1 − 5λ1 )r [3]  5 + 2λ[1] [1] [2] 1 + 20λ3 + 4λ1 − 5λ1 a . The next step is to choose the blow-up indicated in Theorem 4.1 with b = 1. Hence, the new ] (r, μ) where Melnikov functions (15) are N[k] (r, μ) = (1 − r)r −1 N [0] (r, μ) = μ[4] r 3 + μ[3] r 2 + μ[2] r + μ[1] , N 0 0 0 0 [1] (r, μ) = −μ[3] r 4 − (μ[2] − μ[3] − μ[4] )r 3 − (μ[1] − μ[2] − μ[3] )r 2 + (μ[1] + μ[2] )r +μ[1] , N 1 1 1 0 1 1 0 1 0 0 [2] (r, μ) = −(μ[1] + μ[3] )r 4 − (μ[1] + μ[2] − μ[3] − μ[4] )r 3 − (μ[1] − μ[1] − μ[2] − μ[3] )r 2 N 2 1 2 1 1 0 1 2 1 0 + (μ[1] [1] [2] [1] 1 + μ2 + μ0 )r + μ0 , [3] (r, μ) = μ[1] r 6 + μ[1] r 5 − (μ[1] + 2μ[1] + μ[3] )r 4 − (μ[1] + 2μ[1] + μ[2] − μ[3] − μ[4] )r 3 N 3 3 2 3 1 2 3 1 1 0 − (μ[1] [1] [1] [2] [3] 2 [1] [1] [1] [2] [1] 1 − μ2 − μ3 − μ1 − μ0 )r + (μ1 + μ2 + μ3 + μ0 )r + μ0 , and N[3] = N[4] = N [5] . A last linear change of coordinates provides the functions N[k] indicated in the statement and the conclusion on the number of limit cycles follows immediately for k = 0, 1, 2. For k = 3, 4, 5, Descartes’ rule ensures that we have at most five positive zeros of the corresponding N[k] function and it is easy to check that they exist as simple ones.  Acknowledgments This work has been realized thanks to the Spanish Ministerio de Ciencia, Innovación y Uni- versidades - Agencia Estatal de Investigación grants, PID2019-104658GB-I00 and CEX2020- 309 J. Itikawa, R. Oliveira and J. Torregrosa Journal of Differential Equations 309 (2022) 291–310 001084-M; the Generalitat de Catalunya - AGAUR grant 2017SGR1617; the European Com- munity Marie Skłodowska-Curie grant Dynamics-H2020-MSCA-RISE-2017-777911; and the Brazilian grants FAPESP 2019/21181-0 and CNPq 304766/2019-4. We thank professors Carlos H. Grossi and Daniel Peralta-Salas for their helpful discussions during the realization of this paper. References [1] V.I. Arnol’d, Geometrical Methods in the Theory of Ordinary Differential Equations, Grundlehren der Math- ematischen Wissenschaften (Fundamental Principles of Mathematical Science), vol. 250, Springer-Verlag, New York-Berlin, 1983, translated from the Russian by Joseph Szücs, translation edited by Mark Levi. [2] A. Buică, J. Llibre, Averaging methods for finding periodic orbits via Brouwer degree, Bull. Sci. Math. 128 (1) (2004) 7–22. [3] A. Buică, J. Llibre, Limit cycles of a perturbed cubic polynomial differential center, Chaos Solitons Fractals 32 (3) (2007) 1059–1069. [4] F. Dumortier, R. Roussarie, Abelian integrals and limit cycles, J. Differ. Equ. 227 (1) (2006) 116–165. [5] H. Flanders, Differentiation under the integral sign, Am. Math. Mon. 80 (1973) 615–627; correction, Am. Math. Mon. 81 (1974) 145. [6] J. Giné, L.F.S. Gouveia, J. Torregrosa, Lower bounds for the local cyclicity for families of centers, J. Differ. Equ. 275 (2021) 309–331. [7] M. Han, P. Yu, Normal Forms, Melnikov Functions and Bifurcations of Limit Cycles, Applied Mathematical Sci- ences, vol. 181, Springer, London, 2012. [8] D. Hilbert, Mathematical problems, Bull. Am. Math. Soc. 8 (10) (1902) 437–479. [9] Y. Ilyashenko, Centennial history of Hilbert’s 16th problem, Bull. Am. Math. Soc. (N.S.) 39 (3) (2002) 301–354. [10] J. Li, Hilbert’s 16th problem and bifurcations of planar polynomial vector fields, Int. J. Bifurc. Chaos Appl. Sci. Eng. 13 (1) (2003) 47–106. [11] J. Llibre, D.D. Novaes, M.A. Teixeira, On the birth of limit cycles for non-smooth dynamical systems, Bull. Sci. Math. 139 (3) (2015) 229–244. [12] J. Llibre, J.S. Pérez del Río, J.A. Rodríguez, Averaging analysis of a perturbated quadratic center, Nonlinear Anal., Theory Methods Appl. 46 (1) (2001) 45–51. [13] J. Llibre, J. Yu, On the upper bound of the number of limit cycles obtained by the second order averaging method, Dyn. Contin. Discrete Impuls. Syst., Ser. B, Appl. Algorithms 14 (6) (2007) 841–873. [14] D. Marín, J. Villadelprat, Asymptotic expansion of the dulac map and time for unfoldings of hyperbolic saddles: general setting, J. Differ. Equ. 275 (2021) 684–732. [15] H. Poincaré, Mémoire sur les courbes définies par une équation différentielle, J. Math. Pures Appl. (1) 7 (1881) 375–422, (2) 3 (1882) 251–296, (3) 4 (1885) 167–244, (4) 2 (1886) 151–217. [16] L.S. Pontrjagin, Über autoschwingungssysteme, die den hamiltonshen nahe liegen, Phys. Z. Sowjetunion 6 (1–2) (1934) 25–28. [17] R. Prohens, J. Torregrosa, Periodic orbits from second order perturbation via rational trigonometric integrals, Phys- ica D 280/281 (2014) 59–72. 310