Simultaneous Bifurcation of Limit Cycles and Critical Periods

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.

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), 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 aṡ 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, It is straightforward that from this system we can write dr dϕ = a 2 (ϕ)r 2 + a 3 (ϕ)r 3 where we denote by O j (r ) the monomials in r of at least degree j. From here we deduce the solution being ρ = r (0) the initial condition. Evaluating it at 2π we obtain the so-called Poincaré return map: 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, and by integrating we obtain the expression of the time, the so-called period function, as a Taylor series in ρ, this is 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 variables 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 derivative 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 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][20][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 V 2k+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 T l . 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 For this system, the origin is bi-weak of type [1,3], because we can easily find that the first nonzero Lyapunov constant is V 3 = π and the first nonzero coefficient of the period function is T 3 = π, so the origin is not a center because V 3 = 0 and we have T 3 = 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. 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
being a 2 , a 3 , b 2 , b 3 ∈ R, which have a bi-weak [2,4] focus at the origin. (ii) There exist quartic Liénard systems of the form which have a bi-weak [3,6] focus at the origin. (iii) There exist quadratic systems of the form ẋ = −y + a 20 x 2 + a 11 x y + a 02 y 2 , being a 20 , a 11 , a 02 , b 20 , b 11 , b 02 ∈ R, which have a bi-weak [2,4] focus at the origin. (iv) There exist linear plus cubic homogeneous systems of the form ẋ = −y + a 30 x 3 + a 21 x 2 y + a 12 x y 2 + a 03 y 3 , being a 30 , a 21 , a 12 , a 03 , b 30 , b 21 , b 12 , b 03 ∈ 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 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].

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].

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 F 1 = (F 1 1 , F 2 1 ) and F 2 = (F 1 2 , F 2 2 ) in variables (x, y) ∈ K 2 , where K = R or C, as a new vector field which has the form 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 This definition also appears in [13]. The first proof of the next geometrical equivalence was done by Sabatini in [30]. 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,

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 where U i (ϕ) and W i (ϕ) are homogeneous polynomials in sin ϕ and cos ϕ of degree i + 2. Eliminating time and doing the Taylor series expansion in r we obtain The initial value problem for (16) with the initial condition (r , ϕ) = (ρ, 0) has a unique analytic solution which can be expanded as 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 where V m := u m (2π) is the first coefficient which does not vanish. It is a wellknown fact that the first nonidentically zero V m has odd m ( [5,28]). As a consequence, expression (18) can be rewritten as We will denote the V 2k+1 with odd subscript as L k := V 2k+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 Alternatively, this can be written as 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 for some trigonometric polynomials F i (ϕ). Rewriting this equation as and integrating ϕ from 0 to 2π yields the period function 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 T i of the period function are then given by the expression and the first nonzero T l is known as the lth period constant of the system. If we assume now that system (15) depends on some parameters λ ∈ R d , we can follow exactly the same procedure as before, and now we have that both the Lyapunov constants V 2k+1 = u 2k+1 (2π, λ) and the period constants are polynomials in the parameters λ (see [10,28]). Finally, we consider the case with nonzero trace (α = 0). The following result studies 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
where X and Y are polynomials which start at least with cubic monomials, and such that when α = 0 the first non vanishing Lyapunov and period constants are V 3 and T 2 .
Moreover, if α V 3 < 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 where Moreover, if (e 2πα −1) T 1 (α) T 2 < 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 e 2πα −1 have the same sign, it is immediate to see that if α V 3 < 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, 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 (e 2πα −1) T 1 (α) T 2 < 0, there is an extra change of sign in the coefficients of the period function, which implies the bifurcation of an extra critical period.

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ż where N ∈ N is arbitrary, O 2N +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) Proposition 5 The truncated system (26) in polar coordinates takes the form Proof The proof is straightforward by applying the change of variables (r , ϕ) → z = r e i ϕ , w = r e − i ϕ to differential equation (27). Indeed, 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 Then the Lie bracket between Z and U has the form where O 2N +2 (z, w) denotes a sum of monomials of degree at least 2N + 2 in z and w, and are the coefficients of the return map and the period function of system (26), respectively.
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.
and the second component is Notice that we have denoted by p 2 j+1 the coefficient of the Lie bracket with degree 2 j +1, and it has the expression p 2 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 V 2 j+1 and T 2 j , but only under the conditions V 2i+1 = T 2i = 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 V 2 j+1 the jth Lyapunov constant assuming that V 2i+1 = T 2i = 0 for every i < j, and by T 2 j the jth period constant V 2i+1 = T 2i = 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. y)). Then the following equivalences hold: Proof The proof is straightforward by applying the chain rule for two variable functions. For (30), we use that φ 1 (u, v) = x and φ 2 (u, v) = y and the chain rule to rewrite it as 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. y)). Let us denote by G 1 = (G 1 1 , G 2 1 ) and G 2 = (G 1 2 , G 2 2 ) the vector fields F 1 and F 2 , respectively, in variables (x, y) after applying the change of variables (x, y) = φ (u, v). Then the following equivalence between the Lie brackets of F 1 , F 2 and G 1 , G 2 holds: y), and the superindex T denotes the transpose vector.

Theorem 9 Let us consider two vector fields F
Proof Let us start by observing that equivalence (34) can be rewritten in matrix form as We will prove the first component of equivalence (35), since the second one can be analogously seen. Then, our aim is to show that 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 F 1 and F 2 . Let us denote by G y) for i, j = 1, 2. As vector fields G 1 i and G 2 i 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 since vector fields F 1 i and F 2 i correspond respectively tou andv. The first partial derivative we need to find for [G 1 , G 2 ] 1 is ∂G 1 1 /∂ x and, by applying the chain rule, it is We can find the expression for ∂ G 1 1 /∂u by taking the derivative of (37) with respect to u, and we get ∂u .
To find ∂ G 1 1 /∂v we proceed analogously. Now, by substituting these two expressions in (38), we get ∂G 1 1 /∂ x in terms of the derivatives of F 1 and F 2 as we wanted. The same procedure can be applied to ∂G 1 1 /∂ y, ∂G 1 2 /∂ x, and ∂G 1 2 /∂ y, and then we have all the terms of [G 1 , G 2 ] 1 expressed with F 1 and F 2 .
For the sake of compactness, from now on we will use the following notation to write this expression of [G 1 , G 2 ] 1 : 1 (x, y)), 1 (x, y)).
Also, with a slight abuse of notation, we will denote the derivatives 1 (x, y)) simply as ∂ F j i /∂u and ∂ F j i /∂v, respectively. The expression of [G 1 , G 2 ] 1 can then be written as This expression can be expanded, and a long sum of terms in products of F 1 1 , F 2 1 , F 1 2 , F 2 2 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 F 1 1 F 1 2 cancel with each other, as well as those terms with F 2 1 F 2 2 . Now let us focus on the terms with F 1 1 F 2 2 and F 2 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 F 1 1 F 2 2 vanishes in (39). Analogously, by applying Lemma 8, we can see that the coefficient of F 2 1 F 1 2 also vanishes. Now let us see how the Lie bracket of F 1 and F 2 appears in (39). One can see that where we have applied (30) from Lemma 8. Analogously, applying both (30) and (31) from such lemma, we see that the coefficients of , 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 The same procedure can be followed with terms (39), and we see that these terms equal φ 1v [F 1 , F 2 ] 2 . After all this, expression (39) can be rewritten as 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 To prove the relation for the second component, this is we follow an analogous procedure and see that and the proof finishes. Now let us consider the Lie bracket of two planar real vector fields written (abusing notation) in complex coordinates as (F 1 , F 1 ) and (F 2 , 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 F 1 and F 2 , 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 [F 1 , F 2 ] 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 With this relation, the following corollary from Theorem 9 is straightforward. (F 1 , F 1 ) and (F 2 , 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 (φ, φ).

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 V 2m+1 and T 2m for m ≤ k then V 2k+3 = 0 and T 2k+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 V 3 and T 2 by using the Lie bracket method, System The next step is to find V 5 and T 4 assuming that conditions in S 1 , we find V 7 = 0 and T 6 = 0, and this case corresponds to system ẋ = −y + a 2 x 2 , which is an isochronous center as we will see in Proposition 12. System (9) under condition S 1 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 V 3 and T 2 are the same that in the cubic case, (41), so they vanish again for S (1) 2 = {a 3 = 2 3 b 2 a 2 , b 3 = 4 9 a 2 2 + 10 9 b 2 2 }. Now we find the next constants, which are . The next step is to find V 7 and T 6 , which under in V 9 and T 8 we obtain V 9 = 0 and T 8 = 0 in both cases, so the maximal bi-weak type is of the form [k, 2k] is [3,6] and this proves Theorem 1 (ii). We notice that case S 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 + r 20 z 2 + r 11 zw + r 02 w 2 , w = − i w + s 02 z 2 + s 11 zw + s 20 w 2 , being s i j = r i j . The first Lyapunov and period constants are 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 V 2k+1 and T 2k in the Bautin type ideal B k−1 : First we want to check whether V 5 and T 4 belong to B 1 = V 3 , T 2 . The expression of V 5 and T 4 in B 1 is 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 V 7 and However, we can see that actually V 2 7 = 0 in B 2 , i.e. V 7 = 0 in the variety of zeros defined by the elements in B 2 , 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 T 6 = 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 T 2 , T 4 , and T 6 are necessary to characterize the isochronicity property. In fact, Let us make the following observation regarding the quadratic case. It is a wellknown fact that the general quadratic family (11) classically unfolds 3 limit cycles, hence the case V 3 = V 5 = 0 with nonvanishing V 7 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.
(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 + r 30 z 3 + r 21 z 2 w + r 12 zw 2 + r 03 w 3 , w = − i w + s 03 z 3 + s 12 z 2 w + s 21 zw 2 + s 30 w 3 , being s i j = r i j . We will proceed analogously to the quadratic case, by studying V 2k+1 and T 2k in the Bautin ideal B k−1 := V 2m+1 , T 2m m≤k−1 .
The first Lyapunov and period constants are The next step is to find the expressions of V 5 and T 4 simplified with respect to B 1 = V 3 , T 2 , 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

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.

A Liénard Family Starting with an Odd and an Even Degree Monomials
In this subsection we will consider the Liénard family ẋ = −y + a m x m + a n x n + x d P(x), 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 P d (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 a n = 0 we have a center at the origin. In this case the even power controls the isochronicity property, so if a m = 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 a m and a n 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 V n = a n 2 n (n − 1) 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 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 Binomial Theorem, We will consider the system 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 real vector field, by using (14) their Lie bracket can be described only from its first component: where each H q 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 This homogeneous mth degree part vanishes taking The homogeneous nth degree polynomial H n is analogous to (49), only changing m by n. We can take then as for the term corresponding to p = n−1 2 we have that 2 p − n + 1 = 0, and therefore the coefficient of z = 0. Notice that this fact occurs with n because it is an odd number, but not with m which is even.
Let us see the structure of H 2m−1 . After substituting (50), we obtain 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 For the second summation in (53), the equality 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 After these calculations we can finally prove our result. To this end, we will consider (48) assuming (50) and (51), so H m = 0 and H n only has the z(zw) n− 1 2 term. If n < 2m − 1, then the lowest degree homogeneous polynomial is the nth degree one, H n = a n 2 n (n − 1) = 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 H 2m−1 , so the period constant is the imaginary part of (56) and all Lyapunov constants vanish up to this order.

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.
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 and we compute that the only nontrivial solution is , 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) → ( 3 2a 2 x, 3 2a 2 y) for a 2 = 0, which results in Observe that for the case a 2 = 0, the system is isochronous because it corresponds to the linear center. System (58) is isochronous because it commutes with ẋ = x A(x, y), where A(x, y) = −1 − y + 1 2 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 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π.
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 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 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 = 2h 2 y + 3h 2 − 2y + 1, or equivalently y = 3h 2 −z 2 +1 2(1−h 2 ) . After applying this change both to X + (y, h) and the integration limits y ± , the integral becomes 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 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 , 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.

Proposition 13
System (60) with even n, b = n 2 (n+1) 2 a 2 , and a = 0 has a first integral of the form being A(x, y) a polynomial in x, y.
Proof For b = n 2 (n+1) 2 a 2 and a = 0, system (60) can be rewritten as ẋ = −y + (n + 1)x n =: P(x, y), after the rescaling (x, y) → n+1 a 1 n−1 x, n+1 a 1 n−1 y . It can be checked that function 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 and R(x, y) = x 2 + y 2 − 2x n y + x 2n − n+1 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 ∂ H ∂ x = Q(x, y)R(x, y) and ∂ H ∂ y = −P(x, y)R(x, y), so we can integrate this second equation to find the form of such first integral, The first term in (63) is an immediate integral and its result is , taking into account that due to being a first integral we can consider that the integration constant is 0. For the second integral we will perform the change of variables ω = x √ x 2 +y 2 −2x n y+x 2n , so dy = − x ω 2 √ 1−ω 2 dω and the integral can be written as x n Now we can apply a trigonometric change of variables ω = sin φ, and As n is even, so is n − 2 and the integral becomes 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 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 , so joining both terms from (63) the result follows. 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 and compute 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 By considering (64) on a level curve such that H (x, y) = h 2 2 being h > 0, we obtain that y 2 = h 2 −x 2 − b n x 2n , and we can use this relation to perform a change of variables from y to x in (65) as follows: where we have used the symmetry of the integral and x − h , x + h = x h , are the intersections of the level curve with the horizontal axis, i.e. the real solutions of H (x, 0) = h 2 2 or equivalently x 2 + b n 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 As we are not able to find the explicit solution of L(h) := x 2 + b n 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 as for this x h we have Now by substituting (68) in (67) and developing the expression we have where we can expand the denominator as a power series and obtain Notice that on the last equality we have used the formula for the sum of terms in a geometric progression. Finally, this integral becomes = 2π − 2b n h 2n−2 1 0 1 + z 2 + z 4 + · · · + z 2n−2 √ 1 − z 2 dz + O 4n−4 (h).
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, we will see that if b = n 2 (n+1) 2 a 2 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 b = n 2 (n+1) 2 a 2 then the period function is not constant. Let us assume that b = n 2 (n+1) 2 a 2 , and consider the rescaled system (62). According to Proposition 13, this system has a first integral of the form (61). Observe that, due to the fact of being a first integral, the condition ∂ H ∂ x P + ∂ H ∂ y Q = 0 must be satisfied. In our case, when finding the partial derivatives of the first integral (61) we obtain 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 Let us consider system ẋ = x 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. Finally, we have to check that if b = n 2 (n+1) 2 a 2 , then the center is not isochronous. To this end, we will use the Lie bracket method to find the first period constant and check that it only vanishes for b = n 2 (n+1) 2 a 2 . System (60) can be rewritten in complex coordinates asż = i z + a z + w 2 which is actually system (46) choosing a m = a, a n = 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 By Theorem 6, the imaginary part of (72) is the period constant T 2n−2 = b 2 2n−1 (n − 1) 2n n − a 2 2 2n−1 n 2 (n − 1) (n + 1) 2 2n n = b − n 2 (n + 1) 2 a 2 n − 1 2 2n−1 2n n , which only vanishes for b = n 2 (n+1) 2 a 2 .

Observation 15
After presenting the new isochronous family (60) for even n and b = n 2 (n+1) 2 a 2 , a natural question is to wonder if these isochronous centers will 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 perturbation 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. 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.

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: