Qualitative Theory of Dynamical Systems (2022) 21:20 https://doi.org/10.1007/s12346-021-00546-x Simultaneous Bifurcation of Limit Cycles and Critical Periods Regilene D. S. Oliveira1 · Iván Sánchez-Sánchez2 · Joan Torregrosa2,3 Received: 16 April 2021 / Accepted: 12 November 2021 / Published online: 19 December 2021 © The Author(s) 2021 Abstract The present work introduces the problem of simultaneous bifurcation of limit cycles and critical periods for a system of polynomial differential equations in the plane. The simultaneity concept is defined, as well as the idea of bi-weakness in the return map and the period function. Together with the classical methods, we present an approach which uses the Lie bracket to address the simultaneity in some cases. This approach is used to find the bi-weakness of cubic and quartic Liénard systems, the general quadratic family, and the linear plus cubic homogeneous family. We finish with an illustrative example by solving the problem of simultaneous bifurcation of limit cycles and critical periods for the cubic Liénard family. Keywords Limit cycles · Critical periods · Simultaneous bifurcation Mathematics Subject Classification 34C07 · 34C23 · 34C25 1 Introduction A classical problem in the study of qualitative theory of planar differential equations in the plane is the 16th Hilbert Problem, related to the bifurcation of limit cycles (isolated B Joan Torregrosa torre@mat.uab.cat Regilene D. S. Oliveira regilene@icmc.usp.br Iván Sánchez-Sánchez isanchez@mat.uab.cat 1 Departamento de Matemática, Instituto de Ciências Matemáticas e Computação, Universidade de São Paulo, Avenida Trabalhador São Carlense, São Carlos, SP, 400 13566-590, Brazil 2 Departament de Matemàtiques, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Catalonia, Spain 3 Centre de Recerca Matemàtica, Campus de Bellaterra, 08193 Bellaterra, Barcelona, Spain 20 Page 2 of 35 R. D. S. Oliveira et al. periodic orbits) in a polynomial class of fixed degree. A large number of works in this line have been published so far for several polynomial families of differential equations. A different problem that has aroused interest during the last decades is the study of the isochronicity of a system, as well as its bifurcation of critical periods. These problems consist on analyzing the flatness and the oscillations of the period function of the system, respectively. As a matter of fact, the bifurcation of limit cycles and critical periods are analogous in terms of the techniques that can be used to be approached. For this reason, in this paper we suggest the study of the bifurcation of limit cycles and critical periods simultaneously, a problem that to the best of our knowledge has not been formulated yet. Let us consider a real polynomial system of differential equations in the plane whose origin is a nondegenerate monodromic equilibrium point, so the matrix associated to the differential system evaluated at the origin has zero trace and positive determinant. It is a well-known fact that, by a suitable change of coordinates and time rescaling, it can be written in the form  ẋ = αx − y + X (x, y) =: P(x, y), (1) ẏ = x + α y + Y (x, y) =: Q(x, y), where X and Y are polynomials of degree n ≥ 2 which start at least with quadratic monomials and the dot indicates the derivative with respect to the time. We can consider system (1) in complex coordinates (z, w) = (z, z) = (x + i y, x − i y), which will be represented by only one equation as ż = (α + i)z + Z (z, w) = Z(z, w), (2) where Z is a polynomial starting with monomials of at least second degree. For all the results developed in this work, we will consider a fixed transversal section x defined as the positive x-axis, this is the semi-axis {(x, 0), x > 0}. This restriction is necessary for the study of the period function because when the origin of (1) is not a center we can construct an isochronous section, as can be seen in [3]. Then, for example, the bifurcation of critical periods can depend on the transversal section. Consider system (1), with α = 0, in polar coordinates,  ṙ = R(r , ϕ), (3) ϕ̇ = 1 + (r , ϕ). It is straightforward that from this system we can write dr = a2 (ϕ)r 2 + a3 (ϕ)r 3 + O4 (r ), (4) dϕ where we denote by O j (r ) the monomials in r of at least degree j. From here we deduce the solution r (ϕ, ρ) = ρ + v2 (ϕ)ρ 2 + v3 (ϕ)ρ 3 + O4 (ρ), (5) Simultaneous Bifurcation of Limit Cycles and Critical Periods Page 3 of 35 20 being ρ = r (0) the initial condition. Evaluating it at 2π we obtain the so-called Poincaré return map: r (2π, ρ) = ρ + V2k+1 ρ 2k+1 + O2k+2 (ρ). (6) As we will see in Sect. 2, the first nonvanishing term of the displacement map d(ρ) = r (2π, ρ) − ρ has odd degree in ρ. See more details in [28]. Solution (5) can be considered on the second equation of (3) to obtain the expression of the time t in terms of the angle ϕ. Let us invert such equation, dϕ   dt = = 1 + τ1 (ϕ)r + τ2 (ϕ)r 2 + O3 (r ) dϕ, 1 + (r , ϕ) and by integrating we obtain the expression of the time, the so-called period function, as a Taylor series in ρ, this is T (ρ) = 2π + Tl ρ l + Ol+1 (ρ). (7) When α = 0, the origin of (1) is clearly a hyperbolic focus. When α = 0 we are in the center-focus case. Using the expressions of the return and period maps given in (6) and (7), the origin of system (1) satisfies one of the following properties, with respect to the transversal section x : L The origin is an isochronous center, so the system is linearizable and we can say that k = ∞ and l = ∞. C The origin is a center with weakness of order l on the period, so k = ∞ and 1 ≤ l < ∞. W The origin is an isochronous weak focus of order k, so 1 ≤ k < ∞ and l = ∞. B The remaining case is when both the center and isochronicity properties are not kept at the same time, so 1 ≤ k, l < ∞ and we can say that the origin is a bi-weak monodromic equilibrium point of type [k, l]. It is usual to restrict the study of the period function to the class of centers, i.e. systems that remain in type C in the aforementioned canonical form changes of vari- ables classification. The study of the monotonicity or the number of oscillations of such function are difficult problems, see for example [6,11,26,36] and the references therein. There are not so many global studies of the period function for general classes of centers. Gavrilov ([18]) proved the existence of at most one critical period for the Hamiltonian potentials x  + x + ax 2 + bx 3 = 0, a problem started by Chow and Sanders ([8]) in 1986. In 2006, Mañosas and Villadelprat ([25]) proved that the deriva- tive of the period function for Hamiltonian potentials x  + x +ax 3 +bx 5 = 0 has only one zero. Some years later, Grau and Villadelprat ([22]) proved that only two critical periods appear in some cubic homogeneous nonlinearity classes. In those cases, we say that the systems have one and two critical periods, respectively. For centers in the quadratic class, the most relevant study was done by Chicone and Jacobs in 1989 ([7]). Among others, they studied the local problem for the quadratic family, proving that only two critical periods bifurcate from the center equilibrium point. The answer for 20 Page 4 of 35 R. D. S. Oliveira et al. the global problem remains open. The greatest difficulty to deal with is the fact that the outer boundary of the period annulus changes together with the parameters inside a fixed family, see for example [15,33]. Hence, the usual perturbation techniques are not useful and new tools need to be developed ([27]). As we have described, the maximal number of zeros of the derivative of the period function under perturbation (in some fixed class) is known as the criticality of the center. Recently, this problem has been studied for low degree polynomial vector fields in the class of reversible centers, see [31,32]. Similarly to the above problem, we can restrict our analysis to the class of vector fields that remain in type W. This is the case associated to the problem of studying isochronous foci, a problem that was addressed for example by Giné in [19–21]. In this special class, the cyclicity problem is also an interesting problem to be approached, which up to our knowledge is not completely solved even for low degree vector fields. A special family of systems in this class are the so-called rigid (or uniformly isochronous) systems. They satisfy that ϕ̇ = 1 in the usual polar coordinates. Inside this class, quadratics have no limit cycles and there are cubics with at least two ([16]), but there is no answer for the global question about the total number of limit cycles in rigid cubic systems. When α = 0, the origin of system (1) can be a either a center or a weak focus of a certain order. This classical notion of order will be recalled in the next section, where we will generalize it to the bi-weakness property previously defined. This gives birth to the idea of duality in weakness of a nonisochronous focus. In this work we will restrict our analysis to the case in which the transversal section is the horizontal axis. As we have previously commented and as we will see in Sect. 2, the first nonzero term of the displacement map, defined from (6), always has an odd exponent. The coefficient V2k+1 is known as the kth Lyapunov constant. We also know that the first nonzero term of the derivative of the period function (7) is the so-called lth period constant Tl . For centers it corresponds to an even subindex l ([1]) with k = ∞. However, this is not the case for systems which do not have a center at the origin. For instance, let us consider the quadratic system  ẋ = −y + x 2 − 10 9 x y, (8) ẏ = x + x + 4x y − 25 2 9 y . 2 For this system, the origin is bi-weak of type [1, 3], because we can easily find that the first nonzero Lyapunov constant is V3 = π and the first nonzero coefficient of the period function is T3 = π, so the origin is not a center because V3 = 0 and we have T3 = 0. Therefore, if the center property is not kept then the property which states that the first nonzero period constant has even subindex does not hold. The main purpose of this work is to present the problems of simultaneous cyclicity and criticality and bi-weak equilibrium points of type [k, l], as well as some useful tools to deal with them. The first result we present is related to bi-weak [k, l] foci in Liénard, quadratic, and linear plus cubic homogeneous systems, and our aim is then to find the highest [k, l] before the either the center or the isochronicity (or both) occur. Simultaneous Bifurcation of Limit Cycles and Critical Periods Page 5 of 35 20 This problem is introduced in Sect. 2 using the Lie bracket, with a method which only finds bi-weak points of type [m, 2m]. Theorem 1 (i) There exist cubic Liénard systems of the form  ẋ = −y + a2 x 2 + a3 x 3 , (9) ẏ = x + b2 x 2 + b3 x 3 , being a2 , a3 , b2 , b3 ∈ R, which have a bi-weak [2, 4] focus at the origin. (ii) There exist quartic Liénard systems of the form  ẋ = −y + a2 x 2 + a3 x 3 + a4 x 4 , (10) ẏ = x + b2 x 2 + b3 x 3 + b4 x 4 , being a2 , a3 , a4 , b2 , b3 , b4 ∈ R, which have a bi-weak [3, 6] focus at the origin. (iii) There exist quadratic systems of the form  ẋ = −y + a20 x 2 + a11 x y + a02 y 2 , (11) ẏ = x + b20 x 2 + b11 x y + b02 y 2 , being a20 , a11 , a02 , b20 , b11 , b02 ∈ R, which have a bi-weak [2, 4] focus at the origin. (iv) There exist linear plus cubic homogeneous systems of the form  ẋ = −y + a30 x 3 + a21 x 2 y + a12 x y 2 + a03 y 3 , (12) ẏ = x + b30 x 3 + b21 x 2 y + b12 x y 2 + b03 y 3 , being a30 , a21 , a12 , a03 , b30 , b21 , b12 , b03 ∈ R, which have a bi-weak [3, 6] focus at the origin. The second result deals with the simultaneous cyclicity and criticality for cubic Liénard systems. Theorem 2 For the cubic Liénard family (9), adding the trace parameter α as in (1), we have the following properties with respect to the transversal section x : (i) When the origin is a center, either it is isochronous or it has at most 1 critical period. (ii) There at most two limit cycles of small amplitude bifurcating from the origin. Simultaneously, there are at least three critical periods also bifurcating from the origin. (iii) The origin is never an isochronous focus. This paper is structured as follows. In Sect. 2 we introduce the classical methods to find Lyapunov and period constants as well as a method which uses the Lie bracket. We also discuss the pros and cons of each of the approaches. The presented Lie bracket 20 Page 6 of 35 R. D. S. Oliveira et al. method is used in Sect. 3 to study the bi-weakness of Liénard, quadratic, and linear plus cubic homogeneous systems in order to prove Theorem 1. Finally, in Sect. 4 we show some properties of several Liénard systems and give the isochronicity conditions for some cases, and finish by proving Theorem 2. We remark however that the aim of this work is not to provide any isochronicity characterization, not even for an a priori simple class of systems as the Liénard family, a characterization of which has been obtained very recently in [4]. 2 Lyapunov and Period Constants Computation. The Lie Bracket Method In this section we retrieve some classical concepts related with the Taylor developments of the return and period maps. The coefficients of these series, as we will see in Sect. 2.2, allow us to define the Lyapunov and period constants. These values will be computed simultaneously in Sect. 2.3 by using the Lie bracket operator, which is recalled in complex coordinates the following Sect. 2.1 together with the characterization of isochronous centers by using commutators. See more details in [1,30]. 2.1 The Lie Bracket In this subsection we present the Lie bracket, a powerful tool to study the isochronicity of a system and to find its Lyapunov and period constants under certain conditions that we will see. Some parts of this introduction have been directly extracted from a previous work [31]. We define the Lie bracket of two planar vector fields F1 = (F11 , F12 ) and F2 = (F21 , F22 ) in variables (x, y) ∈ K2 , where K = R or C, as a new vector field which has the form  1 ∂ F1 1 ∂ F11 2 ∂ F21 1 ∂ F21 2 [F1 , F2 ] = F + F − F − F , ∂x 2 ∂y 2 ∂x 1 ∂y 1  (13) ∂ F12 1 ∂ F12 2 ∂ F22 1 ∂ F22 2 F + F − F − F . ∂x 2 ∂y 2 ∂x 1 ∂y 1 Let us consider now the Lie bracket in the case of having two complex planar vector fields Z, U, corresponding to two real vector fields. Observe that, in such situation, second components are obtained by complex conjugation of first components, thus both vector fields Z and U and their Lie bracket can be described only from their first components, so in some cases we will simply write ∂Z ∂Z ∂U ∂U [Z, U] = U+ U− Z− Z. (14) ∂z ∂w ∂z ∂w This definition also appears in [13]. The first proof of the next geometrical equivalence was done by Sabatini in [30]. Simultaneous Bifurcation of Limit Cycles and Critical Periods Page 7 of 35 20 Theorem 3 ([1]) Equation (2) with α = 0 has an isochronous center at the origin if and only if there exists ż = U(z, w) = z + O(|z, w|2 ) such that [Z, U] = 0. It is a well-known fact that holomorphic systems are isochronous, see for example [1]. The paper [14] also deals with this problem, and gives ż = i f (z) as a linearizing system. As a first example of application of Theorem 3, we can prove this same result and see that actually U : ż = k f (z) for any k ∈ C is a linearizing system. Indeed, ∂ f (z) ∂ f (z) ∂ (k f (z)) [Z, U] = [ f (z), k f (z)] = k f (z) + k f (z) − f (z) ∂z ∂w ∂z ∂ (k f (z)) ∂ f (z) ∂ f (z) − f (z) = k f (z) + 0 − k f (z) − 0 = 0. ∂w ∂z ∂z 2.2 The Classical Method to Compute Lyapunov and Period Constants We start by presenting the classical method of finding Lyapunov and period constants. Let us write system (1), with α = 0, in polar coordinates by performing the usual change (x, y) = (r cos ϕ, r sin ϕ), and one obtains ⎧ n−1 ⎪ ⎪ ⎨ṙ = Ui (ϕ)r i+1 , i=1 (15) ⎪ ⎪ n−1 ⎩ϕ̇ = 1 + Wi (ϕ)r i , i=1 where Ui (ϕ) and Wi (ϕ) are homogeneous polynomials in sin ϕ and cos ϕ of degree i + 2. Eliminating time and doing the Taylor series expansion in r we obtain ∞ dr = R j (ϕ)r j . (16) dϕ j=2 The initial value problem for (16) with the initial condition (r , ϕ) = (ρ, 0) has a unique analytic solution which can be expanded as ∞ r (ϕ, ρ) = ρ + u j (ϕ)ρ j . (17) j=2 As r (0, ρ) = ρ, it immediately follows that u j (0) = 0 for every j. Let us study the stability near the origin, r = 0, by using r (2π, ρ) = ρ + Vm ρ m + Om+1 (ρ), (18) where Vm := u m (2π ) is the first coefficient which does not vanish. It is a well- known fact that the first nonidentically zero Vm has odd m ([5,28]). As a consequence, expression (18) can be rewritten as 20 Page 8 of 35 R. D. S. Oliveira et al. r (2π, ρ) = ρ + V2k+1 ρ 2k+1 + O2k+2 (ρ). We will denote the V2k+1 with odd subscript as L k := V2k+1 , and these values are the Lyapunov constants. These objects are the key tool to study the center and cyclicity problems of a system of the form (1). Observe that r (2π, ρ) indicates the radius after a whole loop starting in the initial value ρ, and then we define the Poincaré map ∞ (ρ) := ρ + Vj ρ j . j=3 Alternatively, this can be written as ∞ d(ρ) := (ρ) − ρ = Vj ρ j, (19) j=3 which is the so-called displacement map. We recall that we are taking a nondegenerate equilibrium point with zero trace, i.e. α = 0 in (1). Now let us illustrate how the period constants can be found, for which the reader is referred to [24,28,31]. First, we substitute the power series (17) into the second equation of (15), which yields a differential equation of the form ∞ dϕ =1+ Fi (ϕ)ρ i , dt i=1 for some trigonometric polynomials Fi (ϕ). Rewriting this equation as ∞  dϕ dt = ∞ = 1+ i (ϕ)ρ i dϕ 1+ Fi (ϕ)ρ i i=1 i=1 and integrating ϕ from 0 to 2π yields the period function   ∞  ∞   T (ρ) 2π 2π T (ρ)= dt= 1+ i (ϕ)ρ i dϕ = 2π + i (ϕ) dϕ ρi , 0 0 i=1 i=1 0 (20) where the series converges for 0 ≤ ϕ ≤ 2π and sufficiently small ρ ≥ 0. In the center case, (20) can be seen as the lowest period of the trajectory of (1) passing through (x, y) = (ρ, 0) for ρ = 0. The coefficients Ti of the period function are then given by the expression  2π 1 Ti = i (ϕ) dϕ, 2π 0 Simultaneous Bifurcation of Limit Cycles and Critical Periods Page 9 of 35 20 and the first nonzero Tl is known as the lth period constant of the system. If we assume now that system (15) depends on some parameters λ ∈ Rd , we can follow exactly the same procedure as before, and now we have that both the Lyapunov constants V2k+1 = u 2k+1 (2π, λ) and the period constants  2π 1 Tl (λ) = l (ϕ, λ) dϕ (21) 2π 0 are polynomials in the parameters λ (see [10,28]). Finally, we consider the case with nonzero trace (α = 0). The following result stud- ies how the return map and the period function, together with cyclicity and criticality, change when the considered system has nonzero trace. Lemma 4 Let us consider a system of the form (1) written as  ẋ = αx − y + a20 x 2 + a11 x y + a02 y 2 + X (x, y), (22) (x, y), ẏ = x + α y + b20 x 2 + b11 x y + b02 y 2 + Y where   are polynomials which start at least with cubic monomials, and such X and Y that when α = 0 the first non vanishing Lyapunov and period constants are V3 and T2 . (i) The displacement map (19) when α = 0 has the form d(ρ) = (e2π α −1)ρ + O2 (ρ). (23) Moreover, if α V3 < 0 and α is small enough, then a unique limit cycle bifurcates from the origin due to a Hopf bifurcation. (ii) The period function (20) when α = 0 has the form T (ρ) = 2π + (e2π α −1)T1 (α) ρ + O2 (ρ), (24) where −α 3 b20 − α 2 (a20 − b11 ) + α(2a11 − 2b02 − 7b20 ) − 6a02 − 3a20 + 3b11 T1 (α) = . α 4 + 10α 2 + 9 Moreover, if (e2π α −1) T1 (α) T2 < 0 and α is small enough, then a unique critical period bifurcates from the origin. Proof It is a well-known fact and a straightforward computation that the displacement map of system (22) has the form (23). Thus, as α and e2π α −1 have the same sign, it is immediate to see that if α V3 < 0 then d(ρ) has an extra change of sign in the coefficients of ρ, so a unique new positive zero can appear and this implies the unfolding of a limit cycle of small amplitude. See more details in [5,29]. For the case of the period, we can change system (22) to polar coordinates and integrate the angular equation of dϕ/dt as we did above. After doing the computations, 20 Page 10 of 35 R. D. S. Oliveira et al. we obtain that the period function in the case of nonzero trace becomes (24). Notice that in (22) we have only considered the coefficients of the quadratic part because by construction they are the only ones which can actually play a part in the coefficient of ρ in the period function, in the sense that higher degree coefficients would appear in higher degree terms of the period function. Now if (e2π α −1)T1 (α) T2 < 0, there is an extra change of sign in the coefficients of the period function, which implies the bifurcation of an extra critical period.  2.3 Lyapunov and Period Constants via the Lie Bracket In the previous section we presented the classical method to compute Lyapunov and period constants. Such method involves some integrals (equation (20)), which easily become too difficult to be explicitly obtained, so this technique is not useful in many cases. To deal with this inconvenience, a method to find period constants which is based on the Lie bracket tool was introduced in [24] and recently used in [31,32]. However, this technique was only valid when the origin is a center, which is not our case in this paper, since we aim to study cyclicity and criticality simultaneously. Here we present a new approach based on the Lie bracket which will allows to find both Lyapunov and period constants at the same time. This method will provide some valuable advantages, even though it also has its limitations. Let us consider system (2) with α = 0. By applying near the identity changes of variables, as the spirit of normal form transformations, such system can be simplified to N ż = i z + (α2 j+1 + i β2 j+1 )z(zw) j + O2N +3 (z, w), (25) j=1 where N ∈ N is arbitrary, O2N +3 (z, w) denotes a sum of monomials of degree at least 2N + 3 in z and w, and α2 j+1 , β2 j+1 ∈ R (see [1] for more details). Let us consider the N th order truncation of differential system (25) N ż = i z + (α2 j+1 + i β2 j+1 )z(zw) j . (26) j=1 Proposition 5 The truncated system (26) in polar coordinates takes the form ⎧ ⎪ ⎪ N ⎨ṙ = ⎪ α2 j+1r 2 j+1 , j=1 (27) ⎪ ⎪ N ⎪ ⎩ϕ̇ = 1 + β2 j+1r 2 j . j=1 Proof The proof is straightforward by applying the change of variables   (r , ϕ) → z = r ei ϕ , w = r e− i ϕ Simultaneous Bifurcation of Limit Cycles and Critical Periods Page 11 of 35 20 to differential equation (27). Indeed, ⎛ ⎞ ⎛ ⎞ N N ż = ṙ ei ϕ + i r ϕ̇ ei ϕ = ⎝ α2 j+1 r 2 j+1 ⎠ ei ϕ + i r ei ϕ ⎝1 + β2 j+1 r2 j ⎠ j=1 j=1 ⎛ ⎞ ⎛ ⎞ N N N =z⎝ α2 j+1 r 2 j ⎠ + i z ⎝1 + β2 j+1 r 2 j ⎠ = i z + (α2 j+1 + i β2 j+1 )z(zw) j , j=1 j=1 j=1 where we have used the trivial relation r 2 = zw.  The next result shows how we can compute the Lyapunov and period constants of (26) by using the Lie bracket. This result was previously obtained in [24], but its proof is included here for completeness. Theorem 6 Let us denote by Z the vector field (26), and consider system U defined by the differential equation ż = z + ∞ k=2N +2 k l=0 Ak−l,l z k−l wl . Then the Lie bracket between Z and U has the form ⎛ ⎞ N N [Z, U] = ⎝ p2 j+1 z(zw) j + O2N +2 (z, w), p 2 j+1 w(zw) j + O2N +2 (z, w)⎠ , j=1 j=1 where O2N +2 (z, w) denotes a sum of monomials of degree at least 2N + 2 in z and w, and 1 p2 j+1 + p 2 j+1 V2 j+1 = Re( p2 j+1 ) = = α2 j+1 , (28) 2j 4j 1 p2 j+1 − p 2 j+1 T2 j = Im( p2 j+1 ) = = β2 j+1 , (29) 2j 4j i are the coefficients of the return map and the period function of system (26), respec- tively. Proof Observe that, by Proposition 5, the quantities α2 j+1 and β2 j+1 from (26) are in fact the Lyapunov and period constants, respectively, of the system in this normal form. Let us compute the Lie bracket between Z and U. Actually, we only need to find the first component [Z, U]1 of such operation, since its second component is its complex conjugate. 20 Page 12 of 35 R. D. S. Oliveira et al. ⎛ ⎞  N ∞ k [Z , U ] = ⎝i + 1 (α2 j+1 + i β2 j+1 )( j + 1)(zw) j⎠ z+ Ak−l,l z k−l w l j=1 k=2N +2 l=0 ⎛ ⎞  N ∞ k +⎝ (α2 j+1 + i β2 j+1 ) j z j+1 w j−1 ⎠ w+ Ak−l,l z l w k−l j=1 k=2N +2 l=0 ∞ k−1 ⎛ N ⎞ − 1+ Ak−l,l (k − l)z k−l−1 w l ⎝i z + (α2 j+1 + i β2 j+1 )z(zw) j⎠ k=2N +2 l=0 j=1 ∞ k ⎛ N ⎞ − Ak−l,l lz k−l wl−1 ⎝− i w + (α2 j+1 − i β2 j+1 )w(zw) j ⎠ k=2N +2 l=1 j=1 N N = 2 j(α2 j+1 + i β2 j+1 )z(zw) j + O2N +2 (z, w) =: p2 j+1 z(zw) j + O2N +2 (z, w), j=1 j=1 and the second component is N [Z, U]2 = 2 j(α2 j+1 − i β2 j+1 )w(zw) j + O2N +2 (z, w) j=1 N = p 2 j+1 w(zw) j + O2N +2 (z, w). j=1 Notice that we have denoted by p2 j+1 the coefficient of the Lie bracket with degree 2 j +1, and it has the expression p2 j+1 = 2 j(α2 j+1 +i β2 j+1 ). Finally, as we observed that α2 j+1 and β2 j+1 are the Lyapunov and period constants of system (26), formulas (28) and (29) follow.  Remark 7 It is worth remarking that this Lie bracket method does not allow to obtain the general expression of the Lyapunov and period constants V2 j+1 and T2 j , but only under the conditions V2i+1 = T2i = 0 for every i < j. In this sense, and with a slight abuse of notation, when using the Lie bracket method throughout this paper we will denote by V2 j+1 the jth Lyapunov constant assuming that V2i+1 = T2i = 0 for every i < j, and by T2 j the jth period constant V2i+1 = T2i = 0 for every i < j. We have seen that, when having a system in normal form (26), one can use the Lie bracket to find the corresponding Lyapunov and period constants. The last step is to prove that, in fact, the system does not have to be necessarily in such normal form. In other words, we will show that the method is equally valid if a change of variables is performed on the system, and this will allow to apply the Lie bracket method to general systems of the form (2) and not only to those which are in normal form. A similar approach can be also seen in [2]. Before the main result we present the following lemma. Lemma 8 Let φ = (φ1 , φ2 ) : K2 → K2 be a C 2 -diffeomorphism which maps (u, v) ∈ K2 , where K = R or C, to new variables (x, y) = φ(u, v) = (φ1 (u, v), φ2 (u, v)) ∈ Simultaneous Bifurcation of Limit Cycles and Critical Periods Page 13 of 35 20 K2 , and whose inverse is (u, v) = φ −1 (x, y) = (φ1−1 (x, y), φ2−1 (x, y)). Then the following equivalences hold: ∂φ1−1 ∂φ1 −1 ∂φ −1 ∂φ2 −1 (x, y) · (φ (x, y)) + 1 (x, y) · (φ (x, y)) = 1, (30) ∂x ∂u ∂y ∂u ∂φ2−1 ∂φ1 −1 ∂φ −1 ∂φ2 −1 (x, y) · (φ (x, y)) + 2 (x, y) · (φ (x, y)) = 1, (31) ∂x ∂v ∂y ∂v ∂φ1−1 ∂φ1 −1 ∂φ −1 ∂φ2 −1 (x, y) · (φ (x, y)) + 1 (x, y) · (φ (x, y)) = 0, (32) ∂x ∂v ∂y ∂v ∂φ2−1 ∂φ1 −1 ∂φ −1 ∂φ2 −1 (x, y) · (φ (x, y)) + 2 (x, y) · (φ (x, y)) = 0. (33) ∂x ∂u ∂y ∂u Proof The proof is straightforward by applying the chain rule for two variable func- tions. For (30), we use that φ1 (u, v) = x and φ2 (u, v) = y and the chain rule to rewrite it as ∂φ1−1 ∂ x −1 ∂φ −1 ∂ y −1 ∂φ −1 (x, y) ∂u (x, y) · (φ (x, y))+ 1 (x, y) · (φ (x, y))= 1 = =1. ∂x ∂u ∂y ∂u ∂u ∂u For (31), (32), and (33), we follow an analogous procedure and we obtain that they are equivalent to ∂v/∂v = 1, ∂u/∂v = 0, and ∂v/∂u = 0, respectively.  Theorem 9 Let us consider two vector fields F1 = (F11 , F12 ) and F2 = (F21 , F22 ) in variables (u, v) ∈ K2 , where K = R or C. Let φ = (φ1 , φ2 ) : K2 → K2 be a C 2 -diffeomorphism which maps (u, v) ∈ K2 to new variables (x, y) = φ(u, v) = (φ1 (u, v), φ2 (u, v)) ∈ K2 , and whose inverse is (u, v) = φ −1 (x, y) = (φ1−1 (x, y), φ2−1 (x, y)). Let us denote by G 1 = (G 11 , G 21 ) and G 2 = (G 12 , G 22 ) the vector fields F1 and F2 , respectively, in variables (x, y) after applying the change of variables (x, y) = φ(u, v). Then the following equivalence between the Lie brackets of F1 , F2 and G 1 , G 2 holds: [G 1 , G 2 ]T = Jφ (φ −1 (x, y)) · [F1 , F2 ]φT −1 (x,y) , (34) where Jφ (φ −1 (x, y)) is the Jacobian matrix of φ(u, v) evaluated at (u, v) = φ −1 (x, y), [F1 , F2 ]φ −1 (x,y) denotes the Lie bracket between F1 and F2 also evalu- ated at (u, v) = φ −1 (x, y), and the superindex T denotes the transpose vector. Proof Let us start by observing that equivalence (34) can be rewritten in matrix form as ⎛ ⎞   ∂φ1 −1 ∂φ1 −1  [G 1 , G 2 ] 1 ⎜ ∂u (φ (x, y)) ∂v (φ (x, y))⎟ [F1 , F2 ]1 −1 φ (x,y) = ⎝ ∂φ ∂φ2 −1 ⎠· . (35) [G 1 , G 2 ]2 [F1 , F2 ]2φ −1 (x,y) (φ −1 (x, y)) 2 (φ (x, y)) ∂u ∂v 20 Page 14 of 35 R. D. S. Oliveira et al. We will prove the first component of equivalence (35), since the second one can be analogously seen. Then, our aim is to show that ∂φ1 −1 [G 1 , G 2 ]1 = (φ (x, y)) · [F1 , F2 ]1φ −1 (x,y) ∂u ∂φ1 −1 + (φ (x, y)) · [F1 , F2 ]2φ −1 (x,y) . (36) ∂v The idea of the proof is to develop the expression of [G 1 , G 2 ]1 when performing the change of variables φ and to check that it satisfies (36). First we will see how G 1 and G 2 can be expressed in terms of F1 and F2 . Let us denote by G  j (u, v) the vector i j  ◦ φ −1 )(x, y) for i, j = 1, 2. As vector fields G 1 and field such that G i (x, y) = (G j i i G i2 correspond respectively to ẋ and ẏ, we can take the derivative of (x, y) = φ(u, v) with respect to time and apply the chain rule to see that ∂φ j 1 ∂φ j 2 j = G F + F , (37) i ∂u i ∂v i since vector fields Fi1 and Fi2 correspond respectively to u̇ and v̇. The first partial derivative we need to find for [G 1 , G 2 ]1 is ∂G 11 /∂ x and, by applying the chain rule, it is ∂G 11 1 ◦ φ −1 ) ∂(G ∂G 1 ∂u (x, y) = 1 (x, y) = 1 (φ −1 (x, y)) · ∂x ∂x ∂u ∂x 1 ∂G ∂v + 1 (φ −1 (x, y)) · . ∂v ∂x (38) 1 /∂u by taking the derivative of (37) with respect We can find the expression for ∂ G 1 to u, and we get 1   ∂G ∂ ∂φ1 1 ∂φ1 2 ∂ 2 φ1 1 ∂φ1 ∂ F11 ∂ 2 φ1 2 ∂φ1 ∂ F12 1 = F + F = F + + F + . ∂u ∂u ∂u 1 ∂v 1 ∂u 2 1 ∂u ∂u ∂u∂v 1 ∂v ∂u To find ∂ G1 /∂v we proceed analogously. Now, by substituting these two expressions 1 in (38), we get ∂G 11 /∂ x in terms of the derivatives of F1 and F2 as we wanted. The same procedure can be applied to ∂G 11 /∂ y, ∂G 12 /∂ x, and ∂G 12 /∂ y, and then we have all the terms of [G 1 , G 2 ]1 expressed with F1 and F2 . For the sake of compactness, from now on we will use the following notation to write this expression of [G 1 , G 2 ]1 : ∂φ j −1 ∂φ j −1 ∂ 2 φ j −1 φ ju = (φ (x, y)), φ jv = (φ (x, y)), φ juu = (φ (x, y)), ∂u ∂v ∂u 2 ∂ 2 φ j −1 ∂ 2 φ j −1 ∂ 2 φ j −1 φ juv = (φ (x, y)) = (φ (x, y)), φ jvv = (φ (x, y)). ∂u∂v ∂v∂u ∂v 2 Simultaneous Bifurcation of Limit Cycles and Critical Periods Page 15 of 35 20 Also, with a slight abuse of notation, we will denote the derivatives (∂ Fi /∂u)(φ −1 j (x, y)) and (∂ Fi /∂v)(φ −1 (x, y)) simply as ∂ Fi /∂u and ∂ Fi /∂v, respectively. j j j The expression of [G 1 , G 2 ]1 can then be written as ∂G 11 1 ∂G 11 2 ∂G 12 1 ∂G 12 2 [G 1 , G 2 ]1 = G + G − G − G ∂x 2 ∂y 2 ∂x 1 ∂y 1  ∂ F11 ∂ F12 ∂u = φ1uu F1 + φ1u 1 + φ1uv F1 + φ1v 2 ∂u ∂u ∂ x   ∂ F11 ∂ F12 ∂v   + φ1uv F1 + φ1u 1 + φ1vv F1 + φ1v 2 φ1u F21 + φ1v F22 ∂v ∂v ∂ x  ∂ F11 ∂ F12 ∂u + φ1uu F1 + φ1u 1 + φ1uv F1 + φ1v 2 ∂u ∂u ∂ y   ∂ F11 ∂ F12 ∂v   + φ1uv F1 + φ1u 1 + φ1vv F1 + φ1v 2 φ2u F21 + φ2v F22 ∂v ∂v ∂ y  ∂ F 1 ∂ F 2 ∂u − φ1uu F21 + φ1u 2 + φ1uv F22 + φ1v 2 ∂u ∂u ∂ x   ∂ F21 ∂ F22 ∂v   + φ1uv F2 + φ1u 1 + φ1vv F2 + φ1v 2 φ1u F11 + φ1v F12 ∂v ∂v ∂ x  ∂ F 1 ∂ F 2 ∂u − φ1uu F21 + φ1u 2 + φ1uv F22 + φ1v 2 ∂u ∂u ∂ y   ∂ F21 ∂ F22 ∂v   + φ1uv F2 + φ1u 1 + φ1vv F2 + φ1v 2 φ2u F11 + φ2v F12 . ∂v ∂v ∂ y (39) This expression can be expanded, and a long sum of terms in products of F11 , F12 , F21 , F22 and their derivatives is obtained. Even though we will not show the complete expansion here due to length reasons, we will split it into several parts and see what happens in each case. First, it is straightforward to see that the terms with F11 F21 cancel with each other, as well as those terms with F12 F22 . Now let us focus on the terms with F11 F22 and F12 F21 . The coefficient of F11 F22 in (39) is     ∂v ∂v ∂u ∂u φ1v + φ2v − φ1u + φ2u φ1uv ∂x ∂y ∂x ∂y     ∂u ∂u ∂v ∂v + φ1v + φ2v φ1uu − φ1u + φ2u φ1vv . ∂x ∂y ∂x ∂y 20 Page 16 of 35 R. D. S. Oliveira et al. We can show that this expression is zero by applying Lemma 8. In particular, the coefficient of φ1uv equals zero due to (30) and (31), and the coefficients of φ1uu and φ1vv also equal zero due to (32) and (33), respectively. Therefore, the term with F11 F22 vanishes in (39). Analogously, by applying Lemma 8, we can see that the coefficient of F12 F21 also vanishes. Now let us see how the Lie bracket of F1 and F2 appears in (39). One can see that ∂ F11 1 the coefficient of F in this expression is ∂u 2   ∂u ∂u φ1u φ1u + φ2u = φ1u , ∂x ∂y where we have applied (30) from Lemma 8. Analogously, applying both (30) and (31) ∂ F11 2 ∂ F21 1 ∂ F21 2 from such lemma, we see that the coefficients of F2 , F1 , and F in ∂v ∂u ∂v 1 (39) are respectively φ1u , −φ1u , and −φ1u . Therefore, as these are the terms of the first component of the Lie bracket according to (13), we obtain that those terms altogether ∂ F2 ∂ F2 equal φ1u [F1 , F2 ]1 . The same procedure can be followed with terms 1 F21 , 1 F22 , ∂u ∂v ∂ F22 1 ∂ F22 2 F , and F from (39), and we see that these terms equal φ1v [F1 , F2 ]2 . ∂u 1 ∂v 1 After all this, expression (39) can be rewritten as [G 1 , G 2 ]1 = φ1u [F1 , F2 ]1 + φ1v [F1 , F2 ]2    ∂u ∂u ∂ F11 2 ∂ F12 2 ∂ F21 2 ∂ F22 2 + φ1v + φ2v φ1u F + φ1v F − φ1u F − φ1v F ∂x ∂y ∂u 2 ∂u 2 ∂u 1 ∂u 1    ∂v ∂v ∂ F1 ∂ F2 ∂ F1 ∂ F2 + φ1u + φ2u φ1u 1 F21 + φ1v 1 F21 − φ1u 2 F11 − φ1v 2 F11 . ∂x ∂y ∂v ∂v ∂v ∂v (40) Finally, applying (32) and (33) on the first factor in the second and third lines in (40) and seeing that they are zero, we prove relation (36) for the first component of [G 1 , G 2 ]. To prove the relation for the second component, this is ∂φ2 −1 ∂φ2 −1 [G 1 , G 2 ]2 = (φ (x, y)) · [F1 , F2 ]1φ −1 (x,y) + (φ (x, y)) · [F1 , F2 ]2φ −1 (x,y) , ∂u ∂v we follow an analogous procedure and see that [G 1 , G 2 ]2 = φ2u [F1 , F2 ]1 + φ2v [F1 , F2 ]2 , and the proof finishes.  Simultaneous Bifurcation of Limit Cycles and Critical Periods Page 17 of 35 20 Now let us consider the Lie bracket of two planar real vector fields written (abusing notation) in complex coordinates as (F1 , F 1 ) and (F2 , F 2 ), and the Lie bracket of new vector fields (G 1 , G 1 ) and (G 2 , G 2 ) after a change of variables (φ, φ) on F1 and F2 , respectively. As we already stated, these Lie brackets can be described only from their first components, as the second ones are the complex conjugates, so in such case we will denote the fist components of the Lie bracket simply by [F1 , F2 ] and [G 1 , G 2 ], with a slight abuse of notation as in (14). Therefore, relation (34) from the theorem can be written, in this case, as [G 1 , G 2 ] = φu [F1 , F2 ] + φ u [F1 , F2 ]. With this relation, the following corollary from Theorem 9 is straightforward. Corollary 10 Let (F1 , F 1 ) and (F2 , F 2 ) be two planar complex vector fields which correspond to two real vector fiels in the plane, and let (G 1 , G 1 ) and (G 2 , G 2 ), respectively, be the same vector fields after a change of variables (φ, φ). If [F1 , F2 ] = 0 then [G 1 , G 2 ] = 0. 3 Bi-Weakness for Certain Families This section is devoted to prove Theorem 1. The way to do this is to apply the Lie bracket method introduced in Sect. 2, so all the found bi-weak types will be of the form [k, 2k]. Actually, the proved results are maximal in the sense that if we vanish V2m+1 and T2m for m ≤ k then V2k+3 = 0 and T2k+2 = 0. However, this does not mean that such bi-weak type is the maximal for the system, but the maximal which can be found by means of the Lie bracket method, since such method only detects bi-weak types of the form [k, 2k]. (i) Let us start by studying system (9). We find V3 and T2 by using the Lie bracket method, 1 1 V3 = (2b2 a2 − 3a3 ), T2 = (4a 2 + 10b22 − 9b3 ). (41) 4 12 2 (1) System {V3 = 0, T2 = 0} has solution S1 = {a3 = 23 b2 a2 , b3 = 49 a22 + 10 9 b2 }. The 2 (1) next step is to find V5 and T4 assuming that conditions in S1 hold. This gives 5 5 V5 = − (2b2 a23 + 5b23 a2 ), T4 = − (11a22 b22 − 14b24 ), 27 27 (2) and {V5 = 0, T4 = 0} has the solution S1 = {b2 = 0}. Finally, if we compute V7 and (1) (2) T6 under conditions S1 ∪ S1 , we find V7 = 0 and T6 = 0, and this case corresponds to system  ẋ = −y + a2 x 2 , (42) ẏ = x + 49 a22 x 3 , 20 Page 18 of 35 R. D. S. Oliveira et al. which is an isochronous center as we will see in Proposition 12. System (9) under (1) condition S1 has then a bi-weak [2, 4] type, which is the maximal of the form [k, 2k] as we have seen that [3, 6] does not appear. (ii) For the quartic Liénard (10), we proceed analogously. By using the Lie bracket technique we find that V3 and T2 are the same that in the cubic case, (41), so they vanish again for S2(1) = {a3 = 23 b2 a2 , b3 = 49 a22 + 10 9 b2 }. Now we find the next 2 constants, which are 1 V5 = (−20a23 b2 − 50a2 b23 + 27a2 b4 + 90a4 b2 ), 54 1 T4 = (−110a22 b22 − 140b24 + 72a2 a4 + 189b2 b4 ), 54 and {V5 = 0, T4 = 0} has two solutions   (2a) 5 b2 a2 (a 2 + 7b22 ) 10 b2 (8a24 − 35a22 b22 − 70b24 ) S2 = a4 = − 2 2 2 , b4 = 3 4a2 − 35b22 27 4a22 − 35b22 (2b) and S2 = {a2 = 0, b2 = 0}. The next step is to find V7 and T6 , which under conditions S2(1) ∪ S2(2a) are (a) 35 b23 a2 (8a26 + 525a24 b22 − 2520a22 b24 − 4900b26 ) V7 =− , 216 (4a22 − 35b22 )2 (a) 35 b24 (430a26 − 3900a24 b22 + 46431a22 b24 + 56350b26 ) T6 = , 648 (4a22 − 35b22 )2 (1) (2b) and under S2 ∪ S2 , (b) 21 (b) 1 V7 = a4 b4 , T6 = (84a42 + 189b42 ). 16 80 (a) (a) (3a) (b) System {V7 = 0, T6 = 0} has the solution S2 = {b2 = 0}, and {V7 = (b) (3b) 0, T6 = 0} has the solution S2 = {a4 = 0, b4 = 0}. Finally, if we substitute (1) (2a) (3a) (1) (2b) (3b) S2 ∪ S2 ∪ S2 and S2 ∪ S2 ∪ S2 in V9 and T8 we obtain V9 = 0 and T8 = 0 in both cases, so the maximal bi-weak type is of the form [k, 2k] is [3, 6] and this (1) (2a) (3a) proves Theorem 1 (ii). We notice that case S2 ∪ S2 ∪ S2 corresponds to system (1) (2b) (3b) (42) and case S2 ∪ S2 ∪ S2 is the linear center, so both systems are isochronous centers. (iii) Let us study the quadratic case. First, for the sake of simplicity in the obtained expressions, we will consider system (11) in complex coordinates as  ż = i z + r20 z 2 + r11 zw + r02 w 2 , ẇ = − i w + s02 z 2 + s11 zw + s20 w 2 , Simultaneous Bifurcation of Limit Cycles and Critical Periods Page 19 of 35 20 being si j = r i j . The first Lyapunov and period constants are 4 V3 = i(−r11r20 + s11 s20 ), T2 = r02 s02 + 2r11 s11 − r11r20 − s11 s20 . 3 In this case, solving the systems for Lyapunov and period constants equal to zero as we did in Liénard families is cumbersome, as many more solutions are obtained. For this reason, we will study Lyapunov and period constants V2k+1 and T2k in the Bautin type ideal Bk−1 := V3 , T2 , . . . , V2k−1 , T2k−2 . First we want to check whether V5 and T4 belong to B1 = V3 , T2 . The expression of V5 and T4 in B1 is 1 V5 = 2 i(−4r02 r20 s11 + 6r02 r20 s11 2 + 4r02 s11 3 − 4r11 3 s02 − 6r11 2 s02 s20 + 4r11 s02 s20 2 ), 3 1 T4 = −8r02 r20 2 s11 + 36r02 r20 s11 2 − 40r02 s113 − 40r113 s02 − 96r11 2 2 r20 6  +219r11 2 r20 s11 + 36r11 2 s02 s20 − 135r11 2 2 s11 + 12r11r20 2 s20 − 8r11 s02 s20 2 . As these constants are not identically zero, we have a bi-weak [2, 4] type. Let us check now that [3, 6] does not appear, so we find the expressions of V7 and T6 in B2 = V3 , T2 , V5 , T4 , 1  V7 = i −80r02 r11r20 2 2 s11 + 360r02 r11r20 s11 3 − 100r02 r11 s11 4 − 300r11 4 s02 s11 64  +60r11 r20 − 480r11 3 3 r20 s11 + 1095r11 3 2 3 2 r20 s11 − 675r11 s11 , 3 3 1  T6 = 3032r02 r11r20 s11 − 4548r02 r11r20 s11 2 2 3 − 2282r02 r11 s11 4 480 + 750r11 4 s02 s11 + 94178r11 r20 − 249834r11 3 3 r20 s11 + 287559r11 3 2 3 2 r20 s11  −117863r11 s11 − 14832r11 3 3 r20 s20 + 792r11 r20 2 3 s20 . 3 2 However, we can see that actually V72 = 0 in B2 , i.e. V7 = 0 in the variety of zeros defined by the elements in B2 , so the bi-weak [3, 6] type cannot appear as we have a center and in this case [2, 4] is maximal. We notice that there are quadratic centers with T6 = 0, which makes sense as it is a classical result that quadratic centers can have at most 2 oscillations in the period function (see [7]), because in the center case only T2 , T4 , and T6 are necessary to characterize the isochronicity property. In fact, T82 = 0 in B3 = V3 , T2 , V5 , T4 , T6 . Let us make the following observation regarding the quadratic case. It is a well- known fact that the general quadratic family (11) classically unfolds 3 limit cycles, hence the case V3 = V5 = 0 with nonvanishing V7 exists. This is not the case for the simultaneous study of Lyapunov and period constants we have seen here, but as we are adding the conditions for vanishing the corresponding period constants, the fact that the obtained Lyapunov constants vanish at lower orders is not inconsistent. 20 Page 20 of 35 R. D. S. Oliveira et al. (iv) Finally, let us study the bi-weakness of cubic homogeneous nonlinearity family (12). To simplify the expressions we will perform a change to complex coordinates, obtaining  ż = i z + r30 z 3 + r21 z 2 w + r12 zw 2 + r03 w 3 , ẇ = − i w + s03 z 3 + s12 z 2 w + s21 zw 2 + s30 w 3 , being si j = r i j . We will proceed analogously to the quadratic case, by studying V2k+1 and T2k in the Bautin ideal Bk−1 := V2m+1 , T2m m≤k−1 . The first Lyapunov and period constants are V3 = −r21 − s21 , T2 = i(r21 − s21 ). The next step is to find the expressions of V5 and T4 simplified with respect to B1 = V3 , T2 , V5 = −2 i(r12 r30 − s12 s30 ), T4 = 3r03 s03 − 2r12 r30 + 4r12 s12 − 2s12 s30 . Being in B2 = V3 , T2 , V5 , T4  we have 3 V7 = (3r03r30 2 − 8r03r30 s12 − 3r03 s12 2 − 3r12 2 s03 − 8r12 s03 s30 + 3s03 s30 2 ), 8 3 T6 = − i(3r03r30 2 − 16r03r30 s12 + 21r03 s122 − 21r12 2 s03 + 16r12 s03 s30 − 3s03 s30 2 ). 8 As they are not identically zero, this proves the existence of a bi-weak [3, 6] type. The proof finishes due to the fact that V9 = 0 in B3 = V3 , T2 , V5 , T4 , V7 , T6 . 4 Some Results on Liénard Families This section is devoted to the study of Liénard systems and it is divided into three parts. First, we use the Lie bracket method to deduce the structure of Lyapunov and period constants of a Liénard system starting with an odd and an even degree monomials on its first differential equation. Second, we classify the centers and the isochronicity of a Liénard family. Finally, we use the previous results to provide a complete study of the simultaneous bifurcation of limit cycles and critical periods of the cubic Liénard family, which proves Theorem 2. 4.1 A Liénard Family Starting with an Odd and an Even Degree Monomials In this subsection we will consider the Liénard family  ẋ = −y + am x m + an x n + x d P(x), (43) ẏ = x, Simultaneous Bifurcation of Limit Cycles and Critical Periods Page 21 of 35 20 where m and n are even and odd natural numbers, respectively, d = max(m, n) + 1, and P(x) is a polynomial in x. In order to motivate the problem, let us start by considering (43) being Pd (x) = 0. It is a well-known fact that in this classical Liénard family the coefficients corresponding to odd powers control the center property, so if an = 0 we have a center at the origin. In this case the even power controls the isochronicity property, so if am = 0 the system has an isochronous center. A study in this line is presented for example in [10], but in our case we will present our result by using the Lie bracket method introduced in Subsect. 2.3 to find the Lyapunov and period constants. It is worth recalling that different methods can lead to Lyapunov and period constants that differ in a multiplicative constant, but the dependence on parameters am and an is the same and the center and linearizability conditions are also kept. The result is as follows. Theorem 11 For system (43), (i) if n < 2m − 1, then the first nonidentically zero Lyapunov constant is   an n Vn = (n − 1) n−1 , (44) 2n 2 and the system vanishes all its period constants up to order n; (ii) if n = 2m − 1, then the first nonidentically zero Lyapunov and period constants are (44) and   2 am m 2 (m − 1) 2m T2m−2 = − , (45) 22m−1 (m + 1)2 m respectively; (iii) if n > 2m − 1, then the first nonidentically zero period constant is (45), and the system vanishes all its Lyapunov constants constants up to order 2m − 1. Proof Let us start by writing system (43) in complex coordinates. By using the Bino- mial Theorem, ⎧ m m  n n   z+w d ⎪ ⎪ am an ⎨ż = i z + 2m j z m− j w j + 2n k z n−k w k + 2 P(z, w) =: Z (z, w), j=0 k=0 m m  n n   z+w d (46) ⎪ ⎪ am an ⎩ẇ = −iw + 2m j z j w m− j + 2n k z k w n−k + 2 P d (z, w) =: Z (z, w). j=0 k=0 We will consider the system ⎧ m n ⎪ ⎪ ⎨ż = z + (Am−l,l + i Bm−l,l )z m−l wl + (An− p, p + i Bn− p, p )z n− p w p =: U (z, w), l=0 p=0 m n (47) ⎪ ⎪ ⎩ẇ = w + (Am−l,l − i Bm−l,l )z l w m−l + (An− p, p − i Bn− p, p )z p w n− p =: U (z, w). l=0 p=0 To calculate the Lyapunov and period constants we will find the structure of the Lie bracket of (46) and (47). Observe that, due to the fact that they are associated to a 20 Page 22 of 35 R. D. S. Oliveira et al. real vector field, by using (14) their Lie bracket can be described only from its first component: ∂Z ∂Z ∂U ∂U [Z , U ]1 = U+ U− Z− Z = Hm + Hn + H2m−1 + H2n−1 + Om+n−1 (z, w), ∂z ∂w ∂z ∂w (48) where each Hq is an homogeneous qth degree polynomial in z, w. These polynomials have been found, but they are not written here for the sake of brevity. We have that m    am m Hm = (−Bm−l,l + i Am−l,l )(2l − m + 1) + m (m − 1) z m−l wl . (49) 2 l l=0 This homogeneous mth degree part vanishes taking   am m − 1 m Am−l,l = 0, Bm−l,l = for l ∈ {0, . . . , m}. (50) 2m 2l − m + 1 l The homogeneous nth degree polynomial Hn is analogous to (49), only changing m by n. We can take then An− p, p = 0 for p ∈ {0, . . . , n},     an n−1 n n−1 (51) Bn− p, p = n for p ∈ {0, . . . , n}\ , 2 2p − n + 1 p 2 as for the term corresponding to p = n−1 2 we have that 2 p − n + 1 = 0, and therefore n+1 n−1 n−1 the coefficient of z w 2 = z(zw) 2 cannot be vanished; in this case we take 2 B n+1 , n−1 = 0. Notice that this fact occurs with n because it is an odd number, but not 2 2 with m which is even. Let us see the structure of H2m−1 . After substituting (50), we obtain m     2 am m−1 m m H2m−1 = i 2m (l − j)z 2m−l− j−1 w j+l 2 2l − m + 1 l j j,l=0 (52)  − jz w m+l− j m+ j−l−1 − lz m+ j−l w m+l− j−1 . We know by Theorem 6 that the period constant associated to (52) is the coefficient of z(zw)m−1 . Then, after some basic combinatorial computations, this coefficient can be written as ⎛ ⎞ 2 m−1    m  2 am m m 2 j m p2m−1 = i 2m (m − 1) ⎝ − ⎠ . (53) 2 j j +1 2j − m + 1 j j=0 j=0 Simultaneous Bifurcation of Limit Cycles and Critical Periods Page 23 of 35 20 Let us write this expression in a more compact way. We can use the generating functions (1 + x)m and (1 + x)2m to rewrite the first summation in (53). In the relation 2m   m   2 m    2m k  2 m k m m k+ j x = (1 + x)2m = (1 + x)m = x = x , k k k j k=0 k=0 k, j=0 we equate the coefficients of x m−1 , and we obtain m−1        m m 2m m 2m = = . (54) j j +1 m−1 m+1 m j=0 For the second summation in (53), the equality m  2   2j m m(3m + 1) 2m = (55) 2j − m + 1 j (m + 1)2 m j=0 holds. This equality has been obtained with a computer algebra system by means of the Zeilberger’s algorithm. For more details on how to use such method, the reader is referred to [34,35], as well as [23] and the references therein. Notice that this Zeilberger’s algorithm can also be used to justify (54). Adding (54) and (55), one can rewrite (53) as   2 am m 2 (m − 1) 2m p2m−1 = − i 2m−1 . (56) 2 (m + 1)2 m After these calculations we can finally prove our result. To this end, we will consider n−1 (48) assuming (50) and (51), so Hm = 0 and Hn only has the z(zw) 2 term. If n < 2m − 1, then the lowest degree homogeneous polynomial is the nth degree one,   an n n−1 Hn = (n − 1) n−1 z(zw) 2 , 2n 2 and due to Theorem 6 this coefficient is the first nonidentically zero Lyapunov constant (44), while all period constants vanish up to this order. If n = 2m − 1, then the lowest degree homogeneous polynomial in (48) is Hn + n−1 H2m−1 , so the coefficient of z(zw) 2 = z(zw)m−1 is (44) plus (56). Therefore, by Theorem 6 the real part is the Lyapunov constant (44) and the imaginary part is the first nonidentically zero period constant (45). Finally, for n > 2m − 1, an analogous procedure shows that the lowest degree homogeneous polynomial in (48) is H2m−1 , so the period constant is the imaginary part of (56) and all Lyapunov constants vanish up to this order.  20 Page 24 of 35 R. D. S. Oliveira et al. 4.2 Center and Isochronicity Classification of a Liénard Family In this subsection we present, as an application of our approach, a characterization of isochronous centers in a Liénard family. For a general proof we refer the reader to the very recent work of Amel’kin ([4]). We start with the following result. Proposition 12 Let us consider the Liénard system  ẋ = −y + a2 x 2 + a3 x 3 + a4 x 4 + a5 x 5 + a6 x 6 , (57) ẏ = x + b2 x 2 + b3 x 3 + b4 x 4 + b5 x 5 + b6 x 6 . The only isochronous center having the form (57) is (42). Proof We start by finding the first five Lyapunov and period constants of (57) by using the Lie bracket method from Subsect. 2.3. Then we solve the resulting system {V3 = 0, V5 = 0, V7 = 0, V9 = 0, V11 = 0, T2 = 0, T4 = 0, T6 = 0, T8 = 0, T10 = 0}, and we compute that the only nontrivial solution is {a3 = 0, a4 = 0, a5 = 0, b2 = 0, b3 = 49 a22 , b4 = 0, b5 = 0}, which corresponds to system (42). Such system has a center at the origin due to Theorem 1 from [9], which classifies a class of Liénard systems known as composition centers; this center classification can also be found in [17]. To show the isochronicity of (42) we will start by rescaling the system via the change (x, y) → ( 2a32 x, 2a32 y) for a2 = 0, which results in  ẋ = −y + 23 x 2 , (58) ẏ = x + x 3 . Observe that for the case a2 = 0, the system is isochronous because it corresponds to the linear center. System (58) is isochronous because it commutes with  ẋ = x A(x, y),   ẏ = y + 21 x 2 A(x, y), where A(x, y) = −1 − y + 21 x 2 , as the Lie bracket between both systems is 0. We will present an alternative proof for the isochronicity of (58), because we consider that it is interesting as it simply uses a first integral and the system itself. A first integral of (58) is x 4 − 4x 2 y + 4x 2 + 4y 2 H (x, y) = , (59) (x 2 − 2y − 2)2 which satisfies that H (0, 0) = 0. The idea is to prove the isochronicity of (58) by using (59) to find an expression for the level curves γh and check that its integral through a whole loop is 2π. Simultaneous Bifurcation of Limit Cycles and Critical Periods Page 25 of 35 20 Let us consider H (x, y) = h 2 with h > 0 and solve it for x. As it has degree 4 in x we obtain four solutions, two of which are imaginary for values of h close to 0. The other two solutions are ! −2(h 2 − 1)(−h 2 y − h 2 + 2h 2 y + 3h 2 − 2y + 1 + y − 1) X ± (y, h) = ± , h2 − 1 which are real for 0 < h < 1, and they correspond to the level curve on the right and left side of the vertical axis. By solving H (0, y) = h 2 with respect to y, we can find that the intersections of the level curves with the vertical axis for level h 2 are y± = h/(1 ± h), being 0 < h < 1. Considering the second equation in (58), we aim to calculate the integral for the time  T   y+ dy dy T = dt = =2 =: 2I (h), 0 γh ẏ y− X + (y, h) + X + 3 (y, h) and see that we obtain 2π. Here we have used that, due to the symmetry of the problem, we can simply integrate from y− to y+ and check that we obtain I (h) = π, which represents half a loop. To simplify the expression of I (h) we will consider the change of variables z 2 = 2 −z 2 +1 2h 2 y + 3h 2 − 2y + 1, or equivalently y = 3h2(1−h 2 ) . After applying this change both to X + (y, h) and the integration limits y± , the integral becomes  1−h 1 − h2 I (h) = − √ dz. 1+h (z − 2) −h 4 + h 2 z 2 − 2h 2 z + 2h 2 − z 2 + 2z − 1 This integral can be simplified a little more so that the integration limits become ±1. Let us apply the change z = 1 − hw to a new variable w, and we have  √  √ 1 1 − h2 k 1 − h2 I (h) = √ dw = lim √ dw −11 − w 2 (1 + hw) k→1 −k 1 − w 2 (1 + hw)  w=k h+w π  π = lim arctan √ √ = − − = π, k→1 1 − h 2 1 − w 2 w=−k 2 2 for 0 < h < 1, so the result follows.  By studying the Liénard system from the previous result we have come across the Liénard family  ẋ = −y + ax n , (60) ẏ = x + bx 2n−1 , for a, b ∈ R, n ∈ N, being n ≥ 2. Indeed, (42) is a particular case of (60). The centers and isochronicity of this family are classified on Theorem 14. For the proof of this classification result, we first need the following proposition. 20 Page 26 of 35 R. D. S. Oliveira et al. n2 Proposition 13 System (60) with even n, b = (n+1)2 a 2 , and a = 0 has a first integral of the form A(x, y) H (x, y) =   n−1 , (61) x 2 + y 2 − 2x n y + x 2n 2 being A(x, y) a polynomial in x, y. n2 Proof For b = (n+1)2 a2 and a = 0, system (60) can be rewritten as  ẋ = −y + (n + 1)x n =: P(x, y), (62) ẏ = x + n 2 x 2n−1 =: Q(x, y),    n+1  n−1 1  n+1  n−1 1 after the rescaling (x, y) → a x, a y . It can be checked that func- tion F(x, y) = x 2 + y 2 −2x n y + x 2n is an invariant curve of system (62) with cofactor K = 2nx n−1 , and the divergence of the vector field is div(P, Q) = n(n + 1)x n−1 . Now according to Darboux integrability theory (see for instance [12]), for λ satisfying λK = − div(P, Q) we have that R(x, y) = F(x, y)λ is an integrating factor of the system. In our case, λK = − div(P, Q) is λ2nx n−1 = −n(n + 1)x n−1 , so λ = − n+1 2  2 − n+1 and R(x, y) = x + y − 2x y + x 2 n 2n 2 is an integrating factor. Having an integrating factor R(x, y) of the system, we know that a first integral H (x, y) satisfies ∂∂Hx = Q(x, y)R(x, y) and ∂∂Hy = −P(x, y)R(x, y), so we can integrate this second equation to find the form of such first integral,   y − nx n − x n H (x, y) = −P(x, y)R(x, y)dy =   n+1 dy x 2 + y 2 − 2x n y + x 2n 2   y − xn xn =   n+1 dy − n   n+1 dy. x 2 + y 2 − 2x n y + x 2n 2 x 2 + y 2 − 2x n y + x 2n 2 (63) −1 The first term in (63) is an immediate integral and its result is n−1 , (n−1)(x 2 +y 2 −2x n y+x 2n ) 2 taking into account that due to being a first integral we can consider that the integra- tion constant is 0. For the second integral we will perform the change of variables ω= √ 2 2x n , so dy = − 2 √x 2 dω and the integral can be written as x +y −2x y+x2n ω 1−ω   xn ωn−1   n+1 dy = − √ 1 − ω2 dω. x 2 + y 2 − 2x n y + x 2n 2 Simultaneous Bifurcation of Limit Cycles and Critical Periods Page 27 of 35 20 Now we can apply a trigonometric change of variables ω = sin φ, and    ωn−1 sinn−1 φ − √ dω = − ! cos φdφ = − sinn−1 φdφ 1 − ω2 1 − sin2 φ  =− sin φ sinn−2 φdφ. As n is even, so is n − 2 and the integral becomes    n−2    n−2 2 2 − sin φ sin φ 2 dφ = − sin φ 1 − cos2 φ dφ ⎛ ⎞  n−2 2 ⎜ ⎟ = − sin φ ⎝ a j cos2 j φ ⎠ dφ, j=0 for certain coefficients a j , where in the last equality we have used the Binomial Theorem. Then, after swapping the sum and the integral we get n−2 2  n−2 2 n−2 2 aj aj − aj sin φ cos2 j φdφ = cos2 j+1 φ = cos φ cos2 j φ. 2j + 1 2j + 1 j=0 j=0 j=0  j By using cos2 j φ = 1 − sin2 φ and applying the Binomial Theorem again, we have that for new coefficients b j , n−2 n−2 n−2 2 aj 2 ! 2 cos φ cos2 j φ = cos φ b j sin 2j φ= 1 − ω2 b j ω2 j , 2j + 1 j=0 j=0 j=0 after undoing the change ω = sin φ. Finally, we can substitute ω = √ x , x 2 +y 2 −2x n y+x 2n and after some trivial calculations we get n−2 n−2 n−2 ! (y − x n ) 2 2 j (x 2 + y 2 − 2x n y + x 2n ) 2 −j j=0 b j x 2 1 − ω2 bjω 2j = n−1 , j=0 (x 2 + y 2 − 2x n y + x 2n ) 2 so joining both terms from (63) the result follows.  Theorem 14 The Liénard family (60) satisfies that, (i) for odd n, the system is a center if and only if a = 0 –in which case it is a reversible center–, and it is an isochronous center if and only if a = b = 0; (ii) for even n, the system is a reversible center for any a and b, and it is an isochronous n2 center if and only if b = (n+1) 2a . 2 20 Page 28 of 35 R. D. S. Oliveira et al. Proof Let us start with the case of n being odd. If a = 0, it is immediate to see that (60) is a time-reversible center. Consider the function 1 2 1 2 b 2n H (x, y) = x + y + x , (64) 2 2 2n and compute       ∇ H · (ẋ, ẏ) = x + bx 2n−1 , y · −y + ax n , x + bx 2n−1 = ax n+1 1 + bx 2n−2 . We can see then that for a = 0 function H is a first integral of the system, and for a = 0 it is a Lyapunov function of the system. In the latter case, the origin is a focus and the sign of a determines its stability near the origin. We will prove now the isochronicity condition provided that we have a center, i.e. a = 0. If b = 0, (60) becomes the linear center so it is isochronous. Now let us show that if b = 0 then the system is not isochronous, for which we will integrate the second equation in (60) along the level curves γh determined by the first integral (64). The integral to solve is  T   dy dy T = dt = = . (65) 0 γh ẏ γh x + bx 2n−1 2 By considering (64) on a level curve such that H (x, y) = h2 being h > 0, we obtain that y 2 = h 2 −x 2 − nb x 2n , and we can use this relation to perform a change of variables from y to x in (65) as follows:   x h+  xh dy dx dx T = =2 =4 , γh x + bx 2n−1 x h− h 2 − x 2 − nb x 2n 0 h 2 − x 2 − nb x 2n (66) where we have used the symmetry of the integral and x h− , x h+ = x h , are the intersec- 2 tions of the level curve with the horizontal axis, i.e. the real solutions of H (x, 0) = h2 or equivalently x 2 + nb x 2n − h 2 = 0. As x h depends on h, we will perform a second change of variables x = x h z so that the integration limits are constant, so we rewrite (66) as  1 xh T =4 dz. (67) 0 h2 − xh 2 z2 − nb x h 2n z 2n As we are not able to find the explicit solution of L(h) := x 2 + nb x 2n − h 2 = 0 to have an expression of x h , we consider a power series expansion of x h with respect to h. One can check that such series expansion starts Simultaneous Bifurcation of Limit Cycles and Critical Periods Page 29 of 35 20 b 2n−1 xh = h − h + O4n−3 (h), (68) 2n as for this x h we have  2  2n b 2n−1 b b 2n−1 L(h) = h − h + O4n−3 (h) + h− h + O4n−3 (h) − h2 2n n 2n b b  2n  = h 2 − h 2n + O4n−2 (h) + h + O4n−2 (h) − h 2 = O4n−2 (h). n n Now by substituting (68) in (67) and developing the expression we have  1 1− b 2n−2 2n h + O4n−4 (h) T =4   dz, 0 1 − z2 + b n z 2 − z 2n h 2n−2 + O4n−4 (h) where we can expand the denominator as a power series and obtain  1 b 2n−2  1 b z 2 − z 2n 2n−2  T =4 1−h + O4n−4 (h) √ − h + O4n−4 (h) dz 0 2n 1 − z2 2n (1 − z )2 3/2  1    1 b 1 z 2 − z 2n =4 √ − √ + h 2n−2 + O 4n−4 (h) dz 0 1 − z2 2n 1 − z2 (1 − z 2 )3/2  1  1 b 1 + z 2 + z 4 + · · · + z 2n−2 2n−2 =4 √ − √ h + O4n−4 (h) dz. 0 1−z 2 2n 1 − z2 Notice that on the last equality we have used the formula for the sum of terms in a geometric progression. Finally, this integral becomes   1 2b 2n−2 1 1 + z 2 + z 4 + · · · + z 2n−2 dz T =4 √ − h √ dz + O4n−4 (h) 0 1 − z2 n 0 1 − z2  1 2b 1 + z 2 + z 4 + · · · + z 2n−2 = 2π − h 2n−2 √ dz + O4n−4 (h). n 0 1 − z2 As the remaining integral has a strictly positive integrand in the considered interval, the integral is strictly positive and nonzero, so if b = 0 the term of order 2n − 2 in h is nonzero and the center is not isochronous, hence the isochronicity condition is proved. For the case of n being even, it is straightforward to see that (60) is a time-reversible center for any value of the parameters a and b. To prove the isochronicity condition, n2 we will see that if b = (n+1) 2 2 a then there exists a transversal system such that its Lie bracket with the original system vanishes, and for the reciprocal we will check that if n2 b = (n+1) 2 2 a then the period function is not constant. 2 Let us assume that b = (n+1) n 2 a , and consider the rescaled system (62). According 2 to Proposition 13, this system has a first integral of the form (61). Observe that, due 20 Page 30 of 35 R. D. S. Oliveira et al. to the fact of being a first integral, the condition ∂∂Hx P + ∂∂Hy Q = 0 must be satisfied. In our case, when finding the partial derivatives of the first integral (61) we obtain ∂H ∂H    ∂A P+ Q = − n(n − 1)x n A(x, y) + x y − (n + 1)x n+1 ∂x ∂y ∂x   ∂A − x 2 + n 2 x 2n F(x, y), ∂y where F(x, y) = x 2 + y 2 − 2x n y + x 2n is the invariant curve of the system found in the proof of Proposition 13. Therefore, the condition that must be fulfilled is   ∂A   ∂A n(n − 1)x n A(x, y) + x y − (n + 1)x n+1 − x 2 + n 2 x 2n = 0. (69) ∂x ∂y Let us consider system  ẋ = x A(x, y), (70) ẏ = (y + (n − 1)x n ) A(x, y), where A(x, y) is the numerator of the first integral (61) of the system. A straightforward computation shows that the Lie bracket between systems (62) and (70) is exactly the left hand side of equality (69), so by construction it equals 0 and the system is isochronous due to Theorem 3. n2 Finally, we have to check that if b = (n+1) 2 a , then the center is not isochronous. 2 To this end, we will use the Lie bracket method to find the first period constant and n2 check that it only vanishes for b = (n+1) 2 a . System (60) can be rewritten in complex 2 coordinates as  n  2n−1 z+w z+w ż = i z + a +ib , (71) 2 2 which is actually system (46) choosing am = a, an = i b, and P(z, w) = 0, and taking into account that (m, n) in (46) corresponds to (n, 2n − 1) in (71). By applying the Lie bracket method analogously to the proof of Theorem 11, one can see that the first nonzero term is that of degree 2n − 1, whose coefficient is     b 2n − 1 a 2 n 2 (n − 1) 2n p2n−1 = i (2n − 2) − i 2n−1 . (72) 2 2n−1 n−1 2 (n + 1)2 n By Theorem 6, the imaginary part of (72) is the period constant         b 2n a 2 n 2 (n − 1) 2n n2 2 n − 1 2n T2n−2 = (n − 1) − 2n−1 = b − a , 2 2n−1 n 2 (n + 1)2 n (n + 1)2 22n−1 n Simultaneous Bifurcation of Limit Cycles and Critical Periods Page 31 of 35 20 n2 which only vanishes for b = (n+1)2 a2.  Observation 15 After presenting the new isochronous family (60) for even n and n2 b = (n+1) 2 a , a natural question is to wonder if these isochronous centers will 2 provide a high criticality under perturbation –notice that we do not consider odd n because in this case the only isochronous center is the linear one. To inquire into this question, we find the number of critical periods which can bifurcate from the system by considering only linear parts in the perturbative parameters of the period constants, following the ideas presented in previous works [31,32]. In particular, we have performed this study for n = 2, 4, 6, 8, which correspond to systems of degrees 2n − 1 = 3, 7, 11, 15 respectively, and by adding either a time-reversible perturba- tion with respect to the horizontal axis (x, y, t) → (x, −y, −t) or with respect to the vertical axis (x, y, t) → (−x, y, −t), none of which breaks the center property. After computing period constants up to first order and calculating the corresponding ranks, we find the number of critical periods presented in the next chart. We also show the number of critical periods (n 2 + n − 4)/2 obtained in [32] for a class of holomorphic systems also with linear parts, which were found with a reversible perturbation with respect to the horizontal axis. n Degree Critical periods Critical periods n 2 +n−4 2 (reversible respect (reversible respect horizontal axis) vertical axis) 2 3 4 2 4 4 7 20 13 26 6 11 44 30 64 8 15 76 53 118 As we can check, the obtained criticality of the isochronous family (60) for even n found up to linear parts in the period constants is much worse than the one obtained for the holomorphic family in [32] also with linear parts. This leads us believe that family (60) will not provide a high number of oscillations of the period function, so we will not go further in the study of its criticality. 4.3 Limit Cycles and Critical Periods in the Cubic Liénard Family The aim of this subsection is to prove Theorem 2, this is to study the simultaneous bifurcation of limit cycles and critical periods of the cubic Liénard system (9) adding the trace parameter α as in (1). With the tools introduced in Sect. 2 we can find the first Lyapunov and period constants for system (9), assuming that the previous ones vanish: 20 Page 32 of 35 R. D. S. Oliveira et al. π 5π V3 = − (2a2 b2 − 3a3 ), V5 = a2 b2 b3 , 4 12 π π T2 = (4a22 + 10b22 − 9b3 ), T3 = a2 (2a2 b2 − 3a3 ), 12 12 π π T4 = − (3484a22 b22 + 4480b24 + 81a32 ), T5 = − b2 (10a24 + 3193a22 b22 + 4032b24 ). 3456 162 (73) Let us make the following observation. As we already saw for the quadratic system (8), in a system not having a center at the origin it is not generally true that its first period constant must have even subindex. The presented cubic Liénard system is another example of this fact. If we compare V3 and T3 we can observe this: for a center we would have V3 = 0 which means 2a2 b2 − 3a3 = 0, and this implies T3 = 0, but T3 can be nonzero if V3 = 0. We can also see by comparing V5 and T5 that this equality in the factors of V3 and T3 is not a general fact for any pair V2k+1 and T2k+1 . Proof of Theorem 2 (i) Let us start by finding the centers in family (9). To do this we solve the system formed by the two first Lyapunov constants of (9) being equal to 0, this is {V3 = 0, V5 = 0}, and we obtain three solutions: (ẋ, ẏ) = (−y, x + b2 x 2 + b3 x 3 ), (74) (ẋ, ẏ) = (−y + a2 x , x + b3 x ), 2 3 (75)   2 (ẋ, ẏ) = −y + a2 x + a2 x , x + b2 x . 2 3 2 (76) 3 Equations (74) and (75) are time-reversible families with respect to the x and y axis, respectively, and (76) corresponds to a center due to Theorem 1 from [9]. Observe that such theorem is also an alternative proof for (75) having a center at the origin. Therefore, we have seen that for (9) only two Lyapunov constants are necessary to solve the center problem, this is V3 = V5 = 0 implies V2k+1 = 0 for any k ≥ 3, which is actually a well-known fact. Now we will check that (74) unfolds  one critical period. To see this, we find from (73) that for this family T2 = 12 1 π 10b22 − 9b3 , T3 = 0 –as it is a center–, and T4 = − 35 27 π b2 . Then, if b3 = 9 b2 we see that T2 = 0 but T4 can be different from 0, 4 10 2 so 1 critical period unfolds. Now if we vanish T2 and T4 , automatically b2 = b3 = 0 and (74) becomes the linear center, so it is isochronous. We will also see that for (75) and (76), when vanishing T2 an isochronous system is obtained so no critical periods appear, thus  the critical  period obtained from (74) is maximal. For (75) we have that T2 = 12 1 π 4a22 − 9b3 , and it vanishes for b3 = 49 a22 , which corresponds  to system  (42) and it is isochronous due to Proposition 12. For (75), T2 = 12 1 π 4a22 + 10b22 , which only vanishes for a2 = b2 = 0 and in this case the system becomes the linear center, hence is isochronous. This finishes the proof of the statement. (ii) The first part of the statement follows by using the fact that the first two Lyapunov constants in (73) characterize the centers for the Liénard system given in (9), they are independent and the Bautin ideal V3 , V5  is radical. For more details see [29]. Simultaneous Bifurcation of Limit Cycles and Critical Periods Page 33 of 35 20 Taking again (73), the simultaneous bifurcation will follow by looking for a bi-weak equilibrium point of type [2, 4], that is, solving {α = V3 = T2 = T3 = 0} with V5 and T4 nonvanishing. Under these conditions, we get the set S = {a3 = 23 a2 b2 , b3 = 9 a2 + 9 b2 , a2 b2  = 0} and we can check that on S we have 4 2 10 2 5 5 V5 = πa2 b2 (2a22 + 5b22 ), T4 = − π b22 (11a22 + 14b22 ). 54 54 We take α = 0 and (a2 , a3 , b2 , b3 ) → (−3 + εe1 , 6 + εe2 , −3 + εe3 , 14 + εe4 ). Clearly, when ε = 0 we have a bi-weak equilibrium point of type [2, 4]. It can be checked that the linear parts of V3 , T2 , and T3 with respect to ε are V3[1] = T3[1] = 3π ε(2e1 +e2 +2e3 )/4 and T2[1] = −π ε(8e1 +20e3 +3e4 )/4, respectively. Therefore, as they are independent, we have one limit cycle and two critical periods of small amplitude. The proof of this statement finishes by adding the parameter α and using Lemma 4. (iii) We have already shown that when vanishing T2 , T3 , and T4 the system auto- matically becomes a center, so no isochronous foci exist for this system and the result follows.  Acknowledgements This work has been realized thanks to Spanish Agencia Estatal de Investigación grants PID2019-104658GB-I00 and FPU16/04317; the Severo Ochoa and María de Maeztu Program for Centers and Units of Excellence in R&D grant CEX2020-001084-M; the Generalitat de Catalunya - AGAUR grant 2017SGR1617; the European Community grant Dynamics-H2020-MSCA-RISE-2017-777911; and the Brazilian grants Projeto Temático FAPESP 2019/21181-0 and BP-CNPq 304766/2019-4. Funding Open Access Funding provided by Universitat Autònoma de Barcelona. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. References 1. Algaba, A., Freire, E., Gamero, E.: Isochronicity via normal form. Qual. Theory Dyn. Syst. 1(2), 133–156 (2000) 2. Algaba, A., Freire, E., Gamero, E.: Characterizing and computing normal forms using Lie transforms: a survey. Comput. Normal Forms Appl. 8, 449–475 (2001) 3. Algaba, A., Reyes, M.: Characterizing isochronous pairs and computing isochronous sections. J. Math. Anal. Appl. 355(2), 564–576 (2009) 4. Amelkin, V.V.: Positive solution of one conjecture in the theory of polynomial isochronous centers of liénard systems. Differ. Equ. 57(2), 133–138 (2021) 5. Andronov, A. A., Leontovich, E. A., Gordon,I. I., Maı̆er, A. G.: Theory of bifurcations of dynamic systems on a plane. Halsted Press [A division of John Wiley & Sons], New York-Toronto, Ont.; Israel Program for Scientific Translations, Jerusalem-London, 1973. Translated from the Russian 6. Chicone, C.: The monotonicity of the period function for planar Hamiltonian vector fields. J. Differ. Equ. 69(3), 310–321 (1987) 20 Page 34 of 35 R. D. S. Oliveira et al. 7. Chicone, C., Jacobs, M.: Bifurcation of critical periods for plane vector fields. Trans. Amer. Math. Soc. 312(2), 433–486 (1989) 8. Chow, S.-N., Sanders, J.A.: On the number of critical points of the period. J. Differ. Equ. 64(1), 51–66 (1986) 9. Christopher, C.: An algebraic approach to the classification of centers in polynomial Liénard systems. J. Math. Anal. Appl. 229(1), 319–329 (1999) 10. Cima, A., Gasull, A., Mañosa, V., Mañosas, F.: Algebraic properties of the Liapunov and period constants. Rocky Mountain J. Math. 27(2), 471–501 (1997) 11. Cima, A., Mañosas, F., Villadelprat, J.: Isochronicity for several classes of Hamiltonian systems. J. Differ. Equ. 157(2), 373–413 (1999) 12. Dumortier, F., Llibre, J., Artés, J.C.: Qualitative theory of planar differential systems. Universitext. Springer-Verlag, Berlin (2006) 13. Garijo, A., Gasull, A., Jarque, X.: On the period function for a family of complex differential equations. J. Differ. Equ. 224(2), 314–331 (2006) 14. Garijo, A., Gasull, A., Jarque, X.: A note on the period function for certain planar vector fields. J. Differ. Equ. Appl. 16(5–6), 631–645 (2010) 15. Gasull, A., Mañosa, V., Mañosas, F.: Stability of certain planar unbounded polycycles. J. Math. Anal. Appl. 269(1), 332–351 (2002) 16. Gasull, A., Prohens, R., Torregrosa, J.: Limit cycles for rigid cubic systems. J. Math. Anal. Appl. 303(2), 391–404 (2005) 17. Gasull, A., Torregrosa, J.: Center problem for several differential equations via Cherkas method. J. Math. Anal. Appl. 228(2), 322–343 (1998) 18. Gavrilov, L.: Remark on the number of critical points of the period. J. Differ. Equ. 101(1), 58–65 (1993) 19. Giné, J.: Isochronous foci for analytic differential systems. Internat J. Bifur Chaos Appl. Sci. Engrg. 13(6), 1617–1623 (2003) 20. Giné, J., Grau, M.: Characterization of isochronous foci for planar analytic differential systems. Proc. Roy. Soc. Edinburgh Sect. A 135(5), 985–998 (2005) 21. Giné, J., Llibre, J.: A family of isochronous foci with Darboux first integral. Pacific J. Math. 218(2), 343–355 (2005) 22. Grau, M., Villadelprat, J.: Bifurcation of critical periods from Pleshkans isochrones. J. Lond. Math. Soc. (2) 81(1), 142–160 (2010) 23. Koepf, W.: Hypergeometric summation. Universitext. Springer, London, second edition, 2014. An algorithmic approach to summation and special function identities 24. Mañosa, V.: Liapunov and period constants for a planar vector fields. Universitat Autònoma de Barcelona, Bellaterra, Masterthesis (1995). (In Catalan) 25. Mañosas, F., Villadelprat, J.: A note on the critical periods of potential systems. Internat. J. Bifur. Chaos Appl. Sci. Engrg. 16(3), 765–774 (2006) 26. Mardešić, P., Moser-Jauslin, L., Rousseau, C.: Darboux linearization and isochronous centers with a rational first integral. J. Differ. Equ. 134(2), 216–268 (1997) 27. Marín, D., Villadelprat, J.: Asymptotic expansion of the Dulac map and time for unfoldings of hyper- bolic saddles: general setting. J. Differ. Equ. 275, 684–732 (2021) 28. Romanovski, V.G., Shafer, D.S.: The center and cyclicity problems: a computational algebra approach. Birkhäuser Boston Ltd, Boston, MA (2009) 29. Roussarie, R.: Bifurcations of planar vector fields and Hilbert’s sixteenth problem. Modern Birkhäuser Classics. Birkhäuser/Springer, Basel, 1998. [2013] reprint of the 1998 edition [MR1628014] 30. Sabatini, M.: Characterizing isochronous centres by Lie brackets. Differ. Equ. Dynam. Syst. 5(1), 91–99 (1997) 31. Sánchez-Sánchez, I., Torregrosa, J.: Criticality via first order development of the period constants. Nonlinear Anal. 203, 112164 (2021) 32. Sánchez-Sánchez, I., Torregrosa, J.: New lower bounds of the number of critical periods in reversible centers. J. Differ. Equ. 292, 427–460 (2021) 33. Świrszcz, G.: Cyclicity of infinite contour around certain reversible quadratic center. J. Differ. Equ. 154(2), 239–266 (1999) 34. Zeilberger, D.: A fast algorithm for proving terminating hypergeometric identities. Discrete Math. 80(2), 207–211 (1990) 35. Zeilberger, D.: The method of creative telescoping. J. Symbolic Comput. 11(3), 195–204 (1991) Simultaneous Bifurcation of Limit Cycles and Critical Periods Page 35 of 35 20 36. Zhao, Y.: The monotonicity of period function for codimension four quadratic system Q 14 . J. Differ. Equ. 185(1), 370–387 (2002) Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.